00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 #include <utility>
00033 #include <list>
00034
00035 #include "G4StringChipsInterface.hh"
00036 #include "globals.hh"
00037 #include "G4SystemOfUnits.hh"
00038 #include "G4KineticTrackVector.hh"
00039 #include "G4Nucleon.hh"
00040 #include "G4LorentzRotation.hh"
00041
00042 G4StringChipsInterface::G4StringChipsInterface()
00043 {
00044 #ifdef CHIPSdebug
00045 G4cout << "Please enter the energy loss per fermi in GeV"<<G4endl;
00046 G4cin >> theEnergyLossPerFermi;
00047 #endif
00048 #ifdef debug
00049 G4cout<<"G4StringChipsInterface::Constructor is called"<<G4endl;
00050 #endif
00051 theEnergyLossPerFermi = 0.5*GeV;
00052
00053 }
00054
00055 G4HadFinalState* G4StringChipsInterface::
00056 ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& theNucleus)
00057 {
00058 #ifdef debug
00059 G4cout<<"G4StringChipsInterface::ApplyYourself is called"<<G4endl;
00060 #endif
00061 return theModel.ApplyYourself(aTrack, theNucleus);
00062 }
00063
00064 G4ReactionProductVector* G4StringChipsInterface::
00065 Propagate(G4KineticTrackVector* theSecondaries, G4V3DNucleus* theNucleus)
00066 {
00067 #ifdef debug
00068 G4cout<<"G4StringChipsInterface::Propagate is called"<<G4endl;
00069 #endif
00070
00071
00072 if(theSecondaries->size() == 1)
00073 throw G4HadronicException(__FILE__, __LINE__, "G4StringChipsInterface: Only one particle from String models!");
00074
00075
00076 std::pair<G4double, G4double> theImpact = theNucleus->RefetchImpactXandY();
00077 G4double impactX = theImpact.first;
00078 G4double impactY = theImpact.second;
00079 G4double inpactPar2 = impactX*impactX + impactY*impactY;
00080
00081 G4double radius2 = theNucleus->GetNuclearRadius(5*perCent);
00082 radius2 *= radius2;
00083 G4double pathlength = 0;
00084 if(radius2 - inpactPar2>0) pathlength = 2.*std::sqrt(radius2 - inpactPar2);
00085 G4double theEnergyLostInFragmentation = theEnergyLossPerFermi*pathlength/fermi;
00086
00087
00088 std::list<std::pair<G4double, G4KineticTrack *> > theSorted;
00089 std::list<std::pair<G4double, G4KineticTrack *> >::iterator current;
00090 for(unsigned int secondary = 0; secondary<theSecondaries->size(); secondary++)
00091 {
00092 G4LorentzVector a4Mom = theSecondaries->operator[](secondary)->Get4Momentum();
00093
00094 #ifdef CHIPSdebug
00095 G4cout <<"ALL STRING particles "<<theSecondaries->operator[](secondary)->GetDefinition()->GetPDGCharge()<<" "
00096 << theSecondaries->operator[](secondary)->GetDefinition()->GetPDGEncoding()<<" "
00097 << a4Mom <<G4endl;
00098 #endif
00099
00100 G4double toSort = a4Mom.rapidity();
00101 std::pair<G4double, G4KineticTrack *> it;
00102 it.first = toSort;
00103 it.second = theSecondaries->operator[](secondary);
00104 G4bool inserted = false;
00105 for(current = theSorted.begin(); current!=theSorted.end(); current++)
00106 {
00107 if((*current).first > toSort)
00108 {
00109 theSorted.insert(current, it);
00110 inserted = true;
00111 break;
00112 }
00113 }
00114 if(!inserted)
00115 {
00116 theSorted.push_front(it);
00117 }
00118 }
00119
00120 G4LorentzVector proj4Mom(0.,0.,0.,0.);
00121
00122
00123
00124
00125
00126
00127 G4int nD = 0;
00128 G4int nU = 0;
00129 G4int nS = 0;
00130 G4int nAD = 0;
00131 G4int nAU = 0;
00132 G4int nAS = 0;
00133 std::list<std::pair<G4double, G4KineticTrack *> >::iterator firstEscaping = theSorted.begin();
00134 G4double runningEnergy = 0;
00135 G4int particleCount = 0;
00136 G4LorentzVector theLow = (*(theSorted.begin())).second->Get4Momentum();
00137 G4LorentzVector theHigh;
00138
00139 #ifdef CHIPSdebug
00140 G4cout << "CHIPS ENERGY LOST "<<theEnergyLostInFragmentation<<G4endl;
00141 G4cout << "sorted rapidities event start"<<G4endl;
00142 #endif
00143
00144 G4ReactionProductVector * theResult = new G4ReactionProductVector;
00145 G4ReactionProduct * theSec;
00146 G4KineticTrackVector * secondaries;
00147 #ifdef pdebug
00148 G4cout<<"G4StringChipsInterface::Propagate: Absorption"<<G4endl;
00149 #endif
00150
00151
00152 for(current = theSorted.begin(); current!=theSorted.end(); current++)
00153 {
00154 firstEscaping = current;
00155 if((*current).second->GetDefinition()->GetQuarkContent(3)!=0 ||
00156 (*current).second->GetDefinition()->GetAntiQuarkContent(3) !=0)
00157 {
00158 G4KineticTrack * aResult = (*current).second;
00159 G4ParticleDefinition * pdef=aResult->GetDefinition();
00160 secondaries = NULL;
00161 if ( pdef->GetPDGWidth() > 0 && pdef->GetPDGLifeTime() < 5E-17*s )
00162 {
00163 secondaries = aResult->Decay();
00164 }
00165 if ( secondaries == NULL )
00166 {
00167 theSec = new G4ReactionProduct(aResult->GetDefinition());
00168 G4LorentzVector current4Mom = aResult->Get4Momentum();
00169 theSec->SetTotalEnergy(current4Mom.t());
00170 theSec->SetMomentum(current4Mom.vect());
00171 theResult->push_back(theSec);
00172 }
00173 else
00174 {
00175 for (unsigned int aSecondary=0; aSecondary<secondaries->size(); aSecondary++)
00176 {
00177 theSec = new G4ReactionProduct(secondaries->operator[](aSecondary)->GetDefinition());
00178 G4LorentzVector current4Mom = secondaries->operator[](aSecondary)->Get4Momentum();
00179 theSec->SetTotalEnergy(current4Mom.t());
00180 theSec->SetMomentum(current4Mom.vect());
00181 theResult->push_back(theSec);
00182 }
00183 std::for_each(secondaries->begin(), secondaries->end(), DeleteKineticTrack());
00184 delete secondaries;
00185 }
00186 continue;
00187 }
00188 runningEnergy += (*current).second->Get4Momentum().t();
00189 if(runningEnergy > theEnergyLostInFragmentation) break;
00190
00191 #ifdef CHIPSdebug
00192 G4cout <<"ABSORBED STRING particles "<<current->second->GetDefinition()->GetPDGCharge()<<" "
00193 << current->second->GetDefinition()->GetPDGEncoding()<<" "
00194 << current->second->Get4Momentum() <<G4endl;
00195 #endif
00196 #ifdef pdebug
00197 G4cout<<"G4StringChipsInterface::Propagate:C="
00198 <<current->second->GetDefinition()->GetPDGCharge()<<", PDG="
00199 << current->second->GetDefinition()->GetPDGEncoding()<<", 4M="
00200 << current->second->Get4Momentum() <<G4endl;
00201 #endif
00202
00203
00204 particleCount++;
00205 theHigh = (*current).second->Get4Momentum();
00206 proj4Mom += (*current).second->Get4Momentum();
00207
00208 #ifdef CHIPSdebug
00209 G4cout << "sorted rapidities "<<current->second->Get4Momentum().rapidity()<<G4endl;
00210 #endif
00211
00212
00213 nD += (*current).second->GetDefinition()->GetQuarkContent(1);
00214 nU += (*current).second->GetDefinition()->GetQuarkContent(2);
00215 nS += (*current).second->GetDefinition()->GetQuarkContent(3);
00216 nAD += (*current).second->GetDefinition()->GetAntiQuarkContent(1);
00217 nAU += (*current).second->GetDefinition()->GetAntiQuarkContent(2);
00218 nAS += (*current).second->GetDefinition()->GetAntiQuarkContent(3);
00219 }
00220
00221
00222 #ifdef CHIPSdebug
00223 G4cout << "Quark content: d="<<nD<<", u="<<nU<< ", s="<< nS << G4endl;
00224 G4cout << "Anti-quark content: anit-d="<<nAD<<", anti-u="<<nAU<< ", anti-s="<< nAS << G4endl;
00225 #endif
00226
00227 G4QContent theProjectiles(nD, nU, nS, nAD, nAU, nAS);
00228
00229 #ifdef CHIPSdebug
00230 G4cout << "G4QContent is constructed"<<endl;
00231 #endif
00232
00233
00234
00235 theNucleus->StartLoop();
00236 G4Nucleon * aNucleon;
00237 G4int resA = 0;
00238 G4int resZ = 0;
00239 G4ThreeVector hitMomentum(0,0,0);
00240 G4double hitMass = 0;
00241 G4int hitCount = 0;
00242 while((aNucleon = theNucleus->GetNextNucleon()))
00243 {
00244 if(!aNucleon->AreYouHit())
00245 {
00246 resA++;
00247 resZ+=G4int (aNucleon->GetDefinition()->GetPDGCharge());
00248 }
00249 else
00250 {
00251 hitMomentum += aNucleon->GetMomentum().vect();
00252 hitMass += aNucleon->GetMomentum().m();
00253 hitCount ++;
00254 }
00255 }
00256 G4int targetPDGCode = 90000000 + 1000*resZ + (resA-resZ);
00257 G4double targetMass = theNucleus->GetMass();
00258 targetMass -= hitMass;
00259 G4double targetEnergy = std::sqrt(hitMomentum.mag2()+targetMass*targetMass);
00260
00261 G4LorentzVector targ4Mom(-1.*hitMomentum, targetEnergy);
00262
00263
00264 G4int nop = 85;
00265
00266
00267 G4double fractionOfSingleQuasiFreeNucleons = 0.5;
00268 G4double fractionOfPairedQuasiFreeNucleons = 0.05;
00269 G4double clusteringCoefficient = 5.;
00270 G4double temperature = 180.;
00271 G4double halfTheStrangenessOfSee = 0.3;
00272 G4double etaToEtaPrime = 0.3;
00273
00274 G4QNucleus::SetParameters(fractionOfSingleQuasiFreeNucleons,
00275 fractionOfPairedQuasiFreeNucleons,
00276 clusteringCoefficient);
00277 G4Quasmon::SetParameters(temperature,
00278 halfTheStrangenessOfSee,
00279 etaToEtaPrime);
00280
00281 #ifdef CHIPSdebug
00282 G4cout << "G4QNucleus parameters "<< fractionOfSingleQuasiFreeNucleons << " "
00283 << fractionOfPairedQuasiFreeNucleons << " "<< clusteringCoefficient << G4endl;
00284 G4cout << "G4Quasmon parameters "<< temperature << " "<< halfTheStrangenessOfSee << " "
00285 <<etaToEtaPrime << G4endl;
00286 G4cout << "The Target PDG code = "<<targetPDGCode<<G4endl;
00287 G4cout << "The projectile momentum = "<<1./MeV*proj4Mom<<G4endl;
00288 G4cout << "The target momentum = "<<1./MeV*targ4Mom<<G4endl;
00289 #endif
00290
00291
00292
00293
00294 G4QCHIPSWorld::Get()->GetParticles(nop);
00295 G4QHadronVector projHV;
00296
00297 proj4Mom.boost(-1.*targ4Mom.boostVector());
00298
00299 G4LorentzRotation toZ;
00300 toZ.rotateZ(-1*proj4Mom.phi());
00301 toZ.rotateY(-1*proj4Mom.theta());
00302 proj4Mom = toZ*proj4Mom;
00303 G4LorentzRotation toLab(toZ.inverse());
00304
00305 #ifdef CHIPSdebug
00306 G4cout << "a Boosted projectile vector along z"<<proj4Mom<<" "<<proj4Mom.mag()<<G4endl;
00307 #endif
00308
00309 G4QHadron* iH = new G4QHadron(theProjectiles, 1./MeV*proj4Mom);
00310 projHV.push_back(iH);
00311
00312
00313 G4QHadronVector * output = 0;
00314 if (particleCount!=0)
00315 {
00316 G4QEnvironment* pan= new G4QEnvironment(projHV, targetPDGCode);
00317 try
00318 {
00319 output = pan->Fragment();
00320 }
00321 catch(G4HadronicException & aR)
00322 {
00323 G4cerr << "Exception thrown passing through G4ChiralInvariantPhaseSpace "<<G4endl;
00324 G4cerr << " targetPDGCode = "<< targetPDGCode <<G4endl;
00325 G4cerr << " Dumping the information in the pojectile list"<<G4endl;
00326 for(size_t i=0; i< projHV.size(); i++)
00327 {
00328 G4cerr <<" Incoming 4-momentum and PDG code of "<<i<<"'th hadron: "
00329 <<" "<< projHV[i]->Get4Momentum()<<" "<<projHV[i]->GetPDGCode()<<G4endl;
00330 }
00331 throw;
00332 }
00333 std::for_each(projHV.begin(), projHV.end(), DeleteQHadron());
00334 projHV.clear();
00335 delete pan;
00336 }
00337 else output = new G4QHadronVector;
00338
00339
00340 #ifdef CHIPSdebug
00341 G4cout << "NEXT EVENT"<<endl;
00342 #endif
00343
00344
00345 for(current = firstEscaping; current!=theSorted.end(); current++)
00346 {
00347 G4KineticTrack * aResult = (*current).second;
00348 G4ParticleDefinition * pdef=aResult->GetDefinition();
00349 secondaries = NULL;
00350 if ( pdef->GetPDGWidth() > 0 && pdef->GetPDGLifeTime() < 5E-17*s )
00351 {
00352 secondaries = aResult->Decay();
00353 }
00354 if ( secondaries == NULL )
00355 {
00356 theSec = new G4ReactionProduct(aResult->GetDefinition());
00357 G4LorentzVector current4Mom = aResult->Get4Momentum();
00358 current4Mom = toLab*current4Mom;
00359 current4Mom.boost(targ4Mom.boostVector());
00360 theSec->SetTotalEnergy(current4Mom.t());
00361 theSec->SetMomentum(current4Mom.vect());
00362 theResult->push_back(theSec);
00363 }
00364 else
00365 {
00366 for (unsigned int aSecondary=0; aSecondary<secondaries->size(); aSecondary++)
00367 {
00368 theSec = new G4ReactionProduct(secondaries->operator[](aSecondary)->GetDefinition());
00369 G4LorentzVector current4Mom = secondaries->operator[](aSecondary)->Get4Momentum();
00370 current4Mom = toLab*current4Mom;
00371 current4Mom.boost(targ4Mom.boostVector());
00372 theSec->SetTotalEnergy(current4Mom.t());
00373 theSec->SetMomentum(current4Mom.vect());
00374 theResult->push_back(theSec);
00375 }
00376 std::for_each(secondaries->begin(), secondaries->end(), DeleteKineticTrack());
00377 delete secondaries;
00378 }
00379 }
00380 std::for_each(theSecondaries->begin(), theSecondaries->end(), DeleteKineticTrack());
00381 delete theSecondaries;
00382
00383
00384 #ifdef CHIPSdebug
00385 G4cout << "Number of particles from string"<<theResult->size()<<G4endl;
00386 G4cout << "Number of particles from chips"<<output->size()<<G4endl;
00387 #endif
00388
00389 for(unsigned int particle = 0; particle < output->size(); particle++)
00390 {
00391 if(output->operator[](particle)->GetNFragments() != 0)
00392 {
00393 delete output->operator[](particle);
00394 continue;
00395 }
00396
00397 G4int pdgCode = output->operator[](particle)->GetPDGCode();
00398 G4ParticleDefinition * theDefinition;
00399
00400
00401
00402
00403
00404 if(pdgCode>90000000)
00405 {
00406 G4int aZ = (pdgCode-90000000)/1000;
00407 if (aZ>1000) aZ=aZ%1000;
00408 G4int anN = pdgCode-90000000-1000*aZ;
00409 if(anN>1000) anN=anN%1000;
00410 if(pdgCode==91000000) theDefinition = G4Lambda::LambdaDefinition();
00411 else if(pdgCode==92000000) theDefinition = G4Lambda::LambdaDefinition();
00412 else if(pdgCode==93000000) theDefinition = G4Lambda::LambdaDefinition();
00413 else if(pdgCode==94000000) theDefinition = G4Lambda::LambdaDefinition();
00414 else if(pdgCode==95000000) theDefinition = G4Lambda::LambdaDefinition();
00415 else if(pdgCode==96000000) theDefinition = G4Lambda::LambdaDefinition();
00416 else if(pdgCode==97000000) theDefinition = G4Lambda::LambdaDefinition();
00417 else if(pdgCode==98000000) theDefinition = G4Lambda::LambdaDefinition();
00418 else if(aZ == 0 && anN == 1) theDefinition = G4Neutron::Neutron();
00419 else theDefinition = G4ParticleTable::GetParticleTable()->FindIon(aZ,anN+aZ,0,aZ);
00420 }
00421 else
00422 {
00423 theDefinition = G4ParticleTable::GetParticleTable()->FindParticle(output->operator[](particle)->GetPDGCode());
00424 }
00425 #ifdef CHIPSdebug
00426 G4cout << "Particle code produced = "<< pdgCode <<G4endl;
00427 #endif
00428
00429 theSec = new G4ReactionProduct(theDefinition);
00430 G4LorentzVector current4Mom = output->operator[](particle)->Get4Momentum();
00431 current4Mom = toLab*current4Mom;
00432 current4Mom.boost(targ4Mom.boostVector());
00433 theSec->SetTotalEnergy(current4Mom.t());
00434 theSec->SetMomentum(current4Mom.vect());
00435 theResult->push_back(theSec);
00436
00437 #ifdef CHIPSdebug
00438 G4cout <<"CHIPS particles "<<theDefinition->GetPDGCharge()<<" "
00439 << theDefinition->GetPDGEncoding()<<" "
00440 << current4Mom <<G4endl;
00441 #endif
00442
00443 delete output->operator[](particle);
00444 }
00445 delete output;
00446 #ifdef CHIPSdebug
00447 G4cout << "Number of particles"<<theResult->size()<<G4endl;
00448 G4cout << G4endl;
00449
00450 G4cout << "QUASMON preparation info "
00451 << 1./MeV*proj4Mom<<" "
00452 << 1./MeV*targ4Mom<<" "
00453 << nD<<" "<<nU<<" "<<nS<<" "<<nAD<<" "<<nAU<<" "<<nAS<<" "
00454 << hitCount<<" "
00455 << particleCount<<" "
00456 << theLow<<" "
00457 << theHigh<<" "
00458 << G4endl;
00459 #endif
00460
00461 return theResult;
00462 }