101{
105#ifdef debug
106 G4cout<<
"G4QStringChipsParticleLevelInterface::Propagate is called"<<
G4endl;
107#endif
108
109
110 if(theSecondaries->size() == 1)
111 {
118 theFastResult->push_back(theFastSec);
119 return theFastResult;
120
121
122 }
123
124
125
126
133 unsigned int hitCount = 0;
135 {
137 {
138 resA++;
140 }
141 else
142 {
145 hitCount ++;
146 }
147 }
148 G4int targetPDGCode = 90000000 + 1000*resZ + (resA-resZ);
150 if (!resZ)
151 {
152 if (resA>1) targetMass*=resA;
153 }
156 G4double targetEnergy = std::sqrt(hitMomentum.mag2()+targetMass*targetMass);
157
159
160
163 G4double impactY = theImpact.second;
164 G4double impactPar2 = impactX*impactX + impactY*impactY;
165
167 radius2 *= radius2;
169#ifdef pdebug
170 G4cout<<
"G4QStringChipsParticleLevelInterface::Propagate: r="<<std::sqrt(radius2)/fermi
171 <<
", b="<<std::sqrt(impactPar2)/fermi<<
", R="<<theNucleus->
GetOuterRadius()/fermi
172 <<
", b/r="<<std::sqrt(impactPar2/radius2)<<
G4endl;
173#endif
174 if(radius2 - impactPar2>0) pathlength = 2.*std::sqrt(radius2 - impactPar2);
175 G4double theEnergyLostInFragmentation = theEnergyLossPerFermi*pathlength/fermi;
176
177
178 std::list<std::pair<G4double, G4KineticTrack *> > theSorted;
179 std::list<std::pair<G4double, G4KineticTrack *> >::iterator current;
180 for(unsigned int secondary = 0; secondary<theSecondaries->size(); secondary++)
181 {
182 G4LorentzVector a4Mom = theSecondaries->operator[](secondary)->Get4Momentum();
183#ifdef CHIPSdebug
184 G4cout<<
"G4QStringChipsParticleLevelInterface::Propagate: ALL STRING particles "
185 << theSecondaries->operator[](secondary)->GetDefinition()->GetPDGCharge()<<" "
186 << theSecondaries->operator[](secondary)->GetDefinition()->GetPDGEncoding()<<" "
188#endif
189#ifdef pdebug
190 G4cout<<
"G4QStringChipsParticleLevelInterface::Propagate: in C="
191 <<theSecondaries->operator[](secondary)->GetDefinition()->GetPDGCharge()<<",PDG="
192 <<theSecondaries->operator[](secondary)->GetDefinition()->GetPDGEncoding()
193 <<
",4M="<<a4Mom<<
", current nS="<<theSorted.size()<<
G4endl;
194#endif
196 std::pair<G4double, G4KineticTrack *> it;
197 it.first = toSort;
198 it.second = theSecondaries->operator[](secondary);
200 for(current = theSorted.begin(); current!=theSorted.end(); current++)
201 {
202 if((*current).first > toSort)
203 {
204 theSorted.insert(current, it);
205 inserted = true;
206 break;
207 }
208 }
209 if(!inserted) theSorted.push_back(it);
210 }
211
219 std::list<std::pair<G4double,G4KineticTrack*> >::iterator firstEscape=theSorted.begin();
221 G4int particleCount = 0;
222 G4LorentzVector theLow = (*(theSorted.begin())).second->Get4Momentum();
224
225#ifdef CHIPSdebug
226 G4cout<<
"G4QStringChipsParticleLevelInterface::Propagate: CHIPS ENERGY LOST "
227 <<theEnergyLostInFragmentation<<
". Sorted rapidities event start"<<
G4endl;
228#endif
229#ifdef pdebug
230 G4cout<<
"G4QStringChipsParticleLevelInterface::Propagate: total CHIPS energy = "
231 <<theEnergyLostInFragmentation<<". Start rapidity sorting nS="<<theSorted.size()
233#endif
234
236 std::vector<G4QContent> theContents;
237 std::vector<G4LorentzVector*> theMomenta;
241#ifdef pdebug
242 G4cout<<
"G4QStringChipsParticleLevelInterface::Propagate: Absorption nS="
243 <<theSorted.size()<<
G4endl;
244#endif
245 G4bool EscapeExists =
false;
246 for(current = theSorted.begin(); current!=theSorted.end(); current++)
247 {
248#ifdef pdebug
249 G4cout<<
"G4QStringChipsParticleLevelInterface::Propagate: nq="
250 <<(*current).second->GetDefinition()->GetQuarkContent(3)<<", naq="
251 <<(*current).second->GetDefinition()->GetAntiQuarkContent(3)<<", PDG="
252 <<(*current).second->GetDefinition()->GetPDGEncoding()<<",4M="
253 <<(*current).second->Get4Momentum()<<
G4endl;
254#endif
255 firstEscape = current;
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
306 if (baryn>0 && charge>0 && strang<1) runningEnergy-=mProt;
307 else if(baryn>0 && strang<1) runningEnergy-=mNeut;
308 else if(baryn>0) runningEnergy-=mLamb;
309 else if(baryn<0) runningEnergy+=mProt;
310
311#ifdef CHIPSdebug
312 G4cout<<
"G4QStringChipsParticleLevelInterface::Propagate: sorted rapidities "
313 <<(*current).second->Get4Momentum().rapidity()<<
G4endl;
314#endif
315
316#ifdef pdebug
317 G4cout<<
"G4QStringChipsParticleLevelInterface::Propagate: E="<<runningEnergy<<
", EL="
318 <<theEnergyLostInFragmentation<<
G4endl;
319#endif
320 if(runningEnergy > theEnergyLostInFragmentation)
321 {
322 EscapeExists = true;
323 break;
324 }
325#ifdef CHIPSdebug
326 G4cout <<
"G4QStringChipsParticleLevelInterface::Propagate: ABSORBED STRING particles "
327 <<(*current).second->GetDefinition()->GetPDGCharge()<<" "
328 << (*current).second->GetDefinition()->GetPDGEncoding()<<" "
329 << (*current).second->Get4Momentum() <<
G4endl;
330#endif
331#ifdef pdebug
332 G4cout<<
"G4QStringChipsParticleLevelInterface::Propagate:C="
333 <<current->second->GetDefinition()->GetPDGCharge()<<", PDG="
334 <<current->second->GetDefinition()->GetPDGEncoding()<<", 4M="
335 <<current->second->Get4Momentum()<<
G4endl;
336#endif
337
338
339 particleCount++;
340 theHigh = (*current).second->Get4Momentum();
341 proj4Mom = (*current).second->Get4Momentum();
342 proj4Mom.
boost(-1.*targ4Mom.boostVector());
343 nD = (*current).second->GetDefinition()->GetQuarkContent(1);
344 nU = (*current).second->GetDefinition()->GetQuarkContent(2);
345 nS = (*current).second->GetDefinition()->GetQuarkContent(3);
346 nAD = (*current).second->GetDefinition()->GetAntiQuarkContent(1);
347 nAU = (*current).second->GetDefinition()->GetAntiQuarkContent(2);
348 nAS = (*current).second->GetDefinition()->GetAntiQuarkContent(3);
349 G4QContent aProjectile(nD, nU, nS, nAD, nAU, nAS);
350
351#ifdef CHIPSdebug_1
353 G4cout <<
"G4QStringChipsParticleLevelInterface::Propagate: Quark content: d="<<nD
354 <<", u="<<nU<<", s="<<nS<< "Anti-quark content: anit-d="<<nAD<<", anti-u="<<nAU
355 <<", anti-s="<<nAS<<". G4QContent is constructed"<<endl;
356#endif
357
358 theContents.push_back(aProjectile);
360
361#ifdef CHIPSdebug_1
362 G4cout<<
"G4QStringChipsParticleLevelInterface::Propagate: projectile momentum = "
365#endif
366
367 theMomenta.push_back(aVec);
368 }
369 std::vector<G4QContent> theFinalContents;
370 std::vector<G4LorentzVector*> theFinalMomenta;
371
372 for(unsigned int hp = 0; hp<theContents.size(); hp++)
373 {
375 projHV.push_back(aHadron);
376 }
377
378
379 size_t i;
380 for (i=0; i<theFinalMomenta.size(); i++) delete theFinalMomenta[i];
381 for (i=0; i<theMomenta.size(); i++) delete theMomenta[i];
382 theFinalMomenta.clear();
383 theMomenta.clear();
384
386 fractionOfPairedQuasiFreeNucleons,
387 clusteringCoefficient,
388 fusionToExchange);
390
391#ifdef CHIPSdebug
392 G4cout<<
"G4QStringChipsParticleLevelInterface::Propagate: G4QNucleus parameters "
393 <<fractionOfSingleQuasiFreeNucleons<<" "<<fractionOfPairedQuasiFreeNucleons
394 <<
" "<<clusteringCoefficient<<
G4endl;
395 G4cout<<
"G4Quasmon parameters "<<temperature<<
" "<<halfTheStrangenessOfSee<<
" "
396 <<etaToEtaPrime <<
G4endl;
397 G4cout<<
"The Target PDG code = "<<targetPDGCode<<
G4endl;
398 G4cout<<
"The projectile momentum = "<<1./MeV*proj4Mom<<
G4endl;
399 G4cout<<
"The target momentum = "<<1./MeV*targ4Mom<<
G4endl;
400#endif
401
402
404 if (particleCount!=0 && resA!=0)
405 {
406
409#ifdef pdebug
410 G4cout<<
"G4QStringChipsParticleLevelInterface::Propagate: CHIPS fragmentation, rA="
411 <<resA<<
", #AbsPt="<<particleCount<<
G4endl;
412#endif
413 try
414 {
416 }
418 {
419 G4cerr <<
"Exception thrown of G4QStringChipsParticleLevelInterface "<<
G4endl;
421 G4cerr <<
" The projectile momentum = "<<1./MeV*proj4Mom<<
G4endl;
423 G4cerr <<
" Dumping the information in the pojectile list"<<
G4endl;
424 for(i=0; i< projHV.size(); i++)
425 {
426 G4cerr <<
" Incoming 4-momentum and PDG code of "<<i<<
"'th hadron: "
427 <<
" "<< projHV[i]->Get4Momentum()<<
" "<<projHV[i]->GetPDGCode()<<
G4endl;
428 }
429 throw;
430 }
431
432 std::for_each(projHV.begin(), projHV.end(),
DeleteQHadron());
433 projHV.clear();
434 delete pan;
435 }
436 else
437 {
438#ifdef pdebug
439 G4cout<<
"G4QStringChipsParticleLevelInterface::Propagate: NO CHIPS fragmentation, rA="
440 <<resA<<
", #AbsPt="<<particleCount<<
G4endl;
441#endif
443 }
444
445#ifdef CHIPSdebug
446 G4cout <<
"NEXT EVENT, EscapeExists="<<EscapeExists<<
G4endl;
447#endif
448
449
450 if (EscapeExists) for (current = firstEscape; current!=theSorted.end(); current++)
451 {
454 secondaries = NULL;
456 {
457 secondaries = aResult->
Decay();
458 }
459 if ( secondaries == NULL )
460 {
463 current4Mom.
boost(targ4Mom.boostVector());
466#ifdef pdebug
467 G4cout<<
"G4QStringChipsParticleLevelInterface::Propagate: *OUT* QGS stable PDG="
469#endif
470 theResult->push_back(theSec);
471 }
472 else
473 {
474 for (unsigned int aSecondary=0; aSecondary<secondaries->size(); aSecondary++)
475 {
476 theSec=
new G4ReactionProduct(secondaries->operator[](aSecondary)->GetDefinition());
477 G4LorentzVector current4Mom = secondaries->operator[](aSecondary)->Get4Momentum();
478 current4Mom.
boost(targ4Mom.boostVector());
481#ifdef pdebug
482 G4cout<<
"G4QStringChipsParticleLevelInterface::Propagate: *OUT* QGS decay PDG="
483 <<secondaries->operator[](aSecondary)->GetDefinition()->GetPDGEncoding()
484 <<
",4M="<<current4Mom<<
G4endl;
485#endif
486 theResult->push_back(theSec);
487 }
489 delete secondaries;
490 }
491 }
493 delete theSecondaries;
494
495
496 G4int maxParticle=output->size();
497#ifdef CHIPSdebug
498 G4cout <<
"Number of particles from string"<<theResult->size()<<
G4endl;
499 G4cout <<
"Number of particles from chips"<<maxParticle<<
G4endl;
500#endif
501#ifdef pdebug
502 G4cout <<
"Number of particles from QGS="<<theResult->size()<<
G4endl;
503 G4cout <<
"Number of particles from CHIPS="<<maxParticle<<
G4endl;
504#endif
505 if(maxParticle)
for(
G4int particle = 0; particle < maxParticle; particle++)
506 {
507 if(output->operator[](particle)->GetNFragments() != 0)
508 {
509 delete output->operator[](particle);
510 continue;
511 }
512 G4int pdgCode = output->operator[](particle)->GetPDGCode();
513
514
515#ifdef CHIPSdebug
516 G4cerr <<
"PDG code of chips particle = "<<pdgCode<<
G4endl;
517#endif
518
520
521
522
523
524
525 if(pdgCode>90000000)
526 {
527 G4int aZ = (pdgCode-90000000)/1000;
528 if (aZ>1000) aZ=aZ%1000;
529 G4int anN = pdgCode-90000000-1000*aZ;
530 if(anN>1000) anN=anN%1000;
541 }
543
544#ifdef CHIPSdebug
545 G4cout <<
"Particle code produced = "<< pdgCode <<
G4endl;
546#endif
547
548 if(theDefinition)
549 {
551 G4LorentzVector current4Mom = output->operator[](particle)->Get4Momentum();
552 current4Mom.
boost(targ4Mom.boostVector());
555#ifdef pdebug
556 G4cout<<
"G4QStringChipsParticleLevelInterface::Propagate: *OUT* CHIPS PDG="
558#endif
559 theResult->push_back(theSec);
560 }
561 else
562 {
564 G4cerr <<
"Getting unknown pdgCode from chips in ParticleLevelInterface"<<
G4endl;
566 }
567
568#ifdef CHIPSdebug
571 << output->operator[](particle)->Get4Momentum() <<
G4endl;
572#endif
573
574 delete output->operator[](particle);
575 }
576 else
577 {
578 if(resA>0)
579 {
581 if(resA==1)
582 {
584 }
589 theResult->push_back(theSec);
590 if(!resZ && resA>0)
for(
G4int ni=1; ni<resA; ni++)
591 {
595 theResult->push_back(theSec);
596 }
597 }
598 }
599 delete output;
600
601#ifdef CHIPSdebug
602 G4cout <<
"Number of particles"<<theResult->size()<<
G4endl;
604 G4cout <<
"QUASMON preparation info "
605 << 1./MeV*proj4Mom<<" "
606 << 1./MeV*targ4Mom<<" "
607 << nD<<" "<<nU<<" "<<nS<<" "<<nAD<<" "<<nAU<<" "<<nAS<<" "
608 << hitCount<<" "
609 << particleCount<<" "
610 << theLow<<" "
611 << theHigh<<" "
613#endif
614
615 return theResult;
616}
CLHEP::HepLorentzVector G4LorentzVector
std::vector< G4QHadron * > G4QHadronVector
std::vector< G4ReactionProduct * > G4ReactionProductVector
CLHEP::Hep3Vector G4ThreeVector
G4DLLIMPORT std::ostream G4cerr
HepLorentzVector & boost(double, double, double)
G4KineticTrackVector * Decay()
G4ParticleDefinition * GetDefinition() const
const G4LorentzVector & Get4Momentum() const
static G4Lambda * Lambda()
static G4Lambda * LambdaDefinition()
static G4Neutron * Neutron()
virtual G4ParticleDefinition * GetDefinition() const
const G4LorentzVector & GetMomentum() const
G4double GetPDGMass() const
G4int GetPDGEncoding() const
G4int GetQuarkContent(G4int flavor) const
G4double GetPDGWidth() const
G4double GetPDGCharge() const
G4int GetBaryonNumber() const
G4double GetPDGLifeTime() const
G4ParticleDefinition * FindIon(G4int atomicNumber, G4int atomicMass, G4double excitationEnergy)
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
static G4ParticleTable * GetParticleTable()
static G4Proton * Proton()
static G4QCHIPSWorld * Get()
G4QParticleVector * GetParticles(G4int nOfParts=0)
G4QHadronVector * Fragment()
static void SetParameters(G4double fN=.1, G4double fD=.05, G4double cP=4., G4double mR=1., G4double nD=.8 *CLHEP::fermi)
static void SetParameters(G4double temper=180., G4double ssin2g=.3, G4double etaetap=.3)
void SetMomentum(const G4double x, const G4double y, const G4double z)
void SetTotalEnergy(const G4double en)
std::pair< G4double, G4double > RefetchImpactXandY()
virtual G4double GetOuterRadius()=0
virtual G4Nucleon * GetNextNucleon()=0
virtual G4bool StartLoop()=0
virtual G4double GetNuclearRadius()=0