Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4QInelastic.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$
27 //
28 // ---------------- G4QInelastic class -----------------
29 // by Mikhail Kossov, December 2003.
30 // G4QInelastic class of the CHIPS Simulation Branch in GEANT4
31 // ---------------------------------------------------------------
32 // ****************************************************************************************
33 // This Header is a part of the CHIPS physics package (author: M. Kosov)
34 // ****************************************************************************************
35 // Short description: This is a universal class for the incoherent (inelastic)
36 // nuclear interactions in the CHIPS model.
37 // ---------------------------------------------------------------------------
38 //#define debug
39 //#define pdebug
40 //#define pickupd
41 //#define ldebug
42 //#define ppdebug
43 //#define qedebug
44 
45 #include "G4QInelastic.hh"
46 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh"
48 #include "G4HadronicDeprecate.hh"
49 
50 
51 // Initialization of static vectors
52 std::vector<G4int> G4QInelastic::ElementZ; // Z of the element(i) in theLastCalc
53 std::vector<G4double> G4QInelastic::ElProbInMat; // SumProbabilityElements in Material
54 std::vector<std::vector<G4int>*> G4QInelastic::ElIsoN;// N of isotope(j) of Element(i)
55 std::vector<std::vector<G4double>*>G4QInelastic::IsoProbInEl;//SumProbabIsotopes inElementI
56 
57 G4QInelastic::G4QInelastic(const G4String& processName):
58  G4VDiscreteProcess(processName, fHadronic)
59 {
60  G4HadronicDeprecate("G4QInelastic");
61 
62  EnMomConservation = G4LorentzVector(0.,0.,0.,0.);
63  nOfNeutrons = 0;
64 #ifdef debug
65  G4cout<<"G4QInelastic::Constructor is called"<<G4endl;
66 #endif
67  if (verboseLevel>0) G4cout << GetProcessName() << " process is created "<< G4endl;
68  //G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPSWorld (234 part.max)
69  G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio); // Clusterization param's
70  G4Quasmon::SetParameters(Temperature,SSin2Gluons,EtaEtaprime); // Hadronic parameters
71  G4QEnvironment::SetParameters(SolidAngle); // SolAngle of pbar-A secondary mesons capture
72  //@@ Initialize here the G4QuasmonString parameters
73 }
74 
75 G4bool G4QInelastic::manualFlag=false; // If false then standard parameters are used
76 G4double G4QInelastic::Temperature=140.; // Critical Temperature (sensitive at High En)
77 G4double G4QInelastic::SSin2Gluons=0.3; // Supression of s-quarks (in respect to u&d)
78 G4double G4QInelastic::EtaEtaprime=0.3; // Supression of eta mesons (gg->qq/3g->qq)
79 G4double G4QInelastic::freeNuc=.04; // Percentage of free nucleons on the surface
80 G4double G4QInelastic::freeDib=.08; // Percentage of free diBaryons on the surface
81 G4double G4QInelastic::clustProb=4.; // Nuclear clusterization parameter
82 G4double G4QInelastic::mediRatio=1.; // medium/vacuum hadronization ratio
83 //G4int G4QInelastic::nPartCWorld=152;// The#of particles initialized in CHIPS World
84 //G4int G4QInelastic::nPartCWorld=122;// The#of particles initialized in CHIPS World
85 G4int G4QInelastic::nPartCWorld=85; // The#of particles initialized in CHIPS World Red
86 G4double G4QInelastic::SolidAngle=0.5; // Part of Solid Angle to capture (@@A-dep.)
87 G4bool G4QInelastic::EnergyFlux=false; // Flag for Energy Flux use (not MultyQuasmon)
88 G4double G4QInelastic::PiPrThresh=141.4; // Pion Production Threshold for gammas
89 G4double G4QInelastic::M2ShiftVir=20000.;// Shift for M2=-Q2=m_pi^2 of the virtualGamma
90 G4double G4QInelastic::DiNuclMass=1870.; // DoubleNucleon Mass for VirtualNormalization
91 G4double G4QInelastic::photNucBias=1.; // BiasingParameter for photo(e,mu,tau)Nuclear
92 G4double G4QInelastic::weakNucBias=1.; // BiasingParameter for ChargedCurrents(nu,mu)
93 
94 void G4QInelastic::SetManual() {manualFlag=true;}
95 void G4QInelastic::SetStandard() {manualFlag=false;}
96 
97 // Fill the private parameters
98 void G4QInelastic::SetParameters(G4double temper, G4double ssin2g, G4double etaetap,
99  G4double fN, G4double fD, G4double cP, G4double mR,
100  G4int nParCW, G4double solAn, G4bool efFlag,
101  G4double piThresh, G4double mpisq, G4double dinum)
102 {
103  Temperature=temper;
104  SSin2Gluons=ssin2g;
105  EtaEtaprime=etaetap;
106  freeNuc=fN;
107  freeDib=fD;
108  clustProb=cP;
109  mediRatio=mR;
110  nPartCWorld = nParCW;
111  EnergyFlux=efFlag;
112  SolidAngle=solAn;
113  PiPrThresh=piThresh;
114  M2ShiftVir=mpisq;
115  DiNuclMass=dinum;
116  G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World with 234 particles
117  G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio); // Clusterization param's
118  G4Quasmon::SetParameters(Temperature,SSin2Gluons,EtaEtaprime); // Hadronic parameters
119  G4QEnvironment::SetParameters(SolidAngle); // SolAngle of pbar-A secondary mesons capture
120 }
121 
122 void G4QInelastic::SetPhotNucBias(G4double phnB) {photNucBias=phnB;}
123 void G4QInelastic::SetWeakNucBias(G4double ccnB) {weakNucBias=ccnB;}
124 
125 // Destructor
126 
127 G4QInelastic::~G4QInelastic()
128 {
129  // The following is just a copy of what is done in PostStepDoIt every interaction !
130  // The correction is if(IPIE), so just for(...;ip<IPIE;...) does not work ! @@
131  G4int IPIE=IsoProbInEl.size(); // How many old elements?
132  if(IPIE) for(G4int ip=0; ip<IPIE; ++ip) // Clean up the SumProb's of Isotopes (SPI)
133  {
134  std::vector<G4double>* SPI=IsoProbInEl[ip]; // Pointer to the SPI vector
135  SPI->clear();
136  delete SPI;
137  std::vector<G4int>* IsN=ElIsoN[ip]; // Pointer to the N vector
138  IsN->clear();
139  delete IsN;
140  }
141  ElProbInMat.clear(); // Clean up the SumProb's of Elements (SPE)
142  ElementZ.clear(); // Clear the body vector for Z of Elements
143  IsoProbInEl.clear(); // Clear the body vector for SPI
144  ElIsoN.clear(); // Clear the body vector for N of Isotopes
145 }
146 
147 
148 G4LorentzVector G4QInelastic::GetEnegryMomentumConservation()
149 {
150  return EnMomConservation;
151 }
152 
153 G4int G4QInelastic::GetNumberOfNeutronsInTarget()
154 {
155  return nOfNeutrons;
156 }
157 
158 // output of the function must be in units of length! L=1/sig_V,sig_V=SUM(n(j,i)*sig(j,i)),
159 // where n(i,j) is a number of nuclei of the isotop j of the element i in V=1(lengtUnit^3)
160 // ********** All CHIPS cross sections are calculated in the surface units ************
161 G4double G4QInelastic::GetMeanFreePath(const G4Track& aTrack,G4double,G4ForceCondition* Fc)
162 {
163 #ifdef debug
164  G4cout<<"G4QInelastic::GetMeanFreePath: Called Fc="<<*Fc<<G4endl;
165 #endif
166  *Fc = NotForced;
167 #ifdef debug
168  G4cout<<"G4QInelastic::GetMeanFreePath: Before GetDynPart"<<G4endl;
169 #endif
170  const G4DynamicParticle* incidentParticle = aTrack.GetDynamicParticle();
171 #ifdef debug
172  G4cout<<"G4QInelastic::GetMeanFreePath: Before GetDef"<<G4endl;
173 #endif
174  G4ParticleDefinition* incidentParticleDefinition=incidentParticle->GetDefinition();
175  if( !IsApplicable(*incidentParticleDefinition))
176  {
177  G4cout<<"-W-G4QInelastic::GetMeanFreePath called for not implemented particle"<<G4endl;
178  return DBL_MAX;
179  }
180  // Calculate the mean Cross Section for the set of Elements(*Isotopes) in the Material
181  G4double Momentum = incidentParticle->GetTotalMomentum(); // 3-momentum of the Particle
182 #ifdef debug
183  G4cout<<"G4QInelastic::GetMeanFreePath: BeforeGetMaterial"<<G4endl;
184 #endif
185  const G4Material* material = aTrack.GetMaterial(); // Get the current material
186  const G4double* NOfNucPerVolume = material->GetVecNbOfAtomsPerVolume();
187  const G4ElementVector* theElementVector = material->GetElementVector();
188  G4int nE=material->GetNumberOfElements();
189 #ifdef debug
190  G4cout<<"G4QInelastic::GetMeanFreePath:"<<nE<<" Elem's in theMaterial"<<G4endl;
191 #endif
192  //G4bool leptoNuc=false; // By default the reaction is not lepto-nuclear *Growing point*
193  G4VQCrossSection* CSmanager=0;
194  G4VQCrossSection* CSmanager2=0;
195  G4QNeutronCaptureRatio* capMan=0;
196  G4int pPDG=0;
197  G4int pZ = incidentParticleDefinition->GetAtomicNumber();
198  G4int pA = incidentParticleDefinition->GetAtomicMass();
199  if(incidentParticleDefinition == G4Neutron::Neutron())
200  {
201  CSmanager=G4QNeutronNuclearCrossSection::GetPointer();
202  capMan=G4QNeutronCaptureRatio::GetPointer();
203 #ifdef debug
204  G4cout<<"G4QInelastic::GetMeanFreePath: CSmanager is defined for the neutron"<<G4endl;
205 #endif
206  pPDG=2112;
207  }
208  else if(incidentParticleDefinition == G4Proton::Proton())
209  {
210  CSmanager=G4QProtonNuclearCrossSection::GetPointer();
211  pPDG=2212;
212  }
213  else if(incidentParticleDefinition == G4PionMinus::PionMinus())
214  {
215  CSmanager=G4QPionMinusNuclearCrossSection::GetPointer();
216  pPDG=-211;
217  }
218  else if(incidentParticleDefinition == G4PionPlus::PionPlus())
219  {
220  CSmanager=G4QPionPlusNuclearCrossSection::GetPointer();
221  pPDG=211;
222  }
223  else if(incidentParticleDefinition == G4KaonMinus::KaonMinus())
224  {
225  CSmanager=G4QKaonMinusNuclearCrossSection::GetPointer();
226  pPDG=-321;
227  }
228  else if(incidentParticleDefinition == G4KaonPlus::KaonPlus())
229  {
230  CSmanager=G4QKaonPlusNuclearCrossSection::GetPointer();
231  pPDG=321;
232  }
233  else if(incidentParticleDefinition == G4KaonZeroLong::KaonZeroLong() ||
234  incidentParticleDefinition == G4KaonZeroShort::KaonZeroShort() ||
235  incidentParticleDefinition == G4KaonZero::KaonZero() ||
236  incidentParticleDefinition == G4AntiKaonZero::AntiKaonZero() )
237  {
238  CSmanager=G4QKaonZeroNuclearCrossSection::GetPointer();
239  if(G4UniformRand() > 0.5) pPDG= 311;
240  else pPDG=-311;
241  }
242  else if(incidentParticleDefinition == G4Lambda::Lambda())
243  {
244  CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
245  pPDG=3122;
246  }
247  else if(pZ > 0 && pA > 1) // Ions (not implemented yet, should not be used)
248  {
249  G4cout<<"-Warning-G4QInelastic::GetMeanFreePath: G4QInelastic called for ions"<<G4endl;
250  CSmanager=G4QProtonNuclearCrossSection::GetPointer();
251  pPDG=90000000+999*pZ+pA;
252  }
253  else if(incidentParticleDefinition == G4SigmaPlus::SigmaPlus())
254  {
255  CSmanager=G4QHyperonPlusNuclearCrossSection::GetPointer();
256  pPDG=3222;
257  }
258  else if(incidentParticleDefinition == G4SigmaMinus::SigmaMinus())
259  {
260  CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
261  pPDG=3112;
262  }
263  else if(incidentParticleDefinition == G4SigmaZero::SigmaZero())
264  {
265  CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
266  pPDG=3212;
267  }
268  else if(incidentParticleDefinition == G4XiMinus::XiMinus())
269  {
270  CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
271  pPDG=3312;
272  }
273  else if(incidentParticleDefinition == G4XiZero::XiZero())
274  {
275  CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
276  pPDG=3322;
277  }
278  else if(incidentParticleDefinition == G4OmegaMinus::OmegaMinus())
279  {
280  CSmanager=G4QHyperonNuclearCrossSection::GetPointer();
281  pPDG=3334;
282  }
283  else if(incidentParticleDefinition == G4MuonPlus::MuonPlus() ||
284  incidentParticleDefinition == G4MuonMinus::MuonMinus())
285  {
286  CSmanager=G4QMuonNuclearCrossSection::GetPointer();
287  //leptoNuc=true;
288  pPDG=13;
289  }
290  else if(incidentParticleDefinition == G4AntiNeutron::AntiNeutron())
291  {
292  CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
293  pPDG=-2112;
294  }
295  else if(incidentParticleDefinition == G4AntiProton::AntiProton())
296  {
297  CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
298  pPDG=-2212;
299  }
300  else if(incidentParticleDefinition == G4AntiLambda::AntiLambda())
301  {
302  CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
303  pPDG=-3122;
304  }
305  else if(incidentParticleDefinition == G4AntiSigmaPlus::AntiSigmaPlus())
306  {
307  CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
308  pPDG=-3222;
309  }
310  else if(incidentParticleDefinition == G4AntiSigmaMinus::AntiSigmaMinus())
311  {
312  CSmanager=G4QAntiBaryonPlusNuclearCrossSection::GetPointer();
313  pPDG=-3112;
314  }
315  else if(incidentParticleDefinition == G4AntiSigmaZero::AntiSigmaZero())
316  {
317  CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
318  pPDG=-3212;
319  }
320  else if(incidentParticleDefinition == G4AntiXiMinus::AntiXiMinus())
321  {
322  CSmanager=G4QAntiBaryonPlusNuclearCrossSection::GetPointer();
323  pPDG=-3312;
324  }
325  else if(incidentParticleDefinition == G4AntiXiZero::AntiXiZero())
326  {
327  CSmanager=G4QAntiBaryonNuclearCrossSection::GetPointer();
328  pPDG=-3322;
329  }
330  else if(incidentParticleDefinition == G4AntiOmegaMinus::AntiOmegaMinus())
331  {
332  CSmanager=G4QAntiBaryonPlusNuclearCrossSection::GetPointer();
333  pPDG=-3334;
334  }
335  else if(incidentParticleDefinition == G4Gamma::Gamma())
336  {
337  CSmanager=G4QPhotonNuclearCrossSection::GetPointer();
338  pPDG=22;
339  }
340  else if(incidentParticleDefinition == G4Electron::Electron() ||
341  incidentParticleDefinition == G4Positron::Positron())
342  {
343  CSmanager=G4QElectronNuclearCrossSection::GetPointer();
344  //leptoNuc=true;
345  pPDG=11;
346  }
347  else if(incidentParticleDefinition == G4TauPlus::TauPlus() ||
348  incidentParticleDefinition == G4TauMinus::TauMinus())
349  {
350  CSmanager=G4QTauNuclearCrossSection::GetPointer();
351  //leptoNuc=true;
352  pPDG=15;
353  }
354  else if(incidentParticleDefinition == G4NeutrinoMu::NeutrinoMu() )
355  {
356  CSmanager=G4QNuMuNuclearCrossSection::GetPointer();
357  CSmanager2=G4QNuNuNuclearCrossSection::GetPointer();
358  //leptoNuc=true;
359  pPDG=14;
360  }
361  else if(incidentParticleDefinition == G4AntiNeutrinoMu::AntiNeutrinoMu() )
362  {
363  CSmanager=G4QANuMuNuclearCrossSection::GetPointer();
364  CSmanager2=G4QANuANuNuclearCrossSection::GetPointer();
365  //leptoNuc=true;
366  pPDG=-14;
367  }
368  else if(incidentParticleDefinition == G4NeutrinoE::NeutrinoE() )
369  {
370  CSmanager=G4QNuENuclearCrossSection::GetPointer();
371  CSmanager2=G4QNuNuNuclearCrossSection::GetPointer();
372  //leptoNuc=true;
373  pPDG=12;
374  }
375  else if(incidentParticleDefinition == G4AntiNeutrinoE::AntiNeutrinoE() )
376  {
377  CSmanager=G4QANuENuclearCrossSection::GetPointer();
378  CSmanager2=G4QANuANuNuclearCrossSection::GetPointer();
379  //leptoNuc=true;
380  pPDG=-12;
381  }
382  else
383  {
384  G4cout<<"-Warning-G4QInelastic::GetMeanFreePath:Particle "
385  <<incidentParticleDefinition->GetPDGEncoding()<<" isn't supported"<<G4endl;
386  return DBL_MAX;
387  }
388  G4QIsotope* Isotopes = G4QIsotope::Get(); // Pointer to the G4QIsotopes singleton
389  G4double sigma=0.; // Sums over elements for the material
390  G4int IPIE=IsoProbInEl.size(); // How many old elements?
391  if(IPIE) for(G4int ip=0; ip<IPIE; ++ip) // Clean up the SumProb's of Isotopes (SPI)
392  {
393  std::vector<G4double>* SPI=IsoProbInEl[ip]; // Pointer to the SPI vector
394  SPI->clear();
395  delete SPI;
396  std::vector<G4int>* IsN=ElIsoN[ip]; // Pointer to the N vector
397  IsN->clear();
398  delete IsN;
399  }
400  ElProbInMat.clear(); // Clean up the SumProb's of Elements (SPE)
401  ElementZ.clear(); // Clear the body vector for Z of Elements
402  IsoProbInEl.clear(); // Clear the body vector for SPI
403  ElIsoN.clear(); // Clear the body vector for N of Isotopes
404  for(G4int i=0; i<nE; ++i)
405  {
406  G4Element* pElement=(*theElementVector)[i]; // Pointer to the current element
407  G4int Z = static_cast<G4int>(pElement->GetZ()); // Z of the Element
408  ElementZ.push_back(Z); // Remember Z of the Element
409  G4int isoSize=0; // The default for the isoVectorLength is 0
410  G4int indEl=0; // Index of non-trivial element or 0(default)
411  G4IsotopeVector* isoVector=pElement->GetIsotopeVector(); // Get the predefined IsoVect
412  if(isoVector) isoSize=isoVector->size();// Get size of the existing isotopeVector
413 #ifdef debug
414  G4cout<<"G4QInelastic::GetMeanFreePath: isovectorLength="<<isoSize<<G4endl; // Result
415 #endif
416  if(isoSize) // The Element has non-trivial abundance set
417  {
418  indEl=pElement->GetIndex()+1; // Index of the non-trivial element
419  if(!Isotopes->IsDefined(Z,indEl)) // This index is not defined for this Z: define
420  {
421  std::vector<std::pair<G4int,G4double>*>* newAbund =
422  new std::vector<std::pair<G4int,G4double>*>;
423  G4double* abuVector=pElement->GetRelativeAbundanceVector();
424  for(G4int j=0; j<isoSize; j++) // Calculation of abundance vector for isotopes
425  {
426  G4int N=pElement->GetIsotope(j)->GetN()-Z; // N means A=N+Z !
427  if(pElement->GetIsotope(j)->GetZ()!=Z)G4cerr<<"G4QInelastic::GetMeanFreePath"
428  <<": Z="<<pElement->GetIsotope(j)->GetZ()<<"#"<<Z<<G4endl;
429  G4double abund=abuVector[j];
430  std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund);
431 #ifdef debug
432  G4cout<<"G4QInelastic::GetMeanFreePath: p#="<<j<<",N="<<N<<",ab="<<abund<<G4endl;
433 #endif
434  newAbund->push_back(pr);
435  }
436 #ifdef debug
437  G4cout<<"G4QInelastic::GetMeanFreePath: pairVectLength="<<newAbund->size()<<G4endl;
438 #endif
439  indEl=G4QIsotope::Get()->InitElement(Z,indEl,newAbund); // definition of the newInd
440  for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k]; // Cleaning temporary
441  delete newAbund; // Was "new" in the beginning of the name space
442  }
443  }
444  std::vector<std::pair<G4int,G4double>*>* cs= Isotopes->GetCSVector(Z,indEl);//CSPointer
445  std::vector<G4double>* SPI = new std::vector<G4double>; // Pointer to the SPI vector
446  IsoProbInEl.push_back(SPI);
447  std::vector<G4int>* IsN = new std::vector<G4int>; // Pointer to the N vector
448  ElIsoN.push_back(IsN);
449  G4int nIs=cs->size(); // A#Of Isotopes in the Element
450  G4double susi=0.; // sum of CS over isotopes
451 #ifdef debug
452  G4cout<<"G4QInelastic::GetMeanFreePath: Before Loop nIs="<<nIs<<G4endl;
453 #endif
454  if(nIs) for(G4int j=0; j<nIs; j++) // Calculate CS for eachIsotope of El
455  {
456  std::pair<G4int,G4double>* curIs=(*cs)[j]; // A pointer, which is used twice
457  G4int N=curIs->first; // #of Neuterons in the isotope j of El i
458  IsN->push_back(N); // Remember Min N for the Element
459 #ifdef debug
460  G4cout<<"G4QInel::GetMeanFrP: Before CS, P="<<Momentum<<",Z="<<Z<<",N="<<N<<G4endl;
461 #endif
462  if(!pPDG) G4cout<<"-Warning-G4QInelastic::GetMeanFrP: (1) projectile PDG=0"<<G4endl;
463  G4double CSI=CSmanager->GetCrossSection(true,Momentum,Z,N,pPDG);//CS(j,i) for isotope
464  if(CSmanager2)CSI+=CSmanager2->GetCrossSection(true,Momentum,Z,N,pPDG);//CS(j,i)nu,nu
465  if(capMan) CSI*=(1.-capMan->GetRatio(Momentum, Z, N));
466 #ifdef debug
467  G4cout<<"GQC::GMF:X="<<CSI<<",M="<<Momentum<<",Z="<<Z<<",N="<<N<<",P="<<pPDG<<G4endl;
468 #endif
469  curIs->second = CSI;
470  susi+=CSI; // Make a sum per isotopes
471  SPI->push_back(susi); // Remember summed cross-section
472  } // End of temporary initialization of the cross sections in the G4QIsotope singeltone
473  sigma+=Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i];//SUM(MeanCS*NOfNperV)
474  ElProbInMat.push_back(sigma);
475  } // End of LOOP over Elements
476 #ifdef debug
477  G4cout<<"G4QInel::GetMeanFrPa:S="<<sigma<<",e="<<photNucBias<<",w="<<weakNucBias<<G4endl;
478 #endif
479  // Check that cross section is not zero and return the mean free path
480  if(photNucBias!=1.) if(incidentParticleDefinition == G4Gamma::Gamma() ||
481  incidentParticleDefinition == G4MuonPlus::MuonPlus() ||
482  incidentParticleDefinition == G4MuonMinus::MuonMinus() ||
483  incidentParticleDefinition == G4Electron::Electron() ||
484  incidentParticleDefinition == G4Positron::Positron() ||
485  incidentParticleDefinition == G4TauMinus::TauMinus() ||
486  incidentParticleDefinition == G4TauPlus::TauPlus() )
487  sigma*=photNucBias;
488  if(weakNucBias!=1.) if(incidentParticleDefinition==G4NeutrinoE::NeutrinoE() ||
489  incidentParticleDefinition==G4AntiNeutrinoE::AntiNeutrinoE() ||
490  incidentParticleDefinition==G4NeutrinoTau::NeutrinoTau() ||
491  incidentParticleDefinition==G4AntiNeutrinoTau::AntiNeutrinoTau()||
492  incidentParticleDefinition==G4NeutrinoMu::NeutrinoMu() ||
493  incidentParticleDefinition==G4AntiNeutrinoMu::AntiNeutrinoMu() )
494  sigma*=weakNucBias;
495  if(sigma > 0.) return 1./sigma; // Mean path [distance]
496  return DBL_MAX;
497 }
498 
499 
500 G4bool G4QInelastic::IsApplicable(const G4ParticleDefinition& particle)
501 {
502  G4int Z=particle.GetAtomicNumber();
503  G4int A=particle.GetAtomicMass();
504  if (particle == *( G4Neutron::Neutron() )) return true;
505  else if (particle == *( G4Proton::Proton() )) return true;
506  else if (particle == *( G4MuonPlus::MuonPlus() )) return true;
507  else if (particle == *( G4MuonMinus::MuonMinus() )) return true;
508  else if (particle == *( G4Gamma::Gamma() )) return true;
509  else if (particle == *( G4Electron::Electron() )) return true;
510  else if (particle == *( G4Positron::Positron() )) return true;
511  else if (particle == *( G4PionMinus::PionMinus() )) return true;
512  else if (particle == *( G4PionPlus::PionPlus() )) return true;
513  else if (particle == *( G4KaonPlus::KaonPlus() )) return true;
514  else if (particle == *( G4KaonMinus::KaonMinus() )) return true;
515  else if (particle == *( G4KaonZeroLong::KaonZeroLong() )) return true;
516  else if (particle == *( G4KaonZeroShort::KaonZeroShort() )) return true;
517  else if (particle == *( G4Lambda::Lambda() )) return true;
518  else if (particle == *( G4SigmaPlus::SigmaPlus() )) return true;
519  else if (particle == *( G4SigmaMinus::SigmaMinus() )) return true;
520  else if (particle == *( G4SigmaZero::SigmaZero() )) return true;
521  else if (particle == *( G4XiMinus::XiMinus() )) return true;
522  else if (particle == *( G4XiZero::XiZero() )) return true;
523  else if (particle == *( G4OmegaMinus::OmegaMinus() )) return true;
524  else if (particle == *(G4GenericIon::GenericIon()) || (Z > 0 && A > 1)) return true;
525  else if (particle == *( G4AntiNeutron::AntiNeutron() )) return true;
526  else if (particle == *( G4AntiProton::AntiProton() )) return true;
527  else if (particle == *( G4AntiLambda::AntiLambda() )) return true;
528  else if (particle == *( G4AntiSigmaPlus::AntiSigmaPlus() )) return true;
529  else if (particle == *(G4AntiSigmaMinus::AntiSigmaMinus())) return true;
530  else if (particle == *( G4AntiSigmaZero::AntiSigmaZero() )) return true;
531  else if (particle == *( G4AntiXiMinus::AntiXiMinus() )) return true;
532  else if (particle == *( G4AntiXiZero::AntiXiZero() )) return true;
533  else if (particle == *(G4AntiOmegaMinus::AntiOmegaMinus())) return true;
534  else if (particle == *( G4TauPlus::TauPlus() )) return true;
535  else if (particle == *( G4TauMinus::TauMinus() )) return true;
536  else if (particle == *( G4AntiNeutrinoE::AntiNeutrinoE() )) return true;
537  else if (particle == *( G4NeutrinoE::NeutrinoE() )) return true;
538  else if (particle == *(G4AntiNeutrinoMu::AntiNeutrinoMu())) return true;
539  else if (particle == *( G4NeutrinoMu::NeutrinoMu() )) return true;
540  //else if (particle == *(G4AntiNeutrinoTau::AntiNeutrinoTau())) return true;
541  //else if (particle == *( G4NeutrinoTau::NeutrinoTau() )) return true;
542 #ifdef debug
543  G4cout<<"***G4QInelastic::IsApplicable: PDG="<<particle.GetPDGEncoding()<<G4endl;
544 #endif
545  return false;
546 }
547 
548 G4VParticleChange* G4QInelastic::PostStepDoIt(const G4Track& track, const G4Step& step)
549 {
550  static const G4double third = 1./3.;
551  static const G4double me=G4Electron::Electron()->GetPDGMass(); // electron mass
552  static const G4double me2=me*me; // squared electron mass
553  static const G4double mu=G4MuonMinus::MuonMinus()->GetPDGMass(); // muon mass
554  static const G4double mu2=mu*mu; // squared muon mass
555  static const G4double mt=G4TauMinus::TauMinus()->GetPDGMass(); // tau mass
556  static const G4double mt2=mt*mt; // squared tau mass
557  //static const G4double dpi=M_PI+M_PI; // 2*pi (for Phi distr.) ***changed to twopi***
558  static const G4double mNeut= G4QPDGCode(2112).GetMass();
559  static const G4double mNeut2= mNeut*mNeut;
560  static const G4double mProt= G4QPDGCode(2212).GetMass();
561  static const G4double mProt2= mProt*mProt;
562  static const G4double dM=mProt+mNeut; // doubled nucleon mass
563  static const G4double hdM=dM/2.; // M of the "nucleon"
564  static const G4double hdM2=hdM*hdM; // M2 of the "nucleon"
565  static const G4double mLamb= G4QPDGCode(3122).GetMass(); // Mass of Lambda/antiL
566  static const G4double mPi0 = G4QPDGCode(111).GetMass();
567  static const G4double mPi0s= mPi0*mPi0;
568  static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0);// Mass of deuteron
569  //static const G4double mDeut2=mDeut*mDeut; // SquaredMassOfDeuteron
570  static const G4double mTrit= G4QPDGCode(2112).GetNuclMass(1,2,0);// Mass of tritium
571  static const G4double mHel3= G4QPDGCode(2112).GetNuclMass(2,1,0);// Mass of Helium3
572  static const G4double mAlph= G4QPDGCode(2112).GetNuclMass(2,2,0);// Mass of alpha
573  static const G4double mPi = G4QPDGCode(211).GetMass();
574  static const G4double tmPi = mPi+mPi; // Doubled mass of the charged pion
575  static const G4double stmPi= tmPi*tmPi; // Squared Doubled mass of the charged pion
576  static const G4double mPPi = mPi+mProt; // Delta threshold
577  static const G4double mPPi2= mPPi*mPPi; // Delta low threshold for W2
578  //static const G4double mDel2= 1400*1400; // Delta up threshold for W2 (in MeV^2)
579  // Static definitions for electrons (nu,e) -----------------------------------------
580  static const G4double meN = mNeut+me;
581  static const G4double meN2= meN*meN;
582  static const G4double fmeN= 4*mNeut2*me2;
583  static const G4double mesN= mNeut2+me2;
584  static const G4double meP = mProt+me;
585  static const G4double meP2= meP*meP;
586  static const G4double fmeP= 4*mProt2*me2;
587  static const G4double mesP= mProt2+me2;
588  static const G4double medM= me2/dM; // for x limit
589  static const G4double meD = mPPi+me; // Multiperipheral threshold
590  static const G4double meD2= meD*meD;
591  // Static definitions for muons (nu,mu) -----------------------------------------
592  static const G4double muN = mNeut+mu;
593  static const G4double muN2= muN*muN;
594  static const G4double fmuN= 4*mNeut2*mu2;
595  static const G4double musN= mNeut2+mu2;
596  static const G4double muP = mProt+mu;
597  static const G4double muP2= muP*muP; // +
598  static const G4double fmuP= 4*mProt2*mu2; // +
599  static const G4double musP= mProt2+mu2;
600  static const G4double mudM= mu2/dM; // for x limit
601  static const G4double muD = mPPi+mu; // Multiperipheral threshold
602  static const G4double muD2= muD*muD;
603  // Static definitions for muons (nu,nu) -----------------------------------------
604  //static const G4double nuN = mNeut;
605  //static const G4double nuN2= mNeut2;
606  //static const G4double fnuN= 0.;
607  //static const G4double nusN= mNeut2;
608  //static const G4double nuP = mProt;
609  //static const G4double nuP2= mProt2;
610  //static const G4double fnuP= 0.;
611  //static const G4double nusP= mProt2;
612  //static const G4double nudM= 0.; // for x limit
613  //static const G4double nuD = mPPi; // Multiperipheral threshold
614  //static const G4double nuD2= mPPi2;
615  static const G4LorentzVector vacuum4M(0.,0.,0.,0.);
616  //-------------------------------------------------------------------------------------
617  static G4bool CWinit = true; // CHIPS Warld needs to be initted
618  if(CWinit)
619  {
620  CWinit=false;
621  G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part.max)
622  }
623  //-------------------------------------------------------------------------------------
624  const G4DynamicParticle* projHadron = track.GetDynamicParticle();
625  const G4ParticleDefinition* particle=projHadron->GetDefinition();
626 #ifdef debug
627  G4cout<<"G4QInelastic::PostStepDoIt: Before the GetMeanFreePath is called"<<G4endl;
628 #endif
629  G4ForceCondition cond=NotForced;
630  GetMeanFreePath(track, 1., &cond);
631 #ifdef debug
632  G4cout<<"G4QInelastic::PostStepDoIt: After the GetMeanFreePath is called"<<G4endl;
633 #endif
634  G4bool scat=false; // No CHEX in proj scattering
635  G4int scatPDG=0; // Must be filled if true (CHEX)
636  G4LorentzVector proj4M=projHadron->Get4Momentum(); // 4-momentum of the projectile (IU?)
637  G4LorentzVector scat4M=proj4M; // Must be filled if true
638  G4double momentum = projHadron->GetTotalMomentum(); // 3-momentum of the Particle
639  G4double Momentum=proj4M.rho();
640  if(std::fabs(Momentum-momentum)>.001)
641  G4cerr<<"*G4QInelastic::PostStepDoIt: P="<<Momentum<<"#"<<momentum<<G4endl;
642 #ifdef debug
643  G4double mp=proj4M.m();
644  G4cout<<"-->G4QInel::PostStDoIt:*called*, 4M="<<proj4M<<", P="<<Momentum<<"="<<momentum
645  <<",m="<<mp<<G4endl;
646 #endif
647  if (!IsApplicable(*particle)) // Check applicability
648  {
649  G4cerr<<"G4QInelastic::PostStepDoIt:Only gam,e+,e-,mu+,mu-,t+,t-,p are implemented."
650  <<G4endl;
651  return 0;
652  }
653  const G4Material* material = track.GetMaterial(); // Get the current material
654  G4int Z=0;
655  const G4ElementVector* theElementVector = material->GetElementVector();
656  G4int nE=material->GetNumberOfElements();
657 #ifdef debug
658  G4cout<<"G4QInelastic::PostStepDoIt: "<<nE<<" elements in the material."<<G4endl;
659 #endif
660  G4int projPDG=0; // PDG Code prototype for the captured hadron
661  G4int pZ=particle->GetAtomicNumber();
662  G4int pA=particle->GetAtomicMass();
663  if (particle == G4Neutron::Neutron() ) projPDG= 2112;
664  else if (particle == G4Proton::Proton() ) projPDG= 2212;
665  else if (particle == G4PionMinus::PionMinus() ) projPDG= -211;
666  else if (particle == G4PionPlus::PionPlus() ) projPDG= 211;
667  else if (particle == G4KaonPlus::KaonPlus() ) projPDG= 321;
668  else if (particle == G4KaonMinus::KaonMinus() ) projPDG= -321;
669  else if (particle == G4KaonZeroLong::KaonZeroLong() ||
670  particle == G4KaonZeroShort::KaonZeroShort() ||
671  particle == G4KaonZero::KaonZero() ||
672  particle == G4AntiKaonZero::AntiKaonZero() )
673  {
674  if(G4UniformRand() > 0.5) projPDG= 311;
675  else projPDG= -311;
676  }
677  else if ( pZ > 0 && pA > 1 ) projPDG = 90000000+999*pZ+pA;
678  else if (particle == G4Lambda::Lambda() ) projPDG= 3122;
679  else if (particle == G4SigmaPlus::SigmaPlus() ) projPDG= 3222;
680  else if (particle == G4SigmaMinus::SigmaMinus() ) projPDG= 3112;
681  else if (particle == G4SigmaZero::SigmaZero() ) projPDG= 3212;
682  else if (particle == G4XiMinus::XiMinus() ) projPDG= 3312;
683  else if (particle == G4XiZero::XiZero() ) projPDG= 3322;
684  else if (particle == G4OmegaMinus::OmegaMinus() ) projPDG= 3334;
685  else if (particle == G4AntiNeutron::AntiNeutron() ) projPDG=-2112;
686  else if (particle == G4AntiProton::AntiProton() ) projPDG=-2212;
687  else if (particle == G4MuonPlus::MuonPlus() ) projPDG= -13;
688  else if (particle == G4MuonMinus::MuonMinus() ) projPDG= 13;
689  else if (particle == G4Gamma::Gamma() ) projPDG= 22;
690  else if (particle == G4Electron::Electron() ) projPDG= 11;
691  else if (particle == G4Positron::Positron() ) projPDG= -11;
692  else if (particle == G4NeutrinoMu::NeutrinoMu() ) projPDG= 14;
693  else if (particle == G4AntiNeutrinoMu::AntiNeutrinoMu() ) projPDG= -14;
694  else if (particle == G4AntiLambda::AntiLambda() ) projPDG=-3122;
695  else if (particle == G4AntiSigmaPlus::AntiSigmaPlus() ) projPDG=-3222;
696  else if (particle == G4AntiSigmaMinus::AntiSigmaMinus() ) projPDG=-3112;
697  else if (particle == G4AntiSigmaZero::AntiSigmaZero() ) projPDG=-3212;
698  else if (particle == G4AntiXiMinus::AntiXiMinus() ) projPDG=-3312;
699  else if (particle == G4AntiXiZero::AntiXiZero() ) projPDG=-3322;
700  else if (particle == G4AntiOmegaMinus::AntiOmegaMinus() ) projPDG=-3334;
701  else if (particle == G4NeutrinoE::NeutrinoE() ) projPDG= 12;
702  else if (particle == G4AntiNeutrinoE::AntiNeutrinoE() ) projPDG= -12;
703  else if (particle == G4TauPlus::TauPlus() ) projPDG= -15;
704  else if (particle == G4TauMinus::TauMinus() ) projPDG= 15;
705  else if (particle == G4NeutrinoTau::NeutrinoTau() ) projPDG= 16;
706  else if (particle == G4AntiNeutrinoTau::AntiNeutrinoTau()) projPDG= -16;
707  G4int aProjPDG=std::abs(projPDG);
708 #ifdef debug
709  G4int prPDG=particle->GetPDGEncoding();
710  G4cout<<"G4QInelastic::PostStepDoIt: projPDG="<<projPDG<<", stPDG="<<prPDG<<G4endl;
711 #endif
712  if(!projPDG)
713  {
714  G4cerr<<"---Warning---G4QInelastic::PostStepDoIt:Undefined interacting hadron"<<G4endl;
715  return 0;
716  }
717  G4int EPIM=ElProbInMat.size();
718 #ifdef debug
719  G4cout<<"G4QInel::PostStDoIt: m="<<EPIM<<",n="<<nE<<",T="<<ElProbInMat[EPIM-1]<<G4endl;
720 #endif
721  G4int i=0;
722  if(EPIM>1)
723  {
724  G4double rnd = ElProbInMat[EPIM-1]*G4UniformRand();
725  for(i=0; i<nE; ++i)
726  {
727 #ifdef debug
728  G4cout<<"G4QInelastic::PostStepDoIt:E["<<i<<"]="<<ElProbInMat[i]<<",r="<<rnd<<G4endl;
729 #endif
730  if (rnd<ElProbInMat[i]) break;
731  }
732  if(i>=nE) i=nE-1; // Top limit for the Element
733  }
734  G4Element* pElement=(*theElementVector)[i];
735  Z=static_cast<G4int>(pElement->GetZ());
736 #ifdef debug
737  G4cout<<"G4QInelastic::PostStepDoIt: i="<<i<<", Z(element)="<<Z<<G4endl;
738 #endif
739  if(Z<=0)
740  {
741  G4cerr<<"---Warning---G4QInelastic::PostStepDoIt: Element with Z="<<Z<<G4endl;
742  if(Z<0) return 0;
743  }
744  std::vector<G4double>* SPI = IsoProbInEl[i];// Vector of summedProbabilities for isotopes
745  std::vector<G4int>* IsN = ElIsoN[i]; // Vector of "#of neutrons" in the isotope El[i]
746  G4int nofIsot=SPI->size(); // #of isotopes in the element i
747 #ifdef debug
748  G4cout<<"G4QInel::PosStDoIt:n="<<nofIsot<<",T="<<(*SPI)[nofIsot-1]<<G4endl;
749 #endif
750  G4int j=0;
751  if(nofIsot>1)
752  {
753  G4double rndI=(*SPI)[nofIsot-1]*G4UniformRand(); // Randomize the isotop of the Element
754  for(j=0; j<nofIsot; ++j)
755  {
756 #ifdef debug
757  G4cout<<"G4QInelastic::PostStepDoIt: SP["<<j<<"]="<<(*SPI)[j]<<", r="<<rndI<<G4endl;
758 #endif
759  if(rndI < (*SPI)[j]) break;
760  }
761  if(j>=nofIsot) j=nofIsot-1; // Top limit for the isotope
762  }
763  G4int N =(*IsN)[j]; ; // Randomized number of neutrons
764 #ifdef debug
765  G4cout<<"G4QInelastic::PostStepDoIt: Z="<<Z<<", j="<<i<<", N(isotope)="<<N<<G4endl;
766 #endif
767  G4double kinEnergy= projHadron->GetKineticEnergy();
768  G4ParticleMomentum dir = projHadron->GetMomentumDirection();
769  if(projPDG==22 && Z==1 && !N && Momentum<150.)
770  {
771 #ifdef debug
772  G4cerr<<"---LowPHOTON---G4QInelastic::PSDoIt: Called for zero Cross-section"<<G4endl;
773 #endif
774  //Do Nothing Action insead of the reaction
775  aParticleChange.ProposeEnergy(kinEnergy);
776  aParticleChange.ProposeLocalEnergyDeposit(0.);
777  aParticleChange.ProposeMomentumDirection(dir);
778  aParticleChange.ProposeTrackStatus(fAlive);
779  return G4VDiscreteProcess::PostStepDoIt(track,step);
780  }
781  if(N<0)
782  {
783  G4cerr<<"-Warning-G4QInelastic::PostStepDoIt: Isotope with Z="<<Z<<", 0>N="<<N<<G4endl;
784  return 0;
785  }
786  nOfNeutrons=N; // Remember it for the energy-momentum check
787  //G4double am=Z+N;
788  G4int am=Z+N;
789  //G4double dd=0.025;
790  //G4double sr=std::sqrt(am);
791  //G4double dsr=0.01*(sr+sr);
792  //if(dsr<dd)dsr=dd;
793  //G4double medRA=mediRatio*G4QThd(am);
794  G4double medRA=mediRatio;
795  if(manualFlag) G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,medRA); // ManualP
796  //else if(projPDG==-2212)G4QNucleus::SetParameters(1.-dsr-dsr,dd+dd,5.,9.);//aP ClustPars
797  //else if(projPDG==-211) G4QNucleus::SetParameters(.67-dsr,.32-dsr,5.,9.); //Pi-ClustPars
798  //else if(projPDG ==-2212 || projPDG ==-2112)
799  // G4QNucleus::SetParameters(1.-dsr-dsr,dd+dd,5.,1.);//aP ClustPars
800  //else if(projPDG ==-211 ||projPDG == 211 )
801  // G4QNucleus::SetParameters(.67-dsr,.32-dsr,5.,1.); //Pi-ClustPars
802 #ifdef debug
803  G4cout<<"G4QInelastic::PostStepDoIt: N="<<N<<" for element with Z="<<Z<<G4endl;
804 #endif
805  aParticleChange.Initialize(track);
806  G4double weight = track.GetWeight();
807 #ifdef debug
808  G4cout<<"G4QInelastic::PostStepDoIt: weight="<<weight<<G4endl;
809 #endif
810  if(photNucBias!=1.) weight/=photNucBias;
811  else if(weakNucBias!=1.) weight/=weakNucBias;
812  G4double localtime = track.GetGlobalTime();
813 #ifdef debug
814  G4cout<<"G4QInelastic::PostStepDoIt: localtime="<<localtime<<G4endl;
815 #endif
816  G4ThreeVector position = track.GetPosition();
817  G4TouchableHandle trTouchable = track.GetTouchableHandle();
818 #ifdef debug
819  G4cout<<"G4QInelastic::PostStepDoIt: position="<<position<<G4endl;
820 #endif
821  //
822  G4int targPDG=90000000+Z*1000+N; // PDG Code of the target nucleus
823  G4QPDGCode targQPDG(targPDG);
824  G4double tgM=targQPDG.GetMass(); // Target mass
825  G4double tM=tgM; // Target mass (copy to be changed)
826  G4QHadronVector* output=new G4QHadronVector;// Prototype of Fragm Output G4QHadronVector
827  G4double absMom = 0.; // Prototype of absorbed by nucleus Moment
828  G4QHadronVector* leadhs=new G4QHadronVector;// Prototype of QuasmOutput G4QHadronVector
829  G4LorentzVector lead4M(0.,0.,0.,0.); // Prototype of LeadingQ 4-momentum
830 #ifdef debug
831  G4cout<<"G4QInelastic::PostStepDoIt: projPDG="<<aProjPDG<<", targPDG="<<targPDG<<G4endl;
832 #endif
833  //
834  // Leptons with photonuclear
835  // Lepto-nuclear case with the equivalent photon algorithm. @@InFuture + NC (?)
836  //
837  if (aProjPDG == 11 || aProjPDG == 13 || aProjPDG == 15)
838  {
839 #ifdef debug
840  G4cout<<"G4QInelastic::PostStDoIt:startSt="<<aParticleChange.GetTrackStatus()<<G4endl;
841 #endif
842  G4VQCrossSection* CSmanager=G4QElectronNuclearCrossSection::GetPointer();
843  G4double ml=me;
844  G4double ml2=me2;
845  if(aProjPDG== 13)
846  {
847  CSmanager=G4QMuonNuclearCrossSection::GetPointer();
848  ml=mu;
849  ml2=mu2;
850  }
851  if(aProjPDG== 15)
852  {
853  CSmanager=G4QTauNuclearCrossSection::GetPointer();
854  ml=mt;
855  ml2=mt2;
856  }
857  // @@ Probably this is not necessary any more (?)
858  if(!aProjPDG) G4cout<<"-Warning-G4QInelast::PostStepDoIt:(2) projectile PDG=0"<<G4endl;
859  G4double xSec=CSmanager->GetCrossSection(false,Momentum,Z,N,aProjPDG);// Recalculate XS
860  // @@ check a possibility to separate p, n, or alpha (!)
861  if(Z==1 && !N && Momentum<150.) xSec=0.;
862 #ifdef debug
863  G4cout<<"-Forse-G4QInel::PStDoIt:P="<<Momentum<<",PDG="<<projPDG<<",xS="<<xSec<<G4endl;
864 #endif
865  if(xSec <= 0.) // The cross-section is 0 -> Do Nothing
866  {
867 #ifdef debug
868  G4cerr<<"---OUT---G4QInelastic::PSDoIt: Called for zero Cross-section"<<G4endl;
869 #endif
870  //Do Nothing Action insead of the reaction
871  aParticleChange.ProposeEnergy(kinEnergy);
872  aParticleChange.ProposeLocalEnergyDeposit(0.);
873  aParticleChange.ProposeMomentumDirection(dir);
874  aParticleChange.ProposeTrackStatus(fAlive);
875  delete output;
876  return G4VDiscreteProcess::PostStepDoIt(track,step);
877  }
878  G4double photonEnergy = CSmanager->GetExchangeEnergy(); // Energy of EqivExchangePart
879 #ifdef debug
880  G4cout<<"G4QInel::PStDoIt:kE="<<kinEnergy<<",dir="<<dir<<",phE="<<photonEnergy<<G4endl;
881 #endif
882  if( kinEnergy < photonEnergy || photonEnergy < 0.)
883  {
884  //Do Nothing Action insead of the reaction
885  G4cerr<<"--G4QInelastic::PSDoIt: photE="<<photonEnergy<<">leptE="<<kinEnergy<<G4endl;
886  aParticleChange.ProposeEnergy(kinEnergy);
887  aParticleChange.ProposeLocalEnergyDeposit(0.);
888  aParticleChange.ProposeMomentumDirection(dir);
889  aParticleChange.ProposeTrackStatus(fAlive);
890  delete output;
891  return G4VDiscreteProcess::PostStepDoIt(track,step);
892  }
893  G4double photonQ2 = CSmanager->GetExchangeQ2(photonEnergy);// Q2(t) of EqivExchangePart
894  G4double W=photonEnergy-photonQ2/dM;// HadronicEnergyFlow (W-energy) for virtual photon
895  if(tM<999.) W-=mPi0+mPi0s/dM; // Pion production threshold for a nucleon target
896  if(W<0.)
897  {
898  //Do Nothing Action insead of the reaction
899 #ifdef debug
900  G4cout<<"--G4QInelastic::PostStepDoIt:(lN) negative equivalent energy W="<<W<<G4endl;
901 #endif
902  aParticleChange.ProposeEnergy(kinEnergy);
903  aParticleChange.ProposeLocalEnergyDeposit(0.);
904  aParticleChange.ProposeMomentumDirection(dir);
905  aParticleChange.ProposeTrackStatus(fAlive);
906  delete output;
907  return G4VDiscreteProcess::PostStepDoIt(track,step);
908  }
909  // Update G4VParticleChange for the scattered muon
910  G4VQCrossSection* thePhotonData=G4QPhotonNuclearCrossSection::GetPointer();
911  G4double sigNu=thePhotonData->GetCrossSection(true,photonEnergy,Z,N,22);//Integrated XS
912  G4double sigK =thePhotonData->GetCrossSection(true, W, Z, N, 22); // Real XS
913  G4double rndFraction = CSmanager->GetVirtualFactor(photonEnergy, photonQ2);
914  if(sigNu*G4UniformRand()>sigK*rndFraction)
915  {
916  //NothingToDo Action insead of the reaction
917 #ifdef debug
918  G4cout<<"-DoNoth-G4QInelastic::PostStepDoIt: probab. correction - DoNothing"<<G4endl;
919 #endif
920  aParticleChange.ProposeEnergy(kinEnergy);
921  aParticleChange.ProposeLocalEnergyDeposit(0.);
922  aParticleChange.ProposeMomentumDirection(dir);
923  aParticleChange.ProposeTrackStatus(fAlive);
924  delete output;
925  return G4VDiscreteProcess::PostStepDoIt(track,step);
926  }
927  G4double iniE=kinEnergy+ml; // Initial total energy of the lepton
928  G4double finE=iniE-photonEnergy; // Final total energy of the lepton
929 #ifdef pdebug
930  G4cout<<"G4QInelastic::PoStDoIt:E="<<iniE<<",lE="<<finE<<"-"<<ml<<"="<<finE-ml<<G4endl;
931 #endif
932  aParticleChange.ProposeEnergy(finE-ml);
933  if(finE<=ml) // Secondary lepton (e/mu/tau) at rest disappears
934  {
935  aParticleChange.ProposeEnergy(0.);
936  if(aProjPDG== 11) aParticleChange.ProposeTrackStatus(fStopAndKill);
937  else aParticleChange.ProposeTrackStatus(fStopButAlive);
938  aParticleChange.ProposeMomentumDirection(dir);
939  }
940  else aParticleChange.ProposeTrackStatus(fAlive);
941  G4double iniP=std::sqrt(iniE*iniE-ml2); // Initial momentum of the electron
942  G4double finP=std::sqrt(finE*finE-ml2); // Final momentum of the electron
943  G4double cost=(iniE*finE-ml2-photonQ2/2)/iniP/finP; // cos(scat_ang_of_lepton)
944 #ifdef pdebug
945  G4cout<<"G4QI::PSDoIt:Q2="<<photonQ2<<",ct="<<cost<<",Pi="<<iniP<<",Pf="<<finP<<G4endl;
946 #endif
947  if(cost>1.) cost=1.; // To avoid the accuracy of calculation problem
948  if(cost<-1.) cost=-1.; // To avoid the accuracy of calculation problem
949  //
950  // Scatter the lepton ( @@ make the same thing for real photons)
951  // At this point we have photonEnergy and photonQ2 (with notDefinedPhi)->SelectProjPart
952  G4double absEn = G4QThd(am)*GeV; // @@(b) Mean Energy Absorbed by a Nucleus
953  //G4double absEn = std::pow(am,third)*GeV; // @@(b) Mean Energy Absorbed by a Nucleus
954  if(am>1 && absEn < photonEnergy) // --> the absorption of energy can happen
955  //if(absEn < photonEnergy) // --> the absorption of energy can happen
956  {
957  G4double abtEn = absEn+hdM; // @@(b) MeanEnergyAbsorbed by a nucleus (+M_N)
958  G4double abEn2 = abtEn*abtEn; // Squared absorbed Energy + MN
959  G4double abMo2 = abEn2-hdM2; // Squared absorbed Momentum of compound system
960  G4double phEn2 = photonEnergy*photonEnergy;
961  G4double phMo2 = phEn2+photonQ2; // Squared momentum of primary virtual photon
962  G4double phMo = std::sqrt(phMo2); // Momentum of the primary virtual photon
963  absMom = std::sqrt(abMo2); // Absorbed Momentum
964  if(absMom < phMo) // --> the absorption of momentum can happen
965  {
966  G4double dEn = photonEnergy - absEn; // Leading energy
967  G4double dMo = phMo - absMom; // Leading momentum
968  G4double sF = dEn*dEn - dMo*dMo;// s of leading particle
969 #ifdef ppdebug
970  G4cout<<"-PhotoAbsorption-G4QIn::PStDoIt: sF="<<sF<<",phEn="<<photonEnergy<<G4endl;
971 #endif
972  if(sF > stmPi) // --> Leading fragmentation is possible
973  {
974  photonEnergy = absEn; // New value of the photon energy
975  photonQ2=abMo2-absEn*absEn; // New value of the photon Q2
976  absEn = dEn; // Put energy of leading particle to absEn (!)
977  }
978  else absMom=0.; // Flag that nothing has happened
979  }
980  else absMom=0.; // Flag that nothing has happened
981  }
982  // ------------- End of ProjPart selection
983  //
984  // Scattering in respect to the derection of the incident muon is made impicitly:
985  G4ThreeVector ort=dir.orthogonal(); // Not normed orthogonal vector (!) (to dir)
986  G4ThreeVector ortx = ort.unit(); // First unit vector orthogonal to the direction
987  G4ThreeVector orty = dir.cross(ortx);// Second unit vector orthoganal to the direction
988  G4double sint=std::sqrt(1.-cost*cost); // Perpendicular component
989  G4double phi=twopi*G4UniformRand(); // phi of scattered electron
990  G4double sinx=sint*std::sin(phi); // x perpendicular component
991  G4double siny=sint*std::cos(phi); // y perpendicular component
992  G4ThreeVector findir=cost*dir+sinx*ortx+siny*orty;
993  aParticleChange.ProposeMomentumDirection(findir); // new direction for the lepton
994 #ifdef pdebug
995  G4cout<<"G4QInelastic::PostStepDoIt: E="<<aParticleChange.GetEnergy()<<"="<<finE<<"-"
996  <<ml<<", d="<<*aParticleChange.GetMomentumDirection()<<","<<findir<<G4endl;
997 #endif
998  G4ThreeVector photon3M=iniP*dir-finP*findir;// 3D total momentum of photon
999  if(absMom) // Photon must be reduced & LeadingSyst fragmented
1000  {
1001  G4double ptm=photon3M.mag(); // 3M of the virtual photon
1002 #ifdef ppdebug
1003  G4cout<<"-Absorption-G4QInelastic::PostStepDoIt: ph3M="<<photon3M<<", eIn3M="
1004  <<iniP*dir<<", eFin3M="<<finP*findir<<", abs3M="<<absMom<<"<ptm="<<ptm<<G4endl;
1005 #endif
1006  G4ThreeVector lead3M=photon3M*(ptm-absMom)/ptm; // Keep the direction for leading Q
1007  photon3M-=lead3M; // Reduced photon Momentum (photEn already = absEn)
1008  proj4M=G4LorentzVector(lead3M,absEn); // 4-momentum of leading System
1009 #ifdef ppdebug
1010  G4cout<<"-->G4QI::PoStDoIt: new sF="<<proj4M.m2()<<", lead4M="<<proj4M<<G4endl;
1011 #endif
1012  lead4M=proj4M; // Remember 4-mom for the total 4-momentum
1013  G4Quasmon* pan= new G4Quasmon(G4QContent(1,1,0,1,1,0),proj4M);// ---> DELETED -->---+
1014  try // |
1015  { // |
1016  if(leadhs) delete leadhs; // |
1017  G4QNucleus vac(90000000); // |
1018  leadhs=pan->Fragment(vac,1); // DELETED after it is copied to output vector |
1019  } // |
1020  catch (G4QException& error) // |
1021  { // |
1022  G4cerr<<"***G4QInelastic::PostStepDoIt: G4Quasmon Exception is catched"<<G4endl;//|
1023  // G4Exception("G4QInelastic::PostStepDoIt:","72",FatalException,"QuasmonCrash"); //|
1024  G4Exception("G4QInelastic::PostStepDoIt()","HAD_CHPS_0072",
1025  FatalException, "QuasmonCrash");
1026  } // |
1027  delete pan; // Delete the Nuclear Environment <----<---+
1028 #ifdef ppdebug
1029  G4cout<<"G4QInel::PStDoIt: l4M="<<proj4M<<proj4M.m2()<<", N="<<leadhs->size()<<",pt="
1030  <<ptm<<",pa="<<absMom<<",El="<<absEn<<",Pl="<<ptm-absMom<<G4endl;
1031 #endif
1032  }
1033  else
1034  {
1035  G4int qNH=0;
1036  if(leadhs) qNH=leadhs->size();
1037  if(qNH) for(G4int iq=0; iq<qNH; iq++) delete (*leadhs)[iq];
1038  if(leadhs) delete leadhs;
1039  leadhs=0;
1040  }
1041  projPDG=22;
1042  proj4M=G4LorentzVector(photon3M,photonEnergy);
1043 #ifdef debug
1044  G4cout<<"G4QInelastic::PostStDoIt: St="<<aParticleChange.GetTrackStatus()<<", g4m="
1045  <<proj4M<<", lE="<<finE<<", lP="<<finP*findir<<", d="<<findir.mag2()<<G4endl;
1046 #endif
1047 
1048  //
1049  // neutrinoNuclear interactions (nu_e, nu_mu only)
1050  //
1051  }
1052  else if (aProjPDG == 12 || aProjPDG == 14)
1053  {
1054  kinEnergy= projHadron->GetKineticEnergy()/MeV; // Total energy of the neutrino
1055  G4double dKinE=kinEnergy+kinEnergy; // doubled energy for s calculation
1056 #ifdef debug
1057  G4cout<<"G4QInelastic::PostStDoIt: 2*nuEnergy="<<dKinE<<"(MeV), PDG="<<projPDG<<G4endl;
1058 #endif
1059  dir = projHadron->GetMomentumDirection(); // unit vector
1060  G4double ml = mu;
1061  G4double ml2 = mu2;
1062  //G4double mlN = muN;
1063  G4double mlN2= muN2;
1064  G4double fmlN= fmuN;
1065  G4double mlsN= musN;
1066  //G4double mlP = muP;
1067  G4double mlP2= muP2;
1068  G4double fmlP= fmuP;
1069  G4double mlsP= musP;
1070  G4double mldM= mudM;
1071  //G4double mlD = muD;
1072  G4double mlD2= muD2;
1073  if(aProjPDG==12)
1074  {
1075  ml = me;
1076  ml2 = me2;
1077  //mlN = meN;
1078  mlN2= meN2;
1079  fmlN= fmeN;
1080  mlsN= mesN;
1081  //mlP = meP;
1082  mlP2= meP2;
1083  fmlP= fmeP;
1084  mlsP= mesP;
1085  mldM= medM;
1086  //mlD = meD;
1087  mlD2= meD2;
1088  }
1089  G4VQCrossSection* CSmanager =G4QNuMuNuclearCrossSection::GetPointer(); // (nu,l)
1090  G4VQCrossSection* CSmanager2=G4QNuNuNuclearCrossSection::GetPointer(); // (nu,nu)
1091  proj4M=G4LorentzVector(dir*kinEnergy,kinEnergy); // temporary
1092  G4bool nuanu=true;
1093  scatPDG=13; // Prototype = secondary scattered mu-
1094  if(projPDG==-14)
1095  {
1096  nuanu=false; // Anti-neutrino
1097  CSmanager=G4QANuMuNuclearCrossSection::GetPointer(); // (anu,mu+) CC @@ open
1098  CSmanager=G4QANuANuNuclearCrossSection::GetPointer(); // (anu,anu) NC @@ open
1099  scatPDG=-13; // secondary scattered mu+
1100  }
1101  else if(projPDG==12)
1102  {
1103  CSmanager=G4QNuENuclearCrossSection::GetPointer(); // @@ open (only CC is changed)
1104  scatPDG=11; // secondary scattered e-
1105  }
1106  else if(projPDG==-12)
1107  {
1108  nuanu=false; // anti-neutrino
1109  CSmanager=G4QANuENuclearCrossSection::GetPointer(); // (anu,e+) CC @@ open
1110  CSmanager=G4QANuANuNuclearCrossSection::GetPointer(); // (anu,anu) NC @@ open
1111  scatPDG=-11; // secondary scattered e+
1112  }
1113  // @@ Probably this is not necessary any more
1114  if(!projPDG) G4cout<<"-Warning-G4QInelastic::PostStepDoIt:(3)projectile PDG=0"<<G4endl;
1115  G4double xSec1=CSmanager->GetCrossSection(false,Momentum,Z,N,projPDG); //Recalculate XS
1116  G4double xSec2=CSmanager2->GetCrossSection(false,Momentum,Z,N,projPDG);//Recalculate XS
1117  G4double xSec=xSec1+xSec2;
1118  // @@ check a possibility to separate p, n, or alpha (!)
1119  if(xSec <= 0.) // The cross-section = 0 -> Do Nothing
1120  {
1121  G4cerr<<"G4QInelastic::PSDoIt:nuE="<<kinEnergy<<",X1="<<xSec1<<",X2="<<xSec2<<G4endl;
1122  //Do Nothing Action insead of the reaction
1123  aParticleChange.ProposeEnergy(kinEnergy);
1124  aParticleChange.ProposeLocalEnergyDeposit(0.);
1125  aParticleChange.ProposeMomentumDirection(dir);
1126  aParticleChange.ProposeTrackStatus(fAlive);
1127  delete output;
1128  return G4VDiscreteProcess::PostStepDoIt(track,step);
1129  }
1130  G4bool secnu=false;
1131  if(xSec*G4UniformRand()>xSec1) // recover neutrino/antineutrino
1132  {
1133  if(scatPDG>0) scatPDG++;
1134  else scatPDG--;
1135  secnu=true;
1136  }
1137  scat=true; // event with changed scattered projectile
1138  G4double totCS1 = CSmanager->GetLastTOTCS(); // the last total cross section1(isotope?)
1139  G4double totCS2 = CSmanager2->GetLastTOTCS();// the last total cross section2(isotope?)
1140  G4double totCS = totCS1+totCS2; // the last total cross section (isotope?)
1141  if(std::fabs(xSec-totCS*millibarn)/xSec>.0001)
1142  G4cout<<"-Warning-G4QInelastic::PostStepDoIt: xS="<<xSec<<"# CS="<<totCS<<G4endl;
1143  G4double qelCS1 = CSmanager->GetLastQELCS(); // the last quasi-elastic cross section1
1144  G4double qelCS2 = CSmanager2->GetLastQELCS();// the last quasi-elastic cross section2
1145  G4double qelCS = qelCS1+qelCS2; // the last quasi-elastic cross section
1146  if(totCS - qelCS < 0.) // only at low energies
1147  {
1148  totCS = qelCS;
1149  totCS1 = qelCS1;
1150  totCS2 = qelCS2;
1151  }
1152  // make different definitions for neutrino and antineutrino
1153  G4double mIN=mProt; // Just a prototype (for anu, Z=1, N=0)
1154  G4double mOT=mNeut;
1155  G4double OT=mlN2;
1156  //G4double mOT2=mNeut2; // Other formula: sd=s-mlsOT -Not used?-
1157  G4double mlOT=fmlN;
1158  G4double mlsOT=mlsN;
1159  if(secnu)
1160  {
1161  if(am*G4UniformRand()>Z) // Neutron target
1162  {
1163  targPDG-=1; // subtract neutron
1164  projPDG=2112; // neutron is going out
1165  mIN =mNeut;
1166  OT =mNeut2;
1167  //mOT2=mNeut2;
1168  mlOT=0.;
1169  mlsOT=mNeut2;
1170  }
1171  else
1172  {
1173  targPDG-=1000; // subtract neutron
1174  projPDG=2212; // neutron is going out
1175  mOT =mProt;
1176  OT =mProt2;
1177  //mOT2 =mProt2;
1178  mlOT =0.;
1179  mlsOT=mProt2;
1180  }
1181  ml=0.;
1182  ml2=0.;
1183  mldM=0.;
1184  mlD2=mPPi2;
1185  G4QPDGCode temporary_targQPDG(targPDG);
1186  targQPDG = temporary_targQPDG;
1187  G4double rM=targQPDG.GetMass();
1188  mIN=tM-rM; // bounded in-mass of the neutron
1189  tM=rM;
1190  }
1191  else if(nuanu)
1192  {
1193  targPDG-=1; // Neutrino -> subtract neutron
1194  G4QPDGCode temporary_targQPDG(targPDG);
1195  targQPDG = temporary_targQPDG;
1196  G4double rM=targQPDG.GetMass();
1197  mIN=tM-rM; // bounded in-mass of the neutron
1198  tM=rM;
1199  mOT=mProt;
1200  OT=mlP2;
1201  //mOT2=mProt2;
1202  mlOT=fmlP;
1203  mlsOT=mlsP;
1204  projPDG=2212; // proton is going out
1205  }
1206  else
1207  {
1208  if(Z>1||N>0) // Calculate the splitted mass
1209  {
1210  targPDG-=1000; // Anti-Neutrino -> subtract proton
1211  G4QPDGCode temporary_targQPDG(targPDG);
1212  targQPDG = temporary_targQPDG;
1213  G4double rM=targQPDG.GetMass();
1214  mIN=tM-rM; // bounded in-mass of the proton
1215  tM=rM;
1216  }
1217  else
1218  {
1219  targPDG=0;
1220  mIN=tM;
1221  tM=0.;
1222  }
1223  projPDG=2112; // neutron is going out
1224  }
1225  G4double s_value=mIN*(mIN+dKinE); // s_value=(M_cm)^2=m2+2mE (m=targetMass,E=E_nu)
1226 #ifdef debug
1227  G4cout<<"G4QInelastic::PostStDoIt: s="<<s_value<<" >? OT="<<OT<<", mlD2="<<mlD2<<G4endl;
1228 #endif
1229  if(s_value<=OT) // *** Do nothing solution ***
1230  {
1231  //Do NothingToDo Action insead of the reaction (@@ Can we make it common?)
1232  G4cout<<"G4QInelastic::PostStepDoIt: tooSmallFinalMassOfCompound: DoNothing"<<G4endl;
1233  aParticleChange.ProposeEnergy(kinEnergy);
1234  aParticleChange.ProposeLocalEnergyDeposit(0.);
1235  aParticleChange.ProposeMomentumDirection(dir);
1236  aParticleChange.ProposeTrackStatus(fAlive);
1237  delete output;
1238  return G4VDiscreteProcess::PostStepDoIt(track,step);
1239  }
1240 #ifdef debug
1241  G4cout<<"G4QInelastic::PostStDoIt: Stop and kill the projectile neutrino"<<G4endl;
1242 #endif
1243  aParticleChange.ProposeEnergy(0.);
1244  aParticleChange.ProposeTrackStatus(fStopAndKill); // the initial neutrino is killed
1245  // There is no way back from here !
1246  if ( ((secnu || !nuanu || N) && totCS*G4UniformRand() < qelCS) || s_value < mlD2 )
1247  { // Quasi-Elastic interaction
1248  G4double Q2=0.; // Simulate transferred momentum, in MeV^2
1249  if(secnu) Q2=CSmanager2->GetQEL_ExchangeQ2();
1250  else Q2=CSmanager->GetQEL_ExchangeQ2();
1251 #ifdef debug
1252  G4cout<<"G4QInelastic::PostStDoIt:QuasiEl(nu="<<secnu<<"),s="<<s_value<<",Q2="<<Q2<<G4endl;
1253 #endif
1254  //G4double ds=s_value+s_value; // doubled s_value
1255  G4double sqs=std::sqrt(s_value); // M_cm
1256  G4double dsqs=sqs+sqs; // 2*M_cm
1257  G4double p_init=(s_value-mIN*mIN)/dsqs; // initial momentum in CMS (checked MK)
1258  G4double dpi=p_init+p_init; // doubled initial momentum in CMS
1259  G4double sd=s_value-mlsOT; // s_value-ml2-mOT2 (mlsOT=m^2_neut+m^2_lept)
1260  G4double qo2=(sd*sd-mlOT)/dsqs; // squared momentum of secondaries in CMS
1261  G4double qo=std::sqrt(qo2); // momentum of secondaries in CMS
1262  G4double cost=(dpi*std::sqrt(qo2+ml2)-Q2-ml2)/dpi/qo; // cos(theta) in CMS (chck MK)
1263  G4LorentzVector t4M(0.,0.,0.,mIN); // 4mom of the effective target
1264  G4LorentzVector c4M=t4M+proj4M; // 4mom of the compound system
1265  t4M.setT(mOT); // now it is 4mom of the outgoing nucleon
1266  scat4M=G4LorentzVector(0.,0.,0.,ml); // 4mom of the scattered muon
1267  if(!G4QHadron(c4M).RelDecayIn2(scat4M, t4M, proj4M, cost, cost))
1268  {
1269  G4cerr<<"G4QIn::PStD:c4M="<<c4M<<sqs<<",mM="<<ml<<",tM="<<mOT<<",c="<<cost<<G4endl;
1270  // throw G4QException("G4QInelastic::HadronizeQuasm: Can't dec QE nu,lept Compound");
1271  G4Exception("G4QInelastic::PostStepDoIt()", "HAD_CHPS_0000",
1272  FatalException, "Hadronize quasmon: Can't dec QE nu,lept Compound");
1273  }
1274  proj4M=t4M; // 4mom of the new projectile nucleon
1275  }
1276  else // ***** Non Quasi Elastic interaction
1277  {
1278  // Recover the target PDG Code if the quasi-elastic scattering did not happened
1279  if ( (secnu && projPDG == 2212) || (!secnu && projPDG == 2112) ) targPDG+=1;
1280  else if ( (secnu && projPDG == 2112) || (!secnu && projPDG == 2212) ) targPDG+=1000;
1281  G4double Q2=0; // Simulate transferred momentum, in MeV^2
1282  if(secnu) Q2=CSmanager->GetNQE_ExchangeQ2();
1283  else Q2=CSmanager2->GetNQE_ExchangeQ2();
1284 #ifdef debug
1285  G4cout<<"G4QInel::PStDoIt: MultiPeriferal s="<<s_value<<",Q2="<<Q2<<",T="<<targPDG<<G4endl;
1286 #endif
1287  if(secnu) projPDG=CSmanager2->GetExchangePDGCode();// PDG Code of the effective gamma
1288  else projPDG=CSmanager->GetExchangePDGCode(); // PDG Code of the effective pion
1289  //@@ Temporary made only for direct interaction and for N=3 (good for small Q2)
1290  //@@ inFuture use N=GetNPartons and directFraction=GetDirectPart, @@ W2...
1291  G4double r=G4UniformRand();
1292  G4double r1=0.5; // (1-x)
1293  if(r<0.5) r1=std::sqrt(r+r)*(.5+.1579*(r-.5));
1294  else if(r>0.5) r1=1.-std::sqrt(2.-r-r)*(.5+.1579*(.5-r));
1295  G4double xn=1.-mldM/Momentum; // Normalization of (1-x) [x>mldM/Mom]
1296  G4double x1=xn*r1; // (1-x)
1297  G4double x=1.-x1; // x=2k/M
1298  //G4double W2=(hdM2+Q2/x)*x1; // W2 candidate
1299  G4double mx=hdM*x; // Part of the target to interact with
1300  G4double we=Q2/(mx+mx); // transfered energy
1301  if(we>=kinEnergy-ml-.001) we=kinEnergy-ml-.0001; // safety to avoid nan=sqrt(neg)
1302  G4double pot=kinEnergy-we; // energy of the secondary lepton
1303  G4double mlQ2=ml2+Q2;
1304  G4double cost=(pot-mlQ2/dKinE)/std::sqrt(pot*pot-ml2); // LS cos(theta)
1305  if(std::fabs(cost)>1)
1306  {
1307 #ifdef debug
1308  G4cout<<"*G4QInelastic::PostStDoIt: cost="<<cost<<", Q2="<<Q2<<", nu="<<we<<", mx="
1309  <<mx<<", pot="<<pot<<", 2KE="<<dKinE<<G4endl;
1310 #endif
1311  if(cost>1.) cost=1.;
1312  else cost=-1.;
1313  pot=mlQ2/dKinE+dKinE*ml2/mlQ2; // extreme output momentum
1314  }
1315  G4double lEn=std::sqrt(pot*pot+ml2); // Lepton energy
1316  G4double lPl=pot*cost; // Lepton longitudinal momentum
1317  G4double lPt=pot*std::sqrt(1.-cost*cost); // Lepton transverse momentum
1318  std::pair<G4double,G4double> d2d=Random2DDirection(); // Randomize phi
1319  G4double lPx=lPt*d2d.first;
1320  G4double lPy=lPt*d2d.second;
1321  G4ThreeVector vdir=proj4M.vect(); // 3D momentum of the projectile
1322  G4ThreeVector vz= vdir.unit(); // Ort in the direction of the projectile
1323  G4ThreeVector vv= vz.orthogonal(); // Not normed orthogonal vector (!)
1324  G4ThreeVector vx= vv.unit(); // First ort orthogonal to the direction
1325  G4ThreeVector vy= vz.cross(vx); // Second ort orthoganal to the direction
1326  G4ThreeVector lP= lPl*vz+lPx*vx+lPy*vy; // 3D momentum of the scattered lepton
1327  scat4M=G4LorentzVector(lP,lEn); // 4mom of the scattered lepton
1328  proj4M-=scat4M; // 4mom of the W/Z (effective pion/gamma)
1329 #ifdef debug
1330  G4cout<<"G4QInelastic::PostStDoIt: proj4M="<<proj4M<<", ml="<<ml<<G4endl;
1331 #endif
1332  // Check that the en/mom transfer is possible, if not -> elastic
1333  G4int fintPDG=targPDG; // Prototype for the compound nucleus
1334  if(!secnu)
1335  {
1336  if(projPDG<0) fintPDG-= 999;
1337  else fintPDG+= 999;
1338  }
1339  G4double fM=G4QPDGCode(fintPDG).GetMass();// compound nucleus Mass (MeV)
1340  G4double fM2=fM*fM;
1341  G4LorentzVector tg4M=G4LorentzVector(0.,0.,0.,tgM);
1342  G4LorentzVector c4M=tg4M+proj4M;
1343 #ifdef debug
1344  G4cout<<"G4QInelastic::PSDI:fM2="<<fM2<<"<? mc4M="<<c4M.m2()<<",dM="<<fM-tgM<<G4endl;
1345 #endif
1346  if(fM2>=c4M.m2()) // Elastic scattering should be done
1347  {
1348  G4LorentzVector tot4M=tg4M+proj4M+scat4M; // recover the total 4-momentum
1349  s_value=tot4M.m2();
1350  G4double fs=s_value-fM2-ml2;
1351  G4double fMl=fM2*ml2;
1352  G4double hQ2max=(fs*fs/2-fMl-fMl)/s_value; // Maximum possible Q2/2
1353  cost=1.-Q2/hQ2max; // cos(theta) in CMS (use MultProd Q2)
1354 #ifdef debug
1355  G4cout<<"G4QI::PSDI:ct="<<cost<<",Q2="<<Q2<<",hQ2="<<hQ2max<<",4M="<<tot4M<<G4endl;
1356 #endif
1357  G4double acost=std::fabs(cost);
1358  if(acost>1.)
1359  {
1360  if(acost>1.001) G4cout<<"-Warning-G4QInelastic::PostStDoIt: cost="<<cost<<G4endl;
1361  if (cost> 1.) cost= 1.;
1362  else if(cost<-1.) cost=-1.;
1363  }
1364  G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,fM); // 4mom of the recoilNucleus
1365  scat4M=G4LorentzVector(0.,0.,0.,ml); // 4mom of the scatteredLepton
1366  G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-ml)*.01);
1367  if(!G4QHadron(tot4M).RelDecayIn2(scat4M, reco4M, dir4M, cost, cost))
1368  {
1369  G4cerr<<"G4QI::PSDI:t4M="<<tot4M<<",lM="<<ml<<",rM="<<fM<<",cost="<<cost<<G4endl;
1370  //G4Exception("G4QInelastic::PostStepDoIt:","027",FatalException,"ElasticDecay");
1371  }
1372 #ifdef debug
1373  G4cout<<"G4QIn::PStDoI:l4M="<<scat4M<<"+r4M="<<reco4M<<"="<<scat4M+reco4M<<G4endl;
1374 #endif
1375  // ----------------------------------------------------
1376  G4ParticleDefinition* theDefinition=0; // Prototype of a particle for E-Secondaries
1377  // Fill scattered lepton
1378  if (scatPDG==-11) theDefinition = G4Positron::Positron();
1379  else if(scatPDG== 11) theDefinition = G4Electron::Electron();
1380  else if(scatPDG== 13) theDefinition = G4MuonMinus::MuonMinus();
1381  else if(scatPDG==-13) theDefinition = G4MuonPlus::MuonPlus();
1382  //else if(scatPDG== 15) theDefinition = G4TauMinus::TauMinus();
1383  //else if(scatPDG==-15) theDefinition = G4TauPlus::TauPlus();
1384  if (scatPDG==-12) theDefinition = G4AntiNeutrinoE::AntiNeutrinoE();
1385  else if(scatPDG== 12) theDefinition = G4NeutrinoE::NeutrinoE();
1386  else if(scatPDG== 14) theDefinition = G4NeutrinoMu::NeutrinoMu();
1387  else if(scatPDG==-14) theDefinition = G4AntiNeutrinoMu::AntiNeutrinoMu();
1388  //else if(scatPDG== 16) theDefinition = G4NeutrinoTau::NeutrinoTau();
1389  //else if(scatPDG==-16) theDefinition = G4AntiNeutrinoTau::AntiNeutrinoTau();
1390  else G4cout<<"-Warning-G4QInelastic::PostStDoIt: UnknownLepton="<<scatPDG<<G4endl;
1391  G4DynamicParticle* theScL = new G4DynamicParticle(theDefinition,scat4M);
1392  G4Track* scatLep = new G4Track(theScL, localtime, position ); // scattered
1393  scatLep->SetWeight(weight); // weighted
1394  scatLep->SetTouchableHandle(trTouchable); // residual
1395  aParticleChange.AddSecondary(scatLep); // lepton
1396  // Fill residual nucleus
1397  if (fintPDG==90000001) theDefinition = G4Neutron::Neutron(); // neutron
1398  else if(fintPDG==90001000) theDefinition = G4Proton::Proton(); // proton
1399  else // ion
1400  {
1401  G4int fm=static_cast<G4int>(fintPDG/1000000); // Strange part
1402  G4int ZN=fintPDG-1000000*fm;
1403  G4int rZ=static_cast<G4int>(ZN/1000);
1404  G4int rA=ZN-999*rZ;
1405  theDefinition = G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,0,rZ);
1406  }
1407  G4DynamicParticle* theReN = new G4DynamicParticle(theDefinition,reco4M);
1408  G4Track* scatReN = new G4Track(theReN, localtime, position ); // scattered
1409  scatReN->SetWeight(weight); // weighted
1410  scatReN->SetTouchableHandle(trTouchable); // residual
1411  aParticleChange.AddSecondary(scatReN); // nucleus
1412  delete output;
1413  return G4VDiscreteProcess::PostStepDoIt(track, step);
1414  }
1415  } // End of non-quasi-elastic interaction
1416  //
1417  // quasi-elastic (& pickup process) for p+A(Z,N)
1418  //
1419  }
1420  //else if ((projPDG == 2212 || projPDG == 2112) && Z > 0 && N > 0) // Never (in Fragm!)
1421  else if(2>3) // Possibility is closed as quasi-free scattering is made in fragmentation
1422  {
1423  if(momentum<500. && projPDG == 2112) // @@ It's reasonable to add proton capture too !
1424  {
1425  G4LorentzVector tot4M=G4LorentzVector(0.,0.,0.,tgM)+proj4M;
1426  G4double totM2=tot4M.m2();
1427  G4int tZ=Z;
1428  G4int tN=N;
1429  if(projPDG == 2212) tZ++;
1430  else tN++; // @@ Now a neutron is the only alternative
1431  if(momentum<100.) // Try the production threshold (@@ Why 5 MeV?)
1432  {
1433  G4QPDGCode FakeP(2112);
1434  //G4int m1p=tZ-1;
1435  G4int m2n=tN-2;
1436  // @@ For the secondary proton the Coulomb barrier should be taken into account
1437  //G4double mRP= FakeP.GetNuclMass(m1p,tN,0); // Residual to the secondary proton
1438  G4double mR2N= FakeP.GetNuclMass(tZ,m2n,0); // Residual to two secondary neutrons
1439  //G4double mRAl= FakeP.GetNuclMass(tZ-2,m2n,0); // Residual to the secondary alpha
1440  //G4double mR2N= FakeP.GetNuclMass(m1p,tN-1,0); // Residual to the p+n secondaries
1441  G4double minM=mR2N+mNeut+mNeut;
1442  // Compare with other possibilities (sciped for acceleration)
1443  if(totM2 < minM*minM) // DoNothing (under reasonable threshold)
1444  {
1445  aParticleChange.ProposeEnergy(kinEnergy);
1446  aParticleChange.ProposeLocalEnergyDeposit(0.);
1447  aParticleChange.ProposeMomentumDirection(dir);
1448  aParticleChange.ProposeTrackStatus(fAlive);
1449  return G4VDiscreteProcess::PostStepDoIt(track,step);
1450  }
1451  }
1452  }
1453  // Now the quasi-free and PickUp processes should be done:
1454  G4QuasiFreeRatios* qfMan=G4QuasiFreeRatios::GetPointer();
1455  std::pair<G4double,G4double> fief=qfMan->GetRatios(momentum, projPDG, Z, N);
1456  G4double qepart=fief.first*fief.second; // @@ charge exchange is not included @@
1457  // Keep QE part for the quasi-free nucleons
1458  const G4int lCl=3; // The last clProb[lCl]==1. by definition, MUST be increasing
1459  G4double clProb[lCl]={.65,.85,.95};// N/P,D,t/He3,He4, integroProbab for .65,.2,.1,.05
1460  if(qepart>0.45) qepart=.45; // Quasielastic part is too large - shrink
1461  //else qepart/=clProb[0]; // Add corresponding number of 2N, 3N, & 4N clusters
1462  qepart=qepart/clProb[0]-qepart;// Add QE for 2N, 3N, & 4N clusters (N is made in G4QF)
1463  G4double pickup=1.-qepart; // Estimate the rest of the cross-section
1464  G4double thresh=100.;
1465  //if(momentum >thresh) pickup*=50./momentum/std::pow(G4double(Z+N),third);// 50 is Par
1466  if(momentum > thresh) pickup*=50./momentum/G4QThd(Z+N);// 50 is Par
1467  // pickup = 0.; // To exclude the pickup process
1468  if (N) pickup+=qepart;
1469  else pickup =qepart;
1470  G4double rnd=G4UniformRand();
1471 #ifdef ldebug
1472  G4cout<<"-->G4QInelastic::PSD:QE[p("<<proj4M<<")+(Z="<<Z<<",N="<<N<<")]="<<qepart
1473  <<", pickup="<<pickup<<G4endl;
1474 #endif
1475  if(rnd<pickup) // Make a quasi free scattering (out:A-1,h,N) @@ KinLim,CoulBar,PauliBl
1476  {
1477  G4double dmom=91.; // Fermi momentum (proto default for a deuteron)
1478  G4double fmm=287.; // @@ Can be A-dependent
1479  if(Z>1||N>1) dmom=fmm*std::pow(-std::log(G4UniformRand()),third);
1480  // Values to be redefined for quasi-elastic
1481  G4LorentzVector r4M(0.,0.,0.,0.); // Prototype of 4mom of the residual nucleus
1482  G4LorentzVector n4M(0.,0.,0.,0.); // Prototype of 4mom of the quasi-cluster
1483  G4int nPDG=90000001; // Prototype for quasiCluster mass calculation
1484  G4int restPDG=targPDG; // Prototype should be reduced by quasiCluster
1485  G4int rA=Z+N-1; // Prototype for the residualNucl definition
1486  G4int rZ=Z; // residZ: OK for the quasi-free neutron
1487  G4int nA=1; // Prototype for the quasi-cluster definition
1488  G4int nZ=0; // nA=1,nZ=0: OK for the quasi-free neutron
1489  G4double qM=mNeut; // Free mass of the quasi-free cluster
1490  G4int A = Z + N; // Baryon number of the nucleus
1491  if(rnd<qepart)
1492  {
1493  // ===> First decay a nucleus in a nucleon and a residual (A-1) nucleus
1494  // Calculate clusterProbabilities (n,p,d,t,he3,he4 now only, can use UpdateClusters)
1495  G4double base=1.; // Base for randomization (can be reduced by totZ & totN)
1496  G4int max=lCl; // Number of boundaries (can be reduced by totZ & totN)
1497  // Take into account that at least one nucleon must be left !
1498  if(Z<2 || N<2 || A < 6) base = clProb[--max]; // Alpha cluster is impossible
1499  if((Z > 1 && N < 2) || (Z < 2 && N > 1)) base=(clProb[max]+clProb[max-1])/2;// t/He3
1500  if ( (Z < 2 && N < 2) || A < 5) base=clProb[--max];// He3&t clusters are impossible
1501  if(A<3) base=clProb[--max]; // Deuteron cluster is impossible
1502  //G4int cln=0; // Cluster#0 (Default for the selectedNucleon)
1503  G4int cln=1; // Cluster#1 (Default for the selectedDeutron)
1504  //if(max) // Not only nucleons are possible
1505  if(max>1) // Not only deuterons are possible
1506  {
1507  base-=clProb[0]; // Exclude scattering on QF Nucleon
1508  G4double ran=+clProb[0]+base*G4UniformRand(); // Base can be reduced
1509  G4int ic=1; // Start from the smallest cluster boundary
1510  if(max>1) while(ic<max) if(ran>clProb[ic++]) cln=ic;
1511  }
1512  G4ParticleDefinition* theDefinition=0;// Prototype for qfNucleon
1513  G4bool cp1 = cln+2==A; // A=ClusterBN+1 condition
1514  if(!cln || cp1) // Split in nucleon + (A-1) with Fermi momentum
1515  {
1516  G4int nln=0;
1517  if(cln==2) nln=1; // @@ only for cp1: t/He3 choice from A=4
1518  // mass(A)=tM. Calculate masses of A-1 (rM) and mN (mNeut or mProt bounded mass)
1519  if ( ((!cln || cln == 2) && G4UniformRand()*(A-cln) > (N-nln)) ||
1520  ((cln == 3 || cln == 1) && Z > N) )
1521  {
1522  nPDG=90001000; // Update quasi-free nucleon PDGCode to P
1523  nZ=1; // Change charge of the quasiFree nucleon
1524  qM=mProt; // Update quasi-free nucleon mass
1525  rZ--; // Reduce the residual Z
1526  restPDG-=1000; // Reduce the residual PDGCode
1527  }
1528  else restPDG--;
1529  G4LorentzVector t4M(0.,0.,0.,tM); // 4m of the target nucleus to be decayed
1530  G4double rM=G4QPDGCode(restPDG).GetMass();// Mass of the residual nucleus
1531  r4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the residual nucleus
1532  G4double rM2=rM*rM;
1533  G4double nM=std::sqrt(rM2+tM*tM-(tM+tM)*std::sqrt(rM2+dmom*dmom));// M of q-nucleon
1534  n4M=G4LorentzVector(0.,0.,0.,nM); // 4mom of the quasi-nucleon
1535 #ifdef qedebug
1536  G4cout<<"G4QInel::PStDoIt:QE,p="<<dmom<<",tM="<<tM<<",R="<<rM<<",N="<<nM<<G4endl;
1537 #endif
1538  if(!G4QHadron(t4M).DecayIn2(r4M, n4M))
1539  {
1540  //G4cerr<<"G4QInel::PostStDoIt:M="<<tM<<"<r="<<rM<<"+n="<<nM<<"="<<rM+nM<<G4endl;
1541  //G4Exception("G4QInelastic::PostStepDoIt()", "HAD_CHPS_0002",
1542  // FatalException,"Hadronize quasmon: Can't Dec totNuc->QENuc+ResNuc");
1543  G4cerr<<"-W-G4QInel::PoStDoI:M="<<tM<<"<rM="<<rM<<"+nM="<<nM<<"="<<rM+nM<<G4endl;
1544  aParticleChange.ProposeEnergy(kinEnergy);
1545  aParticleChange.ProposeLocalEnergyDeposit(0.);
1546  aParticleChange.ProposeMomentumDirection(dir);
1547  aParticleChange.ProposeTrackStatus(fAlive);
1548  return G4VDiscreteProcess::PostStepDoIt(track,step);
1549  }
1550 #ifdef qedebug
1551  G4cout<<"G4QIn::PStDoIt:QE-N, RA="<<r4M.rho()<<r4M<<",QN="<<n4M.rho()<<n4M<<G4endl;
1552 #endif
1553  if(cp1 && cln) // Quasi-cluster case: swap the output
1554  {
1555  qM=rM; // Scattering will be made on a cluster
1556  nln=nPDG;
1557  nPDG=restPDG;
1558  restPDG=nln;
1559  t4M=n4M;
1560  n4M=r4M;
1561  r4M=t4M;
1562  nln=nZ;
1563  nZ=rZ;
1564  rZ=nln;
1565  nln=nA;
1566  nA=rA;
1567  rA=nln;
1568  }
1569  }
1570  else // Split the cluster (with or w/o "Fermi motion" and "Fermi decay")
1571  {
1572  if(cln==1)
1573  {
1574  nPDG=90001001; // Deuteron
1575  qM=mDeut;
1576  nA=2;
1577  nZ=1;
1578  restPDG-=1001;
1579  // Recalculate dmom
1580  G4ThreeVector sumv=G4ThreeVector(0.,0.,dmom) +
1581  fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection();
1582  dmom=sumv.mag();
1583  }
1584  else if(cln==2)
1585  {
1586  nA=3;
1587  if(G4UniformRand()*(A-2)>(N-1)) // He3
1588  {
1589  nPDG=90002001;
1590  qM=mHel3;
1591  nZ=2;
1592  restPDG-=2001;
1593  }
1594  else // tritium
1595  {
1596  nPDG=90001002;
1597  qM=mTrit;
1598  nZ=1;
1599  restPDG-=1002;
1600  }
1601  // Recalculate dmom
1602  G4ThreeVector sumv=G4ThreeVector(0.,0.,dmom) +
1603  fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection()+
1604  fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection();
1605  dmom=sumv.mag();
1606  }
1607  else
1608  {
1609  nPDG=90002002; // Alpha
1610  qM=mAlph;
1611  nA=4;
1612  nZ=2;
1613  restPDG-=2002;
1614  G4ThreeVector sumv=G4ThreeVector(0.,0.,dmom) +
1615  fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection()+
1616  fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection()+
1617  fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection();
1618  dmom=sumv.mag();
1619  }
1620  rA=A-nA;
1621  rZ=Z-nZ;
1622  // This is a simple case of cluster at rest
1623  //G4double rM=G4QPDGCode(restPDG).GetMass();// Mass of the residual nucleus
1624  //r4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the residual nucleus
1625  //n4M=G4LorentzVector(0.,0.,0.,tM-rM); // 4mom of the quasi-free cluster
1626  // --- End of the "simple case of cluster at rest"
1627  // Make a fake quasi-Fermi distribution for clusters (clusters are not at rest)
1628  G4LorentzVector t4M(0.,0.,0.,tM); // 4m of the target nucleus to be decayed
1629  G4double rM=G4QPDGCode(restPDG).GetMass();// Mass of the residual nucleus
1630  r4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the residual nucleus
1631  G4double rM2=rM*rM;
1632  G4double nM=std::sqrt(rM2+tM*tM-(tM+tM)*std::sqrt(rM2+dmom*dmom));// M of q-cluster
1633  n4M=G4LorentzVector(0.,0.,0.,nM); // 4mom of the quasi-nucleon
1634 #ifdef qedebug
1635  G4cout<<"G4QInel::PStDoIt:QEC, p="<<dmom<<",T="<<tM<<",R="<<rM<<",N="<<nM<<G4endl;
1636 #endif
1637  if(!G4QHadron(t4M).DecayIn2(r4M, n4M))
1638  {
1639  //G4cerr<<"G4QInel::PostStDoIt:M="<<tM<<"<r="<<rM<<"+c="<<nM<<"="<<rM+nM<<G4endl;
1640  //G4Exception("G4QInelastic::PostStepDoIt()", "HAD_CHPS_0003",
1641  // FatalException,"Hadronize quasmon: Can't Dec totNuc->QEClu+ResNuc");
1642  G4cerr<<"-W-G4QInel::PoStDoI:M="<<tM<<"<rM="<<rM<<"+cM="<<nM<<"="<<rM+nM<<G4endl;
1643  aParticleChange.ProposeEnergy(kinEnergy);
1644  aParticleChange.ProposeLocalEnergyDeposit(0.);
1645  aParticleChange.ProposeMomentumDirection(dir);
1646  aParticleChange.ProposeTrackStatus(fAlive);
1647  return G4VDiscreteProcess::PostStepDoIt(track,step);
1648  }
1649  // --- End of the moving cluster implementation ---
1650 #ifdef qedebug
1651  G4cout<<"G4QIn::PStDoIt:QEC, RN="<<r4M.rho()<<r4M<<",QCl="<<n4M.rho()<<n4M<<G4endl;
1652 #endif
1653  }
1654  G4LorentzVector s4M=n4M+proj4M; // Tot 4-momentum for scattering
1655  G4double prjM2 = proj4M.m2(); // Squared mass of the projectile
1656  G4double prjM = std::sqrt(prjM2); // @@ Get from pPDG (?)
1657  G4double minM = prjM+qM; // Min mass sum for the final products
1658  G4double cmM2 =s4M.m2();
1659  if(cmM2>minM*minM)
1660  {
1661 #ifdef qedebug
1662  G4cout<<"G4QInel::PStDoIt:***Enter***,cmM2="<<cmM2<<" > minM2="<<minM*minM<<G4endl;
1663 #endif
1664  // Estimate and randomize charge-exchange with a quasi-free cluster
1665  G4bool chex=false; // Flag of the charge exchange scattering
1666  G4ParticleDefinition* projpt=G4Proton::Proton(); // Prototype, only for chex=true
1667  // This is reserved for the charge-exchange scattering (to help neutron spectras)
1668  //if(cln&&!cp1 &&(projPDG==2212&&rA>rZ || projPDG==2112&&rZ>1))// @@ Use proj chex
1669  if(2>3) // @@ charge exchange is not implemented yet @@
1670  {
1671 #ifdef qedebug
1672  G4cout<<"G4QIn::PStDoIt:-Enter, P="<<projPDG<<",cln="<<cln<<",cp1="<<cp1<<G4endl;
1673 #endif
1674  G4double tprM=mProt;
1675  G4double tprM2=mProt2;
1676  G4int tprPDG=2212;
1677  G4int tresPDG=restPDG+999;
1678  if(projPDG==2212)
1679  {
1680  projpt=G4Neutron::Neutron();
1681  tprM=mNeut;
1682  tprM2=mNeut2;
1683  tprPDG=2112;
1684  tresPDG=restPDG-999;
1685  }
1686  minM=tprM+qM;
1687  G4double efE=(cmM2-tprM2-qM*qM)/(qM+qM);
1688  G4double efP=std::sqrt(efE*efE-tprM2);
1689  G4double chl=qfMan->ChExElCoef(efP*MeV, nZ, nA-nZ, projPDG); // ChEx/Elast(pPDG!)
1690 #ifdef qedebug
1691  G4cout<<"G4QInel::PStDoIt:chl="<<chl<<",P="<<efP<<",nZ="<<nZ<<",nA="<<nA<<G4endl;
1692 #endif
1693  if(chl>0.&&cmM2>minM*minM&&G4UniformRand()<chl/(1.+chl)) // minM is redefined
1694  {
1695  projPDG=tprPDG;
1696  prjM=tprM;
1697  G4double rM=G4QPDGCode(tresPDG).GetMass();// Mass of the residual nucleus
1698  r4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the residual nucleus
1699  n4M=G4LorentzVector(0.,0.,0.,tM-rM); // 4mom of the quasi-free cluster
1700  chex=true; // Confirm charge exchange scattering
1701  }
1702  }
1703  //
1704  std::pair<G4LorentzVector,G4LorentzVector> sctout=qfMan->Scatter(nPDG, n4M,
1705  projPDG, proj4M);
1706 #ifdef qedebug
1707  G4cout<<"G4QInelastic::PStDoIt:QElS,proj="<<prjM<<sctout.second<<",qfCl="<<qM
1708  <<sctout.first<<",chex="<<chex<<",nA="<<nA<<",nZ="<<nZ<<G4endl;
1709 #endif
1710  aParticleChange.ProposeLocalEnergyDeposit(0.); // Everything is in particles
1711  // @@ @@ @@ Coulomb barriers must be checked !! @@ @@ @@ Skip if not
1712  if(chex) // ==> Projectile is changed: fill everything to secondaries
1713  {
1714  aParticleChange.ProposeEnergy(0.); // @@ ??
1715  aParticleChange.ProposeTrackStatus(fStopAndKill); // projectile nucleon is killed
1716  aParticleChange.SetNumberOfSecondaries(3);
1717  G4DynamicParticle* thePrH = new G4DynamicParticle(projpt,sctout.second);
1718  G4Track* scatPrH = new G4Track(thePrH, localtime, position ); // scattered & chex
1719  scatPrH->SetWeight(weight); // weighted
1720  scatPrH->SetTouchableHandle(trTouchable); // projectile
1721  aParticleChange.AddSecondary(scatPrH); // hadron
1722  }
1723  else // ==> The leading particle is filled to the updated projectilee
1724  {
1725  aParticleChange.SetNumberOfSecondaries(2); // @@ if proj=leading
1726  G4double ldT=(sctout.second).e()-prjM; // kin Energy of scat project.
1727  aParticleChange.ProposeEnergy(ldT); // Change the kin Energy
1728  G4ThreeVector ldV=(sctout.second).vect(); // Change momentum direction
1729  aParticleChange.ProposeMomentumDirection(ldV/ldV.mag());
1730  aParticleChange.ProposeTrackStatus(fAlive);
1731  }
1732  // ---------------------------------------------------------
1733  // Fill scattered quasi-free nucleon/fragment
1734  if (nPDG==90000001) theDefinition = G4Neutron::Neutron();
1735  else if(nPDG==90001000) theDefinition = G4Proton::Proton();
1736  else if(nZ>0 && nA>1)
1737  theDefinition = G4ParticleTable::GetParticleTable()->FindIon(nZ,nA,0,nZ);
1738 #ifdef debug
1739  else G4cout<<"-Warning_G4QIn::PSD:scatqfPDG="<<nPDG<<", Z="<<nZ<<",A="<<nA<<G4endl;
1740 #endif
1741  if(nZ>0 && nA>0)
1742  {
1743  G4DynamicParticle* theQFN = new G4DynamicParticle(theDefinition,sctout.first);
1744  G4Track* scatQFN = new G4Track(theQFN, localtime, position ); // scattered
1745  scatQFN->SetWeight(weight); // weighted
1746  scatQFN->SetTouchableHandle(trTouchable); // quasi-free
1747  aParticleChange.AddSecondary(scatQFN); // nucleon/cluster
1748  }
1749  // ----------------------------------------------------
1750  // Fill residual nucleus
1751  if (restPDG==90000001) theDefinition = G4Neutron::Neutron();
1752  else if(restPDG==90001000) theDefinition = G4Proton::Proton();
1753  else if(rZ>0 && rA>1)
1754  theDefinition = G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,0,rZ);
1755 #ifdef debug
1756  else G4cout<<"-Warning_G4QIn::PSD: resPDG="<<restPDG<<",Z="<<rZ<<",A="<<rA<<G4endl;
1757 #endif
1758  if(rZ>0 && rA>0)
1759  {
1760  G4DynamicParticle* theReN = new G4DynamicParticle(theDefinition,r4M);
1761  G4Track* scatReN = new G4Track(theReN, localtime, position ); // scattered
1762  scatReN->SetWeight(weight); // weighted
1763  scatReN->SetTouchableHandle(trTouchable); // residual
1764  aParticleChange.AddSecondary(scatReN); // nucleus
1765  }
1766  delete output;
1767  return G4VDiscreteProcess::PostStepDoIt(track, step);
1768  }
1769 #ifdef qedebug
1770  else G4cout<<"G4QInel::PSD:OUT,M2="<<s4M.m2()<<"<"<<minM*minM<<", N="<<nPDG<<G4endl;
1771 #endif
1772  } // end of proton quasi-elastic (including QE on NF)
1773  else // @@ make cost-condition @@ Pickup process (pickup 1 or 2 n and make d or t)
1774  {
1775  if(projPDG==2212) restPDG--; // Res. nucleus PDG is one neutron less than targ. PDG
1776  else
1777  {
1778  restPDG-=1000; // Res. nucleus PDG is one proton less than targ. PDG
1779  rZ--; // Reduce the charge of the residual nucleus
1780  }
1781  G4double mPUF=mDeut; // Prototype of the mass of the created pickup fragment
1782  G4ParticleDefinition* theDefinition = G4Deuteron::Deuteron(); // Default: make d
1783  //theDefinition = G4ParticleTable::GetParticleTable()->FindIon(nZ,nA,0,nZ);//ion
1784  G4bool din=false; // Flag of picking up 2 neutrons to create t (tritium)(p)
1785  G4bool dip=false; // Flag of picking up 2 protons to create t (tritium) (n)
1786  G4bool pin=false; // Flag of picking up 1 n + 1 p to create He3/t (p/n)
1787  G4bool hin=false; // Flag of PickUp creation of alpha (He4) (p/n)
1788  G4double frn=G4UniformRand();
1789  if(N>2 && frn > 0.95) // .95 is a parameter to pickup 2N (to cteate t or He3)
1790  {
1791  if(projPDG==2212)
1792  {
1793  if(N>1 && G4UniformRand()*(Z+.5*(N-1)) > Z)
1794  {
1795  theDefinition = G4Triton::Triton();
1796  mPUF=mTrit; // The mass of the created pickup fragment is triton
1797  restPDG--; // Res. nucleus PDG is two neutrons less than target PDG
1798  din=true;
1799  }
1800  else
1801  {
1802  theDefinition = G4He3::He3();
1803  mPUF=mHel3; // The mass of the created pickup fragment is Helium3
1804  restPDG-=1000; // Res. nucleus PDG is two neutrons less than target PDG
1805  rZ--; // Reduce the charge of the residual nucleus
1806  pin=true;
1807  }
1808  }
1809  else // Neutron projectile (2112)
1810  {
1811  if(Z>1 && G4UniformRand()*(N+.5*(Z-1)) > N)
1812  {
1813  theDefinition = G4He3::He3();
1814  mPUF=mHel3; // The mass of the created pickup fragment is Helium3
1815  restPDG-=1000; // Res. nucleus PDG is two neutrons less than target PDG
1816  rZ--; // Reduce the charge of the residual nucleus
1817  dip=true;
1818  }
1819  else
1820  {
1821  theDefinition = G4Triton::Triton();
1822  mPUF=mTrit; // The mass of the created pickup fragment is triton
1823  restPDG--; // Res. nucleus PDG is two neutrons less than target PDG
1824  pin=true;
1825  }
1826  }
1827  rA--; // Reduce the baryon number of the residual nucleus
1828  // Recalculate dmom
1829  G4ThreeVector sumv=G4ThreeVector(0.,0.,dmom) +
1830  fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection();
1831  dmom=sumv.mag();
1832  if(Z>1 && frn > 0.99) // .99 is a parameter to pickup 2n1p || 1n2p & create alpha
1833  {
1834  theDefinition = G4Alpha::Alpha();
1835  if((din && projPDG==2212) || (pin && projPDG==2112))
1836  {
1837  restPDG-=1000; // Residual nucleus PDG is reduced to create alpha
1838  rZ--; // Reduce the charge of the residual nucleus
1839  }
1840  else if((pin && projPDG==2212) || (dip && projPDG==2112)) restPDG--;
1841  else G4cout<<"-Warning-G4QIn::PSD: PickUp logic error, proj="<<projPDG<<G4endl;
1842  hin=true;
1843  mPUF=mAlph; // The mass of the created pickup fragment (alpha)
1844  rA--; // Reduce the baryon number of the residual nucleus
1845  // Recalculate dmom
1846  sumv += fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection();
1847  dmom=sumv.mag();
1848  }
1849  }
1850  G4double mPUF2=mPUF*mPUF;
1851  G4LorentzVector t4M(0.,0.,0.,tM); // 4m of the target nucleus to be decayed
1852  G4double rM=G4QPDGCode(restPDG).GetMass();// Mass of the residual nucleus
1853  G4double rM2=rM*rM; // Squared mass of the residual nucleus
1854  G4double nM2=rM2+tM*tM-(tM+tM)*std::sqrt(rM2+dmom*dmom);// M2 of boundQF-nucleon(2N)
1855  if(nM2 < 0) nM2=0.;
1856  G4double den2=(dmom*dmom+nM2); // squared energy of bQF-nucleon
1857  G4double den=std::sqrt(den2); // energy of bQF-nucleon
1858 #ifdef qedebug
1859  G4double nM=std::sqrt(nM2); // M of bQF-nucleon
1860  G4cout<<"G4QInel::PStDoIt:PiU, p="<<dmom<<",tM="<<tM<<", R="<<rM<<",N="<<nM<<G4endl;
1861 #endif
1862  G4double qp=momentum*dmom;
1863  G4double srm=proj4M.e()*den;
1864  G4double mNucl2=mProt2;
1865  if(projPDG == 2112) mNucl2=mNeut2;
1866  G4double cost=(nM2+mNucl2+srm+srm-mPUF2)/(qp+qp);
1867  G4int ict=0;
1868  while(std::fabs(cost)>1. && ict<3)
1869  {
1870  dmom=91.; // Fermi momentum (proDefault for a deuteron)
1871  if(Z>1||N>1) dmom=fmm*std::pow(-std::log(G4UniformRand()),third);
1872  if(din || pin || dip) // Recalculate dmom for the final t/He3
1873  {
1874  G4ThreeVector sumv=G4ThreeVector(0.,0.,dmom) +
1875  fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection();
1876  if(hin)
1877  sumv+=fmm*std::pow(-std::log(G4UniformRand()),third)*G4RandomDirection();
1878  dmom=sumv.mag();
1879  }
1880  nM2=rM2+tM*tM-(tM+tM)*std::sqrt(rM2+dmom*dmom); // M2 of bQF-nucleon
1881  if(nM2 < 0) nM2=0.;
1882  //nM=std::sqrt(nM2); // M of bQF-nucleon --Not used?--
1883  den2=(dmom*dmom+nM2); // squared energy of bQF-nucleon
1884  den=std::sqrt(den2); // energy of bQF-nucleon
1885  qp=momentum*dmom;
1886  srm=proj4M.e()*den;
1887  cost=(nM2+mNucl2+srm+srm-mPUF2)/(qp+qp);
1888  ict++;
1889 #ifdef ldebug
1890  if(ict>2)G4cout<<"G4QInelast::PStDoIt:i="<<ict<<",d="<<dmom<<",ct="<<cost<<G4endl;
1891 #endif
1892  }
1893  if(std::fabs(cost)<=1.)
1894  {// ---- @@ From here can be a MF for QF nucleon extraction (if used by others)
1895  G4ThreeVector vdir = proj4M.vect(); // 3-Vector in the projectile direction
1896  G4ThreeVector vx(0.,0.,1.); // ProtoOrt in theDirection of the projectile
1897  G4ThreeVector vy(0.,1.,0.); // First ProtoOrt orthogonal to the direction
1898  G4ThreeVector vz(1.,0.,0.); // SecondProtoOrt orthoganal to the direction
1899  if(vdir.mag2() > 0.) // the projectile isn't at rest
1900  {
1901  vx = vdir.unit(); // Ort in the direction of the projectile
1902  G4ThreeVector vv= vx.orthogonal(); // Not normed orthogonal vector (!)
1903  vy = vv.unit(); // First ort orthogonal to the proj direction
1904  vz = vx.cross(vy); // Second ort orthoganal to the projDirection
1905  }
1906  // ---- @@ End of possible MF (Similar is in G4QEnvironment)
1907  G4double tdom=dmom*std::sqrt(1-cost*cost);// Transversal component of the Fermi-Mom
1908  G4double phi=twopi*G4UniformRand(); // Phi of the Fermi-Mom
1909  G4ThreeVector vedm=vx*dmom*cost+vy*tdom*std::sin(phi)+vz*tdom*std::cos(phi);// bQFN
1910  G4LorentzVector bqf(vedm,den); // 4-mom of the bQF nucleon
1911  r4M=t4M-bqf; // 4mom of the residual nucleus
1912 #ifdef debug
1913  if(std::fabs(r4M.m()-rM)>.001) G4cout<<"G4QIn::PSD: rM="<<rM<<"#"<<r4M.m()<<G4endl;
1914 #endif
1915  G4LorentzVector f4M=proj4M+bqf; // Prototype of 4-mom of Deuteron
1916 #ifdef debug
1917  if(std::fabs(f4M.m()-mPUF)>.001)G4cout<<"G4QI::PSD:m="<<mPUF<<"#"<<f4M.m()<<G4endl;
1918 #endif
1919  aParticleChange.ProposeEnergy(0.); // @@ ??
1920  aParticleChange.ProposeTrackStatus(fStopAndKill);// projectile nucleon is killed
1921  aParticleChange.SetNumberOfSecondaries(2);
1922  // Fill pickup hadron
1923  G4DynamicParticle* theQFN = new G4DynamicParticle(theDefinition,f4M);
1924  G4Track* scatQFN = new G4Track(theQFN, localtime, position ); // pickup
1925  scatQFN->SetWeight(weight); // weighted
1926  scatQFN->SetTouchableHandle(trTouchable); // nuclear
1927  aParticleChange.AddSecondary(scatQFN); // cluster
1928 #ifdef pickupd
1929  G4cout<<"G4QInelastic::PStDoIt: f="<<theDefinition<<",f4M="<<f4M.m()<<f4M<<G4endl;
1930 #endif
1931  // ----------------------------------------------------
1932  // Fill residual nucleus
1933  if (restPDG==90000001) theDefinition = G4Neutron::Neutron();
1934  else if(restPDG==90001000) theDefinition = G4Proton::Proton();
1935  else theDefinition = G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,0,rZ);//ion
1936  G4DynamicParticle* theReN = new G4DynamicParticle(theDefinition,r4M);
1937  G4Track* scatReN = new G4Track(theReN, localtime, position ); // scattered
1938  scatReN->SetWeight(weight); // weighted
1939  scatReN->SetTouchableHandle(trTouchable); // residual
1940  aParticleChange.AddSecondary(scatReN); // nucleus
1941 #ifdef pickupd
1942  G4cout<<"G4QInelastic::PStDoIt:rZ="<<rZ<<",rA="<<rA<<",r4M="<<r4M.m()<<r4M<<G4endl;
1943 #endif
1944  delete output;
1945  return G4VDiscreteProcess::PostStepDoIt(track, step);
1946  }
1947  }
1948  } // end of the quasi-elastic decoupled process for the neucleon projectile
1949  } // end lepto-nuclear, neutrino-nuclear, proton/neutron decoupled processes
1950  EnMomConservation=proj4M+G4LorentzVector(0.,0.,0.,tM); // Total 4-mom of the reaction
1951  if(absMom) EnMomConservation+=lead4M; // Add E/M of leading System
1952 #ifdef debug
1953  G4cout<<"G4QInelastic::PostStDoIt:before St="<<aParticleChange.GetTrackStatus()<<G4endl;
1954 #endif
1955 
1956  // @@@@@@@@@@@@@@ Temporary for the testing purposes --- Begin
1957  //G4bool elF=false; // Flag of the ellastic scattering is "false" by default
1958  //G4double eWei=1.;
1959  // @@@@@@@@@@@@@@ Temporary for the testing purposes --- End
1960 #ifdef debug
1961  G4cout<<"^^G4QInelastic::PostStepDoIt: projPDG="<<projPDG<<", targPDG="<<targPDG<<G4endl;
1962 #endif
1963  const G4QNucleus targNuc(Z,N); // Define the target nucleus
1964  if(projPDG>9000000)
1965  {
1966  delete output; // Before the output creation
1967  G4QNucleus projNuc(proj4M,projPDG); // Define the projectile nucleus
1968  G4QIonIonCollision IonIon(projNuc, targNuc); // Define deep-inelastic ion-ion reaction
1969 #ifdef debug
1970  G4cout<<"G4QInel::PStDoIt: ProjNuc="<<projPDG<<proj4M<<", TargNuc="<<targPDG<<G4endl;
1971 #endif
1972  try {output = IonIon.Fragment();} // DINR AA interaction products
1973  catch (G4QException& error)
1974  {
1975  G4cerr<<"***G4QInelastic::PostStepDoIt: G4QE Exception is catched in hA"<<G4endl;
1976  // G4Exception("G4QInelastic::PostStepDoIt:","27",FatalException,"CHIPS hA crash");
1977  G4Exception("G4QInelastic::PostStepDoIt()", "HAD_CHPS_0027",
1978  FatalException, "CHIPS hA crash");
1979  }
1980  }
1981  else // --> A projectile hadron
1982  {
1983  // *** Let to produce elastic final states for acceleration ***
1984 #ifdef debug
1985  G4int maxCn=7;
1986 #endif
1987  G4int atCn=0; // Attempts counter
1988  //G4bool inel=false;
1989  //while (inel==false && ++atCn <= maxCn) // To evoid elastic final state
1990  //{
1991 #ifdef debug
1992  G4cout<<"G4QInelastic::PStDoIt: atCn="<<atCn<<" <= maxCn="<<maxCn<<G4endl;
1993 #endif
1994  G4int outN=output->size();
1995  if(outN) // The output is not empty
1996  {
1997  std::for_each(output->begin(), output->end(), DeleteQHadron());
1998  output->clear();
1999  }
2000  delete output; // Before the new output creation
2001  const G4QHadron projH(projPDG,proj4M); // Define the projectile particle
2002 #ifdef debug
2003  G4cout<<"G4QInelastic::PStDoIt: Before Fragmentation"<<G4endl;
2004 #endif
2005  //if(aProjPDG==22 && projPDG==22 && proj4M.m2()<.01 && proj4M.e()<150.)
2006  //{
2007  // G4QHadron* tgH= new G4QHadron(targNuc.GetPDGCode(),targNuc.Get4Momentum());
2008  // G4QHadron* prH= new G4QHadron(projPDG,proj4M);
2009  // output = new G4QHadronVector();
2010  // output->push_back(prH);
2011  // output->push_back(tgH);
2012  //}
2013  //else
2014  //{
2015  G4QFragmentation DINR(targNuc, projH); // Define deep-inel. h-A reaction
2016 #ifdef debug
2017  G4cout<<"G4QInelastic::PStDoIt:Proj="<<projPDG<<proj4M<<",TgNuc="<<targPDG<<G4endl;
2018 #endif
2019  try {output = DINR.Fragment();} // DINR hA fragmentation
2020  catch (G4QException& error)
2021  {
2022  G4cerr<<"***G4QInelastic::PostStepDoIt: G4QE Exception is catched in hA"<<G4endl;
2023  // G4Exception("G4QInelastic::PostStepDoIt:","27",FatalException,"CHIPS hA crash");
2024  G4Exception("G4QInelastic::PostStepDoIt()", "HAD_CHPS_0027",
2025  FatalException, "CHIPS hA crash");
2026  }
2027  //}
2028  outN=output->size();
2029 #ifdef debug
2030  G4cout<<"G4QInelastic::PostStepDoIt: At# "<<atCn<<", nSec="<<outN<<", fPDG="
2031  <<(*output)[0]->GetPDGCode()<<", pPDG="<< projPDG<<G4endl;
2032 #endif
2033  //inel=true; // For while
2034  if(outN < 2)
2035  {
2036  G4cout<<"-Warning-G4QInelastic::PostStepDoIt: nSec="<<outN<<", At# "<<atCn<<G4endl;
2037  //inel=false; // For while
2038  }
2039  //else if(outN==2 && (*output)[0]->GetPDGCode() == projPDG) inel=false; // For while
2040 #ifdef debug
2041  if(atCn==maxCn)G4cout<<"-Warning-G4QI::PostStDoIt:mAt="<<atCn<<" is reached"<<G4endl;
2042 #endif
2043  //}
2044  }
2045  //
2046  // --- the scattered hadron with changed nature can be added here ---
2047  if(scat)
2048  {
2049  G4QHadron* scatHadron = new G4QHadron(scatPDG,scat4M);
2050  output->push_back(scatHadron);
2051  }
2052  G4int qNH=0;
2053  if(leadhs) qNH=leadhs->size();
2054  if(absMom)
2055  {
2056  if(qNH) for(G4int iq=0; iq<qNH; iq++)
2057  {
2058  G4QHadron* loh=(*leadhs)[iq]; // Pointer to the output hadron
2059  output->push_back(loh);
2060  }
2061  if(leadhs) delete leadhs;
2062  leadhs=0;
2063  }
2064  else
2065  {
2066  if(qNH) for(G4int iq=0; iq<qNH; iq++) delete (*leadhs)[iq];
2067  if(leadhs) delete leadhs;
2068  leadhs=0;
2069  }
2070  // ------------- From here the secondaries are filled -------------------------
2071  G4int tNH = output->size(); // A#of hadrons in the output
2072  aParticleChange.SetNumberOfSecondaries(tNH); // @@ tNH can be changed (drop/decay!)
2073  // Now add nuclear fragments
2074 #ifdef debug
2075  G4cout<<"G4QInelastic::PostStepDoIt: "<<tNH<<" particles are generated"<<G4endl;
2076 #endif
2077 #ifdef ppdebug
2078  if(absMom)G4cout<<"G4QInelastic::PostStepDoIt: t="<<tNH<<", q="<<qNH<<G4endl;
2079 #endif
2080  G4int nOut=output->size(); // Real length of the output @@ Temporary
2081  if(tNH==1 && !scat) // @@ Temporary. Find out why it happened!
2082  {
2083  G4cout<<"-Warning-G4QInelastic::PostStepDoIt: 1 secondary! absMom="<<absMom;
2084  if(absMom) G4cout<<", qNH="<<qNH;
2085  G4cout<<", PDG0="<<(*output)[0]->GetPDGCode();
2086  G4cout<<G4endl;
2087  tNH=0;
2088  delete output->operator[](0); // delete the creazy hadron
2089  output->pop_back(); // clean up the output vector
2090  }
2091  if(tNH==2&&2!=nOut) G4cout<<"--Warning--G4QInelastic::PostStepDoIt: 2 # "<<nOut<<G4endl;
2092  // Deal with ParticleChange final state interface to GEANT4 output of the process
2093  //if(tNH==2) for(i=0; i<tNH; i++) // @@ Temporary tNH==2 instead of just tNH
2094  if(tNH) for(i=0; i<tNH; i++) // @@ Temporary tNH==2 instead of just tNH
2095  {
2096  // Note that one still has to take care of Hypernuclei (with Lambda or Sigma inside)
2097  // Hypernucleus mass calculation and ion-table interface upgrade => work for Hisaya @@
2098  // The decau process for hypernuclei must be developed in GEANT4 (change CHIPS body)
2099  G4QHadron* hadr=(*output)[i]; // Pointer to the output hadron
2100  G4int PDGCode = hadr->GetPDGCode();
2101  G4int nFrag = hadr->GetNFragments();
2102 #ifdef pdebug
2103  G4cout<<"G4QInelastic::PostStepDoIt: H#"<<i<<",PDG="<<PDGCode<<",nF="<<nFrag
2104  <<", 4Mom="<<hadr->Get4Momentum()<<G4endl;
2105 #endif
2106  if(nFrag) // Skip intermediate (decayed) hadrons
2107  {
2108 #ifdef debug
2109  G4cout<<"G4QInelastic::PostStepDoIt: Intermediate particle is found i="<<i<<G4endl;
2110 #endif
2111  delete hadr;
2112  continue;
2113  }
2114  G4DynamicParticle* theSec = new G4DynamicParticle;
2115  G4ParticleDefinition* theDefinition;
2116  if (PDGCode==90000001) theDefinition = G4Neutron::Neutron();
2117  else if(PDGCode==90001000) theDefinition = G4Proton::Proton();//While it can be in ions
2118  else if(PDGCode==91000000) theDefinition = G4Lambda::Lambda();
2119  else if(PDGCode==311 || PDGCode==-311)
2120  {
2121  if(G4UniformRand()>.5) theDefinition = G4KaonZeroLong::KaonZeroLong(); // K_L
2122  else theDefinition = G4KaonZeroShort::KaonZeroShort(); // K_S
2123  }
2124  else if(PDGCode==91000999) theDefinition = G4SigmaPlus::SigmaPlus();
2125  else if(PDGCode==90999001) theDefinition = G4SigmaMinus::SigmaMinus();
2126  else if(PDGCode==91999000) theDefinition = G4XiMinus::XiMinus();
2127  else if(PDGCode==91999999) theDefinition = G4XiZero::XiZero();
2128  else if(PDGCode==92998999) theDefinition = G4OmegaMinus::OmegaMinus();
2129  else if(PDGCode==90000000)
2130  {
2131 #ifdef pdebug
2132  G4cout<<"-Warning-G4QInelastic::PostStepDoIt:Vacuum PDG="<<PDGCode
2133  <<", 4Mom="<<hadr->Get4Momentum()<<G4endl;
2134 #endif
2135  theDefinition = G4Gamma::Gamma();
2136  G4double hM2=(hadr->Get4Momentum()).m2();
2137  if(std::fabs(hM2)>.01) G4cout<<"-Warning-G4QInel::PoStDoIt:90000000M2="<<hM2<<G4endl;
2138  else if(hadr->Get4Momentum()==vacuum4M)
2139  {
2140  delete hadr;
2141  continue;
2142  }
2143  }
2144  else if(PDGCode >80000000) // Defines hypernuclei as normal nuclei (N=N+S Correction!)
2145  {
2146  G4int aZ = hadr->GetCharge();
2147  G4int aA = hadr->GetBaryonNumber();
2148 #ifdef pdebug
2149  G4cout<<"G4QInelastic::PostStepDoIt:Ion Z="<<aZ<<", A="<<aA<<G4endl;
2150 #endif
2151  theDefinition = G4ParticleTable::GetParticleTable()->FindIon(aZ,aA,0,aZ);
2152  }
2153  //else theDefinition = G4ParticleTable::GetParticleTable()->FindParticle(PDGCode);
2154  else if(PDGCode==89999998 || PDGCode==89998000 || PDGCode==88000000) // anti-di-baryon
2155  {
2156  G4double rM=mNeut; // Prototype of the baryon mass
2157  G4int rPDG=-2112; // Prototype of the baryon PDG
2158  G4ParticleDefinition* BarDef=G4Neutron::Neutron();// Prototype of theBaryonDefinition
2159  if (PDGCode==89998000)
2160  {
2161  rM=mProt;
2162  rPDG=-2212;
2163  BarDef=G4Proton::Proton();
2164  }
2165  else if(PDGCode==88000000)
2166  {
2167  rM=mLamb;
2168  rPDG=-3122;
2169  BarDef=G4Lambda::Lambda();
2170  }
2171  G4LorentzVector t4M=hadr->Get4Momentum(); // 4m of the di-baryon
2172  G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the 1st anti-baryon
2173  G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the 2nd anti-baryon
2174 #ifdef qedebug
2175  G4cout<<"G4QInel::PStDoIt:AntiDiBar,t4M="<<tM<<",m="<<rM<<",PDG="<<PDGCode<<G4endl;
2176 #endif
2177  if(!G4QHadron(t4M).DecayIn2(f4M, s4M))
2178  {
2179  G4cerr<<"G4QIn::PostStDoIt: ADB, M="<<t4M.m()<<" < 2*rM="<<rM<<" = "<<2*rM<<G4endl;
2180  // throw G4QException("G4QInelastic::HadronizeQuasm:Can't decay anti-dibaryon");
2181  G4Exception("G4QInelastic::PostStepDoIt()", "HAD_CHPS_0004",
2182  FatalException, "Hadronize quasmon: Can't decay anti-dibaryon");
2183  }
2184  // --- End of the moving cluster implementation ---
2185 #ifdef qedebug
2186  G4cout<<"G4QInelastic::PStDoIt: ADB, H1="<<rM<<f4M<<", H2="<<rM<<s4M<<G4endl;
2187 #endif
2188  G4QHadron fH(rPDG,f4M);
2189  hadr->Set4Momentum(f4M);
2190  hadr->SetQPDG(fH.GetQPDG());
2191  theDefinition=BarDef;
2192 #ifdef debug
2193  G4cout<<"G4QInel::PostStDoIt:Anti-DiBar, DecayIn2, h1="<<rPDG<<f4M<<G4endl;
2194 #endif
2195  G4QHadron* sH = new G4QHadron(rPDG,s4M);
2196  output->push_back(sH);
2197  ++tNH;
2198 #ifdef debug
2199  G4cout<<"G4QInel::PostStDoIt:Anti-DiBar, DecayIn2, h2="<<rPDG<<s4M<<G4endl;
2200 #endif
2201  }
2202  else if(PDGCode==90000997 || PDGCode==89997001) // anti-NDelta
2203  {
2204  G4double rM=mNeut; // Prototype of the baryon mass
2205  G4int rPDG=-2112; // Prototype of the baryon PDG
2206  G4double iM=mPi; // Prototype of the pion mass
2207  G4int iPDG= 211; // Prototype of the pion PDG
2208  G4ParticleDefinition* BarDef=G4Neutron::Neutron();// Prototype of theBaryonDefinition
2209  if(PDGCode==90000997)
2210  {
2211  rM=mProt;
2212  rPDG=-2212;
2213  iPDG=-211;
2214  BarDef=G4Proton::Proton();
2215  }
2216  G4LorentzVector t4M=hadr->Get4Momentum(); // 4m of the di-baryon
2217  G4LorentzVector f4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the 1st anti-baryon
2218  G4LorentzVector s4M=G4LorentzVector(0.,0.,0.,rM); // 4mom of the 2nd anti-baryon
2219  G4LorentzVector u4M=G4LorentzVector(0.,0.,0.,iM); // 4mom of the pion
2220 #ifdef qedebug
2221  G4cout<<"G4QInel::PStDoIt:AntiNDelta, t4M="<<tM<<",m="<<rM<<",PDG="<<PDGCode<<G4endl;
2222 #endif
2223  if(!G4QHadron(t4M).DecayIn3(f4M, s4M, u4M))
2224  {
2225  G4cerr<<"G4QIn::PostStDoIt: AND, tM="<<t4M.m()<<" < 2*mB+mPi="<<2*rM+iM<<G4endl;
2226  // throw G4QException("G4QInelastic::HadronizeQuasm:Can't decay anti-NDelta");
2227  G4Exception("G4QInelastic::PostStepDoIt()", "HAD_CHPS_0005",
2228  FatalException, "Hadronize quasmon: Can't decay anti-NDelta");
2229  }
2230  // --- End of the moving cluster implementation ---
2231 #ifdef qedebug
2232  G4cout<<"G4QInel::PStDoI:AND,B1="<<rM<<f4M<<",B2="<<rM<<s4M<<",Pi="<<iM<<u4M<<G4endl;
2233 #endif
2234  G4QHadron fH(rPDG,f4M);
2235  hadr->Set4Momentum(f4M);
2236  hadr->SetQPDG(fH.GetQPDG());
2237  theDefinition=BarDef;
2238 #ifdef debug
2239  G4cout<<"G4QInel::PostStDoIt:Anti-NDelta, DecayIn2, h1="<<rPDG<<f4M<<G4endl;
2240 #endif
2241  G4QHadron* sH = new G4QHadron(rPDG,s4M);
2242  output->push_back(sH);
2243  ++tNH;
2244 #ifdef debug
2245  G4cout<<"G4QInel::PostStDoIt:Anti-NDelta, DecayIn2, h2="<<rPDG<<s4M<<G4endl;
2246 #endif
2247  G4QHadron* uH = new G4QHadron(iPDG,u4M);
2248  output->push_back(uH);
2249  ++tNH;
2250 #ifdef debug
2251  G4cout<<"G4QInel::PostStDoIt:Anti-NDelta, DecayIn2, h2="<<rPDG<<s4M<<G4endl;
2252 #endif
2253  }
2254  else
2255  {
2256 #ifdef pdebug
2257  G4cout<<"G4QInelastic::PostStepDoIt:Define particle with PDG="<<PDGCode<<G4endl;
2258 #endif
2259  theDefinition = G4QPDGToG4Particle::Get()->GetParticleDefinition(PDGCode);
2260 #ifdef pdebug
2261  G4cout<<"G4QInelastic::PostStepDoIt:AfterParticleDefinition PDG="<<PDGCode<<G4endl;
2262 #endif
2263  }
2264  if(!theDefinition)
2265  { // !! Commenting of this print costs days of searching for E/mom nonconservation !!
2266  if(PDGCode!=90000000 || hadr->Get4Momentum()!=vacuum4M)
2267  G4cout<<"---Warning---G4QInelastic::PostStepDoIt: drop PDG="<<PDGCode<<", 4Mom="
2268  <<hadr->Get4Momentum()<<G4endl;
2269  delete hadr;
2270  continue;
2271  }
2272 #ifdef pdebug
2273  G4cout<<"G4QInelastic::PostStepDoIt:Name="<<theDefinition->GetParticleName()<<G4endl;
2274 #endif
2275  theSec->SetDefinition(theDefinition);
2276  G4LorentzVector h4M=hadr->Get4Momentum();
2277  EnMomConservation-=h4M;
2278 #ifdef tdebug
2279  G4cout<<"G4QInelast::PSDI: "<<i<<","<<PDGCode<<h4M<<h4M.m()<<EnMomConservation<<G4endl;
2280 #endif
2281 #ifdef debug
2282  G4cout<<"G4QInelastic::PostStepDoIt:#"<<i<<",PDG="<<PDGCode<<",4M="<<h4M<<G4endl;
2283 #endif
2284  theSec->Set4Momentum(h4M); // ^
2285  delete hadr; // <-----<-----------<-------------<---------------------<---------<-----+
2286 #ifdef debug
2287  G4ThreeVector curD=theSec->GetMomentumDirection(); // ^
2288  G4double curM=theSec->GetMass(); // |
2289  G4double curE=theSec->GetKineticEnergy()+curM; // ^
2290  G4cout<<"G4QInelast::PSDoIt:p="<<curD<<curD.mag()<<",e="<<curE<<",m="<<curM<<G4endl;//|
2291 #endif
2292  G4Track* aNewTrack = new G4Track(theSec, localtime, position ); // ^
2293  aNewTrack->SetWeight(weight); // weighted |
2294  aNewTrack->SetTouchableHandle(trTouchable); // |
2295  aParticleChange.AddSecondary( aNewTrack ); // |
2296 #ifdef debug
2297  G4cout<<"G4QInelastic::PostStepDoIt:#"<<i<<" is done."<<G4endl; // |
2298 #endif
2299  } // |
2300  delete output; // instances of the G4QHadrons from the output are already deleted above +
2301  if(leadhs) // To satisfy Valgrind ( How can that be?)
2302  {
2303  qNH=leadhs->size();
2304  if(qNH) for(G4int iq=0; iq<qNH; iq++) delete (*leadhs)[iq];
2305  delete leadhs;
2306  leadhs=0;
2307  }
2308 #ifdef debug
2309  G4cout<<"G4QInelastic::PostStDoIt: after St="<<aParticleChange.GetTrackStatus()<<G4endl;
2310 #endif
2311  if(aProjPDG!=11 && aProjPDG!=13 && aProjPDG!=15)
2312  aParticleChange.ProposeTrackStatus(fStopAndKill); // Kill the absorbed particle
2313 #ifdef pdebug
2314  G4cout<<"G4QInelastic::PostStepDoIt: E="<<aParticleChange.GetEnergy()
2315  <<", d="<<*aParticleChange.GetMomentumDirection()<<G4endl;
2316 #endif
2317 #ifdef ldebug
2318  G4cout<<"G4QInelastic::PostStepDoIt:*** PostStepDoIt is done ***, P="<<aProjPDG<<", St="
2319  <<aParticleChange.GetTrackStatus()<<G4endl;
2320 #endif
2321  return G4VDiscreteProcess::PostStepDoIt(track, step);
2322 }
2323 
2324 std::pair<G4double,G4double> G4QInelastic::Random2DDirection()
2325 {
2326  G4double sp=0; // sin(phi)
2327  G4double cp=1.; // cos(phi)
2328  G4double r2=2.; // to enter the loop
2329  while(r2>1. || r2<.0001) // pi/4 efficiency
2330  {
2331  G4double sine=G4UniformRand();
2332  G4double cosine=G4UniformRand();
2333  sp=1.-sine-sine;
2334  cp=1.-cosine-cosine;
2335  r2=sp*sp+cp*cp;
2336  }
2337  G4double norm=std::sqrt(r2);
2338  return std::make_pair(sp/norm,cp/norm);
2339 }
Generated on Sat May 25 2013 14:33:50 for Geant4 by   doxygen 1.8.4