TRD/qaRec/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 "AliTRDseedV1.h"
  50 #include "AliTRDcalibDB.h"
  51 #include "AliTRDgeometry.h"
  52
  53 ClassImp(AliTRDpidRefMakerLQ)
  54
  55 //__________________________________________________________________
  56 AliTRDpidRefMakerLQ::AliTRDpidRefMakerLQ()
  57   :AliTRDpidRefMaker("pidRefMakerLQ", "PID(LQ) Reference Maker")
  58   ,fPbin(-1)
  59   ,fSbin(-1)
  60   ,fResults(0x0)
  61 {
  62   //
  63   // AliTRDpidRefMakerLQ default constructor
  64   //
  65 }
  66
  67 //__________________________________________________________________
  68 AliTRDpidRefMakerLQ::~AliTRDpidRefMakerLQ()
  69 {
  70   //
  71   // AliTRDCalPIDQRef destructor
  72   //
  73   if(fResults) {
  74     fResults->Delete();
  75     delete fResults;
  76   }
  77 }
  78
  79 //________________________________________________________________________
  80 void AliTRDpidRefMakerLQ::CreateOutputObjects()
  81 {
  82   // Create histograms
  83   // Called once
  84
  85   AliTRDpidRefMaker::CreateOutputObjects();
  86
  87   // save dE/dx references
  88   TH2 *h2 = 0x0;
  89   for(Int_t ip=AliTRDCalPID::kNMom; ip--;){
  90     TObjArray *arr = new TObjArray(AliPID::kSPECIES);
  91     arr->SetName(Form("Pbin%02d", ip));
  92     for(Int_t is=AliPID::kSPECIES; is--;) {
  93       h2 = new TH2I(Form("h%s%d", AliPID::ParticleShortName(is), ip), Form("%s ref. dEdx @ Pbin[%d]", AliPID::ParticleName(is), ip), 50, 5., 10., 50, 5., 10.);
  94       h2->GetXaxis()->SetTitle("log(dE/dx_{am}) [au]");
  95       h2->GetYaxis()->SetTitle("log(dE/dx_{dr}) [au]");
  96       h2->GetZaxis()->SetTitle("#");
  97       arr->AddAt(h2, is);
  98     }
  99     fContainer->AddAt(arr, 1+ip);
 100   }
 101
 102   fData = new TTree(GetName(), Form("Reference data for %s", GetName()));
 103   fData->Branch("s", &fSbin, "l/b");
 104   fData->Branch("p", &fPbin, "p/b");
 105   fData->Branch("dEdx" , fdEdx , Form("dEdx[%d]/F", GetNslices()));
 106 }
 107
 108
 109 //________________________________________________________________________
 110 Float_t* AliTRDpidRefMakerLQ::CookdEdx(AliTRDseedV1 *trklt)
 111 {
 112   trklt->CookdEdx(AliTRDpidUtil::kLQslices);
 113   const Float_t *dedx = trklt->GetdEdx();
 114   if(dedx[0]+dedx[1] <= 0.) return 0x0;
 115   if(dedx[2] <= 0.) return 0x0;
 116
 117   fdEdx[0] = TMath::Log(dedx[0]+dedx[1]);
 118   fdEdx[1] = TMath::Log(dedx[2]);
 119   return fdEdx;
 120 }
 121
 122 #include "../Cal/AliTRDCalPIDLQ.h"
 123
 124 //__________________________________________________________________
 125 TObject* AliTRDpidRefMakerLQ::GetOCDBEntry(Option_t *opt)
 126 {
 127   TDirectoryFile *d = 0x0;
 128   if(!TFile::Open(Form("TRD.Calib%s.root", GetName()))) return 0x0;
 129   if(!(d=(TDirectoryFile*)gFile->Get(Form("PDF_%s", opt)))) return 0x0;
 130   AliTRDCalPIDLQ *cal = new AliTRDCalPIDLQ("pidLQ", "LQ TRD PID object");
 131   cal->LoadPDF(d);
 132   return cal;
 133 }
 134
 135 //__________________________________________________________________
 136 Bool_t AliTRDpidRefMakerLQ::GetRefFigure(Int_t ifig)
 137 {
 138   if(ifig<0 || ifig>AliTRDCalPID::kNMom-1){
 139     AliError("Ref fig requested outside definition.");
 140     return kFALSE;
 141   }
 142   if(!fResults){
 143     AliError("No results processed.");
 144     return kFALSE;
 145   }
 146   if(!gPad){
 147     AliWarning("Please provide a canvas to draw results.");
 148     return kFALSE;
 149   }
 150
 151   TObjArray *arr = (TObjArray*)fResults->At(ifig);
 152   gPad->Divide(3, 2, 1.e-5, 1.e-5);
 153   TList *l=gPad->GetListOfPrimitives();
 154   for(Int_t is=0; is<AliPID::kSPECIES; is++){
 155     ((TVirtualPad*)l->At(is))->cd();
 156     ((TH2*)arr->At(is))->Draw("cont4z");
 157   }
 158
 159   return kTRUE;
 160 }
 161
 162
 163 //________________________________________________________________________
 164 void AliTRDpidRefMakerLQ::Fill()
 165 {
 166   // get particle type
 167   fSbin = TMath::LocMax(AliPID::kSPECIES, fPID);
 168   // get momentum bin
 169   fPbin = AliTRDpidUtil::GetMomentumBin(fP);
 170   // fill data
 171   fData->Fill();
 172   // fill monitor
 173   ((TH2*)fContainer->At(0))->Fill(fSbin, fPbin);
 174   TH2* h2 = (TH2*)((TObjArray*)fContainer->At(1+fPbin))->At(fSbin);
 175   h2->Fill(fdEdx[0], fdEdx[1]);
 176 }
 177
 178 //________________________________________________________________________
 179 Bool_t AliTRDpidRefMakerLQ::PostProcess()
 180 {
 181 // Analyse merged dedx = f(p) distributions.
 182 //   - select momentum - species bins
 183 //   - rotate to principal components
 184 //   - locally interpolate with TKDPDF
 185 //   - save interpolation to monitoring histograms
 186 //   - write pdf to file for loading to OCDB
 187 //
 188
 189
 190   TFile *fCalib = TFile::Open(Form("TRD.Calib%s.root", GetName()), "update");
 191   fData = dynamic_cast<TTree*>(gFile->Get(GetName()));
 192   if (!fData) {
 193     AliError("Tree not available");
 194     return kFALSE;
 195   }
 196   TDatime d;
 197   TDirectoryFile *pdfs = new TDirectoryFile(Form("PDF_%d", d.GetDate()), "PDFs for LQ TRD-PID", "", gFile);
 198   pdfs->Write();
 199   AliDebug(2, Form("Data[%d]", fData->GetEntries()));
 200
 201   // save a volatile PDF representation
 202   gROOT->cd();
 203   TH2 *h2 = 0x0;
 204   fResults = new TObjArray(AliTRDCalPID::kNMom);
 205   for(Int_t ip=AliTRDCalPID::kNMom; ip--;){
 206     TObjArray *arr = new TObjArray(AliPID::kSPECIES);
 207     arr->SetName(Form("Pbin%02d", ip));
 208     for(Int_t is=AliPID::kSPECIES; is--;) {
 209       h2 = new TH2D(Form("i%s%d", AliPID::ParticleShortName(is), ip), Form("Interpolated %s dEdx @ Pbin[%d]", AliPID::ParticleName(is), ip), 50, -3., 3., 50, -3., 3.);
 210       h2->GetXaxis()->SetTitle("log(dE/dx^{*}_{am}) [au]");
 211       h2->GetYaxis()->SetTitle("log(dE/dx^{*}_{dr}) [au]");
 212       h2->GetZaxis()->SetTitle("#");
 213       arr->AddAt(h2, is);
 214     }
 215     fResults->AddAt(arr, ip);
 216   }
 217   pdfs->cd();
 218
 219   TCanvas *cc = new TCanvas("cc", "", 500, 500);
 220
 221   Float_t *data[] = {0x0, 0x0};
 222   TPrincipal principal(2, "ND");
 223   for(Int_t ip=AliTRDCalPID::kNMom; ip--; ){
 224     for(Int_t is=AliPID::kSPECIES; is--;) {
 225       principal.Clear();
 226       Int_t n = fData->Draw("dEdx[0]:dEdx[1]", Form("p==%d&&s==%d", ip, is), "goff");
 227
 228       // estimate bucket statistics
 229       Int_t nb(kMinBuckets), // number of buckets
 230             ns(Int_t(Float_t(n)/nb));    //statistics/bucket
 231
 232 // if(Float_t(n)/nb < 220.) ns = 200; // 7% stat error
 233 //       else if(Float_t(n)/nb < 420.) ns = 400; // 5% stat error
 234
 235       AliDebug(2, Form("pBin[%d] sBin[%d] N[%d] n[%d] nb[%d]", ip, is, n, ns, nb));
 236       if(ns<Int_t(kMinStat)){
 237         AliWarning(Form("Not enough entries [%d] for %s[%d].", ns, AliPID::ParticleShortName(is), ip));
 238         continue;
 239       }
 240       // allocate storage
 241       data[0] = new Float_t[n];data[1] = new Float_t[n];
 242
 243       // fill and make principal
 244       Double_t *v1 = fData->GetV1(), *v2 = fData->GetV2();
 245       while(n--){
 246         Double_t dedx[] = {v1[n], v2[n]};
 247         principal.AddRow(dedx);
 248       }
 249       principal.MakePrincipals();
 250
 251 //       // calculate covariance ellipse
 252 //       eValues  = principal.GetEigenValues();
 253 //       x0  = 0.;
 254 //       y0  = 0.;
 255 //       rx  = 3.5*sqrt((*eValues)[0]);
 256 //       ry  = 3.5*sqrt((*eValues)[1]);
 257
 258       // rotate to principal axis
 259       const Double_t *xx; Double_t rxy[2];
 260       while((xx = principal.GetRow(++n))){
 261         principal.X2P(xx, rxy);
 262         data[0][n]=rxy[0]; data[1][n]=rxy[1];
 263       }
 264
 265       // build PDF
 266       TKDPDF pdf(n, 2, ns, data);
 267       pdf.SetStore();
 268       pdf.SetAlpha(5.);
 269       //pdf.GetStatus();
 270       Float_t *c, v, ve; Double_t r, e;
 271       for(Int_t in=pdf.GetNTNodes(); in--;){
 272         pdf.GetCOGPoint(in, c, v, ve);
 273         rxy[0] = (Double_t)c[0];rxy[1] = (Double_t)c[1];
 274         pdf.Eval(rxy, r, e, kTRUE);
 275       }
 276       // visual on-line monitoring
 277       pdf.DrawBins(0,1,-5.,5.,-8.,8.);cc->Modified(); cc->Update(); cc->SaveAs(Form("pdf_%s%02d.gif", AliPID::ParticleShortName(is), ip));
 278       cc->SaveAs(Form("%s_%s%02d.gif", GetName(), AliPID::ParticleShortName(is), ip));
 279
 280       // save a discretization of the PDF for result monitoring
 281       TH2 *h2s = (TH2D*)((TObjArray*)fResults->At(ip))->At(is);
 282       TAxis *ax = h2s->GetXaxis(), *ay = h2s->GetYaxis();
 283       for(int ix=1; ix<=ax->GetNbins(); ix++){
 284         rxy[0] = ax->GetBinCenter(ix);
 285         for(int iy=1; iy<=ay->GetNbins(); iy++){
 286           rxy[1] = ay->GetBinCenter(iy);
 287
 288           Double_t r,e;
 289           pdf.Eval(rxy, r, e);
 290           if(r<0. || e/r>.15) continue; // 15% relative error
 291           //printf("x[%2d] x[%2d] r[%f] e[%f]\n", ix, iy, r, e);
 292           h2s->SetBinContent(ix, iy, r);
 293         }
 294       }
 295
 296       // write results to output array
 297       pdf.Write(Form("%s[%d]", AliPID::ParticleShortName(is), ip));
 298
 299       delete [] data[0]; delete [] data[1];
 300     }
 301   }
 302   pdfs->Write();
 303   fCalib->Close(); delete fCalib;
 304
 305   return kTRUE; // testing protection
 306 }
 307
 308