// we suppose the chamber is well calibrated for t_{0} and aligned in
// radial direction.
//
-// Clusters can be radially shifted due to two causes:
+// Clusters can be radially shifted due to three causes:
// - globally shifted - due to residual misalignment/miscalibration(t0)
// - locally shifted - due to different local drift velocity from the mean
// - randomly shifted - due to neighboring (radial direction) clusters
// END_LATEX
//
//
+// Clusters can also be r-phi shifted due to:
+// - wrong PRF or wrong cuts at digits level
+//The following correction is applied :
+// BEGIN_LATEX
+// <#Delta y> = a + b * sin(c*y_{pw})
+// END_LATEX
+
// The Models
//
// Parameterization against total charge
//
// The parameterization of the error in the y direction along track uses
// BEGIN_LATEX
-// #sigma_{y}^{||} = #sigma_{y}^{0} + exp(-a*(x-b))
+// #sigma_{y}^{||} = #sigma_{y}^{0} -a*exp(1/(x-b))
// END_LATEX
//
// with following values for the parameters:
//
// // initialize DB manager
// AliCDBManager *cdb = AliCDBManager::Instance();
-// cdb->SetDefaultStorage("local://$ALICE_ROOT");
+// cdb->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
// cdb->SetRun(0);
// // initialize magnetic field.
-// AliMagFMaps *field = new AliMagFMaps("Maps","Maps", 2, 1., 10., AliMagFMaps::k5kG);
+// AliMagFCheb *field=new AliMagFCheb("Maps","Maps", 2, 1., 10., AliMagFCheb::k5kG);
// AliTracker::SetFieldMap(field, kTRUE);
//
// AliTRDclusterResolution *res = new AliTRDclusterResolution();
// res->SetMCdata();
// res->Load("TRD.TaskClErrParam.root");
-// res->SetExB();
+// res->SetExB();
+// res->SetVisual();
+// //res->SetSaveAs();
+// res->SetProcessCharge(kFALSE);
+// res->SetProcessCenterPad(kFALSE);
+// //res->SetProcessMean(kFALSE);
+// res->SetProcessSigma(kFALSE);
// if(!res->PostProcess()) return;
// new TCanvas;
// res->GetRefFigure(fig);
////////////////////////////////////////////////////////////////////////////
#include "AliTRDclusterResolution.h"
-#include "AliTRDtrackInfo/AliTRDclusterInfo.h"
+#include "info/AliTRDclusterInfo.h"
#include "AliTRDgeometry.h"
#include "AliTRDcalibDB.h"
+#include "AliTRDCommonParam.h"
#include "Cal/AliTRDCalROC.h"
#include "Cal/AliTRDCalDet.h"
#include "AliTracker.h"
#include "AliCDBManager.h"
+#include "TROOT.h"
#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()
- : AliTRDrecoTask("ClErrParam", "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.)
- ,fDet(-1)
+ ////////////
+ ,fLy(0)
+ ,fX(0.)
+ ,fY(0.)
+ ,fZ(0.)
{
- // ideal equidistant binning
- //fAt = new TAxis(kNTB, -0.075, (kNTB-.5)*.15);
-
- // equidistant binning for scaled for X0-MC (track ref)
- fAt = new TAxis(kNTB, -0.075+.088, (kNTB-.5)*.15+.088);
-
- // non-equidistant binning after vdrift correction
-// const Double_t x[kNTB+1] = {
-// 0.000000, 0.185084, 0.400490, 0.519701, 0.554653, 0.653150, 0.805063, 0.990261,
-// 1.179610, 1.356406, 1.524094, 1.685499, 1.843083, 1.997338, 2.148077, 2.298274,
-// 2.448656, 2.598639, 2.747809, 2.896596, 3.045380, 3.195000, 3.340303, 3.461688,
-// 3.540094, 3.702677};
-// fAt = new TAxis(kNTB, x);
-
- fAd = new TAxis(kND, 0., .25);
+ memset(fR, 0, 4*sizeof(Float_t));
+ memset(fP, 0, 4*sizeof(Float_t));
+ // time drift axis
+ fAt = new TAxis(kNTB, 0., kNTB*fgkTimeBinLength);
+
+ // By default register all analysis
+ // The user can switch them off in his steering macro
+ SetProcess(kQRes);
+ SetProcess(kCenter);
+ SetProcess(kMean);
+ SetProcess(kSigm);
}
//_______________________________________________________
AliTRDclusterResolution::~AliTRDclusterResolution()
{
- if(fAd) delete fAd;
+ if(fCanvas) delete fCanvas;
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;
+ TH2 *h2 = 0x0;TH1 *h1 = 0x0;
TGraphErrors *gm(0x0), *gs(0x0), *gp(0x0);
switch(ifig){
case kQRes:
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;
- if(!(h2 = (TH2F*)arr->At(0))) break;
- h2->Draw("lego2");
+ 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;
- gPad->Divide(2, 1);
- gPad->cd(1);
- if(!(h2 = (TH2F*)arr->At(0))) break;
- h2->Draw("lego2fb");
- gPad->cd(2);
- if(!(h2 = (TH2F*)arr->At(1))) break;
- h2->Draw("lego2fb");
+ gPad->Divide(2, 1); l = gPad->GetListOfPrimitives();
+ if(!(h2 = (TH2F*)arr->At(0))) return kFALSE;
+ ((TVirtualPad*)l->At(0))->cd();
+ h1 = h2->ProjectionY("hsx_pyy"); h1->Scale(1.e4/kND); h1->SetMarkerStyle(24);
+ h1->SetYTitle("<#sigma_{x}> [#mum]");
+ h1->GetXaxis()->SetRange(2, kNTB-1); h1->Draw("pc");
+
+ if(!(h2 = (TH2F*)arr->At(1))) return kFALSE;
+ ((TVirtualPad*)l->At(1))->cd();
+ h1 = h2->ProjectionY("hsy_pyy"); h1->Scale(1.e4/kND); h1->SetMarkerStyle(24);
+ h1->SetYTitle("<#sigma_{y}> [#mum]");
+ h1->GetXaxis()->SetRange(2, kNTB-1); h1->Draw("pc");
return kTRUE;
case kMean:
if(!(arr = (TObjArray*)fResults->At(kMean))) break;
- arr->ls();
-/* gPad->Divide(2, 1);
- gPad->cd(1);*/
- if(!(h2 = (TH2F*)arr->At(0))) break;
- h2->Draw("lego2 fb");
-/* gPad->cd(2);
- if(!(h2 = (TH2F*)arr->At(1))) break;
- h2->Draw("lego2fb");*/
+ gPad->Divide(2, 1); l = gPad->GetListOfPrimitives();
+ ((TVirtualPad*)l->At(0))->cd();
+ if(!(gm = (TGraphErrors*)arr->At(0))) return kFALSE;
+ gm->Draw("apl");
+ gm->GetHistogram()->SetXTitle("t_{drift} [#mus]");
+ gm->GetHistogram()->SetYTitle("dx [#mum]");
+
+ ((TVirtualPad*)l->At(1))->cd();
+ if(!(gm = (TGraphErrors*)arr->At(1))) return kFALSE;
+ gm->Draw("apl");
+ gm->GetHistogram()->SetXTitle("t_{drift} [#mus]");
+ gm->GetHistogram()->SetYTitle("dy [#mum]");
+
return kTRUE;
default:
break;
TObjArray* AliTRDclusterResolution::Histos()
{
if(fContainer) return fContainer;
- fContainer = new TObjArray(sizeof(AliTRDclusterResolution::EResultContainers));
+ 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);
- 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);
- 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("h_d%02dt%02d", id, it), Form("d_{wire}(%3.1f-%3.1f)[mm] x_{drift}(%5.2f-%5.2f)[mm]", 10.*fAd->GetBinCenter(id)- 5.*fAd->GetBinWidth(id), 10.*fAd->GetBinCenter(id)+ 5.*fAd->GetBinWidth(id), 10.*fAt->GetBinCenter(it)- 5.*fAt->GetBinWidth(it), 10.*fAt->GetBinCenter(it)+ 5.*fAt->GetBinWidth(it)), 35, -.35, .35, 100, -.5, .5), ih++);
- h2->SetXTitle("tg#phi");
- h2->SetYTitle("#Delta y[cm]");
- h2->SetZTitle("entries");
+ 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");
+ 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);
- 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("h_d%02dt%02d", id, it), Form("d_{wire}(%3.1f-%3.1f)[mm] x_{drift}(%5.2f-%5.2f)[mm]", 10.*fAd->GetBinCenter(id)- 5.*fAd->GetBinWidth(id), 10.*fAd->GetBinCenter(id)+ 5.*fAd->GetBinWidth(id), 10.*fAt->GetBinCenter(it)- 5.*fAt->GetBinWidth(it), 10.*fAt->GetBinCenter(it)+ 5.*fAt->GetBinWidth(it)), 35, -.35, .35, 100, -.5, .5), ih++);
- h2->SetXTitle("tg#phi-h*tg(#theta)");
- h2->SetYTitle("#Delta y[cm]");
- h2->SetZTitle("entries");
+ fContainer->AddAt(arr = new TObjArray(kNTB), kMean);
+ arr->SetName("Systematics");
+ 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(cli->GetDriftLength());
- if(it==0 || it == fAt->GetNbins()+1){
- AliWarning(Form("Drift length %f outside allowed range", cli->GetDriftLength()));
- 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::EResultContainers));
+ 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(2), kCenter);
+ // pad center dependence
+ fResults->AddAt(arr = new TObjArray(AliTRDgeometry::kNlayer+1), kCenter);
arr->SetOwner();
- arr->AddAt(h2 = new TH2F("hYM", "",
- AliTRDgeometry::kNlayer, -.5, AliTRDgeometry::kNlayer-.5, 51, -.51, .51), 0);
- h2->SetXTitle("ly");
- h2->SetYTitle("y [w]");
- h2->SetZTitle("#mu_{x} [mm]");
- arr->AddAt(h2 = (TH2F*)h2->Clone("hYS"), 1);
- h2->SetZTitle("#sigma_{x} [mm]");
-
- fResults->AddAt(arr = new TObjArray(2), kSigm);
- arr->SetOwner();
- arr->AddAt(h2 = new TH2F("hSX", "",
- fAd->GetNbins(), fAd->GetXmin(), fAd->GetXmax(),
- fAt->GetNbins(), fAt->GetXmin(), fAt->GetXmax()), 0);
- h2->SetXTitle("d [mm]");
- h2->SetYTitle("x [mm]");
- h2->SetZTitle("#sigma_{x} [mm]");
- arr->AddAt(h2 = (TH2F*)h2->Clone("hSY"), 1);
- h2->SetZTitle("#sigma_{y} [mm]");
-
- fResults->AddAt(arr = new TObjArray(2), kMean);
- arr->SetOwner();
- arr->AddAt(h2 = (TH2F*)h2->Clone("hDX"), 0);
- h2->SetZTitle("dx [mm]");
- arr->AddAt(h2 = (TH2F*)h2->Clone("hX0"), 1);
- h2->SetZTitle("x0 [mm]");
+ 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);
+ }
+
+
+ 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();
printf("ExB[%e] ExB2[%e]\n", fExB, exb2);
// process resolution dependency on charge
- // ProcessCharge();
+ if(HasProcess(kQRes)) ProcessCharge();
// process resolution dependency on y displacement
- // ProcessCenterPad();
+ if(HasProcess(kCenter)) ProcessCenterPad();
// process resolution dependency on drift legth and drift cell width
- ProcessSigma();
+ if(HasProcess(kSigm)) ProcessSigma();
// process systematic shift on drift legth and drift cell width
- //ProcessMean();
+ if(HasProcess(kMean)) ProcessMean();
return kTRUE;
}
const AliTRDCalDet *fCalVdriftDet = fCalibration->GetVdriftDet();
fVdrift = fCalVdriftDet->GetValue(det) * fCalVdriftROC->GetValue(col, row);
- fExB = fCalibration->GetOmegaTau(fVdrift, -0.1*AliTracker::GetBz());
+ fExB = AliTRDCommonParam::Instance()->GetOmegaTau(fVdrift);
SetBit(kExB);
return kTRUE;
}
+//_______________________________________________________
+void AliTRDclusterResolution::SetVisual()
+{
+ if(fCanvas) return;
+ fCanvas = new TCanvas("clResCanvas", "Cluster Resolution Visualization", 10, 10, 600, 600);
+}
+
//_______________________________________________________
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;
}
Double_t q, n = 0., entries;
ax = h2->GetXaxis();
for(Int_t ix=1; ix<=ax->GetNbins(); ix++){
- Double_t q = TMath::Exp(ax->GetBinCenter(ix));
+ q = TMath::Exp(ax->GetBinCenter(ix));
if(q<20. || q>250.) continue; // ?!
h1 = h2->ProjectionY("py", ix, ix);
// 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;
- TH2I *h2 = 0x0;
- TObjArray *arrg = (TObjArray*)fResults->At(kCenter);
- TH2F *hym = (TH2F*)arrg->At(0);
- TH2F *hys = (TH2F*)arrg->At(1);
- 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);
- if(h1->GetEntries() < 50) continue;
- Adjust(&f, h1);
- h1->Fit(&f, "Q");
-
- // Fill sy = f(y_w)
- hym->Fill(ily, yd, f.GetParameter(1));
- //gym->SetPointError(ip, 0., 10.*f.GetParError(1));
- hys->Fill(ily, yd, f.GetParameter(2));
- //gys->SetPointError(ip, 0., 2.e2*f.GetParameter(2)*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);
+
+ 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);
+
+ //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();
+
+
+ }
+ }
+ 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)*/);
+ }
+ printf("\b},\n");
+ s[il]/=n;
+
+ f.ReleaseParameter(2);
+
+
+ if(!fCanvas) continue;
+ 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");
+
+ // init visualization
+ TGraphErrors *ggs = 0x0;
+ TGraph *line = 0x0;
+ if(fCanvas){
+ ggs = new TGraphErrors();
+ line = new TGraph();
+ line->SetLineColor(kRed);line->SetLineWidth(2);
+ }
+
+ // 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;
- TH1D *h1 = 0x0;
- TH2I *h2 = 0x0;
+ Double_t exb2 = fExB*fExB;
- Double_t d(0.), x(0.), sx, sy;
- 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 = %f[mm]\n", d);
- for(Int_t it=1; it<=fAt->GetNbins(); it++){
- x = fAt->GetBinCenter(it); //[mm]
- printf(" Doing xd = %f[mm]\n", x);
- 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;
+ 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);
- Double_t tgg = (dydx-fExB)/(1.+dydx*fExB);
+
+ 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 = (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, 1.e2*f.GetParameter(2)*f.GetParameter(2), 2.e2*f.GetParameter(2)*f.GetParError(2));
+ gs.AddPoint(&tgg2, s2, s2e);
+
+ if(!ggs) continue;
+ Int_t jp = ggs->GetN();
+ ggs->SetPoint(jp, tgg2, s2);
+ ggs->SetPointError(jp, 0., s2e);
}
if(gs.Eval()) continue;
- sx = TMath::Sqrt(gs.GetParameter(1));
- sy = TMath::Sqrt(gs.GetParameter(0));
- printf("sx[%f] sy[%f] [mm]\n", sx, sy);
- hsx->SetBinContent(id, it, sx);
- //hsx->SetBinError(id, it, sig.GetParError(1));
+ // s^2_x = s0^2_x - x^2*tg^2(a_L)/12
+ 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
+ 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();
+ if(!hFrame){
+ hFrame=new TH1I("hFrame", "", 100, 0., .3);
+ hFrame->SetMinimum(0.);hFrame->SetMaximum(.005);
+ hFrame->SetXTitle("tg^{2}(#phi-#alpha_{L})");
+ hFrame->SetYTitle("#sigma^{2}y[cm^{2}]");
+ hFrame->SetLineColor(1);hFrame->SetLineWidth(1);
+ hFrame->Draw();
+ } else hFrame->Reset();
+ Double_t xx = 0., dxx=.2/50;
+ for(Int_t ip=0;ip<50;ip++){
+ line->SetPoint(ip, xx, gs.GetParameter(0)+xx*gs.GetParameter(1));
+ xx+=dxx;
+ }
+ ggs->Draw("pl"); line->Draw("l");
+ fCanvas->Modified(); fCanvas->Update();
+ if(IsSaveAs()) fCanvas->SaveAs(Form("Figures/ProcessSigma_z[%5.3f]_x[%5.3f].gif", fZ, fX));
+ else gSystem->Sleep(100);
- // s^2_y = s0^2_y + tg^2(a_L) * s^2_x
- hsy->SetBinContent(id, it, sy/* - exb2*sig.GetParameter(1)*/);
- //hsy->SetBinError(id, it, sig.GetParError(0)+exb2*exb2*sig.GetParError(1));
+ printf(" xd=%4.1f[cm] sx=%5.3e[cm] sy=%5.3e[cm]\n", fX, TMath::Sqrt(fR[0]), TMath::Sqrt(fR[1]));
}
}
-
- // fit sy parallel to the drift
- h1 = hsy->ProjectionY("pyy"); h1->Scale(1./hsy->GetNbinsX());
- TF1 fsyD("fsy", "[0]+exp([1]*(x+[2]))", 0.5, 3.5);
- h1->Fit(&fsyD);
-
-
- // fit sx parallel to the drift
- h1 = hsx->ProjectionY("pyy"); h1->Scale(1./hsx->GetNbinsX());
- TF1 fsxD("fsx", "gaus(0)+expo(3)", 0., 3.5);
- h1->Fit(&fsxD);
-
- // fit sx perpendicular to the drift
- Int_t nb = 0;
- TAxis *ay = hsx->GetYaxis();
- TH1D *hx = hsx->ProjectionX("px", 1,1);
- TH1D *hAn = (TH1D*)hx->Clone("hAn"); hAn->Reset();
- for(Int_t ib = 11; ib<24; ib++, nb++){
- hx = hsx->ProjectionX("px", ib, ib);
- for(Int_t id=1; id<=hx->GetNbinsX(); id++)
- hx->SetBinContent(id, hx->GetBinContent(id)-fsxD.Eval(ay->GetBinCenter(ib)));
- hAn->Add(hx);
- }
- hAn->Scale(1./nb);
- hAn->Fit("pol2");
+ 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", Form("%6.4f*[0]+[1]", fExB), -.12, .12);
- Double_t d(0.), x(0.), dx, x0;
+ TF1 line("l", "[0]+[1]*x", -.15, .15);
+ TGraphErrors *gm = new TGraphErrors();
+ TH1 *hFrame=0x0;
+ TH1D *h1 = 0x0; TH3S *h3 =0x0;
TAxis *ax = 0x0;
- TH1D *h1 = 0x0;
- TH2I *h2 = 0x0;
- TCanvas *c=new TCanvas;
-
- TObjArray *arrr = (TObjArray*)fResults->At(kMean);
- TH2F *hdx = (TH2F*)arrr->At(0);
- TH2F *hx0 = (TH2F*)arrr->At(1);
- for(Int_t id=1; id<=fAd->GetNbins(); id++){
- d = fAd->GetBinCenter(id); //[mm]
- printf(" Doing d = %f[mm]\n", d);
- for(Int_t it=1; it<=fAt->GetNbins(); it++){
- x = fAt->GetBinCenter(it); //[mm]
- printf(" Doing xd = %f[mm]\n", x);
- 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;
- }
- Int_t ip = 0;
- 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() < 50) 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);
+
+ 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);
- //printf("\tFit ix[%d] on %s entries [%d]\n", ix, h2->GetName(), (Int_t)h1->GetEntries());
h1->Fit(&f, "QN");
// Fill <Dy> = f(dydx - h*dzdx)
- gm->SetPoint(ip, dydx, 10.*f.GetParameter(1));
- gm->SetPointError(ip, 0., 10.*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");
- c->cd();
- gm->Draw("apl");
- c->Modified(); c->Update(); gSystem->Sleep(500);
- dx = line.GetParameter(1); // = (fVdrift + dVd)*(fT + tCorr);
- x0 = line.GetParameter(0); // = tg(a_L)*dx
- printf(" dx[%f] x0[%f] [mm]\n", dx, x0);
-
- hdx->SetBinContent(id, it, dx);
- hx0->SetBinContent(id, it, x0);
+ 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();
+ if(!hFrame){
+ hFrame=new TH1I("hFrame", "", 100, -.3, .3);
+ hFrame->SetMinimum(-.1);hFrame->SetMaximum(.1);
+ hFrame->SetXTitle("tg#phi-htg#theta");
+ hFrame->SetYTitle("#Deltay[cm]");
+ hFrame->SetLineColor(1);hFrame->SetLineWidth(1);
+ hFrame->Draw();
+ } else hFrame->Reset();
+ gm->Draw("pl"); line.Draw("same");
+ fCanvas->Modified(); fCanvas->Update();
+ if(IsSaveAs()) fCanvas->SaveAs(Form("Figures/ProcessMean_Z[%5.3f]_X[%5.3f].gif", fZ, fX));
+ else gSystem->Sleep(100);
+ printf(" xd=%4.2f[cm] dx=%5.3e[cm] dy=%5.3e[cm]\n", fX, fR[0], fR[2]);
}
}
+
+ // 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");
- TH1D *h=hdx->ProjectionY("py"); h->Scale(1./hdx->GetNbinsX());
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff