,fIdxPlot(0)
,fIdxFrame(0)
,fPtThreshold(1.)
+ ,fDyRange(0.75)
,fDBPDG(NULL)
,fGraphS(NULL)
,fGraphM(NULL)
//
SetNameTitle("TRDresolution", "TRD spatial and momentum resolution");
SetSegmentationLevel();
+ memset(fXcorr, 0, AliTRDgeometry::kNdet*30*sizeof(Float_t));
}
//________________________________________________________
,fIdxPlot(0)
,fIdxFrame(0)
,fPtThreshold(1.)
+ ,fDyRange(0.75)
,fDBPDG(NULL)
,fGraphS(NULL)
,fGraphM(NULL)
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);
}
// Uses functionality defined by AliTRDseedV1.
Double_t t2(tilt*tilt);
+ // exit door until a bug fix is found for AliTRDseedV1::GetCovSqrt
+ return kTRUE;
// rotate along pad
Double_t cc[3];
Int_t sgm[3];
Double_t covR[7], cov[3], dy[2], dz[2];
- Float_t pt, x0, y0, z0, dydx, dzdx, tilt;
+ Float_t pt, x0, y0, z0, dydx, dzdx, tilt(0.);
const AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
AliTRDseedV1 *fTracklet(NULL); TObjArray *clInfoArr(NULL);
// do LINEAR track refit if asked by the user
// it is the user responsibility to check if B=0T
Float_t param[10]; memset(param, 0, 10*sizeof(Float_t));
+ Int_t np(0), nrc(0); AliTrackPoint clusters[300];
if(HasTrackRefit()){
- Int_t np(0); AliTrackPoint clusters[300];
+ Bool_t kPrimary(kFALSE);
for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
if(!fTracklet->IsOK()) continue;
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);
clusters[np].SetXYZ(xyz);
np++;
}
if(fTracklet->GetEstimatedCrossPoint(xcross, zcross)){
clusters[np].SetCharge(tilt);
clusters[np].SetClusterType(1);
+ clusters[np].SetVolumeID(ily);
clusters[np].SetXYZ(xcross, 0., zcross);
np++;
+ nrc++;
}
}
- if(fTracklet->IsPrimary()){
- clusters[np].SetCharge(tilt);
- clusters[np].SetClusterType(1);
- clusters[np].SetXYZ(0., 0., 0.);
- np++;
+ if(fTracklet->IsPrimary()) kPrimary = kTRUE;
+ }
+ if(kPrimary){
+ clusters[np].SetCharge(tilt);
+ clusters[np].SetClusterType(1);
+ clusters[np].SetVolumeID(-1);
+ clusters[np].SetXYZ(0., 0., 0.);
+ np++;
+ }
+ if(!FitTrack(np, clusters, param)) {
+ AliDebug(1, "Linear track Fit failed.");
+ return NULL;
+ }
+ if(HasTrackSelection()){
+ Bool_t kReject(kFALSE);
+ if(fkTrack->GetNumberOfTracklets() != AliTRDgeometry::kNlayer) kReject = kTRUE;
+ if(!kReject && !UseTrack(np, clusters, param)) kReject = kTRUE;
+ if(kReject){
+ AliDebug(1, "Reject track for residuals analysis.");
+ return NULL;
}
}
- FitTrack(np, clusters, param);
}
for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
// retrive the track angle with the chamber
if(HasTrackRefit()){
- dydx = param[3];
+ Float_t par[3];
+ if(!FitTracklet(ily, np, clusters, param, par)) continue;
+ dydx = par[2];//param[3];
dzdx = param[4];
- y0 = param[1] + dydx * (x0 - param[0]);
+ y0 = par[1] + dydx * (x0 - par[0]);//param[1] + dydx * (x0 - param[0]);
z0 = param[2] + dzdx * (x0 - param[0]);
} else {
- y0 = fTracklet->GetYref(0);
+ y0 = fTracklet->GetYfit(0);//fTracklet->GetYref(0);
z0 = fTracklet->GetZref(0);
- dydx = fTracklet->GetYref(1);
+ dydx = fTracklet->GetYfit(1);//fTracklet->GetYref(1);
dzdx = fTracklet->GetZref(1);
}
+ /*printf("RC[%c] Primary[%c]\n"
+ " Fit : y0[%f] z0[%f] dydx[%f] dzdx[%f]\n"
+ " Ref: y0[%f] z0[%f] dydx[%f] dzdx[%f]\n", fTracklet->IsRowCross()?'y':'n', fTracklet->IsPrimary()?'y':'n', y0, z0, dydx, dzdx
+ ,fTracklet->GetYref(0),fTracklet->GetZref(0),fTracklet->GetYref(1),fTracklet->GetZref(1));*/
tilt = fTracklet->GetTilt();
fTracklet->GetCovRef(covR);
Double_t t2(tilt*tilt)
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;
// rotate along pad
dy[1] = cost*(dy[0] - dz[0]*tilt);
dz[1] = cost*(dz[0] + dy[0]*tilt);
- if(pt>fPtThreshold && c->IsInChamber()) ((TH3S*)arr->At(0))->Fill(dydx, dy[1], sgm[fSegmentLevel]);
+ if(pt>fPtThreshold && c->IsInChamber()) ((TH3S*)arr->At(0))->Fill(fTracklet->GetYref(1), dy[0], sgm[fSegmentLevel]);
// tilt rotation of covariance for clusters
Double_t sy2(c->GetSigmaY2()), sz2(c->GetSigmaZ2());
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;
}
}
if(clInfoArr) delete clInfoArr;
+
return (TH3S*)arr->At(0);
}
dy[1]= cost*(dy[0] - dz[0]*tilt);
dz[1]= cost*(dz[0] + dy[0]*tilt);
((TH3S*)arr->At(0))->Fill(phi, dy[1], sgm[fSegmentLevel]+rc*fgkNresYsegm[fSegmentLevel]);
- ((TH3S*)arr->At(2))->Fill(tht, dz[1], rc);
+ if(rc) ((TH2S*)arr->At(2))->Fill(tht, dz[1]);
// compute covariance matrix
fTracklet->GetCovAt(x, cov);
((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dyz[0], dyz[1]);
((TH3S*)arr->At(3))->Fill(tht, dyz[1], rc);
- Double_t dphi((phi-fTracklet->GetYfit(1))/(1-phi*fTracklet->GetYfit(1)));
+ // calculate angular residuals and correct for tilt
+ Double_t phiTrklt = fTracklet->GetYfit(1);
+ phiTrklt += tilt*tht;
+ Double_t dphi((phi-phiTrklt)/(1-phi*phiTrklt));
Double_t dtht((tht-fTracklet->GetZfit(1))/(1-tht*fTracklet->GetZfit(1)));
- ((TH2I*)arr->At(4))->Fill(phi, TMath::ATan(dphi));
+ ((TH3S*)arr->At(4))->Fill(phi, TMath::ATan(dphi), sgm[fSegmentLevel]);
if(DebugLevel()>=1){
UChar_t err(fTracklet->GetErrorMsg());
<< "\n";
}
}
-
-
return (TH2I*)arr->At(0);
}
break;
}
if(!fTracklet || TMath::Abs(x-fTracklet->GetX())>1.e-3){
- AliWarning("Tracklet did not match Track.");
+ AliDebug(1, "Tracklet did not match Track.");
return NULL;
}
Int_t sgm[3];
// PID calculation.
if(track) fkTrack = track;
+ return NULL;
if(!fkTrack){
AliDebug(4, "No Track defined.");
return NULL;
break;
}
if(!fTracklet || TMath::Abs(x-fTracklet->GetX())>1.e-3){
- AliWarning("Tracklet did not match Track position.");
+ AliDebug(1, "Tracklet did not match Track position.");
return NULL;
}
Int_t sgm[3];
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()));
- if(!HasMCdata()){
+ if(!HasMCdata() ||
+ (!fGraphS->At(kMCcluster) || !fGraphM->At(kMCcluster) ||
+ !fGraphS->At(kMCtracklet) || !fGraphM->At(kMCtracklet))){
delete cOut;
return;
}
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);
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, kMCtracklet, 1, 1);
+ if(!GetGraphArray(xy, kMCtracklet, 1, 1)){
+ AliWarning("Failed retrieve tracklet resolution plot.");
+ }
cOut->SaveAs(Form("%s.gif", cOut->GetName()));
delete cOut;
}
//________________________________________________________
-TObjArray* AliTRDresolution::BuildMonitorContainerCluster(const char* name, Bool_t expand)
+TObjArray* AliTRDresolution::BuildMonitorContainerCluster(const char* name, Bool_t expand, Float_t range)
{
// Build performance histograms for AliTRDcluster.vs TRD track or MC
// - y reziduals/pulls
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]);
+ Float_t rr = range<0.?fDyRange:range;
if(!(h = (TH3S*)gROOT->FindObject(hname))){
Int_t nybins=fgkNresYsegm[fSegmentLevel];
if(expand) nybins*=2;
h = new TH3S(hname, htitle,
- 48, -.48, .48, 60, -1.5, 1.5, nybins, -0.5, nybins-0.5);
+ 48, -.48, .48, // phi
+ 60, -rr, rr, // dy
+ 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();
// - y reziduals/pulls
// - z reziduals/pulls
// - phi reziduals
- TObjArray *arr = BuildMonitorContainerCluster(name, expand);
+ TObjArray *arr = BuildMonitorContainerCluster(name, expand, 0.05);
arr->Expand(5);
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);
- if(!(h = (TH3S*)gROOT->FindObject(hname))){
- h = new TH3S(hname, htitle, 50, -1., 1., 100, -1.5, 1.5, 2, -0.5, 1.5);
+ snprintf(hname, 100, "%s_%s_Z", GetNameId(), name);
+ snprintf(htitle, 300, "Z res for \"%s\" @ %s;tg(#theta);#Delta z [cm]", GetNameId(), name);
+ if(!(h = (TH2S*)gROOT->FindObject(hname))){
+ h = new TH2S(hname, htitle, 50, -1., 1., 100, -.05, .05);
} 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);
- if(!(h = (TH2I*)gROOT->FindObject(hname))){
- h = new TH2I(hname, htitle, 21, -.33, .33, 100, -.5, .5);
+ snprintf(hname, 100, "%s_%s_PHI", GetNameId(), name);
+ snprintf(htitle, 300, "#Phi res for \"%s\" @ %s;tg(#phi);#Delta #phi [rad];%s", GetNameId(), name, fgkResYsegmName[fSegmentLevel]);
+ Int_t nsgms=fgkNresYsegm[fSegmentLevel];
+ if(!(h = (TH3S*)gROOT->FindObject(hname))){
+ h = new TH3S(hname, htitle, 48, -.48, .48, 100, -.5, .5, nsgms, -0.5, nsgms-0.5);
} else h->Reset();
arr->AddAt(h, 4);
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);
if(fContainer) return fContainer;
fContainer = new TObjArray(kNviews);
- //fContainer->SetOwner(kTRUE);
+ fContainer->SetOwner(kTRUE);
TH1 *h(NULL);
TObjArray *arr(NULL);
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");
for(Int_t ip(0); ip<np; ip++) x0+=points[ip].GetX();
x0/=Float_t(np);
- Double_t x, y, z, dx, tilt;
+ Double_t x, y, z, dx, tilt(0.);
for(Int_t ip(0); ip<np; ip++){
x = points[ip].GetX(); z = points[ip].GetZ();
dx = x - x0;
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, 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".
+// See function FitTrack for the data stored in the "clusters" array
+
+// The parameters of the straight line fit are stored in the array "param" in the following order :
+// par[0] - x0 reference radial position
+// par[1] - y0 reference r-phi position @ x0
+// par[2] - slope dy/dx
+//
+// Attention :
+// Function should be used to refit tracks for B=0T
+//
+
+ TLinearFitter yfitter(2, "pol1");
+
+ // grep data for tracklet
+ Double_t x0(0.), x[60], y[60], dy[60];
+ Int_t nly(0);
+ for(Int_t ip(0); ip<np; ip++){
+ if(points[ip].GetClusterType()) continue;
+ if(points[ip].GetVolumeID() != ly) continue;
+ Float_t xt(points[ip].GetX())
+ ,yt(param[1] + param[3] * (xt - param[0]));
+ x[nly] = xt;
+ y[nly] = points[ip].GetY();
+ dy[nly]= y[nly]-yt;
+ x0 += xt;
+ nly++;
+ }
+ if(nly<10){
+ 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
+ x0 /= Float_t(nly);
+
+ // find tracklet core
+ Double_t mean(0.), sig(1.e3);
+ AliMathBase::EvaluateUni(nly, dy, mean, sig, 0);
+
+ // simple cluster error parameterization
+ Float_t kSigCut = TMath::Sqrt(5.e-4 + param[3]*param[3]*0.018);
+
+ // fit tracklet core
+ for(Int_t jly(0); jly<nly; jly++){
+ if(TMath::Abs(dy[jly]-mean)>kSigCut) continue;
+ Double_t dx(x[jly]-x0);
+ yfitter.AddPoint(&dx, y[jly], 1.);
+ }
+ if(yfitter.Eval() != 0) return kFALSE;
+ par[0] = x0;
+ par[1] = yfitter.GetParameter(0);
+ par[2] = yfitter.GetParameter(1);
+ return kTRUE;
+}
+
+//____________________________________________________________________
+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
+// The parameters are the same as used ni function FitTrack().
+
+ const Float_t kS(0.6), kM(0.2);
+ TH1S h("h1", "", 100, -5.*kS, 5.*kS);
+ Float_t dy, dz, s, m;
+ for(Int_t ip(0); ip<np; ip++){
+ if(points[ip].GetClusterType()) continue;
+ Float_t x0(points[ip].GetX())
+ ,y0(param[1] + param[3] * (x0 - param[0]))
+ ,z0(param[2] + param[4] * (x0 - param[0]));
+ dy=points[ip].GetY() - y0; h.Fill(dy);
+ dz=points[ip].GetZ() - z0;
+ }
+ TF1 fg("fg", "gaus", -5.*kS, 5.*kS);
+ fg.SetParameter(1, 0.);
+ fg.SetParameter(2, 2.e-2);
+ h.Fit(&fg, "QN");
+ m=fg.GetParameter(1); s=fg.GetParameter(2);
+ if(s>kS || TMath::Abs(m)>kM) return kFALSE;
return kTRUE;
}
};
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