|
|
|
#include <Validation/RecoVertex/src/PrimaryVertexAnalyzer.cc>
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 () | |
Private Member Functions | |
| std::vector< simPrimaryVertex > | getSimPVs (const edm::Handle< edm::HepMCProduct > evtMC, 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 | |
| std::map< std::string, TH1 * > | h |
| std::map< std::string, TDirectory * > | hdir |
| std::string | outputFile_ |
| edm::ESHandle< ParticleDataTable > | pdt |
| std::string | recoTrackProducer_ |
| TFile * | rootFile_ |
| edm::InputTag | simG4_ |
| double | simUnit_ |
| std::vector< std::string > | suffixSample_ |
| bool | verbose_ |
| std::vector< std::string > | vtxSample_ |
Description: simple primary vertex analyzer
Implementation: <Notes on="" implementation>="">
Definition at line 61 of file PrimaryVertexAnalyzer.h.
| PrimaryVertexAnalyzer::PrimaryVertexAnalyzer | ( | const edm::ParameterSet & | iConfig | ) | [explicit] |
Definition at line 49 of file PrimaryVertexAnalyzer.cc.
References spr::find(), edm::ParameterSet::getParameter(), edm::getReleaseVersion(), and edm::ParameterSet::getUntrackedParameter().
{
//now do what ever initialization is needed
simG4_=iConfig.getParameter<edm::InputTag>( "simG4" );
recoTrackProducer_= iConfig.getUntrackedParameter<std::string>("recoTrackProducer");
// 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_;
vtxSample_ = iConfig.getUntrackedParameter<std::vector< std::string > >("vtxSample");
for(std::vector<std::string>::iterator isample = vtxSample_.begin(); isample!=vtxSample_.end(); ++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");
verbose_= iConfig.getUntrackedParameter<bool>("verbose", false);
simUnit_= 1.0; // starting with CMSSW_1_2_x ??
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 77 of file PrimaryVertexAnalyzer.cc.
{
// do anything here that needs to be done at desctruction time
// (e.g. close files, deallocate resources etc.)
delete rootFile_;
/*for(std::map<std::string,TH1*>::const_iterator hist=h.begin(); hist!=h.end(); hist++){
delete hist->second;
}*/
}
| void PrimaryVertexAnalyzer::analyze | ( | const edm::Event & | iEvent, |
| const edm::EventSetup & | iSetup | ||
| ) | [virtual] |
Implements edm::EDAnalyzer.
Definition at line 406 of file PrimaryVertexAnalyzer.cc.
References ChiSquaredProbability(), gather_cfg::cout, runTheMatrix::data, edm::Event::getByLabel(), edm::EventSetup::getData(), h, i, getHLTprescales::index, edm::isnan(), edm::HandleBase::isValid(), j, m, NULL, matplotRender::t, v, and z.
{
Handle<reco::TrackCollection> recTrks;
iEvent.getByLabel(recoTrackProducer_, recTrks);
Handle<SimVertexContainer> simVtxs;
iEvent.getByLabel( simG4_, simVtxs);
Handle<SimTrackContainer> simTrks;
iEvent.getByLabel( simG4_, simTrks);
for (int ivtxSample=0; ivtxSample!= (int)vtxSample_.size(); ++ivtxSample) {
std::string isuffix = suffixSample_[ivtxSample];
Handle<reco::VertexCollection> recVtxs;
iEvent.getByLabel(vtxSample_[ivtxSample], recVtxs);
try{
iSetup.getData(pdt);
}catch(const Exception&){
std::cout << "Some problem occurred with the particle data table. This may not work !." <<std::endl;
}
bool MC=false;
Handle<HepMCProduct> evtMC;
iEvent.getByLabel("generator",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 93 of file PrimaryVertexAnalyzer.cc.
References gather_cfg::cout, and h.
{
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 148 of file PrimaryVertexAnalyzer.cc.
References h, and estimatePileup::hist.
{
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, |
| std::string | suffix = "" |
||
| ) | [private] |
Definition at line 240 of file PrimaryVertexAnalyzer.cc.
References gather_cfg::cout, spr::find(), h, m, NULL, parents, pos, and createPayload::suffix.
{
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 330 of file PrimaryVertexAnalyzer.cc.
References configurableAnalysis::GenParticle, resonance, PrimaryVertexAnalyzer::simPrimaryVertex::simTrackIndex, matplotRender::t, and 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(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 176 of file PrimaryVertexAnalyzer.cc.
| bool PrimaryVertexAnalyzer::isFinalstateParticle | ( | const HepMC::GenParticle * | p | ) | [private] |
Definition at line 172 of file PrimaryVertexAnalyzer.cc.
| bool PrimaryVertexAnalyzer::isResonance | ( | const HepMC::GenParticle * | p | ) | [private] |
Definition at line 166 of file PrimaryVertexAnalyzer.cc.
References abs, and gather_cfg::cout.
| bool PrimaryVertexAnalyzer::matchVertex | ( | const simPrimaryVertex & | vsim, |
| const reco::Vertex & | vrec | ||
| ) | [private] |
Definition at line 161 of file PrimaryVertexAnalyzer.cc.
References PrimaryVertexAnalyzer::simPrimaryVertex::z, and reco::Vertex::z().
| void PrimaryVertexAnalyzer::printRecVtxs | ( | const edm::Handle< reco::VertexCollection > | recVtxs | ) | [private] |
Definition at line 186 of file PrimaryVertexAnalyzer.cc.
References gather_cfg::cout, and v.
{
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 221 of file PrimaryVertexAnalyzer.cc.
References gather_cfg::cout, i, and matplotRender::t.
{
std::cout << " simTrks type, (momentum), vertIndex, genpartIndex" << std::endl;
int i=1;
for(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 205 of file PrimaryVertexAnalyzer.cc.
References gather_cfg::cout, and i.
{
int i=0;
for(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;
}
}
std::map<std::string, TH1*> PrimaryVertexAnalyzer::h [private] |
Definition at line 116 of file PrimaryVertexAnalyzer.h.
std::map<std::string, TDirectory*> PrimaryVertexAnalyzer::hdir [private] |
Definition at line 117 of file PrimaryVertexAnalyzer.h.
std::string PrimaryVertexAnalyzer::outputFile_ [private] |
Definition at line 106 of file PrimaryVertexAnalyzer.h.
edm::ESHandle< ParticleDataTable > PrimaryVertexAnalyzer::pdt [private] |
Definition at line 113 of file PrimaryVertexAnalyzer.h.
std::string PrimaryVertexAnalyzer::recoTrackProducer_ [private] |
Definition at line 105 of file PrimaryVertexAnalyzer.h.
TFile* PrimaryVertexAnalyzer::rootFile_ [private] |
Definition at line 109 of file PrimaryVertexAnalyzer.h.
edm::InputTag PrimaryVertexAnalyzer::simG4_ [private] |
Definition at line 111 of file PrimaryVertexAnalyzer.h.
double PrimaryVertexAnalyzer::simUnit_ [private] |
Definition at line 112 of file PrimaryVertexAnalyzer.h.
std::vector<std::string> PrimaryVertexAnalyzer::suffixSample_ [private] |
Definition at line 108 of file PrimaryVertexAnalyzer.h.
bool PrimaryVertexAnalyzer::verbose_ [private] |
Definition at line 110 of file PrimaryVertexAnalyzer.h.
std::vector<std::string> PrimaryVertexAnalyzer::vtxSample_ [private] |
Definition at line 107 of file PrimaryVertexAnalyzer.h.
1.7.3