Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-12-12 23:53:08 

0001 ///////////////////////////////////////////////////////////////////////////////

0002 // File: HFShower.cc

0003 // Description: Sensitive Detector class for calorimeters

0004 ///////////////////////////////////////////////////////////////////////////////

0005 
0006 #include "SimG4CMS/Calo/interface/HFShower.h"
0007 #include "SimG4CMS/Calo/interface/HFFibreFiducial.h"
0008 
0009 #include "G4NavigationHistory.hh"
0010 #include "G4Step.hh"
0011 #include "G4Track.hh"
0012 #include "G4VSolid.hh"
0013 #include "Randomize.hh"
0014 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0015 #include "CLHEP/Units/GlobalSystemOfUnits.h"
0016 
0017 //#define EDM_ML_DEBUG

0018 
0019 #include <iostream>
0020 
0021 HFShower::HFShower(const std::string &name,
0022                    const HcalDDDSimConstants *hcons,
0023                    const HcalSimulationParameters *hps,
0024                    edm::ParameterSet const &p,
0025                    int chk)
0026     : cherenkov_(p.getParameter<edm::ParameterSet>("HFShower")), fibre_(hcons, hps, p), chkFibre_(chk) {
0027   edm::ParameterSet m_HF = p.getParameter<edm::ParameterSet>("HFShower");
0028   applyFidCut_ = m_HF.getParameter<bool>("ApplyFiducialCut");
0029   edm::ParameterSet m_HF2 = m_HF.getParameter<edm::ParameterSet>("HFShowerBlock");
0030   equalizeTimeShift_ = m_HF2.getParameter<bool>("EqualizeTimeShift");
0031   probMax_ = m_HF2.getParameter<double>("ProbMax");
0032 
0033   edm::LogVerbatim("HFShower") << "HFShower:: Maximum probability cut off " << probMax_ << " Check flag " << chkFibre_
0034                                << " EqualizeTimeShift " << equalizeTimeShift_;
0035 
0036   //Special Geometry parameters

0037   gpar_ = hcons->getGparHF();
0038 }
0039 
0040 std::vector<HFShower::Hit> HFShower::getHits(const G4Step *aStep, double weight) {
0041   std::vector<HFShower::Hit> hits;
0042   int nHit = 0;
0043 
0044   double edep = weight * (aStep->GetTotalEnergyDeposit());
0045 #ifdef EDM_ML_DEBUG
0046   edm::LogVerbatim("HFShower") << "HFShower::getHits: energy " << aStep->GetTotalEnergyDeposit() << " weight " << weight
0047                                << " edep " << edep;
0048 #endif
0049   double stepl = 0.;
0050 
0051   if (aStep->GetTrack()->GetDefinition()->GetPDGCharge() != 0.)
0052     stepl = aStep->GetStepLength();
0053   if ((edep == 0.) || (stepl == 0.)) {
0054 #ifdef EDM_ML_DEBUG
0055     edm::LogVerbatim("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
0056 #endif
0057     return hits;
0058   }
0059   const G4Track *aTrack = aStep->GetTrack();
0060   const G4DynamicParticle *aParticle = aTrack->GetDynamicParticle();
0061 
0062   HFShower::Hit hit;
0063   double energy = aParticle->GetTotalEnergy();
0064   double momentum = aParticle->GetTotalMomentum();
0065   double pBeta = momentum / energy;
0066   double dose = 0.;
0067   int npeDose = 0;
0068 
0069   const G4ThreeVector &momentumDir = aParticle->GetMomentumDirection();
0070   const G4ParticleDefinition *particleDef = aTrack->GetDefinition();
0071 
0072   const G4StepPoint *preStepPoint = aStep->GetPreStepPoint();
0073   const G4ThreeVector &globalPos = preStepPoint->GetPosition();
0074   G4String name = preStepPoint->GetTouchable()->GetSolid(0)->GetName();
0075   //double zv = std::abs(globalPos.z()) - gpar_[4] - 0.5*gpar_[1];

0076   double zv = std::abs(globalPos.z()) - gpar_[4];
0077   G4ThreeVector localPos = G4ThreeVector(globalPos.x(), globalPos.y(), zv);
0078   G4ThreeVector localMom = preStepPoint->GetTouchable()->GetHistory()->GetTopTransform().TransformAxis(momentumDir);
0079   // @@ Here the depth should be changed  (Fibers all long in Geometry!)

0080   int depth = 1;
0081   int npmt = 0;
0082   bool ok = true;
0083   if (zv < 0. || zv > gpar_[1]) {
0084     ok = false;  // beyond the fiber in Z

0085   }
0086   if (ok && applyFidCut_) {
0087     npmt = HFFibreFiducial::PMTNumber(globalPos);
0088     if (npmt <= 0) {
0089       ok = false;
0090     } else if (npmt > 24) {  // a short fibre

0091       if (zv > gpar_[0]) {
0092         depth = 2;
0093       } else {
0094         ok = false;
0095       }
0096     }
0097 #ifdef EDM_ML_DEBUG
0098     edm::LogVerbatim("HFShower") << "HFShower: npmt " << npmt << " zv " << std::abs(globalPos.z()) << ":" << gpar_[4]
0099                                  << ":" << zv << ":" << gpar_[0] << " ok " << ok << " depth " << depth;
0100 #endif
0101   } else {
0102     depth = (preStepPoint->GetTouchable()->GetReplicaNumber(0)) % 10;  // All LONG!

0103   }
0104   G4ThreeVector translation = preStepPoint->GetTouchable()->GetVolume(1)->GetObjectTranslation();
0105 
0106   double u = localMom.x();
0107   double v = localMom.y();
0108   double w = localMom.z();
0109   double zCoor = localPos.z();
0110   double zFibre = (0.5 * gpar_[1] - zCoor - translation.z());
0111   double tSlice = (aStep->GetPostStepPoint()->GetGlobalTime());
0112   double time =
0113       (equalizeTimeShift_) ? (fibre_.tShift(localPos, depth, -1)) : (fibre_.tShift(localPos, depth, chkFibre_));
0114 
0115 #ifdef EDM_ML_DEBUG
0116   edm::LogVerbatim("HFShower") << "HFShower::getHits: in " << name << " Z " << zCoor << "(" << globalPos.z() << ") "
0117                                << zFibre << " Trans " << translation << " Time " << tSlice << " " << time
0118                                << "\n                  Direction " << momentumDir << " Local " << localMom;
0119 #endif
0120   // Here npe should be 0 if there is no fiber (@@ M.K.)

0121   int npe = 1;
0122   std::vector<double> wavelength;
0123   std::vector<double> momz;
0124   if (!applyFidCut_) {  // _____ Tmp close of the cherenkov function

0125     if (ok)
0126       npe = cherenkov_.computeNPE(aStep, particleDef, pBeta, u, v, w, stepl, zFibre, dose, npeDose);
0127     wavelength = cherenkov_.getWL();
0128     momz = cherenkov_.getMom();
0129   }  // ^^^^^ End of Tmp close of the cherenkov function

0130   if (ok && npe > 0) {
0131     for (int i = 0; i < npe; ++i) {
0132       double p = 1.;
0133       if (!applyFidCut_)
0134         p = fibre_.attLength(wavelength[i]);
0135       double r1 = G4UniformRand();
0136       double r2 = G4UniformRand();
0137 #ifdef EDM_ML_DEBUG
0138       edm::LogVerbatim("HFShower") << "HFShower::getHits: " << i << " attenuation " << r1 << ":" << exp(-p * zFibre)
0139                                    << " r2 " << r2 << ":" << probMax_
0140                                    << " Survive: " << (r1 <= exp(-p * zFibre) && r2 <= probMax_);
0141 #endif
0142       if (applyFidCut_ || chkFibre_ < 0 || (r1 <= exp(-p * zFibre) && r2 <= probMax_)) {
0143         hit.depth = depth;
0144         hit.time = tSlice + time;
0145         // Temporary fix

0146         if (!applyFidCut_) {
0147           hit.wavelength = wavelength[i];
0148           hit.momentum = momz[i];
0149         } else {
0150           hit.wavelength = 300.;
0151           hit.momentum = 1.;
0152         }
0153         hit.position = globalPos;
0154         hits.push_back(hit);
0155         nHit++;
0156       }
0157     }
0158   }
0159 
0160 #ifdef EDM_ML_DEBUG
0161   edm::LogVerbatim("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
0162   for (int i = 0; i < nHit; ++i)
0163     edm::LogVerbatim("HFShower") << "HFShower::Hit " << i << " WaveLength " << hits[i].wavelength << " Time "
0164                                  << hits[i].time << " Momentum " << hits[i].momentum << " Position "
0165                                  << hits[i].position;
0166 #endif
0167   return hits;
0168 }
0169 
0170 std::vector<HFShower::Hit> HFShower::getHits(const G4Step *aStep, bool forLibraryProducer, double zoffset) {
0171   std::vector<HFShower::Hit> hits;
0172   int nHit = 0;
0173 
0174   double edep = aStep->GetTotalEnergyDeposit();
0175   double stepl = 0.;
0176 
0177   if (aStep->GetTrack()->GetDefinition()->GetPDGCharge() != 0.)
0178     stepl = aStep->GetStepLength();
0179   if ((edep == 0.) || (stepl == 0.)) {
0180 #ifdef EDM_ML_DEBUG
0181     edm::LogVerbatim("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
0182 #endif
0183     return hits;
0184   }
0185   const G4Track *aTrack = aStep->GetTrack();
0186   const G4DynamicParticle *aParticle = aTrack->GetDynamicParticle();
0187 
0188   HFShower::Hit hit;
0189   double energy = aParticle->GetTotalEnergy();
0190   double momentum = aParticle->GetTotalMomentum();
0191   double pBeta = momentum / energy;
0192   double dose = 0.;
0193   int npeDose = 0;
0194 
0195   const G4ThreeVector &momentumDir = aParticle->GetMomentumDirection();
0196   G4ParticleDefinition *particleDef = aTrack->GetDefinition();
0197 
0198   G4StepPoint *preStepPoint = aStep->GetPreStepPoint();
0199   const G4ThreeVector &globalPos = preStepPoint->GetPosition();
0200   G4String name = preStepPoint->GetTouchable()->GetSolid(0)->GetName();
0201   double zb = std::abs(globalPos.z() - zoffset);  // from beginning of HF

0202   double zv = zb - .5 * gpar_[1];                 // from center of HF

0203   G4ThreeVector localPos = G4ThreeVector(globalPos.x(), globalPos.y(), zv);
0204   G4ThreeVector localMom = preStepPoint->GetTouchable()->GetHistory()->GetTopTransform().TransformAxis(momentumDir);
0205   bool ok = true;
0206   if (zb < 0. || zb > gpar_[1]) {
0207     ok = false;  // beyond the fiber in Z

0208   }
0209   int depth = (preStepPoint->GetTouchable()->GetReplicaNumber(0)) % 10;  // All LONG!

0210   G4ThreeVector translation = preStepPoint->GetTouchable()->GetVolume(1)->GetObjectTranslation();
0211 
0212   double u = localMom.x();
0213   double v = localMom.y();
0214   double w = localMom.z();
0215   double zCoor = localPos.z();
0216   double zFibre = (0.5 * gpar_[1] - zCoor - translation.z());
0217   double tSlice = (aStep->GetPostStepPoint()->GetGlobalTime());
0218 
0219 #ifdef EDM_ML_DEBUG
0220   double time = (equalizeTimeShift_) ? 0 : fibre_.tShift(localPos, depth, chkFibre_);
0221   edm::LogVerbatim("HFShower") << "HFShower::getHits: in " << name << " Z " << zCoor << "(" << globalPos.z() << ") "
0222                                << zFibre << " Trans " << translation << " Time " << tSlice << " " << time
0223                                << "\n                  Direction " << momentumDir << " Local " << localMom;
0224 #endif
0225 
0226   int npe = 0;
0227   std::vector<double> wavelength;
0228   std::vector<double> momz;
0229   if (ok)
0230     npe = cherenkov_.computeNPE(aStep, particleDef, pBeta, u, v, w, stepl, zFibre, dose, npeDose);
0231   wavelength = cherenkov_.getWL();
0232   momz = cherenkov_.getMom();
0233   if (ok && npe > 0) {
0234     for (int i = 0; i < npe; ++i) {
0235       hit.depth = depth;
0236       //      hit.time = tSlice + time;

0237       hit.time = tSlice;  // only shower time, fiber propagation time to be set on reading library

0238       hit.wavelength = wavelength[i];
0239       hit.momentum = momz[i];
0240       hit.position = globalPos;
0241       hits.push_back(hit);
0242       nHit++;
0243     }
0244   }
0245 
0246 #ifdef EDM_ML_DEBUG
0247   edm::LogVerbatim("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
0248   for (int i = 0; i < nHit; ++i)
0249     edm::LogVerbatim("HFShower") << "HFShower::Hit " << i << " WaveLength " << hits[i].wavelength << " Time "
0250                                  << hits[i].time << " Momentum " << hits[i].momentum << " Position "
0251                                  << hits[i].position;
0252 #endif
0253   return hits;
0254 }
0255 
0256 std::vector<HFShower::Hit> HFShower::getHits(const G4Step *aStep, bool forLibrary) {
0257   std::vector<HFShower::Hit> hits;
0258   int nHit = 0;
0259 
0260   double edep = aStep->GetTotalEnergyDeposit();
0261   double stepl = 0.;
0262 
0263   if (aStep->GetTrack()->GetDefinition()->GetPDGCharge() != 0.)
0264     stepl = aStep->GetStepLength();
0265   if ((edep == 0.) || (stepl == 0.)) {
0266 #ifdef EDM_ML_DEBUG
0267     edm::LogVerbatim("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
0268 #endif
0269     return hits;
0270   }
0271 
0272   const G4Track *aTrack = aStep->GetTrack();
0273   const G4DynamicParticle *aParticle = aTrack->GetDynamicParticle();
0274 
0275   HFShower::Hit hit;
0276   double energy = aParticle->GetTotalEnergy();
0277   double momentum = aParticle->GetTotalMomentum();
0278   double pBeta = momentum / energy;
0279   double dose = 0.;
0280   int npeDose = 0;
0281 
0282   const G4ThreeVector &momentumDir = aParticle->GetMomentumDirection();
0283   G4ParticleDefinition *particleDef = aTrack->GetDefinition();
0284 
0285   const G4StepPoint *preStepPoint = aStep->GetPreStepPoint();
0286   const G4ThreeVector &globalPos = preStepPoint->GetPosition();
0287   G4String name = preStepPoint->GetTouchable()->GetSolid(0)->GetName();
0288   double zv = std::abs(globalPos.z()) - gpar_[4] - 0.5 * gpar_[1];
0289   G4ThreeVector localPos = G4ThreeVector(globalPos.x(), globalPos.y(), zv);
0290   G4ThreeVector localMom = preStepPoint->GetTouchable()->GetHistory()->GetTopTransform().TransformAxis(momentumDir);
0291   // @@ Here the depth should be changed (Fibers are all long in Geometry!)

0292   int depth = 1;
0293   int npmt = 0;
0294   bool ok = true;
0295   if (zv < 0 || zv > gpar_[1]) {
0296     ok = false;  // beyond the fiber in Z

0297   }
0298   if (ok && applyFidCut_) {
0299     npmt = HFFibreFiducial::PMTNumber(globalPos);
0300     if (npmt <= 0) {
0301       ok = false;
0302     } else if (npmt > 24) {  // a short fibre

0303       if (zv > gpar_[0]) {
0304         depth = 2;
0305       } else {
0306         ok = false;
0307       }
0308     }
0309 #ifdef EDM_ML_DEBUG
0310     edm::LogVerbatim("HFShower") << "HFShower:getHits(SL): npmt " << npmt << " zv " << std::abs(globalPos.z()) << ":"
0311                                  << gpar_[4] << ":" << zv << ":" << gpar_[0] << " ok " << ok << " depth " << depth;
0312 #endif
0313   } else {
0314     depth = (preStepPoint->GetTouchable()->GetReplicaNumber(0)) % 10;  // All LONG!

0315   }
0316   G4ThreeVector translation = preStepPoint->GetTouchable()->GetVolume(1)->GetObjectTranslation();
0317 
0318   double u = localMom.x();
0319   double v = localMom.y();
0320   double w = localMom.z();
0321   double zCoor = localPos.z();
0322   double zFibre = (0.5 * gpar_[1] - zCoor - translation.z());
0323   double tSlice = (aStep->GetPostStepPoint()->GetGlobalTime());
0324   double time = (equalizeTimeShift_) ? 0 : fibre_.tShift(localPos, depth, chkFibre_);
0325 
0326 #ifdef EDM_ML_DEBUG
0327   edm::LogVerbatim("HFShower") << "HFShower::getHits(SL): in " << name << " Z " << zCoor << "(" << globalPos.z() << ") "
0328                                << zFibre << " Trans " << translation << " Time " << tSlice << " " << time
0329                                << "\n                  Direction " << momentumDir << " Local " << localMom;
0330 #endif
0331   // npe should be 0

0332   int npe = 0;
0333   if (ok)
0334     npe = cherenkov_.computeNPE(aStep, particleDef, pBeta, u, v, w, stepl, zFibre, dose, npeDose);
0335   std::vector<double> wavelength = cherenkov_.getWL();
0336   std::vector<double> momz = cherenkov_.getMom();
0337 
0338   for (int i = 0; i < npe; ++i) {
0339     hit.depth = depth;
0340     hit.time = tSlice + time;
0341     hit.wavelength = wavelength[i];
0342     hit.momentum = momz[i];
0343     hit.position = globalPos;
0344     hits.push_back(hit);
0345     nHit++;
0346   }
0347 
0348   return hits;
0349 }
[ 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!