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