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 <TH1F.h>
  30 #include <TH2F.h>
  31 #include <TFile.h>
  32 #include <TROOT.h>
  33
  34 #include "AliLog.h"
  35 #include "AliPID.h"
  36
  37 #include "AliTRDCalPIDLQ.h"
  38 #include "AliTRDcalibDB.h"
  39
  40 ClassImp(AliTRDCalPIDLQ)
  41
  42 Float_t AliTRDCalPIDLQ::fTrackSegLength[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::LoadReferences(Char_t *refFile)
  79 {
  80   //
  81   // Read the TRD dEdx histograms.
  82   //
  83
  84         Int_t nTimeBins = 22;
  85         // Get number of time bins from CDB
  86         AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
  87         if(!calibration){
  88                 AliWarning(Form("No AliTRDcalibDB available. Using %d time bins.", nTimeBins));
  89         }else{
  90                 if(calibration->GetRun() > -1) nTimeBins = calibration->GetNumberOfTimeBins();
  91                 else AliWarning(Form("Run number not set. Using %d time bins.", nTimeBins));
  92         }
  93
  94
  95   // Read histogram Root file
  96   TFile *histFile = TFile::Open(refFile, "READ");
  97   if (!histFile || !histFile->IsOpen()) {
  98     AliError(Form("Opening TRD histgram file %s failed", refFile));
  99     return kFALSE;
 100   }
 101   gROOT->cd();
 102
 103   // Read histograms
 104   for (Int_t iparticle = 0; iparticle < AliPID::kSPECIES; iparticle++){
 105     for (Int_t imom = 0; imom < kNMom; imom++){
 106                         TH2D* hist = (TH2D*)histFile->Get(Form("h2dEdx%s%d", fPartSymb[iparticle], imom/*, ilength*/))->Clone();
 107                         hist->Scale(1./hist->Integral());
 108                         fModel->AddAt(hist, GetModelID(imom, iparticle, 0));
 109
 110 //                      if (iparticle != AliPID::kElectron && iparticle != AliPID::kPion) continue;
 111 //
 112 //                      TH1F* ht = (TH1F*)histFile->Get(Form("h1MaxTB%s%02d", fPartSymb[iparticle], imom))->Clone();
 113 //                      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));
 114 //                      ht->Scale(1./ht->Integral());
 115 //                      fHistTimeBin->AddAt(ht, ((iparticle==AliPID::kElectron)?0:1)*kNMom + imom);
 116                 }
 117   }
 118
 119   histFile->Close();
 120   delete histFile;
 121
 122   // Number of bins and bin size
 123   //TH1F* hist = (TH1F*) fHistdEdx->At(GetHistID(AliPID::kPion, 1));
 124   //fNbins   = hist->GetNbinsX();
 125   //fBinSize = hist->GetBinWidth(1);
 126
 127   return kTRUE;
 128
 129
 130 }
 131
 132 //_________________________________________________________________________
 133 TObject* AliTRDCalPIDLQ::GetModel(Int_t ip, Int_t iType, Int_t iplane) const          // iType not needed
 134 {
 135   //
 136   // Returns one selected dEdx histogram
 137   //
 138
 139   if (iType < 0 || iType >= AliPID::kSPECIES) return 0x0;
 140         if(ip<0 || ip>= kNMom ) return 0x0;
 141
 142         AliInfo(Form("Retrive dEdx histogram for %s of %5.2f GeV/c", fPartName[iType], fTrackMomentum[ip]));
 143
 144   return fModel->At(GetModelID(ip, iType, iplane));
 145 }
 146
 147 //_________________________________________________________________________
 148 Double_t AliTRDCalPIDLQ::GetProbability(Int_t spec, Float_t mom, Float_t *dedx, Float_t length, Int_t iplane) const
 149 {
 150   //
 151         // Core function of AliTRDCalPID class for calculating the
 152         // likelihood for species "spec" (see AliTRDtrack::kNspecie) of a
 153         // given momentum "mom" and a given dE/dx (slice "dedx") yield per
 154         // layer
 155   //
 156
 157         if (spec < 0 || spec >= AliPID::kSPECIES) return 0.;
 158
 159         //Double_t dedx   = dedx1/fMeanChargeRatio;
 160
 161         // find the interval in momentum and track segment length which applies for this data
 162         Int_t ilength = 1;
 163   while(ilength<kNLength-1 && length>fTrackSegLength[ilength]){
 164                 ilength++;
 165         }
 166         Int_t imom = 1;
 167   while(imom<kNMom-1 && mom>fTrackMomentum[imom]) imom++;
 168
 169         Int_t nbinsx, nbinsy;
 170         TAxis *ax = 0x0, *ay = 0x0;
 171         Double_t LQ1, LQ2;
 172         Double_t mom1 = fTrackMomentum[imom-1], mom2 = fTrackMomentum[imom];
 173         TH2 *hist = 0x0;
 174         if(!(hist = (TH2D*)fModel->At(GetModelID(imom-1, spec, iplane)))){
 175                 AliInfo(Form("Looking for spec(%d) mom(%f) Ex(%f) Ey(%f) length(%f)", spec, mom, dedx[0], dedx[1], length));
 176                 AliError(Form("EHistogram id %d not found in DB.", GetModelID(imom-1, spec, iplane)));
 177                 return 0.;
 178         }
 179         ax = hist->GetXaxis(); nbinsx = ax->GetNbins();
 180         ay = hist->GetYaxis(); nbinsy = ay->GetNbins();
 181         Float_t x = dedx[0]+dedx[1], y = dedx[2];
 182         Bool_t kX = (x < ax->GetBinUpEdge(nbinsx));
 183         Bool_t kY = (y < ay->GetBinUpEdge(nbinsy));
 184         if(kX)
 185                 if(kY) LQ1 = hist->GetBinContent( hist->FindBin(x, y));
 186     //fEstimator->Estimate2D2(hist, x, y);
 187                 else LQ1 = hist->GetBinContent(ax->FindBin(x), nbinsy);
 188         else
 189                 if(kY) LQ1 = hist->GetBinContent(nbinsx, ay->FindBin(y));
 190                 else LQ1 = hist->GetBinContent(nbinsx, nbinsy);
 191
 192
 193         if(!(hist = (TH2D*)fModel->At(GetModelID(imom, 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, spec, iplane)));
 196                 return LQ1;
 197         }
 198         if(kX)
 199                 if(kY) LQ2 = hist->GetBinContent( hist->FindBin(x, y));
 200     //fEstimator->Estimate2D2(hist, x, y);
 201                 else LQ2 = hist->GetBinContent(ax->FindBin(x), nbinsy);
 202         else
 203                 if(kY) LQ2 = hist->GetBinContent(nbinsx, ay->FindBin(y));
 204                 else LQ2 = hist->GetBinContent(nbinsx, nbinsy);
 205
 206         // return interpolation over momentum binning
 207         if(mom < fTrackMomentum[0]) return LQ1;
 208         else if(mom > fTrackMomentum[kNMom-1]) return LQ2;
 209         else return LQ1 + (LQ2 - LQ1)*(mom - mom1)/(mom2 - mom1);
 210
 211 }
 212
 213 //_________________________________________________________________________
 214 void AliTRDCalPIDLQ::Init()
 215 {
 216   //
 217   // Initialization
 218   //
 219
 220   fModel = new TObjArray(AliPID::kSPECIES  * kNMom);
 221   fModel -> SetOwner();
 222
 223 }
 224
 225 //_________________________________________________________________________
 226 Int_t AliTRDCalPIDLQ::GetModelID(Int_t mom, Int_t spec, Int_t) const
 227 {
 228   // returns the ID of the LQ distribution (55 Histos, ID from 1 to 55)
 229
 230         return spec * AliTRDCalPID::kNMom + mom;
 231 }