#include "TObjArray.h"
#include "TAxis.h"
+#include "TF1.h"
+#include "TGraphErrors.h"
#include "TH2I.h"
#include "TMath.h"
AliTRDclusterResolution::AliTRDclusterResolution()
: AliTRDrecoTask("CalibClRes", "Calibrate Cluster Resolution for Tracking")
,fInfo(0x0)
+ ,fResults(0x0)
+ ,fAt(0x0)
+ ,fAd(0x0)
+ ,fExB(0.)
{
+ fAt = new TAxis(kNTB, -0.075, (kNTB-.5)*.15);
+ fAd = new TAxis(kND, 0., .25);
}
//_______________________________________________________
AliTRDclusterResolution::~AliTRDclusterResolution()
{
-
+ if(fAd) delete fAd;
+ if(fAt) delete fAt;
+ if(fResults){
+ fResults->Delete();
+ delete fResults;
+ }
}
//_______________________________________________________
}
//_______________________________________________________
-void AliTRDclusterResolution::GetRefFigure(Int_t /*ifig*/)
+void AliTRDclusterResolution::GetRefFigure(Int_t ifig)
{
-
+ if(!fResults) return /*kFALSE*/;
+
+ TObjArray *arr = 0x0;
+ TH2 *h2 = 0x0;
+ TGraphErrors *gm(0x0), *gs(0x0);
+ switch(ifig){
+ case kQRes:
+ if(!(arr = (TObjArray*)fResults->At(kQRes))) break;
+ if(!(gm = (TGraphErrors*)arr->At(0))) break;
+ if(!(gs = (TGraphErrors*)arr->At(1))) break;
+ gs->Draw("apl");
+ gs->GetHistogram()->SetXTitle("Log(Q) [a.u.]");
+ gs->GetHistogram()->SetYTitle("#sigma_{y}^{2} [mm^{2}] [a.u.]");
+ return /*kTRUE*/;
+ case kYRes:
+ if(!(arr = (TObjArray*)fResults->At(kYRes))) break;
+ if(!(gm = (TGraphErrors*)arr->At(0))) break;
+ if(!(gs = (TGraphErrors*)arr->At(1))) break;
+ gs->Draw("apl");
+ gs->GetHistogram()->SetXTitle("y [w]");
+ gs->GetHistogram()->SetYTitle("#sigma_{y}^{2} [mm^{2}] [a.u.]");
+ return /*kTRUE*/;
+ case kSXRes:
+ if(!(h2 = (TH2F*)fResults->At(kSXRes))) break;
+ h2->Draw("lego2");
+ return /*kTRUE*/;
+ case kSYRes:
+ if(!(h2 = (TH2F*)fResults->At(kSYRes))) break;
+ h2->Draw("lego2");
+ return /*kTRUE*/;
+ default:
+ break;
+ }
+
+ return /*kFALSE*/;
}
//_______________________________________________________
fContainer = new TObjArray(kN+1);
//fContainer->SetOwner(kTRUE);
- TAxis at(kNTB, -0.075, (kNTB-.5)*.15);
- TAxis ad(kND, 0., .25);
Int_t ih = 0;
- for(Int_t id=1; id<=ad.GetNbins(); id++){
- for(Int_t it=1; it<=at.GetNbins(); it++){
- fContainer->AddAt(new TH2I(Form("h_d%02dt%02d", id, it), Form("d_{wire}(%5.2f-%5.2f)[mm] x_{drift}(%5.2f-%5.2f)[mm]", 10.*ad.GetBinCenter(id)- 5.*ad.GetBinWidth(id), 10.*ad.GetBinCenter(id)+ 5.*ad.GetBinWidth(id), 10.*at.GetBinCenter(it)- 5.*at.GetBinWidth(it), 10.*at.GetBinCenter(it)+ 5.*at.GetBinWidth(it)), 30, -.15, .15, 100, -.5, .5), ih++);
+ for(Int_t id=1; id<=fAd->GetNbins(); id++){
+ for(Int_t it=1; it<=fAt->GetNbins(); it++){
+ fContainer->AddAt(new TH2I(Form("h_d%02dt%02d", id, it), Form("d_{wire}(%5.2f-%5.2f)[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)), 30, -.15, .15, 100, -.5, .5), ih++);
}
}
fContainer->AddAt(new TH2I("h_q", "", 50, 2.2, 7.5, 100, -.5, .5), ih++);
+ fContainer->AddAt(new TH2I("h_y", "", 50, -.35, .35, 100, -.5, .5), ih++);
return fContainer;
}
void AliTRDclusterResolution::Exec(Option_t *)
{
Int_t det, t;
- Float_t x, y, z, q, dy, dydx, dzdx, cov[3];
- TAxis at(kNTB, -0.075, (kNTB-.5)*.15); Int_t it = 0;
- TAxis ad(kND, 0., .25); Int_t id = 0;
+ Float_t x, y, z, q, dy, dydx, dzdx, cov[3], covcl[3];
TH2I *h2 = 0x0;
const AliTRDclusterInfo *cli = 0x0;
TIterator *iter=fInfo->MakeIterator();
while((cli=dynamic_cast<AliTRDclusterInfo*>((*iter)()))){
dy = cli->GetResolution();
- it = at.FindBin(cli->GetDriftLength());
- if(it==0 || it == at.GetNbins()+1){
+ 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;
}
- id = ad.FindBin(cli->GetAnisochronity());
- if(id==0 || id == ad.GetNbins()+1){
+ 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()));
continue;
}
}
cli->GetGlobalPosition(y, z, dydx, dzdx, &cov[0]);
- h2->Fill(dydx, dy);
+ h2->Fill(dydx, cov[0]!=0. ? dy/TMath::Sqrt(cov[0]) : dy);
- // resolution as a function of cluster charge
+ // resolution as a function of:
+ // - cluster charge and
+ // - y displacement
// only for phi equal exB
- if(TMath::Abs(dydx)<.01){
- cli->GetCluster(det, x, y, z, q, t);
+ if(TMath::Abs(dydx-fExB)<.01){
+ cli->GetCluster(det, x, y, z, q, t, covcl);
h2 = (TH2I*)fContainer->At(kN);
h2->Fill(TMath::Log(q), dy);
+
+ h2 = (TH2I*)fContainer->At(kN+1);
+ h2->Fill(cli->GetYDisplacement(), dy);
}
}
PostData(0, fContainer);
//_______________________________________________________
Bool_t AliTRDclusterResolution::PostProcess()
{
- return kFALSE;
+ if(!fContainer) return kFALSE;
+
+ TH2 *h2 = 0x0;
+ if(!fResults){
+ TGraphErrors *g = 0x0;
+ fResults = new TObjArray(4);
+ fResults->SetOwner();
+ TObjArray *arr = 0x0;
+ fResults->AddAt(arr = new TObjArray(2), kQRes);
+ arr->SetOwner();
+ arr->AddAt(g = new TGraphErrors(), 0);
+ g->SetLineColor(kBlue); g->SetMarkerColor(kBlue);
+ g->SetMarkerStyle(7);
+ arr->AddAt(g = new TGraphErrors(), 1);
+ g->SetLineColor(kRed); g->SetMarkerColor(kRed);
+ g->SetMarkerStyle(23);
+
+ fResults->AddAt(arr = new TObjArray(2), kYRes);
+ arr->SetOwner();
+ arr->AddAt(g = new TGraphErrors(), 0);
+ g->SetLineColor(kBlue); g->SetMarkerColor(kBlue);
+ g->SetMarkerStyle(7);
+ arr->AddAt(g = new TGraphErrors(), 1);
+ g->SetLineColor(kRed); g->SetMarkerColor(kRed);
+ g->SetMarkerStyle(23);
+
+ fResults->AddAt(h2 = new TH2F("hSX", "",
+ fAd->GetNbins(), fAd->GetXmin(), fAd->GetXmax(),
+ fAt->GetNbins(), fAt->GetXmin(), fAt->GetXmax()), kSXRes);
+ h2->SetXTitle("d [mm]");
+ h2->SetYTitle("x [mm]");
+ h2->SetZTitle("#sigma_{x}^{2} [mm^{2}]");
+ fResults->AddAt(h2 = (TH2F*)h2->Clone("hSY"), kSYRes);
+ h2->SetZTitle("#sigma_{y}^{2} [mm^{2}]");
+ } else {
+ TObject *o = 0x0;
+ TIterator *iter=fResults->MakeIterator();
+ while((o=((*iter)()))) o->Clear();
+ }
+
+ TF1 f("f", "gaus", -.5, .5);
+ TF1 sig("sig", "pol1", 0., .0225);
+
+ TAxis *ax = 0x0;
+ TH1D *h1 = 0x0;
+
+ // process resolution dependency on charge
+ if((h2 = (TH2I*)fContainer->At(kN))) {
+ TObjArray *arr = (TObjArray*)fResults->At(kQRes);
+ TGraphErrors *gqm = (TGraphErrors*)arr->At(0);
+ TGraphErrors *gqs = (TGraphErrors*)arr->At(1);
+ ax = h2->GetXaxis();
+ for(Int_t ix=1; ix<=ax->GetNbins(); ix++){
+ Float_t logq = ax->GetBinCenter(ix);
+ h1 = h2->ProjectionY("py", ix, ix);
+ if(h1->GetEntries() < 50) continue;
+ Adjust(&f, h1);
+ h1->Fit(&f, "Q");
+
+ // Fill sy^2 = f(log(q))
+ Int_t ip = gqm->GetN();
+ gqm->SetPoint(ip, logq, 10.*f.GetParameter(1));
+ gqm->SetPointError(ip, 0., 10.*f.GetParError(1));
+ gqs->SetPoint(ip, logq, 1.e2*f.GetParameter(2)*f.GetParameter(2));
+ gqs->SetPointError(ip, 0., 2.e2*f.GetParameter(2)*f.GetParError(2));
+ }
+ } else AliWarning("Missing dy=f(Q) histo");
+
+ // process resolution dependency on y displacement
+ if((h2 = (TH2I*)fContainer->At(kN+1))) {
+ TObjArray *arr = (TObjArray*)fResults->At(kYRes);
+ TGraphErrors *gym = (TGraphErrors*)arr->At(0);
+ TGraphErrors *gys = (TGraphErrors*)arr->At(1);
+ 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^2 = f(log(q))
+ Int_t ip = gym->GetN();
+ gym->SetPoint(ip, yd, 10.*f.GetParameter(1));
+ gym->SetPointError(ip, 0., 10.*f.GetParError(1));
+ gys->SetPoint(ip, yd, 1.e2*f.GetParameter(2)*f.GetParameter(2));
+ gys->SetPointError(ip, 0., 2.e2*f.GetParameter(2)*f.GetParError(2));
+ }
+ } else AliWarning("Missing dy=f(y_d) histo");
+
+ // process resolution dependency on drift legth and drift cell width
+ TGraphErrors *gm = new TGraphErrors(), *gs = new TGraphErrors();
+ Float_t d(0.), x(0.);
+ TH2F *hsx = (TH2F*)fResults->At(kSXRes);
+ TH2F *hsy = (TH2F*)fResults->At(kSYRes);
+ 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 x=%f[mm]\n", x);
+ Int_t idx = (id-1)*kNTB+it-1;
+ if(!(h2 = (TH2I*)fContainer->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++){
+ Float_t dydx = ax->GetBinCenter(ix) - fExB;
+ if(dydx<0.) continue;
+ h1 = h2->ProjectionY("py", ix, ix);
+ if(h1->GetEntries() < 50) continue;
+ Adjust(&f, h1);
+ h1->Fit(&f, "Q", "", -.2, .2);
+
+ // Fill sy^2 = f(tg_phi^2)
+ gm->SetPoint(ip, dydx, 10.*f.GetParameter(1));
+ gm->SetPointError(ip, 0., 10.*f.GetParError(1));
+ gs->SetPoint(ip, dydx*dydx, 1.e2*f.GetParameter(2)*f.GetParameter(2));
+ gs->SetPointError(ip, 0., 2.e2*f.GetParameter(2)*f.GetParError(2));
+ ip++;
+ }
+ for(;ip<gm->GetN(); ip++){
+ gm->RemovePoint(ip);gs->RemovePoint(ip);
+ }
+ gs->Fit(&sig, "Q");
+ hsx->SetBinContent(id, it, sig.GetParameter(1));
+ hsx->SetBinError(id, it, sig.GetParError(1));
+ hsy->SetBinContent(id, it, sig.GetParameter(0));
+ hsy->SetBinError(id, it, sig.GetParError(0));
+ }
+ }
+ delete gm; delete gs;
+
+ return kTRUE;
}
AliTRDclusterInfo();
Float_t GetAnisochronity() const {return fD;}
- inline void GetCluster(Int_t &det, Float_t &x, Float_t &y, Float_t &z, Float_t &q, Int_t &t) const;
+ inline void GetCluster(Int_t &det, Float_t &x, Float_t &y, Float_t &z, Float_t &q, Int_t &t, Float_t *cov) const;
void GetMC(Int_t &pdg, Int_t &label) const{
pdg = fPdg;
label= fLbl; }
memcpy(cov, fCov, 3*sizeof(Float_t));}
Float_t GetResolution() const {return fdy;}
Float_t GetDriftLength() const {return fXd;}
+ Float_t GetYDisplacement() const {return fYd;}
void Print(Option_t *opt="") const;
void SetAnisochronity(Float_t d) {fD = d;}
- void SetCluster(const AliTRDcluster *c=0x0);
+ void SetCluster(const AliTRDcluster *c, Float_t *cov=0x0);
void SetMC(Int_t pdg, Int_t label){
fPdg = pdg;
fLbl = label;}
private:
UShort_t fDet; // detector
Short_t fPdg; // particle code
- Short_t fLbl; // track label (MC)
- Short_t fLocalTime; // calibrate drift time
+ Short_t fLbl; // track label (MC)
+ Short_t fLocalTime; // calibrate drift time
Float_t fQ; // cluster charge (REC)
Float_t fX; // x coordinate (REC)
Float_t fY; // y coordinate (REC)
+ Float_t fYd; // displacement from pad center y coordinate
Float_t fZ; // z coordinate (REC)
Float_t fdydx; // slope in phi (MC)
Float_t fdzdx; // slope in theta (MC)
Float_t fXd; // drift length
Float_t fYt; // y coordinate (MC)
Float_t fZt; // z coordinate (MC)
- Float_t fCov[3];// covariance matrix in the yz plane (global)
+ Float_t fCov[3];// covariance matrix in the yz plane (track)
+ Float_t fCovCl[3];// covariance matrix in the yz plane (cluster)
Float_t fdy; // difference in y after tilt correction
Float_t fD; // distance to the anode wire
//_________________________________________________
-inline void AliTRDclusterInfo::GetCluster(Int_t &det, Float_t &x, Float_t &y, Float_t &z, Float_t &q, Int_t &t) const
+inline void AliTRDclusterInfo::GetCluster(Int_t &det, Float_t &x, Float_t &y, Float_t &z, Float_t &q, Int_t &t, Float_t *cov) const
{
det = fDet;
x = fX;
z = fZ;
q = fQ;
t = fLocalTime;
+ memcpy(cov, fCovCl, 3*sizeof(Float_t));
}
#endif
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff