musrsim/geant4/LEMuSR/src/TDCheck.cc

//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$//*
//             LOW ENERGY MUON SPIN RELAXATION, ROTATION, RADIATION 
// 
//  ID    :TDCheck.cc , v 1.2
//  AUTHOR: Taofiq PARAISO
//  DATE  : 2005-03-01 10:07
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$//
//
//  &           &&&&&&&&&&                            &&&&&&&         &&&&&&&&
//  &           &                                   &&      &&       &       &&
//  &           &                   &      &       &                &       &&
//  &           &&&&&&&            &      &         &&&&&&         &&&&&&&&  
//  &           &                 &      &&                &      &      &&
//  &           &                &&     &  &     &&      &&      &        &
//  &&&&&&&&&&  &&&&&&&&&&      &  &&&&&    &&    &&&&&&&       &        &&   
//                             &
//                            &
//                           &
//                         &  
//                               TDCHECK   
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$//
#include "TDCheck.hh"

#include "G4SteppingManager.hh"
#include "G4Transform3D.hh"
#include "G4DynamicParticle.hh"
#include "G4UnitsTable.hh"


#include "LEMuSRDetectorConstruction.hh"
#include "LEMuSRPrimaryGeneratorAction.hh"


TDCheck::TDCheck()
{

  BookRoot();
  muon.thickness=0;
  pointer=this  ;
  loop=0;
  oldz=0;
  thk_old=0;
  id=0;
  old_id=0;
  loop=0;
}


TDCheck::~TDCheck()
{
  WriteRoot();
}

TDCheck*  TDCheck::pointer=0;
TDCheck*  TDCheck::GetInstance()
{
  return pointer;
}


void TDCheck::UserSteppingAction(const G4Step* aStep)
  
{ 
  
  if( aStep->GetPreStepPoint()&& aStep->GetPreStepPoint()->GetPhysicalVolume() )
    {
      SetParticleVolumeNames(aStep);

      if(CheckCondition(aStep))
	{	  
	  // Get datas
	  SetPositionMomentum(aStep);
	  SetTimeEnergy(aStep);
	  Update();
	  FillRoot();
#if defined G4UI_USE_ROOT 
	  myFile->Write();
#endif
	}

 
      LoopKiller(aStep);
     
    }

}


G4bool TDCheck::CheckCondition(const G4Step* )
{
  G4bool condition=false;

  if(v_name=="lv_CFOIL" )
    {	  
      if(p_name=="Mu" || p_name == "mu+")
	{
	  condition=true;
	}
    }

  if(v_name=="lv_SAH2" )
    {	  
      if(p_name=="Mu" || p_name == "mu+")
	{
	  condition=true;
	}
    }

  return condition;

}  


void TDCheck::LoopKiller(const G4Step *aStep)
{
  // loop killa
  if(aStep->GetStepLength()<0.1*mm)
    {
      loop++;
      if(loop>20)
	{
	  aStep->GetTrack()->SetTrackStatus(fStopAndKill);
	  loop=0;
	}
    }

  // kill useless particles
  if(aStep->GetTrack()->GetDefinition()->GetParticleName()=="e-"
     ||aStep->GetTrack()->GetDefinition()->GetParticleName()=="gamma")
    {
      //      if(aStep->GetTrack()->GetKineticEnergy()<10.*eV)      
      //	{
	  aStep->GetTrack()->SetTrackStatus(fStopAndKill);
	  //	}
    }
  
}




void TDCheck::SetParticleVolumeNames(const G4Step* aStep)
{
  // NAMES
  p_name = aStep->GetTrack()->GetDefinition()->GetParticleName(); // particle name 
  v_name = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName();
  pv_name = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName();
}


void TDCheck::SetPositionMomentum(const G4Step* aStep)
{
  // POSITION, MOMENTUM
  position = aStep->GetPreStepPoint()->GetPosition(); // position
  momentum = aStep->GetPostStepPoint()->GetMomentumDirection(); // momentum
  momentum_direction = aStep->GetPreStepPoint()->GetMomentumDirection(); // momentum

}


void TDCheck::SetTimeEnergy(const G4Step* aStep)
{
 
  // ENERGY
  kenergy= aStep->GetTrack()->GetDynamicParticle()->GetKineticEnergy(); // position
  
  tenergy= aStep->GetTrack()->GetDynamicParticle()->GetTotalEnergy(); // position
 
  // TIME
  localtime =  (aStep->GetPreStepPoint()->GetLocalTime()+ aStep->GetPreStepPoint()->GetLocalTime())/2.; // time since track creation
  globaltime =  aStep->GetPreStepPoint()->GetGlobalTime();// time since event creation
  proptime =  aStep->GetPreStepPoint()->GetProperTime(); // proper time of the particle
  time = proptime;
  
}

void TDCheck::Update()
{
  muon.localtime = localtime/ns ;
  muon.globaltime= globaltime/ns ;
  muon.proptime   = proptime/ns ;
  
  muon.positionx = LEMuSRPrimaryGeneratorAction::GetPGA()->X/cm; // position.x()/cm;
  muon.positiony =  LEMuSRPrimaryGeneratorAction::GetPGA()->Y/cm; //position.y()/cm;
  muon.positionz = position.z()/cm;
  
  muon.momdirx= LEMuSRPrimaryGeneratorAction::GetPGA()->angle/rad;//momentum_direction.x();
  muon.momdiry= LEMuSRPrimaryGeneratorAction::GetPGA()->radius/cm; //momentum_direction.y();
  muon.momdirz= momentum_direction.z();

  muon.tenergy   = kenergy/keV;
  muon.thickness =  LEMuSRDetectorConstruction::GetInstance()->cfthk;
  
  muon.angle=0;


      muon.index=0;
      muon.ienergy=  LEMuSRPrimaryGeneratorAction::GetPGA()->energy/keV;

      muon.angle=acos(momentum.z())*180/M_PI;
      
      /////////////////
      /*
	G4cout << "particle name     : " << p_name   <<" ;\n " 
	<< "position          : " << position        <<" ;\n " 
	<< "momentum direction: " <<momentum_direction <<" ;\n " 
	<< "time              : " <<  globaltime/ns <<" ;\n " 
	<<std::endl;
      */
      //////////////////   
   
  if( G4RunManager::GetRunManager()->GetCurrentRun()->GetRunID()!=old_id)
    {
      id=id+1;
      old_id= G4RunManager::GetRunManager()->GetCurrentRun()->GetRunID();
    }
  
  if( LEMuSRDetectorConstruction::GetInstance()->cfthk!=thk_old)
    {
      id=0;
      thk_old=LEMuSRDetectorConstruction::GetInstance()->cfthk;
    } 
  muon.id = G4RunManager::GetRunManager()->GetCurrentRun()->GetRunID();

  if(p_name=="Mu")  muon.charge=0;
  else if(p_name=="mu+")  muon.charge=1;


}


//--------------------------------------------------------------------------------//
//                                            ROOT

void TDCheck::BookRoot()
{
  myFile = new TFile("gfoil.root", "RECREATE");

  tree = new TTree ("tree","Muons parameters");

  tree->Branch("muon",&muon.ienergy,"Init_Energy/F:kenergy/F:localtime/F:globaltime:propertime/F:positionx/F:positiony:positionz:init_angle:init_radius:momdirz:Angle/F:Plane/I:foil_thickness/F:runID/I:Charge/I");

}


void TDCheck::FillRoot()
{
  tree->Fill();
}

void TDCheck::WriteRoot()
{
  myFile->Write();
  myFile->Close();

}