source: trunk/source/processes/hadronic/models/parton_string/management/src/G4FTFCrossSection.cc@ 1036

Last change on this file since 1036 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 10.8 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id: G4FTFCrossSection.cc,v 1.2 2007/04/24 10:37:10 gunter Exp $
28// GEANT4 tag $Name: geant4-09-01-patch-02 $
29//
30
31#include "G4FTFCrossSection.hh"
32
33G4FTFCrossSection::G4FTFCrossSection()
34{;}
35
36
37G4FTFCrossSection::~G4FTFCrossSection()
38{;}
39//**********************************************************************************************
40
41G4FTFCrossSection::G4FTFCrossSection(const G4ParticleDefinition * particle, G4double s)
42 {
43 G4int PDGcode = particle->GetPDGEncoding();
44 G4int absPDGcode = std::abs(PDGcode);
45 G4double Elab = (s - 2*0.88*GeV*GeV)/(2*0.939*GeV)/GeV;
46 G4double Plab = std::sqrt(Elab * Elab - 0.88);
47
48 G4double LogPlab = std::log( Plab );
49 G4double sqrLogPlab = LogPlab * LogPlab;
50
51//G4cout<<"G4FTFCrossSection Plab "<<Plab<<G4endl;
52
53 G4int NumberOfTargetProtons = 1; //aNucleus.GetZ(); // ??????????????????????
54 G4int NumberOfTargetNeutrons = 1; //aNucleus.GetN();
55 G4int NumberOfTargetNucleons = NumberOfTargetProtons + NumberOfTargetNeutrons;
56
57 G4double Xtotal, Xelastic;
58
59 if( absPDGcode > 1000 ) //------Projectile is baryon --------
60 {
61 G4double XtotPP = 48.0 + 0. *std::pow(Plab, 0. ) + 0.522*sqrLogPlab - 4.51*LogPlab;
62 G4double XtotPN = 47.3 + 0. *std::pow(Plab, 0. ) + 0.513*sqrLogPlab - 4.27*LogPlab;
63
64 G4double XelPP = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
65 G4double XelPN = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
66
67 Xtotal = ( NumberOfTargetProtons * XtotPP +
68 NumberOfTargetNeutrons * XtotPN ) / NumberOfTargetNucleons;
69 Xelastic = ( NumberOfTargetProtons * XelPP +
70 NumberOfTargetNeutrons * XelPN ) / NumberOfTargetNucleons;
71 }
72 else if( PDGcode == 211 ) //------Projectile is PionPlus -------
73 {
74 G4double XtotPiP = 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab;
75 G4double XtotPiN = 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab;
76
77 G4double XelPiP = 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab;
78 G4double XelPiN = 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab;
79
80 Xtotal = ( NumberOfTargetProtons * XtotPiP +
81 NumberOfTargetNeutrons * XtotPiN ) / NumberOfTargetNucleons;
82 Xelastic = ( NumberOfTargetProtons * XelPiP +
83 NumberOfTargetNeutrons * XelPiN ) / NumberOfTargetNucleons;
84 }
85 else if( PDGcode == -211 ) //------Projectile is PionMinus -------
86 {
87 G4double XtotPiP = 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab;
88 G4double XtotPiN = 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab;
89
90 G4double XelPiP = 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab;
91 G4double XelPiN = 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab;
92
93 Xtotal = ( NumberOfTargetProtons * XtotPiP +
94 NumberOfTargetNeutrons * XtotPiN ) / NumberOfTargetNucleons;
95 Xelastic = ( NumberOfTargetProtons * XelPiP +
96 NumberOfTargetNeutrons * XelPiN ) / NumberOfTargetNucleons;
97 }
98
99 else if( PDGcode == 111 ) //------Projectile is PionZero -------
100 {
101 G4double XtotPiP =(16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab + //Pi+
102 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab)/2; //Pi-
103
104 G4double XtotPiN =(33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab + //Pi+
105 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
106
107 G4double XelPiP =( 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab + //Pi+
108 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
109 G4double XelPiN =( 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab + //Pi+
110 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
111
112 Xtotal = ( NumberOfTargetProtons * XtotPiP +
113 NumberOfTargetNeutrons * XtotPiN ) / NumberOfTargetNucleons;
114 Xelastic = ( NumberOfTargetProtons * XelPiP +
115 NumberOfTargetNeutrons * XelPiN ) / NumberOfTargetNucleons;
116 }
117 else if( PDGcode == 321 ) //------Projectile is KaonPlus -------
118 {
119 G4double XtotKP = 18.1 + 0. *std::pow(Plab, 0. ) + 0.26 *sqrLogPlab - 1.0 *LogPlab;
120 G4double XtotKN = 18.7 + 0. *std::pow(Plab, 0. ) + 0.21 *sqrLogPlab - 0.89*LogPlab;
121
122 G4double XelKP = 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab;
123 G4double XelKN = 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab;
124
125 Xtotal = ( NumberOfTargetProtons * XtotKP +
126 NumberOfTargetNeutrons * XtotKN ) / NumberOfTargetNucleons;
127 Xelastic = ( NumberOfTargetProtons * XelKP +
128 NumberOfTargetNeutrons * XelKN ) / NumberOfTargetNucleons;
129 }
130 else if( PDGcode ==-321 ) //------Projectile is KaonMinus ------
131 {
132 G4double XtotKP = 32.1 + 0. *std::pow(Plab, 0. ) + 0.66 *sqrLogPlab - 5.6 *LogPlab;
133 G4double XtotKN = 25.2 + 0. *std::pow(Plab, 0. ) + 0.38 *sqrLogPlab - 2.9 *LogPlab;
134
135 G4double XelKP = 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab;
136 G4double XelKN = 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab;
137
138 Xtotal = ( NumberOfTargetProtons * XtotKP +
139 NumberOfTargetNeutrons * XtotKN ) / NumberOfTargetNucleons;
140 Xelastic = ( NumberOfTargetProtons * XelKP +
141 NumberOfTargetNeutrons * XelKN ) / NumberOfTargetNucleons;
142 }
143 else if( PDGcode == 311 ) //------Projectile is KaonZero ------
144 {
145 G4double XtotKP =( 18.1 + 0. *std::pow(Plab, 0. ) + 0.26 *sqrLogPlab - 1.0 *LogPlab + //K+
146 32.1 + 0. *std::pow(Plab, 0. ) + 0.66 *sqrLogPlab - 5.6 *LogPlab)/2; //K-
147 G4double XtotKN =( 18.7 + 0. *std::pow(Plab, 0. ) + 0.21 *sqrLogPlab - 0.89*LogPlab + //K+
148 25.2 + 0. *std::pow(Plab, 0. ) + 0.38 *sqrLogPlab - 2.9 *LogPlab)/2; //K-
149
150 G4double XelKP =( 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab + //K+
151 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab)/2; //K-
152 G4double XelKN =( 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab + //K+
153 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab)/2; //K-
154 Xtotal = ( NumberOfTargetProtons * XtotKP +
155 NumberOfTargetNeutrons * XtotKN ) / NumberOfTargetNucleons;
156 Xelastic = ( NumberOfTargetProtons * XelKP +
157 NumberOfTargetNeutrons * XelKN ) / NumberOfTargetNucleons;
158 }
159 else //------Projectile is undefined, Nucleon assumed
160 {
161 G4double XtotPP = 48.0 + 0. *std::pow(Plab, 0. ) + 0.522*sqrLogPlab - 4.51*LogPlab;
162 G4double XtotPN = 47.3 + 0. *std::pow(Plab, 0. ) + 0.513*sqrLogPlab - 4.27*LogPlab;
163
164 G4double XelPP = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
165 G4double XelPN = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
166
167 Xtotal = ( NumberOfTargetProtons * XtotPP +
168 NumberOfTargetNeutrons * XtotPN ) / NumberOfTargetNucleons;
169 Xelastic = ( NumberOfTargetProtons * XelPP +
170 NumberOfTargetNeutrons * XelPN ) / NumberOfTargetNucleons;
171 };
172
173 SetTotalCrossSection(Xtotal);
174 SetElastisCrossSection(Xelastic);
175 SetInelasticCrossSection(Xtotal-Xelastic);
176
177//G4cout<<"G4FTFCrossSection Xt Xel "<<Xtotal<<" "<<Xelastic<<G4endl;
178
179//-----------------------------------------------------------------------------------
180 SetSlope( Xtotal*Xtotal/16./pi/Xelastic/0.3894 ); // Slope parameter of elastic scattering
181 // (GeV/c)^(-2))
182 // Gaussian parametrization of
183 // elastic scattering amplitude assumed
184
185//G4cout<<"G4FTFCrossSection Slope "<<GetSlope()<<G4endl;
186
187//-----------------------------------------------------------------------------------
188 SetGamma0( GetSlope()*Xtotal/10./2./pi );
189//-----------------------------------------------------------------------------------
190
191//G4cout<<"G4FTFCrossSection Out"<<G4endl;
192
193 }
194//**********************************************************************************************
Note: See TracBrowser for help on using the repository browser.