Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-01-16 23:37:22 

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