copy TRD performance train to PWG1
[u/mrichter/AliRoot.git] / PWG1 / TRD / AliTRDpidRefMakerLQ.cxx
CommitLineData
1ee39b3a 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
54ClassImp(AliTRDpidRefMakerLQ)
55
56//__________________________________________________________________
57AliTRDpidRefMakerLQ::AliTRDpidRefMakerLQ()
58 :AliTRDpidRefMaker("PIDrefMakerLQ", "PID(LQ) Reference Maker")
59{
60 //
61 // AliTRDpidRefMakerLQ default constructor
62 //
63}
64
65//__________________________________________________________________
66AliTRDpidRefMakerLQ::~AliTRDpidRefMakerLQ()
67{
68 //
69 // AliTRDCalPIDQRef destructor
70 //
71}
72
73//________________________________________________________________________
74void 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//________________________________________________________________________
99Float_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//__________________________________________________________________
114TObject* 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//__________________________________________________________________
127Bool_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//________________________________________________________________________
152Bool_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