//
SetNameTitle("TRDresolution", "TRD spatial and momentum resolution");
SetSegmentationLevel();
+ memset(fXcorr, 0, AliTRDgeometry::kNdet*30*sizeof(Float_t));
}
//________________________________________________________
InitFunctorList();
SetSegmentationLevel();
+ memset(fXcorr, 0, AliTRDgeometry::kNdet*30*sizeof(Float_t));
DefineOutput(kClToTrk, TObjArray::Class()); // cluster2track
DefineOutput(kClToMC, TObjArray::Class()); // cluster2mc
fCl->Delete();
fMCcl->Delete();
+ if(!HasLoadCorrection()) LoadCorrection("correction.lst");
AliTRDrecoTask::UserExec(opt);
}
AliTRDcluster *c = NULL;
fTracklet->ResetClusterIter(kFALSE);
while((c = fTracklet->PrevCluster())){
- Float_t xyz[3] = {c->GetX(), c->GetY(), c->GetZ()};
+ Float_t xyz[3] = {c->GetX()+fXcorr[c->GetDetector()][c->GetLocalTimeBin()], c->GetY(), c->GetZ()};
clusters[np].SetCharge(tilt);
clusters[np].SetClusterType(0);
clusters[np].SetVolumeID(ily);
AliTRDcluster *c = NULL;
fTracklet->ResetClusterIter(kFALSE);
while((c = fTracklet->PrevCluster())){
- Float_t xc = c->GetX();
+ Float_t xc = c->GetX()+fXcorr[c->GetDetector()][c->GetLocalTimeBin()];
Float_t yc = c->GetY();
Float_t zc = c->GetZ();
Float_t dx = x0 - xc;
Int_t istk = geo->GetStack(c->GetDetector());
AliTRDpadPlane *pp = geo->GetPadPlane(ily, istk);
Float_t row0 = pp->GetRow0();
- Float_t d = row0 - zt + pp->GetAnodeWireOffset();
+ Float_t d = row0 - zt + pp->GetAnodeWireOffset();
d -= ((Int_t)(2 * d)) / 2.0;
if (d > 0.25) d = 0.5 - d;
tt.ResetClusterIter(kFALSE);
while((c = tt.PrevCluster())){
Float_t q = TMath::Abs(c->GetQ());
- x = c->GetX(); y = c->GetY();z = c->GetZ();
+ x = c->GetX()+fXcorr[c->GetDetector()][c->GetLocalTimeBin()]; y = c->GetY();z = c->GetZ();
dx = x0 - x;
ymc= y0 - dx*dydx0;
zmc= z0 - dx*dzdx0;
p=cOut->cd(3);
p->SetRightMargin(0.06);p->SetTopMargin(0.06);
xy[0]=-.5; xy[1]=-0.5; xy[2]=fgkNresYsegm[fSegmentLevel]-.5; xy[3]=2.5;
- GetGraphArray(xy, kCluster, 1, 1);
+ if(!GetGraphArray(xy, kCluster, 1, 1)) return;
p=cOut->cd(4);
p->SetRightMargin(0.16);p->SetTopMargin(0.06);
p=cOut->cd(6);
p->SetRightMargin(0.06);p->SetTopMargin(0.06);
xy[0]=-.5; xy[1]=-0.5; xy[2]=fgkNresYsegm[fSegmentLevel]-.5; xy[3]=2.5;
- GetGraphArray(xy, kTrack, 1, 1);
+ if(!GetGraphArray(xy, kTrack, 1, 1)) return;
cOut->SaveAs(Form("%s.gif", cOut->GetName()));
p=cOut->cd(3);
p->SetRightMargin(0.06);p->SetTopMargin(0.06);
xy[0]=-.5; xy[1]=-0.5; xy[2]=fgkNresYsegm[fSegmentLevel]-.5; xy[3]=2.5;
- GetGraphArray(xy, kMCcluster, 1, 1);
+ if(!GetGraphArray(xy, kMCcluster, 1, 1)) AliWarning("Missing MC cluster plot.");
p=cOut->cd(4);
p->SetRightMargin(0.16);p->SetTopMargin(0.06);
TH1 *h(NULL); char hname[100], htitle[300];
// tracklet resolution/pull in y direction
- sprintf(hname, "%s_%s_Y", GetNameId(), name);
- sprintf(htitle, "Y res for \"%s\" @ %s;tg(#phi);#Delta y [cm];%s", GetNameId(), name, fgkResYsegmName[fSegmentLevel]);
+ snprintf(hname, 100, "%s_%s_Y", GetNameId(), name);
+ snprintf(htitle, 300, "Y res for \"%s\" @ %s;tg(#phi);#Delta y [cm];%s", GetNameId(), name, fgkResYsegmName[fSegmentLevel]);
if(!(h = (TH3S*)gROOT->FindObject(hname))){
Int_t nybins=fgkNresYsegm[fSegmentLevel];
if(expand) nybins*=2;
nybins, -0.5, nybins-0.5);// segment
} else h->Reset();
arr->AddAt(h, 0);
- sprintf(hname, "%s_%s_YZpull", GetNameId(), name);
- sprintf(htitle, "YZ pull for \"%s\" @ %s;%s;#Delta y / #sigma_{y};#Delta z / #sigma_{z}", GetNameId(), name, fgkResYsegmName[fSegmentLevel]);
+ snprintf(hname, 100, "%s_%s_YZpull", GetNameId(), name);
+ snprintf(htitle, 300, "YZ pull for \"%s\" @ %s;%s;#Delta y / #sigma_{y};#Delta z / #sigma_{z}", GetNameId(), name, fgkResYsegmName[fSegmentLevel]);
if(!(h = (TH3S*)gROOT->FindObject(hname))){
h = new TH3S(hname, htitle, fgkNresYsegm[fSegmentLevel], -0.5, fgkNresYsegm[fSegmentLevel]-0.5, 100, -4.5, 4.5, 100, -4.5, 4.5);
} else h->Reset();
TH1 *h(NULL); char hname[100], htitle[300];
// tracklet resolution/pull in z direction
- sprintf(hname, "%s_%s_Z", GetNameId(), name);
- sprintf(htitle, "Z res for \"%s\" @ %s;tg(#theta);#Delta z [cm];row cross", GetNameId(), name);
+ snprintf(hname, 100, "%s_%s_Z", GetNameId(), name);
+ snprintf(htitle, 300, "Z res for \"%s\" @ %s;tg(#theta);#Delta z [cm];row cross", GetNameId(), name);
if(!(h = (TH3S*)gROOT->FindObject(hname))){
h = new TH3S(hname, htitle, 50, -1., 1., 100, -1.5, 1.5, 2, -0.5, 1.5);
} else h->Reset();
arr->AddAt(h, 2);
- sprintf(hname, "%s_%s_Zpull", GetNameId(), name);
- sprintf(htitle, "Z pull for \"%s\" @ %s;tg(#theta);#Delta z / #sigma_{z};row cross", GetNameId(), name);
+ snprintf(hname, 100, "%s_%s_Zpull", GetNameId(), name);
+ snprintf(htitle, 300, "Z pull for \"%s\" @ %s;tg(#theta);#Delta z / #sigma_{z};row cross", GetNameId(), name);
if(!(h = (TH3S*)gROOT->FindObject(hname))){
h = new TH3S(hname, htitle, 50, -1., 1., 100, -5.5, 5.5, 2, -0.5, 1.5);
h->GetZaxis()->SetBinLabel(1, "no RC");
arr->AddAt(h, 3);
// tracklet to track phi resolution
- sprintf(hname, "%s_%s_PHI", GetNameId(), name);
- sprintf(htitle, "#Phi res for \"%s\" @ %s;tg(#phi);#Delta #phi [rad];entries", GetNameId(), name);
+ snprintf(hname, 100, "%s_%s_PHI", GetNameId(), name);
+ snprintf(htitle, 300, "#Phi res for \"%s\" @ %s;tg(#phi);#Delta #phi [rad];entries", GetNameId(), name);
if(!(h = (TH2I*)gROOT->FindObject(hname))){
h = new TH2I(hname, htitle, 21, -.33, .33, 100, -.5, .5);
} else h->Reset();
TAxis *ax(NULL);
// snp pulls
- sprintf(hname, "%s_%s_SNPpull", GetNameId(), name);
- sprintf(htitle, "SNP pull for \"%s\" @ %s;tg(#phi);#Delta snp / #sigma_{snp};entries", GetNameId(), name);
+ snprintf(hname, 100, "%s_%s_SNPpull", GetNameId(), name);
+ snprintf(htitle, 300, "SNP pull for \"%s\" @ %s;tg(#phi);#Delta snp / #sigma_{snp};entries", GetNameId(), name);
if(!(h = (TH2I*)gROOT->FindObject(hname))){
h = new TH2I(hname, htitle, 60, -.3, .3, 100, -4.5, 4.5);
} else h->Reset();
arr->AddAt(h, 5);
// theta resolution
- sprintf(hname, "%s_%s_THT", GetNameId(), name);
- sprintf(htitle, "#Theta res for \"%s\" @ %s;tg(#theta);#Delta #theta [rad];entries", GetNameId(), name);
+ snprintf(hname, 100, "%s_%s_THT", GetNameId(), name);
+ snprintf(htitle, 300, "#Theta res for \"%s\" @ %s;tg(#theta);#Delta #theta [rad];entries", GetNameId(), name);
if(!(h = (TH2I*)gROOT->FindObject(hname))){
h = new TH2I(hname, htitle, 100, -1., 1., 100, -5e-3, 5e-3);
} else h->Reset();
arr->AddAt(h, 6);
// tgl pulls
- sprintf(hname, "%s_%s_TGLpull", GetNameId(), name);
- sprintf(htitle, "TGL pull for \"%s\" @ %s;tg(#theta);#Delta tgl / #sigma_{tgl};entries", GetNameId(), name);
+ snprintf(hname, 100, "%s_%s_TGLpull", GetNameId(), name);
+ snprintf(htitle, 300, "TGL pull for \"%s\" @ %s;tg(#theta);#Delta tgl / #sigma_{tgl};entries", GetNameId(), name);
if(!(h = (TH2I*)gROOT->FindObject(hname))){
h = new TH2I(hname, htitle, 100, -1., 1., 100, -4.5, 4.5);
} else h->Reset();
for(Int_t i=0; i<kNdpt+1; i++,lDPt+=2.e-3) binsDPt[i]=lDPt;
// Pt resolution
- sprintf(hname, "%s_%s_Pt", GetNameId(), name);
- sprintf(htitle, "P_{t} res for \"%s\" @ %s;p_{t} [GeV/c];#Delta p_{t}/p_{t}^{MC};SPECIES", GetNameId(), name);
+ snprintf(hname, 100, "%s_%s_Pt", GetNameId(), name);
+ snprintf(htitle, 300, "#splitline{P_{t} res for}{\"%s\" @ %s};p_{t} [GeV/c];#Delta p_{t}/p_{t}^{MC};SPECIES", GetNameId(), name);
if(!(h = (TH3S*)gROOT->FindObject(hname))){
h = new TH3S(hname, htitle,
kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
} else h->Reset();
arr->AddAt(h, 8);
// 1/Pt pulls
- sprintf(hname, "%s_%s_1Pt", GetNameId(), name);
- sprintf(htitle, "1/P_{t} pull for \"%s\" @ %s;1/p_{t}^{MC} [c/GeV];#Delta(1/p_{t})/#sigma(1/p_{t});SPECIES", GetNameId(), name);
+ snprintf(hname, 100, "%s_%s_1Pt", GetNameId(), name);
+ snprintf(htitle, 300, "#splitline{1/P_{t} pull for}{\"%s\" @ %s};1/p_{t}^{MC} [c/GeV];#Delta(1/p_{t})/#sigma(1/p_{t});SPECIES", GetNameId(), name);
if(!(h = (TH3S*)gROOT->FindObject(hname))){
h = new TH3S(hname, htitle,
kNpt, 0., 2., 100, -4., 4., kNspc, -5.5, 5.5);
} else h->Reset();
arr->AddAt(h, 9);
// P resolution
- sprintf(hname, "%s_%s_P", GetNameId(), name);
- sprintf(htitle, "P res for \"%s\" @ %s;p [GeV/c];#Delta p/p^{MC};SPECIES", GetNameId(), name);
+ snprintf(hname, 100, "%s_%s_P", GetNameId(), name);
+ snprintf(htitle, 300, "P res for \"%s\" @ %s;p [GeV/c];#Delta p/p^{MC};SPECIES", GetNameId(), name);
if(!(h = (TH3S*)gROOT->FindObject(hname))){
h = new TH3S(hname, htitle,
kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
return kTRUE;
}
+//________________________________________________________
+Bool_t AliTRDresolution::Process(TH2* const h2, TGraphErrors **g, Int_t stat)
+{
+// Robust function to process sigma/mean for 2D plot dy(x)
+// For each x bin a gauss fit is performed on the y projection and the range
+// with the minimum chi2/ndf is choosen
+
+ if(!h2) {
+ if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : NULL pointer input container.\n");
+ return kFALSE;
+ }
+ if(!Int_t(h2->GetEntries())){
+ if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : Empty h[%s - %s].\n", h2->GetName(), h2->GetTitle());
+ return kFALSE;
+ }
+ if(!g || !g[0]|| !g[1]) {
+ if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : NULL pointer output container.\n");
+ return kFALSE;
+ }
+ // prepare
+ TAxis *ax(h2->GetXaxis()), *ay(h2->GetYaxis());
+ Float_t ymin(ay->GetXmin()), ymax(ay->GetXmax()), dy(ay->GetBinWidth(1)), y0(0.), y1(0.);
+ TF1 f("f", "gaus", ymin, ymax);
+ Int_t n(0);
+ if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
+ if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
+ TH1D *h(NULL);
+ if((h=(TH1D*)gROOT->FindObject("py"))) delete h;
+ Double_t x(0.), y(0.), ex(0.), ey(0.), sy(0.), esy(0.);
+
+
+ // do actual loop
+ Float_t chi2OverNdf(0.);
+ for(Int_t ix = 1, np=0; ix<=ax->GetNbins(); ix++){
+ x = ax->GetBinCenter(ix); ex= ax->GetBinWidth(ix)*0.288; // w/sqrt(12)
+ ymin = ay->GetXmin(); ymax = ay->GetXmax();
+
+ h = h2->ProjectionY("py", ix, ix);
+ if((n=(Int_t)h->GetEntries())<stat){
+ if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("I-AliTRDresolution::Process() : Low statistics @ x[%f] stat[%d]=%d [%d].\n", x, ix, n, stat);
+ continue;
+ }
+ // looking for a first order mean value
+ f.SetParameter(1, 0.); f.SetParameter(2, 3.e-2);
+ h->Fit(&f, "QNW");
+ chi2OverNdf = f.GetChisquare()/f.GetNDF();
+ printf("x[%f] range[%f %f] chi2[%f] ndf[%d] chi2/ndf[%f]\n", x, ymin, ymax, f.GetChisquare(),f.GetNDF(),chi2OverNdf);
+ y = f.GetParameter(1); y0 = y-4*dy; y1 = y+4*dy;
+ ey = f.GetParError(1);
+ sy = f.GetParameter(2); esy = f.GetParError(2);
+// // looking for the best chi2/ndf value
+// while(ymin<y0 && ymax>y1){
+// f.SetParameter(1, y);
+// f.SetParameter(2, sy);
+// h->Fit(&f, "QNW", "", y0, y1);
+// printf(" range[%f %f] chi2[%f] ndf[%d] chi2/ndf[%f]\n", y0, y1, f.GetChisquare(),f.GetNDF(),f.GetChisquare()/f.GetNDF());
+// if(f.GetChisquare()/f.GetNDF() < Chi2OverNdf){
+// chi2OverNdf = f.GetChisquare()/f.GetNDF();
+// y = f.GetParameter(1); ey = f.GetParError(1);
+// sy = f.GetParameter(2); esy = f.GetParError(2);
+// printf(" set y[%f] sy[%f] chi2/ndf[%f]\n", y, sy, chi2OverNdf);
+// }
+// y0-=dy; y1+=dy;
+// }
+ g[0]->SetPoint(np, x, y);
+ g[0]->SetPointError(np, ex, ey);
+ g[1]->SetPoint(np, x, sy);
+ g[1]->SetPointError(np, ex, esy);
+ np++;
+ }
+ return kTRUE;
+}
+
//________________________________________________________
Bool_t AliTRDresolution::Process(TH2 * const h2, TF1 *f, Float_t k, TGraphErrors **g)
// Do the processing
//
- Char_t pn[10]; sprintf(pn, "p%03d", fIdxPlot);
+ Char_t pn[10]; snprintf(pn, 10, "p%03d", fIdxPlot);
Int_t n = 0;
if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
if(!(gm = (TGraphAsymmErrors*)((TObjArray*)fGraphM->At(plot))->At(0))) return kFALSE;
if(!(gs = (TGraphErrors*)((TObjArray*)fGraphS->At(plot)))) return kFALSE;
- Float_t x, r, mpv, xM, xm;
+ Float_t x(0.), r(0.), mpv(0.), xM(0.), xm(0.);
TAxis *ay = h3->GetYaxis();
for(Int_t iy=1; iy<=h3->GetNbinsY(); iy++){
ay->SetRange(iy, iy);
//
if(np<40){
- if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTrack: Not enough clusters to fit a track [%d].", np);
+ if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTrack: Not enough clusters to fit a track [%d].\n", np);
return kFALSE;
}
TLinearFitter yfitter(2, "pol1"), zfitter(2, "pol1");
Double_t dydx = yfitter.GetParameter(1);
param[0] = x0; param[1] = y0; param[2] = z0; param[3] = dydx; param[4] = dzdx;
- if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>3) printf("D-AliTRDresolution::FitTrack: x0[%f] y0[%f] z0[%f] dydx[%f] dzdx[%f]\n", x0, y0, z0, dydx, dzdx);
+ if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>3) printf("D-AliTRDresolution::FitTrack: x0[%f] y0[%f] z0[%f] dydx[%f] dzdx[%f].\n", x0, y0, z0, dydx, dzdx);
return kTRUE;
}
//____________________________________________________________________
-Bool_t AliTRDresolution::FitTracklet(const Int_t ly, const Int_t np, AliTrackPoint *points, const Float_t param[10], Float_t par[3])
+Bool_t AliTRDresolution::FitTracklet(const Int_t ly, const Int_t np, const AliTrackPoint *points, const Float_t param[10], Float_t par[3])
{
//
// Fit tracklet with a staight line using the coresponding subset of clusters out of the total "np" clusters stored in the array "points".
nly++;
}
if(nly<10){
- if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTracklet: Not enough clusters to fit a tracklet [%d].", nly);
+ if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTracklet: Not enough clusters to fit a tracklet [%d].\n", nly);
return kFALSE;
}
// set radial reference for fit
}
//____________________________________________________________________
-Bool_t AliTRDresolution::UseTrack(const Int_t np, AliTrackPoint *points, Float_t param[10])
+Bool_t AliTRDresolution::UseTrack(const Int_t np, const AliTrackPoint *points, Float_t param[10])
{
//
// Global selection mechanism of tracksbased on cluster to fit residuals
};
memcpy(fAxTitle, lAxTitle, 4*kNprojs*sizeof(Char_t*));
}
+
+#include "TFile.h"
+//________________________________________________________
+Bool_t AliTRDresolution::LoadCorrection(const Char_t *file)
+{
+ if(!file){
+ AliWarning("Use cluster position as in reconstruction.");
+ SetLoadCorrection();
+ return kTRUE;
+ }
+ TDirectory *cwd(gDirectory);
+ TString fileList;
+ FILE *filePtr = fopen(file, "rt");
+ if(!filePtr){
+ AliError(Form("Couldn't open correction list \"%s\". Use cluster position as in reconstruction.", file));
+ SetLoadCorrection();
+ return kFALSE;
+ }
+ TH2 *h2 = new TH2F("h2", ";time [time bins];detector;dx [#mum]", 30, -0.5, 29.5, AliTRDgeometry::kNdet, -0.5, AliTRDgeometry::kNdet-0.5);
+ while(fileList.Gets(filePtr)){
+ if(!TFile::Open(fileList.Data())) {
+ AliWarning(Form("Couldn't open \"%s\"", fileList.Data()));
+ continue;
+ } else AliInfo(Form("\"%s\"", fileList.Data()));
+
+ TTree *tSys = (TTree*)gFile->Get("tSys");
+ h2->SetDirectory(gDirectory); h2->Reset("ICE");
+ tSys->Draw("det:t>>h2", "dx", "goff");
+ for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
+ for(Int_t it(0); it<30; it++) fXcorr[idet][it]+=(1.e-4*h2->GetBinContent(it+1, idet+1));
+ }
+ h2->SetDirectory(cwd);
+ gFile->Close();
+ }
+ cwd->cd();
+
+ if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>=2){
+ for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
+ printf("DET|");for(Int_t it(0); it<30; it++) printf(" tb%02d|", it); printf("\n");
+ for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
+ for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
+ printf("%03d|", idet);
+ for(Int_t it(0); it<30; it++) printf("%+5.0f|", 1.e4*fXcorr[idet][it]);
+ printf("\n");
+ }
+ }
+ SetLoadCorrection();
+ return kTRUE;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff