CMSSW Reference Manual

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