G4StringChipsParticleLevelInterface Class Reference

#include <G4StringChipsParticleLevelInterface.hh>

Inheritance diagram for G4StringChipsParticleLevelInterface:

G4VIntraNuclearTransportModel G4HadronicInteraction

Public Member Functions

 G4StringChipsParticleLevelInterface ()
virtual G4HadFinalStateApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &theNucleus)
virtual G4ReactionProductVectorPropagate (G4KineticTrackVector *theSecondaries, G4V3DNucleus *theNucleus)

Detailed Description

Definition at line 38 of file G4StringChipsParticleLevelInterface.hh.


Constructor & Destructor Documentation

G4StringChipsParticleLevelInterface::G4StringChipsParticleLevelInterface (  ) 

Definition at line 66 of file G4StringChipsParticleLevelInterface.cc.

References G4cin, G4cout, and G4endl.

00067 {
00068 #ifdef debug
00069   G4cout<<"G4StringChipsParticleLevelInterface::Constructor is called"<<G4endl;
00070 #endif
00071   //theEnergyLossPerFermi = 1.*GeV;
00072   theEnergyLossPerFermi = 1.5*GeV;
00073   nop = 152;                               // clusters (A<6)
00074   fractionOfSingleQuasiFreeNucleons = 0.5; // It is A-dependent (C=.85, U=.40) - M.K.
00075   fractionOfPairedQuasiFreeNucleons = 0.05;
00076   clusteringCoefficient = 5.;
00077   temperature = 180.;
00078   halfTheStrangenessOfSee = 0.3; // = s/d = s/u
00079   etaToEtaPrime = 0.3;
00080   fusionToExchange = 1.;
00081   //theInnerCoreDensityCut = 50.;
00082   theInnerCoreDensityCut = 70.;
00083   
00084   if(getenv("ChipsParameterTuning"))
00085   {
00086     G4cout << "Please enter the energy loss per fermi in GeV"<<G4endl;
00087     G4cin >> theEnergyLossPerFermi;
00088     theEnergyLossPerFermi *= GeV;
00089     G4cout << "Please enter nop"<<G4endl;
00090     G4cin >> nop;
00091     G4cout << "Please enter the fractionOfSingleQuasiFreeNucleons"<<G4endl;
00092     G4cin >> fractionOfSingleQuasiFreeNucleons;
00093     G4cout << "Please enter the fractionOfPairedQuasiFreeNucleons"<<G4endl;
00094     G4cin >> fractionOfPairedQuasiFreeNucleons;
00095     G4cout << "Please enter the clusteringCoefficient"<<G4endl;
00096     G4cin >> clusteringCoefficient;
00097     G4cout << "Please enter the temperature"<<G4endl;
00098     G4cin >> temperature;
00099     G4cout << "Please enter halfTheStrangenessOfSee"<<G4endl;
00100     G4cin >> halfTheStrangenessOfSee;
00101     G4cout << "Please enter the etaToEtaPrime"<<G4endl;
00102     G4cin >> etaToEtaPrime;
00103     G4cout << "Please enter the fusionToExchange"<<G4endl;
00104     G4cin >> fusionToExchange;
00105     G4cout << "Please enter the cut-off for calculating the nuclear radius in percent"<<G4endl;
00106     G4cin >> theInnerCoreDensityCut;
00107   }
00108 }


Member Function Documentation

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

Implements G4HadronicInteraction.

Definition at line 111 of file G4StringChipsParticleLevelInterface.cc.

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

00112 {
00113 #ifdef debug
00114   G4cout<<"G4StringChipsParticleLevelInterface::ApplyYourself is called"<<G4endl;
00115 #endif
00116   return theModel.ApplyYourself(aTrack, theNucleus);
00117 }

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

Implements G4VIntraNuclearTransportModel.

Definition at line 120 of file G4StringChipsParticleLevelInterface.cc.

References G4AntiKaonZero::AntiKaonZeroDefinition(), 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(), G4ParticleDefinition::IsShortLived(), G4KaonMinus::KaonMinusDefinition(), G4KaonPlus::KaonPlus(), G4KaonPlus::KaonPlusDefinition(), G4KaonZero::KaonZero(), G4KaonZero::KaonZeroDefinition(), G4Lambda::Lambda(), G4Lambda::LambdaDefinition(), G4Neutron::Neutron(), G4PionMinus::PionMinus(), G4PionMinus::PionMinusDefinition(), G4PionPlus::PionPlusDefinition(), G4Proton::Proton(), G4V3DNucleus::RefetchImpactXandY(), G4ReactionProduct::SetMomentum(), G4Quasmon::SetParameters(), G4QNucleus::SetParameters(), G4ReactionProduct::SetTotalEnergy(), and G4V3DNucleus::StartLoop().

00121 {
00122   static const G4double mProt=G4Proton::Proton()->GetPDGMass();
00123   static const G4double mNeut=G4Neutron::Neutron()->GetPDGMass();
00124   static const G4double mLamb=G4Lambda::Lambda()->GetPDGMass();
00125   static const G4double mKChg=G4KaonPlus::KaonPlus()->GetPDGMass();
00126   static const G4double mKZer=G4KaonZero::KaonZero()->GetPDGMass();
00127   static const G4double mPiCh=G4PionMinus::PionMinus()->GetPDGMass();
00128   static const G4int pcl=4; // clusterization parameter for Energy Flow
00129   static const G4QContent ProtQC(1,2,0,0,0,0);
00130   static const G4QContent NeutQC(2,1,0,0,0,0);
00131   static const G4QContent LambQC(1,1,1,0,0,0);
00132   static const G4QContent KPlsQC(0,1,0,0,0,1);
00133   static const G4QContent KMinQC(0,0,1,0,1,0);
00134   static const G4QContent AKZrQC(1,0,0,0,0,1);
00135   static const G4QContent KZerQC(1,0,0,0,0,1);
00136   static const G4QContent PiMiQC(1,0,0,0,1,0);
00137   static const G4QContent PiPlQC(0,1,0,1,0,0);
00138 #ifdef debug
00139   G4cout<<"G4StringChipsParticleLevelInterface::Propagate is called"<<G4endl;
00140 #endif
00141   // Protection for non physical conditions
00142   
00143   if(theSecondaries->size() == 1) 
00144   {
00145     G4ReactionProductVector* theFastResult = new G4ReactionProductVector;
00146     G4ReactionProduct* theFastSec;
00147     theFastSec = new G4ReactionProduct((*theSecondaries)[0]->GetDefinition());
00148     G4LorentzVector current4Mom = (*theSecondaries)[0]->Get4Momentum();
00149     theFastSec->SetTotalEnergy(current4Mom.t());
00150     theFastSec->SetMomentum(current4Mom.vect());
00151     theFastResult->push_back(theFastSec);
00152     return theFastResult;
00153     //throw G4HadronicException(__FILE__,__LINE__,
00154     //       "G4StringChipsParticleLevelInterface: Only one particle from String models!");
00155   }
00156   
00157   // target properties needed in constructor of quasmon, and for boosting to
00158   // target rest frame
00159   // remove all nucleons already involved in STRING interaction, to make the ResidualTarget
00160   theNucleus->StartLoop();
00161   G4Nucleon * aNucleon;
00162   G4int resA = 0;
00163   G4int resZ = 0;
00164   G4ThreeVector hitMomentum(0,0,0);
00165   G4double hitMass = 0;
00166   unsigned int hitCount = 0;
00167   while((aNucleon = theNucleus->GetNextNucleon()))
00168   {
00169     if(!aNucleon->AreYouHit())
00170     {
00171       resA++;                                                  // Collect A of the ResidNuc
00172       resZ+=G4int (aNucleon->GetDefinition()->GetPDGCharge()); // Collect Z of the ResidNuc
00173     }
00174     else
00175     {
00176       hitMomentum += aNucleon->GetMomentum().vect();           // Sum 3-mom of StringHadr's
00177       hitMass += aNucleon->GetMomentum().m();                  // Sum masses of StringHadrs
00178       hitCount ++;                                             // Calculate STRING hadrons
00179     }
00180   }
00181   G4int targetPDGCode = 90000000 + 1000*resZ + (resA-resZ);    // PDG of theResidualNucleus
00182   G4double targetMass=mNeut;
00183   if (!resZ)                                                   // Nucleus of only neutrons
00184   {
00185     if (resA>1) targetMass*=resA;
00186   }
00187   else targetMass=G4ParticleTable::GetParticleTable()->FindIon(resZ,resA,0,resZ)
00188                                                                             ->GetPDGMass();
00189   G4double targetEnergy = std::sqrt(hitMomentum.mag2()+targetMass*targetMass);
00190   // !! @@ Target should be at rest: hitMomentum=(0,0,0) @@ !! M.K. (go to this system)
00191   G4LorentzVector targ4Mom(-1.*hitMomentum, targetEnergy);
00192   
00193   // Calculate the mean energy lost
00194   std::pair<G4double, G4double> theImpact = theNucleus->RefetchImpactXandY();
00195   G4double impactX = theImpact.first;
00196   G4double impactY = theImpact.second;
00197   G4double impactPar2 = impactX*impactX + impactY*impactY;
00198   G4double radius2 = theNucleus->GetNuclearRadius(theInnerCoreDensityCut*perCent);
00199   //G4double radius2 = theNucleus->GetNuclearRadius(theInnerCoreDensityCut*perCent);
00200   radius2 *= radius2;
00201   G4double pathlength = 0.;
00202 #ifdef ppdebug
00203   G4cout<<"G4StringChipsParticleLevelInterface::Propagate: r="<<std::sqrt(radius2)/fermi
00204         <<", b="<<std::sqrt(impactPar2)/fermi<<", R="<<theNucleus->GetOuterRadius()/fermi
00205         <<", b/r="<<std::sqrt(impactPar2/radius2)<<G4endl; 
00206 #endif
00207   if(radius2 - impactPar2>0) pathlength = 2.*std::sqrt(radius2 - impactPar2);
00208   //pathlength = 27.;                         // *** Only temporary *** Always CHIPS
00209 #ifdef hdebug_SCPLI
00210   toth+=1.;                                 // increment total number of measurements
00211   G4double bfm=std::sqrt(impactPar2)/fermi; // impact parameter
00212   G4double efm=pathlength/fermi;            // energy absorption length
00213   G4int    nbi=static_cast<G4int>(bfm/bhdb);
00214   G4int    nei=static_cast<G4int>(efm/ehde);
00215   if(nbi<nbh) bhis[nbi]++;
00216   else        bover++;
00217   if(nei<nbh) ehis[nei]++;
00218   else        eover++;
00219 #endif
00220   G4double theEnergyLostInFragmentation = theEnergyLossPerFermi*pathlength/fermi;
00221   
00222   // now select all particles in range
00223   std::list<std::pair<G4double, G4KineticTrack *> > theSorted;         // Output
00224   std::list<std::pair<G4double, G4KineticTrack *> >::iterator current; // Input
00225   for(unsigned int secondary = 0; secondary<theSecondaries->size(); secondary++)
00226   {
00227     G4LorentzVector a4Mom = theSecondaries->operator[](secondary)->Get4Momentum();
00228 #ifdef CHIPSdebug
00229     G4cout<<"G4StringChipsParticleLevelInterface::Propagate: ALL STRING particles "
00230           << theSecondaries->operator[](secondary)->GetDefinition()->GetPDGCharge()<<" "
00231           << theSecondaries->operator[](secondary)->GetDefinition()->GetPDGEncoding()<<" "
00232           << a4Mom <<G4endl; 
00233 #endif
00234 #ifdef pdebug
00235     G4cout<<"G4StringChipsParticleLevelInterface::Propagate: in C="
00236           <<theSecondaries->operator[](secondary)->GetDefinition()->GetPDGCharge()<<",PDG="
00237           <<theSecondaries->operator[](secondary)->GetDefinition()->GetPDGEncoding()
00238           <<",4M="<<a4Mom<<", current nS="<<theSorted.size()<<G4endl; 
00239 #endif
00240     G4double toSort = a4Mom.rapidity();          // Rapidity is used for the ordering (?!)
00241     std::pair<G4double, G4KineticTrack *> it;
00242     it.first = toSort;
00243     it.second = theSecondaries->operator[](secondary);
00244     G4bool inserted = false;
00245     for(current = theSorted.begin(); current!=theSorted.end(); current++)
00246     {
00247       if((*current).first > toSort)        // The current is smaller then existing
00248       {
00249         theSorted.insert(current, it);     // It shifts the others up
00250         inserted = true;
00251         break;
00252       }
00253     }
00254     if(!inserted) theSorted.push_back(it); // It is bigger than any previous (@@Check this)
00255   }
00256   
00257   G4LorentzVector proj4Mom(0.,0.,0.,0.);
00258   G4int nD  = 0;
00259   G4int nU  = 0;
00260   G4int nS  = 0;
00261   G4int nAD = 0;
00262   G4int nAU = 0;
00263   G4int nAS = 0;
00264   std::list<std::pair<G4double,G4KineticTrack*> >::iterator firstEscape=theSorted.begin();
00265   G4double runningEnergy = 0;
00266   G4int particleCount = 0;
00267   G4LorentzVector theLow = (*(theSorted.begin())).second->Get4Momentum();
00268   G4LorentzVector theHigh;
00269 
00270 #ifdef CHIPSdebug
00271   G4cout<<"G4StringChipsParticleLevelInterface::Propagate: CHIPS ENERGY LOST "
00272         <<theEnergyLostInFragmentation<<". Sorted rapidities event start"<<G4endl;
00273 #endif
00274 #ifdef pdebug
00275   G4cout<<"G4StringChipsParticleLevelInterface::Propagate: total CHIPS energy = "
00276         <<theEnergyLostInFragmentation<<". Start rapidity sorting nS="<<theSorted.size()
00277         <<G4endl;
00278 #endif
00279 
00280   G4QHadronVector projHV;
00281   std::vector<G4QContent> theContents;
00282   std::vector<G4LorentzVector*> theMomenta;
00283   G4ReactionProductVector* theResult = new G4ReactionProductVector;
00284   G4ReactionProduct* theSec;
00285   G4KineticTrackVector* secondaries;
00286   G4KineticTrackVector* secsec;
00287 #ifdef pdebug
00288   G4cout<<"G4StringChipsParticleLevelInterface::Propagate: Absorption nS="
00289         <<theSorted.size()<<G4endl; 
00290 #endif
00291   G4bool EscapeExists = false;
00292   for(current = theSorted.begin(); current!=theSorted.end(); current++)
00293   {
00294 #ifdef pdebug
00295     G4cout<<"G4StringChipsParticleLevelInterface::Propagate: nq="
00296           <<(*current).second->GetDefinition()->GetQuarkContent(3)<<", naq="
00297           <<(*current).second->GetDefinition()->GetAntiQuarkContent(3)<<", PDG="
00298           <<(*current).second->GetDefinition()->GetPDGEncoding()<<",4M="
00299           <<(*current).second->Get4Momentum()<<G4endl; 
00300 #endif
00301     firstEscape = current;              // Remember to make decays for the rest
00302     G4KineticTrack* aResult = (*current).second;
00303     // This is an old (H.P.) solution, which makes an error in En/Mom conservation
00304     //
00305     // @@ Now it does not include strange particle for the absorption in nuclei (?!) M.K.
00306     //if((*current).second->GetDefinition()->GetQuarkContent(3)!=0 || 
00307     //   (*current).second->GetDefinition()->GetAntiQuarkContent(3) !=0) // Strange quarks
00308     //{
00309     //  G4ParticleDefinition* pdef = aResult->GetDefinition();
00310     //  secondaries = NULL;
00311     //  if ( pdef->GetPDGWidth() > 0 && pdef->GetPDGLifeTime() < 5E-17*s )
00312     //     secondaries = aResult->Decay(); // @@ Decay of only strange resonances (?!) M.K.
00313     //  if ( secondaries == NULL )           // No decay
00314     //  {
00315     //    theSec = new G4ReactionProduct(aResult->GetDefinition());
00316     //    G4LorentzVector current4Mom = aResult->Get4Momentum();
00317     //    current4Mom.boost(targ4Mom.boostVector()); // boost from the targetAtRes system
00318     //    theSec->SetTotalEnergy(current4Mom.t());
00319     //    theSec->SetMomentum(current4Mom.vect());
00320     //    theResult->push_back(theSec);
00321     //  } 
00322     //  else                                 // The decay happened
00323     //  {
00324     //    for (unsigned int aSecondary=0; aSecondary<secondaries->size(); aSecondary++)
00325     //    {
00326     //      theSec = 
00327     //           new G4ReactionProduct(secondaries->operator[](aSecondary)->GetDefinition());
00328     //      G4LorentzVector current4Mom =secondaries->operator[](aSecondary)->Get4Momentum();
00329     //      current4Mom.boost(targ4Mom.boostVector());
00330     //      theSec->SetTotalEnergy(current4Mom.t());
00331     //      theSec->SetMomentum(current4Mom.vect());
00332     //      theResult->push_back(theSec);
00333     //    }
00334     //    std::for_each(secondaries->begin(), secondaries->end(), DeleteKineticTrack());
00335     //    delete secondaries;
00336     //  }
00337     //}
00338     //
00339     //runningEnergy += (*current).second->Get4Momentum().t();
00340     //if((*current).second->GetDefinition() == G4Proton::Proton())
00341     //                                   runningEnergy-=G4Proton::Proton()->GetPDGMass();
00342     //if((*current).second->GetDefinition() == G4Neutron::Neutron())
00343     //                                   runningEnergy-=G4Neutron::Neutron()->GetPDGMass();
00344     //if((*current).second->GetDefinition() == G4Lambda::Lambda())
00345     //                                   runningEnergy-=G4Lambda::Lambda()->GetPDGMass();
00346     //
00347     // New solution starts from here (M.Kossov March 2006) [Strange particles included]
00348     runningEnergy += aResult->Get4Momentum().t();
00349     G4double charge=aResult->GetDefinition()->GetPDGCharge(); // Charge of the particle
00350     G4int strang=aResult->GetDefinition()->GetQuarkContent(3);// Its strangeness
00351     G4int baryn=aResult->GetDefinition()->GetBaryonNumber();  // Its baryon number
00352     if     (baryn>0 && charge>0 && strang<1) runningEnergy-=mProt;  // For positive baryons
00353     else if(baryn>0 && strang<1) runningEnergy-=mNeut;              // For neut/neg baryons
00354     else if(baryn>0) runningEnergy-=mLamb;                          // For strange baryons
00355     else if(baryn<0) runningEnergy+=mProt;                          // For anti-particles
00356     // ------------ End of the new solution
00357 #ifdef CHIPSdebug
00358     G4cout<<"G4StringChipsParticleLevelInterface::Propagate: sorted rapidities "
00359                                     <<(*current).second->Get4Momentum().rapidity()<<G4endl;
00360 #endif
00361 
00362 #ifdef pdebug
00363     G4cout<<"G4StringChipsParticleLevelInterface::Propagate: E="<<runningEnergy<<", EL="
00364                                                     <<theEnergyLostInFragmentation<<G4endl;
00365 #endif
00366     if(runningEnergy > theEnergyLostInFragmentation)
00367     {
00368       EscapeExists = true;
00369       break;
00370     }
00371 #ifdef CHIPSdebug
00372     G4cout <<"G4StringChipsParticleLevelInterface::Propagate: ABSORBED STRING particles "
00373            <<(*current).second->GetDefinition()->GetPDGCharge()<<" "
00374            << (*current).second->GetDefinition()->GetPDGEncoding()<<" "
00375            << (*current).second->Get4Momentum() <<G4endl; 
00376 #endif
00377 #ifdef pdebug
00378     G4cout<<"G4StringChipsParticleLevelInterface::Propagate:C="
00379           <<current->second->GetDefinition()->GetPDGCharge()<<", PDG="
00380           <<current->second->GetDefinition()->GetPDGEncoding()<<", 4M="
00381           <<current->second->Get4Momentum()<<G4endl; 
00382 #endif
00383 
00384     // projectile 4-momentum in target rest frame needed in constructor of QHadron
00385     particleCount++;
00386     theHigh = (*current).second->Get4Momentum(); 
00387     proj4Mom = (*current).second->Get4Momentum(); 
00388     proj4Mom.boost(-1.*targ4Mom.boostVector());   // Back to the system of nucleusAtRest
00389     nD = (*current).second->GetDefinition()->GetQuarkContent(1);
00390     nU = (*current).second->GetDefinition()->GetQuarkContent(2);
00391     nS = (*current).second->GetDefinition()->GetQuarkContent(3);
00392     nAD = (*current).second->GetDefinition()->GetAntiQuarkContent(1);
00393     nAU = (*current).second->GetDefinition()->GetAntiQuarkContent(2);
00394     nAS = (*current).second->GetDefinition()->GetAntiQuarkContent(3);
00395     G4QContent aProjectile(nD, nU, nS, nAD, nAU, nAS);
00396 
00397 #ifdef CHIPSdebug_1
00398     G4cout <<G4endl;
00399     G4cout <<"G4StringChipsParticleLevelInterface::Propagate: Quark content: d="<<nD
00400            <<", u="<<nU<<", s="<<nS<< "Anti-quark content: anit-d="<<nAD<<", anti-u="<<nAU
00401            <<", anti-s="<<nAS<<". G4QContent is constructed"<<endl;
00402 #endif
00403 
00404     theContents.push_back(aProjectile);
00405     G4LorentzVector* aVec = new G4LorentzVector((1./MeV)*proj4Mom); // @@ MeV is basic
00406 
00407 #ifdef CHIPSdebug_1
00408     G4cout<<"G4StringChipsParticleLevelInterface::Propagate: projectile momentum = "
00409           <<*aVec<<G4endl;
00410     G4cout << G4endl;
00411 #endif
00412    
00413     theMomenta.push_back(aVec);
00414   }
00415   std::vector<G4QContent> theFinalContents;
00416   std::vector<G4LorentzVector*> theFinalMomenta;
00417 
00418   // Multiquasmon case: each particle creates a quasmon
00419   //for(unsigned int hp = 0; hp<theContents.size(); hp++)
00420   //{
00421   //  G4QHadron* aHadron = new G4QHadron(theContents[hp], *(theMomenta[hp]) );
00422   //  projHV.push_back(aHadron);
00423   //}
00424   // Energy flow: one Quasmon for each B>0 collection ----------
00425   G4QContent EnFlowQC(0,0,0,0,0,0);
00426   G4LorentzVector EnFlow4M(0.,0.,0.,0.);
00427   //G4bool empty=true;
00428   G4int  barys=0;
00429   G4int  stras=0;
00430   G4int  chars=0;
00431   for(G4int hp = theContents.size()-1; hp>=0; hp--)
00432   {
00433      G4QContent curQC=theContents[hp];
00434      G4int baryn = curQC.GetBaryonNumber();
00435      G4int stran = curQC.GetStrangeness();
00436      G4int charg = curQC.GetCharge();
00437      EnFlowQC += curQC;                        // Keep collecting energy flow
00438      EnFlow4M += *(theMomenta[hp]);
00439      barys += baryn;                           // Collected baryon number
00440      stras += stran;                           // Collected strangeness
00441      chars += charg;                           // Collected charge
00442      //empty = false;
00443   }
00444   if(barys>pcl)                               // Split in two or more parts (to survive!)
00445   {
00446      G4int nprt=(barys-1)/pcl+1;               // Number of parts (pcl=4: 2:5-8,3:9-12...)
00447      G4int curb=barys;
00448      while (nprt>0)
00449      {
00450         nprt--;                                 // One part is going to be created
00451         G4int brnm=pcl;                         // Baryon number of splitting part
00452         curb-=brnm;                             // The residual baryon number
00453         G4double prtM=0.;                       // The resulting GS mass of the part
00454         G4double resM=0.;                       // The resulting GS mass of the residual
00455         G4QContent prtQC(0,0,0,0,0,0);          // The resulting Quark Content of the part
00456         G4int strm=0;                           // Max strangeness per part (stras=0)
00457         if(stras>0) strm=(stras-1)/nprt+1;      // Max strangeness per part (stras>0)
00458         else if(stras<0) strm=(stras+1)/nprt-1; // Max strangeness per part (stras<0)
00459         G4int chgm=0;                           // Max charge per part (chars=0)
00460         if(stras>0) chgm=(chars-1)/nprt+1;      // Max strangeness per part (chars>0)
00461         else if(stras<0) chgm=(chars+1)/nprt-1; // Max strangeness per part (chars<0)
00462         // ---> calculate proposed separated part
00463         //@@ Convert it to a CHIPS function (Which class? G4QH::Conctruct?)
00464         if(!strm)                               // --> The total strangness = 0 (n/p/pi-)
00465         {
00466           if(chgm<0)                            // (n/pi-)
00467           {
00468             prtM=(-chgm)*mPiCh+brnm*mNeut;
00469             prtQC=(-chgm)*PiMiQC+brnm*NeutQC;
00470           }
00471           else                                  // (n/p)
00472           {
00473             prtM=chgm*mProt+(brnm-chgm)*mNeut;
00474             prtQC=chgm*ProtQC+(brnm-chgm)*NeutQC;
00475           }
00476         }
00477         else if(strm>=brnm)                     // ---> BigPositiveStrangeness(L/Pi+/K0/K-)
00478         {
00479           G4int stmb=strm-brnm;
00480           if(chgm<0)                            // (L/K-/K0)
00481           {
00482             prtM=(-chgm)*mKChg+brnm*mLamb+std::abs(stmb+chgm)*mKZer;
00483             prtQC=(-chgm)*KMinQC+brnm*LambQC;
00484             if(stmb>-chgm) prtQC+=(stmb+chgm)*KZerQC;
00485             else if(stmb<-chgm) prtQC+=(-stmb-chgm)*AKZrQC;
00486           }
00487           else                                  // (L/K0/pi+)
00488           {
00489             prtM=chgm*mPiCh+(strm-brnm)*mKZer+brnm*mLamb;
00490             prtQC=chgm*PiPlQC+(strm-brnm)*KZerQC+brnm*LambQC;
00491           }
00492         }
00493         else if(strm>0)                         // ---> PositiveStrangeness<B (L/n/p/Pi+-)
00494         {
00495           G4int bmst=brnm-strm;
00496           if(chgm<0)                            // (L/n/Pi-)
00497           {
00498             prtM=(-chgm)*mPiCh+strm*mLamb+bmst*mNeut;
00499             prtQC=(-chgm)*PiMiQC+strm*LambQC+bmst*NeutQC;
00500           }
00501           else if(chgm>=bmst)                   // (L/p/Pi+)
00502           {
00503             prtM=(chgm-bmst)*mPiCh+strm*mLamb+bmst*mProt;
00504             prtQC=(chgm-bmst)*PiPlQC+strm*LambQC+bmst*ProtQC;
00505           }
00506           else                                  // ch<bmst (L/p/n)
00507           {
00508             prtM=chgm*mProt+strm*mLamb+(bmst-chgm)*mNeut;
00509             prtQC=chgm*ProtQC+strm*LambQC+(bmst-chgm)*NeutQC;
00510           }
00511         }
00512         else                                    // ---> NegativeStrangeness (N/K+/aK0/Pi-)
00513         {
00514           G4int bmst=brnm-strm;
00515           if(chgm>=bmst)                        // (K+/p/Pi+)
00516           {
00517             prtM=(-strm)*mKChg+brnm*mProt+(chgm-bmst)*mPiCh;
00518             prtQC=(-strm)*KPlsQC+brnm*ProtQC+(chgm-bmst)*PiPlQC;
00519           }
00520           else if(chgm>=-strm)                  // (K+/p/n)
00521           {
00522             prtM=(-strm)*mKChg+chgm*mProt+(brnm-chgm)*mNeut;
00523             prtQC=(-strm)*KPlsQC+chgm*ProtQC+(brnm-chgm)*NeutQC;
00524           }
00525           else if(chgm>=0)                      // (K+/aK0/n)
00526           {
00527             prtM=chgm*mKChg+(-chgm-strm)*mKZer+brnm*mNeut;
00528             prtQC=chgm*KPlsQC+(-chgm-strm)*AKZrQC+brnm*NeutQC;
00529           }
00530           else                                  // ch<0 (aK0/n/Pi-)
00531           {
00532             prtM=(-strm)*mKChg+(-chgm)*mPiCh+brnm*mNeut;
00533             prtQC=(-strm)*KPlsQC+(-chgm)*PiMiQC+brnm*NeutQC;
00534           }
00535         }
00536         EnFlowQC-=prtQC;
00537         chgm=chars-chgm;                        // Just to keep the same notation
00538         strm=stras-strm;
00539         brnm=curb;
00540         if(!strm)                               // --> The total strangness = 0 (n/p/pi-)
00541         {
00542           if(chgm<0) resM=(-chgm)*mPiCh+brnm*mNeut;
00543           else       resM=chgm*mProt+(brnm-chgm)*mNeut;
00544         }
00545         else if(strm>=brnm)                     // ---> BigPositiveStrangeness(L/Pi+/K0/K-)
00546         {
00547           G4int stmb=strm-brnm;
00548           if(chgm<0) resM=(-chgm)*mKChg+brnm*mLamb+std::abs(stmb+chgm)*mKZer;
00549           else       resM=chgm*mPiCh+(strm-brnm)*mKZer+brnm*mLamb;
00550         }
00551         else if(strm>0)                         // ---> PositiveStrangeness<B (L/n/p/Pi+-)
00552         {
00553           G4int bmst=brnm-strm;
00554           if     (chgm<0)     resM=(-chgm)*mPiCh+strm*mLamb+bmst*mNeut;
00555           else if(chgm>=bmst) resM=(chgm-bmst)*mPiCh+strm*mLamb+bmst*mProt;
00556           else                resM=chgm*mProt+strm*mLamb+(bmst-chgm)*mNeut;
00557         }
00558         else                                    // ---> NegativeStrangeness (N/K+/aK0/Pi-)
00559         {
00560           G4int bmst=brnm-strm;
00561           if     (chgm>=bmst)  resM=(-strm)*mKChg+brnm*mProt+(chgm-bmst)*mPiCh;
00562           else if(chgm>=-strm) resM=(-strm)*mKChg+chgm*mProt+(brnm-chgm)*mNeut;
00563           else if(chgm>=0)     resM=chgm*mKChg+(-chgm-strm)*mKZer+brnm*mNeut;
00564           else                 resM=(-strm)*mKChg+(-chgm)*mPiCh+brnm*mNeut;
00565         }
00566         G4LorentzVector prt4M=(prtM/(prtM+resM))*EnFlow4M;
00567         EnFlow4M-=prt4M;
00568         EnFlowQC-=prtQC;
00569         G4QHadron* aHadron = new G4QHadron(prtQC, prt4M);
00570         projHV.push_back(aHadron);
00571         if(nprt==1)
00572         {
00573           G4QHadron* fHadron = new G4QHadron(EnFlowQC, EnFlow4M);
00574           projHV.push_back(fHadron);
00575           nprt=0;
00576         }
00577 #ifdef debug
00578         G4cout<<"G4StringChipsParticleLevelInterface::Propagate: nprt="<<nprt<<G4endl;
00579 #endif
00580      } // End of WHILE
00581   }
00582   else
00583   {
00584      G4QHadron* aHadron = new G4QHadron(EnFlowQC, EnFlow4M);
00585      projHV.push_back(aHadron);
00586   }
00587 
00588   // construct the quasmon
00589   size_t i;
00590   for (i=0; i<theFinalMomenta.size(); i++) delete theFinalMomenta[i];
00591   for (i=0; i<theMomenta.size();      i++) delete theMomenta[i];
00592   theFinalMomenta.clear();
00593   theMomenta.clear();
00594 
00595   G4QNucleus::SetParameters(fractionOfSingleQuasiFreeNucleons,
00596                             fractionOfPairedQuasiFreeNucleons,
00597                             clusteringCoefficient, 
00598                             fusionToExchange);
00599   G4Quasmon::SetParameters(temperature, halfTheStrangenessOfSee, etaToEtaPrime);
00600 
00601 #ifdef CHIPSdebug
00602   G4cout<<"G4StringChipsParticleLevelInterface::Propagate: G4QNucleus parameters "
00603         <<fractionOfSingleQuasiFreeNucleons<<" "<<fractionOfPairedQuasiFreeNucleons
00604         <<" "<<clusteringCoefficient<<G4endl;
00605   G4cout<<"G4Quasmon parameters "<<temperature<<" "<<halfTheStrangenessOfSee<<" "
00606         <<etaToEtaPrime << G4endl;
00607   G4cout<<"The Target PDG code = "<<targetPDGCode<<G4endl;
00608   G4cout<<"The projectile momentum = "<<1./MeV*proj4Mom<<G4endl;
00609   G4cout<<"The target momentum = "<<1./MeV*targ4Mom<<G4endl;
00610 #endif
00611 
00612   // now call chips with this info in place
00613   G4QHadronVector* output = 0;
00614   if (particleCount!=0 && resA!=0)
00615   {
00616     //  G4QCHIPSWorld aWorld(nop);              // Create CHIPS World of nop particles
00617     G4QCHIPSWorld::Get()->GetParticles(nop);
00618     G4QEnvironment* pan= new G4QEnvironment(projHV, targetPDGCode);
00619 #ifdef pdebug
00620       G4cout<<"G4StringChipsParticleLevelInterface::Propagate: CHIPS fragmentation, rA="
00621             <<resA<<", #AbsPt="<<particleCount<<G4endl; 
00622 #endif
00623     try
00624     {
00625       output = pan->Fragment();                 // The main fragmentation member function
00626     }
00627     catch(G4HadronicException& aR)
00628     {
00629       G4cerr << "Exception thrown of G4StringChipsParticleLevelInterface "<<G4endl;
00630       G4cerr << " targetPDGCode = "<< targetPDGCode <<G4endl;
00631       G4cerr << " The projectile momentum = "<<1./MeV*proj4Mom<<G4endl;
00632       G4cerr << " The target momentum = "<<1./MeV*targ4Mom<<G4endl<<G4endl;
00633       G4cerr << " Dumping the information in the pojectile list"<<G4endl;
00634       for(i=0; i< projHV.size(); i++)
00635       {
00636         G4cerr <<"  Incoming 4-momentum and PDG code of "<<i<<"'th hadron: "
00637         <<" "<< projHV[i]->Get4Momentum()<<" "<<projHV[i]->GetPDGCode()<<G4endl;
00638       }
00639       throw;
00640     }
00641     delete pan;
00642   }
00643   else
00644   {
00645 #ifdef pdebug
00646     G4cout<<"G4StringChipsParticleLevelInterface::Propagate: NO CHIPS fragmentation, rA="
00647           <<resA<<", #AbsPt="<<particleCount<<G4endl; 
00648 #endif
00649     output = new G4QHadronVector;
00650   }   
00651 
00652   // clean up impinging particles
00653   std::for_each(projHV.begin(), projHV.end(), DeleteQHadron());
00654   projHV.clear();
00655    
00656   // Fill the result.
00657 #ifdef CHIPSdebug
00658   G4cout << "NEXT EVENT, EscapeExists="<<EscapeExists<<G4endl;
00659 #endif
00660 
00661   // first decay and add all escaping particles.
00662   if (EscapeExists) for (current = firstEscape; current!=theSorted.end(); current++)
00663   {
00664     G4KineticTrack* aResult = (*current).second;
00665     G4ParticleDefinition* pdef=aResult->GetDefinition();
00666     secondaries = NULL;
00667     //if(pdef->GetPDGWidth() > 0 && pdef->GetPDGLifeTime() < 5E-17*s ) // HPW version
00668     if ( pdef->IsShortLived() )
00669     {
00670       secondaries = aResult->Decay();
00671       for (unsigned int aSecondary=0; aSecondary<secondaries->size(); aSecondary++)
00672       {
00673         G4KineticTrack* bResult=secondaries->operator[](aSecondary);
00674         G4ParticleDefinition* sdef=bResult->GetDefinition();
00675         if ( sdef->IsShortLived() )
00676         {
00677           secsec = bResult->Decay();
00678           for (unsigned int bSecondary=0; bSecondary<secsec->size(); bSecondary++)
00679           {
00680              G4KineticTrack* cResult=secsec->operator[](bSecondary);
00681              G4ParticleDefinition* cdef=cResult->GetDefinition();
00682              theSec = new G4ReactionProduct(cdef);
00683              G4LorentzVector cur4Mom = cResult->Get4Momentum();
00684              cur4Mom.boost(targ4Mom.boostVector());
00685              theSec->SetTotalEnergy(cur4Mom.t());
00686              theSec->SetMomentum(cur4Mom.vect());
00687 #ifdef trapdebug
00688              if(cdef->GetPDGEncoding()==113) G4cout
00689                 <<"G4StringChipsParticleLevelInterface::Propagate: *Rho0* QGS dec2 PDG="
00690                 <<cdef->GetPDGEncoding()<<",4M="<<cur4Mom<<", grandparPDG= "
00691                 <<pdef->GetPDGEncoding()<<", parPDG= "<<sdef->GetPDGEncoding()<<G4endl; 
00692 #endif
00693 #ifdef pdebug
00694              G4cout<<"G4StringChipsParticleLevelInterface::Propagate: *OUT* QGS dec2 PDG="
00695               <<sdef->GetPDGEncoding()<<",4M="<<cur4Mom<<G4endl; 
00696 #endif
00697              theResult->push_back(theSec);
00698           }
00699           std::for_each(secsec->begin(), secsec->end(), DeleteKineticTrack());
00700           delete secsec;
00701         }
00702         else
00703         {
00704           theSec = new G4ReactionProduct(sdef);
00705           G4LorentzVector current4Mom = bResult->Get4Momentum();
00706           current4Mom.boost(targ4Mom.boostVector());
00707           theSec->SetTotalEnergy(current4Mom.t());
00708           theSec->SetMomentum(current4Mom.vect());
00709 #ifdef trapdebug
00710           if(sdef->GetPDGEncoding()==113)
00711             G4cout<<"G4StringChipsParticleLevelInterface::Propagate:*Rho0* QGS decay PDG="
00712                   <<sdef->GetPDGEncoding()<<",4M="<<current4Mom<<", parentPDG= "
00713                   <<pdef->GetPDGEncoding()<<G4endl; 
00714         //throw G4HadronicException(__FILE__,__LINE__,
00715         //                          "G4StringChipsParticleLevelInterface: Rho0 is found!");
00716 #endif
00717 #ifdef pdebug
00718           G4cout<<"G4StringChipsParticleLevelInterface::Propagate: *OUT* QGS decay PDG="
00719                 <<sdef->GetPDGEncoding()<<",4M="<<current4Mom<<G4endl; 
00720 #endif
00721           theResult->push_back(theSec);
00722         }
00723       }
00724       std::for_each(secondaries->begin(), secondaries->end(), DeleteKineticTrack());
00725       delete secondaries;
00726     }
00727     else
00728     {
00729       theSec = new G4ReactionProduct(aResult->GetDefinition());
00730       G4LorentzVector current4Mom = aResult->Get4Momentum();
00731       current4Mom.boost(targ4Mom.boostVector());
00732       theSec->SetTotalEnergy(current4Mom.t());
00733       theSec->SetMomentum(current4Mom.vect());
00734 #ifdef trapdebug
00735       if(aResult->GetDefinition()->GetPDGEncoding()==113)
00736         G4cout<<"G4StringChipsParticleLevelInterface::Propagate: *OUT* QGS stable PDG="
00737               <<aResult->GetDefinition()->GetPDGEncoding()<<",4M="<<current4Mom<<G4endl; 
00738 #endif
00739 #ifdef pdebug
00740       G4cout<<"G4StringChipsParticleLevelInterface::Propagate: *OUT* QGS stable PDG="
00741             <<aResult->GetDefinition()->GetPDGEncoding()<<",4M="<<current4Mom<<G4endl; 
00742 #endif
00743       theResult->push_back(theSec);
00744     } 
00745   }
00746   std::for_each(theSecondaries->begin(), theSecondaries->end(), DeleteKineticTrack());
00747   delete theSecondaries;
00748     
00749   // now add the quasmon output
00750   G4int maxParticle=output->size();
00751 #ifdef CHIPSdebug
00752   G4cout << "Number of particles from string"<<theResult->size()<<G4endl;
00753   G4cout << "Number of particles from chips"<<maxParticle<<G4endl;
00754 #endif
00755 #ifdef pdebug
00756   G4cout << "Number of particles from QGS="<<theResult->size()<<G4endl;
00757   G4cout << "Number of particles from CHIPS="<<maxParticle<<G4endl;
00758 #endif
00759   if(maxParticle) for(G4int particle = 0; particle < maxParticle; particle++)
00760   {
00761     if(output->operator[](particle)->GetNFragments() != 0) 
00762     {
00763       delete output->operator[](particle);
00764       continue;
00765     }
00766     G4int pdgCode = output->operator[](particle)->GetPDGCode();
00767 
00768 
00769 #ifdef CHIPSdebug
00770     G4cerr << "PDG code of chips particle = "<<pdgCode<<G4endl;
00771 #endif
00772 
00773     G4ParticleDefinition * theDefinition;
00774     // Note that I still have to take care of strange nuclei
00775     // For this I need the mass calculation, and a changed interface
00776     // for ion-table ==> work for Hisaya @@@@@@@
00777     // Then I can sort out the pdgCode. I also need a decau process 
00778     // for strange nuclei; may be another chips interface
00779     if(pdgCode>90000000) 
00780     {
00781       G4int aZ = (pdgCode-90000000)/1000;
00782       if (aZ>1000) aZ=aZ%1000;  // patch for strange nuclei, to be repaired @@@@
00783       G4int anN = pdgCode-90000000-1000*aZ;
00784       if(anN>1000) anN=anN%1000; // patch for strange nuclei, to be repaired @@@@
00785 
00786       if(pdgCode==90000999) theDefinition = G4PionPlus::PionPlusDefinition();
00787       else if(pdgCode==89999001) theDefinition = G4PionMinus::PionMinusDefinition();
00788       else if(pdgCode==90999999) theDefinition = G4KaonZero::KaonZeroDefinition();
00789       else if(pdgCode==90999000) theDefinition = G4KaonMinus::KaonMinusDefinition();
00790       else if(pdgCode==89001000) theDefinition = G4KaonPlus::KaonPlusDefinition();
00791       else if(pdgCode==89000001) theDefinition = G4AntiKaonZero::AntiKaonZeroDefinition();
00792       else if(pdgCode==91000000) theDefinition = G4Lambda::LambdaDefinition();
00793       else if(pdgCode==92000000) theDefinition = G4Lambda::LambdaDefinition(); //NLambd?
00794       else if(pdgCode==93000000) theDefinition = G4Lambda::LambdaDefinition();
00795       else if(pdgCode==94000000) theDefinition = G4Lambda::LambdaDefinition();
00796       else if(pdgCode==95000000) theDefinition = G4Lambda::LambdaDefinition();
00797       else if(pdgCode==96000000) theDefinition = G4Lambda::LambdaDefinition();
00798       else if(pdgCode==97000000) theDefinition = G4Lambda::LambdaDefinition();
00799       else if(pdgCode==98000000) theDefinition = G4Lambda::LambdaDefinition();
00800       else if(aZ == 0 && anN == 1) theDefinition = G4Neutron::Neutron();
00801       else theDefinition = G4ParticleTable::GetParticleTable()->FindIon(aZ,anN+aZ,0,aZ);
00802     }    
00803     else theDefinition = G4ParticleTable::GetParticleTable()->FindParticle(pdgCode);
00804 
00805 #ifdef CHIPSdebug
00806     G4cout << "Particle code produced = "<< pdgCode <<G4endl;
00807 #endif
00808 
00809     if(theDefinition)
00810     {
00811       theSec = new G4ReactionProduct(theDefinition);
00812       G4LorentzVector current4Mom = output->operator[](particle)->Get4Momentum();
00813       current4Mom.boost(targ4Mom.boostVector());
00814       theSec->SetTotalEnergy(current4Mom.t());
00815       theSec->SetMomentum(current4Mom.vect());
00816 #ifdef pdebug
00817       G4cout<<"G4StringChipsParticleLevelInterface::Propagate: *OUT* CHIPS PDG="
00818               <<theDefinition->GetPDGEncoding()<<",4M="<<current4Mom<<G4endl; 
00819 #endif
00820       theResult->push_back(theSec);
00821     }
00822 #ifdef pdebug
00823     else
00824     {
00825       G4cerr <<"G4StringChipsParticleLevelInterface::Propagate: WARNING"<<G4endl;
00826       G4cerr << "Getting unknown pdgCode from chips in ParticleLevelInterface"<<G4endl;
00827       G4cerr << "skipping particle with pdgCode = "<<pdgCode<<G4endl<<G4endl;
00828     }
00829 #endif
00830     
00831 #ifdef CHIPSdebug
00832     G4cout <<"CHIPS particles "<<theDefinition->GetPDGCharge()<<" "
00833            << theDefinition->GetPDGEncoding()<<" "
00834     << output->operator[](particle)->Get4Momentum() <<G4endl; 
00835 #endif
00836 
00837     delete output->operator[](particle);
00838   }
00839   else
00840   {
00841     if(resA>0)
00842     {
00843       G4ParticleDefinition* theDefinition = G4Neutron::Neutron();
00844       if(resA==1) // The residual nucleus at rest must be added to conserve BaryN & Charge
00845       {
00846         if(resZ == 1) theDefinition = G4Proton::Proton();
00847       }
00848       else theDefinition = G4ParticleTable::GetParticleTable()->FindIon(resZ,resA,0,resZ);
00849       theSec = new G4ReactionProduct(theDefinition);
00850       theSec->SetTotalEnergy(theDefinition->GetPDGMass());
00851       theSec->SetMomentum(G4ThreeVector(0.,0.,0.));
00852       theResult->push_back(theSec);
00853       if(!resZ && resA>0) for(G4int ni=1; ni<resA; ni++) 
00854       {
00855         theSec = new G4ReactionProduct(theDefinition);
00856         theSec->SetTotalEnergy(theDefinition->GetPDGMass());
00857         theSec->SetMomentum(G4ThreeVector(0.,0.,0.));
00858         theResult->push_back(theSec);
00859       }
00860     }
00861   }
00862   delete output;
00863 
00864 #ifdef CHIPSdebug
00865   G4cout << "Number of particles"<<theResult->size()<<G4endl;
00866   G4cout << G4endl;
00867   G4cout << "QUASMON preparation info "
00868          << 1./MeV*proj4Mom<<" "
00869   << 1./MeV*targ4Mom<<" "
00870   << nD<<" "<<nU<<" "<<nS<<" "<<nAD<<" "<<nAU<<" "<<nAS<<" "
00871   << hitCount<<" "
00872   << particleCount<<" "
00873   << theLow<<" "
00874   << theHigh<<" "
00875   << G4endl;
00876 #endif
00877 
00878   return theResult;
00879 } 


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