Reference Manual

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