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