LCDG4CalSD.cc

Go to the documentation of this file.

#include "LCDG4CalSD.hh"

// Geant4
// #include "G4String.hh"
// #include "G4Step.hh"
// #include "G4HCofThisEvent.hh"
// #include "G4Track.hh"

#include <map>
using std::map;
#include <fstream>
using std::ifstream;

/*
 * Constructor
 */
LCDG4CalSD::LCDG4CalSD(G4String name)
  : G4VSensitiveDetector(G4String("/lcddet/")+name)
  , HCID(-1)
  , layer0CenterDepth(0)
  , layerThickness(0)
  , debugLevel(0)
{
  // read debugLevel flag from file ./debfile.
  // If no file is present, debugLevel = 0
  // The higher its value, the more output
  ifstream* debfile = new ifstream("debfile");
  if(debfile) {
    *debfile >> debugLevel;
    debfile->close();
    delete debfile;
    debfile = 0;
  }

  // save name of hits collection associated to this Sensitive Detector
  G4String HCname = name+"Collection";
  collectionName.insert( HCname );

  G4cout<< "SensitiveDetectorName = "<< SensitiveDetectorName << G4endl;
}

/*
 * Destructor
 */
LCDG4CalSD::~LCDG4CalSD() {
  // methinks this causes Seg Fault JM062503
  // delete _hitCollection;
}

/*
 * Initialize the hit collection for each event
 */
void LCDG4CalSD::Initialize(G4HCofThisEvent* HCE) {
  // Create a new hit collection for this event
  _hitCollection = 
    new LCDG4CalHitsCollection(SensitiveDetectorName,
                               collectionName[0]);
  // save hit collection id
  if(HCID < 0) HCID = GetCollectionID(0);
  // add hit collection to event
  HCE->AddHitsCollection(HCID,_hitCollection);
}

/*
 * Process the step information, feed it into LCDG4CalHit
 */
G4bool LCDG4CalSD::ProcessHits(G4Step* aStep, G4TouchableHistory* ROhist) {

  G4double edep = aStep->GetTotalEnergyDeposit();
  if(edep==0) return false;

  const G4VPhysicalVolume* physVol
    = aStep->GetPreStepPoint()->GetPhysicalVolume();
  if( ROhist!=NULL ) {
    G4VPhysicalVolume* ROvol = ROhist->GetVolume();
    if(ROvol!=physVol) {
      G4cout<<"ROvol="<< ROvol <<", physVol="<< physVol <<G4endl;
    }
  }

  // some aliases
  unsigned int& sysNo = index[0];
  unsigned int& layno = index[1];

  G4LogicalVolume* logVol=physVol->GetLogicalVolume();
#ifndef TXT_MODE
  // no need to keep track of hits in absorbers
  if (!_mapsens[logVol]) return false;
#endif

  sysNo = _mapsys[logVol];
  layno = _maplayer[logVol];
  barend = _mapbe[logVol];
  north = _mapns[logVol];

  G4Track* ptrk = aStep->GetTrack();
  G4int trkID = ptrk->GetTrackID();
  G4double tdep = ptrk->GetGlobalTime();

  // edep hit position is at middle of preStep and postStep points
  G4ThreeVector pos = 0.5*(aStep->GetPreStepPoint()->GetPosition()+
                           aStep->GetPostStepPoint()->GetPosition() );

  if(debugLevel) {
    G4cout << "\nCheck PhysVol on pre/post step: "
           <<aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName() <<' '
           <<aStep->GetPostStepPoint()->GetPhysicalVolume()->GetName()<<G4endl;
    G4cout << "Hit in layer " << layno << G4endl
           << "hit info (x,y,z)=("<< pos.x() <<';'<< pos.y() <<';'<< pos.z()
           <<"), cylRho="<< pos.rho() <<", th="<< pos.theta() <<", phi="<< pos.phi()
           << G4endl;
  }

  // Geometry dependent calculations
  unsigned int cellIndex = findCell(pos);
  G4ThreeVector cellCenter = getCellCenter(cellIndex);
  if(debugLevel) {
    G4cout<< "cellIndex="<< cellIndex
          <<'/'<< getIndex3(cellIndex)
          <<'/'<< getIndex2(cellIndex)
          <<'/'<< getIndex1(cellIndex) <<'/'<< sysNo
           << " -> cell center info (x,y,z)=("<< cellCenter.x() <<';'<< cellCenter.y() <<';'<< cellCenter.z()
          <<"), cylRho="<< cellCenter.rho() <<", th="<< cellCenter.theta() <<", phi="<< cellCenter.phi() <<' '<<", north="<< north << G4endl;
  }

  // Get a pointer to hit if it exists, or create a new one
  LCDG4CalHit *hit = NULL;
  map<int,LCDG4CalHit*>::iterator pcellhit = _cellhit.find(cellIndex);
  if( pcellhit!=_cellhit.end() ) hit = _cellhit.find(cellIndex)->second;
  if(hit == NULL) {
    if (debugLevel>3) G4cout << "creating new hit cell" << G4endl;
    hit = new LCDG4CalHit();
    hit->SetSystem(sysNo);
    hit->SetBarEnd(barend);
    hit->SetNorth(north);
    hit->SetLayer(layno);

    // set hit position to center of cell
    hit->SetPos(cellCenter);

    // save cell indices (virtual cells)
    hit->SetPhi( getIndex1(cellIndex) );
    hit->SetTheta( getIndex2(cellIndex) );

    // Insert hit into hit collection
    _hitCollection->insert(hit);
    _cellhit.insert( std::pair<G4int,LCDG4CalHit*>( cellIndex, hit) );
  }

  if (_mapsens[logVol]) {
    // hit in sensitive material
    hit->AddMcPart(trkID,edep,tdep);
  } else {
    // hit in absorber
    hit->AddMcPartAbs(edep);
  }
  if(debugLevel>3) {
    G4cout<<"After adding hit: hit(edep,edepa) = "
          << hit->GetEdep() <<' '<< hit->GetEdepAbs() << G4endl;
  }
  //    if ( hit->GetNMC() > 1 ) {
  //      G4cout <<"CalSD: cellid,nmc = "<< layno <<' '<< sphi <<' '<< sthe <<' '<< hit->GetNMC() << G4endl;
  //      for(int i=0; i<hit->GetNMC(); i++) {
  //        G4cout <<"  hit,id,emc = " <<' '<< i <<' '<< hit->GetMCIdx(i) <<' '<< hit->GetEMC(i) << G4endl;
  //      }
  //    }

  return true;
}


/*
 * EndOfEvent method clear cell hit maps for barrel and endcap
 */
void LCDG4CalSD::EndOfEvent(G4HCofThisEvent* HCE) {
  if(HCE) ;  // make compiler happy
  // clear member hit collections.  The hits are owned by the event, 
  // so deleting them here will produce a crash
  _cellhit.clear();
}

/*
 * Add layer no, sysNo, barrel/endcap, N/S, and sensitive information
 * to logical volume maps
 */
G4int LCDG4CalSD::AddLogVolInfo(G4LogicalVolume* g4logvol,
                                G4int startlayerNo,
                                G4int sysNo, G4int be,
                                G4int ns, G4int sens)
{
  _maplayer.insert(std::pair<G4LogicalVolume*,G4int>(g4logvol,startlayerNo));
  _mapsys.insert(std::pair<G4LogicalVolume*,G4int>(g4logvol,sysNo));
  _mapbe.insert(std::pair<G4LogicalVolume*,G4int>(g4logvol,be));
  _mapns.insert(std::pair<G4LogicalVolume*,G4int>(g4logvol,ns));
  _mapsens.insert(std::pair<G4LogicalVolume*,G4int>(g4logvol,sens));
  return _maplayer.size();
}

/*
 * Clear method, not used
 */
void LCDG4CalSD::clear() {
}

/*
 * Draw method, not used
 */
void LCDG4CalSD::DrawAll() {
}

/*
 * Print method, not used
 */
void LCDG4CalSD::PrintAll() {
}

---