source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/interface/src/G4StringChipsParticleLevelInterface.cc

Last change on this file was 1228, checked in by garnier, 14 years ago

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