STEER/AliITSPidParams.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id: */
  17
  18 ///////////////////////////////////////////////////////////////////
  19 //                                                               //
  20 // Implementation of the class to store parameters of ITS        //
  21 // response funcions for dE/dx based PID                         //
  22 // Origin: F.Prino, Torino, prino@to.infn.it                     //
  23 //                                                               //
  24 ///////////////////////////////////////////////////////////////////
  25
  26 #include <TFormula.h>
  27 #include <TNamed.h>
  28 #include <TMath.h>
  29 #include "AliITSPidParams.h"
  30 #include "AliPID.h"
  31
  32 ClassImp(AliITSPidParams)
  33
  34 //______________________________________________________________________
  35 AliITSPidParams::AliITSPidParams():
  36   TNamed("default",""),
  37   fSDDPionMPV(0),
  38   fSDDPionLandauWidth(0),
  39   fSDDPionGaussWidth(0),
  40   fSSDPionMPV(0),
  41   fSSDPionLandauWidth(0),
  42   fSSDPionGaussWidth(0),
  43   fSDDKaonMPV(0),
  44   fSDDKaonLandauWidth(0),
  45   fSDDKaonGaussWidth(0),
  46   fSSDKaonMPV(0),
  47   fSSDKaonLandauWidth(0),
  48   fSSDKaonGaussWidth(0),
  49   fSDDProtMPV(0),
  50   fSDDProtLandauWidth(0),
  51   fSDDProtGaussWidth(0),
  52   fSSDProtMPV(0),
  53   fSSDProtLandauWidth(0),
  54   fSSDProtGaussWidth(0)
  55 {
  56   // default constructor
  57   InitMC();
  58 }
  59 //______________________________________________________________________
  60 AliITSPidParams::AliITSPidParams(Char_t * name):
  61   TNamed(name,""),
  62   fSDDPionMPV(0),
  63   fSDDPionLandauWidth(0),
  64   fSDDPionGaussWidth(0),
  65   fSSDPionMPV(0),
  66   fSSDPionLandauWidth(0),
  67   fSSDPionGaussWidth(0),
  68   fSDDKaonMPV(0),
  69   fSDDKaonLandauWidth(0),
  70   fSDDKaonGaussWidth(0),
  71   fSSDKaonMPV(0),
  72   fSSDKaonLandauWidth(0),
  73   fSSDKaonGaussWidth(0),
  74   fSDDProtMPV(0),
  75   fSDDProtLandauWidth(0),
  76   fSDDProtGaussWidth(0),
  77   fSSDProtMPV(0),
  78   fSSDProtLandauWidth(0),
  79   fSSDProtGaussWidth(0)
  80 {
  81   // standard constructor
  82   InitMC();
  83 }
  84 //______________________________________________________________________
  85 AliITSPidParams::~AliITSPidParams(){
  86   //
  87   if(fSDDPionMPV) delete fSDDPionMPV;
  88   if(fSDDPionLandauWidth) delete fSDDPionLandauWidth;
  89   if(fSDDPionGaussWidth) delete fSDDPionGaussWidth;
  90
  91   if(fSSDPionMPV) delete fSSDPionMPV;
  92   if(fSSDPionLandauWidth) delete fSSDPionLandauWidth;
  93   if(fSSDPionGaussWidth) delete fSSDPionGaussWidth;
  94
  95   if(fSDDKaonMPV) delete fSDDKaonMPV;
  96   if(fSDDKaonLandauWidth) delete fSDDKaonLandauWidth;
  97   if(fSDDKaonGaussWidth) delete fSDDKaonGaussWidth;
  98
  99   if(fSSDKaonMPV) delete fSSDKaonMPV;
 100   if(fSSDKaonLandauWidth) delete fSSDKaonLandauWidth;
 101   if(fSSDKaonGaussWidth) delete fSSDKaonGaussWidth;
 102
 103   if(fSDDProtMPV) delete fSDDProtMPV;
 104   if(fSDDProtLandauWidth) delete fSDDProtLandauWidth;
 105   if(fSDDProtGaussWidth) delete fSDDProtGaussWidth;
 106
 107   if(fSSDProtMPV) delete fSSDProtMPV;
 108   if(fSSDProtLandauWidth) delete fSSDProtLandauWidth;
 109   if(fSSDProtGaussWidth) delete fSSDProtGaussWidth;
 110 }
 111
 112 //______________________________________________________________________
 113 void AliITSPidParams::InitMC(){
 114   // initialize TFormulas to Monte Carlo values (=p-p simulations PYTHIA+GEANT)
 115   // parameter values from LHC10d1
 116
 117   // pions
 118   if(fSDDPionMPV) delete fSDDPionMPV;
 119   fSDDPionMPV=new TFormula("fSDDPionMPV","[0]/(x*x)*TMath::Log(x)+[1]/(x*x*x*x)*TMath::Log(x)+[2]*TMath::Log(x)+[3]");
 120   fSDDPionMPV->SetParameters(-0.690010, 0.002602, 1.185083, 78.454691);
 121
 122   if(fSDDPionLandauWidth) delete fSDDPionLandauWidth;
 123   fSDDPionLandauWidth=new TFormula("fSDDPionLandauWidth","[0]/(x*x)+[1]");
 124   fSDDPionLandauWidth->SetParameters(0.061606, 5.960376);
 125
 126   if(fSDDPionGaussWidth) delete fSDDPionGaussWidth;
 127   fSDDPionGaussWidth=new TFormula("fSDDPionGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
 128   fSDDPionGaussWidth->SetParameters(-0.065307, 7.896339);
 129
 130   if(fSSDPionMPV) delete fSSDPionMPV;
 131   fSSDPionMPV=new TFormula("fSSDPionMPV","[0]/(x*x)*TMath::Log(x)+[1]/(x*x*x*x)*TMath::Log(x)+[2]*TMath::Log(x)+[3]");
 132   fSSDPionMPV->SetParameters(-0.699466, 0.002429, 1.366895, 80.759188);
 133
 134   if(fSSDPionLandauWidth) delete fSSDPionLandauWidth;
 135   fSSDPionLandauWidth=new TFormula("fSSDPionLandauWidth","[0]/(x*x)+[1]");
 136   fSSDPionLandauWidth->SetParameters(0.066319, 5.889438);
 137
 138   if(fSSDPionGaussWidth) delete fSSDPionGaussWidth;
 139   fSSDPionGaussWidth=new TFormula("fSSDPionGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
 140   fSSDPionGaussWidth->SetParameters(-0.077798, 5.903887);
 141
 142   // kaons
 143   if(fSDDKaonMPV) delete fSDDKaonMPV;
 144   fSDDKaonMPV=new TFormula("fSDDKaonMPV","[0]/(x*x)+[1]/(x*x*x*x)*TMath::Log(x)+[2]");
 145   fSDDKaonMPV->SetParameters(15.924230, 0.085357, 73.528107);
 146
 147   if(fSDDKaonLandauWidth) delete fSDDKaonLandauWidth;
 148   fSDDKaonLandauWidth=new TFormula("fSDDKaonLandauWidth","[0]/(x*x)+[1]");
 149   fSDDKaonLandauWidth->SetParameters(1.121062, 5.925409);
 150
 151   if(fSDDKaonGaussWidth) delete fSDDKaonGaussWidth;
 152   fSDDKaonGaussWidth=new TFormula("fSDDKaonGaussWidth","[0]/(x*x)+[1]");
 153   fSDDKaonGaussWidth->SetParameters(2.010609, 5.973445);
 154
 155   if(fSSDKaonMPV) delete fSSDKaonMPV;
 156   fSSDKaonMPV=new TFormula("fSSDKaonMPV","[0]/(x*x)+[1]/(x*x*x*x)*TMath::Log(x)+[2]");
 157   fSSDKaonMPV->SetParameters(15.197250, 0.016714, 76.446132);
 158
 159   if(fSSDKaonLandauWidth) delete fSSDKaonLandauWidth;
 160   fSSDKaonLandauWidth=new TFormula("fSSDKaonLandauWidth","[0]/(x*x)+[1]");
 161   fSSDKaonLandauWidth->SetParameters(1.036749, 6.106413);
 162
 163   if(fSSDKaonGaussWidth) delete fSSDKaonGaussWidth;
 164   fSSDKaonGaussWidth=new TFormula("fSSDKaonGaussWidth","[0]/(x*x)+[1]");
 165   fSSDKaonGaussWidth->SetParameters(2.426498, 3.383779);
 166
 167   // protons
 168   if(fSDDProtMPV) delete fSDDProtMPV;
 169   fSDDProtMPV=new TFormula("fSDDProtMPV","[0]/(x*x)+[1]/(x*x*x*x)*TMath::Log(x)+[2]");
 170   fSDDProtMPV->SetParameters(56.888592, 1.115447, 75.416075);
 171
 172   if(fSDDProtLandauWidth) delete fSDDProtLandauWidth;
 173   fSDDProtLandauWidth=new TFormula("fSDDProtLandauWidth","[0]/(x*x)+[1]");
 174   fSDDProtLandauWidth->SetParameters(6.350805, 4.312568);
 175
 176   if(fSDDProtGaussWidth) delete fSDDProtGaussWidth;
 177   fSDDProtGaussWidth=new TFormula("fSDDProtGaussWidth","[0]/(x*x)+[1]");
 178   fSDDProtGaussWidth->SetParameters(6.556759, 5.953683);
 179
 180   if(fSSDProtMPV) delete fSSDProtMPV;
 181   fSSDProtMPV=new TFormula("fSSDProtMPV","[0]/(x*x)+[1]/(x*x*x*x)*TMath::Log(x)+[2]");
 182   fSSDProtMPV->SetParameters(57.385512, 0.884585, 76.138989);
 183
 184   if(fSSDProtLandauWidth) delete fSSDProtLandauWidth;
 185   fSSDProtLandauWidth=new TFormula("fSSDProtLandauWidth","[0]/(x*x)+[1]");
 186   fSSDProtLandauWidth->SetParameters(6.653282, 3.997930);
 187
 188   if(fSSDProtGaussWidth) delete fSSDProtGaussWidth;
 189   fSSDProtGaussWidth=new TFormula("fSSDProtGaussWidth","[0]/(x*x)+[1]");
 190   fSSDProtGaussWidth->SetParameters(8.203296, 1.491822);
 191
 192 }
 193 //_______________________________________________________________________
 194 Double_t AliITSPidParams::GetLandauGausNormPdgCode(Double_t dedx, Int_t pdgCode, Double_t mom, Int_t lay) const {
 195   // Computes Landau Gauss convolution for given particle specie and given momentum in a given ITS layer
 196   if(TMath::Abs(pdgCode)==211) return GetLandauGausNorm(dedx,AliPID::kPion,mom,lay);
 197   else if(TMath::Abs(pdgCode)==321) return GetLandauGausNorm(dedx,AliPID::kKaon,mom,lay);
 198   else if(TMath::Abs(pdgCode)==2212) return GetLandauGausNorm(dedx,AliPID::kProton,mom,lay);
 199   else return 0.;
 200 }
 201 //_______________________________________________________________________
 202   Double_t AliITSPidParams::GetLandauGausNorm(Double_t dedx, Int_t partType, Double_t mom, Int_t lay) const{
 203   // Computes Landau Gauss convolution for given particle specie and given momentum in a given ITS layer
 204
 205   Double_t par[3];
 206   Bool_t isSet=kFALSE;
 207   if(partType==AliPID::kPion){
 208     if(lay==3 || lay==4){
 209       par[0]=GetSDDPionLandauWidth(mom);
 210       par[1]=GetSDDPionMPV(mom);
 211       par[2]=GetSDDPionGaussWidth(mom);
 212       isSet=kTRUE;
 213     }
 214     else if(lay==5 || lay==6){
 215       par[0]=GetSSDPionLandauWidth(mom);
 216       par[1]=GetSSDPionMPV(mom);
 217       par[2]=GetSSDPionGaussWidth(mom);
 218       isSet=kTRUE;
 219     }
 220   }else if(partType==AliPID::kKaon){
 221     if(lay==3 || lay==4){
 222       par[0]=GetSDDKaonLandauWidth(mom);
 223       par[1]=GetSDDKaonMPV(mom);
 224       par[2]=GetSDDKaonGaussWidth(mom);
 225       isSet=kTRUE;
 226     }
 227     else if(lay==5 || lay==6){
 228       par[0]=GetSSDKaonLandauWidth(mom);
 229       par[1]=GetSSDKaonMPV(mom);
 230       par[2]=GetSSDKaonGaussWidth(mom);
 231       isSet=kTRUE;
 232     }
 233   }else if(partType==AliPID::kProton){
 234     if(lay==3 || lay==4){
 235       par[0]=GetSDDProtLandauWidth(mom);
 236       par[1]=GetSDDProtMPV(mom);
 237       par[2]=GetSDDProtGaussWidth(mom);
 238       isSet=kTRUE;
 239     }
 240     else if(lay==5 || lay==6){
 241       par[0]=GetSSDProtLandauWidth(mom);
 242       par[1]=GetSSDProtMPV(mom);
 243       par[2]=GetSSDProtGaussWidth(mom);
 244       isSet=kTRUE;
 245     }
 246   }
 247   if(!isSet) return 0.;
 248   // Numeric constants
 249   Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
 250   Double_t mpshift  = -0.22278298;       // Landau maximum location
 251   // Control constants
 252   Double_t np = 100.0;      // number of convolution steps
 253   Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas
 254   // Variables
 255   Double_t xx;
 256   Double_t mpc;
 257   Double_t fland;
 258   Double_t sum = 0.0;
 259   Double_t xlow,xupp;
 260   Double_t step = 0.;
 261   Double_t i;
 262
 263   // MP shift correction
 264   mpc = par[1] - mpshift * par[0];
 265   // Range of convolution integral
 266   xlow = dedx - sc * par[2];
 267   xupp = dedx + sc * par[2];
 268   if(np!=0) step = (xupp-xlow) / np;
 269
 270   // Convolution integral of Landau and Gaussian by sum
 271   for(i=1.0; i<=np/2; i++) {
 272     xx = xlow + (i-.5) * step;
 273
 274     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
 275     sum += fland * TMath::Gaus(dedx,xx,par[2]);
 276
 277     xx = xupp - (i-.5) * step;
 278     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
 279     sum += fland * TMath::Gaus(dedx,xx,par[2]);
 280   }
 281
 282   return (step * sum * invsq2pi / par[2]);
 283 }
 284