Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: SimG4Validation.cc
0003 // Description: Main analysis class for Hcal Validation of G4 Hits
0004 ///////////////////////////////////////////////////////////////////////////////
0005 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0006 #include "SimG4Core/Notification/interface/BeginOfJob.h"
0007 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0008 #include "SimG4Core/Notification/interface/EndOfEvent.h"
0009 #include "SimG4Core/Notification/interface/Observer.h"
0010 #include "SimG4Core/Watcher/interface/SimProducer.h"
0011 
0012 // to retreive hits
0013 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0014 #include "DataFormats/Math/interface/GeantUnits.h"
0015 #include "DataFormats/Math/interface/Point3D.h"
0016 #include "SimDataFormats/CaloHit/interface/CaloHit.h"
0017 #include "SimDataFormats/ValidationFormats/interface/PValidationFormats.h"
0018 
0019 #include "SimG4CMS/Calo/interface/CaloG4Hit.h"
0020 #include "SimG4CMS/Calo/interface/CaloG4HitCollection.h"
0021 #include "SimG4CMS/Calo/interface/HCalSD.h"
0022 #include "SimG4CMS/Calo/interface/HcalTestNumberingScheme.h"
0023 
0024 #include "FWCore/Framework/interface/Event.h"
0025 #include "FWCore/Framework/interface/EventSetup.h"
0026 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0027 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0028 
0029 #include "Geometry/HcalCommonData/interface/HcalDDDSimConstants.h"
0030 #include "Geometry/HcalCommonData/interface/HcalNumberingFromDDD.h"
0031 #include "Geometry/Records/interface/HcalSimNumberingRecord.h"
0032 
0033 #include "Validation/HcalHits/interface/SimG4HcalHitCluster.h"
0034 #include "Validation/HcalHits/interface/SimG4HcalHitJetFinder.h"
0035 
0036 #include <CLHEP/Units/SystemOfUnits.h>
0037 
0038 #include "G4HCofThisEvent.hh"
0039 #include "G4SDManager.hh"
0040 #include "G4Step.hh"
0041 #include "G4Track.hh"
0042 #include "G4VProcess.hh"
0043 
0044 #include <iostream>
0045 #include <memory>
0046 #include <string>
0047 #include <vector>
0048 
0049 using namespace geant_units::operators;
0050 using CLHEP::GeV;
0051 using CLHEP::MeV;
0052 
0053 class SimG4HcalValidation : public SimProducer,
0054                             public Observer<const BeginOfRun *>,
0055                             public Observer<const BeginOfEvent *>,
0056                             public Observer<const EndOfEvent *>,
0057                             public Observer<const G4Step *> {
0058 public:
0059   SimG4HcalValidation(const edm::ParameterSet &p);
0060   SimG4HcalValidation(const SimG4HcalValidation &) = delete;  // stop default
0061   const SimG4HcalValidation &operator=(const SimG4HcalValidation &) = delete;
0062   ~SimG4HcalValidation() override;
0063 
0064   void registerConsumes(edm::ConsumesCollector) override;
0065   void produce(edm::Event &, const edm::EventSetup &) override;
0066   void beginRun(edm::EventSetup const &) override;
0067 
0068 private:
0069   void init();
0070 
0071   // observer classes
0072   void update(const BeginOfRun *run) override;
0073   void update(const BeginOfEvent *evt) override;
0074   void update(const G4Step *step) override;
0075   void update(const EndOfEvent *evt) override;
0076 
0077   // jetfinding and analysis-related stuff
0078   void fill(const EndOfEvent *ev);
0079   void layerAnalysis(PHcalValidInfoLayer &);
0080   void nxNAnalysis(PHcalValidInfoNxN &);
0081   void jetAnalysis(PHcalValidInfoJets &);
0082   void fetchHits(PHcalValidInfoLayer &);
0083   void clear();
0084   void collectEnergyRdir(const double, const double);
0085   double getHcalScale(std::string, int) const;
0086 
0087 private:
0088   edm::ESGetToken<HcalDDDSimConstants, HcalSimNumberingRecord> ddconsToken_;
0089 
0090   // Keep parameters to instantiate Jet finder later
0091   std::unique_ptr<SimG4HcalHitJetFinder> jetf;
0092 
0093   // Keep reference to instantiate HcalNumberingFromDDD later
0094   std::unique_ptr<HcalNumberingFromDDD> numberingFromDDD;
0095 
0096   // Keep parameters to instantiate HcalTestNumberingScheme later
0097   std::unique_ptr<HcalTestNumberingScheme> org;
0098 
0099   // Hit cache for cluster analysis
0100   std::vector<CaloHit> hitcache;  // e, eta, phi, time, layer, calo type
0101 
0102   // scale factors :
0103   std::vector<float> scaleHB;
0104   std::vector<float> scaleHE;
0105   std::vector<float> scaleHF;
0106 
0107   // to read from parameter set
0108   std::vector<std::string> names;
0109   double coneSize, ehitThreshold, hhitThreshold;
0110   float timeLowlim, timeUplim, eta0, phi0, jetThreshold;
0111   bool applySampling, hcalOnly;
0112   int infolevel;
0113   std::string labelLayer, labelNxN, labelJets;
0114 
0115   // eta and phi size of windows around eta0, phi0
0116   std::vector<double> dEta;
0117   std::vector<double> dPhi;
0118 
0119   // some private members for ananlysis
0120   unsigned int count;
0121   double edepEB, edepEE, edepHB, edepHE, edepHO;
0122   double edepd[5], edepl[20];
0123   double een, hen, hoen;  // Energy sum in ECAL, HCAL, HO
0124   double vhitec, vhithc, enEcal, enHcal;
0125 };
0126 
0127 SimG4HcalValidation::SimG4HcalValidation(const edm::ParameterSet &p) {
0128   edm::ParameterSet m_Anal = p.getParameter<edm::ParameterSet>("SimG4HcalValidation");
0129   infolevel = m_Anal.getParameter<int>("InfoLevel");
0130   hcalOnly = m_Anal.getParameter<bool>("HcalClusterOnly");
0131   applySampling = m_Anal.getParameter<bool>("HcalSampling");
0132   coneSize = m_Anal.getParameter<double>("ConeSize");
0133   ehitThreshold = m_Anal.getParameter<double>("EcalHitThreshold");
0134   hhitThreshold = m_Anal.getParameter<double>("HcalHitThreshold");
0135   timeLowlim = m_Anal.getParameter<double>("TimeLowLimit");
0136   timeUplim = m_Anal.getParameter<double>("TimeUpLimit");
0137   jetThreshold = m_Anal.getParameter<double>("JetThreshold");
0138   eta0 = m_Anal.getParameter<double>("Eta0");
0139   phi0 = m_Anal.getParameter<double>("Phi0");
0140   names = m_Anal.getParameter<std::vector<std::string>>("Names");
0141   labelLayer = m_Anal.getUntrackedParameter<std::string>("LabelLayerInfo", "HcalInfoLayer");
0142   labelNxN = m_Anal.getUntrackedParameter<std::string>("LabelNxNInfo", "HcalInfoNxN");
0143   labelJets = m_Anal.getUntrackedParameter<std::string>("LabelJetsInfo", "HcalInfoJets");
0144 
0145   produces<PHcalValidInfoLayer>(labelLayer);
0146   if (infolevel > 0)
0147     produces<PHcalValidInfoNxN>(labelNxN);
0148   if (infolevel > 1)
0149     produces<PHcalValidInfoJets>(labelJets);
0150 
0151   edm::LogVerbatim("ValidHcal") << "HcalTestAnalysis:: Initialized as observer of begin/end events and "
0152                                 << "of G4step with Parameter values: \n\tInfoLevel     = " << infolevel
0153                                 << "\n\thcalOnly      = " << hcalOnly << "\n\tapplySampling = " << applySampling
0154                                 << "\n\tconeSize      = " << coneSize << "\n\tehitThreshold = " << ehitThreshold
0155                                 << "\n\thhitThreshold = " << hhitThreshold << "\n\tttimeLowlim   = " << timeLowlim
0156                                 << "\n\tttimeUplim    = " << timeUplim << "\n\teta0          = " << eta0
0157                                 << "\n\tphi0          = " << phi0 << "\nLabels (Layer): " << labelLayer
0158                                 << " (NxN): " << labelNxN << " (Jets): " << labelJets;
0159 
0160   init();
0161 }
0162 
0163 SimG4HcalValidation::~SimG4HcalValidation() {
0164   edm::LogVerbatim("ValidHcal") << "\n -------->  Total number of selected entries"
0165                                 << " : " << count;
0166 }
0167 
0168 void SimG4HcalValidation::registerConsumes(edm::ConsumesCollector cc) {
0169   ddconsToken_ = cc.esConsumes<HcalDDDSimConstants, HcalSimNumberingRecord, edm::Transition::BeginRun>();
0170   edm::LogVerbatim("ValidHcal") << "SimG4HcalValidation::Initialize ESGetToken for HcalDDDSimConstants";
0171 }
0172 
0173 void SimG4HcalValidation::produce(edm::Event &e, const edm::EventSetup &) {
0174   std::unique_ptr<PHcalValidInfoLayer> productLayer(new PHcalValidInfoLayer);
0175   layerAnalysis(*productLayer);
0176   e.put(std::move(productLayer), labelLayer);
0177 
0178   if (infolevel > 0) {
0179     std::unique_ptr<PHcalValidInfoNxN> productNxN(new PHcalValidInfoNxN);
0180     nxNAnalysis(*productNxN);
0181     e.put(std::move(productNxN), labelNxN);
0182   }
0183 
0184   if (infolevel > 1) {
0185     std::unique_ptr<PHcalValidInfoJets> productJets(new PHcalValidInfoJets);
0186     jetAnalysis(*productJets);
0187     e.put(std::move(productJets), labelJets);
0188   }
0189 }
0190 
0191 //==================================================================== per RUN
0192 void SimG4HcalValidation::init() {
0193   float sHB[4] = {117., 117., 178., 217.};
0194   float sHE[3] = {178., 178., 178.};
0195   float sHF[3] = {2.84, 2.09, 0.};  //
0196 
0197   float deta[4] = {0.0435, 0.1305, 0.2175, 0.3045};
0198   float dphi[4] = {0.0436, 0.1309, 0.2182, 0.3054};
0199 
0200   int i = 0;
0201   for (i = 0; i < 4; i++) {
0202     scaleHB.push_back(sHB[i]);
0203   }
0204   for (i = 0; i < 3; i++) {
0205     scaleHE.push_back(sHE[i]);
0206   }
0207   for (i = 0; i < 3; i++) {
0208     scaleHF.push_back(sHF[i]);
0209   }
0210 
0211   // window steps;
0212   for (i = 0; i < 4; i++) {
0213     dEta.push_back(deta[i]);
0214     dPhi.push_back(dphi[i]);
0215   }
0216 
0217   // jetfinder conse size setting
0218   jetf = std::make_unique<SimG4HcalHitJetFinder>(coneSize);
0219 
0220   // counter
0221   count = 0;
0222 }
0223 
0224 void SimG4HcalValidation::beginRun(edm::EventSetup const &es) {
0225   // Numbering From DDD
0226   const HcalDDDSimConstants *hcons = &es.getData(ddconsToken_);
0227   edm::LogVerbatim("ValidHcal") << "HcalTestAnalysis:: Initialise HcalNumberingFromDDD";
0228   numberingFromDDD = std::make_unique<HcalNumberingFromDDD>(hcons);
0229 
0230   // Numbering scheme
0231   org = std::make_unique<HcalTestNumberingScheme>(false);
0232 }
0233 
0234 void SimG4HcalValidation::update(const BeginOfRun *run) {
0235   int irun = (*run)()->GetRunID();
0236 
0237   edm::LogVerbatim("ValidHcal") << " =====> Begin of Run = " << irun;
0238 
0239   std::string sdname = names[0];
0240   G4SDManager *sd = G4SDManager::GetSDMpointerIfExist();
0241   if (sd != nullptr) {
0242     G4VSensitiveDetector *aSD = sd->FindSensitiveDetector(sdname);
0243     if (aSD == nullptr) {
0244       edm::LogWarning("ValidHcal") << "SimG4HcalValidation::beginOfRun: No SD"
0245                                    << " with name " << sdname << " in this "
0246                                    << "Setup";
0247     } else {
0248       HCalSD *theCaloSD = dynamic_cast<HCalSD *>(aSD);
0249       edm::LogVerbatim("ValidHcal") << "SimG4HcalValidation::beginOfRun: Finds SD with name " << theCaloSD->GetName()
0250                                     << " in this Setup";
0251       if (org.get()) {
0252         theCaloSD->setNumberingScheme(org.get());
0253         edm::LogVerbatim("ValidHcal") << "SimG4HcalValidation::beginOfRun: set a new numbering scheme";
0254       }
0255     }
0256   } else {
0257     edm::LogWarning("ValidHcal") << "SimG4HcalValidation::beginOfRun: Could "
0258                                  << "not get SD Manager!";
0259   }
0260 }
0261 
0262 //=================================================================== per EVENT
0263 void SimG4HcalValidation::update(const BeginOfEvent *evt) {
0264   int i = 0;
0265   edepEB = edepEE = edepHB = edepHE = edepHO = 0.;
0266   for (i = 0; i < 5; i++)
0267     edepd[i] = 0.;
0268   for (i = 0; i < 20; i++)
0269     edepl[i] = 0.;
0270   vhitec = vhithc = enEcal = enHcal = 0;
0271   // Cache reset
0272   clear();
0273 
0274   int iev = (*evt)()->GetEventID();
0275   LogDebug("ValidHcal") << "SimG4HcalValidation: =====> Begin of event = " << iev;
0276 }
0277 
0278 //=================================================================== each STEP
0279 void SimG4HcalValidation::update(const G4Step *aStep) {
0280   if (aStep != nullptr) {
0281     G4VPhysicalVolume *curPV = aStep->GetPreStepPoint()->GetPhysicalVolume();
0282     G4String name = curPV->GetName();
0283     name.assign(name, 0, 3);
0284     double edeposit = aStep->GetTotalEnergyDeposit();
0285     int layer = -1, depth = -1;
0286     if (name == "EBR") {
0287       depth = 0;
0288       edepEB += edeposit;
0289     } else if (name == "EFR") {
0290       depth = 0;
0291       edepEE += edeposit;
0292     } else if (name == "HBS") {
0293       layer = (curPV->GetCopyNo() / 10) % 100;
0294       depth = (curPV->GetCopyNo()) % 10 + 1;
0295       if (depth > 0 && depth < 4 && layer >= 0 && layer < 17) {
0296         edepHB += edeposit;
0297       } else {
0298         edm::LogWarning("ValidHcal") << "SimG4HcalValidation:Error " << curPV->GetName() << curPV->GetCopyNo();
0299         depth = -1;
0300         layer = -1;
0301       }
0302     } else if (name == "HES") {
0303       layer = (curPV->GetCopyNo() / 10) % 100;
0304       depth = (curPV->GetCopyNo()) % 10 + 1;
0305       if (depth > 0 && depth < 3 && layer >= 0 && layer < 19) {
0306         edepHE += edeposit;
0307       } else {
0308         edm::LogWarning("ValidHcal") << "SimG4HcalValidation:Error " << curPV->GetName() << curPV->GetCopyNo();
0309         depth = -1;
0310         layer = -1;
0311       }
0312     } else if (name == "HTS") {
0313       layer = (curPV->GetCopyNo() / 10) % 100;
0314       depth = (curPV->GetCopyNo()) % 10 + 1;
0315       if (depth > 3 && depth < 5 && layer >= 17 && layer < 20) {
0316         edepHO += edeposit;
0317       } else {
0318         edm::LogWarning("ValidHcal") << "SimG4HcalValidation:Error " << curPV->GetName() << curPV->GetCopyNo();
0319         depth = -1;
0320         layer = -1;
0321       }
0322     }
0323     if (depth >= 0 && depth < 5)
0324       edepd[depth] += edeposit;
0325     if (layer >= 0 && layer < 20)
0326       edepl[layer] += edeposit;
0327 
0328     if (layer >= 0 && layer < 20) {
0329       LogDebug("ValidHcal") << "SimG4HcalValidation:: G4Step: " << name << " Layer " << std::setw(3) << layer
0330                             << " Depth " << std::setw(2) << depth << " Edep " << std::setw(6) << edeposit / MeV
0331                             << " MeV";
0332     }
0333   }
0334 }
0335 
0336 //================================================================ End of EVENT
0337 void SimG4HcalValidation::update(const EndOfEvent *evt) {
0338   count++;
0339 
0340   // Fill hits cache for jetfinding etc.
0341   fill(evt);
0342   LogDebug("ValidHcal") << "SimG4HcalValidation:: ---  after Fill ";
0343 }
0344 
0345 //---------------------------------------------------
0346 void SimG4HcalValidation::fill(const EndOfEvent *evt) {
0347   LogDebug("ValidHcal") << "SimG4HcalValidation:Fill event " << (*evt)()->GetEventID();
0348 
0349   // access to the G4 hit collections
0350   G4HCofThisEvent *allHC = (*evt)()->GetHCofThisEvent();
0351 
0352   int nhc = 0, j = 0;
0353 
0354   // Hcal
0355   int HCHCid = G4SDManager::GetSDMpointer()->GetCollectionID(names[0]);
0356   CaloG4HitCollection *theHCHC = (CaloG4HitCollection *)allHC->GetHC(HCHCid);
0357   LogDebug("ValidHcal") << "SimG4HcalValidation :: Hit Collection for " << names[0] << " of ID " << HCHCid
0358                         << " is obtained at " << theHCHC;
0359   if (HCHCid >= 0 && theHCHC != nullptr) {
0360     int nhits = theHCHC->entries();
0361     for (j = 0; j < nhits; j++) {
0362       CaloG4Hit *aHit = (*theHCHC)[j];
0363 
0364       double e = aHit->getEnergyDeposit() / GeV;
0365       double time = aHit->getTimeSlice();
0366 
0367       math::XYZPoint pos = aHit->getPosition();
0368       double theta = pos.theta();
0369       double eta = -log(tan(theta * 0.5));
0370       double phi = pos.phi();
0371 
0372       uint32_t unitID = aHit->getUnitID();
0373       int subdet, zside, layer, etaIndex, phiIndex, lay;
0374       org->unpackHcalIndex(unitID, subdet, zside, layer, etaIndex, phiIndex, lay);
0375 
0376       // some logic to separate HO ...
0377       layer--;
0378       std::string det = "HB";
0379       if (subdet == static_cast<int>(HcalForward)) {
0380         det = "HF";
0381         uint16_t depth = aHit->getDepth();
0382         if (depth != 0) {
0383           layer += 2;
0384           lay += 2;
0385         }
0386         if (layer == 1)
0387           vhithc += e;
0388         else if (layer == 0)
0389           vhitec += e;
0390         else
0391           edm::LogVerbatim("ValidHcal") << "SimG4HcalValidation::HitPMT " << subdet << " " << (2 * zside - 1) * etaIndex
0392                                         << " " << phiIndex << " " << layer + 1 << " R " << pos.rho() << " Phi "
0393                                         << convertRadToDeg(phi) << " Edep " << e << " Time " << time;
0394       } else if (subdet == static_cast<int>(HcalEndcap)) {
0395         if (etaIndex <= 20) {
0396           det = "HES";
0397         } else {
0398           det = "HED";
0399         }
0400       }
0401       LogDebug("ValidHcal") << "SimG4HcalValidation::debugFill Hcal " << det << " layer " << std::setw(2) << layer
0402                             << " lay " << std::setw(2) << lay << " time " << std::setw(6) << time << " theta "
0403                             << std::setw(8) << theta << " eta " << std::setw(8) << eta << " phi " << std::setw(8) << phi
0404                             << " e " << std::setw(8) << e << " ID 0x" << std::hex << unitID << " ID dec " << std::dec
0405                             << (int)unitID;
0406 
0407       // if desired, apply sampling factors in HCAL !!!
0408       if (applySampling)
0409         e *= getHcalScale(det, layer);
0410 
0411       //    filter on time & energy
0412       if (time >= timeLowlim && time <= timeUplim && e > hhitThreshold) {
0413         enHcal += e;
0414         CaloHit ahit(subdet, lay, e, eta, phi, time, unitID);
0415         hitcache.push_back(ahit);  // fill cache
0416         ++nhc;
0417       }
0418     }
0419   }
0420   LogDebug("ValidHcal") << "SimG4HcalValidation:: HCAL hits : " << nhc;
0421 
0422   if (!hcalOnly) {  //--------------------------  ECAL hits --------------------
0423     int ndets = names.size();
0424     if (ndets > 3)
0425       ndets = 3;
0426     for (int idty = 1; idty < ndets; idty++) {
0427       std::string det = "EB";
0428       if (idty == 2)
0429         det = "EE";
0430       else if (idty > 2)
0431         det = "ES";
0432 
0433       int nec = 0;
0434       int ECHCid = G4SDManager::GetSDMpointer()->GetCollectionID(names[idty]);
0435       CaloG4HitCollection *theECHC = (CaloG4HitCollection *)allHC->GetHC(ECHCid);
0436       LogDebug("ValidHcal") << "SimG4HcalValidation:: Hit Collection for " << names[idty] << " of ID " << ECHCid
0437                             << " is obtained at " << theECHC;
0438       int theechc_entries = theECHC->entries();
0439       if (ECHCid >= 0 && theECHC != nullptr) {
0440         for (j = 0; j < theechc_entries; j++) {
0441           CaloG4Hit *aHit = (*theECHC)[j];
0442 
0443           double e = aHit->getEnergyDeposit() / GeV;
0444           double time = aHit->getTimeSlice();
0445 
0446           math::XYZPoint pos = aHit->getPosition();
0447           double theta = pos.theta();
0448           double eta = -log(tan(theta / 2.));
0449           double phi = pos.phi();
0450           HcalNumberingFromDDD::HcalID id = numberingFromDDD->unitID(eta, phi, 1, 1);
0451           uint32_t unitID = org->getUnitID(id);
0452           int subdet, zside, layer, ieta, iphi, lay;
0453           org->unpackHcalIndex(unitID, subdet, zside, layer, ieta, iphi, lay);
0454           subdet = idty + 9;
0455           layer = 0;
0456           unitID = org->packHcalIndex(subdet, zside, layer, ieta, iphi, lay);
0457 
0458           //    filter on time & energy
0459           if (time >= timeLowlim && time <= timeUplim && e > ehitThreshold) {
0460             enEcal += e;
0461             CaloHit ahit(subdet, lay, e, eta, phi, time, unitID);
0462             hitcache.push_back(ahit);  // fill cache
0463             ++nec;
0464           }
0465 
0466           LogDebug("ValidHcal") << "SimG4HcalValidation::debugFill Ecal " << det << " layer " << std::setw(2) << layer
0467                                 << " lay " << std::setw(2) << lay << " time " << std::setw(6) << time << " theta "
0468                                 << std::setw(8) << theta << " eta " << std::setw(8) << eta << " phi " << std::setw(8)
0469                                 << phi << " e " << std::setw(8) << e << " ID 0x" << std::hex << unitID << " ID dec "
0470                                 << std::dec << (int)unitID;
0471         }
0472       }
0473 
0474       LogDebug("ValidHcal") << "SimG4HcalValidation:: " << det << " hits : " << nec;
0475     }
0476   }  // end of if(!hcalOnly)
0477 }
0478 
0479 void SimG4HcalValidation::layerAnalysis(PHcalValidInfoLayer &product) {
0480   int i = 0;
0481   LogDebug("ValidHcal") << "\n ===>>> SimG4HcalValidation: Energy deposit "
0482                         << "in MeV\n at EB : " << std::setw(6) << edepEB / MeV << "\n at EE : " << std::setw(6)
0483                         << edepEE / MeV << "\n at HB : " << std::setw(6) << edepHB / MeV
0484                         << "\n at HE : " << std::setw(6) << edepHE / MeV << "\n at HO : " << std::setw(6)
0485                         << edepHO / MeV << "\n ---- SimG4HcalValidation: Energy deposit in";
0486   for (i = 0; i < 5; i++)
0487     LogDebug("ValidHcal") << " Depth " << std::setw(2) << i << " E " << std::setw(8) << edepd[i] / MeV << " MeV";
0488   LogDebug("ValidHcal") << " ---- SimG4HcalValidation: Energy deposit in"
0489                         << "layers";
0490   for (i = 0; i < 20; i++)
0491     LogDebug("ValidHcal") << " Layer " << std::setw(2) << i << " E " << std::setw(8) << edepl[i] / MeV << " MeV";
0492 
0493   product.fillLayers(edepl, edepd, edepHO, edepHB + edepHE, edepEB + edepEE);
0494 
0495   // Hits in HF
0496   product.fillHF(vhitec, vhithc, enEcal, enHcal);
0497   LogDebug("ValidHcal") << "SimG4HcalValidation::HF hits " << vhitec << " in EC and " << vhithc << " in HC\n"
0498                         << "                     HB/HE   " << enEcal << " in EC and " << enHcal << " in HC";
0499 
0500   // Another HCAL hist to porcess and store separately (a bit more complicated)
0501   fetchHits(product);
0502 
0503   LogDebug("ValidHcal") << " layerAnalysis ===> after fetchHits";
0504 }
0505 
0506 //-----------------------------------------------------------------------------
0507 void SimG4HcalValidation::nxNAnalysis(PHcalValidInfoNxN &product) {
0508   std::vector<CaloHit> *hits = &hitcache;
0509   std::vector<CaloHit>::iterator hit_itr;
0510 
0511   LogDebug("ValidHcal") << "SimG4HcalValidation::NxNAnalysis : entrance ";
0512 
0513   collectEnergyRdir(eta0, phi0);  // HCAL and ECAL energy in  SimHitCache
0514                                   // around (eta0,phi0)
0515 
0516   LogDebug("ValidHcal") << " NxNAnalysis : coolectEnergyRdir - Ecal " << een << "   Hcal " << hen;
0517 
0518   double etot = 0.;  // total e deposited     in "cluster"
0519   double ee = 0.;    // ECAL  e deposited     in "cluster"
0520   double he = 0.;    // HCAL  e deposited     in "cluster"
0521   double hoe = 0.;   // HO    e deposited     in "cluster"
0522 
0523   int max = dEta.size();  // 4
0524 
0525   for (hit_itr = hits->begin(); hit_itr < hits->end(); hit_itr++) {
0526     double e = hit_itr->e();
0527     double t = hit_itr->t();
0528     double eta = hit_itr->eta();
0529     double phi = hit_itr->phi();
0530     int type = hit_itr->det();
0531     int layer = hit_itr->layer();
0532 
0533     // NxN calulation
0534 
0535     if (fabs(eta0 - eta) <= dEta[max - 1] && fabs(phi0 - phi) <= dPhi[max - 1]) {
0536       etot += e;
0537       if (type == 10 || type == 11 || type == 12)
0538         ee += e;
0539       if (type == static_cast<int>(HcalBarrel) || type == static_cast<int>(HcalEndcap) ||
0540           type == static_cast<int>(HcalForward)) {
0541         he += e;
0542         if (type == static_cast<int>(HcalBarrel) && layer > 17)
0543           hoe += e;
0544 
0545         // which concrete i-th square ?
0546         for (int i = 0; i < max; i++) {
0547           if ((eta0 - eta) <= dEta[i] && (eta0 - eta) >= -dEta[i] && (phi0 - phi) <= dPhi[i] &&
0548               (phi0 - phi) >= -dPhi[i]) {
0549             LogDebug("ValidHcal") << "SimG4HcalValidation:: nxNAnalysis eta0,"
0550                                   << " phi0 = " << eta0 << " " << phi0 << "   type, layer = " << type << "," << layer
0551                                   << "  eta, phi = " << eta << " " << phi;
0552 
0553             product.fillTProfileNxN(e, i, t);
0554             break;
0555           }
0556         }
0557         // here comes break ...
0558       }
0559     }
0560   }
0561 
0562   product.fillEcollectNxN(ee, he, hoe, etot);
0563   product.fillHvsE(een, hen, hoen, een + hen);
0564 
0565   LogDebug("ValidHcal") << " nxNAnalysis ===> after fillHvsE";
0566 }
0567 
0568 //-----------------------------------------------------------------------------
0569 void SimG4HcalValidation::jetAnalysis(PHcalValidInfoJets &product) {
0570   std::vector<CaloHit> *hhit = &hitcache;
0571 
0572   jetf->setInput(hhit);
0573 
0574   // zeroing cluster list, perfom clustering, fill cluster list & return pntr
0575   std::vector<SimG4HcalHitCluster> *result = jetf->getClusters(hcalOnly);
0576 
0577   std::vector<SimG4HcalHitCluster>::iterator clus_itr;
0578 
0579   LogDebug("ValidHcal") << "\n ---------- Final list of " << (*result).size() << " clusters ---------------";
0580   for (clus_itr = result->begin(); clus_itr < result->end(); clus_itr++)
0581     LogDebug("ValidHcal") << (*clus_itr);
0582 
0583   std::vector<double> enevec, etavec, phivec;
0584 
0585   if (!(*result).empty()) {
0586     sort((*result).begin(), (*result).end());
0587 
0588     clus_itr = result->begin();  // first cluster only
0589     double etac = clus_itr->eta();
0590     double phic = clus_itr->phi();
0591 
0592     double ecal_collect = 0.;  // collect Ecal energy in the cone
0593     if (!hcalOnly) {
0594       ecal_collect = clus_itr->collectEcalEnergyR();
0595     } else {
0596       collectEnergyRdir(etac, phic);
0597       ecal_collect = een;
0598     }
0599     LogDebug("ValidHcal") << " JetAnalysis ===> ecal_collect  = " << ecal_collect;
0600 
0601     // eta-phi deviation of the cluster from nominal (eta0, phi0) values
0602     double dist = jetf->rDist(eta0, phi0, etac, phic);
0603     LogDebug("ValidHcal") << " JetAnalysis ===> eta phi deviation  = " << dist;
0604     product.fillEtaPhiProfileJet(eta0, phi0, etac, phic, dist);
0605 
0606     LogDebug("ValidHcal") << " JetAnalysis ===> after fillEtaPhiProfileJet";
0607 
0608     std::vector<CaloHit> *hits = clus_itr->getHits();
0609     std::vector<CaloHit>::iterator hit_itr;
0610 
0611     double ee = 0., he = 0., hoe = 0., etot = 0.;
0612 
0613     // cycle over all hits in the FIRST cluster
0614     for (hit_itr = hits->begin(); hit_itr < hits->end(); hit_itr++) {
0615       double e = hit_itr->e();
0616       double t = hit_itr->t();
0617       double r = jetf->rDist(&(*clus_itr), &(*hit_itr));
0618 
0619       // energy collection
0620       etot += e;
0621       if (hit_itr->det() == 10 || hit_itr->det() == 11 || hit_itr->det() == 12)
0622         ee += e;
0623       if (hit_itr->det() == static_cast<int>(HcalBarrel) || hit_itr->det() == static_cast<int>(HcalEndcap) ||
0624           hit_itr->det() == static_cast<int>(HcalForward)) {
0625         he += e;
0626         if (hit_itr->det() == static_cast<int>(HcalBarrel) && hit_itr->layer() > 17)
0627           hoe += e;
0628       }
0629 
0630       if (hit_itr->det() == static_cast<int>(HcalBarrel) || hit_itr->det() == static_cast<int>(HcalEndcap) ||
0631           hit_itr->det() == static_cast<int>(HcalForward)) {
0632         product.fillTProfileJet(he, r, t);
0633       }
0634     }
0635 
0636     product.fillEcollectJet(ee, he, hoe, etot);
0637 
0638     LogDebug("ValidHcal") << " JetAnalysis ===> after fillEcollectJet: "
0639                           << "ee/he/hoe/etot " << ee << "/" << he << "/" << hoe << "/" << etot;
0640 
0641     // Loop over clusters
0642     for (clus_itr = result->begin(); clus_itr < result->end(); clus_itr++) {
0643       if (clus_itr->e() > jetThreshold) {
0644         enevec.push_back(clus_itr->e());
0645         etavec.push_back(clus_itr->eta());
0646         phivec.push_back(clus_itr->phi());
0647       }
0648     }
0649     product.fillJets(enevec, etavec, phivec);
0650 
0651     LogDebug("ValidHcal") << " JetAnalysis ===> after fillJets\n"
0652                           << " JetAnalysis ===> (*result).size() " << (*result).size();
0653 
0654     // Di-jet mass
0655     if (etavec.size() > 1) {
0656       if (etavec[0] > -2.5 && etavec[0] < 2.5 && etavec[1] > -2.5 && etavec[1] < 2.5) {
0657         LogDebug("ValidHcal") << " JetAnalysis ===> Di-jet mass enter\n"
0658                               << " JetAnalysis ===> Di-jet vectors ";
0659         for (unsigned int i = 0; i < enevec.size(); i++)
0660           LogDebug("ValidHcal") << "   e, eta, phi = " << enevec[i] << " " << etavec[i] << " " << phivec[i];
0661 
0662         double et0 = enevec[0] / cosh(etavec[0]);
0663         double px0 = et0 * cos(phivec[0]);
0664         double py0 = et0 * sin(phivec[0]);
0665         double pz0 = et0 * sinh(etavec[0]);
0666         double et1 = enevec[1] / cosh(etavec[1]);
0667         double px1 = et1 * cos(phivec[1]);
0668         double py1 = et1 * sin(phivec[1]);
0669         double pz1 = et1 * sinh(etavec[1]);
0670 
0671         double dijetmass2 = (enevec[0] + enevec[1]) * (enevec[0] + enevec[1]) - (px1 + px0) * (px1 + px0) -
0672                             (py1 + py0) * (py1 + py0) - (pz1 + pz0) * (pz1 + pz0);
0673 
0674         LogDebug("ValidHcal") << " JetAnalysis ===> Di-jet massSQ " << dijetmass2;
0675 
0676         double dijetmass;
0677         if (dijetmass2 >= 0.)
0678           dijetmass = sqrt(dijetmass2);
0679         else
0680           dijetmass = -sqrt(-1. * dijetmass2);
0681 
0682         product.fillDiJets(dijetmass);
0683 
0684         LogDebug("ValidHcal") << " JetAnalysis ===> after fillDiJets";
0685       }
0686     }
0687   }
0688 }
0689 
0690 //---------------------------------------------------
0691 void SimG4HcalValidation::fetchHits(PHcalValidInfoLayer &product) {
0692   LogDebug("ValidHcal") << "Enter SimG4HcalValidation::fetchHits with " << hitcache.size() << " hits";
0693   int nHit = hitcache.size();
0694   int hit = 0;
0695   int i;
0696   std::vector<CaloHit>::iterator itr;
0697   std::vector<CaloHit *> lhits(nHit);
0698   for (i = 0, itr = hitcache.begin(); itr != hitcache.end(); i++, itr++) {
0699     uint32_t unitID = itr->id();
0700     int subdet, zside, group, ieta, iphi, lay;
0701     HcalTestNumbering::unpackHcalIndex(unitID, subdet, zside, group, ieta, iphi, lay);
0702     subdet = itr->det();
0703     lay = itr->layer();
0704     group = (subdet & 15) << 20;
0705     group += ((lay - 1) & 31) << 15;
0706     group += (zside & 1) << 14;
0707     group += (ieta & 127) << 7;
0708     group += (iphi & 127);
0709     itr->setId(group);
0710     lhits[i] = &hitcache[i];
0711     LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits:Original " << i << " " << hitcache[i] << "\n"
0712                           << "SimG4HcalValidation::fetchHits:Copied   " << i << " " << *lhits[i];
0713   }
0714   sort(lhits.begin(), lhits.end(), CaloHitIdMore());
0715   std::vector<CaloHit *>::iterator k1, k2;
0716   for (i = 0, k1 = lhits.begin(); k1 != lhits.end(); i++, k1++)
0717     LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits:Sorted " << i << " " << **k1;
0718   int nHits = 0;
0719   for (i = 0, k1 = lhits.begin(); k1 != lhits.end(); i++, k1++) {
0720     double ehit = (**k1).e();
0721     double t = (**k1).t();
0722     uint32_t unitID = (**k1).id();
0723     int jump = 0;
0724     LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits:Start " << i << " U/T/E"
0725                           << " 0x" << std::hex << unitID << std::dec << " " << t << " " << ehit;
0726     for (k2 = k1 + 1; k2 != lhits.end() && (t - (**k2).t()) < 1 && (t - (**k2).t()) > -1 && unitID == (**k2).id();
0727          k2++) {
0728       ehit += (**k2).e();
0729       LogDebug("ValidHcal") << "\t + " << (**k2).e();
0730       jump++;
0731     }
0732     LogDebug("ValidHcal") << "\t = " << ehit << " in " << jump;
0733 
0734     double eta = (*k1)->eta();
0735     double phi = (*k1)->phi();
0736     int lay = ((unitID >> 15) & 31) + 1;
0737     int subdet = (unitID >> 20) & 15;
0738     int zside = (unitID >> 14) & 1;
0739     int ieta = (unitID >> 7) & 127;
0740     int iphi = (unitID) & 127;
0741 
0742     // All hits in cache
0743     product.fillHits(nHits, lay, subdet, eta, phi, ehit, t);
0744     nHits++;
0745 
0746     LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits:Hit " << nHits << " " << i << " ID 0x" << std::hex
0747                           << unitID << "   det " << std::dec << subdet << " " << lay << " " << zside << " " << ieta
0748                           << " " << iphi << " Time " << t << " E " << ehit;
0749 
0750     i += jump;
0751     k1 += jump;
0752   }
0753 
0754   LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits called with " << nHit << " hits"
0755                         << " and writes out " << nHits << '(' << hit << ") hits";
0756 }
0757 //---------------------------------------------------
0758 void SimG4HcalValidation::clear() { hitcache.erase(hitcache.begin(), hitcache.end()); }
0759 
0760 //---------------------------------------------------
0761 void SimG4HcalValidation::collectEnergyRdir(const double eta0, const double phi0) {
0762   std::vector<CaloHit> *hits = &hitcache;
0763   std::vector<CaloHit>::iterator hit_itr;
0764 
0765   double sume = 0., sumh = 0., sumho = 0.;
0766 
0767   for (hit_itr = hits->begin(); hit_itr < hits->end(); hit_itr++) {
0768     double e = hit_itr->e();
0769     double eta = hit_itr->eta();
0770     double phi = hit_itr->phi();
0771 
0772     int type = hit_itr->det();
0773 
0774     double r = jetf->rDist(eta0, phi0, eta, phi);
0775     if (r < coneSize) {
0776       if (type == 10 || type == 11 || type == 12) {
0777         sume += e;
0778       } else {
0779         sumh += e;
0780         if (type == static_cast<int>(HcalBarrel) && hit_itr->layer() > 17)
0781           sumho += e;
0782       }
0783     }
0784   }
0785 
0786   een = sume;
0787   hen = sumh;
0788   hoen = sumho;
0789 }
0790 
0791 //---------------------------------------------------
0792 double SimG4HcalValidation::getHcalScale(std::string det, int layer) const {
0793   double tmp = 0.;
0794 
0795   if (det == "HB") {
0796     tmp = scaleHB[layer];
0797   } else if (det == "HES" || det == "HED") {
0798     tmp = scaleHE[layer];
0799   } else if (det == "HF") {
0800     tmp = scaleHF[layer];
0801   }
0802 
0803   return tmp;
0804 }
0805 
0806 #include "FWCore/Framework/interface/MakerMacros.h"
0807 #include "FWCore/PluginManager/interface/ModuleDef.h"
0808 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0809 
0810 DEFINE_SIMWATCHER(SimG4HcalValidation);
[ 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!