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

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

import all except CVS
File size: 17.1 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 pdebug
29
30#include "G4StringChipsInterface.hh"
31#include "globals.hh"
32#include <utility>
33#include <list>
34#include "G4KineticTrackVector.hh"
35#include "G4Nucleon.hh"
36#include "G4LorentzRotation.hh"
37
38G4StringChipsInterface::G4StringChipsInterface()
39{
40#ifdef CHIPSdebug
41  G4cout << "Please enter the energy loss per fermi in GeV"<<G4endl;
42  G4cin >> theEnergyLossPerFermi;
43#endif
44#ifdef debug
45  G4cout<<"G4StringChipsInterface::Constructor is called"<<G4endl;
46#endif
47  theEnergyLossPerFermi = 0.5*GeV;
48  // theEnergyLossPerFermi = 1.*GeV;
49}
50
51G4HadFinalState* G4StringChipsInterface::
52ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& theNucleus)
53{
54#ifdef debug
55  G4cout<<"G4StringChipsInterface::ApplyYourself is called"<<G4endl;
56#endif
57  return theModel.ApplyYourself(aTrack, theNucleus);
58}
59
60G4ReactionProductVector* G4StringChipsInterface::
61Propagate(G4KineticTrackVector* theSecondaries, G4V3DNucleus* theNucleus)
62{
63#ifdef debug
64  G4cout<<"G4StringChipsInterface::Propagate is called"<<G4endl;
65#endif
66  // Protection for non physical conditions
67 
68  if(theSecondaries->size() == 1) 
69    throw G4HadronicException(__FILE__, __LINE__, "G4StringChipsInterface: Only one particle from String models!");
70 
71  // Calculate the mean energy lost
72  std::pair<G4double, G4double> theImpact = theNucleus->RefetchImpactXandY();
73  G4double impactX = theImpact.first;
74  G4double impactY = theImpact.second;
75  G4double inpactPar2 = impactX*impactX + impactY*impactY;
76 
77  G4double radius2 = theNucleus->GetNuclearRadius(5*perCent);
78  radius2 *= radius2;
79  G4double pathlength = 0;
80  if(radius2 - inpactPar2>0) pathlength = 2.*std::sqrt(radius2 - inpactPar2);
81  G4double theEnergyLostInFragmentation = theEnergyLossPerFermi*pathlength/fermi;
82 
83  // now select all particles in range
84  std::list<std::pair<G4double, G4KineticTrack *> > theSorted;
85  std::list<std::pair<G4double, G4KineticTrack *> >::iterator current;
86  for(unsigned int secondary = 0; secondary<theSecondaries->size(); secondary++)
87  {
88    G4LorentzVector a4Mom = theSecondaries->operator[](secondary)->Get4Momentum();
89
90#ifdef CHIPSdebug
91    G4cout <<"ALL STRING particles "<<theSecondaries->operator[](secondary)->GetDefinition()->GetPDGCharge()<<" "
92           << theSecondaries->operator[](secondary)->GetDefinition()->GetPDGEncoding()<<" "
93           << a4Mom <<G4endl; 
94#endif
95
96    G4double toSort = a4Mom.rapidity();
97    std::pair<G4double, G4KineticTrack *> it;
98    it.first = toSort;
99    it.second = theSecondaries->operator[](secondary);
100    G4bool inserted = false;
101    for(current = theSorted.begin(); current!=theSorted.end(); current++)
102    {
103      if((*current).first > toSort)
104      {
105               theSorted.insert(current, it);
106               inserted = true;
107               break;
108      }
109    }
110    if(!inserted)
111    {
112      theSorted.push_front(it);
113    }
114  }
115 
116  G4LorentzVector proj4Mom(0.,0.,0.,0.);
117  // @@ Use the G4QContent class, which is exactly (nD,nU,nS,nAD,nAU,nAS) !
118  // The G4QContent class is a basic clas in CHIPS (not PDG Code as in GEANT4),
119  // so in CHIPS on can has a hadronic obgect (Quasmon) with any Quark Content.
120  // As a simple extantion for the hadron (which is a special case for Hadron)
121  // there is a clas G4QChipolino, which is a Quasmon, which can decay in two
122  // hadrons. In future the three-hadron class can be added...
123  G4int nD  = 0;
124  G4int nU  = 0;
125  G4int nS  = 0;
126  G4int nAD = 0;
127  G4int nAU = 0;
128  G4int nAS = 0;
129  std::list<std::pair<G4double, G4KineticTrack *> >::iterator firstEscaping = theSorted.begin();
130  G4double runningEnergy = 0;
131  G4int particleCount = 0;
132  G4LorentzVector theLow = (*(theSorted.begin())).second->Get4Momentum();
133  G4LorentzVector theHigh;
134
135#ifdef CHIPSdebug
136  G4cout << "CHIPS ENERGY LOST "<<theEnergyLostInFragmentation<<G4endl;
137  G4cout << "sorted rapidities event start"<<G4endl;
138#endif
139
140  G4ReactionProductVector * theResult = new G4ReactionProductVector;
141  G4ReactionProduct * theSec;
142  G4KineticTrackVector * secondaries;
143#ifdef pdebug
144  G4cout<<"G4StringChipsInterface::Propagate: Absorption"<<G4endl; 
145#endif
146 
147  // first decay and add all escaping particles.
148  for(current = theSorted.begin(); current!=theSorted.end(); current++)
149  {
150    firstEscaping = current;
151    if((*current).second->GetDefinition()->GetQuarkContent(3)!=0 ||
152       (*current).second->GetDefinition()->GetAntiQuarkContent(3) !=0)
153    {
154      G4KineticTrack * aResult = (*current).second;
155      G4ParticleDefinition * pdef=aResult->GetDefinition();
156      secondaries = NULL;
157      if ( pdef->GetPDGWidth() > 0 && pdef->GetPDGLifeTime() < 5E-17*s )
158      {
159        secondaries = aResult->Decay();
160      }
161      if ( secondaries == NULL )
162      {
163        theSec = new G4ReactionProduct(aResult->GetDefinition());
164        G4LorentzVector current4Mom = aResult->Get4Momentum();
165        theSec->SetTotalEnergy(current4Mom.t());
166        theSec->SetMomentum(current4Mom.vect());
167        theResult->push_back(theSec);
168      } 
169      else
170      {
171        for (unsigned int aSecondary=0; aSecondary<secondaries->size(); aSecondary++)
172        {
173          theSec = new G4ReactionProduct(secondaries->operator[](aSecondary)->GetDefinition());
174          G4LorentzVector current4Mom = secondaries->operator[](aSecondary)->Get4Momentum();
175          theSec->SetTotalEnergy(current4Mom.t());
176          theSec->SetMomentum(current4Mom.vect());
177          theResult->push_back(theSec);
178        }
179        std::for_each(secondaries->begin(), secondaries->end(), DeleteKineticTrack());
180        delete secondaries;
181      }
182      continue;
183    }
184    runningEnergy += (*current).second->Get4Momentum().t();
185    if(runningEnergy > theEnergyLostInFragmentation) break;
186   
187#ifdef CHIPSdebug
188    G4cout <<"ABSORBED STRING particles "<<current->second->GetDefinition()->GetPDGCharge()<<" "
189           << current->second->GetDefinition()->GetPDGEncoding()<<" "
190                  << current->second->Get4Momentum() <<G4endl; 
191#endif
192#ifdef pdebug
193    G4cout<<"G4StringChipsInterface::Propagate:C="
194          <<current->second->GetDefinition()->GetPDGCharge()<<", PDG="
195          << current->second->GetDefinition()->GetPDGEncoding()<<", 4M="
196                 << current->second->Get4Momentum() <<G4endl; 
197#endif
198
199    // projectile 4-momentum needed in constructor of quasmon
200    particleCount++;
201    theHigh = (*current).second->Get4Momentum(); 
202    proj4Mom += (*current).second->Get4Momentum(); 
203
204#ifdef CHIPSdebug
205    G4cout << "sorted rapidities "<<current->second->Get4Momentum().rapidity()<<G4endl;
206#endif
207   
208     // projectile quark contents needed for G4QContent construction (@@ ? -> G4QContent class)
209    nD  += (*current).second->GetDefinition()->GetQuarkContent(1);
210    nU  += (*current).second->GetDefinition()->GetQuarkContent(2);
211    nS  += (*current).second->GetDefinition()->GetQuarkContent(3);
212    nAD += (*current).second->GetDefinition()->GetAntiQuarkContent(1);
213    nAU += (*current).second->GetDefinition()->GetAntiQuarkContent(2);
214    nAS += (*current).second->GetDefinition()->GetAntiQuarkContent(3);
215  }
216  // construct G4QContent
217
218#ifdef CHIPSdebug
219  G4cout << "Quark content: d="<<nD<<", u="<<nU<< ", s="<< nS << G4endl;
220  G4cout << "Anti-quark content: anit-d="<<nAD<<", anti-u="<<nAU<< ", anti-s="<< nAS << G4endl;
221#endif
222
223  G4QContent theProjectiles(nD, nU, nS, nAD, nAU, nAS);
224
225#ifdef CHIPSdebug
226  G4cout << "G4QContent is constructed"<<endl;
227#endif
228 
229  // target properties needed in constructor of quasmon
230  // remove all hit nucleons to get Target code
231  theNucleus->StartLoop();
232  G4Nucleon * aNucleon;
233  G4int resA = 0;
234  G4int resZ = 0;
235  G4ThreeVector hitMomentum(0,0,0);
236  G4double hitMass = 0;
237  G4int hitCount = 0;
238  while((aNucleon = theNucleus->GetNextNucleon()))
239  {
240    if(!aNucleon->AreYouHit())
241    {
242      resA++;
243      resZ+=G4int (aNucleon->GetDefinition()->GetPDGCharge());
244    }
245    else
246    {
247      hitMomentum += aNucleon->GetMomentum().vect();
248      hitMass += aNucleon->GetMomentum().m();
249      hitCount ++;
250    }
251  }
252  G4int targetPDGCode = 90000000 + 1000*resZ + (resA-resZ);
253  G4double targetMass = theNucleus->GetMass();
254  targetMass -= hitMass;
255  G4double targetEnergy = std::sqrt(hitMomentum.mag2()+targetMass*targetMass);
256  // !! @@ Target should be at rest: hitMomentum=(0,0,0) @@ !! M.K.
257  G4LorentzVector targ4Mom(-1.*hitMomentum, targetEnergy);
258 
259  // construct the quasmon
260  //  G4int nop = 122; // clusters up to Alpha cluster
261  G4int nop = 152; // V.Ivanchenko set the same parameter as for all CHIPS models
262  G4double fractionOfSingleQuasiFreeNucleons = 0.5; // It is A-dependent (C=.85, U=.40)
263  G4double fractionOfPairedQuasiFreeNucleons = 0.05;
264  G4double clusteringCoefficient = 5.;
265  G4double temperature = 180.;
266  G4double halfTheStrangenessOfSee = 0.3; // = s/d = s/u
267  G4double etaToEtaPrime = 0.3;
268
269  G4QNucleus::SetParameters(fractionOfSingleQuasiFreeNucleons,
270                            fractionOfPairedQuasiFreeNucleons,
271                            clusteringCoefficient);
272  G4Quasmon::SetParameters(temperature,
273                           halfTheStrangenessOfSee,
274                           etaToEtaPrime);
275
276#ifdef CHIPSdebug
277  G4cout << "G4QNucleus parameters "<< fractionOfSingleQuasiFreeNucleons << " "
278         << fractionOfPairedQuasiFreeNucleons << " "<< clusteringCoefficient << G4endl;
279  G4cout << "G4Quasmon parameters "<< temperature << " "<< halfTheStrangenessOfSee << " "
280         <<etaToEtaPrime << G4endl;
281  G4cout << "The Target PDG code = "<<targetPDGCode<<G4endl;
282  G4cout << "The projectile momentum = "<<1./MeV*proj4Mom<<G4endl;
283  G4cout << "The target momentum = "<<1./MeV*targ4Mom<<G4endl;
284#endif
285
286 
287  // Chips expects all in target rest frame, along z.
288  // G4QCHIPSWorld aWorld(nop);              // Create CHIPS World of nop particles
289  G4QCHIPSWorld::Get()->GetParticles(nop);
290  G4QHadronVector projHV;
291  // target rest frame
292  proj4Mom.boost(-1.*targ4Mom.boostVector());
293  // now go along z
294  G4LorentzRotation toZ;
295  toZ.rotateZ(-1*proj4Mom.phi());
296  toZ.rotateY(-1*proj4Mom.theta());
297  proj4Mom = toZ*proj4Mom;
298  G4LorentzRotation toLab(toZ.inverse());
299
300#ifdef CHIPSdebug
301  G4cout << "a Boosted projectile vector along z"<<proj4Mom<<" "<<proj4Mom.mag()<<G4endl;
302#endif
303 
304  G4QHadron* iH = new G4QHadron(theProjectiles, 1./MeV*proj4Mom);
305  projHV.push_back(iH);
306
307  // now call chips with this info in place
308  G4QHadronVector * output = 0;
309  if (particleCount!=0)
310  {
311    G4QEnvironment* pan= new G4QEnvironment(projHV, targetPDGCode);
312    try
313    {
314      output = pan->Fragment();
315    }
316    catch(G4HadronicException & aR)
317    {
318      G4cerr << "Exception thrown passing through G4ChiralInvariantPhaseSpace "<<G4endl;
319      G4cerr << " targetPDGCode = "<< targetPDGCode <<G4endl;
320      G4cerr << " Dumping the information in the pojectile list"<<G4endl;
321      for(size_t i=0; i< projHV.size(); i++)
322      {
323               G4cerr <<"  Incoming 4-momentum and PDG code of "<<i<<"'th hadron: "
324               <<" "<< projHV[i]->Get4Momentum()<<" "<<projHV[i]->GetPDGCode()<<G4endl;
325      }
326      throw;
327    }
328    std::for_each(projHV.begin(), projHV.end(), DeleteQHadron());
329    projHV.clear();
330    delete pan;
331  }
332  else output = new G4QHadronVector;
333   
334  // Fill the result.
335#ifdef CHIPSdebug
336  G4cout << "NEXT EVENT"<<endl;
337#endif
338 
339  // first decay and add all escaping particles.
340  for(current = firstEscaping; current!=theSorted.end(); current++)
341  {
342    G4KineticTrack * aResult = (*current).second;
343    G4ParticleDefinition * pdef=aResult->GetDefinition();
344    secondaries = NULL;
345    if ( pdef->GetPDGWidth() > 0 && pdef->GetPDGLifeTime() < 5E-17*s )
346    {
347      secondaries = aResult->Decay();
348    }
349    if ( secondaries == NULL )
350    {
351      theSec = new G4ReactionProduct(aResult->GetDefinition());
352      G4LorentzVector current4Mom = aResult->Get4Momentum();
353      current4Mom = toLab*current4Mom;
354      current4Mom.boost(targ4Mom.boostVector());
355      theSec->SetTotalEnergy(current4Mom.t());
356      theSec->SetMomentum(current4Mom.vect());
357      theResult->push_back(theSec);
358    } 
359    else
360    {
361      for (unsigned int aSecondary=0; aSecondary<secondaries->size(); aSecondary++)
362      {
363        theSec = new G4ReactionProduct(secondaries->operator[](aSecondary)->GetDefinition());
364        G4LorentzVector current4Mom = secondaries->operator[](aSecondary)->Get4Momentum();
365        current4Mom = toLab*current4Mom;
366        current4Mom.boost(targ4Mom.boostVector());
367        theSec->SetTotalEnergy(current4Mom.t());
368        theSec->SetMomentum(current4Mom.vect());
369        theResult->push_back(theSec);
370      }
371      std::for_each(secondaries->begin(), secondaries->end(), DeleteKineticTrack());
372      delete secondaries;
373    }
374  }
375  std::for_each(theSecondaries->begin(), theSecondaries->end(), DeleteKineticTrack());
376  delete theSecondaries;
377   
378  // now add the quasmon output
379#ifdef CHIPSdebug
380  G4cout << "Number of particles from string"<<theResult->size()<<G4endl;
381  G4cout << "Number of particles from chips"<<output->size()<<G4endl;
382#endif
383
384  for(unsigned int particle = 0; particle < output->size(); particle++)
385  {
386    if(output->operator[](particle)->GetNFragments() != 0) 
387    {
388      delete output->operator[](particle);
389      continue;
390    }
391    theSec = new G4ReactionProduct; 
392    G4int pdgCode = output->operator[](particle)->GetPDGCode();
393    G4ParticleDefinition * theDefinition;
394    // Note that I still have to take care of strange nuclei
395    // For this I need the mass calculation, and a changed interface
396    // for ion-table ==> work for Hisaya @@@@@@@
397    // Then I can sort out the pdgCode. I also need a decau process
398    // for strange nuclei; may be another chips interface
399    if(pdgCode>90000000) 
400    {
401      G4int aZ = (pdgCode-90000000)/1000;
402      if (aZ>1000) aZ=aZ%1000;  // patch for strange nuclei, to be repaired @@@@
403      G4int anN = pdgCode-90000000-1000*aZ;
404      if(anN>1000) anN=anN%1000; // patch for strange nuclei, to be repaired @@@@
405      if(pdgCode==91000000) theDefinition = G4Lambda::LambdaDefinition();
406      else if(pdgCode==92000000) theDefinition = G4Lambda::LambdaDefinition();
407      else if(pdgCode==93000000) theDefinition = G4Lambda::LambdaDefinition();
408      else if(pdgCode==94000000) theDefinition = G4Lambda::LambdaDefinition();
409      else if(pdgCode==95000000) theDefinition = G4Lambda::LambdaDefinition();
410      else if(pdgCode==96000000) theDefinition = G4Lambda::LambdaDefinition();
411      else if(pdgCode==97000000) theDefinition = G4Lambda::LambdaDefinition();
412      else if(pdgCode==98000000) theDefinition = G4Lambda::LambdaDefinition();
413      else if(aZ == 0 && anN == 1) theDefinition = G4Neutron::Neutron();
414      else theDefinition = G4ParticleTable::GetParticleTable()->FindIon(aZ,anN+aZ,0,aZ);
415    }   
416    else
417    {
418      theDefinition = G4ParticleTable::GetParticleTable()->FindParticle(output->operator[](particle)->GetPDGCode());
419    }
420#ifdef CHIPSdebug
421    G4cout << "Particle code produced = "<< pdgCode <<G4endl;
422#endif
423
424    theSec = new G4ReactionProduct(theDefinition);
425    G4LorentzVector current4Mom = output->operator[](particle)->Get4Momentum();
426    current4Mom = toLab*current4Mom;
427    current4Mom.boost(targ4Mom.boostVector());
428    theSec->SetTotalEnergy(current4Mom.t());
429    theSec->SetMomentum(current4Mom.vect());
430    theResult->push_back(theSec);
431   
432#ifdef CHIPSdebug
433    G4cout <<"CHIPS particles "<<theDefinition->GetPDGCharge()<<" "
434           << theDefinition->GetPDGEncoding()<<" "
435           << current4Mom <<G4endl; 
436#endif
437
438    delete output->operator[](particle);
439  }
440  delete output;
441#ifdef CHIPSdebug
442  G4cout << "Number of particles"<<theResult->size()<<G4endl;
443  G4cout << G4endl;
444  // @@ G4QContent has even the out option!
445  G4cout << "QUASMON preparation info "
446         << 1./MeV*proj4Mom<<" "
447         << 1./MeV*targ4Mom<<" "
448         << nD<<" "<<nU<<" "<<nS<<" "<<nAD<<" "<<nAU<<" "<<nAS<<" "
449         << hitCount<<" "
450         << particleCount<<" "
451         << theLow<<" "
452         << theHigh<<" "
453         << G4endl;
454#endif
455
456  return theResult;
457} 
Note: See TracBrowser for help on using the repository browser.