source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/interface/src/G4QStringChipsParticleLevelInterface.cc @ 1197

Last change on this file since 1197 was 1196, checked in by garnier, 15 years ago

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