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

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

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