CMSSW Reference Manual

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