5 #include "TLorentzVector.h"
12 std::vector<math::XYZVector>
makeVecForEventShape(std::vector<pat::Jet>
jets,
bool only6Jets =
true, ROOT::Math::Boost boost = ROOT::Math::Boost(0.,0.,0.)) {
13 std::vector<math::XYZVector>
p;
14 bool doBoost = (boost == ROOT::Math::Boost(0.,0.,0.)) ?
false :
true;
15 for(std::vector<pat::Jet>::const_iterator
jet = jets.begin();
jet != jets.end(); ++
jet){
18 ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > Ljet(
jet->px(),
jet->py(),
jet->pz(),
jet->energy());
19 Vjet =
math::XYZVector(boost(Ljet).Px(), boost(Ljet).Py(), boost(Ljet).Pz());
24 if(only6Jets &&
jet-jets.begin()==5)
break;
92 pts_ = std::vector<double>(0);
94 thetas_ = std::vector<double>(0);
125 dR_ = std::vector<double>(0);
126 dRMass_ = std::vector<double>(0);
134 ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > totalSystem(0.,0.,0.,0.);
136 std::vector< std::pair< double, std::vector<unsigned short> > > dRs(0);
137 std::vector< std::pair< double, std::vector<unsigned short> > > dRs3Jets(0);
139 std::vector< std::pair< double, std::vector<unsigned short> > > M3s(0);
141 std::vector<ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > > fourVectors(0);
143 unsigned short nonBJets = 0;
144 for(std::vector<pat::Jet>::const_iterator
jet = jets.begin();
jet != jets.end(); ++
jet){
149 if(
jet - jets.begin() < 3){
152 if(
jet - jets.begin() == 4 ||
jet - jets.begin() == 5)
155 if(
jet - jets.begin() > 1){
163 TCHP_BJet_Discs_ .push_back(
jet->bDiscriminator(
"trackCountingHighPurBJetTags") );
164 SSVHE_BJet_Discs_ .push_back(
jet->bDiscriminator(
"simpleSecondaryVertexHighEffBJetTags") );
165 SSVHP_BJet_Discs_ .push_back(
jet->bDiscriminator(
"simpleSecondaryVertexHighPurBJetTags") );
166 CSV_BJet_Discs_ .push_back(
jet->bDiscriminator(
"combinedSecondaryVertexBJetTags") );
167 CSVMVA_BJet_Discs_.push_back(
jet->bDiscriminator(
"combinedSecondaryVertexMVABJetTags") );
168 SM_BJet_Discs_ .push_back(
jet->bDiscriminator(
"softMuonBJetTags") );
170 pts_ .push_back(
jet->pt());
171 EtSin2Thetas_.push_back(
jet->et()*
pow(
sin(
jet->theta()),2));
172 thetas_ .push_back( (
jet->theta() >
M_PI/2.)? (
M_PI -
jet->theta()) :
jet->theta() );
174 fourVectors.push_back(
jet->p4());
176 if(
jet->bDiscriminator(
"trackCountingHighEffBJetTags") > 3.3 ) ++
TCHE_Bjets_;
177 if(
jet->bDiscriminator(
"trackCountingHighPurBJetTags") > 3.41 ) ++
TCHP_Bjets_;
178 if(
jet->bDiscriminator(
"simpleSecondaryVertexHighEffBJetTags") > 1.74 ) ++
SSVHE_Bjets_;
179 if(
jet->bDiscriminator(
"simpleSecondaryVertexHighPurBJetTags") > 2.0 ) ++
SSVHP_Bjets_;
180 if(
jet->bDiscriminator(
"combinedSecondaryVertexBJetTags") > 0.75 ) ++
CSV_Bjets_;
181 if(
jet->bDiscriminator(
"combinedSecondaryVertexMVABJetTags") > 0.75 ) ++
CSVMVA_Bjets_;
182 if(
jet->bDiscriminator(
"softMuonBJetTags") > 0.3 ) ++
SM_Bjets_;
184 if(
jet->nConstituents() > 0){
187 etaphiMoments_.push_back(
std::abs(
jet->etaphiMoment()));
188 phiphiMoments_.push_back(
jet->phiphiMoment() );
194 etaetaMomentsLogEt_.push_back(
jet->etaetaMoment() *
log(
jet->et()));
196 phiphiMomentsLogEt_.push_back(
jet->phiphiMoment() *
log(
jet->et()));
202 if(
jet->bDiscriminator(
"trackCountingHighEffBJetTags") < 3.3 &&
jet->bDiscriminator(
"trackCountingHighPurBJetTags") < 1.93 &&
203 jet->bDiscriminator(
"simpleSecondaryVertexHighEffBJetTags") < 1.74 &&
jet->bDiscriminator(
"simpleSecondaryVertexHighPurBJetTags") < 2.0 ){
207 etaetaMomentsNoB_.push_back(
jet->etaetaMoment() );
208 etaphiMomentsNoB_.push_back(
std::abs(
jet->etaphiMoment()));
209 phiphiMomentsNoB_.push_back(
jet->phiphiMoment() );
218 for(std::vector<pat::Jet>::const_iterator jet2 =
jet+1; jet2 != jets.end(); ++jet2){
219 unsigned short comb2A[2] = { (
unsigned short)(
jet-jets.begin()) , (
unsigned short)(jet2-jets.begin()) };
220 std::vector<unsigned short> comb2(comb2A, comb2A +
sizeof(comb2A) /
sizeof(
unsigned short));
221 dRs.push_back( std::make_pair(
deltaR(
jet->phi(),
jet->eta(), jet2->phi(), jet2->eta() ), comb2 ) );
223 for(std::vector<pat::Jet>::const_iterator jet3 = jet2+1; jet3 != jets.end(); ++jet3){
224 unsigned short comb3A[3] = { (
unsigned short)(
jet-jets.begin()) , (
unsigned short)(jet2-jets.begin()) , (
unsigned short)(jet3-jets.begin()) };
225 std::vector<unsigned short> comb3(comb3A, comb3A +
sizeof(comb3A) /
sizeof(
unsigned short));
228 double dR3 =
deltaR( jet2->eta(), jet2->phi(), jet3->eta(), jet3->phi() );
229 dRs3Jets.push_back( std::make_pair( dR1+dR2+dR3, comb3 ) );
230 M3s.push_back( std::make_pair( (
jet->p4() + jet2->p4() + jet3->p4() ).
pt(), comb3 ) );
234 totalSystem +=
jet->p4();
237 ROOT::Math::Boost CoMBoostTotal(totalSystem.BoostToCM());
238 std::vector<reco::LeafCandidate> boostedJets;
240 for(std::vector<pat::Jet>::const_iterator
jet = jets.begin();
jet != jets.end(); ++
jet){
244 EtSin2Theta3jet_ /= ((double)(jets.size()-3));
273 std::sort(dRs3Jets.begin(), dRs3Jets.end());
275 for(std::vector< std::pair<
double, std::vector<unsigned short> > >::const_iterator
dR = dRs.begin();
dR != dRs.end(); ++
dR){
276 dR_.push_back(
dR->first);
277 dRMass_.push_back((jets.at(
dR->second.at(0)).
p4()+jets.at(
dR->second.at(1)).
p4()).mass());
279 ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > wHypo = jets.at(
dR->second.at(0)).
p4()+jets.at(
dR->second.at(1)).
p4();
280 TLorentzVector wHypoHelper(wHypo.Px(), wHypo.Py(), wHypo.Pz(), wHypo.E());
281 wHypoHelper.SetVectM(TVector3(wHypo.Px(), wHypo.Py(), wHypo.Pz()), 80.4);
282 wHypo.SetPxPyPzE(wHypoHelper.Px(), wHypoHelper.Py(), wHypoHelper.Pz(), wHypoHelper.E());
283 ROOT::Math::Boost CoMBoostWHypo(wHypo.BoostToCM());
284 dRAngle_.push_back(ROOT::Math::VectorUtil::Angle(CoMBoostWHypo(jets.at(
dR->second.at(0)).
p4()), CoMBoostWHypo(jets.at(
dR->second.at(1)).
p4())));
287 for(std::vector< std::pair<
double, std::vector<unsigned short> > >::const_iterator
dR = dRs3Jets.begin();
dR != dRs3Jets.end(); ++
dR){
288 dR3Jets_.push_back(
dR->first);
289 dR3JetsMass_.push_back((jets.at(
dR->second.at(0)).
p4()+jets.at(
dR->second.at(1)).
p4()+jets.at(
dR->second.at(2)).
p4()).mass());
292 std::vector< std::pair< double, unsigned short > > massDiff2W;
294 for(std::vector< double >::const_iterator mass = dRMass_.begin(); mass != dRMass_.end(); ++mass){
295 massDiff2W.push_back(std::make_pair(
std::abs((*mass)-80.4), mass - dRMass_.begin()));
298 std::sort(massDiff2W.begin(), massDiff2W.end());
302 for(std::vector< std::pair< double, unsigned short > >::const_iterator
i = massDiff2W.begin();
i != massDiff2W.end(); ++
i){
303 unsigned int mass1 =
i->second;
304 for(std::vector< std::pair< double, unsigned short > >::const_iterator
j =
i + 1;
j != massDiff2W.end(); ++
j){
305 unsigned int mass2 =
j->second;
306 if(dRs.at(mass1).second.at(0) != dRs.at(mass2).second.at(0) && dRs.at(mass1).second.at(0) != dRs.at(mass2).second.at(1) &&
307 dRs.at(mass1).second.at(1) != dRs.at(mass2).second.at(0) && dRs.at(mass1).second.at(1) != dRs.at(mass2).second.at(1)){
311 massDiffMWCands_.push_back(
std::abs(dRMass_.at(mass1)-dRMass_.at(mass2)));
314 if(massDiffMWCands_.size() > 20)
break;
332 std::sort(CSVMVA_BJet_Discs_.begin(), CSVMVA_BJet_Discs_.end());
336 std::sort(etaphiMoments_.begin(), etaphiMoments_.end());
337 std::sort(phiphiMoments_.begin(), phiphiMoments_.end());
339 std::sort(etaetaMomentsLogEt_.begin(), etaetaMomentsLogEt_.end());
340 std::sort(etaphiMomentsLogEt_.begin(), etaphiMomentsLogEt_.end());
341 std::sort(phiphiMomentsLogEt_.begin(), phiphiMomentsLogEt_.end());
343 std::sort(etaetaMomentsMoment_.begin(), etaetaMomentsMoment_.end());
344 std::sort(etaphiMomentsMoment_.begin(), etaphiMomentsMoment_.end());
345 std::sort(phiphiMomentsMoment_.begin(), phiphiMomentsMoment_.end());
347 std::sort(etaetaMomentsMomentLogEt_.begin(), etaetaMomentsMomentLogEt_.end());
348 std::sort(etaphiMomentsMomentLogEt_.begin(), etaphiMomentsMomentLogEt_.end());
349 std::sort(phiphiMomentsMomentLogEt_.begin(), phiphiMomentsMomentLogEt_.end());
351 std::sort(etaetaMomentsNoB_.begin(), etaetaMomentsNoB_.end());
352 std::sort(etaphiMomentsNoB_.begin(), etaphiMomentsNoB_.end());
353 std::sort(phiphiMomentsNoB_.begin(), phiphiMomentsNoB_.end());
356 M3_ = ( jets.at((M3s.back().second.at(0))).p4() + jets.at((M3s.back().second.at(1))).p4() + jets.at((M3s.back().second.at(2))).p4() ).mass();
361 ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > boostedJet = CoMBoostTotal( *
jet );
362 if(
jet - fourVectors.begin() > 1){
367 thetaStars_.push_back((boostedJet.Theta() >
M_PI/2.)? (
M_PI - boostedJet.Theta()) : boostedJet.Theta());
368 EtStars_.push_back(
jet->Et() *
pow(
sin( boostedJet.Theta() ), 2) );
371 theta3jet_ /= (double)fourVectors.size() - 2.;
391 CAll_ = eventshapeAll.
C();
392 DAll_ = eventshapeAll.
D();
403 Thrust thrustAlgo(jets.begin(), jets.end());
406 Thrust thrustAlgoCMS(boostedJets.begin(), boostedJets.end());
std::vector< double > etaetaMoments_
std::vector< double > TCHP_BJet_Discs_
double C(double=2.) const
std::vector< double > dR_
std::vector< double > massDiffMWCands_
double circularity(const unsigned int &numberOfSteps=1000) const
double jets_phiphiMomentNoB_
std::vector< double > etaphiMoments_
std::vector< double > phiphiMomentsMoment_
std::vector< double > phiphiMomentsLogEt_
Class for the calculation of several event shape variables.
std::vector< double > etaphiMomentsLogEt_
std::vector< double > TCHE_BJet_Discs_
std::vector< double > phiphiMomentsNoB_
std::vector< double > dR3JetsMass_
Sin< T >::type sin(const T &t)
std::vector< double > EtStars_
std::vector< double > etaetaMomentsMomentLogEt_
double jets_etaphiMoment_
double pt(unsigned short i) const
std::vector< double > dR3Jets_
std::vector< double > etaphiMomentsMomentLogEt_
std::vector< double > etaetaMomentsLogEt_
double jets_phiphiMomentLogEt_
std::vector< math::XYZVector > makeVecForEventShape(std::vector< pat::Jet > jets, bool only6Jets=true, ROOT::Math::Boost boost=ROOT::Math::Boost(0., 0., 0.))
double jets_phiphiMoment_
std::vector< double > SM_BJet_Discs_
double jets_etaetaMoment_
std::vector< double > EtSin2Thetas_
double sphericity(double=2.) const
std::vector< double > etaetaMomentsNoB_
double deltaR(double eta1, double eta2, double phi1, double phi2)
std::vector< double > thetas_
std::vector< double > pts_
double D(double=2.) const
std::vector< double > dRAngle_
std::vector< double > thetaStars_
XYZVectorD XYZVector
spatial vector with cartesian internal representation
double isotropy(const unsigned int &numberOfSteps=1000) const
double jets_etaetaMomentNoB_
std::vector< double > SSVHP_BJet_Discs_
std::vector< double > CSV_BJet_Discs_
std::vector< double > SSVHE_BJet_Discs_
double jets_etaetaMomentLogEt_
std::vector< double > CSVMVA_BJet_Discs_
std::vector< double > etaphiMomentsNoB_
double jets_etaphiMomentNoB_
std::vector< double > etaphiMomentsMoment_
std::vector< double > etaetaMomentsMoment_
std::vector< double > phiphiMoments_
std::vector< double > phiphiMomentsMomentLogEt_
Power< A, B >::type pow(const A &a, const B &b)
double aplanarity(double=2.) const
double circularityAllCMS_
math::PtEtaPhiELorentzVectorF LorentzVector
std::vector< double > dRMass_
double jets_etaphiMomentLogEt_