29 #include "CLHEP/Units/SystemOfUnits.h"
30 #include "CLHEP/Units/PhysicalConstants.h"
33 #include "HepPDT/ParticleData.hh"
38 #include "CLHEP/Random/RandPoisson.h"
39 #include "CLHEP/Random/RandFlat.h"
40 #include "CLHEP/Random/RandBinomial.h"
41 #include "CLHEP/Random/RandExpZiggurat.h"
42 #include "CLHEP/Random/RandGaussZiggurat.h"
55 const HepPDT::ParticleDataTable* pdt,
64 m_pPAItoolXe(paitoolXe),
65 m_pPAItoolAr(paitoolAr),
66 m_pPAItoolKr(paitoolKr),
67 m_pSimDriftTimeTool(simdrifttool),
70 m_pTimeCorrection(nullptr),
71 m_pElectronicsProcessing(ep),
73 m_pDigConditions(digcond),
74 m_pParticleTable(pdt),
75 m_alreadywarnedagainstpdg0(false),
104 ATH_MSG_FATAL (
"TRT_PAITool for Xenon not defined! no point in continuing!" );
107 ATH_MSG_ERROR (
"TRT_PAITool for Krypton is not defined!!! Xenon TRT_PAITool will be used for Krypton straws!" );
111 ATH_MSG_ERROR (
"TRT_PAITool for Argon is not defined!!! Xenon TRT_PAITool will be used for Argon straws!" );
129 for (
unsigned int k=0;
k<150;
k++) {
130 double dist = 10.0*(
k+0.5);
141 for (
unsigned int iwheel = 0; iwheel <
num->getNEndcapWheels(); ++iwheel)
143 for (
unsigned int iside = 0; iside < 2; ++iside)
145 for (
unsigned int ilayer = 0; ilayer <
num->getNEndcapLayers(iwheel); ++ilayer)
155 ATH_MSG_VERBOSE (
"Failed to retrieve endcap element for (iside,iwheel,ilayer)=("
156 << iside<<
", "<<iwheel<<
", "<<ilayer<<
")." );
165 const double zcoordfrac((
v.z()>0?
v.z():-
v.z())/
v.mag());
166 if (zcoordfrac<0.98 || zcoordfrac > 1.02)
168 ATH_MSG_WARNING (
"Found endcap straw where the assumption that local x-direction"
169 <<
" is parallel to global z-direction is NOT valid."
170 <<
" Drift times will be somewhat off." );
177 for (
unsigned int phi_it = 0; phi_it < 32; phi_it++)
186 const double coordfrac(atan2(
v.x(),
v.y()));
187 if (coordfrac>0.2 || coordfrac < -0.2)
189 ATH_MSG_WARNING (
"Found barrel straw where the assumption that local y-direction"
190 <<
" is along the straw is NOT valid."
191 <<
" Drift times will be somewhat off." );
205 const double& timeOfHit,
206 const double& prex,
const double& prey,
const double& prez,
207 const double& postx,
const double& posty,
const double& postz,
208 std::vector<cluster>& clusterlist,
int strawGasType,
209 CLHEP::HepRandomEngine* rndmEngine,
210 CLHEP::HepRandomEngine* paiRndmEngine)
216 else if (strawGasType==1) { activePAITool =
m_pPAItoolKr; }
217 else if (strawGasType==2) { activePAITool =
m_pPAItoolAr; }
221 const double deltaX(postx - prex);
222 const double deltaY(posty - prey);
223 const double deltaZ(postz - prez);
226 const double meanFreePath(activePAITool->
GetMeanFreePath( scaledKineticEnergy, particleCharge*particleCharge ));
229 const unsigned int numberOfClusters(CLHEP::RandPoisson::shoot(rndmEngine,stepLength / meanFreePath));
233 for (
unsigned int iclus(0); iclus<numberOfClusters; ++iclus)
236 const double lambda(CLHEP::RandFlat::shoot(rndmEngine));
239 double clusE(activePAITool->
GetEnergyTransfer(scaledKineticEnergy, paiRndmEngine));
240 clusterlist.emplace_back(clusE, timeOfHit,
254 bool & alreadyPrintedPDGcodeWarning,
255 double cosmicEventPhase,
257 bool emulationArflag,
258 bool emulationKrflag,
259 CLHEP::HepRandomEngine* rndmEngine,
260 CLHEP::HepRandomEngine* elecProcRndmEngine,
261 CLHEP::HepRandomEngine* elecNoiseRndmEngine,
262 CLHEP::HepRandomEngine* paiRndmEngine)
268 const int hitID((*i)->GetHitID());
274 const bool isECA (region==3);
275 const bool isECB (region==4);
310 double timeOfHit(0.0);
312 const double globalHitTime(
hitTime(*theHit));
313 const double globalTime =
static_cast<double>((*theHit)->GetGlobalTime());
314 const double bunchCrossingTime(globalHitTime - globalTime);
322 const int particleEncoding((*theHit)->GetParticleEncoding());
325 if (particleEncoding == 0)
333 ATH_MSG_WARNING (
"Ignoring sim. particle with pdgcode 0. This warning is only shown once per job" );
352 double ArEmulationScaling_BA = 0.05;
353 double ArEmulationScaling_ECA = 0.20;
354 double ArEmulationScaling_ECB = 0.20;
357 double KrEmulationScaling_BA = 0.20;
358 double KrEmulationScaling_ECA = 0.39;
359 double KrEmulationScaling_ECB = 0.39;
363 double hitx = TRThitGlobalPos[0];
364 double hity = TRThitGlobalPos[1];
365 double hitz = TRThitGlobalPos[2];
366 double hitEta = std::abs(
log(
tan(0.5*atan2(sqrt(hitx*hitx+hity*hity),hitz))));
367 if ( hitEta < 0.5 ) { trEfficiencyBarrel *= ( 0.833333+0.6666667*hitEta*hitEta ); }
369 if ( strawGasType == 0 && emulationArflag ) { trEfficiencyBarrel *= ArEmulationScaling_BA; }
370 if ( strawGasType == 0 && emulationKrflag ) { trEfficiencyBarrel *= KrEmulationScaling_BA; }
371 if ( CLHEP::RandFlat::shoot(rndmEngine) > trEfficiencyBarrel )
continue;
377 if ( strawGasType == 0 && emulationArflag ) { trEfficiencyEndCapA *= ArEmulationScaling_ECA; }
378 if ( strawGasType == 0 && emulationKrflag ) { trEfficiencyEndCapA *= KrEmulationScaling_ECA; }
379 if ( CLHEP::RandFlat::shoot(rndmEngine) > trEfficiencyEndCapA )
continue;
384 if ( strawGasType == 0 && emulationArflag ) { trEfficiencyEndCapB *= ArEmulationScaling_ECB; }
385 if ( strawGasType == 0 && emulationKrflag ) { trEfficiencyEndCapB *= KrEmulationScaling_ECB; }
386 if ( CLHEP::RandFlat::shoot(rndmEngine) > trEfficiencyEndCapB )
continue;
401 m_clusterlist.emplace_back( (*theHit)->GetEnergyDeposit()*
CLHEP::keV, timeOfHit, (*theHit)->GetPostStepX(), (*theHit)->GetPostStepY(), (*theHit)->GetPostStepZ() );
407 double particleCharge(0.);
408 double particleMass(0.);
411 particleCharge =
particle->charge();
412 particleMass =
particle->mass().value();
421 ATH_MSG_WARNING (
"Data for sim. particle with pdgcode "<<particleEncoding
422 <<
" is not a nucleus and could not be retrieved from PartPropSvc. Assuming mass and charge as pion. Please investigate." );
427 particleCharge =
MC::charge(particleEncoding);
433 particleMass = std::abs(
Z*Mp+(
A-
Z)*Mn );
435 if (!alreadyPrintedPDGcodeWarning) {
436 ATH_MSG_WARNING (
"Data for sim. particle with pdgcode "<<particleEncoding
437 <<
" could not be retrieved from PartPropSvc (unexpected ion)."
438 <<
" Please Investigate the PDGTABLE.MeV file."
439 <<
" Calculating mass and charge from pdg code."
440 <<
" The result is: Charge = "<<particleCharge<<
" Mass = "<<particleMass<<
"MeV" );
441 alreadyPrintedPDGcodeWarning =
true;
446 if (!particleCharge) {
continue; }
449 ATH_MSG_WARNING (
"Ignoring ionization from particle with pdg code "<<particleEncoding
450 <<
" since it appears to be a massless charged particle. Please investigate." );
458 const double scaledKineticEnergy(
static_cast<double>((*theHit)->GetKineticEnergy()) * (
CLHEP::proton_mass_c2 / particleMass ));
461 (*theHit)->GetPreStepX(),(*theHit)->GetPreStepY(),(*theHit)->GetPreStepZ(),
462 (*theHit)->GetPostStepX(),(*theHit)->GetPostStepY(),(*theHit)->GetPostStepZ(),
508 double lowthreshold, noiseamplitude;
513 noiseamplitude = 0.0;
522 const std::vector<cluster>&
clusters,
523 std::vector<TRTElectronicsProcessing::Deposit>& deposits,
525 double cosmicEventPhase,
527 CLHEP::HepRandomEngine* rndmEngine)
539 const bool isShort (region==1);
540 const bool isLong (region==2);
553 std::vector<cluster>::const_iterator currentClusterIter(
clusters.begin());
554 const std::vector<cluster>::const_iterator endOfClusterList(
clusters.end());
566 double map_x2(0.),map_y2(0.),map_z2(0.);
567 double effectiveField2(0.);
572 globalPosition[0]=TRThitGlobalPos[0]*
CLHEP::mm;
573 globalPosition[1]=TRThitGlobalPos[1]*
CLHEP::mm;
574 globalPosition[2]=TRThitGlobalPos[2]*
CLHEP::mm;
577 fieldCache.
getField (globalPosition.data(), mField.data());
579 map_x2 = mField.x()*mField.x();
580 map_y2 = mField.y()*mField.y();
581 map_z2 = mField.z()*mField.z();
608 for (;currentClusterIter!=endOfClusterList;++currentClusterIter)
611 const double cluster_x(currentClusterIter->xpos);
612 const double cluster_y(currentClusterIter->ypos);
613 const double cluster_z(this->
setClusterZ(currentClusterIter->zpos, isLong, isShort, isEC));
614 const double cluster_x2(cluster_x*cluster_x);
615 const double cluster_y2(cluster_y*cluster_y);
616 double cluster_r2(cluster_x2+cluster_y2);
619 if (cluster_r2<wire_r2) cluster_r2=wire_r2;
620 if (cluster_r2>straw_r2) cluster_r2=straw_r2;
622 const double cluster_r(std::sqrt(cluster_r2));
623 const double cluster_E(currentClusterIter->energy);
624 const unsigned int nprimaryelectrons(
static_cast<unsigned int>( cluster_E / ionisationPotential + 1.0 ) );
629 double depositEnergy(0.);
633 unsigned int nsurvivingprimaryelectrons =
static_cast<unsigned int>(randBinomial->fire(rndmEngine,nprimaryelectrons) + 0.5);
634 if (nsurvivingprimaryelectrons==0)
continue;
635 const double meanElectronEnergy(ionisationPotential / smearingFactor);
636 for (
unsigned int ielec(0); ielec<nsurvivingprimaryelectrons; ++ielec) {
637 depositEnergy += CLHEP::RandExpZiggurat::shoot(rndmEngine, meanElectronEnergy);
642 const double fluctSigma(sqrt(cluster_E * ionisationPotential * (2 - smearingFactor) / smearingFactor));
644 depositEnergy = CLHEP::RandGaussZiggurat::shoot(rndmEngine, cluster_E, fluctSigma);
645 }
while(depositEnergy<0.0);
660 effectiveField2 = map_z2*cluster_y2/cluster_r2 + map_x2 + map_y2;
669 effectiveField2 = map_z2 + (map_x2+map_y2)*cluster_y2/cluster_r2;
685 double clusterTime(currentClusterIter->time);
695 double timedirect(0.), timereflect(0.);
699 double expdirect(1.0), expreflect(1.0);
713 const unsigned int kdirect =
static_cast<unsigned int>(distdirect/10);
714 const unsigned int kreflect =
static_cast<unsigned int>(distreflect/10);
721 double dt = clusterTime + commondrifttime;
722 deposits.emplace_back(0.5*depositEnergy*expdirect, timedirect+
dt);
723 deposits.emplace_back(0.5*depositEnergy*expreflect, timereflect+
dt);
732 const int mask(0x0000001F);
734 int trtID, ringID, moduleID, layerID, strawID;
735 int wheelID, planeID, sectorID;
740 if ( !(hitID & 0x00200000) ) {
742 strawID = hitID &
mask;
743 hitID >>= word_shift;
744 layerID = hitID &
mask;
745 hitID >>= word_shift;
746 moduleID = hitID &
mask;
747 hitID >>= word_shift;
748 ringID = hitID &
mask;
749 trtID = hitID >> word_shift;
754 const Amg::Vector3D v( (*theHit)->GetPreStepX(),(*theHit)->GetPreStepY(),(*theHit)->GetPreStepZ());
760 strawID = hitID &
mask;
761 hitID >>= word_shift;
762 planeID = hitID &
mask;
763 hitID >>= word_shift;
764 sectorID = hitID &
mask;
765 hitID >>= word_shift;
766 wheelID = hitID &
mask;
767 trtID = hitID >> word_shift;
770 if (trtID == 3) trtID = 0;
775 if ( endcapElement ) {
776 const Amg::Vector3D v( (*theHit)->GetPreStepX(),(*theHit)->GetPreStepY(),(*theHit)->GetPreStepZ());
782 ATH_MSG_WARNING (
"Could not find global coordinate of a straw - drifttime calculation will be inaccurate" );
783 return {0.0,0.0,0.0};
789 double cluster_z(cluster_z_in);
794 const double longBarrelStrawHalfLength(349.315*
CLHEP::mm);
795 const double shortBarrelStrawHalfLength(153.375*
CLHEP::mm);
796 const double EndcapStrawHalfLength(177.150*
CLHEP::mm);
797 if ( isLong && std::abs(cluster_z)>longBarrelStrawHalfLength+30 ) {
798 double d = cluster_z<0 ? cluster_z+longBarrelStrawHalfLength : cluster_z-longBarrelStrawHalfLength;
799 ATH_MSG_WARNING (
"Long barrel straw cluster is outside the active gas volume z = +- 349.315 mm by " <<
d <<
" mm.");
803 if ( isShort && std::abs(cluster_z)>shortBarrelStrawHalfLength+30 ) {
804 double d = cluster_z<0 ? cluster_z+shortBarrelStrawHalfLength : cluster_z-shortBarrelStrawHalfLength;
805 ATH_MSG_WARNING (
"Short barrel straw cluster is outside the active gas volume z = +- 153.375 mm by " <<
d <<
" mm.");
809 if ( isEC && std::abs(cluster_z)>EndcapStrawHalfLength+30 ) {
810 double d = cluster_z<0 ? cluster_z+EndcapStrawHalfLength : cluster_z-EndcapStrawHalfLength;
811 ATH_MSG_WARNING (
"End cap straw cluster is outside the active gas volume z = +- 177.150 mm by " <<
d <<
" mm.");