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: FCALSteppingAction.cc,v 1.7 2006-06-29 16:03:15 gunter Exp $ >> 27 // GEANT4 tag $Name: not supported by cvs2svn $ 26 // 28 // 27 // 29 // 28 30 29 //....oooOO0OOooo........oooOO0OOooo........oo 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 30 //....oooOO0OOooo........oooOO0OOooo........oo 32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 31 33 32 #include <iostream> << 33 << 34 #include "FCALSteppingAction.hh" 34 #include "FCALSteppingAction.hh" 35 #include "G4SteppingManager.hh" 35 #include "G4SteppingManager.hh" 36 36 37 #include "globals.hh" << 38 #include "G4SystemOfUnits.hh" << 39 #include "G4Track.hh" 37 #include "G4Track.hh" 40 #include "G4DynamicParticle.hh" 38 #include "G4DynamicParticle.hh" 41 #include "G4Material.hh" 39 #include "G4Material.hh" 42 40 43 #include "G4LogicalVolume.hh" 41 #include "G4LogicalVolume.hh" 44 #include "G4VPhysicalVolume.hh" 42 #include "G4VPhysicalVolume.hh" 45 #include "G4VTouchable.hh" 43 #include "G4VTouchable.hh" 46 #include "G4TouchableHistory.hh" 44 #include "G4TouchableHistory.hh" 47 45 48 #include "G4Event.hh" 46 #include "G4Event.hh" 49 47 50 #include "G4ThreeVector.hh" 48 #include "G4ThreeVector.hh" 51 49 52 #include "G4ios.hh" 50 #include "G4ios.hh" >> 51 #include <iostream> >> 52 #include "globals.hh" 53 53 54 //....oooOO0OOooo........oooOO0OOooo........oo 54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 55 55 56 FCALSteppingAction::FCALSteppingAction():IDold 56 FCALSteppingAction::FCALSteppingAction():IDold(-1),IDout(-1) 57 {;} 57 {;} 58 58 59 //....oooOO0OOooo........oooOO0OOooo........oo 59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 60 60 61 FCALSteppingAction::~FCALSteppingAction() 61 FCALSteppingAction::~FCALSteppingAction() 62 {;} 62 {;} 63 63 64 //....oooOO0OOooo........oooOO0OOooo........oo 64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 65 65 66 void FCALSteppingAction::UserSteppingAction(co 66 void FCALSteppingAction::UserSteppingAction(const G4Step* astep) 67 { 67 { 68 // Get Edep 68 // Get Edep 69 G4double Edep = astep->GetTotalEnergyDeposi 69 G4double Edep = astep->GetTotalEnergyDeposit(); 70 70 71 // Get Track 71 // Get Track 72 G4Track* aTrack = astep->GetTrack(); 72 G4Track* aTrack = astep->GetTrack(); 73 73 74 // Get Touchable History 74 // Get Touchable History 75 G4TouchableHistory* theTouchable = (G4Touch 75 G4TouchableHistory* theTouchable = (G4TouchableHistory*)(aTrack->GetTouchable()); 76 76 77 // Energy deposit in FCAL1 and FCAL2 77 // Energy deposit in FCAL1 and FCAL2 78 if(Edep != 0.) 78 if(Edep != 0.) 79 { 79 { 80 G4VPhysicalVolume* physVol = theTouchabl 80 G4VPhysicalVolume* physVol = theTouchable->GetVolume(); 81 81 82 if(strcmp(physVol->GetName(),"FCALEmModu 82 if(strcmp(physVol->GetName(),"FCALEmModulePhysical")== 0 || 83 strcmp(physVol->GetName(),"F1LArGapPhysical 83 strcmp(physVol->GetName(),"F1LArGapPhysical") == 0) 84 { 84 { 85 EdepFCALEm = EdepFCALEm + Edep; 85 EdepFCALEm = EdepFCALEm + Edep; 86 }; 86 }; 87 87 88 if( (strcmp(physVol->GetName(), "FCALHad 88 if( (strcmp(physVol->GetName(), "FCALHadModulePhysical") == 0) || 89 (strcmp(physVol->GetName(), "CuPlateAPhysi 89 (strcmp(physVol->GetName(), "CuPlateAPhysical") == 0) || 90 (strcmp(physVol->GetName(), "CuPlateBPhysi 90 (strcmp(physVol->GetName(), "CuPlateBPhysical") == 0) || 91 (strcmp(physVol->GetName(), "WAbsorberPhys 91 (strcmp(physVol->GetName(), "WAbsorberPhysical") == 0) || 92 (strcmp(physVol->GetName(), "F2RodPhysical 92 (strcmp(physVol->GetName(), "F2RodPhysical") == 0) || 93 (strcmp(physVol->GetName(), "F2LArGapPhysi 93 (strcmp(physVol->GetName(), "F2LArGapPhysical") == 0) ) 94 { 94 { 95 EdepFCALHad = EdepFCALHad + Edep; 95 EdepFCALHad = EdepFCALHad + Edep; 96 }; 96 }; 97 }; 97 }; 98 98 99 // Get Tracks properties 99 // Get Tracks properties 100 G4int TrackID = aTrack->GetTrackID(); 100 G4int TrackID = aTrack->GetTrackID(); 101 G4int ParentID = aTrack->GetParentID(); 101 G4int ParentID = aTrack->GetParentID(); 102 // Get Associated particle 102 // Get Associated particle 103 const G4DynamicParticle * aDynamicParticle = 103 const G4DynamicParticle * aDynamicParticle = aTrack->GetDynamicParticle(); 104 G4ParticleDefinition * aParticle = aTrack->G 104 G4ParticleDefinition * aParticle = aTrack->GetDefinition(); 105 G4String ParticleName = aParticle->GetPartic 105 G4String ParticleName = aParticle->GetParticleName(); 106 106 107 IDnow = EventNo + 10000*TrackID+ 100000000*P 107 IDnow = EventNo + 10000*TrackID+ 100000000*ParentID; 108 108 109 if(IDnow != IDold) 109 if(IDnow != IDold) 110 { 110 { 111 IDold = IDnow; 111 IDold = IDnow; 112 112 113 // Get the primary particle 113 // Get the primary particle 114 if(TrackID==1 && ParentID==0 && (aTrack- 114 if(TrackID==1 && ParentID==0 && (aTrack->GetCurrentStepNumber()) == 1) 115 { 115 { 116 PrimaryVertex = aTrack->GetVertexPositi 116 PrimaryVertex = aTrack->GetVertexPosition(); 117 PrimaryDirection = aTrack->GetVertexMoment 117 PrimaryDirection = aTrack->GetVertexMomentumDirection(); 118 118 119 NSecondaries = 1; 119 NSecondaries = 1; 120 Secondaries[NSecondaries][1] = aParticle-> 120 Secondaries[NSecondaries][1] = aParticle->GetPDGEncoding(); 121 Secondaries[NSecondaries][2] = PrimaryVert 121 Secondaries[NSecondaries][2] = PrimaryVertex.x(); 122 Secondaries[NSecondaries][3] = PrimaryVert 122 Secondaries[NSecondaries][3] = PrimaryVertex.y(); 123 Secondaries[NSecondaries][4] = PrimaryVert 123 Secondaries[NSecondaries][4] = PrimaryVertex.z(); 124 Secondaries[NSecondaries][5] = (aDynamicPa 124 Secondaries[NSecondaries][5] = (aDynamicParticle->GetMomentum()).x(); 125 Secondaries[NSecondaries][6] = (aDynamicPa 125 Secondaries[NSecondaries][6] = (aDynamicParticle->GetMomentum()).y(); 126 Secondaries[NSecondaries][7] = (aDynamicPa 126 Secondaries[NSecondaries][7] = (aDynamicParticle->GetMomentum()).z(); 127 Secondaries[NSecondaries][8] = aDynamicPar 127 Secondaries[NSecondaries][8] = aDynamicParticle->GetTotalMomentum(); 128 Secondaries[NSecondaries][9] = aDynamicPar 128 Secondaries[NSecondaries][9] = aDynamicParticle->GetTotalEnergy(); 129 Secondaries[NSecondaries][10] = aDynamicPa 129 Secondaries[NSecondaries][10] = aDynamicParticle->GetKineticEnergy(); 130 130 131 G4cout << " **** Primary : " << EventNo < 131 G4cout << " **** Primary : " << EventNo << G4endl; 132 G4cout << " Vertex : " << PrimaryVertex << 132 G4cout << " Vertex : " << PrimaryVertex << G4endl; 133 } 133 } 134 134 135 135 136 // Get secondaries in air close to the p 136 // Get secondaries in air close to the primary tracks (DCA < 2.mm) 137 G4double DCACut = 2.*mm; 137 G4double DCACut = 2.*mm; 138 G4String Material = aTrack->GetMaterial( 138 G4String Material = aTrack->GetMaterial()->GetName(); 139 G4ThreeVector TrackPos = aTrack->GetVert 139 G4ThreeVector TrackPos = aTrack->GetVertexPosition(); 140 140 141 if(TrackID != 1 && ParentID == 1 && (str 141 if(TrackID != 1 && ParentID == 1 && (strcmp(Material,"Air")==0) && (TrackPos.z() > 135.*cm)) 142 { 142 { 143 SecondaryVertex = aTrack->GetVertexPositio 143 SecondaryVertex = aTrack->GetVertexPosition(); 144 SecondaryDirection = aTrack->GetVertexMome 144 SecondaryDirection = aTrack->GetVertexMomentumDirection(); 145 145 146 // calculate DCA of secondries to primary 146 // calculate DCA of secondries to primary particle 147 Distance = PrimaryVertex - SecondaryVertex 147 Distance = PrimaryVertex - SecondaryVertex ; 148 VectorProduct = PrimaryDirection.cross(Sec 148 VectorProduct = PrimaryDirection.cross(SecondaryDirection); 149 if(VectorProduct == G4ThreeVector() && 149 if(VectorProduct == G4ThreeVector() && 150 PrimaryDirection != G4ThreeVector() && 150 PrimaryDirection != G4ThreeVector() && SecondaryDirection != G4ThreeVector()) 151 { 151 { 152 G4ThreeVector Temp = Distance.cross(Pr 152 G4ThreeVector Temp = Distance.cross(PrimaryDirection); 153 VectorProduct = Temp.cross(PrimaryDire 153 VectorProduct = Temp.cross(PrimaryDirection); 154 }; 154 }; 155 155 156 VectorProductMagnitude = VectorProduct.mag(); 156 VectorProductMagnitude = VectorProduct.mag(); 157 if(VectorProductMagnitude == 0.) 157 if(VectorProductMagnitude == 0.) 158 { 158 { 159 VectorProductNorm = G4ThreeVector(); 159 VectorProductNorm = G4ThreeVector(); 160 } else { 160 } else { 161 VectorProductNorm = (1./VectorProduct. 161 VectorProductNorm = (1./VectorProduct.mag()) * VectorProduct ; 162 }; 162 }; 163 DistOfClosestApproach = Distance * VectorP 163 DistOfClosestApproach = Distance * VectorProductNorm ; 164 164 165 if(std::abs(DistOfClosestApproach) < DCACu 165 if(std::abs(DistOfClosestApproach) < DCACut) 166 { 166 { 167 NSecondaries++; 167 NSecondaries++; 168 Secondaries[0][0] = NSecondaries; 168 Secondaries[0][0] = NSecondaries; 169 Secondaries[NSecondaries][1] = aPartic 169 Secondaries[NSecondaries][1] = aParticle->GetPDGEncoding(); 170 Secondaries[NSecondaries][2] = (aTrack 170 Secondaries[NSecondaries][2] = (aTrack->GetVertexPosition()).x(); 171 Secondaries[NSecondaries][3] = (aTrack 171 Secondaries[NSecondaries][3] = (aTrack->GetVertexPosition()).y(); 172 Secondaries[NSecondaries][4] = (aTrack 172 Secondaries[NSecondaries][4] = (aTrack->GetVertexPosition()).z(); 173 Secondaries[NSecondaries][5] =(aDynami 173 Secondaries[NSecondaries][5] =(aDynamicParticle->GetMomentum()).x(); 174 Secondaries[NSecondaries][6] = (aDynam 174 Secondaries[NSecondaries][6] = (aDynamicParticle->GetMomentum()).y(); 175 Secondaries[NSecondaries][7] = (aDynam 175 Secondaries[NSecondaries][7] = (aDynamicParticle->GetMomentum()).z(); 176 Secondaries[NSecondaries][8] = aDynami 176 Secondaries[NSecondaries][8] = aDynamicParticle->GetTotalMomentum(); 177 Secondaries[NSecondaries][9] = aDynami 177 Secondaries[NSecondaries][9] = aDynamicParticle->GetTotalEnergy(); 178 Secondaries[NSecondaries][10] =aDynami 178 Secondaries[NSecondaries][10] =aDynamicParticle->GetKineticEnergy(); 179 }; 179 }; 180 }; 180 }; 181 }; 181 }; 182 182 183 183 184 // Get the World leaving particle 184 // Get the World leaving particle 185 if(aTrack->GetNextVolume() == 0) { 185 if(aTrack->GetNextVolume() == 0) { 186 if(IDnow != IDout) { 186 if(IDnow != IDout) { 187 IDout = IDnow; 187 IDout = IDnow; 188 188 189 NTracks++; 189 NTracks++; 190 190 191 OutOfWorldTracksData[0][0] = NTracks; 191 OutOfWorldTracksData[0][0] = NTracks; 192 192 193 OutOfWorldTracksData[NTracks][1] = aPart 193 OutOfWorldTracksData[NTracks][1] = aParticle->GetPDGEncoding(); 194 194 195 OutOfWorldTracksData[NTracks][2] = (aTra 195 OutOfWorldTracksData[NTracks][2] = (aTrack->GetVertexPosition()).x(); 196 OutOfWorldTracksData[NTracks][3] = (aTra 196 OutOfWorldTracksData[NTracks][3] = (aTrack->GetVertexPosition()).y(); 197 OutOfWorldTracksData[NTracks][4] = (aTra 197 OutOfWorldTracksData[NTracks][4] = (aTrack->GetVertexPosition()).z(); 198 198 199 OutOfWorldTracksData[NTracks][5] = (aDyn 199 OutOfWorldTracksData[NTracks][5] = (aDynamicParticle->GetMomentum()).x(); 200 OutOfWorldTracksData[NTracks][6] = (aDyn 200 OutOfWorldTracksData[NTracks][6] = (aDynamicParticle->GetMomentum()).y(); 201 OutOfWorldTracksData[NTracks][7] = (aDyn 201 OutOfWorldTracksData[NTracks][7] = (aDynamicParticle->GetMomentum()).z(); 202 202 203 OutOfWorldTracksData[NTracks][8] = aDyna 203 OutOfWorldTracksData[NTracks][8] = aDynamicParticle->GetTotalMomentum(); 204 204 205 OutOfWorldTracksData[NTracks][9] = aDyna 205 OutOfWorldTracksData[NTracks][9] = aDynamicParticle->GetTotalEnergy(); 206 206 207 OutOfWorldTracksData[NTracks][10] = aDyn 207 OutOfWorldTracksData[NTracks][10] = aDynamicParticle->GetKineticEnergy(); 208 }; 208 }; 209 }; 209 }; 210 210 211 211 212 } 212 } 213 213 214 void FCALSteppingAction::initialize(G4int Nev) 214 void FCALSteppingAction::initialize(G4int Nev) { 215 EventNo = Nev; 215 EventNo = Nev; 216 NTracks = 0; 216 NTracks = 0; 217 NSecondaries = 0; 217 NSecondaries = 0; 218 EdepFCALEm = EdepFCALHad = 0.; 218 EdepFCALEm = EdepFCALHad = 0.; 219 219 220 for(G4int i=0; i<6000; i++) 220 for(G4int i=0; i<6000; i++) 221 { 221 { 222 for(G4int j=0; j<11; j++) 222 for(G4int j=0; j<11; j++) 223 { 223 { 224 OutOfWorldTracksData[i][j] = 0.; 224 OutOfWorldTracksData[i][j] = 0.; 225 Secondaries[i][j] = 0.; 225 Secondaries[i][j] = 0.; 226 } 226 } 227 }; 227 }; 228 } 228 } 229 229 230 G4double FCALSteppingAction::GetOutOfWorldTrac 230 G4double FCALSteppingAction::GetOutOfWorldTracks(G4int i, G4int j){ 231 return OutOfWorldTracksData[i][j]; 231 return OutOfWorldTracksData[i][j]; 232 } 232 } 233 233 234 G4double FCALSteppingAction::GetSecondaries(G4 234 G4double FCALSteppingAction::GetSecondaries(G4int i, G4int j){ 235 return Secondaries[i][j]; 235 return Secondaries[i][j]; 236 } 236 } 237 237 238 G4double FCALSteppingAction::GetEdepFCAL(G4Str 238 G4double FCALSteppingAction::GetEdepFCAL(G4String FCAL) { 239 if(strcmp(FCAL,"FCALEm") == 0) { 239 if(strcmp(FCAL,"FCALEm") == 0) { 240 return EdepFCALEm; 240 return EdepFCALEm; 241 } else { 241 } else { 242 if(strcmp(FCAL,"FCALHad") == 0) { 242 if(strcmp(FCAL,"FCALHad") == 0) { 243 return EdepFCALHad;} 243 return EdepFCALHad;} 244 } 244 } 245 return 0.0; 245 return 0.0; 246 } 246 } 247 247 248 248 249 //....oooOO0OOooo........oooOO0OOooo........oo 249 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 250 250 251 251 252 252 253 253