Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:24:12 

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