////////////////////////////////////////////////////////////////////////////
#include "AliTRDclusterResolution.h"
-#include "AliTRDtrackInfo/AliTRDclusterInfo.h"
+#include "info/AliTRDclusterInfo.h"
#include "AliTRDgeometry.h"
#include "AliTRDcalibDB.h"
#include "AliTRDCommonParam.h"
#include "TGraphErrors.h"
#include "TLine.h"
#include "TH2I.h"
+#include "THnSparse.h"
#include "TMath.h"
#include "TLinearFitter.h"
const Float_t AliTRDclusterResolution::fgkTimeBinLength = 1./ AliTRDCommonParam::Instance()->GetSamplingFrequency();
//_______________________________________________________
-AliTRDclusterResolution::AliTRDclusterResolution(const char *name)
- : AliTRDrecoTask(name, "Cluster Error Parametrization")
+AliTRDclusterResolution::AliTRDclusterResolution(const char *name, const char *title)
+ : AliTRDrecoTask(name, title)
,fCanvas(0x0)
,fInfo(0x0)
,fResults(0x0)
TList *l = 0x0;
TObjArray *arr = 0x0;
- TH2 *h2 = 0x0;
+ TH2 *h2 = 0x0;TH1 *h1 = 0x0;
TGraphErrors *gm(0x0), *gs(0x0), *gp(0x0);
switch(ifig){
case kQRes:
case kSigm:
if(!(arr = (TObjArray*)fResults->At(kSigm))) break;
gPad->Divide(2, 1); l = gPad->GetListOfPrimitives();
- for(Int_t ipad = 2; ipad--;){
- if(!(h2 = (TH2F*)arr->At(ipad))) return kFALSE;
- ((TVirtualPad*)l->At(ipad))->cd();
- h2->Draw("lego2fb");
- }
+ if(!(h2 = (TH2F*)arr->At(0))) return kFALSE;
+ ((TVirtualPad*)l->At(0))->cd();
+ h1 = h2->ProjectionY("hsx_pyy"); h1->Scale(1.e4/kND); h1->SetMarkerStyle(24);
+ h1->SetYTitle("<#sigma_{x}> [#mum]");
+ h1->GetXaxis()->SetRange(2, kNTB-1); h1->Draw("pc");
+
+ if(!(h2 = (TH2F*)arr->At(1))) return kFALSE;
+ ((TVirtualPad*)l->At(1))->cd();
+ h1 = h2->ProjectionY("hsy_pyy"); h1->Scale(1.e4/kND); h1->SetMarkerStyle(24);
+ h1->SetYTitle("<#sigma_{y}> [#mum]");
+ h1->GetXaxis()->SetRange(2, kNTB-1); h1->Draw("pc");
return kTRUE;
case kMean:
if(!(arr = (TObjArray*)fResults->At(kMean))) break;
gPad->Divide(2, 1); l = gPad->GetListOfPrimitives();
- for(Int_t ipad = 2; ipad--;){
- if(!(h2 = (TH2F*)arr->At(ipad))) return kFALSE;
- ((TVirtualPad*)l->At(ipad))->cd();
- h2->Draw("lego2fb");
- }
+ ((TVirtualPad*)l->At(0))->cd();
+ if(!(gm = (TGraphErrors*)arr->At(0))) return kFALSE;
+ gm->Draw("apl");
+ gm->GetHistogram()->SetXTitle("t_{drift} [#mus]");
+ gm->GetHistogram()->SetYTitle("dx [#mum]");
+
+ ((TVirtualPad*)l->At(1))->cd();
+ if(!(gm = (TGraphErrors*)arr->At(1))) return kFALSE;
+ gm->Draw("apl");
+ gm->GetHistogram()->SetXTitle("t_{drift} [#mus]");
+ gm->GetHistogram()->SetYTitle("dy [#mum]");
+
return kTRUE;
default:
break;
TH2I *h2 = 0x0;
TObjArray *arr = 0x0;
- fContainer->AddAt(h2 = new TH2I("h_q", "", 50, 2.2, 7.5, 100, -.5, .5), kQRes);
+ fContainer->AddAt(h2 = new TH2I("h_q", "dy=f(q)", 50, 2.2, 7.5, 100, -.5, .5), kQRes);
h2->SetXTitle("log(q) [a.u.]");
h2->SetYTitle("#Delta y[cm]");
h2->SetZTitle("entries");
- fContainer->AddAt(arr = new TObjArray(AliTRDgeometry::kNlayer), kCenter);
- arr->SetName("y(PadWidth)");
- for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
- arr->AddAt(h2 = new TH2I(Form("h_y%d", ily), Form("Ly[%d]", ily), 51, -.51, .51, 100, -.5, .5), ily);
- h2->SetXTitle("y_{w} [w]");
- h2->SetYTitle("#Delta y[cm]");
- h2->SetZTitle("entries");
+ THnSparseS *hn = 0x0;
+ const Int_t nd = 4;
+ // ly x pw dy
+ Int_t nbins[nd] = { 6, 24, 51, 100};
+ Double_t xmin[nd] = {-0.5, -0.5,-0.51,-0.5},
+ xmax[nd] = { 5.5, 23.5, 0.51, 0.5};
+ if(!(hn = (THnSparseS*)gROOT->FindObject("hn"))){
+ hn = new THnSparseS("hn", "dy=f(pw|x,ly)", nd, nbins, xmin, xmax);
+ arr = hn->GetListOfAxes();
+ ((TAxis*)arr->At(0))->SetTitle("layer");
+ ((TAxis*)arr->At(1))->SetTitle("x");
+ ((TAxis*)arr->At(2))->SetTitle("pw");
+ ((TAxis*)arr->At(3))->SetTitle("dy");
}
+ fContainer->AddAt(hn, kCenter);
fContainer->AddAt(arr = new TObjArray(kN), kSigm);
arr->SetName("Resolution");
Int_t det, t;
Float_t x, y, z, q, dy, dydx, dzdx, cov[3], covcl[3];
TH2I *h2 = 0x0;
+ THnSparseS *hn = 0x0;
// define limits around ExB for which x contribution is negligible
const Float_t kDtgPhi = 3.5e-2; //(+- 2 deg)
- TObjArray *arr0 = (TObjArray*)fContainer->At(kCenter);
+ //TObjArray *arr0 = (TObjArray*)fContainer->At(kCenter);
TObjArray *arr1 = (TObjArray*)fContainer->At(kSigm);
TObjArray *arr2 = (TObjArray*)fContainer->At(kMean);
// resolution as a function of y displacement from pad center
// only for phi equal exB
if(TMath::Abs(dydx-fExB) < kDtgPhi){
- h2 = (TH2I*)arr0->At(AliTRDgeometry::GetLayer(det));
- h2->Fill(cli->GetYDisplacement(), dy);
+ hn = (THnSparseS*)fContainer->At(kCenter);
+ Double_t x[]={
+ AliTRDgeometry::GetLayer(det),
+ t,
+ cli->GetYDisplacement(),
+ dy
+ };
+ hn->Fill(x);
}
Int_t it = fAt->FindBin((t+.5)*fgkTimeBinLength);
fAt->GetNbins(), fAt->GetXmin(), fAt->GetXmax());
}
arr->AddAt(h2, 0);
- h2->SetXTitle("d [mm]");
+ h2->SetXTitle("d [cm]");
h2->SetYTitle("t_{drift} [#mus]");
h2->SetZTitle("#sigma_{x} [cm]");
arr->AddAt(h2 = (TH2F*)h2->Clone("hSY"), 1);
h2->SetZTitle("#sigma_{y} [cm]");
- fResults->AddAt(arr = new TObjArray(2), kMean);
+ fResults->AddAt(arr = new TObjArray(4), kMean);
arr->SetOwner();
- arr->AddAt(h2 = (TH2F*)h2->Clone("hDX"), 0);
+ arr->AddAt(g = new TGraphErrors(), 0);
+ g->SetLineColor(kBlue); g->SetMarkerColor(kBlue);
+ g->SetMarkerStyle(24);
+ arr->AddAt(g = new TGraphErrors(), 1);
+ g->SetLineColor(kRed); g->SetMarkerColor(kRed);
+ g->SetMarkerStyle(24);
+ arr->AddAt(h2 = (TH2F*)h2->Clone("hDX"), 2);
h2->SetZTitle("dx [cm]");
- arr->AddAt(h2 = (TH2F*)h2->Clone("hX0"), 1);
- h2->SetZTitle("x0 [cm]");
+ arr->AddAt(h2 = (TH2F*)h2->Clone("hDY"), 3);
+ h2->SetZTitle("dy [cm]");
} else {
TObject *o = 0x0;
TIterator *iter=fResults->MakeIterator();
//_______________________________________________________
void AliTRDclusterResolution::ProcessCenterPad()
{
- TObjArray *arr = (TObjArray*)fContainer->At(kCenter);
- if(!arr) {
- AliWarning("Missing dy=f(y_d) container");
+ THnSparseS *hn = (THnSparseS*)fContainer->At(kCenter);
+ if(!hn) {
+ AliWarning("Missing dy=f(y | x, ly) container");
return;
}
TF1 f("f", "gaus", -.5, .5);
- TAxis *ax = 0x0;
+ TAxis *aly = 0x0, *ax = 0x0, *ay = 0x0;
TH1D *h1 = 0x0, *h11 = 0x0;
- TH2I *h2 = 0x0;
TObjArray *arrg = (TObjArray*)fResults->At(kCenter);
TH2F *hym = (TH2F*)arrg->At(0);
TH2F *hys = (TH2F*)arrg->At(1);
TH2F *hyp = (TH2F*)arrg->At(2);
- for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
- if(!(h2 = (TH2I*)arr->At(ily))) continue;
- ax = h2->GetXaxis();
- for(Int_t ix=1; ix<=ax->GetNbins(); ix++){
- Float_t yd = ax->GetBinCenter(ix);
- h1 = h2->ProjectionY("py", ix, ix);
- Int_t entries = (Int_t)h1->GetEntries();
- if(entries < 50) continue;
- //Adjust(&f, h1);
- h1->Fit(&f, "QN");
-
- // Fill sy = f(y_w)
- hyp->Fill(ily, yd, entries);
- hym->Fill(ily, yd, f.GetParameter(1));
- //hym->SetPointError(ip, 0., f.GetParError(1));
- hys->Fill(ily, yd, f.GetParameter(2));
- //hys->SetPointError(ip, 0., f.GetParError(2));
- }
+
+ aly = hn->GetAxis(0); // layer axis
+ ax = hn->GetAxis(1); // drift length axis
+ ay = hn->GetAxis(2); // y2center axis
+ for(Int_t ily=0; ily<aly->GetNbins(); ily++){
+ aly->SetRangeUser(ily, ily+1);
+ for(Int_t ix=0; ix<ax->GetNbins(); ix++){
+ ax->SetRangeUser(ix, ix+1);
+ for(Int_t iy=0; iy<ay->GetNbins(); ix++){
+ ay->SetRangeUser(iy, iy+1);
+ // finish navigation in the HnSparse
+
+ Double_t yd = ay->GetBinCenter(iy+1);
+ h1 = hn->Projection(3, "O");
+ Int_t entries = (Int_t)h1->GetEntries();
+ if(entries < 50) continue;
+ //Adjust(&f, h1);
+ h1->Fit(&f, "QN");
+
+ // Fill sy = f(y_w)
+ hyp->Fill(ily, yd, entries);
+ hym->Fill(ily, yd, f.GetParameter(1));
+ //hym->SetPointError(ip, 0., f.GetParError(1));
+ hys->Fill(ily, yd, f.GetParameter(2));
+ //hys->SetPointError(ip, 0., f.GetParError(2));
+ }
+ }
}
// POSTPROCESS SPECTRA
TAxis *ax = 0x0;
TH1D *h1 = 0x0;
TH2I *h2 = 0x0;
+ TH1 *hFrame=0x0;
+
// init visualization
TGraphErrors *ggs = 0x0;
- TLine *line = 0x0;
+ TGraph *line = 0x0;
if(fCanvas){
ggs = new TGraphErrors();
- line = new TLine();
+ line = new TGraph();
+ line->SetLineColor(kRed);line->SetLineWidth(2);
}
- Double_t d(0.), x(0.), sx, sy;
+ Double_t d(0.), x(0.), sx, sy, exb2=0.;//fExB*fExB;
TObjArray *arrr = (TObjArray*)fResults->At(kSigm);
TH2F *hsx = (TH2F*)arrr->At(0);
TH2F *hsy = (TH2F*)arrr->At(1);
ax = h2->GetXaxis();
for(Int_t ix=1; ix<=ax->GetNbins(); ix++){
Float_t dydx = ax->GetBinCenter(ix);
+ //if(TMath::Abs(dydx)>0.18) continue;
Double_t tgg = (dydx-fExB)/(1.+dydx*fExB);
Double_t tgg2 = tgg*tgg;
h1 = h2->ProjectionY("py", ix, ix);
Int_t ip = ggs->GetN();
ggs->SetPoint(ip, tgg2, s2);
ggs->SetPointError(ip, 0., s2e);
- //printf("%f, %f, %f, \n", tgg2, s2, s2e);
}
if(gs.Eval()) continue;
- sx = gs.GetParameter(1);
- sy = gs.GetParameter(0);
+
+ // s^2_x = s0^2_x - x^2*tg^2(a_L)/12
+ sx = gs.GetParameter(1)/* - x*x*exb2/12.*/;
+ if(sx<0.) continue;
hsx->SetBinContent(id, it, TMath::Sqrt(sx));
- //hsx->SetBinError(id, it, sig.GetParError(1));
+ //hsx->SetBinError(id, it, .5*gs.GetParError(1)/TMath::Sqrt(sx));
- // s^2_y = s0^2_y + tg^2(a_L) * s^2_x
- hsy->SetBinContent(id, it, TMath::Sqrt(sy)/* - exb2*sig.GetParameter(1)*/);
+ // s^2_y = s0^2_y + tg^2(a_L) * s^2_x
+ // s0^2_y = f(D_L)*x + s_PRF^2
+ sy= gs.GetParameter(0)/*-exb2*sx*/;
+ if(sy <0.) continue;
+ hsy->SetBinContent(id, it, TMath::Sqrt(sy));
//hsy->SetBinError(id, it, sig.GetParError(0)+exb2*exb2*sig.GetParError(1));
- /*if(!fCanvas) */continue;
- fCanvas->cd();
- new(line) TLine(0, sy, .025, sy+.025*sx);
- line->SetLineColor(kRed);line->SetLineWidth(2);
- ggs->Draw("apl"); line->Draw("");
+ if(!fCanvas) continue;
+ fCanvas->cd(); fCanvas->SetLogx(); //fCanvas->SetLogy();
+ if(!hFrame){
+ hFrame=new TH1I("hFrame", "", 100, 0., .3);
+ hFrame->SetMinimum(0.);hFrame->SetMaximum(.005);
+ hFrame->SetXTitle("tg^{2}(#phi-#alpha_{L})");
+ hFrame->SetYTitle("#sigma^{2}y[cm^{2}]");
+ hFrame->SetLineColor(1);hFrame->SetLineWidth(1);
+ hFrame->Draw();
+ } else hFrame->Reset();
+ Double_t xx = 0., dxx=.2/50;
+ for(Int_t ip=0;ip<50;ip++){
+ line->SetPoint(ip, xx, gs.GetParameter(0)+xx*gs.GetParameter(1));
+ xx+=dxx;
+ }
+ ggs->Draw("pl"); line->Draw("l");
fCanvas->Modified(); fCanvas->Update();
if(IsSaveAs()) fCanvas->SaveAs(Form("Figures/ProcessSigma_D%d_T%02d.gif", id, it));
else gSystem->Sleep(100);
- printf(" xd=%4.1f[cm] sx=%5.3e[cm] sy=%5.3e[cm]\n", x, sx, sy);
+ printf(" xd=%4.1f[cm] sx=%5.3e[cm] sy=%5.3e[cm]\n", x, TMath::Sqrt(sx), TMath::Sqrt(sy));
}
}
-// if(ggs) delete ggs; // TODO fix this memory leak !
-// if(line) delete line;
-
- // fit sy parallel to the drift
- TF1 fsyD("fsy", "[0]+[1]*exp(1./(x+[2]))", 0.5, 3.5);
- fsyD.SetParameter(0, .025);
- fsyD.SetParameter(1, -8.e-4);
- fsyD.SetParameter(2, -.25);
- h1 = hsy->ProjectionY("hsy_pyy"); h1->Scale(1./hsy->GetNbinsX());
- for(Int_t ib=1; ib<=h1->GetNbinsX(); ib++) h1->SetBinError(ib, 0.002);
- h1->Fit(&fsyD, "Q", "", .5, 3.5);
- printf(" const Float_t sy0 = %5.3e;\n", fsyD.GetParameter(0));
- printf(" const Float_t sya = %5.3e;\n", fsyD.GetParameter(1));
- printf(" const Float_t syb = %5.3e;\n", fsyD.GetParameter(2));
- if(fCanvas){
- fCanvas->cd();
- h1->Draw("e1"); fsyD.Draw("same");
- fCanvas->Modified(); fCanvas->Update();
- if(IsSaveAs()) fCanvas->SaveAs("Figures/ProcessSigma_SY||.gif");
- else gSystem->Sleep(1000);
- }
-
- // fit sx parallel to the drift
- TF1 fsxD("fsx", "gaus(0)+expo(3)", 0., 3.5);
- h1 = hsx->ProjectionY("hsx_pyy");
- h1->Scale(1./hsx->GetNbinsX());
- for(Int_t ib=1; ib<=h1->GetNbinsX(); ib++) if(h1->GetBinContent(ib)>0.) h1->SetBinError(ib, 0.02);
- h1->Fit(&f, "QN", "", 0., 1.3);
- fsxD.SetParameter(0, 0.);
- fsxD.SetParameter(1, f.GetParameter(1));
- fsxD.SetParameter(2, f.GetParameter(2));
- TF1 expo("fexpo", "expo", 0., 3.5);
- h1->Fit(&expo, "QN");
- fsxD.SetParameter(3, expo.GetParameter(0));
- fsxD.SetParameter(4, expo.GetParameter(1));
- h1->Fit(&fsxD, "Q");
- printf(" const Float_t sxgc = %5.3e;\n", fsxD.GetParameter(0));
- printf(" const Float_t sxgm = %5.3e;\n", fsxD.GetParameter(1));
- printf(" const Float_t sxgs = %5.3e;\n", fsxD.GetParameter(2));
- printf(" const Float_t sxe0 = %5.3e;\n", fsxD.GetParameter(3));
- printf(" const Float_t sxe1 = %5.3e;\n", fsxD.GetParameter(4));
- if(fCanvas){
- fCanvas->cd();
- h1->Draw("e1"); fsxD.Draw("same");
- fCanvas->Modified(); fCanvas->Update();
- if(IsSaveAs()) fCanvas->SaveAs("Figures/ProcessSigma_SX||.gif");
- else gSystem->Sleep(1000);
+ printf(" const Double_t sx[%d][%d]={\n", kNTB-1, kND);
+ for(Int_t iy=1; iy<kNTB; iy++){
+ printf(" {");
+ for(Int_t ix=1; ix<kND; ix++){
+ printf("%5.3e, ", hsx->GetBinContent(ix, iy));
+ }
+ printf("%5.3e}", hsx->GetBinContent(kND, iy));
+ printf("%c\n", iy==(kNTB-1)?' ':',');
}
+ printf(" };\n");
- // fit sx perpendicular to the drift
- Int_t nb = 0;
- TAxis *ay = hsx->GetYaxis();
- TH1D *hx = hsx->ProjectionX("hsx_px", 1,1);
- TH1D *hAn = (TH1D*)hx->Clone("hAn"); hAn->Reset();
- for(Int_t ib = 11; ib<24; ib++, nb++){
- hx = hsx->ProjectionX("hsx_px", ib, ib);
- for(Int_t id=1; id<=hx->GetNbinsX(); id++)
- hx->SetBinContent(id, hx->GetBinContent(id)-fsxD.Eval(ay->GetBinCenter(ib)));
- hAn->Add(hx);
- }
- TF1 fsxA("fsxA", "pol2", 0., 0.25);
- hAn->Scale(1./nb);
- for(Int_t ib=1; ib<=hAn->GetNbinsX(); ib++) hAn->SetBinError(ib, 0.002);
- hAn->Fit(&fsxA, "Q");
- printf(" const Float_t sxd0 = %5.3e;\n", fsxA.GetParameter(0));
- printf(" const Float_t sxd1 = %5.3e;\n", fsxA.GetParameter(1));
- printf(" const Float_t sxd2 = %5.3e;\n", fsxA.GetParameter(2));
- if(fCanvas){
- fCanvas->cd();
- hAn->Draw("e1"); fsxA.Draw("same");
- fCanvas->Modified(); fCanvas->Update();
- if(IsSaveAs()) fCanvas->SaveAs("Figures/ProcessSigma_SXL.gif");
- else gSystem->Sleep(100);
+ printf(" const Double_t sy[%d][%d]={\n", kNTB-1, kND);
+ for(Int_t iy=1; iy<kNTB; iy++){
+ printf(" {");
+ for(Int_t ix=1; ix<kND; ix++){
+ printf("%5.3e, ", hsy->GetBinContent(ix, iy));
+ }
+ printf("%5.3e}", hsy->GetBinContent(kND, iy));
+ printf("%c\n", iy==(kNTB-1)?' ':',');
}
+ printf(" };\n");
+
+ return;
}
//_______________________________________________________
}
TGraphErrors *gm = new TGraphErrors();
TF1 f("f", "gaus", -.5, .5);
- TF1 line("l", Form("%6.4f*[0]+[1]*x", fExB), -.15, .15);
- Double_t d(0.), x(0.), dx, x0;
+ TF1 line("l", "[0]+[1]*x", -.15, .15);
+ Double_t d(0.), x(0.), dx, dy;
TAxis *ax = 0x0;
TH1D *h1 = 0x0;
TH2I *h2 = 0x0;
+ TH1 *hFrame=0x0;
TObjArray *arrr = (TObjArray*)fResults->At(kMean);
- TH2F *hdx = (TH2F*)arrr->At(0);
- TH2F *hx0 = (TH2F*)arrr->At(1);
+ TH2F *hdx = (TH2F*)arrr->At(2);
+ TH2F *hdy = (TH2F*)arrr->At(3);
for(Int_t id=1; id<=fAd->GetNbins(); id++){
d = fAd->GetBinCenter(id); //[mm]
printf(" Doing d = %5.3f [mm]\n", d);
for(Int_t ix=1; ix<=ax->GetNbins(); ix++){
Double_t dydx = ax->GetBinCenter(ix);
h1 = h2->ProjectionY("py", ix, ix);
- if(h1->GetEntries() < 50) continue;
+ if(h1->GetEntries() < 200) continue;
//Adjust(&f, h1);
h1->Fit(&f, "QN");
if(gm->GetN()<4) continue;
gm->Fit(&line, "QN");
- dx = line.GetParameter(1); // = (fVdrift + dVd)*(fT + tCorr);
- x0 = line.GetParameter(0); // = tg(a_L)*dx
+ dx = line.GetParameter(1);
+ Double_t xs = line.GetParameter(0);
+ dy = xs + fExB*dx; // xs = dy - tg(a_L)*dx
hdx->SetBinContent(id, it, dx);
- hx0->SetBinContent(id, it, x0);
+ hdy->SetBinContent(id, it, dy);
if(!fCanvas) continue;
fCanvas->cd();
- gm->Draw("apl"); line.Draw("same");
+ if(!hFrame){
+ hFrame=new TH1I("hFrame", "", 100, -.3, .3);
+ hFrame->SetMinimum(-.1);hFrame->SetMaximum(.1);
+ hFrame->SetXTitle("tg#phi-htg#theta");
+ hFrame->SetYTitle("#Deltay[cm]");
+ hFrame->SetLineColor(1);hFrame->SetLineWidth(1);
+ hFrame->Draw();
+ } else hFrame->Reset();
+ gm->Draw("pl"); line.Draw("same");
fCanvas->Modified(); fCanvas->Update();
if(IsSaveAs()) fCanvas->SaveAs(Form("Figures/ProcessMean_D%d_T%02d.gif", id, it));
else gSystem->Sleep(100);
- printf(" xd=%4.1f[cm] dx=%5.3e[cm] x0=%5.3e[cm]\n", x, dx, x0);
+ printf(" xd=%4.1f[cm] dx=%5.3e[cm] dy=%5.3e[cm] xs=%5.3e[cm]\n", x, dx, dy, xs);
}
}
- TH1D *h=hdx->ProjectionY("hdx_py"); h->Scale(1./hdx->GetNbinsX());
- printf(" const Float_t cx[] = {\n ");
- for(Int_t ib=1; ib<=h->GetNbinsX(); ib++){
- printf("%6.4e, ", h->GetBinContent(ib));
- if(ib%6==0) printf("\n ");
+ // dump to stdout correction map
+ printf(" const Double_t dx[%d][%d]={\n", kNTB-1, kND);
+ for(Int_t iy=1; iy<kNTB; iy++){
+ printf(" {");
+ for(Int_t ix=1; ix<kND; ix++){
+ printf("%+5.3e,", hdx->GetBinContent(ix, iy));
+ }
+ printf("%+5.3e}", hdx->GetBinContent(kND, iy));
+ printf("%c\n", iy==(kNTB-1)?' ':',');
}
printf(" };\n");
+ // Collapse the z direction
+ TAxis *ay = hdx->GetYaxis();
+ // radial systematics
+ TGraphErrors *g = (TGraphErrors*)arrr->At(0);
+ for(Int_t iy = 1; iy<kNTB; iy++){
+ Double_t m=0., rms=0.;
+ for(Int_t ix = 1; ix<=kND; ix++){
+ d = hdx->GetBinContent(ix, iy);
+ m += d; rms += (d*d);
+ }
+ m /= Int_t(kND); rms = TMath::Sqrt(rms/Int_t(kND)-m*m);
+ g->SetPoint(iy-1, ay->GetBinCenter(iy), 1.e4*m);
+ g->SetPointError(iy-1, 0., 1.e4*rms);
+ }
+
+ // r-phi systematics
+ g = (TGraphErrors*)arrr->At(1);
+ for(Int_t iy = 1; iy<kNTB; iy++){
+ Double_t m=0., rms=0.;
+ for(Int_t ix = 1; ix<=kND; ix++){
+ d = hdy->GetBinContent(ix, iy);
+ m += d; rms += (d*d);
+ }
+ m /= Int_t(kND); rms = TMath::Sqrt(rms/Int_t(kND)-m*m);
+ g->SetPoint(iy-1, ay->GetBinCenter(iy), 1.e4*m);
+ g->SetPointError(iy-1, 0., 1.e4*rms);
+ }
// delete gm; TODO memory leak ?
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff