Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-02-23 01:39:02 

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   double etot1 = 0, etot2 = 0;
0487 
0488   // Look for the Hit Collection of HCal
0489   G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
0490   std::string sdName = names[0];
0491   idHC = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
0492   theHC = (CaloG4HitCollection*)allHC->GetHC(idHC);
0493 #ifdef EDM_ML_DEBUG
0494   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Hit Collection for " << sdName << " of ID " << idHC
0495                                 << " is obtained at " << theHC << " with " << theHC->entries() << " entries";
0496 #endif
0497   int thehc_entries = theHC->entries();
0498   if (idHC >= 0 && theHC != nullptr) {
0499     hhits.reserve(theHC->entries());
0500     hhitl.reserve(theHC->entries());
0501     for (j = 0; j < thehc_entries; j++) {
0502       CaloG4Hit* aHit = (*theHC)[j];
0503       double e = aHit->getEnergyDeposit() / CLHEP::GeV;
0504       double time = aHit->getTimeSlice();
0505       math::XYZPoint pos = aHit->getEntry();
0506       unsigned int id = aHit->getUnitID();
0507       double theta = pos.theta();
0508       double eta = -std::log(std::tan(theta * 0.5));
0509       double phi = pos.phi();
0510       int det, z, group, ieta, iphi, layer;
0511       HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, layer);
0512       double jitter = time - timeOfFlight(det, layer, eta);
0513       if (jitter < 0)
0514         jitter = 0;
0515       if (e < 0 || e > 1.)
0516         e = 0;
0517       double escl = e * scale(det, layer);
0518       unsigned int idx = HcalTBNumberingScheme::getUnitID(id, mode);
0519       CaloHit hit(det, layer, escl, eta, phi, jitter, idx);
0520       hhits.push_back(hit);
0521       CaloHit hitl(det, layer, escl, eta, phi, jitter, id);
0522       hhitl.push_back(hitl);
0523       primaries[aHit->getTrackID()] += e;
0524       etot1 += escl;
0525 #ifdef EDM_ML_DEBUG
0526       edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Hcal Hit i/p " << j << "  ID 0x" << std::hex << id << " 0x"
0527                                     << idx << std::dec << " time " << std::setw(6) << time << " " << std::setw(6)
0528                                     << jitter << " theta " << std::setw(8) << theta << " eta " << std::setw(8) << eta
0529                                     << " phi " << std::setw(8) << phi << " e " << std::setw(8) << e << " "
0530                                     << std::setw(8) << escl;
0531 #endif
0532     }
0533   }
0534 
0535   // Add hits in the same channel within same time slice
0536   std::vector<CaloHit>::iterator itr;
0537   int nHit = hhits.size();
0538   std::vector<CaloHit*> hits(nHit);
0539   for (j = 0, itr = hhits.begin(); itr != hhits.end(); j++, itr++) {
0540     hits[j] = &hhits[j];
0541   }
0542   sort(hits.begin(), hits.end(), CaloHitIdMore());
0543   std::vector<CaloHit*>::iterator k1, k2;
0544   int nhit = 0;
0545   for (k1 = hits.begin(); k1 != hits.end(); k1++) {
0546     int det = (**k1).det();
0547     int layer = (**k1).layer();
0548     double ehit = (**k1).e();
0549     double eta = (**k1).eta();
0550     double phi = (**k1).phi();
0551     double jitter = (**k1).t();
0552     uint32_t unitID = (**k1).id();
0553     int jump = 0;
0554     for (k2 = k1 + 1; k2 != hits.end() && std::fabs(jitter - (**k2).t()) < 1 && unitID == (**k2).id(); k2++) {
0555       ehit += (**k2).e();
0556       jump++;
0557     }
0558     nhit++;
0559     CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
0560     hcalHitCache.push_back(hit);
0561     etot2 += ehit;
0562     k1 += jump;
0563 #ifdef EDM_ML_DEBUG
0564     edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Hcal Hit store " << nhit << "  ID 0x" << std::hex << unitID
0565                                   << std::dec << " time " << std::setw(6) << jitter << " eta " << std::setw(8) << eta
0566                                   << " phi " << std::setw(8) << phi << " e " << std::setw(8) << ehit;
0567 #endif
0568   }
0569 #ifdef EDM_ML_DEBUG
0570   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Stores " << nhit << " HCal hits from " << nHit
0571                                 << " input hits E(Hcal) " << etot1 << " " << etot2;
0572 #endif
0573   //Repeat for Hit in each layer (hhits and hhitl sizes are the same)
0574   for (j = 0, itr = hhitl.begin(); itr != hhitl.end(); j++, itr++) {
0575     hits[j] = &hhitl[j];
0576   }
0577   sort(hits.begin(), hits.end(), CaloHitIdMore());
0578   int nhitl = 0;
0579   double etotl = 0;
0580   for (k1 = hits.begin(); k1 != hits.end(); k1++) {
0581     int det = (**k1).det();
0582     int layer = (**k1).layer();
0583     double ehit = (**k1).e();
0584     double eta = (**k1).eta();
0585     double phi = (**k1).phi();
0586     double jitter = (**k1).t();
0587     uint32_t unitID = (**k1).id();
0588     int jump = 0;
0589     for (k2 = k1 + 1; k2 != hits.end() && std::fabs(jitter - (**k2).t()) < 1 && unitID == (**k2).id(); k2++) {
0590       ehit += (**k2).e();
0591       jump++;
0592     }
0593     nhitl++;
0594     CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
0595     hcalHitLayer.push_back(hit);
0596     etotl += ehit;
0597     k1 += jump;
0598 #ifdef EDM_ML_DEBUG
0599     edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Hcal Hit store " << nhitl << "  ID 0x" << std::hex << unitID
0600                                   << std::dec << " time " << std::setw(6) << jitter << " eta " << std::setw(8) << eta
0601                                   << " phi " << std::setw(8) << phi << " e " << std::setw(8) << ehit;
0602 #endif
0603   }
0604 #ifdef EDM_ML_DEBUG
0605   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Stores " << nhitl << " HCal hits from " << nHit
0606                                 << " input hits E(Hcal) " << etot1 << " " << etotl;
0607 #endif
0608   // Look for the Hit Collection of ECal
0609   std::vector<CaloHit> ehits;
0610   sdName = names[1];
0611   idHC = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
0612   theHC = (CaloG4HitCollection*)allHC->GetHC(idHC);
0613   etot1 = etot2 = 0;
0614 #ifdef EDM_ML_DEBUG
0615   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Hit Collection for " << sdName << " of ID " << idHC
0616                                 << " is obtained at " << theHC << " with " << theHC->entries() << " entries";
0617 #endif
0618   if (idHC >= 0 && theHC != nullptr) {
0619     thehc_entries = theHC->entries();
0620     ehits.reserve(theHC->entries());
0621     for (j = 0; j < thehc_entries; j++) {
0622       CaloG4Hit* aHit = (*theHC)[j];
0623       double e = aHit->getEnergyDeposit() / CLHEP::GeV;
0624       double time = aHit->getTimeSlice();
0625       if (e < 0 || e > 100000.)
0626         e = 0;
0627       if (e > 0) {
0628         math::XYZPoint pos = aHit->getEntry();
0629         unsigned int id = aHit->getUnitID();
0630         double theta = pos.theta();
0631         double eta = -std::log(std::tan(theta * 0.5));
0632         double phi = pos.phi();
0633         int det, z, group, ieta, iphi, layer;
0634         HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, layer);
0635         CaloHit hit(det, 0, e, eta, phi, time, id);
0636         ehits.push_back(hit);
0637         primaries[aHit->getTrackID()] += e;
0638         etot1 += e;
0639 #ifdef EDM_ML_DEBUG
0640         edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Ecal Hit i/p " << j << "  ID 0x" << std::hex << id
0641                                       << std::dec << " time " << std::setw(6) << time << " theta " << std::setw(8)
0642                                       << theta << " eta " << std::setw(8) << eta << " phi " << std::setw(8) << phi
0643                                       << " e " << std::setw(8) << e;
0644 #endif
0645       }
0646     }
0647   }
0648 
0649   // Add hits in the same channel within same time slice
0650   nHit = ehits.size();
0651   std::vector<CaloHit*> hite(nHit);
0652   for (j = 0, itr = ehits.begin(); itr != ehits.end(); j++, itr++) {
0653     hite[j] = &ehits[j];
0654   }
0655   sort(hite.begin(), hite.end(), CaloHitIdMore());
0656   nhit = 0;
0657   for (k1 = hite.begin(); k1 != hite.end(); k1++) {
0658     int det = (**k1).det();
0659     int layer = (**k1).layer();
0660     double ehit = (**k1).e();
0661     double eta = (**k1).eta();
0662     double phi = (**k1).phi();
0663     double jitter = (**k1).t();
0664     uint32_t unitID = (**k1).id();
0665     int jump = 0;
0666     for (k2 = k1 + 1; k2 != hite.end() && std::fabs(jitter - (**k2).t()) < 1 && unitID == (**k2).id(); k2++) {
0667       ehit += (**k2).e();
0668       jump++;
0669     }
0670     nhit++;
0671     CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
0672     ecalHitCache.push_back(hit);
0673     etot2 += ehit;
0674     k1 += jump;
0675 #ifdef EDM_ML_DEBUG
0676     edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Ecal Hit store " << nhit << "  ID 0x" << std::hex << unitID
0677                                   << std::dec << " time " << std::setw(6) << jitter << " eta " << std::setw(8) << eta
0678                                   << " phi " << std::setw(8) << phi << " e " << std::setw(8) << ehit;
0679 #endif
0680   }
0681 #ifdef EDM_ML_DEBUG
0682   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Stores " << nhit << " ECal hits from " << nHit
0683                                 << " input hits E(Ecal) " << etot1 << " " << etot2;
0684 #endif
0685   // Find Primary info:
0686   nPrimary = static_cast<int>(primaries.size());
0687   int trackID = 0;
0688   G4PrimaryParticle* thePrim = nullptr;
0689   int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
0690 #ifdef EDM_ML_DEBUG
0691   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Event has " << nvertex << " verteices";
0692 #endif
0693   if (nvertex <= 0)
0694     edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::EndOfEvent ERROR: no vertex found for event " << evNum;
0695   for (int i = 0; i < nvertex; i++) {
0696     G4PrimaryVertex* avertex = (*evt)()->GetPrimaryVertex(i);
0697     if (avertex == nullptr) {
0698       edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::EndOfEvent ERR: pointer to vertex = 0 for event " << evNum;
0699     } else {
0700       edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::Vertex number :" << i << " " << avertex->GetPosition();
0701       int npart = avertex->GetNumberOfParticle();
0702       if (npart == 0)
0703         edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::End Of Event ERR: no primary!";
0704       if (thePrim == nullptr)
0705         thePrim = avertex->GetPrimary(trackID);
0706     }
0707   }
0708 
0709   if (thePrim != nullptr) {
0710     double px = thePrim->GetPx();
0711     double py = thePrim->GetPy();
0712     double pz = thePrim->GetPz();
0713     double p = std::sqrt(pow(px, 2.) + pow(py, 2.) + pow(pz, 2.));
0714     pInit = p / CLHEP::GeV;
0715     if (p == 0)
0716       edm::LogWarning("HcalTBSim") << "HcalTB04Analysis:: EndOfEvent ERR: primary has p=0 ";
0717     else {
0718       double costheta = pz / p;
0719       double theta = acos(std::min(std::max(costheta, -1.), 1.));
0720       etaInit = -std::log(std::tan(theta / 2));
0721       if (px != 0 || py != 0)
0722         phiInit = std::atan2(py, px);
0723     }
0724     particleType = thePrim->GetPDGcode();
0725   } else
0726     edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::EndOfEvent ERR: could not find primary";
0727 }
0728 
0729 void HcalTB04Analysis::qieAnalysis(CLHEP::HepRandomEngine* engine) {
0730   int hittot = hcalHitCache.size();
0731   if (hittot <= 0)
0732     hittot = 1;
0733   std::vector<CaloHit> hits(hittot);
0734   std::vector<int> todo(nTower, 0);
0735 
0736 #ifdef EDM_ML_DEBUG
0737   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::qieAnalysis: Size " << hits.size() << " " << todo.size() << " "
0738                                 << idTower.size() << " " << esimh.size() << " " << eqie.size();
0739 #endif
0740   // Loop over all HCal hits
0741   for (unsigned int k1 = 0; k1 < hcalHitCache.size(); k1++) {
0742     CaloHit hit = hcalHitCache[k1];
0743     uint32_t id = hit.id();
0744     int nhit = 0;
0745     double esim = hit.e();
0746     hits[nhit] = hit;
0747     for (unsigned int k2 = k1 + 1; k2 < hcalHitCache.size(); k2++) {
0748       hit = hcalHitCache[k2];
0749       if (hit.id() == id) {
0750         nhit++;
0751         hits[nhit] = hit;
0752         esim += hit.e();
0753       }
0754     }
0755     k1 += nhit;
0756     nhit++;
0757     std::vector<int> cd = myQie->getCode(nhit, hits, engine);
0758     double eq = myQie->getEnergy(cd);
0759 #ifdef EDM_ML_DEBUG
0760     edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::  ID 0x" << std::hex << id << std::dec << " registers " << esim
0761                                   << " energy from " << nhit << " hits starting with hit # " << k1
0762                                   << " energy with noise " << eq;
0763 #endif
0764     for (int k2 = 0; k2 < nTower; k2++) {
0765       if (id == idTower[k2]) {
0766         todo[k2] = 1;
0767         esimh[k2] = esim;
0768         eqie[k2] = eq;
0769       }
0770     }
0771   }
0772 
0773   // Towers with no hit
0774   for (int k2 = 0; k2 < nTower; k2++) {
0775     if (todo[k2] == 0) {
0776       std::vector<int> cd = myQie->getCode(0, hits, engine);
0777       double eq = myQie->getEnergy(cd);
0778       esimh[k2] = 0;
0779       eqie[k2] = eq;
0780 #ifdef EDM_ML_DEBUG
0781       edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::  ID 0x" << std::hex << idTower[k2] << std::dec
0782                                     << " registers " << esimh[k2] << " energy from hits and energy after QIE analysis "
0783                                     << eqie[k2];
0784 #endif
0785     }
0786   }
0787 }
0788 
0789 void HcalTB04Analysis::xtalAnalysis(CLHEP::HepRandomEngine* engine) {
0790   CLHEP::RandGaussQ randGauss(*engine);
0791 
0792   // Crystal Data
0793   std::vector<int> iok(nCrystal, 0);
0794 #ifdef EDM_ML_DEBUG
0795   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::xtalAnalysis: Size " << iok.size() << " " << idEcal.size() << " "
0796                                 << esime.size() << " " << enois.size();
0797 #endif
0798   for (unsigned int k1 = 0; k1 < ecalHitCache.size(); k1++) {
0799     uint32_t id = ecalHitCache[k1].id();
0800     int nhit = 0;
0801     double esim = ecalHitCache[k1].e();
0802     for (unsigned int k2 = k1 + 1; k2 < ecalHitCache.size(); k2++) {
0803       if (ecalHitCache[k2].id() == id) {
0804         nhit++;
0805         esim += ecalHitCache[k2].e();
0806       }
0807     }
0808     k1 += nhit;
0809     nhit++;
0810     double eq = esim + randGauss.fire(0., ecalNoise);
0811 #ifdef EDM_ML_DEBUG
0812     edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::  ID 0x" << std::hex << id << std::dec << " registers " << esim
0813                                   << " energy from " << nhit << " hits starting with hit # " << k1
0814                                   << " energy with noise " << eq;
0815 #endif
0816     for (int k2 = 0; k2 < nCrystal; k2++) {
0817       if (id == idEcal[k2]) {
0818         iok[k2] = 1;
0819         esime[k2] = esim;
0820         enois[k2] = eq;
0821       }
0822     }
0823   }
0824 
0825   // Crystals with no hit
0826   for (int k2 = 0; k2 < nCrystal; k2++) {
0827     if (iok[k2] == 0) {
0828       esime[k2] = 0;
0829       enois[k2] = randGauss.fire(0., ecalNoise);
0830 #ifdef EDM_ML_DEBUG
0831       edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::  ID 0x" << std::hex << idEcal[k2] << std::dec
0832                                     << " registers " << esime[k2] << " energy from hits and energy from noise "
0833                                     << enois[k2];
0834 #endif
0835     }
0836   }
0837 }
0838 
0839 void HcalTB04Analysis::finalAnalysis() {
0840   //Beam Information
0841   histo->fillPrimary(pInit, etaInit, phiInit);
0842 
0843   // Total Energy
0844   eecals = ehcals = eecalq = ehcalq = 0.;
0845   for (int i = 0; i < nTower; i++) {
0846     ehcals += esimh[i];
0847     ehcalq += eqie[i];
0848   }
0849   for (int i = 0; i < nCrystal; i++) {
0850     eecals += esime[i];
0851     eecalq += enois[i];
0852   }
0853   etots = eecals + ehcals;
0854   etotq = eecalq + ehcalq;
0855 #ifdef EDM_ML_DEBUG
0856   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Energy deposit at Sim Level (Total) " << etots << " (ECal) "
0857                                 << eecals << " (HCal) " << ehcals
0858                                 << "\nHcalTB04Analysis:: Energy deposit at Qie Level (Total) " << etotq << " (ECal) "
0859                                 << eecalq << " (HCal) " << ehcalq;
0860 #endif
0861   histo->fillEdep(etots, eecals, ehcals, etotq, eecalq, ehcalq);
0862 
0863   // Lateral Profile
0864   for (int i = 0; i < 5; i++) {
0865     eseta[i] = 0.;
0866     eqeta[i] = 0.;
0867   }
0868   for (int i = 0; i < 3; i++) {
0869     esphi[i] = 0.;
0870     eqphi[i] = 0.;
0871   }
0872   double e1 = 0, e2 = 0;
0873   unsigned int id;
0874   for (int i = 0; i < nTower; i++) {
0875     int det, z, group, ieta, iphi, layer;
0876     id = idTower[i];
0877     HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, layer);
0878     iphi -= (icphi - 1);
0879     if (icphi > 4) {
0880       if (ieta == 0)
0881         ieta = 2;
0882       else
0883         ieta = -1;
0884     } else {
0885       ieta = ieta - iceta + 2;
0886     }
0887     if (iphi >= 0 && iphi < 3 && ieta >= 0 && ieta < 5) {
0888       eseta[ieta] += esimh[i];
0889       esphi[iphi] += esimh[i];
0890       e1 += esimh[i];
0891       eqeta[ieta] += eqie[i];
0892       eqphi[iphi] += eqie[i];
0893       e2 += eqie[i];
0894     }
0895   }
0896   for (int i = 0; i < 3; i++) {
0897     if (e1 > 0)
0898       esphi[i] /= e1;
0899     if (e2 > 0)
0900       eqphi[i] /= e2;
0901   }
0902   for (int i = 0; i < 5; i++) {
0903     if (e1 > 0)
0904       eseta[i] /= e1;
0905     if (e2 > 0)
0906       eqeta[i] /= e2;
0907   }
0908 #ifdef EDM_ML_DEBUG
0909   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Energy fraction along Eta and Phi (Sim/Qie)";
0910   for (int i = 0; i < 5; i++)
0911     edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: [" << i << "] Eta Sim = " << eseta[i] << " Qie = " << eqeta[i]
0912                                   << " Phi Sim = " << esphi[i] << " Qie = " << eqphi[i];
0913 #endif
0914   histo->fillTrnsProf(eseta, eqeta, esphi, eqphi);
0915 
0916   // Longitudianl profile
0917   for (int i = 0; i < 20; i++) {
0918     eslay[i] = 0.;
0919     eqlay[i] = 0.;
0920   }
0921   e1 = 0;
0922   e2 = 0;
0923   for (int i = 0; i < nTower; i++) {
0924     int det, z, group, ieta, iphi, layer;
0925     id = idTower[i];
0926     HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, layer);
0927     iphi -= (icphi - 1);
0928     layer -= 1;
0929     if (iphi >= 0 && iphi < 3 && layer >= 0 && layer < 20) {
0930       eslay[layer] += esimh[i];
0931       e1 += esimh[i];
0932       eqlay[layer] += eqie[i];
0933       e2 += eqie[i];
0934     }
0935   }
0936   for (int i = 0; i < 20; i++) {
0937     if (e1 > 0)
0938       eslay[i] /= e1;
0939     if (e2 > 0)
0940       eqlay[i] /= e2;
0941   }
0942 #ifdef EDM_ML_DEBUG
0943   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Energy fraction along Layer";
0944   for (int i = 0; i < 20; i++)
0945     edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: [" << i << "] Sim = " << eslay[i] << " Qie = " << eqlay[i];
0946 #endif
0947   histo->fillLongProf(eslay, eqlay);
0948 }
0949 
0950 void HcalTB04Analysis::fillEvent(PHcalTB04Info& product) {
0951   //Setup the ID's
0952   product.setIDs(idHcal, idXtal);
0953 
0954   //Beam Information
0955   product.setPrimary(nPrimary, particleType, pInit, etaInit, phiInit);
0956 
0957   //Energy deposits in the crystals and towers
0958   product.setEdepHcal(esimh, eqie);
0959   product.setEdepHcal(esime, enois);
0960 
0961   // Total Energy
0962   product.setEdep(etots, eecals, ehcals, etotq, eecalq, ehcalq);
0963 
0964   // Lateral Profile
0965   product.setTrnsProf(eseta, eqeta, esphi, eqphi);
0966 
0967   // Longitudianl profile
0968   product.setLongProf(eslay, eqlay);
0969 
0970   //Save Hits
0971   int i, nhit = 0;
0972   std::vector<CaloHit>::iterator itr;
0973   for (i = 0, itr = ecalHitCache.begin(); itr != ecalHitCache.end(); i++, itr++) {
0974     uint32_t id = itr->id();
0975     int det, z, group, ieta, iphi, lay;
0976     HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, lay);
0977     product.saveHit(det, lay, ieta, iphi, itr->e(), itr->t());
0978     nhit++;
0979 #ifdef EDM_ML_DEBUG
0980     edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Save Hit " << std::setw(3) << i + 1 << " ID 0x" << std::hex
0981                                   << group << std::dec << " " << std::setw(2) << det << " " << std::setw(2) << lay
0982                                   << " " << std::setw(1) << z << " " << std::setw(3) << ieta << " " << std::setw(3)
0983                                   << iphi << " T " << std::setw(6) << itr->t() << " E " << std::setw(6) << itr->e();
0984 #endif
0985   }
0986   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Saves " << nhit << " hits from Crystals";
0987   int hit = nhit;
0988   nhit = 0;
0989 
0990   for (i = hit, itr = hcalHitCache.begin(); itr != hcalHitCache.end(); i++, itr++) {
0991     uint32_t id = itr->id();
0992     int det, z, group, ieta, iphi, lay;
0993     HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, lay);
0994     product.saveHit(det, lay, ieta, iphi, itr->e(), itr->t());
0995     nhit++;
0996 #ifdef EDM_ML_DEBUG
0997     edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Save Hit " << std::setw(3) << i + 1 << " ID 0x" << std::hex
0998                                   << group << std::dec << " " << std::setw(2) << det << " " << std::setw(2) << lay
0999                                   << " " << std::setw(1) << z << " " << std::setw(3) << ieta << " " << std::setw(3)
1000                                   << iphi << " T " << std::setw(6) << itr->t() << " E " << std::setw(6) << itr->e();
1001 #endif
1002   }
1003   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Saves " << nhit << " hits from HCal";
1004 
1005   //Vertex associated quantities
1006   product.setVtxPrim(evNum,
1007                      pvType,
1008                      pvPosition.x(),
1009                      pvPosition.y(),
1010                      pvPosition.z(),
1011                      pvUVW.x(),
1012                      pvUVW.y(),
1013                      pvUVW.z(),
1014                      pvMomentum.x(),
1015                      pvMomentum.y(),
1016                      pvMomentum.z());
1017   for (unsigned int i = 0; i < secTrackID.size(); i++) {
1018     product.setVtxSec(
1019         secTrackID[i], secPartID[i], secMomentum[i].x(), secMomentum[i].y(), secMomentum[i].z(), secEkin[i]);
1020   }
1021 }
1022 
1023 void HcalTB04Analysis::clear() {
1024   pvFound = false;
1025   pvType = -2;
1026   pvPosition = G4ThreeVector();
1027   pvMomentum = G4ThreeVector();
1028   pvUVW = G4ThreeVector();
1029   secTrackID.clear();
1030   secPartID.clear();
1031   secMomentum.clear();
1032   secEkin.clear();
1033   shortLivedSecondaries.clear();
1034 
1035   ecalHitCache.erase(ecalHitCache.begin(), ecalHitCache.end());
1036   hcalHitCache.erase(hcalHitCache.begin(), hcalHitCache.end());
1037   hcalHitLayer.erase(hcalHitLayer.begin(), hcalHitLayer.end());
1038   nPrimary = particleType = 0;
1039   pInit = etaInit = phiInit = 0;
1040 
1041   esimh.clear();
1042   eqie.clear();
1043   esimh.reserve(nTower);
1044   eqie.reserve(nTower);
1045   for (int i = 0; i < nTower; i++) {
1046     esimh.push_back(0.);
1047     eqie.push_back(0.);
1048   }
1049   esime.clear();
1050   enois.clear();
1051   esime.reserve(nCrystal);
1052   enois.reserve(nCrystal);
1053   for (int i = 0; i < nCrystal; i++) {
1054     esime.push_back(0.);
1055     enois.push_back(0.);
1056   }
1057 }
1058 
1059 int HcalTB04Analysis::unitID(uint32_t id) {
1060   int det, z, group, ieta, iphi, lay;
1061   HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, lay);
1062   group = (det & 15) << 20;
1063   group += ((lay - 1) & 31) << 15;
1064   group += (z & 1) << 14;
1065   group += (ieta & 127) << 7;
1066   group += (iphi & 127);
1067   return group;
1068 }
1069 
1070 double HcalTB04Analysis::scale(int det, int layer) {
1071   double tmp = 1.;
1072   if (det == static_cast<int>(HcalBarrel)) {
1073     if (layer == 1)
1074       tmp = scaleHB0;
1075     else if (layer == 17)
1076       tmp = scaleHB16;
1077     else if (layer > 17)
1078       tmp = scaleHO;
1079   } else {
1080     if (layer <= 2)
1081       tmp = scaleHE0;
1082   }
1083   return tmp;
1084 }
1085 
1086 double HcalTB04Analysis::timeOfFlight(int det, int layer, double eta) {
1087   double theta = 2.0 * std::atan(std::exp(-eta));
1088   double dist = beamOffset;
1089   if (det == static_cast<int>(HcalBarrel)) {
1090     const double rLay[19] = {1836.0,
1091                              1902.0,
1092                              1962.0,
1093                              2022.0,
1094                              2082.0,
1095                              2142.0,
1096                              2202.0,
1097                              2262.0,
1098                              2322.0,
1099                              2382.0,
1100                              2448.0,
1101                              2514.0,
1102                              2580.0,
1103                              2646.0,
1104                              2712.0,
1105                              2776.0,
1106                              2862.5,
1107                              3847.0,
1108                              4052.0};
1109     if (layer > 0 && layer <= 19)
1110       dist += rLay[layer - 1] * mm / sin(theta);
1111   } else {
1112     const double zLay[19] = {4034.0,
1113                              4032.0,
1114                              4123.0,
1115                              4210.0,
1116                              4297.0,
1117                              4384.0,
1118                              4471.0,
1119                              4558.0,
1120                              4645.0,
1121                              4732.0,
1122                              4819.0,
1123                              4906.0,
1124                              4993.0,
1125                              5080.0,
1126                              5167.0,
1127                              5254.0,
1128                              5341.0,
1129                              5428.0,
1130                              5515.0};
1131     if (layer > 0 && layer <= 19)
1132       dist += zLay[layer - 1] * mm / cos(theta);
1133   }
1134 
1135   double tmp = dist / c_light / ns;
1136 #ifdef EDM_ML_DEBUG
1137   edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::timeOfFlight " << tmp << " for det/lay " << det << " " << layer
1138                                 << " eta/theta " << eta << " " << theta / deg << " dist " << dist;
1139 #endif
1140   return tmp;
1141 }
1142 
1143 DEFINE_SIMWATCHER(HcalTB04Analysis);
[ 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!