Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4QStringChipsParticleLevelInterface Class Reference

#include <G4QStringChipsParticleLevelInterface.hh>

+ Inheritance diagram for G4QStringChipsParticleLevelInterface:

Public Member Functions

 G4QStringChipsParticleLevelInterface ()
 
virtual G4HadFinalStateApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &theNucleus)
 
virtual G4ReactionProductVectorPropagate (G4KineticTrackVector *theSecondaries, G4V3DNucleus *theNucleus)
 
- Public Member Functions inherited from G4VIntraNuclearTransportModel
 G4VIntraNuclearTransportModel (const G4String &modelName="CascadeModel", G4VPreCompoundModel *ptr=0)
 
virtual ~G4VIntraNuclearTransportModel ()
 
virtual G4ReactionProductVectorPropagate (G4KineticTrackVector *theSecondaries, G4V3DNucleus *theNucleus)=0
 
void SetDeExcitation (G4VPreCompoundModel *ptr)
 
void Set3DNucleus (G4V3DNucleus *const value)
 
void SetPrimaryProjectile (const G4HadProjectile &aPrimary)
 
const G4StringGetModelName () const
 
virtual void ModelDescription (std::ostream &outFile) const
 
virtual void PropagateModelDescription (std::ostream &outFile) const
 
- Public Member Functions inherited from G4HadronicInteraction
 G4HadronicInteraction (const G4String &modelName="HadronicModel")
 
virtual ~G4HadronicInteraction ()
 
virtual G4HadFinalStateApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)=0
 
virtual G4double SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
virtual G4bool IsApplicable (const G4HadProjectile &, G4Nucleus &)
 
G4double GetMinEnergy () const
 
G4double GetMinEnergy (const G4Material *aMaterial, const G4Element *anElement) const
 
void SetMinEnergy (G4double anEnergy)
 
void SetMinEnergy (G4double anEnergy, const G4Element *anElement)
 
void SetMinEnergy (G4double anEnergy, const G4Material *aMaterial)
 
G4double GetMaxEnergy () const
 
G4double GetMaxEnergy (const G4Material *aMaterial, const G4Element *anElement) const
 
void SetMaxEnergy (const G4double anEnergy)
 
void SetMaxEnergy (G4double anEnergy, const G4Element *anElement)
 
void SetMaxEnergy (G4double anEnergy, const G4Material *aMaterial)
 
const G4HadronicInteractionGetMyPointer () const
 
G4int GetVerboseLevel () const
 
void SetVerboseLevel (G4int value)
 
const G4StringGetModelName () const
 
void DeActivateFor (const G4Material *aMaterial)
 
void ActivateFor (const G4Material *aMaterial)
 
void DeActivateFor (const G4Element *anElement)
 
void ActivateFor (const G4Element *anElement)
 
G4bool IsBlocked (const G4Material *aMaterial) const
 
G4bool IsBlocked (const G4Element *anElement) const
 
void SetRecoilEnergyThreshold (G4double val)
 
G4double GetRecoilEnergyThreshold () const
 
G4bool operator== (const G4HadronicInteraction &right) const
 
G4bool operator!= (const G4HadronicInteraction &right) const
 
virtual const std::pair< G4double, G4doubleGetFatalEnergyCheckLevels () const
 
virtual std::pair< G4double, G4doubleGetEnergyMomentumCheckLevels () const
 
void SetEnergyMomentumCheckLevels (G4double relativeLevel, G4double absoluteLevel)
 
virtual void ModelDescription (std::ostream &outFile) const
 

Additional Inherited Members

- Protected Member Functions inherited from G4VIntraNuclearTransportModel
G4V3DNucleusGet3DNucleus () const
 
G4VPreCompoundModelGetDeExcitation () const
 
const G4HadProjectileGetPrimaryProjectile () const
 
- Protected Member Functions inherited from G4HadronicInteraction
void SetModelName (const G4String &nam)
 
G4bool IsBlocked () const
 
void Block ()
 
- Protected Attributes inherited from G4VIntraNuclearTransportModel
G4String theTransportModelName
 
G4V3DNucleusthe3DNucleus
 
G4VPreCompoundModeltheDeExcitation
 
const G4HadProjectilethePrimaryProjectile
 
- Protected Attributes inherited from G4HadronicInteraction
G4HadFinalState theParticleChange
 
G4int verboseLevel
 
G4double theMinEnergy
 
G4double theMaxEnergy
 
G4bool isBlocked
 

Detailed Description

Definition at line 35 of file G4QStringChipsParticleLevelInterface.hh.

Constructor & Destructor Documentation

◆ G4QStringChipsParticleLevelInterface()

G4QStringChipsParticleLevelInterface::G4QStringChipsParticleLevelInterface ( )

Definition at line 48 of file G4QStringChipsParticleLevelInterface.cc.

49{
50#ifdef debug
51 G4cout<<"G4QStringChipsParticleLevelInterface::Constructor is called"<<G4endl;
52#endif
53 theEnergyLossPerFermi = 1.*GeV;
54 nop = 152; // clusters (A<6)
55 fractionOfSingleQuasiFreeNucleons = 0.5; // It is A-dependent (C=.85, U=.40)
56 fractionOfPairedQuasiFreeNucleons = 0.05;
57 clusteringCoefficient = 5.;
58 temperature = 180.;
59 halfTheStrangenessOfSee = 0.3; // = s/d = s/u
60 etaToEtaPrime = 0.3;
61 fusionToExchange = 1.;
62 theInnerCoreDensityCut = 50.;
63
64 if(getenv("ChipsParameterTuning"))
65 {
66 G4cout << "Please enter the energy loss per fermi in GeV"<<G4endl;
67 G4cin >> theEnergyLossPerFermi;
68 theEnergyLossPerFermi *= GeV;
69 G4cout << "Please enter nop"<<G4endl;
70 G4cin >> nop;
71 G4cout << "Please enter the fractionOfSingleQuasiFreeNucleons"<<G4endl;
72 G4cin >> fractionOfSingleQuasiFreeNucleons;
73 G4cout << "Please enter the fractionOfPairedQuasiFreeNucleons"<<G4endl;
74 G4cin >> fractionOfPairedQuasiFreeNucleons;
75 G4cout << "Please enter the clusteringCoefficient"<<G4endl;
76 G4cin >> clusteringCoefficient;
77 G4cout << "Please enter the temperature"<<G4endl;
78 G4cin >> temperature;
79 G4cout << "Please enter halfTheStrangenessOfSee"<<G4endl;
80 G4cin >> halfTheStrangenessOfSee;
81 G4cout << "Please enter the etaToEtaPrime"<<G4endl;
82 G4cin >> etaToEtaPrime;
83 G4cout << "Please enter the fusionToExchange"<<G4endl;
84 G4cin >> fusionToExchange;
85 G4cout << "Please enter cut-off for calculating the nuclear radius in percent"<<G4endl;
86 G4cin >> theInnerCoreDensityCut;
87 }
88}
#define G4endl
Definition: G4ios.hh:52
#define G4cin
Definition: G4ios.hh:51
G4DLLIMPORT std::ostream G4cout

Member Function Documentation

◆ ApplyYourself()

G4HadFinalState * G4QStringChipsParticleLevelInterface::ApplyYourself ( const G4HadProjectile aTrack,
G4Nucleus theNucleus 
)
virtual

Implements G4HadronicInteraction.

Definition at line 90 of file G4QStringChipsParticleLevelInterface.cc.

92{
93#ifdef debug
94 G4cout<<"G4QStringChipsParticleLevelInterface::ApplyYourself is called"<<G4endl;
95#endif
96 return theModel.ApplyYourself(aTrack, theNucleus);
97}
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &aTargetNucleus, G4HadFinalState *aChange=0)

◆ Propagate()

G4ReactionProductVector * G4QStringChipsParticleLevelInterface::Propagate ( G4KineticTrackVector theSecondaries,
G4V3DNucleus theNucleus 
)
virtual

Implements G4VIntraNuclearTransportModel.

Definition at line 99 of file G4QStringChipsParticleLevelInterface.cc.

101{
102 static const G4double mProt=G4Proton::Proton()->GetPDGMass();
103 static const G4double mNeut=G4Neutron::Neutron()->GetPDGMass();
104 static const G4double mLamb=G4Lambda::Lambda()->GetPDGMass();
105#ifdef debug
106 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate is called"<<G4endl;
107#endif
108 // Protection for non physical conditions
109
110 if(theSecondaries->size() == 1)
111 {
113 G4ReactionProduct* theFastSec;
114 theFastSec = new G4ReactionProduct((*theSecondaries)[0]->GetDefinition());
115 G4LorentzVector current4Mom = (*theSecondaries)[0]->Get4Momentum();
116 theFastSec->SetTotalEnergy(current4Mom.t());
117 theFastSec->SetMomentum(current4Mom.vect());
118 theFastResult->push_back(theFastSec);
119 return theFastResult;
120 //throw G4HadronicException(__FILE__,__LINE__,
121 // "G4QStringChipsParticleLevelInterface: Only one particle from String models!");
122 }
123
124 // target properties needed in constructor of quasmon, and for boosting to
125 // target rest frame
126 // remove all nucleons already involved in STRING interaction, to make the ResidualTarget
127 theNucleus->StartLoop();
128 G4Nucleon * aNucleon;
129 G4int resA = 0;
130 G4int resZ = 0;
131 G4ThreeVector hitMomentum(0,0,0);
132 G4double hitMass = 0;
133 unsigned int hitCount = 0;
134 while((aNucleon = theNucleus->GetNextNucleon()))
135 {
136 if(!aNucleon->AreYouHit())
137 {
138 resA++; // Collect A of the ResidNuc
139 resZ+=G4int (aNucleon->GetDefinition()->GetPDGCharge()); // Collect Z of the ResidNuc
140 }
141 else
142 {
143 hitMomentum += aNucleon->GetMomentum().vect(); // Sum 3-mom of StringHadr's
144 hitMass += aNucleon->GetMomentum().m(); // Sum masses of StringHadrs
145 hitCount ++; // Calculate STRING hadrons
146 }
147 }
148 G4int targetPDGCode = 90000000 + 1000*resZ + (resA-resZ); // PDG of theResidualNucleus
149 G4double targetMass=mNeut;
150 if (!resZ) // Nucleus of only neutrons
151 {
152 if (resA>1) targetMass*=resA;
153 }
154 else targetMass=G4ParticleTable::GetParticleTable()->FindIon(resZ,resA,0,resZ)
155 ->GetPDGMass();
156 G4double targetEnergy = std::sqrt(hitMomentum.mag2()+targetMass*targetMass);
157 // !! @@ Target should be at rest: hitMomentum=(0,0,0) @@ !! M.K. (go to this system)
158 G4LorentzVector targ4Mom(-1.*hitMomentum, targetEnergy);
159
160 // Calculate the mean energy lost
161 std::pair<G4double, G4double> theImpact = theNucleus->RefetchImpactXandY();
162 G4double impactX = theImpact.first;
163 G4double impactY = theImpact.second;
164 G4double impactPar2 = impactX*impactX + impactY*impactY;
165
166 G4double radius2 = theNucleus->GetNuclearRadius(theInnerCoreDensityCut*perCent);
167 radius2 *= radius2;
168 G4double pathlength = 0.;
169#ifdef pdebug
170 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: r="<<std::sqrt(radius2)/fermi
171 <<", b="<<std::sqrt(impactPar2)/fermi<<", R="<<theNucleus->GetOuterRadius()/fermi
172 <<", b/r="<<std::sqrt(impactPar2/radius2)<<G4endl;
173#endif
174 if(radius2 - impactPar2>0) pathlength = 2.*std::sqrt(radius2 - impactPar2);
175 G4double theEnergyLostInFragmentation = theEnergyLossPerFermi*pathlength/fermi;
176
177 // now select all particles in range
178 std::list<std::pair<G4double, G4KineticTrack *> > theSorted; // Output
179 std::list<std::pair<G4double, G4KineticTrack *> >::iterator current; // Input
180 for(unsigned int secondary = 0; secondary<theSecondaries->size(); secondary++)
181 {
182 G4LorentzVector a4Mom = theSecondaries->operator[](secondary)->Get4Momentum();
183#ifdef CHIPSdebug
184 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: ALL STRING particles "
185 << theSecondaries->operator[](secondary)->GetDefinition()->GetPDGCharge()<<" "
186 << theSecondaries->operator[](secondary)->GetDefinition()->GetPDGEncoding()<<" "
187 << a4Mom <<G4endl;
188#endif
189#ifdef pdebug
190 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: in C="
191 <<theSecondaries->operator[](secondary)->GetDefinition()->GetPDGCharge()<<",PDG="
192 <<theSecondaries->operator[](secondary)->GetDefinition()->GetPDGEncoding()
193 <<",4M="<<a4Mom<<", current nS="<<theSorted.size()<<G4endl;
194#endif
195 G4double toSort = a4Mom.rapidity(); // Rapidity is used for the ordering (?!)
196 std::pair<G4double, G4KineticTrack *> it;
197 it.first = toSort;
198 it.second = theSecondaries->operator[](secondary);
199 G4bool inserted = false;
200 for(current = theSorted.begin(); current!=theSorted.end(); current++)
201 {
202 if((*current).first > toSort) // The current is smaller then existing
203 {
204 theSorted.insert(current, it); // It shifts the others up
205 inserted = true;
206 break;
207 }
208 }
209 if(!inserted) theSorted.push_back(it); // It is bigger than any previous
210 }
211
212 G4LorentzVector proj4Mom(0.,0.,0.,0.);
213 G4int nD = 0;
214 G4int nU = 0;
215 G4int nS = 0;
216 G4int nAD = 0;
217 G4int nAU = 0;
218 G4int nAS = 0;
219 std::list<std::pair<G4double,G4KineticTrack*> >::iterator firstEscape=theSorted.begin();
220 G4double runningEnergy = 0;
221 G4int particleCount = 0;
222 G4LorentzVector theLow = (*(theSorted.begin())).second->Get4Momentum();
223 G4LorentzVector theHigh;
224
225#ifdef CHIPSdebug
226 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: CHIPS ENERGY LOST "
227 <<theEnergyLostInFragmentation<<". Sorted rapidities event start"<<G4endl;
228#endif
229#ifdef pdebug
230 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: total CHIPS energy = "
231 <<theEnergyLostInFragmentation<<". Start rapidity sorting nS="<<theSorted.size()
232 <<G4endl;
233#endif
234
235 G4QHadronVector projHV;
236 std::vector<G4QContent> theContents;
237 std::vector<G4LorentzVector*> theMomenta;
239 G4ReactionProduct* theSec;
240 G4KineticTrackVector* secondaries;
241#ifdef pdebug
242 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: Absorption nS="
243 <<theSorted.size()<<G4endl;
244#endif
245 G4bool EscapeExists = false;
246 for(current = theSorted.begin(); current!=theSorted.end(); current++)
247 {
248#ifdef pdebug
249 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: nq="
250 <<(*current).second->GetDefinition()->GetQuarkContent(3)<<", naq="
251 <<(*current).second->GetDefinition()->GetAntiQuarkContent(3)<<", PDG="
252 <<(*current).second->GetDefinition()->GetPDGEncoding()<<",4M="
253 <<(*current).second->Get4Momentum()<<G4endl;
254#endif
255 firstEscape = current; // Remember to make decays for the rest
256 G4KineticTrack* aResult = (*current).second;
257 // This is an old (H.P.) solution, which makes an error in En/Mom conservation
258 //
259 // @@ Now it does not include strange particle for the absorption in nuclei (?!) M.K.
260 //if((*current).second->GetDefinition()->GetQuarkContent(3)!=0 ||
261 // (*current).second->GetDefinition()->GetAntiQuarkContent(3) !=0) // Strange quarks
262 //{
263 // G4ParticleDefinition* pdef = aResult->GetDefinition();
264 // secondaries = NULL;
265 // if ( pdef->GetPDGWidth() > 0 && pdef->GetPDGLifeTime() < 5E-17*s )
266 // secondaries = aResult->Decay(); // @@ Decay of only strange resonances (?!) M.K.
267 // if ( secondaries == NULL ) // No decay
268 // {
269 // theSec = new G4ReactionProduct(aResult->GetDefinition());
270 // G4LorentzVector current4Mom = aResult->Get4Momentum();
271 // current4Mom.boost(targ4Mom.boostVector()); // boost from the targetAtRes system
272 // theSec->SetTotalEnergy(current4Mom.t());
273 // theSec->SetMomentum(current4Mom.vect());
274 // theResult->push_back(theSec);
275 // }
276 // else // The decay happened
277 // {
278 // for (unsigned int aSecondary=0; aSecondary<secondaries->size(); aSecondary++)
279 // {
280 // theSec =
281 // new G4ReactionProduct(secondaries->operator[](aSecondary)->GetDefinition());
282 // G4LorentzVector current4Mom=secondaries->operator[](aSecondary)->Get4Momentum();
283 // current4Mom.boost(targ4Mom.boostVector());
284 // theSec->SetTotalEnergy(current4Mom.t());
285 // theSec->SetMomentum(current4Mom.vect());
286 // theResult->push_back(theSec);
287 // }
288 // std::for_each(secondaries->begin(), secondaries->end(), DeleteKineticTrack());
289 // delete secondaries;
290 // }
291 //}
292 //
293 //runningEnergy += (*current).second->Get4Momentum().t();
294 //if((*current).second->GetDefinition() == G4Proton::Proton())
295 // runningEnergy-=G4Proton::Proton()->GetPDGMass();
296 //if((*current).second->GetDefinition() == G4Neutron::Neutron())
297 // runningEnergy-=G4Neutron::Neutron()->GetPDGMass();
298 //if((*current).second->GetDefinition() == G4Lambda::Lambda())
299 // runningEnergy-=G4Lambda::Lambda()->GetPDGMass();
300 //
301 // New solution starts from here (M.Kossov March 2006) [Strange particles included]
302 runningEnergy += aResult->Get4Momentum().t();
303 G4double charge=aResult->GetDefinition()->GetPDGCharge(); // Charge of the particle
304 G4int strang=aResult->GetDefinition()->GetQuarkContent(3);// Its strangeness
305 G4int baryn=aResult->GetDefinition()->GetBaryonNumber(); // Its baryon number
306 if (baryn>0 && charge>0 && strang<1) runningEnergy-=mProt; // For positive baryons
307 else if(baryn>0 && strang<1) runningEnergy-=mNeut; // For neut/neg baryons
308 else if(baryn>0) runningEnergy-=mLamb; // For strange baryons
309 else if(baryn<0) runningEnergy+=mProt; // For anti-particles
310 // ------------ End of the new solution
311#ifdef CHIPSdebug
312 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: sorted rapidities "
313 <<(*current).second->Get4Momentum().rapidity()<<G4endl;
314#endif
315
316#ifdef pdebug
317 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: E="<<runningEnergy<<", EL="
318 <<theEnergyLostInFragmentation<<G4endl;
319#endif
320 if(runningEnergy > theEnergyLostInFragmentation)
321 {
322 EscapeExists = true;
323 break;
324 }
325#ifdef CHIPSdebug
326 G4cout <<"G4QStringChipsParticleLevelInterface::Propagate: ABSORBED STRING particles "
327 <<(*current).second->GetDefinition()->GetPDGCharge()<<" "
328 << (*current).second->GetDefinition()->GetPDGEncoding()<<" "
329 << (*current).second->Get4Momentum() <<G4endl;
330#endif
331#ifdef pdebug
332 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate:C="
333 <<current->second->GetDefinition()->GetPDGCharge()<<", PDG="
334 <<current->second->GetDefinition()->GetPDGEncoding()<<", 4M="
335 <<current->second->Get4Momentum()<<G4endl;
336#endif
337
338 // projectile 4-momentum in target rest frame needed in constructor of QHadron
339 particleCount++;
340 theHigh = (*current).second->Get4Momentum();
341 proj4Mom = (*current).second->Get4Momentum();
342 proj4Mom.boost(-1.*targ4Mom.boostVector()); // Back to the system of nucleusAtRest
343 nD = (*current).second->GetDefinition()->GetQuarkContent(1);
344 nU = (*current).second->GetDefinition()->GetQuarkContent(2);
345 nS = (*current).second->GetDefinition()->GetQuarkContent(3);
346 nAD = (*current).second->GetDefinition()->GetAntiQuarkContent(1);
347 nAU = (*current).second->GetDefinition()->GetAntiQuarkContent(2);
348 nAS = (*current).second->GetDefinition()->GetAntiQuarkContent(3);
349 G4QContent aProjectile(nD, nU, nS, nAD, nAU, nAS);
350
351#ifdef CHIPSdebug_1
352 G4cout <<G4endl;
353 G4cout <<"G4QStringChipsParticleLevelInterface::Propagate: Quark content: d="<<nD
354 <<", u="<<nU<<", s="<<nS<< "Anti-quark content: anit-d="<<nAD<<", anti-u="<<nAU
355 <<", anti-s="<<nAS<<". G4QContent is constructed"<<endl;
356#endif
357
358 theContents.push_back(aProjectile);
359 G4LorentzVector* aVec = new G4LorentzVector((1./MeV)*proj4Mom); // @@ MeV is basic
360
361#ifdef CHIPSdebug_1
362 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: projectile momentum = "
363 <<*aVec<<G4endl;
364 G4cout << G4endl;
365#endif
366
367 theMomenta.push_back(aVec);
368 }
369 std::vector<G4QContent> theFinalContents;
370 std::vector<G4LorentzVector*> theFinalMomenta;
371
372 for(unsigned int hp = 0; hp<theContents.size(); hp++)
373 {
374 G4QHadron* aHadron = new G4QHadron(theContents[hp], *(theMomenta[hp]) );
375 projHV.push_back(aHadron);
376 }
377
378 // construct the quasmon
379 size_t i;
380 for (i=0; i<theFinalMomenta.size(); i++) delete theFinalMomenta[i];
381 for (i=0; i<theMomenta.size(); i++) delete theMomenta[i];
382 theFinalMomenta.clear();
383 theMomenta.clear();
384
385 G4QNucleus::SetParameters(fractionOfSingleQuasiFreeNucleons,
386 fractionOfPairedQuasiFreeNucleons,
387 clusteringCoefficient,
388 fusionToExchange);
389 G4Quasmon::SetParameters(temperature, halfTheStrangenessOfSee, etaToEtaPrime);
390
391#ifdef CHIPSdebug
392 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: G4QNucleus parameters "
393 <<fractionOfSingleQuasiFreeNucleons<<" "<<fractionOfPairedQuasiFreeNucleons
394 <<" "<<clusteringCoefficient<<G4endl;
395 G4cout<<"G4Quasmon parameters "<<temperature<<" "<<halfTheStrangenessOfSee<<" "
396 <<etaToEtaPrime << G4endl;
397 G4cout<<"The Target PDG code = "<<targetPDGCode<<G4endl;
398 G4cout<<"The projectile momentum = "<<1./MeV*proj4Mom<<G4endl;
399 G4cout<<"The target momentum = "<<1./MeV*targ4Mom<<G4endl;
400#endif
401
402 // now call chips with this info in place
403 G4QHadronVector* output = 0;
404 if (particleCount!=0 && resA!=0)
405 {
406 // G4QCHIPSWorld aWorld(nop); // Create CHIPS World of nop particles
408 G4QEnvironment* pan= new G4QEnvironment(projHV, targetPDGCode);
409#ifdef pdebug
410 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: CHIPS fragmentation, rA="
411 <<resA<<", #AbsPt="<<particleCount<<G4endl;
412#endif
413 try
414 {
415 output = pan->Fragment(); // The main fragmentation member function
416 }
417 catch(G4HadronicException& aR)
418 {
419 G4cerr << "Exception thrown of G4QStringChipsParticleLevelInterface "<<G4endl;
420 G4cerr << " targetPDGCode = "<< targetPDGCode <<G4endl;
421 G4cerr << " The projectile momentum = "<<1./MeV*proj4Mom<<G4endl;
422 G4cerr << " The target momentum = "<<1./MeV*targ4Mom<<G4endl<<G4endl;
423 G4cerr << " Dumping the information in the pojectile list"<<G4endl;
424 for(i=0; i< projHV.size(); i++)
425 {
426 G4cerr <<" Incoming 4-momentum and PDG code of "<<i<<"'th hadron: "
427 <<" "<< projHV[i]->Get4Momentum()<<" "<<projHV[i]->GetPDGCode()<<G4endl;
428 }
429 throw;
430 }
431 // clean up particles
432 std::for_each(projHV.begin(), projHV.end(), DeleteQHadron());
433 projHV.clear();
434 delete pan;
435 }
436 else
437 {
438#ifdef pdebug
439 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: NO CHIPS fragmentation, rA="
440 <<resA<<", #AbsPt="<<particleCount<<G4endl;
441#endif
442 output = new G4QHadronVector;
443 }
444 // Fill the result.
445#ifdef CHIPSdebug
446 G4cout << "NEXT EVENT, EscapeExists="<<EscapeExists<<G4endl;
447#endif
448
449 // first decay and add all escaping particles.
450 if (EscapeExists) for (current = firstEscape; current!=theSorted.end(); current++)
451 {
452 G4KineticTrack* aResult = (*current).second;
453 G4ParticleDefinition* pdef=aResult->GetDefinition();
454 secondaries = NULL;
455 if(pdef->GetPDGWidth() > 0 && pdef->GetPDGLifeTime() < 5E-17*s )
456 {
457 secondaries = aResult->Decay(); // @@ Uses standard Decay, which is now wrong!
458 }
459 if ( secondaries == NULL )
460 {
461 theSec = new G4ReactionProduct(aResult->GetDefinition());
462 G4LorentzVector current4Mom = aResult->Get4Momentum();
463 current4Mom.boost(targ4Mom.boostVector());
464 theSec->SetTotalEnergy(current4Mom.t());
465 theSec->SetMomentum(current4Mom.vect());
466#ifdef pdebug
467 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: *OUT* QGS stable PDG="
468 <<aResult->GetDefinition()->GetPDGEncoding()<<",4M="<<current4Mom<<G4endl;
469#endif
470 theResult->push_back(theSec);
471 }
472 else
473 {
474 for (unsigned int aSecondary=0; aSecondary<secondaries->size(); aSecondary++)
475 {
476 theSec=new G4ReactionProduct(secondaries->operator[](aSecondary)->GetDefinition());
477 G4LorentzVector current4Mom = secondaries->operator[](aSecondary)->Get4Momentum();
478 current4Mom.boost(targ4Mom.boostVector());
479 theSec->SetTotalEnergy(current4Mom.t());
480 theSec->SetMomentum(current4Mom.vect());
481#ifdef pdebug
482 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: *OUT* QGS decay PDG="
483 <<secondaries->operator[](aSecondary)->GetDefinition()->GetPDGEncoding()
484 <<",4M="<<current4Mom<<G4endl;
485#endif
486 theResult->push_back(theSec);
487 }
488 std::for_each(secondaries->begin(), secondaries->end(), DeleteKineticTrack());
489 delete secondaries;
490 }
491 }
492 std::for_each(theSecondaries->begin(), theSecondaries->end(), DeleteKineticTrack());
493 delete theSecondaries;
494
495 // now add the quasmon output
496 G4int maxParticle=output->size();
497#ifdef CHIPSdebug
498 G4cout << "Number of particles from string"<<theResult->size()<<G4endl;
499 G4cout << "Number of particles from chips"<<maxParticle<<G4endl;
500#endif
501#ifdef pdebug
502 G4cout << "Number of particles from QGS="<<theResult->size()<<G4endl;
503 G4cout << "Number of particles from CHIPS="<<maxParticle<<G4endl;
504#endif
505 if(maxParticle) for(G4int particle = 0; particle < maxParticle; particle++)
506 {
507 if(output->operator[](particle)->GetNFragments() != 0)
508 {
509 delete output->operator[](particle);
510 continue;
511 }
512 G4int pdgCode = output->operator[](particle)->GetPDGCode();
513
514
515#ifdef CHIPSdebug
516 G4cerr << "PDG code of chips particle = "<<pdgCode<<G4endl;
517#endif
518
519 G4ParticleDefinition * theDefinition;
520 // Note that I still have to take care of strange nuclei
521 // For this I need the mass calculation, and a changed interface
522 // for ion-table ==> work for Hisaya @@@@@@@
523 // Then I can sort out the pdgCode. I also need a decau process
524 // for strange nuclei; may be another chips interface
525 if(pdgCode>90000000)
526 {
527 G4int aZ = (pdgCode-90000000)/1000;
528 if (aZ>1000) aZ=aZ%1000; // patch for strange nuclei, to be repaired @@@@
529 G4int anN = pdgCode-90000000-1000*aZ;
530 if(anN>1000) anN=anN%1000; // patch for strange nuclei, to be repaired @@@@
531 if(pdgCode==91000000) theDefinition = G4Lambda::LambdaDefinition();
532 else if(pdgCode==92000000) theDefinition = G4Lambda::LambdaDefinition();
533 else if(pdgCode==93000000) theDefinition = G4Lambda::LambdaDefinition();
534 else if(pdgCode==94000000) theDefinition = G4Lambda::LambdaDefinition();
535 else if(pdgCode==95000000) theDefinition = G4Lambda::LambdaDefinition();
536 else if(pdgCode==96000000) theDefinition = G4Lambda::LambdaDefinition();
537 else if(pdgCode==97000000) theDefinition = G4Lambda::LambdaDefinition();
538 else if(pdgCode==98000000) theDefinition = G4Lambda::LambdaDefinition();
539 else if(aZ == 0 && anN == 1) theDefinition = G4Neutron::Neutron();
540 else theDefinition = G4ParticleTable::GetParticleTable()->FindIon(aZ,anN+aZ,0,aZ);
541 }
542 else theDefinition = G4ParticleTable::GetParticleTable()->FindParticle(pdgCode);
543
544#ifdef CHIPSdebug
545 G4cout << "Particle code produced = "<< pdgCode <<G4endl;
546#endif
547
548 if(theDefinition)
549 {
550 theSec = new G4ReactionProduct(theDefinition);
551 G4LorentzVector current4Mom = output->operator[](particle)->Get4Momentum();
552 current4Mom.boost(targ4Mom.boostVector());
553 theSec->SetTotalEnergy(current4Mom.t());
554 theSec->SetMomentum(current4Mom.vect());
555#ifdef pdebug
556 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: *OUT* CHIPS PDG="
557 <<theDefinition->GetPDGEncoding()<<",4M="<<current4Mom<<G4endl;
558#endif
559 theResult->push_back(theSec);
560 }
561 else
562 {
563 G4cerr << G4endl<<"WARNING: "<<G4endl;
564 G4cerr << "Getting unknown pdgCode from chips in ParticleLevelInterface"<<G4endl;
565 G4cerr << "skipping particle with pdgCode = "<<pdgCode<<G4endl<<G4endl;
566 }
567
568#ifdef CHIPSdebug
569 G4cout <<"CHIPS particles "<<theDefinition->GetPDGCharge()<<" "
570 << theDefinition->GetPDGEncoding()<<" "
571 << output->operator[](particle)->Get4Momentum() <<G4endl;
572#endif
573
574 delete output->operator[](particle);
575 }
576 else
577 {
578 if(resA>0)
579 {
580 G4ParticleDefinition* theDefinition = G4Neutron::Neutron();
581 if(resA==1) // The residual nucleus at rest must be added to conserve BaryN & Charge
582 {
583 if(resZ == 1) theDefinition = G4Proton::Proton();
584 }
585 else theDefinition = G4ParticleTable::GetParticleTable()->FindIon(resZ,resA,0,resZ);
586 theSec = new G4ReactionProduct(theDefinition);
587 theSec->SetTotalEnergy(theDefinition->GetPDGMass());
588 theSec->SetMomentum(G4ThreeVector(0.,0.,0.));
589 theResult->push_back(theSec);
590 if(!resZ && resA>0) for(G4int ni=1; ni<resA; ni++)
591 {
592 theSec = new G4ReactionProduct(theDefinition);
593 theSec->SetTotalEnergy(theDefinition->GetPDGMass());
594 theSec->SetMomentum(G4ThreeVector(0.,0.,0.));
595 theResult->push_back(theSec);
596 }
597 }
598 }
599 delete output;
600
601#ifdef CHIPSdebug
602 G4cout << "Number of particles"<<theResult->size()<<G4endl;
603 G4cout << G4endl;
604 G4cout << "QUASMON preparation info "
605 << 1./MeV*proj4Mom<<" "
606 << 1./MeV*targ4Mom<<" "
607 << nD<<" "<<nU<<" "<<nS<<" "<<nAD<<" "<<nAU<<" "<<nAS<<" "
608 << hitCount<<" "
609 << particleCount<<" "
610 << theLow<<" "
611 << theHigh<<" "
612 << G4endl;
613#endif
614
615 return theResult;
616}
CLHEP::HepLorentzVector G4LorentzVector
std::vector< G4QHadron * > G4QHadronVector
std::vector< G4ReactionProduct * > G4ReactionProductVector
CLHEP::Hep3Vector G4ThreeVector
double G4double
Definition: G4Types.hh:64
int G4int
Definition: G4Types.hh:66
bool G4bool
Definition: G4Types.hh:67
G4DLLIMPORT std::ostream G4cerr
HepLorentzVector & boost(double, double, double)
Hep3Vector vect() const
G4KineticTrackVector * Decay()
G4ParticleDefinition * GetDefinition() const
const G4LorentzVector & Get4Momentum() const
static G4Lambda * Lambda()
Definition: G4Lambda.cc:108
static G4Lambda * LambdaDefinition()
Definition: G4Lambda.cc:103
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4bool AreYouHit() const
Definition: G4Nucleon.hh:97
virtual G4ParticleDefinition * GetDefinition() const
Definition: G4Nucleon.hh:85
const G4LorentzVector & GetMomentum() const
Definition: G4Nucleon.hh:71
G4int GetQuarkContent(G4int flavor) const
G4double GetPDGWidth() const
G4double GetPDGCharge() const
G4double GetPDGLifeTime() const
G4ParticleDefinition * FindIon(G4int atomicNumber, G4int atomicMass, G4double excitationEnergy)
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
static G4ParticleTable * GetParticleTable()
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4QCHIPSWorld * Get()
G4QParticleVector * GetParticles(G4int nOfParts=0)
G4QHadronVector * Fragment()
static void SetParameters(G4double fN=.1, G4double fD=.05, G4double cP=4., G4double mR=1., G4double nD=.8 *CLHEP::fermi)
Definition: G4QNucleus.cc:347
static void SetParameters(G4double temper=180., G4double ssin2g=.3, G4double etaetap=.3)
Definition: G4Quasmon.cc:192
void SetMomentum(const G4double x, const G4double y, const G4double z)
void SetTotalEnergy(const G4double en)
std::pair< G4double, G4double > RefetchImpactXandY()
Definition: G4V3DNucleus.hh:77
virtual G4double GetOuterRadius()=0
virtual G4Nucleon * GetNextNucleon()=0
virtual G4bool StartLoop()=0
virtual G4double GetNuclearRadius()=0

The documentation for this class was generated from the following files: