]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWG1/TRD/AliTRDpidRefMakerLQ.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
e THnSparse structure to store MC residuals
[u/mrichter/AliRoot.git] / PWG1 / TRD / AliTRDpidRefMakerLQ.cxx
... / ...
CommitLineData
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 "../STAT/TKDInterpolator.h"
43#include "AliTRDpidRefMakerLQ.h"
44#include "Cal/AliTRDCalPID.h"
45#include "Cal/AliTRDCalPIDLQ.h"
46#include "AliTRDseedV1.h"
47#include "AliTRDcalibDB.h"
48#include "AliTRDgeometry.h"
49#include "info/AliTRDpidInfo.h"
50
51ClassImp(AliTRDpidRefMakerLQ)
52
53//__________________________________________________________________
54AliTRDpidRefMakerLQ::AliTRDpidRefMakerLQ()
55 : AliTRDpidRefMaker()
56 ,fPDF(NULL)
57{
58 //
59 // AliTRDpidRefMakerLQ default constructor
60 //
61 SetNameTitle("TRDrefMakerLQ", "PID(LQ) Reference Maker");
62}
63
64//__________________________________________________________________
65AliTRDpidRefMakerLQ::AliTRDpidRefMakerLQ(const char *name)
66 : AliTRDpidRefMaker(name, "PID(LQ) Reference Maker")
67 ,fPDF(NULL)
68{
69 //
70 // AliTRDpidRefMakerLQ default constructor
71 //
72 DefineOutput(3, TObjArray::Class());
73}
74
75//__________________________________________________________________
76AliTRDpidRefMakerLQ::~AliTRDpidRefMakerLQ()
77{
78 //
79 // AliTRDCalPIDQRef destructor
80 //
81 if(fPDF){
82 //fPDF->Write("PDF_LQ", TObject::kSingleKey);
83 fPDF->Delete();
84 delete fPDF;
85 }
86}
87
88// //________________________________________________________________________
89void AliTRDpidRefMakerLQ::UserCreateOutputObjects()
90{
91 // Create histograms
92 // Called once
93
94 fContainer = Histos();
95 PostData(1, fContainer);
96
97 //OpenFile(2, "RECREATE");
98 fPDF = new TObjArray(AliTRDCalPIDLQ::GetNrefs());
99 fPDF->SetOwner();fPDF->SetName("PDF_LQ");
100 PostData(3, fPDF);
101}
102
103
104//__________________________________________________________________
105TObjArray* AliTRDpidRefMakerLQ::Histos()
106{
107 // Create histograms
108
109 if(fContainer) return fContainer;
110 fContainer = new TObjArray(AliTRDCalPID::kNMom);
111
112 // save dE/dx references
113 TH1 *h(NULL);
114 for(Int_t ip=AliTRDCalPID::kNMom; ip--;){
115 TObjArray *arr = new TObjArray(2*AliPID::kSPECIES);
116 arr->SetName(Form("Pbin%02d", ip)); arr->SetOwner();
117 for(Int_t is=AliPID::kSPECIES; is--;) {
118 h = new TH1D(Form("h1%s%d", AliPID::ParticleShortName(is), ip), Form("1D %s @ Pbin[%d]", AliPID::ParticleName(is), ip), 50, 7., 12.);
119 h->GetXaxis()->SetTitle("log(dE/dx) [au]");
120 h->GetYaxis()->SetTitle("#");
121 h->SetLineColor(AliTRDCalPIDLQ::GetPartColor(is));
122 arr->AddAt(h, is);
123 }
124 for(Int_t is=AliPID::kSPECIES; is--;) {
125 h = new TH2D(Form("h2%s%d", AliPID::ParticleShortName(is), ip), Form("2D %s @ Pbin[%d]", AliPID::ParticleName(is), ip), 50, 7., 12., 50, 6.5, 11.);
126 h->GetXaxis()->SetTitle("log(dE/dx_{am}) [au]");
127 h->GetYaxis()->SetTitle("log(dE/dx_{dr}) [au]");
128 h->GetZaxis()->SetTitle("#");
129 arr->AddAt(h, AliPID::kSPECIES+is);
130 }
131 fContainer->AddAt(arr, ip);
132 }
133 fNRefFigures=AliTRDCalPID::kNMom;
134 return fContainer;
135}
136
137
138//__________________________________________________________________
139TObject* AliTRDpidRefMakerLQ::GetOCDBEntry(Option_t */*opt*/)
140{
141// Steer loading of OCDB LQ PID
142
143 if(gSystem->AccessPathName(Form("TRD.Calib%s.root", GetName()), kReadPermission)){
144 AliError(Form("File TRD.Calib%s.root not readable", GetName()));
145 return NULL;
146 }
147 AliTRDCalPIDLQ *cal = new AliTRDCalPIDLQ("pidLQ", "LQ TRD PID object");
148 cal->LoadReferences(Form("TRD.Calib%s.root", GetName()));
149 return cal;
150}
151
152//__________________________________________________________________
153Bool_t AliTRDpidRefMakerLQ::GetRefFigure(Int_t ifig)
154{
155// Steering reference picture
156
157 if(ifig<0 || ifig>=fNRefFigures){
158 AliError("Ref fig requested outside definition.");
159 return kFALSE;
160 }
161 if(!gPad){
162 AliWarning("Please provide a canvas to draw results.");
163 return kFALSE;
164 }
165
166 TObjArray *arr(NULL);TList *l(NULL);TH1 *h(NULL);
167 if(!(arr = (TObjArray*)fContainer->At(ifig))){
168 AliError(Form("PDF container @ pBin[%d] missing.", ifig));
169 return kFALSE;
170 }
171 gPad->Divide(5, 2, 1.e-5, 1.e-5);l=gPad->GetListOfPrimitives();
172 for(Int_t is=0; is<AliPID::kSPECIES; is++){
173 ((TVirtualPad*)l->At(is))->cd();
174 if(!(h=(TH1*)arr->At(is))){
175 AliError(Form("1D for %s @ pBin[%d] missing.", AliPID::ParticleShortName(is), ifig));
176 return kFALSE;
177 }
178 h->GetYaxis()->SetRangeUser(0., 1.2);
179 h->Draw("l");
180
181 ((TVirtualPad*)l->At(AliPID::kSPECIES+is))->cd();
182 if(!(h=(TH1*)arr->At(AliPID::kSPECIES+is))){
183 AliError(Form("2D for %s @ pBin[%d] missing.", AliPID::ParticleShortName(is), ifig));
184 return kFALSE;
185 }
186 h->Draw("cont4z");
187 }
188
189 return kTRUE;
190}
191
192
193//________________________________________________________________________
194Bool_t AliTRDpidRefMakerLQ::Load(const Char_t *file, const Char_t *dir)
195{
196// Load tree with data in case of detached PostProcess processing.
197
198 if(gSystem->AccessPathName(file, kReadPermission)){
199 AliError(Form("File %s not readable", file));
200 return kFALSE;
201 }
202 if(!TFile::Open(file)) {
203 AliError(Form("File %s corrupted", file));
204 return kFALSE;
205 }
206 if(!gFile->cd(dir)){
207 AliWarning(Form("Couldn't cd to %s in %s.", dir, file));
208 return kFALSE;
209 }
210 if (!(fData = dynamic_cast<TTree*>(gDirectory->Get("PIDrefMaker")))) {
211 AliError("PIDref Tree not available");
212 return kFALSE;
213 }
214 LinkPIDdata();
215
216 TObjArray *o(NULL);
217/* if(!(o = (TObjArray*)gFile->Get(Form("Moni%s", GetName())))){
218 AliWarning(Form("Monitor container Moni%s not available.", name));
219 return kFALSE;
220 }
221 fContainer = (TObjArray*)o->Clone("monitor");
222 fNRefFigures=AliTRDCalPID::kNMom;
223*/
224 // temporary until new calibration data are being produced
225 Histos();
226
227
228 if(!TFile::Open(Form("TRD.Calib%s.root", GetName()), "UPDATE")){
229 AliError(Form("File TRD.Calib%s.root corrupted", GetName()));
230 return kFALSE;
231 }
232 if(!(o = (TObjArray*)gFile->Get(Form("PDF_LQ")))) {
233 AliInfo("PDF container not available. Create.");
234 fPDF = new TObjArray(AliTRDCalPIDLQ::GetNrefs());
235 fPDF->SetOwner();fPDF->SetName("PDF_LQ");
236 } else fPDF = (TObjArray*)o->Clone("PDF_LQ");
237
238 return kTRUE;
239}
240
241#include "AliAnalysisManager.h"
242//________________________________________________________________________
243void AliTRDpidRefMakerLQ::UserExec(Option_t */*opt*/)
244{
245// Process pid info data array
246// The tasks have also access to the following containers which, for the time being, are not used
247// fTracks = dynamic_cast<TObjArray*>(GetInputData(1))
248// fEvent = dynamic_cast<AliTRDeventInfo*>(GetInputData(2))
249// fV0s = dynamic_cast<TObjArray*>(GetInputData(3))
250
251 if(!(fInfo = dynamic_cast<TObjArray*>(GetInputData(4)))){
252 Int_t ev((Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry());
253 AliDebug(3, Form("Missing pid info container in ev %d", ev));
254 return;
255 }
256
257 AliDebug(2, Form("ENTRIES pid[%d]\n", fInfo->GetEntriesFast()));
258 AliTRDpidInfo *pid(NULL);
259 const AliTRDpidInfo::AliTRDpidData *data(NULL);
260 Char_t s(-1);
261 for(Int_t itrk=fInfo->GetEntriesFast(); itrk--;){
262 if(!(pid=(AliTRDpidInfo*)fInfo->At(itrk))) continue;
263 if((s=pid->GetPID())<0) continue;
264 for(Int_t itrklt=pid->GetNtracklets();itrklt--;){
265 data=pid->GetData(itrklt);
266 Int_t ip(data->Momentum());
267
268 Double_t dedx[] = {0., 0.};
269 if(!AliTRDCalPIDLQ::CookdEdx(data->fdEdx, dedx)) continue;
270 ((TH2*)((TObjArray*)fContainer->At(ip))->At(s+Int_t(AliPID::kSPECIES)))->Fill(dedx[0], dedx[1]);
271 AliTRDCalPIDLQ::CookdEdx(data->fdEdx, dedx, kFALSE);
272 ((TH1*)((TObjArray*)fContainer->At(ip))->At(s))->Fill(dedx[0]+dedx[1]);
273 }
274 }
275}
276
277
278//________________________________________________________________________
279Bool_t AliTRDpidRefMakerLQ::PostProcess()
280{
281// Analyse merged dedx = f(p) distributions.
282// - select momentum - species bins
283// - rotate to principal components
284// - locally interpolate with TKDPDF
285// - save interpolation to monitoring histograms
286// - write pdf to file for loading to OCDB
287//
288
289 TCanvas *fMonitor(NULL);
290 // allocate working storage
291 const Int_t kWS(AliPID::kSPECIES*AliTRDCalPID::kNMom);
292 Float_t *data[2*kWS], *data1D[kWS];
293 for(Int_t i(0); i<kWS; i++){
294 data1D[i] = new Float_t[kMaxStat];
295 data[i] = new Float_t[kMaxStat];
296 data[kWS+i] = new Float_t[kMaxStat];
297 }
298 Int_t ndata[kWS]; memset(ndata, 0, kWS*sizeof(Int_t));
299
300 AliDebug(1, Form("Loading %d entries.", (Int_t)fData->GetEntries()));
301 for(Int_t itrk=0; itrk < fData->GetEntries(); itrk++){
302 if(!(fData->GetEntry(itrk))) continue;
303 Int_t sbin(fPIDdataArray->GetPID());
304 for(Int_t itrklt=fPIDdataArray->GetNtracklets(); itrklt--;){
305 Int_t pbin(fPIDdataArray->GetData(itrklt)->Momentum());
306
307 Double_t dedx[] = {0., 0.},
308 dedx1D[] = {0., 0.};
309 if(!AliTRDCalPIDLQ::CookdEdx(fPIDdataArray->GetData(itrklt)->fdEdx, dedx)) continue;
310 AliTRDCalPIDLQ::CookdEdx(fPIDdataArray->GetData(itrklt)->fdEdx, dedx1D, kFALSE);
311 Int_t idx=AliTRDCalPIDLQ::GetModelID(pbin,sbin);
312 if(ndata[idx]==kMaxStat) continue;
313
314 // store data
315 data1D[idx][ndata[idx]] = dedx1D[0];
316 data[idx][ndata[idx]] = dedx[0];
317 data[idx+kWS][ndata[idx]] = dedx[1];
318 ndata[idx]++;
319 }
320 }
321
322 TKDPDF *pdf(NULL);
323 Int_t in(0); Float_t par[6], *pp(NULL);
324 for(Int_t ip=0;ip<AliTRDCalPID::kNMom;ip++){
325 for(Int_t is=AliPID::kSPECIES; is--;) {
326 // estimate bucket statistics
327 Int_t idx(AliTRDCalPIDLQ::GetModelID(ip,is)),
328 nb(kMinBuckets), // number of buckets
329 ns((Int_t)(((Float_t)(ndata[idx]))/nb)); //statistics/bucket
330
331 AliDebug(2, Form("pBin[%d] sBin[%d] n[%d] ns[%d] nb[%d]", ip, is, ndata[idx], ns, nb));
332 if(ns<Int_t(kMinStat)){
333 AliWarning(Form("Not enough entries [%d] for %s[%d].", ndata[idx], AliPID::ParticleShortName(is), ip));
334 continue;
335 }
336
337 // build helper 1D PDF
338 pdf = new TKDPDF(ndata[idx], 1, ns, &data1D[idx]);
339 pdf->SetCOG(kFALSE);
340 pdf->SetWeights();
341 //pdf->SetStore();
342 pdf->SetAlpha(15.);
343 //pdf.GetStatus();
344 fPDF->AddAt(pdf, idx);
345 in=pdf->GetNTNodes(); pp=&par[0];
346 while(in--){
347 const TKDNodeInfo *nn = pdf->GetNodeInfo(in);
348 nn->GetCOG(pp);
349 Double_t p(par[0]),r,e;
350 pdf->Eval(&p,r,e,1);
351 }
352
353 // build helper 2D PDF
354 Float_t *ldata[2]={data[idx], data[kWS+idx]};
355 pdf = new TKDPDF(ndata[idx], 2, ns, ldata);
356 pdf->SetCOG(kFALSE);
357 pdf->SetWeights();
358 //pdf->SetStore();
359 pdf->SetAlpha(5.);
360 //pdf.GetStatus();
361 fPDF->AddAt(pdf, AliTRDCalPIDLQ::GetModelID(ip,is, 2));
362 in=pdf->GetNTNodes(); pp=&par[0];
363 while(in--){
364 const TKDNodeInfo *nn = pdf->GetNodeInfo(in);
365 nn->GetCOG(pp);
366 Double_t p[] = {par[0], par[1]}, r,e;
367 pdf->Eval(p,r,e,1);
368 }
369/*
370 Int_t nnodes = pdf.GetNTNodes(),
371 nside = Int_t(0.05*nnodes),
372 nzeros = 4*(nside+1);
373 printf("nnodes[%d] nside[%d] nzeros[%d]\n", nnodes, nside, nzeros);
374
375
376 // Build interpolator on the pdf skeleton
377 TKDInterpolator interpolator(2, nnodes+nzeros);
378 for(Int_t in=nnodes; in--;)
379 interpolator.SetNode(in, *pdf.GetNodeInfo(in));
380 TKDNodeInfo *nodes = new TKDNodeInfo[nzeros], *node = &nodes[0];
381 Float_t ax0min, ax0max, ax1min, ax1max;
382 pdf.GetRange(0, ax0min, ax0max); Float_t dx = (ax0max-ax0min)/nside;
383 pdf.GetRange(1, ax1min, ax1max); Float_t dy = (ax1max-ax1min)/nside;
384 printf("x=[%f %f] y[%f %f]\n", ax0min, ax0max, ax1min, ax1max);
385
386 Int_t jn = nnodes;
387 SetZeroes(&interpolator, node, nside, jn, ax0min, dx, ax1min, -dy, 'x');
388 SetZeroes(&interpolator, node, nside, jn, ax1min, dy, ax0max, dx, 'y');
389 SetZeroes(&interpolator, node, nside, jn, ax0max,-dx, ax1max, dy, 'x');
390 SetZeroes(&interpolator, node, nside, jn ,ax1max, -dy, ax0min, -dx, 'y');
391 delete [] nodes;
392 Int_t in=nnodes; Float_t par[6], *pp=&par[0];
393 while(in--){
394 const TKDNodeInfo *nn = interpolator.GetNodeInfo(in);
395 nn->GetCOG(pp);
396 //printf("evaluate for node[%d] @ [%f %f]\n",in, par[0], par[1]);
397 Double_t p[] = {par[0], par[1]}, r,e;
398 interpolator.Eval(p,r,e,1);
399 }
400*/
401
402 // visual on-line monitoring
403 if(HasOnlineMonitor()){
404 if(!fMonitor) fMonitor = new TCanvas("cc", "PDF 2D LQ", 500, 500);
405 pdf->DrawProjection();
406 fMonitor->Modified(); fMonitor->Update();
407 fMonitor->SaveAs(Form("pdf_%s%02d.gif", AliPID::ParticleShortName(is), ip));
408 }
409
410 //fContainer->ls();
411 // save a discretization of the PDF for result monitoring
412 Double_t rxy[]={0.,0.};
413 TH2 *h2s = (TH2D*)((TObjArray*)fContainer->At(ip))->At(AliPID::kSPECIES+is);
414 TAxis *ax = h2s->GetXaxis(), *ay = h2s->GetYaxis();
415 h2s->Clear();
416 for(int ix=1; ix<=ax->GetNbins(); ix++){
417 rxy[0] = ax->GetBinCenter(ix);
418 for(int iy=1; iy<=ay->GetNbins(); iy++){
419 rxy[1] = ay->GetBinCenter(iy);
420
421 Double_t rr,ee;
422 pdf->Eval(rxy, rr, ee, kFALSE);
423 if(rr<0. || ee/rr>.15) continue; // 15% relative error
424 //printf("x[%2d] x[%2d] r[%f] e[%f]\n", ix, iy, rr, ee);
425 h2s->SetBinContent(ix, iy, rr);
426 }
427 }
428
429 if(!(pdf=dynamic_cast<TKDPDF*>(fPDF->At(idx)))){
430 AliWarning(Form("Missing pdf for model id[%d]", idx));
431 continue;
432 }
433 TH1 *h1 = (TH1D*)((TObjArray*)fContainer->At(ip))->At(is);
434 ax = h1->GetXaxis();
435 h1->Clear();
436 for(int ix=1; ix<=ax->GetNbins(); ix++){
437 rxy[0] = ax->GetBinCenter(ix);
438
439 Double_t rr,ee;
440 pdf->Eval(rxy, rr, ee, kFALSE);
441 if(rr<0. || ee/rr>.15) continue; // 15% relative error
442 //printf("x[%2d] x[%2d] r[%f] e[%f]\n", ix, iy, rr, ee);
443 h1->SetBinContent(ix, rr);
444 }
445 }
446 }
447 for(Int_t i(0); i<kWS; i++){
448 delete [] data1D[i];
449 delete [] data[i];
450 delete [] data[i+kWS];
451 }
452 return kTRUE;
453}
454
455
456//__________________________________________________________________
457void AliTRDpidRefMakerLQ::SetZeroes(TKDInterpolator *interpolator, TKDNodeInfo *node, Int_t n, Int_t& idx, Float_t x, Float_t dx, Float_t y, Float_t dy, const Char_t opt)
458{
459// Set extra nodes to ensure boundary conditions
460
461 printf("SetZeroes(%c)\n", opt);
462 Float_t par[6], val[] = {0., 1.};
463 Int_t a[6];
464 if(opt=='x'){
465 a[0]=0; a[1]=1; a[2]=2; a[3]=3; a[4]=4; a[5]=5;
466 } else if(opt=='y'){
467 a[0]=1; a[1]=0; a[2]=4; a[3]=5; a[4]=2; a[5]=3;
468 } else return;
469 Float_t tmp;
470 par[a[1]] = y;
471 par[a[4]] = y; par[a[5]] = y+dy;
472 if(dy<0.){tmp=par[a[4]]; par[a[4]]=par[a[5]]; par[a[5]]=tmp;}
473 for(Int_t in=n; in--; node++, idx++, x+=dx){
474 par[a[0]] = x+.5*dx;
475 par[a[2]] = x; par[a[3]] = x+dx;
476 if(dx<0.){tmp=par[a[2]]; par[a[2]]=par[a[3]]; par[a[3]]=tmp;}
477 node->SetNode(2, par, val);
478 printf("\n\tnode[%d]\n", idx); node->Print();
479 interpolator->SetNode(idx, *node);
480 }
481 par[a[0]] = x;
482 par[a[2]] = x; par[a[3]] = x+dx;
483 if(dx<0.){tmp=par[a[2]]; par[a[2]]=par[a[3]]; par[a[3]]=tmp;}
484 node->SetNode(2, par, val);
485 printf("\n\tnode[%d]\n", idx); node->Print();
486 interpolator->SetNode(idx, *node);node++;idx++;
487}
488