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

Last change on this file since 1343 was 1340, checked in by garnier, 15 years ago

update ti head
File size: 31.2 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.36 2010/09/20 15:50:46 vuzhinsk Exp $
28// GEANT4 tag $Name: geant4-09-03-ref-09 $
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
102 theParticipants.Init(aNucleus.GetA_asInt(),aNucleus.GetZ_asInt());
103
104// ----------- N-mass number Z-charge -------------------------
105
106// --- cms energy
107 G4double s = sqr( theProjectile.GetMass() ) +
108 sqr( G4Proton::Proton()->GetPDGMass() ) +
109 2*theProjectile.GetTotalEnergy()*G4Proton::Proton()->GetPDGMass();
110
111 if( theParameters != 0 ) delete theParameters;
112 theParameters = new G4FTFParameters(theProjectile.GetDefinition(),
113 aNucleus.GetN(),aNucleus.GetZ(),
114 s);
115
116//theParameters->SetProbabilityOfElasticScatt(0.);
117//G4cout<<theParameters->GetProbabilityOfElasticScatt()<<G4endl;
118//G4int Uzhi; G4cin>>Uzhi;
119// To turn on/off (1/0) elastic scattering
120
121}
122
123// ------------------------------------------------------------
124struct DeleteVSplitableHadron { void operator()(G4VSplitableHadron * aH){ delete aH;} };
125
126
127// ------------------------------------------------------------
128G4ExcitedStringVector * G4FTFModel::GetStrings()
129{
130 G4ExcitedStringVector * theStrings(0);
131//G4cout<<"GetString"<<G4endl;
132 theParticipants.GetList(theProjectile,theParameters);
133//G4cout<<"Reggeon"<<G4endl;
134 ReggeonCascade();
135
136 G4bool Success(true);
137 if( PutOnMassShell() )
138 {
139//G4cout<<"PutOn mass Shell OK"<<G4endl;
140 if( ExciteParticipants() )
141 {
142//G4cout<<"Excite partic OK"<<G4endl;
143 theStrings = BuildStrings();
144//G4cout<<"Build String OK"<<G4endl;
145 GetResidualNucleus();
146
147 if( theParameters != 0 )
148 {
149 delete theParameters;
150 theParameters=0;
151 }
152 } else // if( ExciteParticipants() )
153 { Success=false;}
154 } else // if( PutOnMassShell() )
155 { Success=false;}
156
157 if(!Success)
158 {
159 // -------------- Erase the projectile ----------------
160 std::vector<G4VSplitableHadron *> primaries;
161
162 theParticipants.StartLoop(); // restart a loop
163 while ( theParticipants.Next() )
164 {
165 const G4InteractionContent & interaction=theParticipants.GetInteraction();
166 // do not allow for duplicates ...
167 if ( primaries.end() == std::find(primaries.begin(), primaries.end(),
168 interaction.GetProjectile()) )
169 primaries.push_back(interaction.GetProjectile());
170 }
171 std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron());
172 primaries.clear();
173 }
174
175// -------------- Cleaning of the memory --------------
176// -------------- Erase the target nucleons -----------
177 G4VSplitableHadron * aNucleon = 0;
178 for(G4int i=0; i < NumberOfInvolvedNucleon; i++)
179 {
180 aNucleon = TheInvolvedNucleon[i]->GetSplitableHadron();
181 if(aNucleon) delete aNucleon;
182 }
183
184 NumberOfInvolvedNucleon=0;
185
186 return theStrings;
187
188}
189//-------------------------------------------------------------------
190void G4FTFModel::ReggeonCascade()
191{ //--- Implementation of reggeon theory inspired model-------
192 NumberOfInvolvedNucleon=0;
193
194 theParticipants.StartLoop();
195 while (theParticipants.Next())
196 {
197 const G4InteractionContent & collision=theParticipants.GetInteraction();
198 G4Nucleon * TargetNucleon=collision.GetTargetNucleon();
199
200 TheInvolvedNucleon[NumberOfInvolvedNucleon]=TargetNucleon;
201 NumberOfInvolvedNucleon++;
202//G4cout<<"Prim NumberOfInvolvedNucleon "<<NumberOfInvolvedNucleon<<G4endl;
203 G4double XofWoundedNucleon = TargetNucleon->GetPosition().x();
204 G4double YofWoundedNucleon = TargetNucleon->GetPosition().y();
205
206 theParticipants.theNucleus->StartLoop();
207 G4Nucleon * Neighbour(0);
208
209 while ( (Neighbour = theParticipants.theNucleus->GetNextNucleon()) )
210 {
211 if(!Neighbour->AreYouHit())
212 {
213 G4double impact2= sqr(XofWoundedNucleon - Neighbour->GetPosition().x()) +
214 sqr(YofWoundedNucleon - Neighbour->GetPosition().y());
215
216 if(G4UniformRand() < theParameters->GetCofNuclearDestruction()*
217 std::exp(-impact2/theParameters->GetR2ofNuclearDestruction()))
218 { // The neighbour nucleon is involved in the reggeon cascade
219
220 TheInvolvedNucleon[NumberOfInvolvedNucleon]=Neighbour;
221 NumberOfInvolvedNucleon++;
222//G4cout<<"Seco NumberOfInvolvedNucleon "<<NumberOfInvolvedNucleon<<G4endl;
223
224 G4VSplitableHadron *targetSplitable;
225 targetSplitable = new G4DiffractiveSplitableHadron(*Neighbour);
226
227 Neighbour->Hit(targetSplitable);
228 targetSplitable->SetStatus(2);
229 }
230 } // end of if(!Neighbour->AreYouHit())
231 } // end of while (theParticipant.theNucleus->GetNextNucleon())
232 } // end of while (theParticipants.Next())
233
234// ---------------- Calculation of creation time for each target nucleon -----------
235 theParticipants.StartLoop(); // restart a loop
236 theParticipants.Next();
237 G4VSplitableHadron * primary = theParticipants.GetInteraction().GetProjectile();
238 G4double betta_z=primary->Get4Momentum().pz()/primary->Get4Momentum().e();
239 primary->SetTimeOfCreation(0.);
240
241 G4double ZcoordinateOfPreviousCollision(0.);
242 G4double ZcoordinateOfCurrentInteraction(0.);
243 G4double TimeOfPreviousCollision(0.);
244 G4double TimeOfCurrentCollision(0);
245
246 theParticipants.theNucleus->StartLoop();
247 G4Nucleon * aNucleon;
248 G4bool theFirstInvolvedNucleon(true);
249 while ( (aNucleon = theParticipants.theNucleus->GetNextNucleon()) )
250 {
251 if(aNucleon->AreYouHit())
252 {
253 if(theFirstInvolvedNucleon)
254 {
255 ZcoordinateOfPreviousCollision=aNucleon->GetPosition().z();
256 theFirstInvolvedNucleon=false;
257 }
258
259 ZcoordinateOfCurrentInteraction=aNucleon->GetPosition().z();
260 TimeOfCurrentCollision=TimeOfPreviousCollision+
261 (ZcoordinateOfCurrentInteraction-ZcoordinateOfPreviousCollision)/betta_z;
262// It is assumed that the nucleons are ordered on increasing z-coordinate ------------
263 aNucleon->GetSplitableHadron()->SetTimeOfCreation(TimeOfCurrentCollision);
264
265 ZcoordinateOfPreviousCollision=ZcoordinateOfCurrentInteraction;
266 TimeOfPreviousCollision=TimeOfCurrentCollision;
267 } // end of if(aNucleon->AreYouHit())
268 } // end of while (theParticipant.theNucleus->GetNextNucleon())
269//
270// The algorithm can be improved, but it will be more complicated, and will require
271// changes in G4DiffractiveExcitation.cc and G4ElasticHNScattering.cc
272} // Uzhi 26 July 2009
273
274// ------------------------------------------------------------
275G4bool G4FTFModel::PutOnMassShell()
276{
277// -------------- Properties of the projectile ----------------
278 theParticipants.StartLoop(); // restart a loop
279 theParticipants.Next();
280 G4VSplitableHadron * primary = theParticipants.GetInteraction().GetProjectile();
281 G4LorentzVector Pprojectile=primary->Get4Momentum();
282
283//G4cout<<"Pprojectile "<<Pprojectile<<G4endl;
284// To get original projectile particle
285
286 if(Pprojectile.z() < 0.){return false;}
287
288 G4double Mprojectile = Pprojectile.mag();
289 G4double M2projectile = Pprojectile.mag2();
290//-------------------------------------------------------------
291 G4LorentzVector Psum = Pprojectile;
292 G4double SumMasses = Mprojectile + 20.*MeV; // 13.12.09
293 // Separation energy for projectile
294//G4cout<<"SumMasses Pr "<<SumMasses<<G4endl;
295//--------------- Target nucleus ------------------------------
296 G4V3DNucleus *theNucleus = GetWoundedNucleus();
297 G4Nucleon * aNucleon;
298 G4int ResidualMassNumber=theNucleus->GetMassNumber();
299 G4int ResidualCharge =theNucleus->GetCharge();
300 ResidualExcitationEnergy=0.;
301 G4LorentzVector PnuclearResidual(0.,0.,0.,0.);
302
303 G4double ExcitationEnergyPerWoundedNucleon=
304 theParameters->GetExcitationEnergyPerWoundedNucleon();
305
306 theNucleus->StartLoop();
307
308 while ((aNucleon = theNucleus->GetNextNucleon()))
309 {
310 if(aNucleon->AreYouHit())
311 { // Involved nucleons
312 Psum += aNucleon->Get4Momentum();
313 SumMasses += aNucleon->GetDefinition()->GetPDGMass();
314 SumMasses += 20.*MeV; // 13.12.09 Separation energy for a nucleon
315//G4cout<<"SumMasses Tr "<<SumMasses<<G4endl;
316 ResidualMassNumber--;
317 ResidualCharge-=(G4int) aNucleon->GetDefinition()->GetPDGCharge();
318 ResidualExcitationEnergy+=ExcitationEnergyPerWoundedNucleon;
319 }
320 else
321 { // Spectator nucleons
322 PnuclearResidual += aNucleon->Get4Momentum();
323 } // end of if(!aNucleon->AreYouHit())
324 } // end of while (theNucleus->GetNextNucleon())
325
326 Psum += PnuclearResidual;
327//G4cout<<"ResidualCharge ,ResidualMassNumber "<<ResidualCharge<<" "<<ResidualMassNumber<<G4endl;
328 G4double ResidualMass(0.);
329 if(ResidualMassNumber == 0)
330 {
331 ResidualMass=0.;
332 ResidualExcitationEnergy=0.;
333 }
334 else
335 {
336 ResidualMass=G4ParticleTable::GetParticleTable()->GetIonTable()->
337 GetIonMass(ResidualCharge ,ResidualMassNumber);
338 if(ResidualMassNumber == 1) {ResidualExcitationEnergy=0.;}
339 }
340
341// ResidualMass +=ResidualExcitationEnergy; // Will be given after checks
342//G4cout<<"SumMasses End ResidualMass "<<SumMasses<<" "<<ResidualMass<<G4endl;
343 SumMasses += ResidualMass;
344//G4cout<<"SumMasses + ResM "<<SumMasses<<G4endl;
345//G4cout<<"Psum "<<Psum<<G4endl;
346//-------------------------------------------------------------
347
348 G4double SqrtS=Psum.mag();
349 G4double S=Psum.mag2();
350
351//G4cout<<"SqrtS < SumMasses "<<SqrtS<<" "<<SumMasses<<G4endl;
352 if(SqrtS < SumMasses) {return false;} // It is impossible to simulate
353 // after putting nuclear nucleons
354 // on mass-shell
355
356 if(SqrtS < SumMasses+ResidualExcitationEnergy) {ResidualExcitationEnergy=0.;}
357
358 ResidualMass +=ResidualExcitationEnergy;
359 SumMasses +=ResidualExcitationEnergy;
360//G4cout<<"ResidualMass "<<ResidualMass<<" "<<SumMasses<<G4endl;
361//-------------------------------------------------------------
362// Sampling of nucleons what are transfered to delta-isobars --
363 G4int MaxNumberOfDeltas = (int)((SqrtS - SumMasses)/(400.*MeV));
364 G4int NumberOfDeltas(0);
365
366 if(theNucleus->GetMassNumber() != 1)
367 {
368 G4double ProbDeltaIsobar(0.); // 1. *** Can be set if it is needed
369 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
370 {
371 if((G4UniformRand() < ProbDeltaIsobar)&&(NumberOfDeltas < MaxNumberOfDeltas))
372 {
373 NumberOfDeltas++;
374 G4VSplitableHadron * targetSplitable=TheInvolvedNucleon[i]->GetSplitableHadron();
375 SumMasses-=targetSplitable->GetDefinition()->GetPDGMass();
376
377 G4int PDGcode = targetSplitable->GetDefinition()->GetPDGEncoding();
378 G4int newPDGcode = PDGcode/10; newPDGcode=newPDGcode*10+4; // Delta
379 G4ParticleDefinition* ptr =
380 G4ParticleTable::GetParticleTable()->FindParticle(newPDGcode);
381 targetSplitable->SetDefinition(ptr);
382 SumMasses+=targetSplitable->GetDefinition()->GetPDGMass();
383 }
384 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
385 } // end of if(theNucleus.GetMassNumber() != 1)
386//-------------------------------------------------------------
387
388 G4LorentzRotation toCms(-1*Psum.boostVector());
389 G4LorentzVector Ptmp=toCms*Pprojectile;
390 if ( Ptmp.pz() <= 0. )
391 { // "String" moving backwards in CMS, abort collision !!
392 //G4cout << " abort ColliDeleteVSplitableHadronsion!! " << G4endl;
393 return false;
394 }
395
396// toCms.rotateZ(-1*Ptmp.phi()); // Uzhi 5.12.09
397// toCms.rotateY(-1*Ptmp.theta()); // Uzhi 5.12.09
398
399 G4LorentzRotation toLab(toCms.inverse());
400
401//-------------------------------------------------------------
402//------- Ascribing of the involved nucleons Pt and Xminus ----
403 G4double Dcor = theParameters->GetDofNuclearDestruction()/
404 theNucleus->GetMassNumber();
405
406 G4double AveragePt2 = theParameters->GetPt2ofNuclearDestruction();
407 G4double maxPtSquare = theParameters->GetMaxPt2ofNuclearDestruction();
408//G4cout<<"Dcor "<<Dcor<<" AveragePt2 "<<AveragePt2<<G4endl;
409 G4double M2target(0.);
410 G4double WminusTarget(0.);
411 G4double WplusProjectile(0.);
412
413 G4int NumberOfTries(0);
414 G4double ScaleFactor(1.);
415 G4bool OuterSuccess(true);
416 do // while (!OuterSuccess)
417 {
418 OuterSuccess=true;
419
420 do // while (SqrtS < Mprojectile + std::sqrt(M2target))
421 { // while (DecayMomentum < 0.)
422
423 NumberOfTries++;
424//G4cout<<"NumberOfTries "<<NumberOfTries<<G4endl;
425 if(NumberOfTries == 100*(NumberOfTries/100)) // 100
426 { // At large number of tries it would be better to reduce the values
427 ScaleFactor/=2.;
428 Dcor *=ScaleFactor;
429 AveragePt2 *=ScaleFactor;
430 }
431
432 G4ThreeVector PtSum(0.,0.,0.);
433 G4double XminusSum(0.);
434 G4double Xminus(0.);
435 G4bool InerSuccess=true;
436
437 do // while(!InerSuccess);
438 {
439 InerSuccess=true;
440
441 PtSum =G4ThreeVector(0.,0.,0.);
442 XminusSum=0.;
443
444 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
445 {
446 G4Nucleon * aNucleon = TheInvolvedNucleon[i];
447
448 G4ThreeVector tmpPt = GaussianPt(AveragePt2, maxPtSquare);
449 PtSum += tmpPt;
450 G4ThreeVector tmpX=GaussianPt(Dcor*Dcor, 1.);
451 Xminus=tmpX.x();
452 XminusSum+=Xminus;
453
454 G4LorentzVector tmp(tmpPt.x(),tmpPt.y(),Xminus,0.);
455//G4cout<<"Inv i mom "<<i<<" "<<tmp<<G4endl;
456 aNucleon->SetMomentum(tmp);
457 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
458
459//---------------------------------------------------------------------------
460 G4double DeltaX(0.);
461 G4double DeltaY(0.);
462 G4double DeltaXminus(0.);
463
464//G4cout<<"ResidualMassNumber "<<ResidualMassNumber<<" "<<PtSum<<G4endl;
465 if(ResidualMassNumber == 0)
466 {
467 DeltaX = PtSum.x()/NumberOfInvolvedNucleon;
468 DeltaY = PtSum.y()/NumberOfInvolvedNucleon;
469 DeltaXminus = (XminusSum-1.)/NumberOfInvolvedNucleon;
470 }
471 else
472 {
473 DeltaXminus = -1./theNucleus->GetMassNumber();
474 }
475//G4cout<<"Dx y xmin "<<DeltaX<<" "<<DeltaY<<" "<<DeltaXminus<<G4endl;
476 XminusSum=1.;
477 M2target =0.;
478
479 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
480 {
481 G4Nucleon * aNucleon = TheInvolvedNucleon[i];
482
483 Xminus = aNucleon->Get4Momentum().pz() - DeltaXminus;
484 XminusSum-=Xminus;
485//G4cout<<" i X-sum "<<i<<" "<<Xminus<<" "<<XminusSum<<G4endl;
486 if(ResidualMassNumber == 0) // Uzhi 5.07.10
487 {
488 if((Xminus <= 0.) || (Xminus > 1.)) {InerSuccess=false; break;}
489 } else
490 {
491 if((Xminus <= 0.) || (Xminus > 1.) ||
492 (XminusSum <=0.) || (XminusSum > 1.)) {InerSuccess=false; break;}
493 } // Uzhi 5.07.10
494
495 G4double Px=aNucleon->Get4Momentum().px() - DeltaX;
496 G4double Py=aNucleon->Get4Momentum().py() - DeltaY;
497
498 M2target +=(aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass()*
499 aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass() +
500 Px*Px + Py*Py)/Xminus;
501
502 G4LorentzVector tmp(Px,Py,Xminus,0.);
503 aNucleon->SetMomentum(tmp);
504 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
505
506//G4cout<<"Rescale O.K."<<G4endl;
507
508 if(InerSuccess && (ResidualMassNumber != 0))
509 {
510 M2target +=(ResidualMass*ResidualMass + PtSum.mag2())/XminusSum;
511 }
512//G4cout<<"InerSuccess "<<InerSuccess<<G4endl;
513//G4int Uzhi;G4cin>>Uzhi;
514 } while(!InerSuccess);
515 } while (SqrtS < Mprojectile + std::sqrt(M2target));
516//-------------------------------------------------------------
517 G4double DecayMomentum2= S*S+M2projectile*M2projectile+M2target*M2target
518 -2.*S*M2projectile - 2.*S*M2target
519 -2.*M2projectile*M2target;
520
521 WminusTarget=(S-M2projectile+M2target+std::sqrt(DecayMomentum2))/2./SqrtS;
522 WplusProjectile=SqrtS - M2target/WminusTarget;
523//G4cout<<"DecayMomentum2 "<<DecayMomentum2<<G4endl;
524//-------------------------------------------------------------
525 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
526 {
527 G4Nucleon * aNucleon = TheInvolvedNucleon[i];
528 G4LorentzVector tmp=aNucleon->Get4Momentum();
529
530 G4double Mt2 = sqr(tmp.x())+sqr(tmp.y())+
531 aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass()*
532 aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass();
533 G4double Xminus=tmp.z();
534
535 G4double Pz=-WminusTarget*Xminus/2. + Mt2/(2.*WminusTarget*Xminus);
536 G4double E = WminusTarget*Xminus/2. + Mt2/(2.*WminusTarget*Xminus);
537
538 if( E+Pz > WplusProjectile ){OuterSuccess=false; break;}
539 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
540//G4int Uzhi;G4cin>>Uzhi;
541 } while(!OuterSuccess);
542
543//-------------------------------------------------------------
544 G4double Pzprojectile=WplusProjectile/2. - M2projectile/2./WplusProjectile;
545 G4double Eprojectile =WplusProjectile/2. + M2projectile/2./WplusProjectile;
546 Pprojectile.setPz(Pzprojectile); Pprojectile.setE(Eprojectile);
547
548 Pprojectile.transform(toLab); // The work with the projectile
549 primary->Set4Momentum(Pprojectile); // is finished at the moment.
550
551//-------------------------------------------------------------
552 G4ThreeVector Residual3Momentum(0.,0.,1.);
553
554 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
555 {
556 G4Nucleon * aNucleon = TheInvolvedNucleon[i];
557 G4LorentzVector tmp=aNucleon->Get4Momentum();
558 Residual3Momentum-=tmp.vect();
559
560 G4double Mt2 = sqr(tmp.x())+sqr(tmp.y())+
561 aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass()*
562 aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass();
563 G4double Xminus=tmp.z();
564
565 G4double Pz=-WminusTarget*Xminus/2. + Mt2/(2.*WminusTarget*Xminus);
566 G4double E = WminusTarget*Xminus/2. + Mt2/(2.*WminusTarget*Xminus);
567
568 tmp.setPz(Pz);
569 tmp.setE(E);
570
571 tmp.transform(toLab);
572
573 aNucleon->SetMomentum(tmp);
574
575 G4VSplitableHadron * targetSplitable=aNucleon->GetSplitableHadron();
576 targetSplitable->Set4Momentum(tmp);
577
578 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
579
580 G4double Mt2Residual=sqr(ResidualMass) +
581 sqr(Residual3Momentum.x())+sqr(Residual3Momentum.y());
582
583 G4double PzResidual=-WminusTarget*Residual3Momentum.z()/2. +
584 Mt2Residual/(2.*WminusTarget*Residual3Momentum.z());
585 G4double EResidual = WminusTarget*Residual3Momentum.z()/2. +
586 Mt2Residual/(2.*WminusTarget*Residual3Momentum.z());
587
588 Residual4Momentum.setPx(Residual3Momentum.x());
589 Residual4Momentum.setPy(Residual3Momentum.y());
590 Residual4Momentum.setPz(PzResidual);
591 Residual4Momentum.setE(EResidual);
592
593 Residual4Momentum.transform(toLab);
594//-------------------------------------------------------------
595 return true;
596}
597
598// ------------------------------------------------------------
599G4bool G4FTFModel::ExciteParticipants()
600{
601 G4bool Successfull(false);
602// do { // } while (Successfull == false) // Closed 15.12.09
603 Successfull=false;
604 theParticipants.StartLoop();
605
606G4int MaxNumOfInelCollisions=G4int(theParameters->GetMaxNumberOfCollisions());
607G4double NumberOfInel(0.);
608//
609if(MaxNumOfInelCollisions > 0)
610{ // Plab > Pbound, Normal application of FTF is possible
611 G4double ProbMaxNumber=theParameters->GetMaxNumberOfCollisions()-MaxNumOfInelCollisions;
612 if(G4UniformRand() < ProbMaxNumber) {MaxNumOfInelCollisions++;}
613 NumberOfInel=MaxNumOfInelCollisions;
614} else
615{ // Plab < Pbound, Normal application of FTF is impossible, low energy corrections
616 if(theParticipants.theNucleus->GetMassNumber() > 1)
617 {
618 NumberOfInel = theParameters->GetProbOfInteraction();
619 MaxNumOfInelCollisions = 1;
620 } else
621 { // Special case for hadron-nucleon interactions
622 NumberOfInel = 1.;
623 MaxNumOfInelCollisions = 1;
624 }
625} // end of if(MaxNumOfInelCollisions > 0)
626//
627 while (theParticipants.Next())
628 {
629 const G4InteractionContent & collision=theParticipants.GetInteraction();
630
631 G4VSplitableHadron * projectile=collision.GetProjectile();
632 G4VSplitableHadron * target=collision.GetTarget();
633//G4cout<<"ProbabilityOfElasticScatt "<<theParameters->GetProbabilityOfElasticScatt()<<G4endl;
634 if(G4UniformRand()< theParameters->GetProbabilityOfElasticScatt())
635 { // Elastic scattering -------------------------
636//G4cout<<"Elastic FTF"<<G4endl;
637 if(theElastic->ElasticScattering(projectile, target, theParameters))
638 {
639 Successfull = Successfull || true;
640 } else
641 {
642 Successfull = Successfull || false;
643 target->SetStatus(2);
644 }
645 }
646 else
647 { // Inelastic scattering ----------------------
648/*
649 if(theExcitation->ExciteParticipants(projectile, target,
650 theParameters, theElastic))
651 {
652 Successfull = Successfull || true;
653 } else
654 {
655 Successfull = Successfull || false;
656 target->SetStatus(2);
657 }
658*/
659//G4cout<<"InElastic FTF"<<G4endl;
660 if(G4UniformRand()< NumberOfInel/MaxNumOfInelCollisions)
661 {
662 if(theExcitation->ExciteParticipants(projectile, target,
663 theParameters, theElastic))
664 {
665 Successfull = Successfull || true;
666NumberOfInel--;
667 } else
668 {
669 Successfull = Successfull || false;
670 target->SetStatus(2);
671 }
672 } else // If NumOfInel
673 {
674 if(theElastic->ElasticScattering(projectile, target, theParameters))
675 {
676 Successfull = Successfull || true;
677 } else
678 {
679 Successfull = Successfull || false;
680 target->SetStatus(2);
681 }
682 } // end if NumOfInel
683 }
684 } // end of while (theParticipants.Next())
685// } while (Successfull == false); // Closed 15.12.09
686 return Successfull;
687}
688// ------------------------------------------------------------
689G4ExcitedStringVector * G4FTFModel::BuildStrings()
690{
691// Loop over all collisions; find all primaries, and all target ( targets may
692// be duplicate in the List ( to unique G4VSplitableHadrons)
693
694 G4ExcitedStringVector * strings;
695 strings = new G4ExcitedStringVector();
696
697 std::vector<G4VSplitableHadron *> primaries;
698
699 G4ExcitedString * FirstString(0); // If there will be a kink,
700 G4ExcitedString * SecondString(0); // two strings will be produced.
701
702 theParticipants.StartLoop(); // restart a loop
703//
704 while ( theParticipants.Next() )
705 {
706 const G4InteractionContent & interaction=theParticipants.GetInteraction();
707 // do not allow for duplicates ...
708
709 if ( primaries.end() == std::find(primaries.begin(), primaries.end(),
710 interaction.GetProjectile()) )
711 primaries.push_back(interaction.GetProjectile());
712 }
713
714 unsigned int ahadron;
715 for ( ahadron=0; ahadron < primaries.size() ; ahadron++)
716 {
717 G4bool isProjectile(0);
718 if(primaries[ahadron]->GetStatus() == 1) {isProjectile=true; }
719 if(primaries[ahadron]->GetStatus() == 3) {isProjectile=false;}
720
721 FirstString=0; SecondString=0;
722 theExcitation->CreateStrings(primaries[ahadron], isProjectile,
723 FirstString, SecondString,
724 theParameters);
725
726 if(FirstString != 0) strings->push_back(FirstString);
727 if(SecondString != 0) strings->push_back(SecondString);
728 }
729//
730 for (G4int ahadron=0; ahadron < NumberOfInvolvedNucleon ; ahadron++)
731 {
732 if(TheInvolvedNucleon[ahadron]->GetSplitableHadron()->GetStatus() !=0) //== 2)
733 {
734 G4bool isProjectile=false;
735 FirstString=0; SecondString=0;
736 theExcitation->CreateStrings(
737 TheInvolvedNucleon[ahadron]->GetSplitableHadron(),
738 isProjectile,
739 FirstString, SecondString,
740 theParameters);
741 if(FirstString != 0) strings->push_back(FirstString);
742 if(SecondString != 0) strings->push_back(SecondString);
743 }
744 }
745
746 std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron());
747 primaries.clear();
748 return strings;
749}
750// ------------------------------------------------------------
751void G4FTFModel::GetResidualNucleus()
752{ // This method is needed for the correct application of G4PrecompoundModelInterface
753 G4double DeltaExcitationE=ResidualExcitationEnergy/
754 (G4double) NumberOfInvolvedNucleon;
755 G4LorentzVector DeltaPResidualNucleus = Residual4Momentum/
756 (G4double) NumberOfInvolvedNucleon;
757
758 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
759 {
760 G4Nucleon * aNucleon = TheInvolvedNucleon[i];
761// G4LorentzVector tmp=aNucleon->Get4Momentum()-DeltaPResidualNucleus;
762 G4LorentzVector tmp=-DeltaPResidualNucleus;
763 aNucleon->SetMomentum(tmp);
764 aNucleon->SetBindingEnergy(DeltaExcitationE);
765 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
766
767}
768
769// ------------------------------------------------------------
770G4ThreeVector G4FTFModel::GaussianPt(G4double AveragePt2, G4double maxPtSquare) const
771{ // @@ this method is used in FTFModel as well. Should go somewhere common!
772
773 G4double Pt2(0.);
774 if(AveragePt2 <= 0.) {Pt2=0.;}
775 else
776 {
777 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() *
778 (std::exp(-maxPtSquare/AveragePt2)-1.));
779 }
780 G4double Pt=std::sqrt(Pt2);
781 G4double phi=G4UniformRand() * twopi;
782
783 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
784}
Note: See TracBrowser for help on using the repository browser.