CaloHitMakers/src/HcalHitMaker.cc

0001 #include "FastSimulation/CaloHitMakers/interface/HcalHitMaker.h"
0002 #include "FastSimulation/CaloGeometryTools/interface/CaloGeometryHelper.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "DataFormats/HcalDetId/interface/HcalDetId.h"  // PV
0005 #include <algorithm>
0006 #include <cmath>
0007
0008 typedef ROOT::Math::Transform3DPJ::Point Point;
0009
0010 HcalHitMaker::HcalHitMaker(EcalHitMaker& grid, unsigned shower)
0011     : CaloHitMaker(grid.getCalorimeter(),
0012                    DetId::Hcal,
0013                    HcalHitMaker::getSubHcalDet(grid.getFSimTrack()),
0014                    grid.getFSimTrack()->onHcal() ? grid.getFSimTrack()->onHcal() : grid.getFSimTrack()->onVFcal() + 1,
0015                    shower),
0016       myGrid(grid),
0017       myTrack((grid.getFSimTrack())) {
0018   // normalize the direction
0019   ecalEntrance_ = myGrid.ecalEntrance();
0020   particleDirection = myTrack->ecalEntrance().Vect().Unit();
0021   radiusFactor_ = (EMSHOWER) ? moliereRadius : interactionLength;
0022   mapCalculated_ = false;
0023   //std::cout << " Famos HCAL " << grid.getTrack()->onHcal() << " " <<  grid.getTrack()->onVFcal() << " " << showerType << std::endl;
0024   if (EMSHOWER && (abs(grid.getFSimTrack()->type()) != 11 && grid.getFSimTrack()->type() != 22)) {
0025     std::cout << " FamosHcalHitMaker : Strange. The following shower has EM type" << std::endl
0026               << *grid.getFSimTrack() << std::endl;
0027   }
0028 }
0029
0030 bool HcalHitMaker::addHit(double r, double phi, unsigned layer) {
0031   //  std::cout << " FamosHcalHitMaker::addHit - radiusFactor = " << radiusFactor
0032   //        << std::endl;
0033
0034   XYZPoint point(r * radiusFactor_ * std::cos(phi), r * radiusFactor_ * std::sin(phi), 0.);
0035
0036   // Watch out !!!! (Point) is a real point in the MathCore terminology (not a redefined a XYZPoint which
0037   // is actually a XYZVector in the MatchCore terminology). Therefore, the Transform3D is correctly applied
0038   point = locToGlobal_((Point)point);
0039   return addHit(point, layer);
0040 }
0041
0042 bool HcalHitMaker::addHit(const XYZPoint& point, unsigned layer) {
0043   // Temporary nasty hacks to avoid misbehaviour of not-intended-for-that
0044   //  getClosestCell in case of large (eta beyond HF ...)  and in EM showers
0045   if (fabs(point.Z()) > 2000 || fabs(point.X()) > 2000 || fabs(point.Y()) > 2000) {
0046     if (EMSHOWER)
0047       edm::LogWarning("HcalHitMaker") << "received a hit very far from the detector " << point
0048                                       << " coming from an electromagnetic shower. - Ignoring it" << std::endl;
0049     else if (HADSHOWER)
0050       edm::LogWarning("HcalHitMaker") << "received a hit very far from the detector " << point
0051                                       << " coming from a hadron shower. - Ignoring it" << std::endl;
0052     else if (MIP)
0053       edm::LogWarning("HcalHitMaker") << "received a hit very far from the detector " << point
0054                                       << " coming from a muon. - Ignoring it" << std::endl;
0055     return false;
0056   }
0057
0058   double pointeta = fabs(point.eta());
0059   if (pointeta > 5.19)
0060     return false;
0061
0062   //calculate time of flight
0063   double dist = std::sqrt(point.X() * point.X() + point.Y() * point.Y() + point.Z() * point.Z());
0064   double tof = dist / 29.98;  //speed of light
0065
0066   DetId thecellID(myCalorimeter->getClosestCell(point, false, false));
0067
0068   HcalDetId myDetId(thecellID);
0069
0070   //   if ( myDetId.subdetId() == HcalForward ) {
0071   //     std::cout << "HcalHitMaker : " << point.Z() << " " << myDetId.depth()    << std::endl;
0072   //   }
0073
0074   //   std::cout << "BEFORE" << std::endl;
0075   //   std::cout << "HcalHitMaker : subdetId : " << myDetId.subdetId() << std::endl;
0076   //   std::cout << "HcalHitMaker : depth    : " << myDetId.depth()    << std::endl;
0077   //   std::cout << "HcalHitMaker : ieta     : " << myDetId.ieta()     << std::endl;
0078   //   std::cout << "HcalHitMaker : iphi     : " << myDetId.iphi()     << std::endl;
0079   //   std::cout << "HcalHitMaker : spotE    : " << spotEnergy         << std::endl;
0080   //   std::cout << "HcalHitMaker : point.X  : " << point.X()          << std::endl;
0081   //   std::cout << "HcalHitMaker : point.Y  : " << point.Y()          << std::endl;
0082   //   std::cout << "HcalHitMaker : point.Z  : " << point.Z()          << std::endl;
0083
0084   if (myDetId.subdetId() == HcalForward) {
0085     int mylayer = layer;
0086     if (myDetId.depth() == 2) {
0087       mylayer = (int)layer;
0088     } else {
0089       mylayer = 1;
0090     }
0091     HcalDetId myDetId2((HcalSubdetector)myDetId.subdetId(), myDetId.ieta(), myDetId.iphi(), mylayer);
0092     thecellID = myDetId2;
0093     myDetId = myDetId2;
0094   }
0095
0096   if (!thecellID.null() && myDetId.depth() > 0) {
0097     CaloHitID current_id(thecellID.rawId(), tof, 0);  //no track yet
0098
0099     //      std::cout << " FamosHcalHitMaker::addHit - the cell num " << cell
0100     //              << std::endl;
0101
0102     std::map<CaloHitID, float>::iterator cellitr;
0103     cellitr = hitMap_.find(current_id);
0104     if (cellitr == hitMap_.end()) {
0105       hitMap_.insert(std::pair<CaloHitID, float>(current_id, spotEnergy));
0106     } else {
0107       cellitr->second += spotEnergy;
0108     }
0109     return true;
0110   }
0111   return false;
0112 }
0113
0114 bool HcalHitMaker::setDepth(double depth, bool inCm) {
0115   currentDepth_ = depth;
0116   std::vector<CaloSegment>::const_iterator segiterator;
0117   if (inCm) {
0118     segiterator =
0119         find_if(myGrid.getSegments().begin(), myGrid.getSegments().end(), CaloSegment::inSegment(currentDepth_));
0120   } else {
0121     if (EMSHOWER)
0122       segiterator =
0123           find_if(myGrid.getSegments().begin(), myGrid.getSegments().end(), CaloSegment::inX0Segment(currentDepth_));
0124
0125     //Hadron shower
0126     if (HADSHOWER)
0127       segiterator =
0128           find_if(myGrid.getSegments().begin(), myGrid.getSegments().end(), CaloSegment::inL0Segment(currentDepth_));
0129   }
0130
0131   if (segiterator == myGrid.getSegments().end()) {
0132     // Special trick  - As advised by Salavat, no leakage should be simulated
0133     if (depth > myGrid.getSegments().back().sL0Exit()) {
0134       segiterator = find_if(myGrid.getSegments().begin(),
0135                             myGrid.getSegments().end(),
0136                             CaloSegment::inL0Segment(myGrid.getSegments().back().sL0Exit() - 1.));
0137       depth = segiterator->sL0Exit() - 1.;
0138       currentDepth_ = depth;
0139       if (segiterator == myGrid.getSegments().end()) {
0140         std::cout << " Could not go at such depth " << EMSHOWER << "  " << currentDepth_ << std::endl;
0141         std::cout << " Track " << *(myGrid.getFSimTrack()) << std::endl;
0142         return false;
0143       }
0144     } else {
0145       std::cout << " Could not go at such depth " << EMSHOWER << "  " << currentDepth_ << " "
0146                 << myGrid.getSegments().back().sL0Exit() << std::endl;
0147       std::cout << " Track " << *(myGrid.getFSimTrack()) << std::endl;
0148       return false;
0149     }
0150   }
0151
0152   XYZPoint origin;
0153   if (inCm) {
0154     origin = segiterator->positionAtDepthincm(currentDepth_);
0155   } else {
0156     if (EMSHOWER)
0157       origin = segiterator->positionAtDepthinX0(currentDepth_);
0158     if (HADSHOWER)
0159       origin = segiterator->positionAtDepthinL0(currentDepth_);
0160   }
0161   XYZVector zaxis(0, 0, 1);
0162   XYZVector planeVec1 = (zaxis.Cross(particleDirection)).Unit();
0163   locToGlobal_ = Transform3D(Point(0, 0, 0),
0164                              Point(0, 0, 1),
0165                              Point(1, 0, 0),
0166                              (Point)origin,
0167                              (Point)(origin + particleDirection),
0168                              (Point)(origin + planeVec1));
0169   return true;
0170 }