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