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