CMSSW Reference Manual

00001 #include "RecoParticleFlow/Benchmark/interface/PFJetBenchmark.h"
00002 
00003 
00004 // preprocessor macro for booking 1d histos with DQMStore -or- bare Root
00005 #define BOOK1D(name,title,nbinsx,lowx,highx) \
00006   h##name = dbe_ ? dbe_->book1D(#name,title,nbinsx,lowx,highx)->getTH1F() \
00007     : new TH1F(#name,title,nbinsx,lowx,highx)
00008 
00009 // preprocessor macro for booking 2d histos with DQMStore -or- bare Root
00010 #define BOOK2D(name,title,nbinsx,lowx,highx,nbinsy,lowy,highy) \
00011   h##name = dbe_ ? dbe_->book2D(#name,title,nbinsx,lowx,highx,nbinsy,lowy,highy)->getTH2F() \
00012     : new TH2F(#name,title,nbinsx,lowx,highx,nbinsy,lowy,highy)
00013 
00014 //macros for building barrel and endcap histos with one call
00015 #define DBOOK1D(name,title,nbinsx,lowx,highx) \
00016   BOOK1D(B##name,"Barrel "#title,nbinsx,lowx,highx); BOOK1D(E##name,"Endcap "#title,nbinsx,lowx,highx);
00017 #define DBOOK2D(name,title,nbinsx,lowx,highx,nbinsy,lowy,highy) \
00018   BOOK2D(B##name,"Barrel "#title,nbinsx,lowx,highx,nbinsy,lowy,highy); BOOK2D(E##name,"Endcap "#title,nbinsx,lowx,highx,nbinsy,lowy,highy);
00019 
00020 // all versions OK
00021 // preprocesor macro for setting axis titles
00022 #define SETAXES(name,xtitle,ytitle) \
00023   h##name->GetXaxis()->SetTitle(xtitle); h##name->GetYaxis()->SetTitle(ytitle)
00024 
00025 //macro for setting the titles for barrel and endcap together
00026 #define DSETAXES(name,xtitle,ytitle) \
00027   SETAXES(B##name,xtitle,ytitle);SETAXES(E##name,xtitle,ytitle)
00028 /*#define SET2AXES(name,xtitle,ytitle) \
00029   hE##name->GetXaxis()->SetTitle(xtitle); hE##name->GetYaxis()->SetTitle(ytitle);  hB##name->GetXaxis()->SetTitle(xtitle); hB##name->GetYaxis()->SetTitle(ytitle)
00030 */
00031 
00032 #define PT (plotAgainstReco_)?"reconstructed P_{T}" :"generated P_{T}"
00033 
00034 using namespace reco;
00035 using namespace std;
00036 
00037 class MonitorElement;
00038 
00039 PFJetBenchmark::PFJetBenchmark() : file_(0) {}
00040 
00041 PFJetBenchmark::~PFJetBenchmark() {
00042   if(file_) file_->Close();
00043 }
00044 
00045 void PFJetBenchmark::write() {
00046    // Store the DAQ Histograms 
00047   if (outputFile_.size() != 0) {
00048     if (dbe_)
00049           dbe_->save(outputFile_.c_str());
00050     // use bare Root if no DQM (FWLite applications)
00051     else if (file_) {
00052        file_->Write(outputFile_.c_str());
00053        cout << "Benchmark output written to file " << outputFile_.c_str() << endl;
00054        file_->Close();
00055        }
00056   }
00057   else 
00058     cout << "No output file specified ("<<outputFile_<<"). Results will not be saved!" << endl;
00059     
00060 } 
00061 
00062 void PFJetBenchmark::setup(
00063                            string Filename,
00064                            bool debug, 
00065                            bool plotAgainstReco,
00066                            bool onlyTwoJets,
00067                            double deltaRMax, 
00068                            string benchmarkLabel_, 
00069                            double recPt, 
00070                            double maxEta, 
00071                            DQMStore * dbe_store) {
00072   debug_ = debug; 
00073   plotAgainstReco_ = plotAgainstReco;
00074   onlyTwoJets_ = onlyTwoJets;
00075   deltaRMax_ = deltaRMax;
00076   outputFile_=Filename;
00077   recPt_cut = recPt;
00078   maxEta_cut= maxEta;
00079   file_ = NULL;
00080   dbe_ = dbe_store;
00081   // print parameters
00082   cout<< "PFJetBenchmark Setup parameters =============================================="<<endl;
00083   cout << "Filename to write histograms " << Filename<<endl;
00084   cout << "PFJetBenchmark debug " << debug_<< endl;
00085   cout << "plotAgainstReco " << plotAgainstReco_ << endl;
00086   cout << "onlyTwoJets " << onlyTwoJets_ << endl;
00087   cout << "deltaRMax " << deltaRMax << endl;
00088   cout << "benchmarkLabel " << benchmarkLabel_ << endl;
00089   cout << "recPt_cut " << recPt_cut << endl;
00090   cout << "maxEta_cut " << maxEta_cut << endl;
00091   
00092   // Book histogram
00093 
00094   // Establish DQM Store
00095   string path = "PFTask/Benchmarks/"+ benchmarkLabel_ + "/";
00096   if (plotAgainstReco) path += "Reco"; else path += "Gen";
00097   if (dbe_) {
00098     dbe_->setCurrentFolder(path.c_str());
00099   }
00100   else {
00101     file_ = new TFile(outputFile_.c_str(), "recreate");
00102 //    TTree * tr = new TTree("PFTast");
00103 //    tr->Branch("Benchmarks/ParticleFlow")
00104     cout << "Info: DQM is not available to provide data storage service. Using TFile to save histograms. "<<endl;
00105   }
00106   // Jets inclusive  distributions  (Pt > 20 or specified recPt GeV)
00107   char cutString[35];
00108   sprintf(cutString,"Jet multiplicity P_{T}>%4.1f GeV", recPt_cut);
00109   BOOK1D(Njets,cutString,50, 0, 50);
00110 
00111   BOOK1D(jetsPt,"Jets P_{T} Distribution",100, 0, 500);
00112 
00113   sprintf(cutString,"Jets #eta Distribution |#eta|<%4.1f", maxEta_cut);
00114   BOOK1D(jetsEta,cutString,100, -5, 5);
00115         
00116   // delta Pt or E quantities for Barrel
00117   DBOOK1D(RPt,#DeltaP_{T}/P_{T},80,-1,1);
00118   DBOOK1D(RCHE,#DeltaE/E (charged had),80,-2,2);
00119   DBOOK1D(RNHE,#DeltaE/E (neutral had),80,-2,2);
00120   DBOOK1D(RNEE,#DeltaE/E (neutral em),80,-2,2);
00121   DBOOK1D(Rneut,#DeltaE/E (neutral),80,-2,2);
00122   DBOOK2D(RPtvsPt,#DeltaP_{T}/P_{T} vs P_{T},250, 0, 500, 100,-2,2);       //used to be 50 bin for each in x-direction
00123   DBOOK2D(RCHEvsPt,#DeltaE/E (charged had) vs P_{T},40, 0, 500, 80,-2,2);
00124   DBOOK2D(RNHEvsPt,#DeltaE/E (neutral had) vs P_{T},40, 0, 500, 80,-2,2);
00125   DBOOK2D(RNEEvsPt,#DeltaE/E (neutral em) vs P_{T},40, 0, 500, 80,-2,2);
00126   DBOOK2D(RneutvsPt,#DeltaE/E (neutral) vs P_{T},40, 0, 500, 80,-2,2);
00127   DBOOK1D(RPt20_40,#DeltaP_{T}/P_{T},80,-1,1);
00128   DBOOK1D(RPt40_60,#DeltaP_{T}/P_{T},80,-1,1);
00129   DBOOK1D(RPt60_80,#DeltaP_{T}/P_{T},80,-1,1);
00130   DBOOK1D(RPt80_100,#DeltaP_{T}/P_{T},80,-1,1);
00131   DBOOK1D(RPt100_150,#DeltaP_{T}/P_{T},80,-1,1);
00132   DBOOK1D(RPt150_200,#DeltaP_{T}/P_{T},80,-1,1);
00133   DBOOK1D(RPt200_250,#DeltaP_{T}/P_{T},80,-1,1);
00134   DBOOK1D(RPt250_300,#DeltaP_{T}/P_{T},80,-1,1);
00135   DBOOK1D(RPt300_400,#DeltaP_{T}/P_{T},80,-1,1);
00136   DBOOK1D(RPt400_500,#DeltaP_{T}/P_{T},80,-1,1);
00137   DBOOK1D(RPt500_750,#DeltaP_{T}/P_{T},80,-1,1);
00138         
00139  // Set Axis Titles
00140  
00141  // Jets inclusive  distributions  (Pt > 20 GeV)
00142   SETAXES(Njets,"","Multiplicity");
00143   SETAXES(jetsPt, PT, "Number of Events");
00144   SETAXES(jetsEta, "#eta", "Number of Events");
00145   // delta Pt or E quantities for Barrel and Endcap
00146   DSETAXES(RPt, "#DeltaP_{T}/P_{T}", "Events");
00147   DSETAXES(RPt20_40, "#DeltaP_{T}/P_{T}", "Events");
00148   DSETAXES(RPt40_60, "#DeltaP_{T}/P_{T}", "Events");
00149   DSETAXES(RPt60_80, "#DeltaP_{T}/P_{T}", "Events");
00150   DSETAXES(RPt80_100, "#DeltaP_{T}/P_{T}", "Events");
00151   DSETAXES(RPt100_150, "#DeltaP_{T}/P_{T}", "Events");
00152   DSETAXES(RPt150_200, "#DeltaP_{T}/P_{T}", "Events");
00153   DSETAXES(RPt200_250, "#DeltaP_{T}/P_{T}", "Events");
00154   DSETAXES(RPt250_300, "#DeltaP_{T}/P_{T}", "Events");
00155   DSETAXES(RPt300_400, "#DeltaP_{T}/P_{T}", "Events");
00156   DSETAXES(RPt400_500, "#DeltaP_{T}/P_{T}", "Events");
00157   DSETAXES(RPt500_750, "#DeltaP_{T}/P_{T}", "Events");
00158   DSETAXES(RCHE, "#DeltaE/E(charged had)", "Events");
00159   DSETAXES(RNHE, "#DeltaE/E(neutral had)", "Events");
00160   DSETAXES(RNEE, "#DeltaE/E(neutral em)", "Events");
00161   DSETAXES(Rneut, "#DeltaE/E(neutral)", "Events");
00162   DSETAXES(RPtvsPt, PT, "#DeltaP_{T}/P_{T}");
00163   DSETAXES(RCHEvsPt, PT, "#DeltaE/E(charged had)");
00164   DSETAXES(RNHEvsPt, PT, "#DeltaE/E(neutral had)");
00165   DSETAXES(RNEEvsPt, PT, "#DeltaE/E(neutral em)");
00166   DSETAXES(RneutvsPt, PT, "#DeltaE/E(neutral)");
00167 
00168 }
00169 
00170 
00171 void PFJetBenchmark::process(const reco::PFJetCollection& pfJets, const reco::GenJetCollection& genJets) {
00172   // loop over reco  pf  jets
00173   resPtMax_ = 0.;
00174   resChargedHadEnergyMax_ = 0.;
00175   resNeutralHadEnergyMax_ = 0.;
00176   resNeutralEmEnergyMax_ = 0.; 
00177   int NPFJets = 0;
00178         
00179   for(unsigned i=0; i<pfJets.size(); i++) {   
00180 
00181     // Count the number of jets with a larger energy
00182     unsigned highJets = 0;
00183     for(unsigned j=0; j<pfJets.size(); j++) { 
00184       if ( j != i && pfJets[j].pt() > pfJets[i].pt() ) highJets++;
00185     }
00186     if ( onlyTwoJets_ && highJets > 1 ) continue;
00187                 
00188                 
00189     const reco::PFJet& pfj = pfJets[i];
00190     double rec_pt = pfj.pt();
00191     double rec_eta = pfj.eta();
00192     double rec_phi = pfj.phi();
00193     // skip PFjets with pt < recPt_cut GeV
00194     if (rec_pt<recPt_cut and recPt_cut != -1.) continue;
00195     // skip PFjets with eta > maxEta_cut
00196     if (fabs(rec_eta)>maxEta_cut and maxEta_cut != -1.) continue;
00197 
00198     NPFJets++;
00199                 
00200     // fill inclusive PFjet distribution pt > 20 GeV
00201     hNjets->Fill(NPFJets);
00202     hjetsPt->Fill(rec_pt);
00203     hjetsEta->Fill(rec_eta);
00204 
00205     // separate Barrel PFJets from Endcap PFJets
00206     bool Barrel = false;
00207     bool Endcap = false;
00208     if (abs(rec_eta) < 1.4 ) Barrel = true;
00209     if (abs (rec_eta) > 1.6 && abs (rec_eta) < 3. ) Endcap = true;
00210 
00211     // do only barrel for now
00212     //  if(!Barrel) continue;
00213 
00214     // look for the closets gen Jet : truth
00215     const GenJet *truth = algo_->matchByDeltaR(&pfj,&genJets);
00216     if(!truth) continue;   
00217     double deltaR = algo_->deltaR(&pfj, truth);
00218     // check deltaR is small enough
00219     if(deltaR < deltaRMax_ || deltaRMax_ == -1.0 ) {//start case deltaR < deltaRMax
00220       // generate histograms comparing the reco and truth candidate (truth = closest in delta-R) 
00221       // get the quantities to place on the denominator and/or divide by
00222       double pt_denom;
00223       double true_pt = truth->pt();
00224       if (plotAgainstReco_) {pt_denom = rec_pt;}
00225       else {pt_denom = true_pt;}
00226       // get true specific quantities
00227       double true_ChargedHadEnergy;
00228       double true_NeutralHadEnergy;
00229       double true_NeutralEmEnergy;
00230       gettrue (truth, true_ChargedHadEnergy, true_NeutralHadEnergy, true_NeutralEmEnergy);
00231       double true_NeutralEnergy = true_NeutralHadEnergy + true_NeutralEmEnergy;
00232       double rec_ChargedHadEnergy = pfj.chargedHadronEnergy();
00233       double rec_NeutralHadEnergy = pfj.neutralHadronEnergy();
00234       double rec_NeutralEmEnergy = pfj.neutralEmEnergy();
00235       double rec_NeutralEnergy = rec_NeutralHadEnergy + rec_NeutralEmEnergy;
00236       bool plot1 = false;
00237       bool plot2 = false;
00238       bool plot3 = false;
00239       bool plot4 = false;
00240       bool plot5 = false;
00241       double cut1 = 0.0001;
00242       double cut2 = 0.0001;
00243       double cut3 = 0.0001;
00244       double cut4 = 0.0001;
00245       double cut5 = 0.0001;
00246       double resPt =0.;
00247       double resChargedHadEnergy= 0.;
00248       double resNeutralHadEnergy= 0.;
00249       double resNeutralEmEnergy= 0.;
00250       double resNeutralEnergy= 0.;
00251       // get relative delta quantities (protect against division by zero!)
00252       if (true_pt > 0.0001){
00253         resPt = (rec_pt -true_pt)/true_pt ; 
00254         plot1 = true;}
00255       if (true_ChargedHadEnergy > cut1){
00256         resChargedHadEnergy = (rec_ChargedHadEnergy- true_ChargedHadEnergy)/true_ChargedHadEnergy;
00257         plot2 = true;}
00258       if (true_NeutralHadEnergy > cut2){
00259         resNeutralHadEnergy = (rec_NeutralHadEnergy- true_NeutralHadEnergy)/true_NeutralHadEnergy;
00260         plot3 = true;}
00261       else 
00262         if (rec_NeutralHadEnergy > cut3){
00263           resNeutralHadEnergy = (rec_NeutralHadEnergy- true_NeutralHadEnergy)/rec_NeutralHadEnergy;
00264           plot3 = true;}
00265       if (true_NeutralEmEnergy > cut4){
00266         resNeutralEmEnergy = (rec_NeutralEmEnergy- true_NeutralEmEnergy)/true_NeutralEmEnergy;
00267         plot4 = true;}
00268       if (true_NeutralEnergy > cut5){
00269         resNeutralEnergy = (rec_NeutralEnergy- true_NeutralEnergy)/true_NeutralEnergy;
00270         plot5 = true;}
00271       
00272       //double deltaEta = algo_->deltaEta(&pfj, truth);
00273       //double deltaPhi = algo_->deltaPhi(&pfj, truth);
00274       if(abs(resPt) > abs(resPtMax_)) resPtMax_ = resPt;
00275       if(abs(resChargedHadEnergy) > abs(resChargedHadEnergyMax_) ) resChargedHadEnergyMax_ = resChargedHadEnergy;
00276       if(abs(resNeutralHadEnergy) > abs(resNeutralHadEnergyMax_) ) resNeutralHadEnergyMax_ = resNeutralHadEnergy;
00277       if(abs(resNeutralEmEnergy) > abs(resNeutralEmEnergyMax_) ) resNeutralEmEnergyMax_ = resNeutralEmEnergy;
00278       if (debug_) {
00279         cout << i <<"  =========PFJet Pt "<< rec_pt
00280              << " eta " << rec_eta
00281              << " phi " << rec_phi
00282              << " Charged Had Energy " << rec_ChargedHadEnergy
00283              << " Neutral Had Energy " << rec_NeutralHadEnergy
00284              << " Neutral elm Energy " << rec_NeutralEmEnergy << endl;
00285         cout << " matching Gen Jet Pt " << true_pt
00286              << " eta " << truth->eta()
00287              << " phi " << truth->phi()
00288              << " Charged Had Energy " << true_ChargedHadEnergy
00289              << " Neutral Had Energy " << true_NeutralHadEnergy
00290              << " Neutral elm Energy " << true_NeutralEmEnergy << endl;
00291         printPFJet(&pfj);
00292         //      cout<<pfj.print()<<endl;
00293         printGenJet(truth);
00294         //cout <<truth->print()<<endl;
00295                                 
00296         cout << "==============deltaR " << deltaR << "  resPt " << resPt
00297              << " resChargedHadEnergy " << resChargedHadEnergy
00298              << " resNeutralHadEnergy " << resNeutralHadEnergy
00299              << " resNeutralEmEnergy " << resNeutralEmEnergy
00300              << endl;
00301       }
00302                         
00303       // fill histograms for relative delta quantitites of matched jets
00304       // delta Pt or E quantities for Barrel
00305       if (Barrel){
00306         if(plot1) { 
00307           hBRPt->Fill (resPt);
00308           if ( pt_denom >  20. && pt_denom <  40. ) hBRPt20_40->Fill (resPt);
00309           if ( pt_denom >  40. && pt_denom <  60. ) hBRPt40_60->Fill (resPt);
00310           if ( pt_denom >  60. && pt_denom <  80. ) hBRPt60_80->Fill (resPt);
00311           if ( pt_denom >  80. && pt_denom < 100. ) hBRPt80_100->Fill (resPt);
00312           if ( pt_denom > 100. && pt_denom < 150. ) hBRPt100_150->Fill (resPt);
00313           if ( pt_denom > 150. && pt_denom < 200. ) hBRPt150_200->Fill (resPt);
00314           if ( pt_denom > 200. && pt_denom < 250. ) hBRPt200_250->Fill (resPt);
00315           if ( pt_denom > 250. && pt_denom < 300. ) hBRPt250_300->Fill (resPt);
00316           if ( pt_denom > 300. && pt_denom < 400. ) hBRPt300_400->Fill (resPt);
00317           if ( pt_denom > 400. && pt_denom < 500. ) hBRPt400_500->Fill (resPt);
00318           if ( pt_denom > 500. && pt_denom < 750. ) hBRPt500_750->Fill (resPt);
00319         }
00320         if(plot2)hBRCHE->Fill(resChargedHadEnergy);
00321         if(plot3)hBRNHE->Fill(resNeutralHadEnergy);
00322         if(plot4)hBRNEE->Fill(resNeutralEmEnergy);
00323         if(plot5)hBRneut->Fill(resNeutralEnergy);
00324         if(plot1)hBRPtvsPt->Fill(pt_denom, resPt);
00325         if(plot2)hBRCHEvsPt->Fill(pt_denom, resChargedHadEnergy);
00326         if(plot3)hBRNHEvsPt->Fill(pt_denom, resNeutralHadEnergy);
00327         if(plot4)hBRNEEvsPt->Fill(pt_denom, resNeutralEmEnergy);
00328         if(plot5)hBRneutvsPt->Fill(pt_denom, resNeutralEnergy);
00329       }
00330       // delta Pt or E quantities for Endcap
00331       if (Endcap){
00332         if(plot1) {
00333           hERPt->Fill (resPt);
00334           if ( pt_denom >  20. && pt_denom <  40. ) hERPt20_40->Fill (resPt);
00335           if ( pt_denom >  40. && pt_denom <  60. ) hERPt40_60->Fill (resPt);
00336           if ( pt_denom >  60. && pt_denom <  80. ) hERPt60_80->Fill (resPt);
00337           if ( pt_denom >  80. && pt_denom < 100. ) hERPt80_100->Fill (resPt);
00338           if ( pt_denom > 100. && pt_denom < 150. ) hERPt100_150->Fill (resPt);
00339           if ( pt_denom > 150. && pt_denom < 200. ) hERPt150_200->Fill (resPt);
00340           if ( pt_denom > 200. && pt_denom < 250. ) hERPt200_250->Fill (resPt);
00341           if ( pt_denom > 250. && pt_denom < 300. ) hERPt250_300->Fill (resPt);
00342           if ( pt_denom > 300. && pt_denom < 400. ) hERPt300_400->Fill (resPt);
00343           if ( pt_denom > 400. && pt_denom < 500. ) hERPt400_500->Fill (resPt);
00344           if ( pt_denom > 500. && pt_denom < 750. ) hERPt500_750->Fill (resPt);
00345         }
00346         if(plot2)hERCHE->Fill(resChargedHadEnergy);
00347         if(plot3)hERNHE->Fill(resNeutralHadEnergy);
00348         if(plot4)hERNEE->Fill(resNeutralEmEnergy);
00349         if(plot5)hERneut->Fill(resNeutralEnergy);
00350         if(plot1)hERPtvsPt->Fill(pt_denom, resPt);
00351         if(plot2)hERCHEvsPt->Fill(pt_denom, resChargedHadEnergy);
00352         if(plot3)hERNHEvsPt->Fill(pt_denom, resNeutralHadEnergy);
00353         if(plot4)hERNEEvsPt->Fill(pt_denom, resNeutralEmEnergy);
00354         if(plot5)hERneutvsPt->Fill(pt_denom, resNeutralEnergy);
00355       }                                         
00356     } // end case deltaR < deltaRMax
00357                 
00358   } // i loop on pf Jets        
00359 }
00360 
00361 void PFJetBenchmark::gettrue (const reco::GenJet* truth, double& true_ChargedHadEnergy, 
00362                               double& true_NeutralHadEnergy, double& true_NeutralEmEnergy){
00363   std::vector <const GenParticle*> mcparts = truth->getGenConstituents ();
00364   true_NeutralEmEnergy = 0.;
00365   true_ChargedHadEnergy = 0.;
00366   true_NeutralHadEnergy = 0.;
00367   // for each MC particle in turn  
00368   for (unsigned i = 0; i < mcparts.size (); i++) {
00369     const GenParticle* mcpart = mcparts[i];
00370     int PDG = abs( mcpart->pdgId());
00371     double e = mcpart->energy(); 
00372     switch(PDG){  // start PDG switch
00373     case 22: // photon
00374       true_NeutralEmEnergy += e;
00375       break;
00376     case 211: // pi
00377     case 321: // K
00378     case 2212: // p
00379     case 11: //electrons (until recognised)
00380       true_ChargedHadEnergy += e;
00381       break;
00382     case 310: // K_S0
00383     case 130: // K_L0
00384     case 3122: // Lambda0
00385     case 2112: // n0
00386       true_NeutralHadEnergy += e;
00387     default:
00388       break;
00389     }  // end PDG switch                
00390   }  // end loop on constituents.
00391 }
00392 
00393 void PFJetBenchmark::printPFJet(const reco::PFJet* pfj){
00394   cout<<setiosflags(ios::right);
00395   cout<<setiosflags(ios::fixed);
00396   cout<<setprecision(3);
00397 
00398 //formerly read out all the constituents, but now there is an error ??????
00399 //std::vector <const reco::PFCandidate*> pfCandidates = pfj->getPFConstituents ();
00400 
00401   cout << "PFJet  p/px/py/pz/pt: " << pfj->p() << "/" << pfj->px () 
00402        << "/" << pfj->py() << "/" << pfj->pz() << "/" << pfj->pt() << endl
00403        << "    eta/phi: " << pfj->eta () << "/" << pfj->phi () << endl                  
00404        << "    PFJet specific:" << std::endl
00405        << "      charged/neutral hadrons energy: " << pfj->chargedHadronEnergy () << '/' << pfj->neutralHadronEnergy () << endl
00406        << "      charged/neutral em energy: " << pfj->chargedEmEnergy () << '/' << pfj->neutralEmEnergy () << endl
00407        << "      charged muon energy: " << pfj->chargedMuEnergy () << '/' << endl
00408        << "      charged/neutral multiplicity: " << pfj->chargedMultiplicity () << '/' << pfj->neutralMultiplicity () << endl;
00409 /* cout  << "    # of pfCandidates: " << pfCandidates.size() << endl;
00410 
00411 //  vector <PFBlockRef> PFBRef;
00412 // print PFCandidates constituents of the jet
00413   for(unsigned i=0; i<pfCandidates.size(); i++) {
00414     const PFCandidate* pfCand = pfCandidates[i];
00415     cout<<i <<" " << *pfCand << endl;
00416   } // end loop on i (PFCandidates)
00417 */
00418   // print blocks associated to the jet (to be done with new format PFCandiates)
00419         
00420   cout<<resetiosflags(ios::right|ios::fixed);
00421 }
00422 
00423 
00424 void PFJetBenchmark::printGenJet (const reco::GenJet* truth){
00425   std::vector <const GenParticle*> mcparts = truth->getGenConstituents ();
00426   cout << "GenJet p/px/py/pz/pt: " << truth->p() << '/' << truth->px () 
00427        << '/' << truth->py() << '/' << truth->pz() << '/' << truth->pt() << endl
00428        << "    eta/phi: " << truth->eta () << '/' << truth->phi () << endl
00429        << "    # of constituents: " << mcparts.size() << endl;
00430   cout << "    constituents:" << endl;
00431   for (unsigned i = 0; i < mcparts.size (); i++) {
00432     const GenParticle* mcpart = mcparts[i];
00433     cout << "      #" << i << "  PDG code:" << mcpart->pdgId() 
00434          << ", p/pt/eta/phi: " << mcpart->p() << '/' << mcpart->pt() 
00435          << '/' << mcpart->eta() << '/' << mcpart->phi() << endl;       
00436   }    
00437 }
00438
PFJetBenchmark.cc
