| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // $Id: G4NeutrinoNucleusModel.cc 91806 2015-0
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
80
81 const G4double G4NeutrinoNucleusModel::fResMas
82 {2190., 1920., 1700., 1600., 1440., 1232. };
83
84 const G4int G4NeutrinoNucleusModel::fClustNumb
85
86 const G4double G4NeutrinoNucleusModel::fMesMas
87 const G4int G4NeutrinoNucleusModel::fMesPDG
88
89 // const G4double G4NeutrinoNucleusModel::fBar
90 // const G4int G4NeutrinoNucleusModel::fBar
91
92 const G4double G4NeutrinoNucleusModel::fBarMas
93 const G4int G4NeutrinoNucleusModel::fBarPDG
94
95 const G4double G4NeutrinoNucleusModel::fNuMuE
96 115.603, 133.424, 153.991, 177.729, 205.126, 2
97 745.387, 860.289, 992.903, 1145.96, 1322.61, 1
98 4806.1, 5546.97, 6402.04, 7388.91, 8527.92, 98
99 30988.8, 35765.7, 41279, 47642.2, 54986.3, 634
100
101
102 G4double G4NeutrinoNucleusModel::fNuMuXarrayKR
103 G4double G4NeutrinoNucleusModel::fNuMuXdistrKR
104 G4double G4NeutrinoNucleusModel::fNuMuQarrayKR
105 G4double G4NeutrinoNucleusModel::fNuMuQdistrKR
106
107 ///////////////////////////////////////////
108
109 G4NeutrinoNucleusModel::G4NeutrinoNucleusModel
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 = fBrea
121
122 fNuEnergy = fQ2 = fQtransfer = fXsample =
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-,
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
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 b
158
159 fPrecoInterface = new G4GeneratorPrecompou
160
161 if( !fPreCompound )
162 {
163 G4HadronicInteraction* p =
164 G4HadronicInteractionRegistry::Instanc
165 G4VPreCompoundModel* fPreCompound = stat
166
167 if(!fPreCompound)
168 {
169 fPreCompound = new G4PreCompoundModel();
170 }
171 fPrecoInterface->SetDeExcitation(fPreCom
172 }
173 fDeExcitation = GetDeExcitation()->GetExci
174 */
175
176 fDeExcitation = new G4ExcitationHandler();
177 fPreCompound = new G4PreCompoundModel(fDe
178 fPrecoInterface = new G4GeneratorPrecompoun
179 fPrecoInterface->SetDeExcitation(fPreCompoun
180
181 fPDGencoding = 0; // unphysical as default
182 fRecoil = nullptr;
183
184 // Creator model ID
185 fSecID = G4PhysicsModelCatalog::GetModelID(
186 }
187
188
189 G4NeutrinoNucleusModel::~G4NeutrinoNucleusMode
190 {
191 if(fPrecoInterface) delete fPrecoInterface;
192 }
193
194
195 void G4NeutrinoNucleusModel::ModelDescription(
196 {
197
198 outFile << "G4NeutrinoNucleusModel is a ne
199 << "model which uses the standard
200 << "transfer parameterization. Th
201
202 }
203
204 //////////////////////////////////////////////
205
206 G4bool G4NeutrinoNucleusModel::IsApplicable(co
207 G4Nucleus & )
208 {
209 G4bool result = false;
210 G4String pName = aPart.GetDefinition()->GetP
211 G4double energy = aPart.GetTotalEnergy();
212
213 if( pName == "nu_mu" // || pName == "anti_n
214 &&
215 energy > fMinNuEnergy
216 {
217 result = true;
218 }
219
220 return result;
221 }
222
223 //////////////////////////////////////
224
225 G4double G4NeutrinoNucleusModel::SampleXkr(G4d
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] ) br
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
251 G4double e2 = G4Log(fNuMuEnergyLogVector[i
252 G4double e = G4Log(energy);
253
254 if( e2 <= e1) xx = x1 + G4UniformRand()*(x
255 else xx = x1 + (e-e1)*(x2-x1)/(e2
256 }
257 return xx;
258 }
259
260 //////////////////////////////////////////////
261 //
262 // sample X according to prob (xmin,1) at a gi
263
264 G4double G4NeutrinoNucleusModel::GetXkr(G4int
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-
293
294 G4double p2 = fNuMuXdistrKR[iEnergy][i];
295
296 if( p2 <= p1 ) xx = x1 + G4UniformRand()*(
297 else xx = x1 + (prob-p1)*(x2-x1)
298 }
299 return xx;
300 }
301
302 //////////////////////////////////////
303 //
304 // Sample fQtransfer at a given Enu and fX
305
306 G4double G4NeutrinoNucleusModel::SampleQkr( G4
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[i
328 G4double e2 = G4Log(fNuMuEnergyLogVector[i
329 G4double e = G4Log(energy);
330
331 if( e2 <= e1) qq1 = q1 + G4UniformRand()*(
332 else qq1 = q1 + (e-e1)*(q2-q1)/(e
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()*(
355 else qq2 = q1 + (e-e1)*(q2-q1)/(e
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
365
366 G4double G4NeutrinoNucleusModel::GetQkr( G4int
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()*(
399 else qq = q1 + (prob-p1)*(q2-q1)
400 }
401 return qq;
402 }
403
404
405 //////////////////////////////////////////////
406 //
407 // Final meson to theParticleChange
408
409 void G4NeutrinoNucleusModel::FinalMeson( G4Lor
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)
416 {
417 G4ParticleDefinition* pd2 = G4ParticleTabl
418 G4DynamicParticle* dp2 = new G4DynamicP
419 theParticleChange.AddSecondary( dp2, fSecI
420 }
421 else // meson resonances
422 {
423 G4ParticleDefinition* rePart = G4ParticleT
424 FindParticle(pdg);
425 G4KineticTrack ddkt( rePart, 0., G4ThreeVe
426 G4KineticTrackVector* ddktv = ddkt.Decay()
427
428 G4DecayKineticTracks decay( ddktv );
429
430 for( unsigned int i = 0; i < ddktv->size()
431 {
432 G4DynamicParticle * aNew =
433 new G4DynamicParticle( ddktv->operator[]
434 ddktv->operator[]
435
436 // G4cout<<" "<<i<<", "<<aNew->Get
437
438 theParticleChange.AddSecondary( aNew, fS
439 delete ddktv->operator[](i);
440 }
441 delete ddktv;
442 }
443 }
444
445 //////////////////////////////////////////////
446 //
447 // Final barion to theParticleChange, and reco
448
449 void G4NeutrinoNucleusModel::FinalBarion( G4Lo
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()-
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
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); //->AtomicM
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 - p
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
511 {
512 G4ParticleDefinition* pd2 = G4ParticleTabl
513 G4DynamicParticle* dp2 = new G4DynamicP
514 theParticleChange.AddSecondary( dp2, fSecI
515
516 }
517 else // delta resonances
518 {
519 G4ParticleDefinition* rePart = G4ParticleT
520 G4KineticTrack ddkt( rePart, 0., G4ThreeVe
521 G4KineticTrackVector* ddktv = ddkt.Decay()
522
523 G4DecayKineticTracks decay( ddktv );
524
525 for( unsigned int i = 0; i < ddktv->size()
526 {
527 G4DynamicParticle * aNew =
528 new G4DynamicParticle( ddktv->operator[]
529 ddktv->operator[]
530
531 // G4cout<<" "<<i<<", "<<aNew->Get
532
533 theParticleChange.AddSecondary( aNew, fS
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; // (vRe
548
549
550 if( fRecoil ) // proton*?
551 {
552 G4double grM = G4NucleiProperties::GetNucl
553 G4double exE = fLVt.m() - grM;
554 if( exE < 5.*MeV ) exE = 5.*MeV + G4Unifor
555
556 const G4LorentzVector in4v( G4ThreeVector(
557 G4Fragment fragment( A, Z, in4v); // lvTar
558 fragment.SetNumberOfHoles(1);
559 fragment.SetExcEnergyAndMomentum( exE, lvT
560
561 RecoilDeexcitation(fragment);
562 }
563 else // momentum?
564 {
565 theParticleChange.SetLocalEnergyDeposit( f
566 }
567 }
568
569
570 //////////////////////////////////////////////
571 //
572 // Get final particles from excited recoil nuc
573
574 void G4NeutrinoNucleusModel::RecoilDeexcitatio
575 {
576 G4ReactionProductVector* products = fPreComp
577
578 if( products != nullptr )
579 {
580 for( auto & prod : *products ) // prod is
581 {
582 theParticleChange.AddSecondary(new G4Dyn
583
584
585 delete prod;
586 }
587 delete products;
588 }
589 }
590
591 ///////////////////////////////////////////
592 //
593 // Fragmentation of lvX directly to pion and r
594
595 void G4NeutrinoNucleusModel::CoherentPion( G4L
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()->Fi
603
604 // mN = 1.*139.57 + G4UniformRand()*(938. -
605
606 G4ThreeVector vN = lvP.boostVector(), bst(0.
607 // G4double gN = lvP.e()/lvP.m();
608 // G4LorentzVector lvNu(vN*gN*mN, mN*gN);
609 G4LorentzVector lvNu(0.,0.,0., mN); // lvN
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); //->Atom
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
638 // G4cout<<fLVt<<" "<<bst<<G4endl;
639 }
640 lvP.boost(-bst); // 9-3-20
641 fMr = G4ParticleTable::GetParticleTable()->F
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 - p
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-
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
669
670 fLVcpi = lvN;
671
672 G4ParticleDefinition* pd2 = G4ParticleTable:
673 G4DynamicParticle* dp2 = new G4DynamicPar
674 theParticleChange.AddSecondary( dp2, fSecID
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::GetNucl
688 G4double exE = fLVt.m() - grM;
689
690 if( exE < 5.*MeV ) exE = 5.*MeV + G4Unifor
691
692 const G4LorentzVector in4v( G4ThreeVector(
693 G4Fragment fragment( A, Z, in4v); // lvTar
694 fragment.SetNumberOfHoles(1);
695 fragment.SetExcEnergyAndMomentum( exE, lvT
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 G4Dynamic
707 theParticleChange.AddSecondary(aSec, fSe
708 }
709 else theParticleChange.SetLocalEnergyDepos
710 }
711 return;
712 }
713
714
715 //////////////////////////////////////////////
716 //
717 // Excited barion decay to meson and barion,
718 // mass distributions and charge exchange are
719
720 void G4NeutrinoNucleusModel::ClusterDecay( G4L
721 {
722 G4bool finB = false;
723 G4int pdgB(0), i(0), qM(0), qB(0); // p
724 G4double mM(0.), mB(0.), eM(0.), eB(0.), pM(
725 G4double mm1(0.), mm22(0.), M1(0.), M2(0.),
726
727 mX = lvX.m();
728
729 G4double mN = G4ParticleTable::GetParticleT
730 G4double mPi = G4ParticleTable::GetParticleT
731
732 // G4double deltaM = 1.*MeV; // 30.*MeV; /
733 G4double deltaMr[4] = { 0.*MeV, 0.*MeV, 100
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
741 {
742 if( mX >= fBarMass[i] )
743 {
744 pdgB = fBarPDG[i];
745 // mB = G4ParticleTable::GetParticleTa
746 break;
747 }
748 }
749 if( i == fClustNumber || i == fClustNumber-1
750 {
751 if ( qX == 2 || qX == 0) { pdgB = 221
752
753 else if( qX == 1 || qX == -1) { pdgB = 211
754
755 return FinalBarion( lvX, qB, pdgB);
756 }
757 else if( mX < fBarMass[i] + deltaMr[i] || m
758 {
759 finB = true; // final barion -> out
760
761 if ( qX == 1 && pdgB != 2212) pdgB =
762 else if( qX == 0 && pdgB != 2212) pdgB =
763 else if( qX == 0 && pdgB == 2212) pdgB =
764
765 if( finB ) return FinalBarion( lvX, qX, pd
766 }
767 // no barion resonance, try 1->2 decay in CO
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 = 221
777 else if( qX == 1 || qX == -1) { pdgB = 211
778
779 return FinalBarion(lvX, qB, pdgB);
780 }
781 else // try decay -> meson(cluster) + barion
782 {
783 // G4double sigmaM = 50.*MeV; // 100.*MeV;
784 G4double rand = G4UniformRand();
785
786 // mM = mm1*mm22/( mm1 + rand*(mm22 - mm1)
787 // mM = mm1*mm22/sqrt( mm1*mm1 + rand*(mm2
788 // mM = -sigmaM*log( (1.- rand)*exp(-mm22/
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::GetParticleTa
798 break;
799 }
800 }
801 M1 = G4ParticleTable::GetParticleTable()->
802 M2 = mX - mM;
803
804 if( M2 <= M1 ) //
805 {
806 if ( qX == 2 || qX == 0) { pdgB = 2
807 else if( qX == 1 || qX == -1) { pdgB = 2
808
809 return FinalBarion(lvX, qB, pdgB);
810 }
811 mB = M1 + G4UniformRand()*(M2-M1);
812 // mB = -sigmaM*log( (1.- rand)*exp(-M2/si
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
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
858
859 void G4NeutrinoNucleusModel::MesonDecay( G4Lor
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(
864 G4double mm1(0.), mm22(0.), M1(0.), M2(0.),
865
866 mX = lvX.m();
867 Tkin = lvX.e() - mX;
868
869 // if( mX < 1120*MeV && mX > 1020*MeV ) // p
870 if( mX < 1080*MeV && mX > 990*MeV && Tkin <
871 {
872 return FinalMeson( lvX, qB, 333);
873 }
874 G4double mPi = G4ParticleTable::GetParticleT
875
876 G4double deltaMr[4] = { 0.*MeV, 0.*MeV, 100
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
884 {
885 if( mX >= fMesMass[i] )
886 {
887 pdgB = fMesPDG[i];
888 // mB = G4ParticleTable::GetParticleTa
889 break;
890 }
891 }
892 if( i == fClustNumber ) // || i == fClustNum
893 {
894 if ( qX == 1) { pdgB = 211; qB = 1;
895 else if( qX == 0 ) { pdgB = 111; qB = 0;
896 else if( qX == -1) { pdgB = -211; qB = -1;
897
898 return FinalMeson( lvX, qB, pdgB);
899 }
900 else if( mX < fMesMass[i] + deltaMr[i] ) //
901 {
902 finB = true; // final barion -> out
903 pdgB = fMesPDG[i];
904
905 // if ( qX == 1 && pdgB != 2212) pdg
906
907 if( qX == 0 ) pdgB = pdgB - 100;
908 else if( qX == -1 ) pdgB = -pdgB;
909
910 if( finB ) return FinalMeson( lvX, qX, pdg
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;}
922 else if( qX == 0 ) { pdgB = 111; qB = 0;}
923 else if( qX == -1) { pdgB = -211; qB = -1;
924
925 return FinalMeson(lvX, qB, pdgB);
926 }
927 else // try decay -> pion + meson(cluster)
928 {
929 // G4double sigmaM = 50.*MeV; // 100.*MeV;
930 G4double rand = G4UniformRand();
931
932 if ( qX == 1 ) { qM = 1; qB = 0;}
933 else if( qX == 0 ) { qM = -1; qB = 1;} /
934 else if( qX == -1 ) { qM = -1; qB = 0;}
935 /*
936 mM = mPi;
937 if(qM == 0) mM = G4ParticleTable::GetParti
938 pdgM = fMesPDG[fClustNumber-1];
939 */
940 // mm1*mm22/( mm1 + rand*(mm22 - mm1) );
941 // mM = mm1*mm22/sqrt( mm1*mm1 + rand*(mm2
942 // mM = -sigmaM*log( (1.- rand)*exp(-mm22/
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::GetParticleTa
953 break;
954 }
955 }
956 if( i == fClustNumber || i == fClustNumber
957 {
958 if ( qX == 1) { pdgB = 211; qB = 1;
959 else if( qX == 0 ) { pdgB = 111; qB = 0;
960 else if( qX == -1) { pdgB = -211; qB = -
961
962 return FinalMeson( lvX, qB, pdgB);
963 }
964 else if( mX < fMesMass[i] + deltaMr[i] ) /
965 {
966 finB = true; // final barion -> out
967 pdgB = fMesPDG[i];
968
969 // if ( qX == 1 && pdgB != 2212) pdg
970
971 if( qX == 0 ) pdgB = pdgB - 100;
972 else if( qX == -1 ) pdgB = -pdgB;
973
974 if( finB ) return FinalMeson( lvX, qX, p
975 }
976
977 M1 = G4ParticleTable::GetParticleTable()->
978 M2 = mX - mM;
979
980 if( M2 <= M1 ) //
981 {
982 if ( qX == 1) { pdgB = 211; qB = 1;
983 else if( qX == 0 ) { pdgB = 111; qB = 0;
984 else if( qX == -1) { pdgB = -211; qB = -
985
986 return FinalMeson(lvX, qB, pdgB);
987 }
988 mB = M1 + G4UniformRand()*(M2-M1);
989 // mB = -sigmaM*log( (1.- rand)*exp(-M2/si
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 lv
1025
1026 G4double G4NeutrinoNucleusModel::FinalMomentu
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::FermiMomentu
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.;
1071 else if ( Z == 1 && A == 2 ) { kF = 87.*
1072 else if ( Z == 2 && A == 3 ) { kF = 134.
1073 else if ( Z == 6 && A == 12 ) { kF = 221.
1074 else if ( Z == 14 && A == 28 ) { kF = 239.
1075 else if ( Z == 26 && A == 56 ) { kF = 257.
1076 else if ( Z == 82 && A == 208 ) { kF = 265.
1077 else
1078 {
1079 kF = kp*ZpA*( 1 - pow( G4double(A), -t1 )
1080 }
1081 return kF;
1082 }
1083
1084 /////////////////////////////////////////////
1085 //
1086 // sample nucleon momentum of Fermi motion fo
1087
1088 G4double G4NeutrinoNucleusModel::NucleonMomen
1089 {
1090 G4int A = targetNucleus.GetA_asInt();
1091 G4double kF = FermiMomentum( targetNucleus)
1092 G4double mom(0.), kCut = 0.5*GeV; // kCut =
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 ) // 1p1
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
1116 {
1117 G4double eX(10.*MeV), a1(0.), a2(0.), e1(0.
1118 G4int i(0);
1119 const G4int maxBin = 12;
1120
1121 G4double refA[maxBin] = { 2., 6., 12., 16.,
1122
1123 G4double pEx[maxBin] = { 0., 12.2, 10.1, 10
1124
1125 G4double nEx[maxBin] = { 0., 12.2, 10., 10.
1126
1127 G4DataVector dE(12,0.);
1128
1129 if(fP) for( i = 0; i < maxBin; ++i ) dE[i]
1130 else dE[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 mo
1155
1156 G4double G4NeutrinoNucleusModel::GgSampleNM(G
1157 {
1158 f2p2h = false;
1159 G4double /* distr(0.), tail(0.), */ shift(1
1160 G4double kF = FermiMomentum( nucl);
1161 G4double momMax = 2.*kF; // 1.*GeV; // 1.
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(G4doubl
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()*mo
1201 if( mom > 2.*kF ) f2p2h = true;
1202
1203 // mom = 0.;
1204
1205 return mom;
1206 }
1207
1208
1209 ///////////////////////////////////// experim
1210 //
1211 // Return index of nu/anu energy array corres
1212
1213 G4int G4NeutrinoNucleusModel::GetEnergyIndex(
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 line
1233
1234 G4double G4NeutrinoNucleusModel::GetNuMuQeTot
1235 {
1236 G4double ratio(0.);
1237 // GetMinNuMuEnergy()
1238 if( index <= 0 || energy < fNuMuEnergy[0] )
1239 else if (index >= fIndex) ratio = fNuMuQeTo
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::fNuMuE
1260 {
1261 0.112103, 0.117359, 0.123119, 0.129443, 0.1
1262 0.144084, 0.152576, 0.161991, 0.172458, 0.1
1263 0.197171, 0.211801, 0.228261, 0.24684, 0.26
1264 0.291816, 0.319125, 0.350417, 0.386422, 0.4
1265 0.47634, 0.532692, 0.598756, 0.676612, 0.76
1266 0.878812, 1.01062, 1.16963, 1.36271, 1.5987
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::fNuMuQ
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
1281 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0
1282 0.97, 0.96, 0.95, 0.93,
1283 0.917794, 0.850239, 0.780412, 0.709339, 0.6
1284 0.500236, 0.435528, 0.375015, 0.319157, 0.2
1285 0.148627, 0.119008, 0.0940699, 0.0733255, 0
1286 0.0235026, 0.0170486, 0.0122149, 0.00857825
1287 };
1288
1289 /////////////////////////////////////////////
1290 //
1291 // Return index of one pion array correspondi
1292
1293 G4int G4NeutrinoNucleusModel::GetOnePionIndex
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 = fOnePionI
1306 // G4cout<<"eIndex = "<<eIndex<<G4endl;
1307 return eIndex;
1308 }
1309
1310 /////////////////////////////////////////////
1311 //
1312 // nu_mu 1pi/Tot ratio for index-1, index lin
1313
1314 G4double G4NeutrinoNucleusModel::GetNuMuOnePi
1315 {
1316 G4double ratio(0.);
1317
1318 if( index <= 0 || energy < fOnePionEn
1319 else if ( index >= fOnePionIndex ) rat
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::fOnePi
1340 {
1341
1342 0.275314, 0.293652, 0.31729, 0.33409, 0.351
1343 0.504391, 0.537803, 0.588487, 0.627532, 0.6
1344 1.2982, 1.40393, 1.49854, 1.64168, 1.7524,
1345 3.03005, 3.40733, 3.88128, 4.53725, 5.16786
1346 11.5735, 13.1801, 15.2052, 17.5414, 19.7178
1347 53.6065, 63.4668, 73.2147, 85.5593, 99.9854
1348 };
1349
1350
1351 /////////////////////////////////////////////
1352
1353 const G4double G4NeutrinoNucleusModel::fOnePi
1354 {
1355 0.0019357, 0.0189361, 0.0378722, 0.0502758,
1356 0.153787, 0.18308, 0.213996, 0.245358, 0.27
1357 0.328092, 0.313557, 0.304965, 0.292169, 0.2
1358 0.205492, 0.198781, 0.182216, 0.162251, 0.1
1359 0.0755204, 0.0703121, 0.0607066, 0.0554278,
1360 0.0244296, 0.0218526, 0.019121, 0.016477, 0
1361 };
1362
1363 //////////////////// QEratio(Z,A,Enu)
1364
1365 G4double G4NeutrinoNucleusModel::CalculateQEr
1366 {
1367 energy /= GeV;
1368 G4double qerata(0.5), rr(0.), x1(0.), x2(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
1402 else
1403 }
1404 fQEratioA = qerata;
1405
1406 return qerata;
1407 }
1408
1409 const G4double G4NeutrinoNucleusModel::fQEner
1410 {
1411 0.12, 0.1416, 0.167088, 0.197164, 0.232653,
1412 0.62806, 0.741111, 0.874511, 1.03192, 1.217
1413 3.28716, 3.87885, 4.57705, 5.40092, 6.37308
1414 17.2045, 20.3013, 23.9555, 28.2675, 33.3557
1415 90.0454, 106.254, 125.379, 147.947, 174.578
1416 };
1417
1418 const G4double G4NeutrinoNucleusModel::fANeMu
1419 {
1420 1, 1, 1, 1, 1, 1, 1, 0.97506, 0.920938, 0.8
1421 0.52538, 0.461466, 0.405329, 0.356154, 0.31
1422 0.161991, 0.141339, 0.123078, 0.106952, 0.0
1423 0.044193, 0.0378696, 0.0324138, 0.0276955,
1424 0.0105536, 0.00896322, 0.00761004, 0.006458
1425 0.00333961, 0.00283086, 0.00239927
1426 };
1427
1428 const G4double G4NeutrinoNucleusModel::fNeMuQ
1429 {
1430 1, 1, 1, 1, 1, 1, 1, 0.977592, 0.926073, 0.
1431 0.490818, 0.428384, 0.371865, 0.321413, 0.2
1432 0.109456, 0.093514, 0.0798548, 0.0681575, 0
1433 0.0261061, 0.0222231, 0.0189152, 0.0160987,
1434 0.00611714, 0.00520618, 0.00443105, 0.00377
1435 };
1436
1437 //
1438 //
1439 ///////////////////////////
1440