Ignore:
Timestamp:
Nov 25, 2009, 5:13:58 PM (15 years ago)
Author:
garnier
Message:

update CVS release candidate geant4.9.3.01

Location:
trunk/source/processes/hadronic/models/parton_string/diffraction
Files:
15 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/parton_string/diffraction/History

    r962 r1196  
    1 $Id: History,v 1.5 2008/12/09 10:43:47 vuzhinsk Exp $
     1$Id: History,v 1.22 2009/10/29 14:58:52 vuzhinsk Exp $
    22-------------------------------------------------------------------
    33
    … …  
    1515     * Please list in reverse chronological order (last date on top)
    1616     ---------------------------------------------------------------
     1729 October 2009, V. Uzhinsky (hadr-string-diff-V09-02-17)
     18   Warning meassage is erased for using slc5 compilator in FTFModel.cc
     19
     2025 October 2009, V. Uzhinsky (hadr-string-diff-V09-02-16) 
     21   Excitation energy calculation of nuclear residuals is implemenred.
     22   The interface to the G4PrecompoundModelInterface is improved.
     23----------------------------------------------------------
     246 October 2009, V. Uzhinsky (hadr-string-diff-V09-02-15)
     25   Compilation warning are erased.
     26----------------------------------------------------------
     275 October 2009, V. Uzhinsky (hadr-string-diff-V09-02-14)
     28   FTFP with tuned parameters of nuclear de-excitation.
     29
     30-----------------------------------------------------------
     3119 Sept., 2009, V. Uzhinsky (hadr-string-diff-V09-02-13)
     32   I have foggoten to open quasi-elastic scattering!
     33   Now it is openned.
     34----------------------------------------------------------
     3517 Sept., 2009, V. Uzhinsky (hadr-string-diff-V09-02-12)
     36   FTF parameters have been tunned on hN-interactions.
     37   Binary reactions were introduced and checked.
     38   Creation of s-channel Delta's was erased.
     39
     40-----------------------------------------
     416 August 2009, V. Uzhinsky (hadr-string-diff-V09-02-11)
     42  Kaon A interactions are improved in FTFModel.
     43----------------------------------------------------------
     446 August 2009, V. Uzhinsky (hadr-string-diff-V09-02-10)
     45  A warning message at FTFmodel compilation connected with Xminus
     46  initialisation was erased.
     47----------------------------------------------------------
     485 August 2009, V.Uzhinsky (hadr-string-diff-V09-02-09)
     49  Some additional warning found by Gabriela were erased.
     50------------------------------------------------------------
     515 August 2009, V. Uzhinsky (hadr-string-diff-V09-02-08)
     52  Warning meassages were erased at a compilation of FTF model
     53-----------------------------------------------------------
     543 August 2009, V. Uzhinsky (hadr-string-diff-V09-02-07)
     55   FTF model with the reggeon cascading was complited by
     56   formation time for a following coupling with the binary model.
     57-----------------------------------------------------------
     5831 July 2009, V. Uzhinsky (hadr-string-diff-V09-02-06)
     59   Inconsistency in G4FTFModel was erased thank to Gunter.
     60
     61----------------------------------------------------------
     6231 July 2009, V. Uzhinsky (hadr-string-diff-V09-02-05)
     63   The reggeon theory inspired model of duclear desctruction is
     64   implemented for a simulation of low particle cascading.
     65   The code can operate with Precompound model only.
     66-------------------------------------------------------------------
     67
     6817 July V. Uzhinsky        (hadr-string-diff-V09-02-04)
     69   A Status of nuclear nucleon involved in an interaction is introdused.
     70   Status: 0 - spectator, 1 - primary involved nucleon, 2 - a nucleon
     71   involved in the reggeon cascading, 3 - absorbed nucleon.
     72
     73   A connection between a participant nucleon and a nuclear nucleon was
     74   introsuced in G4InteractionContent.
     75
     76   All of these allow to improve FTF model for pion obsorption on nuclei.
     77   These required a lot of changes.
     78---------------------------------------------------------
     79
     8010 July V. Uzhinsky        (hadr-string-diff-V09-02-03)
     81Bug was fixed in G4FTFModel.cc thank to A. Ribon
     82---------------------------------------------------------
     839 July 2009 V. Uzhinsky    (hadr-string-diff-V09-02-02)
     84- Charge-exchange was implemented for pn->np
     85  elastic and inelastic interactions.
     86- Pion absorption by a nucleon was implemented.
     87  See comments marked by Uzhi 07.07.09
     88---------------------------------------------------------
     8919 December 08 V. Uzhinsky (hadr-string-diff-V09-02-01)
     90---------------------------------------------------------
     91 Version of FTF suted for low energies (2-10 GeV/c)
     92
    17939 December 08, V. Uzhinsky (hadr-string-diff-V09-01-04)
    1894- Improvement of delete operators in FTF

  • trunk/source/processes/hadronic/models/parton_string/diffraction/include/G4DiffractiveExcitation.hh

    r962 r1196  
    2525//
    2626//
    27 // $Id: G4DiffractiveExcitation.hh,v 1.2 2008/04/25 14:20:13 vuzhinsk Exp $
     27// $Id: G4DiffractiveExcitation.hh,v 1.4 2009/09/17 18:24:30 vuzhinsk Exp $
    2828
    2929#ifndef G4DiffractiveExcitation_h
    … …  
    4242class G4VSplitableHadron;
    4343class G4ExcitedString;
    44 #include "G4FTFParameters.hh"                            // Uzhi 19.04.08
     44#include "G4FTFParameters.hh"
     45#include "G4ElasticHNScattering.hh"  // Uzhi 3.09.09
    4546#include "G4ThreeVector.hh"
    4647
    … …  
    5051  public:
    5152
    52       G4DiffractiveExcitation();                           // Uzhi
     53      G4DiffractiveExcitation();
    5354      virtual ~G4DiffractiveExcitation();
    5455
    5556      virtual G4bool ExciteParticipants (G4VSplitableHadron *aPartner,
    5657                                         G4VSplitableHadron * bPartner,
     58                                         G4FTFParameters *theParameters,
     59                                         G4ElasticHNScattering *theElastic) const;
     60
     61      virtual void CreateStrings        (G4VSplitableHadron * aHadron,
     62                                         G4bool isProjectile,
     63                                         G4ExcitedString * &FirstString,
     64                                         G4ExcitedString * &SecondString,
    5765                                         G4FTFParameters *theParameters) const;
    58 
    59       virtual G4ExcitedString * String(G4VSplitableHadron * aHadron, G4bool isProjectile) const;
    60 
    6166  private:
    6267
    6368      G4DiffractiveExcitation(const G4DiffractiveExcitation &right);
    6469     
    65       G4ThreeVector GaussianPt(G4double  AveragePt2, G4double maxPtSquare) const; // Uzhi
    66       G4double ChooseP(G4double Pmin, G4double Pmax) const;                       // Uzhi
    67      
     70      G4ThreeVector GaussianPt(G4double  AveragePt2, G4double maxPtSquare) const;
     71      G4double ChooseP(G4double Pmin, G4double Pmax) const;
     72      G4double GetQuarkFractionOfKink(G4double zmin, G4double zmax) const;
     73      void UnpackMeson( G4int IdPDG, G4int &Q1, G4int &Q2) const;  // Uzhi 7.09.09
     74      void UnpackBaryon(G4int IdPDG, G4int &Q1, G4int &Q2, G4int &Q3) const; // Uzhi 7.09.09
     75      G4int NewNucleonId(G4int Q1, G4int Q2, G4int Q3) const; // Uzhi 7.09.09
     76
    6877      const G4DiffractiveExcitation & operator=(const G4DiffractiveExcitation &right);
    6978      int operator==(const G4DiffractiveExcitation &right) const;

  • trunk/source/processes/hadronic/models/parton_string/diffraction/include/G4DiffractiveHHScatterer.hh

    r962 r1196  
    2525//
    2626//
    27 // $Id: G4DiffractiveHHScatterer.hh,v 1.4 2008/04/25 14:20:13 vuzhinsk Exp $
     27// $Id: G4DiffractiveHHScatterer.hh,v 1.7 2009/10/06 10:10:36 vuzhinsk Exp $
    2828
    2929#ifndef G4DiffractiveHHScatterer_h
    … …  
    4444class G4KineticTrack;
    4545#include "G4KineticTrackVector.hh"
    46 #include "G4FTFParameters.hh"                            // Uzhi 21.04.08
     46#include "G4FTFParameters.hh"
     47#include "G4ExcitedString.hh"
    4748
    4849class G4DiffractiveHHScatterer
    … …  
    5253   G4DiffractiveHHScatterer();
    5354   
    54    G4KineticTrackVector * Scatter(const G4KineticTrack & aTrack, const G4KineticTrack & bTrack);
     55//   G4KineticTrackVector * Scatter(const G4KineticTrack & aTrack, const G4KineticTrack & bTrack);
    5556
     57      virtual void CreateStrings() const;
     58/*
     59                                        (G4VSplitableHadron * aHadron,
     60                                         G4bool isProjectile,
     61                                         G4ExcitedString * FirstString,
     62                                         G4ExcitedString * SecondString,
     63                                         G4FTFParameters *theParameters) const;
     64*/
    5665private:
    5766
    5867const G4DiffractiveExcitation * theExcitation;
    5968G4LundStringFragmentation * theStringFragmentation;
    60 G4FTFParameters  *theParameters;                       // Uzhi  21.04.08
     69G4FTFParameters  *theParameters;
    6170};
    6271

  • trunk/source/processes/hadronic/models/parton_string/diffraction/include/G4DiffractiveSplitableHadron.hh

    r1007 r1196  
    2525//
    2626//
    27 // $Id: G4DiffractiveSplitableHadron.hh,v 1.4 2006/06/29 20:54:27 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02 $
     27// $Id: G4DiffractiveSplitableHadron.hh,v 1.5 2009/08/03 13:14:19 vuzhinsk Exp $
     28// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2929//
    3030

  • trunk/source/processes/hadronic/models/parton_string/diffraction/include/G4ElasticHNScattering.hh

    r967 r1196  
    2525//
    2626//
    27 // $Id: G4ElasticHNScattering.hh,v 1.1 2008/03/31 15:34:01 vuzhinsk Exp $
     27// $Id: G4ElasticHNScattering.hh,v 1.2 2009/08/03 13:14:19 vuzhinsk Exp $
    2828
    2929#ifndef G4ElasticHNScattering_h
    … …  
    4242class G4VSplitableHadron;
    4343class G4ExcitedString;
    44 #include "G4FTFParameters.hh"                            // Uzhi 29.03.08
     44#include "G4FTFParameters.hh"
    4545#include "G4ThreeVector.hh"
    4646
    … …  
    5050  public:
    5151
    52       G4ElasticHNScattering();                           // Uzhi
     52      G4ElasticHNScattering();
    5353      virtual ~G4ElasticHNScattering();
    5454
    … …  
    6161      G4ElasticHNScattering(const G4ElasticHNScattering &right);
    6262     
    63       G4ThreeVector GaussianPt(G4double  AveragePt2, G4double maxPtSquare) const;  // Uzhi
     63      G4ThreeVector GaussianPt(G4double  AveragePt2, G4double maxPtSquare) const;
    6464     
    6565      const G4ElasticHNScattering & operator=(const G4ElasticHNScattering &right);

  • trunk/source/processes/hadronic/models/parton_string/diffraction/include/G4FTFModel.hh

    r1007 r1196  
    2525//
    2626//
    27 // $Id: G4FTFModel.hh,v 1.7 2008/04/25 14:20:13 vuzhinsk Exp $
    28 // GEANT4 tag $Name: geant4-09-02 $
     27// $Id: G4FTFModel.hh,v 1.10 2009/10/25 10:50:54 vuzhinsk Exp $
     28// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2929//
    3030// Class Description
    … …  
    5353class G4ExcitedString;
    5454
    55 #include "G4FTFParameters.hh"                            // Uzhi 29.03.08
     55#include "G4FTFParameters.hh"
    5656#include "G4FTFParticipants.hh"
    5757
    … …  
    6464
    6565  public:
    66       G4FTFModel();                                          // Uzhi
    67       G4FTFModel(G4double , G4double , G4double );           // Uzhi
     66      G4FTFModel();
     67      G4FTFModel(G4double , G4double , G4double );
    6868      G4FTFModel(G4DiffractiveExcitation * anExcitation);
    6969      G4FTFModel(const G4FTFModel &right);
    … …  
    8282 
    8383  private:
     84       void ReggeonCascade();
     85       G4bool PutOnMassShell();
    8486       G4bool ExciteParticipants();
    8587       G4ExcitedStringVector * BuildStrings();
     88       void GetResidualNucleus();                  // 23 Oct. 2009
     89       G4ThreeVector GaussianPt(G4double  AveragePt2, G4double maxPtSquare) const;
    8690 
    8791  private:     
    … …  
    8993       G4ReactionProduct theProjectile;       
    9094       G4FTFParticipants theParticipants;
     95       
     96       G4Nucleon * TheInvolvedNucleon[250];
     97       G4int NumberOfInvolvedNucleon;
    9198
    92        G4FTFParameters  *theParameters;        // Uzhi  29.03.08
     99       G4FTFParameters  *theParameters;
    93100       G4DiffractiveExcitation * theExcitation;
    94        G4ElasticHNScattering   * theElastic;   // Uzhi 29.03.08
     101       G4ElasticHNScattering   * theElastic;
    95102
     103       G4LorentzVector Residual4Momentum;
     104       G4double ResidualExcitationEnergy;
    96105
    97106};
    … …  
    105114
    106115#endif
    107 
    108 

  • trunk/source/processes/hadronic/models/parton_string/diffraction/include/G4FTFParameters.hh

    r1007 r1196  
    2727#define G4FTFParameters_h 1
    2828//
    29 // $Id: G4FTFParameters.hh,v 1.2 2008/06/13 12:49:23 vuzhinsk Exp $
    30 // GEANT4 tag $Name: geant4-09-02 $
     29// $Id: G4FTFParameters.hh,v 1.7 2009/10/25 10:50:54 vuzhinsk Exp $
     30// GEANT4 tag $Name: geant4-09-03-cand-01 $
    3131//
    3232#include "G4Proton.hh"
    … …  
    6060
    6161// --------- Set parameters of excitations ---------------------------
    62         void SetProjMinDiffMass(const G4double aValue); // Uzhi 19.04.08
     62        void SetMagQuarkExchange(const G4double aValue);
     63        void SetSlopeQuarkExchange(const G4double aValue);
     64        void SetDeltaProbAtQuarkExchange(const G4double aValue);
     65
     66        void SetProjMinDiffMass(const G4double aValue);
    6367        void SetProjMinNonDiffMass(const G4double aValue);
    6468        void SetProbabilityOfProjDiff(const G4double aValue);
    6569
    66         void SetTarMinDiffMass(const G4double aValue);  // Uzhi 19.04.08
     70        void SetTarMinDiffMass(const G4double aValue);
    6771        void SetTarMinNonDiffMass(const G4double aValue);
    6872        void SetProbabilityOfTarDiff(const G4double aValue);
    … …  
    7074        void SetAveragePt2(const G4double aValue);
    7175
    72    // Perevod mass*GeV Pt2*GeV*GeV
    73 
     76// --------- Set parameters of a string kink --------------------------------
     77        void SetPt2Kink(const G4double aValue);
     78        void SetQuarkProbabilitiesAtGluonSplitUp(const G4double Puubar,
     79                                                 const G4double Pddbar,
     80                                                 const G4double Pssbar );
     81
     82// --------- Set parameters of nuclear destruction--------------------
     83        void SetCofNuclearDestruction(const G4double aValue);
     84        void SetR2ofNuclearDestruction(const G4double aValue);
     85
     86        void SetExcitationEnergyPerWoundedNucleon(const G4double aValue);
     87
     88        void SetDofNuclearDestruction(const G4double aValue);
     89        void SetPt2ofNuclearDestruction(const G4double aValue);
     90        void SetMaxPt2ofNuclearDestruction(const G4double aValue);
     91
     92//--------------------------------------------------------------------
    7493// --------- Get geometrical parameteres -----------------------------
    7594        G4double GetTotalCrossSection();
    … …  
    87106
    88107// --------- Get parameters of excitations ---------------------------
    89         G4double GetProjMinDiffMass(); // Uzhi 19.04.08
     108        G4double GetMagQuarkExchange();
     109        G4double GetSlopeQuarkExchange();
     110        G4double GetDeltaProbAtQuarkExchange();
     111
     112        G4double GetProjMinDiffMass();
    90113        G4double GetProjMinNonDiffMass();
    91114        G4double GetProbabilityOfProjDiff();
    92115
    93         G4double GetTarMinDiffMass();  // Uzhi 19.04.08
     116        G4double GetTarMinDiffMass();
    94117        G4double GetTarMinNonDiffMass();
    95118        G4double GetProbabilityOfTarDiff();
    96119
    97 
    98120        G4double GetAveragePt2();
     121
     122// --------- Get parameters of a string kink --------------------------------
     123        G4double GetPt2Kink();
     124        std::vector<G4double>  GetQuarkProbabilitiesAtGluonSplitUp();
     125
     126// --------- Get parameters of nuclear destruction---------------------
     127        G4double GetCofNuclearDestruction();
     128        G4double GetR2ofNuclearDestruction();
     129
     130        G4double GetExcitationEnergyPerWoundedNucleon();
     131
     132        G4double GetDofNuclearDestruction();
     133        G4double GetPt2ofNuclearDestruction();
     134        G4double GetMaxPt2ofNuclearDestruction();
     135
    99136//  private:
    100137
    … …  
    115152
    116153// --------- Parameters of excitations -------------------------------
    117         G4double ProjMinDiffMass; // Uzhi 19.04.08
     154        G4double MagQuarkExchange;
     155        G4double SlopeQuarkExchange;
     156        G4double DeltaProbAtQuarkExchange;
     157
     158        G4double ProjMinDiffMass;
    118159        G4double ProjMinNonDiffMass;
    119160        G4double ProbabilityOfProjDiff;
    … …  
    124165
    125166        G4double AveragePt2;
     167
     168// ---------- Parameters of kink -------------------------------------
     169        G4double Pt2kink;
     170        std::vector<G4double> QuarkProbabilitiesAtGluonSplitUp;
     171
     172// --------- Parameters of nuclear destruction------------------------
     173        G4double CofNuclearDestruction;         // Cnd of nuclear destruction
     174        G4double R2ofNuclearDestruction;        // R2nd
     175
     176        G4double ExcitationEnergyPerWoundedNucleon;
     177
     178        G4double DofNuclearDestruction;         // D for momentum sampling
     179        G4double Pt2ofNuclearDestruction;       // Pt2
     180        G4double MaxPt2ofNuclearDestruction;    // Max Pt2
     181
    126182};
    127183
    … …  
    162218inline  void G4FTFParameters::SetAvaragePt2ofElasticScattering(const G4double aPt2)
    163219                 {
    164 //G4cout<<"Pt2 El "<<aPt2<<" "<<std::sqrt(aPt2)<<G4endl;
    165 //G4int Uzhi; G4cin>>Uzhi;
    166220AvaragePt2ofElasticScattering = aPt2;}
    167221
    168222// --------- Set parameters of excitations ----------------------------
    169 inline  void G4FTFParameters::SetProjMinDiffMass(const G4double aValue)   // Uzhi 19.04.08
     223inline  void G4FTFParameters::SetMagQuarkExchange(const G4double aValue)
     224             {MagQuarkExchange = aValue;}
     225inline  void G4FTFParameters::SetSlopeQuarkExchange(const G4double aValue)
     226             {SlopeQuarkExchange = aValue;}
     227inline  void G4FTFParameters::SetDeltaProbAtQuarkExchange(const G4double aValue)
     228             {DeltaProbAtQuarkExchange = aValue;}
     229
     230inline  void G4FTFParameters::SetProjMinDiffMass(const G4double aValue)
    170231             {ProjMinDiffMass = aValue*GeV;}
    171232inline  void G4FTFParameters::SetProjMinNonDiffMass(const G4double aValue)
    … …  
    174235             {ProbabilityOfProjDiff = aValue;}
    175236
    176 inline  void G4FTFParameters::SetTarMinDiffMass(const G4double aValue)  // Uzhi 19.04.08
     237inline  void G4FTFParameters::SetTarMinDiffMass(const G4double aValue)
    177238             {TarMinDiffMass = aValue*GeV;}
    178239inline  void G4FTFParameters::SetTarMinNonDiffMass(const G4double aValue)
    … …  
    184245             {AveragePt2 = aValue*GeV*GeV;}
    185246
     247// --------- Set parameters of a string kink --------------------------------
     248inline  void G4FTFParameters::SetPt2Kink(const G4double aValue)
     249             {Pt2kink = aValue;}
     250
     251inline  void G4FTFParameters::SetQuarkProbabilitiesAtGluonSplitUp(
     252                                                 const G4double Puubar,
     253                                                 const G4double Pddbar,
     254                                                 const G4double Pssbar )
     255             {
     256              QuarkProbabilitiesAtGluonSplitUp.push_back(Puubar);
     257              QuarkProbabilitiesAtGluonSplitUp.push_back(Puubar+Pddbar);
     258              QuarkProbabilitiesAtGluonSplitUp.push_back(Puubar+Pddbar+Pssbar);
     259             }
     260
     261// --------- Set parameters of nuclear destruction--------------------
     262inline  void G4FTFParameters::SetCofNuclearDestruction(const G4double aValue)
     263             {CofNuclearDestruction = aValue;}
     264inline  void G4FTFParameters::SetR2ofNuclearDestruction(const G4double aValue)
     265             {R2ofNuclearDestruction = aValue;}
     266
     267inline  void G4FTFParameters::SetExcitationEnergyPerWoundedNucleon(const G4double aValue)
     268             {ExcitationEnergyPerWoundedNucleon = aValue;}
     269
     270inline  void G4FTFParameters::SetDofNuclearDestruction(const G4double aValue)
     271             {DofNuclearDestruction = aValue;}
     272inline  void G4FTFParameters::SetPt2ofNuclearDestruction(const G4double aValue)
     273             {Pt2ofNuclearDestruction =aValue;}
     274inline  void G4FTFParameters::SetMaxPt2ofNuclearDestruction(const G4double aValue)
     275             {MaxPt2ofNuclearDestruction = aValue;}
     276
    186277// --------- Get geometrical parameteres ------------------------------
    187278inline  G4double G4FTFParameters::GetTotalCrossSection()     {return FTFXtotal;}
    … …  
    211302
    212303// --------- Get parameters of excitations ---------------------------
     304inline  G4double G4FTFParameters::GetMagQuarkExchange()       {return MagQuarkExchange;}
     305inline  G4double G4FTFParameters::GetSlopeQuarkExchange()     {return SlopeQuarkExchange;}
     306inline  G4double G4FTFParameters::GetDeltaProbAtQuarkExchange()
     307                                                              {return DeltaProbAtQuarkExchange;}
     308
     309
    213310inline  G4double G4FTFParameters::GetProjMinDiffMass()        {return ProjMinDiffMass;}
    214311inline  G4double G4FTFParameters::GetProjMinNonDiffMass()     {return ProjMinNonDiffMass;}
    … …  
    221318inline  G4double G4FTFParameters::GetAveragePt2()             {return AveragePt2;}
    222319
     320// --------- Get parameters of a string kink --------------------------
     321inline  G4double G4FTFParameters::GetPt2Kink()                {return Pt2kink;}
     322inline  std::vector<G4double>
     323                 G4FTFParameters::GetQuarkProbabilitiesAtGluonSplitUp()
     324                                  {return QuarkProbabilitiesAtGluonSplitUp;}
     325
     326// --------- Get parameters of nuclear destruction---------------------
     327inline  G4double G4FTFParameters::GetCofNuclearDestruction(){return CofNuclearDestruction;}
     328inline  G4double G4FTFParameters::GetR2ofNuclearDestruction(){return R2ofNuclearDestruction;}
     329
     330inline  G4double G4FTFParameters::GetExcitationEnergyPerWoundedNucleon()
     331            {return ExcitationEnergyPerWoundedNucleon;}
     332
     333
     334inline  G4double G4FTFParameters::GetDofNuclearDestruction()
     335                 {return DofNuclearDestruction;}
     336inline  G4double G4FTFParameters::GetPt2ofNuclearDestruction(){return Pt2ofNuclearDestruction;}
     337inline  G4double G4FTFParameters::GetMaxPt2ofNuclearDestruction()
     338                 {return MaxPt2ofNuclearDestruction;}
    223339#endif

  • trunk/source/processes/hadronic/models/parton_string/diffraction/include/G4FTFParticipants.hh

    r1007 r1196  
    2525//
    2626//
    27 // $Id: G4FTFParticipants.hh,v 1.5 2008/03/31 15:34:01 vuzhinsk Exp $
    28 // GEANT4 tag $Name: geant4-09-02 $
     27// $Id: G4FTFParticipants.hh,v 1.6 2009/08/03 13:14:19 vuzhinsk Exp $
     28// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2929//
    3030
    … …  
    4141
    4242#include "G4VParticipants.hh"
    43 #include "G4FTFParameters.hh"   // Uzhi 29.03.08
     43#include "G4FTFParameters.hh"
    4444#include <vector>
    4545#include "G4Nucleon.hh"
    … …  
    6262
    6363      void GetList(const G4ReactionProduct  &thePrimary,
    64                          G4FTFParameters    *theParameters); // Uzhi 29.03.08
     64                         G4FTFParameters    *theParameters);
    6565
    6666      void StartLoop();

  • trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4DiffractiveExcitation.cc

    r962 r1196  
    2525//
    2626//
    27 // $Id: G4DiffractiveExcitation.cc,v 1.7 2008/12/18 13:01:58 gunter Exp $
     27// $Id: G4DiffractiveExcitation.cc,v 1.16 2009/10/06 10:10:36 vuzhinsk Exp $
    2828// ------------------------------------------------------------
    2929//      GEANT 4 class implemetation file
    … …  
    4848
    4949#include "G4DiffractiveExcitation.hh"
     50#include "G4FTFParameters.hh"
     51#include "G4ElasticHNScattering.hh"
     52
    5053#include "G4LorentzRotation.hh"
     54#include "G4RotationMatrix.hh"
    5155#include "G4ThreeVector.hh"
    52 #include "G4ParticleDefinition.hh"
     56#include "G4ParticleDefinition.hh" 
    5357#include "G4VSplitableHadron.hh"
    5458#include "G4ExcitedString.hh"
    55 #include "G4FTFParameters.hh"                            // Uzhi 19.04.08
     59#include "G4ParticleTable.hh"
     60#include "G4Neutron.hh"
     61#include "G4ParticleDefinition.hh"
     62
    5663//#include "G4ios.hh"
    5764//#include "UZHI_diffraction.hh"
    … …  
    6572  ExciteParticipants(G4VSplitableHadron *projectile,
    6673                     G4VSplitableHadron *target,
    67                      G4FTFParameters    *theParameters) const
     74                     G4FTFParameters    *theParameters,
     75                     G4ElasticHNScattering *theElastic) const  // Uzhi 03.09.09
    6876{
    69      G4bool PutOnMassShell=0;
    70 
    7177// -------------------- Projectile parameters -----------------------
    72 
    7378     G4LorentzVector Pprojectile=projectile->Get4Momentum();
     79//G4cout<<"Excite P "<<Pprojectile<<G4endl;
     80
     81     if(Pprojectile.z() < 0.)
     82     {
     83       target->SetStatus(2);
     84       return false;
     85     }
     86
     87     G4double ProjectileRapidity = Pprojectile.rapidity();
     88
     89     G4int    ProjectilePDGcode=projectile->GetDefinition()->GetPDGEncoding();
     90//     G4ParticleDefinition * ProjectileDefinition = projectile->GetDefinition();
     91     G4int    absProjectilePDGcode=std::abs(ProjectilePDGcode);
     92
     93     G4bool PutOnMassShell(false);
    7494//   G4double M0projectile=projectile->GetDefinition()->GetPDGMass(); // With de-excitation
    7595     G4double M0projectile = Pprojectile.mag();                       // Without de-excitation
    76 /*
    77 G4cout<<"ExciteParticipants-------------------"<<G4endl;
    78 G4cout<<"Mom "<<Pprojectile<<" mass "<<M0projectile<<G4endl;
    79 */
     96
    8097     if(M0projectile < projectile->GetDefinition()->GetPDGMass())
    8198     {
    82       PutOnMassShell=1;
     99      PutOnMassShell=true;
    83100      M0projectile=projectile->GetDefinition()->GetPDGMass();
    84101     }
    85102
    86103     G4double M0projectile2 = M0projectile * M0projectile;
    87 
    88      G4int    PDGcode=projectile->GetDefinition()->GetPDGEncoding();
    89      G4int    absPDGcode=std::abs(PDGcode);
    90104
    91105     G4double ProjectileDiffStateMinMass=theParameters->GetProjMinDiffMass();
    92106     G4double ProjectileNonDiffStateMinMass=theParameters->GetProjMinNonDiffMass();
    93107     G4double ProbProjectileDiffraction=theParameters->GetProbabilityOfProjDiff();
    94 /*
    95 G4cout<<ProjectileDiffStateMinMass<<" "<<ProjectileNonDiffStateMinMass<<" "<<ProbProjectileDiffraction<<G4endl;
    96 */
    97 // -------------------- Target paraExciteParticipantsmeters -------------------------
     108
     109// -------------------- Target parameters -------------------------
     110     G4int    TargetPDGcode=target->GetDefinition()->GetPDGEncoding();
     111     G4int    absTargetPDGcode=std::abs(TargetPDGcode);
     112
    98113     G4LorentzVector Ptarget=target->Get4Momentum();
     114
    99115     G4double M0target = Ptarget.mag();
    100116
    101 //G4cout<<"Mom "<<Ptarget<<" mass "<<M0target<<G4endl;
     117//G4cout<<"Excite T "<<Ptarget<<G4endl;
     118//G4cout<<"PDGs "<<ProjectilePDGcode<<" "<<TargetPDGcode<<G4endl;
     119//G4cout<<"Masses "<<M0projectile<<" "<<M0target<<G4endl;
     120
     121     G4double TargetRapidity = Ptarget.rapidity();
    102122
    103123     if(M0target < target->GetDefinition()->GetPDGMass())
    104124     {
    105       PutOnMassShell=1;
     125      PutOnMassShell=true;
    106126      M0target=target->GetDefinition()->GetPDGMass();
    107127     }
    108128
    109      G4double M0target2 = M0target * M0target;             //Ptarget.mag2();
    110                                                            // for AA-inter.
     129     G4double M0target2 = M0target * M0target;
     130 
    111131     G4double TargetDiffStateMinMass=theParameters->GetTarMinDiffMass();   
    112132     G4double TargetNonDiffStateMinMass=theParameters->GetTarMinNonDiffMass();   
    113133     G4double ProbTargetDiffraction=theParameters->GetProbabilityOfTarDiff();
    114 /*
    115 G4cout<<TargetDiffStateMinMass<<" "<<TargetNonDiffStateMinMass<<" "<<ProbTargetDiffraction<<G4endl;
    116 */
     134
    117135     G4double AveragePt2=theParameters->GetAveragePt2();
     136
     137     G4double ProbOfDiffraction=ProbProjectileDiffraction +
     138                                ProbTargetDiffraction;
     139
     140     G4double SumMasses=M0projectile+M0target+200.*MeV;
    118141
    119142// Kinematical properties of the interactions --------------
    … …  
    122145     G4double S=Psum.mag2();
    123146
    124 //G4cout<<" sqrt(s) "<<std::sqrt(S)<<G4endl;
    125 
    126 // ------------------------------------------------------------------
    127 
    128 //ProbProjectileDiffraction=1.;
    129 //ProbTargetDiffraction    =1.;
    130      G4double ProbOfDiffraction=ProbProjectileDiffraction +
    131                                 ProbTargetDiffraction;
    132 
    133      if(ProbOfDiffraction!=0.)
    134      {
    135       ProbProjectileDiffraction/=ProbOfDiffraction;
    136      }
    137      else
    138      {
    139       ProbProjectileDiffraction=0.;
    140      }
    141 //   ProbTargetDiffraction    /=ProbOfDiffraction;
    142 
    143 //G4cout<<"ProbOfDiffraction "<<ProbOfDiffraction<<"ProbProjectileDiffraction "<<ProbProjectileDiffraction<<G4endl;   // Vova
    144 
    145      G4double ProjectileDiffStateMinMass2    = ProjectileDiffStateMinMass    *
    146                                                ProjectileDiffStateMinMass;
    147      G4double ProjectileNonDiffStateMinMass2 = ProjectileNonDiffStateMinMass *
    148                                                ProjectileNonDiffStateMinMass;
    149 
    150      G4double TargetDiffStateMinMass2        = TargetDiffStateMinMass        *
    151                                                TargetDiffStateMinMass;
    152      G4double TargetNonDiffStateMinMass2     = TargetNonDiffStateMinMass     *
    153                                                TargetNonDiffStateMinMass;
    154 
    155147// Transform momenta to cms and then rotate parallel to z axis;
    156 
    157 //         G4LorentzVector Psum;
    158 //         Psum=Pprojectile+Ptarget;
    159 
    160148     G4LorentzRotation toCms(-1*Psum.boostVector());
    161149
    162150     G4LorentzVector Ptmp=toCms*Pprojectile;
    163151     if ( Ptmp.pz() <= 0. )
    164         {
    165            // "String" moving backwards in  CMS, abort collision !!
    166            //G4cout << " abort Collision!! " << G4endl;
    167           return false;
    168         }
     152     {
     153       target->SetStatus(2);
     154   // "String" moving backwards in  CMS, abort collision !!
     155       return false;
     156    }
    169157
    170158     toCms.rotateZ(-1*Ptmp.phi());
    … …  
    176164     Ptarget.transform(toCms);
    177165
    178      G4double Pt2;
    179      G4double ProjMassT2, ProjMassT;
    180      G4double TargMassT2, TargMassT;
    181166     G4double PZcms2, PZcms;
    182      G4double PMinusMin, PMinusMax;
    183 //     G4double PPlusMin , PPlusMax;
    184      G4double TPlusMin , TPlusMax;
    185      G4double PMinusNew, PPlusNew, TPlusNew, TMinusNew;
    186 
    187 //   G4double S=Psum.mag2();
     167
    188168     G4double SqrtS=std::sqrt(S);
    189 
    190      if(absPDGcode > 1000 && SqrtS < 2200*MeV)
    191      {return false;}                         // The model cannot work for
     169//G4cout<<" SqrtS "<<SqrtS<<" 2200 "<<G4endl;
     170     if(absProjectilePDGcode > 1000 && SqrtS < 2300*MeV && SqrtS < SumMasses)
     171     {target->SetStatus(2);  return false;}  // The model cannot work for
    192172                                             // p+p-interactions
    193                                              // at Plab < 1.3 GeV/c.
    194 
    195      if(( absPDGcode == 211 || PDGcode ==  111) && SqrtS < 1600*MeV)
    196      {return false;}                         // The model cannot work for
     173                                             // at Plab < 1.62 GeV/c.
     174
     175     if(( absProjectilePDGcode == 211 || ProjectilePDGcode ==  111) &&
     176        (SqrtS < 1600*MeV) && (SqrtS < SumMasses))
     177     {target->SetStatus(2);  return false;}  // The model cannot work for
    197178                                             // Pi+p-interactions
    198179                                             // at Plab < 1. GeV/c.
    199180
    200      if(( absPDGcode == 321 || PDGcode == -311) && SqrtS < 1600*MeV)
    201      {return false;}                         // The model cannot work for
     181     if(( (absProjectilePDGcode == 321) || (ProjectilePDGcode == -311)   ||
     182          (absProjectilePDGcode == 311) || (absProjectilePDGcode == 130) ||
     183          (absProjectilePDGcode == 310)) && (SqrtS < 1600*MeV) && (SqrtS < SumMasses))
     184     {target->SetStatus(2);  return false;}  // The model cannot work for
    202185                                             // K+p-interactions
    203186                                             // at Plab < ??? GeV/c.  ???
    … …  
    208191
    209192     if(PZcms2 < 0)
    210      {return false;}   // It can be in an interaction with off-shell nuclear nucleon
     193     {target->SetStatus(2);  return false;}   // It can be in an interaction with
     194                                              // off-shell nuclear nucleon
    211195
    212196     PZcms = std::sqrt(PZcms2);
    … …  
    236220
    237221     G4double maxPtSquare; // = PZcms2;
    238 /*
    239 G4cout << "Pprojectile aft boost : " << Pprojectile <<" "<<Pprojectile.mag()<< G4endl;
    240 G4cout << "Ptarget aft boost : " << Ptarget <<" "<<Ptarget.mag()<< G4endl;
    241 G4cout << "cms aft boost : " << (Pprojectile+ Ptarget) << G4endl;
    242 G4cout << " Projectile Xplus / Xminus : " <<
    243             Pprojectile.plus() << " / " << Pprojectile.minus() << G4endl;
    244 G4cout << " Target Xplus / Xminus : " <<           Ptarget.plus() << " / " << Ptarget.minus() << G4endl;
    245 G4cout<<"maxPtSquare "<<maxPtSquare<<G4endl;
    246 */
     222
     223// Charge exchange can be possible for baryons -----------------
     224
     225//G4cout<<1./std::exp(-2.0*(ProjectileRapidity - TargetRapidity))<<G4endl;
     226//G4int Uzhi; G4cin>>Uzhi;
     227// Getting the values needed for exchange ----------------------
     228     G4double MagQuarkExchange        =theParameters->GetMagQuarkExchange();
     229     G4double SlopeQuarkExchange      =theParameters->GetSlopeQuarkExchange();
     230     G4double DeltaProbAtQuarkExchange=theParameters->GetDeltaProbAtQuarkExchange();
     231
     232//     G4double NucleonMass=
     233              (G4ParticleTable::GetParticleTable()->FindParticle(2112))->GetPDGMass();     
     234     G4double DeltaMass=
     235              (G4ParticleTable::GetParticleTable()->FindParticle(2224))->GetPDGMass();
     236
     237// Check for possible quark excjane -----------------------------------
     238     if(G4UniformRand() < MagQuarkExchange*
     239        std::exp(-SlopeQuarkExchange*(ProjectileRapidity - TargetRapidity)))
     240     {   
     241
     242//G4cout<<"Exchange "<<G4endl;
     243
     244      G4int NewProjCode(0), NewTargCode(0);
     245
     246//G4bool StandardExcitation(false); // =================================
     247      G4int ProjQ1(0), ProjQ2(0), ProjQ3(0);
     248
     249//  Projectile unpacking --------------------------
     250      if(absProjectilePDGcode < 1000 )
     251      {    // projectile is meson -----------------
     252       UnpackMeson(ProjectilePDGcode, ProjQ1, ProjQ2); 
     253      } else
     254      {    // projectile is baryon ----------------
     255       UnpackBaryon(ProjectilePDGcode, ProjQ1, ProjQ2, ProjQ3);
     256      } // End of the hadron's unpacking ----------
     257
     258//  Target unpacking ------------------------------
     259      G4int TargQ1(0), TargQ2(0), TargQ3(0);
     260      UnpackBaryon(TargetPDGcode, TargQ1, TargQ2, TargQ3);
     261
     262// Sampling of exchanged quarks -------------------
     263      G4int ProjExchangeQ(0);
     264      G4int TargExchangeQ(0);
     265
     266//G4cout<<"Targ "<<TargQ1<<" "<<TargQ2<<" "<<TargQ3<<G4endl;
     267      if(absProjectilePDGcode < 1000 )
     268      {    // projectile is meson -----------------
     269//G4cout<<"Proj Q1 Q2 "<<ProjQ1<<" "<<ProjQ2<<G4endl;
     270
     271       if(ProjQ1 > 0 ) // ProjQ1 is quark
     272       { 
     273        ProjExchangeQ = ProjQ1;
     274        if(ProjExchangeQ != TargQ1)
     275        {
     276         TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjQ1=TargExchangeQ;
     277        } else
     278        if(ProjExchangeQ != TargQ2)
     279        {
     280         TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjQ1=TargExchangeQ;
     281        } else         
     282        {
     283         TargExchangeQ = TargQ3;  TargQ3=ProjExchangeQ; ProjQ1=TargExchangeQ;
     284        }
     285//G4cout<<" Pr Tr "<<ProjQ1<<" "<<TargQ1<<" "<<TargQ2<<" "<<TargQ3<<G4endl;
     286       } else          // ProjQ2 is quark
     287       { 
     288        ProjExchangeQ = ProjQ2;
     289        if(ProjExchangeQ != TargQ1)
     290        {
     291         TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjQ2=TargExchangeQ;
     292        } else
     293        if(ProjExchangeQ != TargQ2)
     294        {
     295         TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjQ2=TargExchangeQ;
     296        } else
     297        {
     298         TargExchangeQ = TargQ3;  TargQ3=ProjExchangeQ; ProjQ2=TargExchangeQ;
     299        }
     300
     301       } // End of if(ProjQ1 > 0 ) // ProjQ1 is quark
     302
     303       G4int aProjQ1=std::abs(ProjQ1);
     304       G4int aProjQ2=std::abs(ProjQ2);
     305       if(aProjQ1 == aProjQ2)          {NewProjCode = 111;} // Pi0-meson
     306       else  // |ProjQ1| # |ProjQ2|
     307       {
     308        if(aProjQ1 > aProjQ2)          {NewProjCode = aProjQ1*100+aProjQ2*10+1;}
     309        else                           {NewProjCode = aProjQ2*100+aProjQ1*10+1;}
     310       }
     311
     312//       projectile->SetDefinition(
     313//       G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode));
     314
     315       NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3);
     316 
     317       if( (TargQ1 == TargQ2) && (TargQ1 == TargQ3) &&
     318           (SqrtS > M0projectile+DeltaMass))           {NewTargCode +=2;} //Create Delta isobar
     319
     320       else if( target->GetDefinition()->GetPDGiIsospin() == 3 )         //Delta was the target
     321       { if(G4UniformRand() > DeltaProbAtQuarkExchange){NewTargCode +=2;} //Save   Delta isobar
     322         else                                          {}     // De-excite initial Delta isobar
     323       }
     324
     325       else if((G4UniformRand() < DeltaProbAtQuarkExchange) &&         //Nucleon was the target
     326               (SqrtS > M0projectile+DeltaMass))      {NewTargCode +=2;}  //Create Delta isobar
     327       else                                           {}                 //Save initial nucleon
     328
     329       target->SetDefinition(
     330       G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode));
     331//G4cout<<"New PDGs "<<NewProjCode<<" "<<NewTargCode<<G4endl;
     332//G4int Uzhi; G4cin>>Uzhi;
     333       //}
     334      } else
     335      {    // projectile is baryon ----------------
     336G4double Same=0.; //0.3; //0.5;
     337       G4bool ProjDeltaHasCreated(false);
     338       G4bool TargDeltaHasCreated(false);
     339 
     340       G4double Ksi=G4UniformRand();
     341       if(G4UniformRand() < 0.5)     // Sampling exchange quark from proj. or targ.
     342       {   // Sampling exchanged quark from the projectile ---
     343//G4cout<<"Proj Exc"<<G4endl;
     344        if( Ksi < 0.333333 )
     345        {ProjExchangeQ = ProjQ1;}
     346        else if( (0.333333 <= Ksi) && (Ksi < 0.666667))
     347        {ProjExchangeQ = ProjQ2;}
     348        else
     349        {ProjExchangeQ = ProjQ3;}
     350//G4cout<<"ProjExchangeQ "<<ProjExchangeQ<<G4endl;
     351
     352        if((ProjExchangeQ != TargQ1)||(G4UniformRand()<Same)) // Vova
     353        {
     354         TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjExchangeQ=TargExchangeQ;
     355        } else
     356        if((ProjExchangeQ != TargQ2)||(G4UniformRand()<Same)) // Vova
     357        {
     358         TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjExchangeQ=TargExchangeQ;
     359        } else
     360        {
     361         TargExchangeQ = TargQ3;  TargQ3=ProjExchangeQ; ProjExchangeQ=TargExchangeQ;
     362        }
     363//G4cout<<"ProjExchangeQ "<<ProjExchangeQ<<G4endl;
     364        if( Ksi < 0.333333 )
     365        {ProjQ1=ProjExchangeQ;}
     366        else if( (0.333333 <= Ksi) && (Ksi < 0.666667))
     367        {ProjQ2=ProjExchangeQ;}
     368        else
     369        {ProjQ3=ProjExchangeQ;}
     370
     371/*
     372        NewProjCode = NewNucleonId(ProjQ1, ProjQ2, ProjQ3); // *****************************
     373
     374        if((ProjQ1==ProjQ2) && (ProjQ1==ProjQ3)) {NewProjCode +=2; ProjDeltaHasCreated=true;}
     375        else if(projectile->GetDefinition()->GetPDGiIsospin() == 3)// Projectile was Delta
     376        { if(G4UniformRand() > DeltaProbAtQuarkExchange)
     377                                                 {NewProjCode +=2; ProjDeltaHasCreated=true;}
     378          else                                   {NewProjCode +=0; ProjDeltaHasCreated=false;}
     379        }
     380        else                                                       // Projectile was Nucleon
     381        {
     382         if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > DeltaMass+M0target))
     383                                                 {NewProjCode +=2; ProjDeltaHasCreated=true;}
     384         else                                    {NewProjCode +=0; ProjDeltaHasCreated=false;}
     385        }
     386
     387        NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3); // *****************************
     388
     389        if((TargQ1==TargQ2) && (TargQ1==TargQ3)) {NewTargCode +=2; TargDeltaHasCreated=true;} 
     390        else if(target->GetDefinition()->GetPDGiIsospin() == 3)    // Target was Delta
     391        { if(G4UniformRand() > DeltaProbAtQuarkExchange)
     392                                                 {NewTargCode +=2; TargDeltaHasCreated=true;}
     393          else                                   {NewTargCode +=0; TargDeltaHasCreated=false;}
     394        }
     395        else                                                       // Target was Nucleon
     396        {
     397         if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > M0projectile+DeltaMass))
     398                                                 {NewTargCode +=2; TargDeltaHasCreated=true;}
     399         else                                    {NewTargCode +=0; TargDeltaHasCreated=false;}
     400        }         
     401
     402        if((absProjectilePDGcode == NewProjCode) && (absTargetPDGcode == NewTargCode))
     403        { // Nothing was changed! It is not right!?
     404        }
     405*/
     406/*-----------------------------------------------------------------------------------------
     407        if(!ProjDeltaHasCreated && !TargDeltaHasCreated) // No Deltas were created then
     408        { if( G4UniformRand() < 0.5) {NewProjCode +=2; ProjDeltaHasCreated = true;}
     409          else                       {NewTargCode +=2; TargDeltaHasCreated = true;}
     410        }
     411
     412        if(!(ProjDeltaHasCreated && TargDeltaHasCreated))
     413        {
     414         if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS >2600.))
     415         {
     416G4cout<<"2 Deltas "<<G4endl;
     417          if(!ProjDeltaHasCreated)
     418               {NewProjCode +=2; ProjDeltaHasCreated = true;}
     419          else {NewTargCode +=2; TargDeltaHasCreated = true;} // For Delta isobar
     420         }
     421        }
     422*/
     423       } else
     424       {   // Sampling exchanged quark from the target -------
     425//G4cout<<"Targ Exc"<<G4endl;
     426        if( Ksi < 0.333333 )
     427        {TargExchangeQ = TargQ1;}
     428        else if( (0.333333 <= Ksi) && (Ksi < 0.666667))
     429        {TargExchangeQ = TargQ2;}
     430        else
     431        {TargExchangeQ = TargQ3;}
     432//G4cout<<"TargExchangeQ "<<TargExchangeQ<<G4endl;
     433        if((TargExchangeQ != ProjQ1)||(G4UniformRand()<Same)) // Vova
     434        {
     435         ProjExchangeQ = ProjQ1; ProjQ1=TargExchangeQ; TargExchangeQ=ProjExchangeQ;
     436        } else
     437        if((TargExchangeQ != ProjQ2)||(G4UniformRand()<Same)) // Vova
     438        {
     439         ProjExchangeQ = ProjQ2; ProjQ2=TargExchangeQ; TargExchangeQ=ProjExchangeQ;
     440        } else
     441        {
     442         ProjExchangeQ = ProjQ3;  ProjQ3=TargExchangeQ; TargExchangeQ=ProjExchangeQ;
     443        }
     444//G4cout<<"TargExchangeQ "<<TargExchangeQ<<G4endl;
     445        if( Ksi < 0.333333 )
     446        {TargQ1=TargExchangeQ;}
     447        else if( (0.333333 <= Ksi) && (Ksi < 0.666667))
     448        {TargQ2=TargExchangeQ;}
     449        else
     450        {TargQ3=TargExchangeQ;}
     451/*
     452        NewProjCode = NewNucleonId(ProjQ1, ProjQ2, ProjQ3);
     453        if( (ProjQ1 == ProjQ2) && (ProjQ1 == ProjQ3) ) 
     454        {NewProjCode +=2; ProjDeltaHasCreated = true;}
     455
     456        NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3);
     457        if( (TargQ1 == TargQ2) && (TargQ1 == TargQ3) ) 
     458        {NewTargCode +=2; TargDeltaHasCreated = true;}
     459
     460        if(!ProjDeltaHasCreated && !TargDeltaHasCreated)
     461        {NewTargCode +=2; TargDeltaHasCreated = true;}
     462
     463        if(!(ProjDeltaHasCreated && TargDeltaHasCreated))
     464        {
     465         if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS >2600.))
     466         {
     467G4cout<<"2 Deltas "<<G4endl;
     468          if(!ProjDeltaHasCreated)
     469               {NewProjCode +=2; ProjDeltaHasCreated = true;}
     470          else {NewTargCode +=2; TargDeltaHasCreated = true;}
     471         }
     472        }
     473*/
     474       } // End of sampling baryon
     475
     476       NewProjCode = NewNucleonId(ProjQ1, ProjQ2, ProjQ3); // *****************************
     477
     478       if((ProjQ1==ProjQ2) && (ProjQ1==ProjQ3)) {NewProjCode +=2; ProjDeltaHasCreated=true;}
     479       else if(projectile->GetDefinition()->GetPDGiIsospin() == 3)// Projectile was Delta
     480       { if(G4UniformRand() > DeltaProbAtQuarkExchange)
     481                                                {NewProjCode +=2; ProjDeltaHasCreated=true;}
     482         else                                   {NewProjCode +=0; ProjDeltaHasCreated=false;}
     483       }
     484       else                                                       // Projectile was Nucleon
     485       {
     486        if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > DeltaMass+M0target))
     487                                                {NewProjCode +=2; ProjDeltaHasCreated=true;}
     488        else                                    {NewProjCode +=0; ProjDeltaHasCreated=false;}
     489       }
     490
     491       NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3); // *****************************
     492
     493       if((TargQ1==TargQ2) && (TargQ1==TargQ3)) {NewTargCode +=2; TargDeltaHasCreated=true;} 
     494       else if(target->GetDefinition()->GetPDGiIsospin() == 3)    // Target was Delta
     495       { if(G4UniformRand() > DeltaProbAtQuarkExchange)
     496                                                {NewTargCode +=2; TargDeltaHasCreated=true;}
     497         else                                   {NewTargCode +=0; TargDeltaHasCreated=false;}
     498       }
     499       else                                                       // Target was Nucleon
     500       {
     501        if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > M0projectile+DeltaMass))
     502                                                {NewTargCode +=2; TargDeltaHasCreated=true;}
     503        else                                    {NewTargCode +=0; TargDeltaHasCreated=false;}
     504       }         
     505
     506       if((absProjectilePDGcode == NewProjCode) && (absTargetPDGcode == NewTargCode))
     507       { // Nothing was changed! It is not right!?
     508       }
     509// Forming baryons --------------------------------------------------
     510
     511//       if( ProjDeltaHasCreated && TargDeltaHasCreated ) {StandardExcitation = true;}
     512
     513//G4cout<<"New PDG "<<NewProjCode<<" "<<NewTargCode<<G4endl;
     514//G4int Uzhi; G4cin>>Uzhi;
     515      } // End of if projectile is baryon ---------------------------
     516
     517
     518// If we assume that final state hadrons after the charge exchange will be
     519// in the ground states, we have to put ----------------------------------
     520      if(M0projectile <
     521         (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass())
     522      {
     523       M0projectile=
     524         (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass();
     525       M0projectile2 = M0projectile * M0projectile;
     526      }
     527
     528      if(M0target <
     529         (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass())
     530      {
     531       M0target=
     532         (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass();
     533       M0target2 = M0target * M0target;
     534      }
     535
     536      PZcms2=(S*S+M0projectile2*M0projectile2+M0target2*M0target2-
     537             2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2)
     538             /4./S;
     539//G4cout<<"New Masses "<<M0projectile<<" "<<M0target<<G4endl;
     540//G4cout<<"PZcms2 "<<PZcms2<<G4endl;
     541
     542      if(PZcms2 < 0) {return false;}  // It can be if energy is not sufficient for Delta
     543//----------------------------------------------------------
     544      projectile->SetDefinition(
     545                  G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode));
     546
     547      target->SetDefinition(
     548                  G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode));
     549//----------------------------------------------------------
     550      PZcms = std::sqrt(PZcms2);
     551
     552      Pprojectile.setPz( PZcms);
     553      Pprojectile.setE(std::sqrt(M0projectile2+PZcms2));
     554
     555      Ptarget.setPz(    -PZcms);
     556      Ptarget.setE(std::sqrt(M0target2+PZcms2));
     557
     558//G4cout<<"Proj "<<Pprojectile<<" "<<Pprojectile.mag()<<G4endl;
     559//G4cout<<"Targ "<<Ptarget<<" "<<Ptarget.mag()<<G4endl;
     560
     561//      if(!StandardExcitation)
     562      {
     563       Pprojectile.transform(toLab);
     564       Ptarget.transform(toLab);
     565
     566       projectile->SetTimeOfCreation(target->GetTimeOfCreation());
     567       projectile->SetPosition(target->GetPosition());
     568
     569       projectile->Set4Momentum(Pprojectile);
     570       target->Set4Momentum(Ptarget);
     571
     572       G4bool Result= theElastic->ElasticScattering (projectile,target,theParameters);
     573//G4cout<<"1 Delta result "<<Result<<"**********"<<G4endl;
     574       return Result;
     575      }
     576/*
     577else
     578      {
     579       if(M0projectile > ProjectileDiffStateMinMass)
     580       { ProjectileDiffStateMinMass +=200.*MeV; ProjectileNonDiffStateMinMass +=200.*MeV;}
     581
     582       if(M0target > TargetDiffStateMinMass)
     583       { TargetDiffStateMinMass +=200.*MeV; TargetNonDiffStateMinMass +=200.*MeV;}
     584
     585       ProbOfDiffraction         = 1.;
     586       ProbProjectileDiffraction =0.5;
     587       ProbTargetDiffraction     =0.5;
     588G4cout<<"Exc DD "<<M0projectile<<" "<<M0target<<" --------------------"<<G4endl;
     589//       After that standard FTF excitation
     590      }
     591*/
     592     }  // End of charge exchange part ------------------------------
     593
     594// ------------------------------------------------------------------
     595//     G4double ProbOfDiffraction=ProbProjectileDiffraction +
     596//                                ProbTargetDiffraction;
     597//G4cout<<ProbOfDiffraction<<" "<<ProbProjectileDiffraction<<" "<<ProbTargetDiffraction<<G4endl;
     598//G4int Uzhi; G4cin>>Uzhi;
     599     if(ProbOfDiffraction!=0.)
     600     {
     601      ProbProjectileDiffraction/=ProbOfDiffraction;
     602     }
     603     else
     604     {
     605      ProbProjectileDiffraction=0.;
     606     }
     607
     608//ProbOfDiffraction=1.; //0.5; // Vova
     609
     610     G4double ProjectileDiffStateMinMass2    = ProjectileDiffStateMinMass    *
     611                                               ProjectileDiffStateMinMass;
     612     G4double ProjectileNonDiffStateMinMass2 = ProjectileNonDiffStateMinMass *
     613                                               ProjectileNonDiffStateMinMass;
     614
     615     G4double TargetDiffStateMinMass2        = TargetDiffStateMinMass        *
     616                                               TargetDiffStateMinMass;
     617     G4double TargetNonDiffStateMinMass2     = TargetNonDiffStateMinMass     *
     618                                               TargetNonDiffStateMinMass;
     619
     620
     621     G4double Pt2;
     622     G4double ProjMassT2, ProjMassT;
     623     G4double TargMassT2, TargMassT;
     624     G4double PMinusMin, PMinusMax;
     625//   G4double PPlusMin , PPlusMax;
     626     G4double TPlusMin , TPlusMax;
     627     G4double PMinusNew, PPlusNew, TPlusNew, TMinusNew;
     628
    247629     G4LorentzVector Qmomentum;
    248630     G4double Qminus, Qplus;
    249631
    250632     G4int whilecount=0;
    251 //  Choose a process
     633//   Choose a process ---------------------------
    252634
    253635     if(G4UniformRand() < ProbOfDiffraction)
    … …  
    255637        if(G4UniformRand() < ProbProjectileDiffraction)
    256638        { //-------- projectile diffraction ---------------
    257 //G4cout<<"   Projectile diffraction"<<G4endl;
    258 //Uzhi_projectilediffraction++;
     639//G4cout<<"Proj difr"<<G4endl;
    259640         do {
    260641//             Generate pt
    … …  
    266647             {
    267648              Qmomentum=G4LorentzVector(0.,0.,0.,0.);
    268               return false;    //  Ignore this interaction
     649              target->SetStatus(2);  return false;    //  Ignore this interaction
    269650             };
    270651// --------------- Check that the interaction is possible -----------
    … …  
    278659                     2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
    279660                    /4./S;
    280 //G4cout<<" Pt2 Mpt Mtt Pz2 "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
    281 
    282 if(PZcms2 < 0 )
    283 {
    284 /*
    285 G4cout<<"whilecount "<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
    286 G4int Uzhi; G4cin>>Uzhi;
    287 */
    288 return false;
    289 };
     661//G4cout<<"PZcms2 < 0 false "<<PZcms2<<G4endl;
     662             if(PZcms2 < 0 )
     663             {
     664               target->SetStatus(2); 
     665               return false;
     666             }
    290667             maxPtSquare=PZcms2;
    291668
    … …  
    302679                     2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
    303680                    /4./S;
    304 //G4cout<<" Pt2 Mpt Mtt Pz2 "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
    305 
    306 //             if(PZcms2 < 0 ) {PZcms2=0;};
     681
    307682             if(PZcms2 < 0 ) continue;
    308683             PZcms =std::sqrt(PZcms2);
    … …  
    310685             PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
    311686             PMinusMax=SqrtS-TargMassT;
    312 //G4cout<<" SqrtS P+mim max "<<SqrtS<<" "<<PMinusMin<<" "<<PMinusMax<<G4endl;
    313687
    314688             PMinusNew=ChooseP(PMinusMin, PMinusMax);
    … …  
    328702        else
    329703        { // -------------- Target diffraction ----------------
    330 //G4cout<<"   Target difraction"<<G4endl;
    331 //Uzhi_targetdiffraction++;
     704         do {
     705//G4cout<<"Targ difr"<<G4endl;
     706//             Generate pt
     707//             if (whilecount++ >= 500 && (whilecount%100)==0)
     708//               G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping"
     709//               << ", loop count/ maxPtSquare : "
     710//               << whilecount << " / " << maxPtSquare << G4endl;
     711             if (whilecount > 1000 )
     712             {
     713              Qmomentum=G4LorentzVector(0.,0.,0.,0.);
     714              target->SetStatus(2);  return false;    //  Ignore this interaction
     715             };
     716// --------------- Check that the interaction is possible -----------
     717             ProjMassT2=M0projectile2;
     718             ProjMassT =M0projectile;
     719
     720             TargMassT2=TargetDiffStateMinMass2;
     721             TargMassT =TargetDiffStateMinMass;
     722
     723             PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
     724                     2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
     725                    /4./S;
     726
     727             if(PZcms2 < 0 )
     728             {
     729               target->SetStatus(2); 
     730               return false;
     731             }
     732             maxPtSquare=PZcms2;
     733
     734             Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
     735             Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
     736
     737             ProjMassT2=M0projectile2+Pt2;
     738             ProjMassT =std::sqrt(ProjMassT2);
     739
     740             TargMassT2=TargetDiffStateMinMass2+Pt2;
     741             TargMassT =std::sqrt(TargMassT2);
     742
     743             PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
     744                     2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
     745                    /4./S;
     746
     747             if(PZcms2 < 0 ) continue;
     748             PZcms =std::sqrt(PZcms2);
     749
     750             TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
     751             TPlusMax=SqrtS-ProjMassT;
     752
     753             TPlusNew=ChooseP(TPlusMin, TPlusMax);
     754
     755//TPlusNew=TPlusMax;
     756
     757             PPlusNew=SqrtS-TPlusNew;
     758             Qplus=PPlusNew-Pprojectile.plus();
     759             PMinusNew=ProjMassT2/PPlusNew;
     760             Qminus=PMinusNew-Pprojectile.minus();
     761
     762             Qmomentum.setPz( (Qplus-Qminus)/2 );
     763             Qmomentum.setE(  (Qplus+Qminus)/2 );
     764
     765          } while (
     766 ((Pprojectile+Qmomentum).mag2() <  M0projectile2          ) ||  //No without excitation
     767 ((Ptarget    -Qmomentum).mag2() <  TargetDiffStateMinMass2));
     768         }
     769        }
     770        else  //----------- Non-diffraction process ------------
     771        {
     772//G4cout<<"Non difr"<<G4endl;
    332773         do {
    333774//             Generate pt
    … …  
    339780             {
    340781              Qmomentum=G4LorentzVector(0.,0.,0.,0.);
    341               return false;    //  Ignore this interaction
    342              };
    343 // --------------- Check that the interaction is possible -----------
    344              ProjMassT2=M0projectile2;
    345              ProjMassT =M0projectile;
    346 
    347              TargMassT2=TargetDiffStateMinMass2;
    348              TargMassT =TargetDiffStateMinMass;
    349 
    350              PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
    351                      2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
    352                     /4./S;
    353 //G4cout<<" Pt2 Mpt Mtt Pz2 "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
    354 
    355 if(PZcms2 < 0 )
    356 {
    357 /*
    358 G4cout<<"whilecount "<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
    359 G4int Uzhi; G4cin>>Uzhi;
    360 */
    361 return false;
    362 };
    363              maxPtSquare=PZcms2;
    364 
    365              Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
    366              Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
    367 
    368              ProjMassT2=M0projectile2+Pt2;
    369              ProjMassT =std::sqrt(ProjMassT2);
    370 
    371              TargMassT2=TargetDiffStateMinMass2+Pt2;
    372              TargMassT =std::sqrt(TargMassT2);
    373 
    374              PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
    375                      2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
    376                     /4./S;
    377 /*
    378 if(PZcms2 < 0 )
    379 {
    380 G4cout<<"whilecount "<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
    381 G4int Uzhi; G4cin>>Uzhi;
    382 return false;
    383 };
    384 */
    385              if(PZcms2 < 0 ) continue;
    386              PZcms =std::sqrt(PZcms2);
    387 
    388              TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
    389              TPlusMax=SqrtS-ProjMassT;
    390 
    391 //G4cout<<" Tmin max "<<TPlusMin<<" "<<TPlusMax<<G4endl;
    392 
    393              TPlusNew=ChooseP(TPlusMin, TPlusMax);
    394 
    395 //TPlusNew=TPlusMax;
    396 //G4cout<<"T+new "<<TPlusNew<<G4endl;
    397 
    398              PPlusNew=SqrtS-TPlusNew;
    399              Qplus=PPlusNew-Pprojectile.plus();
    400              PMinusNew=ProjMassT2/PPlusNew;
    401              Qminus=PMinusNew-Pprojectile.minus();
    402 
    403              Qmomentum.setPz( (Qplus-Qminus)/2 );
    404              Qmomentum.setE(  (Qplus+Qminus)/2 );
    405 
    406           } while (
    407  ((Pprojectile+Qmomentum).mag2() <  M0projectile2          ) ||  //No without excitation
    408  ((Ptarget    -Qmomentum).mag2() <  TargetDiffStateMinMass2));
    409          }
    410         }
    411         else  //----------- Non-diffraction process ------------
    412         {
    413 //G4cout<<"   Non-difraction"<<G4endl;
    414          do {
    415 //             Generate pt
    416 //             if (whilecount++ >= 500 && (whilecount%100)==0)
    417 //               G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping"
    418 //               << ", loop count/ maxPtSquare : "
    419 //               << whilecount << " / " << maxPtSquare << G4endl;
    420              if (whilecount > 1000 )
    421              {
    422               Qmomentum=G4LorentzVector(0.,0.,0.,0.);
    423               return false;    //  Ignore this interaction
     782              target->SetStatus(2);  return false;    //  Ignore this interaction
    424783             };
    425784// --------------- Check that the interaction is possible -----------
    … …  
    433792                    2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
    434793                   /4./S;
    435 //G4cout<<" Pt2 Mpt Mtt Pz2 "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
    436 
    437 if(PZcms2 < 0 )
    438 {
    439 /*
    440 G4cout<<"whilecount "<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
    441 G4int Uzhi; G4cin>>Uzhi;
    442 */
    443 return false;
    444 };
     794
     795             if(PZcms2 < 0 )
     796             {
     797               target->SetStatus(2); 
     798               return false;
     799             }
    445800             maxPtSquare=PZcms2;
    446801             Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
    … …  
    456811                    2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
    457812                   /4./S;
    458 /*
    459 G4cout<<"ProjectileNonDiffStateMinMass2 "<<ProjectileNonDiffStateMinMass2<<G4endl;
    460 G4cout<<"TargetNonDiffStateMinMass2     "<<TargetNonDiffStateMinMass2<<G4endl;
    461 G4cout<<"Mt "<<ProjMassT<<" "<<TargMassT<<" "<<Pt2<<" "<<PZcms2<<G4endl<<G4endl;
    462 */
    463 //             if(PZcms2 < 0 ) {PZcms2=0;};
     813
    464814             if(PZcms2 < 0 ) continue;
    465815             PZcms =std::sqrt(PZcms2);
    … …  
    469819
    470820             PMinusNew=ChooseP(PMinusMin, PMinusMax);
    471 // PMinusNew=1./sqrt(1./PMinusMin-G4UniformRand()*(1./PMinusMin-1./PMinusMax));
    472 
    473 //G4cout<<"Proj "<<PMinusMin<<" "<<PMinusMax<<" "<<PMinusNew<<G4endl;
    474 
    475 //PMinusNew=PMinusMax; //+++++++++++++++++++++++++++++++++++ Vova
    476821
    477822             Qminus=PMinusNew-Pprojectile.minus();
    478823
    479824             TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
    480 //             TPlusMax=SqrtS-PMinusNew;                      // Vova
    481              TPlusMax=SqrtS-ProjMassT;                        // Vova
     825//           TPlusMax=SqrtS-PMinusNew;                     
     826             TPlusMax=SqrtS-ProjMassT;     
    482827
    483828             TPlusNew=ChooseP(TPlusMin, TPlusMax);
    484 
    485 //G4cout<<"Targ "<<TPlusMin<<" "<<TPlusMax<<" "<<TPlusNew<<G4endl;
    486 //G4cout<<PMinusNew<<" "<<TPlusNew<<G4endl;
    487829
    488830             Qplus=-(TPlusNew-Ptarget.plus());
    … …  
    490832             Qmomentum.setPz( (Qplus-Qminus)/2 );
    491833             Qmomentum.setE(  (Qplus+Qminus)/2 );
    492 /*
    493 G4cout << "Qplus / Qminus " << Qplus << " / " << Qminus<<G4endl;
    494 G4cout << "pt2" << pt2 << G4endl;
    495 G4cout << "Qmomentum " << Qmomentum << G4endl;
    496 G4cout << " Masses (P/T) : " << (Pprojectile+Qmomentum).mag() <<
    497            " / " << (Ptarget-Qmomentum).mag() << G4endl;   // mag()
    498 G4cout<<"Mprojectile "<<std::sqrt(M0projectile2)<<G4endl;
    499 G4cout<<"Mtarget     "<<std::sqrt(M0target2    )<<G4endl;
    500 G4cout<<"ProjectileDiffStateMinMass "<<std::sqrt(ProjectileDiffStateMinMass2)<<G4endl;
    501 G4cout<<"TargetDiffStateMinMass "<<std::sqrt(TargetDiffStateMinMass2)<<G4endl;
    502 */
     834
    503835       } while (
    504836 ((Pprojectile+Qmomentum).mag2() <  ProjectileNonDiffStateMinMass2) || //No double Diffraction
    … …  
    506838         }
    507839
    508 //G4int Uzhiinp; G4cin>>Uzhiinp;    // Vova
    509 
    510840           Pprojectile += Qmomentum;
    511841           Ptarget     -= Qmomentum;
    512 /*
    513 G4cout << "Pprojectile with Q : " << Pprojectile << G4endl;
    514 G4cout << "Ptarget with Q : " << Ptarget << G4endl;
    515 G4cout << "Target        mass  " <<  Ptarget.mag() << G4endl;
    516 G4cout << "Projectile mass  " <<  Pprojectile.mag() << G4endl;
    517 //
    518 //G4cout << "Projectile back: " << toLab * Pprojectile << G4endl;
    519 //G4cout << "Target back: " << toLab * Ptarget << G4endl;
    520 */
    521 //-------------- Flip if projectale moves in backward direction ------------
    522 //G4bool Flip=Pprojectile.pz()< 0.;
    523 
     842
     843//G4cout<<"Masses "<<Pprojectile.mag()<<" "<<Ptarget.mag()<<G4endl;
     844//G4int Uzhi; G4cin>>Uzhi;
    524845
    525846// Transform back and update SplitableHadron Participant.
    … …  
    527848           Ptarget.transform(toLab);
    528849
    529 //G4cout << "Pprojectile with Q M: " << Pprojectile<<" "<<  Pprojectile.mag() << G4endl;
    530 //G4cout << "Ptarget     with Q M: " << Ptarget    <<" "<<  Ptarget.mag()     << G4endl;
    531 //G4cout << "Target      mass  " <<  Ptarget.mag() << G4endl;
    532 //G4cout << "Projectile mass  " <<  Pprojectile.mag() << G4endl;
    533 
    534 /*
    535 if(!Flip){
     850// Calculation of the creation time ---------------------
     851      projectile->SetTimeOfCreation(target->GetTimeOfCreation());
     852      projectile->SetPosition(target->GetPosition());
     853// Creation time and position of target nucleon were determined at
     854// ReggeonCascade() of G4FTFModel
     855// ------------------------------------------------------
     856
    536857           projectile->Set4Momentum(Pprojectile);
    537858           target->Set4Momentum(Ptarget);
    538          }
    539 else     {
    540            G4ParticleDefinition * t_Definition=projectile->GetDefinition();
    541            projectile->SetDefinition(target->GetDefinition());
    542            projectile->Set4Momentum(Ptarget);
    543            target->SetDefinition(t_Definition);
    544            target->Set4Momentum(Pprojectile);
    545           }
    546 */
    547 //
    548 /*
    549 if(G4UniformRand() < 1.) {
    550            G4ParticleDefinition * t_Definition=projectile->GetDefinition();
    551            projectile->SetDefinition(target->GetDefinition());
    552            target->SetDefinition(t_Definition);
    553 }
    554 */ // For flip, for HARP
    555 
    556            G4double ZcoordinateOfCurrentInteraction = target->GetPosition().z();
    557 // It is assumed that nucleon z-coordinates are ordered on increasing -----------
    558 
    559            G4double betta_z=projectile->Get4Momentum().pz()/projectile->Get4Momentum().e();
    560 
    561            G4double ZcoordinateOfPreviousCollision=projectile->GetPosition().z();
    562            if(projectile->GetSoftCollisionCount()==0) {
    563               projectile->SetTimeOfCreation(0.);
    564               target->SetTimeOfCreation(0.);
    565               ZcoordinateOfPreviousCollision=ZcoordinateOfCurrentInteraction;
    566            }
    567            
    568            G4ThreeVector thePosition(projectile->GetPosition().x(),
    569                                      projectile->GetPosition().y(),
    570                                      ZcoordinateOfCurrentInteraction);
    571            projectile->SetPosition(thePosition);
    572 
    573            G4double TimeOfPreviousCollision=projectile->GetTimeOfCreation();
    574            G4double TimeOfCurrentCollision=TimeOfPreviousCollision+
    575                     (ZcoordinateOfCurrentInteraction-ZcoordinateOfPreviousCollision)/betta_z;
    576 
    577            projectile->SetTimeOfCreation(TimeOfCurrentCollision);
    578            target->SetTimeOfCreation(TimeOfCurrentCollision);
    579 
    580            projectile->Set4Momentum(Pprojectile);
    581            target->Set4Momentum(Ptarget);
    582859
    583860           projectile->IncrementCollisionCount(1);
    584861           target->IncrementCollisionCount(1);
    585862
    586 //
    587 //G4cout<<"Out of Excitation --------------------"<<G4endl;
    588 //G4int Uzhiinp; G4cin>>Uzhiinp;    // Vova
    589 
    590863           return true;
    591864}
    592865
    593866// ---------------------------------------------------------------------
    594 G4ExcitedString * G4DiffractiveExcitation::
    595            String(G4VSplitableHadron * hadron, G4bool isProjectile) const
     867//G4ExcitedString * G4DiffractiveExcitation::
     868//           String(G4VSplitableHadron * hadron, G4bool isProjectile) const
     869void G4DiffractiveExcitation::CreateStrings(G4VSplitableHadron * hadron,
     870                                            G4bool isProjectile,
     871                                            G4ExcitedString * &FirstString,
     872                                            G4ExcitedString * &SecondString,
     873                                            G4FTFParameters *theParameters) const
    596874{
    597 
    598 //G4cout<<"G4DiffractiveExcitation::String isProj"<<isProjectile<<G4endl;
    599 
    600875        hadron->SplitUp();
    601876        G4Parton *start= hadron->GetNextParton();
    602877        if ( start==NULL)
    603878        { G4cout << " G4FTFModel::String() Error:No start parton found"<< G4endl;
    604           return NULL;
     879          FirstString=0; SecondString=0;
     880          return;
    605881        }
    606882        G4Parton *end  = hadron->GetNextParton();
    607883        if ( end==NULL)
    608884        { G4cout << " G4FTFModel::String() Error:No end parton found"<< G4endl;
    609           return NULL;
     885          FirstString=0; SecondString=0;
     886          return;
    610887        }
    611 
    612         G4ExcitedString * string;
    613         if ( isProjectile )
    614         {
    615                 string= new G4ExcitedString(end,start, +1);
    616         } else {
    617                 string= new G4ExcitedString(start,end, -1);
    618         }
    619 // Uzhi
    620 //G4cout<<"G4ExcitedString * G4DiffractiveExcitation::String"<<G4endl;
    621 //G4cout<<hadron->GetTimeOfCreation()<<" "<<hadron->GetPosition()/fermi<<G4endl;
    622 
    623         string->SetTimeOfCreation(hadron->GetTimeOfCreation());
    624         string->SetPosition(hadron->GetPosition());
     888        G4LorentzVector Phadron=hadron->Get4Momentum();
     889/*
     890G4cout<<"Create strings had "<<hadron->GetDefinition()->GetParticleName()<<" "<<Phadron<<" "<<Phadron.mag()<<G4endl;
     891G4cout<<"isProjectile "<<isProjectile<<G4endl;
     892G4cout<<"start Q "<<start->GetDefinition()->GetPDGEncoding()<<G4endl;
     893G4cout<<"end   Q "<<end->GetDefinition()->GetPDGEncoding()<<G4endl;
     894*/
     895        G4LorentzVector Pstart(0.,0.,0.,0.);
     896        G4LorentzVector Pend(0.,0.,0.,0.);
     897        G4LorentzVector Pkink(0.,0.,0.,0.);
     898        G4LorentzVector PkinkQ1(0.,0.,0.,0.);
     899        G4LorentzVector PkinkQ2(0.,0.,0.,0.);
     900
     901        G4int PDGcode_startQ = std::abs(start->GetDefinition()->GetPDGEncoding());
     902        G4int PDGcode_endQ   = std::abs(  end->GetDefinition()->GetPDGEncoding());
     903
     904//--------------------------------------------------------------------------------       
     905        G4double Wmin(0.);
     906        if(isProjectile)
     907        {
     908          Wmin=theParameters->GetProjMinDiffMass();
     909        } else
     910        {
     911          Wmin=theParameters->GetTarMinDiffMass();
     912        } // end of if(isProjectile)
     913
     914        G4double W = hadron->Get4Momentum().mag();
     915        G4double W2=W*W;
     916//G4cout<<"W Wmin "<<W<<" "<<Wmin<<G4endl;
     917        G4double Pt(0.), x1(0.), x2(0.), x3(0.);
     918        G4bool Kink=false;
     919
     920        if(W > Wmin)
     921        {                                        // Kink is possible
     922          G4double Pt2kink=theParameters->GetPt2Kink();
     923          Pt = std::sqrt(Pt2kink*(std::pow(W2/16./Pt2kink+1.,G4UniformRand()) - 1.));
     924// Pt=0.;
     925//G4cout<<"Pt2kink Pt Pt2 "<<Pt2kink<<" "<<Pt<<" "<<Pt*Pt<<G4endl;
     926          if(Pt > 500.*MeV)
     927          {
     928             G4double Ymax = std::log(W/2./Pt + std::sqrt(W2/4./Pt/Pt - 1.));
     929             G4double Y=Ymax*(1.- 2.*G4UniformRand());
     930//G4cout<<"Ymax Y "<<Ymax<<" "<<Y<<G4endl;
     931             x1=1.-Pt/W*std::exp( Y);
     932             x3=1.-Pt/W*std::exp(-Y);
     933             x2=2.-x1-x3;
     934//G4cout<<"X1 X2 X3 "<<x1<<" "<<x2<<" "<<x3<<G4endl;
     935             G4double Mass_startQ = 650.*MeV;
     936             if(PDGcode_startQ <  3) Mass_startQ =  325.*MeV;
     937             if(PDGcode_startQ == 3) Mass_startQ =  500.*MeV;
     938             if(PDGcode_startQ == 4) Mass_startQ = 1600.*MeV;
     939
     940             G4double Mass_endQ = 650.*MeV;
     941             if(PDGcode_endQ <  3) Mass_endQ =  325.*MeV;
     942             if(PDGcode_endQ == 3) Mass_endQ =  500.*MeV;
     943             if(PDGcode_endQ == 4) Mass_endQ = 1600.*MeV;
     944
     945             G4double P2_1=W2*x1*x1/4.-Mass_endQ  *Mass_endQ;
     946             G4double P2_3=W2*x3*x3/4.-Mass_startQ*Mass_startQ;
     947     
     948             G4double P2_2=sqr((2.-x1-x3)*W/2.);
     949//G4cout<<"P2_1 P2_2 P2_3 "<<P2_1<<" "<<P2_2<<" "<<P2_3<<G4endl;
     950             if((P2_1 < 0.) || (P2_3 <0.))
     951             { Kink=false;}
     952             else
     953             {
     954               G4double P_1=std::sqrt(P2_1);
     955               G4double P_2=std::sqrt(P2_2);
     956               G4double P_3=std::sqrt(P2_3);
     957//G4cout<<"P_1 P_2 P_3 "<<P_1<<" "<<P_2<<" "<<P_3<<G4endl;
     958               G4double CosT12=(P2_3-P2_1-P2_2)/(2.*P_1*P_2);
     959               G4double CosT13=(P2_2-P2_1-P2_3)/(2.*P_1*P_3);
     960//G4cout<<"CosT12 CosT13 "<<CosT12<<" "<<CosT13<<" "<<std::acos(CosT12)*180./3.14159<<" "<<std::acos(CosT13)*180./3.14159<<G4endl;
     961               if((std::abs(CosT12) >1.) || (std::abs(CosT13) > 1.))
     962               { Kink=false;}
     963               else
     964               {
     965                 Kink=true;
     966                 Pstart.setPx(-Pt); Pstart.setPy(0.); Pstart.setPz(P_3*CosT13);
     967                 Pend.setPx(   0.); Pend.setPy(  0.); Pend.setPz(          P_1);
     968                 Pkink.setPx(  Pt); Pkink.setPy( 0.); Pkink.setPz(  P_2*CosT12);
     969                 Pstart.setE(x3*W/2.);               
     970                 Pkink.setE(Pkink.vect().mag());
     971                 Pend.setE(x1*W/2.);
     972
     973                 G4double XkQ=GetQuarkFractionOfKink(0.,1.);
     974                 if(Pkink.getZ() > 0.)
     975                 {
     976                   if(XkQ > 0.5) {PkinkQ1=XkQ*Pkink;} else {PkinkQ1=(1.-XkQ)*Pkink;}
     977                 } else {
     978                   if(XkQ > 0.5) {PkinkQ1=(1.-XkQ)*Pkink;} else {PkinkQ1=XkQ*Pkink;}
     979                 }
     980
     981                 PkinkQ2=Pkink - PkinkQ1;
     982//------------------------- Minimizing Pt1^2+Pt3^2 ------------------------------
     983/*
     984G4cout<<"Pstart "<<Pstart<<G4endl;
     985G4cout<<"Pkink  "<<Pkink <<G4endl;
     986G4cout<<"Pkink1 "<<PkinkQ1<<G4endl;
     987G4cout<<"Pkink2 "<<PkinkQ2<<G4endl;
     988G4cout<<"Pend   "<<Pend  <<G4endl;
     989*/
     990                 G4double Cos2Psi=(sqr(x1) -sqr(x3)+2.*sqr(x3*CosT13))/
     991                          std::sqrt(sqr(sqr(x1)-sqr(x3)) + sqr(2.*x1*x3*CosT13));
     992                 G4double Psi=std::acos(Cos2Psi);
     993//G4cout<<"Psi "<<Psi<<" "<<Psi*180./twopi<<G4endl;
     994
     995G4LorentzRotation Rotate;
     996if(isProjectile) {Rotate.rotateY(Psi);}
     997else             {Rotate.rotateY(pi-Psi);}                   
     998Rotate.rotateZ(twopi*G4UniformRand());
     999
     1000//G4cout<<"Rotate"<<G4endl;
     1001
     1002Pstart*=Rotate;
     1003Pkink*=Rotate;
     1004PkinkQ1*=Rotate;
     1005PkinkQ2*=Rotate;
     1006Pend*=Rotate;
     1007/*
     1008G4cout<<"Pstart "<<Pstart<<G4endl;
     1009G4cout<<"Pkink1 "<<PkinkQ1 <<G4endl;
     1010G4cout<<"Pkink2 "<<PkinkQ2 <<G4endl;
     1011G4cout<<"Pend   "<<Pend  <<G4endl;
     1012*/
     1013               }
     1014             }      // end of if((P2_1 < 0.) || (P2_3 <0.))
     1015          }         // end of if(Pt > 500.*MeV)
     1016        }           // end of if(W > Wmin) Check for a kink
     1017
     1018//--------------------------------------------------------------------------------
     1019
     1020//G4cout<<"Kink "<<Kink<<G4endl;
     1021
     1022        if(Kink)
     1023        {                                        // Kink is possible
     1024          std::vector<G4double> QuarkProbabilitiesAtGluonSplitUp =
     1025              theParameters->GetQuarkProbabilitiesAtGluonSplitUp();
     1026
     1027          G4int QuarkInGluon(1); G4double Ksi=G4UniformRand();
     1028          for(unsigned int Iq=0; Iq <3; Iq++)
     1029          {
     1030//G4cout<<Iq<<" "<<QuarkProbabilitiesAtGluonSplitUp[Iq]<<G4endl;
     1031if(Ksi > QuarkProbabilitiesAtGluonSplitUp[Iq]) QuarkInGluon++;}
     1032
     1033//G4cout<<"Gquark "<<QuarkInGluon<<G4endl;
     1034          G4Parton * Gquark = new G4Parton(QuarkInGluon);
     1035          G4Parton * Ganti_quark = new G4Parton(-QuarkInGluon);
     1036/*
     1037G4cout<<Gquark->GetDefinition()->GetParticleName()<<" "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->GetDefinition()->GetPDGMass()<<G4endl;
     1038G4cout<<Ganti_quark->GetDefinition()->GetParticleName()<<" "<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->GetDefinition()->GetPDGMass()<<G4endl;
     1039*/
     1040
     1041//G4int Uzhi; G4cin>>Uzhi;
     1042
     1043//-------------------------------------------------------------------------------
     1044/*
     1045DefineMomentumInZ(G4double aLightConeMomentum, G4bool aDirection);
     1046      void Set4Momentum(const G4LorentzVector & aMomentum);
     1047      G4int PDGencoding;
     1048      G4ParticleDefinition * theDefinition;
     1049      G4LorentzVector theMomentum;
     1050      G4ThreeVector   thePosition;
     1051     
     1052      G4int theColour;
     1053      G4double theIsoSpinZ;
     1054      G4double theSpinZ;
     1055     
     1056      G4double theX;
     1057*/
     1058//-------------------------------------------------------------------------------
     1059
     1060//G4cout<<"Phadron "<<Phadron<<" mass "<<Phadron.mag()<<G4endl;
     1061          G4LorentzRotation toCMS(-1*Phadron.boostVector());
     1062//G4cout<<"To lab"<<G4endl;
     1063          G4LorentzRotation toLab(toCMS.inverse());
     1064
     1065          Pstart.transform(toLab);  start->Set4Momentum(Pstart);
     1066          PkinkQ1.transform(toLab);
     1067          PkinkQ2.transform(toLab);
     1068          Pend.transform(toLab);    end->Set4Momentum(Pend);
     1069/*
     1070G4cout<<"Pstart "<<start->GetDefinition()->GetPDGEncoding()<<Pstart<<G4endl;
     1071G4cout<<"Pkink1 "<<PkinkQ1 <<G4endl;
     1072G4cout<<"Pkink2 "<<PkinkQ2 <<G4endl;
     1073G4cout<<"Pend   "<<end->GetDefinition()->GetPDGEncoding()<<Pend  <<G4endl;
     1074*/
     1075// !!!    G4ExcitedString * FirstString(0); G4ExcitedString * SecondString(0);
     1076          G4int absPDGcode=std::abs(hadron->GetDefinition()->GetPDGEncoding());
     1077
     1078//G4cout<<"absPDGcode "<<absPDGcode<<G4endl;
     1079
     1080          if(absPDGcode < 1000)
     1081          {                                // meson
     1082            if ( isProjectile )
     1083            {                              // Projectile
     1084              if(end->GetDefinition()->GetPDGEncoding() > 0 )  // A quark on the end
     1085              {                            // Quark on the end
     1086                FirstString = new G4ExcitedString(end   ,Ganti_quark, +1);
     1087                SecondString= new G4ExcitedString(Gquark,start      ,+1);
     1088                Ganti_quark->Set4Momentum(PkinkQ1);
     1089                Gquark->Set4Momentum(PkinkQ2);
     1090/*
     1091G4cout<<" Proj Meson end Q"<<G4endl;
     1092G4cout<<"First string  ============ "<<G4endl;
     1093G4cout<<"end  Q "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl;
     1094G4cout<<"G antiQ"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl;
     1095G4cout<<"Sum P  "<<(Ganti_quark->Get4Momentum()+end->Get4Momentum())<<G4endl;
     1096G4cout<<"Secondstring   ============ "<<G4endl;
     1097G4cout<<"G    Q "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl;
     1098G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl;
     1099
     1100G4cout<<"Sum P  "<<(Gquark->Get4Momentum()+start->Get4Momentum())<<G4endl;
     1101*/
     1102              } else
     1103              {                            // Anti_Quark on the end
     1104                FirstString = new G4ExcitedString(end        ,Gquark, +1);
     1105                SecondString= new G4ExcitedString(Ganti_quark,start ,+1);
     1106                Gquark->Set4Momentum(PkinkQ1);
     1107                Ganti_quark->Set4Momentum(PkinkQ2);
     1108/*
     1109G4cout<<" Proj Meson end Qbar"<<G4endl;
     1110G4cout<<"First string  ============ "<<G4endl;
     1111G4cout<<"end  Q "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl;
     1112G4cout<<"G     Q"<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl;
     1113G4cout<<"Sum P  "<<(Gquark->Get4Momentum()+end->Get4Momentum())<<G4endl;
     1114G4cout<<"Secondstring   ============ "<<G4endl;
     1115G4cout<<"G antQ "<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl;
     1116G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl;
     1117G4cout<<"Sum P  "<<(Ganti_quark->Get4Momentum()+start->Get4Momentum())<<G4endl;
     1118*/
     1119              }   // end of if(end->GetPDGcode() > 0)
     1120            } else {                      // Target
     1121              if(end->GetDefinition()->GetPDGEncoding() > 0 )  // A quark on the end
     1122              {                           // Quark on the end
     1123                FirstString = new G4ExcitedString(Ganti_quark,end   ,-1);
     1124                SecondString= new G4ExcitedString(start      ,Gquark,-1);
     1125                Ganti_quark->Set4Momentum(PkinkQ2);
     1126                Gquark->Set4Momentum(PkinkQ1);
     1127/*
     1128G4cout<<" Targ Meson end Q"<<G4endl;
     1129G4cout<<"First string   ============ "<<G4endl;
     1130G4cout<<"G antiQ"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl;
     1131G4cout<<"end  Q "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl;
     1132G4cout<<"Sum P  "<<(Ganti_quark->Get4Momentum()+end->Get4Momentum())<<G4endl;
     1133G4cout<<"Secondstring   ============ "<<G4endl;
     1134G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl;
     1135G4cout<<"G    Q "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl;
     1136G4cout<<"Sum P  "<<(Gquark->Get4Momentum()+start->Get4Momentum())<<G4endl;
     1137*/
     1138              } else
     1139              {                            // Anti_Quark on the end
     1140                FirstString = new G4ExcitedString(Gquark,end        ,-1);
     1141                SecondString= new G4ExcitedString(start ,Ganti_quark,-1);
     1142                Gquark->Set4Momentum(PkinkQ2);
     1143                Ganti_quark->Set4Momentum(PkinkQ1);
     1144/*
     1145G4cout<<" Targ Meson end Qbar"<<G4endl;
     1146G4cout<<"First string   ============ "<<G4endl;
     1147G4cout<<"G     Q"<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl;
     1148G4cout<<"end  Q "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl;
     1149G4cout<<"Sum P  "<<(Gquark->Get4Momentum()+end->Get4Momentum())<<G4endl;
     1150G4cout<<"Secondstring   ============ "<<G4endl;
     1151G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl;
     1152G4cout<<"G antQ "<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl;
     1153G4cout<<"Sum P  "<<(Ganti_quark->Get4Momentum()+start->Get4Momentum())<<G4endl;
     1154*/
     1155              }   // end of if(end->GetPDGcode() > 0)
     1156            }     // end of if ( isProjectile )
     1157          } else  // if(absPDGCode < 1000)
     1158          {                             // Baryon/AntiBaryon
     1159            if ( isProjectile )
     1160            {                              // Projectile
     1161              if((end->GetDefinition()->GetParticleType() == "diquarks") &&
     1162                 (end->GetDefinition()->GetPDGEncoding() > 0           )   )
     1163              {                            // DiQuark on the end
     1164                FirstString = new G4ExcitedString(end        ,Gquark, +1);
     1165                SecondString= new G4ExcitedString(Ganti_quark,start ,+1);
     1166                Gquark->Set4Momentum(PkinkQ1);
     1167                Ganti_quark->Set4Momentum(PkinkQ2);
     1168/*
     1169G4cout<<" Proj baryon end QQ"<<G4endl;
     1170G4cout<<"First string   ============ "<<G4endl;
     1171G4cout<<"end OQ "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl;
     1172G4cout<<"G    Q"<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl;
     1173G4cout<<"Sum P  "<<(Gquark->Get4Momentum()+end->Get4Momentum())<<G4endl;
     1174G4cout<<"Secondstring   ============ "<<G4endl;
     1175G4cout<<"G  Qbar"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl;
     1176G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl;
     1177G4cout<<"Sum P  "<<(Ganti_quark->Get4Momentum()+start->Get4Momentum())<<G4endl;
     1178*/
     1179              } else
     1180              {                            // Anti_DiQuark on the end or quark
     1181                FirstString = new G4ExcitedString(end   ,Ganti_quark, +1);
     1182                SecondString= new G4ExcitedString(Gquark,start      ,+1);
     1183                Ganti_quark->Set4Momentum(PkinkQ1);
     1184                Gquark->Set4Momentum(PkinkQ2);
     1185/*
     1186G4cout<<" Proj baryon end Q"<<G4endl;
     1187G4cout<<"First string   ============ "<<G4endl;
     1188G4cout<<"end OQ "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl;
     1189G4cout<<"G antQ"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl;
     1190G4cout<<"Sum P  "<<(Ganti_quark->Get4Momentum()+end->Get4Momentum())<<G4endl;
     1191G4cout<<"Secondstring   ============ "<<G4endl;
     1192G4cout<<"G  Q   "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl;
     1193G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl;
     1194G4cout<<"Sum P  "<<(Gquark->Get4Momentum()+start->Get4Momentum())<<G4endl;
     1195*/
     1196              }   // end of if(end->GetPDGcode() > 0)
     1197            } else {                      // Target
     1198
     1199              if((end->GetDefinition()->GetParticleType() == "diquarks") &&
     1200                 (end->GetDefinition()->GetPDGEncoding() > 0           )   )
     1201              {                            // DiQuark on the end
     1202                FirstString = new G4ExcitedString(Gquark,end        ,-1);
     1203
     1204                SecondString= new G4ExcitedString(start ,Ganti_quark,-1);
     1205                Gquark->Set4Momentum(PkinkQ1);
     1206                Ganti_quark->Set4Momentum(PkinkQ2);
     1207/*
     1208G4cout<<" Targ baryon end QQ"<<G4endl;
     1209G4cout<<"First string   ============ "<<G4endl;
     1210G4cout<<"G  Q   "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl;
     1211G4cout<<"end OQ "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl;
     1212G4cout<<"Sum P  "<<(Gquark->Get4Momentum()+end->Get4Momentum())<<G4endl;
     1213G4cout<<"Secondstring   ============ "<<G4endl;
     1214G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl;
     1215G4cout<<"G  Qbar"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl;
     1216G4cout<<"Sum P  "<<(Ganti_quark->Get4Momentum()+start->Get4Momentum())<<G4endl;
     1217*/
     1218              } else
     1219              {                            // Anti_DiQuark on the end or Q
     1220                FirstString = new G4ExcitedString(Ganti_quark,end   ,-1);
     1221                SecondString= new G4ExcitedString(start      ,Gquark,-1);
     1222                Gquark->Set4Momentum(PkinkQ2);
     1223                Ganti_quark->Set4Momentum(PkinkQ1);
     1224/*
     1225G4cout<<" Targ baryon end Q"<<G4endl;
     1226G4cout<<"First string   ============ "<<G4endl;
     1227G4cout<<"G  Qbar"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl;
     1228G4cout<<"end O "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl;
     1229G4cout<<"Sum P  "<<(Ganti_quark->Get4Momentum()+end->Get4Momentum())<<G4endl;
     1230G4cout<<"Secondstring   ============ "<<G4endl;
     1231G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl;
     1232G4cout<<"G  Q   "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl;
     1233G4cout<<"Sum P  "<<(Gquark->Get4Momentum()+start->Get4Momentum())<<G4endl;
     1234*/
     1235              }   // end of if(end->GetPDGcode() > 0)
     1236            }     // end of if ( isProjectile )
     1237          }  // end of if(absPDGcode < 1000)
     1238
     1239          FirstString->SetTimeOfCreation(hadron->GetTimeOfCreation());
     1240          FirstString->SetPosition(hadron->GetPosition());
     1241
     1242          SecondString->SetTimeOfCreation(hadron->GetTimeOfCreation());
     1243          SecondString->SetPosition(hadron->GetPosition());
     1244
     1245// ------------------------------------------------------------------------- 
     1246        } else                                   // End of kink is possible
     1247        {                                        // Kink is impossible
     1248          if ( isProjectile )
     1249          {
     1250                FirstString= new G4ExcitedString(end,start, +1);
     1251          } else {
     1252                FirstString= new G4ExcitedString(start,end, -1);
     1253          }
     1254
     1255          SecondString=0;
     1256
     1257          FirstString->SetTimeOfCreation(hadron->GetTimeOfCreation());
     1258          FirstString->SetPosition(hadron->GetPosition());
    6251259
    6261260// momenta of string ends
    6271261//
    628         G4double Momentum=hadron->Get4Momentum().vect().mag();
    629         G4double  Plus=hadron->Get4Momentum().e() + Momentum;
    630         G4double Minus=hadron->Get4Momentum().e() - Momentum;
    631 
    632         G4ThreeVector tmp;
    633         if(Momentum > 0.)
    634         {
    635          tmp.set(hadron->Get4Momentum().px(),
    636                  hadron->Get4Momentum().py(),
    637                  hadron->Get4Momentum().pz());
    638          tmp/=Momentum;
    639          }
    640          else
    641          {
    642           tmp.set(0.,0.,1.);
    643          };
    644 
    645         G4LorentzVector Pstart(tmp,0.);
    646         G4LorentzVector   Pend(tmp,0.);
    647 
    648         if(isProjectile)
    649         {
    650          Pstart*=(-1.)*Minus/2.;
    651          Pend  *=(+1.)*Plus /2.;
    652          }
    653         else
    654          {
    655           Pstart*=(+1.)*Plus/2.;
    656           Pend  *=(-1.)*Minus/2.;
    657          };
    658 
    659          Momentum=-Pstart.mag();
    660          Pstart.setT(Momentum);  // It is assumed that quark has m=0.
    661 
    662          Momentum=-Pend.mag();
    663          Pend.setT(Momentum);    // It is assumed that di-quark has m=0.
    664 //
    665 /* Uzhi
    666         G4double ptSquared= hadron->Get4Momentum().perp2();
    667         G4double transverseMassSquared= hadron->Get4Momentum().plus()
    668                                     *   hadron->Get4Momentum().minus();
    669 
    670 
    671         G4double maxAvailMomentumSquared=
    672                  sqr( std::sqrt(transverseMassSquared) - std::sqrt(ptSquared) );
    673 
    674         G4double widthOfPtSquare = 0.25*GeV*GeV;       // Uzhi 11.07 <Pt^2>=0.25 ??????????????????
    675         G4ThreeVector pt=GaussianPt(widthOfPtSquare,maxAvailMomentumSquared);
    676 
    677         G4LorentzVector Pstart(G4LorentzVector(pt,0.));
    678         G4LorentzVector Pend;
    679         Pend.setPx(hadron->Get4Momentum().px() - pt.x());
    680         Pend.setPy(hadron->Get4Momentum().py() - pt.y());
    681 
    682         G4double tm1=hadron->Get4Momentum().minus() +
    683           ( Pend.perp2()-Pstart.perp2() ) / hadron->Get4Momentum().plus();
    684 
    685         G4double tm2= std::sqrt( std::max(0., sqr(tm1) -
    686              4. * Pend.perp2() * hadron->Get4Momentum().minus()
    687               /  hadron->Get4Momentum().plus() ));
    688 
    689         G4int Sign= isProjectile ? -1 : 1;
    690 
    691         G4double endMinus  = 0.5 * (tm1 + Sign*tm2);
    692         G4double startMinus= hadron->Get4Momentum().minus() - endMinus;
    693 
    694         G4double startPlus= Pstart.perp2() /  startMinus;
    695         G4double endPlus  = hadron->Get4Momentum().plus() - startPlus;
    696 
    697         Pstart.setPz(0.5*(startPlus - startMinus));
    698         Pstart.setE(0.5*(startPlus + startMinus));
    699 
    700         Pend.setPz(0.5*(endPlus - endMinus));
    701         Pend.setE(0.5*(endPlus + endMinus));
    702 */ // Uzhi
    703         start->Set4Momentum(Pstart);
    704         end->Set4Momentum(Pend);
    705 /*
    706 G4cout<<"G4DiffractiveExcitation::String hadro"<<hadron->Get4Momentum()<<" "<<hadron->Get4Momentum().mag2()<<G4endl;
    707 
    708 G4cout<<"G4DiffractiveExcitation::String start"<<start->Get4Momentum()<<" "<<start->GetPDGcode()<<G4endl;
    709 
    710 G4cout<<"G4DiffractiveExcitation::String end  "<<  end->Get4Momentum()<<" "<<  end->GetPDGcode()<<G4endl;
    711 G4int Uzhi; G4cin>>Uzhi;
    712 */
     1262          G4double Momentum=hadron->Get4Momentum().vect().mag();
     1263          G4double  Plus=hadron->Get4Momentum().e() + Momentum;
     1264          G4double Minus=hadron->Get4Momentum().e() - Momentum;
     1265
     1266          G4ThreeVector tmp;
     1267          if(Momentum > 0.)
     1268          {
     1269           tmp.set(hadron->Get4Momentum().px(),
     1270                   hadron->Get4Momentum().py(),
     1271                   hadron->Get4Momentum().pz());
     1272           tmp/=Momentum;
     1273          }
     1274          else
     1275          {
     1276           tmp.set(0.,0.,1.);
     1277          }
     1278
     1279          G4LorentzVector Pstart(tmp,0.);
     1280          G4LorentzVector   Pend(tmp,0.);
     1281
     1282          if(isProjectile)
     1283          {
     1284           Pstart*=(-1.)*Minus/2.;
     1285           Pend  *=(+1.)*Plus /2.;
     1286          }
     1287          else
     1288          {
     1289           Pstart*=(+1.)*Plus/2.;
     1290           Pend  *=(-1.)*Minus/2.;
     1291          }
     1292
     1293          Momentum=-Pstart.mag();
     1294          Pstart.setT(Momentum);  // It is assumed that quark has m=0.
     1295
     1296          Momentum=-Pend.mag();
     1297          Pend.setT(Momentum);    // It is assumed that di-quark has m=0.
     1298
     1299          start->Set4Momentum(Pstart);
     1300          end->Set4Momentum(Pend);
     1301          SecondString=0;
     1302        }            // End of kink is impossible
     1303
    7131304#ifdef G4_FTFDEBUG
    714         G4cout << " generated string flavors          "
    715                << start->GetPDGcode() << " / "
    716                << end->GetPDGcode()                    << G4endl;
    717         G4cout << " generated string momenta:   quark "
    718                << start->Get4Momentum() << "mass : "
    719                <<start->Get4Momentum().mag()           << G4endl;
    720         G4cout << " generated string momenta: Diquark "
    721                << end ->Get4Momentum()
    722                << "mass : " <<end->Get4Momentum().mag()<< G4endl;
    723         G4cout << " sum of ends                       " << Pstart+Pend << G4endl;
    724         G4cout << " Original                          " << hadron->Get4Momentum() << G4endl;
     1305          G4cout << " generated string flavors          "
     1306                 << start->GetPDGcode() << " / "
     1307                 << end->GetPDGcode()                    << G4endl;
     1308          G4cout << " generated string momenta:   quark "
     1309                 << start->Get4Momentum() << "mass : "
     1310                 <<start->Get4Momentum().mag()           << G4endl;
     1311          G4cout << " generated string momenta: Diquark "
     1312                 << end ->Get4Momentum()
     1313                 << "mass : " <<end->Get4Momentum().mag()<< G4endl;
     1314          G4cout << " sum of ends                       " << Pstart+Pend << G4endl;
     1315          G4cout << " Original                          " << hadron->Get4Momentum() << G4endl;
    7251316#endif
    7261317
    727         return string;
     1318        return;
     1319   
    7281320}
    7291321
    … …  
    7321324
    7331325// ---------------------------------------------------------------------
    734 G4double G4DiffractiveExcitation::ChooseP(G4double Pmin, G4double Pmax) const // Uzhi
     1326G4double G4DiffractiveExcitation::ChooseP(G4double Pmin, G4double Pmax) const
    7351327{
    7361328// choose an x between Xmin and Xmax with P(x) ~ 1/x
    7371329//  to be improved...
    7381330
    739         G4double range=Pmax-Pmin;                                         // Uzhi
     1331        G4double range=Pmax-Pmin;                   
    7401332
    7411333        if ( Pmin <= 0. || range <=0. )
    … …  
    7451337        }
    7461338
    747         G4double P;
    748 /*                                                                          // Uzhi
    749         do {
    750             x=Xmin + G4UniformRand() * range;
    751         }  while ( Xmin/x < G4UniformRand() );
    752 */                                                                          // Uzhi
    753 
    754         P=Pmin * std::pow(Pmax/Pmin,G4UniformRand());                       // Uzhi
    755 
    756 //debug-hpw     cout << "DiffractiveX "<<x<<G4endl;
     1339        G4double P=Pmin * std::pow(Pmax/Pmin,G4UniformRand());
    7571340        return P;
    7581341}
    … …  
    7601343// ---------------------------------------------------------------------
    7611344G4ThreeVector G4DiffractiveExcitation::GaussianPt(G4double AveragePt2,
    762                                                   G4double maxPtSquare) const // Uzhi
     1345                                                  G4double maxPtSquare) const
    7631346{            //  @@ this method is used in FTFModel as well. Should go somewhere common!
    7641347
    7651348        G4double Pt2;
    766 /*                                                                          // Uzhi
    767         do {
    768             pt2=widthSquare * std::log( G4UniformRand() );
    769         } while ( pt2 > maxPtSquare);
    770 */                                                                          // Uzhi
    771 
    7721349        Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() *
    773                            (std::exp(-maxPtSquare/AveragePt2)-1.));// Uzhi
     1350                           (std::exp(-maxPtSquare/AveragePt2)-1.));
    7741351
    7751352        G4double Pt=std::sqrt(Pt2);
    776 
    7771353        G4double phi=G4UniformRand() * twopi;
    778 
    7791354        return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
    7801355}
    7811356
     1357// ---------------------------------------------------------------------
     1358G4double G4DiffractiveExcitation::GetQuarkFractionOfKink(G4double zmin, G4double zmax) const
     1359    {
     1360       G4double z, yf;
     1361       do {
     1362           z = zmin + G4UniformRand()*(zmax-zmin);
     1363           yf = z*z +sqr(1 - z);       
     1364           }
     1365       while (G4UniformRand() > yf);
     1366       return z;
     1367    }
     1368// ---------------------------------------------------------------------
     1369void G4DiffractiveExcitation::UnpackMeson(const G4int IdPDG, G4int &Q1, G4int &Q2) const // Uzhi 7.09.09
     1370    {
     1371       G4int absIdPDG = std::abs(IdPDG);
     1372       Q1=  absIdPDG/ 100;
     1373       Q2= (absIdPDG %100)/10;
     1374           
     1375       G4int anti= 1 -2 * ( std::max( Q1, Q2 ) % 2 );
     1376
     1377       if (IdPDG < 0 ) anti *=-1;
     1378       Q1 *= anti;
     1379       Q2 *= -1 * anti;
     1380       return;
     1381    }
     1382// ---------------------------------------------------------------------
     1383void G4DiffractiveExcitation::UnpackBaryon(G4int IdPDG,
     1384                              G4int &Q1, G4int &Q2, G4int &Q3) const // Uzhi 7.09.09
     1385    {
     1386       Q1 =  IdPDG           / 1000;
     1387       Q2 = (IdPDG % 1000)  / 100;
     1388       Q3 = (IdPDG % 100)   / 10;
     1389       return;
     1390    }
     1391// ---------------------------------------------------------------------
     1392G4int G4DiffractiveExcitation::NewNucleonId(G4int Q1, G4int Q2, G4int Q3) const // Uzhi 7.09.09
     1393    {
     1394       G4int TmpQ(0);
     1395       if( Q3 > Q2 )
     1396       {
     1397        TmpQ = Q2;
     1398        Q2 = Q3;
     1399        Q3 = TmpQ;
     1400       } else if( Q3 > Q1 )
     1401       {
     1402        TmpQ = Q1;
     1403        Q1 = Q3;
     1404        Q3 = TmpQ;
     1405       }
     1406
     1407       if( Q2 > Q1 )
     1408       {
     1409        TmpQ = Q1;
     1410        Q1 = Q2;
     1411        Q2 = TmpQ;
     1412       }
     1413
     1414       G4int NewCode = Q1*1000 + Q2* 100 + Q3*  10 + 2;
     1415       return NewCode;
     1416    }
    7821417// ---------------------------------------------------------------------
    7831418G4DiffractiveExcitation::G4DiffractiveExcitation(const G4DiffractiveExcitation &)

  • trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4DiffractiveHHScatterer.cc

    r962 r1196  
    3838{}
    3939
    40 // -------------------------------------------------------------------
     40void G4DiffractiveHHScatterer::CreateStrings()
     41/*
     42                                             G4VSplitableHadron * aHadron,
     43                                             G4bool isProjectile,
     44                                             G4ExcitedString * FirstString,
     45                                             G4ExcitedString * SecondString,
     46                                             G4FTFParameters *theParameters)
     47*/
     48const
     49
     50{}
     51
     52/* -------------------------------------------------------------------
    4153G4KineticTrackVector * G4DiffractiveHHScatterer::
    4254Scatter(const G4KineticTrack & aTrack, const G4KineticTrack & bTrack)
    … …  
    7688  return result;
    7789}
     90*/

  • trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4DiffractiveSplitableHadron.cc

    r1007 r1196  
    2525//
    2626//
    27 // $Id: G4DiffractiveSplitableHadron.cc,v 1.7 2008/03/31 15:34:01 vuzhinsk Exp $
    28 // GEANT4 tag $Name: geant4-09-02 $
     27// $Id: G4DiffractiveSplitableHadron.cc,v 1.8 2009/07/31 11:03:00 vuzhinsk Exp $
     28// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2929//
    3030
    … …  
    8787       
    8888        G4int PDGcode=GetDefinition()->GetPDGEncoding();
     89
    8990        G4int stringStart, stringEnd;
    9091        ChooseStringEnds(PDGcode, &stringStart,&stringEnd);
    91        
     92
    9293        Parton[0] = new G4Parton(stringStart);
    9394        Parton[1] = new G4Parton(stringEnd);
    9495        PartonIndex=-1;
    95        
    9696}
    9797

  • trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4ElasticHNScattering.cc

    r968 r1196  
    2525//
    2626//
    27 // $Id: G4ElasticHNScattering.cc,v 1.3 2008/05/19 12:56:36 vuzhinsk Exp $
     27// $Id: G4ElasticHNScattering.cc,v 1.11 2009/10/06 10:10:36 vuzhinsk Exp $
    2828// ------------------------------------------------------------
    2929//      GEANT 4 class implemetation file
    … …  
    4646#include "G4VSplitableHadron.hh"
    4747#include "G4ExcitedString.hh"
    48 #include "G4FTFParameters.hh"                            // Uzhi 29.03.08
     48#include "G4FTFParameters.hh"
    4949//#include "G4ios.hh"
    5050
    … …  
    5959{
    6060//G4cout<<"G4ElasticHNScattering::ElasticScattering"<<G4endl;
    61 
     61// -------------------- Projectile parameters -----------------------------------
    6262           G4LorentzVector Pprojectile=projectile->Get4Momentum();
    63 
    64 // -------------------- Projectile parameters -----------------------------------
    65            G4bool PutOnMassShell=0;
     63//G4cout<<"Elastic P "<<Pprojectile<<G4endl;
     64           if(Pprojectile.z() < 0.)
     65           {
     66            target->SetStatus(2);
     67            return false;
     68           }
     69
     70//           G4double ProjectileRapidity = Pprojectile.rapidity();
     71//           G4int    ProjectilePDGcode=projectile->GetDefinition()->GetPDGEncoding();
     72//           G4ParticleDefinition * ProjectileDefinition = projectile->GetDefinition();
     73
     74           G4bool PutOnMassShell(false);
    6675
    6776           G4double M0projectile = Pprojectile.mag();                       
    68 
    6977           if(M0projectile < projectile->GetDefinition()->GetPDGMass())
    70              {
    71               PutOnMassShell=1;
     78           {
     79              PutOnMassShell=true;
    7280              M0projectile=projectile->GetDefinition()->GetPDGMass();
    73              }
     81           }
    7482
    7583           G4double Mprojectile2 = M0projectile * M0projectile;
    7684
    77 //           G4double AveragePt2=theParameters->GetSlope(); // Uzhi ???
    78 //           AveragePt2 = AveragePt2 * GeV*GeV;
    79 
    8085           G4double AveragePt2=theParameters->GetAvaragePt2ofElasticScattering();
    8186
    8287// -------------------- Target parameters ----------------------------------------------
     88//           G4int    TargetPDGcode=target->GetDefinition()->GetPDGEncoding();
     89
    8390           G4LorentzVector Ptarget=target->Get4Momentum();
    84 
     91//G4cout<<"Elastic T "<<Ptarget<<G4endl;
     92//           G4double TargetRapidity = Ptarget.rapidity();
    8593           G4double M0target = Ptarget.mag();
    86 //G4cout<<" Mp Mt Pt2 "<<M0projectile<<" "<<M0target<<" "<<AveragePt2/GeV/GeV<<G4endl;
    8794
    8895           if(M0target < target->GetDefinition()->GetPDGMass())
    89              {
    90               PutOnMassShell=1;
     96           {
     97              PutOnMassShell=true;
    9198              M0target=target->GetDefinition()->GetPDGMass();
    92              }
     99           }
    93100     
    94            G4double Mtarget2 = M0target * M0target;                      //Ptarget.mag2();
    95                                                                          // for AA-inter.
     101           G4double Mtarget2 = M0target * M0target;                     
     102
    96103// Transform momenta to cms and then rotate parallel to z axis;
    97104
    … …  
    103110           G4LorentzVector Ptmp=toCms*Pprojectile;
    104111
    105            if ( Ptmp.pz() <= 0. )                                 // Uzhi ???
     112           if ( Ptmp.pz() <= 0. )                               
    106113           {
    107114           // "String" moving backwards in  CMS, abort collision !!
    108115           //G4cout << " abort Collision!! " << G4endl;
     116                   target->SetStatus(2);
    109117                   return false;
    110118           }
    … …  
    118126           Ptarget.transform(toCms);
    119127
    120 // ---------------------- Sampling of transfered Pt ------------------------
     128// ---------------------- Putting on mass-on-shell, if needed ------------------------
     129           G4double PZcms2, PZcms;                                         
     130
     131           G4double S=Psum.mag2();                                         
     132//         G4double SqrtS=std::sqrt(S);                                     
     133
     134           PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
     135                                 2*S*Mprojectile2-2*S*Mtarget2-2*Mprojectile2*Mtarget2)/4./S;
     136
     137           if(PZcms2 < 0.)
     138           {  // It can be in an interaction with off-shell nuclear nucleon
     139            if(M0projectile > projectile->GetDefinition()->GetPDGMass())
     140            {  // An attempt to de-excite the projectile
     141               // It is assumed that the target is in the ground state
     142             M0projectile = projectile->GetDefinition()->GetPDGMass();
     143             Mprojectile2=M0projectile*M0projectile;
     144             PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
     145                    2*S*Mprojectile2 - 2*S*Mtarget2 - 2*Mprojectile2*Mtarget2)
     146                    /4./S;
     147             if(PZcms2 < 0.){ return false;} // Non succesful attempt after the de-excitation
     148            }
     149            else // if(M0projectile > projectile->GetDefinition()->GetPDGMass())
     150            {
     151             target->SetStatus(2);                                   // Uzhi 17.07.09
     152             return false;                   // The projectile was not excited,
     153                                             // but the energy was too low to put
     154                                             // the target nucleon on mass-shell
     155            }   // end of if(M0projectile > projectile->GetDefinition()->GetPDGMass())
     156           }    // end of if(PZcms2 < 0.)
     157
     158           PZcms = std::sqrt(PZcms2);
     159
     160           if(PutOnMassShell)
     161           {
     162              if(Pprojectile.z() > 0.)
     163              {
     164                 Pprojectile.setPz( PZcms);
     165                 Ptarget.setPz(    -PZcms);
     166              }
     167              else  // if(Pprojectile.z() > 0.)
     168              {
     169                 Pprojectile.setPz(-PZcms);
     170                 Ptarget.setPz(     PZcms);
     171              };
     172
     173              Pprojectile.setE(std::sqrt(Mprojectile2+
     174                                                      Pprojectile.x()*Pprojectile.x()+
     175                                                      Pprojectile.y()*Pprojectile.y()+
     176                                                      PZcms2));
     177              Ptarget.setE(std::sqrt(    Mtarget2    +
     178                                                      Ptarget.x()*Ptarget.x()+
     179                                                      Ptarget.y()*Ptarget.y()+
     180                                                      PZcms2));
     181           }  // end of if(PutOnMassShell)
     182
     183           G4double maxPtSquare = PZcms2;
     184
     185//----- Charge exchange between the projectile and the target, if possible
     186//        (G4UniformRand() < 0.5*std::exp(-0.5*(ProjectileRapidity - TargetRapidity))))
     187/*
     188           if((ProjectilePDGcode != TargetPDGcode) &&
     189             ((ProjectilePDGcode > 1000) && (TargetPDGcode > 1000)) &&
     190             (G4UniformRand() < 1.0*std::exp(-0.5*(ProjectileRapidity - TargetRapidity))))
     191           {
     192              projectile->SetDefinition(target->GetDefinition());
     193              target->SetDefinition(ProjectileDefinition);
     194           }
     195*/
     196// ------ Now we can calculate the transfered Pt --------------------------
    121197           G4double Pt2;                                                   
    122198           G4double ProjMassT2, ProjMassT;                                 
    123            G4double TargMassT2, TargMassT;                                 
    124            G4double PZcms2, PZcms;                                         
    125 
    126            G4double S=Psum.mag2();                                         
    127 //           G4double SqrtS=std::sqrt(S);                                     
    128 
    129            PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
    130                                  2*S*Mprojectile2-2*S*Mtarget2-2*Mprojectile2*Mtarget2)/4./S;
    131            if(PZcms2 < 0)
    132              {return false;}   // It can be in an interaction with off-shell nuclear nucleon
    133 
    134            PZcms = std::sqrt(PZcms2);
    135 
    136            if(PutOnMassShell)
    137              {
    138               if(Pprojectile.z() > 0.)
    139                 {
    140                  Pprojectile.setPz( PZcms);
    141                  Ptarget.setPz(    -PZcms);
    142                 }
    143               else
    144                  {
    145                  Pprojectile.setPz(-PZcms);
    146                  Ptarget.setPz(     PZcms);
    147                 };
    148 
    149                Pprojectile.setE(std::sqrt(Mprojectile2+
    150                                                        Pprojectile.x()*Pprojectile.x()+
    151                                                        Pprojectile.y()*Pprojectile.y()+
    152                                                        PZcms2));
    153                Ptarget.setE(std::sqrt(    Mtarget2    +
    154                                                        Ptarget.x()*Ptarget.x()+
    155                                                        Ptarget.y()*Ptarget.y()+
    156                                                        PZcms2));
    157              }
    158 
    159            G4double maxPtSquare = PZcms2;
     199           G4double TargMassT2, TargMassT;
     200
    160201           G4LorentzVector Qmomentum;
    161202
    … …  
    163204
    164205           Pt2=G4ThreeVector(Qmomentum.vect()).mag2();                 
    165 //G4cout<<"Pt2 GeV^2 "<<(Pt2)/GeV/GeV<<G4endl;
    166206
    167207           ProjMassT2=Mprojectile2+Pt2;                           
    … …  
    175215                    2.*S*ProjMassT2-2.*S*TargMassT2-                 
    176216                    2.*ProjMassT2*TargMassT2)/4./S;                 
    177            if(PZcms2 < 0 ) {PZcms2=0;};                                 
     217           if(PZcms2 < 0 ) {PZcms2=0;};// to avoid the exactness problem
    178218           PZcms =std::sqrt(PZcms2);                                   
    179219
    180            Pprojectile.setPz( PZcms);  // Uzhi Proj can move backward
    181            Ptarget.setPz(    -PZcms);  // Uzhi Proj can move backward
    182 
    183 //G4cout << "Qplus / Qminus " << Qplus << " / " << Qminus<<G4endl;
    184 //           G4cout << "pt2" << pt2 << G4endl;
    185 //           G4cout << "Qmomentum " << Qmomentum << G4endl;
    186 //           G4cout << " Masses (P/T) : " << (Pprojectile+Qmomentum).mag() <<
    187 //                             " / " << (Ptarget-Qmomentum).mag() << G4endl;
     220           Pprojectile.setPz( PZcms); 
     221           Ptarget.setPz(    -PZcms);
    188222
    189223           Pprojectile += Qmomentum;
    190224           Ptarget     -= Qmomentum;
    191 
    192 //G4cout << "Pprojectile with Q : " << Pprojectile << G4endl;
    193 //G4cout << "Ptarget with Q : " << Ptarget << G4endl;
    194        
    195 //         G4cout << "Projectile back: " << toLab * Pprojectile << G4endl;
    196 //         G4cout << "Target back: " << toLab * Ptarget << G4endl;
    197225
    198226// Transform back and update SplitableHadron Participant.
    199227           Pprojectile.transform(toLab);
    200228           Ptarget.transform(toLab);
    201 
    202 //G4cout << "Pprojectile with Q M: " << Pprojectile<<" "<<  Pprojectile.mag() << G4endl;
    203 //G4cout << "Ptarget     with Q M: " << Ptarget    <<" "<<  Ptarget.mag()     << G4endl;
    204 
    205 //G4cout << "Target      mass  " <<  Ptarget.mag() << G4endl;     
    206                            
    207 //G4cout << "Projectile mass  " <<  Pprojectile.mag() << G4endl;
    208 
    209            G4double ZcoordinateOfCurrentInteraction = target->GetPosition().z();
    210 // It is assumed that nucleon z-coordinates are ordered on increasing -----------
    211 
    212            G4double betta_z=projectile->Get4Momentum().pz()/projectile->Get4Momentum().e();
    213 
    214            G4double ZcoordinateOfPreviousCollision=projectile->GetPosition().z();
    215            if(projectile->GetSoftCollisionCount()==0) {
    216               projectile->SetTimeOfCreation(0.);
    217               target->SetTimeOfCreation(0.);
    218               ZcoordinateOfPreviousCollision=ZcoordinateOfCurrentInteraction;
    219            }
    220            
    221            G4ThreeVector thePosition(projectile->GetPosition().x(),
    222                                      projectile->GetPosition().y(),
    223                                      ZcoordinateOfCurrentInteraction);
    224            projectile->SetPosition(thePosition);
    225 
    226            G4double TimeOfPreviousCollision=projectile->GetTimeOfCreation();
    227            G4double TimeOfCurrentCollision=TimeOfPreviousCollision+
    228                     (ZcoordinateOfCurrentInteraction-ZcoordinateOfPreviousCollision)/betta_z;
    229 
    230            projectile->SetTimeOfCreation(TimeOfCurrentCollision);
    231            target->SetTimeOfCreation(TimeOfCurrentCollision);
     229/*  // Maybe it will be needed for an exact calculations--------------------
     230           G4double TargetMomentum=std::sqrt(Ptarget.x()*Ptarget.x()+
     231                                             Ptarget.y()*Ptarget.y()+
     232                                             Ptarget.z()*Ptarget.z());
     233*/
     234
     235// Calculation of the creation time ---------------------
     236      projectile->SetTimeOfCreation(target->GetTimeOfCreation());
     237      projectile->SetPosition(target->GetPosition());
     238// Creation time and position of target nucleon were determined at
     239// ReggeonCascade() of G4FTFModel
     240// ------------------------------------------------------
    232241
    233242           projectile->Set4Momentum(Pprojectile);

  • trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4FTFModel.cc

    r1007 r1196  
    2525//
    2626//
    27 // $Id: G4FTFModel.cc,v 1.13 2008/12/09 10:40:52 vuzhinsk Exp $
    28 // GEANT4 tag $Name: geant4-09-02 $
     27// $Id: G4FTFModel.cc,v 1.28 2009/10/29 14:55:33 vuzhinsk Exp $
     28// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2929//
    3030
    … …  
    3838
    3939#include "G4FTFModel.hh"
    40 #include "G4FTFParameters.hh"                            // Uzhi 29.03.08
     40#include "G4FTFParameters.hh"
    4141#include "G4FTFParticipants.hh"
     42#include "G4DiffractiveSplitableHadron.hh"
    4243#include "G4InteractionContent.hh"
    4344#include "G4LorentzRotation.hh"
    4445#include "G4ParticleDefinition.hh"
     46#include "G4ParticleTable.hh"
    4547#include "G4ios.hh"
    46 #include <utility>                                        // Uzhi 29.03.08
     48#include <utility>
     49#include "G4IonTable.hh"
    4750
    4851// Class G4FTFModel
    4952
    5053G4FTFModel::G4FTFModel():theExcitation(new G4DiffractiveExcitation()),
    51                          theElastic(new G4ElasticHNScattering()) // Uzhi 29.03.08
     54                         theElastic(new G4ElasticHNScattering())
    5255{
    5356        G4VPartonStringModel::SetThisPointer(this);
    54         theParameters=0;                                         // Uzhi 9.12.08
    55 }
    56 
    57 /*
    58 G4FTFModel::G4FTFModel(G4double , G4double , G4double ):theExcitation(new // Uzhi 9.12.08 G4DiffractiveExcitation())
    59 {
    60         G4VPartonStringModel::SetThisPointer(this);
    61 }
    62 
    63 G4FTFModel::G4FTFModel(G4DiffractiveExcitation * anExcitation)
    64 :
    65 theExcitation(anExcitation)
    66 {
    67         G4VPartonStringModel::SetThisPointer(this);
    68 }
    69 */
     57        theParameters=0;
     58}
    7059
    7160
    7261G4FTFModel::~G4FTFModel()
    7362{
    74    if( theParameters != 0 ) delete theParameters;             // Uzhi 5.12.08
     63   if( theParameters != 0 ) delete theParameters;
    7564// Because FTF model can be called for various particles
    7665// theParameters must be erased at the end of each call.
    77 // Thus the delete is olso in G4FTFModel::GetStrings() method
    78    if( theExcitation != 0 ) delete theExcitation;             // Uzhi 5.12.08
    79    if( theElastic    != 0 ) delete theElastic;                // Uzhi 5.12.08
    80 }
    81 
     66// Thus the delete is also in G4FTFModel::GetStrings() method
     67   if( theExcitation != 0 ) delete theExcitation;
     68   if( theElastic    != 0 ) delete theElastic;
     69}
    8270
    8371const G4FTFModel & G4FTFModel::operator=(const G4FTFModel &)
    … …  
    8775}
    8876
    89 
    9077int G4FTFModel::operator==(const G4FTFModel &right) const
    9178{
    … …  
    10289{
    10390        theProjectile = aProjectile; 
    104 //G4cout<<"G4FTFModel::Init "<<aNucleus.GetN()<<" "<<aNucleus.GetZ()<<G4endl;
    10591        theParticipants.Init(aNucleus.GetN(),aNucleus.GetZ());
    106 // Uzhi N-mass number Z-charge ------------------------- Uzhi 29.03.08
     92// ----------- N-mass number Z-charge -------------------------
    10793
    10894// --- cms energy
    … …  
    11298                     2*theProjectile.GetTotalEnergy()*G4Proton::Proton()->GetPDGMass();
    11399/*
     100G4cout<<"----------------------------------------"<<G4endl;
    114101G4cout << " primary Total E (GeV): " << theProjectile.GetTotalEnergy()/GeV << G4endl;
    115102G4cout << " primary Mass    (GeV): " << theProjectile.GetMass() /GeV << G4endl;
     103G4cout << " primary 3Mom           " << theProjectile.GetMomentum() << G4endl;
     104G4cout << " primary space position " << theProjectile.GetPositionInNucleus() << G4endl;
    116105G4cout << "cms std::sqrt(s) (GeV) = " << std::sqrt(s) / GeV << G4endl;
    117106*/
    118107
    119       if( theParameters != 0 ) delete theParameters;                    // Uzhi 9.12.08
     108      if( theParameters != 0 ) delete theParameters;
    120109      theParameters = new G4FTFParameters(theProjectile.GetDefinition(),
    121110                                          aNucleus.GetN(),aNucleus.GetZ(),
    122                                           s);// ------------------------- Uzhi 19.04.08
    123 //theParameters->SetProbabilityOfElasticScatt(0.); // To turn on/off (1/0) elastic scattering
     111                                          s);
     112//theParameters->SetProbabilityOfElasticScatt(0.);
     113// To turn on/off (1/0) elastic scattering
     114
    124115}
    125116
    … …  
    127118G4ExcitedStringVector * G4FTFModel::GetStrings()
    128119{
    129 //G4cout<<"theParticipants.GetList"<<G4endl;
     120//G4int Uzhi; G4cin>>Uzhi;
     121 
     122//G4cout<<"GetList"<<G4endl;
    130123        theParticipants.GetList(theProjectile,theParameters);
    131 //G4cout<<"ExciteParticipants()"<<G4endl;
    132         if (! ExciteParticipants()) return NULL;;
    133 //G4cout<<"theStrings = BuildStrings()"<<G4endl;
     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;
    134131        G4ExcitedStringVector * theStrings = BuildStrings();
    135 //G4cout<<"Return to theStrings "<<G4endl;
    136         if( theParameters != 0 )                              // Uzhi 9.12.08
    137         {                                                     // Uzhi 9.12.08
    138           delete theParameters;                               // Uzhi 9.12.08
    139           theParameters=0;                                    // Uzhi 9.12.08
    140         }                                                     // Uzhi 9.12.08
     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        }
    141141        return theStrings;
    142142}
    143143
    144144// ------------------------------------------------------------
    145 struct DeleteVSplitableHadron { void operator()(G4VSplitableHadron * aH){delete aH;} };
     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}
    146557
    147558// ------------------------------------------------------------
    148559G4bool G4FTFModel::ExciteParticipants()
    149560{
    150 /*    // Uzhi 29.03.08                     For elastic Scatt.
    151 G4cout<<"  In ExciteParticipants() "<<theParticipants.theInteractions.size()<<G4endl;
    152 G4cout<<" test Params Tot "<<theParameters->GetTotalCrossSection()<<G4endl;
    153 G4cout<<" test Params Ela "<<theParameters->GetElasticCrossSection()<<G4endl;
     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;
    154565       
    155 G4int counter=0;
    156 */   // Uzhi 29.03.08
    157 
    158 
    159 
     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
    160576        while (theParticipants.Next())
    161577        {         
    162578           const G4InteractionContent & collision=theParticipants.GetInteraction();
    163 /*
    164 counter++;
    165 G4cout<<" Inter # "<<counter<<G4endl;
    166 */
     579//
     580//counter++;
     581//G4cout<<" Int num "<<counter<<G4endl;
     582//
    167583           G4VSplitableHadron * projectile=collision.GetProjectile();
    168584           G4VSplitableHadron * target=collision.GetTarget();
    169 
    170 //   // Uzhi 29.03.08
    171            G4bool Successfull;
     585//         G4Nucleon * TargetNucleon=collision.GetTargetNucleon(); // Uzhi 16.07.09
     586// Uzhi 16.07.09 ----------------------------
    172587           if(G4UniformRand()< theParameters->GetProbabilityOfElasticScatt())
    173            {
     588           { //   Elastic scattering -------------------------
    174589//G4cout<<"Elastic"<<G4endl;
    175             Successfull=theElastic->ElasticScattering(projectile, target, theParameters);
     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            };
    176606           }
    177607           else
    178            {
    179 //G4cout<<"Inelastic"<<G4endl;
    180             Successfull=theExcitation->ExciteParticipants(projectile, target, theParameters);
     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            };
    181627           }
    182 //           if(!Successfull)
    183 // // Uzhi 29.03.08
    184 
    185 //         if ( ! theExcitation->ExciteParticipants(projectile, target) )
    186            if(!Successfull)
    187            {
     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;
    188634//           give up, clean up
    189                 std::vector<G4VSplitableHadron *> primaries;
    190                 std::vector<G4VSplitableHadron *> targets;
    191                 theParticipants.StartLoop();    // restart a loop
    192                 while ( theParticipants.Next() )
    193                 {
    194                     const G4InteractionContent & interaction=theParticipants.GetInteraction();
     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();
    195641                         //  do not allow for duplicates ...
    196                     if ( primaries.end() == std::find(primaries.begin(), primaries.end(),
    197                                                       interaction.GetProjectile()) )
    198                         primaries.push_back(interaction.GetProjectile());
    199 
    200                     if ( targets.end()   == std::find(targets.begin(), targets.end(),
     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(),
    201647                                                      interaction.GetTarget()) )
    202                         targets.push_back(interaction.GetTarget());
    203                 }
    204                 std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron());
    205                 primaries.clear();
    206                 std::for_each(targets.begin(), targets.end(), DeleteVSplitableHadron());
    207                 targets.clear();
    208 
    209                 return false;
    210            }  // End of the loop Uzhi
    211         }
     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
    212676        return true;
    213677}
    … …  
    217681// Loop over all collisions; find all primaries, and all target ( targets may
    218682//  be duplicate in the List ( to unique G4VSplitableHadrons)
     683
     684//G4cout<<"In build string"<<G4endl;
    219685
    220686        G4ExcitedStringVector * strings;
    … …  
    223689        std::vector<G4VSplitableHadron *> primaries;
    224690        std::vector<G4VSplitableHadron *> targets;
     691        std::vector<G4Nucleon          *> TargetNucleons;     // Uzhi 16.07.09
    225692       
     693        G4ExcitedString * FirstString(0);     // If there will be a kink,
     694        G4ExcitedString * SecondString(0);    // two strings will be prodused.
     695
    226696        theParticipants.StartLoop();    // restart a loop
     697//G4int InterCount(0); // Uzhi
    227698        while ( theParticipants.Next() )
    228699        {
    229700            const G4InteractionContent & interaction=theParticipants.GetInteraction();
    230701                 //  do not allow for duplicates ...
     702
    231703            if ( primaries.end() == std::find(primaries.begin(), primaries.end(),
    232704                                                interaction.GetProjectile()) )
    … …  
    236708                                                interaction.GetTarget()) )
    237709                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++;
    238716        }
    239717           
    240718       
    241 //      G4cout << "BuildStrings prim/targ " << primaries.size() << " , " <<
    242 //                                             targets.size() << G4endl;
     719//G4cout << "BuildStrings prim/targ/woundN " << primaries.size() << " , " <<targets.size() <<", "<<TargetNucleons.size()<< G4endl;
    243720
    244721        unsigned int ahadron;
    245722// Only for hA-interactions Uzhi -------------------------------------
     723//G4int StringN(0);
     724//G4cout<<"Proj strings -----------------------"<<G4endl;
    246725        for ( ahadron=0; ahadron < primaries.size() ; ahadron++)
    247726        {
    248 //G4ThreeVector aPosition=primaries[ahadron]->GetPosition();
    249 //G4cout<<"Proj Build "<<aPosition<<" "<<primaries[ahadron]->GetTimeOfCreation()<<G4endl;
    250             G4bool isProjectile=true;
    251             strings->push_back(theExcitation->String(primaries[ahadron], isProjectile));
     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;
    252741        }
    253 
     742//G4cout<<"Targ strings ------------------------------"<<G4endl;
    254743        for ( ahadron=0; ahadron < targets.size() ; ahadron++)
    255744        {
    256 //G4ThreeVector aPosition=targets[ahadron]->GetPosition();
    257 //G4cout<<"Targ Build "<<aPosition<<" "<<targets[ahadron]->GetTimeOfCreation()<<G4endl;
    258             G4bool isProjectile=false;
    259             strings->push_back(theExcitation->String(targets[ahadron], isProjectile));
     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            }
    260765        }
     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;
    261808
    262809        std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron());
    263810        primaries.clear();
     811
    264812        std::for_each(targets.begin(), targets.end(), DeleteVSplitableHadron());
    265813        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!
    266838       
    267         return strings;
    268 }
    269 // ------------------------------------------------------------
     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}

  • trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4FTFParameters.cc

    r1007 r1196  
    1 //
    2 // ********************************************************************
     1//********************************************************************
    32// * License and Disclaimer                                           *
    43// *                                                                  *
    … …  
    2524//
    2625//
    27 // $Id: G4FTFParameters.cc,v 1.4 2008/12/18 13:02:00 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4FTFParameters.cc,v 1.11 2009/10/25 10:50:54 vuzhinsk Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2928//
    3029
    3130#include "G4FTFParameters.hh"
     31
     32#include "G4ios.hh"
     33#include <utility>                                        // Uzhi 29.03.08
    3234
    3335G4FTFParameters::G4FTFParameters()
    … …  
    4345                                                   G4double   theZ,
    4446                                                   G4double   s)
    45     {
     47{
    4648    G4int PDGcode = particle->GetPDGEncoding();
    4749    G4int absPDGcode = std::abs(PDGcode);
    48     G4double Elab = (s - 2*0.88*GeV*GeV)/(2*0.939*GeV)/GeV;
    49     G4double Plab = std::sqrt(Elab * Elab - 0.88);
     50    G4double ProjectileMass = particle->GetPDGMass();
     51    G4double TargetMass     = G4Proton::Proton()->GetPDGMass();
     52
     53    G4double Elab = (s - ProjectileMass*ProjectileMass - TargetMass*TargetMass)/
     54                     (2*TargetMass);
     55    G4double Plab = std::sqrt(Elab * Elab - ProjectileMass*ProjectileMass);
    5056
    5157    G4double LogPlab = std::log( Plab );
    5258    G4double sqrLogPlab = LogPlab * LogPlab;
    53 
    54 //G4cout<<"G4FTFParameters Plab "<<Plab<<G4endl;
    5559
    5660    G4int NumberOfTargetProtons  = (G4int) theZ;
    … …  
    144148                           NumberOfTargetNeutrons * XelKN   ) / NumberOfTargetNucleons;
    145149      }
    146     else if( PDGcode == 311 )                      //------Projectile is KaonZero ------
     150    else if((PDGcode == 311) || (PDGcode == 130) || (PDGcode == 310))//Projectile is KaonZero
    147151      {
    148152       G4double XtotKP =( 18.1 +  0. *std::pow(Plab, 0.  ) + 0.26 *sqrLogPlab - 1.0 *LogPlab +   //K+
    … …  
    181185      SetInelasticCrossSection(Xtotal-Xelastic);
    182186
    183 //  // Interactions with elastic ans inelastic collisions
     187//  // Interactions with elastic and inelastic collisions
    184188      SetProbabilityOfElasticScatt(Xtotal, Xelastic);
    185189      SetRadiusOfHNinteractions2(Xtotal/pi/10.);
    … …  
    190194*/ //=======================================================
    191195
    192 //G4cout<<" Rnn "<<Xtotal/pi/10.<<" "<<Xtotal/pi/10.*fermi*fermi<<G4endl;
    193 //G4cout<<"G4FTFParameters Xt Xel MeV "<<Xtotal<<" "<<Xelastic<<" "<<GeV<<G4endl;
    194 
    195196//----------------------------------------------------------------------------------- 
    196197      SetSlope( Xtotal*Xtotal/16./pi/Xelastic/0.3894 ); // Slope parameter of elastic scattering
    197198                                                        //      (GeV/c)^(-2))
    198 //G4cout<<"G4FTFParameters Slope "<<GetSlope()<<G4endl;
    199199//-----------------------------------------------------------------------------------
    200200      SetGamma0( GetSlope()*Xtotal/10./2./pi );
    … …  
    208208           if( absPDGcode > 1000 )                        //------Projectile is baryon --------
    209209             {
     210              SetMagQuarkExchange(3.4); //3.8);
     211              SetSlopeQuarkExchange(1.2);
     212              SetDeltaProbAtQuarkExchange(0.1); //(0.1*4.);
     213
    210214              SetProjMinDiffMass(1.1);                    // GeV
    211215              SetProjMinNonDiffMass(1.1);                 // GeV
    212               SetProbabilityOfProjDiff(0.95*std::pow(s/GeV/GeV,-0.35)); // 40/32 X-dif/X-inel
     216              SetProbabilityOfProjDiff(0.76*std::pow(s/GeV/GeV,-0.35));
    213217
    214218              SetTarMinDiffMass(1.1);                     // GeV
    215219              SetTarMinNonDiffMass(1.1);                  // GeV
    216               SetProbabilityOfTarDiff(0.95*std::pow(s/GeV/GeV,-0.35)); // 40/32 X-dif/X-inel
     220              SetProbabilityOfTarDiff(0.76*std::pow(s/GeV/GeV,-0.35));
    217221
    218222              SetAveragePt2(0.3);                         // GeV^2
    … …  
    220224           else if( absPDGcode == 211 || PDGcode ==  111) //------Projectile is Pion -----------
    221225             {
     226              SetMagQuarkExchange(120.); // 210.
     227              SetSlopeQuarkExchange(2.0);
     228              SetDeltaProbAtQuarkExchange(0.6);
     229
    222230              SetProjMinDiffMass(0.5);                    // GeV
    223231              SetProjMinNonDiffMass(0.3);                 // GeV
    224               SetProbabilityOfProjDiff(0.62*std::pow(s/GeV/GeV,-0.51)); // 40/32 X-dif/X-inel
    225 
     232              SetProbabilityOfProjDiff(0.*0.62*std::pow(s/GeV/GeV,-0.51)); // 40/32 X-dif/X-inel
     233//G4cout<<"Params "<<0.6*0.62*std::pow(s/GeV/GeV,-0.51)<<G4endl;
    226234              SetTarMinDiffMass(1.1);                     // GeV
    227235              SetTarMinNonDiffMass(1.1);                  // GeV
    228               SetProbabilityOfTarDiff(0.62*std::pow(s/GeV/GeV,-0.51)); // 40/32 X-dif/X-inel
     236//G4cout<<"       "<<2.7*0.62*std::pow(s/GeV/GeV,-0.51)<<G4endl; // was 2
     237//G4int Uzhi; G4cin>>Uzhi;
     238              SetProbabilityOfTarDiff(2.*0.62*std::pow(s/GeV/GeV,-0.51)); // 40/32 X-dif/X-inel
    229239
    230240/*
    … …  
    236246              SetAveragePt2(0.3);                         // GeV^2
    237247             }
    238            else if( absPDGcode == 321 || PDGcode == -311) //------Projectile is Kaon -----------
     248           else if( (absPDGcode == 321) || (PDGcode == 311) ||
     249                       (PDGcode == 130) || (PDGcode == 310))  //Projectile is Kaon
    239250             {
     251// Must be corrected, taken from PiN
     252              SetMagQuarkExchange(120.);
     253              SetSlopeQuarkExchange(2.0);
     254              SetDeltaProbAtQuarkExchange(0.6);
     255
    240256              SetProjMinDiffMass(0.7);                    // GeV 1.1
    241257              SetProjMinNonDiffMass(0.7);                 // GeV
    … …  
    251267                                                          //------Nucleon assumed
    252268             {
     269              SetMagQuarkExchange(3.5);
     270              SetSlopeQuarkExchange(1.0);
     271              SetDeltaProbAtQuarkExchange(0.1);
     272
    253273              SetProjMinDiffMass((particle->GetPDGMass()+160.*MeV)/GeV);
    254274              SetProjMinNonDiffMass((particle->GetPDGMass()+160.*MeV)/GeV);
    … …  
    260280
    261281              SetAveragePt2(0.3);                         // GeV^2
    262              };
    263 
    264 
    265 //G4cout<<"G4FTFParameters Out"<<G4endl;
    266 
    267     }
     282             }
     283
     284// ---------- Set parameters of a string kink -------------------------------
     285             SetPt2Kink(6.*GeV*GeV);
     286             G4double Puubar(1./3.), Pddbar(1./3.), Pssbar(1./3.); // SU(3) symmetry
     287//           G4double Puubar(0.41 ), Pddbar(0.41 ), Pssbar(0.18 ); // Broken SU(3) symmetry
     288             SetQuarkProbabilitiesAtGluonSplitUp(Puubar, Pddbar, Pssbar);
     289
     290// --------- Set parameters of nuclear destruction--------------------
     291
     292    if( absPDGcode < 1000 )
     293    {
     294      SetCofNuclearDestruction(1.); //1.0);                 // for meson projectile
     295    } else if( theA > 20. )
     296    {
     297      SetCofNuclearDestruction(0.2); //2);                 // for baryon projectile and heavy target
     298    } else
     299    {
     300      SetCofNuclearDestruction(0.2); //1.0);                 // for baryon projectile and light target
     301    }
     302
     303    SetR2ofNuclearDestruction(1.5*fermi*fermi);
     304
     305    SetExcitationEnergyPerWoundedNucleon(100*MeV);
     306
     307    SetDofNuclearDestruction(0.4);
     308    SetPt2ofNuclearDestruction(0.17*GeV*GeV);
     309    SetMaxPt2ofNuclearDestruction(1.0*GeV*GeV);
     310}
    268311//**********************************************************************************************

  • trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4FTFParticipants.cc

    r1007 r1196  
    2525//
    2626//
    27 // $Id: G4FTFParticipants.cc,v 1.9 2008/06/13 12:49:23 vuzhinsk Exp $
    28 // GEANT4 tag $Name: geant4-09-02 $
     27// $Id: G4FTFParticipants.cc,v 1.15 2009/10/25 10:50:54 vuzhinsk Exp $
     28// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2929//
    3030// ------------------------------------------------------------
    … …  
    8787    xyradius =theNucleus->GetOuterRadius() + deltaxy; // Impact parameter sampling
    8888                                                      // radius
    89     G4bool nucleusNeedsShift = true;
     89//    G4bool nucleusNeedsShift = true;                // Uzhi 20 July 2009
    9090   
    9191    while ( theInteractions.size() == 0 )
    … …  
    9696        G4double impactY = theImpactParameter.second;
    9797
     98        G4ThreeVector thePosition(impactX, impactY, -DBL_MAX);     //Uzhi 24 July 09
     99        primarySplitable->SetPosition(thePosition);                //Uzhi 24 July 09
     100
    98101        theNucleus->StartLoop();
    99102        G4Nucleon * nucleon;
    100103//G4int InterNumber=0;           // Uzhi
     104G4int NucleonNumber(0);        // Uzhi
    101105//while ( (nucleon=theNucleus->GetNextNucleon())&& (InterNumber < 1) ) // Uzhi
    102106        while ( (nucleon=theNucleus->GetNextNucleon()) ) // Uzhi
    103107        {
     108//G4cout<<"Nucl# "<<NucleonNumber<<G4endl; // Vova
    104109           G4double impact2= sqr(impactX - nucleon->GetPosition().x()) +
    105110                             sqr(impactY - nucleon->GetPosition().y());
    … …  
    110115           {
    111116//InterNumber++;
     117/*                                                    // Uzhi 20 July 2009
    112118                if ( nucleusNeedsShift )
    113119                {                       // on the first hit, shift nucleus
    … …  
    117123                     impactY=0;
    118124                }
    119                 G4VSplitableHadron * targetSplitable;
    120                 if ( (targetSplitable=nucleon->GetSplitableHadron()) == NULL )
     125*/                                                    // Uzhi 20 July 2009
     126//G4cout<<"          Interact"<<G4endl;
     127                primarySplitable->SetStatus(1);        // It takes part in the interaction
     128
     129                G4VSplitableHadron * targetSplitable=0;
     130                if ( ! nucleon->AreYouHit() )
    121131                {
     132//G4cout<<"Part "<<nucleon->Get4Momentum()<<G4endl;
    122133                    targetSplitable= new G4DiffractiveSplitableHadron(*nucleon);
    123134                    nucleon->Hit(targetSplitable);
     135                    nucleon->SetBindingEnergy(3.*nucleon->GetBindingEnergy()); // Uzhi 5.10.09
     136                    targetSplitable->SetStatus(1);     // It takes part in the interaction
    124137                }
    125138                G4InteractionContent * aInteraction =
    126139                                       new G4InteractionContent(primarySplitable);
    127140                aInteraction->SetTarget(targetSplitable);
     141                aInteraction->SetTargetNucleon(nucleon);     // Uzhi 16.07.09
    128142                theInteractions.push_back(aInteraction);
    129143           }
     144NucleonNumber++; // Uzhi
    130145        }   
    131 
    132 //      G4cout << "Number of Hit nucleons " << theInteractions.size() //  entries()
     146// // Uzhi
     147//      G4cout << "Number of Hit nucleons " << theInteractions.size()<<G4endl; //  entries()
    133148//              << "\t" << impactX/fermi << "\t"<<impactY/fermi
    134149//              << "\t" << std::sqrt(sqr(impactX)+sqr(impactY))/fermi <<G4endl;
    135 
     150// // Uzhi
    136151    }
    137152   
Note: See TracChangeset for help on using the changeset viewer.