9 #include "GaudiKernel/IRndmGenSvc.h"
10 #include "GaudiKernel/RndmGenerators.h"
22 #include "CLHEP/Units/SystemOfUnits.h"
23 #include "CLHEP/Matrix/Vector.h"
24 #include "CLHEP/Random/RandFlat.h"
25 #include "CLHEP/Random/RandGauss.h"
26 #include "CLHEP/Geometry/Transform3D.h"
29 #include "GaudiKernel/ITHistSvc.h"
42 m_minimumHadOutEnergy(200.*
CLHEP::
MeV),
45 m_hadIntFromX0(false),
46 m_hadIntProbScale(1.0),
47 m_minimumHadInitialEnergy(1000.*
CLHEP::
MeV),
48 m_particleBroker(
"ISF_ParticleBroker",
n),
49 m_truthRecordSvc(
"ISF_TruthRecordSvc",
n),
50 m_rndGenSvc(
"AtDSFMTGenSvc",
n),
52 m_randomEngine(nullptr),
53 m_randomEngineName(
"FatrasRnd"),
54 m_validationMode(false),
56 m_hadIntValidation(false),
57 m_hadIntValidationTreeName(
"FatrasHadronicInteractions"),
58 m_hadIntValidationTreeDescription(
"Validation output from the HadIntProcessorParametric"),
59 m_hadIntValidationTreeFolder(
"/val/FatrasHadronicInteractions"),
60 m_hadIntValidationTree(nullptr),
66 m_hadIntMotherBarcode(0),
69 m_hadIntMotherPhi(0.),
70 m_hadIntMotherEta(0.),
107 if (m_rndGenSvc.retrieve().isFailure()){
109 return StatusCode::FAILURE;
114 m_randomEngine = m_rndGenSvc->GetEngine(m_randomEngineName);
115 if (!m_randomEngine) {
116 ATH_MSG_FATAL(
"Could not get random engine '" << m_randomEngineName <<
"'" );
117 return StatusCode::FAILURE;
121 if (m_particleBroker.retrieve().isFailure()){
123 return StatusCode::FAILURE;
125 if (m_truthRecordSvc.retrieve().isFailure()){
127 return StatusCode::FAILURE;
131 ATH_CHECK( m_validationTool.retrieve( DisableTool{ m_validationTool.empty() || !m_validationMode } ) );
133 if (m_hadIntValidation){
135 ITHistSvc* tHistSvc =
nullptr;
137 if (service(
"THistSvc",tHistSvc).isFailure())
138 ATH_MSG_ERROR(
"initialize() Could not find Hist Service -> Switching ValidationMode Off !" );
141 ATH_MSG_VERBOSE(
"Booking hadronic interaction validation TTree ... " );
144 m_hadIntValidationTree =
new TTree(m_hadIntValidationTreeName.c_str(), m_hadIntValidationTreeDescription.c_str());
147 m_hadIntValidationTree->Branch(
"HadIntPointX" , &m_hadIntPointX,
"hintX/F");
148 m_hadIntValidationTree->Branch(
"HadIntPointY" , &m_hadIntPointY,
"hintY/F");
149 m_hadIntValidationTree->Branch(
"HadIntPointR" , &m_hadIntPointR,
"hintR/F");
150 m_hadIntValidationTree->Branch(
"HadIntPointZ" , &m_hadIntPointZ,
"hintZ/F");
151 m_hadIntValidationTree->Branch(
"HadIntMotherPdg" , &m_hadIntMotherPdg,
"hintMotherPdg/I");
152 m_hadIntValidationTree->Branch(
"HadIntMotherBarcode" , &m_hadIntMotherBarcode,
"hintMotherBarcode/I");
153 m_hadIntValidationTree->Branch(
"HadIntMotherP" , &m_hadIntMotherP,
"hintMotherP/F");
154 m_hadIntValidationTree->Branch(
"HadIntMotherPt" , &m_hadIntMotherPt,
"hintMotherPt/F");
155 m_hadIntValidationTree->Branch(
"HadIntMotherPhi" , &m_hadIntMotherPhi,
"hintMotherPhi/F");
156 m_hadIntValidationTree->Branch(
"HadIntMotherEta" , &m_hadIntMotherEta,
"hintMohterEta/F");
158 m_hadIntValidationTree->Branch(
"HadIntChildren" , &m_hadIntChildren,
"hintcs/I");
159 m_hadIntValidationTree->Branch(
"HadIntChildE " , &m_hadIntChildE,
"hintce/F");
160 m_hadIntValidationTree->Branch(
"HadIntChildPdg" , m_hadIntChildPdg,
"hintChildPdg[hintcs]/I");
161 m_hadIntValidationTree->Branch(
"HadIntChildP" , m_hadIntChildP,
"hintChildP[hintcs]/F");
162 m_hadIntValidationTree->Branch(
"HadIntChildPcms" , m_hadIntChildPcms,
"hintChildPcms[hintcs]/F");
163 m_hadIntValidationTree->Branch(
"HadIntChildTh" , m_hadIntChildTh,
"hintChildTh[hintcs]/F");
164 m_hadIntValidationTree->Branch(
"HadIntChildThc" , m_hadIntChildThc,
"hintChildThc[hintcs]/F");
165 m_hadIntValidationTree->Branch(
"HadIntChildPhi" , m_hadIntChildPhi,
"hintChildPhi[hintcs]/F");
166 m_hadIntValidationTree->Branch(
"HadIntChildEta" , m_hadIntChildEta,
"hintChildEta[hintcs]/F");
167 m_hadIntValidationTree->Branch(
"HadIntChildDeltaPhi", m_hadIntChildDeltaPhi,
"hintChildPhi[hintcs]/F");
168 m_hadIntValidationTree->Branch(
"HadIntChildDeltaEta", m_hadIntChildDeltaEta,
"hintChildEta[hintcs]/F");
171 if ((tHistSvc->regTree(m_hadIntValidationTreeFolder, m_hadIntValidationTree)).isFailure()) {
172 ATH_MSG_ERROR(
"initialize() Could not register the validation Tree -> Switching ValidationMode Off !" );
173 delete m_hadIntValidationTree; m_hadIntValidationTree =
nullptr;
175 ATH_MSG_INFO(
"TTree for Hadronic Interactions validation booked." );
179 return StatusCode::SUCCESS;
186 return StatusCode::SUCCESS;
191 double p,
double ,
double charge,
199 if (extMprop && !m_hadIntFromX0) {
207 ATH_MSG_VERBOSE(
" [ had ] Nuclear Interaction length not available, using appox. from X0" );
214 << (1. -
prob) * 100. * m_hadIntProbScale <<
" %." );
219 if (CLHEP::RandFlat::shoot(m_randomEngine) < (1. -
prob) * m_hadIntProbScale )
return recordHadState(0.,
p,
232 double al =
mat->l0();
243 al *= 1./(1.+
exp(-0.5*(
p-270.)*(
p-270.)/60./60.));
262 double E = sqrt(
p*
p +
m*
m);
265 double multiplicity_max = 0.25 *
E/1000. + 18.;
268 double randx , randy,
arg = 0.;
282 randx = 30. * CLHEP::RandFlat::shoot(m_randomEngine);
283 randy = 1. * CLHEP::RandFlat::shoot(m_randomEngine);
284 arg =
exp(-0.5*( (randx-p1)/p2 +
exp(-(randx-p1)/p2) ) );
285 if (randy < arg && randx>3 && randx<multiplicity_max)
break;
288 randx *= (1.2-0.4*
exp(-0.001*
p));
290 int Npart = (
int)randx;
294 ATH_MSG_VERBOSE(
"[ had ] Number of particles smaller than 3, parameterisation not valid."
295 <<
" Doing Nothing");
300 <<
" with " << Npart <<
" outgoing particles " );
305 ATH_MSG_DEBUG(
"[ had ] create hadronic shower for particle with PDG ID "
306 <<
parent->pdgCode() <<
" and barcode "
313 if (m_hadIntValidationTree){
314 m_hadIntPointX =
vertex.x();
315 m_hadIntPointY =
vertex.y();
316 m_hadIntPointZ =
vertex.z();
317 m_hadIntPointR =
vertex.perp();
318 m_hadIntMotherPdg =
parent->pdgCode();
321 m_hadIntMotherPt =
p*particleDir.perp();
322 m_hadIntMotherPhi = particleDir.phi();
323 m_hadIntMotherEta = particleDir.eta();
325 m_hadIntChildren = 0;
329 <<
E <<
" | " <<
m <<
" | " <<
p <<
" | " );
332 if (m_hadIntValidationTree) m_hadIntValidationTree->Fill();
333 ATH_MSG_VERBOSE(
"[ had ] interaction initiated by a secondary particle, no children saved " );
341 double chargedist = 0.;
342 std::vector<double>
charge(Npart);
343 std::vector<Trk::ParticleHypothesis> childType(Npart);
344 std::vector<double> newm(Npart);
345 std::vector<int> pdgid(Npart);
373 for (
int i=0;
i<Npart;
i++) {
374 chargedist = CLHEP::RandFlat::shoot(m_randomEngine);
375 if (chargedist<pif) {
382 if ( chargedist<2*pif) {
389 if (chargedist<3*pif) {
396 if (chargedist<3*pif+nef) {
403 if (chargedist<3*pif+nef+prf) {
418 for (
int i=1;
i<Npart;
i++) {
420 childType[
i]=childType[0];
451 std::vector<double>
mom(Npart);
452 std::vector<double>
th(Npart);
453 std::vector<double> ph(Npart);
456 double eps = 2./Npart;
457 double rnd = CLHEP::RandFlat::shoot(m_randomEngine);
458 mom[0] = 0.5*
pow(eps,rnd);
459 th[0] = acos( 2*CLHEP::RandFlat::shoot(m_randomEngine)-1.);
460 ph[0] = 2*
M_PI*CLHEP::RandFlat::shoot(m_randomEngine);
466 double ptot =
mom[0];
469 for (
int i=1;
i<Npart-2;
i++) {
472 mom[
i] = ( eps + CLHEP::RandFlat::shoot(m_randomEngine)*(1-eps))*(1-ptot);
473 if (ptemp.mag()<1-ptot) {
474 while ( fabs(ptemp.mag()-
mom[
i])>1-ptot-
mom[
i] ){
475 mom[
i] = ( eps + CLHEP::RandFlat::shoot(m_randomEngine)*(1-eps))*(1-ptot);
479 double p_rem=1-ptot-
mom[
i];
480 double cthmax = fmin(1.,(-ptemp.mag()*ptemp.mag()-
mom[
i]*
mom[
i]+p_rem*p_rem)/2/ptemp.mag()/
mom[
i]);
482 double rnd = CLHEP::RandFlat::shoot(m_randomEngine);
483 theta = acos( (cthmax+1.)*rnd-1.);
484 phi = 2*
M_PI*CLHEP::RandFlat::shoot(m_randomEngine);
486 HepGeom::Vector3D<double> dnewHep = HepGeom::RotateZ3D(ptemp.phi())*HepGeom::RotateY3D(ptemp.theta())*
test;
490 ptemp +=
mom[
i]*dnew;
495 mom[Npart-2] =
pow(eps,CLHEP::RandFlat::shoot(m_randomEngine))*(1-ptot);
496 mom[Npart-1] = 1-ptot-
mom[Npart-2];
498 if (ptemp.mag()<1-ptot) {
499 while (
mom[Npart-1]+
mom[Npart-2]<ptemp.mag()) {
500 mom[Npart-2] =
pow(eps,CLHEP::RandFlat::shoot(m_randomEngine))*(1-ptot);
501 mom[Npart-1] = 1-ptot-
mom[Npart-2];
504 if (ptemp.mag()<fabs(
mom[Npart-1]-
mom[Npart-2]) ) {
505 double diff = ptemp.mag()*CLHEP::RandFlat::shoot(m_randomEngine);
510 double cth =(-ptemp.mag()*ptemp.mag()-
mom[Npart-2]*
mom[Npart-2]+
mom[Npart-1]*
mom[Npart-1])/2/ptemp.mag()/
mom[Npart-2];
511 if (fabs(cth)>1.) cth = (cth>0.) ? 1. : -1.;
514 phi = 2*
M_PI*CLHEP::RandFlat::shoot(m_randomEngine);
516 HepGeom::Vector3D<double> dnewHep = HepGeom::RotateZ3D(ptemp.phi())*HepGeom::RotateY3D(ptemp.theta())*
test;
519 th[Npart-2]=dnew.theta();
520 ph[Npart-2]=dnew.phi();
521 ptemp +=
mom[Npart-2]*dnew;
523 th[Npart-1]=dlast.theta();
524 ph[Npart-1]=dlast.phi();
529 for (
int i=0;
i<Npart;
i++) etot += sqrt(
mom[
i]*
mom[
i]+newm[
i]*newm[
i]);
531 for (
int i=0;
i<Npart;
i++) summ += newm[
i];
535 float scale = sqrt(summ*summ+2*summ*
p+
m*
m)/etot;
537 for (
int i=0;
i<Npart;
i++) {
543 CLHEP::HepLorentzVector bv(
p*particleDir.unit().x(),
p*particleDir.unit().y(),
p*particleDir.unit().z(),sqrt(etot*etot+
p*
p));
544 std::vector<CLHEP::HepLorentzVector> childBoost(Npart);
552 for (
int i=0;
i<Npart;
i++) {
557 CLHEP::HepLorentzVector newp(childP.x(),childP.y(),childP.z(),sqrt(
mom[
i]*
mom[
i]+newm[
i]*newm[
i]));
558 CLHEP::HepLorentzVector childPB = newp.boost(bv.boostVector());
559 childBoost[
i]=childPB;
571 unsigned short numChildren = 0;
573 for (
int i=0;
i<Npart;
i++) {
574 if (pdgid[
i]<10000) {
578 eout += childBoost[
i].e();
582 m_hadIntChildPdg[m_hadIntChildren] = pdgid[
i];
583 m_hadIntChildP[m_hadIntChildren] = childP.mag();
584 m_hadIntChildPcms[m_hadIntChildren] =
mom[
i];
585 m_hadIntChildTh[m_hadIntChildren] = childP.unit().dot(particleDir);
586 m_hadIntChildThc[m_hadIntChildren] =chP.dot(particleDir);
587 m_hadIntChildPhi[m_hadIntChildren] = childP.phi();
588 m_hadIntChildEta[m_hadIntChildren] = childP.eta();
589 double deltaPhi = m_hadIntMotherPhi - m_hadIntChildPhi[m_hadIntChildren];
593 m_hadIntChildDeltaPhi[m_hadIntChildren] =
deltaPhi;
594 m_hadIntChildDeltaEta[m_hadIntChildren] = m_hadIntMotherEta - m_hadIntChildEta[m_hadIntChildren];
598 if (childP.mag()> m_minimumHadOutEnergy) {
616 if (m_validationMode) {
625 ++childrenIt; numChildren++;
638 m_truthRecordSvc->registerTruthIncident( truth);
644 if (m_validationMode && m_validationTool.isEnabled()) {
651 m_hadIntChildE = eout;
653 if (m_hadIntValidationTree) m_hadIntValidationTree->Fill();
655 ATH_MSG_VERBOSE(
"[ had ] it was kinematically possible to create " << gen_part <<
" shower particles " );
677 if (processSecondaries && !ispVec.empty() ) {
678 for (
auto *childParticle : ispVec) {
680 if (!childParticle->getTruthBinding()) {
681 ATH_MSG_ERROR(
"Could not retrieve TruthBinding from child ISFParticle "<< *childParticle);
683 m_particleBroker->push(childParticle,
parent);
717 if (!ispVec.empty() ) {
718 for (
auto *childParticle : ispVec) {
720 if (!childParticle->getTruthBinding()) {
721 ATH_MSG_ERROR(
"Could not retrieve TruthBinding from child ISFParticle "<< *childParticle);
723 m_particleBroker->push(childParticle,
parent);