#include "TObjArray.h"
#include "TAxis.h"
#include "TF1.h"
+#include "TLegend.h"
#include "TGraphErrors.h"
#include "TLine.h"
#include "TH2I.h"
+#include "TH3S.h"
+#include "TTree.h"
#include "TMath.h"
#include "TLinearFitter.h"
#include "TCanvas.h"
#include "TSystem.h"
-
ClassImp(AliTRDclusterResolution)
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)
,fAt(0x0)
- ,fAd(0x0)
,fStatus(0)
,fDet(-1)
,fExB(0.)
,fVdrift(0.)
+ ////////////
+ ,fLy(0)
+ ,fX(0.)
+ ,fY(0.)
+ ,fZ(0.)
{
+ memset(fR, 0, 4*sizeof(Float_t));
+ memset(fP, 0, 4*sizeof(Float_t));
// time drift axis
fAt = new TAxis(kNTB, 0., kNTB*fgkTimeBinLength);
- // z axis spans the drift cell 2.5 mm
- fAd = new TAxis(kND, 0., .25);
// By default register all analysis
// The user can switch them off in his steering macro
AliTRDclusterResolution::~AliTRDclusterResolution()
{
if(fCanvas) delete fCanvas;
- if(fAd) delete fAd;
if(fAt) delete fAt;
if(fResults){
fResults->Delete();
Bool_t AliTRDclusterResolution::GetRefFigure(Int_t ifig)
{
if(!fResults) return kFALSE;
-
+ TLegend *leg = 0x0;
TList *l = 0x0;
TObjArray *arr = 0x0;
TH2 *h2 = 0x0;TH1 *h1 = 0x0;
if(!(gs = (TGraphErrors*)arr->At(1))) break;
if(!(gp = (TGraphErrors*)arr->At(2))) break;
gs->Draw("apl");
- gs->GetHistogram()->GetYaxis()->SetRangeUser(-.01, .6);
+ gs->GetHistogram()->GetYaxis()->SetRangeUser(-50., 700.);
gs->GetHistogram()->SetXTitle("Q [a.u.]");
- gs->GetHistogram()->SetYTitle("#sigma_{y} / #mu_{y} [mm] / freq");
+ gs->GetHistogram()->SetYTitle("#sigma_{y} / #mu_{y} [#mum] / freq");
gm->Draw("pl");
gp->Draw("pl");
return kTRUE;
case kCenter:
if(!(arr = (TObjArray*)fResults->At(kCenter))) break;
- gPad->Divide(3, 1); l = gPad->GetListOfPrimitives();
- for(Int_t ipad = 3; ipad--;){
- if(!(h2 = (TH2F*)arr->At(ipad))) return kFALSE;
- ((TVirtualPad*)l->At(ipad))->cd();
- h2->Draw("lego2fb");
+ gPad->Divide(2, 1); l = gPad->GetListOfPrimitives();
+ ((TVirtualPad*)l->At(0))->cd();
+ ((TTree*)arr->At(0))->Draw("y:x>>h(23, 0.1, 2.4, 51, -.51, .51)",
+ "m[0]*(ly==0&&abs(m[0])<1.e-1)", "colz");
+ ((TVirtualPad*)l->At(1))->cd();
+ leg= new TLegend(.7, .7, .9, .95);
+ leg->SetBorderSize(0); leg->SetFillColor(0); leg->SetFillStyle(0);
+ leg->SetHeader("TRD Plane");
+ for(Int_t il = 1; il<=AliTRDgeometry::kNlayer; il++){
+ if(!(gm = (TGraphErrors*)arr->At(il))) return kFALSE;
+ gm->Draw(il>1?"pc":"apc"); leg->AddEntry(gm, Form("%d", il-1), "pl");
+ if(il>1) continue;
+ gm->GetHistogram()->SetXTitle("t_{drift} [#mus]");
+ gm->GetHistogram()->SetYTitle("#sigma_{y}(x|cen=0) [#mum]");
+ gm->GetHistogram()->GetYaxis()->SetRangeUser(150., 500.);
}
+ leg->Draw();
return kTRUE;
case kSigm:
if(!(arr = (TObjArray*)fResults->At(kSigm))) break;
TObjArray* AliTRDclusterResolution::Histos()
{
if(fContainer) return fContainer;
- fContainer = new TObjArray(sizeof(AliTRDclusterResolution::EResultContainer));
+ fContainer = new TObjArray(kNtasks);
//fContainer->SetOwner(kTRUE);
- TH2I *h2 = 0x0;
+ TH3S *h3 = 0x0;
TObjArray *arr = 0x0;
- fContainer->AddAt(h2 = new TH2I("h_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");
+ fContainer->AddAt(arr = new TObjArray(2*AliTRDgeometry::kNlayer), kCenter);
+ arr->SetName("Center");
+ for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
+ // add resolution plot for each layer
+ if(!(h3=(TH3S*)gROOT->FindObject(Form("hCenResLy%d", il)))){
+ h3 = new TH3S(
+ Form("hCenResLy%d", il),
+ Form(" ly [%d]", il),
+ kNTB, fAt->GetBinLowEdge(1), fAt->GetBinUpEdge(kNTB), // x
+ 51, -.51, .51, // y
+ 60, -.3, .3); // dy
+ h3->SetXTitle("x [#mus]");
+ h3->SetYTitle("y [pw]");
+ h3->SetZTitle("#Delta y[cm]");
+ } h3->Reset();
+ arr->AddAt(h3, il);
+ // add Pull plot for each layer
+ if(!(h3=(TH3S*)gROOT->FindObject(Form("hCenPullLy%d", il)))){
+ h3 = new TH3S(
+ Form("hCenPullLy%d", il),
+ Form(" ly [%d]", il),
+ kNTB, fAt->GetBinLowEdge(1), fAt->GetBinUpEdge(kNTB), // x
+ 51, -.51, .51, // y
+ 60, -4., 4.); // dy
+ h3->SetXTitle("x [#mus]");
+ h3->SetYTitle("y [pw]");
+ h3->SetZTitle("#Delta y/#sigma_{y}");
+ } h3->Reset();
+ arr->AddAt(h3, AliTRDgeometry::kNlayer+il);
}
- fContainer->AddAt(arr = new TObjArray(kN), kSigm);
+ if(!(h3 = (TH3S*)gROOT->FindObject("Charge"))){
+ h3 = new TH3S("Charge", "dy=f(q)", 50, 2.2, 7.5, 60, -.3, .3, 60, -4., 4.);
+ h3->SetXTitle("log(q) [a.u.]");
+ h3->SetYTitle("#Delta y[cm]");
+ h3->SetZTitle("#Delta y/#sigma_{y}");
+ }
+ fContainer->AddAt(h3, kQRes);
+
+ fContainer->AddAt(arr = new TObjArray(kNTB), kSigm);
arr->SetName("Resolution");
- Int_t ih = 0;
- for(Int_t id=1; id<=fAd->GetNbins(); id++){
- for(Int_t it=1; it<=fAt->GetNbins(); it++){
- arr->AddAt(h2 = new TH2I(Form("hr_d%02dt%02d", id, it), Form("d_{wire}(%3.1f-%3.1f)[mm] t_{drift}(%3.1f-%3.1f)[#mus]", 10.*fAd->GetBinLowEdge(id), 10.*fAd->GetBinUpEdge(id), fAt->GetBinLowEdge(it), fAt->GetBinUpEdge(it)), 35, -.35, .35, 100, -.5, .5), ih++);
- h2->SetXTitle("tg#phi");
- h2->SetYTitle("#Delta y[cm]");
- h2->SetZTitle("entries");
+ for(Int_t ix=0; ix<kNTB; ix++){
+ if(!(h3=(TH3S*)gROOT->FindObject(Form("hr_x%02d", ix)))){
+ h3 = new TH3S(
+ Form("hr_x%02d", ix),
+ Form(" t_{drift}(%3.1f-%3.1f)[#mus]", fAt->GetBinLowEdge(ix+1), fAt->GetBinUpEdge(ix+1)),
+ kND, 0., 2.5, // z
+ 35, -.35, .35, // tgp
+ 60, -.3, .3); // dy
+ h3->SetXTitle("z [mm]");
+ h3->SetYTitle("tg#phi");
+ h3->SetZTitle("#Delta y[cm]");
}
+ arr->AddAt(h3, ix);
}
- fContainer->AddAt(arr = new TObjArray(kN), kMean);
+ fContainer->AddAt(arr = new TObjArray(kNTB), kMean);
arr->SetName("Systematics");
- ih = 0;
- for(Int_t id=1; id<=fAd->GetNbins(); id++){
- for(Int_t it=1; it<=fAt->GetNbins(); it++){
- arr->AddAt(h2 = new TH2I(Form("hs_d%02dt%02d", id, it), Form("d_{wire}(%3.1f-%3.1f)[mm] t_{drift}(%3.1f-%3.1f)[#mus]", 10.*fAd->GetBinLowEdge(id), 10.*fAd->GetBinUpEdge(id), fAt->GetBinLowEdge(it), fAt->GetBinUpEdge(it)), 35, -.35, .35, 100, -.5, .5), ih++);
- h2->SetXTitle("tg#phi-h*tg(#theta)");
- h2->SetYTitle("#Delta y[cm]");
- h2->SetZTitle("entries");
+ for(Int_t ix=0; ix<kNTB; ix++){
+ if(!(h3=(TH3S*)gROOT->FindObject(Form("hs_x%02d", ix)))){
+ h3 = new TH3S(
+ Form("hs_x%02d", ix),
+ Form(" t_{drift}(%3.1f-%3.1f)[#mus]", fAt->GetBinLowEdge(ix+1), fAt->GetBinUpEdge(ix+1)),
+ kND, 0., 2.5, // z
+ 35, -.35, .35, // tgp-h tgt
+ 60, -.3, .3); // dy
+ h3->SetXTitle("z [mm]");
+ h3->SetYTitle("tg(#phi) - h*tg(#theta)");
+ h3->SetZTitle("#Delta y[cm]");
}
+ arr->AddAt(h3, ix);
}
return fContainer;
Int_t det, t;
Float_t x, y, z, q, dy, dydx, dzdx, cov[3], covcl[3];
- TH2I *h2 = 0x0;
+ TH3S *h3 = 0x0;
// define limits around ExB for which x contribution is negligible
const Float_t kDtgPhi = 3.5e-2; //(+- 2 deg)
// resolution as a function of cluster charge
// only for phi equal exB
if(TMath::Abs(dydx-fExB) < kDtgPhi){
- h2 = (TH2I*)fContainer->At(kQRes);
- h2->Fill(TMath::Log(q), dy);
+ h3 = (TH3S*)fContainer->At(kQRes);
+ h3->Fill(TMath::Log(q), dy, dy/TMath::Sqrt(covcl[0]));
}
// do not use problematic clusters in resolution analysis
// TODO define limits as calibration aware (gain) !!
if(q<20. || q>250.) continue;
+ x = (t+.5)*fgkTimeBinLength; // conservative approach !!
+
// 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);
+ if(TMath::Abs(dydx-fExB) < kDtgPhi/* &&
+ TMath::Abs(x-0.675)<0.225*/){
+ Int_t ly(AliTRDgeometry::GetLayer(det));
+ h3 = (TH3S*)arr0->At(ly);
+ h3->Fill(x, cli->GetYDisplacement(), dy);
+ h3 = (TH3S*)arr0->At(AliTRDgeometry::kNlayer+ly);
+ h3->Fill(x, cli->GetYDisplacement(), dy/TMath::Sqrt(covcl[0]));
}
- Int_t it = fAt->FindBin((t+.5)*fgkTimeBinLength);
- if(it==0 || it == fAt->GetNbins()+1){
- AliWarning(Form("Drift time %f outside allowed range", t));
- continue;
- }
- Int_t id = fAd->FindBin(cli->GetAnisochronity());
- if(id==0 || id == fAd->GetNbins()+1){
- AliWarning(Form("Distance to anode %f outside allowed range", cli->GetAnisochronity()));
+ Int_t ix = fAt->FindBin(x);
+ if(ix==0 || ix == fAt->GetNbins()+1){
+ AliWarning(Form("Drift time %3.1f outside allowed range", x));
continue;
}
- // calculate index of cluster position in array
- Int_t hid = (id-1)*kNTB+it-1;
// fill histo for resolution (sigma)
- h2 = (TH2I*)arr1->At(hid);
- h2->Fill(dydx, dy);
+ ((TH3S*)arr1->At(ix-1))->Fill(10.*cli->GetAnisochronity(), dydx, dy);
// fill histo for systematic (mean)
- h2 = (TH2I*)arr2->At(hid);
- h2->Fill(dydx-cli->GetTilt()*dzdx, dy);
+ ((TH3S*)arr2->At(ix-1))->Fill(10.*cli->GetAnisochronity(), dydx-cli->GetTilt()*dzdx, dy);
}
PostData(0, fContainer);
}
if(!HasExB()) AliWarning("ExB was not set. Call SetExB() before running the post processing.");
TObjArray *arr = 0x0;
- TH2 *h2 = 0x0;
+ TTree *t=0x0;
if(!fResults){
TGraphErrors *g = 0x0;
- fResults = new TObjArray(sizeof(AliTRDclusterResolution::EResultContainer));
+ fResults = new TObjArray(kNtasks);
fResults->SetOwner();
fResults->AddAt(arr = new TObjArray(3), kQRes);
arr->SetOwner();
g->SetLineColor(kGreen); g->SetMarkerColor(kGreen);
g->SetMarkerStyle(7);
- fResults->AddAt(arr = new TObjArray(3), kCenter);
+ // pad center dependence
+ fResults->AddAt(arr = new TObjArray(AliTRDgeometry::kNlayer+1), kCenter);
arr->SetOwner();
-
- if(!(h2 = (TH2F*)gROOT->FindObject("hYM"))){
- h2 = new TH2F("hYM", "",
- AliTRDgeometry::kNlayer, -.5, AliTRDgeometry::kNlayer-.5, 51, -.51, .51);
- }
- arr->AddAt(h2, 0);
- h2->SetXTitle("ly");
- h2->SetYTitle("y [w]");
- h2->SetZTitle("#mu_{x} [cm]");
- arr->AddAt(h2 = (TH2F*)h2->Clone("hYS"), 1);
- h2->SetZTitle("#sigma_{x} [cm]");
- arr->AddAt(h2 = (TH2F*)h2->Clone("hYP"), 2);
- h2->SetZTitle("entries");
-
- fResults->AddAt(arr = new TObjArray(2), kSigm);
- arr->SetOwner();
- if(!(h2 = (TH2F*)gROOT->FindObject("hSX"))){
- h2 = new TH2F("hSX", "",
- fAd->GetNbins(), fAd->GetXmin(), fAd->GetXmax(),
- fAt->GetNbins(), fAt->GetXmin(), fAt->GetXmax());
+ arr->AddAt(
+ t = new TTree("cent", "dy=f(y,x,ly)"), 0);
+ t->Branch("ly", &fLy, "ly/B");
+ t->Branch("x", &fX, "x/F");
+ t->Branch("y", &fY, "y/F");
+ t->Branch("m", &fR[0], "m[2]/F");
+ t->Branch("s", &fR[2], "s[2]/F");
+ t->Branch("pm", &fP[0], "pm[2]/F");
+ t->Branch("ps", &fP[2], "ps[2]/F");
+ for(Int_t il=1; il<=AliTRDgeometry::kNlayer; il++){
+ arr->AddAt(g = new TGraphErrors(), il);
+ g->SetLineColor(il); g->SetLineStyle(il);
+ g->SetMarkerColor(il);g->SetMarkerStyle(4);
}
- arr->AddAt(h2, 0);
- 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(4), kMean);
- arr->SetOwner();
- 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("hDY"), 3);
- h2->SetZTitle("dy [cm]");
+
+
+ fResults->AddAt(t = new TTree("sigm", "dy=f(dw,x,dydx)"), kSigm);
+ t->Branch("x", &fX, "x/F");
+ t->Branch("z", &fZ, "z/F");
+ t->Branch("sx", &fR[0], "sx[2]/F");
+ t->Branch("sy", &fR[2], "sy[2]/F");
+
+
+ fResults->AddAt(t = new TTree("mean", "dy=f(dw,x,dydx - h dzdx)"), kMean);
+ t->Branch("x", &fX, "x/F");
+ t->Branch("z", &fZ, "z/F");
+ t->Branch("dx", &fR[0], "dx[2]/F");
+ t->Branch("dy", &fR[2], "dy[2]/F");
} else {
TObject *o = 0x0;
TIterator *iter=fResults->MakeIterator();
void AliTRDclusterResolution::ProcessCharge()
{
TH2I *h2 = 0x0;
- if((h2 = (TH2I*)fContainer->At(kQRes))) {
+ if(!(h2 = (TH2I*)fContainer->At(kQRes))) {
AliWarning("Missing dy=f(Q) histo");
return;
}
// Fill sy^2 = f(q)
Int_t ip = gqm->GetN();
- gqm->SetPoint(ip, q, 10.*f.GetParameter(1));
- gqm->SetPointError(ip, 0., 10.*f.GetParError(1));
+ gqm->SetPoint(ip, q, 1.e4*f.GetParameter(1));
+ gqm->SetPointError(ip, 0., 1.e4*f.GetParError(1));
// correct sigma for ExB effect
- gqs->SetPoint(ip, q, 1.e1*f.GetParameter(2)/**f.GetParameter(2)-exb2*sxd2*/);
- gqs->SetPointError(ip, 0., 1.e1*f.GetParError(2)/**f.GetParameter(2)*/);
+ gqs->SetPoint(ip, q, 1.e4*f.GetParameter(2)/**f.GetParameter(2)-exb2*sxd2*/);
+ gqs->SetPointError(ip, 0., 1.e4*f.GetParError(2)/**f.GetParameter(2)*/);
// save probability
n += entries;
for(Int_t ip=gqp->GetN(); ip--;){
gqp->GetPoint(ip, q, entries);
entries/=n;
- gqp->SetPoint(ip, q, entries);
+ gqp->SetPoint(ip, q, 1.e3*entries);
gqs->GetPoint(ip, q, sy);
sm += entries*sy;
}
// error parametrization s(q) = <sy> + b(1/q-1/q0)
TF1 fq("fq", "[0] + [1]/x", 20., 250.);
gqs->Fit(&fq);
- printf("sy(Q) :: sm[%f] b[%f] 1/q0[%f]\n", sm, fq.GetParameter(1), (sm-fq.GetParameter(0))/fq.GetParameter(1));
+ //printf("sm=%f [0]=%f [1]=%f\n", 1.e-4*sm, fq.GetParameter(0), fq.GetParameter(1));
+ printf(" const Float_t sq0inv = %f; // [1/q0]\n", (sm-fq.GetParameter(0))/fq.GetParameter(1));
+ printf(" const Float_t sqb = %f; // [cm]\n", 1.e-4*fq.GetParameter(1));
}
//_______________________________________________________
void AliTRDclusterResolution::ProcessCenterPad()
{
+// Resolution as a function of y displacement from pad center
+// only for phi equal exB and clusters close to cathode wire plane
+// for ideal simulations time bins 4,5 and 6.
+
TObjArray *arr = (TObjArray*)fContainer->At(kCenter);
if(!arr) {
- AliWarning("Missing dy=f(y_d) container");
+ AliWarning("Missing dy=f(y | x, ly) container");
return;
}
+ Float_t s[AliTRDgeometry::kNlayer];
TF1 f("f", "gaus", -.5, .5);
+ TF1 fp("fp", "gaus", -3.5, 3.5);
+
+ TH1D *h1 = 0x0; TH2F *h2 = 0x0; TH3S *h3r=0x0, *h3p=0x0;
+ TObjArray *arrRes = (TObjArray*)fResults->At(kCenter);
+ TTree *t = (TTree*)arrRes->At(0);
+ TGraphErrors *gs = 0x0;
TAxis *ax = 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));
- }
- }
+ printf(" const Float_t lSy[6][24] = {\n {");
+ const Int_t nl = AliTRDgeometry::kNlayer;
+ for(Int_t il=0; il<nl; il++){
+ if(!(h3r = (TH3S*)arr->At(il))) continue;
+ if(!(h3p = (TH3S*)arr->At(nl+il))) continue;
+ gs = (TGraphErrors*)arrRes->At(il+1);
+ fLy = il;
+ for(Int_t ix=1; ix<=h3r->GetXaxis()->GetNbins(); ix++){
+ ax = h3r->GetXaxis(); ax->SetRange(ix, ix);
+ ax = h3p->GetXaxis(); ax->SetRange(ix, ix);
+ fX = ax->GetBinCenter(ix);
+ //printf(" x[%2d]=%4.2f\n", ix, fX);
+ for(Int_t iy=1; iy<=h3r->GetYaxis()->GetNbins(); iy++){
+ ax = h3r->GetYaxis(); ax->SetRange(iy, iy);
+ ax = h3p->GetYaxis(); ax->SetRange(iy, iy);
+ fY = ax->GetBinCenter(iy);
+ //printf(" y[%2d]=%5.2f\n", iy, fY);
+ // finish navigation in the HnSparse
+
+ h1 = (TH1D*)h3r->Project3D("z");
+ Int_t entries = (Int_t)h1->Integral();
+ if(entries < 50) continue;
+ //Adjust(&f, h1);
+ h1->Fit(&f, "QN");
+
+ // Fill sy,my=f(y_w,x,ly)
+ fR[0] = f.GetParameter(1); fR[1] = f.GetParError(1);
+ fR[2] = f.GetParameter(2); fR[3] = f.GetParError(2);
- // POSTPROCESS SPECTRA
- // Found correction for systematic deviation
- TF1 fprf("fprf", "[0]+[1]*sin([2]*x)", -.5, .5);
- fprf.SetParameter(0, 0.);
- fprf.SetParameter(1, 1.1E-2);
- fprf.SetParameter(2, -TMath::PiOver2()/0.25);
- printf(" const Float_t cy[AliTRDgeometry::kNlayer][3] = {\n");
- for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
- h1 = hym->ProjectionY("hym_pyy", ily+1, ily+1);
- // adjust errors
- for(Int_t ib=h1->GetNbinsX(); ib--;) h1->SetBinError(ib, 0.002);
- h1->Fit(&fprf, "Q");
- printf(" {%5.3e, %5.3e, %5.3e},\n", fprf.GetParameter(0), fprf.GetParameter(1), fprf.GetParameter(2));
+ h1 = (TH1D*)h3p->Project3D("z");
+ h1->Fit(&fp, "QN");
+ fP[0] = fp.GetParameter(1); fP[1] = fp.GetParError(1);
+ fP[2] = fp.GetParameter(2); fP[3] = fp.GetParError(2);
- if(!fCanvas) continue;
- fCanvas->cd();
- h1->SetMinimum(-0.02);h1->SetMaximum(0.02);h1->Draw("e1");
- h11 = hyp->ProjectionY("hyp_pyy", ily+1, ily+1);
- h11->Scale(.8/h11->Integral());
- h11->SetLineColor(kBlue); h11->Draw("csame");
- fCanvas->Modified(); fCanvas->Update();
- if(IsSaveAs())
- fCanvas->SaveAs(Form("Figures/ProcessCenterPad_M_Ly%d.gif", ily));
- else gSystem->Sleep(500);
- }
- printf(" };\n");
+ //printf("ly[%d] x[%3.1f] y[%+5.2f] m[%5.3f] s[%5.3f] \n", fLy, fX, fY, fR[0], fR[2]);
+ t->Fill();
- // Parameterization for sigma PRF
- TF1 fgaus("fgaus", "gaus", -.5, .5);
- printf(" const Float_t scy[AliTRDgeometry::kNlayer][4] = {\n");
- for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
- h1 = hys->ProjectionY("hys_pyy", ily+1, ily+1);
- // adjust errors
- for(Int_t ib=h1->GetNbinsX(); ib--;) h1->SetBinError(ib, 0.002);
-
- h1->Fit(&fgaus, "Q");
- printf(" {%5.3e, %5.3e, %5.3e, ", fgaus.GetParameter(0), fgaus.GetParameter(1), fgaus.GetParameter(2));
-
- // calculate mean sigma on the pad center distribution
- Float_t sy = 0.;
- h1 = hyp->ProjectionY("hyp_pyy", ily+1, ily+1);
- for(Int_t ix=1; ix<=h1->GetNbinsX(); ix++){
- sy += fgaus.Eval(h1->GetBinCenter(ix))*h1->GetBinContent(ix);
+
+ }
+ }
+ t->Draw("y:x>>h(24, 0., 2.4, 51, -.51, .51)",
+ Form("s[0]*(ly==%d&&abs(m[0])<1.e-1)", fLy),
+ "goff");
+ h2=(TH2F*)gROOT->FindObject("h");
+ f.FixParameter(1, 0.);
+ Int_t n = h2->GetXaxis()->GetNbins(); s[il]=0.;
+ printf(" {");
+ for(Int_t ix=1; ix<=n; ix++){
+ ax = h2->GetXaxis();
+ fX = ax->GetBinCenter(ix);
+ h1 = h2->ProjectionY("hCenPy", ix, ix);
+ //if((Int_t)h1->Integral() < 1.e-10) continue;
+ h1->Fit(&f, "QN");
+ s[il]+=f.GetParameter(2);
+ printf("%6.4f,%s", f.GetParameter(0), ix%6?" ":"\n ");
+ Int_t jx = gs->GetN();
+ gs->SetPoint(jx, fX, 1.e4*f.GetParameter(0));
+ gs->SetPointError(jx, 0., 0./*f.GetParError(0)*/);
}
- sy /= h1->GetEntries();
- printf("%5.3e},\n", sy);
+ printf("\b},\n");
+ s[il]/=n;
+
+ f.ReleaseParameter(2);
+
if(!fCanvas) continue;
- fCanvas->cd();
- h1->SetMinimum(0.01);h1->SetMaximum(0.04);h1->Draw("e1");
- h11 = hyp->ProjectionY("hyp_pyy", ily+1, ily+1);
- h11->Scale(1./h11->Integral());
- h11->SetLineColor(kBlue); h11->Draw("csame");
- fCanvas->Modified(); fCanvas->Update();
- if(IsSaveAs())
- fCanvas->SaveAs(Form("Figures/ProcessCenterPad_S_Ly%d.gif", ily));
- else gSystem->Sleep(500);
+ h2->Draw("lego2fb");
+ fCanvas->Modified(); fCanvas->Update();
+ if(IsSaveAs()) fCanvas->SaveAs(Form("Figures/ProcessCenter_ly[%d].gif", fLy));
+ else gSystem->Sleep(100);
}
printf(" };\n");
+ printf(" const Float_t lPRF[] = {"
+ "%5.3f, %5.3f, %5.3f, %5.3f, %5.3f, %5.3f};\n",
+ s[0], s[1], s[2], s[3], s[4], s[5]);
}
//_______________________________________________________
{
TObjArray *arr = (TObjArray*)fContainer->At(kSigm);
if(!arr){
- AliWarning("Missing dy=f(x_d, d_wire) container");
+ AliWarning("Missing dy=f(x_d, d_w) container");
return;
}
- TLinearFitter gs(1,"pol1");
- TF1 f("f", "gaus", -.5, .5);
- TAxis *ax = 0x0;
- TH1D *h1 = 0x0;
- TH2I *h2 = 0x0;
- TH1 *hFrame=0x0;
// init visualization
TGraphErrors *ggs = 0x0;
line->SetLineColor(kRed);line->SetLineWidth(2);
}
- 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);
- 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 it=1; it<=fAt->GetNbins(); it++){
- x = fAt->GetBinCenter(it); //[mm]
- Int_t idx = (id-1)*kNTB+it-1;
-
- // retrieve data histogram
- if(!(h2 = (TH2I*)arr->At(idx))) {
- AliWarning(Form("Missing histo at index idx[%3d] [id[%2d] it[%2d]] xd[%f] d[%f]\n", idx, id, it, x, d));
- continue;
- }
+ // init logistic support
+ TF1 f("f", "gaus", -.5, .5);
+ TLinearFitter gs(1,"pol1");
+ TH1 *hFrame=0x0;
+ TH1D *h1 = 0x0; TH3S *h3=0x0;
+ TAxis *ax = 0x0;
+ Double_t exb2 = fExB*fExB;
+
+ TTree *t = (TTree*)fResults->At(kSigm);
+ for(Int_t ix=0; ix<kNTB; ix++){
+ if(!(h3=(TH3S*)arr->At(ix))) continue;
+ fX = fAt->GetBinCenter(ix+1);
+
+ for(Int_t iz=1; iz<=h3->GetXaxis()->GetNbins(); iz++){
+ ax = h3->GetXaxis();
+ ax->SetRange(iz, iz);
+ fZ = ax->GetBinCenter(iz);
+ // reset visualization
if(fCanvas){
new(ggs) TGraphErrors();
ggs->SetMarkerStyle(7);
}
gs.ClearPoints();
- ax = h2->GetXaxis();
- for(Int_t ix=1; ix<=ax->GetNbins(); ix++){
- Float_t dydx = ax->GetBinCenter(ix);
+
+ for(Int_t ip=1; ip<=h3->GetYaxis()->GetNbins(); ip++){
+ ax = h3->GetYaxis();
+ ax->SetRange(ip, ip);
+ Double_t tgl = ax->GetBinCenter(ip);
+ // finish navigation in the HnSparse
+
//if(TMath::Abs(dydx)>0.18) continue;
- Double_t tgg = (dydx-fExB)/(1.+dydx*fExB);
+ Double_t tgg = (tgl-fExB)/(1.+tgl*fExB);
Double_t tgg2 = tgg*tgg;
- h1 = h2->ProjectionY("py", ix, ix);
- if(h1->GetEntries() < 100) continue;
+
+ h1 = (TH1D*)h3->Project3D("z");
+ Int_t entries = (Int_t)h1->Integral();
+ if(entries < 50) continue;
//Adjust(&f, h1);
- //printf("\tFit ix[%d] on %s entries [%d]\n", ix, h2->GetName(), (Int_t)h1->GetEntries());
h1->Fit(&f, "QN");
+
Double_t s2 = f.GetParameter(2)*f.GetParameter(2);
Double_t s2e = 2.*f.GetParameter(2)*f.GetParError(2);
// Fill sy^2 = f(tg^2(phi-a_L))
gs.AddPoint(&tgg2, s2, s2e);
if(!ggs) continue;
- Int_t ip = ggs->GetN();
- ggs->SetPoint(ip, tgg2, s2);
- ggs->SetPointError(ip, 0., s2e);
+ Int_t jp = ggs->GetN();
+ ggs->SetPoint(jp, tgg2, s2);
+ ggs->SetPointError(jp, 0., s2e);
}
if(gs.Eval()) continue;
// 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, .5*gs.GetParError(1)/TMath::Sqrt(sx));
+ fR[0] = gs.GetParameter(1)/* - x*x*exb2/12.*/;
+ if(fR[0]<0.) continue;
+ fR[0] = TMath::Sqrt(fR[0]);
+ fR[1] = .5*gs.GetParError(1)/fR[0];
// 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));
+ fR[2]= gs.GetParameter(0)/*-exb2*sx*/;
+ if(fR[1] <0.) continue;
+ fR[2] = TMath::Sqrt(fR[2]);
+ fR[3] = gs.GetParError(0)+exb2*exb2*gs.GetParError(1);
+ t->Fill();
if(!fCanvas) continue;
fCanvas->cd(); fCanvas->SetLogx(); //fCanvas->SetLogy();
}
ggs->Draw("pl"); line->Draw("l");
fCanvas->Modified(); fCanvas->Update();
- if(IsSaveAs()) fCanvas->SaveAs(Form("Figures/ProcessSigma_D%d_T%02d.gif", id, it));
+ if(IsSaveAs()) fCanvas->SaveAs(Form("Figures/ProcessSigma_z[%5.3f]_x[%5.3f].gif", fZ, fX));
else gSystem->Sleep(100);
- printf(" xd=%4.1f[cm] sx=%5.3e[cm] sy=%5.3e[cm]\n", x, TMath::Sqrt(sx), TMath::Sqrt(sy));
- }
- }
-
- 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(" xd=%4.1f[cm] sx=%5.3e[cm] sy=%5.3e[cm]\n", fX, TMath::Sqrt(fR[0]), TMath::Sqrt(fR[1]));
}
- printf("%5.3e}", hsx->GetBinContent(kND, iy));
- printf("%c\n", iy==(kNTB-1)?' ':',');
}
- printf(" };\n");
-
- 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;
}
{
TObjArray *arr = (TObjArray*)fContainer->At(kMean);
if(!arr){
- AliWarning("Missing dy=f(x_d, d_wire) container");
+ AliWarning("Missing dy=f(x_d, d_w) container");
return;
}
- TGraphErrors *gm = new TGraphErrors();
+
+ // init logistic support
TF1 f("f", "gaus", -.5, .5);
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;
+ TGraphErrors *gm = new TGraphErrors();
TH1 *hFrame=0x0;
+ TH1D *h1 = 0x0; TH3S *h3 =0x0;
+ TAxis *ax = 0x0;
- TObjArray *arrr = (TObjArray*)fResults->At(kMean);
- 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 it=1; it<=fAt->GetNbins(); it++){
- x = fAt->GetBinCenter(it); //[mm]
- Int_t idx = (id-1)*kNTB+it-1;
- if(!(h2 = (TH2I*)arr->At(idx))) {
- AliWarning(Form("Missing histo at index idx[%3d] [id[%2d] it[%2d]] xd[%f] d[%f]\n", idx, id, it, x, d));
- continue;
- }
+ TTree *t = (TTree*)fResults->At(kMean);
+ for(Int_t ix=0; ix<kNTB; ix++){
+ if(!(h3=(TH3S*)arr->At(ix))) continue;
+ fX = fAt->GetBinCenter(ix+1);
+
+ for(Int_t iz=1; iz<=h3->GetXaxis()->GetNbins(); iz++){
+ ax = h3->GetXaxis();
+ ax->SetRange(iz, iz);
+ fZ = ax->GetBinCenter(iz);
+ // reset fitter
new(gm) TGraphErrors();
gm->SetMarkerStyle(7);
- ax = h2->GetXaxis();
- for(Int_t ix=1; ix<=ax->GetNbins(); ix++){
- Double_t dydx = ax->GetBinCenter(ix);
- h1 = h2->ProjectionY("py", ix, ix);
- if(h1->GetEntries() < 200) continue;
+
+ for(Int_t ip=1; ip<=h3->GetYaxis()->GetNbins(); ip++){
+ ax = h3->GetYaxis();
+ ax->SetRange(ip, ip);
+ Double_t tgl = ax->GetBinCenter(ip);
+ // finish navigation in the HnSparse
+
+ h1 = (TH1D*)h3->Project3D("z");
+ Int_t entries = (Int_t)h1->Integral();
+ if(entries < 50) continue;
//Adjust(&f, h1);
h1->Fit(&f, "QN");
// Fill <Dy> = f(dydx - h*dzdx)
- Int_t ip = gm->GetN();
- gm->SetPoint(ip, dydx, f.GetParameter(1));
- gm->SetPointError(ip, 0., f.GetParError(1));
- ip++;
+ Int_t jp = gm->GetN();
+ gm->SetPoint(jp, tgl, f.GetParameter(1));
+ gm->SetPointError(jp, 0., f.GetParError(1));
}
if(gm->GetN()<4) continue;
gm->Fit(&line, "QN");
- 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);
- hdy->SetBinContent(id, it, dy);
+ fR[0] = line.GetParameter(1); // dx
+ fR[1] = line.GetParError(1);
+ fR[2] = line.GetParameter(0) + fExB*fR[0]; // xs = dy - tg(a_L)*dx
+ t->Fill();
if(!fCanvas) continue;
fCanvas->cd();
} 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));
+ if(IsSaveAs()) fCanvas->SaveAs(Form("Figures/ProcessMean_Z[%5.3f]_X[%5.3f].gif", fZ, fX));
else gSystem->Sleep(100);
- printf(" xd=%4.1f[cm] dx=%5.3e[cm] dy=%5.3e[cm] xs=%5.3e[cm]\n", x, dx, dy, xs);
- }
- }
-
- // 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);
+ printf(" xd=%4.2f[cm] dx=%5.3e[cm] dy=%5.3e[cm]\n", fX, fR[0], fR[2]);
}
- 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);
}
+
+ // draw shift results
+ //t->Draw("z:x>>h(24, 0, 2.4, 25, 0, 2.5)", "dx*(abs(dx)<1.e-2)", "lego2fb");
+ //t->Draw("z:x>>h(24, 0, 2.4, 25, 0, 2.5)", "dy*(abs(dx)<1.e-2)", "lego2fb");
- // 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