Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:30:05

0001 // -*- C++ -*-
0002 //
0003 // Package:     HcalTestBeam
0004 // Class  :     HcalTB04Analysis
0005 //
0006 // Implementation:
0007 //     Main analysis class for Hcal Test Beam 2004 Analysis
0008 //
0009 //  Usage: A Simwatcher class and can be activated from Oscarproducer module
0010 //
0011 // Original Author:
0012 //         Created:  Tue May 16 10:14:34 CEST 2006
0013 //
0014 
0015 // system include files
0016 #include <cmath>
0017 #include <iomanip>
0018 #include <iostream>
0019 #include <memory>
0020 #include <vector>
0021 #include <string>
0022 
0023 // user include files
0024 #include "SimG4Core/Notification/interface/Observer.h"
0025 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0026 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0027 #include "SimG4Core/Notification/interface/EndOfEvent.h"
0028 #include "SimG4Core/Watcher/interface/SimProducer.h"
0029 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0030 
0031 // to retreive hits
0032 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0033 #include "SimDataFormats/CaloHit/interface/CaloHit.h"
0034 #include "SimDataFormats/HcalTestBeam/interface/PHcalTB04Info.h"
0035 #include "SimG4CMS/Calo/interface/ECalSD.h"
0036 #include "SimG4CMS/Calo/interface/HCalSD.h"
0037 #include "SimG4CMS/Calo/interface/HcalQie.h"
0038 #include "SimG4CMS/Calo/interface/CaloG4Hit.h"
0039 #include "SimG4CMS/Calo/interface/CaloG4HitCollection.h"
0040 #include "SimG4CMS/Calo/interface/HcalTestNumberingScheme.h"
0041 #include "DataFormats/Math/interface/Point3D.h"
0042 
0043 #include "FWCore/Framework/interface/Event.h"
0044 #include "FWCore/Framework/interface/EventSetup.h"
0045 #include "FWCore/Framework/interface/MakerMacros.h"
0046 #include "FWCore/PluginManager/interface/ModuleDef.h"
0047 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0048 
0049 #include "SimG4CMS/HcalTestBeam/interface/HcalTBNumberingScheme.h"
0050 #include "SimG4CMS/HcalTestBeam/interface/HcalTB04Histo.h"
0051 #include "SimG4CMS/HcalTestBeam/interface/HcalTB04XtalNumberingScheme.h"
0052 
0053 #include "G4SDManager.hh"
0054 #include "G4Step.hh"
0055 #include "G4Track.hh"
0056 #include "G4ThreeVector.hh"
0057 #include "G4VProcess.hh"
0058 #include "G4HCofThisEvent.hh"
0059 
0060 #include <CLHEP/Random/RandGaussQ.h>
0061 #include <CLHEP/Random/Randomize.h>
0062 #include <CLHEP/Units/GlobalSystemOfUnits.h>
0063 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0064 
0065 #include <cstdint>
0066 
0067 //#define EDM_ML_DEBUG
0068 
0069 namespace CLHEP {
0070   class HepRandomEngine;
0071 }
0072 
0073 class HcalTB04Analysis : public SimProducer,
0074                          public Observer<const BeginOfRun*>,
0075                          public Observer<const BeginOfEvent*>,
0076                          public Observer<const EndOfEvent*>,
0077                          public Observer<const G4Step*> {
0078 public:
0079   HcalTB04Analysis(const edm::ParameterSet& p);
0080   HcalTB04Analysis(const HcalTB04Analysis&) = delete;  // stop default
0081   const HcalTB04Analysis& operator=(const HcalTB04Analysis&) = delete;
0082   ~HcalTB04Analysis() override;
0083 
0084   void produce(edm::Event&, const edm::EventSetup&) override;
0085 
0086 private:
0087   void init();
0088 
0089   // observer methods
0090   void update(const BeginOfRun* run) override;
0091   void update(const BeginOfEvent* evt) override;
0092   void update(const G4Step* step) override;
0093   void update(const EndOfEvent* evt) override;
0094 
0095   //User methods
0096   void fillBuffer(const EndOfEvent* evt);
0097   void qieAnalysis(CLHEP::HepRandomEngine*);
0098   void xtalAnalysis(CLHEP::HepRandomEngine*);
0099   void finalAnalysis();
0100   void fillEvent(PHcalTB04Info&);
0101 
0102   void clear();
0103   int unitID(uint32_t id);
0104   double scale(int det, int layer);
0105   double timeOfFlight(int det, int layer, double eta);
0106 
0107 private:
0108   // to read from parameter set
0109   const edm::ParameterSet m_Anal;
0110   const bool hcalOnly;
0111   const int mode, type;
0112   const double ecalNoise, beamOffset;
0113   const double scaleHB0, scaleHB16, scaleHO, scaleHE0;
0114   const std::vector<std::string> names;
0115 
0116   HcalQie* myQie;
0117   HcalTB04Histo* histo;
0118 
0119   int iceta, icphi;
0120   G4RotationMatrix* beamline_RM;
0121 
0122   // Constants for the run
0123   int count;
0124   int nTower, nCrystal;
0125   std::vector<int> idHcal, idXtal;
0126   std::vector<uint32_t> idTower, idEcal;
0127 
0128   // Constants for the event
0129   int nPrimary, particleType;
0130   double pInit, etaInit, phiInit;
0131   std::vector<CaloHit> ecalHitCache;
0132   std::vector<CaloHit> hcalHitCache, hcalHitLayer;
0133   std::vector<double> esimh, eqie, esime, enois;
0134   std::vector<double> eseta, eqeta, esphi, eqphi, eslay, eqlay;
0135   double etots, eecals, ehcals, etotq, eecalq, ehcalq;
0136 
0137   bool pvFound;
0138   int evNum, pvType;
0139   G4ThreeVector pvPosition, pvMomentum, pvUVW;
0140   std::vector<int> secTrackID, secPartID;
0141   std::vector<G4ThreeVector> secMomentum;
0142   std::vector<double> secEkin;
0143   std::vector<int> shortLivedSecondaries;
0144 };
0145 
0146 //
0147 // constructors and destructor
0148 //
0149 
0150 HcalTB04Analysis::HcalTB04Analysis(const edm::ParameterSet& p)
0151     : m_Anal(p.getParameter<edm::ParameterSet>("HcalTB04Analysis")),
0152       hcalOnly(m_Anal.getParameter<bool>("HcalOnly")),
0153       mode(m_Anal.getParameter<int>("Mode")),
0154       type(m_Anal.getParameter<int>("Type")),
0155       ecalNoise(m_Anal.getParameter<double>("EcalNoise")),
0156       beamOffset(-m_Anal.getParameter<double>("BeamPosition") * CLHEP::cm),
0157       scaleHB0(m_Anal.getParameter<double>("ScaleHB0")),
0158       scaleHB16(m_Anal.getParameter<double>("ScaleHB16")),
0159       scaleHO(m_Anal.getParameter<double>("ScaleHO")),
0160       scaleHE0(m_Anal.getParameter<double>("ScaleHE0")),
0161       names(m_Anal.getParameter<std::vector<std::string> >("Names")),
0162       myQie(nullptr),
0163       histo(nullptr) {
0164   double fMinEta = m_Anal.getParameter<double>("MinEta");
0165   double fMaxEta = m_Anal.getParameter<double>("MaxEta");
0166   double fMinPhi = m_Anal.getParameter<double>("MinPhi");
0167   double fMaxPhi = m_Anal.getParameter<double>("MaxPhi");
0168   double beamEta = (fMaxEta + fMinEta) / 2.;
0169   double beamPhi = (fMaxPhi + fMinPhi) / 2.;
0170   double beamThet = 2 * atan(exp(-beamEta));
0171   if (beamPhi < 0)
0172     beamPhi += twopi;
0173   iceta = static_cast<int>(beamEta / 0.087) + 1;
0174   icphi = static_cast<int>(std::fabs(beamPhi) / 0.087) + 5;
0175   if (icphi > 72)
0176     icphi -= 73;
0177 
0178   produces<PHcalTB04Info>();
0179 
0180   beamline_RM = new G4RotationMatrix;
0181   beamline_RM->rotateZ(-beamPhi);
0182   beamline_RM->rotateY(-beamThet);
0183 
0184   edm::LogVerbatim("HcalTBSim")
0185       << "HcalTB04:: Initialised as observer of BeginOf Job/BeginOfRun/BeginOfEvent/G4Step/EndOfEvent with Parameter "
0186          "values:\n \thcalOnly = "
0187       << hcalOnly << "\tecalNoise = " << ecalNoise << "\n\tMode = " << mode << " (0: HB2 Standard; 1:HB2 Segmented)"
0188       << "\tType = " << type << " (0: HB; 1 HE; 2 HB+HE)\n\tbeamOffset = " << beamOffset << "\ticeta = " << iceta
0189       << "\ticphi = " << icphi << "\n\tbeamline_RM = " << *beamline_RM;
0190 
0191   init();
0192 
0193   myQie = new HcalQie(p);
0194   histo = new HcalTB04Histo(m_Anal);
0195 }
0196 
0197 HcalTB04Analysis::~HcalTB04Analysis() {
0198 #ifdef EDM_ML_DEBUG
0199   edm::LogVerbatim("HcalTBSim") << "\n -------->  Total number of selected entries : " << count << "\nPointers:: QIE "
0200                                 << myQie << " Histo " << histo;
0201 #endif
0202   if (myQie) {
0203     delete myQie;
0204     myQie = nullptr;
0205   }
0206   if (histo) {
0207     delete histo;
0208     histo = nullptr;
0209   }
0210 }
0211 
0212 //
0213 // member functions
0214 //
0215 
0216 void HcalTB04Analysis::produce(edm::Event& e, const edm::EventSetup&) {
0217   std::unique_ptr<PHcalTB04Info> product(new PHcalTB04Info);
0218   fillEvent(*product);
0219   e.put(std::move(product));
0220 }
0221 
0222 void HcalTB04Analysis::init() {
0223   idTower = HcalTBNumberingScheme::getUnitIDs(type, mode);
0224   nTower = idTower.size();
0225 #ifdef EDM_ML_DEBUG
0226   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Save information from " << nTower << " HCal towers";
0227 #endif
0228   idHcal.reserve(nTower);
0229   for (int i = 0; i < nTower; i++) {
0230     int id = unitID(idTower[i]);
0231     idHcal.push_back(id);
0232 #ifdef EDM_ML_DEBUG
0233     edm::LogVerbatim("HcalTBSim") << "\tTower[" << i << "] Original " << std::hex << idTower[i] << " Stored "
0234                                   << idHcal[i] << std::dec;
0235 #endif
0236   }
0237 
0238   if (!hcalOnly) {
0239     int det = 10;
0240     uint32_t id1;
0241     nCrystal = 0;
0242     for (int lay = 1; lay < 8; lay++) {
0243       for (int icr = 1; icr < 8; icr++) {
0244         id1 = HcalTestNumbering::packHcalIndex(det, 0, 1, icr, lay, 1);
0245         int id = unitID(id1);
0246         idEcal.push_back(id1);
0247         idXtal.push_back(id);
0248         nCrystal++;
0249       }
0250     }
0251     edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Save information from " << nCrystal << " ECal Crystals";
0252 #ifdef EDM_ML_DEBUG
0253     for (int i = 0; i < nCrystal; i++) {
0254       edm::LogVerbatim("HcalTBSim") << "\tCrystal[" << i << "] Original " << std::hex << idEcal[i] << " Stored "
0255                                     << idXtal[i] << std::dec;
0256     }
0257 #endif
0258   }
0259   // Profile vectors
0260   eseta.reserve(5);
0261   eqeta.reserve(5);
0262   esphi.reserve(3);
0263   eqphi.reserve(3);
0264   eslay.reserve(20);
0265   eqlay.reserve(20);
0266   for (int i = 0; i < 5; i++) {
0267     eseta.push_back(0.);
0268     eqeta.push_back(0.);
0269   }
0270   for (int i = 0; i < 3; i++) {
0271     esphi.push_back(0.);
0272     eqphi.push_back(0.);
0273   }
0274   for (int i = 0; i < 20; i++) {
0275     eslay.push_back(0.);
0276     eqlay.push_back(0.);
0277   }
0278 
0279   // counter
0280   count = 0;
0281   evNum = 0;
0282   clear();
0283 }
0284 
0285 void HcalTB04Analysis::update(const BeginOfRun* run) {
0286   int irun = (*run)()->GetRunID();
0287   edm::LogVerbatim("HcalTBSim") << " =====> Begin of Run = " << irun;
0288 
0289   G4SDManager* sd = G4SDManager::GetSDMpointerIfExist();
0290   if (sd != nullptr) {
0291     std::string sdname = names[0];
0292     G4VSensitiveDetector* aSD = sd->FindSensitiveDetector(sdname);
0293     if (aSD == nullptr) {
0294       edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::beginOfRun: No SD"
0295                                    << " with name " << sdname << " in this "
0296                                    << "Setup";
0297     } else {
0298       HCalSD* theCaloSD = dynamic_cast<HCalSD*>(aSD);
0299       edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::beginOfRun: Finds SD with name " << theCaloSD->GetName()
0300                                     << " in this Setup";
0301       HcalNumberingScheme* org = new HcalTestNumberingScheme(false);
0302       theCaloSD->setNumberingScheme(org);
0303       edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::beginOfRun: set a new numbering scheme";
0304     }
0305     if (!hcalOnly) {
0306       sdname = names[1];
0307       aSD = sd->FindSensitiveDetector(sdname);
0308       if (aSD == nullptr) {
0309         edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::beginOfRun: No SD"
0310                                      << " with name " << sdname << " in this "
0311                                      << "Setup";
0312       } else {
0313         ECalSD* theCaloSD = dynamic_cast<ECalSD*>(aSD);
0314         edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::beginOfRun: Finds SD with name " << theCaloSD->GetName()
0315                                       << " in this Setup";
0316         EcalNumberingScheme* org = new HcalTB04XtalNumberingScheme();
0317         theCaloSD->setNumberingScheme(org);
0318         edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::beginOfRun: set a new numbering scheme";
0319       }
0320     }
0321   } else {
0322     edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::beginOfRun: Could "
0323                                  << "not get SD Manager!";
0324   }
0325 }
0326 
0327 void HcalTB04Analysis::update(const BeginOfEvent* evt) {
0328   clear();
0329 #ifdef EDM_ML_DEBUG
0330   evNum = (*evt)()->GetEventID();
0331   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis: =====> Begin of event = " << evNum;
0332 #endif
0333 }
0334 
0335 void HcalTB04Analysis::update(const G4Step* aStep) {
0336   if (aStep != nullptr) {
0337     //Get Step properties
0338     G4ThreeVector thePreStepPoint = aStep->GetPreStepPoint()->GetPosition();
0339     G4ThreeVector thePostStepPoint;
0340 
0341     // Get Tracks properties
0342     G4Track* aTrack = aStep->GetTrack();
0343     int trackID = aTrack->GetTrackID();
0344     int parentID = aTrack->GetParentID();
0345     const G4ThreeVector& position = aTrack->GetPosition();
0346     G4ThreeVector momentum = aTrack->GetMomentum();
0347     G4String partType = aTrack->GetDefinition()->GetParticleType();
0348     G4String partSubType = aTrack->GetDefinition()->GetParticleSubType();
0349     int partPDGEncoding = aTrack->GetDefinition()->GetPDGEncoding();
0350 #ifdef EDM_ML_DEBUG
0351     bool isPDGStable = aTrack->GetDefinition()->GetPDGStable();
0352 #endif
0353     double pDGlifetime = aTrack->GetDefinition()->GetPDGLifeTime();
0354     double gammaFactor = aStep->GetPreStepPoint()->GetGamma();
0355 
0356     if (!pvFound) {  //search for v1
0357       double stepDeltaEnergy = aStep->GetDeltaEnergy();
0358       double kinEnergy = aTrack->GetKineticEnergy();
0359 
0360       // look for DeltaE > 10% kinEnergy of particle, or particle death - Ek=0
0361       if (trackID == 1 && parentID == 0 && ((kinEnergy == 0.) || (std::fabs(stepDeltaEnergy / kinEnergy) > 0.1))) {
0362         pvType = -1;
0363         if (kinEnergy == 0.) {
0364           pvType = 0;
0365         } else {
0366           if (std::fabs(stepDeltaEnergy / kinEnergy) > 0.1)
0367             pvType = 1;
0368         }
0369         pvFound = true;
0370         pvPosition = position;
0371         pvMomentum = momentum;
0372         // Rotated coord.system:
0373         pvUVW = (*beamline_RM) * (pvPosition);
0374 
0375         //Volume name requires some checks:
0376         G4String thePostPVname = "NoName";
0377         G4StepPoint* thePostPoint = aStep->GetPostStepPoint();
0378         if (thePostPoint) {
0379           thePostStepPoint = thePostPoint->GetPosition();
0380           G4VPhysicalVolume* thePostPV = thePostPoint->GetPhysicalVolume();
0381           if (thePostPV)
0382             thePostPVname = thePostPV->GetName();
0383         }
0384 #ifdef EDM_ML_DEBUG
0385         edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: V1 found at: " << thePostStepPoint
0386                                       << " G4VPhysicalVolume: " << thePostPVname;
0387         edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::fill_v1Pos: Primary Track momentum: " << pvMomentum
0388                                       << " psoition " << pvPosition << " u/v/w " << pvUVW;
0389 #endif
0390       }
0391     } else {
0392       // watch for secondaries originating @v1, including the surviving primary
0393       if ((trackID != 1 && parentID == 1 && (aTrack->GetCurrentStepNumber() == 1) && (thePreStepPoint == pvPosition)) ||
0394           (trackID == 1 && thePreStepPoint == pvPosition)) {
0395 #ifdef EDM_ML_DEBUG
0396         edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::A secondary...  PDG:" << partPDGEncoding
0397                                       << " TrackID:" << trackID << " ParentID:" << parentID
0398                                       << " stable: " << isPDGStable << " Tau: " << pDGlifetime
0399                                       << " cTauGamma=" << c_light * pDGlifetime * gammaFactor * 1000.
0400                                       << "um GammaFactor: " << gammaFactor;
0401 #endif
0402         secTrackID.push_back(trackID);
0403         secPartID.push_back(partPDGEncoding);
0404         secMomentum.push_back(momentum);
0405         secEkin.push_back(aTrack->GetKineticEnergy());
0406 
0407         // Check for short-lived secondaries: cTauGamma<100um
0408         double ctaugamma_um = CLHEP::c_light * pDGlifetime * gammaFactor * 1000.;
0409         if ((ctaugamma_um > 0.) && (ctaugamma_um < 100.)) {  //short-lived secondary
0410           shortLivedSecondaries.push_back(trackID);
0411         } else {  //normal secondary - enter into the V1-calorimetric tree
0412                   //          histos->fill_v1cSec (aTrack);
0413         }
0414       }
0415       // Also watch for tertiary particles coming from
0416       // short-lived secondaries from V1
0417       if (aTrack->GetCurrentStepNumber() == 1) {
0418         if (!shortLivedSecondaries.empty()) {
0419           int pid = parentID;
0420           std::vector<int>::iterator pos1 = shortLivedSecondaries.begin();
0421           std::vector<int>::iterator pos2 = shortLivedSecondaries.end();
0422           std::vector<int>::iterator pos;
0423           for (pos = pos1; pos != pos2; pos++) {
0424             if (*pos == pid) {  //ParentID is on the list of short-lived
0425                                 // secondary
0426 #ifdef EDM_ML_DEBUG
0427               edm::LogVerbatim("HcalTBSim")
0428                   << "HcalTB04Analysis:: A tertiary...  PDG:" << partPDGEncoding << " TrackID:" << trackID
0429                   << " ParentID:" << parentID << " stable: " << isPDGStable << " Tau: " << pDGlifetime
0430                   << " cTauGamma=" << c_light * pDGlifetime * gammaFactor * 1000. << "um GammaFactor: " << gammaFactor;
0431 #endif
0432             }
0433           }
0434         }
0435       }
0436     }
0437   }
0438 }
0439 
0440 void HcalTB04Analysis::update(const EndOfEvent* evt) {
0441   count++;
0442 
0443   //fill the buffer
0444 #ifdef EDM_ML_DEBUG
0445   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::Fill event " << (*evt)()->GetEventID();
0446 #endif
0447   fillBuffer(evt);
0448 
0449   //QIE analysis
0450 #ifdef EDM_ML_DEBUG
0451   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::Do QIE analysis with " << hcalHitCache.size() << " hits";
0452 #endif
0453   CLHEP::HepRandomEngine* engine = G4Random::getTheEngine();
0454   qieAnalysis(engine);
0455 
0456   //Energy in Crystal Matrix
0457   if (!hcalOnly) {
0458 #ifdef EDM_ML_DEBUG
0459     edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::Do Xtal analysis with " << ecalHitCache.size() << " hits";
0460 #endif
0461     xtalAnalysis(engine);
0462   }
0463 
0464   //Final Analysis
0465 #ifdef EDM_ML_DEBUG
0466   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::Final analysis";
0467 #endif
0468   finalAnalysis();
0469 
0470   int iEvt = (*evt)()->GetEventID();
0471   if (iEvt < 10)
0472     edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Event " << iEvt;
0473   else if ((iEvt < 100) && (iEvt % 10 == 0))
0474     edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Event " << iEvt;
0475   else if ((iEvt < 1000) && (iEvt % 100 == 0))
0476     edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Event " << iEvt;
0477   else if ((iEvt < 10000) && (iEvt % 1000 == 0))
0478     edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Event " << iEvt;
0479 }
0480 
0481 void HcalTB04Analysis::fillBuffer(const EndOfEvent* evt) {
0482   std::vector<CaloHit> hhits, hhitl;
0483   int idHC, j;
0484   CaloG4HitCollection* theHC;
0485   std::map<int, float, std::less<int> > primaries;
0486 #ifdef EDM_ML_DEBUG
0487   double etot1 = 0, etot2 = 0;
0488 #endif
0489 
0490   // Look for the Hit Collection of HCal
0491   G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
0492   std::string sdName = names[0];
0493   idHC = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
0494   theHC = (CaloG4HitCollection*)allHC->GetHC(idHC);
0495 #ifdef EDM_ML_DEBUG
0496   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Hit Collection for " << sdName << " of ID " << idHC
0497                                 << " is obtained at " << theHC << " with " << theHC->entries() << " entries";
0498 #endif
0499   int thehc_entries = theHC->entries();
0500   if (idHC >= 0 && theHC != nullptr) {
0501     hhits.reserve(theHC->entries());
0502     hhitl.reserve(theHC->entries());
0503     for (j = 0; j < thehc_entries; j++) {
0504       CaloG4Hit* aHit = (*theHC)[j];
0505       double e = aHit->getEnergyDeposit() / CLHEP::GeV;
0506       double time = aHit->getTimeSlice();
0507       math::XYZPoint pos = aHit->getEntry();
0508       unsigned int id = aHit->getUnitID();
0509       double theta = pos.theta();
0510       double eta = -std::log(std::tan(theta * 0.5));
0511       double phi = pos.phi();
0512       int det, z, group, ieta, iphi, layer;
0513       HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, layer);
0514       double jitter = time - timeOfFlight(det, layer, eta);
0515       if (jitter < 0)
0516         jitter = 0;
0517       if (e < 0 || e > 1.)
0518         e = 0;
0519       double escl = e * scale(det, layer);
0520       unsigned int idx = HcalTBNumberingScheme::getUnitID(id, mode);
0521       CaloHit hit(det, layer, escl, eta, phi, jitter, idx);
0522       hhits.push_back(hit);
0523       CaloHit hitl(det, layer, escl, eta, phi, jitter, id);
0524       hhitl.push_back(hitl);
0525       primaries[aHit->getTrackID()] += e;
0526 #ifdef EDM_ML_DEBUG
0527       etot1 += escl;
0528       edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Hcal Hit i/p " << j << "  ID 0x" << std::hex << id << " 0x"
0529                                     << idx << std::dec << " time " << std::setw(6) << time << " " << std::setw(6)
0530                                     << jitter << " theta " << std::setw(8) << theta << " eta " << std::setw(8) << eta
0531                                     << " phi " << std::setw(8) << phi << " e " << std::setw(8) << e << " "
0532                                     << std::setw(8) << escl;
0533 #endif
0534     }
0535   }
0536 
0537   // Add hits in the same channel within same time slice
0538   std::vector<CaloHit>::iterator itr;
0539   int nHit = hhits.size();
0540   std::vector<CaloHit*> hits(nHit);
0541   for (j = 0, itr = hhits.begin(); itr != hhits.end(); j++, itr++) {
0542     hits[j] = &hhits[j];
0543   }
0544   sort(hits.begin(), hits.end(), CaloHitIdMore());
0545   std::vector<CaloHit*>::iterator k1, k2;
0546 #ifdef EDM_ML_DEBUG
0547   int nhit = 0;
0548 #endif
0549   for (k1 = hits.begin(); k1 != hits.end(); k1++) {
0550     int det = (**k1).det();
0551     int layer = (**k1).layer();
0552     double ehit = (**k1).e();
0553     double eta = (**k1).eta();
0554     double phi = (**k1).phi();
0555     double jitter = (**k1).t();
0556     uint32_t unitID = (**k1).id();
0557     int jump = 0;
0558     for (k2 = k1 + 1; k2 != hits.end() && std::fabs(jitter - (**k2).t()) < 1 && unitID == (**k2).id(); k2++) {
0559       ehit += (**k2).e();
0560       jump++;
0561     }
0562     CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
0563     hcalHitCache.push_back(hit);
0564     k1 += jump;
0565 #ifdef EDM_ML_DEBUG
0566     nhit++;
0567     etot2 += ehit;
0568     edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Hcal Hit store " << nhit << "  ID 0x" << std::hex << unitID
0569                                   << std::dec << " time " << std::setw(6) << jitter << " eta " << std::setw(8) << eta
0570                                   << " phi " << std::setw(8) << phi << " e " << std::setw(8) << ehit;
0571 #endif
0572   }
0573 #ifdef EDM_ML_DEBUG
0574   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Stores " << nhit << " HCal hits from " << nHit
0575                                 << " input hits E(Hcal) " << etot1 << " " << etot2;
0576 #endif
0577   //Repeat for Hit in each layer (hhits and hhitl sizes are the same)
0578   for (j = 0, itr = hhitl.begin(); itr != hhitl.end(); j++, itr++) {
0579     hits[j] = &hhitl[j];
0580   }
0581   sort(hits.begin(), hits.end(), CaloHitIdMore());
0582 #ifdef EDM_ML_DEBUG
0583   int nhitl = 0;
0584   double etotl = 0;
0585 #endif
0586   for (k1 = hits.begin(); k1 != hits.end(); k1++) {
0587     int det = (**k1).det();
0588     int layer = (**k1).layer();
0589     double ehit = (**k1).e();
0590     double eta = (**k1).eta();
0591     double phi = (**k1).phi();
0592     double jitter = (**k1).t();
0593     uint32_t unitID = (**k1).id();
0594     int jump = 0;
0595     for (k2 = k1 + 1; k2 != hits.end() && std::fabs(jitter - (**k2).t()) < 1 && unitID == (**k2).id(); k2++) {
0596       ehit += (**k2).e();
0597       jump++;
0598     }
0599     CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
0600     hcalHitLayer.push_back(hit);
0601     k1 += jump;
0602 #ifdef EDM_ML_DEBUG
0603     nhitl++;
0604     etotl += ehit;
0605     edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Hcal Hit store " << nhitl << "  ID 0x" << std::hex << unitID
0606                                   << std::dec << " time " << std::setw(6) << jitter << " eta " << std::setw(8) << eta
0607                                   << " phi " << std::setw(8) << phi << " e " << std::setw(8) << ehit;
0608 #endif
0609   }
0610 #ifdef EDM_ML_DEBUG
0611   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Stores " << nhitl << " HCal hits from " << nHit
0612                                 << " input hits E(Hcal) " << etot1 << " " << etotl;
0613 #endif
0614   // Look for the Hit Collection of ECal
0615   std::vector<CaloHit> ehits;
0616   sdName = names[1];
0617   idHC = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
0618   theHC = (CaloG4HitCollection*)allHC->GetHC(idHC);
0619 #ifdef EDM_ML_DEBUG
0620   etot1 = etot2 = 0;
0621   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Hit Collection for " << sdName << " of ID " << idHC
0622                                 << " is obtained at " << theHC << " with " << theHC->entries() << " entries";
0623 #endif
0624   if (idHC >= 0 && theHC != nullptr) {
0625     thehc_entries = theHC->entries();
0626     ehits.reserve(theHC->entries());
0627     for (j = 0; j < thehc_entries; j++) {
0628       CaloG4Hit* aHit = (*theHC)[j];
0629       double e = aHit->getEnergyDeposit() / CLHEP::GeV;
0630       double time = aHit->getTimeSlice();
0631       if (e < 0 || e > 100000.)
0632         e = 0;
0633       if (e > 0) {
0634         math::XYZPoint pos = aHit->getEntry();
0635         unsigned int id = aHit->getUnitID();
0636         double theta = pos.theta();
0637         double eta = -std::log(std::tan(theta * 0.5));
0638         double phi = pos.phi();
0639         int det, z, group, ieta, iphi, layer;
0640         HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, layer);
0641         CaloHit hit(det, 0, e, eta, phi, time, id);
0642         ehits.push_back(hit);
0643         primaries[aHit->getTrackID()] += e;
0644 #ifdef EDM_ML_DEBUG
0645         etot1 += e;
0646         edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Ecal Hit i/p " << j << "  ID 0x" << std::hex << id
0647                                       << std::dec << " time " << std::setw(6) << time << " theta " << std::setw(8)
0648                                       << theta << " eta " << std::setw(8) << eta << " phi " << std::setw(8) << phi
0649                                       << " e " << std::setw(8) << e;
0650 #endif
0651       }
0652     }
0653   }
0654 
0655   // Add hits in the same channel within same time slice
0656   nHit = ehits.size();
0657   std::vector<CaloHit*> hite(nHit);
0658   for (j = 0, itr = ehits.begin(); itr != ehits.end(); j++, itr++) {
0659     hite[j] = &ehits[j];
0660   }
0661   sort(hite.begin(), hite.end(), CaloHitIdMore());
0662 #ifdef EDM_ML_DEBUG
0663   nhit = 0;
0664 #endif
0665   for (k1 = hite.begin(); k1 != hite.end(); k1++) {
0666     int det = (**k1).det();
0667     int layer = (**k1).layer();
0668     double ehit = (**k1).e();
0669     double eta = (**k1).eta();
0670     double phi = (**k1).phi();
0671     double jitter = (**k1).t();
0672     uint32_t unitID = (**k1).id();
0673     int jump = 0;
0674     for (k2 = k1 + 1; k2 != hite.end() && std::fabs(jitter - (**k2).t()) < 1 && unitID == (**k2).id(); k2++) {
0675       ehit += (**k2).e();
0676       jump++;
0677     }
0678     CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
0679     ecalHitCache.push_back(hit);
0680     k1 += jump;
0681 #ifdef EDM_ML_DEBUG
0682     nhit++;
0683     etot2 += ehit;
0684     edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Ecal Hit store " << nhit << "  ID 0x" << std::hex << unitID
0685                                   << std::dec << " time " << std::setw(6) << jitter << " eta " << std::setw(8) << eta
0686                                   << " phi " << std::setw(8) << phi << " e " << std::setw(8) << ehit;
0687 #endif
0688   }
0689 #ifdef EDM_ML_DEBUG
0690   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Stores " << nhit << " ECal hits from " << nHit
0691                                 << " input hits E(Ecal) " << etot1 << " " << etot2;
0692 #endif
0693   // Find Primary info:
0694   nPrimary = static_cast<int>(primaries.size());
0695   int trackID = 0;
0696   G4PrimaryParticle* thePrim = nullptr;
0697   int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
0698 #ifdef EDM_ML_DEBUG
0699   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Event has " << nvertex << " verteices";
0700 #endif
0701   if (nvertex <= 0)
0702     edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::EndOfEvent ERROR: no vertex found for event " << evNum;
0703   for (int i = 0; i < nvertex; i++) {
0704     G4PrimaryVertex* avertex = (*evt)()->GetPrimaryVertex(i);
0705     if (avertex == nullptr) {
0706       edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::EndOfEvent ERR: pointer to vertex = 0 for event " << evNum;
0707     } else {
0708       edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::Vertex number :" << i << " " << avertex->GetPosition();
0709       int npart = avertex->GetNumberOfParticle();
0710       if (npart == 0)
0711         edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::End Of Event ERR: no primary!";
0712       if (thePrim == nullptr)
0713         thePrim = avertex->GetPrimary(trackID);
0714     }
0715   }
0716 
0717   if (thePrim != nullptr) {
0718     double px = thePrim->GetPx();
0719     double py = thePrim->GetPy();
0720     double pz = thePrim->GetPz();
0721     double p = std::sqrt(pow(px, 2.) + pow(py, 2.) + pow(pz, 2.));
0722     pInit = p / CLHEP::GeV;
0723     if (p == 0)
0724       edm::LogWarning("HcalTBSim") << "HcalTB04Analysis:: EndOfEvent ERR: primary has p=0 ";
0725     else {
0726       double costheta = pz / p;
0727       double theta = acos(std::min(std::max(costheta, -1.), 1.));
0728       etaInit = -std::log(std::tan(theta / 2));
0729       if (px != 0 || py != 0)
0730         phiInit = std::atan2(py, px);
0731     }
0732     particleType = thePrim->GetPDGcode();
0733   } else
0734     edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::EndOfEvent ERR: could not find primary";
0735 }
0736 
0737 void HcalTB04Analysis::qieAnalysis(CLHEP::HepRandomEngine* engine) {
0738   int hittot = hcalHitCache.size();
0739   if (hittot <= 0)
0740     hittot = 1;
0741   std::vector<CaloHit> hits(hittot);
0742   std::vector<int> todo(nTower, 0);
0743 
0744 #ifdef EDM_ML_DEBUG
0745   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::qieAnalysis: Size " << hits.size() << " " << todo.size() << " "
0746                                 << idTower.size() << " " << esimh.size() << " " << eqie.size();
0747 #endif
0748   // Loop over all HCal hits
0749   for (unsigned int k1 = 0; k1 < hcalHitCache.size(); k1++) {
0750     CaloHit hit = hcalHitCache[k1];
0751     uint32_t id = hit.id();
0752     int nhit = 0;
0753     double esim = hit.e();
0754     hits[nhit] = hit;
0755     for (unsigned int k2 = k1 + 1; k2 < hcalHitCache.size(); k2++) {
0756       hit = hcalHitCache[k2];
0757       if (hit.id() == id) {
0758         nhit++;
0759         hits[nhit] = hit;
0760         esim += hit.e();
0761       }
0762     }
0763     k1 += nhit;
0764     nhit++;
0765     std::vector<int> cd = myQie->getCode(nhit, hits, engine);
0766     double eq = myQie->getEnergy(cd);
0767 #ifdef EDM_ML_DEBUG
0768     edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::  ID 0x" << std::hex << id << std::dec << " registers " << esim
0769                                   << " energy from " << nhit << " hits starting with hit # " << k1
0770                                   << " energy with noise " << eq;
0771 #endif
0772     for (int k2 = 0; k2 < nTower; k2++) {
0773       if (id == idTower[k2]) {
0774         todo[k2] = 1;
0775         esimh[k2] = esim;
0776         eqie[k2] = eq;
0777       }
0778     }
0779   }
0780 
0781   // Towers with no hit
0782   for (int k2 = 0; k2 < nTower; k2++) {
0783     if (todo[k2] == 0) {
0784       std::vector<int> cd = myQie->getCode(0, hits, engine);
0785       double eq = myQie->getEnergy(cd);
0786       esimh[k2] = 0;
0787       eqie[k2] = eq;
0788 #ifdef EDM_ML_DEBUG
0789       edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::  ID 0x" << std::hex << idTower[k2] << std::dec
0790                                     << " registers " << esimh[k2] << " energy from hits and energy after QIE analysis "
0791                                     << eqie[k2];
0792 #endif
0793     }
0794   }
0795 }
0796 
0797 void HcalTB04Analysis::xtalAnalysis(CLHEP::HepRandomEngine* engine) {
0798   CLHEP::RandGaussQ randGauss(*engine);
0799 
0800   // Crystal Data
0801   std::vector<int> iok(nCrystal, 0);
0802 #ifdef EDM_ML_DEBUG
0803   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::xtalAnalysis: Size " << iok.size() << " " << idEcal.size() << " "
0804                                 << esime.size() << " " << enois.size();
0805 #endif
0806   for (unsigned int k1 = 0; k1 < ecalHitCache.size(); k1++) {
0807     uint32_t id = ecalHitCache[k1].id();
0808     int nhit = 0;
0809     double esim = ecalHitCache[k1].e();
0810     for (unsigned int k2 = k1 + 1; k2 < ecalHitCache.size(); k2++) {
0811       if (ecalHitCache[k2].id() == id) {
0812         nhit++;
0813         esim += ecalHitCache[k2].e();
0814       }
0815     }
0816     k1 += nhit;
0817     nhit++;
0818     double eq = esim + randGauss.fire(0., ecalNoise);
0819 #ifdef EDM_ML_DEBUG
0820     edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::  ID 0x" << std::hex << id << std::dec << " registers " << esim
0821                                   << " energy from " << nhit << " hits starting with hit # " << k1
0822                                   << " energy with noise " << eq;
0823 #endif
0824     for (int k2 = 0; k2 < nCrystal; k2++) {
0825       if (id == idEcal[k2]) {
0826         iok[k2] = 1;
0827         esime[k2] = esim;
0828         enois[k2] = eq;
0829       }
0830     }
0831   }
0832 
0833   // Crystals with no hit
0834   for (int k2 = 0; k2 < nCrystal; k2++) {
0835     if (iok[k2] == 0) {
0836       esime[k2] = 0;
0837       enois[k2] = randGauss.fire(0., ecalNoise);
0838 #ifdef EDM_ML_DEBUG
0839       edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::  ID 0x" << std::hex << idEcal[k2] << std::dec
0840                                     << " registers " << esime[k2] << " energy from hits and energy from noise "
0841                                     << enois[k2];
0842 #endif
0843     }
0844   }
0845 }
0846 
0847 void HcalTB04Analysis::finalAnalysis() {
0848   //Beam Information
0849   histo->fillPrimary(pInit, etaInit, phiInit);
0850 
0851   // Total Energy
0852   eecals = ehcals = eecalq = ehcalq = 0.;
0853   for (int i = 0; i < nTower; i++) {
0854     ehcals += esimh[i];
0855     ehcalq += eqie[i];
0856   }
0857   for (int i = 0; i < nCrystal; i++) {
0858     eecals += esime[i];
0859     eecalq += enois[i];
0860   }
0861   etots = eecals + ehcals;
0862   etotq = eecalq + ehcalq;
0863 #ifdef EDM_ML_DEBUG
0864   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Energy deposit at Sim Level (Total) " << etots << " (ECal) "
0865                                 << eecals << " (HCal) " << ehcals
0866                                 << "\nHcalTB04Analysis:: Energy deposit at Qie Level (Total) " << etotq << " (ECal) "
0867                                 << eecalq << " (HCal) " << ehcalq;
0868 #endif
0869   histo->fillEdep(etots, eecals, ehcals, etotq, eecalq, ehcalq);
0870 
0871   // Lateral Profile
0872   for (int i = 0; i < 5; i++) {
0873     eseta[i] = 0.;
0874     eqeta[i] = 0.;
0875   }
0876   for (int i = 0; i < 3; i++) {
0877     esphi[i] = 0.;
0878     eqphi[i] = 0.;
0879   }
0880   double e1 = 0, e2 = 0;
0881   unsigned int id;
0882   for (int i = 0; i < nTower; i++) {
0883     int det, z, group, ieta, iphi, layer;
0884     id = idTower[i];
0885     HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, layer);
0886     iphi -= (icphi - 1);
0887     if (icphi > 4) {
0888       if (ieta == 0)
0889         ieta = 2;
0890       else
0891         ieta = -1;
0892     } else {
0893       ieta = ieta - iceta + 2;
0894     }
0895     if (iphi >= 0 && iphi < 3 && ieta >= 0 && ieta < 5) {
0896       eseta[ieta] += esimh[i];
0897       esphi[iphi] += esimh[i];
0898       e1 += esimh[i];
0899       eqeta[ieta] += eqie[i];
0900       eqphi[iphi] += eqie[i];
0901       e2 += eqie[i];
0902     }
0903   }
0904   for (int i = 0; i < 3; i++) {
0905     if (e1 > 0)
0906       esphi[i] /= e1;
0907     if (e2 > 0)
0908       eqphi[i] /= e2;
0909   }
0910   for (int i = 0; i < 5; i++) {
0911     if (e1 > 0)
0912       eseta[i] /= e1;
0913     if (e2 > 0)
0914       eqeta[i] /= e2;
0915   }
0916 #ifdef EDM_ML_DEBUG
0917   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Energy fraction along Eta and Phi (Sim/Qie)";
0918   for (int i = 0; i < 5; i++)
0919     edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: [" << i << "] Eta Sim = " << eseta[i] << " Qie = " << eqeta[i]
0920                                   << " Phi Sim = " << esphi[i] << " Qie = " << eqphi[i];
0921 #endif
0922   histo->fillTrnsProf(eseta, eqeta, esphi, eqphi);
0923 
0924   // Longitudianl profile
0925   for (int i = 0; i < 20; i++) {
0926     eslay[i] = 0.;
0927     eqlay[i] = 0.;
0928   }
0929   e1 = 0;
0930   e2 = 0;
0931   for (int i = 0; i < nTower; i++) {
0932     int det, z, group, ieta, iphi, layer;
0933     id = idTower[i];
0934     HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, layer);
0935     iphi -= (icphi - 1);
0936     layer -= 1;
0937     if (iphi >= 0 && iphi < 3 && layer >= 0 && layer < 20) {
0938       eslay[layer] += esimh[i];
0939       e1 += esimh[i];
0940       eqlay[layer] += eqie[i];
0941       e2 += eqie[i];
0942     }
0943   }
0944   for (int i = 0; i < 20; i++) {
0945     if (e1 > 0)
0946       eslay[i] /= e1;
0947     if (e2 > 0)
0948       eqlay[i] /= e2;
0949   }
0950 #ifdef EDM_ML_DEBUG
0951   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Energy fraction along Layer";
0952   for (int i = 0; i < 20; i++)
0953     edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: [" << i << "] Sim = " << eslay[i] << " Qie = " << eqlay[i];
0954 #endif
0955   histo->fillLongProf(eslay, eqlay);
0956 }
0957 
0958 void HcalTB04Analysis::fillEvent(PHcalTB04Info& product) {
0959   //Setup the ID's
0960   product.setIDs(idHcal, idXtal);
0961 
0962   //Beam Information
0963   product.setPrimary(nPrimary, particleType, pInit, etaInit, phiInit);
0964 
0965   //Energy deposits in the crystals and towers
0966   product.setEdepHcal(esimh, eqie);
0967   product.setEdepHcal(esime, enois);
0968 
0969   // Total Energy
0970   product.setEdep(etots, eecals, ehcals, etotq, eecalq, ehcalq);
0971 
0972   // Lateral Profile
0973   product.setTrnsProf(eseta, eqeta, esphi, eqphi);
0974 
0975   // Longitudianl profile
0976   product.setLongProf(eslay, eqlay);
0977 
0978   //Save Hits
0979   int nhit = 0;
0980   std::vector<CaloHit>::iterator itr;
0981   for (itr = ecalHitCache.begin(); itr != ecalHitCache.end(); itr++) {
0982     uint32_t id = itr->id();
0983     int det, z, group, ieta, iphi, lay;
0984     HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, lay);
0985     product.saveHit(det, lay, ieta, iphi, itr->e(), itr->t());
0986     nhit++;
0987 #ifdef EDM_ML_DEBUG
0988     edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Save Hit " << std::setw(3) << nhit << " ID 0x" << std::hex
0989                                   << group << std::dec << " " << std::setw(2) << det << " " << std::setw(2) << lay
0990                                   << " " << std::setw(1) << z << " " << std::setw(3) << ieta << " " << std::setw(3)
0991                                   << iphi << " T " << std::setw(6) << itr->t() << " E " << std::setw(6) << itr->e();
0992 #endif
0993   }
0994   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Saves " << nhit << " hits from Crystals";
0995 #ifdef EDM_ML_DEBUG
0996   int nhit0 = nhit;
0997 #endif
0998   nhit = 0;
0999   for (itr = hcalHitCache.begin(); itr != hcalHitCache.end(); itr++) {
1000     uint32_t id = itr->id();
1001     int det, z, group, ieta, iphi, lay;
1002     HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, lay);
1003     product.saveHit(det, lay, ieta, iphi, itr->e(), itr->t());
1004     nhit++;
1005 #ifdef EDM_ML_DEBUG
1006     edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Save Hit " << std::setw(3) << nhit + nhit0 << " ID 0x"
1007                                   << std::hex << group << std::dec << " " << std::setw(2) << det << " " << std::setw(2)
1008                                   << lay << " " << std::setw(1) << z << " " << std::setw(3) << ieta << " "
1009                                   << std::setw(3) << iphi << " T " << std::setw(6) << itr->t() << " E " << std::setw(6)
1010                                   << itr->e();
1011 #endif
1012   }
1013   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Saves " << nhit << " hits from HCal";
1014 
1015   //Vertex associated quantities
1016   product.setVtxPrim(evNum,
1017                      pvType,
1018                      pvPosition.x(),
1019                      pvPosition.y(),
1020                      pvPosition.z(),
1021                      pvUVW.x(),
1022                      pvUVW.y(),
1023                      pvUVW.z(),
1024                      pvMomentum.x(),
1025                      pvMomentum.y(),
1026                      pvMomentum.z());
1027   for (unsigned int i = 0; i < secTrackID.size(); i++) {
1028     product.setVtxSec(
1029         secTrackID[i], secPartID[i], secMomentum[i].x(), secMomentum[i].y(), secMomentum[i].z(), secEkin[i]);
1030   }
1031 }
1032 
1033 void HcalTB04Analysis::clear() {
1034   pvFound = false;
1035   pvType = -2;
1036   pvPosition = G4ThreeVector();
1037   pvMomentum = G4ThreeVector();
1038   pvUVW = G4ThreeVector();
1039   secTrackID.clear();
1040   secPartID.clear();
1041   secMomentum.clear();
1042   secEkin.clear();
1043   shortLivedSecondaries.clear();
1044 
1045   ecalHitCache.erase(ecalHitCache.begin(), ecalHitCache.end());
1046   hcalHitCache.erase(hcalHitCache.begin(), hcalHitCache.end());
1047   hcalHitLayer.erase(hcalHitLayer.begin(), hcalHitLayer.end());
1048   nPrimary = particleType = 0;
1049   pInit = etaInit = phiInit = 0;
1050 
1051   esimh.clear();
1052   eqie.clear();
1053   esimh.reserve(nTower);
1054   eqie.reserve(nTower);
1055   for (int i = 0; i < nTower; i++) {
1056     esimh.push_back(0.);
1057     eqie.push_back(0.);
1058   }
1059   esime.clear();
1060   enois.clear();
1061   esime.reserve(nCrystal);
1062   enois.reserve(nCrystal);
1063   for (int i = 0; i < nCrystal; i++) {
1064     esime.push_back(0.);
1065     enois.push_back(0.);
1066   }
1067 }
1068 
1069 int HcalTB04Analysis::unitID(uint32_t id) {
1070   int det, z, group, ieta, iphi, lay;
1071   HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, lay);
1072   group = (det & 15) << 20;
1073   group += ((lay - 1) & 31) << 15;
1074   group += (z & 1) << 14;
1075   group += (ieta & 127) << 7;
1076   group += (iphi & 127);
1077   return group;
1078 }
1079 
1080 double HcalTB04Analysis::scale(int det, int layer) {
1081   double tmp = 1.;
1082   if (det == static_cast<int>(HcalBarrel)) {
1083     if (layer == 1)
1084       tmp = scaleHB0;
1085     else if (layer == 17)
1086       tmp = scaleHB16;
1087     else if (layer > 17)
1088       tmp = scaleHO;
1089   } else {
1090     if (layer <= 2)
1091       tmp = scaleHE0;
1092   }
1093   return tmp;
1094 }
1095 
1096 double HcalTB04Analysis::timeOfFlight(int det, int layer, double eta) {
1097   double theta = 2.0 * std::atan(std::exp(-eta));
1098   double dist = beamOffset;
1099   if (det == static_cast<int>(HcalBarrel)) {
1100     const double rLay[19] = {1836.0,
1101                              1902.0,
1102                              1962.0,
1103                              2022.0,
1104                              2082.0,
1105                              2142.0,
1106                              2202.0,
1107                              2262.0,
1108                              2322.0,
1109                              2382.0,
1110                              2448.0,
1111                              2514.0,
1112                              2580.0,
1113                              2646.0,
1114                              2712.0,
1115                              2776.0,
1116                              2862.5,
1117                              3847.0,
1118                              4052.0};
1119     if (layer > 0 && layer <= 19)
1120       dist += rLay[layer - 1] * mm / sin(theta);
1121   } else {
1122     const double zLay[19] = {4034.0,
1123                              4032.0,
1124                              4123.0,
1125                              4210.0,
1126                              4297.0,
1127                              4384.0,
1128                              4471.0,
1129                              4558.0,
1130                              4645.0,
1131                              4732.0,
1132                              4819.0,
1133                              4906.0,
1134                              4993.0,
1135                              5080.0,
1136                              5167.0,
1137                              5254.0,
1138                              5341.0,
1139                              5428.0,
1140                              5515.0};
1141     if (layer > 0 && layer <= 19)
1142       dist += zLay[layer - 1] * mm / cos(theta);
1143   }
1144 
1145   double tmp = dist / c_light / ns;
1146 #ifdef EDM_ML_DEBUG
1147   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::timeOfFlight " << tmp << " for det/lay " << det << " " << layer
1148                                 << " eta/theta " << eta << " " << theta / deg << " dist " << dist;
1149 #endif
1150   return tmp;
1151 }
1152 
1153 DEFINE_SIMWATCHER(HcalTB04Analysis);