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
G4NuMuNucleusCcModel.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: 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"
74#include "G4ReactionProductVector.hh"
75
76#include "G4GeneratorPrecompoundInterface.hh"
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
86using namespace std;
87using namespace CLHEP;
88
89#ifdef G4MULTITHREADED
90 G4Mutex G4NuMuNucleusCcModel::numuNucleusModel = G4MUTEX_INITIALIZER;
91#endif
92
93
94G4NuMuNucleusCcModel::G4NuMuNucleusCcModel(const G4String& name)
95 : G4NeutrinoNucleusModel(name)
96{
97 fData = fMaster = false;
98 InitialiseModel();
99}
100
101
102G4NuMuNucleusCcModel::~G4NuMuNucleusCcModel()
103{}
104
105
106void G4NuMuNucleusCcModel::ModelDescription(std::ostream& outFile) const
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}
114
115/////////////////////////////////////////////////////////
116//
117// Read data from G4PARTICLEXSDATA (locally PARTICLEXSDATA)
118
119void G4NuMuNucleusCcModel::InitialiseModel()
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}
214
215/////////////////////////////////////////////////////////
216
217G4bool G4NuMuNucleusCcModel::IsApplicable(const G4HadProjectile & aPart,
218 G4Nucleus & )
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}
233
234/////////////////////////////////////////// ClusterDecay ////////////////////////////////////////////////////////////
235//
236//
237
238G4HadFinalState* G4NuMuNucleusCcModel::ApplyYourself(
239 const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
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}
524
525
526/////////////////////////////////////////////////////////////////////
527////////////////////////////////////////////////////////////////////
528///////////////////////////////////////////////////////////////////
529
530/////////////////////////////////////////////////
531//
532// sample x, then Q2
533
534void G4NuMuNucleusCcModel::SampleLVkr(const G4HadProjectile & aTrack, G4Nucleus& targetNucleus)
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}
731
732//
733//
734///////////////////////////
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::theMuonMinus
G4ParticleDefinition * theMuonMinus
Definition: G4NeutrinoNucleusModel.hh:162
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
G4NuMuNucleusCcModel::SampleLVkr
void SampleLVkr(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
Definition: G4NuMuNucleusCcModel.cc:534
G4NuMuNucleusCcModel::G4NuMuNucleusCcModel
G4NuMuNucleusCcModel(const G4String &name="NuMuNuclCcModel")
Definition: G4NuMuNucleusCcModel.cc:94
G4NuMuNucleusCcModel::IsApplicable
virtual G4bool IsApplicable(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
Definition: G4NuMuNucleusCcModel.cc:217
G4NuMuNucleusCcModel::ApplyYourself
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
Definition: G4NuMuNucleusCcModel.cc:238
G4NuMuNucleusCcModel::ModelDescription
virtual void ModelDescription(std::ostream &) const
Definition: G4NuMuNucleusCcModel.cc:106
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
CLHEP
Definition: DoubConv.h:17