|
|
|
#include <Validation/RecoVertex/src/PrimaryVertexAnalyzer4PU.cc>
Classes | |
| class | SimEvent |
| struct | SimPart |
| class | simPrimaryVertex |
Public Member Functions | |
| virtual void | analyze (const edm::Event &, const edm::EventSetup &) |
| virtual void | beginJob () |
| virtual void | endJob () |
| PrimaryVertexAnalyzer4PU (const edm::ParameterSet &) | |
| ~PrimaryVertexAnalyzer4PU () | |
Private Types | |
| typedef math::XYZTLorentzVector | LorentzVector |
| typedef reco::TrackBase::ParameterVector | ParameterVector |
Private Member Functions | |
| void | add (std::map< std::string, TH1 * > &h, TH1 *hist) |
| void | analyzeVertexCollection (std::map< std::string, TH1 * > &h, const edm::Handle< reco::VertexCollection > recVtxs, const edm::Handle< reco::TrackCollection > recTrks, std::vector< simPrimaryVertex > &simpv, const std::string message="") |
| void | analyzeVertexCollectionTP (std::map< std::string, TH1 * > &h, const edm::Handle< reco::VertexCollection > recVtxs, const edm::Handle< reco::TrackCollection > recTrks, std::vector< SimEvent > &simEvt, const std::string message="") |
| std::map< std::string, TH1 * > | bookVertexHistograms () |
| void | Cumulate (TH1 *h) |
| void | dumpHitInfo (const reco::Track &t) |
| void | Fill (std::map< std::string, TH1 * > &h, std::string s, bool yesno, bool signal) |
| void | Fill (std::map< std::string, TH1 * > &h, std::string s, double x, bool signal) |
| void | Fill (std::map< std::string, TH1 * > &h, std::string s, double x, double y) |
| void | Fill (std::map< std::string, TH1 * > &h, std::string s, double x) |
| void | fillTrackHistos (std::map< std::string, TH1 * > &h, const std::string &ttype, const reco::Track &t, const reco::Vertex *v=NULL) |
| std::vector < PrimaryVertexAnalyzer4PU::SimEvent > | getSimEvents (edm::Handle< TrackingParticleCollection >, edm::Handle< TrackingVertexCollection >, edm::Handle< edm::View< reco::Track > >) |
| std::vector< simPrimaryVertex > | getSimPVs (const edm::Handle< edm::HepMCProduct > evtMC) |
| std::vector< simPrimaryVertex > | getSimPVs (const edm::Handle< edm::HepMCProduct > evt, const edm::Handle< edm::SimVertexContainer > simVtxs, const edm::Handle< edm::SimTrackContainer > simTrks) |
| std::vector < PrimaryVertexAnalyzer4PU::simPrimaryVertex > | getSimPVs (const edm::Handle< TrackingVertexCollection >) |
| std::vector< SimPart > | getSimTrkParameters (edm::Handle< edm::SimTrackContainer > &simTrks, edm::Handle< edm::SimVertexContainer > &simVtcs, double simUnit=1.0) |
| void | getTc (const std::vector< reco::TransientTrack > &, double &, double &, double &, double &, double &) |
| std::vector< edm::RefToBase < reco::Track > > | getTruthMatchedVertexTracks (const reco::Vertex &) |
| void | history (const edm::Handle< edm::View< reco::Track > > &tracks, const size_t trackindex=10000) |
| bool | isCharged (const HepMC::GenParticle *p) |
| bool | isFinalstateParticle (const HepMC::GenParticle *p) |
| bool | isResonance (const HepMC::GenParticle *p) |
| void | matchRecTracksToVertex (simPrimaryVertex &pv, const std::vector< SimPart > &tsim, const edm::Handle< reco::TrackCollection > &recTrks) |
| bool | matchVertex (const simPrimaryVertex &vsim, const reco::Vertex &vrec) |
| std::string | particleString (int) const |
| void | printEventSummary (std::map< std::string, TH1 * > &h, const edm::Handle< reco::VertexCollection > recVtxs, const edm::Handle< reco::TrackCollection > recTrks, std::vector< SimEvent > &simEvt, const std::string message) |
| void | printPVTrks (const edm::Handle< reco::TrackCollection > &recTrks, const edm::Handle< reco::VertexCollection > &recVtxs, std::vector< SimPart > &tsim, std::vector< SimEvent > &simEvt, const bool selectedOnly=true) |
| void | printRecTrks (const edm::Handle< reco::TrackCollection > &recTrks) |
| void | printRecVtxs (const edm::Handle< reco::VertexCollection > recVtxs, std::string title="Reconstructed Vertices") |
| void | printSimTrks (const edm::Handle< edm::SimTrackContainer > simVtrks) |
| void | printSimVtxs (const edm::Handle< edm::SimVertexContainer > simVtxs) |
| int * | supf (std::vector< SimPart > &simtrks, const reco::TrackCollection &trks) |
| bool | truthMatchedTrack (edm::RefToBase< reco::Track >, TrackingParticleRef &) |
| std::string | vertexString (HepMC::GenVertex::particles_in_const_iterator, HepMC::GenVertex::particles_in_const_iterator, HepMC::GenVertex::particles_out_const_iterator, HepMC::GenVertex::particles_out_const_iterator) const |
| std::string | vertexString (TrackingParticleRefVector, TrackingParticleRefVector) const |
Static Private Member Functions | |
| static bool | match (const ParameterVector &a, const ParameterVector &b) |
Private Attributes | |
| TrackAssociatorBase * | associatorByHits_ |
| int | bunchCrossing_ |
| bool | DEBUG_ |
| bool | doMatching_ |
| int | dumpcounter_ |
| bool | dumpPUcandidates_ |
| bool | dumpThisEvent_ |
| int | event_ |
| int | eventcounter_ |
| double | fBfield_ |
| std::map< std::string, TH1 * > | hBS |
| std::map< std::string, TH1 * > | hDA |
| std::map< std::string, TH1 * > | hMVF |
| std::map< std::string, TH1 * > | hnoBS |
| std::map< std::string, TH1 * > | hPIX |
| std::map< std::string, TH1 * > | hsimPV |
| int | luminosityBlock_ |
| math::XYZPoint | myBeamSpot |
| int | ndump_ |
| int | orbitNumber_ |
| std::string | outputFile_ |
| edm::ESHandle< ParticleDataTable > | pdt_ |
| bool | printXBS_ |
| reco::RecoToSimCollection | r2s_ |
| edm::Handle< reco::BeamSpot > | recoBeamSpotHandle_ |
| std::string | recoTrackProducer_ |
| TFile * | rootFile_ |
| int | run_ |
| edm::InputTag | simG4_ |
| double | simUnit_ |
| edm::ESHandle < TransientTrackBuilder > | theB_ |
| TrackFilterForPVFinding | theTrackFilter |
| bool | verbose_ |
| reco::BeamSpot | vertexBeamSpot_ |
| std::vector< std::string > | vtxSample_ |
| double | wxy2_ |
| std::map< double, TrackingParticleRef > | z2tp_ |
| double | zmatch_ |
Description: primary vertex analyzer for events with pile-up
Implementation: <Notes on="" implementation>="">
Definition at line 75 of file PrimaryVertexAnalyzer4PU.h.
typedef math::XYZTLorentzVector PrimaryVertexAnalyzer4PU::LorentzVector [private] |
Definition at line 77 of file PrimaryVertexAnalyzer4PU.h.
typedef reco::TrackBase::ParameterVector PrimaryVertexAnalyzer4PU::ParameterVector [private] |
Definition at line 79 of file PrimaryVertexAnalyzer4PU.h.
| PrimaryVertexAnalyzer4PU::PrimaryVertexAnalyzer4PU | ( | const edm::ParameterSet & | iConfig | ) | [explicit] |
Definition at line 74 of file PrimaryVertexAnalyzer4PU.cc.
References gather_cfg::cout, DEBUG_, doMatching_, dumpcounter_, dumpPUcandidates_, eventcounter_, spr::find(), edm::ParameterSet::getParameter(), edm::getReleaseVersion(), edm::ParameterSet::getUntrackedParameter(), ndump_, outputFile_, printXBS_, recoTrackProducer_, rootFile_, simG4_, simUnit_, verbose_, and zmatch_.
:theTrackFilter(iConfig.getParameter<edm::ParameterSet>("TkFilterParameters")) { //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"); rootFile_ = TFile::Open(outputFile_.c_str(),"RECREATE"); verbose_= iConfig.getUntrackedParameter<bool>("verbose", false); doMatching_= iConfig.getUntrackedParameter<bool>("matching", false); printXBS_= iConfig.getUntrackedParameter<bool>("XBS", 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 } dumpPUcandidates_=iConfig.getUntrackedParameter<bool>("dumpPUcandidates", false); zmatch_=iConfig.getUntrackedParameter<double>("zmatch", 0.0500); cout << "PrimaryVertexAnalyzer4PU: zmatch=" << zmatch_ << endl; eventcounter_=0; dumpcounter_=0; ndump_=10; DEBUG_=false; //DEBUG_=true; }
| PrimaryVertexAnalyzer4PU::~PrimaryVertexAnalyzer4PU | ( | ) |
Definition at line 105 of file PrimaryVertexAnalyzer4PU.cc.
References rootFile_.
{
// do anything here that needs to be done at desctruction time
// (e.g. close files, deallocate resources etc.)
delete rootFile_;
}
| void PrimaryVertexAnalyzer4PU::add | ( | std::map< std::string, TH1 * > & | h, |
| TH1 * | hist | ||
| ) | [inline, private] |
Definition at line 184 of file PrimaryVertexAnalyzer4PU.h.
References estimatePileup::hist.
Referenced by beginJob(), and bookVertexHistograms().
| void PrimaryVertexAnalyzer4PU::analyze | ( | const edm::Event & | iEvent, |
| const edm::EventSetup & | iSetup | ||
| ) | [virtual] |
Implements edm::EDAnalyzer.
Definition at line 1912 of file PrimaryVertexAnalyzer4PU.cc.
References analyzeVertexCollection(), analyzeVertexCollectionTP(), TrackAssociatorBase::associateRecoToSim(), associatorByHits_, align::BeamSpot, reco::BeamSpot::BeamWidthX(), reco::BeamSpot::BeamWidthY(), edm::EventBase::bunchCrossing(), bunchCrossing_, gather_cfg::cout, doMatching_, dumpcounter_, dumpThisEvent_, edm::EventID::event(), event_, eventcounter_, fBfield_, Fill(), edm::EventSetup::get(), edm::Event::getByLabel(), edm::Event::getByType(), edm::EventSetup::getData(), getSimEvents(), getSimPVs(), getSimTrkParameters(), hBS, hDA, reco::TrackBase::highPurity, hnoBS, hsimPV, edm::EventBase::id(), edm::HandleBase::isValid(), edm::EventBase::luminosityBlock(), luminosityBlock_, matchRecTracksToVertex(), myBeamSpot, ndump_, edm::EventBase::orbitNumber(), orbitNumber_, pdt_, funct::pow(), printPVTrks(), printRecVtxs(), printXBS_, edm::ESHandle< T >::product(), edm::Handle< T >::product(), r2s_, recoBeamSpotHandle_, recoTrackProducer_, edm::EventID::run(), run_, simG4_, simUnit_, lumiQTWidget::t, theB_, v, verbose_, vertexBeamSpot_, wxy2_, reco::BeamSpot::x0(), reco::BeamSpot::y0(), z, and reco::BeamSpot::z0().
{
std::vector<simPrimaryVertex> simpv; // a list of primary MC vertices
std::vector<SimPart> tsim;
std::string mcproduct="generator"; // starting with 3_1_0 pre something
eventcounter_++;
run_ = iEvent.id().run();
luminosityBlock_ = iEvent.luminosityBlock();
event_ = iEvent.id().event();
bunchCrossing_ = iEvent.bunchCrossing();
orbitNumber_ = iEvent.orbitNumber();
dumpThisEvent_ = false;
if(verbose_){
std::cout << endl
<< "PrimaryVertexAnalyzer4PU::analyze event counter=" << eventcounter_
<< " Run=" << run_ << " LumiBlock " << luminosityBlock_ << " event " << event_
<< " bx=" << bunchCrossing_ << " orbit=" << orbitNumber_
//<< " selected = " << good
<< std::endl;
}
try{
iSetup.getData(pdt_);
}catch(const Exception&){
std::cout << "Some problem occurred with the particle data table. This may not work !" <<std::endl;
}
Handle<reco::VertexCollection> recVtxs;
bool bnoBS=iEvent.getByLabel("offlinePrimaryVertices", recVtxs);
Handle<reco::VertexCollection> recVtxsBS;
bool bBS=iEvent.getByLabel("offlinePrimaryVerticesWithBS", recVtxsBS);
Handle<reco::VertexCollection> recVtxsDA;
bool bDA=iEvent.getByLabel("offlinePrimaryVerticesDA", recVtxsDA);
// Handle<reco::VertexCollection> recVtxsPIX;
// bool bPIX=iEvent.getByLabel("pixelVertices", recVtxsPIX);
// bPIX=false;
// Handle<reco::VertexCollection> recVtxsMVF;
// bool bMVF=iEvent.getByLabel("offlinePrimaryVerticesMVF", recVtxsMVF);
Handle<reco::TrackCollection> recTrks;
iEvent.getByLabel(recoTrackProducer_, recTrks);
int nhighpurity=0, ntot=0;
for(reco::TrackCollection::const_iterator t=recTrks->begin(); t!=recTrks->end(); ++t){
ntot++;
if(t->quality(reco::TrackBase::highPurity)) nhighpurity++;
}
if(ntot>10) hnoBS["highpurityTrackFraction"]->Fill(float(nhighpurity)/float(ntot));
if((recoTrackProducer_=="generalTracks")&&(nhighpurity<0.25*ntot)){
if(verbose_){ cout << "rejected, " << nhighpurity << " highpurity out of " << ntot << " total tracks "<< endl<< endl;}
return;
}
if(iEvent.getByType(recoBeamSpotHandle_)){
vertexBeamSpot_= *recoBeamSpotHandle_;
wxy2_=pow(vertexBeamSpot_.BeamWidthX(),2)+pow(vertexBeamSpot_.BeamWidthY(),2);
Fill(hsimPV, "xbeam",vertexBeamSpot_.x0()); Fill(hsimPV, "wxbeam",vertexBeamSpot_.BeamWidthX());
Fill(hsimPV, "ybeam",vertexBeamSpot_.y0()); Fill(hsimPV, "wybeam",vertexBeamSpot_.BeamWidthY());
Fill(hsimPV, "zbeam",vertexBeamSpot_.z0());// Fill("wzbeam",vertexBeamSpot_.BeamWidthZ());
}else{
cout << " beamspot not found, using dummy " << endl;
vertexBeamSpot_=reco::BeamSpot();// dummy
}
if(bnoBS){
if((recVtxs->begin()->isValid())&&(recVtxs->begin()->ndof()>1)&&(recVtxs->begin()->ndof()>(0.0*recVtxs->begin()->tracksSize()))){ // was 5 and 0.2
Fill(hnoBS,"xrecBeamvsdxXBS",recVtxs->begin()->xError(),recVtxs->begin()->x()-vertexBeamSpot_.x0());
Fill(hnoBS,"yrecBeamvsdyXBS",recVtxs->begin()->yError(),recVtxs->begin()->y()-vertexBeamSpot_.y0());
if(printXBS_) {
cout << Form("XBS %10d %5d %10d %4d %lu %5.2f %8f %8f %8f %8f %8f %8f",
run_,luminosityBlock_,event_,bunchCrossing_,
(unsigned long)(recVtxs->begin()->tracksSize()), recVtxs->begin()->ndof(),
recVtxs->begin()->x(), recVtxs->begin()->xError(),
recVtxs->begin()->y(), recVtxs->begin()->yError(),
recVtxs->begin()->z(), recVtxs->begin()->zError()
) << endl;
}
}
}
// for the associator
Handle<View<Track> > trackCollectionH;
iEvent.getByLabel(recoTrackProducer_,trackCollectionH);
Handle<HepMCProduct> evtMC;
Handle<SimVertexContainer> simVtxs;
iEvent.getByLabel( simG4_, simVtxs);
Handle<SimTrackContainer> simTrks;
iEvent.getByLabel( simG4_, simTrks);
edm::Handle<TrackingParticleCollection> TPCollectionH ;
edm::Handle<TrackingVertexCollection> TVCollectionH ;
bool gotTP=iEvent.getByLabel("mergedtruth","MergedTrackTruth",TPCollectionH);
bool gotTV=iEvent.getByLabel("mergedtruth","MergedTrackTruth",TVCollectionH);
iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder",theB_);
fBfield_=((*theB_).field()->inTesla(GlobalPoint(0.,0.,0.))).z();
vector<SimEvent> simEvt;
if (gotTP && gotTV ){
edm::ESHandle<TrackAssociatorBase> theHitsAssociator;
iSetup.get<TrackAssociatorRecord>().get("TrackAssociatorByHits",theHitsAssociator);
associatorByHits_ = (TrackAssociatorBase *) theHitsAssociator.product();
r2s_ = associatorByHits_->associateRecoToSim (trackCollectionH,TPCollectionH, &iEvent );
//simEvt=getSimEvents(iEvent, TPCollectionH, TVCollectionH, trackCollectionH);
simEvt=getSimEvents(TPCollectionH, TVCollectionH, trackCollectionH);
if (simEvt.size()==0){
cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
cout << " !!!!!!!!!!!!!!!!!!!!!! got TrackingParticles but no simEvt !!!!!!!!!!!!!!!!!" << endl;
cout << " dumping Tracking particles " << endl;
const TrackingParticleCollection* simTracks = TPCollectionH.product();
for(TrackingParticleCollection::const_iterator it=simTracks->begin(); it!=simTracks->end(); it++){
cout << *it << endl;
}
cout << " dumping Tracking Vertices " << endl;
const TrackingVertexCollection* tpVtxs = TVCollectionH.product();
for(TrackingVertexCollection::const_iterator iv=tpVtxs->begin(); iv!=tpVtxs->end(); iv++){
cout << *iv << endl;
}
if(iEvent.getByLabel(mcproduct,evtMC)){
cout << "dumping simTracks" << endl;
for(SimTrackContainer::const_iterator t=simTrks->begin(); t!=simTrks->end(); ++t){
cout << *t << endl;
}
cout << "dumping simVertexes" << endl;
for(SimVertexContainer::const_iterator vv=simVtxs->begin();
vv!=simVtxs->end();
++vv){
cout << *vv << endl;
}
}else{
cout << "no hepMC" << endl;
}
cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
const HepMC::GenEvent* evt=evtMC->GetEvent();
if(evt){
std::cout << "process id " <<evt->signal_process_id()<<std::endl;
std::cout <<"signal process vertex "<< ( evt->signal_process_vertex() ?
evt->signal_process_vertex()->barcode() : 0 ) <<std::endl;
std::cout <<"number of vertices " << evt->vertices_size() << std::endl;
evt->print();
}else{
cout << "no event in HepMC" << endl;
}
cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
}
}
if(gotTV){
if(verbose_){
cout << "Found Tracking Vertices " << endl;
}
simpv=getSimPVs(TVCollectionH);
}else if(iEvent.getByLabel(mcproduct,evtMC)){
simpv=getSimPVs(evtMC);
if(verbose_){
cout << "Using HepMCProduct " << endl;
std::cout << "simtrks " << simTrks->size() << std::endl;
}
tsim = PrimaryVertexAnalyzer4PU::getSimTrkParameters(simTrks, simVtxs, simUnit_);
}else{
// if(verbose_({cout << "No MC info at all" << endl;}
}
// get the beam spot from the appropriate dummy vertex, if available
for(reco::VertexCollection::const_iterator v=recVtxs->begin();
v!=recVtxs->end(); ++v){
if ((v->ndof()==-3) && (v->chi2()==0) ){
myBeamSpot=math::XYZPoint(v->x(), v->y(), v->z());
}
}
hsimPV["nsimvtx"]->Fill(simpv.size());
for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
vsim!=simpv.end(); vsim++){
if(doMatching_){
matchRecTracksToVertex(*vsim, tsim, recTrks); // hepmc, based on track parameters
}
hsimPV["nbsimtksinvtx"]->Fill(vsim->nGenTrk);
hsimPV["xsim"]->Fill(vsim->x*simUnit_);
hsimPV["ysim"]->Fill(vsim->y*simUnit_);
hsimPV["zsim"]->Fill(vsim->z*simUnit_);
hsimPV["xsim1"]->Fill(vsim->x*simUnit_);
hsimPV["ysim1"]->Fill(vsim->y*simUnit_);
hsimPV["zsim1"]->Fill(vsim->z*simUnit_);
Fill(hsimPV,"xsim2",vsim->x*simUnit_,vsim==simpv.begin());
Fill(hsimPV,"ysim2",vsim->y*simUnit_,vsim==simpv.begin());
Fill(hsimPV,"zsim2",vsim->z*simUnit_,vsim==simpv.begin());
hsimPV["xsim2"]->Fill(vsim->x*simUnit_);
hsimPV["ysim2"]->Fill(vsim->y*simUnit_);
hsimPV["zsim2"]->Fill(vsim->z*simUnit_);
hsimPV["xsim3"]->Fill(vsim->x*simUnit_);
hsimPV["ysim3"]->Fill(vsim->y*simUnit_);
hsimPV["zsim3"]->Fill(vsim->z*simUnit_);
hsimPV["xsimb"]->Fill(vsim->x*simUnit_-vertexBeamSpot_.x0());
hsimPV["ysimb"]->Fill(vsim->y*simUnit_-vertexBeamSpot_.y0());
hsimPV["zsimb"]->Fill(vsim->z*simUnit_-vertexBeamSpot_.z0());
hsimPV["xsimb1"]->Fill(vsim->x*simUnit_-vertexBeamSpot_.x0());
hsimPV["ysimb1"]->Fill(vsim->y*simUnit_-vertexBeamSpot_.y0());
hsimPV["zsimb1"]->Fill(vsim->z*simUnit_-vertexBeamSpot_.z0());
}
if(bnoBS){
analyzeVertexCollection(hnoBS, recVtxs, recTrks, simpv,"noBS");
analyzeVertexCollectionTP(hnoBS, recVtxs, recTrks, simEvt,"noBS");
}
if(bBS){
analyzeVertexCollection(hBS, recVtxsBS, recTrks, simpv,"BS");
analyzeVertexCollectionTP(hBS, recVtxsBS, recTrks, simEvt,"BS");
}
if(bDA){
analyzeVertexCollection(hDA, recVtxsDA, recTrks, simpv,"DA");
analyzeVertexCollectionTP(hDA, recVtxsDA, recTrks, simEvt,"DA");
}
// if(bPIX){
// analyzeVertexCollection(hPIX, recVtxsPIX, recTrks, simpv,"PIX");
// analyzeVertexCollectionTP(hPIX, recVtxsPIX, recTrks, simEvt,"PIX");
// }
if((dumpThisEvent_&& (dumpcounter_<100)) ||(verbose_ && (eventcounter_<ndump_))){
cout << endl << "Event dump" << endl
<< "event counter=" << eventcounter_
<< " Run=" << run_ << " LumiBlock " << luminosityBlock_ << " event " << event_
<< " bx=" << bunchCrossing_ << " orbit=" << orbitNumber_
<< std::endl;
dumpcounter_++;
//evtMC->GetEvent()->print();
//printRecTrks(recTrks); // very verbose !!
// if (bPIX) printRecVtxs(recVtxsPIX,"pixel vertices");
if (bnoBS) {printRecVtxs(recVtxs,"Offline without Beamspot");}
if (bnoBS && (!bDA)){ printPVTrks(recTrks, recVtxs, tsim, simEvt, false);}
if (bBS) printRecVtxs(recVtxsBS,"Offline with Beamspot");
if (bDA) {
printRecVtxs(recVtxsDA,"Offline DA");
printPVTrks(recTrks, recVtxsDA, tsim, simEvt, false);
}
if (dumpcounter_<2){cout << "beamspot " << vertexBeamSpot_ << endl;}
}
if(verbose_){
std::cout << std::endl;
}
}
| void PrimaryVertexAnalyzer4PU::analyzeVertexCollection | ( | std::map< std::string, TH1 * > & | h, |
| const edm::Handle< reco::VertexCollection > | recVtxs, | ||
| const edm::Handle< reco::TrackCollection > | recTrks, | ||
| std::vector< simPrimaryVertex > & | simpv, | ||
| const std::string | message = "" |
||
| ) | [private] |
Definition at line 2758 of file PrimaryVertexAnalyzer4PU.cc.
References reco::BeamSpot::BeamWidthX(), reco::BeamSpot::BeamWidthY(), bunchCrossing_, ChiSquaredProbability(), gather_cfg::cout, data, DEBUG_, dumpPUcandidates_, dumpThisEvent_, event_, eventcounter_, Fill(), fillTrackHistos(), i, getHLTprescales::index, edm::detail::isnan(), edm::HandleBase::isValid(), j, create_public_lumi_plots::log, m, M_PI, nt, NULL, position, funct::pow(), run_, reco::TransientTrack::setBeamSpot(), simUnit_, mathSSE::sqrt(), lumiQTWidget::t, theB_, theTrackFilter, reco::Vertex::tracksSize(), groupFilesInBlocks::tt, v, verbose_, vertexBeamSpot_, w(), x, reco::BeamSpot::x0(), reco::BeamSpot::y0(), z, reco::BeamSpot::z0(), and zmatch_.
Referenced by analyze().
{
//cout <<"PrimaryVertexAnalyzer4PU::analyzeVertexCollection (HepMC), simpvs=" << simpv.size() << endl;
int nrectrks=recTrks->size();
int nrecvtxs=recVtxs->size();
int nseltrks=-1;
reco::TrackCollection selTrks; // selected tracks
reco::TrackCollection lostTrks; // selected but lost tracks (part of dropped clusters)
// extract dummy vertices representing clusters
reco::VertexCollection clusters;
reco::Vertex allSelected;
double cpufit=0;
double cpuclu=0;
for(reco::VertexCollection::const_iterator v=recVtxs->begin();
v!=recVtxs->end(); ++v){
if ( (fabs(v->ndof()+3.)<0.0001) && (v->chi2()<=0) ){
// this dummy vertex is for the full event
allSelected=(*v);
nseltrks=(allSelected.tracksSize());
nrecvtxs--;
cpuclu=-v->chi2();
continue;
}else if( (fabs(v->ndof()+2.)<0.0001) && (v->chi2()==0) ){
// this is a cluster, not a vertex
clusters.push_back(*v);
Fill(h,"cpuvsntrk",(double) v->tracksSize(),fabs(v->y()));
cpufit+=fabs(v->y());
Fill(h,"nclutrkall",(double) v->tracksSize());
Fill(h,"selstat",v->x());
//Fill(h,"nclutrkvtx",);// see below
nrecvtxs--;
}
}
Fill(h,"cpuclu",cpuclu);
Fill(h,"cpufit",cpufit);
Fill(h,"cpucluvsntrk",nrectrks, cpuclu);
if(simpv.size()>0){//this is mc
double dsimrecx=0.;
double dsimrecy=0.;//0.0011;
double dsimrecz=0.;//0.0012;
// vertex matching and efficiency bookkeeping
int nsimtrk=0;
for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
vsim!=simpv.end(); vsim++){
nsimtrk+=vsim->nGenTrk;
// look for a matching reconstructed vertex
vsim->recVtx=NULL;
vsim->cluster=-1;
for(reco::VertexCollection::const_iterator vrec=recVtxs->begin(); vrec!=recVtxs->end(); ++vrec){
if( vrec->isFake() ) {
continue; // skip fake vertices (=beamspot)
cout << "fake vertex" << endl;
}
if( vrec->ndof()<0. )continue; // skip dummy clusters, if any
// 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-dsimrecz)>fabs(vrec->z()-vsim->z-dsimrecz)))
|| (!vsim->recVtx) )
{
vsim->recVtx=&(*vrec);
// find the corresponding cluster
for(unsigned int iclu=0; iclu<clusters.size(); iclu++){
if( fabs(clusters[iclu].position().z()-vrec->position().z()) < 0.001 ){
vsim->cluster=iclu;
vsim->nclutrk=clusters[iclu].position().y();
}
}
}
// the following only works in MC samples without pile-up
if ((simpv.size()==1) && ( fabs(vsim->recVtx->position().z()-vsim->z-dsimrecz)>zmatch_ )){
// now we have a recvertex without a matching simvertex, I would call this fake
// however, the G4 info does not contain pile-up
Fill(h,"fakeVtxZ",vrec->z());
if (vrec->ndof()>=0.5) Fill(h,"fakeVtxZNdofgt05",vrec->z());
if (vrec->ndof()>=2.0) Fill(h,"fakeVtxZNdofgt2",vrec->z());
Fill(h,"fakeVtxNdof",vrec->ndof());
//Fill(h,"fakeVtxNdof1",vrec->ndof());
Fill(h,"fakeVtxNtrk",vrec->tracksSize());
if(vrec->tracksSize()==2){ Fill(h,"fake2trkchi2vsndof",vrec->ndof(),vrec->chi2()); }
if(vrec->tracksSize()==3){ Fill(h,"fake3trkchi2vsndof",vrec->ndof(),vrec->chi2()); }
if(vrec->tracksSize()==4){ Fill(h,"fake4trkchi2vsndof",vrec->ndof(),vrec->chi2()); }
if(vrec->tracksSize()==5){ Fill(h,"fake5trkchi2vsndof",vrec->ndof(),vrec->chi2()); }
}
}
Fill(h,"nsimtrk",float(nsimtrk));
Fill(h,"nsimtrk",float(nsimtrk),vsim==simpv.begin());
Fill(h,"nrecsimtrk",float(vsim->nMatchedTracks));
Fill(h,"nrecnosimtrk",float(nsimtrk-vsim->nMatchedTracks));
// histogram properties of matched vertices
if (vsim->recVtx && ( fabs(vsim->recVtx->z()-vsim->z*simUnit_)<zmatch_ )){
if(verbose_){std::cout <<"primary matched " << message << " " << setw(8) << setprecision(4) << vsim->x << " " << vsim->y << " " << vsim->z << std:: endl;}
Fill(h,"matchedVtxNdof", vsim->recVtx->ndof());
// residuals an pulls with respect to simulated vertex
Fill(h,"resx", vsim->recVtx->x()-vsim->x*simUnit_ );
Fill(h,"resy", vsim->recVtx->y()-vsim->y*simUnit_ );
Fill(h,"resz", vsim->recVtx->z()-vsim->z*simUnit_ );
Fill(h,"resz10", vsim->recVtx->z()-vsim->z*simUnit_ );
Fill(h,"pullx", (vsim->recVtx->x()-vsim->x*simUnit_)/vsim->recVtx->xError() );
Fill(h,"pully", (vsim->recVtx->y()-vsim->y*simUnit_)/vsim->recVtx->yError() );
Fill(h,"pullz", (vsim->recVtx->z()-vsim->z*simUnit_)/vsim->recVtx->zError() );
Fill(h,"resxr", vsim->recVtx->x()-vsim->x*simUnit_-dsimrecx);
Fill(h,"resyr", vsim->recVtx->y()-vsim->y*simUnit_-dsimrecy );
Fill(h,"reszr", vsim->recVtx->z()-vsim->z*simUnit_-dsimrecz);
Fill(h,"pullxr", (vsim->recVtx->x()-vsim->x*simUnit_-dsimrecx)/vsim->recVtx->xError() );
Fill(h,"pullyr", (vsim->recVtx->y()-vsim->y*simUnit_-dsimrecy)/vsim->recVtx->yError() );
Fill(h,"pullzr", (vsim->recVtx->z()-vsim->z*simUnit_-dsimrecz)/vsim->recVtx->zError() );
// efficiency with zmatch within 500 um (or whatever zmatch is)
Fill(h,"eff", 1.);
if(simpv.size()==1){
if (vsim->recVtx==&(*recVtxs->begin())){
Fill(h,"efftag", 1.);
}else{
Fill(h,"efftag", 0.);
cout << "signal vertex not tagged " << message << " " << eventcounter_ << endl;
// call XXClusterizerInZ.vertices(seltrks,3)
}
}
Fill(h,"effvsptsq",vsim->ptsq,1.);
Fill(h,"effvsnsimtrk",vsim->nGenTrk,1.);
Fill(h,"effvsnrectrk",nrectrks,1.);
Fill(h,"effvsnseltrk",nseltrks,1.);
Fill(h,"effvsz",vsim->z*simUnit_,1.);
Fill(h,"effvsz2",vsim->z*simUnit_,1.);
Fill(h,"effvsr",sqrt(vsim->x*vsim->x+vsim->y*vsim->y)*simUnit_,1.);
}else{ // no matching rec vertex found for this simvertex
bool plapper=verbose_ && vsim->nGenTrk;
if(plapper){
// be quiet about invisble vertices
std::cout << "primary not found " << message << " " << eventcounter_ << " x=" <<vsim->x*simUnit_ << " y=" << vsim->y*simUnit_ << " z=" << vsim->z*simUnit_ << " nGenTrk=" << vsim->nGenTrk << std::endl;
}
int mistype=0;
if (vsim->recVtx){
if(plapper){
std::cout << "nearest recvertex at " << vsim->recVtx->z() << " dz=" << vsim->recVtx->z()-vsim->z*simUnit_ << std::endl;
}
if (fabs(vsim->recVtx->z()-vsim->z*simUnit_)<0.2 ){
Fill(h,"effvsz2",vsim->z*simUnit_,1.);
}
if (fabs(vsim->recVtx->z()-vsim->z*simUnit_)<0.5 ){
if(verbose_){std::cout << "type 1, lousy z vertex" << std::endl;}
Fill(h,"zlost1", vsim->z*simUnit_,1.);
mistype=1;
}else{
if(plapper){std::cout << "type 2a no vertex anywhere near" << std::endl;}
mistype=2;
}
}else{// no recVtx at all
mistype=2;
if(plapper){std::cout << "type 2b, no vertex at all" << std::endl;}
}
if(mistype==2){
int selstat=-3;
// no matching vertex found, is there a cluster?
for(unsigned int iclu=0; iclu<clusters.size(); iclu++){
if( fabs(clusters[iclu].position().z()-vsim->z*simUnit_) < 0.1 ){
selstat=int(clusters[iclu].position().x()+0.1);
if(verbose_){std::cout << "matching cluster found with selstat=" << clusters[iclu].position().x() << std::endl;}
}
}
if (selstat==0){
if(plapper){std::cout << "vertex rejected (distance to beam)" << std::endl;}
Fill(h,"zlost3", vsim->z*simUnit_,1.);
}else if(selstat==-1){
if(plapper) {std::cout << "vertex invalid" << std::endl;}
Fill(h,"zlost4", vsim->z*simUnit_,1.);
}else if(selstat==1){
if(plapper){std::cout << "vertex accepted, this cannot be right!!!!!!!!!!" << std::endl;}
}else if(selstat==-2){
if(plapper){std::cout << "dont know what this means !!!!!!!!!!" << std::endl;}
}else if(selstat==-3){
if(plapper){std::cout << "no matching cluster found " << std::endl;}
Fill(h,"zlost2", vsim->z*simUnit_,1.);
}else{
if(plapper){std::cout << "dont know what this means either !!!!!!!!!!" << selstat << std::endl;}
}
}//
Fill(h,"eff", 0.);
if(simpv.size()==1){ Fill(h,"efftag", 0.); }
Fill(h,"effvsptsq",vsim->ptsq,0.);
Fill(h,"effvsnsimtrk",float(vsim->nGenTrk),0.);
Fill(h,"effvsnrectrk",nrectrks,0.);
Fill(h,"effvsnseltrk",nseltrks,0.);
Fill(h,"effvsz",vsim->z*simUnit_,0.);
Fill(h,"effvsr",sqrt(vsim->x*vsim->x+vsim->y*vsim->y)*simUnit_,0.);
} // no recvertex for this simvertex
}
// end of sim/rec matching
// purity of event vertex tags
if (recVtxs->size()>0){
Double_t dz=(*recVtxs->begin()).z() - (*simpv.begin()).z*simUnit_;
Fill(h,"zdistancetag",dz);
Fill(h,"abszdistancetag",fabs(dz));
if( fabs(dz)<zmatch_){
Fill(h,"puritytag",1.);
}else{
// bad tag: the true primary was more than 500 um (or zmatch) away from the tagged primary
Fill(h,"puritytag",0.);
}
}
}else{
//if(verbose_) cout << "PrimaryVertexAnalyzer4PU::analyzeVertexCollection: simPV is empty!" << endl;
}
//******* the following code does not require MC and will/should work for data **********
Fill(h,"bunchCrossing",bunchCrossing_);
if(recTrks->size()>0) Fill(h,"bunchCrossingLogNtk",bunchCrossing_,log(recTrks->size())/log(10.));
// ----------------- reconstructed tracks ------------------------
// the list of selected tracks can only be correct if the selection parameterset and track collection
// is the same that was used for the reconstruction
int nt=0;
for(reco::TrackCollection::const_iterator t=recTrks->begin();
t!=recTrks->end(); ++t){
if((recVtxs->size()>0) && (recVtxs->begin()->isValid())){
fillTrackHistos(h,"all",*t,&(*recVtxs->begin()));
}else{
fillTrackHistos(h,"all",*t);
}
if(recTrks->size()>100) fillTrackHistos(h,"M",*t);
TransientTrack tt = theB_->build(&(*t)); tt.setBeamSpot(vertexBeamSpot_);
if (theTrackFilter(tt)){
selTrks.push_back(*t);
fillTrackHistos(h,"sel",*t);
int foundinvtx=0;
int nvtemp=-1;
for(reco::VertexCollection::const_iterator v=recVtxs->begin();
v!=recVtxs->end(); ++v){
nvtemp++;
if(( v->isFake()) || (v->ndof()<-2) ) break;
for(trackit_t tv=v->tracks_begin(); tv!=v->tracks_end(); tv++ ){
if( ((**tv).vz()==t->vz()&&((**tv).phi()==t->phi())) ) {
foundinvtx++;
}
}
}
if(foundinvtx==0){
fillTrackHistos(h,"sellost",*t);
}else if(foundinvtx>1){
cout << "hmmmm " << foundinvtx << endl;
}
}
nt++;
}
if (nseltrks<0){
nseltrks=selTrks.size();
}else if( ! (nseltrks==(int)selTrks.size()) ){
std::cout << "Warning: inconsistent track selection !" << std::endl;
}
// fill track histograms of vertex tracks
int nrec=0, nrec0=0, nrec8=0, nrec2=0, nrec4=0;
for(reco::VertexCollection::const_iterator v=recVtxs->begin();
v!=recVtxs->end(); ++v){
if (! (v->isFake()) && v->ndof()>0 && v->chi2()>0 ){
nrec++;
if (v->ndof()>0) nrec0++;
if (v->ndof()>8) nrec8++;
if (v->ndof()>2) nrec2++;
if (v->ndof()>4) nrec4++;
for(trackit_t t=v->tracks_begin(); t!=v->tracks_end(); t++){
if(v==recVtxs->begin()){
fillTrackHistos(h,"tagged",**t, &(*v));
}else{
fillTrackHistos(h,"untagged",**t, &(*v));
}
Float_t wt=v->trackWeight(*t);
//dumpHitInfo(**t); cout << " w=" << wt << endl;
Fill(h,"trackWt",wt);
if(wt>0.5){
fillTrackHistos(h,"wgt05",**t, &(*v));
if(v->ndof()>4) fillTrackHistos(h,"ndof4",**t, &(*v));
}else{
fillTrackHistos(h,"wlt05",**t, &(*v));
}
}
}
}
// bachelor tracks (only available through clusters right now)
for(unsigned int iclu=0; iclu<clusters.size(); iclu++){
if (clusters[iclu].tracksSize()==1){
for(trackit_t t = clusters[iclu].tracks_begin();
t!=clusters[iclu].tracks_end(); t++){
fillTrackHistos(h,"bachelor",**t);
}
}
}
// ----------------- reconstructed vertices ------------------------
// event
Fill(h,"szRecVtx",recVtxs->size());
Fill(h,"nclu",clusters.size());
Fill(h,"nseltrk",nseltrks);
Fill(h,"nrectrk",nrectrks);
Fill(h,"nrecvtx",nrec);
Fill(h,"nrecvtx2",nrec2);
Fill(h,"nrecvtx4",nrec4);
Fill(h,"nrecvtx8",nrec8);
if(nrec>0){
Fill(h,"eff0vsntrec",nrectrks,1.);
Fill(h,"eff0vsntsel",nseltrks,1.);
}else{
Fill(h,"eff0vsntrec",nrectrks,0.);
Fill(h,"eff0vsntsel",nseltrks,0.);
if((nseltrks>1)&&(verbose_)){
cout << Form("PrimaryVertexAnalyzer4PU: %s may have lost a vertex %10d %10d %4d / %4d ",message.c_str(),run_, event_, nrectrks,nseltrks) << endl;
dumpThisEvent_=true;
}
}
if(nrec0>0) { Fill(h,"eff0ndof0vsntsel",nseltrks,1.);}else{ Fill(h,"eff0ndof0vsntsel",nseltrks,0.);}
if(nrec2>0) { Fill(h,"eff0ndof2vsntsel",nseltrks,1.);}else{ Fill(h,"eff0ndof2vsntsel",nseltrks,0.);}
if(nrec4>0) { Fill(h,"eff0ndof4vsntsel",nseltrks,1.);}else{ Fill(h,"eff0ndof4vsntsel",nseltrks,0.);}
if(nrec8>0) { Fill(h,"eff0ndof8vsntsel",nseltrks,1.);}else{ Fill(h,"eff0ndof8vsntsel",nseltrks,0.);}
if((nrec>1)&&(DEBUG_)) {
cout << "multivertex event" << endl;
dumpThisEvent_=true;
}
if((nrectrks>10)&&(nseltrks<3)){
cout << "small fraction of selected tracks " << endl;
dumpThisEvent_=true;
}
// properties of events without a vertex
if((nrec==0)||(recVtxs->begin()->isFake())){
Fill(h,"nrectrk0vtx",nrectrks);
Fill(h,"nseltrk0vtx",nseltrks);
Fill(h,"nclu0vtx",clusters.size());
}
// properties of (valid) vertices
double ndof2=-10,ndof1=-10, zndof1=0, zndof2=0;
for(reco::VertexCollection::const_iterator v=recVtxs->begin();
v!=recVtxs->end(); ++v){
if(v->isFake()){ Fill(h,"isFake",1.);}else{ Fill(h,"isFake",0.);}
if(v->isFake()||((v->ndof()<0)&&(v->ndof()>-3))){ Fill(h,"isFake1",1.);}else{ Fill(h,"isFake1",0.);}
if((v->isFake())||(v->ndof()<0)) continue;
if(v->ndof()>ndof1){ ndof2=ndof1; zndof2=zndof1; ndof1=v->ndof(); zndof1=v->position().z();}
else if(v->ndof()>ndof2){ ndof2=v->ndof(); zndof2=v->position().z();}
// some special histogram for two track vertices
if(v->tracksSize()==2){
const TrackBaseRef& t1= *(v->tracks_begin());
const TrackBaseRef& t2=*(v->tracks_begin()+1);
bool os=(t1->charge()*t2->charge()<0);
double dphi=t1->phi()-t2->phi(); if (dphi<0) dphi+=2*M_PI;
double m12=sqrt(pow( sqrt(pow(0.139,2)+pow( t1->p(),2)) +sqrt(pow(0.139,2)+pow( t2->p(),2)) ,2)
-pow(t1->px()+t2->px(),2)
-pow(t1->py()+t2->py(),2)
-pow(t1->pz()+t2->pz(),2)
);
if(os){
Fill(h,"2trkdetaOS",t1->eta()-t2->eta());
Fill(h,"2trkmassOS",m12);
}else{
Fill(h,"2trkdetaSS",t1->eta()-t2->eta());
Fill(h,"2trkmassSS",m12);
}
Fill(h,"2trkdphi",dphi);
Fill(h,"2trkseta",t1->eta()+t2->eta());
if(fabs(dphi-M_PI)<0.1) Fill(h,"2trksetacurl",t1->eta()+t2->eta());
if(fabs(t1->eta()+t2->eta())<0.1) Fill(h,"2trkdphicurl",dphi);
// fill separately for extra vertices
if(v!=recVtxs->begin()){
if(os){
Fill(h,"2trkdetaOSPU",t1->eta()-t2->eta());
Fill(h,"2trkmassOSPU",m12);
}else{
Fill(h,"2trkdetaSSPU",t1->eta()-t2->eta());
Fill(h,"2trkmassSSPU",m12);
}
Fill(h,"2trkdphiPU",dphi);
Fill(h,"2trksetaPU",t1->eta()+t2->eta());
if(fabs(dphi-M_PI)<0.1) Fill(h,"2trksetacurlPU",t1->eta()+t2->eta());
if(fabs(t1->eta()+t2->eta())<0.1) Fill(h,"2trkdphicurlPU",dphi);
}
}
Fill(h,"trkchi2vsndof",v->ndof(),v->chi2());
if(v->ndof()>0){ Fill(h,"trkchi2overndof",v->chi2()/v->ndof()); }
if(v->tracksSize()==2){ Fill(h,"2trkchi2vsndof",v->ndof(),v->chi2()); }
if(v->tracksSize()==3){ Fill(h,"3trkchi2vsndof",v->ndof(),v->chi2()); }
if(v->tracksSize()==4){ Fill(h,"4trkchi2vsndof",v->ndof(),v->chi2()); }
if(v->tracksSize()==5){ Fill(h,"5trkchi2vsndof",v->ndof(),v->chi2()); }
Fill(h,"nbtksinvtx",v->tracksSize());
Fill(h,"nbtksinvtx2",v->tracksSize());
Fill(h,"vtxchi2",v->chi2());
Fill(h,"vtxndf",v->ndof());
Fill(h,"vtxprob",ChiSquaredProbability(v->chi2() ,v->ndof()));
Fill(h,"vtxndfvsntk",v->tracksSize(), v->ndof());
Fill(h,"vtxndfoverntk",v->ndof()/v->tracksSize());
Fill(h,"vtxndf2overntk",(v->ndof()+2)/v->tracksSize());
Fill(h,"zrecvsnt",v->position().z(),float(nrectrks));
if(nrectrks>100){
Fill(h,"zrecNt100",v->position().z());
}
if(v->ndof()>2.0){ // enter only vertices that really contain tracks
Fill(h,"xrec",v->position().x());
Fill(h,"yrec",v->position().y());
Fill(h,"zrec",v->position().z());
Fill(h,"xrec1",v->position().x());
Fill(h,"yrec1",v->position().y());
Fill(h,"zrec1",v->position().z());
Fill(h,"xrec2",v->position().x());
Fill(h,"yrec2",v->position().y());
Fill(h,"zrec2",v->position().z());
Fill(h,"xrec3",v->position().x());
Fill(h,"yrec3",v->position().y());
Fill(h,"zrec3",v->position().z());
Fill(h,"xrecb",v->position().x()-vertexBeamSpot_.x0());
Fill(h,"yrecb",v->position().y()-vertexBeamSpot_.y0());
Fill(h,"zrecb",v->position().z()-vertexBeamSpot_.z0());
Fill(h,"xrecBeam",v->position().x()-vertexBeamSpot_.x0());
Fill(h,"yrecBeam",v->position().y()-vertexBeamSpot_.y0());
Fill(h,"zrecBeam",v->position().z()-vertexBeamSpot_.z0());
Fill(h,"xrecBeamPull",(v->position().x()-vertexBeamSpot_.x0())/sqrt(pow(v->xError(),2)+pow(vertexBeamSpot_.BeamWidthX(),2)));
Fill(h,"yrecBeamPull",(v->position().y()-vertexBeamSpot_.y0())/sqrt(pow(v->yError(),2)+pow(vertexBeamSpot_.BeamWidthY(),2)));
Fill(h,"xrecBeamvsdx",v->xError(),v->position().x()-vertexBeamSpot_.x0());
Fill(h,"yrecBeamvsdy",v->yError(),v->position().y()-vertexBeamSpot_.y0());
Fill(h,"xrecBeamvsdxR2",v->position().x()-vertexBeamSpot_.x0(),v->xError());
Fill(h,"yrecBeamvsdyR2",v->position().y()-vertexBeamSpot_.y0(),v->yError());
Fill(h,"xrecBeam2vsdx2prof",pow(v->xError(),2),pow(v->position().x()-vertexBeamSpot_.x0(),2));
Fill(h,"yrecBeam2vsdy2prof",pow(v->yError(),2),pow(v->position().y()-vertexBeamSpot_.y0(),2));
Fill(h,"xrecBeamvsdx2",pow(v->xError(),2),pow(v->position().x()-vertexBeamSpot_.x0(),2));
Fill(h,"yrecBeamvsdy2",pow(v->yError(),2),pow(v->position().y()-vertexBeamSpot_.y0(),2));
Fill(h,"xrecBeamvsz",v->position().z(),v->position().x()-vertexBeamSpot_.x0());
Fill(h,"yrecBeamvsz",v->position().z(),v->position().y()-vertexBeamSpot_.y0());
Fill(h,"xrecBeamvszprof",v->position().z(),v->position().x()-vertexBeamSpot_.x0());
Fill(h,"yrecBeamvszprof",v->position().z(),v->position().y()-vertexBeamSpot_.y0());
Fill(h,"xrecBeamvsdxprof",v->xError(),v->position().x()-vertexBeamSpot_.x0());
Fill(h,"yrecBeamvsdyprof",v->yError(),v->position().y()-vertexBeamSpot_.y0());
Fill(h,"errx",v->xError());
Fill(h,"erry",v->yError());
Fill(h,"errz",v->zError());
double vxx=v->covariance(0,0);
double vyy=v->covariance(1,1);
double vxy=v->covariance(1,0);
double dv=0.25*(vxx+vyy)*(vxx+vyy-(vxx*vyy-vxy*vxy));
if(dv>0){
double l1=0.5*(vxx+vyy)+sqrt(dv);
Fill(h,"err1",sqrt(l1));
double l2=sqrt(0.5*(vxx+vyy)-sqrt(dv));
if(l2>0) Fill(h,"err2",sqrt(l2));
}
// look at the tagged vertex separately
if (v==recVtxs->begin()){
Fill(h,"nbtksinvtxTag",v->tracksSize());
Fill(h,"nbtksinvtxTag2",v->tracksSize());
Fill(h,"xrectag",v->position().x());
Fill(h,"yrectag",v->position().y());
Fill(h,"zrectag",v->position().z());
}else{
Fill(h,"nbtksinvtxPU",v->tracksSize());
Fill(h,"nbtksinvtxPU2",v->tracksSize());
}
// vertex resolution vs number of tracks
Fill(h,"xresvsntrk",v->tracksSize(),v->xError());
Fill(h,"yresvsntrk",v->tracksSize(),v->yError());
Fill(h,"zresvsntrk",v->tracksSize(),v->zError());
}
// cluster properties (if available)
for(unsigned int iclu=0; iclu<clusters.size(); iclu++){
if( fabs(clusters[iclu].position().z()-v->position().z()) < 0.0001 ){
Fill(h,"nclutrkvtx",clusters[iclu].tracksSize());
}
}
// properties of (valid) neighbour vertices
reco::VertexCollection::const_iterator v1=v; v1++;
for(; v1!=recVtxs->end(); ++v1){
if((v1->isFake())||(v1->ndof()<0)) continue;
Fill(h,"zdiffrec",v->position().z()-v1->position().z());
// if(fabs(v->position().z()-v1->position().z())>1){
// cout << message << " zdiffrec=" << v->position().z()-v1->position().z() << " " << v->ndof() << " " << v1->ndof() << endl;
// dumpThisEvent_=true;
// }
double z0=v->position().z()-vertexBeamSpot_.z0();
double z1=v1->position().z()-vertexBeamSpot_.z0();
Fill(h,"zPUcand",z0); Fill(h,"zPUcand",z1);
Fill(h,"ndofPUcand",v->ndof()); Fill(h,"ndofPUcand",v1->ndof());
Fill(h,"zdiffvsz",z1-z0,0.5*(z1+z0));
if ((v->ndof()>2) && (v1->ndof()>2)){
Fill(h,"zdiffrec2",v->position().z()-v1->position().z());
Fill(h,"zPUcand2",z0);
Fill(h,"zPUcand2",z1);
Fill(h,"ndofPUcand2",v->ndof());
Fill(h,"ndofPUcand2",v1->ndof());
Fill(h,"zvszrec2",z0, z1);
Fill(h,"pzvspz2",TMath::Freq(z0/2.16),TMath::Freq(z1/2.16) );
}
if ((v->ndof()>4) && (v1->ndof()>4)){
Fill(h,"zdiffvsz4",z1-z0,0.5*(z1+z0));
Fill(h,"zdiffrec4",v->position().z()-v1->position().z());
Fill(h,"zvszrec4",z0, z1);
Fill(h,"pzvspz4",TMath::Freq(z0/2.16),TMath::Freq(z1/2.16) );
//cout << "ndof4 pu-candidate " << run_ << " " << event_ << endl ;
if(fabs(z0-z1)>1.0){
Fill(h,"xbeamPUcand",v->position().x()-vertexBeamSpot_.x0());
Fill(h,"ybeamPUcand",v->position().y()-vertexBeamSpot_.y0());
Fill(h,"xbeamPullPUcand",(v->position().x()-vertexBeamSpot_.x0())/v->xError());
Fill(h,"ybeamPullPUcand",(v->position().y()-vertexBeamSpot_.y0())/v->yError());
//Fill(h,"sumwOverNtkPUcand",sumw/v->tracksSize());
//Fill(h,"sumwOverNtkPUcand",sumw/v1->tracksSize());
Fill(h,"ndofOverNtkPUcand",v->ndof()/v->tracksSize());
Fill(h,"ndofOverNtkPUcand",v1->ndof()/v1->tracksSize());
Fill(h,"xbeamPUcand",v1->position().x()-vertexBeamSpot_.x0());
Fill(h,"ybeamPUcand",v1->position().y()-vertexBeamSpot_.y0());
Fill(h,"xbeamPullPUcand",(v1->position().x()-vertexBeamSpot_.x0())/v1->xError());
Fill(h,"ybeamPullPUcand",(v1->position().y()-vertexBeamSpot_.y0())/v1->yError());
Fill(h,"zPUcand4",z0);
Fill(h,"zPUcand4",z1);
Fill(h,"ndofPUcand4",v->ndof());
Fill(h,"ndofPUcand4",v1->ndof());
for(trackit_t t=v->tracks_begin(); t!=v->tracks_end(); t++){ if(v->trackWeight(*t)>0.5) fillTrackHistos(h,"PUcand",**t, &(*v)); }
for(trackit_t t=v1->tracks_begin(); t!=v1->tracks_end(); t++){ if(v1->trackWeight(*t)>0.5) fillTrackHistos(h,"PUcand",**t, &(*v1)); }
}
}
if ((v->ndof()>4) && (v1->ndof()>2) && (v1->ndof()<4)){
for(trackit_t t=v1->tracks_begin(); t!=v1->tracks_end(); t++){ if(v1->trackWeight(*t)>0.5) fillTrackHistos(h,"PUfake",**t, &(*v1)); }
}
if ((v->ndof()>8) && (v1->ndof()>8)){
Fill(h,"zdiffrec8",v->position().z()-v1->position().z());
if(dumpPUcandidates_ && fabs(z0-z1)>1.0){
cout << "PU candidate " << run_ << " " << event_ << " " << message << " zdiff=" <<z0-z1 << endl;
// cerr << "PU candidate " << run_ << " "<< event_ << " " << message << endl;
dumpThisEvent_=true;
}
}
}
// test track links, use reconstructed vertices
bool problem = false;
Fill(h,"nans",1.,std::isnan(v->position().x())*1.);
Fill(h,"nans",2.,std::isnan(v->position().y())*1.);
Fill(h,"nans",3.,std::isnan(v->position().z())*1.);
int index = 3;
for (int i = 0; i != 3; i++) {
for (int j = i; j != 3; j++) {
index++;
Fill(h,"nans",index*1., std::isnan(v->covariance(i, j))*1.);
if (std::isnan(v->covariance(i, j))) problem = true;
// in addition, diagonal element must be positive
if (j == i && v->covariance(i, j) < 0) {
Fill(h,"nans",index*1., 1.);
problem = true;
}
}
}
try {
for(trackit_t t = v->tracks_begin();
t!=v->tracks_end(); t++) {
// illegal charge
if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
Fill(h,"tklinks",0.);
}
else {
Fill(h,"tklinks",1.);
}
}
} catch (...) {
// exception thrown when trying to use linked track
Fill(h,"tklinks",0.);
}
if (problem) {
// analyze track parameter covariance definiteness
double data[25];
try {
int itk = 0;
for(trackit_t 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;
}// catch()
}// if (problem)
} // vertex loop (v)
// 2nd highest ndof
if (ndof2>0){
Fill(h,"ndofnr2",ndof2);
if(fabs(zndof1-zndof2)>1) Fill(h,"ndofnr2d1cm",ndof2);
if(fabs(zndof1-zndof2)>2) Fill(h,"ndofnr2d2cm",ndof2);
}
}
| void PrimaryVertexAnalyzer4PU::analyzeVertexCollectionTP | ( | std::map< std::string, TH1 * > & | h, |
| const edm::Handle< reco::VertexCollection > | recVtxs, | ||
| const edm::Handle< reco::TrackCollection > | recTrks, | ||
| std::vector< SimEvent > & | simEvt, | ||
| const std::string | message = "" |
||
| ) | [private] |
Definition at line 2447 of file PrimaryVertexAnalyzer4PU.cc.
References gather_cfg::cout, eventcounter_, Fill(), getTc(), getTruthMatchedVertexTracks(), create_public_lumi_plots::log, n, printEventSummary(), mathSSE::sqrt(), PrimaryVertexAnalyzer4PU::SimEvent::tk, PrimaryVertexAnalyzer4PU::SimEvent::tkprim, reco::Vertex::tracks_begin(), reco::Vertex::tracks_end(), reco::Vertex::tracksSize(), v, reco::TrackBase::vz(), and zmatch_.
Referenced by analyze().
{
// cout <<"PrimaryVertexAnalyzer4PU::analyzeVertexCollectionTP size=" << simEvt.size() << endl;
if(simEvt.size()==0)return;
printEventSummary(h, recVtxs,recTrks,simEvt,message);
//const int iSignal=0;
EncodedEventId iSignal=simEvt[0].eventId;
Fill(h,"npu0",simEvt.size());
for(vector<SimEvent>::iterator ev=simEvt.begin(); ev!=simEvt.end(); ev++){
Fill(h,"Tc", ev->Tc, ev==simEvt.begin());
Fill(h,"Chisq", ev->chisq, ev==simEvt.begin());
if(ev->chisq>0)Fill(h,"logChisq", log(ev->chisq), ev==simEvt.begin());
Fill(h,"dzmax", ev->dzmax, ev==simEvt.begin());
Fill(h,"dztrim",ev->dztrim,ev==simEvt.begin());
Fill(h,"m4m2", ev->m4m2, ev==simEvt.begin());
if(ev->Tc>0){ Fill(h,"logTc",log(ev->Tc)/log(10.),ev==simEvt.begin());}
for(vector<SimEvent>::iterator ev2=ev+1; ev2!=simEvt.end(); ev2++){
vector<TransientTrack> xt;
if((ev->tkprimsel.size()>0)&&(ev2->tkprimsel.size()>0)&&(ev->tkprimsel.size()+ev2->tkprimsel.size())>1){
xt.insert (xt.end() ,ev->tkprimsel.begin(),ev->tkprimsel.end());
xt.insert (xt.end() ,ev2->tkprimsel.begin(),ev2->tkprimsel.end());
double xTc,xChsq,xDzmax,xDztrim,xm4m2;
getTc(xt, xTc, xChsq, xDzmax, xDztrim,xm4m2);
if(xTc>0){
Fill(h,"xTc",xTc,ev==simEvt.begin());
Fill(h,"logxTc", log(xTc)/log(10),ev==simEvt.begin());
Fill(h,"xChisq", xChsq,ev==simEvt.begin());
if(xChsq>0){Fill(h,"logxChisq", log(xChsq),ev==simEvt.begin());};
Fill(h,"xdzmax", xDzmax,ev==simEvt.begin());
Fill(h,"xdztrim", xDztrim,ev==simEvt.begin());
Fill(h,"xm4m2", xm4m2,ev==simEvt.begin());
}
}
}
}
// --------------------------------------- count actual rec vtxs ----------------------
int nrecvtxs=0;//, nrecvtxs1=0, nrecvtxs2=0;
for(reco::VertexCollection::const_iterator v=recVtxs->begin();
v!=recVtxs->end(); ++v){
if ( (v->isFake()) || (v->ndof()<0) || (v->chi2()<=0) ) continue; // skip clusters
nrecvtxs++;
}
// --------------------------------------- fill the track assignment matrix ----------------------
for(vector<SimEvent>::iterator ev=simEvt.begin(); ev!=simEvt.end(); ev++){
ev->ntInRecVz.clear(); // just in case
ev->zmatch=-99.;
ev->nmatch=0;
for(reco::VertexCollection::const_iterator v=recVtxs->begin();
v!=recVtxs->end(); ++v){
double n=0, wt=0;
for(vector<TransientTrack>::iterator te=ev->tk.begin(); te!=ev->tk.end(); te++){
const reco::Track& RTe=te->track();
for(trackit_t tv=v->tracks_begin(); tv!=v->tracks_end(); tv++){
const reco::Track & RTv=*(tv->get());
if(RTe.vz()==RTv.vz()) {n++; wt+=v->trackWeight(*tv);}
}
}
ev->ntInRecVz[v->z()]=n;
if (n > ev->nmatch){ ev->nmatch=n; ev->zmatch=v->z(); ev->zmatch=v->z(); }
}
}
// call a vertex a fake if for every sim vertex there is another recvertex containing more tracks
// from that sim vertex than the current recvertex
double nfake=0;
for(reco::VertexCollection::const_iterator v=recVtxs->begin();
v!=recVtxs->end(); ++v){
bool matched=false;
for(vector<SimEvent>::iterator ev=simEvt.begin(); ev!=simEvt.end(); ev++){
if ((ev->nmatch>0)&&(ev->zmatch==v->z())){
matched=true;
}
}
if(!matched && !v->isFake()) {
nfake++;
cout << " fake rec vertex at z=" << v->z() << endl;
// some histograms of fake vertex properties here
Fill(h,"unmatchedVtxZ",v->z());
Fill(h,"unmatchedVtxNdof",v->ndof());
}
}
if(nrecvtxs>0){
Fill(h,"unmatchedVtx",nfake);
Fill(h,"unmatchedVtxFrac",nfake/nrecvtxs);
}
// --------------------------------------- match rec to sim ---------------------------------------
for(reco::VertexCollection::const_iterator v=recVtxs->begin();
v!=recVtxs->end(); ++v){
if ( (v->ndof()<0) || (v->chi2()<=0) ) continue; // skip clusters
double nmatch=-1; // highest number of tracks in recvtx v found in any event
EncodedEventId evmatch;
bool matchFound=false;
double nmatchvtx=0; // number of simvtcs contributing to recvtx v
for(vector<SimEvent>::iterator ev=simEvt.begin(); ev!=simEvt.end(); ev++){
double n=0; // number of tracks that are in both, the recvtx v and the event ev
for(vector<TransientTrack>::iterator te=ev->tk.begin(); te!=ev->tk.end(); te++){
const reco::Track& RTe=te->track();
for(trackit_t tv=v->tracks_begin(); tv!=v->tracks_end(); tv++){
const reco::Track & RTv=*(tv->get());
if(RTe.vz()==RTv.vz()){ n++;}
}
}
// find the best match in terms of the highest number of tracks
// from a simvertex in this rec vertex
if (n > nmatch){
nmatch=n;
evmatch=ev->eventId;
matchFound=true;
}
if(n>0){
nmatchvtx++;
}
}
double nmatchany=getTruthMatchedVertexTracks(*v).size();
if (matchFound && (nmatchany>0)){
// highest number of tracks in recvtx matched to (the same) sim vertex
// purity := -----------------------------------------------------------------
// number of truth matched tracks in this recvtx
double purity =nmatch/nmatchany;
Fill(h,"recmatchPurity",purity);
if(v==recVtxs->begin()){
Fill(h,"recmatchPurityTag",purity, (bool)(evmatch==iSignal));
}else{
Fill(h,"recmatchPuritynoTag",purity,(bool)(evmatch==iSignal));
}
}
Fill(h,"recmatchvtxs",nmatchvtx);
if(v==recVtxs->begin()){
Fill(h,"recmatchvtxsTag",nmatchvtx);
}else{
Fill(h,"recmatchvtxsnoTag",nmatchvtx);
}
} // recvtx loop
Fill(h,"nrecv",nrecvtxs);
// --------------------------------------- match sim to rec ---------------------------------------
int npu1=0, npu2=0;
for(vector<SimEvent>::iterator ev=simEvt.begin(); ev!=simEvt.end(); ev++){
if(ev->tk.size()>0) npu1++;
if(ev->tk.size()>1) npu2++;
bool isSignal= ev->eventId==iSignal;
Fill(h,"nRecTrkInSimVtx",(double) ev->tk.size(),isSignal);
Fill(h,"nPrimRecTrkInSimVtx",(double) ev->tkprim.size(),isSignal);
Fill(h,"sumpt2rec",sqrt(ev->sumpt2rec),isSignal);
Fill(h,"sumpt2",sqrt(ev->sumpt2),isSignal);
Fill(h,"sumpt",sqrt(ev->sumpt),isSignal);
double nRecVWithTrk=0; // vertices with tracks from this simvertex
double nmatch=0, ntmatch=0, zmatch=-99;
for(reco::VertexCollection::const_iterator v=recVtxs->begin();
v!=recVtxs->end(); ++v){
if ( (v->ndof()<-1) || (v->chi2()<=0) ) continue; // skip clusters
// count tracks found in both, sim and rec
double n=0, wt=0;
for(vector<TransientTrack>::iterator te=ev->tk.begin(); te!=ev->tk.end(); te++){
const reco::Track& RTe=te->track();
for(trackit_t tv=v->tracks_begin(); tv!=v->tracks_end(); tv++){
const reco::Track & RTv=*(tv->get());
if(RTe.vz()==RTv.vz()) {n++; wt+=v->trackWeight(*tv);}
}
}
if(n>0){ nRecVWithTrk++; }
if (n > nmatch){
nmatch=n; ntmatch=v->tracksSize(); zmatch=v->position().z();
}
}// end of reco vertex loop
// nmatch is the highest number of tracks from this sim vertex found in a single reco-vertex
if(ev->tk.size()>0){ Fill(h,"trkAssignmentEfficiency", nmatch/ev->tk.size(), isSignal); };
if(ev->tkprim.size()>0){ Fill(h,"primtrkAssignmentEfficiency", nmatch/ev->tkprim.size(), isSignal); };
// matched efficiency = efficiency for finding a reco vertex with > 50% of the simvertexs reconstructed tracks
double ntsim=ev->tk.size(); // may be better to use the number of primary tracks here ?
double matchpurity=nmatch/ntmatch;
if(ntsim>0){
Fill(h,"matchVtxFraction",nmatch/ntsim,isSignal);
if(nmatch/ntsim>=0.5){
Fill(h,"matchVtxEfficiency",1.,isSignal);
if(ntsim>1){Fill(h,"matchVtxEfficiency2",1.,isSignal);}
if(matchpurity>0.5){Fill(h,"matchVtxEfficiency5",1.,isSignal);}
}else{
Fill(h,"matchVtxEfficiency",0.,isSignal);
if(ntsim>1){Fill(h,"matchVtxEfficiency2",0.,isSignal);}
Fill(h,"matchVtxEfficiency5",0.,isSignal); // no (matchpurity>5) here !!
if(isSignal){
cout << "Signal vertex not matched " << message << " event=" << eventcounter_ << " nmatch=" << nmatch << " ntsim=" << ntsim << endl;
}
}
} // ntsim >0
if(zmatch>-99){
Fill(h,"matchVtxZ",zmatch-ev->z);
Fill(h,"matchVtxZ",zmatch-ev->z,isSignal);
Fill(h,"matchVtxZCum",fabs(zmatch-ev->z));
Fill(h,"matchVtxZCum",fabs(zmatch-ev->z),isSignal);
}else{
Fill(h,"matchVtxZCum",1.0);
Fill(h,"matchVtxZCum",1.0,isSignal);
}
if(fabs(zmatch-ev->z)<zmatch_){
Fill(h,"matchVtxEfficiencyZ",1.,isSignal);
}else{
Fill(h,"matchVtxEfficiencyZ",0.,isSignal);
}
if(ntsim>0) Fill(h, "matchVtxEfficiencyZ1", fabs(zmatch-ev->z)<zmatch_ , isSignal);
if(ntsim>1) Fill(h, "matchVtxEfficiencyZ2", fabs(zmatch-ev->z)<zmatch_ , isSignal);
Fill(h,"vtxMultiplicity",nRecVWithTrk,isSignal);
// efficiency vs number of tracks, use your favorite definition of efficiency here
//if(nmatch>=0.5*ntmatch){ // purity
if(fabs(zmatch-ev->z)<zmatch_){ // zmatch
Fill(h,"vtxFindingEfficiencyVsNtrk",(double) ev->tk.size(),1.);
if(isSignal){
Fill(h,"vtxFindingEfficiencyVsNtrkSignal",ev->tk.size(),1.);
}else{
Fill(h,"vtxFindingEfficiencyVsNtrkPU",(double) ev->tk.size(),1.);
}
}else{
Fill(h,"vtxFindingEfficiencyVsNtrk",(double) ev->tk.size(),0.);
if(isSignal){
Fill(h,"vtxFindingEfficiencyVsNtrkSignal",(double) ev->tk.size(),1.);
}else{
Fill(h,"vtxFindingEfficiencyVsNtrkPU",(double) ev->tk.size(),1.);
}
}
}
Fill(h,"npu1",npu1);
Fill(h,"npu2",npu2);
Fill(h,"nrecvsnpu",npu1,float(nrecvtxs));
Fill(h,"nrecvsnpu2",npu2,float(nrecvtxs));
// --------------------------------------- sim-signal vs rec-tag ---------------------------------------
SimEvent* ev=&(simEvt[0]);
const reco::Vertex* v=&(*recVtxs->begin());
double n=0;
for(vector<TransientTrack>::iterator te=ev->tk.begin(); te!=ev->tk.end(); te++){
const reco::Track& RTe=te->track();
for(trackit_t tv=v->tracks_begin(); tv!=v->tracks_end(); tv++){
const reco::Track & RTv=*(tv->get());
if(RTe.vz()==RTv.vz()) {n++;}
}
}
cout << "Number of tracks in reco tagvtx " << v->tracksSize() << endl;
cout << "Number of selected tracks in sim event vtx " << ev->tk.size() << " (prim=" << ev->tkprim.size() << ")"<<endl;
cout << "Number of tracks in both " << n << endl;
double ntruthmatched=getTruthMatchedVertexTracks(*v).size();
if (ntruthmatched>0){
cout << "TrackPurity = "<< n/ntruthmatched <<endl;
Fill(h,"TagVtxTrkPurity",n/ntruthmatched);
}
if (ev->tk.size()>0){
cout << "TrackEfficiency = "<< n/ev->tk.size() <<endl;
Fill(h,"TagVtxTrkEfficiency",n/ev->tk.size());
}
}
| void PrimaryVertexAnalyzer4PU::beginJob | ( | void | ) | [virtual] |
Reimplemented from edm::EDAnalyzer.
Definition at line 564 of file PrimaryVertexAnalyzer4PU.cc.
References add(), bookVertexHistograms(), gather_cfg::cout, hBS, hDA, estimatePileup::hist, hnoBS, hsimPV, rootFile_, and simUnit_.
{
std::cout << " PrimaryVertexAnalyzer4PU::beginJob conversion from sim units to rec units is " << simUnit_ << std::endl;
rootFile_->cd();
TDirectory *noBS = rootFile_->mkdir("noBS");
noBS->cd();
hnoBS=bookVertexHistograms();
for(std::map<std::string,TH1*>::const_iterator hist=hnoBS.begin(); hist!=hnoBS.end(); hist++){
hist->second->SetDirectory(noBS);
}
TDirectory *withBS = rootFile_->mkdir("BS");
withBS->cd();
hBS=bookVertexHistograms();
for(std::map<std::string,TH1*>::const_iterator hist=hBS.begin(); hist!=hBS.end(); hist++){
hist->second->SetDirectory(withBS);
}
TDirectory *DA = rootFile_->mkdir("DA");
DA->cd();
hDA=bookVertexHistograms();
for(std::map<std::string,TH1*>::const_iterator hist=hDA.begin(); hist!=hDA.end(); hist++){
hist->second->SetDirectory(DA);
}
// TDirectory *PIX = rootFile_->mkdir("PIX");
// PIX->cd();
// hPIX=bookVertexHistograms();
// for(std::map<std::string,TH1*>::const_iterator hist=hPIX.begin(); hist!=hPIX.end(); hist++){
// hist->second->SetDirectory(PIX);
// }
// TDirectory *MVF = rootFile_->mkdir("MVF");
// MVF->cd();
// hMVF=bookVertexHistograms();
// for(std::map<std::string,TH1*>::const_iterator hist=hMVF.begin(); hist!=hMVF.end(); hist++){
// hist->second->SetDirectory(MVF);
// }
rootFile_->cd();
hsimPV["rapidity"] = new TH1F("rapidity","rapidity ",100,-10., 10.);
hsimPV["chRapidity"] = new TH1F("chRapidity","charged rapidity ",100,-10., 10.);
hsimPV["recRapidity"] = new TH1F("recRapidity","reconstructed rapidity ",100,-10., 10.);
hsimPV["pt"] = new TH1F("pt","pt ",100,0., 20.);
hsimPV["xsim"] = new TH1F("xsim","simulated x",100,-0.01,0.01); // 0.01cm = 100 um
hsimPV["ysim"] = new TH1F("ysim","simulated y",100,-0.01,0.01);
hsimPV["zsim"] = new TH1F("zsim","simulated z",100,-20.,20.);
hsimPV["xsim1"] = new TH1F("xsim1","simulated x",100,-4.,4.);
hsimPV["ysim1"] = new TH1F("ysim1","simulated y",100,-4.,4.);
hsimPV["zsim1"] = new TH1F("zsim1","simulated z",100,-40.,40.);
add(hsimPV, new TH1F("xsim2PU","simulated x (Pile-up)",100,-1.,1.));
add(hsimPV, new TH1F("ysim2PU","simulated y (Pile-up)",100,-1.,1.));
add(hsimPV, new TH1F("zsim2PU","simulated z (Pile-up)",100,-20.,20.));
add(hsimPV, new TH1F("xsim2Signal","simulated x (Signal)",100,-1.,1.));
add(hsimPV, new TH1F("ysim2Signal","simulated y (Signal)",100,-1.,1.));
add(hsimPV, new TH1F("zsim2Signal","simulated z (Signal)",100,-20.,20.));
hsimPV["xsim2"] = new TH1F("xsim2","simulated x",100,-1,1); // 0.01cm = 100 um
hsimPV["ysim2"] = new TH1F("ysim2","simulated y",100,-1,1);
hsimPV["zsim2"] = new TH1F("zsim2","simulated z",100,-20.,20.);
hsimPV["xsim3"] = new TH1F("xsim3","simulated x",100,-0.1,0.1); // 0.01cm = 100 um
hsimPV["ysim3"] = new TH1F("ysim3","simulated y",100,-0.1,0.1);
hsimPV["zsim3"] = new TH1F("zsim3","simulated z",100,-20.,20.);
hsimPV["xsimb"] = new TH1F("xsimb","simulated x",100,-0.01,0.01); // 0.01cm = 100 um
hsimPV["ysimb"] = new TH1F("ysimb","simulated y",100,-0.01,0.01);
hsimPV["zsimb"] = new TH1F("zsimb","simulated z",100,-20.,20.);
hsimPV["xsimb1"] = new TH1F("xsimb1","simulated x",100,-0.1,0.1); // 0.01cm = 100 um
hsimPV["ysimb1"] = new TH1F("ysimb1","simulated y",100,-0.1,0.1);
hsimPV["zsimb1"] = new TH1F("zsimb1","simulated z",100,-20.,20.);
add(hsimPV, new TH1F("xbeam","beamspot x",100,-1.,1.));
add(hsimPV, new TH1F("ybeam","beamspot y",100,-1.,1.));
add(hsimPV, new TH1F("zbeam","beamspot z",100,-1.,1.));
add(hsimPV, new TH1F("wxbeam","beamspot sigma x",100,-1.,1.));
add(hsimPV, new TH1F("wybeam","beamspot sigma y",100,-1.,1.));
add(hsimPV, new TH1F("wzbeam","beamspot sigma z",100,-1.,1.));
hsimPV["xsim2"]->StatOverflows(kTRUE);
hsimPV["ysim2"]->StatOverflows(kTRUE);
hsimPV["zsim2"]->StatOverflows(kTRUE);
hsimPV["xsimb"]->StatOverflows(kTRUE);
hsimPV["ysimb"]->StatOverflows(kTRUE);
hsimPV["zsimb"]->StatOverflows(kTRUE);
hsimPV["nsimvtx"] = new TH1F("nsimvtx","# of simulated vertices", 50, -0.5, 49.5);
// hsimPV["nsimtrk"] = new TH1F("nsimtrk","# of simulated tracks", 100, -0.5, 99.5); // not filled right now, also exists in hBS..
// hsimPV["nsimtrk"]->StatOverflows(kTRUE);
hsimPV["nbsimtksinvtx"]= new TH1F("nbsimtksinvtx","simulated tracks in vertex",100,-0.5,99.5);
hsimPV["nbsimtksinvtx"]->StatOverflows(kTRUE);
}
| std::map< std::string, TH1 * > PrimaryVertexAnalyzer4PU::bookVertexHistograms | ( | ) | [private] |
Definition at line 119 of file PrimaryVertexAnalyzer4PU.cc.
References add(), h, M_PI, and lumiQTWidget::t.
Referenced by beginJob().
{
std::map<std::string, TH1*> h;
// release validation histograms used in DoCompare.C
h["nbtksinvtx"] = new TH1F("nbtksinvtx","reconstructed tracks in vertex",40,-0.5,39.5);
h["nbtksinvtxPU"] = new TH1F("nbtksinvtxPU","reconstructed tracks in vertex",40,-0.5,39.5);
h["nbtksinvtxTag"]= new TH1F("nbtksinvtxTag","reconstructed tracks in vertex",40,-0.5,39.5);
h["resx"] = new TH1F("resx","residual x",100,-0.04,0.04);
h["resy"] = new TH1F("resy","residual y",100,-0.04,0.04);
h["resz"] = new TH1F("resz","residual z",100,-0.1,0.1);
h["resz10"] = new TH1F("resz10","residual z",100,-1.0,1.);
h["pullx"] = new TH1F("pullx","pull x",100,-25.,25.);
h["pully"] = new TH1F("pully","pull y",100,-25.,25.);
h["pullz"] = new TH1F("pullz","pull z",100,-25.,25.);
h["vtxchi2"] = new TH1F("vtxchi2","chi squared",100,0.,100.);
h["vtxndf"] = new TH1F("vtxndf","degrees of freedom",500,0.,100.);
h["vtxndfc"] = new TH1F("vtxndfc","expected 2nd highest ndof",500,0.,100.);
h["vtxndfvsntk"] = new TH2F("vtxndfvsntk","ndof vs #tracks",20,0.,100, 20, 0., 200.);
h["vtxndfoverntk"]= new TH1F("vtxndfoverntk","ndof / #tracks",40,0.,2.);
h["vtxndf2overntk"]= new TH1F("vtxndf2overntk","(ndof+2) / #tracks",40,0.,2.);
h["tklinks"] = new TH1F("tklinks","Usable track links",2,-0.5,1.5);
h["nans"] = new TH1F("nans","Illegal values for x,y,z,xx,xy,xz,yy,yz,zz",9,0.5,9.5);
// raw
add(h, new TH1F("szRecVtx","size of recvtx collection",20, -0.5, 19.5));
add(h, new TH1F("isFake","fake vertex",2, -0.5, 1.5));
add(h, new TH1F("isFake1","fake vertex or ndof<0",2, -0.5, 1.5));
add(h, new TH1F("bunchCrossing","bunchCrossing",4000, 0., 4000.));
add(h, new TH2F("bunchCrossingLogNtk","bunchCrossingLogNtk",4000, 0., 4000.,5,0., 5.));
add(h, new TH1F("highpurityTrackFraction","fraction of high purity tracks",20, 0., 1.));
add(h, new TH2F("trkchi2vsndof","vertices chi2 vs ndof",50, 0., 100., 50, 0., 200.));
add(h, new TH1F("trkchi2overndof","vertices chi2 / ndof",50, 0., 5.));
// two track vertices
add(h,new TH2F("2trkchi2vsndof","two-track vertices chi2 vs ndof",40, 0., 10., 20, 0., 20.));
add(h,new TH1F("2trkmassSS","two-track vertices mass (same sign)",100, 0., 2.));
add(h,new TH1F("2trkmassOS","two-track vertices mass (opposite sign)",100, 0., 2.));
add(h,new TH1F("2trkdphi","two-track vertices delta-phi",360, 0, 2*M_PI));
add(h,new TH1F("2trkseta","two-track vertices sum-eta",50, -2., 2.));
add(h,new TH1F("2trkdphicurl","two-track vertices delta-phi (sum eta<0.1)",360, 0, 2*M_PI));
add(h,new TH1F("2trksetacurl","two-track vertices sum-eta (delta-phi<0.1)",50, -2., 2.));
add(h,new TH1F("2trkdetaOS","two-track vertices delta-eta (same sign)",50, -0.5, 0.5));
add(h,new TH1F("2trkdetaSS","two-track vertices delta-eta (opposite sign)",50, -0.5, 0.5));
// two track PU vertices
add(h,new TH1F("2trkmassSSPU","two-track vertices mass (same sign)",100, 0., 2.));
add(h,new TH1F("2trkmassOSPU","two-track vertices mass (opposite sign)",100, 0., 2.));
add(h,new TH1F("2trkdphiPU","two-track vertices delta-phi",360, 0, 2*M_PI));
add(h,new TH1F("2trksetaPU","two-track vertices sum-eta",50, -2., 2.));
add(h,new TH1F("2trkdphicurlPU","two-track vertices delta-phi (sum eta<0.1)",360, 0, 2*M_PI));
add(h,new TH1F("2trksetacurlPU","two-track vertices sum-eta (delta-phi<0.1)",50, -2., 2.));
add(h,new TH1F("2trkdetaOSPU","two-track vertices delta-eta (same sign)",50, -0.5, 0.5));
add(h,new TH1F("2trkdetaSSPU","two-track vertices delta-eta (opposite sign)",50, -0.5, 0.5));
// three track vertices
add(h,new TH2F("2trkchi2vsndof","two-track vertices chi2 vs ndof",40, 0., 10., 20, 0., 20.));
add(h,new TH2F("3trkchi2vsndof","three-track vertices chi2 vs ndof",40, 0., 10., 20, 0., 20.));
add(h,new TH2F("4trkchi2vsndof","four-track vertices chi2 vs ndof",40, 0., 10., 20, 0., 20.));
add(h,new TH2F("5trkchi2vsndof","five-track vertices chi2 vs ndof",40, 0., 10., 20, 0., 20.));
// same for fakes
add(h,new TH2F("fake2trkchi2vsndof","fake two-track vertices chi2 vs ndof",40, 0., 10., 20, 0., 20.));
add(h,new TH2F("fake3trkchi2vsndof","fake three-track vertices chi2 vs ndof",40, 0., 10., 20, 0., 20.));
add(h,new TH2F("fake4trkchi2vsndof","fake four-track vertices chi2 vs ndof",40, 0., 10., 20, 0., 20.));
add(h,new TH2F("fake5trkchi2vsndof","fake five-track vertices chi2 vs ndof",40, 0., 10., 20, 0., 20.));
h["resxr"] = new TH1F("resxr","relative residual x",100,-0.04,0.04);
h["resyr"] = new TH1F("resyr","relative residual y",100,-0.04,0.04);
h["reszr"] = new TH1F("reszr","relative residual z",100,-0.1,0.1);
h["pullxr"] = new TH1F("pullxr","relative pull x",100,-25.,25.);
h["pullyr"] = new TH1F("pullyr","relative pull y",100,-25.,25.);
h["pullzr"] = new TH1F("pullzr","relative pull z",100,-25.,25.);
h["vtxprob"] = new TH1F("vtxprob","chisquared probability",100,0.,1.);
h["eff"] = new TH1F("eff","efficiency",2, -0.5, 1.5);
h["efftag"] = new TH1F("efftag","efficiency tagged vertex",2, -0.5, 1.5);
h["zdistancetag"] = new TH1F("zdistancetag","z-distance between tagged and generated",100, -0.1, 0.1);
h["abszdistancetag"] = new TH1F("abszdistancetag","z-distance between tagged and generated",1000, 0., 1.0);
h["abszdistancetagcum"] = new TH1F("abszdistancetagcum","z-distance between tagged and generated",1000, 0., 1.0);
h["puritytag"] = new TH1F("puritytag","purity of primary vertex tags",2, -0.5, 1.5);
h["effvsptsq"] = new TProfile("effvsptsq","efficiency vs ptsq",20, 0., 10000., 0, 1.);
h["effvsnsimtrk"] = new TProfile("effvsnsimtrk","efficiency vs # simtracks",50, 0., 50., 0, 1.);
h["effvsnrectrk"] = new TProfile("effvsnrectrk","efficiency vs # rectracks",50, 0., 50., 0, 1.);
h["effvsnseltrk"] = new TProfile("effvsnseltrk","efficiency vs # selected tracks",50, 0., 50., 0, 1.);
h["effvsz"] = new TProfile("effvsz","efficiency vs z",20, -20., 20., 0, 1.);
h["effvsz2"] = new TProfile("effvsz2","efficiency vs z (2mm)",20, -20., 20., 0, 1.);
h["effvsr"] = new TProfile("effvsr","efficiency vs r",20, 0., 1., 0, 1.);
h["xresvsntrk"] = new TProfile("xresvsntrk","xresolution vs # vertex tracks",40, 0., 200., 0, 0.01);
h["yresvsntrk"] = new TProfile("yresvsntrk","yresolution vs # vertex tracks",40, 0., 200., 0, 0.01);
h["zresvsntrk"] = new TProfile("zresvsntrk","zresolution vs # vertex tracks",40, 0., 200., 0, 0.01);
h["cpuvsntrk"] = new TProfile("cpuvsntrk","cpu time vs # of fitted tracks",40, 0., 200., 0, 200.);
h["cpucluvsntrk"] = new TProfile("cpucluvsntrk","clustering cpu time # of tracks",40, 0., 200., 0, 10.);
h["cpufit"] = new TH1F("cpufit","cpu time for fitting",100, 0., 200.);
h["cpuclu"] = new TH1F("cpuclu","cpu time for clustering",100, 0., 200.);
h["nbtksinvtx2"] = new TH1F("nbtksinvtx2","reconstructed tracks in vertex",40,0.,200.);
h["nbtksinvtxPU2"] = new TH1F("nbtksinvtxPU2","reconstructed tracks in vertex",40,0.,200.);
h["nbtksinvtxTag2"] = new TH1F("nbtksinvtxTag2","reconstructed tracks in vertex",40,0.,200.);
h["xrec"] = new TH1F("xrec","reconstructed x",100,-0.1,0.1);
h["yrec"] = new TH1F("yrec","reconstructed y",100,-0.1,0.1);
h["zrec"] = new TH1F("zrec","reconstructed z",100,-20.,20.);
h["err1"] = new TH1F("err1","error 1",100,0.,0.1);
h["err2"] = new TH1F("err2","error 2",100,0.,0.1);
h["errx"] = new TH1F("errx","error x",100,0.,0.1);
h["erry"] = new TH1F("erry","error y",100,0.,0.1);
h["errz"] = new TH1F("errz","error z",100,0.,2.0);
h["errz1"] = new TH1F("errz1","error z",100,0.,0.2);
h["zrecNt100"] = new TH1F("zrecNt100","reconstructed z for high multiplicity events",80,-40.,40.);
add(h, new TH2F("zrecvsnt","reconstructed z vs number of tracks",100,-50., 50., 20, 0., 200.));
add(h, new TH2F("xyrec","reconstructed xy",100, -4., 4., 100, -4., 4.));
h["xrecBeam"] = new TH1F("xrecBeam","reconstructed x - beam x",100,-0.1,0.1);
h["yrecBeam"] = new TH1F("yrecBeam","reconstructed y - beam y",100,-0.1,0.1);
h["zrecBeam"] = new TH1F("zrecBeam","reconstructed z - beam z",100,-20.,20.);
h["xrecBeamvsz"] = new TH2F("xrecBeamvsz","reconstructed x - beam x vs z", 20, -20., 20.,100,-0.1,0.1);
h["yrecBeamvsz"] = new TH2F("yrecBeamvsz","reconstructed y - beam y vs z", 20, -20., 20.,100,-0.1,0.1);
h["xrecBeamvszprof"] = new TProfile("xrecBeamvszprof","reconstructed x - beam x vs z-z0", 20, -20., 20.,-0.1,0.1);
h["yrecBeamvszprof"] = new TProfile("yrecBeamvszprof","reconstructed y - beam y vs z-z0", 20, -20., 20.,-0.1,0.1);
add(h, new TH2F("xrecBeamvsdxXBS","reconstructed x - beam x vs resolution",10,0., 0.02, 100, -0.1,0.1)); // just a test
add(h, new TH2F("yrecBeamvsdyXBS","reconstructed z - beam z vs resolution",10,0., 0.02, 100, -0.1,0.1));
add(h, new TH2F("xrecBeamvsdx","reconstructed x - beam x vs resolution",10,0., 0.02, 100, -0.1,0.1));
add(h, new TH2F("yrecBeamvsdy","reconstructed z - beam z vs resolution",10,0., 0.02, 100, -0.1,0.1));
add(h, new TH2F("xrecBeamvsdxR2","reconstructed x - beam x vs resolution",20,0., 0.04, 100, -0.1,0.1));
add(h, new TH2F("yrecBeamvsdyR2","reconstructed z - beam z vs resolution",20,0., 0.04, 100, -0.1,0.1));
// add(h, new TH2F("xrecBeamvsdx","reconstructed x - beam x vs resolution",100,-0.1,0.1, 10, 0., 0.04));
// add(h, new TH2F("yrecBeamvsdy","reconstructed y - beam y vs resolution",100,-0.1,0.1, 10, 0., 0.04));
h["xrecBeamvsdxprof"] = new TProfile("xrecBeamvsdxprof","reconstructed x - beam x vs resolution",10, 0., 0.04,-0.1,0.1 );
h["yrecBeamvsdyprof"] = new TProfile("yrecBeamvsdyprof","reconstructed y - beam y vs resolution",10, 0., 0.04,-0.1,0.1 );
add(h, new TProfile("xrecBeam2vsdx2prof","reconstructed x - beam x vs resolution",10,0., 0.002, 0., 0.01));
add(h, new TProfile("yrecBeam2vsdy2prof","reconstructed y - beam y vs resolution",10,0., 0.002, 0., 0.01));
add(h,new TH2F("xrecBeamvsdx2","reconstructed x - beam x vs resolution",10,0., 0.002, 100, -0.01, 0.01));
add(h,new TH2F("yrecBeamvsdy2","reconstructed y - beam y vs resolution",10,0., 0.002, 100, -0.01, 0.01));
h["xrecb"] = new TH1F("xrecb","reconstructed x - beam x",100,-0.01,0.01);
h["yrecb"] = new TH1F("yrecb","reconstructed y - beam y",100,-0.01,0.01);
h["zrecb"] = new TH1F("zrecb","reconstructed z - beam z",100,-20.,20.);
h["xrec1"] = new TH1F("xrec1","reconstructed x",100,-4,4);
h["yrec1"] = new TH1F("yrec1","reconstructed y",100,-4,4); // should match the sim histos
h["zrec1"] = new TH1F("zrec1","reconstructed z",100,-80.,80.);
h["xrec2"] = new TH1F("xrec2","reconstructed x",100,-1,1);
h["yrec2"] = new TH1F("yrec2","reconstructed y",100,-1,1);
h["zrec2"] = new TH1F("zrec2","reconstructed z",100,-40.,40.);
h["xrec3"] = new TH1F("xrec3","reconstructed x",100,-0.1,0.1);
h["yrec3"] = new TH1F("yrec3","reconstructed y",100,-0.1,0.1);
h["zrec3"] = new TH1F("zrec3","reconstructed z",100,-20.,20.);
add(h, new TH1F("xrecBeamPull","normalized residuals reconstructed x - beam x",100,-20,20));
add(h, new TH1F("yrecBeamPull","normalized residuals reconstructed y - beam y",100,-20,20));
add(h, new TH1F("zdiffrec","z-distance between vertices",200,-20., 20.));
add(h, new TH1F("zdiffrec2","z-distance between ndof>2 vertices",200,-20., 20.));
add(h, new TH1F("zdiffrec4","z-distance between ndof>4 vertices",200,-20., 20.));
add(h, new TH1F("zdiffrec8","z-distance between ndof>8 vertices",200,-20., 20.));
add(h, new TH2F("zvszrec2","z positions of multiple vertices",200,-20., 20., 200,-20., 20.));
add(h, new TH2F("pzvspz2","prob(z) of multiple vertices",100, 0.,1.,100,0., 1.));
add(h, new TH2F("zvszrec4","z positions of multiple vertices",100,-20., 20., 100,-20., 20.));
add(h, new TH2F("pzvspz4","prob(z) of multiple vertices",100, 0.,1.,100,0., 1.));
add(h, new TH2F("zdiffvsz","z-distance vs z",100,-10.,10.,30,-15.,15.));
add(h, new TH2F("zdiffvsz4","z-distance vs z (ndof>4)",100,-10.,10.,30,-15.,15.));
add(h, new TProfile("eff0vsntrec","efficiency vs # reconstructed tracks",50, 0., 50., 0, 1.));
add(h, new TProfile("eff0vsntsel","efficiency vs # selected tracks",50, 0., 50., 0, 1.));
add(h, new TProfile("eff0ndof0vsntsel","efficiency (ndof>0) vs # selected tracks",50, 0., 50., 0, 1.));
add(h, new TProfile("eff0ndof8vsntsel","efficiency (ndof>8) vs # selected tracks",50, 0., 50., 0, 1.));
add(h, new TProfile("eff0ndof2vsntsel","efficiency (ndof>2) vs # selected tracks",50, 0., 50., 0, 1.));
add(h, new TProfile("eff0ndof4vsntsel","efficiency (ndof>4) vs # selected tracks",50, 0., 50., 0, 1.));
add(h, new TH1F("xbeamPUcand","x - beam of PU candidates (ndof>4)",100,-1., 1.));
add(h, new TH1F("ybeamPUcand","y - beam of PU candidates (ndof>4)",100,-1., 1.));
add(h, new TH1F("xbeamPullPUcand","x - beam pull of PU candidates (ndof>4)",100,-20., 20.));
add(h, new TH1F("ybeamPullPUcand","y - beam pull of PU candidates (ndof>4)",100,-20., 20.));
add(h, new TH1F("ndofOverNtk","ndof / ntk of ndidates (ndof>4)",100,0., 2.));
//add(h, new TH1F("sumwOverNtk","sumw / ntk of ndidates (ndof>4)",100,0., 2.));
add(h, new TH1F("ndofOverNtkPUcand","ndof / ntk of ndidates (ndof>4)",100,0., 2.));
//add(h, new TH1F("sumwOverNtkPUcand","sumw / ntk of ndidates (ndof>4)",100,0., 2.));
add(h, new TH1F("zPUcand","z positions of PU candidates (all)",100,-20., 20.));
add(h, new TH1F("zPUcand2","z positions of PU candidates (ndof>2)",100,-20., 20.));
add(h, new TH1F("zPUcand4","z positions of PU candidates (ndof>4)",100,-20., 20.));
add(h, new TH1F("ndofPUcand","ndof of PU candidates (all)",50,0., 100.));
add(h, new TH1F("ndofPUcand2","ndof of PU candidates (ndof>2)",50,0., 100.));
add(h, new TH1F("ndofPUcand4","ndof of PU candidates (ndof>4)",50,0., 100.));
add(h, new TH1F("ndofnr2","second highest ndof",500,0., 100.));
add(h, new TH1F("ndofnr2d1cm","second highest ndof (dz>1cm)",500,0., 100.));
add(h, new TH1F("ndofnr2d2cm","second highest ndof (dz>2cm)",500,0., 100.));
h["nrecvtx"] = new TH1F("nrecvtx","# of reconstructed vertices", 50, -0.5, 49.5);
h["nrecvtx8"] = new TH1F("nrecvtx8","# of reconstructed vertices with ndof>8", 50, -0.5, 49.5);
h["nrecvtx2"] = new TH1F("nrecvtx2","# of reconstructed vertices with ndof>2", 50, -0.5, 49.5);
h["nrecvtx4"] = new TH1F("nrecvtx4","# of reconstructed vertices with ndof>4", 50, -0.5, 49.5);
// h["nsimvtx"] = new TH1F("nsimvtx","# of simulated vertices", 50, -0.5, 49.5);
h["nrectrk"] = new TH1F("nrectrk","# of reconstructed tracks", 100, -0.5, 99.5);
add(h, new TH1F("nsimtrk","# of simulated tracks", 100, -0.5, 99.5));
add(h, new TH1F("nsimtrkSignal","# of simulated tracks (Signal)", 100, -0.5, 99.5));
add(h, new TH1F("nsimtrkPU","# of simulated tracks (PU)", 100, -0.5, 99.5));
h["nsimtrk"]->StatOverflows(kTRUE);
h["nsimtrkPU"]->StatOverflows(kTRUE);
h["nsimtrkSignal"]->StatOverflows(kTRUE);
h["xrectag"] = new TH1F("xrectag","reconstructed x, signal vtx",100,-0.05,0.05);
h["yrectag"] = new TH1F("yrectag","reconstructed y, signal vtx",100,-0.05,0.05);
h["zrectag"] = new TH1F("zrectag","reconstructed z, signal vtx",100,-20.,20.);
h["nrectrk0vtx"] = new TH1F("nrectrk0vtx","# rec tracks no vertex ",100,-0.5, 99.5);
h["nseltrk0vtx"] = new TH1F("nseltrk0vtx","# rec tracks no vertex ",100,-0.5, 99.5);
h["nrecsimtrk"] = new TH1F("nrecsimtrk","# rec tracks matched to sim tracks in vertex",100,-0.5, 99.5);
h["nrecnosimtrk"] = new TH1F("nrecsimtrk","# rec tracks not matched to sim tracks in vertex",100,-0.5, 99.5);
h["trackAssEffvsPt"] = new TProfile("trackAssEffvsPt","track association efficiency vs pt",20, 0., 100., 0, 1.);
// cluster stuff
h["nseltrk"] = new TH1F("nseltrk","# of reconstructed tracks selected for PV", 100, -0.5, 99.5);
h["nclutrkall"] = new TH1F("nclutrkall","# of reconstructed tracks in clusters", 100, -0.5, 99.5);
h["nclutrkvtx"] = new TH1F("nclutrkvtx","# of reconstructed tracks in clusters of reconstructed vertices", 100, -0.5, 99.5);
h["nclu"] = new TH1F("nclu","# of clusters", 100, -0.5, 99.5);
h["nclu0vtx"] = new TH1F("nclu0vtx","# of clusters in events with no PV", 100, -0.5, 99.5);
h["zlost1"] = new TH1F("zlost1","z of lost vertices (bad z)", 100, -20., 20.);
h["zlost2"] = new TH1F("zlost2","z of lost vertices (no matching cluster)", 100, -20., 20.);
h["zlost3"] = new TH1F("zlost3","z of lost vertices (vertex too far from beam)", 100, -20., 20.);
h["zlost4"] = new TH1F("zlost4","z of lost vertices (invalid vertex)", 100, -20., 20.);
h["selstat"] = new TH1F("selstat","selstat", 5, -2.5, 2.5);
// properties of fake vertices (MC only)_
add(h, new TH1F("fakeVtxZNdofgt05","z of fake vertices with ndof>0.5", 100, -20., 20.));
add(h, new TH1F("fakeVtxZNdofgt2","z of fake vertices with ndof>2", 100, -20., 20.));
add(h, new TH1F("fakeVtxZ","z of fake vertices", 100, -20., 20.));
add(h, new TH1F("fakeVtxNdof","ndof of fake vertices", 500,0., 100.));
add(h,new TH1F("fakeVtxNtrk","number of tracks in fake vertex",20,-0.5, 19.5));
add(h,new TH1F("matchedVtxNdof","ndof of matched vertices", 500,0., 100.));
// histograms of track quality (Data and MC)
string types[] = {"all","sel","bachelor","sellost","wgt05","wlt05","M","tagged","untagged","ndof4","PUcand","PUfake"};
for(int t=0; t<12; t++){
add(h, new TH1F(("rapidity_"+types[t]).c_str(),"rapidity ",100,-3., 3.));
h["z0_"+types[t]] = new TH1F(("z0_"+types[t]).c_str(),"z0 ",80,-40., 40.);
h["phi_"+types[t]] = new TH1F(("phi_"+types[t]).c_str(),"phi ",80,-3.14159, 3.14159);
h["eta_"+types[t]] = new TH1F(("eta_"+types[t]).c_str(),"eta ",80,-4., 4.);
h["pt_"+types[t]] = new TH1F(("pt_"+types[t]).c_str(),"pt ",100,0., 20.);
h["found_"+types[t]] = new TH1F(("found_"+types[t]).c_str(),"found hits",20, 0., 20.);
h["lost_"+types[t]] = new TH1F(("lost_"+types[t]).c_str(),"lost hits",20, 0., 20.);
h["nchi2_"+types[t]] = new TH1F(("nchi2_"+types[t]).c_str(),"normalized track chi2",100, 0., 20.);
h["rstart_"+types[t]] = new TH1F(("rstart_"+types[t]).c_str(),"start radius",100, 0., 20.);
h["tfom_"+types[t]] = new TH1F(("tfom_"+types[t]).c_str(),"track figure of merit",100, 0., 100.);
h["expectedInner_"+types[t]] = new TH1F(("expectedInner_"+types[t]).c_str(),"expected inner hits ",10, 0., 10.);
h["expectedOuter_"+types[t]] = new TH1F(("expectedOuter_"+types[t]).c_str(),"expected outer hits ",10, 0., 10.);
h["logtresxy_"+types[t]] = new TH1F(("logtresxy_"+types[t]).c_str(),"log10(track r-phi resolution/um)",100, 0., 5.);
h["logtresz_"+types[t]] = new TH1F(("logtresz_"+types[t]).c_str(),"log10(track z resolution/um)",100, 0., 5.);
h["tpullxy_"+types[t]] = new TH1F(("tpullxy_"+types[t]).c_str(),"track r-phi pull",100, -10., 10.);
add(h, new TH2F( ("lvseta_"+types[t]).c_str(),"cluster length vs eta",60,-3., 3., 20, 0., 20));
add(h, new TH2F( ("lvstanlambda_"+types[t]).c_str(),"cluster length vs tan lambda",60,-6., 6., 20, 0., 20));
add(h, new TH1F( ("restrkz_"+types[t]).c_str(),"z-residuals (track vs vertex)", 200, -5., 5.));
add(h, new TH2F( ("restrkzvsphi_"+types[t]).c_str(),"z-residuals (track - vertex)", 12,-3.14159,3.14159,100, -5., 5.));
add(h, new TH2F( ("restrkzvseta_"+types[t]).c_str(),"z-residuals (track - vertex)", 12,-3.,3.,200, -5., 5.));
add(h, new TH2F( ("pulltrkzvsphi_"+types[t]).c_str(),"normalized z-residuals (track - vertex)", 12,-3.14159,3.14159,100, -5., 5.));
add(h, new TH2F( ("pulltrkzvseta_"+types[t]).c_str(),"normalized z-residuals (track - vertex)", 12,-3.,3.,100, -5., 5.));
add(h, new TH1F( ("pulltrkz_"+types[t]).c_str(),"normalized z-residuals (track vs vertex)", 100, -5., 5.));
add(h, new TH1F( ("sigmatrkz0_"+types[t]).c_str(),"z-resolution (excluding beam)", 100, 0., 5.));
add(h, new TH1F( ("sigmatrkz_"+types[t]).c_str(),"z-resolution (including beam)", 100,0., 5.));
add(h, new TH1F( ("nbarrelhits_"+types[t]).c_str(),"number of pixel barrel hits", 10, 0., 10.));
add(h, new TH1F( ("nbarrelLayers_"+types[t]).c_str(),"number of pixel barrel layers", 10, 0., 10.));
add(h, new TH1F( ("nPxLayers_"+types[t]).c_str(),"number of pixel layers (barrel+endcap)", 10, 0., 10.));
add(h, new TH1F( ("nSiLayers_"+types[t]).c_str(),"number of Tracker layers ", 20, 0., 20.));
add(h, new TH1F( ("trackAlgo_"+types[t]).c_str(),"track algorithm ", 30, 0., 30.));
add(h, new TH1F( ("trackQuality_"+types[t]).c_str(),"track quality ", 7, -1., 6.));
}
add(h, new TH1F("trackWt","track weight in vertex",100,0., 1.));
h["nrectrk"]->StatOverflows(kTRUE);
h["nrectrk"]->StatOverflows(kTRUE);
h["nrectrk0vtx"]->StatOverflows(kTRUE);
h["nseltrk0vtx"]->StatOverflows(kTRUE);
h["nrecsimtrk"]->StatOverflows(kTRUE);
h["nrecnosimtrk"]->StatOverflows(kTRUE);
h["nseltrk"]->StatOverflows(kTRUE);
h["nbtksinvtx"]->StatOverflows(kTRUE);
h["nbtksinvtxPU"]->StatOverflows(kTRUE);
h["nbtksinvtxTag"]->StatOverflows(kTRUE);
h["nbtksinvtx2"]->StatOverflows(kTRUE);
h["nbtksinvtxPU2"]->StatOverflows(kTRUE);
h["nbtksinvtxTag2"]->StatOverflows(kTRUE);
// pile-up and track assignment related histograms (MC)
h["npu0"] =new TH1F("npu0","Number of simulated vertices",40,0.,40.);
h["npu1"] =new TH1F("npu1","Number of simulated vertices with >0 track",40,0.,40.);
h["npu2"] =new TH1F("npu2","Number of simulated vertices with >1 track",40,0.,40.);
h["nrecv"] =new TH1F("nrecv","# of reconstructed vertices", 40, 0, 40);
add(h,new TH2F("nrecvsnpu","#rec vertices vs number of sim vertices with >0 tracks", 40, 0., 40, 40, 0, 40));
add(h,new TH2F("nrecvsnpu2","#rec vertices vs number of sim vertices with >1 tracks", 40, 0., 40, 40, 0, 40));
add(h,new TH1F("sumpt","sumpt of simulated tracks",100,0.,100.));
add(h,new TH1F("sumptSignal","sumpt of simulated tracks in Signal events",100,0.,200.));
add(h,new TH1F("sumptPU","sumpt of simulated tracks in PU events",100,0.,200.));
add(h,new TH1F("sumpt2rec","sumpt2 of reconstructed and matched tracks",100,0.,100.));
add(h,new TH1F("sumpt2","sumpt2 of simulated tracks",100,0.,100.));
add(h,new TH1F("sumpt2Signal","sumpt2 of simulated tracks in Signal events",100,0.,200.));
add(h,new TH1F("sumpt2PU","sumpt2 of simulated tracks in PU events",100,0.,200.));
add(h,new TH1F("sumpt2rec","sumpt2 of reconstructed and matched tracks",100,0.,100.));
add(h,new TH1F("sumpt2recSignal","sumpt2 of reconstructed and matched tracks in Signal events",100,0.,200.));
add(h,new TH1F("sumpt2recPU","sumpt2 of reconstructed and matched tracks in PU events",100,0.,200.));
add(h,new TH1F("nRecTrkInSimVtx","number of reco tracks matched to sim-vertex", 101, 0., 100));
add(h,new TH1F("nRecTrkInSimVtxSignal","number of reco tracks matched to signal sim-vertex", 101, 0., 100));
add(h,new TH1F("nRecTrkInSimVtxPU","number of reco tracks matched to PU-vertex", 101, 0., 100));
add(h,new TH1F("nPrimRecTrkInSimVtx","number of reco primary tracks matched to sim-vertex", 101, 0., 100));
add(h,new TH1F("nPrimRecTrkInSimVtxSignal","number of reco primary tracks matched to signal sim-vertex", 101, 0., 100));
add(h,new TH1F("nPrimRecTrkInSimVtxPU","number of reco primary tracks matched to PU-vertex", 101, 0., 100));
// obsolete, scheduled for removal
add(h,new TH1F("recPurity","track purity of reconstructed vertices", 101, 0., 1.01));
add(h,new TH1F("recPuritySignal","track purity of reconstructed Signal vertices", 101, 0., 1.01));
add(h,new TH1F("recPurityPU","track purity of reconstructed PU vertices", 101, 0., 1.01));
// end of obsolete
add(h,new TH1F("recmatchPurity","track purity of all vertices", 101, 0., 1.01));
add(h,new TH1F("recmatchPurityTag","track purity of tagged vertices", 101, 0., 1.01));
add(h,new TH1F("recmatchPurityTagSignal","track purity of tagged Signal vertices", 101, 0., 1.01));
add(h,new TH1F("recmatchPurityTagPU","track purity of tagged PU vertices", 101, 0., 1.01));
add(h,new TH1F("recmatchPuritynoTag","track purity of untagged vertices", 101, 0., 1.01));
add(h,new TH1F("recmatchPuritynoTagSignal","track purity of untagged Signal vertices", 101, 0., 1.01));
add(h,new TH1F("recmatchPuritynoTagPU","track purity of untagged PU vertices", 101, 0., 1.01));
add(h,new TH1F("recmatchvtxs","number of sim vertices contributing to recvtx", 10, 0., 10.));
add(h,new TH1F("recmatchvtxsTag","number of sim vertices contributing to recvtx", 10, 0., 10.));
add(h,new TH1F("recmatchvtxsnoTag","number of sim vertices contributing to recvtx", 10, 0., 10.));
//
add(h,new TH1F("trkAssignmentEfficiency", "track to vertex assignment efficiency", 101, 0., 1.01) );
add(h,new TH1F("trkAssignmentEfficiencySignal", "track to signal vertex assignment efficiency", 101, 0., 1.01) );
add(h,new TH1F("trkAssignmentEfficiencyPU", "track to PU vertex assignment efficiency", 101, 0., 1.01) );
add(h,new TH1F("primtrkAssignmentEfficiency", "track to vertex assignment efficiency", 101, 0., 1.01) );
add(h,new TH1F("primtrkAssignmentEfficiencySignal", "track to signal vertex assignment efficiency", 101, 0., 1.01) );
add(h,new TH1F("primtrkAssignmentEfficiencyPU", "track to PU vertex assignment efficiency", 101, 0., 1.01) );
add(h,new TH1F("vtxMultiplicity", "number of rec vertices containg tracks from one true vertex", 10, 0., 10.) );
add(h,new TH1F("vtxMultiplicitySignal", "number of rec vertices containg tracks from the Signal Vertex", 10, 0., 10.) );
add(h,new TH1F("vtxMultiplicityPU", "number of rec vertices containg tracks from a PU Vertex", 10, 0., 10.) );
add(h,new TProfile("vtxFindingEfficiencyVsNtrk","finding efficiency vs number of associated rec tracks",100, 0., 100., 0., 1.) );
add(h,new TProfile("vtxFindingEfficiencyVsNtrkSignal","Signal vertex finding efficiency vs number of associated rec tracks",100, 0., 100., 0., 1.) );
add(h,new TProfile("vtxFindingEfficiencyVsNtrkPU","PU vertex finding efficiency vs number of associated rec tracks",100, 0., 100., 0., 1.) );
add(h,new TH1F("TagVtxTrkPurity","TagVtxTrkPurity",100,0.,1.01));
add(h,new TH1F("TagVtxTrkEfficiency","TagVtxTrkEfficiency",100,0.,1.01));
add(h,new TH1F("matchVtxFraction","fraction of sim vertex tracks found in a recvertex",101,0,1.01));
add(h,new TH1F("matchVtxFractionSignal","fraction of sim vertex tracks found in a recvertex",101,0,1.01));
add(h,new TH1F("matchVtxFractionPU","fraction of sim vertex tracks found in a recvertex",101,0,1.01));
add(h,new TH1F("matchVtxFractionCum","fraction of sim vertex tracks found in a recvertex",101,0,1.01));
add(h,new TH1F("matchVtxFractionCumSignal","fraction of sim vertexs track found in a recvertex",101,0,1.01));
add(h,new TH1F("matchVtxFractionCumPU","fraction of sim vertex tracks found in a recvertex",101,0,1.01));
add(h,new TH1F("matchVtxEfficiency","efficiency for finding matching rec vertex",2,-0.5,1.5));
add(h,new TH1F("matchVtxEfficiencySignal","efficiency for finding matching rec vertex",2,-0.5,1.5));
add(h,new TH1F("matchVtxEfficiencyPU","efficiency for finding matching rec vertex",2,-0.5,1.5));
add(h,new TH1F("matchVtxEfficiency2","efficiency for finding matching rec vertex (nt>1)",2,-0.5,1.5));
add(h,new TH1F("matchVtxEfficiency2Signal","efficiency for finding matching rec vertex (nt>1)",2,-0.5,1.5));
add(h,new TH1F("matchVtxEfficiency2PU","efficiency for finding matching rec vertex (nt>1)",2,-0.5,1.5));
add(h,new TH1F("matchVtxEfficiency5","efficiency for finding matching rec vertex (purity>0.5)",2,-0.5,1.5));
add(h,new TH1F("matchVtxEfficiency5Signal","efficiency for finding matching rec vertex (purity>0.5)",2,-0.5,1.5));
add(h,new TH1F("matchVtxEfficiency5PU","efficiency for finding matching rec vertex (purity>0.5)",2,-0.5,1.5));
add(h,new TH1F("matchVtxEfficiencyZ","efficiency for finding matching rec vertex within 1 mm",2,-0.5,1.5));
add(h,new TH1F("matchVtxEfficiencyZSignal","efficiency for finding matching rec vertex within 1 mm",2,-0.5,1.5));
add(h,new TH1F("matchVtxEfficiencyZPU","efficiency for finding matching rec vertex within 1 mm",2,-0.5,1.5));
add(h,new TH1F("matchVtxEfficiencyZ1","efficiency for finding matching rec vertex within 1 mm (nt>0)",2,-0.5,1.5));
add(h,new TH1F("matchVtxEfficiencyZ1Signal","efficiency for finding matching rec vertex within 1 mm (nt>0)",2,-0.5,1.5));
add(h,new TH1F("matchVtxEfficiencyZ1PU","efficiency for finding matching rec vertex within 1 mm (nt>0)",2,-0.5,1.5));
add(h,new TH1F("matchVtxEfficiencyZ2","efficiency for finding matching rec vertex within 1 mm (nt>1)",2,-0.5,1.5));
add(h,new TH1F("matchVtxEfficiencyZ2Signal","efficiency for finding matching rec vertex within 1 mm (nt>1)",2,-0.5,1.5));
add(h,new TH1F("matchVtxEfficiencyZ2PU","efficiency for finding matching rec vertex within 1 mm (nt>1)",2,-0.5,1.5));
add(h,new TH1F("matchVtxZ","z distance to matched recvtx",100, -0.1, 0.1));
add(h,new TH1F("matchVtxZPU","z distance to matched recvtx",100, -0.1, 0.1));
add(h,new TH1F("matchVtxZSignal","z distance to matched recvtx",100, -0.1, 0.1));
add(h,new TH1F("matchVtxZCum","z distance to matched recvtx",1001,0, 1.01));
add(h,new TH1F("matchVtxZCumSignal","z distance to matched recvtx",1001,0, 1.01));
add(h,new TH1F("matchVtxZCumPU","z distance to matched recvtx",1001,0, 1.01));
add(h, new TH1F("unmatchedVtx","number of unmatched rec vertices (fakes)",10,0.,10.));
add(h, new TH1F("unmatchedVtxFrac","fraction of unmatched rec vertices (fakes)",1000,0.,1.0));
add(h, new TH1F("unmatchedVtxZ","z of unmached rec vertices (fakes)",100,-20., 20.));
add(h, new TH1F("unmatchedVtxNdof","ndof of unmatched rec vertices (fakes)", 500,0., 100.));
add(h,new TH2F("correctlyassigned","pt and eta of correctly assigned tracks", 60, -3., 3., 100, 0, 10.));
add(h,new TH2F("misassigned","pt and eta of mis assigned tracks", 60, -3., 3., 100, 0, 10.));
add(h,new TH1F("ptcat","pt of correctly assigned tracks", 100, 0, 10.));
add(h,new TH1F("etacat","eta of correctly assigned tracks", 60, -3., 3.));
add(h,new TH1F("phicat","phi of correctly assigned tracks", 100, -3.14159, 3.14159));
add(h,new TH1F("dzcat","dz of correctly assigned tracks", 100, 0., 1.));
add(h,new TH1F("ptmis","pt of mis-assigned tracks", 100, 0, 10.));
add(h,new TH1F("etamis","eta of mis-assigned tracks", 60, -3., 3.));
add(h,new TH1F("phimis","phi of mis-assigned tracks",100, -3.14159, 3.14159));
add(h,new TH1F("dzmis","dz of mis-assigned tracks", 100, 0., 1.));
add(h,new TH1F("Tc","Tc computed with Truth matched Reco Tracks",100,0.,20.));
add(h,new TH1F("TcSignal","Tc of signal vertices computed with Truth matched Reco Tracks",100,0.,20.));
add(h,new TH1F("TcPU","Tc of PU vertices computed with Truth matched Reco Tracks",100,0.,20.));
add(h,new TH1F("logTc","log Tc computed with Truth matched Reco Tracks",100,-2.,8.));
add(h,new TH1F("logTcSignal","log Tc of signal vertices computed with Truth matched Reco Tracks",100,-2.,8.));
add(h,new TH1F("logTcPU","log Tc of PU vertices computed with Truth matched Reco Tracks",100,-2.,8.));
add(h,new TH1F("xTc","Tc of merged clusters",100,0.,20.));
add(h,new TH1F("xTcSignal","Tc of signal vertices merged with PU",100,0.,20.));
add(h,new TH1F("xTcPU","Tc of merged PU vertices",100,0.,20.));
add(h,new TH1F("logxTc","log Tc merged vertices",100,-2.,8.));
add(h,new TH1F("logxTcSignal","log Tc of signal vertices merged with PU",100,-2.,8.));
add(h,new TH1F("logxTcPU","log Tc of merged PU vertices ",100,-2.,8.));
add(h,new TH1F("logChisq","Chisq/ntrk computed with Truth matched Reco Tracks",100,-2.,8.));
add(h,new TH1F("logChisqSignal","Chisq/ntrk of signal vertices computed with Truth matched Reco Tracks",100,-2.,8.));
add(h,new TH1F("logChisqPU","Chisq/ntrk of PU vertices computed with Truth matched Reco Tracks",100,-2.,8.));
add(h,new TH1F("logxChisq","Chisq/ntrk of merged clusters",100,-2.,8.));
add(h,new TH1F("logxChisqSignal","Chisq/ntrk of signal vertices merged with PU",100,-2.,8.));
add(h,new TH1F("logxChisqPU","Chisq/ntrk of merged PU vertices",100,-2.,8.));
add(h,new TH1F("Chisq","Chisq/ntrk computed with Truth matched Reco Tracks",100,0.,20.));
add(h,new TH1F("ChisqSignal","Chisq/ntrk of signal vertices computed with Truth matched Reco Tracks",100,0.,20.));
add(h,new TH1F("ChisqPU","Chisq/ntrk of PU vertices computed with Truth matched Reco Tracks",100,0.,20.));
add(h,new TH1F("xChisq","Chisq/ntrk of merged clusters",100,0.,20.));
add(h,new TH1F("xChisqSignal","Chisq/ntrk of signal vertices merged with PU",100,0.,20.));
add(h,new TH1F("xChisqPU","Chisq/ntrk of merged PU vertices",100,0.,20.));
add(h,new TH1F("dzmax","dzmax computed with Truth matched Reco Tracks",100,0.,2.));
add(h,new TH1F("dzmaxSignal","dzmax of signal vertices computed with Truth matched Reco Tracks",100,0.,2.));
add(h,new TH1F("dzmaxPU","dzmax of PU vertices computed with Truth matched Reco Tracks",100,0.,2.));
add(h,new TH1F("xdzmax","dzmax of merged clusters",100,0.,2.));
add(h,new TH1F("xdzmaxSignal","dzmax of signal vertices merged with PU",100,0.,2.));
add(h,new TH1F("xdzmaxPU","dzmax of merged PU vertices",100,0.,2.));
add(h,new TH1F("dztrim","dzmax computed with Truth matched Reco Tracks",100,0.,2.));
add(h,new TH1F("dztrimSignal","dzmax of signal vertices computed with Truth matched Reco Tracks",100,0.,2.));
add(h,new TH1F("dztrimPU","dzmax of PU vertices computed with Truth matched Reco Tracks",100,0.,2.));
add(h,new TH1F("xdztrim","dzmax of merged clusters",100,0.,2.));
add(h,new TH1F("xdztrimSignal","dzmax of signal vertices merged with PU",100,0.,2.));
add(h,new TH1F("xdztrimPU","dzmax of merged PU vertices",100,0.,2.));
add(h,new TH1F("m4m2","m4m2 computed with Truth matched Reco Tracks",100,0.,100.));
add(h,new TH1F("m4m2Signal","m4m2 of signal vertices computed with Truth matched Reco Tracks",100,0.,100.));
add(h,new TH1F("m4m2PU","m4m2 of PU vertices computed with Truth matched Reco Tracks",100,0.,100.));
add(h,new TH1F("xm4m2","m4m2 of merged clusters",100,0.,100.));
add(h,new TH1F("xm4m2Signal","m4m2 of signal vertices merged with PU",100,0.,100.));
add(h,new TH1F("xm4m2PU","m4m2 of merged PU vertices",100,0.,100.));
return h;
}
| void PrimaryVertexAnalyzer4PU::Cumulate | ( | TH1 * | h | ) | [inline, private] |
Definition at line 234 of file PrimaryVertexAnalyzer4PU.h.
References gather_cfg::cout, and funct::integral().
Referenced by endJob().
{
if((h->GetEntries()==0) || (h->Integral()<=0) ){
std::cout << "DEBUG : Cumulate called with empty histogram " << h->GetTitle() << std::endl;
return;
}
std::cout << "DEBUG : cumulating " << h->GetTitle() << std::endl;
try{
h->ComputeIntegral();
Double_t * integral=h->GetIntegral();
h->SetContent(integral);
}catch(...){
std::cout << "DEBUG : an error occurred cumulating " << h->GetTitle() << std::endl;
}
std::cout << "DEBUG : cumulating " << h->GetTitle() << "done " << std::endl;
}
| void PrimaryVertexAnalyzer4PU::dumpHitInfo | ( | const reco::Track & | t | ) | [private] |
Definition at line 1102 of file PrimaryVertexAnalyzer4PU.cc.
References Reference_intrackfit_cff::barrel, gather_cfg::cout, reco::TrackBase::eta(), edm::Ref< C, T, F >::isNonnull(), GeomDetEnumerators::PixelBarrel, reco::TrackBase::pt(), reco::Track::recHitsBegin(), reco::Track::recHitsEnd(), and align::Tracker.
{
// collect some info on hits and clusters
int longesthit=0, nbarrel=0;
cout << Form("%5.2f %5.2f : ",t.pt(),t.eta());
for(trackingRecHit_iterator hit=t.recHitsBegin(); hit!=t.recHitsEnd(); hit++){
if ((**hit).isValid() && (**hit).geographicalId().det() == DetId::Tracker ){
bool barrel = DetId((**hit).geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelBarrel);
//bool endcap = DetId::DetId((**hit).geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelEndcap);
if (barrel){
nbarrel++;
const SiPixelRecHit *pixhit = dynamic_cast<const SiPixelRecHit*>( &(**hit));
edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const& clust = (*pixhit).cluster();
if (clust.isNonnull()) {
cout << Form("%4d",clust->sizeY());
if (clust->sizeY()-longesthit>0) longesthit=clust->sizeY();
}
}
}
}
//cout << endl;
}
| void PrimaryVertexAnalyzer4PU::endJob | ( | void | ) | [virtual] |
Reimplemented from edm::EDAnalyzer.
Definition at line 659 of file PrimaryVertexAnalyzer4PU.cc.
References gather_cfg::cout, Cumulate(), hBS, hDA, estimatePileup::hist, hnoBS, hsimPV, i, AlCaHLTBitMon_ParallelJobs::p, and rootFile_.
{
std::cout << "this is void PrimaryVertexAnalyzer4PU::endJob() " << std::endl;
//cumulate some histos
double sumDA=0,sumBS=0,sumnoBS=0;
// double sumPIX=0,sumMVF=0;
for(int i=101; i>0; i--){
sumDA+=hDA["matchVtxFractionSignal"]->GetBinContent(i)/hDA["matchVtxFractionSignal"]->Integral();
hDA["matchVtxFractionCumSignal"]->SetBinContent(i,sumDA);
sumBS+=hBS["matchVtxFractionSignal"]->GetBinContent(i)/hBS["matchVtxFractionSignal"]->Integral();
hBS["matchVtxFractionCumSignal"]->SetBinContent(i,sumBS);
sumnoBS+=hnoBS["matchVtxFractionSignal"]->GetBinContent(i)/hnoBS["matchVtxFractionSignal"]->Integral();
hnoBS["matchVtxFractionCumSignal"]->SetBinContent(i,sumnoBS);
// sumPIX+=hPIX["matchVtxFractionSignal"]->GetBinContent(i)/hPIX["matchVtxFractionSignal"]->Integral();
// hPIX["matchVtxFractionCumSignal"]->SetBinContent(i,sumPIX);
// sumMVF+=hMVF["matchVtxFractionSignal"]->GetBinContent(i)/hMVF["matchVtxFractionSignal"]->Integral();
// hMVF["matchVtxFractionCumSignal"]->SetBinContent(i,sumMVF);
}
sumDA=0,sumBS=0,sumnoBS=0;
//sumPIX=0,sumMVF=0;
for(int i=1; i<1001; i++){
sumDA+=hDA["abszdistancetag"]->GetBinContent(i);
hDA["abszdistancetagcum"]->SetBinContent(i,sumDA/float(hDA["abszdistancetag"]->GetEntries()));
sumBS+=hBS["abszdistancetag"]->GetBinContent(i);
hBS["abszdistancetagcum"]->SetBinContent(i,sumBS/float(hBS["abszdistancetag"]->GetEntries()));
sumnoBS+=hnoBS["abszdistancetag"]->GetBinContent(i);
hnoBS["abszdistancetagcum"]->SetBinContent(i,sumnoBS/float(hnoBS["abszdistancetag"]->GetEntries()));
// sumPIX+=hPIX["abszdistancetag"]->GetBinContent(i);
// hPIX["abszdistancetagcum"]->SetBinContent(i,sumPIX/float(hPIX["abszdistancetag"]->GetEntries()));
// sumMVF+=hMVF["abszdistancetag"]->GetBinContent(i);
// hMVF["abszdistancetagcum"]->SetBinContent(i,sumMVF/float(hMVF["abszdistancetag"]->GetEntries()));
}
Cumulate(hBS["matchVtxZCum"]); Cumulate(hBS["matchVtxZCumSignal"]); Cumulate(hBS["matchVtxZCumPU"]);
Cumulate(hnoBS["matchVtxZCum"]); Cumulate(hnoBS["matchVtxZCumSignal"]); Cumulate(hnoBS["matchVtxZCumPU"]);
Cumulate(hDA["matchVtxZCum"]); Cumulate(hDA["matchVtxZCumSignal"]); Cumulate(hDA["matchVtxZCumPU"]);
/*
h->ComputeIntegral();
Double_t *integral = h->GetIntegral();
h->SetContent(integral);
*/
// make a reference for ndofnr2
//hDA["vtxndof"]->ComputeIntegral();
//Double_t *integral = hDA["vtxndf"]->GetIntegral();
// h->SetContent(integral);
double p;
for(int i=1; i<501; i++){
if(hDA["vtxndf"]->GetEntries()>0){
p= hDA["vtxndf"]->Integral(i,501)/hDA["vtxndf"]->GetEntries(); hDA["vtxndfc"]->SetBinContent(i,p*hDA["vtxndf"]->GetBinContent(i));
}
if(hBS["vtxndf"]->GetEntries()>0){
p= hBS["vtxndf"]->Integral(i,501)/hBS["vtxndf"]->GetEntries(); hBS["vtxndfc"]->SetBinContent(i,p*hBS["vtxndf"]->GetBinContent(i));
}
if(hnoBS["vtxndf"]->GetEntries()>0){
p=hnoBS["vtxndf"]->Integral(i,501)/hnoBS["vtxndf"]->GetEntries(); hnoBS["vtxndfc"]->SetBinContent(i,p*hnoBS["vtxndf"]->GetBinContent(i));
}
// if(hPIX["vtxndf"]->GetEntries()>0){
// p=hPIX["vtxndf"]->Integral(i,501)/hPIX["vtxndf"]->GetEntries(); hPIX["vtxndfc"]->SetBinContent(i,p*hPIX["vtxndf"]->GetBinContent(i));
// }
}
rootFile_->cd();
for(std::map<std::string,TH1*>::const_iterator hist=hsimPV.begin(); hist!=hsimPV.end(); hist++){
std::cout << "writing " << hist->first << std::endl;
hist->second->Write();
}
rootFile_->Write();
std::cout << "PrimaryVertexAnalyzer4PU::endJob: done" << std::endl;
}
| void PrimaryVertexAnalyzer4PU::Fill | ( | std::map< std::string, TH1 * > & | h, |
| std::string | s, | ||
| double | x | ||
| ) | [inline, private] |
Definition at line 186 of file PrimaryVertexAnalyzer4PU.h.
References gather_cfg::cout, and alignCSCRings::s.
Referenced by analyze(), analyzeVertexCollection(), analyzeVertexCollectionTP(), Fill(), fillTrackHistos(), and printEventSummary().
| void PrimaryVertexAnalyzer4PU::Fill | ( | std::map< std::string, TH1 * > & | h, |
| std::string | s, | ||
| double | x, | ||
| double | y | ||
| ) | [inline, private] |
Definition at line 195 of file PrimaryVertexAnalyzer4PU.h.
References gather_cfg::cout, and alignCSCRings::s.
| void PrimaryVertexAnalyzer4PU::Fill | ( | std::map< std::string, TH1 * > & | h, |
| std::string | s, | ||
| double | x, | ||
| bool | signal | ||
| ) | [inline, private] |
Definition at line 204 of file PrimaryVertexAnalyzer4PU.h.
References gather_cfg::cout, and alignCSCRings::s.
{
if(h.count(s)==0){
std::cout << "Trying to fill non-exiting Histogram named " << s << std::endl;
return;
}
h[s]->Fill(x);
if(signal){
if(h.count(s+"Signal")==0){
std::cout << "Trying to fill non-exiting Histogram named " << s+"Signal" << std::endl;
return;
}
h[s+"Signal"]->Fill(x);
}else{
if(h.count(s+"PU")==0){
std::cout << "Trying to fill non-exiting Histogram named " << s+"PU" << std::endl;
return;
}
h[s+"PU"]->Fill(x);
}
}
| void PrimaryVertexAnalyzer4PU::Fill | ( | std::map< std::string, TH1 * > & | h, |
| std::string | s, | ||
| bool | yesno, | ||
| bool | signal | ||
| ) | [inline, private] |
| void PrimaryVertexAnalyzer4PU::fillTrackHistos | ( | std::map< std::string, TH1 * > & | h, |
| const std::string & | ttype, | ||
| const reco::Track & | t, | ||
| const reco::Vertex * | v = NULL |
||
| ) | [private] |
Definition at line 1003 of file PrimaryVertexAnalyzer4PU.cc.
References reco::TrackBase::algo(), Reference_intrackfit_cff::barrel, funct::cos(), reco::TrackBase::d0Error(), reco::TrackBase::dxy(), reco::TrackBase::dzError(), reco::TrackBase::eta(), fBfield_, Fill(), reco::Track::found(), reco::TrackBase::hitPattern(), reco::Track::innerPosition(), edm::Ref< C, T, F >::isNonnull(), reco::TrackBase::lambda(), create_public_lumi_plots::log, reco::Track::lost(), reco::Vertex::ndof(), reco::TrackBase::normalizedChi2(), NULL, reco::HitPattern::numberOfHits(), reco::TrackBase::phi(), GeomDetEnumerators::PixelBarrel, reco::HitPattern::pixelBarrelLayersWithMeasurement(), reco::HitPattern::pixelLayersWithMeasurement(), reco::Vertex::position(), position, reco::BeamSpot::position(), funct::pow(), reco::TrackBase::pt(), lumiQueryAPI::q, reco::TrackBase::qoverp(), reco::TrackBase::qualityMask(), reco::Track::recHitsBegin(), reco::Track::recHitsEnd(), reco::TransientTrack::setBeamSpot(), funct::sin(), mathSSE::sqrt(), reco::TransientTrack::stateAtBeamLine(), funct::tan(), theB_, reco::TrackBase::theta(), reco::TransientTrack::track(), align::Tracker, reco::TrackBase::trackerExpectedHitsInner(), reco::TrackBase::trackerExpectedHitsOuter(), reco::HitPattern::trackerLayersWithMeasurement(), TrajectoryStateClosestToBeamLine::trackStateAtPCA(), groupFilesInBlocks::tt, vertexBeamSpot_, reco::TrackBase::vx(), reco::TrackBase::vy(), reco::TrackBase::vz(), wxy2_, reco::BeamSpot::x0(), reco::BeamSpot::y0(), and z.
Referenced by analyzeVertexCollection().
{
Fill(h,"rapidity_"+ttype,t.eta());
Fill(h,"z0_"+ttype,t.vz());
Fill(h,"phi_"+ttype,t.phi());
Fill(h,"eta_"+ttype,t.eta());
Fill(h,"pt_"+ttype,t.pt());
Fill(h,"found_"+ttype,t.found());
Fill(h,"lost_"+ttype,t.lost());
Fill(h,"nchi2_"+ttype,t.normalizedChi2());
Fill(h,"rstart_"+ttype,(t.innerPosition()).Rho());
double d0Error=t.d0Error();
double d0=t.dxy(vertexBeamSpot_.position());
if (d0Error>0){
Fill(h,"logtresxy_"+ttype,log(d0Error/0.0001)/log(10.));
Fill(h,"tpullxy_"+ttype,d0/d0Error);
}
//double z0=t.vz();
double dzError=t.dzError();
if(dzError>0){
Fill(h,"logtresz_"+ttype,log(dzError/0.0001)/log(10.));
}
//
Fill(h,"sigmatrkz_"+ttype,sqrt(pow(t.dzError(),2)+wxy2_/pow(tan(t.theta()),2)));
Fill(h,"sigmatrkz0_"+ttype,t.dzError());
// track vs vertex
if((! (v==NULL)) && (v->ndof()>10.)) {
// emulate clusterizer input
//const TransientTrack & tt = theB_->build(&t); wrong !!!!
TransientTrack tt = theB_->build(&t); tt.setBeamSpot(vertexBeamSpot_); // need the setBeamSpot !
double z=(tt.stateAtBeamLine().trackStateAtPCA()).position().z();
double tantheta=tan((tt.stateAtBeamLine().trackStateAtPCA()).momentum().theta());
double dz2= pow(tt.track().dzError(),2)+wxy2_/pow(tantheta,2);
Fill(h,"restrkz_"+ttype,z-v->position().z());
Fill(h,"restrkzvsphi_"+ttype,t.phi(), z-v->position().z());
Fill(h,"restrkzvseta_"+ttype,t.eta(), z-v->position().z());
Fill(h,"pulltrkzvsphi_"+ttype,t.phi(), (z-v->position().z())/sqrt(dz2));
Fill(h,"pulltrkzvseta_"+ttype,t.eta(), (z-v->position().z())/sqrt(dz2));
Fill(h,"pulltrkz_"+ttype,(z-v->position().z())/sqrt(dz2));
double x1=t.vx()-vertexBeamSpot_.x0(); double y1=t.vy()-vertexBeamSpot_.y0();
double kappa=-0.002998*fBfield_*t.qoverp()/cos(t.theta());
double D0=x1*sin(t.phi())-y1*cos(t.phi())-0.5*kappa*(x1*x1+y1*y1);
double q=sqrt(1.-2.*kappa*D0);
double s0=(x1*cos(t.phi())+y1*sin(t.phi()))/q;
// double s1;
if (fabs(kappa*s0)>0.001){
//s1=asin(kappa*s0)/kappa;
}else{
//double ks02=(kappa*s0)*(kappa*s0);
//s1=s0*(1.+ks02/6.+3./40.*ks02*ks02+5./112.*pow(ks02,3));
}
// sp.ddcap=-2.*D0/(1.+q);
//double zdcap=t.vz()-s1/tan(t.theta());
}
//
// collect some info on hits and clusters
Fill(h,"nbarrelLayers_"+ttype,t.hitPattern().pixelBarrelLayersWithMeasurement());
Fill(h,"nPxLayers_"+ttype,t.hitPattern().pixelLayersWithMeasurement());
Fill(h,"nSiLayers_"+ttype,t.hitPattern().trackerLayersWithMeasurement());
Fill(h,"expectedInner_"+ttype,t.trackerExpectedHitsInner().numberOfHits());
Fill(h,"expectedOuter_"+ttype,t.trackerExpectedHitsOuter().numberOfHits());
Fill(h,"trackAlgo_"+ttype,t.algo());
Fill(h,"trackQuality_"+ttype,t.qualityMask());
//
int longesthit=0, nbarrel=0;
for(trackingRecHit_iterator hit=t.recHitsBegin(); hit!=t.recHitsEnd(); hit++){
if ((**hit).isValid() && (**hit).geographicalId().det() == DetId::Tracker ){
bool barrel = DetId((**hit).geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelBarrel);
//bool endcap = DetId::DetId((**hit).geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelEndcap);
if (barrel){
const SiPixelRecHit *pixhit = dynamic_cast<const SiPixelRecHit*>( &(**hit));
edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const& clust = (*pixhit).cluster();
if (clust.isNonnull()) {
nbarrel++;
if (clust->sizeY()-longesthit>0) longesthit=clust->sizeY();
if (clust->sizeY()>20.){
Fill(h,"lvseta_"+ttype,t.eta(), 19.9);
Fill(h,"lvstanlambda_"+ttype,tan(t.lambda()), 19.9);
}else{
Fill(h,"lvseta_"+ttype,t.eta(), float(clust->sizeY()));
Fill(h,"lvstanlambda_"+ttype,tan(t.lambda()), float(clust->sizeY()));
}
}
}
}
}
Fill(h,"nbarrelhits_"+ttype,float(nbarrel));
//-------------------------------------------------------------------
}
| std::vector< PrimaryVertexAnalyzer4PU::SimEvent > PrimaryVertexAnalyzer4PU::getSimEvents | ( | edm::Handle< TrackingParticleCollection > | TPCollectionH, |
| edm::Handle< TrackingVertexCollection > | TVCollectionH, | ||
| edm::Handle< edm::View< reco::Track > > | trackCollectionH | ||
| ) | [private] |
Definition at line 1578 of file PrimaryVertexAnalyzer4PU.cc.
References abs, EncodedEventId::bunchCrossing(), gather_cfg::cout, alignCSCRings::e, event(), EncodedEventId::event(), PrimaryVertexAnalyzer4PU::SimEvent::eventId, TrackingVertex::eventId(), edm::RefToBase< T >::get(), edm::Ref< C, T, F >::get(), i, errorMatrix2Lands_multiChannel::id, TransientVertex::isValid(), TrackingVertex::position(), TransientVertex::position(), funct::pow(), edm::Handle< T >::product(), reco::TransientTrack::setBeamSpot(), edm::View< T >::size(), mathSSE::sqrt(), lumiQTWidget::t, v, AdaptiveVertexFitter::vertex(), PrimaryVertexAnalyzer4PU::SimEvent::x, PV3DBase< T, PVType, FrameType >::x(), x, PrimaryVertexAnalyzer4PU::SimEvent::y, PV3DBase< T, PVType, FrameType >::y(), detailsBasic3DVector::y, PrimaryVertexAnalyzer4PU::SimEvent::z, PV3DBase< T, PVType, FrameType >::z(), and z.
Referenced by analyze().
{
const TrackingParticleCollection* simTracks = TPCollectionH.product();
const View<Track> tC = *(trackCollectionH.product());
vector<SimEvent> simEvt;
map<EncodedEventId, unsigned int> eventIdToEventMap;
map<EncodedEventId, unsigned int>::iterator id;
bool dumpTP=false;
for(TrackingParticleCollection::const_iterator it=simTracks->begin(); it!=simTracks->end(); it++){
if( fabs(it->parentVertex().get()->position().z())>100.) continue; // skip funny entries @ -13900
unsigned int event=0; //note, this is no longer the same as eventId().event()
id=eventIdToEventMap.find(it->eventId());
if (id==eventIdToEventMap.end()){
// new event here
SimEvent e;
e.eventId=it->eventId();
event=simEvt.size();
const TrackingVertex *parentVertex= it->parentVertex().get();
if(DEBUG_){
cout << "getSimEvents: ";
cout << it->eventId().bunchCrossing() << "," << it->eventId().event()
<< " z="<< it->vz() << " "
<< parentVertex->eventId().bunchCrossing() << "," <<parentVertex->eventId().event()
<< " z=" << parentVertex->position().z() << endl;
}
if (it->eventId()==parentVertex->eventId()){
//e.x=it->vx(); e.y=it->vy(); e.z=it->vz();// wrong ???
e.x=parentVertex->position().x();
e.y=parentVertex->position().y();
e.z=parentVertex->position().z();
if(e.z<-100){
dumpTP=true;
}
}else{
e.x=0; e.y=0; e.z=-88.;
}
simEvt.push_back(e);
eventIdToEventMap[e.eventId]=event;
}else{
event=id->second;
}
simEvt[event].tp.push_back(&(*it));
if( (abs(it->eta())<2.5) && (it->charge()!=0) ){
simEvt[event].sumpt2+=pow(it->pt(),2); // should keep track of decays ?
simEvt[event].sumpt+=it->pt();
}
}
if(dumpTP){
for(TrackingParticleCollection::const_iterator it=simTracks->begin(); it!=simTracks->end(); it++){
std::cout << *it << std::endl;
}
}
for(View<Track>::size_type i=0; i<tC.size(); ++i) {
RefToBase<Track> track(trackCollectionH, i);
TrackingParticleRef tpr;
if( truthMatchedTrack(track,tpr)){
if( eventIdToEventMap.find(tpr->eventId())==eventIdToEventMap.end() ){ cout << "Bug in getSimEvents" << endl; break; }
z2tp_[track.get()->vz()]=tpr;
unsigned int event=eventIdToEventMap[tpr->eventId()];
double ipsig=0,ipdist=0;
const TrackingVertex *parentVertex= tpr->parentVertex().get();
double vx=parentVertex->position().x(); // problems with tpr->vz()
double vy=parentVertex->position().y();
double vz=parentVertex->position().z();
double d=sqrt(pow(simEvt[event].x-vx,2)+pow(simEvt[event].y-vy,2)+pow(simEvt[event].z-vz,2))*1.e4;
ipdist=d;
double dxy=track->dxy(vertexBeamSpot_.position());
ipsig=dxy/track->dxyError();
TransientTrack t = theB_->build(tC[i]);
t.setBeamSpot(vertexBeamSpot_);
if (theTrackFilter(t)){
simEvt[event].tk.push_back(t);
if(ipdist<5){simEvt[event].tkprim.push_back(t);}
if(ipsig<5){simEvt[event].tkprimsel.push_back(t);}
}
}
}
AdaptiveVertexFitter theFitter;
cout << "SimEvents " << simEvt.size() << endl;
for(unsigned int i=0; i<simEvt.size(); i++){
if(simEvt[i].tkprim.size()>0){
getTc(simEvt[i].tkprimsel, simEvt[i].Tc, simEvt[i].chisq, simEvt[i].dzmax, simEvt[i].dztrim, simEvt[i].m4m2);
TransientVertex v = theFitter.vertex(simEvt[i].tkprim, vertexBeamSpot_);
if (v.isValid()){
simEvt[i].xfit=v.position().x();
simEvt[i].yfit=v.position().y();
simEvt[i].zfit=v.position().z();
// if (simEvt[i].z<-80.){simEvt[i].z=v.position().z();} // for now
}
}
if(DEBUG_){
cout << i <<" ) nTP=" << simEvt[i].tp.size()
<< " z=" << simEvt[i].z
<< " recTrks=" << simEvt[i].tk.size()
<< " recTrksPrim=" << simEvt[i].tkprim.size()
<< " zfit=" << simEvt[i].zfit
<< endl;
}
}
return simEvt;
}
| std::vector< PrimaryVertexAnalyzer4PU::simPrimaryVertex > PrimaryVertexAnalyzer4PU::getSimPVs | ( | const edm::Handle< edm::HepMCProduct > | evtMC | ) | [private] |
Definition at line 1710 of file PrimaryVertexAnalyzer4PU.cc.
References gather_cfg::cout, DEBUG_, alignCSCRings::e, spr::find(), hsimPV, isCharged(), isFinalstateParticle(), m, NULL, parents, pos, and verbose_.
Referenced by analyze().
{
if(DEBUG_){std::cout << "getSimPVs HepMC " << std::endl;}
std::vector<PrimaryVertexAnalyzer4PU::simPrimaryVertex> simpv;
const HepMC::GenEvent* evt=evtMC->GetEvent();
if (evt) {
// std::cout << "process id " <<evt->signal_process_id()<<std::endl;
// std::cout <<"signal process vertex "<< ( evt->signal_process_vertex() ?
// evt->signal_process_vertex()->barcode() : 0 ) <<std::endl;
// 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();
// if (pos.t()>0) { continue;} // skip secondary vertices, doesn't work for some samples
if (fabs(pos.z())>1000) continue; // skip funny junk vertices
bool hasMotherVertex=false;
//std::cout << "mothers" << 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;}
//std::cout << "\t"; (*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
//std::cout << "this is a new vertex" << sv.x << " " << sv.y << " " << sv.z <<std::endl;
simpv.push_back(sv);
vp=&simpv.back();
// }else{
// std::cout << "this is not a new vertex" << std::endl;
}
// store the gen vertex barcode with this simpv
vp->genVertex.push_back((*vitr)->barcode());
// collect final state descendants and sum up momenta etc
for ( HepMC::GenVertex::particle_iterator
daughter = (*vitr)->particles_begin(HepMC::descendants);
daughter != (*vitr)->particles_end(HepMC::descendants);
++daughter ) {
//std::cout << "checking daughter type " << (*daughter)->pdg_id() << " final :" <<isFinalstateParticle(*daughter) << std::endl;
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();
//std::cout << "adding particle to primary " << m.px()<< " " << m.py() << " " << m.pz() << std::endl;
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());
// count relevant particles
if ( (m.perp()>0.2) && (fabs(m.pseudoRapidity())<2.5) && isCharged( *daughter ) ){
vp->nGenTrk++;
}
hsimPV["rapidity"]->Fill(m.pseudoRapidity());
if( (m.perp()>0.8) && isCharged( *daughter ) ){
hsimPV["chRapidity"]->Fill(m.pseudoRapidity());
}
hsimPV["pt"]->Fill(m.perp());
}//new final state particle for this vertex
}//daughter is a final state particle
}
//idx++;
}
}
if(verbose_){
cout << "------- PrimaryVertexAnalyzer4PU simPVs -------" << endl;
for(std::vector<simPrimaryVertex>::iterator v0=simpv.begin();
v0!=simpv.end(); v0++){
cout << "z=" << v0->z
<< " px=" << v0->ptot.px()
<< " py=" << v0->ptot.py()
<< " pz=" << v0->ptot.pz()
<< " pt2="<< v0->ptsq
<< endl;
}
cout << "-----------------------------------------------" << endl;
}
return simpv;
}
| std::vector<simPrimaryVertex> PrimaryVertexAnalyzer4PU::getSimPVs | ( | const edm::Handle< edm::HepMCProduct > | evt, |
| const edm::Handle< edm::SimVertexContainer > | simVtxs, | ||
| const edm::Handle< edm::SimTrackContainer > | simTrks | ||
| ) | [private] |
| std::vector< PrimaryVertexAnalyzer4PU::simPrimaryVertex > PrimaryVertexAnalyzer4PU::getSimPVs | ( | const edm::Handle< TrackingVertexCollection > | tVC | ) | [private] |
Definition at line 1831 of file PrimaryVertexAnalyzer4PU.cc.
References begin, gather_cfg::cout, DEBUG_, alignCSCRings::e, end, event(), PrimaryVertexAnalyzer4PU::simPrimaryVertex::eventId, NULL, benchmark_cfg::pdgId, position, and v.
{
std::vector<PrimaryVertexAnalyzer4PU::simPrimaryVertex> simpv;
//std::cout <<"number of vertices " << tVC->size() << std::endl;
if(DEBUG_){std::cout << "getSimPVs TrackingVertexCollection " << std::endl;}
for (TrackingVertexCollection::const_iterator v = tVC -> begin(); v != tVC -> end(); ++v) {
if(DEBUG_){
std::cout << (v->eventId()).event() << v -> position() << v->g4Vertices().size() <<" " <<v->genVertices().size() << " t=" <<v->position().t()*1.e12 <<" ==0:" <<(v->position().t()>0) << std::endl;
for( TrackingVertex::g4v_iterator gv=v->g4Vertices_begin(); gv!=v->g4Vertices_end(); gv++){
std::cout << *gv << std::endl;
}
std::cout << "----------" << std::endl;
}
// bool hasMotherVertex=false;
if ((unsigned int) v->eventId().event()<simpv.size()) continue;
//if (v->position().t()>0) { continue;} // skip secondary vertices (obsolete + doesn't work)
if (fabs(v->position().z())>1000) continue; // skip funny junk vertices
// could be a new vertex, check all primaries found so far to avoid multiple entries
const double mm=1.0; // for tracking vertices
simPrimaryVertex sv(v->position().x()*mm,v->position().y()*mm,v->position().z()*mm);
//cout << "sv: " << (v->eventId()).event() << endl;
sv.eventId=v->eventId();
for (TrackingParticleRefVector::iterator iTrack = v->daughterTracks_begin(); iTrack != v->daughterTracks_end(); ++iTrack){
//cout <<((**iTrack).eventId()).event() << " "; // an iterator of Refs, dereference twice. Cool eyh?
sv.eventId=(**iTrack).eventId();
}
//cout <<endl;
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( (sv.eventId==v0->eventId) && (fabs(sv.x-v0->x)<1e-5) && (fabs(sv.y-v0->y)<1e-5) && (fabs(sv.z-v0->z)<1e-5)){
vp=&(*v0);
break;
}
}
if(!vp){
// this is a new vertex
if(DEBUG_){std::cout << "this is a new vertex " << sv.eventId.event() << " " << sv.x << " " << sv.y << " " << sv.z <<std::endl;}
simpv.push_back(sv);
vp=&simpv.back();
}else{
if(DEBUG_){std::cout << "this is not a new vertex" << sv.x << " " << sv.y << " " << sv.z <<std::endl;}
}
// Loop over daughter track(s)
if(DEBUG_){
for (TrackingVertex::tp_iterator iTP = v -> daughterTracks_begin(); iTP != v -> daughterTracks_end(); ++iTP) {
std::cout << " Daughter momentum: " << (*(*iTP)).momentum();
std::cout << " Daughter type " << (*(*iTP)).pdgId();
std::cout << std::endl;
}
}
}
if(DEBUG_){
cout << "------- PrimaryVertexAnalyzer4PU simPVs from TrackingVertices -------" << endl;
for(std::vector<simPrimaryVertex>::iterator v0=simpv.begin();
v0!=simpv.end(); v0++){
cout << "z=" << v0->z << " event=" << v0->eventId.event() << endl;
}
cout << "-----------------------------------------------" << endl;
}
return simpv;
}
| std::vector< PrimaryVertexAnalyzer4PU::SimPart > PrimaryVertexAnalyzer4PU::getSimTrkParameters | ( | edm::Handle< edm::SimTrackContainer > & | simTrks, |
| edm::Handle< edm::SimVertexContainer > & | simVtcs, | ||
| double | simUnit = 1.0 |
||
| ) | [private] |
Definition at line 733 of file PrimaryVertexAnalyzer4PU.cc.
References funct::cos(), gather_cfg::cout, PrimaryVertexAnalyzer4PU::SimPart::ddcap, alignCSCRings::e, fBfield_, reco::TrackBase::i_dsz, reco::TrackBase::i_dxy, reco::TrackBase::i_lambda, reco::TrackBase::i_phi, reco::TrackBase::i_qoverp, M_PI, AlCaHLTBitMon_ParallelJobs::p, PrimaryVertexAnalyzer4PU::SimPart::par, PrimaryVertexAnalyzer4PU::SimPart::pdg, funct::pow(), lumiQueryAPI::q, funct::sin(), mathSSE::sqrt(), lumiQTWidget::t, funct::tan(), PrimaryVertexAnalyzer4PU::SimPart::type, v, verbose_, PrimaryVertexAnalyzer4PU::SimPart::xvtx, PrimaryVertexAnalyzer4PU::SimPart::yvtx, PrimaryVertexAnalyzer4PU::SimPart::zdcap, and PrimaryVertexAnalyzer4PU::SimPart::zvtx.
Referenced by analyze().
{
std::vector<SimPart > tsim;
if(simVtcs->begin()==simVtcs->end()){
if(verbose_){
cout << " PrimaryVertexAnalyzer4PU::getSimTrkParameters no simvtcs" << endl;
}
return tsim;
}
if(verbose_){
cout << " PrimaryVertexAnalyzer4PU::getSimTrkParameters simVtcs n=" << simVtcs->size() << endl;
cout << " PrimaryVertexAnalyzer4PU::getSimTrkParameters 1st position" << setw(8) << setprecision(4) << simVtcs->begin()->position() << endl;
}
double t0=simVtcs->begin()->position().e();
for(edm::SimTrackContainer::const_iterator t=simTrks->begin();
t!=simTrks->end(); ++t){
if (t->noVertex()){
std::cout << "simtrk has no vertex" << std::endl;
}else{
// get the vertex position
//HepLorentzVector v=(*simVtcs)[t->vertIndex()].position();
math::XYZTLorentzVectorD v((*simVtcs)[t->vertIndex()].position().x(),
(*simVtcs)[t->vertIndex()].position().y(),
(*simVtcs)[t->vertIndex()].position().z(),
(*simVtcs)[t->vertIndex()].position().e());
int pdgCode=t->type();
if( pdgCode==-99 ){
// such entries cause crashes, no idea what they are
std::cout << "funny particle skipped , code=" << pdgCode << std::endl;
}else{
double Q=0; //double Q=HepPDT::theTable().getParticleData(pdgCode)->charge();
if ((pdgCode==11)||(pdgCode==13)||(pdgCode==15)||(pdgCode==-211)||(pdgCode==-2212)||(pdgCode==-321)||(pdgCode==-3222)){Q=-1;}
else if((pdgCode==-11)||(pdgCode==-13)||(pdgCode==-15)||(pdgCode==211)||(pdgCode==2212)||(pdgCode==321)||(pdgCode==3222)){Q=1;}
else {
//std::cout << pdgCode << " " <<std::endl;
}
math::XYZTLorentzVectorD p(t->momentum().x(),t->momentum().y(),t->momentum().z(),t->momentum().e());
if ( (Q != 0) && (p.pt()>0.1) && (fabs(t->momentum().eta())<3.0)
&& fabs(v.z()*simUnit<20) && (sqrt(v.x()*v.x()+v.y()*v.y())<10.)){
double x0=v.x()*simUnit;
double y0=v.y()*simUnit;
double z0=v.z()*simUnit;
double kappa=-Q*0.002998*fBfield_/p.pt();
double D0=x0*sin(p.phi())-y0*cos(p.phi())-0.5*kappa*(x0*x0+y0*y0);
double q=sqrt(1.-2.*kappa*D0);
double s0=(x0*cos(p.phi())+y0*sin(p.phi()))/q;
double s1;
if (fabs(kappa*s0)>0.001){
s1=asin(kappa*s0)/kappa;
}else{
double ks02=(kappa*s0)*(kappa*s0);
s1=s0*(1.+ks02/6.+3./40.*ks02*ks02+5./112.*pow(ks02,3));
}
SimPart sp;//ParameterVector par;
sp.par[reco::TrackBase::i_qoverp] = Q/p.P();
sp.par[reco::TrackBase::i_lambda] = M_PI/2.-p.theta();
sp.par[reco::TrackBase::i_phi] = p.phi()-asin(kappa*s0);
sp.par[reco::TrackBase::i_dxy] = -2.*D0/(1.+q);
sp.par[reco::TrackBase::i_dsz] = z0*sin(p.theta())-s1*cos(p.theta());
sp.pdg=pdgCode;
if (v.t()-t0<1e-15){
sp.type=0; // primary
}else{
sp.type=1; //secondary
}
// now get zpca (get perigee wrt beam)
double x1=x0-0.033; double y1=y0-0.; // FIXME how do we get the simulated beam position?
D0=x1*sin(p.phi())-y1*cos(p.phi())-0.5*kappa*(x1*x1+y1*y1);
q=sqrt(1.-2.*kappa*D0);
s0=(x1*cos(p.phi())+y1*sin(p.phi()))/q;
if (fabs(kappa*s0)>0.001){
s1=asin(kappa*s0)/kappa;
}else{
double ks02=(kappa*s0)*(kappa*s0);
s1=s0*(1.+ks02/6.+3./40.*ks02*ks02+5./112.*pow(ks02,3));
}
sp.ddcap=-2.*D0/(1.+q);
sp.zdcap=z0-s1/tan(p.theta());
sp.zvtx=z0;
sp.xvtx=x0;
sp.yvtx=y0;
tsim.push_back(sp);
}
}
}// has vertex
}//for loop
return tsim;
}
| void PrimaryVertexAnalyzer4PU::getTc | ( | const std::vector< reco::TransientTrack > & | , |
| double & | , | ||
| double & | , | ||
| double & | , | ||
| double & | , | ||
| double & | |||
| ) | [private] |
Definition at line 1466 of file PrimaryVertexAnalyzer4PU.cc.
References a, b, SiPixelRawToDigiRegional_cfi::beamSpot, position, funct::pow(), mathSSE::sqrt(), funct::tan(), w(), and z.
Referenced by analyzeVertexCollectionTP().
{
if (tracks.size()<2){ Tc=-1; chsq=-1; dzmax=-1; dztrim=-1; m4m2=-1; return;}
double sumw=0, sumwz=0, sumww=0,sumwwz=0,sumwwzz=0;
double zmin=1e10, zmin1=1e10, zmax1=-1e10, zmax=-1e10;
double m4=0,m3=0,m2=0,m1=0,m0=0;
for(vector<reco::TransientTrack>::const_iterator it=tracks.begin(); it!=tracks.end(); it++){
double tantheta=tan(((*it).stateAtBeamLine().trackStateAtPCA()).momentum().theta());
reco::BeamSpot beamspot=(it->stateAtBeamLine()).beamSpot();
double z=((*it).stateAtBeamLine().trackStateAtPCA()).position().z();
double dz2= pow((*it).track().dzError(),2)+pow(beamspot.BeamWidthX()/tantheta,2);
// t.dz2= pow((*it).track().dzError(),2) + pow(wxy0/tantheta,2) + 1./(1.+exp(pow(t.ip/t.dip,2)-pow(2.)))*pow(ip/tantheta,2);
double w=1./dz2; // take p=1
sumw += w;
sumwz += w*z;
sumwwz += w*w*z;;
sumwwzz += w*w*z*z;
sumww += w*w;
m0 += w;
m1 += w*z;
m2 += w*z*z;
m3 += w*z*z*z;
m4 += w*z*z*z*z;
if(dz2<pow(0.1,2)){
if(z<zmin1){zmin1=z;} if(z<zmin){zmin1=zmin; zmin=z;}
if(z>zmax1){zmax1=z;} if(z>zmax){zmax1=zmax; zmax=z;}
}
}
double z=sumwz/sumw;
double a=sumwwzz-2*z*sumwwz+z*z*sumww;
double b=sumw;
if(tracks.size()>1){
chsq=(m2-m0*z*z)/(tracks.size()-1);
Tc=2.*a/b;
m4m2=sqrt((m4-4*m3*z+6*m2*z*z-3*m1*z*z*z+m0*z*z*z*z)/(m2-2*m1*z+z*z*m0));
}else{
chsq=0;
Tc=0;
m4m2=0;
}
dzmax=zmax-zmin;
dztrim=zmax1-zmin1;// truncated
}
| std::vector< edm::RefToBase< reco::Track > > PrimaryVertexAnalyzer4PU::getTruthMatchedVertexTracks | ( | const reco::Vertex & | v | ) | [private] |
Definition at line 1551 of file PrimaryVertexAnalyzer4PU.cc.
References b, reco::Vertex::tracks_begin(), reco::Vertex::tracks_end(), and truthMatchedTrack().
Referenced by analyzeVertexCollectionTP(), and printEventSummary().
{
std::vector< edm::RefToBase<reco::Track> > b;
TrackingParticleRef tpr;
for(trackit_t tv=v.tracks_begin(); tv!=v.tracks_end(); tv++){
edm::RefToBase<reco::Track> track=*tv;
if (truthMatchedTrack(track, tpr)){
b.push_back(*tv);
}
}
return b;
}
| void PrimaryVertexAnalyzer4PU::history | ( | const edm::Handle< edm::View< reco::Track > > & | tracks, |
| const size_t | trackindex = 10000 |
||
| ) | [private] |
| bool PrimaryVertexAnalyzer4PU::isCharged | ( | const HepMC::GenParticle * | p | ) | [private] |
Definition at line 990 of file PrimaryVertexAnalyzer4PU.cc.
References pdt_.
Referenced by getSimPVs().
| bool PrimaryVertexAnalyzer4PU::isFinalstateParticle | ( | const HepMC::GenParticle * | p | ) | [private] |
Definition at line 985 of file PrimaryVertexAnalyzer4PU.cc.
Referenced by getSimPVs().
| bool PrimaryVertexAnalyzer4PU::isResonance | ( | const HepMC::GenParticle * | p | ) | [private] |
Definition at line 979 of file PrimaryVertexAnalyzer4PU.cc.
References abs, alignCSCRings::e, and pdt_.
| bool PrimaryVertexAnalyzer4PU::match | ( | const ParameterVector & | a, |
| const ParameterVector & | b | ||
| ) | [static, private] |
Definition at line 962 of file PrimaryVertexAnalyzer4PU.cc.
Referenced by matchRecTracksToVertex(), and supf().
{
double dqoverp =a(0)-b(0);
double dlambda =a(1)-b(1);
double dphi =a(2)-b(2);
double dsz =a(4)-b(4);
if (dphi>M_PI){ dphi-=M_2_PI; }else if(dphi<-M_PI){dphi+=M_2_PI;}
// return ( (fabs(dqoverp)<0.2) && (fabs(dlambda)<0.02) && (fabs(dphi)<0.04) && (fabs(dsz)<0.1) );
return ( (fabs(dqoverp)<0.2) && (fabs(dlambda)<0.02) && (fabs(dphi)<0.04) && (fabs(dsz)<1.0) );
}
| void PrimaryVertexAnalyzer4PU::matchRecTracksToVertex | ( | simPrimaryVertex & | pv, |
| const std::vector< SimPart > & | tsim, | ||
| const edm::Handle< reco::TrackCollection > & | recTrks | ||
| ) | [private] |
Definition at line 1418 of file PrimaryVertexAnalyzer4PU.cc.
References gather_cfg::cout, match(), PrimaryVertexAnalyzer4PU::simPrimaryVertex::matchedRecTrackIndex, PrimaryVertexAnalyzer4PU::simPrimaryVertex::nMatchedTracks, AlCaHLTBitMon_ParallelJobs::p, alignCSCRings::s, and lumiQTWidget::t.
Referenced by analyze().
{
// find all recTracks that belong to this simulated vertex (not just the ones that are used in
// matching recVertex)
std::cout << "dump rec tracks: " << std::endl;
int irec=0;
for(reco::TrackCollection::const_iterator t=recTrks->begin();
t!=recTrks->end(); ++t){
reco::TrackBase::ParameterVector p = t->parameters();
std::cout << irec++ << ") " << p << std::endl;
}
std::cout << "match sim tracks: " << std::endl;
pv.matchedRecTrackIndex.clear();
pv.nMatchedTracks=0;
int isim=0;
for(std::vector<SimPart>::const_iterator s=tsim.begin();
s!=tsim.end(); ++s){
std::cout << isim++ << " " << s->par;// << std::endl;
int imatch=-1;
int irec=0;
for(reco::TrackCollection::const_iterator t=recTrks->begin();
t!=recTrks->end(); ++t){
reco::TrackBase::ParameterVector p = t->parameters();
if (match(s->par,p)){ imatch=irec; }
irec++;
}
pv.matchedRecTrackIndex.push_back(imatch);
if(imatch>-1){
pv.nMatchedTracks++;
std::cout << " matched to rec trk" << imatch << std::endl;
}else{
std::cout << " not matched" << std::endl;
}
}
}
| bool PrimaryVertexAnalyzer4PU::matchVertex | ( | const simPrimaryVertex & | vsim, |
| const reco::Vertex & | vrec | ||
| ) | [private] |
Definition at line 974 of file PrimaryVertexAnalyzer4PU.cc.
References simUnit_, reco::Vertex::z(), PrimaryVertexAnalyzer4PU::simPrimaryVertex::z, and zmatch_.
| std::string PrimaryVertexAnalyzer4PU::particleString | ( | int | ) | const [private] |
| void PrimaryVertexAnalyzer4PU::printEventSummary | ( | std::map< std::string, TH1 * > & | h, |
| const edm::Handle< reco::VertexCollection > | recVtxs, | ||
| const edm::Handle< reco::TrackCollection > | recTrks, | ||
| std::vector< SimEvent > & | simEvt, | ||
| const std::string | message | ||
| ) | [private] |
Definition at line 2222 of file PrimaryVertexAnalyzer4PU.cc.
References SiPixelRawToDigiRegional_cfi::beamSpot, gather_cfg::cout, reco::TrackBase::dzError(), reco::TrackBase::eta(), eventcounter_, Fill(), edm::Ref< C, T, F >::get(), getTruthMatchedVertexTracks(), lt, n, reco::TrackBase::phi(), funct::pow(), reco::TrackBase::pt(), findQualityFiles::size, mathSSE::sqrt(), funct::tan(), theTrackFilter, PrimaryVertexAnalyzer4PU::SimEvent::tk, PrimaryVertexAnalyzer4PU::SimEvent::tkprim, reco::Vertex::tracks_begin(), reco::Vertex::tracks_end(), reco::Vertex::trackWeight(), v, reco::TrackBase::vz(), PrimaryVertexAnalyzer4PU::SimEvent::z, z, and z2tp_.
Referenced by analyzeVertexCollectionTP().
{
// make a readable summary of the vertex finding if the TrackingParticles are availabe
if (simEvt.size()==0){return;}
// sort vertices in z ... for nicer printout
vector< pair<double,unsigned int> > zrecv;
for(unsigned int idx=0; idx<recVtxs->size(); idx++){
if ( (recVtxs->at(idx).ndof()<0) || (recVtxs->at(idx).chi2()<=0) ) continue; // skip clusters
zrecv.push_back( make_pair(recVtxs->at(idx).z(),idx) );
}
stable_sort(zrecv.begin(),zrecv.end(),lt);
// same for simulated vertices
vector< pair<double,unsigned int> > zsimv;
for(unsigned int idx=0; idx<simEvt.size(); idx++){
zsimv.push_back(make_pair(simEvt[idx].z, idx));
}
stable_sort(zsimv.begin(), zsimv.end(),lt);
cout << "---------------------------" << endl;
cout << "event counter = " << eventcounter_ << " " << message << endl;
cout << "---------------------------" << endl;
cout << " z[cm] rec --> ";
cout.precision(4);
for(vector< pair<double,unsigned int> >::iterator itrec=zrecv.begin(); itrec!=zrecv.end(); itrec++){
cout << setw(7) << fixed << itrec->first;
if (itrec->second==0){cout << "*" ;}else{cout << " " ;}
}
cout << endl;
cout << " ";
for(vector< pair<double,unsigned int> >::iterator itrec=zrecv.begin(); itrec!=zrecv.end(); itrec++){
cout << setw(7) << fixed << recVtxs->at(itrec->second).tracksSize();
if (itrec->second==0){cout << "*" ;}else{cout << " " ;}
}
cout << " rec tracks" << endl;
cout << " ";
map<unsigned int, int> truthMatchedVertexTracks;
for(vector< pair<double,unsigned int> >::iterator itrec=zrecv.begin(); itrec!=zrecv.end(); itrec++){
truthMatchedVertexTracks[itrec->second]=getTruthMatchedVertexTracks(recVtxs->at(itrec->second)).size();
cout << setw(7) << fixed << truthMatchedVertexTracks[itrec->second];
if (itrec->second==0){cout << "*" ;}else{cout << " " ;}
}
cout << " truth matched " << endl;
cout << "sim ------- trk prim ----" << endl;
map<unsigned int, unsigned int> rvmatch; // reco vertex matched to sim vertex (sim to rec)
map<unsigned int, double > nmatch; // highest number of truth-matched tracks of ev found in a recvtx
map<unsigned int, double > purity; // highest purity of a rec vtx (i.e. highest number of tracks from the same simvtx)
map<unsigned int, double > wpurity; // same for the sum of weights
for(vector< pair<double,unsigned int> >::iterator itrec=zrecv.begin(); itrec!=zrecv.end(); itrec++){
purity[itrec->second]=0.;
wpurity[itrec->second]=0.;
}
for(vector< pair<double,unsigned int> >::iterator itsim=zsimv.begin(); itsim!=zsimv.end(); itsim++){
SimEvent* ev =&(simEvt[itsim->second]);
cout.precision(4);
if (itsim->second==0){
cout << setw(8) << fixed << ev->z << ")*" << setw(5) << ev->tk.size() << setw(5) << ev->tkprim.size() << " | ";
}else{
cout << setw(8) << fixed << ev->z << ") " << setw(5) << ev->tk.size() << setw(5) << ev->tkprim.size() << " | ";
}
nmatch[itsim->second]=0; // highest number of truth-matched tracks of ev found in a recvtx
double matchpurity=0,matchwpurity=0;
for(vector< pair<double,unsigned int> >::iterator itrec=zrecv.begin(); itrec!=zrecv.end(); itrec++){
const reco::Vertex *v = &(recVtxs->at(itrec->second));
// count tracks found in both, sim and rec
double n=0,wt=0;
for(vector<TransientTrack>::iterator te=ev->tk.begin(); te!=ev->tk.end(); te++){
const reco::Track& RTe=te->track();
for(trackit_t tv=v->tracks_begin(); tv!=v->tracks_end(); tv++){
const reco::Track & RTv=*(tv->get());
if(RTe.vz()==RTv.vz()) {n++; wt+=v->trackWeight(*tv);}
}
}
cout << setw(7) << int(n)<< " ";
if (n > nmatch[itsim->second]){
nmatch[itsim->second]=n;
rvmatch[itsim->second]=itrec->second;
matchpurity=n/truthMatchedVertexTracks[itrec->second];
matchwpurity=wt/truthMatchedVertexTracks[itrec->second];
}
if(n > purity[itrec->second]){
purity[itrec->second]=n;
}
if(wt > wpurity[itrec->second]){
wpurity[itrec->second]=wt;
}
}// end of reco vertex loop
cout << " | ";
if (nmatch[itsim->second]>0 ){
if(matchpurity>0.5){
cout << "found ";
}else{
cout << "merged ";
}
cout << " max eff. = " << setw(8) << nmatch[itsim->second]/ev->tk.size() << " p=" << matchpurity << " w=" << matchwpurity << endl;
}else{
if(ev->tk.size()==0){
cout << " invisible" << endl;
}else if (ev->tk.size()==1){
cout << "single track " << endl;
}else{
cout << "lost " << endl;
}
}
}
cout << "---------------------------" << endl;
// the purity of the reconstructed vertex
cout << " purity ";
for(vector< pair<double,unsigned int> >::iterator itrec=zrecv.begin(); itrec!=zrecv.end(); itrec++){
cout << setw(7) << fixed << purity[itrec->second]/truthMatchedVertexTracks[itrec->second];
if (itrec->second==0){cout << "*" ;}else{cout << " " ;}
}
cout << endl;
// // classification of reconstructed vertex fake/real
// cout << " ";
// for(vector< pair<double,unsigned int> >::iterator itrec=zrecv.begin(); itrec!=zrecv.end(); itrec++){
// cout << setw(7) << fixed << purity[itrec->second]/truthMatchedVertexTracks[itrec->second];
// if (itrec->second==0){cout << "*" ;}else{cout << " " ;}
// }
// cout << endl;
cout << "---------------------------" << endl;
// list problematic tracks
for(vector< pair<double,unsigned int> >::iterator itsim=zsimv.begin(); itsim!=zsimv.end(); itsim++){
SimEvent* ev =&(simEvt[itsim->second]);
for(vector<TransientTrack>::iterator te=ev->tk.begin(); te!=ev->tk.end(); te++){
const reco::Track& RTe=te->track();
int ivassign=-1; // will become the index of the vertex to which a track was assigned
for(vector< pair<double,unsigned int> >::iterator itrec=zrecv.begin(); itrec!=zrecv.end(); itrec++){
const reco::Vertex *v = &(recVtxs->at(itrec->second));
for(trackit_t tv=v->tracks_begin(); tv!=v->tracks_end(); tv++){
const reco::Track & RTv=*(tv->get());
if(RTe.vz()==RTv.vz()) {ivassign=itrec->second;}
}
}
double tantheta=tan((te->stateAtBeamLine().trackStateAtPCA()).momentum().theta());
reco::BeamSpot beamspot=(te->stateAtBeamLine()).beamSpot();
//double z=(te->stateAtBeamLine().trackStateAtPCA()).position().z();
double dz2= pow(RTe.dzError(),2)+pow(beamspot.BeamWidthX()/tantheta,2);
if(ivassign==(int)rvmatch[itsim->second]){
Fill(h,"correctlyassigned",RTe.eta(),RTe.pt());
Fill(h,"ptcat",RTe.pt());
Fill(h,"etacat",RTe.eta());
Fill(h,"phicat",RTe.phi());
Fill(h,"dzcat",sqrt(dz2));
}else{
Fill(h,"misassigned",RTe.eta(),RTe.pt());
Fill(h,"ptmis",RTe.pt());
Fill(h,"etamis",RTe.eta());
Fill(h,"phimis",RTe.phi());
Fill(h,"dzmis",sqrt(dz2));
cout << "vertex " << setw(8) << fixed << ev->z;
if (ivassign<0){
cout << " track lost ";
// for some clusterizers there shouldn't be any lost tracks,
// are there differences in the track selection?
}else{
cout << " track misassigned " << setw(8) << fixed << recVtxs->at(ivassign).z();
}
cout << " track z=" << setw(8) << fixed << RTe.vz() << "+/-" << RTe.dzError() << " pt=" << setw(8) << fixed<< RTe.pt() << " eta=" << setw(8) << fixed << RTe.eta()<< " sel=" <<theTrackFilter(*te);
//
//cout << " ztrack=" << te->track().vz();
TrackingParticleRef tpr = z2tp_[te->track().vz()];
double zparent=tpr->parentVertex().get()->position().z();
if(zparent==ev->z) {
cout << " prim";
}else{
cout << " sec ";
}
cout << " id=" << tpr->pdgId();
cout << endl;
//
}
}// next simvertex-track
}//next simvertex
cout << "---------------------------" << endl;
}
| void PrimaryVertexAnalyzer4PU::printPVTrks | ( | const edm::Handle< reco::TrackCollection > & | recTrks, |
| const edm::Handle< reco::VertexCollection > & | recVtxs, | ||
| std::vector< SimPart > & | tsim, | ||
| std::vector< SimEvent > & | simEvt, | ||
| const bool | selectedOnly = true |
||
| ) | [private] |
Definition at line 1283 of file PrimaryVertexAnalyzer4PU.cc.
References reco::BeamSpot::BeamWidthX(), reco::BeamSpot::BeamWidthY(), gather_cfg::cout, Measurement1D::error(), edm::Ref< C, T, F >::get(), reco::TrackBase::highPurity, i, TrackingVertex::position(), position, funct::pow(), reco::TransientTrack::setBeamSpot(), mathSSE::sqrt(), ntuplemaker::status, supf(), lumiQTWidget::t, funct::tan(), theB_, theta(), theTrackFilter, groupFilesInBlocks::tt, v, Measurement1D::value(), vertexBeamSpot_, reco::TrackBase::vz(), and z2tp_.
Referenced by analyze().
{
// make a printout of the tracks selected for PV reconstructions, show matching MC tracks, too
vector<TransientTrack> selTrks;
for(reco::TrackCollection::const_iterator t=recTrks->begin();
t!=recTrks->end(); ++t){
TransientTrack tt = theB_->build(&(*t)); tt.setBeamSpot(vertexBeamSpot_);
if( (!selectedOnly) || (selectedOnly && theTrackFilter(tt))){ selTrks.push_back(tt); }
}
if (selTrks.size()==0) return;
stable_sort(selTrks.begin(), selTrks.end(), recTrackLessZ);
// select tracks like for PV reconstruction and match them to sim tracks
reco::TrackCollection selRecTrks;
for(unsigned int i=0; i<selTrks.size(); i++){ selRecTrks.push_back(selTrks[i].track());}
int* rectosim=supf(tsim, selRecTrks);
// now dump in a format similar to the clusterizer
cout << "printPVTrks" << endl;
cout << "---- z +/- dz 1bet-m ip +/-dip pt phi eta";
if((tsim.size()>0)||(simEvt.size()>0)) {cout << " type pdg zvtx zdca dca zvtx zdca dsz";}
cout << endl;
cout.precision(4);
int isel=0;
for(unsigned int i=0; i<selTrks.size(); i++){
if (selectedOnly || (theTrackFilter(selTrks[i]))) {
cout << setw (3)<< isel;
isel++;
}else{
cout << " ";
}
// is this track in the tracklist of a recvtx ?
int vmatch=-1;
int iv=0;
for(reco::VertexCollection::const_iterator v=recVtxs->begin();
v!=recVtxs->end(); ++v){
if ( (v->ndof()<-1) || (v->chi2()<=0) ) continue; // skip clusters
for(trackit_t tv=v->tracks_begin(); tv!=v->tracks_end(); tv++){
const reco::Track & RTv=*(tv->get());
if(selTrks[i].track().vz()==RTv.vz()) {vmatch=iv;}
}
iv++;
}
double tz=(selTrks[i].stateAtBeamLine().trackStateAtPCA()).position().z();
double tantheta=tan((selTrks[i].stateAtBeamLine().trackStateAtPCA()).momentum().theta());
double tdz0= selTrks[i].track().dzError();
double tdz2= pow(selTrks[i].track().dzError(),2)+ (pow(vertexBeamSpot_.BeamWidthX(),2)+pow(vertexBeamSpot_.BeamWidthY(),2))/pow(tantheta,2);
if(vmatch>-1){
cout << "["<<setw(2)<<vmatch<<"]";
}else{
//int status=theTrackFilter.status(selTrks[i]);
int status=0;
if(status==0){
cout <<" ";
}else{
if(status&0x1){cout << "i";}else{cout << ".";};
if(status&0x2){cout << "c";}else{cout << ".";};
if(status&0x4){cout << "h";}else{cout << ".";};
if(status&0x8){cout << "q";}else{cout << ".";};
}
}
cout << setw (8) << fixed << setprecision(4)<< tz << " +/-" << setw (6)<< sqrt(tdz2);
// track quality and hit information, see DataFormats/TrackReco/interface/HitPattern.h
if(selTrks[i].track().quality(reco::TrackBase::highPurity)){ cout << " *";}else{cout <<" ";}
if(selTrks[i].track().hitPattern().hasValidHitInFirstPixelBarrel()){cout <<"+";}else{cout << "-";}
cout << setw(1) << selTrks[i].track().hitPattern().pixelBarrelLayersWithMeasurement();
cout << setw(1) << selTrks[i].track().hitPattern().pixelEndcapLayersWithMeasurement();
cout << setw(1) << hex << selTrks[i].track().hitPattern().trackerLayersWithMeasurement()-selTrks[i].track().hitPattern().pixelLayersWithMeasurement()<<dec;
cout << "-" << setw(1)<<hex <<selTrks[i].track().trackerExpectedHitsOuter().numberOfHits() << dec;
Measurement1D IP=selTrks[i].stateAtBeamLine().transverseImpactParameter();
cout << setw (8) << IP.value() << "+/-" << setw (6) << IP.error();
if(selTrks[i].track().ptError()<1){
cout << " " << setw(8) << setprecision(2) << selTrks[i].track().pt()*selTrks[i].track().charge();
}else{
cout << " " << setw(7) << setprecision(1) << selTrks[i].track().pt()*selTrks[i].track().charge() << "*";
//cout << "+/-" << setw(6)<< setprecision(2) << selTrks[i].track().ptError();
}
cout << " " << setw(5) << setprecision(2) << selTrks[i].track().phi()
<< " " << setw(5) << setprecision(2) << selTrks[i].track().eta() ;
// print MC info, if available
if(simEvt.size()>0){
reco::Track t=selTrks[i].track();
try{
TrackingParticleRef tpr = z2tp_[t.vz()];
const TrackingVertex *parentVertex= tpr->parentVertex().get();
//double vx=parentVertex->position().x(); // problems with tpr->vz()
//double vy=parentVertex->position().y(); work in progress
double vz=parentVertex->position().z();
cout << " " << tpr->eventId().event();
cout << " " << setw(5) << tpr->pdgId();
cout << " " << setw(8) << setprecision(4) << vz;
}catch (...){
cout << " not matched";
}
}else{
// cout << " " << rectosim[i];
if(rectosim[i]>0){
if(tsim[rectosim[i]].type==0){ cout << " prim " ;}else{cout << " sec ";}
cout << " " << setw(5) << tsim[rectosim[i]].pdg;
cout << " " << setw(8) << setprecision(4) << tsim[rectosim[i]].zvtx;
cout << " " << setw(8) << setprecision(4) << tsim[rectosim[i]].zdcap;
cout << " " << setw(8) << setprecision(4) << tsim[rectosim[i]].ddcap;
double zvtxpull=(tz-tsim[rectosim[i]].zvtx)/sqrt(tdz2);
cout << setw(6) << setprecision(1) << zvtxpull;
double zdcapull=(tz-tsim[rectosim[i]].zdcap)/tdz0;
cout << setw(6) << setprecision(1) << zdcapull;
double dszpull=(selTrks[i].track().dsz()-tsim[rectosim[i]].par[4])/selTrks[i].track().dszError();
cout << setw(6) << setprecision(1) << dszpull;
}
}
cout << endl;
}
delete rectosim;
}
| void PrimaryVertexAnalyzer4PU::printRecTrks | ( | const edm::Handle< reco::TrackCollection > & | recTrks | ) | [private] |
Definition at line 1197 of file PrimaryVertexAnalyzer4PU.cc.
References Reference_intrackfit_cff::barrel, reco::TrackBase::confirmed, gather_cfg::cout, Reference_intrackfit_cff::endcap, reco::TrackBase::goodIterative, reco::TrackBase::highPurity, i, edm::Ref< C, T, F >::isNonnull(), TrajectoryStateClosestToBeamLine::isValid(), reco::TrackBase::loose, reco::HitPattern::numberOfHits(), AlCaHLTBitMon_ParallelJobs::p, GeomDetEnumerators::PixelBarrel, GeomDetEnumerators::PixelEndcap, reco::HitPattern::printHitPattern(), reco::TransientTrack::setBeamSpot(), Measurement1D::significance(), reco::TransientTrack::stateAtBeamLine(), lumiQTWidget::t, funct::tan(), theB_, reco::TrackBase::tight, align::Tracker, TrajectoryStateClosestToBeamLine::transverseImpactParameter(), and vertexBeamSpot_.
{
cout << "printRecTrks" << endl;
int i =1;
for(reco::TrackCollection::const_iterator t=recTrks->begin(); t!=recTrks->end(); ++t){
// reco::TrackBase::ParameterVector par = t->parameters();
cout << endl;
cout << "Track "<<i << " " ; i++;
//enum TrackQuality { undefQuality=-1, loose=0, tight=1, highPurity=2, confirmed=3, goodIterative=4, qualitySize=5};
if( t->quality(reco::TrackBase::loose)){ cout << "loose ";};
if( t->quality(reco::TrackBase::tight)){ cout << "tight ";};
if( t->quality(reco::TrackBase::highPurity)){ cout << "highPurity ";};
if( t->quality(reco::TrackBase::confirmed)){ cout << "confirmed ";};
if( t->quality(reco::TrackBase::goodIterative)){ cout << "goodIterative ";};
cout << endl;
TransientTrack tk = theB_->build(&(*t)); tk.setBeamSpot(vertexBeamSpot_);
double ipsig=0;
if (tk.stateAtBeamLine().isValid()){
ipsig= tk.stateAtBeamLine().transverseImpactParameter().significance();
}else{
ipsig=-1;
}
cout << Form("pt=%8.3f phi=%6.3f eta=%6.3f z=%8.4f dz=%6.4f, ipsig=%6.1f",t->pt(), t->phi(), t->eta(), t->vz(), t->dzError(),ipsig) << endl;
cout << Form(" found=%6d lost=%6d chi2/ndof=%10.3f ",t->found(), t->lost(),t->normalizedChi2())<<endl;
const reco::HitPattern & p= t->hitPattern();
cout << "subdet layers valid lost" << endl;
cout << Form(" barrel %2d %2d %2d",p.pixelBarrelLayersWithMeasurement(),p.numberOfValidPixelBarrelHits(), p.numberOfLostPixelBarrelHits()) << endl;
cout << Form(" fwd %2d %2d %2d",p.pixelEndcapLayersWithMeasurement(),p.numberOfValidPixelEndcapHits(), p.numberOfLostPixelEndcapHits()) << endl;
cout << Form(" pixel %2d %2d %2d",p.pixelLayersWithMeasurement(), p.numberOfValidPixelHits(), p.numberOfLostPixelHits()) << endl;
cout << Form(" tracker %2d %2d %2d",p.trackerLayersWithMeasurement(), p.numberOfValidTrackerHits(), p.numberOfLostTrackerHits()) << endl;
cout << endl;
const reco::HitPattern & pinner= t->trackerExpectedHitsInner();
const reco::HitPattern & pouter= t->trackerExpectedHitsOuter();
int ninner=pinner.numberOfHits();
int nouter=pouter.numberOfHits();
// double d0Error=t.d0Error();
// double d0=t.dxy(myBeamSpot);
//
for(trackingRecHit_iterator hit=t->recHitsBegin(); hit!=t->recHitsEnd(); hit++){
if ((**hit).isValid() && (**hit).geographicalId().det() == DetId::Tracker ){
bool barrel = DetId((**hit).geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelBarrel);
bool endcap = DetId((**hit).geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelEndcap);
if (barrel){
const SiPixelRecHit *pixhit = dynamic_cast<const SiPixelRecHit*>( &(**hit));
edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const& clust = (*pixhit).cluster();
if (clust.isNonnull()) {
cout << Form(" barrel cluster size=%2d charge=%6.1f wx=%2d wy=%2d, expected=%3.1f",clust->size(),clust->charge(),clust->sizeX(),clust->sizeY(),1.+2./fabs(tan(t->theta()))) << endl;
}
}else if(endcap){
const SiPixelRecHit *pixhit = dynamic_cast<const SiPixelRecHit*>( &(**hit));
edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const& clust = (*pixhit).cluster();
if (clust.isNonnull()) {
cout << Form(" endcap cluster size=%2d charge=%6.1f wx=%2d wy=%2d",clust->size(),clust->charge(),clust->sizeX(),clust->sizeY()) << endl;
}
}
}
}
cout << "hitpattern" << endl;
for(int i=0; i<p.numberOfHits(); i++){ p.printHitPattern(i,std::cout); }
cout << "expected inner " << ninner << endl;
for(int i=0; i<pinner.numberOfHits(); i++){ pinner.printHitPattern(i,std::cout); }
cout << "expected outer " << nouter << endl;
for(int i=0; i<pouter.numberOfHits(); i++){ pouter.printHitPattern(i,std::cout); }
}
}
| void PrimaryVertexAnalyzer4PU::printRecVtxs | ( | const edm::Handle< reco::VertexCollection > | recVtxs, |
| std::string | title = "Reconstructed Vertices" |
||
| ) | [private] |
Definition at line 1124 of file PrimaryVertexAnalyzer4PU.cc.
References gather_cfg::cout, and v.
Referenced by analyze().
{
int ivtx=0;
std::cout << std::endl << title << std::endl;
for(reco::VertexCollection::const_iterator v=recVtxs->begin();
v!=recVtxs->end(); ++v){
string vtype=" recvtx ";
if( v->isFake()){
vtype=" fake ";
}else if (v->ndof()==-5){
vtype=" cluster "; // pos=selector[iclu],cputime[iclu],clusterz[iclu]
}else if(v->ndof()==-3){
vtype=" event ";
}
std::cout << "vtx "<< std::setw(3) << std::setfill(' ')<<ivtx++
<< vtype
<< " #trk " << std::fixed << std::setprecision(4) << std::setw(3) << v->tracksSize()
<< " chi2 " << std::fixed << std::setw(4) << std::setprecision(1) << v->chi2()
<< " ndof " << std::fixed << std::setw(5) << std::setprecision(2) << v->ndof()
<< " 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 PrimaryVertexAnalyzer4PU::printSimTrks | ( | const edm::Handle< edm::SimTrackContainer > | simVtrks | ) | [private] |
Definition at line 1176 of file PrimaryVertexAnalyzer4PU.cc.
References gather_cfg::cout, i, and lumiQTWidget::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->eventId().event() << "," << t->eventId().bunchCrossing()
<< (*t)
<< " index="
<< (*t).genpartIndex();
//if (gp) {
// HepMC::GenVertex *gv=gp->production_vertex();
// std::cout << " genvertex =" << (*gv);
//}
std::cout << std::endl;
}
}
| void PrimaryVertexAnalyzer4PU::printSimVtxs | ( | const edm::Handle< edm::SimVertexContainer > | simVtxs | ) | [private] |
Definition at line 1153 of file PrimaryVertexAnalyzer4PU.cc.
References gather_cfg::cout, i, and simUnit_.
{
int i=0;
for(SimVertexContainer::const_iterator vsim=simVtxs->begin();
vsim!=simVtxs->end(); ++vsim){
if ( vsim->position().x()*vsim->position().x()+vsim->position().y()*vsim->position().y() < 1.){
std::cout << i++ << ")" << std::scientific
<< " evtid=" << vsim->eventId().event() << "," << vsim->eventId().bunchCrossing()
<< " 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;
}
}
}
| int * PrimaryVertexAnalyzer4PU::supf | ( | std::vector< SimPart > & | simtrks, |
| const reco::TrackCollection & | trks | ||
| ) | [private] |
Definition at line 831 of file PrimaryVertexAnalyzer4PU.cc.
References trackerHits::c, gather_cfg::cout, create_public_lumi_plots::exp, i, j, gen::k, prof2calltree::l, create_public_lumi_plots::log, M_PI, match(), RPCpg::mu, PrimaryVertexAnalyzer4PU::SimPart::par, alignCSCRings::s, lumiQTWidget::t, and funct::tan().
Referenced by printPVTrks().
{
int nsim=simtrks.size();
int nrec=trks.size();
int *rectosim=new int[nrec]; // pointer to associated simtrk
double** pij=new double*[nrec];
double mu=100.; // initial chi^2 cut-off (5 dofs !)
int nmatch=0;
int i=0;
for(reco::TrackCollection::const_iterator t=trks.begin(); t!=trks.end(); ++t){
pij[i]=new double[nsim];
rectosim[i]=-1;
ParameterVector par = t->parameters();
//reco::TrackBase::CovarianceMatrix V = t->covariance();
reco::TrackBase::CovarianceMatrix S = t->covariance();
S.Invert();
for(int j=0; j<nsim; j++){
simtrks[j].rec=-1;
SimPart s=simtrks[j];
double c=0;
for(int k=0; k<5; k++){
for(int l=0; l<5; l++){
c+=(par(k)-s.par[k])*(par(l)-s.par[l])*S(k,l);
}
}
pij[i][j]=exp(-0.5*c);
// double c0=pow((par[0]-s.par[0])/t->qoverpError(),2)*0.1
// +pow((par[1]-s.par[1])/t->lambdaError(),2)
// +pow((par[2]-s.par[2])/t->phiError(),2)
// +pow((par[3]-s.par[3])/t->dxyError(),2)*0.1;
// +pow((par[4]-s.par[4])/t->dszError(),2)*0.1;
// pij[i][j]=exp(-0.5*c0);
// if( c0 <100 ){
// cout << setw(3) << i << " rec " << setw(6) << par << endl;
// cout << setw(3) << j << " sim " << setw(6) << s.par << " ---> C=" << c << endl;
// cout << " " << setw(6)
// << (par[0]-s.par[0])<< ","
// << (par[1]-s.par[1])<< ","
// << (par[2]-s.par[2])<< ","
// << (par[3]-s.par[3])<< ","
// << (par[4]-s.par[4])
// << " match=" << match(simtrks[j].par, trks.at(i).parameters())
// << endl;
// cout << " " << setw(6)
// << (par[0]-s.par[0])/t->qoverpError() << ","
// << (par[1]-s.par[1])/t->lambdaError() << ","
// << (par[2]-s.par[2])/t->phiError() << ","
// << (par[3]-s.par[3])/t->dxyError() << ","
// << (par[4]-s.par[4])/t->dszError() << " c0=" << c0
// << endl <<endl;
// }
}
i++;
}
for(int k=0; k<nrec; k++){
int imatch=-1; int jmatch=-1;
double pmatch=0;
for(int j=0; j<nsim; j++){
if ((simtrks[j].rec)<0){
double psum=exp(-0.5*mu); //cutoff
for(int i=0; i<nrec; i++){
if (rectosim[i]<0){ psum+=pij[i][j];}
}
for(int i=0; i<nrec; i++){
if ((rectosim[i]<0)&&(pij[i][j]/psum>pmatch)){
pmatch=pij[i][j]/psum;
imatch=i; jmatch=j;
}
}
}
}// sim loop
if((jmatch>=0)||(imatch>=0)){
//std::cout << pmatch << " " << pij[imatch][jmatch] << " match=" <<
// match(simtrks[jmatch].par, trks.at(imatch).parameters()) <<std::endl;
if (pmatch>0.01){
rectosim[imatch]=jmatch;
simtrks[jmatch].rec=imatch;
nmatch++;
}else if (match(simtrks[jmatch].par, trks.at(imatch).parameters())){
// accept it anyway if it matches crudely and relax the cut-off
rectosim[imatch]=jmatch;
simtrks[jmatch].rec=imatch;
nmatch++;
mu=mu*2;
}
}
}
// std::cout << ">>>>>>>>>>>>>>>--------------supf----------------------" << std::endl;
// std::cout <<"nsim=" << nsim << " nrec=" << nrec << " nmatch=" << nmatch << std::endl;
// std::cout << "rec to sim " << std::endl;
// for(int i=0; i<nrec; i++){
// std::cout << i << " ---> " << rectosim[i] << std::endl;
// }
// std::cout << "sim to rec " << std::endl;
// for(int j=0; j<nsim; j++){
// std::cout << j << " ---> " << simtrks[j].rec << std::endl;
// }
std::cout << "unmatched sim " << std::endl;
for(int j=0; j<nsim; j++){
if(simtrks[j].rec<0){
double pt= 1./simtrks[j].par[0]/tan(simtrks[j].par[1]);
if((fabs(pt))>1.){
std::cout << setw(3) << j << setw(8) << simtrks[j].pdg
<< setw(8) << setprecision(4) << " (" << simtrks[j].xvtx << "," << simtrks[j].yvtx << "," << simtrks[j].zvtx << ")"
<< " pt= " << pt
<< " phi=" << simtrks[j].par[2]
<< " eta= " << -log(tan(0.5*(M_PI/2-simtrks[j].par[1])))
<< std::endl;
}
}
}
// std::cout << "<<<<<<<<<<<<<<<--------------supf----------------------" << std::endl;
//delete rectosim; // or return it?
for(int i=0; i<nrec; i++){delete pij[i];}
delete pij;
return rectosim; // caller must delete it
}
| bool PrimaryVertexAnalyzer4PU::truthMatchedTrack | ( | edm::RefToBase< reco::Track > | track, |
| TrackingParticleRef & | tpr | ||
| ) | [private] |
Definition at line 1517 of file PrimaryVertexAnalyzer4PU.cc.
References event(), f, and r2s_.
Referenced by getTruthMatchedVertexTracks().
{
double f=0;
try{
std::vector<std::pair<TrackingParticleRef, double> > tp = r2s_[track];
for (std::vector<std::pair<TrackingParticleRef, double> >::const_iterator it = tp.begin();
it != tp.end(); ++it) {
if (it->second>f){
tpr=it->first;
f=it->second;
}
}
} catch (Exception event) {
// silly way of testing whether track is in r2s_
}
// sanity check on track parameters?
return (f>0.5);
}
| std::string PrimaryVertexAnalyzer4PU::vertexString | ( | TrackingParticleRefVector | , |
| TrackingParticleRefVector | |||
| ) | const [private] |
| std::string PrimaryVertexAnalyzer4PU::vertexString | ( | HepMC::GenVertex::particles_in_const_iterator | , |
| HepMC::GenVertex::particles_in_const_iterator | , | ||
| HepMC::GenVertex::particles_out_const_iterator | , | ||
| HepMC::GenVertex::particles_out_const_iterator | |||
| ) | const [private] |
Definition at line 355 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyze().
int PrimaryVertexAnalyzer4PU::bunchCrossing_ [private] |
Definition at line 342 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyze(), and analyzeVertexCollection().
bool PrimaryVertexAnalyzer4PU::DEBUG_ [private] |
Definition at line 345 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyzeVertexCollection(), getSimPVs(), and PrimaryVertexAnalyzer4PU().
bool PrimaryVertexAnalyzer4PU::doMatching_ [private] |
Definition at line 324 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyze(), and PrimaryVertexAnalyzer4PU().
int PrimaryVertexAnalyzer4PU::dumpcounter_ [private] |
Definition at line 333 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyze(), and PrimaryVertexAnalyzer4PU().
bool PrimaryVertexAnalyzer4PU::dumpPUcandidates_ [private] |
Definition at line 336 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyzeVertexCollection(), and PrimaryVertexAnalyzer4PU().
bool PrimaryVertexAnalyzer4PU::dumpThisEvent_ [private] |
Definition at line 335 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyze(), and analyzeVertexCollection().
int PrimaryVertexAnalyzer4PU::event_ [private] |
Definition at line 341 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyze(), and analyzeVertexCollection().
int PrimaryVertexAnalyzer4PU::eventcounter_ [private] |
Definition at line 332 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyze(), analyzeVertexCollection(), analyzeVertexCollectionTP(), PrimaryVertexAnalyzer4PU(), and printEventSummary().
double PrimaryVertexAnalyzer4PU::fBfield_ [private] |
Definition at line 321 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyze(), fillTrackHistos(), and getSimTrkParameters().
std::map<std::string, TH1*> PrimaryVertexAnalyzer4PU::hBS [private] |
Definition at line 348 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyze(), beginJob(), and endJob().
std::map<std::string, TH1*> PrimaryVertexAnalyzer4PU::hDA [private] |
Definition at line 350 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyze(), beginJob(), and endJob().
std::map<std::string, TH1*> PrimaryVertexAnalyzer4PU::hMVF [private] |
Definition at line 352 of file PrimaryVertexAnalyzer4PU.h.
std::map<std::string, TH1*> PrimaryVertexAnalyzer4PU::hnoBS [private] |
Definition at line 349 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyze(), beginJob(), and endJob().
std::map<std::string, TH1*> PrimaryVertexAnalyzer4PU::hPIX [private] |
Definition at line 351 of file PrimaryVertexAnalyzer4PU.h.
std::map<std::string, TH1*> PrimaryVertexAnalyzer4PU::hsimPV [private] |
Definition at line 353 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyze(), beginJob(), endJob(), and getSimPVs().
int PrimaryVertexAnalyzer4PU::luminosityBlock_ [private] |
Definition at line 340 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyze().
Definition at line 330 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyze().
int PrimaryVertexAnalyzer4PU::ndump_ [private] |
Definition at line 334 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyze(), and PrimaryVertexAnalyzer4PU().
int PrimaryVertexAnalyzer4PU::orbitNumber_ [private] |
Definition at line 343 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyze().
std::string PrimaryVertexAnalyzer4PU::outputFile_ [private] |
Definition at line 319 of file PrimaryVertexAnalyzer4PU.h.
Referenced by PrimaryVertexAnalyzer4PU().
edm::ESHandle< ParticleDataTable > PrimaryVertexAnalyzer4PU::pdt_ [private] |
Definition at line 329 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyze(), isCharged(), and isResonance().
bool PrimaryVertexAnalyzer4PU::printXBS_ [private] |
Definition at line 325 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyze(), and PrimaryVertexAnalyzer4PU().
Definition at line 356 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyze(), and truthMatchedTrack().
Definition at line 362 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyze().
std::string PrimaryVertexAnalyzer4PU::recoTrackProducer_ [private] |
Definition at line 318 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyze(), and PrimaryVertexAnalyzer4PU().
TFile* PrimaryVertexAnalyzer4PU::rootFile_ [private] |
Definition at line 322 of file PrimaryVertexAnalyzer4PU.h.
Referenced by beginJob(), endJob(), PrimaryVertexAnalyzer4PU(), and ~PrimaryVertexAnalyzer4PU().
int PrimaryVertexAnalyzer4PU::run_ [private] |
Definition at line 339 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyze(), and analyzeVertexCollection().
Definition at line 326 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyze(), and PrimaryVertexAnalyzer4PU().
double PrimaryVertexAnalyzer4PU::simUnit_ [private] |
Definition at line 327 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyze(), analyzeVertexCollection(), beginJob(), matchVertex(), PrimaryVertexAnalyzer4PU(), and printSimVtxs().
Definition at line 363 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyze(), analyzeVertexCollection(), fillTrackHistos(), printPVTrks(), and printRecTrks().
Definition at line 359 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyzeVertexCollection(), printEventSummary(), and printPVTrks().
bool PrimaryVertexAnalyzer4PU::verbose_ [private] |
Definition at line 323 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyze(), analyzeVertexCollection(), getSimPVs(), getSimTrkParameters(), and PrimaryVertexAnalyzer4PU().
Definition at line 360 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyze(), analyzeVertexCollection(), fillTrackHistos(), printPVTrks(), and printRecTrks().
std::vector<std::string> PrimaryVertexAnalyzer4PU::vtxSample_ [private] |
Definition at line 320 of file PrimaryVertexAnalyzer4PU.h.
double PrimaryVertexAnalyzer4PU::wxy2_ [private] |
Definition at line 361 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyze(), and fillTrackHistos().
std::map<double, TrackingParticleRef> PrimaryVertexAnalyzer4PU::z2tp_ [private] |
Definition at line 357 of file PrimaryVertexAnalyzer4PU.h.
Referenced by printEventSummary(), and printPVTrks().
double PrimaryVertexAnalyzer4PU::zmatch_ [private] |
Definition at line 328 of file PrimaryVertexAnalyzer4PU.h.
Referenced by analyzeVertexCollection(), analyzeVertexCollectionTP(), matchVertex(), and PrimaryVertexAnalyzer4PU().
1.7.3