#include "TList.h"\r
#include "TMap.h"\r
#include "TH1.h"\r
+#include "TH2.h"\r
+#include "TH3.h"\r
#include "TF1.h"\r
#include "TObjArray.h"\r
#include "TDirectory.h"\r
}\r
Float_t xmin(0.), xmax(0.);\r
TBox *gdet = new TBox();\r
- gdet->SetLineWidth(kBlack);gdet->SetFillColor(kBlack);\r
+ gdet->SetLineColor(kBlack);gdet->SetFillColor(kBlack);\r
Int_t style[] = {0, 3003};\r
for(Int_t idet(0); idet<540; idet++){\r
if(idet%6 != ly) continue;\r
f->SetParLimits(5, 0., 1.e2);\r
f->SetParameter(5, 2.e-1);\r
}\r
+\r
+\r
+//________________________________________________________\r
+void AliTRDrecoTask::SetNormZ(TH2 *h2, Int_t bxmin, Int_t bxmax, Int_t bymin, Int_t bymax, Float_t thr)\r
+{\r
+// Normalize histo content to the mean value in the range specified by bin ranges\r
+// [bxmin, bxmax] on the x axis and [bymin, bymax] on the y axis.\r
+// Optionally a threshold "thr" can be specified to disregard entries with no meaning\r
+\r
+ Float_t s = 0., c=0.; Int_t is(0);\r
+ for(Int_t ix(bxmin); ix<=(bxmax>0?bxmax:(h2->GetXaxis()->GetNbins())); ix++){\r
+ for(Int_t iy(bymin); iy<=(bymax>0?bymax:(h2->GetYaxis()->GetNbins())); iy++){\r
+ if((c = h2->GetBinContent(ix, iy))<thr) continue;\r
+ s += c; is++;\r
+ }\r
+ }\r
+ s/= is?is:1;\r
+ for(Int_t ix(1); ix<=h2->GetXaxis()->GetNbins(); ix++){\r
+ for(Int_t iy(1); iy<=h2->GetYaxis()->GetNbins(); iy++){\r
+ if((c = h2->GetBinContent(ix, iy))<thr) h2->SetBinContent(ix, iy, thr-100);\r
+ else h2->SetBinContent(ix, iy, 100.*(c/s-1.));\r
+ }\r
+ }\r
+}\r
+\r
+//________________________________________________________\r
+void AliTRDrecoTask::SetRangeZ(TH2 *h2, Float_t min, Float_t max, Float_t thr)\r
+{\r
+// Set range on Z axis such to avoid outliers\r
+\r
+ Float_t c(0.), dz(1.e-3*(max-min));\r
+ for(Int_t ix(1); ix<=h2->GetXaxis()->GetNbins(); ix++){\r
+ for(Int_t iy(1); iy<=h2->GetYaxis()->GetNbins(); iy++){\r
+ if((c = h2->GetBinContent(ix, iy))<thr) continue;\r
+ if(c<=min) h2->SetBinContent(ix, iy, min+dz);\r
+ }\r
+ }\r
+ h2->GetZaxis()->SetRangeUser(min, max);\r
+}\r
+\r
+//________________________________________________________\r
+Float_t AliTRDrecoTask::GetMeanStat(TH1 *h, Float_t cut, Option_t *opt)\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
+ 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
+ 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
+ }\r
+ }\r
+ return n>0?sum/n:0.;\r
+}\r
+\r
+//________________________________________________________\r
+AliTRDrecoTask::AliTRDrecoProjection::AliTRDrecoProjection()\r
+ :TNamed()\r
+ ,fH(NULL)\r
+ ,fNrebin(0)\r
+ ,fRebinX(NULL)\r
+ ,fRebinY(NULL)\r
+{\r
+ // constructor\r
+ memset(fAx, 0, 3*sizeof(Int_t));\r
+ memset(fRange, 0, 4*sizeof(Float_t));\r
+}\r
+\r
+//________________________________________________________\r
+AliTRDrecoTask::AliTRDrecoProjection::~AliTRDrecoProjection()\r
+{\r
+ // destructor\r
+ if(fH) delete fH;\r
+}\r
+\r
+//________________________________________________________\r
+void AliTRDrecoTask::AliTRDrecoProjection::Build(const Char_t *n, const Char_t *t, Int_t ix, Int_t iy, Int_t iz, TAxis *aa[])\r
+{\r
+// check and build (if neccessary) projection determined by axis "ix", "iy" and "iz"\r
+ if(!aa[ix] || !aa[iy] || !aa[iz]) return;\r
+ TAxis *ax(aa[ix]), *ay(aa[iy]), *az(aa[iz]);\r
+ // check ax definiton to protect against older versions of the data\r
+ if(ax->GetNbins()<=0 || (ax->GetXmax()-ax->GetXmin())<=0.){\r
+ AliWarning(Form("Wrong definition of axis[%d] \"%s\"[%d](%f %f).", ix, ax->GetTitle(), ax->GetNbins(), ax->GetXmin(), ax->GetXmax()));\r
+ return;\r
+ }\r
+ if(ay->GetNbins()<=0 || (ay->GetXmax()-ay->GetXmin())<=0.){\r
+ AliWarning(Form("Wrong definition of axis[%d] \"%s\"[%d](%f %f).", ix, ay->GetTitle(), ay->GetNbins(), ay->GetXmin(), ay->GetXmax()));\r
+ return;\r
+ }\r
+ if(az->GetNbins()<=0 || (az->GetXmax()-az->GetXmin())<=0.){\r
+ AliWarning(Form("Wrong definition of axis[%d] \"%s\"[%d](%f %f).", ix, az->GetTitle(), az->GetNbins(), az->GetXmin(), az->GetXmax()));\r
+ return;\r
+ }\r
+ SetNameTitle(n,t);\r
+ fH = new TH3I(n, Form("%s;%s;%s;%s", t, ax->GetTitle(), ay->GetTitle(), az->GetTitle()),\r
+ ax->GetNbins(), ax->GetXmin(), ax->GetXmax(),\r
+ ay->GetNbins(), ay->GetXmin(), ay->GetXmax(),\r
+ az->GetNbins(), az->GetXmin(), az->GetXmax());\r
+ fAx[0] = ix; fAx[1] = iy; fAx[2] = iz;\r
+ fRange[0] = az->GetXmin()/3.; fRange[1] = az->GetXmax()/3.;\r
+ AliDebug(2, Form("H3(%s, %s) :: %s[%3d %4.2f %4.2f]%s[%3d %4.2f %4.2f]%s[%3d %4.2f %4.2f]", n, t,\r
+ ax->GetTitle(), ax->GetNbins(), ax->GetXmin(), ax->GetXmax(),\r
+ ay->GetTitle(), ay->GetNbins(), ay->GetXmin(), ay->GetXmax(),\r
+ az->GetTitle(), az->GetNbins(), az->GetXmin(), az->GetXmax()));\r
+}\r
+\r
+//________________________________________________________\r
+AliTRDrecoTask::AliTRDrecoProjection& AliTRDrecoTask::AliTRDrecoProjection::operator=(const AliTRDrecoProjection& rhs)\r
+{\r
+// copy projections\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
+ memcpy(fAx, rhs.fAx, 3*sizeof(Int_t));\r
+ memcpy(fRange, rhs.fRange, 4*sizeof(Float_t));\r
+ if(fH) delete fH;\r
+ if(rhs.fH) fH=(TH3I*)rhs.fH->Clone(Form("%s_CLONE", rhs.fH->GetName()));\r
+ return *this;\r
+}\r
+\r
+//________________________________________________________\r
+AliTRDrecoTask::AliTRDrecoProjection& AliTRDrecoTask::AliTRDrecoProjection::operator+=(const AliTRDrecoProjection& other)\r
+{\r
+// increment projections\r
+ if(!fH || !other.fH) return *this;\r
+ AliDebug(2, Form("%s+=%s [%s+=%s]", GetName(), other.GetName(), fH->GetName(), (other.fH)->GetName()));\r
+ fH->Add(other.fH);\r
+ return *this;\r
+}\r
+\r
+//________________________________________________________\r
+void AliTRDrecoTask::AliTRDrecoProjection::Increment(Int_t bin[], Double_t v)\r
+{\r
+// increment bin with value "v" pointed by general coord in "bin"\r
+ if(!fH) return;\r
+ AliDebug(4, Form(" %s[%2d]", fH->GetName(), Int_t(v)));\r
+ fH->AddBinContent(fH->GetBin(bin[fAx[0]],bin[fAx[1]],bin[fAx[2]]), Int_t(v));\r
+}\r
+\r
+//________________________________________________________\r
+TH2* AliTRDrecoTask::AliTRDrecoProjection::Projection2D(const Int_t nstat, const Int_t ncol, const Int_t mid, Bool_t del)\r
+{\r
+// build the 2D projection and adjust binning\r
+\r
+ const Char_t *title[] = {"Mean", "#mu", "MPV"};\r
+ if(!fH) return NULL;\r
+ TAxis *ax(fH->GetXaxis()), *ay(fH->GetYaxis()), *az(fH->GetZaxis());\r
+ TH2D *h2s(NULL), *hyx(NULL);\r
+ if(!(h2s = (TH2D*)fH->Project3D("yx"))) return NULL;\r
+ // save a copy of the original distribution\r
+ if(!del){\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 nx(h2s->GetNbinsX()), ny(h2s->GetNbinsY());\r
+ delete h2s;\r
+ if(mid<0) return NULL;\r
+\r
+ // start projection\r
+ TH1 *h(NULL); Int_t n(0);\r
+ Float_t dz=(fRange[1]-fRange[1])/ncol;\r
+ TString titlez(az->GetTitle()); TObjArray *tokenTitle(titlez.Tokenize(" "));\r
+ Int_t nt(tokenTitle->GetEntriesFast());\r
+ TH2 *h2 = new TH2F(Form("%s_2D", fH->GetName()),\r
+ Form("%s;%s;%s;%s(%s) %s", fH->GetTitle(), ax->GetTitle(), ay->GetTitle(), title[mid], nt>0?(*tokenTitle)[0]->GetName():"", nt>1?(*tokenTitle)[1]->GetName():""),\r
+ nx, ax->GetXmin(), ax->GetXmax(), ny, ay->GetXmin(), ay->GetXmax());\r
+ h2->SetContour(ncol);\r
+ h2->GetZaxis()->CenterTitle();\r
+ h2->GetZaxis()->SetTitleOffset(1.4);\r
+ h2->GetZaxis()->SetRangeUser(fRange[0], fRange[1]);\r
+ AliDebug(2, Form("%s[%s] nx[%d] ny[%d]", h2->GetName(), h2->GetTitle(), nx, ny));\r
+ for(Int_t iy(0); iy<ny; iy++){\r
+ for(Int_t ix(0); ix<nx; ix++){\r
+ h = fH->ProjectionZ(Form("%s_z", h2->GetName()), ix*dxBin+1, (ix+1)*dxBin+1, iy*dyBin+1, (iy+1)*dyBin+1);\r
+ Int_t ne((Int_t)h->Integral());\r
+ if(ne<nstat/2){\r
+ h2->SetBinContent(ix+1, iy+1, -999);\r
+ h2->SetBinError(ix+1, iy+1, 1.);\r
+ n++;\r
+ }else{\r
+ Float_t v(h->GetMean()), ve(h->GetRMS());\r
+ if(mid==1){\r
+ TF1 fg("fg", "gaus", az->GetXmin(), az->GetXmax());\r
+ fg.SetParameter(0, Float_t(ne)); fg.SetParameter(1, v); fg.SetParameter(2, ve);\r
+ h->Fit(&fg, "WQ");\r
+ v = fg.GetParameter(1); 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, v); fl.SetParameter(2, ve);\r
+ h->Fit(&fl, "WQ");\r
+ v = fl.GetMaximumX(); ve = fl.GetParameter(2);\r
+/* TF1 fgle("gle", "[0]*TMath::Landau(x, [1], [2], 1)*TMath::Exp(-[3]*x/[1])", az->GetXmin(), az->GetXmax());\r
+ fgle.SetParameter(0, fl.GetParameter(0));\r
+ fgle.SetParameter(1, fl.GetParameter(1));\r
+ fgle.SetParameter(2, fl.GetParameter(2));\r
+ fgle.SetParameter(3, 1.);fgle.SetParLimits(3, 0., 5.);\r
+ h->Fit(&fgle, "WQ");\r
+ v = fgle.GetMaximumX(); ve = fgle.GetParameter(2);*/\r
+ }\r
+ if(v<fRange[0]) h2->SetBinContent(ix+1, iy+1, fRange[0]+0.1*dz);\r
+ else h2->SetBinContent(ix+1, iy+1, v);\r
+ h2->SetBinError(ix+1, iy+1, ve);\r
+ }\r
+ }\r
+ }\r
+ if(h) delete h;\r
+ if(n==nx*ny){delete h2; h2=NULL;} // clean empty projections\r
+ return h2;\r
+}\r
+\r
+//________________________________________________________\r
+void AliTRDrecoTask::AliTRDrecoProjection::SetRebinStrategy(Int_t n, Int_t rebx[], Int_t reby[])\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
+}\r
+\r
+\r
+\r
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff