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