250{
255
256
257 static G4bool CWinit =
true;
258 if(CWinit)
259 {
260 CWinit=false;
263 }
264
267#ifdef debug
268 G4cout<<
"G4QDiffraction::PostStepDoIt: Before the GetMeanFreePath is called In4M="
271#endif
274#ifdef debug
275 G4cout<<
"G4QDiffraction::PostStepDoIt: After GetMeanFreePath is called"<<
G4endl;
276#endif
280 if(std::fabs(Momentum-momentum)>.000001)
281 G4cerr<<
"-Warning-G4QDiffraction::PostStepDoIt:P_IU="<<Momentum<<
"#"<<momentum<<
G4endl;
282#ifdef pdebug
283 G4cout<<
"G4QDiffraction::PostStepDoIt: pP(IU)="<<Momentum<<
"="<<momentum
284 <<
",proj4M="<<proj4M<<
", projM="<<proj4M.
m()<<
G4endl;
285#endif
287 {
288 G4cerr<<
"G4QDiffraction::PostStepDoIt: Only NA is implemented."<<
G4endl;
289 return 0;
290 }
295#ifdef debug
296 G4cout<<
"G4QDiffraction::PostStepDoIt: "<<nE<<
" elements in the material."<<
G4endl;
297#endif
299
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330#ifdef debug
332 G4cout<<
"G4QDiffraction::PostStepDoIt: projPDG="<<projPDG<<
", stPDG="<<prPDG<<
G4endl;
333#endif
334 if(!projPDG)
335 {
336 G4cerr<<
"-Warning-G4QDiffraction::PostStepDoIt:UndefProjHadron(PDG=0) ->ret 0"<<
G4endl;
337 return 0;
338 }
339
340
342 if(projPDG==2112) pM=mProt;
343
344 G4int EPIM=ElProbInMat.size();
345#ifdef debug
346 G4cout<<
"G4QDiffra::PostStDoIt: m="<<EPIM<<
",n="<<nE<<
",T="<<ElProbInMat[EPIM-1]<<
G4endl;
347#endif
349 if(EPIM>1)
350 {
352 for(i=0; i<nE; ++i)
353 {
354#ifdef debug
355 G4cout<<
"G4QDiffra::PostStepDoIt: EPM["<<i<<
"]="<<ElProbInMat[i]<<
",r="<<rnd<<
G4endl;
356#endif
357 if (rnd<ElProbInMat[i]) break;
358 }
359 if(i>=nE) i=nE-1;
360 }
361 G4Element* pElement=(*theElementVector)[i];
363#ifdef debug
364 G4cout<<
"G4QDiffraction::PostStepDoIt: i="<<i<<
", Z(element)="<<Z<<
G4endl;
365#endif
366 if(Z<=0)
367 {
368 G4cerr<<
"-Warning-G4QDiffraction::PostStepDoIt: Element with Z="<<Z<<
G4endl;
369 if(Z<0) return 0;
370 }
371 std::vector<G4double>* SPI = IsoProbInEl[i];
372 std::vector<G4int>* IsN = ElIsoN[i];
373 G4int nofIsot=SPI->size();
374#ifdef debug
375 G4cout<<
"G4QDiffraction::PostStepDoIt: nI="<<nofIsot<<
", T="<<(*SPI)[nofIsot-1]<<
G4endl;
376#endif
378 if(nofIsot>1)
379 {
381 for(j=0; j<nofIsot; ++j)
382 {
383#ifdef debug
384 G4cout<<
"G4QDiffraction::PostStepDoIt: SP["<<j<<
"]="<<(*SPI)[j]<<
",r="<<rndI<<
G4endl;
385#endif
386 if(rndI < (*SPI)[j]) break;
387 }
388 if(j>=nofIsot) j=nofIsot-1;
389 }
390 G4int N =(*IsN)[j]; ;
391#ifdef debug
392 G4cout<<
"G4QDiffraction::PostStepDoIt: j="<<i<<
", N(isotope)="<<N<<
", MeV="<<MeV<<
G4endl;
393#endif
394 if(N<0)
395 {
396 G4cerr<<
"-Warning-G4QDiffraction::PostStepDoIt: Isotope Z="<<Z<<
" has 0>N="<<N<<
G4endl;
397 return 0;
398 }
399 nOfNeutrons=N;
400#ifdef debug
401 G4cout<<
"G4QDiffraction::PostStepDoIt: N="<<N<<
" for element with Z="<<Z<<
G4endl;
402#endif
403 if(N<0)
404 {
405 G4cerr<<
"*Warning*G4QDiffraction::PostStepDoIt:Element with N="<<N<<
G4endl;
406 return 0;
407 }
409#ifdef debug
410 G4cout<<
"G4QDiffraction::PostStepDoIt: track is initialized"<<
G4endl;
411#endif
412 G4double weight = track.GetWeight();
413 G4double localtime = track.GetGlobalTime();
415#ifdef debug
416 G4cout<<
"G4QDiffraction::PostStepDoIt: before Touchable extraction"<<
G4endl;
417#endif
419#ifdef debug
420 G4cout<<
"G4QDiffraction::PostStepDoIt: Touchable is extracted"<<
G4endl;
421#endif
422 G4int targPDG=90000000+Z*1000+N;
428#ifdef debug
429 G4cout<<
"G4QDiffraction::PostStepDoIt: tM="<<tM<<
",p4M="<<proj4M<<
",t4M="<<tot4M<<
G4endl;
430#endif
431 EnMomConservation=tot4M;
432
433#ifdef debug
434 G4cout<<
"G4QDiff::PSDI:false,P="<<Momentum<<
",Z="<<Z<<
",N="<<N<<
",PDG="<<projPDG<<
G4endl;
435#endif
436 G4double xSec=CalculateXS(Momentum, Z, N, projPDG);
437#ifdef debug
438 G4cout<<
"G4QDiffra::PSDI:PDG="<<projPDG<<
",P="<<Momentum<<
",XS="<<xSec/millibarn<<
G4endl;
439#endif
440#ifdef nandebug
441 if(xSec>0. || xSec<0. || xSec==0);
442 else G4cout<<
"-Warning-G4QDiffraction::PostSDI: *NAN* xSec="<<xSec/millibarn<<
G4endl;
443#endif
444
445 if(xSec <= 0.)
446 {
447#ifdef pdebug
448 G4cerr<<
"*Warning*G4QDiffraction::PSDoIt:*Zero cross-section* PDG="<<projPDG
449 <<
",tPDG="<<targPDG<<
",P="<<Momentum<<
G4endl;
450#endif
451
456 }
458 if(totCMMass < mPion+pM+tM)
459 {
460#ifdef pdebug
461 G4cerr<<
"*Warning*G4QDiffraction::PSDoIt:*Below Diffraction Threshold* cmM="<<totCMMass
462 <<
">pM="<<pM<<
"+tM="<<tM<<
"+pi0="<<mPion<<
"=="<<pM+tM+mPion<<
G4endl;
463#endif
464
469 }
470
473 G4int nSec=out->size();
478#ifdef pdebug
479 G4cout<<
"G4QDiffraction::PostStepDoIt: =---=found=---= nSecondaries="<<nSec<<
G4endl;
480#endif
481 for(i=0; i<nSec; i++)
482 {
483 difQH = (*out)[i];
494 else if(difPDG== 130 || difPDG==-311 || difPDG==89000001)
496 else if(difPDG== 310 || difPDG== 311 || difPDG==90999999)
518 else
519 {
523 if(S||Z>B||Z<0)
G4cout<<
"-Warning-G4QDif::PoStDoIt:Z="<<Z<<
",A="<<B<<
",S="<<S<<
G4endl;
525#ifdef pdebug
526 G4cout<<
"G4QDiffraction::PoStDoIt:Ion,Z="<<Z<<
",A="<<B<<
",D="<<theDefinition<<
G4endl;
527#endif
528 }
529 if(theDefinition)
530 {
534 EnMomConservation-=dif4M;
540#ifdef pdebug
541 G4cout<<
"G4QDiffraction::PostStepDoIt: Filled 4M="<<dif4M<<
", PDG="<<difPDG<<
G4endl;
542#endif
543 }
544 else G4cout<<
"-Warning-G4QDif::PSDI: Lost PDG="<<difPDG<<
", Z="<<difQH->
GetCharge()
546 delete difQH;
547 }
548 delete out;
549#ifdef debug
550 G4cout<<
"G4QDiffraction::PostStepDoIt:*** PostStepDoIt is done ***"<<
G4endl;
551#endif
553}
CLHEP::HepLorentzVector G4LorentzVector
std::vector< G4QHadron * > G4QHadronVector
static G4AntiLambda * AntiLambda()
static G4AntiNeutron * AntiNeutron()
static G4AntiOmegaMinus * AntiOmegaMinus()
static G4AntiProton * AntiProton()
static G4AntiSigmaMinus * AntiSigmaMinus()
static G4AntiSigmaPlus * AntiSigmaPlus()
static G4AntiSigmaZero * AntiSigmaZero()
static G4AntiXiMinus * AntiXiMinus()
static G4AntiXiZero * AntiXiZero()
const G4ThreeVector & GetMomentumDirection() const
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
G4LorentzVector Get4Momentum() const
void Set4Momentum(const G4LorentzVector &momentum)
static G4Electron * Electron()
static G4KaonMinus * KaonMinus()
static G4KaonPlus * KaonPlus()
static G4KaonZeroLong * KaonZeroLong()
static G4KaonZeroShort * KaonZeroShort()
static G4Lambda * Lambda()
static G4MuonMinus * MuonMinus()
static G4MuonPlus * MuonPlus()
static G4OmegaMinus * OmegaMinus()
void AddSecondary(G4Track *aSecondary)
void ProposeEnergy(G4double finalEnergy)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
virtual void Initialize(const G4Track &)
const G4String & GetParticleType() const
const G4String & GetParticleSubType() const
G4ParticleDefinition * FindIon(G4int atomicNumber, G4int atomicMass, G4double excitationEnergy)
static G4ParticleTable * GetParticleTable()
static G4PionMinus * PionMinus()
static G4PionPlus * PionPlus()
static G4PionZero * PionZero()
static G4Positron * Positron()
static G4QDiffractionRatio * GetPointer()
G4QHadronVector * TargFragment(G4int pPDG, G4LorentzVector p4M, G4int tgZ, G4int tgN)
G4double GetMeanFreePath(const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition)
G4LorentzVector Get4Momentum() const
G4int GetBaryonNumber() const
G4int GetStrangeness() const
static G4SigmaMinus * SigmaMinus()
static G4SigmaPlus * SigmaPlus()
static G4SigmaZero * SigmaZero()
void SetWeight(G4double aValue)
void SetTouchableHandle(const G4TouchableHandle &apValue)
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
void ProposeTrackStatus(G4TrackStatus status)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4ParticleChange aParticleChange
static G4XiMinus * XiMinus()
static G4XiZero * XiZero()