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