si_syst Namespace Reference

Classes
class	ParticleHistory
class	ParticleSIManager
	Manager. More...
class	ParticleSIPropagator
struct	ParticleState
	A simple structure to represent a state object: position, direction, momentum, charge. More...
class	SISystInput
	A simple class to provide the data for the systematic itself: a collection of particles with the relevant history. More...
struct	SIXSecData
	A simple structure to repesent a pair: xsection value and its uncertainty. More...

Typedefs
typedef std::vector< std::pair< AnaTrueParticleB *, std::vector< AnaTrueParticleB * > > >	ParentsAndDaughters
	Basic structures to be able to associate a track with its "relatives".
typedef std::vector< std::pair< AnaTrueParticleB *, std::vector< AnaTrueParticleB * > > >::iterator	ParentsAndDaughtersIt

Enumerations
enum	InteractionType { kNoInteraction = 0, kElastic, kInelastic }
	Enums. More...
enum	MaterialEnum {   kCarbon = 0, kOxygen, kAluminium, kIron,   kAlG10, kWaterSystem, kG10, kFGDGlue,   kG10Roha, kCounter }

Functions
Double_t	GetZoverAMaterial (TGeoNode *)
Double_t	GetIMaterial (TGeoNode *node)
Double_t	GetIElement (Int_t Z)
MaterialEnum	GetMaterialEnum (const std::string &)
	Get the material enum given a material name.
Double_t	BetheBloch (const ParticleState &)
Double_t	BetheBloch (const ParticleState &, TGeoNode *)
void	TakeSmallStep (ParticleState &, const Double_t &)
	Propagate a state in a small step.
void	TakeSmallStep (ParticleState &, const Double_t &, TGeoNode *)

Variables
const TVector3	field = TVector3(0.2, 0., 0.)
	Constants. More...
const Double_t	K = 0.307075
	Bethe-Bloch constants.

Detailed Description

A general class to do various business related to secondary interactions: retrieve particles of intereset and

Enumeration Type Documentation

InteractionType

enum si_syst::InteractionType

Enums.

Interaction types, keep a wide list here, one can add what needed

Definition at line 45 of file SecondaryInteractionSystematic.hxx.

                      {
    kNoInteraction = 0, // particle does not interact in the volume of interest
    kElastic, 
    kInelastic 
  };

MaterialEnum

enum si_syst::MaterialEnum

Enumerate various materials of interest Will be moved to secondary interactions syst...

Definition at line 54 of file SecondaryInteractionSystematic.hxx.

                   {
    kCarbon = 0,  //  probably better to use exact Z values, but want to be able to do simple loops etc      
    kOxygen,        
    kAluminium,     
    kIron,          
    kAlG10,         
    kWaterSystem,
    kG10,
    kFGDGlue,
    kG10Roha, 
    kCounter
  };

Function Documentation

BetheBloch() [1/2]

Double_t si_syst::BetheBloch

(

const ParticleState &

)

Propagation Bethe-Bloch dEdX given particle state + TGeoNode (the particle is currently crossing) to specify material properties

BetheBloch() [2/2]

Double_t si_syst::BetheBloch	(	const ParticleState &	state,
		TGeoNode *	node
	)

Bethe-Bloch dEdX given particle state + TGeoNode (the particle is currently crossing) to specify material properties this is not to read the node if already avaiable

Definition at line 184 of file SecondaryInteractionSystematic.cxx.

                                                                      {
  //********************************************************************


  TGeoVolume   *volume   = node->GetVolume();
  TGeoMaterial *material = volume->GetMaterial();

  // Now, set up the values that we need for the dEdX
  Double_t betagamma  = state.Momentum/state.Mass;
  Double_t energy     = sqrt(state.Momentum*state.Momentum + state.Mass*state.Mass);
  Double_t beta       = state.Momentum/energy;
  Double_t gamma      = betagamma/beta;
  Double_t T_max      = (2*units::mass_electron*pow(betagamma,2))/(1 + 2*gamma*(units::mass_electron/state.Mass) + pow(units::mass_electron/state.Mass,2));

  // Potential
  Double_t I       = GetIMaterial(node)*(1E-6); //Make sure to convert to MeV!

  // The conversion factor converts the density stored in the geometry to g/cm^3
  Double_t rho     = (1000.0/(6.2415e21))*(material->GetDensity());

  Double_t ZoverA = si_syst::GetZoverAMaterial(node);
  //0.1 to convert to mm 
  return 0.1*rho*K*1.*(ZoverA)*(1/pow(beta,2))*((1.0/2.0)*log(2*units::mass_electron*pow(betagamma,2)*T_max/pow(I,2)) - pow(beta,2));

}

GetIElement()

Double_t si_syst::GetIElement

(

Int_t

Z

)

Takes in an atomic number and returns the Mean Ionization Potential for that element, from http://physics.nist.gov/PhysRefData/XrayMassCoef/tab1.html Only filled for elements present in ND280UserDetectorConstruction.cc Uses units of eV. Must be converted to MeV later.

Definition at line 74 of file SecondaryInteractionSystematic.cxx.

                                    {
  //********************************************************************
  switch(Z){    
    case 1:     //Hydrogen
      return 19.2;    
    case 5:      //Boron
      return 76.0;    
    case 6:      //Carbon
      return 78.0;    
    case 7:     //Nitrogen
      return 82.0;
    case 8:     //Oxygen
      return 95.0;    
    case 9:     //Fluorine
      return 115.0;
    case 11:    //Sodium
      return 149.0;    
    case 13:    //Aluminum
      return 166.0;
    case 14:    //Silicon
      return 173.0;    
    case 17:    //Chlorine
      return 174.0;
    case 18:    //Argon
      return 188.0;
    case 22:    //Titanium
      return 233.0;    
    case 26:    //Iron
      return 286.0;    
    case 27:    //Cobalt
      return 297.0;    
    case 29:    //Copper
      return 322.0;    
    case 30:    //Zinc
      return 330.0;    
    case 50:    //Tin
      return 488.0;    
    case 82:    //Lead
      return 823.0;    
    default:    //If unrecognized, just assume FGD scintillator.
      return 68.7;
  }    
}

GetIMaterial()

Double_t si_syst::GetIMaterial

(

TGeoNode *

node

)

Calculates the effective mean excitation energy for Equation taken from Leo p.29. Material properties retrieved from node

Definition at line 149 of file SecondaryInteractionSystematic.cxx.

                                            {
  //********************************************************************

  TGeoVolume   *volume   = node->GetVolume();
  TGeoMaterial *material = volume->GetMaterial();
  TGeoMixture  *mixture  = (TGeoMixture*)material;

  Double_t ZoverA = GetZoverAMaterial(node);

  Double_t lnI = 0;

  //Get the number of elements
  Int_t NElements = mixture->GetNelements();

  //Get the Z values for the material constituents.
  Double_t* Zmixt = mixture->GetZmixt();

  //Get the A values for the material constituents.
  Double_t* Amixt = mixture->GetAmixt();

  //Get the weight by mass of the elements.
  Double_t* Wmixt = mixture->GetWmixt();

  //lnIeff = sum(w_i*(Z_i/A_i)*lnI_i)/(Z/A)
  for(Int_t i = 0; i < NElements; i++){

    lnI += (Wmixt[i]*(Zmixt[i]/Amixt[i])*log(GetIElement((Int_t)Zmixt[i])))/ZoverA;
  }

  return exp(lnI);  
}

GetZoverAMaterial()

Double_t si_syst::GetZoverAMaterial

(

TGeoNode *

node

)

General kinematics methods ZoverA, material properties retrieved from node

Definition at line 119 of file SecondaryInteractionSystematic.cxx.

                                                 { 
  //******************************************************************** 
  TGeoVolume   *volume   = node->GetVolume();
  TGeoMaterial *material = volume->GetMaterial();
  TGeoMixture  *mixture  = (TGeoMixture*)material;

  Double_t ZoverA = 0;

  // Get the number of elements
  Int_t NElements = mixture->GetNelements();

  // Get the Z values for the material constituents.
  Double_t* Zmixt = mixture->GetZmixt();

  // Get the A values for the material constituents.
  Double_t* Amixt = mixture->GetAmixt();

  // Get the weight by mass of the elements.
  Double_t* Wmixt = mixture->GetWmixt();

  // Compute Z/A
  for(Int_t i = 0; i < NElements; i++){

    ZoverA += Wmixt[i]*(Zmixt[i]/Amixt[i]);

  }  
  return ZoverA;
}

TakeSmallStep()

void si_syst::TakeSmallStep	(	ParticleState &	state,
		const Double_t &	stepLength,
		TGeoNode *	node
	)

Propagate a state in a small step provided the node (material properties) and a step this is not to read the node if already avaiable

Definition at line 211 of file SecondaryInteractionSystematic.cxx.

                                                                                           {
  //********************************************************************

  //Find the initial kinetic energy.
  Double_t initEkin = anaUtils::ComputeKineticEnergy(state.Momentum, state.Mass);

  //Find the initial direction at the beginning of the step
  TVector3 initDir = state.Dir;

  Double_t p          = state.Momentum;
  Double_t betagamma  = p / state.Mass;
  Double_t beta       =  p /(initEkin + state.Mass);
  Double_t gamma      = betagamma / beta;

  Double_t dEdx = BetheBloch(state, node);

  //Remember, it is dE/dx * length in material * density of material
  //Used 1.032 g/cm^3, the density of the scintillator bars (from TN-91), as opposed
  //to the full module density which took into account the G10, glue, and fiber.
  //Factor of 0.1 to convert the mm tracklength to cm.
  Double_t finalEkin = initEkin - dEdx*(0.1)*stepLength;

  if (finalEkin <= 0 ){
    state.Dir = TVector3(0., 0., 0.);
    state.Momentum = 0.;
    return;
  }

  //Now, since we're taking small steps, we're assuming that the change in momentum is small, so we'll
  //just use the initial momentum as the feed through into the helix code.
  //The velocity of the pion in m/s)
  //PDG speed of light value: 299 792 458 m/s
  Double_t c = 299792458;
  TVector3 initVelocity = beta*c*initDir;

  //dp/dt = q(v x B), where p = \beta\gamma*mc
  //1 unit of charge is 1.602176e-19 Coulombs
  TVector3 dpdt = 1.*(1.602176e-19)*initVelocity.Cross(field);

  //So the change in momentum is therefore dpdt*dt, where dt = stepLength/initVelocity.Mag()
  //0.001 converts stepLength from mm to m.
  Double_t dt = ((0.001*stepLength)/initVelocity.Mag());
  TVector3 dp = dpdt*dt;

  //So that we're working all in the same units, use gamma*m*initVelocity to get the initial momentum back
  //(i.e. p = gamma*m*v)
  //Convert m_pi from MeV/c^2 to kg with factor of 1.782662e-30
  TVector3 pFinal = gamma*state.Mass*(1.782662e-30)*initVelocity + dp;

  //Now we know the final momentum in kg*m/s.  The Lorentz force only changes the direction and not the
  //magnitude of the momentum.  So, all we care about is the direction.
  state.Dir = pFinal*(1/pFinal.Mag());

  //The final position of a step is going to be
  //the result of going the step length in the direction of the initial momentum.
  //The final momentum is going to be the final kinetic energy converted back to a momentum
  //and in the direction of the final direction.
  // set new position
  state.Pos += stepLength*initDir;

  state.Momentum = sqrt(pow(finalEkin + state.Mass, 2) - pow(state.Mass, 2));

  return; 

}

Variable Documentation

field

const TVector3 si_syst::field = TVector3(0.2, 0., 0.)

Constants.

The magnetic field is 0.2 T in the x direction –> used for stepping –> keep it general

Definition at line 33 of file SecondaryInteractionSystematic.hxx.