Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-08-10 22:49:56

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: ZdcShowerLibrary.cc
0003 // Description: Shower library for the Zero Degree Calorimeter
0004 // E. Garcia June 2008
0005 ///////////////////////////////////////////////////////////////////////////////
0006 
0007 #include "SimG4CMS/Forward/interface/ZdcShowerLibrary.h"
0008 #include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "FWCore/ServiceRegistry/interface/Service.h"
0011 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0012 #include "FWCore/Utilities/interface/Exception.h"
0013 
0014 #include "G4VPhysicalVolume.hh"
0015 #include "G4Step.hh"
0016 #include "G4Track.hh"
0017 #include "Randomize.hh"
0018 #include "CLHEP/Units/GlobalSystemOfUnits.h"
0019 
0020 ZdcShowerLibrary::ZdcShowerLibrary(const std::string& name, edm::ParameterSet const& p) {
0021   edm::ParameterSet m_HS = p.getParameter<edm::ParameterSet>("ZdcShowerLibrary");
0022   verbose = m_HS.getUntrackedParameter<int>("Verbosity", 0);
0023 
0024   npe = 9;  // number of channels or fibers where the energy will be deposited
0025   hits.reserve(npe);
0026 }
0027 
0028 ZdcShowerLibrary::~ZdcShowerLibrary() {}
0029 
0030 std::vector<ZdcShowerLibrary::Hit>& ZdcShowerLibrary::getHits(const G4Step* aStep, bool& ok) {
0031   const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0032   const G4StepPoint* postStepPoint = aStep->GetPostStepPoint();
0033   const G4Track* track = aStep->GetTrack();
0034 
0035   const G4DynamicParticle* aParticle = track->GetDynamicParticle();
0036   const G4ThreeVector& momDir = aParticle->GetMomentumDirection();
0037   double energy = preStepPoint->GetKineticEnergy();
0038   G4ThreeVector hitPoint = preStepPoint->GetPosition();
0039   const G4ThreeVector& hitPointOrig = preStepPoint->GetPosition();
0040   G4int parCode = track->GetDefinition()->GetPDGEncoding();
0041 
0042   hits.clear();
0043 
0044   ok = false;
0045   bool isEM = G4TrackToParticleID::isGammaElectronPositron(parCode);
0046   bool isHad = G4TrackToParticleID::isStableHadronIon(track);
0047   if (!isEM && !isHad)
0048     return hits;
0049   ok = true;
0050 
0051   G4ThreeVector pos;
0052   G4ThreeVector posLocal;
0053   double tSlice = (postStepPoint->GetGlobalTime()) / nanosecond;
0054 
0055   int nHit = 0;
0056   HcalZDCDetId::Section section;
0057   bool side = false;
0058   int channel = 0;
0059   double xx, yy, zz;
0060   double xxlocal, yylocal, zzlocal;
0061 
0062   ZdcShowerLibrary::Hit oneHit;
0063   side = (hitPointOrig.z() > 0.) ? true : false;
0064 
0065   float xWidthEM = std::abs(theXChannelBoundaries[0] - theXChannelBoundaries[1]);
0066   float zWidthEM = std::abs(theZSectionBoundaries[0] - theZSectionBoundaries[1]);
0067   float zWidthHAD = std::abs(theZHadChannelBoundaries[0] - theZHadChannelBoundaries[1]);
0068 
0069   for (int i = 0; i < npe; i++) {
0070     if (i < 5) {
0071       section = HcalZDCDetId::EM;
0072       channel = i + 1;
0073       xxlocal = theXChannelBoundaries[i] + (xWidthEM / 2.);
0074       xx = xxlocal + X0;
0075       yy = 0.0;
0076       yylocal = yy + Y0;
0077       zzlocal = theZSectionBoundaries[0] + (zWidthEM / 2.);
0078       zz = (hitPointOrig.z() > 0.) ? zzlocal + Z0 : zzlocal - Z0;
0079       pos = G4ThreeVector(xx, yy, zz);
0080       posLocal = G4ThreeVector(xxlocal, yylocal, zzlocal);
0081     }
0082     if (i > 4) {
0083       section = HcalZDCDetId::HAD;
0084       channel = i - 4;
0085       xxlocal = 0.0;
0086       xx = xxlocal + X0;
0087       yylocal = 0;
0088       yy = yylocal + Y0;
0089       zzlocal = (hitPointOrig.z() > 0.) ? theZHadChannelBoundaries[i - 5] + (zWidthHAD / 2.)
0090                                         : theZHadChannelBoundaries[i - 5] - (zWidthHAD / 2.);
0091       zz = (hitPointOrig.z() > 0.) ? zzlocal + Z0 : zzlocal - Z0;
0092       pos = G4ThreeVector(xx, yy, zz);
0093       posLocal = G4ThreeVector(xxlocal, yylocal, zzlocal);
0094     }
0095 
0096     oneHit.position = pos;
0097     oneHit.entryLocal = posLocal;
0098     oneHit.depth = channel;
0099     oneHit.time = tSlice;
0100     oneHit.detID = HcalZDCDetId(section, side, channel);
0101 
0102     // Note: coodinates of hit are relative to center of detector (X0,Y0,Z0)
0103     hitPoint.setX(hitPointOrig.x() - X0);
0104     hitPoint.setY(hitPointOrig.y() - Y0);
0105     double setZ = (hitPointOrig.z() > 0.) ? hitPointOrig.z() - Z0 : fabs(hitPointOrig.z()) - Z0;
0106     hitPoint.setZ(setZ);
0107 
0108     int dE = getEnergyFromLibrary(hitPoint, momDir, energy, parCode, section, side, channel);
0109 
0110     if (isEM) {
0111       oneHit.DeEM = dE;
0112       oneHit.DeHad = 0.;
0113     } else {
0114       oneHit.DeEM = 0;
0115       oneHit.DeHad = dE;
0116     }
0117 
0118     hits.push_back(oneHit);
0119 
0120     edm::LogVerbatim("ZdcShower") << "\nZdcShowerLibrary:Generated Hit " << nHit << " orig hit pos " << hitPointOrig
0121                                   << " orig hit pos local coord" << hitPoint << " new position "
0122                                   << (hits[nHit].position) << " Channel " << (hits[nHit].depth) << " side " << side
0123                                   << " Time " << (hits[nHit].time) << " DetectorID " << (hits[nHit].detID)
0124                                   << " Had Energy " << (hits[nHit].DeHad) << " EM Energy  " << (hits[nHit].DeEM)
0125                                   << "\n";
0126     nHit++;
0127   }
0128   return hits;
0129 }
0130 
0131 int ZdcShowerLibrary::getEnergyFromLibrary(const G4ThreeVector& hitPoint,
0132                                            const G4ThreeVector& momDir,
0133                                            double energy,
0134                                            G4int parCode,
0135                                            HcalZDCDetId::Section section,
0136                                            bool side,
0137                                            int channel) {
0138   int nphotons = -1;
0139 
0140   energy = energy / GeV;
0141 
0142   edm::LogVerbatim("ZdcShower") << "\n ZdcShowerLibrary::getEnergyFromLibrary input/output variables:"
0143                                 << " phi: " << 59.2956 * momDir.phi() << " theta: " << 59.2956 * momDir.theta()
0144                                 << " xin : " << hitPoint.x() << " yin : " << hitPoint.y() << " zin : " << hitPoint.z()
0145                                 << " track en: " << energy << "(GeV)"
0146                                 << " section: " << section << " side: " << side << " channel: " << channel
0147                                 << " partID: " << parCode;
0148 
0149   double eav = 0.;
0150   double esig = 0.;
0151   double edis = 0.;
0152 
0153   float xin = hitPoint.x();
0154   float yin = hitPoint.y();
0155   float fact = 0.;
0156 
0157   bool isEM = G4TrackToParticleID::isGammaElectronPositron(parCode);
0158 
0159   if (section == 1 && !isEM) {
0160     if (channel < 5)
0161       if (((theXChannelBoundaries[channel - 1]) < (xin + X0)) && ((xin + X0) <= theXChannelBoundaries[channel]))
0162         fact = 0.18;
0163     if (channel == 5)
0164       if (theXChannelBoundaries[channel - 1] < xin + X0)
0165         fact = 0.18;
0166   }
0167 
0168   if (section == 2 && !isEM) {
0169     if (channel == 1)
0170       fact = 0.34;
0171     if (channel == 2)
0172       fact = 0.24;
0173     if (channel == 3)
0174       fact = 0.17;
0175     if (channel == 4)
0176       fact = 0.07;
0177   }
0178   if (section == 1 && isEM) {
0179     if (channel < 5)
0180       if (((theXChannelBoundaries[channel - 1]) < (xin + X0)) && ((xin + X0) <= theXChannelBoundaries[channel]))
0181         fact = 1.;
0182     if (channel == 5)
0183       if (theXChannelBoundaries[channel - 1] < xin + X0)
0184         fact = 1.0;
0185   }
0186 
0187   //change to cm for parametrization
0188   yin = yin / cm;
0189   xin = xin / cm;
0190 
0191   if (isEM) {
0192     eav = ((((((-0.0002 * xin - 2.0e-13) * xin + 0.0022) * xin + 1.0e-11) * xin - 0.0217) * xin - 3.0e-10) * xin +
0193            1.0028) *
0194           (((0.0001 * yin + 0.0056) * yin + 0.0508) * yin + 1.0) * 300.0 * pow((energy / 300.0), 0.99);  // EM
0195     esig = ((((((0.0005 * xin - 1.0e-12) * xin - 0.0052) * xin + 5.0e-11) * xin + 0.032) * xin - 2.0e-10) * xin + 1.0) *
0196            (((0.0006 * yin + 0.0071) * yin - 0.031) * yin + 1.0) * 30.0 * pow((energy / 300.0), 0.54);  // EM
0197     edis = 1.0;
0198   } else {
0199     eav = ((((((-0.0002 * xin - 2.0e-13) * xin + 0.0022) * xin + 1.0e-11) * xin - 0.0217) * xin - 3.0e-10) * xin +
0200            1.0028) *
0201           (((0.0001 * yin + 0.0056) * yin + 0.0508) * yin + 1.0) * 300.0 * pow((energy / 300.0), 1.12);  // HD
0202     esig = ((((((0.0005 * xin - 1.0e-12) * xin - 0.0052) * xin + 5.0e-11) * xin + 0.032) * xin - 2.0e-10) * xin + 1.0) *
0203            (((0.0006 * yin + 0.0071) * yin - 0.031) * yin + 1.0) * 54.0 * pow((energy / 300.0), 0.93);  //HD
0204     edis = 3.0;
0205   }
0206 
0207   if (eav < 0. || esig < 0.) {
0208     edm::LogVerbatim("ZdcShower") << " Negative everage energy or esigma from parametrization \n"
0209                                   << " xin: " << xin << "(cm)"
0210                                   << " yin: " << yin << "(cm)"
0211                                   << " track en: " << energy << "(GeV)"
0212                                   << " eaverage: " << eav << " (GeV)"
0213                                   << " esigma: " << esig << "  (GeV)"
0214                                   << " edist: " << edis << " (GeV)";
0215     return 0;
0216   }
0217 
0218   // Convert from GeV to MeV for the code
0219   eav = eav * GeV;
0220   esig = esig * GeV;
0221 
0222   while (nphotons == -1 || nphotons > int(eav + 5. * esig))
0223     nphotons = (int)(fact * photonFluctuation(eav, esig, edis));
0224 
0225   edm::LogVerbatim("ZdcShower") << " track en: " << energy << "(GeV)"
0226                                 << " eaverage: " << eav / GeV << " (GeV)"
0227                                 << " esigma: " << esig / GeV << "  (GeV)"
0228                                 << " edist: " << edis << " (GeV)"
0229                                 << " dE hit: " << nphotons / GeV << " (GeV)";
0230 
0231   return nphotons;
0232 }
0233 
0234 int ZdcShowerLibrary::photonFluctuation(double eav, double esig, double edis) {
0235   int nphot = 0;
0236   double efluct = 0.;
0237   if (edis == 1.0)
0238     efluct = eav + esig * CLHEP::RandGaussQ::shoot();
0239   if (edis == 3.0)
0240     efluct = eav + esig * CLHEP::RandLandau::shoot();
0241   nphot = int(efluct);
0242   if (nphot < 0)
0243     nphot = 0;
0244   return nphot;
0245 }