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