CMS 3D CMS Logo
EcalHitMaker.cc
Go to the documentation of this file.
1 #include "DataFormats/DetId/interface/DetId.h"
2 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
3 #include "DataFormats/EcalDetId/interface/EBDetId.h"
4 #include "DataFormats/EcalDetId/interface/EEDetId.h"
5 
6 //#include "Geometry/EcalAlgo/interface/EcalBarrelGeometry.h"
7 //#include "Geometry/EcalAlgo/interface/EcalEndcapGeometry.h"
8 //#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
9 
10 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
11 #include "FastSimulation/CaloHitMakers/interface/EcalHitMaker.h"
12 #include "FastSimulation/CaloGeometryTools/interface/CaloGeometryHelper.h"
13 #include "FastSimulation/CaloGeometryTools/interface/CrystalWindowMap.h"
14 #include "FastSimulation/CaloGeometryTools/interface/CaloDirectionOperations.h"
15 #include "FastSimulation/Event/interface/FSimVertex.h"
16 #include "FastSimulation/Event/interface/FSimTrack.h"
17 #include "FastSimulation/CalorimeterProperties/interface/PreshowerLayer1Properties.h"
18 #include "FastSimulation/CalorimeterProperties/interface/PreshowerLayer2Properties.h"
19 #include "FastSimulation/CalorimeterProperties/interface/HCALProperties.h"
20 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
21 //#include "FastSimulation/Utilities/interface/Histos.h"
22 
23 //#include "Math/GenVector/Transform3D.h"
24 #include "FastSimulation/CaloGeometryTools/interface/Transform3DPJ.h"
25 
26 #include <algorithm>
27 #include <cmath>
28 
29 typedef ROOT::Math::Plane3D::Vector Vector;
30 typedef ROOT::Math::Plane3D::Point Point;
31 typedef ROOT::Math::Transform3DPJ Transform3DR;
32 
33 EcalHitMaker::EcalHitMaker(CaloGeometryHelper* theCalo,
34  const XYZPoint& ecalentrance,
35  const DetId& cell,
36  int onEcal,
37  unsigned size,
38  unsigned showertype,
39  const RandomEngineAndDistribution* engine)
40  : CaloHitMaker(theCalo, DetId::Ecal, ((onEcal == 1) ? EcalBarrel : EcalEndcap), onEcal, showertype),
41  EcalEntrance_(ecalentrance),
42  onEcal_(onEcal),
43  myTrack_(nullptr),
44  random(engine) {
45 #ifdef FAMOSDEBUG
46  myHistos = Histos::instance();
47 #endif
48  // myHistos->debug("Constructeur EcalHitMaker");
49  simulatePreshower_ = true;
50  X0depthoffset_ = 0.;
51  X0PS1_ = 0.;
52  X0PS2_ = 0.;
53  X0PS2EE_ = 0.;
54  X0ECAL_ = 0.;
55  X0EHGAP_ = 0.;
56  X0HCAL_ = 0.;
57  L0PS1_ = 0.;
58  L0PS2_ = 0.;
59  L0PS2EE_ = 0.;
60  L0ECAL_ = 0.;
61  L0EHGAP_ = 0.;
62  L0HCAL_ = 0.;
63  maxX0_ = 0.;
64  totalX0_ = 0;
65  totalL0_ = 0.;
66  pulledPadProbability_ = 1.;
67  outsideWindowEnergy_ = 0.;
68  rearleakage_ = 0.;
69  bfactor_ = 1.;
70  ncrystals_ = 0;
71 
72  doreorg_ = !showertype;
73 
74  hitmaphasbeencalculated_ = false;
75 
76  if (onEcal)
77  myCalorimeter->buildCrystal(cell, pivot_);
78  else
79  pivot_ = Crystal();
80  central_ = onEcal == 1;
81  ecalFirstSegment_ = -1;
82 
83  myCrystalWindowMap_ = nullptr;
84  // In some cases, a "dummy" grid, not based on a cell, can be built. The previous variables
85  // should however be initialized. In such a case onEcal=0
86  if (!onEcal)
87  return;
88 
89  // Same size in eta-phi
90  etasize_ = size;
91  phisize_ = size;
92 
93  // Build the grid
94  // The result is put in CellsWindow and is ordered by distance to the pivot
95  myCalorimeter->getWindow(pivot_.getDetId(), size, size, CellsWindow_);
96 
97  buildGeometry();
98  // std::cout << " Geometry built " << regionOfInterest_.size() << std::endl;
99 
100  truncatedGrid_ = CellsWindow_.size() != (etasize_ * phisize_);
101 
102  // A local vector of corners
103  mycorners.resize(4);
104  corners.resize(4);
105 
106 #ifdef DEBUGGW
107  myHistos->fill("h10", EcalEntrance_.eta(), CellsWindow_.size());
108  if (onEcal == 2) {
109  myHistos->fill("h20", EcalEntrance_.perp(), CellsWindow_.size());
110  if (EcalEntrance_.perp() > 70 && EcalEntrance_.perp() < 80 && CellsWindow_.size() < 35) {
111  std::cout << " Truncated grid " << CellsWindow_.size() << " " << EcalEntrance_.perp() << std::endl;
112  std::cout << " Pivot "
113  << myCalorimeter->getEcalEndcapGeometry()->getGeometry(pivot_.getDetId())->getPosition().perp();
114  std::cout << EEDetId(pivot_.getDetId()) << std::endl;
115 
116  std::cout << " Test getClosestCell " << EcalEntrance_ << std::endl;
117  DetId testcell = myCalorimeter->getClosestCell(EcalEntrance_, true, false);
118  std::cout << " Result " << EEDetId(testcell) << std::endl;
119  std::cout << " Position " << myCalorimeter->getEcalEndcapGeometry()->getGeometry(testcell)->getPosition()
120  << std::endl;
121  }
122  }
123 
124 #endif
125 }
126 
127 EcalHitMaker::~EcalHitMaker() {
128  if (myCrystalWindowMap_ != nullptr) {
129  delete myCrystalWindowMap_;
130  }
131 }
132 
133 bool EcalHitMaker::addHitDepth(double r, double phi, double depth) {
134  // std::cout << " Add hit depth called; Current deph is " << currentdepth_;
135  // std::cout << " Required depth is " << depth << std::endl;
136  depth += X0depthoffset_;
137  double sp(1.);
138  r *= radiusFactor_;
139  CLHEP::Hep2Vector point(r * std::cos(phi), r * std::sin(phi));
140 
141  unsigned xtal = fastInsideCell(point, sp);
142  // if(cellid.isZero()) std::cout << " cell is Zero " << std::endl;
143  // if(xtal<1000)
144  // {
145  // std::cout << "Result " << regionOfInterest_[xtal].getX0Back() << " " ;
146  // std::cout << depth << std::endl;
147  // }
148  // myHistos->fill("h5000",depth);
149  if (xtal < 1000) {
150  // myHistos->fill("h5002",regionOfInterest_[xtal].getX0Back(),depth);
151  // myHistos->fill("h5003",ecalentrance_.eta(),maxX0_);
152  if (regionOfInterest_[xtal].getX0Back() > depth) {
153  hits_[xtal] += spotEnergy;
154  // myHistos->fill("h5005",r);
155  return true;
156  } else {
157  rearleakage_ += spotEnergy;
158  }
159  } else {
160  // std::cout << " Return false " << std::endl;
161  // std::cout << " Add hit depth called; Current deph is " << currentdepth_;
162  // std::cout << " Required depth is " << depth << std::endl;
163  // std::cout << " R = " << r << " " << radiusFactor_ << std::endl;
164  }
165 
166  outsideWindowEnergy_ += spotEnergy;
167  return false;
168 }
169 
170 //bool EcalHitMaker::addHitDepth(double r,double phi,double depth)
171 //{
172 // std::map<unsigned,float>::iterator itcheck=hitMap_.find(pivot_.getDetId().rawId());
173 // if(itcheck==hitMap_.end())
174 // {
175 // hitMap_.insert(std::pair<uint32_t,float>(pivot_.getDetId().rawId(),spotEnergy));
176 // }
177 // else
178 // {
179 // itcheck->second+=spotEnergy;
180 // }
181 // return true;
182 //}
183 
184 bool EcalHitMaker::addHit(double r, double phi, unsigned layer) {
185  // std::cout <<" Addhit " << std::endl;
186  // std::cout << " Before insideCell " << std::endl;
187  double sp(1.);
188  // std::cout << " Trying to add " << r << " " << phi << " " << radiusFactor_ << std::endl;
189  r *= radiusFactor_;
190  CLHEP::Hep2Vector point(r * std::cos(phi), r * std::sin(phi));
191  // std::cout << "point " << point << std::endl;
192  // CellID cellid=insideCell(point,sp);
193  unsigned xtal = fastInsideCell(point, sp);
194  // if(cellid.isZero()) std::cout << " cell is Zero " << std::endl;
195  if (xtal < 1000) {
196  if (sp == 1.)
197  hits_[xtal] += spotEnergy;
198  else
199  hits_[xtal] += (random->flatShoot() < sp) * spotEnergy;
200  return true;
201  }
202 
203  outsideWindowEnergy_ += spotEnergy;
204  // std::cout << " This hit ; r= " << point << " hasn't been added "<<std::endl;
205  // std::cout << " Xtal " << xtal << std::endl;
206  // for(unsigned ip=0;ip<npadsatdepth_;++ip)
207  // {
208  // std::cout << padsatdepth_[ip] << std::endl;
209  // }
210 
211  return false;
212 }
213 
214 // Temporary solution
215 //bool EcalHitMaker::addHit(double r,double phi,unsigned layer)
216 //{
217 // std::map<unsigned,float>::iterator itcheck=hitMap_.find(pivot_.getDetId().rawId());
218 // if(itcheck==hitMap_.end())
219 // {
220 // hitMap_.insert(std::pair<uint32_t,float>(pivot_.getDetId().rawId(),spotEnergy));
221 // }
222 // else
223 // {
224 // itcheck->second+=spotEnergy;
225 // }
226 // return true;
227 //}
228 
229 unsigned EcalHitMaker::fastInsideCell(const CLHEP::Hep2Vector& point, double& sp, bool debug) {
230  // debug = true;
231  bool found = false;
232  unsigned niter = 0;
233  // something clever has to be implemented here
234  unsigned d1, d2;
235  convertIntegerCoordinates(point.x(), point.y(), d1, d2);
236  // std::cout << "Fastinside cell " << point.x() << " " << point.y() << " " << d1 << " "<< d2 << " " << nx_ << " " << ny_ << std::endl;
237  if (d1 >= nx_ || d2 >= ny_) {
238  // std::cout << " Not in the map " <<std::endl;
239  return 9999;
240  }
241  unsigned cell = myCrystalNumberArray_[d1][d2];
242  // We are likely to be lucky
243  // std::cout << " Got the cell " << cell << std::endl;
244  if (validPads_[cell] && padsatdepth_[cell].inside(point)) {
245  // std::cout << " We are lucky " << cell << std::endl;
246  sp = padsatdepth_[cell].survivalProbability();
247  return cell;
248  }
249 
250  // std::cout << "Starting the loop " << std::endl;
251  bool status(true);
252  const std::vector<unsigned>& localCellVector(myCrystalWindowMap_->getCrystalWindow(cell, status));
253  if (status) {
254  unsigned size = localCellVector.size();
255  // std::cout << " Starting from " << EBDetId(regionOfInterest_[cell].getDetId()) << std::endl;
256  // const std::vector<DetId>& neighbours=myCalorimeter->getNeighbours(regionOfInterest_[cell].getDetId());
257  // std::cout << " The neighbours are " << std::endl;
258  // for(unsigned ic=0;ic<neighbours.size(); ++ic)
259  // {
260  // std::cout << EBDetId(neighbours[ic]) << std::endl;
261  // }
262  // std::cout << " Done " << std::endl;
263  for (unsigned ic = 0; ic < 8 && ic < size; ++ic) {
264  unsigned iq = localCellVector[ic];
265  // std::cout << " Testing " << EBDetId(regionOfInterest_[iq].getDetId()) << std::endl; ;
266  // std::cout << " " << iq << std::endl;
267  // std::cout << padsatdepth_[iq] ;
268  if (validPads_[iq] && padsatdepth_[iq].inside(point)) {
269  // std::cout << " Yes " << std::endl;
270  // myHistos->fill("h1000",niter);
271  sp = padsatdepth_[iq].survivalProbability();
272  // std::cout << "Finished the loop " << niter << std::endl;
273  // std::cout << "Inside " << std::endl;
274  return iq;
275  }
276  // std::cout << " Not inside " << std::endl;
277  // std::cout << "No " << std::endl;
278  ++niter;
279  }
280  }
281  if (debug)
282  std::cout << " not found in a quad, let's check the " << ncrackpadsatdepth_ << " cracks " << std::endl;
283  // std::cout << "Finished the loop " << niter << std::endl;
284  // Let's check the cracks
285  // std::cout << " Let's check the cracks " << ncrackpadsatdepth_ << " " << crackpadsatdepth_.size() << std::endl;
286  unsigned iquad = 0;
287  unsigned iquadinside = 999;
288  while (iquad < ncrackpadsatdepth_ && !found) {
289  // std::cout << " Inside the while " << std::endl;
290  if (crackpadsatdepth_[iquad].inside(point)) {
291  iquadinside = iquad;
292  found = true;
293  sp = crackpadsatdepth_[iquad].survivalProbability();
294  }
295  ++iquad;
296  ++niter;
297  }
298  // myHistos->fill("h1002",niter);
299  if (!found && debug)
300  std::cout << " Not found in the cracks " << std::endl;
301  return (found) ? crackpadsatdepth_[iquadinside].getNumber() : 9999;
302 }
303 
304 void EcalHitMaker::setTrackParameters(const XYZNormal& normal, double X0depthoffset, const FSimTrack& theTrack) {
305  // myHistos->debug("setTrackParameters");
306  // std::cout << " Track " << theTrack << std::endl;
307  intersections_.clear();
308  // This is certainly enough
309  intersections_.reserve(50);
310  myTrack_ = &theTrack;
311  normal_ = normal.Unit();
312  X0depthoffset_ = X0depthoffset;
313  cellLine(intersections_);
314  buildSegments(intersections_);
315  // std::cout << " Segments " << segments_.size() << std::endl;
316  // for(unsigned ii=0; ii<segments_.size() ; ++ii)
317  // {
318  // std::cout << segments_[ii] << std::endl;
319  // }
320 
321  // This is only needed in case of electromagnetic showers
322  if (EMSHOWER && onEcal_ && ecalTotalX0() > 0.) {
323  // std::cout << "Total X0 " << ecalTotalX0() << std::endl;
324  for (unsigned ic = 0; ic < ncrystals_; ++ic) {
325  for (unsigned idir = 0; idir < 4; ++idir) {
326  XYZVector norm = regionOfInterest_[ic].exitingNormal(CaloDirectionOperations::Side(idir));
327  regionOfInterest_[ic].crystalNeighbour(idir).setToBeGlued((norm.Dot(normal_) < 0.));
328  }
329  // Now calculate the distance in X0 of the back sides of the crystals
330  // (only for EM showers)
331  if (EMSHOWER) {
332  XYZVector dir = regionOfInterest_[ic].getBackCenter() - segments_[ecalFirstSegment_].entrance();
333  double dist = dir.Dot(normal_);
334  double absciss = dist + segments_[ecalFirstSegment_].sEntrance();
335  std::vector<CaloSegment>::const_iterator segiterator;
336  // First identify the correct segment
337  // std::cout << " Crystal " << ic << regionOfInterest_[ic].getBackCenter() ;
338  // std::cout << " Entrance : " << segments_[ecalFirstSegment_].entrance()<< std::endl;
339  // std::cout << " Looking for the segment " << dist << std::endl;
340  segiterator = find_if(segments_.begin(), segments_.end(), CaloSegment::inSegment(absciss));
341  // std::cout << " Done " << std::endl;
342  if (segiterator == segments_.end()) {
343  // in this case, we won't have any problem. No need to
344  // calculate the real depth.
345  regionOfInterest_[ic].setX0Back(9999);
346  } else {
347  DetId::Detector det(segiterator->whichDetector());
348  if (det != DetId::Ecal) {
349  regionOfInterest_[ic].setX0Back(9999);
350  } else {
351  double x0 = segiterator->x0FromCm(dist);
352  if (x0 < maxX0_)
353  maxX0_ = x0;
354  regionOfInterest_[ic].setX0Back(x0);
355  }
356  }
357  // myHistos->fill("h4000",ecalentrance_.eta(), regionOfInterest_[ic].getX0Back());
358  //}
359  // else
360  //{
361  // // in this case, we won't have any problem. No need to
362  // // calculate the real depth.
363  // regionOfInterest_[ic].setX0Back(9999);
364  //}
365  } //EMSHOWER
366  } // ndir
367  // myHistos->fill("h6000",segments_[ecalFirstSegment_].entrance().eta(),maxX0_);
368  }
369  // std::cout << "Leaving setTrackParameters" << std::endl
370 }
371 
372 void EcalHitMaker::cellLine(std::vector<CaloPoint>& cp) {
373  cp.clear();
374  // if(myTrack->onVFcal()!=2)
375  // {
376  if (!central_ && onEcal_ && simulatePreshower_)
377  preshowerCellLine(cp);
378  if (onEcal_)
379  ecalCellLine(EcalEntrance_, EcalEntrance_ + normal_, cp);
380  // }
381 
382  XYZPoint vertex(myTrack_->vertex().position().Vect());
383 
384  //sort the points by distance (in the ECAL they are not necessarily ordered)
385  XYZVector dir(0., 0., 0.);
386  if (myTrack_->onLayer1()) {
387  vertex = (myTrack_->layer1Entrance().vertex()).Vect();
388  dir = myTrack_->layer1Entrance().Vect().Unit();
389  } else if (myTrack_->onLayer2()) {
390  vertex = (myTrack_->layer2Entrance().vertex()).Vect();
391  dir = myTrack_->layer2Entrance().Vect().Unit();
392  } else if (myTrack_->onEcal()) {
393  vertex = (myTrack_->ecalEntrance().vertex()).Vect();
394  dir = myTrack_->ecalEntrance().Vect().Unit();
395  } else if (myTrack_->onHcal()) {
396  vertex = (myTrack_->hcalEntrance().vertex()).Vect();
397  dir = myTrack_->hcalEntrance().Vect().Unit();
398  } else if (myTrack_->onVFcal() == 2) {
399  vertex = (myTrack_->vfcalEntrance().vertex()).Vect();
400  dir = myTrack_->vfcalEntrance().Vect().Unit();
401  } else {
402  std::cout << " Problem with the grid " << std::endl;
403  }
404 
405  // Move the vertex for distance comparison (5cm)
406  vertex -= 5. * dir;
407  CaloPoint::DistanceToVertex myDistance(vertex);
408  sort(cp.begin(), cp.end(), myDistance);
409 
410  // The intersections with the HCAL shouldn't need to be sorted
411  // with the N.I it is actually a source of problems
412  hcalCellLine(cp);
413 
414  // std::cout << " Intersections ordered by distance to " << vertex << std::endl;
415  //
416  // for (unsigned ic=0;ic<cp.size();++ic)
417  // {
418  // XYZVector t=cp[ic]-vertex;
419  // std::cout << cp[ic] << " " << t.mag() << std::endl;
420  // }
421 }
422 
423 void EcalHitMaker::preshowerCellLine(std::vector<CaloPoint>& cp) const {
424  // FSimEvent& mySimEvent = myEventMgr->simSignal();
425  // FSimTrack myTrack = mySimEvent.track(fsimtrack_);
426  // std::cout << "FsimTrack " << fsimtrack_<< std::endl;
427  // std::cout << " On layer 1 " << myTrack.onLayer1() << std::endl;
428  // std::cout << " preshowerCellLine " << std::endl;
429  if (myTrack_->onLayer1()) {
430  XYZPoint point1 = (myTrack_->layer1Entrance().vertex()).Vect();
431  double phys_eta = myTrack_->layer1Entrance().eta();
432  double cmthickness = myCalorimeter->layer1Properties(1)->thickness(phys_eta);
433 
434  if (cmthickness > 0) {
435  XYZVector dir = myTrack_->layer1Entrance().Vect().Unit();
436  XYZPoint point2 = point1 + dir * cmthickness;
437 
438  CaloPoint cp1(DetId::Ecal, EcalPreshower, 1, point1);
439  CaloPoint cp2(DetId::Ecal, EcalPreshower, 1, point2);
440  cp.push_back(cp1);
441  cp.push_back(cp2);
442  } else {
443  // std::cout << "Track on ECAL " << myTrack.EcalEntrance_().vertex()*0.1<< std::endl;
444  }
445  }
446 
447  // std::cout << " On layer 2 " << myTrack.onLayer2() << std::endl;
448  if (myTrack_->onLayer2()) {
449  XYZPoint point1 = (myTrack_->layer2Entrance().vertex()).Vect();
450  double phys_eta = myTrack_->layer2Entrance().eta();
451  double cmthickness = myCalorimeter->layer2Properties(1)->thickness(phys_eta);
452  if (cmthickness > 0) {
453  XYZVector dir = myTrack_->layer2Entrance().Vect().Unit();
454  XYZPoint point2 = point1 + dir * cmthickness;
455 
456  CaloPoint cp1(DetId::Ecal, EcalPreshower, 2, point1);
457  CaloPoint cp2(DetId::Ecal, EcalPreshower, 2, point2);
458 
459  cp.push_back(cp1);
460  cp.push_back(cp2);
461  } else {
462  // std::cout << "Track on ECAL " << myTrack.EcalEntrance_().vertex()*0.1 << std::endl;
463  }
464  }
465  // std::cout << " Exit preshower CellLine " << std::endl;
466 }
467 
468 void EcalHitMaker::hcalCellLine(std::vector<CaloPoint>& cp) const {
469  // FSimEvent& mySimEvent = myEventMgr->simSignal();
470  // FSimTrack myTrack = mySimEvent.track(fsimtrack_);
471  int onHcal = myTrack_->onHcal();
472 
473  if (onHcal <= 2 && onHcal > 0) {
474  XYZPoint point1 = (myTrack_->hcalEntrance().vertex()).Vect();
475 
476  double eta = point1.eta();
477  // HCAL thickness in cm (assuming that the particle is coming from 000)
478  double thickness = myCalorimeter->hcalProperties(onHcal)->thickness(eta);
479  cp.push_back(CaloPoint(DetId::Hcal, point1));
480  XYZVector dir = myTrack_->hcalEntrance().Vect().Unit();
481  XYZPoint point2 = point1 + dir * thickness;
482 
483  cp.push_back(CaloPoint(DetId::Hcal, point2));
484  }
485  int onVFcal = myTrack_->onVFcal();
486  if (onVFcal == 2) {
487  XYZPoint point1 = (myTrack_->vfcalEntrance().vertex()).Vect();
488  double eta = point1.eta();
489  // HCAL thickness in cm (assuming that the particle is coming from 000)
490  double thickness = myCalorimeter->hcalProperties(3)->thickness(eta);
491  cp.push_back(CaloPoint(DetId::Hcal, point1));
492  XYZVector dir = myTrack_->vfcalEntrance().Vect().Unit();
493  if (thickness > 0) {
494  XYZPoint point2 = point1 + dir * thickness;
495  cp.push_back(CaloPoint(DetId::Hcal, point2));
496  }
497  }
498 }
499 
500 void EcalHitMaker::ecalCellLine(const XYZPoint& a, const XYZPoint& b, std::vector<CaloPoint>& cp) {
501  // std::vector<XYZPoint> corners;
502  // corners.resize(4);
503  unsigned ic = 0;
504  double t;
505  XYZPoint xp;
506  DetId c_entrance, c_exit;
507  bool entrancefound(false), exitfound(false);
508  // std::cout << " Look for intersections " << ncrystals_ << std::endl;
509  // std::cout << " regionOfInterest_ " << truncatedGrid_ << " " << regionOfInterest_.size() << std::endl;
510  // try to determine the number of crystals to test
511  // First determine the incident angle
512  double angle = std::acos(normal_.Dot(regionOfInterest_[0].getAxis().Unit()));
513 
514  // std::cout << " Normal " << normal_<< " Axis " << regionOfInterest_[0].getAxis().Unit() << std::endl;
515  double backdistance = std::sqrt(regionOfInterest_[0].getAxis().mag2()) * std::tan(angle);
516  // 1/2.2cm = 0.45
517  // std::cout << " Angle " << angle << std::endl;
518  // std::cout << " Back distance " << backdistance << std::endl;
519  unsigned ncrystals = (unsigned)(backdistance * 0.45);
520  unsigned highlim = (ncrystals + 4);
521  highlim *= highlim;
522  if (highlim > ncrystals_)
523  highlim = ncrystals_;
524  // unsigned lowlim=(ncrystals>2)? (ncrystals-2):0;
525  // std::cout << " Ncrys " << ncrystals << std::endl;
526 
527  while (ic < ncrystals_ && (ic < highlim || !exitfound)) {
528  // Check front side
529  // if(!entrancefound)
530  {
531  const Plane3D& plan = regionOfInterest_[ic].getFrontPlane();
532  // XYZVector axis1=(plan.Normal());
533  // XYZVector axis2=regionOfInterest_[ic].getFirstEdge();
534  xp = intersect(plan, a, b, t, false);
535  regionOfInterest_[ic].getFrontSide(corners);
536  // CrystalPad pad(9999,onEcal_,corners,regionOfInterest_[ic].getCorner(0),axis1,axis2);
537  // if(pad.globalinside(xp))
538  if (inside3D(corners, xp)) {
539  cp.push_back(CaloPoint(regionOfInterest_[ic].getDetId(), UP, xp));
540  entrancefound = true;
541  c_entrance = regionOfInterest_[ic].getDetId();
542  // myHistos->fill("j12",highlim,ic);
543  }
544  }
545 
546  // check rear side
547  // if(!exitfound)
548  {
549  const Plane3D& plan = regionOfInterest_[ic].getBackPlane();
550  // XYZVector axis1=(plan.Normal());
551  // XYZVector axis2=regionOfInterest_[ic].getFifthEdge();
552  xp = intersect(plan, a, b, t, false);
553  regionOfInterest_[ic].getBackSide(corners);
554  // CrystalPad pad(9999,onEcal_,corners,regionOfInterest_[ic].getCorner(4),axis1,axis2);
555  // if(pad.globalinside(xp))
556  if (inside3D(corners, xp)) {
557  cp.push_back(CaloPoint(regionOfInterest_[ic].getDetId(), DOWN, xp));
558  exitfound = true;
559  c_exit = regionOfInterest_[ic].getDetId();
560  // std::cout << " Crystal : " << ic << std::endl;
561  // myHistos->fill("j10",highlim,ic);
562  }
563  }
564 
565  if (entrancefound && exitfound && c_entrance == c_exit)
566  return;
567  // check lateral sides
568  for (unsigned iside = 0; iside < 4; ++iside) {
569  const Plane3D& plan = regionOfInterest_[ic].getLateralPlane(iside);
570  xp = intersect(plan, a, b, t, false);
571  // XYZVector axis1=(plan.Normal());
572  // XYZVector axis2=regionOfInterest_[ic].getLateralEdge(iside);
573  regionOfInterest_[ic].getLateralSide(iside, corners);
574  // CrystalPad pad(9999,onEcal_,corners,regionOfInterest_[ic].getCorner(iside),axis1,axis2);
575  // if(pad.globalinside(xp))
576  if (inside3D(corners, xp)) {
577  cp.push_back(CaloPoint(regionOfInterest_[ic].getDetId(), CaloDirectionOperations::Side(iside), xp));
578  // std::cout << cp[cp.size()-1] << std::endl;
579  }
580  }
581  // Go to next crystal
582  ++ic;
583  }
584 }
585 
586 void EcalHitMaker::buildSegments(const std::vector<CaloPoint>& cp) {
587  // myHistos->debug();
588  // TimeMe theT("FamosGrid::buildSegments");
589  unsigned size = cp.size();
590  if (size % 2 != 0) {
591  // std::cout << " There is a problem " << std::endl;
592  return;
593  }
594  // myHistos->debug();
595  unsigned nsegments = size / 2;
596  segments_.reserve(nsegments);
597  if (size == 0)
598  return;
599  // curv abs
600  double s = 0.;
601  double sX0 = 0.;
602  double sL0 = 0.;
603 
604  unsigned ncrossedxtals = 0;
605  unsigned is = 0;
606  while (is < nsegments) {
607  if (cp[2 * is].getDetId() != cp[2 * is + 1].getDetId() && cp[2 * is].whichDetector() != DetId::Hcal &&
608  cp[2 * is + 1].whichDetector() != DetId::Hcal) {
609  // std::cout << " Problem with the segments " << std::endl;
610  // std::cout << cp[2*is].whichDetector() << " " << cp[2*is+1].whichDetector() << std::endl;
611  // std::cout << is << " " <<cp[2*is].getDetId().rawId() << std::endl;
612  // std::cout << (2*is+1) << " " <<cp[2*is+1].getDetId().rawId() << std::endl;
613  ++is;
614  continue;
615  }
616 
617  // Check if it is a Preshower segment - Layer 1
618  // One segment per layer, nothing between
619  // myHistos->debug("Just avant Preshower");
620  if (cp[2 * is].whichDetector() == DetId::Ecal && cp[2 * is].whichSubDetector() == EcalPreshower &&
621  cp[2 * is].whichLayer() == 1) {
622  if (cp[2 * is + 1].whichDetector() == DetId::Ecal && cp[2 * is + 1].whichSubDetector() == EcalPreshower &&
623  cp[2 * is + 1].whichLayer() == 1) {
624  CaloSegment preshsegment(cp[2 * is], cp[2 * is + 1], s, sX0, sL0, CaloSegment::PS, myCalorimeter);
625  segments_.push_back(preshsegment);
626  // std::cout << " Added (1-1)" << preshsegment << std::endl;
627  s += preshsegment.length();
628  sX0 += preshsegment.X0length();
629  sL0 += preshsegment.L0length();
630  X0PS1_ += preshsegment.X0length();
631  L0PS1_ += preshsegment.L0length();
632  } else {
633  std::cout << " Strange segment between Preshower1 and " << cp[2 * is + 1].whichDetector();
634  std::cout << std::endl;
635  }
636  ++is;
637  continue;
638  }
639 
640  // Check if it is a Preshower segment - Layer 2
641  // One segment per layer, nothing between
642  if (cp[2 * is].whichDetector() == DetId::Ecal && cp[2 * is].whichSubDetector() == EcalPreshower &&
643  cp[2 * is].whichLayer() == 2) {
644  if (cp[2 * is + 1].whichDetector() == DetId::Ecal && cp[2 * is + 1].whichSubDetector() == EcalPreshower &&
645  cp[2 * is + 1].whichLayer() == 2) {
646  CaloSegment preshsegment(cp[2 * is], cp[2 * is + 1], s, sX0, sL0, CaloSegment::PS, myCalorimeter);
647  segments_.push_back(preshsegment);
648  // std::cout << " Added (1-2)" << preshsegment << std::endl;
649  s += preshsegment.length();
650  sX0 += preshsegment.X0length();
651  sL0 += preshsegment.L0length();
652  X0PS2_ += preshsegment.X0length();
653  L0PS2_ += preshsegment.L0length();
654 
655  // material between preshower and EE
656  if (is < (nsegments - 1) && cp[2 * is + 2].whichDetector() == DetId::Ecal &&
657  cp[2 * is + 2].whichSubDetector() == EcalEndcap) {
658  CaloSegment gapsef(cp[2 * is + 1], cp[2 * is + 2], s, sX0, sL0, CaloSegment::PSEEGAP, myCalorimeter);
659  segments_.push_back(gapsef);
660  s += gapsef.length();
661  sX0 += gapsef.X0length();
662  sL0 += gapsef.L0length();
663  X0PS2EE_ += gapsef.X0length();
664  L0PS2EE_ += gapsef.L0length();
665  // std::cout << " Created a segment " << gapsef.length()<< " " << gapsef.X0length()<< std::endl;
666  }
667  } else {
668  std::cout << " Strange segment between Preshower2 and " << cp[2 * is + 1].whichDetector();
669  std::cout << std::endl;
670  }
671  ++is;
672  continue;
673  }
674  // Now deal with the ECAL
675  // One segment in each crystal. Segment corresponding to cracks/gaps are added
676  // myHistos->debug("Just avant ECAL");
677  if (cp[2 * is].whichDetector() == DetId::Ecal &&
678  (cp[2 * is].whichSubDetector() == EcalBarrel || cp[2 * is].whichSubDetector() == EcalEndcap)) {
679  if (cp[2 * is + 1].whichDetector() == DetId::Ecal &&
680  (cp[2 * is + 1].whichSubDetector() == EcalBarrel || cp[2 * is + 1].whichSubDetector() == EcalEndcap)) {
681  DetId cell2 = cp[2 * is + 1].getDetId();
682  // set the real entrance
683  if (ecalFirstSegment_ < 0)
684  ecalFirstSegment_ = segments_.size();
685 
686  // !! Approximatiom : the first segment is always in a crystal
687  if (cp[2 * is].getDetId() == cell2) {
688  CaloSegment segment(cp[2 * is], cp[2 * is + 1], s, sX0, sL0, CaloSegment::PbWO4, myCalorimeter);
689  segments_.push_back(segment);
690  // std::cout << " Added (2)" << segment << std::endl;
691  s += segment.length();
692  sX0 += segment.X0length();
693  sL0 += segment.L0length();
694  X0ECAL_ += segment.X0length();
695  L0ECAL_ += segment.L0length();
696  ++ncrossedxtals;
697  ++is;
698  } else {
699  std::cout << " One more bug in the segment " << std::endl;
700  ++is;
701  }
702  // Now check if a gap or crack should be added
703  if (is > 0 && is < nsegments) {
704  DetId cell3 = cp[2 * is].getDetId();
705  if (cp[2 * is].whichDetector() != DetId::Hcal) {
706  // Crack inside the ECAL
707  bool bordercrossing = myCalorimeter->borderCrossing(cell2, cell3);
708  CaloSegment cracksegment(cp[2 * is - 1],
709  cp[2 * is],
710  s,
711  sX0,
712  sL0,
713  (bordercrossing) ? CaloSegment::CRACK : CaloSegment::GAP,
714  myCalorimeter);
715  segments_.push_back(cracksegment);
716  s += cracksegment.length();
717  sX0 += cracksegment.X0length();
718  sL0 += cracksegment.L0length();
719  X0ECAL_ += cracksegment.X0length();
720  L0ECAL_ += cracksegment.L0length();
721  // std::cout <<" Added(3) "<< cracksegment << std::endl;
722  } else {
723  // a segment corresponding to ECAL/HCAL transition should be
724  // added here
725  CaloSegment cracksegment(cp[2 * is - 1], cp[2 * is], s, sX0, sL0, CaloSegment::ECALHCALGAP, myCalorimeter);
726  segments_.push_back(cracksegment);
727  s += cracksegment.length();
728  sX0 += cracksegment.X0length();
729  sL0 += cracksegment.L0length();
730  X0EHGAP_ += cracksegment.X0length();
731  L0EHGAP_ += cracksegment.L0length();
732  }
733  }
734  continue;
735  } else {
736  std::cout << " Strange segment between " << cp[2 * is].whichDetector();
737  std::cout << " and " << cp[2 * is + 1].whichDetector() << std::endl;
738  ++is;
739  continue;
740  }
741  }
742  // myHistos->debug("Just avant HCAL");
743  // HCAL
744  if (cp[2 * is].whichDetector() == DetId::Hcal && cp[2 * is + 1].whichDetector() == DetId::Hcal) {
745  CaloSegment segment(cp[2 * is], cp[2 * is + 1], s, sX0, sL0, CaloSegment::HCAL, myCalorimeter);
746  segments_.push_back(segment);
747  s += segment.length();
748  sX0 += segment.X0length();
749  sL0 += segment.L0length();
750  X0HCAL_ += segment.X0length();
751  L0HCAL_ += segment.L0length();
752  // std::cout <<" Added(4) "<< segment << std::endl;
753  ++is;
754  }
755  }
756  // std::cout << " PS1 " << X0PS1_ << " " << L0PS1_ << std::endl;
757  // std::cout << " PS2 " << X0PS2_ << " " << L0PS2_ << std::endl;
758  // std::cout << " ECAL " << X0ECAL_ << " " << L0ECAL_ << std::endl;
759  // std::cout << " HCAL " << X0HCAL_ << " " << L0HCAL_ << std::endl;
760 
761  totalX0_ = X0PS1_ + X0PS2_ + X0PS2EE_ + X0ECAL_ + X0EHGAP_ + X0HCAL_;
762  totalL0_ = L0PS1_ + L0PS2_ + L0PS2EE_ + L0ECAL_ + L0EHGAP_ + L0HCAL_;
763  // myHistos->debug("Just avant le fill");
764 
765 #ifdef DEBUGCELLLINE
766  myHistos->fill("h200", fabs(EcalEntrance_.eta()), X0ECAL_);
767  myHistos->fill("h210", EcalEntrance_.phi(), X0ECAL_);
768  if (X0ECAL_ < 20)
769  myHistos->fill("h212", EcalEntrance_.phi(), X0ECAL_);
770  // if(X0ECAL_<1.)
771  // {
772  // for(unsigned ii=0; ii<segments_.size() ; ++ii)
773  // {
774  // std::cout << segments_[ii] << std::endl;
775  // }
776  // }
777  myHistos->fillByNumber("h30", ncrossedxtals, EcalEntrance_.eta(), X0ECAL_);
778 
779  double zvertex = myTrack_->vertex().position().z();
780 
781  myHistos->fill("h310", EcalEntrance_.eta(), X0ECAL_);
782  if (X0ECAL_ < 22)
783  myHistos->fill("h410", EcalEntrance_.phi());
784  myHistos->fill("h400", zvertex, X0ECAL_);
785 #endif
786  // std::cout << " Finished the segments " << std::endl;
787 }
788 
789 void EcalHitMaker::buildGeometry() {
790  configuredGeometry_ = false;
791  ncrystals_ = CellsWindow_.size();
792  // create the vector with of pads with the appropriate size
793  padsatdepth_.resize(ncrystals_);
794 
795  // This is fully correct in the barrel.
796  ny_ = phisize_;
797  nx_ = ncrystals_ / ny_;
798  std::vector<unsigned> empty;
799  empty.resize(ny_, 0);
800  myCrystalNumberArray_.reserve((unsigned)nx_);
801  for (unsigned inx = 0; inx < (unsigned)nx_; ++inx) {
802  myCrystalNumberArray_.push_back(empty);
803  }
804 
805  hits_.resize(ncrystals_, 0.);
806  regionOfInterest_.clear();
807  regionOfInterest_.resize(ncrystals_);
808  validPads_.resize(ncrystals_);
809  for (unsigned ic = 0; ic < ncrystals_; ++ic) {
810  myCalorimeter->buildCrystal(CellsWindow_[ic], regionOfInterest_[ic]);
811  regionOfInterest_[ic].setNumber(ic);
812  DetIdMap_.insert(std::pair<DetId, unsigned>(CellsWindow_[ic], ic));
813  }
814 
815  // Computes the map of the neighbours
816  myCrystalWindowMap_ = new CrystalWindowMap(myCalorimeter, regionOfInterest_);
817 }
818 
819 // depth is in X0 , L0 (depending on EMSHOWER/HADSHOWER) or in CM if inCM
820 bool EcalHitMaker::getPads(double depth, bool inCm) {
821  //std::cout << " New depth " << depth << std::endl;
822  // The first time, the relationship between crystals must be calculated
823  // but only in the case of EM showers
824 
825  if (EMSHOWER && !configuredGeometry_)
826  configureGeometry();
827 
828  radiusFactor_ = (EMSHOWER) ? moliereRadius * radiusCorrectionFactor_ : interactionLength;
829  detailedShowerTail_ = false;
830  if (EMSHOWER)
831  currentdepth_ = depth + X0depthoffset_;
832  else
833  currentdepth_ = depth;
834 
835  // if(currentdepth_>maxX0_+ps1TotalX0()+ps2TotalX0())
836  // {
837  // currentdepth_=maxX0_+ps1TotalX0()+ps2TotalX0()-1.; // the -1 is for safety
838  // detailedShowerTail_=true;
839  // }
840 
841  // std::cout << " FamosGrid::getQuads " << currentdepth_ << " " << maxX0_ << std::endl;
842 
843  ncrackpadsatdepth_ = 0;
844 
845  xmin_ = ymin_ = 999;
846  xmax_ = ymax_ = -999;
847  double locxmin, locxmax, locymin, locymax;
848 
849  // Get the depth of the pivot
850  std::vector<CaloSegment>::const_iterator segiterator;
851  // First identify the correct segment
852 
853  if (inCm) // centimeter
854  {
855  segiterator = find_if(segments_.begin(), segments_.end(), CaloSegment::inSegment(currentdepth_));
856  } else {
857  // EM shower
858  if (EMSHOWER)
859  segiterator = find_if(segments_.begin(), segments_.end(), CaloSegment::inX0Segment(currentdepth_));
860 
861  //Hadron shower
862  if (HADSHOWER)
863  segiterator = find_if(segments_.begin(), segments_.end(), CaloSegment::inL0Segment(currentdepth_));
864  }
865  if (segiterator == segments_.end()) {
866  std::cout << " FamosGrid: Could not go at such depth " << depth << std::endl;
867  std::cout << " EMSHOWER " << EMSHOWER << std::endl;
868  std::cout << " Track " << *myTrack_ << std::endl;
869  std::cout << " Segments " << segments_.size() << std::endl;
870  for (unsigned ii = 0; ii < segments_.size(); ++ii) {
871  std::cout << segments_[ii] << std::endl;
872  }
873 
874  return false;
875  }
876  // std::cout << *segiterator << std::endl;
877 
878  if (segiterator->whichDetector() != DetId::Ecal) {
879  std::cout << " In " << segiterator->whichDetector() << std::endl;
880  // buildPreshower();
881  return false;
882  }
883 
884  // std::cout << *segiterator << std::endl;
885  // get the position of the origin
886 
887  XYZPoint origin;
888  if (inCm) {
889  origin = segiterator->positionAtDepthincm(currentdepth_);
890  } else {
891  if (EMSHOWER)
892  origin = segiterator->positionAtDepthinX0(currentdepth_);
893  if (HADSHOWER)
894  origin = segiterator->positionAtDepthinL0(currentdepth_);
895  }
896  // std::cout << " currentdepth_ " << currentdepth_ << " " << origin << std::endl;
897  XYZVector newaxis = pivot_.getFirstEdge().Cross(normal_);
898 
899  // std::cout << "Normal " << normal_ << std::endl;
900  // std::cout << " New axis " << newaxis << std::endl;
901 
902  // std::cout << " ncrystals " << ncrystals << std::endl;
903  plan_ = Plane3D((Vector)normal_, (Point)origin);
904 
905  unsigned nquads = 0;
906  double sign = (central_) ? -1. : 1.;
907  Transform3DR trans((Point)origin,
908  (Point)(origin + normal_),
909  (Point)(origin + newaxis),
910  Point(0, 0, 0),
911  Point(0., 0., sign),
912  Point(0., 1., 0.));
913  for (unsigned ic = 0; ic < ncrystals_; ++ic) {
914  // std::cout << " Building geometry for " << regionOfInterest_[ic].getCellID() << std::endl;
915  XYZPoint a, b;
916 
917  // std::cout << " Origin " << origin << std::endl;
918 
919  // std::vector<XYZPoint> corners;
920  // corners.reserve(4);
921  double dummyt;
922  bool hasbeenpulled = false;
923  bool behindback = false;
924  for (unsigned il = 0; il < 4; ++il) {
925  // a is the il-th front corner of the crystal. b is the corresponding rear corner
926  regionOfInterest_[ic].getLateralEdges(il, a, b);
927 
928  // pull the surface if necessary (only in the front of the crystals)
929  XYZPoint aprime = a;
930  if (pulled(origin, normal_, a)) {
931  b = aprime;
932  hasbeenpulled = true;
933  }
934 
935  // compute the intersection.
936  // Check that the intersection is in the [a,b] segment if HADSHOWER
937  // if EMSHOWER the intersection is calculated as if the crystals were infinite
938  XYZPoint xx = (EMSHOWER) ? intersect(plan_, a, b, dummyt, false) : intersect(plan_, a, b, dummyt, true);
939 
940  if (dummyt > 1)
941  behindback = true;
942  // std::cout << " Intersect " << il << " " << a << " " << b << " " << plan_ << " " << xx << std::endl;
943  // check that the intersection actually exists
944  if (xx.mag2() != 0) {
945  corners[il] = xx;
946  }
947  }
948  // std::cout << " ncorners " << corners.size() << std::endl;
949  if (behindback && EMSHOWER)
950  detailedShowerTail_ = true;
951  // If the quad is completly defined. Store it !
952  if (corners.size() == 4) {
953  padsatdepth_[ic] = CrystalPad(ic, corners, trans, bfactor_, !central_);
954  // Parameter to be tuned
955  if (hasbeenpulled)
956  padsatdepth_[ic].setSurvivalProbability(pulledPadProbability_);
957  validPads_[ic] = true;
958  ++nquads;
959  // In principle, this should be done after the quads reorganization. But it would cost one more loop
960  // quadsatdepth_[ic].extrems(locxmin,locxmax,locymin,locymax);
961  // if(locxmin<xmin_) xmin_=locxmin;
962  // if(locymin<ymin_) ymin_=locymin;
963  // if(locxmax>xmax_) xmax_=locxmax;
964  // if(locymax>ymax_) ymax_=locymax;
965  } else {
966  padsatdepth_[ic] = CrystalPad();
967  validPads_[ic] = false;
968  }
969  }
970  // std::cout << " Number of quads " << quadsatdepth_.size() << std::endl;
971  if (doreorg_)
972  reorganizePads();
973  // std::cout << "Finished to reorganize " << std::endl;
974  npadsatdepth_ = nquads;
975  // std::cout << " prepareCellIDMap " << std::endl;
976 
977  // Resize the Quads to allow for some numerical inaccuracy
978  // in the "inside" function
979  for (unsigned ic = 0; ic < ncrystals_; ++ic) {
980  if (!validPads_[ic])
981  continue;
982 
983  if (EMSHOWER)
984  padsatdepth_[ic].resetCorners();
985 
986  padsatdepth_[ic].extrems(locxmin, locxmax, locymin, locymax);
987  if (locxmin < xmin_)
988  xmin_ = locxmin;
989  if (locymin < ymin_)
990  ymin_ = locymin;
991  if (locxmax > xmax_)
992  xmax_ = locxmax;
993  if (locymax > ymax_)
994  ymax_ = locymax;
995  }
996 
997  sizex_ = (xmax_ - xmin_) / nx_;
998  sizey_ = (ymax_ - ymin_) / ny_;
999 
1000  // Make sure that sizex_ and sizey_ are set before running prepareCellIDMap
1001  prepareCrystalNumberArray();
1002 
1003  // std::cout << " Finished prepareCellIDMap " << std::endl;
1004  ncrackpadsatdepth_ = crackpadsatdepth_.size();
1005 
1006  return true;
1007 }
1008 
1009 void EcalHitMaker::configureGeometry() {
1010  configuredGeometry_ = true;
1011  for (unsigned ic = 0; ic < ncrystals_; ++ic) {
1012  // std::cout << " Building " << cellids_[ic] << std::endl;
1013  for (unsigned idir = 0; idir < 8; ++idir) {
1014  unsigned oppdir = CaloDirectionOperations::oppositeDirection(idir);
1015  // Is there something else to do ?
1016  // The relationship with the neighbour may have been set previously.
1017  if (regionOfInterest_[ic].crystalNeighbour(idir).status() >= 0) {
1018  // std::cout << " Nothing to do " << std::endl;
1019  continue;
1020  }
1021 
1022  const DetId& oldcell(regionOfInterest_[ic].getDetId());
1023  CaloDirection dir = CaloDirectionOperations::neighbourDirection(idir);
1024  DetId newcell(oldcell);
1025  if (!myCalorimeter->move(newcell, dir)) {
1026  // no neighbour in this direction
1027  regionOfInterest_[ic].crystalNeighbour(idir).setStatus(-1);
1028  continue;
1029  }
1030  // Determine the number of this neighbour
1031  // std::cout << " The neighbour is " << newcell << std::endl;
1032  std::map<DetId, unsigned>::const_iterator niter(DetIdMap_.find(newcell));
1033  if (niter == DetIdMap_.end()) {
1034  // std::cout << " The neighbour is not in the map " << std::endl;
1035  regionOfInterest_[ic].crystalNeighbour(idir).setStatus(-1);
1036  continue;
1037  }
1038  // Now there is a neighbour
1039  // std::cout << " The neighbour is " << niter->second << " " << cellids_[niter->second] << std::endl;
1040  regionOfInterest_[ic].crystalNeighbour(idir).setNumber(niter->second);
1041  // std::cout << " Managed to set crystalNeighbour " << ic << " " << idir << std::endl;
1042  // std::cout << " Trying " << niter->second << " " << oppdir << std::endl;
1043  regionOfInterest_[niter->second].crystalNeighbour(oppdir).setNumber(ic);
1044  // std::cout << " Crack/gap " << std::endl;
1045  if (myCalorimeter->borderCrossing(oldcell, newcell)) {
1046  regionOfInterest_[ic].crystalNeighbour(idir).setStatus(1);
1047  regionOfInterest_[niter->second].crystalNeighbour(oppdir).setStatus(1);
1048  // std::cout << " Crack ! " << std::endl;
1049  } else {
1050  regionOfInterest_[ic].crystalNeighbour(idir).setStatus(0);
1051  regionOfInterest_[niter->second].crystalNeighbour(oppdir).setStatus(0);
1052  // std::cout << " Gap" << std::endl;
1053  }
1054  }
1055  }
1056  // Magnetic field a la Charlot
1057  double theta = EcalEntrance_.theta();
1058  if (theta > M_PI_2)
1059  theta = M_PI - theta;
1060  bfactor_ = 1. / (1. + 0.133 * theta);
1061  if (myCalorimeter->magneticField() == 0.)
1062  bfactor_ = 1.;
1063 }
1064 
1065 // project fPoint on the plane (original,normal)
1066 bool EcalHitMaker::pulled(const XYZPoint& origin, const XYZNormal& normal, XYZPoint& fPoint) const {
1067  // check if fPoint is behind the origin
1068  double dotproduct = normal.Dot(fPoint - origin);
1069  if (dotproduct <= 0.)
1070  return false;
1071 
1072  //norm of normal is 1
1073  fPoint -= (1 + dotproduct) * normal;
1074  return true;
1075 }
1076 
1077 void EcalHitMaker::prepareCrystalNumberArray() {
1078  for (unsigned iq = 0; iq < npadsatdepth_; ++iq) {
1079  if (!validPads_[iq])
1080  continue;
1081  unsigned d1, d2;
1082  convertIntegerCoordinates(padsatdepth_[iq].center().x(), padsatdepth_[iq].center().y(), d1, d2);
1083  myCrystalNumberArray_[d1][d2] = iq;
1084  }
1085 }
1086 
1087 void EcalHitMaker::convertIntegerCoordinates(double x, double y, unsigned& ix, unsigned& iy) const {
1088  int tix = (int)((x - xmin_) / sizex_);
1089  int tiy = (int)((y - ymin_) / sizey_);
1090  ix = iy = 9999;
1091  if (tix >= 0)
1092  ix = (unsigned)tix;
1093  if (tiy >= 0)
1094  iy = (unsigned)tiy;
1095 }
1096 
1097 const std::map<CaloHitID, float>& EcalHitMaker::getHits() {
1098  if (hitmaphasbeencalculated_)
1099  return hitMap_;
1100  for (unsigned ic = 0; ic < ncrystals_; ++ic) {
1101  //calculate time of flight
1102  float tof = 0.0;
1103  if (onEcal_ == 1 || onEcal_ == 2)
1104  tof =
1105  (myCalorimeter->getEcalGeometry(onEcal_)->getGeometry(regionOfInterest_[ic].getDetId())->getPosition().mag()) /
1106  29.98; //speed of light
1107 
1108  if (onEcal_ == 1) {
1109  CaloHitID current_id(EBDetId(regionOfInterest_[ic].getDetId().rawId()).hashedIndex(), tof, 0); //no track yet
1110  hitMap_.insert(std::pair<CaloHitID, float>(current_id, hits_[ic]));
1111  } else if (onEcal_ == 2) {
1112  CaloHitID current_id(EEDetId(regionOfInterest_[ic].getDetId().rawId()).hashedIndex(), tof, 0); //no track yet
1113  hitMap_.insert(std::pair<CaloHitID, float>(current_id, hits_[ic]));
1114  }
1115  }
1116  hitmaphasbeencalculated_ = true;
1117  return hitMap_;
1118 }
1119 
1121 
1122 void EcalHitMaker::reorganizePads() {
1123  // Some cleaning first
1124  // std::cout << " Starting reorganize " << std::endl;
1125  crackpadsatdepth_.clear();
1126  crackpadsatdepth_.reserve(etasize_ * phisize_);
1127  ncrackpadsatdepth_ = 0;
1128  std::vector<neighbour> gaps;
1129  std::vector<std::vector<neighbour> > cracks;
1130  // cracks.clear();
1131  cracks.resize(ncrystals_);
1132 
1133  for (unsigned iq = 0; iq < ncrystals_; ++iq) {
1134  if (!validPads_[iq])
1135  continue;
1136 
1137  gaps.clear();
1138  // std::cout << padsatdepth_[iq] << std::endl;
1139  //check all the directions
1140  for (unsigned iside = 0; iside < 4; ++iside) {
1141  // std::cout << " To be glued " << iside << " " << regionOfInterest_[iq].crystalNeighbour(iside).toBeGlued() << std::endl;
1142  CaloDirection thisside = CaloDirectionOperations::Side(iside);
1143  if (regionOfInterest_[iq].crystalNeighbour(iside).toBeGlued()) {
1144  // look for the neighbour and check that it exists
1145  int neighbourstatus = regionOfInterest_[iq].crystalNeighbour(iside).status();
1146  if (neighbourstatus < 0)
1147  continue;
1148 
1149  unsigned neighbourNumber = regionOfInterest_[iq].crystalNeighbour(iside).number();
1150  if (!validPads_[neighbourNumber])
1151  continue;
1152  // there is a crack between
1153  if (neighbourstatus == 1) {
1154  // std::cout << " 1 Crack : " << thisside << " " << cellids_[iq]<< " " << cellids_[neighbourNumber] << std::endl;
1155  cracks[iq].push_back(neighbour(thisside, neighbourNumber));
1156  } // else it is a gap
1157  else {
1158  gaps.push_back(neighbour(thisside, neighbourNumber));
1159  }
1160  }
1161  }
1162  // Now lift the gaps
1163  gapsLifting(gaps, iq);
1164  }
1165 
1166  unsigned ncracks = cracks.size();
1167  // std::cout << " Cracks treatment : " << cracks.size() << std::endl;
1168  for (unsigned icrack = 0; icrack < ncracks; ++icrack) {
1169  // std::cout << " Crack number " << crackiter->first << std::endl;
1170  cracksPads(cracks[icrack], icrack);
1171  }
1172 }
1173 
1174 //dir 2 = N,E,W,S
1175 CLHEP::Hep2Vector& EcalHitMaker::correspondingEdge(neighbour& myneighbour, CaloDirection dir2) {
1176  CaloDirection dir = CaloDirectionOperations::oppositeSide(myneighbour.first);
1177  CaloDirection corner = CaloDirectionOperations::add2d(dir, dir2);
1178  // std::cout << "Corresponding Edge " << dir<< " " << dir2 << " " << corner << std::endl;
1179  return padsatdepth_[myneighbour.second].edge(corner);
1180 }
1181 
1182 bool EcalHitMaker::diagonalEdge(unsigned myPad, CaloDirection dir, CLHEP::Hep2Vector& point) {
1183  unsigned idir = CaloDirectionOperations::neighbourDirection(dir);
1184  if (regionOfInterest_[myPad].crystalNeighbour(idir).status() < 0)
1185  return false;
1186  unsigned nneighbour = regionOfInterest_[myPad].crystalNeighbour(idir).number();
1187  if (!validPads_[nneighbour]) {
1188  // std::cout << " Wasn't able to move " << std::endl;
1189  return false;
1190  }
1191  point = padsatdepth_[nneighbour].edge(CaloDirectionOperations::oppositeSide(dir));
1192  return true;
1193 }
1194 
1195 bool EcalHitMaker::unbalancedDirection(const std::vector<neighbour>& dirs,
1196  unsigned& unb,
1197  unsigned& dir1,
1198  unsigned& dir2) {
1199  if (dirs.size() == 1)
1200  return false;
1201  if (dirs.size() % 2 == 0)
1202  return false;
1203  CaloDirection tmp;
1204  tmp = CaloDirectionOperations::add2d(dirs[0].first, dirs[1].first);
1205  if (tmp == NONE) {
1206  unb = 2;
1207  dir1 = 0;
1208  dir2 = 1;
1209  return true;
1210  }
1211  tmp = CaloDirectionOperations::add2d(dirs[0].first, dirs[2].first);
1212  if (tmp == NONE) {
1213  unb = 1;
1214  dir1 = 0;
1215  dir2 = 2;
1216  return true;
1217  }
1218  unb = 0;
1219  dir1 = 1;
1220  dir2 = 2;
1221  return true;
1222 }
1223 
1224 void EcalHitMaker::gapsLifting(std::vector<neighbour>& gaps, unsigned iq) {
1225  // std::cout << " Entering gapsLifting " << std::endl;
1226  CrystalPad& myPad = padsatdepth_[iq];
1227  unsigned ngaps = gaps.size();
1228  constexpr bool debug = false;
1229  if (ngaps == 1) {
1230  if (debug) {
1231  std::cout << " Avant " << ngaps << " " << gaps[0].first << std::endl;
1232  std::cout << myPad << std::endl;
1233  }
1234  if (gaps[0].first == NORTH || gaps[0].first == SOUTH) {
1235  CaloDirection dir1 = CaloDirectionOperations::add2d(gaps[0].first, EAST);
1236  CaloDirection dir2 = CaloDirectionOperations::add2d(gaps[0].first, WEST);
1237  myPad.edge(dir1) = correspondingEdge(gaps[0], EAST);
1238  myPad.edge(dir2) = correspondingEdge(gaps[0], WEST);
1239  } else {
1240  CaloDirection dir1 = CaloDirectionOperations::add2d(gaps[0].first, NORTH);
1241  CaloDirection dir2 = CaloDirectionOperations::add2d(gaps[0].first, SOUTH);
1242  myPad.edge(dir1) = correspondingEdge(gaps[0], NORTH);
1243  myPad.edge(dir2) = correspondingEdge(gaps[0], SOUTH);
1244  }
1245  if (debug) {
1246  std::cout << " Apres " << std::endl;
1247  std::cout << myPad << std::endl;
1248  }
1249  } else if (ngaps == 2) {
1250  if (debug) {
1251  std::cout << " Avant " << ngaps << " " << gaps[0].first << " " << gaps[1].first << std::endl;
1252  std::cout << myPad << std::endl;
1253  std::cout << " Voisin 1 " << (gaps[0].second) << std::endl;
1254  std::cout << " Voisin 2 " << (gaps[1].second) << std::endl;
1255  }
1256  CaloDirection corner0 = CaloDirectionOperations::add2d(gaps[0].first, gaps[1].first);
1257 
1258  CLHEP::Hep2Vector point(0., 0.);
1259  if (corner0 != NONE && diagonalEdge(iq, corner0, point)) {
1260  CaloDirection corner1 =
1261  CaloDirectionOperations::add2d(CaloDirectionOperations::oppositeSide(gaps[0].first), gaps[1].first);
1262  CaloDirection corner2 =
1263  CaloDirectionOperations::add2d(gaps[0].first, CaloDirectionOperations::oppositeSide(gaps[1].first));
1264  myPad.edge(corner0) = point;
1265  myPad.edge(corner1) = correspondingEdge(gaps[1], CaloDirectionOperations::oppositeSide(gaps[0].first));
1266  myPad.edge(corner2) = correspondingEdge(gaps[0], CaloDirectionOperations::oppositeSide(gaps[1].first));
1267  } else if (corner0 == NONE) {
1268  if (gaps[0].first == EAST || gaps[0].first == WEST) {
1269  CaloDirection corner1 = CaloDirectionOperations::add2d(gaps[0].first, NORTH);
1270  CaloDirection corner2 = CaloDirectionOperations::add2d(gaps[0].first, SOUTH);
1271  myPad.edge(corner1) = correspondingEdge(gaps[0], NORTH);
1272  myPad.edge(corner2) = correspondingEdge(gaps[0], SOUTH);
1273 
1274  corner1 = CaloDirectionOperations::add2d(gaps[1].first, NORTH);
1275  corner2 = CaloDirectionOperations::add2d(gaps[1].first, SOUTH);
1276  myPad.edge(corner1) = correspondingEdge(gaps[1], NORTH);
1277  myPad.edge(corner2) = correspondingEdge(gaps[1], SOUTH);
1278  } else {
1279  CaloDirection corner1 = CaloDirectionOperations::add2d(gaps[0].first, EAST);
1280  CaloDirection corner2 = CaloDirectionOperations::add2d(gaps[0].first, WEST);
1281  myPad.edge(corner1) = correspondingEdge(gaps[0], EAST);
1282  myPad.edge(corner2) = correspondingEdge(gaps[0], WEST);
1283 
1284  corner1 = CaloDirectionOperations::add2d(gaps[1].first, EAST);
1285  corner2 = CaloDirectionOperations::add2d(gaps[1].first, WEST);
1286  myPad.edge(corner1) = correspondingEdge(gaps[1], EAST);
1287  myPad.edge(corner2) = correspondingEdge(gaps[1], WEST);
1288  }
1289  }
1290  if (debug) {
1291  std::cout << " Apres " << std::endl;
1292  std::cout << myPad << std::endl;
1293  }
1294  } else if (ngaps == 3) {
1295  // in this case the four corners have to be changed
1296  unsigned iubd, idir1, idir2;
1297  CaloDirection diag;
1298  CLHEP::Hep2Vector point(0., 0.);
1299  // std::cout << " Yes : 3 gaps" << std::endl;
1300  if (unbalancedDirection(gaps, iubd, idir1, idir2)) {
1301  CaloDirection ubd(gaps[iubd].first), dir1(gaps[idir1].first);
1302  CaloDirection dir2(gaps[idir2].first);
1303 
1304  // std::cout << " Avant " << std::endl << myPad << std::endl;
1305  // std::cout << ubd << " " << dir1 << " " << dir2 << std::endl;
1306  diag = CaloDirectionOperations::add2d(ubd, dir1);
1307  if (diagonalEdge(iq, diag, point))
1308  myPad.edge(diag) = point;
1309  diag = CaloDirectionOperations::add2d(ubd, dir2);
1310  if (diagonalEdge(iq, diag, point))
1311  myPad.edge(diag) = point;
1312  CaloDirection oppside = CaloDirectionOperations::oppositeSide(ubd);
1313  myPad.edge(CaloDirectionOperations::add2d(oppside, dir1)) = correspondingEdge(gaps[idir1], oppside);
1314  myPad.edge(CaloDirectionOperations::add2d(oppside, dir2)) = correspondingEdge(gaps[idir2], oppside);
1315  // std::cout << " Apres " << std::endl << myPad << std::endl;
1316  }
1317  } else if (ngaps == 4) {
1318  // std::cout << " Waouh :4 gaps" << std::endl;
1319  // std::cout << " Avant " << std::endl;
1320  // std::cout << myPad<< std::endl;
1321  CLHEP::Hep2Vector point(0., 0.);
1322  if (diagonalEdge(iq, NORTHEAST, point))
1323  myPad.edge(NORTHEAST) = point;
1324  if (diagonalEdge(iq, NORTHWEST, point))
1325  myPad.edge(NORTHWEST) = point;
1326  if (diagonalEdge(iq, SOUTHWEST, point))
1327  myPad.edge(SOUTHWEST) = point;
1328  if (diagonalEdge(iq, SOUTHEAST, point))
1329  myPad.edge(SOUTHEAST) = point;
1330  // std::cout << " Apres " << std::endl;
1331  // std::cout << myPad<< std::endl;
1332  }
1333 }
1334 
1335 void EcalHitMaker::cracksPads(std::vector<neighbour>& cracks, unsigned iq) {
1336  // std::cout << " myPad " << &myPad << std::endl;
1337  unsigned ncracks = cracks.size();
1338  CrystalPad& myPad = padsatdepth_[iq];
1339  for (unsigned ic = 0; ic < ncracks; ++ic) {
1340  // std::vector<CLHEP::Hep2Vector> mycorners;
1341  // mycorners.reserve(4);
1342  switch (cracks[ic].first) {
1343  case NORTH: {
1344  mycorners[0] = (padsatdepth_[cracks[ic].second].edge(SOUTHWEST));
1345  mycorners[1] = (padsatdepth_[cracks[ic].second].edge(SOUTHEAST));
1346  mycorners[2] = (myPad.edge(NORTHEAST));
1347  mycorners[3] = (myPad.edge(NORTHWEST));
1348  } break;
1349  case SOUTH: {
1350  mycorners[0] = (myPad.edge(SOUTHWEST));
1351  mycorners[1] = (myPad.edge(SOUTHEAST));
1352  mycorners[2] = (padsatdepth_[cracks[ic].second].edge(NORTHEAST));
1353  mycorners[3] = (padsatdepth_[cracks[ic].second].edge(NORTHWEST));
1354  } break;
1355  case EAST: {
1356  mycorners[0] = (myPad.edge(NORTHEAST));
1357  mycorners[1] = (padsatdepth_[cracks[ic].second].edge(NORTHWEST));
1358  mycorners[2] = (padsatdepth_[cracks[ic].second].edge(SOUTHWEST));
1359  mycorners[3] = (myPad.edge(SOUTHEAST));
1360  } break;
1361  case WEST: {
1362  mycorners[0] = (padsatdepth_[cracks[ic].second].edge(NORTHEAST));
1363  mycorners[1] = (myPad.edge(NORTHWEST));
1364  mycorners[2] = (myPad.edge(SOUTHWEST));
1365  mycorners[3] = (padsatdepth_[cracks[ic].second].edge(SOUTHEAST));
1366  } break;
1367  default: {
1368  }
1369  }
1370  CrystalPad crackpad(ic, mycorners);
1371  // to be tuned. A simpleconfigurable should be used
1372  crackpad.setSurvivalProbability(crackPadProbability_);
1373  crackpadsatdepth_.push_back(crackpad);
1374  }
1375  // std::cout << " Finished cracksPads " << std::endl;
1376 }
1377 
1378 bool EcalHitMaker::inside3D(const std::vector<XYZPoint>& corners, const XYZPoint& p) const {
1379  // corners and p are in the same plane
1380  // p is inside "corners" if the four crossproducts (corners[i]xcorners[i+1])
1381  // are in the same direction
1382 
1383  XYZVector crossproduct(0., 0., 0.), previouscrossproduct(0., 0., 0.);
1384 
1385  for (unsigned ip = 0; ip < 4; ++ip) {
1386  crossproduct = (corners[ip] - p).Cross(corners[(ip + 1) % 4] - p);
1387  if (ip == 0)
1388  previouscrossproduct = crossproduct;
1389  else if (crossproduct.Dot(previouscrossproduct) < 0.)
1390  return false;
1391  }
1392 
1393  return true;
1394 }
EcalHitMaker::totalL0_
double totalL0_
Definition: EcalHitMaker.h:192
DDAxes::y
CaloDirection
CaloDirection
Codes the local directions in the cell lattice.
Definition: CaloDirection.h:9
L1TDiffHarvesting_cfi.dir2
dir2
Definition: L1TDiffHarvesting_cfi.py:12
EcalHitMaker::outsideWindowEnergy_
double outsideWindowEnergy_
Definition: EcalHitMaker.h:209
EcalHitMaker::correspondingEdge
CLHEP::Hep2Vector & correspondingEdge(neighbour &myneighbour, CaloDirection dir2)
Definition: EcalHitMaker.cc:1175
CaloHitMaker::XYZPoint
math::XYZVector XYZPoint
Definition: CaloHitMaker.h:23
EcalHitMaker::L0PS1_
double L0PS1_
Definition: EcalHitMaker.h:200
EcalHitMaker::simulatePreshower_
bool simulatePreshower_
Definition: EcalHitMaker.h:256
Calorimeter::hcalProperties
const HCALProperties * hcalProperties(int onHcal) const
HCAL properties.
Definition: Calorimeter.cc:88
EcalHitMaker::buildGeometry
void buildGeometry()
Definition: EcalHitMaker.cc:789
Histos::instance
static Histos * instance()
Definition: Histos.cc:15
EcalHitMaker::preshowerCellLine
void preshowerCellLine(std::vector< CaloPoint > &cp) const
Definition: EcalHitMaker.cc:423
FSimTrack::onLayer2
int onLayer2() const
Definition: FSimTrack.h:106
CrystalWindowMap::getCrystalWindow
bool getCrystalWindow(unsigned, std::vector< unsigned > &) const
get the ordered list of the crystals around the crystal given as a first argument
Definition: CrystalWindowMap.cc:43
EcalHitMaker::gapsLifting
void gapsLifting(std::vector< neighbour > &gaps, unsigned iq)
Definition: EcalHitMaker.cc:1224
FSimTrack::onLayer1
int onLayer1() const
Definition: FSimTrack.h:101
CaloSegment::PS
Definition: CaloSegment.h:20
EcalHitMaker::L0EHGAP_
double L0EHGAP_
Definition: EcalHitMaker.h:205
PreshowerLayer1Properties.h
EcalHitMaker::setTrackParameters
void setTrackParameters(const XYZNormal &normal, double X0depthoffset, const FSimTrack &theTrack)
Definition: EcalHitMaker.cc:304
EcalHitMaker::corners
std::vector< XYZPoint > corners
Definition: EcalHitMaker.h:288
mps_update.status
status
Definition: mps_update.py:69
EcalHitMaker::central_
int central_
Definition: EcalHitMaker.h:215
Crystal
Definition: Crystal.h:20
heppy_check.dirs
dirs
Definition: heppy_check.py:26
Calorimeter::layer2Properties
const PreshowerLayer2Properties * layer2Properties(int onLayer2) const
Preshower Layer2 properties.
Definition: Calorimeter.cc:110
RawParticle::vertex
const XYZTLorentzVector & vertex() const
the vertex fourvector
Definition: RawParticle.h:320
EBDetId
Definition: EBDetId.h:17
EcalHitMaker::ecalFirstSegment_
int ecalFirstSegment_
Definition: EcalHitMaker.h:236
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
EcalHitMaker::onEcal_
int onEcal_
Definition: EcalHitMaker.h:216
gather_cfg.cout
cout
Definition: gather_cfg.py:144
FSimTrack::layer1Entrance
const RawParticle & layer1Entrance() const
The particle at Preshower Layer 1.
Definition: FSimTrack.h:137
EcalHitMaker::X0PS2EE_
double X0PS2EE_
Definition: EcalHitMaker.h:196
EcalHitMaker::ymin_
double ymin_
Definition: EcalHitMaker.h:222
EBDetId.h
EEDetId.h
DetId::Hcal
Definition: DetId.h:28
EcalHitMaker::Plane3D
ROOT::Math::Plane3D Plane3D
Definition: EcalHitMaker.h:29
EcalHitMaker::normal_
XYZNormal normal_
Definition: EcalHitMaker.h:214
EcalHitMaker::myCrystalWindowMap_
CrystalWindowMap * myCrystalWindowMap_
Definition: EcalHitMaker.h:233
EcalHitMaker::ymax_
double ymax_
Definition: EcalHitMaker.h:222
EcalHitMaker::xmax_
double xmax_
Definition: EcalHitMaker.h:222
FSimTrack.h
EcalHitMaker::pivot_
Crystal pivot_
Definition: EcalHitMaker.h:212
Transform3DR
ROOT::Math::Transform3DPJ Transform3DR
Definition: EcalHitMaker.cc:31
EcalHitMaker::diagonalEdge
bool diagonalEdge(unsigned myPad, CaloDirection dir, CLHEP::Hep2Vector &point)
Definition: EcalHitMaker.cc:1182
Validation_hcalonly_cfi.sign
sign
Definition: Validation_hcalonly_cfi.py:32
CaloHitMaker::EMSHOWER
bool EMSHOWER
Definition: CaloHitMaker.h:50
EcalHitMaker::buildSegments
void buildSegments(const std::vector< CaloPoint > &cp)
Definition: EcalHitMaker.cc:586
Variation::UP
PreshowerLayer1Properties::thickness
double thickness(double eta) const override
Definition: PreshowerLayer1Properties.cc:11
hgcal_conditions::parameters
Definition: HGCConditions.h:86
CaloHitMaker::interactionLength
double interactionLength
Definition: CaloHitMaker.h:47
FSimTrack::layer2Entrance
const RawParticle & layer2Entrance() const
The particle at Preshower Layer 2.
Definition: FSimTrack.h:140
DDAxes::x
EcalHitMaker::L0PS2EE_
double L0PS2EE_
Definition: EcalHitMaker.h:202
FSimVertex.h
EcalHitMaker::X0depthoffset_
double X0depthoffset_
Definition: EcalHitMaker.h:193
EcalHitMaker::totalX0_
double totalX0_
Definition: EcalHitMaker.h:191
CaloDirectionOperations::add2d
static CaloDirection add2d(const CaloDirection &dir1, const CaloDirection &dir2)
Definition: CaloDirectionOperations.cc:3
EcalHitMaker::doreorg_
bool doreorg_
Definition: EcalHitMaker.h:278
EcalHitMaker::validPads_
std::vector< bool > validPads_
Definition: EcalHitMaker.h:228
Calorimeter::layer1Properties
const PreshowerLayer1Properties * layer1Properties(int onLayer1) const
Preshower Layer1 properties.
Definition: Calorimeter.cc:103
createJobs.tmp
tmp
align.sh
Definition: createJobs.py:716
newFWLiteAna.found
found
Definition: newFWLiteAna.py:118
EcalHitMaker::getPads
bool getPads(double depth, bool inCm=false)
Definition: EcalHitMaker.cc:820
EcalHitMaker::ncrackpadsatdepth_
unsigned ncrackpadsatdepth_
Definition: EcalHitMaker.h:259
EcalHitMaker.h
RandomEngineAndDistribution.h
EcalBarrel
Definition: EcalSubdetector.h:10
EcalHitMaker::padsatdepth_
std::vector< CrystalPad > padsatdepth_
Definition: EcalHitMaker.h:281
EcalHitMaker::configureGeometry
void configureGeometry()
Definition: EcalHitMaker.cc:1009
EcalHitMaker::~EcalHitMaker
~EcalHitMaker() override
Definition: EcalHitMaker.cc:127
CaloDirectionOperations::oppositeDirection
static unsigned oppositeDirection(unsigned iside)
Definition: CaloDirectionOperations.cc:126
EcalHitMaker::myTrack_
const FSimTrack * myTrack_
Definition: EcalHitMaker.h:271
EcalHitMaker::addHit
bool addHit(double r, double phi, unsigned layer=0) override
Definition: EcalHitMaker.cc:184
Vector
ROOT::Math::Plane3D::Vector Vector
Definition: EcalHitMaker.cc:29
Crystal::getDetId
const DetId & getDetId() const
get the DetId
Definition: Crystal.h:46
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
EcalHitMaker::unbalancedDirection
bool unbalancedDirection(const std::vector< neighbour > &dirs, unsigned &unb, unsigned &dir1, unsigned &dir2)
Definition: EcalHitMaker.cc:1195
Calorimeter::getEcalEndcapGeometry
const EcalEndcapGeometry * getEcalEndcapGeometry() const
Definition: Calorimeter.h:52
EcalHitMaker::npadsatdepth_
unsigned npadsatdepth_
Definition: EcalHitMaker.h:260
DetId
Definition: DetId.h:17
EcalHitMaker::X0EHGAP_
double X0EHGAP_
Definition: EcalHitMaker.h:198
EcalHitMaker::X0PS1_
double X0PS1_
Definition: EcalHitMaker.h:194
alignCSCRings.s
s
Definition: alignCSCRings.py:92
CaloSegment::GAP
Definition: CaloSegment.h:20
CrystalPad
Definition: CrystalPad.h:13
CaloGeometryHelper::getClosestCell
DetId getClosestCell(const XYZPoint &point, bool ecal, bool central) const
Definition: CaloGeometryHelper.cc:58
debug
#define debug
Definition: HDRShower.cc:19
CaloSegment::PSEEGAP
Definition: CaloSegment.h:20
EcalHitMaker::radiusFactor_
double radiusFactor_
Definition: EcalHitMaker.h:248
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
EcalHitMaker::hcalCellLine
void hcalCellLine(std::vector< CaloPoint > &cp) const
Definition: EcalHitMaker.cc:468
CaloDirectionOperations::neighbourDirection
static unsigned neighbourDirection(const CaloDirection &side)
unsigned int -> Direction for the neighbours
Definition: CaloDirectionOperations.cc:19
FSimTrack::vfcalEntrance
const RawParticle & vfcalEntrance() const
The particle at VFCAL entrance.
Definition: FSimTrack.h:149
FSimTrack::onHcal
int onHcal() const
Definition: FSimTrack.h:116
EAST
Definition: CaloDirection.h:14
CalorimeterProperties::thickness
virtual double thickness(double eta) const =0
Thickness (in cm) of the homegeneous material as a function of rapidity.
CrystalPad::edge
CLHEP::Hep2Vector & edge(unsigned iside, int n)
access to the corners in direction iside; n=0,1
Definition: CrystalPad.cc:222
NORTH
Definition: CaloDirection.h:18
CaloSegment::PbWO4
Definition: CaloSegment.h:20
EcalHitMaker::sizey_
double sizey_
Definition: EcalHitMaker.h:268
Calorimetry_cff.thickness
thickness
Definition: Calorimetry_cff.py:115
EcalHitMaker::XYZNormal
math::XYZVector XYZNormal
Definition: EcalHitMaker.h:28
PVValHelper::eta
Definition: PVValidationHelpers.h:70
RawParticle::eta
double eta() const
Definition: RawParticle.h:279
FSimTrack::hcalEntrance
const RawParticle & hcalEntrance() const
The particle at HCAL entrance.
Definition: FSimTrack.h:146
CaloSegment::inSegment
Definition: CaloSegment.h:131
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
EcalHitMaker::L0HCAL_
double L0HCAL_
Definition: EcalHitMaker.h:204
SOUTHEAST
Definition: CaloDirection.h:12
CaloSegment::L0length
double L0length() const
length of the segment (in L9)
Definition: CaloSegment.h:47
EcalHitMaker::mycorners
std::vector< CLHEP::Hep2Vector > mycorners
Definition: EcalHitMaker.h:287
CaloDirectionOperations::Side
static CaloDirection Side(unsigned i)
unsigned int -> Side conversion
Definition: CaloDirectionOperations.cc:13
CaloSegment::HCAL
Definition: CaloSegment.h:20
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:150
EcalHitMaker::rearleakage_
double rearleakage_
Definition: EcalHitMaker.h:208
L1TDiffHarvesting_cfi.dir1
dir1
Definition: L1TDiffHarvesting_cfi.py:10
FSimTrack::ecalEntrance
const RawParticle & ecalEntrance() const
The particle at ECAL entrance.
Definition: FSimTrack.h:143
FSimTrack::onVFcal
int onVFcal() const
Definition: FSimTrack.h:121
b
double b
Definition: hdecay.h:118
CaloPoint::DistanceToVertex
Definition: CaloPoint.h:70
LEDCalibrationChannels.depth
depth
Definition: LEDCalibrationChannels.py:65
first
auto first
Definition: CAHitNtupletGeneratorKernelsImpl.h:112
RawParticle::Vect
XYZVector Vect() const
the momentum threevector
Definition: RawParticle.h:323
CaloHitMaker::intersect
static XYZPoint intersect(const Plane3D &p, const XYZPoint &a, const XYZPoint &b, double &t, bool segment, bool debug=false)
Definition: CaloHitMaker.cc:41
EEDetId
Definition: EEDetId.h:14
phase1PixelTopology::layer
constexpr std::array< uint8_t, layerIndexSize > layer
Definition: phase1PixelTopology.h:99
EcalHitMaker::CellsWindow_
std::vector< DetId > CellsWindow_
Definition: EcalHitMaker.h:224
EcalSubdetector.h
EcalEndcap
Definition: EcalSubdetector.h:10
CaloSegment
Definition: CaloSegment.h:15
EcalHitMaker::regionOfInterest_
std::vector< Crystal > regionOfInterest_
Definition: EcalHitMaker.h:225
EcalHitMaker::pulled
bool pulled(const XYZPoint &origin, const XYZNormal &normal, XYZPoint &fPoint) const
Definition: EcalHitMaker.cc:1066
Point
Structure Point Contains parameters of Gaussian fits to DMRs.
Definition: DMRtrends.cc:57
CaloSubdetectorGeometry.h
EcalHitMaker::X0HCAL_
double X0HCAL_
Definition: EcalHitMaker.h:199
EcalHitMaker::configuredGeometry_
bool configuredGeometry_
Definition: EcalHitMaker.h:218
bphysicsOniaDQM_cfi.vertex
vertex
Definition: bphysicsOniaDQM_cfi.py:7
EcalHitMaker::pulledPadProbability_
double pulledPadProbability_
Definition: EcalHitMaker.h:264
EcalHitMaker::bfactor_
double bfactor_
Definition: EcalHitMaker.h:254
a
double a
Definition: hdecay.h:119
EcalHitMaker::cellLine
void cellLine(std::vector< CaloPoint > &cp)
Definition: EcalHitMaker.cc:372
FSimTrack::vertex
const FSimVertex vertex() const
Origin vertex.
PVValHelper::phi
Definition: PVValidationHelpers.h:69
EcalHitMaker::segments_
std::vector< CaloSegment > segments_
Definition: EcalHitMaker.h:276
CaloGeometryHelper.h
CaloSegment::CRACK
Definition: CaloSegment.h:20
jetUpdater_cfi.sort
sort
Definition: jetUpdater_cfi.py:29
EcalHitMaker::detailedShowerTail_
bool detailedShowerTail_
Definition: EcalHitMaker.h:250
CaloGeometryHelper::buildCrystal
void buildCrystal(const DetId &id, Crystal &) const
Definition: CaloGeometryHelper.cc:157
FSimVertex::position
const math::XYZTLorentzVector & position() const
Temporary (until CMSSW moves to Mathcore) - No ! Actually very useful.
Definition: FSimVertex.h:48
EcalHitMaker::X0ECAL_
double X0ECAL_
Definition: EcalHitMaker.h:197
CrystalPad::setSurvivalProbability
void setSurvivalProbability(double val)
Definition: CrystalPad.h:64
EcalHitMaker::DetIdMap_
std::map< DetId, unsigned > DetIdMap_
Definition: EcalHitMaker.h:231
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
EcalHitMaker::reorganizePads
void reorganizePads()
Definition: EcalHitMaker.cc:1122
createfilelist.int
int
Definition: createfilelist.py:10
EcalHitMaker::hits_
std::vector< float > hits_
Definition: EcalHitMaker.h:226
EcalHitMaker::nx_
unsigned nx_
Definition: EcalHitMaker.h:221
SOUTHWEST
Definition: CaloDirection.h:13
EcalHitMaker::inside3D
bool inside3D(const std::vector< XYZPoint > &, const XYZPoint &p) const
Definition: EcalHitMaker.cc:1378
mag2
T mag2() const
The vector magnitude squared. Equivalent to vec.dot(vec)
Definition: Basic3DVectorLD.h:124
M_PI
#define M_PI
Definition: BXVectorInputProducer.cc:49
CaloGeometryHelper::getWindow
void getWindow(const DetId &pivot, int s1, int s2, std::vector< DetId > &) const
Definition: CaloGeometryHelper.cc:148
EcalHitMaker::cracksPads
void cracksPads(std::vector< neighbour > &cracks, unsigned iq)
Definition: EcalHitMaker.cc:1335
EcalHitMaker::ecalCellLine
void ecalCellLine(const XYZPoint &, const XYZPoint &, std::vector< CaloPoint > &cp)
Definition: EcalHitMaker.cc:500
EcalHitMaker::ny_
unsigned ny_
Definition: EcalHitMaker.h:221
hgcalTopologyTester_cfi.cell2
cell2
Definition: hgcalTopologyTester_cfi.py:12
CaloGeometryHelper::borderCrossing
bool borderCrossing(const DetId &, const DetId &) const
Definition: CaloGeometryHelper.cc:370
EcalHitMaker::crackPadProbability_
double crackPadProbability_
Definition: EcalHitMaker.h:266
EcalHitMaker::convertIntegerCoordinates
void convertIntegerCoordinates(double x, double y, unsigned &ix, unsigned &iy) const
Definition: EcalHitMaker.cc:1087
DetId::Detector
Detector
Definition: DetId.h:24
CaloSubdetectorGeometry::getGeometry
virtual std::shared_ptr< const CaloCellGeometry > getGeometry(const DetId &id) const
Get the cell geometry of a given detector id. Should return false if not found.
Definition: CaloSubdetectorGeometry.cc:36
PreshowerLayer2Properties.h
EcalHitMaker::ecalTotalX0
double ecalTotalX0() const
in the ECAL
Definition: EcalHitMaker.h:67
EcalHitMaker::random
const RandomEngineAndDistribution * random
Definition: EcalHitMaker.h:290
EcalHitMaker::intersections_
std::vector< CaloPoint > intersections_
Definition: EcalHitMaker.h:274
DetId::Ecal
Definition: DetId.h:27
ROOT::Math::Transform3DPJ
Definition: Transform3DPJ.h:62
RandomEngineAndDistribution::flatShoot
double flatShoot(double xmin=0.0, double xmax=1.0) const
Definition: RandomEngineAndDistribution.h:27
EcalHitMaker::hitmaphasbeencalculated_
bool hitmaphasbeencalculated_
Definition: EcalHitMaker.h:284
CaloHitMaker::spotEnergy
double spotEnergy
Definition: CaloHitMaker.h:48
HCALProperties.h
CaloHitMaker::myCalorimeter
const CaloGeometryHelper * myCalorimeter
Definition: CaloHitMaker.h:44
EcalHitMaker::crackpadsatdepth_
std::vector< CrystalPad > crackpadsatdepth_
Definition: EcalHitMaker.h:282
EcalHitMaker::addHitDepth
bool addHitDepth(double r, double phi, double depth=-1)
Definition: EcalHitMaker.cc:133
alignCSCRings.r
r
Definition: alignCSCRings.py:93
CaloHitMaker
Definition: CaloHitMaker.h:20
CaloDirectionOperations.h
DDAxes::phi
CaloHitMaker::HADSHOWER
bool HADSHOWER
Definition: CaloHitMaker.h:51
EcalPreshower
Definition: EcalSubdetector.h:10
EcalHitMaker::phisize_
unsigned phisize_
Definition: EcalHitMaker.h:240
Crystal::getFirstEdge
const XYZVector & getFirstEdge() const
Direction of the first edge.
Definition: Crystal.h:42
Transform3DPJ.h
EcalHitMaker::neighbour
std::pair< CaloDirection, unsigned > neighbour
Definition: EcalHitMaker.h:168
EcalHitMaker::prepareCrystalNumberArray
void prepareCrystalNumberArray()
Definition: EcalHitMaker.cc:1077
WEST
Definition: CaloDirection.h:15
CaloSegment::length
double length() const
length of the segment (in cm)
Definition: CaloSegment.h:43
DetId.h
CaloHitMaker::Plane3D
ROOT::Math::Plane3D Plane3D
Definition: CaloHitMaker.h:24
EcalHitMaker::maxX0_
double maxX0_
Definition: EcalHitMaker.h:207
EcalHitMaker::ncrystals_
unsigned ncrystals_
Definition: EcalHitMaker.h:219
relativeConstraints.empty
bool empty
Definition: relativeConstraints.py:46
CaloHitID
Definition: CaloHitID.h:11
CaloHitMaker::hitMap_
std::map< CaloHitID, float > hitMap_
Definition: CaloHitMaker.h:61
NORTHWEST
Definition: CaloDirection.h:17
genVertex_cff.x
x
Definition: genVertex_cff.py:12
Variation::DOWN
CaloSegment::ECALHCALGAP
Definition: CaloSegment.h:20
detailsBasic3DVector::y
float float y
Definition: extBasic3DVector.h:14
CaloGeometryHelper::magneticField
double magneticField() const
Definition: CaloGeometryHelper.h:44
EcalHitMaker::sizex_
double sizex_
Definition: EcalHitMaker.h:268
angle
T angle(T x1, T y1, T z1, T x2, T y2, T z2)
Definition: angle.h:11
NONE
Definition: TkAlStyle.cc:47
TCMET_cfi.corner
corner
Definition: TCMET_cfi.py:34
CaloGeometryHelper::move
bool move(DetId &cell, const CaloDirection &dir, bool fast=true) const
Definition: CaloGeometryHelper.cc:272
PreshowerLayer2Properties::thickness
double thickness(const double eta) const override
Definition: PreshowerLayer2Properties.cc:14
EcalHitMaker::truncatedGrid_
bool truncatedGrid_
Definition: EcalHitMaker.h:242
SOUTH
Definition: CaloDirection.h:11
M_PI_2
#define M_PI_2
Definition: CSCGattiFunction.cc:6
FSimTrack
Definition: FSimTrack.h:30
EcalHitMaker::EcalEntrance_
XYZPoint EcalEntrance_
Definition: EcalHitMaker.h:213
EcalHitMaker::X0PS2_
double X0PS2_
Definition: EcalHitMaker.h:195
EcalHitMaker::radiusCorrectionFactor_
double radiusCorrectionFactor_
Definition: EcalHitMaker.h:246
NORTHEAST
Definition: CaloDirection.h:16
EcalHitMaker::L0PS2_
double L0PS2_
Definition: EcalHitMaker.h:201
point
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5
EcalHitMaker::etasize_
unsigned etasize_
Definition: EcalHitMaker.h:239
EcalHitMaker::EcalHitMaker
EcalHitMaker(CaloGeometryHelper *calo, const XYZPoint &ecalentrance, const DetId &cell, int onEcal, unsigned size, unsigned showertype, const RandomEngineAndDistribution *engine)
Definition: EcalHitMaker.cc:33
FSimTrack::onEcal
int onEcal() const
Definition: FSimTrack.h:111
CaloGeometryHelper
Definition: CaloGeometryHelper.h:18
EcalHitMaker::currentdepth_
double currentdepth_
Definition: EcalHitMaker.h:252
submitPVValidationJobs.t
string t
Definition: submitPVValidationJobs.py:644
Point
ROOT::Math::Plane3D::Point Point
Definition: EcalHitMaker.cc:30
CaloHitMaker::moliereRadius
double moliereRadius
Definition: CaloHitMaker.h:46
EcalHitMaker::xmin_
double xmin_
Definition: EcalHitMaker.h:222
CrystalWindowMap
Definition: CrystalWindowMap.h:23
ecalpyutils::hashedIndex
int hashedIndex(int ieta, int iphi)
Definition: EcalPyUtils.cc:38
CaloSegment::X0length
double X0length() const
length of the segment (in X0)
Definition: CaloSegment.h:45
CaloHitMaker::XYZVector
math::XYZVector XYZVector
Definition: CaloHitMaker.h:22
cuy.ii
ii
Definition: cuy.py:590
EcalHitMaker::plan_
Plane3D plan_
Definition: EcalHitMaker.h:261
Calorimeter::getEcalGeometry
const CaloSubdetectorGeometry * getEcalGeometry(int subdetn) const
Definition: Calorimeter.cc:132
d1
static constexpr float d1
Definition: L1EGammaCrystalsEmulatorProducer.cc:85
geometryCSVtoXML.xx
xx
Definition: geometryCSVtoXML.py:19
EcalHitMaker::L0ECAL_
double L0ECAL_
Definition: EcalHitMaker.h:203
CrystalWindowMap.h
EcalHitMaker::myCrystalNumberArray_
std::vector< std::vector< unsigned > > myCrystalNumberArray_
Definition: EcalHitMaker.h:188
CaloSegment::inL0Segment
Definition: CaloSegment.h:109
EcalHitMaker::fastInsideCell
unsigned fastInsideCell(const CLHEP::Hep2Vector &point, double &sp, bool debug=false)
Definition: EcalHitMaker.cc:229
CaloDirectionOperations::oppositeSide
static CaloDirection oppositeSide(const CaloDirection &side)
Definition: CaloDirectionOperations.cc:86
DeadROC_duringRun.dir
dir
Definition: DeadROC_duringRun.py:23
findQualityFiles.size
size
Write out results.
Definition: findQualityFiles.py:443
EcalHitMaker::getHits
const std::map< CaloHitID, float > & getHits() override
Definition: EcalHitMaker.cc:1097
CaloSegment::inX0Segment
This class is used to determine if a point lies in the segment.
Definition: CaloSegment.h:87
RandomEngineAndDistribution
Definition: RandomEngineAndDistribution.h:18
Generated for CMSSW Reference Manual by doxygen 1.8.16