]>
Commit | Line | Data |
---|---|---|
1ee39b3a | 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) | |
1ee39b3a | 24 | // |
25 | // Origin | |
26 | // Alex Bercuci (A.Bercuci@gsi.de) | |
27 | // | |
28 | /////////////////////////////////////////////////////////////////////////////// | |
29 | ||
b91fdd71 | 30 | #include <TSystem.h> |
1ee39b3a | 31 | #include <TROOT.h> |
b91fdd71 | 32 | #include <TFile.h> |
1ee39b3a | 33 | #include <TTree.h> |
34 | #include <TMath.h> | |
35 | #include <TEventList.h> | |
36 | #include <TH2D.h> | |
37 | #include <TH2I.h> | |
1ee39b3a | 38 | #include <TCanvas.h> |
1ee39b3a | 39 | |
40 | #include "AliLog.h" | |
da27d983 | 41 | #include "../STAT/TKDPDF.h" |
4826d5b1 | 42 | #include "../STAT/TKDInterpolator.h" |
1ee39b3a | 43 | #include "AliTRDpidRefMakerLQ.h" |
da27d983 | 44 | #include "Cal/AliTRDCalPID.h" |
45 | #include "Cal/AliTRDCalPIDLQ.h" | |
1ee39b3a | 46 | #include "AliTRDseedV1.h" |
47 | #include "AliTRDcalibDB.h" | |
48 | #include "AliTRDgeometry.h" | |
4826d5b1 | 49 | #include "info/AliTRDpidInfo.h" |
4860c9db | 50 | #include "AliAnalysisManager.h" |
1ee39b3a | 51 | |
52 | ClassImp(AliTRDpidRefMakerLQ) | |
53 | ||
54 | //__________________________________________________________________ | |
4826d5b1 | 55 | AliTRDpidRefMakerLQ::AliTRDpidRefMakerLQ() |
705f8b0a | 56 | : AliTRDpidRefMaker() |
a5a3321d | 57 | ,fPDF(NULL) |
1ee39b3a | 58 | { |
59 | // | |
60 | // AliTRDpidRefMakerLQ default constructor | |
61 | // | |
4fa7d600 | 62 | SetNameTitle("TRDrefMakerLQ", "PID(LQ) Reference Maker"); |
705f8b0a | 63 | } |
64 | ||
65 | //__________________________________________________________________ | |
66 | AliTRDpidRefMakerLQ::AliTRDpidRefMakerLQ(const char *name) | |
67 | : AliTRDpidRefMaker(name, "PID(LQ) Reference Maker") | |
a5a3321d | 68 | ,fPDF(NULL) |
705f8b0a | 69 | { |
70 | // | |
71 | // AliTRDpidRefMakerLQ default constructor | |
72 | // | |
73 | DefineOutput(3, TObjArray::Class()); | |
1ee39b3a | 74 | } |
75 | ||
76 | //__________________________________________________________________ | |
77 | AliTRDpidRefMakerLQ::~AliTRDpidRefMakerLQ() | |
78 | { | |
79 | // | |
80 | // AliTRDCalPIDQRef destructor | |
81 | // | |
4860c9db | 82 | if (AliAnalysisManager::GetAnalysisManager()->IsProofMode()) return; |
4826d5b1 | 83 | if(fPDF){ |
705f8b0a | 84 | //fPDF->Write("PDF_LQ", TObject::kSingleKey); |
4826d5b1 | 85 | fPDF->Delete(); |
86 | delete fPDF; | |
87 | } | |
1ee39b3a | 88 | } |
89 | ||
b91fdd71 | 90 | // //________________________________________________________________________ |
f8f46e4d | 91 | void AliTRDpidRefMakerLQ::UserCreateOutputObjects() |
1ee39b3a | 92 | { |
93 | // Create histograms | |
94 | // Called once | |
95 | ||
b91fdd71 | 96 | fContainer = Histos(); |
068e2c00 | 97 | PostData(1, fContainer); |
4826d5b1 | 98 | |
a5a3321d | 99 | //OpenFile(2, "RECREATE"); |
092fbd8a | 100 | fPDF = new TObjArray(AliTRDCalPIDLQ::GetNrefs()); |
101 | fPDF->SetOwner();fPDF->SetName("PDF_LQ"); | |
068e2c00 | 102 | PostData(3, fPDF); |
b91fdd71 | 103 | } |
104 | ||
105 | ||
106 | //__________________________________________________________________ | |
107 | TObjArray* AliTRDpidRefMakerLQ::Histos() | |
108 | { | |
109 | // Create histograms | |
110 | ||
a5a3321d | 111 | if(fContainer) return fContainer; |
112 | fContainer = new TObjArray(AliTRDCalPID::kNMom); | |
1ee39b3a | 113 | |
114 | // save dE/dx references | |
092fbd8a | 115 | TH1 *h(NULL); |
1ee39b3a | 116 | for(Int_t ip=AliTRDCalPID::kNMom; ip--;){ |
092fbd8a | 117 | TObjArray *arr = new TObjArray(2*AliPID::kSPECIES); |
b91fdd71 | 118 | arr->SetName(Form("Pbin%02d", ip)); arr->SetOwner(); |
1ee39b3a | 119 | for(Int_t is=AliPID::kSPECIES; is--;) { |
092fbd8a | 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); | |
1ee39b3a | 132 | } |
a5a3321d | 133 | fContainer->AddAt(arr, ip); |
1ee39b3a | 134 | } |
a5a3321d | 135 | fNRefFigures=AliTRDCalPID::kNMom; |
b91fdd71 | 136 | return fContainer; |
1ee39b3a | 137 | } |
138 | ||
1ee39b3a | 139 | |
140 | //__________________________________________________________________ | |
4826d5b1 | 141 | TObject* AliTRDpidRefMakerLQ::GetOCDBEntry(Option_t */*opt*/) |
1ee39b3a | 142 | { |
143 | // Steer loading of OCDB LQ PID | |
144 | ||
4826d5b1 | 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 | } | |
1ee39b3a | 149 | AliTRDCalPIDLQ *cal = new AliTRDCalPIDLQ("pidLQ", "LQ TRD PID object"); |
4826d5b1 | 150 | cal->LoadReferences(Form("TRD.Calib%s.root", GetName())); |
1ee39b3a | 151 | return cal; |
152 | } | |
153 | ||
154 | //__________________________________________________________________ | |
155 | Bool_t AliTRDpidRefMakerLQ::GetRefFigure(Int_t ifig) | |
156 | { | |
157 | // Steering reference picture | |
158 | ||
b91fdd71 | 159 | if(ifig<0 || ifig>=fNRefFigures){ |
1ee39b3a | 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 | ||
092fbd8a | 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)); | |
b91fdd71 | 178 | return kFALSE; |
179 | } | |
092fbd8a | 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)); | |
b91fdd71 | 186 | return kFALSE; |
187 | } | |
092fbd8a | 188 | h->Draw("cont4z"); |
b91fdd71 | 189 | } |
092fbd8a | 190 | |
b91fdd71 | 191 | return kTRUE; |
192 | } | |
193 | ||
194 | ||
195 | //________________________________________________________________________ | |
fe1d1beb | 196 | Bool_t AliTRDpidRefMakerLQ::Load(const Char_t *file, const Char_t *dir) |
b91fdd71 | 197 | { |
64d57299 | 198 | // Load tree with data in case of detached PostProcess processing. |
199 | ||
fe1d1beb | 200 | if(gSystem->AccessPathName(file, kReadPermission)){ |
201 | AliError(Form("File %s not readable", file)); | |
b91fdd71 | 202 | return kFALSE; |
203 | } | |
fe1d1beb | 204 | if(!TFile::Open(file)) { |
205 | AliError(Form("File %s corrupted", file)); | |
b91fdd71 | 206 | return kFALSE; |
207 | } | |
fe1d1beb | 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"); | |
b91fdd71 | 214 | return kFALSE; |
1ee39b3a | 215 | } |
b91fdd71 | 216 | LinkPIDdata(); |
1ee39b3a | 217 | |
b91fdd71 | 218 | TObjArray *o(NULL); |
092fbd8a | 219 | /* if(!(o = (TObjArray*)gFile->Get(Form("Moni%s", GetName())))){ |
b91fdd71 | 220 | AliWarning(Form("Monitor container Moni%s not available.", name)); |
221 | return kFALSE; | |
222 | } | |
223 | fContainer = (TObjArray*)o->Clone("monitor"); | |
4826d5b1 | 224 | fNRefFigures=AliTRDCalPID::kNMom; |
092fbd8a | 225 | */ |
226 | // temporary until new calibration data are being produced | |
227 | Histos(); | |
4826d5b1 | 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 | } | |
092fbd8a | 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"); | |
4826d5b1 | 239 | |
1ee39b3a | 240 | return kTRUE; |
241 | } | |
242 | ||
94b94be0 | 243 | #include "AliAnalysisManager.h" |
b91fdd71 | 244 | //________________________________________________________________________ |
f8f46e4d | 245 | void AliTRDpidRefMakerLQ::UserExec(Option_t */*opt*/) |
b91fdd71 | 246 | { |
94b94be0 | 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 | } | |
87c9fc00 | 258 | |
a5a3321d | 259 | AliDebug(2, Form("ENTRIES pid[%d]\n", fInfo->GetEntriesFast())); |
4826d5b1 | 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()); | |
4826d5b1 | 269 | |
270 | Double_t dedx[] = {0., 0.}; | |
271 | if(!AliTRDCalPIDLQ::CookdEdx(data->fdEdx, dedx)) continue; | |
092fbd8a | 272 | ((TH2*)((TObjArray*)fContainer->At(ip))->At(s+Int_t(AliPID::kSPECIES)))->Fill(dedx[0], dedx[1]); |
273 | AliTRDCalPIDLQ::CookdEdx(data->fdEdx, dedx, kFALSE); | |
a5a3321d | 274 | ((TH1*)((TObjArray*)fContainer->At(ip))->At(s))->Fill(dedx[0]+dedx[1]); |
4826d5b1 | 275 | } |
276 | } | |
b91fdd71 | 277 | } |
278 | ||
279 | ||
1ee39b3a | 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 | ||
674c7d26 | 291 | TCanvas *fMonitor(NULL); |
b91fdd71 | 292 | // allocate working storage |
674c7d26 | 293 | const Int_t kWS(AliPID::kSPECIES*AliTRDCalPID::kNMom); |
092fbd8a | 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 | } | |
674c7d26 | 300 | Int_t ndata[kWS]; memset(ndata, 0, kWS*sizeof(Int_t)); |
301 | ||
ca50a37e | 302 | AliDebug(1, Form("Loading %d entries.", (Int_t)fData->GetEntries())); |
674c7d26 | 303 | for(Int_t itrk=0; itrk < fData->GetEntries(); itrk++){ |
304 | if(!(fData->GetEntry(itrk))) continue; | |
a5a3321d | 305 | Int_t sbin(fPIDdataArray->GetPID()); |
306 | for(Int_t itrklt=fPIDdataArray->GetNtracklets(); itrklt--;){ | |
307 | Int_t pbin(fPIDdataArray->GetData(itrklt)->Momentum()); | |
674c7d26 | 308 | |
092fbd8a | 309 | Double_t dedx[] = {0., 0.}, |
310 | dedx1D[] = {0., 0.}; | |
a5a3321d | 311 | if(!AliTRDCalPIDLQ::CookdEdx(fPIDdataArray->GetData(itrklt)->fdEdx, dedx)) continue; |
312 | AliTRDCalPIDLQ::CookdEdx(fPIDdataArray->GetData(itrklt)->fdEdx, dedx1D, kFALSE); | |
674c7d26 | 313 | Int_t idx=AliTRDCalPIDLQ::GetModelID(pbin,sbin); |
314 | if(ndata[idx]==kMaxStat) continue; | |
315 | ||
316 | // store data | |
092fbd8a | 317 | data1D[idx][ndata[idx]] = dedx1D[0]; |
318 | data[idx][ndata[idx]] = dedx[0]; | |
674c7d26 | 319 | data[idx+kWS][ndata[idx]] = dedx[1]; |
320 | ndata[idx]++; | |
321 | } | |
322 | } | |
092fbd8a | 323 | |
324 | TKDPDF *pdf(NULL); | |
325 | Int_t in(0); Float_t par[6], *pp(NULL); | |
4826d5b1 | 326 | for(Int_t ip=0;ip<AliTRDCalPID::kNMom;ip++){ |
1ee39b3a | 327 | for(Int_t is=AliPID::kSPECIES; is--;) { |
1ee39b3a | 328 | // estimate bucket statistics |
674c7d26 | 329 | Int_t idx(AliTRDCalPIDLQ::GetModelID(ip,is)), |
330 | nb(kMinBuckets), // number of buckets | |
7fe4e88b | 331 | ns((Int_t)(((Float_t)(ndata[idx]))/nb)); //statistics/bucket |
1ee39b3a | 332 | |
674c7d26 | 333 | AliDebug(2, Form("pBin[%d] sBin[%d] n[%d] ns[%d] nb[%d]", ip, is, ndata[idx], ns, nb)); |
1ee39b3a | 334 | if(ns<Int_t(kMinStat)){ |
674c7d26 | 335 | AliWarning(Form("Not enough entries [%d] for %s[%d].", ndata[idx], AliPID::ParticleShortName(is), ip)); |
1ee39b3a | 336 | continue; |
337 | } | |
338 | ||
092fbd8a | 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 | |
674c7d26 | 356 | Float_t *ldata[2]={data[idx], data[kWS+idx]}; |
092fbd8a | 357 | pdf = new TKDPDF(ndata[idx], 2, ns, ldata); |
4826d5b1 | 358 | pdf->SetCOG(kFALSE); |
359 | pdf->SetWeights(); | |
674c7d26 | 360 | //pdf->SetStore(); |
4826d5b1 | 361 | pdf->SetAlpha(5.); |
362 | //pdf.GetStatus(); | |
092fbd8a | 363 | fPDF->AddAt(pdf, AliTRDCalPIDLQ::GetModelID(ip,is, 2)); |
364 | in=pdf->GetNTNodes(); pp=&par[0]; | |
4826d5b1 | 365 | while(in--){ |
366 | const TKDNodeInfo *nn = pdf->GetNodeInfo(in); | |
367 | nn->GetCOG(pp); | |
4826d5b1 | 368 | Double_t p[] = {par[0], par[1]}, r,e; |
369 | pdf->Eval(p,r,e,1); | |
1ee39b3a | 370 | } |
4826d5b1 | 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; | |
674c7d26 | 394 | Int_t in=nnodes; Float_t par[6], *pp=&par[0]; |
4826d5b1 | 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 | */ | |
092fbd8a | 403 | |
4826d5b1 | 404 | // visual on-line monitoring |
674c7d26 | 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 | } | |
1ee39b3a | 411 | |
4826d5b1 | 412 | //fContainer->ls(); |
1ee39b3a | 413 | // save a discretization of the PDF for result monitoring |
4826d5b1 | 414 | Double_t rxy[]={0.,0.}; |
092fbd8a | 415 | TH2 *h2s = (TH2D*)((TObjArray*)fContainer->At(ip))->At(AliPID::kSPECIES+is); |
1ee39b3a | 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; | |
4826d5b1 | 424 | pdf->Eval(rxy, rr, ee, kFALSE); |
1ee39b3a | 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 | } | |
092fbd8a | 430 | |
7fe4e88b | 431 | if(!(pdf=dynamic_cast<TKDPDF*>(fPDF->At(idx)))){ |
432 | AliWarning(Form("Missing pdf for model id[%d]", idx)); | |
433 | continue; | |
434 | } | |
092fbd8a | 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 | } | |
1ee39b3a | 447 | } |
448 | } | |
092fbd8a | 449 | for(Int_t i(0); i<kWS; i++){ |
450 | delete [] data1D[i]; | |
451 | delete [] data[i]; | |
452 | delete [] data[i+kWS]; | |
453 | } | |
4826d5b1 | 454 | return kTRUE; |
1ee39b3a | 455 | } |
456 | ||
457 | ||
4826d5b1 | 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 |