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

Inheritance diagram for PrimaryVertexAnalyzer4PU:

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 ()

Public Member Functions inherited from edm::EDAnalyzer
	EDAnalyzer ()

std::string	workerType () const

virtual	~EDAnalyzer ()

Public Member Functions inherited from edm::EDConsumerBase
	EDConsumerBase ()

ProductHolderIndex	indexFrom (EDGetToken, BranchType, TypeID const &) const

void	itemsMayGet (BranchType, std::vector< ProductHolderIndex > &) const

void	itemsToGet (BranchType, std::vector< ProductHolderIndex > &) const

void	labelsForToken (EDGetToken iToken, Labels &oLabels) const

void	updateLookup (BranchType iBranchType, ProductHolderIndexHelper const &)

virtual	~EDConsumerBase ()

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, double x)

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, bool signal)

void	Fill (std::map< std::string, TH1 * > &h, std::string s, bool yesno, bool signal)

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 (TrackingParticleRefVector, TrackingParticleRefVector) const

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

Static Private Member Functions
static bool	match (const ParameterVector &a, const ParameterVector &b)

Private Attributes
TrackAssociatorBase *	associatorByHits_

edm::InputTag	beamSpot_

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_

Additional Inherited Members
Public Types inherited from edm::EDAnalyzer
typedef EDAnalyzer	ModuleType

typedef WorkerT< EDAnalyzer >	WorkerType

Static Public Member Functions inherited from edm::EDAnalyzer
static const std::string &	baseType ()

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &)

Protected Member Functions inherited from edm::EDAnalyzer
void	callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)

CurrentProcessingContext const *	currentContext () const

Protected Member Functions inherited from edm::EDConsumerBase
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	consumes (edm::InputTag const &tag)

EDGetToken	consumes (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	consumes (TypeToGet const &id, edm::InputTag const &tag)

ConsumesCollector	consumesCollector ()
	Use a ConsumesCollector to gather consumes information from helper functions. More...

template<typename ProductType , BranchType B = InEvent>
void	consumesMany ()

void	consumesMany (const TypeToGet &id)

template<BranchType B>
void	consumesMany (const TypeToGet &id)

template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	mayConsume (edm::InputTag const &tag)

EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

Detailed Description

Description: primary vertex analyzer for events with pile-up

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

Definition at line 76 of file PrimaryVertexAnalyzer4PU.h.

Member Typedef Documentation

typedef math::XYZTLorentzVector PrimaryVertexAnalyzer4PU::LorentzVector

private

Definition at line 78 of file PrimaryVertexAnalyzer4PU.h.

typedef reco::TrackBase::ParameterVector PrimaryVertexAnalyzer4PU::ParameterVector

private

Definition at line 80 of file PrimaryVertexAnalyzer4PU.h.

Constructor & Destructor Documentation

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_, AlCaHLTBitMon_QueryRunRegistry::string, verbose_, and zmatch_.

                                                                               :
   theTrackFilter(iConfig.getParameter<edm::ParameterSet>("TkFilterParameters")),
   beamSpot_(iConfig.getParameter<edm::InputTag>("beamSpot"))
 {
    //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 107 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_;
 }

Member Function Documentation

void PrimaryVertexAnalyzer4PU::add	(	std::map< std::string, TH1 * > &	h,
		TH1 *	hist
	)

inlineprivate

Definition at line 185 of file PrimaryVertexAnalyzer4PU.h.

References estimatePileup::hist.

Referenced by beginJob(), and bookVertexHistograms().

185 { h[hist->GetName()]=hist; hist->StatOverflows(kTRUE);}

estimatePileup.hist

tuple hist

Definition: estimatePileup.py:215

h

The Signals That Services Can Subscribe To This is based on ActivityRegistry h

Helper function to determine trigger accepts.

Definition: Activities.doc:4

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

virtual

Implements edm::EDAnalyzer.

Definition at line 1914 of file PrimaryVertexAnalyzer4PU.cc.

References analyzeVertexCollection(), analyzeVertexCollectionTP(), TrackAssociatorBase::associateRecoToSim(), associatorByHits_, align::BeamSpot, 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::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::Handle< T >::product(), edm::ESHandle< class >::product(), r2s_, recoBeamSpotHandle_, recoTrackProducer_, edm::EventID::run(), run_, simG4_, simUnit_, AlCaHLTBitMon_QueryRunRegistry::string, lumiQTWidget::t, theB_, findQualityFiles::v, verbose_, vertexBeamSpot_, wxy2_, reco::BeamSpot::x0(), reco::BeamSpot::y0(), detailsBasic3DVector::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.getByLabel(beamSpot_, 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("mix","MergedTrackTruth",TPCollectionH);
   bool gotTV=iEvent.getByLabel("mix","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, &iSetup ); 
     //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 2760 of file PrimaryVertexAnalyzer4PU.cc.

References reco::BeamSpot::BeamWidthX(), reco::BeamSpot::BeamWidthY(), bunchCrossing_, reco::TrackBase::charge(), ChiSquaredProbability(), gather_cfg::cout, data, DEBUG_, dumpPUcandidates_, dumpThisEvent_, reco::TrackBase::eta(), event_, eventcounter_, Fill(), fillTrackHistos(), i, getHLTprescales::index, edm::detail::isnan(), edm::HandleBase::isValid(), j, create_public_lumi_plots::log, m, M_PI, nt, NULL, reco::TrackBase::p(), reco::TrackBase::phi(), position, funct::pow(), reco::TrackBase::px(), reco::TrackBase::py(), reco::TrackBase::pz(), run_, reco::TransientTrack::setBeamSpot(), simUnit_, mathSSE::sqrt(), lumiQTWidget::t, reco::t2, theB_, theTrackFilter, reco::Vertex::tracksSize(), groupFilesInBlocks::tt, findQualityFiles::v, verbose_, vertexBeamSpot_, w(), x, reco::BeamSpot::x0(), reco::BeamSpot::y0(), detailsBasic3DVector::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 2449 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(), findQualityFiles::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 566 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 121 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 )

inlineprivate

Definition at line 235 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 1104 of file PrimaryVertexAnalyzer4PU.cc.

References Reference_intrackfit_cff::barrel, gather_cfg::cout, reco::TrackBase::eta(), edm::Ref< C, T, F >::isNonnull(), PixelSubdetector::PixelBarrel, reco::TrackBase::pt(), reco::Track::recHitsBegin(), reco::Track::recHitsEnd(), and DetId::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 661 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
	)

inlineprivate

Definition at line 187 of file PrimaryVertexAnalyzer4PU.h.

References gather_cfg::cout, and alignCSCRings::s.

Referenced by analyze(), analyzeVertexCollection(), analyzeVertexCollectionTP(), Fill(), fillTrackHistos(), and printEventSummary().

                                                               {
     //    cout << "Fill1 " << s << endl;
     if(h.count(s)==0){
       std::cout << "Trying to fill non-exiting Histogram named " << s << std::endl;
       return;
     }
     h[s]->Fill(x);
   }

void PrimaryVertexAnalyzer4PU::Fill	(	std::map< std::string, TH1 * > &	h,
		std::string	s,
		double	x,
		double	y
	)

inlineprivate

Definition at line 196 of file PrimaryVertexAnalyzer4PU.h.

References gather_cfg::cout, and alignCSCRings::s.

                                                                         {
     //    cout << "Fill2 " << s << endl;
     if(h.count(s)==0){
       std::cout << "Trying to fill non-exiting Histogram named " << s << std::endl;
       return;
     }
     h[s]->Fill(x,y);
   }

void PrimaryVertexAnalyzer4PU::Fill	(	std::map< std::string, TH1 * > &	h,
		std::string	s,
		double	x,
		bool	signal
	)

inlineprivate

Definition at line 205 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
	)

inlineprivate

Definition at line 227 of file PrimaryVertexAnalyzer4PU.h.

References Fill().

                                                                              {
     if (yesno){
       Fill(h, s,1.,signal);
     }else{
       Fill(h, s,0.,signal);
     }
   }

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 1005 of file PrimaryVertexAnalyzer4PU.cc.

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 1580 of file PrimaryVertexAnalyzer4PU.cc.

References abs, EncodedEventId::bunchCrossing(), gather_cfg::cout, reco::TrackBase::dxy(), reco::TrackBase::dxyError(), alignCSCRings::e, EncodedEventId::event(), event(), TrackingVertex::eventId(), PrimaryVertexAnalyzer4PU::SimEvent::eventId, edm::RefToBase< T >::get(), edm::Ref< C, T, F >::get(), i, TransientVertex::isValid(), TrackingVertex::position(), TransientVertex::position(), funct::pow(), edm::Handle< T >::product(), reco::TransientTrack::setBeamSpot(), edm::View< T >::size(), mathSSE::sqrt(), lumiQTWidget::t, findQualityFiles::v, AdaptiveVertexFitter::vertex(), reco::TrackBase::vz(), x, PV3DBase< T, PVType, FrameType >::x(), PrimaryVertexAnalyzer4PU::SimEvent::x, detailsBasic3DVector::y, PV3DBase< T, PVType, FrameType >::y(), PrimaryVertexAnalyzer4PU::SimEvent::y, detailsBasic3DVector::z, PV3DBase< T, PVType, FrameType >::z(), and PrimaryVertexAnalyzer4PU::SimEvent::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 1712 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 1833 of file PrimaryVertexAnalyzer4PU.cc.

References begin, gather_cfg::cout, DEBUG_, alignCSCRings::e, end, event(), PrimaryVertexAnalyzer4PU::simPrimaryVertex::eventId, NULL, benchmark_cfg::pdgId, position, and findQualityFiles::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 735 of file PrimaryVertexAnalyzer4PU.cc.

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 > &	tracks,
		double &	Tc,
		double &	chsq,
		double &	dzmax,
		double &	dztrim,
		double &	m4m2
	)

private

Definition at line 1468 of file PrimaryVertexAnalyzer4PU.cc.

References a, b, SiPixelRawToDigiRegional_cfi::beamSpot, position, funct::pow(), mathSSE::sqrt(), funct::tan(), w(), detailsBasic3DVector::z, SiStripMonitorClusterAlca_cfi::zmax, and SiStripMonitorClusterAlca_cfi::zmin.

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 1553 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 992 of file PrimaryVertexAnalyzer4PU.cc.

References pdt_.

Referenced by getSimPVs().

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

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

private

Definition at line 987 of file PrimaryVertexAnalyzer4PU.cc.

Referenced by getSimPVs().

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

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

private

Definition at line 981 of file PrimaryVertexAnalyzer4PU.cc.

References abs, alignCSCRings::e, and pdt_.

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

bool PrimaryVertexAnalyzer4PU::match	(	const ParameterVector &	a,
		const ParameterVector &	b
	)

staticprivate

Definition at line 964 of file PrimaryVertexAnalyzer4PU.cc.

References a, b, and M_PI.

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 1420 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 976 of file PrimaryVertexAnalyzer4PU.cc.

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

                                                     {
   return (fabs(vsim.z*simUnit_-vrec.z())<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 2224 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(), customizeTrackingMonitorSeedNumber::idx, 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(), findQualityFiles::v, reco::TrackBase::vz(), detailsBasic3DVector::z, PrimaryVertexAnalyzer4PU::SimEvent::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 1285 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, listHistos::IP, TrackingVertex::position(), position, funct::pow(), reco::TransientTrack::setBeamSpot(), mathSSE::sqrt(), ntuplemaker::status, supf(), lumiQTWidget::t, funct::tan(), theB_, theta(), theTrackFilter, groupFilesInBlocks::tt, findQualityFiles::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 1199 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, PixelSubdetector::PixelBarrel, PixelSubdetector::PixelEndcap, reco::HitPattern::printHitPattern(), reco::TransientTrack::setBeamSpot(), Measurement1D::significance(), reco::TransientTrack::stateAtBeamLine(), lumiQTWidget::t, funct::tan(), theB_, reco::TrackBase::tight, DetId::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 1126 of file PrimaryVertexAnalyzer4PU.cc.

References gather_cfg::cout, and findQualityFiles::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 1178 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 1155 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 833 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 1519 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