source: trunk/examples/advanced/composite_calorimeter/src/CCaloSD.cc @ 1304

//
// ********************************************************************
// * 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; j<theHC->entries()&&!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 !!" <<G4endl;
    return;
  }

  theHC->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 <<G4endl;
  G4cout << " primary "    << PrimaryID
         << " time slice " << TSliceID 
         << " For Track  " << theTrack->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;
}