G4QStringChipsParticleLevelInterface Class Reference

#include <G4QStringChipsParticleLevelInterface.hh>

Inheritance diagram for G4QStringChipsParticleLevelInterface:

G4VIntraNuclearTransportModel G4HadronicInteraction

Public Member Functions

 G4QStringChipsParticleLevelInterface ()
virtual G4HadFinalState * ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &theNucleus)
virtual G4ReactionProductVector * Propagate (G4KineticTrackVector *theSecondaries, G4V3DNucleus *theNucleus)

Detailed Description

Definition at line 35 of file G4QStringChipsParticleLevelInterface.hh.

Constructor & Destructor Documentation

G4QStringChipsParticleLevelInterface::G4QStringChipsParticleLevelInterface (  ) 

Definition at line 48 of file G4QStringChipsParticleLevelInterface.cc.

References G4cin, G4cout, and G4endl.

00049 {
00050 #ifdef debug
00051   G4cout<<"G4QStringChipsParticleLevelInterface::Constructor is called"<<G4endl;
00052 #endif
00053   theEnergyLossPerFermi = 1.*GeV;
00054   nop = 152;                               // clusters (A<6)
00055   fractionOfSingleQuasiFreeNucleons = 0.5; // It is A-dependent (C=.85, U=.40)
00056   fractionOfPairedQuasiFreeNucleons = 0.05;
00057   clusteringCoefficient = 5.;
00058   temperature = 180.;
00059   halfTheStrangenessOfSee = 0.3; // = s/d = s/u
00060   etaToEtaPrime = 0.3;
00061   fusionToExchange = 1.;
00062   theInnerCoreDensityCut = 50.;
00063   
00064   if(getenv("ChipsParameterTuning"))
00065   {
00066     G4cout << "Please enter the energy loss per fermi in GeV"<<G4endl;
00067     G4cin >> theEnergyLossPerFermi;
00068     theEnergyLossPerFermi *= GeV;
00069     G4cout << "Please enter nop"<<G4endl;
00070     G4cin >> nop;
00071     G4cout << "Please enter the fractionOfSingleQuasiFreeNucleons"<<G4endl;
00072     G4cin >> fractionOfSingleQuasiFreeNucleons;
00073     G4cout << "Please enter the fractionOfPairedQuasiFreeNucleons"<<G4endl;
00074     G4cin >> fractionOfPairedQuasiFreeNucleons;
00075     G4cout << "Please enter the clusteringCoefficient"<<G4endl;
00076     G4cin >> clusteringCoefficient;
00077     G4cout << "Please enter the temperature"<<G4endl;
00078     G4cin >> temperature;
00079     G4cout << "Please enter halfTheStrangenessOfSee"<<G4endl;
00080     G4cin >> halfTheStrangenessOfSee;
00081     G4cout << "Please enter the etaToEtaPrime"<<G4endl;
00082     G4cin >> etaToEtaPrime;
00083     G4cout << "Please enter the fusionToExchange"<<G4endl;
00084     G4cin >> fusionToExchange;
00085     G4cout << "Please enter cut-off for calculating the nuclear radius in percent"<<G4endl;
00086     G4cin >> theInnerCoreDensityCut;
00087   }
00088 }

Member Function Documentation

G4HadFinalState * G4QStringChipsParticleLevelInterface::ApplyYourself ( const G4HadProjectile &  aTrack,
G4Nucleus &  theNucleus 
) [virtual]

Implements G4HadronicInteraction.

Definition at line 91 of file G4QStringChipsParticleLevelInterface.cc.

References G4ChiralInvariantPhaseSpace::ApplyYourself(), G4cout, and G4endl.

00092 {
00093 #ifdef debug
00094   G4cout<<"G4QStringChipsParticleLevelInterface::ApplyYourself is called"<<G4endl;
00095 #endif
00096   return theModel.ApplyYourself(aTrack, theNucleus);
00097 }

G4ReactionProductVector * G4QStringChipsParticleLevelInterface::Propagate ( G4KineticTrackVector *  theSecondaries,
G4V3DNucleus *  theNucleus 
) [virtual]

Implements G4VIntraNuclearTransportModel.

Definition at line 100 of file G4QStringChipsParticleLevelInterface.cc.

References G4Nucleon::AreYouHit(), G4ParticleTable::FindIon(), G4ParticleTable::FindParticle(), G4QEnvironment::Fragment(), G4cerr, G4cout, G4endl, G4QCHIPSWorld::Get(), G4KineticTrack::Get4Momentum(), G4ParticleDefinition::GetBaryonNumber(), G4KineticTrack::GetDefinition(), G4Nucleon::GetDefinition(), G4Nucleon::GetMomentum(), G4V3DNucleus::GetNextNucleon(), G4V3DNucleus::GetNuclearRadius(), G4V3DNucleus::GetOuterRadius(), G4QCHIPSWorld::GetParticles(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGEncoding(), G4ParticleDefinition::GetPDGMass(), G4ParticleDefinition::GetQuarkContent(), G4Lambda::Lambda(), G4Lambda::LambdaDefinition(), G4Neutron::Neutron(), G4Proton::Proton(), G4V3DNucleus::RefetchImpactXandY(), G4ReactionProduct::SetMomentum(), G4Quasmon::SetParameters(), G4QNucleus::SetParameters(), G4ReactionProduct::SetTotalEnergy(), and G4V3DNucleus::StartLoop().

00101 {
00102   static const G4double mProt=G4Proton::Proton()->GetPDGMass();
00103   static const G4double mNeut=G4Neutron::Neutron()->GetPDGMass();
00104   static const G4double mLamb=G4Lambda::Lambda()->GetPDGMass();
00105 #ifdef debug
00106   G4cout<<"G4QStringChipsParticleLevelInterface::Propagate is called"<<G4endl;
00107 #endif
00108   // Protection for non physical conditions
00109   
00110   if(theSecondaries->size() == 1) 
00111   {
00112     G4ReactionProductVector* theFastResult = new G4ReactionProductVector;
00113     G4ReactionProduct* theFastSec;
00114     theFastSec = new G4ReactionProduct((*theSecondaries)[0]->GetDefinition());
00115     G4LorentzVector current4Mom = (*theSecondaries)[0]->Get4Momentum();
00116     theFastSec->SetTotalEnergy(current4Mom.t());
00117     theFastSec->SetMomentum(current4Mom.vect());
00118     theFastResult->push_back(theFastSec);
00119     return theFastResult;
00120     //throw G4HadronicException(__FILE__,__LINE__,
00121     //      "G4QStringChipsParticleLevelInterface: Only one particle from String models!");
00122   }
00123   
00124   // target properties needed in constructor of quasmon, and for boosting to
00125   // target rest frame
00126   // remove all nucleons already involved in STRING interaction, to make the ResidualTarget
00127   theNucleus->StartLoop();
00128   G4Nucleon * aNucleon;
00129   G4int resA = 0;
00130   G4int resZ = 0;
00131   G4ThreeVector hitMomentum(0,0,0);
00132   G4double hitMass = 0;
00133   unsigned int hitCount = 0;
00134   while((aNucleon = theNucleus->GetNextNucleon()))
00135   {
00136     if(!aNucleon->AreYouHit())
00137     {
00138       resA++;                                                  // Collect A of the ResidNuc
00139       resZ+=G4int (aNucleon->GetDefinition()->GetPDGCharge()); // Collect Z of the ResidNuc
00140     }
00141     else
00142     {
00143       hitMomentum += aNucleon->GetMomentum().vect();           // Sum 3-mom of StringHadr's
00144       hitMass += aNucleon->GetMomentum().m();                  // Sum masses of StringHadrs
00145       hitCount ++;                                             // Calculate STRING hadrons
00146     }
00147   }
00148   G4int targetPDGCode = 90000000 + 1000*resZ + (resA-resZ);    // PDG of theResidualNucleus
00149   G4double targetMass=mNeut;
00150   if (!resZ)                                                   // Nucleus of only neutrons
00151   {
00152     if (resA>1) targetMass*=resA;
00153   }
00154   else targetMass=G4ParticleTable::GetParticleTable()->FindIon(resZ,resA,0,resZ)
00155                                                                             ->GetPDGMass();
00156   G4double targetEnergy = std::sqrt(hitMomentum.mag2()+targetMass*targetMass);
00157   // !! @@ Target should be at rest: hitMomentum=(0,0,0) @@ !! M.K. (go to this system)
00158   G4LorentzVector targ4Mom(-1.*hitMomentum, targetEnergy);
00159   
00160   // Calculate the mean energy lost
00161   std::pair<G4double, G4double> theImpact = theNucleus->RefetchImpactXandY();
00162   G4double impactX = theImpact.first;
00163   G4double impactY = theImpact.second;
00164   G4double impactPar2 = impactX*impactX + impactY*impactY;
00165   
00166   G4double radius2 = theNucleus->GetNuclearRadius(theInnerCoreDensityCut*perCent);
00167   radius2 *= radius2;
00168   G4double pathlength = 0.;
00169 #ifdef pdebug
00170   G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: r="<<std::sqrt(radius2)/fermi
00171         <<", b="<<std::sqrt(impactPar2)/fermi<<", R="<<theNucleus->GetOuterRadius()/fermi
00172         <<", b/r="<<std::sqrt(impactPar2/radius2)<<G4endl; 
00173 #endif
00174   if(radius2 - impactPar2>0) pathlength = 2.*std::sqrt(radius2 - impactPar2);
00175   G4double theEnergyLostInFragmentation = theEnergyLossPerFermi*pathlength/fermi;
00176   
00177   // now select all particles in range
00178   std::list<std::pair<G4double, G4KineticTrack *> > theSorted;         // Output
00179   std::list<std::pair<G4double, G4KineticTrack *> >::iterator current; // Input
00180   for(unsigned int secondary = 0; secondary<theSecondaries->size(); secondary++)
00181   {
00182     G4LorentzVector a4Mom = theSecondaries->operator[](secondary)->Get4Momentum();
00183 #ifdef CHIPSdebug
00184     G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: ALL STRING particles "
00185           << theSecondaries->operator[](secondary)->GetDefinition()->GetPDGCharge()<<" "
00186           << theSecondaries->operator[](secondary)->GetDefinition()->GetPDGEncoding()<<" "
00187           << a4Mom <<G4endl; 
00188 #endif
00189 #ifdef pdebug
00190     G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: in C="
00191           <<theSecondaries->operator[](secondary)->GetDefinition()->GetPDGCharge()<<",PDG="
00192           <<theSecondaries->operator[](secondary)->GetDefinition()->GetPDGEncoding()
00193           <<",4M="<<a4Mom<<", current nS="<<theSorted.size()<<G4endl; 
00194 #endif
00195     G4double toSort = a4Mom.rapidity();          // Rapidity is used for the ordering (?!)
00196     std::pair<G4double, G4KineticTrack *> it;
00197     it.first = toSort;
00198     it.second = theSecondaries->operator[](secondary);
00199     G4bool inserted = false;
00200     for(current = theSorted.begin(); current!=theSorted.end(); current++)
00201     {
00202       if((*current).first > toSort)        // The current is smaller then existing
00203       {
00204         theSorted.insert(current, it);     // It shifts the others up
00205         inserted = true;
00206         break;
00207       }
00208     }
00209     if(!inserted) theSorted.push_back(it); // It is bigger than any previous
00210   }
00211   
00212   G4LorentzVector proj4Mom(0.,0.,0.,0.);
00213   G4int nD  = 0;
00214   G4int nU  = 0;
00215   G4int nS  = 0;
00216   G4int nAD = 0;
00217   G4int nAU = 0;
00218   G4int nAS = 0;
00219   std::list<std::pair<G4double,G4KineticTrack*> >::iterator firstEscape=theSorted.begin();
00220   G4double runningEnergy = 0;
00221   G4int particleCount = 0;
00222   G4LorentzVector theLow = (*(theSorted.begin())).second->Get4Momentum();
00223   G4LorentzVector theHigh;
00224 
00225 #ifdef CHIPSdebug
00226   G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: CHIPS ENERGY LOST "
00227         <<theEnergyLostInFragmentation<<". Sorted rapidities event start"<<G4endl;
00228 #endif
00229 #ifdef pdebug
00230   G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: total CHIPS energy = "
00231         <<theEnergyLostInFragmentation<<". Start rapidity sorting nS="<<theSorted.size()
00232         <<G4endl;
00233 #endif
00234 
00235   G4QHadronVector projHV;
00236   std::vector<G4QContent> theContents;
00237   std::vector<G4LorentzVector*> theMomenta;
00238   G4ReactionProductVector* theResult = new G4ReactionProductVector;
00239   G4ReactionProduct* theSec;
00240   G4KineticTrackVector* secondaries;
00241 #ifdef pdebug
00242   G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: Absorption nS="
00243         <<theSorted.size()<<G4endl; 
00244 #endif
00245   G4bool EscapeExists = false;
00246   for(current = theSorted.begin(); current!=theSorted.end(); current++)
00247   {
00248 #ifdef pdebug
00249     G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: nq="
00250           <<(*current).second->GetDefinition()->GetQuarkContent(3)<<", naq="
00251           <<(*current).second->GetDefinition()->GetAntiQuarkContent(3)<<", PDG="
00252           <<(*current).second->GetDefinition()->GetPDGEncoding()<<",4M="
00253           <<(*current).second->Get4Momentum()<<G4endl; 
00254 #endif
00255     firstEscape = current;              // Remember to make decays for the rest
00256     G4KineticTrack* aResult = (*current).second;
00257     // This is an old (H.P.) solution, which makes an error in En/Mom conservation
00258     //
00259     // @@ Now it does not include strange particle for the absorption in nuclei (?!) M.K.
00260     //if((*current).second->GetDefinition()->GetQuarkContent(3)!=0 || 
00261     //   (*current).second->GetDefinition()->GetAntiQuarkContent(3) !=0) // Strange quarks
00262     //{
00263     //  G4ParticleDefinition* pdef = aResult->GetDefinition();
00264     //  secondaries = NULL;
00265     //  if ( pdef->GetPDGWidth() > 0 && pdef->GetPDGLifeTime() < 5E-17*s )
00266     //     secondaries = aResult->Decay(); // @@ Decay of only strange resonances (?!) M.K.
00267     //  if ( secondaries == NULL )           // No decay
00268     //  {
00269     //    theSec = new G4ReactionProduct(aResult->GetDefinition());
00270     //    G4LorentzVector current4Mom = aResult->Get4Momentum();
00271     //    current4Mom.boost(targ4Mom.boostVector()); // boost from the targetAtRes system
00272     //    theSec->SetTotalEnergy(current4Mom.t());
00273     //    theSec->SetMomentum(current4Mom.vect());
00274     //    theResult->push_back(theSec);
00275     //  } 
00276     //  else                                 // The decay happened
00277     //  {
00278     //   for (unsigned int aSecondary=0; aSecondary<secondaries->size(); aSecondary++)
00279     //   {
00280     //     theSec = 
00281     //        new G4ReactionProduct(secondaries->operator[](aSecondary)->GetDefinition());
00282     //     G4LorentzVector current4Mom=secondaries->operator[](aSecondary)->Get4Momentum();
00283     //     current4Mom.boost(targ4Mom.boostVector());
00284     //     theSec->SetTotalEnergy(current4Mom.t());
00285     //     theSec->SetMomentum(current4Mom.vect());
00286     //     theResult->push_back(theSec);
00287     //   }
00288     //   std::for_each(secondaries->begin(), secondaries->end(), DeleteKineticTrack());
00289     //   delete secondaries;
00290     //  }
00291     //}
00292     //
00293     //runningEnergy += (*current).second->Get4Momentum().t();
00294     //if((*current).second->GetDefinition() == G4Proton::Proton())
00295     //                                   runningEnergy-=G4Proton::Proton()->GetPDGMass();
00296     //if((*current).second->GetDefinition() == G4Neutron::Neutron())
00297     //                                   runningEnergy-=G4Neutron::Neutron()->GetPDGMass();
00298     //if((*current).second->GetDefinition() == G4Lambda::Lambda())
00299     //                                   runningEnergy-=G4Lambda::Lambda()->GetPDGMass();
00300     //
00301     // New solution starts from here (M.Kossov March 2006) [Strange particles included]
00302     runningEnergy += aResult->Get4Momentum().t();
00303     G4double charge=aResult->GetDefinition()->GetPDGCharge(); // Charge of the particle
00304     G4int strang=aResult->GetDefinition()->GetQuarkContent(3);// Its strangeness
00305     G4int baryn=aResult->GetDefinition()->GetBaryonNumber();  // Its baryon number
00306     if     (baryn>0 && charge>0 && strang<1) runningEnergy-=mProt;  // For positive baryons
00307     else if(baryn>0 && strang<1) runningEnergy-=mNeut;              // For neut/neg baryons
00308     else if(baryn>0) runningEnergy-=mLamb;                          // For strange baryons
00309     else if(baryn<0) runningEnergy+=mProt;                          // For anti-particles
00310     // ------------ End of the new solution
00311 #ifdef CHIPSdebug
00312     G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: sorted rapidities "
00313                                     <<(*current).second->Get4Momentum().rapidity()<<G4endl;
00314 #endif
00315 
00316 #ifdef pdebug
00317     G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: E="<<runningEnergy<<", EL="
00318                                                     <<theEnergyLostInFragmentation<<G4endl;
00319 #endif
00320     if(runningEnergy > theEnergyLostInFragmentation)
00321     {
00322       EscapeExists = true;
00323       break;
00324     }
00325 #ifdef CHIPSdebug
00326     G4cout <<"G4QStringChipsParticleLevelInterface::Propagate: ABSORBED STRING particles "
00327            <<(*current).second->GetDefinition()->GetPDGCharge()<<" "
00328            << (*current).second->GetDefinition()->GetPDGEncoding()<<" "
00329            << (*current).second->Get4Momentum() <<G4endl; 
00330 #endif
00331 #ifdef pdebug
00332     G4cout<<"G4QStringChipsParticleLevelInterface::Propagate:C="
00333           <<current->second->GetDefinition()->GetPDGCharge()<<", PDG="
00334           <<current->second->GetDefinition()->GetPDGEncoding()<<", 4M="
00335           <<current->second->Get4Momentum()<<G4endl; 
00336 #endif
00337 
00338     // projectile 4-momentum in target rest frame needed in constructor of QHadron
00339     particleCount++;
00340     theHigh = (*current).second->Get4Momentum(); 
00341     proj4Mom = (*current).second->Get4Momentum(); 
00342     proj4Mom.boost(-1.*targ4Mom.boostVector());   // Back to the system of nucleusAtRest
00343     nD = (*current).second->GetDefinition()->GetQuarkContent(1);
00344     nU = (*current).second->GetDefinition()->GetQuarkContent(2);
00345     nS = (*current).second->GetDefinition()->GetQuarkContent(3);
00346     nAD = (*current).second->GetDefinition()->GetAntiQuarkContent(1);
00347     nAU = (*current).second->GetDefinition()->GetAntiQuarkContent(2);
00348     nAS = (*current).second->GetDefinition()->GetAntiQuarkContent(3);
00349     G4QContent aProjectile(nD, nU, nS, nAD, nAU, nAS);
00350 
00351 #ifdef CHIPSdebug_1
00352     G4cout <<G4endl;
00353     G4cout <<"G4QStringChipsParticleLevelInterface::Propagate: Quark content: d="<<nD
00354            <<", u="<<nU<<", s="<<nS<< "Anti-quark content: anit-d="<<nAD<<", anti-u="<<nAU
00355            <<", anti-s="<<nAS<<". G4QContent is constructed"<<endl;
00356 #endif
00357 
00358     theContents.push_back(aProjectile);
00359     G4LorentzVector* aVec = new G4LorentzVector((1./MeV)*proj4Mom); // @@ MeV is basic
00360 
00361 #ifdef CHIPSdebug_1
00362     G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: projectile momentum = "
00363           <<*aVec<<G4endl;
00364     G4cout << G4endl;
00365 #endif
00366    
00367     theMomenta.push_back(aVec);
00368   }
00369   std::vector<G4QContent> theFinalContents;
00370   std::vector<G4LorentzVector*> theFinalMomenta;
00371 
00372   for(unsigned int hp = 0; hp<theContents.size(); hp++)
00373   {
00374     G4QHadron* aHadron = new G4QHadron(theContents[hp], *(theMomenta[hp]) );
00375     projHV.push_back(aHadron);
00376   }
00377  
00378   // construct the quasmon
00379   size_t i;
00380   for (i=0; i<theFinalMomenta.size(); i++) delete theFinalMomenta[i];
00381   for (i=0; i<theMomenta.size();      i++) delete theMomenta[i];
00382   theFinalMomenta.clear();
00383   theMomenta.clear();
00384 
00385   G4QNucleus::SetParameters(fractionOfSingleQuasiFreeNucleons,
00386                             fractionOfPairedQuasiFreeNucleons,
00387                             clusteringCoefficient, 
00388                             fusionToExchange);
00389   G4Quasmon::SetParameters(temperature, halfTheStrangenessOfSee, etaToEtaPrime);
00390 
00391 #ifdef CHIPSdebug
00392   G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: G4QNucleus parameters "
00393         <<fractionOfSingleQuasiFreeNucleons<<" "<<fractionOfPairedQuasiFreeNucleons
00394         <<" "<<clusteringCoefficient<<G4endl;
00395   G4cout<<"G4Quasmon parameters "<<temperature<<" "<<halfTheStrangenessOfSee<<" "
00396         <<etaToEtaPrime << G4endl;
00397   G4cout<<"The Target PDG code = "<<targetPDGCode<<G4endl;
00398   G4cout<<"The projectile momentum = "<<1./MeV*proj4Mom<<G4endl;
00399   G4cout<<"The target momentum = "<<1./MeV*targ4Mom<<G4endl;
00400 #endif
00401 
00402   // now call chips with this info in place
00403   G4QHadronVector* output = 0;
00404   if (particleCount!=0 && resA!=0)
00405   {
00406     //  G4QCHIPSWorld aWorld(nop);              // Create CHIPS World of nop particles
00407     G4QCHIPSWorld::Get()->GetParticles(nop);
00408     G4QEnvironment* pan= new G4QEnvironment(projHV, targetPDGCode);
00409 #ifdef pdebug
00410       G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: CHIPS fragmentation, rA="
00411             <<resA<<", #AbsPt="<<particleCount<<G4endl; 
00412 #endif
00413     try
00414     {
00415       output = pan->Fragment();                 // The main fragmentation member function
00416     }
00417     catch(G4HadronicException& aR)
00418     {
00419       G4cerr << "Exception thrown of G4QStringChipsParticleLevelInterface "<<G4endl;
00420       G4cerr << " targetPDGCode = "<< targetPDGCode <<G4endl;
00421       G4cerr << " The projectile momentum = "<<1./MeV*proj4Mom<<G4endl;
00422       G4cerr << " The target momentum = "<<1./MeV*targ4Mom<<G4endl<<G4endl;
00423       G4cerr << " Dumping the information in the pojectile list"<<G4endl;
00424       for(i=0; i< projHV.size(); i++)
00425       {
00426         G4cerr <<"  Incoming 4-momentum and PDG code of "<<i<<"'th hadron: "
00427         <<" "<< projHV[i]->Get4Momentum()<<" "<<projHV[i]->GetPDGCode()<<G4endl;
00428       }
00429       throw;
00430     }
00431     // clean up particles
00432     std::for_each(projHV.begin(), projHV.end(), DeleteQHadron());
00433     projHV.clear();
00434     delete pan;
00435   }
00436   else
00437   {
00438 #ifdef pdebug
00439     G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: NO CHIPS fragmentation, rA="
00440           <<resA<<", #AbsPt="<<particleCount<<G4endl; 
00441 #endif
00442     output = new G4QHadronVector;
00443   }   
00444   // Fill the result.
00445 #ifdef CHIPSdebug
00446   G4cout << "NEXT EVENT, EscapeExists="<<EscapeExists<<G4endl;
00447 #endif
00448 
00449   // first decay and add all escaping particles.
00450   if (EscapeExists) for (current = firstEscape; current!=theSorted.end(); current++)
00451   {
00452     G4KineticTrack* aResult = (*current).second;
00453     G4ParticleDefinition* pdef=aResult->GetDefinition();
00454     secondaries = NULL;
00455     if(pdef->GetPDGWidth() > 0 && pdef->GetPDGLifeTime() < 5E-17*s )
00456     {
00457       secondaries = aResult->Decay();  // @@ Uses standard Decay, which is now wrong!
00458     }
00459     if ( secondaries == NULL )
00460     {
00461       theSec = new G4ReactionProduct(aResult->GetDefinition());
00462       G4LorentzVector current4Mom = aResult->Get4Momentum();
00463       current4Mom.boost(targ4Mom.boostVector());
00464       theSec->SetTotalEnergy(current4Mom.t());
00465       theSec->SetMomentum(current4Mom.vect());
00466 #ifdef pdebug
00467       G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: *OUT* QGS stable PDG="
00468             <<aResult->GetDefinition()->GetPDGEncoding()<<",4M="<<current4Mom<<G4endl; 
00469 #endif
00470       theResult->push_back(theSec);
00471     } 
00472     else
00473     {
00474       for (unsigned int aSecondary=0; aSecondary<secondaries->size(); aSecondary++)
00475       {
00476         theSec=new G4ReactionProduct(secondaries->operator[](aSecondary)->GetDefinition());
00477         G4LorentzVector current4Mom = secondaries->operator[](aSecondary)->Get4Momentum();
00478         current4Mom.boost(targ4Mom.boostVector());
00479         theSec->SetTotalEnergy(current4Mom.t());
00480         theSec->SetMomentum(current4Mom.vect());
00481 #ifdef pdebug
00482         G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: *OUT* QGS decay PDG="
00483               <<secondaries->operator[](aSecondary)->GetDefinition()->GetPDGEncoding()
00484               <<",4M="<<current4Mom<<G4endl; 
00485 #endif
00486         theResult->push_back(theSec);
00487       }
00488       std::for_each(secondaries->begin(), secondaries->end(), DeleteKineticTrack());
00489       delete secondaries;
00490     }
00491   }
00492   std::for_each(theSecondaries->begin(), theSecondaries->end(), DeleteKineticTrack());
00493   delete theSecondaries;
00494     
00495   // now add the quasmon output
00496   G4int maxParticle=output->size();
00497 #ifdef CHIPSdebug
00498   G4cout << "Number of particles from string"<<theResult->size()<<G4endl;
00499   G4cout << "Number of particles from chips"<<maxParticle<<G4endl;
00500 #endif
00501 #ifdef pdebug
00502   G4cout << "Number of particles from QGS="<<theResult->size()<<G4endl;
00503   G4cout << "Number of particles from CHIPS="<<maxParticle<<G4endl;
00504 #endif
00505   if(maxParticle) for(G4int particle = 0; particle < maxParticle; particle++)
00506   {
00507     if(output->operator[](particle)->GetNFragments() != 0) 
00508     {
00509       delete output->operator[](particle);
00510       continue;
00511     }
00512     G4int pdgCode = output->operator[](particle)->GetPDGCode();
00513 
00514 
00515 #ifdef CHIPSdebug
00516     G4cerr << "PDG code of chips particle = "<<pdgCode<<G4endl;
00517 #endif
00518 
00519     G4ParticleDefinition * theDefinition;
00520     // Note that I still have to take care of strange nuclei
00521     // For this I need the mass calculation, and a changed interface
00522     // for ion-table ==> work for Hisaya @@@@@@@
00523     // Then I can sort out the pdgCode. I also need a decau process 
00524     // for strange nuclei; may be another chips interface
00525     if(pdgCode>90000000) 
00526     {
00527       G4int aZ = (pdgCode-90000000)/1000;
00528       if (aZ>1000) aZ=aZ%1000;  // patch for strange nuclei, to be repaired @@@@
00529       G4int anN = pdgCode-90000000-1000*aZ;
00530       if(anN>1000) anN=anN%1000; // patch for strange nuclei, to be repaired @@@@
00531       if(pdgCode==91000000) theDefinition = G4Lambda::LambdaDefinition();
00532       else if(pdgCode==92000000) theDefinition = G4Lambda::LambdaDefinition();
00533       else if(pdgCode==93000000) theDefinition = G4Lambda::LambdaDefinition();
00534       else if(pdgCode==94000000) theDefinition = G4Lambda::LambdaDefinition();
00535       else if(pdgCode==95000000) theDefinition = G4Lambda::LambdaDefinition();
00536       else if(pdgCode==96000000) theDefinition = G4Lambda::LambdaDefinition();
00537       else if(pdgCode==97000000) theDefinition = G4Lambda::LambdaDefinition();
00538       else if(pdgCode==98000000) theDefinition = G4Lambda::LambdaDefinition();
00539       else if(aZ == 0 && anN == 1) theDefinition = G4Neutron::Neutron();
00540       else theDefinition = G4ParticleTable::GetParticleTable()->FindIon(aZ,anN+aZ,0,aZ);
00541     }    
00542     else theDefinition = G4ParticleTable::GetParticleTable()->FindParticle(pdgCode);
00543 
00544 #ifdef CHIPSdebug
00545     G4cout << "Particle code produced = "<< pdgCode <<G4endl;
00546 #endif
00547 
00548     if(theDefinition)
00549     {
00550       theSec = new G4ReactionProduct(theDefinition);
00551       G4LorentzVector current4Mom = output->operator[](particle)->Get4Momentum();
00552       current4Mom.boost(targ4Mom.boostVector());
00553       theSec->SetTotalEnergy(current4Mom.t());
00554       theSec->SetMomentum(current4Mom.vect());
00555 #ifdef pdebug
00556       G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: *OUT* CHIPS PDG="
00557               <<theDefinition->GetPDGEncoding()<<",4M="<<current4Mom<<G4endl; 
00558 #endif
00559       theResult->push_back(theSec);
00560     }
00561     else
00562     {
00563       G4cerr << G4endl<<"WARNING: "<<G4endl;
00564       G4cerr << "Getting unknown pdgCode from chips in ParticleLevelInterface"<<G4endl;
00565       G4cerr << "skipping particle with pdgCode = "<<pdgCode<<G4endl<<G4endl;
00566     }
00567     
00568 #ifdef CHIPSdebug
00569     G4cout <<"CHIPS particles "<<theDefinition->GetPDGCharge()<<" "
00570            << theDefinition->GetPDGEncoding()<<" "
00571     << output->operator[](particle)->Get4Momentum() <<G4endl; 
00572 #endif
00573 
00574     delete output->operator[](particle);
00575   }
00576   else
00577   {
00578     if(resA>0)
00579     {
00580       G4ParticleDefinition* theDefinition = G4Neutron::Neutron();
00581       if(resA==1) // The residual nucleus at rest must be added to conserve BaryN & Charge
00582       {
00583         if(resZ == 1) theDefinition = G4Proton::Proton();
00584       }
00585       else theDefinition = G4ParticleTable::GetParticleTable()->FindIon(resZ,resA,0,resZ);
00586       theSec = new G4ReactionProduct(theDefinition);
00587       theSec->SetTotalEnergy(theDefinition->GetPDGMass());
00588       theSec->SetMomentum(G4ThreeVector(0.,0.,0.));
00589       theResult->push_back(theSec);
00590       if(!resZ && resA>0) for(G4int ni=1; ni<resA; ni++) 
00591       {
00592         theSec = new G4ReactionProduct(theDefinition);
00593         theSec->SetTotalEnergy(theDefinition->GetPDGMass());
00594         theSec->SetMomentum(G4ThreeVector(0.,0.,0.));
00595         theResult->push_back(theSec);
00596       }
00597     }
00598   }
00599   delete output;
00600 
00601 #ifdef CHIPSdebug
00602   G4cout << "Number of particles"<<theResult->size()<<G4endl;
00603   G4cout << G4endl;
00604   G4cout << "QUASMON preparation info "
00605          << 1./MeV*proj4Mom<<" "
00606   << 1./MeV*targ4Mom<<" "
00607   << nD<<" "<<nU<<" "<<nS<<" "<<nAD<<" "<<nAU<<" "<<nAS<<" "
00608   << hitCount<<" "
00609   << particleCount<<" "
00610   << theLow<<" "
00611   << theHigh<<" "
00612   << G4endl;
00613 #endif
00614 
00615   return theResult;
00616 }

The documentation for this class was generated from the following files:
Generated on Mon May 27 17:53:14 2013 for Geant4 by  doxygen 1.4.7