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