Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4NuTauNucleusCcModel.cc
Go to the documentation of this file.
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: G4NuTauNucleusCcModel.cc 91806 2015-08-06 12:20:45Z gcosmo $
27//
28// Geant4 Header : G4NuTauNucleusCcModel
29//
30// Author : V.Grichine 12.2.19
31//
32
33#include <iostream>
34#include <fstream>
35#include <sstream>
36
37#include "G4NuTauNucleusCcModel.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"
74#include "G4ReactionProductVector.hh"
75
76#include "G4GeneratorPrecompoundInterface.hh"
77#include "G4PreCompoundModel.hh"
78#include "G4ExcitationHandler.hh"
79
80
81#include "G4TauMinus.hh"
82#include "G4TauPlus.hh"
83#include "G4Nucleus.hh"
84#include "G4LorentzVector.hh"
85
86using namespace std;
87using namespace CLHEP;
88
89#ifdef G4MULTITHREADED
90 G4Mutex G4NuTauNucleusCcModel::numuNucleusModel = G4MUTEX_INITIALIZER;
91#endif
92
93
94G4NuTauNucleusCcModel::G4NuTauNucleusCcModel(const G4String& name)
95 : G4NeutrinoNucleusModel(name)
96{
97 fData = fMaster = false;
98 fMtau = 1.77686*GeV;
99 theTauMinus = G4TauMinus::TauMinus();
100 theTauPlus = G4TauPlus::TauPlus();
101 InitialiseModel();
102}
103
104
105G4NuTauNucleusCcModel::~G4NuTauNucleusCcModel()
106{}
107
108
109void G4NuTauNucleusCcModel::ModelDescription(std::ostream& outFile) const
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}
117
118/////////////////////////////////////////////////////////
119//
120// Read data from G4PARTICLEXSDATA (locally PARTICLEXSDATA)
121
122void G4NuTauNucleusCcModel::InitialiseModel()
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}
217
218/////////////////////////////////////////////////////////
219
220G4bool G4NuTauNucleusCcModel::IsApplicable(const G4HadProjectile & aPart,
221 G4Nucleus & )
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}
236
237/////////////////////////////////////////// ClusterDecay ////////////////////////////////////////////////////////////
238//
239//
240
241G4HadFinalState* G4NuTauNucleusCcModel::ApplyYourself(
242 const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
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}
527
528
529/////////////////////////////////////////////////////////////////////
530////////////////////////////////////////////////////////////////////
531///////////////////////////////////////////////////////////////////
532
533/////////////////////////////////////////////////
534//
535// sample x, then Q2
536
537void G4NuTauNucleusCcModel::SampleLVkr(const G4HadProjectile & aTrack, G4Nucleus& targetNucleus)
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}
734
735//
736//
737///////////////////////////
G4FindDataDir
const char * G4FindDataDir(const char *)
Definition: G4FindDataDir.cc:93
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition: G4LorentzVector.hh:33
G4RandomDirection
G4ThreeVector G4RandomDirection()
Definition: G4RandomDirection.hh:58
G4MUTEX_INITIALIZER
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:85
G4MUTEXLOCK
#define G4MUTEXLOCK(mutex)
Definition: G4Threading.hh:251
G4MUTEXUNLOCK
#define G4MUTEXUNLOCK(mutex)
Definition: G4Threading.hh:254
G4Mutex
std::mutex G4Mutex
Definition: G4Threading.hh:81
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:36
G4double
double G4double
Definition: G4Types.hh:83
G4bool
bool G4bool
Definition: G4Types.hh:86
G4int
int G4int
Definition: G4Types.hh:85
Z
const G4int Z[17]
Definition: G4WaterStopping.cc:51
A
const G4double A[17]
Definition: G4WaterStopping.cc:53
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
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::fNuMuQarrayKR
static G4double fNuMuQarrayKR[50][51][51]
Definition: G4NeutrinoNucleusModel.hh:212
G4NeutrinoNucleusModel::fNuMuXarrayKR
static G4double fNuMuXarrayKR[50][51]
Definition: G4NeutrinoNucleusModel.hh:210
G4NeutrinoNucleusModel::GetOnePionIndex
G4int GetOnePionIndex(G4double energy)
Definition: G4NeutrinoNucleusModel.cc:1293
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
G4NeutrinoNucleusModel::GetNuMuOnePionProb
G4double GetNuMuOnePionProb(G4int index, G4double energy)
Definition: G4NeutrinoNucleusModel.cc:1314
G4NeutrinoNucleusModel::fNuMuXdistrKR
static G4double fNuMuXdistrKR[50][50]
Definition: G4NeutrinoNucleusModel.hh:211
G4NeutrinoNucleusModel::fNuMuQdistrKR
static G4double fNuMuQdistrKR[50][51][50]
Definition: G4NeutrinoNucleusModel.hh:213
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::ModelDescription
virtual void ModelDescription(std::ostream &) const
Definition: G4NuTauNucleusCcModel.cc:109
G4NuTauNucleusCcModel::ApplyYourself
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
Definition: G4NuTauNucleusCcModel.cc:241
G4NuTauNucleusCcModel::SampleLVkr
void SampleLVkr(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
Definition: G4NuTauNucleusCcModel.cc:537
G4NuTauNucleusCcModel::IsApplicable
virtual G4bool IsApplicable(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
Definition: G4NuTauNucleusCcModel.cc:220
G4NuTauNucleusCcModel::G4NuTauNucleusCcModel
G4NuTauNucleusCcModel(const G4String &name="NuTauNuclCcModel")
Definition: G4NuTauNucleusCcModel.cc:94
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:357
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:109
G4ParticleTable::FindParticle
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
Definition: G4ParticleTable.cc:586
G4ParticleTable::GetParticleTable
static G4ParticleTable * GetParticleTable()
Definition: G4ParticleTable.cc:87
G4TauMinus::TauMinus
static G4TauMinus * TauMinus()
Definition: G4TauMinus.cc:134
G4TauPlus::TauPlus
static G4TauPlus * TauPlus()
Definition: G4TauPlus.cc:133
CLHEP
Definition: DoubConv.h:17