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