| 101 |
|
useTP_ = iConfig.getUntrackedParameter<bool>("useTP",false); |
| 102 |
|
useAssociator_ = iConfig.getUntrackedParameter<bool>("useAssociator",false); |
| 103 |
|
histoFileName_ = iConfig.getUntrackedParameter<std::string> ("histoFileName"); |
| 104 |
+ |
outputfilename_ = iConfig.getUntrackedParameter<string>("OutputFileName"); |
| 105 |
|
nTrkMin_ = iConfig.getUntrackedParameter<int>("nTrkMin"); |
| 106 |
|
nTrkMax_ = iConfig.getUntrackedParameter<int>("nTrkMax"); |
| 107 |
|
zsigncut_ = iConfig.getUntrackedParameter<double>("zsigncut"); |
| 108 |
|
ptcut_ = iConfig.getUntrackedParameter<double>("ptcut"); |
| 109 |
|
analyze_ = iConfig.getUntrackedParameter<bool>("analyze",false); |
| 110 |
+ |
saventuple_ = iConfig.getUntrackedParameter<bool>("saventuple",false); |
| 111 |
|
bsSrc = iConfig.getParameter< edm::InputTag >("beamSpot"); |
| 112 |
|
reqCluster_ = iConfig.getUntrackedParameter<bool>("reqCluster",false); |
| 113 |
|
|
| 117 |
|
datamode.push_back(0); |
| 118 |
|
datamode.push_back(1); |
| 119 |
|
} |
| 120 |
+ |
|
| 121 |
+ |
if(saventuple_) { |
| 122 |
+ |
file_ = TFile::Open(outputfilename_.c_str(),"RECREATE"); |
| 123 |
+ |
// pvtxtree_ analyzing the pvtxs ootb |
| 124 |
+ |
pvtxtree_ = new TTree("pvtxtree","pvtxtree"); |
| 125 |
+ |
if(!realData_) { |
| 126 |
+ |
pvtxtree_->Branch("nsimPV",&nsimPV_,"nsimPV/I"); |
| 127 |
+ |
pvtxtree_->Branch("nsimTrkPV",&nsimTrkPV_,"nsimTrkPV[nsimPV]/I"); |
| 128 |
+ |
pvtxtree_->Branch("simx",&simx_,"simx[nsimPV]/D"); |
| 129 |
+ |
pvtxtree_->Branch("simy",&simy_,"simy[nsimPV]/D"); |
| 130 |
+ |
pvtxtree_->Branch("simz",&simz_,"simz[nsimPV]/D"); |
| 131 |
+ |
pvtxtree_->Branch("simptsq",&simptsq_,"simptsq[nsimPV]/D"); |
| 132 |
+ |
pvtxtree_->Branch("simDeltaZ",&simDeltaZ_,"simDeltaZ[nsimPV]/D"); |
| 133 |
+ |
} |
| 134 |
+ |
} |
| 135 |
|
|
| 136 |
|
theFile = new TFile(histoFileName_.c_str(), "RECREATE"); |
| 137 |
|
theFile->mkdir("Summary"); |
| 138 |
|
theFile->cd(); |
| 139 |
|
|
| 140 |
+ |
|
| 141 |
|
// Book MC only plots |
| 142 |
|
if (!realData_) { |
| 143 |
|
h_gen = new PVEffHistograms(); |
| 273 |
|
return simpv; |
| 274 |
|
} |
| 275 |
|
|
| 276 |
+ |
|
| 277 |
+ |
// get sim pv from TrackingParticles/TrackingVertex |
| 278 |
+ |
std::vector<PVEffAnalyzer::simPrimaryVertex> PVEffAnalyzer::getSimPVs( const edm::Handle<TrackingVertexCollection> & tVC ) |
| 279 |
+ |
{ |
| 280 |
+ |
std::vector<simPrimaryVertex> simpv; |
| 281 |
+ |
//std::cout <<"number of vertices " << tVC->size() << std::endl; |
| 282 |
+ |
|
| 283 |
+ |
if(verbose_){std::cout << "getSimPVs TrackingVertexCollection " << std::endl;} |
| 284 |
+ |
|
| 285 |
+ |
for (TrackingVertexCollection::const_iterator v = tVC -> begin(); v != tVC -> end(); ++v) { |
| 286 |
+ |
|
| 287 |
+ |
if(verbose_){ |
| 288 |
+ |
std::cout << (v->eventId()).event() << v -> position() << v->g4Vertices().size() <<" " <<v->genVertices().size() << " t=" <<v->position().t()*1.e12 <<" ==0:" <<(v->position().t()>0) << std::endl; |
| 289 |
+ |
for( TrackingVertex::g4v_iterator gv=v->g4Vertices_begin(); gv!=v->g4Vertices_end(); gv++){ |
| 290 |
+ |
std::cout << *gv << std::endl; |
| 291 |
+ |
} |
| 292 |
+ |
std::cout << "----------" << std::endl; |
| 293 |
+ |
} |
| 294 |
+ |
|
| 295 |
+ |
// bool hasMotherVertex=false; |
| 296 |
+ |
if ((unsigned int) v->eventId().event()<simpv.size()) continue; |
| 297 |
+ |
//if (v->position().t()>0) { continue;} // skip secondary vertices (obsolete + doesn't work) |
| 298 |
+ |
if (fabs(v->position().z())>1000) continue; // skip funny junk vertices |
| 299 |
+ |
|
| 300 |
+ |
// could be a new vertex, check all primaries found so far to avoid multiple entries |
| 301 |
+ |
const double mm=1.0; // for tracking vertices |
| 302 |
+ |
simPrimaryVertex sv(v->position().x()*mm,v->position().y()*mm,v->position().z()*mm); |
| 303 |
+ |
|
| 304 |
+ |
for (TrackingParticleRefVector::iterator iTrack = v->daughterTracks_begin(); iTrack != v->daughterTracks_end(); ++iTrack){ |
| 305 |
+ |
//cout <<((**iTrack).eventId()).event() << " "; // an iterator of Refs, dereference twice. Cool eyh? |
| 306 |
+ |
sv.eventId=(**iTrack).eventId(); |
| 307 |
+ |
sv.ptsq += pow((**iTrack).pt(),2); |
| 308 |
+ |
sv.nGenTrk++; |
| 309 |
+ |
} |
| 310 |
+ |
|
| 311 |
+ |
simPrimaryVertex *vp=NULL; // will become non-NULL if a vertex is found and then point to it |
| 312 |
+ |
for(std::vector<simPrimaryVertex>::iterator v0=simpv.begin(); |
| 313 |
+ |
v0!=simpv.end(); v0++){ |
| 314 |
+ |
if( (sv.eventId==v0->eventId) && (fabs(sv.x-v0->x)<1e-5) && (fabs(sv.y-v0->y)<1e-5) && (fabs(sv.z-v0->z)<1e-5)){ |
| 315 |
+ |
vp=&(*v0); |
| 316 |
+ |
break; |
| 317 |
+ |
} |
| 318 |
+ |
} |
| 319 |
+ |
if(!vp){ |
| 320 |
+ |
// this is a new vertex |
| 321 |
+ |
if(verbose_){std::cout << "this is a new vertex " << sv.eventId.event() << " " << sv.x << " " << sv.y << " " << sv.z <<std::endl;} |
| 322 |
+ |
simpv.push_back(sv); |
| 323 |
+ |
vp=&simpv.back(); |
| 324 |
+ |
}else{ |
| 325 |
+ |
if(verbose_){std::cout << "this is not a new vertex" << sv.x << " " << sv.y << " " << sv.z <<std::endl;} |
| 326 |
+ |
} |
| 327 |
+ |
|
| 328 |
+ |
|
| 329 |
+ |
// Loop over daughter track(s) |
| 330 |
+ |
if(verbose_){ |
| 331 |
+ |
for (TrackingVertex::tp_iterator iTP = v -> daughterTracks_begin(); iTP != v -> daughterTracks_end(); ++iTP) { |
| 332 |
+ |
std::cout << " Daughter momentum: " << (*(*iTP)).momentum(); |
| 333 |
+ |
std::cout << " Daughter type " << (*(*iTP)).pdgId(); |
| 334 |
+ |
std::cout << std::endl; |
| 335 |
+ |
} |
| 336 |
+ |
} |
| 337 |
+ |
} |
| 338 |
+ |
if(verbose_){ |
| 339 |
+ |
cout << "------- PVEffAnalyzer simPVs from TrackingVertices -------" << endl; |
| 340 |
+ |
for(std::vector<simPrimaryVertex>::iterator v0=simpv.begin(); |
| 341 |
+ |
v0!=simpv.end(); v0++){ |
| 342 |
+ |
cout << "z=" << v0->z << " event=" << v0->eventId.event() << endl; |
| 343 |
+ |
} |
| 344 |
+ |
cout << "-----------------------------------------------" << endl; |
| 345 |
+ |
} |
| 346 |
+ |
return simpv; |
| 347 |
+ |
} |
| 348 |
+ |
|
| 349 |
+ |
|
| 350 |
+ |
|
| 351 |
+ |
|
| 352 |
|
// ------------ method called to for each event ------------ |
| 353 |
|
void |
| 354 |
|
PVEffAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) |
| 437 |
|
|
| 438 |
|
// ======= Analyze MC efficiency |
| 439 |
|
if (!realData_) { |
| 440 |
< |
bool MC=false; |
| 440 |
> |
|
| 441 |
> |
bool MC=false; // the flag to see if the MC level information is complete |
| 442 |
> |
|
| 443 |
|
Handle<HepMCProduct> evtMC; |
| 444 |
|
iEvent.getByLabel("generator",evtMC); |
| 445 |
< |
if (!evtMC.isValid()) { |
| 446 |
< |
MC=false; |
| 447 |
< |
edm::LogInfo("Debug") << "[PVEffAnalyzer] no HepMCProduct found"<< endl; |
| 448 |
< |
} else { |
| 449 |
< |
edm::LogInfo("Debug") << "[PVEffAnalyzer] generator HepMCProduct found"<< endl; |
| 445 |
> |
|
| 446 |
> |
edm::Handle<TrackingParticleCollection> TPCollectionH ; |
| 447 |
> |
bool gotTP=iEvent.getByLabel("mergedtruth","MergedTrackTruth",TPCollectionH); |
| 448 |
> |
|
| 449 |
> |
edm::Handle<TrackingVertexCollection> TVCollectionH ; |
| 450 |
> |
bool gotTV=iEvent.getByLabel("mergedtruth","MergedTrackTruth",TVCollectionH); |
| 451 |
> |
|
| 452 |
> |
std::vector<simPrimaryVertex> simpv; |
| 453 |
> |
|
| 454 |
> |
if(gotTV){ |
| 455 |
|
MC=true; |
| 456 |
+ |
if(verbose_){ |
| 457 |
+ |
cout << "Found Tracking Vertices " << endl; |
| 458 |
+ |
} |
| 459 |
+ |
simpv=getSimPVs(TVCollectionH); |
| 460 |
|
} |
| 461 |
+ |
else if (evtMC.isValid()) { |
| 462 |
+ |
MC=true; |
| 463 |
+ |
edm::LogInfo("Debug") << "[PVEffAnalyzer] generator HepMCProduct found"<< endl; |
| 464 |
+ |
} else { |
| 465 |
+ |
edm::LogInfo("Debug") << "[PVEffAnalyzer] no HepMCProduct found"<< endl; |
| 466 |
+ |
MC=false; |
| 467 |
+ |
} |
| 468 |
+ |
|
| 469 |
+ |
|
| 470 |
|
if(MC){ |
| 471 |
|
TString suffix = "_mct"; |
| 472 |
|
// make a list of primary vertices: |
| 359 |
– |
std::vector<simPrimaryVertex> simpv; |
| 360 |
– |
simpv=getSimPVs(evtMC,"",ptcut_); |
| 473 |
|
h_gen->Fill1d("nsimPV", simpv.size()); |
| 474 |
|
|
| 475 |
|
int isimpv = 0; |
| 476 |
< |
|
| 477 |
< |
for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin(); |
| 478 |
< |
vsim!=simpv.end(); vsim++, isimpv++){ |
| 479 |
< |
//nsimTrkPV_[isimpv] =vsim->nGenTrk; |
| 480 |
< |
//simx_[isimpv] = vsim->x; |
| 481 |
< |
//simy_[isimpv] = vsim->y; |
| 482 |
< |
//simz_[isimpv] = vsim->y; |
| 483 |
< |
//simptsq_[isimpv] = vsim->ptsq; |
| 484 |
< |
if(verbose_ && simpv.size() > 1 ) |
| 485 |
< |
std::cout<<"Simulated Vertex # " << isimpv << ": ptsq = " << vsim->ptsq<<std::endl; |
| 476 |
> |
nsimPV_ = simpv.size(); |
| 477 |
> |
|
| 478 |
> |
for(int isimpv = 0; isimpv<simpv.size();isimpv++) { |
| 479 |
> |
simPrimaryVertex vsim = simpv.at(isimpv); |
| 480 |
> |
nsimTrkPV_[isimpv] =vsim.nGenTrk; |
| 481 |
> |
simx_[isimpv] = vsim.x; |
| 482 |
> |
simy_[isimpv] = vsim.y; |
| 483 |
> |
simz_[isimpv] = vsim.z; |
| 484 |
> |
simptsq_[isimpv] = vsim.ptsq; |
| 485 |
> |
if(simpv.size() > 1 ) |
| 486 |
> |
if(verbose_) |
| 487 |
> |
std::cout<<"Simulated Vertex # " << isimpv << ": ptsq = " << vsim.ptsq<<std::endl; |
| 488 |
> |
if(isimpv >=1 ) |
| 489 |
> |
simDeltaZ_[isimpv-1] = vsim.z-simz_[0]; |
| 490 |
|
} |
| 491 |
|
|
| 492 |
|
// Just analyze the first vertex in the collection |
| 630 |
|
} |
| 631 |
|
} |
| 632 |
|
// == End of Efficiency Analyzing |
| 633 |
+ |
// == Save Ntuple |
| 634 |
+ |
if(saventuple_) { |
| 635 |
+ |
pvtxtree_->Fill(); |
| 636 |
+ |
} |
| 637 |
|
|
| 638 |
|
} |
| 639 |
|
|
| 648 |
|
void |
| 649 |
|
PVEffAnalyzer::endJob() { |
| 650 |
|
edm::LogInfo("Analysis") << "[endJob] Analyzing PV info" << endl; |
| 651 |
< |
|
| 651 |
> |
|
| 652 |
> |
if(saventuple_) { |
| 653 |
> |
file_->cd(); |
| 654 |
> |
pvtxtree_->Write(); |
| 655 |
> |
file_->Close(); |
| 656 |
> |
} |
| 657 |
> |
|
| 658 |
> |
|
| 659 |
|
for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) { |
| 660 |
|
string suffix; |
| 661 |
|
edm::LogInfo("Debug")<<"datamode = "<< *it<<endl; |
| 668 |
|
if (analyze_) { |
| 669 |
|
MakeEff(h_summary->ReadHisto1D(TString("eff_numer_ntrack"+suffix)), h_summary->ReadHisto1D(TString("eff_denom_ntrack"+suffix)), h_summary->ReadHisto1D(TString("eff_ntrack"+suffix)), false, 1); |
| 670 |
|
MakeEff(h_summary->ReadHisto1D(TString("fakerate_numer_ntrack"+suffix)), h_summary->ReadHisto1D(TString("fakerate_denom_ntrack"+suffix)), h_summary->ReadHisto1D(TString("fakerate_ntrack"+suffix)), false, 1); |
| 544 |
– |
|
| 671 |
|
} |
| 672 |
|
} |
| 673 |
|
} |
| 882 |
|
else |
| 883 |
|
return (vtx.ndof()+3.0)/2.0/double(vtx.tracksSize()); |
| 884 |
|
} |
| 885 |
+ |
|
| 886 |
+ |
|
| 887 |
+ |
void PVEffAnalyzer::SetVarToZero() { |
| 888 |
+ |
nsimPV_ = 0; |
| 889 |
+ |
//simpv |
| 890 |
+ |
for(int i = 0; i < nMaxPVs_; i++) { |
| 891 |
+ |
nsimTrkPV_[i] = 0; |
| 892 |
+ |
simx_[i] = 0; |
| 893 |
+ |
simy_[i] = 0; |
| 894 |
+ |
simz_[i] = 0; |
| 895 |
+ |
simptsq_[i] = 0; |
| 896 |
+ |
if(i>=1) |
| 897 |
+ |
simDeltaZ_[i-1] = 0; |
| 898 |
+ |
} |
| 899 |
+ |
} |
