GFEnergyLossCoulomb.cxx

Go to the documentation of this file.

/* Copyright 2008-2009, Technische Universitaet Muenchen,
   Authors: Christian Hoeppner & Sebastian Neubert

   This file is part of GENFIT.

   GENFIT is free software: you can redistribute it and/or modify
   it under the terms of the GNU Lesser General Public License as published
   by the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.

   GENFIT is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU Lesser General Public License for more details.

   You should have received a copy of the GNU Lesser General Public License
   along with GENFIT.  If not, see <http://www.gnu.org/licenses/>.
*/

#include "larreco/Genfit/GFEnergyLossCoulomb.h"
#include <string>
#include "larreco/Genfit/GFException.h"
#include "math.h"

genf::GFEnergyLossCoulomb::~GFEnergyLossCoulomb(){
}

double genf::GFEnergyLossCoulomb::energyLoss(const double& step,
                                       const double& mom,
                                       const int&    pdg,              
                                       const double& /* matDensity */,
                                       const double& matZ,
                                       const double& /* matA */,
                                       const double& radiationLength,
                                       const double& /* meanExcitationEnergy */,
                                       const bool&   doNoise,
                                             TMatrixT<Double_t>* noise,
                                       const TMatrixT<Double_t>* jacobian,
                                       const TVector3* directionBefore,
                                       const TVector3* directionAfter){

  double charge, mass;
  getParticleParameters(pdg, charge, mass);
  
  const double beta = mom/sqrt(mass*mass+mom*mom);
 
  if (doNoise) {
    if (!noise) throw GFException(std::string(__func__) + ": no noise given!", __LINE__, __FILE__).setFatal();
    if (!jacobian) throw GFException(std::string(__func__) + ": no jacobian given!", __LINE__, __FILE__).setFatal();
    if (!directionBefore) throw GFException(std::string(__func__) + ": no directionBefore given!", __LINE__, __FILE__).setFatal();
    if (!directionAfter) throw GFException(std::string(__func__) + ": no directionAfter given!", __LINE__, __FILE__).setFatal();
    
    // MULTIPLE SCATTERING; calculate sigma^2 
    // PANDA report PV/01-07 eq(43); linear in step length
    double sigma2 = 225.E-6/(beta*beta*mom*mom) * step/radiationLength * matZ/(matZ+1) * log(159.*pow(matZ,-1./3.))/log(287./sqrt(matZ)); // sigma^2 = 225E-6/mom^2 * XX0/beta^2 * Z/(Z+1) * ln(159*Z^(-1/3))/ln(287*Z^(-1/2)
        
    
    // noiseBefore
      TMatrixT<Double_t> noiseBefore(7,7);
      
      // calculate euler angles theta, psi (so that directionBefore' points in z' direction)
      double psi = 0;
      if (fabs((*directionBefore)[1]) < 1E-14) {  // numerical case: arctan(+-inf)=+-PI/2
        if ((*directionBefore)[0] >= 0.) psi = M_PI/2.; 
        else psi = 3.*M_PI/2.;
      }
      else {
        if ((*directionBefore)[1] > 0.) psi = M_PI - atan((*directionBefore)[0]/(*directionBefore)[1]);
        else psi = -atan((*directionBefore)[0]/(*directionBefore)[1]);
      }
      // cache sin and cos
      double sintheta = sqrt(1-(*directionBefore)[2]*(*directionBefore)[2]);  // theta = arccos(directionBefore[2])
      double costheta = (*directionBefore)[2];  
      double sinpsi = sin(psi);
      double cospsi = cos(psi);
   
      // calculate NoiseBefore Matrix R M R^T  
      double noiseBefore34 =  sigma2 * cospsi * sinpsi * sintheta*sintheta; // noiseBefore_ij = noiseBefore_ji
      double noiseBefore35 = -sigma2 * costheta * sinpsi * sintheta;
      double noiseBefore45 =  sigma2 * costheta * cospsi * sintheta;

      noiseBefore[3][3] = sigma2 * (cospsi*cospsi + costheta*costheta - costheta*costheta * cospsi*cospsi);
      noiseBefore[4][3] = noiseBefore34;
      noiseBefore[5][3] = noiseBefore35;

      noiseBefore[3][4] = noiseBefore34;
      noiseBefore[4][4] = sigma2 * (sinpsi*sinpsi + costheta*costheta * cospsi*cospsi);
      noiseBefore[5][4] = noiseBefore45;

      noiseBefore[3][5] = noiseBefore35;
      noiseBefore[4][5] = noiseBefore45;
      noiseBefore[5][5] = sigma2 * sintheta*sintheta;
      
      TMatrixT<Double_t> jacobianT(7,7);
      jacobianT = (*jacobian);
      jacobianT.T();
    
      noiseBefore = jacobianT*noiseBefore*(*jacobian); //propagate
    
    // noiseAfter
      TMatrixT<Double_t> noiseAfter(7,7);
      
      // calculate euler angles theta, psi (so that A' points in z' direction)
      psi = 0;
      if (fabs((*directionAfter)[1]) < 1E-14) {  // numerical case: arctan(+-inf)=+-PI/2
        if ((*directionAfter)[0] >= 0.) psi = M_PI/2.; 
        else psi = 3.*M_PI/2.;
      }
      else {
        if ((*directionAfter)[1] > 0.) psi = M_PI - atan((*directionAfter)[0]/(*directionAfter)[1]);
        else psi = -atan((*directionAfter)[0]/(*directionAfter)[1]);
      }
      // cache sin and cos
      sintheta = sqrt(1-(*directionAfter)[2]*(*directionAfter)[2]);  // theta = arccos(directionAfter[2])
      costheta = (*directionAfter)[2];  
      sinpsi = sin(psi);
      cospsi = cos(psi);
   
      // calculate NoiseAfter Matrix R M R^T  
      double noiseAfter34 =  sigma2 * cospsi * sinpsi * sintheta*sintheta; // noiseAfter_ij = noiseAfter_ji
      double noiseAfter35 = -sigma2 * costheta * sinpsi * sintheta;
      double noiseAfter45 =  sigma2 * costheta * cospsi * sintheta;

      noiseAfter[3][3] = sigma2 * (cospsi*cospsi + costheta*costheta - costheta*costheta * cospsi*cospsi);
      noiseAfter[4][3] = noiseAfter34;
      noiseAfter[5][3] = noiseAfter35;

      noiseAfter[3][4] = noiseAfter34;
      noiseAfter[4][4] = sigma2 * (sinpsi*sinpsi + costheta*costheta * cospsi*cospsi);
      noiseAfter[5][4] = noiseAfter45;

      noiseAfter[3][5] = noiseAfter35;
      noiseAfter[4][5] = noiseAfter45;
      noiseAfter[5][5] = sigma2 * sintheta*sintheta;
      
    //calculate mean of noiseBefore and noiseAfter and update noise
      (*noise) += 0.5*noiseBefore + 0.5*noiseAfter;
  }
  
  return 0.;
}

//ClassImp(GFEnergyLossCoulomb)