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) > 3.3 && std::abs (rec_eta) < 4.7 ) Forward = true;
00270 
00271     // do only barrel for now
00272     //  if(!Barrel) continue;
00273 
00274     // look for the closets gen Jet : truth
00275     const GenJet *truth = algo_->matchByDeltaR(&pfj,&genJets);
00276     if(!truth) continue;   
00277     double deltaR = algo_->deltaR(&pfj, truth);
00278     // check deltaR is small enough
00279     if(deltaR < deltaRMax_ || deltaRMax_ == -1.0 ) {//start case deltaR < deltaRMax
00280       // generate histograms comparing the reco and truth candidate (truth = closest in delta-R) 
00281       // get the quantities to place on the denominator and/or divide by
00282       double pt_denom;
00283       double true_E = truth->p();
00284       double true_pt = truth->pt();
00285       double true_eta = truth->eta();
00286       double true_phi = truth->phi();
00287 
00288       if (plotAgainstReco_) {pt_denom = rec_pt;}
00289       else {pt_denom = true_pt;}
00290       // get true specific quantities
00291       double true_ChargedHadEnergy;
00292       double true_NeutralHadEnergy;
00293       double true_NeutralEmEnergy;
00294       gettrue (truth, true_ChargedHadEnergy, true_NeutralHadEnergy, true_NeutralEmEnergy);
00295       double true_NeutralEnergy = true_NeutralHadEnergy + true_NeutralEmEnergy;
00296       double rec_ChargedHadEnergy = pfj.chargedHadronEnergy();
00297       double rec_NeutralHadEnergy = pfj.neutralHadronEnergy();
00298       double rec_NeutralEmEnergy = pfj.neutralEmEnergy();
00299       double rec_NeutralEnergy = rec_NeutralHadEnergy + rec_NeutralEmEnergy;
00300       double rec_ChargedMultiplicity = pfj.chargedMultiplicity();
00301       std::vector <PFCandidatePtr> constituents = pfj.getPFConstituents ();
00302       std::vector <unsigned int> chMult(7, static_cast<unsigned int>(0)); 
00303       for (unsigned ic = 0; ic < constituents.size (); ++ic) {
00304         if ( constituents[ic]->particleId() > 3 ) continue;
00305         reco::TrackRef trackRef = constituents[ic]->trackRef();
00306         if ( trackRef.isNull() ) {
00307           //std::cout << "Warning in entry " << entry_ 
00308           //        << " : a track with Id " << constituents[ic]->particleId() 
00309           //        << " has no track ref.." << std::endl;
00310           continue;
00311         }
00312         unsigned int iter = 0; 
00313         switch (trackRef->algo()) {
00314         case TrackBase::ctf:
00315         case TrackBase::iter0:
00316           iter = 0;
00317           break;
00318         case TrackBase::iter1:
00319           iter = 1;
00320           break;
00321         case TrackBase::iter2:
00322           iter = 2;
00323           break;
00324         case TrackBase::iter3:
00325           iter = 3;
00326           break;
00327         case TrackBase::iter4:
00328           iter = 4;
00329           break;
00330         case TrackBase::iter5:
00331           iter = 5;
00332           break;
00333         default:
00334           iter = 6;
00335           std::cout << "Warning in entry " << entry_ << " : iter = 6... " << std::endl;
00336           break;
00337         }
00338         ++(chMult[iter]);
00339       }
00340 
00341       bool plot1 = false;
00342       bool plot2 = false;
00343       bool plot3 = false;
00344       bool plot4 = false;
00345       bool plot5 = false;
00346       bool plot6 = false;
00347       bool plot7 = false;
00348       double cut1 = 0.0001;
00349       double cut2 = 0.0001;
00350       double cut3 = 0.0001;
00351       double cut4 = 0.0001;
00352       double cut5 = 0.0001;
00353       double cut6 = 0.0001;
00354       double cut7 = 0.0001;
00355       double resPt =0.;
00356       double resChargedHadEnergy= 0.;
00357       double resNeutralHadEnergy= 0.;
00358       double resNeutralEmEnergy= 0.;
00359       double resNeutralEnergy= 0.;
00360 
00361       double resHCALEnergy = 0.;
00362       double resNEUTEnergy = 0.;
00363       if ( rec_NeutralHadEnergy > cut6 && rec_ChargedHadEnergy < cut1 ) { 
00364         double true_NEUTEnergy = true_NeutralHadEnergy + true_NeutralEmEnergy;
00365         double true_HCALEnergy = true_NEUTEnergy - rec_NeutralEmEnergy;
00366         double rec_NEUTEnergy = rec_NeutralHadEnergy+rec_NeutralEmEnergy; 
00367         double rec_HCALEnergy = rec_NeutralHadEnergy; 
00368         resHCALEnergy = (rec_HCALEnergy-true_HCALEnergy)/rec_HCALEnergy;
00369         resNEUTEnergy = (rec_NEUTEnergy-true_NEUTEnergy)/rec_NEUTEnergy;
00370         if ( rec_NeutralEmEnergy > cut7 ) {
00371           plot6 = true;
00372         } else {
00373           plot7 = true;
00374         }
00375       }
00376 
00377       // get relative delta quantities (protect against division by zero!)
00378       if (true_pt > 0.0001){
00379         resPt = (rec_pt -true_pt)/true_pt ; 
00380         plot1 = true;}
00381       if (true_ChargedHadEnergy > cut1){
00382         resChargedHadEnergy = (rec_ChargedHadEnergy- true_ChargedHadEnergy)/true_ChargedHadEnergy;
00383         plot2 = true;}
00384       if (true_NeutralHadEnergy > cut2){
00385         resNeutralHadEnergy = (rec_NeutralHadEnergy- true_NeutralHadEnergy)/true_NeutralHadEnergy;
00386         plot3 = true;}
00387       else 
00388         if (rec_NeutralHadEnergy > cut3){
00389           resNeutralHadEnergy = (rec_NeutralHadEnergy- true_NeutralHadEnergy)/rec_NeutralHadEnergy;
00390           plot3 = true;}
00391       if (true_NeutralEmEnergy > cut4){
00392         resNeutralEmEnergy = (rec_NeutralEmEnergy- true_NeutralEmEnergy)/true_NeutralEmEnergy;
00393         plot4 = true;}
00394       if (true_NeutralEnergy > cut5){
00395         resNeutralEnergy = (rec_NeutralEnergy- true_NeutralEnergy)/true_NeutralEnergy;
00396         plot5 = true;}
00397       
00398       //double deltaEta = algo_->deltaEta(&pfj, truth);
00399       //double deltaPhi = algo_->deltaPhi(&pfj, truth);
00400 
00401       // Print outliers for further debugging
00402       if ( ( resPt > 0.2 && true_pt > 100. ) || 
00403            ( resPt < -0.5 && true_pt > 100. ) ) {
00404         //if ( ( true_pt > 50. && 
00405         //     ( ( truth->eta()>3.0 && rec_eta-truth->eta() < -0.1 ) || 
00406         //       ( truth->eta()<-3.0 && rec_eta-truth->eta() > 0.1 ) ))) {
00407         std::cout << "Entry " << entry_ 
00408                   << " resPt = " << resPt
00409                   <<" resCharged  " << resChargedHadEnergy
00410                   <<" resNeutralHad  " << resNeutralHadEnergy
00411                   << " resNeutralEm  " << resNeutralEmEnergy
00412                   << " pT (T/R) " << true_pt << "/" << rec_pt 
00413                   << " Eta (T/R) " << truth->eta() << "/" << rec_eta 
00414                   << " Phi (T/R) " << truth->phi() << "/" << rec_phi 
00415                   << std::endl;
00416 
00417         // check overlapping PF jets
00418         const reco::PFJet* pfoj = 0; 
00419         double dRo = 1E9;
00420         for(unsigned j=0; j<pfJets.size(); j++) { 
00421           const reco::PFJet& pfo = pfJets[j];
00422           if ( j != i &&  algo_->deltaR(&pfj,&pfo) < dRo && pfo.pt() > 0.25*pfj.pt()) { 
00423             dRo = algo_->deltaR(&pfj,&pfo);     
00424             pfoj = &pfo;
00425           }
00426         }
00427         
00428         // Check overlapping Gen Jet 
00429         math::XYZTLorentzVector overlappinGenJet(0.,0.,0.,0.);
00430         const reco::GenJet* genoj = 0;
00431         double dRgo = 1E9;
00432         for(unsigned j=0; j<genJets.size(); j++) { 
00433           const reco::GenJet* gjo = &(genJets[j]);
00434           if ( gjo != truth && algo_->deltaR(truth,gjo) < dRgo && gjo->pt() > 0.25*truth->pt() ) { 
00435             dRgo = algo_->deltaR(truth,gjo);
00436             genoj = gjo;
00437           }
00438         }
00439         
00440         if ( dRo < 0.8 && dRgo < 0.8 && algo_->deltaR(genoj,pfoj) < 2.*deltaRMax_ ) 
00441           std::cout << "Excess probably due to overlapping jets (DR = " <<   algo_->deltaR(genoj,pfoj) << "),"
00442                     << " at DeltaR(T/R) = " << dRgo << "/" << dRo  
00443                     << " with pT(T/R) " << genoj->pt() << "/" << pfoj->pt()
00444                     << " and Eta (T/R) " << genoj->eta() << "/" << pfoj->eta()
00445                     << " and Phi (T/R) " << genoj->phi() << "/" << pfoj->phi()
00446                     << std::endl;
00447       }
00448 
00449       if(std::abs(resPt) > std::abs(resPtMax_)) resPtMax_ = resPt;
00450       if(std::abs(resChargedHadEnergy) > std::abs(resChargedHadEnergyMax_) ) resChargedHadEnergyMax_ = resChargedHadEnergy;
00451       if(std::abs(resNeutralHadEnergy) > std::abs(resNeutralHadEnergyMax_) ) resNeutralHadEnergyMax_ = resNeutralHadEnergy;
00452       if(std::abs(resNeutralEmEnergy) > std::abs(resNeutralEmEnergyMax_) ) resNeutralEmEnergyMax_ = resNeutralEmEnergy;
00453       if (debug_) {
00454         cout << i <<"  =========PFJet Pt "<< rec_pt
00455              << " eta " << rec_eta
00456              << " phi " << rec_phi
00457              << " Charged Had Energy " << rec_ChargedHadEnergy
00458              << " Neutral Had Energy " << rec_NeutralHadEnergy
00459              << " Neutral elm Energy " << rec_NeutralEmEnergy << endl;
00460         cout << " matching Gen Jet Pt " << true_pt
00461              << " eta " << truth->eta()
00462              << " phi " << truth->phi()
00463              << " Charged Had Energy " << true_ChargedHadEnergy
00464              << " Neutral Had Energy " << true_NeutralHadEnergy
00465              << " Neutral elm Energy " << true_NeutralEmEnergy << endl;
00466         printPFJet(&pfj);
00467         //      cout<<pfj.print()<<endl;
00468         printGenJet(truth);
00469         //cout <<truth->print()<<endl;
00470                                 
00471         cout << "==============deltaR " << deltaR << "  resPt " << resPt
00472              << " resChargedHadEnergy " << resChargedHadEnergy
00473              << " resNeutralHadEnergy " << resNeutralHadEnergy
00474              << " resNeutralEmEnergy " << resNeutralEmEnergy
00475              << endl;
00476       }
00477                         
00478 
00479       if(plot1) {
00480         if ( rec_eta > 0. ) 
00481           hDEtavsEta->Fill(true_eta,rec_eta-true_eta);
00482         else
00483           hDEtavsEta->Fill(true_eta,-rec_eta+true_eta);
00484         hDPhivsEta->Fill(true_eta,rec_phi-true_phi);
00485 
00486         hRPtvsEta->Fill(true_eta, resPt);
00487         hNCHvsEta->Fill(true_eta, rec_ChargedMultiplicity);
00488         hNCH0vsEta->Fill(true_eta,chMult[0]);
00489         hNCH1vsEta->Fill(true_eta,chMult[1]);
00490         hNCH2vsEta->Fill(true_eta,chMult[2]);
00491         hNCH3vsEta->Fill(true_eta,chMult[3]);
00492         hNCH4vsEta->Fill(true_eta,chMult[4]);
00493         hNCH5vsEta->Fill(true_eta,chMult[5]);
00494         hNCH6vsEta->Fill(true_eta,chMult[6]);
00495       }
00496       if(plot2)hRCHEvsEta->Fill(true_eta, resChargedHadEnergy);
00497       if(plot5)hRNeutvsEta->Fill(true_eta, resNeutralEnergy);
00498       if(plot6) { 
00499         hRHCALvsEta->Fill(true_eta, resHCALEnergy);
00500         hRNEUTvsEta->Fill(true_eta, resNEUTEnergy);
00501       }
00502       if(plot7) {  
00503         hRHONLvsEta->Fill(true_eta, resHCALEnergy);
00504         hRNONLvsEta->Fill(true_eta, resNEUTEnergy);
00505       }
00506 
00507       // fill histograms for relative delta quantitites of matched jets
00508       // delta Pt or E quantities for Barrel
00509       if (Barrel){
00510         if(plot1) { 
00511           hBRPt->Fill (resPt);
00512           if ( pt_denom >  20. && pt_denom <  40. ) hBRPt20_40->Fill (resPt);
00513           if ( pt_denom >  40. && pt_denom <  60. ) hBRPt40_60->Fill (resPt);
00514           if ( pt_denom >  60. && pt_denom <  80. ) hBRPt60_80->Fill (resPt);
00515           if ( pt_denom >  80. && pt_denom < 100. ) hBRPt80_100->Fill (resPt);
00516           if ( pt_denom > 100. && pt_denom < 150. ) hBRPt100_150->Fill (resPt);
00517           if ( pt_denom > 150. && pt_denom < 200. ) hBRPt150_200->Fill (resPt);
00518           if ( pt_denom > 200. && pt_denom < 250. ) hBRPt200_250->Fill (resPt);
00519           if ( pt_denom > 250. && pt_denom < 300. ) hBRPt250_300->Fill (resPt);
00520           if ( pt_denom > 300. && pt_denom < 400. ) hBRPt300_400->Fill (resPt);
00521           if ( pt_denom > 400. && pt_denom < 500. ) hBRPt400_500->Fill (resPt);
00522           if ( pt_denom > 500. && pt_denom < 750. ) hBRPt500_750->Fill (resPt);
00523           if ( pt_denom > 750. && pt_denom < 1250. ) hBRPt750_1250->Fill (resPt);
00524           if ( pt_denom > 1250. && pt_denom < 2000. ) hBRPt1250_2000->Fill (resPt);
00525           if ( pt_denom > 2000. && pt_denom < 5000. ) hBRPt2000_5000->Fill (resPt);
00526           hBNCH->Fill(rec_ChargedMultiplicity);
00527           hBNCH0vsPt->Fill(pt_denom,chMult[0]);
00528           hBNCH1vsPt->Fill(pt_denom,chMult[1]);
00529           hBNCH2vsPt->Fill(pt_denom,chMult[2]);
00530           hBNCH3vsPt->Fill(pt_denom,chMult[3]);
00531           hBNCH4vsPt->Fill(pt_denom,chMult[4]);
00532           hBNCH5vsPt->Fill(pt_denom,chMult[5]);
00533           hBNCH6vsPt->Fill(pt_denom,chMult[6]);
00534           hBNCHvsPt->Fill(pt_denom,rec_ChargedMultiplicity);
00535           if ( rec_eta > 0. ) 
00536             hBDEtavsPt->Fill(pt_denom,rec_eta-true_eta);
00537           else
00538             hBDEtavsPt->Fill(pt_denom,-rec_eta+true_eta);
00539           hBDPhivsPt->Fill(pt_denom,rec_phi-true_phi);
00540         }
00541         if(plot2)hBRCHE->Fill(resChargedHadEnergy);
00542         if(plot3)hBRNHE->Fill(resNeutralHadEnergy);
00543         if(plot4)hBRNEE->Fill(resNeutralEmEnergy);
00544         if(plot5)hBRneut->Fill(resNeutralEnergy);
00545         if(plot1)hBRPtvsPt->Fill(pt_denom, resPt);
00546         if(plot2)hBRCHEvsPt->Fill(pt_denom, resChargedHadEnergy);
00547         if(plot3)hBRNHEvsPt->Fill(pt_denom, resNeutralHadEnergy);
00548         if(plot4)hBRNEEvsPt->Fill(pt_denom, resNeutralEmEnergy);
00549         if(plot5)hBRneutvsPt->Fill(pt_denom, resNeutralEnergy);
00550         if(plot6) { 
00551           hBRHCALvsP->Fill(true_E, resHCALEnergy);
00552           hBRNEUTvsP->Fill(true_E, resNEUTEnergy);
00553         }
00554         if(plot7) { 
00555           hBRHONLvsP->Fill(true_E, resHCALEnergy);
00556           hBRNONLvsP->Fill(true_E, resNEUTEnergy);
00557         }
00558 
00559       }
00560       // delta Pt or E quantities for Endcap
00561       if (Endcap){
00562         if(plot1) {
00563           hERPt->Fill (resPt);
00564           if ( pt_denom >  20. && pt_denom <  40. ) hERPt20_40->Fill (resPt);
00565           if ( pt_denom >  40. && pt_denom <  60. ) hERPt40_60->Fill (resPt);
00566           if ( pt_denom >  60. && pt_denom <  80. ) hERPt60_80->Fill (resPt);
00567           if ( pt_denom >  80. && pt_denom < 100. ) hERPt80_100->Fill (resPt);
00568           if ( pt_denom > 100. && pt_denom < 150. ) hERPt100_150->Fill (resPt);
00569           if ( pt_denom > 150. && pt_denom < 200. ) hERPt150_200->Fill (resPt);
00570           if ( pt_denom > 200. && pt_denom < 250. ) hERPt200_250->Fill (resPt);
00571           if ( pt_denom > 250. && pt_denom < 300. ) hERPt250_300->Fill (resPt);
00572           if ( pt_denom > 300. && pt_denom < 400. ) hERPt300_400->Fill (resPt);
00573           if ( pt_denom > 400. && pt_denom < 500. ) hERPt400_500->Fill (resPt);
00574           if ( pt_denom > 500. && pt_denom < 750. ) hERPt500_750->Fill (resPt);
00575           if ( pt_denom > 750. && pt_denom < 1250. ) hERPt750_1250->Fill (resPt);
00576           if ( pt_denom > 1250. && pt_denom < 2000. ) hERPt1250_2000->Fill (resPt);
00577           if ( pt_denom > 2000. && pt_denom < 5000. ) hERPt2000_5000->Fill (resPt);
00578           hENCH->Fill(rec_ChargedMultiplicity);
00579           hENCHvsPt->Fill(pt_denom,rec_ChargedMultiplicity);
00580           hENCH0vsPt->Fill(pt_denom,chMult[0]);
00581           hENCH1vsPt->Fill(pt_denom,chMult[1]);
00582           hENCH2vsPt->Fill(pt_denom,chMult[2]);
00583           hENCH3vsPt->Fill(pt_denom,chMult[3]);
00584           hENCH4vsPt->Fill(pt_denom,chMult[4]);
00585           hENCH5vsPt->Fill(pt_denom,chMult[5]);
00586           hENCH6vsPt->Fill(pt_denom,chMult[6]);
00587           if ( rec_eta > 0. ) 
00588             hEDEtavsPt->Fill(pt_denom,rec_eta-true_eta);
00589           else
00590             hEDEtavsPt->Fill(pt_denom,-rec_eta+true_eta);
00591           hEDPhivsPt->Fill(pt_denom,rec_phi-true_phi);
00592         }
00593         if(plot2)hERCHE->Fill(resChargedHadEnergy);
00594         if(plot3)hERNHE->Fill(resNeutralHadEnergy);
00595         if(plot4)hERNEE->Fill(resNeutralEmEnergy);
00596         if(plot5)hERneut->Fill(resNeutralEnergy);
00597         if(plot1)hERPtvsPt->Fill(pt_denom, resPt);
00598         if(plot2)hERCHEvsPt->Fill(pt_denom, resChargedHadEnergy);
00599         if(plot3)hERNHEvsPt->Fill(pt_denom, resNeutralHadEnergy);
00600         if(plot4)hERNEEvsPt->Fill(pt_denom, resNeutralEmEnergy);
00601         if(plot5)hERneutvsPt->Fill(pt_denom, resNeutralEnergy);
00602         if(plot6) {
00603           hERHCALvsP->Fill(true_E, resHCALEnergy);
00604           hERNEUTvsP->Fill(true_E, resNEUTEnergy);
00605         }
00606         if(plot7) {
00607           hERHONLvsP->Fill(true_E, resHCALEnergy);
00608           hERNONLvsP->Fill(true_E, resNEUTEnergy);
00609         }
00610       }                                         
00611       // delta Pt or E quantities for Forward
00612       if (Forward){
00613         if(plot1) {
00614           hFRPt->Fill (resPt);
00615           if ( pt_denom >  20. && pt_denom <  40. ) hFRPt20_40->Fill (resPt);
00616           if ( pt_denom >  40. && pt_denom <  60. ) hFRPt40_60->Fill (resPt);
00617           if ( pt_denom >  60. && pt_denom <  80. ) hFRPt60_80->Fill (resPt);
00618           if ( pt_denom >  80. && pt_denom < 100. ) hFRPt80_100->Fill (resPt);
00619           if ( pt_denom > 100. && pt_denom < 150. ) hFRPt100_150->Fill (resPt);
00620           if ( pt_denom > 150. && pt_denom < 200. ) hFRPt150_200->Fill (resPt);
00621           if ( pt_denom > 200. && pt_denom < 250. ) hFRPt200_250->Fill (resPt);
00622           if ( pt_denom > 250. && pt_denom < 300. ) hFRPt250_300->Fill (resPt);
00623           if ( pt_denom > 300. && pt_denom < 400. ) hFRPt300_400->Fill (resPt);
00624           if ( pt_denom > 400. && pt_denom < 500. ) hFRPt400_500->Fill (resPt);
00625           if ( pt_denom > 500. && pt_denom < 750. ) hFRPt500_750->Fill (resPt);
00626           if ( pt_denom > 750. && pt_denom < 1250. ) hFRPt750_1250->Fill (resPt);
00627           if ( pt_denom > 1250. && pt_denom < 2000. ) hFRPt1250_2000->Fill (resPt);
00628           if ( pt_denom > 2000. && pt_denom < 5000. ) hFRPt2000_5000->Fill (resPt);
00629           if ( rec_eta > 0. ) 
00630             hFDEtavsPt->Fill(pt_denom,rec_eta-true_eta);
00631           else
00632             hFDEtavsPt->Fill(pt_denom,-rec_eta+true_eta);
00633           hFDPhivsPt->Fill(pt_denom,rec_phi-true_phi);
00634         }
00635         if(plot2)hFRCHE->Fill(resChargedHadEnergy);
00636         if(plot3)hFRNHE->Fill(resNeutralHadEnergy);
00637         if(plot4)hFRNEE->Fill(resNeutralEmEnergy);
00638         if(plot5)hFRneut->Fill(resNeutralEnergy);
00639         if(plot1)hFRPtvsPt->Fill(pt_denom, resPt);
00640         if(plot2)hFRCHEvsPt->Fill(pt_denom, resChargedHadEnergy);
00641         if(plot3)hFRNHEvsPt->Fill(pt_denom, resNeutralHadEnergy);
00642         if(plot4)hFRNEEvsPt->Fill(pt_denom, resNeutralEmEnergy);
00643         if(plot5)hFRneutvsPt->Fill(pt_denom, resNeutralEnergy);
00644         if(plot6) {
00645           hFRHCALvsP->Fill(true_E, resHCALEnergy);
00646           hFRNEUTvsP->Fill(true_E, resNEUTEnergy);
00647         }
00648         if(plot7) {
00649           hFRHONLvsP->Fill(true_E, resHCALEnergy);
00650           hFRNONLvsP->Fill(true_E, resNEUTEnergy);
00651         }
00652       }                                         
00653     } // end case deltaR < deltaRMax
00654                 
00655   } // i loop on pf Jets        
00656 
00657   // Increment counter
00658   entry_++;
00659 }
00660 
00661 void PFJetBenchmark::gettrue (const reco::GenJet* truth, double& true_ChargedHadEnergy, 
00662                               double& true_NeutralHadEnergy, double& true_NeutralEmEnergy){
00663   std::vector <const GenParticle*> mcparts = truth->getGenConstituents ();
00664   true_NeutralEmEnergy = 0.;
00665   true_ChargedHadEnergy = 0.;
00666   true_NeutralHadEnergy = 0.;
00667   // for each MC particle in turn  
00668   for (unsigned i = 0; i < mcparts.size (); i++) {
00669     const GenParticle* mcpart = mcparts[i];
00670     int PDG = std::abs( mcpart->pdgId());
00671     double e = mcpart->energy(); 
00672     switch(PDG){  // start PDG switch
00673     case 22: // photon
00674       true_NeutralEmEnergy += e;
00675       break;
00676     case 211: // pi
00677     case 321: // K
00678     case 2212: // p
00679     case 11: //electrons (until recognised)
00680       true_ChargedHadEnergy += e;
00681       break;
00682     case 310: // K_S0
00683     case 130: // K_L0
00684     case 3122: // Lambda0
00685     case 2112: // n0
00686       true_NeutralHadEnergy += e;
00687     default:
00688       break;
00689     }  // end PDG switch                
00690   }  // end loop on constituents.
00691 }
00692 
00693 void PFJetBenchmark::printPFJet(const reco::PFJet* pfj){
00694   cout<<setiosflags(ios::right);
00695   cout<<setiosflags(ios::fixed);
00696   cout<<setprecision(3);
00697 
00698   cout << "PFJet  p/px/py/pz/pt: " << pfj->p() << "/" << pfj->px () 
00699        << "/" << pfj->py() << "/" << pfj->pz() << "/" << pfj->pt() << endl
00700        << "    eta/phi: " << pfj->eta () << "/" << pfj->phi () << endl                  
00701        << "    PFJet specific:" << std::endl
00702        << "      charged/neutral hadrons energy: " << pfj->chargedHadronEnergy () << '/' << pfj->neutralHadronEnergy () << endl
00703        << "      charged/neutral em energy: " << pfj->chargedEmEnergy () << '/' << pfj->neutralEmEnergy () << endl
00704        << "      charged muon energy: " << pfj->chargedMuEnergy () << '/' << endl
00705        << "      charged/neutral multiplicity: " << pfj->chargedMultiplicity () << '/' << pfj->neutralMultiplicity () << endl;
00706   
00707   // And print the constituents
00708   std::cout << pfj->print() << std::endl;
00709 
00710   cout<<resetiosflags(ios::right|ios::fixed);
00711 }
00712 
00713 
00714 void PFJetBenchmark::printGenJet (const reco::GenJet* truth){
00715   std::vector <const GenParticle*> mcparts = truth->getGenConstituents ();
00716   cout << "GenJet p/px/py/pz/pt: " << truth->p() << '/' << truth->px () 
00717        << '/' << truth->py() << '/' << truth->pz() << '/' << truth->pt() << endl
00718        << "    eta/phi: " << truth->eta () << '/' << truth->phi () << endl
00719        << "    # of constituents: " << mcparts.size() << endl;
00720   cout << "    constituents:" << endl;
00721   for (unsigned i = 0; i < mcparts.size (); i++) {
00722     const GenParticle* mcpart = mcparts[i];
00723     cout << "      #" << i << "  PDG code:" << mcpart->pdgId() 
00724          << ", p/pt/eta/phi: " << mcpart->p() << '/' << mcpart->pt() 
00725          << '/' << mcpart->eta() << '/' << mcpart->phi() << endl;       
00726   }    
00727 }
00728
/data/refman/pasoursint/CMSSW_4_1_8_patch9/src/RecoParticleFlow/Benchmark/src/PFJetBenchmark.cc
