Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-05-14 22:52:34 

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

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

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

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

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

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

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

0298   int depth = 1;
0299   int npmt = 0;
0300   bool ok = true;
0301   if (zv < 0 || zv > gpar_[1]) {
0302     ok = false;  // beyond the fiber in Z

0303   }
0304   if (ok && applyFidCut_) {
0305     npmt = HFFibreFiducial::PMTNumber(globalPos);
0306     if (npmt <= 0) {
0307       ok = false;
0308     } else if (npmt > 24) {  // a short fibre

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

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

0338   int npe = 0;
0339   if (ok)
0340     npe = cherenkov_.computeNPE(aStep, particleDef, pBeta, u, v, w, stepl, zFibre, dose, npeDose);
0341   std::vector<double> wavelength = cherenkov_.getWL();
0342   std::vector<double> momz = cherenkov_.getMom();
0343 
0344   for (int i = 0; i < npe; ++i) {
0345     hit.depth = depth;
0346     hit.time = tSlice + time;
0347     hit.wavelength = wavelength[i];
0348     hit.momentum = momz[i];
0349     hit.position = globalPos;
0350     hits.push_back(hit);
0351 #ifdef EDM_ML_DEBUG
0352     nHit++;
0353 #endif
0354   }
0355 
0356   return hits;
0357 }
[ 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!