CMS 3D CMS Logo

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