PWG1/TRD/AliTRDpidRefMakerLQ.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: AliTRDpidRefMakerLQ.cxx 34163 2009-08-07 11:28:51Z cblume $ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //
  20 //
  21 //  TRD calibration class for building reference data for PID
  22 //  - 2D reference histograms (responsible A.Bercuci)
  23 //  - 3D reference histograms (not yet implemented) (responsible A.Bercuci)
  24 //  - Neural Network (responsible A.Wilk)
  25 //
  26 //   Origin
  27 //   Alex Bercuci  (A.Bercuci@gsi.de)
  28 //
  29 ///////////////////////////////////////////////////////////////////////////////
  30
  31 #include <TFile.h>
  32 #include <TROOT.h>
  33 #include <TTree.h>
  34 #include <TMath.h>
  35 #include <TEventList.h>
  36 #include <TH2D.h>
  37 #include <TH2I.h>
  38 #include <TPrincipal.h>
  39 //#include <TVector3.h>
  40 //#include <TLinearFitter.h>
  41 #include <TCanvas.h>
  42 //#include <TEllipse.h>
  43 //#include <TMarker.h>
  44
  45 #include "AliLog.h"
  46 #include "../../STAT/TKDPDF.h"
  47 #include "AliTRDpidRefMakerLQ.h"
  48 #include "../Cal/AliTRDCalPID.h"
  49 #include "../Cal/AliTRDCalPIDLQ.h"
  50 #include "AliTRDseedV1.h"
  51 #include "AliTRDcalibDB.h"
  52 #include "AliTRDgeometry.h"
  53
  54 ClassImp(AliTRDpidRefMakerLQ)
  55
  56 //__________________________________________________________________
  57 AliTRDpidRefMakerLQ::AliTRDpidRefMakerLQ()
  58   :AliTRDpidRefMaker("PIDrefMakerLQ", "PID(LQ) Reference Maker")
  59 {
  60   //
  61   // AliTRDpidRefMakerLQ default constructor
  62   //
  63 }
  64
  65 //__________________________________________________________________
  66 AliTRDpidRefMakerLQ::~AliTRDpidRefMakerLQ()
  67 {
  68   //
  69   // AliTRDCalPIDQRef destructor
  70   //
  71 }
  72
  73 //________________________________________________________________________
  74 void AliTRDpidRefMakerLQ::CreateOutputObjects()
  75 {
  76   // Create histograms
  77   // Called once
  78
  79   AliTRDpidRefMaker::CreateOutputObjects();
  80
  81   // save dE/dx references
  82   TH2 *h2 = 0x0;
  83   for(Int_t ip=AliTRDCalPID::kNMom; ip--;){
  84     TObjArray *arr = new TObjArray(AliPID::kSPECIES);
  85     arr->SetName(Form("Pbin%02d", ip));
  86     for(Int_t is=AliPID::kSPECIES; is--;) {
  87       h2 = new TH2D(Form("h%s%d", AliPID::ParticleShortName(is), ip), Form("%s ref. dEdx @ Pbin[%d]", AliPID::ParticleName(is), ip), 50, 5., 10., 50, 5., 10.);
  88       h2->GetXaxis()->SetTitle("log(dE/dx_{am}) [au]");
  89       h2->GetYaxis()->SetTitle("log(dE/dx_{dr}) [au]");
  90       h2->GetZaxis()->SetTitle("#");
  91       arr->AddAt(h2, is);
  92     }
  93     fContainer->AddAt(arr, 1+ip);
  94   }
  95 }
  96
  97 /*
  98 //________________________________________________________________________
  99 Float_t* AliTRDpidRefMakerLQ::CookdEdx(AliTRDseedV1 *trklt)
 100 {
 101 // Fill dEdx array for multidim LQ PID
 102
 103   trklt->CookdEdx(AliTRDpidUtil::kLQslices);
 104   const Float_t *dedx = trklt->GetdEdx();
 105   if(dedx[0]+dedx[1] <= 0.) return 0x0;
 106   if(dedx[2] <= 0.) return 0x0;
 107
 108   fdEdx[0] = TMath::Log(dedx[0]+dedx[1]);
 109   fdEdx[1] = TMath::Log(dedx[2]);
 110   return fdEdx;
 111 }*/
 112
 113 //__________________________________________________________________
 114 TObject* AliTRDpidRefMakerLQ::GetOCDBEntry(Option_t *opt)
 115 {
 116 // Steer loading of OCDB LQ PID
 117
 118   TDirectoryFile *d = 0x0;
 119   if(!TFile::Open(Form("TRD.Calib%s.root", GetName()))) return 0x0;
 120   if(!(d=(TDirectoryFile*)gFile->Get(Form("PDF_%s", opt)))) return 0x0;
 121   AliTRDCalPIDLQ *cal = new AliTRDCalPIDLQ("pidLQ", "LQ TRD PID object");
 122   cal->LoadPDF(d);
 123   return cal;
 124 }
 125
 126 //__________________________________________________________________
 127 Bool_t AliTRDpidRefMakerLQ::GetRefFigure(Int_t ifig)
 128 {
 129 // Steering reference picture
 130
 131   if(ifig<0 || ifig>AliTRDCalPID::kNMom-1){
 132     AliError("Ref fig requested outside definition.");
 133     return kFALSE;
 134   }
 135   if(!gPad){
 136     AliWarning("Please provide a canvas to draw results.");
 137     return kFALSE;
 138   }
 139
 140   TObjArray *arr = (TObjArray*)fContainer->At(ifig);
 141   gPad->Divide(3, 2, 1.e-5, 1.e-5);
 142   TList *l=gPad->GetListOfPrimitives();
 143   for(Int_t is=0; is<AliPID::kSPECIES; is++){
 144     ((TVirtualPad*)l->At(is))->cd();
 145     ((TH2*)arr->At(is))->Draw("cont4z");
 146   }
 147
 148   return kTRUE;
 149 }
 150
 151 //________________________________________________________________________
 152 Bool_t AliTRDpidRefMakerLQ::PostProcess()
 153 {
 154 // Analyse merged dedx = f(p) distributions.
 155 //   - select momentum - species bins
 156 //   - rotate to principal components
 157 //   - locally interpolate with TKDPDF
 158 //   - save interpolation to monitoring histograms
 159 //   - write pdf to file for loading to OCDB
 160 //
 161
 162
 163   TFile *fCalib = TFile::Open(Form("TRD.Calib%s.root", GetName()), "update");
 164   fData = dynamic_cast<TTree*>(gFile->Get(GetName()));
 165   if (!fData) {
 166     AliError("Calibration data not available");
 167     return kFALSE;
 168   }
 169   TObjArray *o = 0x0;
 170   if(!(o = (TObjArray*)gFile->Get(Form("Moni%s", GetName())))){
 171     AliWarning("Missing monitoring container.");
 172     return kFALSE;
 173   }
 174   fContainer = (TObjArray*)o->Clone("monitor");
 175
 176   TDatime d;
 177   TDirectoryFile *pdfs = new TDirectoryFile(Form("PDF_%d", d.GetDate()), "PDFs for LQ TRD-PID", "", gFile);
 178   pdfs->Write();
 179   AliDebug(2, Form("Data[%d]", fData->GetEntries()));
 180   pdfs->cd();
 181
 182   //TCanvas *cc = new TCanvas("cc", "", 500, 500);
 183   LinkPIDdata();
 184   Float_t *data[] = {0x0, 0x0};
 185   // allocate storage
 186   data[0] = new Float_t[kMaxStat];data[1] = new Float_t[kMaxStat];
 187   for(Int_t ip=AliTRDCalPID::kNMom; ip--; ){
 188     for(Int_t is=AliPID::kSPECIES; is--;) {
 189       Int_t n(0); // index of data
 190       for(Int_t itrk=0; (itrk < fData->GetEntries()) && (n<kMaxStat); itrk++){
 191         if(!(fData->GetEntry(itrk))) continue;
 192         if(fPIDbin!=is) continue;
 193         for(Int_t ily=fPIDdataArray->fNtracklets; ily--;){
 194           if((fPIDdataArray->fData[ily].fPLbin & 0xf)!= ip) continue;
 195
 196           Float_t dedx[] = {0., 0.};
 197           for(Int_t islice=AliTRDCalPID::kNSlicesNN; islice--;){
 198             Int_t jslice = islice>kNN2LQtransition;
 199             dedx[jslice]+=fPIDdataArray->fData[ily].fdEdx[islice];
 200           }
 201
 202           // check data integrity
 203           if(dedx[0]<1.e-30) continue;
 204           if(dedx[1]<1.e-30) continue;
 205
 206           // store data
 207           data[0][n] = TMath::Log(dedx[0]);
 208           data[1][n] = TMath::Log(dedx[1]);
 209           n++; if(n==kMaxStat) break;
 210         }
 211       }
 212
 213       // estimate bucket statistics
 214       Int_t nb(kMinBuckets), // number of buckets
 215             ns(Int_t(Float_t(n)/nb));    //statistics/bucket
 216
 217 // if(Float_t(n)/nb < 220.) ns = 200; // 7% stat error
 218 //       else if(Float_t(n)/nb < 420.) ns = 400; // 5% stat error
 219
 220       AliDebug(2, Form("pBin[%d] sBin[%d] n[%d] ns[%d] nb[%d]", ip, is, n, ns, nb));
 221       if(ns<Int_t(kMinStat)){
 222         AliWarning(Form("Not enough entries [%d] for %s[%d].", n, AliPID::ParticleShortName(is), ip));
 223         continue;
 224       }
 225
 226       // build PDF
 227       TKDPDF pdf(n, 2, ns, data);
 228       pdf.SetCOG(kFALSE);
 229       pdf.SetWeights();
 230       pdf.SetStore();
 231       pdf.SetAlpha(5.);
 232       pdf.GetStatus();
 233       Float_t *c, v, ve; Double_t r, e, rxy[2];
 234       for(Int_t in=pdf.GetNTNodes(); in--;){
 235         pdf.GetCOGPoint(in, c, v, ve);
 236         rxy[0] = (Double_t)c[0];rxy[1] = (Double_t)c[1];
 237         pdf.Eval(rxy, r, e, kTRUE);
 238       }
 239 //       // visual on-line monitoring
 240 //       pdf.DrawProjection();cc->Modified(); cc->Update(); cc->SaveAs(Form("pdf_%s%02d.gif", AliPID::ParticleShortName(is), ip));
 241 //       cc->SaveAs(Form("%s_%s%02d.gif", GetName(), AliPID::ParticleShortName(is), ip));
 242
 243       // save a discretization of the PDF for result monitoring
 244       TH2 *h2s = (TH2D*)((TObjArray*)fContainer->At(ip))->At(is);
 245       TAxis *ax = h2s->GetXaxis(), *ay = h2s->GetYaxis();
 246       h2s->Clear();
 247       for(int ix=1; ix<=ax->GetNbins(); ix++){
 248         rxy[0] = ax->GetBinCenter(ix);
 249         for(int iy=1; iy<=ay->GetNbins(); iy++){
 250           rxy[1] = ay->GetBinCenter(iy);
 251
 252           Double_t rr,ee;
 253           pdf.Eval(rxy, rr, ee, kFALSE);
 254           if(rr<0. || ee/rr>.15) continue; // 15% relative error
 255           //printf("x[%2d] x[%2d] r[%f] e[%f]\n", ix, iy, rr, ee);
 256           h2s->SetBinContent(ix, iy, rr);
 257         }
 258       }
 259
 260       // write results to output array
 261       //pdf.GetStatus();
 262       pdf.Write(Form("%s[%d]", AliPID::ParticleShortName(is), ip));
 263     }
 264   }
 265   delete [] data[0]; delete [] data[1];
 266   pdfs->Write();
 267   fCalib->Close(); delete fCalib;
 268
 269   return kTRUE; // testing protection
 270 }
 271
 272