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

Last change on this file since 1198 was 1196, checked in by garnier, 15 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.