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