6f71a78e1fd46687fc11f9a8fa0a2bb114ec61ab
[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   
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
133 void AnaDeltaTree2(TString file/*, TString outDir="."*/)
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
216 void 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);
231   arrMaps->Add(corrAverage); // correction with the average Map
232   arrMaps->Add(corrFluct);   // distortion with the fluctuation Map
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
237   // !!!!! In AliTPCComposedCorrection::GetDistortion MakeInverseIterator is called !!!!
238   // for this reason we have to add the maps in the wrong order
239   
240   AliTPCComposedCorrection *residualDistortion = new AliTPCComposedCorrection(arrMaps, AliTPCComposedCorrection::kQueueResidual);
241   Float_t dummy=0;
242   TVectorD weights(2);
243   weights(1)=+1.;
244   weights(0)=-AliToyMCEventGenerator::GetSCScalingFactor(corrFluct, corrAverage,dummy);
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 }
299
300 void DumpHn(THn *hn, TTreeSRedirector &stream)
301 {
302   TAxis *ar   = hn->GetAxis(0);
303   TAxis *aphi = hn->GetAxis(1);
304   TAxis *az   = hn->GetAxis(2);
305   
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);
317         
318         TH1 *hProj = hn->Projection(3);
319         if (hProj->GetEntries()<1) {
320           delete hProj;
321           continue;
322         }
323         
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);
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 )
330         );
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();
340
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   }
361   
362 }
363
364 //______________________________________________________________________________
365 TVectorD* 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 //______________________________________________________________________________
393 TVectorD* 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 //_____________________________________________________________________________
415 TVectorD* 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;
441 }
442
443 void 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");
467   TGraphErrors *grmean_phi=new TGraphErrors(d->GetSelectedRows(),d->GetV3(),d->GetV1(),0,d->GetV2());
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");
478   TGraphErrors *grmean_r=new TGraphErrors(d->GetSelectedRows(),d->GetV3(),d->GetV1(),0,d->GetV2());
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   
497 }
498
499