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

Inheritance diagram for PrimaryVertexAnalyzer:

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< simPrimaryVertex >	getSimPVs (const edm::Handle< edm::HepMCProduct > &evtMC, const std::string &suffix="")

std::vector< simPrimaryVertex >	getSimPVs (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< ParticleDataTable >	pdt

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_.

   : verbose_( iConfig.getUntrackedParameter<bool>( "verbose", false ) )
   , simUnit_( 1.0 ) // starting with CMSSW_1_2_x ??
   , recoTrackCollectionToken_( consumes< reco::TrackCollection >( edm::InputTag( iConfig.getUntrackedParameter<std::string>( "recoTrackProducer" ) ) ) )
   , edmSimVertexContainerToken_( consumes< edm::SimVertexContainer>( iConfig.getParameter<edm::InputTag>( "simG4" ) ) )
   , edmSimTrackContainerToken_( consumes< edm::SimTrackContainer >( iConfig.getParameter<edm::InputTag>( "simG4" ) ) )
   , edmHepMCProductToken_( consumes< edm::HepMCProduct >( edm::InputTag( std::string( "generator" ) ) ) )
 {
    //now do what ever initialization is needed
   // open output file to store histograms
   outputFile_ = iConfig.getUntrackedParameter<std::string>("outputFile");
   TString tversion(edm::getReleaseVersion());
   tversion = tversion.Remove(0,1);
   tversion = tversion.Remove(tversion.Length()-1,tversion.Length());
   outputFile_ = std::string(tversion)+"_"+outputFile_;
 
   std::vector<std::string> vtxSample = iConfig.getUntrackedParameter<std::vector< std::string > >("vtxSample");
   recoVertexCollectionTokens_.reserve( vtxSample.size() );
   suffixSample_.reserve( vtxSample.size() );
   auto vtxSampleBegin = vtxSample.begin();
   auto vtxSampleEnd = vtxSample.end();
   for( auto isample = vtxSampleBegin; isample != vtxSampleEnd; ++isample ) {
     recoVertexCollectionTokens_.push_back( consumes< reco::VertexCollection >( edm::InputTag( *isample ) ) );
     if ( *isample == "offlinePrimaryVertices" ) suffixSample_.push_back("AVF");
     if ( *isample == "offlinePrimaryVerticesWithBS" ) suffixSample_.push_back("wBS");
   }
   if ( suffixSample_.size() == 0 ) throw cms::Exception("NoVertexSamples") << " no known vertex samples given";
 
   rootFile_ = new TFile(outputFile_.c_str(),"RECREATE");
   if ( (edm::getReleaseVersion()).find("CMSSW_1_1_",0)!=std::string::npos){
     simUnit_=0.1;  // for use in  CMSSW_1_1_1 tutorial
   }
 }

PrimaryVertexAnalyzer::~PrimaryVertexAnalyzer ( )

Definition at line 80 of file PrimaryVertexAnalyzer.cc.

References rootFile_.

 {
     // do anything here that needs to be done at desctruction time
    // (e.g. close files, deallocate resources etc.)
    delete rootFile_;
 }

Member Function Documentation

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

virtual

Implements edm::EDAnalyzer.

Definition at line 402 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.

 {
   
   edm::Handle<reco::TrackCollection> recTrks;
   iEvent.getByToken( recoTrackCollectionToken_, recTrks );
 
   edm::Handle<edm::SimVertexContainer> simVtxs;
   iEvent.getByToken( edmSimVertexContainerToken_, simVtxs );
   
   edm::Handle<edm::SimTrackContainer> simTrks;
   iEvent.getByToken( edmSimTrackContainerToken_, simTrks );
 
   for (int ivtxSample=0; ivtxSample!= (int)recoVertexCollectionTokens_.size(); ++ivtxSample) {
     std::string isuffix = suffixSample_[ivtxSample];
 
     edm::Handle<reco::VertexCollection> recVtxs;
     iEvent.getByToken( recoVertexCollectionTokens_[ ivtxSample ], recVtxs );
 
     try{
       iSetup.getData(pdt);
     }catch(const edm::Exception&){
       std::cout << "Some problem occurred with the particle data table. This may not work !." <<std::endl;
     }
 
     bool MC=false;
     edm::Handle<edm::HepMCProduct> evtMC;
     iEvent.getByToken( edmHepMCProductToken_, evtMC );
     if (!evtMC.isValid()) {
       MC=false;
       if(verbose_){
     std::cout << "no HepMCProduct found"<< std::endl;
       }
     } else {
       if(verbose_){
     std::cout << "generator HepMCProduct found"<< std::endl;
       }
       MC=true;
     }
     
     
     if(verbose_){
       //evtMC->GetEvent()->print();
       printRecVtxs(recVtxs);
       //printSimVtxs(simVtxs);
       //printSimTrks(simTrks);
     }
     
     if(MC){
       // make a list of primary vertices:
       std::vector<simPrimaryVertex> simpv;
       //simpv=getSimPVs(evtMC, simVtxs, simTrks);
       simpv=getSimPVs(evtMC,isuffix);
       
       
       // vertex matching and efficiency bookkeeping
       h["nsimvtx"+isuffix]->Fill(simpv.size());
       int nsimtrk=0;
       for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
       vsim!=simpv.end(); vsim++){
     
     hdir[isuffix]->cd();
     
     h["nbsimtksinvtx"+isuffix]->Fill(vsim->nGenTrk);
     h["xsim"+isuffix]->Fill(vsim->x*simUnit_);
     h["ysim"+isuffix]->Fill(vsim->y*simUnit_);
     h["zsim"+isuffix]->Fill(vsim->z*simUnit_);
     nsimtrk+=vsim->nGenTrk;
     // look for a matching reconstructed vertex
     vsim->recVtx=NULL;
     for(reco::VertexCollection::const_iterator vrec=recVtxs->begin(); 
         vrec!=recVtxs->end(); ++vrec){
       if ( matchVertex(*vsim,*vrec) ){
         // if the matching critera are fulfilled, accept the rec-vertex that is closest in z
         if(    ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z)>fabs(vrec->z()-vsim->z)))
            || (!vsim->recVtx) )
           {
         vsim->recVtx=&(*vrec);
           }
       }
     }
     h["nsimtrk"+isuffix]->Fill(float(nsimtrk));
     
     
     // histogram properties of matched vertices
     if (vsim->recVtx){
       
       if(verbose_){std::cout <<"primary matched " << vsim->x << " " << vsim->y << " " << vsim->z << std:: endl;}
       
       h["resx"+isuffix]->Fill( vsim->recVtx->x()-vsim->x*simUnit_ );
       h["resy"+isuffix]->Fill( vsim->recVtx->y()-vsim->y*simUnit_ );
       h["resz"+isuffix]->Fill( vsim->recVtx->z()-vsim->z*simUnit_ );
       h["pullx"+isuffix]->Fill( (vsim->recVtx->x()-vsim->x*simUnit_)/vsim->recVtx->xError() );
       h["pully"+isuffix]->Fill( (vsim->recVtx->y()-vsim->y*simUnit_)/vsim->recVtx->yError() );
       h["pullz"+isuffix]->Fill( (vsim->recVtx->z()-vsim->z*simUnit_)/vsim->recVtx->zError() );
       h["eff"+isuffix]->Fill( 1.);
       if(simpv.size()==1){
         if (vsim->recVtx==&(*recVtxs->begin())){
           h["efftag"+isuffix]->Fill( 1.); 
         }else{
           h["efftag"+isuffix]->Fill( 0.); 
         }
       }
       h["effvseta"+isuffix]->Fill(vsim->ptot.pseudoRapidity(),1.);
       h["effvsptsq"+isuffix]->Fill(vsim->ptsq,1.);
       h["effvsntrk"+isuffix]->Fill(vsim->nGenTrk,1.);
       h["effvsnrectrk"+isuffix]->Fill(recTrks->size(),1.);
       h["effvsz"+isuffix]->Fill(vsim->z*simUnit_,1.);
       
     }else{  // no rec vertex found for this simvertex
       
       if(verbose_){std::cout <<"primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << std::endl;}
       
       h["eff"+isuffix]->Fill( 0.);
       if(simpv.size()==1){ h["efftag"+isuffix]->Fill( 0.); }
       h["effvseta"+isuffix]->Fill(vsim->ptot.pseudoRapidity(),0.);
       h["effvsptsq"+isuffix]->Fill(vsim->ptsq,0.);
       h["effvsntrk"+isuffix]->Fill(float(vsim->nGenTrk),0.);
       h["effvsnrectrk"+isuffix]->Fill(recTrks->size(),0.);
       h["effvsz"+isuffix]->Fill(vsim->z*simUnit_,0.);
     } // no recvertex for this simvertex
       }//found MC event
       // end of sim/rec matching 
       
       
       // purity of event vertex tags
       if (recVtxs->size()>0 && simpv.size()>0){
     Double_t dz=(*recVtxs->begin()).z() - (*simpv.begin()).z*simUnit_;
     h["zdistancetag"+isuffix]->Fill(dz);
     if( fabs(dz)<0.0500){
       h["puritytag"+isuffix]->Fill(1.);
     }else{
       // bad tag: the true primary was more than 500 um away from the tagged primary
       h["puritytag"+isuffix]->Fill(0.);
     }
       }
       
       
     }// MC event
     
     // test track links, use reconstructed vertices
     
     h["nrecvtx"+isuffix]->Fill(recVtxs->size());
     h["nrectrk"+isuffix]->Fill(recTrks->size());
     h["nbvtx"+isuffix]->Fill(recVtxs->size()*1.);
     if((recVtxs->size()==0)||recVtxs->begin()->isFake()) {h["nrectrk0vtx"+isuffix]->Fill(recTrks->size());}
     
     for(reco::VertexCollection::const_iterator v=recVtxs->begin(); 
     v!=recVtxs->end(); ++v){
 
       if(v->isFake()) continue;
       
       try {
     for(reco::Vertex::trackRef_iterator t = v->tracks_begin(); 
         t!=v->tracks_end(); t++) {
       // illegal charge
       if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
         h["tklinks"+isuffix]->Fill(0.);
       }
       else {
         h["tklinks"+isuffix]->Fill(1.);
       }
     }
       }
       catch (...) {
     // exception thrown when trying to use linked track
     h["tklinks"+isuffix]->Fill(0.);
       }
       
       h["nbtksinvtx"+isuffix]->Fill(v->tracksSize());
       h["vtxchi2"+isuffix]->Fill(v->chi2());
       h["vtxndf"+isuffix]->Fill(v->ndof());
       h["vtxprob"+isuffix]->Fill(ChiSquaredProbability(v->chi2() ,v->ndof()));
       h["xrec"+isuffix]->Fill(v->position().x());
       h["yrec"+isuffix]->Fill(v->position().y());
       h["zrec"+isuffix]->Fill(v->position().z());
       if (v==recVtxs->begin()){
     h["xrectag"+isuffix]->Fill(v->position().x());
     h["yrectag"+isuffix]->Fill(v->position().y());
     h["zrectag"+isuffix]->Fill(v->position().z());
       }
       
       bool problem = false;
       h["nans"+isuffix]->Fill(1.,isnan(v->position().x())*1.);
       h["nans"+isuffix]->Fill(2.,isnan(v->position().y())*1.);
       h["nans"+isuffix]->Fill(3.,isnan(v->position().z())*1.);
       
       int index = 3;
       for (int i = 0; i != 3; i++) {
     for (int j = i; j != 3; j++) {
       index++;
       h["nans"+isuffix]->Fill(index*1., isnan(v->covariance(i, j))*1.);
       if (isnan(v->covariance(i, j))) problem = true;
       // in addition, diagonal element must be positive
       if (j == i && v->covariance(i, j) < 0) {
         h["nans"+isuffix]->Fill(index*1., 1.);
         problem = true;
       }
     }
       }
       
       if (problem) {
     // analyze track parameter covariance definiteness
     double data[25];
     try {
       int itk = 0;
       for(reco::Vertex::trackRef_iterator t = v->tracks_begin(); 
           t!=v->tracks_end(); t++) {
         std::cout << "Track " << itk++ << std::endl;
         int i2 = 0;
         for (int i = 0; i != 5; i++) {
           for (int j = 0; j != 5; j++) {
         data[i2] = (**t).covariance(i, j);
         std::cout << std:: scientific << data[i2] << " ";
         i2++;
           }
           std::cout << std::endl;
         }
         gsl_matrix_view m 
           = gsl_matrix_view_array (data, 5, 5);
         
         gsl_vector *eval = gsl_vector_alloc (5);
         gsl_matrix *evec = gsl_matrix_alloc (5, 5);
         
         gsl_eigen_symmv_workspace * w = 
           gsl_eigen_symmv_alloc (5);
       
         gsl_eigen_symmv (&m.matrix, eval, evec, w);
         
         gsl_eigen_symmv_free (w);
         
         gsl_eigen_symmv_sort (eval, evec, 
                   GSL_EIGEN_SORT_ABS_ASC);
         
         // print sorted eigenvalues
         {
           int i;
           for (i = 0; i < 5; i++) {
         double eval_i 
           = gsl_vector_get (eval, i);
         //gsl_vector_view evec_i 
         //  = gsl_matrix_column (evec, i);
         
         printf ("eigenvalue = %g\n", eval_i);
         //        printf ("eigenvector = \n");
         //        gsl_vector_fprintf (stdout, 
         //                &evec_i.vector, "%g");
           }
         }
       }
     }
     catch (...) {
       // exception thrown when trying to use linked track
       break;
     }
       }// if (problem)
     }
   } // for vertex loop
 }

void PrimaryVertexAnalyzer::beginJob ( void )

virtual

Reimplemented from edm::EDAnalyzer.

Definition at line 91 of file PrimaryVertexAnalyzer.cc.

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

                                     {
   std::cout << " PrimaryVertexAnalyzer::beginJob  conversion from sim units to rec units is " << simUnit_ << std::endl;
 
 
   rootFile_->cd();
   // release validation histograms used in DoCompare.C
   for (std::vector<std::string>::iterator isuffix= suffixSample_.begin(); isuffix!=suffixSample_.end(); isuffix++) {
 
     hdir[*isuffix] = rootFile_->mkdir(TString(*isuffix));
     hdir[*isuffix]->cd();
     h["nbvtx"+ *isuffix]        = new TH1F(TString("nbvtx"+ *isuffix),"Reconstructed Vertices in Event",20,-0.5,19.5);
     h["nbtksinvtx"+ *isuffix]   = new TH1F(TString("nbtksinvtx"+ *isuffix),"Reconstructed Tracks in Vertex",200,-0.5,199.5); 
     h["resx"+ *isuffix]         = new TH1F(TString("resx"+ *isuffix),"Residual X",100,-0.04,0.04);
     h["resy"+ *isuffix]         = new TH1F(TString("resy"+ *isuffix),"Residual Y",100,-0.04,0.04);
     h["resz"+ *isuffix]         = new TH1F(TString("resz"+ *isuffix),"Residual Z",100,-0.1,0.1);
     h["pullx"+ *isuffix]        = new TH1F(TString("pullx"+ *isuffix),"Pull X",100,-25.,25.);
     h["pully"+ *isuffix]        = new TH1F(TString("pully"+ *isuffix),"Pull Y",100,-25.,25.);
     h["pullz"+ *isuffix]        = new TH1F(TString("pullz"+ *isuffix),"Pull Z",100,-25.,25.);
     h["vtxchi2"+ *isuffix]      = new TH1F(TString("vtxchi2"+ *isuffix),"#chi^{2}",100,0.,100.);
     h["vtxndf"+ *isuffix]       = new TH1F(TString("vtxndf"+ *isuffix),"ndof",100,0.,100.);
     h["tklinks"+ *isuffix]      = new TH1F(TString("tklinks"+ *isuffix),"Usable track links",2,-0.5,1.5);
     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);
     // more histograms
     h["vtxprob"+ *isuffix]      = new TH1F(TString("vtxprob"+ *isuffix),"#chi^{2} probability",100,0.,1.);
     h["eff"+ *isuffix]          = new TH1F(TString("eff"+ *isuffix),"efficiency",2, -0.5, 1.5);
     h["efftag"+ *isuffix]       = new TH1F(TString("efftag"+ *isuffix),"efficiency tagged vertex",2, -0.5, 1.5);
     h["effvseta"+ *isuffix]     = new TProfile(TString("effvseta"+ *isuffix),"efficiency vs eta",20, -2.5, 2.5, 0, 1.);
     h["effvsptsq"+ *isuffix]    = new TProfile(TString("effvsptsq"+ *isuffix),"efficiency vs ptsq",20, 0., 10000., 0, 1.);
     h["effvsntrk"+ *isuffix]    = new TProfile(TString("effvsntrk"+ *isuffix),"efficiency vs # simtracks",200, 0., 200., 0, 1.);
     h["effvsnrectrk"+ *isuffix] = new TProfile(TString("effvsnrectrk"+ *isuffix),"efficiency vs # rectracks",200, 0., 200., 0, 1.);
     h["effvsz"+ *isuffix]       = new TProfile(TString("effvsz"+ *isuffix),"efficiency vs z",40, -20., 20., 0, 1.);
     h["nbsimtksinvtx"+ *isuffix]= new TH1F(TString("nbsimtksinvtx"+ *isuffix),"simulated tracks in vertex",100,-0.5,99.5); 
     h["xrec"+ *isuffix]         = new TH1F(TString("xrec"+ *isuffix),"reconstructed x",100,-0.1,0.1);
     h["yrec"+ *isuffix]         = new TH1F(TString("yrec"+ *isuffix),"reconstructed y",100,-0.1,0.1);
     h["zrec"+ *isuffix]         = new TH1F(TString("zrec"+ *isuffix),"reconstructed z",100,-20.,20.);
     h["xsim"+ *isuffix]         = new TH1F(TString("xsim"+ *isuffix),"simulated x",100,-0.1,0.1);
     h["ysim"+ *isuffix]         = new TH1F(TString("ysim"+ *isuffix),"simulated y",100,-0.1,0.1);
     h["zsim"+ *isuffix]         = new TH1F(TString("zsim"+ *isuffix),"simulated z",100,-20.,20.);
     h["nrecvtx"+ *isuffix]      = new TH1F(TString("nrecvtx"+ *isuffix),"# of reconstructed vertices", 50, -0.5, 49.5);
     h["nsimvtx"+ *isuffix]      = new TH1F(TString("nsimvtx"+ *isuffix),"# of simulated vertices", 50, -0.5, 49.5);
     h["nrectrk"+ *isuffix]      = new TH1F(TString("nrectrk"+ *isuffix),"# of reconstructed tracks", 200, -0.5, 199.5);
     h["nsimtrk"+ *isuffix]      = new TH1F(TString("nsimtrk"+ *isuffix),"# of simulated tracks", 200, -0.5, 199.5);
     h["xrectag"+ *isuffix]      = new TH1F(TString("xrectag"+ *isuffix),"reconstructed x, signal vtx",100,-0.1,0.1);
     h["yrectag"+ *isuffix]      = new TH1F(TString("yrectag"+ *isuffix),"reconstructed y, signal vtx",100,-0.1,0.1);
     h["zrectag"+ *isuffix]      = new TH1F(TString("zrectag"+ *isuffix),"reconstructed z, signal vtx",100,-20.,20.);
     h["rapidity"+ *isuffix] = new TH1F(TString("rapidity"+ *isuffix),"rapidity ",100,-10., 10.);
     h["pt"+ *isuffix] = new TH1F(TString("pt"+ *isuffix),"pt ",100,0., 20.);
     h["nrectrk0vtx"+ *isuffix] = new TH1F(TString("nrectrk0vtx"+ *isuffix),"# rec tracks no vertex ",200,0., 200.);
     h["zdistancetag"+ *isuffix] = new TH1F(TString("zdistancetag"+ *isuffix),"z-distance between tagged and generated",100, -0.1, 0.1);
     h["puritytag"+ *isuffix]    = new TH1F(TString("puritytag"+ *isuffix),"purity of primary vertex tags",2, -0.5, 1.5);
   }
 }

void PrimaryVertexAnalyzer::endJob ( void )

virtual

Reimplemented from edm::EDAnalyzer.

Definition at line 145 of file PrimaryVertexAnalyzer.cc.

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

                                    {
   rootFile_->cd();
   // save all histograms created in beginJob()
   for (std::map<std::string, TDirectory*>::const_iterator idir=hdir.begin(); idir!=hdir.end(); idir++){
     idir->second->cd();
     for(std::map<std::string,TH1*>::const_iterator hist=h.begin(); hist!=h.end(); hist++){
       if (TString(hist->first).Contains(idir->first)) hist->second->Write();
     }
   }
 }

std::vector< PrimaryVertexAnalyzer::simPrimaryVertex > PrimaryVertexAnalyzer::getSimPVs	(	const edm::Handle< edm::HepMCProduct > &	evtMC,
		const std::string &	suffix = `""`
	)

private

Definition at line 237 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().

 {
   std::vector<PrimaryVertexAnalyzer::simPrimaryVertex> simpv;
   const HepMC::GenEvent* evt=evtMC->GetEvent();
   if (evt) {
     if(verbose_) std::cout << "process id " <<evt->signal_process_id()<<std::endl;
     if(verbose_) std::cout <<"signal process vertex "<< ( evt->signal_process_vertex() ?
                          evt->signal_process_vertex()->barcode() : 0 )   <<std::endl;
     if(verbose_) std::cout <<"number of vertices " << evt->vertices_size() << std::endl;
 
 
     int idx=0;
     for(HepMC::GenEvent::vertex_const_iterator vitr= evt->vertices_begin();
     vitr != evt->vertices_end(); ++vitr ) 
       { // loop for vertex ...
     HepMC::FourVector pos = (*vitr)->position();
     //HepLorentzVector pos = (*vitr)->position();
 
     bool hasMotherVertex=false;
     if(verbose_) std::cout << "mothers of vertex[" << ++idx << "]: " << std::endl;
     for ( HepMC::GenVertex::particle_iterator
           mother  = (*vitr)->particles_begin(HepMC::parents);
           mother != (*vitr)->particles_end(HepMC::parents);
               ++mother ) {
       HepMC::GenVertex * mv=(*mother)->production_vertex();
       if (mv) {
         hasMotherVertex=true;
         if(!verbose_) break; //if verbose_, print all particles of gen vertices
       }
       if(verbose_)std::cout << "\t";
       if(verbose_)(*mother)->print();
     }
 
     if(hasMotherVertex) {continue;}
 
     // could be a new vertex, check  all primaries found so far to avoid multiple entries
         const double mm=0.1;
     simPrimaryVertex sv(pos.x()*mm,pos.y()*mm,pos.z()*mm);
     simPrimaryVertex *vp=NULL;  // will become non-NULL if a vertex is found and then point to it
     for(std::vector<simPrimaryVertex>::iterator v0=simpv.begin();
         v0!=simpv.end(); v0++){
       if( (fabs(sv.x-v0->x)<1e-5) && (fabs(sv.y-v0->y)<1e-5) && (fabs(sv.z-v0->z)<1e-5)){
         vp=&(*v0);
         break;
       }
     }
 
     if(!vp){
       // this is a new vertex
       if(verbose_)std::cout << "this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << std::endl;
       simpv.push_back(sv);
       vp=&simpv.back();
     }else{
       if(verbose_)std::cout << "this is not new vertex " << sv.x << " " << sv.y << " " << sv.z << std::endl;
     }
     vp->genVertex.push_back((*vitr)->barcode());
     // collect final state descendants
     for ( HepMC::GenVertex::particle_iterator
           daughter  = (*vitr)->particles_begin(HepMC::descendants);
           daughter != (*vitr)->particles_end(HepMC::descendants);
               ++daughter ) {
       if (isFinalstateParticle(*daughter)){ 
         if ( find(vp->finalstateParticles.begin(), vp->finalstateParticles.end(),(*daughter)->barcode())
          == vp->finalstateParticles.end()){
           vp->finalstateParticles.push_back((*daughter)->barcode());
           HepMC::FourVector m=(*daughter)->momentum();
           // the next four lines used to be "vp->ptot+=m;" in the days of CLHEP::HepLorentzVector
           // but adding FourVectors seems not to be foreseen
           vp->ptot.setPx(vp->ptot.px()+m.px());
           vp->ptot.setPy(vp->ptot.py()+m.py());
           vp->ptot.setPz(vp->ptot.pz()+m.pz());
           vp->ptot.setE(vp->ptot.e()+m.e());
           vp->ptsq+=(m.perp())*(m.perp());
           if ( (m.perp()>0.8) && (fabs(m.pseudoRapidity())<2.5) && isCharged( *daughter ) ){
         vp->nGenTrk++;
           }
           
           h["rapidity"+suffix]->Fill(m.pseudoRapidity());
           h["pt"+suffix]->Fill(m.perp());
         }
       }
     }
       }
   }
   return simpv;
 }

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 328 of file PrimaryVertexAnalyzer.cc.

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

 {
    // simvertices don't have enough information to decide, 
    // genVertices don't have the simulated coordinates  ( with VtxSmeared they might)
    // go through simtracks to get the link between simulated and generated vertices
    std::vector<PrimaryVertexAnalyzer::simPrimaryVertex> simpv;
    int idx=0;
    for(edm::SimTrackContainer::const_iterator t=simTrks->begin();
        t!=simTrks->end(); ++t){
      if ( !(t->noVertex()) && !(t->type()==-99) ){
        double ptsq=0;
        bool primary=false;   // something coming directly from the primary vertex
        bool resonance=false; // resonance
        bool track=false;     // undecayed, charged particle
        HepMC::GenParticle* gp=evtMC->GetEvent()->barcode_to_particle( (*t).genpartIndex() );
        if (gp) {
      HepMC::GenVertex * gv=gp->production_vertex();
      if (gv) {
        for ( HepMC::GenVertex::particle_iterator 
                  daughter =  gv->particles_begin(HepMC::descendants);
          daughter != gv->particles_end(HepMC::descendants);
          ++daughter ) {
          if (isFinalstateParticle(*daughter)){
            ptsq+=(*daughter)->momentum().perp()*(*daughter)->momentum().perp();
          }
        }
        primary =  ( gv->position().t()==0);
        //resonance= ( gp->mother() && isResonance(gp->mother()));  // in CLHEP/HepMC days
        // no more mother pointer in the improved HepMC GenParticle
        resonance= ( isResonance(*(gp->production_vertex()->particles_in_const_begin())));
        if (gp->status()==1){
          //track=((pdt->particle(gp->pdg_id()))->charge() != 0);
          track=not isCharged(gp);
        }
      }
        }
 
        const HepMC::FourVector & v=(*simVtxs)[t->vertIndex()].position();
        //const HepLorentzVector & v=(*simVtxs)[t->vertIndex()].position();
        if(primary or resonance){
      {
        // check all primaries found so far to avoid multiple entries
        bool newVertex=true;
        for(std::vector<simPrimaryVertex>::iterator v0=simpv.begin();
            v0!=simpv.end(); v0++){
          if( (fabs(v0->x-v.x())<0.001) && (fabs(v0->y-v.y())<0.001) && (fabs(v0->z-v.z())<0.001) ){
            if (track) {
          v0->simTrackIndex.push_back(idx);
          if (ptsq>(*v0).ptsq){(*v0).ptsq=ptsq;}
            }
            newVertex=false;
          }
        }
        if(newVertex && !resonance){
          simPrimaryVertex anotherVertex(v.x(),v.y(),v.z());
          if (track) anotherVertex.simTrackIndex.push_back(idx);
          anotherVertex.ptsq=ptsq;
          simpv.push_back(anotherVertex);
        }
      }// 
        }
        
      }// simtrack has vertex and valid type
      idx++;
    }//simTrack loop
    return simpv;
 }

bool PrimaryVertexAnalyzer::isCharged ( const HepMC::GenParticle * p )

private

Definition at line 173 of file PrimaryVertexAnalyzer.cc.

References pdt.

Referenced by getSimPVs().

                                                                {
   const ParticleData * part = pdt->particle( p->pdg_id() );
   if (part){
     return part->charge()!=0;
   }else{
     // the new/improved particle table doesn't know anti-particles
     return  pdt->particle( -p->pdg_id() )!=0;
   }
 }

bool PrimaryVertexAnalyzer::isFinalstateParticle ( const HepMC::GenParticle * p )

private

Definition at line 169 of file PrimaryVertexAnalyzer.cc.

Referenced by getSimPVs().

                                                                           {
   return ( !p->end_vertex() && p->status()==1 );
 }

bool PrimaryVertexAnalyzer::isResonance ( const HepMC::GenParticle * p )

private

Definition at line 163 of file PrimaryVertexAnalyzer.cc.

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

Referenced by getSimPVs().

                                                                  {
   double ctau=(pdt->particle( abs(p->pdg_id()) ))->lifetime();
   if(verbose_) std::cout << "isResonance   " << p->pdg_id() << " " << ctau << std::endl;
   return  ctau >0 && ctau <1e-6;
 }

bool PrimaryVertexAnalyzer::matchVertex	(	const simPrimaryVertex &	vsim,
		const reco::Vertex &	vrec
	)

private

Definition at line 158 of file PrimaryVertexAnalyzer.cc.

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

Referenced by analyze().

                                                     {
   return (fabs(vsim.z*simUnit_-vrec.z())<0.0500); // =500um
 }

void PrimaryVertexAnalyzer::printRecVtxs ( const edm::Handle< reco::VertexCollection > & recVtxs )

private

Definition at line 183 of file PrimaryVertexAnalyzer.cc.

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

Referenced by analyze().

                                                                                        {
     int ivtx=0;
     for(reco::VertexCollection::const_iterator v=recVtxs->begin(); 
     v!=recVtxs->end(); ++v){
       std::cout << "Recvtx "<< std::setw(3) << std::setfill(' ')<<ivtx++
         << "#trk " << std::setw(3) << v->tracksSize() 
         << " chi2 " << std::setw(4) << v->chi2() 
         << " ndof " << std::setw(3) << v->ndof() << std::endl 
         << " x "  << std::setw(8) <<std::fixed << std::setprecision(4) << v->x() 
         << " dx " << std::setw(8) << v->xError()<< std::endl
         << " y "  << std::setw(8) << v->y() 
         << " dy " << std::setw(8) << v->yError()<< std::endl
         << " z "  << std::setw(8) << v->z() 
         << " dz " << std::setw(8) << v->zError()
         << std::endl;
     }
 }

void PrimaryVertexAnalyzer::printSimTrks ( const edm::Handle< edm::SimTrackContainer > & simVtrks )

private

Definition at line 218 of file PrimaryVertexAnalyzer.cc.

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

                                                                                        {
   std::cout <<  " simTrks   type, (momentum), vertIndex, genpartIndex"  << std::endl;
   int i=1;
   for(edm::SimTrackContainer::const_iterator t=simTrks->begin();
       t!=simTrks->end(); ++t){
     //HepMC::GenParticle* gp=evtMC->GetEvent()->particle( (*t).genpartIndex() );
     std::cout << i++ << ")" 
           << (*t)
           << " index="
           << (*t).genpartIndex();
     //if (gp) {
     //  HepMC::GenVertex *gv=gp->production_vertex();
     //  std::cout  <<  " genvertex =" << (*gv);
     //}
     std::cout << std::endl;
   }
 }

void PrimaryVertexAnalyzer::printSimVtxs ( const edm::Handle< edm::SimVertexContainer > & simVtxs )

private

Definition at line 202 of file PrimaryVertexAnalyzer.cc.

References gather_cfg::cout, i, and simUnit_.

                                                                                         {
     int i=0;
     for(edm::SimVertexContainer::const_iterator vsim=simVtxs->begin();
     vsim!=simVtxs->end(); ++vsim){
       std::cout << i++ << ")" << std::scientific
                 << " evtid=" << vsim->eventId().event() 
         << " sim x=" << vsim->position().x()*simUnit_
         << " sim y=" << vsim->position().y()*simUnit_
         << " sim z=" << vsim->position().z()*simUnit_
         << " sim t=" << vsim->position().t()
         << " parent=" << vsim->parentIndex() 
         << std::endl;
     }
 }