| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // $Id: G4NeutrinoNucleusModel.cc 91806 2015-0 26 // $Id: G4NeutrinoNucleusModel.cc 91806 2015-08-06 12:20:45Z gcosmo $
27 // 27 //
28 // Geant4 Header : G4NeutrinoNucleusModel 28 // Geant4 Header : G4NeutrinoNucleusModel
29 // 29 //
30 // Author : V.Grichine 12.2.19 30 // Author : V.Grichine 12.2.19
31 // 31 //
32 32
33 #include "G4NeutrinoNucleusModel.hh" 33 #include "G4NeutrinoNucleusModel.hh"
34 34
35 #include "G4SystemOfUnits.hh" 35 #include "G4SystemOfUnits.hh"
36 #include "G4ParticleTable.hh" 36 #include "G4ParticleTable.hh"
37 #include "G4ParticleDefinition.hh" 37 #include "G4ParticleDefinition.hh"
38 #include "G4IonTable.hh" 38 #include "G4IonTable.hh"
39 #include "Randomize.hh" 39 #include "Randomize.hh"
40 #include "G4RandomDirection.hh" 40 #include "G4RandomDirection.hh"
41 41
42 //#include "G4Integrator.hh" 42 //#include "G4Integrator.hh"
43 #include "G4DataVector.hh" 43 #include "G4DataVector.hh"
44 #include "G4PhysicsTable.hh" 44 #include "G4PhysicsTable.hh"
45 45
46 #include "G4CascadeInterface.hh" 46 #include "G4CascadeInterface.hh"
47 // #include "G4BinaryCascade.hh" 47 // #include "G4BinaryCascade.hh"
48 #include "G4TheoFSGenerator.hh" 48 #include "G4TheoFSGenerator.hh"
49 #include "G4GeneratorPrecompoundInterface.hh" 49 #include "G4GeneratorPrecompoundInterface.hh"
50 #include "G4ExcitationHandler.hh" 50 #include "G4ExcitationHandler.hh"
51 #include "G4PreCompoundModel.hh" 51 #include "G4PreCompoundModel.hh"
52 #include "G4LundStringFragmentation.hh" 52 #include "G4LundStringFragmentation.hh"
53 #include "G4ExcitedStringDecay.hh" 53 #include "G4ExcitedStringDecay.hh"
54 #include "G4FTFModel.hh" 54 #include "G4FTFModel.hh"
55 // #include "G4BinaryCascade.hh" 55 // #include "G4BinaryCascade.hh"
56 #include "G4HadFinalState.hh" 56 #include "G4HadFinalState.hh"
57 #include "G4HadSecondary.hh" 57 #include "G4HadSecondary.hh"
58 #include "G4HadronicInteractionRegistry.hh" 58 #include "G4HadronicInteractionRegistry.hh"
59 // #include "G4INCLXXInterface.hh" 59 // #include "G4INCLXXInterface.hh"
60 #include "G4KineticTrack.hh" 60 #include "G4KineticTrack.hh"
61 #include "G4DecayKineticTracks.hh" 61 #include "G4DecayKineticTracks.hh"
62 #include "G4KineticTrackVector.hh" 62 #include "G4KineticTrackVector.hh"
63 #include "G4Fragment.hh" 63 #include "G4Fragment.hh"
64 #include "G4ReactionProductVector.hh" 64 #include "G4ReactionProductVector.hh"
65 65
66 // #include "G4QGSModel.hh" 66 // #include "G4QGSModel.hh"
67 // #include "G4QGSMFragmentation.hh" 67 // #include "G4QGSMFragmentation.hh"
68 // #include "G4QGSParticipants.hh" 68 // #include "G4QGSParticipants.hh"
69 69
70 #include "G4MuonMinus.hh" 70 #include "G4MuonMinus.hh"
71 #include "G4MuonPlus.hh" 71 #include "G4MuonPlus.hh"
72 #include "G4Nucleus.hh" 72 #include "G4Nucleus.hh"
73 #include "G4LorentzVector.hh" 73 #include "G4LorentzVector.hh"
74 #include "G4PhysicsModelCatalog.hh" 74 #include "G4PhysicsModelCatalog.hh"
75 75
76 using namespace std; 76 using namespace std;
77 using namespace CLHEP; 77 using namespace CLHEP;
78 78
79 const G4int G4NeutrinoNucleusModel::fResNumber 79 const G4int G4NeutrinoNucleusModel::fResNumber = 6;
80 80
81 const G4double G4NeutrinoNucleusModel::fResMas 81 const G4double G4NeutrinoNucleusModel::fResMass[6] = // [fResNumber] =
82 {2190., 1920., 1700., 1600., 1440., 1232. }; 82 {2190., 1920., 1700., 1600., 1440., 1232. };
83 83
84 const G4int G4NeutrinoNucleusModel::fClustNumb 84 const G4int G4NeutrinoNucleusModel::fClustNumber = 4;
85 85
86 const G4double G4NeutrinoNucleusModel::fMesMas 86 const G4double G4NeutrinoNucleusModel::fMesMass[4] = {1260., 980., 770., 139.57};
87 const G4int G4NeutrinoNucleusModel::fMesPDG 87 const G4int G4NeutrinoNucleusModel::fMesPDG[4] = {20213, 9000211, 213, 211};
88 88
89 // const G4double G4NeutrinoNucleusModel::fBar 89 // const G4double G4NeutrinoNucleusModel::fBarMass[4] = {1905., 1600., 1232., 939.57};
90 // const G4int G4NeutrinoNucleusModel::fBar 90 // const G4int G4NeutrinoNucleusModel::fBarPDG[4] = {2226, 32224, 2224, 2212};
91 91
92 const G4double G4NeutrinoNucleusModel::fBarMas 92 const G4double G4NeutrinoNucleusModel::fBarMass[4] = {1700., 1600., 1232., 939.57};
93 const G4int G4NeutrinoNucleusModel::fBarPDG 93 const G4int G4NeutrinoNucleusModel::fBarPDG[4] = {12224, 32224, 2224, 2212};
94 94
95 const G4double G4NeutrinoNucleusModel::fNuMuE 95 const G4double G4NeutrinoNucleusModel::fNuMuEnergyLogVector[50] = {
96 115.603, 133.424, 153.991, 177.729, 205.126, 2 96 115.603, 133.424, 153.991, 177.729, 205.126, 236.746, 273.24, 315.361, 363.973, 420.08, 484.836, 559.573, 645.832,
97 745.387, 860.289, 992.903, 1145.96, 1322.61, 1 97 745.387, 860.289, 992.903, 1145.96, 1322.61, 1526.49, 1761.8, 2033.38, 2346.83, 2708.59, 3126.12, 3608.02, 4164.19,
98 4806.1, 5546.97, 6402.04, 7388.91, 8527.92, 98 98 4806.1, 5546.97, 6402.04, 7388.91, 8527.92, 9842.5, 11359.7, 13110.8, 15131.9, 17464.5, 20156.6, 23263.8, 26849.9,
99 30988.8, 35765.7, 41279, 47642.2, 54986.3, 634 99 30988.8, 35765.7, 41279, 47642.2, 54986.3, 63462.4, 73245.2, 84536, 97567.2, 112607, 129966 };
100 100
101 101
102 G4double G4NeutrinoNucleusModel::fNuMuXarrayKR 102 G4double G4NeutrinoNucleusModel::fNuMuXarrayKR[50][51] = {{1.0}};
103 G4double G4NeutrinoNucleusModel::fNuMuXdistrKR 103 G4double G4NeutrinoNucleusModel::fNuMuXdistrKR[50][50] = {{1.0}};
104 G4double G4NeutrinoNucleusModel::fNuMuQarrayKR 104 G4double G4NeutrinoNucleusModel::fNuMuQarrayKR[50][51][51] = {{{1.0}}};
105 G4double G4NeutrinoNucleusModel::fNuMuQdistrKR 105 G4double G4NeutrinoNucleusModel::fNuMuQdistrKR[50][51][50] = {{{1.0}}};
106 106
107 /////////////////////////////////////////// 107 ///////////////////////////////////////////
108 108
109 G4NeutrinoNucleusModel::G4NeutrinoNucleusModel 109 G4NeutrinoNucleusModel::G4NeutrinoNucleusModel(const G4String& name)
110 : G4HadronicInteraction(name), fSecID(-1) 110 : G4HadronicInteraction(name), fSecID(-1)
111 { 111 {
112 SetMinEnergy( 0.0*GeV ); 112 SetMinEnergy( 0.0*GeV );
113 SetMaxEnergy( 100.*TeV ); 113 SetMaxEnergy( 100.*TeV );
114 SetMinEnergy(1.e-6*eV); 114 SetMinEnergy(1.e-6*eV);
115 115
116 fNbin = 50; 116 fNbin = 50;
117 fEindex = fXindex = fQindex = 0; 117 fEindex = fXindex = fQindex = 0;
118 fOnePionIndex = 58; 118 fOnePionIndex = 58;
119 fIndex = 50; 119 fIndex = 50;
120 fCascade = fString = fProton = f2p2h = fBrea 120 fCascade = fString = fProton = f2p2h = fBreak = false;
121 121
122 fNuEnergy = fQ2 = fQtransfer = fXsample = 122 fNuEnergy = fQ2 = fQtransfer = fXsample = fDp = fTr = 0.;
123 fCosTheta = fCosThetaPi = 1.; 123 fCosTheta = fCosThetaPi = 1.;
124 fEmuPi = fW2 = fW2pi = 0.; 124 fEmuPi = fW2 = fW2pi = 0.;
125 125
126 fMu = 105.6583745*MeV; 126 fMu = 105.6583745*MeV;
127 fMpi = 139.57018*MeV; 127 fMpi = 139.57018*MeV;
128 fM1 = 939.5654133*MeV; // for nu_mu -> mu-, 128 fM1 = 939.5654133*MeV; // for nu_mu -> mu-, and n -> p
129 fM2 = 938.2720813*MeV; 129 fM2 = 938.2720813*MeV;
130 130
131 fEmu = fMu; 131 fEmu = fMu;
132 fEx = fM1; 132 fEx = fM1;
133 fMr = 1232.*MeV; 133 fMr = 1232.*MeV;
134 fMt = fM2; // threshold for N*-diffraction 134 fMt = fM2; // threshold for N*-diffraction
135 135
136 fMinNuEnergy = GetMinNuMuEnergy(); 136 fMinNuEnergy = GetMinNuMuEnergy();
137 137
138 fLVh = G4LorentzVector(0.,0.,0.,0.); 138 fLVh = G4LorentzVector(0.,0.,0.,0.);
139 fLVl = G4LorentzVector(0.,0.,0.,0.); 139 fLVl = G4LorentzVector(0.,0.,0.,0.);
140 fLVt = G4LorentzVector(0.,0.,0.,0.); 140 fLVt = G4LorentzVector(0.,0.,0.,0.);
141 fLVcpi = G4LorentzVector(0.,0.,0.,0.); 141 fLVcpi = G4LorentzVector(0.,0.,0.,0.);
142 142
143 fQEratioA = 0.5; // mean value around 1 GeV 143 fQEratioA = 0.5; // mean value around 1 GeV neutrino beams
144 144
145 theMuonMinus = G4MuonMinus::MuonMinus(); 145 theMuonMinus = G4MuonMinus::MuonMinus();
146 theMuonPlus = G4MuonPlus::MuonPlus(); 146 theMuonPlus = G4MuonPlus::MuonPlus();
147 147
148 // PDG2016: sin^2 theta Weinberg 148 // PDG2016: sin^2 theta Weinberg
149 149
150 fSin2tW = 0.23129; // 0.2312; 150 fSin2tW = 0.23129; // 0.2312;
151 151
152 fCutEnergy = 0.; // default value 152 fCutEnergy = 0.; // default value
153 153
154 154
155 /* 155 /*
156 // G4VPreCompoundModel* ptr; 156 // G4VPreCompoundModel* ptr;
157 // reuse existing pre-compound model as in b 157 // reuse existing pre-compound model as in binary cascade
158 158
159 fPrecoInterface = new G4GeneratorPrecompou 159 fPrecoInterface = new G4GeneratorPrecompoundInterface ;
160 160
161 if( !fPreCompound ) 161 if( !fPreCompound )
162 { 162 {
163 G4HadronicInteraction* p = 163 G4HadronicInteraction* p =
164 G4HadronicInteractionRegistry::Instanc 164 G4HadronicInteractionRegistry::Instance()->FindModel("PRECO");
165 G4VPreCompoundModel* fPreCompound = stat 165 G4VPreCompoundModel* fPreCompound = static_cast<G4VPreCompoundModel*>(p);
166 166
167 if(!fPreCompound) 167 if(!fPreCompound)
168 { 168 {
169 fPreCompound = new G4PreCompoundModel(); 169 fPreCompound = new G4PreCompoundModel();
170 } 170 }
171 fPrecoInterface->SetDeExcitation(fPreCom 171 fPrecoInterface->SetDeExcitation(fPreCompound);
172 } 172 }
173 fDeExcitation = GetDeExcitation()->GetExci 173 fDeExcitation = GetDeExcitation()->GetExcitationHandler();
174 */ 174 */
175 175
176 fDeExcitation = new G4ExcitationHandler(); 176 fDeExcitation = new G4ExcitationHandler();
177 fPreCompound = new G4PreCompoundModel(fDe 177 fPreCompound = new G4PreCompoundModel(fDeExcitation);
178 fPrecoInterface = new G4GeneratorPrecompoun 178 fPrecoInterface = new G4GeneratorPrecompoundInterface ;
179 fPrecoInterface->SetDeExcitation(fPreCompoun 179 fPrecoInterface->SetDeExcitation(fPreCompound);
180 180
181 fPDGencoding = 0; // unphysical as default 181 fPDGencoding = 0; // unphysical as default
182 fRecoil = nullptr; 182 fRecoil = nullptr;
183 183
184 // Creator model ID 184 // Creator model ID
185 fSecID = G4PhysicsModelCatalog::GetModelID( 185 fSecID = G4PhysicsModelCatalog::GetModelID( "model_" + GetModelName() );
186 } 186 }
187 187
188 188
189 G4NeutrinoNucleusModel::~G4NeutrinoNucleusMode 189 G4NeutrinoNucleusModel::~G4NeutrinoNucleusModel()
190 { 190 {
191 if(fPrecoInterface) delete fPrecoInterface; 191 if(fPrecoInterface) delete fPrecoInterface;
192 } 192 }
193 193
194 194
195 void G4NeutrinoNucleusModel::ModelDescription( 195 void G4NeutrinoNucleusModel::ModelDescription(std::ostream& outFile) const
196 { 196 {
197 197
198 outFile << "G4NeutrinoNucleusModel is a ne 198 outFile << "G4NeutrinoNucleusModel is a neutrino-nucleus general\n"
199 << "model which uses the standard 199 << "model which uses the standard model \n"
200 << "transfer parameterization. Th 200 << "transfer parameterization. The model is fully relativistic\n";
201 201
202 } 202 }
203 203
204 ////////////////////////////////////////////// 204 /////////////////////////////////////////////////////////
205 205
206 G4bool G4NeutrinoNucleusModel::IsApplicable(co 206 G4bool G4NeutrinoNucleusModel::IsApplicable(const G4HadProjectile & aPart,
207 G4Nucleus & ) 207 G4Nucleus & )
208 { 208 {
209 G4bool result = false; 209 G4bool result = false;
210 G4String pName = aPart.GetDefinition()->GetP 210 G4String pName = aPart.GetDefinition()->GetParticleName();
211 G4double energy = aPart.GetTotalEnergy(); 211 G4double energy = aPart.GetTotalEnergy();
212 212
213 if( pName == "nu_mu" // || pName == "anti_n 213 if( pName == "nu_mu" // || pName == "anti_nu_mu" )
214 && 214 &&
215 energy > fMinNuEnergy 215 energy > fMinNuEnergy )
216 { 216 {
217 result = true; 217 result = true;
218 } 218 }
219 219
220 return result; 220 return result;
221 } 221 }
222 222
223 ////////////////////////////////////// 223 //////////////////////////////////////
224 224
225 G4double G4NeutrinoNucleusModel::SampleXkr(G4d 225 G4double G4NeutrinoNucleusModel::SampleXkr(G4double energy)
226 { 226 {
227 G4int i(0), nBin(50); 227 G4int i(0), nBin(50);
228 G4double xx(0.), prob = G4UniformRand(); 228 G4double xx(0.), prob = G4UniformRand();
229 229
230 for( i = 0; i < nBin; ++i ) 230 for( i = 0; i < nBin; ++i )
231 { 231 {
232 if( energy <= fNuMuEnergyLogVector[i] ) br 232 if( energy <= fNuMuEnergyLogVector[i] ) break;
233 } 233 }
234 if( i <= 0) // E-edge 234 if( i <= 0) // E-edge
235 { 235 {
236 fEindex = 0; 236 fEindex = 0;
237 xx = GetXkr( 0, prob); 237 xx = GetXkr( 0, prob);
238 } 238 }
239 else if ( i >= nBin) 239 else if ( i >= nBin)
240 { 240 {
241 fEindex = nBin-1; 241 fEindex = nBin-1;
242 xx = GetXkr( nBin-1, prob); 242 xx = GetXkr( nBin-1, prob);
243 } 243 }
244 else 244 else
245 { 245 {
246 fEindex = i; 246 fEindex = i;
247 G4double x1 = GetXkr(i-1,prob); 247 G4double x1 = GetXkr(i-1,prob);
248 G4double x2 = GetXkr(i,prob); 248 G4double x2 = GetXkr(i,prob);
249 249
250 G4double e1 = G4Log(fNuMuEnergyLogVector[i 250 G4double e1 = G4Log(fNuMuEnergyLogVector[i-1]);
251 G4double e2 = G4Log(fNuMuEnergyLogVector[i 251 G4double e2 = G4Log(fNuMuEnergyLogVector[i]);
252 G4double e = G4Log(energy); 252 G4double e = G4Log(energy);
253 253
254 if( e2 <= e1) xx = x1 + G4UniformRand()*(x 254 if( e2 <= e1) xx = x1 + G4UniformRand()*(x2-x1);
255 else xx = x1 + (e-e1)*(x2-x1)/(e2 255 else xx = x1 + (e-e1)*(x2-x1)/(e2-e1); // lin in energy log-scale
256 } 256 }
257 return xx; 257 return xx;
258 } 258 }
259 259
260 ////////////////////////////////////////////// 260 //////////////////////////////////////////////
261 // 261 //
262 // sample X according to prob (xmin,1) at a gi 262 // sample X according to prob (xmin,1) at a given energy index iEnergy
263 263
264 G4double G4NeutrinoNucleusModel::GetXkr(G4int 264 G4double G4NeutrinoNucleusModel::GetXkr(G4int iEnergy, G4double prob)
265 { 265 {
266 G4int i(0), nBin=50; 266 G4int i(0), nBin=50;
267 G4double xx(0.); 267 G4double xx(0.);
268 268
269 for( i = 0; i < nBin; ++i ) 269 for( i = 0; i < nBin; ++i )
270 { 270 {
271 if( prob <= fNuMuXdistrKR[iEnergy][i] ) 271 if( prob <= fNuMuXdistrKR[iEnergy][i] )
272 break; 272 break;
273 } 273 }
274 if(i <= 0 ) // X-edge 274 if(i <= 0 ) // X-edge
275 { 275 {
276 fXindex = 0; 276 fXindex = 0;
277 xx = fNuMuXarrayKR[iEnergy][0]; 277 xx = fNuMuXarrayKR[iEnergy][0];
278 } 278 }
279 if ( i >= nBin ) 279 if ( i >= nBin )
280 { 280 {
281 fXindex = nBin; 281 fXindex = nBin;
282 xx = fNuMuXarrayKR[iEnergy][nBin]; 282 xx = fNuMuXarrayKR[iEnergy][nBin];
283 } 283 }
284 else 284 else
285 { 285 {
286 fXindex = i; 286 fXindex = i;
287 G4double x1 = fNuMuXarrayKR[iEnergy][i]; 287 G4double x1 = fNuMuXarrayKR[iEnergy][i];
288 G4double x2 = fNuMuXarrayKR[iEnergy][i+1]; 288 G4double x2 = fNuMuXarrayKR[iEnergy][i+1];
289 289
290 G4double p1 = 0.; 290 G4double p1 = 0.;
291 291
292 if( i > 0 ) p1 = fNuMuXdistrKR[iEnergy][i- 292 if( i > 0 ) p1 = fNuMuXdistrKR[iEnergy][i-1];
293 293
294 G4double p2 = fNuMuXdistrKR[iEnergy][i]; 294 G4double p2 = fNuMuXdistrKR[iEnergy][i];
295 295
296 if( p2 <= p1 ) xx = x1 + G4UniformRand()*( 296 if( p2 <= p1 ) xx = x1 + G4UniformRand()*(x2-x1);
297 else xx = x1 + (prob-p1)*(x2-x1) 297 else xx = x1 + (prob-p1)*(x2-x1)/(p2-p1);
298 } 298 }
299 return xx; 299 return xx;
300 } 300 }
301 301
302 ////////////////////////////////////// 302 //////////////////////////////////////
303 // 303 //
304 // Sample fQtransfer at a given Enu and fX 304 // Sample fQtransfer at a given Enu and fX
305 305
306 G4double G4NeutrinoNucleusModel::SampleQkr( G4 306 G4double G4NeutrinoNucleusModel::SampleQkr( G4double energy, G4double xx)
307 { 307 {
308 G4int nBin(50), iE=fEindex, jX=fXindex; 308 G4int nBin(50), iE=fEindex, jX=fXindex;
309 G4double qq(0.), qq1(0.), qq2(0.); 309 G4double qq(0.), qq1(0.), qq2(0.);
310 G4double prob = G4UniformRand(); 310 G4double prob = G4UniformRand();
311 311
312 // first E 312 // first E
313 313
314 if( iE <= 0 ) 314 if( iE <= 0 )
315 { 315 {
316 qq1 = GetQkr( 0, jX, prob); 316 qq1 = GetQkr( 0, jX, prob);
317 } 317 }
318 else if ( iE >= nBin-1) 318 else if ( iE >= nBin-1)
319 { 319 {
320 qq1 = GetQkr( nBin-1, jX, prob); 320 qq1 = GetQkr( nBin-1, jX, prob);
321 } 321 }
322 else 322 else
323 { 323 {
324 G4double q1 = GetQkr(iE-1,jX, prob); 324 G4double q1 = GetQkr(iE-1,jX, prob);
325 G4double q2 = GetQkr(iE,jX, prob); 325 G4double q2 = GetQkr(iE,jX, prob);
326 326
327 G4double e1 = G4Log(fNuMuEnergyLogVector[i 327 G4double e1 = G4Log(fNuMuEnergyLogVector[iE-1]);
328 G4double e2 = G4Log(fNuMuEnergyLogVector[i 328 G4double e2 = G4Log(fNuMuEnergyLogVector[iE]);
329 G4double e = G4Log(energy); 329 G4double e = G4Log(energy);
330 330
331 if( e2 <= e1) qq1 = q1 + G4UniformRand()*( 331 if( e2 <= e1) qq1 = q1 + G4UniformRand()*(q2-q1);
332 else qq1 = q1 + (e-e1)*(q2-q1)/(e 332 else qq1 = q1 + (e-e1)*(q2-q1)/(e2-e1); // lin in energy log-scale
333 } 333 }
334 334
335 // then X 335 // then X
336 336
337 if( jX <= 0 ) 337 if( jX <= 0 )
338 { 338 {
339 qq2 = GetQkr( iE, 0, prob); 339 qq2 = GetQkr( iE, 0, prob);
340 } 340 }
341 else if ( jX >= nBin) 341 else if ( jX >= nBin)
342 { 342 {
343 qq2 = GetQkr( iE, nBin, prob); 343 qq2 = GetQkr( iE, nBin, prob);
344 } 344 }
345 else 345 else
346 { 346 {
347 G4double q1 = GetQkr(iE,jX-1, prob); 347 G4double q1 = GetQkr(iE,jX-1, prob);
348 G4double q2 = GetQkr(iE,jX, prob); 348 G4double q2 = GetQkr(iE,jX, prob);
349 349
350 G4double e1 = G4Log(fNuMuXarrayKR[iE][jX-1 350 G4double e1 = G4Log(fNuMuXarrayKR[iE][jX-1]);
351 G4double e2 = G4Log(fNuMuXarrayKR[iE][jX]) 351 G4double e2 = G4Log(fNuMuXarrayKR[iE][jX]);
352 G4double e = G4Log(xx); 352 G4double e = G4Log(xx);
353 353
354 if( e2 <= e1) qq2 = q1 + G4UniformRand()*( 354 if( e2 <= e1) qq2 = q1 + G4UniformRand()*(q2-q1);
355 else qq2 = q1 + (e-e1)*(q2-q1)/(e 355 else qq2 = q1 + (e-e1)*(q2-q1)/(e2-e1); // lin in energy log-scale
356 } 356 }
357 qq = 0.5*(qq1+qq2); 357 qq = 0.5*(qq1+qq2);
358 358
359 return qq; 359 return qq;
360 } 360 }
361 361
362 ////////////////////////////////////////////// 362 //////////////////////////////////////////////
363 // 363 //
364 // sample Q according to prob (qmin,qmax) at a 364 // sample Q according to prob (qmin,qmax) at a given energy index iE and X index jX
365 365
366 G4double G4NeutrinoNucleusModel::GetQkr( G4int 366 G4double G4NeutrinoNucleusModel::GetQkr( G4int iE, G4int jX, G4double prob )
367 { 367 {
368 G4int i(0), nBin=50; 368 G4int i(0), nBin=50;
369 G4double qq(0.); 369 G4double qq(0.);
370 370
371 for( i = 0; i < nBin; ++i ) 371 for( i = 0; i < nBin; ++i )
372 { 372 {
373 if( prob <= fNuMuQdistrKR[iE][jX][i] ) 373 if( prob <= fNuMuQdistrKR[iE][jX][i] )
374 break; 374 break;
375 } 375 }
376 if(i <= 0 ) // Q-edge 376 if(i <= 0 ) // Q-edge
377 { 377 {
378 fQindex = 0; 378 fQindex = 0;
379 qq = fNuMuQarrayKR[iE][jX][0]; 379 qq = fNuMuQarrayKR[iE][jX][0];
380 } 380 }
381 if ( i >= nBin ) 381 if ( i >= nBin )
382 { 382 {
383 fQindex = nBin; 383 fQindex = nBin;
384 qq = fNuMuQarrayKR[iE][jX][nBin]; 384 qq = fNuMuQarrayKR[iE][jX][nBin];
385 } 385 }
386 else 386 else
387 { 387 {
388 fQindex = i; 388 fQindex = i;
389 G4double q1 = fNuMuQarrayKR[iE][jX][i]; 389 G4double q1 = fNuMuQarrayKR[iE][jX][i];
390 G4double q2 = fNuMuQarrayKR[iE][jX][i+1]; 390 G4double q2 = fNuMuQarrayKR[iE][jX][i+1];
391 391
392 G4double p1 = 0.; 392 G4double p1 = 0.;
393 393
394 if( i > 0 ) p1 = fNuMuQdistrKR[iE][jX][i-1 394 if( i > 0 ) p1 = fNuMuQdistrKR[iE][jX][i-1];
395 395
396 G4double p2 = fNuMuQdistrKR[iE][jX][i]; 396 G4double p2 = fNuMuQdistrKR[iE][jX][i];
397 397
398 if( p2 <= p1 ) qq = q1 + G4UniformRand()*( 398 if( p2 <= p1 ) qq = q1 + G4UniformRand()*(q2-q1);
399 else qq = q1 + (prob-p1)*(q2-q1) 399 else qq = q1 + (prob-p1)*(q2-q1)/(p2-p1);
400 } 400 }
401 return qq; 401 return qq;
402 } 402 }
403 403
404 404
405 ////////////////////////////////////////////// 405 ///////////////////////////////////////////////////////////
406 // 406 //
407 // Final meson to theParticleChange 407 // Final meson to theParticleChange
408 408
409 void G4NeutrinoNucleusModel::FinalMeson( G4Lor 409 void G4NeutrinoNucleusModel::FinalMeson( G4LorentzVector & lvM, G4int, G4int pdgM) // qM
410 { 410 {
411 G4int pdg = pdgM; 411 G4int pdg = pdgM;
412 // if ( qM == 0 ) pdg = pdgM - 100; 412 // if ( qM == 0 ) pdg = pdgM - 100;
413 // else if ( qM == -1 ) pdg = -pdgM; 413 // else if ( qM == -1 ) pdg = -pdgM;
414 414
415 if( pdg == 211 || pdg == -211 || pdg == 111) 415 if( pdg == 211 || pdg == -211 || pdg == 111) // pions
416 { 416 {
417 G4ParticleDefinition* pd2 = G4ParticleTabl 417 G4ParticleDefinition* pd2 = G4ParticleTable::GetParticleTable()->FindParticle(pdg);
418 G4DynamicParticle* dp2 = new G4DynamicP 418 G4DynamicParticle* dp2 = new G4DynamicParticle( pd2, lvM);
419 theParticleChange.AddSecondary( dp2, fSecI 419 theParticleChange.AddSecondary( dp2, fSecID );
420 } 420 }
421 else // meson resonances 421 else // meson resonances
422 { 422 {
423 G4ParticleDefinition* rePart = G4ParticleT 423 G4ParticleDefinition* rePart = G4ParticleTable::GetParticleTable()->
424 FindParticle(pdg); 424 FindParticle(pdg);
425 G4KineticTrack ddkt( rePart, 0., G4ThreeVe 425 G4KineticTrack ddkt( rePart, 0., G4ThreeVector(0.,0.,0.), lvM);
426 G4KineticTrackVector* ddktv = ddkt.Decay() 426 G4KineticTrackVector* ddktv = ddkt.Decay();
427 427
428 G4DecayKineticTracks decay( ddktv ); 428 G4DecayKineticTracks decay( ddktv );
429 429
430 for( unsigned int i = 0; i < ddktv->size() 430 for( unsigned int i = 0; i < ddktv->size(); i++ ) // add products to partchange
431 { 431 {
432 G4DynamicParticle * aNew = 432 G4DynamicParticle * aNew =
433 new G4DynamicParticle( ddktv->operator[] 433 new G4DynamicParticle( ddktv->operator[](i)->GetDefinition(),
434 ddktv->operator[] 434 ddktv->operator[](i)->Get4Momentum());
435 435
436 // G4cout<<" "<<i<<", "<<aNew->Get 436 // G4cout<<" "<<i<<", "<<aNew->GetDefinition()->GetParticleName()<<", "<<aNew->Get4Momentum()<<G4endl;
437 437
438 theParticleChange.AddSecondary( aNew, fS 438 theParticleChange.AddSecondary( aNew, fSecID );
439 delete ddktv->operator[](i); 439 delete ddktv->operator[](i);
440 } 440 }
441 delete ddktv; 441 delete ddktv;
442 } 442 }
443 } 443 }
444 444
445 ////////////////////////////////////////////// 445 ////////////////////////////////////////////////////////
446 // 446 //
447 // Final barion to theParticleChange, and reco 447 // Final barion to theParticleChange, and recoil nucleus treatment
448 448
449 void G4NeutrinoNucleusModel::FinalBarion( G4Lo 449 void G4NeutrinoNucleusModel::FinalBarion( G4LorentzVector & lvB, G4int, G4int pdgB) // qB
450 { 450 {
451 G4int A(0), Z(0), pdg = pdgB; 451 G4int A(0), Z(0), pdg = pdgB;
452 // G4bool FiNucleon(false); 452 // G4bool FiNucleon(false);
453 453
454 // if ( qB == 1 ) pdg = pdgB - 10; 454 // if ( qB == 1 ) pdg = pdgB - 10;
455 // else if ( qB == 0 ) pdg = pdgB - 110; 455 // else if ( qB == 0 ) pdg = pdgB - 110;
456 // else if ( qB == -1 ) pdg = pdgB - 1110; 456 // else if ( qB == -1 ) pdg = pdgB - 1110;
457 457
458 if( pdg == 2212 || pdg == 2112) 458 if( pdg == 2212 || pdg == 2112)
459 { 459 {
460 fMr = G4ParticleTable::GetParticleTable()- 460 fMr = G4ParticleTable::GetParticleTable()->FindParticle(pdg)->GetPDGMass();
461 // FiNucleon = true; 461 // FiNucleon = true;
462 } 462 }
463 else fMr = lvB.m(); 463 else fMr = lvB.m();
464 464
465 G4ThreeVector bst = fLVt.boostVector(); 465 G4ThreeVector bst = fLVt.boostVector();
466 lvB.boost(-bst); // in fLVt rest system 466 lvB.boost(-bst); // in fLVt rest system
467 467
468 G4double eX = lvB.e(); 468 G4double eX = lvB.e();
469 G4double det(0.), det2(0.), rM(0.), mX = lvB 469 G4double det(0.), det2(0.), rM(0.), mX = lvB.m();
470 G4ThreeVector dX = (lvB.vect()).unit(); 470 G4ThreeVector dX = (lvB.vect()).unit();
471 G4double pX = sqrt(eX*eX-mX*mX); 471 G4double pX = sqrt(eX*eX-mX*mX);
472 472
473 if( fRecoil ) 473 if( fRecoil )
474 { 474 {
475 Z = fRecoil->GetZ_asInt(); 475 Z = fRecoil->GetZ_asInt();
476 A = fRecoil->GetA_asInt(); 476 A = fRecoil->GetA_asInt();
477 rM = fRecoil->AtomicMass(A,Z); //->AtomicM 477 rM = fRecoil->AtomicMass(A,Z); //->AtomicMass(); //
478 rM = fLVt.m(); 478 rM = fLVt.m();
479 } 479 }
480 else // A=0 nu+p 480 else // A=0 nu+p
481 { 481 {
482 A = 0; 482 A = 0;
483 Z = 1; 483 Z = 1;
484 rM = electron_mass_c2; 484 rM = electron_mass_c2;
485 } 485 }
486 // G4cout<<A<<", "; 486 // G4cout<<A<<", ";
487 487
488 G4double sumE = eX + rM; 488 G4double sumE = eX + rM;
489 G4double B = sumE*sumE + rM*rM - fMr*fMr - p 489 G4double B = sumE*sumE + rM*rM - fMr*fMr - pX*pX;
490 G4double a = 4.*(sumE*sumE - pX*pX); 490 G4double a = 4.*(sumE*sumE - pX*pX);
491 G4double b = -4.*B*pX; 491 G4double b = -4.*B*pX;
492 G4double c = 4.*sumE*sumE*rM*rM - B*B; 492 G4double c = 4.*sumE*sumE*rM*rM - B*B;
493 det2 = b*b-4.*a*c; 493 det2 = b*b-4.*a*c;
494 if( det2 <= 0. ) det = 0.; 494 if( det2 <= 0. ) det = 0.;
495 else det = sqrt(det2); 495 else det = sqrt(det2);
496 G4double dP = 0.5*(-b - det )/a; 496 G4double dP = 0.5*(-b - det )/a;
497 497
498 fDp = dP; 498 fDp = dP;
499 499
500 pX -= dP; 500 pX -= dP;
501 501
502 if(pX < 0.) pX = 0.; 502 if(pX < 0.) pX = 0.;
503 503
504 // if( A == 0 ) G4cout<<pX/MeV<<", "; 504 // if( A == 0 ) G4cout<<pX/MeV<<", ";
505 505
506 eX = sqrt( pX*pX + fMr*fMr ); 506 eX = sqrt( pX*pX + fMr*fMr );
507 G4LorentzVector lvN( pX*dX, eX ); 507 G4LorentzVector lvN( pX*dX, eX );
508 lvN.boost(bst); // back to lab 508 lvN.boost(bst); // back to lab
509 509
510 if( pdg == 2212 || pdg == 2112) // nucleons 510 if( pdg == 2212 || pdg == 2112) // nucleons mX >= fMr, dP >= 0
511 { 511 {
512 G4ParticleDefinition* pd2 = G4ParticleTabl 512 G4ParticleDefinition* pd2 = G4ParticleTable::GetParticleTable()->FindParticle(pdg);
513 G4DynamicParticle* dp2 = new G4DynamicP 513 G4DynamicParticle* dp2 = new G4DynamicParticle( pd2, lvN);
514 theParticleChange.AddSecondary( dp2, fSecI 514 theParticleChange.AddSecondary( dp2, fSecID );
515 515
516 } 516 }
517 else // delta resonances 517 else // delta resonances
518 { 518 {
519 G4ParticleDefinition* rePart = G4ParticleT 519 G4ParticleDefinition* rePart = G4ParticleTable::GetParticleTable()->FindParticle(pdg);
520 G4KineticTrack ddkt( rePart, 0., G4ThreeVe 520 G4KineticTrack ddkt( rePart, 0., G4ThreeVector( 0., 0., 0.), lvN);
521 G4KineticTrackVector* ddktv = ddkt.Decay() 521 G4KineticTrackVector* ddktv = ddkt.Decay();
522 522
523 G4DecayKineticTracks decay( ddktv ); 523 G4DecayKineticTracks decay( ddktv );
524 524
525 for( unsigned int i = 0; i < ddktv->size() 525 for( unsigned int i = 0; i < ddktv->size(); i++ ) // add products to partchange
526 { 526 {
527 G4DynamicParticle * aNew = 527 G4DynamicParticle * aNew =
528 new G4DynamicParticle( ddktv->operator[] 528 new G4DynamicParticle( ddktv->operator[](i)->GetDefinition(),
529 ddktv->operator[] 529 ddktv->operator[](i)->Get4Momentum() );
530 530
531 // G4cout<<" "<<i<<", "<<aNew->Get 531 // G4cout<<" "<<i<<", "<<aNew->GetDefinition()->GetParticleName()<<", "<<aNew->Get4Momentum()<<G4endl;
532 532
533 theParticleChange.AddSecondary( aNew, fS 533 theParticleChange.AddSecondary( aNew, fSecID );
534 delete ddktv->operator[](i); 534 delete ddktv->operator[](i);
535 } 535 }
536 delete ddktv; 536 delete ddktv;
537 } 537 }
538 // recoil nucleus 538 // recoil nucleus
539 539
540 G4double eRecoil = sqrt( rM*rM + dP*dP ); 540 G4double eRecoil = sqrt( rM*rM + dP*dP );
541 fTr = eRecoil - rM; 541 fTr = eRecoil - rM;
542 G4ThreeVector vRecoil(dP*dX); 542 G4ThreeVector vRecoil(dP*dX);
543 // dP += G4UniformRand()*10.*MeV; 543 // dP += G4UniformRand()*10.*MeV;
544 G4LorentzVector rec4v(vRecoil, 0.); 544 G4LorentzVector rec4v(vRecoil, 0.);
545 rec4v.boost(bst); // back to lab 545 rec4v.boost(bst); // back to lab
546 fLVt += rec4v; 546 fLVt += rec4v;
547 const G4LorentzVector lvTarg = fLVt; // (vRe 547 const G4LorentzVector lvTarg = fLVt; // (vRecoil, eRecoil);
548 548
549 549
550 if( fRecoil ) // proton*? 550 if( fRecoil ) // proton*?
551 { 551 {
552 G4double grM = G4NucleiProperties::GetNucl 552 G4double grM = G4NucleiProperties::GetNuclearMass(A,Z);
553 G4double exE = fLVt.m() - grM; 553 G4double exE = fLVt.m() - grM;
554 if( exE < 5.*MeV ) exE = 5.*MeV + G4Unifor 554 if( exE < 5.*MeV ) exE = 5.*MeV + G4UniformRand()*10.*MeV;
555 555
556 const G4LorentzVector in4v( G4ThreeVector( 556 const G4LorentzVector in4v( G4ThreeVector( 0., 0., 0.), grM );
557 G4Fragment fragment( A, Z, in4v); // lvTar 557 G4Fragment fragment( A, Z, in4v); // lvTarg );
558 fragment.SetNumberOfHoles(1); 558 fragment.SetNumberOfHoles(1);
559 fragment.SetExcEnergyAndMomentum( exE, lvT 559 fragment.SetExcEnergyAndMomentum( exE, lvTarg );
560 560
561 RecoilDeexcitation(fragment); 561 RecoilDeexcitation(fragment);
562 } 562 }
563 else // momentum? 563 else // momentum?
564 { 564 {
565 theParticleChange.SetLocalEnergyDeposit( f 565 theParticleChange.SetLocalEnergyDeposit( fTr ); // eRecoil );
566 } 566 }
567 } 567 }
568 568
569 569
570 ////////////////////////////////////////////// 570 ///////////////////////////////////////////////////////
571 // 571 //
572 // Get final particles from excited recoil nuc 572 // Get final particles from excited recoil nucleus and write them to theParticleChange, delete the particle vector
573 573
574 void G4NeutrinoNucleusModel::RecoilDeexcitatio 574 void G4NeutrinoNucleusModel::RecoilDeexcitation( G4Fragment& fragment)
575 { 575 {
576 G4ReactionProductVector* products = fPreComp 576 G4ReactionProductVector* products = fPreCompound->DeExcite(fragment);
577 577
578 if( products != nullptr ) 578 if( products != nullptr )
579 { 579 {
580 for( auto & prod : *products ) // prod is 580 for( auto & prod : *products ) // prod is the pointer to final hadronic particle
581 { 581 {
582 theParticleChange.AddSecondary(new G4Dyn 582 theParticleChange.AddSecondary(new G4DynamicParticle( prod->GetDefinition(),
583 583 prod->GetTotalEnergy(),
584 584 prod->GetMomentum() ), fSecID );
585 delete prod; 585 delete prod;
586 } 586 }
587 delete products; 587 delete products;
588 } 588 }
589 } 589 }
590 590
591 /////////////////////////////////////////// 591 ///////////////////////////////////////////
592 // 592 //
593 // Fragmentation of lvX directly to pion and r 593 // Fragmentation of lvX directly to pion and recoil nucleus (A,Z): fLVh + fLVt -> pi + A*
594 594
595 void G4NeutrinoNucleusModel::CoherentPion( G4L 595 void G4NeutrinoNucleusModel::CoherentPion( G4LorentzVector & lvP, G4int pdgP, G4Nucleus & targetNucleus)
596 { 596 {
597 G4int A(0), Z(0), pdg = pdgP; 597 G4int A(0), Z(0), pdg = pdgP;
598 fLVcpi = G4LorentzVector(0.,0.,0.,0.); 598 fLVcpi = G4LorentzVector(0.,0.,0.,0.);
599 599
600 G4double rM(0.), mN(938.), det(0.), det2(0. 600 G4double rM(0.), mN(938.), det(0.), det2(0.);
601 G4double mI(0.); 601 G4double mI(0.);
602 mN = G4ParticleTable::GetParticleTable()->Fi 602 mN = G4ParticleTable::GetParticleTable()->FindParticle(2212)->GetPDGMass(); // *0.85; // *0.9; //
603 603
604 // mN = 1.*139.57 + G4UniformRand()*(938. - 604 // mN = 1.*139.57 + G4UniformRand()*(938. - 1.*139.57);
605 605
606 G4ThreeVector vN = lvP.boostVector(), bst(0. 606 G4ThreeVector vN = lvP.boostVector(), bst(0.,0.,0.);
607 // G4double gN = lvP.e()/lvP.m(); 607 // G4double gN = lvP.e()/lvP.m();
608 // G4LorentzVector lvNu(vN*gN*mN, mN*gN); 608 // G4LorentzVector lvNu(vN*gN*mN, mN*gN);
609 G4LorentzVector lvNu(0.,0.,0., mN); // lvN 609 G4LorentzVector lvNu(0.,0.,0., mN); // lvNu(bst, mN);
610 lvP.boost(-vN); // 9-3-20 610 lvP.boost(-vN); // 9-3-20
611 lvP = lvP - lvNu; // 9-3-20 already 1pi 611 lvP = lvP - lvNu; // 9-3-20 already 1pi
612 lvP.boost(vN); // 9-3-20 612 lvP.boost(vN); // 9-3-20
613 lvNu.boost(vN); // 9-3-20 613 lvNu.boost(vN); // 9-3-20
614 614
615 // G4cout<<vN-lvP.boostVector()<<", "; 615 // G4cout<<vN-lvP.boostVector()<<", ";
616 616
617 Z = targetNucleus.GetZ_asInt(); 617 Z = targetNucleus.GetZ_asInt();
618 A = targetNucleus.GetA_asInt(); 618 A = targetNucleus.GetA_asInt();
619 rM = targetNucleus.AtomicMass(A,Z); //->Atom 619 rM = targetNucleus.AtomicMass(A,Z); //->AtomicMass(); //
620 620
621 // G4cout<<rM<<", "; 621 // G4cout<<rM<<", ";
622 // G4cout<<A<<", "; 622 // G4cout<<A<<", ";
623 623
624 if( A == 1 ) 624 if( A == 1 )
625 { 625 {
626 bst = vN; // lvNu.boostVector(); // 9-3-20 626 bst = vN; // lvNu.boostVector(); // 9-3-20
627 // mI = 0.; // 9-3-20 627 // mI = 0.; // 9-3-20
628 rM = mN; 628 rM = mN;
629 } 629 }
630 else 630 else
631 { 631 {
632 G4Nucleus targ(A-1,Z); 632 G4Nucleus targ(A-1,Z);
633 mI = targ.AtomicMass(A-1,Z); 633 mI = targ.AtomicMass(A-1,Z);
634 G4LorentzVector lvTar(0.,0.,0.,mI); 634 G4LorentzVector lvTar(0.,0.,0.,mI);
635 lvNu = lvNu + lvTar; 635 lvNu = lvNu + lvTar;
636 bst = lvNu.boostVector(); 636 bst = lvNu.boostVector();
637 // bst = fLVt.boostVector(); // to recoil 637 // bst = fLVt.boostVector(); // to recoil rest frame
638 // G4cout<<fLVt<<" "<<bst<<G4endl; 638 // G4cout<<fLVt<<" "<<bst<<G4endl;
639 } 639 }
640 lvP.boost(-bst); // 9-3-20 640 lvP.boost(-bst); // 9-3-20
641 fMr = G4ParticleTable::GetParticleTable()->F 641 fMr = G4ParticleTable::GetParticleTable()->FindParticle(pdg)->GetPDGMass();
642 G4double eX = lvP.e(); 642 G4double eX = lvP.e();
643 G4double mX = lvP.m(); 643 G4double mX = lvP.m();
644 // G4cout<<mX-fMr<<", "; 644 // G4cout<<mX-fMr<<", ";
645 G4ThreeVector dX = (lvP.vect()).unit(); 645 G4ThreeVector dX = (lvP.vect()).unit();
646 // G4cout<<dX<<", "; 646 // G4cout<<dX<<", ";
647 G4double pX = sqrt(eX*eX-mX*mX); 647 G4double pX = sqrt(eX*eX-mX*mX);
648 // G4cout<<pX<<", "; 648 // G4cout<<pX<<", ";
649 G4double sumE = eX + rM; 649 G4double sumE = eX + rM;
650 G4double B = sumE*sumE + rM*rM - fMr*fMr - p 650 G4double B = sumE*sumE + rM*rM - fMr*fMr - pX*pX;
651 G4double a = 4.*(sumE*sumE - pX*pX); 651 G4double a = 4.*(sumE*sumE - pX*pX);
652 G4double b = -4.*B*pX; 652 G4double b = -4.*B*pX;
653 G4double c = 4.*sumE*sumE*rM*rM - B*B; 653 G4double c = 4.*sumE*sumE*rM*rM - B*B;
654 det2 = b*b-4.*a*c; 654 det2 = b*b-4.*a*c;
655 if(det2 > 0.) det = sqrt(det2); 655 if(det2 > 0.) det = sqrt(det2);
656 G4double dP = 0.5*(-b - det )/a; 656 G4double dP = 0.5*(-b - det )/a;
657 657
658 // dP = FinalMomentum( mI, rM, fMr, lvP); 658 // dP = FinalMomentum( mI, rM, fMr, lvP);
659 dP = FinalMomentum( rM, rM, fMr, lvP); // 9- 659 dP = FinalMomentum( rM, rM, fMr, lvP); // 9-3-20
660 660
661 // G4cout<<dP<<", "; 661 // G4cout<<dP<<", ";
662 pX -= dP; 662 pX -= dP;
663 if( pX < 0. ) pX = 0.; 663 if( pX < 0. ) pX = 0.;
664 664
665 eX = sqrt( dP*dP + fMr*fMr ); 665 eX = sqrt( dP*dP + fMr*fMr );
666 G4LorentzVector lvN( dP*dX, eX ); 666 G4LorentzVector lvN( dP*dX, eX );
667 667
668 if( A >= 1 ) lvN.boost(bst); // 9-3-20 back 668 if( A >= 1 ) lvN.boost(bst); // 9-3-20 back to lab
669 669
670 fLVcpi = lvN; 670 fLVcpi = lvN;
671 671
672 G4ParticleDefinition* pd2 = G4ParticleTable: 672 G4ParticleDefinition* pd2 = G4ParticleTable::GetParticleTable()->FindParticle(pdg);
673 G4DynamicParticle* dp2 = new G4DynamicPar 673 G4DynamicParticle* dp2 = new G4DynamicParticle( pd2, lvN);
674 theParticleChange.AddSecondary( dp2, fSecID 674 theParticleChange.AddSecondary( dp2, fSecID ); // coherent pion
675 675
676 // recoil nucleus 676 // recoil nucleus
677 677
678 G4double eRecoil = sqrt( rM*rM + pX*pX ); 678 G4double eRecoil = sqrt( rM*rM + pX*pX );
679 G4ThreeVector vRecoil(pX*dX); 679 G4ThreeVector vRecoil(pX*dX);
680 G4LorentzVector lvTarg1( vRecoil, eRecoil); 680 G4LorentzVector lvTarg1( vRecoil, eRecoil);
681 lvTarg1.boost(bst); 681 lvTarg1.boost(bst);
682 682
683 const G4LorentzVector lvTarg = lvTarg1; 683 const G4LorentzVector lvTarg = lvTarg1;
684 684
685 if( A > 1 ) // recoil target nucleus* 685 if( A > 1 ) // recoil target nucleus*
686 { 686 {
687 G4double grM = G4NucleiProperties::GetNucl 687 G4double grM = G4NucleiProperties::GetNuclearMass(A,Z);
688 G4double exE = fLVt.m() - grM; 688 G4double exE = fLVt.m() - grM;
689 689
690 if( exE < 5.*MeV ) exE = 5.*MeV + G4Unifor 690 if( exE < 5.*MeV ) exE = 5.*MeV + G4UniformRand()*10.*MeV; // vmg???
691 691
692 const G4LorentzVector in4v( G4ThreeVector( 692 const G4LorentzVector in4v( G4ThreeVector( 0., 0., 0.), grM );
693 G4Fragment fragment( A, Z, in4v); // lvTar 693 G4Fragment fragment( A, Z, in4v); // lvTarg );
694 fragment.SetNumberOfHoles(1); 694 fragment.SetNumberOfHoles(1);
695 fragment.SetExcEnergyAndMomentum( exE, lvT 695 fragment.SetExcEnergyAndMomentum( exE, lvTarg );
696 696
697 RecoilDeexcitation(fragment); 697 RecoilDeexcitation(fragment);
698 } 698 }
699 else // recoil target proton 699 else // recoil target proton
700 { 700 {
701 G4double eTkin = eRecoil - rM; 701 G4double eTkin = eRecoil - rM;
702 G4double eTh = 0.01*MeV; // 10.*MeV; 702 G4double eTh = 0.01*MeV; // 10.*MeV;
703 703
704 if( eTkin > eTh ) 704 if( eTkin > eTh )
705 { 705 {
706 G4DynamicParticle * aSec = new G4Dynamic 706 G4DynamicParticle * aSec = new G4DynamicParticle( G4Proton::Proton(), lvTarg);
707 theParticleChange.AddSecondary(aSec, fSe 707 theParticleChange.AddSecondary(aSec, fSecID);
708 } 708 }
709 else theParticleChange.SetLocalEnergyDepos 709 else theParticleChange.SetLocalEnergyDeposit( eTkin );
710 } 710 }
711 return; 711 return;
712 } 712 }
713 713
714 714
715 ////////////////////////////////////////////// 715 ////////////////////////////////////////////////////////////
716 // 716 //
717 // Excited barion decay to meson and barion, 717 // Excited barion decay to meson and barion,
718 // mass distributions and charge exchange are 718 // mass distributions and charge exchange are free parameters
719 719
720 void G4NeutrinoNucleusModel::ClusterDecay( G4L 720 void G4NeutrinoNucleusModel::ClusterDecay( G4LorentzVector & lvX, G4int qX)
721 { 721 {
722 G4bool finB = false; 722 G4bool finB = false;
723 G4int pdgB(0), i(0), qM(0), qB(0); // p 723 G4int pdgB(0), i(0), qM(0), qB(0); // pdgM(0),
724 G4double mM(0.), mB(0.), eM(0.), eB(0.), pM( 724 G4double mM(0.), mB(0.), eM(0.), eB(0.), pM(0.), pB(0.);
725 G4double mm1(0.), mm22(0.), M1(0.), M2(0.), 725 G4double mm1(0.), mm22(0.), M1(0.), M2(0.), mX(0.);
726 726
727 mX = lvX.m(); 727 mX = lvX.m();
728 728
729 G4double mN = G4ParticleTable::GetParticleT 729 G4double mN = G4ParticleTable::GetParticleTable()->FindParticle(2112)->GetPDGMass();
730 G4double mPi = G4ParticleTable::GetParticleT 730 G4double mPi = G4ParticleTable::GetParticleTable()->FindParticle(211)->GetPDGMass();
731 731
732 // G4double deltaM = 1.*MeV; // 30.*MeV; / 732 // G4double deltaM = 1.*MeV; // 30.*MeV; // 10.*MeV; // 100.*MeV; // 20.*MeV; //
733 G4double deltaMr[4] = { 0.*MeV, 0.*MeV, 100 733 G4double deltaMr[4] = { 0.*MeV, 0.*MeV, 100.*MeV, 0.*MeV};
734 734
735 G4ThreeVector dir(0.,0.,0.); 735 G4ThreeVector dir(0.,0.,0.);
736 G4ThreeVector bst(0.,0.,0.); 736 G4ThreeVector bst(0.,0.,0.);
737 G4LorentzVector lvM(0.,0.,0.,0.); 737 G4LorentzVector lvM(0.,0.,0.,0.);
738 G4LorentzVector lvB(0.,0.,0.,0.); 738 G4LorentzVector lvB(0.,0.,0.,0.);
739 739
740 for( i = 0; i < fClustNumber; ++i) // check 740 for( i = 0; i < fClustNumber; ++i) // check resonance
741 { 741 {
742 if( mX >= fBarMass[i] ) 742 if( mX >= fBarMass[i] )
743 { 743 {
744 pdgB = fBarPDG[i]; 744 pdgB = fBarPDG[i];
745 // mB = G4ParticleTable::GetParticleTa 745 // mB = G4ParticleTable::GetParticleTable()->FindParticle(pdgB)->GetPDGMass();
746 break; 746 break;
747 } 747 }
748 } 748 }
749 if( i == fClustNumber || i == fClustNumber-1 749 if( i == fClustNumber || i == fClustNumber-1 ) // low mass, p || n
750 { 750 {
751 if ( qX == 2 || qX == 0) { pdgB = 221 751 if ( qX == 2 || qX == 0) { pdgB = 2212; qB = 1;} // p for 2, 0
752 752
753 else if( qX == 1 || qX == -1) { pdgB = 211 753 else if( qX == 1 || qX == -1) { pdgB = 2112; qB = 0;} // n for 1, -1
754 754
755 return FinalBarion( lvX, qB, pdgB); 755 return FinalBarion( lvX, qB, pdgB);
756 } 756 }
757 else if( mX < fBarMass[i] + deltaMr[i] || m 757 else if( mX < fBarMass[i] + deltaMr[i] || mX < mN + mPi )
758 { 758 {
759 finB = true; // final barion -> out 759 finB = true; // final barion -> out
760 760
761 if ( qX == 1 && pdgB != 2212) pdgB = 761 if ( qX == 1 && pdgB != 2212) pdgB = pdgB - 10;
762 else if( qX == 0 && pdgB != 2212) pdgB = 762 else if( qX == 0 && pdgB != 2212) pdgB = pdgB - 110;
763 else if( qX == 0 && pdgB == 2212) pdgB = 763 else if( qX == 0 && pdgB == 2212) pdgB = pdgB - 100;
764 764
765 if( finB ) return FinalBarion( lvX, qX, pd 765 if( finB ) return FinalBarion( lvX, qX, pdgB ); // out
766 } 766 }
767 // no barion resonance, try 1->2 decay in CO 767 // no barion resonance, try 1->2 decay in COM frame
768 768
769 // try meson mass 769 // try meson mass
770 770
771 mm1 = mPi + 1.*MeV; // pi+ 771 mm1 = mPi + 1.*MeV; // pi+
772 mm22 = mX - mN; // mX-n 772 mm22 = mX - mN; // mX-n
773 773
774 if( mm22 <= mm1 ) // out with p or n 774 if( mm22 <= mm1 ) // out with p or n
775 { 775 {
776 if( qX == 2 || qX == 0) { pdgB = 221 776 if( qX == 2 || qX == 0) { pdgB = 2212; qB = 1;} // p
777 else if( qX == 1 || qX == -1) { pdgB = 211 777 else if( qX == 1 || qX == -1) { pdgB = 2112; qB = 0;} // n
778 778
779 return FinalBarion(lvX, qB, pdgB); 779 return FinalBarion(lvX, qB, pdgB);
780 } 780 }
781 else // try decay -> meson(cluster) + barion 781 else // try decay -> meson(cluster) + barion(cluster)
782 { 782 {
783 // G4double sigmaM = 50.*MeV; // 100.*MeV; 783 // G4double sigmaM = 50.*MeV; // 100.*MeV; // 200.*MeV; // 400.*MeV; // 800.*MeV; //
784 G4double rand = G4UniformRand(); 784 G4double rand = G4UniformRand();
785 785
786 // mM = mm1*mm22/( mm1 + rand*(mm22 - mm1) 786 // mM = mm1*mm22/( mm1 + rand*(mm22 - mm1) );
787 // mM = mm1*mm22/sqrt( mm1*mm1 + rand*(mm2 787 // mM = mm1*mm22/sqrt( mm1*mm1 + rand*(mm22*mm22 - mm1*mm1) );
788 // mM = -sigmaM*log( (1.- rand)*exp(-mm22/ 788 // mM = -sigmaM*log( (1.- rand)*exp(-mm22/sigmaM) + rand*exp(-mm1/sigmaM) );
789 mM = mm1 + rand*(mm22-mm1); 789 mM = mm1 + rand*(mm22-mm1);
790 790
791 791
792 for( i = 0; i < fClustNumber; ++i) 792 for( i = 0; i < fClustNumber; ++i)
793 { 793 {
794 if( mM >= fMesMass[i] ) 794 if( mM >= fMesMass[i] )
795 { 795 {
796 // pdgM = fMesPDG[i]; 796 // pdgM = fMesPDG[i];
797 // mM = G4ParticleTable::GetParticleTa 797 // mM = G4ParticleTable::GetParticleTable()->FindParticle(pdgM)->GetPDGMass();
798 break; 798 break;
799 } 799 }
800 } 800 }
801 M1 = G4ParticleTable::GetParticleTable()-> 801 M1 = G4ParticleTable::GetParticleTable()->FindParticle(2112)->GetPDGMass()+2.*MeV; // n
802 M2 = mX - mM; 802 M2 = mX - mM;
803 803
804 if( M2 <= M1 ) // 804 if( M2 <= M1 ) //
805 { 805 {
806 if ( qX == 2 || qX == 0) { pdgB = 2 806 if ( qX == 2 || qX == 0) { pdgB = 2212; qB = 1;} // p
807 else if( qX == 1 || qX == -1) { pdgB = 2 807 else if( qX == 1 || qX == -1) { pdgB = 2112; qB = 0;} // n
808 808
809 return FinalBarion(lvX, qB, pdgB); 809 return FinalBarion(lvX, qB, pdgB);
810 } 810 }
811 mB = M1 + G4UniformRand()*(M2-M1); 811 mB = M1 + G4UniformRand()*(M2-M1);
812 // mB = -sigmaM*log( (1.- rand)*exp(-M2/si 812 // mB = -sigmaM*log( (1.- rand)*exp(-M2/sigmaM) + rand*exp(-M1/sigmaM) );
813 813
814 bst = lvX.boostVector(); 814 bst = lvX.boostVector();
815 815
816 // dir = G4RandomDirection(); // ??? 816 // dir = G4RandomDirection(); // ???
817 // dir = G4ThreeVector(0.,0.,1.); 817 // dir = G4ThreeVector(0.,0.,1.);
818 dir = bst.orthogonal().unit(); // ?? 818 dir = bst.orthogonal().unit(); // ??
819 // G4double cost = exp(-G4UniformRand()); 819 // G4double cost = exp(-G4UniformRand());
820 // G4double sint = sqrt((1.-cost)*(1.+cost 820 // G4double sint = sqrt((1.-cost)*(1.+cost));
821 // G4double phi = twopi*G4UniformRand(); 821 // G4double phi = twopi*G4UniformRand();
822 // dir = G4ThreeVector(sint*cos(phi), sint 822 // dir = G4ThreeVector(sint*cos(phi), sint*sin(phi), cost);
823 823
824 eM = 0.5*(mX*mX + mM*mM - mB*mB)/mX; 824 eM = 0.5*(mX*mX + mM*mM - mB*mB)/mX;
825 pM = sqrt(eM*eM - mM*mM); 825 pM = sqrt(eM*eM - mM*mM);
826 lvM = G4LorentzVector( pM*dir, eM); 826 lvM = G4LorentzVector( pM*dir, eM);
827 lvM.boost(bst); 827 lvM.boost(bst);
828 828
829 eB = 0.5*(mX*mX + mB*mB - mM*mM)/mX; 829 eB = 0.5*(mX*mX + mB*mB - mM*mM)/mX;
830 pB = sqrt(eB*eB - mB*mB); 830 pB = sqrt(eB*eB - mB*mB);
831 lvB = G4LorentzVector(-pB*dir, eB); 831 lvB = G4LorentzVector(-pB*dir, eB);
832 lvB.boost(bst); 832 lvB.boost(bst);
833 833
834 // G4cout<<mM<<"/"<<mB<<", "; 834 // G4cout<<mM<<"/"<<mB<<", ";
835 835
836 // charge exchange 836 // charge exchange
837 837
838 if ( qX == 2 ) { qM = 1; qB = 1;} 838 if ( qX == 2 ) { qM = 1; qB = 1;}
839 else if( qX == 1 ) { qM = 0; qB = 1;} 839 else if( qX == 1 ) { qM = 0; qB = 1;}
840 else if( qX == 0 ) { qM = 0; qB = 0;} 840 else if( qX == 0 ) { qM = 0; qB = 0;}
841 else if( qX == -1 ) { qM = -1; qB = 0;} 841 else if( qX == -1 ) { qM = -1; qB = 0;}
842 842
843 // if ( qM == 0 ) pdgM = pdgM - 100; 843 // if ( qM == 0 ) pdgM = pdgM - 100;
844 // else if( qM == -1 ) pdgM = -pdgM; 844 // else if( qM == -1 ) pdgM = -pdgM;
845 845
846 MesonDecay( lvM, qM); // pdgM ); // 846 MesonDecay( lvM, qM); // pdgM ); //
847 847
848 // else 848 // else
849 ClusterDecay( lvB, qB ); // continue 849 ClusterDecay( lvB, qB ); // continue
850 } 850 }
851 } 851 }
852 852
853 853
854 ////////////////////////////////////////////// 854 ////////////////////////////////////////////////////////////
855 // 855 //
856 // Excited barion decay to meson and barion, 856 // Excited barion decay to meson and barion,
857 // mass distributions and charge exchange are 857 // mass distributions and charge exchange are free parameters
858 858
859 void G4NeutrinoNucleusModel::MesonDecay( G4Lor 859 void G4NeutrinoNucleusModel::MesonDecay( G4LorentzVector & lvX, G4int qX)
860 { 860 {
861 G4bool finB = false; 861 G4bool finB = false;
862 G4int pdgM(0), pdgB(0), i(0), qM(0), qB(0 862 G4int pdgM(0), pdgB(0), i(0), qM(0), qB(0);
863 G4double mM(0.), mB(0.), eM(0.), eB(0.), pM( 863 G4double mM(0.), mB(0.), eM(0.), eB(0.), pM(0.), pB(0.);
864 G4double mm1(0.), mm22(0.), M1(0.), M2(0.), 864 G4double mm1(0.), mm22(0.), M1(0.), M2(0.), mX(0.), Tkin(0.);
865 865
866 mX = lvX.m(); 866 mX = lvX.m();
867 Tkin = lvX.e() - mX; 867 Tkin = lvX.e() - mX;
868 868
869 // if( mX < 1120*MeV && mX > 1020*MeV ) // p 869 // if( mX < 1120*MeV && mX > 1020*MeV ) // phi(1020)->K+K-
870 if( mX < 1080*MeV && mX > 990*MeV && Tkin < 870 if( mX < 1080*MeV && mX > 990*MeV && Tkin < 600*MeV ) // phi(1020)->K+K-
871 { 871 {
872 return FinalMeson( lvX, qB, 333); 872 return FinalMeson( lvX, qB, 333);
873 } 873 }
874 G4double mPi = G4ParticleTable::GetParticleT 874 G4double mPi = G4ParticleTable::GetParticleTable()->FindParticle(211)->GetPDGMass();
875 875
876 G4double deltaMr[4] = { 0.*MeV, 0.*MeV, 100 876 G4double deltaMr[4] = { 0.*MeV, 0.*MeV, 100.*MeV, 0.*MeV};
877 877
878 G4ThreeVector dir(0.,0.,0.); 878 G4ThreeVector dir(0.,0.,0.);
879 G4ThreeVector bst(0.,0.,0.); 879 G4ThreeVector bst(0.,0.,0.);
880 G4LorentzVector lvM(0.,0.,0.,0.); 880 G4LorentzVector lvM(0.,0.,0.,0.);
881 G4LorentzVector lvB(0.,0.,0.,0.); 881 G4LorentzVector lvB(0.,0.,0.,0.);
882 882
883 for( i = 0; i < fClustNumber; ++i) // check 883 for( i = 0; i < fClustNumber; ++i) // check resonance
884 { 884 {
885 if( mX >= fMesMass[i] ) 885 if( mX >= fMesMass[i] )
886 { 886 {
887 pdgB = fMesPDG[i]; 887 pdgB = fMesPDG[i];
888 // mB = G4ParticleTable::GetParticleTa 888 // mB = G4ParticleTable::GetParticleTable()->FindParticle(pdgB)->GetPDGMass();
889 break; 889 break;
890 } 890 }
891 } 891 }
892 if( i == fClustNumber ) // || i == fClustNum 892 if( i == fClustNumber ) // || i == fClustNumber-1 ) // low mass, p || n
893 { 893 {
894 if ( qX == 1) { pdgB = 211; qB = 1; 894 if ( qX == 1) { pdgB = 211; qB = 1;} // pi+
895 else if( qX == 0 ) { pdgB = 111; qB = 0; 895 else if( qX == 0 ) { pdgB = 111; qB = 0;} // pi0
896 else if( qX == -1) { pdgB = -211; qB = -1; 896 else if( qX == -1) { pdgB = -211; qB = -1;} // pi-
897 897
898 return FinalMeson( lvX, qB, pdgB); 898 return FinalMeson( lvX, qB, pdgB);
899 } 899 }
900 else if( mX < fMesMass[i] + deltaMr[i] ) // 900 else if( mX < fMesMass[i] + deltaMr[i] ) // || mX < mPi + mPi ) //
901 { 901 {
902 finB = true; // final barion -> out 902 finB = true; // final barion -> out
903 pdgB = fMesPDG[i]; 903 pdgB = fMesPDG[i];
904 904
905 // if ( qX == 1 && pdgB != 2212) pdg 905 // if ( qX == 1 && pdgB != 2212) pdgB = pdgB - 10;
906 906
907 if( qX == 0 ) pdgB = pdgB - 100; 907 if( qX == 0 ) pdgB = pdgB - 100;
908 else if( qX == -1 ) pdgB = -pdgB; 908 else if( qX == -1 ) pdgB = -pdgB;
909 909
910 if( finB ) return FinalMeson( lvX, qX, pdg 910 if( finB ) return FinalMeson( lvX, qX, pdgB ); // out
911 } 911 }
912 // no resonance, try 1->2 decay in COM frame 912 // no resonance, try 1->2 decay in COM frame
913 913
914 // try meson 914 // try meson
915 915
916 mm1 = mPi + 1.*MeV; // pi+ 916 mm1 = mPi + 1.*MeV; // pi+
917 mm22 = mX - mPi - 1.*MeV; // mX-n 917 mm22 = mX - mPi - 1.*MeV; // mX-n
918 918
919 if( mm22 <= mm1 ) // out 919 if( mm22 <= mm1 ) // out
920 { 920 {
921 if ( qX == 1) { pdgB = 211; qB = 1;} 921 if ( qX == 1) { pdgB = 211; qB = 1;} // pi+
922 else if( qX == 0 ) { pdgB = 111; qB = 0;} 922 else if( qX == 0 ) { pdgB = 111; qB = 0;} // pi0
923 else if( qX == -1) { pdgB = -211; qB = -1; 923 else if( qX == -1) { pdgB = -211; qB = -1;} // pi-
924 924
925 return FinalMeson(lvX, qB, pdgB); 925 return FinalMeson(lvX, qB, pdgB);
926 } 926 }
927 else // try decay -> pion + meson(cluster) 927 else // try decay -> pion + meson(cluster)
928 { 928 {
929 // G4double sigmaM = 50.*MeV; // 100.*MeV; 929 // G4double sigmaM = 50.*MeV; // 100.*MeV; // 200.*MeV; // 400.*MeV; // 800.*MeV; //
930 G4double rand = G4UniformRand(); 930 G4double rand = G4UniformRand();
931 931
932 if ( qX == 1 ) { qM = 1; qB = 0;} 932 if ( qX == 1 ) { qM = 1; qB = 0;}
933 else if( qX == 0 ) { qM = -1; qB = 1;} / 933 else if( qX == 0 ) { qM = -1; qB = 1;} // { qM = 0; qB = 0;} //
934 else if( qX == -1 ) { qM = -1; qB = 0;} 934 else if( qX == -1 ) { qM = -1; qB = 0;}
935 /* 935 /*
936 mM = mPi; 936 mM = mPi;
937 if(qM == 0) mM = G4ParticleTable::GetParti 937 if(qM == 0) mM = G4ParticleTable::GetParticleTable()->FindParticle(111)->GetPDGMass(); //pi0
938 pdgM = fMesPDG[fClustNumber-1]; 938 pdgM = fMesPDG[fClustNumber-1];
939 */ 939 */
940 // mm1*mm22/( mm1 + rand*(mm22 - mm1) ); 940 // mm1*mm22/( mm1 + rand*(mm22 - mm1) );
941 // mM = mm1*mm22/sqrt( mm1*mm1 + rand*(mm2 941 // mM = mm1*mm22/sqrt( mm1*mm1 + rand*(mm22*mm22 - mm1*mm1) );
942 // mM = -sigmaM*log( (1.- rand)*exp(-mm22/ 942 // mM = -sigmaM*log( (1.- rand)*exp(-mm22/sigmaM) + rand*exp(-mm1/sigmaM) );
943 mM = mm1 + rand*(mm22-mm1); 943 mM = mm1 + rand*(mm22-mm1);
944 // mM = mm1 + 0.9*(mm22-mm1); 944 // mM = mm1 + 0.9*(mm22-mm1);
945 945
946 946
947 for( i = 0; i < fClustNumber; ++i) 947 for( i = 0; i < fClustNumber; ++i)
948 { 948 {
949 if( mM >= fMesMass[i] ) 949 if( mM >= fMesMass[i] )
950 { 950 {
951 pdgM = fMesPDG[i]; 951 pdgM = fMesPDG[i];
952 // mM = G4ParticleTable::GetParticleTa 952 // mM = G4ParticleTable::GetParticleTable()->FindParticle(pdgM)->GetPDGMass();
953 break; 953 break;
954 } 954 }
955 } 955 }
956 if( i == fClustNumber || i == fClustNumber 956 if( i == fClustNumber || i == fClustNumber-1 ) // low mass, p || n
957 { 957 {
958 if ( qX == 1) { pdgB = 211; qB = 1; 958 if ( qX == 1) { pdgB = 211; qB = 1;} // pi+
959 else if( qX == 0 ) { pdgB = 111; qB = 0; 959 else if( qX == 0 ) { pdgB = 111; qB = 0;} // pi0
960 else if( qX == -1) { pdgB = -211; qB = - 960 else if( qX == -1) { pdgB = -211; qB = -1;} // pi-
961 961
962 return FinalMeson( lvX, qB, pdgB); 962 return FinalMeson( lvX, qB, pdgB);
963 } 963 }
964 else if( mX < fMesMass[i] + deltaMr[i] ) / 964 else if( mX < fMesMass[i] + deltaMr[i] ) // || mX < mPi + mPi ) //
965 { 965 {
966 finB = true; // final barion -> out 966 finB = true; // final barion -> out
967 pdgB = fMesPDG[i]; 967 pdgB = fMesPDG[i];
968 968
969 // if ( qX == 1 && pdgB != 2212) pdg 969 // if ( qX == 1 && pdgB != 2212) pdgB = pdgB - 10;
970 970
971 if( qX == 0 ) pdgB = pdgB - 100; 971 if( qX == 0 ) pdgB = pdgB - 100;
972 else if( qX == -1 ) pdgB = -pdgB; 972 else if( qX == -1 ) pdgB = -pdgB;
973 973
974 if( finB ) return FinalMeson( lvX, qX, p 974 if( finB ) return FinalMeson( lvX, qX, pdgB ); // out
975 } 975 }
976 976
977 M1 = G4ParticleTable::GetParticleTable()-> 977 M1 = G4ParticleTable::GetParticleTable()->FindParticle(211)->GetPDGMass()+2.*MeV; // n
978 M2 = mX - mM; 978 M2 = mX - mM;
979 979
980 if( M2 <= M1 ) // 980 if( M2 <= M1 ) //
981 { 981 {
982 if ( qX == 1) { pdgB = 211; qB = 1; 982 if ( qX == 1) { pdgB = 211; qB = 1;} // pi+
983 else if( qX == 0 ) { pdgB = 111; qB = 0; 983 else if( qX == 0 ) { pdgB = 111; qB = 0;} // pi0
984 else if( qX == -1) { pdgB = -211; qB = - 984 else if( qX == -1) { pdgB = -211; qB = -1;} // pi-
985 985
986 return FinalMeson(lvX, qB, pdgB); 986 return FinalMeson(lvX, qB, pdgB);
987 } 987 }
988 mB = M1 + G4UniformRand()*(M2-M1); 988 mB = M1 + G4UniformRand()*(M2-M1);
989 // mB = -sigmaM*log( (1.- rand)*exp(-M2/si 989 // mB = -sigmaM*log( (1.- rand)*exp(-M2/sigmaM) + rand*exp(-M1/sigmaM) );
990 // mB = M1 + 0.9*(M2-M1); 990 // mB = M1 + 0.9*(M2-M1);
991 991
992 bst = lvX.boostVector(); 992 bst = lvX.boostVector();
993 993
994 // dir = G4RandomDirection(); 994 // dir = G4RandomDirection();
995 dir = bst.orthogonal().unit(); 995 dir = bst.orthogonal().unit();
996 996
997 eM = 0.5*(mX*mX + mM*mM - mB*mB)/mX; 997 eM = 0.5*(mX*mX + mM*mM - mB*mB)/mX;
998 pM = sqrt(eM*eM - mM*mM); 998 pM = sqrt(eM*eM - mM*mM);
999 lvM = G4LorentzVector( pM*dir, eM); 999 lvM = G4LorentzVector( pM*dir, eM);
1000 lvM.boost(bst); 1000 lvM.boost(bst);
1001 1001
1002 eB = 0.5*(mX*mX + mB*mB - mM*mM)/mX; 1002 eB = 0.5*(mX*mX + mB*mB - mM*mM)/mX;
1003 pB = sqrt(eB*eB - mB*mB); 1003 pB = sqrt(eB*eB - mB*mB);
1004 lvB = G4LorentzVector(-pB*dir, eB); 1004 lvB = G4LorentzVector(-pB*dir, eB);
1005 lvB.boost(bst); 1005 lvB.boost(bst);
1006 1006
1007 // G4cout<<mM<<"/"<<mB<<", "; 1007 // G4cout<<mM<<"/"<<mB<<", ";
1008 1008
1009 // charge exchange 1009 // charge exchange
1010 1010
1011 // if ( qX == 2 ) { qM = 1; qB = 1; 1011 // if ( qX == 2 ) { qM = 1; qB = 1;}
1012 1012
1013 if ( qM == 0 ) pdgM = pdgM - 100; 1013 if ( qM == 0 ) pdgM = pdgM - 100;
1014 else if( qM == -1 ) pdgM = -pdgM; 1014 else if( qM == -1 ) pdgM = -pdgM;
1015 1015
1016 MesonDecay( lvM, qM ); // 1016 MesonDecay( lvM, qM ); //
1017 1017
1018 MesonDecay( lvB, qB ); // continue 1018 MesonDecay( lvB, qB ); // continue
1019 } 1019 }
1020 } 1020 }
1021 1021
1022 ///////////////////////////////////////////// 1022 ///////////////////////////////////////////////////////////////////////
1023 // 1023 //
1024 // return final momentum x in the reaction lv 1024 // return final momentum x in the reaction lvX + mI -> mF + mP with momenta p-x, x
1025 1025
1026 G4double G4NeutrinoNucleusModel::FinalMomentu 1026 G4double G4NeutrinoNucleusModel::FinalMomentum(G4double mI, G4double mF, G4double mP, G4LorentzVector lvX)
1027 { 1027 {
1028 G4double result(0.), delta(0.); 1028 G4double result(0.), delta(0.);
1029 // G4double mI2 = mI*mI; 1029 // G4double mI2 = mI*mI;
1030 G4double mF2 = mF*mF; 1030 G4double mF2 = mF*mF;
1031 G4double mP2 = mP*mP; 1031 G4double mP2 = mP*mP;
1032 G4double eX = lvX.e(); 1032 G4double eX = lvX.e();
1033 // G4double mX = lvX.m(); 1033 // G4double mX = lvX.m();
1034 G4double pX = lvX.vect().mag(); 1034 G4double pX = lvX.vect().mag();
1035 G4double pX2 = pX*pX; 1035 G4double pX2 = pX*pX;
1036 G4double sI = eX + mI; 1036 G4double sI = eX + mI;
1037 G4double sI2 = sI*sI; 1037 G4double sI2 = sI*sI;
1038 G4double B = sI2 - mF2 -pX2 + mP2; 1038 G4double B = sI2 - mF2 -pX2 + mP2;
1039 G4double B2 = B*B; 1039 G4double B2 = B*B;
1040 G4double a = 4.*(sI2-pX2); 1040 G4double a = 4.*(sI2-pX2);
1041 G4double b = -4.*B*pX; 1041 G4double b = -4.*B*pX;
1042 G4double c = 4.*sI2*mP2 - B2; 1042 G4double c = 4.*sI2*mP2 - B2;
1043 G4double delta2 = b*b -4.*a*c; 1043 G4double delta2 = b*b -4.*a*c;
1044 1044
1045 if( delta2 >= 0. ) delta = sqrt(delta2); 1045 if( delta2 >= 0. ) delta = sqrt(delta2);
1046 1046
1047 result = 0.5*(-b-delta)/a; 1047 result = 0.5*(-b-delta)/a;
1048 // result = 0.5*(-b+delta)/a; 1048 // result = 0.5*(-b+delta)/a;
1049 1049
1050 return result; 1050 return result;
1051 } 1051 }
1052 1052
1053 ///////////////////////////////////////////// 1053 /////////////////////////////////////////////////////////////////
1054 // 1054 //
1055 // 1055 //
1056 1056
1057 G4double G4NeutrinoNucleusModel::FermiMomentu 1057 G4double G4NeutrinoNucleusModel::FermiMomentum( G4Nucleus & targetNucleus)
1058 { 1058 {
1059 G4int Z = targetNucleus.GetZ_asInt(); 1059 G4int Z = targetNucleus.GetZ_asInt();
1060 G4int A = targetNucleus.GetA_asInt(); 1060 G4int A = targetNucleus.GetA_asInt();
1061 1061
1062 G4double kF(250.*MeV); 1062 G4double kF(250.*MeV);
1063 G4double kp = 365.*MeV; 1063 G4double kp = 365.*MeV;
1064 G4double kn = 231.*MeV; 1064 G4double kn = 231.*MeV;
1065 G4double t1 = 0.479; 1065 G4double t1 = 0.479;
1066 G4double t2 = 0.526; 1066 G4double t2 = 0.526;
1067 G4double ZpA = G4double(Z)/G4double(A); 1067 G4double ZpA = G4double(Z)/G4double(A);
1068 G4double NpA = 1. - ZpA; 1068 G4double NpA = 1. - ZpA;
1069 1069
1070 if ( Z == 1 && A == 1 ) { kF = 0.; 1070 if ( Z == 1 && A == 1 ) { kF = 0.; } // hydrogen ???
1071 else if ( Z == 1 && A == 2 ) { kF = 87.* 1071 else if ( Z == 1 && A == 2 ) { kF = 87.*MeV; }
1072 else if ( Z == 2 && A == 3 ) { kF = 134. 1072 else if ( Z == 2 && A == 3 ) { kF = 134.*MeV; }
1073 else if ( Z == 6 && A == 12 ) { kF = 221. 1073 else if ( Z == 6 && A == 12 ) { kF = 221.*MeV; }
1074 else if ( Z == 14 && A == 28 ) { kF = 239. 1074 else if ( Z == 14 && A == 28 ) { kF = 239.*MeV; }
1075 else if ( Z == 26 && A == 56 ) { kF = 257. 1075 else if ( Z == 26 && A == 56 ) { kF = 257.*MeV; }
1076 else if ( Z == 82 && A == 208 ) { kF = 265. 1076 else if ( Z == 82 && A == 208 ) { kF = 265.*MeV; }
1077 else 1077 else
1078 { 1078 {
1079 kF = kp*ZpA*( 1 - pow( G4double(A), -t1 ) 1079 kF = kp*ZpA*( 1 - pow( G4double(A), -t1 ) ) + kn*NpA*( 1 - pow( G4double(A), -t2 ) );
1080 } 1080 }
1081 return kF; 1081 return kF;
1082 } 1082 }
1083 1083
1084 ///////////////////////////////////////////// 1084 /////////////////////////////////////////////////////////////////
1085 // 1085 //
1086 // sample nucleon momentum of Fermi motion fo 1086 // sample nucleon momentum of Fermi motion for 1p1h and 2p2h modes
1087 1087
1088 G4double G4NeutrinoNucleusModel::NucleonMomen 1088 G4double G4NeutrinoNucleusModel::NucleonMomentum( G4Nucleus & targetNucleus)
1089 { 1089 {
1090 G4int A = targetNucleus.GetA_asInt(); 1090 G4int A = targetNucleus.GetA_asInt();
1091 G4double kF = FermiMomentum( targetNucleus) 1091 G4double kF = FermiMomentum( targetNucleus);
1092 G4double mom(0.), kCut = 0.5*GeV; // kCut = 1092 G4double mom(0.), kCut = 0.5*GeV; // kCut = 1.*GeV; // kCut = 2.*GeV; // kCut = 4.*GeV; //
1093 // G4double cof = 2./GeV; 1093 // G4double cof = 2./GeV;
1094 // G4double ksi = kF*kF*cof*cof/pi/pi; 1094 // G4double ksi = kF*kF*cof*cof/pi/pi;
1095 G4double th = 1.; // 1. - 6.*ksi; // 1095 G4double th = 1.; // 1. - 6.*ksi; //
1096 1096
1097 if( G4UniformRand() < th || A < 3 ) // 1p1 1097 if( G4UniformRand() < th || A < 3 ) // 1p1h
1098 { 1098 {
1099 mom = kF*pow( G4UniformRand(), 1./3.); 1099 mom = kF*pow( G4UniformRand(), 1./3.);
1100 } 1100 }
1101 else // 2p2h 1101 else // 2p2h
1102 { 1102 {
1103 mom = kF*kCut; 1103 mom = kF*kCut;
1104 mom /= kCut - G4UniformRand()*(kCut - kF) 1104 mom /= kCut - G4UniformRand()*(kCut - kF);
1105 f2p2h = true; 1105 f2p2h = true;
1106 } 1106 }
1107 return mom; 1107 return mom;
1108 } 1108 }
1109 1109
1110 1110
1111 ///////////////////////////////////////////// 1111 //////////////////////////////////////////////////////
1112 // 1112 //
1113 // Excitation energy according Bodek 1113 // Excitation energy according Bodek
1114 1114
1115 G4double G4NeutrinoNucleusModel::GetEx( G4int 1115 G4double G4NeutrinoNucleusModel::GetEx( G4int A, G4bool fP )
1116 { 1116 {
1117 G4double eX(10.*MeV), a1(0.), a2(0.), e1(0. 1117 G4double eX(10.*MeV), a1(0.), a2(0.), e1(0.), e2(0.), aa = G4double(A);
1118 G4int i(0); 1118 G4int i(0);
1119 const G4int maxBin = 12; 1119 const G4int maxBin = 12;
1120 1120
1121 G4double refA[maxBin] = { 2., 6., 12., 16., 1121 G4double refA[maxBin] = { 2., 6., 12., 16., 27., 28., 40., 50., 56., 58., 197., 208. };
1122 1122
1123 G4double pEx[maxBin] = { 0., 12.2, 10.1, 10 1123 G4double pEx[maxBin] = { 0., 12.2, 10.1, 10.9, 21.6, 12.4, 17.8, 17., 19., 16.8, 19.5, 14.7 };
1124 1124
1125 G4double nEx[maxBin] = { 0., 12.2, 10., 10. 1125 G4double nEx[maxBin] = { 0., 12.2, 10., 10.2, 21.6, 12.4, 21.8, 17., 19., 16.8, 19.5, 16.9 };
1126 1126
1127 G4DataVector dE(12,0.); 1127 G4DataVector dE(12,0.);
1128 1128
1129 if(fP) for( i = 0; i < maxBin; ++i ) dE[i] 1129 if(fP) for( i = 0; i < maxBin; ++i ) dE[i] = pEx[i];
1130 else dE[i] 1130 else dE[i] = nEx[i];
1131 1131
1132 for( i = 0; i < maxBin; ++i ) 1132 for( i = 0; i < maxBin; ++i )
1133 { 1133 {
1134 if( aa <= refA[i] ) break; 1134 if( aa <= refA[i] ) break;
1135 } 1135 }
1136 if( i >= maxBin ) eX = dE[maxBin-1]; 1136 if( i >= maxBin ) eX = dE[maxBin-1];
1137 else if( i <= 0 ) eX = dE[0]; 1137 else if( i <= 0 ) eX = dE[0];
1138 else 1138 else
1139 { 1139 {
1140 a1 = refA[i-1]; 1140 a1 = refA[i-1];
1141 a2 = refA[i]; 1141 a2 = refA[i];
1142 e1 = dE[i-1]; 1142 e1 = dE[i-1];
1143 e2 = dE[i]; 1143 e2 = dE[i];
1144 if (a1 == a2 || e1 == e2 ) eX = dE[i]; 1144 if (a1 == a2 || e1 == e2 ) eX = dE[i];
1145 else eX = e1 + (e2-e1)*(aa-a1)/(a2-a1); 1145 else eX = e1 + (e2-e1)*(aa-a1)/(a2-a1);
1146 } 1146 }
1147 return eX; 1147 return eX;
1148 } 1148 }
1149 1149
1150 1150
1151 1151
1152 ///////////////////////////////////////////// 1152 ///////////////////////////////////////////////////////
1153 // 1153 //
1154 // Two G-function sampling for the nucleon mo 1154 // Two G-function sampling for the nucleon momentum
1155 1155
1156 G4double G4NeutrinoNucleusModel::GgSampleNM(G 1156 G4double G4NeutrinoNucleusModel::GgSampleNM(G4Nucleus & nucl)
1157 { 1157 {
1158 f2p2h = false; 1158 f2p2h = false;
1159 G4double /* distr(0.), tail(0.), */ shift(1 1159 G4double /* distr(0.), tail(0.), */ shift(1.), xx(1.), mom(0.), th(0.1);
1160 G4double kF = FermiMomentum( nucl); 1160 G4double kF = FermiMomentum( nucl);
1161 G4double momMax = 2.*kF; // 1.*GeV; // 1. 1161 G4double momMax = 2.*kF; // 1.*GeV; // 1.*GeV; //
1162 G4double aa = 5.5; 1162 G4double aa = 5.5;
1163 G4double ll = 6.0; // 6.5; // 1163 G4double ll = 6.0; // 6.5; //
1164 1164
1165 G4int A = nucl.GetA_asInt(); 1165 G4int A = nucl.GetA_asInt();
1166 1166
1167 if( A <= 12) th = 0.1; 1167 if( A <= 12) th = 0.1;
1168 else 1168 else
1169 { 1169 {
1170 // th = 0.1/(1.+log(G4double(A)/12.)); 1170 // th = 0.1/(1.+log(G4double(A)/12.));
1171 th = 1.2/( G4double(A) + 1.35*log(G4doubl 1171 th = 1.2/( G4double(A) + 1.35*log(G4double(A)/12.) );
1172 } 1172 }
1173 shift = 0.99; // 0.95; // 1173 shift = 0.99; // 0.95; //
1174 xx = mom/shift/kF; 1174 xx = mom/shift/kF;
1175 1175
1176 G4double rr = G4UniformRand(); 1176 G4double rr = G4UniformRand();
1177 1177
1178 if( rr > th ) 1178 if( rr > th )
1179 { 1179 {
1180 aa = 5.5; 1180 aa = 5.5;
1181 1181
1182 if( A <= 12 ) ll = 6.0; 1182 if( A <= 12 ) ll = 6.0;
1183 else 1183 else
1184 { 1184 {
1185 ll = 6.0 + 1.35*log(G4double(A)/12.); 1185 ll = 6.0 + 1.35*log(G4double(A)/12.);
1186 } 1186 }
1187 xx = RandGamma::shoot(aa,ll); 1187 xx = RandGamma::shoot(aa,ll);
1188 shift = 0.99; 1188 shift = 0.99;
1189 mom = xx*shift*kF; 1189 mom = xx*shift*kF;
1190 } 1190 }
1191 else 1191 else
1192 { 1192 {
1193 f2p2h = true; 1193 f2p2h = true;
1194 aa = 6.5; 1194 aa = 6.5;
1195 ll = 6.5; 1195 ll = 6.5;
1196 xx = RandGamma::shoot(aa,ll); 1196 xx = RandGamma::shoot(aa,ll);
1197 shift = 2.5; 1197 shift = 2.5;
1198 mom = xx*shift*kF; 1198 mom = xx*shift*kF;
1199 } 1199 }
1200 if( mom > momMax ) mom = G4UniformRand()*mo 1200 if( mom > momMax ) mom = G4UniformRand()*momMax;
1201 if( mom > 2.*kF ) f2p2h = true; 1201 if( mom > 2.*kF ) f2p2h = true;
1202 1202
1203 // mom = 0.; 1203 // mom = 0.;
1204 1204
1205 return mom; 1205 return mom;
1206 } 1206 }
1207 1207
1208 1208
1209 ///////////////////////////////////// experim 1209 ///////////////////////////////////// experimental arrays and get functions ////////////////////////////////////////
1210 // 1210 //
1211 // Return index of nu/anu energy array corres 1211 // Return index of nu/anu energy array corresponding to the neutrino energy
1212 1212
1213 G4int G4NeutrinoNucleusModel::GetEnergyIndex( 1213 G4int G4NeutrinoNucleusModel::GetEnergyIndex(G4double energy)
1214 { 1214 {
1215 G4int i, eIndex = 0; 1215 G4int i, eIndex = 0;
1216 1216
1217 for( i = 0; i < fIndex; i++) 1217 for( i = 0; i < fIndex; i++)
1218 { 1218 {
1219 if( energy <= fNuMuEnergy[i]*GeV ) 1219 if( energy <= fNuMuEnergy[i]*GeV )
1220 { 1220 {
1221 eIndex = i; 1221 eIndex = i;
1222 break; 1222 break;
1223 } 1223 }
1224 } 1224 }
1225 if( i >= fIndex ) eIndex = fIndex; 1225 if( i >= fIndex ) eIndex = fIndex;
1226 // G4cout<<"eIndex = "<<eIndex<<G4endl; 1226 // G4cout<<"eIndex = "<<eIndex<<G4endl;
1227 return eIndex; 1227 return eIndex;
1228 } 1228 }
1229 1229
1230 ///////////////////////////////////////////// 1230 /////////////////////////////////////////////////////
1231 // 1231 //
1232 // nu_mu QE/Tot ratio for index-1, index line 1232 // nu_mu QE/Tot ratio for index-1, index linear over energy
1233 1233
1234 G4double G4NeutrinoNucleusModel::GetNuMuQeTot 1234 G4double G4NeutrinoNucleusModel::GetNuMuQeTotRat(G4int index, G4double energy)
1235 { 1235 {
1236 G4double ratio(0.); 1236 G4double ratio(0.);
1237 // GetMinNuMuEnergy() 1237 // GetMinNuMuEnergy()
1238 if( index <= 0 || energy < fNuMuEnergy[0] ) 1238 if( index <= 0 || energy < fNuMuEnergy[0] ) ratio = 0.;
1239 else if (index >= fIndex) ratio = fNuMuQeTo 1239 else if (index >= fIndex) ratio = fNuMuQeTotRat[fIndex-1]*fOnePionEnergy[fIndex-1]*GeV/energy;
1240 else 1240 else
1241 { 1241 {
1242 G4double x1 = fNuMuEnergy[index-1]*GeV; 1242 G4double x1 = fNuMuEnergy[index-1]*GeV;
1243 G4double x2 = fNuMuEnergy[index]*GeV; 1243 G4double x2 = fNuMuEnergy[index]*GeV;
1244 G4double y1 = fNuMuQeTotRat[index-1]; 1244 G4double y1 = fNuMuQeTotRat[index-1];
1245 G4double y2 = fNuMuQeTotRat[index]; 1245 G4double y2 = fNuMuQeTotRat[index];
1246 1246
1247 if(x1 >= x2) return fNuMuQeTotRat[index]; 1247 if(x1 >= x2) return fNuMuQeTotRat[index];
1248 else 1248 else
1249 { 1249 {
1250 G4double angle = (y2-y1)/(x2-x1); 1250 G4double angle = (y2-y1)/(x2-x1);
1251 ratio = y1 + (energy-x1)*angle; 1251 ratio = y1 + (energy-x1)*angle;
1252 } 1252 }
1253 } 1253 }
1254 return ratio; 1254 return ratio;
1255 } 1255 }
1256 1256
1257 ///////////////////////////////////////////// 1257 ////////////////////////////////////////////////////////
1258 1258
1259 const G4double G4NeutrinoNucleusModel::fNuMuE 1259 const G4double G4NeutrinoNucleusModel::fNuMuEnergy[50] =
1260 { 1260 {
1261 0.112103, 0.117359, 0.123119, 0.129443, 0.1 1261 0.112103, 0.117359, 0.123119, 0.129443, 0.136404,
1262 0.144084, 0.152576, 0.161991, 0.172458, 0.1 1262 0.144084, 0.152576, 0.161991, 0.172458, 0.184126,
1263 0.197171, 0.211801, 0.228261, 0.24684, 0.26 1263 0.197171, 0.211801, 0.228261, 0.24684, 0.267887,
1264 0.291816, 0.319125, 0.350417, 0.386422, 0.4 1264 0.291816, 0.319125, 0.350417, 0.386422, 0.428032,
1265 0.47634, 0.532692, 0.598756, 0.676612, 0.76 1265 0.47634, 0.532692, 0.598756, 0.676612, 0.768868,
1266 0.878812, 1.01062, 1.16963, 1.36271, 1.5987 1266 0.878812, 1.01062, 1.16963, 1.36271, 1.59876,
1267 1.88943, 2.25002, 2.70086, 3.26916, 3.99166 1267 1.88943, 2.25002, 2.70086, 3.26916, 3.99166,
1268 4.91843, 6.11836, 7.6872, 9.75942, 12.5259, 1268 4.91843, 6.11836, 7.6872, 9.75942, 12.5259,
1269 16.2605, 21.3615, 28.4141, 38.2903, 52.3062 1269 16.2605, 21.3615, 28.4141, 38.2903, 52.3062,
1270 72.4763, 101.93, 145.6, 211.39, 312.172 1270 72.4763, 101.93, 145.6, 211.39, 312.172
1271 }; 1271 };
1272 1272
1273 ///////////////////////////////////////////// 1273 ////////////////////////////////////////////////////////
1274 1274
1275 const G4double G4NeutrinoNucleusModel::fNuMuQ 1275 const G4double G4NeutrinoNucleusModel::fNuMuQeTotRat[50] =
1276 { 1276 {
1277 // 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1277 // 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
1278 // 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1278 // 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
1279 // 1., 1., 1., 0.982311, 1279 // 1., 1., 1., 0.982311,
1280 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0 1280 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98,
1281 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0 1281 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98,
1282 0.97, 0.96, 0.95, 0.93, 1282 0.97, 0.96, 0.95, 0.93,
1283 0.917794, 0.850239, 0.780412, 0.709339, 0.6 1283 0.917794, 0.850239, 0.780412, 0.709339, 0.638134, 0.568165,
1284 0.500236, 0.435528, 0.375015, 0.319157, 0.2 1284 0.500236, 0.435528, 0.375015, 0.319157, 0.268463, 0.2232, 0.183284,
1285 0.148627, 0.119008, 0.0940699, 0.0733255, 0 1285 0.148627, 0.119008, 0.0940699, 0.0733255, 0.0563819, 0.0427312, 0.0319274,
1286 0.0235026, 0.0170486, 0.0122149, 0.00857825 1286 0.0235026, 0.0170486, 0.0122149, 0.00857825, 0.00594018, 0.00405037
1287 }; 1287 };
1288 1288
1289 ///////////////////////////////////////////// 1289 /////////////////////////////////////////////////////
1290 // 1290 //
1291 // Return index of one pion array correspondi 1291 // Return index of one pion array corresponding to the neutrino energy
1292 1292
1293 G4int G4NeutrinoNucleusModel::GetOnePionIndex 1293 G4int G4NeutrinoNucleusModel::GetOnePionIndex(G4double energy)
1294 { 1294 {
1295 G4int i, eIndex = 0; 1295 G4int i, eIndex = 0;
1296 1296
1297 for( i = 0; i < fOnePionIndex; i++) 1297 for( i = 0; i < fOnePionIndex; i++)
1298 { 1298 {
1299 if( energy <= fOnePionEnergy[i]*GeV ) 1299 if( energy <= fOnePionEnergy[i]*GeV )
1300 { 1300 {
1301 eIndex = i; 1301 eIndex = i;
1302 break; 1302 break;
1303 } 1303 }
1304 } 1304 }
1305 if( i >= fOnePionIndex ) eIndex = fOnePionI 1305 if( i >= fOnePionIndex ) eIndex = fOnePionIndex;
1306 // G4cout<<"eIndex = "<<eIndex<<G4endl; 1306 // G4cout<<"eIndex = "<<eIndex<<G4endl;
1307 return eIndex; 1307 return eIndex;
1308 } 1308 }
1309 1309
1310 ///////////////////////////////////////////// 1310 /////////////////////////////////////////////////////
1311 // 1311 //
1312 // nu_mu 1pi/Tot ratio for index-1, index lin 1312 // nu_mu 1pi/Tot ratio for index-1, index linear over energy
1313 1313
1314 G4double G4NeutrinoNucleusModel::GetNuMuOnePi 1314 G4double G4NeutrinoNucleusModel::GetNuMuOnePionProb(G4int index, G4double energy)
1315 { 1315 {
1316 G4double ratio(0.); 1316 G4double ratio(0.);
1317 1317
1318 if( index <= 0 || energy < fOnePionEn 1318 if( index <= 0 || energy < fOnePionEnergy[0] ) ratio = 0.;
1319 else if ( index >= fOnePionIndex ) rat 1319 else if ( index >= fOnePionIndex ) ratio = fOnePionProb[fOnePionIndex-1]*fOnePionEnergy[fOnePionIndex-1]*GeV/energy;
1320 else 1320 else
1321 { 1321 {
1322 G4double x1 = fOnePionEnergy[index-1]*GeV 1322 G4double x1 = fOnePionEnergy[index-1]*GeV;
1323 G4double x2 = fOnePionEnergy[index]*GeV; 1323 G4double x2 = fOnePionEnergy[index]*GeV;
1324 G4double y1 = fOnePionProb[index-1]; 1324 G4double y1 = fOnePionProb[index-1];
1325 G4double y2 = fOnePionProb[index]; 1325 G4double y2 = fOnePionProb[index];
1326 1326
1327 if( x1 >= x2) return fOnePionProb[index]; 1327 if( x1 >= x2) return fOnePionProb[index];
1328 else 1328 else
1329 { 1329 {
1330 G4double angle = (y2-y1)/(x2-x1); 1330 G4double angle = (y2-y1)/(x2-x1);
1331 ratio = y1 + (energy-x1)*angle; 1331 ratio = y1 + (energy-x1)*angle;
1332 } 1332 }
1333 } 1333 }
1334 return ratio; 1334 return ratio;
1335 } 1335 }
1336 1336
1337 ///////////////////////////////////////////// 1337 ////////////////////////////////////////////////////////////////////////////////////////////////////
1338 1338
1339 const G4double G4NeutrinoNucleusModel::fOnePi 1339 const G4double G4NeutrinoNucleusModel::fOnePionEnergy[58] =
1340 { 1340 {
1341 1341
1342 0.275314, 0.293652, 0.31729, 0.33409, 0.351 1342 0.275314, 0.293652, 0.31729, 0.33409, 0.351746, 0.365629, 0.380041, 0.400165, 0.437941, 0.479237,
1343 0.504391, 0.537803, 0.588487, 0.627532, 0.6 1343 0.504391, 0.537803, 0.588487, 0.627532, 0.686839, 0.791905, 0.878332, 0.987405, 1.08162, 1.16971,
1344 1.2982, 1.40393, 1.49854, 1.64168, 1.7524, 1344 1.2982, 1.40393, 1.49854, 1.64168, 1.7524, 1.87058, 2.02273, 2.15894, 2.3654, 2.55792, 2.73017,
1345 3.03005, 3.40733, 3.88128, 4.53725, 5.16786 1345 3.03005, 3.40733, 3.88128, 4.53725, 5.16786, 5.73439, 6.53106, 7.43879, 8.36214, 9.39965, 10.296,
1346 11.5735, 13.1801, 15.2052, 17.5414, 19.7178 1346 11.5735, 13.1801, 15.2052, 17.5414, 19.7178, 22.7462, 25.9026, 29.4955, 33.5867, 39.2516, 46.4716,
1347 53.6065, 63.4668, 73.2147, 85.5593, 99.9854 1347 53.6065, 63.4668, 73.2147, 85.5593, 99.9854
1348 }; 1348 };
1349 1349
1350 1350
1351 ///////////////////////////////////////////// 1351 ////////////////////////////////////////////////////////////////////////////////////////////////////
1352 1352
1353 const G4double G4NeutrinoNucleusModel::fOnePi 1353 const G4double G4NeutrinoNucleusModel::fOnePionProb[58] =
1354 { 1354 {
1355 0.0019357, 0.0189361, 0.0378722, 0.0502758, 1355 0.0019357, 0.0189361, 0.0378722, 0.0502758, 0.0662559, 0.0754581, 0.0865008, 0.0987275, 0.124112,
1356 0.153787, 0.18308, 0.213996, 0.245358, 0.27 1356 0.153787, 0.18308, 0.213996, 0.245358, 0.274425, 0.301536, 0.326612, 0.338208, 0.337806, 0.335948,
1357 0.328092, 0.313557, 0.304965, 0.292169, 0.2 1357 0.328092, 0.313557, 0.304965, 0.292169, 0.28481, 0.269474, 0.254138, 0.247499, 0.236249, 0.221654,
1358 0.205492, 0.198781, 0.182216, 0.162251, 0.1 1358 0.205492, 0.198781, 0.182216, 0.162251, 0.142878, 0.128631, 0.116001, 0.108435, 0.0974843, 0.082092,
1359 0.0755204, 0.0703121, 0.0607066, 0.0554278, 1359 0.0755204, 0.0703121, 0.0607066, 0.0554278, 0.0480401, 0.0427023, 0.0377123, 0.0323248, 0.0298584,
1360 0.0244296, 0.0218526, 0.019121, 0.016477, 0 1360 0.0244296, 0.0218526, 0.019121, 0.016477, 0.0137309, 0.0137963, 0.0110371, 0.00834028, 0.00686127, 0.00538226
1361 }; 1361 };
1362 1362
1363 //////////////////// QEratio(Z,A,Enu) 1363 //////////////////// QEratio(Z,A,Enu)
1364 1364
1365 G4double G4NeutrinoNucleusModel::CalculateQEr 1365 G4double G4NeutrinoNucleusModel::CalculateQEratioA( G4int Z, G4int A, G4double energy, G4int nepdg)
1366 { 1366 {
1367 energy /= GeV; 1367 energy /= GeV;
1368 G4double qerata(0.5), rr(0.), x1(0.), x2(0. 1368 G4double qerata(0.5), rr(0.), x1(0.), x2(0.), y1(0.), y2(0.), aa(0.);
1369 G4int i(0), N(0); 1369 G4int i(0), N(0);
1370 1370
1371 if( A > Z ) N = A-Z; 1371 if( A > Z ) N = A-Z;
1372 1372
1373 for( i = 0; i < 50; i++) 1373 for( i = 0; i < 50; i++)
1374 { 1374 {
1375 if( fQEnergy[i] >= energy ) break; 1375 if( fQEnergy[i] >= energy ) break;
1376 } 1376 }
1377 if(i <= 0) return 1.; 1377 if(i <= 0) return 1.;
1378 else if (i >= 49) return 0.; 1378 else if (i >= 49) return 0.;
1379 else 1379 else
1380 { 1380 {
1381 x1 = fQEnergy[i-1]; 1381 x1 = fQEnergy[i-1];
1382 x2 = fQEnergy[i]; 1382 x2 = fQEnergy[i];
1383 1383
1384 if( nepdg == 12 || nepdg == 14 ) 1384 if( nepdg == 12 || nepdg == 14 )
1385 { 1385 {
1386 if( x1 >= x2) return fNeMuQEratio[i]; 1386 if( x1 >= x2) return fNeMuQEratio[i];
1387 1387
1388 y1 = fNeMuQEratio[i-1]; 1388 y1 = fNeMuQEratio[i-1];
1389 y2 = fNeMuQEratio[i]; 1389 y2 = fNeMuQEratio[i];
1390 } 1390 }
1391 else 1391 else
1392 { 1392 {
1393 if( x1 >= x2) return fANeMuQEratio[i]; 1393 if( x1 >= x2) return fANeMuQEratio[i];
1394 1394
1395 y1 = fANeMuQEratio[i-1]; 1395 y1 = fANeMuQEratio[i-1];
1396 y2 = fANeMuQEratio[i]; 1396 y2 = fANeMuQEratio[i];
1397 } 1397 }
1398 aa = (y2-y1)/(x2-x1); 1398 aa = (y2-y1)/(x2-x1);
1399 rr = y1 + (energy-x1)*aa; 1399 rr = y1 + (energy-x1)*aa;
1400 1400
1401 if( nepdg == 12 || nepdg == 14 ) qerata 1401 if( nepdg == 12 || nepdg == 14 ) qerata = N*rr/( N*rr + A*( 1 - rr ) );
1402 else 1402 else qerata = Z*rr/( Z*rr + A*( 1 - rr ) );
1403 } 1403 }
1404 fQEratioA = qerata; 1404 fQEratioA = qerata;
1405 1405
1406 return qerata; 1406 return qerata;
1407 } 1407 }
1408 1408
1409 const G4double G4NeutrinoNucleusModel::fQEner 1409 const G4double G4NeutrinoNucleusModel::fQEnergy[50] =
1410 { 1410 {
1411 0.12, 0.1416, 0.167088, 0.197164, 0.232653, 1411 0.12, 0.1416, 0.167088, 0.197164, 0.232653, 0.274531, 0.323946, 0.382257, 0.451063, 0.532254,
1412 0.62806, 0.741111, 0.874511, 1.03192, 1.217 1412 0.62806, 0.741111, 0.874511, 1.03192, 1.21767, 1.43685, 1.69548, 2.00067, 2.36079, 2.78573,
1413 3.28716, 3.87885, 4.57705, 5.40092, 6.37308 1413 3.28716, 3.87885, 4.57705, 5.40092, 6.37308, 7.52024, 8.87388, 10.4712, 12.356, 14.5801,
1414 17.2045, 20.3013, 23.9555, 28.2675, 33.3557 1414 17.2045, 20.3013, 23.9555, 28.2675, 33.3557, 39.3597, 46.4444, 54.8044, 64.6692, 76.3097,
1415 90.0454, 106.254, 125.379, 147.947, 174.578 1415 90.0454, 106.254, 125.379, 147.947, 174.578, 206.002, 243.082, 286.837, 338.468, 399.392
1416 }; 1416 };
1417 1417
1418 const G4double G4NeutrinoNucleusModel::fANeMu 1418 const G4double G4NeutrinoNucleusModel::fANeMuQEratio[50] =
1419 { 1419 {
1420 1, 1, 1, 1, 1, 1, 1, 0.97506, 0.920938, 0.8 1420 1, 1, 1, 1, 1, 1, 1, 0.97506, 0.920938, 0.847671, 0.762973, 0.677684, 0.597685,
1421 0.52538, 0.461466, 0.405329, 0.356154, 0.31 1421 0.52538, 0.461466, 0.405329, 0.356154, 0.312944, 0.274984, 0.241341, 0.211654, 0.185322,
1422 0.161991, 0.141339, 0.123078, 0.106952, 0.0 1422 0.161991, 0.141339, 0.123078, 0.106952, 0.0927909, 0.0803262, 0.0693698, 0.0598207, 0.0514545,
1423 0.044193, 0.0378696, 0.0324138, 0.0276955, 1423 0.044193, 0.0378696, 0.0324138, 0.0276955, 0.0236343, 0.0201497, 0.0171592, 0.014602, 0.0124182,
1424 0.0105536, 0.00896322, 0.00761004, 0.006458 1424 0.0105536, 0.00896322, 0.00761004, 0.00645821, 0.00547859, 0.00464595, 0.00393928,
1425 0.00333961, 0.00283086, 0.00239927 1425 0.00333961, 0.00283086, 0.00239927
1426 }; 1426 };
1427 1427
1428 const G4double G4NeutrinoNucleusModel::fNeMuQ 1428 const G4double G4NeutrinoNucleusModel::fNeMuQEratio[50] =
1429 { 1429 {
1430 1, 1, 1, 1, 1, 1, 1, 0.977592, 0.926073, 0. 1430 1, 1, 1, 1, 1, 1, 1, 0.977592, 0.926073, 0.858783, 0.783874, 0.706868, 0.63113, 0.558681,
1431 0.490818, 0.428384, 0.371865, 0.321413, 0.2 1431 0.490818, 0.428384, 0.371865, 0.321413, 0.276892, 0.237959, 0.204139, 0.1749, 0.149706, 0.128047,
1432 0.109456, 0.093514, 0.0798548, 0.0681575, 0 1432 0.109456, 0.093514, 0.0798548, 0.0681575, 0.0581455, 0.0495804, 0.0422578, 0.036002, 0.0306614,
1433 0.0261061, 0.0222231, 0.0189152, 0.0160987, 1433 0.0261061, 0.0222231, 0.0189152, 0.0160987, 0.0137011, 0.0116604, 0.00992366, 0.00844558, 0.00718766,
1434 0.00611714, 0.00520618, 0.00443105, 0.00377 1434 0.00611714, 0.00520618, 0.00443105, 0.00377158, 0.00321062, 0.0027335, 0.00232774, 0.00198258
1435 }; 1435 };
1436 1436
1437 // 1437 //
1438 // 1438 //
1439 /////////////////////////// 1439 ///////////////////////////
1440 1440