source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/interface/src/G4QCollision.cc @ 836

Last change on this file since 836 was 819, checked in by garnier, 16 years ago

import all except CVS

File size: 82.6 KB
Line 
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: G4QCollision.cc,v 1.24.2.2 2008/05/07 13:47:31 gcosmo Exp $
27// GEANT4 tag $Name: geant4-09-01-patch-02 $
28//
29//      ---------------- G4QCollision class -----------------
30//                 by Mikhail Kossov, December 2003.
31// G4QCollision class of the CHIPS Simulation Branch in GEANT4
32// ---------------------------------------------------------------
33// ****************************************************************************************
34// ********** This CLASS is temporary moved from the photolepton_hadron directory *********
35// ****************************************************************************************
36
37//#define debug
38//#define pdebug
39//#define ppdebug
40//#define qedebug
41
42#include "G4QCollision.hh"
43
44// Initialization of static vectors
45std::vector<G4int> G4QCollision::ElementZ;            // Z of the element(i) in theLastCalc
46std::vector<G4double> G4QCollision::ElProbInMat;      // SumProbabilityElements in Material
47std::vector<std::vector<G4int>*> G4QCollision::ElIsoN;    // N of isotope(j) of Element(i)
48std::vector<std::vector<G4double>*>G4QCollision::IsoProbInEl;//SumProbabIsotopes inElementI
49
50G4QCollision::G4QCollision(const G4String& processName) : G4VDiscreteProcess(processName)
51{
52#ifdef debug
53  G4cout<<"G4QCollision::Constructor is called"<<G4endl;
54#endif
55  if (verboseLevel>0) G4cout << GetProcessName() << " process is created "<< G4endl;
56
57  //G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPSWorld (234 part.max)
58  G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio); // Clusterization param's
59  G4Quasmon::SetParameters(Temperature,SSin2Gluons,EtaEtaprime);  // Hadronic parameters
60  G4QEnvironment::SetParameters(SolidAngle); // SolAngle of pbar-A secondary mesons capture
61  //@@ Initialize here the G4QuasmonString parameters
62}
63
64G4bool   G4QCollision::manualFlag=false; // If false then standard parameters are used
65G4double G4QCollision::Temperature=180.; // Critical Temperature (sensitive at High En)
66G4double G4QCollision::SSin2Gluons=0.3;  // Supression of s-quarks (in respect to u&d)
67G4double G4QCollision::EtaEtaprime=0.3;  // Supression of eta mesons (gg->qq/3g->qq)
68G4double G4QCollision::freeNuc=0.5;      // Percentage of free nucleons on the surface
69G4double G4QCollision::freeDib=0.05;     // Percentage of free diBaryons on the surface
70G4double G4QCollision::clustProb=5.;     // Nuclear clusterization parameter
71G4double G4QCollision::mediRatio=10.;    // medium/vacuum hadronization ratio
72G4int    G4QCollision::nPartCWorld=152;  // The#of particles initialized in CHIPS World
73G4double G4QCollision::SolidAngle=0.5;   // Part of Solid Angle to capture (@@A-dep.)
74G4bool   G4QCollision::EnergyFlux=false; // Flag for Energy Flux use (not MultyQuasmon)
75G4double G4QCollision::PiPrThresh=141.4; // Pion Production Threshold for gammas
76G4double G4QCollision::M2ShiftVir=20000.;// Shift for M2=-Q2=m_pi^2 of the virtualGamma
77G4double G4QCollision::DiNuclMass=1880.; // DoubleNucleon Mass for VirtualNormalization
78G4double G4QCollision::photNucBias=1.;   // BiasingParameter for photo(e,mu,tau)Nuclear
79G4double G4QCollision::weakNucBias=1.;   // BiasingParameter for ChargedCurrents(nu,mu)
80
81void G4QCollision::SetManual()   {manualFlag=true;}
82void G4QCollision::SetStandard() {manualFlag=false;}
83
84// Fill the private parameters
85void G4QCollision::SetParameters(G4double temper, G4double ssin2g, G4double etaetap,
86                                     G4double fN, G4double fD, G4double cP, G4double mR,
87                                     G4int nParCW, G4double solAn, G4bool efFlag,
88                                     G4double piThresh, G4double mpisq, G4double dinum)
89{//  =============================================================================
90  Temperature=temper;
91  SSin2Gluons=ssin2g;
92  EtaEtaprime=etaetap;
93  freeNuc=fN;
94  freeDib=fD;
95  clustProb=cP;
96  mediRatio=mR;
97  nPartCWorld = nParCW;
98  EnergyFlux=efFlag;
99  SolidAngle=solAn;
100  PiPrThresh=piThresh;
101  M2ShiftVir=mpisq;
102  DiNuclMass=dinum;
103  G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World with 234 particles
104  G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio); // Clusterization param's
105  G4Quasmon::SetParameters(Temperature,SSin2Gluons,EtaEtaprime);  // Hadronic parameters
106  G4QEnvironment::SetParameters(SolidAngle); // SolAngle of pbar-A secondary mesons capture
107}
108
109void G4QCollision::SetPhotNucBias(G4double phnB) {photNucBias=phnB;}
110void G4QCollision::SetWeakNucBias(G4double ccnB) {weakNucBias=ccnB;}
111
112// Destructor
113
114G4QCollision::~G4QCollision() {}
115
116
117G4LorentzVector G4QCollision::GetEnegryMomentumConservation()
118{
119  return EnMomConservation;
120}
121
122G4int G4QCollision::GetNumberOfNeutronsInTarget()
123{
124  return nOfNeutrons;
125}
126
127// output of the function must be in units of length! L=1/sig_V,sig_V=SUM(n(j,i)*sig(j,i)),
128// where n(i,j) is a number of nuclei of the isotop j of the element i in V=1(lengtUnit^3)
129// ********** All CHIPS cross sections are calculated in the surface units ************
130G4double G4QCollision::GetMeanFreePath(const G4Track& aTrack,G4double,G4ForceCondition* Fc)
131{
132#ifdef debug
133  G4cout<<"G4QCollision::GetMeanFreePath: Called Fc="<<*Fc<<G4endl;
134#endif
135  *Fc = NotForced;
136#ifdef debug
137  G4cout<<"G4QCollision::GetMeanFreePath: Before GetDynPart"<<G4endl;
138#endif
139  const G4DynamicParticle* incidentParticle = aTrack.GetDynamicParticle();
140#ifdef debug
141  G4cout<<"G4QCollision::GetMeanFreePath: Before GetDef"<<G4endl;
142#endif
143  G4ParticleDefinition* incidentParticleDefinition=incidentParticle->GetDefinition();
144  if( !IsApplicable(*incidentParticleDefinition))
145    G4cout<<"-W-G4QCollision::GetMeanFreePath called for not implemented particle"<<G4endl;
146  // Calculate the mean Cross Section for the set of Elements(*Isotopes) in the Material
147  G4double Momentum = incidentParticle->GetTotalMomentum(); // 3-momentum of the Particle
148#ifdef debug
149  G4cout<<"G4QCollis::GetMeanFreePath: BeforeGetMaterial"<<G4endl;
150#endif
151  const G4Material* material = aTrack.GetMaterial();        // Get the current material
152  const G4double* NOfNucPerVolume = material->GetVecNbOfAtomsPerVolume();
153  const G4ElementVector* theElementVector = material->GetElementVector();
154  G4int nE=material->GetNumberOfElements();
155#ifdef debug
156  G4cout<<"G4QCollision::GetMeanFreePath:"<<nE<<" Elem's in theMaterial"<<G4endl;
157#endif
158  G4bool leptoNuc=false;       // By default the reaction is not lepto-nuclear
159  G4VQCrossSection* CSmanager=0;
160  G4VQCrossSection* CSmanager2=0;
161  G4int pPDG=0;
162  if(incidentParticleDefinition == G4Proton::Proton())
163  {
164    CSmanager=G4QProtonNuclearCrossSection::GetPointer();
165    pPDG=2212;
166  }
167  else if(incidentParticleDefinition == G4Gamma::Gamma())
168  {
169    CSmanager=G4QPhotonNuclearCrossSection::GetPointer();
170    pPDG=22;
171  }
172  else if(incidentParticleDefinition == G4Electron::Electron() ||
173          incidentParticleDefinition == G4Positron::Positron())
174  {
175    CSmanager=G4QElectronNuclearCrossSection::GetPointer();
176    leptoNuc=true;
177    pPDG=11;
178  }
179  else if(incidentParticleDefinition == G4MuonPlus::MuonPlus() ||
180          incidentParticleDefinition == G4MuonMinus::MuonMinus())
181  {
182    CSmanager=G4QMuonNuclearCrossSection::GetPointer();
183    leptoNuc=true;
184    pPDG=13;
185  }
186  else if(incidentParticleDefinition == G4TauPlus::TauPlus() ||
187          incidentParticleDefinition == G4TauMinus::TauMinus())
188  {
189    CSmanager=G4QTauNuclearCrossSection::GetPointer();
190    leptoNuc=true;
191    pPDG=15;
192  }
193  else if(incidentParticleDefinition == G4NeutrinoMu::NeutrinoMu() )
194  {
195    CSmanager=G4QNuMuNuclearCrossSection::GetPointer();
196    CSmanager2=G4QNuNuNuclearCrossSection::GetPointer();
197    leptoNuc=true;
198    pPDG=14;
199  }
200  else if(incidentParticleDefinition == G4AntiNeutrinoMu::AntiNeutrinoMu() )
201  {
202    CSmanager=G4QANuMuNuclearCrossSection::GetPointer();
203    CSmanager2=G4QANuANuNuclearCrossSection::GetPointer();
204    leptoNuc=true;
205    pPDG=-14;
206  }
207  else if(incidentParticleDefinition == G4NeutrinoE::NeutrinoE() )
208  {
209    CSmanager=G4QNuENuclearCrossSection::GetPointer();
210    CSmanager2=G4QNuNuNuclearCrossSection::GetPointer();
211    leptoNuc=true;
212    pPDG=12;
213  }
214  else if(incidentParticleDefinition == G4AntiNeutrinoE::AntiNeutrinoE() )
215  {
216    CSmanager=G4QANuENuclearCrossSection::GetPointer();
217    CSmanager2=G4QANuANuNuclearCrossSection::GetPointer();
218    leptoNuc=true;
219    pPDG=-12;
220  }
221  else G4cout<<"G4QCollision::GetMeanFreePath:Particle isn't implemented in CHIPS"<<G4endl;
222 
223  G4QIsotope* Isotopes = G4QIsotope::Get(); // Pointer to the G4QIsotopes singleton
224  G4double sigma=0.;                        // Sums over elements for the material
225  G4int IPIE=IsoProbInEl.size();            // How many old elements?
226  if(IPIE) for(G4int ip=0; ip<IPIE; ++ip)   // Clean up the SumProb's of Isotopes (SPI)
227  {
228    std::vector<G4double>* SPI=IsoProbInEl[ip]; // Pointer to the SPI vector
229    SPI->clear();
230    delete SPI;
231    std::vector<G4int>* IsN=ElIsoN[ip];     // Pointer to the N vector
232    IsN->clear();
233    delete IsN;
234  }
235  ElProbInMat.clear();                      // Clean up the SumProb's of Elements (SPE)
236  ElementZ.clear();                         // Clear the body vector for Z of Elements
237  IsoProbInEl.clear();                      // Clear the body vector for SPI
238  ElIsoN.clear();                           // Clear the body vector for N of Isotopes
239  for(G4int i=0; i<nE; ++i)
240  {
241    G4Element* pElement=(*theElementVector)[i]; // Pointer to the current element
242    G4int Z = static_cast<G4int>(pElement->GetZ()); // Z of the Element
243    ElementZ.push_back(Z);                  // Remember Z of the Element
244    G4int isoSize=0;                        // The default for the isoVectorLength is 0
245    G4int indEl=0;                          // Index of non-trivial element or 0(default)
246    G4IsotopeVector* isoVector=pElement->GetIsotopeVector(); // Get the predefined IsoVect
247    if(isoVector) isoSize=isoVector->size();// Get size of the existing isotopeVector
248#ifdef debug
249    G4cout<<"G4QCollision::GetMeanFreePath: isovectorLength="<<isoSize<<G4endl; // Result
250#endif
251    if(isoSize)                             // The Element has non-trivial abumdance set
252    {
253      indEl=pElement->GetIndex()+1;         // Index of the non-trivial element
254      if(!Isotopes->IsDefined(Z,indEl))     // This index is not defined for this Z: define
255      {
256        std::vector<std::pair<G4int,G4double>*>* newAbund =
257                                               new std::vector<std::pair<G4int,G4double>*>;
258        G4double* abuVector=pElement->GetRelativeAbundanceVector();
259        for(G4int j=0; j<isoSize; j++)      // Calculation of abundance vector for isotopes
260        {
261          G4int N=pElement->GetIsotope(j)->GetN()-Z; // N means A=N+Z !
262          if(pElement->GetIsotope(j)->GetZ()!=Z)G4cerr<<"G4QCollision::GetMeanFreePath"
263                                                                                                                                                                                                                                                                        <<": Z="<<pElement->GetIsotope(j)->GetZ()<<"#"<<Z<<G4endl;
264          G4double abund=abuVector[j];
265                                                                  std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund);
266#ifdef debug
267          G4cout<<"G4QCollision::GetMeanFreePath: p#="<<j<<",N="<<N<<",ab="<<abund<<G4endl;
268#endif
269          newAbund->push_back(pr);
270                                                  }
271#ifdef debug
272        G4cout<<"G4QCollision::GetMeanFreePath: pairVectLength="<<newAbund->size()<<G4endl;
273#endif
274        indEl=G4QIsotope::Get()->InitElement(Z,indEl,newAbund); // definition of the newInd
275        for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k];   // Cleaning temporary
276        delete newAbund; // Was "new" in the beginning of the name space
277      }
278    }
279    std::vector<std::pair<G4int,G4double>*>* cs= Isotopes->GetCSVector(Z,indEl);//CSPointer
280    std::vector<G4double>* SPI = new std::vector<G4double>; // Pointer to the SPI vector
281    IsoProbInEl.push_back(SPI);
282    std::vector<G4int>* IsN = new std::vector<G4int>; // Pointer to the N vector
283    ElIsoN.push_back(IsN);
284    G4int nIs=cs->size();                   // A#Of Isotopes in the Element
285    G4double susi=0.;                       // sum of CS over isotopes
286    if(nIs) for(G4int j=0; j<nIs; j++)      // Calculate CS for eachIsotope of El
287    {
288      std::pair<G4int,G4double>* curIs=(*cs)[j]; // A pointer, which is used twice
289      G4int N=curIs->first;                 // #of Neuterons in the isotope j of El i
290      IsN->push_back(N);                    // Remember Min N for the Element
291      G4double CSI=CSmanager->GetCrossSection(true,Momentum,Z,N,pPDG);//CS(j,i) for isotope
292      if(CSmanager2)CSI+=CSmanager2->GetCrossSection(true,Momentum,Z,N,pPDG);//CS(j,i)nu,nu
293#ifdef debug
294      G4cout<<"GQC::GMF:X="<<CSI<<",M="<<Momentum<<",Z="<<Z<<",N="<<N<<",P="<<pPDG<<G4endl;
295#endif
296      curIs->second = CSI;
297      susi+=CSI;                            // Make a sum per isotopes
298      SPI->push_back(susi);                 // Remember summed cross-section
299    } // End of temporary initialization of the cross sections in the G4QIsotope singeltone
300    sigma+=Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i];//SUM(MeanCS*NOfNperV)
301    ElProbInMat.push_back(sigma);
302  } // End of LOOP over Elements
303#ifdef debug
304  G4cout<<"G4QCol::GetMeanFrPa: S="<<sigma<<",e="<<photNucBias<<",w="<<weakNucBias<<G4endl;
305#endif
306  // Check that cross section is not zero and return the mean free path
307  if(photNucBias!=1.) if(incidentParticleDefinition == G4Gamma::Gamma()         ||
308                         incidentParticleDefinition == G4MuonPlus::MuonPlus()   ||
309                         incidentParticleDefinition == G4MuonMinus::MuonMinus() ||
310                         incidentParticleDefinition == G4Electron::Electron()   ||
311                         incidentParticleDefinition == G4Positron::Positron()   || 
312                         incidentParticleDefinition == G4TauMinus::TauMinus()   ||
313                         incidentParticleDefinition == G4TauPlus::TauPlus()       )
314                                                                        sigma*=photNucBias;
315  if(weakNucBias!=1.) if(incidentParticleDefinition==G4NeutrinoE::NeutrinoE()            ||
316                         incidentParticleDefinition==G4AntiNeutrinoE::AntiNeutrinoE()    ||
317                         incidentParticleDefinition==G4NeutrinoTau::NeutrinoTau()        ||
318                         incidentParticleDefinition==G4AntiNeutrinoTau::AntiNeutrinoTau()||
319                         incidentParticleDefinition==G4NeutrinoMu::NeutrinoMu()          ||
320                         incidentParticleDefinition==G4AntiNeutrinoMu::AntiNeutrinoMu()   )
321                                                                        sigma*=weakNucBias;
322  if(sigma > 0.) return 1./sigma;                 // Mean path [distance]
323  return DBL_MAX;
324}
325
326
327G4bool G4QCollision::IsApplicable(const G4ParticleDefinition& particle) 
328{
329  if      (particle == *(      G4MuonPlus::MuonPlus()      )) return true;
330  else if (particle == *(     G4MuonMinus::MuonMinus()     )) return true; 
331  else if (particle == *(       G4TauPlus::TauPlus()       )) return true;
332  else if (particle == *(      G4TauMinus::TauMinus()      )) return true;
333  else if (particle == *(      G4Electron::Electron()      )) return true;
334  else if (particle == *(      G4Positron::Positron()      )) return true;
335  else if (particle == *(         G4Gamma::Gamma()         )) return true;
336  else if (particle == *(        G4Proton::Proton()        )) return true;
337  else if (particle == *( G4AntiNeutrinoE::AntiNeutrinoE() )) return true;
338  else if (particle == *(     G4NeutrinoE::NeutrinoE()     )) return true;
339  else if (particle == *(G4AntiNeutrinoMu::AntiNeutrinoMu())) return true;
340  else if (particle == *(    G4NeutrinoMu::NeutrinoMu()    )) return true;
341  //else if (particle == *(       G4Neutron::Neutron()       )) return true;
342  //else if (particle == *(     G4PionMinus::PionMinus()     )) return true;
343  //else if (particle == *(      G4PionPlus::PionPlus()      )) return true;
344  //else if (particle == *(      G4KaonPlus::KaonPlus()      )) return true;
345  //else if (particle == *(     G4KaonMinus::KaonMinus()     )) return true;
346  //else if (particle == *(  G4KaonZeroLong::KaonZeroLong()  )) return true;
347  //else if (particle == *( G4KaonZeroShort::KaonZeroShort() )) return true;
348  //else if (particle == *(        G4Lambda::Lambda()        )) return true;
349  //else if (particle == *(     G4SigmaPlus::SigmaPlus()     )) return true;
350  //else if (particle == *(    G4SigmaMinus::SigmaMinus()    )) return true;
351  //else if (particle == *(     G4SigmaZero::SigmaZero()     )) return true;
352  //else if (particle == *(       G4XiMinus::XiMinus()       )) return true;
353  //else if (particle == *(        G4XiZero::XiZero()        )) return true;
354  //else if (particle == *(    G4OmegaMinus::OmegaMinus()    )) return true;
355  //else if (particle == *(   G4AntiNeutron::AntiNeutron()   )) return true;
356  //else if (particle == *(    G4AntiProton::AntiProton()    )) return true;
357  //else if (particle == *(G4AntiNeutrinoTau::AntiNeutrinoTau())) return true;
358  //else if (particle == *(    G4NeutrinoTau::NeutrinoTau()    )) return true;
359#ifdef debug
360  G4cout<<"***G4QCollision::IsApplicable: PDG="<<particle.GetPDGEncoding()<<G4endl;
361#endif
362  return false;
363}
364
365G4VParticleChange* G4QCollision::PostStepDoIt(const G4Track& track, const G4Step& step)
366{
367  static const G4double third = 1./3.;
368  static const G4double me=G4Electron::Electron()->GetPDGMass();   // electron mass
369  static const G4double me2=me*me;                                 // squared electron mass
370  static const G4double mu=G4MuonMinus::MuonMinus()->GetPDGMass(); // muon mass
371  static const G4double mu2=mu*mu;                                 // squared muon mass
372  static const G4double mt=G4TauMinus::TauMinus()->GetPDGMass();   // tau mass
373  static const G4double mt2=mt*mt;                                 // squared tau mass
374  //static const G4double dpi=M_PI+M_PI;   // 2*pi (for Phi distr.) ***changed to twopi***
375  static const G4double mNeut= G4QPDGCode(2112).GetMass();
376  static const G4double mNeut2= mNeut*mNeut;
377  static const G4double mProt= G4QPDGCode(2212).GetMass();
378  static const G4double mProt2= mProt*mProt;
379  static const G4double dM=mProt+mNeut;                            // doubled nucleon mass
380  static const G4double hdM=dM/2.;                                 // M of the "nucleon"
381  static const G4double hdM2=hdM*hdM;                              // M2 of the "nucleon"
382  static const G4double mPi0 = G4QPDGCode(111).GetMass();
383  static const G4double mPi0s= mPi0*mPi0;
384  static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0);// Mass of deuteron
385  static const G4double mTrit= G4QPDGCode(2112).GetNuclMass(1,2,0);// Mass of tritium
386  static const G4double mHel3= G4QPDGCode(2112).GetNuclMass(2,1,0);// Mass of Helium3
387  static const G4double mAlph= G4QPDGCode(2112).GetNuclMass(2,2,0);// Mass of alpha
388  static const G4double mPi  = G4QPDGCode(211).GetMass();
389  static const G4double tmPi = mPi+mPi;     // Doubled mass of the charged pion
390  static const G4double stmPi= tmPi*tmPi;   // Squared Doubled mass of the charged pion
391  static const G4double mPPi = mPi+mProt;   // Delta threshold
392  static const G4double mPPi2= mPPi*mPPi; // Delta low threshold for W2
393  //static const G4double mDel2= 1400*1400; // Delta up threshold for W2 (in MeV^2)
394  // Static definitions for electrons (nu,e) -----------------------------------------
395  static const G4double meN = mNeut+me;
396  static const G4double meN2= meN*meN;
397  static const G4double fmeN= 4*mNeut2*me2;
398  static const G4double mesN= mNeut2+me2;
399  static const G4double meP = mProt+me;
400  static const G4double meP2= meP*meP;
401  static const G4double fmeP= 4*mProt2*me2;
402  static const G4double mesP= mProt2+me2;
403  static const G4double medM= me2/dM;       // for x limit
404  static const G4double meD = mPPi+me;      // Multiperipheral threshold
405  static const G4double meD2= meD*meD;
406  // Static definitions for muons (nu,mu) -----------------------------------------
407  static const G4double muN = mNeut+mu;
408  static const G4double muN2= muN*muN;
409  static const G4double fmuN= 4*mNeut2*mu2;
410  static const G4double musN= mNeut2+mu2;
411  static const G4double muP = mProt+mu;
412  static const G4double muP2= muP*muP;      // +
413  static const G4double fmuP= 4*mProt2*mu2; // +
414  static const G4double musP= mProt2+mu2;
415  static const G4double mudM= mu2/dM;       // for x limit
416  static const G4double muD = mPPi+mu;      // Multiperipheral threshold
417  static const G4double muD2= muD*muD;
418  // Static definitions for muons (nu,nu) -----------------------------------------
419  //static const G4double nuN = mNeut;
420  //static const G4double nuN2= mNeut2;
421  //static const G4double fnuN= 0.;
422  //static const G4double nusN= mNeut2;
423  //static const G4double nuP = mProt;
424  //static const G4double nuP2= mProt2;
425  //static const G4double fnuP= 0.;
426  //static const G4double nusP= mProt2;
427  //static const G4double nudM= 0.;           // for x limit
428  //static const G4double nuD = mPPi;         // Multiperipheral threshold
429  //static const G4double nuD2= mPPi2;
430  //-------------------------------------------------------------------------------------
431  static G4bool CWinit = true;              // CHIPS Warld needs to be initted
432  if(CWinit)
433                {
434    CWinit=false;
435    G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World (234 part.max)
436  }
437  //-------------------------------------------------------------------------------------
438  const G4DynamicParticle* projHadron = track.GetDynamicParticle();
439  const G4ParticleDefinition* particle=projHadron->GetDefinition();
440#ifdef debug
441  G4cout<<"G4QCollision::PostStepDoIt: Before the GetMeanFreePath is called"<<G4endl;
442#endif
443  G4ForceCondition cond=NotForced;
444  GetMeanFreePath(track, 1., &cond);
445#ifdef debug
446  G4cout<<"G4QCollision::PostStepDoIt: After the GetMeanFreePath is called"<<G4endl;
447#endif
448  G4bool scat=false;                                  // No CHEX in proj scattering
449  G4int  scatPDG=0;                                   // Must be filled if true (CHEX)
450  G4LorentzVector proj4M=projHadron->Get4Momentum();  // 4-momentum of the projectile (IU?)
451  G4LorentzVector scat4M=proj4M;                      // Must be filled if true
452  G4double momentum = projHadron->GetTotalMomentum(); // 3-momentum of the Particle
453  G4double Momentum=proj4M.rho();
454  if(std::fabs(Momentum-momentum)>.001)
455    G4cerr<<"*G4QCollision::PostStepDoIt: P="<<Momentum<<"#"<<momentum<<G4endl;
456#ifdef debug
457  G4double mp=proj4M.m();
458  G4cout<<"G4QCollis::PostStepDoIt:called, P="<<Momentum<<"="<<momentum<<",m="<<mp<<G4endl;
459#endif
460  if (!IsApplicable(*particle))  // Check applicability
461  {
462    G4cerr<<"G4QCollision::PostStepDoIt:Only gam,e+,e-,mu+,mu-,t+,t-,p are implemented."
463          <<G4endl;
464    return 0;
465  }
466  const G4Material* material = track.GetMaterial();      // Get the current material
467  G4int Z=0;
468  const G4ElementVector* theElementVector = material->GetElementVector();
469  G4int nE=material->GetNumberOfElements();
470#ifdef debug
471  G4cout<<"G4QCollision::PostStepDoIt: "<<nE<<" elements in the material."<<G4endl;
472#endif
473  G4int projPDG=0;                           // PDG Code prototype for the captured hadron
474  // Not all these particles are implemented yet (see Is Applicable)
475  if      (particle ==        G4MuonPlus::MuonPlus()       ) projPDG=  -13;
476  else if (particle ==       G4MuonMinus::MuonMinus()      ) projPDG=   13;
477  else if (particle ==      G4NeutrinoMu::NeutrinoMu()     ) projPDG=   14;
478  else if (particle ==  G4AntiNeutrinoMu::AntiNeutrinoMu() ) projPDG=  -14;
479  else if (particle ==        G4Electron::Electron()       ) projPDG=   11;
480  else if (particle ==        G4Positron::Positron()       ) projPDG=  -11;
481  else if (particle ==       G4NeutrinoE::NeutrinoE()      ) projPDG=   12;
482  else if (particle ==   G4AntiNeutrinoE::AntiNeutrinoE()  ) projPDG=  -12;
483  else if (particle ==           G4Gamma::Gamma()          ) projPDG=   22;
484  else if (particle ==          G4Proton::Proton()         ) projPDG= 2212;
485  else if (particle ==         G4Neutron::Neutron()        ) projPDG= 2112;
486  else if (particle ==       G4PionMinus::PionMinus()      ) projPDG= -211;
487  else if (particle ==        G4PionPlus::PionPlus()       ) projPDG=  211;
488  else if (particle ==        G4KaonPlus::KaonPlus()       ) projPDG= 2112;
489  else if (particle ==       G4KaonMinus::KaonMinus()      ) projPDG= -321;
490  else if (particle ==    G4KaonZeroLong::KaonZeroLong()   ) projPDG=  130;
491  else if (particle ==   G4KaonZeroShort::KaonZeroShort()  ) projPDG=  310;
492  else if (particle ==         G4TauPlus::TauPlus()        ) projPDG=  -15;
493  else if (particle ==        G4TauMinus::TauMinus()       ) projPDG=   15;
494  else if (particle ==     G4NeutrinoTau::NeutrinoTau()    ) projPDG=   16;
495  else if (particle == G4AntiNeutrinoTau::AntiNeutrinoTau()) projPDG=  -16;
496  else if (particle ==          G4Lambda::Lambda()         ) projPDG= 3122;
497  else if (particle ==       G4SigmaPlus::SigmaPlus()      ) projPDG= 3222;
498  else if (particle ==      G4SigmaMinus::SigmaMinus()     ) projPDG= 3112;
499  else if (particle ==       G4SigmaZero::SigmaZero()      ) projPDG= 3212;
500  else if (particle ==         G4XiMinus::XiMinus()        ) projPDG= 3312;
501  else if (particle ==          G4XiZero::XiZero()         ) projPDG= 3322;
502  else if (particle ==      G4OmegaMinus::OmegaMinus()     ) projPDG= 3334;
503  else if (particle ==     G4AntiNeutron::AntiNeutron()    ) projPDG=-2112;
504  else if (particle ==      G4AntiProton::AntiProton()     ) projPDG=-2212;
505  G4int aProjPDG=std::abs(projPDG);
506#ifdef debug
507  G4int prPDG=particle->GetPDGEncoding();
508                G4cout<<"G4QCollision::PostStepDoIt: projPDG="<<projPDG<<", stPDG="<<prPDG<<G4endl;
509#endif
510  if(!projPDG)
511  {
512    G4cerr<<"---Warning---G4QCollision::PostStepDoIt:Undefined interacting hadron"<<G4endl;
513    return 0;
514  }
515  G4int EPIM=ElProbInMat.size();
516#ifdef debug
517                G4cout<<"G4QCollis::PostStDoIt: m="<<EPIM<<",n="<<nE<<",T="<<ElProbInMat[EPIM-1]<<G4endl;
518#endif
519  G4int i=0;
520  if(EPIM>1)
521  {
522    G4double rnd = ElProbInMat[EPIM-1]*G4UniformRand();
523    for(i=0; i<nE; ++i)
524                  {
525#ifdef debug
526                                  G4cout<<"G4QCollision::PostStepDoIt:E["<<i<<"]="<<ElProbInMat[i]<<",r="<<rnd<<G4endl;
527#endif
528      if (rnd<ElProbInMat[i]) break;
529    }
530    if(i>=nE) i=nE-1;                        // Top limit for the Element
531  }
532  G4Element* pElement=(*theElementVector)[i];
533  Z=static_cast<G4int>(pElement->GetZ());
534#ifdef debug
535                                G4cout<<"G4QCollision::PostStepDoIt: i="<<i<<", Z(element)="<<Z<<G4endl;
536#endif
537  if(Z<=0)
538  {
539    G4cerr<<"---Warning---G4QCollision::PostStepDoIt: Element with Z="<<Z<<G4endl;
540    if(Z<0) return 0;
541  }
542  std::vector<G4double>* SPI = IsoProbInEl[i];// Vector of summedProbabilities for isotopes
543  std::vector<G4int>* IsN = ElIsoN[i];     // Vector of "#of neutrons" in the isotope El[i]
544  G4int nofIsot=SPI->size();               // #of isotopes in the element i
545#ifdef debug
546                G4cout<<"G4QCollis::PosStDoIt:n="<<nofIsot<<",T="<<(*SPI)[nofIsot-1]<<G4endl;
547#endif
548  G4int j=0;
549  if(nofIsot>1)
550  {
551    G4double rndI=(*SPI)[nofIsot-1]*G4UniformRand(); // Randomize the isotop of the Element
552    for(j=0; j<nofIsot; ++j)
553    {
554#ifdef debug
555                                  G4cout<<"G4QCollision::PostStepDoIt: SP["<<j<<"]="<<(*SPI)[j]<<", r="<<rndI<<G4endl;
556#endif
557      if(rndI < (*SPI)[j]) break;
558    }
559    if(j>=nofIsot) j=nofIsot-1;            // Top limit for the isotope
560  }
561  G4int N =(*IsN)[j]; ;                    // Randomized number of neutrons
562#ifdef debug
563                G4cout<<"G4QCollision::PostStepDoIt: j="<<i<<", N(isotope)="<<N<<G4endl;
564#endif
565  if(N<0)
566  {
567    G4cerr<<"-Warning-G4QCollision::PostStepDoIt: Isotope with Z="<<Z<<", 0>N="<<N<<G4endl;
568    return 0;
569  }
570  nOfNeutrons=N;                           // Remember it for the energy-momentum check
571  G4double dd=0.025;
572  G4double am=Z+N;
573  G4double sr=std::sqrt(am);
574  G4double dsr=0.01*(sr+sr);
575  if(dsr<dd)dsr=dd;
576  if(manualFlag) G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio); // ManualP
577                else if(projPDG==-2212) G4QNucleus::SetParameters(1.-dsr-dsr,dd+dd,5.,10.);//aP ClustPars
578  else if(projPDG==-211)  G4QNucleus::SetParameters(.67-dsr,.32-dsr,5.,9.);//Pi- ClustPars
579#ifdef debug
580  G4cout<<"G4QCollision::PostStepDoIt: N="<<N<<" for element with Z="<<Z<<G4endl;
581#endif
582  if(N<0)
583  {
584    G4cerr<<"---Warning---G4QCollision::PostStepDoIt:Element with N="<<N<< G4endl;
585    return 0;
586  }
587  aParticleChange.Initialize(track);
588  G4double weight = track.GetWeight();
589  if(photNucBias!=1.)      weight/=photNucBias;
590  else if(weakNucBias!=1.) weight/=weakNucBias;
591  G4double localtime = track.GetGlobalTime();
592  G4ThreeVector position = track.GetPosition();
593  G4TouchableHandle trTouchable = track.GetTouchableHandle();
594  //
595  G4int targPDG=90000000+Z*1000+N;            // PDG Code of the target nucleus
596  G4QPDGCode targQPDG(targPDG);
597  G4double tgM=targQPDG.GetMass();            // Target mass
598  G4double tM=tgM;                            // Target mass (copy to be changed)
599  G4QHadronVector* output=new G4QHadronVector;// Prototype of EnvironOutput G4QHadronVector
600  G4double absMom = 0.;                       // Prototype of absorbed by nucleus Moment
601  G4QHadronVector* leadhs=new G4QHadronVector;// Prototype of QuasmOutput G4QHadronVectorum
602  G4LorentzVector lead4M(0.,0.,0.,0.);        // Prototype of LeadingQ 4-momentum
603
604  // 
605  // Leptons with photonuclear
606  // Lepto-nuclear case with the equivalent photon algorithm. @@InFuture + NC (?)
607  //
608  if (aProjPDG == 11 || aProjPDG == 13 || aProjPDG == 15) {
609
610#ifdef debug
611    G4cout<<"G4QCollision::PostStDoIt:startSt="<<aParticleChange.GetTrackStatus()<<G4endl;
612#endif
613    G4double kinEnergy= projHadron->GetKineticEnergy();
614    G4ParticleMomentum dir = projHadron->GetMomentumDirection();
615    G4VQCrossSection* CSmanager=G4QElectronNuclearCrossSection::GetPointer();
616    G4double ml=me;
617    G4double ml2=me2;
618    if(aProjPDG== 13)
619    {
620      CSmanager=G4QMuonNuclearCrossSection::GetPointer();
621      ml=mu;
622      ml2=mu2;
623    }
624    if(aProjPDG== 15)
625    {
626      CSmanager=G4QTauNuclearCrossSection::GetPointer();
627      ml=mt;
628      ml2=mt2;
629    }
630    // @@ Probably this is not necessary any more (?)
631    G4double xSec=CSmanager->GetCrossSection(false,Momentum,Z,N,aProjPDG);// Recalculate XS
632    // @@ check a possibility to separate p, n, or alpha (!)
633    if(xSec <= 0.) // The cross-section is 0 -> Do Nothing
634    {
635#ifdef debug
636      G4cerr<<"---OUT---G4QCollision::PSDoIt: Called for zero Cross-section"<<G4endl;
637#endif
638      //Do Nothing Action insead of the reaction
639      aParticleChange.ProposeEnergy(kinEnergy);
640      aParticleChange.ProposeLocalEnergyDeposit(0.);
641      aParticleChange.ProposeMomentumDirection(dir);
642      aParticleChange.ProposeTrackStatus(fAlive);
643      return G4VDiscreteProcess::PostStepDoIt(track,step);
644    }
645    G4double photonEnergy = CSmanager->GetExchangeEnergy(); // Energy of EqivExchangePart
646#ifdef debug
647                                G4cout<<"G4QCol::PStDoIt: kE="<<kinEnergy<<",dir="<<dir<<",phE="<<photonEnergy<<G4endl;
648#endif
649    if( kinEnergy < photonEnergy || photonEnergy < 0.)
650    {
651      //Do Nothing Action insead of the reaction
652      G4cerr<<"--G4QCollision::PSDoIt: photE="<<photonEnergy<<">leptE="<<kinEnergy<<G4endl;
653      aParticleChange.ProposeEnergy(kinEnergy);
654      aParticleChange.ProposeLocalEnergyDeposit(0.);
655      aParticleChange.ProposeMomentumDirection(dir);
656      aParticleChange.ProposeTrackStatus(fAlive);
657      return G4VDiscreteProcess::PostStepDoIt(track,step);
658    }
659    G4double photonQ2 = CSmanager->GetExchangeQ2(photonEnergy);// Q2(t) of EqivExchangePart
660    G4double W=photonEnergy-photonQ2/dM;// HadronicEnergyFlow (W-energy) for virtual photon
661    if(tM<999.) W-=mPi0+mPi0s/dM;       // Pion production threshold for a nucleon target
662    if(W<0.) 
663    {
664      //Do Nothing Action insead of the reaction
665#ifdef debug
666      G4cout<<"--G4QCollision::PostStepDoIt:(lN) negative equivalent energy W="<<W<<G4endl;
667#endif
668      aParticleChange.ProposeEnergy(kinEnergy);
669      aParticleChange.ProposeLocalEnergyDeposit(0.);
670      aParticleChange.ProposeMomentumDirection(dir);
671      aParticleChange.ProposeTrackStatus(fAlive);
672      return G4VDiscreteProcess::PostStepDoIt(track,step);
673    }
674    // Update G4VParticleChange for the scattered muon
675    G4VQCrossSection* thePhotonData=G4QPhotonNuclearCrossSection::GetPointer();
676    G4double sigNu=thePhotonData->GetCrossSection(true,photonEnergy,Z,N,22);//Integrated XS
677    G4double sigK =thePhotonData->GetCrossSection(true, W, Z, N, 22);       // Real XS
678    G4double rndFraction = CSmanager->GetVirtualFactor(photonEnergy, photonQ2);
679    if(sigNu*G4UniformRand()>sigK*rndFraction) 
680    {
681      //Do NothingToDo Action insead of the reaction
682#ifdef debug
683      G4cout<<"-DoNoth-G4QCollision::PostStepDoIt: probab. correction - DoNothing"<<G4endl;
684#endif
685      aParticleChange.ProposeEnergy(kinEnergy);
686      aParticleChange.ProposeLocalEnergyDeposit(0.);
687      aParticleChange.ProposeMomentumDirection(dir);
688      aParticleChange.ProposeTrackStatus(fAlive);
689      return G4VDiscreteProcess::PostStepDoIt(track,step);
690    }
691    G4double iniE=kinEnergy+ml;          // Initial total energy of the lepton
692    G4double finE=iniE-photonEnergy;     // Final total energy of the lepton
693#ifdef pdebug
694    G4cout<<"G4QCollision::PoStDoIt:E="<<iniE<<",lE="<<finE<<"-"<<ml<<"="<<finE-ml<<G4endl;
695#endif
696    aParticleChange.ProposeEnergy(finE-ml);
697    if(finE<=ml)                         // Secondary lepton (e/mu/tau) at rest disappears
698    {
699      aParticleChange.ProposeEnergy(0.);
700      if(aProjPDG== 11) aParticleChange.ProposeTrackStatus(fStopAndKill);
701      else aParticleChange.ProposeTrackStatus(fStopButAlive);
702      aParticleChange.ProposeMomentumDirection(dir);
703    }
704    else aParticleChange.ProposeTrackStatus(fAlive);
705    G4double iniP=std::sqrt(iniE*iniE-ml2); // Initial momentum of the electron
706    G4double finP=std::sqrt(finE*finE-ml2); // Final momentum of the electron
707    G4double cost=(iniE*finE-ml2-photonQ2/2)/iniP/finP; // cos(scat_ang_of_lepton)
708#ifdef pdebug
709                  G4cout<<"G4QC::PSDoIt:Q2="<<photonQ2<<",ct="<<cost<<",Pi="<<iniP<<",Pf="<<finP<<G4endl;
710#endif
711    if(cost>1.) cost=1.;                 // To avoid the accuracy of calculation problem
712    if(cost<-1.) cost=-1.;               // To avoid the accuracy of calculation problem
713    //
714    // Scatter the lepton ( @@ make the same thing for real photons)
715    // At this point we have photonEnergy and photonQ2 (with notDefinedPhi)->SelectProjPart
716    G4double absEn = std::pow(am,third)*GeV;  // @@(b) Mean Energy Absorbed by a Nucleus
717    //if(am>1 && absEn < photonEnergy)     // --> the absorption of energy can happen
718                                if(absEn < photonEnergy)     // --> the absorption of energy can happen
719    {
720      G4double abtEn = absEn+hdM;        // @@(b) MeanEnergyAbsorbed by a nucleus (+M_N)
721      G4double abEn2 = abtEn*abtEn;      // Squared absorbed Energy + MN
722      G4double abMo2 = abEn2-hdM2;       // Squared absorbed Momentum of compound system
723      G4double phEn2 = photonEnergy*photonEnergy;
724      G4double phMo2 = phEn2+photonQ2;   // Squared momentum of primary virtual photon
725      G4double phMo  = std::sqrt(phMo2); // Momentum of the primary virtual photon
726      absMom         = std::sqrt(abMo2); // Absorbed Momentum
727      if(absMom < phMo)                  // --> the absorption of momentum can happen
728                                  {
729        G4double dEn = photonEnergy - absEn; // Leading energy
730        G4double dMo = phMo - absMom;    // Leading momentum
731        G4double sF  = dEn*dEn - dMo*dMo;// s of leading particle
732#ifdef ppdebug
733                                    G4cout<<"-PhotoAbsorption-G4QCol::PStDoIt:sF="<<sF<<",phEn="<<photonEnergy<<G4endl;
734#endif
735        if(sF > stmPi)                   // --> Leading fragmentation is possible
736                                                                {
737          photonEnergy = absEn;          // New value of the photon energy
738          photonQ2=abMo2-absEn*absEn;    // New value of the photon Q2
739          absEn = dEn;                   // Put energy of leading particle to absEn (!)
740        }
741        else absMom=0.;                  // Flag that nothing has happened
742      }
743      else absMom=0.;                    // Flag that nothing has happened
744    }
745    // ------------- End of ProjPart selection
746    //
747    // Scattering in respect to the derection of the incident muon is made impicitly:
748    G4ThreeVector ort=dir.orthogonal();  // Not normed orthogonal vector (!) (to dir)
749    G4ThreeVector ortx = ort.unit();     // First unit vector orthogonal to the direction
750    G4ThreeVector orty = dir.cross(ortx);// Second unit vector orthoganal to the direction
751    G4double sint=std::sqrt(1.-cost*cost); // Perpendicular component
752    G4double phi=twopi*G4UniformRand();  // phi of scattered electron
753    G4double sinx=sint*std::sin(phi);    // x perpendicular component
754    G4double siny=sint*std::cos(phi);    // y perpendicular component
755    G4ThreeVector findir=cost*dir+sinx*ortx+siny*orty;
756    aParticleChange.ProposeMomentumDirection(findir); // new direction for the lepton
757#ifdef pdebug
758                  G4cout<<"G4QCollision::PostStepDoIt: E="<<aParticleChange.GetEnergy()<<"="<<finE<<"-"
759          <<ml<<", d="<<*aParticleChange.GetMomentumDirection()<<","<<findir<<G4endl;
760#endif
761    G4ThreeVector photon3M=iniP*dir-finP*findir;// 3D total momentum of photon
762    if(absMom)                           // Photon must be reduced & LeadingSyst fragmented
763    {
764      G4double ptm=photon3M.mag();                    // 3M of the virtual photon
765#ifdef ppdebug
766                    G4cout<<"-Absorption-G4QCollision::PostStepDoIt: ph3M="<<photon3M<<", eIn3M="
767            <<iniP*dir<<", eFin3M="<<finP*findir<<", abs3M="<<absMom<<"<ptm="<<ptm<<G4endl;
768#endif
769      G4ThreeVector lead3M=photon3M*(ptm-absMom)/ptm; // Keep the direction for leading Q
770      photon3M-=lead3M; // Reduced photon Momentum (photEn already = absEn)
771      proj4M=G4LorentzVector(lead3M,absEn); // 4-momentum of leading System
772#ifdef ppdebug
773                    G4cout<<"-->G4QC::PoStDoIt: new sF="<<proj4M.m2()<<", lead4M="<<proj4M<<G4endl;
774#endif
775      lead4M=proj4M;                        // Remember 4-mom for the total 4-momentum
776      G4Quasmon* pan= new G4Quasmon(G4QContent(1,1,0,1,1,0),proj4M);// ---> DELETED -->---+
777      try                                                           //                    |
778             {                                                             //                    |
779               delete leadhs;                                              //                    |
780        G4QNucleus vac(90000000);                                   //                    |
781        leadhs=pan->Fragment(vac,1);  // DELETED after it is copied to output vector      |
782      }                                                             //                    |
783      catch (G4QException& error)                                   //                    |
784             {                                                             //                    |
785        G4cerr<<"***G4QCollision::PostStepDoIt: G4Quasmon Exception is catched"<<G4endl;//|
786        G4Exception("G4QCollision::PostStepDoIt:","72",FatalException,"QuasmonCrash");  //|
787      }                                                             //                    |
788      delete pan;                              // Delete the Nuclear Environment <----<---+
789#ifdef ppdebug
790                                                G4cout<<"G4QCol::PStDoIt: l4M="<<proj4M<<proj4M.m2()<<", N="<<leadhs->size()<<",pt="
791            <<ptm<<",pa="<<absMom<<",El="<<absEn<<",Pl="<<ptm-absMom<<G4endl;
792#endif
793    }
794    projPDG=22;
795    proj4M=G4LorentzVector(photon3M,photonEnergy);
796#ifdef debug
797                  G4cout<<"G4QCollision::PostStDoIt: St="<<aParticleChange.GetTrackStatus()<<", g4m="
798          <<proj4M<<", lE="<<finE<<", lP="<<finP*findir<<", d="<<findir.mag2()<<G4endl;
799#endif
800
801  //
802  // neutrinoNuclear interactions (nu_e, nu_mu only)
803  //
804  } else if (aProjPDG == 12 || aProjPDG == 14) {
805
806    G4double kinEnergy= projHadron->GetKineticEnergy()/MeV; // Total energy of the neutrino
807    G4double dKinE=kinEnergy+kinEnergy;  // doubled energy for s calculation
808#ifdef debug
809                  G4cout<<"G4QCollision::PostStDoIt: 2*nuEnergy="<<dKinE<<"(MeV), PDG="<<projPDG<<G4endl;
810#endif
811    G4ParticleMomentum dir = projHadron->GetMomentumDirection(); // unit vector
812    G4double ml  = mu;
813    G4double ml2 =      mu2;
814    //G4double mlN =    muN;
815    G4double mlN2=      muN2;
816    G4double fmlN=      fmuN;
817    G4double mlsN=      musN;
818    //G4double mlP =    muP;
819    G4double mlP2=      muP2;
820    G4double fmlP=      fmuP;
821    G4double mlsP=      musP;
822    G4double mldM=      mudM;
823    //G4double mlD =    muD;
824    G4double mlD2=      muD2;
825    if(aProjPDG==12)
826    {
827      ml  = me;
828      ml2 =     me2;
829      //mlN =   meN;
830      mlN2=     meN2;
831      fmlN=     fmeN;
832      mlsN=     mesN;
833      //mlP =   meP;
834      mlP2=     meP2;
835      fmlP=     fmeP;
836      mlsP=     mesP;
837      mldM=     medM;
838      //mlD =   meD;
839      mlD2=     meD2;
840    }
841    G4VQCrossSection* CSmanager =G4QNuMuNuclearCrossSection::GetPointer(); // (nu,l)
842    G4VQCrossSection* CSmanager2=G4QNuNuNuclearCrossSection::GetPointer(); // (nu,nu)
843    proj4M=G4LorentzVector(dir*kinEnergy,kinEnergy);   // temporary
844    G4bool nuanu=true;
845    scatPDG=13;                          // Prototype = secondary scattered mu-
846    if(projPDG==-14)
847    {
848      nuanu=false;                       // Anti-neutrino
849      CSmanager=G4QANuMuNuclearCrossSection::GetPointer();  // (anu,mu+) CC @@ open
850      CSmanager=G4QANuANuNuclearCrossSection::GetPointer(); // (anu,anu) NC @@ open
851      scatPDG=-13;                       // secondary scattered mu+
852    }
853    else if(projPDG==12)
854    {
855      CSmanager=G4QNuENuclearCrossSection::GetPointer(); // @@ open (only CC is changed)
856      scatPDG=11;                        // secondary scattered e-
857    }
858    else if(projPDG==-12)
859    {
860      nuanu=false;                       // anti-neutrino
861      CSmanager=G4QANuENuclearCrossSection::GetPointer();   // (anu,e+) CC @@ open
862      CSmanager=G4QANuANuNuclearCrossSection::GetPointer(); // (anu,anu) NC @@ open
863      scatPDG=-11;                       // secondary scattered e+
864    }
865    // @@ Probably this is not necessary any more
866    G4double xSec1=CSmanager->GetCrossSection(false,Momentum,Z,N,projPDG); //Recalculate XS
867    G4double xSec2=CSmanager2->GetCrossSection(false,Momentum,Z,N,projPDG);//Recalculate XS
868    G4double xSec=xSec1+xSec2;
869    // @@ check a possibility to separate p, n, or alpha (!)
870    if(xSec <= 0.) // The cross-section = 0 -> Do Nothing
871    {
872      G4cerr<<"G4QCollision::PSDoIt:nuE="<<kinEnergy<<",X1="<<xSec1<<",X2="<<xSec2<<G4endl;
873      //Do Nothing Action insead of the reaction
874      aParticleChange.ProposeEnergy(kinEnergy);
875      aParticleChange.ProposeLocalEnergyDeposit(0.);
876      aParticleChange.ProposeMomentumDirection(dir);
877      aParticleChange.ProposeTrackStatus(fAlive);
878      return G4VDiscreteProcess::PostStepDoIt(track,step);
879    }
880    G4bool secnu=false;
881    if(xSec*G4UniformRand()>xSec1)               // recover neutrino/antineutrino
882                                {
883      if(scatPDG>0) scatPDG++;
884      else          scatPDG--;
885      secnu=true;
886    }
887    scat=true;                                   // event with changed scattered projectile
888    G4double totCS1 = CSmanager->GetLastTOTCS(); // the last total cross section1(isotope?)
889    G4double totCS2 = CSmanager2->GetLastTOTCS();// the last total cross section2(isotope?)
890    G4double totCS  = totCS1+totCS2;             // the last total cross section (isotope?)
891    if(std::fabs(xSec-totCS*millibarn)/xSec>.0001)
892          G4cout<<"-Warning-G4QCollision::PostStepDoIt: xS="<<xSec<<"# CS="<<totCS<<G4endl;
893    G4double qelCS1 = CSmanager->GetLastQELCS(); // the last quasi-elastic cross section1
894    G4double qelCS2 = CSmanager2->GetLastQELCS();// the last quasi-elastic cross section2
895    G4double qelCS  = qelCS1+qelCS2;             // the last quasi-elastic cross section
896    if(totCS - qelCS < 0.)                       // only at low energies
897    {
898      totCS  = qelCS;
899      totCS1 = qelCS1;
900      totCS2 = qelCS2;
901                                }
902    // make different definitions for neutrino and antineutrino
903    G4double mIN=mProt;                          // Just a prototype (for anu, Z=1, N=0)
904    G4double mOT=mNeut;
905    G4double OT=mlN2;
906    G4double mOT2=mNeut2;
907    G4double mlOT=fmlN;
908    G4double mlsOT=mlsN;
909    if(secnu)
910    {
911      if(am*G4UniformRand()>Z)                   // Neutron target
912      {
913        targPDG-=1;                              // subtract neutron
914        projPDG=2112;                            // neutron is going out
915        mIN =mNeut;                     
916        OT  =mNeut2;
917        mOT2=mNeut2;
918        mlOT=0.;
919        mlsOT=mNeut2;
920      }
921      else
922      {
923        targPDG-=1000;                           // subtract neutron
924        projPDG=2212;                            // neutron is going out
925        mOT  =mProt;
926        OT   =mProt2;
927        mOT2 =mProt2;
928        mlOT =0.;
929        mlsOT=mProt2;
930      }
931      ml=0.;
932      ml2=0.;
933      mldM=0.;
934      mlD2=mPPi2;
935      G4QPDGCode targQPDG(targPDG);
936      G4double rM=targQPDG.GetMass();
937      mIN=tM-rM;                                 // bounded in-mass of the neutron
938      tM=rM;
939    }
940    else if(nuanu)
941    {
942      targPDG-=1;                                // Neutrino -> subtract neutron
943      G4QPDGCode targQPDG(targPDG);
944      G4double rM=targQPDG.GetMass();
945      mIN=tM-rM;                                 // bounded in-mass of the neutron
946      tM=rM;
947      mOT=mProt;
948      OT=mlP2;
949      mOT2=mProt2;
950      mlOT=fmlP;
951      mlsOT=mlsP;
952      projPDG=2212;                              // proton is going out
953    }
954    else
955    {
956      if(Z>1||N>0)                               // Calculate the splitted mass
957                                                {
958        targPDG-=1000;                           // Anti-Neutrino -> subtract proton
959        G4QPDGCode targQPDG(targPDG);
960        G4double rM=targQPDG.GetMass();
961        mIN=tM-rM;                               // bounded in-mass of the proton
962        tM=rM;
963      }
964      else
965      {
966        targPDG=0;
967        mIN=tM;
968        tM=0.;
969      }
970      projPDG=2112;                              // neutron is going out
971    }
972    G4double s=mIN*(mIN+dKinE);                  // s=(M_cm)^2=m2+2mE (m=targetMass,E=E_nu)
973#ifdef debug
974                  G4cout<<"G4QCollision::PostStDoIt: s="<<s<<" >? OT="<<OT<<", mlD2="<<mlD2<<G4endl;
975#endif
976    if(s<=OT)                                    // *** Do nothing solution ***
977    {
978      //Do NothingToDo Action insead of the reaction (@@ Can we make it common?)
979      G4cout<<"G4QCollision::PostStepDoIt: tooSmallFinalMassOfCompound: DoNothing"<<G4endl;
980      aParticleChange.ProposeEnergy(kinEnergy);
981      aParticleChange.ProposeLocalEnergyDeposit(0.);
982      aParticleChange.ProposeMomentumDirection(dir);
983      aParticleChange.ProposeTrackStatus(fAlive);
984      return G4VDiscreteProcess::PostStepDoIt(track,step);
985    }
986#ifdef debug
987                G4cout<<"G4QCollision::PostStDoIt: Stop and kill the projectile neutrino"<<G4endl;
988#endif
989    aParticleChange.ProposeEnergy(0.);
990    aParticleChange.ProposeTrackStatus(fStopAndKill); // the initial neutrino is killed
991    // There is no way back from here
992
993    if ( ((secnu || !nuanu || N) && totCS*G4UniformRand() < qelCS) || s < mlD2 ) 
994    {   // Quasi-Elastic interaction
995      G4double Q2=0.;                           // Simulate transferred momentum, in MeV^2
996      if(secnu) Q2=CSmanager2->GetQEL_ExchangeQ2();
997      else      Q2=CSmanager->GetQEL_ExchangeQ2();
998#ifdef debug
999                  G4cout<<"G4QCollision::PostStDoIt:QuasiEl(nu="<<secnu<<"),s="<<s<<",Q2="<<Q2<<G4endl;
1000#endif
1001      //G4double ds=s+s;                          // doubled s
1002      G4double sqs=std::sqrt(s);                // M_cm
1003      G4double dsqs=sqs+sqs;                    // 2*M_cm
1004      G4double pi=(s-mIN*mIN)/dsqs;             // initial momentum in CMS (checked MK)
1005      G4double dpi=pi+pi;                       // doubled initial momentum in CMS
1006      G4double sd=s-mlsOT;                      // s-ml2-mOT2 (mlsOT=m^2_neut+m^2_lept)
1007      G4double qo2=(sd*sd-mlOT)/dsqs;           // squared momentum of secondaries in CMS
1008      G4double qo=std::sqrt(qo2);               // momentum of secondaries in CMS
1009      G4double cost=(dpi*std::sqrt(qo2+ml2)-Q2-ml2)/dpi/qo; // cos(theta) in CMS (chck MK)
1010      G4LorentzVector t4M(0.,0.,0.,mIN);        // 4mom of the effective target
1011      G4LorentzVector c4M=t4M+proj4M;           // 4mom of the compound system
1012      t4M.setT(mOT);                            // now it is 4mom of the outgoing nucleon
1013      scat4M=G4LorentzVector(0.,0.,0.,ml);      // 4mom of the scattered muon
1014      if(!G4QHadron(c4M).RelDecayIn2(scat4M, t4M, proj4M, cost, cost))
1015      {
1016        G4cerr<<"G4QCol::PSD:c4M="<<c4M<<sqs<<",mM="<<ml<<",tM="<<mOT<<",c="<<cost<<G4endl;
1017        throw G4QException("G4QCollision::HadronizeQuasm: Can't dec QE nu,lept Compound");
1018      }
1019      proj4M=t4M;                               // 4mom of the new projectile nucleon
1020    }
1021    else                                        // ***** Non Quasi Elastic interaction
1022    {
1023     
1024      if ( (secnu && projPDG == 2212) || (!secnu && projPDG == 2112) ) {
1025        targPDG+=1;    // Recover target PDG,
1026      } else if ( (secnu && projPDG == 2112) || (!secnu && projPDG == 2212) ) { 
1027        targPDG+=1000; // if not quasiEl
1028      }
1029      G4double Q2=0;                            // Simulate transferred momentum, in MeV^2
1030      if(secnu) Q2=CSmanager->GetNQE_ExchangeQ2();
1031      else      Q2=CSmanager2->GetNQE_ExchangeQ2();
1032#ifdef debug
1033                  G4cout<<"G4QColl::PStDoIt: MultiPeriferal s="<<s<<",Q2="<<Q2<<",T="<<targPDG<<G4endl;
1034#endif
1035      if(secnu) projPDG=CSmanager2->GetExchangePDGCode();// PDG Code of the effective gamma
1036      else      projPDG=CSmanager->GetExchangePDGCode(); // PDG Code of the effective pion
1037      //@@ Temporary made only for direct interaction and for N=3 (good for small Q2)
1038      //@@ inFuture use N=GetNPartons and directFraction=GetDirectPart, @@ W2...
1039      G4double r=G4UniformRand();
1040      G4double r1=0.5;                          // (1-x)
1041      if(r<0.5)      r1=std::sqrt(r+r)*(.5+.1579*(r-.5));
1042      else if(r>0.5) r1=1.-std::sqrt(2.-r-r)*(.5+.1579*(.5-r));
1043      G4double xn=1.-mldM/Momentum;             // Normalization of (1-x) [x>mldM/Mom]
1044      G4double x1=xn*r1;                        // (1-x)
1045      G4double x=1.-x1;                         // x=2k/M
1046      //G4double W2=(hdM2+Q2/x)*x1;               // W2 candidate
1047      G4double mx=hdM*x;                        // Part of the target to interact with
1048      G4double we=Q2/(mx+mx);                   // transfered energy
1049      if(we>=kinEnergy-ml-.001) we=kinEnergy-ml-.0001; // safety to avoid nan=sqrt(neg)
1050      G4double pot=kinEnergy-we;                // energy of the secondary lepton
1051      G4double mlQ2=ml2+Q2;
1052      G4double cost=(pot-mlQ2/dKinE)/std::sqrt(pot*pot-ml2); // LS cos(theta)
1053      if(std::fabs(cost)>1)
1054      {
1055#ifdef debug
1056                    G4cout<<"*G4QCollision::PostStDoIt: cost="<<cost<<", Q2="<<Q2<<", nu="<<we<<", mx="
1057              <<mx<<", pot="<<pot<<", 2KE="<<dKinE<<G4endl;
1058#endif
1059        if(cost>1.) cost=1.;
1060        else        cost=-1.;
1061        pot=mlQ2/dKinE+dKinE*ml2/mlQ2;          // extreme output momentum
1062      }
1063      G4double lEn=std::sqrt(pot*pot+ml2);      // Lepton energy
1064      G4double lPl=pot*cost;                    // Lepton longitudinal momentum
1065      G4double lPt=pot*std::sqrt(1.-cost*cost); // Lepton transverse momentum
1066      std::pair<G4double,G4double> d2d=Random2DDirection(); // Randomize phi
1067      G4double lPx=lPt*d2d.first;
1068      G4double lPy=lPt*d2d.second;
1069      G4ThreeVector vdir=proj4M.vect();         // 3D momentum of the projectile
1070      G4ThreeVector vz= vdir.unit();            // Ort in the direction of the projectile
1071      G4ThreeVector vv= vz.orthogonal();        // Not normed orthogonal vector (!)
1072      G4ThreeVector vx= vv.unit();              // First ort orthogonal to the direction
1073      G4ThreeVector vy= vz.cross(vx);           // Second ort orthoganal to the direction
1074      G4ThreeVector lP= lPl*vz+lPx*vx+lPy*vy;   // 3D momentum of the scattered lepton
1075      scat4M=G4LorentzVector(lP,lEn);           // 4mom of the scattered lepton
1076      proj4M-=scat4M;                           // 4mom of the W/Z (effective pion/gamma)
1077#ifdef debug
1078                  G4cout<<"G4QCollision::PostStDoIt: proj4M="<<proj4M<<", ml="<<ml<<G4endl;
1079#endif
1080      // Check that the en/mom transfer is possible, if not -> elastic
1081      G4int fintPDG=targPDG;                    // Prototype for the compound nucleus
1082      if(!secnu)
1083      {
1084        if(projPDG<0) fintPDG-= 999;
1085        else          fintPDG+= 999;
1086      }
1087      G4double fM=G4QPDGCode(fintPDG).GetMass();// compound nucleus Mass (MeV)
1088      G4double fM2=fM*fM;
1089      G4LorentzVector tg4M=G4LorentzVector(0.,0.,0.,tgM);
1090      G4LorentzVector c4M=tg4M+proj4M;
1091#ifdef debug
1092      G4cout<<"G4QCol::PSDI:fM2="<<fM2<<" <? mc4M="<<c4M.m2()<<",dM="<<fM-tgM<<G4endl;
1093#endif
1094      if(fM2>=c4M.m2())                         // Elastic scattering should be done
1095      {
1096        G4LorentzVector tot4M=tg4M+proj4M+scat4M; // recover the total 4-momentum
1097        s=tot4M.m2();
1098        G4double fs=s-fM2-ml2;
1099        G4double fMl=fM2*ml2;
1100        G4double hQ2max=(fs*fs/2-fMl-fMl)/s;    // Maximum possible Q2/2
1101        G4double cost=1.-Q2/hQ2max;             // cos(theta) in CMS (use MultProd Q2)
1102#ifdef debug
1103        G4cout<<"G4QC::PSDI:ct="<<cost<<",Q2="<<Q2<<",hQ2="<<hQ2max<<",4M="<<tot4M<<G4endl;
1104#endif
1105        G4double acost=std::fabs(cost);
1106        if(acost>1.)
1107        {
1108          if(acost>1.001) G4cout<<"-Warning-G4QCollision::PostStDoIt: cost="<<cost<<G4endl;
1109          if     (cost> 1.) cost= 1.;
1110          else if(cost<-1.) cost=-1.;
1111        }
1112        G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,fM); // 4mom of the recoilNucleus
1113        scat4M=G4LorentzVector(0.,0.,0.,ml); // 4mom of the scatteredLepton
1114        G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-ml)*.01);
1115        if(!G4QHadron(tot4M).RelDecayIn2(scat4M, reco4M, dir4M, cost, cost))
1116        {
1117          G4cerr<<"G4QC::PSDI:t4M="<<tot4M<<",lM="<<ml<<",rM="<<fM<<",cost="<<cost<<G4endl;
1118          //G4Exception("G4QCollision::PostStepDoIt:","027",FatalException,"ElasticDecay");
1119        }
1120#ifdef debug
1121        G4cout<<"G4QCol::PStDoI:l4M="<<scat4M<<"+r4M="<<reco4M<<"="<<scat4M+reco4M<<G4endl;
1122#endif
1123        // ----------------------------------------------------
1124        G4ParticleDefinition* theDefinition=0; // Prototype of a particle for E-Secondaries
1125        // Fill scattered lepton
1126        if     (scatPDG==-11) theDefinition = G4Positron::Positron();
1127        else if(scatPDG== 11) theDefinition = G4Electron::Electron();
1128        else if(scatPDG== 13) theDefinition = G4MuonMinus::MuonMinus();
1129        else if(scatPDG==-13) theDefinition = G4MuonPlus::MuonPlus();
1130        //else if(scatPDG== 15) theDefinition = G4TauMinus::TauMinus();
1131        //else if(scatPDG==-15) theDefinition = G4TauPlus::TauPlus();
1132        if     (scatPDG==-12) theDefinition = G4AntiNeutrinoE::AntiNeutrinoE();
1133        else if(scatPDG== 12) theDefinition = G4NeutrinoE::NeutrinoE();
1134        else if(scatPDG== 14) theDefinition = G4NeutrinoMu::NeutrinoMu();
1135        else if(scatPDG==-14) theDefinition = G4AntiNeutrinoMu::AntiNeutrinoMu();
1136        //else if(scatPDG== 16) theDefinition = G4NeutrinoTau::NeutrinoTau();
1137        //else if(scatPDG==-16) theDefinition = G4AntiNeutrinoTau::AntiNeutrinoTau();
1138        else  G4cout<<"-Warning-G4QCollision::PostStDoIt: UnknownLepton="<<scatPDG<<G4endl;
1139        G4DynamicParticle* theScL = new G4DynamicParticle(theDefinition,scat4M);
1140        G4Track* scatLep = new G4Track(theScL, localtime, position ); //    scattered
1141        scatLep->SetWeight(weight);                                   //    weighted
1142        scatLep->SetTouchableHandle(trTouchable);                     //    residual
1143        aParticleChange.AddSecondary(scatLep);                        //    lepton
1144        // Fill residual nucleus
1145        if     (fintPDG==90000001) theDefinition = G4Neutron::Neutron(); // neutron
1146        else if(fintPDG==90001000) theDefinition = G4Proton::Proton();   // proton
1147        else                                                             // ion
1148        {
1149          G4int fm=static_cast<G4int>(fintPDG/1000000);               // Strange part
1150          G4int ZN=fintPDG-1000000*fm;
1151          G4int rZ=static_cast<G4int>(ZN/1000);
1152          G4int rA=ZN-999*rZ;
1153          theDefinition = G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,0,rZ);
1154        }
1155        G4DynamicParticle* theReN = new G4DynamicParticle(theDefinition,reco4M);
1156        G4Track* scatReN = new G4Track(theReN, localtime, position ); //    scattered
1157        scatReN->SetWeight(weight);                                   //    weighted
1158        scatReN->SetTouchableHandle(trTouchable);                     //    residual
1159        aParticleChange.AddSecondary(scatReN);                        //    nucleus
1160        return G4VDiscreteProcess::PostStepDoIt(track, step);
1161      }
1162    }
1163
1164  //
1165  // quasi-elastic for p+A(Z,N)
1166  //
1167  } else if (aProjPDG == 2212 && Z > 0 && N > 0) {
1168    //else if(2>3)
1169
1170    G4QuasiFreeRatios* qfMan=G4QuasiFreeRatios::GetPointer();
1171    std::pair<G4double,G4double> fief=qfMan->GetRatios(momentum, aProjPDG, Z, N);
1172    G4double qepart=fief.first*fief.second;
1173#ifdef qedebug
1174    G4cout<<"G4QCol::PSD:QE[p("<<proj4M<<")+(Z="<<Z<<",N="<<N<<",)="<<qepart<<G4endl;
1175#endif
1176    if(G4UniformRand()<qepart) // Make a quasi free scattering (out:A-1,h,N) @@ KinLim
1177    {
1178      // First decay a nucleus in a nucleon and a residual (A-1) nucleus
1179      G4double dmom=91.; // Fermi momentum (proto default for a deuteron)
1180      if(Z>1||N>1) dmom=286.2*std::pow(-std::log(G4UniformRand()),third);// p_max=250 MeV/c
1181
1182      // Calculate cluster probabilities (n,p,d,t,he3,he4 now only, can use UpdateClusters)
1183      const G4int lCl=3; // The last clProb[lCl]==1. by definition, MUST be increasing
1184      G4double clProb[lCl]={0.6,0.7,0.8}; // N/P,D,t/He3,Alpha, integrated prob for .6,.1,.1,.2
1185      G4double base=1.;  // Base for randomization (can be reduced by totZ & totN)
1186      G4int max=lCl;   // Number of boundaries (can be reduced by totZ & totN)
1187
1188      // Take into account that at least one nucleon must be left !
1189      // Change max-- to --max - DHW 05/08
1190      G4int A = Z + N;       // Baryon number of the nucleus
1191      if (Z<2 || N<2 || A<5) base = clProb[--max]; // Alpha cluster is impossible
1192      if ( (Z > 1 && N < 2) || (Z < 2 && N > 1) ) 
1193        base=(clProb[max]+clProb[max-1])/2; // t or He3 is impossible
1194
1195      if ( (Z < 2 && N < 2) || A < 4) base=clProb[--max]; // Both He3 and t clusters are impossible
1196
1197      if(A<3)           base=clProb[--max]; // Deuteron cluster is impossible
1198      G4int cln=0;                          // Cluster#0 (Default for the selected nucleon)
1199      if(max)                               // Not only nucleons are possible
1200      //if(2>3)
1201      {
1202        G4double ran=base*G4UniformRand();  // Base can be reduced
1203        G4int ic=0;                         // Start from the smallest cluster boundary
1204        while(ic<max) if(ran>clProb[ic++]) cln=ic;
1205      }
1206      G4ParticleDefinition* theDefinition;  // Prototype for qfNucleon
1207      G4bool cp1 = cln+2==A;                // A=ClusterBN+1 condition
1208      // Values to be defined in the following IF/ELSE
1209      G4LorentzVector r4M(0.,0.,0.,0.);     // Prototype of 4mom of the residual nucleus
1210      G4LorentzVector n4M(0.,0.,0.,0.);     // Prototype of 4mom of the quasi-cluster
1211      G4int nPDG=90000001;                  // Prototype for quasi-cluster mass calculation
1212      G4int restPDG=targPDG;                // Prototype should be reduced by quasi-cluster
1213      G4int rA=Z+N-1;                       // Prototype for the residualNucl definition
1214      G4int rZ=Z;                           // residZ: OK for the quasi-free neutron
1215      G4int nA=1;                           // Prototype for the quasi-cluster definition
1216      G4int nZ=0;                           // nA=1,nZ=0: OK for the quasi-free neutron
1217      G4double qM=mNeut;                    // Free mass of the quasi-free cluster
1218      if(!cln || cp1)                       // Split in nucleon + (A-1) with Fermi momentum
1219      {
1220        G4int nln=0;
1221        if(cln==2) nln=1;                         // @@ only for cp1: t/He3 choice from A=4
1222        // mass(A)=tM. Calculate masses of A-1 (rM) and mN (mNeut or mProt bounded mass)
1223        if ( ((!cln || cln == 2) && G4UniformRand()*(A-cln) > (N-nln)) || 
1224             ((cln == 3 || cln == 1) && Z > N) )
1225        {
1226          nPDG=90001000;                          // Update quasi-free nucleon PDGCode to P
1227          nZ=1;                                   // Change charge of the quasiFree nucleon
1228          qM=mProt;                               // Update quasi-free nucleon mass
1229          rZ--;                                   // Reduce the residual Z
1230          restPDG-=1000;                          // Reduce the residual PDGCode
1231        }
1232        else restPDG--;
1233        G4LorentzVector t4M(0.,0.,0.,tM);         // 4m of the target nucleus to be decayed
1234        G4double rM=G4QPDGCode(restPDG).GetMass();// Mass of the residual nucleus
1235        r4M=G4LorentzVector(0.,0.,0.,rM);         // 4mom of the residual nucleus
1236        G4double rM2=rM*rM;
1237        G4double nM=std::sqrt(rM2+tM*tM-(tM+tM)*std::sqrt(rM2+dmom*dmom));// M of q-nucleon
1238        n4M=G4LorentzVector(0.,0.,0.,nM);         // 4mom of the quasi-nucleon
1239#ifdef qedebug
1240                      G4cout<<"G4QCollis::PStDoIt:QE,p="<<dmom<<",tM="<<tM<<",R="<<rM<<",N="<<nM<<G4endl;
1241#endif
1242        if(!G4QHadron(t4M).DecayIn2(r4M, n4M))
1243        {
1244          G4cerr<<"G4QCol::PostStDoIt: M="<<tM<<"<rM="<<rM<<"+nM="<<nM<<"="<<rM+nM<<G4endl;
1245          throw G4QException("G4QCollision::HadronizeQuasm:Can'tDec totNuc->QENuc+ResNuc");
1246        }
1247#ifdef qedebug
1248                      G4cout<<"G4QCol::PStDoIt:QE-N,RA="<<r4M.rho()<<r4M<<",QN="<<n4M.rho()<<n4M<<G4endl;
1249#endif
1250        if(cp1 && cln)                           // Quasi-cluster case: swap the output
1251        {
1252          qM=rM;                                 // Scattering will be made on a cluster
1253          nln=nPDG;
1254          nPDG=restPDG;
1255          restPDG=nln;
1256          t4M=n4M;
1257          n4M=r4M;
1258          r4M=t4M;
1259          nln=nZ;
1260          nZ=rZ;
1261          rZ=nln;
1262          nln=nA;
1263          nA=rA;
1264          rA=nln;
1265        }
1266      }
1267      else // Split a cluster (w or w/o "Fermi motion" and "Fermi decay")
1268      {
1269        if(cln==1)
1270        {
1271          nPDG=90001001;                  // Deuteron
1272          qM=mDeut;
1273          nA=2;
1274          nZ=1;
1275          restPDG-=1001;
1276        }
1277        else if(cln==2)
1278        {
1279          nA=3;
1280                                                                                if(G4UniformRand()*(A-2)>(N-1)) // He3
1281          {
1282            nPDG=90002001;
1283            qM=mHel3;
1284            nZ=2;
1285            restPDG-=2001;
1286          }
1287          else                            // tritium
1288          {
1289            nPDG=90001002;
1290            qM=mTrit;
1291            nZ=1;
1292            restPDG-=1002;
1293          }
1294        }
1295        else
1296        {
1297          nPDG=90002002;                  // Alpha
1298          qM=mAlph;
1299          nA=4;
1300          nZ=2;
1301          restPDG-=2002;
1302        }
1303        rA=A-nA;
1304        rZ=Z-nZ;
1305        // This is a simple case of cluster at rest
1306        //G4double rM=G4QPDGCode(restPDG).GetMass();// Mass of the residual nucleus
1307        //r4M=G4LorentzVector(0.,0.,0.,rM);         // 4mom of the residual nucleus
1308        //n4M=G4LorentzVector(0.,0.,0.,tM-rM);      // 4mom of the quasi-free cluster
1309        // --- End of the "simple case of cluster at rest"
1310        // Make a fake quasi-Fermi distribution for clusters (clusters are not at rest)
1311        G4LorentzVector t4M(0.,0.,0.,tM);         // 4m of the target nucleus to be decayed
1312        G4double rM=G4QPDGCode(restPDG).GetMass();// Mass of the residual nucleus
1313        r4M=G4LorentzVector(0.,0.,0.,rM);         // 4mom of the residual nucleus
1314        G4double rM2=rM*rM;
1315        G4double nM=std::sqrt(rM2+tM*tM-(tM+tM)*std::sqrt(rM2+dmom*dmom));// M of q-cluster
1316        n4M=G4LorentzVector(0.,0.,0.,nM);         // 4mom of the quasi-nucleon
1317#ifdef qedebug
1318                      G4cout<<"G4QCollis::PStDoIt:QEC,p="<<dmom<<",T="<<tM<<",R="<<rM<<",N="<<nM<<G4endl;
1319#endif
1320        if(!G4QHadron(t4M).DecayIn2(r4M, n4M))
1321        {
1322          G4cerr<<"G4QCol::PostStDoIt: M="<<tM<<"<rM="<<rM<<"+cM="<<nM<<"="<<rM+nM<<G4endl;
1323          throw G4QException("G4QCollision::HadronizeQuasm:Can'tDec totNuc->QEClu+ResNuc");
1324        }
1325        // --- End of the moving cluster implementation ---
1326#ifdef qedebug
1327                      G4cout<<"G4QCol::PStDoIt:QEC,RN="<<r4M.rho()<<r4M<<",QCl="<<n4M.rho()<<n4M<<G4endl;
1328#endif
1329      }
1330      G4LorentzVector s4M=n4M+proj4M;             // Tot 4-momentum for scattering
1331      G4double prjM2 = proj4M.m2();
1332      G4double prjM = std::sqrt(prjM2);           // @@ Get from pPDG (?)
1333      G4double minM = prjM+qM;                    // Min mass sum for the final products
1334      G4double cmM2 =s4M.m2();
1335      if(cmM2>minM*minM)
1336      {
1337#ifdef qedebug
1338                      G4cout<<"G4QCol::PStDoIt:***Enter***,cmM2="<<cmM2<<" > minM2="<<minM*minM<<G4endl;
1339#endif
1340        // Estimate and randomize charge-exchange with quasi-free cluster
1341        G4bool chex=false;                        // Flag of the charge exchange scattering
1342        G4ParticleDefinition* projpt=G4Proton::Proton(); // Prototype, only for chex=true
1343        //if(cln&&!cp1 &&(projPDG==2212&&rA>rZ || projPDG==2112&&rZ>1))// @@ Use proj chex
1344               if(2>3)
1345        {
1346#ifdef qedebug
1347                        G4cout<<"G4QCol::PStDoIt:-Enter,P="<<projPDG<<",cln="<<cln<<",cp1="<<cp1<<G4endl;
1348#endif
1349          G4double tprM=mProt;
1350          G4double tprM2=mProt2;
1351          G4int tprPDG=2212;
1352          G4int tresPDG=restPDG+999;
1353          if(projPDG==2212)
1354          {
1355            projpt=G4Neutron::Neutron();
1356            tprM=mNeut;
1357            tprM2=mNeut2;
1358            tprPDG=2112;
1359            tresPDG=restPDG-999;
1360          }
1361          minM=tprM+qM;
1362          G4double efE=(cmM2-tprM2-qM*qM)/(qM+qM);
1363          G4double efP=std::sqrt(efE*efE-tprM2);
1364          G4double chl=qfMan->ChExElCoef(efP*MeV, nZ, nA-nZ, projPDG); // ChEx/Elast(pPDG!)
1365#ifdef qedebug
1366                        G4cout<<"G4QCol::PStDoIt:chl="<<chl<<",P="<<efP<<",nZ="<<nZ<<",nA="<<nA<<G4endl;
1367#endif
1368          if(chl>0.&&cmM2>minM*minM&&G4UniformRand()<chl/(1.+chl))     // minM is redefined
1369          {
1370            projPDG=tprPDG;
1371            prjM=tprM;
1372            G4double rM=G4QPDGCode(tresPDG).GetMass();// Mass of the residual nucleus
1373            r4M=G4LorentzVector(0.,0.,0.,rM);         // 4mom of the residual nucleus
1374            n4M=G4LorentzVector(0.,0.,0.,tM-rM);      // 4mom of the quasi-free cluster
1375            chex=true;                                // Confirm charge exchange scattering
1376          }
1377        }
1378        //
1379        std::pair<G4LorentzVector,G4LorentzVector> sctout=qfMan->Scatter(nPDG, n4M,
1380                                                                         projPDG, proj4M);
1381#ifdef qedebug
1382                      G4cout<<"G4QCollis::PStDoIt:QElS,proj="<<prjM<<sctout.second<<",qfCl="<<qM
1383              <<sctout.first<<",chex="<<chex<<",nA="<<nA<<",nZ="<<nZ<<G4endl;
1384#endif
1385        aParticleChange.ProposeLocalEnergyDeposit(0.); // Everything is in particles
1386        // @@ @@ @@ Coulomb barriers must be checked !! @@ @@ @@ Skip if not
1387        if(chex)               // ==> Projectile is changed: fill everything to secondaries
1388        {
1389          aParticleChange.ProposeEnergy(0.);                // @@ ??
1390          aParticleChange.ProposeTrackStatus(fStopAndKill); // projectile nucleon is killed
1391          aParticleChange.SetNumberOfSecondaries(3); 
1392          G4DynamicParticle* thePrH = new G4DynamicParticle(projpt,sctout.second);
1393          G4Track* scatPrH = new G4Track(thePrH, localtime, position ); // scattered & chex
1394          scatPrH->SetWeight(weight);                                   //    weighted
1395          scatPrH->SetTouchableHandle(trTouchable);                     //   projectile
1396          aParticleChange.AddSecondary(scatPrH);                        //     hadron
1397        }
1398        else               // ==> The leading particle is filled to the updated projectilee
1399        {
1400          aParticleChange.SetNumberOfSecondaries(2);        // @@ if proj=leading
1401          G4double ldT=(sctout.second).e()-prjM;            // kin Energy of scat project.
1402          aParticleChange.ProposeEnergy(ldT);               // Change the kin Energy
1403          G4ThreeVector ldV=(sctout.second).vect();         // Change momentum direction
1404          aParticleChange.ProposeMomentumDirection(ldV/ldV.mag());
1405          aParticleChange.ProposeTrackStatus(fAlive);
1406        }
1407        // ---------------------------------------------------------
1408        // Fill scattered quasi-free nucleon
1409        if     (nPDG==90000001) theDefinition = G4Neutron::Neutron();
1410        else if(nPDG==90001000) theDefinition = G4Proton::Proton();
1411        else theDefinition = G4ParticleTable::GetParticleTable()->FindIon(nZ,nA,0,nZ);//ion
1412        G4DynamicParticle* theQFN = new G4DynamicParticle(theDefinition,sctout.first);
1413        G4Track* scatQFN = new G4Track(theQFN, localtime, position ); //   scattered
1414        scatQFN->SetWeight(weight);                                   //    weighted
1415        scatQFN->SetTouchableHandle(trTouchable);                     //   quasi-free
1416        aParticleChange.AddSecondary(scatQFN);                        //  nucleon/cluster
1417        // ----------------------------------------------------
1418        // Fill residual nucleus
1419        if     (restPDG==90000001) theDefinition = G4Neutron::Neutron();
1420        else if(restPDG==90001000) theDefinition = G4Proton::Proton();
1421        else theDefinition = G4ParticleTable::GetParticleTable()->FindIon(rZ,rA,0,rZ);//ion
1422        G4DynamicParticle* theReN = new G4DynamicParticle(theDefinition,r4M);
1423        G4Track* scatReN = new G4Track(theReN, localtime, position ); //    scattered
1424        scatReN->SetWeight(weight);                                   //    weighted
1425        scatReN->SetTouchableHandle(trTouchable);                     //    residual
1426        aParticleChange.AddSecondary(scatReN);                        //    nucleus
1427        return G4VDiscreteProcess::PostStepDoIt(track, step);
1428      }
1429#ifdef qedebug
1430      else G4cout<<"G4QCol::PSD: OUT, M2="<<s4M.m2()<<"<"<<minM*minM<<", N="<<nPDG<<G4endl;
1431#endif
1432    }
1433  }  // end lepto-nuclear, neutrino-nuclear, proton quasi-elastic
1434
1435  EnMomConservation=proj4M+G4LorentzVector(0.,0.,0.,tM);    // Total 4-mom of the reaction
1436  if(absMom) EnMomConservation+=lead4M;         // Add E/M of leading System
1437#ifdef debug
1438  G4cout<<"G4QCollision::PostStDoIt:before St="<<aParticleChange.GetTrackStatus()<<G4endl;
1439#endif
1440
1441  // @@@@@@@@@@@@@@ Temporary for the testing purposes --- Begin
1442  //G4bool elF=false; // Flag of the ellastic scattering is "false" by default
1443  //G4double eWei=1.;
1444  // @@@@@@@@@@@@@@ Temporary for the testing purposes --- End 
1445#ifdef debug
1446  G4cout<<"G4QCollision::PostStepDoIt: projPDG="<<projPDG<<", targPDG="<<targPDG<<G4endl;
1447#endif
1448  G4QHadron* pH = new G4QHadron(projPDG,proj4M);                // ---> DELETED -->--  -+
1449  //if(momentum<1000.) // Condition for using G4QEnvironment (not G4QuasmonString)      |
1450  { //                                                                                  |
1451    G4QHadronVector projHV;                                 //                          |
1452    projHV.push_back(pH);                                   // DESTROYED over 2 lines-+ |
1453    G4QEnvironment* pan= new G4QEnvironment(projHV,targPDG);// ---> DELETED --->----+ | |
1454    std::for_each(projHV.begin(), projHV.end(), DeleteQHadron()); // <---<------<---+-+-+
1455    projHV.clear(); // <------------<---------------<-------------------<-----------+-+ .
1456#ifdef debug
1457    G4cout<<"G4QCol::PStDoIt:Proj="<<projPDG<<proj4M<<",Targ="<<targPDG<<G4endl; // |   .
1458#endif
1459    try                                                           //                |   .
1460    {                                                             //                |   .
1461      delete output;                                              //                |   .
1462      output = pan->Fragment();// DESTROYED in the end of the LOOP work space       |   .
1463    }                                                             //                |   .
1464    catch (G4QException& error)//                                                   |   .
1465    {                                                             //                |   .
1466             //#ifdef pdebug
1467      G4cerr<<"***G4QCollision::PostStepDoIt: G4QE Exception is catched"<<G4endl;// |   .
1468             //#endif
1469      G4Exception("G4QCollision::PostStepDoIt:","27",FatalException,"CHIPSCrash");//|   .
1470    }                                                             //                |   .
1471    delete pan;                              // Delete the Nuclear Environment <-<--+   .
1472  } //                                                                                  .
1473  //else             // Use G4QuasmonString                                             .
1474  //{ //                                                                                  ^
1475  //  G4QuasmonString* pan= new G4QuasmonString(pH,false,targPDG,false);//-> DELETED --+  |
1476  //  delete pH;                                                    // --------<-------+--+
1477#ifdef debug
1478  //  G4double mp=G4QPDGCode(projPDG).GetMass();   // Mass of the projectile particle  |
1479  //  G4cout<<"G4QCollision::PostStepDoIt: pPDG="<<projPDG<<", pM="<<mp<<G4endl; //    |
1480#endif
1481  //  //G4int tNH=0;                    // Prototype of the number of secondaries inOut|
1482  //  try                                                           //                 |
1483  //  {                                                             //                 |
1484  //              delete output;                                  //                   |
1485  //    output = pan->Fragment();// DESTROYED in the end of the LOOP work space        |
1486  //    // @@@@@@@@@@@@@@ Temporary for the testing purposes --- Begin                 |
1487  //    //tNH=pan->GetNOfHadrons();     // For the test purposes of the String         |
1488  //    //if(tNH==2)                    // At least 2 hadrons are in the Constr.Output |
1489  //    //{//                                                                          |
1490  //    //  elF=true;                   // Just put a flag for the ellastic Scattering |
1491  //    //  delete output;              // Delete a prototype of dummy G4QHadronVector |
1492  //    //  output = pan->GetHadrons(); // DESTROYED in the end of the LOOP work space |
1493  //    //}//                                                                          |
1494  //    //eWei=pan->GetWeight();        // Just an example for the weight of the event |
1495#ifdef debug
1496  //    //G4cout<<"=====>>G4QCollision::PostStepDoIt: elF="<<elF<<",n="<<tNH<<G4endl;//|
1497#endif
1498  //    // @@@@@@@@@@@@@@ Temporary for the testing purposes --- End                   |
1499  //  }                                                             //                 |
1500  //  catch (G4QException& error)//                                                    |
1501  //  {                                                             //                 |
1502  //    //#ifdef pdebug
1503  //    G4cerr<<"***G4QCollision::PostStepDoIt: GEN Exception is catched"<<G4endl; //  |
1504         //    //#endif
1505  //    G4Exception("G4QCollision::PostStDoIt:","27",FatalException,"QString Excep");//|
1506  //  }                                                             //                 |
1507  //  delete pan;                              // Delete the Nuclear Environment ---<--+
1508  //}
1509  // --- the scattered hadron with changed nature can be added here ---
1510  if(scat)
1511  {
1512    G4QHadron* scatHadron = new G4QHadron(scatPDG,scat4M);
1513    output->push_back(scatHadron);
1514  }
1515  G4int qNH=leadhs->size();
1516  if(absMom)
1517  {
1518    if(qNH) for(G4int iq=0; iq<qNH; iq++)
1519    {
1520      G4QHadron* loh=(*leadhs)[iq];   // Pointer to the output hadron
1521      output->push_back(loh);
1522    }
1523    delete leadhs;
1524  }
1525
1526
1527  // ------------- From here the secondaries are filled -------------------------
1528  G4int tNH = output->size();       // A#of hadrons in the output
1529  aParticleChange.SetNumberOfSecondaries(tNH); 
1530  // Now add nuclear fragments
1531#ifdef debug
1532  G4cout<<"G4QCollision::PostStepDoIt: "<<tNH<<" particles are generated"<<G4endl;
1533#endif
1534#ifdef ppdebug
1535  if(absMom)G4cout<<"G4QCollision::PostStepDoIt: t="<<tNH<<", q="<<qNH<<G4endl;
1536#endif
1537  G4int nOut=output->size();               // Real length of the output @@ Temporary
1538  if(tNH==1 && !scat)                      // @@ Temporary. Find out why it happened!
1539  {
1540    G4cout<<"-Warning-G4QCollision::PostStepDoIt: 1 secondary! absMom="<<absMom;
1541    if(absMom) G4cout<<", qNH="<<qNH;
1542    G4cout<<", PDG0="<<(*output)[0]->GetPDGCode();
1543    G4cout<<G4endl;
1544    tNH=0;
1545    delete output->operator[](0);          // delete the creazy hadron
1546    output->pop_back();                    // clean up the output vector
1547  }
1548  if(tNH==2&&2!=nOut) G4cout<<"--Warning--G4QCollision::PostStepDoIt: 2 # "<<nOut<<G4endl;
1549  // Deal with ParticleChange final state interface to GEANT4 output of the process
1550  //if(tNH==2) for(i=0; i<tNH; i++)          // @@ Temporary tNH==2 instead of just tNH
1551  if(tNH) for(i=0; i<tNH; i++)             // @@ Temporary tNH==2 instead of just tNH
1552  {
1553    // Note that one still has to take care of Hypernuclei (with Lambda or Sigma inside)
1554    // Hypernucleus mass calculation and ion-table interface upgrade => work for Hisaya @@
1555    // The decau process for hypernuclei must be developed in GEANT4 (change CHIPS body)
1556    G4QHadron* hadr=(*output)[i];          // Pointer to the output hadron   
1557    G4int PDGCode = hadr->GetPDGCode();
1558    G4int nFrag   = hadr->GetNFragments();
1559#ifdef pdebug
1560    G4cout<<"G4QCollision::PostStepDoIt: H#"<<i<<",PDG="<<PDGCode<<",nF="<<nFrag
1561          <<", 4Mom="<<hadr->Get4Momentum()<<G4endl;
1562#endif
1563    if(nFrag)                              // Skip intermediate (decayed) hadrons
1564    {
1565#ifdef debug
1566             G4cout<<"G4QCollision::PostStepDoIt: Intermediate particle is found i="<<i<<G4endl;
1567#endif
1568      delete hadr;
1569      continue;
1570    }
1571    G4DynamicParticle* theSec = new G4DynamicParticle;
1572    G4ParticleDefinition* theDefinition;
1573    if     (PDGCode==90000001) theDefinition = G4Neutron::Neutron();
1574    else if(PDGCode==90001000) theDefinition = G4Proton::Proton();//While it can be in ions
1575    else if(PDGCode==91000000) theDefinition = G4Lambda::Lambda();
1576    else if(PDGCode==311 || PDGCode==-311)
1577    {
1578      if(G4UniformRand()>.5) theDefinition = G4KaonZeroLong::KaonZeroLong();   // K_L
1579                                                else                   theDefinition = G4KaonZeroShort::KaonZeroShort(); // K_S
1580    }
1581    else if(PDGCode==91000999) theDefinition = G4SigmaPlus::SigmaPlus();
1582    else if(PDGCode==90999001) theDefinition = G4SigmaMinus::SigmaMinus();
1583    else if(PDGCode==91999000) theDefinition = G4XiMinus::XiMinus();
1584    else if(PDGCode==91999999) theDefinition = G4XiZero::XiZero();
1585    else if(PDGCode==92998999) theDefinition = G4OmegaMinus::OmegaMinus();
1586           else if(PDGCode >80000000) // Defines hypernuclei as normal nuclei (N=N+S Correction!)
1587    {
1588      G4int aZ = hadr->GetCharge();
1589      G4int aA = hadr->GetBaryonNumber();
1590#ifdef pdebug
1591                                                G4cout<<"G4QCollision::PostStepDoIt:Ion Z="<<aZ<<", A="<<aA<<G4endl;
1592#endif
1593      theDefinition = G4ParticleTable::GetParticleTable()->FindIon(aZ,aA,0,aZ);
1594    }
1595    //else theDefinition = G4ParticleTable::GetParticleTable()->FindParticle(PDGCode);
1596    else
1597    {
1598#ifdef pdebug
1599                                                G4cout<<"G4QCollision::PostStepDoIt:Define particle with PDG="<<PDGCode<<G4endl;
1600#endif
1601      theDefinition = G4QPDGToG4Particle::Get()->GetParticleDefinition(PDGCode);
1602#ifdef pdebug
1603                                                G4cout<<"G4QCollision::PostStepDoIt:AfterParticleDefinition PDG="<<PDGCode<<G4endl;
1604#endif
1605    }
1606    if(!theDefinition)
1607    {
1608#ifdef debug
1609      G4cout<<"---Warning---G4QCollision::PostStepDoIt: drop PDG="<<PDGCode<<G4endl;
1610#endif
1611      delete hadr;
1612      continue;
1613    }
1614#ifdef pdebug
1615    G4cout<<"G4QCollision::PostStepDoIt:Name="<<theDefinition->GetParticleName()<<G4endl;
1616#endif
1617    theSec->SetDefinition(theDefinition);
1618    G4LorentzVector h4M=hadr->Get4Momentum();
1619    EnMomConservation-=h4M;
1620#ifdef tdebug
1621    G4cout<<"G4QCollis::PSDI:"<<i<<","<<PDGCode<<h4M<<h4M.m()<<EnMomConservation<<G4endl;
1622#endif
1623#ifdef debug
1624    G4cout<<"G4QCollision::PostStepDoIt:#"<<i<<",PDG="<<PDGCode<<",4M="<<h4M<<G4endl;
1625#endif
1626    theSec->Set4Momentum(h4M); //                                                         ^
1627    delete hadr; // <-----<-----------<-------------<---------------------<---------<-----+
1628#ifdef debug
1629    G4ThreeVector curD=theSec->GetMomentumDirection();              //                    ^
1630    G4double curM=theSec->GetMass();                                //                    |
1631    G4double curE=theSec->GetKineticEnergy()+curM;                  //                    ^
1632    G4cout<<"G4QCollis::PSDoIt:p="<<curD<<curD.mag()<<",e="<<curE<<",m="<<curM<<G4endl;// |
1633#endif
1634    G4Track* aNewTrack = new G4Track(theSec, localtime, position ); //                    ^
1635    aNewTrack->SetWeight(weight);                                   //    weighted        |
1636    aNewTrack->SetTouchableHandle(trTouchable);                     //                    |
1637    aParticleChange.AddSecondary( aNewTrack );                      //                    |
1638#ifdef debug
1639    G4cout<<"G4QCollision::PostStepDoIt:#"<<i<<" is done."<<G4endl; //                    |
1640#endif
1641  } //                                                                                    |
1642  delete output; // instances of the G4QHadrons from the output are already deleted above +
1643#ifdef debug
1644                G4cout<<"G4QCollision::PostStDoIt: after St="<<aParticleChange.GetTrackStatus()<<G4endl;
1645#endif
1646  if(aProjPDG!=11 && aProjPDG!=13 && aProjPDG!=15)
1647    aParticleChange.ProposeTrackStatus(fStopAndKill);        // Kill the absorbed particle
1648  //return &aParticleChange;                               // This is not enough (ClearILL)
1649#ifdef pdebug
1650                  G4cout<<"G4QCollision::PostStepDoIt: E="<<aParticleChange.GetEnergy()
1651          <<", d="<<*aParticleChange.GetMomentumDirection()<<G4endl;
1652#endif
1653#ifdef debug
1654                G4cout<<"G4QCollision::PostStepDoIt:*** PostStepDoIt is done ***, P="<<aProjPDG<<", St="
1655        <<aParticleChange.GetTrackStatus()<<G4endl;
1656#endif
1657  return G4VDiscreteProcess::PostStepDoIt(track, step);
1658}
1659
1660std::pair<G4double,G4double> G4QCollision::Random2DDirection()
1661{
1662  G4double sp=0;              // sin(phi)
1663  G4double cp=1.;             // cos(phi)
1664  G4double r2=2.;             // to enter the loop
1665  while(r2>1. || r2<.0001)    // pi/4 efficiency
1666  {
1667    G4double s=G4UniformRand();
1668    G4double c=G4UniformRand();
1669    sp=1.-s-s;
1670    cp=1.-c-c;
1671    r2=sp*sp+cp*cp;
1672  }
1673  G4double norm=std::sqrt(r2);
1674  return std::make_pair(sp/norm,cp/norm);
1675}
Note: See TracBrowser for help on using the repository browser.