ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4NuMuNucleusCcModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4NuMuNucleusCcModel.cc
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: G4NuMuNucleusCcModel.cc 91806 2015-08-06 12:20:45Z gcosmo $
27 //
28 // Geant4 Header : G4NuMuNucleusCcModel
29 //
30 // Author : V.Grichine 12.2.19
31 //
32 
33 #include <iostream>
34 #include <fstream>
35 #include <sstream>
36 
37 #include "G4NuMuNucleusCcModel.hh"
38 // #include "G4NuMuNuclCcDistrKR.hh"
39 
40 // #include "G4NuMuResQX.hh"
41 
42 #include "G4SystemOfUnits.hh"
43 #include "G4ParticleTable.hh"
44 #include "G4ParticleDefinition.hh"
45 #include "G4IonTable.hh"
46 #include "Randomize.hh"
47 #include "G4RandomDirection.hh"
48 // #include "G4Threading.hh"
49 
50 // #include "G4Integrator.hh"
51 #include "G4DataVector.hh"
52 #include "G4PhysicsTable.hh"
53 /*
54 #include "G4CascadeInterface.hh"
55 // #include "G4BinaryCascade.hh"
56 #include "G4TheoFSGenerator.hh"
57 #include "G4LundStringFragmentation.hh"
58 #include "G4ExcitedStringDecay.hh"
59 #include "G4FTFModel.hh"
60 // #include "G4BinaryCascade.hh"
61 #include "G4HadFinalState.hh"
62 #include "G4HadSecondary.hh"
63 #include "G4HadronicInteractionRegistry.hh"
64 // #include "G4INCLXXInterface.hh"
65 #include "G4QGSModel.hh"
66 #include "G4QGSMFragmentation.hh"
67 #include "G4QGSParticipants.hh"
68 */
69 #include "G4KineticTrack.hh"
70 #include "G4DecayKineticTracks.hh"
71 #include "G4KineticTrackVector.hh"
72 #include "G4Fragment.hh"
73 #include "G4NucleiProperties.hh"
75 
77 #include "G4PreCompoundModel.hh"
78 #include "G4ExcitationHandler.hh"
79 
80 
81 #include "G4MuonMinus.hh"
82 #include "G4MuonPlus.hh"
83 #include "G4Nucleus.hh"
84 #include "G4LorentzVector.hh"
85 
86 using namespace std;
87 using namespace CLHEP;
88 
90 
91 const G4double G4NuMuNucleusCcModel::fResMass[6] = // [fResNumber] =
92  {2190., 1920., 1700., 1600., 1440., 1232. };
93 
95 
96 const G4double G4NuMuNucleusCcModel::fMesMass[4] = {1260., 980., 770., 139.57};
97 const G4int G4NuMuNucleusCcModel::fMesPDG[4] = {20213, 9000211, 213, 211};
98 
99 // const G4double G4NuMuNucleusCcModel::fBarMass[4] = {1905., 1600., 1232., 939.57};
100 // const G4int G4NuMuNucleusCcModel::fBarPDG[4] = {2226, 32224, 2224, 2212};
101 
102 const G4double G4NuMuNucleusCcModel::fBarMass[4] = {1700., 1600., 1232., 939.57};
103 const G4int G4NuMuNucleusCcModel::fBarPDG[4] = {12224, 32224, 2224, 2212};
104 
106 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,
107 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,
108 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,
109 30988.8, 35765.7, 41279, 47642.2, 54986.3, 63462.4, 73245.2, 84536, 97567.2, 112607, 129966 };
110 
111 
114 G4double G4NuMuNucleusCcModel::fNuMuQarrayKR[50][51][51] = {{{1.0}}};
115 G4double G4NuMuNucleusCcModel::fNuMuQdistrKR[50][51][50] = {{{1.0}}};
116 
117 #ifdef G4MULTITHREADED
118  G4Mutex G4NuMuNucleusCcModel::numuNucleusModel = G4MUTEX_INITIALIZER;
119 #endif
120 
121 
123  : G4NeutrinoNucleusModel(name)
124 {
125  fData = fMaster = false;
126  InitialiseModel();
127 }
128 
129 
131 {}
132 
133 
134 void G4NuMuNucleusCcModel::ModelDescription(std::ostream& outFile) const
135 {
136 
137  outFile << "G4NuMuNucleusCcModel is a neutrino-nucleus (charge current) scattering\n"
138  << "model which uses the standard model \n"
139  << "transfer parameterization. The model is fully relativistic\n";
140 
141 }
142 
144 //
145 // Read data from G4PARTICLEXSDATA (locally PARTICLEXSDATA)
146 
148 {
149  G4String pName = "nu_mu";
150 
151  G4int nSize(0), i(0), j(0), k(0);
152 
153  if(!fData)
154  {
155 #ifdef G4MULTITHREADED
156  G4MUTEXLOCK(&numuNucleusModel);
157  if(!fData)
158  {
159 #endif
160  fMaster = true;
161 #ifdef G4MULTITHREADED
162  }
163  G4MUTEXUNLOCK(&numuNucleusModel);
164 #endif
165  }
166 
167  if(fMaster)
168  {
169  char* path = getenv("G4PARTICLEXSDATA");
170  std::ostringstream ost1, ost2, ost3, ost4;
171  ost1 << path << "/" << "neutrino" << "/" << pName << "/xarraycckr";
172 
173  std::ifstream filein1( ost1.str().c_str() );
174 
175  // filein.open("$PARTICLEXSDATA/");
176 
177  filein1>>nSize;
178 
179  for( k = 0; k < fNbin; ++k )
180  {
181  for( i = 0; i <= fNbin; ++i )
182  {
183  filein1 >> fNuMuXarrayKR[k][i];
184  // G4cout<< fNuMuXarrayKR[k][i] << " ";
185  }
186  }
187  // G4cout<<G4endl<<G4endl;
188 
189  ost2 << path << "/" << "neutrino" << "/" << pName << "/xdistrcckr";
190  std::ifstream filein2( ost2.str().c_str() );
191 
192  filein2>>nSize;
193 
194  for( k = 0; k < fNbin; ++k )
195  {
196  for( i = 0; i < fNbin; ++i )
197  {
198  filein2 >> fNuMuXdistrKR[k][i];
199  // G4cout<< fNuMuXdistrKR[k][i] << " ";
200  }
201  }
202  // G4cout<<G4endl<<G4endl;
203 
204  ost3 << path << "/" << "neutrino" << "/" << pName << "/q2arraycckr";
205  std::ifstream filein3( ost3.str().c_str() );
206 
207  filein3>>nSize;
208 
209  for( k = 0; k < fNbin; ++k )
210  {
211  for( i = 0; i <= fNbin; ++i )
212  {
213  for( j = 0; j <= fNbin; ++j )
214  {
215  filein3 >> fNuMuQarrayKR[k][i][j];
216  // G4cout<< fNuMuQarrayKR[k][i][j] << " ";
217  }
218  }
219  }
220  // G4cout<<G4endl<<G4endl;
221 
222  ost4 << path << "/" << "neutrino" << "/" << pName << "/q2distrcckr";
223  std::ifstream filein4( ost4.str().c_str() );
224 
225  filein4>>nSize;
226 
227  for( k = 0; k < fNbin; ++k )
228  {
229  for( i = 0; i <= fNbin; ++i )
230  {
231  for( j = 0; j < fNbin; ++j )
232  {
233  filein4 >> fNuMuQdistrKR[k][i][j];
234  // G4cout<< fNuMuQdistrKR[k][i][j] << " ";
235  }
236  }
237  }
238  fData = true;
239  }
240 }
241 
243 
245  G4Nucleus & targetNucleus)
246 {
247  G4bool result = false;
248  G4String pName = aPart.GetDefinition()->GetParticleName();
249  G4double energy = aPart.GetTotalEnergy();
250 
251  if( pName == "nu_mu" // || pName == "anti_nu_mu" )
252  &&
253  energy > fMinNuEnergy )
254  {
255  result = true;
256  }
257  G4int Z = targetNucleus.GetZ_asInt();
258  Z *= 1;
259 
260  return result;
261 }
262 
264 //
265 //
266 
268  const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
269 {
271  fProton = f2p2h = fBreak = false;
272  const G4HadProjectile* aParticle = &aTrack;
273  G4double energy = aParticle->GetTotalEnergy();
274 
275  G4String pName = aParticle->GetDefinition()->GetParticleName();
276 
277  if( energy < fMinNuEnergy )
278  {
281  return &theParticleChange;
282  }
283  SampleLVkr( aTrack, targetNucleus);
284 
285  if( fBreak == true || fEmu < fMu ) // ~5*10^-6
286  {
287  // G4cout<<"ni, ";
290  return &theParticleChange;
291  }
292 
293  // LVs of initial state
294 
295  G4LorentzVector lvp1 = aParticle->Get4Momentum();
296  G4LorentzVector lvt1( 0., 0., 0., fM1 );
298 
299  // 1-pi by fQtransfer && nu-energy
300  G4LorentzVector lvpip1( 0., 0., 0., mPip );
301  G4LorentzVector lvsum, lv2, lvX;
302  G4ThreeVector eP;
303  G4double cost(1.), sint(0.), phi(0.), muMom(0.), massX2(0.), massX(0.), massR(0.), eCut(0.);
304  G4DynamicParticle* aLept = nullptr; // lepton lv
305 
306  G4int Z = targetNucleus.GetZ_asInt();
307  G4int A = targetNucleus.GetA_asInt();
308  G4double mTarg = targetNucleus.AtomicMass(A,Z);
309  G4int pdgP(0), qB(0);
310  // G4double mSum = G4ParticleTable::GetParticleTable()->FindParticle(2212)->GetPDGMass() + mPip;
311 
312  G4int iPi = GetOnePionIndex(energy);
313  G4double p1pi = GetNuMuOnePionProb( iPi, energy);
314 
315  if( p1pi > G4UniformRand() ) // && fQtransfer < 0.95*GeV ) // mu- & coherent pion + nucleus
316  {
317  // lvsum = lvp1 + lvpip1;
318  lvsum = lvp1 + lvt1;
319  // cost = fCosThetaPi;
320  cost = fCosTheta;
321  sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
323  eP = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost );
324 
325  // muMom = sqrt(fEmuPi*fEmuPi-fMu*fMu);
326  muMom = sqrt(fEmu*fEmu-fMu*fMu);
327 
328  eP *= muMom;
329 
330  // lv2 = G4LorentzVector( eP, fEmuPi );
331  // lv2 = G4LorentzVector( eP, fEmu );
332  lv2 = fLVl;
333 
334  // lvX = lvsum - lv2;
335  lvX = fLVh;
336  massX2 = lvX.m2();
337  massX = lvX.m();
338  massR = fLVt.m();
339 
340  if ( massX2 <= 0. ) // vmg: very rarely ~ (1-4)e-6 due to big Q2/x, to be improved
341  {
344  return &theParticleChange;
345  }
346  fW2 = massX2;
347 
348  if( pName == "nu_mu" ) aLept = new G4DynamicParticle( theMuonMinus, lv2 );
349  else if( pName == "anti_nu_mu") aLept = new G4DynamicParticle( theMuonPlus, lv2 );
350  else
351  {
354  return &theParticleChange;
355  }
356  if( pName == "nu_mu" ) pdgP = 211;
357  else pdgP = -211;
358 
359  // eCut = fMpi + 0.5*(fMpi*fMpi-massX2)/mTarg; // massX -> fMpi
360 
361  eCut = (fMpi + mTarg)*(fMpi + mTarg) - (massX + massR)*(massX + massR);
362  eCut /= 2.*massR;
363  eCut += massX;
364 
365  if ( lvX.e() > eCut ) // && sqrt( GetW2() ) < 1.4*GeV ) //
366  {
367  CoherentPion( lvX, pdgP, targetNucleus);
368  }
369  else
370  {
373  return &theParticleChange;
374  }
376 
377  return &theParticleChange;
378  }
379  else // lepton part in lab
380  {
381  lvsum = lvp1 + lvt1;
382  cost = fCosTheta;
383  sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
385  eP = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost );
386 
387  muMom = sqrt(fEmu*fEmu-fMu*fMu);
388 
389  eP *= muMom;
390 
391  lv2 = G4LorentzVector( eP, fEmu );
392  lv2 = fLVl;
393  lvX = lvsum - lv2;
394  lvX = fLVh;
395  massX2 = lvX.m2();
396 
397  if ( massX2 <= 0. ) // vmg: very rarely ~ (1-4)e-6 due to big Q2/x, to be improved
398  {
401  return &theParticleChange;
402  }
403  fW2 = massX2;
404 
405  if( pName == "nu_mu" ) aLept = new G4DynamicParticle( theMuonMinus, lv2 );
406  else if( pName == "anti_nu_mu") aLept = new G4DynamicParticle( theMuonPlus, lv2 );
407  else
408  {
411  return &theParticleChange;
412  }
414  }
415 
416  // hadron part
417 
418  fRecoil = nullptr;
419  fCascade = false;
420  fString = false;
421 
422  if( A == 1 )
423  {
424  if( pName == "nu_mu" ) qB = 2;
425  else qB = 0;
426 
427  // if( G4UniformRand() > 0.1 ) // > 0.9999 ) // > 0.0001 ) //
428  {
429  ClusterDecay( lvX, qB );
430  }
431  return &theParticleChange;
432  }
433  /*
434  // else
435  {
436  if( pName == "nu_mu" ) pdgP = 211;
437  else pdgP = -211;
438 
439 
440  if ( fQtransfer < 0.95*GeV ) // < 0.35*GeV ) //
441  {
442  if( lvX.m() > mSum ) CoherentPion( lvX, pdgP, targetNucleus);
443  }
444  }
445  return &theParticleChange;
446  }
447  */
448  G4Nucleus recoil;
449  G4double rM(0.), ratio = G4double(Z)/G4double(A);
450 
451  if( ratio > G4UniformRand() ) // proton is excited
452  {
453  fProton = true;
454  recoil = G4Nucleus(A-1,Z-1);
455  fRecoil = &recoil;
456  rM = recoil.AtomicMass(A-1,Z-1);
457 
458  if( pName == "nu_mu" ) // (++) state -> p + pi+
459  {
462  }
463  else // (0) state -> p + pi-, n + pi0
464  {
467  }
468  }
469  else // excited neutron
470  {
471  fProton = false;
472  recoil = G4Nucleus(A-1,Z);
473  fRecoil = &recoil;
474  rM = recoil.AtomicMass(A-1,Z);
475 
476  if( pName == "nu_mu" ) // (+) state -> n + pi+
477  {
480  }
481  else // (-) state -> n + pi-, // n + pi0
482  {
485  }
486  }
487  G4int index = GetEnergyIndex(energy);
488  G4double qeTotRat = GetNuMuQeTotRat(index, energy);
489 
490  G4ThreeVector dX = (lvX.vect()).unit();
491  G4double eX = lvX.e(); // excited nucleon
492  G4double mX = sqrt(massX2);
493  // G4double pX = sqrt( eX*eX - mX*mX );
494  // G4double sumE = eX + rM;
495 
496  if( qeTotRat > G4UniformRand() || mX <= fMt ) // || eX <= 1232.*MeV) // QE
497  {
498  fString = false;
499 
500  if( pName == "nu_mu" )
501  {
502  fPDGencoding = 2212;
504  recoil = G4Nucleus(A-1,Z);
505  fRecoil = &recoil;
506  rM = recoil.AtomicMass(A-1,Z);
507  }
508  else // if( pName == "anti_nu_mu" )
509  {
510  fPDGencoding = 2112;
512  FindParticle(fPDGencoding)->GetPDGMass(); // 939.5654133*MeV;
513  recoil = G4Nucleus(A-1,Z-1);
514  fRecoil = &recoil;
515  rM = recoil.AtomicMass(A-1,Z-1);
516  }
517  // sumE = eX + rM;
518  G4double eTh = fMr + 0.5*(fMr*fMr - mX*mX)/rM;
519 
520  if( eX <= eTh ) // vmg, very rarely out of kinematics
521  {
524  return &theParticleChange;
525  }
526  FinalBarion( fLVh, 0, fPDGencoding ); // p(n)+deexcited recoil
527  }
528  else if ( eX < 95000.*GeV ) // < 25.*GeV) // < 95.*GeV ) // < 2.5*GeV ) //cluster decay
529  {
530  if ( fProton && pName == "nu_mu" ) qB = 2;
531  else if( fProton && pName == "anti_nu_mu" ) qB = 0;
532  else if( !fProton && pName == "nu_mu" ) qB = 1;
533  else if( !fProton && pName == "anti_nu_mu" ) qB = -1;
534 
535  // if( G4UniformRand() > 0.1 )
536  {
537  ClusterDecay( lvX, qB );
538  }
539  // else
540  {
541  if( pName == "nu_mu" ) pdgP = 211;
542  else pdgP = -211;
543 
544  if ( fQtransfer < 0.95*GeV ) // < 0.99*GeV ) //
545  {
546  // if( lvX.m() > mSum ) CoherentPion( lvX, pdgP, targetNucleus);
547  }
548  }
549  }
550  else // string
551  {
552  return &theParticleChange;
553  }
554  return &theParticleChange;
555 }
556 
557 
561 
563 //
564 // sample x, then Q2
565 
566 void G4NuMuNucleusCcModel::SampleLVkr(const G4HadProjectile & aTrack, G4Nucleus& targetNucleus)
567 {
568  fBreak = false;
569  G4int A = targetNucleus.GetA_asInt(), iTer(0), iTerMax(100);
570  G4int Z = targetNucleus.GetZ_asInt();
571  G4double e3(0.), pMu2(0.), pX2(0.), nMom(0.), rM(0.), hM(0.), tM = targetNucleus.AtomicMass(A,Z);
572  G4double Ex(0.), ei(0.), nm2(0.);
573  G4double cost(1.), sint(0.), phi(0.), muMom(0.);
574  G4ThreeVector eP, bst;
575  const G4HadProjectile* aParticle = &aTrack;
576  G4LorentzVector lvp1 = aParticle->Get4Momentum();
577 
578  if( A == 1 ) // hydrogen, no Fermi motion ???
579  {
580  fNuEnergy = aParticle->GetTotalEnergy();
581  iTer = 0;
582 
583  do
584  {
588 
589  if( fXsample > 0. )
590  {
591  fW2 = fM1*fM1 - fQ2 + fQ2/fXsample; // sample excited hadron mass
592  fEmu = fNuEnergy - fQ2/2./fM1/fXsample;
593  }
594  else
595  {
596  fW2 = fM1*fM1;
597  fEmu = fNuEnergy;
598  }
599  e3 = fNuEnergy + fM1 - fEmu;
600 
601  if( e3 < sqrt(fW2) ) G4cout<<"energyX = "<<e3/GeV<<", fW = "<<sqrt(fW2)/GeV<<G4endl;
602 
603  pMu2 = fEmu*fEmu - fMu*fMu;
604  pX2 = e3*e3 - fW2;
605 
606  fCosTheta = fNuEnergy*fNuEnergy + pMu2 - pX2;
607  fCosTheta /= 2.*fNuEnergy*sqrt(pMu2);
608  iTer++;
609  }
610  while( ( abs(fCosTheta) > 1. || fEmu < fMu ) && iTer < iTerMax );
611 
612  if( iTer >= iTerMax ) { fBreak = true; return; }
613 
614  if( abs(fCosTheta) > 1.) // vmg: due to big Q2/x values. To be improved ...
615  {
616  G4cout<<"H2: fCosTheta = "<<fCosTheta<<", fEmu = "<<fEmu<<G4endl;
617  // fCosTheta = -1. + 2.*G4UniformRand();
618  if(fCosTheta < -1.) fCosTheta = -1.;
619  if(fCosTheta > 1.) fCosTheta = 1.;
620  }
621  // LVs
622 
623  G4LorentzVector lvt1 = G4LorentzVector( 0., 0., 0., fM1 );
624  G4LorentzVector lvsum = lvp1 + lvt1;
625 
626  cost = fCosTheta;
627  sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
629  eP = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost );
630  muMom = sqrt(fEmu*fEmu-fMu*fMu);
631  eP *= muMom;
632  fLVl = G4LorentzVector( eP, fEmu );
633 
634  fLVh = lvsum - fLVl;
635  fLVt = G4LorentzVector( 0., 0., 0., 0. ); // no recoil
636  }
637  else // Fermi motion, Q2 in nucleon rest frame
638  {
639  G4Nucleus recoil1( A-1, Z );
640  rM = recoil1.AtomicMass(A-1,Z);
641  do
642  {
643  // nMom = NucleonMomentumBR( targetNucleus ); // BR
644  nMom = GgSampleNM( targetNucleus ); // Gg
645  Ex = GetEx(A-1, fProton);
646  ei = tM - sqrt( (rM + Ex)*(rM + Ex) + nMom*nMom );
647  // ei = 0.5*( tM - s2M - 2*eX );
648 
649  nm2 = ei*ei - nMom*nMom;
650  iTer++;
651  }
652  while( nm2 < 0. && iTer < iTerMax );
653 
654  if( iTer >= iTerMax ) { fBreak = true; return; }
655 
656  G4ThreeVector nMomDir = nMom*G4RandomDirection();
657 
658  if( !f2p2h ) // 1p1h
659  {
660  // hM = tM - rM;
661 
662  fLVt = G4LorentzVector( -nMomDir, sqrt( (rM + Ex)*(rM + Ex) + nMom*nMom ) ); // rM ); //
663  fLVh = G4LorentzVector( nMomDir, ei ); // hM); //
664  }
665  else // 2p2h
666  {
667  G4Nucleus recoil(A-2,Z-1);
668  rM = recoil.AtomicMass(A-2,Z-1)+sqrt(nMom*nMom+fM1*fM1);
669  hM = tM - rM;
670 
671  fLVt = G4LorentzVector( nMomDir, sqrt( rM*rM+nMom*nMom ) );
672  fLVh = G4LorentzVector(-nMomDir, sqrt( hM*hM+nMom*nMom ) );
673  }
674  // G4cout<<hM<<", ";
675  bst = fLVh.boostVector();
676 
677  lvp1.boost(-bst); // -> nucleon rest system, where Q2 transfer is ???
678 
679  fNuEnergy = lvp1.e();
680  G4double mN = fLVh.m();
681  iTer = 0;
682 
683  do
684  {
688 
689  if( fXsample > 0. )
690  {
691  // fW2 = fM1*fM1 - fQ2 + fQ2/fXsample; // sample excited hadron mass
692  fW2 = mN*mN - fQ2 + fQ2/fXsample; // sample excited hadron mass
693  fEmu = fNuEnergy - fQ2/2./fM1/fXsample;
694  }
695  else
696  {
697  fW2 = fM1*fM1;
698  fEmu = fNuEnergy;
699  }
700 
701  // if(fEmu < 0.) G4cout<<"fEmu = "<<fEmu<<" hM = "<<hM<<G4endl;
702 
703  e3 = fNuEnergy + fM1 - fEmu;
704 
705  // if( e3 < sqrt(fW2) ) G4cout<<"energyX = "<<e3/GeV<<", fW = "<<sqrt(fW2)/GeV<<G4endl;
706 
707  pMu2 = fEmu*fEmu - fMu*fMu;
708  pX2 = e3*e3 - fW2;
709 
710  if(pMu2 < 0.)
711  {
712  fBreak = true;
713  return;
714  }
715  fCosTheta = fNuEnergy*fNuEnergy + pMu2 - pX2;
716  fCosTheta /= 2.*fNuEnergy*sqrt(pMu2);
717  iTer++;
718  }
719  while( ( abs(fCosTheta) > 1. || fEmu < fMu ) && iTer < iTerMax );
720 
721  if( iTer >= iTerMax ) { fBreak = true; return; }
722 
723  if( abs(fCosTheta) > 1.) // vmg: due to big Q2/x values. To be improved ...
724  {
725  G4cout<<"FM: fCosTheta = "<<fCosTheta<<", fEmu = "<<fEmu<<G4endl;
726  // fCosTheta = -1. + 2.*G4UniformRand();
727  if(fCosTheta < -1.) fCosTheta = -1.;
728  if(fCosTheta > 1.) fCosTheta = 1.;
729  }
730  // LVs
731  G4LorentzVector lvt1 = G4LorentzVector( 0., 0., 0., mN ); // fM1 );
732  G4LorentzVector lvsum = lvp1 + lvt1;
733 
734  cost = fCosTheta;
735  sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
737  eP = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost );
738  muMom = sqrt(fEmu*fEmu-fMu*fMu);
739  eP *= muMom;
740  fLVl = G4LorentzVector( eP, fEmu );
741  fLVh = lvsum - fLVl;
742  // back to lab system
743  fLVl.boost(bst);
744  fLVh.boost(bst);
745  }
746  //G4cout<<iTer<<", "<<fBreak<<"; ";
747 }
748 
750 
752 {
753  G4int i(0), nBin(50);
754  G4double xx(0.), prob = G4UniformRand();
755 
756  for( i = 0; i < nBin; ++i )
757  {
758  if( energy <= fNuMuEnergyLogVector[i] ) break;
759  }
760  if( i <= 0) // E-edge
761  {
762  fEindex = 0;
763  xx = GetXkr( 0, prob);
764  }
765  else if ( i >= nBin)
766  {
767  fEindex = nBin-1;
768  xx = GetXkr( nBin-1, prob);
769  }
770  else
771  {
772  fEindex = i;
773  G4double x1 = GetXkr(i-1,prob);
774  G4double x2 = GetXkr(i,prob);
775 
778  G4double e = G4Log(energy);
779 
780  if( e2 <= e1) xx = x1 + G4UniformRand()*(x2-x1);
781  else xx = x1 + (e-e1)*(x2-x1)/(e2-e1); // lin in energy log-scale
782  }
783  return xx;
784 }
785 
787 //
788 // sample X according to prob (xmin,1) at a given energy index iEnergy
789 
791 {
792  G4int i(0), nBin=50;
793  G4double xx(0.);
794 
795  for( i = 0; i < nBin; ++i )
796  {
797  if( prob <= fNuMuXdistrKR[iEnergy][i] )
798  break;
799  }
800  if(i <= 0 ) // X-edge
801  {
802  fXindex = 0;
803  xx = fNuMuXarrayKR[iEnergy][0];
804  }
805  if ( i >= nBin )
806  {
807  fXindex = nBin;
808  xx = fNuMuXarrayKR[iEnergy][nBin];
809  }
810  else
811  {
812  fXindex = i;
813  G4double x1 = fNuMuXarrayKR[iEnergy][i];
814  G4double x2 = fNuMuXarrayKR[iEnergy][i+1];
815 
816  G4double p1 = 0.;
817 
818  if( i > 0 ) p1 = fNuMuXdistrKR[iEnergy][i-1];
819 
820  G4double p2 = fNuMuXdistrKR[iEnergy][i];
821 
822  if( p2 <= p1 ) xx = x1 + G4UniformRand()*(x2-x1);
823  else xx = x1 + (prob-p1)*(x2-x1)/(p2-p1);
824  }
825  return xx;
826 }
827 
829 //
830 // Sample fQtransfer at a given Enu and fX
831 
833 {
834  G4int nBin(50), iE=fEindex, jX=fXindex;
835  G4double qq(0.), qq1(0.), qq2(0.);
836  G4double prob = G4UniformRand();
837 
838  // first E
839 
840  if( iE <= 0 )
841  {
842  qq1 = GetQkr( 0, jX, prob);
843  }
844  else if ( iE >= nBin-1)
845  {
846  qq1 = GetQkr( nBin-1, jX, prob);
847  }
848  else
849  {
850  G4double q1 = GetQkr(iE-1,jX, prob);
851  G4double q2 = GetQkr(iE,jX, prob);
852 
855  G4double e = G4Log(energy);
856 
857  if( e2 <= e1) qq1 = q1 + G4UniformRand()*(q2-q1);
858  else qq1 = q1 + (e-e1)*(q2-q1)/(e2-e1); // lin in energy log-scale
859  }
860 
861  // then X
862 
863  if( jX <= 0 )
864  {
865  qq2 = GetQkr( iE, 0, prob);
866  }
867  else if ( jX >= nBin)
868  {
869  qq2 = GetQkr( iE, nBin, prob);
870  }
871  else
872  {
873  G4double q1 = GetQkr(iE,jX-1, prob);
874  G4double q2 = GetQkr(iE,jX, prob);
875 
876  G4double e1 = G4Log(fNuMuXarrayKR[iE][jX-1]);
877  G4double e2 = G4Log(fNuMuXarrayKR[iE][jX]);
878  G4double e = G4Log(xx);
879 
880  if( e2 <= e1) qq2 = q1 + G4UniformRand()*(q2-q1);
881  else qq2 = q1 + (e-e1)*(q2-q1)/(e2-e1); // lin in energy log-scale
882  }
883  qq = 0.5*(qq1+qq2);
884 
885  return qq;
886 }
887 
889 //
890 // sample Q according to prob (qmin,qmax) at a given energy index iE and X index jX
891 
893 {
894  G4int i(0), nBin=50;
895  G4double qq(0.);
896 
897  for( i = 0; i < nBin; ++i )
898  {
899  if( prob <= fNuMuQdistrKR[iE][jX][i] )
900  break;
901  }
902  if(i <= 0 ) // Q-edge
903  {
904  fQindex = 0;
905  qq = fNuMuQarrayKR[iE][jX][0];
906  }
907  if ( i >= nBin )
908  {
909  fQindex = nBin;
910  qq = fNuMuQarrayKR[iE][jX][nBin];
911  }
912  else
913  {
914  fQindex = i;
915  G4double q1 = fNuMuQarrayKR[iE][jX][i];
916  G4double q2 = fNuMuQarrayKR[iE][jX][i+1];
917 
918  G4double p1 = 0.;
919 
920  if( i > 0 ) p1 = fNuMuQdistrKR[iE][jX][i-1];
921 
922  G4double p2 = fNuMuQdistrKR[iE][jX][i];
923 
924  if( p2 <= p1 ) qq = q1 + G4UniformRand()*(q2-q1);
925  else qq = q1 + (prob-p1)*(q2-q1)/(p2-p1);
926  }
927  return qq;
928 }
929 
930 //
931 //