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