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

Last change on this file since 968 was 819, checked in by garnier, 16 years ago

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