ATLAS Offline Software
FullModelHadronicProcess.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include "FullModelHadronicProcess.hh"
6 #include "G4Version.hh"
7 #include "G4ProcessManager.hh"
8 #include "G4ProcessHelper.hh"
9 #include "G4ParticleTable.hh"
10 #include "FullModelReactionDynamics.hh"
11 
12 #if G4VERSION_NUMBER < 1100
13 #include "G4HadReentrentException.hh"
14 #else
15 #include "G4HadronicException.hh"
16 #endif
17 
18 #include "CustomParticle.h"
19 
20 
21 FullModelHadronicProcess::FullModelHadronicProcess(const G4String& processName) :
22  G4VDiscreteProcess(processName,fUserDefined),m_theParticle(0),m_newParticle(0),m_cache(0)
23 {
24  // Instantiating helper class
25  m_theHelper = G4ProcessHelper::Instance();
26 
27 }
28 
29 
30 G4bool FullModelHadronicProcess::IsApplicable(const G4ParticleDefinition& aP)
31 {
32 
33  return m_theHelper->ApplicabilityTester(aP);
34 
35 }
36 
37 G4double FullModelHadronicProcess::GetMicroscopicCrossSection(const G4DynamicParticle *aParticle,
38  const G4Element *anElement,
39  G4double /*aTemp*/)
40 {
41  // Get the cross section for this particle/element combination from the ProcessHelper
42  G4double InclXsec = m_theHelper->GetInclusiveCrossSection(aParticle,anElement);
43  // G4cout<<"Returned cross section from helper was: "<<InclXsec/CLHEP::millibarn<<" millibarn"<<G4endl;
44 
45  // Need to provide Set-methods for these in time
46  G4double HighestEnergyLimit = 10 * CLHEP::TeV ;
47  G4double LowestEnergyLimit = 1 * CLHEP::eV;
48 
49  G4double ParticleEnergy = aParticle->GetKineticEnergy();
50 
51 
52  if (ParticleEnergy > HighestEnergyLimit || ParticleEnergy < LowestEnergyLimit){
53  return 0;
54  } else {
55  // G4cout << "Microscopic Cross Section: "<<InclXsec / CLHEP::millibarn<<" millibarn"<<G4endl;
56  return InclXsec;
57  }
58 
59 }
60 
61 G4double FullModelHadronicProcess::GetMeanFreePath(const G4Track& aTrack, G4double, G4ForceCondition*)
62 {
63  G4Material *aMaterial = aTrack.GetMaterial();
64  const G4DynamicParticle *aParticle = aTrack.GetDynamicParticle();
65  G4double sigma = 0.0;
66 
67  G4int nElements = aMaterial->GetNumberOfElements();
68 
69  const G4double *theAtomicNumDensityVector =
70  aMaterial->GetAtomicNumDensityVector();
71  G4double aTemp = aMaterial->GetTemperature();
72 
73  for( G4int i=0; i<nElements; ++i )
74  {
75  G4double xSection =
76  GetMicroscopicCrossSection( aParticle, (*aMaterial->GetElementVector())[i], aTemp);
77  sigma += theAtomicNumDensityVector[i] * xSection;
78  }
79 
80  return 1.0/sigma;
81 
82 }
83 
84 
85 G4VParticleChange* FullModelHadronicProcess::PostStepDoIt(const G4Track& aTrack,
86  const G4Step& /*aStep*/)
87 {
88  // A little setting up
89  aParticleChange.Initialize(aTrack);
90  // G4DynamicParticle* OrgPart = const_cast<G4DynamicParticle*>(aTrack.GetDynamicParticle());
91  const G4DynamicParticle* IncidentRhadron = aTrack.GetDynamicParticle();
92  CustomParticle* CustomIncident = static_cast<CustomParticle*>(IncidentRhadron->GetDefinition());
93  const G4ThreeVector& aPosition = aTrack.GetPosition();
94  const G4int theIncidentPDG = IncidentRhadron->GetDefinition()->GetPDGEncoding();
95  G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();
96  std::vector<G4ParticleDefinition*> theParticleDefinitions;
97  G4bool IncidentSurvives = false;
98  G4bool TargetSurvives = false;
99  G4Nucleus targetNucleus(aTrack.GetMaterial());
100  G4ParticleDefinition* outgoingRhadron=0;
101  G4ParticleDefinition* outgoingCloud=0;
102  G4ParticleDefinition* outgoingTarget=0;
103 
104 
105  G4ThreeVector p_0 = IncidentRhadron->GetMomentum();
106 
107  // static int n=0;
108 
109  G4double e_kin_0 = IncidentRhadron->GetKineticEnergy();
110 
111  G4DynamicParticle* cloudParticle = new G4DynamicParticle();
112 
113  cloudParticle->SetDefinition(CustomIncident->GetCloud());
114 
115  double scale=cloudParticle->GetDefinition()->GetPDGMass()/IncidentRhadron->GetDefinition()->GetPDGMass();
116 
117  G4LorentzVector cloudMomentum;
118  cloudMomentum.setVectM(IncidentRhadron->GetMomentum()*scale,cloudParticle->GetDefinition()->GetPDGMass());
119  G4LorentzVector gluinoMomentum;
120 
121  gluinoMomentum.setVectM(IncidentRhadron->GetMomentum()*(1.-scale),CustomIncident->GetSpectator()->GetPDGMass());
122 
123  G4LorentzVector FullRhadron4Momentum = IncidentRhadron->Get4Momentum();
124  G4LorentzVector Cloud4Momentum = cloudMomentum;
125 
126  cloudParticle->Set4Momentum(cloudMomentum);
127 
128  G4DynamicParticle* OrgPart = cloudParticle;
129 
130  G4double ek = OrgPart->GetKineticEnergy();
131  G4double amas = OrgPart->GetDefinition()->GetPDGMass();
132  G4double tkin = targetNucleus.Cinema( ek );
133  ek += tkin;
134  OrgPart->SetKineticEnergy( ek );
135  G4double et = ek + amas;
136  G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
137  G4double pp = OrgPart->GetMomentum().mag();
138  if( pp > 0.0 )
139  {
140  G4ThreeVector momentum = OrgPart->GetMomentum();
141  OrgPart->SetMomentum( momentum * (p/pp) );
142  }
143 
144  // calculate black track energies
145  if(ek>0.) {tkin = targetNucleus.EvaporationEffects( ek ); ek -= tkin;}
146  if(ek+gluinoMomentum.e()-gluinoMomentum.m()<=0.1*CLHEP::MeV||ek<=0.) {
147  G4cout<<"Stopping particle...:"<<G4endl;
148  G4cout<<"ek: "<<ek/CLHEP::MeV<<" MeV"<<G4endl;
149  aParticleChange.ProposeTrackStatus( fStopButAlive );
150  // h4->Fill(aPosition.x()/CLHEP::cm);
151  return &aParticleChange;
152  }
153  OrgPart->SetKineticEnergy( ek );
154  et = ek + amas;
155  p = std::sqrt( std::abs((et-amas)*(et+amas)) );
156  pp = OrgPart->GetMomentum().mag();
157 
158  if( pp > 0.0 )
159  {
160  G4ThreeVector momentum = OrgPart->GetMomentum();
161  OrgPart->SetMomentum( momentum * (p/pp) );
162  }
163 
164 
165 
166  //Get the final state particles
167 
168  G4ParticleDefinition* aTarget;
169  ReactionProduct rp = m_theHelper->GetFinalState(aTrack,aTarget);
170 
171  //FIXME: Due to force2to2 always being set to false the while loop
172  //below can never be called, so have commented it out.
173 
174  // G4bool force2to2 = false;
175  // // G4cout<<"Trying to get final state..."<<G4endl;
176  // while(rp.size()!=2 && force2to2){
177  // rp = m_theHelper->GetFinalState(aTrack,aTarget);
178  // }
179 
180  G4int NumberOfSecondaries = rp.size();
181  // G4cout<<"Multiplicity of selected final state: "<<rp.size()<<G4endl;
182 
183  //Getting CMS transforms. Boosting is done at histogram filling
184  G4LorentzVector Target4Momentum;
185  Target4Momentum.setVectM(CLHEP::Hep3Vector(0.,0.,0.),aTarget->GetPDGMass());
186  // Target4Momentum.setVectM(0.,targetNucleus.GetN()*CLHEP::GeV);
187  G4LorentzVector psum_full,psum_cloud;
188  psum_full = FullRhadron4Momentum + Target4Momentum;
189  psum_cloud = Cloud4Momentum + Target4Momentum;
190  G4ThreeVector trafo_full_cms = (-1)*psum_full.boostVector();
191  G4ThreeVector trafo_cloud_cms = (-1)*psum_cloud.boostVector();
192 
193  // OK Let's make some particles :-)
194  // We're not using them for anything yet, I know, but let's make sure the machinery is there
195 
196  for(ReactionProduct::iterator it = rp.begin();
197  it != rp.end();
198  ++it)
199  {
200  G4ParticleDefinition* tempDef = theParticleTable->FindParticle(*it);
201  CustomParticle* tempCust = dynamic_cast<CustomParticle*>(tempDef);
202  if (tempDef==aTarget) TargetSurvives = true;
203 
204  if (tempCust!=0)
205  {
206  outgoingRhadron = tempDef;
207  outgoingCloud=tempCust->GetCloud();
208  if(outgoingCloud == 0)
209  {
210  G4cout << "FullModelHadronicProcess::PostStepDoIt() Definition of outgoing particle cloud for " << tempDef->GetParticleName() << " not available!!" << G4endl;
211  }
212  }
213 
214  if (tempDef==G4Proton::Proton()||tempDef==G4Neutron::Neutron()) outgoingTarget = tempDef;
215  if (tempCust==0&&rp.size()==2) outgoingTarget = tempDef;
216  if (tempDef->GetPDGEncoding()==theIncidentPDG) {
217  IncidentSurvives = true;
218  } else {
219  theParticleDefinitions.push_back(tempDef);
220  }
221  }
222 
223  if (NULL==outgoingRhadron) {
224  //FIXME: What should be done in the case that outgoingRhadron is
225  //still 0 at this point? Throwing an exception for now to prevent
226  //a null-pointer being dereferenced later on in the code.
227  G4ExceptionDescription description;
228  description << "PostStepDoIt: outgoingRHadron is null";
229  G4Exception("FullModelHadronicProcess", "NoOutGoingRHadron", FatalException, description);
230  }
231 
232  if (outgoingTarget==0) outgoingTarget = theParticleTable->FindParticle(rp[1]);
233 
234  // If the incident particle survives it is not a "secondary", although
235  // it would probably be easier to fStopAndKill the track every time.
236  if(IncidentSurvives) NumberOfSecondaries--;
237  aParticleChange.SetNumberOfSecondaries(NumberOfSecondaries);
238 
239 
240  // ADAPTED FROM G4LEPionMinusInelastic::ApplyYourself
241  // It is bloody ugly, but such is the way of cut 'n' paste
242 
243  // Set up the incident
244  const G4HadProjectile* originalIncident = new G4HadProjectile(*OrgPart);//This is where rotation to z-axis is done (Aarrggh!)
245 
246 
247  //Maybe I need to calculate trafo from R-hadron... Bollocks! Labframe does not move. CMS does.
248  G4HadProjectile* org2 = new G4HadProjectile(*OrgPart);
249  G4LorentzRotation trans = org2->GetTrafoToLab();
250  delete org2;
251 
252  G4DynamicParticle *originalTarget = new G4DynamicParticle;
253  originalTarget->SetDefinition(aTarget);
254 
255  G4ReactionProduct targetParticle(aTarget);
256 
257 
258  G4ReactionProduct currentParticle(const_cast<G4ParticleDefinition *>(originalIncident->GetDefinition() ) );
259  currentParticle.SetMomentum( originalIncident->Get4Momentum().vect() );
260  currentParticle.SetKineticEnergy( originalIncident->GetKineticEnergy() );
261 
262 
263  G4ReactionProduct modifiedOriginal = currentParticle;
264 
265  currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
266  targetParticle.SetSide( -1 ); // target always goes in backward hemisphere
267  G4bool incidentHasChanged = false;
268  if (!IncidentSurvives) incidentHasChanged = true; //I wonder if I am supposed to do this...
269  G4bool targetHasChanged = false;
270  if (!TargetSurvives) targetHasChanged = true; //Ditto here
271  G4bool quasiElastic = false;
272  if (rp.size()==2) quasiElastic = true; //Oh well...
273  G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec; // vec will contain the secondary particles
274  G4int vecLen = 0;
275  vec.Initialize( 0 );
276 
277 
278 
279  // I hope my understanding of "secondary" is correct here
280  // I think that it entails only what is not a surviving incident of target
281 
282  for (G4int i = 0; i!=NumberOfSecondaries;i++){
283  if(theParticleDefinitions[i]!=aTarget
284  && theParticleDefinitions[i]!=originalIncident->GetDefinition()
285  && theParticleDefinitions[i]!=outgoingRhadron
286  && theParticleDefinitions[i]!=outgoingTarget)
287  {
288  G4ReactionProduct* pa = new G4ReactionProduct;
289  pa->SetDefinition( theParticleDefinitions[i] );
290  (G4UniformRand() < 0.5) ? pa->SetSide( -1 ) : pa->SetSide( 1 );
291  vec.SetElement( vecLen++, pa );
292  }
293  }
294 
295  // if (incidentHasChanged) currentParticle.SetDefinitionAndUpdateE( outgoingRhadron );
296 
297  if (incidentHasChanged) currentParticle.SetDefinitionAndUpdateE( outgoingCloud );
298  if (incidentHasChanged) modifiedOriginal.SetDefinition( outgoingCloud );//Is this correct? It solves the "free energy" checking in ReactionDynamics
299  if (targetHasChanged) targetParticle.SetDefinitionAndUpdateE( outgoingTarget );
300 
301 
302  CalculateMomenta( vec, vecLen,
303  originalIncident, originalTarget, modifiedOriginal,
304  targetNucleus, currentParticle, targetParticle,
305  incidentHasChanged, targetHasChanged, quasiElastic );
306 
307 
308  G4String cPname = currentParticle.GetDefinition()->GetParticleName();
309 
310 
311  aParticleChange.SetNumberOfSecondaries(vecLen+NumberOfSecondaries);
312  G4double e_kin;
313  G4LorentzVector cloud_p4_new; //Cloud 4-momentum in lab after collision
314 
315  cloud_p4_new.setVectM(currentParticle.GetMomentum(),currentParticle.GetMass());
316  cloud_p4_new *= trans;
317 
318  G4LorentzVector cloud_p4_old_full = Cloud4Momentum; //quark system in CMS BEFORE collision
319  cloud_p4_old_full.boost(trafo_full_cms);
320  G4LorentzVector cloud_p4_old_cloud = Cloud4Momentum; //quark system in cloud CMS BEFORE collision
321  cloud_p4_old_cloud.boost(trafo_cloud_cms);
322  G4LorentzVector cloud_p4_full = cloud_p4_new; //quark system in CMS AFTER collision
323  cloud_p4_full.boost(trafo_full_cms);
324  G4LorentzVector cloud_p4_cloud = cloud_p4_new; //quark system in cloud CMS AFTER collision
325  cloud_p4_cloud.boost(trafo_cloud_cms);
326 
327  G4LorentzVector p_g_cms = gluinoMomentum; //gluino in CMS BEFORE collision
328  p_g_cms.boost(trafo_full_cms);
329 
330  G4LorentzVector p4_new;
331  p4_new.setVectM( cloud_p4_new.v() + gluinoMomentum.v(), outgoingRhadron->GetPDGMass() );
332 
333  G4ThreeVector p_new;
334  p_new = p4_new.vect();
335 
336  aParticleChange.ProposeLocalEnergyDeposit((p4_new-cloud_p4_new-gluinoMomentum).m());
337 
338  if( incidentHasChanged )
339  {
340  G4DynamicParticle* p0 = new G4DynamicParticle;
341  p0->SetDefinition( outgoingRhadron );
342  p0->SetMomentum( p_new );
343 
344  // May need to run SetDefinitionAndUpdateE here...
345  G4Track* Track0 = new G4Track(p0,
346  aTrack.GetGlobalTime(),
347  aPosition);
348  aParticleChange.AddSecondary(Track0);
349 
350  // Because of the above calculations, and the use of mass, there's going to be a lot of squaring and
351  // square rooting to get the total energy. For that reason, we should allow a little buffer before
352  // we freak out over energy non-conservation...
353  if(p0->GetKineticEnergy()>e_kin_0+1.e-9) { // Allow 1 meV of energy non-conservation in an interaction
354  G4cout<<"ALAAAAARM!!! (incident changed from "
355  <<IncidentRhadron->GetDefinition()->GetParticleName()
356  <<" to "<<p0->GetDefinition()->GetParticleName()<<")"<<G4endl;
357  G4cout<<"Energy loss: "<<(e_kin_0-p0->GetKineticEnergy())/CLHEP::GeV<<" GeV (should be positive)"<<G4endl;
358  //Turns out problem is only in 2 -> 3 (Won't fix 2 -> 2 energy deposition)
359  if(rp.size()!=3) G4cout<<"DOUBLE ALAAAAARM!!!"<<G4endl;
360  } /*else {
361  G4cout<<"NO ALAAAAARM!!!"<<G4endl;
362  }*/
363  if(std::abs(p0->GetKineticEnergy()-e_kin_0)>100*CLHEP::GeV) {
364  G4cout<<"Diff. too big"<<G4endl;
365  }
366 
367  aParticleChange.ProposeTrackStatus( fStopAndKill );
368  }
369  else
370  {
371 
372  G4double p = p_new.mag();
373  if( p > DBL_MIN )
374  aParticleChange.ProposeMomentumDirection( p_new.x()/p, p_new.y()/p, p_new.z()/p );
375  else
376  aParticleChange.ProposeMomentumDirection( 1.0, 0.0, 0.0 );
377 
378  G4double aE = sqrt(p*p+(outgoingRhadron->GetPDGMass()*outgoingRhadron->GetPDGMass()) );
379  e_kin = aE - outgoingRhadron->GetPDGMass();
380 
381  if(e_kin>e_kin_0) {
382  G4cout<<"ALAAAAARM!!!"<<G4endl;
383  G4cout<<"Energy loss: "<<(e_kin_0-e_kin)/CLHEP::GeV<<" GeV (should be positive)"<<G4endl;
384  if(rp.size()!=3) G4cout<<"DOUBLE ALAAAAARM!!!"<<G4endl;
385  }
386  if(std::abs(e_kin-e_kin_0)>100*CLHEP::GeV) {
387  G4cout<<"Diff. too big"<<G4endl;
388  }
389 
390  if (std::fabs(aE)<.1*CLHEP::eV) aE=.1*CLHEP::eV;
391  aParticleChange.ProposeEnergy( aE- outgoingRhadron->GetPDGMass() ); //I do not know if this is correct...
392  if(std::abs(e_kin-e_kin_0)>100*CLHEP::GeV) {
393  G4cout<<"Diff. too big"<<G4endl;
394  }
395 
396  }
397 
398  if( targetParticle.GetMass() > 0.0 ) // targetParticle can be eliminated in TwoBody
399  {
400  G4DynamicParticle *p1 = new G4DynamicParticle;
401  p1->SetDefinition( targetParticle.GetDefinition() );
402  // G4cout<<"Target secondary: "<<targetParticle.GetDefinition()->GetParticleName()<<G4endl;
403  G4ThreeVector momentum = targetParticle.GetMomentum();
404  momentum = momentum.rotate(m_cache,m_what);
405  p1->SetMomentum( momentum );
406  p1->SetMomentum( (trans*p1->Get4Momentum()).vect());
407  G4Track* Track1 = new G4Track(p1,
408  aTrack.GetGlobalTime(),
409  aPosition);
410  aParticleChange.AddSecondary(Track1);
411  }
412  G4DynamicParticle *pa;
413 
414 
415  for(int i=0; i<vecLen; ++i )
416  {
417  pa = new G4DynamicParticle();
418  pa->SetDefinition( vec[i]->GetDefinition() );
419  pa->SetMomentum( vec[i]->GetMomentum() );
420  pa->Set4Momentum(trans*(pa->Get4Momentum()));
421  G4ThreeVector pvec = pa->GetMomentum();
422  G4Track* Trackn = new G4Track(pa,
423  aTrack.GetGlobalTime(),
424  aPosition);
425  aParticleChange.AddSecondary(Trackn);
426  delete vec[i];
427  }
428 
429  delete originalIncident;
430  delete originalTarget;
431 
432  //clear interaction length
433  ClearNumberOfInteractionLengthLeft();
434 
435 
436  return &aParticleChange;
437 
438 }
439 
440 
441 void FullModelHadronicProcess::CalculateMomenta(
442  G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
443  G4int &vecLen,
444  const G4HadProjectile *originalIncident, // the original incident particle
445  const G4DynamicParticle *originalTarget,
446  G4ReactionProduct &modifiedOriginal, // Fermi motion and evap. effects included
447  G4Nucleus &targetNucleus,
448  G4ReactionProduct &currentParticle,
449  G4ReactionProduct &targetParticle,
450  G4bool &incidentHasChanged,
451  G4bool &targetHasChanged,
452  G4bool quasiElastic )
453 {
454  FullModelReactionDynamics theReactionDynamics;
455 
456  m_cache = 0;
457  m_what = originalIncident->Get4Momentum().vect();
458 
459  //Commented out like in casqmesmin.F where CALL STPAIR is commented out
460  /*
461  theReactionDynamics.ProduceStrangeParticlePairs( vec, vecLen,
462  modifiedOriginal, originalTarget,
463  currentParticle, targetParticle,
464  incidentHasChanged, targetHasChanged );
465  */
466 
467  if( quasiElastic )
468  {
469  // G4cout<<"We are calling TwoBody..."<<G4endl;
470  theReactionDynamics.TwoBody( vec, vecLen,
471  modifiedOriginal, originalTarget,
472  currentParticle, targetParticle,
473  targetNucleus, targetHasChanged );
474 
475  return;
476  }
477 
478  //If ProduceStrangeParticlePairs is commented out, let's cut this one as well
479  G4ReactionProduct leadingStrangeParticle;
480  G4bool leadFlag = MarkLeadingStrangeParticle( currentParticle,
481  targetParticle,
482  leadingStrangeParticle );
483 
484 
485 
486  //
487  // Note: the number of secondaries can be reduced in GenerateXandPt and TwoCluster
488  //
489  G4bool finishedGenXPt = false;
490  G4bool annihilation = false;
491  if( originalIncident->GetDefinition()->GetPDGEncoding() < 0 &&
492  currentParticle.GetMass() == 0.0 && targetParticle.GetMass() == 0.0 )
493  {
494  // original was an anti-particle and annihilation has taken place
495  annihilation = true;
496  G4double ekcor = 1.0;
497  G4double ek = originalIncident->GetKineticEnergy();
498  G4double ekOrg = ek;
499 
500  const G4double tarmas = originalTarget->GetDefinition()->GetPDGMass();
501  if( ek > 1.0*CLHEP::GeV )ekcor = 1./(ek/CLHEP::GeV);
502  const G4double atomicWeight = targetNucleus.AtomicMass(targetNucleus.GetA_asInt(),targetNucleus.GetZ_asInt()); //targetNucleus.GetN();
503  ek = 2*tarmas + ek*(1.+ekcor/atomicWeight);
504 
505  G4double tkin = targetNucleus.Cinema( ek );
506  ek += tkin;
507  ekOrg += tkin;
508  modifiedOriginal.SetKineticEnergy( ekOrg );
509  }
510 
511  const G4double twsup[] = { 1.0, 0.7, 0.5, 0.3, 0.2, 0.1 };
512  G4double rand1 = G4UniformRand();
513  G4double rand2 = G4UniformRand();
514  if (annihilation || vecLen >= 6 ||
515  (modifiedOriginal.GetKineticEnergy()/CLHEP::GeV >= 1.0 &&
516  ((rand1 < 0.5 &&
517  (originalIncident->GetDefinition() == G4KaonPlus::KaonPlus() ||
518  originalIncident->GetDefinition() == G4KaonMinus::KaonMinus() ||
519  originalIncident->GetDefinition() == G4KaonZeroLong::KaonZeroLong() ||
520  originalIncident->GetDefinition() == G4KaonZeroShort::KaonZeroShort()))
521  || rand2 > twsup[vecLen])))
522  finishedGenXPt =
523  theReactionDynamics.GenerateXandPt(vec, vecLen,
524  modifiedOriginal, originalIncident,
525  currentParticle, targetParticle,
526  targetNucleus, incidentHasChanged,
527  targetHasChanged, leadFlag,
528  leadingStrangeParticle);
529  if( finishedGenXPt )
530  {
531  Rotate(vec, vecLen);
532  return;
533  }
534 
535  G4bool finishedTwoClu = false;
536  if( modifiedOriginal.GetTotalMomentum()/CLHEP::MeV < 1.0 )
537  {
538 
539  for(G4int i=0; i<vecLen; i++) delete vec[i];
540  vecLen = 0;
541  }
542  else
543  {
544 
545  theReactionDynamics.SuppressChargedPions( vec, vecLen,
546  modifiedOriginal, currentParticle,
547  targetNucleus,
548  incidentHasChanged );
549 
550  try
551  {
552  finishedTwoClu = theReactionDynamics.TwoCluster( vec, vecLen,
553  modifiedOriginal, originalIncident,
554  currentParticle, targetParticle,
555  targetNucleus, incidentHasChanged,
556  targetHasChanged, leadFlag,
557  leadingStrangeParticle );
558  }
559 #if G4VERSION_NUMBER < 1100
560  catch(G4HadReentrentException& aC)
561  {
562  aC.Report(G4cout);
563  throw G4HadReentrentException(__FILE__, __LINE__, "Failing to calculate momenta");
564  }
565 #else
566  catch(G4HadronicException& aC)
567  {
568  aC.Report(G4cout);
569  throw G4HadronicException(__FILE__, __LINE__, "Failing to calculate momenta");
570  }
571 #endif
572  }
573  if( finishedTwoClu )
574  {
575  Rotate(vec, vecLen);
576  return;
577  }
578 
579  //
580  // PNBlackTrackEnergy is the kinetic energy available for
581  // proton/neutron black track particles [was enp(1) in fortran code]
582  // DTABlackTrackEnergy is the kinetic energy available for
583  // deuteron/triton/alpha particles [was enp(3) in fortran code]
584  //const G4double pnCutOff = 0.1;
585  //const G4double dtaCutOff = 0.1;
586  //if( (targetNucleus.GetN() >= 1.5)
587  // && !(incidentHasChanged || targetHasChanged)
588  // && (targetNucleus.GetPNBlackTrackEnergy()/CLHEP::MeV <= pnCutOff)
589  // && (targetNucleus.GetDTABlackTrackEnergy()/CLHEP::MeV <= dtaCutOff) )
590  //{
591  // the atomic weight of the target nucleus is >= 1.5 AND
592  // neither the incident nor the target particles have changed AND
593  // there is no kinetic energy available for either proton/neutron
594  // or for deuteron/triton/alpha black track particles
595  // For diffraction scattering on heavy nuclei use elastic routines instead
596  //G4cerr << "*** Error in G4InelasticInteraction::CalculateMomenta" << G4endl;
597  //G4cerr << "*** the elastic scattering would be better here ***" <<G4endl;
598  //}
599  theReactionDynamics.TwoBody( vec, vecLen,
600  modifiedOriginal, originalTarget,
601  currentParticle, targetParticle,
602  targetNucleus, targetHasChanged );
603 }
604 
605 
606 G4bool FullModelHadronicProcess::MarkLeadingStrangeParticle(
607  const G4ReactionProduct &currentParticle,
608  const G4ReactionProduct &targetParticle,
609  G4ReactionProduct &leadParticle )
610 {
611  // the following was in GenerateXandPt and TwoCluster
612  // add a parameter to the GenerateXandPt function telling it about the strange particle
613  //
614  // assumes that the original particle was a strange particle
615  //
616  G4bool lead = false;
617  if( (currentParticle.GetMass() >= G4KaonPlus::KaonPlus()->GetPDGMass()) &&
618  (currentParticle.GetDefinition() != G4Proton::Proton()) &&
619  (currentParticle.GetDefinition() != G4Neutron::Neutron()) )
620  {
621  lead = true;
622  leadParticle = currentParticle; // set lead to the incident particle
623  }
624  else if( (targetParticle.GetMass() >= G4KaonPlus::KaonPlus()->GetPDGMass()) &&
625  (targetParticle.GetDefinition() != G4Proton::Proton()) &&
626  (targetParticle.GetDefinition() != G4Neutron::Neutron()) )
627  {
628  lead = true;
629  leadParticle = targetParticle; // set lead to the target particle
630  }
631  return lead;
632 }
633 
634 void FullModelHadronicProcess::Rotate(G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec, G4int &vecLen)
635 {
636  G4double rotation = 2.*CLHEP::pi*G4UniformRand();
637  m_cache = rotation;
638  G4int i;
639  for( i=0; i<vecLen; ++i )
640  {
641  G4ThreeVector momentum = vec[i]->GetMomentum();
642  momentum = momentum.rotate(rotation, m_what);
643  vec[i]->SetMomentum(momentum);
644  }
645 }
646 
647 const G4DynamicParticle* FullModelHadronicProcess::FindRhadron(G4ParticleChange* aParticleChange)
648 {
649  G4int nsec = aParticleChange->GetNumberOfSecondaries();
650  if (nsec==0) return 0;
651  int i = 0;
652  G4bool found = false;
653  while (i!=nsec && !found){
654  // G4cout<<"Checking "<<aParticleChange->GetSecondary(i)->GetDynamicParticle()->GetDefinition()->GetParticleName()<<G4endl;
655  // if (aParticleChange->GetSecondary(i)->GetDynamicParticle()->GetDefinition()->GetParticleType()=="rhadron") found = true;
656  if (dynamic_cast<CustomParticle*>(aParticleChange->GetSecondary(i)->GetDynamicParticle()->GetDefinition())!=0) found = true;
657  i++;
658  }
659  i--;
660  if(found) return aParticleChange->GetSecondary(i)->GetDynamicParticle();
661  return 0;
662 }
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
et
Extra patterns decribing particle interation process.
GeV
#define GeV
Definition: PhysicsAnalysis/TauID/TauAnalysisTools/Root/HelperFunctions.cxx:18
pdg_comparison.sigma
sigma
Definition: pdg_comparison.py:324
VP1PartSpect::Neutron
@ Neutron
Definition: VP1PartSpectFlags.h:25
TRTCalib_cfilter.p1
p1
Definition: TRTCalib_cfilter.py:130
skel.it
it
Definition: skel.GENtoEVGEN.py:407
CustomParticle::GetSpectator
G4ParticleDefinition * GetSpectator()
Definition: CustomParticle.h:45
vec
std::vector< size_t > vec
Definition: CombinationsGeneratorTest.cxx:9
yodamerge_tmp.scale
scale
Definition: yodamerge_tmp.py:138
python.SystemOfUnits.TeV
float TeV
Definition: SystemOfUnits.py:176
CustomParticle::GetCloud
G4ParticleDefinition * GetCloud()
Definition: CustomParticle.h:43
pi
#define pi
Definition: TileMuonFitter.cxx:65
python.SystemOfUnits.MeV
float MeV
Definition: SystemOfUnits.py:172
CxxUtils::vec
typename vecDetail::vec_typedef< T, N >::type vec
Define a nice alias for the vectorized type.
Definition: vec.h:207
python.SystemOfUnits.eV
float eV
Definition: SystemOfUnits.py:173
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:209
ParticleGun_EoverP_Config.momentum
momentum
Definition: ParticleGun_EoverP_Config.py:63
lumiFormat.i
int i
Definition: lumiFormat.py:85
xAOD::rotation
rotation
Definition: TrackSurface_v1.cxx:15
CustomParticle
Definition: CustomParticle.h:18
checkTriggerxAOD.found
found
Definition: checkTriggerxAOD.py:328
CustomParticle.h
VP1PartSpect::Proton
@ Proton
Definition: VP1PartSpectFlags.h:28
pvec
std::array< fp_t, 2 > pvec
Definition: FPGATrackSimLLPDoubletHoughTransformTool.cxx:9
rp
ReadCards * rp
Definition: IReadCards.cxx:26
TRTCalib_cfilter.p0
p0
Definition: TRTCalib_cfilter.py:129
python.SystemOfUnits.m
float m
Definition: SystemOfUnits.py:106
description
std::string description
glabal timer - how long have I taken so far?
Definition: hcg.cxx:88