File indexing completed on 2021-02-14 13:26:35
0001 #include "DataFormats/DetId/interface/DetId.h"
0002 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0003 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0004 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0005
0006
0007
0008
0009
0010 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0011 #include "FastSimulation/CaloHitMakers/interface/EcalHitMaker.h"
0012 #include "FastSimulation/CaloGeometryTools/interface/CaloGeometryHelper.h"
0013 #include "FastSimulation/CaloGeometryTools/interface/CrystalWindowMap.h"
0014 #include "FastSimulation/CaloGeometryTools/interface/CaloDirectionOperations.h"
0015 #include "FastSimulation/Event/interface/FSimVertex.h"
0016 #include "FastSimulation/Event/interface/FSimTrack.h"
0017 #include "FastSimulation/CalorimeterProperties/interface/PreshowerLayer1Properties.h"
0018 #include "FastSimulation/CalorimeterProperties/interface/PreshowerLayer2Properties.h"
0019 #include "FastSimulation/CalorimeterProperties/interface/HCALProperties.h"
0020 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0021
0022
0023
0024 #include "FastSimulation/CaloGeometryTools/interface/Transform3DPJ.h"
0025
0026 #include <algorithm>
0027 #include <cmath>
0028
0029 typedef ROOT::Math::Plane3D::Vector Vector;
0030 typedef ROOT::Math::Plane3D::Point Point;
0031 typedef ROOT::Math::Transform3DPJ Transform3DR;
0032
0033 EcalHitMaker::EcalHitMaker(CaloGeometryHelper* theCalo,
0034 const XYZPoint& ecalentrance,
0035 const DetId& cell,
0036 int onEcal,
0037 unsigned size,
0038 unsigned showertype,
0039 const RandomEngineAndDistribution* engine)
0040 : CaloHitMaker(theCalo, DetId::Ecal, ((onEcal == 1) ? EcalBarrel : EcalEndcap), onEcal, showertype),
0041 EcalEntrance_(ecalentrance),
0042 onEcal_(onEcal),
0043 myTrack_(nullptr),
0044 random(engine) {
0045 #ifdef FAMOSDEBUG
0046 myHistos = Histos::instance();
0047 #endif
0048
0049 simulatePreshower_ = true;
0050 X0depthoffset_ = 0.;
0051 X0PS1_ = 0.;
0052 X0PS2_ = 0.;
0053 X0PS2EE_ = 0.;
0054 X0ECAL_ = 0.;
0055 X0EHGAP_ = 0.;
0056 X0HCAL_ = 0.;
0057 L0PS1_ = 0.;
0058 L0PS2_ = 0.;
0059 L0PS2EE_ = 0.;
0060 L0ECAL_ = 0.;
0061 L0EHGAP_ = 0.;
0062 L0HCAL_ = 0.;
0063 maxX0_ = 0.;
0064 totalX0_ = 0;
0065 totalL0_ = 0.;
0066 pulledPadProbability_ = 1.;
0067 outsideWindowEnergy_ = 0.;
0068 rearleakage_ = 0.;
0069 bfactor_ = 1.;
0070 ncrystals_ = 0;
0071
0072 doreorg_ = !showertype;
0073
0074 hitmaphasbeencalculated_ = false;
0075
0076 if (onEcal)
0077 myCalorimeter->buildCrystal(cell, pivot_);
0078 else
0079 pivot_ = Crystal();
0080 central_ = onEcal == 1;
0081 ecalFirstSegment_ = -1;
0082
0083 myCrystalWindowMap_ = nullptr;
0084
0085
0086 if (!onEcal)
0087 return;
0088
0089
0090 etasize_ = size;
0091 phisize_ = size;
0092
0093
0094
0095 myCalorimeter->getWindow(pivot_.getDetId(), size, size, CellsWindow_);
0096
0097 buildGeometry();
0098
0099
0100 truncatedGrid_ = CellsWindow_.size() != (etasize_ * phisize_);
0101
0102
0103 mycorners.resize(4);
0104 corners.resize(4);
0105
0106 #ifdef DEBUGGW
0107 myHistos->fill("h10", EcalEntrance_.eta(), CellsWindow_.size());
0108 if (onEcal == 2) {
0109 myHistos->fill("h20", EcalEntrance_.perp(), CellsWindow_.size());
0110 if (EcalEntrance_.perp() > 70 && EcalEntrance_.perp() < 80 && CellsWindow_.size() < 35) {
0111 std::cout << " Truncated grid " << CellsWindow_.size() << " " << EcalEntrance_.perp() << std::endl;
0112 std::cout << " Pivot "
0113 << myCalorimeter->getEcalEndcapGeometry()->getGeometry(pivot_.getDetId())->getPosition().perp();
0114 std::cout << EEDetId(pivot_.getDetId()) << std::endl;
0115
0116 std::cout << " Test getClosestCell " << EcalEntrance_ << std::endl;
0117 DetId testcell = myCalorimeter->getClosestCell(EcalEntrance_, true, false);
0118 std::cout << " Result " << EEDetId(testcell) << std::endl;
0119 std::cout << " Position " << myCalorimeter->getEcalEndcapGeometry()->getGeometry(testcell)->getPosition()
0120 << std::endl;
0121 }
0122 }
0123
0124 #endif
0125 }
0126
0127 EcalHitMaker::~EcalHitMaker() {
0128 if (myCrystalWindowMap_ != nullptr) {
0129 delete myCrystalWindowMap_;
0130 }
0131 }
0132
0133 bool EcalHitMaker::addHitDepth(double r, double phi, double depth) {
0134
0135
0136 depth += X0depthoffset_;
0137 double sp(1.);
0138 r *= radiusFactor_;
0139 CLHEP::Hep2Vector point(r * std::cos(phi), r * std::sin(phi));
0140
0141 unsigned xtal = fastInsideCell(point, sp);
0142
0143
0144
0145
0146
0147
0148
0149 if (xtal < 1000) {
0150
0151
0152 if (regionOfInterest_[xtal].getX0Back() > depth) {
0153 hits_[xtal] += spotEnergy;
0154
0155 return true;
0156 } else {
0157 rearleakage_ += spotEnergy;
0158 }
0159 } else {
0160
0161
0162
0163
0164 }
0165
0166 outsideWindowEnergy_ += spotEnergy;
0167 return false;
0168 }
0169
0170
0171
0172
0173
0174
0175
0176
0177
0178
0179
0180
0181
0182
0183
0184 bool EcalHitMaker::addHit(double r, double phi, unsigned layer) {
0185
0186
0187 double sp(1.);
0188
0189 r *= radiusFactor_;
0190 CLHEP::Hep2Vector point(r * std::cos(phi), r * std::sin(phi));
0191
0192
0193 unsigned xtal = fastInsideCell(point, sp);
0194
0195 if (xtal < 1000) {
0196 if (sp == 1.)
0197 hits_[xtal] += spotEnergy;
0198 else
0199 hits_[xtal] += (random->flatShoot() < sp) * spotEnergy;
0200 return true;
0201 }
0202
0203 outsideWindowEnergy_ += spotEnergy;
0204
0205
0206
0207
0208
0209
0210
0211 return false;
0212 }
0213
0214
0215
0216
0217
0218
0219
0220
0221
0222
0223
0224
0225
0226
0227
0228
0229 unsigned EcalHitMaker::fastInsideCell(const CLHEP::Hep2Vector& point, double& sp, bool debug) {
0230
0231 bool found = false;
0232 unsigned niter = 0;
0233
0234 unsigned d1, d2;
0235 convertIntegerCoordinates(point.x(), point.y(), d1, d2);
0236
0237 if (d1 >= nx_ || d2 >= ny_) {
0238
0239 return 9999;
0240 }
0241 unsigned cell = myCrystalNumberArray_[d1][d2];
0242
0243
0244 if (validPads_[cell] && padsatdepth_[cell].inside(point)) {
0245
0246 sp = padsatdepth_[cell].survivalProbability();
0247 return cell;
0248 }
0249
0250
0251 bool status(true);
0252 const std::vector<unsigned>& localCellVector(myCrystalWindowMap_->getCrystalWindow(cell, status));
0253 if (status) {
0254 unsigned size = localCellVector.size();
0255
0256
0257
0258
0259
0260
0261
0262
0263 for (unsigned ic = 0; ic < 8 && ic < size; ++ic) {
0264 unsigned iq = localCellVector[ic];
0265
0266
0267
0268 if (validPads_[iq] && padsatdepth_[iq].inside(point)) {
0269
0270
0271 sp = padsatdepth_[iq].survivalProbability();
0272
0273
0274 return iq;
0275 }
0276
0277
0278 ++niter;
0279 }
0280 }
0281 if (debug)
0282 std::cout << " not found in a quad, let's check the " << ncrackpadsatdepth_ << " cracks " << std::endl;
0283
0284
0285
0286 unsigned iquad = 0;
0287 unsigned iquadinside = 999;
0288 while (iquad < ncrackpadsatdepth_ && !found) {
0289
0290 if (crackpadsatdepth_[iquad].inside(point)) {
0291 iquadinside = iquad;
0292 found = true;
0293 sp = crackpadsatdepth_[iquad].survivalProbability();
0294 }
0295 ++iquad;
0296 ++niter;
0297 }
0298
0299 if (!found && debug)
0300 std::cout << " Not found in the cracks " << std::endl;
0301 return (found) ? crackpadsatdepth_[iquadinside].getNumber() : 9999;
0302 }
0303
0304 void EcalHitMaker::setTrackParameters(const XYZNormal& normal, double X0depthoffset, const FSimTrack& theTrack) {
0305
0306
0307 intersections_.clear();
0308
0309 intersections_.reserve(50);
0310 myTrack_ = &theTrack;
0311 normal_ = normal.Unit();
0312 X0depthoffset_ = X0depthoffset;
0313 cellLine(intersections_);
0314 buildSegments(intersections_);
0315
0316
0317
0318
0319
0320
0321
0322 if (EMSHOWER && onEcal_ && ecalTotalX0() > 0.) {
0323
0324 for (unsigned ic = 0; ic < ncrystals_; ++ic) {
0325 for (unsigned idir = 0; idir < 4; ++idir) {
0326 XYZVector norm = regionOfInterest_[ic].exitingNormal(CaloDirectionOperations::Side(idir));
0327 regionOfInterest_[ic].crystalNeighbour(idir).setToBeGlued((norm.Dot(normal_) < 0.));
0328 }
0329
0330
0331 if (EMSHOWER) {
0332 XYZVector dir = regionOfInterest_[ic].getBackCenter() - segments_[ecalFirstSegment_].entrance();
0333 double dist = dir.Dot(normal_);
0334 double absciss = dist + segments_[ecalFirstSegment_].sEntrance();
0335 std::vector<CaloSegment>::const_iterator segiterator;
0336
0337
0338
0339
0340 segiterator = find_if(segments_.begin(), segments_.end(), CaloSegment::inSegment(absciss));
0341
0342 if (segiterator == segments_.end()) {
0343
0344
0345 regionOfInterest_[ic].setX0Back(9999);
0346 } else {
0347 DetId::Detector det(segiterator->whichDetector());
0348 if (det != DetId::Ecal) {
0349 regionOfInterest_[ic].setX0Back(9999);
0350 } else {
0351 double x0 = segiterator->x0FromCm(dist);
0352 if (x0 < maxX0_)
0353 maxX0_ = x0;
0354 regionOfInterest_[ic].setX0Back(x0);
0355 }
0356 }
0357
0358
0359
0360
0361
0362
0363
0364
0365 }
0366 }
0367
0368 }
0369
0370 }
0371
0372 void EcalHitMaker::cellLine(std::vector<CaloPoint>& cp) {
0373 cp.clear();
0374
0375
0376 if (!central_ && onEcal_ && simulatePreshower_)
0377 preshowerCellLine(cp);
0378 if (onEcal_)
0379 ecalCellLine(EcalEntrance_, EcalEntrance_ + normal_, cp);
0380
0381
0382 XYZPoint vertex(myTrack_->vertex().position().Vect());
0383
0384
0385 XYZVector dir(0., 0., 0.);
0386 if (myTrack_->onLayer1()) {
0387 vertex = (myTrack_->layer1Entrance().vertex()).Vect();
0388 dir = myTrack_->layer1Entrance().Vect().Unit();
0389 } else if (myTrack_->onLayer2()) {
0390 vertex = (myTrack_->layer2Entrance().vertex()).Vect();
0391 dir = myTrack_->layer2Entrance().Vect().Unit();
0392 } else if (myTrack_->onEcal()) {
0393 vertex = (myTrack_->ecalEntrance().vertex()).Vect();
0394 dir = myTrack_->ecalEntrance().Vect().Unit();
0395 } else if (myTrack_->onHcal()) {
0396 vertex = (myTrack_->hcalEntrance().vertex()).Vect();
0397 dir = myTrack_->hcalEntrance().Vect().Unit();
0398 } else if (myTrack_->onVFcal() == 2) {
0399 vertex = (myTrack_->vfcalEntrance().vertex()).Vect();
0400 dir = myTrack_->vfcalEntrance().Vect().Unit();
0401 } else {
0402 std::cout << " Problem with the grid " << std::endl;
0403 }
0404
0405
0406 vertex -= 5. * dir;
0407 CaloPoint::DistanceToVertex myDistance(vertex);
0408 sort(cp.begin(), cp.end(), myDistance);
0409
0410
0411
0412 hcalCellLine(cp);
0413
0414
0415
0416
0417
0418
0419
0420
0421 }
0422
0423 void EcalHitMaker::preshowerCellLine(std::vector<CaloPoint>& cp) const {
0424
0425
0426
0427
0428
0429 if (myTrack_->onLayer1()) {
0430 XYZPoint point1 = (myTrack_->layer1Entrance().vertex()).Vect();
0431 double phys_eta = myTrack_->layer1Entrance().eta();
0432 double cmthickness = myCalorimeter->layer1Properties(1)->thickness(phys_eta);
0433
0434 if (cmthickness > 0) {
0435 XYZVector dir = myTrack_->layer1Entrance().Vect().Unit();
0436 XYZPoint point2 = point1 + dir * cmthickness;
0437
0438 CaloPoint cp1(DetId::Ecal, EcalPreshower, 1, point1);
0439 CaloPoint cp2(DetId::Ecal, EcalPreshower, 1, point2);
0440 cp.push_back(cp1);
0441 cp.push_back(cp2);
0442 } else {
0443
0444 }
0445 }
0446
0447
0448 if (myTrack_->onLayer2()) {
0449 XYZPoint point1 = (myTrack_->layer2Entrance().vertex()).Vect();
0450 double phys_eta = myTrack_->layer2Entrance().eta();
0451 double cmthickness = myCalorimeter->layer2Properties(1)->thickness(phys_eta);
0452 if (cmthickness > 0) {
0453 XYZVector dir = myTrack_->layer2Entrance().Vect().Unit();
0454 XYZPoint point2 = point1 + dir * cmthickness;
0455
0456 CaloPoint cp1(DetId::Ecal, EcalPreshower, 2, point1);
0457 CaloPoint cp2(DetId::Ecal, EcalPreshower, 2, point2);
0458
0459 cp.push_back(cp1);
0460 cp.push_back(cp2);
0461 } else {
0462
0463 }
0464 }
0465
0466 }
0467
0468 void EcalHitMaker::hcalCellLine(std::vector<CaloPoint>& cp) const {
0469
0470
0471 int onHcal = myTrack_->onHcal();
0472
0473 if (onHcal <= 2 && onHcal > 0) {
0474 XYZPoint point1 = (myTrack_->hcalEntrance().vertex()).Vect();
0475
0476 double eta = point1.eta();
0477
0478 double thickness = myCalorimeter->hcalProperties(onHcal)->thickness(eta);
0479 cp.push_back(CaloPoint(DetId::Hcal, point1));
0480 XYZVector dir = myTrack_->hcalEntrance().Vect().Unit();
0481 XYZPoint point2 = point1 + dir * thickness;
0482
0483 cp.push_back(CaloPoint(DetId::Hcal, point2));
0484 }
0485 int onVFcal = myTrack_->onVFcal();
0486 if (onVFcal == 2) {
0487 XYZPoint point1 = (myTrack_->vfcalEntrance().vertex()).Vect();
0488 double eta = point1.eta();
0489
0490 double thickness = myCalorimeter->hcalProperties(3)->thickness(eta);
0491 cp.push_back(CaloPoint(DetId::Hcal, point1));
0492 XYZVector dir = myTrack_->vfcalEntrance().Vect().Unit();
0493 if (thickness > 0) {
0494 XYZPoint point2 = point1 + dir * thickness;
0495 cp.push_back(CaloPoint(DetId::Hcal, point2));
0496 }
0497 }
0498 }
0499
0500 void EcalHitMaker::ecalCellLine(const XYZPoint& a, const XYZPoint& b, std::vector<CaloPoint>& cp) {
0501
0502
0503 unsigned ic = 0;
0504 double t;
0505 XYZPoint xp;
0506 DetId c_entrance, c_exit;
0507 bool entrancefound(false), exitfound(false);
0508
0509
0510
0511
0512 double angle = std::acos(normal_.Dot(regionOfInterest_[0].getAxis().Unit()));
0513
0514
0515 double backdistance = std::sqrt(regionOfInterest_[0].getAxis().mag2()) * std::tan(angle);
0516
0517
0518
0519 unsigned ncrystals = (unsigned)(backdistance * 0.45);
0520 unsigned highlim = (ncrystals + 4);
0521 highlim *= highlim;
0522 if (highlim > ncrystals_)
0523 highlim = ncrystals_;
0524
0525
0526
0527 while (ic < ncrystals_ && (ic < highlim || !exitfound)) {
0528
0529
0530 {
0531 const Plane3D& plan = regionOfInterest_[ic].getFrontPlane();
0532
0533
0534 xp = intersect(plan, a, b, t, false);
0535 regionOfInterest_[ic].getFrontSide(corners);
0536
0537
0538 if (inside3D(corners, xp)) {
0539 cp.push_back(CaloPoint(regionOfInterest_[ic].getDetId(), UP, xp));
0540 entrancefound = true;
0541 c_entrance = regionOfInterest_[ic].getDetId();
0542
0543 }
0544 }
0545
0546
0547
0548 {
0549 const Plane3D& plan = regionOfInterest_[ic].getBackPlane();
0550
0551
0552 xp = intersect(plan, a, b, t, false);
0553 regionOfInterest_[ic].getBackSide(corners);
0554
0555
0556 if (inside3D(corners, xp)) {
0557 cp.push_back(CaloPoint(regionOfInterest_[ic].getDetId(), DOWN, xp));
0558 exitfound = true;
0559 c_exit = regionOfInterest_[ic].getDetId();
0560
0561
0562 }
0563 }
0564
0565 if (entrancefound && exitfound && c_entrance == c_exit)
0566 return;
0567
0568 for (unsigned iside = 0; iside < 4; ++iside) {
0569 const Plane3D& plan = regionOfInterest_[ic].getLateralPlane(iside);
0570 xp = intersect(plan, a, b, t, false);
0571
0572
0573 regionOfInterest_[ic].getLateralSide(iside, corners);
0574
0575
0576 if (inside3D(corners, xp)) {
0577 cp.push_back(CaloPoint(regionOfInterest_[ic].getDetId(), CaloDirectionOperations::Side(iside), xp));
0578
0579 }
0580 }
0581
0582 ++ic;
0583 }
0584 }
0585
0586 void EcalHitMaker::buildSegments(const std::vector<CaloPoint>& cp) {
0587
0588
0589 unsigned size = cp.size();
0590 if (size % 2 != 0) {
0591
0592 return;
0593 }
0594
0595 unsigned nsegments = size / 2;
0596 segments_.reserve(nsegments);
0597 if (size == 0)
0598 return;
0599
0600 double s = 0.;
0601 double sX0 = 0.;
0602 double sL0 = 0.;
0603
0604 unsigned ncrossedxtals = 0;
0605 unsigned is = 0;
0606 while (is < nsegments) {
0607 if (cp[2 * is].getDetId() != cp[2 * is + 1].getDetId() && cp[2 * is].whichDetector() != DetId::Hcal &&
0608 cp[2 * is + 1].whichDetector() != DetId::Hcal) {
0609
0610
0611
0612
0613 ++is;
0614 continue;
0615 }
0616
0617
0618
0619
0620 if (cp[2 * is].whichDetector() == DetId::Ecal && cp[2 * is].whichSubDetector() == EcalPreshower &&
0621 cp[2 * is].whichLayer() == 1) {
0622 if (cp[2 * is + 1].whichDetector() == DetId::Ecal && cp[2 * is + 1].whichSubDetector() == EcalPreshower &&
0623 cp[2 * is + 1].whichLayer() == 1) {
0624 CaloSegment preshsegment(cp[2 * is], cp[2 * is + 1], s, sX0, sL0, CaloSegment::PS, myCalorimeter);
0625 segments_.push_back(preshsegment);
0626
0627 s += preshsegment.length();
0628 sX0 += preshsegment.X0length();
0629 sL0 += preshsegment.L0length();
0630 X0PS1_ += preshsegment.X0length();
0631 L0PS1_ += preshsegment.L0length();
0632 } else {
0633 std::cout << " Strange segment between Preshower1 and " << cp[2 * is + 1].whichDetector();
0634 std::cout << std::endl;
0635 }
0636 ++is;
0637 continue;
0638 }
0639
0640
0641
0642 if (cp[2 * is].whichDetector() == DetId::Ecal && cp[2 * is].whichSubDetector() == EcalPreshower &&
0643 cp[2 * is].whichLayer() == 2) {
0644 if (cp[2 * is + 1].whichDetector() == DetId::Ecal && cp[2 * is + 1].whichSubDetector() == EcalPreshower &&
0645 cp[2 * is + 1].whichLayer() == 2) {
0646 CaloSegment preshsegment(cp[2 * is], cp[2 * is + 1], s, sX0, sL0, CaloSegment::PS, myCalorimeter);
0647 segments_.push_back(preshsegment);
0648
0649 s += preshsegment.length();
0650 sX0 += preshsegment.X0length();
0651 sL0 += preshsegment.L0length();
0652 X0PS2_ += preshsegment.X0length();
0653 L0PS2_ += preshsegment.L0length();
0654
0655
0656 if (is < (nsegments - 1) && cp[2 * is + 2].whichDetector() == DetId::Ecal &&
0657 cp[2 * is + 2].whichSubDetector() == EcalEndcap) {
0658 CaloSegment gapsef(cp[2 * is + 1], cp[2 * is + 2], s, sX0, sL0, CaloSegment::PSEEGAP, myCalorimeter);
0659 segments_.push_back(gapsef);
0660 s += gapsef.length();
0661 sX0 += gapsef.X0length();
0662 sL0 += gapsef.L0length();
0663 X0PS2EE_ += gapsef.X0length();
0664 L0PS2EE_ += gapsef.L0length();
0665
0666 }
0667 } else {
0668 std::cout << " Strange segment between Preshower2 and " << cp[2 * is + 1].whichDetector();
0669 std::cout << std::endl;
0670 }
0671 ++is;
0672 continue;
0673 }
0674
0675
0676
0677 if (cp[2 * is].whichDetector() == DetId::Ecal &&
0678 (cp[2 * is].whichSubDetector() == EcalBarrel || cp[2 * is].whichSubDetector() == EcalEndcap)) {
0679 if (cp[2 * is + 1].whichDetector() == DetId::Ecal &&
0680 (cp[2 * is + 1].whichSubDetector() == EcalBarrel || cp[2 * is + 1].whichSubDetector() == EcalEndcap)) {
0681 DetId cell2 = cp[2 * is + 1].getDetId();
0682
0683 if (ecalFirstSegment_ < 0)
0684 ecalFirstSegment_ = segments_.size();
0685
0686
0687 if (cp[2 * is].getDetId() == cell2) {
0688 CaloSegment segment(cp[2 * is], cp[2 * is + 1], s, sX0, sL0, CaloSegment::PbWO4, myCalorimeter);
0689 segments_.push_back(segment);
0690
0691 s += segment.length();
0692 sX0 += segment.X0length();
0693 sL0 += segment.L0length();
0694 X0ECAL_ += segment.X0length();
0695 L0ECAL_ += segment.L0length();
0696 ++ncrossedxtals;
0697 ++is;
0698 } else {
0699 std::cout << " One more bug in the segment " << std::endl;
0700 ++is;
0701 }
0702
0703 if (is > 0 && is < nsegments) {
0704 DetId cell3 = cp[2 * is].getDetId();
0705 if (cp[2 * is].whichDetector() != DetId::Hcal) {
0706
0707 bool bordercrossing = myCalorimeter->borderCrossing(cell2, cell3);
0708 CaloSegment cracksegment(cp[2 * is - 1],
0709 cp[2 * is],
0710 s,
0711 sX0,
0712 sL0,
0713 (bordercrossing) ? CaloSegment::CRACK : CaloSegment::GAP,
0714 myCalorimeter);
0715 segments_.push_back(cracksegment);
0716 s += cracksegment.length();
0717 sX0 += cracksegment.X0length();
0718 sL0 += cracksegment.L0length();
0719 X0ECAL_ += cracksegment.X0length();
0720 L0ECAL_ += cracksegment.L0length();
0721
0722 } else {
0723
0724
0725 CaloSegment cracksegment(cp[2 * is - 1], cp[2 * is], s, sX0, sL0, CaloSegment::ECALHCALGAP, myCalorimeter);
0726 segments_.push_back(cracksegment);
0727 s += cracksegment.length();
0728 sX0 += cracksegment.X0length();
0729 sL0 += cracksegment.L0length();
0730 X0EHGAP_ += cracksegment.X0length();
0731 L0EHGAP_ += cracksegment.L0length();
0732 }
0733 }
0734 continue;
0735 } else {
0736 std::cout << " Strange segment between " << cp[2 * is].whichDetector();
0737 std::cout << " and " << cp[2 * is + 1].whichDetector() << std::endl;
0738 ++is;
0739 continue;
0740 }
0741 }
0742
0743
0744 if (cp[2 * is].whichDetector() == DetId::Hcal && cp[2 * is + 1].whichDetector() == DetId::Hcal) {
0745 CaloSegment segment(cp[2 * is], cp[2 * is + 1], s, sX0, sL0, CaloSegment::HCAL, myCalorimeter);
0746 segments_.push_back(segment);
0747 s += segment.length();
0748 sX0 += segment.X0length();
0749 sL0 += segment.L0length();
0750 X0HCAL_ += segment.X0length();
0751 L0HCAL_ += segment.L0length();
0752
0753 ++is;
0754 }
0755 }
0756
0757
0758
0759
0760
0761 totalX0_ = X0PS1_ + X0PS2_ + X0PS2EE_ + X0ECAL_ + X0EHGAP_ + X0HCAL_;
0762 totalL0_ = L0PS1_ + L0PS2_ + L0PS2EE_ + L0ECAL_ + L0EHGAP_ + L0HCAL_;
0763
0764
0765 #ifdef DEBUGCELLLINE
0766 myHistos->fill("h200", fabs(EcalEntrance_.eta()), X0ECAL_);
0767 myHistos->fill("h210", EcalEntrance_.phi(), X0ECAL_);
0768 if (X0ECAL_ < 20)
0769 myHistos->fill("h212", EcalEntrance_.phi(), X0ECAL_);
0770
0771
0772
0773
0774
0775
0776
0777 myHistos->fillByNumber("h30", ncrossedxtals, EcalEntrance_.eta(), X0ECAL_);
0778
0779 double zvertex = myTrack_->vertex().position().z();
0780
0781 myHistos->fill("h310", EcalEntrance_.eta(), X0ECAL_);
0782 if (X0ECAL_ < 22)
0783 myHistos->fill("h410", EcalEntrance_.phi());
0784 myHistos->fill("h400", zvertex, X0ECAL_);
0785 #endif
0786
0787 }
0788
0789 void EcalHitMaker::buildGeometry() {
0790 configuredGeometry_ = false;
0791 ncrystals_ = CellsWindow_.size();
0792
0793 padsatdepth_.resize(ncrystals_);
0794
0795
0796 ny_ = phisize_;
0797 nx_ = ncrystals_ / ny_;
0798 std::vector<unsigned> empty;
0799 empty.resize(ny_, 0);
0800 myCrystalNumberArray_.reserve((unsigned)nx_);
0801 for (unsigned inx = 0; inx < (unsigned)nx_; ++inx) {
0802 myCrystalNumberArray_.push_back(empty);
0803 }
0804
0805 hits_.resize(ncrystals_, 0.);
0806 regionOfInterest_.clear();
0807 regionOfInterest_.resize(ncrystals_);
0808 validPads_.resize(ncrystals_);
0809 for (unsigned ic = 0; ic < ncrystals_; ++ic) {
0810 myCalorimeter->buildCrystal(CellsWindow_[ic], regionOfInterest_[ic]);
0811 regionOfInterest_[ic].setNumber(ic);
0812 DetIdMap_.insert(std::pair<DetId, unsigned>(CellsWindow_[ic], ic));
0813 }
0814
0815
0816 myCrystalWindowMap_ = new CrystalWindowMap(myCalorimeter, regionOfInterest_);
0817 }
0818
0819
0820 bool EcalHitMaker::getPads(double depth, bool inCm) {
0821
0822
0823
0824
0825 if (EMSHOWER && !configuredGeometry_)
0826 configureGeometry();
0827
0828 radiusFactor_ = (EMSHOWER) ? moliereRadius * radiusCorrectionFactor_ : interactionLength;
0829 detailedShowerTail_ = false;
0830 if (EMSHOWER)
0831 currentdepth_ = depth + X0depthoffset_;
0832 else
0833 currentdepth_ = depth;
0834
0835
0836
0837
0838
0839
0840
0841
0842
0843 ncrackpadsatdepth_ = 0;
0844
0845 xmin_ = ymin_ = 999;
0846 xmax_ = ymax_ = -999;
0847 double locxmin, locxmax, locymin, locymax;
0848
0849
0850 std::vector<CaloSegment>::const_iterator segiterator;
0851
0852
0853 if (inCm)
0854 {
0855 segiterator = find_if(segments_.begin(), segments_.end(), CaloSegment::inSegment(currentdepth_));
0856 } else {
0857
0858 if (EMSHOWER)
0859 segiterator = find_if(segments_.begin(), segments_.end(), CaloSegment::inX0Segment(currentdepth_));
0860
0861
0862 if (HADSHOWER)
0863 segiterator = find_if(segments_.begin(), segments_.end(), CaloSegment::inL0Segment(currentdepth_));
0864 }
0865 if (segiterator == segments_.end()) {
0866 std::cout << " FamosGrid: Could not go at such depth " << depth << std::endl;
0867 std::cout << " EMSHOWER " << EMSHOWER << std::endl;
0868 std::cout << " Track " << *myTrack_ << std::endl;
0869 std::cout << " Segments " << segments_.size() << std::endl;
0870 for (unsigned ii = 0; ii < segments_.size(); ++ii) {
0871 std::cout << segments_[ii] << std::endl;
0872 }
0873
0874 return false;
0875 }
0876
0877
0878 if (segiterator->whichDetector() != DetId::Ecal) {
0879 std::cout << " In " << segiterator->whichDetector() << std::endl;
0880
0881 return false;
0882 }
0883
0884
0885
0886
0887 XYZPoint origin;
0888 if (inCm) {
0889 origin = segiterator->positionAtDepthincm(currentdepth_);
0890 } else {
0891 if (EMSHOWER)
0892 origin = segiterator->positionAtDepthinX0(currentdepth_);
0893 if (HADSHOWER)
0894 origin = segiterator->positionAtDepthinL0(currentdepth_);
0895 }
0896
0897 XYZVector newaxis = pivot_.getFirstEdge().Cross(normal_);
0898
0899
0900
0901
0902
0903 plan_ = Plane3D((Vector)normal_, (Point)origin);
0904
0905 unsigned nquads = 0;
0906 double sign = (central_) ? -1. : 1.;
0907 Transform3DR trans((Point)origin,
0908 (Point)(origin + normal_),
0909 (Point)(origin + newaxis),
0910 Point(0, 0, 0),
0911 Point(0., 0., sign),
0912 Point(0., 1., 0.));
0913 for (unsigned ic = 0; ic < ncrystals_; ++ic) {
0914
0915 XYZPoint a, b;
0916
0917
0918
0919
0920
0921 double dummyt;
0922 bool hasbeenpulled = false;
0923 bool behindback = false;
0924 for (unsigned il = 0; il < 4; ++il) {
0925
0926 regionOfInterest_[ic].getLateralEdges(il, a, b);
0927
0928
0929 XYZPoint aprime = a;
0930 if (pulled(origin, normal_, a)) {
0931 b = aprime;
0932 hasbeenpulled = true;
0933 }
0934
0935
0936
0937
0938 XYZPoint xx = (EMSHOWER) ? intersect(plan_, a, b, dummyt, false) : intersect(plan_, a, b, dummyt, true);
0939
0940 if (dummyt > 1)
0941 behindback = true;
0942
0943
0944 if (xx.mag2() != 0) {
0945 corners[il] = xx;
0946 }
0947 }
0948
0949 if (behindback && EMSHOWER)
0950 detailedShowerTail_ = true;
0951
0952 if (corners.size() == 4) {
0953 padsatdepth_[ic] = CrystalPad(ic, corners, trans, bfactor_, !central_);
0954
0955 if (hasbeenpulled)
0956 padsatdepth_[ic].setSurvivalProbability(pulledPadProbability_);
0957 validPads_[ic] = true;
0958 ++nquads;
0959
0960
0961
0962
0963
0964
0965 } else {
0966 padsatdepth_[ic] = CrystalPad();
0967 validPads_[ic] = false;
0968 }
0969 }
0970
0971 if (doreorg_)
0972 reorganizePads();
0973
0974 npadsatdepth_ = nquads;
0975
0976
0977
0978
0979 for (unsigned ic = 0; ic < ncrystals_; ++ic) {
0980 if (!validPads_[ic])
0981 continue;
0982
0983 if (EMSHOWER)
0984 padsatdepth_[ic].resetCorners();
0985
0986 padsatdepth_[ic].extrems(locxmin, locxmax, locymin, locymax);
0987 if (locxmin < xmin_)
0988 xmin_ = locxmin;
0989 if (locymin < ymin_)
0990 ymin_ = locymin;
0991 if (locxmax > xmax_)
0992 xmax_ = locxmax;
0993 if (locymax > ymax_)
0994 ymax_ = locymax;
0995 }
0996
0997 sizex_ = (xmax_ - xmin_) / nx_;
0998 sizey_ = (ymax_ - ymin_) / ny_;
0999
1000
1001 prepareCrystalNumberArray();
1002
1003
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
1013 for (unsigned idir = 0; idir < 8; ++idir) {
1014 unsigned oppdir = CaloDirectionOperations::oppositeDirection(idir);
1015
1016
1017 if (regionOfInterest_[ic].crystalNeighbour(idir).status() >= 0) {
1018
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
1027 regionOfInterest_[ic].crystalNeighbour(idir).setStatus(-1);
1028 continue;
1029 }
1030
1031
1032 std::map<DetId, unsigned>::const_iterator niter(DetIdMap_.find(newcell));
1033 if (niter == DetIdMap_.end()) {
1034
1035 regionOfInterest_[ic].crystalNeighbour(idir).setStatus(-1);
1036 continue;
1037 }
1038
1039
1040 regionOfInterest_[ic].crystalNeighbour(idir).setNumber(niter->second);
1041
1042
1043 regionOfInterest_[niter->second].crystalNeighbour(oppdir).setNumber(ic);
1044
1045 if (myCalorimeter->borderCrossing(oldcell, newcell)) {
1046 regionOfInterest_[ic].crystalNeighbour(idir).setStatus(1);
1047 regionOfInterest_[niter->second].crystalNeighbour(oppdir).setStatus(1);
1048
1049 } else {
1050 regionOfInterest_[ic].crystalNeighbour(idir).setStatus(0);
1051 regionOfInterest_[niter->second].crystalNeighbour(oppdir).setStatus(0);
1052
1053 }
1054 }
1055 }
1056
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
1066 bool EcalHitMaker::pulled(const XYZPoint& origin, const XYZNormal& normal, XYZPoint& fPoint) const {
1067
1068 double dotproduct = normal.Dot(fPoint - origin);
1069 if (dotproduct <= 0.)
1070 return false;
1071
1072
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
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;
1107
1108 if (onEcal_ == 1) {
1109 CaloHitID current_id(EBDetId(regionOfInterest_[ic].getDetId().rawId()).hashedIndex(), tof, 0);
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);
1113 hitMap_.insert(std::pair<CaloHitID, float>(current_id, hits_[ic]));
1114 }
1115 }
1116 hitmaphasbeencalculated_ = true;
1117 return hitMap_;
1118 }
1119
1120
1121
1122 void EcalHitMaker::reorganizePads() {
1123
1124
1125 crackpadsatdepth_.clear();
1126 crackpadsatdepth_.reserve(etasize_ * phisize_);
1127 ncrackpadsatdepth_ = 0;
1128 std::vector<neighbour> gaps;
1129 std::vector<std::vector<neighbour> > cracks;
1130
1131 cracks.resize(ncrystals_);
1132
1133 for (unsigned iq = 0; iq < ncrystals_; ++iq) {
1134 if (!validPads_[iq])
1135 continue;
1136
1137 gaps.clear();
1138
1139
1140 for (unsigned iside = 0; iside < 4; ++iside) {
1141
1142 CaloDirection thisside = CaloDirectionOperations::Side(iside);
1143 if (regionOfInterest_[iq].crystalNeighbour(iside).toBeGlued()) {
1144
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
1153 if (neighbourstatus == 1) {
1154
1155 cracks[iq].push_back(neighbour(thisside, neighbourNumber));
1156 }
1157 else {
1158 gaps.push_back(neighbour(thisside, neighbourNumber));
1159 }
1160 }
1161 }
1162
1163 gapsLifting(gaps, iq);
1164 }
1165
1166 unsigned ncracks = cracks.size();
1167
1168 for (unsigned icrack = 0; icrack < ncracks; ++icrack) {
1169
1170 cracksPads(cracks[icrack], icrack);
1171 }
1172 }
1173
1174
1175 CLHEP::Hep2Vector& EcalHitMaker::correspondingEdge(neighbour& myneighbour, CaloDirection dir2) {
1176 CaloDirection dir = CaloDirectionOperations::oppositeSide(myneighbour.first);
1177 CaloDirection corner = CaloDirectionOperations::add2d(dir, dir2);
1178
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
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
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
1296 unsigned iubd, idir1, idir2;
1297 CaloDirection diag;
1298 CLHEP::Hep2Vector point(0., 0.);
1299
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
1305
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
1316 }
1317 } else if (ngaps == 4) {
1318
1319
1320
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
1331
1332 }
1333 }
1334
1335 void EcalHitMaker::cracksPads(std::vector<neighbour>& cracks, unsigned iq) {
1336
1337 unsigned ncracks = cracks.size();
1338 CrystalPad& myPad = padsatdepth_[iq];
1339 for (unsigned ic = 0; ic < ncracks; ++ic) {
1340
1341
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
1372 crackpad.setSurvivalProbability(crackPadProbability_);
1373 crackpadsatdepth_.push_back(crackpad);
1374 }
1375
1376 }
1377
1378 bool EcalHitMaker::inside3D(const std::vector<XYZPoint>& corners, const XYZPoint& p) const {
1379
1380
1381
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 }