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