MrpcCalibSvc.cxx

Go to the documentation of this file.

//********************************************************
//
//  MRPC Reconstruction Service - Provides Calibration constants
//  Matthias Ullrich - <[email protected]>
//
//  Description:
//     Class Provides Calibration constant, which are use in MRPC Reconstruction:
//     - Time-Amplitude Correlation is used in TofRawDataProviderSvc
//     - Sigma is used in MrpcRec and MrpcDBSRec
//     - Energy-Calibration is used in TofEnergyRec  
//*********************************************************
 
 
#include "GaudiKernel/Kernel.h"
#include "GaudiKernel/IInterface.h"
#include "GaudiKernel/IIncidentSvc.h"
#include "GaudiKernel/Incident.h"
#include "GaudiKernel/IIncidentListener.h"
#include "GaudiKernel/StatusCode.h"
#include "GaudiKernel/SvcFactory.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/DataSvc.h"
 
#include "RawEvent/RawDataUtil.h"
 
#include "MrpcCalibSvc/MrpcCalibSvc.h"
 
 
using namespace std;
 
 
 
 
 
 
 
 
//Time Amplitude Correction extracted in BOSS 6.6.4a
double MrpcCalibSvc::Time_Amplitude_Correction(unsigned int tdcChannel, unsigned int tdc2Channel)
{
 
  double pulse_length=RawDataUtil::TofTime(tdc2Channel-tdcChannel);
 
 
  //Calculate the correction:
  double correction_ta=0;
 
  if(pulse_length < 4.164){correction_ta=-100000.;}   // 4.164 ns = 10 fC; Everything smaller is expected to be noise! ->Deposit charge =0;
  
  //This function has been extracted with a fit on MC data
  if(pulse_length>=4.164 && pulse_length<=14.0)
    {
      correction_ta=278.03; //value of the function below at pulselength=14.0 as it is ~constant.
    }
  else if(pulse_length>14.0 && 15.793)
    {
      correction_ta = 1152674.673-307911.4654*pulse_length+30788.47591*pulse_length*pulse_length-1365.232906*pulse_length*pulse_length*pulse_length+22.64702791*pulse_length*pulse_length*pulse_length*pulse_length;
      if(correction_ta<0)correction_ta=0;
    }
  else if(pulse_length >= 15.793)
    {
      6104.055949-375.6901764*pulse_length; // Linear fit for the last part of the spectra. The x^4 - polynom was not able to describe this part satisfying.
      if(correction_ta<0)correction_ta=0;
    }
 
  return (correction_ta/1000.); 
 
}
 
 
 
 
 
 
// All the folloqing values have been extracted in BOSS 6.6.2
 
//Get the resolution of (t_measured-t_estimated-t_expected-t_tacorrection-t_transitiontime) <--- This should be the correct distribution 
// which should be used for the PID!
//The correction for ta is done in TOF-RawDataProvider, the correction for transitiontime in MRPRREC
// The time --- t_measured ---  stored in the files eventually is  the time corrected for transition time, TA correction and t_estimated.
//DBS is true, when one use the double-sided-readout version of MRPC
 
double MrpcCalibSvc::GetSigma(double momentum, unsigned tof_id,  int particle_species, bool dbs) 
{
  //paricle species: 0=Electron, 1=Muon, 2=Pion, 3=Kaon, 4 = Proton
 
  double sigma=0;
 
  if(particle_species==0) sigma=GetSigmaElectron(momentum,tof_id,dbs);
  if(particle_species==1) sigma=GetSigmaMuon(momentum,tof_id,dbs);
  if(particle_species==2) sigma=GetSigmaPion(momentum,tof_id,dbs);
  if(particle_species==3) sigma=GetSigmaKaon(momentum,tof_id,dbs);
  if(particle_species==4) sigma=GetSigmaProton(momentum,tof_id,dbs);
 
  return sigma;
}
 
 
 
double MrpcCalibSvc::GetSigmaElectron(double momentum, unsigned tof_id, bool dbs) 
{
  double sigma=0.08;
 
  if(!dbs) sigma =0.0864116*exp(-momentum/0.329858)+0.0898002;
  if(dbs)  sigma =0.103841*exp(-momentum/0.252881)+0.0840766;
 
  return sigma;
}
 
 
double MrpcCalibSvc::GetSigmaMuon(double momentum, unsigned tof_id, bool dbs) 
{
  double sigma=0.08;
  
  if(!dbs) sigma=0.165892*exp(-momentum/0.25094)+0.0718477;
  if(dbs)  sigma=0.166412*exp(-momentum/0.203065)+0.0718146;
 
  return sigma;
}
 
 
 
double MrpcCalibSvc::GetSigmaPion(double momentum, unsigned tof_id, bool dbs) 
{
  double sigma=0.08;
 
  if(!dbs) sigma  =0.338622*exp(-momentum/0.187429)+0.0802867;
  if(dbs)  sigma  =0.448576*exp(-momentum/0.147072)+0.078237;
 
  return sigma;
}
 
 
 
double MrpcCalibSvc::GetSigmaKaon(double momentum, unsigned tof_id, bool dbs) 
{
  double sigma=0.08;
 
  if(!dbs) sigma =0.174625*exp(-momentum/0.201226+1.13756+1.10314*momentum*momentum)+0.0752667;
  if(dbs)  sigma =0.581175*exp(-momentum/0.141733+0.581804+1.75623*momentum*momentum)+0.0767948;
     
  return sigma;
}
 
 
 
double MrpcCalibSvc::GetSigmaProton(double momentum, unsigned tof_id, bool dbs) 
{
  double sigma=0.08;
 
  if(!dbs) sigma=0.8125*exp(-momentum/0.276857)+0.0828829;
  if(dbs)  sigma=1.1716*exp(-momentum/0.227311)+0.0823852;
 
  return sigma;
}
 
 
 
 
 
//Calibration value for TofEnergyRec
 
double MrpcCalibSvc::GetEnergyCalibration(double charge_ns, bool neutral , bool neighborhood)
{
 
  double energy=0;
 
  if(neutral==true)
    {
      // The values are from photons simulated with geant4!                                                                                                                                                                        
      if(charge_ns<15.0716) energy = 326.35-10.6284*15.0716-1.4169*15.0716*15.0716-0.0514292*15.0716*15.0716*15.0716+0.00645421*15.0716*15.0716*15.0716*15.0716;//in MeV: Below this value we assume a constant energyloss!        
      else if(charge_ns>15.9131) energy = -7.91205+0.814588*charge_ns; //Over this value we expect the energyloss to be the same as for charged particles!                                                                         
      else energy = 326.35-10.6284*charge_ns-1.4169*charge_ns*charge_ns-0.0514292*charge_ns*charge_ns*charge_ns+0.00645421*charge_ns*charge_ns*charge_ns*charge_ns;// in MeV                                                       
    }
  else //charged tracks, (combined electrons/positrons in geant4)                                                                                                                                                                  
    {
      if(neighborhood==false) energy = -7.91205+0.814588*charge_ns;// in MeV                                                                                                                                                       
      else
    {
      if(15.16>= charge_ns)  energy= -72550.7+14230.9*15.16-930.499*15.16*15.16+20.2818*15.16*15.16*15.16;//Minimal loss from fit!                                                                                             
      if(15.16< charge_ns &&charge_ns<15.7633) energy= -72550.7+14230.9*charge_ns-930.499*charge_ns*charge_ns+20.2818*charge_ns*charge_ns*charge_ns;
      else energy = -7.91205+0.814588*charge_ns; //Else we will use the not neighborhood distribution as our fit will begin to fall                                                                                            
    }
 
    }
 
  return energy;
 
}