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