Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/composite_calorimeter/src/CCaloSD.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 ///////////////////////////////////////////////////////////////////////////////
 27 // File: CCaloSD.cc
 28 // Description: Stores hits of calorimetric type in appropriate container
 29 ///////////////////////////////////////////////////////////////////////////////
 30 
 31 #include "CCaloSD.hh"
 32 #include "G4SystemOfUnits.hh"
 33 #include "G4VProcess.hh"
 34 #include "G4SDManager.hh"
 35 #include "G4VTouchable.hh"
 36 #include "CCalVOrganization.hh"
 37 #include "CCalSDList.hh"
 38 
 39 //#define debug
 40 //#define ddebug
 41  
 42 CCaloSD::CCaloSD(G4String name, CCalVOrganization* numberingScheme):
 43   G4VSensitiveDetector(name), HCID(-1), SDname(name), theHC(0),
 44   CurrentHit(0), theTrack(0), CurrentPV(0), PreviousPV(0), UnitID(0), 
 45   PreviousUnitID(0), PreStepPoint(0), PostStepPoint(0), 
 46   theDescription(numberingScheme) {
 47   
 48   collectionName.insert(name);
 49   
 50   G4cout << "*******************************************************" << G4endl;
 51   G4cout << "*                                                     *" << G4endl;
 52   G4cout << "* Constructing a CCaloSD  with name " << name            << G4endl;
 53   G4cout << "*                                                     *" << G4endl;
 54   G4cout << "*******************************************************" << G4endl;
 55 
 56   CCalSDList::getInstance()->addCalo(name);
 57 }
 58 
 59 CCaloSD::~CCaloSD() {
 60   delete theDescription;
 61 }
 62 
 63 void CCaloSD::Initialize(G4HCofThisEvent*HCE) {
 64 
 65 #ifdef debug
 66   G4cout << "CCaloSD : Initialize called for " << SDname << G4endl;
 67 #endif
 68   //This initialization is performed at the beginning of an event
 69   //------------------------------------------------------------
 70 
 71   theHC = new CCalG4HitCollection(SDname, collectionName[0]);
 72   if (HCID<0) {
 73     HCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
 74   }
 75   HCE->AddHitsCollection( HCID, theHC );
 76 
 77   TSID = -2;
 78   PrimID = -2;
 79   /////PrimaryID = -99;  <--- initialized by StackingAction.
 80 }
 81 
 82 G4bool CCaloSD::ProcessHits(G4Step*aStep, G4TouchableHistory*) {
 83 
 84   if(aStep->GetTotalEnergyDeposit() > CLHEP::eV) {
 85     getStepInfo(aStep);
 86     if (hitExists() == false && EdepositEM+EdepositEHAD>0.f) {
 87       createNewHit();
 88     }
 89   }
 90   return true;
 91 } 
 92 
 93 void CCaloSD::getStepInfo(const G4Step* aStep) {
 94   
 95   PreStepPoint = aStep->GetPreStepPoint(); 
 96   PostStepPoint= aStep->GetPostStepPoint(); 
 97   HitPoint     = PreStepPoint->GetPosition();        
 98 
 99   theTrack = aStep->GetTrack();   
100   CurrentPV= PreStepPoint->GetPhysicalVolume();
101 
102   G4String pname = theTrack->GetDefinition()->GetParticleName();
103   G4double de = aStep->GetTotalEnergyDeposit();
104   //G4cout << "##### Process new step dE= " << de
105   //         << "  " << pname << " inside " << GetName() << G4endl;
106 
107   const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
108   G4double weight = 1.;
109   if (touch->GetVolume(0)->GetName() == "CrystalMatrixCrystal") {
110     weight = curve_LY(PreStepPoint);
111   }
112 
113   if (pname == "e-" || pname == "e+" || pname == "gamma" ){
114     EdepositEM   = weight*de;
115     EdepositEHAD = 0.f;
116   } else {
117     EdepositEM   = 0.f;
118     EdepositEHAD = weight*de;
119   }
120 
121   TSlice = (PostStepPoint) ? PostStepPoint->GetGlobalTime()/nanosecond : 0.;
122   //G4cout << "     W= " << weight << " T= " << TSlice << G4endl;
123   G4int it = (G4int)TSlice;
124   TSliceID = std::min(it, 10000);
125   //G4cout << " tID= " <<  TSliceID << G4endl;
126 
127   UnitID = (theDescription) ? theDescription->GetUnitID(aStep) : 0;
128 }
129 
130 G4bool CCaloSD::hitExists() {
131    
132   if ( CurrentPV==PreviousPV && PrimaryID == PrimID && TSliceID == TSID &&
133        UnitID==PreviousUnitID) {
134     updateHit();
135     return true;
136   }
137    
138   if (PrimaryID != PrimID) { ResetForNewPrimary(); }
139    
140   // look in HC whether a hit with the same primID,UnitID,TSliceID 
141   // exists already:
142    
143   G4bool found = false;
144   for (std::size_t j=0; j<theHC->entries(); ++j) {
145 
146     CCalG4Hit* aPreviousHit = (*theHC)[j];
147     if (aPreviousHit->getTrackID()  == PrimaryID &&
148         aPreviousHit->getTimeSliceID() == TSliceID  &&
149         aPreviousHit->getUnitID()== UnitID       ) {
150       CurrentHit = aPreviousHit;
151       found = true;
152       break;
153     }
154   }          
155 
156   if (found) {
157     updateHit();
158     return true;
159   } else {
160     return false;
161   }    
162 }
163 
164 void CCaloSD::ResetForNewPrimary() {
165   
166   EntrancePoint = SetToLocal(HitPoint);
167   IncidentEnergy = PreStepPoint->GetKineticEnergy();
168 
169 }
170 
171 void CCaloSD::StoreHit(CCalG4Hit* hit){
172 
173   if (PrimID<0) return;
174   if (hit == 0 ) {
175     G4cout << "CCaloSD: hit to be stored is NULL !!" <<G4endl;
176     return;
177   }
178 
179   theHC->insert( hit );
180 }
181 
182 void CCaloSD::createNewHit() {
183 
184 #ifdef debug
185   G4int currentCopyNo = -999;
186   G4int motherCopyNo  = -999;
187   G4TouchableHistory* theTouchable
188     = (G4TouchableHistory*)( theTrack->GetTouchable() );
189   if ( theTouchable ) {
190     currentCopyNo = theTouchable->GetReplicaNumber( 0 );
191     if ( theTouchable->GetHistoryDepth() > 0 ) {
192       motherCopyNo = theTouchable->GetReplicaNumber( 1 );
193     }
194   }
195   G4cout << "CCaloSD createNewHit for"
196          << " PV "     << CurrentPV->GetName()
197          << " PVid = " << currentCopyNo
198          << " MVid = " << motherCopyNo
199          << " Unit "   << UnitID <<G4endl;
200   G4cout << " primary "    << PrimaryID
201          << " time slice " << TSliceID 
202          << " For Track  " << theTrack->GetTrackID()
203          << " which is a " <<  theTrack->GetDefinition()->GetParticleName();
204            
205   if (theTrack->GetTrackID()==1) {
206     G4cout << " of energy "     << theTrack->GetTotalEnergy();
207   } else {
208     G4cout << " daughter of part. " << theTrack->GetParentID();
209   }
210 
211   G4cout  << " and created by " ;
212   if (theTrack->GetCreatorProcess()!=NULL)
213     G4cout << theTrack->GetCreatorProcess()->GetProcessName() ;
214   else 
215     G4cout << "NO process";
216   G4cout << G4endl;
217 #endif          
218     
219   CurrentHit = new CCalG4Hit;
220   CurrentHit->setTrackID(PrimaryID);
221   CurrentHit->setTimeSlice(TSlice);
222   CurrentHit->setUnitID(UnitID);
223   CurrentHit->setEntry(EntrancePoint);
224   CurrentHit->setIncidentEnergy(IncidentEnergy);
225   updateHit();
226   
227   StoreHit(CurrentHit);
228 
229 }         
230 
231 void CCaloSD::updateHit() {
232   if (EdepositEM+EdepositEHAD != 0) {
233     CurrentHit->addEnergyDeposit(EdepositEM,EdepositEHAD);
234 #ifdef debug
235     G4cout << "Energy deposit in Unit " << UnitID << " em " << EdepositEM/MeV
236          << " hadronic " << EdepositEHAD/MeV << " MeV" << G4endl;
237 #endif
238   }
239 
240   // buffer for next steps:
241   TSID = TSliceID;
242   PrimID = PrimaryID;
243   PreviousPV = CurrentPV;
244   PreviousUnitID = UnitID;
245 }
246 
247 G4ThreeVector CCaloSD::SetToLocal(const G4ThreeVector& global) const{
248 
249   G4ThreeVector localPoint;
250   const G4VTouchable*   touch= PreStepPoint->GetTouchable();
251   localPoint=touch->GetHistory()->GetTopTransform().TransformPoint(global);
252   
253   return localPoint;  
254 }
255      
256 void CCaloSD::EndOfEvent(G4HCofThisEvent*) {
257   summarize();
258 }
259      
260 void CCaloSD::summarize() {
261 }
262 
263 void CCaloSD::clear() {
264 } 
265 
266 void CCaloSD::DrawAll() {
267 } 
268 
269 void CCaloSD::PrintAll() {
270   G4cout << "CCaloSD: Collection " << theHC->GetName() << G4endl;
271   theHC->PrintAllHits();
272 } 
273 
274 void CCaloSD::SetOrganization(CCalVOrganization* org){
275 
276   if (theDescription!=0) 
277     delete theDescription;
278   theDescription = org;
279 }
280 
281 G4double CCaloSD::curve_LY(const G4StepPoint* stepPoint) {
282 
283   G4double weight = 1.;
284   G4ThreeVector localPoint = SetToLocal(stepPoint->GetPosition());
285   const G4double crlength = 230.;
286   G4double dapd = 0.5 * crlength - localPoint.z();
287   if (dapd >= -0.1 || dapd <= crlength+0.1) {
288     if (dapd <= 100.)
289       weight = 1.05 - dapd * 0.0005;
290   } else {
291     G4cout << "CCaloSD, light coll curve : wrong distance to APD " << dapd 
292            << " crlength = " << crlength
293            << " z of localPoint = " << localPoint.z() 
294            << " take weight = " << weight << G4endl;
295   }
296 #ifdef ddebug
297   G4cout << "CCaloSD, light coll curve : " << dapd 
298          << " crlength = " << crlength
299          << " z of localPoint = " << localPoint.z() 
300          << " take weight = " << weight << G4endl;
301 #endif
302   return weight;
303 }
304