prepare PID LQ ref maker for production
[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)
1ee39b3a 24//
25// Origin
26// Alex Bercuci (A.Bercuci@gsi.de)
27//
28///////////////////////////////////////////////////////////////////////////////
29
b91fdd71 30#include <TSystem.h>
1ee39b3a 31#include <TROOT.h>
b91fdd71 32#include <TFile.h>
1ee39b3a 33#include <TTree.h>
34#include <TMath.h>
35#include <TEventList.h>
36#include <TH2D.h>
37#include <TH2I.h>
1ee39b3a 38#include <TCanvas.h>
1ee39b3a 39
40#include "AliLog.h"
da27d983 41#include "../STAT/TKDPDF.h"
1ee39b3a 42#include "AliTRDpidRefMakerLQ.h"
da27d983 43#include "Cal/AliTRDCalPID.h"
44#include "Cal/AliTRDCalPIDLQ.h"
1ee39b3a 45#include "AliTRDseedV1.h"
46#include "AliTRDcalibDB.h"
47#include "AliTRDgeometry.h"
48
49ClassImp(AliTRDpidRefMakerLQ)
50
51//__________________________________________________________________
b91fdd71 52AliTRDpidRefMakerLQ::AliTRDpidRefMakerLQ() :
53 AliTRDpidRefMaker("PIDrefMakerLQ", "PID(LQ) Reference Maker")
1ee39b3a 54{
55 //
56 // AliTRDpidRefMakerLQ default constructor
57 //
58}
59
60//__________________________________________________________________
61AliTRDpidRefMakerLQ::~AliTRDpidRefMakerLQ()
62{
63 //
64 // AliTRDCalPIDQRef destructor
65 //
66}
67
b91fdd71 68// //________________________________________________________________________
1ee39b3a 69void AliTRDpidRefMakerLQ::CreateOutputObjects()
70{
71 // Create histograms
72 // Called once
73
b91fdd71 74 //OpenFile(0, "RECREATE");
75 fContainer = Histos();
76}
77
78
79//__________________________________________________________________
80TObjArray* 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()));
1ee39b3a 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);
b91fdd71 94 arr->SetName(Form("Pbin%02d", ip)); arr->SetOwner();
1ee39b3a 95 for(Int_t is=AliPID::kSPECIES; is--;) {
b91fdd71 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.);
1ee39b3a 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 }
b91fdd71 102 fContainer->AddAt(arr, ip);
1ee39b3a 103 }
b91fdd71 104 fNRefFigures=AliTRDCalPID::kNMom;
105 return fContainer;
1ee39b3a 106}
107
1ee39b3a 108
109//__________________________________________________________________
110TObject* 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//__________________________________________________________________
123Bool_t AliTRDpidRefMakerLQ::GetRefFigure(Int_t ifig)
124{
125// Steering reference picture
126
b91fdd71 127 if(ifig<0 || ifig>=fNRefFigures){
1ee39b3a 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
b91fdd71 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//________________________________________________________________________
165Bool_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;
1ee39b3a 179 }
b91fdd71 180 LinkPIDdata();
1ee39b3a 181
b91fdd71 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");
1ee39b3a 188 return kTRUE;
189}
190
191//________________________________________________________________________
b91fdd71 192void 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//________________________________________________________________________
1ee39b3a 201Bool_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
1ee39b3a 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);
b91fdd71 218 // allocate working storage
219 Float_t *data[] = {
220 new Float_t[kMaxStat],
221 new Float_t[kMaxStat]};
1ee39b3a 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();
b91fdd71 302 //fCalib->Close(); delete fCalib;
1ee39b3a 303
304 return kTRUE; // testing protection
305}
306
307