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

Last change on this file since 1350 was 1347, checked in by garnier, 14 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.