TRD/AliTRDprobdist.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 //-----------------------------------------------------------------
  17 // Class for dE/dx and Time Bin of Max. Cluster for Electrons and
  18 // pions in TRD.
  19 // It is instantiated in class AliTRDpidESD for particle identification
  20 // in TRD
  21 // Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de>
  22 //-----------------------------------------------------------------
  23
  24 #include "AliTRDprobdist.h"
  25
  26 ClassImp(AliTRDprobdist)
  27
  28 //_________________________________________________________________________
  29 AliTRDprobdist::AliTRDprobdist(Char_t *responseFile)
  30 {
  31   //
  32   //  The main constructor
  33   //
  34   Char_t *partName[5] = {"electron", "muon", "pion", "kaon", "proton"};
  35   for(Int_t i=0; i<5; i++) {
  36     fpartName[i] = partName[i];
  37   }
  38   // ADC Gain normalization
  39   fADCNorm=1.0;
  40   fNMom=kNMom;
  41   // Set track momenta for which response functions are available
  42   Double_t trackMomentum[kNMom]= {0.6, 0.8, 1, 1.5, 2, 3, 4, 5, 6, 8, 10};
  43   for(Int_t imom=0; imom<fNMom; imom++) {
  44     fTrackMomentum[imom] = trackMomentum[imom];
  45   }
  46   // Read the histogram file
  47   ReadData(responseFile);
  48   //
  49 }
  50
  51 Bool_t AliTRDprobdist::ReadData(Char_t *responseFile)
  52 {
  53   //
  54   // Read the TRD dEdx histograms.
  55   //
  56   // Read histogram Root file
  57   TFile *histFile = new TFile(responseFile);
  58   if (!histFile || !histFile->IsOpen()) {
  59     Error("AliTRDprobdist", "opening TRD histgram file %s failed", responseFile);
  60     return kFALSE;
  61   }
  62
  63   // Read histograms
  64   Char_t text[200];
  65   for (Int_t imom = 0; imom < kNMom; imom++) {
  66     sprintf(text,"h1dEdxEL%01d",imom+1);
  67     fh1dEdxEL[imom] = (TH1F*)histFile->Get(text);
  68     fh1dEdxEL[imom]->Scale(1.0/fh1dEdxEL[imom]->Integral());
  69
  70     sprintf(text,"h1dEdxPI%01d",imom+1);
  71     fh1dEdxPI[imom] = (TH1F*)histFile->Get(text);
  72     fh1dEdxPI[imom]->Scale(1.0/fh1dEdxPI[imom]->Integral());
  73
  74     sprintf(text,"h1dEdxMU%01d",imom+1);
  75     fh1dEdxMU[imom] = (TH1F*)histFile->Get(text);
  76     fh1dEdxMU[imom]->Scale(1.0/fh1dEdxMU[imom]->Integral());
  77
  78     sprintf(text,"h1dEdxKA%01d",imom+1);
  79     fh1dEdxKA[imom] = (TH1F*)histFile->Get(text);
  80     fh1dEdxKA[imom]->Scale(1.0/fh1dEdxKA[imom]->Integral());
  81
  82     sprintf(text,"h1dEdxPR%01d",imom+1);
  83     fh1dEdxPR[imom] = (TH1F*)histFile->Get(text);
  84     fh1dEdxPR[imom]->Scale(1.0/fh1dEdxPR[imom]->Integral());
  85
  86     sprintf(text,"h1MaxTimBinEL%01d",imom+1);
  87     fh1MaxTimBinEL[imom] = (TH1F*)histFile->Get(text);
  88     fh1MaxTimBinEL[imom]->Scale(1.0/fh1MaxTimBinEL[imom]->Integral());
  89
  90     sprintf(text,"h1MaxTimBinPI%01d",imom+1);
  91     fh1MaxTimBinPI[imom] = (TH1F*)histFile->Get(text);
  92     fh1MaxTimBinPI[imom]->Scale(1.0/fh1MaxTimBinPI[imom]->Integral());
  93   }
  94   // Number of bins and bin size
  95   fNbins = fh1dEdxPI[1]->GetNbinsX();
  96   fBinSize = fh1dEdxPI[1]->GetBinWidth(1);
  97   return kTRUE;
  98 }
  99
 100 AliTRDprobdist::AliTRDprobdist(const AliTRDprobdist& pd):TNamed()
 101 {
 102   //
 103   // Copy constructor.
 104   //
 105   for(Int_t i=0; i<5; i++) {
 106     fpartName[i] = pd.fpartName[i];
 107   }
 108   // ADC Gain normalization
 109   fADCNorm=pd.fADCNorm;
 110   fNMom=pd.fNMom;
 111   // Set track momenta for which response functions are available
 112   for(Int_t imom=0; imom<fNMom; imom++) {
 113     fTrackMomentum[imom] = pd.fTrackMomentum[imom];
 114   }
 115
 116   fNbins=pd.fNbins;
 117   fBinSize=pd.fBinSize;
 118
 119   for (Int_t imom = 0; imom < kNMom; imom++) {
 120     fh1dEdxEL[imom] = pd.fh1dEdxEL[imom];
 121     fh1dEdxPI[imom] = pd.fh1dEdxPI[imom];
 122     fh1dEdxMU[imom] = pd.fh1dEdxMU[imom];
 123     fh1dEdxKA[imom] = pd.fh1dEdxKA[imom];
 124     fh1dEdxPR[imom] = pd.fh1dEdxPR[imom];
 125     fh1MaxTimBinEL[imom] = pd.fh1MaxTimBinEL[imom];
 126     fh1MaxTimBinPI[imom] = pd.fh1MaxTimBinPI[imom];
 127   }
 128 }
 129
 130
 131 AliTRDprobdist::~AliTRDprobdist()
 132 {
 133   // Destructor
 134
 135 }
 136
 137
 138 //_________________________________________________________________________
 139 Double_t  AliTRDprobdist::GetMean(Int_t k, Int_t ip) const
 140 {
 141   //
 142   // Gets mean of de/dx dist. of e
 143   printf("Mean for particle = %s and momentum = %.2f is:\n", fpartName[k], fTrackMomentum[ip]);
 144   if(k==0) return fh1dEdxEL[ip]->GetMean();
 145   if(k==1) return fh1dEdxMU[ip]->GetMean();
 146   if(k==2) return fh1dEdxPI[ip]->GetMean();
 147   if(k==3) return fh1dEdxKA[ip]->GetMean();
 148   if(k==4) return fh1dEdxPR[ip]->GetMean();
 149   return fh1dEdxPR[ip]->GetMean();
 150 }
 151
 152 //_________________________________________________________________________
 153 Double_t  AliTRDprobdist::GetNormalization(Int_t k, Int_t ip) const
 154 {
 155   //
 156   // Gets Normalization of de/dx dist. of e
 157
 158   printf("Normalization for particle = %s and momentum = %.2f is:\n",fpartName[k], fTrackMomentum[ip]);
 159   if(k==0) return fh1dEdxEL[ip]->Integral();
 160   if(k==1) return fh1dEdxMU[ip]->Integral();
 161   if(k==2) return fh1dEdxPI[ip]->Integral();
 162   if(k==3) return fh1dEdxKA[ip]->Integral();
 163   if(k==4) return fh1dEdxPR[ip]->Integral();
 164   return fh1dEdxPR[ip]->Integral();
 165 }
 166
 167 TH1F* AliTRDprobdist::GetHistogram(Int_t k, Int_t ip) const
 168 {
 169   //
 170   //
 171   printf("Histogram for particle = %s and momentum = %.2f is:\n", fpartName[k], fTrackMomentum[ip]);
 172   if(k==0) return fh1dEdxEL[ip];
 173   if(k==1) return fh1dEdxMU[ip];
 174   if(k==2) return fh1dEdxPI[ip];
 175   if(k==3) return fh1dEdxKA[ip];
 176   if(k==4) return fh1dEdxPR[ip];
 177   return fh1dEdxPR[ip];
 178 }
 179
 180 //_________________________________________________________________________
 181 Double_t AliTRDprobdist::GetProbability(Int_t k, Double_t mom, Double_t dedx1) const
 182 {
 183   //
 184   // Gets the Probability of having dedx at a given momentum (mom)
 185   // and particle type k (0 for e) and (2 for pi)
 186   // from the precalculated de/dx distributions
 187   Double_t probability = 1.0;
 188   Double_t dedx = dedx1/fADCNorm;
 189   Int_t iEnBin= ((Int_t) (dedx/fBinSize+1));
 190   if(iEnBin > fNbins) iEnBin = fNbins;
 191
 192   ////Electron//////////////////////////
 193   if(k==0){    // electron
 194     Double_t slop;
 195     // Lower limit
 196     if(mom<=fTrackMomentum[0]) {
 197       slop=(fh1dEdxEL[1]->GetBinContent(iEnBin)-fh1dEdxEL[0]->GetBinContent(iEnBin))/(fTrackMomentum[1] - fTrackMomentum[0]);
 198       probability= fh1dEdxEL[0]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[0]);
 199       return probability;
 200     }
 201     // Upper Limit
 202     if(mom>=fTrackMomentum[fNMom-1]) {
 203       slop=(fh1dEdxEL[fNMom-1]->GetBinContent(iEnBin)-fh1dEdxEL[fNMom-2]->GetBinContent(iEnBin))/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]);
 204       probability= fh1dEdxEL[fNMom-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[fNMom-1]);
 205       return probability;
 206     }
 207     // In the range
 208     for(Int_t ip=1; ip<fNMom; ip++)
 209       if((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
 210         slop=(fh1dEdxEL[ip]->GetBinContent(iEnBin)-fh1dEdxEL[ip-1]->GetBinContent(iEnBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]);
 211         // Linear Interpolation
 212         probability= fh1dEdxEL[ip-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[ip-1]);
 213         return probability;
 214       }
 215   }
 216
 217   ////Pion//////////////////////////
 218   if(k==2){    // Pion
 219     Double_t slop;
 220     // Lower limit
 221     if(mom<=fTrackMomentum[0]) {
 222       slop=(fh1dEdxPI[1]->GetBinContent(iEnBin)-fh1dEdxPI[0]->GetBinContent(iEnBin))/(fTrackMomentum[1] - fTrackMomentum[0]);
 223       probability= fh1dEdxPI[0]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[0]);
 224       return probability;
 225     }
 226     // Upper Limit
 227     if(mom>=fTrackMomentum[fNMom-1]) {
 228       slop=(fh1dEdxPI[fNMom-1]->GetBinContent(iEnBin)-fh1dEdxPI[fNMom-2]->GetBinContent(iEnBin))/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]);
 229       probability= fh1dEdxPI[fNMom-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[fNMom-1]);
 230       return probability;
 231     }
 232     // In the range
 233     for(Int_t ip=1; ip<fNMom; ip++)
 234       if((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
 235         slop=(fh1dEdxPI[ip]->GetBinContent(iEnBin)-fh1dEdxPI[ip-1]->GetBinContent(iEnBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]);
 236         // Linear Interpolation
 237         probability= fh1dEdxPI[ip-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[ip-1]);
 238         return probability;
 239       }
 240   }
 241
 242   ////Muon//////////////////////////
 243   if(k==1){    // Muon
 244     Double_t slop;
 245     // Lower limit
 246     if(mom<=fTrackMomentum[0]) {
 247       slop=(fh1dEdxMU[1]->GetBinContent(iEnBin)-fh1dEdxMU[0]->GetBinContent(iEnBin))/(fTrackMomentum[1] - fTrackMomentum[0]);
 248       probability= fh1dEdxMU[0]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[0]);
 249       return probability;
 250     }
 251     // Upper Limit
 252     if(mom>=fTrackMomentum[fNMom-1]) {
 253       slop=(fh1dEdxMU[fNMom-1]->GetBinContent(iEnBin)-fh1dEdxMU[fNMom-2]->GetBinContent(iEnBin))/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]);
 254       probability= fh1dEdxMU[fNMom-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[fNMom-1]);
 255       return probability;
 256     }
 257     // In the range
 258     for(Int_t ip=1; ip<fNMom; ip++)
 259       if((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
 260         slop=(fh1dEdxMU[ip]->GetBinContent(iEnBin)-fh1dEdxMU[ip-1]->GetBinContent(iEnBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]);
 261         // Linear Interpolation
 262         probability= fh1dEdxMU[ip-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[ip-1]);
 263         return probability;
 264       }
 265   }
 266
 267   ////Kaon//////////////////////////
 268   if(k==3){    // Kaon
 269     Double_t slop;
 270     // Lower limit
 271     if(mom<=fTrackMomentum[0]) {
 272       slop=(fh1dEdxKA[1]->GetBinContent(iEnBin)-fh1dEdxKA[0]->GetBinContent(iEnBin))/(fTrackMomentum[1] - fTrackMomentum[0]);
 273       probability= fh1dEdxKA[0]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[0]);
 274       return probability;
 275     }
 276     // Upper Limit
 277     if(mom>=fTrackMomentum[fNMom-1]) {
 278       slop=(fh1dEdxKA[fNMom-1]->GetBinContent(iEnBin)-fh1dEdxKA[fNMom-2]->GetBinContent(iEnBin))/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]);
 279       probability= fh1dEdxKA[fNMom-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[fNMom-1]);
 280       return probability;
 281     }
 282     // In the range
 283     for(Int_t ip=1; ip<fNMom; ip++)
 284       if((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
 285         slop=(fh1dEdxKA[ip]->GetBinContent(iEnBin)-fh1dEdxKA[ip-1]->GetBinContent(iEnBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]);
 286         // Linear Interpolation
 287         probability= fh1dEdxKA[ip-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[ip-1]);
 288         return probability;
 289       }
 290   }
 291
 292   ////Proton//////////////////////////
 293   if(k==4){    // Proton
 294     Double_t slop;
 295     // Lower limit
 296     if(mom<=fTrackMomentum[0]) {
 297       slop=(fh1dEdxPR[1]->GetBinContent(iEnBin)-fh1dEdxPR[0]->GetBinContent(iEnBin))/(fTrackMomentum[1] - fTrackMomentum[0]);
 298       probability= fh1dEdxPR[0]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[0]);
 299       return probability;
 300     }
 301     // Upper Limit
 302     if(mom>=fTrackMomentum[fNMom-1]) {
 303       slop=(fh1dEdxPR[fNMom-1]->GetBinContent(iEnBin)-fh1dEdxPR[fNMom-2]->GetBinContent(iEnBin))/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]);
 304       probability= fh1dEdxPR[fNMom-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[fNMom-1]);
 305       return probability;
 306     }
 307     // In the range
 308     for(Int_t ip=1; ip<fNMom; ip++)
 309       if((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
 310         slop=(fh1dEdxPR[ip]->GetBinContent(iEnBin)-fh1dEdxPR[ip-1]->GetBinContent(iEnBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]);
 311         // Linear Interpolation
 312         probability= fh1dEdxPR[ip-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[ip-1]);
 313         return probability;
 314       }
 315   }
 316
 317   return probability;
 318 }
 319
 320
 321 //_________________________________________________________________________
 322 Double_t AliTRDprobdist::GetProbabilityT(Int_t k, Double_t mom, Int_t timbin) const
 323 {
 324   //
 325   // Gets the Probability of having timbin at a given momentum (mom)
 326   // and particle type k (0 for e) and (2 for pi)
 327   // from the precalculated timbin distributions
 328   Double_t probabilityT = 1.0;
 329   if(timbin<=0) return 0.;
 330   Int_t iTBin=timbin+1;
 331
 332   if(k==0){    // electron
 333     if(mom<=fTrackMomentum[0]) probabilityT = fh1MaxTimBinEL[0]->GetBinContent(iTBin);
 334     if(mom>=fTrackMomentum[fNMom-1]) probabilityT = fh1MaxTimBinEL[fNMom-1]->GetBinContent(iTBin);
 335   }
 336   if(k==1||k==2||k==3||k==4){    // pion
 337     if(mom<=fTrackMomentum[0]) probabilityT = fh1MaxTimBinPI[0]->GetBinContent(iTBin);
 338     if(mom>=fTrackMomentum[fNMom-1]) probabilityT = fh1MaxTimBinPI[fNMom-1]->GetBinContent(iTBin);
 339   }
 340
 341   if(k==0)    // electron
 342   for(Int_t ip=1; ip<fNMom; ip++)
 343   if((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
 344     Double_t slop=(fh1MaxTimBinEL[ip]->GetBinContent(iTBin)-fh1MaxTimBinEL[ip-1]->GetBinContent(iTBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]);
 345     // Linear Interpolation
 346     probabilityT= fh1MaxTimBinEL[ip-1]->GetBinContent(iTBin) + slop*(mom-fTrackMomentum[ip-1]);
 347     return probabilityT;
 348   }
 349
 350   if(k==1||k==2||k==3||k==4)   // pion and other particles
 351   for(Int_t ip=1; ip<fNMom; ip++)
 352   if((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
 353     Double_t slop=(fh1MaxTimBinPI[ip]->GetBinContent(iTBin)-fh1MaxTimBinPI[ip-1]->GetBinContent(iTBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]);
 354     // Linear Interpolation
 355     probabilityT= fh1MaxTimBinPI[ip-1]->GetBinContent(iTBin) + slop*(mom-fTrackMomentum[ip-1]);
 356     return probabilityT;
 357   }
 358   return probabilityT;
 359 }
 360
 361