SimG4HGCalValidation.cc

Go to the documentation of this file.

// File: SimG4HGCalValidation.cc
// Description: Main analysis class for HGCal Validation of G4 Hits

#include "SimG4Core/Notification/interface/BeginOfJob.h"
#include "SimG4Core/Notification/interface/BeginOfRun.h"
#include "SimG4Core/Notification/interface/BeginOfEvent.h"
#include "SimG4Core/Notification/interface/EndOfEvent.h"
#include "SimG4Core/Watcher/interface/SimProducer.h"
#include "SimG4Core/Notification/interface/Observer.h"

// to retreive hits
#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
#include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
#include "DataFormats/Math/interface/Point3D.h"
#include "SimDataFormats/ValidationFormats/interface/PHGCalValidInfo.h"
#include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h"
#include "DataFormats/HcalDetId/interface/HcalTestNumbering.h"
#include "SimG4CMS/Calo/interface/HGCNumberingScheme.h"
#include "SimG4CMS/Calo/interface/HGCalNumberingScheme.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "Geometry/Records/interface/IdealGeometryRecord.h"
#include "Geometry/Records/interface/HcalSimNumberingRecord.h"
#include "Geometry/HcalCommonData/interface/HcalNumberingFromDDD.h"
#include "Geometry/HcalCommonData/interface/HcalDDDSimConstants.h"
#include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h"
#include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"

#include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
#include "FWCore/PluginManager/interface/ModuleDef.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "G4Step.hh"
#include "G4Track.hh"
#include "G4NavigationHistory.hh"
#include "G4TouchableHistory.hh"

#include "CLHEP/Units/GlobalSystemOfUnits.h"
#include "CLHEP/Units/GlobalPhysicalConstants.h"

#include <cmath>
#include <iostream>
#include <iomanip>
#include <memory>
#include <vector>
#include <string>

class SimG4HGCalValidation : public SimProducer,
                 public Observer<const BeginOfJob *>, 
                 public Observer<const BeginOfEvent *>, 
                 public Observer<const G4Step *> {

public:
  SimG4HGCalValidation(const edm::ParameterSet &p);
  ~SimG4HGCalValidation() override;

  void produce(edm::Event&, const edm::EventSetup&) override;

private:
  SimG4HGCalValidation(const SimG4HGCalValidation&); // stop default
  const SimG4HGCalValidation& operator=(const SimG4HGCalValidation&);

  void  init();

  // observer classes
  void update(const BeginOfJob * job) override;
  void update(const BeginOfEvent * evt) override;
  void update(const G4Step * step) override;

  // analysis related class
  void layerAnalysis(PHGCalValidInfo&);
  void clear();

private:
  //Keep reference to instantiate HcalNumberingFromDDD later
  HcalNumberingFromDDD *             numberingFromDDD_;

  //HGCal numbering scheme
  std::vector<HGCNumberingScheme*>   hgcNumbering_;
  std::vector<HGCalNumberingScheme*> hgcalNumbering_;
  
  // to read from parameter set
  std::vector<std::string>           names_;
  std::vector<int>                   types_, detTypes_, subdet_;
  std::string                        labelLayer_;

  // parameters from geometry
  int                                levelT1_, levelT2_;

  // some private members for ananlysis 
  unsigned int                       count_;
  int                                verbosity_;
  double                             edepEE_, edepHEF_, edepHEB_;
  std::vector<double>                hgcEEedep_, hgcHEFedep_, hgcHEBedep_;
  std::vector<unsigned int>          dets_, hgchitDets_, hgchitIndex_;
  std::vector<double>                hgchitX_, hgchitY_, hgchitZ_; 
};

SimG4HGCalValidation::SimG4HGCalValidation(const edm::ParameterSet &p): 
  numberingFromDDD_(nullptr), levelT1_(999), levelT2_(999), count_(0) {

  edm::ParameterSet m_Anal = p.getParameter<edm::ParameterSet>("SimG4HGCalValidation");
  names_        = m_Anal.getParameter<std::vector<std::string> >("Names");
  types_        = m_Anal.getParameter<std::vector<int> >("Types");
  detTypes_     = m_Anal.getParameter<std::vector<int> >("DetTypes");
  labelLayer_   = m_Anal.getParameter<std::string>("LabelLayerInfo");
  verbosity_    = m_Anal.getUntrackedParameter<int>("Verbosity",0);

  produces<PHGCalValidInfo>(labelLayer_);

  if (verbosity_ > 0) {
    edm::LogVerbatim("ValidHGCal") << "HGCalTestAnalysis:: Initialised as "
                   << "observer of begin events and of G4step "
                   << "with Parameter values: \n\tLabel : " 
                   << labelLayer_ << " and with " 
                   << names_.size() << " detectors";
    for (unsigned int k=0; k<names_.size(); ++k)
      edm::LogVerbatim("ValidHGCal") << " [" << k << "] " << names_[k] 
                     << " Type " << types_[k] << " DetType "
                     << detTypes_[k];
  }
} 
   
SimG4HGCalValidation::~SimG4HGCalValidation() {
  delete numberingFromDDD_;
  for (auto number : hgcNumbering_)   delete number;
  for (auto number : hgcalNumbering_) delete number;
}

void SimG4HGCalValidation::produce(edm::Event& e, const edm::EventSetup&) {

  std::unique_ptr<PHGCalValidInfo> productLayer(new PHGCalValidInfo);
  layerAnalysis(*productLayer);
  e.put(std::move(productLayer),labelLayer_);
}

void SimG4HGCalValidation::update(const BeginOfJob * job) {

  const edm::EventSetup* es = (*job)();
  for (unsigned int type=0; type<types_.size(); ++type) {
    int      layers(0);
    int      detType = detTypes_[type];
    G4String nameX = "HGCal";
    if (types_[type] <= 1) {
      if (types_[type] == 0) {
    subdet_.emplace_back(ForwardSubdetector::ForwardEmpty);
    if      (detType == 0) dets_.emplace_back((int)(DetId::HGCalEE));
    else if (detType == 1) dets_.emplace_back((int)(DetId::HGCalHSi));
    else                   dets_.emplace_back((int)(DetId::HGCalHSc));
      } else {
    dets_.push_back((unsigned int)(DetId::Forward));
    if      (detType == 0) subdet_.emplace_back((int)(ForwardSubdetector::HGCEE));
    else if (detType == 1) subdet_.emplace_back((int)(ForwardSubdetector::HGCHEF));
    else                   subdet_.emplace_back((int)(ForwardSubdetector::HGCHEB));
      }
      if        (detType == 0) nameX = "HGCalEESensitive";
      else if   (detType == 1) nameX = "HGCalHESiliconSensitive";
      else                     nameX = "HGCalHEScintillatorSensitive";
      edm::ESHandle<HGCalDDDConstants>    hdc;
      es->get<IdealGeometryRecord>().get(nameX,hdc);
      if (hdc.isValid()) {
    HGCalGeometryMode::GeometryMode m_mode = hdc->geomMode();
    levelT1_ = hdc->levelTop(0); 
    levelT2_ = hdc->levelTop(1); 
    if (m_mode == HGCalGeometryMode::Trapezoid) {
      types_[type] =-1;
      hgcalNumbering_.emplace_back(new HGCalNumberingScheme(*hdc,(DetId::Detector)(dets_[type]),nameX));
    } else if ((m_mode == HGCalGeometryMode::Hexagon) ||
           (m_mode == HGCalGeometryMode::HexagonFull)) {
      types_[type] = 1;
      hgcNumbering_.push_back(new HGCNumberingScheme(*hdc,nameX));
    } else {
      types_[type] = 0;
      hgcalNumbering_.emplace_back(new HGCalNumberingScheme(*hdc,(DetId::Detector)(dets_[type]),nameX));
    }
    layers = hdc->layers(false);
      } else {
    edm::LogError("ValidHGCal") << "Cannot find HGCalDDDConstants for "
                    << nameX;
    throw cms::Exception("Unknown", "ValidHGCal") 
      << "Cannot find HGCalDDDConstants for " << nameX << "\n";
      }
    } else {
      nameX = "HcalEndcap";
      dets_.push_back((unsigned int)(DetId::Hcal));
      subdet_.push_back((int)(HcalSubdetector::HcalEndcap));
      edm::ESHandle<HcalDDDSimConstants>    hdc;
      es->get<HcalSimNumberingRecord>().get(hdc);
      if (hdc.isValid()) {
    numberingFromDDD_ = new HcalNumberingFromDDD(hdc.product());
    layers = 18;
      } else {
    edm::LogError("ValidHGCal") << "Cannot find HcalDDDSimConstant";
    throw cms::Exception("Unknown", "ValidHGCal") 
      << "Cannot find HcalDDDSimConstant\n";
      }
    }
    if        (detType == 0) {
      for (int i=0; i<layers; ++i) hgcEEedep_.push_back(0);
    } else if (detType == 1) {
      for (int i=0; i<layers; ++i) hgcHEFedep_.push_back(0);
    } else                   {
      for (int i=0; i<layers; ++i) hgcHEBedep_.push_back(0);
    }
    if (verbosity_ > 0)
      edm::LogVerbatim("ValidHGCal") << "[" << type << "]: " << nameX 
                     << " det " << dets_[type] << " subdet " 
                     << subdet_[type] << " with " << layers 
                     << " layers";
  }
}

//=================================================================== per EVENT
void SimG4HGCalValidation::update(const BeginOfEvent * evt) {
 
  int iev = (*evt)()->GetEventID();
  if (verbosity_ > 0) 
    edm::LogVerbatim("ValidHGCal") << "SimG4HGCalValidation: =====> Begin "
                   << "event = " << iev;
  
  ++count_;
  edepEE_ = edepHEF_ = edepHEB_ = 0.;

  //HGCal variables 
  for (unsigned int i = 0; i<hgcEEedep_.size();  i++) hgcEEedep_[i]  = 0.;
  for (unsigned int i = 0; i<hgcHEFedep_.size(); i++) hgcHEFedep_[i] = 0.; 
  for (unsigned int i = 0; i<hgcHEBedep_.size(); i++) hgcHEBedep_[i] = 0.; 

  //Cache reset  
  clear();
}

//=================================================================== each STEP
void SimG4HGCalValidation::update(const G4Step * aStep) {

  if (aStep != nullptr) {
    G4VPhysicalVolume* curPV  = aStep->GetPreStepPoint()->GetPhysicalVolume();
    const G4String& name = curPV->GetName();
    G4VSensitiveDetector* curSD = aStep->GetPreStepPoint()->GetSensitiveDetector();
    if (verbosity_ > 1)
      edm::LogVerbatim("ValidHGCal") << "ValidHGCal::Step in " << name <<" at "
                     << aStep->GetPreStepPoint()->GetPosition();

    // Only for Sensitive detector
    if (curSD != nullptr) {
      int type(-1);
      for (unsigned int k=0; k<names_.size(); ++k) {
    if (name.find(names_[k].c_str()) != std::string::npos) {
      type = k; break;
    }
      }
      int detType = (type >= 0) ? detTypes_[type] : -1;
      if (verbosity_ > 0)
    edm::LogVerbatim("ValidHGCal") << "ValidHGCal::In SD " << name 
                       << " type " << type << ":" << detType;

      // Right type of SD
      if (type >= 0) {
    //Get the 32-bit index of the hit cell
    const G4TouchableHistory* touchable = static_cast<const G4TouchableHistory*>(aStep->GetPreStepPoint()->GetTouchable());
    unsigned int index(0);
    int          layer(0);
    G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
    if (types_[type] <= 1) {
      // HGCal
      G4ThreeVector localpos = touchable->GetHistory()->GetTopTransform().TransformPoint(hitPoint);   
      float globalZ = touchable->GetTranslation(0).z();
      int iz(globalZ>0 ? 1 : -1);
      int module(-1), cell(-1);
      if (types_[type] == 1) {
        if (touchable->GetHistoryDepth() == levelT1_) {
          layer  = touchable->GetReplicaNumber(0);
        } else {
          layer  = touchable->GetReplicaNumber(2);
          module = touchable->GetReplicaNumber(1);
          cell   = touchable->GetReplicaNumber(0);
        }
        index = hgcNumbering_[type]->getUnitID((ForwardSubdetector)(subdet_[type]), layer, module, cell, iz, localpos);
      } else {
        if ((touchable->GetHistoryDepth() == levelT1_) || 
        (touchable->GetHistoryDepth() == levelT2_)) {
          layer  = touchable->GetReplicaNumber(0);
        } else {
          layer  = touchable->GetReplicaNumber(3);
          module = touchable->GetReplicaNumber(2);
          cell   = touchable->GetReplicaNumber(1);
        }
        double weight(0);
        index = hgcalNumbering_[type]->getUnitID(layer, module, cell, iz,
                             hitPoint, weight);
      }
      if (verbosity_ > 1)
        edm::LogVerbatim("ValidHGCal") << "HGCal: " << name << " Layer "
                       << layer << " Module " << module
                       << " Cell " << cell;
    } else {
      // Hcal
      int depth = (touchable->GetReplicaNumber(0))%10 + 1;
      int lay   = (touchable->GetReplicaNumber(0)/10)%100 + 1;
      int det   = (touchable->GetReplicaNumber(1))/1000;
      HcalNumberingFromDDD::HcalID tmp = 
        numberingFromDDD_->unitID(det, 
                      math::XYZVectorD(hitPoint.x(),hitPoint.y(),
                               hitPoint.z()), depth, lay);
      index = HcalTestNumbering::packHcalIndex(tmp.subdet,tmp.zside,tmp.depth,tmp.etaR,tmp.phis,tmp.lay);
      layer = tmp.lay;
      if (verbosity_ > 1)
        edm::LogVerbatim("ValidHGCal") << "HCAL: " << det << ":" << depth 
                       << ":" << lay << " o/p " 
                       << tmp.subdet << ":" << tmp.zside 
                       << ":" << tmp.depth << ":" 
                       << tmp.etaR << ":" << tmp.phis 
                       << ":" << tmp.lay << " point "
                       << hitPoint << " " << hitPoint.rho()
                       << ":"  << hitPoint.eta() << ":" 
                       << hitPoint.phi();
    }

    double edeposit = aStep->GetTotalEnergyDeposit();
    if (verbosity_ > 0)
      edm::LogVerbatim("ValidHGCal") << "Layer " << layer << " Index " 
                     << std::hex << index << std::dec
                     << " Edep " << edeposit << " hit at "
                     << hitPoint;
    if        (detType == 0) {
      edepEE_  += edeposit;
      if (layer < (int)(hgcEEedep_.size()))
        hgcEEedep_[layer]  += edeposit;
    } else if (detType == 1) {
      edepHEF_ += edeposit;
      if (layer < (int)(hgcHEFedep_.size())) 
        hgcHEFedep_[layer] += edeposit;
    } else                   {
      edepHEB_ += edeposit;
      if (layer < (int)(hgcHEBedep_.size())) 
        hgcHEBedep_[layer] += edeposit;
    }
    G4String nextVolume("XXX");
    if (aStep->GetTrack()->GetNextVolume()!=nullptr)
      nextVolume = aStep->GetTrack()->GetNextVolume()->GetName();

    if (nextVolume.c_str()!=name.c_str()) { //save hit when it exits cell
      if (std::find(hgchitIndex_.begin(),hgchitIndex_.end(),index) == hgchitIndex_.end()) {
        hgchitDets_.push_back(dets_[type]);
        hgchitIndex_.push_back(index);
        hgchitX_.push_back(hitPoint.x());
        hgchitY_.push_back(hitPoint.y()); 
        hgchitZ_.push_back(hitPoint.z());
      }
    }
      } // it is right type of SD
    } // it is in a SD
  }//if aStep!=NULL
}//end update aStep

//================================================================ End of EVENT

void SimG4HGCalValidation::layerAnalysis(PHGCalValidInfo& product) {

  if (verbosity_ > 0)
    edm::LogVerbatim("ValidHGCal") << "\n ===>>> SimG4HGCalValidation: Energy "
                   << "deposit\n at EE : " << std::setw(6) 
                   << edepEE_/CLHEP::MeV << "\n at HEF: " 
                   << std::setw(6) << edepHEF_/CLHEP::MeV 
                   << "\n at HEB: " << std::setw(6)
                   << edepHEB_/CLHEP::MeV;
  
  //Fill HGC Variables
  product.fillhgcHits(hgchitDets_,hgchitIndex_, hgchitX_, hgchitY_, hgchitZ_);
  product.fillhgcLayers(edepEE_,edepHEF_,edepHEB_,hgcEEedep_, hgcHEFedep_,hgcHEBedep_);
}

//---------------------------------------------------
void SimG4HGCalValidation::clear() {

  hgchitDets_.erase(hgchitDets_.begin(),hgchitDets_.end()); 
  hgchitIndex_.erase(hgchitIndex_.begin(),hgchitIndex_.end()); 
  hgchitX_.erase(hgchitX_.begin(),hgchitX_.end()); 
  hgchitY_.erase(hgchitY_.begin(),hgchitY_.end()); 
  hgchitZ_.erase(hgchitZ_.begin(),hgchitZ_.end());      
}

DEFINE_SIMWATCHER (SimG4HGCalValidation);