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