Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-09-30 01:18:28 

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     : hcalConstant_(hcons), 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   cherenkov_ = std::make_unique<HFCherenkov>(m_HF);
0037   fibre_ = std::make_unique<HFFibre>(name, hcalConstant_, hps, p);
0038 
0039   //Special Geometry parameters

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

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

0085   int depth = 1;
0086   int npmt = 0;
0087   bool ok = true;
0088   if (zv < 0. || zv > gpar_[1]) {
0089     ok = false;  // beyond the fiber in Z

0090   }
0091   if (ok && applyFidCut_) {
0092     npmt = HFFibreFiducial::PMTNumber(globalPos);
0093     if (npmt <= 0) {
0094       ok = false;
0095     } else if (npmt > 24) {  // a short fibre

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

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

0126   int npe = 1;
0127   std::vector<double> wavelength;
0128   std::vector<double> momz;
0129   if (!applyFidCut_) {  // _____ Tmp close of the cherenkov function

0130     if (ok)
0131       npe = cherenkov_->computeNPE(aStep, particleDef, pBeta, u, v, w, stepl, zFibre, dose, npeDose);
0132     wavelength = cherenkov_->getWL();
0133     momz = cherenkov_->getMom();
0134   }  // ^^^^^ End of Tmp close of the cherenkov function

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

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

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

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

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

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

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

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

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

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

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

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

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