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 ////////////////////////////////////////////// 26 /////////////////////////////////////////////////////////////////////////////// 27 // File: CCaloSD.cc 27 // File: CCaloSD.cc 28 // Description: Stores hits of calorimetric ty 28 // Description: Stores hits of calorimetric type in appropriate container 29 ////////////////////////////////////////////// 29 /////////////////////////////////////////////////////////////////////////////// 30 30 31 #include "CCaloSD.hh" 31 #include "CCaloSD.hh" 32 #include "G4SystemOfUnits.hh" 32 #include "G4SystemOfUnits.hh" 33 #include "G4VProcess.hh" 33 #include "G4VProcess.hh" 34 #include "G4SDManager.hh" 34 #include "G4SDManager.hh" 35 #include "G4VTouchable.hh" 35 #include "G4VTouchable.hh" 36 #include "CCalVOrganization.hh" 36 #include "CCalVOrganization.hh" 37 #include "CCalSDList.hh" 37 #include "CCalSDList.hh" 38 38 39 //#define debug 39 //#define debug 40 //#define ddebug 40 //#define ddebug 41 41 42 CCaloSD::CCaloSD(G4String name, CCalVOrganizat 42 CCaloSD::CCaloSD(G4String name, CCalVOrganization* numberingScheme): 43 G4VSensitiveDetector(name), HCID(-1), SDname 43 G4VSensitiveDetector(name), HCID(-1), SDname(name), theHC(0), 44 CurrentHit(0), theTrack(0), CurrentPV(0), Pr 44 CurrentHit(0), theTrack(0), CurrentPV(0), PreviousPV(0), UnitID(0), 45 PreviousUnitID(0), PreStepPoint(0), PostStep 45 PreviousUnitID(0), PreStepPoint(0), PostStepPoint(0), 46 theDescription(numberingScheme) { 46 theDescription(numberingScheme) { 47 47 48 collectionName.insert(name); 48 collectionName.insert(name); 49 49 50 G4cout << "********************************* 50 G4cout << "*******************************************************" << G4endl; 51 G4cout << "* 51 G4cout << "* *" << G4endl; 52 G4cout << "* Constructing a CCaloSD with na << 52 G4cout << "* Constructing a CCaloSD with name " << name << G4endl; 53 G4cout << "* 53 G4cout << "* *" << G4endl; 54 G4cout << "********************************* 54 G4cout << "*******************************************************" << G4endl; 55 55 56 CCalSDList::getInstance()->addCalo(name); 56 CCalSDList::getInstance()->addCalo(name); 57 } 57 } 58 58 >> 59 59 CCaloSD::~CCaloSD() { 60 CCaloSD::~CCaloSD() { 60 delete theDescription; << 61 G4cout << "CCaloSD: " << SDname << " deleted" << G4endl; >> 62 if (theDescription) >> 63 delete theDescription; 61 } 64 } 62 65 >> 66 63 void CCaloSD::Initialize(G4HCofThisEvent*HCE) 67 void CCaloSD::Initialize(G4HCofThisEvent*HCE) { 64 68 65 #ifdef debug 69 #ifdef debug 66 G4cout << "CCaloSD : Initialize called for " 70 G4cout << "CCaloSD : Initialize called for " << SDname << G4endl; 67 #endif 71 #endif 68 //This initialization is performed at the be 72 //This initialization is performed at the beginning of an event 69 //------------------------------------------ 73 //------------------------------------------------------------ 70 74 71 theHC = new CCalG4HitCollection(SDname, coll 75 theHC = new CCalG4HitCollection(SDname, collectionName[0]); 72 if (HCID<0) { << 76 if (HCID<0) 73 HCID = G4SDManager::GetSDMpointer()->GetCo 77 HCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]); 74 } << 75 HCE->AddHitsCollection( HCID, theHC ); 78 HCE->AddHitsCollection( HCID, theHC ); 76 79 77 TSID = -2; 80 TSID = -2; 78 PrimID = -2; 81 PrimID = -2; 79 /////PrimaryID = -99; <--- initialized by S 82 /////PrimaryID = -99; <--- initialized by StackingAction. 80 } 83 } 81 84 82 G4bool CCaloSD::ProcessHits(G4Step*aStep, G4To << 83 85 84 if(aStep->GetTotalEnergyDeposit() > CLHEP::e << 86 G4bool CCaloSD::ProcessHits(G4Step*aStep,G4TouchableHistory*) { 85 getStepInfo(aStep); << 87 86 if (hitExists() == false && EdepositEM+Ede << 88 if (aStep == NULL) return true; 87 createNewHit(); << 89 88 } << 90 getStepInfo(aStep); 89 } << 91 if (hitExists() == false && EdepositEM+EdepositEHAD>0.) >> 92 createNewHit(); >> 93 90 return true; 94 return true; 91 } 95 } 92 96 93 void CCaloSD::getStepInfo(const G4Step* aStep) << 97 >> 98 void CCaloSD::getStepInfo(G4Step* aStep) { 94 99 95 PreStepPoint = aStep->GetPreStepPoint(); 100 PreStepPoint = aStep->GetPreStepPoint(); 96 PostStepPoint= aStep->GetPostStepPoint(); 101 PostStepPoint= aStep->GetPostStepPoint(); 97 HitPoint = PreStepPoint->GetPosition(); << 102 HitPoint = PreStepPoint->GetPosition(); 98 103 99 theTrack = aStep->GetTrack(); 104 theTrack = aStep->GetTrack(); 100 CurrentPV= PreStepPoint->GetPhysicalVolume() 105 CurrentPV= PreStepPoint->GetPhysicalVolume(); 101 106 102 G4String pname = theTrack->GetDefinition()-> << 107 G4String particleType = theTrack->GetDefinition()->GetParticleName(); 103 G4double de = aStep->GetTotalEnergyDeposit() << 104 //G4cout << "##### Process new step dE= " << << 105 // << " " << pname << " inside " << << 106 108 107 const G4VTouchable* touch = aStep->GetPreSte 109 const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable(); 108 G4double weight = 1.; << 110 G4double weight; 109 if (touch->GetVolume(0)->GetName() == "Cryst << 111 if (touch->GetVolume(0)->GetName() == "CrystalMatrixCrystal") 110 weight = curve_LY(PreStepPoint); 112 weight = curve_LY(PreStepPoint); 111 } << 113 else >> 114 weight = 1.; 112 115 113 if (pname == "e-" || pname == "e+" || pname << 116 if (particleType == "e-" || 114 EdepositEM = weight*de; << 117 particleType == "e+" || 115 EdepositEHAD = 0.f; << 118 particleType == "gamma" ){ >> 119 EdepositEM = weight*(aStep->GetTotalEnergyDeposit()); >> 120 EdepositEHAD = 0.; 116 } else { 121 } else { 117 EdepositEM = 0.f; << 122 EdepositEM = 0.; 118 EdepositEHAD = weight*de; << 123 EdepositEHAD = weight*(aStep->GetTotalEnergyDeposit()); 119 } 124 } 120 125 121 TSlice = (PostStepPoint) ? PostStepPoint->Ge << 126 TSlice = PostStepPoint->GetGlobalTime() / nanosecond; 122 //G4cout << " W= " << weight << " T= " < << 127 123 G4int it = (G4int)TSlice; << 128 if ( TSlice > 1.0E9 ) TSliceID = 999999999; 124 TSliceID = std::min(it, 10000); << 129 else TSliceID = (int) TSlice; 125 //G4cout << " tID= " << TSliceID << G4endl; << 126 130 127 UnitID = (theDescription) ? theDescription-> << 131 if (theDescription!=0) >> 132 UnitID = theDescription->GetUnitID(aStep); >> 133 else >> 134 UnitID = 0; >> 135 128 } 136 } 129 137 >> 138 130 G4bool CCaloSD::hitExists() { 139 G4bool CCaloSD::hitExists() { 131 140 >> 141 if (PrimaryID<1) { >> 142 G4cerr << "***** CCaloSD error: PrimaryID = " << PrimaryID >> 143 << " Maybe detector name changed" >> 144 << G4endl; >> 145 } >> 146 >> 147 132 if ( CurrentPV==PreviousPV && PrimaryID == P 148 if ( CurrentPV==PreviousPV && PrimaryID == PrimID && TSliceID == TSID && 133 UnitID==PreviousUnitID) { 149 UnitID==PreviousUnitID) { 134 updateHit(); 150 updateHit(); 135 return true; 151 return true; 136 } 152 } 137 153 138 if (PrimaryID != PrimID) { ResetForNewPrimar << 154 >> 155 if (PrimaryID != PrimID) >> 156 ResetForNewPrimary(); 139 157 140 // look in HC whether a hit with the same pr << 158 141 // exists already: << 159 //look in HC whether a hit with the same primID,UnitID,TSliceID exists already: 142 160 143 G4bool found = false; 161 G4bool found = false; 144 for (std::size_t j=0; j<theHC->entries(); ++ << 162 for (G4int j=0; j<theHC->entries()&&!found; j++) { 145 163 146 CCalG4Hit* aPreviousHit = (*theHC)[j]; 164 CCalG4Hit* aPreviousHit = (*theHC)[j]; 147 if (aPreviousHit->getTrackID() == Primary 165 if (aPreviousHit->getTrackID() == PrimaryID && 148 aPreviousHit->getTimeSliceID() == TSli << 166 aPreviousHit->getTimeSliceID() == TSliceID && 149 aPreviousHit->getUnitID()== UnitID << 167 aPreviousHit->getUnitID()== UnitID ) { 150 CurrentHit = aPreviousHit; 168 CurrentHit = aPreviousHit; 151 found = true; 169 found = true; 152 break; << 153 } 170 } 154 } 171 } 155 172 156 if (found) { 173 if (found) { 157 updateHit(); 174 updateHit(); 158 return true; 175 return true; 159 } else { 176 } else { 160 return false; 177 return false; 161 } 178 } 162 } 179 } 163 180 >> 181 164 void CCaloSD::ResetForNewPrimary() { 182 void CCaloSD::ResetForNewPrimary() { 165 183 166 EntrancePoint = SetToLocal(HitPoint); 184 EntrancePoint = SetToLocal(HitPoint); 167 IncidentEnergy = PreStepPoint->GetKineticEne 185 IncidentEnergy = PreStepPoint->GetKineticEnergy(); 168 186 169 } 187 } 170 188 >> 189 171 void CCaloSD::StoreHit(CCalG4Hit* hit){ 190 void CCaloSD::StoreHit(CCalG4Hit* hit){ 172 191 173 if (PrimID<0) return; 192 if (PrimID<0) return; 174 if (hit == 0 ) { 193 if (hit == 0 ) { 175 G4cout << "CCaloSD: hit to be stored is NU 194 G4cout << "CCaloSD: hit to be stored is NULL !!" <<G4endl; 176 return; 195 return; 177 } 196 } 178 197 179 theHC->insert( hit ); 198 theHC->insert( hit ); 180 } 199 } 181 200 >> 201 182 void CCaloSD::createNewHit() { 202 void CCaloSD::createNewHit() { 183 203 184 #ifdef debug 204 #ifdef debug 185 G4int currentCopyNo = -999; 205 G4int currentCopyNo = -999; 186 G4int motherCopyNo = -999; 206 G4int motherCopyNo = -999; 187 G4TouchableHistory* theTouchable 207 G4TouchableHistory* theTouchable 188 = (G4TouchableHistory*)( theTrack->GetTouc 208 = (G4TouchableHistory*)( theTrack->GetTouchable() ); 189 if ( theTouchable ) { 209 if ( theTouchable ) { 190 currentCopyNo = theTouchable->GetReplicaNu 210 currentCopyNo = theTouchable->GetReplicaNumber( 0 ); 191 if ( theTouchable->GetHistoryDepth() > 0 ) 211 if ( theTouchable->GetHistoryDepth() > 0 ) { 192 motherCopyNo = theTouchable->GetReplicaN 212 motherCopyNo = theTouchable->GetReplicaNumber( 1 ); 193 } 213 } 194 } 214 } 195 G4cout << "CCaloSD createNewHit for" 215 G4cout << "CCaloSD createNewHit for" 196 << " PV " << CurrentPV->GetName() 216 << " PV " << CurrentPV->GetName() 197 << " PVid = " << currentCopyNo 217 << " PVid = " << currentCopyNo 198 << " MVid = " << motherCopyNo 218 << " MVid = " << motherCopyNo 199 << " Unit " << UnitID <<G4endl; 219 << " Unit " << UnitID <<G4endl; 200 G4cout << " primary " << PrimaryID 220 G4cout << " primary " << PrimaryID 201 << " time slice " << TSliceID 221 << " time slice " << TSliceID 202 << " For Track " << theTrack->GetTra 222 << " For Track " << theTrack->GetTrackID() 203 << " which is a " << theTrack->GetDe 223 << " which is a " << theTrack->GetDefinition()->GetParticleName(); 204 << 224 205 if (theTrack->GetTrackID()==1) { 225 if (theTrack->GetTrackID()==1) { 206 G4cout << " of energy " << theTrack->G 226 G4cout << " of energy " << theTrack->GetTotalEnergy(); 207 } else { 227 } else { 208 G4cout << " daughter of part. " << theTrac 228 G4cout << " daughter of part. " << theTrack->GetParentID(); 209 } 229 } 210 230 211 G4cout << " and created by " ; 231 G4cout << " and created by " ; 212 if (theTrack->GetCreatorProcess()!=NULL) 232 if (theTrack->GetCreatorProcess()!=NULL) 213 G4cout << theTrack->GetCreatorProcess()->G 233 G4cout << theTrack->GetCreatorProcess()->GetProcessName() ; 214 else 234 else 215 G4cout << "NO process"; 235 G4cout << "NO process"; 216 G4cout << G4endl; 236 G4cout << G4endl; 217 #endif 237 #endif 218 238 >> 239 219 CurrentHit = new CCalG4Hit; 240 CurrentHit = new CCalG4Hit; 220 CurrentHit->setTrackID(PrimaryID); 241 CurrentHit->setTrackID(PrimaryID); 221 CurrentHit->setTimeSlice(TSlice); 242 CurrentHit->setTimeSlice(TSlice); 222 CurrentHit->setUnitID(UnitID); 243 CurrentHit->setUnitID(UnitID); 223 CurrentHit->setEntry(EntrancePoint); 244 CurrentHit->setEntry(EntrancePoint); 224 CurrentHit->setIncidentEnergy(IncidentEnergy 245 CurrentHit->setIncidentEnergy(IncidentEnergy); 225 updateHit(); 246 updateHit(); 226 247 227 StoreHit(CurrentHit); 248 StoreHit(CurrentHit); 228 249 229 } << 250 } >> 251 230 252 231 void CCaloSD::updateHit() { 253 void CCaloSD::updateHit() { 232 if (EdepositEM+EdepositEHAD != 0) { 254 if (EdepositEM+EdepositEHAD != 0) { 233 CurrentHit->addEnergyDeposit(EdepositEM,Ed 255 CurrentHit->addEnergyDeposit(EdepositEM,EdepositEHAD); 234 #ifdef debug 256 #ifdef debug 235 G4cout << "Energy deposit in Unit " << Uni 257 G4cout << "Energy deposit in Unit " << UnitID << " em " << EdepositEM/MeV 236 << " hadronic " << EdepositEHAD/MeV < << 258 << " hadronic " << EdepositEHAD/MeV << " MeV" << G4endl; 237 #endif 259 #endif 238 } 260 } 239 261 240 // buffer for next steps: 262 // buffer for next steps: 241 TSID = TSliceID; 263 TSID = TSliceID; 242 PrimID = PrimaryID; 264 PrimID = PrimaryID; 243 PreviousPV = CurrentPV; 265 PreviousPV = CurrentPV; 244 PreviousUnitID = UnitID; 266 PreviousUnitID = UnitID; 245 } 267 } 246 268 247 G4ThreeVector CCaloSD::SetToLocal(const G4Thre << 269 >> 270 G4ThreeVector CCaloSD::SetToLocal(G4ThreeVector global){ 248 271 249 G4ThreeVector localPoint; 272 G4ThreeVector localPoint; 250 const G4VTouchable* touch= PreStepPoint->G 273 const G4VTouchable* touch= PreStepPoint->GetTouchable(); 251 localPoint=touch->GetHistory()->GetTopTransf 274 localPoint=touch->GetHistory()->GetTopTransform().TransformPoint(global); 252 275 253 return localPoint; 276 return localPoint; >> 277 254 } 278 } 255 279 >> 280 256 void CCaloSD::EndOfEvent(G4HCofThisEvent*) { 281 void CCaloSD::EndOfEvent(G4HCofThisEvent*) { 257 summarize(); 282 summarize(); 258 } 283 } 259 284 >> 285 260 void CCaloSD::summarize() { 286 void CCaloSD::summarize() { 261 } 287 } 262 288 >> 289 263 void CCaloSD::clear() { 290 void CCaloSD::clear() { 264 } 291 } 265 292 >> 293 266 void CCaloSD::DrawAll() { 294 void CCaloSD::DrawAll() { 267 } 295 } 268 296 >> 297 269 void CCaloSD::PrintAll() { 298 void CCaloSD::PrintAll() { 270 G4cout << "CCaloSD: Collection " << theHC->G 299 G4cout << "CCaloSD: Collection " << theHC->GetName() << G4endl; 271 theHC->PrintAllHits(); 300 theHC->PrintAllHits(); 272 } 301 } 273 302 >> 303 274 void CCaloSD::SetOrganization(CCalVOrganizatio 304 void CCaloSD::SetOrganization(CCalVOrganization* org){ 275 305 276 if (theDescription!=0) 306 if (theDescription!=0) 277 delete theDescription; 307 delete theDescription; 278 theDescription = org; 308 theDescription = org; 279 } 309 } 280 310 281 G4double CCaloSD::curve_LY(const G4StepPoint* << 311 >> 312 G4double CCaloSD::curve_LY(G4StepPoint* stepPoint) { 282 313 283 G4double weight = 1.; 314 G4double weight = 1.; 284 G4ThreeVector localPoint = SetToLocal(stepPo << 315 G4ThreeVector localPoint = SetToLocal(stepPoint->GetPosition()); 285 const G4double crlength = 230.; << 316 G4double crlength = 230.; 286 G4double dapd = 0.5 * crlength - localPoint. 317 G4double dapd = 0.5 * crlength - localPoint.z(); 287 if (dapd >= -0.1 || dapd <= crlength+0.1) { 318 if (dapd >= -0.1 || dapd <= crlength+0.1) { 288 if (dapd <= 100.) 319 if (dapd <= 100.) 289 weight = 1.05 - dapd * 0.0005; 320 weight = 1.05 - dapd * 0.0005; 290 } else { 321 } else { 291 G4cout << "CCaloSD, light coll curve : wro 322 G4cout << "CCaloSD, light coll curve : wrong distance to APD " << dapd 292 << " crlength = " << crlength << 323 << " crlength = " << crlength 293 << " z of localPoint = " << localPo << 324 << " z of localPoint = " << localPoint.z() 294 << " take weight = " << weight << G << 325 << " take weight = " << weight << G4endl; 295 } 326 } 296 #ifdef ddebug 327 #ifdef ddebug 297 G4cout << "CCaloSD, light coll curve : " << 328 G4cout << "CCaloSD, light coll curve : " << dapd 298 << " crlength = " << crlength << 329 << " crlength = " << crlength 299 << " z of localPoint = " << localPoin << 330 << " z of localPoint = " << localPoint.z() 300 << " take weight = " << weight << G4e << 331 << " take weight = " << weight << G4endl; 301 #endif 332 #endif 302 return weight; 333 return weight; 303 } 334 } 304 335