TimeProjectionChamberSensitiveDetector.cpp

#include "TimeProjectionChamberSensitiveDetector.h"

#include "G4Step.hh"
#include "G4VProcess.hh"
//#include "G4SteppingManager.hh"
#include "G4SDManager.hh"
//#include "UserTrackInformation.hh"
#include "DD4hep/DD4hepUnits.h"

TimeProjectionChamberSensitiveDetector::TimeProjectionChamberSensitiveDetector(const std::string& name,
									       dd4hep::Detector& description)
  : DDG4SensitiveDetector(name, description){
  
  G4String CollName1=name+"Collection";
  collectionName.insert(CollName1);
  
  G4String CollName2=name+"SpacePointCollection";
  collectionName.insert(CollName2);
  
  G4String CollName3=name+"LowPtCollection";
  collectionName.insert(CollName3);
}

void TimeProjectionChamberSensitiveDetector::Initialize(G4HCofThisEvent* HCE){
  m_hc = new HitCollection(GetName(), collectionName[0]);
  int HCID = G4SDManager::GetSDMpointer()->GetCollectionID(m_hc);
  HCE->AddHitsCollection( HCID, m_hc ); 

  m_spaceHC = new HitCollection(GetName(), collectionName[1]);
  int HCIDSpace = G4SDManager::GetSDMpointer()->GetCollectionID(m_spaceHC);
  HCE->AddHitsCollection( HCIDSpace, m_spaceHC );

  m_lowPtHC = new HitCollection(GetName(), collectionName[2]);
  int HCIDLowPt = G4SDManager::GetSDMpointer()->GetCollectionID(m_lowPtHC);
  HCE->AddHitsCollection( HCIDLowPt, m_lowPtHC );
  
  CumulativeNumSteps=0;
  
  dEInPadRow = 0.0;
  CumulativeEnergyDeposit = 0.0;
  globalTimeAtPadRingCentre=0.0;
  pathLengthInPadRow=0.0;
  CumulativePathLength=0.0;
  CrossingOfPadRingCentre.SetXYZ(0,0,0);// = dd4hep::Position(0,0,0);
  MomentumAtPadRingCentre.SetXYZ(0,0,0);// = dd4hep::sim::Momentum(0,0,0);
  
  CumulativeMeanPosition.SetXYZ(0,0,0);// = dd4hep::Position(0,0,0);
  CumulativeMeanMomentum.SetXYZ(0,0,0);// = dd4hep::sim::Momentum(0,0,0);
  
  PreviousPostStepPosition.SetXYZ(0,0,0);// = dd4hep::Position(0,0,0);
  CurrentPDGEncoding=0;
  CurrentTrackID=-1;
  CurrentGlobalTime=0.0;
  CurrentCopyNumber=0;
  
  /*
  if( m_lowPtStepLimit ) {
    std::cout << "TimeProjectionChamberSensitiveDetector: TPCLowPtStepLimit ON" << std::endl;
  }
  else {
    std::cout << "TimeProjectionChamberSensitiveDetector: TPCLowPtStepLimit OFF" << std::endl;
  }

  if( m_writeMCParticleForLowPtHits ) {
    std::cout << "TimeProjectionChamberSensitiveDetector: TPCLowPtStoreMCPForHits ON" << std::endl;
  }
  else {
    std::cout << "TimeProjectionChamberSensitiveDetector: TPCLowPtStoreMCPForHits OFF" << std::endl;
  }

  std::cout << "TimeProjectionChamberSensitiveDetector: Threshold for Energy Deposit = " << m_thresholdEnergyDeposit / CLHEP::eV << " eV" << G4endl;
  std::cout << "TimeProjectionChamberSensitiveDetector: TPCLowPtCut = " << m_lowPtCut / CLHEP::MeV << " MeV "<< G4endl;
  std::cout << "TimeProjectionChamberSensitiveDetector: TPCLowPt Max hit separation "<< m_lowPtMaxHitSeparation / CLHEP::mm << " mm" << G4endl;
  */
}

G4bool TimeProjectionChamberSensitiveDetector::ProcessHits(G4Step* step, G4TouchableHistory*){
  // FIXME: 
  // (i) in the following algorithm if a particle "appears" within a pad-ring half and 
  // leaves without passing through the middle of the pad-ring it will not create a hit in 
  // this ring
  // (ii) a particle that crosses the boundry between two pad-ring halves will have the hit 
  // placed on this surface at the last crossing point, and will be assinged the total energy 
  // deposited in the whole pad-ring. This is a possible source of bias for the hit
  
  G4TouchableHandle touchPost = step->GetPostStepPoint()->GetTouchableHandle(); 
  G4TouchableHandle touchPre  = step->GetPreStepPoint()->GetTouchableHandle(); 
  dd4hep::sim::Geant4StepHandler h(step);
  
  if (fabs(h.trackDef()->GetPDGCharge()) < 0.01) return true;
  
  const dd4hep::Position PrePosition   = h.prePos();
  const dd4hep::Position PostPosition  = h.postPos();
  const dd4hep::sim::Momentum Momentum = h.postMom(); //use post step momentum, Mokka history
  
  float ptSQRD = Momentum.X()*Momentum.X()+Momentum.Y()*Momentum.Y();
  if( ptSQRD >= (m_lowPtCut*m_lowPtCut) ){
    // Step finishes at a geometric boundry
    if(step->GetPostStepPoint()->GetStepStatus() == fGeomBoundary){
    //if(h.postStepStatus() == "GeomBoundary"){
      // step within the same pair of upper and lower pad ring halves
      if(int(touchPre->GetCopyNumber()/2)==int(touchPost->GetCopyNumber()/2)) {
        //this step must have ended on the boundry between these two pad ring halfs 
        //record the tracks coordinates at this position 
        //and return

        CrossingOfPadRingCentre.SetXYZ(PostPosition.X(),PostPosition.Y(),PostPosition.Z());
        MomentumAtPadRingCentre   = Momentum;  
        dEInPadRow               += h.deposit();
        globalTimeAtPadRingCentre = h.trkTime();
        pathLengthInPadRow       += h.stepLength();
        
        //	    G4cout << "step must have ended on the boundry between these two pad ring halfs" << G4endl;
        //	    G4cout << "CrossingOfPadRingCentre = "   
        //		   << sqrt( CrossingOfPadRingCentre[0]*CrossingOfPadRingCentre[0]
        //			    +CrossingOfPadRingCentre[1]*CrossingOfPadRingCentre[1] )
        //		   << G4endl;
        
        return true;
      }
      else if(!(CrossingOfPadRingCentre.X()==0.0 && CrossingOfPadRingCentre.Y()==0.0 && CrossingOfPadRingCentre.Z()==0.0)) {
        // the above IF statment is to catch the case where the particle "appears" in this pad-row half volume and 
        // leaves with out crossing the pad-ring centre, as mentioned above
        
        //it is leaving the pad ring couplet
        //write out a hit
        //make sure particle is added to MC list
        //and return
        //	    G4cout << "step must be leaving the pad ring couplet" << G4endl;
        //	    G4cout << "write out hit at " 
        //		   << sqrt( CrossingOfPadRingCentre[0]*CrossingOfPadRingCentre[0]
        //			    +CrossingOfPadRingCentre[1]*CrossingOfPadRingCentre[1] )
        //		   << " " << "dEdx = " << step->GetTotalEnergyDeposit()+dEInPadRow 
        //		   << " " << "step length = " << step->GetStepLength()+pathLengthInPadRow  
        //		   << G4endl;
        
        double dE = h.deposit()+dEInPadRow;
	if ( dE > m_thresholdEnergyDeposit || m_trackingPhysicsListELossOn == false ) {
	  // needed for delta electrons: all particles causing hits have to be saved
	  // TODO: set track flag to save

	  dd4hep::sim::Geant4TrackerHit* hit = new dd4hep::sim::Geant4TrackerHit(h.trkID(),h.trkPdgID(),dE,globalTimeAtPadRingCentre);
	  hit->cellID        = touchPre->GetCopyNumber()/2+1;//getCellID(step);
	  hit->energyDeposit = dE;
	  hit->position      = CrossingOfPadRingCentre;
	  hit->momentum      = MomentumAtPadRingCentre;
	  hit->length        = h.stepLength()+pathLengthInPadRow;
          m_hc->insert(hit);
        }
        
        // zero cumulative variables 
        dEInPadRow                = 0.0;
        globalTimeAtPadRingCentre = 0.0;
        pathLengthInPadRow        = 0.0;
        CrossingOfPadRingCentre.SetXYZ(0.,0.,0.);
        MomentumAtPadRingCentre.SetXYZ(0.,0.,0.);
        return true;
      }
    }
    //case for which the step remains within geometric volume
    //FIXME: need and another IF case to catch particles which Stop within the padring
    else if(step->GetPostStepPoint()->GetStepStatus() != fGeomBoundary){
    //else if(h.postStepStatus() != "GeomBoundary") {
      //the step is not limited by the step length 
      if(step->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName()!="StepLimiter"){
        // if(particle not stoped){ 	    
        //add the dEdx and return 
        //	    G4cout << "Step ended by Physics Process: Add dEdx and carry on" << G4endl;
        dEInPadRow         += h.deposit();
        pathLengthInPadRow += h.stepLength();
	return true;
        //}
        //else{
        //  write out the hit and clear counters
        //}
      }
      else {
	//double position_x = PostPosition[0];
        //double position_y = PostPosition[1];
        //double position_z = PostPosition[2];
        
        const dd4hep::sim::Momentum PreMomentum = h.preMom();
        const dd4hep::sim::Momentum PostMomentum = h.postMom();
        
        //double momentum_x = PostMomentum[0];
        //double momentum_y = PostMomentum[1];
        //double momentum_z = PostMomentum[2];
        
        
        //	    G4cout << "step must have been stopped by the step limiter" << G4endl;
        //	    G4cout << "write out hit at " 
        //		   << sqrt( position_x*position_x
        //			    +position_y*position_y )
        //		   << " " << "dEdx = " << step->GetTotalEnergyDeposit() 
        //		   << " " << "step length = " << step->GetStepLength()  
        //		   << G4endl;
        
        // write out step limited hit 
        // these are just space point hits so do not save the dE, which is set to ZERO
        //	    if ( step->GetTotalEnergyDeposit() > fThresholdEnergyDeposit ) 

        // needed for delta electrons: all particles causing hits have to be saved in the LCIO file
	// TODO: set track flag to be saved
        //UserTrackInformation *info = (UserTrackInformation*)(step->GetTrack()->GetUserInformation());
        //if (info) info->GetTheTrackSummary()->SetToBeSaved();

	dd4hep::sim::Geant4TrackerHit* hit = new dd4hep::sim::Geant4TrackerHit(h.trkID(),h.trkPdgID(),0.0,h.trkTime());
	hit->cellID        = touchPre->GetCopyNumber()/2+1;//getCellID(step);
	hit->energyDeposit = 0.0;
	hit->position      = PostPosition;
	hit->momentum      = PostMomentum;
	hit->length        = h.stepLength();
	m_spaceHC->insert(hit);
	
        // add dE and pathlegth and return
        dEInPadRow += h.deposit();
        pathLengthInPadRow += h.stepLength();
        return true;
      }
    }
  }
  else if(m_lowPtStepLimit) { // low pt tracks will be treated differently as their step length is limited by the special low pt steplimiter
    if ( ( PreviousPostStepPosition - PrePosition ).R() > 1.0e-6 * dd4hep::mm ) {
      // This step does not continue the previous path. Deposit the energy and begin a new Pt hit.
      if (CumulativeEnergyDeposit > m_thresholdEnergyDeposit) {
        DepositLowPtHit();
      }
      else {
        // reset the cumulative variables if the hit has not been deposited.
        // The previous track has ended and the cumulated energy left at the end 
        // was not enough to ionize
        //G4cout << "reset due to new track , discarding " << CumulativeEnergyDeposit / eV << " eV" << std::endl;
        ResetCumulativeVariables();
      }
    }
    else {
      //G4cout << "continuing track" << endl;
    }
    
    CumulateLowPtStep(h);  

    // check whether to deposit the hit
    if( ( CumulativePathLength > m_lowPtMaxHitSeparation )  ) {
      // hit is deposited because the step limit is reached and there is enough energy
      // to ionize
      if ( CumulativeEnergyDeposit > m_thresholdEnergyDeposit) {
        DepositLowPtHit();
      }
      //else {
      //G4cout << "not deposited, energy is " << CumulativeEnergyDeposit/eV << " eV" << std::endl;
      //}
    }
    else { // even if the step lenth has not been reached the hit might
           // be deposited because the particle track ends
      if ( h.post->GetKineticEnergy() == 0 ) {
	// only deposit the hit if the energy is high enough
        if (CumulativeEnergyDeposit > m_thresholdEnergyDeposit) {
	  DepositLowPtHit();
        }
        else { // energy is not enoug to ionize.
               // However, the track has ended and the energy is discarded and not added to the next step
               //G4cout << "reset due to end of track, discarding " << CumulativeEnergyDeposit/eV << " eV" << std::endl;
          ResetCumulativeVariables();
        }
      }
    }
    
    PreviousPostStepPosition = h.postPos();
    
    return true;
  }

  return true;
}

void TimeProjectionChamberSensitiveDetector::EndOfEvent(G4HCofThisEvent* HCE){
  // There might be one low Pt hit which has not been added to the collection.
  if (CumulativeEnergyDeposit > m_thresholdEnergyDeposit) {
    DepositLowPtHit(); 
  }  
}

void TimeProjectionChamberSensitiveDetector::DepositLowPtHit(){
  dd4hep::sim::Geant4TrackerHit* hit = new dd4hep::sim::Geant4TrackerHit(CurrentTrackID,CurrentPDGEncoding,CumulativeEnergyDeposit,CurrentGlobalTime);
  hit->cellID        = CurrentCopyNumber;
  hit->energyDeposit = CumulativeEnergyDeposit; // FIXME: also use the dedx
  hit->position      = CumulativeMeanPosition;
  hit->momentum      = CumulativeMeanMomentum;
  hit->length        = CumulativePathLength;
  m_lowPtHC->insert(hit);
  // reset the cumulative variables after positioning the hit
  ResetCumulativeVariables();
}

void TimeProjectionChamberSensitiveDetector::ResetCumulativeVariables(){
  CumulativeMeanPosition.SetXYZ(0.,0.,0.);// = dd4hep::Position(0.0,0.0,0.0);
  CumulativeMeanMomentum.SetXYZ(0.,0.,0.);// = dd4hep::sim::Momentum(0.0,0.0,0.0);
  CumulativeNumSteps = 0;
  CumulativeEnergyDeposit = 0;
  CumulativePathLength = 0;
}

void TimeProjectionChamberSensitiveDetector::CumulateLowPtStep(dd4hep::sim::Geant4StepHandler& h){
  //dd4hep::Position prePos    = h.prePos();
  //dd4hep::Position postPos   = h.postPos();
  //dd4hep::Position direction = postPos - prePos;
  dd4hep::Position meanPosition  = h.avgPosition();//0.5*(prePos+postPos);
  //double           hit_len   = direction.R();
  dd4hep::Position meanMomentum = 0.5*(h.preMom() + h.postMom());
  
  ++CumulativeNumSteps;    
  CumulativeMeanPosition   = ( (CumulativeMeanPosition*(CumulativeNumSteps-1)) + meanPosition ) / CumulativeNumSteps;
  CumulativeMeanMomentum   = ( (CumulativeMeanMomentum*(CumulativeNumSteps-1)) + meanMomentum ) / CumulativeNumSteps;
  CumulativeEnergyDeposit += h.deposit();
  CumulativePathLength    += h.stepLength();
  CurrentPDGEncoding       = h.trkPdgID();
  CurrentTrackID           = h.trkID();
  CurrentGlobalTime        = h.trkTime();
  CurrentCopyNumber        = h.preTouchable()->GetCopyNumber()/2+1;//getCellID( h.step );
  
  // needed for delta electrons: all particles causing hits have to be saved in the LCIO file 
  if ( CumulativeEnergyDeposit > m_thresholdEnergyDeposit ) {
    // This low Pt hit will eventually be saved, so set the flag to store the particle
    
    // writing the MC Particles can be turned on or off for the lowPt particles
    if(m_writeMCParticleForLowPtHits){
      // TODO: here set write flag for this track 

    }
  }
}