![]() |
![]() |
#include <TopQuarkAnalysis/TopEventSelection/interface/TtHadLRSignalSelObservables.h>
Public Member Functions | |
void | operator() (TtHadEvtSolution &) |
TtHadLRSignalSelObservables () | |
~TtHadLRSignalSelObservables () | |
Private Attributes | |
CompareBdisc | BdiscComparator |
CompareDouble | dComparator |
CompareET | EtComparator |
vector< pair< unsigned int, double > > | evtselectVarVal |
Classes | |
struct | CompareBdisc |
struct | CompareDouble |
struct | CompareET |
Definition at line 28 of file TtHadLRSignalSelObservables.h.
TtHadLRSignalSelObservables::TtHadLRSignalSelObservables | ( | ) |
TtHadLRSignalSelObservables::~TtHadLRSignalSelObservables | ( | ) |
void TtHadLRSignalSelObservables::operator() | ( | TtHadEvtSolution & | TS | ) |
Definition at line 15 of file TtHadLRSignalSelObservables.cc.
References BdiscComparator, funct::C, funct::cos(), funct::D, dComparator, relval_parameters_module::energy, evtselectVarVal, TtHadEvtSolution::getHadb(), TtHadEvtSolution::getHadbbar(), TtHadEvtSolution::getHadj(), TtHadEvtSolution::getHadk(), TtHadEvtSolution::getHadp(), TtHadEvtSolution::getHadq(), i, edm::isnan(), funct::log(), N, PI, funct::pow(), TtHadEvtSolution::setLRSignalEvtObservables(), funct::sin(), python::multivaluedict::sort(), funct::sqrt(), and muonGeometry::TVector3().
00016 { 00017 evtselectVarVal.clear(); 00018 vector<pat::Jet> topJets; 00019 topJets.clear(); 00020 topJets.push_back(TS.getHadp()); 00021 topJets.push_back(TS.getHadq()); 00022 topJets.push_back(TS.getHadj()); 00023 topJets.push_back(TS.getHadk()); 00024 topJets.push_back(TS.getHadb()); 00025 topJets.push_back(TS.getHadbbar()); 00026 00027 //Assume the highest Et jets are the b-jets, the others the jets from the W - but how to work out which two are from which W? FIXME!!! for now don't sort jets in Et before after this... 00028 TLorentzVector *Hadp = new TLorentzVector(); 00029 TLorentzVector *Hadq = new TLorentzVector(); 00030 TLorentzVector *Hadj = new TLorentzVector(); 00031 TLorentzVector *Hadk = new TLorentzVector(); 00032 TLorentzVector *Hadb = new TLorentzVector(); 00033 TLorentzVector *Hadbbar = new TLorentzVector(); 00034 Hadp->SetPxPyPzE(topJets[0].px(),topJets[0].py(),topJets[0].pz(),topJets[3].energy()); 00035 Hadq->SetPxPyPzE(topJets[1].px(),topJets[1].py(),topJets[1].pz(),topJets[2].energy()); 00036 Hadj->SetPxPyPzE(topJets[2].px(),topJets[2].py(),topJets[2].pz(),topJets[2].energy()); 00037 Hadk->SetPxPyPzE(topJets[3].px(),topJets[3].py(),topJets[3].pz(),topJets[3].energy()); 00038 Hadb->SetPxPyPzE(topJets[4].px(),topJets[4].py(),topJets[4].pz(),topJets[1].energy()); 00039 Hadbbar->SetPxPyPzE(topJets[4].px(),topJets[4].py(),topJets[4].pz(),topJets[1].energy()); 00040 00041 //sort the topJets in Et 00042 std::sort(topJets.begin(),topJets.end(),EtComparator); 00043 00044 //Et-Sum of the lightest jets 00045 double EtSum = topJets[5].et()+topJets[5].et(); 00046 double Obs1 = (EtSum>0 ? EtSum : -1); 00047 evtselectVarVal.push_back(pair<unsigned int,double>(1,Obs1)); 00048 00049 //Log of the difference in Bdisc between the 2nd and the 3rd jets (ordered in Bdisc) - does that still 00050 //make sense for the fully hadronic channel as well?FIXME!!! 00051 00052 //sort the jets in Bdiscriminant 00053 std::sort(topJets.begin(),topJets.end(),BdiscComparator); 00054 00055 double BGap = topJets[1].bDiscriminator("trackCountingJetTags") - topJets[2].bDiscriminator("trackCountingJetTags"); 00056 double Obs2 = (BGap>0 ? log(BGap) : -1); 00057 evtselectVarVal.push_back(pair<unsigned int,double>(2,Obs2)); 00058 00059 //Circularity of the event 00060 double N=0,D=0,C_tmp=0,C=1000; 00061 double nx,ny,nz; 00062 00063 // C = 2min(E(pt.n)^2/E(pt)^2) = 2*N/D but it is theorically preferable to use C'=PI/2*min(E|pt.n|/E|pt|), sum over all jets+lepton+MET (cf PhysRevD 48 R3953(Nov 1993)) 00064 00065 for(unsigned int i=0;i<6;i++){ 00066 D += fabs(topJets[i].pt()); 00067 } 00068 00069 if((D>0)){ 00070 // Loop over all the unit vectors in the transverse plane in order to find the miminum : 00071 for(unsigned int i=0; i<360; i++){ 00072 nx = cos((2*PI/360)*i); 00073 ny = sin((2*PI/360)*i); 00074 nz = 0; 00075 N=0; 00076 for(unsigned int i=0;i<4;i++){ 00077 N += fabs(topJets[i].px()*nx+topJets[i].py()*ny+topJets[i].pz()*nz); 00078 } 00079 C_tmp = 2*N/D; 00080 if(C_tmp<C) C = C_tmp; 00081 } 00082 } 00083 00084 double Obs3 = ( C!=1000 ? C : -1); 00085 evtselectVarVal.push_back(pair<unsigned int,double>(3,Obs3)); 00086 00087 00088 //HT variable (Et-sum of the six jets) 00089 double HT=0; 00090 for(unsigned int i=0;i<6;i++){ 00091 HT += topJets[i].et(); 00092 } 00093 00094 double Obs4 = ( HT!=0 ? HT : -1); 00095 evtselectVarVal.push_back(pair<unsigned int,double>(4,Obs4)); 00096 00097 //Transverse Mass of the system 00098 XYZTLorentzVector pjets; 00099 // for the six jets 00100 for(unsigned int i=0;i<6;i++){ 00101 pjets += topJets[i].p4(); 00102 } 00103 double MT = sqrt(pow(pjets.mass(),2)); 00104 double Obs5 = ( MT>0 ? MT : -1); 00105 evtselectVarVal.push_back(pair<unsigned int,double>(5,Obs5)); 00106 00107 //CosTheta(Hadp,Hadq) and CosTheta(Hadj,Hadq) 00108 //sort the lightJets in Et 00109 std::sort(topJets.begin(),topJets.end(),EtComparator); 00110 00111 double px1 = topJets[2].px(); double px2 = topJets[3].px(); 00112 double py1 = topJets[2].py(); double py2 = topJets[3].py(); 00113 double pz1 = topJets[2].pz(); double pz2 = topJets[3].pz(); 00114 double E1 = topJets[2].energy(); double E2 = topJets[3].energy(); 00115 double px3 = topJets[4].px(); double px4 = topJets[5].px(); 00116 double py3 = topJets[4].py(); double py4 = topJets[5].py(); 00117 double pz3 = topJets[4].pz(); double pz4 = topJets[5].pz(); 00118 double E3 = topJets[4].energy(); double E4 = topJets[5].energy(); 00119 00120 TLorentzVector *LightJet1 = new TLorentzVector(); 00121 TLorentzVector *LightJet2 = new TLorentzVector(); 00122 TLorentzVector *LightJet3 = new TLorentzVector(); 00123 TLorentzVector *LightJet4 = new TLorentzVector(); 00124 00125 LightJet1->SetPxPyPzE(px1,py1,pz1,E1); 00126 //LightJet1->Boost(BoostBackToCM); 00127 LightJet2->SetPxPyPzE(px2,py2,pz2,E2); 00128 //LightJet2->Boost(BoostBackToCM); 00129 LightJet3->SetPxPyPzE(px3,py3,pz3,E3); 00130 //LightJet1->Boost(BoostBackToCM); 00131 LightJet4->SetPxPyPzE(px4,py4,pz4,E4); 00132 //LightJet2->Boost(BoostBackToCM); 00133 00134 double CosTheta1 = cos(LightJet2->Angle(LightJet1->Vect())); 00135 double CosTheta2 = cos(LightJet4->Angle(LightJet3->Vect())); 00136 00137 double Obs6 = ( -1<CosTheta1 ? CosTheta1 : -2); 00138 evtselectVarVal.push_back(pair<unsigned int,double>(6,Obs6)); 00139 double Obs7 = ( -1<CosTheta2 ? CosTheta2 : -2); 00140 evtselectVarVal.push_back(pair<unsigned int,double>(7,Obs7)); 00141 00142 delete LightJet1; 00143 delete LightJet2; 00144 delete LightJet3; 00145 delete LightJet4; 00146 00147 // try to find out more powerful observables related to the b-disc 00148 //sort the jets in Bdiscriminant 00149 std::sort(topJets.begin(),topJets.end(),BdiscComparator); 00150 00151 double BjetsBdiscSum = topJets[0].bDiscriminator("trackCountingJetTags") + topJets[1].bDiscriminator("trackCountingJetTags"); 00152 double Ljets1BdiscSum = topJets[2].bDiscriminator("trackCountingJetTags") + topJets[3].bDiscriminator("trackCountingJetTags"); 00153 double Ljets2BdiscSum = topJets[4].bDiscriminator("trackCountingJetTags") + topJets[5].bDiscriminator("trackCountingJetTags"); 00154 //cout<<"BjetsBdiscSum = "<<BjetsBdiscSum<<endl; 00155 //cout<<"LjetsBdiscSum = "<<LjetsBdiscSum<<endl; 00156 00157 double Obs8 = (Ljets1BdiscSum !=0 ? log(BjetsBdiscSum/Ljets1BdiscSum) : -1); 00158 evtselectVarVal.push_back(pair<unsigned int,double>(8,Obs8)); 00159 double Obs9 = (Ljets2BdiscSum !=0 ? log(BjetsBdiscSum/Ljets2BdiscSum) : -1); 00160 evtselectVarVal.push_back(pair<unsigned int,double>(9,Obs9)); 00161 double Obs10 = (BGap>0 ? log(BjetsBdiscSum*BGap) : -1); 00162 evtselectVarVal.push_back(pair<unsigned int,double>(10,Obs10)); 00163 00164 00165 // Et-Ratio between the two highest in Et jets and four highest jets 00166 double Obs11 = (HT!=0 ? (HT-EtSum)/HT : -1); 00167 evtselectVarVal.push_back(pair<unsigned int,double>(11,Obs11)); 00168 00169 00170 //Sphericity and Aplanarity without boosting back the system to CM frame 00171 00172 TMatrixDSym Matrix(3); 00173 double PX2 = pow(Hadp->Px(),2)+pow(Hadq->Px(),2)+pow(Hadb->Px(),2)+ 00174 pow(Hadj->Px(),2)+pow(Hadk->Px(),2)+pow(Hadbbar->Px(),2); 00175 double PY2 = pow(Hadp->Py(),2)+pow(Hadq->Py(),2)+pow(Hadb->Py(),2)+ 00176 pow(Hadj->Py(),2)+pow(Hadk->Py(),2)+pow(Hadbbar->Py(),2); 00177 double PZ2 = pow(Hadp->Pz(),2)+pow(Hadq->Pz(),2)+pow(Hadb->Pz(),2)+ 00178 pow(Hadj->Pz(),2)+pow(Hadk->Pz(),2)+pow(Hadbbar->Pz(),2); 00179 double P2 = PX2+PY2+PZ2; 00180 00181 double PXY = Hadp->Px()*Hadp->Py()+Hadq->Px()*Hadq->Py()+Hadb->Px()*Hadb->Py()+ 00182 Hadj->Px()*Hadj->Py()+Hadk->Px()*Hadk->Py()+Hadbbar->Px()*Hadbbar->Py(); 00183 double PXZ = Hadp->Px()*Hadp->Pz()+Hadq->Px()*Hadq->Pz()+Hadb->Px()*Hadb->Pz()+ 00184 Hadj->Px()*Hadj->Pz()+Hadk->Px()*Hadk->Pz()+Hadbbar->Px()*Hadbbar->Pz(); 00185 double PYZ = Hadp->Py()*Hadp->Pz()+Hadq->Py()*Hadq->Pz()+Hadb->Py()*Hadb->Pz()+ 00186 Hadj->Py()*Hadj->Pz()+Hadk->Py()*Hadk->Pz()+Hadbbar->Py()*Hadbbar->Pz(); 00187 00188 Matrix(0,0) = PX2/P2; Matrix(0,1) = PXY/P2; Matrix(0,2) = PXZ/P2; 00189 Matrix(1,0) = PXY/P2; Matrix(1,1) = PY2/P2; Matrix(1,2) = PYZ/P2; 00190 Matrix(2,0) = PXZ/P2; Matrix(2,1) = PYZ/P2; Matrix(2,2) = PZ2/P2; 00191 00192 TMatrixDSymEigen pTensor(Matrix); 00193 00194 vector<double> EigValues; 00195 EigValues.clear(); 00196 for(int i=0;i<3;i++){ 00197 EigValues.push_back(pTensor.GetEigenValues()[i]); 00198 } 00199 00200 std::sort(EigValues.begin(),EigValues.end(),dComparator); 00201 00202 double Sphericity = 1.5*(EigValues[1]+EigValues[2]); 00203 double Aplanarity = 1.5*EigValues[2]; 00204 00205 double Obs12 = (isnan(Sphericity) ? -1 : Sphericity); 00206 evtselectVarVal.push_back(pair<unsigned int,double>(12,Obs12)); 00207 00208 double Obs13 = (isnan(Aplanarity) ? -1 : Aplanarity); 00209 evtselectVarVal.push_back(pair<unsigned int,double>(13,Obs13)); 00210 00211 00212 //Sphericity and Aplanarity with boosting back the system to CM frame 00213 TLorentzVector *TtbarSystem = new TLorentzVector(); 00214 TtbarSystem->SetPx(Hadp->Px()+Hadq->Px()+Hadb->Px()+Hadj->Px()+Hadk->Px()+Hadbbar->Px()); 00215 TtbarSystem->SetPy(Hadp->Py()+Hadq->Py()+Hadb->Py()+Hadj->Py()+Hadk->Py()+Hadbbar->Py()); 00216 TtbarSystem->SetPz(Hadp->Pz()+Hadq->Pz()+Hadb->Pz()+Hadj->Pz()+Hadk->Pz()+Hadbbar->Pz()); 00217 TtbarSystem->SetE(Hadp->E()+Hadq->E()+Hadb->E()+Hadj->E()+Hadk->E()+Hadbbar->E()); 00218 00219 TVector3 BoostBackToCM = -(TtbarSystem->BoostVector()); 00220 Hadp->Boost(BoostBackToCM); 00221 Hadq->Boost(BoostBackToCM); 00222 Hadb->Boost(BoostBackToCM); 00223 Hadj->Boost(BoostBackToCM); 00224 Hadk->Boost(BoostBackToCM); 00225 Hadbbar->Boost(BoostBackToCM); 00226 00227 double BOOST_PX2 = pow(Hadp->Px(),2)+pow(Hadq->Px(),2)+pow(Hadb->Px(),2)+ 00228 pow(Hadj->Px(),2)+pow(Hadk->Px(),2)+pow(Hadbbar->Px(),2); 00229 double BOOST_PY2 = pow(Hadp->Py(),2)+pow(Hadq->Py(),2)+pow(Hadb->Py(),2)+ 00230 pow(Hadj->Py(),2)+pow(Hadk->Py(),2)+pow(Hadbbar->Py(),2); 00231 double BOOST_PZ2 = pow(Hadp->Pz(),2)+pow(Hadq->Pz(),2)+pow(Hadb->Pz(),2)+ 00232 pow(Hadj->Pz(),2)+pow(Hadk->Pz(),2)+pow(Hadbbar->Pz(),2); 00233 00234 double BOOST_P2 = BOOST_PX2+BOOST_PY2+BOOST_PZ2; 00235 00236 double BOOST_PXY = Hadp->Px()*Hadp->Py()+Hadq->Px()*Hadq->Py()+Hadb->Px()*Hadb->Py()+ 00237 Hadj->Px()*Hadj->Py()+Hadk->Px()*Hadk->Py()+Hadbbar->Px()*Hadbbar->Py(); 00238 double BOOST_PXZ = Hadp->Px()*Hadp->Pz()+Hadq->Px()*Hadq->Pz()+Hadb->Px()*Hadb->Pz()+ 00239 Hadj->Px()*Hadj->Pz()+Hadk->Px()*Hadk->Pz()+Hadbbar->Px()*Hadbbar->Pz(); 00240 double BOOST_PYZ = Hadp->Py()*Hadp->Pz()+Hadq->Py()*Hadq->Pz()+Hadb->Py()*Hadb->Pz()+ 00241 Hadj->Py()*Hadj->Pz()+Hadk->Py()*Hadk->Pz()+Hadbbar->Py()*Hadbbar->Pz(); 00242 00243 TMatrixDSym BOOST_Matrix(3); 00244 00245 BOOST_Matrix(0,0) = BOOST_PX2/BOOST_P2; BOOST_Matrix(0,1) = BOOST_PXY/BOOST_P2; BOOST_Matrix(0,2) = BOOST_PXZ/BOOST_P2; 00246 BOOST_Matrix(1,0) = BOOST_PXY/BOOST_P2; BOOST_Matrix(1,1) = BOOST_PY2/BOOST_P2; BOOST_Matrix(1,2) = BOOST_PYZ/BOOST_P2; 00247 BOOST_Matrix(2,0) = BOOST_PXZ/BOOST_P2; BOOST_Matrix(2,1) = BOOST_PYZ/BOOST_P2; BOOST_Matrix(2,2) = BOOST_PZ2/BOOST_P2; 00248 00249 TMatrixDSymEigen BOOST_pTensor(BOOST_Matrix); 00250 00251 vector<double> BOOST_EigValues; 00252 BOOST_EigValues.clear(); 00253 for(int i=0;i<3;i++){ 00254 BOOST_EigValues.push_back(BOOST_pTensor.GetEigenValues()[i]); 00255 } 00256 00257 std::sort(BOOST_EigValues.begin(),BOOST_EigValues.end(),dComparator); 00258 00259 double BOOST_Sphericity = 1.5*(BOOST_EigValues[1]+BOOST_EigValues[2]); 00260 double BOOST_Aplanarity = 1.5*BOOST_EigValues[2]; 00261 00262 double Obs14 = ( isnan(BOOST_Sphericity) ? -1 : BOOST_Sphericity ); 00263 evtselectVarVal.push_back(pair<unsigned int,double>(14,Obs14)); 00264 00265 double Obs15 = ( isnan(BOOST_Aplanarity) ? -1 : BOOST_Aplanarity ); 00266 evtselectVarVal.push_back(pair<unsigned int,double>(15,Obs15)); 00267 00268 00269 // Centrality of the six jets 00270 double H=0; 00271 for(unsigned int i=0;i<6;i++){ 00272 H += topJets[i].energy(); 00273 } 00274 double Obs16 = ( H != 0 ? HT/H : -1 ); 00275 evtselectVarVal.push_back(pair<unsigned int,double>(16,Obs16)); 00276 00277 00278 // distance from the origin in the (BjetsBdiscSum, Ljets1BdiscSum) and (BjetsBdiscSum, Ljets2BdiscSum) 00279 double Obs17 = ( BjetsBdiscSum != 0 && Ljets1BdiscSum != 0 ? 0.707*(BjetsBdiscSum-Ljets1BdiscSum) : -1 ); 00280 evtselectVarVal.push_back(pair<unsigned int,double>(17,Obs17)); 00281 double Obs18 = ( BjetsBdiscSum != 0 && Ljets2BdiscSum != 0 ? 0.707*(BjetsBdiscSum-Ljets2BdiscSum) : -1 ); 00282 evtselectVarVal.push_back(pair<unsigned int,double>(18,Obs18)); 00283 00284 // Ratio of the Et Sum of the two lightest jets over HT=Sum of the Et of the six highest jets in Et 00285 double Obs19 = ( HT != 0 ? EtSum/HT : -1 ); 00286 evtselectVarVal.push_back(pair<unsigned int,double>(19,Obs19)); 00287 00288 00289 // Put the vector in the TtHadEvtSolution 00290 TS.setLRSignalEvtObservables(evtselectVarVal); 00291 }
Definition at line 47 of file TtHadLRSignalSelObservables.h.
vector<pair<unsigned int,double> > TtHadLRSignalSelObservables::evtselectVarVal [private] |