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