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

Last change on this file since 1201 was 1196, checked in by garnier, 16 years ago

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