CMS 3D CMS Logo

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