| 1 | //
|
|---|
| 2 | // ********************************************************************
|
|---|
| 3 | // * License and Disclaimer *
|
|---|
| 4 | // * *
|
|---|
| 5 | // * The Geant4 software is copyright of the Copyright Holders of *
|
|---|
| 6 | // * the Geant4 Collaboration. It is provided under the terms and *
|
|---|
| 7 | // * conditions of the Geant4 Software License, included in the file *
|
|---|
| 8 | // * LICENSE and available at http://cern.ch/geant4/license . These *
|
|---|
| 9 | // * include a list of copyright holders. *
|
|---|
| 10 | // * *
|
|---|
| 11 | // * Neither the authors of this software system, nor their employing *
|
|---|
| 12 | // * institutes,nor the agencies providing financial support for this *
|
|---|
| 13 | // * work make any representation or warranty, express or implied, *
|
|---|
| 14 | // * regarding this software system or assume any liability for its *
|
|---|
| 15 | // * use. Please see the license in the file LICENSE and URL above *
|
|---|
| 16 | // * for the full disclaimer and the limitation of liability. *
|
|---|
| 17 | // * *
|
|---|
| 18 | // * This code implementation is the result of the scientific and *
|
|---|
| 19 | // * technical work of the GEANT4 collaboration. *
|
|---|
| 20 | // * By using, copying, modifying or distributing the software (or *
|
|---|
| 21 | // * any work based on the software) you agree to acknowledge its *
|
|---|
| 22 | // * use in resulting scientific publications, and indicate your *
|
|---|
| 23 | // * acceptance of all terms of the Geant4 Software license. *
|
|---|
| 24 | // ********************************************************************
|
|---|
| 25 | //
|
|---|
| 26 | //
|
|---|
| 27 | // $Id: TrackingAction.cc,v 1.18 2008/08/28 15:28:04 maire Exp $
|
|---|
| 28 | // GEANT4 tag $Name: geant4-09-04-beta-01 $
|
|---|
| 29 | //
|
|---|
| 30 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|---|
| 31 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|---|
| 32 |
|
|---|
| 33 | #include "TrackingAction.hh"
|
|---|
| 34 |
|
|---|
| 35 | #include "DetectorConstruction.hh"
|
|---|
| 36 | #include "RunAction.hh"
|
|---|
| 37 | #include "EventAction.hh"
|
|---|
| 38 | #include "HistoManager.hh"
|
|---|
| 39 |
|
|---|
| 40 | #include "G4Track.hh"
|
|---|
| 41 |
|
|---|
| 42 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|---|
| 43 |
|
|---|
| 44 | TrackingAction::TrackingAction(DetectorConstruction* DET,RunAction* RA,
|
|---|
| 45 | EventAction* EA, HistoManager* HM)
|
|---|
| 46 | :detector(DET), runaction(RA), eventaction(EA), histoManager(HM)
|
|---|
| 47 | { }
|
|---|
| 48 |
|
|---|
| 49 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|---|
| 50 |
|
|---|
| 51 | void TrackingAction::PreUserTrackingAction(const G4Track* aTrack )
|
|---|
| 52 | {
|
|---|
| 53 | // few initialisations
|
|---|
| 54 | //
|
|---|
| 55 | if (aTrack->GetTrackID() == 1) {
|
|---|
| 56 | xstartAbs = detector->GetxstartAbs();
|
|---|
| 57 | xendAbs = detector->GetxendAbs();
|
|---|
| 58 | primaryCharge = aTrack->GetDefinition()->GetPDGCharge();
|
|---|
| 59 | }
|
|---|
| 60 | }
|
|---|
| 61 |
|
|---|
| 62 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|---|
| 63 |
|
|---|
| 64 | void TrackingAction::PostUserTrackingAction(const G4Track* aTrack)
|
|---|
| 65 | {
|
|---|
| 66 | G4ThreeVector position = aTrack->GetPosition();
|
|---|
| 67 | G4ThreeVector vertex = aTrack->GetVertexPosition();
|
|---|
| 68 | G4double charge = aTrack->GetDefinition()->GetPDGCharge();
|
|---|
| 69 |
|
|---|
| 70 | G4bool transmit = ((position.x() >= xendAbs) && (vertex.x() < xendAbs));
|
|---|
| 71 | G4bool reflect = (position.x() <= xstartAbs);
|
|---|
| 72 | G4bool notabsor = (transmit || reflect);
|
|---|
| 73 |
|
|---|
| 74 | //transmitted + reflected particles counter
|
|---|
| 75 | //
|
|---|
| 76 | G4int flag = 0;
|
|---|
| 77 | if (charge == primaryCharge) flag = 1;
|
|---|
| 78 | if (aTrack->GetTrackID() == 1) flag = 2;
|
|---|
| 79 | if (transmit) eventaction->SetTransmitFlag(flag);
|
|---|
| 80 | if (reflect) eventaction->SetReflectFlag(flag);
|
|---|
| 81 |
|
|---|
| 82 | //
|
|---|
| 83 | //histograms
|
|---|
| 84 | //
|
|---|
| 85 | G4bool charged = (charge != 0.);
|
|---|
| 86 | G4bool neutral = !charged;
|
|---|
| 87 |
|
|---|
| 88 | //energy spectrum at exit
|
|---|
| 89 | //
|
|---|
| 90 | G4int id = 0;
|
|---|
| 91 | if (transmit && charged) id = 10;
|
|---|
| 92 | else if (transmit && neutral) id = 20;
|
|---|
| 93 | else if (reflect && charged) id = 30;
|
|---|
| 94 | else if (reflect && neutral) id = 40;
|
|---|
| 95 |
|
|---|
| 96 | if (id>0)
|
|---|
| 97 | histoManager->FillHisto(id, aTrack->GetKineticEnergy());
|
|---|
| 98 |
|
|---|
| 99 | //energy leakage
|
|---|
| 100 | //
|
|---|
| 101 | if (notabsor) {
|
|---|
| 102 | G4int trackID = aTrack->GetTrackID();
|
|---|
| 103 | G4int index = 0; if (trackID > 1) index = 1; //primary=0, secondaries=1
|
|---|
| 104 | G4double eleak = aTrack->GetKineticEnergy();
|
|---|
| 105 | if ((aTrack->GetDefinition() == G4Positron::Positron()) && (trackID > 1))
|
|---|
| 106 | eleak += 2*electron_mass_c2;
|
|---|
| 107 | runaction->AddEnergyLeak(eleak,index);
|
|---|
| 108 | }
|
|---|
| 109 |
|
|---|
| 110 | //space angle distribution at exit : dN/dOmega
|
|---|
| 111 | //
|
|---|
| 112 | G4ThreeVector direction = aTrack->GetMomentumDirection();
|
|---|
| 113 | id = 0;
|
|---|
| 114 | if (transmit && charged) id = 12;
|
|---|
| 115 | else if (transmit && neutral) id = 22;
|
|---|
| 116 | else if (reflect && charged) id = 32;
|
|---|
| 117 | else if (reflect && neutral) id = 42;
|
|---|
| 118 |
|
|---|
| 119 | if (id>0) {
|
|---|
| 120 | G4double theta = std::acos(direction.x());
|
|---|
| 121 | G4double dteta = histoManager->GetBinWidth(id);
|
|---|
| 122 | G4double unit = histoManager->GetHistoUnit(id);
|
|---|
| 123 | G4double weight = (unit*unit)/(twopi*std::sin(theta)*dteta);
|
|---|
| 124 | histoManager->FillHisto(id,theta,weight);
|
|---|
| 125 | }
|
|---|
| 126 |
|
|---|
| 127 | //energy fluence at exit : dE(MeV)/dOmega
|
|---|
| 128 | //
|
|---|
| 129 | id = 0;
|
|---|
| 130 | if (transmit && charged) id = 11;
|
|---|
| 131 | else if (transmit && neutral) id = 21;
|
|---|
| 132 | else if (reflect && charged) id = 31;
|
|---|
| 133 | else if (reflect && neutral) id = 41;
|
|---|
| 134 |
|
|---|
| 135 | if (id>0) {
|
|---|
| 136 | G4double theta = std::acos(direction.x());
|
|---|
| 137 | G4double dteta = histoManager->GetBinWidth(id);
|
|---|
| 138 | G4double unit = histoManager->GetHistoUnit(id);
|
|---|
| 139 | G4double weight = (unit*unit)/(twopi*std::sin(theta)*dteta);
|
|---|
| 140 | weight *= (aTrack->GetKineticEnergy()/MeV);
|
|---|
| 141 | histoManager->FillHisto(id,theta,weight);
|
|---|
| 142 | }
|
|---|
| 143 |
|
|---|
| 144 | //projected angles distribution at exit
|
|---|
| 145 | //
|
|---|
| 146 | id = 0;
|
|---|
| 147 | if (transmit && charged) id = 13;
|
|---|
| 148 | else if (transmit && neutral) id = 23;
|
|---|
| 149 | else if (reflect && charged) id = 33;
|
|---|
| 150 | else if (reflect && neutral) id = 43;
|
|---|
| 151 |
|
|---|
| 152 | if(id>0) {
|
|---|
| 153 | G4double tet = std::atan(direction.y()/std::fabs(direction.x()));
|
|---|
| 154 | histoManager->FillHisto(id,tet);
|
|---|
| 155 | if (transmit && (flag == 2)) runaction->AddMscProjTheta(tet);
|
|---|
| 156 |
|
|---|
| 157 | tet = std::atan(direction.z()/std::fabs(direction.x()));
|
|---|
| 158 | histoManager->FillHisto(id,tet);
|
|---|
| 159 | if (transmit && (flag == 2)) runaction->AddMscProjTheta(tet);
|
|---|
| 160 | }
|
|---|
| 161 |
|
|---|
| 162 | //projected position and radius at exit
|
|---|
| 163 | //
|
|---|
| 164 | id = 0;
|
|---|
| 165 | if (transmit && charged) id = 14;
|
|---|
| 166 |
|
|---|
| 167 | if (id>0) {
|
|---|
| 168 | G4double y = position.y(), z = position.z();
|
|---|
| 169 | G4double r = std::sqrt(y*y + z*z);
|
|---|
| 170 | histoManager->FillHisto(id, y);
|
|---|
| 171 | histoManager->FillHisto(id, z);
|
|---|
| 172 | histoManager->FillHisto(id+1, r);
|
|---|
| 173 | }
|
|---|
| 174 |
|
|---|
| 175 | //x-vertex of charged secondaries
|
|---|
| 176 | //
|
|---|
| 177 | if ((aTrack->GetParentID() == 1) && charged) {
|
|---|
| 178 | G4double xVertex = (aTrack->GetVertexPosition()).x();
|
|---|
| 179 | histoManager->FillHisto(4, xVertex);
|
|---|
| 180 | if (notabsor) histoManager->FillHisto(5, xVertex);
|
|---|
| 181 | }
|
|---|
| 182 | }
|
|---|
| 183 |
|
|---|
| 184 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|---|
| 185 |
|
|---|