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()
55 : AliTRDpidRefMaker("PIDrefMakerLQ", "PID(LQ) Reference Maker")
59 // AliTRDpidRefMakerLQ default constructor
61 DefineOutput(2, TObjArray::Class());
64 //__________________________________________________________________
65 AliTRDpidRefMakerLQ::~AliTRDpidRefMakerLQ()
68 // AliTRDCalPIDQRef destructor
71 fPDF->Write("PDF_2DLQ", TObject::kSingleKey);
77 // //________________________________________________________________________
78 void AliTRDpidRefMakerLQ::CreateOutputObjects()
83 //OpenFile(0, "RECREATE");
84 fContainer = Histos();
86 OpenFile(2, "RECREATE");
87 fPDF = new TObjArray(AliTRDCalPID::kNMom*AliPID::kSPECIES);
88 fPDF->SetOwner();fPDF->SetName("PDF_2DLQ");
92 //__________________________________________________________________
93 TObjArray* AliTRDpidRefMakerLQ::Histos()
97 if(fContainer) return fContainer;
99 fContainer = new TObjArray(AliTRDCalPID::kNMom);
100 //fContainer->SetOwner(kTRUE);
101 //fContainer->SetName(Form("Moni%s", GetName()));
103 // save dE/dx references
105 for(Int_t ip=AliTRDCalPID::kNMom; ip--;){
106 TObjArray *arr = new TObjArray(AliPID::kSPECIES);
107 arr->SetName(Form("Pbin%02d", ip)); arr->SetOwner();
108 for(Int_t is=AliPID::kSPECIES; is--;) {
109 h2 = new TH2D(Form("h%s%d", AliPID::ParticleShortName(is), ip), Form("%s @ Pbin[%d]", AliPID::ParticleName(is), ip), 50, 7., 12., 50, 6.5, 11.);
110 h2->GetXaxis()->SetTitle("log(dE/dx_{am}) [au]");
111 h2->GetYaxis()->SetTitle("log(dE/dx_{dr}) [au]");
112 h2->GetZaxis()->SetTitle("#");
115 fContainer->AddAt(arr, ip);
117 fNRefFigures=AliTRDCalPID::kNMom;
122 //__________________________________________________________________
123 TObject* AliTRDpidRefMakerLQ::GetOCDBEntry(Option_t */*opt*/)
125 // Steer loading of OCDB LQ PID
127 if(gSystem->AccessPathName(Form("TRD.Calib%s.root", GetName()), kReadPermission)){
128 AliError(Form("File TRD.Calib%s.root not readable", GetName()));
131 AliTRDCalPIDLQ *cal = new AliTRDCalPIDLQ("pidLQ", "LQ TRD PID object");
132 cal->LoadReferences(Form("TRD.Calib%s.root", GetName()));
136 //__________________________________________________________________
137 Bool_t AliTRDpidRefMakerLQ::GetRefFigure(Int_t ifig)
139 // Steering reference picture
141 if(ifig<0 || ifig>=fNRefFigures){
142 AliError("Ref fig requested outside definition.");
146 AliWarning("Please provide a canvas to draw results.");
150 TObjArray *arr(NULL);TList *l(NULL);TH2 *h2(NULL);
153 if(!(h2=(TH2*)fContainer->At(0))){
154 AliError("Abundance Plot missing.");
160 if(!(arr = (TObjArray*)fContainer->At(ifig))){
161 AliError(Form("2D container @ pBin[%d] missing.", ifig-1));
164 gPad->Divide(3, 2, 1.e-5, 1.e-5);l=gPad->GetListOfPrimitives();
165 for(Int_t is=0; is<AliPID::kSPECIES; is++){
166 ((TVirtualPad*)l->At(is))->cd();
167 if(!(h2=(TH2*)arr->At(is))){
168 AliError(Form("2D for %s @ pBin[%d] missing.", AliPID::ParticleShortName(is), ifig-1));
178 //________________________________________________________________________
179 Bool_t AliTRDpidRefMakerLQ::Load(const Char_t */*fname*/)
181 const Char_t *name("PIDrefMaker");
182 if(gSystem->AccessPathName(Form("TRD.Calib%s.root", name), kReadPermission)){
183 AliError(Form("File TRD.Calib%s.root not readable", name));
186 if(!TFile::Open(Form("TRD.Calib%s.root", name))){
187 AliError(Form("File TRD.Calib%s.root corrupted", name));
190 if (!(fData = dynamic_cast<TTree*>(gFile->Get(name)))) {
191 AliError(Form("Tree %s not available", name));
197 if(!(o = (TObjArray*)gFile->Get(Form("Moni%s", GetName())))){
198 AliWarning(Form("Monitor container Moni%s not available.", name));
201 fContainer = (TObjArray*)o->Clone("monitor");
202 fNRefFigures=AliTRDCalPID::kNMom;
205 if(!TFile::Open(Form("TRD.Calib%s.root", GetName()), "UPDATE")){
206 AliError(Form("File TRD.Calib%s.root corrupted", GetName()));
209 if(!(o = (TObjArray*)gFile->Get(Form("PDF_2DLQ")))) {
210 AliInfo("PDF container PDF_2DLQ not available. Create.");
211 fPDF = new TObjArray(AliTRDCalPID::kNMom*AliPID::kSPECIES);
212 fPDF->SetOwner();fPDF->SetName("PDF_2DLQ");
213 } else fPDF = (TObjArray*)o->Clone("PDF_2DLQ");
219 //________________________________________________________________________
220 void AliTRDpidRefMakerLQ::Exec(Option_t */*opt*/)
222 // Load PID data into local data storage
224 AliTRDpidInfo *pid(NULL);
225 const AliTRDpidInfo::AliTRDpidData *data(NULL);
227 for(Int_t itrk=fInfo->GetEntriesFast(); itrk--;){
228 if(!(pid=(AliTRDpidInfo*)fInfo->At(itrk))) continue;
229 if((s=pid->GetPID())<0) continue;
230 for(Int_t itrklt=pid->GetNtracklets();itrklt--;){
231 data=pid->GetData(itrklt);
232 Int_t ip(data->Momentum());
235 Double_t dedx[] = {0., 0.};
236 if(!AliTRDCalPIDLQ::CookdEdx(data->fdEdx, dedx)) continue;
237 ((TH2*)((TObjArray*)fContainer->At(ip))->At(s))->Fill(dedx[0], dedx[1]);
241 PostData(0, fContainer);
246 //________________________________________________________________________
247 Bool_t AliTRDpidRefMakerLQ::PostProcess()
249 // Analyse merged dedx = f(p) distributions.
250 // - select momentum - species bins
251 // - rotate to principal components
252 // - locally interpolate with TKDPDF
253 // - save interpolation to monitoring histograms
254 // - write pdf to file for loading to OCDB
257 TCanvas *cc = new TCanvas("cc", "", 500, 500);
258 // allocate working storage
260 new Float_t[kMaxStat],
261 new Float_t[kMaxStat]};
262 for(Int_t ip=0;ip<AliTRDCalPID::kNMom;ip++){
263 for(Int_t is=AliPID::kSPECIES; is--;) {
264 Int_t n(0); // index of data
265 for(Int_t itrk=0; (itrk < fData->GetEntries()) && (n<kMaxStat); itrk++){
266 if(!(fData->GetEntry(itrk))) continue;
267 if(fPIDbin!=is) continue;
268 for(Int_t ily=fPIDdataArray->fNtracklets; ily--;){
269 if((fPIDdataArray->fData[ily].fPLbin & 0xf)!= ip) continue;
271 Double_t dedx[] = {0., 0.};
272 if(!AliTRDCalPIDLQ::CookdEdx(fPIDdataArray->fData[ily].fdEdx, dedx)) continue;
274 data[0][n] = dedx[0];
275 data[1][n] = dedx[1];
276 n++; if(n==kMaxStat) break;
280 // estimate bucket statistics
281 Int_t nb(kMinBuckets), // number of buckets
282 ns(Int_t(Float_t(n)/nb)); //statistics/bucket
284 // if(Float_t(n)/nb < 220.) ns = 200; // 7% stat error
285 // else if(Float_t(n)/nb < 420.) ns = 400; // 5% stat error
287 AliDebug(2, Form("pBin[%d] sBin[%d] n[%d] ns[%d] nb[%d]", ip, is, n, ns, nb));
288 if(ns<Int_t(kMinStat)){
289 AliWarning(Form("Not enough entries [%d] for %s[%d].", n, AliPID::ParticleShortName(is), ip));
294 TKDPDF *pdf = new TKDPDF(n, 2, ns, data);
300 fPDF->AddAt(pdf, AliTRDCalPIDLQ::GetModelID(ip,is));
301 Int_t in=pdf->GetNTNodes(); Float_t par[6], *pp=&par[0];
303 const TKDNodeInfo *nn = pdf->GetNodeInfo(in);
305 //printf("evaluate for node[%d] @ [%f %f]\n",in, par[0], par[1]);
306 Double_t p[] = {par[0], par[1]}, r,e;
310 Int_t nnodes = pdf.GetNTNodes(),
311 nside = Int_t(0.05*nnodes),
312 nzeros = 4*(nside+1);
313 printf("nnodes[%d] nside[%d] nzeros[%d]\n", nnodes, nside, nzeros);
316 // Build interpolator on the pdf skeleton
317 TKDInterpolator interpolator(2, nnodes+nzeros);
318 for(Int_t in=nnodes; in--;)
319 interpolator.SetNode(in, *pdf.GetNodeInfo(in));
320 TKDNodeInfo *nodes = new TKDNodeInfo[nzeros], *node = &nodes[0];
321 Float_t ax0min, ax0max, ax1min, ax1max;
322 pdf.GetRange(0, ax0min, ax0max); Float_t dx = (ax0max-ax0min)/nside;
323 pdf.GetRange(1, ax1min, ax1max); Float_t dy = (ax1max-ax1min)/nside;
324 printf("x=[%f %f] y[%f %f]\n", ax0min, ax0max, ax1min, ax1max);
327 SetZeroes(&interpolator, node, nside, jn, ax0min, dx, ax1min, -dy, 'x');
328 SetZeroes(&interpolator, node, nside, jn, ax1min, dy, ax0max, dx, 'y');
329 SetZeroes(&interpolator, node, nside, jn, ax0max,-dx, ax1max, dy, 'x');
330 SetZeroes(&interpolator, node, nside, jn ,ax1max, -dy, ax0min, -dx, 'y');
332 Int_t in=interpolator.GetNTNodes(); Float_t par[6], *pp=&par[0];
334 const TKDNodeInfo *nn = interpolator.GetNodeInfo(in);
336 //printf("evaluate for node[%d] @ [%f %f]\n",in, par[0], par[1]);
337 Double_t p[] = {par[0], par[1]}, r,e;
338 interpolator.Eval(p,r,e,1);
341 // visual on-line monitoring
342 pdf->DrawProjection();cc->Modified(); cc->Update(); cc->SaveAs(Form("pdf_%s%02d.gif", AliPID::ParticleShortName(is), ip));
343 cc->SaveAs(Form("%s_%s%02d.gif", GetName(), AliPID::ParticleShortName(is), ip));
347 // save a discretization of the PDF for result monitoring
348 Double_t rxy[]={0.,0.};
349 TH2 *h2s = (TH2D*)((TObjArray*)fContainer->At(ip))->At(is);
350 TAxis *ax = h2s->GetXaxis(), *ay = h2s->GetYaxis();
352 for(int ix=1; ix<=ax->GetNbins(); ix++){
353 rxy[0] = ax->GetBinCenter(ix);
354 for(int iy=1; iy<=ay->GetNbins(); iy++){
355 rxy[1] = ay->GetBinCenter(iy);
358 pdf->Eval(rxy, rr, ee, kFALSE);
359 if(rr<0. || ee/rr>.15) continue; // 15% relative error
360 //printf("x[%2d] x[%2d] r[%f] e[%f]\n", ix, iy, rr, ee);
361 h2s->SetBinContent(ix, iy, rr);
370 //__________________________________________________________________
371 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)
373 // Set extra nodes to ensure boundary conditions
375 printf("SetZeroes(%c)\n", opt);
376 Float_t par[6], val[] = {0., 1.};
379 a[0]=0; a[1]=1; a[2]=2; a[3]=3; a[4]=4; a[5]=5;
381 a[0]=1; a[1]=0; a[2]=4; a[3]=5; a[4]=2; a[5]=3;
385 par[a[4]] = y; par[a[5]] = y+dy;
386 if(dy<0.){tmp=par[a[4]]; par[a[4]]=par[a[5]]; par[a[5]]=tmp;}
387 for(Int_t in=n; in--; node++, idx++, x+=dx){
389 par[a[2]] = x; par[a[3]] = x+dx;
390 if(dx<0.){tmp=par[a[2]]; par[a[2]]=par[a[3]]; par[a[3]]=tmp;}
391 node->SetNode(2, par, val);
392 printf("\n\tnode[%d]\n", idx); node->Print();
393 interpolator->SetNode(idx, *node);
396 par[a[2]] = x; par[a[3]] = x+dx;
397 if(dx<0.){tmp=par[a[2]]; par[a[2]]=par[a[3]]; par[a[3]]=tmp;}
398 node->SetNode(2, par, val);
399 printf("\n\tnode[%d]\n", idx); node->Print();
400 interpolator->SetNode(idx, *node);node++;idx++;