HcalTestBeam/plugins/HcalTB02Analysis.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:     HcalTestBeam
0004 // Class  :     HcalTB02Analysis
0005 //
0006 // Implementation:
0007 //     Main analysis class for Hcal Test Beam 2002 Analysis
0008 //
0009 // Original Author:
0010 //         Created:  Sun May 21 10:14:34 CEST 2006
0011 //
0012
0013 // system include files
0014 #include <cmath>
0015 #include <iostream>
0016 #include <iomanip>
0017 #include <memory>
0018 #include <map>
0019 #include <vector>
0020 #include <string>
0021
0022 // user include files
0023 // to retreive hits
0024 #include "SimDataFormats/HcalTestBeam/interface/HcalTB02HistoClass.h"
0025 #include "SimG4CMS/Calo/interface/CaloG4Hit.h"
0026 #include "SimG4CMS/Calo/interface/CaloG4HitCollection.h"
0027 #include "SimG4Core/Notification/interface/Observer.h"
0028 #include "SimG4Core/Notification/interface/BeginOfJob.h"
0029 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0030 #include "SimG4Core/Notification/interface/EndOfEvent.h"
0031 #include "SimG4Core/Watcher/interface/SimProducer.h"
0032 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0033
0034 #include "FWCore/Framework/interface/Event.h"
0035 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0036 #include "FWCore/Framework/interface/MakerMacros.h"
0037 #include "FWCore/PluginManager/interface/ModuleDef.h"
0038
0039 #include "SimG4CMS/HcalTestBeam/interface/HcalTB02HcalNumberingScheme.h"
0040 #include "SimG4CMS/HcalTestBeam/interface/HcalTB02XtalNumberingScheme.h"
0041 #include "SimG4CMS/HcalTestBeam/interface/HcalTB02Histo.h"
0042
0043 #include "G4HCofThisEvent.hh"
0044 #include "G4SDManager.hh"
0045 #include "G4Step.hh"
0046 #include "G4Track.hh"
0047 #include "G4ThreeVector.hh"
0048 #include "G4VProcess.hh"
0049
0050 #include <CLHEP/Random/RandGaussQ.h>
0051 #include <CLHEP/Random/Randomize.h>
0052 #include <CLHEP/Units/GlobalSystemOfUnits.h>
0053 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0054
0055 //#define EDM_ML_DEBUG
0056
0057 namespace CLHEP {
0058   class HepRandomEngine;
0059 }
0060
0061 class HcalTB02Analysis : public SimProducer, public Observer<const BeginOfEvent*>, public Observer<const EndOfEvent*> {
0062 public:
0063   HcalTB02Analysis(const edm::ParameterSet& p);
0064   HcalTB02Analysis(const HcalTB02Analysis&) = delete;  // stop default
0065   const HcalTB02Analysis& operator=(const HcalTB02Analysis&) = delete;
0066   ~HcalTB02Analysis() override;
0067
0068   void produce(edm::Event&, const edm::EventSetup&) override;
0069
0070 private:
0071   // observer methods
0072   void update(const BeginOfEvent* evt) override;
0073   void update(const EndOfEvent* evt) override;
0074
0075   void fillEvent(HcalTB02HistoClass&);
0076   void clear();
0077   void finish();
0078
0079 private:
0080   // Private Tuples
0081   std::unique_ptr<HcalTB02Histo> histo;
0082
0083   // to read from parameter set
0084   const edm::ParameterSet m_Anal;
0085   const bool hcalOnly;
0086   const std::vector<std::string> names;
0087
0088   //To be saved
0089   std::map<int, float> energyInScints, energyInCrystals;
0090   std::map<int, float> primaries;
0091   int particleType;
0092   double eta, phi, pInit, incidentEnergy;
0093   float SEnergy, E7x7Matrix, E5x5Matrix;
0094   float SEnergyN, E7x7MatrixN, E5x5MatrixN;
0095   int maxTime;
0096   double xIncidentEnergy;
0097   float xSEnergy, xSEnergyN;
0098   float xE3x3Matrix, xE5x5Matrix;
0099   float xE3x3MatrixN, xE5x5MatrixN;
0100 };
0101
0102 //
0103 // constructors and destructor
0104 //
0105
0106 HcalTB02Analysis::HcalTB02Analysis(const edm::ParameterSet& p)
0107     : m_Anal(p.getParameter<edm::ParameterSet>("HcalTB02Analysis")),
0108       hcalOnly(m_Anal.getUntrackedParameter<bool>("HcalClusterOnly", true)),
0109       names(m_Anal.getParameter<std::vector<std::string> >("Names")) {
0110   produces<HcalTB02HistoClass>();
0111
0112   edm::LogVerbatim("HcalTBSim") << "HcalTB02Analysis:: Initialised as observer of "
0113                                 << "BeginOfJob/BeginOfEvent/EndOfEvent with "
0114                                 << "Parameter values:\n \thcalOnly = " << hcalOnly;
0115
0116   histo = std::make_unique<HcalTB02Histo>(m_Anal);
0117 }
0118
0119 HcalTB02Analysis::~HcalTB02Analysis() { finish(); }
0120
0121 //
0122 // member functions
0123 //
0124
0125 void HcalTB02Analysis::produce(edm::Event& e, const edm::EventSetup&) {
0126   std::unique_ptr<HcalTB02HistoClass> product(new HcalTB02HistoClass);
0127   fillEvent(*product);
0128   e.put(std::move(product));
0129 }
0130
0131 void HcalTB02Analysis::update(const BeginOfEvent* evt) {
0132 #ifdef EDM_ML_DEBUG
0133   edm::LogVerbatim("HcalTBSim") << "HcalTB02Analysis: =====> Begin of event = " << (*evt)()->GetEventID();
0134 #endif
0135   clear();
0136 }
0137
0138 void HcalTB02Analysis::update(const EndOfEvent* evt) {
0139   CLHEP::HepRandomEngine* engine = G4Random::getTheEngine();
0140   CLHEP::RandGaussQ randGauss(*engine);
0141
0142   // Look for the Hit Collection
0143 #ifdef EDM_ML_DEBUG
0144   edm::LogVerbatim("HcalTBSim") << "HcalTB02Analysis::Fill event " << (*evt)()->GetEventID();
0145 #endif
0146   // access to the G4 hit collections
0147   G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
0148   int ihit = 0;
0149
0150   // HCAL
0151   std::string sd = names[0];
0152   int HCHCid = G4SDManager::GetSDMpointer()->GetCollectionID(sd);
0153   CaloG4HitCollection* theHCHC = (CaloG4HitCollection*)allHC->GetHC(HCHCid);
0154 #ifdef EDM_ML_DEBUG
0155   edm::LogVerbatim("HcalTBSim") << "HcalTB02Analysis :: Hit Collection for " << sd << " of ID " << HCHCid
0156                                 << " is obtained at " << theHCHC;
0157 #endif
0158   int nentries = 0;
0159   nentries = theHCHC->entries();
0160   if (nentries == 0)
0161     return;
0162
0163   int xentries = 0;
0164   int XTALSid = 0;
0165   CaloG4HitCollection* theXTHC = nullptr;
0166
0167   if (!hcalOnly) {
0168     // XTALS
0169     sd = names[1];
0170     XTALSid = G4SDManager::GetSDMpointer()->GetCollectionID(sd);
0171     theXTHC = (CaloG4HitCollection*)allHC->GetHC(XTALSid);
0172 #ifdef EDM_ML_DEBUG
0173     edm::LogVerbatim("HcalTBSim") << "HcalTB02Analysis :: Hit Collection for " << sd << " of ID " << XTALSid
0174                                   << " is obtained at " << theXTHC;
0175 #endif
0176     xentries = theXTHC->entries();
0177     if (xentries == 0)
0178       return;
0179   }
0180
0181 #ifdef EDM_ML_DEBUG
0182   edm::LogVerbatim("HcalTBSim") << "HcalTB02Analysis :: There are " << nentries << " HCal hits, and" << xentries
0183                                 << " xtal hits";
0184 #endif
0185   float ETot = 0., xETot = 0.;
0186   int scintID = 0, xtalID = 0;
0187
0188   // HCAL
0189   std::unique_ptr<HcalTB02HcalNumberingScheme> org(new HcalTB02HcalNumberingScheme());
0190
0191   if (HCHCid >= 0 && theHCHC != nullptr) {
0192     for (ihit = 0; ihit < nentries; ihit++) {
0193       CaloG4Hit* aHit = (*theHCHC)[ihit];
0194       scintID = aHit->getUnitID();
0195       int layer = org->getlayerID(scintID);
0196       float enEm = aHit->getEM();
0197       float enhad = aHit->getHadr();
0198
0199       if (layer == 0) {
0200         enEm = enEm / 2.;
0201         enhad = enhad / 2.;
0202       }
0203
0204       energyInScints[scintID] += enEm + enhad;
0205       primaries[aHit->getTrackID()] += enEm + enhad;
0206       float time = aHit->getTimeSliceID();
0207       if (time > maxTime)
0208         maxTime = static_cast<int>(time);
0209       histo->fillAllTime(time);
0210     }
0211
0212     //
0213     // Profile
0214     //
0215
0216     float TowerEne[8][18], TowerEneCF[8][18], LayerEne[19], EnRing[100];
0217     for (int iphi = 0; iphi < 8; iphi++) {
0218       for (int jeta = 0; jeta < 18; jeta++) {
0219         TowerEne[iphi][jeta] = 0.;
0220         TowerEneCF[iphi][jeta] = 0.;
0221       }
0222     }
0223
0224     for (int ilayer = 0; ilayer < 19; ilayer++)
0225       LayerEne[ilayer] = 0.;
0226     for (int iring = 0; iring < 100; iring++)
0227       EnRing[iring] = 0.;
0228
0229     for (std::map<int, float>::iterator is = energyInScints.begin(); is != energyInScints.end(); is++) {
0230       ETot = (*is).second;
0231
0232       int layer = org->getlayerID((*is).first);
0233
0234       if ((layer != 17) && (layer != 18)) {
0235         float eta = org->getetaID((*is).first);
0236         float phi = org->getphiID((*is).first);
0237
0238         SEnergy += ETot;
0239         int iphi = static_cast<int>(phi);
0240         int ieta = static_cast<int>(eta);
0241         TowerEne[iphi][ieta] += ETot;
0242
0243         TowerEneCF[iphi][ieta] += ETot * (1. + 0.1 * randGauss.fire());
0244         double dR = 0.08727 * std::sqrt((eta - 8.) * (eta - 8.) + (phi - 3.) * (phi - 3.));
0245         EnRing[static_cast<int>(dR / 0.01)] += ETot;
0246       }
0247
0248       LayerEne[layer] += ETot;
0249     }
0250
0251     for (int ilayer = 0; ilayer < 19; ilayer++) {
0252       histo->fillProfile(ilayer, LayerEne[ilayer] / GeV);
0253     }
0254
0255     for (int iring = 0; iring < 100; iring++)
0256       histo->fillTransProf(0.01 * iring + 0.005, EnRing[iring] / GeV);
0257
0258     for (int iphi = 0; iphi < 8; iphi++) {
0259       for (int jeta = 0; jeta < 18; jeta++) {
0260         SEnergyN += TowerEneCF[iphi][jeta] + 3. * randGauss.fire();  // QGSP
0261
0262         double Rand = 3. * randGauss.fire();  // QGSP
0263
0264         if ((iphi >= 0) && (iphi < 7)) {
0265           if ((jeta >= 5) && (jeta < 12)) {
0266             E7x7Matrix += TowerEne[iphi][jeta];
0267             E7x7MatrixN += TowerEneCF[iphi][jeta] + Rand;
0268
0269             if ((iphi >= 1) && (iphi < 6)) {
0270               if ((jeta >= 6) && (jeta < 11)) {
0271                 E5x5Matrix += TowerEne[iphi][jeta];
0272                 E5x5MatrixN += TowerEneCF[iphi][jeta] + Rand;
0273               }
0274             }
0275           }
0276         }
0277       }
0278     }
0279
0280     //
0281     // Find Primary info:
0282     //
0283     int trackID = 0;
0284     G4PrimaryParticle* thePrim = nullptr;
0285     G4int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
0286 #ifdef EDM_ML_DEBUG
0287     edm::LogVerbatim("HcalTBSim") << "HcalTB02Analysis :: Event has " << nvertex << " vertex";
0288 #endif
0289     if (nvertex == 0)
0290       edm::LogWarning("HcalTBSim") << "HcalTB02Analysis:: End Of Event  "
0291                                    << "ERROR: no vertex";
0292
0293     for (int i = 0; i < nvertex; i++) {
0294       G4PrimaryVertex* avertex = (*evt)()->GetPrimaryVertex(i);
0295       if (avertex == nullptr) {
0296         edm::LogWarning("HcalTBSim") << "HcalTB02Analysis:: End Of Event "
0297                                      << "ERROR: pointer to vertex = 0";
0298       } else {
0299 #ifdef EDM_ML_DEBUG
0300         int npart = avertex->GetNumberOfParticle();
0301         edm::LogVerbatim("HcalTBSim") << "HcalTB02Analysis::Vertex number :" << i << " with " << npart << " particles";
0302 #endif
0303         if (thePrim == nullptr)
0304           thePrim = avertex->GetPrimary(trackID);
0305       }
0306     }
0307
0308     double px = 0., py = 0., pz = 0.;
0309
0310     if (thePrim != nullptr) {
0311       px = thePrim->GetPx();
0312       py = thePrim->GetPy();
0313       pz = thePrim->GetPz();
0314       pInit = std::sqrt(pow(px, 2.) + pow(py, 2.) + pow(pz, 2.));
0315       if (pInit == 0) {
0316         edm::LogWarning("HcalTBSim") << "HcalTB02Analysis:: End Of Event "
0317                                      << " ERROR: primary has p=0 ";
0318       } else {
0319         float costheta = pz / pInit;
0320         float theta = acos(std::min(std::max(costheta, float(-1.)), float(1.)));
0321         eta = -std::log(std::tan(theta / 2));
0322         if (px != 0)
0323           phi = std::atan(py / px);
0324       }
0325       particleType = thePrim->GetPDGcode();
0326     } else {
0327 #ifdef EDM_ML_DEBUG
0328       edm::LogVerbatim("HcalTBSim") << "HcalTB02Analysis:: End Of Event ERROR: could not find primary ";
0329 #endif
0330     }
0331
0332     CaloG4Hit* firstHit = (*theHCHC)[0];
0333     incidentEnergy = (firstHit->getIncidentEnergy() / GeV);
0334
0335   }  // number of Hits > 0
0336
0337   if (!hcalOnly) {
0338     // XTALS
0339
0340     if (XTALSid >= 0 && theXTHC != nullptr) {
0341       for (int xihit = 0; xihit < xentries; xihit++) {
0342         CaloG4Hit* xaHit = (*theXTHC)[xihit];
0343
0344         float xenEm = xaHit->getEM();
0345         float xenhad = xaHit->getHadr();
0346         xtalID = xaHit->getUnitID();
0347
0348         energyInCrystals[xtalID] += xenEm + xenhad;
0349       }
0350
0351       float xCrysEne[7][7];
0352       for (int irow = 0; irow < 7; irow++) {
0353         for (int jcol = 0; jcol < 7; jcol++) {
0354           xCrysEne[irow][jcol] = 0.;
0355         }
0356       }
0357
0358       for (std::map<int, float>::iterator is = energyInCrystals.begin(); is != energyInCrystals.end(); is++) {
0359         int xtalID = (*is).first;
0360         xETot = (*is).second;
0361
0362         int irow = static_cast<int>(xtalID / 100.);
0363         int jcol = static_cast<int>(xtalID - 100. * irow);
0364
0365         xSEnergy += xETot;
0366         xCrysEne[irow][jcol] = xETot;
0367
0368         float dR = std::sqrt(0.01619 * 0.01619 * (jcol - 3) * (jcol - 3) + 0.01606 * 0.01606 * (irow - 3) * (irow - 3));
0369         histo->fillTransProf(dR, xETot * 1.05);
0370
0371         if ((irow > 0) && (irow < 6)) {
0372           if ((jcol > 0) && (jcol < 6)) {
0373             xE5x5Matrix += xCrysEne[irow][jcol];
0374             xE5x5MatrixN += xCrysEne[irow][jcol] + 108.5 * randGauss.fire();
0375
0376             if ((irow > 1) && (irow < 5)) {
0377               if ((jcol > 1) && (jcol < 5)) {
0378                 xE3x3Matrix += xCrysEne[irow][jcol];
0379                 xE3x3MatrixN += xCrysEne[irow][jcol] + 108.5 * randGauss.fire();
0380               }
0381             }
0382           }
0383         }
0384       }
0385
0386       if (!hcalOnly) {
0387         //  assert(theXTHC);
0388         if (theXTHC != nullptr) {
0389           CaloG4Hit* xfirstHit = (*theXTHC)[0];
0390           xIncidentEnergy = xfirstHit->getIncidentEnergy() / GeV;
0391         }
0392       }
0393
0394     }  // number of Hits > 0
0395   }
0396
0397   int iEvt = (*evt)()->GetEventID();
0398   if (iEvt < 10)
0399     edm::LogVerbatim("HcalTBSim") << " Event " << iEvt;
0400   else if ((iEvt < 100) && (iEvt % 10 == 0))
0401     edm::LogVerbatim("HcalTBSim") << " Event " << iEvt;
0402   else if ((iEvt < 1000) && (iEvt % 100 == 0))
0403     edm::LogVerbatim("HcalTBSim") << " Event " << iEvt;
0404   else if ((iEvt < 10000) && (iEvt % 1000 == 0))
0405     edm::LogVerbatim("HcalTBSim") << " Event " << iEvt;
0406 }
0407
0408 void HcalTB02Analysis::fillEvent(HcalTB02HistoClass& product) {
0409   //Beam information
0410   product.set_Nprim(float(primaries.size()));
0411   product.set_partType(particleType);
0412   product.set_Einit(pInit / GeV);
0413   product.set_eta(eta);
0414   product.set_phi(phi);
0415   product.set_Eentry(incidentEnergy);
0416
0417   //Calorimeter energy
0418   product.set_ETot(SEnergy / GeV);
0419   product.set_E7x7(E7x7Matrix / GeV);
0420   product.set_E5x5(E5x5Matrix / GeV);
0421   product.set_ETotN(SEnergyN / GeV);
0422   product.set_E7x7N(E7x7MatrixN / GeV);
0423   product.set_E5x5N(E5x5MatrixN / GeV);
0424   product.set_NUnit(float(energyInScints.size()));
0425   product.set_Ntimesli(float(maxTime));
0426
0427   //crystal information
0428   product.set_xEentry(xIncidentEnergy);
0429   product.set_xNUnit(float(energyInCrystals.size()));
0430   product.set_xETot(xSEnergy / GeV);
0431   product.set_xETotN(xSEnergyN / GeV);
0432   product.set_xE5x5(xE5x5Matrix / GeV);
0433   product.set_xE3x3(xE3x3Matrix / GeV);
0434   product.set_xE5x5N(xE5x5MatrixN / GeV);
0435   product.set_xE3x3N(xE3x3MatrixN / GeV);
0436 }
0437
0438 void HcalTB02Analysis::clear() {
0439   primaries.clear();
0440   particleType = 0;
0441   pInit = eta = phi = incidentEnergy = 0;
0442
0443   SEnergy = E7x7Matrix = E5x5Matrix = SEnergyN = 0;
0444   E7x7MatrixN = E5x5MatrixN = 0;
0445   energyInScints.clear();
0446   maxTime = 0;
0447
0448   xIncidentEnergy = 0;
0449   energyInCrystals.clear();
0450   xSEnergy = xSEnergyN = xE5x5Matrix = xE3x3Matrix = 0;
0451   xE5x5MatrixN = xE3x3MatrixN = 0;
0452 }
0453
0454 void HcalTB02Analysis::finish() {
0455   /*
0456   //Profile
0457   std::ofstream   oFile;
0458   oFile.open("profile.dat");
0459   float st[19] = {0.8,
0460                   0.4,  0.4,  0.4,  0.4,  0.4,
0461                   0.4,  0.4,  0.4,  0.4,  0.4,
0462                   0.4,  0.4,  0.4,  0.4,  0.4,
0463                   0.8,  1.0,  1.0};
0464
0465   //cm of material (brass) in front of scintillator layer i:
0466
0467   float w[19] = {7.45,                         //iron !
0468          6.00, 6.00, 6.00, 6.00, 6.00, //brass
0469          6.00, 6.00, 6.00, 6.00, 6.60, //brass
0470          6.60, 6.60, 6.60, 6.60, 6.60, //brass
0471          8.90, 20.65, 19.5};            //brass,iron !
0472
0473   for (int ilayer = 0; ilayer<19; ilayer++) {
0474
0475     // Histogram mean and sigma calculated from the ROOT histos
0476     edm::LogVerbatim("HcalTBSim") << "Layer number: " << ilayer << " Mean = "
0477                   << histo->getMean(ilayer) << " sigma = "
0478                   << histo->getRMS(ilayer) << " LThick= "
0479                   << w[ilayer] << " SThick= "   << st[ilayer];
0480
0481     oFile << ilayer << " "  << histo->getMean(ilayer) << " "
0482       << histo->getRMS(ilayer)  << " " << w[ilayer] << " " << st[ilayer]
0483       << std::endl;
0484
0485   }
0486   oFile.close();
0487   */
0488 }
0489
0490 DEFINE_SIMWATCHER(HcalTB02Analysis);