Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
39 #include "G4StringChipsParticleLevelInterface.hh"
40 #include "globals.hh"
41 #include "G4SystemOfUnits.hh"
42 #include "G4KineticTrackVector.hh"
43 #include "G4Nucleon.hh"
44 #include "G4Proton.hh"
45 #include "G4Neutron.hh"
46 #include "G4LorentzRotation.hh"
47 #include "G4HadronicException.hh"
48 // #define CHIPSdebug
49 // #define CHIPSdebug_1
50 
51 #ifdef hdebug_SCPLI
52 const G4int G4StringChipsParticleLevelInterface::nbh=200;
53 G4double G4StringChipsParticleLevelInterface::bhmax=20.;
54 G4double G4StringChipsParticleLevelInterface::ehmax=20.;
55 G4double G4StringChipsParticleLevelInterface::bhdb=0.;
56 G4double G4StringChipsParticleLevelInterface::ehde=0.;
57 G4double G4StringChipsParticleLevelInterface::toth=0.;
58 G4int G4StringChipsParticleLevelInterface::bover=0;
59 G4int G4StringChipsParticleLevelInterface::eover=0;
60 G4int* G4StringChipsParticleLevelInterface::bhis =
61  new G4int[G4StringChipsParticleLevelInterface::nbh];
62 G4int* G4StringChipsParticleLevelInterface::ehis =
63  new G4int[G4StringChipsParticleLevelInterface::nbh];
64 #endif
65 
66 G4StringChipsParticleLevelInterface::G4StringChipsParticleLevelInterface()
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 
110 G4HadFinalState* G4StringChipsParticleLevelInterface::
111 ApplyYourself(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 
119 G4ReactionProductVector* G4StringChipsParticleLevelInterface::
120 Propagate(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  {
145  G4ReactionProductVector* theFastResult = new G4ReactionProductVector;
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;
283  G4ReactionProductVector* theResult = new G4ReactionProductVector;
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
617  G4QCHIPSWorld::Get()->GetParticles(nop);
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 }
Generated on Sat May 25 2013 14:33:50 for Geant4 by   doxygen 1.8.4