TRD/AliTRDCalPIDLQ.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 "AliTRDCalPIDLQ.h"
  25 #include <TH1F.h>
  26 #include <TFile.h>
  27
  28 ClassImp(AliTRDCalPIDLQ)
  29
  30 Char_t* AliTRDCalPIDLQ::fpartName[AliPID::kSPECIES] = {"electron", "muon", "pion", "kaon", "proton"};
  31
  32 //_________________________________________________________________________
  33 AliTRDCalPIDLQ::AliTRDCalPIDLQ(): TNamed()
  34 {
  35   //
  36   //  The Default constructor
  37   //
  38
  39   Init();
  40 }
  41
  42 //_________________________________________________________________________
  43 AliTRDCalPIDLQ::AliTRDCalPIDLQ(const Text_t *name, const Text_t *title) : TNamed(name, title)
  44 {
  45   //
  46   //  The main constructor
  47   //
  48
  49   Init();
  50 }
  51
  52 //_____________________________________________________________________________
  53 AliTRDCalPIDLQ::AliTRDCalPIDLQ(const AliTRDCalPIDLQ &c) : TNamed(c)
  54 {
  55   //
  56   // copy constructor
  57   //
  58
  59   Init();
  60
  61   ((AliTRDCalPIDLQ &) c).Copy(*this);
  62 }
  63
  64 //_________________________________________________________________________
  65 AliTRDCalPIDLQ::~AliTRDCalPIDLQ()
  66 {
  67   //
  68   // Destructor
  69   //
  70
  71   CleanUp();
  72 }
  73
  74 //_________________________________________________________________________
  75 void AliTRDCalPIDLQ::CleanUp()
  76 {
  77   if (fHistdEdx)
  78   {
  79     delete fHistdEdx;
  80     fHistdEdx = 0;
  81   }
  82
  83   if (fHistTimeBin)
  84   {
  85     delete fHistTimeBin;
  86     fHistTimeBin = 0;
  87   }
  88
  89   if (fTrackMomentum)
  90   {
  91     delete[] fTrackMomentum;
  92     fTrackMomentum = 0;
  93   }
  94 }
  95
  96 //_____________________________________________________________________________
  97 AliTRDCalPIDLQ &AliTRDCalPIDLQ::operator=(const AliTRDCalPIDLQ &c)
  98 {
  99   //
 100   // Assignment operator
 101   //
 102
 103   if (this != &c) ((AliTRDCalPIDLQ &) c).Copy(*this);
 104   return *this;
 105 }
 106
 107 //_____________________________________________________________________________
 108 void AliTRDCalPIDLQ::Copy(TObject &c) const
 109 {
 110   //
 111   // Copy function
 112   //
 113
 114   AliTRDCalPIDLQ& target = (AliTRDCalPIDLQ &) c;
 115
 116   target.CleanUp();
 117
 118   target.fNMom = fNMom;
 119
 120   target.fTrackMomentum = new Double_t[fNMom];
 121   for (Int_t i=0; i<fNMom; ++i)
 122     target.fTrackMomentum[i] = fTrackMomentum[i];
 123
 124   target.fMeanChargeRatio = fMeanChargeRatio;
 125
 126   target.fNbins = fNbins;
 127   target.fBinSize = fBinSize;
 128
 129   if (fHistdEdx)
 130     target.fHistdEdx = (TObjArray*) fHistdEdx->Clone();
 131
 132   if (fHistTimeBin)
 133     target.fHistTimeBin = (TObjArray*) fHistTimeBin->Clone();
 134
 135   TObject::Copy(c);
 136 }
 137
 138 //_________________________________________________________________________
 139 void AliTRDCalPIDLQ::Init()
 140 {
 141   fNMom = 11;
 142
 143   fTrackMomentum = new Double_t[fNMom];
 144   Double_t trackMomentum[11] = {0.6, 0.8, 1, 1.5, 2, 3, 4, 5, 6, 8, 10};
 145   for (Int_t imom = 0; imom < 11; imom++)
 146     fTrackMomentum[imom] = trackMomentum[imom];
 147
 148   fHistdEdx = new TObjArray(AliPID::kSPECIES * fNMom);
 149   fHistdEdx->SetOwner();
 150   fHistTimeBin = new TObjArray(AliPID::kSPECIES * fNMom);
 151   fHistTimeBin->SetOwner();
 152
 153   // ADC Gain normalization
 154   fMeanChargeRatio=1.0;
 155
 156   // Number of bins and bin size
 157   fNbins = 0;
 158   fBinSize = 0.;
 159 }
 160
 161 //_________________________________________________________________________
 162 Bool_t AliTRDCalPIDLQ::ReadData(Char_t *responseFile)
 163 {
 164   //
 165   // Read the TRD dEdx histograms.
 166   //
 167   // Read histogram Root file
 168   TFile *histFile = new TFile(responseFile, "READ");
 169   if (!histFile || !histFile->IsOpen()) {
 170     Error("AliTRDCalPIDLQ", "opening TRD histgram file %s failed", responseFile);
 171     return kFALSE;
 172   }
 173   gROOT->cd();
 174
 175   // Read histograms
 176   Char_t text[200];
 177   for (Int_t particle = 0; particle < AliPID::kSPECIES; ++particle)
 178   {
 179     Char_t* particleKey = "";
 180     switch (particle)
 181     {
 182       case AliPID::kElectron: particleKey = "EL"; break;
 183       case AliPID::kPion: particleKey = "PI"; break;
 184       case AliPID::kMuon: particleKey = "MU"; break;
 185       case AliPID::kKaon: particleKey = "KA"; break;
 186       case AliPID::kProton: particleKey = "PR"; break;
 187     }
 188
 189     for (Int_t imom = 0; imom < fNMom; imom++)
 190     {
 191       sprintf(text, "h1dEdx%s%01d", particleKey, imom+1);
 192       TH1F* hist = (TH1F*)histFile->Get(text)->Clone();
 193       hist->Scale(1.0/hist->Integral());
 194       fHistdEdx->AddAt(hist, GetHistID(particle, imom));
 195
 196       if (particle == AliPID::kElectron || particle == AliPID::kPion)
 197       {
 198         sprintf(text,"h1MaxTimBin%s%01d", particleKey, imom+1);
 199         TH1F* hist = (TH1F*)histFile->Get(text)->Clone();
 200         hist->Scale(1.0/hist->Integral());
 201         fHistTimeBin->AddAt(hist, GetHistID(particle,imom));
 202       }
 203     }
 204   }
 205
 206   histFile->Close();
 207   delete histFile;
 208
 209   // Number of bins and bin size
 210   TH1F* hist = (TH1F*) fHistdEdx->At(GetHistID(AliPID::kPion, 1));
 211   fNbins   = hist->GetNbinsX();
 212   fBinSize = hist->GetBinWidth(1);
 213
 214   return kTRUE;
 215 }
 216
 217 //_________________________________________________________________________
 218 Double_t  AliTRDCalPIDLQ::GetMean(Int_t k, Int_t ip) const
 219 {
 220   //
 221   // Gets mean of de/dx dist. of e
 222   printf("Mean for particle = %s and momentum = %.2f is:\n", fpartName[k], fTrackMomentum[ip]);
 223   if (k < 0 || k > AliPID::kSPECIES)
 224     return 0;
 225
 226   return ((TH1F*) fHistdEdx->At(GetHistID(k,ip)))->GetMean();
 227 }
 228
 229 //_________________________________________________________________________
 230 Double_t  AliTRDCalPIDLQ::GetNormalization(Int_t k, Int_t ip) const
 231 {
 232   //
 233   // Gets Normalization of de/dx dist. of e
 234
 235   printf("Normalization for particle = %s and momentum = %.2f is:\n",fpartName[k], fTrackMomentum[ip]);
 236   if (k < 0 || k > AliPID::kSPECIES)
 237     return 0;
 238
 239   return ((TH1F*) fHistdEdx->At(GetHistID(k,ip)))->Integral();
 240 }
 241
 242 //_________________________________________________________________________
 243 TH1F* AliTRDCalPIDLQ::GetHistogram(Int_t k, Int_t ip) const
 244 {
 245   //
 246   //
 247   printf("Histogram for particle = %s and momentum = %.2f is:\n", fpartName[k], fTrackMomentum[ip]);
 248   if (k < 0 || k > AliPID::kSPECIES)
 249     return 0;
 250
 251   return (TH1F*) fHistdEdx->At(GetHistID(k,ip));
 252 }
 253
 254 //_________________________________________________________________________
 255 Double_t AliTRDCalPIDLQ::GetProbability(Int_t k, Double_t mom, Double_t dedx1) const
 256 {
 257   //
 258   // Gets the Probability of having dedx at a given momentum (mom)
 259   // and particle type k (0 for e) and (2 for pi)
 260   // from the precalculated de/dx distributions
 261
 262   Double_t dedx = dedx1/fMeanChargeRatio;
 263   Int_t iEnBin= ((Int_t) (dedx/fBinSize+1));
 264   if(iEnBin > fNbins) iEnBin = fNbins;
 265
 266   if (k < 0 || k > AliPID::kSPECIES)
 267     return 1;
 268
 269   TH1F* hist1 = 0;
 270   TH1F* hist2 = 0;
 271   Double_t mom1 = 0;
 272   Double_t mom2 = 0;
 273
 274   // Lower limit
 275   if (mom<=fTrackMomentum[0])
 276   {
 277     hist1 = (TH1F*) fHistdEdx->At(GetHistID(k,1));
 278     hist2 = (TH1F*) fHistdEdx->At(GetHistID(k,0));
 279     mom1 = fTrackMomentum[1];
 280     mom2 = fTrackMomentum[0];
 281   }
 282
 283   // Upper Limit
 284   if(mom>=fTrackMomentum[fNMom-1])
 285   {
 286     hist1 = (TH1F*) fHistdEdx->At(GetHistID(k,fNMom-1));
 287     hist2 = (TH1F*) fHistdEdx->At(GetHistID(k,fNMom-2));
 288     mom1 = fTrackMomentum[fNMom-1];
 289     mom2 = fTrackMomentum[fNMom-2];
 290   }
 291
 292   // In the range
 293   for (Int_t ip=1; ip<fNMom; ip++)
 294   {
 295     if ((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip]))
 296     {
 297       hist1 = (TH1F*) fHistdEdx->At(GetHistID(k,ip));
 298       hist2 = (TH1F*) fHistdEdx->At(GetHistID(k,ip-1));
 299       mom1 = fTrackMomentum[ip];
 300       mom2 = fTrackMomentum[ip-1];
 301     }
 302   }
 303
 304   Double_t slop = (hist1->GetBinContent(iEnBin) - hist2->GetBinContent(iEnBin)) / (mom1 - mom2);
 305   return hist2->GetBinContent(iEnBin) + slop * (mom - mom2);
 306 }
 307
 308 //_________________________________________________________________________
 309 Double_t AliTRDCalPIDLQ::GetProbabilityT(Int_t k, Double_t mom, Int_t timbin) const
 310 {
 311   //
 312   // Gets the Probability of having timbin at a given momentum (mom)
 313   // and particle type k (0 for e) and (2 for pi)
 314   // from the precalculated timbin distributions
 315
 316   if (timbin<=0)
 317     return 0.;
 318   Int_t iTBin=timbin+1;
 319
 320   // everything which is not electron counts as pion for time bin
 321   if (k != AliPID::kElectron)
 322     k = AliPID::kPion;
 323
 324   if (mom<=fTrackMomentum[0])
 325     return ((TH1F*) fHistTimeBin->At(GetHistID(k,0)))->GetBinContent(iTBin);
 326
 327   if (mom>=fTrackMomentum[fNMom-1])
 328     return ((TH1F*) fHistTimeBin->At(GetHistID(k,fNMom-1)))->GetBinContent(iTBin);
 329
 330   for (Int_t ip=1; ip<fNMom; ip++)
 331   {
 332     if ((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip]))
 333     {
 334       Double_t slop=(((TH1F*) fHistTimeBin->At(GetHistID(k,ip)))->GetBinContent(iTBin) - ((TH1F*) fHistTimeBin->At(GetHistID(k,ip-1)))->GetBinContent(iTBin)) / (fTrackMomentum[ip] - fTrackMomentum[ip-1]);
 335       // Linear Interpolation
 336       return ((TH1F*) fHistTimeBin->At(GetHistID(k,ip-1)))->GetBinContent(iTBin) + slop * (mom - fTrackMomentum[ip-1]);
 337     }
 338   }
 339
 340   return -1;
 341 }