ANALYSIS/FIN_PHYS.h

// @(#)fROOT/PREAN:$Name:  $:$Id: FIN_PHYS.h,v 1.5 2007/09/25 14:51:51 Diego_Faso Exp $
// Author: Diego Faso <mailto:faso@to.infn.it>, 2006/11/09

#ifndef FROOT_FIN_PHYS
#define FROOT_FIN_PHYS

// The FINUDA physics namespace //

#include <TString.h>
#include <TMath.h>
#include <TROOT.h>

namespace FIN_PHYS {

  enum E_FndDataTaking {
    DTAK_2003_2004 = 0, 
    DTAK_2006_2007 = 1, 
  };

  enum E_FndFidaVersion {
    FIDAVER_521 = 0, 
    FIDAVER_601 = 1, // version 600 has been skipped
    FIDAVER_603 = 2, // version 602 has been skipped (development version)
    FIDAVER_END = 3, // used for arrays
  };

  enum E_FinPhys_BhabhaPart_ID{ // negative particles at id=0
    FPh_BhaElec_id = 0,
    FPh_BhaPosit_id = 1,
  };

  enum E_FinPhys_HypePart_ID{ // negative particles at id=0
    FPh_HypKmin_id = 0,
    FPh_HypKplu_id = 1,
  };
  
  enum E_Fnd_PID { // summary of Particle ID codes (geant convention)
    FPh_PID_Positron =     2, // -11 (PDG number)  
    FPh_PID_Electron =     3, //  11
    //
    FPh_PID_Muon_plu =     5, // -13  
    FPh_PID_Muon_min =     6, // 13   
    //
    FPh_PID_Pi_zero =  7, //  111   
    FPh_PID_Pi_plu  =  8, //  211   
    FPh_PID_Pi_min  =  9, // -211  
    //
    FPh_PID_Kaon_Long  =    10, // 130   
    FPh_PID_Kaon_Short =    16, // 310,   
    FPh_PID_Kaon_plu   =   11, // 321  
    FPh_PID_Kaon_min   =   12, // -321
    //
    FPh_PID_Neutron  =  13, // 2112  
    FPh_PID_Proton   =  14, // 2212 
    FPh_PID_AProton  =  15, // -2212, 
    FPh_PID_ANeutron =  25, // -2112, 
    //
    FPh_PID_Lambda  = 18, //3122,  
    FPh_PID_ALambda = 26,
    //
    FPh_PID_Sigma_plu   =  19, // 3222,  
    FPh_PID_Sigma_zero  =  20, // 3212,  
    FPh_PID_Sigma_min   =  21, // 3112,  
    FPh_PID_ASigma_plu  =  29, // -3222, 
    FPh_PID_ASigma_zero =  28, // -3212, 
    FPh_PID_ASigma_min  =  27, // -3112, 
    //
    FPh_PID_Gamma    =  1, // 22
    FPh_PID_Deuteron = 45, //
    FPh_PID_Tritium  = 46, //
    FPh_PID_Alpha    = 47, //
    //
  }; 
  
  
  // methods declaration
  TString FidaVer_Name(E_FndFidaVersion ver = FIDAVER_603);
  Int_t   FidaVer_ID(E_FndFidaVersion ver   = FIDAVER_603);

  Double_t TrackRad2Mom(Double_t B=1);
  Double_t TrackMom2Rad(Double_t B=1);
  Double_t TrackRad_FromMom(Double_t mom,Double_t lam,Double_t B=1);

  Double_t GetBhaRate_Lumin_Factor();

  Double_t GetParticleMass(E_Fnd_PID pid); // Gev/c2
  //
  Double_t Sum(const Float_t &mod1,const Float_t &cdx1,const Float_t &cdy1,const Float_t &cdz1,
               const Float_t &mod2,const Float_t &cdx2,const Float_t &cdy2,const Float_t &cdz2,
               Double_t &cdx_res,Double_t &cdy_res,Double_t &cdz_res);
  //
  void ScalarProduct(const Float_t &mod1,const Float_t &cdx1,const Float_t &cdy1,const Float_t &cdz1,
                     const Float_t &mod2,const Float_t &cdx2,const Float_t &cdy2,const Float_t &cdz2,
                     Double_t &module,Double_t &angle_rad);
  //
  Double_t ParticleEnergy_2(E_Fnd_PID pid,const Float_t &mom); // E^2
  Double_t ParticleEnergy(E_Fnd_PID pid,const Float_t &mom);
  //
  Double_t KineticEnergy_2(const Float_t &mom1,const Float_t &cdx1,const Float_t &cdy1,const Float_t &cdz1,
                           const Float_t &mom2,const Float_t &cdx2,const Float_t &cdy2,const Float_t &cdz2);
  Double_t KineticEnergy(const Float_t &mom1,const Float_t &cdx1,const Float_t &cdy1,const Float_t &cdz1,
                         const Float_t &mom2,const Float_t &cdx2,const Float_t &cdy2,const Float_t &cdz2);
  //
  Double_t CenterMassEnergy(E_Fnd_PID pid1, const Float_t &mom1,const Float_t &cdx1,const Float_t &cdy1,const Float_t &cdz1,
                            E_Fnd_PID pid2, const Float_t &mom2,const Float_t &cdx2,const Float_t &cdy2,const Float_t &cdz2);
  
  
  // --- reference system change
  void RefSys_CosDir_To_Spheric(const Double_t &cx,const Double_t &cy,const Double_t &cz,Double_t &theta,Double_t &phi);
  void RefSys_Sphere_To_CosDir(const Double_t &theta,const Double_t &phi,Double_t &cx,Double_t &cy,Double_t &cz);

  TString HolleritToString(Int_t hollerit_in); // by Hyeongryul So (2003)
  
  // GENERAL METHODS USED FOR OFFLINE ANALYSIS //
  //
  Double_t Dedxfunc(Double_t Momentum, Double_t Mass, Double_t Charge, Double_t Zval, Double_t Aval, Double_t Ival);
  Double_t DedxfuncSi(Double_t Momentum, Double_t Mass, Double_t Charge);
  Double_t DedxfuncCh(Double_t Momentum, Double_t Mass, Double_t Charge);
  //
  Double_t Tenevsmom(Double_t Momentum, Double_t Mass);
  //
};

// Implementation will follow:

//_______________________________
inline Int_t FIN_PHYS::FidaVer_ID(E_FndFidaVersion ver){

  return (Int_t) ver;
}

//_______________________________
inline  TString FIN_PHYS::FidaVer_Name(E_FndFidaVersion ver){
  // upgraded everytime a new version of fidarc/mc is released
  switch (ver){
  case FIDAVER_521: return "v_521";
  case FIDAVER_601: return "v_601";
  case FIDAVER_603: return "v_603";
  default: return "";
  }
}

//__________________________
inline Double_t FIN_PHYS::TrackRad2Mom(Double_t B){
  // transv_mom = radius * <this factor>
  //
  // formula: p = 0.00300 * R * B
  //  [p] = GeV/c
  //  [R] = cm   (meters)
  //  [B] = T   (tesla)

  return (Double_t)(0.00300 * B);
}

//__________________________
inline Double_t FIN_PHYS::TrackMom2Rad(Double_t B){
  // radius = transv_mom * <this factor>

  return (Double_t)( 1. / TrackRad2Mom(B) );
}

//__________________________
inline Double_t FIN_PHYS::TrackRad_FromMom(Double_t mom,Double_t lam,Double_t B){

  Double_t mom_transv = mom * TMath::Cos(lam);
  Double_t rad = mom_transv * TrackMom2Rad(B);
  return rad;
}

//__________________________
inline Double_t FIN_PHYS::GetBhaRate_Lumin_Factor(){
  // this value must be multiplied by the BHABHA rate
  //  in order to get the instantaneous luminosity

  return 1.4e+30;
}

//__________________________
inline Double_t FIN_PHYS::GetParticleMass(E_Fnd_PID pid){
  // Get mass from PID code (Gev/c2)
  // returns "-1" in case pid is not supported
  
  switch(pid){
  case FPh_PID_Positron : return 0.00051099906;
  case FPh_PID_Electron : return 0.00051099906;
    //
  case FPh_PID_Muon_plu : return 0.105658389;
  case FPh_PID_Muon_min : return 0.105658389;
    //
  case FPh_PID_Pi_zero : return 0.1349764;
  case FPh_PID_Pi_plu  : return 0.1395700;
  case FPh_PID_Pi_min  : return 0.1395700;
    //
  case FPh_PID_Kaon_Long  : return 0.497672;
  case FPh_PID_Kaon_Short : return 0.497672;
  case FPh_PID_Kaon_plu   : return 0.493677;
  case FPh_PID_Kaon_min   : return 0.493677;
    //
  case FPh_PID_Neutron  : return 0.93956563;
  case FPh_PID_Proton   : return 0.93827231;
  case FPh_PID_AProton  : return 0.93827231;
  case FPh_PID_ANeutron : return 0.93956563;
    //
  case FPh_PID_Lambda   : return 1.115684;
  case FPh_PID_ALambda  : return 1.115684;
    //
  case FPh_PID_Sigma_plu    : return 1.18937;
  case FPh_PID_Sigma_zero   : return 1.19255;
  case FPh_PID_Sigma_min    : return 1.197436;
  case FPh_PID_ASigma_plu   : return 1.197436;
  case FPh_PID_ASigma_zero  : return 1.19255;
  case FPh_PID_ASigma_min   : return 1.18937;
    //
  case FPh_PID_Gamma     : return 0;
  case FPh_PID_Deuteron  : return 1.875613;
  case FPh_PID_Tritium   : return 2.80925;
  case FPh_PID_Alpha     : return 3.727417;
    //
  default:
    gROOT->Warning("FIN_PHYS::GetParticleMass","PID not suported (%d)",pid);
    return -1;
  }
}

//__________________________
inline Double_t FIN_PHYS::Sum(const Float_t &mod1,const Float_t &cdx1,const Float_t &cdy1,const Float_t &cdz1,
                                        const Float_t &mod2,const Float_t &cdx2,const Float_t &cdy2,const Float_t &cdz2,
                                        Double_t &cdx_res,Double_t &cdy_res,Double_t &cdz_res){
  // vectorial sum of two vectors
  // mod[1-2] : vector modules
  // cd*[1-3] : director cosines

  Double_t module = TMath::Sqrt(
                                mod1*mod1 + mod2*mod2 + 
                                2 * (mod1*cdx1*mod2*cdx2+mod1*cdy1*mod2*cdy2+mod1*cdz1*mod2*cdz2)
                                );
  
  cdx_res = (Double_t) cdx1 + cdx2;
  cdy_res = (Double_t) cdy1 + cdy2;
  cdz_res = (Double_t) cdz1 + cdz2;

  return module;

}

//__________________________
inline void FIN_PHYS::ScalarProduct(const Float_t &mod1,const Float_t &cdx1,const Float_t &cdy1,const Float_t &cdz1,
                                        const Float_t &mod2,const Float_t &cdx2,const Float_t &cdy2,const Float_t &cdz2,
                                        Double_t &module,Double_t &angle_rad){
  // mod[1-2] : vector modules
  // cd*[1-3] : director cosines
  
  if(! mod1 || !mod2){
    module = 0;
    angle_rad = -999;
    return;
  }

  module = (Double_t) ( 
                       mod1 * cdx1 * mod2 * cdx2 + 
                       mod1 * cdy1 * mod2 * cdy2 + 
                       mod1 * cdz1 * mod2 * cdz2
                       );
  
  Double_t cos_ang = module / ( (Double_t) (mod1*mod2) );
  if(module != 0){ // can not evaluate the angle between null vectors
    angle_rad = TMath::ACos(cos_ang);
  }
  else angle_rad = -999;

  //  cout << "module sign: " << module << endl;
  module = TMath::Abs(module);
  //  cout << "module abs: " << module << endl << endl;

  return;
}

//__________________________
inline Double_t FIN_PHYS::KineticEnergy_2(const Float_t &mom1,const Float_t &cdx1,const Float_t &cdy1,const Float_t &cdz1,
                                           const Float_t &mom2,const Float_t &cdx2,const Float_t &cdy2,const Float_t &cdz2){
  // compute the (kinetic energy^2) (Gev^2/c^2) of the given system
  // (according to the given cosines angles)
  
  Double_t kin_en = ((mom1*cdx1 + mom2*cdx2) * (mom1*cdx1 + mom2*cdx2) + 
                     (mom1*cdy1 + mom2*cdy2) * (mom1*cdy1 + mom2*cdy2) + 
                     (mom1*cdz1 + mom2*cdz2) * (mom1*cdz1 + mom2*cdz2)
                     );
  return kin_en;
  
}

//__________________________
inline Double_t FIN_PHYS::KineticEnergy(const Float_t &mom1,const Float_t &cdx1,const Float_t &cdy1,const Float_t &cdz1,
                                        const Float_t &mom2,const Float_t &cdx2,const Float_t &cdy2,const Float_t &cdz2){
  // compute the (kinetic energy) (Gev/c) of the given system
  // (according to the given cosines angles)
  
  return TMath::Sqrt( KineticEnergy_2(mom1,cdx1,cdy1,cdz1,mom2,cdx2,cdy2,cdz2) );
}

//__________________________
inline Double_t FIN_PHYS::CenterMassEnergy(E_Fnd_PID pid1,
                                           const Float_t &mom1,const Float_t &cdx1,const Float_t &cdy1,const Float_t &cdz1,
                                           E_Fnd_PID pid2,
                                           const Float_t &mom2,const Float_t &cdx2,const Float_t &cdy2,const Float_t &cdz2){
  // compute the invariant mass (Gev/c2) of the [pid1,pid2] system
  // (according to the given momenta)
  
  Double_t ene1 = ParticleEnergy(pid1,mom1); // take care: E^2
  Double_t ene2 = ParticleEnergy(pid2,mom2); // take care: E^2
  
  if(ene1 == 0 || ene2 ==0) return 0;

  Double_t Imass = TMath::Sqrt(
                               (ene1+ene2)*(ene1+ene2) - 
                               KineticEnergy_2(mom1,cdx1,cdy1,cdz1,mom2,cdx2,cdy2,cdz2)
                               );
  
  return Imass;
  
}

//__________________________
inline Double_t FIN_PHYS::ParticleEnergy_2(E_Fnd_PID pid,const Float_t &mom){
  // compute the energy^2 (GeV^2) of the [pid] particle
  // (according to the given momentum)
  
  Double_t mass = GetParticleMass(pid);
  
  Double_t energy_2 = (Double_t) ( (mass * mass) + (mom * mom) );
  return energy_2;
}

//__________________________
inline Double_t FIN_PHYS::ParticleEnergy(E_Fnd_PID pid,const Float_t &mom){
  // compute the energy (GeV) of the [pid] particle
  // (according to the given momentum)

  return  TMath::Sqrt( ParticleEnergy_2(pid,mom) );
}

//__________________________
inline TString FIN_PHYS::HolleritToString(Int_t hollerit_in){
  // by Hyeongryul So (year 2003)
  // modified to TString by D. Faso (November 2006)
  // needed while reading content of zebra (hollerit code is used)

/*   result[0] = (Char_t) ((hollerit_in >>  0) & 0xFF); */
/*   result[1] = (Char_t) ((hollerit_in >>  8) & 0xFF); */
/*   result[2] = (Char_t) ((hollerit_in >> 16) & 0xFF); */
/*   result[3] = (Char_t) ((hollerit_in >> 24) & 0xFF); */
/*   result[4] = '\0'; */

  TString res = "";
  res += (Char_t) ((hollerit_in >>  0) & 0xFF);
  res += (Char_t) ((hollerit_in >>  8) & 0xFF);
  res += (Char_t) ((hollerit_in >> 16) & 0xFF);
  res += (Char_t) ((hollerit_in >> 24) & 0xFF);
  res.Resize(4);

/*   TString res; */
/*   res.Form("%s%s%s%s", */
/*         (Char_t) ((hollerit_in >>  0) & 0xFF), */
/*         (Char_t) ((hollerit_in >>  8) & 0xFF), */
/*         (Char_t) ((hollerit_in >> 16) & 0xFF), */
/*         (Char_t) ((hollerit_in >> 24) & 0xFF) */
/*         ); */

  return res;
    
}


// --- reference system change
//__________________________
inline void FIN_PHYS::RefSys_CosDir_To_Spheric(const Double_t &cx,const Double_t &cy,const Double_t &cz,Double_t &theta,Double_t &phi){
  // all angles are expressed in rad
  // cx,cy,cz: director-cosines
  // theta: angle to the z-axis (spherical coordinates)  
  // phi: angle to the x-axis in the xy plane (spherical coordinates)  
  
  theta = (Double_t)TMath::ACos(cz);

  if(cx != 0){
    phi = (Double_t) TMath::ATan(cy/cx);
    if(cx < 0) phi = (Double_t)(TMath::Pi()) + phi;
  }
  else{
    Double_t signe = (cy > 0) ? +1 : -1;
    phi = signe * (Double_t)(TMath::Pi() / 2);
  }
  
}

//__________________________
inline void FIN_PHYS::RefSys_Sphere_To_CosDir(const Double_t &theta,const Double_t &phi,Double_t &cx,Double_t &cy,Double_t &cz){
  // all angles are expressed in rad
  // cx,cy,cz: director-cosines
  // theta: angle to the z-axis (spherical coordinates)  
  // phi: angle to the x-axis in the xy plane (spherical coordinates)  
  
  cx = (Double_t) ( TMath::Sin(theta) * TMath::Cos(phi) );
  cy = (Double_t) ( TMath::Sin(theta) * TMath::Sin(phi) );
  cz = (Double_t) ( TMath::Cos(theta) );
}

// GENERAL METHODS USED FOR OFFLINE ANALYSIS //

//__________________________
inline Double_t FIN_PHYS::Dedxfunc(Double_t Momentum, Double_t Mass, Double_t Charge, Double_t Zval, Double_t Aval, Double_t Ival)
{
  // Momentum in GeV/c, Mass in GeV/c^2, Charge in elementary charge 
  Momentum *= 1000.;
  Mass *= 1000. ;
  // Momentum in MeV/c, Mass in MeV/c^2, Charge in elementary charge 
  // Zval Atomic Number of Medium
  // Aval Atomic Mass of Medium 
  // Ival mean Excitation Energy [eV]
  Double_t Energy = TMath::Sqrt(TMath::Power(Momentum,2.) + TMath::Power(Mass,2.));
  Double_t Beta = Momentum / Energy;
  Double_t Beta2 = TMath::Power(Beta,2.);
  Double_t Gamma = Energy / Mass;
  Double_t Gamma2 = TMath::Power(Gamma,2.);
  Double_t Dcons = 0.30707;
  //
  Double_t melect = GetParticleMass(FPh_PID_Electron) * 1000.;
  // Double_t relect = 2.817940285;
  Double_t ZA = Zval/Aval;
  Double_t I2 = TMath::Power(Ival,2.);
  //
  Double_t Bethe = Dcons * ZA * TMath::Power(Charge,2.) / Beta2;
  Double_t TMax = (2*melect*Beta2*Gamma2);
  TMax /= ((1+2*Gamma*melect/Mass)+(TMath::Power((melect/Mass),2)));
  Double_t ArgLog = (2*melect*Beta2*Gamma2*TMax/I2);
  Bethe *= (0.5 * TMath::Log(ArgLog) - Beta2);
  return Bethe;
}

//__________________________
inline Double_t FIN_PHYS::DedxfuncSi(Double_t Momentum, Double_t Mass, Double_t Charge)
{
  // Silicon
  //  Double_t rho = 2.33;
  Double_t Zval = 14.;
  Double_t Aval = 28.0855;
  Double_t Ival = 173.0E-6;
  //
  return Dedxfunc(Momentum, Mass, Charge, Zval, Aval, Ival);
}

//__________________________
inline Double_t FIN_PHYS::DedxfuncCh(Double_t Momentum, Double_t Mass, Double_t Charge)
{
  // Silicon
  //  Double_t rho = 2.33;
  Double_t Zval = 11.6;
  Double_t Aval = 20.2;
  Double_t Ival = 278.9E-6;
  //
  return Dedxfunc(Momentum, Mass, Charge, Zval, Aval, Ival);
}

//__________________________
inline Double_t FIN_PHYS::Tenevsmom(Double_t Momentum, Double_t Mass)
{
  return TMath::Sqrt(TMath::Power(Momentum,2.)+TMath::Power(Mass,2.)) - Mass;
}

//
#endif // FROOT_FIN_PHYS

---