HcalHitMaker.cc

Go to the documentation of this file.

#include "FastSimulation/CaloHitMakers/interface/HcalHitMaker.h"
#include "FastSimulation/CaloGeometryTools/interface/CaloGeometryHelper.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h" // PV
#include <algorithm>
#include <cmath>

typedef ROOT::Math::Transform3DPJ::Point Point;

HcalHitMaker::HcalHitMaker(EcalHitMaker& grid,unsigned shower)
  :CaloHitMaker(grid.getCalorimeter(),DetId::Hcal,HcalHitMaker::getSubHcalDet(grid.getFSimTrack()),
        grid.getFSimTrack()->onHcal()?grid.getFSimTrack()->onHcal():grid.getFSimTrack()->onVFcal()+1,shower),
   myGrid(grid),  myTrack((grid.getFSimTrack()))
{
  // normalize the direction
  ecalEntrance_=myGrid.ecalEntrance();
  particleDirection=myTrack->ecalEntrance().Vect().Unit();
  radiusFactor_=(EMSHOWER)? moliereRadius:interactionLength;
  mapCalculated_=false;
  //std::cout << " Famos HCAL " << grid.getTrack()->onHcal() << " " <<  grid.getTrack()->onVFcal() << " " << showerType << std::endl;
  if(EMSHOWER&&(abs(grid.getFSimTrack()->type())!=11 && grid.getFSimTrack()->type()!=22))
    {
      std::cout << " FamosHcalHitMaker : Strange. The following shower has EM type" << std::endl <<* grid.getFSimTrack() << std::endl;
    }
}

bool 
HcalHitMaker::addHit(double r,double phi,unsigned layer)
{
    //  std::cout << " FamosHcalHitMaker::addHit - radiusFactor = " << radiusFactor
  //        << std::endl;

  XYZPoint point(r*radiusFactor_*std::cos(phi),r*radiusFactor_*std::sin(phi),0.);

  // Watch out !!!! (Point) is a real point in the MathCore terminology (not a redefined a XYZPoint which
  // is actually a XYZVector in the MatchCore terminology). Therefore, the Transform3D is correctly applied
  point = locToGlobal_((Point)point);
  return addHit(point,layer);
}

bool HcalHitMaker::addHit(const XYZPoint& point, unsigned layer)
{
  // Temporary nasty hacks to avoid misbehaviour of not-intended-for-that
  //  getClosestCell in case of large (eta beyond HF ...)  and in EM showers 
  if(fabs(point.Z())>2000 || fabs(point.X())>2000 || fabs(point.Y())>2000) 
    { 
      if(EMSHOWER) 
    edm::LogWarning("HcalHitMaker") << "received a hit very far from the detector " << point << " coming from an electromagnetic shower. - Ignoring it" << std::endl;
      else if(HADSHOWER)
    edm::LogWarning("HcalHitMaker") << "received a hit very far from the detector " << point << " coming from a hadron shower. - Ignoring it" << std::endl;
      else if(MIP)
    edm::LogWarning("HcalHitMaker") << "received a hit very far from the detector " << point << " coming from a muon. - Ignoring it" << std::endl;
      return false; 
    } 


  double pointeta = fabs(point.eta());
  if(pointeta > 5.19) return false; 

  //calculate time of flight
  double dist = std::sqrt(point.X()*point.X() + point.Y()*point.Y() + point.Z()*point.Z());
  double tof = dist/29.98; //speed of light

  DetId thecellID(myCalorimeter->getClosestCell(point,false,false));
  
  HcalDetId myDetId(thecellID);
                                                                                                                                      
//   if ( myDetId.subdetId() == HcalForward ) {
//     std::cout << "HcalHitMaker : " << point.Z() << " " << myDetId.depth()    << std::endl;
//   }
                                                                                                                                      
//   std::cout << "BEFORE" << std::endl;
//   std::cout << "HcalHitMaker : subdetId : " << myDetId.subdetId() << std::endl;
//   std::cout << "HcalHitMaker : depth    : " << myDetId.depth()    << std::endl;
//   std::cout << "HcalHitMaker : ieta     : " << myDetId.ieta()     << std::endl;
//   std::cout << "HcalHitMaker : iphi     : " << myDetId.iphi()     << std::endl;
//   std::cout << "HcalHitMaker : spotE    : " << spotEnergy         << std::endl;
//   std::cout << "HcalHitMaker : point.X  : " << point.X()          << std::endl;
//   std::cout << "HcalHitMaker : point.Y  : " << point.Y()          << std::endl;
//   std::cout << "HcalHitMaker : point.Z  : " << point.Z()          << std::endl;
                                                                                                                                      
  if ( myDetId.subdetId() == HcalForward ) {
    int mylayer = layer;
    if ( myDetId.depth()==2 ) {
      mylayer = (int)layer;
    } else {
      mylayer = 1;
    }
    HcalDetId myDetId2((HcalSubdetector)myDetId.subdetId(),myDetId.ieta(),myDetId.iphi(),mylayer);
    thecellID = myDetId2;
    myDetId = myDetId2;
  }


  
  if(!thecellID.null()  && myDetId.depth()>0)
    {   
      CaloHitID current_id(thecellID.rawId(),tof,0); //no track yet
      
      //      std::cout << " FamosHcalHitMaker::addHit - the cell num " << cell
      //            << std::endl;

      std::map<CaloHitID,float>::iterator cellitr;
      cellitr = hitMap_.find(current_id);
      if(cellitr==hitMap_.end())
    {
      hitMap_.insert(std::pair<CaloHitID,float>(current_id,spotEnergy));
    }
      else
    {
      cellitr->second+=spotEnergy;
    }
      return true;
    }
  return false;
}

bool 
HcalHitMaker::setDepth(double depth,bool inCm)
{
  currentDepth_=depth;
  std::vector<CaloSegment>::const_iterator segiterator;
  if(inCm)
    {
      segiterator = find_if(myGrid.getSegments().begin(),myGrid.getSegments().end(),CaloSegment::inSegment(currentDepth_));
    }
  else
    {
      if(EMSHOWER)
    segiterator = find_if(myGrid.getSegments().begin(),myGrid.getSegments().end(),CaloSegment::inX0Segment(currentDepth_));
      
      //Hadron shower 
      if(HADSHOWER)
    segiterator = find_if(myGrid.getSegments().begin(),myGrid.getSegments().end(),CaloSegment::inL0Segment(currentDepth_));
    }
  
  if(segiterator==myGrid.getSegments().end()) 
    {
      // Special trick  - As advised by Salavat, no leakage should be simulated
      if(depth > myGrid.getSegments().back().sL0Exit())
    {
      segiterator= find_if(myGrid.getSegments().begin(),myGrid.getSegments().end(),CaloSegment::inL0Segment(myGrid.getSegments().back().sL0Exit()-1.));
      depth=segiterator->sL0Exit()-1.;
      currentDepth_=depth;
      if(segiterator==myGrid.getSegments().end())
        {
          std::cout << " Could not go at such depth " << EMSHOWER << "  " << currentDepth_ << std::endl;
          std::cout << " Track " << *(myGrid.getFSimTrack()) << std::endl;
          return false;
        }
    }
      else
    {
      std::cout << " Could not go at such depth " << EMSHOWER << "  " << currentDepth_ << " " << myGrid.getSegments().back().sL0Exit() << std::endl; 
      std::cout << " Track " << *(myGrid.getFSimTrack()) << std::endl; 
      return false; 
    }
    }


  XYZPoint origin;
  if(inCm)
    {
      origin=segiterator->positionAtDepthincm(currentDepth_);
    }
  else
    {
      if(EMSHOWER)
    origin=segiterator->positionAtDepthinX0(currentDepth_);
      if(HADSHOWER)
    origin=segiterator->positionAtDepthinL0(currentDepth_);
    }
  XYZVector zaxis(0,0,1);
  XYZVector planeVec1=(zaxis.Cross(particleDirection)).Unit();
  locToGlobal_=Transform3D(Point(0,0,0),
               Point(0,0,1),
               Point(1,0,0),
               (Point)origin,
               (Point)(origin+particleDirection),
               (Point)(origin+planeVec1));
  return true;
}