]>
Commit | Line | Data |
---|---|---|
1a18cb2b | 1 | #include <TStyle.h> |
2 | #include <TROOT.h> | |
3 | #include <TAxis.h> | |
4 | #include <TF1.h> | |
5 | #include <TFile.h> | |
6 | #include <TH1.h> | |
7 | #include <THn.h> | |
8 | #include <TObjArray.h> | |
9 | #include <TObject.h> | |
10 | #include <TString.h> | |
03de6136 | 11 | #include <TVectorT.h> |
12 | #include <TCanvas.h> | |
13 | #include <TProfile2D.h> | |
1a18cb2b | 14 | #include <TTreeStream.h> |
15 | ||
16 | /* | |
17 | ||
18 | .L $ALICE_ROOT/TPC/Upgrade/macros/AnaDelta.C+g | |
19 | ||
20 | ||
21 | */ | |
03de6136 | 22 | TVectorD* MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax); |
23 | TVectorD* MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax); | |
24 | TVectorD* MakeArbitraryBinning(const char* bins); | |
25 | ||
26 | void DumpHn(THn *hn, TTreeSRedirector &stream); | |
1a18cb2b | 27 | |
28 | void AnaDelta(TString file, TString outDir=".") | |
29 | { | |
30 | // | |
31 | // | |
32 | // | |
33 | ||
34 | gStyle->SetOptFit(); | |
35 | ||
36 | TTreeSRedirector stream(Form("%s/deltas.root",outDir.Data())); | |
37 | gROOT->cd(); | |
38 | ||
39 | TFile f(file); | |
40 | THn *hn=(THn*)f.Get("hn"); | |
41 | ||
03de6136 | 42 | DumpHn(hn, stream); |
43 | ||
44 | delete hn; | |
45 | } | |
46 | ||
47 | ||
48 | void AnaDeltaTree(TString file, TString outFile="deltas_tree.root") | |
49 | { | |
50 | TFile f(file); | |
51 | gROOT->cd(); | |
52 | TTree *t = (TTree*)f.Get("delta"); | |
53 | Float_t soneOverPt=0.; | |
54 | Float_t radius=0.; | |
55 | Float_t trackPhi=0.; | |
56 | Float_t trackY=0.; | |
57 | Float_t trackZ=0.; | |
58 | Float_t resRphi=0.; | |
59 | Double_t trackRes=0.; | |
60 | Float_t pointY=0.; | |
61 | Float_t pointZ=0.; | |
62 | Short_t npTRD=0.; | |
63 | Short_t event=0.; | |
64 | ||
65 | t->SetBranchAddress("soneOverPt" , &soneOverPt); | |
66 | t->SetBranchAddress("r" , &radius); | |
67 | t->SetBranchAddress("trackPhi" , &trackPhi); | |
68 | t->SetBranchAddress("trackY" , &trackY); | |
69 | t->SetBranchAddress("trackZ" , &trackZ); | |
70 | t->SetBranchAddress("resRphi" , &resRphi); | |
71 | t->SetBranchAddress("trackRes" , &trackRes); | |
72 | t->SetBranchAddress("pointY" , &pointY); | |
73 | t->SetBranchAddress("pointZ" , &pointZ); | |
74 | t->SetBranchAddress("npTRD" , &npTRD); | |
75 | t->SetBranchAddress("event" , &event); | |
76 | ||
77 | // make binning | |
78 | TVectorD *vR = MakeLinBinning(10,86.,250.); | |
79 | TVectorD *vPhi = MakeLinBinning(18*8,0.,2*TMath::Pi()); | |
80 | TVectorD *vZ = MakeLinBinning(50,-250.,250.); | |
81 | ||
82 | const Int_t nbins=4; | |
83 | Int_t bins[nbins] = {vR->GetNrows()-1, vPhi->GetNrows()-1, vZ->GetNrows()-1, 80}; | |
84 | // Int_t bins[nbins] = {16, 18*5, 50, 80}; | |
85 | Double_t xmin[nbins] = {86. , 0., -250., -2.}; | |
86 | Double_t xmax[nbins] = {250., 2*TMath::Pi(), 250., 2.}; | |
87 | THnF *hn = new THnF("hn", "hn", nbins, bins, xmin, xmax); | |
88 | ||
89 | hn->GetAxis(0)->Set(vR ->GetNrows()-1, vR ->GetMatrixArray()); | |
90 | hn->GetAxis(1)->Set(vPhi->GetNrows()-1, vPhi->GetMatrixArray()); | |
91 | hn->GetAxis(2)->Set(vZ ->GetNrows()-1, vZ ->GetMatrixArray()); | |
92 | ||
93 | ||
94 | for (Int_t iev=0; iev<t->GetEntries(); ++iev) { | |
95 | t->GetEntry(iev); | |
96 | ||
97 | // cuts | |
98 | // -- on trd | |
99 | if (npTRD<2) continue; | |
100 | Double_t pt=1./TMath::Abs(soneOverPt); | |
101 | if (pt<0.8) continue; | |
102 | ||
103 | Float_t resRphiRandom = resRphi*trackRes; | |
104 | Float_t deviation = trackY+resRphiRandom-pointY; | |
105 | ||
106 | Double_t xx[4]={radius, trackPhi, trackZ ,deviation}; | |
107 | hn->Fill(xx); | |
108 | } | |
109 | ||
110 | // do fits and fill tree | |
111 | TTreeSRedirector stream(outFile.Data()); | |
112 | gROOT->cd(); | |
113 | ||
114 | DumpHn(hn, stream); | |
115 | ||
116 | stream.GetFile()->cd(); | |
117 | hn->Write(); | |
118 | ||
119 | delete hn; | |
120 | delete vR; | |
121 | delete vPhi; | |
122 | delete vZ; | |
123 | } | |
124 | ||
125 | ||
126 | void AnaDeltaTree2(TString file, TString outDir=".") | |
127 | { | |
128 | // | |
129 | // NOTE: not finished | |
130 | // | |
131 | TFile f(file); | |
132 | gROOT->cd(); | |
133 | TTree *t = (TTree*)f.Get("delta"); | |
134 | Float_t soneOverPt=0.; | |
135 | Float_t radius=0.; | |
136 | Float_t trackPhi=0.; | |
137 | Float_t trackY=0.; | |
138 | Float_t trackZ=0.; | |
139 | Float_t resRphi=0.; | |
140 | Float_t trackRes=0.; | |
141 | Float_t pointY=0.; | |
142 | Float_t pointZ=0.; | |
143 | Float_t npTRD=0.; | |
144 | Float_t event=0.; | |
145 | ||
146 | t->SetBranchAddress("soneOverPt" , &soneOverPt); | |
147 | t->SetBranchAddress("r" , &radius); | |
148 | t->SetBranchAddress("trackPhi" , &trackPhi); | |
149 | t->SetBranchAddress("trackY" , &trackY); | |
150 | t->SetBranchAddress("trackZ" , &trackZ); | |
151 | t->SetBranchAddress("resRphi" , &resRphi); | |
152 | t->SetBranchAddress("trackRes" , &trackRes); | |
153 | t->SetBranchAddress("pointY" , &pointY); | |
154 | t->SetBranchAddress("pointZ" , &pointZ); | |
155 | t->SetBranchAddress("npTRD" , &npTRD); | |
156 | t->SetBranchAddress("event" , &event); | |
157 | ||
158 | // make binning | |
159 | TVectorD *vZ = MakeLinBinning(50,-250.,250.); | |
160 | TVectorD *vPhi = MakeLinBinning(18*8,0.,TMath::Pi()); | |
161 | TVectorD *vR = MakeLinBinning(16,86.,250.); | |
162 | ||
163 | TObjArray arrZ(vZ->GetNrows()-1); | |
164 | arrZ.SetOwner(); | |
165 | ||
166 | for (Int_t iev=0; iev<t->GetEntries(); ++iev) { | |
167 | t->GetEntry(iev); | |
168 | ||
169 | // cuts | |
170 | // -- on trd | |
171 | if (npTRD<2) continue; | |
172 | ||
173 | Float_t resRphiRandom=resRphi*trackRes; | |
174 | ||
175 | Int_t binZ = TMath::BinarySearch(vZ->GetNrows(),vZ->GetMatrixArray(),(Double_t)trackZ); | |
176 | Int_t binPhi = TMath::BinarySearch(vPhi->GetNrows(),vPhi->GetMatrixArray(),(Double_t)trackPhi); | |
177 | Int_t binR = TMath::BinarySearch(vR->GetNrows(),vR->GetMatrixArray(),(Double_t)radius); | |
178 | ||
179 | if (binZ<0) binZ=0; | |
180 | if (binPhi<0) binPhi=0; | |
181 | if (binR<0) binR=0; | |
182 | ||
183 | TObjArray *arrPhi=(TObjArray*)arrZ.UncheckedAt(binZ); | |
184 | if (!arrPhi) { | |
185 | arrPhi=new TObjArray(vPhi->GetNrows()-1); | |
186 | arrZ.AddAt(arrPhi,binZ); | |
187 | } | |
188 | ||
189 | TObjArray *arrR=(TObjArray*)arrPhi->UncheckedAt(binPhi); | |
190 | if (!arrR) { | |
191 | arrR=new TObjArray(vR->GetNrows()-1); | |
192 | arrPhi->AddAt(arrR,binPhi); | |
193 | } | |
194 | ||
195 | TH1S *h = (TH1S*)arrR->UncheckedAt(binR); | |
196 | if (!h) { | |
197 | h = new TH1S(Form("h_%02d_%02d_%d02",binZ, binPhi, binR), | |
198 | Form("z,phi,r: %02d,%02d,%d02; #Delta r#phi (cm)",binZ, binPhi, binR), | |
199 | 80, -2., 2.); | |
200 | arrR->AddAt(h, binR); | |
201 | } | |
202 | ||
203 | h->Fill(trackY+resRphiRandom-pointY); | |
204 | } | |
205 | ||
206 | // do fits and fill tree | |
207 | } | |
208 | ||
209 | ||
210 | void DumpHn(THn *hn, TTreeSRedirector &stream) | |
211 | { | |
1a18cb2b | 212 | TAxis *ar = hn->GetAxis(0); |
213 | TAxis *aphi = hn->GetAxis(1); | |
214 | TAxis *az = hn->GetAxis(2); | |
03de6136 | 215 | |
1a18cb2b | 216 | |
217 | for (Int_t iz=0; iz<az->GetNbins(); ++iz) { | |
218 | az->SetRange(iz+1,iz+1); | |
219 | TObjArray arrFits; | |
220 | arrFits.SetName(Form("z_%02d",iz)); | |
221 | arrFits.SetOwner(); | |
222 | ||
223 | for (Int_t ir=0; ir<ar->GetNbins(); ++ir) { | |
224 | ar->SetRange(ir+1,ir+1); | |
225 | for (Int_t iphi=0; iphi<aphi->GetNbins(); ++iphi) { | |
226 | aphi->SetRange(iphi+1,iphi+1); | |
03de6136 | 227 | |
1a18cb2b | 228 | TH1 *hProj = hn->Projection(3); |
229 | if (hProj->GetEntries()<1) { | |
230 | delete hProj; | |
231 | continue; | |
232 | } | |
03de6136 | 233 | |
1a18cb2b | 234 | TF1 fg("fg","gaus",-2,2); |
235 | Double_t cr = ar->GetBinCenter(ir+1); | |
236 | Double_t cphi = aphi->GetBinCenter(iphi+1); | |
237 | Double_t cz = az->GetBinCenter(iz+1); | |
238 | hProj->SetNameTitle(Form("h_%02d_%02d_%d02",iz+1, iphi+1, ir+1), | |
03de6136 | 239 | Form("z,phi,r: %02d,%02d,%d02 (%.2f, %.2f, %.2f)",iz+1,iphi+1,ir+1, cr, cphi, cz ) |
240 | ); | |
1a18cb2b | 241 | hProj->Fit(&fg,"QR"); |
242 | arrFits.Add(hProj); | |
243 | ||
244 | Double_t mean = fg.GetParameter(1); | |
245 | Double_t meanErr = fg.GetParError(1); | |
246 | Double_t sigma = fg.GetParameter(2); | |
247 | Double_t sigmaErr = fg.GetParError(2); | |
248 | Int_t entries = hProj->GetEntries(); | |
249 | Double_t chi2ndf = fg.GetChisquare()/fg.GetNDF(); | |
03de6136 | 250 | |
1a18cb2b | 251 | stream << "d" << |
252 | "ir=" << ir << | |
253 | "iphi=" << iphi << | |
254 | "iz=" << iz << | |
255 | "cr=" << cr << | |
256 | "cphi=" << cphi << | |
257 | "cz=" << cz << | |
258 | "mean=" << mean << | |
259 | "meanErr=" << meanErr << | |
260 | "sigma=" << sigma << | |
261 | "sigmaErr=" << sigmaErr << | |
262 | "entries=" << entries << | |
263 | "chi2ndf=" << chi2ndf << | |
264 | "\n"; | |
265 | } | |
266 | } | |
267 | stream.GetFile()->cd(); | |
268 | arrFits.Write(0x0,TObject::kSingleKey); | |
269 | gROOT->cd(); | |
270 | } | |
03de6136 | 271 | |
272 | } | |
1a18cb2b | 273 | |
03de6136 | 274 | //______________________________________________________________________________ |
275 | TVectorD* MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax) | |
276 | { | |
277 | // | |
278 | // Make logarithmic binning | |
279 | // the user has to delete the array afterwards!!! | |
280 | // | |
281 | ||
282 | //check limits | |
283 | if (xmin<1e-20 || xmax<1e-20){ | |
284 | printf("For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!"); | |
285 | return MakeLinBinning(nbinsX, xmin, xmax); | |
286 | } | |
287 | if (xmax<xmin){ | |
288 | Double_t tmp=xmin; | |
289 | xmin=xmax; | |
290 | xmax=tmp; | |
291 | } | |
292 | TVectorD *binLim=new TVectorD(nbinsX+1); | |
293 | Double_t first=xmin; | |
294 | Double_t last=xmax; | |
295 | Double_t expMax=TMath::Log(last/first); | |
296 | for (Int_t i=0; i<nbinsX+1; ++i){ | |
297 | (*binLim)[i]=first*TMath::Exp(expMax/nbinsX*(Double_t)i); | |
298 | } | |
299 | return binLim; | |
300 | } | |
301 | ||
302 | //______________________________________________________________________________ | |
303 | TVectorD* MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax) | |
304 | { | |
305 | // | |
306 | // Make linear binning | |
307 | // the user has to delete the array afterwards!!! | |
308 | // | |
309 | if (xmax<xmin){ | |
310 | Double_t tmp=xmin; | |
311 | xmin=xmax; | |
312 | xmax=tmp; | |
313 | } | |
314 | TVectorD *binLim=new TVectorD(nbinsX+1); | |
315 | Double_t first=xmin; | |
316 | Double_t last=xmax; | |
317 | Double_t binWidth=(last-first)/nbinsX; | |
318 | for (Int_t i=0; i<nbinsX+1; ++i){ | |
319 | (*binLim)[i]=first+binWidth*(Double_t)i; | |
320 | } | |
321 | return binLim; | |
322 | } | |
323 | ||
324 | //_____________________________________________________________________________ | |
325 | TVectorD* MakeArbitraryBinning(const char* bins) | |
326 | { | |
327 | // | |
328 | // Make arbitrary binning, bins separated by a ',' | |
329 | // | |
330 | TString limits(bins); | |
331 | if (limits.IsNull()){ | |
332 | printf("Bin Limit string is empty, cannot add the variable"); | |
333 | return 0x0; | |
334 | } | |
335 | ||
336 | TObjArray *arr=limits.Tokenize(","); | |
337 | Int_t nLimits=arr->GetEntries(); | |
338 | if (nLimits<2){ | |
339 | printf("Need at leas 2 bin limits, cannot add the variable"); | |
340 | delete arr; | |
341 | return 0x0; | |
342 | } | |
343 | ||
344 | TVectorD *binLimits=new TVectorD(nLimits); | |
345 | for (Int_t iLim=0; iLim<nLimits; ++iLim){ | |
346 | (*binLimits)[iLim]=(static_cast<TObjString*>(arr->At(iLim)))->GetString().Atof(); | |
347 | } | |
348 | ||
349 | delete arr; | |
350 | return binLimits; | |
1a18cb2b | 351 | } |
352 | ||
03de6136 | 353 | void PlotFromTree(TTree *d, TString outDir=".") |
354 | { | |
355 | TCanvas *c=new TCanvas; | |
356 | gStyle->SetOptStat(0); | |
357 | d->SetMarkerStyle(20); | |
358 | d->SetMarkerSize(1); | |
359 | ||
360 | TProfile2D pRZ("pRZ",";z (cm); r(cm)",50,-250,250,10,85,250); | |
361 | d->Draw("entries:cr:cz>>pRZ","","profcolz"); | |
362 | pRZ.GetZaxis()->UnZoom(); | |
363 | c->SaveAs(Form("%s/entries_average.png",outDir.Data())); | |
364 | d->Draw("entries:cr:cz>>pRZ","iphi==2","profcolz"); | |
365 | c->SaveAs(Form("%s/entries_onePhi.png",outDir.Data())); | |
366 | ||
367 | pRZ.SetMaximum(0.04); | |
368 | d->Draw("meanErr:cr:cz>>pRZ","","profcolz"); | |
369 | c->SaveAs(Form("%s/meanErr_average.png",outDir.Data())); | |
370 | d->Draw("meanErr:cr:cz>>pRZ","iphi==2","profcolz"); | |
371 | c->SaveAs(Form("%s/meanErr_onePhi.png",outDir.Data())); | |
372 | ||
373 | ||
374 | d->Draw("mean:cphi:cr","iz==25","colz"); | |
375 | c->SaveAs(Form("%s/mean_oneZ_phi_allR.png",outDir.Data())); | |
376 | d->Draw("mean:meanErr:cphi","iz==25&&ir==2","goff"); | |
377 | TGraphErrors grmean_phi(d->GetSelectedRows(),d->GetV3(),d->GetV1(),0,d->GetV2()); | |
378 | grmean_phi->SetTitle(";#varphi;#LT#Delta r#varphi#GT"); | |
379 | grmean_phi->SetMarkerStyle(20); | |
380 | grmean_phi->SetMarkerSize(1); | |
381 | grmean_phi->Draw("ap"); | |
382 | c->SaveAs(Form("%s/mean_oneZ_phi_oneR.png",outDir.Data())); | |
383 | ||
384 | d->Draw("mean:cr:cphi","iz==25","colz"); | |
385 | c->SaveAs(Form("%s/mean_oneZ_r_allPhi.png",outDir.Data())); | |
386 | ||
387 | d->Draw("mean:meanErr:cr","iz==25&&iphi==2","goff"); | |
388 | TGraphErrors grmean_r(d->GetSelectedRows(),d->GetV3(),d->GetV1(),0,d->GetV2()); | |
389 | grmean_r->SetTitle(";r (cm);#LT#Delta r#varphi#GT"); | |
390 | grmean_r->SetMarkerStyle(20); | |
391 | grmean_r->SetMarkerSize(1); | |
392 | grmean_r->Draw("ap"); | |
393 | c->SaveAs(Form("%s/mean_oneZ_r_onePhi.png",outDir.Data())); | |
394 | ||
395 | ||
396 | d->Draw("meanErr:cphi:cr","iz==25","colz"); | |
397 | c->SaveAs(Form("%s/meanErr_oneZ_phi_allR.png",outDir.Data())); | |
398 | d->Draw("meanErr:cphi","iz==25&&ir==2"); | |
399 | c->SaveAs(Form("%s/meanErr_oneZ_phi_oneR.png",outDir.Data())); | |
400 | ||
401 | d->Draw("meanErr:cr:cphi","iz==25","colz"); | |
402 | c->SaveAs(Form("%s/meanErr_oneZ_r_allPhi.png",outDir.Data())); | |
403 | ||
404 | d->Draw("meanErr:cr","iz==25&&iphi==2"); | |
405 | c->SaveAs(Form("%s/meanErr_oneZ_r_onePhi.png",outDir.Data())); | |
406 | ||
407 | } |