trigger names changed (Markus)
[u/mrichter/AliRoot.git] / PWG1 / TRD / AliTRDpidRefMakerLQ.cxx
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
51 ClassImp(AliTRDpidRefMakerLQ)
52
53 //__________________________________________________________________
54 AliTRDpidRefMakerLQ::AliTRDpidRefMakerLQ() 
55   : AliTRDpidRefMaker()
56   ,fPDF(NULL)
57 {
58   //
59   // AliTRDpidRefMakerLQ default constructor
60   //
61   SetNameTitle("TRDrefMakerLQ", "PID(LQ) Reference Maker");
62 }
63
64 //__________________________________________________________________
65 AliTRDpidRefMakerLQ::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 //__________________________________________________________________
76 AliTRDpidRefMakerLQ::~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 // //________________________________________________________________________
89 void 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 //__________________________________________________________________
105 TObjArray* 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 //__________________________________________________________________
139 TObject* 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 //__________________________________________________________________
153 Bool_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 //________________________________________________________________________
194 Bool_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
242 //________________________________________________________________________
243 void AliTRDpidRefMakerLQ::UserExec(Option_t */*opt*/)
244 {
245 // Load PID data into local data storage
246
247   if(!(fInfo   = dynamic_cast<TObjArray*>(GetInputData(3)))) return;
248
249   AliDebug(2, Form("ENTRIES pid[%d]\n", fInfo->GetEntriesFast()));
250   AliTRDpidInfo *pid(NULL);
251   const AliTRDpidInfo::AliTRDpidData *data(NULL);
252   Char_t s(-1);
253   for(Int_t itrk=fInfo->GetEntriesFast(); itrk--;){
254     if(!(pid=(AliTRDpidInfo*)fInfo->At(itrk))) continue;
255     if((s=pid->GetPID())<0) continue;
256     for(Int_t itrklt=pid->GetNtracklets();itrklt--;){
257       data=pid->GetData(itrklt);
258       Int_t ip(data->Momentum());
259
260       Double_t dedx[] = {0., 0.};
261       if(!AliTRDCalPIDLQ::CookdEdx(data->fdEdx, dedx)) continue;
262       ((TH2*)((TObjArray*)fContainer->At(ip))->At(s+Int_t(AliPID::kSPECIES)))->Fill(dedx[0], dedx[1]);
263       AliTRDCalPIDLQ::CookdEdx(data->fdEdx, dedx, kFALSE);
264       ((TH1*)((TObjArray*)fContainer->At(ip))->At(s))->Fill(dedx[0]+dedx[1]);
265     }
266   }
267 }
268
269
270 //________________________________________________________________________
271 Bool_t AliTRDpidRefMakerLQ::PostProcess()
272 {
273 // Analyse merged dedx = f(p) distributions.
274 //   - select momentum - species bins
275 //   - rotate to principal components
276 //   - locally interpolate with TKDPDF
277 //   - save interpolation to monitoring histograms
278 //   - write pdf to file for loading to OCDB
279 // 
280
281   TCanvas *fMonitor(NULL);
282   // allocate working storage
283   const Int_t kWS(AliPID::kSPECIES*AliTRDCalPID::kNMom);
284   Float_t *data[2*kWS], *data1D[kWS];
285   for(Int_t i(0); i<kWS; i++){ 
286     data1D[i]   = new Float_t[kMaxStat];
287     data[i]     = new Float_t[kMaxStat];
288     data[kWS+i] = new Float_t[kMaxStat];
289   }
290   Int_t ndata[kWS]; memset(ndata, 0, kWS*sizeof(Int_t));
291
292   AliDebug(1, Form("Loading %d entries.", (Int_t)fData->GetEntries()));
293   for(Int_t itrk=0; itrk < fData->GetEntries(); itrk++){
294     if(!(fData->GetEntry(itrk))) continue;
295     Int_t sbin(fPIDdataArray->GetPID());
296     for(Int_t itrklt=fPIDdataArray->GetNtracklets(); itrklt--;){
297       Int_t pbin(fPIDdataArray->GetData(itrklt)->Momentum());
298       
299       Double_t dedx[] = {0., 0.}, 
300                dedx1D[] = {0., 0.};
301       if(!AliTRDCalPIDLQ::CookdEdx(fPIDdataArray->GetData(itrklt)->fdEdx, dedx)) continue;
302       AliTRDCalPIDLQ::CookdEdx(fPIDdataArray->GetData(itrklt)->fdEdx, dedx1D, kFALSE);
303       Int_t idx=AliTRDCalPIDLQ::GetModelID(pbin,sbin);
304       if(ndata[idx]==kMaxStat) continue;
305
306       // store data
307       data1D[idx][ndata[idx]]   = dedx1D[0];
308       data[idx][ndata[idx]]     = dedx[0];
309       data[idx+kWS][ndata[idx]] = dedx[1];
310       ndata[idx]++;
311     }
312   }
313
314   TKDPDF *pdf(NULL);
315   Int_t in(0); Float_t par[6], *pp(NULL);
316   for(Int_t ip=0;ip<AliTRDCalPID::kNMom;ip++){ 
317     for(Int_t is=AliPID::kSPECIES; is--;) {
318       // estimate bucket statistics
319       Int_t idx(AliTRDCalPIDLQ::GetModelID(ip,is)),
320             nb(kMinBuckets), // number of buckets
321         ns((Int_t)(((Float_t)(ndata[idx]))/nb));    //statistics/bucket
322             
323       AliDebug(2, Form("pBin[%d] sBin[%d] n[%d] ns[%d] nb[%d]", ip, is, ndata[idx], ns, nb));
324       if(ns<Int_t(kMinStat)){
325         AliWarning(Form("Not enough entries [%d] for %s[%d].", ndata[idx], AliPID::ParticleShortName(is), ip));
326         continue;
327       }
328
329       // build helper 1D PDF
330       pdf = new TKDPDF(ndata[idx], 1, ns, &data1D[idx]);
331       pdf->SetCOG(kFALSE);
332       pdf->SetWeights();
333       //pdf->SetStore();
334       pdf->SetAlpha(15.);
335       //pdf.GetStatus();
336       fPDF->AddAt(pdf, idx);
337       in=pdf->GetNTNodes(); pp=&par[0];
338       while(in--){
339         const TKDNodeInfo *nn = pdf->GetNodeInfo(in);
340         nn->GetCOG(pp);
341         Double_t p(par[0]),r,e;
342         pdf->Eval(&p,r,e,1);
343       }
344
345       // build helper 2D PDF
346       Float_t *ldata[2]={data[idx], data[kWS+idx]};
347       pdf = new TKDPDF(ndata[idx], 2, ns, ldata);
348       pdf->SetCOG(kFALSE);
349       pdf->SetWeights();
350       //pdf->SetStore();
351       pdf->SetAlpha(5.);
352       //pdf.GetStatus();
353       fPDF->AddAt(pdf, AliTRDCalPIDLQ::GetModelID(ip,is, 2));
354       in=pdf->GetNTNodes(); pp=&par[0];
355       while(in--){
356         const TKDNodeInfo *nn = pdf->GetNodeInfo(in);
357         nn->GetCOG(pp);
358         Double_t p[] = {par[0], par[1]}, r,e;
359         pdf->Eval(p,r,e,1);
360       }
361 /*
362       Int_t nnodes = pdf.GetNTNodes(),
363             nside = Int_t(0.05*nnodes),
364             nzeros = 4*(nside+1);
365       printf("nnodes[%d] nside[%d] nzeros[%d]\n", nnodes, nside, nzeros);
366     
367
368       // Build interpolator on the pdf skeleton
369       TKDInterpolator interpolator(2, nnodes+nzeros); 
370       for(Int_t in=nnodes; in--;)
371         interpolator.SetNode(in, *pdf.GetNodeInfo(in));
372       TKDNodeInfo *nodes = new TKDNodeInfo[nzeros], *node = &nodes[0];
373       Float_t ax0min, ax0max, ax1min, ax1max;
374       pdf.GetRange(0, ax0min, ax0max); Float_t dx = (ax0max-ax0min)/nside;
375       pdf.GetRange(1, ax1min, ax1max); Float_t dy = (ax1max-ax1min)/nside;
376       printf("x=[%f %f] y[%f %f]\n", ax0min, ax0max, ax1min, ax1max);
377
378       Int_t jn = nnodes; 
379       SetZeroes(&interpolator, node, nside, jn, ax0min, dx, ax1min, -dy, 'x');
380       SetZeroes(&interpolator, node, nside, jn, ax1min, dy, ax0max, dx, 'y');
381       SetZeroes(&interpolator, node, nside, jn, ax0max,-dx, ax1max, dy, 'x');
382       SetZeroes(&interpolator, node, nside, jn ,ax1max, -dy, ax0min, -dx, 'y');
383       delete [] nodes;
384       Int_t in=nnodes; Float_t par[6], *pp=&par[0];
385       while(in--){
386         const TKDNodeInfo *nn = interpolator.GetNodeInfo(in);
387         nn->GetCOG(pp);
388         //printf("evaluate for node[%d] @ [%f %f]\n",in, par[0], par[1]);
389         Double_t p[] = {par[0], par[1]}, r,e;
390         interpolator.Eval(p,r,e,1);
391       }
392 */
393
394       // visual on-line monitoring
395       if(HasOnlineMonitor()){
396         if(!fMonitor) fMonitor = new TCanvas("cc", "PDF 2D LQ", 500, 500);
397         pdf->DrawProjection();
398         fMonitor->Modified(); fMonitor->Update(); 
399         fMonitor->SaveAs(Form("pdf_%s%02d.gif", AliPID::ParticleShortName(is), ip));
400       }
401
402       //fContainer->ls();
403       // save a discretization of the PDF for result monitoring
404       Double_t rxy[]={0.,0.};
405       TH2 *h2s = (TH2D*)((TObjArray*)fContainer->At(ip))->At(AliPID::kSPECIES+is);
406       TAxis *ax = h2s->GetXaxis(), *ay = h2s->GetYaxis();
407       h2s->Clear();
408       for(int ix=1; ix<=ax->GetNbins(); ix++){
409         rxy[0] = ax->GetBinCenter(ix);
410         for(int iy=1; iy<=ay->GetNbins(); iy++){
411           rxy[1] = ay->GetBinCenter(iy);
412       
413           Double_t rr,ee;
414           pdf->Eval(rxy, rr, ee, kFALSE);
415           if(rr<0. || ee/rr>.15) continue; // 15% relative error
416           //printf("x[%2d] x[%2d] r[%f] e[%f]\n", ix, iy, rr, ee);
417           h2s->SetBinContent(ix, iy, rr);
418         }
419       }
420
421       pdf=dynamic_cast<TKDPDF*>(fPDF->At(idx));
422       TH1 *h1 = (TH1D*)((TObjArray*)fContainer->At(ip))->At(is);
423       ax = h1->GetXaxis();
424       h1->Clear();
425       for(int ix=1; ix<=ax->GetNbins(); ix++){
426         rxy[0] = ax->GetBinCenter(ix);
427       
428         Double_t rr,ee;
429         pdf->Eval(rxy, rr, ee, kFALSE);
430         if(rr<0. || ee/rr>.15) continue; // 15% relative error
431         //printf("x[%2d] x[%2d] r[%f] e[%f]\n", ix, iy, rr, ee);
432         h1->SetBinContent(ix, rr);
433       }
434     }
435   }
436   for(Int_t i(0); i<kWS; i++){ 
437     delete [] data1D[i];
438     delete [] data[i];
439     delete [] data[i+kWS];
440   }
441   return kTRUE;
442 }
443
444
445 //__________________________________________________________________
446 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)
447 {
448 // Set extra nodes to ensure boundary conditions
449   
450   printf("SetZeroes(%c)\n", opt);
451   Float_t par[6], val[] = {0., 1.};
452   Int_t a[6];
453   if(opt=='x'){
454     a[0]=0; a[1]=1; a[2]=2; a[3]=3; a[4]=4; a[5]=5;
455   } else if(opt=='y'){
456     a[0]=1; a[1]=0; a[2]=4; a[3]=5; a[4]=2; a[5]=3;
457   } else return;
458   Float_t tmp;
459   par[a[1]] = y;
460   par[a[4]] = y; par[a[5]] = y+dy;
461   if(dy<0.){tmp=par[a[4]]; par[a[4]]=par[a[5]]; par[a[5]]=tmp;}
462   for(Int_t in=n; in--; node++, idx++, x+=dx){
463     par[a[0]] = x+.5*dx;
464     par[a[2]] = x;  par[a[3]] = x+dx;
465     if(dx<0.){tmp=par[a[2]]; par[a[2]]=par[a[3]]; par[a[3]]=tmp;}
466     node->SetNode(2, par, val);
467     printf("\n\tnode[%d]\n", idx); node->Print();
468     interpolator->SetNode(idx, *node);
469   }
470   par[a[0]] = x;
471   par[a[2]] = x;  par[a[3]] = x+dx;
472   if(dx<0.){tmp=par[a[2]]; par[a[2]]=par[a[3]]; par[a[3]]=tmp;}
473   node->SetNode(2, par, val);
474   printf("\n\tnode[%d]\n", idx); node->Print();
475   interpolator->SetNode(idx, *node);node++;idx++;
476 }
477