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