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