PWG1/TRD/AliTRDresolution.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: AliTRDresolution.cxx 27496 2008-07-22 08:35:45Z cblume $ */
  17
  18 ////////////////////////////////////////////////////////////////////////////
  19 //                                                                        //
  20 //  TRD tracking resolution                                               //
  21 //
  22 // The class performs resolution and residual studies
  23 // of the TRD tracks for the following quantities :
  24 //   - spatial position (y, [z])
  25 //   - angular (phi) tracklet
  26 //   - momentum at the track level
  27 //
  28 // The class has to be used for regular detector performance checks using the official macros:
  29 //   - $ALICE_ROOT/TRD/qaRec/run.C
  30 //   - $ALICE_ROOT/TRD/qaRec/makeResults.C
  31 //
  32 // For stand alone usage please refer to the following example:
  33 // {
  34 //   gSystem->Load("libANALYSIS.so");
  35 //   gSystem->Load("libTRDqaRec.so");
  36 //   AliTRDresolution *res = new AliTRDresolution();
  37 //   //res->SetMCdata();
  38 //   //res->SetVerbose();
  39 //   //res->SetVisual();
  40 //   res->Load();
  41 //   if(!res->PostProcess()) return;
  42 //   res->GetRefFigure(0);
  43 // }
  44 //
  45 //  Authors:                                                              //
  46 //    Alexandru Bercuci <A.Bercuci@gsi.de>                                //
  47 //    Markus Fasel <M.Fasel@gsi.de>                                       //
  48 //                                                                        //
  49 ////////////////////////////////////////////////////////////////////////////
  50
  51 #include <TSystem.h>
  52 #include <TStyle.h>
  53 #include <TROOT.h>
  54 #include <TObjArray.h>
  55 #include <TH3.h>
  56 #include <TH2.h>
  57 #include <TH1.h>
  58 #include <THnSparse.h>
  59 #include <TF1.h>
  60 #include <TCanvas.h>
  61 #include <TGaxis.h>
  62 #include <TBox.h>
  63 #include <TLegend.h>
  64 #include <TGraphErrors.h>
  65 #include <TGraphAsymmErrors.h>
  66 #include <TLinearFitter.h>
  67 #include <TMath.h>
  68 #include <TMatrixT.h>
  69 #include <TVectorT.h>
  70 #include <TTreeStream.h>
  71 #include <TGeoManager.h>
  72 #include <TDatabasePDG.h>
  73
  74 #include "AliPID.h"
  75 #include "AliLog.h"
  76 #include "AliESDtrack.h"
  77 #include "AliMathBase.h"
  78 #include "AliTrackPointArray.h"
  79
  80 #include "AliTRDresolution.h"
  81 #include "AliTRDgeometry.h"
  82 #include "AliTRDtransform.h"
  83 #include "AliTRDpadPlane.h"
  84 #include "AliTRDcluster.h"
  85 #include "AliTRDseedV1.h"
  86 #include "AliTRDtrackV1.h"
  87 #include "AliTRDReconstructor.h"
  88 #include "AliTRDrecoParam.h"
  89 #include "AliTRDpidUtil.h"
  90 #include "AliTRDinfoGen.h"
  91
  92 #include "info/AliTRDclusterInfo.h"
  93
  94 ClassImp(AliTRDresolution)
  95 //ClassImp(AliTRDresolution::AliTRDresolutionProjection)
  96
  97 Int_t const   AliTRDresolution::fgkNbins[kNdim] = {
  98   Int_t(kNbunchCross)/*bc*/,
  99   180/*phi*/,
 100   50/*eta*/,
 101   50/*dy*/,
 102   40/*dphi*/,
 103   50/*dz*/,
 104   Int_t(kNcharge)*AliPID::kSPECIES+1/*chg*species*/,
 105   kNpt/*pt*/
 106 };  //! no of bins/projection
 107 Double_t const AliTRDresolution::fgkMin[kNdim] = {
 108   -0.5,
 109   -TMath::Pi(),
 110   -1.,
 111   -1.5,
 112   -10.,
 113   -2.5,
 114   -AliPID::kSPECIES-0.5,
 115   -0.5
 116 };    //! low limits for projections
 117 Double_t const AliTRDresolution::fgkMax[kNdim] = {
 118   Int_t(kNbunchCross)-0.5,
 119   TMath::Pi(),
 120   1.,
 121   1.5,
 122   10.,
 123   2.5,
 124   AliPID::kSPECIES+0.5,
 125   kNpt-0.5
 126 };    //! high limits for projections
 127 Char_t const *AliTRDresolution::fgkTitle[kNdim] = {
 128   "bunch cross",
 129   "#phi [rad]",
 130   "#eta",
 131   "#Deltay [cm]",
 132   "#Delta#phi [deg]",
 133   "#Deltaz [cm]",
 134   "chg*spec*rc",
 135   "bin_p_{t}"
 136 };  //! title of projection
 137
 138 Char_t const * AliTRDresolution::fgPerformanceName[kNclasses] = {
 139     "Cluster2Track"
 140     ,"Tracklet2Track"
 141     ,"Tracklet2TRDin"
 142     ,"Cluster2MC"
 143     ,"Tracklet2MC"
 144     ,"TRDin2MC"
 145     ,"TRD2MC"
 146 //    ,"Tracklet2TRDout"
 147 //    ,"TRDout2MC"
 148 };
 149 Float_t AliTRDresolution::fgPtBin[kNpt+1];
 150
 151 //________________________________________________________
 152 AliTRDresolution::AliTRDresolution()
 153   :AliTRDrecoTask()
 154   ,fIdxPlot(0)
 155   ,fIdxFrame(0)
 156   ,fPtThreshold(1.)
 157   ,fDyRange(0.75)
 158   ,fProj(NULL)
 159   ,fDBPDG(NULL)
 160   ,fCl(NULL)
 161   ,fMCcl(NULL)
 162 {
 163   //
 164   // Default constructor
 165   //
 166   SetNameTitle("TRDresolution", "TRD spatial and momentum resolution");
 167   MakePtSegmentation();
 168 }
 169
 170 //________________________________________________________
 171 AliTRDresolution::AliTRDresolution(char* name, Bool_t xchange)
 172   :AliTRDrecoTask(name, "TRD spatial and momentum resolution")
 173   ,fIdxPlot(0)
 174   ,fIdxFrame(0)
 175   ,fPtThreshold(1.)
 176   ,fDyRange(0.75)
 177   ,fProj(NULL)
 178   ,fDBPDG(NULL)
 179   ,fCl(NULL)
 180   ,fMCcl(NULL)
 181 {
 182   //
 183   // Default constructor
 184   //
 185
 186   InitFunctorList();
 187   MakePtSegmentation();
 188   if(xchange){
 189     SetUseExchangeContainers();
 190     DefineOutput(kClToTrk, TObjArray::Class()); // cluster2track
 191     DefineOutput(kClToMC, TObjArray::Class()); // cluster2mc
 192   }
 193 }
 194
 195 //________________________________________________________
 196 AliTRDresolution::~AliTRDresolution()
 197 {
 198   //
 199   // Destructor
 200   //
 201
 202   if(fProj){fProj->Delete(); delete fProj;}
 203   if(fCl){fCl->Delete(); delete fCl;}
 204   if(fMCcl){fMCcl->Delete(); delete fMCcl;}
 205 }
 206
 207
 208 //________________________________________________________
 209 void AliTRDresolution::UserCreateOutputObjects()
 210 {
 211   // spatial resolution
 212
 213   AliTRDrecoTask::UserCreateOutputObjects();
 214   if(UseExchangeContainers()) InitExchangeContainers();
 215 }
 216
 217 //________________________________________________________
 218 void AliTRDresolution::InitExchangeContainers()
 219 {
 220 // Init containers for subsequent tasks (AliTRDclusterResolution)
 221
 222   fCl = new TObjArray(200); fCl->SetOwner(kTRUE);
 223   fMCcl = new TObjArray(); fMCcl->SetOwner(kTRUE);
 224   PostData(kClToTrk, fCl);
 225   PostData(kClToMC, fMCcl);
 226 }
 227
 228 //________________________________________________________
 229 void AliTRDresolution::UserExec(Option_t *opt)
 230 {
 231   //
 232   // Execution part
 233   //
 234
 235   if(fCl) fCl->Delete();
 236   if(fMCcl) fMCcl->Delete();
 237   AliTRDrecoTask::UserExec(opt);
 238 }
 239
 240 //________________________________________________________
 241 Bool_t AliTRDresolution::Pulls(Double_t* /*dyz[2]*/, Double_t* /*cov[3]*/, Double_t /*tilt*/) const
 242 {
 243 // Helper function to calculate pulls in the yz plane
 244 // using proper tilt rotation
 245 // Uses functionality defined by AliTRDseedV1.
 246
 247   return kTRUE;
 248 /*
 249   Double_t t2(tilt*tilt);
 250   // exit door until a bug fix is found for AliTRDseedV1::GetCovSqrt
 251
 252   // rotate along pad
 253   Double_t cc[3];
 254   cc[0] = cov[0] - 2.*tilt*cov[1] + t2*cov[2];
 255   cc[1] = cov[1]*(1.-t2) + tilt*(cov[0] - cov[2]);
 256   cc[2] = t2*cov[0] + 2.*tilt*cov[1] + cov[2];
 257   // do sqrt
 258   Double_t sqr[3]={0., 0., 0.};
 259   if(AliTRDseedV1::GetCovSqrt(cc, sqr)) return kFALSE;
 260   Double_t invsqr[3]={0., 0., 0.};
 261   if(AliTRDseedV1::GetCovInv(sqr, invsqr)<1.e-40) return kFALSE;
 262   Double_t tmp(dyz[0]);
 263   dyz[0] = invsqr[0]*tmp + invsqr[1]*dyz[1];
 264   dyz[1] = invsqr[1]*tmp + invsqr[2]*dyz[1];
 265   return kTRUE;
 266 */
 267 }
 268
 269 //________________________________________________________
 270 TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track)
 271 {
 272   //
 273   // Plot the cluster distributions
 274   //
 275
 276   if(track) fkTrack = track;
 277   if(!fkTrack){
 278     AliDebug(4, "No Track defined.");
 279     return NULL;
 280   }
 281   if(TMath::Abs(fkESD->GetTOFbc()) > 1){
 282     AliDebug(4, Form("Track with BC_index[%d] not used.", fkESD->GetTOFbc()));
 283     return NULL;
 284   }
 285   if(fPt<fPtThreshold){
 286     AliDebug(4, Form("Track with pt[%6.4f] under threshold.", fPt));
 287     return NULL;
 288   }
 289   THnSparse *H(NULL);
 290   if(!fContainer || !(H = ((THnSparse*)fContainer->At(kCluster)))){
 291     AliWarning("No output container defined.");
 292     return NULL;
 293   }
 294
 295   AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
 296   Double_t val[kNdim]; //Float_t exb, vd, t0, s2, dl, dt;
 297   TObjArray     *clInfoArr(NULL);
 298   AliTRDseedV1  *fTracklet(NULL);
 299   AliTRDcluster *c(NULL), *cc(NULL);
 300   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 301     if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
 302     if(!fTracklet->IsOK()) continue;
 303     //fTracklet->GetCalibParam(exb, vd, t0, s2, dl, dt);
 304     val[kBC]  = ily;
 305     val[kPhi] = fPhi;
 306     val[kEta] = fEta;
 307     val[kPt]  = TMath::ATan(fTracklet->GetYref(1))*TMath::RadToDeg();
 308     Float_t corr = 1./TMath::Sqrt(1.+fTracklet->GetYref(1)*fTracklet->GetYref(1)+fTracklet->GetZref(1)*fTracklet->GetZref(1));
 309     Int_t row0(-1);
 310     Float_t padCorr(fTracklet->GetTilt()*fTracklet->GetPadLength());
 311     fTracklet->ResetClusterIter(kTRUE);
 312     while((c = fTracklet->NextCluster())){
 313       Float_t xc(c->GetX()),
 314               q(TMath::Abs(c->GetQ()));
 315       if(row0<0) row0 = c->GetPadRow();
 316
 317       val[kYrez] = c->GetY() + padCorr*(c->GetPadRow() - row0) -fTracklet->GetYat(xc);
 318       val[kPrez] = fTracklet->GetX0()-xc;
 319       val[kZrez] = 0.; Int_t ic(0), tb(c->GetLocalTimeBin());;
 320       if((cc = fTracklet->GetClusters(tb-1))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
 321       if((cc = fTracklet->GetClusters(tb-2))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
 322       if(ic) val[kZrez] /= (ic*q);
 323       val[kSpeciesChgRC]= fTracklet->IsRowCross()?0.:(TMath::Max(q*corr, Float_t(3.)));
 324       H->Fill(val);
 325 /*      // tilt rotation of covariance for clusters
 326       Double_t sy2(c->GetSigmaY2()), sz2(c->GetSigmaZ2());
 327       cov[0] = (sy2+t2*sz2)*corr;
 328       cov[1] = tilt*(sz2 - sy2)*corr;
 329       cov[2] = (t2*sy2+sz2)*corr;
 330       // sum with track covariance
 331       cov[0]+=covR[0]; cov[1]+=covR[1]; cov[2]+=covR[2];
 332       Double_t dyz[2]= {dy[1], dz[1]};
 333       Pulls(dyz, cov, tilt);*/
 334
 335       // Get z-position with respect to anode wire
 336       Float_t yt(fTracklet->GetYref(0)-val[kZrez]*fTracklet->GetYref(1)),
 337               zt(fTracklet->GetZref(0)-val[kZrez]*fTracklet->GetZref(1));
 338       Int_t istk = geo->GetStack(c->GetDetector());
 339       AliTRDpadPlane *pp = geo->GetPadPlane(ily, istk);
 340       Float_t rowZ = pp->GetRow0();
 341       Float_t d  = rowZ - zt + pp->GetAnodeWireOffset();
 342       d -= ((Int_t)(2 * d)) / 2.0;
 343       if (d > 0.25) d  = 0.5 - d;
 344
 345       AliTRDclusterInfo *clInfo(NULL);
 346       clInfo = new AliTRDclusterInfo;
 347       clInfo->SetCluster(c);
 348       //Float_t covF[] = {cov[0], cov[1], cov[2]};
 349       clInfo->SetGlobalPosition(yt, zt, fTracklet->GetYref(1), fTracklet->GetZref(1)/*, covF*/);
 350       clInfo->SetResolution(val[kYrez]);
 351       clInfo->SetAnisochronity(d);
 352       clInfo->SetDriftLength(val[kZrez]);
 353       clInfo->SetTilt(fTracklet->GetTilt());
 354       if(fCl) fCl->Add(clInfo);
 355       //else AliDebug(1, "Cl exchange container missing. Activate by calling \"InitExchangeContainers()\"");
 356
 357       if(DebugLevel()>=2){
 358         if(!clInfoArr){
 359           clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
 360           clInfoArr->SetOwner(kFALSE);
 361         }
 362         clInfoArr->Add(clInfo);
 363       }
 364     }
 365     if(DebugLevel()>=2 && clInfoArr){
 366       ULong_t status = fkESD->GetStatus();
 367       (*DebugStream()) << "cluster"
 368         <<"status="  << status
 369         <<"clInfo.=" << clInfoArr
 370         << "\n";
 371       clInfoArr->Clear();
 372     }
 373   }
 374   if(clInfoArr) delete clInfoArr;
 375
 376   return NULL;//H->Projection(kEta, kPhi);
 377 }
 378
 379
 380 //________________________________________________________
 381 TH1* AliTRDresolution::PlotTracklet(const AliTRDtrackV1 *track)
 382 {
 383 // Plot normalized residuals for tracklets to track.
 384 //
 385 // We start from the result that if X=N(|m|, |Cov|)
 386 // BEGIN_LATEX
 387 // (Cov^{-1})^{1/2}X = N((Cov^{-1})^{1/2}*|m|, |1|)
 388 // END_LATEX
 389 // in our case X=(y_trklt - y_trk z_trklt - z_trk) and |Cov| = |Cov_trklt| + |Cov_trk| at the radial
 390 // reference position.
 391   if(track) fkTrack = track;
 392   if(!fkTrack){
 393     AliDebug(4, "No Track defined.");
 394     return NULL;
 395   }
 396   if(TMath::Abs(fkESD->GetTOFbc())>1){
 397     AliDebug(4, Form("Track with BC_index[%d] not used.", fkESD->GetTOFbc()));
 398     return NULL;
 399   }
 400   THnSparse *H(NULL);
 401   if(!fContainer || !(H = (THnSparse*)fContainer->At(kTracklet))){
 402     AliWarning("No output container defined.");
 403     return NULL;
 404   }
 405 //  return NULL;
 406   Double_t val[kNdim+1];
 407   AliTRDseedV1 *fTracklet(NULL);
 408   for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++){
 409     if(!(fTracklet = fkTrack->GetTracklet(il))) continue;
 410     if(!fTracklet->IsOK()) continue;
 411     val [kBC] = il;
 412     val[kPhi] = fPhi;
 413     val[kEta] = fEta;
 414     val[kSpeciesChgRC]= fSpecies;
 415     val[kPt]  = GetPtBin(fTracklet->GetMomentum());
 416     Double_t dyt(fTracklet->GetYref(0) - fTracklet->GetYfit(0)),
 417              dzt(fTracklet->GetZref(0) - fTracklet->GetZfit(0)),
 418              dydx(fTracklet->GetYfit(1)),
 419              tilt(fTracklet->GetTilt());
 420     // correct for tilt rotation
 421     val[kYrez] = dyt - dzt*tilt;
 422     val[kZrez] = dzt + dyt*tilt;
 423     dydx+= tilt*fTracklet->GetZref(1);
 424     val[kPrez] = TMath::ATan((fTracklet->GetYref(1) - dydx)/(1.+ fTracklet->GetYref(1)*dydx)) * TMath::RadToDeg();
 425     if(fTracklet->IsRowCross()){
 426       val[kSpeciesChgRC]= 0.;
 427 //      val[kPrez] = fkTrack->Charge(); // may be better defined
 428     }/* else {
 429       Float_t exb, vd, t0, s2, dl, dt;
 430       fTracklet->GetCalibParam(exb, vd, t0, s2, dl, dt);
 431       val[kZrez] = TMath::ATan((fTracklet->GetYref(1) - exb)/(1+fTracklet->GetYref(1)*exb));
 432     }*/
 433     val[kNdim] = fTracklet->GetdQdl();
 434     if(DebugLevel()>=1) H->Fill(val);
 435
 436 //     // compute covariance matrix
 437 //     fTracklet->GetCovAt(x, cov);
 438 //     fTracklet->GetCovRef(covR);
 439 //     cov[0] += covR[0]; cov[1] += covR[1]; cov[2] += covR[2];
 440 //     Double_t dyz[2]= {dy[1], dz[1]};
 441 //     Pulls(dyz, cov, tilt);
 442 //     ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dyz[0], dyz[1]);
 443 //     ((TH3S*)arr->At(3))->Fill(tht, dyz[1], rc);
 444
 445     if(DebugLevel()>=3){
 446       Bool_t rc(fTracklet->IsRowCross());
 447       UChar_t err(fTracklet->GetErrorMsg());
 448       Double_t x(fTracklet->GetX()),
 449                pt(fTracklet->GetPt()),
 450                yt(fTracklet->GetYref(0)),
 451                zt(fTracklet->GetZref(0)),
 452                phi(fTracklet->GetYref(1)),
 453                tht(fTracklet->GetZref(1));
 454       Int_t ncl(fTracklet->GetN()),
 455             det(fTracklet->GetDetector());
 456       (*DebugStream()) << "tracklet"
 457         <<"pt="  << pt
 458         <<"x="   << x
 459         <<"yt="  << yt
 460         <<"zt="  << zt
 461         <<"phi=" << phi
 462         <<"tht=" << tht
 463         <<"det=" << det
 464         <<"n="   << ncl
 465         <<"dy0=" << dyt
 466         <<"dz0=" << dzt
 467         <<"dy="  << val[kYrez]
 468         <<"dz="  << val[kZrez]
 469         <<"dphi="<< val[kPrez]
 470         <<"dQ  ="<< val[kNdim]
 471         <<"rc="  << rc
 472         <<"err=" << err
 473         << "\n";
 474     }
 475   }
 476   return NULL;//H->Projection(kEta, kPhi);
 477 }
 478
 479
 480 //________________________________________________________
 481 TH1* AliTRDresolution::PlotTrackIn(const AliTRDtrackV1 *track)
 482 {
 483 // Store resolution/pulls of Kalman before updating with the TRD information
 484 // at the radial position of the first tracklet. The following points are used
 485 // for comparison
 486 //  - the (y,z,snp) of the first TRD tracklet
 487 //  - the (y, z, snp, tgl, pt) of the MC track reference
 488 //
 489 // Additionally the momentum resolution/pulls are calculated for usage in the
 490 // PID calculation.
 491   //printf("AliTRDresolution::PlotTrackIn() :: track[%p]\n", (void*)track);
 492
 493   if(track) fkTrack = track;
 494   if(!fkTrack){
 495     AliDebug(4, "No Track defined.");
 496     return NULL;
 497   }
 498   //fkTrack->Print();
 499   // check container
 500   THnSparseI *H=(THnSparseI*)fContainer->At(kTrackIn);
 501   if(!H){
 502     AliError(Form("Missing container @ %d", Int_t(kTrackIn)));
 503     return NULL;
 504   }
 505   // check input track status
 506   AliExternalTrackParam *tin(NULL);
 507   if(!(tin = fkTrack->GetTrackIn())){
 508     AliError("Track did not entered TRD fiducial volume.");
 509     return NULL;
 510   }
 511   // check first tracklet
 512   AliTRDseedV1 *fTracklet(fkTrack->GetTracklet(0));
 513   if(!fTracklet){
 514     AliDebug(3, "No Tracklet in ly[0]. Skip track.");
 515     return NULL;
 516   }
 517   // check radial position
 518   Double_t x = tin->GetX();
 519   if(TMath::Abs(x-fTracklet->GetX())>1.e-3){
 520     AliDebug(1, Form("Tracklet did not match Track. dx[cm]=%+4.1f", x-fTracklet->GetX()));
 521     return NULL;
 522   }
 523   //printf("USE y[%+f] dydx[%+f]\n", fTracklet->GetYfit(0), fTracklet->GetYfit(1));
 524
 525   Int_t bc(TMath::Abs(fkESD->GetTOFbc())%2);
 526   const Double_t *parR(tin->GetParameter());
 527   Double_t dyt(parR[0] - fTracklet->GetYfit(0)), dzt(parR[1] - fTracklet->GetZfit(0)),
 528             phit(fTracklet->GetYfit(1)),
 529             tilt(fTracklet->GetTilt());
 530
 531   // correct for tilt rotation
 532   Double_t dy  = dyt - dzt*tilt,
 533            dz  = dzt + dyt*tilt;
 534   phit       += tilt*parR[3];
 535   Double_t dphi = TMath::ASin(parR[2])-TMath::ATan(phit);
 536
 537   Double_t val[kNdim];
 538   val[kBC]          = bc;
 539   val[kPhi]         = fPhi;
 540   val[kEta]         = fEta;
 541   val[kSpeciesChgRC]= fTracklet->IsRowCross()?0:fSpecies;
 542   val[kPt]          = GetPtBin(fPt);
 543   val[kYrez]        = dy;
 544   val[kZrez]        = dz;
 545   val[kPrez]        = dphi*TMath::RadToDeg();
 546   H->Fill(val);
 547   if(DebugLevel()>=3){
 548     (*DebugStream()) << "trackIn"
 549       <<"tracklet.="  << fTracklet
 550       <<"trackIn.="   << tin
 551       << "\n";
 552   }
 553
 554   return NULL; // H->Projection(kEta, kPhi);
 555 }
 556
 557 /*
 558 //________________________________________________________
 559 TH1* AliTRDresolution::PlotTrackOut(const AliTRDtrackV1 *track)
 560 {
 561 // Store resolution/pulls of Kalman after last update with the TRD information
 562 // at the radial position of the first tracklet. The following points are used
 563 // for comparison
 564 //  - the (y,z,snp) of the first TRD tracklet
 565 //  - the (y, z, snp, tgl, pt) of the MC track reference
 566 //
 567 // Additionally the momentum resolution/pulls are calculated for usage in the
 568 // PID calculation.
 569
 570   if(track) fkTrack = track;
 571   return NULL;
 572 }
 573 */
 574 //________________________________________________________
 575 TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
 576 {
 577   //
 578   // Plot MC distributions
 579   //
 580
 581   if(!HasMCdata()){
 582     AliDebug(2, "No MC defined. Results will not be available.");
 583     return NULL;
 584   }
 585   if(track) fkTrack = track;
 586   if(!fkTrack){
 587     AliDebug(4, "No Track defined.");
 588     return NULL;
 589   }
 590   Int_t bc(TMath::Abs(fkESD->GetTOFbc()));
 591
 592   THnSparse *H(NULL);
 593   if(!fContainer){
 594     AliWarning("No output container defined.");
 595     return NULL;
 596   }
 597   // retriev track characteristics
 598   Int_t pdg = fkMC->GetPDG(),
 599         sIdx(AliTRDpidUtil::Pdg2Pid(TMath::Abs(pdg))+1), // species index
 600         sign(0),
 601 //        sgm[3],
 602         label(fkMC->GetLabel());
 603 //        fSegmentLevel(0);
 604   if(!fDBPDG) fDBPDG=TDatabasePDG::Instance();
 605   TParticlePDG *ppdg(fDBPDG->GetParticle(pdg));
 606   if(ppdg) sign = ppdg->Charge() > 0. ? 1 : -1;
 607
 608   TH1 *h(NULL);
 609   AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
 610   AliTRDseedV1 *fTracklet(NULL); TObjArray *clInfoArr(NULL);
 611   UChar_t s;
 612   Double_t x, y, z, pt, dydx, dzdx, dzdl;
 613   Float_t pt0, x0, y0, z0, dx, dy, dz, dydx0, dzdx0;
 614   Double_t covR[7]/*, cov[3]*/;
 615
 616 /*  if(DebugLevel()>=3){
 617     // get first detector
 618     Int_t det = -1;
 619     for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 620       if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
 621       det = fTracklet->GetDetector();
 622       break;
 623     }
 624     if(det>=0){
 625       TVectorD X(12), Y(12), Z(12), dX(12), dY(12), dZ(12), vPt(12), dPt(12), budget(12), cCOV(12*15);
 626       Double_t m(-1.);
 627       m = fkTrack->GetMass();
 628       if(fkMC->PropagateKalman(&X, &Y, &Z, &dX, &dY, &dZ, &vPt, &dPt, &budget, &cCOV, m)){
 629         (*DebugStream()) << "MCkalman"
 630           << "pdg=" << pdg
 631           << "det=" << det
 632           << "x="   << &X
 633           << "y="   << &Y
 634           << "z="   << &Z
 635           << "dx="  << &dX
 636           << "dy="  << &dY
 637           << "dz="  << &dZ
 638           << "pt="  << &vPt
 639           << "dpt=" << &dPt
 640           << "bgt=" << &budget
 641           << "cov=" << &cCOV
 642           << "\n";
 643       }
 644     }
 645   }*/
 646   AliTRDcluster *c(NULL);
 647   Double_t val[kNdim+1];
 648   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 649     if(!(fTracklet = fkTrack->GetTracklet(ily)))/* ||
 650        !fTracklet->IsOK())*/ continue;
 651
 652     x= x0 = fTracklet->GetX();
 653     Bool_t rc(fTracklet->IsRowCross()); Float_t eta, phi;
 654     if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, eta, phi, s)) continue;
 655
 656     // MC track position at reference radial position
 657     dx  = x0 - x;
 658     Float_t ymc = y0 - dx*dydx0;
 659     Float_t zmc = z0 - dx*dzdx0;
 660     phi -= TMath::Pi();
 661
 662     val[kBC]  = ily;
 663     val[kPhi] = phi;
 664     val[kEta] = eta;
 665     val[kSpeciesChgRC]= rc?0.:sign*sIdx;
 666     val[kPt]  = GetPtBin(pt0);
 667     Double_t tilt(fTracklet->GetTilt());
 668 //             ,t2(tilt*tilt)
 669 //             ,corr(1./(1. + t2))
 670 //             ,cost(TMath::Sqrt(corr));
 671
 672     AliExternalTrackParam *tin(fkTrack->GetTrackIn());
 673     if(ily==0 && tin){ // trackIn residuals
 674       // check radial position
 675       if(TMath::Abs(tin->GetX()-x)>1.e-3) AliDebug(1, Form("TrackIn radial mismatch. dx[cm]=%+4.1f", tin->GetX()-x));
 676       else{
 677       //  Float_t phi = TMath::ATan2(y0, x0);
 678         val[kBC]          = (bc>=kNbunchCross)?(kNbunchCross-1):bc;
 679         val[kYrez]        = tin->GetY()-ymc;
 680         val[kZrez]        = tin->GetZ()-zmc;
 681         val[kPrez]        = (TMath::ASin(tin->GetSnp())-TMath::ATan(dydx0))*TMath::RadToDeg();
 682         if((H = (THnSparseI*)fContainer->At(kMCtrackIn))) H->Fill(val);
 683       }
 684     }
 685     if(bc>1) break; // do nothing for the rest of TRD objects if satellite bunch
 686
 687     // track residuals
 688     dydx = fTracklet->GetYref(1);
 689     dzdx = fTracklet->GetZref(1);
 690     dzdl = fTracklet->GetTgl();
 691     y  = fTracklet->GetYref(0);
 692     dy = y - ymc;
 693     z  = fTracklet->GetZref(0);
 694     dz = z - zmc;
 695     pt = TMath::Abs(fTracklet->GetPt());
 696     fTracklet->GetCovRef(covR);
 697
 698     val[kYrez] = dy;
 699     val[kPrez] = TMath::ATan((dydx - dydx0)/(1.+ dydx*dydx0))*TMath::RadToDeg();
 700     val[kZrez] = dz;
 701     val[kNdim] = pt/pt0-1.;
 702     if((H = (THnSparse*)fContainer->At(kMCtrack))) H->Fill(val);
 703 /*      // theta resolution/ tgl pulls
 704       Double_t dzdl0 = dzdx0/TMath::Sqrt(1.+dydx0*dydx0),
 705                 dtgl = (dzdl - dzdl0)/(1.- dzdl*dzdl0);
 706       ((TH2I*)arr->At(6))->Fill(dzdl0, TMath::ATan(dtgl));
 707       ((TH2I*)arr->At(7))->Fill(dzdl0, (dzdl - dzdl0)/TMath::Sqrt(covR[4]));
 708       // pt resolution  \\ 1/pt pulls \\ p resolution for PID
 709       Double_t p0 = TMath::Sqrt(1.+ dzdl0*dzdl0)*pt0,
 710               p  = TMath::Sqrt(1.+ dzdl*dzdl)*pt;
 711       ((TH3S*)((TObjArray*)arr->At(8)))->Fill(pt0, pt/pt0-1., sign*sIdx);
 712       ((TH3S*)((TObjArray*)arr->At(9)))->Fill(1./pt0, (1./pt-1./pt0)/TMath::Sqrt(covR[6]), sign*sIdx);
 713       ((TH3S*)((TObjArray*)arr->At(10)))->Fill(p0, p/p0-1., sign*sIdx);*/
 714
 715     // Fill Debug stream for MC track
 716     if(DebugLevel()>=4){
 717       Int_t det(fTracklet->GetDetector());
 718       (*DebugStream()) << "MC"
 719         << "det="     << det
 720         << "pdg="     << pdg
 721         << "sgn="     << sign
 722         << "pt="      << pt0
 723         << "x="       << x0
 724         << "y="       << y0
 725         << "z="       << z0
 726         << "dydx="    << dydx0
 727         << "dzdx="    << dzdx0
 728         << "\n";
 729
 730       // Fill Debug stream for Kalman track
 731       (*DebugStream()) << "MCtrack"
 732         << "pt="      << pt
 733         << "x="       << x
 734         << "y="       << y
 735         << "z="       << z
 736         << "dydx="    << dydx
 737         << "dzdx="    << dzdx
 738         << "s2y="     << covR[0]
 739         << "s2z="     << covR[2]
 740         << "\n";
 741     }
 742
 743     // tracklet residuals
 744     dydx = fTracklet->GetYfit(1) + tilt*dzdx0;
 745     dzdx = fTracklet->GetZfit(1);
 746     y  = fTracklet->GetYfit(0);
 747     dy = y - ymc;
 748     z  = fTracklet->GetZfit(0);
 749     dz = z - zmc;
 750     val[kYrez] = dy - dz*tilt;
 751     val[kPrez] = TMath::ATan((dydx - dydx0)/(1.+ dydx*dydx0))*TMath::RadToDeg();
 752     val[kZrez] = dz + dy*tilt;
 753 //      val[kNdim] = pt/pt0-1.;
 754     if((H = (THnSparse*)fContainer->At(kMCtracklet))) H->Fill(val);
 755
 756
 757     // Fill Debug stream for tracklet
 758     if(DebugLevel()>=4){
 759       Float_t s2y = fTracklet->GetS2Y();
 760       Float_t s2z = fTracklet->GetS2Z();
 761       (*DebugStream()) << "MCtracklet"
 762         << "rc="    << rc
 763         << "x="     << x
 764         << "y="     << y
 765         << "z="     << z
 766         << "dydx="  << dydx
 767         << "s2y="   << s2y
 768         << "s2z="   << s2z
 769         << "\n";
 770     }
 771
 772     AliTRDpadPlane *pp = geo->GetPadPlane(ily, AliTRDgeometry::GetStack(fTracklet->GetDetector()));
 773     Float_t zr0 = pp->GetRow0() + pp->GetAnodeWireOffset();
 774     //Double_t exb = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
 775
 776     H = (THnSparse*)fContainer->At(kMCcluster);
 777     val[kPt]  = TMath::ATan(dydx0)*TMath::RadToDeg();
 778     //Float_t corr = 1./TMath::Sqrt(1.+dydx0*dydx0+dzdx0*dzdx0);
 779     Int_t row0(-1);
 780     Float_t padCorr(tilt*fTracklet->GetPadLength());
 781     fTracklet->ResetClusterIter(kTRUE);
 782     while((c = fTracklet->NextCluster())){
 783       if(row0<0) row0 = c->GetPadRow();
 784       x = c->GetX();//+fXcorr[c->GetDetector()][c->GetLocalTimeBin()];
 785       y = c->GetY()  + padCorr*(c->GetPadRow() - row0);
 786       z = c->GetZ();
 787       dx = x0 - x;
 788       ymc= y0 - dx*dydx0;
 789       zmc= z0 - dx*dzdx0;
 790       dy = y - ymc;
 791       dz = z - zmc;
 792       val[kYrez] = dy - dz*tilt;
 793       val[kPrez] = dx;
 794       val[kZrez] = 0.; AliTRDcluster *cc(NULL); Int_t ic(0), tb(c->GetLocalTimeBin()); Float_t  q(TMath::Abs(c->GetQ()));
 795       if((cc = fTracklet->GetClusters(tb-1))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
 796       if((cc = fTracklet->GetClusters(tb-2))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
 797       if(ic) val[kZrez] /= (ic*q);
 798       if(H) H->Fill(val);
 799
 800
 801       // Fill calibration container
 802       Float_t d = zr0 - zmc;
 803       d -= ((Int_t)(2 * d)) / 2.0;
 804       if (d > 0.25) d  = 0.5 - d;
 805       AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
 806       clInfo->SetCluster(c);
 807       clInfo->SetMC(pdg, label);
 808       clInfo->SetGlobalPosition(ymc, zmc, dydx0, dzdx0);
 809       clInfo->SetResolution(dy);
 810       clInfo->SetAnisochronity(d);
 811       clInfo->SetDriftLength(dx);
 812       clInfo->SetTilt(tilt);
 813       if(fMCcl) fMCcl->Add(clInfo);
 814       else AliDebug(1, "MCcl exchange container missing. Activate by calling \"InitExchangeContainers()\"");
 815       if(DebugLevel()>=5){
 816         if(!clInfoArr){
 817           clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
 818           clInfoArr->SetOwner(kFALSE);
 819         }
 820         clInfoArr->Add(clInfo);
 821       }
 822     }
 823     // Fill Debug Tree
 824     if(DebugLevel()>=5 && clInfoArr){
 825       (*DebugStream()) << "MCcluster"
 826         <<"clInfo.=" << clInfoArr
 827         << "\n";
 828       clInfoArr->Clear();
 829     }
 830   }
 831   if(clInfoArr) delete clInfoArr;
 832   return h;
 833 }
 834
 835
 836 //__________________________________________________________________________
 837 Int_t AliTRDresolution::GetPtBin(Float_t pt)
 838 {
 839 // Find pt bin according to local pt segmentation
 840   Int_t ipt(-1);
 841   while(ipt<AliTRDresolution::kNpt){
 842     if(pt<fgPtBin[ipt+1]) break;
 843     ipt++;
 844   }
 845   return ipt;
 846 }
 847
 848 //________________________________________________________
 849 Float_t AliTRDresolution::GetMeanStat(TH1 *h, Float_t cut, Option_t *opt)
 850 {
 851 // return mean number of entries/bin of histogram "h"
 852 // if option "opt" is given the following values are accepted:
 853 //   "<" : consider only entries less than "cut"
 854 //   ">" : consider only entries greater than "cut"
 855
 856   //Int_t dim(h->GetDimension());
 857   Int_t nbx(h->GetNbinsX()), nby(h->GetNbinsY()), nbz(h->GetNbinsZ());
 858   Double_t sum(0.); Int_t n(0);
 859   for(Int_t ix(1); ix<=nbx; ix++)
 860     for(Int_t iy(1); iy<=nby; iy++)
 861       for(Int_t iz(1); iz<=nbz; iz++){
 862         if(strcmp(opt, "")==0){sum += h->GetBinContent(ix, iy, iz); n++;}
 863         else{
 864           if(strcmp(opt, "<")==0) {
 865             if(h->GetBinContent(ix, iy, iz)<cut) {sum += h->GetBinContent(ix, iy, iz); n++;}
 866           } else if(strcmp(opt, ">")==0){
 867             if(h->GetBinContent(ix, iy, iz)>cut) {sum += h->GetBinContent(ix, iy, iz); n++;}
 868           } else {sum += h->GetBinContent(ix, iy, iz); n++;}
 869         }
 870       }
 871   return n>0?sum/n:0.;
 872 }
 873
 874 //________________________________________________________
 875 Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
 876 {
 877   //
 878   // Get the reference figures
 879   //
 880
 881   if(!gPad){
 882     AliWarning("Please provide a canvas to draw results.");
 883     return kFALSE;
 884   }
 885 /*  Int_t selection[100], n(0), selStart(0); //
 886   Int_t ly0(0), dly(5);
 887   TList *l(NULL); TVirtualPad *pad(NULL); */
 888   switch(ifig){
 889   case 0:
 890     break;
 891   }
 892   AliWarning(Form("Reference plot [%d] missing result", ifig));
 893   return kFALSE;
 894 }
 895
 896
 897 //________________________________________________________
 898 void AliTRDresolution::MakePtSegmentation(Float_t pt0, Float_t dpt)
 899 {
 900 // Build pt segments
 901   for(Int_t j(0); j<=kNpt; j++){
 902     pt0+=(TMath::Exp(j*j*dpt)-1.);
 903     fgPtBin[j]=pt0;
 904   }
 905 }
 906
 907 //________________________________________________________
 908 void AliTRDresolution::MakeSummary()
 909 {
 910 // Build summary plots
 911
 912   if(!fProj){
 913     AliError("Missing results");
 914     return;
 915   }
 916   TVirtualPad *p(NULL); TCanvas *cOut(NULL);
 917   TObjArray *arr(NULL); TH2 *h2(NULL);
 918
 919   // cluster resolution
 920   // define palette
 921   gStyle->SetPalette(1);
 922   const Int_t nClViews(11);
 923   const Char_t *vClName[nClViews] = {"ClY", "ClYn", "ClYp", "ClQn", "ClQp", "ClYXTCp", "ClYXTCn", "ClYXPh", "ClYXPh", "ClY", "ClYn"};
 924   const UChar_t vClOpt[nClViews] = {1, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0};
 925   const Int_t nTrkltViews(10);
 926   const Char_t *vTrkltName[nTrkltViews] = {"TrkltY", "TrkltYn", "TrkltYp", "TrkltPhn", "TrkltPhp", "TrkltZ", "TrkltQn", "TrkltQp", "TrkltPn", "TrkltPp"};
 927   const Int_t nTrkInViews(6);
 928   const Char_t *vTrkInName[nTrkInViews] = {"TrkInY", "TrkInYn", "TrkInYp", "TrkInZ", "TrkInPhn", "TrkInPhp"};
 929   const Int_t nTrkViews(10);
 930   const Char_t *vTrkName[nTrkViews] = {"TrkY", "TrkYn", "TrkYp", "TrkPhn", "TrkPhp", "TrkZ", "TrkQn", "TrkQp", "TrkPn", "TrkPp"};
 931   const Char_t *typName[] = {"", "MC"};
 932
 933   for(Int_t ityp(0); ityp<(HasMCdata()?2:1); ityp++){
 934     if((arr = (TObjArray*)fProj->At(ityp?kMCcluster:kCluster))){
 935       for(Int_t iview(0); iview<nClViews; iview++){
 936         cOut = new TCanvas(Form("TRDsummary%s_%sCl%02d", GetName(), typName[ityp], iview), "Cluster Resolution", 1024, 768);
 937         cOut->Divide(3,2, 1.e-5, 1.e-5);
 938         Int_t nplot(0);
 939         for(Int_t iplot(0); iplot<6; iplot++){
 940           p=cOut->cd(iplot+1);    p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
 941           if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s%d_2D", typName[ityp], vClName[iview], iplot)))) continue;
 942           nplot++;
 943           if(vClOpt[iview]==0) h2->Draw("colz");
 944           else if(vClOpt[iview]==1) DrawSigma(h2, 1.e4, 2.e2, 6.5e2, "#sigma(#Deltay) [#mum]");
 945         }
 946         if(nplot) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
 947         else delete cOut;
 948       }
 949     }
 950     // tracklet systematic
 951     if((arr = (TObjArray*)fProj->At(ityp?kMCtracklet:kTracklet))){
 952       for(Int_t iview(0); iview<nTrkltViews; iview++){
 953         cOut = new TCanvas(Form("TRDsummary%s_%sTrklt%02d", GetName(), typName[ityp], iview), "Tracklet Resolution", 1024, 768);
 954         cOut->Divide(3,2, 1.e-5, 1.e-5);
 955         Int_t nplot(0);
 956         for(Int_t iplot(0); iplot<6; iplot++){
 957           p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
 958           if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s%d_2D", typName[ityp], vTrkltName[iview], iplot)))) continue;
 959           h2->Draw("colz"); nplot++;
 960         }
 961         if(nplot) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
 962         else delete cOut;
 963       }
 964     }
 965     // trackIn systematic
 966     if((arr = (TObjArray*)fProj->At(ityp?kMCtrackIn:kTrackIn))){
 967       cOut = new TCanvas(Form("TRDsummary%s_%sTrkIn", GetName(), typName[ityp]), "Track IN Resolution", 1024, 768);
 968       cOut->Divide(3,2, 1.e-5, 1.e-5);
 969       Int_t nplot(0);
 970       for(Int_t iplot(0); iplot<nTrkInViews; iplot++){
 971         p=cOut->cd(iplot+1);    p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
 972           if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s_2D", typName[ityp], vTrkInName[iplot])))) continue;
 973           h2->Draw("colz"); nplot++;
 974       }
 975       if(nplot) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
 976       else delete cOut;
 977     }
 978   }
 979   // track MC systematic
 980   if((arr = (TObjArray*)fProj->At(kMCtrack))) {
 981     for(Int_t iview(0); iview<nTrkViews; iview++){
 982       cOut = new TCanvas(Form("TRDsummary%s_MCTrk%02d", GetName(), iview), "Track Resolution", 1024, 768);
 983       cOut->Divide(3,2, 1.e-5, 1.e-5);
 984       Int_t nplot(0);
 985       for(Int_t iplot(0); iplot<6; iplot++){
 986         p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
 987         if(!(h2 = (TH2*)arr->FindObject(Form("H%s%d_2D", vTrkName[iview], iplot)))) continue;
 988         h2->Draw("colz"); nplot++;
 989       }
 990       if(nplot) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
 991       else delete cOut;
 992     }
 993   }
 994
 995
 996   gStyle->SetPalette(1);
 997 }
 998
 999 //________________________________________________________
1000 void AliTRDresolution::DrawSigma(TH2 *h2, Float_t scale, Float_t m, Float_t M, const Char_t *title)
1001 {
1002   // Draw error bars scaled with "scale" instead of content values
1003   //use range [m,M] if limits are specified
1004
1005   if(!h2) return;
1006   TH2 *h2e = (TH2F*)h2->Clone(Form("%s_E", h2->GetName()));
1007   h2e->SetContour(10);
1008   if(M>m) h2e->GetZaxis()->SetRangeUser(m, M);
1009   if(title) h2e->GetZaxis()->SetTitle(title);
1010
1011   for(Int_t ix(1); ix<=h2->GetNbinsX(); ix++){
1012     for(Int_t iy(1); iy<=h2->GetNbinsY(); iy++){
1013       if(h2->GetBinContent(ix, iy)<-100.) continue;
1014       Float_t v(scale*h2->GetBinError(ix, iy));
1015       if(M>m && v<m) v=m+TMath::Abs((M-m)*1.e-3);
1016       h2e->SetBinContent(ix, iy, v);
1017     }
1018   }
1019   h2e->Draw("colz");
1020 }
1021
1022 //________________________________________________________
1023 void AliTRDresolution::GetRange(TH2 *h2, Char_t mod, Float_t *range)
1024 {
1025 // Returns the range of the bulk of data in histogram h2. Removes outliers.
1026 // The "range" vector should be initialized with 2 elements
1027 // Option "mod" can be any of
1028 //   - 0 : gaussian like distribution
1029 //   - 1 : tailed distribution
1030
1031   Int_t nx(h2->GetNbinsX())
1032       , ny(h2->GetNbinsY())
1033       , n(nx*ny);
1034   Double_t *data=new Double_t[n];
1035   for(Int_t ix(1), in(0); ix<=nx; ix++){
1036     for(Int_t iy(1); iy<=ny; iy++)
1037       data[in++] = h2->GetBinContent(ix, iy);
1038   }
1039   Double_t mean, sigm;
1040   AliMathBase::EvaluateUni(n, data, mean, sigm, Int_t(n*.8));
1041
1042   range[0]=mean-3.*sigm; range[1]=mean+3.*sigm;
1043   if(mod==1) range[0]=TMath::Max(Float_t(1.e-3), range[0]);
1044   AliDebug(2, Form("h[%s] range0[%f %f]", h2->GetName(), range[0], range[1]));
1045   TH1S h1("h1SF0", "", 100, range[0], range[1]);
1046   h1.FillN(n,data,0);
1047   delete [] data;
1048
1049   switch(mod){
1050   case 0:// gaussian distribution
1051   {
1052     TF1 fg("fg", "gaus", mean-3.*sigm, mean+3.*sigm);
1053     h1.Fit(&fg, "QN");
1054     mean=fg.GetParameter(1); sigm=fg.GetParameter(2);
1055     range[0] = mean-2.5*sigm;range[1] = mean+2.5*sigm;
1056     AliDebug(2, Form("     rangeG[%f %f]", range[0], range[1]));
1057     break;
1058   }
1059   case 1:// tailed distribution
1060   {
1061     Int_t bmax(h1.GetMaximumBin());
1062     Int_t jBinMin(1), jBinMax(100);
1063     for(Int_t ibin(bmax); ibin--;){
1064       if(h1.GetBinContent(ibin)<1.){
1065         jBinMin=ibin; break;
1066       }
1067     }
1068     for(Int_t ibin(bmax); ibin++;){
1069       if(h1.GetBinContent(ibin)<1.){
1070         jBinMax=ibin; break;
1071       }
1072     }
1073     range[0]=h1.GetBinCenter(jBinMin); range[1]=h1.GetBinCenter(jBinMax);
1074     AliDebug(2, Form("     rangeT[%f %f]", range[0], range[1]));
1075     break;
1076   }
1077   }
1078
1079   return;
1080 }
1081
1082
1083 //________________________________________________________
1084 Bool_t AliTRDresolution::MakeProjectionCluster(Bool_t mc)
1085 {
1086 // Analyse cluster
1087   const Int_t kNcontours(9);
1088   const Int_t kNstat(300);
1089   Int_t cidx=mc?kMCcluster:kCluster;
1090   if(fProj && fProj->At(cidx)) return kTRUE;
1091   if(!fContainer){
1092     AliError("Missing data container.");
1093     return kFALSE;
1094   }
1095   THnSparse *H(NULL);
1096   if(!(H = (THnSparse*)fContainer->At(cidx))){
1097     AliError(Form("Missing/Wrong data @ %d.", cidx));
1098     return kFALSE;
1099   }
1100   Int_t ndim(H->GetNdimensions());
1101   Int_t coord[kNdim]; memset(coord, 0, sizeof(Int_t) * kNdim); Double_t v = 0.;
1102   TAxis *aa[kNdim], *as(NULL), *apt(NULL); memset(aa, 0, sizeof(TAxis*) * kNdim);
1103   for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1104   if(ndim > kPt) apt = H->GetAxis(kPt);
1105   if(ndim > kSpeciesChgRC) as  = H->GetAxis(kSpeciesChgRC);
1106   // build list of projections
1107   const Int_t nsel(12), npsel(5);
1108   // define rebinning strategy
1109   const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
1110   AliTRDresolutionProjection hp[kClNproj], *php[nsel][npsel]; memset(php, 0, nsel*npsel*sizeof(AliTRDresolutionProjection*));
1111   Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1112   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1113     isel++; // new selection
1114     hp[ih].Build(Form("H%sClY%d", mc?"MC":"", ily), Form("Clusters :: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1115     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1116       php[isel][np[isel]++] = &hp[ih++];
1117     hp[ih].Build(Form("H%sClYn%d", mc?"MC":"", ily), Form("Clusters[-]:: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1118     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1119       php[isel][np[isel]++] = &hp[ih++];
1120     hp[ih].Build(Form("H%sClQn%d", mc?"MC":"", ily), Form("Clusters[-]:: Charge distribution ly%d", ily), kEta, kPhi, kSpeciesChgRC, aa);
1121     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1122     hp[ih].SetShowRange(20., 40.);
1123       php[isel][np[isel]++] = &hp[ih++];
1124     hp[ih].Build(Form("H%sClYXTCn%d", mc?"MC":"", ily), Form("Clusters[-]:: r-#phi(x,TC) residuals ly%d", ily), kPrez, kZrez, kYrez, aa);
1125 //    hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1126       php[isel][np[isel]++] = &hp[ih++];
1127     hp[ih].Build(Form("H%sClYXPh%d", mc?"MC":"", ily), Form("Clusters :: r-#phi(x,#Phi) residuals ly%d", ily), kPrez, kPt, kYrez, aa);
1128 //    hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1129       php[isel][np[isel]++] = &hp[ih++];
1130     isel++; // new selection
1131       php[isel][np[isel]++] = &hp[ih-5]; // relink HClY
1132     hp[ih].Build(Form("H%sClYp%d", mc?"MC":"", ily), Form("Clusters[+]:: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1133     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1134       php[isel][np[isel]++] = &hp[ih++];
1135     hp[ih].Build(Form("H%sClQp%d", mc?"MC":"", ily), Form("Clusters[+]:: Charge distribution ly%d", ily), kEta, kPhi, kSpeciesChgRC, aa);
1136     hp[ih].SetShowRange(20., 40.);
1137     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1138       php[isel][np[isel]++] = &hp[ih++];
1139     hp[ih].Build(Form("H%sClYXTCp%d", mc?"MC":"", ily), Form("Clusters[+]:: r-#phi(x,TC) residuals ly%d", ily), kPrez, kZrez, kYrez, aa);
1140 //    hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1141       php[isel][np[isel]++] = &hp[ih++];
1142       php[isel][np[isel]++] = &hp[ih-4]; // relink HClYXPh
1143   }
1144
1145   Int_t ly(0), ch(0), rcBin(as?as->FindBin(0.):-1), chBin(apt?apt->FindBin(0.):-1);
1146   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1147     v = H->GetBinContent(ib, coord); if(v<1.) continue;
1148     ly = coord[kBC]-1;
1149     // RC selection
1150     if(rcBin>0 && coord[kSpeciesChgRC] == rcBin) continue;
1151
1152     // charge selection
1153     ch = 0; // [-] track
1154     if(chBin>0 && coord[kPt] > chBin) ch = 1;  // [+] track
1155
1156     isel = ly*2+ch;
1157     for(Int_t jh(0); jh<np[isel]; jh++) php[isel][jh]->Increment(coord, v);
1158   }
1159
1160   if(!fProj){
1161     AliInfo("Building array of projections ...");
1162     fProj = new TObjArray(kNclasses); fProj->SetOwner(kTRUE);
1163   }
1164   TObjArray *arr(NULL);
1165   fProj->AddAt(arr = new TObjArray(kClNproj), cidx);
1166
1167   TH2 *h2(NULL);
1168   for(; ih--; ){
1169     if(!hp[ih].fH) continue;
1170     Int_t mid(1), nstat(kNstat);
1171     if(strchr(hp[ih].fH->GetName(), 'Q')){ mid=2; nstat=300;}
1172     if(!(h2 = hp[ih].Projection2D(nstat, kNcontours, mid))) continue;
1173     arr->AddAt(h2, ih);
1174   }
1175
1176   return kTRUE;
1177 }
1178
1179 //________________________________________________________
1180 Bool_t AliTRDresolution::MakeProjectionTracklet(Bool_t mc)
1181 {
1182 // Analyse tracklet
1183   const Int_t kNcontours(9);
1184   const Int_t kNstat(100);
1185   Int_t cidx=mc?kMCtracklet:kTracklet;
1186   if(fProj && fProj->At(cidx)) return kTRUE;
1187   if(!fContainer){
1188     AliError("Missing data container.");
1189     return kFALSE;
1190   }
1191   THnSparse *H(NULL);
1192   if(!(H = (THnSparse*)fContainer->At(cidx))){
1193     AliError(Form("Missing/Wrong data @ %d.", cidx));
1194     return kFALSE;
1195   }
1196   Int_t ndim(H->GetNdimensions());
1197   Int_t coord[kNdim+1]; memset(coord, 0, sizeof(Int_t) * (kNdim+1)); Double_t v = 0.;
1198   TAxis *aa[kNdim+1], *as(NULL); memset(aa, 0, sizeof(TAxis*) * (kNdim+1));
1199   for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1200   if(ndim > kSpeciesChgRC) as = H->GetAxis(kSpeciesChgRC);
1201   // build list of projections
1202   const Int_t nsel(18), npsel(6);
1203   // define rebinning strategy
1204   const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
1205   AliTRDresolutionProjection hp[kTrkltNproj], *php[nsel][npsel]; memset(php, 0, nsel*npsel*sizeof(AliTRDresolutionProjection*));
1206   Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1207   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1208     isel++; // new selection
1209     hp[ih].Build(Form("H%sTrkltY%d", mc?"MC":"", ily), Form("Tracklets   :: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1210     hp[ih].SetShowRange(-0.03,0.03);
1211     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1212       php[isel][np[isel]++] = &hp[ih++];
1213     hp[ih].Build(Form("H%sTrkltYn%d", mc?"MC":"", ily), Form("Tracklets[-]:: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1214     hp[ih].SetShowRange(-0.03,0.03);
1215     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1216       php[isel][np[isel]++] = &hp[ih++];
1217     hp[ih].Build(Form("H%sTrkltPhn%d", mc?"MC":"", ily), Form("Tracklets[-]:: #Delta#phi residuals ly%d", ily), kEta, kPhi, kPrez, aa);
1218     hp[ih].SetShowRange(-0.5,0.5);
1219     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1220       php[isel][np[isel]++] = &hp[ih++];
1221     hp[ih].Build(Form("H%sTrkltPn%d", mc?"MC":"", ily), Form("Tracklets[-]:: Momentum distribution ly%d", ily), kEta, kPhi, kPt, aa);
1222     hp[ih].SetShowRange(6.,12.);
1223     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1224       php[isel][np[isel]++] = &hp[ih++];
1225     hp[ih].Build(Form("H%sTrkltYPn%d", mc?"MC":"", ily), Form("Tracklets[-]:: r-#phi/p_{t} residuals ly%d", ily), kPt, kPhi, kYrez, aa);
1226       php[isel][np[isel]++] = &hp[ih++];
1227     hp[ih].Build(Form("H%sTrkltQn%d", mc?"MC":"", ily), Form("Tracklets[-]:: dQdl ly%d", ily), kEta, kPhi, kNdim, aa);
1228     hp[ih].SetShowRange(700.,1100.);
1229     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1230       php[isel][np[isel]++] = &hp[ih++];
1231     isel++; // new selection
1232     php[isel][np[isel]++] = &hp[ih-6]; // relink first histo
1233     hp[ih].Build(Form("H%sTrkltYp%d", mc?"MC":"", ily), Form("Tracklets[+]:: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1234     hp[ih].SetShowRange(-0.03,0.03);
1235     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1236       php[isel][np[isel]++] = &hp[ih++];
1237     hp[ih].Build(Form("H%sTrkltPhp%d", mc?"MC":"", ily), Form("Tracklets[+]:: #Delta#phi residuals ly%d", ily), kEta, kPhi, kPrez, aa);
1238     hp[ih].SetShowRange(-0.5,0.5);
1239     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1240       php[isel][np[isel]++] = &hp[ih++];
1241     hp[ih].Build(Form("H%sTrkltPp%d", mc?"MC":"", ily), Form("Tracklets[+]:: Momentum distribution ly%d", ily), kEta, kPhi, kPt, aa);
1242     hp[ih].SetShowRange(6.,12.);
1243     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1244       php[isel][np[isel]++] = &hp[ih++];
1245     hp[ih].Build(Form("H%sTrkltYPp%d", mc?"MC":"", ily), Form("Tracklets[+]:: r-#phi/p_{t} residuals ly%d", ily), kPt, kPhi, kYrez, aa);
1246       php[isel][np[isel]++] = &hp[ih++];
1247     hp[ih].Build(Form("H%sTrkltQp%d", mc?"MC":"", ily), Form("Tracklets[+]:: dQdl ly%d", ily), kEta, kPhi, kNdim, aa);
1248     hp[ih].SetShowRange(700.,1100.);
1249     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1250       php[isel][np[isel]++] = &hp[ih++];
1251     isel++; // new selection
1252     hp[ih].Build(Form("H%sTrkltZ%d", mc?"MC":"", ily), Form("Tracklets[RC]:: z residuals ly%d", ily), kEta, kPhi, kZrez, aa);
1253     hp[ih].SetShowRange(-0.1,0.1);
1254     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1255       php[isel][np[isel]++] = &hp[ih++];
1256   }
1257
1258   Int_t ly(0), ch(0), rcBin(as?as->FindBin(0.):-1);
1259   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1260     v = H->GetBinContent(ib, coord);
1261     if(v<1.) continue;
1262     ly = coord[kBC]-1; // layer selection
1263     // charge selection
1264     ch = 0; // [-] track
1265     if(rcBin>0){ // debug mode in which species are also saved
1266       if(coord[kSpeciesChgRC] > rcBin) ch = 1;  // [+] track
1267       else if(coord[kSpeciesChgRC] == rcBin) ch = 2;  // [RC] track
1268     }
1269     isel = ly*3+ch;
1270     for(Int_t jh(0); jh<np[isel]; jh++) php[isel][jh]->Increment(coord, v);
1271   }
1272
1273   if(!fProj){
1274     AliInfo("Building array of projections ...");
1275     fProj = new TObjArray(kNclasses); fProj->SetOwner(kTRUE);
1276   }
1277   TObjArray *arr(NULL);
1278   fProj->AddAt(arr = new TObjArray(kTrkltNproj), cidx);
1279
1280   TH2 *h2(NULL);
1281   for(; ih--; ){
1282     if(!hp[ih].fH) continue;
1283     Int_t mid(0), nstat(kNstat);
1284     if(strchr(hp[ih].fH->GetName(), 'Q')){ mid=2; nstat=200;}
1285     if(!(h2 = hp[ih].Projection2D(nstat, kNcontours, mid))) continue;
1286     arr->AddAt(h2, ih);
1287   }
1288   return kTRUE;
1289 }
1290
1291 //________________________________________________________
1292 Bool_t AliTRDresolution::MakeProjectionTrackIn(Bool_t mc)
1293 {
1294 // Analyse track in
1295
1296   const Int_t kNcontours(9);
1297   const Int_t kNstat(30);
1298   Int_t cidx=mc?kMCtrackIn:kTrackIn;
1299   if(fProj && fProj->At(cidx)) return kTRUE;
1300   if(!fContainer){
1301     AliError("Missing data container.");
1302     return kFALSE;
1303   }
1304   THnSparse *H(NULL);
1305   if(!(H = (THnSparse*)fContainer->At(cidx))){
1306     AliError(Form("Missing/Wrong data @ %d.", Int_t(cidx)));
1307     return kFALSE;
1308   }
1309
1310   Int_t coord[kNdim]; memset(coord, 0, sizeof(Int_t) * kNdim); Double_t v = 0.;
1311   Int_t ndim(H->GetNdimensions());
1312   TAxis *aa[kNdim+1], *as(NULL); memset(aa, 0, sizeof(TAxis*) * (kNdim+1));
1313   for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1314   if(ndim > kSpeciesChgRC) as = H->GetAxis(kSpeciesChgRC);
1315   // build list of projections
1316   const Int_t nsel(3), npsel(4);
1317   // define rebinning strategy
1318   const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
1319   const Int_t nPtPhi(2); Int_t rebinPtPhiX[nEtaPhi] = {1, 1}, rebinPtPhiY[nEtaPhi] = {2, 5};
1320   AliTRDresolutionProjection hp[kTrkInNproj], *php[nsel][npsel]; memset(php, 0, nsel*npsel*sizeof(AliTRDresolutionProjection*));
1321   Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1322   // define list of projections
1323   isel++;  // negative tracks
1324   hp[ih].Build(Form("H%sTrkInY", mc?"MC":""), "TrackIn :: r-#phi residuals", kEta, kPhi, kYrez, aa);
1325   hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1326     php[isel][np[isel]++] = &hp[ih++];
1327   hp[ih].Build(Form("H%sTrkInYn", mc?"MC":""), "TrackIn[-]:: r-#phi residuals", kEta, kPhi, kYrez, aa);
1328   hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1329     php[isel][np[isel]++] = &hp[ih++];
1330   hp[ih].Build(Form("H%sTrkInPhn", mc?"MC":""), "TrackIn[-]:: #Delta#phi residuals", kEta, kPhi, kPrez, aa);
1331   hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1332     php[isel][np[isel]++] = &hp[ih++];
1333   hp[ih].Build(Form("H%sTrkInYPn", mc?"MC":""), "TrackIn[-]:: r-#phi/p_{t} residuals", kPt, kPhi, kYrez, aa);
1334   hp[ih].SetRebinStrategy(nPtPhi, rebinPtPhiX, rebinPtPhiY);
1335     php[isel][np[isel]++] = &hp[ih++];
1336   isel++; // positive tracks
1337   php[isel][np[isel]++] = &hp[ih-4]; // relink first histo
1338   hp[ih].Build(Form("H%sTrkInYp", mc?"MC":""), "TrackIn[+]:: r-#phi residuals", kEta, kPhi, kYrez, aa);
1339   hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1340     php[isel][np[isel]++] = &hp[ih++];
1341   hp[ih].Build(Form("H%sTrkInPhp", mc?"MC":""), "TrackIn[+]:: #Delta#phi residuals", kEta, kPhi, kPrez, aa);
1342   hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1343     php[isel][np[isel]++] = &hp[ih++];
1344   hp[ih].Build(Form("H%sTrkInYPp", mc?"MC":""), "TrackIn[+]:: r-#phi/p_{t} residuals", kPt, kPhi, kYrez, aa);
1345   hp[ih].SetRebinStrategy(nPtPhi, rebinPtPhiX, rebinPtPhiY);
1346     php[isel][np[isel]++] = &hp[ih++];
1347   isel++; // RC tracks
1348   hp[ih].Build(Form("H%sTrkInZ", mc?"MC":""), "TrackIn[RC]:: z residuals", kEta, kPhi, kZrez, aa);
1349   hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1350     php[isel][np[isel]++] = &hp[ih++];
1351
1352   // fill projections
1353   Int_t ch(0), rcBin(as?as->FindBin(0.):-1);
1354   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1355     v = H->GetBinContent(ib, coord);
1356     if(v<1.) continue;
1357     if(coord[kBC]>1) continue; // bunch cross cut
1358     // charge selection
1359     ch = 0; // [-] track
1360     if(rcBin>0){ // debug mode in which species are also saved
1361       if(coord[kSpeciesChgRC] > rcBin) ch = 1;  // [+] track
1362       else if(coord[kSpeciesChgRC] == rcBin) ch = 2;  // [RC] track
1363     }
1364     for(Int_t jh(0); jh<np[ch]; jh++) php[ch][jh]->Increment(coord, v);
1365   }
1366   if(!fProj){
1367     AliInfo("Building array of projections ...");
1368     fProj = new TObjArray(kNclasses); fProj->SetOwner(kTRUE);
1369   }
1370   TObjArray *arr(NULL);
1371   fProj->AddAt(arr = new TObjArray(kTrkInNproj), cidx);
1372
1373   TH2 *h2(NULL);
1374   for(; ih--; ){
1375     if(!hp[ih].fH) continue;
1376     if(!(h2 = hp[ih].Projection2D(kNstat, kNcontours))) continue;
1377     arr->AddAt(h2, ih);
1378   }
1379   return kTRUE;
1380 }
1381
1382
1383 //________________________________________________________
1384 Bool_t AliTRDresolution::MakeProjectionTrack()
1385 {
1386 // Analyse tracklet
1387   const Int_t kNcontours(9);
1388   const Int_t kNstat(100);
1389   Int_t cidx(kMCtrack);
1390   if(fProj && fProj->At(cidx)) return kTRUE;
1391   if(!fContainer){
1392     AliError("Missing data container.");
1393     return kFALSE;
1394   }
1395   THnSparse *H(NULL);
1396   if(!(H = (THnSparse*)fContainer->At(cidx))){
1397     AliError(Form("Missing/Wrong data @ %d.", cidx));
1398     return kFALSE;
1399   }
1400   Int_t ndim(H->GetNdimensions());
1401   Int_t coord[kNdim+1]; memset(coord, 0, sizeof(Int_t) * (kNdim+1)); Double_t v = 0.;
1402   TAxis *aa[kNdim+1], *as(NULL); memset(aa, 0, sizeof(TAxis*) * (kNdim+1));
1403   for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1404   if(ndim > kSpeciesChgRC) as = H->GetAxis(kSpeciesChgRC);
1405   // build list of projections
1406   const Int_t nsel(18), npsel(6);
1407   // define rebinning strategy
1408   const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
1409   AliTRDresolutionProjection hp[kTrkNproj], *php[nsel][npsel]; memset(php, 0, nsel*npsel*sizeof(AliTRDresolutionProjection*));
1410   Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1411   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1412     isel++; // new selection
1413     hp[ih].Build(Form("HTrkY%d", ily), Form("Tracks   :: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1414     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1415       php[isel][np[isel]++] = &hp[ih++];
1416     hp[ih].Build(Form("HTrkYn%d", ily), Form("Tracks[-]:: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1417     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1418       php[isel][np[isel]++] = &hp[ih++];
1419     hp[ih].Build(Form("HTrkPhn%d", ily), Form("Tracks[-]:: #Delta#phi residuals ly%d", ily), kEta, kPhi, kPrez, aa);
1420     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1421       php[isel][np[isel]++] = &hp[ih++];
1422     hp[ih].Build(Form("HTrkPn%d", ily), Form("Tracks[-]:: Momentum distribution ly%d", ily), kEta, kPhi, kPt, aa);
1423     hp[ih].SetShowRange(6.,12.);
1424     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1425       php[isel][np[isel]++] = &hp[ih++];
1426     hp[ih].Build(Form("HTrkYPn%d", ily), Form("Tracks[-]:: r-#phi/p_{t} residuals ly%d", ily), kPt, kPhi, kYrez, aa);
1427       php[isel][np[isel]++] = &hp[ih++];
1428     hp[ih].Build(Form("HTrkQn%d", ily), Form("Tracks[-]:: #Deltap_{t}/p_{t} ly%d", ily), kEta, kPhi, kNdim, aa);
1429     //hp[ih].SetShowRange(700.,1100.);
1430     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1431       php[isel][np[isel]++] = &hp[ih++];
1432     isel++; // new selection
1433     php[isel][np[isel]++] = &hp[ih-6]; // relink first histo
1434     hp[ih].Build(Form("HTrkYp%d", ily), Form("Tracks[+]:: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1435     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1436       php[isel][np[isel]++] = &hp[ih++];
1437     hp[ih].Build(Form("HTrkPhp%d", ily), Form("Tracks[+]:: #Delta#phi residuals ly%d", ily), kEta, kPhi, kPrez, aa);
1438     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1439       php[isel][np[isel]++] = &hp[ih++];
1440     hp[ih].Build(Form("HTrkPp%d", ily), Form("Tracks[+]:: Momentum distribution ly%d", ily), kEta, kPhi, kPt, aa);
1441     hp[ih].SetShowRange(6.,12.);
1442     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1443       php[isel][np[isel]++] = &hp[ih++];
1444     hp[ih].Build(Form("HTrkYPp%d", ily), Form("Tracks[+]:: r-#phi/p_{t} residuals ly%d", ily), kPt, kPhi, kYrez, aa);
1445       php[isel][np[isel]++] = &hp[ih++];
1446     hp[ih].Build(Form("HTrkQp%d", ily), Form("Tracks[+]:: #Deltap_{t}/p_{t} ly%d", ily), kEta, kPhi, kNdim, aa);
1447     //hp[ih].SetShowRange(700.,1100.);
1448     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1449       php[isel][np[isel]++] = &hp[ih++];
1450     isel++; // new selection
1451     hp[ih].Build(Form("HTrkZ%d", ily), Form("Tracks[RC]:: z residuals ly%d", ily), kEta, kPhi, kZrez, aa);
1452     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1453       php[isel][np[isel]++] = &hp[ih++];
1454   }
1455
1456   Int_t ly(0), ch(0), rcBin(as?as->FindBin(0.):-1);
1457   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1458     v = H->GetBinContent(ib, coord);
1459     if(v<1.) continue;
1460     ly = coord[kBC]-1; // layer selection
1461     // charge selection
1462     ch = 0; // [-] track
1463     if(rcBin>0){ // debug mode in which species are also saved
1464       if(coord[kSpeciesChgRC] > rcBin) ch = 1;  // [+] track
1465       else if(coord[kSpeciesChgRC] == rcBin) ch = 2;  // [RC] track
1466     }
1467     isel = ly*3+ch;
1468     for(Int_t jh(0); jh<np[isel]; jh++) php[isel][jh]->Increment(coord, v);
1469   }
1470
1471   if(!fProj){
1472     AliInfo("Building array of projections ...");
1473     fProj = new TObjArray(kNclasses); fProj->SetOwner(kTRUE);
1474   }
1475   TObjArray *arr(NULL);
1476   fProj->AddAt(arr = new TObjArray(kTrkNproj), cidx);
1477
1478   TH2 *h2(NULL);
1479   for(; ih--; ){
1480     if(!hp[ih].fH) continue;
1481     if(!(h2 = hp[ih].Projection2D(kNstat, kNcontours))) continue;
1482     arr->AddAt(h2, ih);
1483   }
1484   return kTRUE;
1485 }
1486
1487 //________________________________________________________
1488 Bool_t AliTRDresolution::PostProcess()
1489 {
1490 // Fit, Project, Combine, Extract values from the containers filled during execution
1491
1492   if (!fContainer) {
1493     AliError("ERROR: list not available");
1494     return kFALSE;
1495   }
1496
1497   //PROCESS EXPERIMENTAL DISTRIBUTIONS
1498   // Clusters residuals
1499   if(!MakeProjectionCluster()) return kFALSE;
1500   fNRefFigures = 3;
1501   // Tracklet residual/pulls
1502   if(!MakeProjectionTracklet()) return kFALSE;
1503   fNRefFigures = 7;
1504   // TRDin residual/pulls
1505   if(!MakeProjectionTrackIn()) return kFALSE;
1506   fNRefFigures = 11;
1507
1508   if(!HasMCdata()) return kTRUE;
1509   //PROCESS MC RESIDUAL DISTRIBUTIONS
1510
1511   // CLUSTER Y RESOLUTION/PULLS
1512   if(!MakeProjectionCluster(kTRUE)) return kFALSE;
1513   fNRefFigures = 17;
1514
1515   // TRACKLET RESOLUTION/PULLS
1516   if(!MakeProjectionTracklet(kTRUE)) return kFALSE;
1517   fNRefFigures = 21;
1518
1519   // TRACK RESOLUTION/PULLS
1520   if(!MakeProjectionTrack()) return kFALSE;
1521   fNRefFigures+=16;
1522
1523   // TRACK TRDin RESOLUTION/PULLS
1524   if(!MakeProjectionTrackIn(kTRUE)) return kFALSE;
1525   fNRefFigures+=8;
1526
1527   return kTRUE;
1528 }
1529
1530
1531 //________________________________________________________
1532 void AliTRDresolution::Terminate(Option_t *opt)
1533 {
1534   AliTRDrecoTask::Terminate(opt);
1535   if(HasPostProcess()) PostProcess();
1536 }
1537
1538 //________________________________________________________
1539 void AliTRDresolution::AdjustF1(TH1 *h, TF1 *f)
1540 {
1541 // Helper function to avoid duplication of code
1542 // Make first guesses on the fit parameters
1543
1544   // find the intial parameters of the fit !! (thanks George)
1545   Int_t nbinsy = Int_t(.5*h->GetNbinsX());
1546   Double_t sum = 0.;
1547   for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
1548   f->SetParLimits(0, 0., 3.*sum);
1549   f->SetParameter(0, .9*sum);
1550   Double_t rms = h->GetRMS();
1551   f->SetParLimits(1, -rms, rms);
1552   f->SetParameter(1, h->GetMean());
1553
1554   f->SetParLimits(2, 0., 2.*rms);
1555   f->SetParameter(2, rms);
1556   if(f->GetNpar() <= 4) return;
1557
1558   f->SetParLimits(3, 0., sum);
1559   f->SetParameter(3, .1*sum);
1560
1561   f->SetParLimits(4, -.3, .3);
1562   f->SetParameter(4, 0.);
1563
1564   f->SetParLimits(5, 0., 1.e2);
1565   f->SetParameter(5, 2.e-1);
1566 }
1567
1568 //________________________________________________________
1569 TObjArray* AliTRDresolution::BuildMonitorContainerCluster(const char* name, Bool_t expand, Float_t range)
1570 {
1571 // Build performance histograms for AliTRDcluster.vs TRD track or MC
1572 //  - y reziduals/pulls
1573
1574   TObjArray *arr = new TObjArray(2);
1575   arr->SetName(name); arr->SetOwner();
1576   TH1 *h(NULL); char hname[100], htitle[300];
1577
1578   // tracklet resolution/pull in y direction
1579   snprintf(hname, 100, "%s_%s_Y", GetNameId(), name);
1580   snprintf(htitle, 300, "Y res for \"%s\" @ %s;tg(#phi);#Delta y [cm];%s", GetNameId(), name, "Detector");
1581   Float_t rr = range<0.?fDyRange:range;
1582   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1583     Int_t nybins=50;
1584     if(expand) nybins*=2;
1585     h = new TH3S(hname, htitle,
1586                  48, -.48, .48,            // phi
1587                  60, -rr, rr,              // dy
1588                  nybins, -0.5, nybins-0.5);// segment
1589   } else h->Reset();
1590   arr->AddAt(h, 0);
1591   snprintf(hname, 100, "%s_%s_YZpull", GetNameId(), name);
1592   snprintf(htitle, 300, "YZ pull for \"%s\" @ %s;%s;#Delta y  / #sigma_{y};#Delta z  / #sigma_{z}", GetNameId(), name, "Detector");
1593   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1594     h = new TH3S(hname, htitle, 540, -0.5, 540-0.5, 100, -4.5, 4.5, 100, -4.5, 4.5);
1595   } else h->Reset();
1596   arr->AddAt(h, 1);
1597
1598   return arr;
1599 }
1600
1601 //________________________________________________________
1602 TObjArray* AliTRDresolution::BuildMonitorContainerTracklet(const char* name, Bool_t expand)
1603 {
1604 // Build performance histograms for AliExternalTrackParam.vs TRD tracklet
1605 //  - y reziduals/pulls
1606 //  - z reziduals/pulls
1607 //  - phi reziduals
1608   TObjArray *arr = BuildMonitorContainerCluster(name, expand, 0.05);
1609   arr->Expand(5);
1610   TH1 *h(NULL); char hname[100], htitle[300];
1611
1612   // tracklet resolution/pull in z direction
1613   snprintf(hname, 100, "%s_%s_Z", GetNameId(), name);
1614   snprintf(htitle, 300, "Z res for \"%s\" @ %s;tg(#theta);#Delta z [cm]", GetNameId(), name);
1615   if(!(h = (TH2S*)gROOT->FindObject(hname))){
1616     h = new TH2S(hname, htitle, 50, -1., 1., 100, -.05, .05);
1617   } else h->Reset();
1618   arr->AddAt(h, 2);
1619   snprintf(hname, 100, "%s_%s_Zpull", GetNameId(), name);
1620   snprintf(htitle, 300, "Z pull for \"%s\" @ %s;tg(#theta);#Delta z  / #sigma_{z};row cross", GetNameId(), name);
1621   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1622     h = new TH3S(hname, htitle, 50, -1., 1., 100, -5.5, 5.5, 2, -0.5, 1.5);
1623     h->GetZaxis()->SetBinLabel(1, "no RC");
1624     h->GetZaxis()->SetBinLabel(2, "RC");
1625   } else h->Reset();
1626   arr->AddAt(h, 3);
1627
1628   // tracklet to track phi resolution
1629   snprintf(hname, 100, "%s_%s_PHI", GetNameId(), name);
1630   snprintf(htitle, 300, "#Phi res for \"%s\" @ %s;tg(#phi);#Delta #phi [rad];%s", GetNameId(), name, "Detector");
1631   Int_t nsgms=540;
1632   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1633     h = new TH3S(hname, htitle, 48, -.48, .48, 100, -.5, .5, nsgms, -0.5, nsgms-0.5);
1634   } else h->Reset();
1635   arr->AddAt(h, 4);
1636
1637   return arr;
1638 }
1639
1640 //________________________________________________________
1641 TObjArray* AliTRDresolution::BuildMonitorContainerTrack(const char* name)
1642 {
1643 // Build performance histograms for AliExternalTrackParam.vs MC
1644 //  - y resolution/pulls
1645 //  - z resolution/pulls
1646 //  - phi resolution, snp pulls
1647 //  - theta resolution, tgl pulls
1648 //  - pt resolution, 1/pt pulls, p resolution
1649
1650   TObjArray *arr = BuildMonitorContainerTracklet(name);
1651   arr->Expand(11);
1652   TH1 *h(NULL); char hname[100], htitle[300];
1653   //TAxis *ax(NULL);
1654
1655   // snp pulls
1656   snprintf(hname, 100, "%s_%s_SNPpull", GetNameId(), name);
1657   snprintf(htitle, 300, "SNP pull for \"%s\" @ %s;tg(#phi);#Delta snp  / #sigma_{snp};entries", GetNameId(), name);
1658   if(!(h = (TH2I*)gROOT->FindObject(hname))){
1659     h = new TH2I(hname, htitle, 60, -.3, .3, 100, -4.5, 4.5);
1660   } else h->Reset();
1661   arr->AddAt(h, 5);
1662
1663   // theta resolution
1664   snprintf(hname, 100, "%s_%s_THT", GetNameId(), name);
1665   snprintf(htitle, 300, "#Theta res for \"%s\" @ %s;tg(#theta);#Delta #theta [rad];entries", GetNameId(), name);
1666   if(!(h = (TH2I*)gROOT->FindObject(hname))){
1667     h = new TH2I(hname, htitle, 100, -1., 1., 100, -5e-3, 5e-3);
1668   } else h->Reset();
1669   arr->AddAt(h, 6);
1670   // tgl pulls
1671   snprintf(hname, 100, "%s_%s_TGLpull", GetNameId(), name);
1672   snprintf(htitle, 300, "TGL pull for \"%s\" @ %s;tg(#theta);#Delta tgl  / #sigma_{tgl};entries", GetNameId(), name);
1673   if(!(h = (TH2I*)gROOT->FindObject(hname))){
1674     h = new TH2I(hname, htitle, 100, -1., 1., 100, -4.5, 4.5);
1675   } else h->Reset();
1676   arr->AddAt(h, 7);
1677
1678   const Int_t kNdpt(150);
1679   const Int_t kNspc = 2*AliPID::kSPECIES+1;
1680   Float_t lPt=0.1, lDPt=-.1, lSpc=-5.5;
1681   Float_t binsPt[kNpt+1], binsSpc[kNspc+1], binsDPt[kNdpt+1];
1682   for(Int_t i=0;i<kNpt+1; i++,lPt=TMath::Exp(i*.15)-1.) binsPt[i]=lPt;
1683   for(Int_t i=0; i<kNspc+1; i++,lSpc+=1.) binsSpc[i]=lSpc;
1684   for(Int_t i=0; i<kNdpt+1; i++,lDPt+=2.e-3) binsDPt[i]=lDPt;
1685
1686   // Pt resolution
1687   snprintf(hname, 100, "%s_%s_Pt", GetNameId(), name);
1688   snprintf(htitle, 300, "#splitline{P_{t} res for}{\"%s\" @ %s};p_{t} [GeV/c];#Delta p_{t}/p_{t}^{MC};SPECIES", GetNameId(), name);
1689   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1690     h = new TH3S(hname, htitle,
1691                  kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
1692     //ax = h->GetZaxis();
1693     //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
1694   } else h->Reset();
1695   arr->AddAt(h, 8);
1696   // 1/Pt pulls
1697   snprintf(hname, 100, "%s_%s_1Pt", GetNameId(), name);
1698   snprintf(htitle, 300, "#splitline{1/P_{t} pull for}{\"%s\" @ %s};1/p_{t}^{MC} [c/GeV];#Delta(1/p_{t})/#sigma(1/p_{t});SPECIES", GetNameId(), name);
1699   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1700     h = new TH3S(hname, htitle,
1701                  kNpt, 0., 2., 100, -4., 4., kNspc, -5.5, 5.5);
1702     //ax = h->GetZaxis();
1703     //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
1704   } else h->Reset();
1705   arr->AddAt(h, 9);
1706   // P resolution
1707   snprintf(hname, 100, "%s_%s_P", GetNameId(), name);
1708   snprintf(htitle, 300, "P res for \"%s\" @ %s;p [GeV/c];#Delta p/p^{MC};SPECIES", GetNameId(), name);
1709   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1710     h = new TH3S(hname, htitle,
1711                  kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
1712     //ax = h->GetZaxis();
1713     //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
1714   } else h->Reset();
1715   arr->AddAt(h, 10);
1716
1717   return arr;
1718 }
1719
1720
1721 //________________________________________________________
1722 TObjArray* AliTRDresolution::Histos()
1723 {
1724   //
1725   // Define histograms
1726   //
1727
1728   if(fContainer) return fContainer;
1729
1730   fContainer  = new TObjArray(kNclasses); fContainer->SetOwner(kTRUE);
1731   THnSparse *H(NULL);
1732   const Int_t nhn(100); Char_t hn[nhn]; TString st;
1733
1734   //++++++++++++++++++++++
1735   // cluster to tracklet residuals/pulls
1736   snprintf(hn, nhn, "h%s", fgPerformanceName[kCluster]);
1737   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1738     const Char_t *clTitle[kNdim] = {"layer", fgkTitle[kPhi], fgkTitle[kEta], fgkTitle[kYrez], "#Deltax [cm]", "Q</Q", "Q/angle", "#Phi [deg]"};
1739     const Int_t clNbins[kNdim]   = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], fgkNbins[kEta], fgkNbins[kYrez], 45, 10, 30, 15};
1740     const Double_t clMin[kNdim]  = {-0.5, fgkMin[kPhi], fgkMin[kEta], fgkMin[kYrez]/10., -.5, 0.1, -2., -45},
1741                    clMax[kNdim]  = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], fgkMax[kEta], fgkMax[kYrez]/10., 4., 2.1, 118., 45};
1742     st = "cluster spatial&charge resolution;";
1743     // define minimum info to be saved in non debug mode
1744     Int_t ndim=DebugLevel()>=1?kNdim:4;
1745     for(Int_t idim(0); idim<ndim; idim++){ st += clTitle[idim]; st+=";";}
1746     H = new THnSparseI(hn, st.Data(), ndim, clNbins, clMin, clMax);
1747   } else H->Reset();
1748   fContainer->AddAt(H, kCluster);
1749   //++++++++++++++++++++++
1750   // tracklet to TRD track
1751   snprintf(hn, nhn, "h%s", fgPerformanceName[kTracklet]);
1752   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1753     Char_t *trTitle[kNdim+1]; memcpy(trTitle, fgkTitle, kNdim*sizeof(Char_t*));
1754     Int_t trNbins[kNdim+1]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t));
1755     Double_t trMin[kNdim+1]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
1756     Double_t trMax[kNdim+1]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
1757     // set specific fields
1758     trMin[kYrez] = -0.45; trMax[kYrez] = 0.45;
1759     trMin[kPrez] = -4.5; trMax[kPrez] = 4.5;
1760     trMin[kZrez] = -1.5; trMax[kZrez] = 1.5;
1761     trTitle[kBC]=StrDup("layer"); trNbins[kBC] = AliTRDgeometry::kNlayer; trMin[kBC] = -0.5; trMax[kBC] = AliTRDgeometry::kNlayer-0.5;
1762     trTitle[kNdim]=StrDup("dq/dl [a.u.]"); trNbins[kNdim] = 30; trMin[kNdim] = 100.; trMax[kNdim] = 3100;
1763
1764     st = "tracklet spatial&charge resolution;";
1765     // define minimum info to be saved in non debug mode
1766     Int_t ndim=DebugLevel()>=1?(kNdim+1):4;
1767     for(Int_t idim(0); idim<ndim; idim++){ st += trTitle[idim]; st+=";";}
1768     H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax);
1769   } else H->Reset();
1770   fContainer->AddAt(H, kTracklet);
1771   //++++++++++++++++++++++
1772   // tracklet to TRDin
1773   snprintf(hn, nhn, "h%s", fgPerformanceName[kTrackIn]);
1774   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1775     st = "r-#phi/z/angular residuals @ TRD entry;";
1776     // define minimum info to be saved in non debug mode
1777     Int_t ndim=DebugLevel()>=1?kNdim:7;
1778     for(Int_t idim(0); idim<ndim; idim++){ st += fgkTitle[idim]; st+=";";}
1779     H = new THnSparseI(hn, st.Data(), ndim, fgkNbins, fgkMin, fgkMax);
1780   } else H->Reset();
1781   fContainer->AddAt(H, kTrackIn);
1782   // tracklet to TRDout
1783 //  fContainer->AddAt(BuildMonitorContainerTracklet("TrkOUT"), kTrackOut);
1784
1785
1786   // Resolution histos
1787   if(!HasMCdata()) return fContainer;
1788
1789   //++++++++++++++++++++++
1790   // cluster to TrackRef residuals/pulls
1791   snprintf(hn, nhn, "h%s", fgPerformanceName[kMCcluster]);
1792   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1793     const Char_t *clTitle[kNdim] = {"layer", fgkTitle[kPhi], fgkTitle[kEta], fgkTitle[kYrez], "#Deltax [cm]", "Q</Q", fgkTitle[kSpeciesChgRC], "#Phi [deg]"};
1794     const Int_t clNbins[kNdim]   = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], fgkNbins[kEta], fgkNbins[kYrez], 20, 10, fgkNbins[kSpeciesChgRC], 15};
1795     const Double_t clMin[kNdim]  = {-0.5, fgkMin[kPhi], fgkMin[kEta], fgkMin[kYrez]/10., 0., 0.1, fgkMin[kSpeciesChgRC], -45},
1796                    clMax[kNdim]  = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], fgkMax[kEta], fgkMax[kYrez]/10., 4., 2.1, fgkMax[kSpeciesChgRC], 45};
1797     st = "MC cluster spatial resolution;";
1798     // define minimum info to be saved in non debug mode
1799     Int_t ndim=DebugLevel()>=1?kNdim:4;
1800     for(Int_t idim(0); idim<ndim; idim++){ st += clTitle[idim]; st+=";";}
1801     H = new THnSparseI(hn, st.Data(), ndim, clNbins, clMin, clMax);
1802   } else H->Reset();
1803   fContainer->AddAt(H, kMCcluster);
1804   //++++++++++++++++++++++
1805   // tracklet to TrackRef
1806   snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtracklet]);
1807   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1808     Char_t *trTitle[kNdim]; memcpy(trTitle, fgkTitle, kNdim*sizeof(Char_t*));
1809     Int_t trNbins[kNdim]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t));
1810     Double_t trMin[kNdim]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
1811     Double_t trMax[kNdim]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
1812     // set specific fields
1813     trTitle[kBC]=StrDup("layer"); trNbins[kBC] = AliTRDgeometry::kNlayer; trMin[kBC] = -0.5; trMax[kBC] = AliTRDgeometry::kNlayer-0.5;
1814     trMin[kYrez] = -0.54; trMax[kYrez] = -trMin[kYrez];
1815     trMin[kPrez] = -4.5; trMax[kPrez] = -trMin[kPrez];
1816     trMin[kZrez] = -1.5; trMax[kZrez] = -trMin[kZrez];
1817
1818     st = "MC tracklet spatial resolution;";
1819     // define minimum info to be saved in non debug mode
1820     Int_t ndim=DebugLevel()>=1?kNdim:4;
1821     for(Int_t idim(0); idim<ndim; idim++){ st += trTitle[idim]; st+=";";}
1822     H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax);
1823   } else H->Reset();
1824   fContainer->AddAt(H, kMCtracklet);
1825   //++++++++++++++++++++++
1826   // TRDin to TrackRef
1827   snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtrackIn]);
1828   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1829     st = "MC r-#phi/z/angular residuals @ TRD entry;";
1830     // set specific fields
1831     Double_t trMin[kNdim]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
1832     Double_t trMax[kNdim]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
1833     trMin[kYrez] = -0.54; trMax[kYrez] = -trMin[kYrez];
1834     trMin[kPrez] = -2.4; trMax[kPrez] = -trMin[kPrez];
1835     trMin[kZrez] = -0.9; trMax[kZrez] = -trMin[kZrez];
1836     // define minimum info to be saved in non debug mode
1837     Int_t ndim=DebugLevel()>=1?kNdim:7;
1838     for(Int_t idim(0); idim<ndim; idim++){ st += fgkTitle[idim]; st+=";";}
1839     H = new THnSparseI(hn, st.Data(), ndim, fgkNbins, trMin, trMax);
1840   } else H->Reset();
1841   fContainer->AddAt(H, kMCtrackIn);
1842   //++++++++++++++++++++++
1843   // track to TrackRef
1844   snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtrack]);
1845   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1846     Char_t *trTitle[kNdim+1]; memcpy(trTitle, fgkTitle, kNdim*sizeof(Char_t*));
1847     Int_t trNbins[kNdim+1]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t));
1848     Double_t trMin[kNdim+1]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
1849     Double_t trMax[kNdim+1]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
1850     // set specific fields
1851     trTitle[kBC]=StrDup("layer"); trNbins[kBC] = AliTRDgeometry::kNlayer; trMin[kBC] = -0.5; trMax[kBC] = AliTRDgeometry::kNlayer-0.5;
1852     trTitle[kNdim]=StrDup("#Deltap_{t}/p_{t} [%]"); trNbins[kNdim] = 30; trMin[kNdim] = -15.; trMax[kNdim] = 15.;
1853     trMin[kYrez] = -0.9; trMax[kYrez] = -trMin[kYrez];
1854     trMin[kPrez] = -1.5; trMax[kPrez] = -trMin[kPrez];
1855     trMin[kZrez] = -0.9; trMax[kZrez] = -trMin[kZrez];
1856
1857     st = "MC track spatial&p_{t} resolution;";
1858     // define minimum info to be saved in non debug mode
1859     Int_t ndim=DebugLevel()>=1?(kNdim+1):4;
1860     for(Int_t idim(0); idim<ndim; idim++){ st += trTitle[idim]; st+=";";}
1861     H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax);
1862   } else H->Reset();
1863   fContainer->AddAt(H, kMCtrack);
1864
1865 //   // cluster resolution
1866 //   fContainer->AddAt(BuildMonitorContainerCluster("MCcl"),  kMCcluster);
1867 //   // track resolution
1868 //   TObjArray *arr(NULL);
1869 //   fContainer->AddAt(arr = new TObjArray(AliTRDgeometry::kNlayer), kMCtrack);
1870 //   arr->SetName("MCtrk");
1871 //   for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++) arr->AddAt(BuildMonitorContainerTrack(Form("MCtrk_Ly%d", il)), il);
1872 //   // TRDin TRACK RESOLUTION
1873 //   fContainer->AddAt(H, kMCtrackIn);
1874 //   // TRDout TRACK RESOLUTION
1875 //   fContainer->AddAt(BuildMonitorContainerTrack("MCtrkOUT"), kMCtrackOut);
1876
1877   return fContainer;
1878 }
1879
1880 //________________________________________________________
1881 Bool_t AliTRDresolution::Process(TH2* const h2, TGraphErrors **g, Int_t stat)
1882 {
1883 // Robust function to process sigma/mean for 2D plot dy(x)
1884 // For each x bin a gauss fit is performed on the y projection and the range
1885 // with the minimum chi2/ndf is choosen
1886
1887   if(!h2) {
1888     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : NULL pointer input container.\n");
1889     return kFALSE;
1890   }
1891   if(!Int_t(h2->GetEntries())){
1892     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : Empty h[%s - %s].\n", h2->GetName(), h2->GetTitle());
1893     return kFALSE;
1894   }
1895   if(!g || !g[0]|| !g[1]) {
1896     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : NULL pointer output container.\n");
1897     return kFALSE;
1898   }
1899   // prepare
1900   TAxis *ax(h2->GetXaxis()), *ay(h2->GetYaxis());
1901   Float_t ymin(ay->GetXmin()), ymax(ay->GetXmax()), dy(ay->GetBinWidth(1)), y0(0.), y1(0.);
1902   TF1 f("f", "gaus", ymin, ymax);
1903   Int_t n(0);
1904   if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
1905   if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
1906   TH1D *h(NULL);
1907   if((h=(TH1D*)gROOT->FindObject("py"))) delete h;
1908   Double_t x(0.), y(0.), ex(0.), ey(0.), sy(0.), esy(0.);
1909
1910
1911   // do actual loop
1912   Float_t chi2OverNdf(0.);
1913   for(Int_t ix = 1, np=0; ix<=ax->GetNbins(); ix++){
1914     x = ax->GetBinCenter(ix); ex= ax->GetBinWidth(ix)*0.288; // w/sqrt(12)
1915     ymin = ay->GetXmin(); ymax = ay->GetXmax();
1916
1917     h = h2->ProjectionY("py", ix, ix);
1918     if((n=(Int_t)h->GetEntries())<stat){
1919       if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("I-AliTRDresolution::Process() : Low statistics @ x[%f] stat[%d]=%d [%d].\n", x, ix, n, stat);
1920       continue;
1921     }
1922     // looking for a first order mean value
1923     f.SetParameter(1, 0.); f.SetParameter(2, 3.e-2);
1924     h->Fit(&f, "QNW");
1925     chi2OverNdf = f.GetChisquare()/f.GetNDF();
1926     printf("x[%f] range[%f %f] chi2[%f] ndf[%d] chi2/ndf[%f]\n", x, ymin, ymax, f.GetChisquare(),f.GetNDF(),chi2OverNdf);
1927     y = f.GetParameter(1); y0 = y-4*dy; y1 = y+4*dy;
1928     ey  = f.GetParError(1);
1929     sy = f.GetParameter(2); esy = f.GetParError(2);
1930 //     // looking for the best chi2/ndf value
1931 //     while(ymin<y0 && ymax>y1){
1932 //       f.SetParameter(1, y);
1933 //       f.SetParameter(2, sy);
1934 //       h->Fit(&f, "QNW", "", y0, y1);
1935 //       printf("   range[%f %f] chi2[%f] ndf[%d] chi2/ndf[%f]\n", y0, y1, f.GetChisquare(),f.GetNDF(),f.GetChisquare()/f.GetNDF());
1936 //       if(f.GetChisquare()/f.GetNDF() < Chi2OverNdf){
1937 //         chi2OverNdf = f.GetChisquare()/f.GetNDF();
1938 //         y  = f.GetParameter(1); ey  = f.GetParError(1);
1939 //         sy = f.GetParameter(2); esy = f.GetParError(2);
1940 //         printf("    set y[%f] sy[%f] chi2/ndf[%f]\n", y, sy, chi2OverNdf);
1941 //       }
1942 //       y0-=dy; y1+=dy;
1943 //     }
1944     g[0]->SetPoint(np, x, y);
1945     g[0]->SetPointError(np, ex, ey);
1946     g[1]->SetPoint(np, x, sy);
1947     g[1]->SetPointError(np, ex, esy);
1948     np++;
1949   }
1950   return kTRUE;
1951 }
1952
1953
1954 //________________________________________________________
1955 Bool_t AliTRDresolution::Process(TH2 * const h2, TF1 *f, Float_t k, TGraphErrors **g)
1956 {
1957   //
1958   // Do the processing
1959   //
1960
1961   Char_t pn[10]; snprintf(pn, 10, "p%03d", fIdxPlot);
1962   Int_t n = 0;
1963   if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
1964   if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
1965   if(Int_t(h2->GetEntries())){
1966     AliDebug(4, Form("%s: g[%s %s]", pn, g[0]->GetName(), g[0]->GetTitle()));
1967   } else {
1968     AliDebug(2, Form("%s: g[%s %s]: Missing entries.", pn, g[0]->GetName(), g[0]->GetTitle()));
1969     fIdxPlot++;
1970     return kTRUE;
1971   }
1972
1973   const Int_t kINTEGRAL=1;
1974   for(Int_t ibin = 0; ibin < Int_t(h2->GetNbinsX()/kINTEGRAL); ibin++){
1975     Int_t abin(ibin*kINTEGRAL+1),bbin(abin+kINTEGRAL-1),mbin(abin+Int_t(kINTEGRAL/2));
1976     Double_t x = h2->GetXaxis()->GetBinCenter(mbin);
1977     TH1D *h = h2->ProjectionY(pn, abin, bbin);
1978     if((n=(Int_t)h->GetEntries())<150){
1979       AliDebug(4, Form("  x[%f] range[%d %d] stat[%d] low statistics !", x, abin, bbin, n));
1980       continue;
1981     }
1982     h->Fit(f, "QN");
1983     Int_t ip = g[0]->GetN();
1984     AliDebug(4, Form("  x_%d[%f] range[%d %d] stat[%d] M[%f] Sgm[%f]", ip, x, abin, bbin, n, f->GetParameter(1), f->GetParameter(2)));
1985     g[0]->SetPoint(ip, x, k*f->GetParameter(1));
1986     g[0]->SetPointError(ip, 0., k*f->GetParError(1));
1987     g[1]->SetPoint(ip, x, k*f->GetParameter(2));
1988     g[1]->SetPointError(ip, 0., k*f->GetParError(2));
1989 /*
1990     g[0]->SetPoint(ip, x, k*h->GetMean());
1991     g[0]->SetPointError(ip, 0., k*h->GetMeanError());
1992     g[1]->SetPoint(ip, x, k*h->GetRMS());
1993     g[1]->SetPointError(ip, 0., k*h->GetRMSError());*/
1994   }
1995   fIdxPlot++;
1996   return kTRUE;
1997 }
1998
1999
2000 //____________________________________________________________________
2001 Bool_t AliTRDresolution::FitTrack(const Int_t np, AliTrackPoint *points, Float_t param[10])
2002 {
2003 //
2004 // Fit track with a staight line using the "np" clusters stored in the array "points".
2005 // The following particularities are stored in the clusters from points:
2006 //   1. pad tilt as cluster charge
2007 //   2. pad row cross or vertex constrain as fake cluster with cluster type 1
2008 // The parameters of the straight line fit are stored in the array "param" in the following order :
2009 //     param[0] - x0 reference radial position
2010 //     param[1] - y0 reference r-phi position @ x0
2011 //     param[2] - z0 reference z position @ x0
2012 //     param[3] - slope dy/dx
2013 //     param[4] - slope dz/dx
2014 //
2015 // Attention :
2016 // Function should be used to refit tracks for B=0T
2017 //
2018
2019   if(np<40){
2020     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTrack: Not enough clusters to fit a track [%d].\n", np);
2021     return kFALSE;
2022   }
2023   TLinearFitter yfitter(2, "pol1"), zfitter(2, "pol1");
2024
2025   Double_t x0(0.);
2026   for(Int_t ip(0); ip<np; ip++) x0+=points[ip].GetX();
2027   x0/=Float_t(np);
2028
2029   Double_t x, y, z, dx, tilt(0.);
2030   for(Int_t ip(0); ip<np; ip++){
2031     x = points[ip].GetX(); z = points[ip].GetZ();
2032     dx = x - x0;
2033     zfitter.AddPoint(&dx, z, points[ip].GetClusterType()?1.e-3:1.);
2034   }
2035   if(zfitter.Eval() != 0) return kFALSE;
2036
2037   Double_t z0    = zfitter.GetParameter(0);
2038   Double_t dzdx  = zfitter.GetParameter(1);
2039   for(Int_t ip(0); ip<np; ip++){
2040     if(points[ip].GetClusterType()) continue;
2041     x    = points[ip].GetX();
2042     dx   = x - x0;
2043     y    = points[ip].GetY();
2044     z    = points[ip].GetZ();
2045     tilt = points[ip].GetCharge();
2046     y -= tilt*(-dzdx*dx + z - z0);
2047     Float_t xyz[3] = {x, y, z}; points[ip].SetXYZ(xyz);
2048     yfitter.AddPoint(&dx, y, 1.);
2049   }
2050   if(yfitter.Eval() != 0) return kFALSE;
2051   Double_t y0   = yfitter.GetParameter(0);
2052   Double_t dydx = yfitter.GetParameter(1);
2053
2054   param[0] = x0; param[1] = y0; param[2] = z0; param[3] = dydx; param[4] = dzdx;
2055   if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>3) printf("D-AliTRDresolution::FitTrack: x0[%f] y0[%f] z0[%f] dydx[%f] dzdx[%f].\n", x0, y0, z0, dydx, dzdx);
2056   return kTRUE;
2057 }
2058
2059 //____________________________________________________________________
2060 Bool_t AliTRDresolution::FitTracklet(const Int_t ly, const Int_t np, const AliTrackPoint *points, const Float_t param[10], Float_t par[3])
2061 {
2062 //
2063 // Fit tracklet with a staight line using the coresponding subset of clusters out of the total "np" clusters stored in the array "points".
2064 // See function FitTrack for the data stored in the "clusters" array
2065
2066 // The parameters of the straight line fit are stored in the array "param" in the following order :
2067 //     par[0] - x0 reference radial position
2068 //     par[1] - y0 reference r-phi position @ x0
2069 //     par[2] - slope dy/dx
2070 //
2071 // Attention :
2072 // Function should be used to refit tracks for B=0T
2073 //
2074
2075   TLinearFitter yfitter(2, "pol1");
2076
2077   // grep data for tracklet
2078   Double_t x0(0.), x[60], y[60], dy[60];
2079   Int_t nly(0);
2080   for(Int_t ip(0); ip<np; ip++){
2081     if(points[ip].GetClusterType()) continue;
2082     if(points[ip].GetVolumeID() != ly) continue;
2083     Float_t xt(points[ip].GetX())
2084            ,yt(param[1] + param[3] * (xt - param[0]));
2085     x[nly] = xt;
2086     y[nly] = points[ip].GetY();
2087     dy[nly]= y[nly]-yt;
2088     x0    += xt;
2089     nly++;
2090   }
2091   if(nly<10){
2092     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTracklet: Not enough clusters to fit a tracklet [%d].\n", nly);
2093     return kFALSE;
2094   }
2095   // set radial reference for fit
2096   x0 /= Float_t(nly);
2097
2098   // find tracklet core
2099   Double_t mean(0.), sig(1.e3);
2100   AliMathBase::EvaluateUni(nly, dy, mean, sig, 0);
2101
2102   // simple cluster error parameterization
2103   Float_t kSigCut = TMath::Sqrt(5.e-4 + param[3]*param[3]*0.018);
2104
2105   // fit tracklet core
2106   for(Int_t jly(0); jly<nly; jly++){
2107     if(TMath::Abs(dy[jly]-mean)>kSigCut) continue;
2108     Double_t dx(x[jly]-x0);
2109     yfitter.AddPoint(&dx, y[jly], 1.);
2110   }
2111   if(yfitter.Eval() != 0) return kFALSE;
2112   par[0] = x0;
2113   par[1] = yfitter.GetParameter(0);
2114   par[2] = yfitter.GetParameter(1);
2115   return kTRUE;
2116 }
2117
2118 //____________________________________________________________________
2119 Bool_t AliTRDresolution::UseTrack(const Int_t np, const AliTrackPoint *points, Float_t param[10])
2120 {
2121 //
2122 // Global selection mechanism of tracksbased on cluster to fit residuals
2123 // The parameters are the same as used ni function FitTrack().
2124
2125   const Float_t kS(0.6), kM(0.2);
2126   TH1S h("h1", "", 100, -5.*kS, 5.*kS);
2127   Float_t dy, dz, s, m;
2128   for(Int_t ip(0); ip<np; ip++){
2129     if(points[ip].GetClusterType()) continue;
2130     Float_t x0(points[ip].GetX())
2131           ,y0(param[1] + param[3] * (x0 - param[0]))
2132           ,z0(param[2] + param[4] * (x0 - param[0]));
2133     dy=points[ip].GetY() - y0; h.Fill(dy);
2134     dz=points[ip].GetZ() - z0;
2135   }
2136   TF1 fg("fg", "gaus", -5.*kS, 5.*kS);
2137   fg.SetParameter(1, 0.);
2138   fg.SetParameter(2, 2.e-2);
2139   h.Fit(&fg, "QN");
2140   m=fg.GetParameter(1); s=fg.GetParameter(2);
2141   if(s>kS || TMath::Abs(m)>kM) return kFALSE;
2142   return kTRUE;
2143 }
2144
2145 //________________________________________________________
2146 void AliTRDresolution::GetLandauMpvFwhm(TF1 * const f, Float_t &mpv, Float_t &xm, Float_t &xM)
2147 {
2148   //
2149   // Get the most probable value and the full width half mean
2150   // of a Landau distribution
2151   //
2152
2153   const Float_t dx = 1.;
2154   mpv = f->GetParameter(1);
2155   Float_t fx, max = f->Eval(mpv);
2156
2157   xm = mpv - dx;
2158   while((fx = f->Eval(xm))>.5*max){
2159     if(fx>max){
2160       max = fx;
2161       mpv = xm;
2162     }
2163     xm -= dx;
2164   }
2165
2166   xM += 2*(mpv - xm);
2167   while((fx = f->Eval(xM))>.5*max) xM += dx;
2168 }
2169
2170
2171 // #include "TFile.h"
2172 // //________________________________________________________
2173 // Bool_t AliTRDresolution::LoadCorrection(const Char_t *file)
2174 // {
2175 //   if(!file){
2176 //     AliWarning("Use cluster position as in reconstruction.");
2177 //     SetLoadCorrection();
2178 //     return kTRUE;
2179 //   }
2180 //   TDirectory *cwd(gDirectory);
2181 //   TString fileList;
2182 //   FILE *filePtr = fopen(file, "rt");
2183 //   if(!filePtr){
2184 //     AliWarning(Form("Couldn't open correction list \"%s\". Use cluster position as in reconstruction.", file));
2185 //     SetLoadCorrection();
2186 //     return kFALSE;
2187 //   }
2188 //   TH2 *h2 = new TH2F("h2", ";time [time bins];detector;dx [#mum]", 30, -0.5, 29.5, AliTRDgeometry::kNdet, -0.5, AliTRDgeometry::kNdet-0.5);
2189 //   while(fileList.Gets(filePtr)){
2190 //     if(!TFile::Open(fileList.Data())) {
2191 //       AliWarning(Form("Couldn't open \"%s\"", fileList.Data()));
2192 //       continue;
2193 //     } else AliInfo(Form("\"%s\"", fileList.Data()));
2194 //
2195 //     TTree *tSys = (TTree*)gFile->Get("tSys");
2196 //     h2->SetDirectory(gDirectory); h2->Reset("ICE");
2197 //     tSys->Draw("det:t>>h2", "dx", "goff");
2198 //     for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
2199 //       for(Int_t it(0); it<30; it++) fXcorr[idet][it]+=(1.e-4*h2->GetBinContent(it+1, idet+1));
2200 //     }
2201 //     h2->SetDirectory(cwd);
2202 //     gFile->Close();
2203 //   }
2204 //   cwd->cd();
2205 //
2206 //   if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>=2){
2207 //     for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
2208 //     printf("DET|");for(Int_t it(0); it<30; it++) printf(" tb%02d|", it); printf("\n");
2209 //     for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
2210 //     FILE *fout = fopen("TRD.ClusterCorrection.txt", "at");
2211 //     fprintf(fout, "  static const Double_t dx[AliTRDgeometry::kNdet][30]={\n");
2212 //     for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
2213 //       printf("%03d|", idet);
2214 //       fprintf(fout, "    {");
2215 //       for(Int_t it(0); it<30; it++){
2216 //         printf("%+5.0f|", 1.e4*fXcorr[idet][it]);
2217 //         fprintf(fout, "%+6.4f%c", fXcorr[idet][it], it==29?' ':',');
2218 //       }
2219 //       printf("\n");
2220 //       fprintf(fout, "}%c\n", idet==AliTRDgeometry::kNdet-1?' ':',');
2221 //     }
2222 //     fprintf(fout, "  };\n");
2223 //   }
2224 //   SetLoadCorrection();
2225 //   return kTRUE;
2226 // }
2227
2228 //________________________________________________________
2229 AliTRDresolution::AliTRDresolutionProjection::AliTRDresolutionProjection()
2230   :fH(NULL)
2231   ,fNrebin(0)
2232   ,fRebinX(NULL)
2233   ,fRebinY(NULL)
2234 {
2235   // constructor
2236   memset(fAx, 0, 3*sizeof(Int_t));
2237   memset(fRange, 0, 4*sizeof(Float_t));
2238 }
2239
2240 //________________________________________________________
2241 AliTRDresolution::AliTRDresolutionProjection::~AliTRDresolutionProjection()
2242 {
2243   // destructor
2244   if(fH) delete fH;
2245 }
2246
2247 //________________________________________________________
2248 void AliTRDresolution::AliTRDresolutionProjection::Build(const Char_t *n, const Char_t *t, Int_t ix, Int_t iy, Int_t iz, TAxis *aa[])
2249 {
2250 // check and build (if neccessary) projection determined by axis "ix", "iy" and "iz"
2251   if(!aa[ix] || !aa[iy] || !aa[iz]) return;
2252   TAxis *ax(aa[ix]), *ay(aa[iy]), *az(aa[iz]);
2253   fH = new TH3I(n, Form("%s;%s;%s;%s", t, ax->GetTitle(), ay->GetTitle(), az->GetTitle()),
2254     ax->GetNbins(), ax->GetXmin(), ax->GetXmax(),
2255     ay->GetNbins(), ay->GetXmin(), ay->GetXmax(),
2256     az->GetNbins(), az->GetXmin(), az->GetXmax());
2257   fAx[0] = ix; fAx[1] = iy; fAx[2] = iz;
2258   fRange[0] = az->GetXmin()/3.; fRange[1] = az->GetXmax()/3.;
2259 }
2260
2261 //________________________________________________________
2262 void AliTRDresolution::AliTRDresolutionProjection::Increment(Int_t bin[], Double_t v)
2263 {
2264 // increment bin with value "v" pointed by general coord in "bin"
2265   if(!fH) return;
2266   fH->AddBinContent(
2267         fH->GetBin(bin[fAx[0]],bin[fAx[1]],bin[fAx[2]]), v);
2268 }
2269
2270 //________________________________________________________
2271 TH2* AliTRDresolution::AliTRDresolutionProjection::Projection2D(const Int_t nstat, const Int_t ncol, const Int_t mid)
2272 {
2273 // build the 2D projection and adjust binning
2274
2275   const Char_t *title[] = {"Mean", "#mu", "MPV"};
2276   if(!fH) return NULL;
2277   TAxis *ax(fH->GetXaxis()), *ay(fH->GetYaxis()), *az(fH->GetZaxis());
2278   TH2 *h2s = (TH2*)fH->Project3D("yx");
2279   Int_t irebin(0), dxBin(1), dyBin(1);
2280   while(irebin<fNrebin && (AliTRDresolution::GetMeanStat(h2s, .5, ">")<nstat)){
2281     h2s->Rebin2D(fRebinX[irebin], fRebinY[irebin]);
2282     dxBin*=fRebinX[irebin];dyBin*=fRebinY[irebin];
2283     irebin++;
2284   }
2285   Int_t nx(h2s->GetNbinsX()), ny(h2s->GetNbinsY());
2286   if(h2s) delete h2s;
2287
2288   // start projection
2289   TH1 *h(NULL);
2290   Float_t dz=(fRange[1]-fRange[1])/ncol;
2291   TString titlez(az->GetTitle()); TObjArray *tokenTitle(titlez.Tokenize(" "));
2292   TH2 *h2 = new TH2F(Form("%s_2D", fH->GetName()),
2293             Form("%s;%s;%s;%s(%s) %s", fH->GetTitle(), ax->GetTitle(), ay->GetTitle(), title[mid], (*tokenTitle)[0]->GetName(), tokenTitle->GetEntriesFast()>1?(*tokenTitle)[1]->GetName():""),
2294             nx, ax->GetXmin(), ax->GetXmax(), ny, ay->GetXmin(), ay->GetXmax());
2295   h2->SetContour(ncol);
2296   h2->GetZaxis()->CenterTitle();
2297   h2->GetZaxis()->SetRangeUser(fRange[0], fRange[1]);
2298   //printf("%s[%s] nx[%d] ny[%d]\n", h2->GetName(), h2->GetTitle(), nx, ny);
2299   for(Int_t iy(0); iy<ny; iy++){
2300     for(Int_t ix(0); ix<nx; ix++){
2301       h = fH->ProjectionZ(Form("%s_z", h2->GetName()), ix*dxBin+1, (ix+1)*dxBin+1, iy*dyBin+1, (iy+1)*dyBin+1);
2302       Int_t ne((Int_t)h->Integral());
2303       if(ne<nstat/2){
2304         h2->SetBinContent(ix+1, iy+1, -999);
2305         h2->SetBinError(ix+1, iy+1, 1.);
2306       }else{
2307         Float_t v(h->GetMean()), ve(h->GetRMS());
2308         if(mid==1){
2309           TF1 fg("fg", "gaus", az->GetXmin(), az->GetXmax());
2310           fg.SetParameter(0, Float_t(ne)); fg.SetParameter(1, v); fg.SetParameter(2, ve);
2311           h->Fit(&fg, "WQ");
2312           v = fg.GetParameter(1); ve = fg.GetParameter(2);
2313         } else if (mid==2) {
2314           TF1 fl("fl", "landau", az->GetXmin(), az->GetXmax());
2315           fl.SetParameter(0, Float_t(ne)); fl.SetParameter(1, v); fl.SetParameter(2, ve);
2316           h->Fit(&fl, "WQ");
2317           v = fl.GetMaximumX(); ve = fl.GetParameter(2);
2318 /*          TF1 fgle("gle", "[0]*TMath::Landau(x, [1], [2], 1)*TMath::Exp(-[3]*x/[1])", az->GetXmin(), az->GetXmax());
2319           fgle.SetParameter(0, fl.GetParameter(0));
2320           fgle.SetParameter(1, fl.GetParameter(1));
2321           fgle.SetParameter(2, fl.GetParameter(2));
2322           fgle.SetParameter(3, 1.);fgle.SetParLimits(3, 0., 5.);
2323           h->Fit(&fgle, "WQ");
2324           v = fgle.GetMaximumX(); ve = fgle.GetParameter(2);*/
2325         }
2326         if(v<fRange[0]) h2->SetBinContent(ix+1, iy+1, fRange[0]+0.1*dz);
2327         else h2->SetBinContent(ix+1, iy+1, v);
2328         h2->SetBinError(ix+1, iy+1, ve);
2329       }
2330     }
2331   }
2332   if(h) delete h;
2333   return h2;
2334 }
2335
2336 void AliTRDresolution::AliTRDresolutionProjection::SetRebinStrategy(Int_t n, Int_t rebx[], Int_t reby[])
2337 {
2338 // define rebinning strategy for this projection
2339   fNrebin = n;
2340   fRebinX = new Int_t[n]; memcpy(fRebinX, rebx, n*sizeof(Int_t));
2341   fRebinY = new Int_t[n]; memcpy(fRebinY, reby, n*sizeof(Int_t));
2342 }
2343
2344