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