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 
8 //PAT
9 //FastTrackerCluster
11 
12 // SiTracker Gaussian Smearing
14 
15 // SiPixel Gaussian Smearing
17 
18 // SiStripGaussianSmearing
20 
21 // Geometry and magnetic field
29 
30 // Data Formats
33 
34 // Framework
43 
44 // Numbering scheme
49 
50 //For Pileup events
52 
53 // Random engine
55 
56 // topology
57 
58 // the rec hit matcher
59 #include "GSRecHitMatcher.h"
60 
61 // STL
62 #include <memory>
63 
64 // ROOT
65 #include <TFile.h>
66 #include <TH1F.h>
67 
68 
69 //#define FAMOS_DEBUG
70 
72  edm::ParameterSet const& conf)
73  : pset_(conf)
74 {
75  thePixelDataFile = 0;
82 
83  random = 0;
84 
85 
86 #ifdef FAMOS_DEBUG
87  std::cout << "SiTrackerGaussianSmearingRecHitConverter instantiated" << std::endl;
88 #endif
89 
90  // Initialize the random number generator service
92  if ( ! rng.isAvailable() ) {
93  throw cms::Exception("Configuration")
94  << "SiTrackerGaussianSmearingRecHitConverter requires the RandomGeneratorService\n"
95  "which is not present in the configuration file.\n"
96  "You must add the service in the configuration file\n"
97  "or remove the module that requires it";
98  }
99 
100  random = new RandomEngine(&(*rng));
101  //PAT
102  produces<FastTrackerClusterCollection>("TrackerClusters");
103 
104  produces<SiTrackerGSRecHit2DCollection>("TrackerGSRecHits");
105  produces<SiTrackerGSMatchedRecHit2DCollection>("TrackerGSMatchedRecHits");
106 
107  //--- PSimHit Containers
108  // trackerContainers.clear();
109  // trackerContainers = conf.getParameter<std::vector<edm::InputTag> >("ROUList");
110  inputSimHits = conf.getParameter<edm::InputTag>("InputSimHits");
111  //--- delta rays p cut [GeV/c] to filter PSimHits with p>
112  deltaRaysPCut = conf.getParameter<double>("DeltaRaysMomentumCut");
113 
114  //--- switch to have RecHit == PSimHit
115  trackingPSimHits = conf.getParameter<bool>("trackingPSimHits");
116  if(trackingPSimHits) std::cout << "### trackingPSimHits chosen " << trackingPSimHits << std::endl;
117 
118  // switch on/off matching
119  doMatching = conf.getParameter<bool>("doRecHitMatching");
120 
121  // disable/enable dead channels
122  doDisableChannels = conf.getParameter<bool>("killDeadChannels");
123 
124  // Switch between old (ORCA) and new (CMSSW) pixel parameterization
125  useCMSSWPixelParameterization = conf.getParameter<bool>("UseCMSSWPixelParametrization");
126 #ifdef FAMOS_DEBUG
127  std::cout << (useCMSSWPixelParameterization? "CMSSW" : "ORCA") << " pixel parametrization chosen in config file." << std::endl;
128 #endif
129 
130  //Clusters
131  GevPerElectron = conf.getParameter<double>("GevPerElectron");
132  ElectronsPerADC = conf.getParameter<double>("ElectronsPerADC");
133 
134  //
135  // TIB
136  localPositionResolution_TIB1x = conf.getParameter<double>("TIB1x");
137  localPositionResolution_TIB1y = conf.getParameter<double>("TIB1y");
138  localPositionResolution_TIB2x = conf.getParameter<double>("TIB2x");
139  localPositionResolution_TIB2y = conf.getParameter<double>("TIB2y");
140  localPositionResolution_TIB3x = conf.getParameter<double>("TIB3x");
141  localPositionResolution_TIB3y = conf.getParameter<double>("TIB3y");
142  localPositionResolution_TIB4x = conf.getParameter<double>("TIB4x");
143  localPositionResolution_TIB4y = conf.getParameter<double>("TIB4y");
144  //
145  // TID
146  localPositionResolution_TID1x = conf.getParameter<double>("TID1x");
147  localPositionResolution_TID1y = conf.getParameter<double>("TID1y");
148  localPositionResolution_TID2x = conf.getParameter<double>("TID2x");
149  localPositionResolution_TID2y = conf.getParameter<double>("TID2y");
150  localPositionResolution_TID3x = conf.getParameter<double>("TID3x");
151  localPositionResolution_TID3y = conf.getParameter<double>("TID3y");
152  //
153  // TOB
154  localPositionResolution_TOB1x = conf.getParameter<double>("TOB1x");
155  localPositionResolution_TOB1y = conf.getParameter<double>("TOB1y");
156  localPositionResolution_TOB2x = conf.getParameter<double>("TOB2x");
157  localPositionResolution_TOB2y = conf.getParameter<double>("TOB2y");
158  localPositionResolution_TOB3x = conf.getParameter<double>("TOB3x");
159  localPositionResolution_TOB3y = conf.getParameter<double>("TOB3y");
160  localPositionResolution_TOB4x = conf.getParameter<double>("TOB4x");
161  localPositionResolution_TOB4y = conf.getParameter<double>("TOB4y");
162  localPositionResolution_TOB5x = conf.getParameter<double>("TOB5x");
163  localPositionResolution_TOB5y = conf.getParameter<double>("TOB5y");
164  localPositionResolution_TOB6x = conf.getParameter<double>("TOB6x");
165  localPositionResolution_TOB6y = conf.getParameter<double>("TOB6y");
166  //
167  // TEC
168  localPositionResolution_TEC1x = conf.getParameter<double>("TEC1x");
169  localPositionResolution_TEC1y = conf.getParameter<double>("TEC1y");
170  localPositionResolution_TEC2x = conf.getParameter<double>("TEC2x");
171  localPositionResolution_TEC2y = conf.getParameter<double>("TEC2y");
172  localPositionResolution_TEC3x = conf.getParameter<double>("TEC3x");
173  localPositionResolution_TEC3y = conf.getParameter<double>("TEC3y");
174  localPositionResolution_TEC4x = conf.getParameter<double>("TEC4x");
175  localPositionResolution_TEC4y = conf.getParameter<double>("TEC4y");
176  localPositionResolution_TEC5x = conf.getParameter<double>("TEC5x");
177  localPositionResolution_TEC5y = conf.getParameter<double>("TEC5y");
178  localPositionResolution_TEC6x = conf.getParameter<double>("TEC6x");
179  localPositionResolution_TEC6y = conf.getParameter<double>("TEC6y");
180  localPositionResolution_TEC7x = conf.getParameter<double>("TEC7x");
181  localPositionResolution_TEC7y = conf.getParameter<double>("TEC7y");
182  //
183  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
184  //
185 #ifdef FAMOS_DEBUG
186  std::cout << "RecHit local position error set to" << "\n"
187  << "\tTIB1\tx = " << localPositionResolution_TIB1x
188  << " cm\ty = " << localPositionResolution_TIB1y << " cm" << "\n"
189  << "\tTIB2\tx = " << localPositionResolution_TIB2x
190  << " cm\ty = " << localPositionResolution_TIB2y << " cm" << "\n"
191  << "\tTIB3\tx = " << localPositionResolution_TIB3x
192  << " cm\ty = " << localPositionResolution_TIB3y << " cm" << "\n"
193  << "\tTIB4\tx = " << localPositionResolution_TIB4x
194  << " cm\ty = " << localPositionResolution_TIB4y << " cm" << "\n"
195  << "\tTID1\tx = " << localPositionResolution_TID1x
196  << " cm\ty = " << localPositionResolution_TID1y << " cm" << "\n"
197  << "\tTID2\tx = " << localPositionResolution_TID2x
198  << " cm\ty = " << localPositionResolution_TID2y << " cm" << "\n"
199  << "\tTID3\tx = " << localPositionResolution_TID3x
200  << " cm\ty = " << localPositionResolution_TID3y << " cm" << "\n"
201  << "\tTOB1\tx = " << localPositionResolution_TOB1x
202  << " cm\ty = " << localPositionResolution_TOB1y << " cm" << "\n"
203  << "\tTOB2\tx = " << localPositionResolution_TOB2x
204  << " cm\ty = " << localPositionResolution_TOB2y << " cm" << "\n"
205  << "\tTOB3\tx = " << localPositionResolution_TOB3x
206  << " cm\ty = " << localPositionResolution_TOB3y << " cm" << "\n"
207  << "\tTOB4\tx = " << localPositionResolution_TOB4x
208  << " cm\ty = " << localPositionResolution_TOB4y << " cm" << "\n"
209  << "\tTOB5\tx = " << localPositionResolution_TOB5x
210  << " cm\ty = " << localPositionResolution_TOB5y << " cm" << "\n"
211  << "\tTOB6\tx = " << localPositionResolution_TOB6x
212  << " cm\ty = " << localPositionResolution_TOB6y << " cm" << "\n"
213  << "\tTEC1\tx = " << localPositionResolution_TEC1x
214  << " cm\ty = " << localPositionResolution_TEC1y << " cm" << "\n"
215  << "\tTEC2\tx = " << localPositionResolution_TEC2x
216  << " cm\ty = " << localPositionResolution_TEC2y << " cm" << "\n"
217  << "\tTEC3\tx = " << localPositionResolution_TEC3x
218  << " cm\ty = " << localPositionResolution_TEC3y << " cm" << "\n"
219  << "\tTEC4\tx = " << localPositionResolution_TEC4x
220  << " cm\ty = " << localPositionResolution_TEC4y << " cm" << "\n"
221  << "\tTEC5\tx = " << localPositionResolution_TEC5x
222  << " cm\ty = " << localPositionResolution_TEC5y << " cm" << "\n"
223  << "\tTEC6\tx = " << localPositionResolution_TEC6x
224  << " cm\ty = " << localPositionResolution_TEC6y << " cm" << "\n"
225  << "\tTEC7\tx = " << localPositionResolution_TEC7x
226  << " cm\ty = " << localPositionResolution_TEC7y << " cm" << "\n"
227  << "\tAll:\tz = " << localPositionResolution_z << " cm"
228  << std::endl;
229 #endif
230 
231  //--- Number of histograms for alpha/beta barrel/forward multiplicity
232  if(useCMSSWPixelParameterization) {
233  nAlphaBarrel = conf.getParameter<int>("AlphaBarrelMultiplicityNew");
234  nBetaBarrel = conf.getParameter<int>("BetaBarrelMultiplicityNew");
235  nAlphaForward = conf.getParameter<int>("AlphaForwardMultiplicityNew");
236  nBetaForward = conf.getParameter<int>("BetaForwardMultiplicityNew");
237  } else {
238  nAlphaBarrel = conf.getParameter<int>("AlphaBarrelMultiplicity");
239  nBetaBarrel = conf.getParameter<int>("BetaBarrelMultiplicity");
240  nAlphaForward = conf.getParameter<int>("AlphaForwardMultiplicity");
241  nBetaForward = conf.getParameter<int>("BetaForwardMultiplicity");
242  }
243 #ifdef FAMOS_DEBUG
244  std::cout << "Pixel maximum multiplicity set to "
245  << "\nBarrel" << "\talpha " << nAlphaBarrel
246  << "\tbeta " << nBetaBarrel
247  << "\nForward" << "\talpha " << nAlphaForward
248  << "\tbeta " << nBetaForward
249  << std::endl;
250 #endif
251 
252  // Resolution Barrel
253  if(useCMSSWPixelParameterization) {
254  resAlphaBarrel_binMin = conf.getParameter<double>("AlphaBarrel_BinMinNew" );
255  resAlphaBarrel_binWidth = conf.getParameter<double>("AlphaBarrel_BinWidthNew");
256  resAlphaBarrel_binN = conf.getParameter<int>( "AlphaBarrel_BinNNew" );
257  resBetaBarrel_binMin = conf.getParameter<double>("BetaBarrel_BinMinNew" );
258  resBetaBarrel_binWidth = conf.getParameter<double>("BetaBarrel_BinWidthNew" );
259  resBetaBarrel_binN = conf.getParameter<int>( "BetaBarrel_BinNNew" );
260  } else {
261  resAlphaBarrel_binMin = conf.getParameter<double>("AlphaBarrel_BinMin" );
262  resAlphaBarrel_binWidth = conf.getParameter<double>("AlphaBarrel_BinWidth");
263  resAlphaBarrel_binN = conf.getParameter<int>( "AlphaBarrel_BinN" );
264  resBetaBarrel_binMin = conf.getParameter<double>("BetaBarrel_BinMin" );
265  resBetaBarrel_binWidth = conf.getParameter<double>("BetaBarrel_BinWidth" );
266  resBetaBarrel_binN = conf.getParameter<int>( "BetaBarrel_BinN" );
267  }
268 
269  // Resolution Forward
270  if(useCMSSWPixelParameterization) {
271  resAlphaForward_binMin = conf.getParameter<double>("AlphaForward_BinMinNew" );
272  resAlphaForward_binWidth = conf.getParameter<double>("AlphaForward_BinWidthNew" );
273  resAlphaForward_binN = conf.getParameter<int>( "AlphaForward_BinNNew" );
274  resBetaForward_binMin = conf.getParameter<double>("BetaForward_BinMinNew" );
275  resBetaForward_binWidth = conf.getParameter<double>("BetaForward_BinWidthNew" );
276  resBetaForward_binN = conf.getParameter<int>( "BetaForward_BinNNew" );
277  } else {
278  resAlphaForward_binMin = conf.getParameter<double>("AlphaForward_BinMin" );
279  resAlphaForward_binWidth = conf.getParameter<double>("AlphaForward_BinWidth" );
280  resAlphaForward_binN = conf.getParameter<int>( "AlphaForward_BinN" );
281  resBetaForward_binMin = conf.getParameter<double>("BetaForward_BinMin" );
282  resBetaForward_binWidth = conf.getParameter<double>("BetaForward_BinWidth" );
283  resBetaForward_binN = conf.getParameter<int>( "BetaForward_BinN" );
284  }
285 
286  // Hit Finding Probability
287  theHitFindingProbability_PXB = conf.getParameter<double>("HitFindingProbability_PXB" );
288  theHitFindingProbability_PXF = conf.getParameter<double>("HitFindingProbability_PXF" );
289  theHitFindingProbability_TIB1 = conf.getParameter<double>("HitFindingProbability_TIB1");
290  theHitFindingProbability_TIB2 = conf.getParameter<double>("HitFindingProbability_TIB2");
291  theHitFindingProbability_TIB3 = conf.getParameter<double>("HitFindingProbability_TIB3");
292  theHitFindingProbability_TIB4 = conf.getParameter<double>("HitFindingProbability_TIB4");
293  theHitFindingProbability_TID1 = conf.getParameter<double>("HitFindingProbability_TID1");
294  theHitFindingProbability_TID2 = conf.getParameter<double>("HitFindingProbability_TID2");
295  theHitFindingProbability_TID3 = conf.getParameter<double>("HitFindingProbability_TID3");
296  theHitFindingProbability_TOB1 = conf.getParameter<double>("HitFindingProbability_TOB1");
297  theHitFindingProbability_TOB2 = conf.getParameter<double>("HitFindingProbability_TOB2");
298  theHitFindingProbability_TOB3 = conf.getParameter<double>("HitFindingProbability_TOB3");
299  theHitFindingProbability_TOB4 = conf.getParameter<double>("HitFindingProbability_TOB4");
300  theHitFindingProbability_TOB5 = conf.getParameter<double>("HitFindingProbability_TOB5");
301  theHitFindingProbability_TOB6 = conf.getParameter<double>("HitFindingProbability_TOB6");
302  theHitFindingProbability_TEC1 = conf.getParameter<double>("HitFindingProbability_TEC1");
303  theHitFindingProbability_TEC2 = conf.getParameter<double>("HitFindingProbability_TEC2");
304  theHitFindingProbability_TEC3 = conf.getParameter<double>("HitFindingProbability_TEC3");
305  theHitFindingProbability_TEC4 = conf.getParameter<double>("HitFindingProbability_TEC4");
306  theHitFindingProbability_TEC5 = conf.getParameter<double>("HitFindingProbability_TEC5");
307  theHitFindingProbability_TEC6 = conf.getParameter<double>("HitFindingProbability_TEC6");
308  theHitFindingProbability_TEC7 = conf.getParameter<double>("HitFindingProbability_TEC7");
309  //
310 #ifdef FAMOS_DEBUG
311  std::cout << "RecHit finding probability set to" << "\n"
312  << "\tPXB = " << theHitFindingProbability_PXB << "\n"
313  << "\tPXF = " << theHitFindingProbability_PXF << "\n"
314  << "\tTIB1 = " << theHitFindingProbability_TIB1 << "\n"
315  << "\tTIB2 = " << theHitFindingProbability_TIB2 << "\n"
316  << "\tTIB3 = " << theHitFindingProbability_TIB3 << "\n"
317  << "\tTIB4 = " << theHitFindingProbability_TIB4 << "\n"
318  << "\tTID1 = " << theHitFindingProbability_TID1 << "\n"
319  << "\tTID2 = " << theHitFindingProbability_TID2 << "\n"
320  << "\tTID3 = " << theHitFindingProbability_TID3 << "\n"
321  << "\tTOB1 = " << theHitFindingProbability_TOB1 << "\n"
322  << "\tTOB2 = " << theHitFindingProbability_TOB2 << "\n"
323  << "\tTOB3 = " << theHitFindingProbability_TOB3 << "\n"
324  << "\tTOB4 = " << theHitFindingProbability_TOB4 << "\n"
325  << "\tTOB5 = " << theHitFindingProbability_TOB5 << "\n"
326  << "\tTOB6 = " << theHitFindingProbability_TOB6 << "\n"
327  << "\tTEC1 = " << theHitFindingProbability_TEC1 << "\n"
328  << "\tTEC2 = " << theHitFindingProbability_TEC2 << "\n"
329  << "\tTEC3 = " << theHitFindingProbability_TEC3 << "\n"
330  << "\tTEC4 = " << theHitFindingProbability_TEC4 << "\n"
331  << "\tTEC5 = " << theHitFindingProbability_TEC5 << "\n"
332  << "\tTEC6 = " << theHitFindingProbability_TEC6 << "\n"
333  << "\tTEC7 = " << theHitFindingProbability_TEC7 << "\n"
334  << std::endl;
335 #endif
336 
337  // Initialize the si strip error parametrization
340 
341  // Initialization of pixel parameterization posponed to beginRun(), since it depends on the magnetic field
342 
343 }
344 
345 
347  // load multiplicity cumulative probabilities
348  // root files
349  thePixelDataFile = new TFile ( edm::FileInPath( thePixelMultiplicityFileName ).fullPath().c_str() , "READ" );
350  thePixelBarrelResolutionFile = new TFile ( edm::FileInPath( thePixelBarrelResolutionFileName ).fullPath().c_str() , "READ" );
351  thePixelForwardResolutionFile = new TFile ( edm::FileInPath( thePixelForwardResolutionFileName ).fullPath().c_str() , "READ" );
352  //
353 
354  // alpha barrel
356  nAlphaBarrel ,
357  std::string("hist_alpha_barrel") ,
359  //
360  // beta barrel
362  nBetaBarrel ,
363  std::string("hist_beta_barrel") ,
365  //
366  // alpha forward
368  nAlphaForward ,
369  std::string("hist_alpha_forward") ,
371  //
372  // beta forward
374  nBetaForward ,
375  std::string("hist_beta_forward") ,
377 
378  // Load also big pixel data if CMSSW parametrization is on
379  // They are pushed back into the vectors after the normal pixels data:
380  // [0, ..., (size/2)-1] -> Normal pixels
381  // [size/2, ..., size-1] -> Big pixels
383  // alpha barrel
385  nAlphaBarrel ,
386  std::string("hist_alpha_barrel_big") ,
388  true );
389  //
390  // beta barrel
392  nBetaBarrel ,
393  std::string("hist_beta_barrel_big") ,
395  true );
396  //
397  // alpha forward
399  nAlphaForward ,
400  std::string("hist_alpha_forward_big") ,
402  true );
403  //
404  // beta forward
406  nBetaForward ,
407  std::string("hist_beta_forward_big") ,
409  true );
410  }
411  //
412 }
413 
415  TFile* pixelDataFile,
416  unsigned int nMultiplicity,
417  std::string histName,
418  std::vector<TH1F*>& theMultiplicityCumulativeProbabilities,
419  bool bigPixels)
420 {
421 
422  std::string histName_i = histName + "_%u"; // needed to open histograms with a for
423  if(!bigPixels)
424  theMultiplicityCumulativeProbabilities.clear();
425  //
426  // What's this vector? Not needed - MG
427 // std::vector<double> mult; // vector with fixed multiplicity
428  for(unsigned int i = 0; i<nMultiplicity; ++i) {
429  TH1F addHist = *((TH1F*) pixelDataFile->Get( Form( histName_i.c_str() ,i+1 )));
430  if(i==0) {
431  theMultiplicityCumulativeProbabilities.push_back( new TH1F(addHist) );
432  } else {
433  TH1F sumHist;
434  if(bigPixels)
435  sumHist = *(theMultiplicityCumulativeProbabilities[nMultiplicity+i-1]);
436  else
437  sumHist = *(theMultiplicityCumulativeProbabilities[i-1]);
438  sumHist.Add(&addHist);
439  theMultiplicityCumulativeProbabilities.push_back( new TH1F(sumHist) );
440  }
441  }
442 
443  // Logger
444 #ifdef FAMOS_DEBUG
445  const unsigned int maxMult = theMultiplicityCumulativeProbabilities.size();
446  unsigned int iMult, multSize;
448  if(bigPixels) {
449  iMult = maxMult / 2;
450  multSize = maxMult ;
451  } else {
452  iMult = 0;
453  multSize = maxMult;
454  }
455  } else {
456  iMult = 0;
457  multSize = maxMult ;
458  }
459  std::cout << " Multiplicity cumulated probability " << histName << std::endl;
460  for(/* void */; iMult<multSize; ++iMult) {
461  for(int iBin = 1; iBin<=theMultiplicityCumulativeProbabilities[iMult]->GetNbinsX(); ++iBin) {
462  std::cout
463  << " Multiplicity " << iMult+1
464  << " bin " << iBin
465  << " low edge = "
466  << theMultiplicityCumulativeProbabilities[iMult]->GetBinLowEdge(iBin)
467  << " prob = "
468  << (theMultiplicityCumulativeProbabilities[iMult])->GetBinContent(iBin)
469  // remember in ROOT bin starts from 1 (0 underflow, nBin+1 overflow)
470  << std::endl;
471  }
472  }
473 #endif
474 
475 }
476 
477 // Destructor
483 
490 
492 
493  if(random) delete random;
494 
495 }
496 
497 void
499 {
500 
501  // Initialize the Tracker Geometry
502  edm::ESHandle<TrackerGeometry> theGeometry;
503  es.get<TrackerDigiGeometryRecord> ().get (theGeometry);
504  geometry = &(*theGeometry);
505 
506  edm::ESHandle<TrackerGeometry> theMisAlignedGeometry;
507  es.get<TrackerDigiGeometryRecord>().get("MisAligned",theMisAlignedGeometry);
508  misAlignedGeometry = &(*theMisAlignedGeometry);
509 
510  const MagneticField* magfield;
512  es.get<IdealMagneticFieldRecord>().get(magField);
513  magfield=&(*magField);
514  GlobalPoint center(0.0, 0.0, 0.0);
515  double magFieldAtCenter = magfield->inTesla(center).mag();
516 
517  // For new parameterization: select multiplicity and resolution files according to magnetic field
519  if(magFieldAtCenter > 3.9) {
520  thePixelMultiplicityFileName = pset_.getParameter<std::string>( "PixelMultiplicityFile40T");
521  thePixelBarrelResolutionFileName = pset_.getParameter<std::string>( "PixelBarrelResolutionFile40T");
522  thePixelForwardResolutionFileName = pset_.getParameter<std::string>( "PixelForwardResolutionFile40T");
523  } else {
524  thePixelMultiplicityFileName = pset_.getParameter<std::string>( "PixelMultiplicityFile38T");
525  thePixelBarrelResolutionFileName = pset_.getParameter<std::string>( "PixelBarrelResolutionFile38T");
526  thePixelForwardResolutionFileName = pset_.getParameter<std::string>( "PixelForwardResolutionFile38T");
527  }
528  } else {
529  thePixelMultiplicityFileName = pset_.getParameter<std::string>( "PixelMultiplicityFile" );
530  thePixelBarrelResolutionFileName = pset_.getParameter<std::string>( "PixelBarrelResolutionFile");
531  thePixelForwardResolutionFileName = pset_.getParameter<std::string>( "PixelForwardResolutionFile");
532  }
533 
534 
535  // Reading the list of dead pixel modules from DB:
536  edm::ESHandle<SiPixelQuality> siPixelBadModule;
537  es.get<SiPixelQualityRcd>().get(siPixelBadModule);
539  if (doDisableChannels) {
540  disabledModules = new std::vector<SiPixelQuality::disabledModuleType> ( siPixelBadModule->getBadComponentList() );
542  size_t numberOfRecoverableModules = 0;
543  for (size_t id=0;id<numberOfDisabledModules;id++) {
545  // errortype "whole" = int 0 in DB //
546  // errortype "tbmA" = int 1 in DB //
547  // errortype "tbmB" = int 2 in DB //
548  // errortype "none" = int 3 in DB //
550  if ( (*disabledModules)[id-numberOfRecoverableModules].errorType != 0 ){
551  // Disable only the modules totally in error:
552  disabledModules->erase(disabledModules->begin()+id-numberOfRecoverableModules);
553  numberOfRecoverableModules++;
554  }
555  }
556  numberOfDisabledModules = disabledModules->size();
557  }
558 
559 
560 
561 #ifdef FAMOS_DEBUG
562  std::cout << "Pixel multiplicity data are taken from file " << thePixelMultiplicityFileName << std::endl;
563 
564  std::cout << "Pixel maximum multiplicity set to "
565  << "\nBarrel" << "\talpha " << nAlphaBarrel
566  << "\tbeta " << nBetaBarrel
567  << "\nForward" << "\talpha " << nAlphaForward
568  << "\tbeta " << nBetaForward
569  << std::endl;
570 
571  std::cout << "Barrel Pixel resolution data are taken from file "
573  << "Alpha bin min = " << resAlphaBarrel_binMin
574  << "\twidth = " << resAlphaBarrel_binWidth
575  << "\tbins = " << resAlphaBarrel_binN
576  << "\n"
577  << " Beta bin min = " << resBetaBarrel_binMin
578  << "\twidth = " << resBetaBarrel_binWidth
579  << "\tbins = " << resBetaBarrel_binN
580  << std::endl;
581 
582  std::cout << "Forward Pixel resolution data are taken from file "
584  << "Alpha bin min = " << resAlphaForward_binMin
585  << "\twidth = " << resAlphaForward_binWidth
586  << "\tbins = " << resAlphaForward_binN
587  << "\n"
588  << " Beta bin min = " << resBetaForward_binMin
589  << "\twidth = " << resBetaForward_binWidth
590  << "\tbins = " << resBetaForward_binN
591  << std::endl;
592 #endif
593  //
594 
595  //
596  // load pixel data
597  loadPixelData();
598  //
599 
600  // Initialize and open relevant files for the pixel barrel error parametrization
603  pset_,
605  random);
606  // Initialize and open relevant files for the pixel forward error parametrization
609  pset_,
611  random);
612 }
613 
615 {
616 
617  //Retrieve tracker topology from geometry
619  es.get<IdealGeometryRecord>().get(tTopoHand);
620  const TrackerTopology *tTopo=tTopoHand.product();
621 
622 
623  // Step 0: Declare Ref and RefProd
625 
626  // Step A: Get Inputs (PSimHit's)
627  /*
628  edm::Handle<CrossingFrame<PSimHit> > cf_simhit;
629  std::vector<const CrossingFrame<PSimHit> *> cf_simhitvec;
630  for(uint32_t i=0; i<trackerContainers.size(); i++){
631  e.getByLabel(trackerContainers[i], cf_simhit);
632  cf_simhitvec.push_back(cf_simhit.product());
633  }
634  std::auto_ptr<MixCollection<PSimHit> > allTrackerHits(new MixCollection<PSimHit>(cf_simhitvec));
635  */
636 
637  edm::Handle<edm::PSimHitContainer> allTrackerHits_handle;
638  e.getByLabel(inputSimHits,allTrackerHits_handle);
639  const edm::PSimHitContainer& allTrackerHits=*allTrackerHits_handle;
640 
641  // Step B: create temporary RecHit collection and fill it with Gaussian smeared RecHit's
642 
643  std::map<unsigned, edm::OwnVector<SiTrackerGSRecHit2D> > temporaryRecHits;
644  std::map<unsigned, edm::OwnVector<FastTrackerCluster> > theClusters ;
645 
646  smearHits( allTrackerHits, temporaryRecHits, theClusters, tTopo);
647 
648  // Step C: match rechits on stereo layers
649  std::map<unsigned, edm::OwnVector<SiTrackerGSMatchedRecHit2D> > temporaryMatchedRecHits ;
650  if(doMatching) matchHits( temporaryRecHits, temporaryMatchedRecHits);//, allTrackerHits);
651 
652  // Step D: from the temporary RecHit collection, create the real one.
653  std::auto_ptr<SiTrackerGSRecHit2DCollection> recHitCollection(new SiTrackerGSRecHit2DCollection);
654  std::auto_ptr<SiTrackerGSMatchedRecHit2DCollection> recHitCollectionMatched(new SiTrackerGSMatchedRecHit2DCollection);
655  if(doMatching){
656  loadMatchedRecHits(temporaryMatchedRecHits, *recHitCollectionMatched);
657  loadRecHits(temporaryRecHits, *recHitCollection);
658  }
659  else {
660  //might need to have a "matched" hit collection containing the simple hits
661  loadRecHits(temporaryRecHits, *recHitCollection);
662  }
663 
664 
665 
666  //std::cout << "****** TrackerGSRecHits hits are =\t" << (*recHitCollection).size()<<std::endl;
667  //std::cout << "****** TrackerGSRecHitsMatched hits are =\t" << (*recHitCollectionMatched).size()<< std::endl;
668 
669 
670 
671  // Step E: write output to file
672  e.put(recHitCollection,"TrackerGSRecHits");
673  e.put(recHitCollectionMatched,"TrackerGSMatchedRecHits");
674 
675  //STEP F: write clusters
676  std::auto_ptr<FastTrackerClusterCollection> clusterCollection(new FastTrackerClusterCollection);
677  loadClusters(theClusters, *clusterCollection);
678  e.put(clusterCollection,"TrackerClusters");
679 }
680 
681 
682 
683 //void SiTrackerGaussianSmearingRecHitConverter::smearHits(MixCollection<PSimHit>& input,
685  std::map<unsigned, edm::OwnVector<SiTrackerGSRecHit2D> >& temporaryRecHits,
686  std::map<unsigned, edm::OwnVector<FastTrackerCluster> >& theClusters,
687  const TrackerTopology *tTopo)
688 {
689 
690  int numberOfPSimHits = 0;
691 
692 
693  // MixCollection<PSimHit>::iterator isim = input.begin();
694  // MixCollection<PSimHit>::iterator lastSimHit = input.end();
695  edm::PSimHitContainer::const_iterator isim = input.begin();
696  edm::PSimHitContainer::const_iterator lastSimHit = input.end();
697 
698 
699 
700  correspondingSimHit.resize(input.size());
703  LocalError inflatedError;
704 
705  int simHitCounter = -1;
706  int recHitCounter = 0;
707 
708 
709  // loop on PSimHits
710 
711  for ( ; isim != lastSimHit; ++isim ) {
712  ++simHitCounter;
713 
714  DetId det((*isim).detUnitId());
715  unsigned trackID = (*isim).trackId();
716  uint32_t eeID = (*isim).eventId().rawId(); //get the rawId of the eeid for pileup treatment
717 
718 
719 
720 
722 
723  // unsigned int subdetId = det.subdetId();
724  // unsigned int disk = tTopo->pxfDisk(det);
725  // unsigned int side = tTopo->pxfSide(det);
726  // std::cout<< " Pixel Forward Disk Number : "<< disk << " Side : "<<side<<std::endl;
727  // if(subdetId==1 || subdetId==2) std::cout<<" Pixel GSRecHits "<<std::endl;
728  // else if(subdetId==3|| subdetId==4 || subdetId==5 || subdetId == 6) std::cout<<" Strip GSRecHits "<<std::endl;
729 
730  bool isBad = false;
731  unsigned int geoId = det.rawId();
732  for (size_t id=0;id<numberOfDisabledModules;id++) {
733  if(geoId==(*disabledModules)[id].DetID){
734  // Already selected in the beginRun() the ones with errorType = 0
735  // if((*disabledModules)[id].errorType == 0) isBad = true;
736  isBad = true;
737  break;
738  }
739  }
740  if(isBad) continue;
741 
742 
743  ++numberOfPSimHits;
744  // gaussian smearing
745  unsigned int alphaMult = 0;
746  unsigned int betaMult = 0;
747  bool isCreated = gaussianSmearing(*isim, position, error, alphaMult, betaMult, tTopo);
748  // std::cout << "Error as simulated " << error.xx() << " " << error.xy() << " " << error.yy() << std::endl;
749  //
750  unsigned int subdet = det.subdetId();
751 
752  //
753  if(isCreated) {
754  // double dist = theDet->surface().toGlobal((*isim).localPosition()).mag2();
755  // create RecHit
756  // Fill the temporary RecHit on the current DetId collection
757  // temporaryRecHits[trackID].push_back(
758 
759  // Fix by P.Azzi, A.Schmidt:
760  // if this hit is on a glued detector with radial topology,
761  // the x-coordinate needs to be translated correctly to y=0.
762  // this is necessary as intermediate step for the matching of hits
763 
764  // do it only if we do rec hit matching
765 // if(doMatching){
766 // StripSubdetector specDetId=StripSubdetector(det);
767 
768 // if(specDetId.glued()) if(subdet==6 || subdet==4){ // check if on double layer in TEC or TID
769 
770 // const GeomDetUnit* theDetUnit = geometry->idToDetUnit(det);
771 // const StripTopology& topol=(const StripTopology&)theDetUnit->type().topology();
772 
773 // // check if radial topology
774 // if(dynamic_cast <const RadialStripTopology*>(&topol)){
775 
776 // const RadialStripTopology *rtopol = dynamic_cast <const RadialStripTopology*>(&topol);
777 
778 // // translate to measurement coordinates
779 // MeasurementPoint mpoint=rtopol->measurementPosition(position);
780 // // set y=0
781 // MeasurementPoint mpoint0=MeasurementPoint(mpoint.x(),0.);
782 // // translate back to local coordinates with y=0
783 // LocalPoint lpoint0 = rtopol->localPosition(mpoint0);
784 // position = Local3DPoint(lpoint0.x(),0.,0.);
785 
786 // }
787 // } // end if( specDetId.glued()...
788 // } // end if(doMatching)
789 
790  // set y=0 in single-sided strip detectors
791  if(doMatching && (subdet>2)){
792  StripSubdetector specDetId=StripSubdetector(det);
793  if( !specDetId.glued() ){ // this is a single-sided module
794  position = Local3DPoint(position.x(),0.,0.); // set y to 0 on single-sided modules
795  }
796  }
797  else{ if(subdet>2) position = Local3DPoint(position.x(),0.,0.); } // no matching, set y=0 on strips
798 
799  // Inflate errors in case of geometry misalignment
800  const GeomDet* theMADet = misAlignedGeometry->idToDet(det);
801  if ( theMADet->alignmentPositionError() != 0 ) {
802  LocalError lape =
804  theMADet->surface() );
805  error = LocalError ( error.xx()+lape.xx(),
806  error.xy()+lape.xy(),
807  error.yy()+lape.yy() );
808  }
809 
810  float chargeADC = (*isim).energyLoss()/(GevPerElectron * ElectronsPerADC);
811 
812  //create cluster
813  if(subdet > 2) theClusters[trackID].push_back(
814  new FastTrackerCluster(LocalPoint(position.x(), 0.0, 0.0), error, det,
815  simHitCounter, trackID,
816  eeID,
817  //(*isim).energyLoss())
818  chargeADC)
819  );
820  else theClusters[trackID].push_back(
821  new FastTrackerCluster(position, error, det,
822  simHitCounter, trackID,
823  eeID,
824  //(*isim).energyLoss())
825  chargeADC)
826  );
827 
828  // std::cout << "CLUSTER for simhit " << simHitCounter << "\t energy loss = " <<chargeADC << std::endl;
829 
830  // std::cout << "Error as reconstructed " << error.xx() << " " << error.xy() << " " << error.yy() << std::endl;
831 
832  // create rechit
833  temporaryRecHits[trackID].push_back(
834  new SiTrackerGSRecHit2D(position, error, det,
835  simHitCounter, trackID,
836  eeID,
837  ClusterRef(FastTrackerClusterRefProd, simHitCounter),
838  alphaMult, betaMult)
839  );
840 
841  // This a correpondence map between RecHits and SimHits
842  // (for later use in matchHits)
843  correspondingSimHit[recHitCounter++] = isim;
844 
845 
846  } // end if(isCreated)
847 
848  } // end loop on PSimHits
849 
850 
851 }
852 
855  LocalError& error,
856  unsigned& alphaMult,
857  unsigned& betaMult,
858  const TrackerTopology *tTopo)
859 {
860 
861  // A few caracteritics of the detid the SimHit belongs to.
862  unsigned int subdet = DetId(simHit.detUnitId()).subdetId();
863  unsigned int detid = DetId(simHit.detUnitId()).rawId();
864  const GeomDetUnit* theDetUnit = geometry->idToDetUnit((DetId)simHit.detUnitId());
865  const BoundPlane& theDetPlane = theDetUnit->surface();
866  const Bounds& theBounds = theDetPlane.bounds();
867  double boundX = theBounds.width()/2.;
868  double boundY = theBounds.length()/2.;
869 
870 #ifdef FAMOS_DEBUG
871  std::cout << "\tSubdetector " << subdet
872  << " rawid " << detid
873  << std::endl;
874 #endif
875  if(trackingPSimHits) {
876  // z is fixed for all detectors, in case of errors resolution is fixed also for x and y to 1 um (zero)
877  // The Matrix is the Covariance Matrix, sigma^2 on diagonal!!!
879  0.0 ,
880  localPositionResolution_z * localPositionResolution_z );
881  //
882  // starting from PSimHit local position
883  position = simHit.localPosition();
884 #ifdef FAMOS_DEBUG
885  std::cout << " Tracking PSimHit position set to " << position;
886 #endif
887  return true; // RecHit == PSimHit with 100% hit finding efficiency
888  }
889  //
890 
891  // hit finding probability --> RecHit will be created if and only if hitFindingProbability <= theHitFindingProbability_###
892  double hitFindingProbability = random->flatShoot();
893 #ifdef FAMOS_DEBUG
894  std::cout << " Hit finding probability draw: " << hitFindingProbability << std::endl;;
895 #endif
896 
897  switch (subdet) {
898  // Pixel Barrel
899  case 1:
900  {
901 #ifdef FAMOS_DEBUG
902 
903  unsigned int theLayer = tTopo->pxbLayer(detid);
904  std::cout << "\tPixel Barrel Layer " << theLayer << std::endl;
905 #endif
906  if( hitFindingProbability > theHitFindingProbability_PXB ) return false;
907  // Hit smearing
908  const PixelGeomDetUnit* pixelDetUnit = dynamic_cast<const PixelGeomDetUnit*>(theDetUnit);
909  thePixelBarrelParametrization->smearHit(simHit, pixelDetUnit, boundX, boundY);
914  return true;
915  break;
916  }
917  // Pixel Forward
918  case 2:
919  {
920 #ifdef FAMOS_DEBUG
921 
922  unsigned int theDisk = tTopo->pxfDisk(detid);
923  std::cout << "\tPixel Forward Disk " << theDisk << std::endl;
924 #endif
925  if( hitFindingProbability > theHitFindingProbability_PXF ) return false;
926  // Hit smearing
927  const PixelGeomDetUnit* pixelDetUnit = dynamic_cast<const PixelGeomDetUnit*>(theDetUnit);
928  thePixelEndcapParametrization->smearHit(simHit, pixelDetUnit, boundX, boundY);
933  return true;
934  break;
935  }
936  // TIB
937  case 3:
938  {
939 
940  unsigned int theLayer = tTopo->tibLayer(detid);
941 #ifdef FAMOS_DEBUG
942  std::cout << "\tTIB Layer " << theLayer << std::endl;
943 #endif
944  //
945  double resolutionX, resolutionY, resolutionZ;
946  resolutionZ = localPositionResolution_z;
947 
948  switch (theLayer) {
949  case 1:
950  {
951  resolutionX = localPositionResolution_TIB1x;
952  resolutionY = localPositionResolution_TIB1y;
953  if( hitFindingProbability > theHitFindingProbability_TIB1 ) return false;
954  break;
955  }
956  case 2:
957  {
958  resolutionX = localPositionResolution_TIB2x;
959  resolutionY = localPositionResolution_TIB2y;
960  if( hitFindingProbability > theHitFindingProbability_TIB2 ) return false;
961  break;
962  }
963  case 3:
964  {
965  resolutionX = localPositionResolution_TIB3x;
966  resolutionY = localPositionResolution_TIB3y;
967  if( hitFindingProbability > theHitFindingProbability_TIB3 ) return false;
968  break;
969  }
970  case 4:
971  {
972  resolutionX = localPositionResolution_TIB4x;
973  resolutionY = localPositionResolution_TIB4y;
974  if( hitFindingProbability > theHitFindingProbability_TIB4 ) return false;
975  break;
976  }
977  default:
978  {
979  edm::LogError ("SiTrackerGaussianSmearingRecHits")
980  << "\tTIB Layer not valid " << theLayer << std::endl;
981  return false;
982  break;
983  }
984  }
985 
986  // Gaussian smearing
987  theSiStripErrorParametrization->smearHit(simHit, resolutionX, resolutionY, resolutionZ, boundX, boundY);
990  alphaMult = 0;
991  betaMult = 0;
992  return true;
993  break;
994  } // TIB
995 
996  // TID
997  case 4:
998  {
999 
1000  unsigned int theRing = tTopo->tidRing(detid);
1001  double resolutionFactorY =
1002  1. - simHit.localPosition().y() / theDetPlane.position().perp();
1003 
1004 #ifdef FAMOS_DEBUG
1005  std::cout << "\tTID Ring " << theRing << std::endl;
1006 #endif
1007  double resolutionX, resolutionY, resolutionZ;
1008  resolutionZ = localPositionResolution_z;
1009 
1010  switch (theRing) {
1011  case 1:
1012  {
1013  resolutionX = localPositionResolution_TID1x * resolutionFactorY;
1014  resolutionY = localPositionResolution_TID1y;
1015  if( hitFindingProbability > theHitFindingProbability_TID1 ) return false;
1016  break;
1017  }
1018  case 2:
1019  {
1020  resolutionX = localPositionResolution_TID2x * resolutionFactorY;
1021  resolutionY = localPositionResolution_TID2y;
1022  if( hitFindingProbability > theHitFindingProbability_TID2 ) return false;
1023  break;
1024  }
1025  case 3:
1026  {
1027  resolutionX = localPositionResolution_TID3x * resolutionFactorY;
1028  resolutionY = localPositionResolution_TID3y;
1029  if( hitFindingProbability > theHitFindingProbability_TID3 ) return false;
1030  break;
1031  }
1032  default:
1033  {
1034  edm::LogError ("SiTrackerGaussianSmearingRecHits")
1035  << "\tTID Ring not valid " << theRing << std::endl;
1036  return false;
1037  break;
1038  }
1039  }
1040 
1041  boundX *= resolutionFactorY;
1042 
1043  theSiStripErrorParametrization->smearHit(simHit, resolutionX, resolutionY, resolutionZ, boundX, boundY);
1046  alphaMult = 0;
1047  betaMult = 0;
1048  return true;
1049  break;
1050  } // TID
1051 
1052  // TOB
1053  case 5:
1054  {
1055 
1056  unsigned int theLayer = tTopo->tobLayer(detid);
1057 #ifdef FAMOS_DEBUG
1058  std::cout << "\tTOB Layer " << theLayer << std::endl;
1059 #endif
1060  double resolutionX, resolutionY, resolutionZ;
1061  resolutionZ = localPositionResolution_z;
1062 
1063  switch (theLayer) {
1064  case 1:
1065  {
1066  resolutionX = localPositionResolution_TOB1x;
1067  resolutionY = localPositionResolution_TOB1y;
1068  if( hitFindingProbability > theHitFindingProbability_TOB1 ) return false;
1069  break;
1070  }
1071  case 2:
1072  {
1073  resolutionX = localPositionResolution_TOB2x;
1074  resolutionY = localPositionResolution_TOB2y;
1075  if( hitFindingProbability > theHitFindingProbability_TOB2 ) return false;
1076  break;
1077  }
1078  case 3:
1079  {
1080  resolutionX = localPositionResolution_TOB3x;
1081  resolutionY = localPositionResolution_TOB3y;
1082  if( hitFindingProbability > theHitFindingProbability_TOB3 ) return false;
1083  break;
1084  }
1085  case 4:
1086  {
1087  resolutionX = localPositionResolution_TOB4x;
1088  resolutionY = localPositionResolution_TOB4y;
1089  if( hitFindingProbability > theHitFindingProbability_TOB4 ) return false;
1090  break;
1091  }
1092  case 5:
1093  {
1094  resolutionX = localPositionResolution_TOB5x;
1095  resolutionY = localPositionResolution_TOB5y;
1096  if( hitFindingProbability > theHitFindingProbability_TOB5 ) return false;
1097  break;
1098  }
1099  case 6:
1100  {
1101  resolutionX = localPositionResolution_TOB6x;
1102  resolutionY = localPositionResolution_TOB6y;
1103  if( hitFindingProbability > theHitFindingProbability_TOB6 ) return false;
1104  break;
1105  }
1106  default:
1107  {
1108  edm::LogError ("SiTrackerGaussianSmearingRecHits")
1109  << "\tTOB Layer not valid " << theLayer << std::endl;
1110  return false;
1111  break;
1112  }
1113  }
1114  theSiStripErrorParametrization->smearHit(simHit, resolutionX, resolutionY, resolutionZ, boundX, boundY);
1117  alphaMult = 0;
1118  betaMult = 0;
1119  return true;
1120  break;
1121  } // TOB
1122 
1123  // TEC
1124  case 6:
1125  {
1126 
1127  unsigned int theRing = tTopo->tecRing(detid);
1128  double resolutionFactorY =
1129  1. - simHit.localPosition().y() / theDetPlane.position().perp();
1130 
1131 #ifdef FAMOS_DEBUG
1132  std::cout << "\tTEC Ring " << theRing << std::endl;
1133 #endif
1134  double resolutionX, resolutionY, resolutionZ;
1136 
1137  switch (theRing) {
1138  case 1:
1139  {
1140  resolutionX = localPositionResolution_TEC1x * resolutionFactorY;
1141  resolutionY = localPositionResolution_TEC1y;
1142  if( hitFindingProbability > theHitFindingProbability_TEC1 ) return false;
1143  break;
1144  }
1145  case 2:
1146  {
1147  resolutionX = localPositionResolution_TEC2x * resolutionFactorY;
1148  resolutionY = localPositionResolution_TEC2y;
1149  if( hitFindingProbability > theHitFindingProbability_TEC2 ) return false;
1150  break;
1151  }
1152  case 3:
1153  {
1154  resolutionX = localPositionResolution_TEC3x * resolutionFactorY;
1155  resolutionY = localPositionResolution_TEC3y;
1156  if( hitFindingProbability > theHitFindingProbability_TEC3 ) return false;
1157  break;
1158  }
1159  case 4:
1160  {
1161  resolutionX = localPositionResolution_TEC4x * resolutionFactorY;
1162  resolutionY = localPositionResolution_TEC4y;
1163  if( hitFindingProbability > theHitFindingProbability_TEC4 ) return false;
1164  break;
1165  }
1166  case 5:
1167  {
1168  resolutionX = localPositionResolution_TEC5x * resolutionFactorY;
1169  resolutionY = localPositionResolution_TEC5y;
1170  if( hitFindingProbability > theHitFindingProbability_TEC5 ) return false;
1171  break;
1172  }
1173  case 6:
1174  {
1175  resolutionX = localPositionResolution_TEC6x * resolutionFactorY;
1176  resolutionY = localPositionResolution_TEC6y;
1177  if( hitFindingProbability > theHitFindingProbability_TEC6 ) return false;
1178  break;
1179  }
1180  case 7:
1181  {
1182  resolutionX = localPositionResolution_TEC7x * resolutionFactorY;
1183  resolutionY = localPositionResolution_TEC7y;
1184  if( hitFindingProbability > theHitFindingProbability_TEC7 ) return false;
1185  break;
1186  }
1187  default:
1188  {
1189  edm::LogError ("SiTrackerGaussianSmearingRecHits")
1190  << "\tTEC Ring not valid " << theRing << std::endl;
1191  return false;
1192  break;
1193  }
1194  }
1195 
1196  boundX *= resolutionFactorY;
1197  theSiStripErrorParametrization->smearHit(simHit, resolutionX, resolutionY, resolutionZ, boundX, boundY);
1200  alphaMult = 0;
1201  betaMult = 0;
1202  return true;
1203  break;
1204  } // TEC
1205 
1206  default:
1207  {
1208  edm::LogError ("SiTrackerGaussianSmearingRecHits") << "\tTracker subdetector not valid " << subdet << std::endl;
1209  return false;
1210  break;
1211  }
1212 
1213  } // subdetector case
1214  //
1215 }
1216 
1217 
1218 
1219 void
1221  std::map<unsigned,edm::OwnVector<FastTrackerCluster> >& theClusterMap,
1222  FastTrackerClusterCollection& theClusterCollection) const
1223 {
1224  std::map<unsigned,edm::OwnVector<FastTrackerCluster> >::const_iterator
1225  it = theClusterMap.begin();
1226  std::map<unsigned,edm::OwnVector<FastTrackerCluster> >::const_iterator
1227  lastCluster = theClusterMap.end();
1228 
1229  for( ; it != lastCluster ; ++it ) {
1230  theClusterCollection.put(it->first,(it->second).begin(),(it->second).end());
1231  }
1232 }
1233 
1234 
1235 
1236 void
1238  std::map<unsigned,edm::OwnVector<SiTrackerGSRecHit2D> >& theRecHits,
1239  SiTrackerGSRecHit2DCollection& theRecHitCollection) const
1240 {
1241  std::map<unsigned,edm::OwnVector<SiTrackerGSRecHit2D> >::const_iterator
1242  it = theRecHits.begin();
1243  std::map<unsigned,edm::OwnVector<SiTrackerGSRecHit2D> >::const_iterator
1244  lastRecHit = theRecHits.end();
1245 
1246  for( ; it != lastRecHit ; ++it ) {
1247  theRecHitCollection.put(it->first,(it->second).begin(),(it->second).end());
1248  }
1249 
1250 }
1251 
1252 void
1254  std::map<unsigned,edm::OwnVector<SiTrackerGSMatchedRecHit2D> >& theRecHits,
1255  SiTrackerGSMatchedRecHit2DCollection& theRecHitCollection) const
1256 {
1257  std::map<unsigned,edm::OwnVector<SiTrackerGSMatchedRecHit2D> >::const_iterator
1258  it = theRecHits.begin();
1259  std::map<unsigned,edm::OwnVector<SiTrackerGSMatchedRecHit2D> >::const_iterator
1260  lastRecHit = theRecHits.end();
1261 
1262  for( ; it != lastRecHit ; ++it ) {
1263  theRecHitCollection.put(it->first,(it->second).begin(),(it->second).end());
1264  }
1265 
1266 }
1267 
1268 
1269 void
1271  std::map<unsigned, edm::OwnVector<SiTrackerGSRecHit2D> >& theRecHits,
1272  std::map<unsigned, edm::OwnVector<SiTrackerGSMatchedRecHit2D> >& matchedMap) {//,
1273  // MixCollection<PSimHit>& simhits ) {
1274 // const edm::PSimHitContainer& simhits ) { // not used in the function??? (AG)
1275 
1276 
1277  std::map<unsigned, edm::OwnVector<SiTrackerGSRecHit2D> >::iterator it = theRecHits.begin();
1278  std::map<unsigned, edm::OwnVector<SiTrackerGSRecHit2D> >::iterator lastTrack = theRecHits.end();
1279 
1280 
1281  int recHitCounter = 0;
1282 
1283 
1284  //loop over single sided tracker RecHit
1285  for( ; it != lastTrack; ++it ) {
1286 
1287  edm::OwnVector<SiTrackerGSRecHit2D>::const_iterator rit = it->second.begin();
1288  edm::OwnVector<SiTrackerGSRecHit2D>::const_iterator firstRecHit = it->second.begin();
1289  edm::OwnVector<SiTrackerGSRecHit2D>::const_iterator lastRecHit = it->second.end();
1290 
1291  //loop over rechits in track
1292  for ( ; rit != lastRecHit; ++rit,++recHitCounter){
1293 
1294  DetId detid = rit->geographicalId();
1295  unsigned int subdet = detid.subdetId();
1296  // if in the strip (subdet>2)
1297  if(subdet>2){
1298 
1299  StripSubdetector specDetId=StripSubdetector(detid);
1300 
1301  // if this is on a glued, then place only one hit in vector
1302  if(specDetId.glued()){
1303 
1304  // get the track direction from the simhit
1305  LocalVector simtrackdir = correspondingSimHit[recHitCounter]->localDirection();
1306 
1307  const GluedGeomDet* gluedDet = (const GluedGeomDet*)geometry->idToDet(DetId(specDetId.glued()));
1308  const StripGeomDetUnit* stripdet =(StripGeomDetUnit*) gluedDet->stereoDet();
1309 
1310  // global direction of track
1311  GlobalVector globaldir= stripdet->surface().toGlobal(simtrackdir);
1312  LocalVector gluedsimtrackdir=gluedDet->surface().toLocal(globaldir);
1313 
1314  // get partner layer, it is the next one or previous one in the vector
1318  partnerNext++; partnerPrev--;
1319 
1320  // check if this hit is on a stereo layer (== the second layer of a double sided module)
1321  if( specDetId.stereo() ) {
1322 
1323  int partnersFound = 0;
1324  // check next one in vector
1325  // safety check first
1326  if(partnerNext != it->second.end() )
1327  if( StripSubdetector( partnerNext->geographicalId() ).partnerDetId() == detid.rawId() ) {
1328  partner= partnerNext;
1329  partnersFound++;
1330  }
1331  // check prevoius one in vector
1332  if( rit != it->second.begin())
1333  if( StripSubdetector( partnerPrev->geographicalId() ).partnerDetId() == detid.rawId() ) {
1334  partnersFound++;
1335  partner= partnerPrev;
1336  }
1337 
1338 
1339  // in case partner has not been found this way, need to loop over all rechits in track to be sure
1340  // (rare cases fortunately)
1341  if(partnersFound==0){
1342  for(edm::OwnVector<SiTrackerGSRecHit2D>::const_iterator iter = firstRecHit; iter != lastRecHit; ++iter ){
1343  if( StripSubdetector( iter->geographicalId() ).partnerDetId() == detid.rawId()){
1344  partnersFound++;
1345  partner = iter;
1346  }
1347  }
1348  }
1349 
1350  if(partnersFound == 1) {
1351  SiTrackerGSMatchedRecHit2D * theMatchedHit =GSRecHitMatcher().match( &(*partner), &(*rit), gluedDet , gluedsimtrackdir );
1352 
1353 
1354  // std::cout << "Matched hit: isMatched =\t" << theMatchedHit->isMatched() << ", "
1355  // << theMatchedHit->monoHit() << ", " << theMatchedHit->stereoHit() << std::endl;
1356 
1357  matchedMap[it->first].push_back( theMatchedHit );
1358  }
1359  else{
1360  // no partner to match, place projected one in map
1361  SiTrackerGSMatchedRecHit2D * theProjectedHit = GSRecHitMatcher().projectOnly( &(*rit), geometry->idToDet(detid),gluedDet, gluedsimtrackdir );
1362  matchedMap[it->first].push_back( theProjectedHit );
1363  //there is no partner here
1364  }
1365  } // end if stereo
1366  else { // we are on a mono layer
1367  // usually this hit is already matched, but not if stereo hit is missing (rare cases)
1368  // check if stereo hit is missing
1369  int partnersFound = 0;
1370  // check next one in vector
1371  // safety check first
1372  if(partnerNext != it->second.end() )
1373  if( StripSubdetector( partnerNext->geographicalId() ).partnerDetId() == detid.rawId() ) {
1374  partnersFound++;
1375  }
1376  // check prevoius one in vector
1377  if( rit != it->second.begin())
1378  if( StripSubdetector( partnerPrev->geographicalId() ).partnerDetId() == detid.rawId() ) {
1379  partnersFound++;
1380  }
1381 
1382  if(partnersFound==0){ // no partner hit found this way, need to iterate over all hits to be sure (rare cases)
1383  for(edm::OwnVector<SiTrackerGSRecHit2D>::const_iterator iter = firstRecHit; iter != lastRecHit; ++iter ){
1384  if( StripSubdetector( iter->geographicalId() ).partnerDetId() == detid.rawId()){
1385  partnersFound++;
1386  }
1387  }
1388  }
1389 
1390 
1391  if(partnersFound==0){ // no partner hit found
1392  // no partner to match, place projected one one in map
1393  SiTrackerGSMatchedRecHit2D * theProjectedHit =
1394  GSRecHitMatcher().projectOnly( &(*rit), geometry->idToDet(detid),gluedDet, gluedsimtrackdir );
1395 
1396  //std::cout << "Projected hit: isMatched =\t" << theProjectedHit->isMatched() << ", "
1397  // << theProjectedHit->monoHit() << ", " << theProjectedHit->stereoHit() << std::endl;
1398 
1399  matchedMap[it->first].push_back( theProjectedHit );
1400  }
1401  } // end we are on a a mono layer
1402  } // end if glued
1403  // else matchedMap[it->first].push_back( rit->clone() ); // if not glued place the original one in vector
1404  else{ //need to copy the original in a "matched" type rechit
1405 
1406  SiTrackerGSMatchedRecHit2D* rit_copy = new SiTrackerGSMatchedRecHit2D(rit->localPosition(), rit->localPositionError(),
1407  rit->geographicalId(),
1408  rit->simhitId(), rit->simtrackId(), rit->eeId(),
1409  rit->cluster(),
1410  rit->simMultX(), rit->simMultY());
1411  //std::cout << "Simple hit hit: isMatched =\t" << rit_copy->isMatched() << ", "
1412  // << rit_copy->monoHit() << ", " << rit_copy->stereoHit() << std::endl;
1413 
1414  matchedMap[it->first].push_back( rit_copy ); // if not strip place the original one in vector (making it into a matched)
1415  }
1416 
1417  }// end if strip
1418  // else matchedMap[it->first].push_back( rit->clone() ); // if not strip place the original one in vector
1419  else { //need to copy the original in a "matched" type rechit
1420 
1421  SiTrackerGSMatchedRecHit2D* rit_copy = new SiTrackerGSMatchedRecHit2D(rit->localPosition(), rit->localPositionError(),
1422  rit->geographicalId(),
1423  rit->simhitId(), rit->simtrackId(), rit->eeId(),
1424  rit->cluster(),
1425  rit->simMultX(), rit->simMultY());
1426 
1427  //std::cout << "Simple hit hit: isMatched =\t" << rit_copy->isMatched() << ", "
1428  // << rit_copy->monoHit() << ", " << rit_copy->stereoHit() << std::endl;
1429  matchedMap[it->first].push_back( rit_copy ); // if not strip place the original one in vector (makining it into a matched)
1430  }
1431 
1432  } // end loop over rechits
1433 
1434  }// end loop over tracks
1435 }
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:114
T getParameter(std::string const &) const
static GlobalError transform(const LocalError &le, const Surface &surf)
int i
Definition: DBlmapReader.cc:9
float xx() const
Definition: LocalError.h:24
void loadRecHits(std::map< unsigned, edm::OwnVector< SiTrackerGSRecHit2D > > &theRecHits, SiTrackerGSRecHit2DCollection &theRecHitCollection) const
SiStripGaussianSmearingRecHitConverterAlgorithm * theSiStripErrorParametrization
virtual float length() const =0
unsigned int tibLayer(const DetId &id) const
unsigned int tidRing(const DetId &id) const
SiPixelGaussianSmearingRecHitConverterAlgorithm * thePixelEndcapParametrization
void smearHit(const PSimHit &simHit, double localPositionResolutionX, double localPositionResolutionY, double localPositionResolutionZ, double boundX, double boundY)
SiTrackerGSMatchedRecHit2D * match(const SiTrackerGSRecHit2D *monoRH, const SiTrackerGSRecHit2D *stereoRH, const GluedGeomDet *gluedDet, LocalVector &trackdirection) const
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
void loadClusters(std::map< unsigned, edm::OwnVector< FastTrackerCluster > > &theClusterMap, FastTrackerClusterCollection &theClusterCollection) const
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:35
static int position[TOTALCHAMBERS][3]
Definition: ReadPGInfo.cc:509
dictionary map
Definition: Association.py:205
uint32_t rawId() const
get the raw id
Definition: DetId.h:45
float xy() const
Definition: LocalError.h:25
unsigned int glued() const
glued
T mag() const
Definition: PV3DBase.h:67
void smearHits(const edm::PSimHitContainer &input, std::map< unsigned, edm::OwnVector< SiTrackerGSRecHit2D > > &theRecHits, std::map< unsigned, edm::OwnVector< FastTrackerCluster > > &theClusters, const TrackerTopology *tTopo)
float yy() const
Definition: LocalError.h:26
Local3DPoint localPosition() const
Definition: PSimHit.h:44
bool gaussianSmearing(const PSimHit &simHit, Local3DPoint &position, LocalError &error, unsigned &alphaMult, unsigned &betaMult, const TrackerTopology *tTopo)
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:94
LocalPoint toLocal(const GlobalPoint &gp) const
void loadMatchedRecHits(std::map< unsigned, edm::OwnVector< SiTrackerGSMatchedRecHit2D > > &theRecHits, SiTrackerGSMatchedRecHit2DCollection &theRecHitCollection) const
bool isAvailable() const
Definition: Service.h:47
#define end
Definition: vmac.h:38
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:39
SiPixelGaussianSmearingRecHitConverterAlgorithm * thePixelBarrelParametrization
std::vector< SiPixelQuality::disabledModuleType > * disabledModules
virtual const GeomDet * idToDet(DetId) const
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:361
RefProd< PROD > getRefBeforePut()
Definition: Event.h:106
tuple conf
Definition: dbtoconf.py:185
void put(ID id, CI begin, CI end)
insert an object range with specified identifier
Definition: RangeMap.h:119
unsigned int pxbLayer(const DetId &id) const
Point3DBase< float, LocalTag > Local3DPoint
Definition: LocalPoint.h:9
Definition: DetId.h:20
unsigned int stereo() const
stereo
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:55
SiTrackerGSMatchedRecHit2D * projectOnly(const SiTrackerGSRecHit2D *monoRH, const GeomDet *monoDet, const GluedGeomDet *gluedDet, LocalVector &ldir) const
AlignmentPositionError * alignmentPositionError() const
Return pointer to alignment errors.
Definition: GeomDet.h:88
T const * product() const
Definition: ESHandle.h:62
void matchHits(std::map< unsigned, edm::OwnVector< SiTrackerGSRecHit2D > > &theRecHits, std::map< unsigned, edm::OwnVector< SiTrackerGSMatchedRecHit2D > > &matchedMap)
const GlobalError & globalError() const
double flatShoot(double xmin=0.0, double xmax=1.0) const
Definition: RandomEngine.h:30
virtual const GeomDetUnit * idToDetUnit(DetId) const
Return the pointer to the GeomDetUnit corresponding to a given DetId.
Local3DPoint LocalPoint
Definition: LocalPoint.h:11
void smearHit(const PSimHit &simHit, const PixelGeomDetUnit *detUnit, const double boundX, const double boundY)
tuple cout
Definition: gather_cfg.py:121
std::vector< PSimHit > PSimHitContainer
Definition: Bounds.h:22
T x() const
Definition: PV3DBase.h:62
virtual float width() const =0
const GeomDetUnit * stereoDet() const
Definition: GluedGeomDet.h:21
unsigned int detUnitId() const
Definition: PSimHit.h:93
unsigned int tobLayer(const DetId &id) const
Definition: Run.h:36