#include <G4StringChipsParticleLevelInterface.hh>
Inheritance diagram for G4StringChipsParticleLevelInterface:
Public Member Functions | |
G4StringChipsParticleLevelInterface () | |
virtual G4HadFinalState * | ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &theNucleus) |
virtual G4ReactionProductVector * | Propagate (G4KineticTrackVector *theSecondaries, G4V3DNucleus *theNucleus) |
Definition at line 38 of file G4StringChipsParticleLevelInterface.hh.
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 }
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 }