Back to home page

Project CMSSW displayed by LXR
 
 

    

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 //#include "Geometry/EcalAlgo/interface/EcalBarrelGeometry.h"
0007 //#include "Geometry/EcalAlgo/interface/EcalEndcapGeometry.h"
0008 //#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
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 //#include "FastSimulation/Utilities/interface/Histos.h"
0022 
0023 //#include "Math/GenVector/Transform3D.h"
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   //  myHistos->debug("Constructeur EcalHitMaker");
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   // In some cases, a "dummy" grid, not based on a cell, can be built. The previous variables
0085   // should however be initialized. In such a case onEcal=0
0086   if (!onEcal)
0087     return;
0088 
0089   // Same size in eta-phi
0090   etasize_ = size;
0091   phisize_ = size;
0092 
0093   // Build the grid
0094   // The result is put in CellsWindow and is ordered by distance to the pivot
0095   myCalorimeter->getWindow(pivot_.getDetId(), size, size, CellsWindow_);
0096 
0097   buildGeometry();
0098   //  std::cout << " Geometry built " << regionOfInterest_.size() << std::endl;
0099 
0100   truncatedGrid_ = CellsWindow_.size() != (etasize_ * phisize_);
0101 
0102   // A local vector of corners
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   //  std::cout << " Add hit depth called; Current deph is "  << currentdepth_;
0135   //  std::cout << " Required depth is " << depth << std::endl;
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   //  if(cellid.isZero()) std::cout << " cell is Zero " << std::endl;
0143   //  if(xtal<1000)
0144   //    {
0145   //      std::cout << "Result " << regionOfInterest_[xtal].getX0Back() << " " ;
0146   //      std::cout << depth << std::endl;
0147   //    }
0148   //  myHistos->fill("h5000",depth);
0149   if (xtal < 1000) {
0150     //      myHistos->fill("h5002",regionOfInterest_[xtal].getX0Back(),depth);
0151     //      myHistos->fill("h5003",ecalentrance_.eta(),maxX0_);
0152     if (regionOfInterest_[xtal].getX0Back() > depth) {
0153       hits_[xtal] += spotEnergy;
0154       //      myHistos->fill("h5005",r);
0155       return true;
0156     } else {
0157       rearleakage_ += spotEnergy;
0158     }
0159   } else {
0160     //      std::cout << " Return false " << std::endl;
0161     //      std::cout << " Add hit depth called; Current deph is "  << currentdepth_;
0162     //      std::cout << " Required depth is " << depth << std::endl;
0163     //      std::cout << " R = " << r << " " << radiusFactor_ << std::endl;
0164   }
0165 
0166   outsideWindowEnergy_ += spotEnergy;
0167   return false;
0168 }
0169 
0170 //bool EcalHitMaker::addHitDepth(double r,double phi,double depth)
0171 //{
0172 //  std::map<unsigned,float>::iterator itcheck=hitMap_.find(pivot_.getDetId().rawId());
0173 //  if(itcheck==hitMap_.end())
0174 //    {
0175 //      hitMap_.insert(std::pair<uint32_t,float>(pivot_.getDetId().rawId(),spotEnergy));
0176 //    }
0177 //  else
0178 //    {
0179 //      itcheck->second+=spotEnergy;
0180 //    }
0181 //  return true;
0182 //}
0183 
0184 bool EcalHitMaker::addHit(double r, double phi, unsigned layer) {
0185   //  std::cout <<" Addhit " << std::endl;
0186   //  std::cout << " Before insideCell " << std::endl;
0187   double sp(1.);
0188   //  std::cout << " Trying to add " << r << " " << phi << " " << radiusFactor_ << std::endl;
0189   r *= radiusFactor_;
0190   CLHEP::Hep2Vector point(r * std::cos(phi), r * std::sin(phi));
0191   //  std::cout << "point " << point << std::endl;
0192   //  CellID cellid=insideCell(point,sp);
0193   unsigned xtal = fastInsideCell(point, sp);
0194   //  if(cellid.isZero()) std::cout << " cell is Zero " << std::endl;
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   //  std::cout << " This hit ; r= " << point << " hasn't been added "<<std::endl;
0205   //  std::cout << " Xtal " << xtal << std::endl;
0206   //  for(unsigned ip=0;ip<npadsatdepth_;++ip)
0207   //    {
0208   //      std::cout << padsatdepth_[ip] << std::endl;
0209   //    }
0210 
0211   return false;
0212 }
0213 
0214 // Temporary solution
0215 //bool EcalHitMaker::addHit(double r,double phi,unsigned layer)
0216 //{
0217 //  std::map<unsigned,float>::iterator itcheck=hitMap_.find(pivot_.getDetId().rawId());
0218 //  if(itcheck==hitMap_.end())
0219 //    {
0220 //      hitMap_.insert(std::pair<uint32_t,float>(pivot_.getDetId().rawId(),spotEnergy));
0221 //    }
0222 //  else
0223 //    {
0224 //      itcheck->second+=spotEnergy;
0225 //    }
0226 //  return true;
0227 //}
0228 
0229 unsigned EcalHitMaker::fastInsideCell(const CLHEP::Hep2Vector& point, double& sp, bool debug) {
0230   //  debug = true;
0231   bool found = false;
0232   unsigned niter = 0;
0233   // something clever has to be implemented here
0234   unsigned d1, d2;
0235   convertIntegerCoordinates(point.x(), point.y(), d1, d2);
0236   //  std::cout << "Fastinside cell " << point.x() << " " <<  point.y() << " " << d1 << " "<< d2 << " " << nx_ << " " << ny_ << std::endl;
0237   if (d1 >= nx_ || d2 >= ny_) {
0238     //      std::cout << " Not in the map " <<std::endl;
0239     return 9999;
0240   }
0241   unsigned cell = myCrystalNumberArray_[d1][d2];
0242   // We are likely to be lucky
0243   //  std::cout << " Got the cell " << cell << std::endl;
0244   if (validPads_[cell] && padsatdepth_[cell].inside(point)) {
0245     //      std::cout << " We are lucky " << cell << std::endl;
0246     sp = padsatdepth_[cell].survivalProbability();
0247     return cell;
0248   }
0249 
0250   //  std::cout << "Starting the loop " << std::endl;
0251   bool status(true);
0252   const std::vector<unsigned>& localCellVector(myCrystalWindowMap_->getCrystalWindow(cell, status));
0253   if (status) {
0254     unsigned size = localCellVector.size();
0255     //      std::cout << " Starting from " << EBDetId(regionOfInterest_[cell].getDetId()) << std::endl;
0256     //      const std::vector<DetId>& neighbours=myCalorimeter->getNeighbours(regionOfInterest_[cell].getDetId());
0257     //      std::cout << " The neighbours are " << std::endl;
0258     //      for(unsigned ic=0;ic<neighbours.size(); ++ic)
0259     //  {
0260     //    std::cout << EBDetId(neighbours[ic]) << std::endl;
0261     //  }
0262     //      std::cout << " Done " << std::endl;
0263     for (unsigned ic = 0; ic < 8 && ic < size; ++ic) {
0264       unsigned iq = localCellVector[ic];
0265       //      std::cout << " Testing " << EBDetId(regionOfInterest_[iq].getDetId()) << std::endl; ;
0266       //      std::cout << " " << iq << std::endl;
0267       //      std::cout << padsatdepth_[iq] ;
0268       if (validPads_[iq] && padsatdepth_[iq].inside(point)) {
0269         //        std::cout << " Yes " << std::endl;
0270         //        myHistos->fill("h1000",niter);
0271         sp = padsatdepth_[iq].survivalProbability();
0272         //        std::cout << "Finished the loop " << niter << std::endl;
0273         //        std::cout << "Inside " << std::endl;
0274         return iq;
0275       }
0276       //      std::cout << " Not inside " << std::endl;
0277       //      std::cout << "No " << std::endl;
0278       ++niter;
0279     }
0280   }
0281   if (debug)
0282     std::cout << " not found in a quad, let's check the " << ncrackpadsatdepth_ << " cracks " << std::endl;
0283   //  std::cout << "Finished the loop " << niter << std::endl;
0284   // Let's check the cracks
0285   //  std::cout << " Let's check the cracks " << ncrackpadsatdepth_ << " " << crackpadsatdepth_.size() << std::endl;
0286   unsigned iquad = 0;
0287   unsigned iquadinside = 999;
0288   while (iquad < ncrackpadsatdepth_ && !found) {
0289     //      std::cout << " Inside the while " << std::endl;
0290     if (crackpadsatdepth_[iquad].inside(point)) {
0291       iquadinside = iquad;
0292       found = true;
0293       sp = crackpadsatdepth_[iquad].survivalProbability();
0294     }
0295     ++iquad;
0296     ++niter;
0297   }
0298   //  myHistos->fill("h1002",niter);
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   //  myHistos->debug("setTrackParameters");
0306   //  std::cout << " Track " << theTrack << std::endl;
0307   intersections_.clear();
0308   // This is certainly enough
0309   intersections_.reserve(50);
0310   myTrack_ = &theTrack;
0311   normal_ = normal.Unit();
0312   X0depthoffset_ = X0depthoffset;
0313   cellLine(intersections_);
0314   buildSegments(intersections_);
0315   //  std::cout << " Segments " << segments_.size() << std::endl;
0316   //  for(unsigned ii=0; ii<segments_.size() ; ++ii)
0317   //    {
0318   //      std::cout << segments_[ii] << std::endl;
0319   //    }
0320 
0321   // This is only needed in case of electromagnetic showers
0322   if (EMSHOWER && onEcal_ && ecalTotalX0() > 0.) {
0323     //      std::cout << "Total X0  " << ecalTotalX0() << std::endl;
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       // Now calculate the distance in X0 of the back sides of the crystals
0330       // (only for EM showers)
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         // First identify the correct segment
0337         //      std::cout << " Crystal " << ic  << regionOfInterest_[ic].getBackCenter() ;
0338         //      std::cout << " Entrance : " << segments_[ecalFirstSegment_].entrance()<< std::endl;
0339         //      std::cout << " Looking for the segment " << dist << std::endl;
0340         segiterator = find_if(segments_.begin(), segments_.end(), CaloSegment::inSegment(absciss));
0341         //      std::cout << " Done " << std::endl;
0342         if (segiterator == segments_.end()) {
0343           // in this case, we won't have any problem. No need to
0344           // calculate the real depth.
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         //  myHistos->fill("h4000",ecalentrance_.eta(), regionOfInterest_[ic].getX0Back());
0358         //}
0359         //      else
0360         //{
0361         //   // in this case, we won't have any problem. No need to
0362         //  // calculate the real depth.
0363         //  regionOfInterest_[ic].setX0Back(9999);
0364         //}
0365       }  //EMSHOWER
0366     }    // ndir
0367     //      myHistos->fill("h6000",segments_[ecalFirstSegment_].entrance().eta(),maxX0_);
0368   }
0369   //  std::cout << "Leaving setTrackParameters" << std::endl
0370 }
0371 
0372 void EcalHitMaker::cellLine(std::vector<CaloPoint>& cp) {
0373   cp.clear();
0374   //  if(myTrack->onVFcal()!=2)
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   //sort the points by distance (in the ECAL they are not necessarily ordered)
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   // Move the vertex for distance comparison (5cm)
0406   vertex -= 5. * dir;
0407   CaloPoint::DistanceToVertex myDistance(vertex);
0408   sort(cp.begin(), cp.end(), myDistance);
0409 
0410   // The intersections with the HCAL shouldn't need to be sorted
0411   // with the N.I it is actually a source of problems
0412   hcalCellLine(cp);
0413 
0414   //  std::cout << " Intersections ordered by distance to " << vertex << std::endl;
0415   //
0416   //  for (unsigned ic=0;ic<cp.size();++ic)
0417   //    {
0418   //      XYZVector t=cp[ic]-vertex;
0419   //      std::cout << cp[ic] << " " << t.mag() << std::endl;
0420   //    }
0421 }
0422 
0423 void EcalHitMaker::preshowerCellLine(std::vector<CaloPoint>& cp) const {
0424   //  FSimEvent& mySimEvent = myEventMgr->simSignal();
0425   //  FSimTrack myTrack = mySimEvent.track(fsimtrack_);
0426   //  std::cout << "FsimTrack " << fsimtrack_<< std::endl;
0427   //  std::cout << " On layer 1 " << myTrack.onLayer1() << std::endl;
0428   // std::cout << " preshowerCellLine " << std::endl;
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       //      std::cout << "Track on ECAL " << myTrack.EcalEntrance_().vertex()*0.1<< std::endl;
0444     }
0445   }
0446 
0447   //  std::cout << " On layer 2 " << myTrack.onLayer2() << std::endl;
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       //      std::cout << "Track on ECAL " << myTrack.EcalEntrance_().vertex()*0.1 << std::endl;
0463     }
0464   }
0465   //  std::cout << " Exit preshower CellLine " << std::endl;
0466 }
0467 
0468 void EcalHitMaker::hcalCellLine(std::vector<CaloPoint>& cp) const {
0469   //  FSimEvent& mySimEvent = myEventMgr->simSignal();
0470   //  FSimTrack myTrack = mySimEvent.track(fsimtrack_);
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     // HCAL thickness in cm (assuming that the particle is coming from 000)
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     // HCAL thickness in cm (assuming that the particle is coming from 000)
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   //  std::vector<XYZPoint> corners;
0502   //  corners.resize(4);
0503   unsigned ic = 0;
0504   double t;
0505   XYZPoint xp;
0506   DetId c_entrance, c_exit;
0507   bool entrancefound(false), exitfound(false);
0508   //  std::cout << " Look for intersections " << ncrystals_ << std::endl;
0509   //  std::cout << " regionOfInterest_ " << truncatedGrid_ << " " << regionOfInterest_.size() << std::endl;
0510   // try to determine the number of crystals to test
0511   // First determine the incident angle
0512   double angle = std::acos(normal_.Dot(regionOfInterest_[0].getAxis().Unit()));
0513 
0514   //  std::cout << " Normal " << normal_<< " Axis " << regionOfInterest_[0].getAxis().Unit() << std::endl;
0515   double backdistance = std::sqrt(regionOfInterest_[0].getAxis().mag2()) * std::tan(angle);
0516   // 1/2.2cm = 0.45
0517   //   std::cout << " Angle " << angle << std::endl;
0518   //   std::cout << " Back distance " << backdistance << std::endl;
0519   unsigned ncrystals = (unsigned)(backdistance * 0.45);
0520   unsigned highlim = (ncrystals + 4);
0521   highlim *= highlim;
0522   if (highlim > ncrystals_)
0523     highlim = ncrystals_;
0524   //  unsigned lowlim=(ncrystals>2)? (ncrystals-2):0;
0525   //  std::cout << " Ncrys " << ncrystals << std::endl;
0526 
0527   while (ic < ncrystals_ && (ic < highlim || !exitfound)) {
0528     // Check front side
0529     //      if(!entrancefound)
0530     {
0531       const Plane3D& plan = regionOfInterest_[ic].getFrontPlane();
0532       //      XYZVector axis1=(plan.Normal());
0533       //      XYZVector axis2=regionOfInterest_[ic].getFirstEdge();
0534       xp = intersect(plan, a, b, t, false);
0535       regionOfInterest_[ic].getFrontSide(corners);
0536       //      CrystalPad pad(9999,onEcal_,corners,regionOfInterest_[ic].getCorner(0),axis1,axis2);
0537       //      if(pad.globalinside(xp))
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         //      myHistos->fill("j12",highlim,ic);
0543       }
0544     }
0545 
0546     // check rear side
0547     //  if(!exitfound)
0548     {
0549       const Plane3D& plan = regionOfInterest_[ic].getBackPlane();
0550       //      XYZVector axis1=(plan.Normal());
0551       //      XYZVector axis2=regionOfInterest_[ic].getFifthEdge();
0552       xp = intersect(plan, a, b, t, false);
0553       regionOfInterest_[ic].getBackSide(corners);
0554       //      CrystalPad pad(9999,onEcal_,corners,regionOfInterest_[ic].getCorner(4),axis1,axis2);
0555       //      if(pad.globalinside(xp))
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         //        std::cout << " Crystal : " << ic << std::endl;
0561         //        myHistos->fill("j10",highlim,ic);
0562       }
0563     }
0564 
0565     if (entrancefound && exitfound && c_entrance == c_exit)
0566       return;
0567     // check lateral sides
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       //      XYZVector axis1=(plan.Normal());
0572       //      XYZVector axis2=regionOfInterest_[ic].getLateralEdge(iside);
0573       regionOfInterest_[ic].getLateralSide(iside, corners);
0574       //      CrystalPad pad(9999,onEcal_,corners,regionOfInterest_[ic].getCorner(iside),axis1,axis2);
0575       //      if(pad.globalinside(xp))
0576       if (inside3D(corners, xp)) {
0577         cp.push_back(CaloPoint(regionOfInterest_[ic].getDetId(), CaloDirectionOperations::Side(iside), xp));
0578         //        std::cout << cp[cp.size()-1] << std::endl;
0579       }
0580     }
0581     // Go to next crystal
0582     ++ic;
0583   }
0584 }
0585 
0586 void EcalHitMaker::buildSegments(const std::vector<CaloPoint>& cp) {
0587   //  myHistos->debug();
0588   //  TimeMe theT("FamosGrid::buildSegments");
0589   unsigned size = cp.size();
0590   if (size % 2 != 0) {
0591     //      std::cout << " There is a problem " << std::endl;
0592     return;
0593   }
0594   //  myHistos->debug();
0595   unsigned nsegments = size / 2;
0596   segments_.reserve(nsegments);
0597   if (size == 0)
0598     return;
0599   // curv abs
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       //      std::cout << " Problem with the segments " << std::endl;
0610       //      std::cout << cp[2*is].whichDetector() << " " << cp[2*is+1].whichDetector() << std::endl;
0611       //      std::cout << is << " " <<cp[2*is].getDetId().rawId() << std::endl;
0612       //      std::cout << (2*is+1) << " " <<cp[2*is+1].getDetId().rawId() << std::endl;
0613       ++is;
0614       continue;
0615     }
0616 
0617     // Check if it is a Preshower segment - Layer 1
0618     // One segment per layer, nothing between
0619     //      myHistos->debug("Just avant Preshower");
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         //        std::cout << " Added (1-1)" << preshsegment << std::endl;
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     // Check if it is a Preshower segment - Layer 2
0641     // One segment per layer, nothing between
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         //        std::cout << " Added (1-2)" << preshsegment << std::endl;
0649         s += preshsegment.length();
0650         sX0 += preshsegment.X0length();
0651         sL0 += preshsegment.L0length();
0652         X0PS2_ += preshsegment.X0length();
0653         L0PS2_ += preshsegment.L0length();
0654 
0655         // material between preshower and EE
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           //          std::cout << " Created  a segment " << gapsef.length()<< " " << gapsef.X0length()<< std::endl;
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     // Now deal with the ECAL
0675     // One segment in each crystal. Segment corresponding to cracks/gaps are added
0676     //      myHistos->debug("Just avant ECAL");
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         // set the real entrance
0683         if (ecalFirstSegment_ < 0)
0684           ecalFirstSegment_ = segments_.size();
0685 
0686         // !! Approximatiom : the first segment is always in a crystal
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           // std::cout << " Added (2)" << segment << std::endl;
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         // Now check if a gap or crack should be added
0703         if (is > 0 && is < nsegments) {
0704           DetId cell3 = cp[2 * is].getDetId();
0705           if (cp[2 * is].whichDetector() != DetId::Hcal) {
0706             // Crack inside the ECAL
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             //   std::cout <<" Added(3) "<< cracksegment << std::endl;
0722           } else {
0723             // a segment corresponding to ECAL/HCAL transition should be
0724             // added here
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     //      myHistos->debug("Just avant HCAL");
0743     // HCAL
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       //      std::cout <<" Added(4) "<< segment << std::endl;
0753       ++is;
0754     }
0755   }
0756   //  std::cout << " PS1 " << X0PS1_ << " " << L0PS1_ << std::endl;
0757   //  std::cout << " PS2 " << X0PS2_ << " " << L0PS2_ << std::endl;
0758   //  std::cout << " ECAL " << X0ECAL_ << " " << L0ECAL_ << std::endl;
0759   //  std::cout << " HCAL " << X0HCAL_ << " " << L0HCAL_ << std::endl;
0760 
0761   totalX0_ = X0PS1_ + X0PS2_ + X0PS2EE_ + X0ECAL_ + X0EHGAP_ + X0HCAL_;
0762   totalL0_ = L0PS1_ + L0PS2_ + L0PS2EE_ + L0ECAL_ + L0EHGAP_ + L0HCAL_;
0763   //  myHistos->debug("Just avant le fill");
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   //  if(X0ECAL_<1.)
0771   //    {
0772   //      for(unsigned ii=0; ii<segments_.size() ; ++ii)
0773   //    {
0774   //      std::cout << segments_[ii] << std::endl;
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   //  std::cout << " Finished the segments " << std::endl;
0787 }
0788 
0789 void EcalHitMaker::buildGeometry() {
0790   configuredGeometry_ = false;
0791   ncrystals_ = CellsWindow_.size();
0792   // create the vector with of pads with the appropriate size
0793   padsatdepth_.resize(ncrystals_);
0794 
0795   // This is fully correct in the barrel.
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   // Computes the map of the neighbours
0816   myCrystalWindowMap_ = new CrystalWindowMap(myCalorimeter, regionOfInterest_);
0817 }
0818 
0819 // depth is in X0 , L0 (depending on EMSHOWER/HADSHOWER) or in CM if inCM
0820 bool EcalHitMaker::getPads(double depth, bool inCm) {
0821   //std::cout << " New depth " << depth << std::endl;
0822   // The first time, the relationship between crystals must be calculated
0823   // but only in the case of EM showers
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   //  if(currentdepth_>maxX0_+ps1TotalX0()+ps2TotalX0())
0836   //    {
0837   //      currentdepth_=maxX0_+ps1TotalX0()+ps2TotalX0()-1.; // the -1 is for safety
0838   //      detailedShowerTail_=true;
0839   //    }
0840 
0841   //  std::cout << " FamosGrid::getQuads " << currentdepth_ << " " << maxX0_ << std::endl;
0842 
0843   ncrackpadsatdepth_ = 0;
0844 
0845   xmin_ = ymin_ = 999;
0846   xmax_ = ymax_ = -999;
0847   double locxmin, locxmax, locymin, locymax;
0848 
0849   // Get the depth of the pivot
0850   std::vector<CaloSegment>::const_iterator segiterator;
0851   // First identify the correct segment
0852 
0853   if (inCm)  // centimeter
0854   {
0855     segiterator = find_if(segments_.begin(), segments_.end(), CaloSegment::inSegment(currentdepth_));
0856   } else {
0857     // EM shower
0858     if (EMSHOWER)
0859       segiterator = find_if(segments_.begin(), segments_.end(), CaloSegment::inX0Segment(currentdepth_));
0860 
0861     //Hadron shower
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   //  std::cout << *segiterator << std::endl;
0877 
0878   if (segiterator->whichDetector() != DetId::Ecal) {
0879     std::cout << " In  " << segiterator->whichDetector() << std::endl;
0880     //      buildPreshower();
0881     return false;
0882   }
0883 
0884   //  std::cout << *segiterator << std::endl;
0885   // get the position of the origin
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   //  std::cout << " currentdepth_ " << currentdepth_ << " " << origin << std::endl;
0897   XYZVector newaxis = pivot_.getFirstEdge().Cross(normal_);
0898 
0899   //  std::cout  << "Normal " << normal_ << std::endl;
0900   //  std::cout << " New axis " << newaxis << std::endl;
0901 
0902   //  std::cout << " ncrystals  " << ncrystals << std::endl;
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     //      std::cout << " Building geometry for " << regionOfInterest_[ic].getCellID() << std::endl;
0915     XYZPoint a, b;
0916 
0917     //      std::cout << " Origin " << origin << std::endl;
0918 
0919     //      std::vector<XYZPoint> corners;
0920     //      corners.reserve(4);
0921     double dummyt;
0922     bool hasbeenpulled = false;
0923     bool behindback = false;
0924     for (unsigned il = 0; il < 4; ++il) {
0925       // a is the il-th front corner of the crystal. b is the corresponding rear corner
0926       regionOfInterest_[ic].getLateralEdges(il, a, b);
0927 
0928       // pull the surface if necessary (only in the front of the crystals)
0929       XYZPoint aprime = a;
0930       if (pulled(origin, normal_, a)) {
0931         b = aprime;
0932         hasbeenpulled = true;
0933       }
0934 
0935       // compute the intersection.
0936       // Check that the intersection is in the [a,b] segment  if HADSHOWER
0937       // if EMSHOWER the intersection is calculated as if the crystals were infinite
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       //      std::cout << " Intersect " << il << " " << a << " " << b << " " << plan_ << " " << xx << std::endl;
0943       // check that the intersection actually exists
0944       if (xx.mag2() != 0) {
0945         corners[il] = xx;
0946       }
0947     }
0948     //      std::cout << " ncorners " << corners.size() << std::endl;
0949     if (behindback && EMSHOWER)
0950       detailedShowerTail_ = true;
0951     // If the quad is completly defined. Store it !
0952     if (corners.size() == 4) {
0953       padsatdepth_[ic] = CrystalPad(ic, corners, trans, bfactor_, !central_);
0954       // Parameter to be tuned
0955       if (hasbeenpulled)
0956         padsatdepth_[ic].setSurvivalProbability(pulledPadProbability_);
0957       validPads_[ic] = true;
0958       ++nquads;
0959       // In principle, this should be done after the quads reorganization. But it would cost one more loop
0960       //  quadsatdepth_[ic].extrems(locxmin,locxmax,locymin,locymax);
0961       //  if(locxmin<xmin_) xmin_=locxmin;
0962       //  if(locymin<ymin_) ymin_=locymin;
0963       //  if(locxmax>xmax_) xmax_=locxmax;
0964       //  if(locymax>ymax_) ymax_=locymax;
0965     } else {
0966       padsatdepth_[ic] = CrystalPad();
0967       validPads_[ic] = false;
0968     }
0969   }
0970   //  std::cout << " Number of quads " << quadsatdepth_.size() << std::endl;
0971   if (doreorg_)
0972     reorganizePads();
0973   //  std::cout << "Finished to reorganize " << std::endl;
0974   npadsatdepth_ = nquads;
0975   //  std::cout << " prepareCellIDMap " << std::endl;
0976 
0977   // Resize the Quads to allow for some numerical inaccuracy
0978   // in the "inside" function
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   // 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 
1120 ///////////////////////////// GAPS/CRACKS TREATMENT////////////
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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!