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 ) {