Exercise #3 -- Sensitive Detectors

Introduction

This tutorial has been updated for use at the October 2003 Fermilab Geant4 Tutorial, but is also designed for use by anyone else who wants to learn defining geometrical setup in Geant4 application.

To start with this try out, you should start with provided copy TryOut03.tar.gz. Once you succeed to compile/link, just try it to see how it runs. 

  $ $G4WORKDIR/bin/Linux-g++/A01app
After some initialization, a prompt Idle> should appear. Then, adopt a UI command for proper visualization, and then execute an event. 
  Idle> /run/beamOn 1
You should get this drawing.

On the display, you should get a summary of this event similar to this

An event summary 
>>> Event 0 >>> Simulation truth : e+ (-4.61596,0,989.293)
Hodoscope 1 has 1 hits.
  Hodoscope[7] 4.98684 (nsec)
Drift Chamber 1 has 5 hits.
  Layer[0] : time 6.67102 (nsec) --- local (x,y) -8.96713, 0.263021
  Layer[1] : time 8.33886 (nsec) --- local (x,y) -11.1916, -0.0772413
  Layer[2] : time 10.0067 (nsec) --- local (x,y) -13.3434, -0.591037
  Layer[3] : time 11.6745 (nsec) --- local (x,y) -15.5842, -1.03141
  Layer[4] : time 13.3424 (nsec) --- local (x,y) -17.7618, -1.50829

As you see, you do not have any output for detectors in the second lever arm. In this exercise, we are working for this.

Setting sensitivity to Hodoscope and Drift Chamber

We already have sensitive detector classes for hodoscope and drift chamber. So, for these detectors in the second lever arm, what you have to do are
  1. Set sensitive detector class objects to corresponding logical volumes. For drift chamber, logical volume which should be sensitive is not the chamber volume but the virtual wire plane.
  2. Add print-out routine for them in EndOfEventAction.

A part of A01DetectorConstruction.cc 
G4VPhysicalVolume* A01DetectorConstruction::Construct()
{
  // All managed (deleted) by SDManager
  G4VSensitiveDetector* hodoscope1;
  G4VSensitiveDetector* chamber1;
  G4VSensitiveDetector* hodoscope2;
  G4VSensitiveDetector* chamber2;

  ConstructMaterials();

     ..... snipped .....

  // sensitive detectors -----------------------------------------------------
  G4SDManager* SDman = G4SDManager::GetSDMpointer();
  G4String SDname;

  hodoscope1 = new A01Hodoscope(SDname="/hodoscope1");
  SDman->AddNewDetector(hodoscope1);
  hodoscope1Logical->SetSensitiveDetector(hodoscope1);
  hodoscope2 = new A01Hodoscope(SDname="/hodoscope2");
  SDman->AddNewDetector(hodoscope2);
  hodoscope2Logical->SetSensitiveDetector(hodoscope2);

  chamber1 = new A01DriftChamber(SDname="/chamber1");
  SDman->AddNewDetector(chamber1);
  wirePlane1Logical->SetSensitiveDetector(chamber1);
  chamber2 = new A01DriftChamber(SDname="/chamber2");
  SDman->AddNewDetector(chamber2);
  wirePlane2Logical->SetSensitiveDetector(chamber2);

     ..... snipped .....

Then you proceed to touching to A01EventAction.cc. At first, you have to get collection ID's for hodoscope and drift chamber in the second lever arm. To get them, you have to ask them with the full collection names (detector name plus collection name) for them. Please note that these collection ID's are valid through the program unless you dynamically delete sensitive detectors.

Constructor of A01EventAction.cc 
A01EventAction::A01EventAction()
{
  G4String colName;
  G4SDManager* SDman = G4SDManager::GetSDMpointer();
  HHC1ID = SDman->GetCollectionID(colName="hodoscope1/hodoscopeColl");
  DHC1ID = SDman->GetCollectionID(colName="chamber1/driftChamberColl");
  HHC2ID = SDman->GetCollectionID(colName="hodoscope2/hodoscopeColl");
  DHC2ID = SDman->GetCollectionID(colName="chamber2/driftChamberColl");
  verboseLevel = 1;
  messenger = new A01EventActionMessenger(this);
}

To access to the hits objects, you have to follow the following procedures.

  1. Get G4HCofThisEvent object from G4Event object.
  2. Get G4VHitsCollection object from G4HCofThisEvent and cast it to the concrete class.
  3. Get number of hits stored in the collection and loop over all stored hits.

EndOfEventAction method of A01EventAction.cc 
void A01EventAction::EndOfEventAction(const G4Event* evt)
{
  G4HCofThisEvent * HCE = evt->GetHCofThisEvent();
  A01HodoscopeHitsCollection* HHC1 = 0;
  A01DriftChamberHitsCollection* DHC1 = 0;
  A01HodoscopeHitsCollection* HHC2 = 0;
  A01DriftChamberHitsCollection* DHC2 = 0;
  if(HCE)
  {
    HHC1 = (A01HodoscopeHitsCollection*)(HCE->GetHC(HHC1ID));
    DHC1 = (A01DriftChamberHitsCollection*)(HCE->GetHC(DHC1ID));
    HHC2 = (A01HodoscopeHitsCollection*)(HCE->GetHC(HHC2ID));
    DHC2 = (A01DriftChamberHitsCollection*)(HCE->GetHC(DHC2ID));
  }

  // Diagnostics

  if (verboseLevel==0 || evt->GetEventID() % verboseLevel != 0) return;

  G4PrimaryParticle* primary = evt->GetPrimaryVertex(0)->GetPrimary(0);
  G4cout << G4endl
         << ">>> Event " << evt->GetEventID() << " >>> Simulation truth : "
         << primary->GetG4code()->GetParticleName()
         << " " << primary->GetMomentum() << G4endl;

  if(HHC1)
  {
    int n_hit = HHC1->entries();
    G4cout << "Hodoscope 1 has " << n_hit << " hits." << G4endl;
    for(int i1=0;i1<n_hit;i1++)
    {
      A01HodoscopeHit* aHit = (*HHC1)[i1];
      aHit->Print();
    }
  }
  if(HHC2)
  {
    int n_hit = HHC2->entries();
    G4cout << "Hodoscope 2 has " << n_hit << " hits." << G4endl;
    for(int i1=0;i1<n_hit;i1++)
    {
      A01HodoscopeHit* aHit = (*HHC2)[i1];
      aHit->Print();
    }
  }
  if(DHC1)
  {
    int n_hit = DHC1->entries();
    G4cout << "Drift Chamber 1 has " << n_hit << " hits." << G4endl;
    for(int i2=0;i2<5;i2++)
    {
      for(int i1=0;i1<n_hit;i1++)
      {
        A01DriftChamberHit* aHit = (*DHC1)[i1];
        if(aHit->GetLayerID()==i2) aHit->Print();
      }
    }
  }
  if(DHC2)
  {
    int n_hit = DHC2->entries();
    G4cout << "Drift Chamber 2 has " << n_hit << " hits." << G4endl;
    for(int i2=0;i2<5;i2++)
    {
      for(int i1=0;i1<n_hit;i1++)
      {
        A01DriftChamberHit* aHit = (*DHC2)[i1];
        if(aHit->GetLayerID()==i2) aHit->Print();
      }
    }
  }
}

Now, complie/link and execute your code. On the display, you should get a summary of this event similar to this

An event summary 
>>> Event 0 >>> Simulation truth : e+ (-4.61596,0,989.293)
Hodoscope 1 has 1 hits.
  Hodoscope[7] 4.98684 (nsec)
Hodoscope 2 has 1 hits.
  Hodoscope[8] 43.1263 (nsec)
Drift Chamber 1 has 5 hits.
  Layer[0] : time 6.67102 (nsec) --- local (x,y) -8.96713, 0.263021
  Layer[1] : time 8.33886 (nsec) --- local (x,y) -11.1916, -0.0772413
  Layer[2] : time 10.0067 (nsec) --- local (x,y) -13.3434, -0.591037
  Layer[3] : time 11.6745 (nsec) --- local (x,y) -15.5842, -1.03141
  Layer[4] : time 13.3424 (nsec) --- local (x,y) -17.7618, -1.50829
Drift Chamber 2 has 5 hits.
  Layer[0] : time 34.7853 (nsec) --- local (x,y) -201.441, -8.35232
  Layer[1] : time 36.4567 (nsec) --- local (x,y) -234.029, -9.05564
  Layer[2] : time 38.1281 (nsec) --- local (x,y) -266.918, -9.89591
  Layer[3] : time 39.7996 (nsec) --- local (x,y) -299.961, -11.0535
  Layer[4] : time 41.4711 (nsec) --- local (x,y) -333.273, -11.9274

Now you get hits from hodoscope and drift chamber in the second lever arm.

Sensitivity of EM and hadronic calorimeter

This time, you have to start with implementing a sensitive detector class and hit class of the hadronic calorimeter. The header file and skeleton of implementation is provided. Let's look at A01HadCalorimeterHit.hh

A01HadCalorimeterHit.hh 
#ifndef A01HadCalorimeterHit_h
#define A01HadCalorimeterHit_h 1

#include "G4VHit.hh"
#include "G4THitsCollection.hh"

     ..... snipped .....

class A01HadCalorimeterHit : public G4VHit
{

     ..... snipped .....

  private:
      G4int columnID;
      G4int rowID;
      G4double edep;

     ..... snipped .....

  public:
      inline void SetColumnID(G4int z) { columnID = z; }
      inline G4int GetColumnID() const { return columnID; }
      inline void SetRowID(G4int z) { rowID = z; }
      inline G4int GetRowID() const { return rowID; }
      inline void SetEdep(G4double de) { edep = de; }
      inline void AddEdep(G4double de) { edep += de; }
      inline G4double GetEdep() const { return edep; }
};

typedef G4THitsCollection<A01HadCalorimeterHit> A01HadCalorimeterHitsCollection;

     ..... snipped .....

As you see in the code above, A01HadCalorimeterHit is derived from G4VHit and A01HadCalorimeterHitsCollection is an instantiation of G4THitsCollection template class. Taking above class into account, let's start implementing each method of A01HadCaloriemter.cc.

In the constructor, the name of hits collection must be defined.

Constructor of A01HadCalorimeterHit 
A01HadCalorimeter::A01HadCalorimeter(G4String name)
:G4VSensitiveDetector(name)
{
  G4String HCname;
  collectionName.insert(HCname="HadCalorimeterColl");
  HCID = -1;
}

Initialize() method is invoked at the beginning of each event. Here you have to instantiate a hits collection object and set it to G4HCofThisEvent object. Also, since it is a collection for caloriemter hits, hits for all calorimeter cells are instantiated here with empty energy deposition and inserted to the collection. For the case of tracker, you can instantiate empty collection without any hits. Collection ID is stable for all runs in this example. Thus, the collection ID is gotten only once.

Initialize method of A01HadCalorimeterHit 
void A01HadCalorimeter::Initialize(G4HCofThisEvent*HCE)
{
  hitsCollection = new A01HadCalorimeterHitsCollection
                   (SensitiveDetectorName,collectionName[0]);
  if(HCID<0)
  { HCID = G4SDManager::GetSDMpointer()->GetCollectionID(hitsCollection); }
  HCE->AddHitsCollection(HCID,hitsCollection);

  // fill calorimeter hits with zero energy deposition
  for(int iColumn=0;iColumn<10;iColumn++)
  for(int iRow=0;iRow<2;iRow++)
  {
    A01HadCalorimeterHit* aHit = new A01HadCalorimeterHit();
    hitsCollection->insert( aHit );
  }
}

Now, in ProcessHits() method, you should add the energy deposition along the step to the hit corresponds to the cell where the step takes place. The row and column number of the cell are taken from second and third grandparents' copy numbers. (Note that the copy number of the current volume is always 0 for this case, and the copy number of its parent volume means the layer number.)

In the if block with "// check if it is first touch" is meant for storing transformation and rotation matrix of the cell for the sake of drawing the hit.

ProcessHits method of A01HadCalorimeterHit 
G4bool A01HadCalorimeter::ProcessHits(G4Step*aStep,G4TouchableHistory* /*ROhist*/)
{
  G4double edep = aStep->GetTotalEnergyDeposit();
  if(edep==0.) return true;

  G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
  G4TouchableHistory* theTouchable
    = (G4TouchableHistory*)(preStepPoint->GetTouchable());
  G4VPhysicalVolume* theCellPhysical = theTouchable->GetVolume(2);
  G4int rowNo = theCellPhysical->GetCopyNo();
  G4VPhysicalVolume* theColumnPhysical = theTouchable->GetVolume(3);
  G4int columnNo = theColumnPhysical->GetCopyNo();
  G4int hitID = 2*columnNo+rowNo;
  A01HadCalorimeterHit* aHit = (*hitsCollection)[hitID];

  // add energy deposition
  aHit->AddEdep(edep);

  // check if it is first touch
  if(aHit->GetColumnID()<0)
  {
    aHit->SetColumnID(columnNo);
    aHit->SetRowID(rowNo);
    G4int depth = theTouchable->GetHistory()->GetDepth();
    G4AffineTransform aTrans = theTouchable->GetHistory()->GetTransform(depth-2);
    aTrans.Invert();
    aHit->SetRot(aTrans.NetRotation());
    aHit->SetPos(aTrans.NetTranslation());
  }

  return true;
}

Now you finalized implementing A01HadCalorimeter. Nest is using it and A01EmCalorimeter (which is provided) in A01DetectorConstruction. This implementation is exactly the same way as for hodoscope and drift chamber.

A part of A01DetectorConstruction.cc 
#include "A01EmCalorimeter.hh"
#include "A01HadCalorimeter.hh"

G4VPhysicalVolume* A01DetectorConstruction::Construct()
{
  // All managed (deleted) by SDManager
  G4VSensitiveDetector* hodoscope1;
  G4VSensitiveDetector* chamber1;
  G4VSensitiveDetector* hodoscope2;
  G4VSensitiveDetector* chamber2;
  G4VSensitiveDetector* EMcalorimeter;
  G4VSensitiveDetector* HadCalorimeter;

  ConstructMaterials();

     ..... snipped .....

  // sensitive detectors -----------------------------------------------------
  G4SDManager* SDman = G4SDManager::GetSDMpointer();
  G4String SDname;

  hodoscope1 = new A01Hodoscope(SDname="/hodoscope1");
  SDman->AddNewDetector(hodoscope1);
  hodoscope1Logical->SetSensitiveDetector(hodoscope1);
  hodoscope2 = new A01Hodoscope(SDname="/hodoscope2");
  SDman->AddNewDetector(hodoscope2);
  hodoscope2Logical->SetSensitiveDetector(hodoscope2);

  chamber1 = new A01DriftChamber(SDname="/chamber1");
  SDman->AddNewDetector(chamber1);
  wirePlane1Logical->SetSensitiveDetector(chamber1);
  chamber2 = new A01DriftChamber(SDname="/chamber2");
  SDman->AddNewDetector(chamber2);
  wirePlane2Logical->SetSensitiveDetector(chamber2);

  EMcalorimeter = new A01EmCalorimeter(SDname="/EMcalorimeter");
  SDman->AddNewDetector(EMcalorimeter);
  cellLogical->SetSensitiveDetector(EMcalorimeter);

  HadCalorimeter = new A01HadCalorimeter(SDname="/HadCalorimeter");
  SDman->AddNewDetector(HadCalorimeter);
  HadCalScintiLogical->SetSensitiveDetector(HadCalorimeter);

     ..... snipped .....

Finally, touch to A01EventAction for accessing to corresponding hits.

A01EventAction.cc 
#include "A01EmCalorimeterHit.hh"
#include "A01HadCalorimeterHit.hh"

A01EventAction::A01EventAction()
{
  G4String colName;
  G4SDManager* SDman = G4SDManager::GetSDMpointer();
  HHC1ID = SDman->GetCollectionID(colName="hodoscope1/hodoscopeColl");
  DHC1ID = SDman->GetCollectionID(colName="chamber1/driftChamberColl");
  HHC2ID = SDman->GetCollectionID(colName="hodoscope2/hodoscopeColl");
  DHC2ID = SDman->GetCollectionID(colName="chamber2/driftChamberColl");
  ECHCID = SDman->GetCollectionID(colName="EMcalorimeter/EMcalorimeterColl");
  HCHCID = SDman->GetCollectionID(colName="HadCalorimeter/HadCalorimeterColl");

  verboseLevel = 1;
  messenger = new A01EventActionMessenger(this);

}

     ..... snipped .....

void A01EventAction::EndOfEventAction(const G4Event* evt)
{
  G4HCofThisEvent * HCE = evt->GetHCofThisEvent();
  A01HodoscopeHitsCollection* HHC1 = 0;
  A01DriftChamberHitsCollection* DHC1 = 0;
  A01HodoscopeHitsCollection* HHC2 = 0;
  A01DriftChamberHitsCollection* DHC2 = 0;
  A01EmCalorimeterHitsCollection* ECHC = 0;
  A01HadCalorimeterHitsCollection* HCHC = 0;
  if(HCE)
  {
    HHC1 = (A01HodoscopeHitsCollection*)(HCE->GetHC(HHC1ID));
    DHC1 = (A01DriftChamberHitsCollection*)(HCE->GetHC(DHC1ID));
    HHC2 = (A01HodoscopeHitsCollection*)(HCE->GetHC(HHC2ID));
    DHC2 = (A01DriftChamberHitsCollection*)(HCE->GetHC(DHC2ID));
    ECHC = (A01EmCalorimeterHitsCollection*)(HCE->GetHC(ECHCID));
    HCHC = (A01HadCalorimeterHitsCollection*)(HCE->GetHC(HCHCID));
  }

      ..... snipped .....

  if(DHC2)
  {
    int n_hit = DHC2->entries();
    G4cout << "Drift Chamber 2 has " << n_hit << " hits." << G4endl;
    for(int i2=0;i2<5;i2++)
    {
      for(int i1=0;i1<n_hit;i1++)
      {
        A01DriftChamberHit* aHit = (*DHC2)[i1];
        if(aHit->GetLayerID()==i2) aHit->Print();
      }
    }
  }
  if(ECHC)
  {
    int iHit = 0;
    double totalE = 0.;
    for(int i1=0;i1<80;i1++)
    {
      A01EmCalorimeterHit* aHit = (*ECHC)[i1];
      double eDep = aHit->GetEdep();
      if(eDep>0.)
      {
        iHit++;
        totalE += eDep;
      }
    }
    G4cout << "EM Calorimeter has " << iHit << " hits. Total Edep is "
           << totalE/MeV << " (MeV)" << G4endl;
  }
  if(HCHC)
  {
    int iHit = 0;
    double totalE = 0.;
    for(int i1=0;i1<20;i1++)
    {
      A01HadCalorimeterHit* aHit = (*HCHC)[i1];
      double eDep = aHit->GetEdep();
      if(eDep>0.)
      {
        iHit++;
        totalE += eDep;
      }
    }
    G4cout << "Hadron Calorimeter has " << iHit << " hits. Total Edep is "
           << totalE/MeV << " (MeV)" << G4endl;
  }
}

Now, you can try shooting hadron to see hits in EM and hadronic calorimeters. After compiling and linking, execute the program. After some initialization, a prompt Idle> should appear. Then, adopt a UI command for proper visualization, and then execute an event. 

  Idle> /gun/particle proton
  Idle> /mydet/momentum 10 GeV
  Idle> /mydet/fieldValue 10 tesla
  Idle> /run/beamOn 1
You should get these event summary and a drawing.

An event summary 
>>> Event 0 >>> Simulation truth : proton (-46.6088,0,9989.19)
Hodoscope 1 has 1 hits.
  Hodoscope[7] 5.00879 (nsec)
Hodoscope 2 has 1 hits.
  Hodoscope[9] 43.311 (nsec)
Drift Chamber 1 has 5 hits.
  Layer[0] : time 6.70039 (nsec) --- local (x,y) -9.33916, 0.035515
  Layer[1] : time 8.37558 (nsec) --- local (x,y) -11.8348, 0.0592039
  Layer[2] : time 10.0508 (nsec) --- local (x,y) -14.2634, 0.0896775
  Layer[3] : time 11.7259 (nsec) --- local (x,y) -16.6732, 0.107581
  Layer[4] : time 13.4011 (nsec) --- local (x,y) -19.0977, 0.157995
Drift Chamber 2 has 5 hits.
  Layer[0] : time 34.9371 (nsec) --- local (x,y) -187.033, 0.796012
  Layer[1] : time 36.6152 (nsec) --- local (x,y) -216.479, 0.856811
  Layer[2] : time 38.2932 (nsec) --- local (x,y) -245.938, 0.902882
  Layer[3] : time 39.9713 (nsec) --- local (x,y) -275.416, 0.917685
  Layer[4] : time 41.6494 (nsec) --- local (x,y) -304.883, 0.913767
EM Calorimeter has 18 hits. Total Edep is 2489.86 (MeV)
Hadron Calorimeter has 6 hits. Total Edep is 129.766 (MeV)

Tutorial by:
Makoto Asai

11 October 2003