TRD/qaRec/AliTRDclusterResolution.cxx

   1 #include "AliTRDclusterResolution.h"
   2 #include "AliTRDtrackInfo/AliTRDclusterInfo.h"
   3
   4 #include "AliLog.h"
   5
   6 #include "TObjArray.h"
   7 #include "TAxis.h"
   8 #include "TF1.h"
   9 #include "TGraphErrors.h"
  10 #include "TH2I.h"
  11 #include "TMath.h"
  12
  13 ClassImp(AliTRDclusterResolution)
  14
  15
  16 //_______________________________________________________
  17 AliTRDclusterResolution::AliTRDclusterResolution()
  18   : AliTRDrecoTask("CalibClRes", "Calibrate Cluster Resolution for Tracking")
  19   ,fInfo(0x0)
  20   ,fResults(0x0)
  21   ,fAt(0x0)
  22   ,fAd(0x0)
  23   ,fExB(0.)
  24 {
  25   fAt = new TAxis(kNTB, -0.075, (kNTB-.5)*.15);
  26   fAd = new TAxis(kND, 0., .25);
  27 }
  28
  29 //_______________________________________________________
  30 AliTRDclusterResolution::~AliTRDclusterResolution()
  31 {
  32   if(fAd) delete fAd;
  33   if(fAt) delete fAt;
  34   if(fResults){
  35     fResults->Delete();
  36     delete fResults;
  37   }
  38 }
  39
  40 //_______________________________________________________
  41 void AliTRDclusterResolution::ConnectInputData(Option_t *)
  42 {
  43   fInfo = dynamic_cast<TObjArray *>(GetInputData(0));
  44 }
  45
  46 //_______________________________________________________
  47 void AliTRDclusterResolution::CreateOutputObjects()
  48 {
  49   OpenFile(0, "RECREATE");
  50   fContainer = Histos();
  51 }
  52
  53 //_______________________________________________________
  54 void AliTRDclusterResolution::GetRefFigure(Int_t ifig)
  55 {
  56   if(!fResults) return /*kFALSE*/;
  57
  58   TObjArray *arr = 0x0;
  59   TH2 *h2 = 0x0;
  60   TGraphErrors *gm(0x0), *gs(0x0);
  61   switch(ifig){
  62   case kQRes:
  63     if(!(arr = (TObjArray*)fResults->At(kQRes))) break;
  64     if(!(gm = (TGraphErrors*)arr->At(0))) break;
  65     if(!(gs = (TGraphErrors*)arr->At(1))) break;
  66     gs->Draw("apl");
  67     gs->GetHistogram()->SetXTitle("Log(Q) [a.u.]");
  68     gs->GetHistogram()->SetYTitle("#sigma_{y}^{2} [mm^{2}] [a.u.]");
  69     return /*kTRUE*/;
  70   case kYRes:
  71     if(!(arr = (TObjArray*)fResults->At(kYRes))) break;
  72     if(!(gm = (TGraphErrors*)arr->At(0))) break;
  73     if(!(gs = (TGraphErrors*)arr->At(1))) break;
  74     gs->Draw("apl");
  75     gs->GetHistogram()->SetXTitle("y [w]");
  76     gs->GetHistogram()->SetYTitle("#sigma_{y}^{2} [mm^{2}] [a.u.]");
  77     return /*kTRUE*/;
  78   case kSXRes:
  79     if(!(h2 = (TH2F*)fResults->At(kSXRes))) break;
  80     h2->Draw("lego2");
  81     return /*kTRUE*/;
  82   case kSYRes:
  83     if(!(h2 = (TH2F*)fResults->At(kSYRes))) break;
  84     h2->Draw("lego2");
  85     return /*kTRUE*/;
  86   default:
  87     break;
  88   }
  89
  90   return /*kFALSE*/;
  91 }
  92
  93 //_______________________________________________________
  94 TObjArray* AliTRDclusterResolution::Histos()
  95 {
  96   if(fContainer) return fContainer;
  97   fContainer = new TObjArray(kN+1);
  98   //fContainer->SetOwner(kTRUE);
  99
 100   Int_t ih = 0;
 101   for(Int_t id=1; id<=fAd->GetNbins(); id++){
 102     for(Int_t it=1; it<=fAt->GetNbins(); it++){
 103       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++);
 104     }
 105   }
 106   fContainer->AddAt(new TH2I("h_q", "", 50, 2.2, 7.5, 100, -.5, .5), ih++);
 107   fContainer->AddAt(new TH2I("h_y", "", 50, -.35, .35, 100, -.5, .5), ih++);
 108   return fContainer;
 109 }
 110
 111 //_______________________________________________________
 112 void AliTRDclusterResolution::Exec(Option_t *)
 113 {
 114   Int_t det, t;
 115   Float_t x, y, z, q, dy, dydx, dzdx, cov[3], covcl[3];
 116   TH2I *h2 = 0x0;
 117   const AliTRDclusterInfo *cli = 0x0;
 118   TIterator *iter=fInfo->MakeIterator();
 119   while((cli=dynamic_cast<AliTRDclusterInfo*>((*iter)()))){
 120     dy = cli->GetResolution();
 121     Int_t it = fAt->FindBin(cli->GetDriftLength());
 122     if(it==0 || it == fAt->GetNbins()+1){
 123       AliWarning(Form("Drift length %f outside allowed range", cli->GetDriftLength()));
 124       continue;
 125     }
 126     Int_t id = fAd->FindBin(cli->GetAnisochronity());
 127     if(id==0 || id == fAd->GetNbins()+1){
 128       AliWarning(Form("Distance to anode %f outside allowed range", cli->GetAnisochronity()));
 129       continue;
 130     }
 131     if(!(h2 = (TH2I*)fContainer->At((id-1)*kNTB+it-1))){
 132       AliWarning(Form("Missing histo at index idx[%3d] [id[%2d] it[%2d]] xd[%f] d[%f]\n", (id-1)*kNTB+it-1, id, it, cli->GetDriftLength(), cli->GetAnisochronity()));
 133       continue;
 134     }
 135
 136     cli->GetGlobalPosition(y, z, dydx, dzdx, &cov[0]);
 137     h2->Fill(dydx, cov[0]!=0. ? dy/TMath::Sqrt(cov[0]) : dy);
 138
 139     // resolution as a function of:
 140     //   - cluster charge and
 141     //   - y displacement
 142     // only for phi equal exB
 143     if(TMath::Abs(dydx-fExB)<.01){
 144       cli->GetCluster(det, x, y, z, q, t, covcl);
 145       h2 = (TH2I*)fContainer->At(kN);
 146       h2->Fill(TMath::Log(q), dy);
 147
 148       h2 = (TH2I*)fContainer->At(kN+1);
 149       h2->Fill(cli->GetYDisplacement(), dy);
 150     }
 151   }
 152   PostData(0, fContainer);
 153 }
 154
 155
 156 //_______________________________________________________
 157 Bool_t AliTRDclusterResolution::PostProcess()
 158 {
 159   if(!fContainer) return kFALSE;
 160
 161   TH2 *h2 = 0x0;
 162   if(!fResults){
 163     TGraphErrors *g = 0x0;
 164     fResults = new TObjArray(4);
 165     fResults->SetOwner();
 166     TObjArray *arr = 0x0;
 167     fResults->AddAt(arr = new TObjArray(2), kQRes);
 168     arr->SetOwner();
 169     arr->AddAt(g = new TGraphErrors(), 0);
 170     g->SetLineColor(kBlue); g->SetMarkerColor(kBlue);
 171     g->SetMarkerStyle(7);
 172     arr->AddAt(g = new TGraphErrors(), 1);
 173     g->SetLineColor(kRed); g->SetMarkerColor(kRed);
 174     g->SetMarkerStyle(23);
 175
 176     fResults->AddAt(arr = new TObjArray(2), kYRes);
 177     arr->SetOwner();
 178     arr->AddAt(g = new TGraphErrors(), 0);
 179     g->SetLineColor(kBlue); g->SetMarkerColor(kBlue);
 180     g->SetMarkerStyle(7);
 181     arr->AddAt(g = new TGraphErrors(), 1);
 182     g->SetLineColor(kRed); g->SetMarkerColor(kRed);
 183     g->SetMarkerStyle(23);
 184
 185     fResults->AddAt(h2 = new TH2F("hSX", "",
 186       fAd->GetNbins(), fAd->GetXmin(), fAd->GetXmax(),
 187       fAt->GetNbins(), fAt->GetXmin(), fAt->GetXmax()), kSXRes);
 188     h2->SetXTitle("d [mm]");
 189     h2->SetYTitle("x [mm]");
 190     h2->SetZTitle("#sigma_{x}^{2} [mm^{2}]");
 191     fResults->AddAt(h2 = (TH2F*)h2->Clone("hSY"), kSYRes);
 192     h2->SetZTitle("#sigma_{y}^{2} [mm^{2}]");
 193   } else {
 194     TObject *o = 0x0;
 195     TIterator *iter=fResults->MakeIterator();
 196     while((o=((*iter)()))) o->Clear();
 197   }
 198
 199   TF1 f("f", "gaus", -.5, .5);
 200   TF1 sig("sig", "pol1", 0., .0225);
 201
 202   TAxis *ax = 0x0;
 203   TH1D *h1 = 0x0;
 204
 205   // process resolution dependency on charge
 206   if((h2 = (TH2I*)fContainer->At(kN))) {
 207     TObjArray *arr = (TObjArray*)fResults->At(kQRes);
 208     TGraphErrors *gqm = (TGraphErrors*)arr->At(0);
 209     TGraphErrors *gqs = (TGraphErrors*)arr->At(1);
 210     ax = h2->GetXaxis();
 211     for(Int_t ix=1; ix<=ax->GetNbins(); ix++){
 212       Float_t logq = ax->GetBinCenter(ix);
 213       h1 = h2->ProjectionY("py", ix, ix);
 214       if(h1->GetEntries() < 50) continue;
 215       Adjust(&f, h1);
 216       h1->Fit(&f, "Q");
 217
 218       // Fill sy^2 = f(log(q))
 219       Int_t ip = gqm->GetN();
 220       gqm->SetPoint(ip, logq, 10.*f.GetParameter(1));
 221       gqm->SetPointError(ip, 0., 10.*f.GetParError(1));
 222       gqs->SetPoint(ip, logq, 1.e2*f.GetParameter(2)*f.GetParameter(2));
 223       gqs->SetPointError(ip, 0., 2.e2*f.GetParameter(2)*f.GetParError(2));
 224     }
 225   } else AliWarning("Missing dy=f(Q) histo");
 226
 227   // process resolution dependency on y displacement
 228   if((h2 = (TH2I*)fContainer->At(kN+1))) {
 229     TObjArray *arr = (TObjArray*)fResults->At(kYRes);
 230     TGraphErrors *gym = (TGraphErrors*)arr->At(0);
 231     TGraphErrors *gys = (TGraphErrors*)arr->At(1);
 232     ax = h2->GetXaxis();
 233     for(Int_t ix=1; ix<=ax->GetNbins(); ix++){
 234       Float_t yd = ax->GetBinCenter(ix);
 235       h1 = h2->ProjectionY("py", ix, ix);
 236       if(h1->GetEntries() < 50) continue;
 237       Adjust(&f, h1);
 238       h1->Fit(&f, "Q");
 239
 240       // Fill sy^2 = f(log(q))
 241       Int_t ip = gym->GetN();
 242       gym->SetPoint(ip, yd, 10.*f.GetParameter(1));
 243       gym->SetPointError(ip, 0., 10.*f.GetParError(1));
 244       gys->SetPoint(ip, yd, 1.e2*f.GetParameter(2)*f.GetParameter(2));
 245       gys->SetPointError(ip, 0., 2.e2*f.GetParameter(2)*f.GetParError(2));
 246     }
 247   } else AliWarning("Missing dy=f(y_d) histo");
 248
 249   // process resolution dependency on drift legth and drift cell width
 250   TGraphErrors *gm = new TGraphErrors(), *gs = new TGraphErrors();
 251   Float_t d(0.), x(0.);
 252   TH2F *hsx = (TH2F*)fResults->At(kSXRes);
 253   TH2F *hsy = (TH2F*)fResults->At(kSYRes);
 254   for(Int_t id=1; id<=fAd->GetNbins(); id++){
 255     d = fAd->GetBinCenter(id); //[mm]
 256     printf("Doing d=%f[mm]\n", d);
 257     for(Int_t it=1; it<=fAt->GetNbins(); it++){
 258       x = fAt->GetBinCenter(it); //[mm]
 259       printf("Doing x=%f[mm]\n", x);
 260       Int_t idx = (id-1)*kNTB+it-1;
 261       if(!(h2 = (TH2I*)fContainer->At(idx))) {
 262         AliWarning(Form("Missing histo at index idx[%3d] [id[%2d] it[%2d]] xd[%f] d[%f]\n", idx, id, it, x, d));
 263         continue;
 264       }
 265
 266       Int_t ip = 0;
 267       ax = h2->GetXaxis();
 268       for(Int_t ix=1; ix<=ax->GetNbins(); ix++){
 269         Float_t dydx = ax->GetBinCenter(ix) - fExB;
 270         if(dydx<0.) continue;
 271         h1 = h2->ProjectionY("py", ix, ix);
 272         if(h1->GetEntries() < 50) continue;
 273         Adjust(&f, h1);
 274         h1->Fit(&f, "Q", "", -.2, .2);
 275
 276         // Fill sy^2 = f(tg_phi^2)
 277         gm->SetPoint(ip, dydx, 10.*f.GetParameter(1));
 278         gm->SetPointError(ip, 0., 10.*f.GetParError(1));
 279         gs->SetPoint(ip, dydx*dydx, 1.e2*f.GetParameter(2)*f.GetParameter(2));
 280         gs->SetPointError(ip, 0., 2.e2*f.GetParameter(2)*f.GetParError(2));
 281         ip++;
 282       }
 283       for(;ip<gm->GetN(); ip++){
 284         gm->RemovePoint(ip);gs->RemovePoint(ip);
 285       }
 286       gs->Fit(&sig, "Q");
 287       hsx->SetBinContent(id, it, sig.GetParameter(1));
 288       hsx->SetBinError(id, it, sig.GetParError(1));
 289       hsy->SetBinContent(id, it, sig.GetParameter(0));
 290       hsy->SetBinError(id, it, sig.GetParError(0));
 291     }
 292   }
 293   delete gm; delete gs;
 294
 295   return kTRUE;
 296 }
 297
 298