#include <SiPixelHitEfficiencySource.h>

Inheritance diagram for SiPixelHitEfficiencySource:

Public Member Functions
virtual void	analyze (const edm::Event &, const edm::EventSetup &)

virtual void	bookHistograms (DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override

virtual void	dqmBeginRun (const edm::Run &r, edm::EventSetup const &iSetup)

virtual void	fillClusterProbability (int, int, bool, double)

	SiPixelHitEfficiencySource (const edm::ParameterSet &)

	~SiPixelHitEfficiencySource ()

Public Member Functions inherited from DQMEDAnalyzer
virtual void	beginRun (edm::Run const &, edm::EventSetup const &) final

virtual void	beginStream (edm::StreamID id) final

	DQMEDAnalyzer (void)

virtual void	endLuminosityBlockSummary (edm::LuminosityBlock const &, edm::EventSetup const &, dqmDetails::NoCache *) const final

virtual void	endRunSummary (edm::Run const &, edm::EventSetup const &, dqmDetails::NoCache *) const final

uint32_t	streamId () const

Public Member Functions inherited from edm::stream::EDAnalyzer< edm::RunSummaryCache< dqmDetails::NoCache >, edm::LuminosityBlockSummaryCache< dqmDetails::NoCache > >
	EDAnalyzer ()=default

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

	EDAnalyzerBase ()

ModuleDescription const &	moduleDescription () const

virtual	~EDAnalyzerBase ()

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

	EDConsumerBase ()

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

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

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

std::vector < ProductHolderIndexAndSkipBit > const &	itemsToGetFromEvent () const

void	labelsForToken (EDGetToken iToken, Labels &oLabels) const

void	modulesDependentUpon (std::string const &iProcessName, std::string const &iModuleLabel, bool iPrint, std::vector< char const * > &oModuleLabels) const

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

bool	registeredToConsume (ProductHolderIndex, bool, BranchType) const

bool	registeredToConsumeMany (TypeID const &, BranchType) const

void	updateLookup (BranchType iBranchType, ProductHolderIndexHelper const &)

virtual	~EDConsumerBase ()

Private Attributes
bool	applyEdgeCut_

bool	bladeOn

edm::EDGetTokenT < edmNew::DetSetVector < SiPixelCluster > >	clusterCollectionToken_

bool	debug_

bool	diskOn

bool	firstRun

bool	isUpgrade

bool	ladOn

bool	layOn

edm::EDGetTokenT < MeasurementTrackerEvent >	measurementTrackerEventToken_

MonitorElement *	meClusterProbabilityD1_Minus_

MonitorElement *	meClusterProbabilityD1_Plus_

MonitorElement *	meClusterProbabilityD2_Minus_

MonitorElement *	meClusterProbabilityD2_Plus_

MonitorElement *	meClusterProbabilityL1_Minus_

MonitorElement *	meClusterProbabilityL1_Plus_

MonitorElement *	meClusterProbabilityL2_Minus_

MonitorElement *	meClusterProbabilityL2_Plus_

MonitorElement *	meClusterProbabilityL3_Minus_

MonitorElement *	meClusterProbabilityL3_Plus_

bool	modOn

int	nmissing

double	nSigma_EdgeCut_

int	nvalid

int	nvtx_

bool	phiOn

edm::ParameterSet	pSet_

bool	ringOn

edm::InputTag	src_

std::map< uint32_t, SiPixelHitEfficiencyModule * >	theSiPixelStructure

edm::EDGetTokenT < TrajTrackAssociationCollection >	tracksrc_

edm::EDGetTokenT < reco::VertexCollection >	vertexCollectionToken_

double	vtxchi2_

double	vtxD0_

double	vtxndof_

int	vtxntrk_

std::string	vtxsrc_

double	vtxX_

double	vtxY_

double	vtxZ_

Additional Inherited Members
Public Types inherited from edm::stream::EDAnalyzer< edm::RunSummaryCache< dqmDetails::NoCache >, edm::LuminosityBlockSummaryCache< dqmDetails::NoCache > >
typedef CacheContexts< T...>	CacheTypes

typedef CacheTypes::GlobalCache	GlobalCache

typedef AbilityChecker< T...>	HasAbility

typedef CacheTypes::LuminosityBlockCache	LuminosityBlockCache

typedef LuminosityBlockContextT < LuminosityBlockCache, RunCache, GlobalCache >	LuminosityBlockContext

typedef CacheTypes::LuminosityBlockSummaryCache	LuminosityBlockSummaryCache

typedef CacheTypes::RunCache	RunCache

typedef RunContextT< RunCache, GlobalCache >	RunContext

typedef CacheTypes::RunSummaryCache	RunSummaryCache

Public Types inherited from edm::stream::EDAnalyzerBase
typedef EDAnalyzerAdaptorBase	ModuleType

Public Types inherited from edm::EDConsumerBase
typedef ProductLabels	Labels

Static Public Member Functions inherited from DQMEDAnalyzer
static std::shared_ptr < dqmDetails::NoCache >	globalBeginLuminosityBlockSummary (edm::LuminosityBlock const &, edm::EventSetup const &, LuminosityBlockContext const *)

static std::shared_ptr < dqmDetails::NoCache >	globalBeginRunSummary (edm::Run const &, edm::EventSetup const &, RunContext const *)

static void	globalEndLuminosityBlockSummary (edm::LuminosityBlock const &, edm::EventSetup const &, LuminosityBlockContext const , dqmDetails::NoCache )

static void	globalEndRunSummary (edm::Run const &, edm::EventSetup const &, RunContext const , dqmDetails::NoCache )

Static Public Member Functions inherited from edm::stream::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 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

Definition at line 44 of file SiPixelHitEfficiencySource.h.

Constructor & Destructor Documentation

SiPixelHitEfficiencySource::SiPixelHitEfficiencySource ( const edm::ParameterSet & pSet )

explicit

Definition at line 73 of file SiPixelHitEfficiencySource.cc.

References applyEdgeCut_, bladeOn, clusterCollectionToken_, debug_, diskOn, firstRun, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), ladOn, layOn, measurementTrackerEventToken_, modOn, nSigma_EdgeCut_, phiOn, pSet_, ringOn, AlCaHLTBitMon_QueryRunRegistry::string, tracksrc_, vertexCollectionToken_, and vtxsrc_.

                                                                                   :
   pSet_(pSet),
   modOn( pSet.getUntrackedParameter<bool>("modOn",true) ),
   ladOn( pSet.getUntrackedParameter<bool>("ladOn",false) ), 
   layOn( pSet.getUntrackedParameter<bool>("layOn",false) ), 
   phiOn( pSet.getUntrackedParameter<bool>("phiOn",false) ), 
   ringOn( pSet.getUntrackedParameter<bool>("ringOn",false) ), 
   bladeOn( pSet.getUntrackedParameter<bool>("bladeOn",false) ), 
   diskOn( pSet.getUntrackedParameter<bool>("diskOn",false) ), 
   isUpgrade( pSet.getUntrackedParameter<bool>("isUpgrade",false) )
   //updateEfficiencies( pSet.getUntrackedParameter<bool>("updateEfficiencies",false) )
 { 
    pSet_ = pSet; 
    debug_ = pSet_.getUntrackedParameter<bool>("debug", false); 
    applyEdgeCut_ = pSet_.getUntrackedParameter<bool>("applyEdgeCut");
    nSigma_EdgeCut_ = pSet_.getUntrackedParameter<double>("nSigma_EdgeCut");
    vtxsrc_= pSet_.getUntrackedParameter<std::string>("vtxsrc","offlinePrimaryVertices");
    vertexCollectionToken_ = consumes<reco::VertexCollection>(vtxsrc_);
    tracksrc_ = consumes<TrajTrackAssociationCollection>(pSet_.getParameter<edm::InputTag>("trajectoryInput"));
    clusterCollectionToken_ = consumes<edmNew::DetSetVector<SiPixelCluster> >(std::string("siPixelClusters"));
 
    measurementTrackerEventToken_ = consumes<MeasurementTrackerEvent>(std::string("MeasurementTrackerEvent"));
 
    firstRun = true;
    
    LogInfo("PixelDQM") << "SiPixelHitEfficiencySource constructor" << endl;
    LogInfo ("PixelDQM") << "Mod/Lad/Lay/Phi " << modOn << "/" << ladOn << "/" 
             << layOn << "/" << phiOn << std::endl;
    LogInfo ("PixelDQM") << "Blade/Disk/Ring" << bladeOn << "/" << diskOn << "/" 
             << ringOn << std::endl;
 }

SiPixelHitEfficiencySource::~SiPixelHitEfficiencySource ( )

Definition at line 106 of file SiPixelHitEfficiencySource.cc.

References theSiPixelStructure.

                                                         {
   LogInfo("PixelDQM") << "SiPixelHitEfficiencySource destructor" << endl;
 
   std::map<uint32_t,SiPixelHitEfficiencyModule*>::iterator struct_iter;
   for (struct_iter = theSiPixelStructure.begin() ; struct_iter != theSiPixelStructure.end() ; struct_iter++){
     delete struct_iter->second;
     struct_iter->second = 0;
   }
 }

Member Function Documentation

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

virtual

Implements edm::stream::EDAnalyzerBase.

Definition at line 261 of file SiPixelHitEfficiencySource.cc.

References funct::abs(), edm::AssociationMap< Tag >::begin(), edmNew::DetSetVector< T >::begin(), bladeOn, clusterCollectionToken_, SiPixelRecHit::clusterProbability(), edm::AssociationMap< edm::OneToOne< std::vector< Trajectory >, reco::TrackCollection, unsigned short > >::const_iterator, gather_cfg::cout, debug_, PixelEndcapName::diskName(), diskOn, edm::AssociationMap< Tag >::end(), edmNew::DetSetVector< T >::end(), HLT_25ns10e33_v2_cff::estimator, f, fillClusterProbability(), edm::EventSetup::get(), edm::Event::getByToken(), PixelClusterParameterEstimator::getParameters(), TrajectoryStateOnSurface::globalPosition(), PixelEndcapName::halfCylinder(), i, edmNew::DetSetVector< T >::id(), TrackerGeometry::idToDetUnit(), PixelBarrelName::isHalfModule(), isUpgrade, edm::ESHandleBase::isValid(), TrajectoryStateOnSurface::isValid(), edm::HandleBase::isValid(), edm::Ref< C, T, F >::key(), PXBDetId::ladder(), ladOn, PixelBarrelName::layerName(), layOn, TrajectoryStateOnSurface::localParameters(), TrajectoryStateOnSurface::localPosition(), match(), LayerMeasurements::measurements(), measurementTrackerEventToken_, TrackingRecHit::missing, modOn, PXBDetId::module(), nmissing, nvalid, nvtx_, oppositeToMomentum, AlCaHLTBitMon_ParallelJobs::p, PixelEndcapName::pannelName(), phiOn, PixelSubdetector::PixelBarrel, PixelSubdetector::PixelEndcap, PixelEndcapName::plaquetteName(), reco::HaloData::plus, edm::Handle< T >::product(), edm::ESHandle< class >::product(), TrajectoryMeasurement::recHit(), ringOn, Propagator::setPropagationDirection(), PixelBarrelName::shell(), edm::AssociationMap< Tag >::size(), mathSSE::sqrt(), DetId::subdetId(), StripSubdetector::TEC, theSiPixelStructure, StripSubdetector::TIB, StripSubdetector::TID, StripSubdetector::TOB, DetId::Tracker, patCandidatesForDimuonsSequences_cff::tracker, tracksrc_, TrackingRecHit::valid, vertexCollectionToken_, HLT_25ns10e33_v2_cff::vertices, vtxchi2_, vtxD0_, vtxndof_, vtxntrk_, vtxX_, vtxY_, vtxZ_, x, PV3DBase< T, PVType, FrameType >::x(), y, PV3DBase< T, PVType, FrameType >::y(), z, and PV3DBase< T, PVType, FrameType >::z().

                                                                                             {
 
   edm::ESHandle<TrackerTopology> tTopoHandle;
   iSetup.get<TrackerTopologyRcd>().get(tTopoHandle);
   const TrackerTopology *pTT = tTopoHandle.product();
 
   edm::Handle<reco::VertexCollection> vertexCollectionHandle;
   iEvent.getByToken( vertexCollectionToken_, vertexCollectionHandle );
   if(!vertexCollectionHandle.isValid()) return;
   nvtx_=0;
   vtxntrk_=-9999;
   vtxD0_=-9999.; vtxX_=-9999.; vtxY_=-9999.; vtxZ_=-9999.; vtxndof_=-9999.; vtxchi2_=-9999.;
   const reco::VertexCollection & vertices = *vertexCollectionHandle.product();
   reco::VertexCollection::const_iterator bestVtx=vertices.end();
   for(reco::VertexCollection::const_iterator it=vertices.begin(); it!=vertices.end(); ++it){
     if(!it->isValid()) continue;
     if(vtxntrk_==-9999 ||
        vtxntrk_<int(it->tracksSize()) ||
        (vtxntrk_==int(it->tracksSize()) && fabs(vtxZ_)>fabs(it->z()))){
       vtxntrk_=it->tracksSize();
       vtxD0_=it->position().rho();
       vtxX_=it->x();
       vtxY_=it->y();
       vtxZ_=it->z();
       vtxndof_=it->ndof();
       vtxchi2_=it->chi2();
       bestVtx=it;
     }
     if(fabs(it->z())<=20. && fabs(it->position().rho())<=2. && it->ndof()>4) nvtx_++;
   }
   if(nvtx_<1) return;
     
   //get the map
   edm::Handle<TrajTrackAssociationCollection> match;
   iEvent.getByToken( tracksrc_, match );
   const TrajTrackAssociationCollection ttac = *(match.product());
 
   if(debug_){
     std::cout << "+++ NEW EVENT +++"<< std::endl;
     std::cout << "Map entries \t : " << ttac.size() << std::endl;
   }
 
   std::set<SiPixelCluster> clusterSet;
   //  TrajectoryStateCombiner tsoscomb;
   //define variables for extrapolation
   int extrapolateFrom_ = 2;
   int extrapolateTo_ = 1;
   float maxlxmatch_=0.2;
   float maxlymatch_=0.2;
   bool keepOriginalMissingHit_=true;
   ESHandle<MeasurementTracker> measurementTrackerHandle;
 
   iSetup.get<CkfComponentsRecord>().get(measurementTrackerHandle);
   
   edm::ESHandle<Chi2MeasurementEstimatorBase> est;
   iSetup.get<TrackingComponentsRecord>().get("Chi2",est);
   edm::Handle<MeasurementTrackerEvent> measurementTrackerEventHandle;
   iEvent.getByToken(measurementTrackerEventToken_, measurementTrackerEventHandle);
   edm::ESHandle<Propagator> prop;
   iSetup.get<TrackingComponentsRecord>().get("PropagatorWithMaterial",prop);
   Propagator* thePropagator = prop.product()->clone();
   //determines direction of the propagator => inward
   if (extrapolateFrom_>=extrapolateTo_) {
     thePropagator->setPropagationDirection(oppositeToMomentum);
   }
   TrajectoryStateCombiner trajStateComb;
   bool debug_ = false;
   
   //Loop over track collection
   for(TrajTrackAssociationCollection::const_iterator it =  ttac.begin();it !=  ttac.end(); ++it){
     //define vector to save extrapolated tracks
     std::vector<TrajectoryMeasurement> expTrajMeasurements;
 
     const edm::Ref<std::vector<Trajectory> > traj_iterator = it->key;  
     // Trajectory Map, extract Trajectory for this track
     reco::TrackRef trackref = it->val;
     //tracks++;
     
     bool isBpixtrack = false, isFpixtrack = false;
     int nStripHits=0; int L1hits=0; int L2hits=0; int L3hits=0; int L4hits=0; int D1hits=0; int D2hits=0; int D3hits=0;
     std::vector<TrajectoryMeasurement> tmeasColl =traj_iterator->measurements();
     std::vector<TrajectoryMeasurement>::const_iterator tmeasIt;
     //loop on measurements to find out what kind of hits there are
     for(tmeasIt = tmeasColl.begin();tmeasIt!=tmeasColl.end();tmeasIt++){
       //if(! tmeasIt->updatedState().isValid()) continue; NOT NECESSARY (I HOPE)
       TransientTrackingRecHit::ConstRecHitPointer testhit = tmeasIt->recHit();
       if(testhit->geographicalId().det() != DetId::Tracker) continue; 
       uint testSubDetID = (testhit->geographicalId().subdetId()); 
       const DetId & hit_detId = testhit->geographicalId();
       int hit_layer = 0;
       int hit_ladder = 0;
       int hit_mod = 0;
       int hit_disk = 0;
       
       if(testSubDetID==PixelSubdetector::PixelBarrel){
         isBpixtrack = true;
     hit_layer = PixelBarrelName(hit_detId,pTT,isUpgrade).layerName();
     
     hit_ladder = PXBDetId(hit_detId).ladder();
     hit_mod = PXBDetId(hit_detId).module();
 
     if(hit_layer==1) L1hits++;
     if(hit_layer==2) L2hits++;
     if(hit_layer==3) L3hits++;
     if(hit_layer==4) L4hits++;
       }
       if(testSubDetID==PixelSubdetector::PixelEndcap){
         isFpixtrack = true;
         hit_disk = PixelEndcapName(hit_detId,pTT,isUpgrade).diskName();
         
     if(hit_disk==1) D1hits++;
     if(hit_disk==2) D2hits++;
         if(hit_disk==3) D3hits++;
       }
       if(testSubDetID==StripSubdetector::TIB) nStripHits++;
       if(testSubDetID==StripSubdetector::TOB) nStripHits++;
       if(testSubDetID==StripSubdetector::TID) nStripHits++;
       if(testSubDetID==StripSubdetector::TEC) nStripHits++;
       //check if last valid hit is in Layer 2 or Disk 1
 
       bool lastValidL2 = false;      
       if ((testSubDetID == PixelSubdetector::PixelBarrel && hit_layer == extrapolateFrom_) 
       || (testSubDetID == PixelSubdetector::PixelEndcap && hit_disk == 1)) {
     if (testhit->isValid()) {
       if (tmeasIt == tmeasColl.end()-1) {
         lastValidL2=true;
       } else {
         tmeasIt++;
         TransientTrackingRecHit::ConstRecHitPointer nextRecHit = tmeasIt->recHit();
         uint nextSubDetID = (nextRecHit->geographicalId().subdetId()); 
         int nextlayer = PixelBarrelName(nextRecHit->geographicalId()).layerName();
         if (nextSubDetID == PixelSubdetector::PixelBarrel && nextlayer==extrapolateTo_ ) {
           lastValidL2=true; //&& !nextRecHit->isValid()) lastValidL2=true;
         }
         tmeasIt--;
       }
     }
       }//end check last valid layer
       if (lastValidL2) {
     std::vector< const BarrelDetLayer*> pxbLayers = measurementTrackerHandle->geometricSearchTracker()->pixelBarrelLayers();
     const DetLayer* pxb1 = pxbLayers[extrapolateTo_-1];
     const MeasurementEstimator* estimator = est.product();
     const LayerMeasurements* theLayerMeasurements =    new LayerMeasurements(*measurementTrackerHandle, *measurementTrackerEventHandle);  
     const TrajectoryStateOnSurface tsosPXB2 = tmeasIt->updatedState();
     expTrajMeasurements = theLayerMeasurements->measurements(*pxb1, tsosPXB2, *thePropagator, *estimator);
     delete theLayerMeasurements;
     if ( !expTrajMeasurements.empty()) {
       for(uint p=0; p<expTrajMeasurements.size();p++){
         TrajectoryMeasurement pxb1TM(expTrajMeasurements[p]);
         auto pxb1Hit = pxb1TM.recHit();
         //remove hits with rawID == 0
         if(pxb1Hit->geographicalId().rawId()==0){
           expTrajMeasurements.erase(expTrajMeasurements.begin()+p);
           continue;
         }
       }
     }
     //
       }
       //check if extrapolated hit to layer 1 one matches the original hit
       TrajectoryStateOnSurface chkPredTrajState=trajStateComb(tmeasIt->forwardPredictedState(), tmeasIt->backwardPredictedState());
       if (!chkPredTrajState.isValid()) continue;
       float chkx=chkPredTrajState.globalPosition().x();
       float chky=chkPredTrajState.globalPosition().y();
       float chkz=chkPredTrajState.globalPosition().z();
       LocalPoint chklp=chkPredTrajState.localPosition();
       if (testSubDetID == PixelSubdetector::PixelBarrel && hit_layer == extrapolateTo_) {
     // Here we will drop the extrapolated hits if there is a hit and use that hit
     vector<int > imatches;
     size_t imatch=0;
     float glmatch=9999.;
     for (size_t iexp=0; iexp<expTrajMeasurements.size(); iexp++) {
       const DetId & exphit_detId = expTrajMeasurements[iexp].recHit()->geographicalId();
       int exphit_ladder = PXBDetId(exphit_detId).ladder();
       int exphit_mod = PXBDetId(exphit_detId).module();
       int dladder = abs( exphit_ladder - hit_ladder ); 
       if (dladder > 10) dladder = 20 - dladder;
       int dmodule = abs( exphit_mod - hit_mod);
       if (dladder != 0 || dmodule != 0) {
         continue;
       }
 
       TrajectoryStateOnSurface predTrajState=expTrajMeasurements[iexp].updatedState();
       float x=predTrajState.globalPosition().x();
       float y=predTrajState.globalPosition().y();
       float z=predTrajState.globalPosition().z();
       float dxyz=sqrt((chkx-x)*(chkx-x)+(chky-y)*(chky-y)+(chkz-z)*(chkz-z));
 
       if (dxyz<=glmatch) {
         glmatch=dxyz;
         imatch=iexp;
         imatches.push_back(int(imatch));
       }
 
     } // found the propagated traj best matching the hit in data
 
     float lxmatch = 9999.0;
     float lymatch = 9999.0;
     if(!expTrajMeasurements.empty()){
       if (glmatch<9999.) { // if there is any propagated trajectory for this hit
         const DetId & matchhit_detId = expTrajMeasurements[imatch].recHit()->geographicalId();
       
         int matchhit_ladder = PXBDetId(matchhit_detId).ladder();
         int dladder = abs(matchhit_ladder-hit_ladder);
         if (dladder > 10) dladder = 20 - dladder;
         LocalPoint lp = expTrajMeasurements[imatch].updatedState().localPosition();
         lxmatch=fabs(lp.x() - chklp.x());
         lymatch=fabs(lp.y() - chklp.y());
       }
       if (lxmatch < maxlxmatch_ && lymatch < maxlymatch_) {
       
         if (testhit->getType()!=TrackingRecHit::missing || keepOriginalMissingHit_) {
           expTrajMeasurements.erase(expTrajMeasurements.begin()+imatch);
         }
         
       }
       
     } //expected trajectory measurment not empty
       }
     }//loop on trajectory measurments tmeasColl
     
     //if an extrapolated hit was found but not matched to an exisitng L1 hit then push the hit back into the collection
     //now keep the first one that is left
     if(!expTrajMeasurements.empty()){
       for (size_t f=0; f<expTrajMeasurements.size(); f++) {
     TrajectoryMeasurement AddHit=expTrajMeasurements[f];
     if (AddHit.recHit()->getType()==TrackingRecHit::missing){
       tmeasColl.push_back(AddHit);
       isBpixtrack = true;
 
     }
     
       }
     }
 
     if(isBpixtrack || isFpixtrack){
       if(trackref->pt()<0.6 ||
          nStripHits<8 ||
      fabs(trackref->dxy(bestVtx->position()))>0.05 ||
      fabs(trackref->dz(bestVtx->position()))>0.5) continue;
     
       if(debug_){
         std::cout << "isBpixtrack : " << isBpixtrack << std::endl;
         std::cout << "isFpixtrack : " << isFpixtrack << std::endl;
       }
       //std::cout<<"This tracks has so many hits: "<<tmeasColl.size()<<std::endl;
       for(std::vector<TrajectoryMeasurement>::const_iterator tmeasIt = tmeasColl.begin(); tmeasIt!=tmeasColl.end(); tmeasIt++){   
     //if(! tmeasIt->updatedState().isValid()) continue; 
     TrajectoryStateOnSurface tsos = tmeasIt->updatedState();
 
     TransientTrackingRecHit::ConstRecHitPointer hit = tmeasIt->recHit();
     if(hit->geographicalId().det() != DetId::Tracker )
       continue; 
     else {
       
       //      //residual
           const DetId & hit_detId = hit->geographicalId();
       //uint IntRawDetID = (hit_detId.rawId());
       uint IntSubDetID = (hit_detId.subdetId());
       
       if(IntSubDetID == 0 ){
         if(debug_) std::cout << "NO IntSubDetID\n";
         continue;
       }
       if(IntSubDetID!=PixelSubdetector::PixelBarrel && IntSubDetID!=PixelSubdetector::PixelEndcap)
         continue;
       
       int disk=0; int layer=0; int panel=0; int module=0; bool isHalfModule=false;
 
       const SiPixelRecHit* pixhit = dynamic_cast<const SiPixelRecHit*>(hit->hit());
 
       if(IntSubDetID==PixelSubdetector::PixelBarrel){ // it's a BPIX hit
         layer = PixelBarrelName(hit_detId,pTT,isUpgrade).layerName();
         isHalfModule = PixelBarrelName(hit_detId,pTT,isUpgrade).isHalfModule();
 
         if (hit->isValid()){ //fill the cluster probability in barrel
           bool plus=true;
           if ((PixelBarrelName(hit_detId,pTT,isUpgrade).shell()== PixelBarrelName::Shell::mO) || (PixelBarrelName(hit_detId,pTT,isUpgrade).shell()==PixelBarrelName::Shell::mI)) plus=false;       
           double clusterProbability= pixhit->clusterProbability(0);      
           if (clusterProbability!=0) fillClusterProbability(layer,0,plus,log10(clusterProbability));          
         }
 
       }else if(IntSubDetID==PixelSubdetector::PixelEndcap){ // it's an FPIX hit
         disk = PixelEndcapName(hit_detId,pTT,isUpgrade).diskName();
         panel = PixelEndcapName(hit_detId,pTT,isUpgrade).pannelName();
         module = PixelEndcapName(hit_detId,pTT,isUpgrade).plaquetteName();
 
         if (hit->isValid()){
           bool plus=true;
           if ((PixelEndcapName(hit_detId,pTT,isUpgrade).halfCylinder()== PixelEndcapName::HalfCylinder::mO) || (PixelEndcapName(hit_detId,pTT,isUpgrade).halfCylinder()==PixelEndcapName::HalfCylinder::mI)) plus=false;
           double clusterProbability= pixhit->clusterProbability(0);
           if (clusterProbability!=0) fillClusterProbability(0,disk,plus,log10(clusterProbability));
         }
           }
               
       if(layer==1){
         if(fabs(trackref->dxy(bestVtx->position()))>0.01 ||
            fabs(trackref->dz(bestVtx->position()))>0.1) continue;
         if(!(L2hits>0&&L3hits>0) && !(L2hits>0&&D1hits>0) && !(D1hits>0&&D2hits>0)) continue;
       }else if(layer==2){
         if(fabs(trackref->dxy(bestVtx->position()))>0.02 ||
            fabs(trackref->dz(bestVtx->position()))>0.1) continue;
         if(!(L1hits>0&&L3hits>0) && !(L1hits>0&&D1hits>0)) continue;
       }else if(layer==3){
         if(fabs(trackref->dxy(bestVtx->position()))>0.02 ||
            fabs(trackref->dz(bestVtx->position()))>0.1) continue;
         if(!(L1hits>0&&L2hits>0)) continue;
       }else if(layer==4){
         if(fabs(trackref->dxy(bestVtx->position()))>0.02 ||
            fabs(trackref->dz(bestVtx->position()))>0.1) continue;
       }else if(disk==1){
         if(fabs(trackref->dxy(bestVtx->position()))>0.05 ||
            fabs(trackref->dz(bestVtx->position()))>0.5) continue;
         if(!(L1hits>0&&D2hits>0) && !(L2hits>0&&D2hits>0)) continue;
       }else if(disk==2){
         if(fabs(trackref->dxy(bestVtx->position()))>0.05 ||
            fabs(trackref->dz(bestVtx->position()))>0.5) continue;
         if(!(L1hits>0&&D1hits>0)) continue;
       }else if(disk==3){
         if(fabs(trackref->dxy(bestVtx->position()))>0.05 ||
            fabs(trackref->dz(bestVtx->position()))>0.5) continue;
       }
       
       //check whether hit is valid or missing using track algo flag
           bool isHitValid   =hit->hit()->getType()==TrackingRecHit::valid;
           bool isHitMissing =hit->hit()->getType()==TrackingRecHit::missing;
       
       if(debug_) std::cout << "the hit is persistent\n";
      
       std::map<uint32_t, SiPixelHitEfficiencyModule*>::iterator pxd = theSiPixelStructure.find((*hit).geographicalId().rawId());
 
       // calculate alpha and beta from cluster position
       LocalTrajectoryParameters ltp = tsos.localParameters();
       //LocalVector localDir = ltp.momentum()/ltp.momentum().mag();
           
       //*************** Edge cut ********************
       double lx=tsos.localPosition().x();
       double ly=tsos.localPosition().y();
 
       if(fabs(lx)>0.55 || fabs(ly)>3.0) continue;
 
       bool passedFiducial=true;
 
       // Module fiducials:
       if(IntSubDetID==PixelSubdetector::PixelBarrel && fabs(ly)>=3.1) passedFiducial=false;
       if(IntSubDetID==PixelSubdetector::PixelEndcap &&
          !((panel==1 &&
         ((module==1 && fabs(ly)<0.7) ||
          ((module==2 && fabs(ly)<1.1) &&
           !(disk==-1 && ly>0.8 && lx>0.2) &&
           !(disk==1 && ly<-0.7 && lx>0.2) &&
           !(disk==2 && ly<-0.8)) ||
          ((module==3 && fabs(ly)<1.5) &&
           !(disk==-2 && lx>0.1 && ly>1.0) &&
           !(disk==2 && lx>0.1 && ly<-1.0)) ||
          ((module==4 && fabs(ly)<1.9) &&
           !(disk==-2 && ly>1.5) &&
           !(disk==2 && ly<-1.5)))) ||
            (panel==2 &&
         ((module==1 && fabs(ly)<0.7) ||
          (module==2 && fabs(ly)<1.2 &&
           !(disk>0 && ly>1.1) &&
           !(disk<0 && ly<-1.1)) ||
          (module==3 && fabs(ly)<1.6 &&
           !(disk>0 && ly>1.5) &&
           !(disk<0 && ly<-1.5)))))) passedFiducial=false;
       if(IntSubDetID==PixelSubdetector::PixelEndcap &&
          ((panel==1 && (module==1 || (module>=3 && abs(disk)==1))) ||
           (panel==2 && ((module==1 && abs(disk)==2) ||
                 (module==3 && abs(disk)==1))))) passedFiducial=false;
       // ROC fiducials:
       double ly_mod = fabs(ly);
       if(IntSubDetID==PixelSubdetector::PixelEndcap && (panel+module)%2==1) ly_mod=ly_mod+0.405;
       float d_rocedge = fabs(fmod(ly_mod,0.81)-0.405);
       if(d_rocedge<=0.0625) passedFiducial=false;
       if(!( (IntSubDetID==PixelSubdetector::PixelBarrel &&
          ((!isHalfModule && fabs(lx)<0.6) ||
           (isHalfModule && lx>-0.3 && lx<0.2))) ||
         (IntSubDetID==PixelSubdetector::PixelEndcap &&
          ((panel==1 &&
            ((module==1 && fabs(lx)<0.2) ||
             (module==2 &&
              ((fabs(lx)<0.55 && abs(disk)==1) ||
               (lx>-0.5 && lx<0.2 && disk==-2) ||
               (lx>-0.5 && lx<0.0 && disk==2))) ||
             (module==3 && lx>-0.6 && lx<0.5) ||
             (module==4 && lx>-0.3 && lx<0.15))) ||
           (panel==2 &&
            ((module==1 && fabs(lx)<0.6) ||
             (module==2 &&
              ((fabs(lx)<0.55 && abs(disk)==1) ||
               (lx>-0.6 && lx<0.5 && abs(disk)==2))) ||
             (module==3 && fabs(lx)<0.5))))))) passedFiducial=false;
       if(((IntSubDetID==PixelSubdetector::PixelBarrel && !isHalfModule) || 
           (IntSubDetID==PixelSubdetector::PixelEndcap && !(panel==1 && (module==1 || module==4)))) &&
          fabs(lx)<0.06) passedFiducial=false;
            
            
       //*************** find closest clusters ********************
       float dx_cl[2]; float dy_cl[2]; dx_cl[0]=dx_cl[1]=dy_cl[0]=dy_cl[1]=-9999.;
       ESHandle<PixelClusterParameterEstimator> cpEstimator;
       iSetup.get<TkPixelCPERecord>().get("PixelCPEGeneric", cpEstimator);
       if(cpEstimator.isValid()){
         const PixelClusterParameterEstimator &cpe(*cpEstimator);
         edm::ESHandle<TrackerGeometry> tracker;
         iSetup.get<TrackerDigiGeometryRecord>().get(tracker);
         if(tracker.isValid()){
           const TrackerGeometry *tkgeom=&(*tracker);
           edm::Handle<edmNew::DetSetVector<SiPixelCluster> > clusterCollectionHandle;
           iEvent.getByToken( clusterCollectionToken_, clusterCollectionHandle );
           if(clusterCollectionHandle.isValid()){
         const edmNew::DetSetVector<SiPixelCluster>& clusterCollection=*clusterCollectionHandle;
         edmNew::DetSetVector<SiPixelCluster>::const_iterator itClusterSet=clusterCollection.begin();
         float minD[2]; minD[0]=minD[1]=10000.;
         for( ; itClusterSet!=clusterCollection.end(); itClusterSet++){
           DetId detId(itClusterSet->id());
           if(detId.rawId()!=hit->geographicalId().rawId()) continue;
           //unsigned int sdId=detId.subdetId();
           const PixelGeomDetUnit *pixdet=(const PixelGeomDetUnit*) tkgeom->idToDetUnit(detId);
           edmNew::DetSet<SiPixelCluster>::const_iterator itCluster=itClusterSet->begin();
           for( ; itCluster!=itClusterSet->end(); ++itCluster){
             LocalPoint lp(itCluster->x(), itCluster->y(), 0.);
             PixelClusterParameterEstimator::ReturnType params=cpe.getParameters(*itCluster,*pixdet);
             lp=std::get<0>(params);
             float D = sqrt((lp.x()-lx)*(lp.x()-lx)+(lp.y()-ly)*(lp.y()-ly));
             if(D<minD[0]){
               minD[1]=minD[0];
               dx_cl[1]=dx_cl[0];
               dy_cl[1]=dy_cl[0];
               minD[0]=D;
               dx_cl[0]=lp.x();
               dy_cl[0]=lp.y();
             }else if(D<minD[1]){
               minD[1]=D;
               dx_cl[1]=lp.x();
               dy_cl[1]=lp.y();
             }
           }  // loop on clusterSets 
         } // loop on clusterCollection
         for(size_t i=0; i<2; i++){
           if(minD[i]<9999.){
             dx_cl[i]=fabs(dx_cl[i]-lx);
             dy_cl[i]=fabs(dy_cl[i]-ly);
           }
         }
           } // valid clusterCollectionHandle
         } // valid tracker
       } // valid cpEstimator
       // distance of hit from closest cluster!
       float d_cl[2]; d_cl[0]=d_cl[1]=-9999.;
       if(dx_cl[0]!=-9999. && dy_cl[0]!=-9999.) d_cl[0]=sqrt(dx_cl[0]*dx_cl[0]+dy_cl[0]*dy_cl[0]);
       if(dx_cl[1]!=-9999. && dy_cl[1]!=-9999.) d_cl[1]=sqrt(dx_cl[1]*dx_cl[1]+dy_cl[1]*dy_cl[1]);
       if(isHitMissing && (d_cl[0]<0.05 || d_cl[1]<0.05)){ isHitMissing=0; isHitValid=1; }         
           
       if(debug_){
         std::cout << "Ready to add hit in histogram:\n";
         //std::cout << "detid: "<<hit_detId<<std::endl;
         std::cout << "isHitValid: "<<isHitValid<<std::endl;
         std::cout << "isHitMissing: "<<isHitMissing<<std::endl;
         //std::cout << "passedEdgeCut: "<<passedFiducial<<std::endl;        
       }    
       
       if (nStripHits<11) continue; //Efficiency plots are filled with hits on tracks that have at least 11 Strip hits 
     
       if(pxd!=theSiPixelStructure.end() && isHitValid && passedFiducial)
         ++nvalid;
       if(pxd!=theSiPixelStructure.end() && isHitMissing && passedFiducial)
         ++nmissing;
         
       if (pxd!=theSiPixelStructure.end() && passedFiducial && (isHitValid || isHitMissing))
         (*pxd).second->fill(pTT,ltp, isHitValid, modOn, ladOn, layOn, phiOn, bladeOn, diskOn, ringOn);  
 
       //}//end if (persistent hit exists and is pixel hit)
       
     }//end of else 
     
     
       }//end for (all traj measurements of pixeltrack)
     }//end if (is pixeltrack)
     else
       if(debug_) std::cout << "no pixeltrack:\n";
     
   }//end loop on map entries
 }

void SiPixelHitEfficiencySource::bookHistograms	(	DQMStore::IBooker &	iBooker,
		edm::Run const &	iRun,
		edm::EventSetup const &	iSetup
	)

overridevirtual

Implements DQMEDAnalyzer.

Definition at line 187 of file SiPixelHitEfficiencySource.cc.

                                                                                                                            {
 
   // book residual histograms in theSiPixelFolder - one (x,y) pair of histograms per det
   SiPixelFolderOrganizer theSiPixelFolder(false);
   for (std::map<uint32_t, SiPixelHitEfficiencyModule*>::iterator pxd = theSiPixelStructure.begin(); 
        pxd!=theSiPixelStructure.end(); pxd++) {
 
     if(modOn){
       if (theSiPixelFolder.setModuleFolder(iBooker,(*pxd).first,0,isUpgrade)) (*pxd).second->book(pSet_,iSetup,iBooker,0,isUpgrade);
       else throw cms::Exception("LogicError") << "SiPixelHitEfficiencySource Folder Creation Failed! "; 
     }
     if(ladOn){
       if (theSiPixelFolder.setModuleFolder(iBooker,(*pxd).first,1,isUpgrade)) (*pxd).second->book(pSet_,iSetup,iBooker,1,isUpgrade);
       else throw cms::Exception("LogicError") << "SiPixelHitEfficiencySource ladder Folder Creation Failed! "; 
     }
     if(layOn){
       if (theSiPixelFolder.setModuleFolder(iBooker,(*pxd).first,2,isUpgrade)) (*pxd).second->book(pSet_,iSetup,iBooker,2,isUpgrade);
       else throw cms::Exception("LogicError") << "SiPixelHitEfficiencySource layer Folder Creation Failed! "; 
     }
     if(phiOn){
       if (theSiPixelFolder.setModuleFolder(iBooker,(*pxd).first,3,isUpgrade)) (*pxd).second->book(pSet_,iSetup,iBooker,3,isUpgrade);
       else throw cms::Exception("LogicError") << "SiPixelHitEfficiencySource phi Folder Creation Failed! "; 
     }
     if(bladeOn){
       if (theSiPixelFolder.setModuleFolder(iBooker,(*pxd).first,4,isUpgrade)) (*pxd).second->book(pSet_,iSetup,iBooker,4,isUpgrade);
       else throw cms::Exception("LogicError") << "SiPixelHitEfficiencySource Blade Folder Creation Failed! "; 
     }
     if(diskOn){
       if (theSiPixelFolder.setModuleFolder(iBooker,(*pxd).first,5,isUpgrade)) (*pxd).second->book(pSet_,iSetup,iBooker,5,isUpgrade);
       else throw cms::Exception("LogicError") << "SiPixelHitEfficiencySource Disk Folder Creation Failed! "; 
     }
     if(ringOn){
       if (theSiPixelFolder.setModuleFolder(iBooker,(*pxd).first,6,isUpgrade)) (*pxd).second->book(pSet_,iSetup,iBooker,6,isUpgrade);
       else throw cms::Exception("LogicError") << "SiPixelHitEfficiencySource Ring Folder Creation Failed! "; 
     }
   }
 
   //book cluster probability histos for Barrel and Endcap
   iBooker.setCurrentFolder("Pixel/Barrel");
 
   meClusterProbabilityL1_Plus_  = iBooker.book1D("ClusterProbabilityXY_Layer1_Plus","ClusterProbabilityXY_Layer1_Plus",500,-14,0.1);
   meClusterProbabilityL1_Plus_->setAxisTitle("Log(ClusterProbability)",1);
 
   meClusterProbabilityL1_Minus_ = iBooker.book1D("ClusterProbabilityXY_Layer1_Minus","ClusterProbabilityXY_Layer1_Minus",500,-14,0.1);
   meClusterProbabilityL1_Minus_->setAxisTitle("Log(ClusterProbability)",1);
 
   meClusterProbabilityL2_Plus_  = iBooker.book1D("ClusterProbabilityXY_Layer2_Plus","ClusterProbabilityXY_Layer2_Plus",500,-14,0.1);
   meClusterProbabilityL2_Plus_ ->setAxisTitle("Log(ClusterProbability)",1);
 
   meClusterProbabilityL2_Minus_ = iBooker.book1D("ClusterProbabilityXY_Layer2_Minus","ClusterProbabilityXY_Layer2_Minus",500,-14,0.1);
   meClusterProbabilityL2_Minus_ ->setAxisTitle("Log(ClusterProbability)",1);
 
   meClusterProbabilityL3_Plus_  = iBooker.book1D("ClusterProbabilityXY_Layer3_Plus","ClusterProbabilityXY_Layer3_Plus",500,-14,0.1);
   meClusterProbabilityL3_Plus_->setAxisTitle("Log(ClusterProbability)",1);
  
   meClusterProbabilityL3_Minus_ = iBooker.book1D("ClusterProbabilityXY_Layer3_Minus","ClusterProbabilityXY_Layer3_Minus",500,-14,0.1);
   meClusterProbabilityL3_Minus_->setAxisTitle("Log(ClusterProbability)",1);
 
   iBooker.setCurrentFolder("Pixel/Endcap");
 
   meClusterProbabilityD1_Plus_  = iBooker.book1D("ClusterProbabilityXY_Disk1_Plus","ClusterProbabilityXY_Disk1_Plus",500,-14,0.1);
   meClusterProbabilityD1_Plus_ ->setAxisTitle("Log(ClusterProbability)",1);
 
   meClusterProbabilityD1_Minus_ = iBooker.book1D("ClusterProbabilityXY_Disk1_Minus","ClusterProbabilityXY_Disk1_Minus",500,-14,0.1);
   meClusterProbabilityD1_Minus_->setAxisTitle("Log(ClusterProbability)",1);
 
   meClusterProbabilityD2_Plus_  = iBooker.book1D("ClusterProbabilityXY_Disk2_Plus","ClusterProbabilityXY_Disk2_Plus",500,-14,0.1);
   meClusterProbabilityD2_Plus_ ->setAxisTitle("Log(ClusterProbability)",1);
 
   meClusterProbabilityD2_Minus_ = iBooker.book1D("ClusterProbabilityXY_Disk2_Minus","ClusterProbabilityXY_Disk2_Minus",500,-14,0.1);
   meClusterProbabilityD2_Minus_->setAxisTitle("Log(ClusterProbability)",1);
 
 }

void SiPixelHitEfficiencySource::dqmBeginRun	(	const edm::Run &	r,
		edm::EventSetup const &	iSetup
	)

virtual

Reimplemented from DQMEDAnalyzer.

Definition at line 152 of file SiPixelHitEfficiencySource.cc.

References debug_, firstRun, edm::EventSetup::get(), python.rootplot.argparse::module, nmissing, nvalid, and theSiPixelStructure.

                                                                                          {
   LogInfo("PixelDQM") << "SiPixelHitEfficiencySource beginRun()" << endl;
   
   if(firstRun){
     // retrieve TrackerGeometry for pixel dets
   
     nvalid=0;
     nmissing=0;
   
     firstRun = false;
   }
 
   edm::ESHandle<TrackerGeometry> TG;
   iSetup.get<TrackerDigiGeometryRecord>().get(TG);
   if (debug_) LogVerbatim("PixelDQM") << "TrackerGeometry "<< &(*TG) <<" size is "<< TG->dets().size() << endl;
  
   // build theSiPixelStructure with the pixel barrel and endcap dets from TrackerGeometry
   for (TrackerGeometry::DetContainer::const_iterator pxb = TG->detsPXB().begin();  
        pxb!=TG->detsPXB().end(); pxb++) {
     if (dynamic_cast<PixelGeomDetUnit const *>((*pxb))!=0) {
       SiPixelHitEfficiencyModule* module = new SiPixelHitEfficiencyModule((*pxb)->geographicalId().rawId());
       theSiPixelStructure.insert(pair<uint32_t, SiPixelHitEfficiencyModule*>((*pxb)->geographicalId().rawId(), module));
     }
   }
   for (TrackerGeometry::DetContainer::const_iterator pxf = TG->detsPXF().begin(); 
        pxf!=TG->detsPXF().end(); pxf++) {
     if (dynamic_cast<PixelGeomDetUnit const *>((*pxf))!=0) {
       SiPixelHitEfficiencyModule* module = new SiPixelHitEfficiencyModule((*pxf)->geographicalId().rawId());
       theSiPixelStructure.insert(pair<uint32_t, SiPixelHitEfficiencyModule*>((*pxf)->geographicalId().rawId(), module));
     }
   }
   LogInfo("PixelDQM") << "SiPixelStructure size is " << theSiPixelStructure.size() << endl;
   
 }

void SiPixelHitEfficiencySource::fillClusterProbability	(	int	layer,
		int	disk,
		bool	plus,
		double	probability
	)

virtual

Definition at line 116 of file SiPixelHitEfficiencySource.cc.

References MonitorElement::Fill(), meClusterProbabilityD1_Minus_, meClusterProbabilityD1_Plus_, meClusterProbabilityD2_Minus_, meClusterProbabilityD2_Plus_, meClusterProbabilityL1_Minus_, meClusterProbabilityL1_Plus_, meClusterProbabilityL2_Minus_, meClusterProbabilityL2_Plus_, meClusterProbabilityL3_Minus_, and meClusterProbabilityL3_Plus_.

Referenced by analyze().

                                                                                                          {
  
   //barrel
   if (layer!=0){
     if (layer==1){
       if (plus)  meClusterProbabilityL1_Plus_->Fill(probability);
       else  meClusterProbabilityL1_Minus_->Fill(probability);
     }
 
     else if (layer==2){
       if (plus)  meClusterProbabilityL2_Plus_->Fill(probability);
       else  meClusterProbabilityL2_Minus_->Fill(probability);
     }
 
     else if (layer==3){
       if (plus)  meClusterProbabilityL3_Plus_->Fill(probability);
       else  meClusterProbabilityL3_Minus_->Fill(probability);
     }
   }
   //Endcap
   if (disk!=0){
     if (disk==1){
       if (plus)  meClusterProbabilityD1_Plus_->Fill(probability);
       else  meClusterProbabilityD1_Minus_->Fill(probability);
     }
     if (disk==2){
       if (plus)  meClusterProbabilityD2_Plus_->Fill(probability);
       else  meClusterProbabilityD2_Minus_->Fill(probability);
     }      
   }    
 }