// // ******************************************************************** // * License and Disclaimer * // * * // * The Geant4 software is copyright of the Copyright Holders of * // * the Geant4 Collaboration. It is provided under the terms and * // * conditions of the Geant4 Software License, included in the file * // * LICENSE and available at http://cern.ch/geant4/license . These * // * include a list of copyright holders. * // * * // * Neither the authors of this software system, nor their employing * // * institutes,nor the agencies providing financial support for this * // * work make any representation or warranty, express or implied, * // * regarding this software system or assume any liability for its * // * use. Please see the license in the file LICENSE and URL above * // * for the full disclaimer and the limitation of liability. * // * * // * This code implementation is the result of the scientific and * // * technical work of the GEANT4 collaboration. * // * By using, copying, modifying or distributing the software (or * // * any work based on the software) you agree to acknowledge its * // * use in resulting scientific publications, and indicate your * // * acceptance of all terms of the Geant4 Software license. * // ******************************************************************** // /////////////////////////////////////////////////////////////////////////////// // File: CCaloSD.cc // Description: Stores hits of calorimetric type in appropriate container /////////////////////////////////////////////////////////////////////////////// #include "CCaloSD.hh" #include "G4VProcess.hh" #include "G4SDManager.hh" #include "G4VTouchable.hh" #include "CCalVOrganization.hh" #include "CCalSDList.hh" //#define debug //#define ddebug CCaloSD::CCaloSD(G4String name, CCalVOrganization* numberingScheme): G4VSensitiveDetector(name), HCID(-1), SDname(name), theHC(0), CurrentHit(0), theTrack(0), CurrentPV(0), PreviousPV(0), UnitID(0), PreviousUnitID(0), PreStepPoint(0), PostStepPoint(0), theDescription(numberingScheme) { collectionName.insert(name); G4cout << "*******************************************************" << G4endl; G4cout << "* *" << G4endl; G4cout << "* Constructing a CCaloSD with name " << name << G4endl; G4cout << "* *" << G4endl; G4cout << "*******************************************************" << G4endl; CCalSDList::getInstance()->addCalo(name); } CCaloSD::~CCaloSD() { G4cout << "CCaloSD: " << SDname << " deleted" << G4endl; if (theDescription) delete theDescription; } void CCaloSD::Initialize(G4HCofThisEvent*HCE) { #ifdef debug G4cout << "CCaloSD : Initialize called for " << SDname << G4endl; #endif //This initialization is performed at the beginning of an event //------------------------------------------------------------ theHC = new CCalG4HitCollection(SDname, collectionName[0]); if (HCID<0) HCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]); HCE->AddHitsCollection( HCID, theHC ); TSID = -2; PrimID = -2; /////PrimaryID = -99; <--- initialized by StackingAction. } G4bool CCaloSD::ProcessHits(G4Step*aStep,G4TouchableHistory*) { if (aStep == NULL) return true; getStepInfo(aStep); if (hitExists() == false && EdepositEM+EdepositEHAD>0.) createNewHit(); return true; } void CCaloSD::getStepInfo(G4Step* aStep) { PreStepPoint = aStep->GetPreStepPoint(); PostStepPoint= aStep->GetPostStepPoint(); HitPoint = PreStepPoint->GetPosition(); theTrack = aStep->GetTrack(); CurrentPV= PreStepPoint->GetPhysicalVolume(); G4String particleType = theTrack->GetDefinition()->GetParticleName(); const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable(); G4double weight; if (touch->GetVolume(0)->GetName() == "CrystalMatrixCrystal") weight = curve_LY(PreStepPoint); else weight = 1.; if (particleType == "e-" || particleType == "e+" || particleType == "gamma" ){ EdepositEM = weight*(aStep->GetTotalEnergyDeposit()); EdepositEHAD = 0.; } else { EdepositEM = 0.; EdepositEHAD = weight*(aStep->GetTotalEnergyDeposit()); } TSlice = (PostStepPoint->GetGlobalTime() )/nanosecond; TSliceID = (int) TSlice; if (theDescription!=0) UnitID = theDescription->GetUnitID(aStep); else UnitID = 0; } G4bool CCaloSD::hitExists() { if (PrimaryID<1) { G4cerr << "***** CCaloSD error: PrimaryID = " << PrimaryID << " Maybe detector name changed" << G4endl; } if ( CurrentPV==PreviousPV && PrimaryID == PrimID && TSliceID == TSID && UnitID==PreviousUnitID) { updateHit(); return true; } if (PrimaryID != PrimID) ResetForNewPrimary(); //look in HC whether a hit with the same primID,UnitID,TSliceID exists already: G4bool found = false; for (G4int j=0; jentries()&&!found; j++) { CCalG4Hit* aPreviousHit = (*theHC)[j]; if (aPreviousHit->getTrackID() == PrimaryID && aPreviousHit->getTimeSliceID() == TSliceID && aPreviousHit->getUnitID()== UnitID ) { CurrentHit = aPreviousHit; found = true; } } if (found) { updateHit(); return true; } else { return false; } } void CCaloSD::ResetForNewPrimary() { EntrancePoint = SetToLocal(HitPoint); IncidentEnergy = PreStepPoint->GetKineticEnergy(); } void CCaloSD::StoreHit(CCalG4Hit* hit){ if (PrimID<0) return; if (hit == 0 ) { G4cout << "CCaloSD: hit to be stored is NULL !!" <insert( hit ); } void CCaloSD::createNewHit() { #ifdef debug G4int currentCopyNo = -999; G4int motherCopyNo = -999; G4TouchableHistory* theTouchable = (G4TouchableHistory*)( theTrack->GetTouchable() ); if ( theTouchable ) { currentCopyNo = theTouchable->GetReplicaNumber( 0 ); if ( theTouchable->GetHistoryDepth() > 0 ) { motherCopyNo = theTouchable->GetReplicaNumber( 1 ); } } G4cout << "CCaloSD createNewHit for" << " PV " << CurrentPV->GetName() << " PVid = " << currentCopyNo << " MVid = " << motherCopyNo << " Unit " << UnitID <GetTrackID() << " which is a " << theTrack->GetDefinition()->GetParticleName(); if (theTrack->GetTrackID()==1) { G4cout << " of energy " << theTrack->GetTotalEnergy(); } else { G4cout << " daughter of part. " << theTrack->GetParentID(); } G4cout << " and created by " ; if (theTrack->GetCreatorProcess()!=NULL) G4cout << theTrack->GetCreatorProcess()->GetProcessName() ; else G4cout << "NO process"; G4cout << G4endl; #endif CurrentHit = new CCalG4Hit; CurrentHit->setTrackID(PrimaryID); CurrentHit->setTimeSlice(TSlice); CurrentHit->setUnitID(UnitID); CurrentHit->setEntry(EntrancePoint); CurrentHit->setIncidentEnergy(IncidentEnergy); updateHit(); StoreHit(CurrentHit); } void CCaloSD::updateHit() { if (EdepositEM+EdepositEHAD != 0) { CurrentHit->addEnergyDeposit(EdepositEM,EdepositEHAD); #ifdef debug G4cout << "Energy deposit in Unit " << UnitID << " em " << EdepositEM/MeV << " hadronic " << EdepositEHAD/MeV << " MeV" << G4endl; #endif } // buffer for next steps: TSID = TSliceID; PrimID = PrimaryID; PreviousPV = CurrentPV; PreviousUnitID = UnitID; } G4ThreeVector CCaloSD::SetToLocal(G4ThreeVector global){ G4ThreeVector localPoint; const G4VTouchable* touch= PreStepPoint->GetTouchable(); localPoint=touch->GetHistory()->GetTopTransform().TransformPoint(global); return localPoint; } void CCaloSD::EndOfEvent(G4HCofThisEvent*) { summarize(); } void CCaloSD::summarize() { } void CCaloSD::clear() { } void CCaloSD::DrawAll() { } void CCaloSD::PrintAll() { G4cout << "CCaloSD: Collection " << theHC->GetName() << G4endl; theHC->PrintAllHits(); } void CCaloSD::SetOrganization(CCalVOrganization* org){ if (theDescription!=0) delete theDescription; theDescription = org; } G4double CCaloSD::curve_LY(G4StepPoint* stepPoint) { G4double weight = 1.; G4ThreeVector localPoint = SetToLocal(stepPoint->GetPosition()); G4double crlength = 230.; G4double dapd = 0.5 * crlength - localPoint.z(); if (dapd >= -0.1 || dapd <= crlength+0.1) { if (dapd <= 100.) weight = 1.05 - dapd * 0.0005; } else { G4cout << "CCaloSD, light coll curve : wrong distance to APD " << dapd << " crlength = " << crlength << " z of localPoint = " << localPoint.z() << " take weight = " << weight << G4endl; } #ifdef ddebug G4cout << "CCaloSD, light coll curve : " << dapd << " crlength = " << crlength << " z of localPoint = " << localPoint.z() << " take weight = " << weight << G4endl; #endif return weight; }