composite_calorimeter/src/CCaloSD.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 ///////////////////////////////////////////////////////////////////////////////
0027 // File: CCaloSD.cc
0028 // Description: Stores hits of calorimetric type in appropriate container
0029 ///////////////////////////////////////////////////////////////////////////////
0030
0031 #include "CCaloSD.hh"
0032 #include "G4SystemOfUnits.hh"
0033 #include "G4VProcess.hh"
0034 #include "G4SDManager.hh"
0035 #include "G4VTouchable.hh"
0036 #include "CCalVOrganization.hh"
0037 #include "CCalSDList.hh"
0038
0039 //#define debug
0040 //#define ddebug
0041
0042 CCaloSD::CCaloSD(G4String name, CCalVOrganization* numberingScheme):
0043   G4VSensitiveDetector(name), HCID(-1), SDname(name), theHC(0),
0044   CurrentHit(0), theTrack(0), CurrentPV(0), PreviousPV(0), UnitID(0),
0045   PreviousUnitID(0), PreStepPoint(0), PostStepPoint(0),
0046   theDescription(numberingScheme) {
0047
0048   collectionName.insert(name);
0049
0050   G4cout << "*******************************************************" << G4endl;
0051   G4cout << "*                                                     *" << G4endl;
0052   G4cout << "* Constructing a CCaloSD  with name " << name            << G4endl;
0053   G4cout << "*                                                     *" << G4endl;
0054   G4cout << "*******************************************************" << G4endl;
0055
0056   CCalSDList::getInstance()->addCalo(name);
0057 }
0058
0059 CCaloSD::~CCaloSD() {
0060   delete theDescription;
0061 }
0062
0063 void CCaloSD::Initialize(G4HCofThisEvent*HCE) {
0064
0065 #ifdef debug
0066   G4cout << "CCaloSD : Initialize called for " << SDname << G4endl;
0067 #endif
0068   //This initialization is performed at the beginning of an event
0069   //------------------------------------------------------------
0070
0071   theHC = new CCalG4HitCollection(SDname, collectionName[0]);
0072   if (HCID<0) {
0073     HCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
0074   }
0075   HCE->AddHitsCollection( HCID, theHC );
0076
0077   TSID = -2;
0078   PrimID = -2;
0079   /////PrimaryID = -99;  <--- initialized by StackingAction.
0080 }
0081
0082 G4bool CCaloSD::ProcessHits(G4Step*aStep, G4TouchableHistory*) {
0083
0084   if(aStep->GetTotalEnergyDeposit() > CLHEP::eV) {
0085     getStepInfo(aStep);
0086     if (hitExists() == false && EdepositEM+EdepositEHAD>0.f) {
0087       createNewHit();
0088     }
0089   }
0090   return true;
0091 }
0092
0093 void CCaloSD::getStepInfo(const G4Step* aStep) {
0094
0095   PreStepPoint = aStep->GetPreStepPoint();
0096   PostStepPoint= aStep->GetPostStepPoint();
0097   HitPoint     = PreStepPoint->GetPosition();
0098
0099   theTrack = aStep->GetTrack();
0100   CurrentPV= PreStepPoint->GetPhysicalVolume();
0101
0102   G4String pname = theTrack->GetDefinition()->GetParticleName();
0103   G4double de = aStep->GetTotalEnergyDeposit();
0104   //G4cout << "##### Process new step dE= " << de
0105   //         << "  " << pname << " inside " << GetName() << G4endl;
0106
0107   const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
0108   G4double weight = 1.;
0109   if (touch->GetVolume(0)->GetName() == "CrystalMatrixCrystal") {
0110     weight = curve_LY(PreStepPoint);
0111   }
0112
0113   if (pname == "e-" || pname == "e+" || pname == "gamma" ){
0114     EdepositEM   = weight*de;
0115     EdepositEHAD = 0.f;
0116   } else {
0117     EdepositEM   = 0.f;
0118     EdepositEHAD = weight*de;
0119   }
0120
0121   TSlice = (PostStepPoint) ? PostStepPoint->GetGlobalTime()/nanosecond : 0.;
0122   //G4cout << "     W= " << weight << " T= " << TSlice << G4endl;
0123   G4int it = (G4int)TSlice;
0124   TSliceID = std::min(it, 10000);
0125   //G4cout << " tID= " <<  TSliceID << G4endl;
0126
0127   UnitID = (theDescription) ? theDescription->GetUnitID(aStep) : 0;
0128 }
0129
0130 G4bool CCaloSD::hitExists() {
0131
0132   if ( CurrentPV==PreviousPV && PrimaryID == PrimID && TSliceID == TSID &&
0133        UnitID==PreviousUnitID) {
0134     updateHit();
0135     return true;
0136   }
0137
0138   if (PrimaryID != PrimID) { ResetForNewPrimary(); }
0139
0140   // look in HC whether a hit with the same primID,UnitID,TSliceID
0141   // exists already:
0142
0143   G4bool found = false;
0144   for (std::size_t j=0; j<theHC->entries(); ++j) {
0145
0146     CCalG4Hit* aPreviousHit = (*theHC)[j];
0147     if (aPreviousHit->getTrackID()  == PrimaryID &&
0148         aPreviousHit->getTimeSliceID() == TSliceID  &&
0149         aPreviousHit->getUnitID()== UnitID       ) {
0150       CurrentHit = aPreviousHit;
0151       found = true;
0152       break;
0153     }
0154   }
0155
0156   if (found) {
0157     updateHit();
0158     return true;
0159   } else {
0160     return false;
0161   }
0162 }
0163
0164 void CCaloSD::ResetForNewPrimary() {
0165
0166   EntrancePoint = SetToLocal(HitPoint);
0167   IncidentEnergy = PreStepPoint->GetKineticEnergy();
0168
0169 }
0170
0171 void CCaloSD::StoreHit(CCalG4Hit* hit){
0172
0173   if (PrimID<0) return;
0174   if (hit == 0 ) {
0175     G4cout << "CCaloSD: hit to be stored is NULL !!" <<G4endl;
0176     return;
0177   }
0178
0179   theHC->insert( hit );
0180 }
0181
0182 void CCaloSD::createNewHit() {
0183
0184 #ifdef debug
0185   G4int currentCopyNo = -999;
0186   G4int motherCopyNo  = -999;
0187   G4TouchableHistory* theTouchable
0188     = (G4TouchableHistory*)( theTrack->GetTouchable() );
0189   if ( theTouchable ) {
0190     currentCopyNo = theTouchable->GetReplicaNumber( 0 );
0191     if ( theTouchable->GetHistoryDepth() > 0 ) {
0192       motherCopyNo = theTouchable->GetReplicaNumber( 1 );
0193     }
0194   }
0195   G4cout << "CCaloSD createNewHit for"
0196          << " PV "     << CurrentPV->GetName()
0197          << " PVid = " << currentCopyNo
0198          << " MVid = " << motherCopyNo
0199          << " Unit "   << UnitID <<G4endl;
0200   G4cout << " primary "    << PrimaryID
0201          << " time slice " << TSliceID
0202          << " For Track  " << theTrack->GetTrackID()
0203          << " which is a " <<  theTrack->GetDefinition()->GetParticleName();
0204
0205   if (theTrack->GetTrackID()==1) {
0206     G4cout << " of energy "     << theTrack->GetTotalEnergy();
0207   } else {
0208     G4cout << " daughter of part. " << theTrack->GetParentID();
0209   }
0210
0211   G4cout  << " and created by " ;
0212   if (theTrack->GetCreatorProcess()!=NULL)
0213     G4cout << theTrack->GetCreatorProcess()->GetProcessName() ;
0214   else
0215     G4cout << "NO process";
0216   G4cout << G4endl;
0217 #endif
0218
0219   CurrentHit = new CCalG4Hit;
0220   CurrentHit->setTrackID(PrimaryID);
0221   CurrentHit->setTimeSlice(TSlice);
0222   CurrentHit->setUnitID(UnitID);
0223   CurrentHit->setEntry(EntrancePoint);
0224   CurrentHit->setIncidentEnergy(IncidentEnergy);
0225   updateHit();
0226
0227   StoreHit(CurrentHit);
0228
0229 }
0230
0231 void CCaloSD::updateHit() {
0232   if (EdepositEM+EdepositEHAD != 0) {
0233     CurrentHit->addEnergyDeposit(EdepositEM,EdepositEHAD);
0234 #ifdef debug
0235     G4cout << "Energy deposit in Unit " << UnitID << " em " << EdepositEM/MeV
0236          << " hadronic " << EdepositEHAD/MeV << " MeV" << G4endl;
0237 #endif
0238   }
0239
0240   // buffer for next steps:
0241   TSID = TSliceID;
0242   PrimID = PrimaryID;
0243   PreviousPV = CurrentPV;
0244   PreviousUnitID = UnitID;
0245 }
0246
0247 G4ThreeVector CCaloSD::SetToLocal(const G4ThreeVector& global) const{
0248
0249   G4ThreeVector localPoint;
0250   const G4VTouchable*   touch= PreStepPoint->GetTouchable();
0251   localPoint=touch->GetHistory()->GetTopTransform().TransformPoint(global);
0252
0253   return localPoint;
0254 }
0255
0256 void CCaloSD::EndOfEvent(G4HCofThisEvent*) {
0257   summarize();
0258 }
0259
0260 void CCaloSD::summarize() {
0261 }
0262
0263 void CCaloSD::clear() {
0264 }
0265
0266 void CCaloSD::DrawAll() {
0267 }
0268
0269 void CCaloSD::PrintAll() {
0270   G4cout << "CCaloSD: Collection " << theHC->GetName() << G4endl;
0271   theHC->PrintAllHits();
0272 }
0273
0274 void CCaloSD::SetOrganization(CCalVOrganization* org){
0275
0276   if (theDescription!=0)
0277     delete theDescription;
0278   theDescription = org;
0279 }
0280
0281 G4double CCaloSD::curve_LY(const G4StepPoint* stepPoint) {
0282
0283   G4double weight = 1.;
0284   G4ThreeVector localPoint = SetToLocal(stepPoint->GetPosition());
0285   const G4double crlength = 230.;
0286   G4double dapd = 0.5 * crlength - localPoint.z();
0287   if (dapd >= -0.1 || dapd <= crlength+0.1) {
0288     if (dapd <= 100.)
0289       weight = 1.05 - dapd * 0.0005;
0290   } else {
0291     G4cout << "CCaloSD, light coll curve : wrong distance to APD " << dapd
0292            << " crlength = " << crlength
0293            << " z of localPoint = " << localPoint.z()
0294            << " take weight = " << weight << G4endl;
0295   }
0296 #ifdef ddebug
0297   G4cout << "CCaloSD, light coll curve : " << dapd
0298          << " crlength = " << crlength
0299          << " z of localPoint = " << localPoint.z()
0300          << " take weight = " << weight << G4endl;
0301 #endif
0302   return weight;
0303 }