source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/interface/src/G4QGluonString.cc @ 1055

Last change on this file since 1055 was 1055, checked in by garnier, 15 years ago

maj sur la beta de geant 4.9.3
File size: 29.9 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: G4QGluonString.cc,v 1.4 2009/02/23 09:49:24 mkossov Exp $
27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
28//
29//      ---------------- G4QGluonString class -----------------
30//                 by Mikhail Kossov, December 2003.
31// G4QGluonString class of the CHIPS Simulation Branch in GEANT4
32// ---------------------------------------------------------------
33// Short description: CHIPS object for the Quark-Gluon String
34// ---------------------------------------------------------------
35
36//#define debug
37//#define pdebug
38//#define ppdebug
39
40#include "G4QGluonString.hh"
41
42// Initialization of static vectors
43std::vector<G4int> G4QGluonString::ElementZ;       // Z of the element(i) in theLastCalc
44std::vector<G4double> G4QGluonString::ElProbInMat; // SumProbabilityElements in Material
45std::vector<std::vector<G4int>*> G4QGluonString::ElIsoN; // N of isotope(j) of Element(i)
46std::vector<std::vector<G4double>*>G4QGluonString::IsoProbInEl;//SumProbabIsotopeInElementI
47
48G4QGluonString::G4QGluonString(const G4String& processName):
49 G4VDiscreteProcess(processName, fHadronic)
50{
51#ifdef debug
52  G4cout<<"G4QGluonString::Constructor is called"<<G4endl;
53#endif
54  if (verboseLevel>0) G4cout<<GetProcessName()<<" process is created by CHIPS"<<G4endl;
55  SetProcessSubType(fHadronInelastic);
56  G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World with 234 particles
57  G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio); // Clusterization param's
58  G4Quasmon::SetParameters(Temperature,SSin2Gluons,EtaEtaprime);  // Hadronic parameters
59  G4QEnvironment::SetParameters(SolidAngle); // SolAngle of pbar-A secondary mesons capture
60  //@@ Initialize here other parameters
61}
62
63G4bool   G4QGluonString::manualFlag=false; // If false then standard parameters are used
64G4double G4QGluonString::Temperature=180.; // Critical Temperature (sensitive at High En)
65G4double G4QGluonString::SSin2Gluons=0.3;  // Supression of s-quarks (in respect to u&d)
66G4double G4QGluonString::EtaEtaprime=0.3;  // Supression of eta mesons (gg->qq/3g->qq)
67G4double G4QGluonString::freeNuc=0.5;      // Percentage of free nucleons on the surface
68G4double G4QGluonString::freeDib=0.05;     // Percentage of free diBaryons on the surface
69G4double G4QGluonString::clustProb=5.;     // Nuclear clusterization parameter
70G4double G4QGluonString::mediRatio=10.;    // medium/vacuum hadronization ratio
71G4int    G4QGluonString::nPartCWorld=152;  // The#of particles initialized in CHIPS World
72G4double G4QGluonString::SolidAngle=0.5;   // Part of Solid Angle to capture (@@A-dep.)
73G4bool   G4QGluonString::EnergyFlux=false; // Flag for Energy Flux use (not MultyQuasmon)
74G4double G4QGluonString::PiPrThresh=141.4; // Pion Production Threshold for gammas
75G4double G4QGluonString::M2ShiftVir=20000.;// Shift for M2=-Q2=m_pi^2 of the virtualGamma
76G4double G4QGluonString::DiNuclMass=1880.; // DoubleNucleon Mass for VirtualNormalization
77
78void G4QGluonString::SetManual()   {manualFlag=true;}
79void G4QGluonString::SetStandard() {manualFlag=false;}
80
81// Fill the private parameters
82void G4QGluonString::SetParameters(G4double temper, G4double ssin2g, G4double etaetap,
83                                     G4double fN, G4double fD, G4double cP, G4double mR,
84                                     G4int nParCW, G4double solAn, G4bool efFlag,
85                                     G4double piThresh, G4double mpisq, G4double dinum)
86{//  =============================================================================
87  Temperature=temper;
88  SSin2Gluons=ssin2g;
89  EtaEtaprime=etaetap;
90  freeNuc=fN;
91  freeDib=fD;
92  clustProb=cP;
93  mediRatio=mR;
94  nPartCWorld = nParCW;
95  EnergyFlux=efFlag;
96  SolidAngle=solAn;
97  PiPrThresh=piThresh;
98  M2ShiftVir=mpisq;
99  DiNuclMass=dinum;
100  G4QCHIPSWorld::Get()->GetParticles(nPartCWorld); // Create CHIPS World with 234 particles
101  G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio); // Clusterization param's
102  G4Quasmon::SetParameters(Temperature,SSin2Gluons,EtaEtaprime);  // Hadronic parameters
103  G4QEnvironment::SetParameters(SolidAngle); // SolAngle of pbar-A secondary mesons capture
104}
105
106// Destructor
107
108G4QGluonString::~G4QGluonString() {}
109
110// Internal E/P conservation 4-Mom & the selected number of neutrons in the Element
111
112G4LorentzVector G4QGluonString::GetEnegryMomentumConservation()
113{
114  return EnMomConservation;
115}
116
117G4int G4QGluonString::GetNumberOfNeutronsInTarget()
118{
119  return nOfNeutrons;
120}
121
122// output of the function must be in units of length! L=1/sig_V,sig_V=SUM(n(j,i)*sig(j,i)),
123// where n(i,j) is a number of nuclei of the isotop j of the element i in V=1(lengtUnit^3)
124// ********** All CHIPS cross sections are calculated in the surface units ************
125G4double G4QGluonString::GetMeanFreePath(const G4Track& aTrack,
126                                         G4double, G4ForceCondition* Fc)
127{
128#ifdef debug
129  G4cout<<"G4QGluonString::GetMeanFreePath: Called Fc="<<*Fc<<G4endl;
130#endif
131  *Fc = NotForced;
132#ifdef debug
133  G4cout<<"G4QGluonString::GetMeanFreePath: Before GetDynPart"<<G4endl;
134#endif
135  const G4DynamicParticle* incidentParticle = aTrack.GetDynamicParticle();
136#ifdef debug
137  G4cout<<"G4QGluonString::GetMeanFreePath: Before GetDef"<<G4endl;
138#endif
139  G4ParticleDefinition* incidentParticleDefinition=incidentParticle->GetDefinition();
140  if( !IsApplicable(*incidentParticleDefinition))
141    G4cout<<"-W-G4QGluonString::GetMeanFreePath called for NotImplementedParticle"<<G4endl;
142  // Calculate the mean Cross Section for the set of Elements(*Isotopes) in the Material
143  G4double Momentum = incidentParticle->GetTotalMomentum(); // 3-momentum of the Particle
144#ifdef debug
145  G4cout<<"G4QCollis::GetMeanFreePath: BeforeGetMaterial"<<G4endl;
146#endif
147  const G4Material* material = aTrack.GetMaterial();        // Get the current material
148  const G4double* NOfNucPerVolume = material->GetVecNbOfAtomsPerVolume();
149  const G4ElementVector* theElementVector = material->GetElementVector();
150  G4int nE=material->GetNumberOfElements();
151#ifdef debug
152  G4cout<<"G4QGluonString::GetMeanFreePath:"<<nE<<" Elem's in theMaterial"<<G4endl;
153#endif
154  G4VQCrossSection* CSmanager=0;
155  G4int pPDG=0;
156  if(incidentParticleDefinition == G4Proton::Proton())
157  {
158    CSmanager=G4QProtonNuclearCrossSection::GetPointer();
159    pPDG=2212;
160  }                            //@@ Make cross-section mahnagers for other mesons & baryons
161  else
162  {
163    G4cerr<<"***G4QGluonString::GetMeanFreePath: Particle isn't implemented"<<G4endl;
164    G4Exception("G4QGluonString::PostStepDoIt:","72",FatalException,"BadProjectile");
165  }
166 
167  G4QIsotope* Isotopes = G4QIsotope::Get(); // Pointer to the G4QIsotopes singleton
168  G4double sigma=0.;                        // Sums over elements for the material
169  G4int IPIE=IsoProbInEl.size();            // How many old elements?
170  if(IPIE) for(G4int ip=0; ip<IPIE; ++ip)   // Clean up the SumProb's of Isotopes (SPI)
171  {
172    std::vector<G4double>* SPI=IsoProbInEl[ip]; // Pointer to the SPI vector
173    SPI->clear();
174    delete SPI;
175    std::vector<G4int>* IsN=ElIsoN[ip];     // Pointer to the N vector
176    IsN->clear();
177    delete IsN;
178  }
179  ElProbInMat.clear();                      // Clean up the SumProb's of Elements (SPE)
180  ElementZ.clear();                         // Clear the body vector for Z of Elements
181  IsoProbInEl.clear();                      // Clear the body vector for SPI
182  ElIsoN.clear();                           // Clear the body vector for N of Isotopes
183  for(G4int i=0; i<nE; ++i)
184  {
185    G4Element* pElement=(*theElementVector)[i]; // Pointer to the current element
186    G4int Z = static_cast<G4int>(pElement->GetZ()); // Z of the Element
187    ElementZ.push_back(Z);                  // Remember Z of the Element
188    G4int isoSize=0;                        // The default for the isoVectorLength is 0
189    G4int indEl=0;                          // Index of non-trivial element or 0(default)
190    G4IsotopeVector* isoVector=pElement->GetIsotopeVector(); // Get the predefined IsoVect
191    if(isoVector) isoSize=isoVector->size();// Get size of the existing isotopeVector
192#ifdef debug
193    G4cout<<"G4QGluonString::GetMeanFreePath: isovectorLength="<<isoSize<<G4endl; // Result
194#endif
195    if(isoSize)                             // The Element has non-trivial abumdance set
196    {
197      indEl=pElement->GetIndex();           // Index of the non-trivial element
198      if(!Isotopes->IsDefined(Z,indEl))     // This index is not defined for this Z: define
199      {
200        std::vector<std::pair<G4int,G4double>*>* newAbund =
201                                               new std::vector<std::pair<G4int,G4double>*>;
202        G4double* abuVector=pElement->GetRelativeAbundanceVector();
203        for(G4int j=0; j<isoSize; j++)      // Calculation of abundance vector for isotopes
204        {
205          G4int N=pElement->GetIsotope(j)->GetN()-Z; // N means A=N+Z !
206          if(pElement->GetIsotope(j)->GetZ()!=Z)G4cerr<<"G4QCaptureAtRest::GetMeanFreePath"
207                                 <<": Z="<<pElement->GetIsotope(j)->GetZ()<<"#"<<Z<<G4endl;
208          G4double abund=abuVector[j];
209          std::pair<G4int,G4double>* pr= new std::pair<G4int,G4double>(N,abund);
210#ifdef debug
211          G4cout<<"G4QGluonString::GetMeanFreePath:p#"<<j<<",N="<<N<<",ab="<<abund<<G4endl;
212#endif
213          newAbund->push_back(pr);
214        }
215#ifdef debug
216        G4cout<<"G4QGluonString::PostStepDoIt:pairVectorLength="<<newAbund->size()<<G4endl;
217#endif
218        indEl=G4QIsotope::Get()->InitElement(Z,indEl,newAbund); // definition of the newInd
219        for(G4int k=0; k<isoSize; k++) delete (*newAbund)[k];   // Cleaning temporary
220        delete newAbund; // Was "new" in the beginning of the name space
221      }
222    }
223    std::vector<std::pair<G4int,G4double>*>* cs= Isotopes->GetCSVector(Z,indEl);//CSPointer
224    std::vector<G4double>* SPI = new std::vector<G4double>; // Pointer to the SPI vector
225    IsoProbInEl.push_back(SPI);
226    std::vector<G4int>* IsN = new std::vector<G4int>; // Pointer to the N vector
227    ElIsoN.push_back(IsN);
228    G4int nIs=cs->size();                   // A#Of Isotopes in the Element
229    G4double susi=0.;                       // sum of CS over isotopes
230    if(nIs) for(G4int j=0; j<nIs; j++)      // Calculate CS for eachIsotope of El
231    {
232      std::pair<G4int,G4double>* curIs=(*cs)[j]; // A pointer, which is used twice
233      G4int N=curIs->first;                 // #of Neuterons in the isotope j of El i
234      IsN->push_back(N);                    // Remember Min N for the Element
235      G4double CSI=CSmanager->GetCrossSection(true,Momentum,Z,N,pPDG);//CS(j,i) for isotope
236#ifdef debug
237      G4cout<<"GQC::GMF:X="<<CSI<<",M="<<Momentum<<",Z="<<Z<<",N="<<N<<",P="<<pPDG<<G4endl;
238#endif
239      curIs->second = CSI;                  // Remenber the calculated cross-section
240      susi+=CSI;                            // Make a sum per isotopes
241      SPI->push_back(susi);                 // Remember summed cross-section
242    } // End of temporary initialization of the cross sections in the G4QIsotope singeltone
243    sigma+=Isotopes->GetMeanCrossSection(Z,indEl)*NOfNucPerVolume[i];//SUM(MeanCS*NOfNperV)
244    ElProbInMat.push_back(sigma);
245  } // End of LOOP over Elements
246#ifdef debug
247  G4cout<<"G4QCol::GetMeanFrPa: S="<<sigma<<",e="<<photNucBias<<",w="<<weakNucBias<<G4endl;
248#endif
249  // Check that cross section is not zero and return the mean free path
250  if(sigma > 0.) return 1./sigma;                 // Mean path [distance]
251  return DBL_MAX;                                 // If Sigma=0, return max value for PATH
252}
253
254// Check applicability of the process
255G4bool G4QGluonString::IsApplicable(const G4ParticleDefinition& particle) 
256{
257  if      (particle == *(        G4Proton::Proton()        )) return true;
258  //else if (particle == *(       G4Neutron::Neutron()       )) return true;
259  //else if (particle == *(     G4PionMinus::PionMinus()     )) return true;
260  //else if (particle == *(      G4PionPlus::PionPlus()      )) return true;
261  //else if (particle == *(      G4KaonPlus::KaonPlus()      )) return true;
262  //else if (particle == *(     G4KaonMinus::KaonMinus()     )) return true;
263  //else if (particle == *(  G4KaonZeroLong::KaonZeroLong()  )) return true;
264  //else if (particle == *( G4KaonZeroShort::KaonZeroShort() )) return true;
265  //else if (particle == *(        G4Lambda::Lambda()        )) return true;
266  //else if (particle == *(     G4SigmaPlus::SigmaPlus()     )) return true;
267  //else if (particle == *(    G4SigmaMinus::SigmaMinus()    )) return true;
268  //else if (particle == *(     G4SigmaZero::SigmaZero()     )) return true;
269  //else if (particle == *(       G4XiMinus::XiMinus()       )) return true;
270  //else if (particle == *(        G4XiZero::XiZero()        )) return true;
271  //else if (particle == *(    G4OmegaMinus::OmegaMinus()    )) return true;
272  //else if (particle == *(   G4AntiNeutron::AntiNeutron()   )) return true;
273  //else if (particle == *(    G4AntiProton::AntiProton()    )) return true;
274  //else if (particle == *(    G4AntiLambda::AntiLambda()    )) return true;
275  //else if (particle == *( G4AntiSigmaPlus::AntiSigmaPlus() )) return true;
276  //else if (particle == *(G4AntiSigmaMinus::AntiSigmaMinus())) return true;
277  //else if (particle == *( G4AntiSigmaZero::AntiSigmaZero() )) return true;
278  //else if (particle == *(   G4AntiXiMinus::AntiXiMinus()   )) return true;
279  //else if (particle == *(    G4AntiXiZero::AntiXiZero()    )) return true;
280  //else if (particle == *(G4AntiOmegaMinus::AntiOmegaMinus())) return true;
281  else G4cerr<<"***G4QGluonString::IsApplicable: PDG?="<<particle.GetPDGEncoding()<<G4endl;
282#ifdef debug
283  G4cout<<"***G4QGluonString::IsApplicable: PDG="<<particle.GetPDGEncoding()<<G4endl;
284#endif
285  return false;
286}
287
288G4VParticleChange* G4QGluonString::PostStepDoIt(const G4Track& track, const G4Step& step)
289{
290  //static const G4double dpi=M_PI+M_PI;   // 2*pi (for Phi distr.) ***changed to twopi***
291  //static const G4double mNeut= G4QPDGCode(2112).GetMass();
292  //static const G4double mNeut2= mNeut*mNeut;          // Squared neutron mass
293  //static const G4double mProt= G4QPDGCode(2212).GetMass();
294  //static const G4double mProt2= mProt*mProt;          // Squared Proton mass
295  //static const G4double dM=mProt+mNeut;               // doubled nucleon mass
296  //static const G4double hdM=dM/2.;                    // M of the "nucleon"
297  //static const G4double hdM2=hdM*hdM;                 // M2 of the "nucleon"
298  //static const G4double mPi0 = G4QPDGCode(111).GetMass();
299  //static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0);
300  //static const G4double mPi  = G4QPDGCode(211).GetMass();
301  //static const G4double tmPi = mPi+mPi;               // DoubledMass of the charged pion
302  //static const G4double stmPi= tmPi*tmPi;             // SquareDoubledMass of ChargedPion
303  //static const G4double mPPi = mPi+mProt;             // Delta threshold
304  //static const G4double mPPi2= mPPi*mPPi;             // Delta low threshold for W2
305  //-------------------------------------------------------------------------------------
306  const G4DynamicParticle* projHadron = track.GetDynamicParticle();
307  const G4ParticleDefinition* particle=projHadron->GetDefinition();
308#ifdef debug
309  G4cout<<"G4QGluonString::PostStepDoIt: Before the GetMeanFreePath is called"<<G4endl;
310#endif
311  G4ForceCondition cond=NotForced;
312  GetMeanFreePath(track, 1., &cond);                  // Just to check that still sig>0
313#ifdef debug
314  G4cout<<"G4QGluonString::PostStepDoIt: After the GetMeanFreePath is called"<<G4endl;
315#endif
316  G4LorentzVector proj4M=projHadron->Get4Momentum();
317  G4double momentum = projHadron->GetTotalMomentum(); // 3-momentum of the Particle
318  G4double Momentum=proj4M.rho();
319  if(std::fabs(Momentum-momentum)>.001)
320    G4cerr<<"G4QGluonString::PostStepDoIt: P="<<Momentum<<"="<<momentum<<G4endl;
321#ifdef debug
322  G4double mp=proj4M.m();
323  G4cout<<"G4QGluonString::PostStepDoIt is called, P="<<Momentum<<"="<<momentum<<G4endl;
324#endif
325  if (!IsApplicable(*particle))  // Check applicability
326  {
327    G4cerr<<"G4QGluonString::PostStepDoIt:Only gam,e+,e-,mu+,mu-,t+,t-,p are implemented."
328          <<G4endl;
329    return 0;
330  }
331  const G4Material* material = track.GetMaterial();      // Get the current material
332  G4int Z=0;
333  const G4ElementVector* theElementVector = material->GetElementVector();
334  G4int nE=material->GetNumberOfElements();
335#ifdef debug
336  G4cout<<"G4QGluonString::PostStepDoIt: "<<nE<<" elements in the material."<<G4endl;
337#endif
338  G4int projPDG=0;                           // PDG Code prototype for the captured hadron
339  // Not all these particles are implemented yet (see Is Applicable)
340  if      (particle ==          G4Proton::Proton()         ) projPDG= 2212;
341  //else if (particle ==         G4Neutron::Neutron()        ) projPDG= 2112;
342  //else if (particle ==       G4PionMinus::PionMinus()      ) projPDG= -211;
343  //else if (particle ==        G4PionPlus::PionPlus()       ) projPDG=  211;
344  //else if (particle ==        G4KaonPlus::KaonPlus()       ) projPDG= 2112;
345  //else if (particle ==       G4KaonMinus::KaonMinus()      ) projPDG= -321;
346  //else if (particle ==    G4KaonZeroLong::KaonZeroLong()   ) projPDG=  130;
347  //else if (particle ==   G4KaonZeroShort::KaonZeroShort()  ) projPDG=  310;
348  //else if (particle ==         G4TauPlus::TauPlus()        ) projPDG=  -15;
349  //else if (particle ==        G4TauMinus::TauMinus()       ) projPDG=   15;
350  //else if (particle ==     G4NeutrinoTau::NeutrinoTau()    ) projPDG=   16;
351  //else if (particle == G4AntiNeutrinoTau::AntiNeutrinoTau()) projPDG=  -16;
352  //else if (particle ==          G4Lambda::Lambda()         ) projPDG= 3122;
353  //else if (particle ==       G4SigmaPlus::SigmaPlus()      ) projPDG= 3222;
354  //else if (particle ==      G4SigmaMinus::SigmaMinus()     ) projPDG= 3112;
355  //else if (particle ==       G4SigmaZero::SigmaZero()      ) projPDG= 3212;
356  //else if (particle ==         G4XiMinus::XiMinus()        ) projPDG= 3312;
357  //else if (particle ==          G4XiZero::XiZero()         ) projPDG= 3322;
358  //else if (particle ==      G4OmegaMinus::OmegaMinus()     ) projPDG= 3334;
359  //else if (particle ==     G4AntiNeutron::AntiNeutron()    ) projPDG=-2112;
360  //else if (particle ==      G4AntiProton::AntiProton()     ) projPDG=-2212;
361  //else if (particle ==      G4AntiLambda::AntiLambda()     ) projPDG=-3122;
362  //else if (particle ==   G4AntiSigmaPlus::AntiSigmaPlus()  ) projPDG=-3222;
363  //else if (particle ==  G4AntiSigmaMinus::AntiSigmaMinus() ) projPDG=-3112;
364  //else if (particle ==   G4AntiSigmaZero::AntiSigmaZero()  ) projPDG=-3212;
365  //else if (particle ==     G4AntiXiMinus::AntiXiMinus()    ) projPDG=-3312;
366  //else if (particle ==      G4AntiXiZero::AntiXiZero()     ) projPDG=-3322;
367  //else if (particle ==  G4AntiOmegaMinus::AntiOmegaMinus() ) projPDG=-3334;
368  //G4int aProjPDG=std::abs(projPDG);
369#ifdef debug
370  G4int prPDG=particle->GetPDGEncoding();
371  G4cout<<"G4QGluonString::PostStepDoIt: projPDG="<<projPDG<<", stPDG="<<prPDG<<G4endl;
372#endif
373  if(!projPDG)
374  {
375    G4cerr<<"--Warning--G4QGluonString::PostStepDoIt:Undefined interacting hadron"<<G4endl;
376    return 0;
377  }
378  G4int EPIM=ElProbInMat.size();
379#ifdef debug
380  G4cout<<"G4QCollis::PostStDoIt: m="<<EPIM<<",n="<<nE<<",T="<<ElProbInMat[EPIM-1]<<G4endl;
381#endif
382  G4int i=0;
383  if(EPIM>1)
384  {
385    G4double rnd = ElProbInMat[EPIM-1]*G4UniformRand();
386    for(i=0; i<nE; ++i)
387    {
388#ifdef debug
389      G4cout<<"G4QGluonString::PostStDoIt:E["<<i<<"]="<<ElProbInMat[i]<<",r="<<rnd<<G4endl;
390#endif
391      if (rnd<ElProbInMat[i]) break;
392    }
393    if(i>=nE) i=nE-1;                        // Top limit for the Element
394  }
395  G4Element* pElement=(*theElementVector)[i];
396  Z=static_cast<G4int>(pElement->GetZ());
397#ifdef debug
398    G4cout<<"G4QGluonString::PostStepDoIt: i="<<i<<", Z(element)="<<Z<<G4endl;
399#endif
400  if(Z<=0)
401  {
402    G4cerr<<"---Warning---G4QGluonString::PostStepDoIt: Element with Z="<<Z<<G4endl;
403    if(Z<0) return 0;
404  }
405  std::vector<G4double>* SPI = IsoProbInEl[i];// Vector of summedProbabilities for isotopes
406  std::vector<G4int>* IsN = ElIsoN[i];     // Vector of "#of neutrons" in the isotope El[i]
407  G4int nofIsot=SPI->size();               // #of isotopes in the element i
408#ifdef debug
409  G4cout<<"G4QCollis::PosStDoIt:n="<<nofIsot<<",T="<<(*SPI)[nofIsot-1]<<G4endl;
410#endif
411  G4int j=0;
412  if(nofIsot>1)
413  {
414    G4double rndI=(*SPI)[nofIsot-1]*G4UniformRand(); // Randomize the isotop of the Element
415    for(j=0; j<nofIsot; ++j)
416    {
417#ifdef debug
418      G4cout<<"G4QGluonString::PostStepDoIt:SP["<<j<<"]="<<(*SPI)[j]<<",r="<<rndI<<G4endl;
419#endif
420      if(rndI < (*SPI)[j]) break;
421    }
422    if(j>=nofIsot) j=nofIsot-1;            // Top limit for the isotope
423  }
424  G4int N =(*IsN)[j]; ;                    // Randomized number of neutrons
425#ifdef debug
426  G4cout<<"G4QGluonString::PostStepDoIt: j="<<i<<", N(isotope)="<<N<<G4endl;
427#endif
428  if(N<0)
429  {
430    G4cerr<<"-Warning-G4QGluonString::PostStepDoIt:Isotope with N="<<Z<<"<0,Z="<<Z<<G4endl;
431    return 0;
432  }
433  nOfNeutrons=N;                           // Remember it for the energy-momentum check
434  G4double dd=0.025;
435  G4double am=Z+N;
436  G4double sr=std::sqrt(am);
437  G4double dsr=0.01*(sr+sr);
438  if(dsr<dd)dsr=dd;
439  if(manualFlag) G4QNucleus::SetParameters(freeNuc,freeDib,clustProb,mediRatio);// ManualPa
440  //else if(projPDG==-2212) G4QNucleus::SetParameters(1.-dsr-dsr,dd+dd,5.,10.);//aP CluPars
441  //else if(projPDG==-211)  G4QNucleus::SetParameters(.67-dsr,.32-dsr,5.,9.); //Pi- CluPars
442#ifdef debug
443  G4cout<<"G4QGluonString::PostStepDoIt: N="<<N<<" for element with Z="<<Z<<G4endl;
444#endif
445  if(N<0)
446  {
447    G4cerr<<"---Warning---G4QGluonString::PostStepDoIt:Element with N="<<N<< G4endl;
448    return 0;
449  }
450  aParticleChange.Initialize(track);
451  G4double localtime = track.GetGlobalTime();
452  G4ThreeVector position = track.GetPosition();
453  G4TouchableHandle trTouchable = track.GetTouchableHandle();
454  //
455  G4int targPDG=90000000+Z*1000+N;            // PDG Code of the target nucleus
456  //    =========================
457  G4QPDGCode targQPDG(targPDG);
458  G4double tM=targQPDG.GetMass();
459  G4QHadronVector* output=new G4QHadronVector;// Prototype of EnvironOutput G4QHadronVector
460  G4double absMom = 0.;                       // Prototype of absorbed by nucleus Moment
461  G4QHadronVector* leadhs=new G4QHadronVector;// Prototype of QuasmOutput G4QHadronVectorum
462  G4LorentzVector lead4M(0.,0.,0.,0.);        // Prototype of LeadingQ 4-momentum
463  EnMomConservation=proj4M+G4LorentzVector(0.,0.,0.,tM);    // Total 4-mom of the reaction
464  if(absMom) EnMomConservation+=lead4M;       // Add E/M of leading System
465#ifdef debug
466  G4cout<<"G4QGluonString::PostStepDoIt: projPDG="<<projPDG<<", targPDG="<<targPDG<<G4endl;
467#endif
468  //G4QHadron* pH = new G4QHadron(projPDG,proj4M);
469  // @@@@@@@@@@@@@@@@@@@@@@@@@@ Hrere the CHIPS_QGS must be implemented @@@@@@@@@@@@@@@
470  G4int qNH=leadhs->size();
471  if(absMom)
472  {
473    if(qNH) for(G4int iq=0; iq<qNH; iq++)
474    {
475      G4QHadron* loh=(*leadhs)[iq];   // Pointer to the output hadron
476      output->push_back(loh);
477    }
478    delete leadhs;
479  }
480  // ------------- From here the secondaries are filled -------------------------
481  G4int tNH = output->size();       // A#of hadrons in the output
482  aParticleChange.SetNumberOfSecondaries(tNH); 
483  // Now add nuclear fragments
484#ifdef debug
485  G4cout<<"G4QGluonString::PostStepDoIt: "<<tNH<<" particles are generated"<<G4endl;
486#endif
487#ifdef ppdebug
488  if(absMom)G4cout<<"G4QGluonString::PostStepDoIt: t="<<tNH<<", q="<<qNH<<G4endl;
489#endif
490  G4int nOut=output->size();               // Real length of the output @@ Temporary
491  if(tNH==1)                               // @@ Temporary. Find out why it happened!
492  {
493    G4cout<<"-Warning-G4QGluonString::PostStepDoIt: 1 secondary! absMom="<<absMom;
494    if(absMom) G4cout<<", qNH="<<qNH;
495    G4cout<<", PDG0="<<(*output)[0]->GetPDGCode();
496    G4cout<<G4endl;
497    tNH=0;
498    delete output->operator[](0);          // delete the creazy hadron
499    output->pop_back();                    // clean up the output vector
500  }
501  if(tNH==2&&2!=nOut) G4cout<<"--Warning--G4QGluonString::PostStepDoIt:2 # "<<nOut<<G4endl;
502  // Deal with ParticleChange final state interface to GEANT4 output of the process
503  //if(tNH==2) for(i=0; i<tNH; i++) // @@ Temporary tNH==2 instead of just tNH
504  if(tNH) for(i=0; i<tNH; i++) // @@ Temporary tNH==2 instead of just tNH
505  {
506    // Note that one still has to take care of Hypernuclei (with Lambda or Sigma inside)
507    // Hypernucleus mass calculation and ion-table interface upgrade => work for Hisaya @@
508    // The decau process for hypernuclei must be developed in GEANT4 (change CHIPS body)
509    G4QHadron* hadr=(*output)[i];          // Pointer to the output hadron   
510    G4int PDGCode = hadr->GetPDGCode();
511    G4int nFrag   = hadr->GetNFragments();
512#ifdef pdebug
513    G4cout<<"G4QGluonString::AtRestDoIt: H#"<<i<<",PDG="<<PDGCode<<",nF="<<nFrag<<G4endl;
514#endif
515    if(nFrag)                // Skip intermediate (decayed) hadrons
516    {
517#ifdef debug
518      G4cout<<"G4QGluonString::PostStepDoIt: Intermediate particle is found i="<<i<<G4endl;
519#endif
520      delete hadr;
521      continue;
522    }
523    G4DynamicParticle* theSec = new G4DynamicParticle; 
524    G4ParticleDefinition* theDefinition;
525    if     (PDGCode==90000001) theDefinition = G4Neutron::Neutron();
526    else if(PDGCode==90001000) theDefinition = G4Proton::Proton();//While it can be in ions
527    else if(PDGCode==91000000) theDefinition = G4Lambda::Lambda();
528    else if(PDGCode==311 || PDGCode==-311)
529    {
530      if(G4UniformRand()>.5) theDefinition = G4KaonZeroLong::KaonZeroLong();   // K_L
531      else                   theDefinition = G4KaonZeroShort::KaonZeroShort(); // K_S
532    }
533    else if(PDGCode==91000999) theDefinition = G4SigmaPlus::SigmaPlus();
534    else if(PDGCode==90999001) theDefinition = G4SigmaMinus::SigmaMinus();
535    else if(PDGCode==91999000) theDefinition = G4XiMinus::XiMinus();
536    else if(PDGCode==91999999) theDefinition = G4XiZero::XiZero();
537    else if(PDGCode==92998999) theDefinition = G4OmegaMinus::OmegaMinus();
538    else if(PDGCode >80000000) // Defines hypernuclei as normal nuclei (N=N+S Correction!)
539    {
540      G4int aZ = hadr->GetCharge();
541      G4int aA = hadr->GetBaryonNumber();
542#ifdef pdebug
543      G4cout<<"G4QGluonString::AtRestDoIt:Ion Z="<<aZ<<", A="<<aA<<G4endl;
544#endif
545      theDefinition = G4ParticleTable::GetParticleTable()->FindIon(aZ,aA,0,aZ);
546    }
547    //else theDefinition = G4ParticleTable::GetParticleTable()->FindParticle(PDGCode);
548    else
549    {
550#ifdef pdebug
551      G4cout<<"G4QGluonString::PostStepDoIt:Define particle with PDG="<<PDGCode<<G4endl;
552#endif
553      theDefinition = G4QPDGToG4Particle::Get()->GetParticleDefinition(PDGCode);
554#ifdef pdebug
555      G4cout<<"G4QGluonString::PostStepDoIt:AfterParticleDefinition PDG="<<PDGCode<<G4endl;
556#endif
557    }
558    if(!theDefinition)
559    {
560#ifdef debug
561      G4cout<<"---Warning---G4QGluonString::PostStepDoIt: drop PDG="<<PDGCode<<G4endl;
562#endif
563      delete hadr;
564      continue;
565    }
566#ifdef pdebug
567    G4cout<<"G4QGluonString::PostStepDoIt:Name="<<theDefinition->GetParticleName()<<G4endl;
568#endif
569    theSec->SetDefinition(theDefinition);
570    G4LorentzVector h4M=hadr->Get4Momentum();
571    EnMomConservation-=h4M;
572#ifdef tdebug
573    G4cout<<"G4QCollis::PSDI:"<<i<<","<<PDGCode<<h4M<<h4M.m()<<EnMomConservation<<G4endl;
574#endif
575#ifdef debug
576    G4cout<<"G4QGluonString::PostStepDoIt:#"<<i<<",PDG="<<PDGCode<<",4M="<<h4M<<G4endl;
577#endif
578    theSec->Set4Momentum(h4M); //                                                         ^
579    delete hadr; // <-----<-----------<-------------<---------------------<---------<-----+
580#ifdef debug
581    G4ThreeVector curD=theSec->GetMomentumDirection();              //                    ^
582    G4double curM=theSec->GetMass();                                //                    |
583    G4double curE=theSec->GetKineticEnergy()+curM;                  //                    ^
584    G4cout<<"G4QCollis::PSDoIt:p="<<curD<<curD.mag()<<",e="<<curE<<",m="<<curM<<G4endl;// |
585#endif
586    G4Track* aNewTrack = new G4Track(theSec, localtime, position ); //                    ^
587    aNewTrack->SetTouchableHandle(trTouchable);                     //                    |
588    aParticleChange.AddSecondary( aNewTrack );                      //                    |
589#ifdef debug
590    G4cout<<"G4QGluonString::PostStepDoIt:#"<<i<<" is done"<<G4endl;//                    |
591#endif
592  } //                                                                                    |
593  delete output; // instances of the G4QHadrons from the output are already deleted above +
594#ifdef debug
595  G4cout<<"G4QGluonString::PostStDoIt: afterSt="<<aParticleChange.GetTrackStatus()<<G4endl;
596#endif
597  aParticleChange.ProposeTrackStatus(fStopAndKill);        // Kill the absorbed particle
598#ifdef debug
599  G4cout<<"G4QGluonString::PostStepDoIt:*** PostStepDoIt is done ***, P="<<aProjPDG
600        <<", St="<<aParticleChange.GetTrackStatus()<<G4endl;
601#endif
602  return G4VDiscreteProcess::PostStepDoIt(track, step);
603}
Note: See TracBrowser for help on using the repository browser.