Reference Manual

Go to the documentation of this file.
00001 #include "RecoParticleFlow/Benchmark/interface/PFJetBenchmark.h"
00002 #include "DataFormats/TrackReco/interface/Track.h"
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); BOOK1D(F##name,"Forward "#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); BOOK2D(F##name,"Forward "#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);SETAXES(F##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 #define P (plotAgainstReco_)?"generated P" :"generated P"
00034 
00035 using namespace reco;
00036 using namespace std;
00037 
00038 class MonitorElement;
00039 
00040 PFJetBenchmark::PFJetBenchmark() : file_(0), entry_(0) {}
00041 
00042 PFJetBenchmark::~PFJetBenchmark() {
00043   if(file_) file_->Close();
00044 }
00045 
00046 void PFJetBenchmark::write() {
00047    // Store the DAQ Histograms 
00048   if (outputFile_.size() != 0) {
00049     if (dbe_)
00050           dbe_->save(outputFile_.c_str());
00051     // use bare Root if no DQM (FWLite applications)
00052     else if (file_) {
00053        file_->Write(outputFile_.c_str());
00054        cout << "Benchmark output written to file " << outputFile_.c_str() << endl;
00055        file_->Close();
00056        }
00057   }
00058   else 
00059     cout << "No output file specified ("<<outputFile_<<"). Results will not be saved!" << endl;
00060     
00061 } 
00062 
00063 void PFJetBenchmark::setup(
00064                            string Filename,
00065                            bool debug, 
00066                            bool plotAgainstReco,
00067                            bool onlyTwoJets,
00068                            double deltaRMax, 
00069                            string benchmarkLabel_, 
00070                            double recPt, 
00071                            double maxEta, 
00072                            DQMStore * dbe_store) {
00073   debug_ = debug; 
00074   plotAgainstReco_ = plotAgainstReco;
00075   onlyTwoJets_ = onlyTwoJets;
00076   deltaRMax_ = deltaRMax;
00077   outputFile_=Filename;
00078   recPt_cut = recPt;
00079   maxEta_cut= maxEta;
00080   file_ = NULL;
00081   dbe_ = dbe_store;
00082   // print parameters
00083   cout<< "PFJetBenchmark Setup parameters =============================================="<<endl;
00084   cout << "Filename to write histograms " << Filename<<endl;
00085   cout << "PFJetBenchmark debug " << debug_<< endl;
00086   cout << "plotAgainstReco " << plotAgainstReco_ << endl;
00087   cout << "onlyTwoJets " << onlyTwoJets_ << endl;
00088   cout << "deltaRMax " << deltaRMax << endl;
00089   cout << "benchmarkLabel " << benchmarkLabel_ << endl;
00090   cout << "recPt_cut " << recPt_cut << endl;
00091   cout << "maxEta_cut " << maxEta_cut << endl;
00092   
00093   // Book histogram
00094 
00095   // Establish DQM Store
00096   string path = "PFTask/Benchmarks/"+ benchmarkLabel_ + "/";
00097   if (plotAgainstReco) path += "Reco"; else path += "Gen";
00098   if (dbe_) {
00099     dbe_->setCurrentFolder(path.c_str());
00100   }
00101   else {
00102     file_ = new TFile(outputFile_.c_str(), "recreate");
00103 //    TTree * tr = new TTree("PFTast");
00104 //    tr->Branch("Benchmarks/ParticleFlow")
00105     cout << "Info: DQM is not available to provide data storage service. Using TFile to save histograms. "<<endl;
00106   }
00107   // Jets inclusive  distributions  (Pt > 20 or specified recPt GeV)
00108   char cutString[35];
00109   sprintf(cutString,"Jet multiplicity P_{T}>%4.1f GeV", recPt_cut);
00110   BOOK1D(Njets,cutString,50, 0, 50);
00111 
00112   BOOK1D(jetsPt,"Jets P_{T} Distribution",100, 0, 500);
00113 
00114   sprintf(cutString,"Jets #eta Distribution |#eta|<%4.1f", maxEta_cut);
00115   BOOK1D(jetsEta,cutString,100, -5, 5);
00116         
00117   BOOK2D(RPtvsEta,"#DeltaP_{T}/P_{T} vs #eta",200, -5., 5., 200,-1,1); 
00118   BOOK2D(RNeutvsEta,"R_{Neutral} vs #eta",200, -5., 5., 200,-1,1); 
00119   BOOK2D(RNEUTvsEta,"R_{HCAL+ECAL} vs #eta",200, -5., 5., 200,-1,1); 
00120   BOOK2D(RNONLvsEta,"R_{HCAL+ECAL - Hcal Only} vs #eta",200, -5., 5., 200,-1,1); 
00121   BOOK2D(RHCALvsEta,"R_{HCAL} vs #eta",200, -5., 5., 200,-1,1); 
00122   BOOK2D(RHONLvsEta,"R_{HCAL only} vs #eta",200, -5., 5., 200,-1,1); 
00123   BOOK2D(RCHEvsEta,"R_{Charged} vs #eta",200, -5., 5., 200,-1,1); 
00124   BOOK2D(NCHvsEta,"N_{Charged} vs #eta",200, -5., 5., 200,0.,200.);
00125   BOOK2D(NCH0vsEta,"N_{Charged} vs #eta, iter 0",200, -5., 5., 200,0.,200.);
00126   BOOK2D(NCH1vsEta,"N_{Charged} vs #eta, iter 1",200, -5., 5., 200,0.,200.);
00127   BOOK2D(NCH2vsEta,"N_{Charged} vs #eta, iter 2",200, -5., 5., 200,0.,200.);
00128   BOOK2D(NCH3vsEta,"N_{Charged} vs #eta, iter 3",200, -5., 5., 200,0.,200.);
00129   BOOK2D(NCH4vsEta,"N_{Charged} vs #eta, iter 4",200, -5., 5., 200,0.,200.);
00130   BOOK2D(NCH5vsEta,"N_{Charged} vs #eta, iter 5",200, -5., 5., 200,0.,200.);
00131   BOOK2D(NCH6vsEta,"N_{Charged} vs #eta, iter 6",200, -5., 5., 200,0.,200.);
00132   // delta Pt or E quantities for Barrel
00133   DBOOK1D(RPt,#DeltaP_{T}/P_{T},80,-1,1);
00134   DBOOK1D(RCHE,#DeltaE/E (charged had),80,-2,2);
00135   DBOOK1D(RNHE,#DeltaE/E (neutral had),80,-2,2);
00136   DBOOK1D(RNEE,#DeltaE/E (neutral em),80,-2,2);
00137   DBOOK1D(Rneut,#DeltaE/E (neutral),80,-2,2);
00138   DBOOK1D(NCH, #N_{charged},200,0.,200.);
00139   DBOOK2D(RPtvsPt,#DeltaP_{T}/P_{T} vs P_{T},250, 0, 500, 200,-1,1);       //used to be 50 bin for each in x-direction
00140   DBOOK2D(RCHEvsPt,#DeltaE/E (charged had) vs P_{T},250, 0, 500, 120,-1,2);
00141   DBOOK2D(RNHEvsPt,#DeltaE/E (neutral had) vs P_{T},250, 0, 500, 120,-1,2);
00142   DBOOK2D(RNEEvsPt,#DeltaE/E (neutral em) vs P_{T},250, 0, 500, 120,-1,2);
00143   DBOOK2D(RneutvsPt,#DeltaE/E (neutral) vs P_{T},250, 0, 500, 120,-1,2);
00144   DBOOK2D(NCHvsPt,N_{charged} vs P_{T},250,0,500,200,0.,200.);
00145   DBOOK2D(NCH0vsPt, N_{charged} vs P_{T} iter 0,250,0,500,200,0.,200.);
00146   DBOOK2D(NCH1vsPt, N_{charged} vs P_{T} iter 1,250,0,500,200,0.,200.);
00147   DBOOK2D(NCH2vsPt, N_{charged} vs P_{T} iter 2,250,0,500,200,0.,200.);
00148   DBOOK2D(NCH3vsPt, N_{charged} vs P_{T} iter 3,250,0,500,200,0.,200.);
00149   DBOOK2D(NCH4vsPt, N_{charged} vs P_{T} iter 4,250,0,500,200,0.,200.);
00150   DBOOK2D(NCH5vsPt, N_{charged} vs P_{T} iter 5,250,0,500,200,0.,200.);
00151   DBOOK2D(NCH6vsPt, N_{charged} vs P_{T} iter 6,250,0,500,200,0.,200.);
00152   
00153 
00154   DBOOK2D(RNEUTvsP,#DeltaE/E (ECAL+HCAL) vs P,250, 0, 1000, 150,-1.5,1.5);
00155   DBOOK2D(RNONLvsP,#DeltaE/E (ECAL+HCAL-only) vs P,250, 0, 1000, 150,-1.5,1.5);
00156   DBOOK2D(RHCALvsP,#DeltaE/E (HCAL) vs P,250, 0, 1000, 150,-1.5,1.5);
00157   DBOOK2D(RHONLvsP,#DeltaE/E (HCAL only) vs P,250, 0, 1000, 150,-1.5,1.5);
00158   DBOOK1D(RPt20_40,#DeltaP_{T}/P_{T},80,-1,1);
00159   DBOOK1D(RPt40_60,#DeltaP_{T}/P_{T},80,-1,1);
00160   DBOOK1D(RPt60_80,#DeltaP_{T}/P_{T},80,-1,1);
00161   DBOOK1D(RPt80_100,#DeltaP_{T}/P_{T},80,-1,1);
00162   DBOOK1D(RPt100_150,#DeltaP_{T}/P_{T},80,-1,1);
00163   DBOOK1D(RPt150_200,#DeltaP_{T}/P_{T},80,-1,1);
00164   DBOOK1D(RPt200_250,#DeltaP_{T}/P_{T},80,-1,1);
00165   DBOOK1D(RPt250_300,#DeltaP_{T}/P_{T},80,-1,1);
00166   DBOOK1D(RPt300_400,#DeltaP_{T}/P_{T},160,-1,1);
00167   DBOOK1D(RPt400_500,#DeltaP_{T}/P_{T},160,-1,1);
00168   DBOOK1D(RPt500_750,#DeltaP_{T}/P_{T},160,-1,1);
00169   DBOOK1D(RPt750_1250,#DeltaP_{T}/P_{T},160,-1,1);
00170   DBOOK1D(RPt1250_2000,#DeltaP_{T}/P_{T},160,-1,1);
00171   DBOOK1D(RPt2000_5000,#DeltaP_{T}/P_{T},160,-1,1);
00172 
00173   DBOOK2D(DEtavsPt,#Delta#eta vs P_{T},1000,0,2000,500,-0.5,0.5);
00174   DBOOK2D(DPhivsPt,#Delta#phi vs P_{T},1000,0,2000,500,-0.5,0.5);
00175   BOOK2D(DEtavsEta,"#Delta#eta vs P_{T}",1000,-5.,+5.,500,-0.5,0.5);
00176   BOOK2D(DPhivsEta,"#Delta#phi vs P_{T}",1000,-5.,+5.,500,-0.5,0.5);
00177         
00178  // Set Axis Titles
00179  
00180  // Jets inclusive  distributions  (Pt > 20 GeV)
00181   SETAXES(Njets,"","Multiplicity");
00182   SETAXES(jetsPt, PT, "Number of Events");
00183   SETAXES(jetsEta, "#eta", "Number of Events");
00184   SETAXES(RNeutvsEta, "#eta", "#DeltaE/E (Neutral)");
00185   SETAXES(RNEUTvsEta, "#eta", "#DeltaE/E (ECAL+HCAL)");
00186   SETAXES(RNONLvsEta, "#eta", "#DeltaE/E (ECAL+HCAL-only)");
00187   SETAXES(RHCALvsEta, "#eta", "#DeltaE/E (HCAL)");
00188   SETAXES(RHONLvsEta, "#eta", "#DeltaE/E (HCAL Only)");
00189   SETAXES(RCHEvsEta, "#eta", "#DeltaE/E (Charged)");
00190   SETAXES(RPtvsEta, "#eta", "#DeltaP_{T}/P_{T}");
00191   SETAXES(DEtavsEta, "#eta", "#Delta#eta");
00192   SETAXES(DPhivsEta,"#eta", "#Delta#phi");
00193   // delta Pt or E quantities for Barrel and Endcap
00194   DSETAXES(RPt, "#DeltaP_{T}/P_{T}", "Events");
00195   DSETAXES(RPt20_40, "#DeltaP_{T}/P_{T}", "Events");
00196   DSETAXES(RPt40_60, "#DeltaP_{T}/P_{T}", "Events");
00197   DSETAXES(RPt60_80, "#DeltaP_{T}/P_{T}", "Events");
00198   DSETAXES(RPt80_100, "#DeltaP_{T}/P_{T}", "Events");
00199   DSETAXES(RPt100_150, "#DeltaP_{T}/P_{T}", "Events");
00200   DSETAXES(RPt150_200, "#DeltaP_{T}/P_{T}", "Events");
00201   DSETAXES(RPt200_250, "#DeltaP_{T}/P_{T}", "Events");
00202   DSETAXES(RPt250_300, "#DeltaP_{T}/P_{T}", "Events");
00203   DSETAXES(RPt300_400, "#DeltaP_{T}/P_{T}", "Events");
00204   DSETAXES(RPt400_500, "#DeltaP_{T}/P_{T}", "Events");
00205   DSETAXES(RPt500_750, "#DeltaP_{T}/P_{T}", "Events");
00206   DSETAXES(RPt750_1250, "#DeltaP_{T}/P_{T}", "Events");
00207   DSETAXES(RPt1250_2000, "#DeltaP_{T}/P_{T}", "Events");
00208   DSETAXES(RPt2000_5000, "#DeltaP_{T}/P_{T}", "Events");
00209   DSETAXES(RCHE, "#DeltaE/E(charged had)", "Events");
00210   DSETAXES(RNHE, "#DeltaE/E(neutral had)", "Events");
00211   DSETAXES(RNEE, "#DeltaE/E(neutral em)", "Events");
00212   DSETAXES(Rneut, "#DeltaE/E(neutral)", "Events");
00213   DSETAXES(RPtvsPt, PT, "#DeltaP_{T}/P_{T}");
00214   DSETAXES(RCHEvsPt, PT, "#DeltaE/E(charged had)");
00215   DSETAXES(RNHEvsPt, PT, "#DeltaE/E(neutral had)");
00216   DSETAXES(RNEEvsPt, PT, "#DeltaE/E(neutral em)");
00217   DSETAXES(RneutvsPt, PT, "#DeltaE/E(neutral)");
00218   DSETAXES(RHCALvsP, P, "#DeltaE/E(HCAL)");
00219   DSETAXES(RHONLvsP, P, "#DeltaE/E(HCAL-only)");
00220   DSETAXES(RNEUTvsP, P, "#DeltaE/E(ECAL+HCAL)");
00221   DSETAXES(RNONLvsP, P, "#DeltaE/E(ECAL+HCAL-only)");
00222   DSETAXES(DEtavsPt, PT, "#Delta#eta");
00223   DSETAXES(DPhivsPt, PT, "#Delta#phi");
00224 
00225 }
00226 
00227 
00228 void PFJetBenchmark::process(const reco::PFJetCollection& pfJets, const reco::GenJetCollection& genJets) {
00229   // loop over reco  pf  jets
00230   resPtMax_ = 0.;
00231   resChargedHadEnergyMax_ = 0.;
00232   resNeutralHadEnergyMax_ = 0.;
00233   resNeutralEmEnergyMax_ = 0.; 
00234   int NPFJets = 0;
00235         
00236   for(unsigned i=0; i<pfJets.size(); i++) {   
00237 
00238     // Count the number of jets with a larger energy
00239     unsigned highJets = 0;
00240     for(unsigned j=0; j<pfJets.size(); j++) { 
00241       if ( j != i && pfJets[j].pt() > pfJets[i].pt() ) highJets++;
00242     }
00243     if ( onlyTwoJets_ && highJets > 1 ) continue;
00244                 
00245                 
00246     const reco::PFJet& pfj = pfJets[i];
00247     double rec_pt = pfj.pt();
00248     double rec_eta = pfj.eta();
00249     double rec_phi = pfj.phi();
00250 
00251     // skip PFjets with pt < recPt_cut GeV
00252     if (rec_pt<recPt_cut and recPt_cut != -1.) continue;
00253     // skip PFjets with eta > maxEta_cut
00254     if (fabs(rec_eta)>maxEta_cut and maxEta_cut != -1.) continue;
00255 
00256     NPFJets++;
00257                 
00258     // fill inclusive PFjet distribution pt > 20 GeV
00259     hNjets->Fill(NPFJets);
00260     hjetsPt->Fill(rec_pt);
00261     hjetsEta->Fill(rec_eta);
00262 
00263     // separate Barrel PFJets from Endcap PFJets
00264     bool Barrel = false;
00265     bool Endcap = false;
00266     bool Forward = false;
00267     if (std::abs(rec_eta) < 1.4 ) Barrel = true;
00268     if (std::abs (rec_eta) > 1.6 && std::abs (rec_eta) < 2.4 ) Endcap = true;
00269     if (std::abs (rec_eta) > 2.5 && std::abs (rec_eta) < 2.9 ) Forward = true;
00270     if (std::abs (rec_eta) > 3.1 && std::abs (rec_eta) < 4.7 ) Forward = true;
00271 
00272     // do only barrel for now
00273     //  if(!Barrel) continue;
00274 
00275     // look for the closets gen Jet : truth
00276     const GenJet *truth = algo_->matchByDeltaR(&pfj,&genJets);
00277     if(!truth) continue;   
00278     double deltaR = algo_->deltaR(&pfj, truth);
00279     // check deltaR is small enough
00280     if(deltaR < deltaRMax_ || (abs(rec_eta)>2.5 && deltaR < 0.2) || deltaRMax_ == -1.0 ) {//start case deltaR < deltaRMax
00281 
00282       // generate histograms comparing the reco and truth candidate (truth = closest in delta-R) 
00283       // get the quantities to place on the denominator and/or divide by
00284       double pt_denom;
00285       double true_E = truth->p();
00286       double true_pt = truth->pt();
00287       double true_eta = truth->eta();
00288       double true_phi = truth->phi();
00289 
00290       if (plotAgainstReco_) {pt_denom = rec_pt;}
00291       else {pt_denom = true_pt;}
00292       // get true specific quantities
00293       double true_ChargedHadEnergy;
00294       double true_NeutralHadEnergy;
00295       double true_NeutralEmEnergy;
00296       gettrue (truth, true_ChargedHadEnergy, true_NeutralHadEnergy, true_NeutralEmEnergy);
00297       double true_NeutralEnergy = true_NeutralHadEnergy + true_NeutralEmEnergy;
00298       double rec_ChargedHadEnergy = pfj.chargedHadronEnergy();
00299       double rec_NeutralHadEnergy = pfj.neutralHadronEnergy();
00300       double rec_NeutralEmEnergy = pfj.neutralEmEnergy();
00301       double rec_NeutralEnergy = rec_NeutralHadEnergy + rec_NeutralEmEnergy;
00302       double rec_ChargedMultiplicity = pfj.chargedMultiplicity();
00303       std::vector <PFCandidatePtr> constituents = pfj.getPFConstituents ();
00304       std::vector <unsigned int> chMult(7, static_cast<unsigned int>(0)); 
00305       for (unsigned ic = 0; ic < constituents.size (); ++ic) {
00306         if ( constituents[ic]->particleId() > 3 ) continue;
00307         reco::TrackRef trackRef = constituents[ic]->trackRef();
00308         if ( trackRef.isNull() ) {
00309           //std::cout << "Warning in entry " << entry_ 
00310           //        << " : a track with Id " << constituents[ic]->particleId() 
00311           //        << " has no track ref.." << std::endl;
00312           continue;
00313         }
00314         unsigned int iter = 0; 
00315         switch (trackRef->algo()) {
00316         case TrackBase::ctf:
00317         case TrackBase::iter0:
00318           iter = 0;
00319           break;
00320         case TrackBase::iter1:
00321           iter = 1;
00322           break;
00323         case TrackBase::iter2:
00324           iter = 2;
00325           break;
00326         case TrackBase::iter3:
00327           iter = 3;
00328           break;
00329         case TrackBase::iter4:
00330           iter = 4;
00331           break;
00332         case TrackBase::iter5:
00333           iter = 5;
00334           break;
00335         default:
00336           iter = 6;
00337           std::cout << "Warning in entry " << entry_ << " : iter = 6... " << std::endl;
00338           break;
00339         }
00340         ++(chMult[iter]);
00341       }
00342 
00343       bool plot1 = false;
00344       bool plot2 = false;
00345       bool plot3 = false;
00346       bool plot4 = false;
00347       bool plot5 = false;
00348       bool plot6 = false;
00349       bool plot7 = false;
00350       double cut1 = 0.0001;
00351       double cut2 = 0.0001;
00352       double cut3 = 0.0001;
00353       double cut4 = 0.0001;
00354       double cut5 = 0.0001;
00355       double cut6 = 0.0001;
00356       double cut7 = 0.0001;
00357       double resPt =0.;
00358       double resChargedHadEnergy= 0.;
00359       double resNeutralHadEnergy= 0.;
00360       double resNeutralEmEnergy= 0.;
00361       double resNeutralEnergy= 0.;
00362 
00363       double resHCALEnergy = 0.;
00364       double resNEUTEnergy = 0.;
00365       if ( rec_NeutralHadEnergy > cut6 && rec_ChargedHadEnergy < cut1 ) { 
00366         double true_NEUTEnergy = true_NeutralHadEnergy + true_NeutralEmEnergy;
00367         double true_HCALEnergy = true_NEUTEnergy - rec_NeutralEmEnergy;
00368         double rec_NEUTEnergy = rec_NeutralHadEnergy+rec_NeutralEmEnergy; 
00369         double rec_HCALEnergy = rec_NeutralHadEnergy; 
00370         resHCALEnergy = (rec_HCALEnergy-true_HCALEnergy)/rec_HCALEnergy;
00371         resNEUTEnergy = (rec_NEUTEnergy-true_NEUTEnergy)/rec_NEUTEnergy;
00372         if ( rec_NeutralEmEnergy > cut7 ) {
00373           plot6 = true;
00374         } else {
00375           plot7 = true;
00376         }
00377       }
00378 
00379       // get relative delta quantities (protect against division by zero!)
00380       if (true_pt > 0.0001){
00381         resPt = (rec_pt -true_pt)/true_pt ; 
00382         plot1 = true;}
00383       if (true_ChargedHadEnergy > cut1){
00384         resChargedHadEnergy = (rec_ChargedHadEnergy- true_ChargedHadEnergy)/true_ChargedHadEnergy;
00385         plot2 = true;}
00386       if (true_NeutralHadEnergy > cut2){
00387         resNeutralHadEnergy = (rec_NeutralHadEnergy- true_NeutralHadEnergy)/true_NeutralHadEnergy;
00388         plot3 = true;}
00389       else 
00390         if (rec_NeutralHadEnergy > cut3){
00391           resNeutralHadEnergy = (rec_NeutralHadEnergy- true_NeutralHadEnergy)/rec_NeutralHadEnergy;
00392           plot3 = true;}
00393       if (true_NeutralEmEnergy > cut4){
00394         resNeutralEmEnergy = (rec_NeutralEmEnergy- true_NeutralEmEnergy)/true_NeutralEmEnergy;
00395         plot4 = true;}
00396       if (true_NeutralEnergy > cut5){
00397         resNeutralEnergy = (rec_NeutralEnergy- true_NeutralEnergy)/true_NeutralEnergy;
00398         plot5 = true;}
00399       
00400       //double deltaEta = algo_->deltaEta(&pfj, truth);
00401       //double deltaPhi = algo_->deltaPhi(&pfj, truth);
00402 
00403       // Print outliers for further debugging
00404       if ( ( resPt > 0.2 && true_pt > 100. ) || 
00405            ( resPt < -0.5 && true_pt > 100. ) ) {
00406         //if ( ( true_pt > 50. && 
00407         //     ( ( truth->eta()>3.0 && rec_eta-truth->eta() < -0.1 ) || 
00408         //       ( truth->eta()<-3.0 && rec_eta-truth->eta() > 0.1 ) ))) {
00409         std::cout << "Entry " << entry_ 
00410                   << " resPt = " << resPt
00411                   <<" resCharged  " << resChargedHadEnergy
00412                   <<" resNeutralHad  " << resNeutralHadEnergy
00413                   << " resNeutralEm  " << resNeutralEmEnergy
00414                   << " pT (T/R) " << true_pt << "/" << rec_pt 
00415                   << " Eta (T/R) " << truth->eta() << "/" << rec_eta 
00416                   << " Phi (T/R) " << truth->phi() << "/" << rec_phi 
00417                   << std::endl;
00418 
00419         // check overlapping PF jets
00420         const reco::PFJet* pfoj = 0; 
00421         double dRo = 1E9;
00422         for(unsigned j=0; j<pfJets.size(); j++) { 
00423           const reco::PFJet& pfo = pfJets[j];
00424           if ( j != i &&  algo_->deltaR(&pfj,&pfo) < dRo && pfo.pt() > 0.25*pfj.pt()) { 
00425             dRo = algo_->deltaR(&pfj,&pfo);     
00426             pfoj = &pfo;
00427           }
00428         }
00429         
00430         // Check overlapping Gen Jet 
00431         math::XYZTLorentzVector overlappinGenJet(0.,0.,0.,0.);
00432         const reco::GenJet* genoj = 0;
00433         double dRgo = 1E9;
00434         for(unsigned j=0; j<genJets.size(); j++) { 
00435           const reco::GenJet* gjo = &(genJets[j]);
00436           if ( gjo != truth && algo_->deltaR(truth,gjo) < dRgo && gjo->pt() > 0.25*truth->pt() ) { 
00437             dRgo = algo_->deltaR(truth,gjo);
00438             genoj = gjo;
00439           }
00440         }
00441         
00442         if ( dRo < 0.8 && dRgo < 0.8 && algo_->deltaR(genoj,pfoj) < 2.*deltaRMax_ ) 
00443           std::cout << "Excess probably due to overlapping jets (DR = " <<   algo_->deltaR(genoj,pfoj) << "),"
00444                     << " at DeltaR(T/R) = " << dRgo << "/" << dRo  
00445                     << " with pT(T/R) " << genoj->pt() << "/" << pfoj->pt()
00446                     << " and Eta (T/R) " << genoj->eta() << "/" << pfoj->eta()
00447                     << " and Phi (T/R) " << genoj->phi() << "/" << pfoj->phi()
00448                     << std::endl;
00449       }
00450 
00451       if(std::abs(resPt) > std::abs(resPtMax_)) resPtMax_ = resPt;
00452       if(std::abs(resChargedHadEnergy) > std::abs(resChargedHadEnergyMax_) ) resChargedHadEnergyMax_ = resChargedHadEnergy;
00453       if(std::abs(resNeutralHadEnergy) > std::abs(resNeutralHadEnergyMax_) ) resNeutralHadEnergyMax_ = resNeutralHadEnergy;
00454       if(std::abs(resNeutralEmEnergy) > std::abs(resNeutralEmEnergyMax_) ) resNeutralEmEnergyMax_ = resNeutralEmEnergy;
00455       if (debug_) {
00456         cout << i <<"  =========PFJet Pt "<< rec_pt
00457              << " eta " << rec_eta
00458              << " phi " << rec_phi
00459              << " Charged Had Energy " << rec_ChargedHadEnergy
00460              << " Neutral Had Energy " << rec_NeutralHadEnergy
00461              << " Neutral elm Energy " << rec_NeutralEmEnergy << endl;
00462         cout << " matching Gen Jet Pt " << true_pt
00463              << " eta " << truth->eta()
00464              << " phi " << truth->phi()
00465              << " Charged Had Energy " << true_ChargedHadEnergy
00466              << " Neutral Had Energy " << true_NeutralHadEnergy
00467              << " Neutral elm Energy " << true_NeutralEmEnergy << endl;
00468         printPFJet(&pfj);
00469         //      cout<<pfj.print()<<endl;
00470         printGenJet(truth);
00471         //cout <<truth->print()<<endl;
00472                                 
00473         cout << "==============deltaR " << deltaR << "  resPt " << resPt
00474              << " resChargedHadEnergy " << resChargedHadEnergy
00475              << " resNeutralHadEnergy " << resNeutralHadEnergy
00476              << " resNeutralEmEnergy " << resNeutralEmEnergy
00477              << endl;
00478       }
00479                         
00480 
00481       if(plot1) {
00482         if ( rec_eta > 0. ) 
00483           hDEtavsEta->Fill(true_eta,rec_eta-true_eta);
00484         else
00485           hDEtavsEta->Fill(true_eta,-rec_eta+true_eta);
00486         hDPhivsEta->Fill(true_eta,rec_phi-true_phi);
00487 
00488         hRPtvsEta->Fill(true_eta, resPt);
00489         hNCHvsEta->Fill(true_eta, rec_ChargedMultiplicity);
00490         hNCH0vsEta->Fill(true_eta,chMult[0]);
00491         hNCH1vsEta->Fill(true_eta,chMult[1]);
00492         hNCH2vsEta->Fill(true_eta,chMult[2]);
00493         hNCH3vsEta->Fill(true_eta,chMult[3]);
00494         hNCH4vsEta->Fill(true_eta,chMult[4]);
00495         hNCH5vsEta->Fill(true_eta,chMult[5]);
00496         hNCH6vsEta->Fill(true_eta,chMult[6]);
00497       }
00498       if(plot2)hRCHEvsEta->Fill(true_eta, resChargedHadEnergy);
00499       if(plot5)hRNeutvsEta->Fill(true_eta, resNeutralEnergy);
00500       if(plot6) { 
00501         hRHCALvsEta->Fill(true_eta, resHCALEnergy);
00502         hRNEUTvsEta->Fill(true_eta, resNEUTEnergy);
00503       }
00504       if(plot7) {  
00505         hRHONLvsEta->Fill(true_eta, resHCALEnergy);
00506         hRNONLvsEta->Fill(true_eta, resNEUTEnergy);
00507       }
00508 
00509       // fill histograms for relative delta quantitites of matched jets
00510       // delta Pt or E quantities for Barrel
00511       if (Barrel){
00512         if(plot1) { 
00513           hBRPt->Fill (resPt);
00514           if ( pt_denom >  20. && pt_denom <  40. ) hBRPt20_40->Fill (resPt);
00515           if ( pt_denom >  40. && pt_denom <  60. ) hBRPt40_60->Fill (resPt);
00516           if ( pt_denom >  60. && pt_denom <  80. ) hBRPt60_80->Fill (resPt);
00517           if ( pt_denom >  80. && pt_denom < 100. ) hBRPt80_100->Fill (resPt);
00518           if ( pt_denom > 100. && pt_denom < 150. ) hBRPt100_150->Fill (resPt);
00519           if ( pt_denom > 150. && pt_denom < 200. ) hBRPt150_200->Fill (resPt);
00520           if ( pt_denom > 200. && pt_denom < 250. ) hBRPt200_250->Fill (resPt);
00521           if ( pt_denom > 250. && pt_denom < 300. ) hBRPt250_300->Fill (resPt);
00522           if ( pt_denom > 300. && pt_denom < 400. ) hBRPt300_400->Fill (resPt);
00523           if ( pt_denom > 400. && pt_denom < 500. ) hBRPt400_500->Fill (resPt);
00524           if ( pt_denom > 500. && pt_denom < 750. ) hBRPt500_750->Fill (resPt);
00525           if ( pt_denom > 750. && pt_denom < 1250. ) hBRPt750_1250->Fill (resPt);
00526           if ( pt_denom > 1250. && pt_denom < 2000. ) hBRPt1250_2000->Fill (resPt);
00527           if ( pt_denom > 2000. && pt_denom < 5000. ) hBRPt2000_5000->Fill (resPt);
00528           hBNCH->Fill(rec_ChargedMultiplicity);
00529           hBNCH0vsPt->Fill(pt_denom,chMult[0]);
00530           hBNCH1vsPt->Fill(pt_denom,chMult[1]);
00531           hBNCH2vsPt->Fill(pt_denom,chMult[2]);
00532           hBNCH3vsPt->Fill(pt_denom,chMult[3]);
00533           hBNCH4vsPt->Fill(pt_denom,chMult[4]);
00534           hBNCH5vsPt->Fill(pt_denom,chMult[5]);
00535           hBNCH6vsPt->Fill(pt_denom,chMult[6]);
00536           hBNCHvsPt->Fill(pt_denom,rec_ChargedMultiplicity);
00537           if ( rec_eta > 0. ) 
00538             hBDEtavsPt->Fill(pt_denom,rec_eta-true_eta);
00539           else
00540             hBDEtavsPt->Fill(pt_denom,-rec_eta+true_eta);
00541           hBDPhivsPt->Fill(pt_denom,rec_phi-true_phi);
00542         }
00543         if(plot2)hBRCHE->Fill(resChargedHadEnergy);
00544         if(plot3)hBRNHE->Fill(resNeutralHadEnergy);
00545         if(plot4)hBRNEE->Fill(resNeutralEmEnergy);
00546         if(plot5)hBRneut->Fill(resNeutralEnergy);
00547         if(plot1)hBRPtvsPt->Fill(pt_denom, resPt);
00548         if(plot2)hBRCHEvsPt->Fill(pt_denom, resChargedHadEnergy);
00549         if(plot3)hBRNHEvsPt->Fill(pt_denom, resNeutralHadEnergy);
00550         if(plot4)hBRNEEvsPt->Fill(pt_denom, resNeutralEmEnergy);
00551         if(plot5)hBRneutvsPt->Fill(pt_denom, resNeutralEnergy);
00552         if(plot6) { 
00553           hBRHCALvsP->Fill(true_E, resHCALEnergy);
00554           hBRNEUTvsP->Fill(true_E, resNEUTEnergy);
00555         }
00556         if(plot7) { 
00557           hBRHONLvsP->Fill(true_E, resHCALEnergy);
00558           hBRNONLvsP->Fill(true_E, resNEUTEnergy);
00559         }
00560 
00561       }
00562       // delta Pt or E quantities for Endcap
00563       if (Endcap){
00564         if(plot1) {
00565           hERPt->Fill (resPt);
00566           if ( pt_denom >  20. && pt_denom <  40. ) hERPt20_40->Fill (resPt);
00567           if ( pt_denom >  40. && pt_denom <  60. ) hERPt40_60->Fill (resPt);
00568           if ( pt_denom >  60. && pt_denom <  80. ) hERPt60_80->Fill (resPt);
00569           if ( pt_denom >  80. && pt_denom < 100. ) hERPt80_100->Fill (resPt);
00570           if ( pt_denom > 100. && pt_denom < 150. ) hERPt100_150->Fill (resPt);
00571           if ( pt_denom > 150. && pt_denom < 200. ) hERPt150_200->Fill (resPt);
00572           if ( pt_denom > 200. && pt_denom < 250. ) hERPt200_250->Fill (resPt);
00573           if ( pt_denom > 250. && pt_denom < 300. ) hERPt250_300->Fill (resPt);
00574           if ( pt_denom > 300. && pt_denom < 400. ) hERPt300_400->Fill (resPt);
00575           if ( pt_denom > 400. && pt_denom < 500. ) hERPt400_500->Fill (resPt);
00576           if ( pt_denom > 500. && pt_denom < 750. ) hERPt500_750->Fill (resPt);
00577           if ( pt_denom > 750. && pt_denom < 1250. ) hERPt750_1250->Fill (resPt);
00578           if ( pt_denom > 1250. && pt_denom < 2000. ) hERPt1250_2000->Fill (resPt);
00579           if ( pt_denom > 2000. && pt_denom < 5000. ) hERPt2000_5000->Fill (resPt);
00580           hENCH->Fill(rec_ChargedMultiplicity);
00581           hENCHvsPt->Fill(pt_denom,rec_ChargedMultiplicity);
00582           hENCH0vsPt->Fill(pt_denom,chMult[0]);
00583           hENCH1vsPt->Fill(pt_denom,chMult[1]);
00584           hENCH2vsPt->Fill(pt_denom,chMult[2]);
00585           hENCH3vsPt->Fill(pt_denom,chMult[3]);
00586           hENCH4vsPt->Fill(pt_denom,chMult[4]);
00587           hENCH5vsPt->Fill(pt_denom,chMult[5]);
00588           hENCH6vsPt->Fill(pt_denom,chMult[6]);
00589           if ( rec_eta > 0. ) 
00590             hEDEtavsPt->Fill(pt_denom,rec_eta-true_eta);
00591           else
00592             hEDEtavsPt->Fill(pt_denom,-rec_eta+true_eta);
00593           hEDPhivsPt->Fill(pt_denom,rec_phi-true_phi);
00594         }
00595         if(plot2)hERCHE->Fill(resChargedHadEnergy);
00596         if(plot3)hERNHE->Fill(resNeutralHadEnergy);
00597         if(plot4)hERNEE->Fill(resNeutralEmEnergy);
00598         if(plot5)hERneut->Fill(resNeutralEnergy);
00599         if(plot1)hERPtvsPt->Fill(pt_denom, resPt);
00600         if(plot2)hERCHEvsPt->Fill(pt_denom, resChargedHadEnergy);
00601         if(plot3)hERNHEvsPt->Fill(pt_denom, resNeutralHadEnergy);
00602         if(plot4)hERNEEvsPt->Fill(pt_denom, resNeutralEmEnergy);
00603         if(plot5)hERneutvsPt->Fill(pt_denom, resNeutralEnergy);
00604         if(plot6) {
00605           hERHCALvsP->Fill(true_E, resHCALEnergy);
00606           hERNEUTvsP->Fill(true_E, resNEUTEnergy);
00607         }
00608         if(plot7) {
00609           hERHONLvsP->Fill(true_E, resHCALEnergy);
00610           hERNONLvsP->Fill(true_E, resNEUTEnergy);
00611         }
00612       }                                         
00613       // delta Pt or E quantities for Forward
00614       if (Forward){
00615         if(plot1) {
00616           hFRPt->Fill (resPt);
00617           if ( pt_denom >  20. && pt_denom <  40. ) hFRPt20_40->Fill (resPt);
00618           if ( pt_denom >  40. && pt_denom <  60. ) hFRPt40_60->Fill (resPt);
00619           if ( pt_denom >  60. && pt_denom <  80. ) hFRPt60_80->Fill (resPt);
00620           if ( pt_denom >  80. && pt_denom < 100. ) hFRPt80_100->Fill (resPt);
00621           if ( pt_denom > 100. && pt_denom < 150. ) hFRPt100_150->Fill (resPt);
00622           if ( pt_denom > 150. && pt_denom < 200. ) hFRPt150_200->Fill (resPt);
00623           if ( pt_denom > 200. && pt_denom < 250. ) hFRPt200_250->Fill (resPt);
00624           if ( pt_denom > 250. && pt_denom < 300. ) hFRPt250_300->Fill (resPt);
00625           if ( pt_denom > 300. && pt_denom < 400. ) hFRPt300_400->Fill (resPt);
00626           if ( pt_denom > 400. && pt_denom < 500. ) hFRPt400_500->Fill (resPt);
00627           if ( pt_denom > 500. && pt_denom < 750. ) hFRPt500_750->Fill (resPt);
00628           if ( pt_denom > 750. && pt_denom < 1250. ) hFRPt750_1250->Fill (resPt);
00629           if ( pt_denom > 1250. && pt_denom < 2000. ) hFRPt1250_2000->Fill (resPt);
00630           if ( pt_denom > 2000. && pt_denom < 5000. ) hFRPt2000_5000->Fill (resPt);
00631           if ( rec_eta > 0. ) 
00632             hFDEtavsPt->Fill(pt_denom,rec_eta-true_eta);
00633           else
00634             hFDEtavsPt->Fill(pt_denom,-rec_eta+true_eta);
00635           hFDPhivsPt->Fill(pt_denom,rec_phi-true_phi);
00636         }
00637         if(plot2)hFRCHE->Fill(resChargedHadEnergy);
00638         if(plot3)hFRNHE->Fill(resNeutralHadEnergy);
00639         if(plot4)hFRNEE->Fill(resNeutralEmEnergy);
00640         if(plot5)hFRneut->Fill(resNeutralEnergy);
00641         if(plot1)hFRPtvsPt->Fill(pt_denom, resPt);
00642         if(plot2)hFRCHEvsPt->Fill(pt_denom, resChargedHadEnergy);
00643         if(plot3)hFRNHEvsPt->Fill(pt_denom, resNeutralHadEnergy);
00644         if(plot4)hFRNEEvsPt->Fill(pt_denom, resNeutralEmEnergy);
00645         if(plot5)hFRneutvsPt->Fill(pt_denom, resNeutralEnergy);
00646         if(plot6) {
00647           hFRHCALvsP->Fill(true_E, resHCALEnergy);
00648           hFRNEUTvsP->Fill(true_E, resNEUTEnergy);
00649         }
00650         if(plot7) {
00651           hFRHONLvsP->Fill(true_E, resHCALEnergy);
00652           hFRNONLvsP->Fill(true_E, resNEUTEnergy);
00653         }
00654       }                                         
00655     } // end case deltaR < deltaRMax
00656                 
00657   } // i loop on pf Jets        
00658 
00659   // Increment counter
00660   entry_++;
00661 }
00662 
00663 void PFJetBenchmark::gettrue (const reco::GenJet* truth, double& true_ChargedHadEnergy, 
00664                               double& true_NeutralHadEnergy, double& true_NeutralEmEnergy){
00665   std::vector <const GenParticle*> mcparts = truth->getGenConstituents ();
00666   true_NeutralEmEnergy = 0.;
00667   true_ChargedHadEnergy = 0.;
00668   true_NeutralHadEnergy = 0.;
00669   // for each MC particle in turn  
00670   for (unsigned i = 0; i < mcparts.size (); i++) {
00671     const GenParticle* mcpart = mcparts[i];
00672     int PDG = std::abs( mcpart->pdgId());
00673     double e = mcpart->energy(); 
00674     switch(PDG){  // start PDG switch
00675     case 22: // photon
00676       true_NeutralEmEnergy += e;
00677       break;
00678     case 211: // pi
00679     case 321: // K
00680     case 2212: // p
00681     case 11: //electrons (until recognised)
00682       true_ChargedHadEnergy += e;
00683       break;
00684     case 310: // K_S0
00685     case 130: // K_L0
00686     case 3122: // Lambda0
00687     case 2112: // n0
00688       true_NeutralHadEnergy += e;
00689     default:
00690       break;
00691     }  // end PDG switch                
00692   }  // end loop on constituents.
00693 }
00694 
00695 void PFJetBenchmark::printPFJet(const reco::PFJet* pfj){
00696   cout<<setiosflags(ios::right);
00697   cout<<setiosflags(ios::fixed);
00698   cout<<setprecision(3);
00699 
00700   cout << "PFJet  p/px/py/pz/pt: " << pfj->p() << "/" << pfj->px () 
00701        << "/" << pfj->py() << "/" << pfj->pz() << "/" << pfj->pt() << endl
00702        << "    eta/phi: " << pfj->eta () << "/" << pfj->phi () << endl                  
00703        << "    PFJet specific:" << std::endl
00704        << "      charged/neutral hadrons energy: " << pfj->chargedHadronEnergy () << '/' << pfj->neutralHadronEnergy () << endl
00705        << "      charged/neutral em energy: " << pfj->chargedEmEnergy () << '/' << pfj->neutralEmEnergy () << endl
00706        << "      charged muon energy: " << pfj->chargedMuEnergy () << '/' << endl
00707        << "      charged/neutral multiplicity: " << pfj->chargedMultiplicity () << '/' << pfj->neutralMultiplicity () << endl;
00708   
00709   // And print the constituents
00710   std::cout << pfj->print() << std::endl;
00711 
00712   cout<<resetiosflags(ios::right|ios::fixed);
00713 }
00714 
00715 
00716 void PFJetBenchmark::printGenJet (const reco::GenJet* truth){
00717   std::vector <const GenParticle*> mcparts = truth->getGenConstituents ();
00718   cout << "GenJet p/px/py/pz/pt: " << truth->p() << '/' << truth->px () 
00719        << '/' << truth->py() << '/' << truth->pz() << '/' << truth->pt() << endl
00720        << "    eta/phi: " << truth->eta () << '/' << truth->phi () << endl
00721        << "    # of constituents: " << mcparts.size() << endl;
00722   cout << "    constituents:" << endl;
00723   for (unsigned i = 0; i < mcparts.size (); i++) {
00724     const GenParticle* mcpart = mcparts[i];
00725     cout << "      #" << i << "  PDG code:" << mcpart->pdgId() 
00726          << ", p/pt/eta/phi: " << mcpart->p() << '/' << mcpart->pt() 
00727          << '/' << mcpart->eta() << '/' << mcpart->phi() << endl;       
00728   }    
00729 }
00730
/data/refman/pasoursint/CMSSW_4_2_9_HLT1_bphpatch4/src/RecoParticleFlow/Benchmark/src/PFJetBenchmark.cc
