#include <Alignment/OfflineValidation/plugins/PrimaryVertexValidation.cc>

Inheritance diagram for PrimaryVertexValidation:

Public Member Functions
	PrimaryVertexValidation (const edm::ParameterSet &)

	~PrimaryVertexValidation () override

Public Member Functions inherited from edm::one::EDAnalyzer< edm::one::SharedResources >
	EDAnalyzer ()=default

SerialTaskQueue *	globalLuminosityBlocksQueue () final

SerialTaskQueue *	globalRunsQueue () final

bool	wantsGlobalLuminosityBlocks () const final

bool	wantsGlobalRuns () const final

Public Member Functions inherited from edm::one::EDAnalyzerBase
void	callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)

	EDAnalyzerBase ()

ModuleDescription const &	moduleDescription () const

bool	wantsStreamLuminosityBlocks () const

bool	wantsStreamRuns () const

	~EDAnalyzerBase () override

Public Member Functions inherited from edm::EDConsumerBase
std::vector< ConsumesInfo >	consumesInfo () const

void	convertCurrentProcessAlias (std::string const &processName)
	Convert "@currentProcess" in InputTag process names to the actual current process name. More...

	EDConsumerBase ()

	EDConsumerBase (EDConsumerBase const &)=delete

	EDConsumerBase (EDConsumerBase &&)=default

ESProxyIndex const *	esGetTokenIndices (edm::Transition iTrans) const

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

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

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

std::vector< ProductResolverIndexAndSkipBit > const &	itemsToGetFrom (BranchType iType) const

void	labelsForToken (EDGetToken iToken, Labels &oLabels) const

void	modulesWhoseProductsAreConsumed (std::vector< ModuleDescription const * > &modules, ProductRegistry const &preg, std::map< std::string, ModuleDescription const * > const &labelsToDesc, std::string const &processName) const

EDConsumerBase const &	operator= (EDConsumerBase const &)=delete

EDConsumerBase &	operator= (EDConsumerBase &&)=default

bool	registeredToConsume (ProductResolverIndex, bool, BranchType) const

bool	registeredToConsumeMany (TypeID const &, BranchType) const

ProductResolverIndexAndSkipBit	uncheckedIndexFrom (EDGetToken) const

void	updateLookup (BranchType iBranchType, ProductResolverIndexHelper const &, bool iPrefetchMayGet)

void	updateLookup (eventsetup::ESRecordsToProxyIndices const &)

virtual	~EDConsumerBase () noexcept(false)

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

void	analyze (const edm::Event &, const edm::EventSetup &) override

void	beginJob () override

std::vector< TH1F * >	bookResidualsHistogram (const TFileDirectory &dir, unsigned int theNOfBins, PVValHelper::residualType resType, PVValHelper::plotVariable varType, bool isNormalized=false)

std::map< std::string, TH1 * >	bookVertexHistograms (const TFileDirectory &dir)

void	endJob () override

void	fill (std::map< std::string, TH1 * > &h, const std::string &s, double x)

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

void	fillByIndex (std::vector< TH1F * > &h, unsigned int index, double x, std::string tag="")

void	fillMap (TH2F trendMap, TH1F residualsMapPlot[100][100], PVValHelper::estimator fitPar_)

void	fillTrackHistos (std::map< std::string, TH1 * > &h, const std::string &ttype, const reco::TransientTrack *tt, const reco::Vertex &v, const reco::BeamSpot &beamSpot, double fBfield)

void	fillTrendPlot (TH1F trendPlot, TH1F residualsPlot[100], PVValHelper::estimator fitPar_, const std::string &var_)

void	fillTrendPlotByIndex (TH1F trendPlot, std::vector< TH1F > &h, PVValHelper::estimator fitPar_, PVValHelper::plotVariable plotVar=PVValHelper::END_OF_PLOTS)

std::pair< Measurement1D, Measurement1D >	fitResiduals (TH1 *hist)

Measurement1D	getMAD (TH1F *histo)

Measurement1D	getMedian (TH1F *histo)

std::tuple< std::string, std::string, std::string >	getTypeString (PVValHelper::residualType type)

std::tuple< std::string, std::string, std::string >	getVarString (PVValHelper::plotVariable var)

bool	hasFirstLayerPixelHits (const reco::TransientTrack &track)

bool	isBFieldConsistentWithMode (const edm::EventSetup &iSetup) const

bool	isHit2D (const TrackingRecHit &hit) const

bool	passesTrackCuts (const reco::Track &track, const reco::Vertex &vertex, const std::string &qualityString_, double dxyErrMax_, double dzErrMax_, double ptErrMax_)

std::pair< bool, bool >	pixelHitsCheck (const reco::TransientTrack &track)

void	SetVarToZero ()

void	shrinkHistVectorToFit (std::vector< TH1F * > &h, unsigned int desired_size)

double	square (double x)

Static Private Member Functions
static bool	vtxSort (const reco::Vertex &a, const reco::Vertex &b)

Private Attributes
std::vector< TH1F * >	a_d3DEtaResiduals

std::vector< TH1F * >	a_d3DPhiResiduals

TH1F *	a_d3DResidualsMap [nMaxBins_][nMaxBins_]

std::vector< TH1F * >	a_dxEtaResiduals

std::vector< TH1F * >	a_dxPhiResiduals

TH1F *	a_dxyBiasResidualsMap [nMaxBins_][nMaxBins_]

std::vector< TH1F * >	a_dxyEtaBiasResiduals

TH1F *	a_dxyEtaMADBiasTrend

TH1F *	a_dxyEtaMADTrend

TH1F *	a_dxyEtaMeanBiasTrend

TH1F *	a_dxyEtaMeanTrend

TH1F *	a_dxyEtaMedianBiasTrend

TH1F *	a_dxyEtaMedianTrend

std::vector< TH1F * >	a_dxyEtaResiduals

TH1F *	a_dxyEtaWidthBiasTrend

TH1F *	a_dxyEtaWidthTrend

TH1F *	a_dxyladderMeanTrend

TH1F *	a_dxyladderWidthTrend

TH2F *	a_dxyMeanBiasMap

TH2F *	a_dxyMeanMap

TH1F *	a_dxymodZMeanTrend

TH1F *	a_dxymodZWidthTrend

std::vector< TH1F * >	a_dxyPhiBiasResiduals

TH1F *	a_dxyPhiMADBiasTrend

TH1F *	a_dxyPhiMADTrend

TH1F *	a_dxyPhiMeanBiasTrend

TH1F *	a_dxyPhiMeanTrend

TH1F *	a_dxyPhiMedianBiasTrend

TH1F *	a_dxyPhiMedianTrend

std::vector< TH1F * >	a_dxyPhiResiduals

TH1F *	a_dxyPhiWidthBiasTrend

TH1F *	a_dxyPhiWidthTrend

TH1F *	a_dxypTCentralMeanTrend

TH1F *	a_dxypTCentralWidthTrend

TH1F *	a_dxypTMeanTrend

TH1F *	a_dxypTWidthTrend

TH1F *	a_dxyResidualsMap [nMaxBins_][nMaxBins_]

TH2F *	a_dxyVsEta

TH2F *	a_dxyVsPhi

TH2F *	a_dxyWidthBiasMap

TH2F *	a_dxyWidthMap

std::vector< TH1F * >	a_dyEtaResiduals

std::vector< TH1F * >	a_dyPhiResiduals

TH1F *	a_dzBiasResidualsMap [nMaxBins_][nMaxBins_]

std::vector< TH1F * >	a_dzEtaBiasResiduals

TH1F *	a_dzEtaMADBiasTrend

TH1F *	a_dzEtaMADTrend

TH1F *	a_dzEtaMeanBiasTrend

TH1F *	a_dzEtaMeanTrend

TH1F *	a_dzEtaMedianBiasTrend

TH1F *	a_dzEtaMedianTrend

std::vector< TH1F * >	a_dzEtaResiduals

TH1F *	a_dzEtaWidthBiasTrend

TH1F *	a_dzEtaWidthTrend

TH1F *	a_dzladderMeanTrend

TH1F *	a_dzladderWidthTrend

TH2F *	a_dzMeanBiasMap

TH2F *	a_dzMeanMap

TH1F *	a_dzmodZMeanTrend

TH1F *	a_dzmodZWidthTrend

std::vector< TH1F * >	a_dzPhiBiasResiduals

TH1F *	a_dzPhiMADBiasTrend

TH1F *	a_dzPhiMADTrend

TH1F *	a_dzPhiMeanBiasTrend

TH1F *	a_dzPhiMeanTrend

TH1F *	a_dzPhiMedianBiasTrend

TH1F *	a_dzPhiMedianTrend

std::vector< TH1F * >	a_dzPhiResiduals

TH1F *	a_dzPhiWidthBiasTrend

TH1F *	a_dzPhiWidthTrend

TH1F *	a_dzpTCentralMeanTrend

TH1F *	a_dzpTCentralWidthTrend

TH1F *	a_dzpTMeanTrend

TH1F *	a_dzpTWidthTrend

TH1F *	a_dzResidualsMap [nMaxBins_][nMaxBins_]

TH2F *	a_dzVsEta

TH2F *	a_dzVsPhi

TH2F *	a_dzWidthBiasMap

TH2F *	a_dzWidthMap

std::vector< TH1F * >	a_IP2DEtaResiduals

std::vector< TH1F * >	a_IP2DPhiResiduals

std::vector< TH1F * >	a_IP3DEtaResiduals

std::vector< TH1F * >	a_IP3DPhiResiduals

std::vector< TH1F * >	a_reszEtaResiduals

std::vector< TH1F * >	a_reszPhiResiduals

bool	askFirstLayerHit_

double	Beamdxdz_

double	Beamsigmaz_

double	BeamWidthX_

double	BeamWidthY_

double	BSx0_

double	BSy0_

double	BSz0_

int	charge_ [nMaxtracks_]

double	chi2_ [nMaxtracks_]

double	chi2ndof_ [nMaxtracks_]

float	chi2normUnbiasedVertex_ [nMaxtracks_]

float	chi2ProbUnbiasedVertex_ [nMaxtracks_]

float	chi2UnbiasedVertex_ [nMaxtracks_]

double	d3DErrorFromMyVertex_ [nMaxtracks_]

double	d3DFromMyVertex_ [nMaxtracks_]

bool	debug_

bool	doBPix_

bool	doFPix_

float	DOFUnbiasedVertex_ [nMaxtracks_]

double	dxy_ [nMaxtracks_]

double	dxyBs_ [nMaxtracks_]

double	dxyErrorFromMyVertex_ [nMaxtracks_]

double	dxyFromMyVertex_ [nMaxtracks_]

double	dz_ [nMaxtracks_]

double	dzBs_ [nMaxtracks_]

double	dzErrorFromMyVertex_ [nMaxtracks_]

double	dzFromMyVertex_ [nMaxtracks_]

double	eta_ [nMaxtracks_]

double	etaOfProbe_

unsigned int	EventNumber_

bool	forceBeamSpotContraint_

edm::Service< TFileService >	fs

TH2F *	h2_probeEtaPhi_

TH2F *	h2_probeEtaPt_

TH1F *	h_Beamsigmaz

TH1F *	h_BeamWidthX

TH1F *	h_BeamWidthY

TH1F *	h_BSx0

TH1F *	h_BSy0

TH1F *	h_BSz0

std::vector< TH1F * >	h_dxy_Central_pT_

std::vector< TH1F * >	h_dxy_ladder_

std::vector< TH1F * >	h_dxy_ladderNoOverlap_

std::vector< TH1F * >	h_dxy_ladderOverlap_

std::vector< TH1F * >	h_dxy_modZ_

std::vector< TH1F * >	h_dxy_pT_

std::vector< TH1F * >	h_dz_Central_pT_

std::vector< TH1F * >	h_dz_ladder_

std::vector< TH1F * >	h_dz_modZ_

std::vector< TH1F * >	h_dz_pT_

TH1F *	h_etaMax

TH1F *	h_fitVtxChi2_

TH1F *	h_fitVtxChi2ndf_

TH1F *	h_fitVtxChi2Prob_

TH1F *	h_fitVtxNdof_

TH1F *	h_fitVtxNtracks_

TH1F *	h_fitVtxTrackAverageWeight_

TH1F *	h_fitVtxTrackWeights_

TH1F *	h_lumiFromConfig

TH1F *	h_nbins

TH1F *	h_nClus

TH1F *	h_nLadders

TH1F *	h_nOfflineVertices

std::vector< TH1F * >	h_norm_dxy_Central_pT_

std::vector< TH1F * >	h_norm_dxy_ladder_

std::vector< TH1F * >	h_norm_dxy_modZ_

std::vector< TH1F * >	h_norm_dxy_pT_

std::vector< TH1F * >	h_norm_dz_Central_pT_

std::vector< TH1F * >	h_norm_dz_ladder_

std::vector< TH1F * >	h_norm_dz_modZ_

std::vector< TH1F * >	h_norm_dz_pT_

TH1F *	h_nTracks

TH1F *	h_probeCharge_

TH1F *	h_probeChi2_

TH1F *	h_probed0RefitV_

TH1F *	h_probed3DRefitV_

TH1F *	h_probedxyRecoV_

TH1F *	h_probedxyRefitV_

TH1F *	h_probedzRecoV_

TH1F *	h_probedzRefitV_

TH1F *	h_probeEta_

TH1I *	h_probeHasBPixL1Overlap_

TH1F *	h_probeHits1D_

TH1F *	h_probeHits2D_

TH1F *	h_probeHits_

TH1F *	h_probeHitsInBPIX_

TH1F *	h_probeHitsInFPIX_

TH1F *	h_probeHitsInTEC_

TH1F *	h_probeHitsInTIB_

TH1F *	h_probeHitsInTID_

TH1F *	h_probeHitsInTOB_

TH1F *	h_probeL1Ladder_

TH1F *	h_probeL1Module_

TH1F *	h_probeNormChi2_

TH1F *	h_probeP_

TH1F *	h_probePhi_

TH1F *	h_probePt_

TH1F *	h_probePtRebin_

TH1F *	h_probeQoverP_

TH1F *	h_probeRecoVSigXY_

TH1F *	h_probeRecoVSigZ_

TH1F *	h_probeRefitVLogSig3D_

TH1F *	h_probeRefitVSig3D_

TH1F *	h_probeRefitVSigResZ_

TH1F *	h_probeRefitVSigXY_

TH1F *	h_probeRefitVSigZ_

TH1F *	h_probereszRefitV_

TH1F *	h_probesignIP2DRefitV_

TH1F *	h_probez0RefitV_

TH1F *	h_pTinfo

TH1F *	h_recoVtxChi2ndf_

TH1F *	h_recoVtxChi2Prob_

TH1F *	h_recoVtxNtracks_

TH1F *	h_recoVtxSumPt_

TH1I *	h_runFromConfig

TH1I *	h_runFromEvent

TH1F *	h_runNumber

TH1F *	h_xErrOfflineVertex

TH1F *	h_xOfflineVertex

TH1F *	h_yErrOfflineVertex

TH1F *	h_yOfflineVertex

TH1F *	h_zErrOfflineVertex

TH1F *	h_zOfflineVertex

int	hasRecVertex_ [nMaxtracks_]

std::map< std::string, TH1 * >	hDA

double	intLumi_

double	IP3DsigFromMyVertex_ [nMaxtracks_]

double	IPLsigFromMyVertex_ [nMaxtracks_]

double	IPTsigFromMyVertex_ [nMaxtracks_]

int	isGoodTrack_ [nMaxtracks_]

int	isHighPurity_ [nMaxtracks_]

bool	isPhase1_

bool	lightNtupleSwitch_

unsigned int	LuminosityBlockNumber_

TFileDirectory	MADTrendsDir

const double	maxPt_

TFileDirectory	Mean2DMapsDir

TFileDirectory	MeanTrendsDir

TFileDirectory	MedianTrendsDir

const double	minPt_

std::array< float, nPtBins_+1 >	mypT_bins_

std::vector< TH1F * >	n_d3DEtaResiduals

std::vector< TH1F * >	n_d3DPhiResiduals

TH1F *	n_d3DResidualsMap [nMaxBins_][nMaxBins_]

TH1F *	n_dxyBiasResidualsMap [nMaxBins_][nMaxBins_]

std::vector< TH1F * >	n_dxyEtaBiasResiduals

TH1F *	n_dxyEtaMADBiasTrend

TH1F *	n_dxyEtaMADTrend

TH1F *	n_dxyEtaMeanBiasTrend

TH1F *	n_dxyEtaMeanTrend

TH1F *	n_dxyEtaMedianBiasTrend

TH1F *	n_dxyEtaMedianTrend

std::vector< TH1F * >	n_dxyEtaResiduals

TH1F *	n_dxyEtaWidthBiasTrend

TH1F *	n_dxyEtaWidthTrend

TH1F *	n_dxyladderMeanTrend

TH1F *	n_dxyladderWidthTrend

TH2F *	n_dxyMeanBiasMap

TH2F *	n_dxyMeanMap

TH1F *	n_dxymodZMeanTrend

TH1F *	n_dxymodZWidthTrend

std::vector< TH1F * >	n_dxyPhiBiasResiduals

TH1F *	n_dxyPhiMADBiasTrend

TH1F *	n_dxyPhiMADTrend

TH1F *	n_dxyPhiMeanBiasTrend

TH1F *	n_dxyPhiMeanTrend

TH1F *	n_dxyPhiMedianBiasTrend

TH1F *	n_dxyPhiMedianTrend

std::vector< TH1F * >	n_dxyPhiResiduals

TH1F *	n_dxyPhiWidthBiasTrend

TH1F *	n_dxyPhiWidthTrend

TH1F *	n_dxypTCentralMeanTrend

TH1F *	n_dxypTCentralWidthTrend

TH1F *	n_dxypTMeanTrend

TH1F *	n_dxypTWidthTrend

TH1F *	n_dxyResidualsMap [nMaxBins_][nMaxBins_]

TH2F *	n_dxyVsEta

TH2F *	n_dxyVsPhi

TH2F *	n_dxyWidthBiasMap

TH2F *	n_dxyWidthMap

TH1F *	n_dzBiasResidualsMap [nMaxBins_][nMaxBins_]

std::vector< TH1F * >	n_dzEtaBiasResiduals

TH1F *	n_dzEtaMADBiasTrend

TH1F *	n_dzEtaMADTrend

TH1F *	n_dzEtaMeanBiasTrend

TH1F *	n_dzEtaMeanTrend

TH1F *	n_dzEtaMedianBiasTrend

TH1F *	n_dzEtaMedianTrend

std::vector< TH1F * >	n_dzEtaResiduals

TH1F *	n_dzEtaWidthBiasTrend

TH1F *	n_dzEtaWidthTrend

TH1F *	n_dzladderMeanTrend

TH1F *	n_dzladderWidthTrend

TH2F *	n_dzMeanBiasMap

TH2F *	n_dzMeanMap

TH1F *	n_dzmodZMeanTrend

TH1F *	n_dzmodZWidthTrend

std::vector< TH1F * >	n_dzPhiBiasResiduals

TH1F *	n_dzPhiMADBiasTrend

TH1F *	n_dzPhiMADTrend

TH1F *	n_dzPhiMeanBiasTrend

TH1F *	n_dzPhiMeanTrend

TH1F *	n_dzPhiMedianBiasTrend

TH1F *	n_dzPhiMedianTrend

std::vector< TH1F * >	n_dzPhiResiduals

TH1F *	n_dzPhiWidthBiasTrend

TH1F *	n_dzPhiWidthTrend

TH1F *	n_dzpTCentralMeanTrend

TH1F *	n_dzpTCentralWidthTrend

TH1F *	n_dzpTMeanTrend

TH1F *	n_dzpTWidthTrend

TH1F *	n_dzResidualsMap [nMaxBins_][nMaxBins_]

TH2F *	n_dzVsEta

TH2F *	n_dzVsPhi

TH2F *	n_dzWidthBiasMap

TH2F *	n_dzWidthMap

std::vector< TH1F * >	n_IP2DEtaResiduals

std::vector< TH1F * >	n_IP2DPhiResiduals

std::vector< TH1F * >	n_IP3DEtaResiduals

std::vector< TH1F * >	n_IP3DPhiResiduals

std::vector< TH1F * >	n_reszEtaResiduals

std::vector< TH1F * >	n_reszPhiResiduals

int	nBins_

int	nClus_

int	Nevt_

int	nhits1D_ [nMaxtracks_]

int	nhits2D_ [nMaxtracks_]

int	nhits_ [nMaxtracks_]

int	nhitsBPIX_ [nMaxtracks_]

int	nhitsFPIX_ [nMaxtracks_]

double	nHitsOfProbe_

int	nhitsTEC_ [nMaxtracks_]

int	nhitsTIB_ [nMaxtracks_]

int	nhitsTID_ [nMaxtracks_]

int	nhitsTOB_ [nMaxtracks_]

unsigned int	nLadders_ = 20

unsigned int	nModZ_ = 8

int	nOfflineVertices_

int	nTracks_

int	nTracksPerClus_

double	p_ [nMaxtracks_]

double	phi_ [nMaxtracks_]

double	pOfProbe_

double	pt_ [nMaxtracks_]

double	ptOfProbe_

double	qoverp_ [nMaxtracks_]

TTree *	rootTree_

bool	runControl_

std::vector< unsigned int >	runControlNumbers_

unsigned int	RunNumber_

bool	storeNtuple_

float	sumOfWeightsUnbiasedVertex_ [nMaxtracks_]

edm::EDGetTokenT< reco::BeamSpot >	theBeamspotToken

edm::ParameterSet	theConfig

PVValHelper::histodetails	theDetails_

double	theta_ [nMaxtracks_]

TrackClusterizerInZ *	theTrackClusterizer_

edm::EDGetTokenT< reco::TrackCollection >	theTrackCollectionToken

TrackFilterForPVFindingBase *	theTrackFilter_

edm::EDGetTokenT< reco::VertexCollection >	theVertexCollectionToken

int	tracksUsedForVertexing_ [nMaxtracks_]

bool	useTracksFromRecoVtx_

double	vertexZMax_

TFileDirectory	Width2DMapsDir

TFileDirectory	WidthTrendsDir

double	wxy2_

double	xErrOfflineVertex_

double	xOfflineVertex_

double	xPCA_ [nMaxtracks_]

double	xUnbiasedVertex_ [nMaxtracks_]

double	yErrOfflineVertex_

double	yOfflineVertex_

double	yPCA_ [nMaxtracks_]

double	yUnbiasedVertex_ [nMaxtracks_]

double	zErrOfflineVertex_

double	zOfflineVertex_

double	zPCA_ [nMaxtracks_]

double	zUnbiasedVertex_ [nMaxtracks_]

Static Private Attributes
static const int	cmToum = 10000

static const int	nMaxBins_ = 100

static const int	nMaxtracks_ = 1000

static const int	nPtBins_ = 48

Additional Inherited Members
Public Types inherited from edm::one::EDAnalyzerBase
typedef EDAnalyzerBase	ModuleType

Public Types inherited from edm::EDConsumerBase
typedef ProductLabels	Labels

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

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &descriptions)

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 ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto	esConsumes ()

template<typename ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto	esConsumes (ESInputTag const &tag)

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: Validate alignment constants using unbiased vertex residuals

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

Definition at line 67 of file PrimaryVertexValidation.h.

Constructor & Destructor Documentation

PrimaryVertexValidation::PrimaryVertexValidation ( const edm::ParameterSet & iConfig )

explicit

Definition at line 69 of file PrimaryVertexValidation.cc.

References PVValHelper::d3D, debug_, PVValHelper::dx, PVValHelper::dxy, PVValHelper::END_OF_PLOTS, PVValHelper::END_OF_TYPES, PVValHelper::eta, etaOfProbe_, PVValHelper::generateBins(), PVValHelper::histodetails::getHigh(), PVValHelper::histodetails::getLow(), edm::ParameterSet::getParameter(), PVValHelper::getTypeString(), edm::ParameterSet::getUntrackedParameter(), PVValHelper::getVarString(), PVValHelper::histodetails::histobins, mps_fire::i, TFileService::kSharedResource, PVValHelper::ladder, maxPt_, minPt_, PVValHelper::modZ, mypT_bins_, nBins_, PVValHelper::norm_d3D, PVValHelper::phi, PVValHelper::pT, PVValHelper::pTCentral, PVValHelper::histodetails::range, runControlNumbers_, PVValHelper::histodetails::setMap(), AlCaHLTBitMon_QueryRunRegistry::string, theBeamspotToken, theDetails_, theTrackClusterizer_, theTrackCollectionToken, theTrackFilter_, theVertexCollectionToken, and PVValHelper::histodetails::trendbins.

                                                                               :
   storeNtuple_(iConfig.getParameter<bool>("storeNtuple")),
   lightNtupleSwitch_(iConfig.getParameter<bool>("isLightNtuple")),
   useTracksFromRecoVtx_(iConfig.getParameter<bool>("useTracksFromRecoVtx")),
   vertexZMax_(iConfig.getUntrackedParameter<double>("vertexZMax",99.)),
   intLumi_(iConfig.getUntrackedParameter<double>("intLumi",0.)),
   askFirstLayerHit_(iConfig.getParameter<bool>("askFirstLayerHit")),
   doBPix_(iConfig.getUntrackedParameter<bool>("doBPix",true)),
   doFPix_(iConfig.getUntrackedParameter<bool>("doFPix",true)),
   ptOfProbe_(iConfig.getUntrackedParameter<double>("probePt",0.)),
   pOfProbe_(iConfig.getUntrackedParameter<double>("probeP",0.)),
   etaOfProbe_(iConfig.getUntrackedParameter<double>("probeEta",2.4)),
   nHitsOfProbe_(iConfig.getUntrackedParameter<double>("probeNHits",0.)),
   nBins_(iConfig.getUntrackedParameter<int>("numberOfBins",24)),
   minPt_(iConfig.getUntrackedParameter<double>("minPt",1.)),
   maxPt_(iConfig.getUntrackedParameter<double>("maxPt",20.)),
   debug_(iConfig.getParameter<bool>("Debug")),
   runControl_(iConfig.getUntrackedParameter<bool>("runControl",false)),
   forceBeamSpotContraint_(iConfig.getUntrackedParameter<bool>("forceBeamSpot",false))
 {
   
   // now do what ever initialization is needed
   // initialize phase space boundaries
   
   usesResource(TFileService::kSharedResource);
   
   std::vector<unsigned int> defaultRuns;
   defaultRuns.push_back(0);
   runControlNumbers_ = iConfig.getUntrackedParameter<std::vector<unsigned int> >("runControlNumber",defaultRuns);
 
   edm::InputTag TrackCollectionTag_ = iConfig.getParameter<edm::InputTag>("TrackCollectionTag");
   theTrackCollectionToken = consumes<reco::TrackCollection>(TrackCollectionTag_);
 
   edm::InputTag VertexCollectionTag_ = iConfig.getParameter<edm::InputTag>("VertexCollectionTag");
   theVertexCollectionToken = consumes<reco::VertexCollection>(VertexCollectionTag_);
 
   edm::InputTag BeamspotTag_ = edm::InputTag("offlineBeamSpot");
   theBeamspotToken = consumes<reco::BeamSpot>(BeamspotTag_);
 
   // select and configure the track filter 
   theTrackFilter_= new TrackFilterForPVFinding(iConfig.getParameter<edm::ParameterSet>("TkFilterParameters") );
   // select and configure the track clusterizer  
   std::string clusteringAlgorithm=iConfig.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<std::string>("algorithm");
   if (clusteringAlgorithm=="gap"){
     theTrackClusterizer_ = new GapClusterizerInZ(iConfig.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkGapClusParameters"));
   }else if(clusteringAlgorithm=="DA"){
     theTrackClusterizer_ = new DAClusterizerInZ(iConfig.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkDAClusParameters"));
     // provide the vectorized version of the clusterizer, if supported by the build
   } else if(clusteringAlgorithm=="DA_vect") {
     theTrackClusterizer_ = new DAClusterizerInZ_vect(iConfig.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkDAClusParameters"));
   } else {
     throw VertexException("PrimaryVertexProducerAlgorithm: unknown clustering algorithm: " + clusteringAlgorithm);  
   }
 
   theDetails_.histobins = 500;
   theDetails_.setMap(PVValHelper::dxy,PVValHelper::phi      ,-2000.,2000.); 
   theDetails_.setMap(PVValHelper::dxy,PVValHelper::eta      ,-3000.,3000.);
   theDetails_.setMap(PVValHelper::dxy,PVValHelper::pT       ,-1000.,1000.);
   theDetails_.setMap(PVValHelper::dxy,PVValHelper::pTCentral,-1000.,1000.);
   theDetails_.setMap(PVValHelper::dxy,PVValHelper::ladder   ,-1000.,1000.);
   theDetails_.setMap(PVValHelper::dxy,PVValHelper::modZ     ,-1000.,1000.);
 
   for (int i = PVValHelper::phi; i < PVValHelper::END_OF_PLOTS; i++ ){
     for (int j = PVValHelper::dx; j < PVValHelper::END_OF_TYPES; j++ ){
 
       auto plot_index = static_cast<PVValHelper::plotVariable>(i);
       auto res_index  = static_cast<PVValHelper::residualType>(j);
       
       if(debug_){
     edm::LogInfo("PrimaryVertexValidation")<<"==> "<<std::get<0>(PVValHelper::getTypeString(res_index)) << " "<< std::setw(10)<< std::get<0>(PVValHelper::getVarString(plot_index))<<std::endl;
       }
       if(res_index!=PVValHelper::d3D && res_index!=PVValHelper::norm_d3D)
     theDetails_.setMap(res_index,plot_index,theDetails_.getLow(PVValHelper::dxy,plot_index),theDetails_.getHigh(PVValHelper::dxy,plot_index));
       else
     theDetails_.setMap(res_index,plot_index,0.,theDetails_.getHigh(PVValHelper::dxy,plot_index));
     }
   }
 
   edm::LogVerbatim("PrimaryVertexValidation") <<"######################################";
   for (const auto & it : theDetails_.range){
     edm::LogVerbatim("PrimaryVertexValidation")<< "|" <<std::setw(10) << std::get<0>(PVValHelper::getTypeString(it.first.first)) << "|" << std::setw(10)<< std::get<0>(PVValHelper::getVarString(it.first.second)) << "| (" << std::setw(5)<< it.second.first << ";" <<std::setw(5)<< it.second.second << ") |"<<std::endl;
   }
 
   theDetails_.trendbins[PVValHelper::phi] = PVValHelper::generateBins(nBins_+1,-180.,360.);
   theDetails_.trendbins[PVValHelper::eta] = PVValHelper::generateBins(nBins_+1,-etaOfProbe_,2*etaOfProbe_);
 
   if(debug_){
     edm::LogVerbatim("PrimaryVertexValidation") << "etaBins: ";
     for (auto ieta: theDetails_.trendbins[PVValHelper::eta]) {
       edm::LogVerbatim("PrimaryVertexValidation") << ieta << " ";
     }
     edm::LogVerbatim("PrimaryVertexValidation") << "\n";
     
     edm::LogVerbatim("PrimaryVertexValidation") << "phiBins: ";
     for (auto iphi: theDetails_.trendbins[PVValHelper::phi]) {
       edm::LogVerbatim("PrimaryVertexValidation") << iphi << " ";
     }
     edm::LogVerbatim("PrimaryVertexValidation") << "\n";
   }
   
   // create the bins of the pT-binned distributions 
 
   mypT_bins_ = PVValHelper::makeLogBins<float,nPtBins_>(minPt_,maxPt_);
 
   std::string toOutput="";
   for (auto ptbin: mypT_bins_){
     toOutput+=" ";
     toOutput+=std::to_string(ptbin);
     toOutput+=",";
   }
   
   edm::LogVerbatim("PrimaryVertexValidation") <<"######################################\n";
   edm::LogVerbatim("PrimaryVertexValidation") <<"The pT binning is: [" << toOutput << "] \n";
 
 }

PrimaryVertexValidation::~PrimaryVertexValidation ( )

override

Definition at line 186 of file PrimaryVertexValidation.cc.

References theTrackClusterizer_, and theTrackFilter_.

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

Member Function Documentation

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

private

Referenced by counter.Counter::register(), SequenceTypes.Task::remove(), and SequenceTypes.Task::replace().

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

overrideprivate

Definition at line 201 of file PrimaryVertexValidation.cc.

 {
 
   using namespace std;
   using namespace reco;
   using namespace IPTools;
 
   if (!isBFieldConsistentWithMode(iSetup)) {
     edm::LogWarning("PrimaryVertexValidation") << "*********************************************************************************\n" 
                            << "* The configuration (ptOfProbe > " << ptOfProbe_ << "GeV) is not correctly set for current value of magnetic field \n" 
                            << "* Switching it to 0. !!! \n"
                            << "*********************************************************************************"<< std::endl;
     ptOfProbe_=0.;
   }
 
   if(nBins_!=24 && debug_){ 
     edm::LogInfo("PrimaryVertexValidation")<<"Using: "<<nBins_<<" bins plots";
   }
   
   bool passesRunControl = false;
 
   if(runControl_){
     for(const auto & runControlNumber : runControlNumbers_){
       if(iEvent.eventAuxiliary().run() == runControlNumber){ 
     if (debug_){
       edm::LogInfo("PrimaryVertexValidation")<<" run number: "<<iEvent.eventAuxiliary().run()<<" keeping run:"<<runControlNumber;
     }
     passesRunControl = true;
     break;
       }
     }
     if (!passesRunControl) return;
   }
 
   Nevt_++;  
 
   //=======================================================
   // Initialize Root-tuple variables
   //=======================================================
 
   SetVarToZero();
  
   //=======================================================
   // Retrieve the Magnetic Field information
   //=======================================================
 
   edm::ESHandle<MagneticField> theMGField;
   iSetup.get<IdealMagneticFieldRecord>().get( theMGField );
 
   //=======================================================
   // Retrieve the Tracking Geometry information
   //=======================================================
 
   edm::ESHandle<GlobalTrackingGeometry> theTrackingGeometry;
   iSetup.get<GlobalTrackingGeometryRecord>().get( theTrackingGeometry );
 
   //=======================================================
   // Retrieve geometry information
   //=======================================================
  
   edm::LogInfo("read tracker geometry...");
   edm::ESHandle<TrackerGeometry> pDD;
   iSetup.get<TrackerDigiGeometryRecord>().get( pDD );
   edm::LogInfo("tracker geometry read")<<"There are: "<< pDD->dets().size() <<" detectors";
 
   // switch on the phase1 
   if( (pDD->isThere(GeomDetEnumerators::P1PXB)) || 
       (pDD->isThere(GeomDetEnumerators::P1PXEC)) ) {
     isPhase1_ = true;
     nLadders_ = 12;
                         
     if(h_dxy_ladderOverlap_.size()!=nLadders_){
 
       PVValHelper::shrinkHistVectorToFit(h_dxy_ladderOverlap_,nLadders_);
       PVValHelper::shrinkHistVectorToFit(h_dxy_ladderNoOverlap_,nLadders_);
       PVValHelper::shrinkHistVectorToFit(h_dxy_ladder_,nLadders_);    
       PVValHelper::shrinkHistVectorToFit(h_dz_ladder_,nLadders_);     
       PVValHelper::shrinkHistVectorToFit(h_norm_dxy_ladder_,nLadders_);  
       PVValHelper::shrinkHistVectorToFit(h_norm_dz_ladder_,nLadders_);   
       
       if (debug_){
     edm::LogInfo("PrimaryVertexValidation")<<"checking size:"<<h_dxy_ladder_.size()<<std::endl;
       }
     }
 
     if (debug_){
       edm::LogInfo("PrimaryVertexValidation")<<" pixel phase1 setup, nLadders: "<<nLadders_;
     }
     
   } else {
     isPhase1_ = false;
     nLadders_ = 20;
     if (debug_){
       edm::LogInfo("PrimaryVertexValidation")<<" pixel phase0 setup, nLadders: "<<nLadders_;
     }
   }
 
   if(isPhase1_){
     etaOfProbe_ = std::min(etaOfProbe_,2.7);
   } else {
     etaOfProbe_ = std::min(etaOfProbe_,2.5);
   }
 
   if(h_etaMax->GetEntries()==0.){
     h_etaMax->SetBinContent(1.,etaOfProbe_);
     h_nbins->SetBinContent(1.,nBins_);
     h_nLadders->SetBinContent(1.,nLadders_);
     h_pTinfo->SetBinContent(1.,mypT_bins_.size());
     h_pTinfo->SetBinContent(2.,minPt_);
     h_pTinfo->SetBinContent(3.,maxPt_);
   }
 
   //=======================================================
   // Retrieve the Transient Track Builder information
   //=======================================================
 
   edm::ESHandle<TransientTrackBuilder> theB_;
   iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder",theB_);
   double fBfield_=((*theB_).field()->inTesla(GlobalPoint(0.,0.,0.))).z();
 
   //=======================================================
   // Retrieve the Track information
   //=======================================================
   
   edm::Handle<TrackCollection>  trackCollectionHandle;
   iEvent.getByToken(theTrackCollectionToken, trackCollectionHandle);
   if(!trackCollectionHandle.isValid()) return;
   auto const & tracks = *trackCollectionHandle;
 
   //=======================================================
   // Retrieve tracker topology from geometry
   //=======================================================
 
   edm::ESHandle<TrackerTopology> tTopoHandle;
   iSetup.get<TrackerTopologyRcd>().get(tTopoHandle);
   const TrackerTopology* const tTopo = tTopoHandle.product();
 
   //=======================================================
   // Retrieve offline vartex information (only for reco)
   //=======================================================
  
   //edm::Handle<VertexCollection> vertices;
   edm::Handle<std::vector<Vertex> > vertices;
 
   try {
     iEvent.getByToken(theVertexCollectionToken, vertices);
   } catch ( cms::Exception& er ) {
     LogTrace("PrimaryVertexValidation")<<"caught std::exception "<<er.what()<<std::endl;  
   }
 
   std::vector<Vertex> vsorted = *(vertices);
   // sort the vertices by number of tracks in descending order
   // use chi2 as tiebreaker
   std::sort( vsorted.begin(), vsorted.end(), PrimaryVertexValidation::vtxSort );
   
   // skip events with no PV, this should not happen
 
   if( vsorted.empty()) return;
 
   // skip events failing vertex cut
   if( std::abs(vsorted[0].z()) > vertexZMax_ ) return; 
   
   if ( vsorted[0].isValid() ) {
     xOfflineVertex_ = (vsorted)[0].x();
     yOfflineVertex_ = (vsorted)[0].y();
     zOfflineVertex_ = (vsorted)[0].z();
 
     xErrOfflineVertex_ = (vsorted)[0].xError();
     yErrOfflineVertex_ = (vsorted)[0].yError();
     zErrOfflineVertex_ = (vsorted)[0].zError();
   }
 
   h_xOfflineVertex->Fill(xOfflineVertex_);    
   h_yOfflineVertex->Fill(yOfflineVertex_);     
   h_zOfflineVertex->Fill(zOfflineVertex_);     
   h_xErrOfflineVertex->Fill(xErrOfflineVertex_);  
   h_yErrOfflineVertex->Fill(yErrOfflineVertex_);  
   h_zErrOfflineVertex->Fill(zErrOfflineVertex_);  
 
   unsigned int vertexCollectionSize = vsorted.size();
   int nvvertex = 0;
   
   for (unsigned int i=0; i<vertexCollectionSize; i++) {
     const Vertex& vertex = vsorted.at(i);
     if (vertex.isValid()) nvvertex++;
   }
 
   nOfflineVertices_ = nvvertex;
   h_nOfflineVertices->Fill(nvvertex);
 
   if ( !vsorted.empty() && useTracksFromRecoVtx_ ) {
    
     double sumpt    = 0;
     size_t ntracks  = 0;
     double chi2ndf  = 0.; 
     double chi2prob = 0.;
 
     if (!vsorted.at(0).isFake()) {
       
       Vertex pv = vsorted.at(0);
       
       ntracks  = pv.tracksSize();
       chi2ndf  = pv.normalizedChi2();
       chi2prob = TMath::Prob(pv.chi2(),(int)pv.ndof());
       
       h_recoVtxNtracks_->Fill(ntracks);   
       h_recoVtxChi2ndf_->Fill(chi2ndf);    
       h_recoVtxChi2Prob_->Fill(chi2prob);   
       
       for (Vertex::trackRef_iterator itrk = pv.tracks_begin();itrk != pv.tracks_end(); ++itrk) {
     double pt = (**itrk).pt();
     sumpt += pt*pt;
     
     const math::XYZPoint myVertex(pv.position().x(),pv.position().y(),pv.position().z());
     
     double dxyRes = (**itrk).dxy(myVertex);
     double dzRes  = (**itrk).dz(myVertex);
     
     double dxy_err = (**itrk).dxyError();
     double dz_err  = (**itrk).dzError();
     
     float trackphi = ((**itrk).phi())*(180/M_PI);
     float tracketa = (**itrk).eta();
     
     for(int i=0; i<nBins_; i++){
       
       float phiF = theDetails_.trendbins[PVValHelper::phi][i];
       float phiL = theDetails_.trendbins[PVValHelper::phi][i+1];
       
       float etaF = theDetails_.trendbins[PVValHelper::eta][i];
       float etaL = theDetails_.trendbins[PVValHelper::eta][i+1];
 
       if(tracketa >= etaF && tracketa < etaL ){
 
         PVValHelper::fillByIndex(a_dxyEtaBiasResiduals,i,dxyRes*cmToum);
         PVValHelper::fillByIndex(a_dzEtaBiasResiduals,i,dzRes*cmToum); 
         PVValHelper::fillByIndex(n_dxyEtaBiasResiduals,i,(dxyRes)/dxy_err);
         PVValHelper::fillByIndex(n_dzEtaBiasResiduals,i,(dzRes)/dz_err);
 
       }
 
       if(trackphi >= phiF && trackphi < phiL ){ 
 
         PVValHelper::fillByIndex(a_dxyPhiBiasResiduals,i,dxyRes*cmToum);
         PVValHelper::fillByIndex(a_dzPhiBiasResiduals,i,dzRes*cmToum); 
         PVValHelper::fillByIndex(n_dxyPhiBiasResiduals,i,(dxyRes)/dxy_err);
         PVValHelper::fillByIndex(n_dzPhiBiasResiduals,i,(dzRes)/dz_err); 
         
         for(int j=0; j<nBins_; j++){
           
           float etaJ = theDetails_.trendbins[PVValHelper::eta][j];
           float etaK = theDetails_.trendbins[PVValHelper::eta][j+1];
           
           if(tracketa >= etaJ && tracketa < etaK ){
         
         a_dxyBiasResidualsMap[i][j]->Fill(dxyRes*cmToum); 
         a_dzBiasResidualsMap[i][j]->Fill(dzRes*cmToum);   
         
         n_dxyBiasResidualsMap[i][j]->Fill((dxyRes)/dxy_err); 
         n_dzBiasResidualsMap[i][j]->Fill((dzRes)/dz_err);  
         
           }
         }
       }     
     }
       }
       
       h_recoVtxSumPt_->Fill(sumpt);   
 
     }
   }
 
   //=======================================================
   // Retrieve Beamspot information
   //=======================================================
 
   BeamSpot beamSpot;
   edm::Handle<BeamSpot> beamSpotHandle;
   iEvent.getByToken(theBeamspotToken, beamSpotHandle);
     
   if ( beamSpotHandle.isValid() ) {
     beamSpot = *beamSpotHandle;
     BSx0_    = beamSpot.x0();
     BSy0_    = beamSpot.y0();
     BSz0_    = beamSpot.z0();
     Beamsigmaz_ = beamSpot.sigmaZ();    
     Beamdxdz_   = beamSpot.dxdz();       
     BeamWidthX_ = beamSpot.BeamWidthX();
     BeamWidthY_ = beamSpot.BeamWidthY();
 
     wxy2_=TMath::Power(BeamWidthX_,2)+TMath::Power(BeamWidthY_,2);
 
   } else {
     edm::LogWarning("PrimaryVertexValidation")<<"No BeamSpot found!";
   }
 
   h_BSx0->Fill(BSx0_);      
   h_BSy0->Fill(BSy0_);     
   h_BSz0->Fill(BSz0_);      
   h_Beamsigmaz->Fill(Beamsigmaz_);
   h_BeamWidthX->Fill(BeamWidthX_);
   h_BeamWidthY->Fill(BeamWidthY_);
 
   if(debug_)
     edm::LogInfo("PrimaryVertexValidation")<<"Beamspot x:" <<BSx0_<<" y:"<<BSy0_<<" z:"<<BSz0_;
    
   //=======================================================
   // Starts here ananlysis
   //=======================================================
   
   RunNumber_=iEvent.eventAuxiliary().run();
   h_runNumber->Fill(RunNumber_);
 
   if(h_runFromEvent->GetEntries()==0){
     h_runFromEvent->SetBinContent(1,RunNumber_);
   }
 
   LuminosityBlockNumber_=iEvent.eventAuxiliary().luminosityBlock();
   EventNumber_=iEvent.eventAuxiliary().id().event();
     
   if(debug_)
     edm::LogInfo("PrimaryVertexValidation")<<" looping over "<<trackCollectionHandle->size()<< "tracks";
 
   h_nTracks->Fill(trackCollectionHandle->size()); 
 
   //======================================================
   // Interface RECO tracks to vertex reconstruction
   //======================================================
  
   std::vector<TransientTrack> t_tks;
   for (const auto & track : tracks){
     TransientTrack tt = theB_->build(&(track));  
     tt.setBeamSpot(beamSpot);
     t_tks.push_back(tt);
   
   }
   
   if(debug_) {
     edm::LogInfo("PrimaryVertexValidation") << "Found: " << t_tks.size() << " reconstructed tracks";
   }
   
   //======================================================
   // select the tracks
   //======================================================
 
   std::vector<TransientTrack> seltks = theTrackFilter_->select(t_tks);
     
   //======================================================
   // clusterize tracks in Z
   //======================================================
 
   vector< vector<TransientTrack> > clusters = theTrackClusterizer_->clusterize(seltks);
   
   if (debug_){
     edm::LogInfo("PrimaryVertexValidation")<<" looping over: "<< clusters.size() << " clusters  from " << t_tks.size() << " selected tracks";
   }
   
   nClus_=clusters.size();  
   h_nClus->Fill(nClus_);
 
   //======================================================
   // Starts loop on clusters 
   //======================================================
   for (const auto & iclus : clusters){
 
     nTracksPerClus_=0;
 
     unsigned int i=0;
     for(const auto & theTTrack : iclus)
       {
     i++;
 
     if ( nTracks_ >= nMaxtracks_ ) {
       edm::LogError("PrimaryVertexValidation")<<" Warning - Number of tracks: " << nTracks_ << " , greater than " << nMaxtracks_;
       continue;
     }
     
     const Track & theTrack = theTTrack.track();
 
     pt_[nTracks_]       = theTrack.pt();
     p_[nTracks_]        = theTrack.p();
     nhits_[nTracks_]    = theTrack.numberOfValidHits();
     eta_[nTracks_]      = theTrack.eta();
     theta_[nTracks_]    = theTrack.theta();
     phi_[nTracks_]      = theTrack.phi();
     chi2_[nTracks_]     = theTrack.chi2();
     chi2ndof_[nTracks_] = theTrack.normalizedChi2();
     charge_[nTracks_]   = theTrack.charge();
     qoverp_[nTracks_]   = theTrack.qoverp();
     dz_[nTracks_]       = theTrack.dz();
     dxy_[nTracks_]      = theTrack.dxy();
     
     TrackBase::TrackQuality _trackQuality = TrackBase::qualityByName("highPurity");
     isHighPurity_[nTracks_] = theTrack.quality(_trackQuality);
     
     math::XYZPoint point(BSx0_,BSy0_,BSz0_);
     dxyBs_[nTracks_]    = theTrack.dxy(point);
     dzBs_[nTracks_]     = theTrack.dz(point);
 
     xPCA_[nTracks_]     = theTrack.vertex().x();
     yPCA_[nTracks_]     = theTrack.vertex().y();
     zPCA_[nTracks_]     = theTrack.vertex().z(); 
     
     //=======================================================
     // Retrieve rechit information
     //=======================================================  
     
     const reco::HitPattern& hits = theTrack.hitPattern();
     
     int nRecHit1D=0;
     int nRecHit2D=0;
     int nhitinTIB  = hits.numberOfValidStripTIBHits(); 
     int nhitinTOB  = hits.numberOfValidStripTOBHits(); 
     int nhitinTID  = hits.numberOfValidStripTIDHits(); 
     int nhitinTEC  = hits.numberOfValidStripTECHits();
     int nhitinBPIX = hits.numberOfValidPixelBarrelHits();
     int nhitinFPIX = hits.numberOfValidPixelEndcapHits();
     for (trackingRecHit_iterator iHit = theTTrack.recHitsBegin(); iHit != theTTrack.recHitsEnd(); ++iHit) {
       if((*iHit)->isValid()) {  
         
         if (this->isHit2D(**iHit)) {++nRecHit2D;}
         else {++nRecHit1D; }
       }
     } 
     
     nhits1D_[nTracks_]     = nRecHit1D;
     nhits2D_[nTracks_]     = nRecHit2D;
     nhitsBPIX_[nTracks_]   = nhitinBPIX;
     nhitsFPIX_[nTracks_]   = nhitinFPIX;
     nhitsTIB_[nTracks_]    = nhitinTIB;
     nhitsTID_[nTracks_]    = nhitinTID;
     nhitsTOB_[nTracks_]    = nhitinTOB;
     nhitsTEC_[nTracks_]    = nhitinTEC;
     
     //=======================================================
     // Good tracks for vertexing selection
     //=======================================================  
 
     bool pass = true;
     if(askFirstLayerHit_) pass = this->hasFirstLayerPixelHits(theTTrack);
     if (pass 
         && (theTrack.pt() >=ptOfProbe_) 
         && std::abs(theTrack.eta()) <= etaOfProbe_ 
         && (theTrack.numberOfValidHits())>=nHitsOfProbe_
         && (theTrack.p()) >= pOfProbe_ ){
       isGoodTrack_[nTracks_]=1;
     }
       
     //=======================================================
     // Fit unbiased vertex
     //=======================================================
     
     vector<TransientTrack> theFinalTracks;
     theFinalTracks.clear();
       
     for (const auto & tk : iclus) {
 
       pass = this->hasFirstLayerPixelHits(tk);
       if (pass){
         if( tk == theTTrack ) continue;
         else {
           theFinalTracks.push_back(tk);
         }
       }
     }
     
     if(theFinalTracks.size() > 1){
         
       if(debug_)
         edm::LogInfo("PrimaryVertexValidation")<<"Transient Track Collection size: "<<theFinalTracks.size();
       try{
 
         //AdaptiveVertexFitter* theFitter = new AdaptiveVertexFitter;
         auto theFitter = std::unique_ptr<VertexFitter<5> >( new AdaptiveVertexFitter());
         TransientVertex theFittedVertex;
 
         if(forceBeamSpotContraint_){
           theFittedVertex = theFitter->vertex(theFinalTracks,beamSpot);  // if you want the beam constraint
         } else {
           theFittedVertex = theFitter->vertex(theFinalTracks);
         }
 
 
 
         double totalTrackWeights=0;
         if(theFittedVertex.isValid ()){
 
           
           if(theFittedVertex.hasTrackWeight()){
         for(const auto & theFinalTrack : theFinalTracks){
           sumOfWeightsUnbiasedVertex_[nTracks_] += theFittedVertex.trackWeight(theFinalTrack);
           totalTrackWeights+= theFittedVertex.trackWeight(theFinalTrack);
           h_fitVtxTrackWeights_->Fill(theFittedVertex.trackWeight(theFinalTrack));
         }
           }
           
           h_fitVtxTrackAverageWeight_->Fill(totalTrackWeights/theFinalTracks.size());
 
           const math::XYZPoint theRecoVertex(xOfflineVertex_,yOfflineVertex_,zOfflineVertex_);      
           const math::XYZPoint myVertex(theFittedVertex.position().x(),theFittedVertex.position().y(),theFittedVertex.position().z());
 
           const Vertex vertex = theFittedVertex;
           fillTrackHistos(hDA,"all",&theTTrack,vertex,beamSpot,fBfield_);
 
           hasRecVertex_[nTracks_]    = 1;
           xUnbiasedVertex_[nTracks_] = theFittedVertex.position().x();
           yUnbiasedVertex_[nTracks_] = theFittedVertex.position().y();
           zUnbiasedVertex_[nTracks_] = theFittedVertex.position().z();
           
           chi2normUnbiasedVertex_[nTracks_] = theFittedVertex.normalisedChiSquared();
           chi2UnbiasedVertex_[nTracks_]     = theFittedVertex.totalChiSquared();
           DOFUnbiasedVertex_[nTracks_]      = theFittedVertex.degreesOfFreedom();   
           chi2ProbUnbiasedVertex_[nTracks_] = TMath::Prob(theFittedVertex.totalChiSquared(),(int)theFittedVertex.degreesOfFreedom());
           tracksUsedForVertexing_[nTracks_] = theFinalTracks.size();
 
           h_fitVtxNtracks_->Fill(theFinalTracks.size());        
           h_fitVtxChi2_->Fill(theFittedVertex.totalChiSquared());
           h_fitVtxNdof_->Fill(theFittedVertex.degreesOfFreedom());
           h_fitVtxChi2ndf_->Fill(theFittedVertex.normalisedChiSquared());        
           h_fitVtxChi2Prob_->Fill(TMath::Prob(theFittedVertex.totalChiSquared(),(int)theFittedVertex.degreesOfFreedom()));       
                 
           // from my Vertex
           double dxyFromMyVertex = theTrack.dxy(myVertex);
           double dzFromMyVertex  = theTrack.dz(myVertex);
 
           GlobalPoint vert(theFittedVertex.position().x(),theFittedVertex.position().y(),theFittedVertex.position().z());
 
           //FreeTrajectoryState theTrackNearVertex = theTTrack.trajectoryStateClosestToPoint(vert).theState();
           //double dz_err = sqrt(theFittedVertex.positionError().czz() + theTrackNearVertex.cartesianError().position().czz());      
           //double dz_err = hypot(theTrack.dzError(),theFittedVertex.positionError().czz());
           
           double dz_err = sqrt(std::pow(theTrack.dzError(), 2) + theFittedVertex.positionError().czz());
 
 
           // PV2D 
           std::pair<bool,Measurement1D> s_ip2dpv = signedTransverseImpactParameter(theTTrack,
                                                GlobalVector(theTrack.px(),
                                                     theTrack.py(),
                                                     theTrack.pz()),
                                                theFittedVertex);            
           
           double s_ip2dpv_corr = s_ip2dpv.second.value();
           double s_ip2dpv_err  = s_ip2dpv.second.error();
           
           // PV3D
           std::pair<bool, Measurement1D> s_ip3dpv = signedImpactParameter3D(theTTrack,          
                                         GlobalVector(theTrack.px(),  
                                                  theTrack.py(),  
                                                  theTrack.pz()), 
                                         theFittedVertex);            
           
           double s_ip3dpv_corr = s_ip3dpv.second.value();
           double s_ip3dpv_err  = s_ip3dpv.second.error();
 
           // PV3D absolute
           std::pair<bool,Measurement1D> ip3dpv = absoluteImpactParameter3D(theTTrack,theFittedVertex);
           double ip3d_corr = ip3dpv.second.value(); 
           double ip3d_err  = ip3dpv.second.error(); 
           
           // with respect to any specified vertex, such as primary vertex
           TrajectoryStateClosestToPoint traj = (theTTrack).trajectoryStateClosestToPoint(vert);
 
           GlobalPoint refPoint = traj.position();
           GlobalPoint cPToVtx  = traj.theState().position();
 
           float my_dx = refPoint.x() - myVertex.x();
           float my_dy = refPoint.y() - myVertex.y();
 
           float my_dx2 = cPToVtx.x() - myVertex.x();
           float my_dy2 = cPToVtx.y() - myVertex.y();
 
           float my_dxy = std::sqrt(my_dx*my_dx + my_dy*my_dy);
 
           double d0 = traj.perigeeParameters().transverseImpactParameter();
           //double d0_error = traj.perigeeError().transverseImpactParameterError();
           double z0 = traj.perigeeParameters().longitudinalImpactParameter();
           double z0_error = traj.perigeeError().longitudinalImpactParameterError();
 
           if(debug_){
         edm::LogInfo("PrimaryVertexValidation")<< "my_dx:"   << my_dx  
                                << " my_dy:"  << my_dy  
                                << " my_dxy:" << my_dxy
                                << " my_dx2:" << my_dx2  
                                << " my_dy2:" << my_dy2
                                << " d0: "    << d0
                                << " dxyFromVtx:" << dxyFromMyVertex << "\n"
                                << " ============================== "<< "\n"
                                << "diff1:"    << std::abs(d0) - std::abs(my_dxy) << "\n"
                                << "diff2:"    << std::abs(d0) - std::abs(dxyFromMyVertex) << "\n"
                                << "diff3:"    << (my_dx - my_dx2) << " " << (my_dy - my_dy2) << "\n" 
                                << std::endl;    
           }
 
           // define IPs
          
           dxyFromMyVertex_[nTracks_]      = dxyFromMyVertex; 
           dxyErrorFromMyVertex_[nTracks_] = s_ip2dpv_err;
           IPTsigFromMyVertex_[nTracks_]   = dxyFromMyVertex/s_ip2dpv_err;
  
           dzFromMyVertex_[nTracks_]       = dzFromMyVertex;  
           dzErrorFromMyVertex_[nTracks_]  = dz_err;
           IPLsigFromMyVertex_[nTracks_]   = dzFromMyVertex/dz_err;
 
           d3DFromMyVertex_[nTracks_]      = ip3d_corr;
           d3DErrorFromMyVertex_[nTracks_] = ip3d_err;
           IP3DsigFromMyVertex_[nTracks_]  = (ip3d_corr/ip3d_err);
                   
           // fill directly the histograms of residuals
           
           float trackphi = (theTrack.phi())*(180./M_PI);
           float tracketa = theTrack.eta();
           float trackpt  = theTrack.pt();
           float trackp   = theTrack.p();
           float tracknhits = theTrack.numberOfValidHits();
 
           // determine the module number and ladder
 
           int ladder_num = -1.;
           int module_num = -1.;
           int L1BPixHitCount = 0;
 
           for(auto const& hit : theTrack.recHits()) {
         const DetId& detId = hit->geographicalId();
         unsigned int subid = detId.subdetId();
         
         if(hit->isValid() && ( subid == PixelSubdetector::PixelBarrel ) ) {
           int layer = tTopo->pxbLayer(detId);
           if(layer==1){
             L1BPixHitCount+=1;
             ladder_num = tTopo->pxbLadder(detId);    
             module_num = tTopo->pxbModule(detId);
           }
         }
           }
 
           h_probeL1Ladder_->Fill(ladder_num);
           h_probeL1Module_->Fill(module_num);
           h_probeHasBPixL1Overlap_->Fill(L1BPixHitCount);
 
           // filling the pT-binned distributions
 
           for(int ipTBin=0; ipTBin<nPtBins_; ipTBin++){
         
         float pTF = mypT_bins_[ipTBin];
         float pTL = mypT_bins_[ipTBin+1];
         
         if(debug_)
           edm::LogInfo("PrimaryVertexValidation")<<"ipTBin:"<<ipTBin<< " "<<mypT_bins_[ipTBin]<< " < pT < "<<mypT_bins_[ipTBin+1]<<std::endl;
         
         if( std::abs(tracketa)<1.5 && (trackpt >= pTF && trackpt < pTL) ){
           
           if(debug_)
             edm::LogInfo("PrimaryVertexValidation")<<"passes this cut: "<<mypT_bins_[ipTBin]<<std::endl;
           PVValHelper::fillByIndex(h_dxy_pT_,ipTBin,dxyFromMyVertex*cmToum);
           PVValHelper::fillByIndex(h_dz_pT_,ipTBin,dzFromMyVertex*cmToum);
           PVValHelper::fillByIndex(h_norm_dxy_pT_,ipTBin,dxyFromMyVertex/s_ip2dpv_err);
           PVValHelper::fillByIndex(h_norm_dz_pT_,ipTBin,dzFromMyVertex/dz_err);
           
           if(std::abs(tracketa)<1.){
             
             if(debug_)
               edm::LogInfo("PrimaryVertexValidation")<<"passes tight eta cut: "<<mypT_bins_[ipTBin]<<std::endl;
             PVValHelper::fillByIndex(h_dxy_Central_pT_,ipTBin,dxyFromMyVertex*cmToum);
             PVValHelper::fillByIndex(h_dz_Central_pT_,ipTBin,dzFromMyVertex*cmToum);
             PVValHelper::fillByIndex(h_norm_dxy_Central_pT_,ipTBin,dxyFromMyVertex/s_ip2dpv_err);
             PVValHelper::fillByIndex(h_norm_dz_Central_pT_,ipTBin,dzFromMyVertex/dz_err);
           }
         }
           }
           
           // checks on the probe track quality
           if(trackpt >= ptOfProbe_ 
          && std::abs(tracketa)<= etaOfProbe_ 
          && tracknhits>=nHitsOfProbe_
          && trackp >= pOfProbe_){
 
         std::pair<bool,bool> pixelOcc = pixelHitsCheck((theTTrack));
 
         if(debug_){
           if(pixelOcc.first == true)
             edm::LogInfo("PrimaryVertexValidation")<<"has BPIx hits"<<std::endl;
           if(pixelOcc.second == true)
             edm::LogInfo("PrimaryVertexValidation")<<"has FPix hits"<<std::endl;
         }         
 
         if(!doBPix_ && (pixelOcc.first == true))  continue;
         if(!doFPix_ && (pixelOcc.second == true)) continue;
         
         fillTrackHistos(hDA,"sel",&(theTTrack),vertex,beamSpot,fBfield_);
 
         // probe checks
         h_probePt_->Fill(theTrack.pt());
         h_probePtRebin_->Fill(theTrack.pt());
         h_probeP_->Fill(theTrack.p());
         h_probeEta_->Fill(theTrack.eta());
         h_probePhi_->Fill(theTrack.phi());
         h2_probeEtaPhi_->Fill(theTrack.eta(),theTrack.phi());
         h2_probeEtaPt_->Fill(theTrack.eta(),theTrack.pt());
 
         h_probeChi2_->Fill(theTrack.chi2());
         h_probeNormChi2_->Fill(theTrack.normalizedChi2());
         h_probeCharge_->Fill(theTrack.charge());
         h_probeQoverP_->Fill(theTrack.qoverp());
         h_probeHits_->Fill(theTrack.numberOfValidHits());       
         h_probeHits1D_->Fill(nRecHit1D);
         h_probeHits2D_->Fill(nRecHit2D);
         h_probeHitsInTIB_->Fill(nhitinBPIX);
         h_probeHitsInTOB_->Fill(nhitinFPIX);
         h_probeHitsInTID_->Fill(nhitinTIB);
         h_probeHitsInTEC_->Fill(nhitinTID);
         h_probeHitsInBPIX_->Fill(nhitinTOB);
         h_probeHitsInFPIX_->Fill(nhitinTEC);
         
         float dxyRecoV = theTrack.dz(theRecoVertex);
         float dzRecoV  = theTrack.dxy(theRecoVertex);
         float dxysigmaRecoV = TMath::Sqrt(theTrack.d0Error()*theTrack.d0Error()+xErrOfflineVertex_*yErrOfflineVertex_);
         float dzsigmaRecoV  = TMath::Sqrt(theTrack.dzError()*theTrack.dzError()+zErrOfflineVertex_*zErrOfflineVertex_);
 
         double zTrack=(theTTrack.stateAtBeamLine().trackStateAtPCA()).position().z();
         double zVertex=theFittedVertex.position().z();
         double tantheta=tan((theTTrack.stateAtBeamLine().trackStateAtPCA()).momentum().theta());
 
         double dz2= pow(theTrack.dzError(),2)+wxy2_/pow(tantheta,2);
         double restrkz   = zTrack-zVertex;
         double pulltrkz  = (zTrack-zVertex)/TMath::Sqrt(dz2);
 
         h_probedxyRecoV_->Fill(dxyRecoV);
         h_probedzRecoV_->Fill(dzRecoV);
     
         h_probedzRefitV_->Fill(dxyFromMyVertex);
         h_probedxyRefitV_->Fill(dzFromMyVertex);
 
         h_probed0RefitV_->Fill(d0);
         h_probez0RefitV_->Fill(z0);
 
         h_probesignIP2DRefitV_->Fill(s_ip2dpv_corr);
         h_probed3DRefitV_->Fill(ip3d_corr);
         h_probereszRefitV_->Fill(restrkz);
 
         h_probeRecoVSigZ_->Fill(dzRecoV/dzsigmaRecoV);
         h_probeRecoVSigXY_->Fill(dxyRecoV/dxysigmaRecoV);   
         h_probeRefitVSigZ_->Fill(dzFromMyVertex/dz_err);
         h_probeRefitVSigXY_->Fill(dxyFromMyVertex/s_ip2dpv_err); 
         h_probeRefitVSig3D_->Fill(ip3d_corr/ip3d_err);
         h_probeRefitVLogSig3D_->Fill(log10(ip3d_corr/ip3d_err));
         h_probeRefitVSigResZ_->Fill(pulltrkz);
         
         a_dxyVsPhi->Fill(trackphi,dxyFromMyVertex*cmToum);
         a_dzVsPhi->Fill(trackphi,z0*cmToum);  
         n_dxyVsPhi->Fill(trackphi,dxyFromMyVertex/s_ip2dpv_err); 
         n_dzVsPhi->Fill(trackphi,z0/z0_error);
   
         a_dxyVsEta->Fill(tracketa,dxyFromMyVertex*cmToum); 
         a_dzVsEta->Fill(tracketa,z0*cmToum);  
         n_dxyVsEta->Fill(tracketa,dxyFromMyVertex/s_ip2dpv_err); 
         n_dzVsEta->Fill(tracketa,z0/z0_error); 
 
         if( ladder_num > 0 && module_num > 0 ) {
           
           LogDebug("PrimaryVertexValidation")<<" ladder_num"<<ladder_num <<" module_num"<<module_num <<std::endl;
 
           PVValHelper::fillByIndex(h_dxy_modZ_,module_num-1,dxyFromMyVertex*cmToum);
           PVValHelper::fillByIndex(h_dz_modZ_,module_num-1,dzFromMyVertex*cmToum);    
           PVValHelper::fillByIndex(h_norm_dxy_modZ_,module_num-1,dxyFromMyVertex/s_ip2dpv_err);   
           PVValHelper::fillByIndex(h_norm_dz_modZ_,module_num-1,dzFromMyVertex/dz_err);   
           
           PVValHelper::fillByIndex(h_dxy_ladder_,ladder_num-1,dxyFromMyVertex*cmToum);    
 
           LogDebug("PrimaryVertexValidation")<<"h_dxy_ladder size:" <<h_dxy_ladder_.size() << std::endl;
 
           if(L1BPixHitCount==1){
             PVValHelper::fillByIndex(h_dxy_ladderNoOverlap_,ladder_num-1,dxyFromMyVertex*cmToum);     
           } else {
             PVValHelper::fillByIndex(h_dxy_ladderOverlap_,ladder_num-1,dxyFromMyVertex*cmToum);   
           }
 
           PVValHelper::fillByIndex(h_dz_ladder_,ladder_num-1,dzFromMyVertex*cmToum);      
           PVValHelper::fillByIndex(h_norm_dxy_ladder_,ladder_num-1,dxyFromMyVertex/s_ip2dpv_err);  
           PVValHelper::fillByIndex(h_norm_dz_ladder_,ladder_num-1,dzFromMyVertex/dz_err);   
           
         }
 
         // filling the binned distributions
         for(int i=0; i<nBins_; i++){
         
           float phiF = theDetails_.trendbins[PVValHelper::phi][i];
           float phiL = theDetails_.trendbins[PVValHelper::phi][i+1];
 
           float etaF = theDetails_.trendbins[PVValHelper::eta][i];
           float etaL = theDetails_.trendbins[PVValHelper::eta][i+1];
   
           if(tracketa >= etaF && tracketa < etaL ){
 
             PVValHelper::fillByIndex(a_dxyEtaResiduals,i,dxyFromMyVertex*cmToum,"1");
             PVValHelper::fillByIndex(a_dxEtaResiduals,i,my_dx*cmToum,"2");
             PVValHelper::fillByIndex(a_dyEtaResiduals,i,my_dy*cmToum,"3");
             PVValHelper::fillByIndex(a_dzEtaResiduals,i,dzFromMyVertex*cmToum,"4");   
             PVValHelper::fillByIndex(n_dxyEtaResiduals,i,dxyFromMyVertex/s_ip2dpv_err,"5");
             PVValHelper::fillByIndex(n_dzEtaResiduals,i,dzFromMyVertex/dz_err,"6");     
             PVValHelper::fillByIndex(a_IP2DEtaResiduals,i,s_ip2dpv_corr*cmToum,"7");
             PVValHelper::fillByIndex(n_IP2DEtaResiduals,i,s_ip2dpv_corr/s_ip2dpv_err,"8");
             PVValHelper::fillByIndex(a_reszEtaResiduals,i,restrkz*cmToum,"9");
             PVValHelper::fillByIndex(n_reszEtaResiduals,i,pulltrkz,"10");
             PVValHelper::fillByIndex(a_d3DEtaResiduals,i,ip3d_corr*cmToum,"11");   
             PVValHelper::fillByIndex(n_d3DEtaResiduals,i,ip3d_corr/ip3d_err,"12");
             PVValHelper::fillByIndex(a_IP3DEtaResiduals,i,s_ip3dpv_corr*cmToum,"13");
             PVValHelper::fillByIndex(n_IP3DEtaResiduals,i,s_ip3dpv_corr/s_ip3dpv_err,"14");
 
           }
               
           if(trackphi >= phiF && trackphi < phiL ){
 
             PVValHelper::fillByIndex(a_dxyPhiResiduals,i,dxyFromMyVertex*cmToum,"15");
             PVValHelper::fillByIndex(a_dxPhiResiduals,i,my_dx*cmToum,"16");
             PVValHelper::fillByIndex(a_dyPhiResiduals,i,my_dy*cmToum,"17");
             PVValHelper::fillByIndex(a_dzPhiResiduals,i,dzFromMyVertex*cmToum,"18"); 
             PVValHelper::fillByIndex(n_dxyPhiResiduals,i,dxyFromMyVertex/s_ip2dpv_err,"19");
             PVValHelper::fillByIndex(n_dzPhiResiduals,i,dzFromMyVertex/dz_err,"20"); 
             PVValHelper::fillByIndex(a_IP2DPhiResiduals,i,s_ip2dpv_corr*cmToum,"21");
             PVValHelper::fillByIndex(n_IP2DPhiResiduals,i,s_ip2dpv_corr/s_ip2dpv_err,"22"); 
             PVValHelper::fillByIndex(a_reszPhiResiduals,i,restrkz*cmToum,"23");
             PVValHelper::fillByIndex(n_reszPhiResiduals,i,pulltrkz,"24");
             PVValHelper::fillByIndex(a_d3DPhiResiduals,i,ip3d_corr*cmToum,"25");   
             PVValHelper::fillByIndex(n_d3DPhiResiduals,i,ip3d_corr/ip3d_err,"26");
             PVValHelper::fillByIndex(a_IP3DPhiResiduals,i,s_ip3dpv_corr*cmToum,"27");
             PVValHelper::fillByIndex(n_IP3DPhiResiduals,i,s_ip3dpv_corr/s_ip3dpv_err,"28");
 
             for(int j=0; j<nBins_; j++){
 
               float etaJ = theDetails_.trendbins[PVValHelper::eta][j];
               float etaK = theDetails_.trendbins[PVValHelper::eta][j+1];
 
               if(tracketa >= etaJ && tracketa < etaK ){
             a_dxyResidualsMap[i][j]->Fill(dxyFromMyVertex*cmToum); 
             a_dzResidualsMap[i][j]->Fill(dzFromMyVertex*cmToum);        
             n_dxyResidualsMap[i][j]->Fill(dxyFromMyVertex/s_ip2dpv_err); 
             n_dzResidualsMap[i][j]->Fill(dzFromMyVertex/dz_err);  
             a_d3DResidualsMap[i][j]->Fill(ip3d_corr*cmToum);   
             n_d3DResidualsMap[i][j]->Fill(ip3d_corr/ip3d_err); 
 
               }
             }
           }     
         }
           }
               
           if(debug_){
         edm::LogInfo("PrimaryVertexValidation")<<" myVertex.x()= "<<myVertex.x()<<"\n"
                                <<" myVertex.y()= "<<myVertex.y()<<" \n"
                                <<" myVertex.z()= "<<myVertex.z()<<" \n"
                                <<" theTrack.dz(myVertex)= "<<theTrack.dz(myVertex)<<" \n"
                                <<" zPCA -myVertex.z() = "<<(theTrack.vertex().z() -myVertex.z());
         
           }// ends if debug_
         } // ends if the fitted vertex is Valid
 
         //delete theFitter;
 
       }  catch ( cms::Exception& er ) {
         LogTrace("PrimaryVertexValidation")<<"caught std::exception "<<er.what()<<std::endl;
       }
         
     } //ends if theFinalTracks.size() > 2
     
     else {
       if(debug_)
         edm::LogInfo("PrimaryVertexValidation")<<"Not enough tracks to make a vertex.  Returns no vertex info";
     }
       
     ++nTracks_;  
     ++nTracksPerClus_;
 
     if(debug_)
       edm::LogInfo("PrimaryVertexValidation")<<"Track "<<i<<" : pT = "<<theTrack.pt();
       
       }// for loop on tracks
   } // for loop on track clusters
   
   // Fill the TTree if needed
 
   if(storeNtuple_){  
     rootTree_->Fill();
   }
 
   
 } 

void PrimaryVertexValidation::beginJob ( void )

overrideprivatevirtual

Reimplemented from edm::one::EDAnalyzerBase.

Definition at line 1156 of file PrimaryVertexValidation.cc.

References a_d3DEtaResiduals, a_d3DPhiResiduals, a_d3DResidualsMap, a_dxEtaResiduals, a_dxPhiResiduals, a_dxyBiasResidualsMap, a_dxyEtaBiasResiduals, a_dxyEtaMADBiasTrend, a_dxyEtaMADTrend, a_dxyEtaMeanBiasTrend, a_dxyEtaMeanTrend, a_dxyEtaMedianBiasTrend, a_dxyEtaMedianTrend, a_dxyEtaResiduals, a_dxyEtaWidthBiasTrend, a_dxyEtaWidthTrend, a_dxyMeanBiasMap, a_dxyMeanMap, a_dxyPhiBiasResiduals, a_dxyPhiMADBiasTrend, a_dxyPhiMADTrend, a_dxyPhiMeanBiasTrend, a_dxyPhiMeanTrend, a_dxyPhiMedianBiasTrend, a_dxyPhiMedianTrend, a_dxyPhiResiduals, a_dxyPhiWidthBiasTrend, a_dxyPhiWidthTrend, a_dxypTCentralMeanTrend, a_dxypTCentralWidthTrend, a_dxypTMeanTrend, a_dxypTWidthTrend, a_dxyResidualsMap, a_dxyVsEta, a_dxyVsPhi, a_dxyWidthBiasMap, a_dxyWidthMap, a_dyEtaResiduals, a_dyPhiResiduals, a_dzBiasResidualsMap, a_dzEtaBiasResiduals, a_dzEtaMADBiasTrend, a_dzEtaMADTrend, a_dzEtaMeanBiasTrend, a_dzEtaMeanTrend, a_dzEtaMedianBiasTrend, a_dzEtaMedianTrend, a_dzEtaResiduals, a_dzEtaWidthBiasTrend, a_dzEtaWidthTrend, a_dzMeanBiasMap, a_dzMeanMap, a_dzPhiBiasResiduals, a_dzPhiMADBiasTrend, a_dzPhiMADTrend, a_dzPhiMeanBiasTrend, a_dzPhiMeanTrend, a_dzPhiMedianBiasTrend, a_dzPhiMedianTrend, a_dzPhiResiduals, a_dzPhiWidthBiasTrend, a_dzPhiWidthTrend, a_dzpTCentralMeanTrend, a_dzpTCentralWidthTrend, a_dzpTMeanTrend, a_dzpTWidthTrend, a_dzResidualsMap, a_dzVsEta, a_dzVsPhi, a_dzWidthBiasMap, a_dzWidthMap, a_IP2DEtaResiduals, a_IP2DPhiResiduals, a_IP3DEtaResiduals, a_IP3DPhiResiduals, a_reszEtaResiduals, a_reszPhiResiduals, Beamdxdz_, Beamsigmaz_, BeamWidthX_, BeamWidthY_, bookResidualsHistogram(), bookVertexHistograms(), BSx0_, BSy0_, BSz0_, charge_, chi2_, chi2ndof_, chi2normUnbiasedVertex_, chi2ProbUnbiasedVertex_, chi2UnbiasedVertex_, PVValHelper::d3D, d3DFromMyVertex_, DOFUnbiasedVertex_, PVValHelper::dx, PVValHelper::dxy, dxy_, dxyBs_, dxyErrorFromMyVertex_, dxyFromMyVertex_, PVValHelper::dy, PVValHelper::dz, dz_, dzBs_, dzErrorFromMyVertex_, dzFromMyVertex_, PVValHelper::eta, eta_, etaOfProbe_, EventNumber_, fs, PVValHelper::histodetails::getHigh(), h2_probeEtaPhi_, h2_probeEtaPt_, h_Beamsigmaz, h_BeamWidthX, h_BeamWidthY, h_BSx0, h_BSy0, h_BSz0, h_dxy_Central_pT_, h_dxy_ladder_, h_dxy_ladderNoOverlap_, h_dxy_ladderOverlap_, h_dxy_modZ_, h_dxy_pT_, h_dz_Central_pT_, h_dz_ladder_, h_dz_modZ_, h_dz_pT_, h_etaMax, h_fitVtxChi2_, h_fitVtxChi2ndf_, h_fitVtxChi2Prob_, h_fitVtxNdof_, h_fitVtxNtracks_, h_fitVtxTrackAverageWeight_, h_fitVtxTrackWeights_, h_lumiFromConfig, h_nbins, h_nClus, h_nLadders, h_nOfflineVertices, h_norm_dxy_Central_pT_, h_norm_dxy_ladder_, h_norm_dxy_modZ_, h_norm_dxy_pT_, h_norm_dz_Central_pT_, h_norm_dz_ladder_, h_norm_dz_modZ_, h_norm_dz_pT_, h_nTracks, h_probeCharge_, h_probeChi2_, h_probed0RefitV_, h_probed3DRefitV_, h_probedxyRecoV_, h_probedxyRefitV_, h_probedzRecoV_, h_probedzRefitV_, h_probeEta_, h_probeHasBPixL1Overlap_, h_probeHits1D_, h_probeHits2D_, h_probeHits_, h_probeHitsInBPIX_, h_probeHitsInFPIX_, h_probeHitsInTEC_, h_probeHitsInTIB_, h_probeHitsInTID_, h_probeHitsInTOB_, h_probeL1Ladder_, h_probeL1Module_, h_probeNormChi2_, h_probeP_, h_probePhi_, h_probePt_, h_probePtRebin_, h_probeQoverP_, h_probeRecoVSigXY_, h_probeRecoVSigZ_, h_probeRefitVLogSig3D_, h_probeRefitVSig3D_, h_probeRefitVSigResZ_, h_probeRefitVSigXY_, h_probeRefitVSigZ_, h_probereszRefitV_, h_probesignIP2DRefitV_, h_probez0RefitV_, h_pTinfo, h_recoVtxChi2ndf_, h_recoVtxChi2Prob_, h_recoVtxNtracks_, h_recoVtxSumPt_, h_runFromConfig, h_runFromEvent, h_runNumber, h_xErrOfflineVertex, h_xOfflineVertex, h_yErrOfflineVertex, h_yOfflineVertex, h_zErrOfflineVertex, h_zOfflineVertex, hasRecVertex_, hDA, PVValHelper::histodetails::histobins, mps_fire::i, intLumi_, PVValHelper::IP2D, PVValHelper::IP3D, IP3DsigFromMyVertex_, IPLsigFromMyVertex_, IPTsigFromMyVertex_, isGoodTrack_, isHighPurity_, PVValHelper::ladder, lightNtupleSwitch_, LuminosityBlockNumber_, M_PI, MADTrendsDir, TFileDirectory::make(), TFileService::make(), Mean2DMapsDir, MeanTrendsDir, MedianTrendsDir, TFileService::mkdir(), PVValHelper::modZ, mypT_bins_, n_d3DEtaResiduals, n_d3DPhiResiduals, n_d3DResidualsMap, n_dxyBiasResidualsMap, n_dxyEtaBiasResiduals, n_dxyEtaMADBiasTrend, n_dxyEtaMADTrend, n_dxyEtaMeanBiasTrend, n_dxyEtaMeanTrend, n_dxyEtaMedianBiasTrend, n_dxyEtaMedianTrend, n_dxyEtaResiduals, n_dxyEtaWidthBiasTrend, n_dxyEtaWidthTrend, n_dxyMeanBiasMap, n_dxyMeanMap, n_dxyPhiBiasResiduals, n_dxyPhiMADBiasTrend, n_dxyPhiMADTrend, n_dxyPhiMeanBiasTrend, n_dxyPhiMeanTrend, n_dxyPhiMedianBiasTrend, n_dxyPhiMedianTrend, n_dxyPhiResiduals, n_dxyPhiWidthBiasTrend, n_dxyPhiWidthTrend, n_dxypTCentralMeanTrend, n_dxypTCentralWidthTrend, n_dxypTMeanTrend, n_dxypTWidthTrend, n_dxyResidualsMap, n_dxyVsEta, n_dxyVsPhi, n_dxyWidthBiasMap, n_dxyWidthMap, n_dzBiasResidualsMap, n_dzEtaBiasResiduals, n_dzEtaMADBiasTrend, n_dzEtaMADTrend, n_dzEtaMeanBiasTrend, n_dzEtaMeanTrend, n_dzEtaMedianBiasTrend, n_dzEtaMedianTrend, n_dzEtaResiduals, n_dzEtaWidthBiasTrend, n_dzEtaWidthTrend, n_dzMeanBiasMap, n_dzMeanMap, n_dzPhiBiasResiduals, n_dzPhiMADBiasTrend, n_dzPhiMADTrend, n_dzPhiMeanBiasTrend, n_dzPhiMeanTrend, n_dzPhiMedianBiasTrend, n_dzPhiMedianTrend, n_dzPhiResiduals, n_dzPhiWidthBiasTrend, n_dzPhiWidthTrend, n_dzpTCentralMeanTrend, n_dzpTCentralWidthTrend, n_dzpTMeanTrend, n_dzpTWidthTrend, n_dzResidualsMap, n_dzVsEta, n_dzVsPhi, n_dzWidthBiasMap, n_dzWidthMap, n_IP2DEtaResiduals, n_IP2DPhiResiduals, n_IP3DEtaResiduals, n_IP3DPhiResiduals, n_reszEtaResiduals, n_reszPhiResiduals, nBins_, nClus_, Nevt_, nhits1D_, nhits2D_, nhits_, nhitsBPIX_, nhitsFPIX_, nhitsTEC_, nhitsTIB_, nhitsTID_, nhitsTOB_, nLadders_, nOfflineVertices_, PVValHelper::norm_d3D, PVValHelper::norm_dxy, PVValHelper::norm_dz, PVValHelper::norm_IP2D, PVValHelper::norm_IP3D, PVValHelper::norm_resz, nPtBins_, nTracks_, nTracksPerClus_, p_, PVValHelper::phi, phi_, PVValHelper::pT, pt_, PVValHelper::pTCentral, qoverp_, alignCSCRings::r, PVValHelper::resz, rootTree_, runControlNumbers_, RunNumber_, sumOfWeightsUnbiasedVertex_, theDetails_, theta_, tracksUsedForVertexing_, PVValHelper::histodetails::trendbins, useTracksFromRecoVtx_, Width2DMapsDir, WidthTrendsDir, xOfflineVertex_, xPCA_, xUnbiasedVertex_, yOfflineVertex_, yPCA_, yUnbiasedVertex_, zOfflineVertex_, zPCA_, and zUnbiasedVertex_.

 {
   edm::LogInfo("PrimaryVertexValidation") 
     <<"######################################\n"
     <<"Begin Job \n" 
     <<"######################################";
 
   // Define TTree for output
   Nevt_    = 0;
   
   //  rootFile_ = new TFile(filename_.c_str(),"recreate");
   rootTree_ = fs->make<TTree>("tree","PV Validation tree");
  
   // Track Paramters 
 
   if(lightNtupleSwitch_){ 
 
     rootTree_->Branch("EventNumber",&EventNumber_,"EventNumber/i");
     rootTree_->Branch("RunNumber",&RunNumber_,"RunNumber/i");
     rootTree_->Branch("LuminosityBlockNumber",&LuminosityBlockNumber_,"LuminosityBlockNumber/i");
     rootTree_->Branch("nOfflineVertices",&nOfflineVertices_,"nOfflineVertices/I");
     rootTree_->Branch("nTracks",&nTracks_,"nTracks/I");
     rootTree_->Branch("phi",&phi_,"phi[nTracks]/D");
     rootTree_->Branch("eta",&eta_,"eta[nTracks]/D");
     rootTree_->Branch("pt",&pt_,"pt[nTracks]/D");
     rootTree_->Branch("dxyFromMyVertex",&dxyFromMyVertex_,"dxyFromMyVertex[nTracks]/D");
     rootTree_->Branch("dzFromMyVertex",&dzFromMyVertex_,"dzFromMyVertex[nTracks]/D"); 
     rootTree_->Branch("d3DFromMyVertex",&d3DFromMyVertex_,"d3DFromMyVertex[nTracks]/D"); 
     rootTree_->Branch("IPTsigFromMyVertex",&IPTsigFromMyVertex_,"IPTsigFromMyVertex_[nTracks]/D");    
     rootTree_->Branch("IPLsigFromMyVertex",&IPLsigFromMyVertex_,"IPLsigFromMyVertex_[nTracks]/D"); 
     rootTree_->Branch("IP3DsigFromMyVertex",&IP3DsigFromMyVertex_,"IP3DsigFromMyVertex_[nTracks]/D"); 
     rootTree_->Branch("hasRecVertex",&hasRecVertex_,"hasRecVertex[nTracks]/I");
     rootTree_->Branch("isGoodTrack",&isGoodTrack_,"isGoodTrack[nTracks]/I");
     rootTree_->Branch("isHighPurity",&isHighPurity_,"isHighPurity_[nTracks]/I");
     
   } else {
     
     rootTree_->Branch("nTracks",&nTracks_,"nTracks/I");
     rootTree_->Branch("nTracksPerClus",&nTracksPerClus_,"nTracksPerClus/I");
     rootTree_->Branch("nClus",&nClus_,"nClus/I");
     rootTree_->Branch("xOfflineVertex",&xOfflineVertex_,"xOfflineVertex/D");
     rootTree_->Branch("yOfflineVertex",&yOfflineVertex_,"yOfflineVertex/D");
     rootTree_->Branch("zOfflineVertex",&zOfflineVertex_,"zOfflineVertex/D");
     rootTree_->Branch("BSx0",&BSx0_,"BSx0/D");
     rootTree_->Branch("BSy0",&BSy0_,"BSy0/D");
     rootTree_->Branch("BSz0",&BSz0_,"BSz0/D");
     rootTree_->Branch("Beamsigmaz",&Beamsigmaz_,"Beamsigmaz/D");
     rootTree_->Branch("Beamdxdz",&Beamdxdz_,"Beamdxdz/D");
     rootTree_->Branch("BeamWidthX",&BeamWidthX_,"BeamWidthX/D");
     rootTree_->Branch("BeamWidthY",&BeamWidthY_,"BeamWidthY/D");
     rootTree_->Branch("pt",&pt_,"pt[nTracks]/D");
     rootTree_->Branch("p",&p_,"p[nTracks]/D");
     rootTree_->Branch("nhits",&nhits_,"nhits[nTracks]/I");
     rootTree_->Branch("nhits1D",&nhits1D_,"nhits1D[nTracks]/I");
     rootTree_->Branch("nhits2D",&nhits2D_,"nhits2D[nTracks]/I");
     rootTree_->Branch("nhitsBPIX",&nhitsBPIX_,"nhitsBPIX[nTracks]/I");
     rootTree_->Branch("nhitsFPIX",&nhitsFPIX_,"nhitsFPIX[nTracks]/I");
     rootTree_->Branch("nhitsTIB",&nhitsTIB_,"nhitsTIB[nTracks]/I");
     rootTree_->Branch("nhitsTID",&nhitsTID_,"nhitsTID[nTracks]/I");
     rootTree_->Branch("nhitsTOB",&nhitsTOB_,"nhitsTOB[nTracks]/I");
     rootTree_->Branch("nhitsTEC",&nhitsTEC_,"nhitsTEC[nTracks]/I");
     rootTree_->Branch("eta",&eta_,"eta[nTracks]/D");
     rootTree_->Branch("theta",&theta_,"theta[nTracks]/D");
     rootTree_->Branch("phi",&phi_,"phi[nTracks]/D");
     rootTree_->Branch("chi2",&chi2_,"chi2[nTracks]/D");
     rootTree_->Branch("chi2ndof",&chi2ndof_,"chi2ndof[nTracks]/D");
     rootTree_->Branch("charge",&charge_,"charge[nTracks]/I");
     rootTree_->Branch("qoverp",&qoverp_,"qoverp[nTracks]/D");
     rootTree_->Branch("dz",&dz_,"dz[nTracks]/D");
     rootTree_->Branch("dxy",&dxy_,"dxy[nTracks]/D");
     rootTree_->Branch("dzBs",&dzBs_,"dzBs[nTracks]/D");
     rootTree_->Branch("dxyBs",&dxyBs_,"dxyBs[nTracks]/D");
     rootTree_->Branch("xPCA",&xPCA_,"xPCA[nTracks]/D");
     rootTree_->Branch("yPCA",&yPCA_,"yPCA[nTracks]/D");
     rootTree_->Branch("zPCA",&zPCA_,"zPCA[nTracks]/D");
     rootTree_->Branch("xUnbiasedVertex",&xUnbiasedVertex_,"xUnbiasedVertex[nTracks]/D");
     rootTree_->Branch("yUnbiasedVertex",&yUnbiasedVertex_,"yUnbiasedVertex[nTracks]/D");
     rootTree_->Branch("zUnbiasedVertex",&zUnbiasedVertex_,"zUnbiasedVertex[nTracks]/D");
     rootTree_->Branch("chi2normUnbiasedVertex",&chi2normUnbiasedVertex_,"chi2normUnbiasedVertex[nTracks]/F");
     rootTree_->Branch("chi2UnbiasedVertex",&chi2UnbiasedVertex_,"chi2UnbiasedVertex[nTracks]/F");
     rootTree_->Branch("DOFUnbiasedVertex",&DOFUnbiasedVertex_," DOFUnbiasedVertex[nTracks]/F");
     rootTree_->Branch("chi2ProbUnbiasedVertex",&chi2ProbUnbiasedVertex_,"chi2ProbUnbiasedVertex[nTracks]/F");
     rootTree_->Branch("sumOfWeightsUnbiasedVertex",&sumOfWeightsUnbiasedVertex_,"sumOfWeightsUnbiasedVertex[nTracks]/F");
     rootTree_->Branch("tracksUsedForVertexing",&tracksUsedForVertexing_,"tracksUsedForVertexing[nTracks]/I");
     rootTree_->Branch("dxyFromMyVertex",&dxyFromMyVertex_,"dxyFromMyVertex[nTracks]/D");
     rootTree_->Branch("dzFromMyVertex",&dzFromMyVertex_,"dzFromMyVertex[nTracks]/D");
     rootTree_->Branch("dxyErrorFromMyVertex",&dxyErrorFromMyVertex_,"dxyErrorFromMyVertex_[nTracks]/D"); 
     rootTree_->Branch("dzErrorFromMyVertex",&dzErrorFromMyVertex_,"dzErrorFromMyVertex_[nTracks]/D");  
     rootTree_->Branch("IPTsigFromMyVertex",&IPTsigFromMyVertex_,"IPTsigFromMyVertex_[nTracks]/D");    
     rootTree_->Branch("IPLsigFromMyVertex",&IPLsigFromMyVertex_,"IPLsigFromMyVertex_[nTracks]/D"); 
     rootTree_->Branch("hasRecVertex",&hasRecVertex_,"hasRecVertex[nTracks]/I");
     rootTree_->Branch("isGoodTrack",&isGoodTrack_,"isGoodTrack[nTracks]/I");   
   }
 
   // event histograms
   TFileDirectory EventFeatures = fs->mkdir("EventFeatures");
 
   TH1F::SetDefaultSumw2(kTRUE);
 
   h_lumiFromConfig    = EventFeatures.make<TH1F>("h_lumiFromConfig","luminosity from config;;luminosity of present run",1,-0.5,0.5);
   h_lumiFromConfig->SetBinContent(1,intLumi_);
 
   h_runFromConfig     = EventFeatures.make<TH1I>("h_runFromConfig","run number from config;;run number (from configuration)",
                          runControlNumbers_.size(),0.,runControlNumbers_.size());
   for(const auto r : runControlNumbers_){
     h_runFromConfig->SetBinContent(r+1, r);
   }
   
   h_runFromEvent      = EventFeatures.make<TH1I>("h_runFromEvent","run number from config;;run number (from event)",1,-0.5,0.5);
   h_nTracks           = EventFeatures.make<TH1F>("h_nTracks","number of tracks per event;n_{tracks}/event;n_{events}",300,-0.5,299.5);       
   h_nClus             = EventFeatures.make<TH1F>("h_nClus","number of track clusters;n_{clusters}/event;n_{events}",50,-0.5,49.5);       
   h_nOfflineVertices  = EventFeatures.make<TH1F>("h_nOfflineVertices","number of offline reconstructed vertices;n_{vertices}/event;n_{events}",50,-0.5,49.5);  
   h_runNumber         = EventFeatures.make<TH1F>("h_runNumber","run number;run number;n_{events}",100000,250000.,350000.);       
   h_xOfflineVertex    = EventFeatures.make<TH1F>("h_xOfflineVertex","x-coordinate of offline vertex;x_{vertex};n_{events}",100,-0.1,0.1);    
   h_yOfflineVertex    = EventFeatures.make<TH1F>("h_yOfflineVertex","y-coordinate of offline vertex;y_{vertex};n_{events}",100,-0.1,0.1);    
   h_zOfflineVertex    = EventFeatures.make<TH1F>("h_zOfflineVertex","z-coordinate of offline vertex;z_{vertex};n_{events}",100,-30.,30.);    
   h_xErrOfflineVertex = EventFeatures.make<TH1F>("h_xErrOfflineVertex","x-coordinate error of offline vertex;err_{x}^{vtx};n_{events}",100,0.,0.01); 
   h_yErrOfflineVertex = EventFeatures.make<TH1F>("h_yErrOfflineVertex","y-coordinate error of offline vertex;err_{y}^{vtx};n_{events}",100,0.,0.01); 
   h_zErrOfflineVertex = EventFeatures.make<TH1F>("h_zErrOfflineVertex","z-coordinate error of offline vertex;err_{z}^{vtx};n_{events}",100,0.,10.); 
   h_BSx0              = EventFeatures.make<TH1F>("h_BSx0","x-coordinate of reco beamspot;x^{BS}_{0};n_{events}",100,-0.1,0.1);           
   h_BSy0              = EventFeatures.make<TH1F>("h_BSy0","y-coordinate of reco beamspot;y^{BS}_{0};n_{events}",100,-0.1,0.1);           
   h_BSz0              = EventFeatures.make<TH1F>("h_BSz0","z-coordinate of reco beamspot;z^{BS}_{0};n_{events}",100,-1.,1.);             
   h_Beamsigmaz        = EventFeatures.make<TH1F>("h_Beamsigmaz","z-coordinate beam width;#sigma_{Z}^{beam};n_{events}",100,0.,1.);               
   h_BeamWidthX        = EventFeatures.make<TH1F>("h_BeamWidthX","x-coordinate beam width;#sigma_{X}^{beam};n_{events}",100,0.,0.01);         
   h_BeamWidthY        = EventFeatures.make<TH1F>("h_BeamWidthY","y-coordinate beam width;#sigma_{Y}^{beam};n_{events}",100,0.,0.01);        
 
   h_etaMax            = EventFeatures.make<TH1F>("etaMax","etaMax",1,-0.5,0.5);
   h_pTinfo            = EventFeatures.make<TH1F>("pTinfo","pTinfo",3,-1.5,1.5);
   h_pTinfo->GetXaxis()->SetBinLabel(1,"n. bins");
   h_pTinfo->GetXaxis()->SetBinLabel(2,"pT min");
   h_pTinfo->GetXaxis()->SetBinLabel(3,"pT max");
 
   h_nbins             = EventFeatures.make<TH1F>("nbins","nbins",1,-0.5,0.5);
   h_nLadders          = EventFeatures.make<TH1F>("nladders","n. ladders",1,-0.5,0.5);
 
   // probe track histograms
   TFileDirectory ProbeFeatures = fs->mkdir("ProbeTrackFeatures");
 
   h_probePt_         = ProbeFeatures.make<TH1F>("h_probePt","p_{T} of probe track;track p_{T} (GeV); tracks",100,0.,50.);   
   h_probePtRebin_    = ProbeFeatures.make<TH1F>("h_probePtRebin","p_{T} of probe track;track p_{T} (GeV); tracks",mypT_bins_.size()-1,mypT_bins_.data());   
   h_probeP_          = ProbeFeatures.make<TH1F>("h_probeP","momentum of probe track;track p (GeV); tracks",100,0.,100.);   
   h_probeEta_        = ProbeFeatures.make<TH1F>("h_probeEta","#eta of the probe track;track #eta;tracks",54,-2.8,2.8);  
   h_probePhi_        = ProbeFeatures.make<TH1F>("h_probePhi","#phi of probe track;track #phi (rad);tracks",100,-3.15,3.15);  
 
   h2_probeEtaPhi_    = ProbeFeatures.make<TH2F>("h2_probeEtaPhi","probe track #phi vs #eta;#eta of probe track;track #phi of probe track (rad); tracks",54,-2.8,2.8,100,-3.15,3.15);  
   h2_probeEtaPt_     = ProbeFeatures.make<TH2F>("h2_probeEtaPt","probe track p_{T} vs #eta;#eta of probe track;track p_{T} (GeV); tracks",54,-2.8,2.8,100,0.,50.);    
 
   h_probeChi2_       = ProbeFeatures.make<TH1F>("h_probeChi2","#chi^{2} of probe track;track #chi^{2}; tracks",100,0.,100.); 
   h_probeNormChi2_   = ProbeFeatures.make<TH1F>("h_probeNormChi2"," normalized #chi^{2} of probe track;track #chi^{2}/ndof; tracks",100,0.,10.);
   h_probeCharge_     = ProbeFeatures.make<TH1F>("h_probeCharge","charge of profe track;track charge Q;tracks",3,-1.5,1.5);
   h_probeQoverP_     = ProbeFeatures.make<TH1F>("h_probeQoverP","q/p of probe track; track Q/p (GeV^{-1});tracks",200,-1.,1.);
   h_probedzRecoV_    = ProbeFeatures.make<TH1F>("h_probedzRecoV","d_{z}(V_{offline}) of probe track;track d_{z}(V_{off}) (cm);tracks",200,-1.,1.);  
   h_probedxyRecoV_   = ProbeFeatures.make<TH1F>("h_probedxyRecoV","d_{xy}(V_{offline}) of probe track;track d_{xy}(V_{off}) (cm);tracks",200,-1.,1.);      
   h_probedzRefitV_   = ProbeFeatures.make<TH1F>("h_probedzRefitV","d_{z}(V_{refit}) of probe track;track d_{z}(V_{fit}) (cm);tracks",200,-0.5,0.5);
   h_probesignIP2DRefitV_ = ProbeFeatures.make<TH1F>("h_probesignIPRefitV","ip_{2D}(V_{refit}) of probe track;track ip_{2D}(V_{fit}) (cm);tracks",200,-1.,1.);
   h_probedxyRefitV_  = ProbeFeatures.make<TH1F>("h_probedxyRefitV","d_{xy}(V_{refit}) of probe track;track d_{xy}(V_{fit}) (cm);tracks",200,-0.5,0.5); 
 
   h_probez0RefitV_   = ProbeFeatures.make<TH1F>("h_probez0RefitV","z_{0}(V_{refit}) of probe track;track z_{0}(V_{fit}) (cm);tracks",200,-1.,1.);
   h_probed0RefitV_   = ProbeFeatures.make<TH1F>("h_probed0RefitV","d_{0}(V_{refit}) of probe track;track d_{0}(V_{fit}) (cm);tracks",200,-1.,1.);
 
   h_probed3DRefitV_  = ProbeFeatures.make<TH1F>("h_probed3DRefitV","d_{3D}(V_{refit}) of probe track;track d_{3D}(V_{fit}) (cm);tracks",200,0.,1.); 
   h_probereszRefitV_ = ProbeFeatures.make<TH1F>("h_probeReszRefitV","z_{track} -z_{V_{refit}};track res_{z}(V_{refit}) (cm);tracks",200,-1.,1.); 
 
   h_probeRecoVSigZ_  = ProbeFeatures.make<TH1F>("h_probeRecoVSigZ"  ,"Longitudinal DCA Significance (reco);d_{z}(V_{off})/#sigma_{dz};tracks",100,-8,8);
   h_probeRecoVSigXY_ = ProbeFeatures.make<TH1F>("h_probeRecoVSigXY" ,"Transverse DCA Significance (reco);d_{xy}(V_{off})/#sigma_{dxy};tracks",100,-8,8);
   h_probeRefitVSigZ_ = ProbeFeatures.make<TH1F>("h_probeRefitVSigZ" ,"Longitudinal DCA Significance (refit);d_{z}(V_{fit})/#sigma_{dz};tracks",100,-8,8);
   h_probeRefitVSigXY_= ProbeFeatures.make<TH1F>("h_probeRefitVSigXY","Transverse DCA Significance (refit);d_{xy}(V_{fit})/#sigma_{dxy};tracks",100,-8,8);
   h_probeRefitVSig3D_= ProbeFeatures.make<TH1F>("h_probeRefitVSig3D","3D DCA Significance (refit);d_{3D}/#sigma_{3D};tracks",100,0.,20.); 
   h_probeRefitVLogSig3D_ = ProbeFeatures.make<TH1F>("h_probeRefitVLogSig3D","log_{10}(3D DCA-Significance) (refit);log_{10}(d_{3D}/#sigma_{3D});tracks",100,-5.,4.); 
   h_probeRefitVSigResZ_ = ProbeFeatures.make<TH1F>("h_probeRefitVSigResZ" ,"Longitudinal residual significance (refit);(z_{track} -z_{V_{fit}})/#sigma_{res_{z}};tracks",100,-8,8);
 
   h_probeHits_       = ProbeFeatures.make<TH1F>("h_probeNRechits"    ,"N_{hits}     ;N_{hits}    ;tracks",40,-0.5,39.5);
   h_probeHits1D_     = ProbeFeatures.make<TH1F>("h_probeNRechits1D"  ,"N_{hits} 1D  ;N_{hits} 1D ;tracks",40,-0.5,39.5);
   h_probeHits2D_     = ProbeFeatures.make<TH1F>("h_probeNRechits2D"  ,"N_{hits} 2D  ;N_{hits} 2D ;tracks",40,-0.5,39.5);
   h_probeHitsInTIB_  = ProbeFeatures.make<TH1F>("h_probeNRechitsTIB" ,"N_{hits} TIB ;N_{hits} TIB;tracks",40,-0.5,39.5);
   h_probeHitsInTOB_  = ProbeFeatures.make<TH1F>("h_probeNRechitsTOB" ,"N_{hits} TOB ;N_{hits} TOB;tracks",40,-0.5,39.5);
   h_probeHitsInTID_  = ProbeFeatures.make<TH1F>("h_probeNRechitsTID" ,"N_{hits} TID ;N_{hits} TID;tracks",40,-0.5,39.5);
   h_probeHitsInTEC_  = ProbeFeatures.make<TH1F>("h_probeNRechitsTEC" ,"N_{hits} TEC ;N_{hits} TEC;tracks",40,-0.5,39.5);
   h_probeHitsInBPIX_ = ProbeFeatures.make<TH1F>("h_probeNRechitsBPIX","N_{hits} BPIX;N_{hits} BPIX;tracks",40,-0.5,39.5);
   h_probeHitsInFPIX_ = ProbeFeatures.make<TH1F>("h_probeNRechitsFPIX","N_{hits} FPIX;N_{hits} FPIX;tracks",40,-0.5,39.5);
 
   h_probeL1Ladder_         = ProbeFeatures.make<TH1F>("h_probeL1Ladder","Ladder number (L1 hit); ladder number",22,-1.5,20.5); 
   h_probeL1Module_         = ProbeFeatures.make<TH1F>("h_probeL1Module","Module number (L1 hit); module number",10,-1.5,8.5);
   h_probeHasBPixL1Overlap_ = ProbeFeatures.make<TH1I>("h_probeHasBPixL1Overlap","n. hits in L1;n. L1-BPix hits;tracks",5,0,5);
 
   // refit vertex features
   TFileDirectory RefitVertexFeatures = fs->mkdir("RefitVertexFeatures");
   h_fitVtxNtracks_          = RefitVertexFeatures.make<TH1F>("h_fitVtxNtracks"  ,"N_{trks} used in vertex fit;N^{fit}_{tracks};vertices"        ,100,-0.5,99.5);
   h_fitVtxNdof_             = RefitVertexFeatures.make<TH1F>("h_fitVtxNdof"     ,"N_{DOF} of vertex fit;N_{DOF} of refit vertex;vertices"          ,100,-0.5,99.5);
   h_fitVtxChi2_             = RefitVertexFeatures.make<TH1F>("h_fitVtxChi2"     ,"#chi^{2} of vertex fit;vertex #chi^{2};vertices"            ,100,-0.5,99.5);
   h_fitVtxChi2ndf_          = RefitVertexFeatures.make<TH1F>("h_fitVtxChi2ndf"  ,"#chi^{2}/ndf of vertex fit;vertex #chi^{2}/ndf;vertices"    ,100,-0.5,9.5);
   h_fitVtxChi2Prob_         = RefitVertexFeatures.make<TH1F>("h_fitVtxChi2Prob" ,"Prob(#chi^{2},ndf) of vertex fit;Prob(#chi^{2},ndf);vertices",40,0.,1.);
   h_fitVtxTrackWeights_     = RefitVertexFeatures.make<TH1F>("h_fitVtxTrackWeights","track weights associated to track;track weights;tracks",40,0.,1.);
   h_fitVtxTrackAverageWeight_ = RefitVertexFeatures.make<TH1F>("h_fitVtxTrackAverageWeight_","average track weight per vertex;#LT track weight #GT;vertices",40,0.,1.);
 
   if(useTracksFromRecoVtx_) {
 
     TFileDirectory RecoVertexFeatures = fs->mkdir("RecoVertexFeatures");
     h_recoVtxNtracks_          = RecoVertexFeatures.make<TH1F>("h_recoVtxNtracks"  ,"N^{vtx}_{trks};N^{vtx}_{trks};vertices"        ,100,-0.5,99.5);
     h_recoVtxChi2ndf_          = RecoVertexFeatures.make<TH1F>("h_recoVtxChi2ndf"  ,"#chi^{2}/ndf vtx;#chi^{2}/ndf vtx;vertices"    ,10,-0.5,9.5);
     h_recoVtxChi2Prob_         = RecoVertexFeatures.make<TH1F>("h_recoVtxChi2Prob" ,"Prob(#chi^{2},ndf);Prob(#chi^{2},ndf);vertices",40,0.,1.);
     h_recoVtxSumPt_            = RecoVertexFeatures.make<TH1F>("h_recoVtxSumPt"    ,"Sum(p^{trks}_{T});Sum(p^{trks}_{T});vertices"  ,100,0.,200.);
     
   }
 
 
   TFileDirectory DA = fs->mkdir("DA");
   //DA.cd();
   hDA=bookVertexHistograms(DA);
   //for(std::map<std::string,TH1*>::const_iterator hist=hDA.begin(); hist!=hDA.end(); hist++){
   //hist->second->SetDirectory(DA);
   // DA.make<TH1F>(hist->second);
   // }
 
   // initialize the residuals histograms 
 
   const float dxymax_phi = theDetails_.getHigh(PVValHelper::dxy,PVValHelper::phi);
   const float dzmax_phi  = theDetails_.getHigh(PVValHelper::dz ,PVValHelper::eta); 
   const float dxymax_eta = theDetails_.getHigh(PVValHelper::dxy,PVValHelper::phi);
   const float dzmax_eta  = theDetails_.getHigh(PVValHelper::dz, PVValHelper::eta);
   //const float d3Dmax_phi = theDetails_.getHigh(PVValHelper::d3D,PVValHelper::phi);
   const float d3Dmax_eta = theDetails_.getHigh(PVValHelper::d3D,PVValHelper::eta);
 
   //  
   // Unbiased track-to-vertex residuals
   // The vertex is refit without the probe track
   //
 
   //    _   _            _      _         ___        _    _           _    
   //   /_\ | |__ ___ ___| |_  _| |_ ___  | _ \___ __(_)__| |_  _ __ _| |___
   //  / _ \| '_ (_-</ _ \ | || |  _/ -_) |   / -_|_-< / _` | || / _` | (_-<
   // /_/ \_\_.__/__/\___/_|\_,_|\__\___| |_|_\___/__/_\__,_|\_,_\__,_|_/__/
   //
   
   TFileDirectory AbsTransPhiRes  = fs->mkdir("Abs_Transv_Phi_Residuals");
   a_dxyPhiResiduals  = bookResidualsHistogram(AbsTransPhiRes,nBins_,PVValHelper::dxy,PVValHelper::phi);
   a_dxPhiResiduals   = bookResidualsHistogram(AbsTransPhiRes,nBins_,PVValHelper::dx,PVValHelper::phi);
   a_dyPhiResiduals   = bookResidualsHistogram(AbsTransPhiRes,nBins_,PVValHelper::dy,PVValHelper::phi);  
   a_IP2DPhiResiduals = bookResidualsHistogram(AbsTransPhiRes,nBins_,PVValHelper::IP2D,PVValHelper::phi);
 
   TFileDirectory AbsTransEtaRes  = fs->mkdir("Abs_Transv_Eta_Residuals");
   a_dxyEtaResiduals  = bookResidualsHistogram(AbsTransEtaRes,nBins_,PVValHelper::dxy,PVValHelper::eta);
   a_dxEtaResiduals   = bookResidualsHistogram(AbsTransEtaRes,nBins_,PVValHelper::dx,PVValHelper::eta);
   a_dyEtaResiduals   = bookResidualsHistogram(AbsTransEtaRes,nBins_,PVValHelper::dy,PVValHelper::eta);
   a_IP2DEtaResiduals = bookResidualsHistogram(AbsTransEtaRes,nBins_,PVValHelper::IP2D,PVValHelper::eta);
                               
   TFileDirectory AbsLongPhiRes   = fs->mkdir("Abs_Long_Phi_Residuals");
   a_dzPhiResiduals   = bookResidualsHistogram(AbsLongPhiRes,nBins_,PVValHelper::dz,PVValHelper::phi);
   a_reszPhiResiduals = bookResidualsHistogram(AbsLongPhiRes,nBins_,PVValHelper::resz,PVValHelper::phi);
 
   TFileDirectory AbsLongEtaRes   = fs->mkdir("Abs_Long_Eta_Residuals");
   a_dzEtaResiduals   = bookResidualsHistogram(AbsLongEtaRes,nBins_,PVValHelper::dz,PVValHelper::eta);
   a_reszEtaResiduals = bookResidualsHistogram(AbsLongEtaRes,nBins_,PVValHelper::resz,PVValHelper::eta);
                   
   TFileDirectory Abs3DPhiRes     = fs->mkdir("Abs_3D_Phi_Residuals");
   a_d3DPhiResiduals  = bookResidualsHistogram(Abs3DPhiRes,nBins_,PVValHelper::d3D,PVValHelper::phi);
   a_IP3DPhiResiduals = bookResidualsHistogram(Abs3DPhiRes,nBins_,PVValHelper::IP3D,PVValHelper::phi);
 
   TFileDirectory Abs3DEtaRes     = fs->mkdir("Abs_3D_Eta_Residuals");
   a_d3DEtaResiduals  = bookResidualsHistogram(Abs3DEtaRes ,nBins_,PVValHelper::d3D,PVValHelper::eta);
   a_IP3DEtaResiduals = bookResidualsHistogram(Abs3DEtaRes ,nBins_,PVValHelper::IP3D,PVValHelper::eta);
 
   TFileDirectory NormTransPhiRes = fs->mkdir("Norm_Transv_Phi_Residuals");
   n_dxyPhiResiduals  = bookResidualsHistogram(NormTransPhiRes,nBins_,PVValHelper::norm_dxy,PVValHelper::phi,true);
   n_IP2DPhiResiduals = bookResidualsHistogram(NormTransPhiRes,nBins_,PVValHelper::norm_IP2D,PVValHelper::phi,true);
 
   TFileDirectory NormTransEtaRes = fs->mkdir("Norm_Transv_Eta_Residuals");
   n_dxyEtaResiduals  = bookResidualsHistogram(NormTransEtaRes,nBins_,PVValHelper::norm_dxy,PVValHelper::eta,true);
   n_IP2DEtaResiduals = bookResidualsHistogram(NormTransEtaRes,nBins_,PVValHelper::norm_IP2D,PVValHelper::eta,true);
                               
   TFileDirectory NormLongPhiRes  = fs->mkdir("Norm_Long_Phi_Residuals");
   n_dzPhiResiduals   = bookResidualsHistogram(NormLongPhiRes,nBins_,PVValHelper::norm_dz,PVValHelper::phi,true);
   n_reszPhiResiduals = bookResidualsHistogram(NormLongPhiRes,nBins_,PVValHelper::norm_resz,PVValHelper::phi,true);
 
   TFileDirectory NormLongEtaRes  = fs->mkdir("Norm_Long_Eta_Residuals");
   n_dzEtaResiduals   = bookResidualsHistogram(NormLongEtaRes,nBins_,PVValHelper::norm_dz,PVValHelper::eta,true);
   n_reszEtaResiduals = bookResidualsHistogram(NormLongEtaRes,nBins_,PVValHelper::norm_resz,PVValHelper::eta,true);
 
   TFileDirectory Norm3DPhiRes    = fs->mkdir("Norm_3D_Phi_Residuals");
   n_d3DPhiResiduals  = bookResidualsHistogram(Norm3DPhiRes,nBins_,PVValHelper::norm_d3D,PVValHelper::phi,true);
   n_IP3DPhiResiduals = bookResidualsHistogram(Norm3DPhiRes,nBins_,PVValHelper::norm_IP3D,PVValHelper::phi,true);
   
   TFileDirectory Norm3DEtaRes    = fs->mkdir("Norm_3D_Eta_Residuals");
   n_d3DEtaResiduals  = bookResidualsHistogram(Norm3DEtaRes,nBins_,PVValHelper::norm_d3D,PVValHelper::eta,true);
   n_IP3DEtaResiduals = bookResidualsHistogram(Norm3DEtaRes,nBins_,PVValHelper::norm_IP3D,PVValHelper::eta,true);
 
   TFileDirectory AbsDoubleDiffRes   = fs->mkdir("Abs_DoubleDiffResiduals");
   TFileDirectory NormDoubleDiffRes  = fs->mkdir("Norm_DoubleDiffResiduals");
   
   // book residuals vs pT histograms
   
   TFileDirectory AbsTranspTRes  = fs->mkdir("Abs_Transv_pT_Residuals"); 
   h_dxy_pT_      = bookResidualsHistogram(AbsTranspTRes,nPtBins_,PVValHelper::dxy,PVValHelper::pT);        
 
   TFileDirectory AbsLongpTRes   = fs->mkdir("Abs_Long_pT_Residuals"); 
   h_dz_pT_       = bookResidualsHistogram(AbsLongpTRes,nPtBins_,PVValHelper::dz,PVValHelper::pT);              
 
   TFileDirectory NormTranspTRes = fs->mkdir("Norm_Transv_pT_Residuals"); 
   h_norm_dxy_pT_ = bookResidualsHistogram(NormTranspTRes,nPtBins_,PVValHelper::norm_dxy,PVValHelper::pT,true);         
 
   TFileDirectory NormLongpTRes  = fs->mkdir("Norm_Long_pT_Residuals"); 
   h_norm_dz_pT_  = bookResidualsHistogram(NormLongpTRes,nPtBins_,PVValHelper::norm_dz,PVValHelper::pT,true);           
 
   // book residuals vs pT histograms in central region (|eta|<1.0)
                
   TFileDirectory AbsTranspTCentralRes  = fs->mkdir("Abs_Transv_pTCentral_Residuals"); 
   h_dxy_Central_pT_ = bookResidualsHistogram(AbsTranspTCentralRes,nPtBins_,PVValHelper::dxy,PVValHelper::pTCentral);       
 
   TFileDirectory AbsLongpTCentralRes   = fs->mkdir("Abs_Long_pTCentral_Residuals"); 
   h_dz_Central_pT_  = bookResidualsHistogram(AbsLongpTCentralRes,nPtBins_,PVValHelper::dz,PVValHelper::pTCentral);         
 
   TFileDirectory NormTranspTCentralRes = fs->mkdir("Norm_Transv_pTCentral_Residuals"); 
   h_norm_dxy_Central_pT_ = bookResidualsHistogram(NormTranspTCentralRes,nPtBins_,PVValHelper::norm_dxy,PVValHelper::pTCentral,true);  
 
   TFileDirectory NormLongpTCentralRes  = fs->mkdir("Norm_Long_pTCentral_Residuals"); 
   h_norm_dz_Central_pT_  = bookResidualsHistogram(NormLongpTCentralRes,nPtBins_,PVValHelper::norm_dz,PVValHelper::pTCentral,true);   
 
   // book residuals vs module number
   
   TFileDirectory AbsTransModZRes  = fs->mkdir("Abs_Transv_modZ_Residuals"); 
   h_dxy_modZ_      = bookResidualsHistogram(AbsTransModZRes,8,PVValHelper::dxy,PVValHelper::modZ);         
 
   TFileDirectory AbsLongModZRes   = fs->mkdir("Abs_Long_modZ_Residuals"); 
   h_dz_modZ_       = bookResidualsHistogram(AbsLongModZRes,8,PVValHelper::dz,PVValHelper::modZ);               
 
 
   //  _  _                    _ _           _   ___        _    _           _    
   // | \| |___ _ _ _ __  __ _| (_)______ __| | | _ \___ __(_)__| |_  _ __ _| |___
   // | .` / _ \ '_| '  \/ _` | | |_ / -_) _` | |   / -_|_-< / _` | || / _` | (_-<
   // |_|\_\___/_| |_|_|_\__,_|_|_/__\___\__,_| |_|_\___/__/_\__,_|\_,_\__,_|_/__/
   //
   
   TFileDirectory NormTransModZRes = fs->mkdir("Norm_Transv_modZ_Residuals"); 
   h_norm_dxy_modZ_ = bookResidualsHistogram(NormTransModZRes,8,PVValHelper::norm_dxy,PVValHelper::modZ,true);          
 
   TFileDirectory NormLongModZRes  = fs->mkdir("Norm_Long_modZ_Residuals"); 
   h_norm_dz_modZ_  = bookResidualsHistogram(NormLongModZRes,8,PVValHelper::norm_dz,PVValHelper::modZ,true);        
 
   TFileDirectory AbsTransLadderRes  = fs->mkdir("Abs_Transv_ladder_Residuals"); 
   h_dxy_ladder_ = bookResidualsHistogram(AbsTransLadderRes,nLadders_,PVValHelper::dxy,PVValHelper::ladder);
 
   TFileDirectory AbsTransLadderResOverlap  = fs->mkdir("Abs_Transv_ladderOverlap_Residuals"); 
   h_dxy_ladderOverlap_ = bookResidualsHistogram(AbsTransLadderResOverlap,nLadders_,PVValHelper::dxy,PVValHelper::ladder);       
 
   TFileDirectory AbsTransLadderResNoOverlap  = fs->mkdir("Abs_Transv_ladderNoOverlap_Residuals"); 
   h_dxy_ladderNoOverlap_ = bookResidualsHistogram(AbsTransLadderResNoOverlap,nLadders_,PVValHelper::dxy,PVValHelper::ladder);       
 
   TFileDirectory AbsLongLadderRes   = fs->mkdir("Abs_Long_ladder_Residuals"); 
   h_dz_ladder_  = bookResidualsHistogram(AbsLongLadderRes,nLadders_,PVValHelper::dz,PVValHelper::ladder);          
 
   TFileDirectory NormTransLadderRes = fs->mkdir("Norm_Transv_ladder_Residuals"); 
   h_norm_dxy_ladder_ = bookResidualsHistogram(NormTransLadderRes,nLadders_,PVValHelper::norm_dxy,PVValHelper::ladder,true);  
 
   TFileDirectory NormLongLadderRes  = fs->mkdir("Norm_Long_ladder_Residuals"); 
   h_norm_dz_ladder_  = bookResidualsHistogram(NormLongLadderRes,nLadders_,PVValHelper::norm_dz,PVValHelper::ladder,true);   
 
 
   // book residuals as function of phi and eta
 
   for (int i=0; i<nBins_; ++i ) {
 
     float phiF = theDetails_.trendbins[PVValHelper::phi][i];
     float phiL = theDetails_.trendbins[PVValHelper::phi][i+1];
                                                                                     
     //  ___           _    _     ___  _  __  __   ___        _    _           _    
     // |   \ ___ _  _| |__| |___|   \(_)/ _|/ _| | _ \___ __(_)__| |_  _ __ _| |___
     // | |) / _ \ || | '_ \ / -_) |) | |  _|  _| |   / -_|_-< / _` | || / _` | (_-<
     // |___/\___/\_,_|_.__/_\___|___/|_|_| |_|   |_|_\___/__/_\__,_|\_,_\__,_|_/__/
     
     for ( int j=0; j<nBins_; ++j ) {
 
       float etaF = theDetails_.trendbins[PVValHelper::eta][j];
       float etaL = theDetails_.trendbins[PVValHelper::eta][j+1];
 
       a_dxyResidualsMap[i][j] = AbsDoubleDiffRes.make<TH1F>(Form("histo_dxy_eta_plot%i_phi_plot%i",i,j),
                                 Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy};tracks",etaF,etaL,phiF,phiL),
                                 theDetails_.histobins,-dzmax_eta,dzmax_eta);
       
       a_dzResidualsMap[i][j]  = AbsDoubleDiffRes.make<TH1F>(Form("histo_dz_eta_plot%i_phi_plot%i",i,j),
                                 Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z};tracks",etaF,etaL,phiF,phiL),
                                 theDetails_.histobins,-dzmax_eta,dzmax_eta);
       
       a_d3DResidualsMap[i][j] = AbsDoubleDiffRes.make<TH1F>(Form("histo_d3D_eta_plot%i_phi_plot%i",i,j),
                                 Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{3D};tracks",etaF,etaL,phiF,phiL),
                                 theDetails_.histobins,0.,d3Dmax_eta);
       
       n_dxyResidualsMap[i][j] = NormDoubleDiffRes.make<TH1F>(Form("histo_norm_dxy_eta_plot%i_phi_plot%i",i,j),
                                  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy}/#sigma_{d_{xy}};tracks",etaF,etaL,phiF,phiL),
                                  theDetails_.histobins,-dzmax_eta/100,dzmax_eta/100);
 
       n_dzResidualsMap[i][j]  = NormDoubleDiffRes.make<TH1F>(Form("histo_norm_dz_eta_plot%i_phi_plot%i",i,j),
                                  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z}/#sigma_{d_{z}};tracks",etaF,etaL,phiF,phiL),
                                  theDetails_.histobins,-dzmax_eta/100,dzmax_eta/100);
 
       n_d3DResidualsMap[i][j] = NormDoubleDiffRes.make<TH1F>(Form("histo_norm_d3D_eta_plot%i_phi_plot%i",i,j),
                                  Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{3D}/#sigma_{d_{3D}};tracks",etaF,etaL,phiF,phiL),
                                  theDetails_.histobins,0.,d3Dmax_eta);
 
     }
   }
 
   // declaration of the directories
   
   TFileDirectory BiasVsParameter = fs->mkdir("BiasVsParameter");
 
   a_dxyVsPhi = BiasVsParameter.make<TH2F>("h2_dxy_vs_phi","d_{xy} vs track #phi;track #phi [rad];track d_{xy}(PV) [#mum]",
                       nBins_,-M_PI,M_PI,theDetails_.histobins,-dxymax_phi,dxymax_phi); 
  
   a_dzVsPhi  = BiasVsParameter.make<TH2F>("h2_dz_vs_phi","d_{z} vs track #phi;track #phi [rad];track d_{z}(PV) [#mum]",
                       nBins_,-M_PI,M_PI,theDetails_.histobins,-dzmax_phi,dzmax_phi);   
                
   n_dxyVsPhi = BiasVsParameter.make<TH2F>("h2_n_dxy_vs_phi","d_{xy}/#sigma_{d_{xy}} vs track #phi;track #phi [rad];track d_{xy}(PV)/#sigma_{d_{xy}}",
                       nBins_,-M_PI,M_PI,theDetails_.histobins,-dxymax_phi/100.,dxymax_phi/100.); 
   
   n_dzVsPhi  = BiasVsParameter.make<TH2F>("h2_n_dz_vs_phi","d_{z}/#sigma_{d_{z}} vs track #phi;track #phi [rad];track d_{z}(PV)/#sigma_{d_{z}}",
                       nBins_,-M_PI,M_PI,theDetails_.histobins,-dzmax_phi/100.,dzmax_phi/100.);   
                
   a_dxyVsEta = BiasVsParameter.make<TH2F>("h2_dxy_vs_eta","d_{xy} vs track #eta;track #eta;track d_{xy}(PV) [#mum]",
                       nBins_,-etaOfProbe_,etaOfProbe_,theDetails_.histobins,-dxymax_eta,dzmax_eta);
   
   a_dzVsEta  = BiasVsParameter.make<TH2F>("h2_dz_vs_eta","d_{z} vs track #eta;track #eta;track d_{z}(PV) [#mum]",
                       nBins_,-etaOfProbe_,etaOfProbe_,theDetails_.histobins,-dzmax_eta,dzmax_eta);   
                
   n_dxyVsEta = BiasVsParameter.make<TH2F>("h2_n_dxy_vs_eta","d_{xy}/#sigma_{d_{xy}} vs track #eta;track #eta;track d_{xy}(PV)/#sigma_{d_{xy}}",
                       nBins_,-etaOfProbe_,etaOfProbe_,theDetails_.histobins,-dxymax_eta/100.,dxymax_eta/100.);  
 
   n_dzVsEta  = BiasVsParameter.make<TH2F>("h2_n_dz_vs_eta","d_{z}/#sigma_{d_{z}} vs track #eta;track #eta;track d_{z}(PV)/#sigma_{d_{z}}",
                       nBins_,-etaOfProbe_,etaOfProbe_,theDetails_.histobins,-dzmax_eta/100.,dzmax_eta/100.);   
 
   MeanTrendsDir   = fs->mkdir("MeanTrends");
   WidthTrendsDir  = fs->mkdir("WidthTrends");
   MedianTrendsDir = fs->mkdir("MedianTrends");
   MADTrendsDir    = fs->mkdir("MADTrends");
 
   Mean2DMapsDir   = fs->mkdir("MeanMaps");
   Width2DMapsDir  = fs->mkdir("WidthMaps");
 
   double highedge=nBins_-0.5;
   double lowedge=-0.5;
 
   // means and widths from the fit
 
   a_dxyPhiMeanTrend  = MeanTrendsDir.make<TH1F> ("means_dxy_phi",
                          "#LT d_{xy} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy} #GT [#mum]",
                          nBins_,lowedge,highedge); 
   
   a_dxyPhiWidthTrend = WidthTrendsDir.make<TH1F>("widths_dxy_phi",
                          "#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{xy}} [#mum]",
                          nBins_,lowedge,highedge);
   
   a_dzPhiMeanTrend   = MeanTrendsDir.make<TH1F> ("means_dz_phi",
                          "#LT d_{z} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z} #GT [#mum]",
                          nBins_,lowedge,highedge); 
 
   a_dzPhiWidthTrend  = WidthTrendsDir.make<TH1F>("widths_dz_phi","#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{z}} [#mum]",
                          nBins_,lowedge,highedge);
               
   a_dxyEtaMeanTrend  = MeanTrendsDir.make<TH1F> ("means_dxy_eta",
                          "#LT d_{xy} #GT vs #eta sector;#eta (sector);#LT d_{xy} #GT [#mum]",
                          nBins_,lowedge,highedge);
 
   a_dxyEtaWidthTrend = WidthTrendsDir.make<TH1F>("widths_dxy_eta",
                          "#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{xy}} [#mum]",
                          nBins_,lowedge,highedge);
   
   a_dzEtaMeanTrend   = MeanTrendsDir.make<TH1F> ("means_dz_eta",
                          "#LT d_{z} #GT vs #eta sector;#eta (sector);#LT d_{z} #GT [#mum]"
                          ,nBins_,lowedge,highedge); 
   
   a_dzEtaWidthTrend  = WidthTrendsDir.make<TH1F>("widths_dz_eta",
                          "#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{z}} [#mum]",
                          nBins_,lowedge,highedge);
               
   n_dxyPhiMeanTrend  = MeanTrendsDir.make<TH1F> ("norm_means_dxy_phi",
                          "#LT d_{xy}/#sigma_{d_{xy}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy}/#sigma_{d_{xy}} #GT",
                          nBins_,lowedge,highedge);
   
   n_dxyPhiWidthTrend = WidthTrendsDir.make<TH1F>("norm_widths_dxy_phi",
                          "width(d_{xy}/#sigma_{d_{xy}}) vs #phi sector;#varphi (sector) [degrees]; width(d_{xy}/#sigma_{d_{xy}})",
                          nBins_,lowedge,highedge);
   
   n_dzPhiMeanTrend   = MeanTrendsDir.make<TH1F> ("norm_means_dz_phi",
                          "#LT d_{z}/#sigma_{d_{z}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z}/#sigma_{d_{z}} #GT",
                          nBins_,lowedge,highedge); 
   
   n_dzPhiWidthTrend  = WidthTrendsDir.make<TH1F>("norm_widths_dz_phi",
                          "width(d_{z}/#sigma_{d_{z}}) vs #phi sector;#varphi (sector) [degrees];width(d_{z}/#sigma_{d_{z}})",
                          nBins_,lowedge,highedge);
                                             
   n_dxyEtaMeanTrend  = MeanTrendsDir.make<TH1F> ("norm_means_dxy_eta",
                          "#LT d_{xy}/#sigma_{d_{xy}} #GT vs #eta sector;#eta (sector);#LT d_{xy}/#sigma_{d_{z}} #GT",
                          nBins_,lowedge,highedge);
   
   n_dxyEtaWidthTrend = WidthTrendsDir.make<TH1F>("norm_widths_dxy_eta",
                          "width(d_{xy}/#sigma_{d_{xy}}) vs #eta sector;#eta (sector);width(d_{xy}/#sigma_{d_{z}})",
                          nBins_,lowedge,highedge);
 
   n_dzEtaMeanTrend   = MeanTrendsDir.make<TH1F> ("norm_means_dz_eta",
                          "#LT d_{z}/#sigma_{d_{z}} #GT vs #eta sector;#eta (sector);#LT d_{z}/#sigma_{d_{z}} #GT",
                          nBins_,lowedge,highedge);  
   
   n_dzEtaWidthTrend  = WidthTrendsDir.make<TH1F>("norm_widths_dz_eta",
                          "width(d_{z}/#sigma_{d_{z}}) vs #eta sector;#eta (sector);width(d_{z}/#sigma_{d_{z}})",
                          nBins_,lowedge,highedge);                        
   
   // means and widhts vs pT and pTCentral
   
   a_dxypTMeanTrend  = MeanTrendsDir.make<TH1F> ("means_dxy_pT",
                         "#LT d_{xy} #GT vs pT;p_{T} [GeV];#LT d_{xy} #GT [#mum]",
                         mypT_bins_.size()-1,mypT_bins_.data()); 
   
   a_dxypTWidthTrend = WidthTrendsDir.make<TH1F>("widths_dxy_pT",
                         "#sigma_{d_{xy}} vs pT;p_{T} [GeV];#sigma_{d_{xy}} [#mum]",
                         mypT_bins_.size()-1,mypT_bins_.data());
   
   a_dzpTMeanTrend   = MeanTrendsDir.make<TH1F> ("means_dz_pT",
                         "#LT d_{z} #GT vs pT;p_{T} [GeV];#LT d_{z} #GT [#mum]",
                         mypT_bins_.size()-1,mypT_bins_.data()); 
   
   a_dzpTWidthTrend  = WidthTrendsDir.make<TH1F>("widths_dz_pT","#sigma_{d_{z}} vs pT;p_{T} [GeV];#sigma_{d_{z}} [#mum]",
                         mypT_bins_.size()-1,mypT_bins_.data());
   
   
   n_dxypTMeanTrend  = MeanTrendsDir.make<TH1F> ("norm_means_dxy_pT",
                         "#LT d_{xy}/#sigma_{d_{xy}} #GT vs pT;p_{T} [GeV];#LT d_{xy}/#sigma_{d_{xy}} #GT",
                         mypT_bins_.size()-1,mypT_bins_.data());
   
   n_dxypTWidthTrend = WidthTrendsDir.make<TH1F>("norm_widths_dxy_pT",
                         "width(d_{xy}/#sigma_{d_{xy}}) vs pT;p_{T} [GeV]; width(d_{xy}/#sigma_{d_{xy}})",
                         mypT_bins_.size()-1,mypT_bins_.data());
   
   n_dzpTMeanTrend   = MeanTrendsDir.make<TH1F> ("norm_means_dz_pT",
                         "#LT d_{z}/#sigma_{d_{z}} #GT vs pT;p_{T} [GeV];#LT d_{z}/#sigma_{d_{z}} #GT",
                         mypT_bins_.size()-1,mypT_bins_.data()); 
   
   n_dzpTWidthTrend  = WidthTrendsDir.make<TH1F>("norm_widths_dz_pT",
                         "width(d_{z}/#sigma_{d_{z}}) vs pT;p_{T} [GeV];width(d_{z}/#sigma_{d_{z}})",
                         mypT_bins_.size()-1,mypT_bins_.data());
 
 
   a_dxypTCentralMeanTrend  = MeanTrendsDir.make<TH1F> ("means_dxy_pTCentral",
                                "#LT d_{xy} #GT vs p_{T};p_{T}(|#eta|<1.) [GeV];#LT d_{xy} #GT [#mum]",
                                mypT_bins_.size()-1,mypT_bins_.data());
   
   a_dxypTCentralWidthTrend = WidthTrendsDir.make<TH1F>("widths_dxy_pTCentral",
                                "#sigma_{d_{xy}} vs p_{T};p_{T}(|#eta|<1.) [GeV];#sigma_{d_{xy}} [#mum]",
                                mypT_bins_.size()-1,mypT_bins_.data());
   
   a_dzpTCentralMeanTrend   = MeanTrendsDir.make<TH1F> ("means_dz_pTCentral",
                                "#LT d_{z} #GT vs p_{T};p_{T}(|#eta|<1.) [GeV];#LT d_{z} #GT [#mum]"
                                ,mypT_bins_.size()-1,mypT_bins_.data()); 
   
   a_dzpTCentralWidthTrend  = WidthTrendsDir.make<TH1F>("widths_dz_pTCentral",
                                "#sigma_{d_{z}} vs p_{T};p_{T}(|#eta|<1.) [GeV];#sigma_{d_{z}} [#mum]",
                                mypT_bins_.size()-1,mypT_bins_.data());
   
   n_dxypTCentralMeanTrend  = MeanTrendsDir.make<TH1F> ("norm_means_dxy_pTCentral",
                                "#LT d_{xy}/#sigma_{d_{xy}} #GT vs p_{T};p_{T}(|#eta|<1.) [GeV];#LT d_{xy}/#sigma_{d_{z}} #GT",
                                mypT_bins_.size()-1,mypT_bins_.data());
   
   n_dxypTCentralWidthTrend = WidthTrendsDir.make<TH1F>("norm_widths_dxy_pTCentral",
                                "width(d_{xy}/#sigma_{d_{xy}}) vs p_{T};p_{T}(|#eta|<1.) [GeV];width(d_{xy}/#sigma_{d_{z}})",
                                mypT_bins_.size()-1,mypT_bins_.data());
   
   n_dzpTCentralMeanTrend   = MeanTrendsDir.make<TH1F> ("norm_means_dz_pTCentral",
                                "#LT d_{z}/#sigma_{d_{z}} #GT vs p_{T};p_{T}(|#eta|<1.) [GeV];#LT d_{z}/#sigma_{d_{z}} #GT",
                                mypT_bins_.size()-1,mypT_bins_.data());  
   
   n_dzpTCentralWidthTrend  = WidthTrendsDir.make<TH1F>("norm_widths_dz_pTCentral",
                                "width(d_{z}/#sigma_{d_{z}}) vs p_{T};p_{T}(|#eta|<1.) [GeV];width(d_{z}/#sigma_{d_{z}})",
                                mypT_bins_.size()-1,mypT_bins_.data()); 
 
   // 2D maps
 
   a_dxyMeanMap       =  Mean2DMapsDir.make<TH2F>  ("means_dxy_map",
                            "#LT d_{xy} #GT map;#eta (sector);#varphi (sector) [degrees]",
                            nBins_,lowedge,highedge,nBins_,lowedge,highedge);
   
   a_dzMeanMap        =  Mean2DMapsDir.make<TH2F>  ("means_dz_map",
                            "#LT d_{z} #GT map;#eta (sector);#varphi (sector) [degrees]",
                            nBins_,lowedge,highedge,nBins_,lowedge,highedge);
              
   n_dxyMeanMap       =  Mean2DMapsDir.make<TH2F>  ("norm_means_dxy_map",
                            "#LT d_{xy}/#sigma_{d_{xy}} #GT map;#eta (sector);#varphi (sector) [degrees]",
                            nBins_,lowedge,highedge,nBins_,lowedge,highedge);
 
   n_dzMeanMap        =  Mean2DMapsDir.make<TH2F>  ("norm_means_dz_map",
                            "#LT d_{z}/#sigma_{d_{z}} #GT map;#eta (sector);#varphi (sector) [degrees]",
                            nBins_,lowedge,highedge,nBins_,lowedge,highedge);
              
   a_dxyWidthMap      =  Width2DMapsDir.make<TH2F> ("widths_dxy_map",
                            "#sigma_{d_{xy}} map;#eta (sector);#varphi (sector) [degrees]",
                            nBins_,lowedge,highedge,nBins_,lowedge,highedge);
 
   a_dzWidthMap       =  Width2DMapsDir.make<TH2F> ("widths_dz_map",
                            "#sigma_{d_{z}} map;#eta (sector);#varphi (sector) [degrees]",
                            nBins_,lowedge,highedge,nBins_,lowedge,highedge);
              
   n_dxyWidthMap      =  Width2DMapsDir.make<TH2F> ("norm_widths_dxy_map",
                            "width(d_{xy}/#sigma_{d_{xy}}) map;#eta (sector);#varphi (sector) [degrees]",
                            nBins_,lowedge,highedge,nBins_,lowedge,highedge);
   
   n_dzWidthMap       =  Width2DMapsDir.make<TH2F> ("norm_widths_dz_map",
                            "width(d_{z}/#sigma_{d_{z}}) map;#eta (sector);#varphi (sector) [degrees]",
                            nBins_,lowedge,highedge,nBins_,lowedge,highedge);
 
   // medians and MADs
 
   a_dxyPhiMedianTrend = MedianTrendsDir.make<TH1F>("medians_dxy_phi",
                            "Median of d_{xy} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}) [#mum]",
                            nBins_,lowedge,highedge);        
          
   a_dxyPhiMADTrend    = MADTrendsDir.make<TH1F>   ("MADs_dxy_phi",
                            "Median absolute deviation of d_{xy} vs #phi sector;#varphi (sector) [degrees];MAD(d_{xy}) [#mum]",
                            nBins_,lowedge,highedge);
                     
   a_dzPhiMedianTrend  = MedianTrendsDir.make<TH1F>("medians_dz_phi",
                            "Median of d_{z} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}) [#mum]",
                            nBins_,lowedge,highedge); 
                     
   a_dzPhiMADTrend     = MADTrendsDir.make<TH1F>   ("MADs_dz_phi",
                            "Median absolute deviation of d_{z} vs #phi sector;#varphi (sector) [degrees];MAD(d_{z}) [#mum]",
                            nBins_,lowedge,highedge);                    
   
   a_dxyEtaMedianTrend = MedianTrendsDir.make<TH1F>("medians_dxy_eta",
                            "Median of d_{xy} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}) [#mum]",
                            nBins_,lowedge,highedge);        
                   
   a_dxyEtaMADTrend    = MADTrendsDir.make<TH1F>   ("MADs_dxy_eta",
                            "Median absolute deviation of d_{xy} vs #eta sector;#eta (sector);MAD(d_{xy}) [#mum]",
                            nBins_,lowedge,highedge);    
                     
   a_dzEtaMedianTrend  = MedianTrendsDir.make<TH1F>("medians_dz_eta",
                            "Median of d_{z} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}) [#mum]",
                            nBins_,lowedge,highedge);    
                       
   a_dzEtaMADTrend     = MADTrendsDir.make<TH1F>   ("MADs_dz_eta",
                            "Median absolute deviation of d_{z} vs #eta sector;#eta (sector);MAD(d_{z}) [#mum]",
                            nBins_,lowedge,highedge);                          
   
   n_dxyPhiMedianTrend = MedianTrendsDir.make<TH1F>("norm_medians_dxy_phi",
                            "Median of d_{xy}/#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}/#sigma_{d_{xy}})",
                            nBins_,lowedge,highedge); 
 
   n_dxyPhiMADTrend    = MADTrendsDir.make<TH1F>   ("norm_MADs_dxy_phi",
                            "Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees]; MAD(d_{xy}/#sigma_{d_{xy}})",
                            nBins_,lowedge,highedge);   
 
   n_dzPhiMedianTrend  = MedianTrendsDir.make<TH1F>("norm_medians_dz_phi",
                            "Median of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}/#sigma_{d_{z}})",
                            nBins_,lowedge,highedge);  
     
   n_dzPhiMADTrend     = MADTrendsDir.make<TH1F>   ("norm_MADs_dz_phi",
                            "Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];MAD(d_{z}/#sigma_{d_{z}})",
                            nBins_,lowedge,highedge);        
   
   n_dxyEtaMedianTrend = MedianTrendsDir.make<TH1F>("norm_medians_dxy_eta",
                            "Median of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}/#sigma_{d_{z}})",
                            nBins_,lowedge,highedge);
             
   n_dxyEtaMADTrend    = MADTrendsDir.make<TH1F>   ("norm_MADs_dxy_eta",
                            "Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);MAD(d_{xy}/#sigma_{d_{z}})",
                            nBins_,lowedge,highedge);
             
   n_dzEtaMedianTrend  = MedianTrendsDir.make<TH1F>("norm_medians_dz_eta",
                            "Median of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}/#sigma_{d_{z}})",
                            nBins_,lowedge,highedge);    
         
   n_dzEtaMADTrend     = MADTrendsDir.make<TH1F>   ("norm_MADs_dz_eta",
                            "Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);MAD(d_{z}/#sigma_{d_{z}})",
                            nBins_,lowedge,highedge);                      
 
 
   //  
   // plots of biased residuals
   // The vertex includes the probe track
   //
 
   if (useTracksFromRecoVtx_){
 
     TFileDirectory AbsTransPhiBiasRes  = fs->mkdir("Abs_Transv_Phi_BiasResiduals");
     a_dxyPhiBiasResiduals  = bookResidualsHistogram(AbsTransPhiBiasRes,nBins_,PVValHelper::dxy,PVValHelper::phi);
 
     TFileDirectory AbsTransEtaBiasRes  = fs->mkdir("Abs_Transv_Eta_BiasResiduals");
     a_dxyEtaBiasResiduals  = bookResidualsHistogram(AbsTransEtaBiasRes,nBins_,PVValHelper::dxy,PVValHelper::eta);
 
     TFileDirectory AbsLongPhiBiasRes   = fs->mkdir("Abs_Long_Phi_BiasResiduals");
     a_dzPhiBiasResiduals  = bookResidualsHistogram(AbsLongPhiBiasRes,nBins_,PVValHelper::dz,PVValHelper::phi);
                             
     TFileDirectory AbsLongEtaBiasRes   = fs->mkdir("Abs_Long_Eta_BiasResiduals");
     a_dzEtaBiasResiduals  = bookResidualsHistogram(AbsLongEtaBiasRes,nBins_,PVValHelper::dz,PVValHelper::eta);
     
     TFileDirectory NormTransPhiBiasRes = fs->mkdir("Norm_Transv_Phi_BiasResiduals");
     n_dxyPhiBiasResiduals  = bookResidualsHistogram(NormTransPhiBiasRes,nBins_,PVValHelper::dxy,PVValHelper::phi);
     
     TFileDirectory NormTransEtaBiasRes = fs->mkdir("Norm_Transv_Eta_BiasResiduals");
     n_dxyEtaBiasResiduals  = bookResidualsHistogram(NormTransEtaBiasRes,nBins_,PVValHelper::dxy,PVValHelper::eta);
     
     TFileDirectory NormLongPhiBiasRes  = fs->mkdir("Norm_Long_Phi_BiasResiduals");
     n_dzPhiBiasResiduals  = bookResidualsHistogram(NormLongPhiBiasRes,nBins_,PVValHelper::dz,PVValHelper::phi);
 
     TFileDirectory NormLongEtaBiasRes  = fs->mkdir("Norm_Long_Eta_BiasResiduals");
     n_dzEtaBiasResiduals  = bookResidualsHistogram(NormLongEtaBiasRes,nBins_,PVValHelper::dz,PVValHelper::eta);
     
     TFileDirectory AbsDoubleDiffBiasRes   = fs->mkdir("Abs_DoubleDiffBiasResiduals");
     TFileDirectory NormDoubleDiffBiasRes  = fs->mkdir("Norm_DoubleDiffBiasResiduals");
     
     for ( int i=0; i<nBins_; ++i ) {
          
       float phiF = theDetails_.trendbins[PVValHelper::phi][i];
       float phiL = theDetails_.trendbins[PVValHelper::phi][i+1];
       
       for ( int j=0; j<nBins_; ++j ) {
 
     float etaF = theDetails_.trendbins[PVValHelper::eta][j];
     float etaL = theDetails_.trendbins[PVValHelper::eta][j+1];
 
     a_dxyBiasResidualsMap[i][j] = AbsDoubleDiffBiasRes.make<TH1F>(Form("histo_dxy_eta_plot%i_phi_plot%i",i,j),
                                       Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy} [#mum];tracks",etaF,etaL,phiF,phiL),
                                       theDetails_.histobins,-dzmax_eta,dzmax_eta);
     
     a_dzBiasResidualsMap[i][j]  = AbsDoubleDiffBiasRes.make<TH1F>(Form("histo_dxy_eta_plot%i_phi_plot%i",i,j),
                                       Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z} [#mum];tracks",etaF,etaL,phiF,phiL),
                                       theDetails_.histobins,-dzmax_eta,dzmax_eta);
     
     n_dxyBiasResidualsMap[i][j] = NormDoubleDiffBiasRes.make<TH1F>(Form("histo_norm_dxy_eta_plot%i_phi_plot%i",i,j),
                                        Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy}/#sigma_{d_{xy}};tracks",etaF,etaL,phiF,phiL),
                                        theDetails_.histobins,-dzmax_eta/100,dzmax_eta/100);
 
     n_dzBiasResidualsMap[i][j]  = NormDoubleDiffBiasRes.make<TH1F>(Form("histo_norm_dxy_eta_plot%i_phi_plot%i",i,j),
                                        Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z}/#sigma_{d_{z}};tracks",etaF,etaL,phiF,phiL),
                                        theDetails_.histobins,-dzmax_eta/100,dzmax_eta/100);
       }
     }
     
     // declaration of the directories
 
     TFileDirectory MeanBiasTrendsDir   = fs->mkdir("MeanBiasTrends");
     TFileDirectory WidthBiasTrendsDir  = fs->mkdir("WidthBiasTrends");
     TFileDirectory MedianBiasTrendsDir = fs->mkdir("MedianBiasTrends");
     TFileDirectory MADBiasTrendsDir    = fs->mkdir("MADBiasTrends");
     
     TFileDirectory Mean2DBiasMapsDir   = fs->mkdir("MeanBiasMaps");
     TFileDirectory Width2DBiasMapsDir  = fs->mkdir("WidthBiasMaps");
     
     // means and widths from the fit
     
     a_dxyPhiMeanBiasTrend  = MeanBiasTrendsDir.make<TH1F> ("means_dxy_phi",
                                "#LT d_{xy} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy} #GT [#mum]",
                                nBins_,lowedge,highedge); 
 
     a_dxyPhiWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("widths_dxy_phi",
                                "#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{xy}} [#mum]",
                                nBins_,lowedge,highedge);
 
     a_dzPhiMeanBiasTrend   = MeanBiasTrendsDir.make<TH1F> ("means_dz_phi",
                                "#LT d_{z} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z} #GT [#mum]",
                                nBins_,lowedge,highedge); 
 
     a_dzPhiWidthBiasTrend  = WidthBiasTrendsDir.make<TH1F>("widths_dz_phi",
                                "#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{z}} [#mum]",
                                nBins_,lowedge,highedge);
     
     a_dxyEtaMeanBiasTrend  = MeanBiasTrendsDir.make<TH1F> ("means_dxy_eta",
                                "#LT d_{xy} #GT vs #eta sector;#eta (sector);#LT d_{xy} #GT [#mum]",
                                nBins_,lowedge,highedge);
 
     a_dxyEtaWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("widths_dxy_eta",
                                "#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{xy}} [#mum]",
                                nBins_,lowedge,highedge);
 
     a_dzEtaMeanBiasTrend   = MeanBiasTrendsDir.make<TH1F> ("means_dz_eta",
                                "#LT d_{z} #GT vs #eta sector;#eta (sector);#LT d_{z} #GT [#mum]",
                                nBins_,lowedge,highedge); 
 
     a_dzEtaWidthBiasTrend  = WidthBiasTrendsDir.make<TH1F>("widths_dz_eta",
                                "#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{z}} [#mum]",
                                nBins_,lowedge,highedge);
     
     n_dxyPhiMeanBiasTrend  = MeanBiasTrendsDir.make<TH1F> ("norm_means_dxy_phi",
                                "#LT d_{xy}/#sigma_{d_{xy}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy}/#sigma_{d_{xy}} #GT",
                                nBins_,lowedge,highedge);
 
     n_dxyPhiWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("norm_widths_dxy_phi",
                                "width(d_{xy}/#sigma_{d_{xy}}) vs #phi sector;#varphi (sector) [degrees]; width(d_{xy}/#sigma_{d_{xy}})",
                                nBins_,lowedge,highedge);
 
     n_dzPhiMeanBiasTrend   = MeanBiasTrendsDir.make<TH1F> ("norm_means_dz_phi",
                                "#LT d_{z}/#sigma_{d_{z}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z}/#sigma_{d_{z}} #GT",
                                nBins_,lowedge,highedge); 
 
     n_dzPhiWidthBiasTrend  = WidthBiasTrendsDir.make<TH1F>("norm_widths_dz_phi",
                                "width(d_{z}/#sigma_{d_{z}}) vs #phi sector;#varphi (sector) [degrees];width(d_{z}/#sigma_{d_{z}})",
                                nBins_,lowedge,highedge);
     
     n_dxyEtaMeanBiasTrend  = MeanBiasTrendsDir.make<TH1F> ("norm_means_dxy_eta",
                                "#LT d_{xy}/#sigma_{d_{xy}} #GT vs #eta sector;#eta (sector);#LT d_{xy}/#sigma_{d_{z}} #GT",
                                nBins_,lowedge,highedge);
 
     n_dxyEtaWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("norm_widths_dxy_eta",
                                "width(d_{xy}/#sigma_{d_{xy}}) vs #eta sector;#eta (sector);width(d_{xy}/#sigma_{d_{z}})",
                                nBins_,lowedge,highedge);
 
     n_dzEtaMeanBiasTrend   = MeanBiasTrendsDir.make<TH1F> ("norm_means_dz_eta",
                                "#LT d_{z}/#sigma_{d_{z}} #GT vs #eta sector;#eta (sector);#LT d_{z}/#sigma_{d_{z}} #GT",
                                nBins_,lowedge,highedge);  
 
     n_dzEtaWidthBiasTrend  = WidthBiasTrendsDir.make<TH1F>("norm_widths_dz_eta",
                                "width(d_{z}/#sigma_{d_{z}}) vs #eta sector;#eta (sector);width(d_{z}/#sigma_{d_{z}})",
                                nBins_,lowedge,highedge);                        
     
     // 2D maps
     
     a_dxyMeanBiasMap       =  Mean2DBiasMapsDir.make<TH2F>  ("means_dxy_map",
                                  "#LT d_{xy} #GT map;#eta (sector);#varphi (sector) [degrees]",
                                  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
 
     a_dzMeanBiasMap        =  Mean2DBiasMapsDir.make<TH2F>  ("means_dz_map",
                                  "#LT d_{z} #GT map;#eta (sector);#varphi (sector) [degrees]",
                                  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
     
     n_dxyMeanBiasMap       =  Mean2DBiasMapsDir.make<TH2F>  ("norm_means_dxy_map",
                                  "#LT d_{xy}/#sigma_{d_{xy}} #GT map;#eta (sector);#varphi (sector) [degrees]",
                                  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
 
     n_dzMeanBiasMap        =  Mean2DBiasMapsDir.make<TH2F>  ("norm_means_dz_map",
                                  "#LT d_{z}/#sigma_{d_{z}} #GT map;#eta (sector);#varphi (sector) [degrees]",
                                  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
     
     a_dxyWidthBiasMap      =  Width2DBiasMapsDir.make<TH2F> ("widths_dxy_map",
                                  "#sigma_{d_{xy}} map;#eta (sector);#varphi (sector) [degrees]",
                                  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
 
     a_dzWidthBiasMap       =  Width2DBiasMapsDir.make<TH2F> ("widths_dz_map",
                                  "#sigma_{d_{z}} map;#eta (sector);#varphi (sector) [degrees]",
                                  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
     
     n_dxyWidthBiasMap      =  Width2DBiasMapsDir.make<TH2F> ("norm_widths_dxy_map",
                                  "width(d_{xy}/#sigma_{d_{xy}}) map;#eta (sector);#varphi (sector) [degrees]",
                                  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
 
     n_dzWidthBiasMap       =  Width2DBiasMapsDir.make<TH2F> ("norm_widths_dz_map",
                                  "width(d_{z}/#sigma_{d_{z}}) map;#eta (sector);#varphi (sector) [degrees]",
                                  nBins_,lowedge,highedge,nBins_,lowedge,highedge);
     
     // medians and MADs
     
     a_dxyPhiMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("medians_dxy_phi",
                                  "Median of d_{xy} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}) [#mum]",
                                  nBins_,lowedge,highedge);  
              
     a_dxyPhiMADBiasTrend    = MADBiasTrendsDir.make<TH1F>   ("MADs_dxy_phi",
                                  "Median absolute deviation of d_{xy} vs #phi sector;#varphi (sector) [degrees];MAD(d_{xy}) [#mum]",
                                  nBins_,lowedge,highedge);  
                 
     a_dzPhiMedianBiasTrend  = MedianBiasTrendsDir.make<TH1F>("medians_dz_phi",
                                  "Median of d_{z} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}) [#mum]",
                                  nBins_,lowedge,highedge);      
             
     a_dzPhiMADBiasTrend     = MADBiasTrendsDir.make<TH1F>   ("MADs_dz_phi",
                                  "Median absolute deviation of d_{z} vs #phi sector;#varphi (sector) [degrees];MAD(d_{z}) [#mum]",
                                  nBins_,lowedge,highedge);                  
     
     a_dxyEtaMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("medians_dxy_eta",
                                  "Median of d_{xy} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}) [#mum]",
                                  nBins_,lowedge,highedge);  
                       
     a_dxyEtaMADBiasTrend    = MADBiasTrendsDir.make<TH1F>   ("MADs_dxy_eta",
                                  "Median absolute deviation of d_{xy} vs #eta sector;#eta (sector);MAD(d_{xy}) [#mum]",
                                  nBins_,lowedge,highedge);      
                 
     a_dzEtaMedianBiasTrend  = MedianBiasTrendsDir.make<TH1F>("medians_dz_eta",
                                  "Median of d_{z} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}) [#mum]",
                                  nBins_,lowedge,highedge);  
                       
     a_dzEtaMADBiasTrend     = MADBiasTrendsDir.make<TH1F>   ("MADs_dz_eta",
                                  "Median absolute deviation of d_{z} vs #eta sector;#eta (sector);MAD(d_{z}) [#mum]",
                                  nBins_,lowedge,highedge);                        
     
     n_dxyPhiMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("norm_medians_dxy_phi",
                                  "Median of d_{xy}/#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}/#sigma_{d_{xy}})",
                                  nBins_,lowedge,highedge); 
 
     n_dxyPhiMADBiasTrend    = MADBiasTrendsDir.make<TH1F>   ("norm_MADs_dxy_phi",
                                  "Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees]; MAD(d_{xy}/#sigma_{d_{xy}})",
                                  nBins_,lowedge,highedge); 
   
     n_dzPhiMedianBiasTrend  = MedianBiasTrendsDir.make<TH1F>("norm_medians_dz_phi",
                                  "Median of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}/#sigma_{d_{z}})",
                                  nBins_,lowedge,highedge); 
      
     n_dzPhiMADBiasTrend     = MADBiasTrendsDir.make<TH1F>   ("norm_MADs_dz_phi",
                                  "Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];MAD(d_{z}/#sigma_{d_{z}})",
                                  nBins_,lowedge,highedge);      
     
     n_dxyEtaMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("norm_medians_dxy_eta",

Public Member Functions

Private Member Functions

Static Private Member Functions

Private Attributes

Static Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation