Reference Manual

Go to the documentation of this file.
00001 #include "TopQuarkAnalysis/TopEventSelection/interface/TtSemiLRSignalSelObservables.h"
00002 
00003 TtSemiLRSignalSelObservables::TtSemiLRSignalSelObservables()
00004 {
00005 }
00006 
00007 TtSemiLRSignalSelObservables::~TtSemiLRSignalSelObservables()
00008 {
00009 }
00010 
00011 void TtSemiLRSignalSelObservables::operator() (TtSemiEvtSolution &TS,const std::vector<pat::Jet>& SelectedTopJets)
00012 {
00013   evtselectVarVal.clear();
00014   
00015   // choice the B-tag algorithm
00016   const char *BtagAlgo = "trackCountingJetTags";
00017   
00018   // activate the "debug mode"
00019   bool DEBUG = false;
00020   
00021   if(DEBUG) std::cout<<"---------- Start calculating the LR observables ----------"<<std::endl;
00022   
00023   
00024   std::vector<pat::Jet> TopJets;
00025   TopJets.clear();
00026   
00027   for(unsigned int i = 0 ; i<SelectedTopJets.size() ; i++)
00028     {
00029       TopJets.push_back(SelectedTopJets[i]);
00030     }
00031   
00032   //sort the TopJets in Et
00033   std::sort(TopJets.begin(),TopJets.end(),EtComparator);
00034   
00035   TLorentzVector * Hadp = new TLorentzVector();
00036   Hadp->SetPxPyPzE(TopJets[3].px(),TopJets[3].py(),TopJets[3].pz(),TopJets[3].energy());
00037   
00038   TLorentzVector * Hadq = new TLorentzVector();
00039   Hadq->SetPxPyPzE(TopJets[2].px(),TopJets[2].py(),TopJets[2].pz(),TopJets[2].energy());
00040   
00041   TLorentzVector * Hadb = new TLorentzVector();
00042   Hadb->SetPxPyPzE(TopJets[1].px(),TopJets[1].py(),TopJets[1].pz(),TopJets[1].energy());
00043   
00044   TLorentzVector * Lepb = new TLorentzVector();
00045   Lepb->SetPxPyPzE(TopJets[0].px(),TopJets[0].py(),TopJets[0].pz(),TopJets[0].energy());
00046   
00047   TLorentzVector * Lept = new TLorentzVector();
00048   if(TS.getDecay() == "muon") Lept->SetPxPyPzE(TS.getMuon().px(),TS.getMuon().py(),TS.getMuon().pz(),TS.getMuon().energy());
00049   else if(TS.getDecay() == "electron") Lept->SetPxPyPzE(TS.getElectron().px(),TS.getElectron().py(),TS.getElectron().pz(),TS.getElectron().energy());
00050   
00051   TLorentzVector *Lepn = new TLorentzVector();
00052   Lepn->SetPxPyPzE(TS.getNeutrino().px(),TS.getNeutrino().py(),TS.getNeutrino().pz(),TS.getNeutrino().energy());
00053   
00054   // Calculation of the pz of the neutrino due to W-mass constraint
00055   
00056   MEzCalculator *Mez = new MEzCalculator();
00057   Mez->SetMET(TS.getNeutrino());
00058   if(TS.getDecay() == "muon") Mez->SetLepton(TS.getMuon());
00059   else Mez->SetLepton(TS.getElectron(), false);
00060   double MEZ = Mez->Calculate();
00061   Lepn->SetPz(MEZ);
00062   
00063   // Pt of the lepton
00064   
00065   double Obs1 = Lept->Pt();
00066   evtselectVarVal.push_back(std::pair<unsigned int,double>(1,Obs1));
00067   if(DEBUG) std::cout<<"------ LR observable 1 "<<Obs1<<" calculated ------"<<std::endl;
00068   
00069   // Missing transverse energy
00070   
00071   double Obs2 = TS.getNeutrino().et();
00072   evtselectVarVal.push_back(std::pair<unsigned int,double>(2,Obs2));    
00073   if(DEBUG) std::cout<<"------ LR observable 2 "<<Obs2<<" calculated ------"<<std::endl;
00074   
00075   //HT variable (Et-sum of the four jets)
00076   
00077   double HT=0;
00078   for(unsigned int i=0;i<4;i++)
00079     {
00080       HT += TopJets[i].et();
00081     }
00082   
00083   double Obs3 = HT;
00084   evtselectVarVal.push_back(std::pair<unsigned int,double>(3,Obs3));
00085   if(DEBUG) std::cout<<"------ LR observable 3 "<<Obs3<<" calculated ------"<<std::endl;
00086   
00087   //Et-Sum of the lightest jets
00088   
00089   double EtSum = TopJets[2].et()+TopJets[3].et();
00090   double Obs4 = EtSum;
00091   evtselectVarVal.push_back(std::pair<unsigned int,double>(4,Obs4));
00092   if(DEBUG) std::cout<<"------ LR observable 4 "<<Obs4<<" calculated ------"<<std::endl;
00093   
00094   // Et-Ratio between the two lowest jets in Et and four highest jets
00095   
00096   double Obs5 = EtSum/HT;
00097   evtselectVarVal.push_back(std::pair<unsigned int,double>(5,Obs5));
00098   if(DEBUG) std::cout<<"------ LR observable 5 "<<Obs5<<" calculated ------"<<std::endl;
00099   
00100   // Et-Ratio between the two highest jets in Et and four highest jets
00101   
00102   double Obs6 = (HT-EtSum)/HT;
00103   evtselectVarVal.push_back(std::pair<unsigned int,double>(6,Obs6));
00104   if(DEBUG) std::cout<<"------ LR observable 6 "<<Obs6<<" calculated ------"<<std::endl;
00105   
00106   // Transverse Mass of the system
00107   
00108   TLorentzVector TtbarSystem = (*Hadp)+(*Hadq)+(*Hadb)+(*Lepb)+(*Lept)+(*Lepn);
00109   double MT = TtbarSystem.Mt();
00110   double Obs7 = MT;
00111   evtselectVarVal.push_back(std::pair<unsigned int,double>(7,Obs7));
00112   if(DEBUG) std::cout<<"------ LR observable 7 "<<Obs7<<" calculated ------"<<std::endl;
00113   
00114   // Observables related to the b-disc (jets ordered in Bdisc)
00115   
00116   //sort the TopJets in Bdiscriminant
00117   std::sort(TopJets.begin(),TopJets.end(),BdiscComparator);
00118   
00119   double Obs8;
00120   double Obs9;  
00121   double Obs10;
00122   double Obs11; 
00123   
00124   // Difference in bdisc between the 2nd and the 3rd jets
00125   double BGap = TopJets[1].bDiscriminator(BtagAlgo) - TopJets[2].bDiscriminator(BtagAlgo);
00126   
00127   // Sum of the bdisc of the two highest/lowest jets
00128   double BjetsBdiscSum = TopJets[0].bDiscriminator(BtagAlgo) + TopJets[1].bDiscriminator(BtagAlgo);
00129   double LjetsBdiscSum = TopJets[2].bDiscriminator(BtagAlgo) + TopJets[3].bDiscriminator(BtagAlgo);
00130   
00131   Obs8  = (TopJets[2].bDiscriminator(BtagAlgo) > -10 ? log(BGap) : -5);
00132   if(DEBUG) std::cout<<"------ LR observable 8 "<<Obs8<<" calculated ------"<<std::endl;
00133   Obs9  = (BjetsBdiscSum*BGap); 
00134   if(DEBUG) std::cout<<"------ LR observable 9 "<<Obs9<<" calculated ------"<<std::endl;
00135   Obs10 = (BjetsBdiscSum/LjetsBdiscSum);
00136   if(DEBUG) std::cout<<"------ LR observable 10 "<<Obs10<<" calculated ------"<<std::endl;
00137   // Distance from the origin in the (BjetsBdiscSum, LjetsBdiscSum) plane
00138   Obs11 = 0.707*((BjetsBdiscSum+LjetsBdiscSum)/2 +10);  
00139   if(DEBUG) std::cout<<"------ LR observable 11 "<<Obs11<<" calculated ------"<<std::endl;
00140   
00141   
00142   evtselectVarVal.push_back(std::pair<unsigned int,double>(8,Obs8));
00143   evtselectVarVal.push_back(std::pair<unsigned int,double>(9,Obs9));
00144   evtselectVarVal.push_back(std::pair<unsigned int,double>(10,Obs10));
00145   evtselectVarVal.push_back(std::pair<unsigned int,double>(11,Obs11));
00146   
00147   //sort the TopJets in Et
00148   std::sort(TopJets.begin(),TopJets.end(),EtComparator);
00149   
00150   // Circularity of the event
00151   
00152   double N=0,D=0,C_tmp=0,C=10000;
00153   double N_NoNu=0,D_NoNu=0,C_NoNu_tmp=0,C_NoNu=10000;
00154   double nx,ny,nz;
00155   
00156   // 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))
00157   
00158   for(unsigned int i=0;i<4;i++)
00159     {
00160       D += fabs(TopJets[i].pt());
00161     }
00162   // modified the 26th of September to calculate also the circularity without the neutrino contribution
00163   D += fabs(Lept->Pt());
00164   D_NoNu = D;
00165   D += fabs(Lepn->Pt());
00166   
00167   if((D>0))
00168     {
00169       
00170       // Loop over all the unit vectors in the transverse plane in order to find the miminum : 
00171       for(unsigned int i=0; i<720; i++)
00172         {
00173           
00174           nx = cos((2*PI/720)*i);
00175           ny = sin((2*PI/720)*i);
00176           nz = 0;
00177           N=0;
00178           
00179           for(unsigned int i=0;i<4;i++)
00180             {
00181               N += fabs(TopJets[i].px()*nx+TopJets[i].py()*ny+TopJets[i].pz()*nz);
00182             }
00183           // modified the 26th of September to calculate also the circularity without the neutrino contribution
00184           N += fabs(Lept->Px()*nx+Lept->Py()*ny+Lept->Pz()*nz);
00185           N_NoNu = N;
00186           N += fabs(Lepn->Px()*nx+Lepn->Py()*ny+Lepn->Pz()*nz);
00187           
00188           C_tmp = 2*N/D;
00189           
00190           // modified the 26th of September to calculate also the circularity without the neutrino contribution
00191           C_NoNu_tmp = 2*N_NoNu/D_NoNu;
00192           
00193           if(C_tmp<C) C = C_tmp;
00194           
00195           // modified the 26th of September to calculate also the circularity without the neutrino contribution
00196           if(C_NoNu_tmp<C_NoNu) C_NoNu = C_NoNu_tmp;
00197         }
00198     }
00199   
00200   double Obs12 = ( C!=10000 ? C : 0);
00201   evtselectVarVal.push_back(std::pair<unsigned int,double>(12,Obs12));
00202   if(DEBUG) std::cout<<"------ LR observable 12 "<<Obs12<<" calculated ------"<<std::endl;
00203   
00204   // Circularity of the event without neutrino
00205   
00206   double Obs13 = ( C_NoNu != 10000 ? C_NoNu : 0);
00207   evtselectVarVal.push_back(std::pair<unsigned int,double>(13,Obs13));
00208   if(DEBUG) std::cout<<"------ LR observable 13 "<<Obs13<<" calculated ------"<<std::endl;
00209   
00210   // Centrality of the four highest jets
00211   
00212   double H=0;
00213   for(unsigned int i=0;i<4;i++)
00214     {
00215       H += TopJets[i].energy();
00216     }
00217   
00218   double Obs14 = HT/H;
00219   evtselectVarVal.push_back(std::pair<unsigned int,double>(14,Obs14));
00220   if(DEBUG) std::cout<<"------ LR observable 14 "<<Obs14<<" calculated ------"<<std::endl;
00221   
00222   // Sphericity and Aplanarity without boosting back the system to CM frame
00223   
00224   TMatrixDSym Matrix(3);
00225   
00226   double PX2 = std::pow(Hadp->Px(),2)+std::pow(Hadq->Px(),2)+std::pow(Hadb->Px(),2)+std::pow(Lepb->Px(),2)+std::pow(Lept->Px(),2)+std::pow(Lepn->Px(),2);
00227   double PY2 = std::pow(Hadp->Py(),2)+std::pow(Hadq->Py(),2)+std::pow(Hadb->Py(),2)+std::pow(Lepb->Py(),2)+std::pow(Lept->Py(),2)+std::pow(Lepn->Py(),2);
00228   double PZ2 = std::pow(Hadp->Pz(),2)+std::pow(Hadq->Pz(),2)+std::pow(Hadb->Pz(),2)+std::pow(Lepb->Pz(),2)+std::pow(Lept->Pz(),2)+std::pow(Lepn->Pz(),2);
00229   
00230   double P2  = PX2+PY2+PZ2;
00231   
00232   double PXY = Hadp->Px()*Hadp->Py()+Hadq->Px()*Hadq->Py()+Hadb->Px()*Hadb->Py()+Lepb->Px()*Lepb->Py()+Lept->Px()*Lept->Py()+Lepn->Px()*Lepn->Py();
00233   double PXZ = Hadp->Px()*Hadp->Pz()+Hadq->Px()*Hadq->Pz()+Hadb->Px()*Hadb->Pz()+Lepb->Px()*Lepb->Pz()+Lept->Px()*Lept->Pz()+Lepn->Px()*Lepn->Pz();
00234   double PYZ = Hadp->Py()*Hadp->Pz()+Hadq->Py()*Hadq->Pz()+Hadb->Py()*Hadb->Pz()+Lepb->Py()*Lepb->Pz()+Lept->Py()*Lept->Pz()+Lepn->Py()*Lepn->Pz();
00235   
00236   Matrix(0,0) = PX2/P2; Matrix(0,1) = PXY/P2; Matrix(0,2) = PXZ/P2;
00237   Matrix(1,0) = PXY/P2; Matrix(1,1) = PY2/P2; Matrix(1,2) = PYZ/P2;
00238   Matrix(2,0) = PXZ/P2; Matrix(2,1) = PYZ/P2; Matrix(2,2) = PZ2/P2;
00239   
00240   TMatrixDSymEigen pTensor(Matrix);
00241   
00242   std::vector<double> EigValues;
00243   EigValues.clear();
00244   for(int i=0;i<3;i++)
00245     {
00246       EigValues.push_back(pTensor.GetEigenValues()[i]);
00247     }
00248   
00249   std::sort(EigValues.begin(),EigValues.end(),dComparator);
00250   
00251   double Sphericity = 1.5*(EigValues[1]+EigValues[2]);
00252   double Aplanarity = 1.5*EigValues[2];
00253   
00254   double Obs15 = (isnan(Sphericity) ? 0 : Sphericity);
00255   evtselectVarVal.push_back(std::pair<unsigned int,double>(15,Obs15));
00256   if(DEBUG) std::cout<<"------ LR observable 15 "<<Obs15<<" calculated ------"<<std::endl;
00257   
00258   double Obs16 = (isnan(Aplanarity) ? 0 : Aplanarity);
00259   evtselectVarVal.push_back(std::pair<unsigned int,double>(16,Obs16));
00260   if(DEBUG) std::cout<<"------ LR observable 16 "<<Obs16<<" calculated ------"<<std::endl;
00261   
00262   
00263   // Sphericity and Aplanarity with boosting back the system to CM frame        
00264   
00265   TVector3 BoostBackToCM = -(TtbarSystem.BoostVector());
00266   Hadp->Boost(BoostBackToCM);
00267   Hadq->Boost(BoostBackToCM);
00268   Hadb->Boost(BoostBackToCM);
00269   Lepb->Boost(BoostBackToCM);
00270   Lept->Boost(BoostBackToCM);
00271   Lepn->Boost(BoostBackToCM);
00272   
00273   
00274   double BOOST_PX2 = std::pow(Hadp->Px(),2)+std::pow(Hadq->Px(),2)+std::pow(Hadb->Px(),2)+std::pow(Lepb->Px(),2)+std::pow(Lept->Px(),2)+std::pow(Lepn->Px(),2);
00275   double BOOST_PY2 = std::pow(Hadp->Py(),2)+std::pow(Hadq->Py(),2)+std::pow(Hadb->Py(),2)+std::pow(Lepb->Py(),2)+std::pow(Lept->Py(),2)+std::pow(Lepn->Py(),2);
00276   double BOOST_PZ2 = std::pow(Hadp->Pz(),2)+std::pow(Hadq->Pz(),2)+std::pow(Hadb->Pz(),2)+std::pow(Lepb->Pz(),2)+std::pow(Lept->Pz(),2)+std::pow(Lepn->Pz(),2);
00277   
00278   double BOOST_P2  = BOOST_PX2+BOOST_PY2+BOOST_PZ2;
00279   
00280   double BOOST_PXY = Hadp->Px()*Hadp->Py()+Hadq->Px()*Hadq->Py()+Hadb->Px()*Hadb->Py()+Lepb->Px()*Lepb->Py()+Lept->Px()*Lept->Py()+Lepn->Px()*Lepn->Py();
00281   double BOOST_PXZ = Hadp->Px()*Hadp->Pz()+Hadq->Px()*Hadq->Pz()+Hadb->Px()*Hadb->Pz()+Lepb->Px()*Lepb->Pz()+Lept->Px()*Lept->Pz()+Lepn->Px()*Lepn->Pz();
00282   double BOOST_PYZ = Hadp->Py()*Hadp->Pz()+Hadq->Py()*Hadq->Pz()+Hadb->Py()*Hadb->Pz()+Lepb->Py()*Lepb->Pz()+Lept->Py()*Lept->Pz()+Lepn->Py()*Lepn->Pz();
00283   
00284   TMatrixDSym BOOST_Matrix(3);
00285   
00286   BOOST_Matrix(0,0) = BOOST_PX2/BOOST_P2; BOOST_Matrix(0,1) = BOOST_PXY/BOOST_P2; BOOST_Matrix(0,2) = BOOST_PXZ/BOOST_P2;
00287   BOOST_Matrix(1,0) = BOOST_PXY/BOOST_P2; BOOST_Matrix(1,1) = BOOST_PY2/BOOST_P2; BOOST_Matrix(1,2) = BOOST_PYZ/BOOST_P2;
00288   BOOST_Matrix(2,0) = BOOST_PXZ/BOOST_P2; BOOST_Matrix(2,1) = BOOST_PYZ/BOOST_P2; BOOST_Matrix(2,2) = BOOST_PZ2/BOOST_P2;
00289   
00290   TMatrixDSymEigen BOOST_pTensor(BOOST_Matrix);
00291   
00292   std::vector<double> BOOST_EigValues;
00293   BOOST_EigValues.clear();
00294   for(int i=0;i<3;i++)
00295     {
00296       BOOST_EigValues.push_back(BOOST_pTensor.GetEigenValues()[i]);
00297     }
00298   
00299   std::sort(BOOST_EigValues.begin(),BOOST_EigValues.end(),dComparator);
00300   
00301   double BOOST_Sphericity = 1.5*(BOOST_EigValues[1]+BOOST_EigValues[2]);
00302   double BOOST_Aplanarity = 1.5*BOOST_EigValues[2];
00303   
00304   double Obs17 = ( isnan(BOOST_Sphericity) ? 0 : BOOST_Sphericity );
00305   evtselectVarVal.push_back(std::pair<unsigned int,double>(17,Obs17));
00306   if(DEBUG) std::cout<<"------ LR observable 17 "<<Obs17<<" calculated ------"<<std::endl;
00307   
00308   double Obs18 = ( isnan(BOOST_Aplanarity) ? 0 : BOOST_Aplanarity );
00309   evtselectVarVal.push_back(std::pair<unsigned int,double>(18,Obs18));
00310   if(DEBUG) std::cout<<"------ LR observable 18 "<<Obs18<<" calculated ------"<<std::endl;
00311   
00312   //ratio between ET of the fifth jet and HT
00313   
00314   double Obs19 = (TopJets.size() > 4) ? TopJets[4].et()/HT : 1.0;
00315   evtselectVarVal.push_back(std::pair<unsigned int,double>(19,Obs19));
00316   if(DEBUG) std::cout<<"------ LR observable 19 "<<Obs19<<" calculated ------"<<std::endl;
00317   
00318   // HT variable calculated with all the jets in the event.
00319   
00320   double HT_alljets = 0;
00321   double  H_alljets = 0;
00322   for(unsigned int i=0;i<TopJets.size();i++)
00323     {
00324       HT_alljets += TopJets[i].et();
00325       H_alljets  += TopJets[i].energy();
00326     }
00327   double Obs20 = HT_alljets;
00328   evtselectVarVal.push_back(std::pair<unsigned int,double>(20,Obs20));
00329   if(DEBUG) std::cout<<"------ LR observable 20 "<<Obs20<<" calculated ------"<<std::endl;
00330   
00331   // HT3 = HT calculated with all jets except the two leading jets
00332   
00333   double HT3_alljets = 0;
00334   for(unsigned int i=2;i<TopJets.size();i++)
00335     {
00336       HT3_alljets += TopJets[i].et();
00337     }
00338   double Obs21 = HT3_alljets;
00339   evtselectVarVal.push_back(std::pair<unsigned int,double>(21,Obs21));
00340   if(DEBUG) std::cout<<"------ LR observable 21 "<<Obs21<<" calculated ------"<<std::endl;
00341   
00342   // ET0, ratio of the Et of the leading and HT_alljets
00343   
00344   double ET0 = TopJets[0].et()/HT_alljets;
00345   double Obs22 = ET0;
00346   evtselectVarVal.push_back(std::pair<unsigned int,double>(22,Obs22));
00347   if(DEBUG) std::cout<<"------ LR observable 22 "<<Obs22<<" calculated ------"<<std::endl;
00348   
00349   // Centrality of the event computed with all jets
00350   
00351   double Obs23 = ( (H_alljets>0) ? HT_alljets/H_alljets : 0);
00352   evtselectVarVal.push_back(std::pair<unsigned int,double>(23,Obs23));
00353   if(DEBUG) std::cout<<"------ LR observable 23 "<<Obs23<<" calculated ------"<<std::endl;
00354   
00355   //Fox-Wolfram momenta (1st to 6th), modified for hadron collider and using a Legendre polynomials expansion
00356   
00357   double FW_momentum_0=0, FW_momentum_1=0, FW_momentum_2=0, FW_momentum_3=0, FW_momentum_4=0, FW_momentum_5=0, FW_momentum_6=0;
00358   
00359   for(unsigned int i=0;i<TopJets.size();i++)
00360     {
00361       for(unsigned int j=0;j<TopJets.size();j++)
00362         {
00363           double ET_ij_over_ETSum2= TopJets[i].et()*TopJets[j].et()/(std::pow(HT_alljets,2));
00364           double cosTheta_ij = (TopJets[i].px()*TopJets[j].px()+
00365                                 TopJets[i].py()*TopJets[j].py()+
00366                                 TopJets[i].pz()*TopJets[j].pz())
00367             /(TopJets[i].p4().R()*TopJets[j].p4().R());
00368           FW_momentum_0 += ET_ij_over_ETSum2;
00369           FW_momentum_1 += ET_ij_over_ETSum2 * cosTheta_ij;
00370           FW_momentum_2 += ET_ij_over_ETSum2 * 0.5   * (  3*std::pow(cosTheta_ij,2)- 1);
00371           FW_momentum_3 += ET_ij_over_ETSum2 * 0.5   * (  5*std::pow(cosTheta_ij,3)-  3*cosTheta_ij);
00372           FW_momentum_4 += ET_ij_over_ETSum2 * 0.125 * ( 35*std::pow(cosTheta_ij,4)- 30*std::pow(cosTheta_ij,2)+3);
00373           FW_momentum_5 += ET_ij_over_ETSum2 * 0.125 * ( 63*std::pow(cosTheta_ij,5)- 70*std::pow(cosTheta_ij,3)+15*cosTheta_ij);
00374           FW_momentum_6 += ET_ij_over_ETSum2 * 0.0625* (231*std::pow(cosTheta_ij,6)-315*std::pow(cosTheta_ij,4)+105*std::pow(cosTheta_ij,2)-5);
00375         }
00376     }
00377   
00378   double Obs24 = FW_momentum_0;
00379   evtselectVarVal.push_back(std::pair<unsigned int,double>(24,Obs24));
00380   if(DEBUG) std::cout<<"------ LR observable 24 "<<Obs24<<" calculated ------"<<std::endl;
00381   double Obs25 = FW_momentum_1;
00382   evtselectVarVal.push_back(std::pair<unsigned int,double>(25,Obs25));
00383   if(DEBUG) std::cout<<"------ LR observable 25 "<<Obs25<<" calculated ------"<<std::endl;
00384   double Obs26 = FW_momentum_2;
00385   evtselectVarVal.push_back(std::pair<unsigned int,double>(26,Obs26));
00386   if(DEBUG) std::cout<<"------ LR observable 26 "<<Obs26<<" calculated ------"<<std::endl;
00387   double Obs27 = FW_momentum_3;
00388   evtselectVarVal.push_back(std::pair<unsigned int,double>(27,Obs27));
00389   if(DEBUG) std::cout<<"------ LR observable 27 "<<Obs27<<" calculated ------"<<std::endl;
00390   double Obs28 = FW_momentum_4;
00391   evtselectVarVal.push_back(std::pair<unsigned int,double>(28,Obs28));
00392   if(DEBUG) std::cout<<"------ LR observable 28 "<<Obs28<<" calculated ------"<<std::endl;
00393   double Obs29 = FW_momentum_5;
00394   evtselectVarVal.push_back(std::pair<unsigned int,double>(29,Obs29));
00395   if(DEBUG) std::cout<<"------ LR observable 29 "<<Obs29<<" calculated ------"<<std::endl;
00396   double Obs30 = FW_momentum_6;
00397   evtselectVarVal.push_back(std::pair<unsigned int,double>(30,Obs30));
00398   if(DEBUG) std::cout<<"------ LR observable 30 "<<Obs30<<" calculated ------"<<std::endl;
00399   
00400   // Thrust, thrust major and thrust minor
00401   
00402   TVector3 n(0,0,0), n_tmp(0,0,0);
00403   
00404   double Thrust_D = 0,      Thrust_N = 0;
00405   double Thrust = -1,       Thrust_tmp =0;
00406   double Thrust_Major = -1, Thrust_Major_tmp =0;
00407   double Thrust_Minor = -1;
00408   
00409   for(unsigned int i=0;i<TopJets.size();i++)
00410     {
00411       Thrust_D += TopJets[i].p();
00412     }
00413   
00414   Thrust_D += Lept->P();
00415   
00416   if(Thrust_D>0)
00417     {
00418       // Calculation of the thrust :
00419       for(unsigned int i=0;i<720;i++)
00420         {
00421           for(unsigned int j=0;j<360;j++)
00422             {
00423               n_tmp.SetX(cos((2*PI/720)*i)*sin((PI/360)*j));
00424               n_tmp.SetY(sin((2*PI/720)*i)*sin((PI/360)*j));
00425               n_tmp.SetZ(cos((PI/360)*j));
00426               
00427               for(unsigned int k=0;k<TopJets.size();k++)
00428                 {
00429                   Thrust_N += fabs(TopJets[k].px()*(n_tmp.x())+TopJets[k].py()*(n_tmp.y())+TopJets[k].pz()*(n_tmp.z()));
00430                 }
00431               Thrust_N += fabs(Lept->Px()*(n_tmp.x())+Lept->Py()*(n_tmp.y())+Lept->Pz()*(n_tmp.z()));
00432               
00433               Thrust_tmp = Thrust_N/Thrust_D;
00434               
00435               Thrust_N = 0;
00436               if(Thrust_tmp > Thrust)
00437                 {
00438                   Thrust = Thrust_tmp;
00439                   n.SetXYZ(n_tmp.x(),n_tmp.y(),n_tmp.z());
00440                 }       
00441             }
00442         }
00443       
00444       // Calculation of the thrust major :
00445       TVector3 nT = n.Orthogonal();
00446       nT = nT.Unit();
00447       for(unsigned int i=0;i<720;i++)
00448         {
00449           nT.Rotate((2*PI/720)*i,n);
00450           for(unsigned int j=0;j<TopJets.size();j++)
00451             {
00452               Thrust_N += fabs(TopJets[j].px()*(nT.x())+TopJets[j].py()*(nT.y())+TopJets[j].pz()*(nT.z()));
00453             }
00454           Thrust_N += fabs(Lept->Px()*nT.x()+Lept->Py()*(nT.y())+Lept->Pz()*(nT.z()));
00455           
00456           Thrust_Major_tmp = Thrust_N/Thrust_D;
00457           Thrust_N = 0;
00458           
00459           if(Thrust_Major_tmp > Thrust_Major)
00460             {
00461               Thrust_Major = Thrust_Major_tmp;
00462             }
00463         }
00464       
00465       // Calculation of the thrust minor :
00466       
00467       TVector3 nMinor = nT.Cross(n);
00468       nMinor = nMinor.Unit();
00469       
00470       for(unsigned int i=0;i<TopJets.size();i++)
00471         {
00472           Thrust_N += fabs(TopJets[i].px()*(nMinor.x())+TopJets[i].py()*(nMinor.y())+TopJets[i].pz()*(nMinor.z()));
00473         }
00474       Thrust_N += fabs(Lept->Px()*nMinor.x()+Lept->Py()*(nMinor.y())+Lept->Pz()*(nMinor.z()));
00475       
00476       Thrust_Minor =  Thrust_N/Thrust_D;
00477     }
00478   
00479   double Obs31 = Thrust;
00480   evtselectVarVal.push_back(std::pair<unsigned int,double>(31,Obs31));
00481   if(DEBUG) std::cout<<"------ LR observable 31 "<<Obs31<<" calculated ------"<<std::endl;
00482   double Obs32 = Thrust_Major;
00483   evtselectVarVal.push_back(std::pair<unsigned int,double>(32,Obs32));
00484   if(DEBUG) std::cout<<"------ LR observable 32 "<<Obs32<<" calculated ------"<<std::endl;
00485   double Obs33 = Thrust_Minor;
00486   evtselectVarVal.push_back(std::pair<unsigned int,double>(33,Obs33));
00487   if(DEBUG) std::cout<<"------ LR observable 33 "<<Obs33<<" calculated ------"<<std::endl;
00488   
00489   // Oblateness
00490   
00491   double Obs34 = Thrust_Major-Thrust_Minor;
00492   evtselectVarVal.push_back(std::pair<unsigned int,double>(34,Obs34));
00493   if(DEBUG) std::cout<<"------ LR observable 34 "<<Obs34<<" calculated ------"<<std::endl;
00494   
00495   // Sphericity and Aplanarity without boosting back the system to CM frame and without neutrino
00496   
00497   TMatrixDSym Matrix_NoNu(3);
00498   
00499   double PX2_NoNu = std::pow(Hadp->Px(),2)+std::pow(Hadq->Px(),2)+std::pow(Hadb->Px(),2)+std::pow(Lepb->Px(),2)+std::pow(Lept->Px(),2);
00500   double PY2_NoNu = std::pow(Hadp->Py(),2)+std::pow(Hadq->Py(),2)+std::pow(Hadb->Py(),2)+std::pow(Lepb->Py(),2)+std::pow(Lept->Py(),2);
00501   double PZ2_NoNu = std::pow(Hadp->Pz(),2)+std::pow(Hadq->Pz(),2)+std::pow(Hadb->Pz(),2)+std::pow(Lepb->Pz(),2)+std::pow(Lept->Pz(),2);
00502   
00503   double P2_NoNu  = PX2_NoNu+PY2_NoNu+PZ2_NoNu;
00504   
00505   double PXY_NoNu = Hadp->Px()*Hadp->Py()+Hadq->Px()*Hadq->Py()+Hadb->Px()*Hadb->Py()+Lepb->Px()*Lepb->Py()+Lept->Px()*Lept->Py();
00506   double PXZ_NoNu = Hadp->Px()*Hadp->Pz()+Hadq->Px()*Hadq->Pz()+Hadb->Px()*Hadb->Pz()+Lepb->Px()*Lepb->Pz()+Lept->Px()*Lept->Pz();
00507   double PYZ_NoNu = Hadp->Py()*Hadp->Pz()+Hadq->Py()*Hadq->Pz()+Hadb->Py()*Hadb->Pz()+Lepb->Py()*Lepb->Pz()+Lept->Py()*Lept->Pz();
00508   
00509   Matrix_NoNu(0,0) = PX2_NoNu/P2_NoNu; Matrix_NoNu(0,1) = PXY_NoNu/P2_NoNu; Matrix_NoNu(0,2) = PXZ_NoNu/P2_NoNu;
00510   Matrix_NoNu(1,0) = PXY_NoNu/P2_NoNu; Matrix_NoNu(1,1) = PY2_NoNu/P2_NoNu; Matrix_NoNu(1,2) = PYZ_NoNu/P2_NoNu;
00511   Matrix_NoNu(2,0) = PXZ_NoNu/P2_NoNu; Matrix_NoNu(2,1) = PYZ_NoNu/P2_NoNu; Matrix_NoNu(2,2) = PZ2_NoNu/P2_NoNu;
00512   
00513   TMatrixDSymEigen pTensor_NoNu(Matrix_NoNu);
00514   
00515   std::vector<double> EigValues_NoNu;
00516   EigValues_NoNu.clear();
00517   for(int i=0;i<3;i++)
00518     {
00519       EigValues_NoNu.push_back(pTensor_NoNu.GetEigenValues()[i]);
00520     }
00521   
00522   std::sort(EigValues_NoNu.begin(),EigValues_NoNu.end(),dComparator);
00523   
00524   double Sphericity_NoNu = 1.5*(EigValues_NoNu[1]+EigValues_NoNu[2]);
00525   double Aplanarity_NoNu = 1.5*EigValues_NoNu[2];
00526   
00527   double Obs35 = (isnan(Sphericity_NoNu) ? 0 : Sphericity_NoNu);
00528   evtselectVarVal.push_back(std::pair<unsigned int,double>(35,Obs35));
00529   if(DEBUG) std::cout<<"------ LR observable 35 "<<Obs35<<" calculated ------"<<std::endl;
00530   
00531   double Obs36 = (isnan(Aplanarity_NoNu) ? 0 : Aplanarity_NoNu);
00532   evtselectVarVal.push_back(std::pair<unsigned int,double>(36,Obs36));
00533   if(DEBUG) std::cout<<"------ LR observable 36 "<<Obs36<<" calculated ------"<<std::endl;
00534   
00535   // Sphericity and Aplanarity without boosting back the system to CM frame and without neutrino or lepton
00536   
00537   TMatrixDSym Matrix_NoNuNoLep(3);
00538   
00539   double PX2_NoNuNoLep = std::pow(Hadp->Px(),2)+std::pow(Hadq->Px(),2)+std::pow(Hadb->Px(),2)+std::pow(Lepb->Px(),2);
00540   double PY2_NoNuNoLep = std::pow(Hadp->Py(),2)+std::pow(Hadq->Py(),2)+std::pow(Hadb->Py(),2)+std::pow(Lepb->Py(),2);
00541   double PZ2_NoNuNoLep = std::pow(Hadp->Pz(),2)+std::pow(Hadq->Pz(),2)+std::pow(Hadb->Pz(),2)+std::pow(Lepb->Pz(),2);
00542   
00543   double P2_NoNuNoLep  = PX2_NoNuNoLep+PY2_NoNuNoLep+PZ2_NoNuNoLep;
00544   
00545   double PXY_NoNuNoLep = Hadp->Px()*Hadp->Py()+Hadq->Px()*Hadq->Py()+Hadb->Px()*Hadb->Py()+Lepb->Px()*Lepb->Py();
00546   double PXZ_NoNuNoLep = Hadp->Px()*Hadp->Pz()+Hadq->Px()*Hadq->Pz()+Hadb->Px()*Hadb->Pz()+Lepb->Px()*Lepb->Pz();
00547   double PYZ_NoNuNoLep = Hadp->Py()*Hadp->Pz()+Hadq->Py()*Hadq->Pz()+Hadb->Py()*Hadb->Pz()+Lepb->Py()*Lepb->Pz();
00548   
00549   Matrix_NoNuNoLep(0,0) = PX2_NoNuNoLep/P2_NoNuNoLep; Matrix_NoNuNoLep(0,1) = PXY_NoNuNoLep/P2_NoNuNoLep; Matrix_NoNuNoLep(0,2) = PXZ_NoNuNoLep/P2_NoNuNoLep;
00550   Matrix_NoNuNoLep(1,0) = PXY_NoNuNoLep/P2_NoNuNoLep; Matrix_NoNuNoLep(1,1) = PY2_NoNuNoLep/P2_NoNuNoLep; Matrix_NoNuNoLep(1,2) = PYZ_NoNuNoLep/P2_NoNuNoLep;
00551   Matrix_NoNuNoLep(2,0) = PXZ_NoNuNoLep/P2_NoNuNoLep; Matrix_NoNuNoLep(2,1) = PYZ_NoNuNoLep/P2_NoNuNoLep; Matrix_NoNuNoLep(2,2) = PZ2_NoNuNoLep/P2_NoNuNoLep;
00552   
00553   TMatrixDSymEigen pTensor_NoNuNoLep(Matrix_NoNuNoLep);
00554   
00555   std::vector<double> EigValues_NoNuNoLep;
00556   EigValues_NoNuNoLep.clear();
00557   for(int i=0;i<3;i++)
00558     {
00559       EigValues_NoNuNoLep.push_back(pTensor_NoNuNoLep.GetEigenValues()[i]);
00560     }
00561   
00562   std::sort(EigValues_NoNuNoLep.begin(),EigValues_NoNuNoLep.end(),dComparator);
00563   
00564   double Sphericity_NoNuNoLep = 1.5*(EigValues_NoNuNoLep[1]+EigValues_NoNuNoLep[2]);
00565   double Aplanarity_NoNuNoLep = 1.5*EigValues_NoNuNoLep[2];
00566   
00567   double Obs37 = (isnan(Sphericity_NoNuNoLep) ? 0 : Sphericity_NoNuNoLep);
00568   evtselectVarVal.push_back(std::pair<unsigned int,double>(37,Obs37));
00569   if(DEBUG) std::cout<<"------ LR observable 37 "<<Obs37<<" calculated ------"<<std::endl;
00570   
00571   double Obs38 = (isnan(Aplanarity_NoNuNoLep) ? 0 : Aplanarity_NoNuNoLep);
00572   evtselectVarVal.push_back(std::pair<unsigned int,double>(38,Obs38));
00573   if(DEBUG) std::cout<<"------ LR observable 38 "<<Obs38<<" calculated ------"<<std::endl;
00574   
00575   
00576   // Put the vector in the TtSemiEvtSolution
00577   TS.setLRSignalEvtObservables(evtselectVarVal);
00578   if(DEBUG) std::cout<<"------  Observable values stored in the event  ------"<<std::endl;
00579   
00580   delete Hadp;
00581   if(DEBUG) std::cout<<"------     Pointer to Hadp deleted             ------"<<std::endl;
00582   delete Hadq;
00583   if(DEBUG) std::cout<<"------     Pointer to Hadq deleted             ------"<<std::endl;
00584   delete Hadb;
00585   if(DEBUG) std::cout<<"------     Pointer to Hadb deleted             ------"<<std::endl;
00586   delete Lepb;
00587   if(DEBUG) std::cout<<"------     Pointer to Lepb deleted      ------"<<std::endl;
00588   delete Lept;
00589   if(DEBUG) std::cout<<"------     Pointer to Lepn deleted      ------"<<std::endl;
00590   delete Lepn;
00591   if(DEBUG) std::cout<<"------     Pointer to Mez deleted       ------"<<std::endl;
00592   delete Mez;
00593   
00594   if(DEBUG) std::cout<<"------------ LR observables calculated -----------"<<std::endl;  
00595 }
/data/refman/pasoursint/CMSSW_5_3_10_patch2/src/TopQuarkAnalysis/TopEventSelection/src/TtSemiLRSignalSelObservables.cc
