CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Classes | Public Member Functions | Private Member Functions | Private Attributes
PrimaryVertexAnalyzer Class Reference

#include <Validation/RecoVertex/src/PrimaryVertexAnalyzer.cc>

Inheritance diagram for PrimaryVertexAnalyzer:
edm::EDAnalyzer edm::EDConsumerBase

Classes

class  simPrimaryVertex
 

Public Member Functions

virtual void analyze (const edm::Event &, const edm::EventSetup &)
 
virtual void beginJob ()
 
virtual void endJob ()
 
 PrimaryVertexAnalyzer (const edm::ParameterSet &)
 
 ~PrimaryVertexAnalyzer ()
 
- Public Member Functions inherited from edm::EDAnalyzer
void callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)
 
 EDAnalyzer ()
 
ModuleDescription const & moduleDescription () const
 
std::string workerType () const
 
virtual ~EDAnalyzer ()
 
- Public Member Functions inherited from edm::EDConsumerBase
 EDConsumerBase ()
 
ProductHolderIndexAndSkipBit indexFrom (EDGetToken, BranchType, TypeID const &) const
 
void itemsMayGet (BranchType, std::vector< ProductHolderIndexAndSkipBit > &) const
 
void itemsToGet (BranchType, std::vector< ProductHolderIndexAndSkipBit > &) const
 
std::vector
< ProductHolderIndexAndSkipBit >
const & 
itemsToGetFromEvent () const
 
void labelsForToken (EDGetToken iToken, Labels &oLabels) const
 
void modulesDependentUpon (const std::string &iProcessName, std::vector< const char * > &oModuleLabels) const
 
bool registeredToConsume (ProductHolderIndex, bool, BranchType) const
 
bool registeredToConsumeMany (TypeID const &, BranchType) const
 
void updateLookup (BranchType iBranchType, ProductHolderIndexHelper const &)
 
virtual ~EDConsumerBase ()
 

Private Member Functions

std::vector< simPrimaryVertexgetSimPVs (const edm::Handle< edm::HepMCProduct > &evtMC, const std::string &suffix="")
 
std::vector< simPrimaryVertexgetSimPVs (const edm::Handle< edm::HepMCProduct > &evt, const edm::Handle< edm::SimVertexContainer > &simVtxs, const edm::Handle< edm::SimTrackContainer > &simTrks)
 
bool isCharged (const HepMC::GenParticle *p)
 
bool isFinalstateParticle (const HepMC::GenParticle *p)
 
bool isResonance (const HepMC::GenParticle *p)
 
bool matchVertex (const simPrimaryVertex &vsim, const reco::Vertex &vrec)
 
void printRecVtxs (const edm::Handle< reco::VertexCollection > &recVtxs)
 
void printSimTrks (const edm::Handle< edm::SimTrackContainer > &simVtrks)
 
void printSimVtxs (const edm::Handle< edm::SimVertexContainer > &simVtxs)
 

Private Attributes

edm::EDGetTokenT
< edm::HepMCProduct
edmHepMCProductToken_
 
edm::EDGetTokenT
< edm::SimTrackContainer
edmSimTrackContainerToken_
 
edm::EDGetTokenT
< edm::SimVertexContainer
edmSimVertexContainerToken_
 
std::map< std::string, TH1 * > h
 
std::map< std::string,
TDirectory * > 
hdir
 
std::string outputFile_
 
edm::ESHandle< ParticleDataTablepdt
 
edm::EDGetTokenT
< reco::TrackCollection
recoTrackCollectionToken_
 
std::vector< edm::EDGetTokenT
< reco::VertexCollection > > 
recoVertexCollectionTokens_
 
TFile * rootFile_
 
double simUnit_
 
std::vector< std::string > suffixSample_
 
bool verbose_
 

Additional Inherited Members

- Public Types inherited from edm::EDAnalyzer
typedef EDAnalyzer ModuleType
 
- Static Public Member Functions inherited from edm::EDAnalyzer
static const std::string & baseType ()
 
static void fillDescriptions (ConfigurationDescriptions &descriptions)
 
static void prevalidate (ConfigurationDescriptions &)
 
- Protected Member Functions inherited from edm::EDConsumerBase
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType > consumes (edm::InputTag const &tag)
 
EDGetToken consumes (const TypeToGet &id, edm::InputTag const &tag)
 
template<BranchType B>
EDGetToken consumes (TypeToGet const &id, edm::InputTag const &tag)
 
ConsumesCollector consumesCollector ()
 Use a ConsumesCollector to gather consumes information from helper functions. More...
 
template<typename ProductType , BranchType B = InEvent>
void consumesMany ()
 
void consumesMany (const TypeToGet &id)
 
template<BranchType B>
void consumesMany (const TypeToGet &id)
 
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType > mayConsume (edm::InputTag const &tag)
 
EDGetToken mayConsume (const TypeToGet &id, edm::InputTag const &tag)
 
template<BranchType B>
EDGetToken mayConsume (const TypeToGet &id, edm::InputTag const &tag)
 

Detailed Description

Description: simple primary vertex analyzer

Implementation: <Notes on="" implementation>="">

Definition at line 54 of file PrimaryVertexAnalyzer.h.

Constructor & Destructor Documentation

PrimaryVertexAnalyzer::PrimaryVertexAnalyzer ( const edm::ParameterSet iConfig)
explicit

Definition at line 45 of file PrimaryVertexAnalyzer.cc.

References edm::getReleaseVersion(), edm::ParameterSet::getUntrackedParameter(), outputFile_, recoVertexCollectionTokens_, rootFile_, simUnit_, AlCaHLTBitMon_QueryRunRegistry::string, and suffixSample_.

46  : verbose_( iConfig.getUntrackedParameter<bool>( "verbose", false ) )
47  , simUnit_( 1.0 ) // starting with CMSSW_1_2_x ??
48  , recoTrackCollectionToken_( consumes< reco::TrackCollection >( edm::InputTag( iConfig.getUntrackedParameter<std::string>( "recoTrackProducer" ) ) ) )
49  , edmSimVertexContainerToken_( consumes< edm::SimVertexContainer>( iConfig.getParameter<edm::InputTag>( "simG4" ) ) )
50  , edmSimTrackContainerToken_( consumes< edm::SimTrackContainer >( iConfig.getParameter<edm::InputTag>( "simG4" ) ) )
51  , edmHepMCProductToken_( consumes< edm::HepMCProduct >( edm::InputTag( std::string( "generator" ) ) ) )
52 {
53  //now do what ever initialization is needed
54  // open output file to store histograms
55  outputFile_ = iConfig.getUntrackedParameter<std::string>("outputFile");
56  TString tversion(edm::getReleaseVersion());
57  tversion = tversion.Remove(0,1);
58  tversion = tversion.Remove(tversion.Length()-1,tversion.Length());
59  outputFile_ = std::string(tversion)+"_"+outputFile_;
60 
61  std::vector<std::string> vtxSample = iConfig.getUntrackedParameter<std::vector< std::string > >("vtxSample");
62  recoVertexCollectionTokens_.reserve( vtxSample.size() );
63  suffixSample_.reserve( vtxSample.size() );
64  auto vtxSampleBegin = vtxSample.begin();
65  auto vtxSampleEnd = vtxSample.end();
66  for( auto isample = vtxSampleBegin; isample != vtxSampleEnd; ++isample ) {
67  recoVertexCollectionTokens_.push_back( consumes< reco::VertexCollection >( edm::InputTag( *isample ) ) );
68  if ( *isample == "offlinePrimaryVertices" ) suffixSample_.push_back("AVF");
69  if ( *isample == "offlinePrimaryVerticesWithBS" ) suffixSample_.push_back("wBS");
70  }
71  if ( suffixSample_.size() == 0 ) throw cms::Exception("NoVertexSamples") << " no known vertex samples given";
72 
73  rootFile_ = new TFile(outputFile_.c_str(),"RECREATE");
74  if ( (edm::getReleaseVersion()).find("CMSSW_1_1_",0)!=std::string::npos){
75  simUnit_=0.1; // for use in CMSSW_1_1_1 tutorial
76  }
77 }
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
edm::EDGetTokenT< edm::SimTrackContainer > edmSimTrackContainerToken_
std::string getReleaseVersion()
edm::EDGetTokenT< edm::HepMCProduct > edmHepMCProductToken_
std::vector< edm::EDGetTokenT< reco::VertexCollection > > recoVertexCollectionTokens_
edm::EDGetTokenT< reco::TrackCollection > recoTrackCollectionToken_
std::vector< std::string > suffixSample_
edm::EDGetTokenT< edm::SimVertexContainer > edmSimVertexContainerToken_
PrimaryVertexAnalyzer::~PrimaryVertexAnalyzer ( )

Definition at line 80 of file PrimaryVertexAnalyzer.cc.

References rootFile_.

81 {
82  // do anything here that needs to be done at desctruction time
83  // (e.g. close files, deallocate resources etc.)
84  delete rootFile_;
85  /*for(std::map<std::string,TH1*>::const_iterator hist=h.begin(); hist!=h.end(); hist++){
86  delete hist->second;
87  }*/
88 
89 }

Member Function Documentation

void PrimaryVertexAnalyzer::analyze ( const edm::Event iEvent,
const edm::EventSetup iSetup 
)
virtual

Implements edm::EDAnalyzer.

Definition at line 410 of file PrimaryVertexAnalyzer.cc.

References ChiSquaredProbability(), gather_cfg::cout, data, edmHepMCProductToken_, edmSimTrackContainerToken_, edmSimVertexContainerToken_, edm::Event::getByToken(), edm::EventSetup::getData(), getSimPVs(), h, hdir, i, cmsHarvester::index, edm::detail::isnan(), edm::HandleBase::isValid(), j, m, matchVertex(), NULL, pdt, printRecVtxs(), recoTrackCollectionToken_, recoVertexCollectionTokens_, simUnit_, AlCaHLTBitMon_QueryRunRegistry::string, suffixSample_, edmStreamStallGrapher::t, findQualityFiles::v, verbose_, w, and detailsBasic3DVector::z.

411 {
412 
414  iEvent.getByToken( recoTrackCollectionToken_, recTrks );
415 
417  iEvent.getByToken( edmSimVertexContainerToken_, simVtxs );
418 
420  iEvent.getByToken( edmSimTrackContainerToken_, simTrks );
421 
422  for (int ivtxSample=0; ivtxSample!= (int)recoVertexCollectionTokens_.size(); ++ivtxSample) {
423  std::string isuffix = suffixSample_[ivtxSample];
424 
426  iEvent.getByToken( recoVertexCollectionTokens_[ ivtxSample ], recVtxs );
427 
428  try{
429  iSetup.getData(pdt);
430  }catch(const edm::Exception&){
431  std::cout << "Some problem occurred with the particle data table. This may not work !." <<std::endl;
432  }
433 
434  bool MC=false;
436  iEvent.getByToken( edmHepMCProductToken_, evtMC );
437  if (!evtMC.isValid()) {
438  MC=false;
439  if(verbose_){
440  std::cout << "no HepMCProduct found"<< std::endl;
441  }
442  } else {
443  if(verbose_){
444  std::cout << "generator HepMCProduct found"<< std::endl;
445  }
446  MC=true;
447  }
448 
449 
450  if(verbose_){
451  //evtMC->GetEvent()->print();
452  printRecVtxs(recVtxs);
453  //printSimVtxs(simVtxs);
454  //printSimTrks(simTrks);
455  }
456 
457  if(MC){
458  // make a list of primary vertices:
459  std::vector<simPrimaryVertex> simpv;
460  //simpv=getSimPVs(evtMC, simVtxs, simTrks);
461  simpv=getSimPVs(evtMC,isuffix);
462 
463 
464  // vertex matching and efficiency bookkeeping
465  h["nsimvtx"+isuffix]->Fill(simpv.size());
466  int nsimtrk=0;
467  for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
468  vsim!=simpv.end(); vsim++){
469 
470  hdir[isuffix]->cd();
471 
472  h["nbsimtksinvtx"+isuffix]->Fill(vsim->nGenTrk);
473  h["xsim"+isuffix]->Fill(vsim->x*simUnit_);
474  h["ysim"+isuffix]->Fill(vsim->y*simUnit_);
475  h["zsim"+isuffix]->Fill(vsim->z*simUnit_);
476  nsimtrk+=vsim->nGenTrk;
477  // look for a matching reconstructed vertex
478  vsim->recVtx=NULL;
479  for(reco::VertexCollection::const_iterator vrec=recVtxs->begin();
480  vrec!=recVtxs->end(); ++vrec){
481  if ( matchVertex(*vsim,*vrec) ){
482  // if the matching critera are fulfilled, accept the rec-vertex that is closest in z
483  if( ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z)>fabs(vrec->z()-vsim->z)))
484  || (!vsim->recVtx) )
485  {
486  vsim->recVtx=&(*vrec);
487  }
488  }
489  }
490  h["nsimtrk"+isuffix]->Fill(float(nsimtrk));
491 
492 
493  // histogram properties of matched vertices
494  if (vsim->recVtx){
495 
496  if(verbose_){std::cout <<"primary matched " << vsim->x << " " << vsim->y << " " << vsim->z << std:: endl;}
497 
498  h["resx"+isuffix]->Fill( vsim->recVtx->x()-vsim->x*simUnit_ );
499  h["resy"+isuffix]->Fill( vsim->recVtx->y()-vsim->y*simUnit_ );
500  h["resz"+isuffix]->Fill( vsim->recVtx->z()-vsim->z*simUnit_ );
501  h["pullx"+isuffix]->Fill( (vsim->recVtx->x()-vsim->x*simUnit_)/vsim->recVtx->xError() );
502  h["pully"+isuffix]->Fill( (vsim->recVtx->y()-vsim->y*simUnit_)/vsim->recVtx->yError() );
503  h["pullz"+isuffix]->Fill( (vsim->recVtx->z()-vsim->z*simUnit_)/vsim->recVtx->zError() );
504  h["eff"+isuffix]->Fill( 1.);
505  if(simpv.size()==1){
506  if (vsim->recVtx==&(*recVtxs->begin())){
507  h["efftag"+isuffix]->Fill( 1.);
508  }else{
509  h["efftag"+isuffix]->Fill( 0.);
510  }
511  }
512  h["effvseta"+isuffix]->Fill(vsim->ptot.pseudoRapidity(),1.);
513  h["effvsptsq"+isuffix]->Fill(vsim->ptsq,1.);
514  h["effvsntrk"+isuffix]->Fill(vsim->nGenTrk,1.);
515  h["effvsnrectrk"+isuffix]->Fill(recTrks->size(),1.);
516  h["effvsz"+isuffix]->Fill(vsim->z*simUnit_,1.);
517 
518  }else{ // no rec vertex found for this simvertex
519 
520  if(verbose_){std::cout <<"primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << std::endl;}
521 
522  h["eff"+isuffix]->Fill( 0.);
523  if(simpv.size()==1){ h["efftag"+isuffix]->Fill( 0.); }
524  h["effvseta"+isuffix]->Fill(vsim->ptot.pseudoRapidity(),0.);
525  h["effvsptsq"+isuffix]->Fill(vsim->ptsq,0.);
526  h["effvsntrk"+isuffix]->Fill(float(vsim->nGenTrk),0.);
527  h["effvsnrectrk"+isuffix]->Fill(recTrks->size(),0.);
528  h["effvsz"+isuffix]->Fill(vsim->z*simUnit_,0.);
529  } // no recvertex for this simvertex
530  }//found MC event
531  // end of sim/rec matching
532 
533 
534  // purity of event vertex tags
535  if (recVtxs->size()>0 && simpv.size()>0){
536  Double_t dz=(*recVtxs->begin()).z() - (*simpv.begin()).z*simUnit_;
537  h["zdistancetag"+isuffix]->Fill(dz);
538  if( fabs(dz)<0.0500){
539  h["puritytag"+isuffix]->Fill(1.);
540  }else{
541  // bad tag: the true primary was more than 500 um away from the tagged primary
542  h["puritytag"+isuffix]->Fill(0.);
543  }
544  }
545 
546 
547  }// MC event
548 
549  // test track links, use reconstructed vertices
550 
551  h["nrecvtx"+isuffix]->Fill(recVtxs->size());
552  h["nrectrk"+isuffix]->Fill(recTrks->size());
553  h["nbvtx"+isuffix]->Fill(recVtxs->size()*1.);
554  if((recVtxs->size()==0)||recVtxs->begin()->isFake()) {h["nrectrk0vtx"+isuffix]->Fill(recTrks->size());}
555 
556  for(reco::VertexCollection::const_iterator v=recVtxs->begin();
557  v!=recVtxs->end(); ++v){
558 
559  if(v->isFake()) continue;
560 
561  try {
562  for(reco::Vertex::trackRef_iterator t = v->tracks_begin();
563  t!=v->tracks_end(); t++) {
564  // illegal charge
565  if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
566  h["tklinks"+isuffix]->Fill(0.);
567  }
568  else {
569  h["tklinks"+isuffix]->Fill(1.);
570  }
571  }
572  }
573  catch (...) {
574  // exception thrown when trying to use linked track
575  h["tklinks"+isuffix]->Fill(0.);
576  }
577 
578  h["nbtksinvtx"+isuffix]->Fill(v->tracksSize());
579  h["vtxchi2"+isuffix]->Fill(v->chi2());
580  h["vtxndf"+isuffix]->Fill(v->ndof());
581  h["vtxprob"+isuffix]->Fill(ChiSquaredProbability(v->chi2() ,v->ndof()));
582  h["xrec"+isuffix]->Fill(v->position().x());
583  h["yrec"+isuffix]->Fill(v->position().y());
584  h["zrec"+isuffix]->Fill(v->position().z());
585  if (v==recVtxs->begin()){
586  h["xrectag"+isuffix]->Fill(v->position().x());
587  h["yrectag"+isuffix]->Fill(v->position().y());
588  h["zrectag"+isuffix]->Fill(v->position().z());
589  }
590 
591  bool problem = false;
592  h["nans"+isuffix]->Fill(1.,isnan(v->position().x())*1.);
593  h["nans"+isuffix]->Fill(2.,isnan(v->position().y())*1.);
594  h["nans"+isuffix]->Fill(3.,isnan(v->position().z())*1.);
595 
596  int index = 3;
597  for (int i = 0; i != 3; i++) {
598  for (int j = i; j != 3; j++) {
599  index++;
600  h["nans"+isuffix]->Fill(index*1., isnan(v->covariance(i, j))*1.);
601  if (isnan(v->covariance(i, j))) problem = true;
602  // in addition, diagonal element must be positive
603  if (j == i && v->covariance(i, j) < 0) {
604  h["nans"+isuffix]->Fill(index*1., 1.);
605  problem = true;
606  }
607  }
608  }
609 
610  if (problem) {
611  // analyze track parameter covariance definiteness
612  double data[25];
613  try {
614  int itk = 0;
615  for(reco::Vertex::trackRef_iterator t = v->tracks_begin();
616  t!=v->tracks_end(); t++) {
617  std::cout << "Track " << itk++ << std::endl;
618  int i2 = 0;
619  for (int i = 0; i != 5; i++) {
620  for (int j = 0; j != 5; j++) {
621  data[i2] = (**t).covariance(i, j);
622  std::cout << std:: scientific << data[i2] << " ";
623  i2++;
624  }
625  std::cout << std::endl;
626  }
627  gsl_matrix_view m
628  = gsl_matrix_view_array (data, 5, 5);
629 
630  gsl_vector *eval = gsl_vector_alloc (5);
631  gsl_matrix *evec = gsl_matrix_alloc (5, 5);
632 
633  gsl_eigen_symmv_workspace * w =
634  gsl_eigen_symmv_alloc (5);
635 
636  gsl_eigen_symmv (&m.matrix, eval, evec, w);
637 
638  gsl_eigen_symmv_free (w);
639 
640  gsl_eigen_symmv_sort (eval, evec,
641  GSL_EIGEN_SORT_ABS_ASC);
642 
643  // print sorted eigenvalues
644  {
645  int i;
646  for (i = 0; i < 5; i++) {
647  double eval_i
648  = gsl_vector_get (eval, i);
649  //gsl_vector_view evec_i
650  // = gsl_matrix_column (evec, i);
651 
652  printf ("eigenvalue = %g\n", eval_i);
653  // printf ("eigenvector = \n");
654  // gsl_vector_fprintf (stdout,
655  // &evec_i.vector, "%g");
656  }
657  }
658  }
659  }
660  catch (...) {
661  // exception thrown when trying to use linked track
662  break;
663  }
664  }// if (problem)
665  }
666  } // for vertex loop
667 }
int i
Definition: DBlmapReader.cc:9
void printRecVtxs(const edm::Handle< reco::VertexCollection > &recVtxs)
const double w
Definition: UKUtility.cc:23
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:434
edm::EDGetTokenT< edm::SimTrackContainer > edmSimTrackContainerToken_
std::map< std::string, TDirectory * > hdir
#define NULL
Definition: scimark2.h:8
float float float z
void getData(T &iHolder) const
Definition: EventSetup.h:67
edm::ESHandle< ParticleDataTable > pdt
bool isnan(float x)
Definition: math.h:13
int j
Definition: DBlmapReader.cc:9
float ChiSquaredProbability(double chiSquared, double nrDOF)
bool isValid() const
Definition: HandleBase.h:76
edm::EDGetTokenT< edm::HepMCProduct > edmHepMCProductToken_
std::map< std::string, TH1 * > h
std::vector< simPrimaryVertex > getSimPVs(const edm::Handle< edm::HepMCProduct > &evtMC, const std::string &suffix="")
std::vector< edm::EDGetTokenT< reco::VertexCollection > > recoVertexCollectionTokens_
edm::EDGetTokenT< reco::TrackCollection > recoTrackCollectionToken_
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
std::vector< TrackBaseRef >::const_iterator trackRef_iterator
The iteratator for the vector&lt;TrackRef&gt;
Definition: Vertex.h:37
tuple cout
Definition: gather_cfg.py:121
std::vector< std::string > suffixSample_
bool matchVertex(const simPrimaryVertex &vsim, const reco::Vertex &vrec)
edm::EDGetTokenT< edm::SimVertexContainer > edmSimVertexContainerToken_
void PrimaryVertexAnalyzer::beginJob ( void  )
virtual

Reimplemented from edm::EDAnalyzer.

Definition at line 96 of file PrimaryVertexAnalyzer.cc.

References gather_cfg::cout, h, hdir, rootFile_, simUnit_, and suffixSample_.

96  {
97  std::cout << " PrimaryVertexAnalyzer::beginJob conversion from sim units to rec units is " << simUnit_ << std::endl;
98 
99 
100  rootFile_->cd();
101  // release validation histograms used in DoCompare.C
102  for (std::vector<std::string>::iterator isuffix= suffixSample_.begin(); isuffix!=suffixSample_.end(); isuffix++) {
103 
104  hdir[*isuffix] = rootFile_->mkdir(TString(*isuffix));
105  hdir[*isuffix]->cd();
106  h["nbvtx"+ *isuffix] = new TH1F(TString("nbvtx"+ *isuffix),"Reconstructed Vertices in Event",20,-0.5,19.5);
107  h["nbtksinvtx"+ *isuffix] = new TH1F(TString("nbtksinvtx"+ *isuffix),"Reconstructed Tracks in Vertex",200,-0.5,199.5);
108  h["resx"+ *isuffix] = new TH1F(TString("resx"+ *isuffix),"Residual X",100,-0.04,0.04);
109  h["resy"+ *isuffix] = new TH1F(TString("resy"+ *isuffix),"Residual Y",100,-0.04,0.04);
110  h["resz"+ *isuffix] = new TH1F(TString("resz"+ *isuffix),"Residual Z",100,-0.1,0.1);
111  h["pullx"+ *isuffix] = new TH1F(TString("pullx"+ *isuffix),"Pull X",100,-25.,25.);
112  h["pully"+ *isuffix] = new TH1F(TString("pully"+ *isuffix),"Pull Y",100,-25.,25.);
113  h["pullz"+ *isuffix] = new TH1F(TString("pullz"+ *isuffix),"Pull Z",100,-25.,25.);
114  h["vtxchi2"+ *isuffix] = new TH1F(TString("vtxchi2"+ *isuffix),"#chi^{2}",100,0.,100.);
115  h["vtxndf"+ *isuffix] = new TH1F(TString("vtxndf"+ *isuffix),"ndof",100,0.,100.);
116  h["tklinks"+ *isuffix] = new TH1F(TString("tklinks"+ *isuffix),"Usable track links",2,-0.5,1.5);
117  h["nans"+ *isuffix] = new TH1F(TString("nans"+ *isuffix),"Illegal values for x,y,z,xx,xy,xz,yy,yz,zz",9,0.5,9.5);
118  // more histograms
119  h["vtxprob"+ *isuffix] = new TH1F(TString("vtxprob"+ *isuffix),"#chi^{2} probability",100,0.,1.);
120  h["eff"+ *isuffix] = new TH1F(TString("eff"+ *isuffix),"efficiency",2, -0.5, 1.5);
121  h["efftag"+ *isuffix] = new TH1F(TString("efftag"+ *isuffix),"efficiency tagged vertex",2, -0.5, 1.5);
122  h["effvseta"+ *isuffix] = new TProfile(TString("effvseta"+ *isuffix),"efficiency vs eta",20, -2.5, 2.5, 0, 1.);
123  h["effvsptsq"+ *isuffix] = new TProfile(TString("effvsptsq"+ *isuffix),"efficiency vs ptsq",20, 0., 10000., 0, 1.);
124  h["effvsntrk"+ *isuffix] = new TProfile(TString("effvsntrk"+ *isuffix),"efficiency vs # simtracks",200, 0., 200., 0, 1.);
125  h["effvsnrectrk"+ *isuffix] = new TProfile(TString("effvsnrectrk"+ *isuffix),"efficiency vs # rectracks",200, 0., 200., 0, 1.);
126  h["effvsz"+ *isuffix] = new TProfile(TString("effvsz"+ *isuffix),"efficiency vs z",40, -20., 20., 0, 1.);
127  h["nbsimtksinvtx"+ *isuffix]= new TH1F(TString("nbsimtksinvtx"+ *isuffix),"simulated tracks in vertex",100,-0.5,99.5);
128  h["xrec"+ *isuffix] = new TH1F(TString("xrec"+ *isuffix),"reconstructed x",100,-0.1,0.1);
129  h["yrec"+ *isuffix] = new TH1F(TString("yrec"+ *isuffix),"reconstructed y",100,-0.1,0.1);
130  h["zrec"+ *isuffix] = new TH1F(TString("zrec"+ *isuffix),"reconstructed z",100,-20.,20.);
131  h["xsim"+ *isuffix] = new TH1F(TString("xsim"+ *isuffix),"simulated x",100,-0.1,0.1);
132  h["ysim"+ *isuffix] = new TH1F(TString("ysim"+ *isuffix),"simulated y",100,-0.1,0.1);
133  h["zsim"+ *isuffix] = new TH1F(TString("zsim"+ *isuffix),"simulated z",100,-20.,20.);
134  h["nrecvtx"+ *isuffix] = new TH1F(TString("nrecvtx"+ *isuffix),"# of reconstructed vertices", 50, -0.5, 49.5);
135  h["nsimvtx"+ *isuffix] = new TH1F(TString("nsimvtx"+ *isuffix),"# of simulated vertices", 50, -0.5, 49.5);
136  h["nrectrk"+ *isuffix] = new TH1F(TString("nrectrk"+ *isuffix),"# of reconstructed tracks", 200, -0.5, 199.5);
137  h["nsimtrk"+ *isuffix] = new TH1F(TString("nsimtrk"+ *isuffix),"# of simulated tracks", 200, -0.5, 199.5);
138  h["xrectag"+ *isuffix] = new TH1F(TString("xrectag"+ *isuffix),"reconstructed x, signal vtx",100,-0.1,0.1);
139  h["yrectag"+ *isuffix] = new TH1F(TString("yrectag"+ *isuffix),"reconstructed y, signal vtx",100,-0.1,0.1);
140  h["zrectag"+ *isuffix] = new TH1F(TString("zrectag"+ *isuffix),"reconstructed z, signal vtx",100,-20.,20.);
141  h["rapidity"+ *isuffix] = new TH1F(TString("rapidity"+ *isuffix),"rapidity ",100,-10., 10.);
142  h["pt"+ *isuffix] = new TH1F(TString("pt"+ *isuffix),"pt ",100,0., 20.);
143  h["nrectrk0vtx"+ *isuffix] = new TH1F(TString("nrectrk0vtx"+ *isuffix),"# rec tracks no vertex ",200,0., 200.);
144  h["zdistancetag"+ *isuffix] = new TH1F(TString("zdistancetag"+ *isuffix),"z-distance between tagged and generated",100, -0.1, 0.1);
145  h["puritytag"+ *isuffix] = new TH1F(TString("puritytag"+ *isuffix),"purity of primary vertex tags",2, -0.5, 1.5);
146  }
147 
148 }
std::map< std::string, TDirectory * > hdir
std::map< std::string, TH1 * > h
tuple cout
Definition: gather_cfg.py:121
std::vector< std::string > suffixSample_
void PrimaryVertexAnalyzer::endJob ( void  )
virtual

Reimplemented from edm::EDAnalyzer.

Definition at line 151 of file PrimaryVertexAnalyzer.cc.

References h, hdir, estimatePileup::hist, and rootFile_.

151  {
152  rootFile_->cd();
153  // save all histograms created in beginJob()
154  for (std::map<std::string, TDirectory*>::const_iterator idir=hdir.begin(); idir!=hdir.end(); idir++){
155  idir->second->cd();
156  for(std::map<std::string,TH1*>::const_iterator hist=h.begin(); hist!=h.end(); hist++){
157  if (TString(hist->first).Contains(idir->first)) hist->second->Write();
158  }
159  }
160 }
std::map< std::string, TDirectory * > hdir
std::map< std::string, TH1 * > h
std::vector< PrimaryVertexAnalyzer::simPrimaryVertex > PrimaryVertexAnalyzer::getSimPVs ( const edm::Handle< edm::HepMCProduct > &  evtMC,
const std::string &  suffix = "" 
)
private

Definition at line 243 of file PrimaryVertexAnalyzer.cc.

References gather_cfg::cout, alignCSCRings::e, spr::find(), h, customizeTrackingMonitorSeedNumber::idx, isCharged(), isFinalstateParticle(), m, NULL, parents, createPayload::suffix, and verbose_.

Referenced by analyze().

246 {
247  std::vector<PrimaryVertexAnalyzer::simPrimaryVertex> simpv;
248  const HepMC::GenEvent* evt=evtMC->GetEvent();
249  if (evt) {
250  if(verbose_) std::cout << "process id " <<evt->signal_process_id()<<std::endl;
251  if(verbose_) std::cout <<"signal process vertex "<< ( evt->signal_process_vertex() ?
252  evt->signal_process_vertex()->barcode() : 0 ) <<std::endl;
253  if(verbose_) std::cout <<"number of vertices " << evt->vertices_size() << std::endl;
254 
255 
256  int idx=0;
257  for(HepMC::GenEvent::vertex_const_iterator vitr= evt->vertices_begin();
258  vitr != evt->vertices_end(); ++vitr )
259  { // loop for vertex ...
260  HepMC::FourVector pos = (*vitr)->position();
261  //HepLorentzVector pos = (*vitr)->position();
262 
263  bool hasMotherVertex=false;
264  if(verbose_) std::cout << "mothers of vertex[" << ++idx << "]: " << std::endl;
265  for ( HepMC::GenVertex::particle_iterator
266  mother = (*vitr)->particles_begin(HepMC::parents);
267  mother != (*vitr)->particles_end(HepMC::parents);
268  ++mother ) {
269  HepMC::GenVertex * mv=(*mother)->production_vertex();
270  if (mv) {
271  hasMotherVertex=true;
272  if(!verbose_) break; //if verbose_, print all particles of gen vertices
273  }
274  if(verbose_)std::cout << "\t";
275  if(verbose_)(*mother)->print();
276  }
277 
278  if(hasMotherVertex) {continue;}
279 
280  // could be a new vertex, check all primaries found so far to avoid multiple entries
281  const double mm=0.1;
282  simPrimaryVertex sv(pos.x()*mm,pos.y()*mm,pos.z()*mm);
283  simPrimaryVertex *vp=NULL; // will become non-NULL if a vertex is found and then point to it
284  for(std::vector<simPrimaryVertex>::iterator v0=simpv.begin();
285  v0!=simpv.end(); v0++){
286  if( (fabs(sv.x-v0->x)<1e-5) && (fabs(sv.y-v0->y)<1e-5) && (fabs(sv.z-v0->z)<1e-5)){
287  vp=&(*v0);
288  break;
289  }
290  }
291 
292  if(!vp){
293  // this is a new vertex
294  if(verbose_)std::cout << "this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << std::endl;
295  simpv.push_back(sv);
296  vp=&simpv.back();
297  }else{
298  if(verbose_)std::cout << "this is not new vertex " << sv.x << " " << sv.y << " " << sv.z << std::endl;
299  }
300  vp->genVertex.push_back((*vitr)->barcode());
301  // collect final state descendants
302  for ( HepMC::GenVertex::particle_iterator
303  daughter = (*vitr)->particles_begin(HepMC::descendants);
304  daughter != (*vitr)->particles_end(HepMC::descendants);
305  ++daughter ) {
306  if (isFinalstateParticle(*daughter)){
307  if ( find(vp->finalstateParticles.begin(), vp->finalstateParticles.end(),(*daughter)->barcode())
308  == vp->finalstateParticles.end()){
309  vp->finalstateParticles.push_back((*daughter)->barcode());
310  HepMC::FourVector m=(*daughter)->momentum();
311  // the next four lines used to be "vp->ptot+=m;" in the days of CLHEP::HepLorentzVector
312  // but adding FourVectors seems not to be foreseen
313  vp->ptot.setPx(vp->ptot.px()+m.px());
314  vp->ptot.setPy(vp->ptot.py()+m.py());
315  vp->ptot.setPz(vp->ptot.pz()+m.pz());
316  vp->ptot.setE(vp->ptot.e()+m.e());
317  vp->ptsq+=(m.perp())*(m.perp());
318  if ( (m.perp()>0.8) && (fabs(m.pseudoRapidity())<2.5) && isCharged( *daughter ) ){
319  vp->nGenTrk++;
320  }
321 
322  h["rapidity"+suffix]->Fill(m.pseudoRapidity());
323  h["pt"+suffix]->Fill(m.perp());
324  }
325  }
326  }
327  }
328  }
329  return simpv;
330 }
TPRegexp parents
Definition: eve_filter.cc:24
bool isCharged(const HepMC::GenParticle *p)
bool isFinalstateParticle(const HepMC::GenParticle *p)
#define NULL
Definition: scimark2.h:8
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:7
tuple idx
DEBUGGING if hasattr(process,&quot;trackMonIterativeTracking2012&quot;): print &quot;trackMonIterativeTracking2012 D...
std::map< std::string, TH1 * > h
tuple cout
Definition: gather_cfg.py:121
std::vector< PrimaryVertexAnalyzer::simPrimaryVertex > PrimaryVertexAnalyzer::getSimPVs ( const edm::Handle< edm::HepMCProduct > &  evt,
const edm::Handle< edm::SimVertexContainer > &  simVtxs,
const edm::Handle< edm::SimTrackContainer > &  simTrks 
)
private

Definition at line 334 of file PrimaryVertexAnalyzer.cc.

References configurableAnalysis::GenParticle, customizeTrackingMonitorSeedNumber::idx, isCharged(), isFinalstateParticle(), isResonance(), or, PrimaryVertexAnalyzer::simPrimaryVertex::simTrackIndex, edmStreamStallGrapher::t, and findQualityFiles::v.

338 {
339  // simvertices don't have enough information to decide,
340  // genVertices don't have the simulated coordinates ( with VtxSmeared they might)
341  // go through simtracks to get the link between simulated and generated vertices
342  std::vector<PrimaryVertexAnalyzer::simPrimaryVertex> simpv;
343  int idx=0;
344  for(edm::SimTrackContainer::const_iterator t=simTrks->begin();
345  t!=simTrks->end(); ++t){
346  if ( !(t->noVertex()) && !(t->type()==-99) ){
347  double ptsq=0;
348  bool primary=false; // something coming directly from the primary vertex
349  bool resonance=false; // resonance
350  bool track=false; // undecayed, charged particle
351  HepMC::GenParticle* gp=evtMC->GetEvent()->barcode_to_particle( (*t).genpartIndex() );
352  if (gp) {
353  HepMC::GenVertex * gv=gp->production_vertex();
354  if (gv) {
355  for ( HepMC::GenVertex::particle_iterator
356  daughter = gv->particles_begin(HepMC::descendants);
357  daughter != gv->particles_end(HepMC::descendants);
358  ++daughter ) {
359  if (isFinalstateParticle(*daughter)){
360  ptsq+=(*daughter)->momentum().perp()*(*daughter)->momentum().perp();
361  }
362  }
363  primary = ( gv->position().t()==0);
364  //resonance= ( gp->mother() && isResonance(gp->mother())); // in CLHEP/HepMC days
365  // no more mother pointer in the improved HepMC GenParticle
366  resonance= ( isResonance(*(gp->production_vertex()->particles_in_const_begin())));
367  if (gp->status()==1){
368  //track=((pdt->particle(gp->pdg_id()))->charge() != 0);
369  track=not isCharged(gp);
370  }
371  }
372  }
373 
374  const HepMC::FourVector & v=(*simVtxs)[t->vertIndex()].position();
375  //const HepLorentzVector & v=(*simVtxs)[t->vertIndex()].position();
376  if(primary or resonance){
377  {
378  // check all primaries found so far to avoid multiple entries
379  bool newVertex=true;
380  for(std::vector<simPrimaryVertex>::iterator v0=simpv.begin();
381  v0!=simpv.end(); v0++){
382  if( (fabs(v0->x-v.x())<0.001) && (fabs(v0->y-v.y())<0.001) && (fabs(v0->z-v.z())<0.001) ){
383  if (track) {
384  v0->simTrackIndex.push_back(idx);
385  if (ptsq>(*v0).ptsq){(*v0).ptsq=ptsq;}
386  }
387  newVertex=false;
388  }
389  }
390  if(newVertex && !resonance){
391  simPrimaryVertex anotherVertex(v.x(),v.y(),v.z());
392  if (track) anotherVertex.simTrackIndex.push_back(idx);
393  anotherVertex.ptsq=ptsq;
394  simpv.push_back(anotherVertex);
395  }
396  }//
397  }
398 
399  }// simtrack has vertex and valid type
400  idx++;
401  }//simTrack loop
402  return simpv;
403 }
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e< void, edm::EventIDconst &, edm::Timestampconst & > We also list in braces which AR_WATCH_USING_METHOD_ is used for those or
Definition: Activities.doc:12
bool isCharged(const HepMC::GenParticle *p)
bool isFinalstateParticle(const HepMC::GenParticle *p)
bool isResonance(const HepMC::GenParticle *p)
tuple idx
DEBUGGING if hasattr(process,&quot;trackMonIterativeTracking2012&quot;): print &quot;trackMonIterativeTracking2012 D...
bool PrimaryVertexAnalyzer::isCharged ( const HepMC::GenParticle *  p)
private

Definition at line 179 of file PrimaryVertexAnalyzer.cc.

References pdt.

Referenced by getSimPVs().

179  {
180  const ParticleData * part = pdt->particle( p->pdg_id() );
181  if (part){
182  return part->charge()!=0;
183  }else{
184  // the new/improved particle table doesn't know anti-particles
185  return pdt->particle( -p->pdg_id() )!=0;
186  }
187 }
edm::ESHandle< ParticleDataTable > pdt
HepPDT::ParticleData ParticleData
part
Definition: HCALResponse.h:20
bool PrimaryVertexAnalyzer::isFinalstateParticle ( const HepMC::GenParticle *  p)
private

Definition at line 175 of file PrimaryVertexAnalyzer.cc.

Referenced by getSimPVs().

175  {
176  return ( !p->end_vertex() && p->status()==1 );
177 }
bool PrimaryVertexAnalyzer::isResonance ( const HepMC::GenParticle *  p)
private

Definition at line 169 of file PrimaryVertexAnalyzer.cc.

References funct::abs(), gather_cfg::cout, alignCSCRings::e, pdt, and verbose_.

Referenced by getSimPVs().

169  {
170  double ctau=(pdt->particle( abs(p->pdg_id()) ))->lifetime();
171  if(verbose_) std::cout << "isResonance " << p->pdg_id() << " " << ctau << std::endl;
172  return ctau >0 && ctau <1e-6;
173 }
edm::ESHandle< ParticleDataTable > pdt
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
tuple cout
Definition: gather_cfg.py:121
bool PrimaryVertexAnalyzer::matchVertex ( const simPrimaryVertex vsim,
const reco::Vertex vrec 
)
private

Definition at line 164 of file PrimaryVertexAnalyzer.cc.

References simUnit_, PrimaryVertexAnalyzer::simPrimaryVertex::z, and reco::Vertex::z().

Referenced by analyze().

165  {
166  return (fabs(vsim.z*simUnit_-vrec.z())<0.0500); // =500um
167 }
double z() const
y coordinate
Definition: Vertex.h:112
void PrimaryVertexAnalyzer::printRecVtxs ( const edm::Handle< reco::VertexCollection > &  recVtxs)
private

Definition at line 189 of file PrimaryVertexAnalyzer.cc.

References gather_cfg::cout, and findQualityFiles::v.

Referenced by analyze().

189  {
190  int ivtx=0;
191  for(reco::VertexCollection::const_iterator v=recVtxs->begin();
192  v!=recVtxs->end(); ++v){
193  std::cout << "Recvtx "<< std::setw(3) << std::setfill(' ')<<ivtx++
194  << "#trk " << std::setw(3) << v->tracksSize()
195  << " chi2 " << std::setw(4) << v->chi2()
196  << " ndof " << std::setw(3) << v->ndof() << std::endl
197  << " x " << std::setw(8) <<std::fixed << std::setprecision(4) << v->x()
198  << " dx " << std::setw(8) << v->xError()<< std::endl
199  << " y " << std::setw(8) << v->y()
200  << " dy " << std::setw(8) << v->yError()<< std::endl
201  << " z " << std::setw(8) << v->z()
202  << " dz " << std::setw(8) << v->zError()
203  << std::endl;
204  }
205 }
tuple cout
Definition: gather_cfg.py:121
void PrimaryVertexAnalyzer::printSimTrks ( const edm::Handle< edm::SimTrackContainer > &  simVtrks)
private

Definition at line 224 of file PrimaryVertexAnalyzer.cc.

References gather_cfg::cout, i, and edmStreamStallGrapher::t.

224  {
225  std::cout << " simTrks type, (momentum), vertIndex, genpartIndex" << std::endl;
226  int i=1;
227  for(edm::SimTrackContainer::const_iterator t=simTrks->begin();
228  t!=simTrks->end(); ++t){
229  //HepMC::GenParticle* gp=evtMC->GetEvent()->particle( (*t).genpartIndex() );
230  std::cout << i++ << ")"
231  << (*t)
232  << " index="
233  << (*t).genpartIndex();
234  //if (gp) {
235  // HepMC::GenVertex *gv=gp->production_vertex();
236  // std::cout << " genvertex =" << (*gv);
237  //}
238  std::cout << std::endl;
239  }
240 }
int i
Definition: DBlmapReader.cc:9
tuple cout
Definition: gather_cfg.py:121
void PrimaryVertexAnalyzer::printSimVtxs ( const edm::Handle< edm::SimVertexContainer > &  simVtxs)
private

Definition at line 208 of file PrimaryVertexAnalyzer.cc.

References gather_cfg::cout, i, and simUnit_.

208  {
209  int i=0;
210  for(edm::SimVertexContainer::const_iterator vsim=simVtxs->begin();
211  vsim!=simVtxs->end(); ++vsim){
212  std::cout << i++ << ")" << std::scientific
213  << " evtid=" << vsim->eventId().event()
214  << " sim x=" << vsim->position().x()*simUnit_
215  << " sim y=" << vsim->position().y()*simUnit_
216  << " sim z=" << vsim->position().z()*simUnit_
217  << " sim t=" << vsim->position().t()
218  << " parent=" << vsim->parentIndex()
219  << std::endl;
220  }
221 }
int i
Definition: DBlmapReader.cc:9
tuple cout
Definition: gather_cfg.py:121

Member Data Documentation

edm::EDGetTokenT< edm::HepMCProduct > PrimaryVertexAnalyzer::edmHepMCProductToken_
private

Definition at line 104 of file PrimaryVertexAnalyzer.h.

Referenced by analyze().

edm::EDGetTokenT< edm::SimTrackContainer > PrimaryVertexAnalyzer::edmSimTrackContainerToken_
private

Definition at line 103 of file PrimaryVertexAnalyzer.h.

Referenced by analyze().

edm::EDGetTokenT< edm::SimVertexContainer > PrimaryVertexAnalyzer::edmSimVertexContainerToken_
private

Definition at line 102 of file PrimaryVertexAnalyzer.h.

Referenced by analyze().

std::map<std::string, TH1*> PrimaryVertexAnalyzer::h
private

Definition at line 109 of file PrimaryVertexAnalyzer.h.

Referenced by analyze(), beginJob(), endJob(), and getSimPVs().

std::map<std::string, TDirectory*> PrimaryVertexAnalyzer::hdir
private

Definition at line 110 of file PrimaryVertexAnalyzer.h.

Referenced by analyze(), beginJob(), and endJob().

std::string PrimaryVertexAnalyzer::outputFile_
private

Definition at line 106 of file PrimaryVertexAnalyzer.h.

Referenced by PrimaryVertexAnalyzer().

edm::ESHandle< ParticleDataTable > PrimaryVertexAnalyzer::pdt
private

Definition at line 100 of file PrimaryVertexAnalyzer.h.

Referenced by analyze(), isCharged(), and isResonance().

edm::EDGetTokenT< reco::TrackCollection > PrimaryVertexAnalyzer::recoTrackCollectionToken_
private

Definition at line 101 of file PrimaryVertexAnalyzer.h.

Referenced by analyze().

std::vector< edm::EDGetTokenT< reco::VertexCollection > > PrimaryVertexAnalyzer::recoVertexCollectionTokens_
private

Definition at line 105 of file PrimaryVertexAnalyzer.h.

Referenced by analyze(), and PrimaryVertexAnalyzer().

TFile* PrimaryVertexAnalyzer::rootFile_
private
double PrimaryVertexAnalyzer::simUnit_
private
std::vector<std::string> PrimaryVertexAnalyzer::suffixSample_
private

Definition at line 107 of file PrimaryVertexAnalyzer.h.

Referenced by analyze(), beginJob(), and PrimaryVertexAnalyzer().

bool PrimaryVertexAnalyzer::verbose_
private

Definition at line 98 of file PrimaryVertexAnalyzer.h.

Referenced by analyze(), getSimPVs(), and isResonance().