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