prepare PID LQ ref maker for production
[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 "AliTRDpidRefMakerLQ.h"
43 #include "Cal/AliTRDCalPID.h"
44 #include "Cal/AliTRDCalPIDLQ.h"
45 #include "AliTRDseedV1.h"
46 #include "AliTRDcalibDB.h"
47 #include "AliTRDgeometry.h"
48
49 ClassImp(AliTRDpidRefMakerLQ)
50
51 //__________________________________________________________________
52 AliTRDpidRefMakerLQ::AliTRDpidRefMakerLQ() :
53   AliTRDpidRefMaker("PIDrefMakerLQ", "PID(LQ) Reference Maker")
54 {
55   //
56   // AliTRDpidRefMakerLQ default constructor
57   //
58 }
59
60 //__________________________________________________________________
61 AliTRDpidRefMakerLQ::~AliTRDpidRefMakerLQ()
62 {
63   //
64   // AliTRDCalPIDQRef destructor
65   //
66 }
67
68 // //________________________________________________________________________
69 void AliTRDpidRefMakerLQ::CreateOutputObjects()
70 {
71   // Create histograms
72   // Called once
73
74   //OpenFile(0, "RECREATE");
75   fContainer = Histos();
76 }
77
78
79 //__________________________________________________________________
80 TObjArray* AliTRDpidRefMakerLQ::Histos()
81 {
82   // Create histograms
83
84   if(fContainer) return fContainer;
85
86   fContainer = new TObjArray(AliTRDCalPID::kNMom);
87   //fContainer->SetOwner(kTRUE);
88   //fContainer->SetName(Form("Moni%s", GetName()));
89
90   // save dE/dx references
91   TH2 *h2 = 0x0;
92   for(Int_t ip=AliTRDCalPID::kNMom; ip--;){ 
93     TObjArray *arr = new TObjArray(AliPID::kSPECIES);
94     arr->SetName(Form("Pbin%02d", ip)); arr->SetOwner();
95     for(Int_t is=AliPID::kSPECIES; is--;) {
96       h2 = new TH2D(Form("h%s%d", AliPID::ParticleShortName(is), ip), Form("%s @ Pbin[%d]", AliPID::ParticleName(is), ip), 50, 5., 10., 50, 5., 10.);
97       h2->GetXaxis()->SetTitle("log(dE/dx_{am}) [au]");
98       h2->GetYaxis()->SetTitle("log(dE/dx_{dr}) [au]");
99       h2->GetZaxis()->SetTitle("#");
100       arr->AddAt(h2, is);
101     }
102     fContainer->AddAt(arr, ip);
103   }
104   fNRefFigures=AliTRDCalPID::kNMom;
105   return fContainer;
106 }
107
108
109 //__________________________________________________________________
110 TObject* AliTRDpidRefMakerLQ::GetOCDBEntry(Option_t *opt)
111 {
112 // Steer loading of OCDB LQ PID
113
114   TDirectoryFile *d = 0x0;
115   if(!TFile::Open(Form("TRD.Calib%s.root", GetName()))) return 0x0;
116   if(!(d=(TDirectoryFile*)gFile->Get(Form("PDF_%s", opt)))) return 0x0;
117   AliTRDCalPIDLQ *cal = new AliTRDCalPIDLQ("pidLQ", "LQ TRD PID object");
118   cal->LoadPDF(d);
119   return cal;
120 }
121
122 //__________________________________________________________________
123 Bool_t AliTRDpidRefMakerLQ::GetRefFigure(Int_t ifig)
124 {
125 // Steering reference picture
126
127   if(ifig<0 || ifig>=fNRefFigures){ 
128     AliError("Ref fig requested outside definition.");
129     return kFALSE;
130   }
131   if(!gPad){
132     AliWarning("Please provide a canvas to draw results.");
133     return kFALSE;
134   }
135
136   TObjArray *arr(NULL);TList *l(NULL);TH2 *h2(NULL);
137   switch(ifig){
138   case 0: // PDG plot
139     if(!(h2=(TH2*)fContainer->At(0))){
140       AliError("Abundance Plot missing.");
141       return kFALSE;
142     }
143     h2->Draw("cont4z");
144     break;
145   default:
146     if(!(arr = (TObjArray*)fContainer->At(ifig))){
147       AliError(Form("2D container @ pBin[%d] missing.", ifig-1));
148       return kFALSE;
149     }
150     gPad->Divide(3, 2, 1.e-5, 1.e-5);l=gPad->GetListOfPrimitives(); 
151     for(Int_t is=0; is<AliPID::kSPECIES; is++){
152       ((TVirtualPad*)l->At(is))->cd();
153       if(!(h2=(TH2*)arr->At(is))){
154         AliError(Form("2D for %s @ pBin[%d] missing.", AliPID::ParticleShortName(is), ifig-1));
155         return kFALSE;
156       }
157       h2->Draw("cont4z");
158     }
159   }
160   return kTRUE;
161 }
162
163
164 //________________________________________________________________________
165 Bool_t AliTRDpidRefMakerLQ::Load(const Char_t */*fname*/)
166 {
167   const Char_t *name("PIDrefMaker");
168   if(gSystem->AccessPathName(Form("TRD.Calib%s.root", name), kReadPermission)){
169     AliError(Form("File TRD.Calib%s.root not readable", name));
170     return kFALSE;
171   }
172   if(!TFile::Open(Form("TRD.Calib%s.root", name))){
173     AliError(Form("File TRD.Calib%s.root corrupted", name));
174     return kFALSE;
175   }
176   if (!(fData = dynamic_cast<TTree*>(gFile->Get(name)))) {
177     AliError(Form("Tree %s not available", name));
178     return kFALSE;
179   }
180   LinkPIDdata();
181
182   TObjArray *o(NULL);
183   if(!(o = (TObjArray*)gFile->Get(Form("Moni%s", name)))){
184     AliWarning(Form("Monitor container Moni%s not available.", name));
185     return kFALSE;
186   }
187   fContainer = (TObjArray*)o->Clone("monitor");
188   return kTRUE;
189 }
190
191 //________________________________________________________________________
192 void AliTRDpidRefMakerLQ::Exec(Option_t */*opt*/)
193 {
194 // Mock up function to load PID data into local data storage
195   AliInfo(Form("fInfo[%d]\n", fInfo->GetEntriesFast()));
196   PostData(0, fContainer);
197 }
198
199
200 //________________________________________________________________________
201 Bool_t AliTRDpidRefMakerLQ::PostProcess()
202 {
203 // Analyse merged dedx = f(p) distributions.
204 //   - select momentum - species bins
205 //   - rotate to principal components
206 //   - locally interpolate with TKDPDF
207 //   - save interpolation to monitoring histograms
208 //   - write pdf to file for loading to OCDB
209 // 
210
211   TDatime d;
212   TDirectoryFile *pdfs = new TDirectoryFile(Form("PDF_%d", d.GetDate()), "PDFs for LQ TRD-PID", "", gFile);
213   pdfs->Write();
214   AliDebug(2, Form("Data[%d]", fData->GetEntries()));
215   pdfs->cd();
216
217   //TCanvas *cc = new TCanvas("cc", "", 500, 500);
218   // allocate working storage
219   Float_t *data[] = {
220     new Float_t[kMaxStat], 
221     new Float_t[kMaxStat]};
222   for(Int_t ip=AliTRDCalPID::kNMom; ip--; ){ 
223     for(Int_t is=AliPID::kSPECIES; is--;) {
224       Int_t n(0); // index of data
225       for(Int_t itrk=0; (itrk < fData->GetEntries()) && (n<kMaxStat); itrk++){
226         if(!(fData->GetEntry(itrk))) continue;
227         if(fPIDbin!=is) continue;
228         for(Int_t ily=fPIDdataArray->fNtracklets; ily--;){
229           if((fPIDdataArray->fData[ily].fPLbin & 0xf)!= ip) continue;
230           
231           Float_t dedx[] = {0., 0.};
232           for(Int_t islice=AliTRDCalPID::kNSlicesNN; islice--;){
233             Int_t jslice = islice>kNN2LQtransition;
234             dedx[jslice]+=fPIDdataArray->fData[ily].fdEdx[islice];
235           }
236           
237           // check data integrity
238           if(dedx[0]<1.e-30) continue;
239           if(dedx[1]<1.e-30) continue;
240
241           // store data
242           data[0][n] = TMath::Log(dedx[0]);
243           data[1][n] = TMath::Log(dedx[1]);
244           n++; if(n==kMaxStat) break;
245         }
246       }
247
248       // estimate bucket statistics
249       Int_t nb(kMinBuckets), // number of buckets
250             ns(Int_t(Float_t(n)/nb));    //statistics/bucket
251             
252 // if(Float_t(n)/nb < 220.) ns = 200; // 7% stat error
253 //       else if(Float_t(n)/nb < 420.) ns = 400; // 5% stat error
254
255       AliDebug(2, Form("pBin[%d] sBin[%d] n[%d] ns[%d] nb[%d]", ip, is, n, ns, nb));
256       if(ns<Int_t(kMinStat)){
257         AliWarning(Form("Not enough entries [%d] for %s[%d].", n, AliPID::ParticleShortName(is), ip));
258         continue;
259       }
260
261       // build PDF
262       TKDPDF pdf(n, 2, ns, data);
263       pdf.SetCOG(kFALSE);
264       pdf.SetWeights();
265       pdf.SetStore();
266       pdf.SetAlpha(5.);
267       pdf.GetStatus();
268       Float_t *c, v, ve; Double_t r, e, rxy[2];
269       for(Int_t in=pdf.GetNTNodes(); in--;){
270         pdf.GetCOGPoint(in, c, v, ve);
271         rxy[0] = (Double_t)c[0];rxy[1] = (Double_t)c[1];
272         pdf.Eval(rxy, r, e, kTRUE);
273       }
274 //       // visual on-line monitoring
275 //       pdf.DrawProjection();cc->Modified(); cc->Update(); cc->SaveAs(Form("pdf_%s%02d.gif", AliPID::ParticleShortName(is), ip));
276 //       cc->SaveAs(Form("%s_%s%02d.gif", GetName(), AliPID::ParticleShortName(is), ip));
277
278       // save a discretization of the PDF for result monitoring
279       TH2 *h2s = (TH2D*)((TObjArray*)fContainer->At(ip))->At(is);
280       TAxis *ax = h2s->GetXaxis(), *ay = h2s->GetYaxis();
281       h2s->Clear();
282       for(int ix=1; ix<=ax->GetNbins(); ix++){
283         rxy[0] = ax->GetBinCenter(ix);
284         for(int iy=1; iy<=ay->GetNbins(); iy++){
285           rxy[1] = ay->GetBinCenter(iy);
286       
287           Double_t rr,ee;
288           pdf.Eval(rxy, rr, ee, kFALSE);
289           if(rr<0. || ee/rr>.15) continue; // 15% relative error
290           //printf("x[%2d] x[%2d] r[%f] e[%f]\n", ix, iy, rr, ee);
291           h2s->SetBinContent(ix, iy, rr);
292         }
293       }
294
295       // write results to output array
296       //pdf.GetStatus();
297       pdf.Write(Form("%s[%d]", AliPID::ParticleShortName(is), ip));
298     }
299   }
300   delete [] data[0]; delete [] data[1];
301   pdfs->Write();
302   //fCalib->Close(); delete fCalib;
303
304   return kTRUE; // testing protection
305 }
306
307