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