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