Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
G4StringChipsParticleLevelInterface.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// Short description: Interface of QGSC to CHIPS (SoftHadrons)
27// ----------------------------------------------------------------
28//
29
30//#define debug
31//#define trapdebug
32//#define pdebug
33//#define ppdebug
34
35#include <utility>
36#include <list>
37#include <vector>
38
40#include "globals.hh"
41#include "G4SystemOfUnits.hh"
43#include "G4Nucleon.hh"
44#include "G4Proton.hh"
45#include "G4Neutron.hh"
46#include "G4LorentzRotation.hh"
48// #define CHIPSdebug
49// #define CHIPSdebug_1
50
51#ifdef hdebug_SCPLI
52const G4int G4StringChipsParticleLevelInterface::nbh=200;
53G4double G4StringChipsParticleLevelInterface::bhmax=20.;
54G4double G4StringChipsParticleLevelInterface::ehmax=20.;
55G4double G4StringChipsParticleLevelInterface::bhdb=0.;
56G4double G4StringChipsParticleLevelInterface::ehde=0.;
57G4double G4StringChipsParticleLevelInterface::toth=0.;
58G4int G4StringChipsParticleLevelInterface::bover=0;
59G4int G4StringChipsParticleLevelInterface::eover=0;
60G4int* G4StringChipsParticleLevelInterface::bhis =
61 new G4int[G4StringChipsParticleLevelInterface::nbh];
62G4int* G4StringChipsParticleLevelInterface::ehis =
63 new G4int[G4StringChipsParticleLevelInterface::nbh];
64#endif
65
67{
68#ifdef debug
69 G4cout<<"G4StringChipsParticleLevelInterface::Constructor is called"<<G4endl;
70#endif
71 //theEnergyLossPerFermi = 1.*GeV;
72 theEnergyLossPerFermi = 1.5*GeV;
73 nop = 152; // clusters (A<6)
74 fractionOfSingleQuasiFreeNucleons = 0.5; // It is A-dependent (C=.85, U=.40) - M.K.
75 fractionOfPairedQuasiFreeNucleons = 0.05;
76 clusteringCoefficient = 5.;
77 temperature = 180.;
78 halfTheStrangenessOfSee = 0.3; // = s/d = s/u
79 etaToEtaPrime = 0.3;
80 fusionToExchange = 1.;
81 //theInnerCoreDensityCut = 50.;
82 theInnerCoreDensityCut = 70.;
83
84 if(getenv("ChipsParameterTuning"))
85 {
86 G4cout << "Please enter the energy loss per fermi in GeV"<<G4endl;
87 G4cin >> theEnergyLossPerFermi;
88 theEnergyLossPerFermi *= GeV;
89 G4cout << "Please enter nop"<<G4endl;
90 G4cin >> nop;
91 G4cout << "Please enter the fractionOfSingleQuasiFreeNucleons"<<G4endl;
92 G4cin >> fractionOfSingleQuasiFreeNucleons;
93 G4cout << "Please enter the fractionOfPairedQuasiFreeNucleons"<<G4endl;
94 G4cin >> fractionOfPairedQuasiFreeNucleons;
95 G4cout << "Please enter the clusteringCoefficient"<<G4endl;
96 G4cin >> clusteringCoefficient;
97 G4cout << "Please enter the temperature"<<G4endl;
98 G4cin >> temperature;
99 G4cout << "Please enter halfTheStrangenessOfSee"<<G4endl;
100 G4cin >> halfTheStrangenessOfSee;
101 G4cout << "Please enter the etaToEtaPrime"<<G4endl;
102 G4cin >> etaToEtaPrime;
103 G4cout << "Please enter the fusionToExchange"<<G4endl;
104 G4cin >> fusionToExchange;
105 G4cout << "Please enter the cut-off for calculating the nuclear radius in percent"<<G4endl;
106 G4cin >> theInnerCoreDensityCut;
107 }
108}
109
111ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& theNucleus)
112{
113#ifdef debug
114 G4cout<<"G4StringChipsParticleLevelInterface::ApplyYourself is called"<<G4endl;
115#endif
116 return theModel.ApplyYourself(aTrack, theNucleus);
117}
118
120Propagate(G4KineticTrackVector* theSecondaries, G4V3DNucleus* theNucleus)
121{
122 static const G4double mProt=G4Proton::Proton()->GetPDGMass();
123 static const G4double mNeut=G4Neutron::Neutron()->GetPDGMass();
124 static const G4double mLamb=G4Lambda::Lambda()->GetPDGMass();
125 static const G4double mKChg=G4KaonPlus::KaonPlus()->GetPDGMass();
126 static const G4double mKZer=G4KaonZero::KaonZero()->GetPDGMass();
127 static const G4double mPiCh=G4PionMinus::PionMinus()->GetPDGMass();
128 static const G4int pcl=4; // clusterization parameter for Energy Flow
129 static const G4QContent ProtQC(1,2,0,0,0,0);
130 static const G4QContent NeutQC(2,1,0,0,0,0);
131 static const G4QContent LambQC(1,1,1,0,0,0);
132 static const G4QContent KPlsQC(0,1,0,0,0,1);
133 static const G4QContent KMinQC(0,0,1,0,1,0);
134 static const G4QContent AKZrQC(1,0,0,0,0,1);
135 static const G4QContent KZerQC(1,0,0,0,0,1);
136 static const G4QContent PiMiQC(1,0,0,0,1,0);
137 static const G4QContent PiPlQC(0,1,0,1,0,0);
138#ifdef debug
139 G4cout<<"G4StringChipsParticleLevelInterface::Propagate is called"<<G4endl;
140#endif
141 // Protection for non physical conditions
142
143 if(theSecondaries->size() == 1)
144 {
146 G4ReactionProduct* theFastSec;
147 theFastSec = new G4ReactionProduct((*theSecondaries)[0]->GetDefinition());
148 G4LorentzVector current4Mom = (*theSecondaries)[0]->Get4Momentum();
149 theFastSec->SetTotalEnergy(current4Mom.t());
150 theFastSec->SetMomentum(current4Mom.vect());
151 theFastResult->push_back(theFastSec);
152 return theFastResult;
153 //throw G4HadronicException(__FILE__,__LINE__,
154 // "G4StringChipsParticleLevelInterface: Only one particle from String models!");
155 }
156
157 // target properties needed in constructor of quasmon, and for boosting to
158 // target rest frame
159 // remove all nucleons already involved in STRING interaction, to make the ResidualTarget
160 theNucleus->StartLoop();
161 G4Nucleon * aNucleon;
162 G4int resA = 0;
163 G4int resZ = 0;
164 G4ThreeVector hitMomentum(0,0,0);
165 G4double hitMass = 0;
166 unsigned int hitCount = 0;
167 while((aNucleon = theNucleus->GetNextNucleon()))
168 {
169 if(!aNucleon->AreYouHit())
170 {
171 resA++; // Collect A of the ResidNuc
172 resZ+=G4int (aNucleon->GetDefinition()->GetPDGCharge()); // Collect Z of the ResidNuc
173 }
174 else
175 {
176 hitMomentum += aNucleon->GetMomentum().vect(); // Sum 3-mom of StringHadr's
177 hitMass += aNucleon->GetMomentum().m(); // Sum masses of StringHadrs
178 hitCount ++; // Calculate STRING hadrons
179 }
180 }
181 G4int targetPDGCode = 90000000 + 1000*resZ + (resA-resZ); // PDG of theResidualNucleus
182 G4double targetMass=mNeut;
183 if (!resZ) // Nucleus of only neutrons
184 {
185 if (resA>1) targetMass*=resA;
186 }
187 else targetMass=G4ParticleTable::GetParticleTable()->FindIon(resZ,resA,0,resZ)
188 ->GetPDGMass();
189 G4double targetEnergy = std::sqrt(hitMomentum.mag2()+targetMass*targetMass);
190 // !! @@ Target should be at rest: hitMomentum=(0,0,0) @@ !! M.K. (go to this system)
191 G4LorentzVector targ4Mom(-1.*hitMomentum, targetEnergy);
192
193 // Calculate the mean energy lost
194 std::pair<G4double, G4double> theImpact = theNucleus->RefetchImpactXandY();
195 G4double impactX = theImpact.first;
196 G4double impactY = theImpact.second;
197 G4double impactPar2 = impactX*impactX + impactY*impactY;
198 G4double radius2 = theNucleus->GetNuclearRadius(theInnerCoreDensityCut*perCent);
199 //G4double radius2 = theNucleus->GetNuclearRadius(theInnerCoreDensityCut*perCent);
200 radius2 *= radius2;
201 G4double pathlength = 0.;
202#ifdef ppdebug
203 G4cout<<"G4StringChipsParticleLevelInterface::Propagate: r="<<std::sqrt(radius2)/fermi
204 <<", b="<<std::sqrt(impactPar2)/fermi<<", R="<<theNucleus->GetOuterRadius()/fermi
205 <<", b/r="<<std::sqrt(impactPar2/radius2)<<G4endl;
206#endif
207 if(radius2 - impactPar2>0) pathlength = 2.*std::sqrt(radius2 - impactPar2);
208 //pathlength = 27.; // *** Only temporary *** Always CHIPS
209#ifdef hdebug_SCPLI
210 toth+=1.; // increment total number of measurements
211 G4double bfm=std::sqrt(impactPar2)/fermi; // impact parameter
212 G4double efm=pathlength/fermi; // energy absorption length
213 G4int nbi=static_cast<G4int>(bfm/bhdb);
214 G4int nei=static_cast<G4int>(efm/ehde);
215 if(nbi<nbh) bhis[nbi]++;
216 else bover++;
217 if(nei<nbh) ehis[nei]++;
218 else eover++;
219#endif
220 G4double theEnergyLostInFragmentation = theEnergyLossPerFermi*pathlength/fermi;
221
222 // now select all particles in range
223 std::list<std::pair<G4double, G4KineticTrack *> > theSorted; // Output
224 std::list<std::pair<G4double, G4KineticTrack *> >::iterator current; // Input
225 for(unsigned int secondary = 0; secondary<theSecondaries->size(); secondary++)
226 {
227 G4LorentzVector a4Mom = theSecondaries->operator[](secondary)->Get4Momentum();
228#ifdef CHIPSdebug
229 G4cout<<"G4StringChipsParticleLevelInterface::Propagate: ALL STRING particles "
230 << theSecondaries->operator[](secondary)->GetDefinition()->GetPDGCharge()<<" "
231 << theSecondaries->operator[](secondary)->GetDefinition()->GetPDGEncoding()<<" "
232 << a4Mom <<G4endl;
233#endif
234#ifdef pdebug
235 G4cout<<"G4StringChipsParticleLevelInterface::Propagate: in C="
236 <<theSecondaries->operator[](secondary)->GetDefinition()->GetPDGCharge()<<",PDG="
237 <<theSecondaries->operator[](secondary)->GetDefinition()->GetPDGEncoding()
238 <<",4M="<<a4Mom<<", current nS="<<theSorted.size()<<G4endl;
239#endif
240 G4double toSort = a4Mom.rapidity(); // Rapidity is used for the ordering (?!)
241 std::pair<G4double, G4KineticTrack *> it;
242 it.first = toSort;
243 it.second = theSecondaries->operator[](secondary);
244 G4bool inserted = false;
245 for(current = theSorted.begin(); current!=theSorted.end(); current++)
246 {
247 if((*current).first > toSort) // The current is smaller then existing
248 {
249 theSorted.insert(current, it); // It shifts the others up
250 inserted = true;
251 break;
252 }
253 }
254 if(!inserted) theSorted.push_back(it); // It is bigger than any previous (@@Check this)
255 }
256
257 G4LorentzVector proj4Mom(0.,0.,0.,0.);
258 G4int nD = 0;
259 G4int nU = 0;
260 G4int nS = 0;
261 G4int nAD = 0;
262 G4int nAU = 0;
263 G4int nAS = 0;
264 std::list<std::pair<G4double,G4KineticTrack*> >::iterator firstEscape=theSorted.begin();
265 G4double runningEnergy = 0;
266 G4int particleCount = 0;
267 G4LorentzVector theLow = (*(theSorted.begin())).second->Get4Momentum();
268 G4LorentzVector theHigh;
269
270#ifdef CHIPSdebug
271 G4cout<<"G4StringChipsParticleLevelInterface::Propagate: CHIPS ENERGY LOST "
272 <<theEnergyLostInFragmentation<<". Sorted rapidities event start"<<G4endl;
273#endif
274#ifdef pdebug
275 G4cout<<"G4StringChipsParticleLevelInterface::Propagate: total CHIPS energy = "
276 <<theEnergyLostInFragmentation<<". Start rapidity sorting nS="<<theSorted.size()
277 <<G4endl;
278#endif
279
280 G4QHadronVector projHV;
281 std::vector<G4QContent> theContents;
282 std::vector<G4LorentzVector*> theMomenta;
284 G4ReactionProduct* theSec;
285 G4KineticTrackVector* secondaries;
286 G4KineticTrackVector* secsec;
287#ifdef pdebug
288 G4cout<<"G4StringChipsParticleLevelInterface::Propagate: Absorption nS="
289 <<theSorted.size()<<G4endl;
290#endif
291 G4bool EscapeExists = false;
292 for(current = theSorted.begin(); current!=theSorted.end(); current++)
293 {
294#ifdef pdebug
295 G4cout<<"G4StringChipsParticleLevelInterface::Propagate: nq="
296 <<(*current).second->GetDefinition()->GetQuarkContent(3)<<", naq="
297 <<(*current).second->GetDefinition()->GetAntiQuarkContent(3)<<", PDG="
298 <<(*current).second->GetDefinition()->GetPDGEncoding()<<",4M="
299 <<(*current).second->Get4Momentum()<<G4endl;
300#endif
301 firstEscape = current; // Remember to make decays for the rest
302 G4KineticTrack* aResult = (*current).second;
303 // This is an old (H.P.) solution, which makes an error in En/Mom conservation
304 //
305 // @@ Now it does not include strange particle for the absorption in nuclei (?!) M.K.
306 //if((*current).second->GetDefinition()->GetQuarkContent(3)!=0 ||
307 // (*current).second->GetDefinition()->GetAntiQuarkContent(3) !=0) // Strange quarks
308 //{
309 // G4ParticleDefinition* pdef = aResult->GetDefinition();
310 // secondaries = NULL;
311 // if ( pdef->GetPDGWidth() > 0 && pdef->GetPDGLifeTime() < 5E-17*s )
312 // secondaries = aResult->Decay(); // @@ Decay of only strange resonances (?!) M.K.
313 // if ( secondaries == NULL ) // No decay
314 // {
315 // theSec = new G4ReactionProduct(aResult->GetDefinition());
316 // G4LorentzVector current4Mom = aResult->Get4Momentum();
317 // current4Mom.boost(targ4Mom.boostVector()); // boost from the targetAtRes system
318 // theSec->SetTotalEnergy(current4Mom.t());
319 // theSec->SetMomentum(current4Mom.vect());
320 // theResult->push_back(theSec);
321 // }
322 // else // The decay happened
323 // {
324 // for (unsigned int aSecondary=0; aSecondary<secondaries->size(); aSecondary++)
325 // {
326 // theSec =
327 // new G4ReactionProduct(secondaries->operator[](aSecondary)->GetDefinition());
328 // G4LorentzVector current4Mom =secondaries->operator[](aSecondary)->Get4Momentum();
329 // current4Mom.boost(targ4Mom.boostVector());
330 // theSec->SetTotalEnergy(current4Mom.t());
331 // theSec->SetMomentum(current4Mom.vect());
332 // theResult->push_back(theSec);
333 // }
334 // std::for_each(secondaries->begin(), secondaries->end(), DeleteKineticTrack());
335 // delete secondaries;
336 // }
337 //}
338 //
339 //runningEnergy += (*current).second->Get4Momentum().t();
340 //if((*current).second->GetDefinition() == G4Proton::Proton())
341 // runningEnergy-=G4Proton::Proton()->GetPDGMass();
342 //if((*current).second->GetDefinition() == G4Neutron::Neutron())
343 // runningEnergy-=G4Neutron::Neutron()->GetPDGMass();
344 //if((*current).second->GetDefinition() == G4Lambda::Lambda())
345 // runningEnergy-=G4Lambda::Lambda()->GetPDGMass();
346 //
347 // New solution starts from here (M.Kossov March 2006) [Strange particles included]
348 runningEnergy += aResult->Get4Momentum().t();
349 G4double charge=aResult->GetDefinition()->GetPDGCharge(); // Charge of the particle
350 G4int strang=aResult->GetDefinition()->GetQuarkContent(3);// Its strangeness
351 G4int baryn=aResult->GetDefinition()->GetBaryonNumber(); // Its baryon number
352 if (baryn>0 && charge>0 && strang<1) runningEnergy-=mProt; // For positive baryons
353 else if(baryn>0 && strang<1) runningEnergy-=mNeut; // For neut/neg baryons
354 else if(baryn>0) runningEnergy-=mLamb; // For strange baryons
355 else if(baryn<0) runningEnergy+=mProt; // For anti-particles
356 // ------------ End of the new solution
357#ifdef CHIPSdebug
358 G4cout<<"G4StringChipsParticleLevelInterface::Propagate: sorted rapidities "
359 <<(*current).second->Get4Momentum().rapidity()<<G4endl;
360#endif
361
362#ifdef pdebug
363 G4cout<<"G4StringChipsParticleLevelInterface::Propagate: E="<<runningEnergy<<", EL="
364 <<theEnergyLostInFragmentation<<G4endl;
365#endif
366 if(runningEnergy > theEnergyLostInFragmentation)
367 {
368 EscapeExists = true;
369 break;
370 }
371#ifdef CHIPSdebug
372 G4cout <<"G4StringChipsParticleLevelInterface::Propagate: ABSORBED STRING particles "
373 <<(*current).second->GetDefinition()->GetPDGCharge()<<" "
374 << (*current).second->GetDefinition()->GetPDGEncoding()<<" "
375 << (*current).second->Get4Momentum() <<G4endl;
376#endif
377#ifdef pdebug
378 G4cout<<"G4StringChipsParticleLevelInterface::Propagate:C="
379 <<current->second->GetDefinition()->GetPDGCharge()<<", PDG="
380 <<current->second->GetDefinition()->GetPDGEncoding()<<", 4M="
381 <<current->second->Get4Momentum()<<G4endl;
382#endif
383
384 // projectile 4-momentum in target rest frame needed in constructor of QHadron
385 particleCount++;
386 theHigh = (*current).second->Get4Momentum();
387 proj4Mom = (*current).second->Get4Momentum();
388 proj4Mom.boost(-1.*targ4Mom.boostVector()); // Back to the system of nucleusAtRest
389 nD = (*current).second->GetDefinition()->GetQuarkContent(1);
390 nU = (*current).second->GetDefinition()->GetQuarkContent(2);
391 nS = (*current).second->GetDefinition()->GetQuarkContent(3);
392 nAD = (*current).second->GetDefinition()->GetAntiQuarkContent(1);
393 nAU = (*current).second->GetDefinition()->GetAntiQuarkContent(2);
394 nAS = (*current).second->GetDefinition()->GetAntiQuarkContent(3);
395 G4QContent aProjectile(nD, nU, nS, nAD, nAU, nAS);
396
397#ifdef CHIPSdebug_1
398 G4cout <<G4endl;
399 G4cout <<"G4StringChipsParticleLevelInterface::Propagate: Quark content: d="<<nD
400 <<", u="<<nU<<", s="<<nS<< "Anti-quark content: anit-d="<<nAD<<", anti-u="<<nAU
401 <<", anti-s="<<nAS<<". G4QContent is constructed"<<endl;
402#endif
403
404 theContents.push_back(aProjectile);
405 G4LorentzVector* aVec = new G4LorentzVector((1./MeV)*proj4Mom); // @@ MeV is basic
406
407#ifdef CHIPSdebug_1
408 G4cout<<"G4StringChipsParticleLevelInterface::Propagate: projectile momentum = "
409 <<*aVec<<G4endl;
410 G4cout << G4endl;
411#endif
412
413 theMomenta.push_back(aVec);
414 }
415 std::vector<G4QContent> theFinalContents;
416 std::vector<G4LorentzVector*> theFinalMomenta;
417
418 // Multiquasmon case: each particle creates a quasmon
419 //for(unsigned int hp = 0; hp<theContents.size(); hp++)
420 //{
421 // G4QHadron* aHadron = new G4QHadron(theContents[hp], *(theMomenta[hp]) );
422 // projHV.push_back(aHadron);
423 //}
424 // Energy flow: one Quasmon for each B>0 collection ----------
425 G4QContent EnFlowQC(0,0,0,0,0,0);
426 G4LorentzVector EnFlow4M(0.,0.,0.,0.);
427 //G4bool empty=true;
428 G4int barys=0;
429 G4int stras=0;
430 G4int chars=0;
431 for(G4int hp = theContents.size()-1; hp>=0; hp--)
432 {
433 G4QContent curQC=theContents[hp];
434 G4int baryn = curQC.GetBaryonNumber();
435 G4int stran = curQC.GetStrangeness();
436 G4int charg = curQC.GetCharge();
437 EnFlowQC += curQC; // Keep collecting energy flow
438 EnFlow4M += *(theMomenta[hp]);
439 barys += baryn; // Collected baryon number
440 stras += stran; // Collected strangeness
441 chars += charg; // Collected charge
442 //empty = false;
443 }
444 if(barys>pcl) // Split in two or more parts (to survive!)
445 {
446 G4int nprt=(barys-1)/pcl+1; // Number of parts (pcl=4: 2:5-8,3:9-12...)
447 G4int curb=barys;
448 while (nprt>0)
449 {
450 nprt--; // One part is going to be created
451 G4int brnm=pcl; // Baryon number of splitting part
452 curb-=brnm; // The residual baryon number
453 G4double prtM=0.; // The resulting GS mass of the part
454 G4double resM=0.; // The resulting GS mass of the residual
455 G4QContent prtQC(0,0,0,0,0,0); // The resulting Quark Content of the part
456 G4int strm=0; // Max strangeness per part (stras=0)
457 if(stras>0) strm=(stras-1)/nprt+1; // Max strangeness per part (stras>0)
458 else if(stras<0) strm=(stras+1)/nprt-1; // Max strangeness per part (stras<0)
459 G4int chgm=0; // Max charge per part (chars=0)
460 if(stras>0) chgm=(chars-1)/nprt+1; // Max strangeness per part (chars>0)
461 else if(stras<0) chgm=(chars+1)/nprt-1; // Max strangeness per part (chars<0)
462 // ---> calculate proposed separated part
463 //@@ Convert it to a CHIPS function (Which class? G4QH::Conctruct?)
464 if(!strm) // --> The total strangness = 0 (n/p/pi-)
465 {
466 if(chgm<0) // (n/pi-)
467 {
468 prtM=(-chgm)*mPiCh+brnm*mNeut;
469 prtQC=(-chgm)*PiMiQC+brnm*NeutQC;
470 }
471 else // (n/p)
472 {
473 prtM=chgm*mProt+(brnm-chgm)*mNeut;
474 prtQC=chgm*ProtQC+(brnm-chgm)*NeutQC;
475 }
476 }
477 else if(strm>=brnm) // ---> BigPositiveStrangeness(L/Pi+/K0/K-)
478 {
479 G4int stmb=strm-brnm;
480 if(chgm<0) // (L/K-/K0)
481 {
482 prtM=(-chgm)*mKChg+brnm*mLamb+std::abs(stmb+chgm)*mKZer;
483 prtQC=(-chgm)*KMinQC+brnm*LambQC;
484 if(stmb>-chgm) prtQC+=(stmb+chgm)*KZerQC;
485 else if(stmb<-chgm) prtQC+=(-stmb-chgm)*AKZrQC;
486 }
487 else // (L/K0/pi+)
488 {
489 prtM=chgm*mPiCh+(strm-brnm)*mKZer+brnm*mLamb;
490 prtQC=chgm*PiPlQC+(strm-brnm)*KZerQC+brnm*LambQC;
491 }
492 }
493 else if(strm>0) // ---> PositiveStrangeness<B (L/n/p/Pi+-)
494 {
495 G4int bmst=brnm-strm;
496 if(chgm<0) // (L/n/Pi-)
497 {
498 prtM=(-chgm)*mPiCh+strm*mLamb+bmst*mNeut;
499 prtQC=(-chgm)*PiMiQC+strm*LambQC+bmst*NeutQC;
500 }
501 else if(chgm>=bmst) // (L/p/Pi+)
502 {
503 prtM=(chgm-bmst)*mPiCh+strm*mLamb+bmst*mProt;
504 prtQC=(chgm-bmst)*PiPlQC+strm*LambQC+bmst*ProtQC;
505 }
506 else // ch<bmst (L/p/n)
507 {
508 prtM=chgm*mProt+strm*mLamb+(bmst-chgm)*mNeut;
509 prtQC=chgm*ProtQC+strm*LambQC+(bmst-chgm)*NeutQC;
510 }
511 }
512 else // ---> NegativeStrangeness (N/K+/aK0/Pi-)
513 {
514 G4int bmst=brnm-strm;
515 if(chgm>=bmst) // (K+/p/Pi+)
516 {
517 prtM=(-strm)*mKChg+brnm*mProt+(chgm-bmst)*mPiCh;
518 prtQC=(-strm)*KPlsQC+brnm*ProtQC+(chgm-bmst)*PiPlQC;
519 }
520 else if(chgm>=-strm) // (K+/p/n)
521 {
522 prtM=(-strm)*mKChg+chgm*mProt+(brnm-chgm)*mNeut;
523 prtQC=(-strm)*KPlsQC+chgm*ProtQC+(brnm-chgm)*NeutQC;
524 }
525 else if(chgm>=0) // (K+/aK0/n)
526 {
527 prtM=chgm*mKChg+(-chgm-strm)*mKZer+brnm*mNeut;
528 prtQC=chgm*KPlsQC+(-chgm-strm)*AKZrQC+brnm*NeutQC;
529 }
530 else // ch<0 (aK0/n/Pi-)
531 {
532 prtM=(-strm)*mKChg+(-chgm)*mPiCh+brnm*mNeut;
533 prtQC=(-strm)*KPlsQC+(-chgm)*PiMiQC+brnm*NeutQC;
534 }
535 }
536 EnFlowQC-=prtQC;
537 chgm=chars-chgm; // Just to keep the same notation
538 strm=stras-strm;
539 brnm=curb;
540 if(!strm) // --> The total strangness = 0 (n/p/pi-)
541 {
542 if(chgm<0) resM=(-chgm)*mPiCh+brnm*mNeut;
543 else resM=chgm*mProt+(brnm-chgm)*mNeut;
544 }
545 else if(strm>=brnm) // ---> BigPositiveStrangeness(L/Pi+/K0/K-)
546 {
547 G4int stmb=strm-brnm;
548 if(chgm<0) resM=(-chgm)*mKChg+brnm*mLamb+std::abs(stmb+chgm)*mKZer;
549 else resM=chgm*mPiCh+(strm-brnm)*mKZer+brnm*mLamb;
550 }
551 else if(strm>0) // ---> PositiveStrangeness<B (L/n/p/Pi+-)
552 {
553 G4int bmst=brnm-strm;
554 if (chgm<0) resM=(-chgm)*mPiCh+strm*mLamb+bmst*mNeut;
555 else if(chgm>=bmst) resM=(chgm-bmst)*mPiCh+strm*mLamb+bmst*mProt;
556 else resM=chgm*mProt+strm*mLamb+(bmst-chgm)*mNeut;
557 }
558 else // ---> NegativeStrangeness (N/K+/aK0/Pi-)
559 {
560 G4int bmst=brnm-strm;
561 if (chgm>=bmst) resM=(-strm)*mKChg+brnm*mProt+(chgm-bmst)*mPiCh;
562 else if(chgm>=-strm) resM=(-strm)*mKChg+chgm*mProt+(brnm-chgm)*mNeut;
563 else if(chgm>=0) resM=chgm*mKChg+(-chgm-strm)*mKZer+brnm*mNeut;
564 else resM=(-strm)*mKChg+(-chgm)*mPiCh+brnm*mNeut;
565 }
566 G4LorentzVector prt4M=(prtM/(prtM+resM))*EnFlow4M;
567 EnFlow4M-=prt4M;
568 EnFlowQC-=prtQC;
569 G4QHadron* aHadron = new G4QHadron(prtQC, prt4M);
570 projHV.push_back(aHadron);
571 if(nprt==1)
572 {
573 G4QHadron* fHadron = new G4QHadron(EnFlowQC, EnFlow4M);
574 projHV.push_back(fHadron);
575 nprt=0;
576 }
577#ifdef debug
578 G4cout<<"G4StringChipsParticleLevelInterface::Propagate: nprt="<<nprt<<G4endl;
579#endif
580 } // End of WHILE
581 }
582 else
583 {
584 G4QHadron* aHadron = new G4QHadron(EnFlowQC, EnFlow4M);
585 projHV.push_back(aHadron);
586 }
587
588 // construct the quasmon
589 size_t i;
590 for (i=0; i<theFinalMomenta.size(); i++) delete theFinalMomenta[i];
591 for (i=0; i<theMomenta.size(); i++) delete theMomenta[i];
592 theFinalMomenta.clear();
593 theMomenta.clear();
594
595 G4QNucleus::SetParameters(fractionOfSingleQuasiFreeNucleons,
596 fractionOfPairedQuasiFreeNucleons,
597 clusteringCoefficient,
598 fusionToExchange);
599 G4Quasmon::SetParameters(temperature, halfTheStrangenessOfSee, etaToEtaPrime);
600
601#ifdef CHIPSdebug
602 G4cout<<"G4StringChipsParticleLevelInterface::Propagate: G4QNucleus parameters "
603 <<fractionOfSingleQuasiFreeNucleons<<" "<<fractionOfPairedQuasiFreeNucleons
604 <<" "<<clusteringCoefficient<<G4endl;
605 G4cout<<"G4Quasmon parameters "<<temperature<<" "<<halfTheStrangenessOfSee<<" "
606 <<etaToEtaPrime << G4endl;
607 G4cout<<"The Target PDG code = "<<targetPDGCode<<G4endl;
608 G4cout<<"The projectile momentum = "<<1./MeV*proj4Mom<<G4endl;
609 G4cout<<"The target momentum = "<<1./MeV*targ4Mom<<G4endl;
610#endif
611
612 // now call chips with this info in place
613 G4QHadronVector* output = 0;
614 if (particleCount!=0 && resA!=0)
615 {
616 // G4QCHIPSWorld aWorld(nop); // Create CHIPS World of nop particles
618 G4QEnvironment* pan= new G4QEnvironment(projHV, targetPDGCode);
619#ifdef pdebug
620 G4cout<<"G4StringChipsParticleLevelInterface::Propagate: CHIPS fragmentation, rA="
621 <<resA<<", #AbsPt="<<particleCount<<G4endl;
622#endif
623 try
624 {
625 output = pan->Fragment(); // The main fragmentation member function
626 }
627 catch(G4HadronicException& aR)
628 {
629 G4cerr << "Exception thrown of G4StringChipsParticleLevelInterface "<<G4endl;
630 G4cerr << " targetPDGCode = "<< targetPDGCode <<G4endl;
631 G4cerr << " The projectile momentum = "<<1./MeV*proj4Mom<<G4endl;
632 G4cerr << " The target momentum = "<<1./MeV*targ4Mom<<G4endl<<G4endl;
633 G4cerr << " Dumping the information in the pojectile list"<<G4endl;
634 for(i=0; i< projHV.size(); i++)
635 {
636 G4cerr <<" Incoming 4-momentum and PDG code of "<<i<<"'th hadron: "
637 <<" "<< projHV[i]->Get4Momentum()<<" "<<projHV[i]->GetPDGCode()<<G4endl;
638 }
639 throw;
640 }
641 delete pan;
642 }
643 else
644 {
645#ifdef pdebug
646 G4cout<<"G4StringChipsParticleLevelInterface::Propagate: NO CHIPS fragmentation, rA="
647 <<resA<<", #AbsPt="<<particleCount<<G4endl;
648#endif
649 output = new G4QHadronVector;
650 }
651
652 // clean up impinging particles
653 std::for_each(projHV.begin(), projHV.end(), DeleteQHadron());
654 projHV.clear();
655
656 // Fill the result.
657#ifdef CHIPSdebug
658 G4cout << "NEXT EVENT, EscapeExists="<<EscapeExists<<G4endl;
659#endif
660
661 // first decay and add all escaping particles.
662 if (EscapeExists) for (current = firstEscape; current!=theSorted.end(); current++)
663 {
664 G4KineticTrack* aResult = (*current).second;
665 G4ParticleDefinition* pdef=aResult->GetDefinition();
666 secondaries = NULL;
667 //if(pdef->GetPDGWidth() > 0 && pdef->GetPDGLifeTime() < 5E-17*s ) // HPW version
668 if ( pdef->IsShortLived() )
669 {
670 secondaries = aResult->Decay();
671 for (unsigned int aSecondary=0; aSecondary<secondaries->size(); aSecondary++)
672 {
673 G4KineticTrack* bResult=secondaries->operator[](aSecondary);
674 G4ParticleDefinition* sdef=bResult->GetDefinition();
675 if ( sdef->IsShortLived() )
676 {
677 secsec = bResult->Decay();
678 for (unsigned int bSecondary=0; bSecondary<secsec->size(); bSecondary++)
679 {
680 G4KineticTrack* cResult=secsec->operator[](bSecondary);
681 G4ParticleDefinition* cdef=cResult->GetDefinition();
682 theSec = new G4ReactionProduct(cdef);
683 G4LorentzVector cur4Mom = cResult->Get4Momentum();
684 cur4Mom.boost(targ4Mom.boostVector());
685 theSec->SetTotalEnergy(cur4Mom.t());
686 theSec->SetMomentum(cur4Mom.vect());
687#ifdef trapdebug
688 if(cdef->GetPDGEncoding()==113) G4cout
689 <<"G4StringChipsParticleLevelInterface::Propagate: *Rho0* QGS dec2 PDG="
690 <<cdef->GetPDGEncoding()<<",4M="<<cur4Mom<<", grandparPDG= "
691 <<pdef->GetPDGEncoding()<<", parPDG= "<<sdef->GetPDGEncoding()<<G4endl;
692#endif
693#ifdef pdebug
694 G4cout<<"G4StringChipsParticleLevelInterface::Propagate: *OUT* QGS dec2 PDG="
695 <<sdef->GetPDGEncoding()<<",4M="<<cur4Mom<<G4endl;
696#endif
697 theResult->push_back(theSec);
698 }
699 std::for_each(secsec->begin(), secsec->end(), DeleteKineticTrack());
700 delete secsec;
701 }
702 else
703 {
704 theSec = new G4ReactionProduct(sdef);
705 G4LorentzVector current4Mom = bResult->Get4Momentum();
706 current4Mom.boost(targ4Mom.boostVector());
707 theSec->SetTotalEnergy(current4Mom.t());
708 theSec->SetMomentum(current4Mom.vect());
709#ifdef trapdebug
710 if(sdef->GetPDGEncoding()==113)
711 G4cout<<"G4StringChipsParticleLevelInterface::Propagate:*Rho0* QGS decay PDG="
712 <<sdef->GetPDGEncoding()<<",4M="<<current4Mom<<", parentPDG= "
713 <<pdef->GetPDGEncoding()<<G4endl;
714 //throw G4HadronicException(__FILE__,__LINE__,
715 // "G4StringChipsParticleLevelInterface: Rho0 is found!");
716#endif
717#ifdef pdebug
718 G4cout<<"G4StringChipsParticleLevelInterface::Propagate: *OUT* QGS decay PDG="
719 <<sdef->GetPDGEncoding()<<",4M="<<current4Mom<<G4endl;
720#endif
721 theResult->push_back(theSec);
722 }
723 }
724 std::for_each(secondaries->begin(), secondaries->end(), DeleteKineticTrack());
725 delete secondaries;
726 }
727 else
728 {
729 theSec = new G4ReactionProduct(aResult->GetDefinition());
730 G4LorentzVector current4Mom = aResult->Get4Momentum();
731 current4Mom.boost(targ4Mom.boostVector());
732 theSec->SetTotalEnergy(current4Mom.t());
733 theSec->SetMomentum(current4Mom.vect());
734#ifdef trapdebug
735 if(aResult->GetDefinition()->GetPDGEncoding()==113)
736 G4cout<<"G4StringChipsParticleLevelInterface::Propagate: *OUT* QGS stable PDG="
737 <<aResult->GetDefinition()->GetPDGEncoding()<<",4M="<<current4Mom<<G4endl;
738#endif
739#ifdef pdebug
740 G4cout<<"G4StringChipsParticleLevelInterface::Propagate: *OUT* QGS stable PDG="
741 <<aResult->GetDefinition()->GetPDGEncoding()<<",4M="<<current4Mom<<G4endl;
742#endif
743 theResult->push_back(theSec);
744 }
745 }
746 std::for_each(theSecondaries->begin(), theSecondaries->end(), DeleteKineticTrack());
747 delete theSecondaries;
748
749 // now add the quasmon output
750 G4int maxParticle=output->size();
751#ifdef CHIPSdebug
752 G4cout << "Number of particles from string"<<theResult->size()<<G4endl;
753 G4cout << "Number of particles from chips"<<maxParticle<<G4endl;
754#endif
755#ifdef pdebug
756 G4cout << "Number of particles from QGS="<<theResult->size()<<G4endl;
757 G4cout << "Number of particles from CHIPS="<<maxParticle<<G4endl;
758#endif
759 if(maxParticle) for(G4int particle = 0; particle < maxParticle; particle++)
760 {
761 if(output->operator[](particle)->GetNFragments() != 0)
762 {
763 delete output->operator[](particle);
764 continue;
765 }
766 G4int pdgCode = output->operator[](particle)->GetPDGCode();
767
768
769#ifdef CHIPSdebug
770 G4cerr << "PDG code of chips particle = "<<pdgCode<<G4endl;
771#endif
772
773 G4ParticleDefinition * theDefinition;
774 // Note that I still have to take care of strange nuclei
775 // For this I need the mass calculation, and a changed interface
776 // for ion-table ==> work for Hisaya @@@@@@@
777 // Then I can sort out the pdgCode. I also need a decau process
778 // for strange nuclei; may be another chips interface
779 if(pdgCode>90000000)
780 {
781 G4int aZ = (pdgCode-90000000)/1000;
782 if (aZ>1000) aZ=aZ%1000; // patch for strange nuclei, to be repaired @@@@
783 G4int anN = pdgCode-90000000-1000*aZ;
784 if(anN>1000) anN=anN%1000; // patch for strange nuclei, to be repaired @@@@
785
786 if(pdgCode==90000999) theDefinition = G4PionPlus::PionPlusDefinition();
787 else if(pdgCode==89999001) theDefinition = G4PionMinus::PionMinusDefinition();
788 else if(pdgCode==90999999) theDefinition = G4KaonZero::KaonZeroDefinition();
789 else if(pdgCode==90999000) theDefinition = G4KaonMinus::KaonMinusDefinition();
790 else if(pdgCode==89001000) theDefinition = G4KaonPlus::KaonPlusDefinition();
791 else if(pdgCode==89000001) theDefinition = G4AntiKaonZero::AntiKaonZeroDefinition();
792 else if(pdgCode==91000000) theDefinition = G4Lambda::LambdaDefinition();
793 else if(pdgCode==92000000) theDefinition = G4Lambda::LambdaDefinition(); //NLambd?
794 else if(pdgCode==93000000) theDefinition = G4Lambda::LambdaDefinition();
795 else if(pdgCode==94000000) theDefinition = G4Lambda::LambdaDefinition();
796 else if(pdgCode==95000000) theDefinition = G4Lambda::LambdaDefinition();
797 else if(pdgCode==96000000) theDefinition = G4Lambda::LambdaDefinition();
798 else if(pdgCode==97000000) theDefinition = G4Lambda::LambdaDefinition();
799 else if(pdgCode==98000000) theDefinition = G4Lambda::LambdaDefinition();
800 else if(aZ == 0 && anN == 1) theDefinition = G4Neutron::Neutron();
801 else theDefinition = G4ParticleTable::GetParticleTable()->FindIon(aZ,anN+aZ,0,aZ);
802 }
803 else theDefinition = G4ParticleTable::GetParticleTable()->FindParticle(pdgCode);
804
805#ifdef CHIPSdebug
806 G4cout << "Particle code produced = "<< pdgCode <<G4endl;
807#endif
808
809 if(theDefinition)
810 {
811 theSec = new G4ReactionProduct(theDefinition);
812 G4LorentzVector current4Mom = output->operator[](particle)->Get4Momentum();
813 current4Mom.boost(targ4Mom.boostVector());
814 theSec->SetTotalEnergy(current4Mom.t());
815 theSec->SetMomentum(current4Mom.vect());
816#ifdef pdebug
817 G4cout<<"G4StringChipsParticleLevelInterface::Propagate: *OUT* CHIPS PDG="
818 <<theDefinition->GetPDGEncoding()<<",4M="<<current4Mom<<G4endl;
819#endif
820 theResult->push_back(theSec);
821 }
822#ifdef pdebug
823 else
824 {
825 G4cerr <<"G4StringChipsParticleLevelInterface::Propagate: WARNING"<<G4endl;
826 G4cerr << "Getting unknown pdgCode from chips in ParticleLevelInterface"<<G4endl;
827 G4cerr << "skipping particle with pdgCode = "<<pdgCode<<G4endl<<G4endl;
828 }
829#endif
830
831#ifdef CHIPSdebug
832 G4cout <<"CHIPS particles "<<theDefinition->GetPDGCharge()<<" "
833 << theDefinition->GetPDGEncoding()<<" "
834 << output->operator[](particle)->Get4Momentum() <<G4endl;
835#endif
836
837 delete output->operator[](particle);
838 }
839 else
840 {
841 if(resA>0)
842 {
843 G4ParticleDefinition* theDefinition = G4Neutron::Neutron();
844 if(resA==1) // The residual nucleus at rest must be added to conserve BaryN & Charge
845 {
846 if(resZ == 1) theDefinition = G4Proton::Proton();
847 }
848 else theDefinition = G4ParticleTable::GetParticleTable()->FindIon(resZ,resA,0,resZ);
849 theSec = new G4ReactionProduct(theDefinition);
850 theSec->SetTotalEnergy(theDefinition->GetPDGMass());
851 theSec->SetMomentum(G4ThreeVector(0.,0.,0.));
852 theResult->push_back(theSec);
853 if(!resZ && resA>0) for(G4int ni=1; ni<resA; ni++)
854 {
855 theSec = new G4ReactionProduct(theDefinition);
856 theSec->SetTotalEnergy(theDefinition->GetPDGMass());
857 theSec->SetMomentum(G4ThreeVector(0.,0.,0.));
858 theResult->push_back(theSec);
859 }
860 }
861 }
862 delete output;
863
864#ifdef CHIPSdebug
865 G4cout << "Number of particles"<<theResult->size()<<G4endl;
866 G4cout << G4endl;
867 G4cout << "QUASMON preparation info "
868 << 1./MeV*proj4Mom<<" "
869 << 1./MeV*targ4Mom<<" "
870 << nD<<" "<<nU<<" "<<nS<<" "<<nAD<<" "<<nAU<<" "<<nAS<<" "
871 << hitCount<<" "
872 << particleCount<<" "
873 << theLow<<" "
874 << theHigh<<" "
875 << G4endl;
876#endif
877
878 return theResult;
879}
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
#define G4endl
Definition: G4ios.hh:52
#define G4cin
Definition: G4ios.hh:51
G4DLLIMPORT std::ostream G4cerr
G4DLLIMPORT std::ostream G4cout
double mag2() const
Hep3Vector boostVector() const
HepLorentzVector & boost(double, double, double)
Hep3Vector vect() const
static G4AntiKaonZero * AntiKaonZeroDefinition()
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &aTargetNucleus, G4HadFinalState *aChange=0)
static G4KaonMinus * KaonMinusDefinition()
Definition: G4KaonMinus.cc:108
static G4KaonPlus * KaonPlus()
Definition: G4KaonPlus.cc:113
static G4KaonPlus * KaonPlusDefinition()
Definition: G4KaonPlus.cc:108
static G4KaonZero * KaonZero()
Definition: G4KaonZero.cc:104
static G4KaonZero * KaonZeroDefinition()
Definition: G4KaonZero.cc:99
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 GetPDGCharge() const
G4ParticleDefinition * FindIon(G4int atomicNumber, G4int atomicMass, G4double excitationEnergy)
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
static G4ParticleTable * GetParticleTable()
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
static G4PionMinus * PionMinusDefinition()
Definition: G4PionMinus.cc:93
static G4PionPlus * PionPlusDefinition()
Definition: G4PionPlus.cc:93
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4QCHIPSWorld * Get()
G4QParticleVector * GetParticles(G4int nOfParts=0)
G4int GetCharge() const
Definition: G4QContent.cc:1159
G4int GetBaryonNumber() const
Definition: G4QContent.cc:1182
G4int GetStrangeness() const
Definition: G4QContent.hh:184
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)
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &theNucleus)
virtual G4ReactionProductVector * Propagate(G4KineticTrackVector *theSecondaries, G4V3DNucleus *theNucleus)
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