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

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

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