CMS 3D CMS Logo

/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/RecoParticleFlow/PFRootEvent/src/PFRootEventManager.cc

Go to the documentation of this file.
00001 #include "DataFormats/Common/interface/OrphanHandle.h"
00002 #include "DataFormats/Provenance/interface/ProductID.h"
00003 #include "DataFormats/Common/interface/RefToPtr.h"
00004 
00005 #include "DataFormats/Math/interface/LorentzVector.h"
00006 #include "DataFormats/Math/interface/Point3D.h"
00007 
00008 #include "DataFormats/ParticleFlowReco/interface/PFLayer.h"
00009 
00010 #include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
00011 #include "DataFormats/ParticleFlowReco/interface/PFBlockElement.h"
00012 
00013 #include "DataFormats/FWLite/interface/ChainEvent.h"
00014 
00015 #include "RecoParticleFlow/PFClusterProducer/interface/PFClusterAlgo.h"
00016 #include "RecoParticleFlow/PFProducer/interface/PFGeometry.h"
00017 
00018 #include "RecoParticleFlow/PFRootEvent/interface/PFRootEventManager.h"
00019 
00020 #include "RecoParticleFlow/PFRootEvent/interface/IO.h"
00021 
00022 #include "RecoParticleFlow/PFRootEvent/interface/PFJetAlgorithm.h" 
00023 #include "RecoParticleFlow/PFRootEvent/interface/JetMaker.h"
00024 
00025 #include "RecoParticleFlow/PFRootEvent/interface/Utils.h" 
00026 #include "RecoParticleFlow/PFRootEvent/interface/EventColin.h" 
00027 #include "RecoParticleFlow/PFRootEvent/interface/METManager.h"
00028 
00029 #include "RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibration.h"
00030 #include "RecoParticleFlow/PFClusterTools/interface/PFSCEnergyCalibration.h"
00031 #include "RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibrationHF.h"
00032 #include "RecoParticleFlow/PFClusterTools/interface/PFClusterCalibration.h"
00033 #include "RecoParticleFlow/PFClusterTools/interface/PFEnergyResolution.h"
00034 
00035 #include "FWCore/ServiceRegistry/interface/Service.h" 
00036 #include "DQMServices/Core/interface/DQMStore.h"
00037 #include "DQMServices/Core/interface/MonitorElement.h"
00038 
00039 #include "FWCore/FWLite/interface/AutoLibraryLoader.h"
00040 
00041 #include "TFile.h"
00042 #include "TTree.h"
00043 #include "TBranch.h"
00044 #include "TCutG.h"
00045 #include "TVector3.h"
00046 #include "TROOT.h"
00047 
00048 #include <iostream>
00049 #include <vector>
00050 #include <cstdlib>
00051 
00052 using namespace std;
00053 // using namespace boost;
00054 using namespace reco;
00055 
00056 using boost::shared_ptr;
00057 
00058 PFRootEventManager::PFRootEventManager() {}
00059 
00060 
00061 
00062 PFRootEventManager::PFRootEventManager(const char* file)
00063   : 
00064   iEvent_(0),
00065   options_(0),
00066   ev_(0),
00067   tree_(0),
00068   outTree_(0),
00069   outEvent_(0),
00070   //   clusters_(new reco::PFClusterCollection),
00071   eventAuxiliary_( new edm::EventAuxiliary ),
00072   clustersECAL_(new reco::PFClusterCollection),
00073   clustersHCAL_(new reco::PFClusterCollection),
00074   clustersHFEM_(new reco::PFClusterCollection),
00075   clustersHFHAD_(new reco::PFClusterCollection),
00076   clustersPS_(new reco::PFClusterCollection),
00077   pfBlocks_(new reco::PFBlockCollection),
00078   pfCandidates_(new reco::PFCandidateCollection),
00079   //pfJets_(new reco::PFJetCollection),
00080   outFile_(0),
00081   calibFile_(0)
00082 {
00083   
00084   
00085   //   iEvent_=0;
00086   h_deltaETvisible_MCEHT_ 
00087     = new TH1F("h_deltaETvisible_MCEHT","Jet Et difference CaloTowers-MC"
00088                ,1000,-50.,50.);
00089   h_deltaETvisible_MCPF_  
00090     = new TH1F("h_deltaETvisible_MCPF" ,"Jet Et difference ParticleFlow-MC"
00091                ,1000,-50.,50.);
00092 
00093   readOptions(file, true, true);
00094  
00095   initializeEventInformation();
00096        
00097   //   maxERecHitEcal_ = -1;
00098   //   maxERecHitHcal_ = -1;
00099 
00100 }
00101 
00102 
00103 void PFRootEventManager::initializeEventInformation() {
00104 
00105   unsigned int nev = ev_->size();
00106   for ( unsigned int entry = 0; entry < nev; ++entry ) { 
00107     ev_->to(entry);
00108     const edm::EventBase& iEv = *ev_;
00109     mapEventToEntry_[iEv.id().run()][iEv.id().luminosityBlock()][iEv.id().event()] = entry;
00110   }
00111 
00112   cout<<"Number of events: "<< nev
00113       <<" starting with event: "<<mapEventToEntry_.begin()->first<<endl;
00114 }
00115 
00116 
00117 void PFRootEventManager::reset() { 
00118 
00119   if(outEvent_) {
00120     outEvent_->reset();
00121     outTree_->GetBranch("event")->SetAddress(&outEvent_);
00122   }  
00123 
00124   if ( ev_ && ev_->isValid() ) 
00125     ev_->getTFile()->cd();
00126 }
00127 
00128 
00129 
00130 void PFRootEventManager::readOptions(const char* file, 
00131                                      bool refresh, 
00132                                      bool reconnect) {
00133                                      
00134   readSpecificOptions(file);
00135   
00136   cout<<"calling PFRootEventManager::readOptions"<<endl;
00137    
00138 
00139   reset();
00140   
00141   PFGeometry pfGeometry; // initialize geometry
00142   
00143   // cout<<"reading options "<<endl;
00144 
00145   try {
00146     if( !options_ )
00147       options_ = new IO(file);
00148     else if( refresh) {
00149       delete options_;
00150       options_ = new IO(file);
00151     }
00152   }
00153   catch( const string& err ) {
00154     cout<<err<<endl;
00155     return;
00156   }
00157 
00158 
00159   debug_ = false; 
00160   options_->GetOpt("rootevent", "debug", debug_);
00161 
00162   
00163   // output text file for calibration
00164   string calibfilename;
00165   options_->GetOpt("calib","outfile",calibfilename);
00166   if (!calibfilename.empty()) { 
00167     calibFile_ = new std::ofstream(calibfilename.c_str());
00168     std::cout << "Calib file name " << calibfilename << " " << calibfilename.c_str() << std::endl;
00169   }
00170 
00171   // output root file   ------------------------------------------
00172 
00173   
00174   if(!outFile_) {
00175     string outfilename;
00176     options_->GetOpt("root","outfile", outfilename);
00177     bool doOutTree = false;
00178     options_->GetOpt("root","outtree", doOutTree);
00179     if(doOutTree) {
00180       if(!outfilename.empty() ) {
00181         outFile_ = TFile::Open(outfilename.c_str(), "recreate");
00182         
00183         outFile_->cd();
00184         //options_->GetOpt("root","outtree", doOutTree);
00185         //if(doOutTree) {
00186         // cout<<"do tree"<<endl;
00187         outEvent_ = new EventColin();
00188         outTree_ = new TTree("Eff","");
00189         outTree_->Branch("event","EventColin", &outEvent_,32000,2);
00190       }
00191       // cout<<"don't do tree"<<endl;
00192     }
00193   }
00194   // PFJet benchmark options and output jetfile to be open before input file!!!--
00195 
00196   doPFJetBenchmark_ = false;
00197   options_->GetOpt("pfjet_benchmark", "on/off", doPFJetBenchmark_);
00198   
00199   if (doPFJetBenchmark_) {
00200     string outjetfilename;
00201     options_->GetOpt("pfjet_benchmark", "outjetfile", outjetfilename);
00202         
00203     bool pfjBenchmarkDebug;
00204     options_->GetOpt("pfjet_benchmark", "debug", pfjBenchmarkDebug);
00205     
00206     bool plotAgainstReco=0;
00207     options_->GetOpt("pfjet_benchmark", "plotAgainstReco", plotAgainstReco);
00208     
00209     bool onlyTwoJets=1;
00210     options_->GetOpt("pfjet_benchmark", "onlyTwoJets", onlyTwoJets);
00211     
00212     double deltaRMax=0.1;
00213     options_->GetOpt("pfjet_benchmark", "deltaRMax", deltaRMax);
00214 
00215     fastsim_=true;
00216     options_->GetOpt("Simulation","Fast",fastsim_);
00217  
00218     PFJetBenchmark_.setup( outjetfilename, 
00219                            pfjBenchmarkDebug,
00220                            plotAgainstReco,
00221                            onlyTwoJets,
00222                            deltaRMax );
00223   }
00224 
00225 // PFMET benchmark options and output jetfile to be open before input file!!!--
00226 
00227   doPFMETBenchmark_ = false;
00228   options_->GetOpt("pfmet_benchmark", "on/off", doPFMETBenchmark_);
00229   
00230   if (doPFMETBenchmark_) {
00231     //COLIN : looks like this benchmark is not written in the standard output file. Any particular reason for it? 
00232     
00233     doMet_ = false;
00234     options_->GetOpt("MET", "on/off", doMet_);
00235 
00236     JECinCaloMet_ = false;
00237     options_->GetOpt("pfmet_benchmark", "JECinCaloMET", JECinCaloMet_);
00238 
00239     std::string outmetfilename;
00240     options_->GetOpt("pfmet_benchmark", "outmetfile", outmetfilename);
00241 
00242     // define here the various benchmark comparison
00243     metManager_.reset( new METManager(outmetfilename) );
00244     metManager_->addGenBenchmark("PF");
00245     metManager_->addGenBenchmark("Calo");
00246     if ( doMet_ ) metManager_->addGenBenchmark("recompPF");
00247     if (JECinCaloMet_) metManager_->addGenBenchmark("corrCalo");
00248 
00249     bool pfmetBenchmarkDebug;
00250     options_->GetOpt("pfmet_benchmark", "debug", pfmetBenchmarkDebug);
00251         
00252     MET1cut = 10.0;
00253     options_->GetOpt("pfmet_benchmark", "truemetcut", MET1cut);
00254     
00255     DeltaMETcut = 30.0;
00256     options_->GetOpt("pfmet_benchmark", "deltametcut", DeltaMETcut);
00257     
00258     DeltaPhicut = 0.8;
00259     options_->GetOpt("pfmet_benchmark", "deltaphicut", DeltaPhicut);
00260     
00261 
00262     std::vector<unsigned int> vIgnoreParticlesIDs;
00263     options_->GetOpt("pfmet_benchmark", "trueMetIgnoreParticlesIDs", vIgnoreParticlesIDs);
00264     //std::cout << "FL: vIgnoreParticlesIDs.size() = " << vIgnoreParticlesIDs.size() << std::endl;
00265     //std::cout << "FL: first = " << vIgnoreParticlesIDs[0] << std::endl;
00266     metManager_->SetIgnoreParticlesIDs(&vIgnoreParticlesIDs);
00267 
00268     std::vector<unsigned int> trueMetSpecificIdCut;
00269     std::vector<double> trueMetSpecificEtaCut;
00270     options_->GetOpt("pfmet_benchmark", "trueMetSpecificIdCut", trueMetSpecificIdCut);
00271     options_->GetOpt("pfmet_benchmark", "trueMetSpecificEtaCut", trueMetSpecificEtaCut);
00272     if (trueMetSpecificIdCut.size()!=trueMetSpecificEtaCut.size()) std::cout << "Warning: PFRootEventManager: trueMetSpecificIdCut.size()!=trueMetSpecificEtaCut.size()" << std::endl;
00273     else metManager_->SetSpecificIdCut(&trueMetSpecificIdCut,&trueMetSpecificEtaCut);
00274 
00275   }
00276 
00277   doPFCandidateBenchmark_ = true;
00278   options_->GetOpt("pfCandidate_benchmark", "on/off", doPFCandidateBenchmark_); 
00279   if (doPFCandidateBenchmark_) {    
00280     cout<<"+++ Setting PFCandidate benchmark"<<endl;
00281     TDirectory* dir = outFile_->mkdir("DQMData");
00282     dir = dir->mkdir("PFTask");    
00283     dir = dir->mkdir("particleFlowManager");
00284     pfCandidateManager_.setDirectory( dir );
00285 
00286     float dRMax = 0.5; 
00287     options_->GetOpt("pfCandidate_benchmark", "dRMax", dRMax); 
00288     float ptMin = 2; 
00289     options_->GetOpt("pfCandidate_benchmark", "ptMin", ptMin); 
00290     bool matchCharge = true; 
00291     options_->GetOpt("pfCandidate_benchmark", "matchCharge", matchCharge); 
00292     int mode = static_cast<int>(Benchmark::DEFAULT);
00293     options_->GetOpt("pfCandidate_benchmark", "mode", mode); 
00294     
00295     pfCandidateManager_.setParameters( dRMax, matchCharge, 
00296                                        static_cast<Benchmark::Mode>(mode));
00297     pfCandidateManager_.setRange( ptMin, 10e5, -4.5, 4.5, -10, 10);
00298     pfCandidateManager_.setup();
00299     //COLIN need to set the subdirectory.  
00300     cout<<"+++ Done "<<endl;
00301   }
00302   // Addition to have DQM histograms : by S. Dutta   
00303   doPFDQM_ = true;
00304   options_->GetOpt("pfDQM_monitor", "on/off", doPFDQM_); 
00305 
00306   if (doPFDQM_) {
00307     cout<<"+++ Setting PFDQM Monitoring " <<endl;
00308     string dqmfilename;
00309     dqmFile_ = 0;
00310     options_->GetOpt("pfDQM_monitor","DQMFilename", dqmfilename);
00311     dqmFile_ = TFile::Open(dqmfilename.c_str(), "recreate");
00312 
00313     TDirectory* dir = dqmFile_->mkdir("DQMData");
00314     TDirectory* dir1 = dir->mkdir("PFJetValidation");
00315     TDirectory* dir2 = dir->mkdir("PFMETValidation");
00316     pfJetMonitor_.setDirectory( dir1 );
00317     pfMETMonitor_.setDirectory( dir2 );
00318     float dRMax = 0.5;
00319     options_->GetOpt("pfCandidate_benchmark", "dRMax", dRMax);
00320     float ptMin = 2;
00321     options_->GetOpt("pfCandidate_benchmark", "ptMin", ptMin);
00322     bool matchCharge = true;
00323     options_->GetOpt("pfCandidate_benchmark", "matchCharge", matchCharge);
00324     int mode = static_cast<int>(Benchmark::DEFAULT);
00325     options_->GetOpt("pfCandidate_benchmark", "mode", mode);
00326 
00327     pfJetMonitor_.setParameters( dRMax, matchCharge, static_cast<Benchmark::Mode>(mode),
00328                                  ptMin, 10e5, -4.5, 4.5, -10.0, 10.0, false);
00329     pfJetMonitor_.setup();
00330 
00331     pfMETMonitor_.setParameters( static_cast<Benchmark::Mode>(mode),
00332                                  ptMin, 10e5, -4.5, 4.5, -10.0, 10.0, false);
00333     pfMETMonitor_.setup();
00334   }
00335 //-----------------------------------------------
00336 
00337 
00338   std::string outputFile0;
00339   TFile* file0 = 0;
00340   TH2F* hBNeighbour0 = 0;
00341   TH2F* hENeighbour0 = 0;
00342   options_->GetOpt("clustering", "ECAL", outputFile0);
00343   if(!outputFile0.empty() ) {
00344     file0 = TFile::Open(outputFile0.c_str(),"recreate");
00345     hBNeighbour0 = new TH2F("BNeighbour0","f_{Neighbours} vs 1/E_{Seed} (ECAL Barrel)",500, 0., 0.5, 1000,0.,1.);
00346     hENeighbour0 = new TH2F("ENeighbour0","f_{Neighbours} vs 1/E_{Seed} (ECAL Endcap)",500, 0., 0.2, 1000,0.,1.);
00347   }
00348 
00349   std::string outputFile1;
00350   TFile* file1 = 0;
00351   TH2F* hBNeighbour1 = 0;
00352   TH2F* hENeighbour1 = 0;
00353   options_->GetOpt("clustering", "HCAL", outputFile1);
00354   if(!outputFile1.empty() ) {
00355     file1 = TFile::Open(outputFile1.c_str(),"recreate");
00356     hBNeighbour1 = new TH2F("BNeighbour1","f_{Neighbours} vs 1/E_{Seed} (HCAL Barrel)",500, 0., 0.05, 400,0.,1.);
00357     hENeighbour1 = new TH2F("ENeighbour1","f_{Neighbours} vs 1/E_{Seed} (HCAL Endcap)",500, 0., 0.05, 400,0.,1.);
00358   }
00359 
00360   std::string outputFile2;
00361   TFile* file2 = 0;
00362   TH2F* hBNeighbour2 = 0;
00363   TH2F* hENeighbour2 = 0;
00364   options_->GetOpt("clustering", "HFEM", outputFile2);
00365   if(!outputFile2.empty() ) {
00366     file2 = TFile::Open(outputFile2.c_str(),"recreate");
00367     hBNeighbour2 = new TH2F("BNeighbour2","f_{Neighbours} vs 1/E_{Seed} (HFEM Barrel)",500, 0., 0.02, 400,0.,1.);
00368     hENeighbour2 = new TH2F("ENeighbour2","f_{Neighbours} vs 1/E_{Seed} (HFEM Endcap)",500, 0., 0.02, 400,0.,1.);
00369   }
00370 
00371   std::string outputFile3;
00372   TFile* file3 = 0;
00373   TH2F* hBNeighbour3 = 0;
00374   TH2F* hENeighbour3 = 0;
00375   options_->GetOpt("clustering", "HFHAD", outputFile3);
00376   if(!outputFile3.empty() ) {
00377     file3 = TFile::Open(outputFile3.c_str(),"recreate");
00378     hBNeighbour3 = new TH2F("BNeighbour3","f_{Neighbours} vs 1/E_{Seed} (HFHAD Barrel)",500, 0., 0.02, 400,0.,1.);
00379     hENeighbour3 = new TH2F("ENeighbour3","f_{Neighbours} vs 1/E_{Seed} (HFHAD Endcap)",500, 0., 0.02, 400,0.,1.);
00380   }
00381 
00382   std::string outputFile4;
00383   TFile* file4 = 0;
00384   TH2F* hBNeighbour4 = 0;
00385   TH2F* hENeighbour4 = 0;
00386   options_->GetOpt("clustering", "Preshower", outputFile4);
00387   if(!outputFile4.empty() ) {
00388     file4 = TFile::Open(outputFile4.c_str(),"recreate");
00389     hBNeighbour4 = new TH2F("BNeighbour4","f_{Neighbours} vs 1/E_{Seed} (Preshower Barrel)",200, 0., 1000., 400,0.,1.);
00390     hENeighbour4 = new TH2F("ENeighbour4","f_{Neighbours} vs 1/E_{Seed} (Preshower Endcap)",200, 0., 1000., 400,0.,1.);
00391   }
00392 
00393   // input root file --------------------------------------------
00394 
00395   if( reconnect )
00396     connect(); 
00397 
00398   // filter --------------------------------------------------------------
00399 
00400   filterNParticles_ = 0;
00401   options_->GetOpt("filter", "nparticles", filterNParticles_);
00402   
00403   filterHadronicTaus_ = true;
00404   options_->GetOpt("filter", "hadronic_taus", filterHadronicTaus_);
00405   
00406   filterTaus_.clear();
00407   options_->GetOpt("filter", "taus", filterTaus_);
00408   if( !filterTaus_.empty() &&
00409       filterTaus_.size() != 2 ) {
00410     cerr<<"PFRootEventManager::ReadOptions, bad filter/taus option."<<endl
00411         <<"please use : "<<endl
00412         <<"\tfilter taus n_charged n_neutral"<<endl;
00413     filterTaus_.clear();
00414   }
00415   
00416   
00417   // clustering parameters -----------------------------------------------
00418 
00419   doClustering_ = true;
00420   //options_->GetOpt("clustering", "on/off", doClustering_);
00421   
00422   bool clusteringDebug = false;
00423   options_->GetOpt("clustering", "debug", clusteringDebug );
00424 
00425   findRecHitNeighbours_ = true;
00426   options_->GetOpt("clustering", "findRecHitNeighbours", 
00427                    findRecHitNeighbours_);
00428   
00429   double threshEcalBarrel = 0.1;
00430   options_->GetOpt("clustering", "thresh_Ecal_Barrel", threshEcalBarrel);
00431   
00432   double threshPtEcalBarrel = 0.0;
00433   options_->GetOpt("clustering", "thresh_Pt_Ecal_Barrel", threshPtEcalBarrel);
00434   
00435   double threshSeedEcalBarrel = 0.3;
00436   options_->GetOpt("clustering", "thresh_Seed_Ecal_Barrel", 
00437                    threshSeedEcalBarrel);
00438 
00439   double threshPtSeedEcalBarrel = 0.0;
00440   options_->GetOpt("clustering", "thresh_Pt_Seed_Ecal_Barrel", 
00441                    threshPtSeedEcalBarrel);
00442 
00443   double threshCleanEcalBarrel = 1E5;
00444   options_->GetOpt("clustering", "thresh_Clean_Ecal_Barrel", 
00445                    threshCleanEcalBarrel);
00446 
00447   std::vector<double> minS4S1CleanEcalBarrel;
00448   options_->GetOpt("clustering", "minS4S1_Clean_Ecal_Barrel", 
00449                    minS4S1CleanEcalBarrel);
00450 
00451   double threshDoubleSpikeEcalBarrel = 10.;
00452   options_->GetOpt("clustering", "thresh_DoubleSpike_Ecal_Barrel", 
00453                    threshDoubleSpikeEcalBarrel);
00454 
00455   double minS6S2DoubleSpikeEcalBarrel = 0.04;
00456   options_->GetOpt("clustering", "minS6S2_DoubleSpike_Ecal_Barrel", 
00457                    minS6S2DoubleSpikeEcalBarrel);
00458 
00459   double threshEcalEndcap = 0.2;
00460   options_->GetOpt("clustering", "thresh_Ecal_Endcap", threshEcalEndcap);
00461 
00462   double threshPtEcalEndcap = 0.0;
00463   options_->GetOpt("clustering", "thresh_Pt_Ecal_Endcap", threshPtEcalEndcap);
00464 
00465   double threshSeedEcalEndcap = 0.8;
00466   options_->GetOpt("clustering", "thresh_Seed_Ecal_Endcap",
00467                    threshSeedEcalEndcap);
00468 
00469   double threshPtSeedEcalEndcap = 0.0;
00470   options_->GetOpt("clustering", "thresh_Pt_Seed_Ecal_Endcap",
00471                    threshPtSeedEcalEndcap);
00472 
00473   double threshCleanEcalEndcap = 1E5;
00474   options_->GetOpt("clustering", "thresh_Clean_Ecal_Endcap", 
00475                    threshCleanEcalEndcap);
00476 
00477   std::vector<double> minS4S1CleanEcalEndcap;
00478   options_->GetOpt("clustering", "minS4S1_Clean_Ecal_Endcap", 
00479                    minS4S1CleanEcalEndcap);
00480 
00481   double threshDoubleSpikeEcalEndcap = 1E9;
00482   options_->GetOpt("clustering", "thresh_DoubleSpike_Ecal_Endcap", 
00483                    threshDoubleSpikeEcalEndcap);
00484 
00485   double minS6S2DoubleSpikeEcalEndcap = -1.;
00486   options_->GetOpt("clustering", "minS6S2_DoubleSpike_Ecal_Endcap", 
00487                    minS6S2DoubleSpikeEcalEndcap);
00488 
00489   double showerSigmaEcal = 3;  
00490   options_->GetOpt("clustering", "shower_Sigma_Ecal",
00491                    showerSigmaEcal);
00492 
00493   int nNeighboursEcal = 4;
00494   options_->GetOpt("clustering", "neighbours_Ecal", nNeighboursEcal);
00495   
00496   int posCalcNCrystalEcal = -1;
00497   options_->GetOpt("clustering", "posCalc_nCrystal_Ecal", 
00498                    posCalcNCrystalEcal);
00499 
00500   double posCalcP1Ecal 
00501     = threshEcalBarrel<threshEcalEndcap ? threshEcalBarrel:threshEcalEndcap;
00502 //   options_->GetOpt("clustering", "posCalc_p1_Ecal", 
00503 //                    posCalcP1Ecal);
00504   
00505   bool useCornerCellsEcal = false;
00506   options_->GetOpt("clustering", "useCornerCells_Ecal",
00507                    useCornerCellsEcal);
00508 
00509   clusterAlgoECAL_.setHistos(file0,hBNeighbour0,hENeighbour0);
00510 
00511   clusterAlgoECAL_.setThreshBarrel( threshEcalBarrel );
00512   clusterAlgoECAL_.setThreshSeedBarrel( threshSeedEcalBarrel );
00513   
00514   clusterAlgoECAL_.setThreshPtBarrel( threshPtEcalBarrel );
00515   clusterAlgoECAL_.setThreshPtSeedBarrel( threshPtSeedEcalBarrel );
00516   
00517   clusterAlgoECAL_.setThreshCleanBarrel(threshCleanEcalBarrel);
00518   clusterAlgoECAL_.setS4S1CleanBarrel(minS4S1CleanEcalBarrel);
00519 
00520   clusterAlgoECAL_.setThreshDoubleSpikeBarrel( threshDoubleSpikeEcalBarrel );
00521   clusterAlgoECAL_.setS6S2DoubleSpikeBarrel( minS6S2DoubleSpikeEcalBarrel );
00522 
00523   clusterAlgoECAL_.setThreshEndcap( threshEcalEndcap );
00524   clusterAlgoECAL_.setThreshSeedEndcap( threshSeedEcalEndcap );
00525 
00526   clusterAlgoECAL_.setThreshPtEndcap( threshPtEcalEndcap );
00527   clusterAlgoECAL_.setThreshPtSeedEndcap( threshPtSeedEcalEndcap );
00528 
00529   clusterAlgoECAL_.setThreshCleanEndcap(threshCleanEcalEndcap);
00530   clusterAlgoECAL_.setS4S1CleanEndcap(minS4S1CleanEcalEndcap);
00531 
00532   clusterAlgoECAL_.setThreshDoubleSpikeEndcap( threshDoubleSpikeEcalEndcap );
00533   clusterAlgoECAL_.setS6S2DoubleSpikeEndcap( minS6S2DoubleSpikeEcalEndcap );
00534 
00535   clusterAlgoECAL_.setNNeighbours( nNeighboursEcal );
00536   clusterAlgoECAL_.setShowerSigma( showerSigmaEcal );
00537 
00538   clusterAlgoECAL_.setPosCalcNCrystal( posCalcNCrystalEcal );
00539   clusterAlgoECAL_.setPosCalcP1( posCalcP1Ecal );
00540 
00541   clusterAlgoECAL_.setUseCornerCells( useCornerCellsEcal );
00542 
00543   clusterAlgoECAL_.enableDebugging( clusteringDebug ); 
00544 
00545   int dcormode = 0;
00546   options_->GetOpt("clustering", "depthCor_Mode", dcormode);
00547   
00548   double dcora = -1;
00549   options_->GetOpt("clustering", "depthCor_A", dcora);
00550   double dcorb = -1;
00551   options_->GetOpt("clustering", "depthCor_B", dcorb);
00552   double dcorap = -1;
00553   options_->GetOpt("clustering", "depthCor_A_preshower", dcorap);
00554   double dcorbp = -1;
00555   options_->GetOpt("clustering", "depthCor_B_preshower", dcorbp);
00556 
00557   //   if( dcormode > 0 && 
00558   //       dcora > -0.5 && 
00559   //       dcorb > -0.5 && 
00560   //       dcorap > -0.5 && 
00561   //       dcorbp > -0.5 ) {
00562 
00563   //     cout<<"set depth correction "
00564   //    <<dcormode<<" "<<dcora<<" "<<dcorb<<" "<<dcorap<<" "<<dcorbp<<endl;
00565   reco::PFCluster::setDepthCorParameters( dcormode, 
00566                                           dcora, dcorb, 
00567                                           dcorap, dcorbp);
00568   //   }
00569   //   else {
00570   //     reco::PFCluster::setDepthCorParameters( -1, 
00571   //                                        0,0 , 
00572   //                                        0,0 );
00573   //   }
00574 
00575   
00576 
00577   double threshHcalBarrel = 0.8;
00578   options_->GetOpt("clustering", "thresh_Hcal_Barrel", threshHcalBarrel);
00579   
00580   double threshPtHcalBarrel = 0.0;
00581   options_->GetOpt("clustering", "thresh_Pt_Hcal_Barrel", threshPtHcalBarrel);
00582   
00583   double threshSeedHcalBarrel = 1.4;
00584   options_->GetOpt("clustering", "thresh_Seed_Hcal_Barrel", 
00585                    threshSeedHcalBarrel);
00586 
00587   double threshPtSeedHcalBarrel = 0.0;
00588   options_->GetOpt("clustering", "thresh_Pt_Seed_Hcal_Barrel", 
00589                    threshPtSeedHcalBarrel);
00590 
00591   double threshCleanHcalBarrel = 1E5;
00592   options_->GetOpt("clustering", "thresh_Clean_Hcal_Barrel", 
00593                    threshCleanHcalBarrel);
00594 
00595   std::vector<double> minS4S1CleanHcalBarrel;
00596   options_->GetOpt("clustering", "minS4S1_Clean_Hcal_Barrel", 
00597                    minS4S1CleanHcalBarrel);
00598 
00599   double threshHcalEndcap = 0.8;
00600   options_->GetOpt("clustering", "thresh_Hcal_Endcap", threshHcalEndcap);
00601 
00602   double threshPtHcalEndcap = 0.0;
00603   options_->GetOpt("clustering", "thresh_Pt_Hcal_Endcap", threshPtHcalEndcap);
00604 
00605   double threshSeedHcalEndcap = 1.4;
00606   options_->GetOpt("clustering", "thresh_Seed_Hcal_Endcap",
00607                    threshSeedHcalEndcap);
00608 
00609   double threshPtSeedHcalEndcap = 0.0;
00610   options_->GetOpt("clustering", "thresh_Pt_Seed_Hcal_Endcap",
00611                    threshPtSeedHcalEndcap);
00612 
00613   double threshCleanHcalEndcap = 1E5;
00614   options_->GetOpt("clustering", "thresh_Clean_Hcal_Endcap", 
00615                    threshCleanHcalEndcap);
00616 
00617   std::vector<double> minS4S1CleanHcalEndcap;
00618   options_->GetOpt("clustering", "minS4S1_Clean_Hcal_Endcap", 
00619                    minS4S1CleanHcalEndcap);
00620 
00621   double showerSigmaHcal    = 15;
00622   options_->GetOpt("clustering", "shower_Sigma_Hcal",
00623                    showerSigmaHcal);
00624  
00625   int nNeighboursHcal = 4;
00626   options_->GetOpt("clustering", "neighbours_Hcal", nNeighboursHcal);
00627 
00628   int posCalcNCrystalHcal = 5;
00629   options_->GetOpt("clustering", "posCalc_nCrystal_Hcal",
00630                    posCalcNCrystalHcal);
00631 
00632   bool useCornerCellsHcal = false;
00633   options_->GetOpt("clustering", "useCornerCells_Hcal",
00634                    useCornerCellsHcal);
00635 
00636   bool cleanRBXandHPDs = false;
00637   options_->GetOpt("clustering", "cleanRBXandHPDs_Hcal",
00638                    cleanRBXandHPDs);
00639 
00640   double posCalcP1Hcal 
00641     = threshHcalBarrel<threshHcalEndcap ? threshHcalBarrel:threshHcalEndcap;
00642 //   options_->GetOpt("clustering", "posCalc_p1_Hcal", 
00643 //                    posCalcP1Hcal);
00644 
00645 
00646   clusterAlgoHCAL_.setHistos(file1,hBNeighbour1,hENeighbour1);
00647 
00648 
00649   clusterAlgoHCAL_.setThreshBarrel( threshHcalBarrel );
00650   clusterAlgoHCAL_.setThreshSeedBarrel( threshSeedHcalBarrel );
00651   
00652   clusterAlgoHCAL_.setThreshPtBarrel( threshPtHcalBarrel );
00653   clusterAlgoHCAL_.setThreshPtSeedBarrel( threshPtSeedHcalBarrel );
00654   
00655   clusterAlgoHCAL_.setThreshCleanBarrel(threshCleanHcalBarrel);
00656   clusterAlgoHCAL_.setS4S1CleanBarrel(minS4S1CleanHcalBarrel);
00657 
00658   clusterAlgoHCAL_.setThreshEndcap( threshHcalEndcap );
00659   clusterAlgoHCAL_.setThreshSeedEndcap( threshSeedHcalEndcap );
00660 
00661   clusterAlgoHCAL_.setThreshPtEndcap( threshPtHcalEndcap );
00662   clusterAlgoHCAL_.setThreshPtSeedEndcap( threshPtSeedHcalEndcap );
00663 
00664   clusterAlgoHCAL_.setThreshCleanEndcap(threshCleanHcalEndcap);
00665   clusterAlgoHCAL_.setS4S1CleanEndcap(minS4S1CleanHcalEndcap);
00666 
00667   clusterAlgoHCAL_.setNNeighbours( nNeighboursHcal );
00668   clusterAlgoHCAL_.setShowerSigma( showerSigmaHcal );
00669 
00670   clusterAlgoHCAL_.setPosCalcNCrystal( posCalcNCrystalHcal );
00671   clusterAlgoHCAL_.setPosCalcP1( posCalcP1Hcal );
00672 
00673   clusterAlgoHCAL_.setUseCornerCells( useCornerCellsHcal );
00674   clusterAlgoHCAL_.setCleanRBXandHPDs( cleanRBXandHPDs );
00675 
00676   clusterAlgoHCAL_.enableDebugging( clusteringDebug ); 
00677 
00678 
00679   // clustering HF EM 
00680 
00681   double threshHFEM = 0.;
00682   options_->GetOpt("clustering", "thresh_HFEM", threshHFEM);
00683   
00684   double threshPtHFEM = 0.;
00685   options_->GetOpt("clustering", "thresh_Pt_HFEM", threshPtHFEM);
00686   
00687   double threshSeedHFEM = 0.001;
00688   options_->GetOpt("clustering", "thresh_Seed_HFEM", 
00689                    threshSeedHFEM);
00690   
00691   double threshPtSeedHFEM = 0.0;
00692   options_->GetOpt("clustering", "thresh_Pt_Seed_HFEM", 
00693                    threshPtSeedHFEM);
00694   
00695   double threshCleanHFEM = 1E5;
00696   options_->GetOpt("clustering", "thresh_Clean_HFEM", 
00697                    threshCleanHFEM);
00698 
00699   std::vector<double> minS4S1CleanHFEM;
00700   options_->GetOpt("clustering", "minS4S1_Clean_HFEM", 
00701                    minS4S1CleanHFEM);
00702 
00703   double showerSigmaHFEM    = 0.1;
00704   options_->GetOpt("clustering", "shower_Sigma_HFEM",
00705                    showerSigmaHFEM);
00706  
00707   int nNeighboursHFEM = 4;
00708   options_->GetOpt("clustering", "neighbours_HFEM", nNeighboursHFEM);
00709 
00710   int posCalcNCrystalHFEM = -1;
00711   options_->GetOpt("clustering", "posCalc_nCrystal_HFEM",
00712                    posCalcNCrystalHFEM);
00713 
00714   bool useCornerCellsHFEM = false;
00715   options_->GetOpt("clustering", "useCornerCells_HFEM",
00716                    useCornerCellsHFEM);
00717 
00718   double posCalcP1HFEM = threshHFEM;
00719 //   options_->GetOpt("clustering", "posCalc_p1_HFEM", 
00720 //                    posCalcP1HFEM);
00721 
00722 
00723   clusterAlgoHFEM_.setHistos(file2,hBNeighbour2,hENeighbour2);
00724 
00725   clusterAlgoHFEM_.setThreshEndcap( threshHFEM );
00726   clusterAlgoHFEM_.setThreshSeedEndcap( threshSeedHFEM );
00727 
00728   clusterAlgoHFEM_.setThreshPtEndcap( threshPtHFEM );
00729   clusterAlgoHFEM_.setThreshPtSeedEndcap( threshPtSeedHFEM );
00730 
00731   clusterAlgoHFEM_.setThreshCleanEndcap(threshCleanHFEM);
00732   clusterAlgoHFEM_.setS4S1CleanEndcap(minS4S1CleanHFEM);
00733 
00734   clusterAlgoHFEM_.setNNeighbours( nNeighboursHFEM );
00735   clusterAlgoHFEM_.setShowerSigma( showerSigmaHFEM );
00736 
00737   clusterAlgoHFEM_.setPosCalcNCrystal( posCalcNCrystalHFEM );
00738   clusterAlgoHFEM_.setPosCalcP1( posCalcP1HFEM );
00739 
00740   clusterAlgoHFEM_.setUseCornerCells( useCornerCellsHFEM );
00741 
00742   clusterAlgoHFEM_.enableDebugging( clusteringDebug ); 
00743 
00744   // clustering HFHAD 
00745 
00746   double threshHFHAD = 0.;
00747   options_->GetOpt("clustering", "thresh_HFHAD", threshHFHAD);
00748   
00749   double threshPtHFHAD = 0.;
00750   options_->GetOpt("clustering", "thresh_Pt_HFHAD", threshPtHFHAD);
00751   
00752   double threshSeedHFHAD = 0.001;
00753   options_->GetOpt("clustering", "thresh_Seed_HFHAD", 
00754                    threshSeedHFHAD);
00755   
00756   double threshPtSeedHFHAD = 0.0;
00757   options_->GetOpt("clustering", "thresh_Pt_Seed_HFHAD", 
00758                    threshPtSeedHFHAD);
00759   
00760   double threshCleanHFHAD = 1E5;
00761   options_->GetOpt("clustering", "thresh_Clean_HFHAD", 
00762                    threshCleanHFHAD);
00763 
00764   std::vector<double> minS4S1CleanHFHAD;
00765   options_->GetOpt("clustering", "minS4S1_Clean_HFHAD", 
00766                    minS4S1CleanHFHAD);
00767 
00768   double showerSigmaHFHAD    = 0.1;
00769   options_->GetOpt("clustering", "shower_Sigma_HFHAD",
00770                    showerSigmaHFHAD);
00771  
00772   int nNeighboursHFHAD = 4;
00773   options_->GetOpt("clustering", "neighbours_HFHAD", nNeighboursHFHAD);
00774 
00775   int posCalcNCrystalHFHAD = -1;
00776   options_->GetOpt("clustering", "posCalc_nCrystal_HFHAD",
00777                    posCalcNCrystalHFHAD);
00778 
00779   bool useCornerCellsHFHAD = false;
00780   options_->GetOpt("clustering", "useCornerCells_HFHAD",
00781                    useCornerCellsHFHAD);
00782 
00783   double posCalcP1HFHAD = threshHFHAD;
00784 //   options_->GetOpt("clustering", "posCalc_p1_HFHAD", 
00785 //                    posCalcP1HFHAD);
00786 
00787 
00788   clusterAlgoHFHAD_.setHistos(file3,hBNeighbour3,hENeighbour3);
00789 
00790   clusterAlgoHFHAD_.setThreshEndcap( threshHFHAD );
00791   clusterAlgoHFHAD_.setThreshSeedEndcap( threshSeedHFHAD );
00792 
00793   clusterAlgoHFHAD_.setThreshPtEndcap( threshPtHFHAD );
00794   clusterAlgoHFHAD_.setThreshPtSeedEndcap( threshPtSeedHFHAD );
00795 
00796   clusterAlgoHFHAD_.setThreshCleanEndcap(threshCleanHFHAD);
00797   clusterAlgoHFHAD_.setS4S1CleanEndcap(minS4S1CleanHFHAD);
00798 
00799   clusterAlgoHFHAD_.setNNeighbours( nNeighboursHFHAD );
00800   clusterAlgoHFHAD_.setShowerSigma( showerSigmaHFHAD );
00801 
00802   clusterAlgoHFHAD_.setPosCalcNCrystal( posCalcNCrystalHFHAD );
00803   clusterAlgoHFHAD_.setPosCalcP1( posCalcP1HFHAD );
00804 
00805   clusterAlgoHFHAD_.setUseCornerCells( useCornerCellsHFHAD );
00806 
00807   clusterAlgoHFHAD_.enableDebugging( clusteringDebug ); 
00808 
00809   // clustering preshower
00810 
00811   double threshPS = 0.0001;
00812   options_->GetOpt("clustering", "thresh_PS", threshPS);
00813   
00814   double threshPtPS = 0.0;
00815   options_->GetOpt("clustering", "thresh_Pt_PS", threshPtPS);
00816   
00817   double threshSeedPS = 0.001;
00818   options_->GetOpt("clustering", "thresh_Seed_PS", 
00819                    threshSeedPS);
00820   
00821   double threshPtSeedPS = 0.0;
00822   options_->GetOpt("clustering", "thresh_Pt_Seed_PS", threshPtSeedPS);
00823   
00824   double threshCleanPS = 1E5;
00825   options_->GetOpt("clustering", "thresh_Clean_PS", threshCleanPS);
00826 
00827   std::vector<double> minS4S1CleanPS;
00828   options_->GetOpt("clustering", "minS4S1_Clean_PS", minS4S1CleanPS);
00829 
00830   //Comment Michel: PSBarrel shall be removed?
00831   double threshPSBarrel     = threshPS;
00832   double threshSeedPSBarrel = threshSeedPS;
00833 
00834   double threshPtPSBarrel     = threshPtPS;
00835   double threshPtSeedPSBarrel = threshPtSeedPS;
00836 
00837   double threshCleanPSBarrel = threshCleanPS;
00838   std::vector<double> minS4S1CleanPSBarrel = minS4S1CleanPS;
00839 
00840   double threshPSEndcap     = threshPS;
00841   double threshSeedPSEndcap = threshSeedPS;
00842 
00843   double threshPtPSEndcap     = threshPtPS;
00844   double threshPtSeedPSEndcap = threshPtSeedPS;
00845 
00846   double threshCleanPSEndcap = threshCleanPS;
00847   std::vector<double> minS4S1CleanPSEndcap = minS4S1CleanPS;
00848 
00849   double showerSigmaPS    = 0.1;
00850   options_->GetOpt("clustering", "shower_Sigma_PS",
00851                    showerSigmaPS);
00852  
00853   int nNeighboursPS = 4;
00854   options_->GetOpt("clustering", "neighbours_PS", nNeighboursPS);
00855 
00856   int posCalcNCrystalPS = -1;
00857   options_->GetOpt("clustering", "posCalc_nCrystal_PS",
00858                    posCalcNCrystalPS);
00859 
00860   bool useCornerCellsPS = false;
00861   options_->GetOpt("clustering", "useCornerCells_PS",
00862                    useCornerCellsPS);
00863 
00864   double posCalcP1PS = threshPS;
00865 //   options_->GetOpt("clustering", "posCalc_p1_PS", 
00866 //                    posCalcP1PS);
00867 
00868 
00869   clusterAlgoPS_.setHistos(file4,hBNeighbour4,hENeighbour4);
00870 
00871 
00872 
00873   clusterAlgoPS_.setThreshBarrel( threshPSBarrel );
00874   clusterAlgoPS_.setThreshSeedBarrel( threshSeedPSBarrel );
00875   
00876   clusterAlgoPS_.setThreshPtBarrel( threshPtPSBarrel );
00877   clusterAlgoPS_.setThreshPtSeedBarrel( threshPtSeedPSBarrel );
00878   
00879   clusterAlgoPS_.setThreshCleanBarrel(threshCleanPSBarrel);
00880   clusterAlgoPS_.setS4S1CleanBarrel(minS4S1CleanPSBarrel);
00881 
00882   clusterAlgoPS_.setThreshEndcap( threshPSEndcap );
00883   clusterAlgoPS_.setThreshSeedEndcap( threshSeedPSEndcap );
00884 
00885   clusterAlgoPS_.setThreshPtEndcap( threshPtPSEndcap );
00886   clusterAlgoPS_.setThreshPtSeedEndcap( threshPtSeedPSEndcap );
00887 
00888   clusterAlgoPS_.setThreshCleanEndcap(threshCleanPSEndcap);
00889   clusterAlgoPS_.setS4S1CleanEndcap(minS4S1CleanPSEndcap);
00890 
00891   clusterAlgoPS_.setNNeighbours( nNeighboursPS );
00892   clusterAlgoPS_.setShowerSigma( showerSigmaPS );
00893 
00894   clusterAlgoPS_.setPosCalcNCrystal( posCalcNCrystalPS );
00895   clusterAlgoPS_.setPosCalcP1( posCalcP1PS );
00896 
00897   clusterAlgoPS_.setUseCornerCells( useCornerCellsPS );
00898 
00899   clusterAlgoPS_.enableDebugging( clusteringDebug ); 
00900 
00901   // options for particle flow ---------------------------------------------
00902 
00903 
00904   doParticleFlow_ = true;
00905   doCompare_ = false;
00906   options_->GetOpt("particle_flow", "on/off", doParticleFlow_);  
00907   options_->GetOpt("particle_flow", "comparison", doCompare_);  
00908 
00909   std::vector<double> DPtovPtCut;
00910   std::vector<unsigned> NHitCut;
00911   bool useIterTracking;
00912   int nuclearInteractionsPurity;
00913   options_->GetOpt("particle_flow", "DPtoverPt_Cut", DPtovPtCut);
00914   options_->GetOpt("particle_flow", "NHit_Cut", NHitCut);
00915   options_->GetOpt("particle_flow", "useIterTracking", useIterTracking);
00916   options_->GetOpt("particle_flow", "nuclearInteractionsPurity", nuclearInteractionsPurity);
00917   
00918   
00919   std::vector<double> PhotonSelectionCuts;
00920   options_->GetOpt("particle_flow","useEGPhotons",useEGPhotons_);
00921   options_->GetOpt("particle_flow","photonSelection", PhotonSelectionCuts);
00922 
00923   try {
00924     pfBlockAlgo_.setParameters( DPtovPtCut, 
00925                                 NHitCut,
00926                                 useIterTracking,
00927                                 useConvBremPFRecTracks_,
00928                                 nuclearInteractionsPurity,
00929                                 useEGPhotons_,
00930                                 PhotonSelectionCuts); 
00931   }  
00932   catch( std::exception& err ) {
00933     cerr<<"exception setting PFBlockAlgo parameters: "
00934         <<err.what()<<". terminating."<<endl;
00935     delete this;
00936     exit(1);
00937   }
00938   
00939 
00940   bool blockAlgoDebug = false;
00941   options_->GetOpt("blockAlgo", "debug",  blockAlgoDebug);  
00942   pfBlockAlgo_.setDebug( blockAlgoDebug );
00943 
00944   bool AlgoDebug = false;
00945   options_->GetOpt("PFAlgo", "debug",  AlgoDebug);  
00946   pfAlgo_.setDebug( AlgoDebug );
00947 
00948   // read PFCluster calibration parameters
00949   boost::shared_ptr<PFEnergyCalibration> 
00950     calibration( new PFEnergyCalibration() );
00951   calibration_ = calibration;
00952 
00953   bool usePFSCEleCalib;
00954   std::vector<double>  calibPFSCEle_Fbrem_barrel; 
00955   std::vector<double>  calibPFSCEle_Fbrem_endcap;
00956   std::vector<double>  calibPFSCEle_barrel;
00957   std::vector<double>  calibPFSCEle_endcap;
00958   options_->GetOpt("particle_flow","usePFSCEleCalib",usePFSCEleCalib);
00959   options_->GetOpt("particle_flow","calibPFSCEle_Fbrem_barrel",calibPFSCEle_Fbrem_barrel);
00960   options_->GetOpt("particle_flow","calibPFSCEle_Fbrem_endcap",calibPFSCEle_Fbrem_endcap);
00961   options_->GetOpt("particle_flow","calibPFSCEle_barrel",calibPFSCEle_barrel);
00962   options_->GetOpt("particle_flow","calibPFSCEle_endcap",calibPFSCEle_endcap);
00963   boost::shared_ptr<PFSCEnergyCalibration>  
00964     thePFSCEnergyCalibration (new PFSCEnergyCalibration(calibPFSCEle_Fbrem_barrel,calibPFSCEle_Fbrem_endcap,
00965                                                         calibPFSCEle_barrel,calibPFSCEle_endcap ));
00966   
00967   bool useEGammaSupercluster;
00968   double sumEtEcalIsoForEgammaSC_barrel;
00969   double sumEtEcalIsoForEgammaSC_endcap;
00970   double coneEcalIsoForEgammaSC;
00971   double sumPtTrackIsoForEgammaSC_barrel;
00972   double sumPtTrackIsoForEgammaSC_endcap;
00973   unsigned int nTrackIsoForEgammaSC;
00974   double coneTrackIsoForEgammaSC;
00975   bool useEGammaElectrons;
00976   options_->GetOpt("particle_flow","useEGammaSupercluster",useEGammaSupercluster);
00977   options_->GetOpt("particle_flow","sumEtEcalIsoForEgammaSC_barrel",sumEtEcalIsoForEgammaSC_barrel);
00978   options_->GetOpt("particle_flow","sumEtEcalIsoForEgammaSC_endcap",sumEtEcalIsoForEgammaSC_endcap);
00979   options_->GetOpt("particle_flow","coneEcalIsoForEgammaSC",coneEcalIsoForEgammaSC);
00980   options_->GetOpt("particle_flow","sumPtTrackIsoForEgammaSC_barrel",sumPtTrackIsoForEgammaSC_barrel);
00981   options_->GetOpt("particle_flow","sumPtTrackIsoForEgammaSC_endcap",sumPtTrackIsoForEgammaSC_endcap);
00982   options_->GetOpt("particle_flow","nTrackIsoForEgammaSC",nTrackIsoForEgammaSC);
00983   options_->GetOpt("particle_flow","coneTrackIsoForEgammaSC",coneTrackIsoForEgammaSC);
00984   options_->GetOpt("particle_flow","useEGammaElectrons",useEGammaElectrons);
00985 
00986   //--ab: get calibration factors for HF:
00987   bool calibHF_use = false;
00988   std::vector<double>  calibHF_eta_step;
00989   std::vector<double>  calibHF_a_EMonly;
00990   std::vector<double>  calibHF_b_HADonly;
00991   std::vector<double>  calibHF_a_EMHAD;
00992   std::vector<double>  calibHF_b_EMHAD;
00993 
00994   options_->GetOpt("particle_flow","calib_calibHF_use",calibHF_use);
00995   options_->GetOpt("particle_flow","calib_calibHF_eta_step",calibHF_eta_step);
00996   options_->GetOpt("particle_flow","calib_calibHF_a_EMonly",calibHF_a_EMonly);
00997   options_->GetOpt("particle_flow","calib_calibHF_b_HADonly",calibHF_b_HADonly);
00998   options_->GetOpt("particle_flow","calib_calibHF_a_EMHAD",calibHF_a_EMHAD);
00999   options_->GetOpt("particle_flow","calib_calibHF_b_EMHAD",calibHF_b_EMHAD);
01000 
01001   boost::shared_ptr<PFEnergyCalibrationHF>  thepfEnergyCalibrationHF
01002     ( new PFEnergyCalibrationHF(calibHF_use,calibHF_eta_step,calibHF_a_EMonly,calibHF_b_HADonly,calibHF_a_EMHAD,calibHF_b_EMHAD) ) ;
01003 
01004   thepfEnergyCalibrationHF_ = thepfEnergyCalibrationHF;
01005 
01006 
01007   //----------------------------------------
01008   double nSigmaECAL = 99999;
01009   options_->GetOpt("particle_flow", "nsigma_ECAL", nSigmaECAL);
01010   double nSigmaHCAL = 99999;
01011   options_->GetOpt("particle_flow", "nsigma_HCAL", nSigmaHCAL);
01012 
01013   try {
01014     pfAlgo_.setParameters( nSigmaECAL, nSigmaHCAL, 
01015                            calibration, thepfEnergyCalibrationHF_);
01016   }
01017   catch( std::exception& err ) {
01018     cerr<<"exception setting PFAlgo parameters: "
01019         <<err.what()<<". terminating."<<endl;
01020     delete this;
01021     exit(1);
01022   }
01023 
01024   std::vector<double> muonHCAL;
01025   std::vector<double> muonECAL;
01026   options_->GetOpt("particle_flow", "muon_HCAL", muonHCAL);
01027   options_->GetOpt("particle_flow", "muon_ECAL", muonECAL);
01028   assert ( muonHCAL.size() == 2 && muonECAL.size() == 2 );
01029 
01030   double nSigmaTRACK = 3.0;
01031   options_->GetOpt("particle_flow", "nsigma_TRACK", nSigmaTRACK);
01032 
01033   double ptError = 1.0;
01034   options_->GetOpt("particle_flow", "pt_error", ptError);
01035   
01036   std::vector<double> factors45;
01037   options_->GetOpt("particle_flow", "factors_45", factors45);
01038   assert ( factors45.size() == 2 );
01039   
01040   bool usePFMuonMomAssign = false;
01041   options_->GetOpt("particle_flow", "usePFMuonMomAssign", usePFMuonMomAssign);
01042  
01043   try { 
01044     pfAlgo_.setPFMuonAndFakeParameters(muonHCAL,
01045                                        muonECAL,
01046                                        nSigmaTRACK,
01047                                        ptError,
01048                                        factors45,
01049                                        usePFMuonMomAssign);
01050   }
01051   catch( std::exception& err ) {
01052     cerr<<"exception setting PFAlgo Muon and Fake parameters: "
01053         <<err.what()<<". terminating."<<endl;
01054     delete this;
01055     exit(1);
01056   }
01057   
01058   bool postHFCleaning = true;
01059   options_->GetOpt("particle_flow", "postHFCleaning", postHFCleaning);
01060   double minHFCleaningPt = 5.;
01061   options_->GetOpt("particle_flow", "minHFCleaningPt", minHFCleaningPt);
01062   double minSignificance = 2.5;
01063   options_->GetOpt("particle_flow", "minSignificance", minSignificance);
01064   double maxSignificance = 2.5;
01065   options_->GetOpt("particle_flow", "maxSignificance", maxSignificance);
01066   double minSignificanceReduction = 1.4;
01067   options_->GetOpt("particle_flow", "minSignificanceReduction", minSignificanceReduction);
01068   double maxDeltaPhiPt = 7.0;
01069   options_->GetOpt("particle_flow", "maxDeltaPhiPt", maxDeltaPhiPt);
01070   double minDeltaMet = 0.4;
01071   options_->GetOpt("particle_flow", "minDeltaMet", minDeltaMet);
01072 
01073   // Set post HF cleaning muon parameters
01074   try { 
01075     pfAlgo_.setPostHFCleaningParameters(postHFCleaning,
01076                                         minHFCleaningPt,
01077                                         minSignificance,
01078                                         maxSignificance,
01079                                         minSignificanceReduction,
01080                                         maxDeltaPhiPt,
01081                                         minDeltaMet);
01082   }
01083   catch( std::exception& err ) {
01084     cerr<<"exception setting post HF cleaning parameters: "
01085         <<err.what()<<". terminating."<<endl;
01086     delete this;
01087     exit(1);
01088   }
01089   
01090   useAtHLT = false;
01091   options_->GetOpt("particle_flow", "useAtHLT", useAtHLT);
01092   cout<<"use HLT tracking "<<useAtHLT<<endl;
01093 
01094 
01095   bool usePFElectrons = false;   // set true to use PFElectrons
01096   options_->GetOpt("particle_flow", "usePFElectrons", usePFElectrons);
01097   cout<<"use PFElectrons "<<usePFElectrons<<endl;
01098 
01099   if( usePFElectrons ) { 
01100     // PFElectrons options -----------------------------
01101     double mvaEleCut = -1.;  // if = -1. get all the pre-id electrons
01102     options_->GetOpt("particle_flow", "electron_mvaCut", mvaEleCut);
01103 
01104     bool applyCrackCorrections=true;
01105     options_->GetOpt("particle_flow","electron_crackCorrection",applyCrackCorrections);
01106 
01107     string mvaWeightFileEleID = "";
01108     options_->GetOpt("particle_flow", "electronID_mvaWeightFile", 
01109                      mvaWeightFileEleID);
01110     mvaWeightFileEleID = expand(mvaWeightFileEleID);
01111     
01112     try { 
01113       pfAlgo_.setPFEleParameters(mvaEleCut,
01114                                  mvaWeightFileEleID,
01115                                  usePFElectrons,
01116                                  thePFSCEnergyCalibration,
01117                                  calibration,
01118                                  sumEtEcalIsoForEgammaSC_barrel,
01119                                  sumEtEcalIsoForEgammaSC_endcap,
01120                                  coneEcalIsoForEgammaSC,
01121                                  sumPtTrackIsoForEgammaSC_barrel,
01122                                  sumPtTrackIsoForEgammaSC_endcap,
01123                                  nTrackIsoForEgammaSC,
01124                                  coneTrackIsoForEgammaSC,
01125                                  applyCrackCorrections,
01126                                  usePFSCEleCalib,
01127                                  useEGammaElectrons,
01128                                  useEGammaSupercluster);
01129     }
01130     catch( std::exception& err ) {
01131       cerr<<"exception setting PFAlgo Electron parameters: "
01132           <<err.what()<<". terminating."<<endl;
01133       delete this;
01134       exit(1);
01135     }
01136   }
01137 
01138   bool usePFPhotons = true;   
01139   string mvaWeightFileConvID = "";
01140   double mvaConvCut=-1.;
01141   options_->GetOpt("particle_flow", "usePFPhotons", usePFPhotons);
01142   options_->GetOpt("particle_flow", "conv_mvaCut", mvaConvCut);
01143   options_->GetOpt("particle_flow", "convID_mvaWeightFile", mvaWeightFileConvID);  
01144   // cout<<"use PFPhotons "<<usePFPhotons<<endl;
01145 
01146   if( usePFPhotons ) { 
01147     // PFPhoton options -----------------------------
01148     
01149     try { 
01150       pfAlgo_.setPFPhotonParameters(usePFPhotons,
01151                                     mvaWeightFileConvID,
01152                                     mvaConvCut,
01153                                     calibration);
01154     }
01155 
01156     catch( std::exception& err ) {
01157       cerr<<"exception setting PFAlgo Photon parameters: "
01158           <<err.what()<<". terminating."<<endl;
01159       delete this;
01160       exit(1);
01161     }
01162   }
01163 
01164 
01165 
01166   bool rejectTracks_Bad = true;
01167   bool rejectTracks_Step45 = true;
01168   bool usePFConversions = false;   // set true to use PFConversions
01169   bool usePFNuclearInteractions = false;
01170   bool usePFV0s = false;
01171 
01172 
01173   double dptRel_DispVtx = 10;
01174   
01175 
01176   options_->GetOpt("particle_flow", "usePFConversions", usePFConversions);
01177   options_->GetOpt("particle_flow", "usePFV0s", usePFV0s);
01178   options_->GetOpt("particle_flow", "usePFNuclearInteractions", usePFNuclearInteractions);
01179   options_->GetOpt("particle_flow", "dptRel_DispVtx",  dptRel_DispVtx);
01180 
01181   try { 
01182     pfAlgo_.setDisplacedVerticesParameters(rejectTracks_Bad,
01183                                            rejectTracks_Step45,
01184                                            usePFNuclearInteractions,
01185                                            usePFConversions,
01186                                            usePFV0s,
01187                                            dptRel_DispVtx);
01188 
01189   }
01190   catch( std::exception& err ) {
01191     cerr<<"exception setting PFAlgo displaced vertex parameters: "
01192         <<err.what()<<". terminating."<<endl;
01193     delete this;
01194     exit(1);
01195   }
01196 
01197   bool bCorrect = false;
01198   bool bCalibPrimary = false;
01199   double dptRel_PrimaryTrack = 0;
01200   double dptRel_MergedTrack = 0;
01201   double ptErrorSecondary = 0;
01202   vector<double> nuclCalibFactors;
01203 
01204   options_->GetOpt("particle_flow", "bCorrect", bCorrect);
01205   options_->GetOpt("particle_flow", "bCalibPrimary", bCalibPrimary);
01206   options_->GetOpt("particle_flow", "dptRel_PrimaryTrack", dptRel_PrimaryTrack);
01207   options_->GetOpt("particle_flow", "dptRel_MergedTrack", dptRel_MergedTrack);
01208   options_->GetOpt("particle_flow", "ptErrorSecondary", ptErrorSecondary);
01209   options_->GetOpt("particle_flow", "nuclCalibFactors", nuclCalibFactors);
01210 
01211   try { 
01212     pfAlgo_.setCandConnectorParameters(bCorrect, bCalibPrimary, dptRel_PrimaryTrack, dptRel_MergedTrack, ptErrorSecondary, nuclCalibFactors);
01213   }
01214   catch( std::exception& err ) {
01215     cerr<<"exception setting PFAlgo cand connector parameters: "
01216         <<err.what()<<". terminating."<<endl;
01217     delete this;
01218     exit(1);
01219   }
01220 
01221 
01222 
01223 
01224   int    algo = 2;
01225   options_->GetOpt("particle_flow", "algorithm", algo);
01226 
01227   pfAlgo_.setAlgo( algo );
01228   //   pfAlgoOther_.setAlgo( 1 );
01229 
01230 
01231   // jets options ---------------------------------
01232 
01233   doJets_ = false;
01234   options_->GetOpt("jets", "on/off", doJets_);
01235 
01236   jetsDebug_ = false;
01237   options_->GetOpt("jets", "debug", jetsDebug_);
01238 
01239   jetAlgoType_=3; //FastJet as Default
01240   options_->GetOpt("jets", "algo", jetAlgoType_);
01241 
01242   double mEtInputCut = 0.5;
01243   options_->GetOpt("jets", "EtInputCut",  mEtInputCut);           
01244 
01245   double mEInputCut = 0.;
01246   options_->GetOpt("jets", "EInputCut",  mEInputCut);  
01247 
01248   double seedThreshold  = 1.0;
01249   options_->GetOpt("jets", "seedThreshold", seedThreshold);
01250 
01251   double coneRadius = 0.5;
01252   options_->GetOpt("jets", "coneRadius", coneRadius);             
01253 
01254   double coneAreaFraction= 1.0;
01255   options_->GetOpt("jets", "coneAreaFraction",  coneAreaFraction);   
01256 
01257   int maxPairSize=2;
01258   options_->GetOpt("jets", "maxPairSize",  maxPairSize);  
01259 
01260   int maxIterations=100;
01261   options_->GetOpt("jets", "maxIterations",  maxIterations);      
01262 
01263   double overlapThreshold  = 0.75;
01264   options_->GetOpt("jets", "overlapThreshold", overlapThreshold);
01265 
01266   double ptMin = 10.;
01267   options_->GetOpt("jets", "ptMin",  ptMin);      
01268 
01269   double rparam = 1.0;
01270   options_->GetOpt("jets", "rParam",  rparam);    
01271  
01272   jetMaker_.setmEtInputCut (mEtInputCut);
01273   jetMaker_.setmEInputCut(mEInputCut); 
01274   jetMaker_.setSeedThreshold(seedThreshold); 
01275   jetMaker_.setConeRadius(coneRadius);
01276   jetMaker_.setConeAreaFraction(coneAreaFraction);
01277   jetMaker_.setMaxPairSize(maxPairSize);
01278   jetMaker_.setMaxIterations(maxIterations) ;
01279   jetMaker_.setOverlapThreshold(overlapThreshold) ;
01280   jetMaker_.setPtMin (ptMin);
01281   jetMaker_.setRParam (rparam);
01282   jetMaker_.setDebug(jetsDebug_);
01283   jetMaker_.updateParameter();
01284   cout <<"Opt: doJets? " << doJets_  <<endl; 
01285   cout <<"Opt: jetsDebug " << jetsDebug_  <<endl; 
01286   cout <<"Opt: algoType " << jetAlgoType_  <<endl; 
01287   cout <<"----------------------------------" << endl;
01288 
01289 
01290   // tau benchmark options ---------------------------------
01291 
01292   doTauBenchmark_ = false;
01293   options_->GetOpt("tau_benchmark", "on/off", doTauBenchmark_);
01294   
01295   if (doTauBenchmark_) {
01296     double coneAngle = 0.5;
01297     options_->GetOpt("tau_benchmark", "cone_angle", coneAngle);
01298     
01299     double seedEt    = 0.4;
01300     options_->GetOpt("tau_benchmark", "seed_et", seedEt);
01301     
01302     double coneMerge = 100.0;
01303     options_->GetOpt("tau_benchmark", "cone_merge", coneMerge);
01304     
01305     options_->GetOpt("tau_benchmark", "debug", tauBenchmarkDebug_);
01306 
01307     // cout<<"jets debug "<<jetsDebug_<<endl;
01308     
01309     if( tauBenchmarkDebug_ ) {
01310       cout << "Tau Benchmark Options : ";
01311       cout << "Angle=" << coneAngle << " seedEt=" << seedEt 
01312            << " Merge=" << coneMerge << endl;
01313     }
01314 
01315     jetAlgo_.SetConeAngle(coneAngle);
01316     jetAlgo_.SetSeedEt(seedEt);
01317     jetAlgo_.SetConeMerge(coneMerge);   
01318   }
01319 
01320 
01321 
01322   // print flags -------------
01323 
01324   printRecHits_ = false;
01325   printRecHitsEMin_ = 0.;
01326   options_->GetOpt("print", "rechits", printRecHits_ );
01327   options_->GetOpt("print", "rechits_emin", printRecHitsEMin_ );
01328   
01329   printClusters_ = false;
01330   printClustersEMin_ = 0.;
01331   options_->GetOpt("print", "clusters", printClusters_ );
01332   options_->GetOpt("print", "clusters_emin", printClustersEMin_ );
01333 
01334   printPFBlocks_ = false;
01335   options_->GetOpt("print", "PFBlocks", printPFBlocks_ );
01336   
01337   printPFCandidates_ = true;
01338   printPFCandidatesPtMin_ = 0.;
01339   options_->GetOpt("print", "PFCandidates", printPFCandidates_ );
01340   options_->GetOpt("print", "PFCandidates_ptmin", printPFCandidatesPtMin_ );
01341   
01342   printPFJets_ = true;
01343   printPFJetsPtMin_ = 0.;
01344   options_->GetOpt("print", "jets", printPFJets_ );
01345   options_->GetOpt("print", "jets_ptmin", printPFJetsPtMin_ );
01346  
01347   printSimParticles_ = true;
01348   printSimParticlesPtMin_ = 0.;
01349   options_->GetOpt("print", "simParticles", printSimParticles_ );
01350   options_->GetOpt("print", "simParticles_ptmin", printSimParticlesPtMin_ );
01351 
01352   printGenParticles_ = true;
01353   printGenParticlesPtMin_ = 0.;
01354   options_->GetOpt("print", "genParticles", printGenParticles_ );
01355   options_->GetOpt("print", "genParticles_ptmin", printGenParticlesPtMin_ );
01356 
01357   //MCTruthMatching Tool set to false by default
01358   //can only be used with fastsim and the UnFoldedMode set to true
01359   //when generating the simulated file
01360   printMCTruthMatching_ = false; 
01361   options_->GetOpt("print", "mctruthmatching", printMCTruthMatching_ );  
01362 
01363 
01364   verbosity_ = VERBOSE;
01365   options_->GetOpt("print", "verbosity", verbosity_ );
01366   cout<<"verbosity : "<<verbosity_<<endl;
01367 
01368 
01369 
01370   
01371 
01372 }
01373 
01374 void PFRootEventManager::connect( const char* infilename ) {
01375 
01376   cout<<"Opening input root files"<<endl;
01377 
01378   options_->GetOpt("root","file", inFileNames_);
01379   
01380 
01381 
01382   try {
01383     AutoLibraryLoader::enable();
01384   }
01385   catch(string& err) {
01386     cout<<err<<endl;
01387   }
01388 
01389   ev_ = new fwlite::ChainEvent(inFileNames_);
01390 
01391 
01392   if ( !ev_ || !ev_->isValid() ) { 
01393     cout << "The rootfile(s) " << endl;
01394     for ( unsigned int ifile=0; ifile<inFileNames_.size(); ++ifile ) 
01395       std::cout << " - " << inFileNames_[ifile] << std::endl;
01396     cout << " is (are) not valid file(s) to open" << endl;
01397     return;
01398   } else { 
01399     cout << "The rootfile(s) : " << endl;
01400     for ( unsigned int ifile=0; ifile<inFileNames_.size(); ++ifile ) 
01401       std::cout << " - " << inFileNames_[ifile] << std::endl;
01402     cout<<" are opened with " << ev_->size() << " events." <<endl;
01403   }
01404   
01405   // hits branches ----------------------------------------------
01406   std::string rechitsECALtagname;
01407   options_->GetOpt("root","rechits_ECAL_inputTag", rechitsECALtagname);
01408   rechitsECALTag_ = edm::InputTag(rechitsECALtagname);
01409 
01410   std::string rechitsHCALtagname;
01411   options_->GetOpt("root","rechits_HCAL_inputTag", rechitsHCALtagname);
01412   rechitsHCALTag_ = edm::InputTag(rechitsHCALtagname);
01413 
01414   std::string rechitsHFEMtagname;
01415   options_->GetOpt("root","rechits_HFEM_inputTag", rechitsHFEMtagname);
01416   rechitsHFEMTag_ = edm::InputTag(rechitsHFEMtagname);
01417 
01418   std::string rechitsHFHADtagname;
01419   options_->GetOpt("root","rechits_HFHAD_inputTag", rechitsHFHADtagname);
01420   rechitsHFHADTag_ = edm::InputTag(rechitsHFHADtagname);
01421 
01422   std::vector<string> rechitsCLEANEDtagnames;
01423   options_->GetOpt("root","rechits_CLEANED_inputTags", rechitsCLEANEDtagnames);
01424   for ( unsigned tags = 0; tags<rechitsCLEANEDtagnames.size(); ++tags )
01425     rechitsCLEANEDTags_.push_back(edm::InputTag(rechitsCLEANEDtagnames[tags]));
01426   rechitsCLEANEDV_.resize(rechitsCLEANEDTags_.size());
01427   rechitsCLEANEDHandles_.resize(rechitsCLEANEDTags_.size());
01428 
01429 
01430   // Tracks branches
01431   std::string rechitsPStagname;
01432   options_->GetOpt("root","rechits_PS_inputTag", rechitsPStagname);
01433   rechitsPSTag_ = edm::InputTag(rechitsPStagname);
01434 
01435   std::string recTrackstagname;
01436   options_->GetOpt("root","recTracks_inputTag", recTrackstagname);
01437   recTracksTag_ = edm::InputTag(recTrackstagname);
01438 
01439   std::string displacedRecTrackstagname;
01440   options_->GetOpt("root","displacedRecTracks_inputTag", displacedRecTrackstagname);
01441   displacedRecTracksTag_ = edm::InputTag(displacedRecTrackstagname);
01442 
01443   std::string primaryVerticestagname;
01444   options_->GetOpt("root","primaryVertices_inputTag", primaryVerticestagname);
01445   primaryVerticesTag_ = edm::InputTag(primaryVerticestagname);
01446 
01447   std::string stdTrackstagname;
01448   options_->GetOpt("root","stdTracks_inputTag", stdTrackstagname);
01449   stdTracksTag_ = edm::InputTag(stdTrackstagname);
01450 
01451   std::string gsfrecTrackstagname;
01452   options_->GetOpt("root","gsfrecTracks_inputTag", gsfrecTrackstagname);
01453   gsfrecTracksTag_ = edm::InputTag(gsfrecTrackstagname);
01454 
01455   useConvBremGsfTracks_ = false;
01456   options_->GetOpt("particle_flow", "useConvBremGsfTracks", useConvBremGsfTracks_);
01457   if ( useConvBremGsfTracks_ ) { 
01458     std::string convBremGsfrecTrackstagname;
01459     options_->GetOpt("root","convBremGsfrecTracks_inputTag", convBremGsfrecTrackstagname);
01460     convBremGsfrecTracksTag_ = edm::InputTag(convBremGsfrecTrackstagname);
01461   }
01462 
01463   useConvBremPFRecTracks_ = false;
01464   options_->GetOpt("particle_flow", "useConvBremPFRecTracks", useConvBremPFRecTracks_);
01465 
01466 
01467   // muons branch
01468   std::string muonstagname;
01469   options_->GetOpt("root","muon_inputTag", muonstagname);
01470   muonsTag_ = edm::InputTag(muonstagname);
01471 
01472   // conversion
01473   usePFConversions_=false;
01474   options_->GetOpt("particle_flow", "usePFConversions", usePFConversions_);
01475   if( usePFConversions_ ) {
01476     std::string conversiontagname;
01477     options_->GetOpt("root","conversion_inputTag", conversiontagname);
01478     conversionTag_ = edm::InputTag(conversiontagname);
01479   }
01480 
01481   // V0
01482   usePFV0s_=false;
01483   options_->GetOpt("particle_flow", "usePFV0s", usePFV0s_);
01484   if( usePFV0s_ ) {
01485     std::string v0tagname;
01486     options_->GetOpt("root","V0_inputTag", v0tagname);
01487     v0Tag_ = edm::InputTag(v0tagname);
01488   }
01489 
01490   // Photons
01491   std::string photontagname;
01492   options_->GetOpt("root","Photon_inputTag",photontagname);
01493   photonTag_ = edm::InputTag(photontagname);
01494 
01495  //Displaced Vertices
01496   usePFNuclearInteractions_=false;
01497   options_->GetOpt("particle_flow", "usePFNuclearInteractions", usePFNuclearInteractions_);
01498   if( usePFNuclearInteractions_ ) {
01499     std::string pfNuclearTrackerVertextagname;
01500     options_->GetOpt("root","PFDisplacedVertex_inputTag", pfNuclearTrackerVertextagname);
01501     pfNuclearTrackerVertexTag_ = edm::InputTag(pfNuclearTrackerVertextagname);
01502   }
01503 
01504   std::string trueParticlestagname;
01505   options_->GetOpt("root","trueParticles_inputTag", trueParticlestagname);
01506   trueParticlesTag_ = edm::InputTag(trueParticlestagname);
01507 
01508   std::string MCTruthtagname;
01509   options_->GetOpt("root","MCTruth_inputTag", MCTruthtagname);
01510   MCTruthTag_ = edm::InputTag(MCTruthtagname);
01511 
01512   std::string caloTowerstagname;
01513   options_->GetOpt("root","caloTowers_inputTag", caloTowerstagname);
01514   caloTowersTag_ = edm::InputTag(caloTowerstagname);
01515 
01516   std::string genJetstagname;
01517   options_->GetOpt("root","genJets_inputTag", genJetstagname);
01518   genJetsTag_ = edm::InputTag(genJetstagname);
01519 
01520   
01521   std::string genParticlesforMETtagname;
01522   options_->GetOpt("root","genParticlesforMET_inputTag", genParticlesforMETtagname);
01523   genParticlesforMETTag_ = edm::InputTag(genParticlesforMETtagname);
01524 
01525   std::string genParticlesforJetstagname;
01526   options_->GetOpt("root","genParticlesforJets_inputTag", genParticlesforJetstagname);
01527   genParticlesforJetsTag_ = edm::InputTag(genParticlesforJetstagname);
01528 
01529   // PF candidates 
01530   std::string pfCandidatetagname;
01531   options_->GetOpt("root","particleFlowCand_inputTag", pfCandidatetagname);
01532   pfCandidateTag_ = edm::InputTag(pfCandidatetagname);
01533 
01534   std::string caloJetstagname;
01535   options_->GetOpt("root","CaloJets_inputTag", caloJetstagname);
01536   caloJetsTag_ = edm::InputTag(caloJetstagname);
01537 
01538   std::string corrcaloJetstagname;
01539   options_->GetOpt("root","corrCaloJets_inputTag", corrcaloJetstagname);
01540   corrcaloJetsTag_ = edm::InputTag(corrcaloJetstagname);
01541 
01542   std::string pfJetstagname;
01543   options_->GetOpt("root","PFJets_inputTag", pfJetstagname);
01544   pfJetsTag_ = edm::InputTag(pfJetstagname);
01545 
01546   std::string pfMetstagname;
01547   options_->GetOpt("root","PFMET_inputTag", pfMetstagname);
01548   pfMetsTag_ = edm::InputTag(pfMetstagname);
01549 
01550   std::string caloMetstagname;
01551   options_->GetOpt("root","CaloMET_inputTag", caloMetstagname);
01552   caloMetsTag_ = edm::InputTag(caloMetstagname);
01553 
01554   std::string tcMetstagname;
01555   options_->GetOpt("root","TCMET_inputTag", tcMetstagname);
01556   tcMetsTag_ = edm::InputTag(tcMetstagname);
01557 
01558 }
01559 
01560 
01561 
01562 
01563 
01564 PFRootEventManager::~PFRootEventManager() {
01565 
01566   if(outFile_) {
01567     outFile_->Close();
01568   }
01569 
01570   if(outEvent_) delete outEvent_;
01571 
01572   delete options_;
01573 
01574 }
01575 
01576 
01577 void PFRootEventManager::write() {
01578 
01579   if(doPFJetBenchmark_) PFJetBenchmark_.write();
01580   if(doPFMETBenchmark_) metManager_->write();
01581   clusterAlgoECAL_.write();
01582   clusterAlgoHCAL_.write();
01583   clusterAlgoPS_.write();
01584   clusterAlgoHFEM_.write();
01585   clusterAlgoHFHAD_.write();
01586   
01587   // Addition to have DQM histograms : by S. Dutta                                                                                                       
01588   if (doPFDQM_) {
01589     cout << " Writing DQM root file " << endl;
01590     pfJetMonitor_.write();
01591     pfMETMonitor_.write();
01592     dqmFile_->Write();
01593   }
01594   //-----------------------------------------------                                                                                                           
01595   if(outFile_) {
01596     outFile_->Write();
01597 //     outFile_->cd(); 
01598 //     // write histos here
01599 //     cout<<"writing output to "<<outFile_->GetName()<<endl;
01600 //     h_deltaETvisible_MCEHT_->Write();
01601 //     h_deltaETvisible_MCPF_->Write();
01602 //     if(outTree_) outTree_->Write();
01603 //     if(doPFCandidateBenchmark_) pfCandidateBenchmark_.write();
01604   }
01605 }
01606 
01607 
01608 int PFRootEventManager::eventToEntry(int run, int lumi, int event) const {
01609   
01610   RunsMap::const_iterator iR = mapEventToEntry_.find( run );
01611   if( iR != mapEventToEntry_.end() ) {
01612     LumisMap::const_iterator iL = iR->second.find( lumi );
01613     if( iL != iR->second.end() ) {
01614       EventToEntry::const_iterator iE = iL->second.find( event );
01615       if( iE != iL->second.end() ) {
01616         return iE->second;
01617       }  
01618       else {
01619         cout<<"event "<<event<<" not found in run "<<run<<", lumi "<<lumi<<endl;
01620       }
01621     }
01622     else {
01623       cout<<"lumi "<<lumi<<" not found in run "<<run<<endl;
01624     }
01625   }
01626   else{
01627     cout<<"run "<<run<<" not found"<<endl;
01628   }
01629   return -1;    
01630 }
01631 
01632 bool PFRootEventManager::processEvent(int run, int lumi, int event) {
01633 
01634   int entry = eventToEntry(run, lumi, event);
01635   if( entry < 0 ) {
01636     cout<<"event "<<event<<" is not present, sorry."<<endl;
01637     return false;
01638   }
01639   else
01640     return processEntry( entry ); 
01641 } 
01642 
01643 
01644 bool PFRootEventManager::processEntry(int entry) {
01645 
01646   reset();
01647 
01648   iEvent_ = entry;
01649  
01650   bool exists = ev_->to(entry);
01651   if ( !exists ) { 
01652     std::cout << "Entry " << entry << " does not exist " << std::endl; 
01653     return false;
01654   }
01655   const edm::EventBase& iEvent = *ev_;
01656 
01657   if( outEvent_ ) outEvent_->setNumber(entry);
01658 
01659   if(verbosity_ == VERBOSE  || 
01660      //entry < 10000 ||
01661      (entry < 100 && entry%10 == 0) || 
01662      (entry < 1000 && entry%100 == 0) || 
01663      entry%1000 == 0 ) 
01664     cout<<"process entry "<< entry 
01665         <<", run "<<iEvent.id().run()
01666         <<", lumi "<<iEvent.id().luminosityBlock()      
01667         <<", event:"<<iEvent.id().event()
01668         << endl;
01669 
01670   //ev_->getTFile()->cd();
01671   bool goodevent =  readFromSimulation(entry);
01672 
01673   /* 
01674   std::cout << "Rechits cleaned : " << std::endl;
01675   for(unsigned i=0; i<rechitsCLEANED_.size(); i++) {
01676     string seedstatus = "    ";
01677     printRecHit(rechitsCLEANED_[i], i, seedstatus.c_str());
01678   }
01679   */
01680 
01681   if(verbosity_ == VERBOSE ) {
01682     cout<<"number of vertices             : "<<primaryVertices_.size()<<endl;
01683     cout<<"number of recTracks            : "<<recTracks_.size()<<endl;
01684     cout<<"number of gsfrecTracks         : "<<gsfrecTracks_.size()<<endl;
01685     cout<<"number of convBremGsfrecTracks : "<<convBremGsfrecTracks_.size()<<endl;
01686     cout<<"number of muons                : "<<muons_.size()<<endl;
01687     cout<<"number of displaced vertices   : "<<pfNuclearTrackerVertex_.size()<<endl;
01688     cout<<"number of conversions          : "<<conversion_.size()<<endl;
01689     cout<<"number of v0                   : "<<v0_.size()<<endl;
01690     cout<<"number of stdTracks            : "<<stdTracks_.size()<<endl;
01691     cout<<"number of true particles       : "<<trueParticles_.size()<<endl;
01692     cout<<"number of ECAL rechits         : "<<rechitsECAL_.size()<<endl;
01693     cout<<"number of HCAL rechits         : "<<rechitsHCAL_.size()<<endl;
01694     cout<<"number of HFEM rechits         : "<<rechitsHFEM_.size()<<endl;
01695     cout<<"number of HFHAD rechits        : "<<rechitsHFHAD_.size()<<endl;
01696     cout<<"number of HF Cleaned rechits   : "<<rechitsCLEANED_.size()<<endl;
01697     cout<<"number of PS rechits           : "<<rechitsPS_.size()<<endl;
01698   }  
01699 
01700   if( doClustering_ ) clustering(); 
01701   else if( verbosity_ == VERBOSE )
01702     cout<<"clustering is OFF - clusters come from the input file"<<endl; 
01703 
01704   if(verbosity_ == VERBOSE ) {
01705     if(clustersECAL_.get() ) {
01706       cout<<"number of ECAL clusters : "<<clustersECAL_->size()<<endl;
01707     }
01708     if(clustersHCAL_.get() ) {
01709       cout<<"number of HCAL clusters : "<<clustersHCAL_->size()<<endl;
01710     }
01711     if(clustersHFEM_.get() ) {
01712       cout<<"number of HFEM clusters : "<<clustersHFEM_->size()<<endl;
01713     }
01714     if(clustersHFHAD_.get() ) {
01715       cout<<"number of HFHAD clusters : "<<clustersHFHAD_->size()<<endl;
01716     }
01717     if(clustersPS_.get() ) {
01718       cout<<"number of PS clusters : "<<clustersPS_->size()<<endl;
01719     }
01720   }
01721 
01722   if(doParticleFlow_) { 
01723     particleFlow();
01724 
01725     if (doCompare_) pfCandCompare(entry);
01726   }
01727 
01728   if(doJets_) {
01729     reconstructGenJets();
01730     reconstructCaloJets();
01731     reconstructPFJets();
01732   }    
01733         
01734  
01735   // call print() in verbose mode
01736   if( verbosity_ == VERBOSE ) print();
01737   
01738   //COLIN the PFJet and PFMET benchmarks are very messy. 
01739   //COLIN    compare with the filling of the PFCandidateBenchmark, which is one line. 
01740   
01741   goodevent = eventAccepted(); 
01742 
01743   // evaluate PFJet Benchmark 
01744   if(doPFJetBenchmark_) { // start PFJet Benchmark
01745 
01746     PFJetBenchmark_.process(pfJets_, genJets_);
01747     double resPt = PFJetBenchmark_.resPtMax();
01748     double resChargedHadEnergy = PFJetBenchmark_.resChargedHadEnergyMax();
01749     double resNeutralHadEnergy = PFJetBenchmark_.resNeutralHadEnergyMax();
01750     double resNeutralEmEnergy = PFJetBenchmark_.resNeutralEmEnergyMax();
01751           
01752     if( verbosity_ == VERBOSE ){ //start debug print
01753 
01754       cout << " =====================PFJetBenchmark =================" << endl;
01755       cout<<"Resol Pt max "<<resPt
01756           <<" resChargedHadEnergy Max " << resChargedHadEnergy
01757           <<" resNeutralHadEnergy Max " << resNeutralHadEnergy
01758           << " resNeutralEmEnergy Max "<< resNeutralEmEnergy << endl;
01759     } // end debug print
01760 
01761     // PJ : printout for bad events (selected by the "if")
01762     /*
01763     if ( fabs(resPt) > 0.4 )
01764       std::cout << "'" << iEvent.id().run() << ":" << iEvent.id().event() << "-" 
01765                 << iEvent.id().run() << ":" << iEvent.id().event() << "'," << std::endl;
01766     */
01767     if ( resPt < -1. ) { 
01768       cout << " =====================PFJetBenchmark =================" << endl;
01769       cout<<"process entry "<< entry << endl;
01770       cout<<"Resol Pt max "<<resPt
01771           <<" resChargedHadEnergy Max " << resChargedHadEnergy
01772           <<" resNeutralHadEnergy Max " << resNeutralHadEnergy
01773           << " resNeutralEmEnergy Max "<< resNeutralEmEnergy 
01774           << " Jet pt " << genJets_[0].pt() << endl;
01775       // return true;
01776     } else { 
01777       // return false;
01778     }
01779     //   if (resNeutralEmEnergy>0.5) return true;
01780     //   else return false;
01781   }// end PFJet Benchmark
01782 
01783   // Addition to have DQM histograms : by S. Dutta 
01784   reco::MET reComputedMet_;    
01785   reco::MET computedGenMet_;
01786   //-----------------------------------------------
01787 
01788   //COLIN would  be nice to move this long code to a separate function. 
01789   // is it necessary to re-set everything at each event?? 
01790   if(doPFMETBenchmark_) { // start PFMet Benchmark
01791 
01792     // Fill here the various met benchmarks
01793     // pfMET vs GenMET
01794     metManager_->setMET1(&genParticlesCMSSW_);
01795     metManager_->setMET2(&pfMetsCMSSW_[0]);
01796     metManager_->FillHisto("PF");
01797     // cout events in tail
01798     metManager_->coutTailEvents(entry,DeltaMETcut,DeltaPhicut, MET1cut);
01799 
01800     // caloMET vs GenMET
01801     metManager_->setMET2(&caloMetsCMSSW_[0]);
01802     metManager_->FillHisto("Calo");
01803 
01804     if ( doMet_ ) { 
01805       // recomputed pfMET vs GenMET
01806       metManager_->setMET2(*pfCandidates_);
01807       metManager_->FillHisto("recompPF");
01808       metManager_->coutTailEvents(entry,DeltaMETcut,DeltaPhicut, MET1cut);
01809     }
01810 
01811     if (JECinCaloMet_)
01812     {
01813       // corrCaloMET vs GenMET
01814       metManager_->setMET2(&caloMetsCMSSW_[0]);
01815       metManager_->propagateJECtoMET2(caloJetsCMSSW_, corrcaloJetsCMSSW_);
01816       metManager_->FillHisto("corrCalo");
01817     }
01818   }// end PFMET Benchmark
01819 
01820   if( goodevent && doPFCandidateBenchmark_ ) {
01821     pfCandidateManager_.fill( *pfCandidates_, genParticlesCMSSW_);
01822   }
01823   
01824   // Addition to have DQM histograms : by S. Dutta                                                                                                          
01825   if( goodevent && doPFDQM_ ) {
01826     float deltaMin, deltaMax;
01827     pfJetMonitor_.fill( pfJets_, genJets_, deltaMin, deltaMax);
01828     if (doPFMETBenchmark_) {
01829       pfMETMonitor_.fillOne( reComputedMet_, computedGenMet_, deltaMin, deltaMax);
01830     }
01831   }
01832   //-----------------------------------------------                                                                    
01833   // evaluate tau Benchmark   
01834   if( goodevent && doTauBenchmark_) { // start tau Benchmark
01835     double deltaEt = 0.;
01836     deltaEt  = tauBenchmark( *pfCandidates_ ); 
01837     if( verbosity_ == VERBOSE ) cout<<"delta E_t ="<<deltaEt <<endl;
01838     //      cout<<"delta E_t ="<<deltaEt<<" delta E_t Other ="<<deltaEt1<<endl;
01839 
01840 
01841     //   if( deltaEt>0.4 ) {
01842     //     cout<<deltaEt<<endl;
01843     //     return true;
01844     //   }  
01845     //   else return false;
01846 
01847   
01848   } // end tau Benchmark
01849 
01850   if(goodevent && outTree_) 
01851     outTree_->Fill();
01852 
01853   if(calibFile_)
01854     printMCCalib(*calibFile_);
01855   
01856   return goodevent;
01857 
01858 }
01859 
01860 
01861 
01862 bool PFRootEventManager::eventAccepted() const {
01863   // return highPtJet(10); 
01864   //return highPtPFCandidate( 10, PFCandidate::h ); 
01865   return true;
01866 } 
01867 
01868 
01869 bool PFRootEventManager::highPtJet(double ptMin) const {
01870   for( unsigned i=0; i<pfJets_.size(); ++i) {
01871     if( pfJets_[i].pt() > ptMin ) return true;
01872   }
01873   return false;
01874 }
01875 
01876 bool PFRootEventManager::highPtPFCandidate( double ptMin, 
01877                                             PFCandidate::ParticleType type) const {
01878   for( unsigned i=0; i<pfCandidates_->size(); ++i) {
01879 
01880     const PFCandidate& pfc = (*pfCandidates_)[i];
01881     if(type!= PFCandidate::X &&  
01882        pfc.particleId() != type ) continue;
01883     if( pfc.pt() > ptMin ) return true;
01884   }
01885   return false;
01886 }
01887 
01888 
01889 bool PFRootEventManager::readFromSimulation(int entry) {
01890 
01891   if (verbosity_ == VERBOSE ) {
01892     cout <<"start reading from simulation"<<endl;
01893   }
01894 
01895 
01896   // if(!tree_) return false;
01897   
01898   const edm::EventBase& iEvent = *ev_;
01899   
01900 
01901   bool foundstdTracks = iEvent.getByLabel(stdTracksTag_,stdTracksHandle_);
01902   if ( foundstdTracks ) { 
01903     stdTracks_ = *stdTracksHandle_;
01904     // cout << "Found " << stdTracks_.size() << " standard tracks" << endl;
01905   } else { 
01906     cerr<<"PFRootEventManager::ProcessEntry : stdTracks Collection not found : "
01907         <<entry << " " << stdTracksTag_<<endl;
01908   }
01909 
01910   bool foundMCTruth = iEvent.getByLabel(MCTruthTag_,MCTruthHandle_);
01911   if ( foundMCTruth ) { 
01912     MCTruth_ = *MCTruthHandle_;
01913     // cout << "Found MC truth" << endl;
01914   } else { 
01915     // cerr<<"PFRootEventManager::ProcessEntry : MCTruth Collection not found : "
01916     //    <<entry << " " << MCTruthTag_<<endl;
01917   }
01918 
01919   bool foundTP = iEvent.getByLabel(trueParticlesTag_,trueParticlesHandle_);
01920   if ( foundTP ) { 
01921     trueParticles_ = *trueParticlesHandle_;
01922     // cout << "Found " << trueParticles_.size() << " true particles" << endl;
01923   } else { 
01924     //cerr<<"PFRootEventManager::ProcessEntry : trueParticles Collection not found : "
01925     //    <<entry << " " << trueParticlesTag_<<endl;
01926   }
01927 
01928   bool foundECAL = iEvent.getByLabel(rechitsECALTag_,rechitsECALHandle_);
01929   if ( foundECAL ) { 
01930     rechitsECAL_ = *rechitsECALHandle_;
01931     // cout << "Found " << rechitsECAL_.size() << " ECAL rechits" << endl;
01932   } else { 
01933     cerr<<"PFRootEventManager::ProcessEntry : rechitsECAL Collection not found : "
01934         <<entry << " " << rechitsECALTag_<<endl;
01935   }
01936 
01937   bool foundHCAL = iEvent.getByLabel(rechitsHCALTag_,rechitsHCALHandle_);
01938   if ( foundHCAL ) { 
01939     rechitsHCAL_ = *rechitsHCALHandle_;
01940     // cout << "Found " << rechitsHCAL_.size() << " HCAL rechits" << endl;
01941   } else { 
01942     cerr<<"PFRootEventManager::ProcessEntry : rechitsHCAL Collection not found : "
01943         <<entry << " " << rechitsHCALTag_<<endl;
01944   }
01945 
01946   bool foundHFEM = iEvent.getByLabel(rechitsHFEMTag_,rechitsHFEMHandle_);
01947   if ( foundHFEM ) { 
01948     rechitsHFEM_ = *rechitsHFEMHandle_;
01949     // cout << "Found " << rechitsHFEM_.size() << " HFEM rechits" << endl;
01950   } else { 
01951     cerr<<"PFRootEventManager::ProcessEntry : rechitsHFEM Collection not found : "
01952         <<entry << " " << rechitsHFEMTag_<<endl;
01953   }
01954 
01955   bool foundHFHAD = iEvent.getByLabel(rechitsHFHADTag_,rechitsHFHADHandle_);
01956   if ( foundHFHAD ) { 
01957     rechitsHFHAD_ = *rechitsHFHADHandle_;
01958     // cout << "Found " << rechitsHFHAD_.size() << " HFHAD rechits" << endl;
01959   } else { 
01960     cerr<<"PFRootEventManager::ProcessEntry : rechitsHFHAD Collection not found : "
01961         <<entry << " " << rechitsHFHADTag_<<endl;
01962   }
01963 
01964   for ( unsigned i=0; i<rechitsCLEANEDTags_.size(); ++i ) { 
01965     bool foundCLEANED = iEvent.getByLabel(rechitsCLEANEDTags_[i],
01966                                           rechitsCLEANEDHandles_[i]);
01967     if ( foundCLEANED ) { 
01968       rechitsCLEANEDV_[i] = *(rechitsCLEANEDHandles_[i]);
01969       // cout << "Found " << rechitsCLEANEDV_[i].size() << " CLEANED rechits" << endl;
01970     } else { 
01971       cerr<<"PFRootEventManager::ProcessEntry : rechitsCLEANED Collection not found : "
01972           <<entry << " " << rechitsCLEANEDTags_[i]<<endl;
01973     }
01974 
01975   }
01976 
01977   bool foundPS = iEvent.getByLabel(rechitsPSTag_,rechitsPSHandle_);
01978   if ( foundPS ) { 
01979     rechitsPS_ = *rechitsPSHandle_;
01980     // cout << "Found " << rechitsPS_.size() << " PS rechits" << endl;
01981   } else { 
01982     cerr<<"PFRootEventManager::ProcessEntry : rechitsPS Collection not found : "
01983         <<entry << " " << rechitsPSTag_<<endl;
01984   }
01985 
01986   bool foundCT = iEvent.getByLabel(caloTowersTag_,caloTowersHandle_);
01987   if ( foundCT ) { 
01988     caloTowers_ = *caloTowersHandle_;
01989     // cout << "Found " << caloTowers_.size() << " calo Towers" << endl;
01990   } else { 
01991     cerr<<"PFRootEventManager::ProcessEntry : caloTowers Collection not found : "
01992         <<entry << " " << caloTowersTag_<<endl;
01993   }
01994 
01995   bool foundPV = iEvent.getByLabel(primaryVerticesTag_,primaryVerticesHandle_);
01996   if ( foundPV ) { 
01997     primaryVertices_ = *primaryVerticesHandle_;
01998     // cout << "Found " << primaryVertices_.size() << " primary vertices" << endl;
01999   } else { 
02000     cerr<<"PFRootEventManager::ProcessEntry : primaryVertices Collection not found : "
02001         <<entry << " " << primaryVerticesTag_<<endl;
02002   }
02003 
02004   bool foundPFV = iEvent.getByLabel(pfNuclearTrackerVertexTag_,pfNuclearTrackerVertexHandle_);
02005   if ( foundPFV ) { 
02006     pfNuclearTrackerVertex_ = *pfNuclearTrackerVertexHandle_;
02007     // cout << "Found " << pfNuclearTrackerVertex_.size() << " secondary PF vertices" << endl;
02008   } else if ( usePFNuclearInteractions_ ) { 
02009     cerr<<"PFRootEventManager::ProcessEntry : pfNuclearTrackerVertex Collection not found : "
02010         <<entry << " " << pfNuclearTrackerVertexTag_<<endl;
02011   }
02012 
02013   bool foundrecTracks = iEvent.getByLabel(recTracksTag_,recTracksHandle_);
02014   if ( foundrecTracks ) { 
02015     recTracks_ = *recTracksHandle_;
02016     // cout << "Found " << recTracks_.size() << " PFRecTracks" << endl;
02017   } else { 
02018     cerr<<"PFRootEventManager::ProcessEntry : recTracks Collection not found : "
02019         <<entry << " " << recTracksTag_<<endl;
02020   }
02021 
02022   bool founddisplacedRecTracks = iEvent.getByLabel(displacedRecTracksTag_,displacedRecTracksHandle_);
02023   if ( founddisplacedRecTracks ) { 
02024     displacedRecTracks_ = *displacedRecTracksHandle_;
02025     // cout << "Found " << displacedRecTracks_.size() << " PFRecTracks" << endl;
02026   } else { 
02027     cerr<<"PFRootEventManager::ProcessEntry : displacedRecTracks Collection not found : "
02028         <<entry << " " << displacedRecTracksTag_<<endl;
02029   }
02030 
02031 
02032   bool foundgsfrecTracks = iEvent.getByLabel(gsfrecTracksTag_,gsfrecTracksHandle_);
02033   if ( foundgsfrecTracks ) { 
02034     gsfrecTracks_ = *gsfrecTracksHandle_;
02035     // cout << "Found " << gsfrecTracks_.size() << " GsfPFRecTracks" << endl;
02036   } else { 
02037     cerr<<"PFRootEventManager::ProcessEntry : gsfrecTracks Collection not found : "
02038         <<entry << " " << gsfrecTracksTag_<<endl;
02039   }
02040 
02041   bool foundconvBremGsfrecTracks = iEvent.getByLabel(convBremGsfrecTracksTag_,convBremGsfrecTracksHandle_);
02042   if ( foundconvBremGsfrecTracks ) { 
02043     convBremGsfrecTracks_ = *convBremGsfrecTracksHandle_;
02044     // cout << "Found " << convBremGsfrecTracks_.size() << " ConvBremGsfPFRecTracks" << endl;
02045   } else if ( useConvBremGsfTracks_ ) { 
02046     cerr<<"PFRootEventManager::ProcessEntry : convBremGsfrecTracks Collection not found : "
02047         <<entry << " " << convBremGsfrecTracksTag_<<endl;
02048   }
02049 
02050   bool foundmuons = iEvent.getByLabel(muonsTag_,muonsHandle_);
02051   if ( foundmuons ) { 
02052     muons_ = *muonsHandle_;
02053     /*
02054     cout << "Found " << muons_.size() << " muons" << endl;
02055     for ( unsigned imu=0; imu<muons_.size(); ++imu ) { 
02056       std::cout << " Muon " << imu << ":" << std::endl;
02057       reco::MuonRef muonref( &muons_, imu );
02058       PFMuonAlgo::printMuonProperties(muonref);
02059     }
02060     */
02061   } else { 
02062     cerr<<"PFRootEventManager::ProcessEntry : muons Collection not found : "
02063         <<entry << " " << muonsTag_<<endl;
02064   }
02065 
02066   bool foundconversion = iEvent.getByLabel(conversionTag_,conversionHandle_);
02067   if ( foundconversion ) { 
02068     conversion_ = *conversionHandle_;
02069     // cout << "Found " << conversion_.size() << " conversion" << endl;
02070   } else if ( usePFConversions_ ) { 
02071     cerr<<"PFRootEventManager::ProcessEntry : conversion Collection not found : "
02072         <<entry << " " << conversionTag_<<endl;
02073   }
02074 
02075 
02076 
02077   bool foundv0 = iEvent.getByLabel(v0Tag_,v0Handle_);
02078   if ( foundv0 ) { 
02079     v0_ = *v0Handle_;
02080     // cout << "Found " << v0_.size() << " v0" << endl;
02081   } else if ( usePFV0s_ ) { 
02082     cerr<<"PFRootEventManager::ProcessEntry : v0 Collection not found : "
02083         <<entry << " " << v0Tag_<<endl;
02084   }
02085 
02086   if(useEGPhotons_) {
02087     bool foundPhotons = iEvent.getByLabel(photonTag_,photonHandle_);
02088     if ( foundPhotons) {
02089       photons_ = *photonHandle_;    
02090     } else {
02091       cerr <<"PFRootEventManager::ProcessEntry : photon collection not found : " 
02092            << entry << " " << photonTag_ << endl;
02093     }
02094   }
02095 
02096   bool foundgenJets = iEvent.getByLabel(genJetsTag_,genJetsHandle_);
02097   if ( foundgenJets ) { 
02098     genJetsCMSSW_ = *genJetsHandle_;
02099     // cout << "Found " << genJetsCMSSW_.size() << " genJets" << endl;
02100   } else { 
02101     //cerr<<"PFRootEventManager::ProcessEntry : genJets Collection not found : "
02102     //    <<entry << " " << genJetsTag_<<endl;
02103   }
02104 
02105   bool foundgenParticlesforJets = iEvent.getByLabel(genParticlesforJetsTag_,genParticlesforJetsHandle_);
02106   if ( foundgenParticlesforJets ) { 
02107     genParticlesforJets_ = *genParticlesforJetsHandle_;
02108     // cout << "Found " << genParticlesforJets_.size() << " genParticlesforJets" << endl;
02109   } else { 
02110     //cerr<<"PFRootEventManager::ProcessEntry : genParticlesforJets Collection not found : "
02111     //    <<entry << " " << genParticlesforJetsTag_<<endl;
02112   }
02113 
02114   bool foundgenParticlesforMET = iEvent.getByLabel(genParticlesforMETTag_,genParticlesforMETHandle_);
02115   if ( foundgenParticlesforMET ) { 
02116     genParticlesCMSSW_ = *genParticlesforMETHandle_;
02117     // cout << "Found " << genParticlesCMSSW_.size() << " genParticlesforMET" << endl;
02118   } else { 
02119     //cerr<<"PFRootEventManager::ProcessEntry : genParticlesforMET Collection not found : "
02120     //    <<entry << " " << genParticlesforMETTag_<<endl;
02121   }
02122 
02123   bool foundcaloJets = iEvent.getByLabel(caloJetsTag_,caloJetsHandle_);
02124   if ( foundcaloJets ) { 
02125     caloJetsCMSSW_ = *caloJetsHandle_;
02126     // cout << "Found " << caloJetsCMSSW_.size() << " caloJets" << endl;
02127   } else { 
02128     cerr<<"PFRootEventManager::ProcessEntry : caloJets Collection not found : "
02129         <<entry << " " << caloJetsTag_<<endl;
02130   }
02131 
02132   bool foundcorrcaloJets = iEvent.getByLabel(corrcaloJetsTag_,corrcaloJetsHandle_);
02133   if ( foundcorrcaloJets ) { 
02134     corrcaloJetsCMSSW_ = *corrcaloJetsHandle_;
02135     // cout << "Found " << corrcaloJetsCMSSW_.size() << " corrcaloJets" << endl;
02136   } else { 
02137     //cerr<<"PFRootEventManager::ProcessEntry : corrcaloJets Collection not found : "
02138     //    <<entry << " " << corrcaloJetsTag_<<endl;
02139   }
02140 
02141   bool foundpfJets = iEvent.getByLabel(pfJetsTag_,pfJetsHandle_);
02142   if ( foundpfJets ) { 
02143     pfJetsCMSSW_ = *pfJetsHandle_;
02144     // cout << "Found " << pfJetsCMSSW_.size() << " PFJets" << endl;
02145   } else { 
02146     cerr<<"PFRootEventManager::ProcessEntry : PFJets Collection not found : "
02147         <<entry << " " << pfJetsTag_<<endl;
02148   }
02149 
02150   bool foundpfCands = iEvent.getByLabel(pfCandidateTag_,pfCandidateHandle_);
02151   if ( foundpfCands ) { 
02152     pfCandCMSSW_ = *pfCandidateHandle_;
02153     // cout << "Found " << pfCandCMSSW_.size() << " PFCandidates" << endl;
02154   } else { 
02155     cerr<<"PFRootEventManager::ProcessEntry : PFCandidate Collection not found : "
02156         <<entry << " " << pfCandidateTag_<<endl;
02157   }
02158 
02159   bool foundpfMets = iEvent.getByLabel(pfMetsTag_,pfMetsHandle_);
02160   if ( foundpfMets ) { 
02161     pfMetsCMSSW_ = *pfMetsHandle_;
02162     //cout << "Found " << pfMetsCMSSW_.size() << " PFMets" << endl;
02163   } else { 
02164     cerr<<"PFRootEventManager::ProcessEntry : PFMets Collection not found : "
02165         <<entry << " " << pfMetsTag_<<endl;
02166   }
02167 
02168   bool foundtcMets = iEvent.getByLabel(tcMetsTag_,tcMetsHandle_);
02169   if ( foundtcMets ) { 
02170     tcMetsCMSSW_ = *tcMetsHandle_;
02171     //cout << "Found " << tcMetsCMSSW_.size() << " TCMets" << endl;
02172   } else { 
02173     cerr<<"TCRootEventManager::ProcessEntry : TCMets Collection not found : "
02174         <<entry << " " << tcMetsTag_<<endl;
02175   }
02176 
02177   bool foundcaloMets = iEvent.getByLabel(caloMetsTag_,caloMetsHandle_);
02178   if ( foundcaloMets ) { 
02179     caloMetsCMSSW_ = *caloMetsHandle_;
02180     //cout << "Found " << caloMetsCMSSW_.size() << " CALOMets" << endl;
02181   } else { 
02182     cerr<<"CALORootEventManager::ProcessEntry : CALOMets Collection not found : "
02183         <<entry << " " << caloMetsTag_<<endl;
02184   }
02185 
02186   // now can use the tree
02187 
02188   bool goodevent = true;
02189   if(trueParticles_.size() ) {
02190     // this is a filter to select single particle events.
02191     if(filterNParticles_ && doTauBenchmark_ &&
02192        trueParticles_.size() != filterNParticles_ ) {
02193       cout << "PFRootEventManager : event discarded Nparticles="
02194            <<filterNParticles_<< endl; 
02195       goodevent = false;
02196     }
02197     if(goodevent && doTauBenchmark_ && filterHadronicTaus_ && !isHadronicTau() ) {
02198       cout << "PFRootEventManager : leptonic tau discarded " << endl; 
02199       goodevent =  false;
02200     }
02201     if( goodevent && doTauBenchmark_ && !filterTaus_.empty() 
02202         && !countChargedAndPhotons() ) {
02203       assert( filterTaus_.size() == 2 );
02204       cout <<"PFRootEventManager : tau discarded: "
02205            <<"number of charged and neutral particles different from "
02206            <<filterTaus_[0]<<","<<filterTaus_[1]<<endl;
02207       goodevent =  false;      
02208     } 
02209     
02210     if(goodevent)
02211       fillOutEventWithSimParticles( trueParticles_ );
02212 
02213   }
02214   
02215   //   if(caloTowersBranch_) {
02216   //     if(goodevent)
02217   //       fillOutEventWithCaloTowers( caloTowers_ );
02218   //   } 
02219 
02220   if(rechitsECAL_.size()) {
02221     PreprocessRecHits( rechitsECAL_ , findRecHitNeighbours_);
02222   }
02223   if(rechitsHCAL_.size()) {
02224     PreprocessRecHits( rechitsHCAL_ , findRecHitNeighbours_);
02225   }
02226   if(rechitsHFEM_.size()) {
02227     PreprocessRecHits( rechitsHFEM_ , findRecHitNeighbours_);
02228   }
02229   if(rechitsHFHAD_.size()) {
02230     PreprocessRecHits( rechitsHFHAD_ , findRecHitNeighbours_);
02231   }
02232   rechitsCLEANED_.clear();
02233   for ( unsigned i=0; i<rechitsCLEANEDV_.size(); ++i ) { 
02234     if(rechitsCLEANEDV_[i].size()) {
02235       PreprocessRecHits( rechitsCLEANEDV_[i] , false);
02236       for ( unsigned j=0; j<rechitsCLEANEDV_[i].size(); ++j ) { 
02237         rechitsCLEANED_.push_back( (rechitsCLEANEDV_[i])[j] );
02238       }
02239     }
02240   }
02241 
02242   if(rechitsPS_.size()) {
02243     PreprocessRecHits( rechitsPS_ , findRecHitNeighbours_);
02244   }
02245 
02246   if ( recTracks_.size() ) { 
02247     PreprocessRecTracks( recTracks_);
02248   }
02249 
02250   if ( displacedRecTracks_.size() ) { 
02251     //   cout << "preprocessing rec tracks" << endl;
02252     PreprocessRecTracks( displacedRecTracks_);
02253   }
02254 
02255 
02256   if(gsfrecTracks_.size()) {
02257     PreprocessRecTracks( gsfrecTracks_);
02258   }
02259    
02260   if(convBremGsfrecTracks_.size()) {
02261     PreprocessRecTracks( convBremGsfrecTracks_);
02262   }
02263 
02264   return goodevent;
02265 }
02266 
02267 
02268 bool PFRootEventManager::isHadronicTau() const {
02269 
02270   for ( unsigned i=0;  i < trueParticles_.size(); i++) {
02271     const reco::PFSimParticle& ptc = trueParticles_[i];
02272     const std::vector<int>& ptcdaughters = ptc.daughterIds();
02273     if (std::abs(ptc.pdgCode()) == 15) {
02274       for ( unsigned int dapt=0; dapt < ptcdaughters.size(); ++dapt) {
02275         
02276         const reco::PFSimParticle& daughter 
02277           = trueParticles_[ptcdaughters[dapt]];
02278         
02279 
02280         int pdgdaugther = daughter.pdgCode();
02281         int abspdgdaughter = std::abs(pdgdaugther);
02282 
02283 
02284         if (abspdgdaughter == 11 || 
02285             abspdgdaughter == 13) { 
02286           return false; 
02287         }//electron or muons?
02288       }//loop daughter
02289     }//tau
02290   }//loop particles
02291 
02292 
02293   return true;
02294 }
02295 
02296 
02297 
02298 bool PFRootEventManager::countChargedAndPhotons() const {
02299   
02300   int nPhoton = 0;
02301   int nCharged = 0;
02302   
02303   for ( unsigned i=0;  i < trueParticles_.size(); i++) {
02304     const reco::PFSimParticle& ptc = trueParticles_[i];
02305    
02306     const std::vector<int>& daughters = ptc.daughterIds();
02307 
02308     // if the particle decays before ECAL, we do not want to 
02309     // consider it.
02310     if(!daughters.empty() ) continue; 
02311 
02312     double charge = ptc.charge();
02313     double pdgCode = ptc.pdgCode();
02314     
02315     if( std::abs(charge)>1e-9) 
02316       nCharged++;
02317     else if( pdgCode==22 )
02318       nPhoton++;
02319   }    
02320 
02321   //   const HepMC::GenEvent* myGenEvent = MCTruth_.GetEvent();
02322   //   if(!myGenEvent) {
02323   //     cerr<<"impossible to filter on the number of charged and "
02324   //    <<"neutral particles without the HepMCProduct. "
02325   //    <<"Please check that the branch edmHepMCProduct_*_*_* is found"<<endl;
02326   //     exit(1);
02327   //   }
02328   
02329   //   for ( HepMC::GenEvent::particle_const_iterator 
02330   //      piter  = myGenEvent->particles_begin();
02331   //    piter != myGenEvent->particles_end(); 
02332   //    ++piter ) {
02333     
02334   //     const HepMC::GenParticle* p = *piter;
02335   //     int partId = p->pdg_id();Long64_t lines = T->ReadFile("mon_fichier","i:j:k:x:y:z");
02336     
02337   // //     pdgTable_->GetParticle( partId )->Print();
02338        
02339   //     int charge = chargeValue(partId);
02340   //     cout<<partId <<" "<<charge/3.<<endl;
02341 
02342   //     if(charge) 
02343   //       nCharged++;
02344   //     else 
02345   //       nNeutral++;
02346   //   }
02347   
02348   if( nCharged == filterTaus_[0] && 
02349       nPhoton == filterTaus_[1]  )
02350     return true;
02351   else 
02352     return false;
02353 }
02354 
02355 
02356 
02357 int PFRootEventManager::chargeValue(const int& Id) const {
02358 
02359   
02360   //...Purpose: to give three times the charge for a particle/parton.
02361 
02362   //      ID     = particle ID
02363   //      hepchg = particle charge times 3
02364 
02365   int kqa,kq1,kq2,kq3,kqj,irt,kqx,kqn;
02366   int hepchg;
02367 
02368 
02369   int ichg[109]={-1,2,-1,2,-1,2,-1,2,0,0,-3,0,-3,0,-3,0,
02370                  -3,0,0,0,0,0,0,3,0,0,0,0,0,0,3,0,3,6,0,0,3,6,0,0,-1,2,-1,2,-1,2,0,0,0,0,
02371                  -3,0,-3,0,-3,0,0,0,0,0,-1,2,-1,2,-1,2,0,0,0,0,
02372                  -3,0,-3,0,-3,0,3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
02373 
02374 
02375   //...Initial values. Simple case of direct readout.
02376   hepchg=0;
02377   kqa=std::abs(Id);
02378   kqn=kqa/1000000000%10;
02379   kqx=kqa/1000000%10;
02380   kq3=kqa/1000%10;
02381   kq2=kqa/100%10;
02382   kq1=kqa/10%10;
02383   kqj=kqa%10;
02384   irt=kqa%10000;
02385 
02386   //...illegal or ion
02387   //...set ion charge to zero - not enough information
02388   if(kqa==0 || kqa >= 10000000) {
02389 
02390     if(kqn==1) {hepchg=0;}
02391   }
02392   //... direct translation
02393   else if(kqa<=100) {hepchg = ichg[kqa-1];}
02394   //... deuteron or tritium
02395   else if(kqa==100 || kqa==101) {hepchg = -3;}
02396   //... alpha or He3
02397   else if(kqa==102 || kqa==104) {hepchg = -6;}
02398   //... KS and KL (and undefined)
02399   else if(kqj == 0) {hepchg = 0;}
02400   //C... direct translation
02401   else if(kqx>0 && irt<100)
02402     {
02403       hepchg = ichg[irt-1];
02404       if(kqa==1000017 || kqa==1000018) {hepchg = 0;}
02405       if(kqa==1000034 || kqa==1000052) {hepchg = 0;}
02406       if(kqa==1000053 || kqa==1000054) {hepchg = 0;}
02407       if(kqa==5100061 || kqa==5100062) {hepchg = 6;}
02408     }
02409   //...Construction from quark content for heavy meson, diquark, baryon.
02410   //...Mesons.
02411   else if(kq3==0)
02412     {
02413       hepchg = ichg[kq2-1]-ichg[kq1-1];
02414       //...Strange or beauty mesons.
02415       if((kq2==3) || (kq2==5)) {hepchg = ichg[kq1-1]-ichg[kq2-1];}
02416     }
02417   else if(kq1 == 0) {
02418     //...Diquarks.
02419     hepchg = ichg[kq3-1] + ichg[kq2-1];
02420   }
02421 
02422   else{
02423     //...Baryons
02424     hepchg = ichg[kq3-1]+ichg[kq2-1]+ichg[kq1-1];
02425   }
02426 
02427   //... fix sign of charge
02428   if(Id<0 && hepchg!=0) {hepchg = -1*hepchg;}
02429 
02430   // cout << hepchg<< endl;
02431   return hepchg;
02432 }
02433 
02434 
02435 
02436 void 
02437 PFRootEventManager::PreprocessRecTracks(reco::PFRecTrackCollection& recTracks) {  
02438   for( unsigned i=0; i<recTracks.size(); ++i ) {     
02439     recTracks[i].calculatePositionREP();
02440   }
02441 }
02442 
02443 void 
02444 PFRootEventManager::PreprocessRecTracks(reco::GsfPFRecTrackCollection& recTracks) {  
02445   for( unsigned i=0; i<recTracks.size(); ++i ) {     
02446     recTracks[i].calculatePositionREP();
02447     recTracks[i].calculateBremPositionREP();
02448   }
02449 }
02450 
02451 
02452  
02453 void 
02454 PFRootEventManager::PreprocessRecHits(reco::PFRecHitCollection& rechits, 
02455                                       bool findNeighbours) {
02456   
02457  
02458   map<unsigned, unsigned> detId2index;
02459 
02460   for(unsigned i=0; i<rechits.size(); i++) { 
02461     rechits[i].calculatePositionREP();
02462     
02463     if(findNeighbours) 
02464       detId2index.insert( make_pair(rechits[i].detId(), i) );
02465   }
02466   
02467   if(findNeighbours) {
02468     for(unsigned i=0; i<rechits.size(); i++) { 
02469       setRecHitNeigbours( rechits[i], detId2index );
02470     }
02471   }
02472 }
02473 
02474 
02475 void PFRootEventManager::setRecHitNeigbours
02476 ( reco::PFRecHit& rh, 
02477   const map<unsigned, unsigned>& detId2index ) {
02478 
02479   rh.clearNeighbours();
02480 
02481   vector<unsigned> neighbours4DetId = rh.neighboursIds4();
02482   vector<unsigned> neighbours8DetId = rh.neighboursIds8();
02483   
02484   for( unsigned i=0; i<neighbours4DetId.size(); i++) {
02485     unsigned detId = neighbours4DetId[i];
02486     //     cout<<"finding n for detId "<<detId<<endl;
02487     const map<unsigned, unsigned>::const_iterator& it = detId2index.find(detId);
02488     
02489     if(it != detId2index.end() ) {
02490       //       cout<<"found n index "<<it->second<<endl;
02491       rh.add4Neighbour( it->second );
02492     }
02493   }
02494 
02495   for( unsigned i=0; i<neighbours8DetId.size(); i++) {
02496     unsigned detId = neighbours8DetId[i];
02497     //     cout<<"finding n for detId "<<detId<<endl;
02498     const map<unsigned, unsigned>::const_iterator& it = detId2index.find(detId);
02499     
02500     if(it != detId2index.end() ) {
02501       //       cout<<"found n index "<<it->second<<endl;
02502       rh.add8Neighbour( it->second );
02503     }
02504   }
02505 
02506   
02507 }
02508 
02509 
02510 void PFRootEventManager::clustering() {
02511 
02512   if (verbosity_ == VERBOSE ) {
02513     cout <<"start clustering"<<endl;
02514   }
02515   
02516   vector<bool> mask;
02517   // ECAL clustering -------------------------------------------
02518 
02519   fillRecHitMask( mask, rechitsECAL_ );
02520   clusterAlgoECAL_.doClustering( rechitsECAL_, mask );
02521   clustersECAL_ = clusterAlgoECAL_.clusters();
02522 
02523   assert(clustersECAL_.get() );
02524 
02525   fillOutEventWithClusters( *clustersECAL_ );
02526 
02527   // HCAL clustering -------------------------------------------
02528 
02529   fillRecHitMask( mask, rechitsHCAL_ );
02530   clusterAlgoHCAL_.doClustering( rechitsHCAL_, mask );
02531   clustersHCAL_ = clusterAlgoHCAL_.clusters();
02532 
02533   fillOutEventWithClusters( *clustersHCAL_ );
02534 
02535   // HF clustering -------------------------------------------
02536 
02537   fillRecHitMask( mask, rechitsHFEM_ );
02538   clusterAlgoHFEM_.doClustering( rechitsHFEM_, mask );
02539   clustersHFEM_ = clusterAlgoHFEM_.clusters();
02540   
02541   fillRecHitMask( mask, rechitsHFHAD_ );
02542   clusterAlgoHFHAD_.doClustering( rechitsHFHAD_, mask );
02543   clustersHFHAD_ = clusterAlgoHFHAD_.clusters();
02544   
02545 
02546   // PS clustering -------------------------------------------
02547 
02548   fillRecHitMask( mask, rechitsPS_ );
02549   clusterAlgoPS_.doClustering( rechitsPS_, mask );
02550   clustersPS_ = clusterAlgoPS_.clusters();
02551 
02552   fillOutEventWithClusters( *clustersPS_ );
02553   
02554 }
02555 
02556 
02557 
02558 void 
02559 PFRootEventManager::fillOutEventWithClusters(const reco::PFClusterCollection& 
02560                                              clusters) {
02561 
02562   if(!outEvent_) return;
02563   
02564   for(unsigned i=0; i<clusters.size(); i++) {
02565     EventColin::Cluster cluster;
02566     cluster.eta = clusters[i].position().Eta();
02567     cluster.phi = clusters[i].position().Phi();
02568     cluster.e = clusters[i].energy();
02569     cluster.layer = clusters[i].layer();
02570     cluster.type = 1;
02571 
02572     reco::PFTrajectoryPoint::LayerType tpLayer = 
02573       reco::PFTrajectoryPoint::NLayers;
02574     switch( clusters[i].layer() ) {
02575     case PFLayer::ECAL_BARREL:
02576     case PFLayer::ECAL_ENDCAP:
02577       tpLayer = reco::PFTrajectoryPoint::ECALEntrance;
02578       break;
02579     case PFLayer::HCAL_BARREL1:
02580     case PFLayer::HCAL_ENDCAP:
02581       tpLayer = reco::PFTrajectoryPoint::HCALEntrance;
02582       break;
02583     default:
02584       break;
02585     }
02586     if(tpLayer < reco::PFTrajectoryPoint::NLayers) {
02587       try {
02588         double peta = -10;
02589         double phi = -10;
02590         double pe = -10;
02591         
02592         const reco::PFSimParticle& ptc 
02593           = closestParticle( tpLayer, 
02594                              cluster.eta, cluster.phi, 
02595                              peta, phi, pe );
02596 
02597         
02598         cluster.particle.eta = peta;
02599         cluster.particle.phi = phi;
02600         cluster.particle.e = pe;
02601         cluster.particle.pdgCode = ptc.pdgCode();
02602         
02603         
02604       }
02605       catch( std::exception& err ) {
02606         // cerr<<err.what()<<endl;
02607       } 
02608     }
02609 
02610     outEvent_->addCluster(cluster);
02611   }   
02612 }
02613 
02614 
02615 void 
02616 PFRootEventManager::fillOutEventWithSimParticles(const reco::PFSimParticleCollection& trueParticles ) {
02617 
02618   if(!outEvent_) return;
02619   
02620   for ( unsigned i=0;  i < trueParticles.size(); i++) {
02621 
02622     const reco::PFSimParticle& ptc = trueParticles[i];
02623 
02624     unsigned ntrajpoints = ptc.nTrajectoryPoints();
02625     
02626     if(ptc.daughterIds().empty() ) { // stable
02627       reco::PFTrajectoryPoint::LayerType ecalEntrance 
02628         = reco::PFTrajectoryPoint::ECALEntrance;
02629 
02630       if(ntrajpoints == 3) { 
02631         // old format for PFSimCandidates. 
02632         // in this case, the PFSimCandidate which does not decay 
02633         // before ECAL has 3 points: initial, ecal entrance, hcal entrance
02634         ecalEntrance = static_cast<reco::PFTrajectoryPoint::LayerType>(1);
02635       }
02636       // else continue; // endcap case we do not care;
02637 
02638       const reco::PFTrajectoryPoint& tpatecal 
02639         = ptc.extrapolatedPoint( ecalEntrance );
02640         
02641       EventColin::Particle outptc;
02642       outptc.eta = tpatecal.position().Eta();
02643       outptc.phi = tpatecal.position().Phi();    
02644       outptc.e = tpatecal.momentum().E();
02645       outptc.pdgCode = ptc.pdgCode();
02646     
02647       
02648       outEvent_->addParticle(outptc);
02649     }  
02650   }   
02651 }      
02652 
02653 void 
02654 PFRootEventManager::fillOutEventWithPFCandidates(const reco::PFCandidateCollection& pfCandidates ) {
02655 
02656   if(!outEvent_) return;
02657   
02658   for ( unsigned i=0;  i < pfCandidates.size(); i++) {
02659 
02660     const reco::PFCandidate& candidate = pfCandidates[i];
02661     
02662     EventColin::Particle outptc;
02663     outptc.eta = candidate.eta();
02664     outptc.phi = candidate.phi();    
02665     outptc.e = candidate.energy();
02666     outptc.pdgCode = candidate.particleId();
02667     
02668     
02669     outEvent_->addCandidate(outptc);  
02670   }   
02671 }      
02672 
02673 
02674 void 
02675 PFRootEventManager::fillOutEventWithCaloTowers( const CaloTowerCollection& cts ){
02676 
02677   if(!outEvent_) return;
02678   
02679   for ( unsigned i=0;  i < cts.size(); i++) {
02680 
02681     const CaloTower& ct = cts[i];
02682     
02683     EventColin::CaloTower outct;
02684     outct.e  = ct.energy();
02685     outct.ee = ct.emEnergy();
02686     outct.eh = ct.hadEnergy();
02687 
02688     outEvent_->addCaloTower( outct );
02689   }
02690 }
02691 
02692 
02693 void 
02694 PFRootEventManager::fillOutEventWithBlocks( const reco::PFBlockCollection& 
02695                                             blocks ) {
02696 
02697   if(!outEvent_) return;
02698   
02699   for ( unsigned i=0;  i < blocks.size(); i++) {
02700 
02701     //    const reco::PFBlock& block = blocks[i];
02702     
02703     EventColin::Block outblock;
02704  
02705     outEvent_->addBlock( outblock );
02706   }
02707 }
02708 
02709 
02710 
02711 void PFRootEventManager::particleFlow() {
02712   
02713   if (verbosity_ == VERBOSE ) {
02714     cout <<"start particle flow"<<endl;
02715   }
02716 
02717 
02718   if( debug_) {
02719     cout<<"PFRootEventManager::particleFlow start"<<endl;
02720     //     cout<<"number of elements in memory: "
02721     //  <<reco::PFBlockElement::instanceCounter()<<endl;
02722   }
02723 
02724 
02725   edm::OrphanHandle< reco::PFRecTrackCollection > trackh( &recTracks_, 
02726                                                           edm::ProductID(1) );  
02727   
02728   edm::OrphanHandle< reco::PFRecTrackCollection > displacedtrackh( &displacedRecTracks_, 
02729                                                           edm::ProductID(77) );  
02730 
02731   edm::OrphanHandle< reco::PFClusterCollection > ecalh( clustersECAL_.get(), 
02732                                                         edm::ProductID(2) );  
02733   
02734   edm::OrphanHandle< reco::PFClusterCollection > hcalh( clustersHCAL_.get(), 
02735                                                         edm::ProductID(3) );  
02736 
02737   edm::OrphanHandle< reco::PFClusterCollection > hfemh( clustersHFEM_.get(), 
02738                                                         edm::ProductID(31) );  
02739 
02740   edm::OrphanHandle< reco::PFClusterCollection > hfhadh( clustersHFHAD_.get(), 
02741                                                         edm::ProductID(32) );  
02742 
02743   edm::OrphanHandle< reco::PFClusterCollection > psh( clustersPS_.get(), 
02744                                                       edm::ProductID(4) );   
02745   
02746   edm::OrphanHandle< reco::GsfPFRecTrackCollection > gsftrackh( &gsfrecTracks_, 
02747                                                                 edm::ProductID(5) );  
02748   
02749   edm::OrphanHandle< reco::MuonCollection > muonh( &muons_, 
02750                                                    edm::ProductID(6) );
02751 
02752   edm::OrphanHandle< reco::PFDisplacedTrackerVertexCollection > nuclearh( &pfNuclearTrackerVertex_, 
02753                                                           edm::ProductID(7) );
02754 
02755 
02756   //recoPFRecTracks_pfNuclearTrackerVertex__TEST.
02757 
02758   edm::OrphanHandle< reco::PFConversionCollection > convh( &conversion_, 
02759                                                            edm::ProductID(8) );
02760 
02761   edm::OrphanHandle< reco::PFV0Collection > v0( &v0_, 
02762                                                 edm::ProductID(9) );
02763 
02764 
02765   edm::OrphanHandle< reco::VertexCollection > vertexh( &primaryVertices_, 
02766                                                        edm::ProductID(10) );  
02767 
02768   edm::OrphanHandle< reco::GsfPFRecTrackCollection > convBremGsftrackh( &convBremGsfrecTracks_, 
02769                                                                         edm::ProductID(11) );  
02770 
02771   edm::OrphanHandle< reco::PhotonCollection > photonh( &photons_, edm::ProductID(12) ) ;
02772   
02773   vector<bool> trackMask;
02774   fillTrackMask( trackMask, recTracks_ );
02775   vector<bool> gsftrackMask;
02776   fillTrackMask( gsftrackMask, gsfrecTracks_ );
02777   vector<bool> ecalMask;
02778   fillClusterMask( ecalMask, *clustersECAL_ );
02779   vector<bool> hcalMask;
02780   fillClusterMask( hcalMask, *clustersHCAL_ );
02781   vector<bool> hfemMask;
02782   fillClusterMask( hfemMask, *clustersHFEM_ );
02783   vector<bool> hfhadMask;
02784   fillClusterMask( hfhadMask, *clustersHFHAD_ );
02785   vector<bool> psMask;
02786   fillClusterMask( psMask, *clustersPS_ );
02787   vector<bool> photonMask;
02788   fillPhotonMask( photonMask, photons_ );
02789 
02790   if ( !useAtHLT )
02791     pfBlockAlgo_.setInput( trackh, gsftrackh, convBremGsftrackh,
02792                            muonh, nuclearh, displacedtrackh, convh, v0,
02793                            ecalh, hcalh, hfemh, hfhadh, psh, photonh,
02794                            trackMask,gsftrackMask,
02795                            ecalMask, hcalMask, hfemMask, hfhadMask, psMask,photonMask );
02796   else    
02797     pfBlockAlgo_.setInput( trackh, muonh, ecalh, hcalh, hfemh, hfhadh, psh,
02798                            trackMask, ecalMask, hcalMask, psMask );
02799 
02800   pfBlockAlgo_.findBlocks();
02801   
02802   if( debug_) cout<<pfBlockAlgo_<<endl;
02803 
02804   pfBlocks_ = pfBlockAlgo_.transferBlocks();
02805 
02806   pfAlgo_.setPFVertexParameters(true, primaryVertices_); 
02807 
02808   pfAlgo_.reconstructParticles( *pfBlocks_.get() );
02809   //   pfAlgoOther_.reconstructParticles( blockh );
02810 
02811   pfAlgo_.postMuonCleaning(muonsHandle_, *vertexh);
02812 
02813   pfAlgo_.checkCleaning( rechitsCLEANED_ );
02814 
02815   if( debug_) cout<< pfAlgo_<<endl;
02816   pfCandidates_ = pfAlgo_.transferCandidates();
02817   //   pfCandidatesOther_ = pfAlgoOther_.transferCandidates();
02818   
02819   fillOutEventWithPFCandidates( *pfCandidates_ );
02820 
02821   if( debug_) cout<<"PFRootEventManager::particleFlow stop"<<endl;
02822 }
02823 
02824 void PFRootEventManager::pfCandCompare(int entry) {
02825 
02826   bool differentSize = pfCandCMSSW_.size() != pfCandidates_->size();
02827   if ( differentSize ) { 
02828     cout << "+++WARNING+++ PFCandidate size changed for entry " 
02829          << entry << " !" << endl
02830          << " - original size : " << pfCandCMSSW_.size() << endl 
02831          << " - current  size : " << pfCandidates_->size() << endl;
02832   } else { 
02833     for(unsigned i=0; i<pfCandidates_->size(); i++) {
02834       double deltaE = (*pfCandidates_)[i].energy()-pfCandCMSSW_[i].energy();
02835       double deltaEta = (*pfCandidates_)[i].eta()-pfCandCMSSW_[i].eta();
02836       double deltaPhi = (*pfCandidates_)[i].phi()-pfCandCMSSW_[i].phi();
02837       if ( fabs(deltaE) > 1E-5 ||
02838            fabs(deltaEta) > 1E-9 ||
02839            fabs(deltaPhi) > 1E-9 ) { 
02840         cout << "+++WARNING+++ PFCandidate " << i 
02841              << " changed  for entry " << entry << " ! " << endl 
02842              << " - Original : " << pfCandCMSSW_[i] << endl
02843              << " - Current  : " << (*pfCandidates_)[i] << endl
02844              << " DeltaE   = : " << deltaE << endl
02845              << " DeltaEta = : " << deltaEta << endl
02846              << " DeltaPhi = : " << deltaPhi << endl << endl;
02847       }
02848     }
02849   }
02850 
02851 }
02852 
02853 
02854 void PFRootEventManager::reconstructGenJets() {
02855 
02856   if (verbosity_ == VERBOSE || jetsDebug_) {
02857     cout<<endl;
02858     cout<<"start reconstruct GenJets  --- "<<endl;
02859     cout<< " input gen particles for jet: all neutrinos removed ; muons present" << endl;
02860   }
02861 
02862   genJets_.clear();
02863   genParticlesforJetsPtrs_.clear();
02864 
02865   if ( !genParticlesforJets_.size() ) return;
02866 
02867   for(unsigned i=0; i<genParticlesforJets_.size(); i++) {
02868 
02869     const reco::GenParticle&    genPart = *(genParticlesforJets_[i]);
02870 
02871     // remove all muons/neutrinos for PFJet studies
02872     //    if (reco::isNeutrino( genPart ) || reco::isMuon( genPart )) continue;
02873     //    remove all neutrinos for PFJet studies
02874     if (reco::isNeutrino( genPart )) continue;
02875     // Work-around a bug in the pythia di-jet gun.
02876     if (std::abs(genPart.pdgId())<7 || std::abs(genPart.pdgId())==21 ) continue;
02877 
02878     if (jetsDebug_ ) {
02879       cout << "      #" << i << "  PDG code:" << genPart.pdgId() 
02880            << " status " << genPart.status()
02881            << ", p/pt/eta/phi: " << genPart.p() << '/' << genPart.pt() 
02882            << '/' << genPart.eta() << '/' << genPart.phi() << endl;
02883     }
02884     
02885     genParticlesforJetsPtrs_.push_back( refToPtr(genParticlesforJets_[i]) );
02886   }
02887   
02888   vector<ProtoJet> protoJets;
02889   reconstructFWLiteJets(genParticlesforJetsPtrs_, protoJets );
02890 
02891 
02892   // Convert Protojets to GenJets
02893   int ijet = 0;
02894   typedef vector <ProtoJet>::const_iterator IPJ;
02895   for  (IPJ ipj = protoJets.begin(); ipj != protoJets.end (); ipj++) {
02896     const ProtoJet& protojet = *ipj;
02897     const ProtoJet::Constituents& constituents = protojet.getTowerList();
02898           
02899     reco::Jet::Point vertex (0,0,0); // do not have true vertex yet, use default
02900     GenJet::Specific specific;
02901     JetMaker::makeSpecific(constituents, &specific);
02902     // constructor without constituents
02903     GenJet newJet (protojet.p4(), vertex, specific);
02904           
02905     // last step is to copy the constituents into the jet (new jet definition since 18X)
02906     // namespace reco {
02907     //class Jet : public CompositeRefBaseCandidate {
02908     // public:
02909     //  typedef reco::CandidateBaseRefVector Constituents;
02910           
02911     ProtoJet::Constituents::const_iterator constituent = constituents.begin();
02912     for (; constituent != constituents.end(); ++constituent) {
02913       // find index of this ProtoJet constituent in the overall collection PFconstit
02914       // see class IndexedCandidate in JetRecoTypes.h
02915       uint index = constituent->index();
02916       newJet.addDaughter( genParticlesforJetsPtrs_[index] );
02917     }  // end loop on ProtoJet constituents
02918     // last step: copy ProtoJet Variables into Jet
02919     newJet.setJetArea(protojet.jetArea()); 
02920     newJet.setPileup(protojet.pileup());
02921     newJet.setNPasses(protojet.nPasses());
02922     ++ijet;
02923     if (jetsDebug_ ) cout<<" gen jet "<<ijet<<" "<<newJet.print()<<endl;
02924     genJets_.push_back (newJet);
02925           
02926   } // end loop on protojets iterator IPJ
02927   
02928 }
02929 
02930 void PFRootEventManager::reconstructCaloJets() {
02931 
02932   if (verbosity_ == VERBOSE || jetsDebug_ ) {
02933     cout<<endl;
02934     cout<<"start reconstruct CaloJets --- "<<endl;
02935   }
02936   caloJets_.clear();
02937   caloTowersPtrs_.clear();
02938 
02939   for( unsigned i=0; i<caloTowers_.size(); i++) {
02940     reco::CandidatePtr candPtr( &caloTowers_, i );
02941     caloTowersPtrs_.push_back( candPtr );
02942   }
02943  
02944   reconstructFWLiteJets( caloTowersPtrs_, caloJets_ );
02945 
02946   if (jetsDebug_ ) {
02947     for(unsigned ipj=0; ipj<caloJets_.size(); ipj++) {
02948       const ProtoJet& protojet = caloJets_[ipj];      
02949       cout<<" calo jet "<<ipj<<" "<<protojet.pt() <<endl;
02950     }
02951   }
02952 
02953 }
02954 
02955 
02956 void PFRootEventManager::reconstructPFJets() {
02957 
02958   if (verbosity_ == VERBOSE || jetsDebug_) {
02959     cout<<endl;
02960     cout<<"start reconstruct PF Jets --- "<<endl;
02961   }
02962   pfJets_.clear();
02963   pfCandidatesPtrs_.clear();
02964         
02965   for( unsigned i=0; i<pfCandidates_->size(); i++) {
02966     reco::CandidatePtr candPtr( pfCandidates_.get(), i );
02967     pfCandidatesPtrs_.push_back( candPtr );
02968   }
02969 
02970   vector<ProtoJet> protoJets;
02971   reconstructFWLiteJets(pfCandidatesPtrs_, protoJets );
02972 
02973   // Convert Protojets to PFJets
02974 
02975   int ijet = 0;
02976   typedef vector <ProtoJet>::const_iterator IPJ;
02977   for  (IPJ ipj = protoJets.begin(); ipj != protoJets.end (); ipj++) {
02978     const ProtoJet& protojet = *ipj;
02979     const ProtoJet::Constituents& constituents = protojet.getTowerList();
02980         
02981     reco::Jet::Point vertex (0,0,0); // do not have true vertex yet, use default
02982     PFJet::Specific specific;
02983     JetMaker::makeSpecific(constituents, &specific);
02984     // constructor without constituents
02985     PFJet newJet (protojet.p4(), vertex, specific);
02986         
02987     // last step is to copy the constituents into the jet (new jet definition since 18X)
02988     // namespace reco {
02989     //class Jet : public CompositeRefBaseCandidate {
02990     // public:
02991     //  typedef reco::CandidateBaseRefVector Constituents;
02992         
02993     ProtoJet::Constituents::const_iterator constituent = constituents.begin();
02994     for (; constituent != constituents.end(); ++constituent) {
02995       // find index of this ProtoJet constituent in the overall collection PFconstit
02996       // see class IndexedCandidate in JetRecoTypes.h
02997       uint index = constituent->index();
02998       newJet.addDaughter(pfCandidatesPtrs_[index]);
02999     }  // end loop on ProtoJet constituents
03000     // last step: copy ProtoJet Variables into Jet
03001     newJet.setJetArea(protojet.jetArea()); 
03002     newJet.setPileup(protojet.pileup());
03003     newJet.setNPasses(protojet.nPasses());
03004     ++ijet;
03005     if (jetsDebug_ )  cout<<" PF jet "<<ijet<<" "<<newJet.print()<<endl;
03006     pfJets_.push_back (newJet);
03007         
03008   } // end loop on protojets iterator IPJ
03009 
03010 }
03011 
03012 void 
03013 PFRootEventManager::reconstructFWLiteJets(const reco::CandidatePtrVector& Candidates, vector<ProtoJet>& output ) {
03014 
03015   // cout<<"!!! Make FWLite Jets  "<<endl;  
03016   JetReco::InputCollection input;
03017   // vector<ProtoJet> output;
03018   jetMaker_.applyCuts (Candidates, &input); 
03019   if (jetAlgoType_==1) {// ICone 
03021     jetMaker_.makeIterativeConeJets(input, &output);
03022   }
03023   if (jetAlgoType_==2) {// MCone
03024     jetMaker_.makeMidpointJets(input, &output);
03025   }     
03026   if (jetAlgoType_==3) {// Fastjet
03027     jetMaker_.makeFastJets(input, &output);  
03028   }
03029   if((jetAlgoType_>3)||(jetAlgoType_<0)) {
03030     cout<<"Unknown Jet Algo ! " <<jetAlgoType_ << endl;
03031   }
03032   //if (jetsDebug_) cout<<"Proto Jet Size " <<output.size()<<endl;
03033 
03034 }
03035 
03036 
03037 
03039 double 
03040 PFRootEventManager::tauBenchmark( const reco::PFCandidateCollection& candidates) {
03041   //std::cout << "building jets from MC particles," 
03042   //    << "PF particles and caloTowers" << std::endl;
03043   
03044   //initialize Jets Reconstruction
03045   jetAlgo_.Clear();
03046 
03047   //COLIN The following comment is not really adequate, 
03048   // since partTOTMC is not an action..
03049   // one should say what this variable is for.
03050   // see my comment later 
03051   //MAKING TRUE PARTICLE JETS
03052 //   TLorentzVector partTOTMC;
03053 
03054   // colin: the following is not necessary
03055   // since the lorentz vectors are initialized to 0,0,0,0. 
03056   // partTOTMC.SetPxPyPzE(0.0, 0.0, 0.0, 0.0);
03057 
03058   //MAKING JETS WITH TAU DAUGHTERS
03059   //Colin: this vector vectPART is not necessary !!
03060   //it was just an efficient copy of trueparticles_.....
03061 //   vector<reco::PFSimParticle> vectPART;
03062 //   for ( unsigned i=0;  i < trueParticles_.size(); i++) {
03063 //     const reco::PFSimParticle& ptc = trueParticles_[i];
03064 //     vectPART.push_back(ptc);
03065 //   }//loop
03066 
03067 
03068   //COLIN one must not loop on trueParticles_ to find taus. 
03069   //the code was giving wrong results on non single tau events. 
03070 
03071   // first check that this is a single tau event. 
03072 
03073   TLorentzVector partTOTMC;
03074   bool tauFound = false;
03075   bool tooManyTaus = false;
03076   if (fastsim_){
03077 
03078     for ( unsigned i=0;  i < trueParticles_.size(); i++) {
03079       const reco::PFSimParticle& ptc = trueParticles_[i];
03080       if (std::abs(ptc.pdgCode()) == 15) {
03081         // this is a tau
03082         if( i ) tooManyTaus = true;
03083         else tauFound=true;
03084       }
03085     }
03086     
03087     if(!tauFound || tooManyTaus ) {
03088       // cerr<<"PFRootEventManager::tauBenchmark : not a single tau event"<<endl;
03089       return -9999;
03090     }
03091     
03092     // loop on the daugthers of the tau
03093     const std::vector<int>& ptcdaughters = trueParticles_[0].daughterIds();
03094     
03095     // will contain the sum of the lorentz vectors of the visible daughters
03096     // of the tau.
03097     
03098     
03099     for ( unsigned int dapt=0; dapt < ptcdaughters.size(); ++dapt) {
03100       
03101       const reco::PFTrajectoryPoint& tpatvtx 
03102         = trueParticles_[ptcdaughters[dapt]].trajectoryPoint(0);
03103       TLorentzVector partMC;
03104       partMC.SetPxPyPzE(tpatvtx.momentum().Px(),
03105                         tpatvtx.momentum().Py(),
03106                         tpatvtx.momentum().Pz(),
03107                         tpatvtx.momentum().E());
03108       
03109       partTOTMC += partMC;
03110       if (tauBenchmarkDebug_) {
03111         //pdgcode
03112         int pdgcode =  trueParticles_[ptcdaughters[dapt]].pdgCode();
03113         cout << pdgcode << endl;
03114         cout << tpatvtx << endl;
03115         cout << partMC.Px() << " " << partMC.Py() << " " 
03116              << partMC.Pz() << " " << partMC.E()
03117              << " PT=" 
03118              << sqrt(partMC.Px()*partMC.Px()+partMC.Py()*partMC.Py()) 
03119              << endl;
03120       }//debug
03121     }//loop daughter
03122   }else{
03123 
03124     uint itau=0;
03125     const HepMC::GenEvent* myGenEvent = MCTruth_.GetEvent();
03126     for ( HepMC::GenEvent::particle_const_iterator 
03127             piter  = myGenEvent->particles_begin();
03128           piter != myGenEvent->particles_end(); 
03129           ++piter ) {
03130       
03131     
03132       if (std::abs((*piter)->pdg_id())==15){
03133         itau++;
03134         tauFound=true;
03135         for ( HepMC::GenVertex::particles_out_const_iterator bp =
03136                 (*piter)->end_vertex()->particles_out_const_begin();
03137               bp != (*piter)->end_vertex()->particles_out_const_end(); ++bp ) {
03138           uint nuId=std::abs((*bp)->pdg_id());
03139           bool isNeutrino=(nuId==12)||(nuId==14)||(nuId==16);
03140           if (!isNeutrino){
03141             
03142 
03143             TLorentzVector partMC;
03144             partMC.SetPxPyPzE((*bp)->momentum().x(),
03145                               (*bp)->momentum().y(),
03146                               (*bp)->momentum().z(),
03147                               (*bp)->momentum().e());
03148             partTOTMC += partMC;
03149           }
03150         }
03151       }
03152     }
03153     if (itau>1) tooManyTaus=true;
03154 
03155     if(!tauFound || tooManyTaus ) {
03156       cerr<<"PFRootEventManager::tauBenchmark : not a single tau event"<<endl;
03157       return -9999;
03158     }
03159   }
03160 
03161 
03162 
03163 
03164 
03165 
03166 
03167   EventColin::Jet jetmc;
03168 
03169   jetmc.eta = partTOTMC.Eta();
03170   jetmc.phi = partTOTMC.Phi();
03171   jetmc.et = partTOTMC.Et();
03172   jetmc.e = partTOTMC.E();
03173   
03174   if(outEvent_) outEvent_->addJetMC( jetmc );
03175 
03176   /*
03177   //MC JETS RECONSTRUCTION (visible energy)
03178   for ( unsigned i=0;  i < trueParticles_.size(); i++) {
03179   const reco::PFSimParticle& ptc = trueParticles_[i];
03180   const std::vector<int>& ptcdaughters = ptc.daughterIds();
03181     
03182   //PARTICULE NOT DISINTEGRATING BEFORE ECAL
03183   if(ptcdaughters.size() != 0) continue;
03184     
03185   //TAKE INFO AT VERTEX //////////////////////////////////////////////////
03186   const reco::PFTrajectoryPoint& tpatvtx = ptc.trajectoryPoint(0);
03187   TLorentzVector partMC;
03188   partMC.SetPxPyPzE(tpatvtx.momentum().Px(),
03189   tpatvtx.momentum().Py(),
03190   tpatvtx.momentum().Pz(),
03191   tpatvtx.momentum().E());
03192     
03193   partTOTMC += partMC;
03194   if (tauBenchmarkDebug_) {
03195   //pdgcode
03196   int pdgcode = ptc.pdgCode();
03197   cout << pdgcode << endl;
03198   cout << tpatvtx << endl;
03199   cout << partMC.Px() << " " << partMC.Py() << " " 
03200   << partMC.Pz() << " " << partMC.E() 
03201   << " PT=" 
03202   << sqrt(partMC.Px()*partMC.Px()+partMC.Py()*partMC.Py()) 
03203   << endl;
03204   }//debug?
03205   }//loop true particles
03206   */
03207   if (tauBenchmarkDebug_) {
03208     cout << " ET Vector=" << partTOTMC.Et() 
03209          << " " << partTOTMC.Eta() 
03210          << " " << partTOTMC.Phi() << endl; cout << endl;
03211   }//debug
03212 
03214   //CALO TOWER JETS (ECAL+HCAL Towers)
03215   //cout << endl;  
03216   //cout << "THERE ARE " << caloTowers_.size() << " CALO TOWERS" << endl;
03217 
03218   vector<TLorentzVector> allcalotowers;
03219   //   vector<double>         allemenergy;
03220   //   vector<double>         allhadenergy;
03221   double threshCaloTowers = 1E-10;
03222   for ( unsigned int i = 0; i < caloTowers_.size(); ++i) {
03223     
03224     if(caloTowers_[i].energy() < threshCaloTowers) {
03225       //     cout<<"skipping calotower"<<endl;
03226       continue;
03227     }
03228 
03229     TLorentzVector caloT;
03230     TVector3 pepr( caloTowers_[i].eta(),
03231                    caloTowers_[i].phi(),
03232                    caloTowers_[i].energy());
03233     TVector3 pxyz = Utils::VectorEPRtoXYZ( pepr );
03234     caloT.SetPxPyPzE(pxyz.X(),pxyz.Y(),pxyz.Z(),caloTowers_[i].energy());
03235     allcalotowers.push_back(caloT);
03236     //     allemenergy.push_back( caloTowers_[i].emEnergy() );
03237     //     allhadenergy.push_back( caloTowers_[i].hadEnergy() );
03238   }//loop calo towers
03239   if ( tauBenchmarkDebug_)  
03240     cout << " RETRIEVED " << allcalotowers.size() 
03241          << " CALOTOWER 4-VECTORS " << endl;
03242   
03243   //ECAL+HCAL tower jets computation
03244   jetAlgo_.Clear();
03245   const vector< PFJetAlgorithm::Jet >&  caloTjets 
03246     = jetAlgo_.FindJets( &allcalotowers );
03247   
03248   //cout << caloTjets.size() << " CaloTower Jets found" << endl;
03249   double JetEHTETmax = 0.0;
03250   for ( unsigned i = 0; i < caloTjets.size(); i++) {
03251     TLorentzVector jetmom = caloTjets[i].GetMomentum();
03252     double jetcalo_pt = sqrt(jetmom.Px()*jetmom.Px()+jetmom.Py()*jetmom.Py());
03253     double jetcalo_et = jetmom.Et();
03254 
03255 
03256 
03257     EventColin::Jet jet;
03258     jet.eta = jetmom.Eta();
03259     jet.phi = jetmom.Phi();
03260     jet.et  = jetmom.Et();
03261     jet.e   = jetmom.E();
03262     
03263     const vector<int>& indexes = caloTjets[i].GetIndexes();
03264     for( unsigned ii=0; ii<indexes.size(); ii++){
03265       jet.ee   +=  caloTowers_[ indexes[ii] ].emEnergy();
03266       jet.eh   +=  caloTowers_[ indexes[ii] ].hadEnergy();
03267       jet.ete   +=  caloTowers_[ indexes[ii] ].emEt();
03268       jet.eth   +=  caloTowers_[ indexes[ii] ].hadEt();
03269     }
03270     
03271     if(outEvent_) outEvent_->addJetEHT( jet );
03272 
03273     if ( tauBenchmarkDebug_) {
03274       cout << " ECAL+HCAL jet : " << caloTjets[i] << endl;
03275       cout << jetmom.Px() << " " << jetmom.Py() << " " 
03276            << jetmom.Pz() << " " << jetmom.E() 
03277            << " PT=" << jetcalo_pt << endl;
03278     }//debug
03279 
03280     if (jetcalo_et >= JetEHTETmax) 
03281       JetEHTETmax = jetcalo_et;
03282   }//loop MCjets
03283 
03285   //PARTICLE FLOW JETS
03286   vector<TLorentzVector> allrecparticles;
03287   //   if ( tauBenchmarkDebug_) {
03288   //     cout << endl;
03289   //     cout << " THERE ARE " << pfBlocks_.size() << " EFLOW BLOCKS" << endl;
03290   //   }//debug
03291 
03292   //   for ( unsigned iefb = 0; iefb < pfBlocks_.size(); iefb++) {
03293   //       const std::vector< PFBlockParticle >& recparticles 
03294   //    = pfBlocks_[iefb].particles();
03295 
03296   
03297   
03298   for(unsigned i=0; i<candidates.size(); i++) {
03299   
03300     //       if (tauBenchmarkDebug_) 
03301     //  cout << " there are " << recparticles.size() 
03302     //       << " particle in this block" << endl;
03303     
03304     const reco::PFCandidate& candidate = candidates[i];
03305 
03306     if (tauBenchmarkDebug_) {
03307       cout << i << " " << candidate << endl;
03308       int type = candidate.particleId();
03309       cout << " type= " << type << " " << candidate.charge() 
03310            << endl;
03311     }//debug
03312 
03313     const math::XYZTLorentzVector& PFpart = candidate.p4();
03314     
03315     TLorentzVector partRec(PFpart.Px(), 
03316                            PFpart.Py(), 
03317                            PFpart.Pz(),
03318                            PFpart.E());
03319     
03320     //loading 4-vectors of Rec particles
03321     allrecparticles.push_back( partRec );
03322 
03323   }//loop on candidates
03324   
03325 
03326   if (tauBenchmarkDebug_) 
03327     cout << " THERE ARE " << allrecparticles.size() 
03328          << " RECONSTRUCTED 4-VECTORS" << endl;
03329 
03330   jetAlgo_.Clear();
03331   const vector< PFJetAlgorithm::Jet >&  PFjets 
03332     = jetAlgo_.FindJets( &allrecparticles );
03333 
03334   if (tauBenchmarkDebug_) 
03335     cout << PFjets.size() << " PF Jets found" << endl;
03336   double JetPFETmax = 0.0;
03337   for ( unsigned i = 0; i < PFjets.size(); i++) {
03338     TLorentzVector jetmom = PFjets[i].GetMomentum();
03339     double jetpf_pt = sqrt(jetmom.Px()*jetmom.Px()+jetmom.Py()*jetmom.Py());
03340     double jetpf_et = jetmom.Et();
03341 
03342     EventColin::Jet jet;
03343     jet.eta = jetmom.Eta();
03344     jet.phi = jetmom.Phi();
03345     jet.et  = jetmom.Et();
03346     jet.e   = jetmom.E();
03347 
03348     if(outEvent_) outEvent_->addJetPF( jet );
03349 
03350     if (tauBenchmarkDebug_) {
03351       cout <<" Rec jet : "<< PFjets[i] <<endl;
03352       cout << jetmom.Px() << " " << jetmom.Py() << " " 
03353            << jetmom.Pz() << " " << jetmom.E() 
03354            << " PT=" << jetpf_pt << " eta="<< jetmom.Eta() 
03355            << " Phi=" << jetmom.Phi() << endl;
03356       cout << "------------------------------------------------" << endl;
03357     }//debug
03358     
03359     if (jetpf_et >= JetPFETmax)  
03360       JetPFETmax = jetpf_et;
03361   }//loop PF jets
03362 
03363   //fill histos
03364 
03365   double deltaEtEHT = JetEHTETmax - partTOTMC.Et();
03366   h_deltaETvisible_MCEHT_->Fill(deltaEtEHT);
03367   
03368   double deltaEt = JetPFETmax - partTOTMC.Et();
03369   h_deltaETvisible_MCPF_ ->Fill(deltaEt);
03370 
03371   if (verbosity_ == VERBOSE ) {
03372     cout << "tau benchmark E_T(PF) - E_T(true) = " << deltaEt << endl;
03373   }
03374 
03375   return deltaEt/partTOTMC.Et();
03376 }//Makejets
03377 
03378 
03379 
03380 
03381 
03382 /*
03383 
03384 void PFRootEventManager::lookForGenParticle(unsigned barcode) {
03385   
03386 const HepMC::GenEvent* event = MCTruth_.GetEvent();
03387 if(!event) {
03388 cerr<<"no GenEvent"<<endl;
03389 return;
03390 }
03391   
03392 const HepMC::GenParticle* particle = event->barcode_to_particle(barcode);
03393 if(!particle) {
03394 cerr<<"no particle with barcode "<<barcode<<endl;
03395 return;
03396 }
03397 
03398 math::XYZTLorentzVector momentum(particle->momentum().px(),
03399 particle->momentum().py(),
03400 particle->momentum().pz(),
03401 particle->momentum().e());
03402 
03403 double eta = momentum.Eta();
03404 double phi = momentum.phi();
03405 
03406 double phisize = 0.05;
03407 double etasize = 0.05;
03408   
03409 double etagate = displayZoomFactor_ * etasize;
03410 double phigate = displayZoomFactor_ * phisize;
03411   
03412 if(displayHist_[EPE]) {
03413 displayHist_[EPE]->GetXaxis()->SetRangeUser(eta-etagate, eta+etagate);
03414 displayHist_[EPE]->GetYaxis()->SetRangeUser(phi-phigate, phi+phigate);
03415 }
03416 if(displayHist_[EPH]) {
03417 displayHist_[EPH]->GetXaxis()->SetRangeUser(eta-etagate, eta+etagate);
03418 displayHist_[EPH]->GetYaxis()->SetRangeUser(phi-phigate, phi+phigate);
03419 }
03420   
03421 updateDisplay();
03422 
03423 }
03424 */
03425 
03426 
03427 
03428 string PFRootEventManager::expand(const string& oldString) const {
03429 
03430   string newString = oldString;
03431  
03432   string dollar = "$";
03433   string slash  = "/";
03434   
03435   // protection necessary or segv !!
03436   int dollarPos = newString.find(dollar,0);
03437   if( dollarPos == -1 ) return oldString;
03438 
03439   int    lengh  = newString.find(slash,0) - newString.find(dollar,0) + 1;
03440   string env_variable =
03441     newString.substr( ( newString.find(dollar,0) + 1 ), lengh -2);
03442   // the env var could be defined between { }
03443   int begin = env_variable.find_first_of("{");
03444   int end = env_variable.find_last_of("}");
03445   
03446   // cout << "var=" << env_variable << begin<<" "<<end<< endl;
03447   
03448 
03449   env_variable = env_variable.substr( begin+1, end-1 );
03450   // cout << "var=" << env_variable <<endl;
03451 
03452 
03453   // cerr<<"call getenv "<<endl;
03454   char* directory = getenv( env_variable.c_str() );
03455 
03456   if(!directory) {
03457     cerr<<"please define environment variable $"<<env_variable<<endl;
03458     delete this;
03459     exit(1);
03460   }
03461   string sdir = directory;
03462   sdir += "/";
03463 
03464   newString.replace( 0, lengh , sdir);
03465 
03466   if (verbosity_ == VERBOSE ) {
03467     cout << "expand " <<oldString<<" to "<< newString << endl;
03468   }
03469 
03470   return newString;
03471 }
03472 
03473 
03474 void 
03475 PFRootEventManager::printMCCalib(ofstream& out) const {
03476 
03477   if(!out) return;
03478   // if (!out.is_open()) return;
03479 
03480   // Use only for one PFSimParticle/GenParticles
03481   const HepMC::GenEvent* myGenEvent = MCTruth_.GetEvent();
03482   if(!myGenEvent) return;
03483   int nGen = 0;
03484   for ( HepMC::GenEvent::particle_const_iterator 
03485           piter  = myGenEvent->particles_begin();
03486           piter != myGenEvent->particles_end(); 
03487         ++piter ) nGen++;
03488   int nSim = trueParticles_.size();
03489   if ( nGen != 1 || nSim != 1 ) return;
03490 
03491   // One GenJet 
03492   if ( genJets_.size() != 1 ) return;
03493   double true_E = genJets_[0].p();
03494   double true_eta = genJets_[0].eta();
03495   double true_phi = genJets_[0].phi();
03496 
03497   // One particle-flow jet
03498   // if ( pfJets_.size() != 1 ) return;
03499   double rec_ECALEnergy = 0.;
03500   double rec_HCALEnergy = 0.;
03501   double deltaRMin = 999.;
03502   unsigned int theJet = 0;
03503   for ( unsigned int ijet=0; ijet<pfJets_.size(); ++ijet ) { 
03504     double rec_ECAL = pfJets_[ijet].neutralEmEnergy();
03505     double rec_HCAL = pfJets_[ijet].neutralHadronEnergy();
03506     double rec_eta = pfJets_[0].eta();
03507     double rec_phi = pfJets_[0].phi();
03508     double deltaR = std::sqrt( (rec_eta-true_eta)*(rec_eta-true_eta)
03509                              + (rec_phi-true_phi)*(rec_phi-true_phi) ); 
03510     if ( deltaR < deltaRMin ) { 
03511       deltaRMin = deltaR;
03512       rec_ECALEnergy = rec_ECAL;
03513       rec_HCALEnergy = rec_HCAL;
03514     }
03515   }
03516   if ( deltaRMin > 0.1 ) return;
03517   
03518   std::vector < PFCandidatePtr > constituents = pfJets_[theJet].getPFConstituents ();
03519   double pat_ECALEnergy = 0.;
03520   double pat_HCALEnergy = 0.;
03521   for (unsigned ic = 0; ic < constituents.size (); ++ic) {
03522     if ( constituents[ic]->particleId() < 4 ) continue;
03523     if ( constituents[ic]->particleId() == 4 ) 
03524       pat_ECALEnergy += constituents[ic]->rawEcalEnergy();
03525     else if ( constituents[ic]->particleId() == 5 ) 
03526       pat_HCALEnergy += constituents[ic]->rawHcalEnergy();
03527   }
03528 
03529   out << true_eta << " " << true_phi << " " << true_E 
03530       << " " <<  rec_ECALEnergy << " " << rec_HCALEnergy
03531       << " " <<  pat_ECALEnergy << " " << pat_HCALEnergy
03532       << " " << deltaRMin << std::endl;
03533 }
03534 
03535 void  PFRootEventManager::print(ostream& out,int maxNLines ) const {
03536 
03537   if(!out) return;
03538 
03539   //If McTruthMatching print a detailed list 
03540   //of matching between simparticles and PFCandidates
03541   //MCTruth Matching vectors.
03542   std::vector< std::list <simMatch> > candSimMatchTrack;
03543   std::vector< std::list <simMatch> >  candSimMatchEcal;  
03544   if( printMCTruthMatching_){
03545     mcTruthMatching( std::cout,
03546                      *pfCandidates_,
03547                      candSimMatchTrack,
03548                      candSimMatchEcal);
03549   }
03550 
03551 
03552   if( printRecHits_ ) {
03553     out<<"ECAL RecHits ==============================================="<<endl;
03554     printRecHits(rechitsECAL_, clusterAlgoECAL_, out );             out<<endl;
03555     out<<"HCAL RecHits ==============================================="<<endl;
03556     printRecHits(rechitsHCAL_, clusterAlgoHCAL_, out );             out<<endl;
03557     out<<"HFEM RecHits ==============================================="<<endl;
03558     printRecHits(rechitsHFEM_, clusterAlgoHFEM_, out );             out<<endl;
03559     out<<"HFHAD RecHits =============================================="<<endl;
03560     printRecHits(rechitsHFHAD_, clusterAlgoHFHAD_, out );           out<<endl;
03561     out<<"PS RecHits ================================================="<<endl;
03562     printRecHits(rechitsPS_, clusterAlgoPS_, out );                 out<<endl;
03563   }
03564 
03565   if( printClusters_ ) {
03566     out<<"ECAL Clusters ============================================="<<endl;
03567     printClusters( *clustersECAL_, out);                           out<<endl;
03568     out<<"HCAL Clusters ============================================="<<endl;
03569     printClusters( *clustersHCAL_, out);                           out<<endl;
03570     out<<"HFEM Clusters ============================================="<<endl;
03571     printClusters( *clustersHFEM_, out);                           out<<endl;
03572     out<<"HFHAD Clusters ============================================"<<endl;
03573     printClusters( *clustersHFHAD_, out);                          out<<endl;
03574     out<<"PS Clusters   ============================================="<<endl;
03575     printClusters( *clustersPS_, out);                             out<<endl;
03576   }
03577   bool printTracks = true;
03578   if( printTracks) {
03579     
03580   }
03581   if( printPFBlocks_ ) {
03582     out<<"Particle Flow Blocks ======================================"<<endl;
03583     for(unsigned i=0; i<pfBlocks_->size(); i++) {
03584       out<<i<<" "<<(*pfBlocks_)[i]<<endl;
03585     }    
03586     out<<endl;
03587   }
03588   if(printPFCandidates_) {
03589     out<<"Particle Flow Candidates =================================="<<endl;
03590     double mex = 0.;
03591     double mey = 0.;
03592     for(unsigned i=0; i<pfCandidates_->size(); i++) {
03593       const PFCandidate& pfc = (*pfCandidates_)[i];
03594       mex -= pfc.px();
03595       mey -= pfc.py();
03596       if(pfc.pt()>printPFCandidatesPtMin_)
03597       out<<i<<" " <<(*pfCandidates_)[i]<<endl;
03598     }    
03599     out<<endl;
03600     out<< "MEX, MEY, MET ============================================" <<endl 
03601        << mex << " " << mey << " " << sqrt(mex*mex+mey*mey);
03602     out<<endl;
03603     out<<endl;
03604 
03605     //print a detailed list of PFSimParticles matching
03606     //the PFCandiates
03607     if(printMCTruthMatching_){
03608       cout<<"MCTruthMatching Results"<<endl;
03609       for(unsigned icand=0; icand<pfCandidates_->size(); 
03610           icand++) {
03611         out <<icand<<" " <<(*pfCandidates_)[icand]<<endl;
03612         out << "is matching:" << endl;
03613 
03614         //tracking
03615         ITM it_t    = candSimMatchTrack[icand].begin();
03616         ITM itend_t = candSimMatchTrack[icand].end();
03617         for(;it_t!=itend_t;++it_t){
03618           unsigned simid = it_t->second;
03619           out << "\tSimParticle " << trueParticles_[simid]
03620               <<endl;
03621           out << "\t\tthrough Track matching pTrectrack=" 
03622               << it_t->first << " GeV" << endl;
03623         }//loop simparticles
03624 
03625         ITM it_e    = candSimMatchEcal[icand].begin();
03626         ITM itend_e = candSimMatchEcal[icand].end();
03627         for(;it_e!=itend_e;++it_e){
03628           unsigned simid = it_e->second;
03629           out << "\tSimParticle " << trueParticles_[simid]
03630               << endl; 
03631           out << "\t\tsimparticle contributing to a total of " 
03632               << it_e->first
03633               << " GeV of its ECAL cluster"
03634               << endl;  
03635         }//loop simparticles
03636         cout<<"________________"<<endl;
03637       }//loop candidates 
03638     }
03639   }
03640   if(printPFJets_) {
03641     out<<"Jets  ====================================================="<<endl;
03642     out<<"Particle Flow: "<<endl;
03643     for(unsigned i=0; i<pfJets_.size(); i++) {      
03644       if (pfJets_[i].pt() > printPFJetsPtMin_ )
03645         out<<i<<pfJets_[i].print()<<endl;
03646     }    
03647     out<<endl;
03648     out<<"Generated: "<<endl;
03649     for(unsigned i=0; i<genJets_.size(); i++) {
03650       if (genJets_[i].pt() > printPFJetsPtMin_ )
03651         out<<i<<genJets_[i].print()<<endl;
03652       // <<" invisible energy = "<<genJets_[i].invisibleEnergy()<<endl;
03653     }        
03654     out<<endl;
03655     out<<"Calo: "<<endl;
03656     for(unsigned i=0; i<caloJets_.size(); i++) {      
03657       out<<"pt = "<<caloJets_[i].pt()<<endl;
03658     }        
03659     out<<endl;  
03660   }
03661   if( printSimParticles_ ) {
03662     out<<"Sim Particles  ==========================================="<<endl;
03663 
03664     for(unsigned i=0; i<trueParticles_.size(); i++) {
03665       if( trackInsideGCut( trueParticles_[i]) ){ 
03666 
03667         const reco::PFSimParticle& ptc = trueParticles_[i];
03668 
03669         // get trajectory at start point
03670         const reco::PFTrajectoryPoint& tp0 = ptc.extrapolatedPoint( 0 );
03671 
03672         if(tp0.momentum().pt()>printSimParticlesPtMin_)
03673           out<<"\t"<<trueParticles_[i]<<endl;
03674       }
03675     }   
03676  
03677     //print a detailed list of PFSimParticles matching
03678     //the PFCandiates
03679     if(printMCTruthMatching_) {
03680       cout<<"MCTruthMatching Results"<<endl;
03681       for ( unsigned i=0;  i < trueParticles_.size(); i++) {
03682         cout << "==== Particle Simulated " << i << endl;
03683         const reco::PFSimParticle& ptc = trueParticles_[i];
03684         out <<i<<" "<<trueParticles_[i]<<endl;
03685         
03686         if(!ptc.daughterIds().empty()){
03687           cout << "Look at the desintegration products" << endl;
03688           cout << endl;
03689           continue;
03690         }
03691         
03692         //TRACKING
03693         if(ptc.rectrackId() != 99999){
03694           cout << "matching pfCandidate (trough tracking): " << endl;
03695           for( unsigned icand=0; icand<pfCandidates_->size()
03696                  ; icand++ ) 
03697             {
03698               ITM it    = candSimMatchTrack[icand].begin();
03699               ITM itend = candSimMatchTrack[icand].end();
03700               for(;it!=itend;++it)
03701                 if( i == it->second ){
03702                   out<<icand<<" "<<(*pfCandidates_)[icand]<<endl;
03703                   cout << endl;
03704                 }
03705             }//loop candidate
03706         }//trackmatch
03707         
03708         //CALORIMETRY
03709         vector<unsigned> rechitSimIDs  
03710           = ptc.recHitContrib();
03711         vector<double>   rechitSimFrac 
03712           = ptc.recHitContribFrac();
03713         //cout << "Number of rechits contrib =" << rechitSimIDs.size() << endl;
03714         if( !rechitSimIDs.size() ) continue; //no rechit
03715         
03716         cout << "matching pfCandidate (through ECAL): " << endl;
03717         
03718         //look at total ECAL desposition:
03719         double totalEcalE = 0.0;
03720         for(unsigned irh=0; irh<rechitsECAL_.size();++irh)
03721           for ( unsigned isimrh=0;  isimrh < rechitSimIDs.size(); 
03722                 isimrh++ )
03723             if(rechitSimIDs[isimrh] == rechitsECAL_[irh].detId())
03724               totalEcalE += (rechitsECAL_[irh].energy()*rechitSimFrac[isimrh]/100.0);
03725         cout << "For info, this particle deposits E=" << totalEcalE 
03726              << "(GeV) in the ECAL" << endl;
03727         
03728         for( unsigned icand=0; icand<pfCandidates_->size()
03729                ; icand++ ) 
03730           {
03731             ITM it    = candSimMatchEcal[icand].begin();
03732             ITM itend = candSimMatchEcal[icand].end();
03733             for(;it!=itend;++it)
03734               if( i == it->second )
03735                 out<<icand<<" "<<it->first<<"GeV "<<(*pfCandidates_)[icand]<<endl;        
03736           }//loop candidate
03737         cout << endl;      
03738       }//loop particles  
03739     }//mctruthmatching
03740 
03741   }
03742 
03743   
03744   if ( printGenParticles_ ) { 
03745     printGenParticles(out,maxNLines);
03746   }
03747 }
03748 
03749 
03750 void
03751 PFRootEventManager::printGenParticles(std::ostream& out,
03752                                       int maxNLines) const {
03753                                  
03754                                  
03755   const HepMC::GenEvent* myGenEvent = MCTruth_.GetEvent();
03756   if(!myGenEvent) return;
03757 
03758   out<<"GenParticles ==========================================="<<endl;
03759 
03760   std::cout << "Id  Gen Name       eta    phi     pT     E    Vtx1   " 
03761             << " x      y      z   " 
03762             << "Moth  Vtx2  eta   phi     R      Z   Da1  Da2 Ecal?" 
03763             << std::endl;
03764 
03765   int nLines = 0;
03766   for ( HepMC::GenEvent::particle_const_iterator 
03767           piter  = myGenEvent->particles_begin();
03768         piter != myGenEvent->particles_end(); 
03769         ++piter ) {
03770     
03771     if( nLines == maxNLines) break;
03772     nLines++;
03773     
03774     HepMC::GenParticle* p = *piter;
03775     /* */
03776     int partId = p->pdg_id();
03777 
03778     // We have here a subset of particles only. 
03779     // To be filled according to the needs.
03780     /*switch(partId) {
03781       case    1: { name = "d"; break; } 
03782       case    2: { name = "u"; break; } 
03783       case    3: { name = "s"; break; } 
03784       case    4: { name = "c"; break; } 
03785       case    5: { name = "b"; break; } 
03786       case    6: { name = "t"; break; } 
03787       case   -1: { name = "~d"; break; } 
03788       case   -2: { name = "~u"; break; } 
03789       case   -3: { name = "~s"; break; } 
03790       case   -4: { name = "~c"; break; } 
03791       case   -5: { name = "~b"; break; } 
03792       case   -6: { name = "~t"; break; } 
03793       case   11: { name = "e-"; break; }
03794       case  -11: { name = "e+"; break; }
03795       case   12: { name = "nu_e"; break; }
03796       case  -12: { name = "~nu_e"; break; }
03797       case   13: { name = "mu-"; break; }
03798       case  -13: { name = "mu+"; break; }
03799       case   14: { name = "nu_mu"; break; }
03800       case  -14: { name = "~nu_mu"; break; }
03801       case   15: { name = "tau-"; break; }
03802       case  -15: { name = "tau+"; break; }
03803       case   16: { name = "nu_tau"; break; }
03804       case  -16: { name = "~nu_tau"; break; }
03805       case   21: { name = "gluon"; break; }
03806       case   22: { name = "gamma"; break; }
03807       case   23: { name = "Z0"; break; }
03808       case   24: { name = "W+"; break; }
03809       case   25: { name = "H0"; break; }
03810       case  -24: { name = "W-"; break; }
03811       case  111: { name = "pi0"; break; }
03812       case  113: { name = "rho0"; break; }
03813       case  223: { name = "omega"; break; }
03814       case  333: { name = "phi"; break; }
03815       case  443: { name = "J/psi"; break; }
03816       case  553: { name = "Upsilon"; break; }
03817       case  130: { name = "K0L"; break; }
03818       case  211: { name = "pi+"; break; }
03819       case -211: { name = "pi-"; break; }
03820       case  213: { name = "rho+"; break; }
03821       case -213: { name = "rho-"; break; }
03822       case  221: { name = "eta"; break; }
03823       case  331: { name = "eta'"; break; }
03824       case  441: { name = "etac"; break; }
03825       case  551: { name = "etab"; break; }
03826       case  310: { name = "K0S"; break; }
03827       case  311: { name = "K0"; break; }
03828       case -311: { name = "Kbar0"; break; }
03829       case  321: { name = "K+"; break; }
03830       case -321: { name = "K-"; break; }
03831       case  411: { name = "D+"; break; }
03832       case -411: { name = "D-"; break; }
03833       case  421: { name = "D0"; break; }
03834       case  431: { name = "Ds_+"; break; }
03835       case -431: { name = "Ds_-"; break; }
03836       case  511: { name = "B0"; break; }
03837       case  521: { name = "B+"; break; }
03838       case -521: { name = "B-"; break; }
03839       case  531: { name = "Bs_0"; break; }
03840       case  541: { name = "Bc_+"; break; }
03841       case -541: { name = "Bc_+"; break; }
03842       case  313: { name = "K*0"; break; }
03843       case -313: { name = "K*bar0"; break; }
03844       case  323: { name = "K*+"; break; }
03845       case -323: { name = "K*-"; break; }
03846       case  413: { name = "D*+"; break; }
03847       case -413: { name = "D*-"; break; }
03848       case  423: { name = "D*0"; break; }
03849       case  513: { name = "B*0"; break; }
03850       case  523: { name = "B*+"; break; }
03851       case -523: { name = "B*-"; break; }
03852       case  533: { name = "B*_s0"; break; }
03853       case  543: { name = "B*_c+"; break; }
03854       case -543: { name = "B*_c-"; break; }
03855       case  1114: { name = "Delta-"; break; }
03856       case -1114: { name = "Deltabar+"; break; }
03857       case -2112: { name = "nbar0"; break; }
03858       case  2112: { name = "n"; break; }
03859       case  2114: { name = "Delta0"; break; }
03860       case -2114: { name = "Deltabar0"; break; }
03861       case  3122: { name = "Lambda0"; break; }
03862       case -3122: { name = "Lambdabar0"; break; }
03863       case  3112: { name = "Sigma-"; break; }
03864       case -3112: { name = "Sigmabar+"; break; }
03865       case  3212: { name = "Sigma0"; break; }
03866       case -3212: { name = "Sigmabar0"; break; }
03867       case  3214: { name = "Sigma*0"; break; }
03868       case -3214: { name = "Sigma*bar0"; break; }
03869       case  3222: { name = "Sigma+"; break; }
03870       case -3222: { name = "Sigmabar-"; break; }
03871       case  2212: { name = "p"; break; }
03872       case -2212: { name = "~p"; break; }
03873       case -2214: { name = "Delta-"; break; }
03874       case  2214: { name = "Delta+"; break; }
03875       case -2224: { name = "Deltabar--"; break; }
03876       case  2224: { name = "Delta++"; break; }
03877       default: { 
03878       name = "unknown"; 
03879       cout << "Unknown code : " << partId << endl;
03880       }   
03881       }
03882     */
03883     std::string latexString;
03884     std::string name = getGenParticleName(partId,latexString);
03885 
03886     math::XYZTLorentzVector momentum1(p->momentum().px(),
03887                                       p->momentum().py(),
03888                                       p->momentum().pz(),
03889                                       p->momentum().e() );
03890 
03891     if(momentum1.pt()<printGenParticlesPtMin_) continue;
03892 
03893     int vertexId1 = 0;
03894 
03895     if ( !p->production_vertex() && p->pdg_id() == 2212 ) continue;
03896 
03897     math::XYZVector vertex1;
03898     vertexId1 = -1;
03899 
03900     if(p->production_vertex() ) {
03901       vertex1.SetCoordinates( p->production_vertex()->position().x()/10.,
03902                               p->production_vertex()->position().y()/10.,
03903                               p->production_vertex()->position().z()/10. );
03904       vertexId1 = p->production_vertex()->barcode();
03905     }
03906 
03907     out.setf(std::ios::fixed, std::ios::floatfield);
03908     out.setf(std::ios::right, std::ios::adjustfield);
03909     
03910     out << std::setw(4) << p->barcode() << " " 
03911         << name;
03912     
03913     for(unsigned int k=0;k<11-name.length() && k<12; k++) out << " ";  
03914     
03915     double eta = momentum1.eta();
03916     if ( eta > +10. ) eta = +10.;
03917     if ( eta < -10. ) eta = -10.;
03918     
03919     out << std::setw(6) << std::setprecision(2) << eta << " " 
03920         << std::setw(6) << std::setprecision(2) << momentum1.phi() << " " 
03921         << std::setw(7) << std::setprecision(2) << momentum1.pt() << " " 
03922         << std::setw(7) << std::setprecision(2) << momentum1.e() << " " 
03923         << std::setw(4) << vertexId1 << " " 
03924         << std::setw(6) << std::setprecision(1) << vertex1.x() << " " 
03925         << std::setw(6) << std::setprecision(1) << vertex1.y() << " " 
03926         << std::setw(6) << std::setprecision(1) << vertex1.z() << " ";
03927 
03928 
03929     if( p->production_vertex() ) {
03930       if ( p->production_vertex()->particles_in_size() ) {
03931         const HepMC::GenParticle* mother = 
03932           *(p->production_vertex()->particles_in_const_begin());
03933         
03934         out << std::setw(4) << mother->barcode() << " ";
03935       }
03936       else 
03937         out << "     " ;
03938     }    
03939 
03940     if ( p->end_vertex() ) {  
03941       math::XYZTLorentzVector vertex2(p->end_vertex()->position().x()/10.,
03942                                       p->end_vertex()->position().y()/10.,
03943                                       p->end_vertex()->position().z()/10.,
03944                                       p->end_vertex()->position().t()/10.);
03945       int vertexId2 = p->end_vertex()->barcode();
03946 
03947       std::vector<const HepMC::GenParticle*> children;
03948       HepMC::GenVertex::particles_out_const_iterator firstDaughterIt = 
03949         p->end_vertex()->particles_out_const_begin();
03950       HepMC::GenVertex::particles_out_const_iterator lastDaughterIt = 
03951         p->end_vertex()->particles_out_const_end();
03952       for ( ; firstDaughterIt != lastDaughterIt ; ++firstDaughterIt ) {
03953         children.push_back(*firstDaughterIt);
03954       }      
03955 
03956       out << std::setw(4) << vertexId2 << " "
03957           << std::setw(6) << std::setprecision(2) << vertex2.eta() << " " 
03958           << std::setw(6) << std::setprecision(2) << vertex2.phi() << " " 
03959           << std::setw(5) << std::setprecision(1) << vertex2.pt() << " " 
03960           << std::setw(6) << std::setprecision(1) << vertex2.z() << " ";
03961 
03962       for ( unsigned id=0; id<children.size(); ++id )
03963         out << std::setw(4) << children[id]->barcode() << " ";
03964     }
03965     out << std::endl;
03966   }
03967 }
03968 
03969 void PFRootEventManager::printRecHits(const reco::PFRecHitCollection& rechits, const PFClusterAlgo& clusterAlgo, ostream& out) const{
03970 
03971     for(unsigned i=0; i<rechits.size(); i++) {
03972       string seedstatus = "    ";
03973       if(clusterAlgo.isSeed(i) ) 
03974         seedstatus = "SEED";
03975       printRecHit(rechits[i], i, seedstatus.c_str(), out);
03976     }
03977     return;
03978 }
03979 
03980 void  PFRootEventManager::printRecHit(const reco::PFRecHit& rh,
03981                                       unsigned index,  
03982                                       const char* seedstatus,
03983                                       ostream& out) const {
03984 
03985   if(!out) return;
03986   double eta = rh.position().Eta();
03987   double phi = rh.position().Phi();
03988   double energy = rh.energy();
03989 
03990   if(energy<printRecHitsEMin_)  return;
03991 
03992   TCutG* cutg = (TCutG*) gROOT->FindObject("CUTG");
03993   if( !cutg || cutg->IsInside( eta, phi ) ) 
03994     out<<index<<"\t"<<seedstatus<<" "<<rh<<endl; 
03995 }
03996 
03997 void PFRootEventManager::printClusters(const reco::PFClusterCollection& clusters,
03998                                        ostream& out ) const {  
03999   for(unsigned i=0; i<clusters.size(); i++) {
04000     printCluster(clusters[i], out);
04001   }
04002   return;
04003 }
04004 
04005 void  PFRootEventManager::printCluster(const reco::PFCluster& cluster,
04006                                        ostream& out ) const {
04007   
04008   if(!out) return;
04009 
04010   double eta = cluster.position().Eta();
04011   double phi = cluster.position().Phi();
04012   double energy = cluster.energy();
04013 
04014   if(energy<printClustersEMin_)  return;
04015 
04016   TCutG* cutg = (TCutG*) gROOT->FindObject("CUTG");
04017   if( !cutg || cutg->IsInside( eta, phi ) ) 
04018     out<<cluster<<endl;
04019 }
04020 
04021 
04022 
04023 
04024 
04025 bool PFRootEventManager::trackInsideGCut( const reco::PFTrack& track ) const {
04026 
04027   TCutG* cutg = (TCutG*) gROOT->FindObject("CUTG");
04028   if(!cutg) return true;
04029   
04030   const vector< reco::PFTrajectoryPoint >& points = track.trajectoryPoints();
04031   
04032   for( unsigned i=0; i<points.size(); i++) {
04033     if( ! points[i].isValid() ) continue;
04034     
04035     const math::XYZPoint& pos = points[i].position();
04036     if( cutg->IsInside( pos.Eta(), pos.Phi() ) ) return true;
04037   }
04038 
04039   // no point inside cut
04040   return false;
04041 }
04042 
04043 
04044 void  
04045 PFRootEventManager::fillRecHitMask( vector<bool>& mask, 
04046                                     const reco::PFRecHitCollection& rechits ) 
04047   const {
04048 
04049   TCutG* cutg = (TCutG*) gROOT->FindObject("CUTG");
04050   if(!cutg) {
04051     mask.resize( rechits.size(), true);
04052     return;
04053   }
04054 
04055   mask.clear();
04056   mask.reserve( rechits.size() );
04057   for(unsigned i=0; i<rechits.size(); i++) {
04058     
04059     double eta = rechits[i].position().Eta();
04060     double phi = rechits[i].position().Phi();
04061 
04062     if( cutg->IsInside( eta, phi ) )
04063       mask.push_back( true );
04064     else 
04065       mask.push_back( false );   
04066   }
04067 }
04068 
04069 void  
04070 PFRootEventManager::fillClusterMask(vector<bool>& mask, 
04071                                     const reco::PFClusterCollection& clusters) 
04072   const {
04073   
04074   TCutG* cutg = (TCutG*) gROOT->FindObject("CUTG");
04075   if(!cutg) return;
04076 
04077   mask.clear();
04078   mask.reserve( clusters.size() );
04079   for(unsigned i=0; i<clusters.size(); i++) {
04080     
04081     double eta = clusters[i].position().Eta();
04082     double phi = clusters[i].position().Phi();
04083 
04084     if( cutg->IsInside( eta, phi ) )
04085       mask.push_back( true );
04086     else 
04087       mask.push_back( false );   
04088   }
04089 }
04090 
04091 void  
04092 PFRootEventManager::fillTrackMask(vector<bool>& mask, 
04093                                   const reco::PFRecTrackCollection& tracks) 
04094   const {
04095   
04096   TCutG* cutg = (TCutG*) gROOT->FindObject("CUTG");
04097   if(!cutg) return;
04098 
04099   mask.clear();
04100   mask.reserve( tracks.size() );
04101   for(unsigned i=0; i<tracks.size(); i++) {
04102     if( trackInsideGCut( tracks[i] ) )
04103       mask.push_back( true );
04104     else 
04105       mask.push_back( false );   
04106   }
04107 }
04108 
04109 void  
04110 PFRootEventManager::fillPhotonMask(vector<bool>& mask, 
04111                                   const reco::PhotonCollection& photons) 
04112   const {
04113   
04114   TCutG* cutg = (TCutG*) gROOT->FindObject("CUTG");
04115   if(!cutg) return;
04116 
04117   mask.clear();
04118   mask.reserve( photons.size() );
04119   for(unsigned i=0; i<photons.size(); i++) {
04120     double eta = photons[i].caloPosition().Eta();
04121     double phi = photons[i].caloPosition().Phi();
04122     if( cutg->IsInside( eta, phi ) )
04123       mask.push_back( true );
04124     else 
04125       mask.push_back( false );   
04126   }
04127 }
04128 
04129 
04130 void  
04131 PFRootEventManager::fillTrackMask(vector<bool>& mask, 
04132                                   const reco::GsfPFRecTrackCollection& tracks) 
04133   const {
04134   
04135   TCutG* cutg = (TCutG*) gROOT->FindObject("CUTG");
04136   if(!cutg) return;
04137 
04138   mask.clear();
04139   mask.reserve( tracks.size() );
04140   for(unsigned i=0; i<tracks.size(); i++) {
04141     if( trackInsideGCut( tracks[i] ) )
04142       mask.push_back( true );
04143     else 
04144       mask.push_back( false );   
04145   }
04146 }
04147 
04148 
04149 const reco::PFSimParticle&
04150 PFRootEventManager::closestParticle( reco::PFTrajectoryPoint::LayerType layer, 
04151                                      double eta, double phi,
04152                                      double& peta, double& pphi, double& pe) 
04153   const {
04154   
04155 
04156   if( trueParticles_.empty() ) {
04157     string err  = "PFRootEventManager::closestParticle : ";
04158     err        += "vector of PFSimParticles is empty";
04159     throw std::length_error( err.c_str() );
04160   }
04161 
04162   double mindist2 = 99999999;
04163   unsigned iClosest=0;
04164   for(unsigned i=0; i<trueParticles_.size(); i++) {
04165     
04166     const reco::PFSimParticle& ptc = trueParticles_[i];
04167 
04168     // protection for old version of the PFSimParticle 
04169     // dataformats. 
04170     if( layer >= reco::PFTrajectoryPoint::NLayers ||
04171         ptc.nTrajectoryMeasurements() + layer >= 
04172         ptc.nTrajectoryPoints() ) {
04173       continue;
04174     }
04175 
04176     const reco::PFTrajectoryPoint& tp
04177       = ptc.extrapolatedPoint( layer );
04178 
04179     peta = tp.position().Eta();
04180     pphi = tp.position().Phi();
04181     pe = tp.momentum().E();
04182 
04183     double deta = peta - eta;
04184     double dphi = pphi - phi;
04185 
04186     double dist2 = deta*deta + dphi*dphi;
04187 
04188     if(dist2<mindist2) {
04189       mindist2 = dist2;
04190       iClosest = i;
04191     }
04192   }
04193 
04194   return trueParticles_[iClosest];
04195 }
04196 
04197 
04198 
04199 //-----------------------------------------------------------
04200 void 
04201 PFRootEventManager::readCMSSWJets() {
04202 
04203   cout<<"CMSSW Gen jets : size = " <<  genJetsCMSSW_.size() << endl;
04204   for ( unsigned i = 0; i < genJetsCMSSW_.size(); i++) {
04205     cout<<"Gen jet Et : " <<  genJetsCMSSW_[i].et() << endl;
04206   }
04207   cout<<"CMSSW PF jets : size = " <<  pfJetsCMSSW_.size() << endl;
04208   for ( unsigned i = 0; i < pfJetsCMSSW_.size(); i++) {
04209     cout<<"PF jet Et : " <<  pfJetsCMSSW_[i].et() << endl;
04210   }
04211   cout<<"CMSSW Calo jets : size = " <<  caloJetsCMSSW_.size() << endl;
04212   for ( unsigned i = 0; i < caloJetsCMSSW_.size(); i++) {
04213     cout<<"Calo jet Et : " << caloJetsCMSSW_[i].et() << endl;
04214   }
04215 }
04216 //________________________________________________________________
04217 std::string PFRootEventManager::getGenParticleName(int partId, std::string &latexString) const
04218 {
04219   std::string  name;
04220   switch(partId) {
04221   case    1: { name = "d";latexString="d"; break; } 
04222   case    2: { name = "u";latexString="u";break; } 
04223   case    3: { name = "s";latexString="s" ;break; } 
04224   case    4: { name = "c";latexString="c" ; break; } 
04225   case    5: { name = "b";latexString="b" ; break; } 
04226   case    6: { name = "t";latexString="t" ; break; } 
04227   case   -1: { name = "~d";latexString="#bar{d}" ; break; } 
04228   case   -2: { name = "~u";latexString="#bar{u}" ; break; } 
04229   case   -3: { name = "~s";latexString="#bar{s}" ; break; } 
04230   case   -4: { name = "~c";latexString="#bar{c}" ; break; } 
04231   case   -5: { name = "~b";latexString="#bar{b}" ; break; } 
04232   case   -6: { name = "~t";latexString="#bar{t}" ; break; } 
04233   case   11: { name = "e-";latexString=name ; break; }
04234   case  -11: { name = "e+";latexString=name ; break; }
04235   case   12: { name = "nu_e";latexString="#nu_{e}" ; break; }
04236   case  -12: { name = "~nu_e";latexString="#bar{#nu}_{e}" ; break; }
04237   case   13: { name = "mu-";latexString="#mu-" ; break; }
04238   case  -13: { name = "mu+";latexString="#mu+" ; break; }
04239   case   14: { name = "nu_mu";latexString="#nu_{mu}" ; break; }
04240   case  -14: { name = "~nu_mu";latexString="#bar{#nu}_{#mu}"; break; }
04241   case   15: { name = "tau-";latexString="#tau^{-}" ; break; }
04242   case  -15: { name = "tau+";latexString="#tau^{+}" ; break; }
04243   case   16: { name = "nu_tau";latexString="#nu_{#tau}" ; break; }
04244   case  -16: { name = "~nu_tau";latexString="#bar{#nu}_{#tau}"; break; }
04245   case   21: { name = "gluon";latexString= name; break; }
04246   case   22: { name = "gamma";latexString= "#gamma"; break; }
04247   case   23: { name = "Z0";latexString="Z^{0}" ; break; }
04248   case   24: { name = "W+";latexString="W^{+}" ; break; }
04249   case   25: { name = "H0";latexString=name ; break; }
04250   case  -24: { name = "W-";latexString="W^{-}" ; break; }
04251   case  111: { name = "pi0";latexString="#pi^{0}" ; break; }
04252   case  113: { name = "rho0";latexString="#rho^{0}" ; break; }
04253   case  223: { name = "omega";latexString="#omega" ; break; }
04254   case  333: { name = "phi";latexString= "#phi"; break; }
04255   case  443: { name = "J/psi";latexString="J/#psi" ; break; }
04256   case  553: { name = "Upsilon";latexString="#Upsilon" ; break; }
04257   case  130: { name = "K0L";latexString=name ; break; }
04258   case  211: { name = "pi+";latexString="#pi^{+}" ; break; }
04259   case -211: { name = "pi-";latexString="#pi^{-}" ; break; }
04260   case  213: { name = "rho+";latexString="#rho^{+}" ; break; }
04261   case -213: { name = "rho-";latexString="#rho^{-}" ; break; }
04262   case  221: { name = "eta";latexString="#eta" ; break; }
04263   case  331: { name = "eta'";latexString="#eta'" ; break; }
04264   case  441: { name = "etac";latexString="#eta_{c}" ; break; }
04265   case  551: { name = "etab";latexString= "#eta_{b}"; break; }
04266   case  310: { name = "K0S";latexString=name ; break; }
04267   case  311: { name = "K0";latexString="K^{0}" ; break; }
04268   case -311: { name = "Kbar0";latexString="#bar{#Kappa}^{0}" ; break; }
04269   case  321: { name = "K+";latexString= "K^{+}"; break; }
04270   case -321: { name = "K-";latexString="K^{-}"; break; }
04271   case  411: { name = "D+";latexString="D^{+}" ; break; }
04272   case -411: { name = "D-";latexString="D^{-}"; break; }
04273   case  421: { name = "D0";latexString=name ; break; }
04274   case  431: { name = "Ds_+";latexString="Ds_{+}" ; break; }
04275   case -431: { name = "Ds_-";latexString="Ds_{-}" ; break; }
04276   case  511: { name = "B0";latexString= name; break; }
04277   case  521: { name = "B+";latexString="B^{+}" ; break; }
04278   case -521: { name = "B-";latexString="B^{-}" ; break; }
04279   case  531: { name = "Bs_0";latexString="Bs_{0}" ; break; }
04280   case  541: { name = "Bc_+";latexString="Bc_{+}" ; break; }
04281   case -541: { name = "Bc_+";latexString="Bc_{+}" ; break; }
04282   case  313: { name = "K*0";latexString="K^{*0}" ; break; }
04283   case -313: { name = "K*bar0";latexString="#bar{K}^{*0}" ; break; }
04284   case  323: { name = "K*+";latexString="#K^{*+}"; break; }
04285   case -323: { name = "K*-";latexString="#K^{*-}" ; break; }
04286   case  413: { name = "D*+";latexString= "D^{*+}"; break; }
04287   case -413: { name = "D*-";latexString= "D^{*-}" ; break; }
04288   case  423: { name = "D*0";latexString="D^{*0}" ; break; }
04289   case  513: { name = "B*0";latexString="B^{*0}" ; break; }
04290   case  523: { name = "B*+";latexString="B^{*+}" ; break; }
04291   case -523: { name = "B*-";latexString="B^{*-}" ; break; }
04292   case  533: { name = "B*_s0";latexString="B^{*}_{s0}" ; break; }
04293   case  543: { name = "B*_c+";latexString= "B^{*}_{c+}"; break; }
04294   case -543: { name = "B*_c-";latexString= "B^{*}_{c-}"; break; }
04295   case  1114: { name = "Delta-";latexString="#Delta^{-}" ; break; }
04296   case -1114: { name = "Deltabar+";latexString="#bar{#Delta}^{+}" ; break; }
04297   case -2112: { name = "nbar0";latexString="{bar}n^{0}" ; break; }
04298   case  2112: { name = "n"; latexString=name ;break;}
04299   case  2114: { name = "Delta0"; latexString="#Delta^{0}" ;break; }
04300   case -2114: { name = "Deltabar0"; latexString="#bar{#Delta}^{0}" ;break; }
04301   case  3122: { name = "Lambda0";latexString= "#Lambda^{0}"; break; }
04302   case -3122: { name = "Lambdabar0";latexString="#bar{#Lambda}^{0}" ; break; }
04303   case  3112: { name = "Sigma-"; latexString="#Sigma" ;break; }
04304   case -3112: { name = "Sigmabar+"; latexString="#bar{#Sigma}^{+}" ;break; }
04305   case  3212: { name = "Sigma0";latexString="#Sigma^{0}" ; break; }
04306   case -3212: { name = "Sigmabar0";latexString="#bar{#Sigma}^{0}" ; break; }
04307   case  3214: { name = "Sigma*0"; latexString="#Sigma^{*0}" ;break; }
04308   case -3214: { name = "Sigma*bar0";latexString="#bar{#Sigma}^{*0}" ; break; }
04309   case  3222: { name = "Sigma+"; latexString="#Sigma^{+}" ;break; }
04310   case -3222: { name = "Sigmabar-"; latexString="#bar{#Sigma}^{-}";break; }
04311   case  2212: { name = "p";latexString=name ; break; }
04312   case -2212: { name = "~p";latexString="#bar{p}" ; break; }
04313   case -2214: { name = "Delta-";latexString="#Delta^{-}" ; break; }
04314   case  2214: { name = "Delta+";latexString="#Delta^{+}" ; break; }
04315   case -2224: { name = "Deltabar--"; latexString="#bar{#Delta}^{--}" ;break; }
04316   case  2224: { name = "Delta++"; latexString= "#Delta^{++}";break; }
04317   default:
04318     {
04319       name = "unknown"; 
04320       cout << "Unknown code : " << partId << endl;
04321       break;
04322     } 
04323                 
04324                   
04325   }
04326   return name;  
04327 
04328 }
04329 
04330 //_____________________________________________________________________________
04331 void PFRootEventManager::mcTruthMatching( std::ostream& out,
04332                                           const reco::PFCandidateCollection& candidates,
04333                                           std::vector< std::list <simMatch> >& candSimMatchTrack,
04334                                           std::vector< std::list <simMatch> >& candSimMatchEcal) const
04335 {
04336   
04337   if(!out) return;
04338   out << endl;
04339   out << "Running Monte Carlo Truth Matching Tool" << endl;
04340   out << endl;
04341 
04342   //resize matching vectors
04343   candSimMatchTrack.resize(candidates.size());
04344   candSimMatchEcal.resize(candidates.size());
04345 
04346   for(unsigned i=0; i<candidates.size(); i++) {
04347     const reco::PFCandidate& pfCand = candidates[i];
04348     
04349     //Matching with ECAL clusters
04350     if (verbosity_ == VERBOSE ) {
04351       out <<i<<" " <<(*pfCandidates_)[i]<<endl;
04352       out << "is matching:" << endl;
04353     }
04354     
04355     PFCandidate::ElementsInBlocks eleInBlocks 
04356       = pfCand.elementsInBlocks();
04357 
04358     for(unsigned iel=0; iel<eleInBlocks.size(); ++iel) {
04359       PFBlockRef blockRef   = eleInBlocks[iel].first;
04360       unsigned indexInBlock = eleInBlocks[iel].second;
04361       
04362       //Retrieving elements of the block
04363       const reco::PFBlock& blockh 
04364         = *blockRef;
04365       const edm::OwnVector< reco::PFBlockElement >& 
04366         elements_h = blockh.elements();
04367       
04368       reco::PFBlockElement::Type type 
04369         = elements_h[ indexInBlock ].type();   
04370 //       cout <<"(" << blockRef.key() << "|" <<indexInBlock <<"|" 
04371 //         << elements_h[ indexInBlock ].type() << ")," << endl;
04372       
04373       //TRACK=================================
04374       if(type == reco::PFBlockElement::TRACK){
04375         const reco::PFRecTrackRef trackref 
04376           = elements_h[ indexInBlock ].trackRefPF();
04377         assert( !trackref.isNull() );     
04378         const reco::PFRecTrack& track = *trackref; 
04379         const reco::TrackRef trkREF = track.trackRef();
04380         unsigned rtrkID = track.trackId();
04381 
04382         //looking for the matching charged simulated particle:
04383         for ( unsigned isim=0;  isim < trueParticles_.size(); isim++) {
04384           const reco::PFSimParticle& ptc = trueParticles_[isim];
04385           unsigned trackIDM = ptc.rectrackId();
04386           if(trackIDM != 99999 
04387              && trackIDM == rtrkID){
04388 
04389             if (verbosity_ == VERBOSE ) 
04390               out << "\tSimParticle " << isim 
04391                   << " through Track matching pTrectrack=" 
04392                   << trkREF->pt() << " GeV" << endl;     
04393             
04394             //store info
04395             std::pair<double, unsigned> simtrackmatch
04396               = make_pair(trkREF->pt(),trackIDM);
04397             candSimMatchTrack[i].push_back(simtrackmatch);
04398           }//match
04399         }//loop simparticles 
04400         
04401       }//TRACK
04402 
04403       //ECAL=================================
04404       if(type == reco::PFBlockElement::ECAL)
04405         {
04406           const reco::PFClusterRef clusterref 
04407             = elements_h[ indexInBlock ].clusterRef();
04408           assert( !clusterref.isNull() );         
04409           const reco::PFCluster& cluster = *clusterref; 
04410           
04411           const std::vector< reco::PFRecHitFraction >& 
04412             fracs = cluster.recHitFractions();  
04413 
04414 //        cout << "This is an ecal cluster of energy " 
04415 //             << cluster.energy() << endl;
04416           vector<unsigned> simpID;
04417           vector<double>   simpEC(trueParticles_.size(),0.0);     
04418           vector<unsigned> simpCN(trueParticles_.size(),0);      
04419           for(unsigned int rhit = 0; rhit < fracs.size(); ++rhit){
04420             
04421             const reco::PFRecHitRef& rh = fracs[rhit].recHitRef();
04422             if(rh.isNull()) continue;
04423             const reco::PFRecHit& rechit_cluster = *rh;
04424 //          cout << rhit << " ID=" << rechit_cluster.detId() 
04425 //               << " E=" << rechit_cluster.energy() 
04426 //               << " fraction=" << fracs[rhit].fraction() << " ";
04427             
04428             //loop on sim particules
04429 //          cout << "coming from sim particles: ";
04430             for ( unsigned isim=0;  isim < trueParticles_.size(); isim++) {
04431               const reco::PFSimParticle& ptc = trueParticles_[isim];
04432               
04433               vector<unsigned> rechitSimIDs  
04434                 = ptc.recHitContrib();
04435               vector<double>   rechitSimFrac 
04436                 = ptc.recHitContribFrac();
04437               //cout << "Number of rechits contrib =" << rechitSimIDs.size() << endl;
04438               if( !rechitSimIDs.size() ) continue; //no rechit
04439                                                                        
04440               for ( unsigned isimrh=0;  isimrh < rechitSimIDs.size(); isimrh++) {
04441                 if( rechitSimIDs[isimrh] == rechit_cluster.detId() ){
04442                   
04443                   bool takenalready = false;
04444                   for(unsigned iss = 0; iss < simpID.size(); ++iss)
04445                     if(simpID[iss] == isim) takenalready = true;
04446                   if(!takenalready) simpID.push_back(isim);
04447                   
04448                   simpEC[isim] += 
04449                     ((rechit_cluster.energy()*rechitSimFrac[isimrh])/100.0)
04450                     *fracs[rhit].fraction();
04451                   
04452                   simpCN[isim]++; //counting rechits
04453 
04454 //                cout << isim << " with contribution of =" 
04455 //                     << rechitSimFrac[isimrh] << "%, "; 
04456                 }//match rechit
04457               }//loop sim rechit
04458             }//loop sim particules
04459 //          cout << endl;
04460           }//loop cand rechit 
04461 
04462           for(unsigned is=0; is < simpID.size(); ++is)
04463             {
04464               double frac_of_cluster 
04465                 = (simpEC[simpID[is]]/cluster.energy())*100.0;
04466               
04467               //store info
04468               std::pair<double, unsigned> simecalmatch
04469                 = make_pair(simpEC[simpID[is]],simpID[is]);
04470               candSimMatchEcal[i].push_back(simecalmatch);
04471               
04472               if (verbosity_ == VERBOSE ) {
04473                 out << "\tSimParticle " << simpID[is] 
04474                     << " through ECAL matching Epfcluster=" 
04475                     << cluster.energy() 
04476                     << " GeV with N=" << simpCN[simpID[is]]
04477                     << " rechits in common "
04478                     << endl; 
04479                 out << "\t\tsimparticle contributing to a total of " 
04480                     << simpEC[simpID[is]]
04481                     << " GeV of this cluster (" 
04482                     <<  frac_of_cluster << "%) " 
04483                     << endl;
04484               }
04485             }//loop particle matched
04486         }//ECAL clusters
04487 
04488     }//loop elements
04489 
04490     if (verbosity_ == VERBOSE )
04491       cout << "===============================================================" 
04492            << endl;
04493 
04494   }//loop pfCandidates_
04495 
04496   if (verbosity_ == VERBOSE ){
04497 
04498     cout << "=================================================================="
04499          << endl;
04500     cout << "SimParticles" << endl;
04501     
04502     //loop simulated particles  
04503     for ( unsigned i=0;  i < trueParticles_.size(); i++) {
04504       cout << "==== Particle Simulated " << i << endl;
04505       const reco::PFSimParticle& ptc = trueParticles_[i];
04506       out <<i<<" "<<trueParticles_[i]<<endl;
04507 
04508       if(!ptc.daughterIds().empty()){
04509         cout << "Look at the desintegration products" << endl;
04510         cout << endl;
04511         continue;
04512       }
04513       
04514       //TRACKING
04515       if(ptc.rectrackId() != 99999){
04516         cout << "matching pfCandidate (trough tracking): " << endl;
04517         for( unsigned icand=0; icand<candidates.size(); icand++ ) 
04518           {
04519             ITM it    = candSimMatchTrack[icand].begin();
04520             ITM itend = candSimMatchTrack[icand].end();
04521             for(;it!=itend;++it)
04522               if( i == it->second ){
04523                 out<<icand<<" "<<(*pfCandidates_)[icand]<<endl;
04524                 cout << endl;
04525               }
04526           }//loop candidate
04527       }//trackmatch
04528       
04529       
04530       //CALORIMETRY
04531       vector<unsigned> rechitSimIDs  
04532         = ptc.recHitContrib();
04533       vector<double>   rechitSimFrac 
04534         = ptc.recHitContribFrac();
04535       //cout << "Number of rechits contrib =" << rechitSimIDs.size() << endl;
04536       if( !rechitSimIDs.size() ) continue; //no rechit
04537       
04538       cout << "matching pfCandidate (through ECAL): " << endl;
04539       
04540       //look at total ECAL desposition:
04541       double totalEcalE = 0.0;
04542       for(unsigned irh=0; irh<rechitsECAL_.size();++irh)
04543         for ( unsigned isimrh=0;  isimrh < rechitSimIDs.size(); 
04544               isimrh++ )
04545           if(rechitSimIDs[isimrh] == rechitsECAL_[irh].detId())
04546             totalEcalE += (rechitsECAL_[irh].energy()*rechitSimFrac[isimrh]/100.0);
04547       cout << "For info, this particle deposits E=" << totalEcalE 
04548            << "(GeV) in the ECAL" << endl;
04549       
04550       for( unsigned icand=0; icand<candidates.size(); icand++ ) 
04551         {
04552           ITM it    = candSimMatchEcal[icand].begin();
04553           ITM itend = candSimMatchEcal[icand].end();
04554           for(;it!=itend;++it)
04555             if( i == it->second )
04556               out<<icand<<" "<<it->first<<"GeV "<<(*pfCandidates_)[icand]<<endl;          
04557         }//loop candidate
04558       cout << endl;
04559     }//loop particles  
04560   }//verbose
04561 
04562 }//mctruthmatching
04563 //_____________________________________________________________________________
04564 
04565 edm::InputTag 
04566 PFRootEventManager::stringToTag(const std::vector< std::string >& tagname) { 
04567 
04568   if ( tagname.size() == 1 ) 
04569     return edm::InputTag(tagname[0]);
04570 
04571   else if ( tagname.size() == 2 ) 
04572     return edm::InputTag(tagname[0], tagname[1]);
04573 
04574   else if ( tagname.size() == 3 ) 
04575     return tagname[2] == '*' ? 
04576       edm::InputTag(tagname[0], tagname[1]) :
04577       edm::InputTag(tagname[0], tagname[1], tagname[2]);
04578   else {
04579     cout << "Invalid tag name with " << tagname.size() << " strings "<< endl;
04580     return edm::InputTag();
04581   }
04582   
04583 }