TRD/qaRec/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 "info/AliTRDclusterInfo.h"
 172 #include "AliTRDgeometry.h"
 173 #include "AliTRDcalibDB.h"
 174 #include "AliTRDCommonParam.h"
 175 #include "Cal/AliTRDCalROC.h"
 176 #include "Cal/AliTRDCalDet.h"
 177
 178 #include "AliLog.h"
 179 #include "AliTracker.h"
 180 #include "AliCDBManager.h"
 181
 182 #include "TROOT.h"
 183 #include "TObjArray.h"
 184 #include "TAxis.h"
 185 #include "TF1.h"
 186 #include "TLegend.h"
 187 #include "TGraphErrors.h"
 188 #include "TLine.h"
 189 #include "TH2I.h"
 190 #include "TH3S.h"
 191 #include "TTree.h"
 192 #include "TMath.h"
 193 #include "TLinearFitter.h"
 194
 195 #include "TCanvas.h"
 196 #include "TSystem.h"
 197
 198 ClassImp(AliTRDclusterResolution)
 199
 200 const Float_t AliTRDclusterResolution::fgkTimeBinLength = 1./ AliTRDCommonParam::Instance()->GetSamplingFrequency();
 201 //_______________________________________________________
 202 AliTRDclusterResolution::AliTRDclusterResolution(const char *name, const char *title)
 203   : AliTRDrecoTask(name, title)
 204   ,fCanvas(0x0)
 205   ,fInfo(0x0)
 206   ,fResults(0x0)
 207   ,fAt(0x0)
 208   ,fStatus(0)
 209   ,fDet(-1)
 210   ,fExB(0.)
 211   ,fVdrift(0.)
 212   ////////////
 213   ,fLy(0)
 214   ,fX(0.)
 215   ,fY(0.)
 216   ,fZ(0.)
 217 {
 218   memset(fR, 0, 4*sizeof(Float_t));
 219   memset(fP, 0, 4*sizeof(Float_t));
 220   // time drift axis
 221   fAt = new TAxis(kNTB, 0., kNTB*fgkTimeBinLength);
 222
 223   // By default register all analysis
 224   // The user can switch them off in his steering macro
 225   SetProcess(kQRes);
 226   SetProcess(kCenter);
 227   SetProcess(kMean);
 228   SetProcess(kSigm);
 229 }
 230
 231 //_______________________________________________________
 232 AliTRDclusterResolution::~AliTRDclusterResolution()
 233 {
 234   if(fCanvas) delete fCanvas;
 235   if(fAt) delete fAt;
 236   if(fResults){
 237     fResults->Delete();
 238     delete fResults;
 239   }
 240 }
 241
 242 //_______________________________________________________
 243 void AliTRDclusterResolution::ConnectInputData(Option_t *)
 244 {
 245   fInfo = dynamic_cast<TObjArray *>(GetInputData(0));
 246 }
 247
 248 //_______________________________________________________
 249 void AliTRDclusterResolution::CreateOutputObjects()
 250 {
 251   OpenFile(0, "RECREATE");
 252   fContainer = Histos();
 253 }
 254
 255 //_______________________________________________________
 256 Bool_t AliTRDclusterResolution::GetRefFigure(Int_t ifig)
 257 {
 258   if(!fResults) return kFALSE;
 259   TLegend *leg = 0x0;
 260   TList *l = 0x0;
 261   TObjArray *arr = 0x0;
 262   TH2 *h2 = 0x0;TH1 *h1 = 0x0;
 263   TGraphErrors *gm(0x0), *gs(0x0), *gp(0x0);
 264   switch(ifig){
 265   case kQRes:
 266     if(!(arr = (TObjArray*)fResults->At(kQRes))) break;
 267     if(!(gm = (TGraphErrors*)arr->At(0))) break;
 268     if(!(gs = (TGraphErrors*)arr->At(1))) break;
 269     if(!(gp = (TGraphErrors*)arr->At(2))) break;
 270     gs->Draw("apl");
 271     gs->GetHistogram()->GetYaxis()->SetRangeUser(-50., 700.);
 272     gs->GetHistogram()->SetXTitle("Q [a.u.]");
 273     gs->GetHistogram()->SetYTitle("#sigma_{y} / #mu_{y} [#mum] / freq");
 274     gm->Draw("pl");
 275     gp->Draw("pl");
 276     return kTRUE;
 277   case kCenter:
 278     if(!(arr = (TObjArray*)fResults->At(kCenter))) break;
 279     gPad->Divide(2, 1); l = gPad->GetListOfPrimitives();
 280     ((TVirtualPad*)l->At(0))->cd();
 281     ((TTree*)arr->At(0))->Draw("y:x>>h(23, 0.1, 2.4, 51, -.51, .51)",
 282             "m[0]*(ly==0&&abs(m[0])<1.e-1)", "colz");
 283     ((TVirtualPad*)l->At(1))->cd();
 284     leg= new TLegend(.7, .7, .9, .95);
 285     leg->SetBorderSize(0); leg->SetFillColor(0); leg->SetFillStyle(0);
 286     leg->SetHeader("TRD Plane");
 287     for(Int_t il = 1; il<=AliTRDgeometry::kNlayer; il++){
 288       if(!(gm = (TGraphErrors*)arr->At(il))) return kFALSE;
 289       gm->Draw(il>1?"pc":"apc"); leg->AddEntry(gm, Form("%d", il-1), "pl");
 290       if(il>1) continue;
 291       gm->GetHistogram()->SetXTitle("t_{drift} [#mus]");
 292       gm->GetHistogram()->SetYTitle("#sigma_{y}(x|cen=0) [#mum]");
 293       gm->GetHistogram()->GetYaxis()->SetRangeUser(150., 500.);
 294     }
 295     leg->Draw();
 296     return kTRUE;
 297   case kSigm:
 298     if(!(arr = (TObjArray*)fResults->At(kSigm))) break;
 299     gPad->Divide(2, 1); l = gPad->GetListOfPrimitives();
 300     if(!(h2 = (TH2F*)arr->At(0))) return kFALSE;
 301     ((TVirtualPad*)l->At(0))->cd();
 302     h1 = h2->ProjectionY("hsx_pyy"); h1->Scale(1.e4/kND); h1->SetMarkerStyle(24);
 303     h1->SetYTitle("<#sigma_{x}> [#mum]");
 304     h1->GetXaxis()->SetRange(2, kNTB-1); h1->Draw("pc");
 305
 306     if(!(h2 = (TH2F*)arr->At(1))) return kFALSE;
 307     ((TVirtualPad*)l->At(1))->cd();
 308     h1 = h2->ProjectionY("hsy_pyy"); h1->Scale(1.e4/kND); h1->SetMarkerStyle(24);
 309     h1->SetYTitle("<#sigma_{y}> [#mum]");
 310     h1->GetXaxis()->SetRange(2, kNTB-1); h1->Draw("pc");
 311     return kTRUE;
 312   case kMean:
 313     if(!(arr = (TObjArray*)fResults->At(kMean))) break;
 314     gPad->Divide(2, 1);  l = gPad->GetListOfPrimitives();
 315     ((TVirtualPad*)l->At(0))->cd();
 316     if(!(gm = (TGraphErrors*)arr->At(0))) return kFALSE;
 317     gm->Draw("apl");
 318     gm->GetHistogram()->SetXTitle("t_{drift} [#mus]");
 319     gm->GetHistogram()->SetYTitle("dx [#mum]");
 320
 321     ((TVirtualPad*)l->At(1))->cd();
 322     if(!(gm = (TGraphErrors*)arr->At(1))) return kFALSE;
 323     gm->Draw("apl");
 324     gm->GetHistogram()->SetXTitle("t_{drift} [#mus]");
 325     gm->GetHistogram()->SetYTitle("dy [#mum]");
 326
 327     return kTRUE;
 328   default:
 329     break;
 330   }
 331   AliWarning("No container/data found.");
 332   return kFALSE;
 333 }
 334
 335 //_______________________________________________________
 336 TObjArray* AliTRDclusterResolution::Histos()
 337 {
 338   if(fContainer) return fContainer;
 339   fContainer = new TObjArray(kNtasks);
 340   //fContainer->SetOwner(kTRUE);
 341
 342   TH3S *h3 = 0x0;
 343   TObjArray *arr = 0x0;
 344
 345   fContainer->AddAt(arr = new TObjArray(2*AliTRDgeometry::kNlayer), kCenter);
 346   arr->SetName("Center");
 347   for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
 348     // add resolution plot for each layer
 349     if(!(h3=(TH3S*)gROOT->FindObject(Form("hCenResLy%d", il)))){
 350       h3 = new TH3S(
 351         Form("hCenResLy%d", il),
 352         Form(" ly [%d]", il),
 353         kNTB, fAt->GetBinLowEdge(1), fAt->GetBinUpEdge(kNTB),   // x
 354         51, -.51, .51, // y
 355         60, -.3, .3); // dy
 356       h3->SetXTitle("x [#mus]");
 357       h3->SetYTitle("y [pw]");
 358       h3->SetZTitle("#Delta y[cm]");
 359     } h3->Reset();
 360     arr->AddAt(h3, il);
 361     // add Pull plot for each layer
 362     if(!(h3=(TH3S*)gROOT->FindObject(Form("hCenPullLy%d", il)))){
 363       h3 = new TH3S(
 364         Form("hCenPullLy%d", il),
 365         Form(" ly [%d]", il),
 366         kNTB, fAt->GetBinLowEdge(1), fAt->GetBinUpEdge(kNTB),   // x
 367         51, -.51, .51, // y
 368         60, -4., 4.); // dy
 369       h3->SetXTitle("x [#mus]");
 370       h3->SetYTitle("y [pw]");
 371       h3->SetZTitle("#Delta y/#sigma_{y}");
 372     } h3->Reset();
 373     arr->AddAt(h3, AliTRDgeometry::kNlayer+il);
 374   }
 375
 376   if(!(h3 = (TH3S*)gROOT->FindObject("Charge"))){
 377     h3 = new TH3S("Charge", "dy=f(q)", 50, 2.2, 7.5, 60, -.3, .3, 60, -4., 4.);
 378     h3->SetXTitle("log(q) [a.u.]");
 379     h3->SetYTitle("#Delta y[cm]");
 380     h3->SetZTitle("#Delta y/#sigma_{y}");
 381   }
 382   fContainer->AddAt(h3, kQRes);
 383
 384   fContainer->AddAt(arr = new TObjArray(kNTB), kSigm);
 385   arr->SetName("Resolution");
 386   for(Int_t ix=0; ix<kNTB; ix++){
 387     if(!(h3=(TH3S*)gROOT->FindObject(Form("hr_x%02d", ix)))){
 388       h3 = new TH3S(
 389         Form("hr_x%02d", ix),
 390         Form(" t_{drift}(%3.1f-%3.1f)[#mus]", fAt->GetBinLowEdge(ix+1), fAt->GetBinUpEdge(ix+1)),
 391         kND, 0., 2.5,   // z
 392         35, -.35, .35, // tgp
 393         60, -.3, .3); // dy
 394       h3->SetXTitle("z [mm]");
 395       h3->SetYTitle("tg#phi");
 396       h3->SetZTitle("#Delta y[cm]");
 397     }
 398     arr->AddAt(h3, ix);
 399   }
 400
 401   fContainer->AddAt(arr = new TObjArray(kNTB), kMean);
 402   arr->SetName("Systematics");
 403   for(Int_t ix=0; ix<kNTB; ix++){
 404     if(!(h3=(TH3S*)gROOT->FindObject(Form("hs_x%02d", ix)))){
 405       h3 = new TH3S(
 406         Form("hs_x%02d", ix),
 407         Form(" t_{drift}(%3.1f-%3.1f)[#mus]", fAt->GetBinLowEdge(ix+1), fAt->GetBinUpEdge(ix+1)),
 408         kND, 0., 2.5,   // z
 409         35, -.35, .35, // tgp-h tgt
 410         60, -.3, .3); // dy
 411       h3->SetXTitle("z [mm]");
 412       h3->SetYTitle("tg(#phi) - h*tg(#theta)");
 413       h3->SetZTitle("#Delta y[cm]");
 414     }
 415     arr->AddAt(h3, ix);
 416   }
 417
 418   return fContainer;
 419 }
 420
 421 //_______________________________________________________
 422 void AliTRDclusterResolution::Exec(Option_t *)
 423 {
 424   if(!HasExB()) AliWarning("ExB was not set. Call SetExB() before running the task.");
 425
 426   Int_t det, t;
 427   Float_t x, y, z, q, dy, dydx, dzdx, cov[3], covcl[3];
 428   TH3S *h3 = 0x0;
 429
 430   // define limits around ExB for which x contribution is negligible
 431   const Float_t kDtgPhi = 3.5e-2; //(+- 2 deg)
 432
 433   TObjArray *arr0 = (TObjArray*)fContainer->At(kCenter);
 434   TObjArray *arr1 = (TObjArray*)fContainer->At(kSigm);
 435   TObjArray *arr2 = (TObjArray*)fContainer->At(kMean);
 436
 437   const AliTRDclusterInfo *cli = 0x0;
 438   TIterator *iter=fInfo->MakeIterator();
 439   while((cli=dynamic_cast<AliTRDclusterInfo*>((*iter)()))){
 440     cli->GetCluster(det, x, y, z, q, t, covcl);
 441     if(fDet>=0 && fDet!=det) continue;
 442
 443     dy = cli->GetResolution();
 444     cli->GetGlobalPosition(y, z, dydx, dzdx, &cov[0]);
 445
 446     // resolution as a function of cluster charge
 447     // only for phi equal exB
 448     if(TMath::Abs(dydx-fExB) < kDtgPhi){
 449       h3 = (TH3S*)fContainer->At(kQRes);
 450       h3->Fill(TMath::Log(q), dy, dy/TMath::Sqrt(covcl[0]));
 451     }
 452
 453     // do not use problematic clusters in resolution analysis
 454     // TODO define limits as calibration aware (gain) !!
 455     if(q<20. || q>250.) continue;
 456
 457     x = (t+.5)*fgkTimeBinLength; // conservative approach !!
 458
 459     // resolution as a function of y displacement from pad center
 460     // only for phi equal exB
 461     if(TMath::Abs(dydx-fExB) < kDtgPhi/* &&
 462        TMath::Abs(x-0.675)<0.225*/){
 463       Int_t ly(AliTRDgeometry::GetLayer(det));
 464       h3 = (TH3S*)arr0->At(ly);
 465       h3->Fill(x, cli->GetYDisplacement(), dy);
 466       h3 = (TH3S*)arr0->At(AliTRDgeometry::kNlayer+ly);
 467       h3->Fill(x, cli->GetYDisplacement(), dy/TMath::Sqrt(covcl[0]));
 468     }
 469
 470     Int_t ix = fAt->FindBin(x);
 471     if(ix==0 || ix == fAt->GetNbins()+1){
 472       AliWarning(Form("Drift time %3.1f outside allowed range", x));
 473       continue;
 474     }
 475
 476     // fill histo for resolution (sigma)
 477     ((TH3S*)arr1->At(ix-1))->Fill(10.*cli->GetAnisochronity(), dydx, dy);
 478
 479     // fill histo for systematic (mean)
 480     ((TH3S*)arr2->At(ix-1))->Fill(10.*cli->GetAnisochronity(), dydx-cli->GetTilt()*dzdx, dy);
 481   }
 482   PostData(0, fContainer);
 483 }
 484
 485
 486 //_______________________________________________________
 487 Bool_t AliTRDclusterResolution::PostProcess()
 488 {
 489   if(!fContainer) return kFALSE;
 490   if(!HasExB()) AliWarning("ExB was not set. Call SetExB() before running the post processing.");
 491
 492   TObjArray *arr = 0x0;
 493   TTree *t=0x0;
 494   if(!fResults){
 495     TGraphErrors *g = 0x0;
 496     fResults = new TObjArray(kNtasks);
 497     fResults->SetOwner();
 498     fResults->AddAt(arr = new TObjArray(3), kQRes);
 499     arr->SetOwner();
 500     arr->AddAt(g = new TGraphErrors(), 0);
 501     g->SetLineColor(kBlue); g->SetMarkerColor(kBlue);
 502     g->SetMarkerStyle(7);
 503     arr->AddAt(g = new TGraphErrors(), 1);
 504     g->SetLineColor(kRed); g->SetMarkerColor(kRed);
 505     g->SetMarkerStyle(23);
 506     arr->AddAt(g = new TGraphErrors(), 2);
 507     g->SetLineColor(kGreen); g->SetMarkerColor(kGreen);
 508     g->SetMarkerStyle(7);
 509
 510     // pad center dependence
 511     fResults->AddAt(arr = new TObjArray(AliTRDgeometry::kNlayer+1), kCenter);
 512     arr->SetOwner();
 513     arr->AddAt(
 514     t = new TTree("cent", "dy=f(y,x,ly)"), 0);
 515     t->Branch("ly", &fLy, "ly/B");
 516     t->Branch("x", &fX, "x/F");
 517     t->Branch("y", &fY, "y/F");
 518     t->Branch("m", &fR[0], "m[2]/F");
 519     t->Branch("s", &fR[2], "s[2]/F");
 520     t->Branch("pm", &fP[0], "pm[2]/F");
 521     t->Branch("ps", &fP[2], "ps[2]/F");
 522     for(Int_t il=1; il<=AliTRDgeometry::kNlayer; il++){
 523       arr->AddAt(g = new TGraphErrors(), il);
 524       g->SetLineColor(il); g->SetLineStyle(il);
 525       g->SetMarkerColor(il);g->SetMarkerStyle(4);
 526     }
 527
 528
 529     fResults->AddAt(t = new TTree("sigm", "dy=f(dw,x,dydx)"), kSigm);
 530     t->Branch("x", &fX, "x/F");
 531     t->Branch("z", &fZ, "z/F");
 532     t->Branch("sx", &fR[0], "sx[2]/F");
 533     t->Branch("sy", &fR[2], "sy[2]/F");
 534
 535
 536     fResults->AddAt(t = new TTree("mean", "dy=f(dw,x,dydx - h dzdx)"), kMean);
 537     t->Branch("x", &fX, "x/F");
 538     t->Branch("z", &fZ, "z/F");
 539     t->Branch("dx", &fR[0], "dx[2]/F");
 540     t->Branch("dy", &fR[2], "dy[2]/F");
 541   } else {
 542     TObject *o = 0x0;
 543     TIterator *iter=fResults->MakeIterator();
 544     while((o=((*iter)()))) o->Clear(); // maybe it is wrong but we should never reach this point
 545   }
 546
 547   // precalculated value of tg^2(alpha_L)
 548   Double_t exb2 = fExB*fExB;
 549   // square of the mean value of sigma drift length.
 550   // has to come from previous calibration
 551   //Double_t sxd2 = 1.;// [mm^2]
 552
 553   printf("ExB[%e] ExB2[%e]\n", fExB, exb2);
 554
 555   // process resolution dependency on charge
 556   if(HasProcess(kQRes)) ProcessCharge();
 557
 558   // process resolution dependency on y displacement
 559   if(HasProcess(kCenter)) ProcessCenterPad();
 560
 561   // process resolution dependency on drift legth and drift cell width
 562   if(HasProcess(kSigm)) ProcessSigma();
 563
 564   // process systematic shift on drift legth and drift cell width
 565   if(HasProcess(kMean)) ProcessMean();
 566
 567   return kTRUE;
 568 }
 569
 570 //_______________________________________________________
 571 Bool_t AliTRDclusterResolution::SetExB(Int_t det, Int_t col, Int_t row)
 572 {
 573   // check OCDB
 574   AliCDBManager *cdb = AliCDBManager::Instance();
 575   if(cdb->GetRun() < 0){
 576     AliError("OCDB manager not properly initialized");
 577     return kFALSE;
 578   }
 579
 580   // check magnetic field
 581   if(TMath::Abs(AliTracker::GetBz()) < 1.e-10){
 582     AliWarning("B=0. Magnetic field may not be initialized. Continue if you know what you are doing !");
 583   }
 584
 585   // set reference detector if any
 586   if(det>=0 && det<AliTRDgeometry::kNdet) fDet = det;
 587   else det = 0;
 588
 589   AliTRDcalibDB *fCalibration  = AliTRDcalibDB::Instance();
 590   AliTRDCalROC  *fCalVdriftROC = fCalibration->GetVdriftROC(det);
 591   const AliTRDCalDet  *fCalVdriftDet = fCalibration->GetVdriftDet();
 592
 593   fVdrift = fCalVdriftDet->GetValue(det) * fCalVdriftROC->GetValue(col, row);
 594   fExB   = AliTRDCommonParam::Instance()->GetOmegaTau(fVdrift);
 595   SetBit(kExB);
 596   return kTRUE;
 597 }
 598
 599 //_______________________________________________________
 600 void AliTRDclusterResolution::SetVisual()
 601 {
 602   if(fCanvas) return;
 603   fCanvas = new TCanvas("clResCanvas", "Cluster Resolution Visualization", 10, 10, 600, 600);
 604 }
 605
 606 //_______________________________________________________
 607 void AliTRDclusterResolution::ProcessCharge()
 608 {
 609   TH2I *h2 = 0x0;
 610   if(!(h2 = (TH2I*)fContainer->At(kQRes))) {
 611     AliWarning("Missing dy=f(Q) histo");
 612     return;
 613   }
 614   TF1 f("f", "gaus", -.5, .5);
 615   TAxis *ax = 0x0;
 616   TH1D *h1 = 0x0;
 617
 618   TObjArray *arr = (TObjArray*)fResults->At(kQRes);
 619   TGraphErrors *gqm = (TGraphErrors*)arr->At(0);
 620   TGraphErrors *gqs = (TGraphErrors*)arr->At(1);
 621   TGraphErrors *gqp = (TGraphErrors*)arr->At(2);
 622   Double_t q, n = 0., entries;
 623   ax = h2->GetXaxis();
 624   for(Int_t ix=1; ix<=ax->GetNbins(); ix++){
 625     q = TMath::Exp(ax->GetBinCenter(ix));
 626     if(q<20. || q>250.) continue; // ?!
 627
 628     h1 = h2->ProjectionY("py", ix, ix);
 629     entries = h1->GetEntries();
 630     if(entries < 50) continue;
 631     Adjust(&f, h1);
 632     h1->Fit(&f, "Q");
 633
 634     // Fill sy^2 = f(q)
 635     Int_t ip = gqm->GetN();
 636     gqm->SetPoint(ip, q, 1.e4*f.GetParameter(1));
 637     gqm->SetPointError(ip, 0., 1.e4*f.GetParError(1));
 638
 639     // correct sigma for ExB effect
 640     gqs->SetPoint(ip, q, 1.e4*f.GetParameter(2)/**f.GetParameter(2)-exb2*sxd2*/);
 641     gqs->SetPointError(ip, 0., 1.e4*f.GetParError(2)/**f.GetParameter(2)*/);
 642
 643     // save probability
 644     n += entries;
 645     gqp->SetPoint(ip, q, entries);
 646     gqp->SetPointError(ip, 0., 0./*TMath::Sqrt(entries)*/);
 647   }
 648
 649   // normalize probability and get mean sy
 650   Double_t sm = 0., sy;
 651   for(Int_t ip=gqp->GetN(); ip--;){
 652     gqp->GetPoint(ip, q, entries);
 653     entries/=n;
 654     gqp->SetPoint(ip, q, 1.e3*entries);
 655     gqs->GetPoint(ip, q, sy);
 656     sm += entries*sy;
 657   }
 658
 659   // error parametrization s(q) = <sy> + b(1/q-1/q0)
 660   TF1 fq("fq", "[0] + [1]/x", 20., 250.);
 661   gqs->Fit(&fq);
 662   //printf("sm=%f [0]=%f [1]=%f\n", 1.e-4*sm, fq.GetParameter(0), fq.GetParameter(1));
 663   printf("  const Float_t sq0inv = %f; // [1/q0]\n", (sm-fq.GetParameter(0))/fq.GetParameter(1));
 664   printf("  const Float_t sqb    = %f; // [cm]\n", 1.e-4*fq.GetParameter(1));
 665 }
 666
 667 //_______________________________________________________
 668 void AliTRDclusterResolution::ProcessCenterPad()
 669 {
 670 // Resolution as a function of y displacement from pad center
 671 // only for phi equal exB and clusters close to cathode wire plane
 672 // for ideal simulations time bins 4,5 and 6.
 673
 674   TObjArray *arr = (TObjArray*)fContainer->At(kCenter);
 675   if(!arr) {
 676     AliWarning("Missing dy=f(y | x, ly) container");
 677     return;
 678   }
 679   Float_t s[AliTRDgeometry::kNlayer];
 680   TF1 f("f", "gaus", -.5, .5);
 681   TF1 fp("fp", "gaus", -3.5, 3.5);
 682
 683   TH1D *h1 = 0x0; TH2F *h2 = 0x0; TH3S *h3r=0x0, *h3p=0x0;
 684   TObjArray *arrRes = (TObjArray*)fResults->At(kCenter);
 685   TTree *t = (TTree*)arrRes->At(0);
 686   TGraphErrors *gs = 0x0;
 687   TAxis *ax = 0x0;
 688
 689   printf("  const Float_t lSy[6][24] = {\n      {");
 690   const Int_t nl = AliTRDgeometry::kNlayer;
 691   for(Int_t il=0; il<nl; il++){
 692     if(!(h3r = (TH3S*)arr->At(il))) continue;
 693     if(!(h3p = (TH3S*)arr->At(nl+il))) continue;
 694     gs = (TGraphErrors*)arrRes->At(il+1);
 695     fLy = il;
 696     for(Int_t ix=1; ix<=h3r->GetXaxis()->GetNbins(); ix++){
 697       ax = h3r->GetXaxis(); ax->SetRange(ix, ix);
 698       ax = h3p->GetXaxis(); ax->SetRange(ix, ix);
 699       fX  = ax->GetBinCenter(ix);
 700       //printf("  x[%2d]=%4.2f\n", ix, fX);
 701       for(Int_t iy=1; iy<=h3r->GetYaxis()->GetNbins(); iy++){
 702         ax = h3r->GetYaxis(); ax->SetRange(iy, iy);
 703         ax = h3p->GetYaxis(); ax->SetRange(iy, iy);
 704         fY  = ax->GetBinCenter(iy);
 705         //printf("    y[%2d]=%5.2f\n", iy, fY);
 706         // finish navigation in the HnSparse
 707
 708         h1 = (TH1D*)h3r->Project3D("z");
 709         Int_t entries = (Int_t)h1->Integral();
 710         if(entries < 50) continue;
 711         //Adjust(&f, h1);
 712         h1->Fit(&f, "QN");
 713
 714         // Fill sy,my=f(y_w,x,ly)
 715         fR[0] = f.GetParameter(1); fR[1] = f.GetParError(1);
 716         fR[2] = f.GetParameter(2); fR[3] = f.GetParError(2);
 717
 718         h1 = (TH1D*)h3p->Project3D("z");
 719         h1->Fit(&fp, "QN");
 720         fP[0] = fp.GetParameter(1); fP[1] = fp.GetParError(1);
 721         fP[2] = fp.GetParameter(2); fP[3] = fp.GetParError(2);
 722
 723         //printf("ly[%d] x[%3.1f] y[%+5.2f] m[%5.3f] s[%5.3f] \n", fLy, fX, fY, fR[0], fR[2]);
 724         t->Fill();
 725
 726
 727       }
 728     }
 729     t->Draw("y:x>>h(24, 0., 2.4, 51, -.51, .51)",
 730             Form("s[0]*(ly==%d&&abs(m[0])<1.e-1)", fLy),
 731             "goff");
 732     h2=(TH2F*)gROOT->FindObject("h");
 733     f.FixParameter(1, 0.);
 734     Int_t n = h2->GetXaxis()->GetNbins(); s[il]=0.;
 735     printf("    {");
 736     for(Int_t ix=1; ix<=n; ix++){
 737       ax = h2->GetXaxis();
 738       fX  = ax->GetBinCenter(ix);
 739       h1 = h2->ProjectionY("hCenPy", ix, ix);
 740       //if((Int_t)h1->Integral() < 1.e-10) continue;
 741       h1->Fit(&f, "QN");
 742       s[il]+=f.GetParameter(2);
 743       printf("%6.4f,%s", f.GetParameter(0), ix%6?" ":"\n     ");
 744       Int_t jx = gs->GetN();
 745       gs->SetPoint(jx, fX, 1.e4*f.GetParameter(0));
 746       gs->SetPointError(jx, 0., 0./*f.GetParError(0)*/);
 747     }
 748     printf("\b},\n");
 749     s[il]/=n;
 750
 751     f.ReleaseParameter(2);
 752
 753
 754     if(!fCanvas) continue;
 755     h2->Draw("lego2fb");
 756     fCanvas->Modified(); fCanvas->Update();
 757     if(IsSaveAs()) fCanvas->SaveAs(Form("Figures/ProcessCenter_ly[%d].gif", fLy));
 758     else gSystem->Sleep(100);
 759   }
 760   printf("  };\n");
 761   printf("  const Float_t lPRF[] = {"
 762     "%5.3f, %5.3f, %5.3f, %5.3f, %5.3f, %5.3f};\n",
 763     s[0], s[1], s[2], s[3], s[4], s[5]);
 764 }
 765
 766 //_______________________________________________________
 767 void AliTRDclusterResolution::ProcessSigma()
 768 {
 769   TObjArray *arr = (TObjArray*)fContainer->At(kSigm);
 770   if(!arr){
 771     AliWarning("Missing dy=f(x_d, d_w) container");
 772     return;
 773   }
 774
 775   // init visualization
 776   TGraphErrors *ggs = 0x0;
 777   TGraph *line = 0x0;
 778   if(fCanvas){
 779     ggs = new TGraphErrors();
 780     line = new TGraph();
 781     line->SetLineColor(kRed);line->SetLineWidth(2);
 782   }
 783
 784   // init logistic support
 785   TF1 f("f", "gaus", -.5, .5);
 786   TLinearFitter gs(1,"pol1");
 787   TH1 *hFrame=0x0;
 788   TH1D *h1 = 0x0; TH3S *h3=0x0;
 789   TAxis *ax = 0x0;
 790   Double_t exb2 = fExB*fExB;
 791
 792   TTree *t = (TTree*)fResults->At(kSigm);
 793   for(Int_t ix=0; ix<kNTB; ix++){
 794     if(!(h3=(TH3S*)arr->At(ix))) continue;
 795     fX = fAt->GetBinCenter(ix+1);
 796
 797     for(Int_t iz=1; iz<=h3->GetXaxis()->GetNbins(); iz++){
 798       ax = h3->GetXaxis();
 799       ax->SetRange(iz, iz);
 800       fZ = ax->GetBinCenter(iz);
 801
 802       // reset visualization
 803       if(fCanvas){
 804         new(ggs) TGraphErrors();
 805         ggs->SetMarkerStyle(7);
 806       }
 807       gs.ClearPoints();
 808
 809       for(Int_t ip=1; ip<=h3->GetYaxis()->GetNbins(); ip++){
 810         ax = h3->GetYaxis();
 811         ax->SetRange(ip, ip);
 812         Double_t tgl = ax->GetBinCenter(ip);
 813         // finish navigation in the HnSparse
 814
 815         //if(TMath::Abs(dydx)>0.18) continue;
 816         Double_t tgg = (tgl-fExB)/(1.+tgl*fExB);
 817         Double_t tgg2 = tgg*tgg;
 818
 819         h1 = (TH1D*)h3->Project3D("z");
 820         Int_t entries = (Int_t)h1->Integral();
 821         if(entries < 50) continue;
 822         //Adjust(&f, h1);
 823         h1->Fit(&f, "QN");
 824
 825         Double_t s2  = f.GetParameter(2)*f.GetParameter(2);
 826         Double_t s2e = 2.*f.GetParameter(2)*f.GetParError(2);
 827         // Fill sy^2 = f(tg^2(phi-a_L))
 828         gs.AddPoint(&tgg2, s2, s2e);
 829
 830         if(!ggs) continue;
 831         Int_t jp = ggs->GetN();
 832         ggs->SetPoint(jp, tgg2, s2);
 833         ggs->SetPointError(jp, 0., s2e);
 834       }
 835       if(gs.Eval()) continue;
 836
 837       // s^2_x = s0^2_x - x^2*tg^2(a_L)/12
 838       fR[0] = gs.GetParameter(1)/* - x*x*exb2/12.*/;
 839       if(fR[0]<0.) continue;
 840       fR[0] = TMath::Sqrt(fR[0]);
 841       fR[1] = .5*gs.GetParError(1)/fR[0];
 842
 843       // s^2_y  = s0^2_y + tg^2(a_L) * s^2_x
 844       // s0^2_y = f(D_L)*x + s_PRF^2
 845       fR[2]= gs.GetParameter(0)/*-exb2*sx*/;
 846       if(fR[1] <0.) continue;
 847       fR[2] = TMath::Sqrt(fR[2]);
 848       fR[3] = gs.GetParError(0)+exb2*exb2*gs.GetParError(1);
 849       t->Fill();
 850
 851       if(!fCanvas) continue;
 852       fCanvas->cd(); fCanvas->SetLogx(); //fCanvas->SetLogy();
 853       if(!hFrame){
 854         hFrame=new TH1I("hFrame", "", 100, 0., .3);
 855         hFrame->SetMinimum(0.);hFrame->SetMaximum(.005);
 856         hFrame->SetXTitle("tg^{2}(#phi-#alpha_{L})");
 857         hFrame->SetYTitle("#sigma^{2}y[cm^{2}]");
 858         hFrame->SetLineColor(1);hFrame->SetLineWidth(1);
 859         hFrame->Draw();
 860       } else hFrame->Reset();
 861       Double_t xx = 0., dxx=.2/50;
 862       for(Int_t ip=0;ip<50;ip++){
 863         line->SetPoint(ip, xx,  gs.GetParameter(0)+xx*gs.GetParameter(1));
 864         xx+=dxx;
 865       }
 866       ggs->Draw("pl"); line->Draw("l");
 867       fCanvas->Modified(); fCanvas->Update();
 868       if(IsSaveAs()) fCanvas->SaveAs(Form("Figures/ProcessSigma_z[%5.3f]_x[%5.3f].gif", fZ, fX));
 869       else gSystem->Sleep(100);
 870
 871       printf("    xd=%4.1f[cm] sx=%5.3e[cm] sy=%5.3e[cm]\n", fX, TMath::Sqrt(fR[0]), TMath::Sqrt(fR[1]));
 872     }
 873   }
 874   return;
 875 }
 876
 877 //_______________________________________________________
 878 void AliTRDclusterResolution::ProcessMean()
 879 {
 880   TObjArray *arr = (TObjArray*)fContainer->At(kMean);
 881   if(!arr){
 882     AliWarning("Missing dy=f(x_d, d_w) container");
 883     return;
 884   }
 885
 886   // init logistic support
 887   TF1 f("f", "gaus", -.5, .5);
 888   TF1 line("l", "[0]+[1]*x", -.15, .15);
 889   TGraphErrors *gm = new TGraphErrors();
 890   TH1 *hFrame=0x0;
 891   TH1D *h1 = 0x0; TH3S *h3 =0x0;
 892   TAxis *ax = 0x0;
 893
 894   TTree *t = (TTree*)fResults->At(kMean);
 895   for(Int_t ix=0; ix<kNTB; ix++){
 896     if(!(h3=(TH3S*)arr->At(ix))) continue;
 897     fX = fAt->GetBinCenter(ix+1);
 898
 899     for(Int_t iz=1; iz<=h3->GetXaxis()->GetNbins(); iz++){
 900       ax = h3->GetXaxis();
 901       ax->SetRange(iz, iz);
 902       fZ = ax->GetBinCenter(iz);
 903
 904       // reset fitter
 905       new(gm) TGraphErrors();
 906       gm->SetMarkerStyle(7);
 907
 908       for(Int_t ip=1; ip<=h3->GetYaxis()->GetNbins(); ip++){
 909         ax = h3->GetYaxis();
 910         ax->SetRange(ip, ip);
 911         Double_t tgl = ax->GetBinCenter(ip);
 912         // finish navigation in the HnSparse
 913
 914         h1 = (TH1D*)h3->Project3D("z");
 915         Int_t entries = (Int_t)h1->Integral();
 916         if(entries < 50) continue;
 917         //Adjust(&f, h1);
 918         h1->Fit(&f, "QN");
 919
 920         // Fill <Dy> = f(dydx - h*dzdx)
 921         Int_t jp = gm->GetN();
 922         gm->SetPoint(jp, tgl, f.GetParameter(1));
 923         gm->SetPointError(jp, 0., f.GetParError(1));
 924       }
 925       if(gm->GetN()<4) continue;
 926
 927       gm->Fit(&line, "QN");
 928       fR[0] = line.GetParameter(1); // dx
 929       fR[1] = line.GetParError(1);
 930       fR[2] = line.GetParameter(0) + fExB*fR[0]; // xs = dy - tg(a_L)*dx
 931       t->Fill();
 932
 933       if(!fCanvas) continue;
 934       fCanvas->cd();
 935       if(!hFrame){
 936         hFrame=new TH1I("hFrame", "", 100, -.3, .3);
 937         hFrame->SetMinimum(-.1);hFrame->SetMaximum(.1);
 938         hFrame->SetXTitle("tg#phi-htg#theta");
 939         hFrame->SetYTitle("#Deltay[cm]");
 940         hFrame->SetLineColor(1);hFrame->SetLineWidth(1);
 941         hFrame->Draw();
 942       } else hFrame->Reset();
 943       gm->Draw("pl"); line.Draw("same");
 944       fCanvas->Modified(); fCanvas->Update();
 945       if(IsSaveAs()) fCanvas->SaveAs(Form("Figures/ProcessMean_Z[%5.3f]_X[%5.3f].gif", fZ, fX));
 946       else gSystem->Sleep(100);
 947       printf("    xd=%4.2f[cm] dx=%5.3e[cm] dy=%5.3e[cm]\n", fX, fR[0], fR[2]);
 948     }
 949   }
 950
 951   // draw shift results
 952   //t->Draw("z:x>>h(24, 0, 2.4, 25, 0, 2.5)", "dx*(abs(dx)<1.e-2)", "lego2fb");
 953   //t->Draw("z:x>>h(24, 0, 2.4, 25, 0, 2.5)", "dy*(abs(dx)<1.e-2)", "lego2fb");
 954
 955 }