Back to home page

Project CMSSW displayed by LXR
 
 

    

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

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