Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4StringChipsInterface.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 // !!! Was used in QBBC PL, NOW it is not. Must be absolete !!!
27 // =----------------------------------------------------------=
28 
29 //#define debug
30 //#define pdebug
31 
32 #include <utility>
33 #include <list>
34 
35 #include "G4StringChipsInterface.hh"
36 #include "globals.hh"
37 #include "G4SystemOfUnits.hh"
38 #include "G4KineticTrackVector.hh"
39 #include "G4Nucleon.hh"
40 #include "G4LorentzRotation.hh"
41 
42 G4StringChipsInterface::G4StringChipsInterface()
43 {
44 #ifdef CHIPSdebug
45  G4cout << "Please enter the energy loss per fermi in GeV"<<G4endl;
46  G4cin >> theEnergyLossPerFermi;
47 #endif
48 #ifdef debug
49  G4cout<<"G4StringChipsInterface::Constructor is called"<<G4endl;
50 #endif
51  theEnergyLossPerFermi = 0.5*GeV;
52  // theEnergyLossPerFermi = 1.*GeV;
53 }
54 
55 G4HadFinalState* G4StringChipsInterface::
56 ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& theNucleus)
57 {
58 #ifdef debug
59  G4cout<<"G4StringChipsInterface::ApplyYourself is called"<<G4endl;
60 #endif
61  return theModel.ApplyYourself(aTrack, theNucleus);
62 }
63 
64 G4ReactionProductVector* G4StringChipsInterface::
65 Propagate(G4KineticTrackVector* theSecondaries, G4V3DNucleus* theNucleus)
66 {
67 #ifdef debug
68  G4cout<<"G4StringChipsInterface::Propagate is called"<<G4endl;
69 #endif
70  // Protection for non physical conditions
71 
72  if(theSecondaries->size() == 1)
73  throw G4HadronicException(__FILE__, __LINE__, "G4StringChipsInterface: Only one particle from String models!");
74 
75  // Calculate the mean energy lost
76  std::pair<G4double, G4double> theImpact = theNucleus->RefetchImpactXandY();
77  G4double impactX = theImpact.first;
78  G4double impactY = theImpact.second;
79  G4double inpactPar2 = impactX*impactX + impactY*impactY;
80 
81  G4double radius2 = theNucleus->GetNuclearRadius(5*perCent);
82  radius2 *= radius2;
83  G4double pathlength = 0;
84  if(radius2 - inpactPar2>0) pathlength = 2.*std::sqrt(radius2 - inpactPar2);
85  G4double theEnergyLostInFragmentation = theEnergyLossPerFermi*pathlength/fermi;
86 
87  // now select all particles in range
88  std::list<std::pair<G4double, G4KineticTrack *> > theSorted;
89  std::list<std::pair<G4double, G4KineticTrack *> >::iterator current;
90  for(unsigned int secondary = 0; secondary<theSecondaries->size(); secondary++)
91  {
92  G4LorentzVector a4Mom = theSecondaries->operator[](secondary)->Get4Momentum();
93 
94 #ifdef CHIPSdebug
95  G4cout <<"ALL STRING particles "<<theSecondaries->operator[](secondary)->GetDefinition()->GetPDGCharge()<<" "
96  << theSecondaries->operator[](secondary)->GetDefinition()->GetPDGEncoding()<<" "
97  << a4Mom <<G4endl;
98 #endif
99 
100  G4double toSort = a4Mom.rapidity();
101  std::pair<G4double, G4KineticTrack *> it;
102  it.first = toSort;
103  it.second = theSecondaries->operator[](secondary);
104  G4bool inserted = false;
105  for(current = theSorted.begin(); current!=theSorted.end(); current++)
106  {
107  if((*current).first > toSort)
108  {
109  theSorted.insert(current, it);
110  inserted = true;
111  break;
112  }
113  }
114  if(!inserted)
115  {
116  theSorted.push_front(it);
117  }
118  }
119 
120  G4LorentzVector proj4Mom(0.,0.,0.,0.);
121  // @@ Use the G4QContent class, which is exactly (nD,nU,nS,nAD,nAU,nAS) !
122  // The G4QContent class is a basic clas in CHIPS (not PDG Code as in GEANT4),
123  // so in CHIPS on can has a hadronic obgect (Quasmon) with any Quark Content.
124  // As a simple extantion for the hadron (which is a special case for Hadron)
125  // there is a clas G4QChipolino, which is a Quasmon, which can decay in two
126  // hadrons. In future the three-hadron class can be added...
127  G4int nD = 0;
128  G4int nU = 0;
129  G4int nS = 0;
130  G4int nAD = 0;
131  G4int nAU = 0;
132  G4int nAS = 0;
133  std::list<std::pair<G4double, G4KineticTrack *> >::iterator firstEscaping = theSorted.begin();
134  G4double runningEnergy = 0;
135  G4int particleCount = 0;
136  G4LorentzVector theLow = (*(theSorted.begin())).second->Get4Momentum();
137  G4LorentzVector theHigh;
138 
139 #ifdef CHIPSdebug
140  G4cout << "CHIPS ENERGY LOST "<<theEnergyLostInFragmentation<<G4endl;
141  G4cout << "sorted rapidities event start"<<G4endl;
142 #endif
143 
144  G4ReactionProductVector * theResult = new G4ReactionProductVector;
145  G4ReactionProduct * theSec;
146  G4KineticTrackVector * secondaries;
147 #ifdef pdebug
148  G4cout<<"G4StringChipsInterface::Propagate: Absorption"<<G4endl;
149 #endif
150 
151  // first decay and add all escaping particles.
152  for(current = theSorted.begin(); current!=theSorted.end(); current++)
153  {
154  firstEscaping = current;
155  if((*current).second->GetDefinition()->GetQuarkContent(3)!=0 ||
156  (*current).second->GetDefinition()->GetAntiQuarkContent(3) !=0)
157  {
158  G4KineticTrack * aResult = (*current).second;
159  G4ParticleDefinition * pdef=aResult->GetDefinition();
160  secondaries = NULL;
161  if ( pdef->GetPDGWidth() > 0 && pdef->GetPDGLifeTime() < 5E-17*s )
162  {
163  secondaries = aResult->Decay();
164  }
165  if ( secondaries == NULL )
166  {
167  theSec = new G4ReactionProduct(aResult->GetDefinition());
168  G4LorentzVector current4Mom = aResult->Get4Momentum();
169  theSec->SetTotalEnergy(current4Mom.t());
170  theSec->SetMomentum(current4Mom.vect());
171  theResult->push_back(theSec);
172  }
173  else
174  {
175  for (unsigned int aSecondary=0; aSecondary<secondaries->size(); aSecondary++)
176  {
177  theSec = new G4ReactionProduct(secondaries->operator[](aSecondary)->GetDefinition());
178  G4LorentzVector current4Mom = secondaries->operator[](aSecondary)->Get4Momentum();
179  theSec->SetTotalEnergy(current4Mom.t());
180  theSec->SetMomentum(current4Mom.vect());
181  theResult->push_back(theSec);
182  }
183  std::for_each(secondaries->begin(), secondaries->end(), DeleteKineticTrack());
184  delete secondaries;
185  }
186  continue;
187  }
188  runningEnergy += (*current).second->Get4Momentum().t();
189  if(runningEnergy > theEnergyLostInFragmentation) break;
190 
191 #ifdef CHIPSdebug
192  G4cout <<"ABSORBED STRING particles "<<current->second->GetDefinition()->GetPDGCharge()<<" "
193  << current->second->GetDefinition()->GetPDGEncoding()<<" "
194  << current->second->Get4Momentum() <<G4endl;
195 #endif
196 #ifdef pdebug
197  G4cout<<"G4StringChipsInterface::Propagate:C="
198  <<current->second->GetDefinition()->GetPDGCharge()<<", PDG="
199  << current->second->GetDefinition()->GetPDGEncoding()<<", 4M="
200  << current->second->Get4Momentum() <<G4endl;
201 #endif
202 
203  // projectile 4-momentum needed in constructor of quasmon
204  particleCount++;
205  theHigh = (*current).second->Get4Momentum();
206  proj4Mom += (*current).second->Get4Momentum();
207 
208 #ifdef CHIPSdebug
209  G4cout << "sorted rapidities "<<current->second->Get4Momentum().rapidity()<<G4endl;
210 #endif
211 
212  // projectile quark contents needed for G4QContent construction (@@ ? -> G4QContent class)
213  nD += (*current).second->GetDefinition()->GetQuarkContent(1);
214  nU += (*current).second->GetDefinition()->GetQuarkContent(2);
215  nS += (*current).second->GetDefinition()->GetQuarkContent(3);
216  nAD += (*current).second->GetDefinition()->GetAntiQuarkContent(1);
217  nAU += (*current).second->GetDefinition()->GetAntiQuarkContent(2);
218  nAS += (*current).second->GetDefinition()->GetAntiQuarkContent(3);
219  }
220  // construct G4QContent
221 
222 #ifdef CHIPSdebug
223  G4cout << "Quark content: d="<<nD<<", u="<<nU<< ", s="<< nS << G4endl;
224  G4cout << "Anti-quark content: anit-d="<<nAD<<", anti-u="<<nAU<< ", anti-s="<< nAS << G4endl;
225 #endif
226 
227  G4QContent theProjectiles(nD, nU, nS, nAD, nAU, nAS);
228 
229 #ifdef CHIPSdebug
230  G4cout << "G4QContent is constructed"<<endl;
231 #endif
232 
233  // target properties needed in constructor of quasmon
234  // remove all hit nucleons to get Target code
235  theNucleus->StartLoop();
236  G4Nucleon * aNucleon;
237  G4int resA = 0;
238  G4int resZ = 0;
239  G4ThreeVector hitMomentum(0,0,0);
240  G4double hitMass = 0;
241  G4int hitCount = 0;
242  while((aNucleon = theNucleus->GetNextNucleon()))
243  {
244  if(!aNucleon->AreYouHit())
245  {
246  resA++;
247  resZ+=G4int (aNucleon->GetDefinition()->GetPDGCharge());
248  }
249  else
250  {
251  hitMomentum += aNucleon->GetMomentum().vect();
252  hitMass += aNucleon->GetMomentum().m();
253  hitCount ++;
254  }
255  }
256  G4int targetPDGCode = 90000000 + 1000*resZ + (resA-resZ);
257  G4double targetMass = theNucleus->GetMass();
258  targetMass -= hitMass;
259  G4double targetEnergy = std::sqrt(hitMomentum.mag2()+targetMass*targetMass);
260  // !! @@ Target should be at rest: hitMomentum=(0,0,0) @@ !! M.K.
261  G4LorentzVector targ4Mom(-1.*hitMomentum, targetEnergy);
262 
263  // construct the quasmon
264  G4int nop = 85; // clusters up to Alpha cluster (Reduced)
265  // G4int nop = 122; // clusters up to Alpha cluster
266  // G4int nop = 152; // not reduced upto Li6
267  G4double fractionOfSingleQuasiFreeNucleons = 0.5; // It is A-dependent (C=.85, U=.40)
268  G4double fractionOfPairedQuasiFreeNucleons = 0.05;
269  G4double clusteringCoefficient = 5.;
270  G4double temperature = 180.;
271  G4double halfTheStrangenessOfSee = 0.3; // = s/d = s/u
272  G4double etaToEtaPrime = 0.3;
273 
274  G4QNucleus::SetParameters(fractionOfSingleQuasiFreeNucleons,
275  fractionOfPairedQuasiFreeNucleons,
276  clusteringCoefficient);
277  G4Quasmon::SetParameters(temperature,
278  halfTheStrangenessOfSee,
279  etaToEtaPrime);
280 
281 #ifdef CHIPSdebug
282  G4cout << "G4QNucleus parameters "<< fractionOfSingleQuasiFreeNucleons << " "
283  << fractionOfPairedQuasiFreeNucleons << " "<< clusteringCoefficient << G4endl;
284  G4cout << "G4Quasmon parameters "<< temperature << " "<< halfTheStrangenessOfSee << " "
285  <<etaToEtaPrime << G4endl;
286  G4cout << "The Target PDG code = "<<targetPDGCode<<G4endl;
287  G4cout << "The projectile momentum = "<<1./MeV*proj4Mom<<G4endl;
288  G4cout << "The target momentum = "<<1./MeV*targ4Mom<<G4endl;
289 #endif
290 
291 
292  // Chips expects all in target rest frame, along z.
293  // G4QCHIPSWorld aWorld(nop); // Create CHIPS World of nop particles
294  G4QCHIPSWorld::Get()->GetParticles(nop);
295  G4QHadronVector projHV;
296  // target rest frame
297  proj4Mom.boost(-1.*targ4Mom.boostVector());
298  // now go along z
299  G4LorentzRotation toZ;
300  toZ.rotateZ(-1*proj4Mom.phi());
301  toZ.rotateY(-1*proj4Mom.theta());
302  proj4Mom = toZ*proj4Mom;
303  G4LorentzRotation toLab(toZ.inverse());
304 
305 #ifdef CHIPSdebug
306  G4cout << "a Boosted projectile vector along z"<<proj4Mom<<" "<<proj4Mom.mag()<<G4endl;
307 #endif
308 
309  G4QHadron* iH = new G4QHadron(theProjectiles, 1./MeV*proj4Mom);
310  projHV.push_back(iH);
311 
312  // now call chips with this info in place
313  G4QHadronVector * output = 0;
314  if (particleCount!=0)
315  {
316  G4QEnvironment* pan= new G4QEnvironment(projHV, targetPDGCode);
317  try
318  {
319  output = pan->Fragment();
320  }
321  catch(G4HadronicException & aR)
322  {
323  G4cerr << "Exception thrown passing through G4ChiralInvariantPhaseSpace "<<G4endl;
324  G4cerr << " targetPDGCode = "<< targetPDGCode <<G4endl;
325  G4cerr << " Dumping the information in the pojectile list"<<G4endl;
326  for(size_t i=0; i< projHV.size(); i++)
327  {
328  G4cerr <<" Incoming 4-momentum and PDG code of "<<i<<"'th hadron: "
329  <<" "<< projHV[i]->Get4Momentum()<<" "<<projHV[i]->GetPDGCode()<<G4endl;
330  }
331  throw;
332  }
333  std::for_each(projHV.begin(), projHV.end(), DeleteQHadron());
334  projHV.clear();
335  delete pan;
336  }
337  else output = new G4QHadronVector;
338 
339  // Fill the result.
340 #ifdef CHIPSdebug
341  G4cout << "NEXT EVENT"<<endl;
342 #endif
343 
344  // first decay and add all escaping particles.
345  for(current = firstEscaping; current!=theSorted.end(); current++)
346  {
347  G4KineticTrack * aResult = (*current).second;
348  G4ParticleDefinition * pdef=aResult->GetDefinition();
349  secondaries = NULL;
350  if ( pdef->GetPDGWidth() > 0 && pdef->GetPDGLifeTime() < 5E-17*s )
351  {
352  secondaries = aResult->Decay();
353  }
354  if ( secondaries == NULL )
355  {
356  theSec = new G4ReactionProduct(aResult->GetDefinition());
357  G4LorentzVector current4Mom = aResult->Get4Momentum();
358  current4Mom = toLab*current4Mom;
359  current4Mom.boost(targ4Mom.boostVector());
360  theSec->SetTotalEnergy(current4Mom.t());
361  theSec->SetMomentum(current4Mom.vect());
362  theResult->push_back(theSec);
363  }
364  else
365  {
366  for (unsigned int aSecondary=0; aSecondary<secondaries->size(); aSecondary++)
367  {
368  theSec = new G4ReactionProduct(secondaries->operator[](aSecondary)->GetDefinition());
369  G4LorentzVector current4Mom = secondaries->operator[](aSecondary)->Get4Momentum();
370  current4Mom = toLab*current4Mom;
371  current4Mom.boost(targ4Mom.boostVector());
372  theSec->SetTotalEnergy(current4Mom.t());
373  theSec->SetMomentum(current4Mom.vect());
374  theResult->push_back(theSec);
375  }
376  std::for_each(secondaries->begin(), secondaries->end(), DeleteKineticTrack());
377  delete secondaries;
378  }
379  }
380  std::for_each(theSecondaries->begin(), theSecondaries->end(), DeleteKineticTrack());
381  delete theSecondaries;
382 
383  // now add the quasmon output
384 #ifdef CHIPSdebug
385  G4cout << "Number of particles from string"<<theResult->size()<<G4endl;
386  G4cout << "Number of particles from chips"<<output->size()<<G4endl;
387 #endif
388 
389  for(unsigned int particle = 0; particle < output->size(); particle++)
390  {
391  if(output->operator[](particle)->GetNFragments() != 0)
392  {
393  delete output->operator[](particle);
394  continue;
395  }
396  // theSec = new G4ReactionProduct; // JA - not used, and memory leaked (Coverity)
397  G4int pdgCode = output->operator[](particle)->GetPDGCode();
398  G4ParticleDefinition * theDefinition;
399  // Note that I still have to take care of strange nuclei
400  // For this I need the mass calculation, and a changed interface
401  // for ion-table ==> work for Hisaya @@@@@@@
402  // Then I can sort out the pdgCode. I also need a decau process
403  // for strange nuclei; may be another chips interface
404  if(pdgCode>90000000)
405  {
406  G4int aZ = (pdgCode-90000000)/1000;
407  if (aZ>1000) aZ=aZ%1000; // patch for strange nuclei, to be repaired @@@@
408  G4int anN = pdgCode-90000000-1000*aZ;
409  if(anN>1000) anN=anN%1000; // patch for strange nuclei, to be repaired @@@@
410  if(pdgCode==91000000) theDefinition = G4Lambda::LambdaDefinition();
411  else if(pdgCode==92000000) theDefinition = G4Lambda::LambdaDefinition();
412  else if(pdgCode==93000000) theDefinition = G4Lambda::LambdaDefinition();
413  else if(pdgCode==94000000) theDefinition = G4Lambda::LambdaDefinition();
414  else if(pdgCode==95000000) theDefinition = G4Lambda::LambdaDefinition();
415  else if(pdgCode==96000000) theDefinition = G4Lambda::LambdaDefinition();
416  else if(pdgCode==97000000) theDefinition = G4Lambda::LambdaDefinition();
417  else if(pdgCode==98000000) theDefinition = G4Lambda::LambdaDefinition();
418  else if(aZ == 0 && anN == 1) theDefinition = G4Neutron::Neutron();
419  else theDefinition = G4ParticleTable::GetParticleTable()->FindIon(aZ,anN+aZ,0,aZ);
420  }
421  else
422  {
423  theDefinition = G4ParticleTable::GetParticleTable()->FindParticle(output->operator[](particle)->GetPDGCode());
424  }
425 #ifdef CHIPSdebug
426  G4cout << "Particle code produced = "<< pdgCode <<G4endl;
427 #endif
428 
429  theSec = new G4ReactionProduct(theDefinition);
430  G4LorentzVector current4Mom = output->operator[](particle)->Get4Momentum();
431  current4Mom = toLab*current4Mom;
432  current4Mom.boost(targ4Mom.boostVector());
433  theSec->SetTotalEnergy(current4Mom.t());
434  theSec->SetMomentum(current4Mom.vect());
435  theResult->push_back(theSec);
436 
437 #ifdef CHIPSdebug
438  G4cout <<"CHIPS particles "<<theDefinition->GetPDGCharge()<<" "
439  << theDefinition->GetPDGEncoding()<<" "
440  << current4Mom <<G4endl;
441 #endif
442 
443  delete output->operator[](particle);
444  }
445  delete output;
446 #ifdef CHIPSdebug
447  G4cout << "Number of particles"<<theResult->size()<<G4endl;
448  G4cout << G4endl;
449  // @@ G4QContent has even the out option!
450  G4cout << "QUASMON preparation info "
451  << 1./MeV*proj4Mom<<" "
452  << 1./MeV*targ4Mom<<" "
453  << nD<<" "<<nU<<" "<<nS<<" "<<nAD<<" "<<nAU<<" "<<nAS<<" "
454  << hitCount<<" "
455  << particleCount<<" "
456  << theLow<<" "
457  << theHigh<<" "
458  << G4endl;
459 #endif
460 
461  return theResult;
462 }
Generated on Sat May 25 2013 14:33:50 for Geant4 by   doxygen 1.8.4