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   TParticle* particle = stack->Particle(label);
1157   if(!particle) return;
1158
1159   Int_t motherPdg = -1;
1160   TParticle* mother = 0;
1161
1162   //TParticle* prim_mother = AlidNdPtHelper::FindPrimaryMother(stack,label);
1163   Int_t motherLabel = particle->GetMother(0);
1164   if(motherLabel>0) mother = stack->Particle(motherLabel);
1165   if(mother) motherPdg = TMath::Abs(mother->GetPdgCode()); // take abs for visualisation only
1166   //Int_t mech = particle->GetUniqueID(); // production mechanism
1167
1168   if(!particle->GetPDG()) return;
1169   Double_t gq = particle->GetPDG()->Charge()/3.0; // Charge units |e|/3
1170   if(TMath::Abs(gq)<0.001) return;
1171   Float_t gpt = particle->Pt();
1172   Float_t geta = particle->Eta();
1173
1174   Double_t dpt=0;
1175   //printf("pt %f, gpt %f \n",pt,gpt);
1176   if(gpt) dpt = (pt-gpt)/gpt;
1177   Double_t deta = (eta-geta);
1178
1179   // fill histograms
1180   if(trackObj == AlidNdPtHelper::kRecTracks)  //RecTracks???
1181   {
1182     Double_t vTrackPtCorrelationMatrix[4]={pt,gpt,geta,centralityF};
1183     fTrackPtCorrelationMatrix->Fill(vTrackPtCorrelationMatrix);
1184
1185     Double_t vRecMCTrackHist1[5]={gpt,geta,dpt,deta,centralityF};
1186     fRecMCTrackHist1->Fill(vRecMCTrackHist1);
1187   }
1188 }
1189
1190 //_____________________________________________________________________________
1191 void AlidNdPtAnalysisPbPb::FillHistograms(AliStack *const stack, Int_t label, AlidNdPtHelper::TrackObject trackObj, Float_t centralityF)
1192 {
1193
1194   // Fill MC histograms
1195   if(!stack) return;
1196
1197   TParticle* particle = stack->Particle(label);
1198   if(!particle) return;
1199
1200   Int_t motherPdg = -1;
1201   TParticle* mother = 0;
1202
1203   //TParticle* prim_mother = AlidNdPtHelper::FindPrimaryMother(stack,label);
1204   Int_t motherLabel = particle->GetMother(0);
1205   if(motherLabel>0) mother = stack->Particle(motherLabel);
1206   if(mother) motherPdg = TMath::Abs(mother->GetPdgCode()); // take abs for visualisation only
1207   Int_t mech = particle->GetUniqueID(); // production mechanism
1208
1209   if(!particle->GetPDG()) return;
1210   Double_t gq = particle->GetPDG()->Charge()/3.0; // Charge units |e|/3
1211   if(TMath::Abs(gq) < 0.001) return;
1212
1213   Float_t gpt = particle->Pt();
1214   //Float_t qgpt = particle->Pt() * gq;
1215   Float_t geta = particle->Eta();
1216   Float_t gphi = particle->Phi();
1217   //Float_t gpz = particle->Pz();
1218
1219   Bool_t prim = stack->IsPhysicalPrimary(label);
1220   //Float_t vx = particle->Vx(); Float_t vy = particle->Vy(); Float_t vz = particle->Vz();
1221
1222   Int_t pid = AlidNdPtHelper::ConvertPdgToPid(particle);
1223
1224   //if(prim&&pid==5) printf("pdgcode %d, production mech %d \n",particle->GetPdgCode(),mech);
1225   //if(!prim) printf("motherPdg %d, particle %d, production mech %d\n",motherPdg, particle->GetPdgCode(),mech);
1226
1227   //
1228   // fill histogram
1229   //
1230   Double_t vMCTrackHist1[4] = {gpt,geta,gphi,centralityF};
1231   fMCTrackHist1[trackObj]->Fill(vMCTrackHist1);
1232
1233   Double_t vMCPrimTrackHist1[6] = {gpt,geta,pid,mech,motherPdg,centralityF};
1234   Double_t vMCPrimTrackHist2[4] = {TMath::Abs(particle->GetPdgCode()),mech,motherPdg,centralityF};
1235
1236   //if(prim && AliPWG0Helper::IsPrimaryCharged(particle, label)) fMCPrimTrackHist1[trackObj]->Fill(vMCPrimTrackHist1);
1237
1238   if(prim) {
1239     fMCPrimTrackHist1[trackObj]->Fill(vMCPrimTrackHist1);
1240     if(pid == 5) fMCPrimTrackHist2[trackObj]->Fill(vMCPrimTrackHist2);
1241   }
1242   else {
1243     fMCSecTrackHist1[trackObj]->Fill(vMCPrimTrackHist1);
1244   }
1245
1246 }
1247
1248 //_____________________________________________________________________________
1249 Long64_t AlidNdPtAnalysisPbPb::Merge(TCollection* const list)
1250 {
1251   // Merge list of objects (needed by PROOF)
1252
1253   //  init if not done already
1254   if (!fIsInit) { Init(); }
1255
1256   if (!list)
1257   return 0;
1258
1259   if (list->IsEmpty())
1260   return 1;
1261
1262   TIterator* iter = list->MakeIterator();
1263   TObject* obj = 0;
1264
1265   //
1266   //TList *collPhysSelection = new TList;
1267
1268   // collection of generated histograms
1269
1270   Int_t count=0;
1271   while((obj = iter->Next()) != 0) {
1272     AlidNdPtAnalysisPbPb* entry = dynamic_cast<AlidNdPtAnalysisPbPb*>(obj);
1273     if (entry == 0) continue;
1274
1275     // physics selection
1276     //printf("entry->GetPhysicsTriggerSelection() %p \n", entry->GetPhysicsTriggerSelection());
1277     //AliPhysicsSelection *physSel = entry->GetPhysicsTriggerSelection();
1278     //if( physSel ) collPhysSelection->Add(physSel);
1279
1280     //
1281     fTrackPtCorrelationMatrix->Add(entry->fTrackPtCorrelationMatrix);
1282
1283     //
1284     fGenEventMatrix->Add(entry->fGenEventMatrix);
1285
1286     fTriggerEventMatrix->Add(entry->fTriggerEventMatrix);
1287
1288     fRecEventMatrix->Add(entry->fRecEventMatrix);
1289     //
1290     fGenTrackEventMatrix->Add(entry->fGenTrackEventMatrix);
1291
1292     fTriggerTrackEventMatrix->Add(entry->fTriggerTrackEventMatrix);
1293
1294     fRecTrackEventMatrix->Add(entry->fRecTrackEventMatrix);
1295
1296     //
1297     fGenTrackMatrix->Add(entry->fGenTrackMatrix);
1298     fGenPrimTrackMatrix->Add(entry->fGenPrimTrackMatrix);
1299     fRecPrimTrackMatrix->Add(entry->fRecPrimTrackMatrix);
1300     //
1301     fRecTrackMatrix->Add(entry->fRecTrackMatrix);
1302     fRecSecTrackMatrix->Add(entry->fRecSecTrackMatrix);
1303     //
1304     fRecMultTrackMatrix->Add(entry->fRecMultTrackMatrix);
1305
1306     //
1307     // control analysis histograms
1308     //
1309     fMCEventHist1->Add(entry->fMCEventHist1);
1310     fRecEventHist1->Add(entry->fRecEventHist1);
1311     fRecEventHist2->Add(entry->fRecEventHist2);
1312     fRecMCEventHist1->Add(entry->fRecMCEventHist1);
1313     fRecMCEventHist2->Add(entry->fRecMCEventHist2);
1314
1315
1316     for(Int_t i=0; i<AlidNdPtHelper::kCutSteps; i++) {
1317       fMCTrackHist1[i]->Add(entry->fMCTrackHist1[i]);
1318
1319       fMCPrimTrackHist1[i]->Add(entry->fMCPrimTrackHist1[i]);
1320       fMCPrimTrackHist2[i]->Add(entry->fMCPrimTrackHist2[i]);
1321       fMCSecTrackHist1[i]->Add(entry->fMCSecTrackHist1[i]);
1322
1323       fRecTrackHist1[i]->Add(entry->fRecTrackHist1[i]);
1324       fRecTrackHist2[i]->Add(entry->fRecTrackHist2[i]);
1325       fRecTrackMultHist1[i]->Add(entry->fRecTrackMultHist1[i]);
1326     }
1327     fRecMCTrackHist1->Add(entry->fRecMCTrackHist1);
1328     fMCMultRecTrackHist1->Add(entry->fMCMultRecTrackHist1);
1329     fRecTrackHist3->Add(entry->fRecTrackHist3);
1330
1331   count++;
1332   }
1333
1334   //AliPhysicsSelection *trigSelection = GetPhysicsTriggerSelection();
1335   //if( trigSelection ) trigSelection->Merge(collPhysSelection);
1336   //if(collPhysSelection) delete collPhysSelection;
1337
1338 return count;
1339 }
1340
1341
1342
1343 //_____________________________________________________________________________
1344 void AlidNdPtAnalysisPbPb::Analyse()
1345 {
1346
1347   //  init if not done already
1348   if (!fIsInit) { Init(); }
1349
1350   // Analyse histograms
1351   //
1352   TH1::AddDirectory(kFALSE);
1353   TH1 *h=0, *h1=0, *h2=0, *h2c = 0;
1354   THnSparse *hs=0;
1355   TH2 *h2D=0;
1356
1357   char name[256];
1358   TObjArray *aFolderObj = new TObjArray;
1359
1360   //
1361   // LHC backgraund in all and 0-bins
1362   //
1363   //AliPhysicsSelection *trigSel = GetPhysicsTriggerSelection();
1364   //trigSel->SaveHistograms("physics_selection");
1365
1366   //
1367   // Reconstructed event vertex
1368   //
1369   h = fRecEventHist1->Projection(0);
1370   h->SetName("Xv");
1371   aFolderObj->Add(h);
1372
1373   h = fRecEventHist1->Projection(1);
1374   h->SetName("Yv");
1375   aFolderObj->Add(h);
1376
1377   h = fRecEventHist1->Projection(2);
1378   h->SetName("Zv");
1379   aFolderObj->Add(h);
1380
1381   //
1382   // multiplicity
1383   //
1384   h = fRecEventHist2->Projection(1);
1385   h->SetName("multMB");
1386   aFolderObj->Add(h);
1387
1388   h2D = fRecEventHist2->Projection(0,1);
1389   h2D->SetName("Zv_vs_multiplicity_MB");
1390   aFolderObj->Add(h2D);
1391
1392   //
1393   // reconstructed pt histograms
1394   //
1395   h = fRecTrackHist1[0]->Projection(0);
1396   h->Scale(1.,"width");
1397   h->SetName("pt_all_ch");
1398   aFolderObj->Add(h);
1399
1400   h = fRecTrackHist1[1]->Projection(0);
1401   h->Scale(1.,"width");
1402   h->SetName("pt_acc");
1403   aFolderObj->Add(h);
1404
1405   h = fRecTrackHist1[2]->Projection(0);
1406   h->Scale(1.,"width");
1407   h->SetName("pt_rec");
1408   aFolderObj->Add(h);
1409
1410   //
1411   // reconstructed eta histograms
1412   //
1413   h = fRecTrackHist1[0]->Projection(1);
1414   h->SetName("eta_all_ch");
1415   aFolderObj->Add(h);
1416
1417   h = fRecTrackHist1[1]->Projection(1);
1418   h->SetName("eta_acc");
1419   aFolderObj->Add(h);
1420
1421   h = fRecTrackHist1[2]->Projection(1);
1422   h->SetName("eta_rec");
1423   aFolderObj->Add(h);
1424
1425   //
1426   // reconstructed phi histograms
1427   //
1428   h = fRecTrackHist1[0]->Projection(2);
1429   h->SetName("phi_all_ch");
1430   aFolderObj->Add(h);
1431
1432   h = fRecTrackHist1[1]->Projection(2);
1433   h->SetName("phi_acc");
1434   aFolderObj->Add(h);
1435
1436   h = fRecTrackHist1[2]->Projection(2);
1437   h->SetName("phi_rec");
1438   aFolderObj->Add(h);
1439
1440   //
1441   // reconstructed eta:pt histograms
1442   //
1443   h2D = fRecTrackHist1[0]->Projection(1,0);
1444   h2D->SetName("pt_eta_all_ch");
1445   aFolderObj->Add(h2D);
1446
1447   h2D = fRecTrackHist1[1]->Projection(1,0);
1448   h2D->SetName("pt_eta_acc");
1449   aFolderObj->Add(h2D);
1450
1451   h2D = fRecTrackHist1[2]->Projection(1,0);
1452   h2D->SetName("pt_eta_rec");
1453   aFolderObj->Add(h2D);
1454
1455   //
1456   // reconstructed phi:pt histograms
1457   //
1458   h2D = fRecTrackHist1[0]->Projection(2,0);
1459   h2D->SetName("pt_phi_all_ch");
1460   aFolderObj->Add(h2D);
1461
1462   h2D = fRecTrackHist1[1]->Projection(2,0);
1463   h2D->SetName("pt_phi_acc");
1464   aFolderObj->Add(h2D);
1465
1466   h2D = fRecTrackHist1[2]->Projection(2,0);
1467   h2D->SetName("pt_phi_rec");
1468   aFolderObj->Add(h2D);
1469
1470   //
1471   // reconstructed phi:eta histograms
1472   //
1473   h2D = fRecTrackHist1[0]->Projection(2,1);
1474   h2D->SetName("eta_phi_all_ch");
1475   aFolderObj->Add(h2D);
1476
1477   h2D = fRecTrackHist1[1]->Projection(2,1);
1478   h2D->SetName("eta_phi_acc");
1479   aFolderObj->Add(h2D);
1480
1481   h2D = fRecTrackHist1[2]->Projection(2,1);
1482   h2D->SetName("eta_phi_rec");
1483   aFolderObj->Add(h2D);
1484
1485   //
1486   // reconstructed nClust, chi2 vs pt, eta, phi
1487   //
1488   if(fHistogramsOn) {
1489
1490     h2D = fRecTrackHist3->Projection(0,1);
1491     h2D->SetName("nClust_chi2_rec");
1492     aFolderObj->Add(h2D);
1493
1494     h2D = fRecTrackHist3->Projection(0,2);
1495     h2D->SetName("nClust_pt_rec");
1496     aFolderObj->Add(h2D);
1497
1498     h2D = fRecTrackHist3->Projection(0,3);
1499     h2D->SetName("nClust_eta_rec");
1500     aFolderObj->Add(h2D);
1501
1502     h2D = fRecTrackHist3->Projection(0,4);
1503     h2D->SetName("nClust_phi_rec");
1504     aFolderObj->Add(h2D);
1505
1506     h2D = fRecTrackHist3->Projection(1,2);
1507     h2D->SetName("chi2_pt_rec");
1508     aFolderObj->Add(h2D);
1509
1510     h2D = fRecTrackHist3->Projection(1,3);
1511     h2D->SetName("chi2_eta_rec");
1512     aFolderObj->Add(h2D);
1513
1514     h2D = fRecTrackHist3->Projection(1,4);
1515     h2D->SetName("chi2_phi_rec");
1516     aFolderObj->Add(h2D);
1517
1518   }
1519
1520   //
1521   // calculate corrections for empty events
1522   // with multMB==0
1523   //
1524
1525   //
1526   // normalised zv to generate zv for triggered events
1527   //
1528   h = fRecEventHist2->Projection(0);
1529   if( h->Integral() ) h->Scale(1./h->Integral());
1530   h->SetName("zv_distribution_norm");
1531   aFolderObj->Add(h);
1532
1533   //
1534   // MC available
1535   //
1536   if(IsUseMCInfo()) {
1537
1538   //
1539   // Event vertex resolution
1540   //
1541   h2D = fRecMCEventHist2->Projection(0,2);
1542   h2D->SetName("DeltaXv_vs_mult");
1543   aFolderObj->Add(h2D);
1544
1545   h2D = fRecMCEventHist2->Projection(1,2);
1546   h2D->SetName("DeltaZv_vs_mult");
1547   aFolderObj->Add(h2D);
1548
1549   //
1550   // normalised zv to get trigger/trigger+vertex event differences
1551   // F(zv) = E_trig(zv,0)/Int(E_trig(zv,0) / Sum(E_trigvtx(zv,n))/Sum(Int(E_trigvtx(zv,n))dzv)
1552   //
1553   fTriggerEventMatrix->GetAxis(1)->SetRangeUser(0.,0.);
1554   h = fTriggerEventMatrix->Projection(0);
1555   h2D = fTriggerEventMatrix->Projection(0,1);
1556   if(h2D->Integral()) h->Scale(1./h2D->Integral());
1557
1558   h1 = fRecEventMatrix->Projection(0);
1559   h2D = fRecEventMatrix->Projection(0,1);
1560   if(h2D->Integral()) h1->Scale(1./h2D->Integral());
1561
1562   h->Divide(h1);
1563   h->SetName("zv_empty_events_norm");
1564   aFolderObj->Add(h);
1565
1566   fTriggerEventMatrix->GetAxis(1)->SetRange(1,fTriggerEventMatrix->GetAxis(1)->GetNbins());
1567
1568   //
1569   // rec. vs true track pt correlation matrix
1570   //
1571   hs = (THnSparse*)fTrackPtCorrelationMatrix->Clone("track_pt_correlation_matrix");
1572   aFolderObj->Add(hs);
1573
1574   //
1575   // trigger efficiency for INEL
1576   //
1577   h = AlidNdPtHelper::GenerateCorrMatrix(fTriggerEventMatrix->Projection(0),fGenEventMatrix->Projection(0),"zv_trig_INEL_eff_matrix");
1578   aFolderObj->Add(h);
1579
1580
1581   //
1582   // trigger bias correction (MB to INEL)
1583   //
1584   hs = AlidNdPtHelper::GenerateCorrMatrix(fGenEventMatrix,fTriggerEventMatrix,"zv_mult_trig_MBtoInel_corr_matrix");
1585   aFolderObj->Add(hs);
1586
1587   h = AlidNdPtHelper::GenerateCorrMatrix(fGenEventMatrix->Projection(0),fTriggerEventMatrix->Projection(0),"zv_trig_MBtoInel_corr_matrix");
1588   aFolderObj->Add(h);
1589
1590   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenEventMatrix->Projection(0,1),fTriggerEventMatrix->Projection(0,1),"zv_mult_trig_MBtoInel_corr_matrix_2D");
1591   aFolderObj->Add(h2D);
1592
1593
1594   h = AlidNdPtHelper::GenerateCorrMatrix(fGenEventMatrix->Projection(1),fTriggerEventMatrix->Projection(1),"mult_trig_MBtoInel_corr_matrix");
1595   aFolderObj->Add(h);
1596
1597
1598   //
1599   // event vertex reconstruction correction (MB)
1600   //
1601   hs = AlidNdPtHelper::GenerateCorrMatrix(fTriggerEventMatrix,fRecEventMatrix,"zv_mult_event_corr_matrix");
1602   aFolderObj->Add(hs);
1603
1604   h2D = AlidNdPtHelper::GenerateCorrMatrix(fTriggerEventMatrix->Projection(0,1),fRecEventMatrix->Projection(0,1),"zv_mult_event_corr_matrix_2D");
1605   aFolderObj->Add(h2D);
1606
1607
1608   h = AlidNdPtHelper::GenerateCorrMatrix(fTriggerEventMatrix->Projection(1),fRecEventMatrix->Projection(1),"mult_event_corr_matrix");
1609   aFolderObj->Add(h);
1610
1611
1612   h = AlidNdPtHelper::GenerateCorrMatrix(fTriggerEventMatrix->Projection(0),fRecEventMatrix->Projection(0),"zv_event_corr_matrix");
1613   aFolderObj->Add(h);
1614
1615
1616
1617   //
1618   // track-event trigger bias correction (MB to INEL)
1619   //
1620   hs = AlidNdPtHelper::GenerateCorrMatrix(fGenTrackEventMatrix,fTriggerTrackEventMatrix,"zv_pt_eta_track_trig_MBtoInel_corr_matrix");
1621   aFolderObj->Add(hs);
1622
1623   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenTrackEventMatrix->Projection(1,2),fTriggerTrackEventMatrix->Projection(1,2),"eta_pt_track_trig_MBtoInel_corr_matrix");
1624   aFolderObj->Add(h2D);
1625
1626   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenTrackEventMatrix->Projection(1,0),fTriggerTrackEventMatrix->Projection(1,0),"pt_zv_track_trig_MBtoInel_corr_matrix");
1627   aFolderObj->Add(h2D);
1628
1629   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenTrackEventMatrix->Projection(2,0),fTriggerTrackEventMatrix->Projection(2,0),"zv_eta_track_trig_MBtoInel_corr_matrix");
1630   aFolderObj->Add(h2D);
1631
1632   // efficiency
1633
1634   h = AlidNdPtHelper::GenerateCorrMatrix(fTriggerTrackEventMatrix->Projection(1),fGenTrackEventMatrix->Projection(1),"pt_track_trig_MBtoInel_eff_matrix");
1635   aFolderObj->Add(h);
1636
1637
1638   //
1639   // track-event vertex reconstruction correction (MB)
1640   //
1641   hs = AlidNdPtHelper::GenerateCorrMatrix(fTriggerTrackEventMatrix,fRecTrackEventMatrix,"zv_pt_eta_track_event_corr_matrix");
1642   aFolderObj->Add(hs);
1643
1644   h2D = AlidNdPtHelper::GenerateCorrMatrix(fTriggerTrackEventMatrix->Projection(1,2),fRecTrackEventMatrix->Projection(1,2),"eta_pt_track_event_corr_matrix");
1645   aFolderObj->Add(h2D);
1646
1647   h2D = AlidNdPtHelper::GenerateCorrMatrix(fTriggerTrackEventMatrix->Projection(1,0),fRecTrackEventMatrix->Projection(1,0),"pt_zv_track_event_corr_matrix");
1648   aFolderObj->Add(h2D);
1649
1650   h2D = AlidNdPtHelper::GenerateCorrMatrix(fTriggerTrackEventMatrix->Projection(2,0),fRecTrackEventMatrix->Projection(2,0),"zv_eta_track_event_corr_matrix");
1651   aFolderObj->Add(h2D);
1652
1653   // efficiency
1654
1655   h = AlidNdPtHelper::GenerateCorrMatrix(fRecTrackEventMatrix->Projection(1),fTriggerTrackEventMatrix->Projection(1),"pt_track_event_eff_matrix");
1656   aFolderObj->Add(h);
1657
1658
1659   //
1660   // track rec. efficiency correction
1661   //
1662   hs = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix,fRecPrimTrackMatrix,"zv_pt_eta_track_corr_matrix");
1663   aFolderObj->Add(hs);
1664
1665   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(1,2),fRecPrimTrackMatrix->Projection(1,2),"eta_pt_track_corr_matrix");
1666   aFolderObj->Add(h2D);
1667
1668   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(1,0),fRecPrimTrackMatrix->Projection(1,0),"pt_zv_track_corr_matrix");
1669   aFolderObj->Add(h2D);
1670
1671   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(2,0),fRecPrimTrackMatrix->Projection(2,0),"zv_eta_track_corr_matrix");
1672   aFolderObj->Add(h2D);
1673
1674
1675   h = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(0),fRecPrimTrackMatrix->Projection(0),"zv_track_corr_matrix");
1676   aFolderObj->Add(h);
1677
1678   h = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(1),fRecPrimTrackMatrix->Projection(1),"pt_track_corr_matrix");
1679   aFolderObj->Add(h);
1680
1681   // efficiency
1682
1683   h = AlidNdPtHelper::GenerateCorrMatrix(fRecPrimTrackMatrix->Projection(1), fGenPrimTrackMatrix->Projection(1),"pt_track_eff_matrix");
1684   aFolderObj->Add(h);
1685
1686   h = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(2),fRecPrimTrackMatrix->Projection(2),"eta_track_corr_matrix");
1687   aFolderObj->Add(h);
1688
1689   //
1690   // secondary track contamination correction
1691   //
1692   //hs = AlidNdPtHelper::GenerateContCorrMatrix(fRecSecTrackMatrix,fRecTrackMatrix,"zv_pt_eta_track_cont_matrix");
1693   hs = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix,fRecTrackMatrix,"zv_pt_eta_track_cont_matrix");
1694   aFolderObj->Add(hs);
1695
1696   h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(1,2),fRecTrackMatrix->Projection(1,2),"eta_pt_track_cont_matrix");
1697   aFolderObj->Add(h2D);
1698
1699   h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(1,0),fRecTrackMatrix->Projection(1,0),"pt_zv_track_cont_matrix");
1700   aFolderObj->Add(h2D);
1701
1702   h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(2,0),fRecTrackMatrix->Projection(2,0),"zv_eta_track_cont_matrix");
1703   aFolderObj->Add(h2D);
1704
1705   h = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(0),fRecTrackMatrix->Projection(0),"zv_track_cont_matrix");
1706   aFolderObj->Add(h);
1707
1708
1709   h = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(1),fRecTrackMatrix->Projection(1),"pt_track_cont_matrix");
1710   aFolderObj->Add(h);
1711
1712   h = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(2),fRecTrackMatrix->Projection(2),"eta_track_cont_matrix");
1713   aFolderObj->Add(h);
1714
1715   //
1716   // multiple track reconstruction correction
1717   //
1718   //hs = AlidNdPtHelper::GenerateContCorrMatrix(fRecMultTrackMatrix,fRecTrackMatrix,"zv_pt_eta_mult_track_cont_matrix");
1719   hs = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix,fRecTrackMatrix,"zv_pt_eta_mult_track_cont_matrix");
1720   aFolderObj->Add(hs);
1721
1722   h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(1,2),fRecTrackMatrix->Projection(1,2),"eta_pt_mult_track_cont_matrix");
1723   aFolderObj->Add(h2D);
1724
1725   h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(1,0),fRecTrackMatrix->Projection(1,0),"pt_zv_mult_track_cont_matrix");
1726   aFolderObj->Add(h2D);
1727
1728   h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(2,0),fRecTrackMatrix->Projection(2,0),"zv_eta_mult_track_cont_matrix");
1729   aFolderObj->Add(h2D);
1730
1731   h = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(0),fRecTrackMatrix->Projection(0),"zv_mult_track_cont_matrix");
1732   aFolderObj->Add(h);
1733
1734   h = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(1),fRecTrackMatrix->Projection(1),"pt_mult_track_cont_matrix");
1735   aFolderObj->Add(h);
1736
1737   h = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(2),fRecTrackMatrix->Projection(2),"eta_mult_track_cont_matrix");
1738   aFolderObj->Add(h);
1739
1740   //
1741   // Control histograms
1742   //
1743
1744   if(fHistogramsOn) {
1745
1746   // Efficiency electrons, muons, pions, kaons, protons, all
1747   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(1,1);
1748   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(1,1);
1749   h1 = fMCPrimTrackHist1[1]->Projection(0);
1750   h2 = fMCPrimTrackHist1[2]->Projection(0);
1751   h2c = (TH1D *)h2->Clone();
1752   h2c->Divide(h1);
1753   h2c->SetName("eff_pt_electrons");
1754   aFolderObj->Add(h2c);
1755
1756   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(2,2);
1757   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(2,2);
1758   h1 = fMCPrimTrackHist1[1]->Projection(0);
1759   h2 = fMCPrimTrackHist1[2]->Projection(0);
1760   h2c = (TH1D *)h2->Clone();
1761   h2c->Divide(h1);
1762   h2c->SetName("eff_pt_muons");
1763   aFolderObj->Add(h2c);
1764
1765   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(3,3);
1766   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(3,3);
1767   h1 = fMCPrimTrackHist1[1]->Projection(0);
1768   h2 = fMCPrimTrackHist1[2]->Projection(0);
1769   h2c = (TH1D *)h2->Clone();
1770   h2c->Divide(h1);
1771   h2c->SetName("eff_pt_pions");
1772   aFolderObj->Add(h2c);
1773
1774   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(4,4);
1775   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(4,4);
1776   h1 = fMCPrimTrackHist1[1]->Projection(0);
1777   h2 = fMCPrimTrackHist1[2]->Projection(0);
1778   h2c = (TH1D *)h2->Clone();
1779   h2c->Divide(h1);
1780   h2c->SetName("eff_pt_kaons");
1781   aFolderObj->Add(h2c);
1782
1783   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(5,5);
1784   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(5,5);
1785   h1 = fMCPrimTrackHist1[1]->Projection(0);
1786   h2 = fMCPrimTrackHist1[2]->Projection(0);
1787   h2c = (TH1D *)h2->Clone();
1788   h2c->Divide(h1);
1789   h2c->SetName("eff_pt_protons");
1790   aFolderObj->Add(h2c);
1791
1792   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(1,5);
1793   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(1,5);
1794   h1 = fMCPrimTrackHist1[1]->Projection(0);
1795   h2 = fMCPrimTrackHist1[2]->Projection(0);
1796   h2c = (TH1D *)h2->Clone();
1797   h2c->Divide(h1);
1798   h2c->SetName("eff_pt_selected");
1799   aFolderObj->Add(h2c);
1800
1801   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(1,6);
1802   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(1,6);
1803   h1 = fMCPrimTrackHist1[1]->Projection(0);
1804   h2 = fMCPrimTrackHist1[2]->Projection(0);
1805   h2c = (TH1D *)h2->Clone();
1806   h2c->Divide(h1);
1807   h2c->SetName("eff_pt_all");
1808   aFolderObj->Add(h2c);
1809
1810   fMCPrimTrackHist1[1]->GetAxis(1)->SetRange(1,fMCPrimTrackHist1[1]->GetAxis(1)->GetNbins());
1811   fMCPrimTrackHist1[2]->GetAxis(1)->SetRange(1,fMCPrimTrackHist1[2]->GetAxis(1)->GetNbins());
1812
1813   //  pt spetra
1814   // - rec, primaries, secondaries
1815   // - primaries (pid)
1816   // - secondaries (pid)
1817   // - secondaries (mech)
1818   // - secondaries (mother)
1819   //
1820
1821   TH1D *mcPtAccall = fMCTrackHist1[1]->Projection(0);
1822   mcPtAccall->SetName("mc_pt_acc_all");
1823   aFolderObj->Add(mcPtAccall);
1824
1825   TH1D *mcPtAccprim = fMCPrimTrackHist1[1]->Projection(0);
1826   mcPtAccprim->SetName("mc_pt_acc_prim");
1827   aFolderObj->Add(mcPtAccprim);
1828
1829   TH1D *mcPtRecall = fMCTrackHist1[2]->Projection(0);
1830   mcPtRecall->SetName("mc_pt_rec_all");
1831   aFolderObj->Add(mcPtRecall);
1832
1833   TH1D *mcPtRecprim = fMCPrimTrackHist1[2]->Projection(0);
1834   mcPtRecprim->SetName("mc_pt_rec_prim");
1835   aFolderObj->Add(mcPtRecprim);
1836
1837   TH1D *mcPtRecsec = fMCSecTrackHist1[2]->Projection(0);
1838   mcPtRecsec->SetName("mc_pt_rec_sec");
1839   aFolderObj->Add(mcPtRecsec);
1840
1841   for(Int_t i = 0; i<6; i++)
1842   {
1843     sprintf(name,"mc_pt_rec_prim_pid_%d",i);
1844     fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(i+1,i+1);
1845     h = fMCPrimTrackHist1[2]->Projection(0);
1846     h->SetName(name);
1847     aFolderObj->Add(h);
1848
1849     sprintf(name,"mc_pt_rec_sec_pid_%d",i);
1850     fMCSecTrackHist1[2]->GetAxis(2)->SetRange(i+1,i+1);
1851     h = fMCSecTrackHist1[2]->Projection(0);
1852     h->SetName(name);
1853     aFolderObj->Add(h);
1854
1855     // production mechanisms for given pid
1856     fMCSecTrackHist1[2]->GetAxis(2)->SetRange(i+1,i+1);
1857
1858     for(Int_t j=0; j<20; j++) {
1859       if(j == 4) {
1860         // decay
1861
1862         sprintf(name,"mc_pt_rec_sec_pid_%d_decay",i);
1863         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1864         h = fMCSecTrackHist1[2]->Projection(0);
1865         h->SetName(name);
1866         aFolderObj->Add(h);
1867
1868         sprintf(name,"mc_eta_rec_sec_pid_%d_decay",i);
1869         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1870         h = fMCSecTrackHist1[2]->Projection(1);
1871         h->SetName(name);
1872         aFolderObj->Add(h);
1873
1874         sprintf(name,"mc_mother_rec_sec_pid_%d_decay",i);
1875         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1876         h = fMCSecTrackHist1[2]->Projection(4);
1877         h->SetName(name);
1878         aFolderObj->Add(h);
1879
1880       } else if (j == 5) {
1881        // conversion
1882
1883         sprintf(name,"mc_pt_rec_sec_pid_%d_conv",i);
1884         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1885         h = fMCSecTrackHist1[2]->Projection(0);
1886         h->SetName(name);
1887         aFolderObj->Add(h);
1888
1889         sprintf(name,"mc_eta_rec_sec_pid_%d_conv",i);
1890         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1891         h = fMCSecTrackHist1[2]->Projection(1);
1892         h->SetName(name);
1893         aFolderObj->Add(h);
1894
1895         sprintf(name,"mc_mother_rec_sec_pid_%d_conv",i);
1896         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1897         h = fMCSecTrackHist1[2]->Projection(4);
1898         h->SetName(name);
1899         aFolderObj->Add(h);
1900
1901      } else if (j == 13) {
1902        // mat
1903
1904         sprintf(name,"mc_pt_rec_sec_pid_%d_mat",i);
1905         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1906         h = fMCSecTrackHist1[2]->Projection(0);
1907         h->SetName(name);
1908         aFolderObj->Add(h);
1909
1910         sprintf(name,"mc_eta_rec_sec_pid_%d_mat",i);
1911         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1912         h = fMCSecTrackHist1[2]->Projection(1);
1913         h->SetName(name);
1914         aFolderObj->Add(h);
1915
1916         sprintf(name,"mc_eta_mother_rec_sec_pid_%d_mat",i);
1917         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1918         h = fMCSecTrackHist1[2]->Projection(4,1);
1919         h->SetName(name);
1920         aFolderObj->Add(h);
1921
1922         sprintf(name,"mc_mother_rec_sec_pid_%d_mat",i);
1923         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1924         h = fMCSecTrackHist1[2]->Projection(4);
1925         h->SetName(name);
1926         aFolderObj->Add(h);
1927
1928         sprintf(name,"mc_pt_mother_rec_sec_pid_%d_mat",i);
1929         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1930         h = fMCSecTrackHist1[2]->Projection(4,0);
1931         h->SetName(name);
1932         aFolderObj->Add(h);
1933
1934      } else {
1935        continue;
1936      }
1937    }
1938
1939   }
1940   } // end fHistogramOn
1941
1942   //
1943   //  resolution histograms
1944   //  only for reconstructed tracks
1945   //
1946
1947   TH2F *h2F=0;
1948   TCanvas * c = new TCanvas("resol","resol");
1949   c->cd();
1950
1951   //
1952   fRecMCTrackHist1->GetAxis(1)->SetRangeUser(-0.8,0.79);
1953
1954   h2F = (TH2F*)fRecMCTrackHist1->Projection(2,0);
1955   h = AlidNdPtHelper::MakeResol(h2F,1,0,kTRUE,10);
1956   h->SetXTitle("p_{tmc} (GeV/c)");
1957   h->SetYTitle("(p_{t}-p_{tmc})/p_{tmc} resolution");
1958   h->Draw();
1959   h->SetName("pt_resolution_vs_mcpt");
1960   aFolderObj->Add(h);
1961
1962   h2F = (TH2F*)fRecMCTrackHist1->Projection(2,0);
1963   h = AlidNdPtHelper::MakeResol(h2F,1,1,kTRUE,10);
1964   h->SetXTitle("p_{tmc} (GeV/c)");
1965   h->SetYTitle("(p_{t}-p_{tmc})/p_{tmc} mean");
1966   h->Draw();
1967   h->SetName("dpt_mean_vs_mcpt");
1968   aFolderObj->Add(h);
1969
1970   //
1971   h2F = (TH2F*)fRecMCTrackHist1->Projection(3,0);
1972   h = AlidNdPtHelper::MakeResol(h2F,1,0,kTRUE,10);
1973   h->SetXTitle("p_{tmc} (GeV/c)");
1974   h->SetYTitle("(#eta-#eta_{mc}) resolution");
1975   h->Draw();
1976   h->SetName("eta_resolution_vs_mcpt");
1977   aFolderObj->Add(h);
1978
1979   h2F = (TH2F*)fRecMCTrackHist1->Projection(3,0);
1980   h = AlidNdPtHelper::MakeResol(h2F,1,1,kTRUE,10);
1981   h->SetXTitle("p_{tmc} (GeV/c)");
1982   h->SetYTitle("(#eta-mc#eta) mean");
1983   h->Draw();
1984   h->SetName("deta_mean_vs_mcpt");
1985   aFolderObj->Add(h);
1986
1987   //
1988   fRecMCTrackHist1->GetAxis(1)->SetRange(1,fRecMCTrackHist1->GetAxis(1)->GetNbins());
1989
1990   h2F = (TH2F*)fRecMCTrackHist1->Projection(2,1);
1991   h = AlidNdPtHelper::MakeResol(h2F,1,0,kTRUE,10);
1992   h->SetXTitle("#eta_{mc}");
1993   h->SetYTitle("(p_{t}-p_{tmc})/p_{tmc} resolution");
1994   h->Draw();
1995   h->SetName("pt_resolution_vs_mceta");
1996   aFolderObj->Add(h);
1997
1998   h2F = (TH2F*)fRecMCTrackHist1->Projection(2,1);
1999   h = AlidNdPtHelper::MakeResol(h2F,1,1,kTRUE,10);
2000   h->SetXTitle("#eta_{mc}");
2001   h->SetYTitle("(p_{t}-p_{tmc})/p_{tmc} mean");
2002   h->Draw();
2003   h->SetName("dpt_mean_vs_mceta");
2004   aFolderObj->Add(h);
2005
2006   //
2007   h2F = (TH2F*)fRecMCTrackHist1->Projection(3,1);
2008   h = AlidNdPtHelper::MakeResol(h2F,1,0,kTRUE,10);
2009   h->SetXTitle("#eta_{mc}");
2010   h->SetYTitle("(#eta-#eta_{mc}) resolution");
2011   h->Draw();
2012   h->SetName("eta_resolution_vs_mceta");
2013   aFolderObj->Add(h);
2014
2015   h2F = (TH2F*)fRecMCTrackHist1->Projection(3,1);
2016   h = AlidNdPtHelper::MakeResol(h2F,1,1,kTRUE,10);
2017   h->SetXTitle("#eta_{mc}");
2018   h->SetYTitle("(#eta-mc#eta) mean");
2019   h->Draw();
2020   h->SetName("deta_mean_vs_mceta");
2021   aFolderObj->Add(h);
2022
2023   fRecMCTrackHist1->GetAxis(0)->SetRange(1,fRecMCTrackHist1->GetAxis(0)->GetNbins());
2024
2025   } // end use MC info
2026
2027   // export objects to analysis folder
2028   fAnalysisFolder = ExportToFolder(aFolderObj);
2029
2030   // delete only TObjArray
2031   if(aFolderObj) delete aFolderObj;
2032 }
2033
2034 //_____________________________________________________________________________
2035 TFolder* AlidNdPtAnalysisPbPb::ExportToFolder(TObjArray * const array)
2036 {
2037   // recreate folder avery time and export objects to new one
2038   //
2039   AlidNdPtAnalysisPbPb * comp=this;
2040   TFolder *folder = comp->GetAnalysisFolder();
2041
2042   TString name, title;
2043   TFolder *newFolder = 0;
2044   Int_t i = 0;
2045   Int_t size = array->GetSize();
2046
2047   if(folder) {
2048      // get name and title from old folder
2049      name = folder->GetName();
2050      title = folder->GetTitle();
2051
2052          // delete old one
2053      delete folder;
2054
2055          // create new one
2056      newFolder = CreateFolder(name.Data(),title.Data());
2057      newFolder->SetOwner();
2058
2059          // add objects to folder
2060      while(i < size) {
2061            newFolder->Add(array->At(i));
2062            i++;
2063          }
2064   }
2065
2066 return newFolder;
2067 }
2068
2069 //_____________________________________________________________________________
2070 TFolder* AlidNdPtAnalysisPbPb::CreateFolder(TString name,TString title) {
2071 // create folder for analysed histograms
2072 //
2073 TFolder *folder = 0;
2074   folder = new TFolder(name.Data(),title.Data());
2075
2076   return folder;
2077 }