9 #include "G4Electron.hh"
10 #include "G4Positron.hh"
11 #include "G4PionPlus.hh"
12 #include "G4PionMinus.hh"
13 #include "G4PionZero.hh"
14 #include "G4KaonPlus.hh"
15 #include "G4KaonMinus.hh"
16 #include "G4KaonZeroShort.hh"
17 #include "G4KaonZeroLong.hh"
18 #include "G4KaonZero.hh"
19 #include "G4AntiKaonZero.hh"
20 #include "G4MuonMinus.hh"
21 #include "G4MuonPlus.hh"
22 #include "G4Neutron.hh"
23 #include "G4AntiNeutron.hh"
24 #include "G4Proton.hh"
25 #include "G4AntiProton.hh"
26 #include "G4Lambda.hh"
27 #include "G4AntiLambda.hh"
28 #include "G4SigmaPlus.hh"
29 #include "G4AntiSigmaPlus.hh"
30 #include "G4SigmaMinus.hh"
31 #include "G4AntiSigmaMinus.hh"
32 #include "G4SigmaZero.hh"
33 #include "G4AntiSigmaZero.hh"
34 #include "G4XiZero.hh"
35 #include "G4AntiXiZero.hh"
36 #include "G4XiMinus.hh"
37 #include "G4AntiXiMinus.hh"
38 #include "G4OmegaMinus.hh"
39 #include "G4AntiOmegaMinus.hh"
42 #include "G4Deuteron.hh"
43 #include "G4Triton.hh"
45 #include "G4Geantino.hh"
47 #include "G4ParticleTable.hh"
48 #include "G4IonTable.hh"
49 #include "G4VProcess.hh"
51 #include "G4StepPoint.hh"
52 #include "G4VSensitiveDetector.hh"
53 #include "Randomize.hh"
409 m_tgpSiA =
new TGraph(fpSiA.c_str());
411 m_tgpSiB =
new TGraph(fpSiB.c_str());
413 m_tgnSiA =
new TGraph(fnSiA.c_str());
415 m_tgnSiB =
new TGraph(fnSiB.c_str());
417 m_tgnSiC =
new TGraph(fnSiC.c_str());
419 m_tgpiSi =
new TGraph(fpiSi.c_str());
421 m_tgeSi =
new TGraph(feSi.c_str());
439 m_tgHn =
new TGraph(fHn.c_str() ,
"%lg %*lg %*lg %*lg %*lg %*lg %lg");
441 m_tgHg =
new TGraph(fHg.c_str() ,
"%lg %*lg %*lg %*lg %*lg %*lg %lg");
443 m_tgHp =
new TGraph(fHp.c_str() ,
"%lg %*lg %*lg %*lg %*lg %*lg %lg");
445 m_tgHem =
new TGraph(fHem.c_str(),
"%lg %*lg %*lg %lg");
447 m_tgHep =
new TGraph(fHep.c_str(),
"%lg %*lg %*lg %lg");
449 m_tgHmm =
new TGraph(fHmm.c_str(),
"%lg %*lg %*lg %lg");
451 m_tgHmp =
new TGraph(fHmp.c_str(),
"%lg %*lg %*lg %lg");
453 m_tgHpm =
new TGraph(fHpm.c_str(),
"%lg %*lg %*lg %lg");
455 m_tgHpp =
new TGraph(fHpp.c_str(),
"%lg %*lg %*lg %lg");
457 m_tgHa =
new TGraph(fHa.c_str() ,
"%lg %*lg %*lg %lg");
462 bool goodMaterial(
false);
473 if (
m_config.
materials[
i].compare(aStep->GetTrack()->GetMaterial()->GetName()) == 0) {
482 const std::vector<const G4Track*>* tv = aStep->GetSecondaryInCurrentStep();
487 for (
unsigned int is=0;is<tv->size();is++) {
489 const G4ParticleDefinition *
pd = (*tv)[is]->GetParticleDefinition ();
490 if (
pd == G4Triton::Definition() || (
pd->GetParticleTable()->GetIonTable()->IsIon(
pd) && !
pd->GetPDGStable() ) ) {
491 G4String svname(
"unknown");
492 G4String svmaterial(
"unknown");
493 if ( (*tv)[is]->GetVolume() ) {
494 svname = (*tv)[is]->GetVolume()->GetName();
495 if ( (*tv)[is]->GetVolume()->GetLogicalVolume() ) {
496 if ( (*tv)[is]->GetVolume()->GetLogicalVolume()->GetMaterial() ) {
497 svmaterial = (*tv)[is]->GetVolume()->GetLogicalVolume()->GetMaterial()->GetName();
503 <<
" Proc.: " << (*tv)[is]->GetCreatorProcess()->GetProcessType() <<
" Mother: ";
505 m_activationFile << aStep->GetTrack()->GetDefinition()->GetParticleName() <<
" (x,y,z,t): (";
527 if (aStep->GetTrack()->GetDefinition()==G4Gamma::Definition()){
529 }
else if (aStep->GetTrack()->GetDefinition()==G4Proton::Definition()){
531 }
else if (aStep->GetTrack()->GetDefinition()==G4PionPlus::Definition() ||
532 aStep->GetTrack()->GetDefinition()==G4PionMinus::Definition()){
534 }
else if(aStep->GetTrack()->GetDefinition()==G4MuonPlus::Definition() ||
535 aStep->GetTrack()->GetDefinition()==G4MuonMinus::Definition()){
537 }
else if(aStep->GetTrack()->GetDefinition()==G4Electron::Definition() ||
538 aStep->GetTrack()->GetDefinition()==G4Positron::Definition()){
539 if (aStep->GetTrack()->GetKineticEnergy()>0.5){
544 }
else if(aStep->GetTrack()->GetDefinition()==G4Neutron::Definition()){
545 if (aStep->GetTrack()->GetKineticEnergy()>10.){
550 }
else if (aStep->GetTrack()->GetDefinition()==G4KaonPlus::Definition() ||
551 aStep->GetTrack()->GetDefinition()==G4KaonMinus::Definition() ||
552 aStep->GetTrack()->GetDefinition()==G4KaonZeroShort::Definition() ||
553 aStep->GetTrack()->GetDefinition()==G4KaonZeroLong::Definition() ||
554 aStep->GetTrack()->GetDefinition()==G4KaonZero::Definition() ||
555 aStep->GetTrack()->GetDefinition()==G4AntiKaonZero::Definition() ||
556 aStep->GetTrack()->GetDefinition()==G4PionZero::Definition()){
559 }
else if (aStep->GetTrack()->GetDefinition()==G4AntiProton::Definition() ||
560 aStep->GetTrack()->GetDefinition()==G4AntiNeutron::Definition() ||
561 aStep->GetTrack()->GetDefinition()==G4Lambda::Definition() ||
562 aStep->GetTrack()->GetDefinition()==G4AntiLambda::Definition() ||
563 aStep->GetTrack()->GetDefinition()==G4SigmaPlus::Definition() ||
564 aStep->GetTrack()->GetDefinition()==G4AntiSigmaPlus::Definition() ||
565 aStep->GetTrack()->GetDefinition()==G4SigmaMinus::Definition() ||
566 aStep->GetTrack()->GetDefinition()==G4AntiSigmaMinus::Definition() ||
567 aStep->GetTrack()->GetDefinition()==G4SigmaZero::Definition() ||
568 aStep->GetTrack()->GetDefinition()==G4AntiSigmaZero::Definition() ||
569 aStep->GetTrack()->GetDefinition()==G4XiMinus::Definition() ||
570 aStep->GetTrack()->GetDefinition()==G4AntiXiMinus::Definition() ||
571 aStep->GetTrack()->GetDefinition()==G4XiZero::Definition() ||
572 aStep->GetTrack()->GetDefinition()==G4AntiXiZero::Definition() ||
573 aStep->GetTrack()->GetDefinition()==G4OmegaMinus::Definition() ||
574 aStep->GetTrack()->GetDefinition()==G4AntiOmegaMinus::Definition() ||
575 aStep->GetTrack()->GetDefinition()==G4Alpha::Definition() ||
576 aStep->GetTrack()->GetDefinition()==G4Deuteron::Definition() ||
577 aStep->GetTrack()->GetDefinition()==G4Triton::Definition() ||
578 aStep->GetTrack()->GetDefinition()==G4He3::Definition()){
581 }
else if (aStep->GetTrack()->GetDefinition()==G4Geantino::Definition()) {
583 }
else if (!aStep->GetTrack()->GetDefinition()->GetPDGCharge())
return;
585 if ( pdgid == 10 && aStep->GetTrack()->GetKineticEnergy() > 10 ) {
587 if ( aStep->GetTrack()->GetDefinition()->GetAtomicNumber() > 1 ) {
594 if ( pdgid == 0 || pdgid == 3 ||
595 pdgid == 1 || pdgid == 2 || pdgid == 4 || pdgid == 5 || pdgid == 6 || pdgid == 7 || pdgid == 8 || pdgid == 9 ||
596 aStep->GetTotalEnergyDeposit() > 0 || pdgid == 999) {
598 double charge = aStep->GetTrack()->GetDefinition()->GetPDGCharge();
600 double absq = std::fabs(
charge);
604 const G4Material * theMaterial = aStep->GetTrack()->GetMaterial();
607 G4StepPoint* pre_step_point = aStep->GetPreStepPoint();
608 G4StepPoint* post_step_point = aStep->GetPostStepPoint();
609 const G4ThreeVector& startPoint = pre_step_point->GetPosition();
610 const G4ThreeVector& endPoint = post_step_point->GetPosition();
611 G4ThreeVector
p = (startPoint + endPoint) * 0.5;
616 double timeOfFlight = (pre_step_point->GetGlobalTime() +
617 post_step_point->GetGlobalTime()) * 0.5;
627 double x0 = aStep->GetPreStepPoint()->GetPosition().x()*0.1;
628 double x1 = aStep->GetPostStepPoint()->GetPosition().x()*0.1;
629 double y0 = aStep->GetPreStepPoint()->GetPosition().y()*0.1;
630 double y1 = aStep->GetPostStepPoint()->GetPosition().y()*0.1;
631 double z0 = aStep->GetPreStepPoint()->GetPosition().z()*0.1;
632 double z1 = aStep->GetPostStepPoint()->GetPosition().z()*0.1;
637 unsigned int nStep =
l/dl0+1;
638 double dx = (
x1-x0)/nStep;
639 double dy = (
y1-y0)/nStep;
640 double dz = (z1-
z0)/nStep;
645 double eKin = aStep->GetTrack()->GetKineticEnergy();
646 double theta = aStep->GetTrack()->GetMomentumDirection().theta()*180./
M_PI;
650 double dphi = (aStep->GetTrack()->GetMomentumDirection().phi()-
p.phi())*180./
M_PI;
651 while ( dphi > 360 ) dphi -= 360.;
652 while ( dphi < 0 ) dphi += 360.;
656 if ( pdgid == 1 || pdgid == 9 ) {
675 else if ( pdgid == 2 || pdgid == 8 ) {
688 else if ( pdgid == 6 || pdgid == 7 ) {
695 else if ( eKin < 800 ) {
702 else if ( pdgid == 4 || pdgid == 5 ) {
715 else if ( pdgid == 0 ) {
720 else if ( pdgid == 3 ) {
731 double dE_TOT = aStep->GetTotalEnergyDeposit()/nStep;
732 double dE_NIEL = aStep->GetNonIonizingEnergyDeposit()/nStep;
733 double dE_ION = dE_TOT-dE_NIEL;
734 double dH_T = 1
e-12*wHt;
736 for(
unsigned int i=0;
i<nStep;
i++) {
738 double subpos = G4UniformRand();
739 double abszorz =
z0+dz*(
i+subpos);
743 double rr = sqrt(
pow(x0+
dx*(
i+subpos),2)+
745 double pphi = atan2(y0+
dy*(
i+subpos),x0+
dx*(
i+subpos))*180/
M_PI;
750 int vBinZoomSpecn = -1;
751 int vBinFullSpecn = -1;
752 int vBinZoomSpeco = -1;
753 int vBinFullSpeco = -1;
754 int vBinZoomTime = -1;
755 int vBinFullTime = -1;
756 int vBinThetaFullSpecn = -1;
757 int vBinThetaFullSpeco = -1;
773 (pdgid == 6 || pdgid == 7)) {
779 pdgid != 6 && pdgid != 7) {
800 (pdgid == 6 || pdgid == 7)) {
813 pdgid != 6 && pdgid != 7) {
836 double phi_mapped = pphi;
838 if (phi_mapped < 0) {
839 phi_mapped = 360 + phi_mapped;
855 if ( goodMaterial && vBinZoom >=0 ) {
856 if ( pdgid == 999 ) {
859 for (
size_t ie=0;
ie<theMaterial->GetNumberOfElements();
ie++ ) {
860 const G4Element * theElement = theMaterial->GetElement (
ie);
861 const double mFrac = theMaterial->GetFractionVector()[
ie];
862 const int zElem = theElement->GetZ();
872 for (
size_t ie=0;
ie<theMaterial->GetNumberOfElements();
ie++ ) {
873 const G4Element * theElement = theMaterial->GetElement (
ie);
874 const double mFrac = theMaterial->GetFractionVector()[
ie];
875 const int zElem = theElement->GetZ();
883 if ( goodMaterial && vBinZoomTime >=0 ) {
889 if ( goodMaterial && vBinFull >=0 ) {
890 if ( pdgid == 999 ) {
893 for (
size_t ie=0;
ie<theMaterial->GetNumberOfElements();
ie++ ) {
894 const G4Element * theElement = theMaterial->GetElement (
ie);
895 const double mFrac = theMaterial->GetFractionVector()[
ie];
896 const int zElem = theElement->GetZ();
906 for (
size_t ie=0;
ie<theMaterial->GetNumberOfElements();
ie++ ) {
907 const G4Element * theElement = theMaterial->GetElement (
ie);
908 const double mFrac = theMaterial->GetFractionVector()[
ie];
909 const int zElem = theElement->GetZ();
917 if ( goodMaterial && vBinFullTime >=0 ) {
923 if ( goodMaterial && vBin3d >=0 ) {
924 if ( pdgid == 999 ) {
934 if ( goodMaterial && (pdgid == 1 || pdgid == 2 || pdgid == 4 || pdgid == 5 || pdgid == 6 || pdgid == 7 || pdgid == 8 || pdgid == 9 )) {
937 if ( vBinZoom >=0 ) {
940 if ( vBinFull >=0 ) {
948 if ( eKin > 20 && (pdgid == 1 || pdgid == 2 || pdgid == 6 || pdgid == 7 || pdgid == 8 || pdgid == 9) ) {
949 if ( vBinZoom >=0 ) {
952 if ( vBinFull >=0 ) {
960 if ( eKin > 0.1 && (pdgid == 6 || pdgid == 7 ) ) {
961 if ( vBinZoom >=0 ) {
964 if ( vBinFull >=0 ) {
972 if ( absq > 0 && (pdgid == 1 || pdgid == 2 || pdgid == 8 || pdgid == 9 ) ) {
973 if ( vBinZoom >=0 ) {
976 if ( vBinFull >=0 ) {
984 if ( pdgid == 6 || pdgid == 7 ) {
985 if ( vBinZoomSpecn >=0 ) {
988 if ( vBinFullSpecn >=0 ) {
991 if ( vBinThetaFullSpecn >=0 ) {
997 if ( goodMaterial && (pdgid < 6 || pdgid > 7 )){
999 if ( vBinZoomSpeco >=0 ) {
1003 else if ( pdgid == 1 ) {
1006 else if ( pdgid == 2 ) {
1009 else if ( pdgid == 3 ) {
1012 else if ( pdgid == 4 || pdgid == 5 ) {
1019 if ( vBinFullSpeco >=0 ) {
1023 else if ( pdgid == 1 ) {
1026 else if ( pdgid == 2 ) {
1029 else if ( pdgid == 3 ) {
1032 else if ( pdgid == 4 || pdgid == 5 ) {
1039 if ( vBinThetaFullSpeco >=0 ) {
1043 else if ( pdgid == 1 ) {
1046 else if ( pdgid == 2 ) {
1049 else if ( pdgid == 3 ) {
1052 else if ( pdgid == 4 || pdgid == 5 ) {
1068 if ( (eKin > 1 && (pdgid == 6 || pdgid == 7)) || pdgid == 999) {
1082 if ( vBinZoom >=0 ) {
1085 if ( vBinFull >=0 ) {
1091 if ( goodMaterial ) {
1095 if ( vBinZoom >=0 ) {
1098 if ( vBinFull >=0 ) {