Sensitive Detector implementation of the CSI detector. More...

#include <SensitiveDetector.h>

Inheritance diagram for SensitiveDetector:

Public Member Functions
	SensitiveDetector ()
	Constructor.

	~SensitiveDetector ()
	Destructor.

int	saveSimHit (const G4int cellId, const G4int trackID, const G4int pid, const G4double tof, const G4double edep, G4ThreeVector mom, G4ThreeVector pos, G4ThreeVector PosCell, G4ThreeVector VecCell)
	Save ECLSimHits of the event into datastore.

Static Public Member Functions
static const std::map< std::string, RelationArray::EConsolidationAction > &	getMCParticleRelations ()
	Return a list of all registered Relations with MCParticles.

static void	setActive (bool activeStatus)
	Enable/Disable all Sensitive Detectors.

static void	registerMCParticleRelation (const std::string &name, RelationArray::EConsolidationAction ignoreAction=RelationArray::c_negativeWeight)
	Register an relation involving MCParticles.

static void	registerMCParticleRelation (const RelationArray &relation, RelationArray::EConsolidationAction ignoreAction=RelationArray::c_negativeWeight)
	Overload to make it easer to register MCParticle relations.

Protected Member Functions
bool	step (G4Step aStep, G4TouchableHistory ) override
	Step processing method.

Private Member Functions
virtual bool	ProcessHits (G4Step aStep, G4TouchableHistory aROhist)
	Check if recording hits is enabled and if so call step() and set the correct MCParticle flag.

Private Attributes
int	m_hitNum
	members of SensitiveDetector

int	m_EvnetNumber
	The current number of created hits in an event.

int	m_oldEvnetNumber
	The current number of created hits in an event.

int	m_trackID
	The current number of created hits in an event.

G4ThreeVector	m_startPos
	track id

G4ThreeVector	m_endPos
	Position of prestep.

G4ThreeVector	m_WightedPos
	Position of poststep.

G4ThreeVector	m_momentum
	Wighted step Position.

double	m_startTime
	momentum of track

double	m_endTime
	global time

double	m_WightedTime
	global time

double	m_startEnergy
	global time

double	m_energyDeposit
	particle energy at the entrance in volume

double	m_trackLength
	energy deposited in volume

int	CsiHitIndex [18][80]
	length of the track in the volume

int	iECLCell
	Hit index of StoreArray.

int	TimeIndex
	Hit Energy of StoreArray.

double	local_pos
	Hit Time of StoreArray.

double	T_ave
	position alongthe vector of crystal axis

int	firstcall
	flight time to diode sensor

int	m_phiID
	flag of first call

int	m_thetaID
	The current phi ID in an event.

int	m_cellID
	The current theta ID in an event.

Const::EDetector	m_subdetector
	Subdetector the class belongs to.

Static Private Attributes
static std::map< std::string, RelationArray::EConsolidationAction >	s_mcRelations
	Static set holding all relations which have to be updated at the end of the Event.

static bool	s_active
	Static bool which indicates wether recording of hits is enabled.

Detailed Description

Sensitive Detector implementation of the CSI detector.

Definition at line 22 of file SensitiveDetector.h.

Constructor & Destructor Documentation

◆ SensitiveDetector()

SensitiveDetector ( )

Constructor.

Definition at line 29 of file SensitiveDetector.cc.

                                        :
      Simulation::SensitiveDetectorBase("CsiSensitiveDetector", Const::invalidDetector)
    {
      m_hitNum = 0;
      m_EvnetNumber = 0;
      m_oldEvnetNumber = 0;
      m_trackID = 0;
      m_startTime = 0;
      m_endTime = 0;
      m_WightedTime = 0;
      m_startEnergy = 0;
      m_energyDeposit = 0;
      m_trackLength = 0;
      iECLCell = 0;
      TimeIndex = 0;
      local_pos = 0;
      T_ave = 0;
      firstcall = 0;
      m_phiID = 0;
      m_thetaID = 0;
      m_cellID = 0;
 
 
      // CsiHitIndex: array of 18 (Nb. crystals) x 80 (max time index (tof [ns]/100))
      for (int i = 0; i < 18; i++) {for (int j = 0; j < 80; j++) { CsiHitIndex[i][j] = 0;}}
 
      //Make sure all collections are registered
      StoreArray<MCParticle>   mcParticles;
      StoreArray<CsiSimHit>  simHits;
      StoreArray<CsiHit> csiHits;
 
      RelationArray simHitRel(mcParticles, simHits);
      RelationArray csiHitRel(mcParticles, csiHits);
 
      //Register all collections we want to modify and require those we want to use
      mcParticles.registerInDataStore();
      simHits.registerInDataStore();
      csiHits.registerInDataStore();
      simHitRel.registerInDataStore();
      csiHitRel.registerInDataStore();
 
      //Register the Relation so that the TrackIDs get replaced by the actual
      //MCParticle indices after simulating the events. This is needed as
      //secondary particles might not be stored so everything relating to those
      //particles will be attributed to the last saved mother particle
      registerMCParticleRelation(simHitRel);
      registerMCParticleRelation(csiHitRel);
    }

◆ ~SensitiveDetector()

~SensitiveDetector ( )

Destructor.

Definition at line 78 of file SensitiveDetector.cc.

    {
 
    }

Member Function Documentation

◆ getMCParticleRelations()

static const std::map< std::string, RelationArray::EConsolidationAction > & getMCParticleRelations ( )

inlinestaticinherited

Return a list of all registered Relations with MCParticles.

Definition at line 42 of file SensitiveDetectorBase.h.

42{ return s_mcRelations; }

Belle2::Simulation::SensitiveDetectorBase::s_mcRelations

static std::map< std::string, RelationArray::EConsolidationAction > s_mcRelations

Static set holding all relations which have to be updated at the end of the Event.

Definition: SensitiveDetectorBase.h:87

◆ ProcessHits()

bool ProcessHits	(	G4Step *	aStep,
		G4TouchableHistory *	aROhist
	)

inlineprivatevirtualinherited

Check if recording hits is enabled and if so call step() and set the correct MCParticle flag.

Called by Geant4 for each step inside the sensitive volumes attached

Definition at line 94 of file SensitiveDetectorBase.h.

    {
      if (!s_active) return false;
      bool result = step(aStep, aROhist);
      // Do not include hits from invalid detector (beast,teastbeam, etc.)
      if (result && (m_subdetector != Const::invalidDetector)) TrackInfo::getInfo(*aStep->GetTrack()).addSeenInDetector(m_subdetector);
      return result;
    }

◆ registerMCParticleRelation() [1/2]

static void registerMCParticleRelation	(	const RelationArray &	relation,
		RelationArray::EConsolidationAction	ignoreAction = `RelationArray::c_negativeWeight`
	)

inlinestaticinherited

Overload to make it easer to register MCParticle relations.

Parameters

relation	RelationArray to register
ignoreAction

Definition at line 66 of file SensitiveDetectorBase.h.

67 { registerMCParticleRelation(relation.getName(), ignoreAction); }

◆ registerMCParticleRelation() [2/2]

void registerMCParticleRelation	(	const std::string &	name,
		RelationArray::EConsolidationAction	ignoreAction = `RelationArray::c_negativeWeight`
	)

staticinherited

Register an relation involving MCParticles.

All Relations which point from an MCParticle to something have to be registered with addMCParticleRelation() because the index of the MCParticles might change at the end of the event. During simulation, the TrackID should be used as index of the MCParticle

Parameters

name	Name of the relation to register
ignoreAction

Definition at line 22 of file SensitiveDetectorBase.cc.

    {
      std::pair<std::map<std::string, RelationArray::EConsolidationAction>::iterator, bool> insert = s_mcRelations.insert(std::make_pair(
            name, ignoreAction));
      //If the relation already exists and the ignoreAction is different we do have a problem
      if (!insert.second && insert.first->second != ignoreAction) {
        B2FATAL("MCParticle Relation " << name << " already registered with different ignore action.");
      }
    }

◆ saveSimHit()

int saveSimHit	(	const G4int	cellId,
		const G4int	trackID,
		const G4int	pid,
		const G4double	tof,
		const G4double	edep,
		G4ThreeVector	mom,
		G4ThreeVector	pos,
		G4ThreeVector	PosCell,
		G4ThreeVector	VecCell
	)

Save ECLSimHits of the event into datastore.

Parameters

cellId	the channel identifier
trackID	the track identifier
pid	the pdg code of the particle
tof	the time of flight of the particle
edep	the energy deposited in the crystal by the particle
mom	the momentum of the particle
pos	TO BE COMPLETED
PosCell	TO BE COMPLETED
VecCell	TO BE COMPLETED

Returns: int the number of hits saved in this event

Definition at line 153 of file SensitiveDetector.cc.

    {
 
      //Get the datastore arrays
      StoreArray<MCParticle>   mcParticles;
      StoreArray<CsiSimHit>    simHits;
 
      RelationArray simHitRel(mcParticles, simHits);
 
      // convert Hep3Vectors to ROOT::Math::XYZVector
      ROOT::Math::XYZVector momentum(mom.getX(), mom.getY(), mom.getZ());
      ROOT::Math::XYZVector position(pos.getX(), pos.getY(), pos.getZ());
 
      //Append
      simHits.appendNew(cellId, trackID, pid, tof / CLHEP::ns, edep / CLHEP::GeV,
                        momentum * (1.0 / CLHEP::GeV), position * (1.0 / CLHEP::cm));
 
      int simhitNumber = simHits.getEntries() - 1;
      simHitRel.add(trackID, simhitNumber);
      B2DEBUG(150, "Saving CsiSimHit Number: " << simhitNumber);
 
      // Store CsiHit: total energy deposited per hit, with average time.
 
      // Old names:
      /*
      StoreArray<CsiHit> eclHitArray;
      RelationArray eclHitRel(mcParticles, eclHitArray);
      */
 
      // New names:
 
      StoreArray<CsiHit> csiHits;
 
      RelationArray csiHitRel(mcParticles, csiHits);
 
 
      StoreObjPtr<EventMetaData> eventMetaDataPtr;
      int m_currentEvnetNumber = eventMetaDataPtr->getEvent();
 
 
      if (firstcall == 0 || m_currentEvnetNumber != m_oldEvnetNumber) {
        m_oldEvnetNumber = m_currentEvnetNumber;
        for (int iCSICell = 0; iCSICell < 18; iCSICell++) {
          for (int  TimeInd = 0; TimeInd < 80; TimeInd++) {
            CsiHitIndex[iCSICell][TimeInd] = -1;
          }
        }
        firstcall++;
      }
 
 
      if (m_currentEvnetNumber == m_oldEvnetNumber) {
        if ((tof / CLHEP::ns) < 8000) {
          TimeIndex = (int)((tof / CLHEP::ns) / 100);
 
          // local_pos is the distance between the hit and the PIN-diode end of the crystal.
          local_pos = (15. - (pos  - PosCell) * VecCell);
 
          // From the ECL package T_ave =  6.05 + 0.0749 * local_pos - 0.00112 * local_pos * local_pos + (tof / CLHEP::ns)  ;
 
          //For CsI: we remove the offset (I hypothetize it if from photo-sensor delay)
          //         we remove the acceleration of the light
          //         we correct the velocity term using n_CsI=1.80
          T_ave = 0.0600 * local_pos + (tof / CLHEP::ns)  ;
 
 
          double E_cell = (edep / CLHEP::GeV);
          if (CsiHitIndex[cellId][TimeIndex] == -1) {
            //m_hitNum = csiHits->GetLast() + 1;
            //new(csiHits->AddrAt(m_hitNum)) CSIHit();
            csiHits.appendNew();
            m_hitNum = csiHits.getEntries() - 1;
            CsiHitIndex[cellId][TimeIndex] = m_hitNum;
            csiHits[m_hitNum]->setCellId(cellId);
            csiHits[m_hitNum]->setEnergyDep(E_cell);
            csiHits[m_hitNum]->setTimeAve(T_ave);
            csiHits[m_hitNum]->setTimeVar(0);
            csiHitRel.add(trackID, m_hitNum);
 
          } else {
            m_hitNum = CsiHitIndex[cellId][TimeIndex];
            double old_edep = csiHits[m_hitNum]->getEnergyDep();
            double old_TimeAve = csiHits[m_hitNum]->getTimeAve();
            double old_TimeVar = csiHits[m_hitNum]->getTimeVar();
 
            //double new_edep = E_cell + old_edep;
            //double new_TimeAve = (old_edep * old_TimeAve + E_cell * T_ave) / new_edep;
 
            double temp = E_cell + old_edep;
            double delta = T_ave - old_TimeAve;
            double R = delta * E_cell / temp;
            double new_TimeAve = old_TimeAve + R;
            double new_TimeVar = old_TimeVar + old_edep * delta * R;
            double new_edep = temp;
 
            //Algo from Wiki:
            /*
              for x, weight in dataWeightPairs:  # Alternatively "for x, weight in zip(data, weights):"
                  temp = weight + sumweight
                delta = x - mean
                R = delta * weight / temp
                mean = mean + R
                M2 = M2 + sumweight * delta * R  # Alternatively, "M2 = M2 + weight * delta * (x−mean)"
                sumweight = temp
 
              variance_n = M2/sumweight
            */
 
            csiHits[m_hitNum]->setEnergyDep(new_edep);
            csiHits[m_hitNum]->setTimeAve(new_TimeAve);
            csiHits[m_hitNum]->setTimeVar(new_TimeVar);
            csiHitRel.add(trackID, m_hitNum);
 
          }
        }
      } else { B2ERROR("m_currentEvnetNumber  ERROR: m_oldEvnetNumber==m_oldEvnetNumber"); }
 
      return (simhitNumber);
    }//saveSimHit

◆ setActive()

static void setActive ( bool activeStatus )

inlinestaticinherited

Enable/Disable all Sensitive Detectors.

By default, all sensitive detectors won't create hits to make it possible to use the Geant4 Navigator for non-simulation purposes. Only during simulation the sensitive detectors will be enabled to record hits

Parameters

activeStatus bool to indicate wether hits should be recorded

Definition at line 50 of file SensitiveDetectorBase.h.

50{ s_active = activeStatus; }

◆ step()

bool step	(	G4Step *	aStep,
		G4TouchableHistory *
	)

overrideprotectedvirtual

Step processing method.

Parameters

aStep the G4Step with the current step information

Returns: true if a Hit has been created, false if the hit was ignored

Implements SensitiveDetectorBase.

Definition at line 85 of file SensitiveDetector.cc.

    {
 
      const G4StepPoint& preStep  = *aStep->GetPreStepPoint();
      const G4StepPoint& postStep = * aStep->GetPostStepPoint();
 
      G4Track& track  = * aStep->GetTrack();
      if (m_trackID != track.GetTrackID()) {
        //TrackID changed, store track informations
        m_trackID = track.GetTrackID();
        //Get momentum
        m_momentum = preStep.GetMomentum() ;
        //Get energy
        m_startEnergy =  preStep.GetKineticEnergy() ;
        //Reset energy deposit;
        m_energyDeposit = 0;
        //Reset Wighted Time;
        m_WightedTime = 0;
        //Reset m_WightedPos;
        m_WightedPos.set(0, 0, 0);
 
      }
      //Update energy deposit
      m_energyDeposit += aStep->GetTotalEnergyDeposit() ;
 
      m_startTime = preStep.GetGlobalTime();
      m_endTime = postStep.GetGlobalTime();
      m_WightedTime += (m_startTime + m_endTime) / 2 * (aStep->GetTotalEnergyDeposit());
 
      m_startPos =  preStep.GetPosition();
      m_endPos = postStep.GetPosition();
 
      G4ThreeVector position(0.5 * (m_startPos + m_endPos));
      m_WightedPos += position * (aStep->GetTotalEnergyDeposit());
 
      //Save Hit if track leaves volume or is killed
      if (track.GetNextVolume() != track.GetVolume() || track.GetTrackStatus() >= fStopAndKill) {
        int pdgCode = track.GetDefinition()->GetPDGEncoding();
 
        const G4VPhysicalVolume& v = * track.GetVolume();
 
        //Column z of the rotation matrix R is the result of R * zhat.
        // This is what we want since G4Trap are created aligned along z
        G4RotationMatrix rotCell = * v.GetObjectRotation();
        //        G4RotationMatrix rotCell = * v.GetRotation();
        G4ThreeVector vecCell = rotCell.colZ();
        G4ThreeVector posCell = v.GetTranslation();
 
        // Get layer ID
        if (v.GetName().find("Crystal") != std::string::npos) {
          CsiGeometryPar* csip = CsiGeometryPar::Instance();
          m_cellID = csip->CsiVolNameToCellID(v.GetName());
 
          double dTotalEnergy = 1.0 / m_energyDeposit; //avoid the error  no match for 'operator/'
          if (m_energyDeposit > 0.) {
            saveSimHit(m_cellID, m_trackID, pdgCode, m_WightedTime * dTotalEnergy,
                       m_energyDeposit, m_momentum, m_WightedPos * dTotalEnergy,
                       posCell, vecCell);
          }
        }
 
        //Reset TrackID
        m_trackID = 0;
      }
 
      return true;
    }

Member Data Documentation

◆ CsiHitIndex

int CsiHitIndex[18][80]

private

length of the track in the volume

Definition at line 83 of file SensitiveDetector.h.

◆ firstcall

int firstcall

private

flight time to diode sensor

Definition at line 90 of file SensitiveDetector.h.

◆ iECLCell

int iECLCell

private

Hit index of StoreArray.

Definition at line 84 of file SensitiveDetector.h.

◆ local_pos

double local_pos

private

Hit Time of StoreArray.

Definition at line 88 of file SensitiveDetector.h.

◆ m_cellID

int m_cellID

private

The current theta ID in an event.

Used to fill the DataStore ECL array

Definition at line 93 of file SensitiveDetector.h.

◆ m_endPos

G4ThreeVector m_endPos

private

Position of prestep.

Definition at line 73 of file SensitiveDetector.h.

◆ m_endTime

double m_endTime

private

global time

Definition at line 77 of file SensitiveDetector.h.

◆ m_energyDeposit

double m_energyDeposit

private

particle energy at the entrance in volume

Definition at line 80 of file SensitiveDetector.h.

◆ m_EvnetNumber

int m_EvnetNumber

private

The current number of created hits in an event.

Used to fill the DataStore ECLHit.

Definition at line 69 of file SensitiveDetector.h.

◆ m_hitNum

int m_hitNum

private

members of SensitiveDetector

Definition at line 68 of file SensitiveDetector.h.

◆ m_momentum

G4ThreeVector m_momentum

private

Wighted step Position.

Definition at line 75 of file SensitiveDetector.h.

◆ m_oldEvnetNumber

int m_oldEvnetNumber

private

The current number of created hits in an event.

Used to fill the DataStore ECL EB array.

Definition at line 70 of file SensitiveDetector.h.

◆ m_phiID

int m_phiID

private

flag of first call

Definition at line 91 of file SensitiveDetector.h.

◆ m_startEnergy

double m_startEnergy

private

global time

Definition at line 79 of file SensitiveDetector.h.

◆ m_startPos

G4ThreeVector m_startPos

private

track id

Definition at line 72 of file SensitiveDetector.h.

◆ m_startTime

double m_startTime

private

momentum of track

Definition at line 76 of file SensitiveDetector.h.

◆ m_subdetector

Const::EDetector m_subdetector

privateinherited

Subdetector the class belongs to.

Definition at line 91 of file SensitiveDetectorBase.h.

◆ m_thetaID

int m_thetaID

private

The current phi ID in an event.

Used to fill the DataStore ECL array

Definition at line 92 of file SensitiveDetector.h.

◆ m_trackID

int m_trackID

private

The current number of created hits in an event.

Used to fill the DataStore

Definition at line 71 of file SensitiveDetector.h.

◆ m_trackLength

double m_trackLength

private

energy deposited in volume

Definition at line 81 of file SensitiveDetector.h.

◆ m_WightedPos

G4ThreeVector m_WightedPos

private

Position of poststep.

Definition at line 74 of file SensitiveDetector.h.

◆ m_WightedTime

double m_WightedTime

private

global time

Definition at line 78 of file SensitiveDetector.h.

◆ s_active

bool s_active

staticprivateinherited

Static bool which indicates wether recording of hits is enabled.

Definition at line 89 of file SensitiveDetectorBase.h.

◆ s_mcRelations

map< string, RelationArray::EConsolidationAction > s_mcRelations

staticprivateinherited

Static set holding all relations which have to be updated at the end of the Event.

Definition at line 87 of file SensitiveDetectorBase.h.

◆ T_ave

double T_ave

private

position alongthe vector of crystal axis

Definition at line 89 of file SensitiveDetector.h.

◆ TimeIndex

int TimeIndex

private

Hit Energy of StoreArray.

Definition at line 85 of file SensitiveDetector.h.

The documentation for this class was generated from the following files:

beast/csi/simulation/include/SensitiveDetector.h
beast/csi/simulation/src/SensitiveDetector.cc

Public Member Functions

Static Public Member Functions

Protected Member Functions

Private Member Functions

Private Attributes

Static Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ SensitiveDetector()

◆ ~SensitiveDetector()

Member Function Documentation

◆ getMCParticleRelations()

◆ ProcessHits()

◆ registerMCParticleRelation() [1/2]

◆ registerMCParticleRelation() [2/2]

◆ saveSimHit()

◆ setActive()

◆ step()

Member Data Documentation

◆ CsiHitIndex

◆ firstcall

◆ iECLCell

◆ local_pos

◆ m_cellID

◆ m_endPos

◆ m_endTime

◆ m_energyDeposit

◆ m_EvnetNumber

◆ m_hitNum

◆ m_momentum

◆ m_oldEvnetNumber

◆ m_phiID

◆ m_startEnergy

◆ m_startPos

◆ m_startTime

◆ m_subdetector

◆ m_thetaID

◆ m_trackID

◆ m_trackLength

◆ m_WightedPos

◆ m_WightedTime

◆ s_active

◆ s_mcRelations

◆ T_ave

◆ TimeIndex