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