Changeset 1196 for trunk/source/processes/hadronic/models/parton_string
- Timestamp:
- Nov 25, 2009, 5:13:58 PM (15 years ago)
- Location:
- trunk/source/processes/hadronic/models/parton_string
- Files:
-
- 49 edited
-
History (modified) (1 diff)
-
diffraction/History (modified) (2 diffs)
-
diffraction/include/G4DiffractiveExcitation.hh (modified) (3 diffs)
-
diffraction/include/G4DiffractiveHHScatterer.hh (modified) (3 diffs)
-
diffraction/include/G4DiffractiveSplitableHadron.hh (modified) (1 diff)
-
diffraction/include/G4ElasticHNScattering.hh (modified) (4 diffs)
-
diffraction/include/G4FTFModel.hh (modified) (6 diffs)
-
diffraction/include/G4FTFParameters.hh (modified) (11 diffs)
-
diffraction/include/G4FTFParticipants.hh (modified) (3 diffs)
-
diffraction/src/G4DiffractiveExcitation.cc (modified) (23 diffs)
-
diffraction/src/G4DiffractiveHHScatterer.cc (modified) (2 diffs)
-
diffraction/src/G4DiffractiveSplitableHadron.cc (modified) (2 diffs)
-
diffraction/src/G4ElasticHNScattering.cc (modified) (7 diffs)
-
diffraction/src/G4FTFModel.cc (modified) (9 diffs)
-
diffraction/src/G4FTFParameters.cc (modified) (11 diffs)
-
diffraction/src/G4FTFParticipants.cc (modified) (5 diffs)
-
hadronization/History (modified) (3 diffs)
-
hadronization/include/G4ExcitedStringDecay.hh (modified) (3 diffs)
-
hadronization/include/G4FragmentingString.hh (modified) (1 diff)
-
hadronization/include/G4HadronBuilder.hh (modified) (1 diff)
-
hadronization/include/G4LundStringFragmentation.hh (modified) (1 diff)
-
hadronization/include/G4QGSMFragmentation.hh (modified) (1 diff)
-
hadronization/include/G4VKinkyStringDecay.hh (modified) (1 diff)
-
hadronization/include/G4VLongitudinalStringDecay.hh (modified) (1 diff)
-
hadronization/src/G4ExcitedStringDecay.cc (modified) (1 diff)
-
hadronization/src/G4HadronBuilder.cc (modified) (1 diff)
-
hadronization/src/G4LundStringFragmentation.cc (modified) (7 diffs)
-
hadronization/src/G4QGSMFragmentation.cc (modified) (1 diff)
-
hadronization/src/G4VKinkyStringDecay.cc (modified) (1 diff)
-
hadronization/src/G4VLongitudinalStringDecay.cc (modified) (5 diffs)
-
management/History (modified) (2 diffs)
-
management/include/G4EventGenerator.hh (modified) (1 diff)
-
management/include/G4InteractionCode.hh (modified) (1 diff)
-
management/include/G4InteractionContent.hh (modified) (5 diffs)
-
management/include/G4PomeronCrossSection.hh (modified) (1 diff)
-
management/include/G4StringModel.hh (modified) (1 diff)
-
management/include/G4VParticipants.hh (modified) (3 diffs)
-
management/include/G4VPartonStringModel.hh (modified) (1 diff)
-
management/include/G4VSplitableHadron.hh (modified) (4 diffs)
-
management/include/G4VStringFragmentation.hh (modified) (1 diff)
-
management/include/G4VertexCode.hh (modified) (1 diff)
-
management/src/G4EventGenerator.cc (modified) (1 diff)
-
management/src/G4InteractionContent.cc (modified) (1 diff)
-
management/src/G4PomeronCrossSection.cc (modified) (1 diff)
-
management/src/G4StringModel.cc (modified) (1 diff)
-
management/src/G4VParticipants.cc (modified) (1 diff)
-
management/src/G4VPartonStringModel.cc (modified) (4 diffs)
-
management/src/G4VSplitableHadron.cc (modified) (6 diffs)
-
management/src/G4VStringFragmentation.cc (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/parton_string/History
r962 r1196 14 14 * Please list in reverse chronological order (last date on top) 15 15 --------------------------------------------------------------- 16 17 19 November 2009, Gunter Folger (hadr-prtn-V09-02-09) 18 - G4VParticipants.hh: use modified name of method to sort nucleons in G4V3DNucleus 19 20 29 October 2009, V. Uzhinsky (hadr-prtn-V09-02-08) 21 Warning meassage is erased in FTFmodel 22 ----------------------------------------------------- 23 25 October 2009, V Uzhinsky (hadr-prtn-V09-02-07) 24 Excitation energy calculation was added in FTF. 25 ----------------------------------------------------- 26 27 6 October 2009, V Uzhinsky (hadr-prtn-V09-02-06) 28 Compilation warning are erased in FTF. 29 30 ----------------------------------------------------- 31 5 October 2009, V. Uzhinsky (hadr-prtn-V09-02-05) 32 FTFP with tuned parameters of nuclear de-excitation. 33 ------------------------------------------------------ 34 6 August 2009, V. Uzhinsky (hadr-prtn-V09-02-04) 35 A warning message at FTFmodel compilation connected with Xminus 36 initialisation was erased. 37 ------------------------------------------------------ 38 5 August 2009, V. Uzhinsky (hadr-prtn-V09-02-03) 39 Some warning were erased at FTFModel compilation. 40 ---------------------------------------------------- 41 3 August 2009, V. Uzhinsky (hadr-prtn-V09-02-02) 42 FTF model was complited by the reggeon cascade model. 43 ----------------------------------------------------------- 44 17 July 2009 V. Uzhinsky (hadr-prtn-V09-02-01) 45 Improvement of FTF model: 46 A Status of nuclear nucleon involved in an interaction is introdused. 47 Status: 0 - spectator, 1 - involved nucleon, 2 - absorbed nucleon 48 (G4VSplitableHadron) 49 50 A connection between a participant nucleon and a nuclear nucleon was 51 introsuced in G4InteractionContent. 52 53 All of these allow to improve FTF model for pion obsorption on nuclei. 54 These required a lot of changes in FTF. 55 ------------------------------------------------- 56 10 July 2009 V. Uzhinsky (hadr-prtn-V09-02-00) 57 Introduction the right tag number. 58 ------------------------------------------------- 59 9 July 2009 V. Uzhinsky (hadr-prtn-V09-01-03) 60 Charge-exchange was implemented for pn->np in elastic and inelastic 61 interactions simulated by FTF. Pion absorption by a nucleon was 62 implemented also. tag - hadr-string-diff-V09-02-02 63 64 Bug was fixed in G4VLongitudinalStringDecay.cc at calculation of 65 formation/creation time, c_light was inserted. Due to this string 66 tension parameter was set to the usual value (1 GeV/fm) in 67 G4LundStringFragmentation.cc. tag - had-hadronization-V09-02-02 68 69 New field was added in G4VSplitableHadron class (G4 bool Activation) 70 and corresponding methods to operatite with it. It was needed for an 71 absorption of meson in nuclear collision generated by FTF. 72 tag - had-partonstring-mgt-V08-02-02 73 ----------------------------------------------- 74 16 75 23 June 2008 V. Uzhinsky (hadr-prtn-V09-01-02) 17 76 G4QGSMFragmentation.cc and G4LundStringFragmentation.cc were -
trunk/source/processes/hadronic/models/parton_string/diffraction/History
r962 r1196 1 $Id: History,v 1. 5 2008/12/09 10:43:47vuzhinsk Exp $1 $Id: History,v 1.22 2009/10/29 14:58:52 vuzhinsk Exp $ 2 2 ------------------------------------------------------------------- 3 3 … … 15 15 * Please list in reverse chronological order (last date on top) 16 16 --------------------------------------------------------------- 17 29 October 2009, V. Uzhinsky (hadr-string-diff-V09-02-17) 18 Warning meassage is erased for using slc5 compilator in FTFModel.cc 19 20 25 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 ---------------------------------------------------------- 24 6 October 2009, V. Uzhinsky (hadr-string-diff-V09-02-15) 25 Compilation warning are erased. 26 ---------------------------------------------------------- 27 5 October 2009, V. Uzhinsky (hadr-string-diff-V09-02-14) 28 FTFP with tuned parameters of nuclear de-excitation. 29 30 ----------------------------------------------------------- 31 19 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 ---------------------------------------------------------- 35 17 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 ----------------------------------------- 41 6 August 2009, V. Uzhinsky (hadr-string-diff-V09-02-11) 42 Kaon A interactions are improved in FTFModel. 43 ---------------------------------------------------------- 44 6 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 ---------------------------------------------------------- 48 5 August 2009, V.Uzhinsky (hadr-string-diff-V09-02-09) 49 Some additional warning found by Gabriela were erased. 50 ------------------------------------------------------------ 51 5 August 2009, V. Uzhinsky (hadr-string-diff-V09-02-08) 52 Warning meassages were erased at a compilation of FTF model 53 ----------------------------------------------------------- 54 3 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 ----------------------------------------------------------- 58 31 July 2009, V. Uzhinsky (hadr-string-diff-V09-02-06) 59 Inconsistency in G4FTFModel was erased thank to Gunter. 60 61 ---------------------------------------------------------- 62 31 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 68 17 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 80 10 July V. Uzhinsky (hadr-string-diff-V09-02-03) 81 Bug was fixed in G4FTFModel.cc thank to A. Ribon 82 --------------------------------------------------------- 83 9 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 --------------------------------------------------------- 89 19 December 08 V. Uzhinsky (hadr-string-diff-V09-02-01) 90 --------------------------------------------------------- 91 Version of FTF suted for low energies (2-10 GeV/c) 92 17 93 9 December 08, V. Uzhinsky (hadr-string-diff-V09-01-04) 18 94 - Improvement of delete operators in FTF -
trunk/source/processes/hadronic/models/parton_string/diffraction/include/G4DiffractiveExcitation.hh
r962 r1196 25 25 // 26 26 // 27 // $Id: G4DiffractiveExcitation.hh,v 1. 2 2008/04/25 14:20:13vuzhinsk Exp $27 // $Id: G4DiffractiveExcitation.hh,v 1.4 2009/09/17 18:24:30 vuzhinsk Exp $ 28 28 29 29 #ifndef G4DiffractiveExcitation_h … … 42 42 class G4VSplitableHadron; 43 43 class G4ExcitedString; 44 #include "G4FTFParameters.hh" // Uzhi 19.04.08 44 #include "G4FTFParameters.hh" 45 #include "G4ElasticHNScattering.hh" // Uzhi 3.09.09 45 46 #include "G4ThreeVector.hh" 46 47 … … 50 51 public: 51 52 52 G4DiffractiveExcitation(); // Uzhi53 G4DiffractiveExcitation(); 53 54 virtual ~G4DiffractiveExcitation(); 54 55 55 56 virtual G4bool ExciteParticipants (G4VSplitableHadron *aPartner, 56 57 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, 57 65 G4FTFParameters *theParameters) const; 58 59 virtual G4ExcitedString * String(G4VSplitableHadron * aHadron, G4bool isProjectile) const;60 61 66 private: 62 67 63 68 G4DiffractiveExcitation(const G4DiffractiveExcitation &right); 64 69 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 68 77 const G4DiffractiveExcitation & operator=(const G4DiffractiveExcitation &right); 69 78 int operator==(const G4DiffractiveExcitation &right) const; -
trunk/source/processes/hadronic/models/parton_string/diffraction/include/G4DiffractiveHHScatterer.hh
r962 r1196 25 25 // 26 26 // 27 // $Id: G4DiffractiveHHScatterer.hh,v 1. 4 2008/04/25 14:20:13vuzhinsk Exp $27 // $Id: G4DiffractiveHHScatterer.hh,v 1.7 2009/10/06 10:10:36 vuzhinsk Exp $ 28 28 29 29 #ifndef G4DiffractiveHHScatterer_h … … 44 44 class G4KineticTrack; 45 45 #include "G4KineticTrackVector.hh" 46 #include "G4FTFParameters.hh" // Uzhi 21.04.08 46 #include "G4FTFParameters.hh" 47 #include "G4ExcitedString.hh" 47 48 48 49 class G4DiffractiveHHScatterer … … 52 53 G4DiffractiveHHScatterer(); 53 54 54 G4KineticTrackVector * Scatter(const G4KineticTrack & aTrack, const G4KineticTrack & bTrack);55 // G4KineticTrackVector * Scatter(const G4KineticTrack & aTrack, const G4KineticTrack & bTrack); 55 56 57 virtual void CreateStrings() const; 58 /* 59 (G4VSplitableHadron * aHadron, 60 G4bool isProjectile, 61 G4ExcitedString * FirstString, 62 G4ExcitedString * SecondString, 63 G4FTFParameters *theParameters) const; 64 */ 56 65 private: 57 66 58 67 const G4DiffractiveExcitation * theExcitation; 59 68 G4LundStringFragmentation * theStringFragmentation; 60 G4FTFParameters *theParameters; // Uzhi 21.04.0869 G4FTFParameters *theParameters; 61 70 }; 62 71 -
trunk/source/processes/hadronic/models/parton_string/diffraction/include/G4DiffractiveSplitableHadron.hh
r1007 r1196 25 25 // 26 26 // 27 // $Id: G4DiffractiveSplitableHadron.hh,v 1. 4 2006/06/29 20:54:27 gunterExp $28 // GEANT4 tag $Name: geant4-09-0 2$27 // $Id: G4DiffractiveSplitableHadron.hh,v 1.5 2009/08/03 13:14:19 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 -
trunk/source/processes/hadronic/models/parton_string/diffraction/include/G4ElasticHNScattering.hh
r967 r1196 25 25 // 26 26 // 27 // $Id: G4ElasticHNScattering.hh,v 1. 1 2008/03/31 15:34:01vuzhinsk Exp $27 // $Id: G4ElasticHNScattering.hh,v 1.2 2009/08/03 13:14:19 vuzhinsk Exp $ 28 28 29 29 #ifndef G4ElasticHNScattering_h … … 42 42 class G4VSplitableHadron; 43 43 class G4ExcitedString; 44 #include "G4FTFParameters.hh" // Uzhi 29.03.0844 #include "G4FTFParameters.hh" 45 45 #include "G4ThreeVector.hh" 46 46 … … 50 50 public: 51 51 52 G4ElasticHNScattering(); // Uzhi52 G4ElasticHNScattering(); 53 53 virtual ~G4ElasticHNScattering(); 54 54 … … 61 61 G4ElasticHNScattering(const G4ElasticHNScattering &right); 62 62 63 G4ThreeVector GaussianPt(G4double AveragePt2, G4double maxPtSquare) const; // Uzhi63 G4ThreeVector GaussianPt(G4double AveragePt2, G4double maxPtSquare) const; 64 64 65 65 const G4ElasticHNScattering & operator=(const G4ElasticHNScattering &right); -
trunk/source/processes/hadronic/models/parton_string/diffraction/include/G4FTFModel.hh
r1007 r1196 25 25 // 26 26 // 27 // $Id: G4FTFModel.hh,v 1. 7 2008/04/25 14:20:13vuzhinsk Exp $28 // GEANT4 tag $Name: geant4-09-0 2$27 // $Id: G4FTFModel.hh,v 1.10 2009/10/25 10:50:54 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 // Class Description … … 53 53 class G4ExcitedString; 54 54 55 #include "G4FTFParameters.hh" // Uzhi 29.03.0855 #include "G4FTFParameters.hh" 56 56 #include "G4FTFParticipants.hh" 57 57 … … 64 64 65 65 public: 66 G4FTFModel(); // Uzhi67 G4FTFModel(G4double , G4double , G4double ); // Uzhi66 G4FTFModel(); 67 G4FTFModel(G4double , G4double , G4double ); 68 68 G4FTFModel(G4DiffractiveExcitation * anExcitation); 69 69 G4FTFModel(const G4FTFModel &right); … … 82 82 83 83 private: 84 void ReggeonCascade(); 85 G4bool PutOnMassShell(); 84 86 G4bool ExciteParticipants(); 85 87 G4ExcitedStringVector * BuildStrings(); 88 void GetResidualNucleus(); // 23 Oct. 2009 89 G4ThreeVector GaussianPt(G4double AveragePt2, G4double maxPtSquare) const; 86 90 87 91 private: … … 89 93 G4ReactionProduct theProjectile; 90 94 G4FTFParticipants theParticipants; 95 96 G4Nucleon * TheInvolvedNucleon[250]; 97 G4int NumberOfInvolvedNucleon; 91 98 92 G4FTFParameters *theParameters; // Uzhi 29.03.0899 G4FTFParameters *theParameters; 93 100 G4DiffractiveExcitation * theExcitation; 94 G4ElasticHNScattering * theElastic; // Uzhi 29.03.08101 G4ElasticHNScattering * theElastic; 95 102 103 G4LorentzVector Residual4Momentum; 104 G4double ResidualExcitationEnergy; 96 105 97 106 }; … … 105 114 106 115 #endif 107 108 -
trunk/source/processes/hadronic/models/parton_string/diffraction/include/G4FTFParameters.hh
r1007 r1196 27 27 #define G4FTFParameters_h 1 28 28 // 29 // $Id: G4FTFParameters.hh,v 1. 2 2008/06/13 12:49:23vuzhinsk Exp $30 // GEANT4 tag $Name: geant4-09-0 2$29 // $Id: G4FTFParameters.hh,v 1.7 2009/10/25 10:50:54 vuzhinsk Exp $ 30 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 31 31 // 32 32 #include "G4Proton.hh" … … 60 60 61 61 // --------- 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); 63 67 void SetProjMinNonDiffMass(const G4double aValue); 64 68 void SetProbabilityOfProjDiff(const G4double aValue); 65 69 66 void SetTarMinDiffMass(const G4double aValue); // Uzhi 19.04.0870 void SetTarMinDiffMass(const G4double aValue); 67 71 void SetTarMinNonDiffMass(const G4double aValue); 68 72 void SetProbabilityOfTarDiff(const G4double aValue); … … 70 74 void SetAveragePt2(const G4double aValue); 71 75 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 //-------------------------------------------------------------------- 74 93 // --------- Get geometrical parameteres ----------------------------- 75 94 G4double GetTotalCrossSection(); … … 87 106 88 107 // --------- Get parameters of excitations --------------------------- 89 G4double GetProjMinDiffMass(); // Uzhi 19.04.08 108 G4double GetMagQuarkExchange(); 109 G4double GetSlopeQuarkExchange(); 110 G4double GetDeltaProbAtQuarkExchange(); 111 112 G4double GetProjMinDiffMass(); 90 113 G4double GetProjMinNonDiffMass(); 91 114 G4double GetProbabilityOfProjDiff(); 92 115 93 G4double GetTarMinDiffMass(); // Uzhi 19.04.08116 G4double GetTarMinDiffMass(); 94 117 G4double GetTarMinNonDiffMass(); 95 118 G4double GetProbabilityOfTarDiff(); 96 119 97 98 120 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 99 136 // private: 100 137 … … 115 152 116 153 // --------- Parameters of excitations ------------------------------- 117 G4double ProjMinDiffMass; // Uzhi 19.04.08 154 G4double MagQuarkExchange; 155 G4double SlopeQuarkExchange; 156 G4double DeltaProbAtQuarkExchange; 157 158 G4double ProjMinDiffMass; 118 159 G4double ProjMinNonDiffMass; 119 160 G4double ProbabilityOfProjDiff; … … 124 165 125 166 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 126 182 }; 127 183 … … 162 218 inline void G4FTFParameters::SetAvaragePt2ofElasticScattering(const G4double aPt2) 163 219 { 164 //G4cout<<"Pt2 El "<<aPt2<<" "<<std::sqrt(aPt2)<<G4endl;165 //G4int Uzhi; G4cin>>Uzhi;166 220 AvaragePt2ofElasticScattering = aPt2;} 167 221 168 222 // --------- Set parameters of excitations ---------------------------- 169 inline void G4FTFParameters::SetProjMinDiffMass(const G4double aValue) // Uzhi 19.04.08 223 inline void G4FTFParameters::SetMagQuarkExchange(const G4double aValue) 224 {MagQuarkExchange = aValue;} 225 inline void G4FTFParameters::SetSlopeQuarkExchange(const G4double aValue) 226 {SlopeQuarkExchange = aValue;} 227 inline void G4FTFParameters::SetDeltaProbAtQuarkExchange(const G4double aValue) 228 {DeltaProbAtQuarkExchange = aValue;} 229 230 inline void G4FTFParameters::SetProjMinDiffMass(const G4double aValue) 170 231 {ProjMinDiffMass = aValue*GeV;} 171 232 inline void G4FTFParameters::SetProjMinNonDiffMass(const G4double aValue) … … 174 235 {ProbabilityOfProjDiff = aValue;} 175 236 176 inline void G4FTFParameters::SetTarMinDiffMass(const G4double aValue) // Uzhi 19.04.08237 inline void G4FTFParameters::SetTarMinDiffMass(const G4double aValue) 177 238 {TarMinDiffMass = aValue*GeV;} 178 239 inline void G4FTFParameters::SetTarMinNonDiffMass(const G4double aValue) … … 184 245 {AveragePt2 = aValue*GeV*GeV;} 185 246 247 // --------- Set parameters of a string kink -------------------------------- 248 inline void G4FTFParameters::SetPt2Kink(const G4double aValue) 249 {Pt2kink = aValue;} 250 251 inline 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-------------------- 262 inline void G4FTFParameters::SetCofNuclearDestruction(const G4double aValue) 263 {CofNuclearDestruction = aValue;} 264 inline void G4FTFParameters::SetR2ofNuclearDestruction(const G4double aValue) 265 {R2ofNuclearDestruction = aValue;} 266 267 inline void G4FTFParameters::SetExcitationEnergyPerWoundedNucleon(const G4double aValue) 268 {ExcitationEnergyPerWoundedNucleon = aValue;} 269 270 inline void G4FTFParameters::SetDofNuclearDestruction(const G4double aValue) 271 {DofNuclearDestruction = aValue;} 272 inline void G4FTFParameters::SetPt2ofNuclearDestruction(const G4double aValue) 273 {Pt2ofNuclearDestruction =aValue;} 274 inline void G4FTFParameters::SetMaxPt2ofNuclearDestruction(const G4double aValue) 275 {MaxPt2ofNuclearDestruction = aValue;} 276 186 277 // --------- Get geometrical parameteres ------------------------------ 187 278 inline G4double G4FTFParameters::GetTotalCrossSection() {return FTFXtotal;} … … 211 302 212 303 // --------- Get parameters of excitations --------------------------- 304 inline G4double G4FTFParameters::GetMagQuarkExchange() {return MagQuarkExchange;} 305 inline G4double G4FTFParameters::GetSlopeQuarkExchange() {return SlopeQuarkExchange;} 306 inline G4double G4FTFParameters::GetDeltaProbAtQuarkExchange() 307 {return DeltaProbAtQuarkExchange;} 308 309 213 310 inline G4double G4FTFParameters::GetProjMinDiffMass() {return ProjMinDiffMass;} 214 311 inline G4double G4FTFParameters::GetProjMinNonDiffMass() {return ProjMinNonDiffMass;} … … 221 318 inline G4double G4FTFParameters::GetAveragePt2() {return AveragePt2;} 222 319 320 // --------- Get parameters of a string kink -------------------------- 321 inline G4double G4FTFParameters::GetPt2Kink() {return Pt2kink;} 322 inline std::vector<G4double> 323 G4FTFParameters::GetQuarkProbabilitiesAtGluonSplitUp() 324 {return QuarkProbabilitiesAtGluonSplitUp;} 325 326 // --------- Get parameters of nuclear destruction--------------------- 327 inline G4double G4FTFParameters::GetCofNuclearDestruction(){return CofNuclearDestruction;} 328 inline G4double G4FTFParameters::GetR2ofNuclearDestruction(){return R2ofNuclearDestruction;} 329 330 inline G4double G4FTFParameters::GetExcitationEnergyPerWoundedNucleon() 331 {return ExcitationEnergyPerWoundedNucleon;} 332 333 334 inline G4double G4FTFParameters::GetDofNuclearDestruction() 335 {return DofNuclearDestruction;} 336 inline G4double G4FTFParameters::GetPt2ofNuclearDestruction(){return Pt2ofNuclearDestruction;} 337 inline G4double G4FTFParameters::GetMaxPt2ofNuclearDestruction() 338 {return MaxPt2ofNuclearDestruction;} 223 339 #endif -
trunk/source/processes/hadronic/models/parton_string/diffraction/include/G4FTFParticipants.hh
r1007 r1196 25 25 // 26 26 // 27 // $Id: G4FTFParticipants.hh,v 1. 5 2008/03/31 15:34:01vuzhinsk Exp $28 // GEANT4 tag $Name: geant4-09-0 2$27 // $Id: G4FTFParticipants.hh,v 1.6 2009/08/03 13:14:19 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 … … 41 41 42 42 #include "G4VParticipants.hh" 43 #include "G4FTFParameters.hh" // Uzhi 29.03.0843 #include "G4FTFParameters.hh" 44 44 #include <vector> 45 45 #include "G4Nucleon.hh" … … 62 62 63 63 void GetList(const G4ReactionProduct &thePrimary, 64 G4FTFParameters *theParameters); // Uzhi 29.03.0864 G4FTFParameters *theParameters); 65 65 66 66 void StartLoop(); -
trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4DiffractiveExcitation.cc
r962 r1196 25 25 // 26 26 // 27 // $Id: G4DiffractiveExcitation.cc,v 1. 7 2008/12/18 13:01:58 gunterExp $27 // $Id: G4DiffractiveExcitation.cc,v 1.16 2009/10/06 10:10:36 vuzhinsk Exp $ 28 28 // ------------------------------------------------------------ 29 29 // GEANT 4 class implemetation file … … 48 48 49 49 #include "G4DiffractiveExcitation.hh" 50 #include "G4FTFParameters.hh" 51 #include "G4ElasticHNScattering.hh" 52 50 53 #include "G4LorentzRotation.hh" 54 #include "G4RotationMatrix.hh" 51 55 #include "G4ThreeVector.hh" 52 #include "G4ParticleDefinition.hh" 56 #include "G4ParticleDefinition.hh" 53 57 #include "G4VSplitableHadron.hh" 54 58 #include "G4ExcitedString.hh" 55 #include "G4FTFParameters.hh" // Uzhi 19.04.08 59 #include "G4ParticleTable.hh" 60 #include "G4Neutron.hh" 61 #include "G4ParticleDefinition.hh" 62 56 63 //#include "G4ios.hh" 57 64 //#include "UZHI_diffraction.hh" … … 65 72 ExciteParticipants(G4VSplitableHadron *projectile, 66 73 G4VSplitableHadron *target, 67 G4FTFParameters *theParameters) const 74 G4FTFParameters *theParameters, 75 G4ElasticHNScattering *theElastic) const // Uzhi 03.09.09 68 76 { 69 G4bool PutOnMassShell=0;70 71 77 // -------------------- Projectile parameters ----------------------- 72 73 78 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); 74 94 // G4double M0projectile=projectile->GetDefinition()->GetPDGMass(); // With de-excitation 75 95 G4double M0projectile = Pprojectile.mag(); // Without de-excitation 76 /* 77 G4cout<<"ExciteParticipants-------------------"<<G4endl; 78 G4cout<<"Mom "<<Pprojectile<<" mass "<<M0projectile<<G4endl; 79 */ 96 80 97 if(M0projectile < projectile->GetDefinition()->GetPDGMass()) 81 98 { 82 PutOnMassShell= 1;99 PutOnMassShell=true; 83 100 M0projectile=projectile->GetDefinition()->GetPDGMass(); 84 101 } 85 102 86 103 G4double M0projectile2 = M0projectile * M0projectile; 87 88 G4int PDGcode=projectile->GetDefinition()->GetPDGEncoding();89 G4int absPDGcode=std::abs(PDGcode);90 104 91 105 G4double ProjectileDiffStateMinMass=theParameters->GetProjMinDiffMass(); 92 106 G4double ProjectileNonDiffStateMinMass=theParameters->GetProjMinNonDiffMass(); 93 107 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 98 113 G4LorentzVector Ptarget=target->Get4Momentum(); 114 99 115 G4double M0target = Ptarget.mag(); 100 116 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(); 102 122 103 123 if(M0target < target->GetDefinition()->GetPDGMass()) 104 124 { 105 PutOnMassShell= 1;125 PutOnMassShell=true; 106 126 M0target=target->GetDefinition()->GetPDGMass(); 107 127 } 108 128 109 G4double M0target2 = M0target * M0target; //Ptarget.mag2();110 // for AA-inter.129 G4double M0target2 = M0target * M0target; 130 111 131 G4double TargetDiffStateMinMass=theParameters->GetTarMinDiffMass(); 112 132 G4double TargetNonDiffStateMinMass=theParameters->GetTarMinNonDiffMass(); 113 133 G4double ProbTargetDiffraction=theParameters->GetProbabilityOfTarDiff(); 114 /* 115 G4cout<<TargetDiffStateMinMass<<" "<<TargetNonDiffStateMinMass<<" "<<ProbTargetDiffraction<<G4endl; 116 */ 134 117 135 G4double AveragePt2=theParameters->GetAveragePt2(); 136 137 G4double ProbOfDiffraction=ProbProjectileDiffraction + 138 ProbTargetDiffraction; 139 140 G4double SumMasses=M0projectile+M0target+200.*MeV; 118 141 119 142 // Kinematical properties of the interactions -------------- … … 122 145 G4double S=Psum.mag2(); 123 146 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 else138 {139 ProbProjectileDiffraction=0.;140 }141 // ProbTargetDiffraction /=ProbOfDiffraction;142 143 //G4cout<<"ProbOfDiffraction "<<ProbOfDiffraction<<"ProbProjectileDiffraction "<<ProbProjectileDiffraction<<G4endl; // Vova144 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 155 147 // Transform momenta to cms and then rotate parallel to z axis; 156 157 // G4LorentzVector Psum;158 // Psum=Pprojectile+Ptarget;159 160 148 G4LorentzRotation toCms(-1*Psum.boostVector()); 161 149 162 150 G4LorentzVector Ptmp=toCms*Pprojectile; 163 151 if ( Ptmp.pz() <= 0. ) 164 165 // "String" moving backwards in CMS, abort collision !! 166 //G4cout << " abort Collision!! " << G4endl;167 168 152 { 153 target->SetStatus(2); 154 // "String" moving backwards in CMS, abort collision !! 155 return false; 156 } 169 157 170 158 toCms.rotateZ(-1*Ptmp.phi()); … … 176 164 Ptarget.transform(toCms); 177 165 178 G4double Pt2;179 G4double ProjMassT2, ProjMassT;180 G4double TargMassT2, TargMassT;181 166 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 188 168 G4double SqrtS=std::sqrt(S); 189 190 if(absP DGcode > 1000 && SqrtS < 2200*MeV)191 { return false;}// The model cannot work for169 //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 192 172 // 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 197 178 // Pi+p-interactions 198 179 // at Plab < 1. GeV/c. 199 180 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 202 185 // K+p-interactions 203 186 // at Plab < ??? GeV/c. ??? … … 208 191 209 192 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 211 195 212 196 PZcms = std::sqrt(PZcms2); … … 236 220 237 221 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 ---------------- 336 G4double 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 { 416 G4cout<<"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 { 467 G4cout<<"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 /* 577 else 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; 588 G4cout<<"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 247 629 G4LorentzVector Qmomentum; 248 630 G4double Qminus, Qplus; 249 631 250 632 G4int whilecount=0; 251 // Choose a process633 // Choose a process --------------------------- 252 634 253 635 if(G4UniformRand() < ProbOfDiffraction) … … 255 637 if(G4UniformRand() < ProbProjectileDiffraction) 256 638 { //-------- projectile diffraction --------------- 257 //G4cout<<" Projectile diffraction"<<G4endl; 258 //Uzhi_projectilediffraction++; 639 //G4cout<<"Proj difr"<<G4endl; 259 640 do { 260 641 // Generate pt … … 266 647 { 267 648 Qmomentum=G4LorentzVector(0.,0.,0.,0.); 268 return false; // Ignore this interaction649 target->SetStatus(2); return false; // Ignore this interaction 269 650 }; 270 651 // --------------- Check that the interaction is possible ----------- … … 278 659 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) 279 660 /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 } 290 667 maxPtSquare=PZcms2; 291 668 … … 302 679 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) 303 680 /4./S; 304 //G4cout<<" Pt2 Mpt Mtt Pz2 "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl; 305 306 // if(PZcms2 < 0 ) {PZcms2=0;}; 681 307 682 if(PZcms2 < 0 ) continue; 308 683 PZcms =std::sqrt(PZcms2); … … 310 685 PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms; 311 686 PMinusMax=SqrtS-TargMassT; 312 //G4cout<<" SqrtS P+mim max "<<SqrtS<<" "<<PMinusMin<<" "<<PMinusMax<<G4endl;313 687 314 688 PMinusNew=ChooseP(PMinusMin, PMinusMax); … … 328 702 else 329 703 { // -------------- 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; 332 773 do { 333 774 // Generate pt … … 339 780 { 340 781 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 424 783 }; 425 784 // --------------- Check that the interaction is possible ----------- … … 433 792 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) 434 793 /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 } 445 800 maxPtSquare=PZcms2; 446 801 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0); … … 456 811 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) 457 812 /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 464 814 if(PZcms2 < 0 ) continue; 465 815 PZcms =std::sqrt(PZcms2); … … 469 819 470 820 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; //+++++++++++++++++++++++++++++++++++ Vova476 821 477 822 Qminus=PMinusNew-Pprojectile.minus(); 478 823 479 824 TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms; 480 // TPlusMax=SqrtS-PMinusNew; // Vova481 TPlusMax=SqrtS-ProjMassT; // Vova825 // TPlusMax=SqrtS-PMinusNew; 826 TPlusMax=SqrtS-ProjMassT; 482 827 483 828 TPlusNew=ChooseP(TPlusMin, TPlusMax); 484 485 //G4cout<<"Targ "<<TPlusMin<<" "<<TPlusMax<<" "<<TPlusNew<<G4endl;486 //G4cout<<PMinusNew<<" "<<TPlusNew<<G4endl;487 829 488 830 Qplus=-(TPlusNew-Ptarget.plus()); … … 490 832 Qmomentum.setPz( (Qplus-Qminus)/2 ); 491 833 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 503 835 } while ( 504 836 ((Pprojectile+Qmomentum).mag2() < ProjectileNonDiffStateMinMass2) || //No double Diffraction … … 506 838 } 507 839 508 //G4int Uzhiinp; G4cin>>Uzhiinp; // Vova509 510 840 Pprojectile += Qmomentum; 511 841 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; 524 845 525 846 // Transform back and update SplitableHadron Participant. … … 527 848 Ptarget.transform(toLab); 528 849 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 536 857 projectile->Set4Momentum(Pprojectile); 537 858 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 HARP555 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);582 859 583 860 projectile->IncrementCollisionCount(1); 584 861 target->IncrementCollisionCount(1); 585 862 586 //587 //G4cout<<"Out of Excitation --------------------"<<G4endl;588 //G4int Uzhiinp; G4cin>>Uzhiinp; // Vova589 590 863 return true; 591 864 } 592 865 593 866 // --------------------------------------------------------------------- 594 G4ExcitedString * G4DiffractiveExcitation:: 595 String(G4VSplitableHadron * hadron, G4bool isProjectile) const 867 //G4ExcitedString * G4DiffractiveExcitation:: 868 // String(G4VSplitableHadron * hadron, G4bool isProjectile) const 869 void G4DiffractiveExcitation::CreateStrings(G4VSplitableHadron * hadron, 870 G4bool isProjectile, 871 G4ExcitedString * &FirstString, 872 G4ExcitedString * &SecondString, 873 G4FTFParameters *theParameters) const 596 874 { 597 598 //G4cout<<"G4DiffractiveExcitation::String isProj"<<isProjectile<<G4endl;599 600 875 hadron->SplitUp(); 601 876 G4Parton *start= hadron->GetNextParton(); 602 877 if ( start==NULL) 603 878 { G4cout << " G4FTFModel::String() Error:No start parton found"<< G4endl; 604 return NULL; 879 FirstString=0; SecondString=0; 880 return; 605 881 } 606 882 G4Parton *end = hadron->GetNextParton(); 607 883 if ( end==NULL) 608 884 { G4cout << " G4FTFModel::String() Error:No end parton found"<< G4endl; 609 return NULL; 885 FirstString=0; SecondString=0; 886 return; 610 887 } 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 /* 890 G4cout<<"Create strings had "<<hadron->GetDefinition()->GetParticleName()<<" "<<Phadron<<" "<<Phadron.mag()<<G4endl; 891 G4cout<<"isProjectile "<<isProjectile<<G4endl; 892 G4cout<<"start Q "<<start->GetDefinition()->GetPDGEncoding()<<G4endl; 893 G4cout<<"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 /* 984 G4cout<<"Pstart "<<Pstart<<G4endl; 985 G4cout<<"Pkink "<<Pkink <<G4endl; 986 G4cout<<"Pkink1 "<<PkinkQ1<<G4endl; 987 G4cout<<"Pkink2 "<<PkinkQ2<<G4endl; 988 G4cout<<"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 995 G4LorentzRotation Rotate; 996 if(isProjectile) {Rotate.rotateY(Psi);} 997 else {Rotate.rotateY(pi-Psi);} 998 Rotate.rotateZ(twopi*G4UniformRand()); 999 1000 //G4cout<<"Rotate"<<G4endl; 1001 1002 Pstart*=Rotate; 1003 Pkink*=Rotate; 1004 PkinkQ1*=Rotate; 1005 PkinkQ2*=Rotate; 1006 Pend*=Rotate; 1007 /* 1008 G4cout<<"Pstart "<<Pstart<<G4endl; 1009 G4cout<<"Pkink1 "<<PkinkQ1 <<G4endl; 1010 G4cout<<"Pkink2 "<<PkinkQ2 <<G4endl; 1011 G4cout<<"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; 1031 if(Ksi > QuarkProbabilitiesAtGluonSplitUp[Iq]) QuarkInGluon++;} 1032 1033 //G4cout<<"Gquark "<<QuarkInGluon<<G4endl; 1034 G4Parton * Gquark = new G4Parton(QuarkInGluon); 1035 G4Parton * Ganti_quark = new G4Parton(-QuarkInGluon); 1036 /* 1037 G4cout<<Gquark->GetDefinition()->GetParticleName()<<" "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->GetDefinition()->GetPDGMass()<<G4endl; 1038 G4cout<<Ganti_quark->GetDefinition()->GetParticleName()<<" "<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->GetDefinition()->GetPDGMass()<<G4endl; 1039 */ 1040 1041 //G4int Uzhi; G4cin>>Uzhi; 1042 1043 //------------------------------------------------------------------------------- 1044 /* 1045 DefineMomentumInZ(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 /* 1070 G4cout<<"Pstart "<<start->GetDefinition()->GetPDGEncoding()<<Pstart<<G4endl; 1071 G4cout<<"Pkink1 "<<PkinkQ1 <<G4endl; 1072 G4cout<<"Pkink2 "<<PkinkQ2 <<G4endl; 1073 G4cout<<"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 /* 1091 G4cout<<" Proj Meson end Q"<<G4endl; 1092 G4cout<<"First string ============ "<<G4endl; 1093 G4cout<<"end Q "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl; 1094 G4cout<<"G antiQ"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl; 1095 G4cout<<"Sum P "<<(Ganti_quark->Get4Momentum()+end->Get4Momentum())<<G4endl; 1096 G4cout<<"Secondstring ============ "<<G4endl; 1097 G4cout<<"G Q "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl; 1098 G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl; 1099 1100 G4cout<<"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 /* 1109 G4cout<<" Proj Meson end Qbar"<<G4endl; 1110 G4cout<<"First string ============ "<<G4endl; 1111 G4cout<<"end Q "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl; 1112 G4cout<<"G Q"<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl; 1113 G4cout<<"Sum P "<<(Gquark->Get4Momentum()+end->Get4Momentum())<<G4endl; 1114 G4cout<<"Secondstring ============ "<<G4endl; 1115 G4cout<<"G antQ "<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl; 1116 G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl; 1117 G4cout<<"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 /* 1128 G4cout<<" Targ Meson end Q"<<G4endl; 1129 G4cout<<"First string ============ "<<G4endl; 1130 G4cout<<"G antiQ"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl; 1131 G4cout<<"end Q "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl; 1132 G4cout<<"Sum P "<<(Ganti_quark->Get4Momentum()+end->Get4Momentum())<<G4endl; 1133 G4cout<<"Secondstring ============ "<<G4endl; 1134 G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl; 1135 G4cout<<"G Q "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl; 1136 G4cout<<"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 /* 1145 G4cout<<" Targ Meson end Qbar"<<G4endl; 1146 G4cout<<"First string ============ "<<G4endl; 1147 G4cout<<"G Q"<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl; 1148 G4cout<<"end Q "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl; 1149 G4cout<<"Sum P "<<(Gquark->Get4Momentum()+end->Get4Momentum())<<G4endl; 1150 G4cout<<"Secondstring ============ "<<G4endl; 1151 G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl; 1152 G4cout<<"G antQ "<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl; 1153 G4cout<<"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 /* 1169 G4cout<<" Proj baryon end QQ"<<G4endl; 1170 G4cout<<"First string ============ "<<G4endl; 1171 G4cout<<"end OQ "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl; 1172 G4cout<<"G Q"<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl; 1173 G4cout<<"Sum P "<<(Gquark->Get4Momentum()+end->Get4Momentum())<<G4endl; 1174 G4cout<<"Secondstring ============ "<<G4endl; 1175 G4cout<<"G Qbar"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl; 1176 G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl; 1177 G4cout<<"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 /* 1186 G4cout<<" Proj baryon end Q"<<G4endl; 1187 G4cout<<"First string ============ "<<G4endl; 1188 G4cout<<"end OQ "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl; 1189 G4cout<<"G antQ"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl; 1190 G4cout<<"Sum P "<<(Ganti_quark->Get4Momentum()+end->Get4Momentum())<<G4endl; 1191 G4cout<<"Secondstring ============ "<<G4endl; 1192 G4cout<<"G Q "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl; 1193 G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl; 1194 G4cout<<"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 /* 1208 G4cout<<" Targ baryon end QQ"<<G4endl; 1209 G4cout<<"First string ============ "<<G4endl; 1210 G4cout<<"G Q "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl; 1211 G4cout<<"end OQ "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl; 1212 G4cout<<"Sum P "<<(Gquark->Get4Momentum()+end->Get4Momentum())<<G4endl; 1213 G4cout<<"Secondstring ============ "<<G4endl; 1214 G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl; 1215 G4cout<<"G Qbar"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl; 1216 G4cout<<"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 /* 1225 G4cout<<" Targ baryon end Q"<<G4endl; 1226 G4cout<<"First string ============ "<<G4endl; 1227 G4cout<<"G Qbar"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl; 1228 G4cout<<"end O "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl; 1229 G4cout<<"Sum P "<<(Ganti_quark->Get4Momentum()+end->Get4Momentum())<<G4endl; 1230 G4cout<<"Secondstring ============ "<<G4endl; 1231 G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl; 1232 G4cout<<"G Q "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl; 1233 G4cout<<"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()); 625 1259 626 1260 // momenta of string ends 627 1261 // 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 713 1304 #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; 725 1316 #endif 726 1317 727 return string; 1318 return; 1319 728 1320 } 729 1321 … … 732 1324 733 1325 // --------------------------------------------------------------------- 734 G4double G4DiffractiveExcitation::ChooseP(G4double Pmin, G4double Pmax) const // Uzhi1326 G4double G4DiffractiveExcitation::ChooseP(G4double Pmin, G4double Pmax) const 735 1327 { 736 1328 // choose an x between Xmin and Xmax with P(x) ~ 1/x 737 1329 // to be improved... 738 1330 739 G4double range=Pmax-Pmin; // Uzhi1331 G4double range=Pmax-Pmin; 740 1332 741 1333 if ( Pmin <= 0. || range <=0. ) … … 745 1337 } 746 1338 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()); 757 1340 return P; 758 1341 } … … 760 1343 // --------------------------------------------------------------------- 761 1344 G4ThreeVector G4DiffractiveExcitation::GaussianPt(G4double AveragePt2, 762 G4double maxPtSquare) const // Uzhi1345 G4double maxPtSquare) const 763 1346 { // @@ this method is used in FTFModel as well. Should go somewhere common! 764 1347 765 1348 G4double Pt2; 766 /* // Uzhi767 do {768 pt2=widthSquare * std::log( G4UniformRand() );769 } while ( pt2 > maxPtSquare);770 */ // Uzhi771 772 1349 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * 773 (std::exp(-maxPtSquare/AveragePt2)-1.)); // Uzhi1350 (std::exp(-maxPtSquare/AveragePt2)-1.)); 774 1351 775 1352 G4double Pt=std::sqrt(Pt2); 776 777 1353 G4double phi=G4UniformRand() * twopi; 778 779 1354 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.); 780 1355 } 781 1356 1357 // --------------------------------------------------------------------- 1358 G4double 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 // --------------------------------------------------------------------- 1369 void 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 // --------------------------------------------------------------------- 1383 void 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 // --------------------------------------------------------------------- 1392 G4int 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 } 782 1417 // --------------------------------------------------------------------- 783 1418 G4DiffractiveExcitation::G4DiffractiveExcitation(const G4DiffractiveExcitation &) -
trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4DiffractiveHHScatterer.cc
r962 r1196 38 38 {} 39 39 40 // ------------------------------------------------------------------- 40 void G4DiffractiveHHScatterer::CreateStrings() 41 /* 42 G4VSplitableHadron * aHadron, 43 G4bool isProjectile, 44 G4ExcitedString * FirstString, 45 G4ExcitedString * SecondString, 46 G4FTFParameters *theParameters) 47 */ 48 const 49 50 {} 51 52 /* ------------------------------------------------------------------- 41 53 G4KineticTrackVector * G4DiffractiveHHScatterer:: 42 54 Scatter(const G4KineticTrack & aTrack, const G4KineticTrack & bTrack) … … 76 88 return result; 77 89 } 90 */ -
trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4DiffractiveSplitableHadron.cc
r1007 r1196 25 25 // 26 26 // 27 // $Id: G4DiffractiveSplitableHadron.cc,v 1. 7 2008/03/31 15:34:01vuzhinsk Exp $28 // GEANT4 tag $Name: geant4-09-0 2$27 // $Id: G4DiffractiveSplitableHadron.cc,v 1.8 2009/07/31 11:03:00 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 … … 87 87 88 88 G4int PDGcode=GetDefinition()->GetPDGEncoding(); 89 89 90 G4int stringStart, stringEnd; 90 91 ChooseStringEnds(PDGcode, &stringStart,&stringEnd); 91 92 92 93 Parton[0] = new G4Parton(stringStart); 93 94 Parton[1] = new G4Parton(stringEnd); 94 95 PartonIndex=-1; 95 96 96 } 97 97 -
trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4ElasticHNScattering.cc
r968 r1196 25 25 // 26 26 // 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 $ 28 28 // ------------------------------------------------------------ 29 29 // GEANT 4 class implemetation file … … 46 46 #include "G4VSplitableHadron.hh" 47 47 #include "G4ExcitedString.hh" 48 #include "G4FTFParameters.hh" // Uzhi 29.03.0848 #include "G4FTFParameters.hh" 49 49 //#include "G4ios.hh" 50 50 … … 59 59 { 60 60 //G4cout<<"G4ElasticHNScattering::ElasticScattering"<<G4endl; 61 61 // -------------------- Projectile parameters ----------------------------------- 62 62 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); 66 75 67 76 G4double M0projectile = Pprojectile.mag(); 68 69 77 if(M0projectile < projectile->GetDefinition()->GetPDGMass()) 70 71 PutOnMassShell= 1;78 { 79 PutOnMassShell=true; 72 80 M0projectile=projectile->GetDefinition()->GetPDGMass(); 73 81 } 74 82 75 83 G4double Mprojectile2 = M0projectile * M0projectile; 76 84 77 // G4double AveragePt2=theParameters->GetSlope(); // Uzhi ???78 // AveragePt2 = AveragePt2 * GeV*GeV;79 80 85 G4double AveragePt2=theParameters->GetAvaragePt2ofElasticScattering(); 81 86 82 87 // -------------------- Target parameters ---------------------------------------------- 88 // G4int TargetPDGcode=target->GetDefinition()->GetPDGEncoding(); 89 83 90 G4LorentzVector Ptarget=target->Get4Momentum(); 84 91 //G4cout<<"Elastic T "<<Ptarget<<G4endl; 92 // G4double TargetRapidity = Ptarget.rapidity(); 85 93 G4double M0target = Ptarget.mag(); 86 //G4cout<<" Mp Mt Pt2 "<<M0projectile<<" "<<M0target<<" "<<AveragePt2/GeV/GeV<<G4endl;87 94 88 95 if(M0target < target->GetDefinition()->GetPDGMass()) 89 90 PutOnMassShell= 1;96 { 97 PutOnMassShell=true; 91 98 M0target=target->GetDefinition()->GetPDGMass(); 92 99 } 93 100 94 G4double Mtarget2 = M0target * M0target; //Ptarget.mag2();95 // for AA-inter. 101 G4double Mtarget2 = M0target * M0target; 102 96 103 // Transform momenta to cms and then rotate parallel to z axis; 97 104 … … 103 110 G4LorentzVector Ptmp=toCms*Pprojectile; 104 111 105 if ( Ptmp.pz() <= 0. ) // Uzhi ???112 if ( Ptmp.pz() <= 0. ) 106 113 { 107 114 // "String" moving backwards in CMS, abort collision !! 108 115 //G4cout << " abort Collision!! " << G4endl; 116 target->SetStatus(2); 109 117 return false; 110 118 } … … 118 126 Ptarget.transform(toCms); 119 127 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 -------------------------- 121 197 G4double Pt2; 122 198 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 160 201 G4LorentzVector Qmomentum; 161 202 … … 163 204 164 205 Pt2=G4ThreeVector(Qmomentum.vect()).mag2(); 165 //G4cout<<"Pt2 GeV^2 "<<(Pt2)/GeV/GeV<<G4endl;166 206 167 207 ProjMassT2=Mprojectile2+Pt2; … … 175 215 2.*S*ProjMassT2-2.*S*TargMassT2- 176 216 2.*ProjMassT2*TargMassT2)/4./S; 177 if(PZcms2 < 0 ) {PZcms2=0;}; 217 if(PZcms2 < 0 ) {PZcms2=0;};// to avoid the exactness problem 178 218 PZcms =std::sqrt(PZcms2); 179 219 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); 188 222 189 223 Pprojectile += Qmomentum; 190 224 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;197 225 198 226 // Transform back and update SplitableHadron Participant. 199 227 Pprojectile.transform(toLab); 200 228 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 // ------------------------------------------------------ 232 241 233 242 projectile->Set4Momentum(Pprojectile); -
trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4FTFModel.cc
r1007 r1196 25 25 // 26 26 // 27 // $Id: G4FTFModel.cc,v 1. 13 2008/12/09 10:40:52vuzhinsk Exp $28 // GEANT4 tag $Name: geant4-09-0 2$27 // $Id: G4FTFModel.cc,v 1.28 2009/10/29 14:55:33 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 … … 38 38 39 39 #include "G4FTFModel.hh" 40 #include "G4FTFParameters.hh" // Uzhi 29.03.0840 #include "G4FTFParameters.hh" 41 41 #include "G4FTFParticipants.hh" 42 #include "G4DiffractiveSplitableHadron.hh" 42 43 #include "G4InteractionContent.hh" 43 44 #include "G4LorentzRotation.hh" 44 45 #include "G4ParticleDefinition.hh" 46 #include "G4ParticleTable.hh" 45 47 #include "G4ios.hh" 46 #include <utility> // Uzhi 29.03.08 48 #include <utility> 49 #include "G4IonTable.hh" 47 50 48 51 // Class G4FTFModel 49 52 50 53 G4FTFModel::G4FTFModel():theExcitation(new G4DiffractiveExcitation()), 51 theElastic(new G4ElasticHNScattering()) // Uzhi 29.03.0854 theElastic(new G4ElasticHNScattering()) 52 55 { 53 56 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 } 70 59 71 60 72 61 G4FTFModel::~G4FTFModel() 73 62 { 74 if( theParameters != 0 ) delete theParameters; // Uzhi 5.12.0863 if( theParameters != 0 ) delete theParameters; 75 64 // Because FTF model can be called for various particles 76 65 // 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 } 82 70 83 71 const G4FTFModel & G4FTFModel::operator=(const G4FTFModel &) … … 87 75 } 88 76 89 90 77 int G4FTFModel::operator==(const G4FTFModel &right) const 91 78 { … … 102 89 { 103 90 theProjectile = aProjectile; 104 //G4cout<<"G4FTFModel::Init "<<aNucleus.GetN()<<" "<<aNucleus.GetZ()<<G4endl;105 91 theParticipants.Init(aNucleus.GetN(),aNucleus.GetZ()); 106 // Uzhi N-mass number Z-charge ------------------------- Uzhi 29.03.0892 // ----------- N-mass number Z-charge ------------------------- 107 93 108 94 // --- cms energy … … 112 98 2*theProjectile.GetTotalEnergy()*G4Proton::Proton()->GetPDGMass(); 113 99 /* 100 G4cout<<"----------------------------------------"<<G4endl; 114 101 G4cout << " primary Total E (GeV): " << theProjectile.GetTotalEnergy()/GeV << G4endl; 115 102 G4cout << " primary Mass (GeV): " << theProjectile.GetMass() /GeV << G4endl; 103 G4cout << " primary 3Mom " << theProjectile.GetMomentum() << G4endl; 104 G4cout << " primary space position " << theProjectile.GetPositionInNucleus() << G4endl; 116 105 G4cout << "cms std::sqrt(s) (GeV) = " << std::sqrt(s) / GeV << G4endl; 117 106 */ 118 107 119 if( theParameters != 0 ) delete theParameters; // Uzhi 9.12.08108 if( theParameters != 0 ) delete theParameters; 120 109 theParameters = new G4FTFParameters(theProjectile.GetDefinition(), 121 110 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 124 115 } 125 116 … … 127 118 G4ExcitedStringVector * G4FTFModel::GetStrings() 128 119 { 129 //G4cout<<"theParticipants.GetList"<<G4endl; 120 //G4int Uzhi; G4cin>>Uzhi; 121 122 //G4cout<<"GetList"<<G4endl; 130 123 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; 134 131 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 } 141 141 return theStrings; 142 142 } 143 143 144 144 // ------------------------------------------------------------ 145 struct DeleteVSplitableHadron { void operator()(G4VSplitableHadron * aH){delete aH;} }; 145 struct DeleteVSplitableHadron { void operator()(G4VSplitableHadron * aH){ delete aH;} }; 146 147 // ------------------------------------------------------------ 148 void 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 // ------------------------------------------------------------ 239 G4bool 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 } 146 557 147 558 // ------------------------------------------------------------ 148 559 G4bool G4FTFModel::ExciteParticipants() 149 560 { 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; 154 565 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 160 576 while (theParticipants.Next()) 161 577 { 162 578 const G4InteractionContent & collision=theParticipants.GetInteraction(); 163 / *164 counter++;165 G4cout<<" Inter #"<<counter<<G4endl;166 */579 // 580 //counter++; 581 //G4cout<<" Int num "<<counter<<G4endl; 582 // 167 583 G4VSplitableHadron * projectile=collision.GetProjectile(); 168 584 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 ---------------------------- 172 587 if(G4UniformRand()< theParameters->GetProbabilityOfElasticScatt()) 173 { 588 { // Elastic scattering ------------------------- 174 589 //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 }; 176 606 } 177 607 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 }; 181 627 } 182 // if(!Successfull) 183 // // Uzhi 29.03.08184 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; 188 634 // give up, clean up 189 190 std::vector<G4VSplitableHadron *> targets; 191 192 193 194 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(); 195 641 // do not allow for duplicates ... 196 197 198 199 200 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(), 201 647 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 212 676 return true; 213 677 } … … 217 681 // Loop over all collisions; find all primaries, and all target ( targets may 218 682 // be duplicate in the List ( to unique G4VSplitableHadrons) 683 684 //G4cout<<"In build string"<<G4endl; 219 685 220 686 G4ExcitedStringVector * strings; … … 223 689 std::vector<G4VSplitableHadron *> primaries; 224 690 std::vector<G4VSplitableHadron *> targets; 691 std::vector<G4Nucleon *> TargetNucleons; // Uzhi 16.07.09 225 692 693 G4ExcitedString * FirstString(0); // If there will be a kink, 694 G4ExcitedString * SecondString(0); // two strings will be prodused. 695 226 696 theParticipants.StartLoop(); // restart a loop 697 //G4int InterCount(0); // Uzhi 227 698 while ( theParticipants.Next() ) 228 699 { 229 700 const G4InteractionContent & interaction=theParticipants.GetInteraction(); 230 701 // do not allow for duplicates ... 702 231 703 if ( primaries.end() == std::find(primaries.begin(), primaries.end(), 232 704 interaction.GetProjectile()) ) … … 236 708 interaction.GetTarget()) ) 237 709 targets.push_back(interaction.GetTarget()); 710 711 if ( TargetNucleons.end() == std::find(TargetNucleons.begin(), 712 TargetNucleons.end(), 713 interaction.GetTargetNucleon()) ) 714 TargetNucleons.push_back(interaction.GetTargetNucleon()); 715 //InterCount++; 238 716 } 239 717 240 718 241 // G4cout << "BuildStrings prim/targ " << primaries.size() << " , " << 242 // targets.size() << G4endl; 719 //G4cout << "BuildStrings prim/targ/woundN " << primaries.size() << " , " <<targets.size() <<", "<<TargetNucleons.size()<< G4endl; 243 720 244 721 unsigned int ahadron; 245 722 // Only for hA-interactions Uzhi ------------------------------------- 723 //G4int StringN(0); 724 //G4cout<<"Proj strings -----------------------"<<G4endl; 246 725 for ( ahadron=0; ahadron < primaries.size() ; ahadron++) 247 726 { 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; 252 741 } 253 742 //G4cout<<"Targ strings ------------------------------"<<G4endl; 254 743 for ( ahadron=0; ahadron < targets.size() ; ahadron++) 255 744 { 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 } 260 765 } 766 767 //G4cout<<"Proj + Targ string "<<strings->size()<<G4endl; 768 //G4cout<<"Involv strings NumberOfInvolvedNucleon "<<NumberOfInvolvedNucleon<<G4endl; 769 for (G4int ahadron=0; ahadron < NumberOfInvolvedNucleon ; ahadron++) 770 { 771 /* 772 G4cout<<"ahadron "<<ahadron<<" Status "<< 773 TheInvolvedNucleon[ahadron]->GetSplitableHadron()->GetStatus()<< 774 TheInvolvedNucleon[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 /* 795 else 796 { 797 G4cout<<"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; 261 808 262 809 std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron()); 263 810 primaries.clear(); 811 264 812 std::for_each(targets.begin(), targets.end(), DeleteVSplitableHadron()); 265 813 targets.clear(); 814 815 return strings; 816 } 817 // ------------------------------------------------------------ 818 void 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 // ------------------------------------------------------------ 836 G4ThreeVector G4FTFModel::GaussianPt(G4double AveragePt2, G4double maxPtSquare) const 837 { // @@ this method is used in FTFModel as well. Should go somewhere common! 266 838 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 //******************************************************************** 3 2 // * License and Disclaimer * 4 3 // * * … … 25 24 // 26 25 // 27 // $Id: G4FTFParameters.cc,v 1. 4 2008/12/18 13:02:00 gunterExp $28 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4FTFParameters.cc,v 1.11 2009/10/25 10:50:54 vuzhinsk Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 28 // 30 29 31 30 #include "G4FTFParameters.hh" 31 32 #include "G4ios.hh" 33 #include <utility> // Uzhi 29.03.08 32 34 33 35 G4FTFParameters::G4FTFParameters() … … 43 45 G4double theZ, 44 46 G4double s) 45 47 { 46 48 G4int PDGcode = particle->GetPDGEncoding(); 47 49 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); 50 56 51 57 G4double LogPlab = std::log( Plab ); 52 58 G4double sqrLogPlab = LogPlab * LogPlab; 53 54 //G4cout<<"G4FTFParameters Plab "<<Plab<<G4endl;55 59 56 60 G4int NumberOfTargetProtons = (G4int) theZ; … … 144 148 NumberOfTargetNeutrons * XelKN ) / NumberOfTargetNucleons; 145 149 } 146 else if( PDGcode == 311 ) //------Projectile is KaonZero ------150 else if((PDGcode == 311) || (PDGcode == 130) || (PDGcode == 310))//Projectile is KaonZero 147 151 { 148 152 G4double XtotKP =( 18.1 + 0. *std::pow(Plab, 0. ) + 0.26 *sqrLogPlab - 1.0 *LogPlab + //K+ … … 181 185 SetInelasticCrossSection(Xtotal-Xelastic); 182 186 183 // // Interactions with elastic an sinelastic collisions187 // // Interactions with elastic and inelastic collisions 184 188 SetProbabilityOfElasticScatt(Xtotal, Xelastic); 185 189 SetRadiusOfHNinteractions2(Xtotal/pi/10.); … … 190 194 */ //======================================================= 191 195 192 //G4cout<<" Rnn "<<Xtotal/pi/10.<<" "<<Xtotal/pi/10.*fermi*fermi<<G4endl;193 //G4cout<<"G4FTFParameters Xt Xel MeV "<<Xtotal<<" "<<Xelastic<<" "<<GeV<<G4endl;194 195 196 //----------------------------------------------------------------------------------- 196 197 SetSlope( Xtotal*Xtotal/16./pi/Xelastic/0.3894 ); // Slope parameter of elastic scattering 197 198 // (GeV/c)^(-2)) 198 //G4cout<<"G4FTFParameters Slope "<<GetSlope()<<G4endl;199 199 //----------------------------------------------------------------------------------- 200 200 SetGamma0( GetSlope()*Xtotal/10./2./pi ); … … 208 208 if( absPDGcode > 1000 ) //------Projectile is baryon -------- 209 209 { 210 SetMagQuarkExchange(3.4); //3.8); 211 SetSlopeQuarkExchange(1.2); 212 SetDeltaProbAtQuarkExchange(0.1); //(0.1*4.); 213 210 214 SetProjMinDiffMass(1.1); // GeV 211 215 SetProjMinNonDiffMass(1.1); // GeV 212 SetProbabilityOfProjDiff(0. 95*std::pow(s/GeV/GeV,-0.35)); // 40/32 X-dif/X-inel216 SetProbabilityOfProjDiff(0.76*std::pow(s/GeV/GeV,-0.35)); 213 217 214 218 SetTarMinDiffMass(1.1); // GeV 215 219 SetTarMinNonDiffMass(1.1); // GeV 216 SetProbabilityOfTarDiff(0. 95*std::pow(s/GeV/GeV,-0.35)); // 40/32 X-dif/X-inel220 SetProbabilityOfTarDiff(0.76*std::pow(s/GeV/GeV,-0.35)); 217 221 218 222 SetAveragePt2(0.3); // GeV^2 … … 220 224 else if( absPDGcode == 211 || PDGcode == 111) //------Projectile is Pion ----------- 221 225 { 226 SetMagQuarkExchange(120.); // 210. 227 SetSlopeQuarkExchange(2.0); 228 SetDeltaProbAtQuarkExchange(0.6); 229 222 230 SetProjMinDiffMass(0.5); // GeV 223 231 SetProjMinNonDiffMass(0.3); // GeV 224 SetProbabilityOfProjDiff(0. 62*std::pow(s/GeV/GeV,-0.51)); // 40/32 X-dif/X-inel225 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; 226 234 SetTarMinDiffMass(1.1); // GeV 227 235 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 229 239 230 240 /* … … 236 246 SetAveragePt2(0.3); // GeV^2 237 247 } 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 239 250 { 251 // Must be corrected, taken from PiN 252 SetMagQuarkExchange(120.); 253 SetSlopeQuarkExchange(2.0); 254 SetDeltaProbAtQuarkExchange(0.6); 255 240 256 SetProjMinDiffMass(0.7); // GeV 1.1 241 257 SetProjMinNonDiffMass(0.7); // GeV … … 251 267 //------Nucleon assumed 252 268 { 269 SetMagQuarkExchange(3.5); 270 SetSlopeQuarkExchange(1.0); 271 SetDeltaProbAtQuarkExchange(0.1); 272 253 273 SetProjMinDiffMass((particle->GetPDGMass()+160.*MeV)/GeV); 254 274 SetProjMinNonDiffMass((particle->GetPDGMass()+160.*MeV)/GeV); … … 260 280 261 281 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 } 268 311 //********************************************************************************************** -
trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4FTFParticipants.cc
r1007 r1196 25 25 // 26 26 // 27 // $Id: G4FTFParticipants.cc,v 1. 9 2008/06/13 12:49:23vuzhinsk Exp $28 // GEANT4 tag $Name: geant4-09-0 2$27 // $Id: G4FTFParticipants.cc,v 1.15 2009/10/25 10:50:54 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 // ------------------------------------------------------------ … … 87 87 xyradius =theNucleus->GetOuterRadius() + deltaxy; // Impact parameter sampling 88 88 // radius 89 G4bool nucleusNeedsShift = true; 89 // G4bool nucleusNeedsShift = true; // Uzhi 20 July 2009 90 90 91 91 while ( theInteractions.size() == 0 ) … … 96 96 G4double impactY = theImpactParameter.second; 97 97 98 G4ThreeVector thePosition(impactX, impactY, -DBL_MAX); //Uzhi 24 July 09 99 primarySplitable->SetPosition(thePosition); //Uzhi 24 July 09 100 98 101 theNucleus->StartLoop(); 99 102 G4Nucleon * nucleon; 100 103 //G4int InterNumber=0; // Uzhi 104 G4int NucleonNumber(0); // Uzhi 101 105 //while ( (nucleon=theNucleus->GetNextNucleon())&& (InterNumber < 1) ) // Uzhi 102 106 while ( (nucleon=theNucleus->GetNextNucleon()) ) // Uzhi 103 107 { 108 //G4cout<<"Nucl# "<<NucleonNumber<<G4endl; // Vova 104 109 G4double impact2= sqr(impactX - nucleon->GetPosition().x()) + 105 110 sqr(impactY - nucleon->GetPosition().y()); … … 110 115 { 111 116 //InterNumber++; 117 /* // Uzhi 20 July 2009 112 118 if ( nucleusNeedsShift ) 113 119 { // on the first hit, shift nucleus … … 117 123 impactY=0; 118 124 } 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() ) 121 131 { 132 //G4cout<<"Part "<<nucleon->Get4Momentum()<<G4endl; 122 133 targetSplitable= new G4DiffractiveSplitableHadron(*nucleon); 123 134 nucleon->Hit(targetSplitable); 135 nucleon->SetBindingEnergy(3.*nucleon->GetBindingEnergy()); // Uzhi 5.10.09 136 targetSplitable->SetStatus(1); // It takes part in the interaction 124 137 } 125 138 G4InteractionContent * aInteraction = 126 139 new G4InteractionContent(primarySplitable); 127 140 aInteraction->SetTarget(targetSplitable); 141 aInteraction->SetTargetNucleon(nucleon); // Uzhi 16.07.09 128 142 theInteractions.push_back(aInteraction); 129 143 } 144 NucleonNumber++; // Uzhi 130 145 } 131 132 // G4cout << "Number of Hit nucleons " << theInteractions.size() // entries()146 // // Uzhi 147 // G4cout << "Number of Hit nucleons " << theInteractions.size()<<G4endl; // entries() 133 148 // << "\t" << impactX/fermi << "\t"<<impactY/fermi 134 149 // << "\t" << std::sqrt(sqr(impactX)+sqr(impactY))/fermi <<G4endl; 135 150 // // Uzhi 136 151 } 137 152 -
trunk/source/processes/hadronic/models/parton_string/hadronization/History
r1055 r1196 1 $Id: History,v 1. 7 2009/05/22 16:36:52 gunterExp $1 $Id: History,v 1.10 2009/07/17 12:30:45 vuzhinsk Exp $ 2 2 ------------------------------------------------------------------- 3 3 … … 15 15 * Please list in reverse chronological order (last date on top) 16 16 --------------------------------------------------------------- 17 17-July-2009 V. Uzhinsky (had-hadronization-V09-02-03) 18 Some small optimization are done. 19 --------------------------------------------------------------------- 20 9-July-2009 V. Uzhinsky (had-hadronization-V09-02-02) 21 --------------------------------------------------------------------- 22 - Bug was fixed in G4VLongitudinalStringDecay.cc at calculation of 23 formation/creation time, c_light was inserted. 24 - String tension parameter was set to the usual value (1 GeV/fm) in 25 G4LundStringFragmentation.cc 17 26 27 -------------------------------------------------------------- 18 28 22-May-2009 Gunter Folger (had-hadronization-V09-02-01) 19 29 - remove temporary workaround - fix is now in QGS … … 23 33 implemented; affected .hh and .cc for new optional argument of max Pt 24 34 to SampleQuarkPt() 35 25 36 26 37 18-May-2009 Gunter Folger (had-hadronization-V09-02-00) -
trunk/source/processes/hadronic/models/parton_string/hadronization/include/G4ExcitedStringDecay.hh
r1007 r1196 25 25 // 26 26 // 27 // $Id: G4ExcitedStringDecay.hh,v 1. 7 2007/05/03 22:06:17 gunterExp $28 // GEANT4 tag $Name: geant4-09-0 2$27 // $Id: G4ExcitedStringDecay.hh,v 1.10 2009/10/05 12:52:48 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 #ifndef G4ExcitedStringDecay_h … … 81 81 G4LorentzVector KTsecondaries; 82 82 G4bool NeedEnergyCorrector=false; 83 83 //G4cout<<"Number of strings "<<theStrings->size()<<G4endl; // Vova 84 84 for ( unsigned int astring=0; astring < theStrings->size(); astring++) 85 85 { 86 //G4cout<<"String# "<<astring<<" 4mom "<<theStrings->operator[](astring)->Get4Momentum()<<G4endl; // Vova 86 87 KTsum+= theStrings->operator[](astring)->Get4Momentum(); 88 //G4cout<<"KTsum "<<KTsum<<G4endl; 87 89 if( !(KTsum.e()<1) && !(KTsum.e()>-1) ) 88 90 { … … 105 107 continue; 106 108 } 107 109 110 //G4cout<<" prod had "<<generatedKineticTracks->size()<<G4endl; // Vova 108 111 G4LorentzVector KTsum1; 109 112 for ( unsigned int aTrack=0; aTrack<generatedKineticTracks->size();aTrack++) -
trunk/source/processes/hadronic/models/parton_string/hadronization/include/G4FragmentingString.hh
r1007 r1196 26 26 // 27 27 // $Id: G4FragmentingString.hh,v 1.4 2007/12/20 15:38:06 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-0 2$28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 -
trunk/source/processes/hadronic/models/parton_string/hadronization/include/G4HadronBuilder.hh
r1055 r1196 26 26 // 27 27 // $Id: G4HadronBuilder.hh,v 1.5 2009/05/18 09:43:40 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 // ----------------------------------------------------------------------------- -
trunk/source/processes/hadronic/models/parton_string/hadronization/include/G4LundStringFragmentation.hh
r1007 r1196 26 26 // 27 27 // $Id: G4LundStringFragmentation.hh,v 1.6 2008/04/25 14:20:14 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-0 2$ Maxim Komogorov28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ Maxim Komogorov 29 29 // 30 30 // ----------------------------------------------------------------------------- -
trunk/source/processes/hadronic/models/parton_string/hadronization/include/G4QGSMFragmentation.hh
r1007 r1196 26 26 // 27 27 // $Id: G4QGSMFragmentation.hh,v 1.5 2007/12/20 15:38:07 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-0 2$28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 // ----------------------------------------------------------------------------- -
trunk/source/processes/hadronic/models/parton_string/hadronization/include/G4VKinkyStringDecay.hh
r1007 r1196 26 26 // 27 27 // $Id: G4VKinkyStringDecay.hh,v 1.3 2006/06/29 20:54:53 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-0 2$28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // Maxim Komogorov 30 30 // -
trunk/source/processes/hadronic/models/parton_string/hadronization/include/G4VLongitudinalStringDecay.hh
r1055 r1196 26 26 // 27 27 // $Id: G4VLongitudinalStringDecay.hh,v 1.7 2009/05/22 16:35:47 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // Maxim Komogorov 30 30 // -
trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4ExcitedStringDecay.cc
r1055 r1196 148 148 Beta = TotalCollisionMom.boostVector(); 149 149 Output->Boost(Beta); 150 /* // Uzhi 151 G4cout<<"Number of produced hadrons // correct E and P "<<Output->size()<<G4endl; // Uzhi 152 for(cHadron = 0; cHadron < Output->size(); cHadron++) 153 { 154 G4cout<<cHadron<<" "<<Output->operator[](cHadron)->Get4Momentum()<<" "<<Output->operator[](cHadron)->Get4Momentum().mag()<<Output->operator[](cHadron)->GetDefinition()->GetParticleName()<<G4endl; 155 } 156 */ // Uzhi 150 157 return success; 151 158 } -
trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4HadronBuilder.cc
r1055 r1196 26 26 // 27 27 // $Id: G4HadronBuilder.cc,v 1.10 2009/05/22 16:34:31 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 // ----------------------------------------------------------------------------- -
trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4LundStringFragmentation.cc
r1055 r1196 25 25 // 26 26 // 27 // $Id: G4LundStringFragmentation.cc,v 1.1 4 2009/05/22 16:36:46 gunterExp $28 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $ 1.827 // $Id: G4LundStringFragmentation.cc,v 1.18 2009/10/05 12:52:48 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 1.8 29 29 // 30 30 // ----------------------------------------------------------------------------- … … 58 58 SmoothParam = 0.2; 59 59 60 SetStringTensionParameter(0.25); // Uzhi 20 June 08 60 // SetStringTensionParameter(0.25); // Uzhi 20 June 08 61 SetStringTensionParameter(1.); // Uzhi 20 June 09 61 62 } 62 63 … … 196 197 LeftVector->operator[](0)->SetPosition(aPosition); 197 198 */ 198 //G4cout<<"Single hadron "<<LeftVector->operator[](0)->Get Position()<<" "<<LeftVector->operator[](0)->GetFormationTime()<<G4endl;199 //G4cout<<"Single hadron "<<LeftVector->operator[](0)->Get4Momentum()<<" "<<LeftVector->operator[](0)->Get4Momentum().mag()<<G4endl; 199 200 } else { // 2 hadrons created from qq-qqbar are stored 200 201 LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation()); … … 293 294 G4double TimeOftheStringCreation=theString.GetTimeOfCreation(); 294 295 G4ThreeVector PositionOftheStringCreation(theString.GetPosition()); 295 296 // Uzhi 11.07.09 297 //G4cout<<" String T Pos"<<TimeOftheStringCreation<<' '<<PositionOftheStringCreation<<G4endl; 296 298 /* // For large formation time open * 297 299 G4double TimeOftheStringCreation=theString.GetTimeOfCreation()+100*fermi; … … 313 315 } 314 316 */ 317 318 //G4cout<<"Lund Frag # hadrons"<<LeftVector->size()<<G4endl; // Uzhi 11.07.09 315 319 for(size_t C1 = 0; C1 < LeftVector->size(); C1++) 316 320 { … … 322 326 G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime()); 323 327 Momentum = toObserverFrame*Coordinate; 324 Hadron->SetFormationTime(TimeOftheStringCreation+Momentum.e()); 328 Hadron->SetFormationTime(TimeOftheStringCreation+Momentum.e() // Uzhi 11.07.09 329 -fermi/c_light); // Uzhi 11.07.09 325 330 G4ThreeVector aPosition(Momentum.vect()); 326 331 // Hadron->SetPosition(theString.GetPosition()+aPosition); 327 332 Hadron->SetPosition(PositionOftheStringCreation+aPosition); 328 333 //G4cout<<"Hadron "<<C1<<" "<<Hadron->GetPosition()/fermi<<" "<<Hadron->GetFormationTime()/fermi<<G4endl; 334 /* // Uzhi 11.07.09 335 G4cout<<C1<<' '<<Hadron->GetDefinition()->GetParticleName()<<G4endl; 336 G4cout<<Hadron->GetDefinition()->GetPDGMass()<<' ' 337 <<Hadron->Get4Momentum()<<G4endl; 338 G4cout<<Hadron->GetFormationTime()<<' ' 339 <<Hadron->GetPosition()<<' ' 340 <<Hadron->GetPosition().z()/fermi<<G4endl; 341 */ // Uzhi 329 342 }; 330 343 … … 519 532 // << " pt2max= " << pt2max ; 520 533 521 Pt=SampleQuarkPt(s qrt(pt2max)); Pt.setZ(0); G4double Pt2=Pt.mag2();534 Pt=SampleQuarkPt(std::sqrt(pt2max)); Pt.setZ(0); G4double Pt2=Pt.mag2(); 522 535 523 536 -
trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4QGSMFragmentation.cc
r1007 r1196 26 26 // 27 27 // $Id: G4QGSMFragmentation.cc,v 1.9 2008/06/23 08:35:55 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-0 2$28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 // ----------------------------------------------------------------------------- -
trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4VKinkyStringDecay.cc
r1007 r1196 26 26 // 27 27 // $Id: G4VKinkyStringDecay.cc,v 1.4 2008/04/25 14:20:14 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-0 2$28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // Maxim Komogorov 30 30 // -
trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4VLongitudinalStringDecay.cc
r1055 r1196 25 25 // 26 26 // 27 // $Id: G4VLongitudinalStringDecay.cc,v 1.1 4 2009/05/22 16:35:47 gunterExp $28 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // $Id: G4VLongitudinalStringDecay.cc,v 1.18 2009/08/03 13:21:25 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 // ----------------------------------------------------------------------------- … … 71 71 72 72 // Changable Parameters below. 73 74 73 SigmaQT = 0.5 * GeV; 75 74 … … 493 492 } else { 494 493 // sample in limited range 495 Pt = -std::log(CLHEP::RandFlat::shoot( exp(-sqr(ptMax)/sqr(SigmaQT)), 1.));494 Pt = -std::log(CLHEP::RandFlat::shoot(std::exp(-sqr(ptMax)/sqr(SigmaQT)), 1.)); 496 495 } 497 496 Pt = SigmaQT * std::sqrt(Pt); … … 505 504 { 506 505 507 //`yo-yo` formation time506 // `yo-yo` formation time 508 507 // const G4double kappa = 1.0 * GeV/fermi/4.; // Uzhi String tension 1.06.08 509 508 G4double kappa = GetStringTensionParameter(); 510 509 //G4cout<<"Kappa "<<kappa<<G4endl; // Uzhi 20.06.08 511 510 //G4int Uzhi; G4cin>>Uzhi; // Uzhi 20.06.08 511 //G4cout<<"Number of hadrons "<<Hadrons->size()<<G4endl; // Vova 512 512 for(size_t c1 = 0; c1 < Hadrons->size(); c1++) 513 513 { … … 521 521 G4double HadronE = Hadrons->operator[](c1)->Get4Momentum().e(); 522 522 G4double HadronPz = Hadrons->operator[](c1)->Get4Momentum().pz(); 523 Hadrons->operator[](c1)->SetFormationTime((theInitialStringMass - 2.*SumPz + HadronE - HadronPz)/(2.*kappa)); 524 G4ThreeVector aPosition(0, 0, (theInitialStringMass - 2.*SumE - HadronE + HadronPz)/(2.*kappa)); 523 Hadrons->operator[](c1)->SetFormationTime( 524 (theInitialStringMass - 2.*SumPz + HadronE - HadronPz)/(2.*kappa)/c_light); //Uzhi1.07.09c_light 525 526 G4ThreeVector aPosition(0, 0, 527 (theInitialStringMass - 2.*SumE - HadronE + HadronPz)/(2.*kappa)); 525 528 Hadrons->operator[](c1)->SetPosition(aPosition); 529 530 /* 531 G4cout<<"kappa "<<kappa<<G4endl; 532 G4cout<<c1<<" Partic time, position Old "<< 533 (theInitialStringMass - 2.*SumPz + HadronE - HadronPz)/(2.*kappa)<<' '<< 534 (theInitialStringMass - 2.*SumE - HadronE + HadronPz)/(2.*kappa)<<G4endl; 535 536 G4cout<<c1<<" Partic time, position New "<< 537 (theInitialStringMass - 2.*SumPz + HadronE - HadronPz)/(2.*kappa)/c_light<<' '<< 538 (theInitialStringMass - 2.*SumE - HadronE + HadronPz)/(2.*kappa)<<G4endl; 539 540 G4cout<<"fermi "<<fermi<<" 1/fermi "<<1./fermi<<' '<<"1*fermi/c "<<1.*fermi/c_light<<G4endl; 541 G4int Uzhi; G4cin>>Uzhi; // Uzhi 20.06.08 542 */ 526 543 } 527 544 } -
trunk/source/processes/hadronic/models/parton_string/management/History
r962 r1196 1 $Id: History,v 1. 4 2008/04/01 08:20:25vuzhinsk Exp $1 $Id: History,v 1.7 2009/07/17 12:41:20 vuzhinsk Exp $ 2 2 ------------------------------------------------------------------- 3 3 … … 15 15 * Please list in reverse chronological order (last date on top) 16 16 --------------------------------------------------------------- 17 17-July-2009 V. Uzhinsky (had-partonstring-mgt-V09-02-01) 18 A Status of nuclear nucleon involved in an interaction is introdused. 19 Status: 0 - spectator, 1 - involved nucleon, 2 - absorbed nucleon 20 (G4VSplitableHadron) 21 22 A connection between a participant nucleon and a nuclear nucleon was 23 introsuced in G4InteractionContent. 24 25 10-July-2009 V. Uzhinsky (had-partonstring-mgt-V09-02-00) 26 Introduction the right tag number. 27 28 9-July-2009 V. Uzhinsky (had-partonstring-mgt-V08-02-02) 29 - New field was added in G4VSplitableHadron class (G4 bool Activation) and 30 corresponding methods to operatite with it. It was needed for an 31 absorption of meson in nuclear collision generated by FTF. 17 32 18 33 24-Apr 2007 Gunter Folger (had-partonstring-mgt-V08-02-01) -
trunk/source/processes/hadronic/models/parton_string/management/include/G4EventGenerator.hh
r1007 r1196 26 26 // 27 27 // $Id: G4EventGenerator.hh,v 1.3 2006/06/29 20:55:13 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-0 2$28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 #ifndef G4EventGenerator_h -
trunk/source/processes/hadronic/models/parton_string/management/include/G4InteractionCode.hh
r1007 r1196 26 26 // 27 27 // $Id: G4InteractionCode.hh,v 1.3 2006/06/29 20:55:15 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-0 2$28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 #ifndef G4InteractionCode_h -
trunk/source/processes/hadronic/models/parton_string/management/include/G4InteractionContent.hh
r1007 r1196 25 25 // 26 26 // 27 // $Id: G4InteractionContent.hh,v 1. 4 2007/01/24 10:28:54 gunterExp $28 // GEANT4 tag $Name: geant4-09-0 2$27 // $Id: G4InteractionContent.hh,v 1.5 2009/07/17 12:36:41 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 … … 43 43 44 44 #include "G4VSplitableHadron.hh" 45 45 #include "G4Nucleon.hh" // Uzhi 16.07.09 46 46 class G4InteractionContent 47 47 { … … 58 58 G4VSplitableHadron * GetProjectile() const ; 59 59 G4VSplitableHadron * GetTarget() const; 60 61 void SetTargetNucleon(G4Nucleon * aNucleon); // Uzhi 16.07.09 62 G4Nucleon * GetTargetNucleon() const; // Uzhi 16.07.09 60 63 61 64 void SetTarget(G4VSplitableHadron *aTarget); … … 86 89 G4VSplitableHadron * theTarget; 87 90 G4VSplitableHadron * theProjectile; 91 G4Nucleon * theTargetNucleon; 88 92 89 93 G4int theNumberOfHard; … … 107 111 { 108 112 theTarget = aTarget; 113 } 114 115 inline void G4InteractionContent::SetTargetNucleon(G4Nucleon * aNucleon) // Uzhi 16.07.09 116 { 117 theTargetNucleon = aNucleon; 118 } 119 120 inline G4Nucleon * G4InteractionContent::GetTargetNucleon() const // Uzhi 16.07.09 121 { 122 return theTargetNucleon; 109 123 } 110 124 -
trunk/source/processes/hadronic/models/parton_string/management/include/G4PomeronCrossSection.hh
r1007 r1196 28 28 // 29 29 // $Id: G4PomeronCrossSection.hh,v 1.3 2006/06/29 20:55:19 gunter Exp $ 30 // GEANT4 tag $Name: geant4-09-0 2$30 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 31 31 // 32 32 #include "G4Proton.hh" -
trunk/source/processes/hadronic/models/parton_string/management/include/G4StringModel.hh
r1007 r1196 26 26 // 27 27 // $Id: G4StringModel.hh,v 1.3 2006/06/29 20:55:23 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-0 2$28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 #ifndef G4StringModel_h -
trunk/source/processes/hadronic/models/parton_string/management/include/G4VParticipants.hh
r1007 r1196 25 25 // 26 26 // 27 // $Id: G4VParticipants.hh,v 1. 4 2008/05/19 13:03:20 vuzhinskExp $28 // GEANT4 tag $Name: geant4-09-0 2$27 // $Id: G4VParticipants.hh,v 1.6 2009/11/19 14:23:09 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 … … 65 65 66 66 67 protected: 67 // protected: // Uzhi 26 July 09 68 68 69 69 … … 90 90 if ( theNucleus == NULL ) theNucleus = new G4Fancy3DNucleus(); 91 91 theNucleus->Init(theA, theZ); 92 theNucleus->SortNucleonsIn Z(); // Uzhi 16.05.08 Sorting of nucleon-Z92 theNucleus->SortNucleonsIncZ(); 93 93 } 94 94 -
trunk/source/processes/hadronic/models/parton_string/management/include/G4VPartonStringModel.hh
r1007 r1196 26 26 // 27 27 // $Id: G4VPartonStringModel.hh,v 1.3 2006/06/29 20:55:27 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-0 2$28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 #ifndef G4VPartonStringModel_h -
trunk/source/processes/hadronic/models/parton_string/management/include/G4VSplitableHadron.hh
r1007 r1196 25 25 // 26 26 // 27 // $Id: G4VSplitableHadron.hh,v 1. 4 2008/05/19 13:03:20vuzhinsk Exp $28 // GEANT4 tag $Name: geant4-09-0 2$27 // $Id: G4VSplitableHadron.hh,v 1.6 2009/07/17 12:36:41 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 … … 81 81 const G4ThreeVector & GetPosition() const; 82 82 83 void SetStatus(const G4int aStatus); // Uzhi 17.07.09 84 G4int GetStatus(); // Uzhi 17.07.09 85 83 86 virtual void SplitUp() = 0; 84 87 virtual G4Parton * GetNextParton() = 0 ; … … 106 109 G4int theCollisionCount; 107 110 111 G4int Status; // Uzhi 17.07.09 108 112 G4bool isSplit; 109 113 … … 165 169 } 166 170 171 inline void G4VSplitableHadron::SetStatus(G4int aStatus) // Uzhi 17.07.09 172 { 173 Status=aStatus; 174 } 175 176 inline G4int G4VSplitableHadron::GetStatus() // Uzhi 17.07.09 177 { 178 return Status; 179 } 167 180 168 181 -
trunk/source/processes/hadronic/models/parton_string/management/include/G4VStringFragmentation.hh
r1007 r1196 26 26 // 27 27 // $Id: G4VStringFragmentation.hh,v 1.3 2006/06/29 20:55:31 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-0 2$28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 #ifndef G4VStringFragmentation_h -
trunk/source/processes/hadronic/models/parton_string/management/include/G4VertexCode.hh
r1007 r1196 26 26 // 27 27 // $Id: G4VertexCode.hh,v 1.3 2006/06/29 20:55:33 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-0 2$28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 #ifndef G4VertexCode_h -
trunk/source/processes/hadronic/models/parton_string/management/src/G4EventGenerator.cc
r1007 r1196 26 26 // 27 27 // $Id: G4EventGenerator.cc,v 1.4 2006/06/29 20:55:37 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-0 2$28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 // G4EventGenerator -
trunk/source/processes/hadronic/models/parton_string/management/src/G4InteractionContent.cc
r1007 r1196 26 26 // 27 27 // $Id: G4InteractionContent.cc,v 1.4 2006/06/29 20:55:39 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-0 2$28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 // ------------------------------------------------------------ -
trunk/source/processes/hadronic/models/parton_string/management/src/G4PomeronCrossSection.cc
r1007 r1196 26 26 // 27 27 // $Id: G4PomeronCrossSection.cc,v 1.6 2006/11/07 12:51:39 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-0 2$28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 -
trunk/source/processes/hadronic/models/parton_string/management/src/G4StringModel.cc
r1007 r1196 26 26 // 27 27 // $Id: G4StringModel.cc,v 1.4 2006/06/29 20:55:45 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-0 2$28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 // G4StringModel -
trunk/source/processes/hadronic/models/parton_string/management/src/G4VParticipants.cc
r1007 r1196 26 26 // 27 27 // $Id: G4VParticipants.cc,v 1.3 2006/06/29 20:55:47 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-0 2$28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 // ------------------------------------------------------------ -
trunk/source/processes/hadronic/models/parton_string/management/src/G4VPartonStringModel.cc
r1007 r1196 25 25 // 26 26 // 27 // $Id: G4VPartonStringModel.cc,v 1. 5 2007/01/24 10:29:30 gunterExp $28 // GEANT4 tag $Name: geant4-09-0 2$27 // $Id: G4VPartonStringModel.cc,v 1.6 2009/10/05 12:52:48 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 //// ------------------------------------------------------------ … … 99 99 throw G4HadronicException(__FILE__, __LINE__, "G4VPartonStringModel::Scatter(): fails to generate strings"); 100 100 } 101 102 101 theThis->Init(theNucleus,thePrimary); 103 102 strings = GetStrings(); … … 106 105 G4KineticTrackVector * theResult = 0; 107 106 G4double stringEnergy(0); 107 108 108 for ( unsigned int astring=0; astring < strings->size(); astring++) 109 109 { … … 137 137 138 138 theResult = stringFragmentationModel->FragmentStrings(strings); 139 /* 140 G4cout<<"Size "<<theResult->size()<<G4endl; 141 for ( unsigned int i=0; i < theResult->size(); i++) 142 { 143 144 G4cout<<(*theResult)[i]->Get4Momentum()<<" "<<(*theResult)[i]->Get4Momentum().mag()<<G4endl;; 145 } 146 */ 139 147 std::for_each(strings->begin(), strings->end(), DeleteString() ); 140 148 delete strings; -
trunk/source/processes/hadronic/models/parton_string/management/src/G4VSplitableHadron.cc
r1007 r1196 25 25 // 26 26 // 27 // $Id: G4VSplitableHadron.cc,v 1. 5 2008/05/19 13:03:20vuzhinsk Exp $28 // GEANT4 tag $Name: geant4-09-0 2$27 // $Id: G4VSplitableHadron.cc,v 1.7 2009/07/17 12:36:41 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 … … 42 42 43 43 G4VSplitableHadron::G4VSplitableHadron() 44 : theDefinition(NULL), TimeOfCreation(0.), theCollisionCount(0), isSplit(false) // Uzhi 8.05.08 44 : theDefinition(NULL), TimeOfCreation(0.), theCollisionCount(0), // Uzhi 8.05.08 45 Status(0), isSplit(false) // Uzhi 17.07.09 45 46 { 46 47 } 47 48 48 49 G4VSplitableHadron::G4VSplitableHadron(const G4ReactionProduct & aPrimary) 49 : TimeOfCreation(0.), theCollisionCount(0), isSplit(false) // Uzhi 8.05.08 50 : TimeOfCreation(0.), theCollisionCount(0), // Uzhi 8.05.08 51 Status(0), isSplit(false) // Uzhi 17.07.09 50 52 { 51 53 theDefinition=aPrimary.GetDefinition(); … … 56 58 G4VSplitableHadron::G4VSplitableHadron(const G4Nucleon & aNucleon) 57 59 { 58 TimeOfCreation = 0.; // Uzhi 8.05.0860 TimeOfCreation = 0.; // Uzhi 8.05.08 59 61 theCollisionCount= 0; 60 62 isSplit = false; … … 62 64 the4Momentum =aNucleon.GetMomentum(); 63 65 thePosition =aNucleon.GetPosition(); 66 Status = 0; // Uzhi 17.07.09 64 67 } 65 68 … … 72 75 the4Momentum =aNucleon->Get4Momentum(); 73 76 thePosition =aNucleon->GetPosition(); 77 Status = 0; // Uzhi 17.07.09 74 78 } 75 79 … … 82 86 the4Momentum = right.Get4Momentum(); 83 87 thePosition = right.GetPosition(); 88 Status = 0; // Uzhi 17.07.09 84 89 } 85 90 -
trunk/source/processes/hadronic/models/parton_string/management/src/G4VStringFragmentation.cc
r1007 r1196 26 26 // 27 27 // $Id: G4VStringFragmentation.cc,v 1.4 2006/06/29 20:55:53 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-0 2$28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 // G4VStringFragmentation
Note: See TracChangeset
for help on using the changeset viewer.
