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 #include <vector>
00035
00036 #include "G4QStringChipsParticleLevelInterface.hh"
00037 #include "globals.hh"
00038 #include "G4SystemOfUnits.hh"
00039 #include "G4KineticTrackVector.hh"
00040 #include "G4Nucleon.hh"
00041 #include "G4Proton.hh"
00042 #include "G4Neutron.hh"
00043 #include "G4LorentzRotation.hh"
00044 #include "G4HadronicException.hh"
00045
00046
00047
00048 G4QStringChipsParticleLevelInterface::G4QStringChipsParticleLevelInterface()
00049 {
00050 #ifdef debug
00051 G4cout<<"G4QStringChipsParticleLevelInterface::Constructor is called"<<G4endl;
00052 #endif
00053 theEnergyLossPerFermi = 1.*GeV;
00054 nop = 152;
00055 fractionOfSingleQuasiFreeNucleons = 0.5;
00056 fractionOfPairedQuasiFreeNucleons = 0.05;
00057 clusteringCoefficient = 5.;
00058 temperature = 180.;
00059 halfTheStrangenessOfSee = 0.3;
00060 etaToEtaPrime = 0.3;
00061 fusionToExchange = 1.;
00062 theInnerCoreDensityCut = 50.;
00063
00064 if(getenv("ChipsParameterTuning"))
00065 {
00066 G4cout << "Please enter the energy loss per fermi in GeV"<<G4endl;
00067 G4cin >> theEnergyLossPerFermi;
00068 theEnergyLossPerFermi *= GeV;
00069 G4cout << "Please enter nop"<<G4endl;
00070 G4cin >> nop;
00071 G4cout << "Please enter the fractionOfSingleQuasiFreeNucleons"<<G4endl;
00072 G4cin >> fractionOfSingleQuasiFreeNucleons;
00073 G4cout << "Please enter the fractionOfPairedQuasiFreeNucleons"<<G4endl;
00074 G4cin >> fractionOfPairedQuasiFreeNucleons;
00075 G4cout << "Please enter the clusteringCoefficient"<<G4endl;
00076 G4cin >> clusteringCoefficient;
00077 G4cout << "Please enter the temperature"<<G4endl;
00078 G4cin >> temperature;
00079 G4cout << "Please enter halfTheStrangenessOfSee"<<G4endl;
00080 G4cin >> halfTheStrangenessOfSee;
00081 G4cout << "Please enter the etaToEtaPrime"<<G4endl;
00082 G4cin >> etaToEtaPrime;
00083 G4cout << "Please enter the fusionToExchange"<<G4endl;
00084 G4cin >> fusionToExchange;
00085 G4cout << "Please enter cut-off for calculating the nuclear radius in percent"<<G4endl;
00086 G4cin >> theInnerCoreDensityCut;
00087 }
00088 }
00089
00090 G4HadFinalState* G4QStringChipsParticleLevelInterface::
00091 ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& theNucleus)
00092 {
00093 #ifdef debug
00094 G4cout<<"G4QStringChipsParticleLevelInterface::ApplyYourself is called"<<G4endl;
00095 #endif
00096 return theModel.ApplyYourself(aTrack, theNucleus);
00097 }
00098
00099 G4ReactionProductVector* G4QStringChipsParticleLevelInterface::
00100 Propagate(G4KineticTrackVector* theSecondaries, G4V3DNucleus* theNucleus)
00101 {
00102 static const G4double mProt=G4Proton::Proton()->GetPDGMass();
00103 static const G4double mNeut=G4Neutron::Neutron()->GetPDGMass();
00104 static const G4double mLamb=G4Lambda::Lambda()->GetPDGMass();
00105 #ifdef debug
00106 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate is called"<<G4endl;
00107 #endif
00108
00109
00110 if(theSecondaries->size() == 1)
00111 {
00112 G4ReactionProductVector* theFastResult = new G4ReactionProductVector;
00113 G4ReactionProduct* theFastSec;
00114 theFastSec = new G4ReactionProduct((*theSecondaries)[0]->GetDefinition());
00115 G4LorentzVector current4Mom = (*theSecondaries)[0]->Get4Momentum();
00116 theFastSec->SetTotalEnergy(current4Mom.t());
00117 theFastSec->SetMomentum(current4Mom.vect());
00118 theFastResult->push_back(theFastSec);
00119 return theFastResult;
00120
00121
00122 }
00123
00124
00125
00126
00127 theNucleus->StartLoop();
00128 G4Nucleon * aNucleon;
00129 G4int resA = 0;
00130 G4int resZ = 0;
00131 G4ThreeVector hitMomentum(0,0,0);
00132 G4double hitMass = 0;
00133 unsigned int hitCount = 0;
00134 while((aNucleon = theNucleus->GetNextNucleon()))
00135 {
00136 if(!aNucleon->AreYouHit())
00137 {
00138 resA++;
00139 resZ+=G4int (aNucleon->GetDefinition()->GetPDGCharge());
00140 }
00141 else
00142 {
00143 hitMomentum += aNucleon->GetMomentum().vect();
00144 hitMass += aNucleon->GetMomentum().m();
00145 hitCount ++;
00146 }
00147 }
00148 G4int targetPDGCode = 90000000 + 1000*resZ + (resA-resZ);
00149 G4double targetMass=mNeut;
00150 if (!resZ)
00151 {
00152 if (resA>1) targetMass*=resA;
00153 }
00154 else targetMass=G4ParticleTable::GetParticleTable()->FindIon(resZ,resA,0,resZ)
00155 ->GetPDGMass();
00156 G4double targetEnergy = std::sqrt(hitMomentum.mag2()+targetMass*targetMass);
00157
00158 G4LorentzVector targ4Mom(-1.*hitMomentum, targetEnergy);
00159
00160
00161 std::pair<G4double, G4double> theImpact = theNucleus->RefetchImpactXandY();
00162 G4double impactX = theImpact.first;
00163 G4double impactY = theImpact.second;
00164 G4double impactPar2 = impactX*impactX + impactY*impactY;
00165
00166 G4double radius2 = theNucleus->GetNuclearRadius(theInnerCoreDensityCut*perCent);
00167 radius2 *= radius2;
00168 G4double pathlength = 0.;
00169 #ifdef pdebug
00170 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: r="<<std::sqrt(radius2)/fermi
00171 <<", b="<<std::sqrt(impactPar2)/fermi<<", R="<<theNucleus->GetOuterRadius()/fermi
00172 <<", b/r="<<std::sqrt(impactPar2/radius2)<<G4endl;
00173 #endif
00174 if(radius2 - impactPar2>0) pathlength = 2.*std::sqrt(radius2 - impactPar2);
00175 G4double theEnergyLostInFragmentation = theEnergyLossPerFermi*pathlength/fermi;
00176
00177
00178 std::list<std::pair<G4double, G4KineticTrack *> > theSorted;
00179 std::list<std::pair<G4double, G4KineticTrack *> >::iterator current;
00180 for(unsigned int secondary = 0; secondary<theSecondaries->size(); secondary++)
00181 {
00182 G4LorentzVector a4Mom = theSecondaries->operator[](secondary)->Get4Momentum();
00183 #ifdef CHIPSdebug
00184 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: ALL STRING particles "
00185 << theSecondaries->operator[](secondary)->GetDefinition()->GetPDGCharge()<<" "
00186 << theSecondaries->operator[](secondary)->GetDefinition()->GetPDGEncoding()<<" "
00187 << a4Mom <<G4endl;
00188 #endif
00189 #ifdef pdebug
00190 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: in C="
00191 <<theSecondaries->operator[](secondary)->GetDefinition()->GetPDGCharge()<<",PDG="
00192 <<theSecondaries->operator[](secondary)->GetDefinition()->GetPDGEncoding()
00193 <<",4M="<<a4Mom<<", current nS="<<theSorted.size()<<G4endl;
00194 #endif
00195 G4double toSort = a4Mom.rapidity();
00196 std::pair<G4double, G4KineticTrack *> it;
00197 it.first = toSort;
00198 it.second = theSecondaries->operator[](secondary);
00199 G4bool inserted = false;
00200 for(current = theSorted.begin(); current!=theSorted.end(); current++)
00201 {
00202 if((*current).first > toSort)
00203 {
00204 theSorted.insert(current, it);
00205 inserted = true;
00206 break;
00207 }
00208 }
00209 if(!inserted) theSorted.push_back(it);
00210 }
00211
00212 G4LorentzVector proj4Mom(0.,0.,0.,0.);
00213 G4int nD = 0;
00214 G4int nU = 0;
00215 G4int nS = 0;
00216 G4int nAD = 0;
00217 G4int nAU = 0;
00218 G4int nAS = 0;
00219 std::list<std::pair<G4double,G4KineticTrack*> >::iterator firstEscape=theSorted.begin();
00220 G4double runningEnergy = 0;
00221 G4int particleCount = 0;
00222 G4LorentzVector theLow = (*(theSorted.begin())).second->Get4Momentum();
00223 G4LorentzVector theHigh;
00224
00225 #ifdef CHIPSdebug
00226 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: CHIPS ENERGY LOST "
00227 <<theEnergyLostInFragmentation<<". Sorted rapidities event start"<<G4endl;
00228 #endif
00229 #ifdef pdebug
00230 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: total CHIPS energy = "
00231 <<theEnergyLostInFragmentation<<". Start rapidity sorting nS="<<theSorted.size()
00232 <<G4endl;
00233 #endif
00234
00235 G4QHadronVector projHV;
00236 std::vector<G4QContent> theContents;
00237 std::vector<G4LorentzVector*> theMomenta;
00238 G4ReactionProductVector* theResult = new G4ReactionProductVector;
00239 G4ReactionProduct* theSec;
00240 G4KineticTrackVector* secondaries;
00241 #ifdef pdebug
00242 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: Absorption nS="
00243 <<theSorted.size()<<G4endl;
00244 #endif
00245 G4bool EscapeExists = false;
00246 for(current = theSorted.begin(); current!=theSorted.end(); current++)
00247 {
00248 #ifdef pdebug
00249 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: nq="
00250 <<(*current).second->GetDefinition()->GetQuarkContent(3)<<", naq="
00251 <<(*current).second->GetDefinition()->GetAntiQuarkContent(3)<<", PDG="
00252 <<(*current).second->GetDefinition()->GetPDGEncoding()<<",4M="
00253 <<(*current).second->Get4Momentum()<<G4endl;
00254 #endif
00255 firstEscape = current;
00256 G4KineticTrack* aResult = (*current).second;
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302 runningEnergy += aResult->Get4Momentum().t();
00303 G4double charge=aResult->GetDefinition()->GetPDGCharge();
00304 G4int strang=aResult->GetDefinition()->GetQuarkContent(3);
00305 G4int baryn=aResult->GetDefinition()->GetBaryonNumber();
00306 if (baryn>0 && charge>0 && strang<1) runningEnergy-=mProt;
00307 else if(baryn>0 && strang<1) runningEnergy-=mNeut;
00308 else if(baryn>0) runningEnergy-=mLamb;
00309 else if(baryn<0) runningEnergy+=mProt;
00310
00311 #ifdef CHIPSdebug
00312 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: sorted rapidities "
00313 <<(*current).second->Get4Momentum().rapidity()<<G4endl;
00314 #endif
00315
00316 #ifdef pdebug
00317 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: E="<<runningEnergy<<", EL="
00318 <<theEnergyLostInFragmentation<<G4endl;
00319 #endif
00320 if(runningEnergy > theEnergyLostInFragmentation)
00321 {
00322 EscapeExists = true;
00323 break;
00324 }
00325 #ifdef CHIPSdebug
00326 G4cout <<"G4QStringChipsParticleLevelInterface::Propagate: ABSORBED STRING particles "
00327 <<(*current).second->GetDefinition()->GetPDGCharge()<<" "
00328 << (*current).second->GetDefinition()->GetPDGEncoding()<<" "
00329 << (*current).second->Get4Momentum() <<G4endl;
00330 #endif
00331 #ifdef pdebug
00332 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate:C="
00333 <<current->second->GetDefinition()->GetPDGCharge()<<", PDG="
00334 <<current->second->GetDefinition()->GetPDGEncoding()<<", 4M="
00335 <<current->second->Get4Momentum()<<G4endl;
00336 #endif
00337
00338
00339 particleCount++;
00340 theHigh = (*current).second->Get4Momentum();
00341 proj4Mom = (*current).second->Get4Momentum();
00342 proj4Mom.boost(-1.*targ4Mom.boostVector());
00343 nD = (*current).second->GetDefinition()->GetQuarkContent(1);
00344 nU = (*current).second->GetDefinition()->GetQuarkContent(2);
00345 nS = (*current).second->GetDefinition()->GetQuarkContent(3);
00346 nAD = (*current).second->GetDefinition()->GetAntiQuarkContent(1);
00347 nAU = (*current).second->GetDefinition()->GetAntiQuarkContent(2);
00348 nAS = (*current).second->GetDefinition()->GetAntiQuarkContent(3);
00349 G4QContent aProjectile(nD, nU, nS, nAD, nAU, nAS);
00350
00351 #ifdef CHIPSdebug_1
00352 G4cout <<G4endl;
00353 G4cout <<"G4QStringChipsParticleLevelInterface::Propagate: Quark content: d="<<nD
00354 <<", u="<<nU<<", s="<<nS<< "Anti-quark content: anit-d="<<nAD<<", anti-u="<<nAU
00355 <<", anti-s="<<nAS<<". G4QContent is constructed"<<endl;
00356 #endif
00357
00358 theContents.push_back(aProjectile);
00359 G4LorentzVector* aVec = new G4LorentzVector((1./MeV)*proj4Mom);
00360
00361 #ifdef CHIPSdebug_1
00362 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: projectile momentum = "
00363 <<*aVec<<G4endl;
00364 G4cout << G4endl;
00365 #endif
00366
00367 theMomenta.push_back(aVec);
00368 }
00369 std::vector<G4QContent> theFinalContents;
00370 std::vector<G4LorentzVector*> theFinalMomenta;
00371
00372 for(unsigned int hp = 0; hp<theContents.size(); hp++)
00373 {
00374 G4QHadron* aHadron = new G4QHadron(theContents[hp], *(theMomenta[hp]) );
00375 projHV.push_back(aHadron);
00376 }
00377
00378
00379 size_t i;
00380 for (i=0; i<theFinalMomenta.size(); i++) delete theFinalMomenta[i];
00381 for (i=0; i<theMomenta.size(); i++) delete theMomenta[i];
00382 theFinalMomenta.clear();
00383 theMomenta.clear();
00384
00385 G4QNucleus::SetParameters(fractionOfSingleQuasiFreeNucleons,
00386 fractionOfPairedQuasiFreeNucleons,
00387 clusteringCoefficient,
00388 fusionToExchange);
00389 G4Quasmon::SetParameters(temperature, halfTheStrangenessOfSee, etaToEtaPrime);
00390
00391 #ifdef CHIPSdebug
00392 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: G4QNucleus parameters "
00393 <<fractionOfSingleQuasiFreeNucleons<<" "<<fractionOfPairedQuasiFreeNucleons
00394 <<" "<<clusteringCoefficient<<G4endl;
00395 G4cout<<"G4Quasmon parameters "<<temperature<<" "<<halfTheStrangenessOfSee<<" "
00396 <<etaToEtaPrime << G4endl;
00397 G4cout<<"The Target PDG code = "<<targetPDGCode<<G4endl;
00398 G4cout<<"The projectile momentum = "<<1./MeV*proj4Mom<<G4endl;
00399 G4cout<<"The target momentum = "<<1./MeV*targ4Mom<<G4endl;
00400 #endif
00401
00402
00403 G4QHadronVector* output = 0;
00404 if (particleCount!=0 && resA!=0)
00405 {
00406
00407 G4QCHIPSWorld::Get()->GetParticles(nop);
00408 G4QEnvironment* pan= new G4QEnvironment(projHV, targetPDGCode);
00409 #ifdef pdebug
00410 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: CHIPS fragmentation, rA="
00411 <<resA<<", #AbsPt="<<particleCount<<G4endl;
00412 #endif
00413 try
00414 {
00415 output = pan->Fragment();
00416 }
00417 catch(G4HadronicException& aR)
00418 {
00419 G4cerr << "Exception thrown of G4QStringChipsParticleLevelInterface "<<G4endl;
00420 G4cerr << " targetPDGCode = "<< targetPDGCode <<G4endl;
00421 G4cerr << " The projectile momentum = "<<1./MeV*proj4Mom<<G4endl;
00422 G4cerr << " The target momentum = "<<1./MeV*targ4Mom<<G4endl<<G4endl;
00423 G4cerr << " Dumping the information in the pojectile list"<<G4endl;
00424 for(i=0; i< projHV.size(); i++)
00425 {
00426 G4cerr <<" Incoming 4-momentum and PDG code of "<<i<<"'th hadron: "
00427 <<" "<< projHV[i]->Get4Momentum()<<" "<<projHV[i]->GetPDGCode()<<G4endl;
00428 }
00429 throw;
00430 }
00431
00432 std::for_each(projHV.begin(), projHV.end(), DeleteQHadron());
00433 projHV.clear();
00434 delete pan;
00435 }
00436 else
00437 {
00438 #ifdef pdebug
00439 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: NO CHIPS fragmentation, rA="
00440 <<resA<<", #AbsPt="<<particleCount<<G4endl;
00441 #endif
00442 output = new G4QHadronVector;
00443 }
00444
00445 #ifdef CHIPSdebug
00446 G4cout << "NEXT EVENT, EscapeExists="<<EscapeExists<<G4endl;
00447 #endif
00448
00449
00450 if (EscapeExists) for (current = firstEscape; current!=theSorted.end(); current++)
00451 {
00452 G4KineticTrack* aResult = (*current).second;
00453 G4ParticleDefinition* pdef=aResult->GetDefinition();
00454 secondaries = NULL;
00455 if(pdef->GetPDGWidth() > 0 && pdef->GetPDGLifeTime() < 5E-17*s )
00456 {
00457 secondaries = aResult->Decay();
00458 }
00459 if ( secondaries == NULL )
00460 {
00461 theSec = new G4ReactionProduct(aResult->GetDefinition());
00462 G4LorentzVector current4Mom = aResult->Get4Momentum();
00463 current4Mom.boost(targ4Mom.boostVector());
00464 theSec->SetTotalEnergy(current4Mom.t());
00465 theSec->SetMomentum(current4Mom.vect());
00466 #ifdef pdebug
00467 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: *OUT* QGS stable PDG="
00468 <<aResult->GetDefinition()->GetPDGEncoding()<<",4M="<<current4Mom<<G4endl;
00469 #endif
00470 theResult->push_back(theSec);
00471 }
00472 else
00473 {
00474 for (unsigned int aSecondary=0; aSecondary<secondaries->size(); aSecondary++)
00475 {
00476 theSec=new G4ReactionProduct(secondaries->operator[](aSecondary)->GetDefinition());
00477 G4LorentzVector current4Mom = secondaries->operator[](aSecondary)->Get4Momentum();
00478 current4Mom.boost(targ4Mom.boostVector());
00479 theSec->SetTotalEnergy(current4Mom.t());
00480 theSec->SetMomentum(current4Mom.vect());
00481 #ifdef pdebug
00482 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: *OUT* QGS decay PDG="
00483 <<secondaries->operator[](aSecondary)->GetDefinition()->GetPDGEncoding()
00484 <<",4M="<<current4Mom<<G4endl;
00485 #endif
00486 theResult->push_back(theSec);
00487 }
00488 std::for_each(secondaries->begin(), secondaries->end(), DeleteKineticTrack());
00489 delete secondaries;
00490 }
00491 }
00492 std::for_each(theSecondaries->begin(), theSecondaries->end(), DeleteKineticTrack());
00493 delete theSecondaries;
00494
00495
00496 G4int maxParticle=output->size();
00497 #ifdef CHIPSdebug
00498 G4cout << "Number of particles from string"<<theResult->size()<<G4endl;
00499 G4cout << "Number of particles from chips"<<maxParticle<<G4endl;
00500 #endif
00501 #ifdef pdebug
00502 G4cout << "Number of particles from QGS="<<theResult->size()<<G4endl;
00503 G4cout << "Number of particles from CHIPS="<<maxParticle<<G4endl;
00504 #endif
00505 if(maxParticle) for(G4int particle = 0; particle < maxParticle; particle++)
00506 {
00507 if(output->operator[](particle)->GetNFragments() != 0)
00508 {
00509 delete output->operator[](particle);
00510 continue;
00511 }
00512 G4int pdgCode = output->operator[](particle)->GetPDGCode();
00513
00514
00515 #ifdef CHIPSdebug
00516 G4cerr << "PDG code of chips particle = "<<pdgCode<<G4endl;
00517 #endif
00518
00519 G4ParticleDefinition * theDefinition;
00520
00521
00522
00523
00524
00525 if(pdgCode>90000000)
00526 {
00527 G4int aZ = (pdgCode-90000000)/1000;
00528 if (aZ>1000) aZ=aZ%1000;
00529 G4int anN = pdgCode-90000000-1000*aZ;
00530 if(anN>1000) anN=anN%1000;
00531 if(pdgCode==91000000) theDefinition = G4Lambda::LambdaDefinition();
00532 else if(pdgCode==92000000) theDefinition = G4Lambda::LambdaDefinition();
00533 else if(pdgCode==93000000) theDefinition = G4Lambda::LambdaDefinition();
00534 else if(pdgCode==94000000) theDefinition = G4Lambda::LambdaDefinition();
00535 else if(pdgCode==95000000) theDefinition = G4Lambda::LambdaDefinition();
00536 else if(pdgCode==96000000) theDefinition = G4Lambda::LambdaDefinition();
00537 else if(pdgCode==97000000) theDefinition = G4Lambda::LambdaDefinition();
00538 else if(pdgCode==98000000) theDefinition = G4Lambda::LambdaDefinition();
00539 else if(aZ == 0 && anN == 1) theDefinition = G4Neutron::Neutron();
00540 else theDefinition = G4ParticleTable::GetParticleTable()->FindIon(aZ,anN+aZ,0,aZ);
00541 }
00542 else theDefinition = G4ParticleTable::GetParticleTable()->FindParticle(pdgCode);
00543
00544 #ifdef CHIPSdebug
00545 G4cout << "Particle code produced = "<< pdgCode <<G4endl;
00546 #endif
00547
00548 if(theDefinition)
00549 {
00550 theSec = new G4ReactionProduct(theDefinition);
00551 G4LorentzVector current4Mom = output->operator[](particle)->Get4Momentum();
00552 current4Mom.boost(targ4Mom.boostVector());
00553 theSec->SetTotalEnergy(current4Mom.t());
00554 theSec->SetMomentum(current4Mom.vect());
00555 #ifdef pdebug
00556 G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: *OUT* CHIPS PDG="
00557 <<theDefinition->GetPDGEncoding()<<",4M="<<current4Mom<<G4endl;
00558 #endif
00559 theResult->push_back(theSec);
00560 }
00561 else
00562 {
00563 G4cerr << G4endl<<"WARNING: "<<G4endl;
00564 G4cerr << "Getting unknown pdgCode from chips in ParticleLevelInterface"<<G4endl;
00565 G4cerr << "skipping particle with pdgCode = "<<pdgCode<<G4endl<<G4endl;
00566 }
00567
00568 #ifdef CHIPSdebug
00569 G4cout <<"CHIPS particles "<<theDefinition->GetPDGCharge()<<" "
00570 << theDefinition->GetPDGEncoding()<<" "
00571 << output->operator[](particle)->Get4Momentum() <<G4endl;
00572 #endif
00573
00574 delete output->operator[](particle);
00575 }
00576 else
00577 {
00578 if(resA>0)
00579 {
00580 G4ParticleDefinition* theDefinition = G4Neutron::Neutron();
00581 if(resA==1)
00582 {
00583 if(resZ == 1) theDefinition = G4Proton::Proton();
00584 }
00585 else theDefinition = G4ParticleTable::GetParticleTable()->FindIon(resZ,resA,0,resZ);
00586 theSec = new G4ReactionProduct(theDefinition);
00587 theSec->SetTotalEnergy(theDefinition->GetPDGMass());
00588 theSec->SetMomentum(G4ThreeVector(0.,0.,0.));
00589 theResult->push_back(theSec);
00590 if(!resZ && resA>0) for(G4int ni=1; ni<resA; ni++)
00591 {
00592 theSec = new G4ReactionProduct(theDefinition);
00593 theSec->SetTotalEnergy(theDefinition->GetPDGMass());
00594 theSec->SetMomentum(G4ThreeVector(0.,0.,0.));
00595 theResult->push_back(theSec);
00596 }
00597 }
00598 }
00599 delete output;
00600
00601 #ifdef CHIPSdebug
00602 G4cout << "Number of particles"<<theResult->size()<<G4endl;
00603 G4cout << G4endl;
00604 G4cout << "QUASMON preparation info "
00605 << 1./MeV*proj4Mom<<" "
00606 << 1./MeV*targ4Mom<<" "
00607 << nD<<" "<<nU<<" "<<nS<<" "<<nAD<<" "<<nAU<<" "<<nAS<<" "
00608 << hitCount<<" "
00609 << particleCount<<" "
00610 << theLow<<" "
00611 << theHigh<<" "
00612 << G4endl;
00613 #endif
00614
00615 return theResult;
00616 }