Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/electromagnetic/TestEm5/src/TrackingAction.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /examples/extended/electromagnetic/TestEm5/src/TrackingAction.cc (Version 11.3.0) and /examples/extended/electromagnetic/TestEm5/src/TrackingAction.cc (Version 10.0.p3)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 /// \file electromagnetic/TestEm5/src/Tracking     26 /// \file electromagnetic/TestEm5/src/TrackingAction.cc
 27 /// \brief Implementation of the TrackingActio     27 /// \brief Implementation of the TrackingAction class
 28 //                                                 28 //
 29 //                                                 29 //
                                                   >>  30 // $Id: TrackingAction.cc 76464 2013-11-11 10:22:56Z gcosmo $
 30 //                                                 31 //
 31 //....oooOO0OOooo........oooOO0OOooo........oo     32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 32 //....oooOO0OOooo........oooOO0OOooo........oo     33 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 33                                                    34 
 34 #include "TrackingAction.hh"                       35 #include "TrackingAction.hh"
 35                                                    36 
 36 #include "DetectorConstruction.hh"                 37 #include "DetectorConstruction.hh"
                                                   >>  38 #include "Run.hh"
 37 #include "EventAction.hh"                          39 #include "EventAction.hh"
 38 #include "HistoManager.hh"                         40 #include "HistoManager.hh"
 39 #include "Run.hh"                              << 
 40                                                    41 
 41 #include "G4PhysicalConstants.hh"              << 
 42 #include "G4RunManager.hh"                         42 #include "G4RunManager.hh"
 43 #include "G4SystemOfUnits.hh"                  << 
 44 #include "G4Track.hh"                              43 #include "G4Track.hh"
                                                   >>  44 #include "G4PhysicalConstants.hh"
                                                   >>  45 #include "G4SystemOfUnits.hh"
 45                                                    46 
 46 //....oooOO0OOooo........oooOO0OOooo........oo     47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 47                                                    48 
 48 TrackingAction::TrackingAction(DetectorConstru <<  49 TrackingAction::TrackingAction(DetectorConstruction* DET, EventAction* EA)
 49   : fDetector(det), fEventAction(event)        <<  50 :G4UserTrackingAction(),fDetector(DET), fEventAction(EA)
 50 {}                                             <<  51 { }
 51                                                <<  52  
 52 //....oooOO0OOooo........oooOO0OOooo........oo     53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 53                                                    54 
 54 void TrackingAction::PreUserTrackingAction(con <<  55 void TrackingAction::PreUserTrackingAction(const G4Track* aTrack )
 55 {                                                  56 {
 56   // few initialisations                           57   // few initialisations
 57   //                                               58   //
 58   if (aTrack->GetTrackID() == 1) {                 59   if (aTrack->GetTrackID() == 1) {
 59     fXstartAbs = fDetector->GetxstartAbs();        60     fXstartAbs = fDetector->GetxstartAbs();
 60     fXendAbs = fDetector->GetxendAbs();        <<  61     fXendAbs   = fDetector->GetxendAbs();
 61     fPrimaryCharge = aTrack->GetDefinition()->     62     fPrimaryCharge = aTrack->GetDefinition()->GetPDGCharge();
 62     fDirX = aTrack->GetMomentumDirection().x() << 
 63   }                                                63   }
 64 }                                                  64 }
 65                                                    65 
 66 //....oooOO0OOooo........oooOO0OOooo........oo     66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 67                                                    67 
 68 void TrackingAction::PostUserTrackingAction(co     68 void TrackingAction::PostUserTrackingAction(const G4Track* aTrack)
 69 {                                                  69 {
 70   G4AnalysisManager* analysisManager = G4Analy     70   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
 71                                                <<  71     
 72   Run* run = static_cast<Run*>(G4RunManager::G <<  72   Run* run = static_cast<Run*>(
                                                   >>  73              G4RunManager::GetRunManager()->GetNonConstCurrentRun()); 
 73                                                    74 
 74   G4ThreeVector position = aTrack->GetPosition     75   G4ThreeVector position = aTrack->GetPosition();
 75   G4ThreeVector vertex = aTrack->GetVertexPosi <<  76   G4ThreeVector vertex   = aTrack->GetVertexPosition();  
 76   G4double charge = aTrack->GetDefinition()->G <<  77   G4double charge        = aTrack->GetDefinition()->GetPDGCharge();
 77   G4double energy = aTrack->GetKineticEnergy() <<  78 
 78                                                <<  79   G4bool transmit = ((position.x() >= fXendAbs) && (vertex.x() < fXendAbs));
 79   G4bool transmit =                            <<  80   G4bool reflect  =  (position.x() <= fXstartAbs);
 80     (((fDirX >= 0.0 && position.x() >= fXendAb << 
 81      && energy > 0.0);                         << 
 82   G4bool reflect =                             << 
 83     (((fDirX >= 0.0 && position.x() <= fXstart << 
 84      && energy > 0.0);                         << 
 85   G4bool notabsor = (transmit || reflect);         81   G4bool notabsor = (transmit || reflect);
 86                                                    82 
 87   // transmitted + reflected particles counter <<  83   //transmitted + reflected particles counter
 88   //                                               84   //
 89   G4int flag = 0;                                  85   G4int flag = 0;
 90   if (charge == fPrimaryCharge) flag = 1;      <<  86   if (charge == fPrimaryCharge)  flag = 1;
 91   if (aTrack->GetTrackID() == 1) flag = 2;         87   if (aTrack->GetTrackID() == 1) flag = 2;
 92   if (transmit) fEventAction->SetTransmitFlag(     88   if (transmit) fEventAction->SetTransmitFlag(flag);
 93   if (reflect) fEventAction->SetReflectFlag(fl <<  89   if (reflect)  fEventAction->SetReflectFlag(flag);
 94                                                <<  90       
 95   //                                               91   //
 96   // histograms                                <<  92   //histograms
 97   //                                               93   //
 98   G4bool charged = (charge != 0.);             <<  94   G4bool charged  = (charge != 0.);
 99   G4bool neutral = !charged;                       95   G4bool neutral = !charged;
100                                                    96 
101   // energy spectrum at exit                   <<  97   //energy spectrum at exit
102   // zero energy charged particle are not take << 
103   //                                               98   //
104   G4int id = 0;                                <<  99   G4int id = 0;  
105   if (transmit && charged)                     << 100        if (transmit && charged) id = 10;
106     id = 10;                                   << 101   else if (transmit && neutral) id = 20;
107   else if (transmit && neutral)                << 102   else if (reflect  && charged) id = 30;
108     id = 20;                                   << 103   else if (reflect  && neutral) id = 40;
109   else if (reflect && charged)                 << 104 
110     id = 30;                                   << 105   if (id>0)   
111   else if (reflect && neutral)                 << 106   analysisManager->FillH1(id, aTrack->GetKineticEnergy());
112     id = 40;                                   << 107     
113                                                << 108   //energy leakage
114   if (id > 0) {                                << 
115     analysisManager->FillH1(id, energy);       << 
116   }                                            << 
117                                                << 
118   // energy leakage                            << 
119   //                                              109   //
120   if (notabsor) {                                 110   if (notabsor) {
121     G4int trackID = aTrack->GetTrackID();         111     G4int trackID = aTrack->GetTrackID();
122     G4int index = 0;                           << 112     G4int index = 0; if (trackID > 1) index = 1;    //primary=0, secondaries=1
123     if (trackID > 1) index = 1;  // primary=0, << 
124     G4double eleak = aTrack->GetKineticEnergy(    113     G4double eleak = aTrack->GetKineticEnergy();
125     if ((aTrack->GetDefinition() == G4Positron    114     if ((aTrack->GetDefinition() == G4Positron::Positron()) && (trackID > 1))
126       eleak += 2 * electron_mass_c2;           << 115       eleak += 2*electron_mass_c2;
127     run->AddEnergyLeak(eleak, index);          << 116     run->AddEnergyLeak(eleak,index);  
128   }                                               117   }
129                                                   118 
130   // space angle distribution at exit : dN/dOm << 119   //space angle distribution at exit : dN/dOmega
131   //                                              120   //
132   G4ThreeVector direction = aTrack->GetMomentu << 121   G4ThreeVector direction = aTrack->GetMomentumDirection();    
133   id = 0;                                      << 122   id = 0;   
134   if (transmit && charged)                     << 123        if (transmit && charged) id = 12;
135     id = 12;                                   << 124   else if (transmit && neutral) id = 22;
136   else if (transmit && neutral)                << 125   else if (reflect  && charged) id = 32;
137     id = 22;                                   << 126   else if (reflect  && neutral) id = 42;
138   else if (reflect && charged)                 << 
139     id = 32;                                   << 
140   else if (reflect && neutral)                 << 
141     id = 42;                                   << 
142                                                   127 
143   if (id > 0 && std::abs(direction.x()) < 1.0) << 128   if (id>0) {
144     G4double theta = std::acos(direction.x()); << 129     G4double theta  = std::acos(direction.x());
145     if (theta > 0.0) {                            130     if (theta > 0.0) {
146       G4double dteta = analysisManager->GetH1W << 131       G4double dteta  = analysisManager->GetH1Width(id);
147       G4double unit = analysisManager->GetH1Un << 132       G4double unit   = analysisManager->GetH1Unit(id);    
148       G4double weight = unit / (twopi * std::s << 133       G4double weight = (unit*unit)/(twopi*std::sin(theta)*dteta);
149       /*                                          134       /*
150       G4cout << "theta, dteta, unit, weight: " << 135       G4cout << "theta, dteta, unit, weight: " 
151              << theta << "  "                  << 136              << theta << "  "   
152              << dteta << "  "                  << 137              << dteta << "  "   
153              << unit << "  "                   << 138              << unit << "  "   
154              << weight << G4endl;              << 139              << weight << G4endl;   
155       */                                          140       */
156       analysisManager->FillH1(id, theta, weigh << 141       analysisManager->FillH1(id,theta,weight);
157     }                                             142     }
158   }                                               143   }
159                                                << 144   
160   // energy fluence at exit : dE(MeV)/dOmega   << 145   //energy fluence at exit : dE(MeV)/dOmega
161   //                                              146   //
162   id = 0;                                      << 147   id = 0;  
163   if (transmit && charged)                     << 148   if (transmit && charged) id = 11;
164     id = 11;                                   << 149   else if (transmit && neutral) id = 21;
165   else if (reflect && charged)                 << 150   else if (reflect  && charged) id = 31;
166     id = 31;                                   << 151   else if (reflect  && neutral) id = 41;
167   else if (transmit && neutral)                << 
168     id = 21;                                   << 
169   else if (reflect && neutral)                 << 
170     id = 41;                                   << 
171                                                   152 
172   if (id > 0 && std::abs(direction.x()) < 1.0) << 153   if (id>0) {
173     G4double theta = std::acos(direction.x()); << 154     G4double theta  = std::acos(direction.x());
174     if (theta > 0.0) {                            155     if (theta > 0.0) {
175       G4double dteta = analysisManager->GetH1W << 156       G4double dteta  = analysisManager->GetH1Width(id);
176       G4double unit = analysisManager->GetH1Un << 157       G4double unit   = analysisManager->GetH1Unit(id);    
177       G4double weight = unit / (twopi * std::s << 158       G4double weight = (unit*unit)/(twopi*std::sin(theta)*dteta);
178       weight *= (aTrack->GetKineticEnergy() /  << 159       weight *= (aTrack->GetKineticEnergy()/MeV); 
179       analysisManager->FillH1(id, theta, weigh << 160       analysisManager->FillH1(id,theta,weight);
180     }                                          << 161     } 
181   }                                            << 162   }
182                                                << 163   
183   // projected angles distribution at exit     << 164   //projected angles distribution at exit
184   //                                           << 165   //
185   id = 0;                                      << 166   id = 0;   
186   if (transmit && charged)                     << 167   if (transmit && charged) id = 13;
187     id = 13;                                   << 168   else if (transmit && neutral) id = 23;
188   else if (transmit && neutral)                << 169   else if (reflect  && charged) id = 33;
189     id = 23;                                   << 170   else if (reflect  && neutral) id = 43;
190   else if (reflect && charged)                 << 
191     id = 33;                                   << 
192   else if (reflect && neutral)                 << 
193     id = 43;                                   << 
194                                                   171 
195   if (id > 0) {                                << 172   if (id>0) {
196     if (direction.x() != 0.0) {                   173     if (direction.x() != 0.0) {
197       G4double tet = std::atan(direction.y() / << 174       G4double tet = std::atan(direction.y()/std::fabs(direction.x()));
198       analysisManager->FillH1(id, tet);        << 175       analysisManager->FillH1(id,tet);
199       if (transmit && (flag == 2)) run->AddMsc    176       if (transmit && (flag == 2)) run->AddMscProjTheta(tet);
200                                                   177 
201       tet = std::atan(direction.z() / std::fab << 178       tet = std::atan(direction.z()/std::fabs(direction.x()));
202       analysisManager->FillH1(id, tet);        << 179       analysisManager->FillH1(id,tet);
203       if (transmit && (flag == 2)) run->AddMsc    180       if (transmit && (flag == 2)) run->AddMscProjTheta(tet);
204     }                                             181     }
205   }                                               182   }
206                                                   183 
207   // projected position and radius at exit     << 184   //projected position and radius at exit
208   //                                              185   //
209   if (transmit && energy > 0.0) {              << 186   id = 0;   
                                                   >> 187   if (transmit && charged) id = 14;
                                                   >> 188   
                                                   >> 189   if (id>0) {
210     G4double y = position.y(), z = position.z(    190     G4double y = position.y(), z = position.z();
211     G4double r = std::sqrt(y * y + z * z);     << 191     G4double r = std::sqrt(y*y + z*z);
212     analysisManager->FillH1(14, y);            << 192     analysisManager->FillH1(id,   y);
213     analysisManager->FillH1(14, z);            << 193     analysisManager->FillH1(id,   z);
214     analysisManager->FillH1(15, r);            << 194     analysisManager->FillH1(id+1, r);
215   }                                               195   }
216                                                << 196   
217   // x-vertex of charged secondaries           << 197   //x-vertex of charged secondaries
218   //                                              198   //
219   if ((aTrack->GetParentID() == 1) && charged)    199   if ((aTrack->GetParentID() == 1) && charged) {
220     G4double xVertex = (aTrack->GetVertexPosit    200     G4double xVertex = (aTrack->GetVertexPosition()).x();
221     analysisManager->FillH1(6, xVertex);          201     analysisManager->FillH1(6, xVertex);
222     if (notabsor) analysisManager->FillH1(7, x << 202     if (notabsor) analysisManager->FillH1(7, xVertex); 
223   }                                               203   }
224 }                                                 204 }
225                                                   205 
226 //....oooOO0OOooo........oooOO0OOooo........oo    206 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >> 207 
227                                                   208