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