source: trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4FTFModel.cc@ 1315

Last change on this file since 1315 was 1228, checked in by garnier, 16 years ago

update geant4.9.3 tag
File size: 27.7 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id: G4FTFModel.cc,v 1.34 2009/12/15 19:14:31 vuzhinsk Exp $
28// GEANT4 tag $Name: geant4-09-03 $
29//
30
31// ------------------------------------------------------------
32// GEANT 4 class implementation file
33//
34// ---------------- G4FTFModel ----------------
35// by Gunter Folger, May 1998.
36// class implementing the excitation in the FTF Parton String Model
37// ------------------------------------------------------------
38
39#include "G4FTFModel.hh"
40#include "G4FTFParameters.hh"
41#include "G4FTFParticipants.hh"
42#include "G4DiffractiveSplitableHadron.hh"
43#include "G4InteractionContent.hh"
44#include "G4LorentzRotation.hh"
45#include "G4ParticleDefinition.hh"
46#include "G4ParticleTable.hh"
47#include "G4ios.hh"
48#include <utility>
49#include "G4IonTable.hh"
50
51// Class G4FTFModel
52
53G4FTFModel::G4FTFModel():theExcitation(new G4DiffractiveExcitation()),
54 theElastic(new G4ElasticHNScattering())
55{
56 G4VPartonStringModel::SetThisPointer(this);
57 theParameters=0;
58 NumberOfInvolvedNucleon=0;
59}
60
61
62G4FTFModel::~G4FTFModel()
63{
64// Because FTF model can be called for various particles
65// theParameters must be erased at the end of each call.
66// Thus the delete is also in G4FTFModel::GetStrings() method
67 if( theParameters != 0 ) delete theParameters;
68 if( theExcitation != 0 ) delete theExcitation;
69 if( theElastic != 0 ) delete theElastic;
70
71 if( NumberOfInvolvedNucleon != 0)
72 {
73 for(G4int i=0; i < NumberOfInvolvedNucleon; i++)
74 {
75 G4VSplitableHadron * aNucleon = TheInvolvedNucleon[i]->GetSplitableHadron();
76 if(aNucleon) delete aNucleon;
77 }
78 }
79}
80
81const G4FTFModel & G4FTFModel::operator=(const G4FTFModel &)
82{
83 throw G4HadronicException(__FILE__, __LINE__, "G4FTFModel::operator= is not meant to be accessed ");
84 return *this;
85}
86
87int G4FTFModel::operator==(const G4FTFModel &right) const
88{
89 return this==&right;
90}
91
92int G4FTFModel::operator!=(const G4FTFModel &right) const
93{
94 return this!=&right;
95}
96
97// ------------------------------------------------------------
98void G4FTFModel::Init(const G4Nucleus & aNucleus, const G4DynamicParticle & aProjectile)
99{
100 theProjectile = aProjectile;
101 theParticipants.Init(aNucleus.GetN(),aNucleus.GetZ());
102// ----------- N-mass number Z-charge -------------------------
103
104// --- cms energy
105 G4double s = sqr( theProjectile.GetMass() ) +
106 sqr( G4Proton::Proton()->GetPDGMass() ) +
107 2*theProjectile.GetTotalEnergy()*G4Proton::Proton()->GetPDGMass();
108
109 if( theParameters != 0 ) delete theParameters;
110 theParameters = new G4FTFParameters(theProjectile.GetDefinition(),
111 aNucleus.GetN(),aNucleus.GetZ(),
112 s);
113
114//theParameters->SetProbabilityOfElasticScatt(0.);
115// To turn on/off (1/0) elastic scattering
116
117}
118
119// ------------------------------------------------------------
120struct DeleteVSplitableHadron { void operator()(G4VSplitableHadron * aH){ delete aH;} };
121
122
123// ------------------------------------------------------------
124G4ExcitedStringVector * G4FTFModel::GetStrings()
125{
126 G4ExcitedStringVector * theStrings(0);
127
128 theParticipants.GetList(theProjectile,theParameters);
129
130 ReggeonCascade();
131
132 G4bool Success(true);
133 if( PutOnMassShell() )
134 {
135 if( ExciteParticipants() )
136 {
137 theStrings = BuildStrings();
138
139 GetResidualNucleus();
140
141 if( theParameters != 0 )
142 {
143 delete theParameters;
144 theParameters=0;
145 }
146 } else // if( ExciteParticipants() )
147 { Success=false;}
148 } else // if( PutOnMassShell() )
149 { Success=false;}
150
151 if(!Success)
152 {
153 // -------------- Erase the projectile ----------------
154 std::vector<G4VSplitableHadron *> primaries;
155
156 theParticipants.StartLoop(); // restart a loop
157 while ( theParticipants.Next() )
158 {
159 const G4InteractionContent & interaction=theParticipants.GetInteraction();
160 // do not allow for duplicates ...
161 if ( primaries.end() == std::find(primaries.begin(), primaries.end(),
162 interaction.GetProjectile()) )
163 primaries.push_back(interaction.GetProjectile());
164 }
165 std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron());
166 primaries.clear();
167 }
168
169// -------------- Cleaning of the memory --------------
170// -------------- Erase the target nucleons -----------
171 G4VSplitableHadron * aNucleon = 0;
172 for(G4int i=0; i < NumberOfInvolvedNucleon; i++)
173 {
174 aNucleon = TheInvolvedNucleon[i]->GetSplitableHadron();
175 if(aNucleon) delete aNucleon;
176 }
177
178 NumberOfInvolvedNucleon=0;
179
180 return theStrings;
181
182}
183//-------------------------------------------------------------------
184void G4FTFModel::ReggeonCascade()
185{ //--- Implementation of reggeon theory inspired model-------
186 NumberOfInvolvedNucleon=0;
187
188 theParticipants.StartLoop();
189 while (theParticipants.Next())
190 {
191 const G4InteractionContent & collision=theParticipants.GetInteraction();
192 G4Nucleon * TargetNucleon=collision.GetTargetNucleon();
193
194 TheInvolvedNucleon[NumberOfInvolvedNucleon]=TargetNucleon;
195 NumberOfInvolvedNucleon++;
196
197 G4double XofWoundedNucleon = TargetNucleon->GetPosition().x();
198 G4double YofWoundedNucleon = TargetNucleon->GetPosition().y();
199
200 theParticipants.theNucleus->StartLoop();
201 G4Nucleon * Neighbour(0);
202
203 while ( (Neighbour = theParticipants.theNucleus->GetNextNucleon()) )
204 {
205 if(!Neighbour->AreYouHit())
206 {
207 G4double impact2= sqr(XofWoundedNucleon - Neighbour->GetPosition().x()) +
208 sqr(YofWoundedNucleon - Neighbour->GetPosition().y());
209
210 if(G4UniformRand() < theParameters->GetCofNuclearDestruction()*
211 std::exp(-impact2/theParameters->GetR2ofNuclearDestruction()))
212 { // The neighbour nucleon is involved in the reggeon cascade
213
214 TheInvolvedNucleon[NumberOfInvolvedNucleon]=Neighbour;
215 NumberOfInvolvedNucleon++;
216
217 G4VSplitableHadron *targetSplitable;
218 targetSplitable = new G4DiffractiveSplitableHadron(*Neighbour);
219
220 Neighbour->Hit(targetSplitable);
221 targetSplitable->SetStatus(2);
222 }
223 } // end of if(!Neighbour->AreYouHit())
224 } // end of while (theParticipant.theNucleus->GetNextNucleon())
225 } // end of while (theParticipants.Next())
226
227// ---------------- Calculation of creation time for each target nucleon -----------
228 theParticipants.StartLoop(); // restart a loop
229 theParticipants.Next();
230 G4VSplitableHadron * primary = theParticipants.GetInteraction().GetProjectile();
231 G4double betta_z=primary->Get4Momentum().pz()/primary->Get4Momentum().e();
232 primary->SetTimeOfCreation(0.);
233
234 G4double ZcoordinateOfPreviousCollision(0.);
235 G4double ZcoordinateOfCurrentInteraction(0.);
236 G4double TimeOfPreviousCollision(0.);
237 G4double TimeOfCurrentCollision(0);
238
239 theParticipants.theNucleus->StartLoop();
240 G4Nucleon * aNucleon;
241 G4bool theFirstInvolvedNucleon(true);
242 while ( (aNucleon = theParticipants.theNucleus->GetNextNucleon()) )
243 {
244 if(aNucleon->AreYouHit())
245 {
246 if(theFirstInvolvedNucleon)
247 {
248 ZcoordinateOfPreviousCollision=aNucleon->GetPosition().z();
249 theFirstInvolvedNucleon=false;
250 }
251
252 ZcoordinateOfCurrentInteraction=aNucleon->GetPosition().z();
253 TimeOfCurrentCollision=TimeOfPreviousCollision+
254 (ZcoordinateOfCurrentInteraction-ZcoordinateOfPreviousCollision)/betta_z;
255// It is assumed that the nucleons are ordered on increasing z-coordinate ------------
256 aNucleon->GetSplitableHadron()->SetTimeOfCreation(TimeOfCurrentCollision);
257
258 ZcoordinateOfPreviousCollision=ZcoordinateOfCurrentInteraction;
259 TimeOfPreviousCollision=TimeOfCurrentCollision;
260 } // end of if(aNucleon->AreYouHit())
261 } // end of while (theParticipant.theNucleus->GetNextNucleon())
262//
263// The algorithm can be improved, but it will be more complicated, and will require
264// changes in G4DiffractiveExcitation.cc and G4ElasticHNScattering.cc
265} // Uzhi 26 July 2009
266
267// ------------------------------------------------------------
268G4bool G4FTFModel::PutOnMassShell()
269{
270// -------------- Properties of the projectile ----------------
271 theParticipants.StartLoop(); // restart a loop
272 theParticipants.Next();
273 G4VSplitableHadron * primary = theParticipants.GetInteraction().GetProjectile();
274 G4LorentzVector Pprojectile=primary->Get4Momentum();
275
276// To get original projectile particle
277
278 if(Pprojectile.z() < 0.){return false;}
279
280 G4double Mprojectile = Pprojectile.mag();
281 G4double M2projectile = Pprojectile.mag2();
282//-------------------------------------------------------------
283 G4LorentzVector Psum = Pprojectile;
284 G4double SumMasses = Mprojectile + 20.*MeV; // 13.12.09
285 // Separation energy for projectile
286
287//--------------- Target nucleus ------------------------------
288 G4V3DNucleus *theNucleus = GetWoundedNucleus();
289 G4Nucleon * aNucleon;
290 G4int ResidualMassNumber=theNucleus->GetMassNumber();
291 G4int ResidualCharge =theNucleus->GetCharge();
292 ResidualExcitationEnergy=0.;
293 G4LorentzVector PnuclearResidual(0.,0.,0.,0.);
294
295 G4double ExcitationEnergyPerWoundedNucleon=
296 theParameters->GetExcitationEnergyPerWoundedNucleon();
297
298 theNucleus->StartLoop();
299
300 while ((aNucleon = theNucleus->GetNextNucleon()))
301 {
302 if(aNucleon->AreYouHit())
303 { // Involved nucleons
304 Psum += aNucleon->Get4Momentum();
305 SumMasses += aNucleon->GetDefinition()->GetPDGMass();
306 SumMasses += 20.*MeV; // 13.12.09 Separation energy for a nucleon
307 ResidualMassNumber--;
308 ResidualCharge-=(G4int) aNucleon->GetDefinition()->GetPDGCharge();
309 ResidualExcitationEnergy+=ExcitationEnergyPerWoundedNucleon;
310 }
311 else
312 { // Spectator nucleons
313 PnuclearResidual += aNucleon->Get4Momentum();
314 } // end of if(!aNucleon->AreYouHit())
315 } // end of while (theNucleus->GetNextNucleon())
316
317 Psum += PnuclearResidual;
318
319 G4double ResidualMass(0.);
320 if(ResidualMassNumber == 0)
321 {
322 ResidualMass=0.;
323 ResidualExcitationEnergy=0.;
324 }
325 else
326 {
327 ResidualMass=G4ParticleTable::GetParticleTable()->GetIonTable()->
328 GetIonMass(ResidualCharge ,ResidualMassNumber);
329 if(ResidualMassNumber == 1) {ResidualExcitationEnergy=0.;}
330 }
331
332// ResidualMass +=ResidualExcitationEnergy; // Will be given after checks
333 SumMasses += ResidualMass;
334
335//-------------------------------------------------------------
336
337 G4double SqrtS=Psum.mag();
338 G4double S=Psum.mag2();
339
340 if(SqrtS < SumMasses) {return false;} // It is impossible to simulate
341 // after putting nuclear nucleons
342 // on mass-shell
343
344 if(SqrtS < SumMasses+ResidualExcitationEnergy) {ResidualExcitationEnergy=0.;}
345
346 ResidualMass +=ResidualExcitationEnergy;
347 SumMasses +=ResidualExcitationEnergy;
348
349//-------------------------------------------------------------
350// Sampling of nucleons what are transfered to delta-isobars --
351 G4int MaxNumberOfDeltas = (int)((SqrtS - SumMasses)/(400.*MeV));
352 G4int NumberOfDeltas(0);
353
354 if(theNucleus->GetMassNumber() != 1)
355 {
356 G4double ProbDeltaIsobar(0.); // 1. *** Can be set if it is needed
357 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
358 {
359 if((G4UniformRand() < ProbDeltaIsobar)&&(NumberOfDeltas < MaxNumberOfDeltas))
360 {
361 NumberOfDeltas++;
362 G4VSplitableHadron * targetSplitable=TheInvolvedNucleon[i]->GetSplitableHadron();
363 SumMasses-=targetSplitable->GetDefinition()->GetPDGMass();
364
365 G4int PDGcode = targetSplitable->GetDefinition()->GetPDGEncoding();
366 G4int newPDGcode = PDGcode/10; newPDGcode=newPDGcode*10+4; // Delta
367 G4ParticleDefinition* ptr =
368 G4ParticleTable::GetParticleTable()->FindParticle(newPDGcode);
369 targetSplitable->SetDefinition(ptr);
370 SumMasses+=targetSplitable->GetDefinition()->GetPDGMass();
371 }
372 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
373 } // end of if(theNucleus.GetMassNumber() != 1)
374//-------------------------------------------------------------
375 G4LorentzRotation toCms(-1*Psum.boostVector());
376 G4LorentzVector Ptmp=toCms*Pprojectile;
377 if ( Ptmp.pz() <= 0. )
378 { // "String" moving backwards in CMS, abort collision !!
379 //G4cout << " abort ColliDeleteVSplitableHadronsion!! " << G4endl;
380 return false;
381 }
382
383// toCms.rotateZ(-1*Ptmp.phi()); // Uzhi 5.12.09
384// toCms.rotateY(-1*Ptmp.theta()); // Uzhi 5.12.09
385
386 G4LorentzRotation toLab(toCms.inverse());
387
388//-------------------------------------------------------------
389//------- Ascribing of the involved nucleons Pt and Xminus ----
390 G4double Dcor = theParameters->GetDofNuclearDestruction()/
391 theNucleus->GetMassNumber();
392
393 G4double AveragePt2 = theParameters->GetPt2ofNuclearDestruction();
394 G4double maxPtSquare = theParameters->GetMaxPt2ofNuclearDestruction();
395
396 G4double M2target(0.);
397 G4double WminusTarget(0.);
398 G4double WplusProjectile(0.);
399
400 G4int NumberOfTries(0);
401 G4double ScaleFactor(1.);
402 G4bool OuterSuccess(true);
403 do // while (!OuterSuccess)
404 {
405 OuterSuccess=true;
406
407 do // while (SqrtS < Mprojectile + std::sqrt(M2target))
408 { // while (DecayMomentum < 0.)
409
410 NumberOfTries++;
411 if(NumberOfTries == 100*(NumberOfTries/100)) // 100
412 { // At large number of tries it would be better to reduce the values
413 ScaleFactor/=2.;
414 Dcor *=ScaleFactor;
415 AveragePt2 *=ScaleFactor;
416 }
417
418 G4ThreeVector PtSum(0.,0.,0.);
419 G4double XminusSum(0.);
420 G4double Xminus(0.);
421 G4bool InerSuccess=true;
422
423 do // while(!InerSuccess);
424 {
425 InerSuccess=true;
426
427 PtSum =G4ThreeVector(0.,0.,0.);
428 XminusSum=0.;
429
430 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
431 {
432 G4Nucleon * aNucleon = TheInvolvedNucleon[i];
433
434 G4ThreeVector tmpPt = GaussianPt(AveragePt2, maxPtSquare);
435 PtSum += tmpPt;
436 G4ThreeVector tmpX=GaussianPt(Dcor*Dcor, 1.);
437 Xminus=tmpX.x();
438 XminusSum+=Xminus;
439
440 G4LorentzVector tmp(tmpPt.x(),tmpPt.y(),Xminus,0.);
441 aNucleon->SetMomentum(tmp);
442 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
443
444//---------------------------------------------------------------------------
445 G4double DeltaX(0.);
446 G4double DeltaY(0.);
447 G4double DeltaXminus(0.);
448
449 if(ResidualMassNumber == 0)
450 {
451 DeltaX = PtSum.x()/NumberOfInvolvedNucleon;
452 DeltaY = PtSum.y()/NumberOfInvolvedNucleon;
453 DeltaXminus = (XminusSum-1.)/NumberOfInvolvedNucleon;
454 }
455 else
456 {
457 DeltaXminus = -1./theNucleus->GetMassNumber();
458 }
459
460 XminusSum=1.;
461 M2target =0.;
462
463 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
464 {
465 G4Nucleon * aNucleon = TheInvolvedNucleon[i];
466
467 Xminus = aNucleon->Get4Momentum().pz() - DeltaXminus;
468 XminusSum-=Xminus;
469
470 if((Xminus <= 0.) || (Xminus > 1.) ||
471 (XminusSum <=0.) || (XminusSum > 1.)) {InerSuccess=false; break;}
472
473 G4double Px=aNucleon->Get4Momentum().px() - DeltaX;
474 G4double Py=aNucleon->Get4Momentum().py() - DeltaY;
475
476 M2target +=(aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass()*
477 aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass() +
478 Px*Px + Py*Py)/Xminus;
479
480 G4LorentzVector tmp(Px,Py,Xminus,0.);
481 aNucleon->SetMomentum(tmp);
482 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
483
484 if(InerSuccess && (ResidualMassNumber != 0))
485 {
486 M2target +=(ResidualMass*ResidualMass + PtSum.mag2())/XminusSum;
487 }
488 } while(!InerSuccess);
489 } while (SqrtS < Mprojectile + std::sqrt(M2target));
490//-------------------------------------------------------------
491 G4double DecayMomentum2= S*S+M2projectile*M2projectile+M2target*M2target
492 -2.*S*M2projectile - 2.*S*M2target
493 -2.*M2projectile*M2target;
494
495 WminusTarget=(S-M2projectile+M2target+std::sqrt(DecayMomentum2))/2./SqrtS;
496 WplusProjectile=SqrtS - M2target/WminusTarget;
497//-------------------------------------------------------------
498 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
499 {
500 G4Nucleon * aNucleon = TheInvolvedNucleon[i];
501 G4LorentzVector tmp=aNucleon->Get4Momentum();
502
503 G4double Mt2 = sqr(tmp.x())+sqr(tmp.y())+
504 aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass()*
505 aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass();
506 G4double Xminus=tmp.z();
507
508 G4double Pz=-WminusTarget*Xminus/2. + Mt2/(2.*WminusTarget*Xminus);
509 G4double E = WminusTarget*Xminus/2. + Mt2/(2.*WminusTarget*Xminus);
510
511 if( E+Pz > WplusProjectile ){OuterSuccess=false; break;}
512 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
513 } while(!OuterSuccess);
514
515//-------------------------------------------------------------
516 G4double Pzprojectile=WplusProjectile/2. - M2projectile/2./WplusProjectile;
517 G4double Eprojectile =WplusProjectile/2. + M2projectile/2./WplusProjectile;
518 Pprojectile.setPz(Pzprojectile); Pprojectile.setE(Eprojectile);
519
520 Pprojectile.transform(toLab); // The work with the projectile
521 primary->Set4Momentum(Pprojectile); // is finished at the moment.
522
523//-------------------------------------------------------------
524 G4ThreeVector Residual3Momentum(0.,0.,1.);
525
526 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
527 {
528 G4Nucleon * aNucleon = TheInvolvedNucleon[i];
529 G4LorentzVector tmp=aNucleon->Get4Momentum();
530 Residual3Momentum-=tmp.vect();
531
532 G4double Mt2 = sqr(tmp.x())+sqr(tmp.y())+
533 aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass()*
534 aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass();
535 G4double Xminus=tmp.z();
536
537 G4double Pz=-WminusTarget*Xminus/2. + Mt2/(2.*WminusTarget*Xminus);
538 G4double E = WminusTarget*Xminus/2. + Mt2/(2.*WminusTarget*Xminus);
539
540 tmp.setPz(Pz);
541 tmp.setE(E);
542
543 tmp.transform(toLab);
544
545 aNucleon->SetMomentum(tmp);
546
547 G4VSplitableHadron * targetSplitable=aNucleon->GetSplitableHadron();
548 targetSplitable->Set4Momentum(tmp);
549
550 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
551
552 G4double Mt2Residual=sqr(ResidualMass) +
553 sqr(Residual3Momentum.x())+sqr(Residual3Momentum.y());
554
555 G4double PzResidual=-WminusTarget*Residual3Momentum.z()/2. +
556 Mt2Residual/(2.*WminusTarget*Residual3Momentum.z());
557 G4double EResidual = WminusTarget*Residual3Momentum.z()/2. +
558 Mt2Residual/(2.*WminusTarget*Residual3Momentum.z());
559
560 Residual4Momentum.setPx(Residual3Momentum.x());
561 Residual4Momentum.setPy(Residual3Momentum.y());
562 Residual4Momentum.setPz(PzResidual);
563 Residual4Momentum.setE(EResidual);
564
565 Residual4Momentum.transform(toLab);
566//-------------------------------------------------------------
567 return true;
568}
569
570// ------------------------------------------------------------
571G4bool G4FTFModel::ExciteParticipants()
572{
573 G4bool Successfull(false);
574// do { // } while (Successfull == false) // Closed 15.12.09
575 Successfull=false;
576 theParticipants.StartLoop();
577 while (theParticipants.Next())
578 {
579 const G4InteractionContent & collision=theParticipants.GetInteraction();
580
581 G4VSplitableHadron * projectile=collision.GetProjectile();
582 G4VSplitableHadron * target=collision.GetTarget();
583
584 if(G4UniformRand()< theParameters->GetProbabilityOfElasticScatt())
585 { // Elastic scattering -------------------------
586 if(theElastic->ElasticScattering(projectile, target, theParameters))
587 {
588 Successfull = Successfull || true;
589 } else
590 {
591 Successfull = Successfull || false;
592 target->SetStatus(2);
593 }
594 }
595 else
596 { // Inelastic scattering ----------------------
597 if(theExcitation->ExciteParticipants(projectile, target,
598 theParameters, theElastic))
599 {
600 Successfull = Successfull || true;
601 } else
602 {
603 Successfull = Successfull || false;
604 target->SetStatus(2);
605 }
606 }
607 } // end of while (theParticipants.Next())
608// } while (Successfull == false); // Closed 15.12.09
609 return Successfull;
610}
611// ------------------------------------------------------------
612G4ExcitedStringVector * G4FTFModel::BuildStrings()
613{
614// Loop over all collisions; find all primaries, and all target ( targets may
615// be duplicate in the List ( to unique G4VSplitableHadrons)
616
617 G4ExcitedStringVector * strings;
618 strings = new G4ExcitedStringVector();
619
620 std::vector<G4VSplitableHadron *> primaries;
621
622 G4ExcitedString * FirstString(0); // If there will be a kink,
623 G4ExcitedString * SecondString(0); // two strings will be produced.
624
625 theParticipants.StartLoop(); // restart a loop
626
627 while ( theParticipants.Next() )
628 {
629 const G4InteractionContent & interaction=theParticipants.GetInteraction();
630 // do not allow for duplicates ...
631
632 if ( primaries.end() == std::find(primaries.begin(), primaries.end(),
633 interaction.GetProjectile()) )
634 primaries.push_back(interaction.GetProjectile());
635 }
636
637 unsigned int ahadron;
638 for ( ahadron=0; ahadron < primaries.size() ; ahadron++)
639 {
640 G4bool isProjectile(0);
641 if(primaries[ahadron]->GetStatus() == 1) {isProjectile=true; }
642 if(primaries[ahadron]->GetStatus() == 3) {isProjectile=false;}
643
644 FirstString=0; SecondString=0;
645 theExcitation->CreateStrings(primaries[ahadron], isProjectile,
646 FirstString, SecondString,
647 theParameters);
648
649 if(FirstString != 0) strings->push_back(FirstString);
650 if(SecondString != 0) strings->push_back(SecondString);
651 }
652
653 for (G4int ahadron=0; ahadron < NumberOfInvolvedNucleon ; ahadron++)
654 {
655 if(TheInvolvedNucleon[ahadron]->GetSplitableHadron()->GetStatus() !=0) //== 2)
656 {
657 G4bool isProjectile=false;
658 FirstString=0; SecondString=0;
659 theExcitation->CreateStrings(
660 TheInvolvedNucleon[ahadron]->GetSplitableHadron(),
661 isProjectile,
662 FirstString, SecondString,
663 theParameters);
664 if(FirstString != 0) strings->push_back(FirstString);
665 if(SecondString != 0) strings->push_back(SecondString);
666 }
667 }
668
669 std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron());
670 primaries.clear();
671
672 return strings;
673}
674// ------------------------------------------------------------
675void G4FTFModel::GetResidualNucleus()
676{ // This method is needed for the correct application of G4PrecompoundModelInterface
677 G4double DeltaExcitationE=ResidualExcitationEnergy/
678 (G4double) NumberOfInvolvedNucleon;
679 G4LorentzVector DeltaPResidualNucleus = Residual4Momentum/
680 (G4double) NumberOfInvolvedNucleon;
681
682 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
683 {
684 G4Nucleon * aNucleon = TheInvolvedNucleon[i];
685// G4LorentzVector tmp=aNucleon->Get4Momentum()-DeltaPResidualNucleus;
686 G4LorentzVector tmp=-DeltaPResidualNucleus;
687 aNucleon->SetMomentum(tmp);
688 aNucleon->SetBindingEnergy(DeltaExcitationE);
689 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
690
691}
692
693// ------------------------------------------------------------
694G4ThreeVector G4FTFModel::GaussianPt(G4double AveragePt2, G4double maxPtSquare) const
695{ // @@ this method is used in FTFModel as well. Should go somewhere common!
696
697 G4double Pt2(0.);
698 if(AveragePt2 <= 0.) {Pt2=0.;}
699 else
700 {
701 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() *
702 (std::exp(-maxPtSquare/AveragePt2)-1.));
703 }
704 G4double Pt=std::sqrt(Pt2);
705 G4double phi=G4UniformRand() * twopi;
706
707 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
708}
Note: See TracBrowser for help on using the repository browser.