STCyclotron/src/STCyclotronSensitiveTarget.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //  Author: F. Poignant, floriane.poignant@gmail.com
0027 //
0028 // file STCyclotronSensitiveTarget.cc
0029 //
0030 #include "STCyclotronRun.hh"
0031 #include "STCyclotronSensitiveTarget.hh"
0032
0033 #include "G4AnalysisManager.hh"
0034 #include "G4RunManager.hh"
0035 #include "G4HCofThisEvent.hh"
0036 #include "G4UnitsTable.hh"
0037 #include "G4Step.hh"
0038 #include "G4SteppingManager.hh"
0039 #include "G4ThreeVector.hh"
0040 #include "G4SDManager.hh"
0041 #include "G4ios.hh"
0042 #include "G4SystemOfUnits.hh"
0043 #include "G4ThreeVector.hh"
0044 #include "G4Track.hh"
0045 #include "G4ParticleDefinition.hh"
0046 #include "G4DecayTable.hh"
0047 #include "G4VDecayChannel.hh"
0048 #include "G4TrackVector.hh"
0049 #include "G4VProcess.hh"
0050 #include "G4Tubs.hh"
0051 #include <map>
0052
0053 STCyclotronSensitiveTarget::STCyclotronSensitiveTarget(G4String name,
0054                                STCyclotronDetectorConstruction* det)
0055   : G4VSensitiveDetector(name),
0056     fDet(det)
0057 {
0058   fTempTrack = 0;
0059   fTempTrack1 = 0;
0060   fTempEnergy = 0.;
0061   fTempVector = G4ThreeVector(0.,0.,0.);
0062   fTrack=0;
0063 }
0064
0065 STCyclotronSensitiveTarget::~STCyclotronSensitiveTarget()
0066 {
0067  delete fTrack;
0068  }
0069
0070 G4bool STCyclotronSensitiveTarget::ProcessHits(G4Step* aStep, G4TouchableHistory*)
0071 {
0072
0073   STCyclotronRun* fRun = static_cast<STCyclotronRun*>(G4RunManager::GetRunManager()->GetNonConstCurrentRun());
0074   fTrack = aStep->GetTrack();
0075
0076   auto analysisManager = G4AnalysisManager::Instance();
0077
0078   //----------------------------------------------
0079   //        Volume info
0080   //----------------------------------------------
0081   G4double targetHalfDiameter= (fDet->GetTargetDiameter())/2.;
0082
0083   //----------------------------------------------
0084   //       Step information
0085   //----------------------------------------------
0086
0087   G4double edep = aStep->GetTotalEnergyDeposit();
0088   G4double energy = aStep->GetPreStepPoint()->GetKineticEnergy();
0089   G4ThreeVector momentumDirection = aStep->GetPreStepPoint()->GetMomentumDirection();
0090   G4ThreeVector vectorPosition = aStep->GetPreStepPoint()->GetPosition();
0091
0092   //----------------------------------------------
0093   //        Track
0094   //----------------------------------------------
0095
0096   G4ParticleDefinition* thePartDef = fTrack->GetDefinition();
0097   G4String partType= fTrack->GetDefinition()->GetParticleType();
0098   G4String name = fTrack->GetDefinition()->GetParticleName();
0099   G4double timeLife = fTrack->GetDefinition()->GetPDGLifeTime(); //<---
0100   const G4VProcess* process = fTrack->GetCreatorProcess();
0101
0102   //----------------------------------------------
0103   //    Collect general information concerning all of the particles
0104   //    Collect energy deposition ; separe decay case to beam case
0105   //----------------------------------------------
0106
0107   fRun->EnergyDepositionTarget(edep);
0108
0109
0110   //----------------------------------------------
0111   //Collect information about protons and deuterons
0112   //----------------------------------------------
0113
0114   if(name == "proton" || name == "deuteron")
0115     {
0116
0117       if(fTrack->GetTrackID()!=fTempTrack && (momentumDirection.getZ()>0.) &&
0118      vectorPosition.getX()<targetHalfDiameter &&
0119      vectorPosition.getX()>-targetHalfDiameter &&
0120      vectorPosition.getY()<targetHalfDiameter &&
0121      vectorPosition.getY()>-targetHalfDiameter)
0122     {
0123       analysisManager->FillH2(0,vectorPosition.getX(),vectorPosition.getY());
0124       analysisManager->FillH1(0,energy);
0125       fRun->CountParticlesTarget();
0126       fTempTrack = fTrack->GetTrackID();
0127     }
0128
0129       if(fTempTrack1 == 0)
0130     {
0131       fTempTrack1 = fTrack->GetTrackID();
0132     }
0133
0134       if(fTrack->GetTrackID()!=fTempTrack1 && (momentumDirection.getZ()>0.) &&
0135      fTempVector.getX()<targetHalfDiameter &&
0136      fTempVector.getX()>-targetHalfDiameter &&
0137      fTempVector.getY()<targetHalfDiameter &&
0138      fTempVector.getY()>-targetHalfDiameter )
0139     {
0140
0141       analysisManager->FillH2(4,fTempVector.getX(),fTempVector.getY());
0142       analysisManager->FillH1(2,fTempEnergy);
0143       fTempTrack1 = fTrack->GetTrackID();
0144     }
0145
0146       fTempVector =  aStep->GetPostStepPoint()->GetPosition();      //vectorPosition;
0147       fTempEnergy = aStep->GetPostStepPoint()->GetKineticEnergy();  //energy;
0148
0149       analysisManager->FillH2(3,vectorPosition.getZ(),energy);
0150
0151     }
0152
0153   //----------------------------------------------
0154   //    Store ID for particles that are
0155   //    not protons/electrons or deuterons
0156   //----------------------------------------------
0157
0158   if((name != "proton") && (name != "e-") && (name != "deuteron"))
0159     {
0160       fRun->StoreIsotopeID(fTrack->GetTrackID(),name);
0161     }
0162
0163   //----------------------------------------------
0164   //  Collect of information for unstable isotopes
0165   //  generated from an interaction with the target
0166   //----------------------------------------------
0167
0168   if (name!="deuteron")
0169     {
0170       if (( partType == "nucleus") &&  !(thePartDef->GetPDGStable()) && (fTrack->GetCurrentStepNumber()==1) && timeLife!=0.)
0171     {
0172       //G4cout << "Saving unstable particles ..." << G4endl;
0173       G4int Z=thePartDef->GetAtomicNumber();
0174       G4int A=thePartDef->GetAtomicMass();
0175       analysisManager->FillH2(2,Z,A);
0176
0177       //----------------------------------------------
0178       //   isotopes count
0179       //----------------------------------------------
0180
0181       fRun->PrimaryIsotopeCountTarget(name,timeLife);
0182       analysisManager->FillH1(4,fTrack->GetPosition().getZ());
0183       //particle that created the nucleus
0184       std::map<G4int,G4String> parentID = fRun->GetIsotopeID();
0185       G4String nameParent = parentID[fTrack->GetParentID()];
0186       fRun->ParticleParent(name, process->GetProcessName());
0187
0188       //G4cout << name << " : " << process->GetProcessName() << " with track ID " << fTrack->GetTrackID() << " and step ID " << fTrack->GetCurrentStepNumber() << G4endl;
0189
0190     }
0191     }
0192
0193   //----------------------------------------------
0194   //   Collect of information for stable isotopes
0195   //   generated from an interaction with the target
0196   //----------------------------------------------
0197
0198   if (name!="deuteron")
0199     {
0200       if (( partType == "nucleus") &&  (thePartDef->GetPDGStable()) && (process->GetProcessName() != "RadioactiveDecay")  && (fTrack->GetCurrentStepNumber()==1) )
0201     {
0202       //----------------------------------------------
0203       //    isotopes count
0204       //----------------------------------------------
0205       fRun->CountStableIsotopes(name);
0206     }
0207     }
0208
0209
0210   //----------------------------------------------
0211   //   Collect unstable isotopes from decay
0212   //----------------------------------------------
0213
0214   if (( partType == "nucleus") &&  !(thePartDef->GetPDGStable()) && (process->GetProcessName() == "RadioactiveDecay")  && (fTrack->GetCurrentStepNumber()==1) && timeLife!=0)
0215     {
0216       std::map<G4int,G4String>::iterator itbis;
0217       std::map<G4int,G4String> parentID = fRun->GetIsotopeID();
0218       G4String nameParent = parentID[fTrack->GetParentID()];
0219       fRun->DecayIsotopeCountTarget(name,nameParent,timeLife);
0220     }
0221
0222   //----------------------------------------------
0223   //   Collect any other particles emitted
0224   //----------------------------------------------
0225
0226   if((partType!="nucleus")&&(name!="proton")&&(name!="deuteron"))
0227     {
0228       fRun->ParticleCountTarget(name);
0229       //Condition so the particle will be counted for only one step
0230       if((fTrack->GetCurrentStepNumber()==1))
0231     {
0232       if(process->GetProcessName() != "RadioactiveDecay")
0233         {
0234           if(name=="e+"){
0235         analysisManager->FillH1(5,energy);
0236           }
0237           if(name=="e-"){
0238         analysisManager->FillH1(6,energy);
0239           }
0240           if(name=="gamma"){
0241         analysisManager->FillH1(7,energy);
0242           }
0243           if(name=="neutron"){
0244         analysisManager->FillH1(8,energy);
0245           }
0246         }
0247
0248       if(process->GetProcessName() == "RadioactiveDecay")
0249         {
0250           if(name=="e+"){
0251         analysisManager->FillH1(9,energy);
0252           }
0253           if(name=="e-"){
0254         analysisManager->FillH1(10,energy);
0255           }
0256           if(name=="gamma"){
0257         analysisManager->FillH1(11,energy);
0258           }
0259           if(name=="neutron"){
0260         analysisManager->FillH1(12,energy);
0261           }
0262           if(name=="nu_e"){
0263         analysisManager->FillH1(13,energy);
0264           }
0265           if(name=="anti_nu_e"){
0266         analysisManager->FillH1(14,energy);
0267           }
0268         }
0269     }
0270     }
0271
0272
0273   fRun->SetTargetVolume(fDet->GetTargetVolume());
0274   fRun->SetTargetThickness(fDet->GetTargetThickness());
0275   fRun->SetTargetDiameter(fDet->GetTargetDiameter());
0276
0277   return true;
0278
0279 }
0280
0281