Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:11:13 

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   // something clever has to be implemented here
0233   unsigned d1, d2;
0234   convertIntegerCoordinates(point.x(), point.y(), d1, d2);
0235   //  std::cout << "Fastinside cell " << point.x() << " " <<  point.y() << " " << d1 << " "<< d2 << " " << nx_ << " " << ny_ << std::endl;
0236   if (d1 >= nx_ || d2 >= ny_) {
0237     //      std::cout << " Not in the map " <<std::endl;
0238     return 9999;
0239   }
0240   unsigned cell = myCrystalNumberArray_[d1][d2];
0241   // We are likely to be lucky
0242   //  std::cout << " Got the cell " << cell << std::endl;
0243   if (validPads_[cell] && padsatdepth_[cell].inside(point)) {
0244     //      std::cout << " We are lucky " << cell << std::endl;
0245     sp = padsatdepth_[cell].survivalProbability();
0246     return cell;
0247   }
0248 
0249   //  std::cout << "Starting the loop " << std::endl;
0250   bool status(true);
0251   const std::vector<unsigned>& localCellVector(myCrystalWindowMap_->getCrystalWindow(cell, status));
0252   if (status) {
0253     unsigned size = localCellVector.size();
0254     //      std::cout << " Starting from " << EBDetId(regionOfInterest_[cell].getDetId()) << std::endl;
0255     //      const std::vector<DetId>& neighbours=myCalorimeter->getNeighbours(regionOfInterest_[cell].getDetId());
0256     //      std::cout << " The neighbours are " << std::endl;
0257     //      for(unsigned ic=0;ic<neighbours.size(); ++ic)
0258     //  {
0259     //    std::cout << EBDetId(neighbours[ic]) << std::endl;
0260     //  }
0261     //      std::cout << " Done " << std::endl;
0262     for (unsigned ic = 0; ic < 8 && ic < size; ++ic) {
0263       unsigned iq = localCellVector[ic];
0264       //      std::cout << " Testing " << EBDetId(regionOfInterest_[iq].getDetId()) << std::endl; ;
0265       //      std::cout << " " << iq << std::endl;
0266       //      std::cout << padsatdepth_[iq] ;
0267       if (validPads_[iq] && padsatdepth_[iq].inside(point)) {
0268         //        std::cout << " Yes " << std::endl;
0269         //        myHistos->fill("h1000",niter);
0270         sp = padsatdepth_[iq].survivalProbability();
0271         //        std::cout << "Finished the loop " << niter << std::endl;
0272         //        std::cout << "Inside " << std::endl;
0273         return iq;
0274       }
0275       //      std::cout << " Not inside " << std::endl;
0276       //      std::cout << "No " << std::endl;
0277     }
0278   }
0279   if (debug)
0280     std::cout << " not found in a quad, let's check the " << ncrackpadsatdepth_ << " cracks " << std::endl;
0281   //  std::cout << "Finished the loop " << niter << std::endl;
0282   // Let's check the cracks
0283   //  std::cout << " Let's check the cracks " << ncrackpadsatdepth_ << " " << crackpadsatdepth_.size() << std::endl;
0284   unsigned iquad = 0;
0285   unsigned iquadinside = 999;
0286   while (iquad < ncrackpadsatdepth_ && !found) {
0287     //      std::cout << " Inside the while " << std::endl;
0288     if (crackpadsatdepth_[iquad].inside(point)) {
0289       iquadinside = iquad;
0290       found = true;
0291       sp = crackpadsatdepth_[iquad].survivalProbability();
0292     }
0293     ++iquad;
0294   }
0295   //  myHistos->fill("h1002",niter);
0296   if (!found && debug)
0297     std::cout << " Not found in the cracks " << std::endl;
0298   return (found) ? crackpadsatdepth_[iquadinside].getNumber() : 9999;
0299 }
0300 
0301 void EcalHitMaker::setTrackParameters(const XYZNormal& normal, double X0depthoffset, const FSimTrack& theTrack) {
0302   //  myHistos->debug("setTrackParameters");
0303   //  std::cout << " Track " << theTrack << std::endl;
0304   intersections_.clear();
0305   // This is certainly enough
0306   intersections_.reserve(50);
0307   myTrack_ = &theTrack;
0308   normal_ = normal.Unit();
0309   X0depthoffset_ = X0depthoffset;
0310   cellLine(intersections_);
0311   buildSegments(intersections_);
0312   //  std::cout << " Segments " << segments_.size() << std::endl;
0313   //  for(unsigned ii=0; ii<segments_.size() ; ++ii)
0314   //    {
0315   //      std::cout << segments_[ii] << std::endl;
0316   //    }
0317 
0318   // This is only needed in case of electromagnetic showers
0319   if (EMSHOWER && onEcal_ && ecalTotalX0() > 0.) {
0320     //      std::cout << "Total X0  " << ecalTotalX0() << std::endl;
0321     for (unsigned ic = 0; ic < ncrystals_; ++ic) {
0322       for (unsigned idir = 0; idir < 4; ++idir) {
0323         XYZVector norm = regionOfInterest_[ic].exitingNormal(CaloDirectionOperations::Side(idir));
0324         regionOfInterest_[ic].crystalNeighbour(idir).setToBeGlued((norm.Dot(normal_) < 0.));
0325       }
0326       // Now calculate the distance in X0 of the back sides of the crystals
0327       // (only for EM showers)
0328       if (EMSHOWER) {
0329         XYZVector dir = regionOfInterest_[ic].getBackCenter() - segments_[ecalFirstSegment_].entrance();
0330         double dist = dir.Dot(normal_);
0331         double absciss = dist + segments_[ecalFirstSegment_].sEntrance();
0332         std::vector<CaloSegment>::const_iterator segiterator;
0333         // First identify the correct segment
0334         //      std::cout << " Crystal " << ic  << regionOfInterest_[ic].getBackCenter() ;
0335         //      std::cout << " Entrance : " << segments_[ecalFirstSegment_].entrance()<< std::endl;
0336         //      std::cout << " Looking for the segment " << dist << std::endl;
0337         segiterator = find_if(segments_.begin(), segments_.end(), CaloSegment::inSegment(absciss));
0338         //      std::cout << " Done " << std::endl;
0339         if (segiterator == segments_.end()) {
0340           // in this case, we won't have any problem. No need to
0341           // calculate the real depth.
0342           regionOfInterest_[ic].setX0Back(9999);
0343         } else {
0344           DetId::Detector det(segiterator->whichDetector());
0345           if (det != DetId::Ecal) {
0346             regionOfInterest_[ic].setX0Back(9999);
0347           } else {
0348             double x0 = segiterator->x0FromCm(dist);
0349             if (x0 < maxX0_)
0350               maxX0_ = x0;
0351             regionOfInterest_[ic].setX0Back(x0);
0352           }
0353         }
0354         //  myHistos->fill("h4000",ecalentrance_.eta(), regionOfInterest_[ic].getX0Back());
0355         //}
0356         //      else
0357         //{
0358         //   // in this case, we won't have any problem. No need to
0359         //  // calculate the real depth.
0360         //  regionOfInterest_[ic].setX0Back(9999);
0361         //}
0362       }  //EMSHOWER
0363     }    // ndir
0364     //      myHistos->fill("h6000",segments_[ecalFirstSegment_].entrance().eta(),maxX0_);
0365   }
0366   //  std::cout << "Leaving setTrackParameters" << std::endl
0367 }
0368 
0369 void EcalHitMaker::cellLine(std::vector<CaloPoint>& cp) {
0370   cp.clear();
0371   //  if(myTrack->onVFcal()!=2)
0372   //    {
0373   if (!central_ && onEcal_ && simulatePreshower_)
0374     preshowerCellLine(cp);
0375   if (onEcal_)
0376     ecalCellLine(EcalEntrance_, EcalEntrance_ + normal_, cp);
0377   //    }
0378 
0379   XYZPoint vertex(myTrack_->vertex().position().Vect());
0380 
0381   //sort the points by distance (in the ECAL they are not necessarily ordered)
0382   XYZVector dir(0., 0., 0.);
0383   if (myTrack_->onLayer1()) {
0384     vertex = (myTrack_->layer1Entrance().vertex()).Vect();
0385     dir = myTrack_->layer1Entrance().Vect().Unit();
0386   } else if (myTrack_->onLayer2()) {
0387     vertex = (myTrack_->layer2Entrance().vertex()).Vect();
0388     dir = myTrack_->layer2Entrance().Vect().Unit();
0389   } else if (myTrack_->onEcal()) {
0390     vertex = (myTrack_->ecalEntrance().vertex()).Vect();
0391     dir = myTrack_->ecalEntrance().Vect().Unit();
0392   } else if (myTrack_->onHcal()) {
0393     vertex = (myTrack_->hcalEntrance().vertex()).Vect();
0394     dir = myTrack_->hcalEntrance().Vect().Unit();
0395   } else if (myTrack_->onVFcal() == 2) {
0396     vertex = (myTrack_->vfcalEntrance().vertex()).Vect();
0397     dir = myTrack_->vfcalEntrance().Vect().Unit();
0398   } else {
0399     std::cout << " Problem with the grid " << std::endl;
0400   }
0401 
0402   // Move the vertex for distance comparison (5cm)
0403   vertex -= 5. * dir;
0404   CaloPoint::DistanceToVertex myDistance(vertex);
0405   sort(cp.begin(), cp.end(), myDistance);
0406 
0407   // The intersections with the HCAL shouldn't need to be sorted
0408   // with the N.I it is actually a source of problems
0409   hcalCellLine(cp);
0410 
0411   //  std::cout << " Intersections ordered by distance to " << vertex << std::endl;
0412   //
0413   //  for (unsigned ic=0;ic<cp.size();++ic)
0414   //    {
0415   //      XYZVector t=cp[ic]-vertex;
0416   //      std::cout << cp[ic] << " " << t.mag() << std::endl;
0417   //    }
0418 }
0419 
0420 void EcalHitMaker::preshowerCellLine(std::vector<CaloPoint>& cp) const {
0421   //  FSimEvent& mySimEvent = myEventMgr->simSignal();
0422   //  FSimTrack myTrack = mySimEvent.track(fsimtrack_);
0423   //  std::cout << "FsimTrack " << fsimtrack_<< std::endl;
0424   //  std::cout << " On layer 1 " << myTrack.onLayer1() << std::endl;
0425   // std::cout << " preshowerCellLine " << std::endl;
0426   if (myTrack_->onLayer1()) {
0427     XYZPoint point1 = (myTrack_->layer1Entrance().vertex()).Vect();
0428     double phys_eta = myTrack_->layer1Entrance().eta();
0429     double cmthickness = myCalorimeter->layer1Properties(1)->thickness(phys_eta);
0430 
0431     if (cmthickness > 0) {
0432       XYZVector dir = myTrack_->layer1Entrance().Vect().Unit();
0433       XYZPoint point2 = point1 + dir * cmthickness;
0434 
0435       CaloPoint cp1(DetId::Ecal, EcalPreshower, 1, point1);
0436       CaloPoint cp2(DetId::Ecal, EcalPreshower, 1, point2);
0437       cp.push_back(cp1);
0438       cp.push_back(cp2);
0439     } else {
0440       //      std::cout << "Track on ECAL " << myTrack.EcalEntrance_().vertex()*0.1<< std::endl;
0441     }
0442   }
0443 
0444   //  std::cout << " On layer 2 " << myTrack.onLayer2() << std::endl;
0445   if (myTrack_->onLayer2()) {
0446     XYZPoint point1 = (myTrack_->layer2Entrance().vertex()).Vect();
0447     double phys_eta = myTrack_->layer2Entrance().eta();
0448     double cmthickness = myCalorimeter->layer2Properties(1)->thickness(phys_eta);
0449     if (cmthickness > 0) {
0450       XYZVector dir = myTrack_->layer2Entrance().Vect().Unit();
0451       XYZPoint point2 = point1 + dir * cmthickness;
0452 
0453       CaloPoint cp1(DetId::Ecal, EcalPreshower, 2, point1);
0454       CaloPoint cp2(DetId::Ecal, EcalPreshower, 2, point2);
0455 
0456       cp.push_back(cp1);
0457       cp.push_back(cp2);
0458     } else {
0459       //      std::cout << "Track on ECAL " << myTrack.EcalEntrance_().vertex()*0.1 << std::endl;
0460     }
0461   }
0462   //  std::cout << " Exit preshower CellLine " << std::endl;
0463 }
0464 
0465 void EcalHitMaker::hcalCellLine(std::vector<CaloPoint>& cp) const {
0466   //  FSimEvent& mySimEvent = myEventMgr->simSignal();
0467   //  FSimTrack myTrack = mySimEvent.track(fsimtrack_);
0468   int onHcal = myTrack_->onHcal();
0469 
0470   if (onHcal <= 2 && onHcal > 0) {
0471     XYZPoint point1 = (myTrack_->hcalEntrance().vertex()).Vect();
0472 
0473     double eta = point1.eta();
0474     // HCAL thickness in cm (assuming that the particle is coming from 000)
0475     double thickness = myCalorimeter->hcalProperties(onHcal)->thickness(eta);
0476     cp.push_back(CaloPoint(DetId::Hcal, point1));
0477     XYZVector dir = myTrack_->hcalEntrance().Vect().Unit();
0478     XYZPoint point2 = point1 + dir * thickness;
0479 
0480     cp.push_back(CaloPoint(DetId::Hcal, point2));
0481   }
0482   int onVFcal = myTrack_->onVFcal();
0483   if (onVFcal == 2) {
0484     XYZPoint point1 = (myTrack_->vfcalEntrance().vertex()).Vect();
0485     double eta = point1.eta();
0486     // HCAL thickness in cm (assuming that the particle is coming from 000)
0487     double thickness = myCalorimeter->hcalProperties(3)->thickness(eta);
0488     cp.push_back(CaloPoint(DetId::Hcal, point1));
0489     XYZVector dir = myTrack_->vfcalEntrance().Vect().Unit();
0490     if (thickness > 0) {
0491       XYZPoint point2 = point1 + dir * thickness;
0492       cp.push_back(CaloPoint(DetId::Hcal, point2));
0493     }
0494   }
0495 }
0496 
0497 void EcalHitMaker::ecalCellLine(const XYZPoint& a, const XYZPoint& b, std::vector<CaloPoint>& cp) {
0498   //  std::vector<XYZPoint> corners;
0499   //  corners.resize(4);
0500   unsigned ic = 0;
0501   double t;
0502   XYZPoint xp;
0503   DetId c_entrance, c_exit;
0504   bool entrancefound(false), exitfound(false);
0505   //  std::cout << " Look for intersections " << ncrystals_ << std::endl;
0506   //  std::cout << " regionOfInterest_ " << truncatedGrid_ << " " << regionOfInterest_.size() << std::endl;
0507   // try to determine the number of crystals to test
0508   // First determine the incident angle
0509   double angle = std::acos(normal_.Dot(regionOfInterest_[0].getAxis().Unit()));
0510 
0511   //  std::cout << " Normal " << normal_<< " Axis " << regionOfInterest_[0].getAxis().Unit() << std::endl;
0512   double backdistance = std::sqrt(regionOfInterest_[0].getAxis().mag2()) * std::tan(angle);
0513   // 1/2.2cm = 0.45
0514   //   std::cout << " Angle " << angle << std::endl;
0515   //   std::cout << " Back distance " << backdistance << std::endl;
0516   unsigned ncrystals = (unsigned)(backdistance * 0.45);
0517   unsigned highlim = (ncrystals + 4);
0518   highlim *= highlim;
0519   if (highlim > ncrystals_)
0520     highlim = ncrystals_;
0521   //  unsigned lowlim=(ncrystals>2)? (ncrystals-2):0;
0522   //  std::cout << " Ncrys " << ncrystals << std::endl;
0523 
0524   while (ic < ncrystals_ && (ic < highlim || !exitfound)) {
0525     // Check front side
0526     //      if(!entrancefound)
0527     {
0528       const Plane3D& plan = regionOfInterest_[ic].getFrontPlane();
0529       //      XYZVector axis1=(plan.Normal());
0530       //      XYZVector axis2=regionOfInterest_[ic].getFirstEdge();
0531       xp = intersect(plan, a, b, t, false);
0532       regionOfInterest_[ic].getFrontSide(corners);
0533       //      CrystalPad pad(9999,onEcal_,corners,regionOfInterest_[ic].getCorner(0),axis1,axis2);
0534       //      if(pad.globalinside(xp))
0535       if (inside3D(corners, xp)) {
0536         cp.push_back(CaloPoint(regionOfInterest_[ic].getDetId(), UP, xp));
0537         entrancefound = true;
0538         c_entrance = regionOfInterest_[ic].getDetId();
0539         //      myHistos->fill("j12",highlim,ic);
0540       }
0541     }
0542 
0543     // check rear side
0544     //  if(!exitfound)
0545     {
0546       const Plane3D& plan = regionOfInterest_[ic].getBackPlane();
0547       //      XYZVector axis1=(plan.Normal());
0548       //      XYZVector axis2=regionOfInterest_[ic].getFifthEdge();
0549       xp = intersect(plan, a, b, t, false);
0550       regionOfInterest_[ic].getBackSide(corners);
0551       //      CrystalPad pad(9999,onEcal_,corners,regionOfInterest_[ic].getCorner(4),axis1,axis2);
0552       //      if(pad.globalinside(xp))
0553       if (inside3D(corners, xp)) {
0554         cp.push_back(CaloPoint(regionOfInterest_[ic].getDetId(), DOWN, xp));
0555         exitfound = true;
0556         c_exit = regionOfInterest_[ic].getDetId();
0557         //        std::cout << " Crystal : " << ic << std::endl;
0558         //        myHistos->fill("j10",highlim,ic);
0559       }
0560     }
0561 
0562     if (entrancefound && exitfound && c_entrance == c_exit)
0563       return;
0564     // check lateral sides
0565     for (unsigned iside = 0; iside < 4; ++iside) {
0566       const Plane3D& plan = regionOfInterest_[ic].getLateralPlane(iside);
0567       xp = intersect(plan, a, b, t, false);
0568       //      XYZVector axis1=(plan.Normal());
0569       //      XYZVector axis2=regionOfInterest_[ic].getLateralEdge(iside);
0570       regionOfInterest_[ic].getLateralSide(iside, corners);
0571       //      CrystalPad pad(9999,onEcal_,corners,regionOfInterest_[ic].getCorner(iside),axis1,axis2);
0572       //      if(pad.globalinside(xp))
0573       if (inside3D(corners, xp)) {
0574         cp.push_back(CaloPoint(regionOfInterest_[ic].getDetId(), CaloDirectionOperations::Side(iside), xp));
0575         //        std::cout << cp[cp.size()-1] << std::endl;
0576       }
0577     }
0578     // Go to next crystal
0579     ++ic;
0580   }
0581 }
0582 
0583 void EcalHitMaker::buildSegments(const std::vector<CaloPoint>& cp) {
0584   //  myHistos->debug();
0585   //  TimeMe theT("FamosGrid::buildSegments");
0586   unsigned size = cp.size();
0587   if (size % 2 != 0) {
0588     //      std::cout << " There is a problem " << std::endl;
0589     return;
0590   }
0591   //  myHistos->debug();
0592   unsigned nsegments = size / 2;
0593   segments_.reserve(nsegments);
0594   if (size == 0)
0595     return;
0596   // curv abs
0597   double s = 0.;
0598   double sX0 = 0.;
0599   double sL0 = 0.;
0600 
0601 #ifdef DEBUGCELLLINE
0602   unsigned ncrossedxtals = 0;
0603 #endif
0604   unsigned is = 0;
0605   while (is < nsegments) {
0606     if (cp[2 * is].getDetId() != cp[2 * is + 1].getDetId() && cp[2 * is].whichDetector() != DetId::Hcal &&
0607         cp[2 * is + 1].whichDetector() != DetId::Hcal) {
0608       //      std::cout << " Problem with the segments " << std::endl;
0609       //      std::cout << cp[2*is].whichDetector() << " " << cp[2*is+1].whichDetector() << std::endl;
0610       //      std::cout << is << " " <<cp[2*is].getDetId().rawId() << std::endl;
0611       //      std::cout << (2*is+1) << " " <<cp[2*is+1].getDetId().rawId() << std::endl;
0612       ++is;
0613       continue;
0614     }
0615 
0616     // Check if it is a Preshower segment - Layer 1
0617     // One segment per layer, nothing between
0618     //      myHistos->debug("Just avant Preshower");
0619     if (cp[2 * is].whichDetector() == DetId::Ecal && cp[2 * is].whichSubDetector() == EcalPreshower &&
0620         cp[2 * is].whichLayer() == 1) {
0621       if (cp[2 * is + 1].whichDetector() == DetId::Ecal && cp[2 * is + 1].whichSubDetector() == EcalPreshower &&
0622           cp[2 * is + 1].whichLayer() == 1) {
0623         CaloSegment preshsegment(cp[2 * is], cp[2 * is + 1], s, sX0, sL0, CaloSegment::PS, myCalorimeter);
0624         segments_.push_back(preshsegment);
0625         //        std::cout << " Added (1-1)" << preshsegment << std::endl;
0626         s += preshsegment.length();
0627         sX0 += preshsegment.X0length();
0628         sL0 += preshsegment.L0length();
0629         X0PS1_ += preshsegment.X0length();
0630         L0PS1_ += preshsegment.L0length();
0631       } else {
0632         std::cout << " Strange segment between Preshower1 and " << cp[2 * is + 1].whichDetector();
0633         std::cout << std::endl;
0634       }
0635       ++is;
0636       continue;
0637     }
0638 
0639     // Check if it is a Preshower segment - Layer 2
0640     // One segment per layer, nothing between
0641     if (cp[2 * is].whichDetector() == DetId::Ecal && cp[2 * is].whichSubDetector() == EcalPreshower &&
0642         cp[2 * is].whichLayer() == 2) {
0643       if (cp[2 * is + 1].whichDetector() == DetId::Ecal && cp[2 * is + 1].whichSubDetector() == EcalPreshower &&
0644           cp[2 * is + 1].whichLayer() == 2) {
0645         CaloSegment preshsegment(cp[2 * is], cp[2 * is + 1], s, sX0, sL0, CaloSegment::PS, myCalorimeter);
0646         segments_.push_back(preshsegment);
0647         //        std::cout << " Added (1-2)" << preshsegment << std::endl;
0648         s += preshsegment.length();
0649         sX0 += preshsegment.X0length();
0650         sL0 += preshsegment.L0length();
0651         X0PS2_ += preshsegment.X0length();
0652         L0PS2_ += preshsegment.L0length();
0653 
0654         // material between preshower and EE
0655         if (is < (nsegments - 1) && cp[2 * is + 2].whichDetector() == DetId::Ecal &&
0656             cp[2 * is + 2].whichSubDetector() == EcalEndcap) {
0657           CaloSegment gapsef(cp[2 * is + 1], cp[2 * is + 2], s, sX0, sL0, CaloSegment::PSEEGAP, myCalorimeter);
0658           segments_.push_back(gapsef);
0659           s += gapsef.length();
0660           sX0 += gapsef.X0length();
0661           sL0 += gapsef.L0length();
0662           X0PS2EE_ += gapsef.X0length();
0663           L0PS2EE_ += gapsef.L0length();
0664           //          std::cout << " Created  a segment " << gapsef.length()<< " " << gapsef.X0length()<< std::endl;
0665         }
0666       } else {
0667         std::cout << " Strange segment between Preshower2 and " << cp[2 * is + 1].whichDetector();
0668         std::cout << std::endl;
0669       }
0670       ++is;
0671       continue;
0672     }
0673     // Now deal with the ECAL
0674     // One segment in each crystal. Segment corresponding to cracks/gaps are added
0675     //      myHistos->debug("Just avant ECAL");
0676     if (cp[2 * is].whichDetector() == DetId::Ecal &&
0677         (cp[2 * is].whichSubDetector() == EcalBarrel || cp[2 * is].whichSubDetector() == EcalEndcap)) {
0678       if (cp[2 * is + 1].whichDetector() == DetId::Ecal &&
0679           (cp[2 * is + 1].whichSubDetector() == EcalBarrel || cp[2 * is + 1].whichSubDetector() == EcalEndcap)) {
0680         DetId cell2 = cp[2 * is + 1].getDetId();
0681         // set the real entrance
0682         if (ecalFirstSegment_ < 0)
0683           ecalFirstSegment_ = segments_.size();
0684 
0685         // !! Approximatiom : the first segment is always in a crystal
0686         if (cp[2 * is].getDetId() == cell2) {
0687           CaloSegment segment(cp[2 * is], cp[2 * is + 1], s, sX0, sL0, CaloSegment::PbWO4, myCalorimeter);
0688           segments_.push_back(segment);
0689           // std::cout << " Added (2)" << segment << std::endl;
0690           s += segment.length();
0691           sX0 += segment.X0length();
0692           sL0 += segment.L0length();
0693           X0ECAL_ += segment.X0length();
0694           L0ECAL_ += segment.L0length();
0695 #ifdef DEBUGCELLLINE
0696           ++ncrossedxtals;
0697 #endif
0698           ++is;
0699         } else {
0700           std::cout << " One more bug in the segment " << std::endl;
0701           ++is;
0702         }
0703         // Now check if a gap or crack should be added
0704         if (is > 0 && is < nsegments) {
0705           DetId cell3 = cp[2 * is].getDetId();
0706           if (cp[2 * is].whichDetector() != DetId::Hcal) {
0707             // Crack inside the ECAL
0708             bool bordercrossing = myCalorimeter->borderCrossing(cell2, cell3);
0709             CaloSegment cracksegment(cp[2 * is - 1],
0710                                      cp[2 * is],
0711                                      s,
0712                                      sX0,
0713                                      sL0,
0714                                      (bordercrossing) ? CaloSegment::CRACK : CaloSegment::GAP,
0715                                      myCalorimeter);
0716             segments_.push_back(cracksegment);
0717             s += cracksegment.length();
0718             sX0 += cracksegment.X0length();
0719             sL0 += cracksegment.L0length();
0720             X0ECAL_ += cracksegment.X0length();
0721             L0ECAL_ += cracksegment.L0length();
0722             //   std::cout <<" Added(3) "<< cracksegment << std::endl;
0723           } else {
0724             // a segment corresponding to ECAL/HCAL transition should be
0725             // added here
0726             CaloSegment cracksegment(cp[2 * is - 1], cp[2 * is], s, sX0, sL0, CaloSegment::ECALHCALGAP, myCalorimeter);
0727             segments_.push_back(cracksegment);
0728             s += cracksegment.length();
0729             sX0 += cracksegment.X0length();
0730             sL0 += cracksegment.L0length();
0731             X0EHGAP_ += cracksegment.X0length();
0732             L0EHGAP_ += cracksegment.L0length();
0733           }
0734         }
0735         continue;
0736       } else {
0737         std::cout << " Strange segment between " << cp[2 * is].whichDetector();
0738         std::cout << " and " << cp[2 * is + 1].whichDetector() << std::endl;
0739         ++is;
0740         continue;
0741       }
0742     }
0743     //      myHistos->debug("Just avant HCAL");
0744     // HCAL
0745     if (cp[2 * is].whichDetector() == DetId::Hcal && cp[2 * is + 1].whichDetector() == DetId::Hcal) {
0746       CaloSegment segment(cp[2 * is], cp[2 * is + 1], s, sX0, sL0, CaloSegment::HCAL, myCalorimeter);
0747       segments_.push_back(segment);
0748       s += segment.length();
0749       sX0 += segment.X0length();
0750       sL0 += segment.L0length();
0751       X0HCAL_ += segment.X0length();
0752       L0HCAL_ += segment.L0length();
0753       //      std::cout <<" Added(4) "<< segment << std::endl;
0754       ++is;
0755     }
0756   }
0757   //  std::cout << " PS1 " << X0PS1_ << " " << L0PS1_ << std::endl;
0758   //  std::cout << " PS2 " << X0PS2_ << " " << L0PS2_ << std::endl;
0759   //  std::cout << " ECAL " << X0ECAL_ << " " << L0ECAL_ << std::endl;
0760   //  std::cout << " HCAL " << X0HCAL_ << " " << L0HCAL_ << std::endl;
0761 
0762   totalX0_ = X0PS1_ + X0PS2_ + X0PS2EE_ + X0ECAL_ + X0EHGAP_ + X0HCAL_;
0763   totalL0_ = L0PS1_ + L0PS2_ + L0PS2EE_ + L0ECAL_ + L0EHGAP_ + L0HCAL_;
0764   //  myHistos->debug("Just avant le fill");
0765 
0766 #ifdef DEBUGCELLLINE
0767   myHistos->fill("h200", fabs(EcalEntrance_.eta()), X0ECAL_);
0768   myHistos->fill("h210", EcalEntrance_.phi(), X0ECAL_);
0769   if (X0ECAL_ < 20)
0770     myHistos->fill("h212", EcalEntrance_.phi(), X0ECAL_);
0771   //  if(X0ECAL_<1.)
0772   //    {
0773   //      for(unsigned ii=0; ii<segments_.size() ; ++ii)
0774   //    {
0775   //      std::cout << segments_[ii] << std::endl;
0776   //    }
0777   //    }
0778   myHistos->fillByNumber("h30", ncrossedxtals, EcalEntrance_.eta(), X0ECAL_);
0779 
0780   double zvertex = myTrack_->vertex().position().z();
0781 
0782   myHistos->fill("h310", EcalEntrance_.eta(), X0ECAL_);
0783   if (X0ECAL_ < 22)
0784     myHistos->fill("h410", EcalEntrance_.phi());
0785   myHistos->fill("h400", zvertex, X0ECAL_);
0786 #endif
0787   //  std::cout << " Finished the segments " << std::endl;
0788 }
0789 
0790 void EcalHitMaker::buildGeometry() {
0791   configuredGeometry_ = false;
0792   ncrystals_ = CellsWindow_.size();
0793   // create the vector with of pads with the appropriate size
0794   padsatdepth_.resize(ncrystals_);
0795 
0796   // This is fully correct in the barrel.
0797   ny_ = phisize_;
0798   nx_ = ncrystals_ / ny_;
0799   std::vector<unsigned> empty;
0800   empty.resize(ny_, 0);
0801   myCrystalNumberArray_.reserve((unsigned)nx_);
0802   for (unsigned inx = 0; inx < (unsigned)nx_; ++inx) {
0803     myCrystalNumberArray_.push_back(empty);
0804   }
0805 
0806   hits_.resize(ncrystals_, 0.);
0807   regionOfInterest_.clear();
0808   regionOfInterest_.resize(ncrystals_);
0809   validPads_.resize(ncrystals_);
0810   for (unsigned ic = 0; ic < ncrystals_; ++ic) {
0811     myCalorimeter->buildCrystal(CellsWindow_[ic], regionOfInterest_[ic]);
0812     regionOfInterest_[ic].setNumber(ic);
0813     DetIdMap_.insert(std::pair<DetId, unsigned>(CellsWindow_[ic], ic));
0814   }
0815 
0816   // Computes the map of the neighbours
0817   myCrystalWindowMap_ = new CrystalWindowMap(myCalorimeter, regionOfInterest_);
0818 }
0819 
0820 // depth is in X0 , L0 (depending on EMSHOWER/HADSHOWER) or in CM if inCM
0821 bool EcalHitMaker::getPads(double depth, bool inCm) {
0822   //std::cout << " New depth " << depth << std::endl;
0823   // The first time, the relationship between crystals must be calculated
0824   // but only in the case of EM showers
0825 
0826   if (EMSHOWER && !configuredGeometry_)
0827     configureGeometry();
0828 
0829   radiusFactor_ = (EMSHOWER) ? moliereRadius * radiusCorrectionFactor_ : interactionLength;
0830   detailedShowerTail_ = false;
0831   if (EMSHOWER)
0832     currentdepth_ = depth + X0depthoffset_;
0833   else
0834     currentdepth_ = depth;
0835 
0836   //  if(currentdepth_>maxX0_+ps1TotalX0()+ps2TotalX0())
0837   //    {
0838   //      currentdepth_=maxX0_+ps1TotalX0()+ps2TotalX0()-1.; // the -1 is for safety
0839   //      detailedShowerTail_=true;
0840   //    }
0841 
0842   //  std::cout << " FamosGrid::getQuads " << currentdepth_ << " " << maxX0_ << std::endl;
0843 
0844   ncrackpadsatdepth_ = 0;
0845 
0846   xmin_ = ymin_ = 999;
0847   xmax_ = ymax_ = -999;
0848   double locxmin, locxmax, locymin, locymax;
0849 
0850   // Get the depth of the pivot
0851   std::vector<CaloSegment>::const_iterator segiterator;
0852   // First identify the correct segment
0853 
0854   if (inCm)  // centimeter
0855   {
0856     segiterator = find_if(segments_.begin(), segments_.end(), CaloSegment::inSegment(currentdepth_));
0857   } else {
0858     // EM shower
0859     if (EMSHOWER)
0860       segiterator = find_if(segments_.begin(), segments_.end(), CaloSegment::inX0Segment(currentdepth_));
0861 
0862     //Hadron shower
0863     if (HADSHOWER)
0864       segiterator = find_if(segments_.begin(), segments_.end(), CaloSegment::inL0Segment(currentdepth_));
0865   }
0866   if (segiterator == segments_.end()) {
0867     std::cout << " FamosGrid: Could not go at such depth " << depth << std::endl;
0868     std::cout << " EMSHOWER " << EMSHOWER << std::endl;
0869     std::cout << " Track " << *myTrack_ << std::endl;
0870     std::cout << " Segments " << segments_.size() << std::endl;
0871     for (unsigned ii = 0; ii < segments_.size(); ++ii) {
0872       std::cout << segments_[ii] << std::endl;
0873     }
0874 
0875     return false;
0876   }
0877   //  std::cout << *segiterator << std::endl;
0878 
0879   if (segiterator->whichDetector() != DetId::Ecal) {
0880     std::cout << " In  " << segiterator->whichDetector() << std::endl;
0881     //      buildPreshower();
0882     return false;
0883   }
0884 
0885   //  std::cout << *segiterator << std::endl;
0886   // get the position of the origin
0887 
0888   XYZPoint origin;
0889   if (inCm) {
0890     origin = segiterator->positionAtDepthincm(currentdepth_);
0891   } else {
0892     if (EMSHOWER)
0893       origin = segiterator->positionAtDepthinX0(currentdepth_);
0894     if (HADSHOWER)
0895       origin = segiterator->positionAtDepthinL0(currentdepth_);
0896   }
0897   //  std::cout << " currentdepth_ " << currentdepth_ << " " << origin << std::endl;
0898   XYZVector newaxis = pivot_.getFirstEdge().Cross(normal_);
0899 
0900   //  std::cout  << "Normal " << normal_ << std::endl;
0901   //  std::cout << " New axis " << newaxis << std::endl;
0902 
0903   //  std::cout << " ncrystals  " << ncrystals << std::endl;
0904   plan_ = Plane3D((Vector)normal_, (Point)origin);
0905 
0906   unsigned nquads = 0;
0907   double sign = (central_) ? -1. : 1.;
0908   Transform3DR trans((Point)origin,
0909                      (Point)(origin + normal_),
0910                      (Point)(origin + newaxis),
0911                      Point(0, 0, 0),
0912                      Point(0., 0., sign),
0913                      Point(0., 1., 0.));
0914   for (unsigned ic = 0; ic < ncrystals_; ++ic) {
0915     //      std::cout << " Building geometry for " << regionOfInterest_[ic].getCellID() << std::endl;
0916     XYZPoint a, b;
0917 
0918     //      std::cout << " Origin " << origin << std::endl;
0919 
0920     //      std::vector<XYZPoint> corners;
0921     //      corners.reserve(4);
0922     double dummyt;
0923     bool hasbeenpulled = false;
0924     bool behindback = false;
0925     for (unsigned il = 0; il < 4; ++il) {
0926       // a is the il-th front corner of the crystal. b is the corresponding rear corner
0927       regionOfInterest_[ic].getLateralEdges(il, a, b);
0928 
0929       // pull the surface if necessary (only in the front of the crystals)
0930       XYZPoint aprime = a;
0931       if (pulled(origin, normal_, a)) {
0932         b = aprime;
0933         hasbeenpulled = true;
0934       }
0935 
0936       // compute the intersection.
0937       // Check that the intersection is in the [a,b] segment  if HADSHOWER
0938       // if EMSHOWER the intersection is calculated as if the crystals were infinite
0939       XYZPoint xx = (EMSHOWER) ? intersect(plan_, a, b, dummyt, false) : intersect(plan_, a, b, dummyt, true);
0940 
0941       if (dummyt > 1)
0942         behindback = true;
0943       //      std::cout << " Intersect " << il << " " << a << " " << b << " " << plan_ << " " << xx << std::endl;
0944       // check that the intersection actually exists
0945       if (xx.mag2() != 0) {
0946         corners[il] = xx;
0947       }
0948     }
0949     //      std::cout << " ncorners " << corners.size() << std::endl;
0950     if (behindback && EMSHOWER)
0951       detailedShowerTail_ = true;
0952     // If the quad is completly defined. Store it !
0953     if (corners.size() == 4) {
0954       padsatdepth_[ic] = CrystalPad(ic, corners, trans, bfactor_, !central_);
0955       // Parameter to be tuned
0956       if (hasbeenpulled)
0957         padsatdepth_[ic].setSurvivalProbability(pulledPadProbability_);
0958       validPads_[ic] = true;
0959       ++nquads;
0960       // In principle, this should be done after the quads reorganization. But it would cost one more loop
0961       //  quadsatdepth_[ic].extrems(locxmin,locxmax,locymin,locymax);
0962       //  if(locxmin<xmin_) xmin_=locxmin;
0963       //  if(locymin<ymin_) ymin_=locymin;
0964       //  if(locxmax>xmax_) xmax_=locxmax;
0965       //  if(locymax>ymax_) ymax_=locymax;
0966     } else {
0967       padsatdepth_[ic] = CrystalPad();
0968       validPads_[ic] = false;
0969     }
0970   }
0971   //  std::cout << " Number of quads " << quadsatdepth_.size() << std::endl;
0972   if (doreorg_)
0973     reorganizePads();
0974   //  std::cout << "Finished to reorganize " << std::endl;
0975   npadsatdepth_ = nquads;
0976   //  std::cout << " prepareCellIDMap " << std::endl;
0977 
0978   // Resize the Quads to allow for some numerical inaccuracy
0979   // in the "inside" function
0980   for (unsigned ic = 0; ic < ncrystals_; ++ic) {
0981     if (!validPads_[ic])
0982       continue;
0983 
0984     if (EMSHOWER)
0985       padsatdepth_[ic].resetCorners();
0986 
0987     padsatdepth_[ic].extrems(locxmin, locxmax, locymin, locymax);
0988     if (locxmin < xmin_)
0989       xmin_ = locxmin;
0990     if (locymin < ymin_)
0991       ymin_ = locymin;
0992     if (locxmax > xmax_)
0993       xmax_ = locxmax;
0994     if (locymax > ymax_)
0995       ymax_ = locymax;
0996   }
0997 
0998   sizex_ = (xmax_ - xmin_) / nx_;
0999   sizey_ = (ymax_ - ymin_) / ny_;
1000 
1001   // Make sure that sizex_ and sizey_ are set before running prepareCellIDMap
1002   prepareCrystalNumberArray();
1003 
1004   //  std::cout << " Finished  prepareCellIDMap " << std::endl;
1005   ncrackpadsatdepth_ = crackpadsatdepth_.size();
1006 
1007   return true;
1008 }
1009 
1010 void EcalHitMaker::configureGeometry() {
1011   configuredGeometry_ = true;
1012   for (unsigned ic = 0; ic < ncrystals_; ++ic) {
1013     //      std::cout << " Building " << cellids_[ic] << std::endl;
1014     for (unsigned idir = 0; idir < 8; ++idir) {
1015       unsigned oppdir = CaloDirectionOperations::oppositeDirection(idir);
1016       // Is there something else to do ?
1017       // The relationship with the neighbour may have been set previously.
1018       if (regionOfInterest_[ic].crystalNeighbour(idir).status() >= 0) {
1019         //        std::cout << " Nothing to do " << std::endl;
1020         continue;
1021       }
1022 
1023       const DetId& oldcell(regionOfInterest_[ic].getDetId());
1024       CaloDirection dir = CaloDirectionOperations::neighbourDirection(idir);
1025       DetId newcell(oldcell);
1026       if (!myCalorimeter->move(newcell, dir)) {
1027         // no neighbour in this direction
1028         regionOfInterest_[ic].crystalNeighbour(idir).setStatus(-1);
1029         continue;
1030       }
1031       // Determine the number of this neighbour
1032       //      std::cout << " The neighbour is " << newcell << std::endl;
1033       std::map<DetId, unsigned>::const_iterator niter(DetIdMap_.find(newcell));
1034       if (niter == DetIdMap_.end()) {
1035         //        std::cout << " The neighbour is not in the map " << std::endl;
1036         regionOfInterest_[ic].crystalNeighbour(idir).setStatus(-1);
1037         continue;
1038       }
1039       // Now there is a neighbour
1040       //      std::cout << " The neighbour is " << niter->second << " " << cellids_[niter->second] << std::endl;
1041       regionOfInterest_[ic].crystalNeighbour(idir).setNumber(niter->second);
1042       //      std::cout << " Managed to set crystalNeighbour " << ic << " " << idir << std::endl;
1043       //      std::cout << " Trying  " << niter->second << " " << oppdir << std::endl;
1044       regionOfInterest_[niter->second].crystalNeighbour(oppdir).setNumber(ic);
1045       //      std::cout << " Crack/gap " << std::endl;
1046       if (myCalorimeter->borderCrossing(oldcell, newcell)) {
1047         regionOfInterest_[ic].crystalNeighbour(idir).setStatus(1);
1048         regionOfInterest_[niter->second].crystalNeighbour(oppdir).setStatus(1);
1049         //        std::cout << " Crack ! " << std::endl;
1050       } else {
1051         regionOfInterest_[ic].crystalNeighbour(idir).setStatus(0);
1052         regionOfInterest_[niter->second].crystalNeighbour(oppdir).setStatus(0);
1053         //        std::cout << " Gap" << std::endl;
1054       }
1055     }
1056   }
1057   // Magnetic field a la Charlot
1058   double theta = EcalEntrance_.theta();
1059   if (theta > M_PI_2)
1060     theta = M_PI - theta;
1061   bfactor_ = 1. / (1. + 0.133 * theta);
1062   if (myCalorimeter->magneticField() == 0.)
1063     bfactor_ = 1.;
1064 }
1065 
1066 // project fPoint on the plane (original,normal)
1067 bool EcalHitMaker::pulled(const XYZPoint& origin, const XYZNormal& normal, XYZPoint& fPoint) const {
1068   // check if fPoint is behind the origin
1069   double dotproduct = normal.Dot(fPoint - origin);
1070   if (dotproduct <= 0.)
1071     return false;
1072 
1073   //norm of normal is 1
1074   fPoint -= (1 + dotproduct) * normal;
1075   return true;
1076 }
1077 
1078 void EcalHitMaker::prepareCrystalNumberArray() {
1079   for (unsigned iq = 0; iq < npadsatdepth_; ++iq) {
1080     if (!validPads_[iq])
1081       continue;
1082     unsigned d1, d2;
1083     convertIntegerCoordinates(padsatdepth_[iq].center().x(), padsatdepth_[iq].center().y(), d1, d2);
1084     myCrystalNumberArray_[d1][d2] = iq;
1085   }
1086 }
1087 
1088 void EcalHitMaker::convertIntegerCoordinates(double x, double y, unsigned& ix, unsigned& iy) const {
1089   int tix = (int)((x - xmin_) / sizex_);
1090   int tiy = (int)((y - ymin_) / sizey_);
1091   ix = iy = 9999;
1092   if (tix >= 0)
1093     ix = (unsigned)tix;
1094   if (tiy >= 0)
1095     iy = (unsigned)tiy;
1096 }
1097 
1098 const std::map<CaloHitID, float>& EcalHitMaker::getHits() {
1099   if (hitmaphasbeencalculated_)
1100     return hitMap_;
1101   for (unsigned ic = 0; ic < ncrystals_; ++ic) {
1102     //calculate time of flight
1103     float tof = 0.0;
1104     if (onEcal_ == 1 || onEcal_ == 2)
1105       tof =
1106           (myCalorimeter->getEcalGeometry(onEcal_)->getGeometry(regionOfInterest_[ic].getDetId())->getPosition().mag()) /
1107           29.98;  //speed of light
1108 
1109     if (onEcal_ == 1) {
1110       CaloHitID current_id(EBDetId(regionOfInterest_[ic].getDetId().rawId()).hashedIndex(), tof, 0);  //no track yet
1111       hitMap_.insert(std::pair<CaloHitID, float>(current_id, hits_[ic]));
1112     } else if (onEcal_ == 2) {
1113       CaloHitID current_id(EEDetId(regionOfInterest_[ic].getDetId().rawId()).hashedIndex(), tof, 0);  //no track yet
1114       hitMap_.insert(std::pair<CaloHitID, float>(current_id, hits_[ic]));
1115     }
1116   }
1117   hitmaphasbeencalculated_ = true;
1118   return hitMap_;
1119 }
1120 
1121 ///////////////////////////// GAPS/CRACKS TREATMENT////////////
1122 
1123 void EcalHitMaker::reorganizePads() {
1124   // Some cleaning first
1125   //  std::cout << " Starting reorganize " << std::endl;
1126   crackpadsatdepth_.clear();
1127   crackpadsatdepth_.reserve(etasize_ * phisize_);
1128   ncrackpadsatdepth_ = 0;
1129   std::vector<neighbour> gaps;
1130   std::vector<std::vector<neighbour> > cracks;
1131   //  cracks.clear();
1132   cracks.resize(ncrystals_);
1133 
1134   for (unsigned iq = 0; iq < ncrystals_; ++iq) {
1135     if (!validPads_[iq])
1136       continue;
1137 
1138     gaps.clear();
1139     //   std::cout << padsatdepth_[iq] << std::endl;
1140     //check all the directions
1141     for (unsigned iside = 0; iside < 4; ++iside) {
1142       //      std::cout << " To be glued " << iside << " " << regionOfInterest_[iq].crystalNeighbour(iside).toBeGlued() << std::endl;
1143       CaloDirection thisside = CaloDirectionOperations::Side(iside);
1144       if (regionOfInterest_[iq].crystalNeighbour(iside).toBeGlued()) {
1145         // look for the neighbour and check that it exists
1146         int neighbourstatus = regionOfInterest_[iq].crystalNeighbour(iside).status();
1147         if (neighbourstatus < 0)
1148           continue;
1149 
1150         unsigned neighbourNumber = regionOfInterest_[iq].crystalNeighbour(iside).number();
1151         if (!validPads_[neighbourNumber])
1152           continue;
1153         // there is a crack between
1154         if (neighbourstatus == 1) {
1155           //          std::cout << " 1 Crack : " << thisside << " " << cellids_[iq]<< " " << cellids_[neighbourNumber] << std::endl;
1156           cracks[iq].push_back(neighbour(thisside, neighbourNumber));
1157         }  // else it is a gap
1158         else {
1159           gaps.push_back(neighbour(thisside, neighbourNumber));
1160         }
1161       }
1162     }
1163     // Now lift the gaps
1164     gapsLifting(gaps, iq);
1165   }
1166 
1167   unsigned ncracks = cracks.size();
1168   //  std::cout << " Cracks treatment : " << cracks.size() << std::endl;
1169   for (unsigned icrack = 0; icrack < ncracks; ++icrack) {
1170     //      std::cout << " Crack number " << crackiter->first << std::endl;
1171     cracksPads(cracks[icrack], icrack);
1172   }
1173 }
1174 
1175 //dir 2 = N,E,W,S
1176 CLHEP::Hep2Vector& EcalHitMaker::correspondingEdge(neighbour& myneighbour, CaloDirection dir2) {
1177   CaloDirection dir = CaloDirectionOperations::oppositeSide(myneighbour.first);
1178   CaloDirection corner = CaloDirectionOperations::add2d(dir, dir2);
1179   //  std::cout << "Corresponding Edge " << dir<< " " << dir2 << " " << corner << std::endl;
1180   return padsatdepth_[myneighbour.second].edge(corner);
1181 }
1182 
1183 bool EcalHitMaker::diagonalEdge(unsigned myPad, CaloDirection dir, CLHEP::Hep2Vector& point) {
1184   unsigned idir = CaloDirectionOperations::neighbourDirection(dir);
1185   if (regionOfInterest_[myPad].crystalNeighbour(idir).status() < 0)
1186     return false;
1187   unsigned nneighbour = regionOfInterest_[myPad].crystalNeighbour(idir).number();
1188   if (!validPads_[nneighbour]) {
1189     //      std::cout << " Wasn't able to move " << std::endl;
1190     return false;
1191   }
1192   point = padsatdepth_[nneighbour].edge(CaloDirectionOperations::oppositeSide(dir));
1193   return true;
1194 }
1195 
1196 bool EcalHitMaker::unbalancedDirection(const std::vector<neighbour>& dirs,
1197                                        unsigned& unb,
1198                                        unsigned& dir1,
1199                                        unsigned& dir2) {
1200   if (dirs.size() == 1)
1201     return false;
1202   if (dirs.size() % 2 == 0)
1203     return false;
1204   CaloDirection tmp;
1205   tmp = CaloDirectionOperations::add2d(dirs[0].first, dirs[1].first);
1206   if (tmp == NONE) {
1207     unb = 2;
1208     dir1 = 0;
1209     dir2 = 1;
1210     return true;
1211   }
1212   tmp = CaloDirectionOperations::add2d(dirs[0].first, dirs[2].first);
1213   if (tmp == NONE) {
1214     unb = 1;
1215     dir1 = 0;
1216     dir2 = 2;
1217     return true;
1218   }
1219   unb = 0;
1220   dir1 = 1;
1221   dir2 = 2;
1222   return true;
1223 }
1224 
1225 void EcalHitMaker::gapsLifting(std::vector<neighbour>& gaps, unsigned iq) {
1226   //  std::cout << " Entering gapsLifting "  << std::endl;
1227   CrystalPad& myPad = padsatdepth_[iq];
1228   unsigned ngaps = gaps.size();
1229   constexpr bool debug = false;
1230   if (ngaps == 1) {
1231     if (debug) {
1232       std::cout << " Avant " << ngaps << " " << gaps[0].first << std::endl;
1233       std::cout << myPad << std::endl;
1234     }
1235     if (gaps[0].first == NORTH || gaps[0].first == SOUTH) {
1236       CaloDirection dir1 = CaloDirectionOperations::add2d(gaps[0].first, EAST);
1237       CaloDirection dir2 = CaloDirectionOperations::add2d(gaps[0].first, WEST);
1238       myPad.edge(dir1) = correspondingEdge(gaps[0], EAST);
1239       myPad.edge(dir2) = correspondingEdge(gaps[0], WEST);
1240     } else {
1241       CaloDirection dir1 = CaloDirectionOperations::add2d(gaps[0].first, NORTH);
1242       CaloDirection dir2 = CaloDirectionOperations::add2d(gaps[0].first, SOUTH);
1243       myPad.edge(dir1) = correspondingEdge(gaps[0], NORTH);
1244       myPad.edge(dir2) = correspondingEdge(gaps[0], SOUTH);
1245     }
1246     if (debug) {
1247       std::cout << " Apres " << std::endl;
1248       std::cout << myPad << std::endl;
1249     }
1250   } else if (ngaps == 2) {
1251     if (debug) {
1252       std::cout << " Avant " << ngaps << " " << gaps[0].first << " " << gaps[1].first << std::endl;
1253       std::cout << myPad << std::endl;
1254       std::cout << " Voisin 1 " << (gaps[0].second) << std::endl;
1255       std::cout << " Voisin 2 " << (gaps[1].second) << std::endl;
1256     }
1257     CaloDirection corner0 = CaloDirectionOperations::add2d(gaps[0].first, gaps[1].first);
1258 
1259     CLHEP::Hep2Vector point(0., 0.);
1260     if (corner0 != NONE && diagonalEdge(iq, corner0, point)) {
1261       CaloDirection corner1 =
1262           CaloDirectionOperations::add2d(CaloDirectionOperations::oppositeSide(gaps[0].first), gaps[1].first);
1263       CaloDirection corner2 =
1264           CaloDirectionOperations::add2d(gaps[0].first, CaloDirectionOperations::oppositeSide(gaps[1].first));
1265       myPad.edge(corner0) = point;
1266       myPad.edge(corner1) = correspondingEdge(gaps[1], CaloDirectionOperations::oppositeSide(gaps[0].first));
1267       myPad.edge(corner2) = correspondingEdge(gaps[0], CaloDirectionOperations::oppositeSide(gaps[1].first));
1268     } else if (corner0 == NONE) {
1269       if (gaps[0].first == EAST || gaps[0].first == WEST) {
1270         CaloDirection corner1 = CaloDirectionOperations::add2d(gaps[0].first, NORTH);
1271         CaloDirection corner2 = CaloDirectionOperations::add2d(gaps[0].first, SOUTH);
1272         myPad.edge(corner1) = correspondingEdge(gaps[0], NORTH);
1273         myPad.edge(corner2) = correspondingEdge(gaps[0], SOUTH);
1274 
1275         corner1 = CaloDirectionOperations::add2d(gaps[1].first, NORTH);
1276         corner2 = CaloDirectionOperations::add2d(gaps[1].first, SOUTH);
1277         myPad.edge(corner1) = correspondingEdge(gaps[1], NORTH);
1278         myPad.edge(corner2) = correspondingEdge(gaps[1], SOUTH);
1279       } else {
1280         CaloDirection corner1 = CaloDirectionOperations::add2d(gaps[0].first, EAST);
1281         CaloDirection corner2 = CaloDirectionOperations::add2d(gaps[0].first, WEST);
1282         myPad.edge(corner1) = correspondingEdge(gaps[0], EAST);
1283         myPad.edge(corner2) = correspondingEdge(gaps[0], WEST);
1284 
1285         corner1 = CaloDirectionOperations::add2d(gaps[1].first, EAST);
1286         corner2 = CaloDirectionOperations::add2d(gaps[1].first, WEST);
1287         myPad.edge(corner1) = correspondingEdge(gaps[1], EAST);
1288         myPad.edge(corner2) = correspondingEdge(gaps[1], WEST);
1289       }
1290     }
1291     if (debug) {
1292       std::cout << " Apres " << std::endl;
1293       std::cout << myPad << std::endl;
1294     }
1295   } else if (ngaps == 3) {
1296     // in this case the four corners have to be changed
1297     unsigned iubd, idir1, idir2;
1298     CaloDirection diag;
1299     CLHEP::Hep2Vector point(0., 0.);
1300     //    std::cout << " Yes : 3 gaps" << std::endl;
1301     if (unbalancedDirection(gaps, iubd, idir1, idir2)) {
1302       CaloDirection ubd(gaps[iubd].first), dir1(gaps[idir1].first);
1303       CaloDirection dir2(gaps[idir2].first);
1304 
1305       //          std::cout << " Avant " << std::endl << myPad << std::endl;
1306       //          std::cout << ubd << " " << dir1 << " " << dir2 << std::endl;
1307       diag = CaloDirectionOperations::add2d(ubd, dir1);
1308       if (diagonalEdge(iq, diag, point))
1309         myPad.edge(diag) = point;
1310       diag = CaloDirectionOperations::add2d(ubd, dir2);
1311       if (diagonalEdge(iq, diag, point))
1312         myPad.edge(diag) = point;
1313       CaloDirection oppside = CaloDirectionOperations::oppositeSide(ubd);
1314       myPad.edge(CaloDirectionOperations::add2d(oppside, dir1)) = correspondingEdge(gaps[idir1], oppside);
1315       myPad.edge(CaloDirectionOperations::add2d(oppside, dir2)) = correspondingEdge(gaps[idir2], oppside);
1316       //          std::cout << " Apres " << std::endl << myPad << std::endl;
1317     }
1318   } else if (ngaps == 4) {
1319     //      std::cout << " Waouh :4 gaps" << std::endl;
1320     //      std::cout << " Avant " << std::endl;
1321     //      std::cout << myPad<< std::endl;
1322     CLHEP::Hep2Vector point(0., 0.);
1323     if (diagonalEdge(iq, NORTHEAST, point))
1324       myPad.edge(NORTHEAST) = point;
1325     if (diagonalEdge(iq, NORTHWEST, point))
1326       myPad.edge(NORTHWEST) = point;
1327     if (diagonalEdge(iq, SOUTHWEST, point))
1328       myPad.edge(SOUTHWEST) = point;
1329     if (diagonalEdge(iq, SOUTHEAST, point))
1330       myPad.edge(SOUTHEAST) = point;
1331     //      std::cout << " Apres " << std::endl;
1332     //      std::cout << myPad<< std::endl;
1333   }
1334 }
1335 
1336 void EcalHitMaker::cracksPads(std::vector<neighbour>& cracks, unsigned iq) {
1337   //  std::cout << " myPad " << &myPad << std::endl;
1338   unsigned ncracks = cracks.size();
1339   CrystalPad& myPad = padsatdepth_[iq];
1340   for (unsigned ic = 0; ic < ncracks; ++ic) {
1341     //      std::vector<CLHEP::Hep2Vector> mycorners;
1342     //      mycorners.reserve(4);
1343     switch (cracks[ic].first) {
1344       case NORTH: {
1345         mycorners[0] = (padsatdepth_[cracks[ic].second].edge(SOUTHWEST));
1346         mycorners[1] = (padsatdepth_[cracks[ic].second].edge(SOUTHEAST));
1347         mycorners[2] = (myPad.edge(NORTHEAST));
1348         mycorners[3] = (myPad.edge(NORTHWEST));
1349       } break;
1350       case SOUTH: {
1351         mycorners[0] = (myPad.edge(SOUTHWEST));
1352         mycorners[1] = (myPad.edge(SOUTHEAST));
1353         mycorners[2] = (padsatdepth_[cracks[ic].second].edge(NORTHEAST));
1354         mycorners[3] = (padsatdepth_[cracks[ic].second].edge(NORTHWEST));
1355       } break;
1356       case EAST: {
1357         mycorners[0] = (myPad.edge(NORTHEAST));
1358         mycorners[1] = (padsatdepth_[cracks[ic].second].edge(NORTHWEST));
1359         mycorners[2] = (padsatdepth_[cracks[ic].second].edge(SOUTHWEST));
1360         mycorners[3] = (myPad.edge(SOUTHEAST));
1361       } break;
1362       case WEST: {
1363         mycorners[0] = (padsatdepth_[cracks[ic].second].edge(NORTHEAST));
1364         mycorners[1] = (myPad.edge(NORTHWEST));
1365         mycorners[2] = (myPad.edge(SOUTHWEST));
1366         mycorners[3] = (padsatdepth_[cracks[ic].second].edge(SOUTHEAST));
1367       } break;
1368       default: {
1369       }
1370     }
1371     CrystalPad crackpad(ic, mycorners);
1372     // to be tuned. A simpleconfigurable should be used
1373     crackpad.setSurvivalProbability(crackPadProbability_);
1374     crackpadsatdepth_.push_back(crackpad);
1375   }
1376   //  std::cout << " Finished cracksPads " << std::endl;
1377 }
1378 
1379 bool EcalHitMaker::inside3D(const std::vector<XYZPoint>& corners, const XYZPoint& p) const {
1380   // corners and p are in the same plane
1381   // p is inside "corners" if the four crossproducts (corners[i]xcorners[i+1])
1382   // are in the same direction
1383 
1384   XYZVector crossproduct(0., 0., 0.), previouscrossproduct(0., 0., 0.);
1385 
1386   for (unsigned ip = 0; ip < 4; ++ip) {
1387     crossproduct = (corners[ip] - p).Cross(corners[(ip + 1) % 4] - p);
1388     if (ip == 0)
1389       previouscrossproduct = crossproduct;
1390     else if (crossproduct.Dot(previouscrossproduct) < 0.)
1391       return false;
1392   }
1393 
1394   return true;
1395 }
[ 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!