Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-05-10 02:21:21 

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