TRD/Cal/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 /* $Id$ */
  17
  18 //////////////////////////////////////////////////////////////////////////
  19 //                                                                      //
  20 // Container for the distributions of dE/dx and the time bin of the     //
  21 // max. cluster for electrons and pions                                 //
  22 //                                                                      //
  23 // Authors:                                                             //
  24 //   Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de>               //
  25 //   Alex Bercuci <a.bercuci@gsi.de>                                    //
  26 //                                                                      //
  27 //////////////////////////////////////////////////////////////////////////
  28
  29 #include <TH2F.h>
  30 #include <TFile.h>
  31 #include <TROOT.h>
  32
  33 #include "AliLog.h"
  34 #include "AliPID.h"
  35
  36 #include "../../STAT/TKDPDF.h"
  37 #include "AliTRDCalPIDLQ.h"
  38 #include "AliTRDcalibDB.h"
  39
  40 ClassImp(AliTRDCalPIDLQ)
  41
  42 Float_t AliTRDCalPIDLQ::fgTrackSegLength[kNLength] = { 3.7, 3.9, 4.2, 5.0 };
  43
  44 //_________________________________________________________________________
  45 AliTRDCalPIDLQ::AliTRDCalPIDLQ()
  46   :AliTRDCalPID("pid", "LQ PID references for TRD")
  47 {
  48   //
  49   //  The Default constructor
  50   //
  51
  52   //Init();
  53
  54 }
  55
  56 //_________________________________________________________________________
  57 AliTRDCalPIDLQ::AliTRDCalPIDLQ(const Text_t *name, const Text_t *title)
  58   :AliTRDCalPID(name,title)
  59 {
  60   //
  61   //  The main constructor
  62   //
  63
  64   Init();
  65
  66 }
  67
  68 //_________________________________________________________________________
  69 AliTRDCalPIDLQ::~AliTRDCalPIDLQ()
  70 {
  71   //
  72   // Destructor
  73   //
  74
  75 }
  76
  77 //_________________________________________________________________________
  78 Bool_t AliTRDCalPIDLQ::LoadPDF(TDirectoryFile *d)
  79 {
  80   // Read histograms
  81   TKDPDF *pdf = 0x0;
  82   for (Int_t is=0; is < AliPID::kSPECIES; is++){
  83     for (Int_t ip = 0; ip < kNMom; ip++){
  84       if(!(pdf = (TKDPDF*)d->Get(Form("%s[%d]", AliPID::ParticleShortName(is), ip)))){
  85         AliWarning(Form("Reference for %s[%d] missing.", AliPID::ParticleShortName(is), ip));
  86         continue;
  87       }
  88       fModel->AddAt(pdf->Clone(), GetModelID(ip, is, 0));
  89     }
  90   }
  91   return kTRUE;
  92 }
  93   //
  94 //_________________________________________________________________________
  95 Bool_t AliTRDCalPIDLQ::LoadReferences(Char_t *refFile)
  96 {
  97   //
  98   // Read the TRD dEdx histograms.
  99   //
 100
 101   Int_t nTimeBins = 22;
 102   // Get number of time bins from CDB
 103   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
 104   if(!calibration){
 105     AliWarning(Form("No AliTRDcalibDB available. Using %d time bins.", nTimeBins));
 106   }else{
 107     if(calibration->GetRun() > -1) nTimeBins = calibration->GetNumberOfTimeBins();
 108     else AliWarning(Form("Run number not set. Using %d time bins.", nTimeBins));
 109   }
 110
 111
 112   // Read histogram Root file
 113   TFile *histFile = TFile::Open(refFile, "READ");
 114   if (!histFile || !histFile->IsOpen()) {
 115     AliError(Form("Opening TRD histgram file %s failed", refFile));
 116     return kFALSE;
 117   }
 118   gROOT->cd();
 119
 120   // Read histograms
 121   for (Int_t iparticle = 0; iparticle < AliPID::kSPECIES; iparticle++){
 122     for (Int_t imom = 0; imom < kNMom; imom++){
 123       TH2D* hist = (TH2D*)histFile->Get(Form("h2dEdx%s%d", fPartSymb[iparticle], imom/*, ilength*/))->Clone();
 124       hist->Scale(1./hist->Integral());
 125       fModel->AddAt(hist, GetModelID(imom, iparticle, 0));
 126
 127 //                      if (iparticle != AliPID::kElectron && iparticle != AliPID::kPion) continue;
 128 //
 129 //                      TH1F* ht = (TH1F*)histFile->Get(Form("h1MaxTB%s%02d", fPartSymb[iparticle], imom))->Clone();
 130 //                      if(ht->GetNbinsX() != nTimeBins) AliWarning(Form("The number of time bins %d defined in h1MaxTB%s%02d differs from calibration value of %d. This may lead to erroneous results.", ht->GetNbinsX(), fPartSymb[iparticle], imom, nTimeBins));
 131 //                      ht->Scale(1./ht->Integral());
 132 //                      fHistTimeBin->AddAt(ht, ((iparticle==AliPID::kElectron)?0:1)*kNMom + imom);
 133     }
 134   }
 135
 136   histFile->Close();
 137   delete histFile;
 138
 139   // Number of bins and bin size
 140   //TH1F* hist = (TH1F*) fHistdEdx->At(GetHistID(AliPID::kPion, 1));
 141   //fNbins   = hist->GetNbinsX();
 142   //fBinSize = hist->GetBinWidth(1);
 143
 144   return kTRUE;
 145
 146
 147 }
 148
 149 //_________________________________________________________________________
 150 TObject* AliTRDCalPIDLQ::GetModel(Int_t ip, Int_t iType, Int_t iplane) const          // iType not needed
 151 {
 152   //
 153   // Returns one selected dEdx histogram
 154   //
 155
 156   if (iType < 0 || iType >= AliPID::kSPECIES) return 0x0;
 157   if(ip<0 || ip>= kNMom ) return 0x0;
 158
 159   AliInfo(Form("Retrive dEdx histogram for %s of %5.2f GeV/c", fPartName[iType], fgTrackMomentum[ip]));
 160
 161   return fModel->At(GetModelID(ip, iType, iplane));
 162 }
 163
 164 //_________________________________________________________________________
 165 Double_t AliTRDCalPIDLQ::GetProbability(Int_t spec, Float_t mom
 166                                       , const Float_t * const dedx
 167                                       , Float_t length, Int_t iplane) const
 168 {
 169   //
 170   // Core function of AliTRDCalPID class for calculating the
 171   // likelihood for species "spec" (see AliTRDtrack::kNspecie) of a
 172   // given momentum "mom" and a given dE/dx (slice "dedx") yield per
 173   // layer
 174   //
 175
 176   if (spec < 0 || spec >= AliPID::kSPECIES) return 0.;
 177
 178   //Double_t dedx   = dedx1/fMeanChargeRatio;
 179
 180   // find the interval in momentum and track segment length which applies for this data
 181   Int_t ilength = 1;
 182   while(ilength<kNLength-1 && length>fgTrackSegLength[ilength]){
 183     ilength++;
 184   }
 185   Int_t imom = 1;
 186   while(imom<kNMom-1 && mom>fgTrackMomentum[imom]) imom++;
 187
 188   Int_t nbinsx, nbinsy;
 189   TAxis *ax = 0x0, *ay = 0x0;
 190   Double_t lq1, lq2;
 191   Double_t mom1 = fgTrackMomentum[imom-1], mom2 = fgTrackMomentum[imom];
 192   TH2 *hist = 0x0;
 193   if(!(hist = (TH2D*)fModel->At(GetModelID(imom-1, spec, iplane)))){
 194     AliInfo(Form("Looking for spec(%d) mom(%f) Ex(%f) Ey(%f) length(%f)", spec, mom, dedx[0], dedx[1], length));
 195     AliError(Form("EHistogram id %d not found in DB.", GetModelID(imom-1, spec, iplane)));
 196     return 0.;
 197   }
 198   ax = hist->GetXaxis(); nbinsx = ax->GetNbins();
 199   ay = hist->GetYaxis(); nbinsy = ay->GetNbins();
 200   Float_t x = dedx[0]+dedx[1], y = dedx[2];
 201         Bool_t kX = (x < ax->GetBinUpEdge(nbinsx));
 202   Bool_t kY = (y < ay->GetBinUpEdge(nbinsy));
 203   if(kX)
 204     if(kY) lq1 = hist->GetBinContent( hist->FindBin(x, y));
 205     //fEstimator->Estimate2D2(hist, x, y);
 206     else lq1 = hist->GetBinContent(ax->FindBin(x), nbinsy);
 207   else
 208     if(kY) lq1 = hist->GetBinContent(nbinsx, ay->FindBin(y));
 209     else lq1 = hist->GetBinContent(nbinsx, nbinsy);
 210
 211
 212   if(!(hist = (TH2D*)fModel->At(GetModelID(imom, spec, iplane)))){
 213     AliInfo(Form("Looking for spec(%d) mom(%f) Ex(%f) Ey(%f) length(%f)", spec, mom, dedx[0], dedx[1], length));
 214     AliError(Form("EHistogram id %d not found in DB.", GetModelID(imom, spec, iplane)));
 215     return lq1;
 216   }
 217   if(kX)
 218     if(kY) lq2 = hist->GetBinContent( hist->FindBin(x, y));
 219     //fEstimator->Estimate2D2(hist, x, y);
 220     else lq2 = hist->GetBinContent(ax->FindBin(x), nbinsy);
 221   else
 222     if(kY) lq2 = hist->GetBinContent(nbinsx, ay->FindBin(y));
 223     else lq2 = hist->GetBinContent(nbinsx, nbinsy);
 224
 225         // return interpolation over momentum binning
 226         if(mom < fgTrackMomentum[0]) return lq1;
 227         else if(mom > fgTrackMomentum[kNMom-1]) return lq2;
 228         else return lq1 + (lq2 - lq1)*(mom - mom1)/(mom2 - mom1);
 229
 230 }
 231
 232 //_________________________________________________________________________
 233 void AliTRDCalPIDLQ::Init()
 234 {
 235   //
 236   // Initialization
 237   //
 238
 239   fModel = new TObjArray(AliPID::kSPECIES  * kNMom);
 240   fModel -> SetOwner();
 241
 242 }
 243
 244 //_________________________________________________________________________
 245 Int_t AliTRDCalPIDLQ::GetModelID(Int_t mom, Int_t spec, Int_t /*ii*/) const
 246 {
 247   // returns the ID of the LQ distribution (55 Histos, ID from 1 to 55)
 248
 249   return spec * AliTRDCalPID::kNMom + mom;
 250 }