test
CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
SiTrackerGaussianSmearingRecHitConverter.cc
Go to the documentation of this file.
1 
9 // fast tracker recHits
14 
15 // SiTracker Gaussian Smearing
17 
18 // SiPixel Gaussian Smearing
20 
21 // SiStripGaussianSmearing
23 
24 // Geometry and magnetic field
32 
33 // Data Formats
36 
37 // Framework
44 
45 // Numbering scheme
48 
49 // Random engine
51 
52 // topology
53 
54 // STL
55 //#include <memory>
56 
57 // ROOT
58 #include <TFile.h>
59 #include <TH1F.h>
60 
61 
62 //#define FAMOS_DEBUG
63 
65  edm::ParameterSet const& conf)
66  : pset_(conf)
67 {
68  thePixelDataFile = 0;
75 
76 #ifdef FAMOS_DEBUG
77  std::cout << "SiTrackerGaussianSmearingRecHitConverter instantiated" << std::endl;
78 #endif
79 
80  produces<FastTrackerRecHitCollection>();
81  produces<FastTrackerRecHitRefCollection>("simHit2RecHitMap");
82 
83  //--- PSimHit Containers
84  // trackerContainers.clear();
85  // trackerContainers = conf.getParameter<std::vector<edm::InputTag> >("ROUList");
86  simHitLabel = conf.getParameter<edm::InputTag>("InputSimHits");
87  simHitToken = consumes<edm::PSimHitContainer>(simHitLabel);
88  //--- delta rays p cut [GeV/c] to filter PSimHits with p>
89  deltaRaysPCut = conf.getParameter<double>("DeltaRaysMomentumCut");
90 
91  //--- switch to have RecHit == PSimHit
92  trackingPSimHits = conf.getParameter<bool>("trackingPSimHits");
93  if(trackingPSimHits) std::cout << "### trackingPSimHits chosen " << trackingPSimHits << std::endl;
94 
95  // disable/enable dead channels
96  doDisableChannels = conf.getParameter<bool>("killDeadChannels");
97 
98  // Switch between old (ORCA) and new (CMSSW) pixel parameterization
99  useCMSSWPixelParameterization = conf.getParameter<bool>("UseCMSSWPixelParametrization");
100 #ifdef FAMOS_DEBUG
101  std::cout << (useCMSSWPixelParameterization? "CMSSW" : "ORCA") << " pixel parametrization chosen in config file." << std::endl;
102 #endif
103 
104  //
105  // TIB
106  localPositionResolution_TIB1x = conf.getParameter<double>("TIB1x");
107  localPositionResolution_TIB1y = conf.getParameter<double>("TIB1y");
108  localPositionResolution_TIB2x = conf.getParameter<double>("TIB2x");
109  localPositionResolution_TIB2y = conf.getParameter<double>("TIB2y");
110  localPositionResolution_TIB3x = conf.getParameter<double>("TIB3x");
111  localPositionResolution_TIB3y = conf.getParameter<double>("TIB3y");
112  localPositionResolution_TIB4x = conf.getParameter<double>("TIB4x");
113  localPositionResolution_TIB4y = conf.getParameter<double>("TIB4y");
114  //
115  // TID
116  localPositionResolution_TID1x = conf.getParameter<double>("TID1x");
117  localPositionResolution_TID1y = conf.getParameter<double>("TID1y");
118  localPositionResolution_TID2x = conf.getParameter<double>("TID2x");
119  localPositionResolution_TID2y = conf.getParameter<double>("TID2y");
120  localPositionResolution_TID3x = conf.getParameter<double>("TID3x");
121  localPositionResolution_TID3y = conf.getParameter<double>("TID3y");
122  //
123  // TOB
124  localPositionResolution_TOB1x = conf.getParameter<double>("TOB1x");
125  localPositionResolution_TOB1y = conf.getParameter<double>("TOB1y");
126  localPositionResolution_TOB2x = conf.getParameter<double>("TOB2x");
127  localPositionResolution_TOB2y = conf.getParameter<double>("TOB2y");
128  localPositionResolution_TOB3x = conf.getParameter<double>("TOB3x");
129  localPositionResolution_TOB3y = conf.getParameter<double>("TOB3y");
130  localPositionResolution_TOB4x = conf.getParameter<double>("TOB4x");
131  localPositionResolution_TOB4y = conf.getParameter<double>("TOB4y");
132  localPositionResolution_TOB5x = conf.getParameter<double>("TOB5x");
133  localPositionResolution_TOB5y = conf.getParameter<double>("TOB5y");
134  localPositionResolution_TOB6x = conf.getParameter<double>("TOB6x");
135  localPositionResolution_TOB6y = conf.getParameter<double>("TOB6y");
136  //
137  // TEC
138  localPositionResolution_TEC1x = conf.getParameter<double>("TEC1x");
139  localPositionResolution_TEC1y = conf.getParameter<double>("TEC1y");
140  localPositionResolution_TEC2x = conf.getParameter<double>("TEC2x");
141  localPositionResolution_TEC2y = conf.getParameter<double>("TEC2y");
142  localPositionResolution_TEC3x = conf.getParameter<double>("TEC3x");
143  localPositionResolution_TEC3y = conf.getParameter<double>("TEC3y");
144  localPositionResolution_TEC4x = conf.getParameter<double>("TEC4x");
145  localPositionResolution_TEC4y = conf.getParameter<double>("TEC4y");
146  localPositionResolution_TEC5x = conf.getParameter<double>("TEC5x");
147  localPositionResolution_TEC5y = conf.getParameter<double>("TEC5y");
148  localPositionResolution_TEC6x = conf.getParameter<double>("TEC6x");
149  localPositionResolution_TEC6y = conf.getParameter<double>("TEC6y");
150  localPositionResolution_TEC7x = conf.getParameter<double>("TEC7x");
151  localPositionResolution_TEC7y = conf.getParameter<double>("TEC7y");
152  //
153  localPositionResolution_z = 0.0001; // not to be changed, set to minimum (1 um), Kalman Filter will crash if errors are exactly 0, setting 1 um means 0
154  //
155 #ifdef FAMOS_DEBUG
156  std::cout << "RecHit local position error set to" << "\n"
157  << "\tTIB1\tx = " << localPositionResolution_TIB1x
158  << " cm\ty = " << localPositionResolution_TIB1y << " cm" << "\n"
159  << "\tTIB2\tx = " << localPositionResolution_TIB2x
160  << " cm\ty = " << localPositionResolution_TIB2y << " cm" << "\n"
161  << "\tTIB3\tx = " << localPositionResolution_TIB3x
162  << " cm\ty = " << localPositionResolution_TIB3y << " cm" << "\n"
163  << "\tTIB4\tx = " << localPositionResolution_TIB4x
164  << " cm\ty = " << localPositionResolution_TIB4y << " cm" << "\n"
165  << "\tTID1\tx = " << localPositionResolution_TID1x
166  << " cm\ty = " << localPositionResolution_TID1y << " cm" << "\n"
167  << "\tTID2\tx = " << localPositionResolution_TID2x
168  << " cm\ty = " << localPositionResolution_TID2y << " cm" << "\n"
169  << "\tTID3\tx = " << localPositionResolution_TID3x
170  << " cm\ty = " << localPositionResolution_TID3y << " cm" << "\n"
171  << "\tTOB1\tx = " << localPositionResolution_TOB1x
172  << " cm\ty = " << localPositionResolution_TOB1y << " cm" << "\n"
173  << "\tTOB2\tx = " << localPositionResolution_TOB2x
174  << " cm\ty = " << localPositionResolution_TOB2y << " cm" << "\n"
175  << "\tTOB3\tx = " << localPositionResolution_TOB3x
176  << " cm\ty = " << localPositionResolution_TOB3y << " cm" << "\n"
177  << "\tTOB4\tx = " << localPositionResolution_TOB4x
178  << " cm\ty = " << localPositionResolution_TOB4y << " cm" << "\n"
179  << "\tTOB5\tx = " << localPositionResolution_TOB5x
180  << " cm\ty = " << localPositionResolution_TOB5y << " cm" << "\n"
181  << "\tTOB6\tx = " << localPositionResolution_TOB6x
182  << " cm\ty = " << localPositionResolution_TOB6y << " cm" << "\n"
183  << "\tTEC1\tx = " << localPositionResolution_TEC1x
184  << " cm\ty = " << localPositionResolution_TEC1y << " cm" << "\n"
185  << "\tTEC2\tx = " << localPositionResolution_TEC2x
186  << " cm\ty = " << localPositionResolution_TEC2y << " cm" << "\n"
187  << "\tTEC3\tx = " << localPositionResolution_TEC3x
188  << " cm\ty = " << localPositionResolution_TEC3y << " cm" << "\n"
189  << "\tTEC4\tx = " << localPositionResolution_TEC4x
190  << " cm\ty = " << localPositionResolution_TEC4y << " cm" << "\n"
191  << "\tTEC5\tx = " << localPositionResolution_TEC5x
192  << " cm\ty = " << localPositionResolution_TEC5y << " cm" << "\n"
193  << "\tTEC6\tx = " << localPositionResolution_TEC6x
194  << " cm\ty = " << localPositionResolution_TEC6y << " cm" << "\n"
195  << "\tTEC7\tx = " << localPositionResolution_TEC7x
196  << " cm\ty = " << localPositionResolution_TEC7y << " cm" << "\n"
197  << "\tAll:\tz = " << localPositionResolution_z << " cm"
198  << std::endl;
199 #endif
200 
201  //--- Number of histograms for alpha/beta barrel/forward multiplicity
202  if(useCMSSWPixelParameterization) {
203  nAlphaBarrel = conf.getParameter<int>("AlphaBarrelMultiplicityNew");
204  nBetaBarrel = conf.getParameter<int>("BetaBarrelMultiplicityNew");
205  nAlphaForward = conf.getParameter<int>("AlphaForwardMultiplicityNew");
206  nBetaForward = conf.getParameter<int>("BetaForwardMultiplicityNew");
207  } else {
208  nAlphaBarrel = conf.getParameter<int>("AlphaBarrelMultiplicity");
209  nBetaBarrel = conf.getParameter<int>("BetaBarrelMultiplicity");
210  nAlphaForward = conf.getParameter<int>("AlphaForwardMultiplicity");
211  nBetaForward = conf.getParameter<int>("BetaForwardMultiplicity");
212  }
213 #ifdef FAMOS_DEBUG
214  std::cout << "Pixel maximum multiplicity set to "
215  << "\nBarrel" << "\talpha " << nAlphaBarrel
216  << "\tbeta " << nBetaBarrel
217  << "\nForward" << "\talpha " << nAlphaForward
218  << "\tbeta " << nBetaForward
219  << std::endl;
220 #endif
221 
222  // Resolution Barrel
223  if(useCMSSWPixelParameterization) {
224  resAlphaBarrel_binMin = conf.getParameter<double>("AlphaBarrel_BinMinNew" );
225  resAlphaBarrel_binWidth = conf.getParameter<double>("AlphaBarrel_BinWidthNew");
226  resAlphaBarrel_binN = conf.getParameter<int>( "AlphaBarrel_BinNNew" );
227  resBetaBarrel_binMin = conf.getParameter<double>("BetaBarrel_BinMinNew" );
228  resBetaBarrel_binWidth = conf.getParameter<double>("BetaBarrel_BinWidthNew" );
229  resBetaBarrel_binN = conf.getParameter<int>( "BetaBarrel_BinNNew" );
230  } else {
231  resAlphaBarrel_binMin = conf.getParameter<double>("AlphaBarrel_BinMin" );
232  resAlphaBarrel_binWidth = conf.getParameter<double>("AlphaBarrel_BinWidth");
233  resAlphaBarrel_binN = conf.getParameter<int>( "AlphaBarrel_BinN" );
234  resBetaBarrel_binMin = conf.getParameter<double>("BetaBarrel_BinMin" );
235  resBetaBarrel_binWidth = conf.getParameter<double>("BetaBarrel_BinWidth" );
236  resBetaBarrel_binN = conf.getParameter<int>( "BetaBarrel_BinN" );
237  }
238 
239  // Resolution Forward
240  if(useCMSSWPixelParameterization) {
241  resAlphaForward_binMin = conf.getParameter<double>("AlphaForward_BinMinNew" );
242  resAlphaForward_binWidth = conf.getParameter<double>("AlphaForward_BinWidthNew" );
243  resAlphaForward_binN = conf.getParameter<int>( "AlphaForward_BinNNew" );
244  resBetaForward_binMin = conf.getParameter<double>("BetaForward_BinMinNew" );
245  resBetaForward_binWidth = conf.getParameter<double>("BetaForward_BinWidthNew" );
246  resBetaForward_binN = conf.getParameter<int>( "BetaForward_BinNNew" );
247  } else {
248  resAlphaForward_binMin = conf.getParameter<double>("AlphaForward_BinMin" );
249  resAlphaForward_binWidth = conf.getParameter<double>("AlphaForward_BinWidth" );
250  resAlphaForward_binN = conf.getParameter<int>( "AlphaForward_BinN" );
251  resBetaForward_binMin = conf.getParameter<double>("BetaForward_BinMin" );
252  resBetaForward_binWidth = conf.getParameter<double>("BetaForward_BinWidth" );
253  resBetaForward_binN = conf.getParameter<int>( "BetaForward_BinN" );
254  }
255 
256  // Hit Finding Probability
257  theHitFindingProbability_PXB = conf.getParameter<double>("HitFindingProbability_PXB" );
258  theHitFindingProbability_PXF = conf.getParameter<double>("HitFindingProbability_PXF" );
259  theHitFindingProbability_TIB1 = conf.getParameter<double>("HitFindingProbability_TIB1");
260  theHitFindingProbability_TIB2 = conf.getParameter<double>("HitFindingProbability_TIB2");
261  theHitFindingProbability_TIB3 = conf.getParameter<double>("HitFindingProbability_TIB3");
262  theHitFindingProbability_TIB4 = conf.getParameter<double>("HitFindingProbability_TIB4");
263  theHitFindingProbability_TID1 = conf.getParameter<double>("HitFindingProbability_TID1");
264  theHitFindingProbability_TID2 = conf.getParameter<double>("HitFindingProbability_TID2");
265  theHitFindingProbability_TID3 = conf.getParameter<double>("HitFindingProbability_TID3");
266  theHitFindingProbability_TOB1 = conf.getParameter<double>("HitFindingProbability_TOB1");
267  theHitFindingProbability_TOB2 = conf.getParameter<double>("HitFindingProbability_TOB2");
268  theHitFindingProbability_TOB3 = conf.getParameter<double>("HitFindingProbability_TOB3");
269  theHitFindingProbability_TOB4 = conf.getParameter<double>("HitFindingProbability_TOB4");
270  theHitFindingProbability_TOB5 = conf.getParameter<double>("HitFindingProbability_TOB5");
271  theHitFindingProbability_TOB6 = conf.getParameter<double>("HitFindingProbability_TOB6");
272  theHitFindingProbability_TEC1 = conf.getParameter<double>("HitFindingProbability_TEC1");
273  theHitFindingProbability_TEC2 = conf.getParameter<double>("HitFindingProbability_TEC2");
274  theHitFindingProbability_TEC3 = conf.getParameter<double>("HitFindingProbability_TEC3");
275  theHitFindingProbability_TEC4 = conf.getParameter<double>("HitFindingProbability_TEC4");
276  theHitFindingProbability_TEC5 = conf.getParameter<double>("HitFindingProbability_TEC5");
277  theHitFindingProbability_TEC6 = conf.getParameter<double>("HitFindingProbability_TEC6");
278  theHitFindingProbability_TEC7 = conf.getParameter<double>("HitFindingProbability_TEC7");
279  //
280 #ifdef FAMOS_DEBUG
281  std::cout << "RecHit finding probability set to" << "\n"
282  << "\tPXB = " << theHitFindingProbability_PXB << "\n"
283  << "\tPXF = " << theHitFindingProbability_PXF << "\n"
284  << "\tTIB1 = " << theHitFindingProbability_TIB1 << "\n"
285  << "\tTIB2 = " << theHitFindingProbability_TIB2 << "\n"
286  << "\tTIB3 = " << theHitFindingProbability_TIB3 << "\n"
287  << "\tTIB4 = " << theHitFindingProbability_TIB4 << "\n"
288  << "\tTID1 = " << theHitFindingProbability_TID1 << "\n"
289  << "\tTID2 = " << theHitFindingProbability_TID2 << "\n"
290  << "\tTID3 = " << theHitFindingProbability_TID3 << "\n"
291  << "\tTOB1 = " << theHitFindingProbability_TOB1 << "\n"
292  << "\tTOB2 = " << theHitFindingProbability_TOB2 << "\n"
293  << "\tTOB3 = " << theHitFindingProbability_TOB3 << "\n"
294  << "\tTOB4 = " << theHitFindingProbability_TOB4 << "\n"
295  << "\tTOB5 = " << theHitFindingProbability_TOB5 << "\n"
296  << "\tTOB6 = " << theHitFindingProbability_TOB6 << "\n"
297  << "\tTEC1 = " << theHitFindingProbability_TEC1 << "\n"
298  << "\tTEC2 = " << theHitFindingProbability_TEC2 << "\n"
299  << "\tTEC3 = " << theHitFindingProbability_TEC3 << "\n"
300  << "\tTEC4 = " << theHitFindingProbability_TEC4 << "\n"
301  << "\tTEC5 = " << theHitFindingProbability_TEC5 << "\n"
302  << "\tTEC6 = " << theHitFindingProbability_TEC6 << "\n"
303  << "\tTEC7 = " << theHitFindingProbability_TEC7 << "\n"
304  << std::endl;
305 #endif
306 
307  // Initialize the si strip error parametrization
310 
311  // Initialization of pixel parameterization posponed to beginRun(), since it depends on the magnetic field
312 
313 }
314 
315 
317  // load multiplicity cumulative probabilities
318  // root files
319  thePixelDataFile = new TFile ( edm::FileInPath( thePixelMultiplicityFileName ).fullPath().c_str() , "READ" );
322  //
323 
324  // alpha barrel
326  nAlphaBarrel ,
327  std::string("hist_alpha_barrel") ,
329  //
330  // beta barrel
332  nBetaBarrel ,
333  std::string("hist_beta_barrel") ,
335  //
336  // alpha forward
338  nAlphaForward ,
339  std::string("hist_alpha_forward") ,
341  //
342  // beta forward
344  nBetaForward ,
345  std::string("hist_beta_forward") ,
347 
348  // Load also big pixel data if CMSSW parametrization is on
349  // They are pushed back into the vectors after the normal pixels data:
350  // [0, ..., (size/2)-1] -> Normal pixels
351  // [size/2, ..., size-1] -> Big pixels
353  // alpha barrel
355  nAlphaBarrel ,
356  std::string("hist_alpha_barrel_big") ,
358  true );
359  //
360  // beta barrel
362  nBetaBarrel ,
363  std::string("hist_beta_barrel_big") ,
365  true );
366  //
367  // alpha forward
369  nAlphaForward ,
370  std::string("hist_alpha_forward_big") ,
372  true );
373  //
374  // beta forward
376  nBetaForward ,
377  std::string("hist_beta_forward_big") ,
379  true );
380  }
381  //
382 }
383 
385  TFile* pixelDataFile,
386  unsigned int nMultiplicity,
387  std::string histName,
388  std::vector<TH1F*>& theMultiplicityCumulativeProbabilities,
389  bool bigPixels)
390 {
391 
392  std::string histName_i = histName + "_%u"; // needed to open histograms with a for
393  if(!bigPixels)
394  theMultiplicityCumulativeProbabilities.clear();
395  //
396  // What's this vector? Not needed - MG
397 // std::vector<double> mult; // vector with fixed multiplicity
398  for(unsigned int i = 0; i<nMultiplicity; ++i) {
399  TH1F addHist = *((TH1F*) pixelDataFile->Get( Form( histName_i.c_str() ,i+1 )));
400  if(i==0) {
401  theMultiplicityCumulativeProbabilities.push_back( new TH1F(addHist) );
402  } else {
403  TH1F sumHist;
404  if(bigPixels)
405  sumHist = *(theMultiplicityCumulativeProbabilities[nMultiplicity+i-1]);
406  else
407  sumHist = *(theMultiplicityCumulativeProbabilities[i-1]);
408  sumHist.Add(&addHist);
409  theMultiplicityCumulativeProbabilities.push_back( new TH1F(sumHist) );
410  }
411  }
412 
413  // Logger
414 #ifdef FAMOS_DEBUG
415  const unsigned int maxMult = theMultiplicityCumulativeProbabilities.size();
416  unsigned int iMult, multSize;
418  if(bigPixels) {
419  iMult = maxMult / 2;
420  multSize = maxMult ;
421  } else {
422  iMult = 0;
423  multSize = maxMult;
424  }
425  } else {
426  iMult = 0;
427  multSize = maxMult ;
428  }
429  std::cout << " Multiplicity cumulated probability " << histName << std::endl;
430  for(/* void */; iMult<multSize; ++iMult) {
431  for(int iBin = 1; iBin<=theMultiplicityCumulativeProbabilities[iMult]->GetNbinsX(); ++iBin) {
432  std::cout
433  << " Multiplicity " << iMult+1
434  << " bin " << iBin
435  << " low edge = "
436  << theMultiplicityCumulativeProbabilities[iMult]->GetBinLowEdge(iBin)
437  << " prob = "
438  << (theMultiplicityCumulativeProbabilities[iMult])->GetBinContent(iBin)
439  // remember in ROOT bin starts from 1 (0 underflow, nBin+1 overflow)
440  << std::endl;
441  }
442  }
443 #endif
444 
445 }
446 
447 // Destructor
453 
460 
462 }
463 
464 void
466 {
467 
468  // Initialize the Tracker Geometry
469  edm::ESHandle<TrackerGeometry> theGeometry;
470  es.get<TrackerDigiGeometryRecord> ().get (theGeometry);
471  geometry = &(*theGeometry);
472 
473  edm::ESHandle<TrackerGeometry> theMisAlignedGeometry;
474  es.get<TrackerDigiGeometryRecord>().get("MisAligned",theMisAlignedGeometry);
475  misAlignedGeometry = &(*theMisAlignedGeometry);
476 
477  const MagneticField* magfield;
479  es.get<IdealMagneticFieldRecord>().get(magField);
480  magfield=&(*magField);
481  GlobalPoint center(0.0, 0.0, 0.0);
482  double magFieldAtCenter = magfield->inTesla(center).mag();
483 
484  // For new parameterization: select multiplicity and resolution files according to magnetic field
486  if(magFieldAtCenter > 3.9) {
487  thePixelMultiplicityFileName = pset_.getParameter<std::string>( "PixelMultiplicityFile40T");
488  thePixelBarrelResolutionFileName = pset_.getParameter<std::string>( "PixelBarrelResolutionFile40T");
489  thePixelForwardResolutionFileName = pset_.getParameter<std::string>( "PixelForwardResolutionFile40T");
490  } else {
491  thePixelMultiplicityFileName = pset_.getParameter<std::string>( "PixelMultiplicityFile38T");
492  thePixelBarrelResolutionFileName = pset_.getParameter<std::string>( "PixelBarrelResolutionFile38T");
493  thePixelForwardResolutionFileName = pset_.getParameter<std::string>( "PixelForwardResolutionFile38T");
494  }
495  } else {
496  thePixelMultiplicityFileName = pset_.getParameter<std::string>( "PixelMultiplicityFile" );
497  thePixelBarrelResolutionFileName = pset_.getParameter<std::string>( "PixelBarrelResolutionFile");
498  thePixelForwardResolutionFileName = pset_.getParameter<std::string>( "PixelForwardResolutionFile");
499  }
500 
501 
502  // Reading the list of dead pixel modules from DB:
503  edm::ESHandle<SiPixelQuality> siPixelBadModule;
504  es.get<SiPixelQualityRcd>().get(siPixelBadModule);
506  if (doDisableChannels) {
507  disabledModules = new std::vector<SiPixelQuality::disabledModuleType> ( siPixelBadModule->getBadComponentList() );
509  size_t numberOfRecoverableModules = 0;
510  for (size_t id=0;id<numberOfDisabledModules;id++) {
512  // errortype "whole" = int 0 in DB //
513  // errortype "tbmA" = int 1 in DB //
514  // errortype "tbmB" = int 2 in DB //
515  // errortype "none" = int 3 in DB //
517  if ( (*disabledModules)[id-numberOfRecoverableModules].errorType != 0 ){
518  // Disable only the modules totally in error:
519  disabledModules->erase(disabledModules->begin()+id-numberOfRecoverableModules);
520  numberOfRecoverableModules++;
521  }
522  }
523  numberOfDisabledModules = disabledModules->size();
524  }
525 
526 
527 
528 #ifdef FAMOS_DEBUG
529  std::cout << "Pixel multiplicity data are taken from file " << thePixelMultiplicityFileName << std::endl;
530 
531  std::cout << "Pixel maximum multiplicity set to "
532  << "\nBarrel" << "\talpha " << nAlphaBarrel
533  << "\tbeta " << nBetaBarrel
534  << "\nForward" << "\talpha " << nAlphaForward
535  << "\tbeta " << nBetaForward
536  << std::endl;
537 
538  std::cout << "Barrel Pixel resolution data are taken from file "
540  << "Alpha bin min = " << resAlphaBarrel_binMin
541  << "\twidth = " << resAlphaBarrel_binWidth
542  << "\tbins = " << resAlphaBarrel_binN
543  << "\n"
544  << " Beta bin min = " << resBetaBarrel_binMin
545  << "\twidth = " << resBetaBarrel_binWidth
546  << "\tbins = " << resBetaBarrel_binN
547  << std::endl;
548 
549  std::cout << "Forward Pixel resolution data are taken from file "
551  << "Alpha bin min = " << resAlphaForward_binMin
552  << "\twidth = " << resAlphaForward_binWidth
553  << "\tbins = " << resAlphaForward_binN
554  << "\n"
555  << " Beta bin min = " << resBetaForward_binMin
556  << "\twidth = " << resBetaForward_binWidth
557  << "\tbins = " << resBetaForward_binN
558  << std::endl;
559 #endif
560  //
561 
562  //
563  // load pixel data
564  loadPixelData();
565  //
566 
567  // Initialize and open relevant files for the pixel barrel error parametrization
570  pset_,
572  // Initialize and open relevant files for the pixel forward error parametrization
575  pset_,
577 }
578 
580 {
582 
583  //Retrieve tracker topology from geometry
585  es.get<TrackerTopologyRcd>().get(tTopoHand);
586  const TrackerTopology *tTopo=tTopoHand.product();
587 
588  // input: simHits
589  edm::Handle<edm::PSimHitContainer> allTrackerHits_handle;
590  e.getByToken(simHitToken,allTrackerHits_handle);
591  const edm::PSimHitContainer& allTrackerHits=*allTrackerHits_handle;
592 
593  // output: recHits
594  std::unique_ptr<FastTrackerRecHitCollection> output_recHits(new FastTrackerRecHitCollection);
595  output_recHits->reserve(allTrackerHits.size());
596 
597  // output: map simHit -> recHit
598  // by default, each simHit is associated to a null ref
600  std::unique_ptr<FastTrackerRecHitRefCollection> output_recHitRefs(new FastTrackerRecHitRefCollection(allTrackerHits.size(),FastTrackerRecHitRef()));
601 
602  // loop on PSimHits
603  for (unsigned simHitCounter = 0;simHitCounter<allTrackerHits.size();++simHitCounter) {
604 
605  const PSimHit & simHit = allTrackerHits[simHitCounter];
606 
607  // skip hits on bad modules
608  DetId det(simHit.detUnitId());
609  bool isBad = false;
610  unsigned int geoId = det.rawId();
611  for (size_t id=0;id<numberOfDisabledModules;id++) {
612  if(geoId==(*disabledModules)[id].DetID){
613  isBad = true;
614  break;
615  }
616  }
617  if(isBad) continue;
618 
619  // smear
622  bool isCreated = smear(simHit, position, error,tTopo, &random);
623  unsigned int subdet = det.subdetId();
624 
625  if(isCreated) {
626 
627  // Inflate errors in case of geometry misaligniment
628  // (still needed! what done in constructor of BaseTrackerRecHit is not effective ad geometry is not missaligned)
629  auto theMADet = misAlignedGeometry->idToDet(det);
630  auto const & lape = theMADet->localAlignmentError();
631  if ( lape.valid() )
632  error = LocalError ( error.xx()+lape.xx(),
633  error.xy()+lape.xy(),
634  error.yy()+lape.yy() );
635 
636  // insert rechit in rechit collection
637  std::auto_ptr<FastSingleTrackerRecHit> recHit (new FastSingleTrackerRecHit(position, error,
638  *geometry->idToDetUnit(det),
639  subdet > 2
642  recHit->addSimTrackId(simHit.trackId());
643  recHit->setId(simHitCounter);
644  output_recHits->push_back(recHit);
645 
646 
647  // update map simHit->recHit
648  (*output_recHitRefs)[simHitCounter] = FastTrackerRecHitRef(output_recHits_refProd,output_recHits->size()-1);
649  } // end if(isCreated)
650 
651  } // end loop on PSimHits
652 
653 
654  // put products in event
655  e.put(std::move(output_recHits));
656  e.put(std::move(output_recHitRefs),"simHit2RecHitMap");
657 
658 }
659 
660 
661 
664  LocalError& error,
665  const TrackerTopology *tTopo,
667 {
668 
669  // A few caracteritics of the detid the SimHit belongs to.
670  unsigned int subdet = DetId(simHit.detUnitId()).subdetId();
671  unsigned int detid = DetId(simHit.detUnitId()).rawId();
672  const GeomDetUnit* theDetUnit = geometry->idToDetUnit((DetId)simHit.detUnitId());
673  const BoundPlane& theDetPlane = theDetUnit->surface();
674  const Bounds& theBounds = theDetPlane.bounds();
675  double boundX = theBounds.width()/2.;
676  double boundY = theBounds.length()/2.;
677 
678 #ifdef FAMOS_DEBUG
679  std::cout << "\tSubdetector " << subdet
680  << " rawid " << detid
681  << std::endl;
682 #endif
683  if(trackingPSimHits) {
684  // z is fixed for all detectors, in case of errors resolution is fixed also for x and y to 1 um (zero)
685  // The Matrix is the Covariance Matrix, sigma^2 on diagonal!!!
687  0.0 ,
688  localPositionResolution_z * localPositionResolution_z );
689  //
690  // starting from PSimHit local position
691  position = simHit.localPosition();
692 #ifdef FAMOS_DEBUG
693  std::cout << " Tracking PSimHit position set to " << position;
694 #endif
695  return true; // RecHit == PSimHit with 100% hit finding efficiency
696  }
697  //
698 
699  // hit finding probability --> RecHit will be created if and only if hitFindingProbability <= theHitFindingProbability_###
700  double hitFindingProbability = random->flatShoot();
701 #ifdef FAMOS_DEBUG
702  std::cout << " Hit finding probability draw: " << hitFindingProbability << std::endl;;
703 #endif
704 
705  switch (subdet) {
706  // Pixel Barrel
707  case 1:
708  {
709 #ifdef FAMOS_DEBUG
710 
711  unsigned int theLayer = tTopo->pxbLayer(detid);
712  std::cout << "\tPixel Barrel Layer " << theLayer << std::endl;
713 #endif
714  if( hitFindingProbability > theHitFindingProbability_PXB ) return false;
715  // Hit smearing
716  const PixelGeomDetUnit* pixelDetUnit = dynamic_cast<const PixelGeomDetUnit*>(theDetUnit);
717  thePixelBarrelParametrization->smearHit(simHit, pixelDetUnit, boundX, boundY, random);
720  return true;
721  break;
722  }
723  // Pixel Forward
724  case 2:
725  {
726 #ifdef FAMOS_DEBUG
727 
728  unsigned int theDisk = tTopo->pxfDisk(detid);
729  std::cout << "\tPixel Forward Disk " << theDisk << std::endl;
730 #endif
731  if( hitFindingProbability > theHitFindingProbability_PXF ) return false;
732  // Hit smearing
733  const PixelGeomDetUnit* pixelDetUnit = dynamic_cast<const PixelGeomDetUnit*>(theDetUnit);
734  thePixelEndcapParametrization->smearHit(simHit, pixelDetUnit, boundX, boundY, random);
737  return true;
738  break;
739  }
740  // TIB
741  case 3:
742  {
743 
744  unsigned int theLayer = tTopo->tibLayer(detid);
745 #ifdef FAMOS_DEBUG
746  std::cout << "\tTIB Layer " << theLayer << std::endl;
747 #endif
748  //
749  double resolutionX, resolutionY, resolutionZ;
750  resolutionZ = localPositionResolution_z;
751 
752  switch (theLayer) {
753  case 1:
754  {
755  resolutionX = localPositionResolution_TIB1x;
756  resolutionY = localPositionResolution_TIB1y;
757  if( hitFindingProbability > theHitFindingProbability_TIB1 ) return false;
758  break;
759  }
760  case 2:
761  {
762  resolutionX = localPositionResolution_TIB2x;
763  resolutionY = localPositionResolution_TIB2y;
764  if( hitFindingProbability > theHitFindingProbability_TIB2 ) return false;
765  break;
766  }
767  case 3:
768  {
769  resolutionX = localPositionResolution_TIB3x;
770  resolutionY = localPositionResolution_TIB3y;
771  if( hitFindingProbability > theHitFindingProbability_TIB3 ) return false;
772  break;
773  }
774  case 4:
775  {
776  resolutionX = localPositionResolution_TIB4x;
777  resolutionY = localPositionResolution_TIB4y;
778  if( hitFindingProbability > theHitFindingProbability_TIB4 ) return false;
779  break;
780  }
781  default:
782  {
783  edm::LogError ("SiTrackerGaussianSmearingRecHits")
784  << "\tTIB Layer not valid " << theLayer << std::endl;
785  return false;
786  break;
787  }
788  }
789 
790  // Gaussian smearing
791  theSiStripErrorParametrization->smearHit(simHit, resolutionX, resolutionY, resolutionZ, boundX, boundY, random);
794  return true;
795  break;
796  } // TIB
797 
798  // TID
799  case 4:
800  {
801 
802  unsigned int theRing = tTopo->tidRing(detid);
803  double resolutionFactorY =
804  1. - simHit.localPosition().y() / theDetPlane.position().perp();
805 
806 #ifdef FAMOS_DEBUG
807  std::cout << "\tTID Ring " << theRing << std::endl;
808 #endif
809  double resolutionX, resolutionY, resolutionZ;
810  resolutionZ = localPositionResolution_z;
811 
812  switch (theRing) {
813  case 1:
814  {
815  resolutionX = localPositionResolution_TID1x * resolutionFactorY;
816  resolutionY = localPositionResolution_TID1y;
817  if( hitFindingProbability > theHitFindingProbability_TID1 ) return false;
818  break;
819  }
820  case 2:
821  {
822  resolutionX = localPositionResolution_TID2x * resolutionFactorY;
823  resolutionY = localPositionResolution_TID2y;
824  if( hitFindingProbability > theHitFindingProbability_TID2 ) return false;
825  break;
826  }
827  case 3:
828  {
829  resolutionX = localPositionResolution_TID3x * resolutionFactorY;
830  resolutionY = localPositionResolution_TID3y;
831  if( hitFindingProbability > theHitFindingProbability_TID3 ) return false;
832  break;
833  }
834  default:
835  {
836  edm::LogError ("SiTrackerGaussianSmearingRecHits")
837  << "\tTID Ring not valid " << theRing << std::endl;
838  return false;
839  break;
840  }
841  }
842 
843  boundX *= resolutionFactorY;
844 
845  theSiStripErrorParametrization->smearHit(simHit, resolutionX, resolutionY, resolutionZ, boundX, boundY, random);
848  return true;
849  break;
850  } // TID
851 
852  // TOB
853  case 5:
854  {
855 
856  unsigned int theLayer = tTopo->tobLayer(detid);
857 #ifdef FAMOS_DEBUG
858  std::cout << "\tTOB Layer " << theLayer << std::endl;
859 #endif
860  double resolutionX, resolutionY, resolutionZ;
861  resolutionZ = localPositionResolution_z;
862 
863  switch (theLayer) {
864  case 1:
865  {
866  resolutionX = localPositionResolution_TOB1x;
867  resolutionY = localPositionResolution_TOB1y;
868  if( hitFindingProbability > theHitFindingProbability_TOB1 ) return false;
869  break;
870  }
871  case 2:
872  {
873  resolutionX = localPositionResolution_TOB2x;
874  resolutionY = localPositionResolution_TOB2y;
875  if( hitFindingProbability > theHitFindingProbability_TOB2 ) return false;
876  break;
877  }
878  case 3:
879  {
880  resolutionX = localPositionResolution_TOB3x;
881  resolutionY = localPositionResolution_TOB3y;
882  if( hitFindingProbability > theHitFindingProbability_TOB3 ) return false;
883  break;
884  }
885  case 4:
886  {
887  resolutionX = localPositionResolution_TOB4x;
888  resolutionY = localPositionResolution_TOB4y;
889  if( hitFindingProbability > theHitFindingProbability_TOB4 ) return false;
890  break;
891  }
892  case 5:
893  {
894  resolutionX = localPositionResolution_TOB5x;
895  resolutionY = localPositionResolution_TOB5y;
896  if( hitFindingProbability > theHitFindingProbability_TOB5 ) return false;
897  break;
898  }
899  case 6:
900  {
901  resolutionX = localPositionResolution_TOB6x;
902  resolutionY = localPositionResolution_TOB6y;
903  if( hitFindingProbability > theHitFindingProbability_TOB6 ) return false;
904  break;
905  }
906  default:
907  {
908  edm::LogError ("SiTrackerGaussianSmearingRecHits")
909  << "\tTOB Layer not valid " << theLayer << std::endl;
910  return false;
911  break;
912  }
913  }
914  theSiStripErrorParametrization->smearHit(simHit, resolutionX, resolutionY, resolutionZ, boundX, boundY, random);
917  return true;
918  break;
919  } // TOB
920 
921  // TEC
922  case 6:
923  {
924 
925  unsigned int theRing = tTopo->tecRing(detid);
926  double resolutionFactorY =
927  1. - simHit.localPosition().y() / theDetPlane.position().perp();
928 
929 #ifdef FAMOS_DEBUG
930  std::cout << "\tTEC Ring " << theRing << std::endl;
931 #endif
932  double resolutionX, resolutionY, resolutionZ;
934 
935  switch (theRing) {
936  case 1:
937  {
938  resolutionX = localPositionResolution_TEC1x * resolutionFactorY;
939  resolutionY = localPositionResolution_TEC1y;
940  if( hitFindingProbability > theHitFindingProbability_TEC1 ) return false;
941  break;
942  }
943  case 2:
944  {
945  resolutionX = localPositionResolution_TEC2x * resolutionFactorY;
946  resolutionY = localPositionResolution_TEC2y;
947  if( hitFindingProbability > theHitFindingProbability_TEC2 ) return false;
948  break;
949  }
950  case 3:
951  {
952  resolutionX = localPositionResolution_TEC3x * resolutionFactorY;
953  resolutionY = localPositionResolution_TEC3y;
954  if( hitFindingProbability > theHitFindingProbability_TEC3 ) return false;
955  break;
956  }
957  case 4:
958  {
959  resolutionX = localPositionResolution_TEC4x * resolutionFactorY;
960  resolutionY = localPositionResolution_TEC4y;
961  if( hitFindingProbability > theHitFindingProbability_TEC4 ) return false;
962  break;
963  }
964  case 5:
965  {
966  resolutionX = localPositionResolution_TEC5x * resolutionFactorY;
967  resolutionY = localPositionResolution_TEC5y;
968  if( hitFindingProbability > theHitFindingProbability_TEC5 ) return false;
969  break;
970  }
971  case 6:
972  {
973  resolutionX = localPositionResolution_TEC6x * resolutionFactorY;
974  resolutionY = localPositionResolution_TEC6y;
975  if( hitFindingProbability > theHitFindingProbability_TEC6 ) return false;
976  break;
977  }
978  case 7:
979  {
980  resolutionX = localPositionResolution_TEC7x * resolutionFactorY;
981  resolutionY = localPositionResolution_TEC7y;
982  if( hitFindingProbability > theHitFindingProbability_TEC7 ) return false;
983  break;
984  }
985  default:
986  {
987  edm::LogError ("SiTrackerGaussianSmearingRecHits")
988  << "\tTEC Ring not valid " << theRing << std::endl;
989  return false;
990  break;
991  }
992  }
993 
994  boundX *= resolutionFactorY;
995  theSiStripErrorParametrization->smearHit(simHit, resolutionX, resolutionY, resolutionZ, boundX, boundY, random);
998  return true;
999  break;
1000  } // TEC
1001 
1002  default:
1003  {
1004  edm::LogError ("SiTrackerGaussianSmearingRecHits") << "\tTracker subdetector not valid " << subdet << std::endl;
1005  return false;
1006  break;
1007  }
1008 
1009  } // subdetector case
1010  //
1011 }
T getParameter(std::string const &) const
virtual const TrackerGeomDet * idToDetUnit(DetId) const
Return the pointer to the GeomDetUnit corresponding to a given DetId.
int i
Definition: DBlmapReader.cc:9
float xx() const
Definition: LocalError.h:24
SiStripGaussianSmearingRecHitConverterAlgorithm * theSiStripErrorParametrization
virtual float length() const =0
unsigned int tibLayer(const DetId &id) const
unsigned int tidRing(const DetId &id) const
SiPixelGaussianSmearingRecHitConverterAlgorithm * thePixelEndcapParametrization
double flatShoot(double xmin=0.0, double xmax=1.0) const
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:464
virtual GlobalVector inTesla(const GlobalPoint &gp) const =0
Field value ad specified global point, in Tesla.
unsigned int pxfDisk(const DetId &id) const
unsigned int tecRing(const DetId &id) const
ring id
tuple magfield
Definition: HLT_ES_cff.py:2311
T y() const
Definition: PV3DBase.h:63
TRandom random
Definition: MVATrainer.cc:138
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:40
edm::Ref< FastTrackerRecHitCollection > FastTrackerRecHitRef
uint32_t rawId() const
get the raw id
Definition: DetId.h:43
float xy() const
Definition: LocalError.h:25
T mag() const
Definition: PV3DBase.h:67
float yy() const
Definition: LocalError.h:26
Local3DPoint localPosition() const
Definition: PSimHit.h:44
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:120
def move
Definition: eostools.py:510
bool smear(const PSimHit &simHit, Local3DPoint &position, LocalError &error, const TrackerTopology *tTopo, RandomEngineAndDistribution const *)
SiPixelGaussianSmearingRecHitConverterAlgorithm * thePixelBarrelParametrization
std::vector< SiPixelQuality::disabledModuleType > * disabledModules
RefProd< PROD > getRefBeforePut()
Definition: Event.h:140
tuple conf
Definition: dbtoconf.py:185
unsigned int pxbLayer(const DetId &id) const
Definition: DetId.h:18
virtual void produce(edm::Event &e, const edm::EventSetup &c) override
virtual void beginRun(edm::Run const &run, const edm::EventSetup &es) override
const T & get() const
Definition: EventSetup.h:56
void smearHit(const PSimHit &simHit, const PixelGeomDetUnit *detUnit, const double boundX, const double boundY, RandomEngineAndDistribution const *)
std::vector< FastTrackerRecHitRef > FastTrackerRecHitRefCollection
static int position[264][3]
Definition: ReadPGInfo.cc:509
StreamID streamID() const
Definition: Event.h:79
unsigned int trackId() const
Definition: PSimHit.h:102
tuple cout
Definition: gather_cfg.py:121
std::vector< PSimHit > PSimHitContainer
Definition: Bounds.h:22
virtual float width() const =0
LocalError const & localAlignmentError() const
Return local alligment error.
unsigned int detUnitId() const
Definition: PSimHit.h:93
unsigned int tobLayer(const DetId &id) const
Definition: Run.h:43
void smearHit(const PSimHit &simHit, double localPositionResolutionX, double localPositionResolutionY, double localPositionResolutionZ, double boundX, double boundY, RandomEngineAndDistribution const *)
virtual const TrackerGeomDet * idToDet(DetId) const