1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 /* $Id: AliTRDpidRefMakerLQ.cxx 34163 2009-08-07 11:28:51Z cblume $ */
18 ///////////////////////////////////////////////////////////////////////////////
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)
26 // Alex Bercuci (A.Bercuci@gsi.de)
28 ///////////////////////////////////////////////////////////////////////////////
35 #include <TEventList.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"
51 ClassImp(AliTRDpidRefMakerLQ)
53 //__________________________________________________________________
54 AliTRDpidRefMakerLQ::AliTRDpidRefMakerLQ()
59 // AliTRDpidRefMakerLQ default constructor
61 SetNameTitle("refMakerLQ", "PID(LQ) Reference Maker");
64 //__________________________________________________________________
65 AliTRDpidRefMakerLQ::AliTRDpidRefMakerLQ(const char *name)
66 : AliTRDpidRefMaker(name, "PID(LQ) Reference Maker")
70 // AliTRDpidRefMakerLQ default constructor
72 DefineOutput(3, TObjArray::Class());
75 //__________________________________________________________________
76 AliTRDpidRefMakerLQ::~AliTRDpidRefMakerLQ()
79 // AliTRDCalPIDQRef destructor
82 //fPDF->Write("PDF_LQ", TObject::kSingleKey);
88 // //________________________________________________________________________
89 void AliTRDpidRefMakerLQ::UserCreateOutputObjects()
94 AliTRDpidRefMaker::UserCreateOutputObjects();
95 fContainer = Histos();
97 OpenFile(2, "RECREATE");
98 fPDF = new TObjArray(AliTRDCalPIDLQ::GetNrefs());
99 fPDF->SetOwner();fPDF->SetName("PDF_LQ");
103 //__________________________________________________________________
104 TObjArray* AliTRDpidRefMakerLQ::Histos()
109 AliError("Monitor histograms missing.");
112 Int_t n(fContainer->GetEntriesFast());
114 // save dE/dx references
116 for(Int_t ip=AliTRDCalPID::kNMom; ip--;){
117 TObjArray *arr = new TObjArray(2*AliPID::kSPECIES);
118 arr->SetName(Form("Pbin%02d", ip)); arr->SetOwner();
119 for(Int_t is=AliPID::kSPECIES; is--;) {
120 h = new TH1D(Form("h1%s%d", AliPID::ParticleShortName(is), ip), Form("1D %s @ Pbin[%d]", AliPID::ParticleName(is), ip), 50, 7., 12.);
121 h->GetXaxis()->SetTitle("log(dE/dx) [au]");
122 h->GetYaxis()->SetTitle("#");
123 h->SetLineColor(AliTRDCalPIDLQ::GetPartColor(is));
126 for(Int_t is=AliPID::kSPECIES; is--;) {
127 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.);
128 h->GetXaxis()->SetTitle("log(dE/dx_{am}) [au]");
129 h->GetYaxis()->SetTitle("log(dE/dx_{dr}) [au]");
130 h->GetZaxis()->SetTitle("#");
131 arr->AddAt(h, AliPID::kSPECIES+is);
133 fContainer->AddAt(arr, n+ip);
135 fNRefFigures+=AliTRDCalPID::kNMom;
140 //__________________________________________________________________
141 TObject* AliTRDpidRefMakerLQ::GetOCDBEntry(Option_t */*opt*/)
143 // Steer loading of OCDB LQ PID
145 if(gSystem->AccessPathName(Form("TRD.Calib%s.root", GetName()), kReadPermission)){
146 AliError(Form("File TRD.Calib%s.root not readable", GetName()));
149 AliTRDCalPIDLQ *cal = new AliTRDCalPIDLQ("pidLQ", "LQ TRD PID object");
150 cal->LoadReferences(Form("TRD.Calib%s.root", GetName()));
154 //__________________________________________________________________
155 Bool_t AliTRDpidRefMakerLQ::GetRefFigure(Int_t ifig)
157 // Steering reference picture
159 if(ifig<0 || ifig>=fNRefFigures){
160 AliError("Ref fig requested outside definition.");
164 AliWarning("Please provide a canvas to draw results.");
168 TObjArray *arr(NULL);TList *l(NULL);TH1 *h(NULL);
169 if(!(arr = (TObjArray*)fContainer->At(ifig))){
170 AliError(Form("PDF container @ pBin[%d] missing.", ifig));
173 gPad->Divide(5, 2, 1.e-5, 1.e-5);l=gPad->GetListOfPrimitives();
174 for(Int_t is=0; is<AliPID::kSPECIES; is++){
175 ((TVirtualPad*)l->At(is))->cd();
176 if(!(h=(TH1*)arr->At(is))){
177 AliError(Form("1D for %s @ pBin[%d] missing.", AliPID::ParticleShortName(is), ifig));
180 h->GetYaxis()->SetRangeUser(0., 1.2);
183 ((TVirtualPad*)l->At(AliPID::kSPECIES+is))->cd();
184 if(!(h=(TH1*)arr->At(AliPID::kSPECIES+is))){
185 AliError(Form("2D for %s @ pBin[%d] missing.", AliPID::ParticleShortName(is), ifig));
195 //________________________________________________________________________
196 Bool_t AliTRDpidRefMakerLQ::Load(const Char_t */*fname*/)
198 const Char_t *name("PIDrefMaker");
199 if(gSystem->AccessPathName(Form("TRD.Calib%s.root", name), kReadPermission)){
200 AliError(Form("File TRD.Calib%s.root not readable", name));
203 if(!TFile::Open(Form("TRD.Calib%s.root", name))){
204 AliError(Form("File TRD.Calib%s.root corrupted", name));
207 if (!(fData = dynamic_cast<TTree*>(gFile->Get(name)))) {
208 AliError(Form("Tree %s not available", name));
214 /* if(!(o = (TObjArray*)gFile->Get(Form("Moni%s", GetName())))){
215 AliWarning(Form("Monitor container Moni%s not available.", name));
218 fContainer = (TObjArray*)o->Clone("monitor");
219 fNRefFigures=AliTRDCalPID::kNMom;
221 // temporary until new calibration data are being produced
225 if(!TFile::Open(Form("TRD.Calib%s.root", GetName()), "UPDATE")){
226 AliError(Form("File TRD.Calib%s.root corrupted", GetName()));
229 if(!(o = (TObjArray*)gFile->Get(Form("PDF_LQ")))) {
230 AliInfo("PDF container not available. Create.");
231 fPDF = new TObjArray(AliTRDCalPIDLQ::GetNrefs());
232 fPDF->SetOwner();fPDF->SetName("PDF_LQ");
233 } else fPDF = (TObjArray*)o->Clone("PDF_LQ");
239 //________________________________________________________________________
240 void AliTRDpidRefMakerLQ::UserExec(Option_t */*opt*/)
242 // Load PID data into local data storage
244 AliTRDpidInfo *pid(NULL);
245 const AliTRDpidInfo::AliTRDpidData *data(NULL);
247 for(Int_t itrk=fInfo->GetEntriesFast(); itrk--;){
248 if(!(pid=(AliTRDpidInfo*)fInfo->At(itrk))) continue;
249 if((s=pid->GetPID())<0) continue;
250 for(Int_t itrklt=pid->GetNtracklets();itrklt--;){
251 data=pid->GetData(itrklt);
252 Int_t ip(data->Momentum());
255 Double_t dedx[] = {0., 0.};
256 if(!AliTRDCalPIDLQ::CookdEdx(data->fdEdx, dedx)) continue;
257 ((TH2*)((TObjArray*)fContainer->At(ip))->At(s+Int_t(AliPID::kSPECIES)))->Fill(dedx[0], dedx[1]);
258 AliTRDCalPIDLQ::CookdEdx(data->fdEdx, dedx, kFALSE);
259 ((TH1*)((TObjArray*)fContainer->At(ip))->At(s))->Fill(dedx[0]);
263 PostData(1, fContainer);
268 //________________________________________________________________________
269 Bool_t AliTRDpidRefMakerLQ::PostProcess()
271 // Analyse merged dedx = f(p) distributions.
272 // - select momentum - species bins
273 // - rotate to principal components
274 // - locally interpolate with TKDPDF
275 // - save interpolation to monitoring histograms
276 // - write pdf to file for loading to OCDB
279 TCanvas *fMonitor(NULL);
280 // allocate working storage
281 const Int_t kWS(AliPID::kSPECIES*AliTRDCalPID::kNMom);
282 Float_t *data[2*kWS], *data1D[kWS];
283 for(Int_t i(0); i<kWS; i++){
284 data1D[i] = new Float_t[kMaxStat];
285 data[i] = new Float_t[kMaxStat];
286 data[kWS+i] = new Float_t[kMaxStat];
288 Int_t ndata[kWS]; memset(ndata, 0, kWS*sizeof(Int_t));
290 AliDebug(1, Form("Loading data[%d]", fData->GetEntries()));
291 for(Int_t itrk=0; itrk < fData->GetEntries(); itrk++){
292 if(!(fData->GetEntry(itrk))) continue;
294 for(Int_t ily=fPIDdataArray->fNtracklets; ily--;){
295 Int_t pbin(fPIDdataArray->fData[ily].fPLbin & 0xf);
297 Double_t dedx[] = {0., 0.},
299 if(!AliTRDCalPIDLQ::CookdEdx(fPIDdataArray->fData[ily].fdEdx, dedx)) continue;
300 AliTRDCalPIDLQ::CookdEdx(fPIDdataArray->fData[ily].fdEdx, dedx1D, kFALSE);
301 Int_t idx=AliTRDCalPIDLQ::GetModelID(pbin,sbin);
302 if(ndata[idx]==kMaxStat) continue;
305 data1D[idx][ndata[idx]] = dedx1D[0];
306 data[idx][ndata[idx]] = dedx[0];
307 data[idx+kWS][ndata[idx]] = dedx[1];
313 Int_t in(0); Float_t par[6], *pp(NULL);
314 for(Int_t ip=0;ip<AliTRDCalPID::kNMom;ip++){
315 for(Int_t is=AliPID::kSPECIES; is--;) {
316 // estimate bucket statistics
317 Int_t idx(AliTRDCalPIDLQ::GetModelID(ip,is)),
318 nb(kMinBuckets), // number of buckets
319 ns((Int_t)(((Float_t)(ndata[idx]))/nb)); //statistics/bucket
321 AliDebug(2, Form("pBin[%d] sBin[%d] n[%d] ns[%d] nb[%d]", ip, is, ndata[idx], ns, nb));
322 if(ns<Int_t(kMinStat)){
323 AliWarning(Form("Not enough entries [%d] for %s[%d].", ndata[idx], AliPID::ParticleShortName(is), ip));
327 // build helper 1D PDF
328 pdf = new TKDPDF(ndata[idx], 1, ns, &data1D[idx]);
334 fPDF->AddAt(pdf, idx);
335 in=pdf->GetNTNodes(); pp=&par[0];
337 const TKDNodeInfo *nn = pdf->GetNodeInfo(in);
339 Double_t p(par[0]),r,e;
343 // build helper 2D PDF
344 Float_t *ldata[2]={data[idx], data[kWS+idx]};
345 pdf = new TKDPDF(ndata[idx], 2, ns, ldata);
351 fPDF->AddAt(pdf, AliTRDCalPIDLQ::GetModelID(ip,is, 2));
352 in=pdf->GetNTNodes(); pp=&par[0];
354 const TKDNodeInfo *nn = pdf->GetNodeInfo(in);
356 Double_t p[] = {par[0], par[1]}, r,e;
360 Int_t nnodes = pdf.GetNTNodes(),
361 nside = Int_t(0.05*nnodes),
362 nzeros = 4*(nside+1);
363 printf("nnodes[%d] nside[%d] nzeros[%d]\n", nnodes, nside, nzeros);
366 // Build interpolator on the pdf skeleton
367 TKDInterpolator interpolator(2, nnodes+nzeros);
368 for(Int_t in=nnodes; in--;)
369 interpolator.SetNode(in, *pdf.GetNodeInfo(in));
370 TKDNodeInfo *nodes = new TKDNodeInfo[nzeros], *node = &nodes[0];
371 Float_t ax0min, ax0max, ax1min, ax1max;
372 pdf.GetRange(0, ax0min, ax0max); Float_t dx = (ax0max-ax0min)/nside;
373 pdf.GetRange(1, ax1min, ax1max); Float_t dy = (ax1max-ax1min)/nside;
374 printf("x=[%f %f] y[%f %f]\n", ax0min, ax0max, ax1min, ax1max);
377 SetZeroes(&interpolator, node, nside, jn, ax0min, dx, ax1min, -dy, 'x');
378 SetZeroes(&interpolator, node, nside, jn, ax1min, dy, ax0max, dx, 'y');
379 SetZeroes(&interpolator, node, nside, jn, ax0max,-dx, ax1max, dy, 'x');
380 SetZeroes(&interpolator, node, nside, jn ,ax1max, -dy, ax0min, -dx, 'y');
382 Int_t in=nnodes; Float_t par[6], *pp=&par[0];
384 const TKDNodeInfo *nn = interpolator.GetNodeInfo(in);
386 //printf("evaluate for node[%d] @ [%f %f]\n",in, par[0], par[1]);
387 Double_t p[] = {par[0], par[1]}, r,e;
388 interpolator.Eval(p,r,e,1);
392 // visual on-line monitoring
393 if(HasOnlineMonitor()){
394 if(!fMonitor) fMonitor = new TCanvas("cc", "PDF 2D LQ", 500, 500);
395 pdf->DrawProjection();
396 fMonitor->Modified(); fMonitor->Update();
397 fMonitor->SaveAs(Form("pdf_%s%02d.gif", AliPID::ParticleShortName(is), ip));
401 // save a discretization of the PDF for result monitoring
402 Double_t rxy[]={0.,0.};
403 TH2 *h2s = (TH2D*)((TObjArray*)fContainer->At(ip))->At(AliPID::kSPECIES+is);
404 TAxis *ax = h2s->GetXaxis(), *ay = h2s->GetYaxis();
406 for(int ix=1; ix<=ax->GetNbins(); ix++){
407 rxy[0] = ax->GetBinCenter(ix);
408 for(int iy=1; iy<=ay->GetNbins(); iy++){
409 rxy[1] = ay->GetBinCenter(iy);
412 pdf->Eval(rxy, rr, ee, kFALSE);
413 if(rr<0. || ee/rr>.15) continue; // 15% relative error
414 //printf("x[%2d] x[%2d] r[%f] e[%f]\n", ix, iy, rr, ee);
415 h2s->SetBinContent(ix, iy, rr);
419 pdf=dynamic_cast<TKDPDF*>(fPDF->At(idx));
420 TH1 *h1 = (TH1D*)((TObjArray*)fContainer->At(ip))->At(is);
423 for(int ix=1; ix<=ax->GetNbins(); ix++){
424 rxy[0] = ax->GetBinCenter(ix);
427 pdf->Eval(rxy, rr, ee, kFALSE);
428 if(rr<0. || ee/rr>.15) continue; // 15% relative error
429 //printf("x[%2d] x[%2d] r[%f] e[%f]\n", ix, iy, rr, ee);
430 h1->SetBinContent(ix, rr);
434 for(Int_t i(0); i<kWS; i++){
437 delete [] data[i+kWS];
443 //__________________________________________________________________
444 void 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)
446 // Set extra nodes to ensure boundary conditions
448 printf("SetZeroes(%c)\n", opt);
449 Float_t par[6], val[] = {0., 1.};
452 a[0]=0; a[1]=1; a[2]=2; a[3]=3; a[4]=4; a[5]=5;
454 a[0]=1; a[1]=0; a[2]=4; a[3]=5; a[4]=2; a[5]=3;
458 par[a[4]] = y; par[a[5]] = y+dy;
459 if(dy<0.){tmp=par[a[4]]; par[a[4]]=par[a[5]]; par[a[5]]=tmp;}
460 for(Int_t in=n; in--; node++, idx++, x+=dx){
462 par[a[2]] = x; par[a[3]] = x+dx;
463 if(dx<0.){tmp=par[a[2]]; par[a[2]]=par[a[3]]; par[a[3]]=tmp;}
464 node->SetNode(2, par, val);
465 printf("\n\tnode[%d]\n", idx); node->Print();
466 interpolator->SetNode(idx, *node);
469 par[a[2]] = x; par[a[3]] = x+dx;
470 if(dx<0.){tmp=par[a[2]]; par[a[2]]=par[a[3]]; par[a[3]]=tmp;}
471 node->SetNode(2, par, val);
472 printf("\n\tnode[%d]\n", idx); node->Print();
473 interpolator->SetNode(idx, *node);node++;idx++;