29 #include "CLHEP/Units/SystemOfUnits.h"
30 #include "CLHEP/Units/PhysicalConstants.h"
33 #include "HepPDT/ParticleData.hh"
36 #include "CLHEP/Random/RandPoisson.h"
37 #include "CLHEP/Random/RandFlat.h"
38 #include "CLHEP/Random/RandBinomial.h"
39 #include "CLHEP/Random/RandExpZiggurat.h"
40 #include "CLHEP/Random/RandGaussZiggurat.h"
53 const HepPDT::ParticleDataTable* pdt,
62 m_pPAItoolXe(paitoolXe),
63 m_pPAItoolAr(paitoolAr),
64 m_pPAItoolKr(paitoolKr),
65 m_pSimDriftTimeTool(simdrifttool),
68 m_pTimeCorrection(nullptr),
69 m_pElectronicsProcessing(ep),
71 m_pDigConditions(digcond),
72 m_pParticleTable(pdt),
73 m_alreadywarnedagainstpdg0(false),
102 ATH_MSG_FATAL (
"TRT_PAITool for Xenon not defined! no point in continuing!" );
105 ATH_MSG_ERROR (
"TRT_PAITool for Krypton is not defined!!! Xenon TRT_PAITool will be used for Krypton straws!" );
109 ATH_MSG_ERROR (
"TRT_PAITool for Argon is not defined!!! Xenon TRT_PAITool will be used for Argon straws!" );
127 for (
unsigned int k=0;
k<150;
k++) {
128 double dist = 10.0*(
k+0.5);
139 for (
unsigned int iwheel = 0; iwheel <
num->getNEndcapWheels(); ++iwheel)
141 for (
unsigned int iside = 0; iside < 2; ++iside)
143 for (
unsigned int ilayer = 0; ilayer <
num->getNEndcapLayers(iwheel); ++ilayer)
153 ATH_MSG_VERBOSE (
"Failed to retrieve endcap element for (iside,iwheel,ilayer)=("
154 << iside<<
", "<<iwheel<<
", "<<ilayer<<
")." );
163 const double zcoordfrac((
v.z()>0?
v.z():-
v.z())/
v.mag());
164 if (zcoordfrac<0.98 || zcoordfrac > 1.02)
166 ATH_MSG_WARNING (
"Found endcap straw where the assumption that local x-direction"
167 <<
" is parallel to global z-direction is NOT valid."
168 <<
" Drift times will be somewhat off." );
175 for (
unsigned int phi_it = 0; phi_it < 32; phi_it++)
184 const double coordfrac(atan2(
v.x(),
v.y()));
185 if (coordfrac>0.2 || coordfrac < -0.2)
187 ATH_MSG_WARNING (
"Found barrel straw where the assumption that local y-direction"
188 <<
" is along the straw is NOT valid."
189 <<
" Drift times will be somewhat off." );
203 const double& timeOfHit,
204 const double& prex,
const double& prey,
const double& prez,
205 const double& postx,
const double& posty,
const double& postz,
206 std::vector<cluster>& clusterlist,
int strawGasType,
207 CLHEP::HepRandomEngine* rndmEngine,
208 CLHEP::HepRandomEngine* paiRndmEngine)
214 else if (strawGasType==1) { activePAITool =
m_pPAItoolKr; }
215 else if (strawGasType==2) { activePAITool =
m_pPAItoolAr; }
219 const double deltaX(postx - prex);
220 const double deltaY(posty - prey);
221 const double deltaZ(postz - prez);
224 const double meanFreePath(activePAITool->
GetMeanFreePath( scaledKineticEnergy, particleCharge*particleCharge ));
227 const unsigned int numberOfClusters(CLHEP::RandPoisson::shoot(rndmEngine,stepLength / meanFreePath));
231 for (
unsigned int iclus(0); iclus<numberOfClusters; ++iclus)
234 const double lambda(CLHEP::RandFlat::shoot(rndmEngine));
237 double clusE(activePAITool->
GetEnergyTransfer(scaledKineticEnergy, paiRndmEngine));
238 clusterlist.emplace_back(clusE, timeOfHit,
252 bool & alreadyPrintedPDGcodeWarning,
253 double cosmicEventPhase,
255 bool emulationArflag,
256 bool emulationKrflag,
257 CLHEP::HepRandomEngine* rndmEngine,
258 CLHEP::HepRandomEngine* elecProcRndmEngine,
259 CLHEP::HepRandomEngine* elecNoiseRndmEngine,
260 CLHEP::HepRandomEngine* paiRndmEngine)
266 const int hitID((*i)->GetHitID());
272 const bool isECA (region==3);
273 const bool isECB (region==4);
308 double timeOfHit(0.0);
310 const double globalHitTime(
hitTime(*theHit));
311 const double globalTime =
static_cast<double>((*theHit)->GetGlobalTime());
312 const double bunchCrossingTime(globalHitTime - globalTime);
320 const int particleEncoding((*theHit)->GetParticleEncoding());
323 if (particleEncoding == 0)
331 ATH_MSG_WARNING (
"Ignoring sim. particle with pdgcode 0. This warning is only shown once per job" );
338 if (particleEncoding == 22)
351 double ArEmulationScaling_BA = 0.05;
352 double ArEmulationScaling_ECA = 0.20;
353 double ArEmulationScaling_ECB = 0.20;
356 double KrEmulationScaling_BA = 0.20;
357 double KrEmulationScaling_ECA = 0.39;
358 double KrEmulationScaling_ECB = 0.39;
362 double hitx = TRThitGlobalPos[0];
363 double hity = TRThitGlobalPos[1];
364 double hitz = TRThitGlobalPos[2];
365 double hitEta = std::abs(
log(
tan(0.5*atan2(sqrt(hitx*hitx+hity*hity),hitz))));
366 if ( hitEta < 0.5 ) { trEfficiencyBarrel *= ( 0.833333+0.6666667*hitEta*hitEta ); }
368 if ( strawGasType == 0 && emulationArflag ) { trEfficiencyBarrel *= ArEmulationScaling_BA; }
369 if ( strawGasType == 0 && emulationKrflag ) { trEfficiencyBarrel *= KrEmulationScaling_BA; }
370 if ( CLHEP::RandFlat::shoot(rndmEngine) > trEfficiencyBarrel )
continue;
376 if ( strawGasType == 0 && emulationArflag ) { trEfficiencyEndCapA *= ArEmulationScaling_ECA; }
377 if ( strawGasType == 0 && emulationKrflag ) { trEfficiencyEndCapA *= KrEmulationScaling_ECA; }
378 if ( CLHEP::RandFlat::shoot(rndmEngine) > trEfficiencyEndCapA )
continue;
383 if ( strawGasType == 0 && emulationArflag ) { trEfficiencyEndCapB *= ArEmulationScaling_ECB; }
384 if ( strawGasType == 0 && emulationKrflag ) { trEfficiencyEndCapB *= KrEmulationScaling_ECB; }
385 if ( CLHEP::RandFlat::shoot(rndmEngine) > trEfficiencyEndCapB )
continue;
398 else if ( ((abs(particleEncoding)/100000==41) && (abs(particleEncoding)/10000000==0)) ||
399 ((
static_cast<int>(abs(particleEncoding)/10000000) == 1) &&
400 (
static_cast<int>(abs(particleEncoding)/100000) == 100) &&
401 (
static_cast<int>((abs(particleEncoding))-10000000)/100>10)) )
403 m_clusterlist.emplace_back((*theHit)->GetEnergyDeposit()*
CLHEP::keV, timeOfHit, (*theHit)->GetPostStepX(), (*theHit)->GetPostStepY(), (*theHit)->GetPostStepZ()
410 double particleCharge(0.);
411 double particleMass(0.);
415 particleCharge =
particle->charge();
416 particleMass =
particle->mass().value();
417 if ((
static_cast<int>(abs(particleEncoding)/10000000) == 1) && (
static_cast<int>(abs(particleEncoding)/100000)==100))
419 particleCharge = (particleEncoding>0 ? 1. : -1.) *(((abs(particleEncoding) / 100000.0) - 100.0) * 1000.0);
421 else if ((
static_cast<int>(abs(particleEncoding)/10000000) == 2) && (
static_cast<int>(abs(particleEncoding)/100000)==200))
423 particleCharge = (particleEncoding>0 ? 1. : -1.) *((
double)((abs(particleEncoding) / 1000) % 100) / (
double)((abs(particleEncoding) / 10) % 100));
428 const int number_of_digits(
static_cast<int>(log10((
double)abs(particleEncoding))+1.));
429 if (number_of_digits != 10)
431 ATH_MSG_ERROR (
"Data for sim. particle with pdgcode "<<particleEncoding
432 <<
" does not have 10 digits and could not be retrieved from PartPropSvc. Assuming mass and charge as pion." );
436 else if (( number_of_digits == 10 ) && (
static_cast<int>(abs(particleEncoding)/100000000)!=10) )
438 ATH_MSG_ERROR (
"Data for sim. particle with pdgcode "<<particleEncoding
439 <<
" has 10 digits, could not be retrieved from PartPropSvc, and is inconsistent with ion pdg convention (+/-10LZZZAAAI)."
440 <<
" Assuming mass and charge as pion." );
444 else if ((number_of_digits==10) && (
static_cast<int>(abs(particleEncoding)/100000000)==10))
446 const int A(
static_cast<int>((((particleEncoding)%1000000000)%10000)/10.));
447 const int Z(
static_cast<int>((((particleEncoding)%10000000)-(((particleEncoding)%1000000000)%10000))/10000.));
452 particleCharge = (particleEncoding>0 ? 1. : -1.) *
static_cast<double>(
Z);
453 particleMass = std::abs(
Z*Mp+(
A-
Z)*Mn );
455 if (!alreadyPrintedPDGcodeWarning)
457 ATH_MSG_WARNING (
"Data for sim. particle with pdgcode "<<particleEncoding
458 <<
" could not be retrieved from PartPropSvc (unexpected ion)."
459 <<
" Calculating mass and charge from pdg code. "
460 <<
" The result is: Charge = "<<particleCharge<<
" Mass = "<<particleMass<<
"MeV" );
461 alreadyPrintedPDGcodeWarning =
true;
472 ATH_MSG_WARNING (
"Ignoring ionization from sim. particle with pdgcode "<<particleEncoding
473 <<
" since it appears to be a massless charged particle." );
481 const double scaledKineticEnergy(
static_cast<double>((*theHit)->GetKineticEnergy()) * (
CLHEP::proton_mass_c2 / particleMass ));
484 (*theHit)->GetPreStepX(),(*theHit)->GetPreStepY(),(*theHit)->GetPreStepZ(),
485 (*theHit)->GetPostStepX(),(*theHit)->GetPostStepY(),(*theHit)->GetPostStepZ(),
531 double lowthreshold, noiseamplitude;
536 noiseamplitude = 0.0;
545 const std::vector<cluster>&
clusters,
546 std::vector<TRTElectronicsProcessing::Deposit>& deposits,
548 double cosmicEventPhase,
550 CLHEP::HepRandomEngine* rndmEngine)
562 const bool isShort (region==1);
563 const bool isLong (region==2);
576 std::vector<cluster>::const_iterator currentClusterIter(
clusters.begin());
577 const std::vector<cluster>::const_iterator endOfClusterList(
clusters.end());
589 double map_x2(0.),map_y2(0.),map_z2(0.);
590 double effectiveField2(0.);
595 globalPosition[0]=TRThitGlobalPos[0]*
CLHEP::mm;
596 globalPosition[1]=TRThitGlobalPos[1]*
CLHEP::mm;
597 globalPosition[2]=TRThitGlobalPos[2]*
CLHEP::mm;
600 fieldCache.
getField (globalPosition.data(), mField.data());
602 map_x2 = mField.x()*mField.x();
603 map_y2 = mField.y()*mField.y();
604 map_z2 = mField.z()*mField.z();
631 for (;currentClusterIter!=endOfClusterList;++currentClusterIter)
634 const double cluster_x(currentClusterIter->xpos);
635 const double cluster_y(currentClusterIter->ypos);
636 const double cluster_z(this->
setClusterZ(currentClusterIter->zpos, isLong, isShort, isEC));
637 const double cluster_x2(cluster_x*cluster_x);
638 const double cluster_y2(cluster_y*cluster_y);
639 double cluster_r2(cluster_x2+cluster_y2);
642 if (cluster_r2<wire_r2) cluster_r2=wire_r2;
643 if (cluster_r2>straw_r2) cluster_r2=straw_r2;
645 const double cluster_r(std::sqrt(cluster_r2));
646 const double cluster_E(currentClusterIter->energy);
647 const unsigned int nprimaryelectrons(
static_cast<unsigned int>( cluster_E / ionisationPotential + 1.0 ) );
652 double depositEnergy(0.);
656 unsigned int nsurvivingprimaryelectrons =
static_cast<unsigned int>(randBinomial->fire(rndmEngine,nprimaryelectrons) + 0.5);
657 if (nsurvivingprimaryelectrons==0)
continue;
658 const double meanElectronEnergy(ionisationPotential / smearingFactor);
659 for (
unsigned int ielec(0); ielec<nsurvivingprimaryelectrons; ++ielec) {
660 depositEnergy += CLHEP::RandExpZiggurat::shoot(rndmEngine, meanElectronEnergy);
665 const double fluctSigma(sqrt(cluster_E * ionisationPotential * (2 - smearingFactor) / smearingFactor));
667 depositEnergy = CLHEP::RandGaussZiggurat::shoot(rndmEngine, cluster_E, fluctSigma);
668 }
while(depositEnergy<0.0);
683 effectiveField2 = map_z2*cluster_y2/cluster_r2 + map_x2 + map_y2;
692 effectiveField2 = map_z2 + (map_x2+map_y2)*cluster_y2/cluster_r2;
708 double clusterTime(currentClusterIter->time);
718 double timedirect(0.), timereflect(0.);
722 double expdirect(1.0), expreflect(1.0);
736 const unsigned int kdirect =
static_cast<unsigned int>(distdirect/10);
737 const unsigned int kreflect =
static_cast<unsigned int>(distreflect/10);
744 double dt = clusterTime + commondrifttime;
745 deposits.emplace_back(0.5*depositEnergy*expdirect, timedirect+
dt);
746 deposits.emplace_back(0.5*depositEnergy*expreflect, timereflect+
dt);
755 const int mask(0x0000001F);
757 int trtID, ringID, moduleID, layerID, strawID;
758 int wheelID, planeID, sectorID;
763 if ( !(hitID & 0x00200000) ) {
765 strawID = hitID &
mask;
766 hitID >>= word_shift;
767 layerID = hitID &
mask;
768 hitID >>= word_shift;
769 moduleID = hitID &
mask;
770 hitID >>= word_shift;
771 ringID = hitID &
mask;
772 trtID = hitID >> word_shift;
777 const Amg::Vector3D v( (*theHit)->GetPreStepX(),(*theHit)->GetPreStepY(),(*theHit)->GetPreStepZ());
783 strawID = hitID &
mask;
784 hitID >>= word_shift;
785 planeID = hitID &
mask;
786 hitID >>= word_shift;
787 sectorID = hitID &
mask;
788 hitID >>= word_shift;
789 wheelID = hitID &
mask;
790 trtID = hitID >> word_shift;
793 if (trtID == 3) trtID = 0;
798 if ( endcapElement ) {
799 const Amg::Vector3D v( (*theHit)->GetPreStepX(),(*theHit)->GetPreStepY(),(*theHit)->GetPreStepZ());
805 ATH_MSG_WARNING (
"Could not find global coordinate of a straw - drifttime calculation will be inaccurate" );
806 return {0.0,0.0,0.0};
812 double cluster_z(cluster_z_in);
817 const double longBarrelStrawHalfLength(349.315*
CLHEP::mm);
818 const double shortBarrelStrawHalfLength(153.375*
CLHEP::mm);
819 const double EndcapStrawHalfLength(177.150*
CLHEP::mm);
820 if ( isLong && std::abs(cluster_z)>longBarrelStrawHalfLength+30 ) {
821 double d = cluster_z<0 ? cluster_z+longBarrelStrawHalfLength : cluster_z-longBarrelStrawHalfLength;
822 ATH_MSG_WARNING (
"Long barrel straw cluster is outside the active gas volume z = +- 349.315 mm by " <<
d <<
" mm.");
826 if ( isShort && std::abs(cluster_z)>shortBarrelStrawHalfLength+30 ) {
827 double d = cluster_z<0 ? cluster_z+shortBarrelStrawHalfLength : cluster_z-shortBarrelStrawHalfLength;
828 ATH_MSG_WARNING (
"Short barrel straw cluster is outside the active gas volume z = +- 153.375 mm by " <<
d <<
" mm.");
832 if ( isEC && std::abs(cluster_z)>EndcapStrawHalfLength+30 ) {
833 double d = cluster_z<0 ? cluster_z+EndcapStrawHalfLength : cluster_z-EndcapStrawHalfLength;
834 ATH_MSG_WARNING (
"End cap straw cluster is outside the active gas volume z = +- 177.150 mm by " <<
d <<
" mm.");