o first try to create a LUT from the residual distortions
[u/mrichter/AliRoot.git] / TPC / Upgrade / macros / AnaDelta.C
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>
11 #include <TVectorT.h>
12 #include <TCanvas.h>
13 #include <TProfile2D.h>
14 #include <TGraphErrors.h>
15 #include <TTreeStream.h>
16
17 #include <AliExternalTrackParam.h>
18 #include <AliTPCComposedCorrection.h>
19 #include <AliTPCCorrectionLookupTable.h>
20
21 #include <AliToyMCEventGenerator.h>
22
23 /*
24
25 .L $ALICE_ROOT/TPC/Upgrade/macros/AnaDelta.C+g
26
27
28 */
29 TVectorD* MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax);
30 TVectorD* MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax);
31 TVectorD* MakeArbitraryBinning(const char* bins);
32
33 void DumpHn(THn *hn, TTreeSRedirector &stream);
34
35 void 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
49   DumpHn(hn, stream);
50
51   delete hn;
52 }
53
54
55 void 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   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)");
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
137 void AnaDeltaTree2(TString file/*, TString outDir="."*/)
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
220 void 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);
235   arrMaps->Add(corrAverage); // correction with the average Map
236   arrMaps->Add(corrFluct);   // distortion with the fluctuation Map
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
241   // !!!!! In AliTPCComposedCorrection::GetDistortion MakeInverseIterator is called !!!!
242   // for this reason we have to add the maps in the wrong order
243   
244   AliTPCComposedCorrection *residualDistortion = new AliTPCComposedCorrection(arrMaps, AliTPCComposedCorrection::kQueueResidual);
245   Float_t dummy=0;
246   TVectorD weights(2);
247   weights(0)=+1.;
248   weights(1)=-AliToyMCEventGenerator::GetSCScalingFactor(corrFluct, corrAverage,dummy);
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 }
303
304 void DumpHn(THn *hn, TTreeSRedirector &stream)
305 {
306   TAxis *ar   = hn->GetAxis(0);
307   TAxis *aphi = hn->GetAxis(1);
308   TAxis *az   = hn->GetAxis(2);
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));
320   
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);
332         
333         TH1 *hProj = hn->Projection(3);
334         if (hProj->GetEntries()<1) {
335           delete hProj;
336           continue;
337         }
338         
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);
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 )
345         );
346         hProj->Fit(&fg,"LMQR");
347         arrFits.Add(hProj);
348
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);
353         Int_t    entries  = hProj->GetEntries();
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         
360         stream << "d" <<
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  <<
377         "\n";
378
379         Double_t x[nbins]={cr, cphi, cz};
380         if (meanErr<0.3) hnRes.Fill(x,mean);
381       }
382     }
383     stream.GetFile()->cd();
384     arrFits.Write(0x0,TObject::kSingleKey);
385     gROOT->cd();
386   }
387
388   stream.GetFile()->cd();
389   hnRes.Write();
390   gROOT->cd();
391 }
392
393 //______________________________________________________________________________
394 TVectorD* 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 //______________________________________________________________________________
422 TVectorD* 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 //_____________________________________________________________________________
444 TVectorD* 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;
470 }
471
472 void 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");
496   TGraphErrors *grmean_phi=new TGraphErrors(d->GetSelectedRows(),d->GetV3(),d->GetV1(),0,d->GetV2());
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");
507   TGraphErrors *grmean_r=new TGraphErrors(d->GetSelectedRows(),d->GetV3(),d->GetV1(),0,d->GetV2());
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   
526 }
527
528