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