Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
G4NuTauNucleusCcModel Class Reference

#include <G4NuTauNucleusCcModel.hh>

+ Inheritance diagram for G4NuTauNucleusCcModel:

Public Member Functions

 G4NuTauNucleusCcModel (const G4String &name="NuTauNuclCcModel")
 
virtual ~G4NuTauNucleusCcModel ()
 
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 G4NuTauNucleusCcModel.hh.

Constructor & Destructor Documentation

◆ G4NuTauNucleusCcModel()

G4NuTauNucleusCcModel::G4NuTauNucleusCcModel ( const G4String & name = "NuTauNuclCcModel")

Definition at line 94 of file G4NuTauNucleusCcModel.cc.

95 : G4NeutrinoNucleusModel(name)
96{
97 fData = fMaster = false;
98 fMtau = 1.77686*GeV;
99 theTauMinus = G4TauMinus::TauMinus();
100 theTauPlus = G4TauPlus::TauPlus();
101 InitialiseModel();
102}
G4NeutrinoNucleusModel::G4NeutrinoNucleusModel
G4NeutrinoNucleusModel(const G4String &name="neutrino-nucleus")
Definition G4NeutrinoNucleusModel.cc:109
G4TauMinus::TauMinus
static G4TauMinus * TauMinus()
Definition G4TauMinus.cc:133
G4TauPlus::TauPlus
static G4TauPlus * TauPlus()
Definition G4TauPlus.cc:133

◆ ~G4NuTauNucleusCcModel()

G4NuTauNucleusCcModel::~G4NuTauNucleusCcModel ( )
virtual

Definition at line 105 of file G4NuTauNucleusCcModel.cc.

106{}

Member Function Documentation

◆ ApplyYourself()

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

Implements G4NeutrinoNucleusModel.

Definition at line 241 of file G4NuTauNucleusCcModel.cc.

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

Definition at line 76 of file G4NuTauNucleusCcModel.hh.

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

◆ InitialiseModel()

void G4NuTauNucleusCcModel::InitialiseModel ( )
virtual

Reimplemented from G4HadronicInteraction.

Definition at line 122 of file G4NuTauNucleusCcModel.cc.

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

◆ IsApplicable()

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

Reimplemented from G4NeutrinoNucleusModel.

Definition at line 220 of file G4NuTauNucleusCcModel.cc.

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

◆ ModelDescription()

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

Reimplemented from G4NeutrinoNucleusModel.

Definition at line 109 of file G4NuTauNucleusCcModel.cc.

110{
111
112 outFile << "G4NuTauNucleusCcModel is a tau neutrino-nucleus (charge current) scattering\n"
113 << "model which uses the standard model \n"
114 << "transfer parameterization. The model is fully relativistic\n";
115
116}

◆ SampleLVkr()

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

Definition at line 537 of file G4NuTauNucleusCcModel.cc.

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

Definition at line 78 of file G4NuTauNucleusCcModel.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: