Geant4 11.3.0
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | List of all members
G4NuMuNucleusCcModel Class Reference

#include <G4NuMuNucleusCcModel.hh>

+ Inheritance diagram for G4NuMuNucleusCcModel:

Public Member Functions

 G4NuMuNucleusCcModel (const G4String &name="NuMuNuclCcModel")
 
virtual ~G4NuMuNucleusCcModel ()
 
virtual void InitialiseModel ()
 
virtual G4bool IsApplicable (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
 
virtual G4HadFinalState * ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
 
void SampleLVkr (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
 
G4double GetMinNuMuEnergy ()
 
G4double ThresholdEnergy (G4double mI, G4double mF, G4double mP)
 
virtual void ModelDescription (std::ostream &) const
 
- Public Member Functions inherited from G4NeutrinoNucleusModel
 G4NeutrinoNucleusModel (const G4String &name="neutrino-nucleus")
 
virtual ~G4NeutrinoNucleusModel ()
 
G4double SampleXkr (G4double energy)
 
G4double GetXkr (G4int iEnergy, G4double prob)
 
G4double SampleQkr (G4double energy, G4double xx)
 
G4double GetQkr (G4int iE, G4int jX, G4double prob)
 
void ClusterDecay (G4LorentzVector &lvX, G4int qX)
 
void MesonDecay (G4LorentzVector &lvX, G4int qX)
 
void FinalBarion (G4LorentzVector &lvB, G4int qB, G4int pdgB)
 
void RecoilDeexcitation (G4Fragment &fragment)
 
void FinalMeson (G4LorentzVector &lvM, G4int qM, G4int pdgM)
 
void CoherentPion (G4LorentzVector &lvP, G4int pdgP, G4Nucleus &targetNucleus)
 
void SetCutEnergy (G4double ec)
 
G4double GetCutEnergy ()
 
G4double GetNuEnergy ()
 
G4double GetQtransfer ()
 
G4double GetQ2 ()
 
G4double GetXsample ()
 
G4int GetPDGencoding ()
 
G4bool GetCascade ()
 
G4bool GetString ()
 
G4double GetCosTheta ()
 
G4double GetEmu ()
 
G4double GetEx ()
 
G4double GetMuMass ()
 
G4double GetW2 ()
 
G4double GetM1 ()
 
G4double GetMr ()
 
G4double GetTr ()
 
G4double GetDp ()
 
G4bool GetfBreak ()
 
G4bool GetfCascade ()
 
G4bool GetfString ()
 
G4LorentzVector GetLVl ()
 
G4LorentzVector GetLVh ()
 
G4LorentzVector GetLVt ()
 
G4LorentzVector GetLVcpi ()
 
G4double GetMinNuMuEnergy ()
 
G4double ThresholdEnergy (G4double mI, G4double mF, G4double mP)
 
G4double GetQEratioA ()
 
void SetQEratioA (G4double qea)
 
G4double FinalMomentum (G4double mI, G4double mF, G4double mP, G4LorentzVector lvX)
 
G4double FermiMomentum (G4Nucleus &targetNucleus)
 
G4double NucleonMomentum (G4Nucleus &targetNucleus)
 
G4double GetEx (G4int A, G4bool fP)
 
G4double GgSampleNM (G4Nucleus &nucl)
 
G4int GetEnergyIndex (G4double energy)
 
G4double GetNuMuQeTotRat (G4int index, G4double energy)
 
G4int GetOnePionIndex (G4double energy)
 
G4double GetNuMuOnePionProb (G4int index, G4double energy)
 
G4double CalculateQEratioA (G4int Z, G4int A, G4double energy, G4int nepdg)
 
- Public Member Functions inherited from G4HadronicInteraction
 G4HadronicInteraction (const G4String &modelName="HadronicModel")
 
virtual ~G4HadronicInteraction ()
 
virtual G4double SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
G4double GetMinEnergy () const
 
G4double GetMinEnergy (const G4Material *aMaterial, const G4Element *anElement) const
 
void SetMinEnergy (G4double anEnergy)
 
void SetMinEnergy (G4double anEnergy, const G4Element *anElement)
 
void SetMinEnergy (G4double anEnergy, const G4Material *aMaterial)
 
G4double GetMaxEnergy () const
 
G4double GetMaxEnergy (const G4Material *aMaterial, const G4Element *anElement) const
 
void SetMaxEnergy (const G4double anEnergy)
 
void SetMaxEnergy (G4double anEnergy, const G4Element *anElement)
 
void SetMaxEnergy (G4double anEnergy, const G4Material *aMaterial)
 
G4int GetVerboseLevel () const
 
void SetVerboseLevel (G4int value)
 
const G4String & GetModelName () const
 
void DeActivateFor (const G4Material *aMaterial)
 
void ActivateFor (const G4Material *aMaterial)
 
void DeActivateFor (const G4Element *anElement)
 
void ActivateFor (const G4Element *anElement)
 
G4bool IsBlocked (const G4Material *aMaterial) const
 
G4bool IsBlocked (const G4Element *anElement) const
 
void SetRecoilEnergyThreshold (G4double val)
 
G4double GetRecoilEnergyThreshold () const
 
virtual const std::pair< G4double, G4double > GetFatalEnergyCheckLevels () const
 
virtual std::pair< G4double, G4double > GetEnergyMomentumCheckLevels () const
 
void SetEnergyMomentumCheckLevels (G4double relativeLevel, G4double absoluteLevel)
 
virtual void BuildPhysicsTable (const G4ParticleDefinition &)
 
 G4HadronicInteraction (const G4HadronicInteraction &right)=delete
 
const G4HadronicInteraction & operator= (const G4HadronicInteraction &right)=delete
 
G4bool operator== (const G4HadronicInteraction &right) const =delete
 
G4bool operator!= (const G4HadronicInteraction &right) const =delete
 

Additional Inherited Members

- Protected Member Functions inherited from G4HadronicInteraction
void SetModelName (const G4String &nam)
 
G4bool IsBlocked () const
 
void Block ()
 
- Protected Attributes inherited from G4NeutrinoNucleusModel
G4ParticleDefinition * theMuonMinus
 
G4ParticleDefinition * theMuonPlus
 
G4double fSin2tW
 
G4double fCutEnergy
 
G4int fNbin
 
G4int fIndex
 
G4int fEindex
 
G4int fXindex
 
G4int fQindex
 
G4int fOnePionIndex
 
G4int fPDGencoding
 
G4bool fCascade
 
G4bool fString
 
G4bool fProton
 
G4bool f2p2h
 
G4bool fBreak
 
G4double fNuEnergy
 
G4double fQ2
 
G4double fQtransfer
 
G4double fXsample
 
G4double fM1
 
G4double fM2
 
G4double fMt
 
G4double fMu
 
G4double fW2
 
G4double fMpi
 
G4double fW2pi
 
G4double fMinNuEnergy
 
G4double fDp
 
G4double fTr
 
G4double fEmu
 
G4double fEmuPi
 
G4double fEx
 
G4double fMr
 
G4double fCosTheta
 
G4double fCosThetaPi
 
G4double fQEratioA
 
G4LorentzVector fLVh
 
G4LorentzVector fLVl
 
G4LorentzVector fLVt
 
G4LorentzVector fLVcpi
 
G4GeneratorPrecompoundInterface * fPrecoInterface
 
G4PreCompoundModel * fPreCompound
 
G4ExcitationHandler * fDeExcitation
 
G4Nucleus * fRecoil
 
G4int fSecID
 
- Protected Attributes inherited from G4HadronicInteraction
G4HadFinalState theParticleChange
 
G4int verboseLevel
 
G4double theMinEnergy
 
G4double theMaxEnergy
 
G4bool isBlocked
 
- Static Protected Attributes inherited from G4NeutrinoNucleusModel
static const G4int fResNumber = 6
 
static const G4double fResMass [6]
 
static const G4int fClustNumber = 4
 
static const G4double fMesMass [4] = {1260., 980., 770., 139.57}
 
static const G4int fMesPDG [4] = {20213, 9000211, 213, 211}
 
static const G4double fBarMass [4] = {1700., 1600., 1232., 939.57}
 
static const G4int fBarPDG [4] = {12224, 32224, 2224, 2212}
 
static const G4double fNuMuResQ [50][50]
 
static const G4double fNuMuEnergy [50]
 
static const G4double fNuMuQeTotRat [50]
 
static const G4double fOnePionEnergy [58]
 
static const G4double fOnePionProb [58]
 
static const G4double fNuMuEnergyLogVector [50]
 
static G4double fNuMuXarrayKR [50][51] = {{1.0}}
 
static G4double fNuMuXdistrKR [50][50] = {{1.0}}
 
static G4double fNuMuQarrayKR [50][51][51] = {{{1.0}}}
 
static G4double fNuMuQdistrKR [50][51][50] = {{{1.0}}}
 
static const G4double fQEnergy [50]
 
static const G4double fANeMuQEratio [50]
 
static const G4double fNeMuQEratio [50]
 

Detailed Description

Definition at line 55 of file G4NuMuNucleusCcModel.hh.

Constructor & Destructor Documentation

◆ G4NuMuNucleusCcModel()

G4NuMuNucleusCcModel::G4NuMuNucleusCcModel ( const G4String & name = "NuMuNuclCcModel")

Definition at line 94 of file G4NuMuNucleusCcModel.cc.

95 : G4NeutrinoNucleusModel(name)
96{
97 fData = fMaster = false;
98 InitialiseModel();
99}
G4NeutrinoNucleusModel::G4NeutrinoNucleusModel
G4NeutrinoNucleusModel(const G4String &name="neutrino-nucleus")
Definition G4NeutrinoNucleusModel.cc:109

◆ ~G4NuMuNucleusCcModel()

G4NuMuNucleusCcModel::~G4NuMuNucleusCcModel ( )
virtual

Definition at line 102 of file G4NuMuNucleusCcModel.cc.

103{}

Member Function Documentation

◆ ApplyYourself()

G4HadFinalState * G4NuMuNucleusCcModel::ApplyYourself ( const G4HadProjectile & aTrack,
G4Nucleus & targetNucleus )
virtual

Implements G4NeutrinoNucleusModel.

Definition at line 238 of file G4NuMuNucleusCcModel.cc.

240{
241 theParticleChange.Clear();
242 fProton = f2p2h = fBreak = false;
243 fCascade = fString = false;
244 fLVh = fLVl = fLVt = fLVcpi = G4LorentzVector(0.,0.,0.,0.);
245
246 const G4HadProjectile* aParticle = &aTrack;
247 G4double energy = aParticle->GetTotalEnergy();
248
249 G4String pName = aParticle->GetDefinition()->GetParticleName();
250
251 if( energy < fMinNuEnergy )
252 {
253 theParticleChange.SetEnergyChange(energy);
254 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
255 return &theParticleChange;
256 }
257
258 SampleLVkr( aTrack, targetNucleus);
259
260 if( fBreak == true || fEmu < fMu ) // ~5*10^-6
261 {
262 // G4cout<<"ni, ";
263 theParticleChange.SetEnergyChange(energy);
264 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
265 return &theParticleChange;
266 }
267
268 // LVs of initial state
269
270 G4LorentzVector lvp1 = aParticle->Get4Momentum();
271 G4LorentzVector lvt1( 0., 0., 0., fM1 );
272 G4double mPip = G4ParticleTable::GetParticleTable()->FindParticle(211)->GetPDGMass();
273
274 // 1-pi by fQtransfer && nu-energy
275 G4LorentzVector lvpip1( 0., 0., 0., mPip );
276 G4LorentzVector lvsum, lv2, lvX;
277 G4ThreeVector eP;
278 G4double cost(1.), sint(0.), phi(0.), muMom(0.), massX2(0.), massX(0.), massR(0.), eCut(0.);
279 G4DynamicParticle* aLept = nullptr; // lepton lv
280
281 G4int Z = targetNucleus.GetZ_asInt();
282 G4int A = targetNucleus.GetA_asInt();
283 G4double mTarg = targetNucleus.AtomicMass(A,Z);
284 G4int pdgP(0), qB(0);
285 // G4double mSum = G4ParticleTable::GetParticleTable()->FindParticle(2212)->GetPDGMass() + mPip;
286
287 G4int iPi = GetOnePionIndex(energy);
288 G4double p1pi = GetNuMuOnePionProb( iPi, energy);
289
290 if( p1pi > G4UniformRand() && fCosTheta > 0.9 ) // && fQtransfer < 0.95*GeV ) // mu- & coherent pion + nucleus
291 {
292 // lvsum = lvp1 + lvpip1;
293 lvsum = lvp1 + lvt1;
294 // cost = fCosThetaPi;
295 cost = fCosTheta;
296 sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
297 phi = G4UniformRand()*CLHEP::twopi;
298 eP = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost );
299
300 // muMom = sqrt(fEmuPi*fEmuPi-fMu*fMu);
301 muMom = sqrt(fEmu*fEmu-fMu*fMu);
302
303 eP *= muMom;
304
305 // lv2 = G4LorentzVector( eP, fEmuPi );
306 // lv2 = G4LorentzVector( eP, fEmu );
307 lv2 = fLVl;
308
309 // lvX = lvsum - lv2;
310 lvX = fLVh;
311 massX2 = lvX.m2();
312 massX = lvX.m();
313 massR = fLVt.m();
314
315 if ( massX2 <= 0. ) // vmg: very rarely ~ (1-4)e-6 due to big Q2/x, to be improved
316 {
317 fCascade = true;
318 theParticleChange.SetEnergyChange(energy);
319 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
320 return &theParticleChange;
321 }
322 fW2 = massX2;
323
324 if( pName == "nu_mu" ) aLept = new G4DynamicParticle( theMuonMinus, lv2 );
325 // else if( pName == "anti_nu_mu") aLept = new G4DynamicParticle( theMuonPlus, lv2 );
326 else
327 {
328 theParticleChange.SetEnergyChange(energy);
329 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
330 return &theParticleChange;
331 }
332 if( pName == "nu_mu" ) pdgP = 211;
333 // else pdgP = -211;
334 // eCut = fMpi + 0.5*(fMpi*fMpi-massX2)/mTarg; // massX -> fMpi
335
336 if( A > 1 )
337 {
338 eCut = (fMpi + mTarg)*(fMpi + mTarg) - (massX + massR)*(massX + massR);
339 eCut /= 2.*massR;
340 eCut += massX;
341 }
342 else eCut = fM1 + fMpi;
343
344 if ( lvX.e() > eCut ) // && sqrt( GetW2() ) < 1.4*GeV ) //
345 {
346 CoherentPion( lvX, pdgP, targetNucleus);
347 }
348 else
349 {
350 fCascade = true;
351 theParticleChange.SetEnergyChange(energy);
352 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
353 return &theParticleChange;
354 }
355 theParticleChange.AddSecondary( aLept, fSecID );
356
357 return &theParticleChange;
358 }
359 else // lepton part in lab
360 {
361 lvsum = lvp1 + lvt1;
362 cost = fCosTheta;
363 sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
364 phi = G4UniformRand()*CLHEP::twopi;
365 eP = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost );
366
367 muMom = sqrt(fEmu*fEmu-fMu*fMu);
368
369 eP *= muMom;
370
371 lv2 = G4LorentzVector( eP, fEmu );
372 lv2 = fLVl;
373 lvX = lvsum - lv2;
374 lvX = fLVh;
375 massX2 = lvX.m2();
376
377 if ( massX2 <= 0. ) // vmg: very rarely ~ (1-4)e-6 due to big Q2/x, to be improved
378 {
379 fCascade = true;
380 theParticleChange.SetEnergyChange(energy);
381 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
382 return &theParticleChange;
383 }
384 fW2 = massX2;
385
386 if( pName == "nu_mu" ) aLept = new G4DynamicParticle( theMuonMinus, lv2 );
387 // else if( pName == "anti_nu_mu") aLept = new G4DynamicParticle( theMuonPlus, lv2 );
388 else
389 {
390 theParticleChange.SetEnergyChange(energy);
391 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
392 return &theParticleChange;
393 }
394 theParticleChange.AddSecondary( aLept, fSecID );
395 }
396
397 // hadron part
398
399 fRecoil = nullptr;
400
401 if( A == 1 )
402 {
403 if( pName == "nu_mu" ) qB = 2;
404 // else qB = 0;
405
406 // if( G4UniformRand() > 0.1 ) // > 0.9999 ) // > 0.0001 ) //
407 {
408 ClusterDecay( lvX, qB );
409 }
410 return &theParticleChange;
411 }
412 /*
413 // else
414 {
415 if( pName == "nu_mu" ) pdgP = 211;
416 else pdgP = -211;
417
418
419 if ( fQtransfer < 0.95*GeV ) // < 0.35*GeV ) //
420 {
421 if( lvX.m() > mSum ) CoherentPion( lvX, pdgP, targetNucleus);
422 }
423 }
424 return &theParticleChange;
425 }
426 */
427 G4Nucleus recoil;
428 G4double rM(0.), ratio = G4double(Z)/G4double(A);
429
430 if( ratio > G4UniformRand() ) // proton is excited
431 {
432 fProton = true;
433 recoil = G4Nucleus(A-1,Z-1);
434 fRecoil = &recoil;
435 rM = recoil.AtomicMass(A-1,Z-1);
436
437 if( pName == "nu_mu" ) // (++) state -> p + pi+
438 {
439 fMt = G4ParticleTable::GetParticleTable()->FindParticle(2212)->GetPDGMass()
440 + G4ParticleTable::GetParticleTable()->FindParticle(211)->GetPDGMass();
441 }
442 else // (0) state -> p + pi-, n + pi0
443 {
444 // fMt = G4ParticleTable::GetParticleTable()->FindParticle(2212)->GetPDGMass()
445 // + G4ParticleTable::GetParticleTable()->FindParticle(-211)->GetPDGMass();
446 }
447 }
448 else // excited neutron
449 {
450 fProton = false;
451 recoil = G4Nucleus(A-1,Z);
452 fRecoil = &recoil;
453 rM = recoil.AtomicMass(A-1,Z);
454
455 if( pName == "nu_mu" ) // (+) state -> n + pi+
456 {
457 fMt = G4ParticleTable::GetParticleTable()->FindParticle(2112)->GetPDGMass()
458 + G4ParticleTable::GetParticleTable()->FindParticle(211)->GetPDGMass();
459 }
460 else // (-) state -> n + pi-, // n + pi0
461 {
462 // fMt = G4ParticleTable::GetParticleTable()->FindParticle(2112)->GetPDGMass()
463 // + G4ParticleTable::GetParticleTable()->FindParticle(-211)->GetPDGMass();
464 }
465 }
466 // G4int index = GetEnergyIndex(energy);
467 G4int nepdg = aParticle->GetDefinition()->GetPDGEncoding();
468
469 G4double qeTotRat; // = GetNuMuQeTotRat(index, energy);
470 qeTotRat = CalculateQEratioA( Z, A, energy, nepdg);
471
472 G4ThreeVector dX = (lvX.vect()).unit();
473 G4double eX = lvX.e(); // excited nucleon
474 G4double mX = sqrt(massX2);
475 // G4double pX = sqrt( eX*eX - mX*mX );
476 // G4double sumE = eX + rM;
477
478 if( qeTotRat > G4UniformRand() || mX <= fMt ) // || eX <= 1232.*MeV) // QE
479 {
480 fString = false;
481
482 if( fProton )
483 {
484 fPDGencoding = 2212;
485 fMr = proton_mass_c2;
486 recoil = G4Nucleus(A-1,Z-1);
487 fRecoil = &recoil;
488 rM = recoil.AtomicMass(A-1,Z-1);
489 }
490 else // if( pName == "anti_nu_mu" )
491 {
492 fPDGencoding = 2112;
493 fMr = G4ParticleTable::GetParticleTable()->
494 FindParticle(fPDGencoding)->GetPDGMass(); // 939.5654133*MeV;
495 recoil = G4Nucleus(A-1,Z);
496 fRecoil = &recoil;
497 rM = recoil.AtomicMass(A-1,Z);
498 }
499 // sumE = eX + rM;
500 G4double eTh = fMr + 0.5*(fMr*fMr - mX*mX)/rM;
501
502 if( eX <= eTh ) // vmg, very rarely out of kinematics
503 {
504 fString = true;
505 theParticleChange.SetEnergyChange(energy);
506 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
507 return &theParticleChange;
508 }
509 // FinalBarion( fLVh, 0, fPDGencoding ); // p(n)+deexcited recoil
510 FinalBarion( lvX, 0, fPDGencoding ); // p(n)+deexcited recoil
511 }
512 else // if ( eX < 9500000.*GeV ) // < 25.*GeV) // < 95.*GeV ) // < 2.5*GeV ) //cluster decay
513 {
514 if ( fProton && pName == "nu_mu" ) qB = 2;
515 // else if( fProton && pName == "anti_nu_mu" ) qB = 0;
516 else if( !fProton && pName == "nu_mu" ) qB = 1;
517 // else if( !fProton && pName == "anti_nu_mu" ) qB = -1;
518
519
520 ClusterDecay( lvX, qB );
521 }
522 return &theParticleChange;
523}
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition G4LorentzVector.hh:33
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition G4ThreeVector.hh:36
G4double
double G4double
Definition G4Types.hh:83
G4int
int G4int
Definition G4Types.hh:85
A
const G4double A[17]
Definition G4WaterStopping.cc:53
G4UniformRand
#define G4UniformRand()
Definition Randomize.hh:52
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
CLHEP::HepLorentzVector::vect
Hep3Vector vect() const
G4HadProjectile::GetDefinition
const G4ParticleDefinition * GetDefinition() const
Definition G4HadProjectile.hh:86
G4HadProjectile::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition G4HadProjectile.hh:91
G4HadProjectile::GetTotalEnergy
G4double GetTotalEnergy() const
Definition G4HadProjectile.hh:106
G4NeutrinoNucleusModel::CoherentPion
void CoherentPion(G4LorentzVector &lvP, G4int pdgP, G4Nucleus &targetNucleus)
Definition G4NeutrinoNucleusModel.cc:595
G4NeutrinoNucleusModel::GetOnePionIndex
G4int GetOnePionIndex(G4double energy)
Definition G4NeutrinoNucleusModel.cc:1293
G4NeutrinoNucleusModel::theMuonMinus
G4ParticleDefinition * theMuonMinus
Definition G4NeutrinoNucleusModel.hh:162
G4NeutrinoNucleusModel::GetNuMuOnePionProb
G4double GetNuMuOnePionProb(G4int index, G4double energy)
Definition G4NeutrinoNucleusModel.cc:1314
G4NeutrinoNucleusModel::CalculateQEratioA
G4double CalculateQEratioA(G4int Z, G4int A, G4double energy, G4int nepdg)
Definition G4NeutrinoNucleusModel.cc:1365
G4NeutrinoNucleusModel::ClusterDecay
void ClusterDecay(G4LorentzVector &lvX, G4int qX)
Definition G4NeutrinoNucleusModel.cc:720
G4NeutrinoNucleusModel::FinalBarion
void FinalBarion(G4LorentzVector &lvB, G4int qB, G4int pdgB)
Definition G4NeutrinoNucleusModel.cc:449
G4NuMuNucleusCcModel::SampleLVkr
void SampleLVkr(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
Definition G4NuMuNucleusCcModel.cc:534
G4Nucleus::GetA_asInt
G4int GetA_asInt() const
Definition G4Nucleus.hh:99
G4Nucleus::GetZ_asInt
G4int GetZ_asInt() const
Definition G4Nucleus.hh:105
G4Nucleus::AtomicMass
G4double AtomicMass(const G4double A, const G4double Z, const G4int numberOfLambdas=0) const
Definition G4Nucleus.cc:369
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition G4ParticleDefinition.hh:96
G4ParticleTable::FindParticle
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
Definition G4ParticleTable.cc:498
G4ParticleTable::GetParticleTable
static G4ParticleTable * GetParticleTable()
Definition G4ParticleTable.cc:82
G4INCL::KinematicsUtils::energy
G4double energy(const ThreeVector &p, const G4double m)
Definition G4INCLKinematicsUtils.cc:185

◆ GetMinNuMuEnergy()

G4double G4NuMuNucleusCcModel::GetMinNuMuEnergy ( )
inline

Definition at line 76 of file G4NuMuNucleusCcModel.hh.

76{ return fMu + 0.5*fMu*fMu/fM1 + 4.*CLHEP::MeV; }; // kinematics + accuracy for sqrts

◆ InitialiseModel()

void G4NuMuNucleusCcModel::InitialiseModel ( )
virtual

Reimplemented from G4HadronicInteraction.

Definition at line 119 of file G4NuMuNucleusCcModel.cc.

120{
121 G4String pName = "nu_mu";
122 // G4String pName = "anti_nu_mu";
123
124 G4int nSize(0), i(0), j(0), k(0);
125
126 if(!fData)
127 {
128#ifdef G4MULTITHREADED
129 G4MUTEXLOCK(&numuNucleusModel);
130 if(!fData)
131 {
132#endif
133 fMaster = true;
134#ifdef G4MULTITHREADED
135 }
136 G4MUTEXUNLOCK(&numuNucleusModel);
137#endif
138 }
139
140 if(fMaster)
141 {
142 const char* path = G4FindDataDir("G4PARTICLEXSDATA");
143 std::ostringstream ost1, ost2, ost3, ost4;
144 ost1 << path << "/" << "neutrino" << "/" << pName << "/xarraycckr";
145
146 std::ifstream filein1( ost1.str().c_str() );
147
148 // filein.open("$PARTICLEXSDATA/");
149
150 filein1>>nSize;
151
152 for( k = 0; k < fNbin; ++k )
153 {
154 for( i = 0; i <= fNbin; ++i )
155 {
156 filein1 >> fNuMuXarrayKR[k][i];
157 // G4cout<< fNuMuXarrayKR[k][i] << " ";
158 }
159 }
160 // G4cout<<G4endl<<G4endl;
161
162 ost2 << path << "/" << "neutrino" << "/" << pName << "/xdistrcckr";
163 std::ifstream filein2( ost2.str().c_str() );
164
165 filein2>>nSize;
166
167 for( k = 0; k < fNbin; ++k )
168 {
169 for( i = 0; i < fNbin; ++i )
170 {
171 filein2 >> fNuMuXdistrKR[k][i];
172 // G4cout<< fNuMuXdistrKR[k][i] << " ";
173 }
174 }
175 // G4cout<<G4endl<<G4endl;
176
177 ost3 << path << "/" << "neutrino" << "/" << pName << "/q2arraycckr";
178 std::ifstream filein3( ost3.str().c_str() );
179
180 filein3>>nSize;
181
182 for( k = 0; k < fNbin; ++k )
183 {
184 for( i = 0; i <= fNbin; ++i )
185 {
186 for( j = 0; j <= fNbin; ++j )
187 {
188 filein3 >> fNuMuQarrayKR[k][i][j];
189 // G4cout<< fNuMuQarrayKR[k][i][j] << " ";
190 }
191 }
192 }
193 // G4cout<<G4endl<<G4endl;
194
195 ost4 << path << "/" << "neutrino" << "/" << pName << "/q2distrcckr";
196 std::ifstream filein4( ost4.str().c_str() );
197
198 filein4>>nSize;
199
200 for( k = 0; k < fNbin; ++k )
201 {
202 for( i = 0; i <= fNbin; ++i )
203 {
204 for( j = 0; j < fNbin; ++j )
205 {
206 filein4 >> fNuMuQdistrKR[k][i][j];
207 // G4cout<< fNuMuQdistrKR[k][i][j] << " ";
208 }
209 }
210 }
211 fData = true;
212 }
213}
G4FindDataDir
const char * G4FindDataDir(const char *)
Definition G4FindDataDir.cc:94
G4MUTEXLOCK
#define G4MUTEXLOCK(mutex)
Definition G4Threading.hh:235
G4MUTEXUNLOCK
#define G4MUTEXUNLOCK(mutex)
Definition G4Threading.hh:238
G4NeutrinoNucleusModel::fNuMuQarrayKR
static G4double fNuMuQarrayKR[50][51][51]
Definition G4NeutrinoNucleusModel.hh:104
G4NeutrinoNucleusModel::fNuMuXarrayKR
static G4double fNuMuXarrayKR[50][51]
Definition G4NeutrinoNucleusModel.hh:102
G4NeutrinoNucleusModel::fNuMuXdistrKR
static G4double fNuMuXdistrKR[50][50]
Definition G4NeutrinoNucleusModel.hh:103
G4NeutrinoNucleusModel::fNuMuQdistrKR
static G4double fNuMuQdistrKR[50][51][50]
Definition G4NeutrinoNucleusModel.hh:105

Referenced by G4NuMuNucleusCcModel().

◆ IsApplicable()

G4bool G4NuMuNucleusCcModel::IsApplicable ( const G4HadProjectile & aTrack,
G4Nucleus & targetNucleus )
virtual

Reimplemented from G4NeutrinoNucleusModel.

Definition at line 217 of file G4NuMuNucleusCcModel.cc.

219{
220 G4bool result = false;
221 G4String pName = aPart.GetDefinition()->GetParticleName();
222 G4double energy = aPart.GetTotalEnergy();
223
224 if( pName == "nu_mu" // || pName == "anti_nu_mu" )
225 &&
226 energy > fMinNuEnergy )
227 {
228 result = true;
229 }
230
231 return result;
232}
G4bool
bool G4bool
Definition G4Types.hh:86

◆ ModelDescription()

void G4NuMuNucleusCcModel::ModelDescription ( std::ostream & outFile) const
virtual

Reimplemented from G4NeutrinoNucleusModel.

Definition at line 106 of file G4NuMuNucleusCcModel.cc.

107{
108
109 outFile << "G4NuMuNucleusCcModel is a neutrino-nucleus (charge current) scattering\n"
110 << "model which uses the standard model \n"
111 << "transfer parameterization. The model is fully relativistic\n";
112
113}

◆ SampleLVkr()

void G4NuMuNucleusCcModel::SampleLVkr ( const G4HadProjectile & aTrack,
G4Nucleus & targetNucleus )

Definition at line 534 of file G4NuMuNucleusCcModel.cc.

535{
536 fBreak = false;
537 G4int A = targetNucleus.GetA_asInt(), iTer(0), iTerMax(100);
538 G4int Z = targetNucleus.GetZ_asInt();
539 G4double e3(0.), pMu2(0.), pX2(0.), nMom(0.), rM(0.), hM(0.), tM = targetNucleus.AtomicMass(A,Z);
540 G4double Ex(0.), ei(0.), nm2(0.);
541 G4double cost(1.), sint(0.), phi(0.), muMom(0.);
542 G4ThreeVector eP, bst;
543 const G4HadProjectile* aParticle = &aTrack;
544 G4LorentzVector lvp1 = aParticle->Get4Momentum();
545
546 if( A == 1 ) // hydrogen, no Fermi motion ???
547 {
548 fNuEnergy = aParticle->GetTotalEnergy();
549 iTer = 0;
550
551 do
552 {
553 fXsample = SampleXkr(fNuEnergy);
554 fQtransfer = SampleQkr(fNuEnergy, fXsample);
555 fQ2 = fQtransfer*fQtransfer;
556
557 if( fXsample > 0. )
558 {
559 fW2 = fM1*fM1 - fQ2 + fQ2/fXsample; // sample excited hadron mass
560 fEmu = fNuEnergy - fQ2/2./fM1/fXsample;
561 }
562 else
563 {
564 fW2 = fM1*fM1;
565 fEmu = fNuEnergy;
566 }
567 e3 = fNuEnergy + fM1 - fEmu;
568
569 if( e3 < sqrt(fW2) ) G4cout<<"energyX = "<<e3/GeV<<", fW = "<<sqrt(fW2)/GeV<<G4endl;
570
571 pMu2 = fEmu*fEmu - fMu*fMu;
572
573 if(pMu2 < 0.) { fBreak = true; return; }
574
575 pX2 = e3*e3 - fW2;
576
577 fCosTheta = fNuEnergy*fNuEnergy + pMu2 - pX2;
578 fCosTheta /= 2.*fNuEnergy*sqrt(pMu2);
579 iTer++;
580 }
581 while( ( abs(fCosTheta) > 1. || fEmu < fMu ) && iTer < iTerMax );
582
583 if( iTer >= iTerMax ) { fBreak = true; return; }
584
585 if( abs(fCosTheta) > 1.) // vmg: due to big Q2/x values. To be improved ...
586 {
587 G4cout<<"H2: fCosTheta = "<<fCosTheta<<", fEmu = "<<fEmu<<G4endl;
588 // fCosTheta = -1. + 2.*G4UniformRand();
589 if(fCosTheta < -1.) fCosTheta = -1.;
590 if(fCosTheta > 1.) fCosTheta = 1.;
591 }
592 // LVs
593
594 G4LorentzVector lvt1 = G4LorentzVector( 0., 0., 0., fM1 );
595 G4LorentzVector lvsum = lvp1 + lvt1;
596
597 cost = fCosTheta;
598 sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
599 phi = G4UniformRand()*CLHEP::twopi;
600 eP = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost );
601 muMom = sqrt(fEmu*fEmu-fMu*fMu);
602 eP *= muMom;
603 fLVl = G4LorentzVector( eP, fEmu );
604
605 fLVh = lvsum - fLVl;
606 fLVt = G4LorentzVector( 0., 0., 0., 0. ); // no recoil
607 }
608 else // Fermi motion, Q2 in nucleon rest frame
609 {
610 G4Nucleus recoil1( A-1, Z );
611 rM = recoil1.AtomicMass(A-1,Z);
612 do
613 {
614 // nMom = NucleonMomentumBR( targetNucleus ); // BR
615 nMom = GgSampleNM( targetNucleus ); // Gg
616 Ex = GetEx(A-1, fProton);
617 ei = tM - sqrt( (rM + Ex)*(rM + Ex) + nMom*nMom );
618 // ei = 0.5*( tM - s2M - 2*eX );
619
620 nm2 = ei*ei - nMom*nMom;
621 iTer++;
622 }
623 while( nm2 < 0. && iTer < iTerMax );
624
625 if( iTer >= iTerMax ) { fBreak = true; return; }
626
627 G4ThreeVector nMomDir = nMom*G4RandomDirection();
628
629 if( !f2p2h || A < 3 ) // 1p1h
630 {
631 // hM = tM - rM;
632
633 fLVt = G4LorentzVector( -nMomDir, sqrt( (rM + Ex)*(rM + Ex) + nMom*nMom ) ); // rM ); //
634 fLVh = G4LorentzVector( nMomDir, ei ); // hM); //
635 }
636 else // 2p2h
637 {
638 G4Nucleus recoil(A-2,Z-1);
639 rM = recoil.AtomicMass(A-2,Z-1)+sqrt(nMom*nMom+fM1*fM1);
640 hM = tM - rM;
641
642 fLVt = G4LorentzVector( nMomDir, sqrt( rM*rM+nMom*nMom ) );
643 fLVh = G4LorentzVector(-nMomDir, sqrt( hM*hM+nMom*nMom ) );
644 }
645 // G4cout<<hM<<", ";
646 // bst = fLVh.boostVector();
647
648 // lvp1.boost(-bst); // -> nucleon rest system, where Q2 transfer is ???
649
650 fNuEnergy = lvp1.e();
651 // G4double mN = fLVh.m(); // better mN = fM1 !? vmg
652 iTer = 0;
653
654 do // no FM!?, 5.4.20 vmg
655 {
656 fXsample = SampleXkr(fNuEnergy);
657 fQtransfer = SampleQkr(fNuEnergy, fXsample);
658 fQ2 = fQtransfer*fQtransfer;
659
660 // G4double mR = mN + fM1*(A-1.)*std::exp(-2.0*fQtransfer/mN); // recoil mass in+el
661
662 if( fXsample > 0. )
663 {
664 fW2 = fM1*fM1 - fQ2 + fQ2/fXsample; // sample excited hadron mass
665
666 // fW2 = mN*mN - fQ2 + fQ2/fXsample; // sample excited hadron mass
667 // fEmu = fNuEnergy - fQ2/2./mR/fXsample; // fM1->mN
668
669 fEmu = fNuEnergy - fQ2/2./fM1/fXsample; // fM1->mN
670 }
671 else
672 {
673 // fW2 = mN*mN;
674
675 fW2 = fM1*fM1;
676 fEmu = fNuEnergy;
677 }
678 // if(fEmu < 0.) G4cout<<"fEmu = "<<fEmu<<" hM = "<<hM<<G4endl;
679 // e3 = fNuEnergy + mR - fEmu;
680
681 e3 = fNuEnergy + fM1 - fEmu;
682
683 // if( e3 < sqrt(fW2) ) G4cout<<"energyX = "<<e3/GeV<<", fW = "<<sqrt(fW2)/GeV<<G4endl;
684
685 pMu2 = fEmu*fEmu - fMu*fMu;
686 pX2 = e3*e3 - fW2;
687
688 if(pMu2 < 0.) { fBreak = true; return; }
689
690 fCosTheta = fNuEnergy*fNuEnergy + pMu2 - pX2;
691 fCosTheta /= 2.*fNuEnergy*sqrt(pMu2);
692 iTer++;
693 }
694 while( ( abs(fCosTheta) > 1. || fEmu < fMu ) && iTer < iTerMax );
695
696 if( iTer >= iTerMax ) { fBreak = true; return; }
697
698 if( abs(fCosTheta) > 1.) // vmg: due to big Q2/x values. To be improved ...
699 {
700 G4cout<<"FM: fCosTheta = "<<fCosTheta<<", fEmu = "<<fEmu<<G4endl;
701 // fCosTheta = -1. + 2.*G4UniformRand();
702 if( fCosTheta < -1.) fCosTheta = -1.;
703 if( fCosTheta > 1.) fCosTheta = 1.;
704 }
705 // LVs
706 // G4LorentzVector lvt1 = G4LorentzVector( 0., 0., 0., mN ); // fM1 );
707
708 G4LorentzVector lvt1 = G4LorentzVector( 0., 0., 0., fM1 ); // fM1 );
709 G4LorentzVector lvsum = lvp1 + lvt1;
710
711 cost = fCosTheta;
712 sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
713 phi = G4UniformRand()*CLHEP::twopi;
714 eP = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost );
715 muMom = sqrt(fEmu*fEmu-fMu*fMu);
716 eP *= muMom;
717 fLVl = G4LorentzVector( eP, fEmu );
718 fLVh = lvsum - fLVl;
719
720 // if( fLVh.e() < mN || fLVh.m2() < 0.) { fBreak = true; return; }
721
722 if( fLVh.e() < fM1 || fLVh.m2() < 0.) { fBreak = true; return; }
723
724 // back to lab system
725
726 // fLVl.boost(bst);
727 // fLVh.boost(bst);
728 }
729 //G4cout<<iTer<<", "<<fBreak<<"; ";
730}
G4RandomDirection
G4ThreeVector G4RandomDirection()
Definition G4RandomDirection.hh:58
G4endl
#define G4endl
Definition G4ios.hh:67
G4cout
G4GLOB_DLL std::ostream G4cout
G4NeutrinoNucleusModel::SampleXkr
G4double SampleXkr(G4double energy)
Definition G4NeutrinoNucleusModel.cc:225
G4NeutrinoNucleusModel::SampleQkr
G4double SampleQkr(G4double energy, G4double xx)
Definition G4NeutrinoNucleusModel.cc:306
G4NeutrinoNucleusModel::GgSampleNM
G4double GgSampleNM(G4Nucleus &nucl)
Definition G4NeutrinoNucleusModel.cc:1156

Referenced by ApplyYourself().

◆ ThresholdEnergy()

G4double G4NuMuNucleusCcModel::ThresholdEnergy ( G4double mI,
G4double mF,
G4double mP )
inline

Definition at line 78 of file G4NuMuNucleusCcModel.hh.

79 {
80 G4double w = std::sqrt(fW2);
81 return w + 0.5*( (mP+mF)*(mP+mF)-(w+mI)*(w+mI) )/mI;
82 };
The documentation for this class was generated from the following files: