PWG1/TRD/AliTRDclusterResolution.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercialf purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 /* $Id: AliTRDclusterResolution.cxx */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //                                                                           //
  20 //  TRD cluster error parameterization                                        //
  21 //                                                                           //
  22 // This class is designed to produce the reference plots for a detailed study//
  23 // and parameterization of TRD cluster errors. The following effects are taken//
  24 // into account :                                                            //
  25 //   - dependence with the total charge of the cluster                       //
  26 //   - dependence with the distance from the center pad. This is monitored
  27 // for each layer individually since the pad size varies with layer
  28 //   - dependence with the drift length - here the influence of anisochronity
  29 // and diffusion are searched
  30 //   - dependence with the distance to the anode wire - anisochronity effects
  31 //   - dependence with track angle (for y resolution)
  32 // The correlation between effects is taken into account.
  33 //
  34 // Since magnetic field plays a very important role in the TRD measurement
  35 // the ExB correction is forced by the setter function SetExB(Int_t). The
  36 // argument is the detector index, if none is specified all will be
  37 // considered.
  38 //
  39 // Two cases are of big importance.
  40 //   - comparison with MC
  41 //   - comparison with Kalman fit. In this case the covariance matrix of the
  42 // Kalman fit are needed.
  43 //
  44 // The functionalities implemented in this class are based on the storage
  45 // class AliTRDclusterInfo.
  46 //
  47 // The Method
  48 // ----------
  49 //
  50 // The method to disentangle s_y and s_x is based on the relation (see also fig.)
  51 // BEGIN_LATEX
  52 // #sigma^{2} = #sigma^{2}_{y} + tg^{2}(#alpha_{L})*#sigma^{2}_{x_{d}} + tg^{2}(#phi-#alpha_{L})*(#sigma^{2}_{x_{d}}+#sigma^{2}_{x_{c}})
  53 // END_LATEX
  54 // with
  55 // BEGIN_LATEX
  56 // #sigma^{2}_{x_{c}} #approx 0
  57 // END_LATEX
  58 // we suppose the chamber is well calibrated for t_{0} and aligned in
  59 // radial direction.
  60 //
  61 // Clusters can be radially shifted due to three causes:
  62 //   - globally shifted - due to residual misalignment/miscalibration(t0)
  63 //   - locally shifted - due to different local drift velocity from the mean
  64 //   - randomly shifted - due to neighboring (radial direction) clusters
  65 // charge induced by asymmetry of the TRF.
  66 //
  67 // We estimate this effects by the relations:
  68 // BEGIN_LATEX
  69 // #mu_{y} = tg(#alpha_{L})*#Delta x_{d}(...) + tg(#phi-#alpha_{L})*(#Delta x_{c}(...) + #Delta x_{d}(...))
  70 // END_LATEX
  71 // where
  72 // BEGIN_LATEX
  73 // #Delta x_{d}(...) = (<v_{d}> + #delta v_{d}(x_{d}, d)) * (t + t^{*}(Q))
  74 // END_LATEX
  75 // and we specified explicitely the variation of drift velocity parallel
  76 // with the track (x_{d}) and perpendicular to it due to anisochronity (d).
  77 //
  78 // For estimating the contribution from asymmetry of TRF the following
  79 // parameterization is being used
  80 // BEGIN_LATEX
  81 // t^{*}(Q) = #delta_{0} * #frac{Q_{t+1} - Q_{t-1}}{Q_{t-1} + Q_{t} + Q_{t+1}}
  82 // END_LATEX
  83 //
  84 //
  85 // Clusters can also be r-phi shifted due to:
  86 //   - wrong PRF or wrong cuts at digits level
  87 //The following correction is applied :
  88 // BEGIN_LATEX
  89 // <#Delta y> = a + b * sin(c*y_{pw})
  90 // END_LATEX
  91
  92 // The Models
  93 //
  94 //   Parameterization against total charge
  95 //
  96 // Obtained for B=0T at phi=0. All other effects integrated out.
  97 // BEGIN_LATEX
  98 // #sigma^{2}_{y}(Q) = #sigma^{2}_{y}(...) + b(#frac{1}{Q} - #frac{1}{Q_{0}})
  99 // END_LATEX
 100 // For B diff 0T the error of the average ExB correction error has to be subtracted !!
 101 //
 102 //   Parameterization Sx
 103 //
 104 // The parameterization of the error in the x direction can be written as
 105 // BEGIN_LATEX
 106 // #sigma_{x} = #sigma_{x}^{||} + #sigma_{x}^{#perp}
 107 // END_LATEX
 108 //
 109 // where the parallel component is given mainly by the TRF width while
 110 // the perpendicular component by the anisochronity. The model employed for
 111 // the parallel is gaus(0)+expo(3) with the following parameters
 112 // 1  C   5.49018e-01   1.23854e+00   3.84540e-04  -8.21084e-06
 113 // 2  M   7.82999e-01   6.22531e-01   2.71272e-04  -6.88485e-05
 114 // 3  S   2.74451e-01   1.13815e+00   2.90667e-04   1.13493e-05
 115 // 4  E1  2.53596e-01   1.08646e+00   9.95591e-05  -2.11625e-05
 116 // 5  E2 -2.40078e-02   4.26520e-01   4.67153e-05  -2.35392e-04
 117 //
 118 // and perpendicular to the track is pol2 with the parameters
 119 //
 120 // Par_0 = 0.190676 +/- 0.41785
 121 // Par_1 = -3.9269  +/- 7.49862
 122 // Par_2 = 14.7851  +/- 27.8012
 123 //
 124 //   Parameterization Sy
 125 //
 126 // The parameterization of the error in the y direction along track uses
 127 // BEGIN_LATEX
 128 // #sigma_{y}^{||} = #sigma_{y}^{0} -a*exp(1/(x-b))
 129 // END_LATEX
 130 //
 131 // with following values for the parameters:
 132 // 1  sy0 2.60967e-01   2.99652e-03   7.82902e-06  -1.89636e-04
 133 // 2  a  -7.68941e+00   1.87883e+00   3.84539e-04   9.38268e-07
 134 // 3  b  -3.41160e-01   7.72850e-02   1.63231e-05   2.51602e-05
 135 //
 136 //==========================================================================
 137 // Example how to retrive reference plots from the task
 138 // void steerClErrParam(Int_t fig=0)
 139 // {
 140 //   gSystem->Load("libANALYSIS.so");
 141 //   gSystem->Load("libTRDqaRec.so");
 142 //
 143 //   // initialize DB manager
 144 //   AliCDBManager *cdb = AliCDBManager::Instance();
 145 //   cdb->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
 146 //   cdb->SetRun(0);
 147 //   // initialize magnetic field.
 148 //   AliMagFCheb *field=new AliMagFCheb("Maps","Maps", 2, 1., 10., AliMagFCheb::k5kG);
 149 //   AliTracker::SetFieldMap(field, kTRUE);
 150 //
 151 //   AliTRDclusterResolution *res = new AliTRDclusterResolution();
 152 //   res->SetMCdata();
 153 //   res->Load("TRD.TaskClErrParam.root");
 154 //   res->SetExB();
 155 //   res->SetVisual();
 156 //   //res->SetSaveAs();
 157 //   res->SetProcessCharge(kFALSE);
 158 //   res->SetProcessCenterPad(kFALSE);
 159 //   //res->SetProcessMean(kFALSE);
 160 //   res->SetProcessSigma(kFALSE);
 161 //   if(!res->PostProcess()) return;
 162 //   new TCanvas;
 163 //   res->GetRefFigure(fig);
 164 // }
 165 //
 166 //  Authors:                                                              //
 167 //    Alexandru Bercuci <A.Bercuci@gsi.de>                                //
 168 ////////////////////////////////////////////////////////////////////////////
 169
 170 #include "AliTRDclusterResolution.h"
 171 #include "AliTRDresolution.h"
 172 #include "AliTRDinfoGen.h"
 173 #include "info/AliTRDclusterInfo.h"
 174
 175 #include "AliTRDcalibDB.h"
 176 #include "Cal/AliTRDCalROC.h"
 177 #include "Cal/AliTRDCalDet.h"
 178 #include "AliTRDCommonParam.h"
 179 #include "AliTRDgeometry.h"
 180 #include "AliTRDpadPlane.h"
 181 #include "AliTRDcluster.h"
 182 #include "AliTRDseedV1.h"
 183
 184 #include "AliESDEvent.h"
 185 #include "AliCDBManager.h"
 186
 187 #include "TROOT.h"
 188 #include "TSystem.h"
 189 #include "TMath.h"
 190 #include "TLinearFitter.h"
 191 #include "TGeoGlobalMagField.h"
 192 #include <TGeoMatrix.h>
 193 #include "TObjArray.h"
 194 #include "TTree.h"
 195 #include "TH2I.h"
 196 #include "TH3S.h"
 197 #include "TAxis.h"
 198 #include "TF1.h"
 199 #include "TCanvas.h"
 200 #include "TLegend.h"
 201 #include "TGraphErrors.h"
 202 #include "TLine.h"
 203
 204 ClassImp(AliTRDclusterResolution)
 205
 206 const Float_t AliTRDclusterResolution::fgkTimeBinLength = 1./ AliTRDCommonParam::Instance()->GetSamplingFrequency();
 207 //_______________________________________________________
 208 AliTRDclusterResolution::AliTRDclusterResolution()
 209   : AliTRDrecoTask()
 210   ,fCanvas(NULL)
 211   ,fInfo(NULL)
 212   ,fResults(NULL)
 213   ,fSubTaskMap(0)
 214   ,fUseCalib(7)
 215   ,fDet(-1)
 216   ,fCol(-1)
 217   ,fRow(-1)
 218   ,fExB(0.)
 219   ,fDt(0.)
 220   ,fDl(0.)
 221   ,fVdrift(1.5)
 222   ,fT0(0.)
 223   ,fGain(1.)
 224   ,fXch(0.)
 225   ,fZch(0.)
 226   ,fH(0.)
 227   ,fDyRange(1.3)
 228   ,fLy(0)
 229   ,fT(0.)
 230   ,fX(0.)
 231   ,fY(0.)
 232   ,fZ(0.)
 233 {
 234 // Constructor
 235   SetNameTitle("ClErrCalib", "Cluster Error Parameterization");
 236   memset(fR, 0, 4*sizeof(Float_t));
 237   memset(fP, 0, 4*sizeof(Float_t));
 238 }
 239
 240 //_______________________________________________________
 241 AliTRDclusterResolution::AliTRDclusterResolution(const char *name)
 242   : AliTRDrecoTask(name, "Cluster Error Parameterization")
 243   ,fCanvas(NULL)
 244   ,fInfo(NULL)
 245   ,fResults(NULL)
 246   ,fSubTaskMap(0)
 247   ,fUseCalib(7)
 248   ,fDet(-1)
 249   ,fCol(-1)
 250   ,fRow(-1)
 251   ,fExB(0.)
 252   ,fDt(0.)
 253   ,fDl(0.)
 254   ,fVdrift(1.5)
 255   ,fT0(0.)
 256   ,fGain(1.)
 257   ,fXch(0.)
 258   ,fZch(0.)
 259   ,fH(0.)
 260   ,fDyRange(1.3)
 261   ,fLy(0)
 262   ,fT(0.)
 263   ,fX(0.)
 264   ,fY(0.)
 265   ,fZ(0.)
 266 {
 267 // Constructor
 268
 269   memset(fR, 0, 4*sizeof(Float_t));
 270   memset(fP, 0, 4*sizeof(Float_t));
 271
 272   // By default register all analysis
 273   // The user can switch them off in his steering macro
 274   SetProcess(kYRes);
 275   SetProcess(kYSys);
 276   SetProcess(kMean);
 277   SetProcess(kSigm);
 278 }
 279
 280 //_______________________________________________________
 281 AliTRDclusterResolution::~AliTRDclusterResolution()
 282 {
 283 // Destructor
 284
 285   if(fCanvas) delete fCanvas;
 286   if(fResults){
 287     fResults->Delete();
 288     delete fResults;
 289   }
 290 }
 291
 292 //_______________________________________________________
 293 void AliTRDclusterResolution::UserCreateOutputObjects()
 294 {
 295 /*  fContainer = Histos();
 296   PostData(1, fContainer);*/
 297 }
 298
 299 //_______________________________________________________
 300 Bool_t AliTRDclusterResolution::GetRefFigure(Int_t ifig)
 301 {
 302 // Steering function to retrieve performance plots
 303
 304   if(!fResults) return kFALSE;
 305   TLegend *leg = NULL;
 306   TList *l = NULL;
 307   TObjArray *arr = NULL;
 308   TTree *t = NULL;
 309   TH2 *h2 = NULL;TH1 *h1 = NULL;
 310   TGraphErrors *gm(NULL), *gs(NULL), *gp(NULL);
 311   switch(ifig){
 312   case kYRes:
 313     if(!(arr = (TObjArray*)fResults->At(kYRes))) break;
 314     if(!(gm = (TGraphErrors*)arr->At(0))) break;
 315     if(!(gs = (TGraphErrors*)arr->At(1))) break;
 316     if(!(gp = (TGraphErrors*)arr->At(2))) break;
 317     leg= new TLegend(.7, .7, .9, .95);
 318     leg->SetBorderSize(0); leg->SetFillColor(0); leg->SetFillStyle(0);
 319     gs->Draw("apl"); leg->AddEntry(gs, "Sigma / Resolution", "pl");
 320     gs->GetHistogram()->GetYaxis()->SetRangeUser(-50., 700.);
 321     gs->GetHistogram()->SetXTitle("Q [a.u.]");
 322     gs->GetHistogram()->SetYTitle("y - x tg(#alpha_{L}) [#mum]");
 323     gm->Draw("pl");leg->AddEntry(gm, "Mean / Systematics", "pl");
 324     gp->Draw("pl");leg->AddEntry(gp, "Abundance / Probability", "pl");
 325     leg->Draw();
 326     return kTRUE;
 327   case kYSys:
 328     if(!(arr = (TObjArray*)fResults->At(kYSys))) break;
 329     gPad->Divide(2, 1); l = gPad->GetListOfPrimitives();
 330     ((TVirtualPad*)l->At(0))->cd();
 331     ((TTree*)arr->At(0))->Draw(Form("y:t>>h(%d, -0.5, %f, 51, -.51, .51)",AliTRDseedV1::kNtb, AliTRDseedV1::kNtb-0.5),
 332             "m[0]*(ly==0&&abs(m[0])<1.e-1)", "colz");
 333     ((TVirtualPad*)l->At(1))->cd();
 334     leg= new TLegend(.7, .7, .9, .95);
 335     leg->SetBorderSize(0); leg->SetFillColor(0); leg->SetFillStyle(0);
 336     leg->SetHeader("TRD Plane");
 337     for(Int_t il = 1; il<=AliTRDgeometry::kNlayer; il++){
 338       if(!(gm = (TGraphErrors*)arr->At(il))) return kFALSE;
 339       gm->Draw(il>1?"pc":"apc"); leg->AddEntry(gm, Form("%d", il-1), "pl");
 340       if(il>1) continue;
 341       gm->GetHistogram()->SetXTitle("t_{drift} [tb]");
 342       gm->GetHistogram()->SetYTitle("#sigma_{y}(x|cen=0) [#mum]");
 343       gm->GetHistogram()->GetYaxis()->SetRangeUser(150., 500.);
 344     }
 345     leg->Draw();
 346     return kTRUE;
 347   case kSigm:
 348     if(!(t = (TTree*)fResults->At(kSigm))) break;
 349     t->Draw("z:t>>h2x(23, 0.1, 2.4, 25, 0., 2.5)","sx*(1)", "lego2fb");
 350     h2 = (TH2F*)gROOT->FindObject("h2x");
 351     printf("  const Double_t sx[24][25]={\n");
 352     for(Int_t ix=1; ix<=h2->GetNbinsX(); ix++){
 353       printf("    {");
 354       for(Int_t iy=1; iy<h2->GetNbinsY(); iy++){
 355         printf("%6.4f ", h2->GetBinContent(ix, iy));
 356       }
 357       printf("%6.4f},\n", h2->GetBinContent(ix, h2->GetNbinsY()));
 358     }
 359     printf("  };\n");
 360     gPad->Divide(2, 1, 1.e-5, 1.e-5); l = gPad->GetListOfPrimitives();
 361     ((TVirtualPad*)l->At(0))->cd();
 362     h1 = h2->ProjectionX("hsx_pxx"); h1->Scale(1.e4/kND); h1->SetMarkerStyle(24);
 363     h1->SetYTitle("<#sigma_{x}> [#mum]");
 364     h1->SetXTitle("t_{drift} [#mus]");
 365     h1->GetXaxis()->SetRange(2, AliTRDseedV1::kNtb-1); h1->Draw("pc");
 366
 367     t->Draw("z:t>>h2y(23, 0.1, 2.4, 25, 0., 2.5)","sy*(1)", "lego2fb");
 368     h2 = (TH2F*)gROOT->FindObject("h2y");
 369     printf("  const Double_t sy[24][25]={\n");
 370     for(Int_t ix=1; ix<=h2->GetNbinsX(); ix++){
 371       printf("    {");
 372       for(Int_t iy=1; iy<h2->GetNbinsY(); iy++){
 373         printf("%6.4f ", h2->GetBinContent(ix, iy));
 374       }
 375       printf("%6.4f},\n", h2->GetBinContent(ix, h2->GetNbinsY()));
 376     }
 377     printf("  };\n");
 378     ((TVirtualPad*)l->At(1))->cd();
 379     h1 = h2->ProjectionX("hsy_pxx"); h1->Scale(1.e4/kND); h1->SetMarkerStyle(24);
 380     h1->SetYTitle("<#sigma_{y}> [#mum]");
 381     h1->SetXTitle("t_{drift} [#mus]");
 382     h1->GetXaxis()->SetRange(2, AliTRDseedV1::kNtb-1); h1->Draw("pc");
 383     return kTRUE;
 384   case kMean:
 385     if(!(t = (TTree*)fResults->At(kMean))) break;
 386     if(!t->Draw(Form("z:t>>h2x(%d, -0.5, %3.1f, %d, 0., 2.5)",
 387       AliTRDseedV1::kNtb, AliTRDseedV1::kNtb-0.5, kND),
 388       "dx*(1)", "goff")) break;
 389     h2 = (TH2F*)gROOT->FindObject("h2x");
 390     printf("  const Double_t dx[%d][%d]={\n", AliTRDseedV1::kNtb, kND);
 391     for(Int_t ix=1; ix<=h2->GetNbinsX(); ix++){
 392       printf("    {");
 393       for(Int_t iy=1; iy<h2->GetNbinsY(); iy++){
 394         printf("%+6.4e, ", h2->GetBinContent(ix, iy));
 395       }
 396       printf("%+6.4e},\n", h2->GetBinContent(ix, h2->GetNbinsY()));
 397     }
 398     printf("  };\n");
 399     gPad->Divide(2, 2, 1.e-5, 1.e-5); l = gPad->GetListOfPrimitives();
 400     ((TVirtualPad*)l->At(0))->cd();
 401     h2->Draw("lego2fb");
 402     ((TVirtualPad*)l->At(2))->cd();
 403     h1 = h2->ProjectionX("hdx_pxx"); h1->Scale(1.e4/kND); h1->SetMarkerStyle(24);
 404     h1->SetYTitle("<#deltax> [#mum]");
 405     h1->SetXTitle("t_{drift} [tb]");
 406     //h1->GetXaxis()->SetRange(2, AliTRDseedV1::kNtb-1);
 407     h1->Draw("pc");
 408
 409     if(!t->Draw(Form("z:t>>h2y(%d, -0.5, %3.1f, %d, 0., 2.5)",
 410       AliTRDseedV1::kNtb, AliTRDseedV1::kNtb-0.5, kND),
 411       "dy*(1)", "goff")) break;
 412     h2 = (TH2F*)gROOT->FindObject("h2y");
 413     printf("  const Double_t dy[%d][%d]={\n", AliTRDseedV1::kNtb, kND);
 414     for(Int_t ix=1; ix<=h2->GetNbinsX(); ix++){
 415       printf("    {");
 416       for(Int_t iy=1; iy<h2->GetNbinsY(); iy++){
 417         printf("%+6.4e ", h2->GetBinContent(ix, iy));
 418       }
 419       printf("%+6.4e},\n", h2->GetBinContent(ix, h2->GetNbinsY()));
 420     }
 421     printf("  };\n");
 422     ((TVirtualPad*)l->At(1))->cd();
 423     h2->Draw("lego2fb");
 424     ((TVirtualPad*)l->At(3))->cd();
 425     h1 = h2->ProjectionX("hdy_pxx"); h1->Scale(1.e4/kND); h1->SetMarkerStyle(24);
 426     h1->SetYTitle("<#deltay> [#mum]");
 427     h1->SetXTitle("t_{drift} [tb]");
 428     //h1->GetXaxis()->SetRange(2, AliTRDseedV1::kNtb-1);
 429     h1->Draw("pc");
 430
 431     return kTRUE;
 432   default:
 433     break;
 434   }
 435   AliWarning("No container/data found.");
 436   return kFALSE;
 437 }
 438
 439 //_______________________________________________________
 440 TObjArray* AliTRDclusterResolution::Histos()
 441 {
 442 // Retrieve histograms array if already build or build it
 443
 444   if(fContainer) return fContainer;
 445   fContainer = new TObjArray(kNtasks);
 446   //fContainer->SetOwner(kTRUE);
 447
 448   TH3S *h3(NULL);TH2I *h2(NULL);
 449   TObjArray *arr(NULL);
 450   if(!HasGlobalPosition() && !LoadGlobalChamberPosition()) return NULL;
 451   Float_t tgt(fZch/fXch), htgt(fH*tgt);
 452
 453   // SYSTEMATIC PLOTS
 454   fContainer->AddAt(arr = new TObjArray(3), kYSys);
 455   arr->SetName("SysY");
 456   // systematic plot on pw and q (dydx=ExB+h*dzdx)
 457   if(!(h3=(TH3S*)gROOT->FindObject(Form("Sys%s%03d", (HasMCdata()?"MC":"") ,fDet)))) {
 458     h3 = new TH3S(
 459       Form("Sys%s%03d", (HasMCdata()?"MC":""),fDet),
 460       Form(" Det[%d] Col[%d] Row[%d];log q [a.u.];#deltay [pw];#Delta y[cm]", fDet, fCol, fRow),
 461       45,   2., 6.5,            // log(q) [a.u.]
 462       25, -.51, .51,            // y [pw]
 463       60, -fDyRange, fDyRange); // dy [cm]
 464   } h3->Reset();
 465   arr->AddAt(h3, 0);
 466   // systematic plot on tb (only for dydx = h*tgt + exb and MPV q)
 467   if(!(h2 = (TH2I*)gROOT->FindObject(Form("SysTb%s%03d", (HasMCdata()?"MC":""), fDet)))){
 468     h2 = new TH2I(Form("SysTb%s%03d", (HasMCdata()?"MC":""), fDet),
 469     Form(" Det[%d] Col[%d] Row[%d];t [time bin];#Delta y[cm]", fDet, fCol, fRow),
 470     AliTRDseedV1::kNtb, -.5, AliTRDseedV1::kNtb-0.5,   // t [tb]
 471     60, -fDyRange, fDyRange);                          // dy [cm]
 472   } h2->Reset();
 473   arr->AddAt(h2, 1);
 474   // systematic plot on tgp and tb (for MPV q)
 475   if(!(h3=(TH3S*)gROOT->FindObject(Form("SysTbTgp%s%03d", (HasMCdata()?"MC":""), fDet)))){
 476     h3 = new TH3S(
 477       Form("SysTbTgp%s%03d", (HasMCdata()?"MC":""), fDet),
 478       Form(" Det[%d];t [time bin];tg(#phi) - h*tg(#theta) %s;#Delta y[cm]", fDet, fExB>1.e-5?"- tg(#alpha_{L})":""),
 479       AliTRDseedV1::kNtb, -.5, AliTRDseedV1::kNtb-0.5,   // t [tb]
 480       36, fExB-.18, fExB+.18,                            // tgp-h tgt-tg(aL)
 481       60, -fDyRange, fDyRange);                          // dy
 482   } h3->Reset();
 483   arr->AddAt(h3, 2);
 484
 485   //  RESOLUTION/PULLS PLOTS
 486   fContainer->AddAt(arr = new TObjArray(6), kYRes);
 487   arr->SetName("ResY");
 488   // resolution plot on pw and q (for dydx=0 && B=0)
 489   if(!(h3=(TH3S*)gROOT->FindObject(Form("Res%s%03d", (HasMCdata()?"MC":"") ,fDet)))) {
 490     h3 = new TH3S(
 491       Form("Res%s%03d", (HasMCdata()?"MC":""),fDet),
 492       Form(" Det[%d] Col[%d] Row[%d];log q [a.u];#deltay [pw];#Delta y[cm]", fDet, fCol, fRow),
 493       45,   2., 6.5,            // log(q) [a.u]
 494       25, -.51, .51,            // y [pw]
 495       60, -fDyRange, fDyRange); // dy
 496   } h3->Reset();
 497   arr->AddAt(h3, 0);
 498   // Pull plot on pw and q (for dydx=0 && B=0)
 499   if(!(h3=(TH3S*)gROOT->FindObject(Form("Pull%s%03d", (HasMCdata()?"MC":""), fDet)))){
 500     h3 = new TH3S(
 501       Form("Pull%s%03d", (HasMCdata()?"MC":""), fDet),
 502       Form(" Det[%d] Col[%d] Row[%d];log q [a.u.];#deltay [pw];#Delta y[cm]/#sigma_{y}", fDet, fCol, fRow),
 503       45,   2., 6.5,            // log(q) [a.u]
 504       25, -.51, .51,            // y [pw]
 505       60, -4., 4.);             // dy/sy
 506   } h3->Reset();
 507   arr->AddAt(h3, 1);
 508   // resolution/pull plot on tb (for dydx=0 && B=0 && MPV q)
 509   if(!(h3 = (TH3S*)gROOT->FindObject(Form("ResPullTb%s%03d", (HasMCdata()?"MC":""), fDet)))){
 510     h3 = new TH3S(Form("ResPullTb%s%03d", (HasMCdata()?"MC":""), fDet),
 511     Form(" Det[%d] Col[%d] Row[%d];t [time bin];#Delta y[cm];#Delta y/#sigma_{y}", fDet, fCol, fRow),
 512     AliTRDseedV1::kNtb, -.5, AliTRDseedV1::kNtb-0.5,   // t [tb]
 513     60, -fDyRange, fDyRange,                           // dy [cm]
 514     60, -4., 4.);                                      // dy/sy
 515   } h3->Reset();
 516   arr->AddAt(h3, 2);
 517   // resolution plot on pw and q (for dydx=0 && B=0) np = 2
 518   if(!(h3=(TH3S*)gROOT->FindObject(Form("Res2%s%03d", (HasMCdata()?"MC":"") ,fDet)))) {
 519     h3 = new TH3S(
 520       Form("Res2%s%03d", (HasMCdata()?"MC":""),fDet),
 521       Form(" Det[%d] Col[%d] Row[%d];log q [a.u];#deltay [pw];#Delta y[cm]", fDet, fCol, fRow),
 522       45,   2., 6.5,            // log(q) [a.u]
 523       25, -.51, .51,            // y [pw]
 524       60, -fDyRange, fDyRange); // dy
 525   } h3->Reset();
 526   arr->AddAt(h3, 3);
 527   // resolution plot on pw and q (for dydx=0 && B=0) np = 4
 528   if(!(h3=(TH3S*)gROOT->FindObject(Form("Res4%s%03d", (HasMCdata()?"MC":"") ,fDet)))) {
 529     h3 = new TH3S(
 530       Form("Res4%s%03d", (HasMCdata()?"MC":""),fDet),
 531       Form(" Det[%d] Col[%d] Row[%d];log q [a.u];#deltay [pw];#Delta y[cm]", fDet, fCol, fRow),
 532       45,   2., 6.5,            // log(q) [a.u]
 533       25, -.51, .51,            // y [pw]
 534       60, -fDyRange, fDyRange); // dy
 535   } h3->Reset();
 536   arr->AddAt(h3, 4);
 537   // systemtic plot of tb on pw and q (for dydx=0 && B=0)
 538   if(!(h3=(TH3S*)gROOT->FindObject(Form("SysTbPwQ%s%03d", (HasMCdata()?"MC":"") ,fDet)))) {
 539     h3 = new TH3S(
 540       Form("SysTbPwQ%s%03d", (HasMCdata()?"MC":""),fDet),
 541       Form(" Det[%d] Col[%d] Row[%d];log q [a.u];#deltay [pw];t [time bin]", fDet, fCol, fRow),
 542       45,   2., 6.5,            // log(q) [a.u]
 543       25, -.51, .51,            // y [pw]
 544       AliTRDseedV1::kNtb, -.5, AliTRDseedV1::kNtb-0.5);   // t [tb]
 545   } h3->Reset();
 546   arr->AddAt(h3, 5);
 547
 548
 549
 550   fContainer->AddAt(arr = new TObjArray(AliTRDseedV1::kNtb), kSigm);
 551   arr->SetName("Resolution");
 552   for(Int_t it=0; it<AliTRDseedV1::kNtb; it++){
 553     if(!(h3=(TH3S*)gROOT->FindObject(Form("hr%s%03d_t%02d", (HasMCdata()?"MC":""), fDet, it)))){
 554       h3 = new TH3S(
 555         Form("hr%s%03d_t%02d", (HasMCdata()?"MC":""), fDet, it),
 556         Form(" Det[%d] t_{drift}(%2d)[bin];h*tg(#theta);tg(#phi);#Delta y[cm]", fDet, it),
 557         35, htgt-0.0035, htgt+0.0035, // h*tgt
 558         36, fExB-.18, fExB+.18,       // tgp
 559         60, -fDyRange, fDyRange);     // dy
 560     } h3->Reset();
 561     arr->AddAt(h3, it);
 562   }
 563   return fContainer;
 564 }
 565
 566 //_______________________________________________________
 567 void AliTRDclusterResolution::UserExec(Option_t *)
 568 {
 569 // Fill container histograms
 570
 571   if(!IsCalibrated()){
 572     LoadCalibration();
 573     if(!IsCalibrated()){
 574       AliFatal("Loading the calibration settings failed. Check OCDB access.");
 575       return;
 576     }
 577   }
 578   if(!fContainer){
 579     fContainer = Histos();
 580     PostData(1, fContainer);
 581   }
 582   fInfo = dynamic_cast<TObjArray *>(GetInputData(1));
 583   AliDebug(2, Form("Clusters[%d]", fInfo->GetEntriesFast()));
 584
 585   Int_t det, t, np;
 586   Float_t x, y, z, q, dy, dydx, dzdx, cov[3], covcl[3];
 587   TH3S *h3(NULL); TH2I *h2(NULL);
 588
 589   // define limits around ExB for which x contribution is negligible
 590   const Float_t kAroundZero = 3.5e-2; //(+- 2 deg)
 591
 592   TObjArray *arr0 = (TObjArray*)fContainer->At(kYSys);
 593   TObjArray *arr1 = (TObjArray*)fContainer->At(kYRes);
 594   TObjArray *arr2 = (TObjArray*)fContainer->At(kSigm);
 595
 596   const AliTRDclusterInfo *cli = NULL;
 597   TIterator *iter=fInfo->MakeIterator();
 598   while((cli=dynamic_cast<AliTRDclusterInfo*>((*iter)()))){
 599     if((np = cli->GetNpads())>4) continue;
 600     cli->GetCluster(det, x, y, z, q, t, covcl);
 601
 602     // select cluster according to detector region if specified
 603     if(fDet>=0 && fDet!=det) continue;
 604     if(fCol>=0 && fRow>=0){
 605       Int_t c,r;
 606       cli->GetCenterPad(c, r);
 607       if(TMath::Abs(fCol-c) > 5) continue;
 608       if(TMath::Abs(fRow-r) > 2) continue;
 609     }
 610     dy = cli->GetResolution();
 611     AliDebug(4, Form("det[%d] tb[%2d] q[%4.0f Log[%6.4f]] np[%d] dy[%7.2f][um] ypull[%5.2f]", det, t, q, TMath::Log(q), np, 1.e4*dy, dy/TMath::Sqrt(covcl[0])));
 612
 613     cli->GetGlobalPosition(y, z, dydx, dzdx, &cov[0]);
 614     Float_t pw(cli->GetYDisplacement());
 615
 616     // systematics as a function of pw and log(q)
 617     // only for dydx = exB + h*dzdx
 618     if(TMath::Abs(dydx-fExB-fH*dzdx) < kAroundZero){
 619       h3 = (TH3S*)arr0->At(0);
 620       h3->Fill(TMath::Log(q), pw, dy);
 621     }
 622     // resolution/pull as a function of pw and log(q)
 623     // only for dydx = 0, ExB=0
 624     if(TMath::Abs(fExB) < kAroundZero &&
 625        TMath::Abs(dydx) < kAroundZero){
 626       switch(np){
 627       case 3: // MPV np
 628         h3 = (TH3S*)arr1->At(0);
 629         h3->Fill(TMath::Log(q), pw, dy);
 630         h3 = (TH3S*)arr1->At(5);
 631         h3->Fill(TMath::Log(q), pw, t);
 632         break;
 633       case 2: // Min np
 634         h3 = (TH3S*)arr1->At(3);
 635         h3->Fill(TMath::Log(q), pw, dy);
 636         break;
 637       case 4: // Max np
 638         h3 = (TH3S*)arr1->At(4);
 639         h3->Fill(TMath::Log(q), pw, dy);
 640         break;
 641       }
 642       h3 = (TH3S*)arr1->At(1);
 643       h3->Fill(TMath::Log(q), pw, dy/TMath::Sqrt(covcl[0]));
 644     }
 645
 646     // do not use problematic clusters in resolution analysis
 647     // TODO define limits as calibration aware (gain) !!
 648     //if(!AcceptableGain(fGain)) continue;
 649     if(q<20. || q>250.) continue;
 650
 651     // systematic as a function of time bin
 652     // only for dydx = exB + h*dzdx and MPV q
 653     if(TMath::Abs(dydx-fExB-fH*dzdx) < kAroundZero){
 654       h2 = (TH2I*)arr0->At(1);
 655       h2->Fill(t, dy);
 656     }
 657     // systematic as function of tb and tgp
 658     // only for MPV q
 659     h3 = (TH3S*)arr0->At(2);
 660     h3->Fill(t, dydx, dy);
 661
 662     // resolution/pull as a function of time bin
 663     // only for dydx = 0, ExB=0 and MPV q
 664     if(TMath::Abs(fExB) < kAroundZero &&
 665        TMath::Abs(dydx) < kAroundZero){
 666       h3 = (TH3S*)arr1->At(2);
 667       h3->Fill(t, dy, dy/TMath::Sqrt(covcl[0]));
 668     }
 669
 670     // resolution as function of tb, tgp and h*tgt
 671     // only for MPV q
 672     ((TH3S*)arr2->At(t))->Fill(fH*dzdx, dydx, dy);
 673   }
 674 }
 675
 676
 677 //_______________________________________________________
 678 Bool_t AliTRDclusterResolution::PostProcess()
 679 {
 680 // Steer processing of various cluster resolution dependences :
 681 //
 682 //   - process resolution dependency cluster charge
 683 //   if(HasProcess(kYRes)) ProcessCharge();
 684 //   - process resolution dependency on y displacement
 685 //   if(HasProcess(kYSys)) ProcessCenterPad();
 686 //   - process resolution dependency on drift legth and drift cell width
 687 //   if(HasProcess(kSigm)) ProcessSigma();
 688 //   - process systematic shift on drift legth and drift cell width
 689 //   if(HasProcess(kMean)) ProcessMean();
 690
 691   if(!fContainer) return kFALSE;
 692   if(!IsCalibrated()){
 693     AliError("Not calibrated instance.");
 694     return kFALSE;
 695   }
 696   TObjArray *arr = NULL;
 697   TTree *t=NULL;
 698   if(!fResults){
 699     TGraphErrors *g = NULL;
 700     fResults = new TObjArray(kNtasks);
 701     fResults->SetOwner();
 702     fResults->AddAt(arr = new TObjArray(3), kYRes);
 703     arr->SetOwner();
 704     arr->AddAt(g = new TGraphErrors(), 0);
 705     g->SetLineColor(kBlue); g->SetMarkerColor(kBlue);
 706     g->SetMarkerStyle(7);
 707     arr->AddAt(g = new TGraphErrors(), 1);
 708     g->SetLineColor(kRed); g->SetMarkerColor(kRed);
 709     g->SetMarkerStyle(23);
 710     arr->AddAt(g = new TGraphErrors(), 2);
 711     g->SetLineColor(kGreen); g->SetMarkerColor(kGreen);
 712     g->SetMarkerStyle(7);
 713
 714     // pad center dependence
 715     fResults->AddAt(arr = new TObjArray(AliTRDgeometry::kNlayer+1), kYSys);
 716     arr->SetOwner();
 717     arr->AddAt(
 718     t = new TTree("cent", "dy=f(y,x,ly)"), 0);
 719     t->Branch("ly", &fLy, "ly/B");
 720     t->Branch("t", &fT, "t/F");
 721     t->Branch("y", &fY, "y/F");
 722     t->Branch("m", &fR[0], "m[2]/F");
 723     t->Branch("s", &fR[2], "s[2]/F");
 724     t->Branch("pm", &fP[0], "pm[2]/F");
 725     t->Branch("ps", &fP[2], "ps[2]/F");
 726     for(Int_t il=1; il<=AliTRDgeometry::kNlayer; il++){
 727       arr->AddAt(g = new TGraphErrors(), il);
 728       g->SetLineColor(il); g->SetLineStyle(il);
 729       g->SetMarkerColor(il);g->SetMarkerStyle(4);
 730     }
 731
 732
 733     fResults->AddAt(t = new TTree("sigm", "dy=f(dw,x,dydx)"), kSigm);
 734     t->Branch("t", &fT, "t/F");
 735     t->Branch("x", &fX, "x/F");
 736     t->Branch("z", &fZ, "z/F");
 737     t->Branch("sx", &fR[0], "sx[2]/F");
 738     t->Branch("sy", &fR[2], "sy[2]/F");
 739
 740
 741     fResults->AddAt(t = new TTree("mean", "dy=f(dw,x,dydx - h dzdx)"), kMean);
 742     t->Branch("t", &fT, "t/F");
 743     t->Branch("x", &fX, "x/F");
 744     t->Branch("z", &fZ, "z/F");
 745     t->Branch("dx", &fR[0], "dx[2]/F");
 746     t->Branch("dy", &fR[2], "dy[2]/F");
 747   } else {
 748     TObject *o = NULL;
 749     TIterator *iter=fResults->MakeIterator();
 750     while((o=((*iter)()))) o->Clear(); // maybe it is wrong but we should never reach this point
 751   }
 752
 753   // process resolution dependency on charge
 754   if(HasProcess(kYRes)) ProcessCharge();
 755
 756   // process resolution dependency on y displacement
 757   if(HasProcess(kYSys)) ProcessNormalTracks();
 758
 759   // process resolution dependency on drift legth and drift cell width
 760   if(HasProcess(kSigm)) ProcessSigma();
 761
 762   // process systematic shift on drift legth and drift cell width
 763   if(HasProcess(kMean)) ProcessMean();
 764
 765   return kTRUE;
 766 }
 767
 768 //_______________________________________________________
 769 Bool_t AliTRDclusterResolution::LoadCalibration()
 770 {
 771 // Retrieve calibration parameters from OCDB, drift velocity and t0 for the detector region specified by
 772 // a previous call to AliTRDclusterResolution::SetCalibrationRegion().
 773
 774   AliCDBManager *cdb = AliCDBManager::Instance(); // check access OCDB
 775   if(cdb->GetRun() < 0){
 776     AliError("OCDB manager not properly initialized");
 777     return kFALSE;
 778   }
 779   // check magnetic field
 780   if(!TGeoGlobalMagField::Instance() || !TGeoGlobalMagField::Instance()->IsLocked()){
 781     AliError("Magnetic field not available.");
 782     return kFALSE;
 783   }
 784
 785   AliTRDcalibDB *fCalibration  = AliTRDcalibDB::Instance();
 786   AliTRDCalROC  *fCalVdriftROC(fCalibration->GetVdriftROC(fDet>=0?fDet:0))
 787                ,*fCalT0ROC(fCalibration->GetT0ROC(fDet>=0?fDet:0));
 788   const AliTRDCalDet  *fCalVdriftDet = fCalibration->GetVdriftDet();
 789   const AliTRDCalDet  *fCalT0Det = fCalibration->GetT0Det();
 790
 791   if(IsUsingCalibParam(kVdrift)){
 792     fVdrift = fCalVdriftDet->GetValue(fDet>=0?fDet:0);
 793     if(fCol>=0 && fRow>=0) fVdrift*= fCalVdriftROC->GetValue(fCol, fRow);
 794   }
 795   fExB    = AliTRDCommonParam::Instance()->GetOmegaTau(fVdrift);
 796   AliTRDCommonParam::Instance()->GetDiffCoeff(fDt, fDl, fVdrift);
 797   if(IsUsingCalibParam(kT0)){
 798     fT0     = fCalT0Det->GetValue(fDet>=0?fDet:0);
 799     if(fCol>=0 && fRow>=0) fT0 *= fCalT0ROC->GetValue(fCol, fRow);
 800   }
 801   if(IsUsingCalibParam(kGain)) fGain = (fCol>=0 && fRow>=0)?fCalibration-> GetGainFactor(fDet, fCol, fRow):fCalibration-> GetGainFactorAverage(fDet);
 802
 803   SetBit(kCalibrated);
 804
 805   AliInfo(Form("Calibration D[%3d] Col[%3d] Row[%2d] : \n   t0[%5.3f] vd[%5.3f] gain[%5.3f] ExB[%f] DiffT[%f] DiffL[%f]", fDet, fCol, fRow, fT0, fVdrift, fGain, fExB, fDt, fDl));
 806
 807   return kTRUE;
 808 }
 809
 810 //_______________________________________________________
 811 Bool_t AliTRDclusterResolution::LoadGlobalChamberPosition()
 812 {
 813 // Retrieve global chamber position from alignment
 814 // a previous call to AliTRDclusterResolution::SetCalibrationRegion() is mandatory.
 815
 816   TGeoHMatrix *matrix(NULL);
 817   Double_t loc[] = {0., 0., 0.}, glb[] = {0., 0., 0.};
 818   AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
 819   if(!(matrix= geo->GetClusterMatrix(fDet))) {
 820     AliFatal(Form("Could not get transformation matrix for detector %d.", fDet));
 821     return kFALSE;
 822   }
 823   matrix->LocalToMaster(loc, glb);
 824   fXch = glb[0]+ AliTRDgeometry::AnodePos()-.5*AliTRDgeometry::AmThick() - AliTRDgeometry::DrThick();
 825   AliTRDpadPlane *pp(geo->GetPadPlane(fDet));
 826   fH = TMath::Tan(pp->GetTiltingAngle()*TMath::DegToRad());
 827   fZch = 0.;
 828   if(fRow>=0){
 829     fZch = pp->GetRowPos(fRow)+0.5*pp->GetLengthIPad();
 830   }else{
 831     Int_t  nrows(pp->GetNrows());
 832     Float_t zmax(pp->GetRow0()),
 833             zmin(zmax - 2 * pp->GetLengthOPad()
 834                 - (nrows-2) * pp->GetLengthIPad()
 835                 - (nrows-1) * pp->GetRowSpacing());
 836     fZch = 0.5*(zmin+zmax);
 837   }
 838
 839   AliInfo(Form("Global pos. D[%3d] Col[%3d] Row[%2d] : \n   x[%6.2f] z[%6.2f] h[%+6.2f].", fDet, fCol, fRow, fXch, fZch, fH));
 840   SetBit(kGlobal);
 841   return kTRUE;
 842 }
 843
 844 //_______________________________________________________
 845 void AliTRDclusterResolution::SetCalibrationRegion(Int_t det, Int_t col, Int_t row)
 846 {
 847 // Set calibration region in terms of detector and pad.
 848 // By default detector 0 mean values are considered.
 849
 850   if(det>=0 && det<AliTRDgeometry::kNdet){
 851     fDet = det;
 852     if(col>=0 && row>=0){
 853       // check pad col/row for detector
 854       AliTRDgeometry *geo = AliTRDinfoGen::Geometry();
 855       AliTRDpadPlane *pp(geo->GetPadPlane(fDet));
 856       if(fCol>=pp->GetNcols() ||
 857         fRow>=pp->GetNrows()){
 858         AliWarning(Form("Pad coordinates col[%d] or row[%d] incorrect for det[%d].\nLimits are max col[%d] max row[%d]. Reset to default", fCol, fRow, fDet, pp->GetNcols(), pp->GetNrows()));
 859         fCol = -1; fRow=-1;
 860       } else {
 861         fCol = col;
 862         fRow = row;
 863       }
 864     }
 865   } else {
 866     AliFatal(Form("Detector index outside range [0 %d].", AliTRDgeometry::kNdet));
 867   }
 868   return;
 869 }
 870
 871 //_______________________________________________________
 872 void AliTRDclusterResolution::SetVisual()
 873 {
 874   if(fCanvas) return;
 875   fCanvas = new TCanvas("clResCanvas", "Cluster Resolution Visualization", 10, 10, 600, 600);
 876 }
 877
 878 //_______________________________________________________
 879 void AliTRDclusterResolution::ProcessCharge()
 880 {
 881 // Resolution as a function of cluster charge.
 882 //
 883 // As described in the function ProcessCenterPad() the error parameterization for clusters for phi = a_L can be
 884 // written as:
 885 // BEGIN_LATEX
 886 // #sigma_{y}^{2} = #sigma_{y}^{2}|_{B=0} + tg^{2}(#alpha_{L})*#sigma_{x}^{2}
 887 // END_LATEX
 888 // with the contribution in case of B=0 given by:
 889 // BEGIN_LATEX
 890 // #sigma_{y}|_{B=0} = #sigma_{diff}*Gauss(0, s_{ly}) + #delta_{#sigma}(q)
 891 // END_LATEX
 892 // which further can be simplified to:
 893 // BEGIN_LATEX
 894 // <#sigma_{y}|_{B=0}>(q) = <#sigma_{y}> + #delta_{#sigma}(q)
 895 // <#sigma_{y}> = #int{f(q)#sigma_{y}dq}
 896 // END_LATEX
 897 // The results for s_y and f(q) are displayed below:
 898 //Begin_Html
 899 //<img src="TRD/clusterQerror.gif">
 900 //End_Html
 901 // The function has to extended to accomodate gain calibration scalling and errors.
 902 //
 903 // Author
 904 // Alexandru Bercuci <A.Bercuci@gsi.de>
 905
 906
 907
 908   TObjArray *arr(NULL);
 909   if(!(arr = (TObjArray*)fContainer->At(kYSys))) {
 910     AliError("Missing systematic container");
 911     return;
 912   }
 913   TH3S *h3s(NULL);
 914   if(!(h3s = (TH3S*)arr->At(0))){
 915     AliError("Missing systematic histo");
 916     return;
 917   }
 918   // PROCESS SYSTEMATIC
 919   Float_t tmin(6.5), tmax(20.5), tmed(0.5*(tmin+tmax));
 920   TGraphErrors *g[2]; TH1 *h(NULL);
 921   g[0] = new TGraphErrors();
 922   g[0]->SetMarkerStyle(24);g[0]->SetMarkerColor(kBlue);g[0]->SetLineColor(kBlue);
 923   g[1] = new TGraphErrors();
 924   g[1]->SetMarkerStyle(24);g[1]->SetMarkerColor(kRed);g[1]->SetLineColor(kRed);
 925   // define model for systematic shift vs pw
 926   TF1 fm("fm", "[0]+[1]*sin(x*[2])", -.45,.45);
 927   fm.SetParameter(0, 0.); fm.SetParameter(1, 1.e-2); fm.FixParameter(2, TMath::TwoPi());
 928   fm.SetParNames("#deltay", "#pm#delta", "2*#pi");
 929   h3s->GetXaxis()->SetRange(tmin, tmax);
 930   if(!AliTRDresolution::Process((TH2*)h3s->Project3D("zy"), g))return;
 931   g[0]->Fit(&fm, "QR");
 932   if(fCanvas){
 933     g[0]->Draw("apl");
 934     fCanvas->Modified(); fCanvas->Update();
 935     h = g[0]->GetHistogram();
 936     h->SetTitle(fm.GetTitle());
 937     h->GetXaxis()->SetTitle("pw");h->GetXaxis()->CenterTitle();
 938     h->GetYaxis()->SetTitle("#Delta y[cm]");h->GetYaxis()->CenterTitle();
 939     if(IsSaveAs()) fCanvas->SaveAs(Form("D%03d_SysNormTrack_pw.gif", fDet));
 940     else gSystem->Sleep(100);
 941   }
 942
 943   // define model for systematic shift vs tb
 944   TF1 fx("fx", "[0]+0.1*[1]*(x-[2])", tmin, tmax);
 945   fx.SetParNames("#deltay", "#deltay/t", "<t>");
 946   fx.FixParameter(2, tmed);
 947   h3s->GetXaxis()->UnZoom();
 948   if(!AliTRDresolution::Process((TH2*)h3s->Project3D("zx"), g)) return;
 949   g[0]->Fit(&fx, "Q", "", tmin, tmax);
 950   if(fCanvas){
 951     g[0]->Draw("apl");
 952     fCanvas->Modified(); fCanvas->Update();
 953     h = g[0]->GetHistogram();
 954     h->SetTitle(fx.GetTitle());
 955     h->GetXaxis()->SetTitle("t [tb]");h->GetXaxis()->CenterTitle();
 956     h->GetYaxis()->SetTitle("#Delta y[cm]");h->GetYaxis()->CenterTitle();
 957     if(IsSaveAs()) fCanvas->SaveAs(Form("D%03d_SysNormTrack_tb.gif", fDet));
 958     else gSystem->Sleep(100);
 959   }
 960
 961   TH3S *h3(NULL);
 962   if(!(h3 = (TH3S*)fContainer->At(kYRes))) {
 963     AliWarning("Missing dy=f(Q) histo");
 964     return;
 965   }
 966   TF1 f("f", "gaus", -.5, .5);
 967   TAxis *ax(NULL);
 968   TH1 *h1(NULL);
 969
 970   // compute mean error on x
 971   Double_t s2x = 0.;
 972   for(Int_t ix=5; ix<AliTRDseedV1::kNtb; ix++){
 973     // retrieve error on the drift length
 974     s2x += AliTRDcluster::GetSX(ix);
 975   }
 976   s2x /= (AliTRDseedV1::kNtb-5); s2x *= s2x;
 977   //Double_t exb2 = fExB*fExB;
 978
 979   arr = (TObjArray*)fResults->At(kYRes);
 980   TGraphErrors *gqm = (TGraphErrors*)arr->At(0);
 981   TGraphErrors *gqs = (TGraphErrors*)arr->At(1);
 982   TGraphErrors *gqp = (TGraphErrors*)arr->At(2);
 983   Double_t q, n = 0., entries;
 984   ax = h3->GetXaxis();
 985   for(Int_t ix=1; ix<=ax->GetNbins(); ix++){
 986     q = TMath::Exp(ax->GetBinCenter(ix));
 987     ax->SetRange(ix, ix);
 988     h1 = h3->Project3D("y");
 989     entries = h1->GetEntries();
 990     if(entries < 150) continue;
 991     h1->Fit(&f, "Q");
 992
 993     // Fill sy^2 = f(q)
 994     Int_t ip = gqm->GetN();
 995     gqm->SetPoint(ip, q, 1.e4*f.GetParameter(1));
 996     gqm->SetPointError(ip, 0., 1.e4*f.GetParError(1));
 997
 998     // correct sigma for ExB effect
 999     gqs->SetPoint(ip, q, 1.e4*f.GetParameter(2)/**f.GetParameter(2)-exb2*s2x)*/);
1000     gqs->SetPointError(ip, 0., 1.e4*f.GetParError(2)/**f.GetParameter(2)*/);
1001
1002     // save probability
1003     n += entries;
1004     gqp->SetPoint(ip, q, entries);
1005     gqp->SetPointError(ip, 0., 0./*TMath::Sqrt(entries)*/);
1006   }
1007
1008   // normalize probability and get mean sy
1009   Double_t sm = 0., sy;
1010   for(Int_t ip=gqp->GetN(); ip--;){
1011     gqp->GetPoint(ip, q, entries);
1012     entries/=n;
1013     gqp->SetPoint(ip, q, 1.e4*entries);
1014     gqs->GetPoint(ip, q, sy);
1015     sm += entries*sy;
1016   }
1017
1018   // error parametrization s(q) = <sy> + b(1/q-1/q0)
1019   TF1 fq("fq", "[0] + [1]/x", 20., 250.);
1020   gqs->Fit(&fq/*, "W"*/);
1021   printf("sm=%f [0]=%f [1]=%f\n", 1.e-4*sm, fq.GetParameter(0), fq.GetParameter(1));
1022   printf("  const Float_t sq0inv = %f; // [1/q0]\n", (sm-fq.GetParameter(0))/fq.GetParameter(1));
1023   printf("  const Float_t sqb    = %f; // [cm]\n", 1.e-4*fq.GetParameter(1));
1024 }
1025
1026 //_______________________________________________________
1027 Bool_t AliTRDclusterResolution::ProcessNormalTracks()
1028 {
1029 // Resolution as a function of y displacement from pad center and drift length.
1030 //
1031 // Since the error parameterization of cluster r-phi position can be written as (see AliTRDcluster::SetSigmaY2()):
1032 // BEGIN_LATEX
1033 // #sigma_{y}^{2} = (#sigma_{diff}*Gauss(0, s_{ly}) + #delta_{#sigma}(q))^{2} + tg^{2}(#alpha_{L})*#sigma_{x}^{2} + tg^{2}(#phi-#alpha_{L})*#sigma_{x}^{2}+[tg(#phi-#alpha_{L})*tg(#alpha_{L})*x]^{2}/12
1034 // END_LATEX
1035 // one can see that for phi = a_L one gets the following expression:
1036 // BEGIN_LATEX
1037 // #sigma_{y}^{2} = #sigma_{y}^{2}|_{B=0} + tg^{2}(#alpha_{L})*#sigma_{x}^{2}
1038 // END_LATEX
1039 // where we have explicitely marked the remaining term in case of absence of magnetic field. Thus one can use the
1040 // previous equation to estimate s_y for B=0 and than by comparing in magnetic field conditions one can get the s_x.
1041 // This is a simplified method to determine the error parameterization for s_x and s_y as compared to the one
1042 // implemented in ProcessSigma(). For more details on cluster error parameterization please see also
1043 // AliTRDcluster::SetSigmaY2()
1044 //
1045 // The representation of dy=f(y_cen, x_drift| layer) can be also used to estimate the systematic shift in the r-phi
1046 // coordinate resulting from imperfection in the cluster shape parameterization. From the expresion of the shift derived
1047 // in ProcessMean() with phi=exb one gets:
1048 // BEGIN_LATEX
1049 // <#Delta y>= <#delta x> * (tg(#alpha_{L})-h*dz/dx) + <#delta y - #delta x * tg(#alpha_{L})>
1050 // <#Delta y>(y_{cen})= -h*<#delta x>(x_{drift}, q_{cl}) * dz/dx + #delta y(y_{cen}, ...)
1051 // END_LATEX
1052 // where all dependences are made explicit. This last expression can be used in two ways:
1053 //   - by average on the dz/dx we can determine directly dy (the method implemented here)
1054 //   - by plotting as a function of dzdx one can determine both dx and dy components in an independent method.
1055 //Begin_Html
1056 //<img src="TRD/clusterYcorr.gif">
1057 //End_Html
1058 // Author
1059 // Alexandru Bercuci <A.Bercuci@gsi.de>
1060
1061   TObjArray *arr(NULL);
1062   TH3S *h3r(NULL), *h3t(NULL);
1063   if(!(arr= (TObjArray*)fContainer->At(kYRes))) {
1064     AliError("Missing resolution container");
1065     return kFALSE;
1066   }
1067   if(!(h3r = (TH3S*)arr->At(0))){
1068     AliError("Missing resolution pw/q histo");
1069     return kFALSE;
1070   } else if(!(Int_t)h3r->GetEntries()){
1071     AliError("Empty resolution pw/q histo");
1072     return kFALSE;
1073   }
1074   if(!(h3t = (TH3S*)arr->At(2))){
1075     AliError("Missing resolution t histo");
1076     return kFALSE;
1077   } else if(!(Int_t)h3t->GetEntries()){
1078     AliError("Empty resolution t histo");
1079     return kFALSE;
1080   }
1081
1082   // local variables
1083   Double_t x(0.), y(0.), ex(0.), ey(0.);
1084   Float_t tmin(6.5), tmax(20.5), tmed(0.5*(tmin+tmax));
1085   TGraphErrors *g[2]; TH1 *h(NULL);
1086   g[0] = new TGraphErrors();
1087   g[0]->SetMarkerStyle(24);g[0]->SetMarkerColor(kBlue);g[0]->SetLineColor(kBlue);
1088   g[1] = new TGraphErrors();
1089   g[1]->SetMarkerStyle(24);g[1]->SetMarkerColor(kRed);g[1]->SetLineColor(kRed);
1090
1091   // PROCESS RESOLUTION VS TB
1092   TF1 fsx("fsx", "[0]*[0]+[1]*[1]*[2]*0.1*(x-[3])", tmin, tmax);
1093   fsx.SetParNames("#sqrt{<#sigma^{2}(prf, q)>}(t_{med})", "D_{T}", "v_{drift}", "t_{med}");
1094   fsx.FixParameter(1, fDt);
1095   fsx.SetParameter(2, fVdrift);
1096   fsx.FixParameter(3, tmed);
1097   if(!AliTRDresolution::Process((TH2*)h3r->Project3D("yx"), g)) return kFALSE;
1098   for(Int_t ip(0); ip<g[1]->GetN(); ip++){
1099     g[1]->GetPoint(ip, x, y);ex = g[1]->GetErrorX(ip); ey = g[1]->GetErrorY(ip);
1100     g[1]->SetPoint(ip, x, y*y);g[1]->SetPointError(ip, ex, 2*y*ey);
1101   }
1102   g[1]->Fit(&fsx, "Q", "", tmin, tmax);
1103   if(fCanvas){
1104     g[1]->Draw("apl");
1105     fCanvas->Modified(); fCanvas->Update();
1106     h = g[1]->GetHistogram();
1107     h->SetTitle(fsx.GetTitle());
1108     h->GetXaxis()->SetTitle("t [tb]");h->GetXaxis()->CenterTitle();
1109     h->GetYaxis()->SetTitle("#sigma^{2} (y) [cm^{2}]");h->GetYaxis()->CenterTitle();
1110     if(IsSaveAs()) fCanvas->SaveAs(Form("D%03d_ResNormTrack_tb.gif", fDet));
1111     else gSystem->Sleep(100);
1112   }
1113
1114   // define model for resolution vs pw
1115   TF1 fg("fg", "gaus", -.5, .5); fg.FixParameter(1, 0.);
1116   TF1 fs("fs", "[0]*[0]*exp(-1*(x/[1])**2)+[2]", -.5, .5);
1117   fs.SetParNames("<#sigma^{max}(q,prf)>_{q}", "#sigma(pw)", "D_{T}^{2}*<x>");
1118   h3r->GetXaxis()->SetRange(tmin, tmax);
1119   if(!AliTRDresolution::Process((TH2*)h3r->Project3D("zy"), g, 200)) return kFALSE;
1120   for(Int_t ip(0); ip<g[1]->GetN(); ip++){
1121     g[1]->GetPoint(ip, x, y); ex = g[1]->GetErrorX(ip); ey = g[1]->GetErrorY(ip);
1122     g[1]->SetPoint(ip, x, y*y);g[1]->SetPointError(ip, ex, 2.*y*ey);
1123   }
1124   g[1]->Fit(&fg, "QR");
1125   fs.SetParameter(0, TMath::Sqrt(fg.GetParameter(0)));
1126   fs.SetParameter(1, fg.GetParameter(2));
1127   Float_t sdiff(fDt*fDt*fsx.GetParameter(2)*tmed*0.1);
1128   fs.SetParameter(2, sdiff);
1129   fs.SetParLimits(2, 0.1*sdiff, 1.9*sdiff);
1130   g[1]->Fit(&fs, "QR");
1131   if(fCanvas){
1132     g[1]->Draw("apl");
1133     fCanvas->Modified(); fCanvas->Update();
1134     h = g[1]->GetHistogram();
1135     h->SetTitle(fs.GetTitle());
1136     h->GetXaxis()->SetTitle("pw");h->GetXaxis()->CenterTitle();
1137     h->GetYaxis()->SetTitle("#sigma^{2} (y) [cm^{2}]");h->GetYaxis()->CenterTitle();
1138     if(IsSaveAs()) fCanvas->SaveAs(Form("D%03d_ResNormTrack_pw.gif", fDet));
1139     else gSystem->Sleep(100);
1140   }
1141
1142   AliDebug(2, Form("<s(q,prf)>[mum] = %7.3f", 1.e4*TMath::Sqrt(fsx.Eval(0.))));
1143   AliDebug(2, Form("<s(q)>[mum] = %7.3f", 1.e4*TMath::Sqrt(fs.Eval(-0.5)-fs.GetParameter(2))));
1144   AliDebug(2, Form("<s(x)>[mum] = %7.3f(prf) %7.3f(diff)", 1.e4*TMath::Sqrt(fs.GetParameter(2)), 1.e4*TMath::Sqrt(sdiff)));
1145
1146   // define model for resolution vs q
1147   TF1 fq("fq", "[0]*[0]*exp(-1*[1]*(x-[2])**2)+[2]", 2.5, 5.5);
1148   fq.SetParNames("<#sigma^{max}(q,prf)>_{prf}", "slope","mean", "D_{T}^{2}*<x>");
1149   if(!AliTRDresolution::Process((TH2*)h3t->Project3D("yx"), g)) return kFALSE;
1150   for(Int_t ip(0); ip<g[1]->GetN(); ip++){
1151     g[1]->GetPoint(ip, x, y); ex = g[1]->GetErrorX(ip); ey = g[1]->GetErrorY(ip);
1152     g[1]->SetPoint(ip, x, y*y);g[1]->SetPointError(ip, ex, 2.*y*ey);
1153   }
1154   fq.SetParameter(0, 8.e-2); fq.SetParLimits(0, 0., 1.);
1155   fq.SetParameter(1, 1.); //fq.SetParLimits(1, -1., 0.);
1156   fq.SetParameter(3, sdiff); fq.SetParLimits(3, 0.1*sdiff, 1.9*sdiff);
1157   g[1]->Fit(&fq, "QR");
1158 //   AliDebug(2, Form("<sq>[mum] = %7.3f", 1.e4*TMath::Sqrt(fs.Eval(-0.5)-fs.GetParameter(2)));
1159 //   AliDebug(2, Form("<sx>[mum] = %7.3f(prf) %7.3f(diff)", 1.e4*TMath::Sqrt(fs.Eval(-0.5)-fs.GetParameter(2)), 1.e4*TMath::Sqrt(sdiff)));
1160   if(fCanvas){
1161     g[1]->Draw("apl");
1162     fCanvas->Modified(); fCanvas->Update();
1163     h = g[1]->GetHistogram();
1164     h->SetTitle(fs.GetTitle());
1165     h->GetXaxis()->SetTitle("log(q) [a.u.]");h->GetXaxis()->CenterTitle();
1166     h->GetYaxis()->SetTitle("#sigma^{2} (y) [cm^{2}]");h->GetYaxis()->CenterTitle();
1167     if(IsSaveAs()) fCanvas->SaveAs(Form("D%03d_ResNormTrack_q.gif", fDet));
1168     else gSystem->Sleep(100);
1169   }
1170   return kTRUE;
1171 }
1172
1173 //_______________________________________________________
1174 void AliTRDclusterResolution::ProcessSigma()
1175 {
1176 // As the r-phi coordinate is the only one which is measured by the TRD detector we have to rely on it to
1177 // estimate both the radial (x) and r-phi (y) errors. This method is based on the following assumptions.
1178 // The measured error in the y direction is the sum of the intrinsic contribution of the r-phi measurement
1179 // with the contribution of the radial measurement - because x is not a parameter of Alice track model (Kalman).
1180 // BEGIN_LATEX
1181 // #sigma^{2}|_{y} = #sigma^{2}_{y*} + #sigma^{2}_{x*}
1182 // END_LATEX
1183 // In the general case
1184 // BEGIN_LATEX
1185 // #sigma^{2}_{y*} = #sigma^{2}_{y} + tg^{2}(#alpha_{L})#sigma^{2}_{x_{drift}}
1186 // #sigma^{2}_{x*} = tg^{2}(#phi - #alpha_{L})*(#sigma^{2}_{x_{drift}} + #sigma^{2}_{x_{0}} + tg^{2}(#alpha_{L})*x^{2}/12)
1187 // END_LATEX
1188 // where we have explicitely show the lorentz angle correction on y and the projection of radial component on the y
1189 // direction through the track angle in the bending plane (phi). Also we have shown that the radial component in the
1190 // last equation has twp terms, the drift and the misalignment (x_0). For ideal geometry or known misalignment one
1191 // can solve the equation
1192 // BEGIN_LATEX
1193 // #sigma^{2}|_{y} = tg^{2}(#phi - #alpha_{L})*(#sigma^{2}_{x} + tg^{2}(#alpha_{L})*x^{2}/12)+ [#sigma^{2}_{y} + tg^{2}(#alpha_{L})#sigma^{2}_{x}]
1194 // END_LATEX
1195 // by fitting a straight line:
1196 // BEGIN_LATEX
1197 // #sigma^{2}|_{y} = a(x_{cl}, z_{cl}) * tg^{2}(#phi - #alpha_{L}) + b(x_{cl}, z_{cl})
1198 // END_LATEX
1199 // the error parameterization will be given by:
1200 // BEGIN_LATEX
1201 // #sigma_{x} (x_{cl}, z_{cl}) = #sqrt{a(x_{cl}, z_{cl}) - tg^{2}(#alpha_{L})*x^{2}/12}
1202 // #sigma_{y} (x_{cl}, z_{cl}) = #sqrt{b(x_{cl}, z_{cl}) - #sigma^{2}_{x} (x_{cl}, z_{cl}) * tg^{2}(#alpha_{L})}
1203 // END_LATEX
1204 // Below there is an example of such dependency.
1205 //Begin_Html
1206 //<img src="TRD/clusterSigmaMethod.gif">
1207 //End_Html
1208 //
1209 // The error parameterization obtained by this method are implemented in the functions AliTRDcluster::GetSX() and
1210 // AliTRDcluster::GetSYdrift(). For an independent method to determine s_y as a function of drift length check the
1211 // function ProcessCenterPad(). One has to keep in mind that while this method return the mean s_y over the distance
1212 // to pad center distribution the other method returns the *STANDARD* value at center=0 (maximum). To recover the
1213 // standard value one has to solve the obvious equation:
1214 // BEGIN_LATEX
1215 // #sigma_{y}^{STANDARD} = #frac{<#sigma_{y}>}{#int{s exp(s^{2}/#sigma) ds}}
1216 // END_LATEX
1217 // with "<s_y>" being the value calculated here and "sigma" the width of the s_y distribution calculated in
1218 // ProcessCenterPad().
1219 //
1220 // Author
1221 // Alexandru Bercuci <A.Bercuci@gsi.de>
1222
1223   TObjArray *arr = (TObjArray*)fContainer->At(kSigm);
1224   if(!arr){
1225     AliWarning("Missing dy=f(x_d, d_w) container");
1226     return;
1227   }
1228
1229   // init visualization
1230   TGraphErrors *ggs = NULL;
1231   TGraph *line = NULL;
1232   if(fCanvas){
1233     ggs = new TGraphErrors();
1234     line = new TGraph();
1235     line->SetLineColor(kRed);line->SetLineWidth(2);
1236   }
1237
1238   // init logistic support
1239   TF1 f("f", "gaus", -.5, .5);
1240   TLinearFitter gs(1,"pol1");
1241   TH1 *hFrame=NULL;
1242   TH1D *h1 = NULL; TH3S *h3=NULL;
1243   TAxis *ax = NULL;
1244   Double_t exb2 = fExB*fExB;
1245   AliTRDcluster c;
1246   TTree *t = (TTree*)fResults->At(kSigm);
1247   for(Int_t ix=0; ix<AliTRDseedV1::kNtb; ix++){
1248     if(!(h3=(TH3S*)arr->At(ix))) continue;
1249     c.SetPadTime(ix);
1250     fX = c.GetXloc(fT0, fVdrift);
1251     fT = c.GetLocalTimeBin(); // ideal
1252     printf(" pad time[%d] local[%f]\n", ix, fT);
1253     for(Int_t iz=1; iz<=h3->GetXaxis()->GetNbins(); iz++){
1254       ax = h3->GetXaxis();
1255       ax->SetRange(iz, iz);
1256       fZ = ax->GetBinCenter(iz);
1257
1258       // reset visualization
1259       if(fCanvas){
1260         new(ggs) TGraphErrors();
1261         ggs->SetMarkerStyle(7);
1262       }
1263       gs.ClearPoints();
1264
1265       for(Int_t ip=1; ip<=h3->GetYaxis()->GetNbins(); ip++){
1266         ax = h3->GetYaxis();
1267         ax->SetRange(ip, ip);
1268         Double_t tgl = ax->GetBinCenter(ip);
1269         // finish navigation in the HnSparse
1270
1271         //if(TMath::Abs(dydx)>0.18) continue;
1272         Double_t tgg = (tgl-fExB)/(1.+tgl*fExB);
1273         Double_t tgg2 = tgg*tgg;
1274
1275         h1 = (TH1D*)h3->Project3D("z");
1276         Int_t entries = (Int_t)h1->Integral();
1277         if(entries < 50) continue;
1278         //Adjust(&f, h1);
1279         h1->Fit(&f, "QN");
1280
1281         Double_t s2  = f.GetParameter(2)*f.GetParameter(2);
1282         Double_t s2e = 2.*f.GetParameter(2)*f.GetParError(2);
1283         // Fill sy^2 = f(tg^2(phi-a_L))
1284         gs.AddPoint(&tgg2, s2, s2e);
1285
1286         if(!ggs) continue;
1287         Int_t jp = ggs->GetN();
1288         ggs->SetPoint(jp, tgg2, s2);
1289         ggs->SetPointError(jp, 0., s2e);
1290       }
1291       // TODO here a more robust fit method has to be provided
1292       // for which lower boundaries on the parameters have to
1293       // be imposed. Unfortunately the Minuit fit does not work
1294       // for the TGraph in the case of B not 0.
1295       if(gs.Eval()) continue;
1296
1297       fR[0] = gs.GetParameter(1) - fX*fX*exb2/12.;
1298       AliDebug(3, Form("  s2x+x2=%f ang=%f s2x=%f", gs.GetParameter(1), fX*fX*exb2/12., fR[0]));
1299       fR[0] = TMath::Max(fR[0], Float_t(4.e-4));
1300
1301       // s^2_y  = s0^2_y + tg^2(a_L) * s^2_x
1302       // s0^2_y = f(D_L)*x + s_PRF^2
1303       fR[2]= gs.GetParameter(0)-exb2*fR[0];
1304       AliDebug(3, Form("  s2y+s2x=%f s2y=%f", fR[0], fR[2]));
1305       fR[2] = TMath::Max(fR[2], Float_t(2.5e-5));
1306       fR[0] = TMath::Sqrt(fR[0]);
1307       fR[1] = .5*gs.GetParError(1)/fR[0];
1308       fR[2] = TMath::Sqrt(fR[2]);
1309       fR[3] = gs.GetParError(0)+exb2*exb2*gs.GetParError(1);
1310       t->Fill();
1311       AliDebug(2, Form("xd=%4.2f[cm] sx=%6.1f[um] sy=%5.1f[um]", fX, 1.e4*fR[0], 1.e4*fR[2]));
1312
1313       if(!fCanvas) continue;
1314       fCanvas->cd(); fCanvas->SetLogx(); //fCanvas->SetLogy();
1315       if(!hFrame){
1316         fCanvas->SetMargin(0.15, 0.01, 0.1, 0.01);
1317         hFrame=new TH1I("hFrame", "", 100, 0., .3);
1318         hFrame->SetMinimum(0.);hFrame->SetMaximum(.005);
1319         hFrame->SetXTitle("tg^{2}(#phi-#alpha_{L})");
1320         hFrame->SetYTitle("#sigma^{2}y[cm^{2}]");
1321         hFrame->GetYaxis()->SetTitleOffset(2.);
1322         hFrame->SetLineColor(1);hFrame->SetLineWidth(1);
1323         hFrame->Draw();
1324       } else hFrame->Reset();
1325       Double_t xx = 0., dxx=.2/50;
1326       for(Int_t ip=0;ip<50;ip++){
1327         line->SetPoint(ip, xx,  gs.GetParameter(0)+xx*gs.GetParameter(1));
1328         xx+=dxx;
1329       }
1330       ggs->Draw("pl"); line->Draw("l");
1331       fCanvas->Modified(); fCanvas->Update();
1332       if(IsSaveAs()) fCanvas->SaveAs(Form("Figures/ProcessSigma_z[%5.3f]_x[%5.3f].gif", fZ, fX));
1333       else gSystem->Sleep(100);
1334     }
1335   }
1336   return;
1337 }
1338
1339 //_______________________________________________________
1340 void AliTRDclusterResolution::ProcessMean()
1341 {
1342 // By this method the cluster shift in r-phi and radial directions can be estimated by comparing with the MC.
1343 // The resolution of the cluster corrected for pad tilt with respect to MC in the r-phi (measuring) plane can be
1344 // expressed by:
1345 // BEGIN_LATEX
1346 // #Delta y=w - y_{MC}(x_{cl})
1347 // w = y_{cl}^{'} + h*(z_{MC}(x_{cl})-z_{cl})
1348 // y_{MC}(x_{cl}) = y_{0} - dy/dx*x_{cl}
1349 // z_{MC}(x_{cl}) = z_{0} - dz/dx*x_{cl}
1350 // y_{cl}^{'} = y_{cl}-x_{cl}*tg(#alpha_{L})
1351 // END_LATEX
1352 // where x_cl is the drift length attached to a cluster, y_cl is the r-phi coordinate of the cluster measured by
1353 // charge sharing on adjacent pads and y_0 and z_0 are MC reference points (as example the track references at
1354 // entrance/exit of a chamber). If we suppose that both r-phi (y) and radial (x) coordinate of the clusters are
1355 // affected by errors we can write
1356 // BEGIN_LATEX
1357 // x_{cl} = x_{cl}^{*} + #delta x
1358 // y_{cl} = y_{cl}^{*} + #delta y
1359 // END_LATEX
1360 // where the starred components are the corrected values. Thus by definition the following quantity
1361 // BEGIN_LATEX
1362 // #Delta y^{*}= w^{*} - y_{MC}(x_{cl}^{*})
1363 // END_LATEX
1364 // has 0 average over all dependency. Using this decomposition we can write:
1365 // BEGIN_LATEX
1366 // <#Delta y>=<#Delta y^{*}> + <#delta x * (dy/dx-h*dz/dx) + #delta y - #delta x * tg(#alpha_{L})>
1367 // END_LATEX
1368 // which can be transformed to the following linear dependence:
1369 // BEGIN_LATEX
1370 // <#Delta y>= <#delta x> * (dy/dx-h*dz/dx) + <#delta y - #delta x * tg(#alpha_{L})>
1371 // END_LATEX
1372 // if expressed as function of dy/dx-h*dz/dx. Furtheremore this expression can be plotted for various clusters
1373 // i.e. we can explicitely introduce the diffusion (x_cl) and drift cell - anisochronity (z_cl) dependences. From
1374 // plotting this dependence and linear fitting it with:
1375 // BEGIN_LATEX
1376 // <#Delta y>= a(x_{cl}, z_{cl}) * (dy/dx-h*dz/dx) + b(x_{cl}, z_{cl})
1377 // END_LATEX
1378 // the systematic shifts will be given by:
1379 // BEGIN_LATEX
1380 // #delta x (x_{cl}, z_{cl}) = a(x_{cl}, z_{cl})
1381 // #delta y (x_{cl}, z_{cl}) = b(x_{cl}, z_{cl}) + a(x_{cl}, z_{cl}) * tg(#alpha_{L})
1382 // END_LATEX
1383 // Below there is an example of such dependency.
1384 //Begin_Html
1385 //<img src="TRD/clusterShiftMethod.gif">
1386 //End_Html
1387 //
1388 // The occurance of the radial shift is due to the following conditions
1389 //   - the approximation of a constant drift velocity over the drift length (larger drift velocities close to
1390 //     cathode wire plane)
1391 //   - the superposition of charge tails in the amplification region (first clusters appear to be located at the
1392 //     anode wire)
1393 //   - the superposition of charge tails in the drift region (shift towards anode wire)
1394 //   - diffusion effects which convolute with the TRF thus enlarging it
1395 //   - approximate knowledge of the TRF (approximate measuring in test beam conditions)
1396 //
1397 // The occurance of the r-phi shift is due to the following conditions
1398 //   - approximate model for cluster shape (LUT)
1399 //   - rounding-up problems
1400 //
1401 // The numerical results for ideal simulations for the radial and r-phi shifts are displayed below and used
1402 // for the cluster reconstruction (see the functions AliTRDcluster::GetXcorr() and AliTRDcluster::GetYcorr()).
1403 //Begin_Html
1404 //<img src="TRD/clusterShiftX.gif">
1405 //<img src="TRD/clusterShiftY.gif">
1406 //End_Html
1407 // More details can be found in the presentation given during the TRD
1408 // software meeting at the end of 2008 and beginning of year 2009, published on indico.cern.ch.
1409 //
1410 // Author
1411 // Alexandru Bercuci <A.Bercuci@gsi.de>
1412
1413
1414
1415   TObjArray *arr = (TObjArray*)fContainer->At(kMean);
1416   if(!arr){
1417     AliWarning("Missing dy=f(x_d, d_w) container");
1418     return;
1419   }
1420
1421   // init logistic support
1422   TF1 f("f", "gaus", -.5, .5);
1423   TF1 line("l", "[0]+[1]*x", -.15, .15);
1424   TGraphErrors *gm = new TGraphErrors();
1425   TH1 *hFrame=NULL;
1426   TH1D *h1 = NULL; TH3S *h3 =NULL;
1427   TAxis *ax = NULL;
1428
1429   AliDebug(1, Form("Calibrate for Det[%3d] t0[%5.3f] vd[%5.3f]", fDet, fT0, fVdrift));
1430
1431   AliTRDcluster c;
1432   TTree *t = (TTree*)fResults->At(kMean);
1433   for(Int_t ix=0; ix<AliTRDseedV1::kNtb; ix++){
1434     if(!(h3=(TH3S*)arr->At(ix))) continue;
1435     c.SetPadTime(ix);
1436     fX = c.GetXloc(fT0, fVdrift);
1437     fT = c.GetLocalTimeBin();
1438     for(Int_t iz=1; iz<=h3->GetXaxis()->GetNbins(); iz++){
1439       ax = h3->GetXaxis();
1440       ax->SetRange(iz, iz);
1441       fZ = ax->GetBinCenter(iz);
1442
1443       // reset fitter
1444       new(gm) TGraphErrors();
1445       gm->SetMarkerStyle(7);
1446
1447       for(Int_t ip=1; ip<=h3->GetYaxis()->GetNbins(); ip++){
1448         ax = h3->GetYaxis();
1449         ax->SetRange(ip, ip);
1450         Double_t tgl = ax->GetBinCenter(ip);
1451         // finish navigation in the HnSparse
1452
1453         h1 = (TH1D*)h3->Project3D("z");
1454         Int_t entries = (Int_t)h1->Integral();
1455         if(entries < 50) continue;
1456         //Adjust(&f, h1);
1457         h1->Fit(&f, "QN");
1458
1459         // Fill <Dy> = f(dydx - h*dzdx)
1460         Int_t jp = gm->GetN();
1461         gm->SetPoint(jp, tgl, f.GetParameter(1));
1462         gm->SetPointError(jp, 0., f.GetParError(1));
1463       }
1464       if(gm->GetN()<10) continue;
1465
1466       gm->Fit(&line, "QN");
1467       fR[0] = line.GetParameter(1); // dx
1468       fR[1] = line.GetParError(1);
1469       fR[2] = line.GetParameter(0) + fExB*fR[0]; // xs = dy - tg(a_L)*dx
1470       t->Fill();
1471       AliDebug(2, Form("tb[%02d] xd=%4.2f[cm] dx=%6.2f[um] dy=%6.2f[um]", ix, fX, 1.e4*fR[0], 1.e4*fR[2]));
1472       if(!fCanvas) continue;
1473
1474       fCanvas->cd();
1475       if(!hFrame){
1476         fCanvas->SetMargin(0.1, 0.02, 0.1, 0.01);
1477         hFrame=new TH1I("hFrame", "", 100, -.3, .3);
1478         hFrame->SetMinimum(-.1);hFrame->SetMaximum(.1);
1479         hFrame->SetXTitle("tg#phi-htg#theta");
1480         hFrame->SetYTitle("#Delta y[cm]");
1481         hFrame->GetYaxis()->SetTitleOffset(1.5);
1482         hFrame->SetLineColor(1);hFrame->SetLineWidth(1);
1483         hFrame->Draw();
1484       } else hFrame->Reset();
1485       gm->Draw("pl"); line.Draw("same");
1486       fCanvas->Modified(); fCanvas->Update();
1487       if(IsSaveAs()) fCanvas->SaveAs(Form("Figures/ProcessMean_Z[%5.3f]_TB[%02d].gif", fZ, ix));
1488       else gSystem->Sleep(100);
1489     }
1490   }
1491 }