\r
ClassImp(AliTRDrecoTask)\r
\r
-Float_t AliTRDrecoTask::fgPt[AliTRDrecoTask::fgNPt] = {0.};\r
+Float_t AliTRDrecoTask::fgPt[AliTRDrecoTask::fgNPt+1] = {0.};\r
TTreeSRedirector* AliTRDrecoTask::fgDebugStream(NULL);\r
+TH1* AliTRDrecoTask::fgProjector(NULL);\r
//_______________________________________________________\r
AliTRDrecoTask::AliTRDrecoTask()\r
: AliAnalysisTaskSE()\r
// Generic task destructor\r
\r
AliDebug(2, Form(" Ending task %s[%s]", GetName(), GetTitle()));\r
- if(fgDebugStream){ \r
+ if(fgDebugStream){\r
delete fgDebugStream;\r
fgDebugStream = NULL;\r
}\r
\r
+ if(fgProjector){\r
+ delete fgProjector;\r
+ fgProjector = NULL;\r
+ }\r
+\r
if(fPlotFuncList){\r
fPlotFuncList->Delete();\r
delete fPlotFuncList;\r
{\r
switch(fNpt){\r
case fgNPt:\r
- fgPt[0]=0.5;\r
+ fgPt[0]=0.3;\r
for(Int_t j(1); j<=fgNPt; j++) fgPt[j]=fgPt[j-1]+(TMath::Exp(j*j*2.e-3)-1.);\r
AliDebug(2, "Using debug momentum segmentation");\r
break;\r
}\r
if(!gFile->cd(dir)){\r
AliWarning(Form("Couldn't cd to %s in %s.", dir, file));\r
+ gFile->Close();\r
return kFALSE;\r
}\r
- TObjArray *o = NULL;\r
- if(!(o = (TObjArray*)gDirectory->Get(GetName()))){\r
+ if(!(fContainer = (TObjArray*)gDirectory->Get(GetName()))){\r
AliWarning("Missing histogram container.");\r
+ gFile->Close();\r
return kFALSE;\r
}\r
- fContainer = (TObjArray*)o->Clone(GetName());\r
gFile->Close();\r
return kTRUE;\r
}\r
}\r
if(!gFile->cd(dir)){\r
AliWarning(Form("Couldn't cd to %s in %s.", dir, file));\r
+ gFile->Close();\r
return kFALSE;\r
}\r
TObjArray *info = NULL;\r
if(!(info = (TObjArray*)gDirectory->Get("TRDinfoGen"))){\r
AliWarning("Missing TRDinfoGen container.");\r
+ gFile->Close();\r
return kFALSE;\r
}\r
- TObjArray *dets = (TObjArray*)info->FindObject("Chambers Status");\r
- if(!dets){\r
+ fDets = (TObjArray*)((TObjArray*)info->FindObject("Chambers Status"))->Clone();\r
+ if(!fDets){\r
if(!info->At(4) || strcmp("TObjArray", info->At(4)->IsA()->GetName())) AliError("Looking for old style chamber status map. Failed.");\r
else {\r
AliWarning("Looking for old style chamber status map.");\r
- TObjArray *vdets = (TObjArray*)info->At(4);\r
- fDetsV = (TObjArray*)vdets->Clone();\r
+ fDetsV = (TObjArray*)((TObjArray*)info->At(4))->Clone();\r
}\r
- } else fDets = (TObjArray*)dets->Clone();\r
+ }\r
gFile->Close();\r
+ info->Delete(); delete info;\r
return kTRUE;\r
}\r
\r
if(AliTRDgeometry::GetLayer(ci->GetDetector()) != ly) continue;\r
ci->Draw(opt);\r
}\r
- \r
- Float_t dsm = TMath::TwoPi()/AliTRDgeometry::kNsector;\r
- Float_t xmed=0.;\r
- if(strcmp(opt, "pad")==0) xmed=38.;\r
- TLatex *sm = new TLatex(); sm->SetTextAlign(22);sm->SetTextColor(kBlack); sm->SetTextFont(32);sm->SetTextSize(0.03);\r
- for(Int_t is(0); is<AliTRDgeometry::kNsector; is++) sm->DrawLatex(xmed, -TMath::Pi()+(is+0.5)*dsm, Form("%02d", is>=9?(is-9):(is+9)));\r
}\r
\r
//_______________________________________________________\r
}\r
\r
//________________________________________________________\r
-Float_t AliTRDrecoTask::GetMeanStat(TH1 *h, Float_t cut, Option_t *opt)\r
+Float_t AliTRDrecoTask::GetMeanStat(TH1 *h, Float_t cut, Int_t opt, Float_t *sigma)\r
{\r
-// return mean number of entries/bin of histogram "h"\r
-// if option "opt" is given the following values are accepted:\r
-// "<" : consider only entries less than "cut"\r
-// ">" : consider only entries greater than "cut"\r
-\r
- //Int_t dim(h->GetDimension());\r
+// Return mean number of entries/bin of histogram "h".\r
+// If optionally sigma is allocated than it is also filled with sigma paramter of the gauss fit \r
+//\r
+// Option "opt" is given the following values are accepted:\r
+// -1 : consider only entries less than "cut"\r
+// 1 : consider only entries greater than "cut"\r
+// 0 : no "cut" [dafault]\r
+// Error codes\r
+// -999. : statistics too low [20]\r
+// -998. : fit failed\r
+\r
+ const Int_t kvd(100000);\r
+ Float_t v[kvd];\r
Int_t nbx(h->GetNbinsX()), nby(h->GetNbinsY()), nbz(h->GetNbinsZ());\r
- Double_t sum(0.); Int_t n(0);\r
- for(Int_t ix(1); ix<=nbx; ix++)\r
- for(Int_t iy(1); iy<=nby; iy++)\r
+ Int_t nv(0); Float_t xmin(1.e5), xmax(-xmin);\r
+ for(Int_t ix(1); ix<=nbx; ix++){\r
+ for(Int_t iy(1); iy<=nby; iy++){\r
for(Int_t iz(1); iz<=nbz; iz++){\r
- if(strcmp(opt, "")==0){sum += h->GetBinContent(ix, iy, iz); n++;}\r
- else{\r
- if(strcmp(opt, "<")==0) {\r
- if(h->GetBinContent(ix, iy, iz)<cut) {sum += h->GetBinContent(ix, iy, iz); n++;}\r
- } else if(strcmp(opt, ">")==0){\r
- if(h->GetBinContent(ix, iy, iz)>cut) {sum += h->GetBinContent(ix, iy, iz); n++;}\r
- } else {sum += h->GetBinContent(ix, iy, iz); n++;}\r
+ Float_t c = h->GetBinContent(ix, iy, iz);\r
+ if(opt*(c-cut)<0.) continue;\r
+ v[nv++] = c;\r
+ if(c<xmin) xmin = c;\r
+ if(c>xmax) xmax = c;\r
+ if(nv==kvd){\r
+ printf("W - AliTRDrecoTask::GetMeanStat() :: Unreliable results for %s[%s]. Statical allocation exceeded.\n", h->GetName(), h->GetTitle());\r
+ break;\r
}\r
}\r
- return n>0?sum/n:0.;\r
+ if(nv==kvd) break;\r
+ }\r
+ if(nv==kvd) break;\r
+ }\r
+ if(nv<10){\r
+ //printf("W - AliTRDrecoTask::GetMeanStat() :: Failed for %s[%s]. Statical undefined [%d].\n", h->GetName(), h->GetTitle(), nv);\r
+ return -999.;\r
+ }\r
+ if(fgProjector) delete fgProjector;\r
+ fgProjector = new TH1F("hProjector", "", 20, 0.5*(3*xmin-xmax), 0.5*(3*xmax - xmin));\r
+ for(Int_t iv(0); iv<nv; iv++) fgProjector->Fill(v[iv]);\r
+ TF1 f("f", "gaus", xmin, xmax);\r
+ f.SetParameter(0, fgProjector->Integral());\r
+ f.SetParameter(1, fgProjector->GetMean()); f.SetParLimits(1, xmin, xmax);\r
+ f.SetParameter(2, fgProjector->GetRMS());\r
+ if(fgProjector->Fit(&f, "WQ0", "goff")) return -998.;\r
+ if(sigma) *sigma = f.GetParameter(2);\r
+ return f.GetParameter(1);\r
+}\r
+\r
+//________________________________________________________\r
+Int_t AliTRDrecoTask::Rebin(TH2 *h, Int_t n, Int_t rebinX[], Int_t rebinY[], Int_t nstat)\r
+{\r
+// Rebin histo "h" according to "rebinning" strategy "n" steps such to obtain mean statistics per bin over "nstat"\r
+\r
+ Int_t irebin(0);\r
+ while(irebin<n && GetMeanStat(h, .5, 1)<nstat){\r
+ h->Rebin2D(rebinX[irebin], rebinY[irebin]);\r
+ irebin++;\r
+ }\r
+ return irebin;\r
}\r
\r
//________________________________________________________\r
:TNamed()\r
,fH(NULL)\r
,fNrebin(0)\r
- ,fRebinX(NULL)\r
- ,fRebinY(NULL)\r
{\r
// constructor\r
+ fRebin[0] = NULL;fRebin[1] = NULL;\r
memset(fAx, 0, 3*sizeof(Int_t));\r
memset(fRange, 0, 4*sizeof(Float_t));\r
}\r
{\r
// destructor\r
if(fH) delete fH;\r
+ if(fRebin[0]) delete [] fRebin[0];\r
+ if(fRebin[1]) delete [] fRebin[1];\r
}\r
\r
//________________________________________________________\r
if(this == &rhs) return *this;\r
\r
TNamed::operator=(rhs);\r
- if(fNrebin){fNrebin=0; delete [] fRebinX; delete [] fRebinY;}\r
- if(rhs.fNrebin) SetRebinStrategy(rhs.fNrebin, rhs.fRebinX, rhs.fRebinY);\r
+ if(fNrebin){fNrebin=0; delete [] fRebin[0]; delete [] fRebin[1];}\r
+ if(rhs.fNrebin) SetRebinStrategy(rhs.fNrebin, rhs.fRebin[0], rhs.fRebin[1]);\r
memcpy(fAx, rhs.fAx, 3*sizeof(Int_t));\r
memcpy(fRange, rhs.fRange, 4*sizeof(Float_t));\r
if(fH) delete fH;\r
}\r
\r
//________________________________________________________\r
-Double_t AliTRDrecoTask::AliTRDrecoProjection::GetTrendValue(const Int_t mid, Double_t *m, Double_t *s) const\r
+Double_t AliTRDrecoTask::AliTRDrecoProjection::GetTrendValue(const Int_t mid, Double_t *e, Double_t *s, Double_t *se) const\r
{\r
// Return result of fitting the main distribution (represented on the z axis) with the function selected\r
// "mid". Optionally return the Mean and RMS of the distribution pointing to "m" and "s"\r
Int_t ne((Int_t)h1s->Integral());\r
if(ne<30){\r
AliDebug(1, Form("Statistics too low[%2d] in %s", ne, GetName()));\r
+ delete h1s;\r
return -999.;\r
}\r
TAxis *az(h1s->GetXaxis());\r
- Float_t vm(h1s->GetMean()), v(vm), ve(h1s->GetRMS());\r
+ Float_t mn(h1s->GetMean()), rms(h1s->GetRMS()),\r
+ v(mn), // main trending value (mean, mu, MPV)\r
+ ve(rms),// dispersion (RMS, sigma, landau 2nd param)\r
+ ev(0.), // error on v\r
+ eve(0.);// error on ve\r
if(mid==1){\r
TF1 fg("fg", "gaus", az->GetXmin(), az->GetXmax());\r
- fg.SetParameter(0, Float_t(ne)); fg.SetParameter(1, vm); fg.SetParameter(2, ve);\r
+ fg.SetParameter(0, Float_t(ne)); fg.SetParameter(1, mn); fg.SetParameter(2, rms);\r
h1s->Fit(&fg, "WQ0");\r
- v = fg.GetParameter(1);\r
+ v = fg.GetParameter(1); ev = fg.GetParError(1);\r
+ ve= fg.GetParameter(2);\r
} else if (mid==2) {\r
TF1 fl("fl", "landau", az->GetXmin(), az->GetXmax());\r
- fl.SetParameter(0, Float_t(ne)); fl.SetParameter(1, vm); fl.SetParameter(2, ve);\r
+ fl.SetParameter(0, Float_t(ne)); fl.SetParameter(1, mn); fl.SetParameter(2, rms);\r
h1s->Fit(&fl, "WQ0");\r
- v = fl.GetMaximumX();\r
+ v = fl.GetMaximumX(); ev = fl.GetParError(1);\r
+ ve= fl.GetParameter(2);\r
}\r
- if(m) *m = vm;\r
- if(s) *s = ve;\r
- AliDebug(2, Form("%s[%d]:: %f {%f %f} Entries[%d]", fH->GetName(), mid, v, m?(*m):0., s?(*s):0., (Int_t)h1s->Integral()));\r
+ if(e) *e = ev;\r
+ if(s) *s = ve;\r
+ if(se) *se = eve;\r
+ AliDebug(2, Form("[%d] %s(%4d) = M{%f+-%f} S{%f+-%f}", mid, fH->GetName(), (Int_t)h1s->Integral(), v, ev, ve, eve));\r
\r
+ delete h1s;\r
return v;\r
}\r
\r
//________________________________________________________\r
-TH2* AliTRDrecoTask::AliTRDrecoProjection::Projection2Dbin(Int_t bin)\r
+TH2* AliTRDrecoTask::AliTRDrecoProjection::Projection2Dbin(Int_t bin, Bool_t mc)\r
{\r
// dumb 2D projection for bin including under/over flow. Default all [bin==-1]\r
\r
ax->GetNbins(), ax->GetXmin(), ax->GetXmax(),\r
ay->GetNbins(), ay->GetXmin(), ay->GetXmax());\r
else h2 = new TH2F(Form("%s%d_2D", fH->GetName(), bin),\r
- Form("%s | #it{%4.2f<=p_{t}[GeV/c]<%4.2f};%s;%s;Entries", fH->GetTitle(),\r
- bin?fgPt[bin-1]:0., bin==nbins?9.99:fgPt[bin], ax->GetTitle(), ay->GetTitle()),\r
+ Form("%s | #it{%4.2f<=p_{t}^{%s}[GeV/c]<%4.2f};%s;%s;Entries", fH->GetTitle(),\r
+ bin?fgPt[bin-1]:0., mc?"MC":"", bin>nbins?99.99:fgPt[bin], ax->GetTitle(), ay->GetTitle()),\r
ax->GetNbins(), ax->GetXmin(), ax->GetXmax(),\r
ay->GetNbins(), ay->GetXmin(), ay->GetXmax());\r
- printf("build %s\n", h2->GetName());\r
for(Int_t ix(1); ix<=ax->GetNbins(); ix++){\r
for(Int_t iy(1); iy<=ay->GetNbins(); iy++){\r
Int_t ibin = h2->GetBin(ix, iy);\r
hyx = (TH2D*)h2s->Clone();\r
hyx->SetName(Form("%sEn", fH->GetName()));\r
}\r
- Int_t irebin(0), dxBin(1), dyBin(1);\r
- while(irebin<fNrebin && (AliTRDrecoTask::GetMeanStat(h2s, .5, ">")<nstat)){\r
- h2s->Rebin2D(fRebinX[irebin], fRebinY[irebin]);\r
- dxBin*=fRebinX[irebin];dyBin*=fRebinY[irebin];\r
- irebin++;\r
- }\r
+ Int_t irebin(Rebin(h2s, fNrebin, fRebin[0], fRebin[1], nstat)), dxBin(1), dyBin(1);\r
+ for(Int_t ir(0); ir<irebin; ir++){dxBin*=fRebin[0][ir]; dyBin*=fRebin[1][ir];}\r
Int_t nx(h2s->GetNbinsX()), ny(h2s->GetNbinsY());\r
delete h2s;\r
if(mid<0) return NULL;\r
{\r
// define rebinning strategy for this projection\r
fNrebin = n;\r
- fRebinX = new Int_t[n]; memcpy(fRebinX, rebx, n*sizeof(Int_t));\r
- fRebinY = new Int_t[n]; memcpy(fRebinY, reby, n*sizeof(Int_t));\r
+ fRebin[0] = new Int_t[n]; memcpy(fRebin[0], rebx, n*sizeof(Int_t));\r
+ fRebin[1] = new Int_t[n]; memcpy(fRebin[1], reby, n*sizeof(Int_t));\r
}\r
\r
//________________________________________________________\r
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff