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