Reference Manual

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