CMSSW Reference Manual

00001 #include "TopQuarkAnalysis/TopEventSelection/interface/TtHadLRSignalSelObservables.h"
00002 
00003 using namespace reco;
00004 using namespace std;
00005 using namespace math;
00006 
00007 TtHadLRSignalSelObservables::TtHadLRSignalSelObservables()
00008 {
00009 }
00010 
00011 TtHadLRSignalSelObservables::~TtHadLRSignalSelObservables()
00012 {
00013 }
00014 
00015 void TtHadLRSignalSelObservables::operator() (TtHadEvtSolution &TS)
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 }
TtHadLRSignalSelObservables.cc
