Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0002 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0003 #include "SimG4Core/Notification/interface/EndOfEvent.h"
0004 #include "SimG4Core/Notification/interface/Observer.h"
0005 #include "SimG4Core/Watcher/interface/SimProducer.h"
0006 
0007 // to retreive hits
0008 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0009 #include "DataFormats/Math/interface/Point3D.h"
0010 #include "SimDataFormats/CaloHit/interface/CaloHit.h"
0011 #include "SimDataFormats/CaloTest/interface/HcalTestHistoClass.h"
0012 
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/Framework/interface/EventSetup.h"
0015 
0016 #include "Geometry/HcalCommonData/interface/HcalNumberingFromDDD.h"
0017 #include "Geometry/HcalCommonData/interface/HcalDDDSimConstants.h"
0018 #include "Geometry/Records/interface/HcalSimNumberingRecord.h"
0019 
0020 #include "SimG4CMS/Calo/interface/CaloG4Hit.h"
0021 #include "SimG4CMS/Calo/interface/CaloG4HitCollection.h"
0022 #include "SimG4CMS/Calo/interface/HCalSD.h"
0023 #include "SimG4CMS/Calo/interface/HcalQie.h"
0024 #include "SimG4CMS/Calo/interface/HcalTestNumberingScheme.h"
0025 
0026 #include "G4SDManager.hh"
0027 #include "G4Step.hh"
0028 #include "G4Track.hh"
0029 #include "G4VProcess.hh"
0030 #include "G4HCofThisEvent.hh"
0031 
0032 #include <CLHEP/Random/Randomize.h>
0033 #include <CLHEP/Units/SystemOfUnits.h>
0034 #include <CLHEP/Units/PhysicalConstants.h>
0035 
0036 #include <cmath>
0037 #include <iomanip>
0038 #include <iostream>
0039 #include <string>
0040 #include <vector>
0041 
0042 using CLHEP::c_light;
0043 using CLHEP::deg;
0044 using CLHEP::GeV;
0045 using CLHEP::MeV;
0046 using CLHEP::mm;
0047 using CLHEP::ns;
0048 
0049 class HcalTestAnalysis : public SimProducer,
0050                          public Observer<const BeginOfRun*>,
0051                          public Observer<const BeginOfEvent*>,
0052                          public Observer<const EndOfEvent*>,
0053                          public Observer<const G4Step*> {
0054 public:
0055   HcalTestAnalysis(const edm::ParameterSet& p);
0056   ~HcalTestAnalysis() override;
0057 
0058   void registerConsumes(edm::ConsumesCollector) override;
0059   void produce(edm::Event&, const edm::EventSetup&) override;
0060   void beginRun(edm::EventSetup const&) override;
0061 
0062 private:
0063   // observer classes
0064   void update(const BeginOfRun* run) override;
0065   void update(const BeginOfEvent* evt) override;
0066   void update(const EndOfEvent* evt) override;
0067   void update(const G4Step* step) override;
0068 
0069   // analysis-related stuff
0070   std::vector<int> layerGrouping(int);
0071   std::vector<int> towersToAdd(int centre, int nadd);
0072   void fill(const EndOfEvent* ev);
0073   void qieAnalysis(CLHEP::HepRandomEngine*);
0074   void layerAnalysis();
0075   double timeOfFlight(int det, int layer, double eta);
0076 
0077 private:
0078   // Qie Analysis
0079   std::unique_ptr<HcalQie> myqie_;
0080   int addTower_;
0081 
0082   // Private Tuples
0083   std::unique_ptr<HcalTestHistoClass> tuples_;
0084 
0085   // to read from ParameterSet
0086   const edm::ParameterSet m_Anal;
0087   double eta0_, phi0_;
0088   const int laygroup_, centralTower_;
0089   const std::vector<std::string> names_;
0090   //Keep parameters to instantiate HcalTestHistoClass later
0091   const std::string fileName_;
0092 
0093   // Numbering scheme
0094   edm::ESGetToken<HcalDDDSimConstants, HcalSimNumberingRecord> ddconsToken_;
0095   std::unique_ptr<HcalNumberingFromDDD> numberingFromDDD_;
0096   const HcalDDDSimConstants* hcons_;
0097   std::unique_ptr<HcalTestNumberingScheme> org_;
0098 
0099   // Hits for qie analysis
0100   std::vector<CaloHit> caloHitCache_;
0101   std::vector<int> group_, tower_;
0102   int nGroup_, nTower_;
0103 
0104   // some private members for ananlysis
0105   unsigned int count_;
0106   double edepEB_, edepEE_, edepHB_, edepHE_;
0107   double edepHO_, edepl_[20];
0108   double mudist_[20];  // Distance of muon from central part
0109 };
0110 
0111 HcalTestAnalysis::HcalTestAnalysis(const edm::ParameterSet& p)
0112     : addTower_(3),
0113       m_Anal(p.getParameter<edm::ParameterSet>("HcalTestAnalysis")),
0114       eta0_(m_Anal.getParameter<double>("Eta0")),
0115       phi0_(m_Anal.getParameter<double>("Phi0")),
0116       laygroup_(m_Anal.getParameter<int>("LayerGrouping")),
0117       centralTower_(m_Anal.getParameter<int>("CentralTower")),
0118       names_(m_Anal.getParameter<std::vector<std::string> >("Names")),
0119       fileName_(m_Anal.getParameter<std::string>("FileName")),
0120       hcons_(nullptr) {
0121   org_.reset(nullptr);
0122 
0123   tuples_.reset(nullptr);
0124   numberingFromDDD_.reset(nullptr);
0125   edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: Initialised as observer of begin/end events"
0126                               << " and of G4step";
0127 
0128   count_ = 0;
0129   group_ = layerGrouping(laygroup_);
0130   nGroup_ = 0;
0131   for (unsigned int i = 0; i < group_.size(); i++)
0132     if (group_[i] > nGroup_)
0133       nGroup_ = group_[i];
0134   tower_ = towersToAdd(centralTower_, addTower_);
0135   nTower_ = tower_.size() / 2;
0136 
0137   edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: initialised for " << nGroup_ << " Longitudinal groups and "
0138                               << nTower_ << " towers";
0139 
0140   // qie
0141   myqie_ = std::make_unique<HcalQie>(p);
0142 
0143   produces<HcalTestHistoClass>();
0144 }
0145 
0146 HcalTestAnalysis::~HcalTestAnalysis() {
0147   edm::LogVerbatim("HcalSim") << "HcalTestAnalysis: -------->  Total number of selected entries : " << count_;
0148   edm::LogVerbatim("HcalSim") << "HcalTestAnalysis: Pointers:: Numbering Scheme " << org_.get();
0149 }
0150 
0151 void HcalTestAnalysis::registerConsumes(edm::ConsumesCollector cc) {
0152   ddconsToken_ = cc.esConsumes<HcalDDDSimConstants, HcalSimNumberingRecord, edm::Transition::BeginRun>();
0153   edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::Initialize ESGetToken for HcalDDDSimConstants";
0154 }
0155 
0156 void HcalTestAnalysis::produce(edm::Event& e, const edm::EventSetup&) {
0157   e.put(std::move(tuples_));
0158   tuples_.reset(nullptr);
0159 }
0160 
0161 std::vector<int> HcalTestAnalysis::layerGrouping(int group) {
0162   std::vector<int> temp(19);
0163   if (group <= 1) {
0164     int grp[19] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2};
0165     for (int i = 0; i < 19; i++)
0166       temp[i] = grp[i];
0167   } else if (group == 2) {
0168     int grp[19] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19};
0169     for (int i = 0; i < 19; i++)
0170       temp[i] = grp[i];
0171   } else if (group == 3) {
0172     int grp[19] = {1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7};
0173     for (int i = 0; i < 19; i++)
0174       temp[i] = grp[i];
0175   } else if (group == 4) {
0176     int grp[19] = {1, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7};
0177     for (int i = 0; i < 19; i++)
0178       temp[i] = grp[i];
0179   } else {
0180     int grp[19] = {1, 1, 2, 3, 3, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 7, 7};
0181     for (int i = 0; i < 19; i++)
0182       temp[i] = grp[i];
0183   }
0184 
0185   edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: Layer Grouping ";
0186   for (int i = 0; i < 19; i++)
0187     edm::LogVerbatim("HcalSim") << "HcalTestAnalysis: Group[" << i << "] = " << temp[i];
0188   return temp;
0189 }
0190 
0191 std::vector<int> HcalTestAnalysis::towersToAdd(int centre, int nadd) {
0192   int etac = (centre / 100) % 100;
0193   int phic = (centre % 100);
0194   int etamin, etamax, phimin, phimax;
0195   if (etac > 0) {
0196     etamin = etac - nadd;
0197     etamax = etac + nadd;
0198   } else {
0199     etamin = etac;
0200     etamax = etac;
0201   }
0202   if (phic > 0) {
0203     phimin = phic - nadd;
0204     phimax = phic + nadd;
0205   } else {
0206     phimin = phic;
0207     phimax = phic;
0208   }
0209 
0210   int nbuf, kount = 0;
0211   nbuf = (etamax - etamin + 1) * (phimax - phimin + 1);
0212   std::vector<int> temp(2 * nbuf);
0213   for (int eta = etamin; eta <= etamax; eta++) {
0214     for (int phi = phimin; phi <= phimax; phi++) {
0215       temp[kount] = (eta * 100 + phi);
0216       temp[kount + nbuf] = std::max(abs(eta - etac), abs(phi - phic));
0217       kount++;
0218     }
0219   }
0220 
0221   edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: Towers to be considered for Central " << centre << " and " << nadd
0222                               << " on either side";
0223   for (int i = 0; i < nbuf; i++)
0224     edm::LogVerbatim("HcalSim") << "HcalTestAnalysis: Tower[" << std::setw(3) << i << "] " << temp[i] << " "
0225                                 << temp[nbuf + i];
0226   return temp;
0227 }
0228 
0229 //==================================================================== per JOB
0230 
0231 void HcalTestAnalysis::beginRun(edm::EventSetup const& es) {
0232   // Numbering From DDD
0233   hcons_ = &es.getData(ddconsToken_);
0234   edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: Initialise HcalNumberingFromDDD for " << names_[0];
0235   numberingFromDDD_ = std::make_unique<HcalNumberingFromDDD>(hcons_);
0236 
0237   // Numbering scheme
0238   org_ = std::make_unique<HcalTestNumberingScheme>(false);
0239 }
0240 
0241 //==================================================================== per RUN
0242 void HcalTestAnalysis::update(const BeginOfRun* run) {
0243   int irun = (*run)()->GetRunID();
0244   edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: Begin of Run = " << irun;
0245 
0246   bool loop = true, eta = true, phi = true;
0247   int etac = (centralTower_ / 100) % 100;
0248   if (etac == 0) {
0249     etac = 1;
0250     eta = false;
0251   }
0252   int phic = (centralTower_ % 100);
0253   if (phic == 0) {
0254     phic = 1;
0255     phi = false;
0256   }
0257   int idet = static_cast<int>(HcalBarrel);
0258   while (loop) {
0259     HcalCellType::HcalCell tmp = hcons_->cell(idet, 1, 1, etac, phic);
0260     if (tmp.ok) {
0261       if (eta)
0262         eta0_ = tmp.eta;
0263       if (phi)
0264         phi0_ = tmp.phi;
0265       loop = false;
0266     } else if (idet == static_cast<int>(HcalBarrel)) {
0267       idet = static_cast<int>(HcalEndcap);
0268     } else if (idet == static_cast<int>(HcalEndcap)) {
0269       idet = static_cast<int>(HcalForward);
0270     } else {
0271       loop = false;
0272     }
0273   }
0274 
0275   edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: Central Tower " << centralTower_
0276                               << " corresponds to eta0 = " << eta0_ << " phi0 = " << phi0_;
0277 
0278   std::string sdname = names_[0];
0279   G4SDManager* sd = G4SDManager::GetSDMpointerIfExist();
0280   if (sd != nullptr) {
0281     G4VSensitiveDetector* aSD = sd->FindSensitiveDetector(sdname);
0282     if (aSD == nullptr) {
0283       edm::LogWarning("HcalSim") << "HcalTestAnalysis::beginOfRun: No SD with "
0284                                  << "name " << sdname << " in this Setup";
0285     } else {
0286       HCalSD* theCaloSD = dynamic_cast<HCalSD*>(aSD);
0287       edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::beginOfRun: Finds SD with name " << theCaloSD->GetName()
0288                                   << " in this Setup";
0289       if (org_.get()) {
0290         theCaloSD->setNumberingScheme(org_.get());
0291         edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::beginOfRun: set a new numbering scheme";
0292       }
0293     }
0294   } else {
0295     edm::LogWarning("HcalSim") << "HcalTestAnalysis::beginOfRun: Could not get"
0296                                << " SD Manager!";
0297   }
0298 }
0299 
0300 //=================================================================== per EVENT
0301 void HcalTestAnalysis::update(const BeginOfEvent* evt) {
0302   // create tuple object
0303   tuples_ = std::make_unique<HcalTestHistoClass>();
0304   // Reset counters
0305   tuples_->setCounters();
0306 
0307   int i = 0;
0308   edepEB_ = edepEE_ = edepHB_ = edepHE_ = edepHO_ = 0.;
0309   for (i = 0; i < 20; i++)
0310     edepl_[i] = 0.;
0311   for (i = 0; i < 20; i++)
0312     mudist_[i] = -1.;
0313 
0314   edm::LogVerbatim("HcalSim") << "HcalTestAnalysis: Begin of event = " << (*evt)()->GetEventID();
0315 }
0316 
0317 //=================================================================== each STEP
0318 void HcalTestAnalysis::update(const G4Step* aStep) {
0319   if (aStep != nullptr) {
0320     G4VPhysicalVolume* curPV = aStep->GetPreStepPoint()->GetPhysicalVolume();
0321     G4String name = curPV->GetName();
0322     name.assign(name, 0, 3);
0323     double edeposit = aStep->GetTotalEnergyDeposit();
0324     int layer = -1;
0325     if (name == "EBR") {
0326       edepEB_ += edeposit;
0327     } else if (name == "EFR") {
0328       edepEE_ += edeposit;
0329     } else if (name == "HBS") {
0330       layer = (curPV->GetCopyNo() / 10) % 100;
0331       if (layer >= 0 && layer < 17) {
0332         edepHB_ += edeposit;
0333       } else {
0334         edm::LogWarning("HcalSim") << "HcalTestAnalysis::Error in HB " << curPV->GetName() << curPV->GetCopyNo();
0335         layer = -1;
0336       }
0337     } else if (name == "HES") {
0338       layer = (curPV->GetCopyNo() / 10) % 100;
0339       if (layer >= 0 && layer < 19) {
0340         edepHE_ += edeposit;
0341       } else {
0342         edm::LogWarning("HcalSim") << "HcalTestAnalysis::Error in HE " << curPV->GetName() << curPV->GetCopyNo();
0343         layer = -1;
0344       }
0345     } else if (name == "HTS") {
0346       layer = (curPV->GetCopyNo() / 10) % 100;
0347       if (layer >= 17 && layer < 20) {
0348         edepHO_ += edeposit;
0349       } else {
0350         edm::LogWarning("HcalSim") << "HcalTestAnalysis::Error in HO " << curPV->GetName() << curPV->GetCopyNo();
0351         layer = -1;
0352       }
0353     }
0354     if (layer >= 0 && layer < 20) {
0355       edepl_[layer] += edeposit;
0356 
0357       // Calculate the distance if it is a muon
0358       G4String part = aStep->GetTrack()->GetDefinition()->GetParticleName();
0359       if ((part == "mu-" || part == "mu+") && mudist_[layer] < 0) {
0360         math::XYZPoint pos(aStep->GetPreStepPoint()->GetPosition().x(),
0361                            aStep->GetPreStepPoint()->GetPosition().y(),
0362                            aStep->GetPreStepPoint()->GetPosition().z());
0363         double theta = pos.theta();
0364         double eta = -log(tan(theta * 0.5));
0365         double phi = pos.phi();
0366         double dist = sqrt((eta - eta0_) * (eta - eta0_) + (phi - phi0_) * (phi - phi0_));
0367         mudist_[layer] = dist * std::sqrt(pos.perp2());
0368       }
0369     }
0370 
0371     if (layer >= 0 && layer < 20) {
0372       edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: G4Step: " << name << " Layer " << std::setw(3) << layer
0373                                   << " Edep " << std::setw(6) << edeposit / MeV << " MeV";
0374     }
0375   } else {
0376     edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: G4Step: Null Step";
0377   }
0378 }
0379 
0380 //================================================================ End of EVENT
0381 void HcalTestAnalysis::update(const EndOfEvent* evt) {
0382   ++count_;
0383   // Fill event input
0384   fill(evt);
0385   edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: ---  after Fill";
0386 
0387   // Qie analysis
0388   CLHEP::HepRandomEngine* engine = G4Random::getTheEngine();
0389   qieAnalysis(engine);
0390   edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: ---  after QieAnalysis";
0391 
0392   // Layers tuples filling
0393   layerAnalysis();
0394   edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: ---  after LayerAnalysis";
0395 }
0396 
0397 //---------------------------------------------------
0398 void HcalTestAnalysis::fill(const EndOfEvent* evt) {
0399   edm::LogVerbatim("HcalSim") << "HcalTestAnalysis: Fill event " << (*evt)()->GetEventID();
0400 
0401   // access to the G4 hit collections
0402   G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
0403 
0404   int nhc = 0, neb = 0, nef = 0, j = 0;
0405   caloHitCache_.erase(caloHitCache_.begin(), caloHitCache_.end());
0406 
0407   // Hcal
0408   int HCHCid = G4SDManager::GetSDMpointer()->GetCollectionID(names_[0]);
0409   CaloG4HitCollection* theHCHC = (CaloG4HitCollection*)allHC->GetHC(HCHCid);
0410   edm::LogVerbatim("HcalSim") << "HcalTestAnalysis :: Hit Collection for " << names_[0] << " of ID " << HCHCid
0411                               << " is obtained at " << theHCHC;
0412   int hchc_entries = theHCHC->entries();
0413   if (HCHCid >= 0 && theHCHC != nullptr) {
0414     for (j = 0; j < hchc_entries; j++) {
0415       CaloG4Hit* aHit = (*theHCHC)[j];
0416 
0417       double e = aHit->getEnergyDeposit() / GeV;
0418       double time = aHit->getTimeSlice();
0419 
0420       math::XYZPoint pos = aHit->getPosition();
0421       double theta = pos.theta();
0422       double eta = -log(tan(theta * 0.5));
0423       double phi = pos.phi();
0424 
0425       uint32_t unitID = aHit->getUnitID();
0426       int subdet, zside, layer, etaIndex, phiIndex, lay;
0427       org_->unpackHcalIndex(unitID, subdet, zside, layer, etaIndex, phiIndex, lay);
0428       double jitter = time - timeOfFlight(subdet, lay, eta);
0429       if (jitter < 0)
0430         jitter = 0;
0431       CaloHit hit(subdet, lay, e, eta, phi, jitter, unitID);
0432       caloHitCache_.push_back(hit);
0433       nhc++;
0434 
0435       std::string det = "HB";
0436       if (subdet == static_cast<int>(HcalForward)) {
0437         det = "HF";
0438       } else if (subdet == static_cast<int>(HcalEndcap)) {
0439         if (etaIndex <= 20) {
0440           det = "HES";
0441         } else {
0442           det = "HED";
0443         }
0444       }
0445 
0446       edm::LogVerbatim("HcalSim") << "HcalTest: " << det << "  layer " << std::setw(2) << layer << " time "
0447                                   << std::setw(6) << time << " theta " << std::setw(8) << theta << " eta "
0448                                   << std::setw(8) << eta << " phi " << std::setw(8) << phi << " e " << std::setw(8)
0449                                   << e;
0450     }
0451   }
0452   edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::HCAL hits : " << nhc;
0453 
0454   // EB
0455   int EBHCid = G4SDManager::GetSDMpointer()->GetCollectionID(names_[1]);
0456   CaloG4HitCollection* theEBHC = (CaloG4HitCollection*)allHC->GetHC(EBHCid);
0457   edm::LogVerbatim("HcalSim") << "HcalTestAnalysis :: Hit Collection for " << names_[1] << " of ID " << EBHCid
0458                               << " is obtained at " << theEBHC;
0459   int ebhc_entries = theEBHC->entries();
0460   if (EBHCid >= 0 && theEBHC != nullptr) {
0461     for (j = 0; j < ebhc_entries; j++) {
0462       CaloG4Hit* aHit = (*theEBHC)[j];
0463 
0464       double e = aHit->getEnergyDeposit() / GeV;
0465       double time = aHit->getTimeSlice();
0466       std::string det = "EB";
0467 
0468       math::XYZPoint pos = aHit->getPosition();
0469       double theta = pos.theta();
0470       double eta = -log(tan(theta / 2.));
0471       double phi = pos.phi();
0472 
0473       HcalNumberingFromDDD::HcalID id = numberingFromDDD_->unitID(eta, phi, 1, 1);
0474       uint32_t unitID = org_->getUnitID(id);
0475       int subdet, zside, layer, ieta, iphi, lay;
0476       org_->unpackHcalIndex(unitID, subdet, zside, layer, ieta, iphi, lay);
0477       subdet = 10;
0478       layer = 0;
0479       unitID = org_->packHcalIndex(subdet, zside, layer, ieta, iphi, lay);
0480       CaloHit hit(subdet, lay, e, eta, phi, time, unitID);
0481       caloHitCache_.push_back(hit);
0482       neb++;
0483       edm::LogVerbatim("HcalSim") << "HcalTest: " << det << "  layer " << std::setw(2) << layer << " time "
0484                                   << std::setw(6) << time << " theta " << std::setw(8) << theta << " eta "
0485                                   << std::setw(8) << eta << " phi " << std::setw(8) << phi << " e " << std::setw(8)
0486                                   << e;
0487     }
0488   }
0489   edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::EB hits : " << neb;
0490 
0491   // EE
0492   int EEHCid = G4SDManager::GetSDMpointer()->GetCollectionID(names_[2]);
0493   CaloG4HitCollection* theEEHC = (CaloG4HitCollection*)allHC->GetHC(EEHCid);
0494   edm::LogVerbatim("HcalSim") << "HcalTestAnalysis :: Hit Collection for " << names_[2] << " of ID " << EEHCid
0495                               << " is obtained at " << theEEHC;
0496   int eehc_entries = theEEHC->entries();
0497   if (EEHCid >= 0 && theEEHC != nullptr) {
0498     for (j = 0; j < eehc_entries; j++) {
0499       CaloG4Hit* aHit = (*theEEHC)[j];
0500 
0501       double e = aHit->getEnergyDeposit() / GeV;
0502       double time = aHit->getTimeSlice();
0503       std::string det = "EE";
0504 
0505       math::XYZPoint pos = aHit->getPosition();
0506       double theta = pos.theta();
0507       double eta = -log(tan(theta / 2.));
0508       double phi = pos.phi();
0509 
0510       HcalNumberingFromDDD::HcalID id = numberingFromDDD_->unitID(eta, phi, 1, 1);
0511       uint32_t unitID = org_->getUnitID(id);
0512       int subdet, zside, layer, ieta, iphi, lay;
0513       org_->unpackHcalIndex(unitID, subdet, zside, layer, ieta, iphi, lay);
0514       subdet = 11;
0515       layer = 0;
0516       unitID = org_->packHcalIndex(subdet, zside, layer, ieta, iphi, lay);
0517       CaloHit hit(subdet, lay, e, eta, phi, time, unitID);
0518       caloHitCache_.push_back(hit);
0519       nef++;
0520       edm::LogVerbatim("HcalSim") << "HcalTest: " << det << "  layer " << std::setw(2) << layer << " time "
0521                                   << std::setw(6) << time << " theta " << std::setw(8) << theta << " eta "
0522                                   << std::setw(8) << eta << " phi " << std::setw(8) << phi << " e " << std::setw(8)
0523                                   << e;
0524     }
0525   }
0526   edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::EE hits : " << nef;
0527 }
0528 
0529 //-----------------------------------------------------------------------------
0530 void HcalTestAnalysis::qieAnalysis(CLHEP::HepRandomEngine* engine) {
0531   //Fill tuple with hit information
0532   int hittot = caloHitCache_.size();
0533   tuples_->fillHits(caloHitCache_);
0534 
0535   //Get the index of the central tower
0536   HcalNumberingFromDDD::HcalID id = numberingFromDDD_->unitID(eta0_, phi0_, 1, 1);
0537   uint32_t unitID = org_->getUnitID(id);
0538   int subdet, zside, layer, ieta, iphi, lay;
0539   org_->unpackHcalIndex(unitID, subdet, zside, layer, ieta, iphi, lay);
0540   int laymax = 0;
0541   std::string det = "Unknown";
0542   if (subdet == static_cast<int>(HcalBarrel)) {
0543     laymax = 4;
0544     det = "HB";
0545   } else if (subdet == static_cast<int>(HcalEndcap)) {
0546     laymax = 2;
0547     det = "HES";
0548   }
0549   edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::Qie: " << det << " Eta " << ieta << " Phi " << iphi << " Laymax "
0550                               << laymax << " Hits " << hittot;
0551 
0552   if (laymax > 0 && hittot > 0) {
0553     std::vector<CaloHit> hits(hittot);
0554     std::vector<double> eqielay(80, 0.0), esimlay(80, 0.0), esimtot(4, 0.0);
0555     std::vector<double> eqietow(200, 0.0), esimtow(200, 0.0), eqietot(4, 0.0);
0556     int etac = (centralTower_ / 100) % 100;
0557     int phic = (centralTower_ % 100);
0558 
0559     for (int layr = 0; layr < nGroup_; layr++) {
0560       /*
0561       int layx, layy=20;
0562       for (int i=0; i<20; i++) 
0563     if (group_[i] == layr+1 && i < layy) layy = i+1;
0564       if (subdet == static_cast<int>(HcalBarrel)) {
0565     if      (layy < 2)    layx = 0;
0566     else if (layy < 17)   layx = 1;
0567     else if (layy == 17)  layx = 2;
0568     else                  layx = 3;
0569       } else {
0570     if      (layy < 2)    layx = 0;
0571     else                  layx = 1;
0572       }
0573       */
0574       for (int it = 0; it < nTower_; it++) {
0575         int nhit = 0;
0576         double esim = 0;
0577         for (int k1 = 0; k1 < hittot; k1++) {
0578           CaloHit hit = caloHitCache_[k1];
0579           int subdetc = hit.det();
0580           int layer = hit.layer();
0581           int group = 0;
0582           if (layer > 0 && layer < 20)
0583             group = group_[layer];
0584           if (subdetc == subdet && group == layr + 1) {
0585             int zsidec, ietac, iphic, idx;
0586             unitID = hit.id();
0587             org_->unpackHcalIndex(unitID, subdetc, zsidec, layer, ietac, iphic, lay);
0588             if (etac > 0 && phic > 0) {
0589               idx = ietac * 100 + iphic;
0590             } else if (etac > 0) {
0591               idx = ietac * 100;
0592             } else if (phic > 0) {
0593               idx = iphic;
0594             } else {
0595               idx = 0;
0596             }
0597             if (zsidec == zside && idx == tower_[it]) {
0598               hits[nhit] = hit;
0599               edm::LogVerbatim("HcalSim") << "HcalTest: Hit " << nhit << " " << hit;
0600               nhit++;
0601               esim += hit.e();
0602             }
0603           }
0604         }
0605 
0606         std::vector<int> cd = myqie_->getCode(nhit, hits, engine);
0607         double eqie = myqie_->getEnergy(cd);
0608 
0609         edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::Qie: Energy in layer " << layr << " Sim " << esim
0610                                     << " After QIE " << eqie;
0611         for (int i = 0; i < 4; i++) {
0612           if (tower_[nTower_ + it] <= i) {
0613             esimtot[i] += esim;
0614             eqietot[i] += eqie;
0615             esimlay[20 * i + layr] += esim;
0616             eqielay[20 * i + layr] += eqie;
0617             esimtow[50 * i + it] += esim;
0618             eqietow[50 * i + it] += eqie;
0619           }
0620         }
0621       }
0622     }
0623     edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::Qie: Total energy " << esimtot[3] << " (SimHit) " << eqietot[3]
0624                                 << " (After QIE)";
0625 
0626     std::vector<double> latphi(10);
0627     int nt = 2 * addTower_ + 1;
0628     for (int it = 0; it < nt; it++)
0629       latphi[it] = it - addTower_;
0630     for (int i = 0; i < 4; i++) {
0631       double scals = 1, scalq = 1;
0632       std::vector<double> latfs(10, 0.), latfq(10, 0.), longs(20), longq(20);
0633       if (esimtot[i] > 0)
0634         scals = 1. / esimtot[i];
0635       if (eqietot[i] > 0)
0636         scalq = 1. / eqietot[i];
0637       for (int it = 0; it < nTower_; it++) {
0638         int phib = it % nt;
0639         latfs[phib] += scals * esimtow[50 * i + it];
0640         latfq[phib] += scalq * eqietow[50 * i + it];
0641       }
0642       for (int layr = 0; layr <= nGroup_; layr++) {
0643         longs[layr] = scals * esimlay[20 * i + layr];
0644         longq[layr] = scalq * eqielay[20 * i + layr];
0645       }
0646       tuples_->fillQie(i, esimtot[i], eqietot[i], nGroup_, longs, longq, nt, latphi, latfs, latfq);
0647     }
0648   }
0649 }
0650 
0651 //---------------------------------------------------
0652 void HcalTestAnalysis::layerAnalysis() {
0653   int i = 0;
0654   edm::LogVerbatim("HcalSim") << "\n ===>>> HcalTestAnalysis: Energy deposit "
0655                               << "\n at EB : " << std::setw(6) << edepEB_ / MeV << "\n at EE : " << std::setw(6)
0656                               << edepEE_ / MeV << "\n at HB : " << std::setw(6) << edepHB_ / MeV
0657                               << "\n at HE : " << std::setw(6) << edepHE_ / MeV << "\n at HO : " << std::setw(6)
0658                               << edepHO_ / MeV << "\n ---- HcalTestAnalysis: Energy deposit in Layers";
0659   for (i = 0; i < 20; i++)
0660     edm::LogVerbatim("HcalSim") << " Layer " << std::setw(2) << i << " E " << std::setw(8) << edepl_[i] / MeV << " MeV";
0661 
0662   tuples_->fillLayers(edepl_, edepHO_, edepHB_ + edepHE_, mudist_);
0663 }
0664 
0665 //---------------------------------------------------
0666 double HcalTestAnalysis::timeOfFlight(int det, int layer, double eta) {
0667   double theta = 2.0 * atan(exp(-eta));
0668   double dist = 0.;
0669   if (det == static_cast<int>(HcalBarrel)) {
0670     const double rLay[19] = {1836.0,
0671                              1902.0,
0672                              1962.0,
0673                              2022.0,
0674                              2082.0,
0675                              2142.0,
0676                              2202.0,
0677                              2262.0,
0678                              2322.0,
0679                              2382.0,
0680                              2448.0,
0681                              2514.0,
0682                              2580.0,
0683                              2646.0,
0684                              2712.0,
0685                              2776.0,
0686                              2862.5,
0687                              3847.0,
0688                              4052.0};
0689     if (layer > 0 && layer < 20)
0690       dist += rLay[layer - 1] * mm / sin(theta);
0691   } else {
0692     const double zLay[19] = {4034.0,
0693                              4032.0,
0694                              4123.0,
0695                              4210.0,
0696                              4297.0,
0697                              4384.0,
0698                              4471.0,
0699                              4558.0,
0700                              4645.0,
0701                              4732.0,
0702                              4819.0,
0703                              4906.0,
0704                              4993.0,
0705                              5080.0,
0706                              5167.0,
0707                              5254.0,
0708                              5341.0,
0709                              5428.0,
0710                              5515.0};
0711     if (layer > 0 && layer < 20)
0712       dist += zLay[layer - 1] * mm / cos(theta);
0713   }
0714   double tmp = dist / c_light / ns;
0715   edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::timeOfFlight " << tmp << " for det/lay " << det << " " << layer
0716                               << " eta/theta " << eta << " " << theta / deg << " dist " << dist;
0717   return tmp;
0718 }
0719 
0720 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0721 #include "FWCore/PluginManager/interface/ModuleDef.h"
0722 
0723 DEFINE_SIMWATCHER(HcalTestAnalysis);
[ 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!