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