d8c0702d9fe230891bc33ac73e80adfbc255adcf
[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 /*
18
19 .L $ALICE_ROOT/TPC/Upgrade/macros/AnaDelta.C+g
20
21
22 */
23 TVectorD* MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax);
24 TVectorD* MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax);
25 TVectorD* MakeArbitraryBinning(const char* bins);
26
27 void DumpHn(THn *hn, TTreeSRedirector &stream);
28
29 void 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
43   DumpHn(hn, stream);
44
45   delete hn;
46 }
47
48
49 void 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
127 void AnaDeltaTree2(TString file/*, TString outDir="."*/)
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
211 void DumpHn(THn *hn, TTreeSRedirector &stream)
212 {
213   TAxis *ar   = hn->GetAxis(0);
214   TAxis *aphi = hn->GetAxis(1);
215   TAxis *az   = hn->GetAxis(2);
216   
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);
228         
229         TH1 *hProj = hn->Projection(3);
230         if (hProj->GetEntries()<1) {
231           delete hProj;
232           continue;
233         }
234         
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),
240                             Form("z,phi,r: %02d,%02d,%d02 (%.2f, %.2f, %.2f)",iz+1,iphi+1,ir+1, cr, cphi, cz )
241         );
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();
251
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   }
272   
273 }
274
275 //______________________________________________________________________________
276 TVectorD* 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 //______________________________________________________________________________
304 TVectorD* 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 //_____________________________________________________________________________
326 TVectorD* 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;
352 }
353
354 void 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");
378   TGraphErrors *grmean_phi=new TGraphErrors(d->GetSelectedRows(),d->GetV3(),d->GetV1(),0,d->GetV2());
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");
389   TGraphErrors *grmean_r=new TGraphErrors(d->GetSelectedRows(),d->GetV3(),d->GetV1(),0,d->GetV2());
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 }