PWG0/dNdPt/AlidNdPtAnalysisPbPb.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-commercial 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 // AlidNdPtAnalysisPbPb class.
  17 //
  18 // a. functionality:
  19 // - fills analysis control histograms
  20 // - fills generic correction matrices
  21 // - generates correction matrices
  22 //
  23 // b. data members:
  24 // - generic correction matrices
  25 // - control histograms
  26 //
  27 // Author: J.Otwinowski 04/11/2008
  28 // last change: 2011-04-04 by M.Knichel
  29 //------------------------------------------------------------------------------
  30
  31 #include "TH1.h"
  32 #include "TH2.h"
  33 #include "TCanvas.h"
  34 #include "THnSparse.h"
  35
  36 #include "AliHeader.h"
  37 #include "AliGenEventHeader.h"
  38 #include "AliInputEventHandler.h"
  39 #include "AliAnalysisManager.h"
  40 #include "AliStack.h"
  41 #include "AliESDEvent.h"
  42 #include "AliMCEvent.h"
  43 #include "AliESDtrackCuts.h"
  44 #include "AliLog.h"
  45 #include "AliMultiplicity.h"
  46 #include "AliTracker.h"
  47
  48 #include "AliCentrality.h"
  49
  50 #include "AlidNdPtEventCuts.h"
  51 #include "AlidNdPtAcceptanceCuts.h"
  52 #include "AliPhysicsSelection.h"
  53 #include "AliTriggerAnalysis.h"
  54
  55 #include "AliPWG0Helper.h"
  56 #include "AlidNdPtHelper.h"
  57 #include "AlidNdPtAnalysisPbPb.h"
  58
  59
  60 using namespace std;
  61
  62 ClassImp(AlidNdPtAnalysisPbPb)
  63
  64 //_____________________________________________________________________________
  65   AlidNdPtAnalysisPbPb::AlidNdPtAnalysisPbPb(): AlidNdPt(),
  66   fAnalysisFolder(0),
  67   fHistogramsOn(kFALSE),
  68
  69   // rec. track pt vs true track pt correlation matrix
  70   fTrackPtCorrelationMatrix(0),
  71
  72   // event level correction
  73   fGenEventMatrix(0),
  74   fTriggerEventMatrix(0),
  75   fRecEventMatrix(0),
  76
  77   //
  78   // track-event level correction
  79   //
  80   fGenTrackEventMatrix(0),
  81
  82   fTriggerTrackEventMatrix(0),
  83
  84   fRecTrackEventMatrix(0),
  85
  86   // track rec. efficiency correction (fRecPrimTrackMatrix / fGenPrimTrackMatrix)
  87   fGenTrackMatrix(0),
  88   fGenPrimTrackMatrix(0),
  89   fRecPrimTrackMatrix(0),
  90
  91   // secondary track contamination correction (fRecSecTrackMatrix / fRecTrackMatrix)
  92   fRecTrackMatrix(0),
  93   fRecSecTrackMatrix(0),
  94
  95   // multiple rec. track contamination corrections (fRecMultTrackMatrix / fRecTrackMatrix)
  96   fRecMultTrackMatrix(0),
  97
  98   // event control histograms
  99   fMCEventHist1(0),
 100   fRecEventHist1(0),
 101   fRecEventHist2(0),
 102   fRecMCEventHist1(0),
 103   fRecMCEventHist2(0),
 104
 105   // rec. pt and eta resolution w.r.t MC
 106   fRecMCTrackHist1(0),
 107
 108   //multple reconstructed tracks
 109   fMCMultRecTrackHist1(0),
 110
 111   // rec. track control histograms
 112   fRecTrackHist3(0),
 113
 114   fTriggerAnalysis(0),
 115
 116   fCentralityEstimator(0),
 117
 118   fMultNbins(0),
 119   fPtNbins(0),
 120   fPtCorrNbins(0),
 121   fEtaNbins(0),
 122   fZvNbins(0),
 123   fCentralityNbins(0),
 124   fBinsMult(0),
 125   fBinsPt(0),
 126   fBinsPtCorr(0),
 127   fBinsEta(0),
 128   fBinsZv(0),
 129   fBinsCentrality(0),
 130
 131   fIsInit(kFALSE)
 132 {
 133   // default constructor
 134   for(Int_t i=0; i<AlidNdPtHelper::kCutSteps; i++) {
 135     fMCTrackHist1[i]=0;
 136     fMCPrimTrackHist1[i]=0;
 137     fMCPrimTrackHist2[i]=0;
 138     fMCSecTrackHist1[i]=0;
 139     fRecTrackHist1[i]=0;
 140     fRecTrackHist2[i]=0;
 141     fRecTrackMultHist1[i]=0;
 142   }
 143   //Init();
 144   SetCentralityEstimator();
 145 }
 146
 147 //_____________________________________________________________________________
 148 AlidNdPtAnalysisPbPb::AlidNdPtAnalysisPbPb(Char_t* name, Char_t* title): AlidNdPt(name,title),
 149   fAnalysisFolder(0),
 150   fHistogramsOn(kFALSE),
 151
 152   // rec. track pt vs true track pt correlation matrix
 153   fTrackPtCorrelationMatrix(0),
 154
 155   // event level correction
 156   fGenEventMatrix(0),
 157   fTriggerEventMatrix(0),
 158   fRecEventMatrix(0),
 159
 160   //
 161   // track-event level correction
 162   //
 163   fGenTrackEventMatrix(0),
 164
 165   fTriggerTrackEventMatrix(0),
 166
 167   fRecTrackEventMatrix(0),
 168
 169   // track rec. efficiency correction (fRecPrimTrackMatrix / fGenPrimTrackMatrix)
 170   fGenTrackMatrix(0),
 171   fGenPrimTrackMatrix(0),
 172   fRecPrimTrackMatrix(0),
 173
 174   // secondary track contamination correction (fRecTrackMatrix - fRecSecTrackMatrix)
 175   fRecTrackMatrix(0),
 176   fRecSecTrackMatrix(0),
 177
 178   // multiple rec. track contamination corrections (fRecTrackMatrix - fRecMultTrackMatrix)
 179   fRecMultTrackMatrix(0),
 180
 181
 182
 183
 184   // event control histograms
 185   fMCEventHist1(0),
 186   fRecEventHist1(0),
 187   fRecEventHist2(0),
 188   fRecMCEventHist1(0),
 189   fRecMCEventHist2(0),
 190
 191   // rec. pt and eta resolution w.r.t MC
 192   fRecMCTrackHist1(0),
 193
 194   //multiple reconstructed tracks
 195   fMCMultRecTrackHist1(0),
 196
 197   // rec. track control histograms
 198   fRecTrackHist3(0),
 199
 200   fTriggerAnalysis(0),
 201
 202   fCentralityEstimator(0),
 203
 204   fMultNbins(0),
 205   fPtNbins(0),
 206   fPtCorrNbins(0),
 207   fEtaNbins(0),
 208   fZvNbins(0),
 209   fCentralityNbins(0),
 210   fBinsMult(0),
 211   fBinsPt(0),
 212   fBinsPtCorr(0),
 213   fBinsEta(0),
 214   fBinsZv(0),
 215   fBinsCentrality(0),
 216
 217   fIsInit(kFALSE)
 218
 219 {
 220   //
 221   // constructor
 222   //
 223   for(Int_t i=0; i<AlidNdPtHelper::kCutSteps; i++) {
 224     fMCTrackHist1[i]=0;
 225     fMCPrimTrackHist1[i]=0;
 226     fMCPrimTrackHist2[i]=0;
 227     fMCSecTrackHist1[i]=0;
 228     fRecTrackHist1[i]=0;
 229     fRecTrackHist2[i]=0;
 230     fRecTrackMultHist1[i]=0;
 231   }
 232
 233  // Init();
 234  SetCentralityEstimator();
 235 }
 236
 237 //_____________________________________________________________________________
 238 AlidNdPtAnalysisPbPb::~AlidNdPtAnalysisPbPb() {
 239   //
 240   // destructor
 241   //
 242   if(fTrackPtCorrelationMatrix) delete fTrackPtCorrelationMatrix; fTrackPtCorrelationMatrix=0;
 243   //
 244   if(fGenEventMatrix) delete fGenEventMatrix; fGenEventMatrix=0;
 245   if(fTriggerEventMatrix) delete fTriggerEventMatrix; fTriggerEventMatrix=0;
 246   if(fRecEventMatrix) delete fRecEventMatrix; fRecEventMatrix=0;
 247
 248   //
 249   if(fGenTrackEventMatrix) delete fGenTrackEventMatrix; fGenTrackEventMatrix=0;
 250   if(fTriggerEventMatrix) delete fTriggerEventMatrix; fTriggerEventMatrix=0;
 251   if(fRecTrackEventMatrix) delete fRecTrackEventMatrix; fRecTrackEventMatrix=0;
 252
 253   //
 254   if(fGenTrackMatrix) delete fGenTrackMatrix; fGenTrackMatrix=0;
 255   if(fGenPrimTrackMatrix) delete fGenPrimTrackMatrix; fGenPrimTrackMatrix=0;
 256   if(fRecPrimTrackMatrix) delete fRecPrimTrackMatrix; fRecPrimTrackMatrix=0;
 257   //
 258   if(fRecTrackMatrix) delete fRecTrackMatrix; fRecTrackMatrix=0;
 259   if(fRecSecTrackMatrix) delete fRecSecTrackMatrix; fRecSecTrackMatrix=0;
 260   //
 261   if(fRecMultTrackMatrix) delete fRecMultTrackMatrix; fRecMultTrackMatrix=0;
 262   //
 263   // Control histograms
 264   //
 265   if(fMCEventHist1) delete fMCEventHist1; fMCEventHist1=0;
 266   if(fRecEventHist1) delete fRecEventHist1; fRecEventHist1=0;
 267   if(fRecEventHist2) delete fRecEventHist2; fRecEventHist2=0;
 268   if(fRecMCEventHist1) delete fRecMCEventHist1; fRecMCEventHist1=0;
 269   if(fRecMCEventHist2) delete fRecMCEventHist2; fRecMCEventHist2=0;
 270
 271   //
 272   for(Int_t i=0; i<AlidNdPtHelper::kCutSteps; i++) {
 273     if(fMCTrackHist1[i]) delete fMCTrackHist1[i]; fMCTrackHist1[i]=0;
 274     if(fMCPrimTrackHist1[i]) delete fMCPrimTrackHist1[i]; fMCPrimTrackHist1[i]=0;
 275     if(fMCPrimTrackHist2[i]) delete fMCPrimTrackHist2[i]; fMCPrimTrackHist2[i]=0;
 276     if(fMCSecTrackHist1[i]) delete fMCSecTrackHist1[i]; fMCSecTrackHist1[i]=0;
 277     if(fRecTrackHist1[i]) delete fRecTrackHist1[i]; fRecTrackHist1[i]=0;
 278     if(fRecTrackHist2[i]) delete fRecTrackHist2[i]; fRecTrackHist2[i]=0;
 279     if(fRecTrackMultHist1[i]) delete fRecTrackMultHist1[i]; fRecTrackMultHist1[i]=0;
 280   }
 281   if(fRecMCTrackHist1) delete fRecMCTrackHist1; fRecMCTrackHist1=0;
 282   if(fMCMultRecTrackHist1) delete fMCMultRecTrackHist1; fMCMultRecTrackHist1=0;
 283   if(fRecTrackHist3) delete fRecTrackHist3; fRecTrackHist3=0;
 284   //
 285   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
 286
 287   if (fTriggerAnalysis) delete fTriggerAnalysis;  fTriggerAnalysis = 0;
 288 }
 289
 290 //_____________________________________________________________________________
 291 void AlidNdPtAnalysisPbPb::Init() {
 292
 293   Int_t multNbins = 47;
 294   Int_t ptNbinsTrackEventCorr = 37;
 295   Int_t ptNbins = 68;
 296   Int_t etaNbins = 30;
 297   Int_t zvNbins = 12;
 298   Int_t centralityBins = 11;
 299
 300   Double_t binsMultDefault[48] = {
 301                    -0.5, 0.5 , 1.5 , 2.5 , 3.5 , 4.5 , 5.5 , 6.5 , 7.5 , 8.5,
 302                     9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5,
 303                     19.5, 20.5, 30.5, 40.5 , 50.5 , 60.5 , 70.5 , 80.5 , 90.5 , 100.5,
 304                     200.5, 300.5, 400.5, 500.5, 600.5, 700.5, 800.5, 900.5, 1000.5, 2000.5,
 305                     3000.5, 4000.5, 5000.5, 6000.5, 7000.5, 8000.5, 9000.5, 10000.5 }; // forPbPb
 306   Double_t binsPtTrackEventCorrDefault[38] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.4,2.6,3.0,4.0,20.0,50.0};
 307   Double_t binsPtDefault[69] = {
 308         0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45,
 309         0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95,
 310         1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9,
 311         2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8,
 312         4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 8.0, 9.0, 10.0,
 313         11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0, 20.0, 22.0, 24.0,
 314         26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 40.0, 45.0, 50.0 };
 315
 316   Double_t binsEtaDefault[31] = {
 317         -1.5, -1.4, -1.3, -1.2, -1.1, -1.0, -0.9, -0.8, -0.7, -0.6,
 318         -0.5, -0.4, -0.3, -0.2, -0.1,  0.0,  0.1,  0.2,  0.3,  0.4,
 319          0.5,  0.6,  0.7,  0.8,  0.9,  1.0,  1.1,  1.2,  1.3,  1.4,
 320          1.5};
 321
 322   Double_t binsZvDefault[13] = {-30.,-25.,-20.,-15.,-10.,-5.,0.,5.,10.,15.,20.,25.,30.};
 323
 324   Double_t binsCentralityDefault[12] = {0., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.};
 325
 326   Double_t* binsMult = binsMultDefault;
 327   Double_t* binsPtTrackEventCorr = binsPtTrackEventCorrDefault;
 328   Double_t* binsPt = binsPtDefault;
 329   Double_t* binsEta = binsEtaDefault;
 330   Double_t* binsZv = binsZvDefault;
 331   Double_t* binsCentrality = binsCentralityDefault;
 332
 333   if (fMultNbins > 0) { multNbins = fMultNbins; binsMult = fBinsMult; }
 334   if (fPtNbins > 0) { ptNbins  = fPtNbins; binsPt = fBinsPt; }
 335   if (fPtCorrNbins > 0) { ptNbinsTrackEventCorr  = fPtCorrNbins; binsPtTrackEventCorr = fBinsPtCorr;  }
 336   if (fEtaNbins > 0) { etaNbins  = fEtaNbins; binsEta = fBinsEta; }
 337   if (fZvNbins > 0) { zvNbins  = fZvNbins; binsZv = fBinsZv; }
 338   if (fCentralityNbins > 0) { centralityBins  = fCentralityNbins; binsCentrality = fBinsCentrality; }
 339
 340   Int_t binsTrackEventCorrMatrix[4]={zvNbins,ptNbinsTrackEventCorr,etaNbins,centralityBins};
 341   Int_t binsTrackEvent[4]={zvNbins,ptNbins,etaNbins,centralityBins};
 342   Int_t binsTrackPtCorrelationMatrix[4]={ptNbins,ptNbins,etaNbins,centralityBins};
 343
 344
 345   fTrackPtCorrelationMatrix = new THnSparseF("fTrackPtCorrelationMatrix","Pt:mcPt:mcEta:Centrality",4,binsTrackPtCorrelationMatrix);
 346   fTrackPtCorrelationMatrix->SetBinEdges(0,binsPt);
 347   fTrackPtCorrelationMatrix->SetBinEdges(1,binsPt);
 348   fTrackPtCorrelationMatrix->SetBinEdges(2,binsEta);
 349   fTrackPtCorrelationMatrix->SetBinEdges(3,binsCentrality);
 350   fTrackPtCorrelationMatrix->GetAxis(0)->SetTitle("Pt (GeV/c)");
 351   fTrackPtCorrelationMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 352   fTrackPtCorrelationMatrix->GetAxis(2)->SetTitle("mcEta");
 353   fTrackPtCorrelationMatrix->GetAxis(3)->SetTitle("Centrality");
 354   fTrackPtCorrelationMatrix->Sumw2();
 355
 356   //
 357   // Efficiency and contamination correction matrices
 358   //
 359   Int_t binsEventMatrix[3]={zvNbins,multNbins,centralityBins};
 360   Double_t minEventMatrix[3]={-30.,-0.5,0.};
 361   Double_t maxEventMatrix[3]={30.,10000.5,100.};
 362
 363   fGenEventMatrix = new THnSparseF("fGenEventMatrix","mcZv:multMB:Centrality",3,binsEventMatrix,minEventMatrix,maxEventMatrix);
 364   fGenEventMatrix->SetBinEdges(0,binsZv);
 365   fGenEventMatrix->SetBinEdges(1,binsMult);
 366   fGenEventMatrix->SetBinEdges(2,binsCentrality);
 367   fGenEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 368   fGenEventMatrix->GetAxis(1)->SetTitle("multiplicity MB");
 369   fGenEventMatrix->GetAxis(2)->SetTitle("Centrality");
 370   fGenEventMatrix->Sumw2();
 371   //
 372   fTriggerEventMatrix = new THnSparseF("fTriggerEventMatrix","mcZv:multMB:Centrality",3,binsEventMatrix,minEventMatrix,maxEventMatrix);
 373   fTriggerEventMatrix->SetBinEdges(0,binsZv);
 374   fTriggerEventMatrix->SetBinEdges(1,binsMult);
 375   fTriggerEventMatrix->SetBinEdges(2,binsCentrality);
 376   fTriggerEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 377   fTriggerEventMatrix->GetAxis(1)->SetTitle("multiplicity MB");
 378   fTriggerEventMatrix->GetAxis(2)->SetTitle("Centrality");
 379   fTriggerEventMatrix->Sumw2();
 380   //
 381   fRecEventMatrix = new THnSparseF("fRecEventMatrix","mcZv:multMB:Centrality",3,binsEventMatrix,minEventMatrix,maxEventMatrix);
 382   fRecEventMatrix->SetBinEdges(0,binsZv);
 383   fRecEventMatrix->SetBinEdges(1,binsMult);
 384   fRecEventMatrix->SetBinEdges(2,binsCentrality);
 385   fRecEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 386   fRecEventMatrix->GetAxis(1)->SetTitle("multiplicity MB");
 387   fRecEventMatrix->GetAxis(2)->SetTitle("Centrality");
 388   fRecEventMatrix->Sumw2();
 389
 390   //
 391   // track to event corrections
 392   //
 393   fGenTrackEventMatrix = new THnSparseF("fGenTrackEventMatrix","mcZv:mcPt:mcEta:Centrality",4,binsTrackEventCorrMatrix);
 394   fGenTrackEventMatrix->SetBinEdges(0,binsZv);
 395   fGenTrackEventMatrix->SetBinEdges(1,binsPtTrackEventCorr);
 396   fGenTrackEventMatrix->SetBinEdges(2,binsEta);
 397   fGenTrackEventMatrix->SetBinEdges(3,binsCentrality);
 398   fGenTrackEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 399   fGenTrackEventMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 400   fGenTrackEventMatrix->GetAxis(2)->SetTitle("mcEta");
 401   fGenTrackEventMatrix->GetAxis(3)->SetTitle("Centrality");
 402   fGenTrackEventMatrix->Sumw2();
 403   //
 404   fTriggerTrackEventMatrix = new THnSparseF("fTriggerTrackEventMatrix","mcZv:mcPt:mcEta:Centrality",4,binsTrackEventCorrMatrix);
 405   fTriggerTrackEventMatrix->SetBinEdges(0,binsZv);
 406   fTriggerTrackEventMatrix->SetBinEdges(1,binsPtTrackEventCorr);
 407   fTriggerTrackEventMatrix->SetBinEdges(2,binsEta);
 408   fTriggerTrackEventMatrix->SetBinEdges(3,binsCentrality);
 409   fTriggerTrackEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 410   fTriggerTrackEventMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 411   fTriggerTrackEventMatrix->GetAxis(2)->SetTitle("mcEta");
 412   fTriggerTrackEventMatrix->GetAxis(3)->SetTitle("Centrality");
 413   fTriggerTrackEventMatrix->Sumw2();
 414   //
 415   fRecTrackEventMatrix = new THnSparseF("fRecTrackEventMatrix","mcZv:mcPt:mcEta:Centrality",4,binsTrackEventCorrMatrix);
 416   fRecTrackEventMatrix->SetBinEdges(0,binsZv);
 417   fRecTrackEventMatrix->SetBinEdges(1,binsPtTrackEventCorr);
 418   fRecTrackEventMatrix->SetBinEdges(2,binsEta);
 419   fRecTrackEventMatrix->SetBinEdges(3,binsCentrality);
 420   fRecTrackEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 421   fRecTrackEventMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 422   fRecTrackEventMatrix->GetAxis(2)->SetTitle("mcEta");
 423   fRecTrackEventMatrix->GetAxis(3)->SetTitle("Centrality");
 424   fRecTrackEventMatrix->Sumw2();
 425
 426   //
 427   // tracks correction matrices
 428   //
 429   fGenTrackMatrix = new THnSparseF("fGenTrackMatrix","mcZv:mcPt:mcEta:Centrality",4,binsTrackEventCorrMatrix);
 430   fGenTrackMatrix->SetBinEdges(0,binsZv);
 431   fGenTrackMatrix->SetBinEdges(1,binsPtTrackEventCorr);
 432   fGenTrackMatrix->SetBinEdges(2,binsEta);
 433   fGenTrackMatrix->SetBinEdges(3,binsCentrality);
 434   fGenTrackMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 435   fGenTrackMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 436   fGenTrackMatrix->GetAxis(2)->SetTitle("mcEta");
 437   fGenTrackMatrix->GetAxis(3)->SetTitle("Centrality");
 438   fGenTrackMatrix->Sumw2();
 439
 440   fGenPrimTrackMatrix = new THnSparseF("fGenPrimTrackMatrix","mcZv:mcPt:mcEta:Centrality",4,binsTrackEventCorrMatrix);
 441   fGenPrimTrackMatrix->SetBinEdges(0,binsZv);
 442   fGenPrimTrackMatrix->SetBinEdges(1,binsPtTrackEventCorr);
 443   fGenPrimTrackMatrix->SetBinEdges(2,binsEta);
 444   fGenPrimTrackMatrix->SetBinEdges(3,binsCentrality);
 445   fGenPrimTrackMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 446   fGenPrimTrackMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 447   fGenPrimTrackMatrix->GetAxis(2)->SetTitle("mcEta");
 448   fGenPrimTrackMatrix->GetAxis(3)->SetTitle("Centrality");
 449   fGenPrimTrackMatrix->Sumw2();
 450
 451
 452   fRecPrimTrackMatrix = new THnSparseF("fRecPrimTrackMatrix","mcZv:mcPt:mcEta:Centrality",4,binsTrackEventCorrMatrix);
 453   fRecPrimTrackMatrix->SetBinEdges(0,binsZv);
 454   fRecPrimTrackMatrix->SetBinEdges(1,binsPtTrackEventCorr);
 455   fRecPrimTrackMatrix->SetBinEdges(2,binsEta);
 456   fRecPrimTrackMatrix->SetBinEdges(3,binsCentrality);
 457   fRecPrimTrackMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 458   fRecPrimTrackMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 459   fRecPrimTrackMatrix->GetAxis(2)->SetTitle("mcEta");
 460   fRecPrimTrackMatrix->GetAxis(3)->SetTitle("Centrality");
 461   fRecPrimTrackMatrix->Sumw2();
 462
 463   //
 464   fRecTrackMatrix = new THnSparseF("fRecTrackMatrix","mcZv:mcPt:mcEta:Centrality",4,binsTrackEventCorrMatrix);
 465   fRecTrackMatrix->SetBinEdges(0,binsZv);
 466   fRecTrackMatrix->SetBinEdges(1,binsPtTrackEventCorr);
 467   fRecTrackMatrix->SetBinEdges(2,binsEta);
 468   fRecTrackMatrix->SetBinEdges(3,binsCentrality);
 469   fRecTrackMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 470   fRecTrackMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 471   fRecTrackMatrix->GetAxis(2)->SetTitle("mcEta");
 472   fRecTrackMatrix->GetAxis(3)->SetTitle("Centrality");
 473   fRecTrackMatrix->Sumw2();
 474
 475   fRecSecTrackMatrix = new THnSparseF("fRecSecTrackMatrix","mcZv:mcPt:mcEta:Centrality",4,binsTrackEventCorrMatrix);
 476   fRecSecTrackMatrix->SetBinEdges(0,binsZv);
 477   fRecSecTrackMatrix->SetBinEdges(1,binsPtTrackEventCorr);
 478   fRecSecTrackMatrix->SetBinEdges(2,binsEta);
 479   fRecSecTrackMatrix->SetBinEdges(3,binsCentrality);
 480   fRecSecTrackMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 481   fRecSecTrackMatrix->GetAxis(1)->SetTitle("Pt (GeV/c)");
 482   fRecSecTrackMatrix->GetAxis(2)->SetTitle("Eta");
 483   fRecSecTrackMatrix->GetAxis(3)->SetTitle("Centrality");
 484   fRecSecTrackMatrix->Sumw2();
 485
 486   //
 487   fRecMultTrackMatrix = new THnSparseF("fRecMultTrackMatrix","mcZv:mcPt:mcEta:Centrality",4,binsTrackEventCorrMatrix);
 488   fRecMultTrackMatrix->SetBinEdges(0,binsZv);
 489   fRecMultTrackMatrix->SetBinEdges(1,binsPtTrackEventCorr);
 490   fRecMultTrackMatrix->SetBinEdges(2,binsEta);
 491   fRecMultTrackMatrix->SetBinEdges(3,binsCentrality);
 492   fRecMultTrackMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 493   fRecMultTrackMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 494   fRecMultTrackMatrix->GetAxis(2)->SetTitle("mcEta");
 495   fRecMultTrackMatrix->GetAxis(3)->SetTitle("Centrality");
 496   fRecMultTrackMatrix->Sumw2();
 497
 498   //
 499   // Control analysis histograms
 500   //
 501   Int_t binsMCEventHist1[4]={100,100,zvNbins,centralityBins};
 502   Double_t minMCEventHist1[4]={-0.1,-0.1,-30.,0.};
 503   Double_t maxMCEventHist1[4]={0.1,0.1,30.,100.};
 504   fMCEventHist1 = new THnSparseF("fMCEventHist1","mcXv:mcYv:mcZv:Centrality",4,binsMCEventHist1,minMCEventHist1,maxMCEventHist1);
 505   fMCEventHist1->SetBinEdges(2,binsZv);
 506   fMCEventHist1->SetBinEdges(3,binsCentrality);
 507   fMCEventHist1->GetAxis(0)->SetTitle("mcXv (cm)");
 508   fMCEventHist1->GetAxis(1)->SetTitle("mcYv (cm)");
 509   fMCEventHist1->GetAxis(2)->SetTitle("mcZv (cm)");
 510   fMCEventHist1->GetAxis(3)->SetTitle("Centrality");
 511   fMCEventHist1->Sumw2();
 512
 513   //
 514   Int_t binsRecEventHist1[4]={100,100,zvNbins,centralityBins};
 515   Double_t minRecEventHist1[4]={-3.,-3.,-30.,0.};
 516   Double_t maxRecEventHist1[4]={3.,3.,30.,100.};
 517
 518   fRecEventHist1 = new THnSparseF("fRecEventHist1","Xv:Yv:Zv:Centrality",4,binsRecEventHist1,minRecEventHist1,maxRecEventHist1);
 519   fRecEventHist1->SetBinEdges(2,binsZv);
 520   fRecEventHist1->SetBinEdges(3,binsCentrality);
 521   fRecEventHist1->GetAxis(0)->SetTitle("Xv (cm)");
 522   fRecEventHist1->GetAxis(1)->SetTitle("Yv (cm)");
 523   fRecEventHist1->GetAxis(2)->SetTitle("Zv (cm)");
 524   fRecEventHist1->GetAxis(3)->SetTitle("Centrality");
 525   fRecEventHist1->Sumw2();
 526
 527   //
 528   Int_t binsRecEventHist2[3]={zvNbins, multNbins, centralityBins};
 529   Double_t minRecEventHist2[3]={-30., -0.5, 0.};
 530   Double_t maxRecEventHist2[3]={30., 10000.5, 100.};
 531
 532   fRecEventHist2 = new THnSparseF("fRecEventHist2","Zv:multMB:Centrality",3,binsRecEventHist2,minRecEventHist2,maxRecEventHist2);
 533   fRecEventHist2->SetBinEdges(0,binsZv);
 534   fRecEventHist2->SetBinEdges(1,binsMult);
 535   fRecEventHist2->SetBinEdges(2,binsCentrality);
 536   fRecEventHist2->GetAxis(0)->SetTitle("Zv (cm)");
 537   fRecEventHist2->GetAxis(1)->SetTitle("multiplicity MB");
 538   fRecEventHist2->GetAxis(2)->SetTitle("Centrality");
 539   fRecEventHist2->Sumw2();
 540
 541   //
 542   Double_t kFact = 0.1;
 543   Int_t binsRecMCEventHist1[4]={100,100,100, centralityBins};
 544   Double_t minRecMCEventHist1[4]={-10.0*kFact,-10.0*kFact,-10.0*kFact, 0.};
 545   Double_t maxRecMCEventHist1[4]={10.0*kFact,10.0*kFact,10.0*kFact, 100.};
 546
 547   fRecMCEventHist1 = new THnSparseF("fRecMCEventHist1","Xv-mcXv:Yv-mcYv:Zv-mcZv:Centrality",4,binsRecMCEventHist1,minRecMCEventHist1,maxRecMCEventHist1);
 548   fRecMCEventHist1->SetBinEdges(3,binsCentrality);
 549   fRecMCEventHist1->GetAxis(0)->SetTitle("Xv-mcXv (cm)");
 550   fRecMCEventHist1->GetAxis(1)->SetTitle("Yv-mcYv (cm)");
 551   fRecMCEventHist1->GetAxis(2)->SetTitle("Zv-mcZv (cm)");
 552   fRecMCEventHist1->GetAxis(3)->SetTitle("Centrality");
 553
 554   fRecMCEventHist1->Sumw2();
 555
 556   //
 557   Int_t binsRecMCEventHist2[4]={100,100,multNbins, centralityBins};
 558   Double_t minRecMCEventHist2[4]={-10.0*kFact,-10.0*kFact,-0.5, 0.};
 559   Double_t maxRecMCEventHist2[4]={10.0*kFact,10.0*kFact,10000.5, 100.};
 560
 561   fRecMCEventHist2 = new THnSparseF("fRecMCEventHist2","Xv-mcXv:Zv-mcZv:mult:Centrality",4,binsRecMCEventHist2,minRecMCEventHist2,maxRecMCEventHist2);
 562   fRecMCEventHist2->SetBinEdges(2,binsMult);
 563   fRecMCEventHist2->SetBinEdges(3,binsCentrality);
 564   fRecMCEventHist2->GetAxis(0)->SetTitle("Xv-mcXv (cm)");
 565   fRecMCEventHist2->GetAxis(1)->SetTitle("Zv-mcZv (cm)");
 566   fRecMCEventHist2->GetAxis(2)->SetTitle("multiplicity");
 567   fRecMCEventHist2->GetAxis(3)->SetTitle("Centrality");
 568   fRecMCEventHist2->Sumw2();
 569   //
 570   char name[256];
 571   char title[256];
 572   for(Int_t i=0; i<AlidNdPtHelper::kCutSteps; i++)
 573   {
 574   // THnSparse track histograms
 575
 576   Int_t binsMCTrackHist1[4]=  {ptNbins, etaNbins, 90, centralityBins};
 577   Double_t minMCTrackHist1[4]={0.,-1.5,0., 0.};
 578   Double_t maxMCTrackHist1[4]={50,1.5,2.*TMath::Pi(), 100.};
 579   sprintf(name,"fMCTrackHist1_%d",i);
 580   sprintf(title,"mcPt:mcEta:mcPhi:Centrality");
 581
 582   fMCTrackHist1[i] = new THnSparseF(name,title,4,binsMCTrackHist1,minMCTrackHist1,maxMCTrackHist1);
 583   fMCTrackHist1[i]->SetBinEdges(0,binsPt);
 584   fMCTrackHist1[i]->SetBinEdges(1,binsEta);
 585   fMCTrackHist1[i]->SetBinEdges(3,binsCentrality);
 586   fMCTrackHist1[i]->GetAxis(0)->SetTitle("mcPt (GeV/c)");
 587   fMCTrackHist1[i]->GetAxis(1)->SetTitle("mcEta");
 588   fMCTrackHist1[i]->GetAxis(2)->SetTitle("mcPhi (rad)");
 589   fMCTrackHist1[i]->GetAxis(3)->SetTitle("Centrality");
 590   fMCTrackHist1[i]->Sumw2();
 591
 592   Int_t binsMCPrimTrackHist1[6]=  {ptNbins,etaNbins,6,20,4000, centralityBins};
 593   Double_t minMCPrimTrackHist1[6]={0.,-1.5,0.,0.,0., 0.};
 594   Double_t maxMCPrimTrackHist1[6]={50.,1.5,6.,20.,4000., 100.};
 595   sprintf(name,"fMCPrimTrackHist1_%d",i);
 596   sprintf(title,"mcPt:mcEta:pid:mech:mother:Centrality");
 597
 598   fMCPrimTrackHist1[i] = new THnSparseF(name,title,6,binsMCPrimTrackHist1,minMCPrimTrackHist1,maxMCPrimTrackHist1);
 599   fMCPrimTrackHist1[i]->SetBinEdges(0,binsPt);
 600   fMCPrimTrackHist1[i]->SetBinEdges(1,binsEta);
 601   fMCPrimTrackHist1[i]->SetBinEdges(5,binsCentrality);
 602   fMCPrimTrackHist1[i]->GetAxis(0)->SetTitle("mcPt (GeV/c)");
 603   fMCPrimTrackHist1[i]->GetAxis(1)->SetTitle("mcEta");
 604   fMCPrimTrackHist1[i]->GetAxis(2)->SetTitle("pid");
 605   fMCPrimTrackHist1[i]->GetAxis(3)->SetTitle("mech");
 606   fMCPrimTrackHist1[i]->GetAxis(4)->SetTitle("mother");
 607   fMCPrimTrackHist1[i]->GetAxis(5)->SetTitle("Centrality");
 608   fMCPrimTrackHist1[i]->Sumw2();
 609
 610   Int_t binsMCPrimTrackHist2[4]=  {4000,20,4000,centralityBins};
 611   Double_t minMCPrimTrackHist2[4]={0.,0.,0., 0.};
 612   Double_t maxMCPrimTrackHist2[4]={4000.,20.,4000., 100.};
 613   sprintf(name,"fMCPrimTrackHist2_%d",i);
 614   sprintf(title,"pdg:mech:mother:Centrality");
 615
 616   fMCPrimTrackHist2[i] = new THnSparseF(name,title,4,binsMCPrimTrackHist2,minMCPrimTrackHist2,maxMCPrimTrackHist2);
 617   fMCPrimTrackHist2[i]->SetBinEdges(3,binsCentrality);
 618   fMCPrimTrackHist2[i]->GetAxis(0)->SetTitle("pdg");
 619   fMCPrimTrackHist2[i]->GetAxis(1)->SetTitle("mech");
 620   fMCPrimTrackHist2[i]->GetAxis(2)->SetTitle("mother");
 621   fMCPrimTrackHist2[i]->GetAxis(3)->SetTitle("Centrality");
 622   fMCPrimTrackHist2[i]->Sumw2();
 623
 624   Int_t binsMCSecTrackHist1[6]=  {ptNbins,etaNbins,6,20,4000, centralityBins};
 625   Double_t minMCSecTrackHist1[6]={0.,-1.5,0.,0.,0., 0.};
 626   Double_t maxMCSecTrackHist1[6]={50.,1.5,6.,20.,4000., 100.};
 627   sprintf(name,"fMCSecTrackHist1_%d",i);
 628   sprintf(title,"mcPt:mcEta:pid:mech:mother:Centrality");
 629
 630   fMCSecTrackHist1[i] = new THnSparseF(name,title,6,binsMCSecTrackHist1,minMCSecTrackHist1,maxMCSecTrackHist1);
 631   fMCSecTrackHist1[i]->SetBinEdges(0,binsPt);
 632   fMCSecTrackHist1[i]->SetBinEdges(1,binsEta);
 633   fMCSecTrackHist1[i]->SetBinEdges(5,binsCentrality);
 634   fMCSecTrackHist1[i]->GetAxis(0)->SetTitle("mcPt (GeV/c)");
 635   fMCSecTrackHist1[i]->GetAxis(1)->SetTitle("mcEta");
 636   fMCSecTrackHist1[i]->GetAxis(2)->SetTitle("pid");
 637   fMCSecTrackHist1[i]->GetAxis(3)->SetTitle("mech");
 638   fMCSecTrackHist1[i]->GetAxis(4)->SetTitle("mother");
 639   fMCSecTrackHist1[i]->GetAxis(5)->SetTitle("Centrality");
 640   fMCSecTrackHist1[i]->Sumw2();
 641
 642   Int_t binsRecTrackHist1[4]={ptNbins,etaNbins,90, centralityBins};
 643   Double_t minRecTrackHist1[4]={0.,-1.5,0., 0.};
 644   Double_t maxRecTrackHist1[4]={50.,1.5,2.*TMath::Pi(), 100.};
 645   sprintf(name,"fRecTrackHist1_%d",i);
 646   sprintf(title,"Pt:Eta:Phi:Centrality");
 647   fRecTrackHist1[i] = new THnSparseF(name,title,4,binsRecTrackHist1,minRecTrackHist1,maxRecTrackHist1);
 648   fRecTrackHist1[i]->SetBinEdges(0,binsPt);
 649   fRecTrackHist1[i]->SetBinEdges(1,binsEta);
 650   fRecTrackHist1[i]->SetBinEdges(3,binsCentrality);
 651   fRecTrackHist1[i]->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
 652   fRecTrackHist1[i]->GetAxis(1)->SetTitle("#eta");
 653   fRecTrackHist1[i]->GetAxis(2)->SetTitle("#phi (rad)");
 654   fRecTrackHist1[i]->GetAxis(3)->SetTitle("Centrality");
 655   fRecTrackHist1[i]->Sumw2();
 656
 657   sprintf(name,"fRecTrackHist2_%d",i);
 658   sprintf(title,"Zv:Pt:Eta:Centrality");
 659   fRecTrackHist2[i] = new THnSparseF(name,title,4,binsTrackEvent);
 660   fRecTrackHist2[i]->SetBinEdges(0,binsZv);
 661   fRecTrackHist2[i]->SetBinEdges(1,binsPt);
 662   fRecTrackHist2[i]->SetBinEdges(2,binsEta);
 663   fRecTrackHist2[i]->SetBinEdges(3,binsCentrality);
 664   fRecTrackHist2[i]->GetAxis(0)->SetTitle("Zv (cm)");
 665   fRecTrackHist2[i]->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
 666   fRecTrackHist2[i]->GetAxis(2)->SetTitle("#eta");
 667   fRecTrackHist2[i]->GetAxis(3)->SetTitle("Centrality");
 668   fRecTrackHist2[i]->Sumw2();
 669
 670   //
 671   Int_t binsRecTrackMultHist1[3]={ptNbins,multNbins, centralityBins};
 672   Double_t minRecTrackMultHist1[3]={0.,-0.5, -0.};
 673   Double_t maxRecTrackMultHist1[3]={50.,10000.5, 100.};
 674   sprintf(name,"fRecTrackMultHist_%d",i);
 675   sprintf(title,"Pt:Mult:Centrality");
 676   fRecTrackMultHist1[i] = new THnSparseF(name,title,3,binsRecTrackMultHist1,minRecTrackMultHist1,maxRecTrackMultHist1);
 677   fRecTrackMultHist1[i]->SetBinEdges(0,binsPt);
 678   fRecTrackMultHist1[i]->SetBinEdges(1,binsMult);
 679   fRecTrackMultHist1[i]->SetBinEdges(2,binsCentrality);
 680   fRecTrackMultHist1[i]->GetAxis(0)->SetTitle("Pt (GeV/c)");
 681   fRecTrackMultHist1[i]->GetAxis(1)->SetTitle("multiplicity");
 682   fRecTrackMultHist1[i]->GetAxis(2)->SetTitle("Centrality");
 683   fRecTrackMultHist1[i]->Sumw2();
 684   }
 685
 686   Int_t binsRecMCTrackHist1[5] = {ptNbins,etaNbins,100,100, centralityBins};
 687   Double_t minRecMCTrackHist1[5]={0.,-1.5,-0.5,-0.5, 0.};
 688   Double_t maxRecMCTrackHist1[5]={50.,1.5,0.5,0.5, 100.};
 689
 690   sprintf(name,"fRecMCTrackHist1");
 691   sprintf(title,"mcPt:mcEta:(Pt-mcPt)/mcPt:(Eta-mcEta):Centrality");
 692   fRecMCTrackHist1 = new THnSparseF(name,title,5,binsRecMCTrackHist1,minRecMCTrackHist1,maxRecMCTrackHist1);
 693   fRecMCTrackHist1->SetBinEdges(0,binsPt);
 694   fRecMCTrackHist1->SetBinEdges(1,binsEta);
 695   fRecMCTrackHist1->SetBinEdges(4,binsCentrality);
 696   fRecMCTrackHist1->GetAxis(0)->SetTitle("mcPt (GeV/c)");
 697   fRecMCTrackHist1->GetAxis(1)->SetTitle("mcEta");
 698   fRecMCTrackHist1->GetAxis(2)->SetTitle("(Pt-mcPt)/mcPt");
 699   fRecMCTrackHist1->GetAxis(3)->SetTitle("Eta-mcEta");
 700   fRecMCTrackHist1->GetAxis(4)->SetTitle("Centrality");
 701
 702   Int_t binsMCMultRecTrackHist1[4] = {ptNbins,etaNbins,6, centralityBins};
 703   Double_t minMCMultRecTrackHist1[4]={0.,-1.5,0., 0.};
 704   Double_t maxMCMultRecTrackHist1[4]={50.,1.5,6., 100.};
 705   sprintf(name,"fMCMultRecTrackHist1");
 706   sprintf(title,"mcPt:mcEta:pid:Centrality");
 707   fMCMultRecTrackHist1 = new THnSparseF(name,title,4,binsMCMultRecTrackHist1,minMCMultRecTrackHist1,maxMCMultRecTrackHist1);
 708   fMCMultRecTrackHist1->SetBinEdges(0,binsPt);
 709   fMCMultRecTrackHist1->SetBinEdges(1,binsEta);
 710   fMCMultRecTrackHist1->SetBinEdges(3,binsCentrality);
 711   fMCMultRecTrackHist1->GetAxis(0)->SetTitle("mcPt (GeV/c)");
 712   fMCMultRecTrackHist1->GetAxis(1)->SetTitle("mcEta");
 713   fMCMultRecTrackHist1->GetAxis(2)->SetTitle("pid");
 714   fMCMultRecTrackHist1->GetAxis(3)->SetTitle("Centrality");
 715
 716   //nClust:chi2PerClust:pt:eta:phi:Centrality
 717   Int_t binsRecTrackHist3[6]={160,100,ptNbins,etaNbins,90, centralityBins};
 718   Double_t minRecTrackHist3[6]={0., 0., 0., -1.5, 0., 0.};
 719   Double_t maxRecRecTrackHist3[6]={160.,10., 50., 1.5, 2.*TMath::Pi(), 100.};
 720
 721   fRecTrackHist3 = new THnSparseF("fRecTrackHist3","nClust:chi2PerClust:pt:eta:phi:Centrality",6,binsRecTrackHist3,minRecTrackHist3,maxRecRecTrackHist3);
 722   fRecTrackHist3->SetBinEdges(2,binsPt);
 723   fRecTrackHist3->SetBinEdges(3,binsEta);
 724   fRecTrackHist3->SetBinEdges(5,binsCentrality);
 725   fRecTrackHist3->GetAxis(0)->SetTitle("nClust");
 726   fRecTrackHist3->GetAxis(1)->SetTitle("chi2PerClust");
 727   fRecTrackHist3->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
 728   fRecTrackHist3->GetAxis(3)->SetTitle("#eta");
 729   fRecTrackHist3->GetAxis(4)->SetTitle("#phi (rad)");
 730   fRecTrackHist3->GetAxis(5)->SetTitle("Centrality");
 731   fRecTrackHist3->Sumw2();
 732
 733
 734   // init folder
 735   fAnalysisFolder = CreateFolder("folderdNdPt","Analysis dNdPt Folder");
 736
 737   // init trigger analysis (for zdc cut)
 738   fTriggerAnalysis = new AliTriggerAnalysis;
 739
 740   // set init flag
 741   fIsInit = kTRUE;
 742 }
 743
 744 //_____________________________________________________________________________
 745 void AlidNdPtAnalysisPbPb::Process(AliESDEvent *const esdEvent, AliMCEvent *const mcEvent)
 746 {
 747   //  init if not done already
 748   if (!fIsInit) { Init(); }
 749   //
 750   // Process real and/or simulated events
 751   //
 752   if(!esdEvent) {
 753     AliDebug(AliLog::kError, "esdEvent not available");
 754     return;
 755   }
 756
 757   // get selection cuts
 758   AlidNdPtEventCuts *evtCuts = GetEventCuts();
 759   AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts();
 760   AlidNdPtAcceptanceCuts *recCuts = GetRecAcceptanceCuts();
 761   AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
 762
 763   if(!evtCuts || !accCuts  || !esdTrackCuts) {
 764     AliDebug(AliLog::kError, "cuts not available");
 765     return;
 766   }
 767   if (0 == recCuts) { recCuts = accCuts;}
 768
 769   // trigger selection
 770   Bool_t isEventTriggered = kTRUE;
 771   AliPhysicsSelection *physicsSelection = NULL;
 772   AliTriggerAnalysis* triggerAnalysis = NULL;
 773
 774   //
 775   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
 776   if (!inputHandler)
 777   {
 778     Printf("ERROR: Could not receive input handler");
 779     return;
 780   }
 781
 782   if(evtCuts->IsTriggerRequired())
 783   {
 784     // always MB
 785     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
 786
 787     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
 788     if(!physicsSelection) return;
 789     //SetPhysicsTriggerSelection(physicsSelection);
 790
 791     if (isEventTriggered) {
 792       // set trigger (V0AND)
 793       triggerAnalysis = physicsSelection->GetTriggerAnalysis();
 794       if(!triggerAnalysis) return;
 795       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
 796     }
 797   }
 798
 799   /*
 800   // trigger selection
 801   Bool_t isEventTriggered = kTRUE;
 802   AliPhysicsSelection *trigSel = NULL;
 803   AliTriggerAnalysis *trigAna = NULL;
 804
 805   if(evtCuts->IsTriggerRequired())
 806   {
 807     //
 808     trigSel = GetPhysicsTriggerSelection();
 809     if(!trigSel) {
 810       printf("cannot get trigSel \n");
 811       return;
 812     }
 813
 814     //
 815     if(IsUseMCInfo())
 816     {
 817       trigSel->SetAnalyzeMC();
 818
 819
 820         isEventTriggered = trigSel->IsCollisionCandidate(esdEvent);
 821
 822         if(GetTrigger() == AliTriggerAnalysis::kV0AND)
 823         {
 824           trigAna = trigSel->GetTriggerAnalysis();
 825           if(!trigAna)
 826             return;
 827
 828           isEventTriggered = trigAna->IsOfflineTriggerFired(esdEvent, GetTrigger());
 829         }//if(GetTrigger() == AliTriggerAnalysis::kV0AND)
 830      }//if(IsUseMCInfo())
 831   }//if(evtCuts->IsTriggerRequired())
 832   */
 833
 834   // centrality determination
 835   Float_t centralityF = 0;
 836   AliCentrality *esdCentrality = esdEvent->GetCentrality();
 837   centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
 838   if (centralityF == 0.) {
 839     centralityF = 100.;
 840   }
 841
 842
 843   // use MC information
 844   AliHeader* header = 0;
 845   AliGenEventHeader* genHeader = 0;
 846   AliStack* stack = 0;
 847   TArrayF vtxMC(3);
 848
 849   Int_t multMCTrueTracks = 0;
 850   if(IsUseMCInfo())
 851   {
 852     //
 853     if(!mcEvent) {
 854       AliDebug(AliLog::kError, "mcEvent not available");
 855       return;
 856     }
 857     // get MC event header
 858     header = mcEvent->Header();
 859     if (!header) {
 860       AliDebug(AliLog::kError, "Header not available");
 861       return;
 862     }
 863     // MC particle stack
 864     stack = mcEvent->Stack();
 865     if (!stack) {
 866       AliDebug(AliLog::kError, "Stack not available");
 867       return;
 868     }
 869
 870     // get MC vertex
 871     genHeader = header->GenEventHeader();
 872     if (!genHeader) {
 873       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
 874       return;
 875     }
 876     genHeader->PrimaryVertex(vtxMC);
 877
 878     Double_t vMCEventHist1[4]={vtxMC[0],vtxMC[1],vtxMC[2],centralityF};
 879     fMCEventHist1->Fill(vMCEventHist1);
 880
 881     // multipliticy of all MC primary tracks
 882     // in Zv, pt and eta ranges)
 883     multMCTrueTracks = AlidNdPtHelper::GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
 884
 885   } // end bUseMC
 886
 887   // get reconstructed vertex
 888   const AliESDVertex* vtxESD = 0;
 889   if(evtCuts->IsRecVertexRequired())
 890   {
 891     if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {
 892       vtxESD = esdEvent->GetPrimaryVertexTPC();
 893     }
 894     else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {
 895       vtxESD = esdEvent->GetPrimaryVertexTracks();
 896     }
 897     else {
 898         return;
 899     }
 900   }
 901
 902   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
 903   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
 904   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
 905
 906   // vertex contributors
 907   Int_t multMBTracks = 0;
 908   if(GetAnalysisMode() == AlidNdPtHelper::kTPC)
 909   {
 910      if(vtxESD->GetStatus()) {
 911          multMBTracks = vtxESD->GetNContributors();
 912      }
 913   }
 914   else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {
 915      if(vtxESD->GetStatus()) {
 916          multMBTracks = vtxESD->GetNContributors();
 917      }
 918   }
 919   else {
 920     AliDebug(AliLog::kError, Form("Found analysis type %d", GetAnalysisMode()));
 921     return;
 922   }
 923
 924   TObjArray *allChargedTracks=0;
 925   //Int_t multAll=0, multAcc=0, multRec=0;
 926   Int_t multAll=0, multRec=0;
 927   Int_t *labelsAll=0, *labelsAcc=0, *labelsRec=0;
 928
 929   // check event cuts
 930   if(isEventOK && isEventTriggered)
 931   {
 932
 933     // get all charged tracks
 934     allChargedTracks = AlidNdPtHelper::GetAllChargedTracks(esdEvent,GetAnalysisMode());
 935     if(!allChargedTracks) return;
 936
 937     Int_t entries = allChargedTracks->GetEntries();
 938
 939     labelsAll = new Int_t[entries];
 940     labelsAcc = new Int_t[entries];
 941     labelsRec = new Int_t[entries];
 942     for(Int_t i=0; i<entries;++i)
 943     {
 944       AliESDtrack *track = (AliESDtrack*)allChargedTracks->At(i);
 945
 946       if(!track) continue;
 947       if(track->Charge()==0) continue;
 948
 949
 950       // only postive charged
 951       if(GetParticleMode() == AlidNdPtHelper::kPlus && track->Charge() < 0)
 952         continue;
 953
 954       // only negative charged
 955       if(GetParticleMode() == AlidNdPtHelper::kMinus && track->Charge() > 0)
 956         continue;
 957
 958       //
 959       Double_t values[4] = {vtxESD->GetZv(),track->Pt(),track->Eta(), centralityF};
 960
 961       fRecTrackHist2[AlidNdPtHelper::kAllTracks]->Fill(values);
 962       FillHistograms(track,stack,AlidNdPtHelper::kAllTracks,centralityF);
 963       labelsAll[multAll] = TMath::Abs(track->GetLabel());
 964
 965       multAll++;
 966       if(esdTrackCuts->AcceptTrack(track) && accCuts->AcceptTrack(track) && recCuts->AcceptTrack(track)) {
 967
 968          fRecTrackHist2[AlidNdPtHelper::kRecTracks]->Fill(values);
 969          FillHistograms(track,stack,AlidNdPtHelper::kRecTracks,centralityF);
 970          labelsRec[multRec] = TMath::Abs(track->GetLabel());
 971
 972          multRec++;
 973
 974       }
 975      }//loop over entries
 976
 977
 978      // fill track multiplicity histograms
 979      //FillHistograms(allChargedTracks,labelsAll,multAll,labelsAcc,multAcc,labelsRec,multRec);
 980
 981      Double_t vRecEventHist1[4] = {vtxESD->GetXv(),vtxESD->GetYv(),vtxESD->GetZv(),centralityF};
 982      fRecEventHist1->Fill(vRecEventHist1);
 983
 984      Double_t vRecEventHist2[3] = {vtxESD->GetZv(),multMBTracks,centralityF};
 985      fRecEventHist2->Fill(vRecEventHist2);
 986
 987    } // triggered and event vertex
 988
 989    if(IsUseMCInfo())
 990    {
 991
 992      //Double_t vMultTrueEventMatrix[2] = { multRec, multMCTrueTracks };
 993
 994      //
 995      // event level corrections (zv,N_MB)
 996      //
 997      // all inelastic
 998      Double_t vEventMatrix[3] = {vtxMC[2],multMBTracks,centralityF};
 999      fGenEventMatrix->Fill(vEventMatrix);
1000      if(isEventTriggered) fTriggerEventMatrix->Fill(vEventMatrix);
1001      if(isEventOK && isEventTriggered) fRecEventMatrix->Fill(vEventMatrix);
1002
1003      //
1004      // track-event level corrections (zv,pt,eta)
1005      //
1006      for (Int_t iMc = 0; iMc < stack->GetNtrack(); ++iMc)
1007      {
1008        TParticle* particle = stack->Particle(iMc);
1009        if (!particle)
1010        continue;
1011
1012        // only charged particles
1013        if(!particle->GetPDG()) continue;
1014        Double_t charge = particle->GetPDG()->Charge()/3.;
1015        if ( TMath::Abs(charge) < 0.001 )
1016         continue;
1017
1018        // only postive charged
1019        if(GetParticleMode() == AlidNdPtHelper::kPlus && charge < 0.)
1020         continue;
1021
1022        // only negative charged
1023        if(GetParticleMode() == AlidNdPtHelper::kMinus && charge > 0.)
1024        continue;
1025
1026        // physical primary
1027        Bool_t prim = stack->IsPhysicalPrimary(iMc);
1028        if(!prim) continue;
1029
1030        // checked accepted
1031        if(accCuts->AcceptTrack(particle))
1032        {
1033         Double_t vTrackEventMatrix[4] = {vtxMC[2], particle->Pt(), particle->Eta(), centralityF};
1034         fGenTrackEventMatrix->Fill(vTrackEventMatrix);
1035
1036         if(!isEventTriggered) continue;
1037         fTriggerTrackEventMatrix->Fill(vTrackEventMatrix);
1038
1039         if(!isEventOK) continue;
1040         fRecTrackEventMatrix->Fill(vTrackEventMatrix);
1041
1042        }// if(accCuts->AcceptTrack(particle))
1043      }// for (Int_t iMc = 0; iMc < stack->GetNtrack(); ++iMc)
1044
1045      //
1046      // track-level corrections (zv,pt,eta)
1047      //
1048      if(isEventOK && isEventTriggered)
1049      {
1050
1051        // fill MC and rec event control histograms
1052        if(fHistogramsOn) {
1053          Double_t vRecMCEventHist1[4] = {vtxESD->GetXv()-vtxMC[0],vtxESD->GetYv()-vtxMC[1],vtxESD->GetZv()-vtxMC[2], centralityF};
1054          fRecMCEventHist1->Fill(vRecMCEventHist1);//
1055
1056          Double_t vRecMCEventHist2[4] = {vtxESD->GetXv()-vtxMC[0],vtxESD->GetZv()-vtxMC[2],multMBTracks, centralityF};
1057          fRecMCEventHist2->Fill(vRecMCEventHist2);
1058
1059        }//
1060
1061        //
1062        // MC histograms for track efficiency studies
1063        //
1064        for (Int_t iMc = 0; iMc < stack->GetNtrack(); ++iMc)
1065        {
1066          TParticle* particle = stack->Particle(iMc);
1067          if (!particle)
1068          continue;
1069
1070          Double_t vTrackMatrix[4] = {vtxMC[2],particle->Pt(),particle->Eta(),centralityF};
1071
1072          // only charged particles
1073          if(!particle->GetPDG()) continue;
1074          Double_t charge = particle->GetPDG()->Charge()/3.;
1075          if (TMath::Abs(charge) < 0.001)
1076          continue;
1077
1078          // only postive charged
1079          if(GetParticleMode() == AlidNdPtHelper::kPlus && charge < 0.)
1080          continue;
1081
1082          // only negative charged
1083          if(GetParticleMode() == AlidNdPtHelper::kMinus && charge > 0.)
1084          continue;
1085
1086          // physical primary
1087          Bool_t prim = stack->IsPhysicalPrimary(iMc);
1088
1089          // check accepted
1090          if(accCuts->AcceptTrack(particle))
1091          {
1092
1093            if( AlidNdPtHelper::IsPrimaryParticle(stack, iMc, GetParticleMode()) )
1094              fGenPrimTrackMatrix->Fill(vTrackMatrix);
1095
1096            // fill control histograms
1097            if(fHistogramsOn)
1098              FillHistograms(stack,iMc,AlidNdPtHelper::kAccTracks, centralityF);
1099
1100            // check multiple found tracks
1101            Int_t multCount = 0;
1102            for(Int_t iRec=0; iRec<multRec; ++iRec)
1103            {
1104              if(iMc == labelsRec[iRec])
1105              {
1106                multCount++;
1107                if(multCount>1)
1108                {
1109                  fRecMultTrackMatrix->Fill(vTrackMatrix);
1110
1111                  // fill control histogram
1112                  if(fHistogramsOn) {
1113                    Int_t pid = AlidNdPtHelper::ConvertPdgToPid(particle);
1114                    Double_t vMCMultRecTrackHist1[4] = {particle->Pt(), particle->Eta(), pid, centralityF};
1115                    fMCMultRecTrackHist1->Fill(vMCMultRecTrackHist1);
1116                  }
1117                }
1118              }
1119            }
1120
1121            // check reconstructed
1122            for(Int_t iRec=0; iRec<multRec; ++iRec)
1123            {
1124              if(iMc == labelsRec[iRec])
1125              {
1126                fRecTrackMatrix->Fill(vTrackMatrix);
1127
1128                if( AlidNdPtHelper::IsPrimaryParticle(stack, iMc, GetParticleMode()) ) {
1129                  fRecPrimTrackMatrix->Fill(vTrackMatrix);
1130                }
1131                if(!prim) fRecSecTrackMatrix->Fill(vTrackMatrix);
1132
1133                // fill control histograms
1134                if(fHistogramsOn)
1135                  FillHistograms(stack,iMc,AlidNdPtHelper::kRecTracks, centralityF);
1136
1137                break;
1138              }//if(iMc == labelsRec[iRec])
1139            }//reco tracks
1140          }//accepted tracks
1141        }//stack loop
1142      }//is triggered
1143    } // end bUseMC
1144
1145   if(allChargedTracks) delete allChargedTracks; allChargedTracks = 0;
1146   if(labelsAll) delete [] labelsAll; labelsAll = 0;
1147   if(labelsAcc) delete [] labelsAcc; labelsAcc = 0;
1148   if(labelsRec) delete [] labelsRec; labelsRec = 0;
1149
1150   if(!evtCuts->IsRecVertexRequired() && vtxESD != NULL) delete vtxESD;
1151
1152 }
1153
1154 //_____________________________________________________________________________
1155 void AlidNdPtAnalysisPbPb::FillHistograms(TObjArray *const allChargedTracks,Int_t *const labelsAll,Int_t multAll,Int_t *const labelsAcc,Int_t multAcc,Int_t *const labelsRec,Int_t multRec, Float_t centralityF) {
1156  // multiplicity  histograms
1157
1158
1159   if(!allChargedTracks) return;
1160   if(!labelsAll) return;
1161   if(!labelsAcc) return;
1162   if(!labelsRec) return;
1163
1164   Int_t entries = allChargedTracks->GetEntries();
1165   for(Int_t i=0; i<entries; ++i)
1166   {
1167      AliESDtrack *track = (AliESDtrack*)allChargedTracks->At(i);
1168      if(!track) continue;
1169      if(track->Charge() == 0) continue;
1170
1171      Int_t label = TMath::Abs(track->GetLabel());
1172      for(Int_t iAll=0; iAll<multAll; ++iAll) {
1173        if(label == labelsAll[iAll]) {
1174          Double_t v1[3] = {track->Pt(), multAll, centralityF};
1175          fRecTrackMultHist1[AlidNdPtHelper::kAllTracks]->Fill(v1);
1176        }
1177      }
1178      for(Int_t iAcc=0; iAcc<multAcc; ++iAcc) {
1179        if(label == labelsAcc[iAcc]) {
1180          Double_t v2[3] = {track->Pt(), multAcc, centralityF};
1181          fRecTrackMultHist1[AlidNdPtHelper::kAccTracks]->Fill(v2);
1182        }
1183      }
1184      for(Int_t iRec=0; iRec<multRec; ++iRec) {
1185        if(label == labelsRec[iRec]) {
1186          Double_t v3[3] = {track->Pt(), multRec, centralityF};
1187          fRecTrackMultHist1[AlidNdPtHelper::kRecTracks]->Fill(v3);
1188        }//out
1189      }
1190   }
1191 }
1192
1193 //_____________________________________________________________________________
1194 void AlidNdPtAnalysisPbPb::FillHistograms(AliESDtrack *const esdTrack, AliStack *const stack, AlidNdPtHelper::TrackObject trackObj, Float_t centralityF)
1195 {
1196
1197   //
1198   // Fill ESD track and MC histograms
1199   //
1200   if(!esdTrack) return;
1201
1202   Float_t q = esdTrack->Charge();
1203   if(TMath::Abs(q) < 0.001) return;
1204
1205   Float_t pt = esdTrack->Pt();
1206   Float_t eta = esdTrack->Eta();
1207   Float_t phi = esdTrack->Phi();
1208
1209   Float_t dca[2], bCov[3];
1210   esdTrack->GetImpactParameters(dca,bCov);
1211
1212   Int_t nClust = esdTrack->GetTPCclusters(0);
1213   Float_t chi2PerCluster = 0.;
1214   if(nClust>0.) chi2PerCluster = esdTrack->GetTPCchi2()/Float_t(nClust);
1215
1216
1217   // fill histograms
1218   Double_t values[4] = {pt,eta,phi,centralityF};
1219   fRecTrackHist1[trackObj]->Fill(values);
1220
1221   //
1222   // Fill rec vs MC information
1223   //
1224   if(!stack) return;
1225
1226   Int_t label = TMath::Abs(esdTrack->GetLabel());
1227   //if(label == 0) return;
1228
1229   TParticle* particle = stack->Particle(label);
1230   if(!particle) return;
1231
1232   Int_t motherPdg = -1;
1233   TParticle* mother = 0;
1234
1235   //TParticle* prim_mother = AlidNdPtHelper::FindPrimaryMother(stack,label);
1236   Int_t motherLabel = particle->GetMother(0);
1237   if(motherLabel>0) mother = stack->Particle(motherLabel);
1238   if(mother) motherPdg = TMath::Abs(mother->GetPdgCode()); // take abs for visualisation only
1239   //Int_t mech = particle->GetUniqueID(); // production mechanism
1240
1241   if(!particle->GetPDG()) return;
1242   Double_t gq = particle->GetPDG()->Charge()/3.0; // Charge units |e|/3
1243   if(TMath::Abs(gq)<0.001) return;
1244   Float_t gpt = particle->Pt();
1245   Float_t geta = particle->Eta();
1246
1247   Double_t dpt=0;
1248   //printf("pt %f, gpt %f \n",pt,gpt);
1249   if(gpt) dpt = (pt-gpt)/gpt;
1250   Double_t deta = (eta-geta);
1251
1252   // fill histograms
1253   if(trackObj == AlidNdPtHelper::kRecTracks)  //RecTracks???
1254   {
1255     Double_t vTrackPtCorrelationMatrix[4]={pt,gpt,geta,centralityF};
1256     fTrackPtCorrelationMatrix->Fill(vTrackPtCorrelationMatrix);
1257
1258     Double_t vRecMCTrackHist1[5]={gpt,geta,dpt,deta,centralityF};
1259     fRecMCTrackHist1->Fill(vRecMCTrackHist1);
1260   }
1261 }
1262
1263 //_____________________________________________________________________________
1264 void AlidNdPtAnalysisPbPb::FillHistograms(AliStack *const stack, Int_t label, AlidNdPtHelper::TrackObject trackObj, Float_t centralityF)
1265 {
1266
1267   // Fill MC histograms
1268   if(!stack) return;
1269
1270   TParticle* particle = stack->Particle(label);
1271   if(!particle) return;
1272
1273   Int_t motherPdg = -1;
1274   TParticle* mother = 0;
1275
1276   //TParticle* prim_mother = AlidNdPtHelper::FindPrimaryMother(stack,label);
1277   Int_t motherLabel = particle->GetMother(0);
1278   if(motherLabel>0) mother = stack->Particle(motherLabel);
1279   if(mother) motherPdg = TMath::Abs(mother->GetPdgCode()); // take abs for visualisation only
1280   Int_t mech = particle->GetUniqueID(); // production mechanism
1281
1282   if(!particle->GetPDG()) return;
1283   Double_t gq = particle->GetPDG()->Charge()/3.0; // Charge units |e|/3
1284   if(TMath::Abs(gq) < 0.001) return;
1285
1286   Float_t gpt = particle->Pt();
1287   //Float_t qgpt = particle->Pt() * gq;
1288   Float_t geta = particle->Eta();
1289   Float_t gphi = particle->Phi();
1290   //Float_t gpz = particle->Pz();
1291
1292   Bool_t prim = stack->IsPhysicalPrimary(label);
1293   //Float_t vx = particle->Vx(); Float_t vy = particle->Vy(); Float_t vz = particle->Vz();
1294
1295   Int_t pid = AlidNdPtHelper::ConvertPdgToPid(particle);
1296
1297   //if(prim&&pid==5) printf("pdgcode %d, production mech %d \n",particle->GetPdgCode(),mech);
1298   //if(!prim) printf("motherPdg %d, particle %d, production mech %d\n",motherPdg, particle->GetPdgCode(),mech);
1299
1300   //
1301   // fill histogram
1302   //
1303   Double_t vMCTrackHist1[4] = {gpt,geta,gphi,centralityF};
1304   fMCTrackHist1[trackObj]->Fill(vMCTrackHist1);
1305
1306   Double_t vMCPrimTrackHist1[6] = {gpt,geta,pid,mech,motherPdg,centralityF};
1307   Double_t vMCPrimTrackHist2[4] = {TMath::Abs(particle->GetPdgCode()),mech,motherPdg,centralityF};
1308
1309   //if(prim && AliPWG0Helper::IsPrimaryCharged(particle, label)) fMCPrimTrackHist1[trackObj]->Fill(vMCPrimTrackHist1);
1310
1311   if(prim) {
1312     fMCPrimTrackHist1[trackObj]->Fill(vMCPrimTrackHist1);
1313     if(pid == 5) fMCPrimTrackHist2[trackObj]->Fill(vMCPrimTrackHist2);
1314   }
1315   else {
1316     fMCSecTrackHist1[trackObj]->Fill(vMCPrimTrackHist1);
1317   }
1318
1319 }
1320
1321 //_____________________________________________________________________________
1322 Long64_t AlidNdPtAnalysisPbPb::Merge(TCollection* const list)
1323 {
1324   // Merge list of objects (needed by PROOF)
1325
1326   //  init if not done already
1327   if (!fIsInit) { Init(); }
1328
1329   if (!list)
1330   return 0;
1331
1332   if (list->IsEmpty())
1333   return 1;
1334
1335   TIterator* iter = list->MakeIterator();
1336   TObject* obj = 0;
1337
1338   //
1339   //TList *collPhysSelection = new TList;
1340
1341   // collection of generated histograms
1342
1343   Int_t count=0;
1344   while((obj = iter->Next()) != 0) {
1345     AlidNdPtAnalysisPbPb* entry = dynamic_cast<AlidNdPtAnalysisPbPb*>(obj);
1346     if (entry == 0) continue;
1347
1348     // physics selection
1349     //printf("entry->GetPhysicsTriggerSelection() %p \n", entry->GetPhysicsTriggerSelection());
1350     //AliPhysicsSelection *physSel = entry->GetPhysicsTriggerSelection();
1351     //if( physSel ) collPhysSelection->Add(physSel);
1352
1353     //
1354     fTrackPtCorrelationMatrix->Add(entry->fTrackPtCorrelationMatrix);
1355
1356     //
1357     fGenEventMatrix->Add(entry->fGenEventMatrix);
1358
1359     fTriggerEventMatrix->Add(entry->fTriggerEventMatrix);
1360
1361     fRecEventMatrix->Add(entry->fRecEventMatrix);
1362     //
1363     fGenTrackEventMatrix->Add(entry->fGenTrackEventMatrix);
1364
1365     fTriggerTrackEventMatrix->Add(entry->fTriggerTrackEventMatrix);
1366
1367     fRecTrackEventMatrix->Add(entry->fRecTrackEventMatrix);
1368
1369     //
1370     fGenTrackMatrix->Add(entry->fGenTrackMatrix);
1371     fGenPrimTrackMatrix->Add(entry->fGenPrimTrackMatrix);
1372     fRecPrimTrackMatrix->Add(entry->fRecPrimTrackMatrix);
1373     //
1374     fRecTrackMatrix->Add(entry->fRecTrackMatrix);
1375     fRecSecTrackMatrix->Add(entry->fRecSecTrackMatrix);
1376     //
1377     fRecMultTrackMatrix->Add(entry->fRecMultTrackMatrix);
1378
1379     //
1380     // control analysis histograms
1381     //
1382     fMCEventHist1->Add(entry->fMCEventHist1);
1383     fRecEventHist1->Add(entry->fRecEventHist1);
1384     fRecEventHist2->Add(entry->fRecEventHist2);
1385     fRecMCEventHist1->Add(entry->fRecMCEventHist1);
1386     fRecMCEventHist2->Add(entry->fRecMCEventHist2);
1387
1388
1389     for(Int_t i=0; i<AlidNdPtHelper::kCutSteps; i++) {
1390       fMCTrackHist1[i]->Add(entry->fMCTrackHist1[i]);
1391
1392       fMCPrimTrackHist1[i]->Add(entry->fMCPrimTrackHist1[i]);
1393       fMCPrimTrackHist2[i]->Add(entry->fMCPrimTrackHist2[i]);
1394       fMCSecTrackHist1[i]->Add(entry->fMCSecTrackHist1[i]);
1395
1396       fRecTrackHist1[i]->Add(entry->fRecTrackHist1[i]);
1397       fRecTrackHist2[i]->Add(entry->fRecTrackHist2[i]);
1398       fRecTrackMultHist1[i]->Add(entry->fRecTrackMultHist1[i]);
1399     }
1400     fRecMCTrackHist1->Add(entry->fRecMCTrackHist1);
1401     fMCMultRecTrackHist1->Add(entry->fMCMultRecTrackHist1);
1402     fRecTrackHist3->Add(entry->fRecTrackHist3);
1403
1404   count++;
1405   }
1406
1407   //AliPhysicsSelection *trigSelection = GetPhysicsTriggerSelection();
1408   //if( trigSelection ) trigSelection->Merge(collPhysSelection);
1409   //if(collPhysSelection) delete collPhysSelection;
1410
1411 return count;
1412 }
1413
1414
1415
1416 //_____________________________________________________________________________
1417 void AlidNdPtAnalysisPbPb::Analyse()
1418 {
1419
1420   //  init if not done already
1421   if (!fIsInit) { Init(); }
1422
1423   // Analyse histograms
1424   //
1425   TH1::AddDirectory(kFALSE);
1426   TH1 *h=0, *h1=0, *h2=0, *h2c = 0;
1427   THnSparse *hs=0;
1428   TH2 *h2D=0;
1429
1430   char name[256];
1431   TObjArray *aFolderObj = new TObjArray;
1432
1433   //
1434   // LHC backgraund in all and 0-bins
1435   //
1436   //AliPhysicsSelection *trigSel = GetPhysicsTriggerSelection();
1437   //trigSel->SaveHistograms("physics_selection");
1438
1439   //
1440   // Reconstructed event vertex
1441   //
1442   h = fRecEventHist1->Projection(0);
1443   h->SetName("Xv");
1444   aFolderObj->Add(h);
1445
1446   h = fRecEventHist1->Projection(1);
1447   h->SetName("Yv");
1448   aFolderObj->Add(h);
1449
1450   h = fRecEventHist1->Projection(2);
1451   h->SetName("Zv");
1452   aFolderObj->Add(h);
1453
1454   //
1455   // multiplicity
1456   //
1457   h = fRecEventHist2->Projection(1);
1458   h->SetName("multMB");
1459   aFolderObj->Add(h);
1460
1461   h2D = fRecEventHist2->Projection(0,1);
1462   h2D->SetName("Zv_vs_multiplicity_MB");
1463   aFolderObj->Add(h2D);
1464
1465   //
1466   // reconstructed pt histograms
1467   //
1468   h = fRecTrackHist1[0]->Projection(0);
1469   h->Scale(1.,"width");
1470   h->SetName("pt_all_ch");
1471   aFolderObj->Add(h);
1472
1473   h = fRecTrackHist1[1]->Projection(0);
1474   h->Scale(1.,"width");
1475   h->SetName("pt_acc");
1476   aFolderObj->Add(h);
1477
1478   h = fRecTrackHist1[2]->Projection(0);
1479   h->Scale(1.,"width");
1480   h->SetName("pt_rec");
1481   aFolderObj->Add(h);
1482
1483   //
1484   // reconstructed eta histograms
1485   //
1486   h = fRecTrackHist1[0]->Projection(1);
1487   h->SetName("eta_all_ch");
1488   aFolderObj->Add(h);
1489
1490   h = fRecTrackHist1[1]->Projection(1);
1491   h->SetName("eta_acc");
1492   aFolderObj->Add(h);
1493
1494   h = fRecTrackHist1[2]->Projection(1);
1495   h->SetName("eta_rec");
1496   aFolderObj->Add(h);
1497
1498   //
1499   // reconstructed phi histograms
1500   //
1501   h = fRecTrackHist1[0]->Projection(2);
1502   h->SetName("phi_all_ch");
1503   aFolderObj->Add(h);
1504
1505   h = fRecTrackHist1[1]->Projection(2);
1506   h->SetName("phi_acc");
1507   aFolderObj->Add(h);
1508
1509   h = fRecTrackHist1[2]->Projection(2);
1510   h->SetName("phi_rec");
1511   aFolderObj->Add(h);
1512
1513   //
1514   // reconstructed eta:pt histograms
1515   //
1516   h2D = fRecTrackHist1[0]->Projection(1,0);
1517   h2D->SetName("pt_eta_all_ch");
1518   aFolderObj->Add(h2D);
1519
1520   h2D = fRecTrackHist1[1]->Projection(1,0);
1521   h2D->SetName("pt_eta_acc");
1522   aFolderObj->Add(h2D);
1523
1524   h2D = fRecTrackHist1[2]->Projection(1,0);
1525   h2D->SetName("pt_eta_rec");
1526   aFolderObj->Add(h2D);
1527
1528   //
1529   // reconstructed phi:pt histograms
1530   //
1531   h2D = fRecTrackHist1[0]->Projection(2,0);
1532   h2D->SetName("pt_phi_all_ch");
1533   aFolderObj->Add(h2D);
1534
1535   h2D = fRecTrackHist1[1]->Projection(2,0);
1536   h2D->SetName("pt_phi_acc");
1537   aFolderObj->Add(h2D);
1538
1539   h2D = fRecTrackHist1[2]->Projection(2,0);
1540   h2D->SetName("pt_phi_rec");
1541   aFolderObj->Add(h2D);
1542
1543   //
1544   // reconstructed phi:eta histograms
1545   //
1546   h2D = fRecTrackHist1[0]->Projection(2,1);
1547   h2D->SetName("eta_phi_all_ch");
1548   aFolderObj->Add(h2D);
1549
1550   h2D = fRecTrackHist1[1]->Projection(2,1);
1551   h2D->SetName("eta_phi_acc");
1552   aFolderObj->Add(h2D);
1553
1554   h2D = fRecTrackHist1[2]->Projection(2,1);
1555   h2D->SetName("eta_phi_rec");
1556   aFolderObj->Add(h2D);
1557
1558   //
1559   // reconstructed nClust, chi2 vs pt, eta, phi
1560   //
1561   if(fHistogramsOn) {
1562
1563     h2D = fRecTrackHist3->Projection(0,1);
1564     h2D->SetName("nClust_chi2_rec");
1565     aFolderObj->Add(h2D);
1566
1567     h2D = fRecTrackHist3->Projection(0,2);
1568     h2D->SetName("nClust_pt_rec");
1569     aFolderObj->Add(h2D);
1570
1571     h2D = fRecTrackHist3->Projection(0,3);
1572     h2D->SetName("nClust_eta_rec");
1573     aFolderObj->Add(h2D);
1574
1575     h2D = fRecTrackHist3->Projection(0,4);
1576     h2D->SetName("nClust_phi_rec");
1577     aFolderObj->Add(h2D);
1578
1579     h2D = fRecTrackHist3->Projection(1,2);
1580     h2D->SetName("chi2_pt_rec");
1581     aFolderObj->Add(h2D);
1582
1583     h2D = fRecTrackHist3->Projection(1,3);
1584     h2D->SetName("chi2_eta_rec");
1585     aFolderObj->Add(h2D);
1586
1587     h2D = fRecTrackHist3->Projection(1,4);
1588     h2D->SetName("chi2_phi_rec");
1589     aFolderObj->Add(h2D);
1590
1591   }
1592
1593   //
1594   // calculate corrections for empty events
1595   // with multMB==0
1596   //
1597
1598   //
1599   // normalised zv to generate zv for triggered events
1600   //
1601   h = fRecEventHist2->Projection(0);
1602   if( h->Integral() ) h->Scale(1./h->Integral());
1603   h->SetName("zv_distribution_norm");
1604   aFolderObj->Add(h);
1605
1606   //
1607   // MC available
1608   //
1609   if(IsUseMCInfo()) {
1610
1611   //
1612   // Event vertex resolution
1613   //
1614   h2D = fRecMCEventHist2->Projection(0,2);
1615   h2D->SetName("DeltaXv_vs_mult");
1616   aFolderObj->Add(h2D);
1617
1618   h2D = fRecMCEventHist2->Projection(1,2);
1619   h2D->SetName("DeltaZv_vs_mult");
1620   aFolderObj->Add(h2D);
1621
1622   //
1623   // normalised zv to get trigger/trigger+vertex event differences
1624   // F(zv) = E_trig(zv,0)/Int(E_trig(zv,0) / Sum(E_trigvtx(zv,n))/Sum(Int(E_trigvtx(zv,n))dzv)
1625   //
1626   fTriggerEventMatrix->GetAxis(1)->SetRangeUser(0.,0.);
1627   h = fTriggerEventMatrix->Projection(0);
1628   h2D = fTriggerEventMatrix->Projection(0,1);
1629   if(h2D->Integral()) h->Scale(1./h2D->Integral());
1630
1631   h1 = fRecEventMatrix->Projection(0);
1632   h2D = fRecEventMatrix->Projection(0,1);
1633   if(h2D->Integral()) h1->Scale(1./h2D->Integral());
1634
1635   h->Divide(h1);
1636   h->SetName("zv_empty_events_norm");
1637   aFolderObj->Add(h);
1638
1639   fTriggerEventMatrix->GetAxis(1)->SetRange(1,fTriggerEventMatrix->GetAxis(1)->GetNbins());
1640
1641   //
1642   // rec. vs true track pt correlation matrix
1643   //
1644   hs = (THnSparse*)fTrackPtCorrelationMatrix->Clone("track_pt_correlation_matrix");
1645   aFolderObj->Add(hs);
1646
1647   //
1648   // trigger efficiency for INEL
1649   //
1650   h = AlidNdPtHelper::GenerateCorrMatrix(fTriggerEventMatrix->Projection(0),fGenEventMatrix->Projection(0),"zv_trig_INEL_eff_matrix");
1651   aFolderObj->Add(h);
1652
1653
1654   //
1655   // trigger bias correction (MB to INEL)
1656   //
1657   hs = AlidNdPtHelper::GenerateCorrMatrix(fGenEventMatrix,fTriggerEventMatrix,"zv_mult_trig_MBtoInel_corr_matrix");
1658   aFolderObj->Add(hs);
1659
1660   h = AlidNdPtHelper::GenerateCorrMatrix(fGenEventMatrix->Projection(0),fTriggerEventMatrix->Projection(0),"zv_trig_MBtoInel_corr_matrix");
1661   aFolderObj->Add(h);
1662
1663   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenEventMatrix->Projection(0,1),fTriggerEventMatrix->Projection(0,1),"zv_mult_trig_MBtoInel_corr_matrix_2D");
1664   aFolderObj->Add(h2D);
1665
1666
1667   h = AlidNdPtHelper::GenerateCorrMatrix(fGenEventMatrix->Projection(1),fTriggerEventMatrix->Projection(1),"mult_trig_MBtoInel_corr_matrix");
1668   aFolderObj->Add(h);
1669
1670
1671   //
1672   // event vertex reconstruction correction (MB)
1673   //
1674   hs = AlidNdPtHelper::GenerateCorrMatrix(fTriggerEventMatrix,fRecEventMatrix,"zv_mult_event_corr_matrix");
1675   aFolderObj->Add(hs);
1676
1677   h2D = AlidNdPtHelper::GenerateCorrMatrix(fTriggerEventMatrix->Projection(0,1),fRecEventMatrix->Projection(0,1),"zv_mult_event_corr_matrix_2D");
1678   aFolderObj->Add(h2D);
1679
1680
1681   h = AlidNdPtHelper::GenerateCorrMatrix(fTriggerEventMatrix->Projection(1),fRecEventMatrix->Projection(1),"mult_event_corr_matrix");
1682   aFolderObj->Add(h);
1683
1684
1685   h = AlidNdPtHelper::GenerateCorrMatrix(fTriggerEventMatrix->Projection(0),fRecEventMatrix->Projection(0),"zv_event_corr_matrix");
1686   aFolderObj->Add(h);
1687
1688
1689
1690   //
1691   // track-event trigger bias correction (MB to INEL)
1692   //
1693   hs = AlidNdPtHelper::GenerateCorrMatrix(fGenTrackEventMatrix,fTriggerTrackEventMatrix,"zv_pt_eta_track_trig_MBtoInel_corr_matrix");
1694   aFolderObj->Add(hs);
1695
1696   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenTrackEventMatrix->Projection(1,2),fTriggerTrackEventMatrix->Projection(1,2),"eta_pt_track_trig_MBtoInel_corr_matrix");
1697   aFolderObj->Add(h2D);
1698
1699   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenTrackEventMatrix->Projection(1,0),fTriggerTrackEventMatrix->Projection(1,0),"pt_zv_track_trig_MBtoInel_corr_matrix");
1700   aFolderObj->Add(h2D);
1701
1702   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenTrackEventMatrix->Projection(2,0),fTriggerTrackEventMatrix->Projection(2,0),"zv_eta_track_trig_MBtoInel_corr_matrix");
1703   aFolderObj->Add(h2D);
1704
1705   // efficiency
1706
1707   h = AlidNdPtHelper::GenerateCorrMatrix(fTriggerTrackEventMatrix->Projection(1),fGenTrackEventMatrix->Projection(1),"pt_track_trig_MBtoInel_eff_matrix");
1708   aFolderObj->Add(h);
1709
1710
1711   //
1712   // track-event vertex reconstruction correction (MB)
1713   //
1714   hs = AlidNdPtHelper::GenerateCorrMatrix(fTriggerTrackEventMatrix,fRecTrackEventMatrix,"zv_pt_eta_track_event_corr_matrix");
1715   aFolderObj->Add(hs);
1716
1717   h2D = AlidNdPtHelper::GenerateCorrMatrix(fTriggerTrackEventMatrix->Projection(1,2),fRecTrackEventMatrix->Projection(1,2),"eta_pt_track_event_corr_matrix");
1718   aFolderObj->Add(h2D);
1719
1720   h2D = AlidNdPtHelper::GenerateCorrMatrix(fTriggerTrackEventMatrix->Projection(1,0),fRecTrackEventMatrix->Projection(1,0),"pt_zv_track_event_corr_matrix");
1721   aFolderObj->Add(h2D);
1722
1723   h2D = AlidNdPtHelper::GenerateCorrMatrix(fTriggerTrackEventMatrix->Projection(2,0),fRecTrackEventMatrix->Projection(2,0),"zv_eta_track_event_corr_matrix");
1724   aFolderObj->Add(h2D);
1725
1726   // efficiency
1727
1728   h = AlidNdPtHelper::GenerateCorrMatrix(fRecTrackEventMatrix->Projection(1),fTriggerTrackEventMatrix->Projection(1),"pt_track_event_eff_matrix");
1729   aFolderObj->Add(h);
1730
1731
1732   //
1733   // track rec. efficiency correction
1734   //
1735   hs = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix,fRecPrimTrackMatrix,"zv_pt_eta_track_corr_matrix");
1736   aFolderObj->Add(hs);
1737
1738   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(1,2),fRecPrimTrackMatrix->Projection(1,2),"eta_pt_track_corr_matrix");
1739   aFolderObj->Add(h2D);
1740
1741   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(1,0),fRecPrimTrackMatrix->Projection(1,0),"pt_zv_track_corr_matrix");
1742   aFolderObj->Add(h2D);
1743
1744   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(2,0),fRecPrimTrackMatrix->Projection(2,0),"zv_eta_track_corr_matrix");
1745   aFolderObj->Add(h2D);
1746
1747
1748   h = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(0),fRecPrimTrackMatrix->Projection(0),"zv_track_corr_matrix");
1749   aFolderObj->Add(h);
1750
1751   h = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(1),fRecPrimTrackMatrix->Projection(1),"pt_track_corr_matrix");
1752   aFolderObj->Add(h);
1753
1754   // efficiency
1755
1756   h = AlidNdPtHelper::GenerateCorrMatrix(fRecPrimTrackMatrix->Projection(1), fGenPrimTrackMatrix->Projection(1),"pt_track_eff_matrix");
1757   aFolderObj->Add(h);
1758
1759   h = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(2),fRecPrimTrackMatrix->Projection(2),"eta_track_corr_matrix");
1760   aFolderObj->Add(h);
1761
1762   //
1763   // secondary track contamination correction
1764   //
1765   //hs = AlidNdPtHelper::GenerateContCorrMatrix(fRecSecTrackMatrix,fRecTrackMatrix,"zv_pt_eta_track_cont_matrix");
1766   hs = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix,fRecTrackMatrix,"zv_pt_eta_track_cont_matrix");
1767   aFolderObj->Add(hs);
1768
1769   h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(1,2),fRecTrackMatrix->Projection(1,2),"eta_pt_track_cont_matrix");
1770   aFolderObj->Add(h2D);
1771
1772   h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(1,0),fRecTrackMatrix->Projection(1,0),"pt_zv_track_cont_matrix");
1773   aFolderObj->Add(h2D);
1774
1775   h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(2,0),fRecTrackMatrix->Projection(2,0),"zv_eta_track_cont_matrix");
1776   aFolderObj->Add(h2D);
1777
1778   h = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(0),fRecTrackMatrix->Projection(0),"zv_track_cont_matrix");
1779   aFolderObj->Add(h);
1780
1781
1782   h = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(1),fRecTrackMatrix->Projection(1),"pt_track_cont_matrix");
1783   aFolderObj->Add(h);
1784
1785   h = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(2),fRecTrackMatrix->Projection(2),"eta_track_cont_matrix");
1786   aFolderObj->Add(h);
1787
1788   //
1789   // multiple track reconstruction correction
1790   //
1791   //hs = AlidNdPtHelper::GenerateContCorrMatrix(fRecMultTrackMatrix,fRecTrackMatrix,"zv_pt_eta_mult_track_cont_matrix");
1792   hs = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix,fRecTrackMatrix,"zv_pt_eta_mult_track_cont_matrix");
1793   aFolderObj->Add(hs);
1794
1795   h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(1,2),fRecTrackMatrix->Projection(1,2),"eta_pt_mult_track_cont_matrix");
1796   aFolderObj->Add(h2D);
1797
1798   h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(1,0),fRecTrackMatrix->Projection(1,0),"pt_zv_mult_track_cont_matrix");
1799   aFolderObj->Add(h2D);
1800
1801   h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(2,0),fRecTrackMatrix->Projection(2,0),"zv_eta_mult_track_cont_matrix");
1802   aFolderObj->Add(h2D);
1803
1804   h = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(0),fRecTrackMatrix->Projection(0),"zv_mult_track_cont_matrix");
1805   aFolderObj->Add(h);
1806
1807   h = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(1),fRecTrackMatrix->Projection(1),"pt_mult_track_cont_matrix");
1808   aFolderObj->Add(h);
1809
1810   h = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(2),fRecTrackMatrix->Projection(2),"eta_mult_track_cont_matrix");
1811   aFolderObj->Add(h);
1812
1813   //
1814   // Control histograms
1815   //
1816
1817   if(fHistogramsOn) {
1818
1819   // Efficiency electrons, muons, pions, kaons, protons, all
1820   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(1,1);
1821   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(1,1);
1822   h1 = fMCPrimTrackHist1[1]->Projection(0);
1823   h2 = fMCPrimTrackHist1[2]->Projection(0);
1824   h2c = (TH1D *)h2->Clone();
1825   h2c->Divide(h1);
1826   h2c->SetName("eff_pt_electrons");
1827   aFolderObj->Add(h2c);
1828
1829   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(2,2);
1830   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(2,2);
1831   h1 = fMCPrimTrackHist1[1]->Projection(0);
1832   h2 = fMCPrimTrackHist1[2]->Projection(0);
1833   h2c = (TH1D *)h2->Clone();
1834   h2c->Divide(h1);
1835   h2c->SetName("eff_pt_muons");
1836   aFolderObj->Add(h2c);
1837
1838   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(3,3);
1839   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(3,3);
1840   h1 = fMCPrimTrackHist1[1]->Projection(0);
1841   h2 = fMCPrimTrackHist1[2]->Projection(0);
1842   h2c = (TH1D *)h2->Clone();
1843   h2c->Divide(h1);
1844   h2c->SetName("eff_pt_pions");
1845   aFolderObj->Add(h2c);
1846
1847   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(4,4);
1848   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(4,4);
1849   h1 = fMCPrimTrackHist1[1]->Projection(0);
1850   h2 = fMCPrimTrackHist1[2]->Projection(0);
1851   h2c = (TH1D *)h2->Clone();
1852   h2c->Divide(h1);
1853   h2c->SetName("eff_pt_kaons");
1854   aFolderObj->Add(h2c);
1855
1856   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(5,5);
1857   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(5,5);
1858   h1 = fMCPrimTrackHist1[1]->Projection(0);
1859   h2 = fMCPrimTrackHist1[2]->Projection(0);
1860   h2c = (TH1D *)h2->Clone();
1861   h2c->Divide(h1);
1862   h2c->SetName("eff_pt_protons");
1863   aFolderObj->Add(h2c);
1864
1865   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(1,5);
1866   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(1,5);
1867   h1 = fMCPrimTrackHist1[1]->Projection(0);
1868   h2 = fMCPrimTrackHist1[2]->Projection(0);
1869   h2c = (TH1D *)h2->Clone();
1870   h2c->Divide(h1);
1871   h2c->SetName("eff_pt_selected");
1872   aFolderObj->Add(h2c);
1873
1874   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(1,6);
1875   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(1,6);
1876   h1 = fMCPrimTrackHist1[1]->Projection(0);
1877   h2 = fMCPrimTrackHist1[2]->Projection(0);
1878   h2c = (TH1D *)h2->Clone();
1879   h2c->Divide(h1);
1880   h2c->SetName("eff_pt_all");
1881   aFolderObj->Add(h2c);
1882
1883   fMCPrimTrackHist1[1]->GetAxis(1)->SetRange(1,fMCPrimTrackHist1[1]->GetAxis(1)->GetNbins());
1884   fMCPrimTrackHist1[2]->GetAxis(1)->SetRange(1,fMCPrimTrackHist1[2]->GetAxis(1)->GetNbins());
1885
1886   //  pt spetra
1887   // - rec, primaries, secondaries
1888   // - primaries (pid)
1889   // - secondaries (pid)
1890   // - secondaries (mech)
1891   // - secondaries (mother)
1892   //
1893
1894   TH1D *mcPtAccall = fMCTrackHist1[1]->Projection(0);
1895   mcPtAccall->SetName("mc_pt_acc_all");
1896   aFolderObj->Add(mcPtAccall);
1897
1898   TH1D *mcPtAccprim = fMCPrimTrackHist1[1]->Projection(0);
1899   mcPtAccprim->SetName("mc_pt_acc_prim");
1900   aFolderObj->Add(mcPtAccprim);
1901
1902   TH1D *mcPtRecall = fMCTrackHist1[2]->Projection(0);
1903   mcPtRecall->SetName("mc_pt_rec_all");
1904   aFolderObj->Add(mcPtRecall);
1905
1906   TH1D *mcPtRecprim = fMCPrimTrackHist1[2]->Projection(0);
1907   mcPtRecprim->SetName("mc_pt_rec_prim");
1908   aFolderObj->Add(mcPtRecprim);
1909
1910   TH1D *mcPtRecsec = fMCSecTrackHist1[2]->Projection(0);
1911   mcPtRecsec->SetName("mc_pt_rec_sec");
1912   aFolderObj->Add(mcPtRecsec);
1913
1914   for(Int_t i = 0; i<6; i++)
1915   {
1916     sprintf(name,"mc_pt_rec_prim_pid_%d",i);
1917     fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(i+1,i+1);
1918     h = fMCPrimTrackHist1[2]->Projection(0);
1919     h->SetName(name);
1920     aFolderObj->Add(h);
1921
1922     sprintf(name,"mc_pt_rec_sec_pid_%d",i);
1923     fMCSecTrackHist1[2]->GetAxis(2)->SetRange(i+1,i+1);
1924     h = fMCSecTrackHist1[2]->Projection(0);
1925     h->SetName(name);
1926     aFolderObj->Add(h);
1927
1928     // production mechanisms for given pid
1929     fMCSecTrackHist1[2]->GetAxis(2)->SetRange(i+1,i+1);
1930
1931     for(Int_t j=0; j<20; j++) {
1932       if(j == 4) {
1933         // decay
1934
1935         sprintf(name,"mc_pt_rec_sec_pid_%d_decay",i);
1936         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1937         h = fMCSecTrackHist1[2]->Projection(0);
1938         h->SetName(name);
1939         aFolderObj->Add(h);
1940
1941         sprintf(name,"mc_eta_rec_sec_pid_%d_decay",i);
1942         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1943         h = fMCSecTrackHist1[2]->Projection(1);
1944         h->SetName(name);
1945         aFolderObj->Add(h);
1946
1947         sprintf(name,"mc_mother_rec_sec_pid_%d_decay",i);
1948         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1949         h = fMCSecTrackHist1[2]->Projection(4);
1950         h->SetName(name);
1951         aFolderObj->Add(h);
1952
1953       } else if (j == 5) {
1954        // conversion
1955
1956         sprintf(name,"mc_pt_rec_sec_pid_%d_conv",i);
1957         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1958         h = fMCSecTrackHist1[2]->Projection(0);
1959         h->SetName(name);
1960         aFolderObj->Add(h);
1961
1962         sprintf(name,"mc_eta_rec_sec_pid_%d_conv",i);
1963         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1964         h = fMCSecTrackHist1[2]->Projection(1);
1965         h->SetName(name);
1966         aFolderObj->Add(h);
1967
1968         sprintf(name,"mc_mother_rec_sec_pid_%d_conv",i);
1969         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1970         h = fMCSecTrackHist1[2]->Projection(4);
1971         h->SetName(name);
1972         aFolderObj->Add(h);
1973
1974      } else if (j == 13) {
1975        // mat
1976
1977         sprintf(name,"mc_pt_rec_sec_pid_%d_mat",i);
1978         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1979         h = fMCSecTrackHist1[2]->Projection(0);
1980         h->SetName(name);
1981         aFolderObj->Add(h);
1982
1983         sprintf(name,"mc_eta_rec_sec_pid_%d_mat",i);
1984         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1985         h = fMCSecTrackHist1[2]->Projection(1);
1986         h->SetName(name);
1987         aFolderObj->Add(h);
1988
1989         sprintf(name,"mc_eta_mother_rec_sec_pid_%d_mat",i);
1990         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1991         h = fMCSecTrackHist1[2]->Projection(4,1);
1992         h->SetName(name);
1993         aFolderObj->Add(h);
1994
1995         sprintf(name,"mc_mother_rec_sec_pid_%d_mat",i);
1996         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1997         h = fMCSecTrackHist1[2]->Projection(4);
1998         h->SetName(name);
1999         aFolderObj->Add(h);
2000
2001         sprintf(name,"mc_pt_mother_rec_sec_pid_%d_mat",i);
2002         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
2003         h = fMCSecTrackHist1[2]->Projection(4,0);
2004         h->SetName(name);
2005         aFolderObj->Add(h);
2006
2007      } else {
2008        continue;
2009      }
2010    }
2011
2012   }
2013   } // end fHistogramOn
2014
2015   //
2016   //  resolution histograms
2017   //  only for reconstructed tracks
2018   //
2019
2020   TH2F *h2F=0;
2021   TCanvas * c = new TCanvas("resol","resol");
2022   c->cd();
2023
2024   //
2025   fRecMCTrackHist1->GetAxis(1)->SetRangeUser(-0.8,0.79);
2026
2027   h2F = (TH2F*)fRecMCTrackHist1->Projection(2,0);
2028   h = AlidNdPtHelper::MakeResol(h2F,1,0,kTRUE,10);
2029   h->SetXTitle("p_{tmc} (GeV/c)");
2030   h->SetYTitle("(p_{t}-p_{tmc})/p_{tmc} resolution");
2031   h->Draw();
2032   h->SetName("pt_resolution_vs_mcpt");
2033   aFolderObj->Add(h);
2034
2035   h2F = (TH2F*)fRecMCTrackHist1->Projection(2,0);
2036   h = AlidNdPtHelper::MakeResol(h2F,1,1,kTRUE,10);
2037   h->SetXTitle("p_{tmc} (GeV/c)");
2038   h->SetYTitle("(p_{t}-p_{tmc})/p_{tmc} mean");
2039   h->Draw();
2040   h->SetName("dpt_mean_vs_mcpt");
2041   aFolderObj->Add(h);
2042
2043   //
2044   h2F = (TH2F*)fRecMCTrackHist1->Projection(3,0);
2045   h = AlidNdPtHelper::MakeResol(h2F,1,0,kTRUE,10);
2046   h->SetXTitle("p_{tmc} (GeV/c)");
2047   h->SetYTitle("(#eta-#eta_{mc}) resolution");
2048   h->Draw();
2049   h->SetName("eta_resolution_vs_mcpt");
2050   aFolderObj->Add(h);
2051
2052   h2F = (TH2F*)fRecMCTrackHist1->Projection(3,0);
2053   h = AlidNdPtHelper::MakeResol(h2F,1,1,kTRUE,10);
2054   h->SetXTitle("p_{tmc} (GeV/c)");
2055   h->SetYTitle("(#eta-mc#eta) mean");
2056   h->Draw();
2057   h->SetName("deta_mean_vs_mcpt");
2058   aFolderObj->Add(h);
2059
2060   //
2061   fRecMCTrackHist1->GetAxis(1)->SetRange(1,fRecMCTrackHist1->GetAxis(1)->GetNbins());
2062
2063   h2F = (TH2F*)fRecMCTrackHist1->Projection(2,1);
2064   h = AlidNdPtHelper::MakeResol(h2F,1,0,kTRUE,10);
2065   h->SetXTitle("#eta_{mc}");
2066   h->SetYTitle("(p_{t}-p_{tmc})/p_{tmc} resolution");
2067   h->Draw();
2068   h->SetName("pt_resolution_vs_mceta");
2069   aFolderObj->Add(h);
2070
2071   h2F = (TH2F*)fRecMCTrackHist1->Projection(2,1);
2072   h = AlidNdPtHelper::MakeResol(h2F,1,1,kTRUE,10);
2073   h->SetXTitle("#eta_{mc}");
2074   h->SetYTitle("(p_{t}-p_{tmc})/p_{tmc} mean");
2075   h->Draw();
2076   h->SetName("dpt_mean_vs_mceta");
2077   aFolderObj->Add(h);
2078
2079   //
2080   h2F = (TH2F*)fRecMCTrackHist1->Projection(3,1);
2081   h = AlidNdPtHelper::MakeResol(h2F,1,0,kTRUE,10);
2082   h->SetXTitle("#eta_{mc}");
2083   h->SetYTitle("(#eta-#eta_{mc}) resolution");
2084   h->Draw();
2085   h->SetName("eta_resolution_vs_mceta");
2086   aFolderObj->Add(h);
2087
2088   h2F = (TH2F*)fRecMCTrackHist1->Projection(3,1);
2089   h = AlidNdPtHelper::MakeResol(h2F,1,1,kTRUE,10);
2090   h->SetXTitle("#eta_{mc}");
2091   h->SetYTitle("(#eta-mc#eta) mean");
2092   h->Draw();
2093   h->SetName("deta_mean_vs_mceta");
2094   aFolderObj->Add(h);
2095
2096   fRecMCTrackHist1->GetAxis(0)->SetRange(1,fRecMCTrackHist1->GetAxis(0)->GetNbins());
2097
2098   } // end use MC info
2099
2100   // export objects to analysis folder
2101   fAnalysisFolder = ExportToFolder(aFolderObj);
2102
2103   // delete only TObjArray
2104   if(aFolderObj) delete aFolderObj;
2105 }
2106
2107 //_____________________________________________________________________________
2108 TFolder* AlidNdPtAnalysisPbPb::ExportToFolder(TObjArray * const array)
2109 {
2110   // recreate folder avery time and export objects to new one
2111   //
2112   AlidNdPtAnalysisPbPb * comp=this;
2113   TFolder *folder = comp->GetAnalysisFolder();
2114
2115   TString name, title;
2116   TFolder *newFolder = 0;
2117   Int_t i = 0;
2118   Int_t size = array->GetSize();
2119
2120   if(folder) {
2121      // get name and title from old folder
2122      name = folder->GetName();
2123      title = folder->GetTitle();
2124
2125          // delete old one
2126      delete folder;
2127
2128          // create new one
2129      newFolder = CreateFolder(name.Data(),title.Data());
2130      newFolder->SetOwner();
2131
2132          // add objects to folder
2133      while(i < size) {
2134            newFolder->Add(array->At(i));
2135            i++;
2136          }
2137   }
2138
2139 return newFolder;
2140 }
2141
2142 //_____________________________________________________________________________
2143 TFolder* AlidNdPtAnalysisPbPb::CreateFolder(TString name,TString title) {
2144 // create folder for analysed histograms
2145 //
2146 TFolder *folder = 0;
2147   folder = new TFolder(name.Data(),title.Data());
2148
2149   return folder;
2150 }