CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/TopQuarkAnalysis/TopJetCombination/src/TtSemiLRJetCombObservables.cc

Go to the documentation of this file.
00001 //
00002 // Author:  Jan Heyninck
00003 // Created: Tue Apr  3 17:33:23 PDT 2007
00004 //
00005 // $Id: TtSemiLRJetCombObservables.cc,v 1.15 2010/03/25 09:22:18 snaumann Exp $
00006 //
00007 #include "TopQuarkAnalysis/TopJetCombination/interface/TtSemiLRJetCombObservables.h"
00008 #include "DataFormats/PatCandidates/interface/Jet.h"
00009 #include "DataFormats/Math/interface/deltaR.h"
00010 #include "DataFormats/Math/interface/deltaPhi.h"
00011 #include "DataFormats/Common/interface/Handle.h"
00012 #include "AnalysisDataFormats/TopObjects/interface/TtGenEvent.h"
00013 
00014 // constructor with path; default should not be used
00015 TtSemiLRJetCombObservables::TtSemiLRJetCombObservables() {}
00016 
00017 
00018 // destructor
00019 TtSemiLRJetCombObservables::~TtSemiLRJetCombObservables() {}
00020 
00021 std::vector< TtSemiLRJetCombObservables::IntBoolPair >
00022 TtSemiLRJetCombObservables::operator() (TtSemiEvtSolution &solution, const edm::Event & iEvent, bool matchOnly)
00023 { 
00024   bool debug=false;
00025   
00026   evtselectVarVal.clear();
00027   evtselectVarMatch.clear();
00028  
00029   // Check whether the objects are matched:
00030   bool matchHadt = false;
00031   bool matchLept = false;
00032   bool matchHadW = false;
00033   bool matchLepW = false;
00034   bool matchHadb = false;
00035   bool matchLepb = false;
00036   bool matchHadp = false;
00037   bool matchHadq = false;  
00038   bool matchHadpq = false;
00039   bool matchHadqp = false;
00040   bool matchLepl = false;
00041   bool matchLepn = false;
00042   
00043   if(debug) std::cout << "== start matching the objects " << std::endl;
00044   
00045   try {
00046     
00047     if(debug) std::cout << "== start trying " << std::endl;
00048     double drLepl=0, drLepn=0, drHadb=0, drLepb=0, drHadp=0, drHadq=0, drHadpq=0, drHadqp=0, drHadt=0, drLept=0, drLepW=0, drHadW=0;
00049     
00050     edm::Handle<TtGenEvent> genEvent;
00051     iEvent.getByLabel ("genEvt",genEvent);
00052     if(debug) std::cout << "== found genEvent " << std::endl;
00053     
00054     if (genEvent->isSemiLeptonic() && genEvent->numberOfBQuarks() == 2) {
00055 
00056       //if(debug) std::cout << "== genEvent->quarkFromAntiTop() " << genEvent->quarkFromAntiTop()->pt() << std::endl;
00057       if(debug) std::cout << "== genEvent->isSemiLeptonic() && genEvent->numberOfBQuarks() == 2 " << std::endl;
00058       if(debug) std::cout << "== solution.getDecay() == " <<solution.getDecay()<< std::endl;     
00059       if(debug) std::cout << "== solution.getRecLepm().pt() = " <<solution.getRecLepm().pt()  << std::endl;
00060       //if(debug) if(solution.getGenLepl() == 0) std::cout << "solution.getGenLepl() == NULL" << std::endl;
00061       if(debug) std::cout << "== *(solution.getGenLept())" << solution.getGenLept()->pt() << std::endl;
00062       if(debug) std::cout << "== *(solution.getGenLepl())" << solution.getGenLepl()->pt() << std::endl;
00063       // std::cout << "Semilepton:\n";
00064       // Match the lepton by deltaR
00065       if (solution.getDecay() == "muon")     drLepl = DeltaR<pat::Particle, reco::GenParticle>()(solution.getRecLepm(), *(solution.getGenLepl()));
00066       if(debug) std::cout << "== matching lepton " << std::endl;
00067       if (solution.getDecay() == "electron") drLepl = DeltaR<pat::Particle, reco::GenParticle>()(solution.getRecLepe(), *(solution.getGenLepl()));
00068       matchLepl = (drLepl < 0.3);
00069       
00070       if(debug) std::cout << "== lepton is matched " << std::endl;
00071       // Match the neutrino by deltaR
00072       drLepn = DeltaR<pat::Particle, reco::GenParticle>()(solution.getRecLepn(), *(solution.getGenLepn()));
00073       matchLepn = (drLepn < 0.3);
00074 
00075       // Match the hadronic b by deltaR
00076       drHadb = DeltaR<pat::Jet, reco::GenParticle>()(solution.getRecHadb(), *(solution.getGenHadb()));
00077       matchHadb = (drHadb < 0.3);
00078 
00079       // Match the hadronicleptonic b by deltaR
00080       drLepb = DeltaR<pat::Jet, reco::GenParticle>()(solution.getRecLepb(), *(solution.getGenLepb()));
00081       matchLepb = (drLepb < 0.3);
00082 
00083       // Match the hadronic p by deltaR
00084       drHadp = DeltaR<pat::Jet, reco::GenParticle>()(solution.getRecHadp(), *(solution.getGenHadp()));
00085       matchHadp = (drHadp < 0.3);
00086       
00087       // Match the hadronic pq by deltaR
00088       drHadpq = DeltaR<pat::Jet, reco::GenParticle>()(solution.getRecHadp(), *(solution.getGenHadq()));
00089       matchHadpq = (drHadpq < 0.3);
00090     
00091       // Match the hadronic q by deltaR
00092       drHadq = DeltaR<pat::Jet, reco::GenParticle>()(solution.getRecHadq(), *(solution.getGenHadq()));
00093       matchHadq = (drHadq < 0.3);      
00094 
00095       // Match the hadronic qp by deltaR
00096       drHadqp = DeltaR<pat::Jet, reco::GenParticle>()(solution.getRecHadq(), *(solution.getGenHadp()));
00097       matchHadqp = (drHadqp < 0.3);  
00098 
00099       // Match the hadronic W by deltaR
00100       drHadW = DeltaR<reco::Particle, reco::GenParticle>()(solution.getRecHadW(), *(solution.getGenHadW()));
00101       matchHadW = (drHadW < 0.3);    
00102 
00103       // Match the leptonic W by deltaR
00104       drLepW = DeltaR<reco::Particle, reco::GenParticle>()(solution.getRecLepW(), *(solution.getGenLepW()));
00105       matchLepW = (drLepW < 0.3);  
00106      
00107       // Match the hadronic t by deltaR
00108       drHadt = DeltaR<reco::Particle, reco::GenParticle>()(solution.getRecHadt(), *(solution.getGenHadt()));
00109       matchHadt = (drHadt < 0.3);    
00110 
00111       // Match the leptonic t by deltaR
00112       drLept = DeltaR<reco::Particle, reco::GenParticle>()(solution.getRecLept(), *(solution.getGenLept()));
00113       matchLept = (drLept < 0.3);   
00114     }
00115   } catch (...){std::cout << "Exception\n";}
00116 
00117   if(debug) std::cout << "== objects matched" <<std::endl;
00118 
00119   edm::Handle<std::vector<pat::Jet> > jets;
00120   iEvent.getByLabel(jetSource_, jets);
00121 
00122   if(debug) std::cout << "== start calculating observables" << std::endl;
00123 
00124 
00125   //obs1 : pt(had top) 
00126   double AverageTop =((solution.getHadb().p4()+solution.getHadq().p4()+solution.getHadp().p4()).pt()+(solution.getLepb().p4()+solution.getHadq().p4()+solution.getHadp().p4()).pt()+(solution.getHadb().p4()+solution.getLepb().p4()+solution.getHadp().p4()).pt()+(solution.getHadb().p4()+solution.getHadq().p4()+solution.getLepb().p4()).pt())/4.;
00127   double Obs1 = ((solution.getHadb().p4()+solution.getHadq().p4()+solution.getHadp().p4()).pt())/AverageTop;
00128   evtselectVarVal.push_back(IntDblPair(1,Obs1));
00129   evtselectVarMatch.push_back(IntBoolPair(1, ((matchHadq&&matchHadp)||(matchHadpq&&matchHadqp))&&matchHadb)); 
00130 
00131   if(debug) std::cout << "== observable 1 " << Obs1 << std::endl;
00132  
00133   //obs2 : (pt_b1 + pt_b2)/(sum jetpt)
00134   double Obs2 = (solution.getHadb().pt()+solution.getLepb().pt())/(solution.getHadp().pt()+solution.getHadq().pt());
00135   evtselectVarVal.push_back(IntDblPair(2,Obs2));
00136   evtselectVarMatch.push_back(IntBoolPair(2,((matchHadp&&matchHadq)||(matchHadpq&&matchHadqp))&&matchHadb&&matchLepb)); 
00137 
00138   if(debug) std::cout << "== observable 2 " << Obs2 << std::endl;
00139 
00140   //obs3: delta R between lep b and lepton 
00141   double Obs3 = -999;
00142   if (solution.getDecay() == "muon")     Obs3 = ROOT::Math::VectorUtil::DeltaR( solution.getLepb().p4(),solution.getRecLepm().p4() );
00143   if (solution.getDecay() == "electron") Obs3 = ROOT::Math::VectorUtil::DeltaR( solution.getLepb().p4(),solution.getRecLepe().p4() );
00144   evtselectVarVal.push_back(IntDblPair(3,Obs3));
00145   evtselectVarMatch.push_back(IntBoolPair(3,matchLepb&&matchLepl)); 
00146 
00147   if(debug) std::cout << "== observable 3 " << Obs3 << std::endl;
00148   
00149    //obs4 : del R ( had b, had W)
00150   double Obs4 = ROOT::Math::VectorUtil::DeltaR( solution.getHadb().p4(), solution.getHadq().p4()+solution.getHadp().p4() );
00151   evtselectVarVal.push_back(IntDblPair(4,Obs4));  
00152   evtselectVarMatch.push_back(IntBoolPair(4,matchHadb&&((matchHadp&&matchHadp)||(matchHadpq&&matchHadqp)))); 
00153 
00154   if(debug) std::cout << "== observable 4 " << Obs4 << std::endl;
00155  
00156   //obs5 : del R between light quarkssolution.getHadp().phi(
00157   double Obs5 = ROOT::Math::VectorUtil::DeltaR( solution.getHadq().p4(),solution.getHadp().p4() );
00158   evtselectVarVal.push_back(IntDblPair(5,Obs5)); 
00159   evtselectVarMatch.push_back(IntBoolPair(5,(matchHadp&&matchHadq)||(matchHadpq&&matchHadqp))); 
00160 
00161   if(debug) std::cout << "== observable 5 " << Obs5 << std::endl;
00162 
00163   //obs6 : b-tagging information
00164   double Obs6 = 0;
00165   if ( fabs(solution.getHadb().bDiscriminator("trackCountingJetTags") +10) < 0.0001 || fabs(solution.getLepb().bDiscriminator("trackCountingJetTags") +10)< 0.0001 ){
00166     Obs6 = -10.;
00167   } else {
00168     Obs6 = (solution.getHadb().bDiscriminator("trackCountingJetTags")+solution.getLepb().bDiscriminator("trackCountingJetTags"));
00169   }
00170   evtselectVarVal.push_back(IntDblPair(6,Obs6)); 
00171   evtselectVarMatch.push_back(IntBoolPair(6,1)); 
00172  
00173   if(debug) std::cout << "== observable 6 " << Obs6 << std::endl;
00174 
00175   //obs7 : chi2 value of kinematical fit with W-mass constraint
00176   double Obs7 =0;
00177   if(solution.getProbChi2() <0){Obs7 = -0;} else { Obs7 = log10(solution.getProbChi2()+.00001);}
00178   evtselectVarVal.push_back(IntDblPair(7,Obs7)); 
00179   evtselectVarMatch.push_back(IntBoolPair(7,((matchHadp&&matchHadq)||(matchHadpq&&matchHadqp)))); 
00180 
00181   if(debug) std::cout << "== observable 7 " << Obs7 << std::endl;
00182 
00183   //obs8(=7+1)
00184   double Obs8 =  solution.getCalHadt().p4().pt();
00185   evtselectVarVal.push_back(IntDblPair(8,Obs8));
00186   evtselectVarMatch.push_back(IntBoolPair(8, matchHadb&&((matchHadp&&matchHadq)||(matchHadpq&&matchHadqp))));
00187 
00188   if(debug) std::cout << "== observable 8 " << Obs8 << std::endl;
00189 
00190   //obs9
00191   double Obs9  = fabs(solution.getCalHadt().p4().eta());
00192   evtselectVarVal.push_back(IntDblPair(9,Obs9));
00193   evtselectVarMatch.push_back(IntBoolPair(9, matchHadb&&((matchHadp&&matchHadq)||(matchHadpq&&matchHadqp))));
00194 
00195   if(debug) std::cout << "== observable 9 " << Obs9 << std::endl;
00196 
00197   //obs10  
00198   double Obs10  = solution.getCalHadt().p4().theta();
00199   evtselectVarVal.push_back(IntDblPair(10,Obs10));
00200   evtselectVarMatch.push_back(IntBoolPair(10, matchHadb&&((matchHadp&&matchHadq)||(matchHadpq&&matchHadqp))));
00201 
00202   if(debug) std::cout << "== observable 10 " << Obs10 << std::endl;
00203  
00204   //obs11
00205   double Obs11  = solution.getCalHadW().p4().pt();
00206   evtselectVarVal.push_back(IntDblPair(11,Obs11));
00207   evtselectVarMatch.push_back(IntBoolPair(11, (matchHadp&&matchHadq)||(matchHadpq&&matchHadqp)));
00208   
00209   if(debug) std::cout << "== observable 11 " << Obs11 << std::endl;
00210 
00211   //obs12
00212   double Obs12  = fabs(solution.getCalHadW().p4().eta());
00213   evtselectVarVal.push_back(IntDblPair(12,Obs12));
00214   evtselectVarMatch.push_back(IntBoolPair(12, (matchHadp&&matchHadq)||(matchHadpq&&matchHadqp)));
00215 
00216   if(debug) std::cout << "== observable 12 " << Obs12 << std::endl;
00217   
00218   //obs13
00219   double Obs13  = solution.getCalHadW().p4().theta();
00220   evtselectVarVal.push_back(IntDblPair(13,Obs13));
00221   evtselectVarMatch.push_back(IntBoolPair(13, (matchHadp&&matchHadq)||(matchHadpq&&matchHadqp)));
00222 
00223   if(debug) std::cout << "== observable 13 " << Obs13 << std::endl;
00224         
00225   //obs14
00226   double Obs14  = solution.getCalHadb().p4().pt();
00227   evtselectVarVal.push_back(IntDblPair(14,Obs14));
00228   evtselectVarMatch.push_back(IntBoolPair(14, matchHadb));
00229 
00230   if(debug) std::cout << "== observable 14 " << Obs14 << std::endl;
00231             
00232   //obs15
00233   double Obs15  = fabs(solution.getCalHadb().p4().eta());
00234   evtselectVarVal.push_back(IntDblPair(15,Obs15));
00235   evtselectVarMatch.push_back(IntBoolPair(15, matchHadb));
00236 
00237   if(debug) std::cout << "== observable 15 " << Obs15 << std::endl;
00238             
00239   //obs16
00240   double Obs16  = solution.getCalHadb().p4().theta();
00241   evtselectVarVal.push_back(IntDblPair(16,Obs16));
00242   evtselectVarMatch.push_back(IntBoolPair(16, matchHadb));
00243 
00244   if(debug) std::cout << "== observable 16 " << Obs16 << std::endl;
00245             
00246   //obs17
00247   double Obs17  = solution.getCalHadp().p4().pt();
00248   evtselectVarVal.push_back(IntDblPair(17,Obs17));
00249   evtselectVarMatch.push_back(IntBoolPair(17, matchHadp));
00250 
00251   if(debug) std::cout << "== observable 17 " << Obs17 << std::endl;
00252             
00253   //obs18
00254   double Obs18  = fabs(solution.getCalHadp().p4().eta());
00255   evtselectVarVal.push_back(IntDblPair(18,Obs18));
00256   evtselectVarMatch.push_back(IntBoolPair(18, matchHadp));
00257 
00258   if(debug) std::cout << "== observable 18 " << Obs18 << std::endl;
00259             
00260   //obs19
00261   double Obs19  = solution.getCalHadp().p4().theta();
00262   evtselectVarVal.push_back(IntDblPair(19,Obs19));
00263   evtselectVarMatch.push_back(IntBoolPair(19, matchHadp));
00264 
00265   if(debug) std::cout << "== observable 19 " << Obs19 << std::endl;
00266             
00267   //obs20
00268   double Obs20  = solution.getCalHadq().p4().pt();
00269   evtselectVarVal.push_back(IntDblPair(20,Obs20));
00270   evtselectVarMatch.push_back(IntBoolPair(20, matchHadq));
00271 
00272   if(debug) std::cout << "== observable 20 " << Obs20 << std::endl;
00273             
00274   //obs21
00275   double Obs21  = fabs(solution.getCalHadq().p4().eta());
00276   evtselectVarVal.push_back(IntDblPair(21,Obs21));
00277   evtselectVarMatch.push_back(IntBoolPair(21, matchHadq));
00278 
00279   if(debug) std::cout << "== observable 21 " << Obs21 << std::endl;
00280             
00281   //obs22
00282   double Obs22  = solution.getCalHadq().p4().theta();
00283   evtselectVarVal.push_back(IntDblPair(22,Obs22));
00284   evtselectVarMatch.push_back(IntBoolPair(22, matchHadq));
00285 
00286   if(debug) std::cout << "== observable 22 " << Obs22 << std::endl;
00287             
00288   //obs23
00289   double Obs23  = solution.getCalLept().p4().pt();
00290   evtselectVarVal.push_back(IntDblPair(23,Obs23));
00291   evtselectVarMatch.push_back(IntBoolPair(23, matchLepl&&matchLepn&&matchLepb));
00292 
00293   if(debug) std::cout << "== observable 23 " << Obs23 << std::endl;
00294             
00295   //obs24
00296   double Obs24  = fabs(solution.getCalLept().p4().eta());
00297   evtselectVarVal.push_back(IntDblPair(24,Obs24));
00298   evtselectVarMatch.push_back(IntBoolPair(24, matchLepl&&matchLepn&&matchLepb));
00299 
00300   if(debug) std::cout << "== observable 24 " << Obs24 << std::endl;
00301             
00302   //obs25
00303   double Obs25  = solution.getCalLept().p4().theta();
00304   evtselectVarVal.push_back(IntDblPair(25,Obs25));
00305   evtselectVarMatch.push_back(IntBoolPair(25, matchLepl&&matchLepn&&matchLepb));
00306 
00307   if(debug) std::cout << "== observable 25 " << Obs25 << std::endl;
00308             
00309   //obs26
00310   double Obs26  = solution.getRecLepW().p4().pt(); 
00311   evtselectVarVal.push_back(IntDblPair(26,Obs26));
00312   evtselectVarMatch.push_back(IntBoolPair(26, matchLepl&&matchLepn));
00313 
00314   if(debug) std::cout << "== observable 26 " << Obs26 << std::endl;
00315 
00316   //obs27
00317   double Obs27  = fabs(solution.getRecLepW().p4().eta());  
00318   evtselectVarVal.push_back(IntDblPair(27,Obs27));
00319   evtselectVarMatch.push_back(IntBoolPair(27, matchLepl&&matchLepn));
00320 
00321   if(debug) std::cout << "== observable 27 " << Obs27 << std::endl;
00322             
00323   //obs28
00324   double Obs28  = solution.getRecLepW().p4().theta(); 
00325   evtselectVarVal.push_back(IntDblPair(28,Obs28));
00326   evtselectVarMatch.push_back(IntBoolPair(28, matchLepl&&matchLepn));
00327 
00328   if(debug) std::cout << "== observable 28 " << Obs28 << std::endl;
00329   
00330   //obs29 
00331   double Obs29  = solution.getCalLepb().p4().pt();
00332   evtselectVarVal.push_back(IntDblPair(29,Obs29));
00333   evtselectVarMatch.push_back(IntBoolPair(29, matchLepb));
00334 
00335   if(debug) std::cout << "== observable 29 " << Obs29 << std::endl;
00336   
00337   //obs30
00338   double Obs30  = fabs(solution.getCalLepb().p4().eta());
00339   evtselectVarVal.push_back(IntDblPair(30,Obs30));
00340   evtselectVarMatch.push_back(IntBoolPair(30, matchLepb));
00341 
00342   if(debug) std::cout << "== observable 30 " << Obs30 << std::endl;
00343   
00344   //obs31
00345   double Obs31  = solution.getCalLepb().p4().theta();
00346   evtselectVarVal.push_back(IntDblPair(31,Obs31));
00347   evtselectVarMatch.push_back(IntBoolPair(31, matchLepb));
00348 
00349   if(debug) std::cout << "== observable 31 " << Obs31 << std::endl;
00350   
00351   //obs32
00352   double Obs32 = -999.;
00353   if (solution.getDecay() == "muon") Obs32 = solution.getRecLepm().p4().pt(); 
00354   if (solution.getDecay() == "electron") Obs32 = solution.getRecLepe().p4().pt(); 
00355   evtselectVarVal.push_back(IntDblPair(32,Obs32));
00356   evtselectVarMatch.push_back(IntBoolPair(32, matchLepl));
00357 
00358   if(debug) std::cout << "== observable 32 " << Obs32 << std::endl;
00359    
00360   //obs33
00361   double Obs33 = -999.;
00362   if (solution.getDecay() == "muon") Obs33 = fabs(solution.getRecLepm().p4().eta());
00363   if (solution.getDecay() == "electron") Obs33 = fabs(solution.getRecLepe().p4().eta());
00364   evtselectVarVal.push_back(IntDblPair(33,Obs33));
00365   evtselectVarMatch.push_back(IntBoolPair(33, matchLepl));
00366 
00367   if(debug) std::cout << "== observable 33 " << Obs33 << std::endl;
00368  
00369   //obs34
00370   double Obs34 = -999.;
00371   if (solution.getDecay() == "muon") Obs34 = fabs(solution.getRecLepm().p4().theta());
00372   if (solution.getDecay() == "electron") Obs34 = fabs(solution.getRecLepe().p4().theta()); 
00373   evtselectVarVal.push_back(IntDblPair(34,Obs34));
00374   evtselectVarMatch.push_back(IntBoolPair(34, matchLepl));
00375 
00376   if(debug) std::cout << "== observable 34 " << Obs34 << std::endl;
00377   
00378   //obs35
00379   double Obs35  = solution.getFitLepn().p4().pt();
00380   evtselectVarVal.push_back(IntDblPair(35,Obs35));
00381   evtselectVarMatch.push_back(IntBoolPair(35, matchLepn));
00382 
00383   if(debug) std::cout << "== observable 35 " << Obs35 << std::endl;
00384   
00385   //obs36
00386   double Obs36  = fabs(solution.getFitLepn().p4().eta());
00387   evtselectVarVal.push_back(IntDblPair(36,Obs36));
00388   evtselectVarMatch.push_back(IntBoolPair(36, matchLepn));
00389 
00390   if(debug) std::cout << "== observable 36 " << Obs36 << std::endl;
00391   
00392   //obs37
00393   double Obs37  = solution.getFitLepn().p4().theta();
00394   evtselectVarVal.push_back(IntDblPair(37,Obs37));
00395   evtselectVarMatch.push_back(IntBoolPair(37, matchLepn));
00396 
00397   if(debug) std::cout << "== observable 37 " << Obs37 << std::endl;
00398   
00399   // 2 particle kinematics
00400   //obs38 
00401   double Obs38  = fabs(solution.getCalHadW().p4().phi()- solution.getRecLepW().p4().phi());
00402   if (Obs38 > 3.1415927)  Obs38 =  2*3.1415927 - Obs31;
00403   if (Obs38 < -3.1415927) Obs38 = -2*3.1415927 - Obs31;
00404   evtselectVarVal.push_back(IntDblPair(38,Obs38));
00405   evtselectVarMatch.push_back(IntBoolPair(38, matchLepl&&matchLepn&&((matchHadp&&matchHadq)||(matchHadpq&&matchHadqp))));
00406 
00407   if(debug) std::cout << "== observable 38 " << Obs38 << std::endl;
00408   
00409   //obs39
00410   double Obs39  = fabs(solution.getCalHadW().p4().eta()- solution.getRecLepW().p4().eta());
00411   evtselectVarVal.push_back(IntDblPair(39,Obs39));
00412   evtselectVarMatch.push_back(IntBoolPair(39, matchLepl&&matchLepn&&((matchHadp&&matchHadq)||(matchHadpq&&matchHadqp))));
00413 
00414   if(debug) std::cout << "== observable 39 " << Obs39 << std::endl;
00415   
00416   //obs40
00417   double Obs40  = fabs(solution.getCalHadW().p4().theta()- solution.getRecLepW().p4().theta());
00418   evtselectVarVal.push_back(IntDblPair(40,Obs40));
00419   evtselectVarMatch.push_back(IntBoolPair(40, matchLepl&&matchLepn&&((matchHadp&&matchHadq)||(matchHadpq&&matchHadqp))));
00420 
00421   if(debug) std::cout << "== observable 40 " << Obs40 << std::endl;
00422   
00423   //obs41
00424   double Obs41  = ROOT::Math::VectorUtil::DeltaR(solution.getCalHadW().p4(), solution.getRecLepW().p4());
00425   evtselectVarVal.push_back(IntDblPair(41,Obs41));
00426   evtselectVarMatch.push_back(IntBoolPair(41, matchLepl&&matchLepn&&((matchHadp&&matchHadq)||(matchHadpq&&matchHadqp))));
00427 
00428   if(debug) std::cout << "== observable 41 " << Obs41 << std::endl;
00429   
00430   //obs42
00431   double Obs42  = fabs(solution.getCalHadb().p4().phi()- solution.getCalLepb().p4().phi());
00432   if (Obs42 > 3.1415927)  Obs42 =  2*3.1415927 - Obs42;
00433   if (Obs42 < -3.1415927) Obs42 = -2*3.1415927 - Obs42;
00434   evtselectVarVal.push_back(IntDblPair(42,Obs42));
00435   evtselectVarMatch.push_back(IntBoolPair(42, matchHadb&&matchLepb));
00436 
00437   if(debug) std::cout << "== observable 42 " << Obs42 << std::endl;
00438   
00439   //obs43
00440   double Obs43  = fabs(solution.getCalHadb().p4().eta()- solution.getCalLepb().p4().eta());
00441   evtselectVarVal.push_back(IntDblPair(43,Obs43));
00442   evtselectVarMatch.push_back(IntBoolPair(43, matchHadb&&matchLepb));
00443 
00444   if(debug) std::cout << "== observable 43 " << Obs43 << std::endl;
00445   
00446   //obs44
00447   double Obs44 = fabs(solution.getCalHadb().p4().theta()- solution.getCalLepb().p4().theta());
00448   evtselectVarVal.push_back(IntDblPair(44,Obs44));
00449   evtselectVarMatch.push_back(IntBoolPair(44, matchHadb&&matchLepb));
00450 
00451   if(debug) std::cout << "== observable 44 " << Obs44 << std::endl;
00452   
00453   //obs45
00454   double Obs45  = ROOT::Math::VectorUtil::DeltaR(solution.getCalHadb().p4(), solution.getCalLepb().p4());
00455   evtselectVarVal.push_back(IntDblPair(45,Obs45));
00456   evtselectVarMatch.push_back(IntBoolPair(45, matchHadb&&matchLepb));
00457 
00458   if(debug) std::cout << "== observable 45 " << Obs45 << std::endl;
00459   
00460   //obs46
00461   double Obs46  = fabs(solution.getCalHadb().p4().phi()- solution.getCalHadW().p4().phi());
00462   if (Obs46 > 3.1415927)  Obs46 =  2*3.1415927 - Obs46;
00463   if (Obs46 < -3.1415927) Obs46 = -2*3.1415927 - Obs46;
00464   evtselectVarVal.push_back(IntDblPair(46,Obs46));
00465   evtselectVarMatch.push_back(IntBoolPair(46, matchHadb&&((matchHadq&&matchHadp)||(matchHadpq&&matchHadqp))));
00466 
00467   if(debug) std::cout << "== observable 46 " << Obs46 << std::endl;
00468 
00469   //obs47
00470   double Obs47  = fabs(solution.getCalHadb().p4().eta()- solution.getCalHadW().p4().eta());
00471   evtselectVarVal.push_back(IntDblPair(47,Obs47));
00472   evtselectVarMatch.push_back(IntBoolPair(47, matchHadb&&((matchHadp&&matchHadq)||(matchHadpq&&matchHadqp))));
00473 
00474   if(debug) std::cout << "== observable 47 " << Obs47 << std::endl;
00475   
00476   //obs48
00477   double Obs48  = fabs(solution.getCalHadb().p4().theta()- solution.getCalHadW().p4().theta());
00478   evtselectVarVal.push_back(IntDblPair(48,Obs48));
00479   evtselectVarMatch.push_back(IntBoolPair(48, matchHadb&&((matchHadp&&matchHadq)||(matchHadpq&&matchHadqp))));
00480 
00481   if(debug) std::cout << "== observable 48 " << Obs48 << std::endl;
00482   
00483   //obs49
00484   double Obs49  = ROOT::Math::VectorUtil::DeltaR(solution.getCalHadb().p4(), solution.getCalHadW().p4());
00485   evtselectVarVal.push_back(IntDblPair(49,Obs49));
00486   evtselectVarMatch.push_back(IntBoolPair(49, matchHadb&&((matchHadp&&matchHadq)||(matchHadpq&&matchHadqp))));
00487 
00488   if(debug) std::cout << "== observable 49 " << Obs49 << std::endl;
00489   
00490   //obs50
00491   double Obs50  = fabs(solution.getCalLepb().p4().phi()- solution.getRecLepW().p4().phi());
00492   if (Obs50 > 3.1415927)  Obs50 =  2*3.1415927 - Obs50;
00493   if (Obs50 < -3.1415927) Obs50 = -2*3.1415927 - Obs50;
00494   evtselectVarVal.push_back(IntDblPair(50,Obs50));
00495   evtselectVarMatch.push_back(IntBoolPair(50, matchLepb&&matchLepn&&matchLepl));
00496 
00497   if(debug) std::cout << "== observable 50 " << Obs50 << std::endl;
00498   
00499   //obs51
00500   double Obs51  = fabs(solution.getCalLepb().p4().eta()- solution.getRecLepW().p4().eta()); 
00501   evtselectVarVal.push_back(IntDblPair(51,Obs51));
00502   evtselectVarMatch.push_back(IntBoolPair(51, matchLepb&&matchLepn&&matchLepl));
00503 
00504   if(debug) std::cout << "== observable 51 " << Obs51 << std::endl;
00505   
00506   //obs52
00507   double Obs52  = fabs(solution.getCalLepb().p4().theta()- solution.getRecLepW().p4().theta());
00508   evtselectVarVal.push_back(IntDblPair(52,Obs52));
00509   evtselectVarMatch.push_back(IntBoolPair(52, matchLepb&&matchLepn&&matchLepl));
00510 
00511   if(debug) std::cout << "== observable 52 " << Obs52 << std::endl;
00512   
00513   //obs53
00514   double Obs53  = ROOT::Math::VectorUtil::DeltaR(solution.getCalLepb().p4(), solution.getRecLepW().p4());
00515   evtselectVarVal.push_back(IntDblPair(53,Obs53));
00516   evtselectVarMatch.push_back(IntBoolPair(53, matchLepb&&matchLepn&&matchLepl));
00517 
00518   if(debug) std::cout << "== observable 53 " << Obs53 << std::endl;
00519   
00520   //obs54
00521   double Obs54 = fabs(solution.getCalHadp().p4().phi()- solution.getCalHadq().p4().phi());
00522   if (Obs54 > 3.1415927)  Obs54 =  2*3.1415927 - Obs54;
00523   if (Obs54 < -3.1415927) Obs54 = -2*3.1415927 - Obs54;
00524   evtselectVarVal.push_back(IntDblPair(54,Obs54));
00525   evtselectVarMatch.push_back(IntBoolPair(54, (matchHadp&&matchHadq)||(matchHadpq&&matchHadqp)));
00526 
00527   if(debug) std::cout << "== observable 54 " << Obs54 << std::endl;
00528 
00529   //obs55
00530   double Obs55 = fabs(solution.getCalHadp().p4().eta()- solution.getCalHadq().p4().eta());
00531   evtselectVarVal.push_back(IntDblPair(55,Obs55));
00532   evtselectVarMatch.push_back(IntBoolPair(55, (matchHadp&&matchHadq)||(matchHadpq&&matchHadqp)));
00533 
00534   if(debug) std::cout << "== observable 55 " << Obs55 << std::endl;
00535   
00536   //obs56
00537   double Obs56  = fabs(solution.getCalHadp().p4().theta()- solution.getCalHadq().p4().theta());
00538   evtselectVarVal.push_back(IntDblPair(56,Obs56));
00539   evtselectVarMatch.push_back(IntBoolPair(56, (matchHadp&&matchHadq)||(matchHadpq&&matchHadqp)));
00540 
00541   if(debug) std::cout << "== observable 56 " << Obs56 << std::endl;
00542   
00543   //obs57
00544   double Obs57  = ROOT::Math::VectorUtil::DeltaR(solution.getCalHadp().p4(), solution.getCalHadq().p4());
00545   evtselectVarVal.push_back(IntDblPair(57,Obs57));
00546   evtselectVarMatch.push_back(IntBoolPair(57, (matchHadp&&matchHadq)||(matchHadpq&&matchHadqp)));
00547 
00548   if(debug) std::cout << "== observable 57 " << Obs57 << std::endl;
00549   
00550   //obs58
00551   double Obs58 = -999.;
00552   if (solution.getDecay() == "muon") Obs58  = fabs(solution.getRecLepm().p4().phi()- solution.getRecLepn().p4().phi()); 
00553   if (solution.getDecay() == "electron") Obs58 = fabs(solution.getRecLepe().p4().phi()- solution.getRecLepn().p4().phi()); 
00554   if (Obs58 > 3.1415927)  Obs58 =  2*3.1415927 - Obs58;
00555   if (Obs58 < -3.1415927) Obs58 = -2*3.1415927 - Obs58;
00556   evtselectVarVal.push_back(IntDblPair(58,Obs58));
00557   evtselectVarMatch.push_back(IntBoolPair(58, matchLepl&&matchLepn));
00558 
00559   if(debug) std::cout << "== observable 58 " << Obs58 << std::endl;
00560  
00561   //obs59
00562   double Obs59 = -999.;
00563   if (solution.getDecay() == "muon") Obs59 = fabs(solution.getRecLepm().p4().eta()- solution.getRecLepn().p4().eta()); 
00564   if (solution.getDecay() == "electron") Obs59 = fabs(solution.getRecLepe().p4().eta()- solution.getRecLepn().p4().eta());  
00565   evtselectVarVal.push_back(IntDblPair(59,Obs59));
00566   evtselectVarMatch.push_back(IntBoolPair(59, matchLepl&&matchLepn));
00567 
00568   if(debug) std::cout << "== observable 59 " << Obs59 << std::endl;
00569   
00570   //obs60
00571   double Obs60 = -999.;
00572   if (solution.getDecay() == "muon") Obs60  = fabs(solution.getRecLepm().p4().theta()- solution.getRecLepn().p4().theta());
00573   if (solution.getDecay() == "electron") Obs60  = fabs(solution.getRecLepe().p4().theta()- solution.getRecLepn().p4().theta());
00574   evtselectVarVal.push_back(IntDblPair(60,Obs60));
00575   evtselectVarMatch.push_back(IntBoolPair(60, matchLepl&&matchLepn));
00576 
00577   if(debug) std::cout << "== observable 60 " << Obs60 << std::endl;
00578   
00579   //obs61 
00580   double Obs61 = -999.; 
00581   if (solution.getDecay() == "muon") Obs61 = ROOT::Math::VectorUtil::DeltaR(solution.getRecLepm().p4(), solution.getRecLepn().p4());
00582   if (solution.getDecay() == "electron") Obs61 = ROOT::Math::VectorUtil::DeltaR(solution.getRecLepe().p4(), solution.getRecLepn().p4());
00583   evtselectVarVal.push_back(IntDblPair(61,Obs61));
00584   evtselectVarMatch.push_back(IntBoolPair(61, matchLepl&&matchLepn));
00585 
00586   if(debug) std::cout << "== observable 61 " << Obs61 << std::endl;
00587   
00588   // miscellaneous event variables
00589    
00590 
00591   //obs62
00592   double Obs62  = ((jets->size() > 4 && (*jets)[3].p4().Et() > 0.00001) ? (*jets)[4].p4().Et() / (*jets)[3].p4().Et() : 1.0);
00593   //double Obs62 = 1;
00594   evtselectVarVal.push_back(IntDblPair(62,Obs62));
00595   evtselectVarMatch.push_back(IntBoolPair(62, ((matchHadp&&matchHadq)||(matchHadpq&&matchHadqp))&&matchHadb&&matchLepb));
00596 
00597   if(debug) std::cout << "== observable 62 " << Obs62 << std::endl;
00598   
00599   float calJetsSumEt = 0;
00600   for (unsigned int i = 4; i < jets->size(); i++) {
00601     calJetsSumEt += (*jets)[i].p4().Et();
00602   }
00603   
00604   //obs63
00605   double Obs63_den = (jets->size() > 4) ? ((*jets)[0].p4().Et()+(*jets)[1].p4().Et()+(*jets)[2].p4().Et()+(*jets)[3].p4().Et()+(*jets)[4].p4().Et()) : 0.0;
00606   double Obs63  = (Obs63_den > 0.00001) ? calJetsSumEt / Obs63_den : 1.0;
00607   //double Obs63 =1;  
00608   evtselectVarVal.push_back(IntDblPair(63,Obs63));
00609   evtselectVarMatch.push_back(IntBoolPair(63, ((matchHadp&&matchHadq)||(matchHadpq&&matchHadqp))&&matchHadb&&matchLepb));
00610 
00611   if(debug) std::cout << "== observable 63 " << Obs63 << std::endl;
00612   
00613   //obs64
00614   double Obs64  = solution.getProbChi2();
00615   evtselectVarVal.push_back(IntDblPair(64,Obs64));
00616   evtselectVarMatch.push_back(IntBoolPair(64, ((matchHadp&&matchHadq)||(matchHadpq&&matchHadqp))&&matchHadb&&matchLepb));
00617 
00618   if(debug) std::cout << "== observable 64 " << Obs64 << std::endl;
00619   
00620   //obs65
00621   double Obs65  = solution.getFitHadt().p4().mass() - solution.getCalHadt().p4().mass();
00622   evtselectVarVal.push_back(IntDblPair(65,Obs65));
00623   evtselectVarMatch.push_back(IntBoolPair(65, ((matchHadp&&matchHadq)||(matchHadpq&&matchHadqp))&&matchHadb&&matchLepb)); 
00624 
00625   if(debug) std::cout << "== observable 65 " << Obs65 << std::endl;
00626   
00627   //obs66
00628   double Obs66  = solution.getFitLept().p4().mass() - solution.getCalLept().p4().mass();
00629   evtselectVarVal.push_back(IntDblPair(66,Obs66));
00630   evtselectVarMatch.push_back(IntBoolPair(66, ((matchHadp&&matchHadq)||(matchHadpq&&matchHadqp))&&matchHadb&&matchLepb)); 
00631   
00632   if(debug) std::cout << "observables calculated" << std::endl;
00633 
00634   if (!matchOnly) solution.setLRJetCombObservables(evtselectVarVal);
00635   return evtselectVarMatch;
00636 }