TestEm5/src/TrackingAction.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 /// \file TrackingAction.cc
0027 /// \brief Implementation of the TrackingAction class
0028
0029 #include "TrackingAction.hh"
0030
0031 #include "DetectorConstruction.hh"
0032 #include "EventAction.hh"
0033 #include "HistoManager.hh"
0034 #include "Run.hh"
0035
0036 #include "G4PhysicalConstants.hh"
0037 #include "G4RunManager.hh"
0038 #include "G4SystemOfUnits.hh"
0039 #include "G4Track.hh"
0040
0041 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0042
0043 TrackingAction::TrackingAction(DetectorConstruction* det, EventAction* event)
0044   : fDetector(det), fEventAction(event)
0045 {}
0046
0047 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0048
0049 void TrackingAction::PreUserTrackingAction(const G4Track* aTrack)
0050 {
0051   // few initialisations
0052   //
0053   if (aTrack->GetTrackID() == 1) {
0054     fXstartAbs = fDetector->GetxstartAbs();
0055     fXendAbs = fDetector->GetxendAbs();
0056     fPrimaryCharge = aTrack->GetDefinition()->GetPDGCharge();
0057     fDirX = aTrack->GetMomentumDirection().x();
0058   }
0059 }
0060
0061 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0062
0063 void TrackingAction::PostUserTrackingAction(const G4Track* aTrack)
0064 {
0065   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0066
0067   Run* run = static_cast<Run*>(G4RunManager::GetRunManager()->GetNonConstCurrentRun());
0068
0069   G4ThreeVector position = aTrack->GetPosition();
0070   G4ThreeVector vertex = aTrack->GetVertexPosition();
0071   G4double charge = aTrack->GetDefinition()->GetPDGCharge();
0072   G4double energy = aTrack->GetKineticEnergy();
0073
0074   G4bool transmit =
0075     (((fDirX >= 0.0 && position.x() >= fXendAbs) || (fDirX < 0.0 && position.x() <= fXstartAbs))
0076      && energy > 0.0);
0077   G4bool reflect =
0078     (((fDirX >= 0.0 && position.x() <= fXstartAbs) || (fDirX < 0.0 && position.x() >= fXendAbs))
0079      && energy > 0.0);
0080   G4bool notabsor = (transmit || reflect);
0081
0082   // transmitted + reflected particles counter
0083   //
0084   G4int flag = 0;
0085   if (charge == fPrimaryCharge) flag = 1;
0086   if (aTrack->GetTrackID() == 1) flag = 2;
0087   if (transmit) fEventAction->SetTransmitFlag(flag);
0088   if (reflect) fEventAction->SetReflectFlag(flag);
0089
0090   //
0091   // histograms
0092   //
0093   G4bool charged = (charge != 0.);
0094   G4bool neutral = !charged;
0095
0096   // energy spectrum at exit
0097   // zero energy charged particle are not taken into account
0098   //
0099   G4int id = 0;
0100   if (transmit && charged)
0101     id = 10;
0102   else if (transmit && neutral)
0103     id = 20;
0104   else if (reflect && charged)
0105     id = 30;
0106   else if (reflect && neutral)
0107     id = 40;
0108
0109   if (id > 0) {
0110     analysisManager->FillH1(id, energy);
0111   }
0112
0113   // energy leakage
0114   //
0115   if (notabsor) {
0116     G4int trackID = aTrack->GetTrackID();
0117     G4int index = 0;
0118     if (trackID > 1) index = 1;  // primary=0, secondaries=1
0119     G4double eleak = aTrack->GetKineticEnergy();
0120     if ((aTrack->GetDefinition() == G4Positron::Positron()) && (trackID > 1))
0121       eleak += 2 * electron_mass_c2;
0122     run->AddEnergyLeak(eleak, index);
0123   }
0124
0125   // space angle distribution at exit : dN/dOmega
0126   //
0127   G4ThreeVector direction = aTrack->GetMomentumDirection();
0128   id = 0;
0129   if (transmit && charged)
0130     id = 12;
0131   else if (transmit && neutral)
0132     id = 22;
0133   else if (reflect && charged)
0134     id = 32;
0135   else if (reflect && neutral)
0136     id = 42;
0137
0138   if (id > 0 && std::abs(direction.x()) < 1.0) {
0139     G4double theta = std::acos(direction.x());
0140     if (theta > 0.0) {
0141       G4double dteta = analysisManager->GetH1Width(id);
0142       G4double unit = analysisManager->GetH1Unit(id);
0143       G4double weight = unit / (twopi * std::sin(theta) * dteta);
0144       /*
0145       G4cout << "theta, dteta, unit, weight: "
0146              << theta << "  "
0147              << dteta << "  "
0148              << unit << "  "
0149              << weight << G4endl;
0150       */
0151       analysisManager->FillH1(id, theta, weight);
0152     }
0153   }
0154
0155   // energy fluence at exit : dE(MeV)/dOmega
0156   //
0157   id = 0;
0158   if (transmit && charged)
0159     id = 11;
0160   else if (reflect && charged)
0161     id = 31;
0162   else if (transmit && neutral)
0163     id = 21;
0164   else if (reflect && neutral)
0165     id = 41;
0166
0167   if (id > 0 && std::abs(direction.x()) < 1.0) {
0168     G4double theta = std::acos(direction.x());
0169     if (theta > 0.0) {
0170       G4double dteta = analysisManager->GetH1Width(id);
0171       G4double unit = analysisManager->GetH1Unit(id);
0172       G4double weight = unit / (twopi * std::sin(theta) * dteta);
0173       weight *= (aTrack->GetKineticEnergy() / MeV);
0174       analysisManager->FillH1(id, theta, weight);
0175     }
0176   }
0177
0178   // projected angles distribution at exit
0179   //
0180   id = 0;
0181   if (transmit && charged)
0182     id = 13;
0183   else if (transmit && neutral)
0184     id = 23;
0185   else if (reflect && charged)
0186     id = 33;
0187   else if (reflect && neutral)
0188     id = 43;
0189
0190   if (id > 0) {
0191     if (direction.x() != 0.0) {
0192       G4double tet = std::atan(direction.y() / std::fabs(direction.x()));
0193       analysisManager->FillH1(id, tet);
0194       if (transmit && (flag == 2)) run->AddMscProjTheta(tet);
0195
0196       tet = std::atan(direction.z() / std::fabs(direction.x()));
0197       analysisManager->FillH1(id, tet);
0198       if (transmit && (flag == 2)) run->AddMscProjTheta(tet);
0199     }
0200   }
0201
0202   // projected position and radius at exit
0203   //
0204   if (transmit && energy > 0.0) {
0205     G4double y = position.y(), z = position.z();
0206     G4double r = std::sqrt(y * y + z * z);
0207     analysisManager->FillH1(14, y);
0208     analysisManager->FillH1(14, z);
0209     analysisManager->FillH1(15, r);
0210   }
0211
0212   // x-vertex of charged secondaries
0213   //
0214   if ((aTrack->GetParentID() == 1) && charged) {
0215     G4double xVertex = (aTrack->GetVertexPosition()).x();
0216     analysisManager->FillH1(6, xVertex);
0217     if (notabsor) analysisManager->FillH1(7, xVertex);
0218   }
0219 }
0220
0221 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......