1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 /* $Id: AliTRDpidRefMakerLQ.cxx 34163 2009-08-07 11:28:51Z cblume $ */
18 ///////////////////////////////////////////////////////////////////////////////
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)
27 // Alex Bercuci (A.Bercuci@gsi.de)
29 ///////////////////////////////////////////////////////////////////////////////
35 #include <TEventList.h>
38 #include <TPrincipal.h>
40 #include <TLinearFitter.h>
46 #include "../../STAT/TKDPDF.h"
47 #include "AliTRDpidRefMakerLQ.h"
48 #include "../Cal/AliTRDCalPID.h"
49 #include "AliTRDseedV1.h"
50 #include "AliTRDcalibDB.h"
51 #include "AliTRDgeometry.h"
53 ClassImp(AliTRDpidRefMakerLQ)
55 //__________________________________________________________________
56 AliTRDpidRefMakerLQ::AliTRDpidRefMakerLQ()
57 :AliTRDpidRefMaker("PidRefMakerLQ", "PID(LQ) Reference Maker")
63 // AliTRDpidRefMakerLQ default constructor
67 //__________________________________________________________________
68 AliTRDpidRefMakerLQ::~AliTRDpidRefMakerLQ()
71 // AliTRDCalPIDQRef destructor
79 //________________________________________________________________________
80 void AliTRDpidRefMakerLQ::CreateOutputObjects()
85 AliTRDpidRefMaker::CreateOutputObjects();
87 // save dE/dx references
89 for(Int_t ip=AliTRDCalPID::kNMom; ip--;){
90 TObjArray *arr = new TObjArray(AliPID::kSPECIES);
91 arr->SetName(Form("Pbin%02d", ip));
92 for(Int_t is=AliPID::kSPECIES; is--;) {
93 h2 = new TH2I(Form("h%s%d", AliPID::ParticleShortName(is), ip), Form("%s ref. dEdx @ Pbin[%d]", AliPID::ParticleName(is), ip), 50, 5., 10., 50, 5., 10.);
94 h2->GetXaxis()->SetTitle("log(dE/dx_{am}) [au]");
95 h2->GetYaxis()->SetTitle("log(dE/dx_{dr}) [au]");
96 h2->GetZaxis()->SetTitle("#");
99 fContainer->AddAt(arr, 1+ip);
102 fData = new TTree(GetName(), Form("Reference data for %s", GetName()));
103 fData->Branch("s", &fSbin, "l/b");
104 fData->Branch("p", &fPbin, "p/b");
105 fData->Branch("dEdx" , fdEdx , Form("dEdx[%d]/F", GetNslices()));
109 //________________________________________________________________________
110 Float_t* AliTRDpidRefMakerLQ::CookdEdx(AliTRDseedV1 *trklt)
112 trklt->CookdEdx(AliTRDpidUtil::kLQslices);
113 const Float_t *dedx = trklt->GetdEdx();
114 if(dedx[0]+dedx[1] <= 0.) return 0x0;
115 if(dedx[2] <= 0.) return 0x0;
117 fdEdx[0] = TMath::Log(dedx[0]+dedx[1]);
118 fdEdx[1] = TMath::Log(dedx[2]);
122 #include "../Cal/AliTRDCalPIDLQ.h"
124 //__________________________________________________________________
125 TObject* AliTRDpidRefMakerLQ::GetOCDBEntry(Option_t *opt)
127 TDirectoryFile *d = 0x0;
128 if(!TFile::Open(Form("TRD.Calib%s.root", GetName()))) return 0x0;
129 if(!(d=(TDirectoryFile*)gFile->Get(Form("PDF_%s", opt)))) return 0x0;
130 AliTRDCalPIDLQ *cal = new AliTRDCalPIDLQ("pidLQ", "LQ TRD PID object");
135 //__________________________________________________________________
136 Bool_t AliTRDpidRefMakerLQ::GetRefFigure(Int_t ifig)
138 if(ifig<0 || ifig>AliTRDCalPID::kNMom-1){
139 AliError("Ref fig requested outside definition.");
143 AliError("No results processed.");
147 AliWarning("Please provide a canvas to draw results.");
151 TObjArray *arr = (TObjArray*)fResults->At(ifig);
152 gPad->Divide(3, 2, 1.e-5, 1.e-5);
153 TList *l=gPad->GetListOfPrimitives();
154 for(Int_t is=0; is<AliPID::kSPECIES; is++){
155 ((TVirtualPad*)l->At(is))->cd();
156 ((TH2*)arr->At(is))->Draw("cont4z");
163 //________________________________________________________________________
164 void AliTRDpidRefMakerLQ::Fill()
167 fSbin = TMath::LocMax(AliPID::kSPECIES, fPID);
169 fPbin = AliTRDpidUtil::GetMomentumBin(fP);
173 ((TH2*)fContainer->At(0))->Fill(fSbin, fPbin);
174 TH2* h2 = (TH2*)((TObjArray*)fContainer->At(1+fPbin))->At(fSbin);
175 h2->Fill(fdEdx[0], fdEdx[1]);
178 //________________________________________________________________________
179 Bool_t AliTRDpidRefMakerLQ::PostProcess()
181 // Analyse merged dedx = f(p) distributions.
182 // - select momentum - species bins
183 // - rotate to principal components
184 // - locally interpolate with TKDPDF
185 // - save interpolation to monitoring histograms
186 // - write pdf to file for loading to OCDB
190 TFile *fCalib = TFile::Open(Form("TRD.Calib%s.root", GetName()), "update");
191 fData = dynamic_cast<TTree*>(gFile->Get(GetName()));
193 AliError("Tree not available");
197 TDirectoryFile *pdfs = new TDirectoryFile(Form("PDF_%d", d.GetDate()), "PDFs for LQ TRD-PID", "", gFile);
199 AliDebug(2, Form("Data[%d]", fData->GetEntries()));
201 // save a volatile PDF representation
204 fResults = new TObjArray(AliTRDCalPID::kNMom);
205 for(Int_t ip=AliTRDCalPID::kNMom; ip--;){
206 TObjArray *arr = new TObjArray(AliPID::kSPECIES);
207 arr->SetName(Form("Pbin%02d", ip));
208 for(Int_t is=AliPID::kSPECIES; is--;) {
209 h2 = new TH2D(Form("i%s%d", AliPID::ParticleShortName(is), ip), Form("Interpolated %s dEdx @ Pbin[%d]", AliPID::ParticleName(is), ip), 50, -3., 3., 50, -3., 3.);
210 h2->GetXaxis()->SetTitle("log(dE/dx^{*}_{am}) [au]");
211 h2->GetYaxis()->SetTitle("log(dE/dx^{*}_{dr}) [au]");
212 h2->GetZaxis()->SetTitle("#");
215 fResults->AddAt(arr, ip);
219 TCanvas *cc = new TCanvas("cc", "", 500, 500);
221 Float_t *data[] = {0x0, 0x0};
222 TPrincipal principal(2, "ND");
223 for(Int_t ip=AliTRDCalPID::kNMom; ip--; ){
224 for(Int_t is=AliPID::kSPECIES; is--;) {
226 Int_t n = fData->Draw("dEdx[0]:dEdx[1]", Form("p==%d&&s==%d", ip, is), "goff");
228 // estimate bucket statistics
229 Int_t nb(kMinBuckets), // number of buckets
230 ns(Int_t(Float_t(n)/nb)); //statistics/bucket
232 // if(Float_t(n)/nb < 220.) ns = 200; // 7% stat error
233 // else if(Float_t(n)/nb < 420.) ns = 400; // 5% stat error
235 AliDebug(2, Form("pBin[%d] sBin[%d] N[%d] n[%d] nb[%d]", ip, is, n, ns, nb));
236 if(ns<Int_t(kMinStat)){
237 AliWarning(Form("Not enough entries [%d] for %s[%d].", ns, AliPID::ParticleShortName(is), ip));
241 data[0] = new Float_t[n];data[1] = new Float_t[n];
243 // fill and make principal
244 Double_t *v1 = fData->GetV1(), *v2 = fData->GetV2();
246 Double_t dedx[] = {v1[n], v2[n]};
247 principal.AddRow(dedx);
249 principal.MakePrincipals();
251 // // calculate covariance ellipse
252 // eValues = principal.GetEigenValues();
255 // rx = 3.5*sqrt((*eValues)[0]);
256 // ry = 3.5*sqrt((*eValues)[1]);
258 // rotate to principal axis
259 const Double_t *xx; Double_t rxy[2];
260 while((xx = principal.GetRow(++n))){
261 principal.X2P(xx, rxy);
262 data[0][n]=rxy[0]; data[1][n]=rxy[1];
266 TKDPDF pdf(n, 2, ns, data);
270 Float_t *c, v, ve; Double_t r, e;
271 for(Int_t in=pdf.GetNTNodes(); in--;){
272 pdf.GetCOGPoint(in, c, v, ve);
273 rxy[0] = (Double_t)c[0];rxy[1] = (Double_t)c[1];
274 pdf.Eval(rxy, r, e, kTRUE);
276 // visual on-line monitoring
277 pdf.DrawBins(0,1,-5.,5.,-8.,8.);cc->Modified(); cc->Update(); cc->SaveAs(Form("pdf_%s%02d.gif", AliPID::ParticleShortName(is), ip));
279 // save a discretization of the PDF for result monitoring
280 TH2 *h2s = (TH2D*)((TObjArray*)fResults->At(ip))->At(is);
281 TAxis *ax = h2s->GetXaxis(), *ay = h2s->GetYaxis();
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);
289 if(r<0. || e/r>.15) continue; // 15% relative error
290 //printf("x[%2d] x[%2d] r[%f] e[%f]\n", ix, iy, r, e);
291 h2s->SetBinContent(ix, iy, r);
295 // write results to output array
296 pdf.Write(Form("%s[%d]", AliPID::ParticleShortName(is), ip));
298 delete [] data[0]; delete [] data[1];
302 fCalib->Close(); delete fCalib;
304 return kTRUE; // testing protection