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