CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
SiPixelHitEfficiencySource.cc
Go to the documentation of this file.
1 // Package: SiPixelMonitorTrack
2 // Class: SiPixelHitEfficiencySource
3 //
4 // class SiPixelHitEfficiencyModule SiPixelHitEfficiencyModule.cc
5 // DQM/SiPixelMonitorTrack/src/SiPixelHitEfficiencyModule.cc
6 //
7 // Description: SiPixel hit efficiency data quality monitoring modules
8 // Implementation: prototype -> improved -> never final - end of the 1st step
9 //
10 // Original Authors: Romain Rougny & Luca Mucibello
11 // Created: Mar Nov 10 13:29:00 CET 2009
12 
13 
14 #include <iostream>
15 #include <map>
16 #include <string>
17 #include <vector>
18 #include <utility>
19 
22 
25 
34 
38 
40 
43 //#include "RecoLocalTracker/SiPixelRecHits/plugins/PixelCPEGenericESProducer.h"
44 
55 
59 
67 
68 using namespace std;
69 using namespace edm;
70 
71 
73  pSet_(pSet),
74  modOn( pSet.getUntrackedParameter<bool>("modOn",true) ),
75  ladOn( pSet.getUntrackedParameter<bool>("ladOn",false) ),
76  layOn( pSet.getUntrackedParameter<bool>("layOn",false) ),
77  phiOn( pSet.getUntrackedParameter<bool>("phiOn",false) ),
78  ringOn( pSet.getUntrackedParameter<bool>("ringOn",false) ),
79  bladeOn( pSet.getUntrackedParameter<bool>("bladeOn",false) ),
80  diskOn( pSet.getUntrackedParameter<bool>("diskOn",false) ),
81  isUpgrade( pSet.getUntrackedParameter<bool>("isUpgrade",false) )
82  //updateEfficiencies( pSet.getUntrackedParameter<bool>("updateEfficiencies",false) )
83  {
84  pSet_ = pSet;
85  debug_ = pSet_.getUntrackedParameter<bool>("debug", false);
86  applyEdgeCut_ = pSet_.getUntrackedParameter<bool>("applyEdgeCut");
87  nSigma_EdgeCut_ = pSet_.getUntrackedParameter<double>("nSigma_EdgeCut");
88  vertexCollectionToken_ = consumes<reco::VertexCollection>(std::string("offlinePrimaryVertices"));
89  tracksrc_ = consumes<TrajTrackAssociationCollection>(pSet_.getParameter<edm::InputTag>("trajectoryInput"));
90  clusterCollectionToken_ = consumes<edmNew::DetSetVector<SiPixelCluster> >(std::string("siPixelClusters"));
91 
92  measurementTrackerEventToken_ = consumes<MeasurementTrackerEvent>(std::string("MeasurementTrackerEvent"));
93 
94  firstRun = true;
95 
96  LogInfo("PixelDQM") << "SiPixelHitEfficiencySource constructor" << endl;
97  LogInfo ("PixelDQM") << "Mod/Lad/Lay/Phi " << modOn << "/" << ladOn << "/"
98  << layOn << "/" << phiOn << std::endl;
99  LogInfo ("PixelDQM") << "Blade/Disk/Ring" << bladeOn << "/" << diskOn << "/"
100  << ringOn << std::endl;
101 }
102 
103 
105  LogInfo("PixelDQM") << "SiPixelHitEfficiencySource destructor" << endl;
106 
107  std::map<uint32_t,SiPixelHitEfficiencyModule*>::iterator struct_iter;
108  for (struct_iter = theSiPixelStructure.begin() ; struct_iter != theSiPixelStructure.end() ; struct_iter++){
109  delete struct_iter->second;
110  struct_iter->second = 0;
111  }
112 }
113 
114 
116  LogInfo("PixelDQM") << "SiPixelHitEfficiencySource beginRun()" << endl;
117 
118  if(firstRun){
119  // retrieve TrackerGeometry for pixel dets
120 
121  nvalid=0;
122  nmissing=0;
123 
124  firstRun = false;
125  }
126 
128  iSetup.get<TrackerDigiGeometryRecord>().get(TG);
129  if (debug_) LogVerbatim("PixelDQM") << "TrackerGeometry "<< &(*TG) <<" size is "<< TG->dets().size() << endl;
130 
131  // build theSiPixelStructure with the pixel barrel and endcap dets from TrackerGeometry
132  for (TrackerGeometry::DetContainer::const_iterator pxb = TG->detsPXB().begin();
133  pxb!=TG->detsPXB().end(); pxb++) {
134  if (dynamic_cast<PixelGeomDetUnit const *>((*pxb))!=0) {
135  SiPixelHitEfficiencyModule* module = new SiPixelHitEfficiencyModule((*pxb)->geographicalId().rawId());
136  theSiPixelStructure.insert(pair<uint32_t, SiPixelHitEfficiencyModule*>((*pxb)->geographicalId().rawId(), module));
137  }
138  }
139  for (TrackerGeometry::DetContainer::const_iterator pxf = TG->detsPXF().begin();
140  pxf!=TG->detsPXF().end(); pxf++) {
141  if (dynamic_cast<PixelGeomDetUnit const *>((*pxf))!=0) {
142  SiPixelHitEfficiencyModule* module = new SiPixelHitEfficiencyModule((*pxf)->geographicalId().rawId());
143  theSiPixelStructure.insert(pair<uint32_t, SiPixelHitEfficiencyModule*>((*pxf)->geographicalId().rawId(), module));
144  }
145  }
146  LogInfo("PixelDQM") << "SiPixelStructure size is " << theSiPixelStructure.size() << endl;
147 
148 }
149 
151 
152  // book residual histograms in theSiPixelFolder - one (x,y) pair of histograms per det
153  SiPixelFolderOrganizer theSiPixelFolder(false);
154  for (std::map<uint32_t, SiPixelHitEfficiencyModule*>::iterator pxd = theSiPixelStructure.begin();
155  pxd!=theSiPixelStructure.end(); pxd++) {
156 
157  if(modOn){
158  if (theSiPixelFolder.setModuleFolder(iBooker,(*pxd).first,0,isUpgrade)) (*pxd).second->book(pSet_,iSetup,iBooker,0,isUpgrade);
159  else throw cms::Exception("LogicError") << "SiPixelHitEfficiencySource Folder Creation Failed! ";
160  }
161  if(ladOn){
162  if (theSiPixelFolder.setModuleFolder(iBooker,(*pxd).first,1,isUpgrade)) (*pxd).second->book(pSet_,iSetup,iBooker,1,isUpgrade);
163  else throw cms::Exception("LogicError") << "SiPixelHitEfficiencySource ladder Folder Creation Failed! ";
164  }
165  if(layOn){
166  if (theSiPixelFolder.setModuleFolder(iBooker,(*pxd).first,2,isUpgrade)) (*pxd).second->book(pSet_,iSetup,iBooker,2,isUpgrade);
167  else throw cms::Exception("LogicError") << "SiPixelHitEfficiencySource layer Folder Creation Failed! ";
168  }
169  if(phiOn){
170  if (theSiPixelFolder.setModuleFolder(iBooker,(*pxd).first,3,isUpgrade)) (*pxd).second->book(pSet_,iSetup,iBooker,3,isUpgrade);
171  else throw cms::Exception("LogicError") << "SiPixelHitEfficiencySource phi Folder Creation Failed! ";
172  }
173  if(bladeOn){
174  if (theSiPixelFolder.setModuleFolder(iBooker,(*pxd).first,4,isUpgrade)) (*pxd).second->book(pSet_,iSetup,iBooker,4,isUpgrade);
175  else throw cms::Exception("LogicError") << "SiPixelHitEfficiencySource Blade Folder Creation Failed! ";
176  }
177  if(diskOn){
178  if (theSiPixelFolder.setModuleFolder(iBooker,(*pxd).first,5,isUpgrade)) (*pxd).second->book(pSet_,iSetup,iBooker,5,isUpgrade);
179  else throw cms::Exception("LogicError") << "SiPixelHitEfficiencySource Disk Folder Creation Failed! ";
180  }
181  if(ringOn){
182  if (theSiPixelFolder.setModuleFolder(iBooker,(*pxd).first,6,isUpgrade)) (*pxd).second->book(pSet_,iSetup,iBooker,6,isUpgrade);
183  else throw cms::Exception("LogicError") << "SiPixelHitEfficiencySource Ring Folder Creation Failed! ";
184  }
185  }
186 
187 
188 }
189 
191 
192  edm::ESHandle<TrackerTopology> tTopoHandle;
193  iSetup.get<TrackerTopologyRcd>().get(tTopoHandle);
194  const TrackerTopology *pTT = tTopoHandle.product();
195 
196  edm::Handle<reco::VertexCollection> vertexCollectionHandle;
197  iEvent.getByToken( vertexCollectionToken_, vertexCollectionHandle );
198  if(!vertexCollectionHandle.isValid()) return;
199  nvtx_=0;
200  vtxntrk_=-9999;
201  vtxD0_=-9999.; vtxX_=-9999.; vtxY_=-9999.; vtxZ_=-9999.; vtxndof_=-9999.; vtxchi2_=-9999.;
202  const reco::VertexCollection & vertices = *vertexCollectionHandle.product();
203  reco::VertexCollection::const_iterator bestVtx=vertices.end();
204  for(reco::VertexCollection::const_iterator it=vertices.begin(); it!=vertices.end(); ++it){
205  if(!it->isValid()) continue;
206  if(vtxntrk_==-9999 ||
207  vtxntrk_<int(it->tracksSize()) ||
208  (vtxntrk_==int(it->tracksSize()) && fabs(vtxZ_)>fabs(it->z()))){
209  vtxntrk_=it->tracksSize();
210  vtxD0_=it->position().rho();
211  vtxX_=it->x();
212  vtxY_=it->y();
213  vtxZ_=it->z();
214  vtxndof_=it->ndof();
215  vtxchi2_=it->chi2();
216  bestVtx=it;
217  }
218  if(fabs(it->z())<=20. && fabs(it->position().rho())<=2. && it->ndof()>4) nvtx_++;
219  }
220  if(nvtx_<1) return;
221 
222  //get the map
224  iEvent.getByToken( tracksrc_, match );
225  const TrajTrackAssociationCollection ttac = *(match.product());
226 
227  if(debug_){
228  std::cout << "+++ NEW EVENT +++"<< std::endl;
229  std::cout << "Map entries \t : " << ttac.size() << std::endl;
230  }
231 
232  std::set<SiPixelCluster> clusterSet;
233  // TrajectoryStateCombiner tsoscomb;
234  //define variables for extrapolation
235  int extrapolateFrom_ = 2;
236  int extrapolateTo_ = 1;
237  float maxlxmatch_=0.2;
238  float maxlymatch_=0.2;
239  bool keepOriginalMissingHit_=true;
240  ESHandle<MeasurementTracker> measurementTrackerHandle;
241 
242  iSetup.get<CkfComponentsRecord>().get(measurementTrackerHandle);
243 
245  iSetup.get<TrackingComponentsRecord>().get("Chi2",est);
246  edm::Handle<MeasurementTrackerEvent> measurementTrackerEventHandle;
247  iEvent.getByToken(measurementTrackerEventToken_, measurementTrackerEventHandle);
249  iSetup.get<TrackingComponentsRecord>().get("PropagatorWithMaterial",prop);
250  Propagator* thePropagator = prop.product()->clone();
251  //determines direction of the propagator => inward
252  if (extrapolateFrom_>=extrapolateTo_) {
254  }
255  TrajectoryStateCombiner trajStateComb;
256  bool debug_ = false;
257 
258  //Loop over track collection
259  for(TrajTrackAssociationCollection::const_iterator it = ttac.begin();it != ttac.end(); ++it){
260  //define vector to save extrapolated tracks
261  std::vector<TrajectoryMeasurement> expTrajMeasurements;
262 
263  const edm::Ref<std::vector<Trajectory> > traj_iterator = it->key;
264  // Trajectory Map, extract Trajectory for this track
265  reco::TrackRef trackref = it->val;
266  //tracks++;
267 
268  bool isBpixtrack = false, isFpixtrack = false;
269  int nStripHits=0; int L1hits=0; int L2hits=0; int L3hits=0; int L4hits=0; int D1hits=0; int D2hits=0; int D3hits=0;
270  std::vector<TrajectoryMeasurement> tmeasColl =traj_iterator->measurements();
271  std::vector<TrajectoryMeasurement>::const_iterator tmeasIt;
272  //loop on measurements to find out what kind of hits there are
273  for(tmeasIt = tmeasColl.begin();tmeasIt!=tmeasColl.end();tmeasIt++){
274  //if(! tmeasIt->updatedState().isValid()) continue; NOT NECESSARY (I HOPE)
275  TransientTrackingRecHit::ConstRecHitPointer testhit = tmeasIt->recHit();
276  if(testhit->geographicalId().det() != DetId::Tracker) continue;
277  uint testSubDetID = (testhit->geographicalId().subdetId());
278  const DetId & hit_detId = testhit->geographicalId();
279  int hit_layer = 0;
280  int hit_ladder = 0;
281  int hit_mod = 0;
282  int hit_disk = 0;
283 
284  if(testSubDetID==PixelSubdetector::PixelBarrel){
285  isBpixtrack = true;
286  hit_layer = PixelBarrelName(hit_detId,pTT,isUpgrade).layerName();
287 
288  hit_ladder = PXBDetId(hit_detId).ladder();
289  hit_mod = PXBDetId(hit_detId).module();
290 
291  if(hit_layer==1) L1hits++;
292  if(hit_layer==2) L2hits++;
293  if(hit_layer==3) L3hits++;
294  if(hit_layer==4) L4hits++;
295  }
296  if(testSubDetID==PixelSubdetector::PixelEndcap){
297  isFpixtrack = true;
298  hit_disk = PixelEndcapName(hit_detId,pTT,isUpgrade).diskName();
299 
300  if(hit_disk==1) D1hits++;
301  if(hit_disk==2) D2hits++;
302  if(hit_disk==3) D3hits++;
303  }
304  if(testSubDetID==StripSubdetector::TIB) nStripHits++;
305  if(testSubDetID==StripSubdetector::TOB) nStripHits++;
306  if(testSubDetID==StripSubdetector::TID) nStripHits++;
307  if(testSubDetID==StripSubdetector::TEC) nStripHits++;
308  //check if last valid hit is in Layer 2 or Disk 1
309 
310  bool lastValidL2 = false;
311  if ((testSubDetID == PixelSubdetector::PixelBarrel && hit_layer == extrapolateFrom_)
312  || (testSubDetID == PixelSubdetector::PixelEndcap && hit_disk == 1)) {
313  if (testhit->isValid()) {
314  if (tmeasIt == tmeasColl.end()-1) {
315  lastValidL2=true;
316  } else {
317  tmeasIt++;
318  TransientTrackingRecHit::ConstRecHitPointer nextRecHit = tmeasIt->recHit();
319  uint nextSubDetID = (nextRecHit->geographicalId().subdetId());
320  int nextlayer = PixelBarrelName(nextRecHit->geographicalId()).layerName();
321  if (nextSubDetID == PixelSubdetector::PixelBarrel && nextlayer==extrapolateTo_ ) {
322  lastValidL2=true; //&& !nextRecHit->isValid()) lastValidL2=true;
323  }
324  tmeasIt--;
325  }
326  }
327  }//end check last valid layer
328  if (lastValidL2) {
329  std::vector< const BarrelDetLayer*> pxbLayers = measurementTrackerHandle->geometricSearchTracker()->pixelBarrelLayers();
330  const DetLayer* pxb1 = pxbLayers[extrapolateTo_-1];
331  const MeasurementEstimator* estimator = est.product();
332  const LayerMeasurements* theLayerMeasurements = new LayerMeasurements(*measurementTrackerHandle, *measurementTrackerEventHandle);
333  const TrajectoryStateOnSurface tsosPXB2 = tmeasIt->updatedState();
334  expTrajMeasurements = theLayerMeasurements->measurements(*pxb1, tsosPXB2, *thePropagator, *estimator);
335  delete theLayerMeasurements;
336  if ( !expTrajMeasurements.empty()) {
337  for(uint p=0; p<expTrajMeasurements.size();p++){
338  TrajectoryMeasurement pxb1TM(expTrajMeasurements[p]);
339  auto pxb1Hit = pxb1TM.recHit();
340  //remove hits with rawID == 0
341  if(pxb1Hit->geographicalId().rawId()==0){
342  expTrajMeasurements.erase(expTrajMeasurements.begin()+p);
343  continue;
344  }
345  }
346  }
347  //
348  }
349  //check if extrapolated hit to layer 1 one matches the original hit
350  TrajectoryStateOnSurface chkPredTrajState=trajStateComb(tmeasIt->forwardPredictedState(), tmeasIt->backwardPredictedState());
351  float chkx=chkPredTrajState.globalPosition().x();
352  float chky=chkPredTrajState.globalPosition().y();
353  float chkz=chkPredTrajState.globalPosition().z();
354  LocalPoint chklp=chkPredTrajState.localPosition();
355  if (testSubDetID == PixelSubdetector::PixelBarrel && hit_layer == extrapolateTo_) {
356  // Here we will drop the extrapolated hits if there is a hit and use that hit
357  vector<int > imatches;
358  size_t imatch=0;
359  float glmatch=9999.;
360  for (size_t iexp=0; iexp<expTrajMeasurements.size(); iexp++) {
361  const DetId & exphit_detId = expTrajMeasurements[iexp].recHit()->geographicalId();
362  int exphit_ladder = PXBDetId(exphit_detId).ladder();
363  int exphit_mod = PXBDetId(exphit_detId).module();
364  int dladder = abs( exphit_ladder - hit_ladder );
365  if (dladder > 10) dladder = 20 - dladder;
366  int dmodule = abs( exphit_mod - hit_mod);
367  if (dladder != 0 || dmodule != 0) {
368  continue;
369  }
370 
371  TrajectoryStateOnSurface predTrajState=expTrajMeasurements[iexp].updatedState();
372  float x=predTrajState.globalPosition().x();
373 float y=predTrajState.globalPosition().y();
374  float z=predTrajState.globalPosition().z();
375  float dxyz=sqrt((chkx-x)*(chkx-x)+(chky-y)*(chky-y)+(chkz-z)*(chkz-z));
376 
377  if (dxyz<=glmatch) {
378  glmatch=dxyz;
379  imatch=iexp;
380  imatches.push_back(int(imatch));
381  }
382 
383  } // found the propagated traj best matching the hit in data
384 
385  float lxmatch = 9999.0;
386  float lymatch = 9999.0;
387  if(!expTrajMeasurements.empty()){
388  if (glmatch<9999.) { // if there is any propagated trajectory for this hit
389  const DetId & matchhit_detId = expTrajMeasurements[imatch].recHit()->geographicalId();
390 
391  int matchhit_ladder = PXBDetId(matchhit_detId).ladder();
392  int dladder = abs(matchhit_ladder-hit_ladder);
393  if (dladder > 10) dladder = 20 - dladder;
394  LocalPoint lp = expTrajMeasurements[imatch].updatedState().localPosition();
395  lxmatch=fabs(lp.x() - chklp.x());
396  lymatch=fabs(lp.y() - chklp.y());
397  }
398  if (lxmatch < maxlxmatch_ && lymatch < maxlymatch_) {
399 
400  if (testhit->getType()!=TrackingRecHit::missing || keepOriginalMissingHit_) {
401  expTrajMeasurements.erase(expTrajMeasurements.begin()+imatch);
402  }
403 
404  }
405 
406  } //expected trajectory measurment not empty
407  }
408  }//loop on trajectory measurments tmeasColl
409 
410  //if an extrapolated hit was found but not matched to an exisitng L1 hit then push the hit back into the collection
411  //now keep the first one that is left
412  if(!expTrajMeasurements.empty()){
413  for (size_t f=0; f<expTrajMeasurements.size(); f++) {
414  TrajectoryMeasurement AddHit=expTrajMeasurements[f];
415  if (AddHit.recHit()->getType()==TrackingRecHit::missing){
416  tmeasColl.push_back(AddHit);
417  isBpixtrack = true;
418 
419  }
420 
421  }
422 
423 
424 
425  }
426  if(isBpixtrack || isFpixtrack){
427  if(trackref->pt()<0.6 ||
428  nStripHits<11 ||
429  fabs(trackref->dxy(bestVtx->position()))>0.05 ||
430  fabs(trackref->dz(bestVtx->position()))>0.5) continue;
431 
432  if(debug_){
433  std::cout << "isBpixtrack : " << isBpixtrack << std::endl;
434  std::cout << "isFpixtrack : " << isFpixtrack << std::endl;
435  }
436  //std::cout<<"This tracks has so many hits: "<<tmeasColl.size()<<std::endl;
437  for(std::vector<TrajectoryMeasurement>::const_iterator tmeasIt = tmeasColl.begin(); tmeasIt!=tmeasColl.end(); tmeasIt++){
438  //if(! tmeasIt->updatedState().isValid()) continue;
439  TrajectoryStateOnSurface tsos = tmeasIt->updatedState();
440 
442  if(hit->geographicalId().det() != DetId::Tracker )
443  continue;
444  else {
445 
446 // //residual
447  const DetId & hit_detId = hit->geographicalId();
448  //uint IntRawDetID = (hit_detId.rawId());
449  uint IntSubDetID = (hit_detId.subdetId());
450 
451  if(IntSubDetID == 0 ){
452  if(debug_) std::cout << "NO IntSubDetID\n";
453  continue;
454  }
455  if(IntSubDetID!=PixelSubdetector::PixelBarrel && IntSubDetID!=PixelSubdetector::PixelEndcap)
456  continue;
457 
458  int disk=0; int layer=0; int panel=0; int module=0; bool isHalfModule=false;
459  if(IntSubDetID==PixelSubdetector::PixelBarrel){ // it's a BPIX hit
460  layer = PixelBarrelName(hit_detId,pTT,isUpgrade).layerName();
461  isHalfModule = PixelBarrelName(hit_detId,pTT,isUpgrade).isHalfModule();
462  }else if(IntSubDetID==PixelSubdetector::PixelEndcap){ // it's an FPIX hit
463  disk = PixelEndcapName(hit_detId,pTT,isUpgrade).diskName();
464  panel = PixelEndcapName(hit_detId,pTT,isUpgrade).pannelName();
465  module = PixelEndcapName(hit_detId,pTT,isUpgrade).plaquetteName();
466  }
467 
468  if(layer==1){
469  if(fabs(trackref->dxy(bestVtx->position()))>0.01 ||
470  fabs(trackref->dz(bestVtx->position()))>0.1) continue;
471  if(!(L2hits>0&&L3hits>0&&L4hits>0) && !(L2hits>0&&D1hits>0&&D2hits) && !(D1hits>0&&D2hits>0&&D3hits>0)) continue;
472  }else if(layer==2){
473  if(fabs(trackref->dxy(bestVtx->position()))>0.02 ||
474  fabs(trackref->dz(bestVtx->position()))>0.1) continue;
475  if(!(L1hits>0&&L3hits>0&&L4hits>0) && !(L1hits>0&&L3hits>0&&D1hits>0) && !(L1hits>0&&D1hits>0&&D2hits>0)) continue;
476  }else if(layer==3){
477  if(fabs(trackref->dxy(bestVtx->position()))>0.02 ||
478  fabs(trackref->dz(bestVtx->position()))>0.1) continue;
479  if(!(L1hits>0&&L2hits>0&&L4hits>0) && !(L1hits>0&&L2hits>0&&D1hits>0)) continue;
480  }else if(layer==4){
481  if(fabs(trackref->dxy(bestVtx->position()))>0.02 ||
482  fabs(trackref->dz(bestVtx->position()))>0.1) continue;
483  if(!(L1hits>0&&L2hits>0&&L3hits>0)) continue;
484  }else if(disk==1){
485  if(fabs(trackref->dxy(bestVtx->position()))>0.05 ||
486  fabs(trackref->dz(bestVtx->position()))>0.5) continue;
487  if(!(L1hits>0&&L2hits>0&&D2hits>0) && !(L1hits>0&&D2hits>0&&D3hits>0) && !(L2hits>0&&D2hits>0&&D3hits>0)) continue;
488  }else if(disk==2){
489  if(fabs(trackref->dxy(bestVtx->position()))>0.05 ||
490  fabs(trackref->dz(bestVtx->position()))>0.5) continue;
491  if(!(L1hits>0&&L2hits>0&&D1hits>0) && !(L1hits>0&&D1hits>0&&D3hits>0) && !(L2hits>0&&D1hits>0&&D3hits>0)) continue;
492  }else if(disk==3){
493  if(fabs(trackref->dxy(bestVtx->position()))>0.05 ||
494  fabs(trackref->dz(bestVtx->position()))>0.5) continue;
495  if(!(L1hits>0&&D1hits>0&&D2hits>0) && !(L2hits>0&&D1hits>0&&D2hits>0)) continue;
496  }
497 
498  //check whether hit is valid or missing using track algo flag
499  bool isHitValid =hit->hit()->getType()==TrackingRecHit::valid;
500  bool isHitMissing =hit->hit()->getType()==TrackingRecHit::missing;
501 
502  if(debug_) std::cout << "the hit is persistent\n";
503 
504  std::map<uint32_t, SiPixelHitEfficiencyModule*>::iterator pxd = theSiPixelStructure.find((*hit).geographicalId().rawId());
505 
506  // calculate alpha and beta from cluster position
508  //LocalVector localDir = ltp.momentum()/ltp.momentum().mag();
509 
510  //*************** Edge cut ********************
511  double lx=tsos.localPosition().x();
512  double ly=tsos.localPosition().y();
513 
514  if(fabs(lx)>0.55 || fabs(ly)>3.0) continue;
515 
516  bool passedFiducial=true;
517  // Module fiducials:
518  if(IntSubDetID==PixelSubdetector::PixelBarrel && fabs(ly)>=3.1) passedFiducial=false;
519  if(IntSubDetID==PixelSubdetector::PixelEndcap &&
520  !((panel==1 &&
521  ((module==1 && fabs(ly)<0.7) ||
522  ((module==2 && fabs(ly)<1.1) &&
523  !(disk==-1 && ly>0.8 && lx>0.2) &&
524  !(disk==1 && ly<-0.7 && lx>0.2) &&
525  !(disk==2 && ly<-0.8)) ||
526  ((module==3 && fabs(ly)<1.5) &&
527  !(disk==-2 && lx>0.1 && ly>1.0) &&
528  !(disk==2 && lx>0.1 && ly<-1.0)) ||
529  ((module==4 && fabs(ly)<1.9) &&
530  !(disk==-2 && ly>1.5) &&
531  !(disk==2 && ly<-1.5)))) ||
532  (panel==2 &&
533  ((module==1 && fabs(ly)<0.7) ||
534  (module==2 && fabs(ly)<1.2 &&
535  !(disk>0 && ly>1.1) &&
536  !(disk<0 && ly<-1.1)) ||
537  (module==3 && fabs(ly)<1.6 &&
538  !(disk>0 && ly>1.5) &&
539  !(disk<0 && ly<-1.5)))))) passedFiducial=false;
540  if(IntSubDetID==PixelSubdetector::PixelEndcap &&
541  ((panel==1 && (module==1 || (module>=3 && abs(disk)==1))) ||
542  (panel==2 && ((module==1 && abs(disk)==2) ||
543  (module==3 && abs(disk)==1))))) passedFiducial=false;
544  // ROC fiducials:
545  double ly_mod = fabs(ly);
546  if(IntSubDetID==PixelSubdetector::PixelEndcap && (panel+module)%2==1) ly_mod=ly_mod+0.405;
547  float d_rocedge = fabs(fmod(ly_mod,0.81)-0.405);
548  if(d_rocedge<=0.0625) passedFiducial=false;
549  if(!( (IntSubDetID==PixelSubdetector::PixelBarrel &&
550  ((!isHalfModule && fabs(lx)<0.6) ||
551  (isHalfModule && lx>-0.3 && lx<0.2))) ||
552  (IntSubDetID==PixelSubdetector::PixelEndcap &&
553  ((panel==1 &&
554  ((module==1 && fabs(lx)<0.2) ||
555  (module==2 &&
556  ((fabs(lx)<0.55 && abs(disk)==1) ||
557  (lx>-0.5 && lx<0.2 && disk==-2) ||
558  (lx>-0.5 && lx<0.0 && disk==2))) ||
559  (module==3 && lx>-0.6 && lx<0.5) ||
560  (module==4 && lx>-0.3 && lx<0.15))) ||
561  (panel==2 &&
562  ((module==1 && fabs(lx)<0.6) ||
563  (module==2 &&
564  ((fabs(lx)<0.55 && abs(disk)==1) ||
565  (lx>-0.6 && lx<0.5 && abs(disk)==2))) ||
566  (module==3 && fabs(lx)<0.5))))))) passedFiducial=false;
567  if(((IntSubDetID==PixelSubdetector::PixelBarrel && !isHalfModule) ||
568  (IntSubDetID==PixelSubdetector::PixelEndcap && !(panel==1 && (module==1 || module==4)))) &&
569  fabs(lx)<0.06) passedFiducial=false;
570 
571 
572  //*************** find closest clusters ********************
573  float dx_cl[2]; float dy_cl[2]; dx_cl[0]=dx_cl[1]=dy_cl[0]=dy_cl[1]=-9999.;
575  iSetup.get<TkPixelCPERecord>().get("PixelCPEGeneric", cpEstimator);
576  if(cpEstimator.isValid()){
577  const PixelClusterParameterEstimator &cpe(*cpEstimator);
579  iSetup.get<TrackerDigiGeometryRecord>().get(tracker);
580  if(tracker.isValid()){
581  const TrackerGeometry *tkgeom=&(*tracker);
582  edm::Handle<edmNew::DetSetVector<SiPixelCluster> > clusterCollectionHandle;
583  iEvent.getByToken( clusterCollectionToken_, clusterCollectionHandle );
584  if(clusterCollectionHandle.isValid()){
585  const edmNew::DetSetVector<SiPixelCluster>& clusterCollection=*clusterCollectionHandle;
586  edmNew::DetSetVector<SiPixelCluster>::const_iterator itClusterSet=clusterCollection.begin();
587  float minD[2]; minD[0]=minD[1]=10000.;
588  for( ; itClusterSet!=clusterCollection.end(); itClusterSet++){
589  DetId detId(itClusterSet->id());
590  if(detId.rawId()!=hit->geographicalId().rawId()) continue;
591  //unsigned int sdId=detId.subdetId();
592  const PixelGeomDetUnit *pixdet=(const PixelGeomDetUnit*) tkgeom->idToDetUnit(detId);
593  edmNew::DetSet<SiPixelCluster>::const_iterator itCluster=itClusterSet->begin();
594  for( ; itCluster!=itClusterSet->end(); ++itCluster){
595  LocalPoint lp(itCluster->x(), itCluster->y(), 0.);
596  PixelClusterParameterEstimator::ReturnType params=cpe.getParameters(*itCluster,*pixdet);
597  lp=std::get<0>(params);
598  float D = sqrt((lp.x()-lx)*(lp.x()-lx)+(lp.y()-ly)*(lp.y()-ly));
599  if(D<minD[0]){
600  minD[1]=minD[0];
601  dx_cl[1]=dx_cl[0];
602  dy_cl[1]=dy_cl[0];
603  minD[0]=D;
604  dx_cl[0]=lp.x();
605  dy_cl[0]=lp.y();
606  }else if(D<minD[1]){
607  minD[1]=D;
608  dx_cl[1]=lp.x();
609  dy_cl[1]=lp.y();
610  }
611  } // loop on clusterSets
612  } // loop on clusterCollection
613  for(size_t i=0; i<2; i++){
614  if(minD[i]<9999.){
615  dx_cl[i]=fabs(dx_cl[i]-lx);
616  dy_cl[i]=fabs(dy_cl[i]-ly);
617  }
618  }
619  } // valid clusterCollectionHandle
620  } // valid tracker
621  } // valid cpEstimator
622  // distance of hit from closest cluster!
623  float d_cl[2]; d_cl[0]=d_cl[1]=-9999.;
624  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]);
625  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]);
626  if(isHitMissing && (d_cl[0]<0.05 || d_cl[1]<0.05)){ isHitMissing=0; isHitValid=1; }
627 
628 
629  if(debug_){
630  std::cout << "Ready to add hit in histogram:\n";
631  //std::cout << "detid: "<<hit_detId<<std::endl;
632  std::cout << "isHitValid: "<<isHitValid<<std::endl;
633  std::cout << "isHitMissing: "<<isHitMissing<<std::endl;
634  //std::cout << "passedEdgeCut: "<<passedFiducial<<std::endl;
635  }
636 
637 
638  if(pxd!=theSiPixelStructure.end() && isHitValid && passedFiducial)
639  ++nvalid;
640  if(pxd!=theSiPixelStructure.end() && isHitMissing && passedFiducial)
641  ++nmissing;
642 
643  if (pxd!=theSiPixelStructure.end() && passedFiducial && (isHitValid || isHitMissing))
644  (*pxd).second->fill(pTT,ltp, isHitValid, modOn, ladOn, layOn, phiOn, bladeOn, diskOn, ringOn);
645 
646  //}//end if (persistent hit exists and is pixel hit)
647 
648  }//end of else
649 
650 
651  }//end for (all traj measurements of pixeltrack)
652  }//end if (is pixeltrack)
653  else
654  if(debug_) std::cout << "no pixeltrack:\n";
655 
656  }//end loop on map entries
657 }
658 
659 
660 DEFINE_FWK_MODULE(SiPixelHitEfficiencySource); // define this as a plug-in
int plaquetteName() const
plaquetteId (in pannel)
edm::EDGetTokenT< MeasurementTrackerEvent > measurementTrackerEventToken_
T getParameter(std::string const &) const
virtual const TrackerGeomDet * idToDetUnit(DetId) const
Return the pointer to the GeomDetUnit corresponding to a given DetId.
T getUntrackedParameter(std::string const &, T const &) const
std::vector< TrajectoryMeasurement > measurements(const DetLayer &layer, const TrajectoryStateOnSurface &startingState, const Propagator &prop, const MeasurementEstimator &est) const
int i
Definition: DBlmapReader.cc:9
virtual void setPropagationDirection(PropagationDirection dir)
Definition: Propagator.h:144
boost::transform_iterator< IterHelp, const_IdIter > const_iterator
const_iterator end(bool update=false) const
ConstRecHitPointer const & recHit() const
SiPixelHitEfficiencySource(const edm::ParameterSet &)
const LocalTrajectoryParameters & localParameters() const
const_iterator end() const
last iterator over the map (read only)
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:457
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
virtual void analyze(const edm::Event &, const edm::EventSetup &)
edm::EDGetTokenT< edmNew::DetSetVector< SiPixelCluster > > clusterCollectionToken_
std::map< uint32_t, SiPixelHitEfficiencyModule * > theSiPixelStructure
T y() const
Definition: PV3DBase.h:63
GlobalPoint globalPosition() const
std::vector< Vertex > VertexCollection
collection of Vertex objects
Definition: VertexFwd.h:9
unsigned int ladder() const
ladder id
Definition: PXBDetId.h:39
data_type const * const_iterator
Definition: DetSetNew.h:30
key_type key() const
Accessor for product key.
Definition: Ref.h:264
id_type id(size_t cell) const
virtual void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override
edm::EDGetTokenT< reco::VertexCollection > vertexCollectionToken_
virtual ReturnType getParameters(const SiPixelCluster &cl, const GeomDetUnit &det) const =0
int iEvent
Definition: GenABIO.cc:230
bool isHalfModule() const
full or half module
std::shared_ptr< TrackingRecHit const > ConstRecHitPointer
T sqrt(T t)
Definition: SSEVec.h:48
bool setModuleFolder(const uint32_t &rawdetid=0, int type=0, bool isUpgrade=false)
Set folder name for a module or plaquette.
edm::EDGetTokenT< TrajTrackAssociationCollection > tracksrc_
T z() const
Definition: PV3DBase.h:64
unsigned int module() const
det id
Definition: PXBDetId.h:43
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double f[11][100]
bool isValid() const
Definition: HandleBase.h:75
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:37
Definition: DetId.h:18
DecomposeProduct< arg, typename Div::arg > D
Definition: Factorize.h:150
size_type size() const
map size
T const * product() const
Definition: Handle.h:81
const T & get() const
Definition: EventSetup.h:55
int layerName() const
layer id
T const * product() const
Definition: ESHandle.h:86
virtual void dqmBeginRun(const edm::Run &r, edm::EventSetup const &iSetup)
int pannelName() const
pannel id
tuple cout
Definition: gather_cfg.py:121
int diskName() const
disk id
std::pair< typename Association::data_type::first_type, double > match(Reference key, Association association, bool bestMatchByMaxValue)
Generic matching function.
Definition: Utils.h:10
const_iterator begin() const
first iterator over the map (read only)
volatile std::atomic< bool > shutdown_flag false
std::tuple< LocalPoint, LocalError, SiPixelRecHitQuality::QualWordType > ReturnType
bool isValid() const
Definition: ESHandle.h:47
T x() const
Definition: PV3DBase.h:62
Definition: vlib.h:208
const_iterator begin(bool update=false) const
Definition: Run.h:41