]>
Commit | Line | Data |
---|---|---|
f4277607 | 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$ */ | |
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 <TROOT.h> | |
32 | #include <TFile.h> | |
33 | #include <TTree.h> | |
34 | #include <TH2D.h> | |
35 | #include <TH2I.h> | |
36 | #include <TH3D.h> | |
37 | #include <TParticle.h> | |
38 | #include <TParticle.h> | |
39 | #include <TPrincipal.h> | |
40 | #include <TVector3.h> | |
41 | #include <TLinearFitter.h> | |
42 | #include <TVectorT.h> | |
43 | #include <TCanvas.h> | |
44 | #include <TEllipse.h> | |
45 | #include <TMarker.h> | |
46 | ||
47 | #include "AliLog.h" | |
48 | #include "AliPID.h" | |
49 | #include "AliESD.h" | |
50 | #include "AliRun.h" | |
51 | #include "AliRunLoader.h" | |
52 | #include "AliStack.h" | |
53 | #include "AliHeader.h" | |
54 | #include "AliGenEventHeader.h" | |
55 | #include "AliESDtrack.h" | |
56 | ||
57 | #include "AliTRDCalPIDRefMaker.h" | |
720a0a16 | 58 | #include "AliTRDCalPID.h" |
f4277607 | 59 | #include "AliTRDcalibDB.h" |
60 | #include "AliTRDgeometry.h" | |
61 | #include "AliTRDtrack.h" | |
62 | ||
63 | #include <vector> | |
64 | ||
65 | ClassImp(AliTRDCalPIDRefMaker) | |
66 | ||
67 | TLinearFitter *AliTRDCalPIDRefMaker::fFitter2D2 = 0x0; | |
68 | TLinearFitter *AliTRDCalPIDRefMaker::fFitter2D1 = 0x0; | |
69 | ||
70 | //__________________________________________________________________ | |
71 | AliTRDCalPIDRefMaker::AliTRDCalPIDRefMaker() | |
72 | :TObject() | |
73 | { | |
74 | // | |
75 | // AliTRDCalPIDRefMaker default constructor | |
76 | // | |
77 | ||
78 | // histogram settings | |
79 | for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){ | |
6c86e4be | 80 | fH2dEdx[ispec] = 0x0; |
f4277607 | 81 | fPrinc[ispec] = 0x0; |
82 | } | |
83 | } | |
84 | ||
85 | //__________________________________________________________________ | |
86 | AliTRDCalPIDRefMaker::AliTRDCalPIDRefMaker(const AliTRDCalPIDRefMaker &ref) | |
87 | :TObject() | |
88 | { | |
89 | // | |
90 | // AliTRDCalPIDRefMaker copy constructor | |
91 | // | |
92 | ||
93 | for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){ | |
6c86e4be | 94 | if(ref.fH2dEdx[ispec]){ |
95 | fH2dEdx[ispec] = new TH2D((TH2D&)*(ref.fH2dEdx[ispec])); | |
96 | } else fH2dEdx[ispec] = 0x0; | |
f4277607 | 97 | fPrinc[ispec] = 0x0; |
98 | } | |
99 | } | |
100 | ||
101 | //__________________________________________________________________ | |
102 | AliTRDCalPIDRefMaker::~AliTRDCalPIDRefMaker() | |
103 | { | |
104 | // | |
105 | // AliTRDCalPIDQRef destructor | |
106 | // | |
107 | ||
108 | for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){ | |
6c86e4be | 109 | if(fH2dEdx[ispec]) delete fH2dEdx[ispec]; |
f4277607 | 110 | if(fPrinc[ispec]) delete fPrinc[ispec]; |
111 | } | |
112 | if(fFitter2D1){ delete fFitter2D1; fFitter2D1 = 0x0;} | |
113 | if(fFitter2D2){ delete fFitter2D2; fFitter2D2 = 0x0;} | |
114 | ||
115 | } | |
116 | ||
117 | //__________________________________________________________________ | |
118 | AliTRDCalPIDRefMaker& AliTRDCalPIDRefMaker::operator=(const AliTRDCalPIDRefMaker &ref) | |
119 | { | |
120 | // | |
121 | // Assignment operator | |
122 | // | |
123 | ||
124 | for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){ | |
6c86e4be | 125 | if(ref.fH2dEdx[ispec]) (ref.fH2dEdx[ispec])->Copy(*fH2dEdx[ispec]); |
f4277607 | 126 | fPrinc[ispec] = 0x0; |
127 | } | |
128 | return *this; | |
129 | } | |
130 | ||
131 | ||
132 | //__________________________________________________________________ | |
133 | Bool_t AliTRDCalPIDRefMaker::BuildLQReferences(Char_t *File, Char_t *dir) | |
134 | { | |
135 | // Build, Fill and write to file the histograms used for PID. | |
136 | // The simulations are looked in the | |
137 | // directories with the general form Form("p%3.1f", momentum) | |
138 | // starting from dir (default .). Here momentum belongs to the list | |
6c86e4be | 139 | // of known momentum to PID (see trackMomentum). |
f4277607 | 140 | // The output histograms are |
141 | // written to the file "File" in cwd (default | |
142 | // TRDPIDHistograms.root). In order to build a DB entry | |
143 | // consider running $ALICE_ROOT/Cal/AliTRDpidCDB.C | |
144 | // | |
145 | // Author: | |
146 | // Alex Bercuci (A.Bercuci@gsi.de) | |
147 | ||
148 | Int_t partCode[AliPID::kSPECIES] = | |
149 | {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton}; | |
150 | ||
151 | // check and retrive number of directories in the production | |
152 | Int_t nBatches; | |
153 | if(!(nBatches = CheckProdDirTree(dir))) return kFALSE; | |
154 | ||
155 | ||
156 | // Number of Time bins | |
157 | Int_t nTimeBins; | |
158 | AliTRDcalibDB *calibration = AliTRDcalibDB::Instance(); | |
159 | if(!calibration){ | |
160 | AliError("No AliTRDcalibDB available."); | |
161 | return kFALSE; | |
162 | } else { | |
163 | nTimeBins = calibration->GetNumberOfTimeBins(); | |
164 | // if (calibration->GetRun() > -1) nTimeBins = calibration->GetNumberOfTimeBins(); | |
165 | // else { | |
166 | // AliError("No run number set."); | |
167 | // return kFALSE; | |
168 | // } | |
169 | } | |
170 | ||
171 | // Build PID reference/working objects | |
6c86e4be | 172 | fH1TB[0] = new TH1F("h1TB_0", "", nTimeBins, -.5, nTimeBins-.5); |
173 | fH1TB[0]->SetLineColor(4); | |
174 | fH1TB[1] = new TH1F("h1TB_1", "", nTimeBins, -.5, nTimeBins-.5); | |
175 | fH1TB[1]->SetLineColor(2); | |
f4277607 | 176 | for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) if(!fPrinc[ispec]) fPrinc[ispec] = new TPrincipal(2, "ND"); |
177 | ||
178 | // Print statistics header | |
179 | Int_t nPart[AliPID::kSPECIES], nTotPart; | |
180 | printf("P[GeV/c] "); | |
44dbae42 | 181 | for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) printf(" %s[%%] ", AliTRDCalPID::GetPartSymb(ispec)); |
f4277607 | 182 | printf("\n-----------------------------------------------\n"); |
183 | ||
720a0a16 | 184 | Float_t trackMomentum[AliTRDCalPID::kNMom]; |
185 | //Float_t trackSegLength[AliTRDCalPID::kNLength]; | |
186 | for(int i=0; i<AliTRDCalPID::kNMom; i++) trackMomentum[i] = AliTRDCalPID::GetMomentum(i); | |
f4277607 | 187 | AliRunLoader *fRunLoader = 0x0; |
188 | TFile *esdFile = 0x0; | |
189 | TTree *esdTree = 0x0; | |
190 | AliESD *esd = 0x0; | |
191 | AliESDtrack *esdTrack = 0x0; | |
192 | ||
193 | // | |
194 | // Momentum loop | |
720a0a16 | 195 | for (Int_t imom = 0; imom < AliTRDCalPID::kNMom; imom++) { |
f4277607 | 196 | Reset(); |
197 | ||
198 | // init statistics for momentum | |
199 | nTotPart = 0; | |
200 | for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) nPart[ispec] = 0; | |
201 | ||
202 | // loop over production directories | |
203 | for(Int_t ibatch = 0; ibatch<nBatches; ibatch++){ | |
204 | // open run loader and load gAlice, kinematics and header | |
6c86e4be | 205 | fRunLoader = AliRunLoader::Open(Form("%s/%3.1fGeV/%03d/galice.root", dir, trackMomentum[imom], ibatch)); |
f4277607 | 206 | if (!fRunLoader) { |
6c86e4be | 207 | AliError(Form("Getting run loader for momentum %3.1f GeV/c batch %03d failed.", trackMomentum[imom], ibatch)); |
f4277607 | 208 | return kFALSE; |
209 | } | |
6c86e4be | 210 | TString s; s.Form("%s/%3.1fGeV/%03d/", dir, trackMomentum[imom], ibatch); |
f4277607 | 211 | fRunLoader->SetDirName(s); |
212 | fRunLoader->LoadgAlice(); | |
213 | gAlice = fRunLoader->GetAliRun(); | |
214 | if (!gAlice) { | |
6c86e4be | 215 | AliError(Form("galice object not found for momentum %3.1f GeV/c batch %03d.", trackMomentum[imom], ibatch)); |
f4277607 | 216 | return kFALSE; |
217 | } | |
218 | fRunLoader->LoadKinematics(); | |
219 | fRunLoader->LoadHeader(); | |
220 | ||
221 | // open the ESD file | |
6c86e4be | 222 | esdFile = TFile::Open(Form("%s/%3.1fGeV/%03d/AliESDs.root", dir, trackMomentum[imom], ibatch)); |
f4277607 | 223 | if (!esdFile || esdFile->IsZombie()) { |
6c86e4be | 224 | AliError(Form("Opening ESD file failed for momentum %3.1f GeV/c batch %03d.", trackMomentum[imom], ibatch)); |
f4277607 | 225 | return kFALSE; |
226 | } | |
227 | esdTree = (TTree*)esdFile->Get("esdTree"); | |
228 | if (!esdTree) { | |
6c86e4be | 229 | AliError(Form("ESD tree not found for momentum %3.1f GeV/c batch %03d.", trackMomentum[imom], ibatch)); |
f4277607 | 230 | return kFALSE; |
231 | } | |
232 | esd = new AliESD; | |
233 | esdTree->SetBranchAddress("ESD", &esd); | |
234 | ||
235 | ||
236 | // Event loop | |
237 | for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) { | |
238 | ||
239 | // Load MC info | |
240 | fRunLoader->GetEvent(iEvent); | |
241 | AliStack* stack = gAlice->Stack(); | |
242 | TArrayF vertex(3); | |
243 | fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex); | |
244 | ||
245 | // Load event summary data | |
246 | esdTree->GetEvent(iEvent); | |
247 | if (!esd) { | |
6c86e4be | 248 | AliWarning(Form("ESD object not found for event %d. (@ momentum %3.1f GeV/c, batch %03d)", iEvent, trackMomentum[imom], ibatch)); |
f4277607 | 249 | continue; |
250 | } | |
251 | ||
252 | // Track loop | |
253 | for(Int_t iTrack=0; iTrack<esd->GetNumberOfTracks(); iTrack++){ | |
254 | esdTrack = esd->GetTrack(iTrack); | |
255 | ||
256 | //if(!AliTRDpidESD::CheckTrack(esdTrack)) continue; | |
257 | ||
258 | if((esdTrack->GetStatus() & AliESDtrack::kITSrefit) == 0) continue; | |
259 | if(esdTrack->GetConstrainedChi2() > 1E9) continue; | |
260 | if ((esdTrack->GetStatus() & AliESDtrack::kESDpid) == 0) continue; | |
261 | if (esdTrack->GetTRDsignal() == 0.) continue; | |
262 | ||
263 | // read MC info | |
264 | Int_t label = esdTrack->GetLabel(); | |
265 | if(label<0) continue; | |
266 | if (label > stack->GetNtrack()) continue; // background | |
267 | TParticle* particle = stack->Particle(label); | |
268 | if(!particle){ | |
6c86e4be | 269 | AliWarning(Form("Retriving particle with index %d from AliStack failed. [@ momentum %3.1f batch %03d event %d track %d]", label, trackMomentum[imom], ibatch, iEvent, iTrack)); |
f4277607 | 270 | continue; |
271 | } | |
272 | if(particle->Pt() < 1.E-3) continue; | |
273 | // if (TMath::Abs(particle->Eta()) > 0.3) continue; | |
274 | TVector3 dVertex(particle->Vx() - vertex[0], | |
275 | particle->Vy() - vertex[1], | |
276 | particle->Vz() - vertex[2]); | |
277 | if (dVertex.Mag() > 1.E-4){ | |
278 | //AliInfo(Form("Particle with index %d generated too far from vertex. Skip from analysis. Details follows. [@ event %d track %d]", label, iEvent, iTrack)); | |
279 | //particle->Print(); | |
280 | continue; | |
281 | } | |
282 | Int_t iGen = -1; | |
283 | for (Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) | |
284 | if(TMath::Abs(particle->GetPdgCode()) == partCode[ispec]){ | |
285 | iGen = ispec; | |
286 | break; | |
287 | } | |
288 | if(iGen<0) continue; | |
289 | ||
290 | nPart[iGen]++; nTotPart++; | |
291 | ||
6c86e4be | 292 | Float_t mom; |
293 | //Float_t length; | |
f4277607 | 294 | Double_t dedx[AliTRDtrack::kNslice], dEdx; |
295 | Int_t timebin; | |
296 | for (Int_t iPlane=0; iPlane<AliTRDgeometry::kNplan; iPlane++){ | |
297 | // read data for track segment | |
298 | for(int iSlice=0; iSlice<AliTRDtrack::kNslice; iSlice++) | |
299 | dedx[iSlice] = esdTrack->GetTRDsignals(iPlane, iSlice); | |
300 | dEdx = esdTrack->GetTRDsignals(iPlane, -1); | |
301 | timebin = esdTrack->GetTRDTimBin(iPlane); | |
302 | ||
303 | // check data | |
304 | if ((dEdx <= 0.) || (timebin <= -1.)) continue; | |
305 | ||
306 | // retrive kinematic info for this track segment | |
307 | //if(!AliTRDpidESD::RecalculateTrackSegmentKine(esdTrack, iPlane, mom, length)) continue; | |
308 | mom = esdTrack->GetOuterParam()->GetP(); | |
309 | ||
310 | // find segment length and momentum bin | |
311 | Int_t jmom = 1, refMom = -1; | |
720a0a16 | 312 | while(jmom<AliTRDCalPID::kNMom-1 && mom>trackMomentum[jmom]) jmom++; |
6c86e4be | 313 | if(TMath::Abs(trackMomentum[jmom-1] - mom) < trackMomentum[jmom-1] * .2) refMom = jmom-1; |
314 | else if(TMath::Abs(trackMomentum[jmom] - mom) < trackMomentum[jmom] * .2) refMom = jmom; | |
f4277607 | 315 | if(refMom<0){ |
6c86e4be | 316 | AliInfo(Form("Momentum at plane %d entrance not in momentum window. [@ momentum %3.1f batch %03d event %d track %d]", iPlane, trackMomentum[imom], ibatch, iEvent, iTrack)); |
f4277607 | 317 | continue; |
318 | } | |
720a0a16 | 319 | /*while(jleng<AliTRDCalPID::kNLength-1 && length>trackSegLength[jleng]) jleng++;*/ |
f4277607 | 320 | |
321 | // this track segment has fulfilled all requierments | |
322 | //nPlanePID++; | |
323 | ||
324 | if(dedx[0] > 0. && dedx[1] > 0.){ | |
325 | dedx[0] = log(dedx[0]); dedx[1] = log(dedx[1]); | |
326 | fPrinc[iGen]->AddRow(dedx); | |
327 | } | |
6c86e4be | 328 | fH1TB[(iGen>0)?1:0]->Fill(timebin); |
f4277607 | 329 | } // end plane loop |
330 | } // end track loop | |
331 | } // end events loop | |
332 | ||
333 | delete esd; esd = 0x0; | |
334 | esdFile->Close(); | |
335 | delete esdFile; esdFile = 0x0; | |
336 | ||
337 | fRunLoader->UnloadHeader(); | |
338 | fRunLoader->UnloadKinematics(); | |
339 | delete fRunLoader; fRunLoader = 0x0; | |
340 | } // end directory loop | |
341 | ||
342 | // use data to prepare references | |
343 | Prepare2D(); | |
344 | // save references | |
345 | SaveReferences(imom, File); | |
346 | ||
347 | ||
348 | // print momentum statistics | |
6c86e4be | 349 | printf(" %3.1f ", trackMomentum[imom]); |
f4277607 | 350 | for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) printf(" %5.2f ", 100.*nPart[ispec]/nTotPart); |
351 | printf("\n"); | |
352 | } // end momentum loop | |
353 | ||
354 | TFile *fSave = 0x0; | |
355 | TListIter it((TList*)gROOT->GetListOfFiles()); | |
356 | while((fSave=(TFile*)it.Next())) | |
357 | if(strcmp(File, fSave->GetName())==0) break; | |
358 | ||
359 | fSave->cd(); | |
360 | fSave->Close(); | |
361 | delete fSave; | |
362 | ||
363 | return kTRUE; | |
364 | } | |
365 | ||
366 | //__________________________________________________________________ | |
6c86e4be | 367 | Bool_t AliTRDCalPIDRefMaker::BuildNNReferences(Char_t* /*File*/, Char_t* /*dir*/) |
f4277607 | 368 | { |
369 | return kTRUE; | |
370 | } | |
371 | ||
372 | //__________________________________________________________________ | |
373 | Double_t AliTRDCalPIDRefMaker::Estimate2D2(TH2 *h, Float_t &x, Float_t &y) | |
374 | { | |
375 | // | |
376 | // Linear interpolation of data point with a parabolic expresion using | |
377 | // the logarithm of the bin content in a close neighbourhood. It is | |
378 | // assumed that the bin content of h is in number of events ! | |
379 | // | |
380 | // Observation: | |
381 | // This function have to be used when the statistics of each bin is | |
382 | // sufficient and all bins are populated. For cases of sparse data | |
383 | // please refere to Estimate2D1(). | |
384 | // | |
385 | // Author : Alex Bercuci (A.Bercuci@gsi.de) | |
386 | // | |
387 | ||
388 | if(!h){ | |
389 | AliErrorGeneral("AliTRDCalPIDRefMaker::Estimate2D2()", "No histogram defined."); | |
390 | return 0.; | |
391 | } | |
392 | ||
393 | TAxis *ax = h->GetXaxis(), *ay = h->GetYaxis(); | |
394 | Int_t binx = ax->FindBin(x); | |
395 | Int_t biny = ay->FindBin(y); | |
396 | Int_t nbinsx = ax->GetNbins(); | |
397 | Int_t nbinsy = ay->GetNbins(); | |
398 | Double_t p[2]; | |
399 | Double_t entries; | |
400 | ||
401 | // late construction of fitter | |
402 | if(!fFitter2D2) fFitter2D2 = new TLinearFitter(6, "1++x++y++x*x++y*y++x*y"); | |
403 | ||
404 | fFitter2D2->ClearPoints(); | |
405 | Int_t npoints=0; | |
406 | Int_t binx0, binx1, biny0, biny1; | |
407 | for(int bin=0; bin<5; bin++){ | |
408 | binx0 = TMath::Max(1, binx-bin); | |
409 | binx1 = TMath::Min(nbinsx, binx+bin); | |
410 | for(int ibin=binx0; ibin<=binx1; ibin++){ | |
411 | biny0 = TMath::Max(1, biny-bin); | |
412 | biny1 = TMath::Min(nbinsy, biny+bin); | |
413 | for(int jbin=biny0; jbin<=biny1; jbin++){ | |
414 | if(ibin != binx0 && ibin != binx1 && jbin != biny0 && jbin != biny1) continue; | |
415 | if((entries = h->GetBinContent(ibin, jbin)) == 0.) continue; | |
416 | p[0] = ax->GetBinCenter(ibin); | |
417 | p[1] = ay->GetBinCenter(jbin); | |
418 | fFitter2D2->AddPoint(p, log(entries), 1./sqrt(entries)); | |
419 | npoints++; | |
420 | } | |
421 | } | |
422 | if(npoints>=25) break; | |
423 | } | |
424 | if(fFitter2D2->Eval() == 1){ | |
425 | printf("<I2> x = %9.4f y = %9.4f\n", x, y); | |
426 | printf("\tbinx %d biny %d\n", binx, biny); | |
427 | printf("\tpoints %d\n", npoints); | |
428 | ||
429 | return 0.; | |
430 | } | |
431 | TVectorD vec(6); | |
432 | fFitter2D2->GetParameters(vec); | |
433 | Double_t result = vec[0] + x*vec[1] + y*vec[2] + x*x*vec[3] + y*y*vec[4] + x*y*vec[5]; | |
434 | return exp(result); | |
435 | ||
436 | } | |
437 | ||
438 | //__________________________________________________________________ | |
439 | Double_t AliTRDCalPIDRefMaker::Estimate2D1(TH2 *h, Float_t &x, Float_t &y | |
440 | , Float_t &dCT, Float_t &rmin | |
441 | , Float_t &rmax) | |
442 | { | |
443 | // | |
444 | // Linear interpolation of data point with a plane using | |
445 | // the logarithm of the bin content in the area defined by the | |
446 | // d(cos(phi)) and dr=(rmin, rmax). It is assumed that the bin content | |
447 | // of h is number of events ! | |
448 | // | |
449 | ||
450 | if(!h){ | |
451 | AliErrorGeneral("AliTRDCalPIDRefMaker::Estimate2D1()", "No histogram defined."); | |
452 | return 0.; | |
453 | } | |
454 | ||
455 | TAxis *ax = h->GetXaxis(), *ay = h->GetYaxis(); | |
456 | // Int_t binx = ax->FindBin(x); | |
457 | // Int_t biny = ay->FindBin(y); | |
458 | Int_t nbinsx = ax->GetNbins(); | |
459 | Int_t nbinsy = ay->GetNbins(); | |
460 | Double_t p[2]; | |
461 | Double_t entries; | |
462 | Double_t rxy = sqrt(x*x + y*y), rpxy; | |
463 | ||
464 | // late construction of fitter | |
465 | if(!fFitter2D1) fFitter2D1 = new TLinearFitter(3, "1++x++y"); | |
466 | ||
467 | fFitter2D1->ClearPoints(); | |
468 | Int_t npoints=0; | |
469 | for(int ibin=1; ibin<=nbinsx; ibin++){ | |
470 | for(int jbin=1; jbin<=nbinsy; jbin++){ | |
471 | if((entries = h->GetBinContent(ibin, jbin)) == 0.) continue; | |
472 | p[0] = ax->GetBinCenter(ibin); | |
473 | p[1] = ay->GetBinCenter(jbin); | |
474 | rpxy = sqrt(p[0]*p[0] + p[1]*p[1]); | |
475 | if((x*p[0] + y*p[1])/rxy/rpxy < dCT) continue; | |
476 | if(rpxy<rmin || rpxy > rmax) continue; | |
477 | ||
478 | fFitter2D1->AddPoint(p, log(entries), 1./sqrt(entries)); | |
479 | npoints++; | |
480 | } | |
481 | } | |
482 | if(npoints<15) return 0.; | |
483 | if(fFitter2D1->Eval() == 1){ | |
484 | printf("<O2> x = %9.4f y = %9.4f\n", x, y); | |
485 | printf("\tpoints %d\n", npoints); | |
486 | return 0.; | |
487 | } | |
488 | TVectorD vec(3); | |
489 | fFitter2D1->GetParameters(vec); | |
490 | Double_t result = vec[0] + x*vec[1] + y*vec[2]; | |
491 | return exp(result); | |
492 | } | |
493 | ||
494 | //__________________________________________________________________ | |
495 | // Double_t AliTRDCalPIDRefMaker::Estimate3D2(TH3 *h, Float_t &x, Float_t &y, Float_t &z) | |
496 | // { | |
497 | // // Author Alex Bercuci (A.Bercuci@gsi.de) | |
498 | // return 0.; | |
499 | // } | |
500 | ||
501 | ||
502 | ||
503 | ///////////// Private functions /////////////////////////////////// | |
504 | ||
505 | //__________________________________________________________________ | |
506 | Int_t AliTRDCalPIDRefMaker::CheckProdDirTree(Char_t *dir) | |
507 | { | |
508 | // Scan directory tree for momenta. Returns the smallest number of | |
509 | // batches found in all directories or 0 if one momentum is missing. | |
510 | ||
511 | const char *pwd = gSystem->pwd(); | |
512 | ||
513 | if(!gSystem->ChangeDirectory(dir)){ | |
514 | AliError(Form("Couldn't access production root directory %s.", dir)); | |
515 | return 0; | |
516 | } | |
517 | ||
518 | Int_t iDir; | |
519 | Int_t nDir = Int_t(1.E6); | |
720a0a16 | 520 | for(int imom=0; imom<AliTRDCalPID::kNMom; imom++){ |
521 | if(!gSystem->ChangeDirectory(Form("%3.1fGeV", AliTRDCalPID::GetMomentum(imom)))){ | |
522 | AliError(Form("Couldn't find data for momentum %3.1f GeV/c.", AliTRDCalPID::GetMomentum(imom))); | |
f4277607 | 523 | return 0; |
524 | } | |
525 | ||
526 | iDir = 0; | |
527 | while(gSystem->ChangeDirectory(Form("%03d", iDir))){ | |
528 | iDir++; | |
529 | gSystem->ChangeDirectory(".."); | |
530 | } | |
531 | if(iDir < nDir) nDir = iDir; | |
532 | gSystem->ChangeDirectory(dir); | |
533 | } | |
534 | ||
535 | gSystem->ChangeDirectory(pwd); | |
536 | ||
537 | return nDir; | |
538 | } | |
539 | ||
540 | ||
541 | //__________________________________________________________________ | |
542 | void AliTRDCalPIDRefMaker::Prepare2D() | |
543 | { | |
544 | // | |
545 | // Prepares the 2-dimensional reference histograms | |
546 | // | |
547 | ||
548 | // histogram settings | |
549 | Float_t xmin, xmax, ymin, ymax; | |
550 | Int_t nbinsx, nbinsy, nBinsSector, nSectors; | |
551 | const Int_t color[] = {4, 3, 2, 7, 6}; | |
552 | for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){ | |
553 | // check PCA data | |
554 | if(!fPrinc[ispec]){ | |
44dbae42 | 555 | AliError(Form("No data defined for %s.", AliTRDCalPID::GetPartName(ispec))); |
f4277607 | 556 | return; |
557 | } | |
558 | // build reference histograms | |
559 | nBinsSector = 10; | |
560 | nSectors = 20; | |
561 | nbinsx = nbinsy = nBinsSector * nSectors; | |
562 | xmin = ymin = 0.; | |
563 | xmax = 8000.; ymax = 6000.; | |
6c86e4be | 564 | if(!fH2dEdx[ispec]){ |
44dbae42 | 565 | fH2dEdx[ispec] = new TH2D(Form("h2%s", AliTRDCalPID::GetPartSymb(ispec)), "", nbinsx, xmin, xmax, nbinsy, ymin, ymax); |
6c86e4be | 566 | fH2dEdx[ispec]->SetLineColor(color[ispec]); |
f4277607 | 567 | } |
568 | } | |
569 | ||
570 | // build transformed rotated histograms | |
571 | nbinsx = nbinsy = 100; | |
572 | xmin = ymin = -6.; | |
573 | xmax = 10.; ymax = 6.; | |
574 | TH2I *hProj = 0x0; | |
575 | if(!(hProj = (TH2I*)gDirectory->Get("hProj"))) hProj = new TH2I("hProj", "", nbinsx, xmin, xmax, nbinsy, ymin, ymax); | |
576 | else hProj->Reset(); | |
577 | ||
578 | printf("Doing interpolation and invertion ... "); fflush(stdout); | |
579 | Bool_t kDebugPlot = kTRUE; | |
580 | //Bool_t kDebugPrint = kFALSE; | |
581 | ||
582 | TCanvas *c = 0x0; | |
583 | TEllipse *ellipse = 0x0; | |
584 | TMarker *mark = 0x0; | |
585 | if(kDebugPlot){ | |
586 | c=new TCanvas("c2", "Interpolation 2D", 10, 10, 500, 500); | |
587 | ellipse = new TEllipse(); | |
588 | ellipse->SetFillStyle(0); ellipse->SetLineColor(2); | |
589 | mark = new TMarker(); | |
590 | mark->SetMarkerColor(2); mark->SetMarkerSize(2); mark->SetMarkerStyle(2); | |
591 | } | |
592 | ||
593 | // define observable variables | |
594 | TAxis *ax, *ay; | |
595 | Double_t xy[2], lxy[2], rxy[2]; | |
596 | Double_t estimate, position; | |
597 | const TVectorD *eValues; | |
598 | ||
599 | // define radial sectors | |
600 | const Int_t nPhi = 36; | |
601 | const Float_t dPhi = TMath::TwoPi()/nPhi; | |
6c86e4be | 602 | //const Float_t dPhiRange = .1; |
f4277607 | 603 | Int_t nPoints[nPhi], nFitPoints, binStart, binStop; |
604 | TLinearFitter refsFitter[nPhi], refsLongFitter(6, "1++x++y++x*x++y*y++x*y"); | |
605 | Float_t fFitterRange[nPhi]; | |
606 | Bool_t kFitterStatus[nPhi]; | |
607 | for(int iphi=0; iphi<nPhi; iphi++){ | |
608 | refsFitter[iphi].SetDim(3); | |
609 | refsFitter[iphi].SetFormula("1++x++y");//++x*x++y*y++x*y"); | |
610 | fFitterRange[iphi] = .8; | |
611 | kFitterStatus[iphi] = kFALSE; | |
612 | } | |
613 | std::vector<UShort_t> storeX[nPhi], storeY[nPhi]; | |
614 | ||
615 | // define working variables | |
6c86e4be | 616 | Float_t x0, y0, rx, ry; |
617 | //Float_t rc, rmin, rmax, dr, dCT; | |
f4277607 | 618 | Double_t Phi, r; |
619 | Int_t iPhi; | |
620 | Double_t entries; | |
621 | for(int ispec=0; ispec<5; ispec++){ | |
622 | hProj->Reset(); | |
623 | for(int iphi=0; iphi<nPhi; iphi++){ | |
624 | nPoints[iphi] = 0; | |
625 | refsFitter[iphi].ClearPoints(); | |
626 | fFitterRange[iphi] = .8; | |
627 | kFitterStatus[iphi] = kFALSE; | |
628 | storeX[iphi].clear(); | |
629 | storeY[iphi].clear(); | |
630 | } | |
631 | ||
632 | // calculate covariance ellipse | |
633 | fPrinc[ispec]->MakePrincipals(); | |
634 | eValues = fPrinc[ispec]->GetEigenValues(); | |
635 | x0 = 0.; | |
636 | y0 = 0.; | |
637 | rx = 3.5*sqrt((*eValues)[0]); | |
638 | ry = 3.5*sqrt((*eValues)[1]); | |
639 | ||
640 | // rotate to principal axis | |
641 | Int_t irow = 0; | |
642 | const Double_t *xx; | |
643 | printf("filling for spec %d ...\n", ispec); | |
644 | while((xx = fPrinc[ispec]->GetRow(irow++))){ | |
645 | fPrinc[ispec]->X2P(xx, rxy); | |
646 | hProj->Fill(rxy[0], rxy[1]); | |
647 | } | |
648 | ||
649 | ||
650 | // // debug plot | |
651 | // if(kDebugPlot){ | |
652 | // hProj->Draw(); | |
653 | // ellipse->DrawEllipse(x0, y0, rx, ry, 0., 360., 0.); | |
654 | // mark->DrawMarker(x0, y0); | |
655 | // gPad->Modified(); gPad->Update(); | |
656 | // } | |
657 | ||
658 | // define radial sectors | |
659 | ax=hProj->GetXaxis(); | |
660 | ay=hProj->GetYaxis(); | |
661 | for(int ibin=1; ibin<=ax->GetNbins(); ibin++){ | |
662 | rxy[0] = ax->GetBinCenter(ibin); | |
663 | for(int jbin=1; jbin<=ay->GetNbins(); jbin++){ | |
664 | rxy[1] = ay->GetBinCenter(jbin); | |
665 | ||
666 | if((entries = hProj->GetBinContent(ibin, jbin)) == 0) continue; | |
667 | ||
668 | position = rxy[0]*rxy[0]/rx/rx + rxy[1]*rxy[1]/ry/ry; | |
669 | if(position < 1.) continue; | |
670 | ||
671 | r = sqrt(rxy[0]*rxy[0] + rxy[1]*rxy[1]); | |
672 | Phi = ((rxy[1] > 0.) ? 1. : -1.) * TMath::ACos(rxy[0]/r); // [-pi, pi] | |
673 | iPhi = nPhi/2 + Int_t(Phi/dPhi) - ((Phi/dPhi > 0.) ? 0 : 1); | |
674 | ||
675 | refsFitter[iPhi].AddPoint(rxy, log(entries), 1./sqrt(entries)); | |
676 | nPoints[iPhi]++; | |
677 | } | |
678 | } | |
679 | ||
680 | // do interpolation | |
6c86e4be | 681 | ax=fH2dEdx[ispec]->GetXaxis(); |
682 | ay=fH2dEdx[ispec]->GetYaxis(); | |
f4277607 | 683 | for(int ibin=1; ibin<=ax->GetNbins(); ibin++){ |
684 | xy[0] = ax->GetBinCenter(ibin); | |
685 | lxy[0] = log(xy[0]); | |
686 | for(int jbin=1; jbin<=ay->GetNbins(); jbin++){ | |
687 | xy[1] = ay->GetBinCenter(jbin); | |
688 | lxy[1] = log(xy[1]); | |
689 | ||
690 | // rotate to PCA | |
691 | fPrinc[ispec]->X2P(lxy, rxy); | |
692 | ||
693 | // calculate border of covariance ellipse | |
694 | position = rxy[0]*rxy[0]/rx/rx + rxy[1]*rxy[1]/ry/ry; | |
695 | ||
696 | // interpolation inside the covariance ellipse | |
697 | if(position < 1.){ | |
b3d727cf | 698 | Float_t xTemp = rxy[0]; |
699 | Float_t yTemp = rxy[1]; | |
700 | estimate = Estimate2D2((TH2 *) hProj, xTemp, yTemp); | |
701 | rxy[0] = xTemp; | |
702 | rxy[1] = yTemp; | |
703 | // estimate = Estimate2D2((TH2 *) hProj, (Float_t)rxy[0], (Float_t)rxy[1]); | |
6c86e4be | 704 | fH2dEdx[ispec]->SetBinContent(ibin, jbin, estimate/xy[0]/xy[1]); |
f4277607 | 705 | } else { // interpolation outside the covariance ellipse |
706 | r = sqrt(rxy[0]*rxy[0] + rxy[1]*rxy[1]); | |
707 | Phi = ((rxy[1] > 0.) ? 1. : -1.) * TMath::ACos(rxy[0]/r); // [-pi, pi] | |
708 | iPhi = nPhi/2 + Int_t(Phi/dPhi) - ((Phi/dPhi > 0.) ? 0 : 1); | |
709 | ||
710 | storeX[iPhi].push_back(ibin); | |
711 | storeY[iPhi].push_back(jbin); | |
712 | } | |
713 | } | |
714 | } | |
715 | ||
716 | // Fill outliers | |
717 | // Radial fit on transformed rotated | |
718 | TVectorD vec(3); | |
719 | Int_t xbin, ybin; | |
720 | for(int iphi=0; iphi<nPhi; iphi++){ | |
721 | Phi = iphi * dPhi - TMath::Pi(); | |
6ab776a4 | 722 | if(TMath::Abs(TMath::Abs(Phi)-TMath::Pi()) < 100.*TMath::DegToRad()) continue; |
f4277607 | 723 | |
724 | ||
725 | refsFitter[iphi].Eval(); | |
726 | refsFitter[iphi].GetParameters(vec); | |
6c86e4be | 727 | for(UInt_t ipoint=0; ipoint<storeX[iphi].size(); ipoint++){ |
f4277607 | 728 | xbin = storeX[iphi].at(ipoint); |
729 | ybin = storeY[iphi].at(ipoint); | |
730 | xy[0] = ax->GetBinCenter(xbin); lxy[0] = log(xy[0]); | |
731 | xy[1] = ay->GetBinCenter(ybin); lxy[1] = log(xy[1]); | |
732 | ||
733 | // rotate to PCA | |
734 | fPrinc[ispec]->X2P(lxy, rxy); | |
735 | estimate = exp(vec[0] + rxy[0]*vec[1] + rxy[1]*vec[2]);//+ rxy[0]*rxy[0]*vec[3]+ rxy[1]*rxy[1]*vec[4]+ rxy[0]*rxy[1]*vec[5]); | |
6c86e4be | 736 | fH2dEdx[ispec]->SetBinContent(xbin, ybin, estimate/xy[0]/xy[1]); |
f4277607 | 737 | } |
738 | } | |
739 | ||
740 | // Longitudinal fit on ref histo in y direction | |
741 | for(int isector=1; isector<nSectors; isector++){ | |
742 | // define sectors | |
743 | binStart = ((isector>0)?isector-1:isector)*nBinsSector + 1; | |
744 | binStop = (isector+1)*nBinsSector; | |
745 | ||
746 | nPoints[0] = 0; | |
747 | refsLongFitter.ClearPoints(); | |
748 | storeX[0].clear(); | |
749 | storeY[0].clear(); | |
750 | ||
751 | for(int ibin=1; ibin<=ax->GetNbins(); ibin++){ | |
752 | xy[0] = ax->GetBinCenter(ibin); | |
753 | for(int jbin=binStart; jbin<=binStop; jbin++){ | |
754 | xy[1] = ay->GetBinCenter(jbin); | |
6c86e4be | 755 | if((entries = fH2dEdx[ispec]->GetBinContent(ibin, jbin)) > 0.){ |
f4277607 | 756 | refsLongFitter.AddPoint(xy, log(entries), 1.e-7); |
757 | nPoints[0]++; | |
758 | } else { | |
759 | storeX[0].push_back(ibin); | |
760 | storeY[0].push_back(jbin); | |
761 | } | |
762 | } | |
763 | if(nPoints[0] >= ((isector==0)?60:420)) break; | |
764 | } | |
765 | ||
766 | ||
767 | if(refsLongFitter.Eval() == 1){ | |
768 | printf("Error in Y sector %d\n", isector); | |
769 | continue; | |
770 | } | |
771 | refsLongFitter.GetParameters(vec); | |
6c86e4be | 772 | for(UInt_t ipoint=0; ipoint<storeX[0].size(); ipoint++){ |
f4277607 | 773 | xbin = storeX[0].at(ipoint); |
774 | ybin = storeY[0].at(ipoint); | |
775 | xy[0] = ax->GetBinCenter(xbin); | |
776 | xy[1] = ay->GetBinCenter(ybin); | |
777 | ||
778 | estimate = vec[0] + xy[0]*vec[1] + xy[1]*vec[2] + xy[0]*xy[0]*vec[3]+ xy[1]*xy[1]*vec[4]+ xy[0]*xy[1]*vec[5]; | |
6c86e4be | 779 | fH2dEdx[ispec]->SetBinContent(xbin, ybin, exp(estimate)); |
f4277607 | 780 | } |
781 | } | |
782 | ||
783 | // Longitudinal fit on ref histo in x direction | |
784 | for(int isector=1; isector<nSectors; isector++){ | |
785 | binStart = ((isector>0)?isector-1:isector)*nBinsSector + 1; | |
786 | binStop = (isector+1)*nBinsSector; | |
787 | ||
788 | nPoints[0] = 0; | |
789 | nFitPoints = 0; | |
790 | refsLongFitter.ClearPoints(); | |
791 | storeX[0].clear(); | |
792 | storeY[0].clear(); | |
793 | ||
794 | for(int jbin=1; jbin<=ay->GetNbins(); jbin++){ | |
795 | xy[1] = ay->GetBinCenter(jbin); | |
796 | for(int ibin=binStart; ibin<=binStop; ibin++){ | |
797 | xy[0] = ax->GetBinCenter(ibin); | |
6c86e4be | 798 | if((entries = fH2dEdx[ispec]->GetBinContent(ibin, jbin)) > 0.){ |
f4277607 | 799 | refsLongFitter.AddPoint(xy, log(entries), 1.e-7); |
800 | nPoints[0]++; | |
801 | } else { | |
802 | storeX[0].push_back(ibin); | |
803 | storeY[0].push_back(jbin); | |
804 | nFitPoints++; | |
805 | } | |
806 | } | |
807 | if(nPoints[0] >= ((isector==0)?60:420)) break; | |
808 | } | |
809 | if(nFitPoints == 0) continue; | |
810 | ||
811 | ||
812 | if(refsLongFitter.Eval() == 1){ | |
813 | printf("Error in X sector %d\n", isector); | |
814 | continue; | |
815 | } | |
816 | refsLongFitter.GetParameters(vec); | |
6c86e4be | 817 | for(UInt_t ipoint=0; ipoint<storeX[0].size(); ipoint++){ |
f4277607 | 818 | xbin = storeX[0].at(ipoint); |
819 | ybin = storeY[0].at(ipoint); | |
820 | xy[0] = ax->GetBinCenter(xbin); | |
821 | xy[1] = ay->GetBinCenter(ybin); | |
822 | ||
823 | estimate = vec[0] + xy[0]*vec[1] + xy[1]*vec[2] + xy[0]*xy[0]*vec[3]+ xy[1]*xy[1]*vec[4]+ xy[0]*xy[1]*vec[5]; | |
6c86e4be | 824 | fH2dEdx[ispec]->SetBinContent(xbin, ybin, exp(estimate)); |
f4277607 | 825 | } |
826 | } | |
827 | ||
828 | ||
6c86e4be | 829 | fH2dEdx[ispec]->Scale(1./fH2dEdx[ispec]->Integral()); |
f4277607 | 830 | |
831 | // debug plot | |
832 | if(kDebugPlot){ | |
6c86e4be | 833 | fH2dEdx[ispec]->Draw("cont1"); gPad->SetLogz(); |
f4277607 | 834 | gPad->Modified(); gPad->Update(); |
835 | } | |
836 | } | |
837 | AliInfo("Finish interpolation."); | |
838 | } | |
839 | ||
840 | //__________________________________________________________________ | |
841 | void AliTRDCalPIDRefMaker::Reset() | |
842 | { | |
843 | // | |
844 | // Reset reference histograms | |
845 | // | |
846 | ||
847 | for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){ | |
6c86e4be | 848 | if(fH2dEdx[ispec]) fH2dEdx[ispec]->Reset(); |
f4277607 | 849 | fPrinc[ispec]->Clear(); |
850 | } | |
6c86e4be | 851 | if(fH1TB[0]) fH1TB[0]->Reset(); |
852 | if(fH1TB[1]) fH1TB[1]->Reset(); | |
f4277607 | 853 | } |
854 | ||
855 | //__________________________________________________________________ | |
856 | void AliTRDCalPIDRefMaker::SaveReferences(const Int_t mom, const char *fn) | |
857 | { | |
858 | // | |
859 | // Save the reference histograms | |
860 | // | |
861 | ||
862 | TFile *fSave = 0x0; | |
863 | TListIter it((TList*)gROOT->GetListOfFiles()); | |
864 | Bool_t kFOUND = kFALSE; | |
865 | TDirectory *pwd = gDirectory; | |
866 | while((fSave=(TFile*)it.Next())) | |
867 | if(strcmp(fn, fSave->GetName())==0){ | |
868 | kFOUND = kTRUE; | |
869 | break; | |
870 | } | |
871 | if(!kFOUND) fSave = TFile::Open(fn, "RECREATE"); | |
872 | fSave->cd(); | |
873 | ||
720a0a16 | 874 | Float_t fmom = AliTRDCalPID::GetMomentum(mom); |
f4277607 | 875 | |
876 | // save dE/dx references | |
877 | TH2 *h2 = 0x0; | |
878 | for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){ | |
44dbae42 | 879 | h2 = (TH2D*)fH2dEdx[ispec]->Clone(Form("h2dEdx%s%d", AliTRDCalPID::GetPartSymb(ispec), mom)); |
880 | h2->SetTitle(Form("2D dEdx for particle %s @ %d", AliTRDCalPID::GetPartName(ispec), mom)); | |
f4277607 | 881 | h2->GetXaxis()->SetTitle("dE/dx_{TRD}^{amplif} [au]"); |
882 | h2->GetYaxis()->SetTitle("dE/dx_{TRD}^{drift} [au]"); | |
883 | h2->GetZaxis()->SetTitle("Entries"); | |
884 | h2->Write(); | |
885 | } | |
886 | ||
887 | // save maximum time bin references | |
888 | TH1 *h1 = 0x0; | |
6c86e4be | 889 | h1 = (TH1F*)fH1TB[0]->Clone(Form("h1MaxTBEL%02d", mom)); |
f4277607 | 890 | h1->SetTitle(Form("Maximum Time Bin distribution for electrons @ %4.1f GeV", fmom)); |
891 | h1->GetXaxis()->SetTitle("time [100 ns]"); | |
892 | h1->GetYaxis()->SetTitle("Probability"); | |
893 | h1->Write(); | |
894 | ||
6c86e4be | 895 | h1 = (TH1F*)fH1TB[1]->Clone(Form("h1MaxTBPI%02d", mom)); |
f4277607 | 896 | h1->SetTitle(Form("Maximum Time Bin distribution for pions @ %4.1f GeV", fmom)); |
897 | h1->GetXaxis()->SetTitle("time [100 ns]"); | |
898 | h1->GetYaxis()->SetTitle("Probability"); | |
899 | h1->Write(); | |
900 | ||
901 | ||
902 | pwd->cd(); | |
903 | } | |
904 |