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   if(!(fInfo = dynamic_cast<TObjArray *>(GetInputData(1)))){
 583     AliError("Cluster array missing.");
 584     return;
 585   }
 586   AliDebug(2, Form("Clusters[%d]", fInfo->GetEntriesFast()));
 587
 588   Int_t det, t, np;
 589   Float_t x, y, z, q, dy, dydx, dzdx, cov[3], covcl[3];
 590   TH3S *h3(NULL); TH2I *h2(NULL);
 591
 592   // define limits around ExB for which x contribution is negligible
 593   const Float_t kAroundZero = 3.5e-2; //(+- 2 deg)
 594
 595   TObjArray *arr0 = (TObjArray*)fContainer->At(kYSys);
 596   TObjArray *arr1 = (TObjArray*)fContainer->At(kYRes);
 597   TObjArray *arr2 = (TObjArray*)fContainer->At(kSigm);
 598
 599   const AliTRDclusterInfo *cli = NULL;
 600   TIterator *iter=fInfo->MakeIterator();
 601   while((cli=dynamic_cast<AliTRDclusterInfo*>((*iter)()))){
 602     if((np = cli->GetNpads())>4) continue;
 603     cli->GetCluster(det, x, y, z, q, t, covcl);
 604
 605     // select cluster according to detector region if specified
 606     if(fDet>=0 && fDet!=det) continue;
 607     if(fCol>=0 && fRow>=0){
 608       Int_t c,r;
 609       cli->GetCenterPad(c, r);
 610       if(TMath::Abs(fCol-c) > 5) continue;
 611       if(TMath::Abs(fRow-r) > 2) continue;
 612     }
 613     dy = cli->GetResolution();
 614     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])));
 615
 616     cli->GetGlobalPosition(y, z, dydx, dzdx, &cov[0]);
 617     Float_t pw(cli->GetYDisplacement());
 618
 619     // systematics as a function of pw and log(q)
 620     // only for dydx = exB + h*dzdx
 621     if(TMath::Abs(dydx-fExB-fH*dzdx) < kAroundZero){
 622       h3 = (TH3S*)arr0->At(0);
 623       h3->Fill(TMath::Log(q), pw, dy);
 624     }
 625     // resolution/pull as a function of pw and log(q)
 626     // only for dydx = 0, ExB=0
 627     if(TMath::Abs(fExB) < kAroundZero &&
 628        TMath::Abs(dydx) < kAroundZero &&
 629        t>5 && t<24 ){
 630       switch(np){
 631       case 3: // MPV np
 632         h3 = (TH3S*)arr1->At(0);
 633         h3->Fill(TMath::Log(q), pw, dy);
 634         h3 = (TH3S*)arr1->At(5);
 635         h3->Fill(TMath::Log(q), pw, t);
 636         break;
 637       case 2: // Min np
 638         h3 = (TH3S*)arr1->At(3);
 639         h3->Fill(TMath::Log(q), pw, dy);
 640         break;
 641       case 4: // Max np
 642         h3 = (TH3S*)arr1->At(4);
 643         h3->Fill(TMath::Log(q), pw, dy);
 644         break;
 645       }
 646       h3 = (TH3S*)arr1->At(1);
 647       h3->Fill(TMath::Log(q), pw, dy/TMath::Sqrt(covcl[0]));
 648     }
 649
 650     // do not use problematic clusters in resolution analysis
 651     // TODO define limits as calibration aware (gain) !!
 652     //if(!AcceptableGain(fGain)) continue;
 653     if(q<20. || q>250.) continue;
 654
 655     // systematic as a function of time bin
 656     // only for dydx = exB + h*dzdx and MPV q
 657     if(TMath::Abs(dydx-fExB-fH*dzdx) < kAroundZero){
 658       h2 = (TH2I*)arr0->At(1);
 659       h2->Fill(t, dy);
 660     }
 661     // systematic as function of tb and tgp
 662     // only for MPV q
 663     h3 = (TH3S*)arr0->At(2);
 664     h3->Fill(t, dydx, dy);
 665
 666     // resolution/pull as a function of time bin
 667     // only for dydx = 0, ExB=0 and MPV q
 668     if(TMath::Abs(fExB) < kAroundZero &&
 669        TMath::Abs(dydx) < kAroundZero){
 670       h3 = (TH3S*)arr1->At(2);
 671       h3->Fill(t, dy, dy/TMath::Sqrt(covcl[0]));
 672     }
 673
 674     // resolution as function of tb, tgp and h*tgt
 675     // only for MPV q
 676     ((TH3S*)arr2->At(t))->Fill(fH*dzdx, dydx, dy);
 677   }
 678 }
 679
 680
 681 //_______________________________________________________
 682 Bool_t AliTRDclusterResolution::PostProcess()
 683 {
 684 // Steer processing of various cluster resolution dependences :
 685 //
 686 //   - process resolution dependency cluster charge
 687 //   if(HasProcess(kYRes)) ProcessCharge();
 688 //   - process resolution dependency on y displacement
 689 //   if(HasProcess(kYSys)) ProcessCenterPad();
 690 //   - process resolution dependency on drift legth and drift cell width
 691 //   if(HasProcess(kSigm)) ProcessSigma();
 692 //   - process systematic shift on drift legth and drift cell width
 693 //   if(HasProcess(kMean)) ProcessMean();
 694
 695   if(!fContainer) return kFALSE;
 696   if(!IsCalibrated()){
 697     AliError("Not calibrated instance.");
 698     return kFALSE;
 699   }
 700   TObjArray *arr = NULL;
 701   TTree *t=NULL;
 702   if(!fResults){
 703     TGraphErrors *g = NULL;
 704     fResults = new TObjArray(kNtasks);
 705     fResults->SetOwner();
 706     fResults->AddAt(arr = new TObjArray(3), kYRes);
 707     arr->SetOwner();
 708     arr->AddAt(g = new TGraphErrors(), 0);
 709     g->SetLineColor(kBlue); g->SetMarkerColor(kBlue);
 710     g->SetMarkerStyle(7);
 711     arr->AddAt(g = new TGraphErrors(), 1);
 712     g->SetLineColor(kRed); g->SetMarkerColor(kRed);
 713     g->SetMarkerStyle(23);
 714     arr->AddAt(g = new TGraphErrors(), 2);
 715     g->SetLineColor(kGreen); g->SetMarkerColor(kGreen);
 716     g->SetMarkerStyle(7);
 717
 718     // pad center dependence
 719     fResults->AddAt(arr = new TObjArray(AliTRDgeometry::kNlayer+1), kYSys);
 720     arr->SetOwner();
 721     arr->AddAt(
 722     t = new TTree("cent", "dy=f(y,x,ly)"), 0);
 723     t->Branch("ly", &fLy, "ly/B");
 724     t->Branch("t", &fT, "t/F");
 725     t->Branch("y", &fY, "y/F");
 726     t->Branch("m", &fR[0], "m[2]/F");
 727     t->Branch("s", &fR[2], "s[2]/F");
 728     t->Branch("pm", &fP[0], "pm[2]/F");
 729     t->Branch("ps", &fP[2], "ps[2]/F");
 730     for(Int_t il=1; il<=AliTRDgeometry::kNlayer; il++){
 731       arr->AddAt(g = new TGraphErrors(), il);
 732       g->SetLineColor(il); g->SetLineStyle(il);
 733       g->SetMarkerColor(il);g->SetMarkerStyle(4);
 734     }
 735
 736
 737     fResults->AddAt(t = new TTree("sigm", "dy=f(dw,x,dydx)"), kSigm);
 738     t->Branch("t", &fT, "t/F");
 739     t->Branch("x", &fX, "x/F");
 740     t->Branch("z", &fZ, "z/F");
 741     t->Branch("sx", &fR[0], "sx[2]/F");
 742     t->Branch("sy", &fR[2], "sy[2]/F");
 743
 744
 745     fResults->AddAt(t = new TTree("mean", "dy=f(dw,x,dydx - h dzdx)"), kMean);
 746     t->Branch("t", &fT, "t/F");
 747     t->Branch("x", &fX, "x/F");
 748     t->Branch("z", &fZ, "z/F");
 749     t->Branch("dx", &fR[0], "dx[2]/F");
 750     t->Branch("dy", &fR[2], "dy[2]/F");
 751   } else {
 752     TObject *o = NULL;
 753     TIterator *iter=fResults->MakeIterator();
 754     while((o=((*iter)()))) o->Clear(); // maybe it is wrong but we should never reach this point
 755   }
 756
 757   // process resolution dependency on charge
 758   if(HasProcess(kYRes)) ProcessCharge();
 759
 760   // process resolution dependency on y displacement
 761   if(HasProcess(kYSys)) ProcessNormalTracks();
 762
 763   // process resolution dependency on drift legth and drift cell width
 764   if(HasProcess(kSigm)) ProcessSigma();
 765
 766   // process systematic shift on drift legth and drift cell width
 767   if(HasProcess(kMean)) ProcessMean();
 768
 769   return kTRUE;
 770 }
 771
 772 //_______________________________________________________
 773 Bool_t AliTRDclusterResolution::LoadCalibration()
 774 {
 775 // Retrieve calibration parameters from OCDB, drift velocity and t0 for the detector region specified by
 776 // a previous call to AliTRDclusterResolution::SetCalibrationRegion().
 777
 778   AliCDBManager *cdb = AliCDBManager::Instance(); // check access OCDB
 779   if(cdb->GetRun() < 0){
 780     AliError("OCDB manager not properly initialized");
 781     return kFALSE;
 782   }
 783   // check magnetic field
 784   if(!TGeoGlobalMagField::Instance() || !TGeoGlobalMagField::Instance()->IsLocked()){
 785     AliError("Magnetic field not available.");
 786     return kFALSE;
 787   }
 788
 789   AliTRDcalibDB *fCalibration  = AliTRDcalibDB::Instance();
 790   AliTRDCalROC  *fCalVdriftROC(fCalibration->GetVdriftROC(fDet>=0?fDet:0))
 791                ,*fCalT0ROC(fCalibration->GetT0ROC(fDet>=0?fDet:0));
 792   const AliTRDCalDet  *fCalVdriftDet = fCalibration->GetVdriftDet();
 793   const AliTRDCalDet  *fCalT0Det = fCalibration->GetT0Det();
 794
 795   if(IsUsingCalibParam(kVdrift)){
 796     fVdrift = fCalVdriftDet->GetValue(fDet>=0?fDet:0);
 797     if(fCol>=0 && fRow>=0) fVdrift*= fCalVdriftROC->GetValue(fCol, fRow);
 798   }
 799   fExB    = AliTRDCommonParam::Instance()->GetOmegaTau(fVdrift);
 800   AliTRDCommonParam::Instance()->GetDiffCoeff(fDt, fDl, fVdrift);
 801   if(IsUsingCalibParam(kT0)){
 802     fT0     = fCalT0Det->GetValue(fDet>=0?fDet:0);
 803     if(fCol>=0 && fRow>=0) fT0 *= fCalT0ROC->GetValue(fCol, fRow);
 804   }
 805   if(IsUsingCalibParam(kGain)) fGain = (fCol>=0 && fRow>=0)?fCalibration-> GetGainFactor(fDet, fCol, fRow):fCalibration-> GetGainFactorAverage(fDet);
 806
 807   SetBit(kCalibrated);
 808
 809   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));
 810
 811   return kTRUE;
 812 }
 813
 814 //_______________________________________________________
 815 Bool_t AliTRDclusterResolution::LoadGlobalChamberPosition()
 816 {
 817 // Retrieve global chamber position from alignment
 818 // a previous call to AliTRDclusterResolution::SetCalibrationRegion() is mandatory.
 819
 820   TGeoHMatrix *matrix(NULL);
 821   Double_t loc[] = {0., 0., 0.}, glb[] = {0., 0., 0.};
 822   AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
 823   if(!(matrix= geo->GetClusterMatrix(fDet))) {
 824     AliFatal(Form("Could not get transformation matrix for detector %d.", fDet));
 825     return kFALSE;
 826   }
 827   matrix->LocalToMaster(loc, glb);
 828   fXch = glb[0]+ AliTRDgeometry::AnodePos()-.5*AliTRDgeometry::AmThick() - AliTRDgeometry::DrThick();
 829   AliTRDpadPlane *pp(geo->GetPadPlane(fDet));
 830   fH = TMath::Tan(pp->GetTiltingAngle()*TMath::DegToRad());
 831   fZch = 0.;
 832   if(fRow>=0){
 833     fZch = pp->GetRowPos(fRow)+0.5*pp->GetLengthIPad();
 834   }else{
 835     Int_t  nrows(pp->GetNrows());
 836     Float_t zmax(pp->GetRow0()),
 837             zmin(zmax - 2 * pp->GetLengthOPad()
 838                 - (nrows-2) * pp->GetLengthIPad()
 839                 - (nrows-1) * pp->GetRowSpacing());
 840     fZch = 0.5*(zmin+zmax);
 841   }
 842
 843   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));
 844   SetBit(kGlobal);
 845   return kTRUE;
 846 }
 847
 848 //_______________________________________________________
 849 void AliTRDclusterResolution::SetCalibrationRegion(Int_t det, Int_t col, Int_t row)
 850 {
 851 // Set calibration region in terms of detector and pad.
 852 // By default detector 0 mean values are considered.
 853
 854   if(det>=0 && det<AliTRDgeometry::kNdet){
 855     fDet = det;
 856     if(col>=0 && row>=0){
 857       // check pad col/row for detector
 858       AliTRDgeometry *geo = AliTRDinfoGen::Geometry();
 859       AliTRDpadPlane *pp(geo->GetPadPlane(fDet));
 860       if(fCol>=pp->GetNcols() ||
 861         fRow>=pp->GetNrows()){
 862         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()));
 863         fCol = -1; fRow=-1;
 864       } else {
 865         fCol = col;
 866         fRow = row;
 867       }
 868     }
 869   } else {
 870     AliFatal(Form("Detector index outside range [0 %d].", AliTRDgeometry::kNdet));
 871   }
 872   return;
 873 }
 874
 875 //_______________________________________________________
 876 void AliTRDclusterResolution::SetVisual()
 877 {
 878   if(fCanvas) return;
 879   fCanvas = new TCanvas("clResCanvas", "Cluster Resolution Visualization", 10, 10, 600, 600);
 880 }
 881
 882 //_______________________________________________________
 883 void AliTRDclusterResolution::ProcessCharge()
 884 {
 885 // Resolution as a function of cluster charge.
 886 //
 887 // As described in the function ProcessCenterPad() the error parameterization for clusters for phi = a_L can be
 888 // written as:
 889 // BEGIN_LATEX
 890 // #sigma_{y}^{2} = #sigma_{y}^{2}|_{B=0} + tg^{2}(#alpha_{L})*#sigma_{x}^{2}
 891 // END_LATEX
 892 // with the contribution in case of B=0 given by:
 893 // BEGIN_LATEX
 894 // #sigma_{y}|_{B=0} = #sigma_{diff}*Gauss(0, s_{ly}) + #delta_{#sigma}(q)
 895 // END_LATEX
 896 // which further can be simplified to:
 897 // BEGIN_LATEX
 898 // <#sigma_{y}|_{B=0}>(q) = <#sigma_{y}> + #delta_{#sigma}(q)
 899 // <#sigma_{y}> = #int{f(q)#sigma_{y}dq}
 900 // END_LATEX
 901 // The results for s_y and f(q) are displayed below:
 902 //Begin_Html
 903 //<img src="TRD/clusterQerror.gif">
 904 //End_Html
 905 // The function has to extended to accomodate gain calibration scalling and errors.
 906 //
 907 // Author
 908 // Alexandru Bercuci <A.Bercuci@gsi.de>
 909
 910
 911
 912   TObjArray *arr(NULL);
 913   if(!(arr = (TObjArray*)fContainer->At(kYSys))) {
 914     AliError("Missing systematic container");
 915     return;
 916   }
 917   TH3S *h3s(NULL);
 918   if(!(h3s = (TH3S*)arr->At(0))){
 919     AliError("Missing systematic histo");
 920     return;
 921   }
 922   // PROCESS SYSTEMATIC
 923   Float_t tmin(6.5), tmax(20.5), tmed(0.5*(tmin+tmax));
 924   TGraphErrors *g[2]; TH1 *h(NULL);
 925   g[0] = new TGraphErrors();
 926   g[0]->SetMarkerStyle(24);g[0]->SetMarkerColor(kBlue);g[0]->SetLineColor(kBlue);
 927   g[1] = new TGraphErrors();
 928   g[1]->SetMarkerStyle(24);g[1]->SetMarkerColor(kRed);g[1]->SetLineColor(kRed);
 929   // define model for systematic shift vs pw
 930   TF1 fm("fm", "[0]+[1]*sin(x*[2])", -.45,.45);
 931   fm.SetParameter(0, 0.); fm.SetParameter(1, 1.e-2); fm.FixParameter(2, TMath::TwoPi());
 932   fm.SetParNames("#deltay", "#pm#delta", "2*#pi");
 933   h3s->GetXaxis()->SetRange(tmin, tmax);
 934   if(!AliTRDresolution::Process((TH2*)h3s->Project3D("zy"), g))return;
 935   g[0]->Fit(&fm, "QR");
 936   if(fCanvas){
 937     g[0]->Draw("apl");
 938     fCanvas->Modified(); fCanvas->Update();
 939     h = g[0]->GetHistogram();
 940     h->SetTitle(fm.GetTitle());
 941     h->GetXaxis()->SetTitle("pw");h->GetXaxis()->CenterTitle();
 942     h->GetYaxis()->SetTitle("#Delta y[cm]");h->GetYaxis()->CenterTitle();
 943     if(IsSaveAs()) fCanvas->SaveAs(Form("D%03d_SysNormTrack_pw.gif", fDet));
 944     else gSystem->Sleep(100);
 945   }
 946
 947   // define model for systematic shift vs tb
 948   TF1 fx("fx", "[0]+0.1*[1]*(x-[2])", tmin, tmax);
 949   fx.SetParNames("#deltay", "#deltay/t", "<t>");
 950   fx.FixParameter(2, tmed);
 951   h3s->GetXaxis()->UnZoom();
 952   if(!AliTRDresolution::Process((TH2*)h3s->Project3D("zx"), g)) return;
 953   g[0]->Fit(&fx, "Q", "", tmin, tmax);
 954   if(fCanvas){
 955     g[0]->Draw("apl");
 956     fCanvas->Modified(); fCanvas->Update();
 957     h = g[0]->GetHistogram();
 958     h->SetTitle(fx.GetTitle());
 959     h->GetXaxis()->SetTitle("t [tb]");h->GetXaxis()->CenterTitle();
 960     h->GetYaxis()->SetTitle("#Delta y[cm]");h->GetYaxis()->CenterTitle();
 961     if(IsSaveAs()) fCanvas->SaveAs(Form("D%03d_SysNormTrack_tb.gif", fDet));
 962     else gSystem->Sleep(100);
 963   }
 964
 965   TH3S *h3(NULL);
 966   if(!(h3 = (TH3S*)fContainer->At(kYRes))) {
 967     AliWarning("Missing dy=f(Q) histo");
 968     return;
 969   }
 970   TF1 f("f", "gaus", -.5, .5);
 971   TAxis *ax(NULL);
 972   TH1 *h1(NULL);
 973
 974   // compute mean error on x
 975   Double_t s2x = 0.;
 976   for(Int_t ix=5; ix<AliTRDseedV1::kNtb; ix++){
 977     // retrieve error on the drift length
 978     s2x += AliTRDcluster::GetSX(ix);
 979   }
 980   s2x /= (AliTRDseedV1::kNtb-5); s2x *= s2x;
 981   //Double_t exb2 = fExB*fExB;
 982
 983   arr = (TObjArray*)fResults->At(kYRes);
 984   TGraphErrors *gqm = (TGraphErrors*)arr->At(0);
 985   TGraphErrors *gqs = (TGraphErrors*)arr->At(1);
 986   TGraphErrors *gqp = (TGraphErrors*)arr->At(2);
 987   Double_t q, n = 0., entries;
 988   ax = h3->GetXaxis();
 989   for(Int_t ix=1; ix<=ax->GetNbins(); ix++){
 990     q = TMath::Exp(ax->GetBinCenter(ix));
 991     ax->SetRange(ix, ix);
 992     h1 = h3->Project3D("y");
 993     entries = h1->GetEntries();
 994     if(entries < 150) continue;
 995     h1->Fit(&f, "Q");
 996
 997     // Fill sy^2 = f(q)
 998     Int_t ip = gqm->GetN();
 999     gqm->SetPoint(ip, q, 1.e4*f.GetParameter(1));
1000     gqm->SetPointError(ip, 0., 1.e4*f.GetParError(1));
1001
1002     // correct sigma for ExB effect
1003     gqs->SetPoint(ip, q, 1.e4*f.GetParameter(2)/**f.GetParameter(2)-exb2*s2x)*/);
1004     gqs->SetPointError(ip, 0., 1.e4*f.GetParError(2)/**f.GetParameter(2)*/);
1005
1006     // save probability
1007     n += entries;
1008     gqp->SetPoint(ip, q, entries);
1009     gqp->SetPointError(ip, 0., 0./*TMath::Sqrt(entries)*/);
1010   }
1011
1012   // normalize probability and get mean sy
1013   Double_t sm = 0., sy;
1014   for(Int_t ip=gqp->GetN(); ip--;){
1015     gqp->GetPoint(ip, q, entries);
1016     entries/=n;
1017     gqp->SetPoint(ip, q, 1.e4*entries);
1018     gqs->GetPoint(ip, q, sy);
1019     sm += entries*sy;
1020   }
1021
1022   // error parametrization s(q) = <sy> + b(1/q-1/q0)
1023   TF1 fq("fq", "[0] + [1]/x", 20., 250.);
1024   gqs->Fit(&fq/*, "W"*/);
1025   printf("sm=%f [0]=%f [1]=%f\n", 1.e-4*sm, fq.GetParameter(0), fq.GetParameter(1));
1026   printf("  const Float_t sq0inv = %f; // [1/q0]\n", (sm-fq.GetParameter(0))/fq.GetParameter(1));
1027   printf("  const Float_t sqb    = %f; // [cm]\n", 1.e-4*fq.GetParameter(1));
1028 }
1029
1030 //_______________________________________________________
1031 Bool_t AliTRDclusterResolution::ProcessNormalTracks()
1032 {
1033 // Resolution as a function of y displacement from pad center and drift length.
1034 //
1035 // Since the error parameterization of cluster r-phi position can be written as (see AliTRDcluster::SetSigmaY2()):
1036 // BEGIN_LATEX
1037 // #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
1038 // END_LATEX
1039 // one can see that for phi = a_L one gets the following expression:
1040 // BEGIN_LATEX
1041 // #sigma_{y}^{2} = #sigma_{y}^{2}|_{B=0} + tg^{2}(#alpha_{L})*#sigma_{x}^{2}
1042 // END_LATEX
1043 // where we have explicitely marked the remaining term in case of absence of magnetic field. Thus one can use the
1044 // previous equation to estimate s_y for B=0 and than by comparing in magnetic field conditions one can get the s_x.
1045 // This is a simplified method to determine the error parameterization for s_x and s_y as compared to the one
1046 // implemented in ProcessSigma(). For more details on cluster error parameterization please see also
1047 // AliTRDcluster::SetSigmaY2()
1048 //
1049 // The representation of dy=f(y_cen, x_drift| layer) can be also used to estimate the systematic shift in the r-phi
1050 // coordinate resulting from imperfection in the cluster shape parameterization. From the expresion of the shift derived
1051 // in ProcessMean() with phi=exb one gets:
1052 // BEGIN_LATEX
1053 // <#Delta y>= <#delta x> * (tg(#alpha_{L})-h*dz/dx) + <#delta y - #delta x * tg(#alpha_{L})>
1054 // <#Delta y>(y_{cen})= -h*<#delta x>(x_{drift}, q_{cl}) * dz/dx + #delta y(y_{cen}, ...)
1055 // END_LATEX
1056 // where all dependences are made explicit. This last expression can be used in two ways:
1057 //   - by average on the dz/dx we can determine directly dy (the method implemented here)
1058 //   - by plotting as a function of dzdx one can determine both dx and dy components in an independent method.
1059 //Begin_Html
1060 //<img src="TRD/clusterYcorr.gif">
1061 //End_Html
1062 // Author
1063 // Alexandru Bercuci <A.Bercuci@gsi.de>
1064
1065   TObjArray *arr(NULL);
1066   TH3S *h3r(NULL), *h3t(NULL);
1067   if(!(arr= (TObjArray*)fContainer->At(kYRes))) {
1068     AliError("Missing resolution container");
1069     return kFALSE;
1070   }
1071   if(!(h3r = (TH3S*)arr->At(0))){
1072     AliError("Missing resolution pw/q histo");
1073     return kFALSE;
1074   } else if(!(Int_t)h3r->GetEntries()){
1075     AliError("Empty resolution pw/q histo");
1076     return kFALSE;
1077   }
1078   if(!(h3t = (TH3S*)arr->At(2))){
1079     AliError("Missing resolution t histo");
1080     return kFALSE;
1081   } else if(!(Int_t)h3t->GetEntries()){
1082     AliError("Empty resolution t histo");
1083     return kFALSE;
1084   }
1085
1086   // local variables
1087   Double_t x(0.), y(0.), ex(0.), ey(0.);
1088   Float_t tmin(6.5), tmax(20.5), tmed(0.5*(tmin+tmax));
1089   TGraphErrors *g[2]; TH1 *h(NULL);
1090   g[0] = new TGraphErrors();
1091   g[0]->SetMarkerStyle(24);g[0]->SetMarkerColor(kBlue);g[0]->SetLineColor(kBlue);
1092   g[1] = new TGraphErrors();
1093   g[1]->SetMarkerStyle(24);g[1]->SetMarkerColor(kRed);g[1]->SetLineColor(kRed);
1094
1095   // PROCESS RESOLUTION VS TB
1096   TF1 fsx("fsx", "[0]*[0]+[1]*[1]*[2]*0.1*(x-[3])", tmin, tmax);
1097   fsx.SetParNames("#sqrt{<#sigma^{2}(prf, q)>}(t_{med})", "D_{T}", "v_{drift}", "t_{med}");
1098   fsx.FixParameter(1, fDt);
1099   fsx.SetParameter(2, fVdrift);
1100   fsx.FixParameter(3, tmed);
1101   if(!AliTRDresolution::Process((TH2*)h3r->Project3D("yx"), g)) return kFALSE;
1102   for(Int_t ip(0); ip<g[1]->GetN(); ip++){
1103     g[1]->GetPoint(ip, x, y);ex = g[1]->GetErrorX(ip); ey = g[1]->GetErrorY(ip);
1104     g[1]->SetPoint(ip, x, y*y);g[1]->SetPointError(ip, ex, 2*y*ey);
1105   }
1106   g[1]->Fit(&fsx, "Q", "", tmin, tmax);
1107   if(fCanvas){
1108     g[1]->Draw("apl");
1109     fCanvas->Modified(); fCanvas->Update();
1110     h = g[1]->GetHistogram();
1111     h->SetTitle(fsx.GetTitle());
1112     h->GetXaxis()->SetTitle("t [tb]");h->GetXaxis()->CenterTitle();
1113     h->GetYaxis()->SetTitle("#sigma^{2} (y) [cm^{2}]");h->GetYaxis()->CenterTitle();
1114     if(IsSaveAs()) fCanvas->SaveAs(Form("D%03d_ResNormTrack_tb.gif", fDet));
1115     else gSystem->Sleep(100);
1116   }
1117
1118   // define model for resolution vs pw
1119   TF1 fg("fg", "gaus", -.5, .5); fg.FixParameter(1, 0.);
1120   TF1 fs("fs", "[0]*[0]*exp(-1*(x/[1])**2)+[2]", -.5, .5);
1121   fs.SetParNames("<#sigma^{max}(q,prf)>_{q}", "#sigma(pw)", "D_{T}^{2}*<x>");
1122   h3r->GetXaxis()->SetRange(tmin, tmax);
1123   if(!AliTRDresolution::Process((TH2*)h3r->Project3D("zy"), g, 200)) return kFALSE;
1124   for(Int_t ip(0); ip<g[1]->GetN(); ip++){
1125     g[1]->GetPoint(ip, x, y); ex = g[1]->GetErrorX(ip); ey = g[1]->GetErrorY(ip);
1126     g[1]->SetPoint(ip, x, y*y);g[1]->SetPointError(ip, ex, 2.*y*ey);
1127   }
1128   g[1]->Fit(&fg, "QR");
1129   fs.SetParameter(0, TMath::Sqrt(fg.GetParameter(0)));
1130   fs.SetParameter(1, fg.GetParameter(2));
1131   Float_t sdiff(fDt*fDt*fsx.GetParameter(2)*tmed*0.1);
1132   fs.SetParameter(2, sdiff);
1133   fs.SetParLimits(2, 0.1*sdiff, 1.9*sdiff);
1134   g[1]->Fit(&fs, "QR");
1135   if(fCanvas){
1136     g[1]->Draw("apl");
1137     fCanvas->Modified(); fCanvas->Update();
1138     h = g[1]->GetHistogram();
1139     h->SetTitle(fs.GetTitle());
1140     h->GetXaxis()->SetTitle("pw");h->GetXaxis()->CenterTitle();
1141     h->GetYaxis()->SetTitle("#sigma^{2} (y) [cm^{2}]");h->GetYaxis()->CenterTitle();
1142     if(IsSaveAs()) fCanvas->SaveAs(Form("D%03d_ResNormTrack_pw.gif", fDet));
1143     else gSystem->Sleep(100);
1144   }
1145
1146   AliDebug(2, Form("<s(q,prf)>[mum] = %7.3f", 1.e4*TMath::Sqrt(fsx.Eval(0.))));
1147   AliDebug(2, Form("<s(q)>[mum] = %7.3f", 1.e4*TMath::Sqrt(fs.Eval(-0.5)-fs.GetParameter(2))));
1148   AliDebug(2, Form("<s(x)>[mum] = %7.3f(prf) %7.3f(diff)", 1.e4*TMath::Sqrt(fs.GetParameter(2)), 1.e4*TMath::Sqrt(sdiff)));
1149
1150   // define model for resolution vs q
1151   TF1 fq("fq", "[0]*[0]*exp(-1*[1]*(x-[2])**2)+[2]", 2.5, 5.5);
1152   fq.SetParNames("<#sigma^{max}(q,prf)>_{prf}", "slope","mean", "D_{T}^{2}*<x>");
1153   if(!AliTRDresolution::Process((TH2*)h3t->Project3D("yx"), g)) return kFALSE;
1154   for(Int_t ip(0); ip<g[1]->GetN(); ip++){
1155     g[1]->GetPoint(ip, x, y); ex = g[1]->GetErrorX(ip); ey = g[1]->GetErrorY(ip);
1156     g[1]->SetPoint(ip, x, y*y);g[1]->SetPointError(ip, ex, 2.*y*ey);
1157   }
1158   fq.SetParameter(0, 8.e-2); fq.SetParLimits(0, 0., 1.);
1159   fq.SetParameter(1, 1.); //fq.SetParLimits(1, -1., 0.);
1160   fq.SetParameter(3, sdiff); fq.SetParLimits(3, 0.1*sdiff, 1.9*sdiff);
1161   g[1]->Fit(&fq, "QR");
1162 //   AliDebug(2, Form("<sq>[mum] = %7.3f", 1.e4*TMath::Sqrt(fs.Eval(-0.5)-fs.GetParameter(2)));
1163 //   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)));
1164   if(fCanvas){
1165     g[1]->Draw("apl");
1166     fCanvas->Modified(); fCanvas->Update();
1167     h = g[1]->GetHistogram();
1168     h->SetTitle(fs.GetTitle());
1169     h->GetXaxis()->SetTitle("log(q) [a.u.]");h->GetXaxis()->CenterTitle();
1170     h->GetYaxis()->SetTitle("#sigma^{2} (y) [cm^{2}]");h->GetYaxis()->CenterTitle();
1171     if(IsSaveAs()) fCanvas->SaveAs(Form("D%03d_ResNormTrack_q.gif", fDet));
1172     else gSystem->Sleep(100);
1173   }
1174   return kTRUE;
1175 }
1176
1177 //_______________________________________________________
1178 void AliTRDclusterResolution::ProcessSigma()
1179 {
1180 // As the r-phi coordinate is the only one which is measured by the TRD detector we have to rely on it to
1181 // estimate both the radial (x) and r-phi (y) errors. This method is based on the following assumptions.
1182 // The measured error in the y direction is the sum of the intrinsic contribution of the r-phi measurement
1183 // with the contribution of the radial measurement - because x is not a parameter of Alice track model (Kalman).
1184 // BEGIN_LATEX
1185 // #sigma^{2}|_{y} = #sigma^{2}_{y*} + #sigma^{2}_{x*}
1186 // END_LATEX
1187 // In the general case
1188 // BEGIN_LATEX
1189 // #sigma^{2}_{y*} = #sigma^{2}_{y} + tg^{2}(#alpha_{L})#sigma^{2}_{x_{drift}}
1190 // #sigma^{2}_{x*} = tg^{2}(#phi - #alpha_{L})*(#sigma^{2}_{x_{drift}} + #sigma^{2}_{x_{0}} + tg^{2}(#alpha_{L})*x^{2}/12)
1191 // END_LATEX
1192 // where we have explicitely show the lorentz angle correction on y and the projection of radial component on the y
1193 // direction through the track angle in the bending plane (phi). Also we have shown that the radial component in the
1194 // last equation has twp terms, the drift and the misalignment (x_0). For ideal geometry or known misalignment one
1195 // can solve the equation
1196 // BEGIN_LATEX
1197 // #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}]
1198 // END_LATEX
1199 // by fitting a straight line:
1200 // BEGIN_LATEX
1201 // #sigma^{2}|_{y} = a(x_{cl}, z_{cl}) * tg^{2}(#phi - #alpha_{L}) + b(x_{cl}, z_{cl})
1202 // END_LATEX
1203 // the error parameterization will be given by:
1204 // BEGIN_LATEX
1205 // #sigma_{x} (x_{cl}, z_{cl}) = #sqrt{a(x_{cl}, z_{cl}) - tg^{2}(#alpha_{L})*x^{2}/12}
1206 // #sigma_{y} (x_{cl}, z_{cl}) = #sqrt{b(x_{cl}, z_{cl}) - #sigma^{2}_{x} (x_{cl}, z_{cl}) * tg^{2}(#alpha_{L})}
1207 // END_LATEX
1208 // Below there is an example of such dependency.
1209 //Begin_Html
1210 //<img src="TRD/clusterSigmaMethod.gif">
1211 //End_Html
1212 //
1213 // The error parameterization obtained by this method are implemented in the functions AliTRDcluster::GetSX() and
1214 // AliTRDcluster::GetSYdrift(). For an independent method to determine s_y as a function of drift length check the
1215 // function ProcessCenterPad(). One has to keep in mind that while this method return the mean s_y over the distance
1216 // to pad center distribution the other method returns the *STANDARD* value at center=0 (maximum). To recover the
1217 // standard value one has to solve the obvious equation:
1218 // BEGIN_LATEX
1219 // #sigma_{y}^{STANDARD} = #frac{<#sigma_{y}>}{#int{s exp(s^{2}/#sigma) ds}}
1220 // END_LATEX
1221 // with "<s_y>" being the value calculated here and "sigma" the width of the s_y distribution calculated in
1222 // ProcessCenterPad().
1223 //
1224 // Author
1225 // Alexandru Bercuci <A.Bercuci@gsi.de>
1226
1227   TObjArray *arr = (TObjArray*)fContainer->At(kSigm);
1228   if(!arr){
1229     AliWarning("Missing dy=f(x_d, d_w) container");
1230     return;
1231   }
1232
1233   // init visualization
1234   TGraphErrors *ggs = NULL;
1235   TGraph *line = NULL;
1236   if(fCanvas){
1237     ggs = new TGraphErrors();
1238     line = new TGraph();
1239     line->SetLineColor(kRed);line->SetLineWidth(2);
1240   }
1241
1242   // init logistic support
1243   TF1 f("f", "gaus", -.5, .5);
1244   TLinearFitter gs(1,"pol1");
1245   TH1 *hFrame=NULL;
1246   TH1D *h1 = NULL; TH3S *h3=NULL;
1247   TAxis *ax = NULL;
1248   Double_t exb2 = fExB*fExB;
1249   AliTRDcluster c;
1250   TTree *t = (TTree*)fResults->At(kSigm);
1251   for(Int_t ix=0; ix<AliTRDseedV1::kNtb; ix++){
1252     if(!(h3=(TH3S*)arr->At(ix))) continue;
1253     c.SetPadTime(ix);
1254     fX = c.GetXloc(fT0, fVdrift);
1255     fT = c.GetLocalTimeBin(); // ideal
1256     printf(" pad time[%d] local[%f]\n", ix, fT);
1257     for(Int_t iz=1; iz<=h3->GetXaxis()->GetNbins(); iz++){
1258       ax = h3->GetXaxis();
1259       ax->SetRange(iz, iz);
1260       fZ = ax->GetBinCenter(iz);
1261
1262       // reset visualization
1263       if(fCanvas){
1264         new(ggs) TGraphErrors();
1265         ggs->SetMarkerStyle(7);
1266       }
1267       gs.ClearPoints();
1268
1269       for(Int_t ip=1; ip<=h3->GetYaxis()->GetNbins(); ip++){
1270         ax = h3->GetYaxis();
1271         ax->SetRange(ip, ip);
1272         Double_t tgl = ax->GetBinCenter(ip);
1273         // finish navigation in the HnSparse
1274
1275         //if(TMath::Abs(dydx)>0.18) continue;
1276         Double_t tgg = (tgl-fExB)/(1.+tgl*fExB);
1277         Double_t tgg2 = tgg*tgg;
1278
1279         h1 = (TH1D*)h3->Project3D("z");
1280         Int_t entries = (Int_t)h1->Integral();
1281         if(entries < 50) continue;
1282         //Adjust(&f, h1);
1283         h1->Fit(&f, "QN");
1284
1285         Double_t s2  = f.GetParameter(2)*f.GetParameter(2);
1286         Double_t s2e = 2.*f.GetParameter(2)*f.GetParError(2);
1287         // Fill sy^2 = f(tg^2(phi-a_L))
1288         gs.AddPoint(&tgg2, s2, s2e);
1289
1290         if(!ggs) continue;
1291         Int_t jp = ggs->GetN();
1292         ggs->SetPoint(jp, tgg2, s2);
1293         ggs->SetPointError(jp, 0., s2e);
1294       }
1295       // TODO here a more robust fit method has to be provided
1296       // for which lower boundaries on the parameters have to
1297       // be imposed. Unfortunately the Minuit fit does not work
1298       // for the TGraph in the case of B not 0.
1299       if(gs.Eval()) continue;
1300
1301       fR[0] = gs.GetParameter(1) - fX*fX*exb2/12.;
1302       AliDebug(3, Form("  s2x+x2=%f ang=%f s2x=%f", gs.GetParameter(1), fX*fX*exb2/12., fR[0]));
1303       fR[0] = TMath::Max(fR[0], Float_t(4.e-4));
1304
1305       // s^2_y  = s0^2_y + tg^2(a_L) * s^2_x
1306       // s0^2_y = f(D_L)*x + s_PRF^2
1307       fR[2]= gs.GetParameter(0)-exb2*fR[0];
1308       AliDebug(3, Form("  s2y+s2x=%f s2y=%f", fR[0], fR[2]));
1309       fR[2] = TMath::Max(fR[2], Float_t(2.5e-5));
1310       fR[0] = TMath::Sqrt(fR[0]);
1311       fR[1] = .5*gs.GetParError(1)/fR[0];
1312       fR[2] = TMath::Sqrt(fR[2]);
1313       fR[3] = gs.GetParError(0)+exb2*exb2*gs.GetParError(1);
1314       t->Fill();
1315       AliDebug(2, Form("xd=%4.2f[cm] sx=%6.1f[um] sy=%5.1f[um]", fX, 1.e4*fR[0], 1.e4*fR[2]));
1316
1317       if(!fCanvas) continue;
1318       fCanvas->cd(); fCanvas->SetLogx(); //fCanvas->SetLogy();
1319       if(!hFrame){
1320         fCanvas->SetMargin(0.15, 0.01, 0.1, 0.01);
1321         hFrame=new TH1I("hFrame", "", 100, 0., .3);
1322         hFrame->SetMinimum(0.);hFrame->SetMaximum(.005);
1323         hFrame->SetXTitle("tg^{2}(#phi-#alpha_{L})");
1324         hFrame->SetYTitle("#sigma^{2}y[cm^{2}]");
1325         hFrame->GetYaxis()->SetTitleOffset(2.);
1326         hFrame->SetLineColor(1);hFrame->SetLineWidth(1);
1327         hFrame->Draw();
1328       } else hFrame->Reset();
1329       Double_t xx = 0., dxx=.2/50;
1330       for(Int_t ip=0;ip<50;ip++){
1331         line->SetPoint(ip, xx,  gs.GetParameter(0)+xx*gs.GetParameter(1));
1332         xx+=dxx;
1333       }
1334       ggs->Draw("pl"); line->Draw("l");
1335       fCanvas->Modified(); fCanvas->Update();
1336       if(IsSaveAs()) fCanvas->SaveAs(Form("Figures/ProcessSigma_z[%5.3f]_x[%5.3f].gif", fZ, fX));
1337       else gSystem->Sleep(100);
1338     }
1339   }
1340   return;
1341 }
1342
1343 //_______________________________________________________
1344 void AliTRDclusterResolution::ProcessMean()
1345 {
1346 // By this method the cluster shift in r-phi and radial directions can be estimated by comparing with the MC.
1347 // The resolution of the cluster corrected for pad tilt with respect to MC in the r-phi (measuring) plane can be
1348 // expressed by:
1349 // BEGIN_LATEX
1350 // #Delta y=w - y_{MC}(x_{cl})
1351 // w = y_{cl}^{'} + h*(z_{MC}(x_{cl})-z_{cl})
1352 // y_{MC}(x_{cl}) = y_{0} - dy/dx*x_{cl}
1353 // z_{MC}(x_{cl}) = z_{0} - dz/dx*x_{cl}
1354 // y_{cl}^{'} = y_{cl}-x_{cl}*tg(#alpha_{L})
1355 // END_LATEX
1356 // where x_cl is the drift length attached to a cluster, y_cl is the r-phi coordinate of the cluster measured by
1357 // charge sharing on adjacent pads and y_0 and z_0 are MC reference points (as example the track references at
1358 // entrance/exit of a chamber). If we suppose that both r-phi (y) and radial (x) coordinate of the clusters are
1359 // affected by errors we can write
1360 // BEGIN_LATEX
1361 // x_{cl} = x_{cl}^{*} + #delta x
1362 // y_{cl} = y_{cl}^{*} + #delta y
1363 // END_LATEX
1364 // where the starred components are the corrected values. Thus by definition the following quantity
1365 // BEGIN_LATEX
1366 // #Delta y^{*}= w^{*} - y_{MC}(x_{cl}^{*})
1367 // END_LATEX
1368 // has 0 average over all dependency. Using this decomposition we can write:
1369 // BEGIN_LATEX
1370 // <#Delta y>=<#Delta y^{*}> + <#delta x * (dy/dx-h*dz/dx) + #delta y - #delta x * tg(#alpha_{L})>
1371 // END_LATEX
1372 // which can be transformed to the following linear dependence:
1373 // BEGIN_LATEX
1374 // <#Delta y>= <#delta x> * (dy/dx-h*dz/dx) + <#delta y - #delta x * tg(#alpha_{L})>
1375 // END_LATEX
1376 // if expressed as function of dy/dx-h*dz/dx. Furtheremore this expression can be plotted for various clusters
1377 // i.e. we can explicitely introduce the diffusion (x_cl) and drift cell - anisochronity (z_cl) dependences. From
1378 // plotting this dependence and linear fitting it with:
1379 // BEGIN_LATEX
1380 // <#Delta y>= a(x_{cl}, z_{cl}) * (dy/dx-h*dz/dx) + b(x_{cl}, z_{cl})
1381 // END_LATEX
1382 // the systematic shifts will be given by:
1383 // BEGIN_LATEX
1384 // #delta x (x_{cl}, z_{cl}) = a(x_{cl}, z_{cl})
1385 // #delta y (x_{cl}, z_{cl}) = b(x_{cl}, z_{cl}) + a(x_{cl}, z_{cl}) * tg(#alpha_{L})
1386 // END_LATEX
1387 // Below there is an example of such dependency.
1388 //Begin_Html
1389 //<img src="TRD/clusterShiftMethod.gif">
1390 //End_Html
1391 //
1392 // The occurance of the radial shift is due to the following conditions
1393 //   - the approximation of a constant drift velocity over the drift length (larger drift velocities close to
1394 //     cathode wire plane)
1395 //   - the superposition of charge tails in the amplification region (first clusters appear to be located at the
1396 //     anode wire)
1397 //   - the superposition of charge tails in the drift region (shift towards anode wire)
1398 //   - diffusion effects which convolute with the TRF thus enlarging it
1399 //   - approximate knowledge of the TRF (approximate measuring in test beam conditions)
1400 //
1401 // The occurance of the r-phi shift is due to the following conditions
1402 //   - approximate model for cluster shape (LUT)
1403 //   - rounding-up problems
1404 //
1405 // The numerical results for ideal simulations for the radial and r-phi shifts are displayed below and used
1406 // for the cluster reconstruction (see the functions AliTRDcluster::GetXcorr() and AliTRDcluster::GetYcorr()).
1407 //Begin_Html
1408 //<img src="TRD/clusterShiftX.gif">
1409 //<img src="TRD/clusterShiftY.gif">
1410 //End_Html
1411 // More details can be found in the presentation given during the TRD
1412 // software meeting at the end of 2008 and beginning of year 2009, published on indico.cern.ch.
1413 //
1414 // Author
1415 // Alexandru Bercuci <A.Bercuci@gsi.de>
1416
1417
1418
1419   TObjArray *arr = (TObjArray*)fContainer->At(kMean);
1420   if(!arr){
1421     AliWarning("Missing dy=f(x_d, d_w) container");
1422     return;
1423   }
1424
1425   // init logistic support
1426   TF1 f("f", "gaus", -.5, .5);
1427   TF1 line("l", "[0]+[1]*x", -.15, .15);
1428   TGraphErrors *gm = new TGraphErrors();
1429   TH1 *hFrame=NULL;
1430   TH1D *h1 = NULL; TH3S *h3 =NULL;
1431   TAxis *ax = NULL;
1432
1433   AliDebug(1, Form("Calibrate for Det[%3d] t0[%5.3f] vd[%5.3f]", fDet, fT0, fVdrift));
1434
1435   AliTRDcluster c;
1436   TTree *t = (TTree*)fResults->At(kMean);
1437   for(Int_t ix=0; ix<AliTRDseedV1::kNtb; ix++){
1438     if(!(h3=(TH3S*)arr->At(ix))) continue;
1439     c.SetPadTime(ix);
1440     fX = c.GetXloc(fT0, fVdrift);
1441     fT = c.GetLocalTimeBin();
1442     for(Int_t iz=1; iz<=h3->GetXaxis()->GetNbins(); iz++){
1443       ax = h3->GetXaxis();
1444       ax->SetRange(iz, iz);
1445       fZ = ax->GetBinCenter(iz);
1446
1447       // reset fitter
1448       new(gm) TGraphErrors();
1449       gm->SetMarkerStyle(7);
1450
1451       for(Int_t ip=1; ip<=h3->GetYaxis()->GetNbins(); ip++){
1452         ax = h3->GetYaxis();
1453         ax->SetRange(ip, ip);
1454         Double_t tgl = ax->GetBinCenter(ip);
1455         // finish navigation in the HnSparse
1456
1457         h1 = (TH1D*)h3->Project3D("z");
1458         Int_t entries = (Int_t)h1->Integral();
1459         if(entries < 50) continue;
1460         //Adjust(&f, h1);
1461         h1->Fit(&f, "QN");
1462
1463         // Fill <Dy> = f(dydx - h*dzdx)
1464         Int_t jp = gm->GetN();
1465         gm->SetPoint(jp, tgl, f.GetParameter(1));
1466         gm->SetPointError(jp, 0., f.GetParError(1));
1467       }
1468       if(gm->GetN()<10) continue;
1469
1470       gm->Fit(&line, "QN");
1471       fR[0] = line.GetParameter(1); // dx
1472       fR[1] = line.GetParError(1);
1473       fR[2] = line.GetParameter(0) + fExB*fR[0]; // xs = dy - tg(a_L)*dx
1474       t->Fill();
1475       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]));
1476       if(!fCanvas) continue;
1477
1478       fCanvas->cd();
1479       if(!hFrame){
1480         fCanvas->SetMargin(0.1, 0.02, 0.1, 0.01);
1481         hFrame=new TH1I("hFrame", "", 100, -.3, .3);
1482         hFrame->SetMinimum(-.1);hFrame->SetMaximum(.1);
1483         hFrame->SetXTitle("tg#phi-htg#theta");
1484         hFrame->SetYTitle("#Delta y[cm]");
1485         hFrame->GetYaxis()->SetTitleOffset(1.5);
1486         hFrame->SetLineColor(1);hFrame->SetLineWidth(1);
1487         hFrame->Draw();
1488       } else hFrame->Reset();
1489       gm->Draw("pl"); line.Draw("same");
1490       fCanvas->Modified(); fCanvas->Update();
1491       if(IsSaveAs()) fCanvas->SaveAs(Form("Figures/ProcessMean_Z[%5.3f]_TB[%02d].gif", fZ, ix));
1492       else gSystem->Sleep(100);
1493     }
1494   }
1495 }