Geant4 Cross Reference |
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 // $Id: SteppingAction.cc 73008 2013-08-15 08:43:13Z gcosmo $ 26 // 27 // 27 //....oooOO0OOooo........oooOO0OOooo........oo 28 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 28 //....oooOO0OOooo........oooOO0OOooo........oo 29 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 29 30 30 #include "SteppingAction.hh" 31 #include "SteppingAction.hh" 31 32 32 #include "DetectorConstruction.hh" 33 #include "DetectorConstruction.hh" >> 34 #include "RunAction.hh" >> 35 #include "PrimaryGeneratorAction.hh" 33 #include "EventAction.hh" 36 #include "EventAction.hh" 34 37 35 #include "G4Step.hh" 38 #include "G4Step.hh" 36 39 >> 40 #include "G4Geantino.hh" >> 41 37 //....oooOO0OOooo........oooOO0OOooo........oo 42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 38 43 39 SteppingAction::SteppingAction(DetectorConstru 44 SteppingAction::SteppingAction(DetectorConstruction* det, EventAction* evt) 40 :G4UserSteppingAction(),detector(det),eventAct 45 :G4UserSteppingAction(),detector(det),eventAct(evt) 41 { 46 { 42 first = true; 47 first = true; 43 lvol_world = lvol_module = lvol_layer = lvol 48 lvol_world = lvol_module = lvol_layer = lvol_fiber = 0; 44 } 49 } 45 50 46 //....oooOO0OOooo........oooOO0OOooo........oo 51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 47 52 48 SteppingAction::~SteppingAction() 53 SteppingAction::~SteppingAction() 49 {} 54 {} 50 55 51 //....oooOO0OOooo........oooOO0OOooo........oo 56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 52 57 53 void SteppingAction::UserSteppingAction(const 58 void SteppingAction::UserSteppingAction(const G4Step* step ) 54 { 59 { 55 //some initialisation 60 //some initialisation 56 // 61 // 57 if (first) { 62 if (first) { 58 lvol_world = detector->GetLvolWorld(); 63 lvol_world = detector->GetLvolWorld(); 59 lvol_module = detector->GetLvolModule(); 64 lvol_module = detector->GetLvolModule(); 60 lvol_layer = detector->GetLvolLayer(); 65 lvol_layer = detector->GetLvolLayer(); 61 lvol_fiber = detector->GetLvolFiber(); << 66 lvol_fiber = detector->GetLvolFiber(); >> 67 >> 68 calorThickness = detector->GetCalorThickness(); >> 69 calorSizeYZ = detector->GetCalorSizeYZ(); >> 70 moduleThickness = detector->GetModuleThickness(); >> 71 d1Pixel = detector->GetD1Pixel(); >> 72 d2Pixel = detector->GetD2Pixel(); >> 73 n1pxl = detector->GetN1Pixels(); >> 74 n2pxl = detector->GetN2Pixels(); >> 75 n1shift = detector->GetN1Shift(); >> 76 62 first = false; 77 first = false; 63 } 78 } >> 79 >> 80 //locate point in geometry >> 81 // >> 82 G4TouchableHandle touch1 = step->GetPreStepPoint()->GetTouchableHandle(); >> 83 G4LogicalVolume* lvol = touch1->GetVolume()->GetLogicalVolume(); >> 84 >> 85 //if world, return >> 86 // >> 87 if (lvol == lvol_world) return; 64 88 >> 89 //sum nb of radiation length of calorimeter with geantino >> 90 // >> 91 G4ParticleDefinition* particle = step->GetTrack()->GetDefinition(); >> 92 if (particle == G4Geantino::Geantino()) { >> 93 G4double radl = lvol->GetMaterial()->GetRadlen(); >> 94 G4double stepl = step->GetStepLength(); >> 95 eventAct->SumNbRadLength(stepl/radl); >> 96 } >> 97 65 //if no edep, return 98 //if no edep, return 66 // 99 // 67 G4double edep = step->GetTotalEnergyDeposit() 100 G4double edep = step->GetTotalEnergyDeposit(); 68 ///if (edep == 0.) return; << 101 if (edep == 0.) return; 69 102 70 //locate point in geometry << 103 //locate position and compute pixel number 71 // 104 // 72 G4int iModule = 0; << 105 G4ThreeVector point1 = step->GetPreStepPoint()->GetPosition(); 73 G4int iLayer = 0; << 106 G4int i1Module = (int) ((point1.x() + 0.5*calorThickness)/moduleThickness); 74 G4int iFiber = 0; << 107 G4int i1Pixel = (int) ((point1.x() + 0.5*calorThickness)/d1Pixel); 75 << 108 if (i1Pixel < 0 ) i1Pixel = 0; 76 G4TouchableHandle touch1 = step->GetPreStepPo << 109 if (i1Pixel >= n1pxl) i1Pixel = n1pxl - 1; 77 G4LogicalVolume* lvol = touch1->GetVolume()-> << 110 G4double point1yz = point1.y(); >> 111 if (i1Module%2 != 0) point1yz = point1.z(); >> 112 G4int i2Pixel = (int) ((point1yz + 0.5*calorSizeYZ)/d2Pixel); >> 113 if (i2Pixel < 0 ) i2Pixel = 0; >> 114 if (i2Pixel >= n2pxl) i2Pixel = n2pxl - 1; 78 115 79 if (lvol == lvol_world) return; << 116 G4int iPixel = i1Pixel*n1shift + i2Pixel; 80 else if (lvol == lvol_module) { iModule = tou << 117 81 else if (lvol == lvol_layer) { iLayer = tou << 118 // sum total energy deposit 82 iModule = tou << 83 else if (lvol == lvol_fiber) { iFiber = tou << 84 iLayer = tou << 85 iModule = tou << 86 << 87 // sum edep << 88 // 119 // 89 eventAct->SumDeStep(iModule, iLayer, iFiber, << 120 eventAct->SumTotalEnergy(iPixel, edep); >> 121 >> 122 //in fiber ? >> 123 // >> 124 if (lvol == lvol_fiber) { >> 125 eventAct->SumVisibleEnergy(iPixel, edep); >> 126 } 90 } 127 } 91 128 92 //....oooOO0OOooo........oooOO0OOooo........oo 129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 93 130 94 G4double SteppingAction::BirksAttenuation(cons 131 G4double SteppingAction::BirksAttenuation(const G4Step* aStep) 95 { 132 { 96 //Example of Birk attenuation law in organic 133 //Example of Birk attenuation law in organic scintillators. 97 //adapted from Geant3 PHYS337. See MIN 80 (19 134 //adapted from Geant3 PHYS337. See MIN 80 (1970) 239-244 98 // 135 // 99 const G4Material* material = aStep->GetTrack( << 136 G4Material* material = aStep->GetTrack()->GetMaterial(); 100 G4double birk1 = material->GetIonisatio 137 G4double birk1 = material->GetIonisation()->GetBirksConstant(); 101 G4double destep = aStep->GetTotalEnergyD 138 G4double destep = aStep->GetTotalEnergyDeposit(); 102 G4double stepl = aStep->GetStepLength() 139 G4double stepl = aStep->GetStepLength(); 103 G4double charge = aStep->GetTrack()->Get 140 G4double charge = aStep->GetTrack()->GetDefinition()->GetPDGCharge(); 104 // 141 // 105 G4double response = destep; 142 G4double response = destep; 106 if (birk1*destep*stepl*charge != 0.) 143 if (birk1*destep*stepl*charge != 0.) 107 { 144 { 108 response = destep/(1. + birk1*destep/step 145 response = destep/(1. + birk1*destep/stepl); 109 } 146 } 110 return response; 147 return response; 111 } 148 } 112 149 113 //....oooOO0OOooo........oooOO0OOooo........oo 150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 114 151 115 152