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

Last change on this file since 1357 was 1347, checked in by garnier, 15 years ago

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