PWGLF/SPECTRA/ChargedHadrons/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   snprintf(name,256,"fMCTrackHist1_%d",i);
 555   snprintf(title,256,"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   snprintf(name,256,"fMCPrimTrackHist1_%d",i);
 571   snprintf(title,256,"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   snprintf(name,256,"fMCPrimTrackHist2_%d",i);
 589   snprintf(title,256,"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   snprintf(name,256,"fMCSecTrackHist1_%d",i);
 603   snprintf(title,256,"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   snprintf(name,256,"fRecTrackHist1_%d",i);
 621   snprintf(title,256,"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   snprintf(name,256,"fRecTrackHist2_%d",i);
 633   snprintf(title,256,"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   snprintf(name,256,"fRecTrackMultHist_%d",i);
 650   snprintf(title,256,"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   snprintf(name,256,"fRecMCTrackHist1");
 666   snprintf(title,256,"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   snprintf(name,256,"fMCMultRecTrackHist1");
 681   snprintf(title,256,"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     isEventTriggered = inputHandler->IsEventSelected() & GetTriggerMask();
 762
 763     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
 764     if(!physicsSelection) return;
 765     //SetPhysicsTriggerSelection(physicsSelection);
 766
 767     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
 768       // set trigger (V0AND)
 769       triggerAnalysis = physicsSelection->GetTriggerAnalysis();
 770       if(!triggerAnalysis) return;
 771       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
 772     }
 773   }
 774
 775
 776   // centrality determination
 777   Float_t centralityF = -1.;
 778   AliCentrality *esdCentrality = esdEvent->GetCentrality();
 779   centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
 780
 781   // use MC information
 782   AliHeader* header = 0;
 783   AliGenEventHeader* genHeader = 0;
 784   AliStack* stack = 0;
 785   TArrayF vtxMC(3);
 786
 787   Int_t multMCTrueTracks = 0;
 788   if(IsUseMCInfo())
 789   {
 790     //
 791     if(!mcEvent) {
 792       AliDebug(AliLog::kError, "mcEvent not available");
 793       return;
 794     }
 795     // get MC event header
 796     header = mcEvent->Header();
 797     if (!header) {
 798       AliDebug(AliLog::kError, "Header not available");
 799       return;
 800     }
 801     // MC particle stack
 802     stack = mcEvent->Stack();
 803     if (!stack) {
 804       AliDebug(AliLog::kError, "Stack not available");
 805       return;
 806     }
 807
 808     // get MC vertex
 809     genHeader = header->GenEventHeader();
 810     if (!genHeader) {
 811       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
 812       return;
 813     }
 814     genHeader->PrimaryVertex(vtxMC);
 815
 816     Double_t vMCEventHist1[4]={vtxMC[0],vtxMC[1],vtxMC[2],centralityF};
 817     fMCEventHist1->Fill(vMCEventHist1);
 818
 819     // multipliticy of all MC primary tracks
 820     // in Zv, pt and eta ranges)
 821     multMCTrueTracks = AlidNdPtHelper::GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
 822
 823   } // end bUseMC
 824
 825   // get reconstructed vertex
 826   const AliESDVertex* vtxESD = 0;
 827   if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {
 828      vtxESD = esdEvent->GetPrimaryVertexTPC();
 829   }
 830   else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {
 831      vtxESD = esdEvent->GetPrimaryVertexTracks();
 832   }
 833   else {
 834      return;
 835   }
 836
 837   if(!vtxESD) return;
 838
 839   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
 840   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
 841   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
 842
 843   // vertex contributors
 844   Int_t multMBTracks = 0;
 845   if(GetAnalysisMode() == AlidNdPtHelper::kTPC || GetAnalysisMode() == AlidNdPtHelper::kTPCITS)
 846   {
 847      if(vtxESD->GetStatus()) {
 848          multMBTracks = vtxESD->GetNContributors();
 849      }
 850   }
 851   else {
 852     AliDebug(AliLog::kError, Form("Found analysis type %d", GetAnalysisMode()));
 853     return;
 854   }
 855
 856   TObjArray *allChargedTracks=0;
 857   //Int_t multAll=0, multAcc=0, multRec=0;
 858   Int_t multAll=0, multRec=0;
 859   Int_t *labelsAll=0, *labelsAcc=0, *labelsRec=0;
 860
 861   // check event cuts
 862   if(isEventOK && isEventTriggered)
 863   {
 864
 865     // get all charged tracks
 866     allChargedTracks = AlidNdPtHelper::GetAllChargedTracks(esdEvent,GetAnalysisMode());
 867     if(!allChargedTracks) return;
 868
 869     Int_t entries = allChargedTracks->GetEntries();
 870
 871     labelsAll = new Int_t[entries];
 872     labelsAcc = new Int_t[entries];
 873     labelsRec = new Int_t[entries];
 874     for(Int_t i=0; i<entries;++i)
 875     {
 876       AliESDtrack *track = (AliESDtrack*)allChargedTracks->At(i);
 877
 878       if(!track) continue;
 879       if(track->Charge()==0) continue;
 880
 881
 882       // only postive charged
 883       if(GetParticleMode() == AlidNdPtHelper::kPlus && track->Charge() < 0)
 884         continue;
 885
 886       // only negative charged
 887       if(GetParticleMode() == AlidNdPtHelper::kMinus && track->Charge() > 0)
 888         continue;
 889
 890       //
 891       Double_t values[4] = {vtxESD->GetZv(),track->Pt(),track->Eta(), centralityF};
 892
 893       fRecTrackHist2[AlidNdPtHelper::kAllTracks]->Fill(values);
 894       FillHistograms(track,stack,AlidNdPtHelper::kAllTracks,centralityF);
 895       labelsAll[multAll] = TMath::Abs(track->GetLabel());
 896
 897       multAll++;
 898       if(esdTrackCuts->AcceptTrack(track) && accCuts->AcceptTrack(track) && recCuts->AcceptTrackLocalTPC(track)) {
 899
 900          fRecTrackHist2[AlidNdPtHelper::kRecTracks]->Fill(values);
 901          FillHistograms(track,stack,AlidNdPtHelper::kRecTracks,centralityF);
 902          labelsRec[multRec] = TMath::Abs(track->GetLabel());
 903
 904          multRec++;
 905
 906       }
 907      }//loop over entries
 908
 909
 910      // fill track multiplicity histograms
 911      //FillHistograms(allChargedTracks,labelsAll,multAll,labelsAcc,multAcc,labelsRec,multRec);
 912
 913      Double_t vRecEventHist1[4] = {vtxESD->GetXv(),vtxESD->GetYv(),vtxESD->GetZv(),centralityF};
 914      fRecEventHist1->Fill(vRecEventHist1);
 915
 916      Double_t vRecEventHist2[3] = {vtxESD->GetZv(),multMBTracks,centralityF};
 917      fRecEventHist2->Fill(vRecEventHist2);
 918
 919    } // triggered and event vertex
 920
 921    if(IsUseMCInfo())
 922    {
 923      //
 924      // event level corrections (zv,N_MB)
 925      //
 926      // all inelastic
 927      Double_t vEventMatrix[3] = {vtxMC[2],multMBTracks,centralityF};
 928      fGenEventMatrix->Fill(vEventMatrix);
 929      if(isEventTriggered) fTriggerEventMatrix->Fill(vEventMatrix);
 930      if(isEventOK && isEventTriggered) fRecEventMatrix->Fill(vEventMatrix);
 931
 932      //
 933      // track-event level corrections (zv,pt,eta)
 934      //
 935      for (Int_t iMc = 0; iMc < stack->GetNtrack(); ++iMc)
 936      {
 937        TParticle* particle = stack->Particle(iMc);
 938        if (!particle)
 939        continue;
 940
 941        // only charged particles
 942        if(!particle->GetPDG()) continue;
 943        Double_t charge = particle->GetPDG()->Charge()/3.;
 944        if ( TMath::Abs(charge) < 0.001 )
 945         continue;
 946
 947        // only postive charged
 948        if(GetParticleMode() == AlidNdPtHelper::kPlus && charge < 0.)
 949         continue;
 950
 951        // only negative charged
 952        if(GetParticleMode() == AlidNdPtHelper::kMinus && charge > 0.)
 953        continue;
 954
 955        // physical primary
 956        Bool_t prim = stack->IsPhysicalPrimary(iMc);
 957        if(!prim) continue;
 958
 959        // checked accepted
 960        if(accCuts->AcceptTrack(particle))
 961        {
 962         Double_t vTrackEventMatrix[4] = {vtxMC[2], particle->Pt(), particle->Eta(), centralityF};
 963         fGenTrackEventMatrix->Fill(vTrackEventMatrix);
 964
 965         if(!isEventTriggered) continue;
 966         fTriggerTrackEventMatrix->Fill(vTrackEventMatrix);
 967
 968         if(!isEventOK) continue;
 969         fRecTrackEventMatrix->Fill(vTrackEventMatrix);
 970
 971        }// if(accCuts->AcceptTrack(particle))
 972      }// for (Int_t iMc = 0; iMc < stack->GetNtrack(); ++iMc)
 973
 974      //
 975      // track-level corrections (zv,pt,eta)
 976      //
 977      if(isEventOK && isEventTriggered)
 978      {
 979
 980        // fill MC and rec event control histograms
 981        if(fHistogramsOn) {
 982          Double_t vRecMCEventHist1[4] = {vtxESD->GetXv()-vtxMC[0],vtxESD->GetYv()-vtxMC[1],vtxESD->GetZv()-vtxMC[2], centralityF};
 983          fRecMCEventHist1->Fill(vRecMCEventHist1);//
 984
 985          Double_t vRecMCEventHist2[4] = {vtxESD->GetXv()-vtxMC[0],vtxESD->GetZv()-vtxMC[2],multMBTracks, centralityF};
 986          fRecMCEventHist2->Fill(vRecMCEventHist2);
 987
 988        }//
 989
 990        //
 991        // MC histograms for track efficiency studies
 992        //
 993        for (Int_t iMc = 0; iMc < stack->GetNtrack(); ++iMc)
 994        {
 995          TParticle* particle = stack->Particle(iMc);
 996          if (!particle)
 997          continue;
 998
 999          Double_t vTrackMatrix[4] = {vtxMC[2],particle->Pt(),particle->Eta(),centralityF};
1000
1001          // only charged particles
1002          if(!particle->GetPDG()) continue;
1003          Double_t charge = particle->GetPDG()->Charge()/3.;
1004          if (TMath::Abs(charge) < 0.001)
1005          continue;
1006
1007          // only postive charged
1008          if(GetParticleMode() == AlidNdPtHelper::kPlus && charge < 0.)
1009          continue;
1010
1011          // only negative charged
1012          if(GetParticleMode() == AlidNdPtHelper::kMinus && charge > 0.)
1013          continue;
1014
1015          // physical primary
1016          Bool_t prim = stack->IsPhysicalPrimary(iMc);
1017
1018          // check accepted
1019          if(accCuts->AcceptTrack(particle))
1020          {
1021
1022            if( AlidNdPtHelper::IsPrimaryParticle(stack, iMc, GetParticleMode()) )
1023              fGenPrimTrackMatrix->Fill(vTrackMatrix);
1024
1025            // fill control histograms
1026            if(fHistogramsOn)
1027              FillHistograms(stack,iMc,AlidNdPtHelper::kAccTracks, centralityF);
1028
1029            // check multiple found tracks
1030            Int_t multCount = 0;
1031            for(Int_t iRec=0; iRec<multRec; ++iRec)
1032            {
1033              if(iMc == labelsRec[iRec])
1034              {
1035                multCount++;
1036                if(multCount>1)
1037                {
1038                  fRecMultTrackMatrix->Fill(vTrackMatrix);
1039
1040                  // fill control histogram
1041                  if(fHistogramsOn) {
1042                    Int_t pid = AlidNdPtHelper::ConvertPdgToPid(particle);
1043                    Double_t vMCMultRecTrackHist1[4] = {particle->Pt(), particle->Eta(), pid, centralityF};
1044                    fMCMultRecTrackHist1->Fill(vMCMultRecTrackHist1);
1045                  }
1046                }
1047              }
1048            }
1049
1050            // check reconstructed
1051            for(Int_t iRec=0; iRec<multRec; ++iRec)
1052            {
1053              if(iMc == labelsRec[iRec])
1054              {
1055                fRecTrackMatrix->Fill(vTrackMatrix);
1056
1057                if( AlidNdPtHelper::IsPrimaryParticle(stack, iMc, GetParticleMode()) ) {
1058                  fRecPrimTrackMatrix->Fill(vTrackMatrix);
1059                }
1060                if(!prim) fRecSecTrackMatrix->Fill(vTrackMatrix);
1061
1062                // fill control histograms
1063                if(fHistogramsOn)
1064                  FillHistograms(stack,iMc,AlidNdPtHelper::kRecTracks, centralityF);
1065
1066                break;
1067              }//if(iMc == labelsRec[iRec])
1068            }//reco tracks
1069          }//accepted tracks
1070        }//stack loop
1071      }//is triggered
1072    } // end bUseMC
1073
1074   if(allChargedTracks) delete allChargedTracks; allChargedTracks = 0;
1075   if(labelsAll) delete [] labelsAll; labelsAll = 0;
1076   if(labelsAcc) delete [] labelsAcc; labelsAcc = 0;
1077   if(labelsRec) delete [] labelsRec; labelsRec = 0;
1078
1079 }
1080
1081 //_____________________________________________________________________________
1082 void AlidNdPtAnalysisPbPb::FillHistograms(TObjArray *const allChargedTracks,Int_t *const labelsAll,Int_t multAll,Int_t *const labelsAcc,Int_t multAcc,Int_t *const labelsRec,Int_t multRec, Float_t centralityF) {
1083  // multiplicity  histograms
1084
1085
1086   if(!allChargedTracks) return;
1087   if(!labelsAll) return;
1088   if(!labelsAcc) return;
1089   if(!labelsRec) return;
1090
1091   Int_t entries = allChargedTracks->GetEntries();
1092   for(Int_t i=0; i<entries; ++i)
1093   {
1094      AliESDtrack *track = (AliESDtrack*)allChargedTracks->At(i);
1095      if(!track) continue;
1096      if(track->Charge() == 0) continue;
1097
1098      Int_t label = TMath::Abs(track->GetLabel());
1099      for(Int_t iAll=0; iAll<multAll; ++iAll) {
1100        if(label == labelsAll[iAll]) {
1101          Double_t v1[3] = {track->Pt(), multAll, centralityF};
1102          fRecTrackMultHist1[AlidNdPtHelper::kAllTracks]->Fill(v1);
1103        }
1104      }
1105      for(Int_t iAcc=0; iAcc<multAcc; ++iAcc) {
1106        if(label == labelsAcc[iAcc]) {
1107          Double_t v2[3] = {track->Pt(), multAcc, centralityF};
1108          fRecTrackMultHist1[AlidNdPtHelper::kAccTracks]->Fill(v2);
1109        }
1110      }
1111      for(Int_t iRec=0; iRec<multRec; ++iRec) {
1112        if(label == labelsRec[iRec]) {
1113          Double_t v3[3] = {track->Pt(), multRec, centralityF};
1114          fRecTrackMultHist1[AlidNdPtHelper::kRecTracks]->Fill(v3);
1115        }//out
1116      }
1117   }
1118 }
1119
1120 //_____________________________________________________________________________
1121 void AlidNdPtAnalysisPbPb::FillHistograms(AliESDtrack *const esdTrack, AliStack *const stack, AlidNdPtHelper::TrackObject trackObj, Float_t centralityF)
1122 {
1123
1124   //
1125   // Fill ESD track and MC histograms
1126   //
1127   if(!esdTrack) return;
1128
1129   Float_t q = esdTrack->Charge();
1130   if(TMath::Abs(q) < 0.001) return;
1131
1132   Float_t pt = esdTrack->Pt();
1133   Float_t eta = esdTrack->Eta();
1134   Float_t phi = esdTrack->Phi();
1135
1136   Float_t dca[2], bCov[3];
1137   esdTrack->GetImpactParameters(dca,bCov);
1138
1139   Int_t nClust = esdTrack->GetTPCclusters(0);
1140   Float_t chi2PerCluster = 0.;
1141   if(nClust>0.) chi2PerCluster = esdTrack->GetTPCchi2()/Float_t(nClust);
1142
1143
1144   // fill histograms
1145   Double_t values[4] = {pt,eta,phi,centralityF};
1146   fRecTrackHist1[trackObj]->Fill(values);
1147
1148   //
1149   // Fill rec vs MC information
1150   //
1151   if(!stack) return;
1152
1153   Int_t label = TMath::Abs(esdTrack->GetLabel());
1154   //if(label == 0) return;
1155
1156   if(label > stack->GetNtrack()) return;
1157   TParticle* particle = stack->Particle(label);
1158   if(!particle) return;
1159
1160   Int_t motherPdg = -1;
1161   TParticle* mother = 0;
1162
1163   //TParticle* prim_mother = AlidNdPtHelper::FindPrimaryMother(stack,label);
1164   Int_t motherLabel = particle->GetMother(0);
1165   if(motherLabel>0) mother = stack->Particle(motherLabel);
1166   if(mother) motherPdg = TMath::Abs(mother->GetPdgCode()); // take abs for visualisation only
1167   //Int_t mech = particle->GetUniqueID(); // production mechanism
1168
1169   if(!particle->GetPDG()) return;
1170   Double_t gq = particle->GetPDG()->Charge()/3.0; // Charge units |e|/3
1171   if(TMath::Abs(gq)<0.001) return;
1172   Float_t gpt = particle->Pt();
1173   Float_t geta = particle->Eta();
1174
1175   Double_t dpt=0;
1176   //printf("pt %f, gpt %f \n",pt,gpt);
1177   if(gpt) dpt = (pt-gpt)/gpt;
1178   Double_t deta = (eta-geta);
1179
1180   // fill histograms
1181   if(trackObj == AlidNdPtHelper::kRecTracks)  //RecTracks???
1182   {
1183     Double_t vTrackPtCorrelationMatrix[4]={pt,gpt,geta,centralityF};
1184     fTrackPtCorrelationMatrix->Fill(vTrackPtCorrelationMatrix);
1185
1186     Double_t vRecMCTrackHist1[5]={gpt,geta,dpt,deta,centralityF};
1187     fRecMCTrackHist1->Fill(vRecMCTrackHist1);
1188   }
1189 }
1190
1191 //_____________________________________________________________________________
1192 void AlidNdPtAnalysisPbPb::FillHistograms(AliStack *const stack, Int_t label, AlidNdPtHelper::TrackObject trackObj, Float_t centralityF)
1193 {
1194
1195   // Fill MC histograms
1196   if(!stack) return;
1197
1198   if(label > stack->GetNtrack()) return;
1199   TParticle* particle = stack->Particle(label);
1200   if(!particle) return;
1201
1202   Int_t motherPdg = -1;
1203   TParticle* mother = 0;
1204
1205   //TParticle* prim_mother = AlidNdPtHelper::FindPrimaryMother(stack,label);
1206   Int_t motherLabel = particle->GetMother(0);
1207   if(motherLabel>0) mother = stack->Particle(motherLabel);
1208   if(mother) motherPdg = TMath::Abs(mother->GetPdgCode()); // take abs for visualisation only
1209   Int_t mech = particle->GetUniqueID(); // production mechanism
1210
1211   if(!particle->GetPDG()) return;
1212   Double_t gq = particle->GetPDG()->Charge()/3.0; // Charge units |e|/3
1213   if(TMath::Abs(gq) < 0.001) return;
1214
1215   Float_t gpt = particle->Pt();
1216   //Float_t qgpt = particle->Pt() * gq;
1217   Float_t geta = particle->Eta();
1218   Float_t gphi = particle->Phi();
1219   //Float_t gpz = particle->Pz();
1220
1221   Bool_t prim = stack->IsPhysicalPrimary(label);
1222   //Float_t vx = particle->Vx(); Float_t vy = particle->Vy(); Float_t vz = particle->Vz();
1223
1224   Int_t pid = AlidNdPtHelper::ConvertPdgToPid(particle);
1225
1226   //if(prim&&pid==5) printf("pdgcode %d, production mech %d \n",particle->GetPdgCode(),mech);
1227   //if(!prim) printf("motherPdg %d, particle %d, production mech %d\n",motherPdg, particle->GetPdgCode(),mech);
1228
1229   //
1230   // fill histogram
1231   //
1232   Double_t vMCTrackHist1[4] = {gpt,geta,gphi,centralityF};
1233   fMCTrackHist1[trackObj]->Fill(vMCTrackHist1);
1234
1235   Double_t vMCPrimTrackHist1[6] = {gpt,geta,pid,mech,motherPdg,centralityF};
1236   Double_t vMCPrimTrackHist2[4] = {TMath::Abs(particle->GetPdgCode()),mech,motherPdg,centralityF};
1237
1238   //if(prim && AliPWG0Helper::IsPrimaryCharged(particle, label)) fMCPrimTrackHist1[trackObj]->Fill(vMCPrimTrackHist1);
1239
1240   if(prim) {
1241     fMCPrimTrackHist1[trackObj]->Fill(vMCPrimTrackHist1);
1242     if(pid == 5) fMCPrimTrackHist2[trackObj]->Fill(vMCPrimTrackHist2);
1243   }
1244   else {
1245     fMCSecTrackHist1[trackObj]->Fill(vMCPrimTrackHist1);
1246   }
1247
1248 }
1249
1250 //_____________________________________________________________________________
1251 Long64_t AlidNdPtAnalysisPbPb::Merge(TCollection* const list)
1252 {
1253   // Merge list of objects (needed by PROOF)
1254
1255   //  init if not done already
1256   if (!fIsInit) { Init(); }
1257
1258   if (!list)
1259   return 0;
1260
1261   if (list->IsEmpty())
1262   return 1;
1263
1264   TIterator* iter = list->MakeIterator();
1265   TObject* obj = 0;
1266
1267   //
1268   //TList *collPhysSelection = new TList;
1269
1270   // collection of generated histograms
1271
1272   Int_t count=0;
1273   while((obj = iter->Next()) != 0) {
1274     AlidNdPtAnalysisPbPb* entry = dynamic_cast<AlidNdPtAnalysisPbPb*>(obj);
1275     if (entry == 0) continue;
1276
1277     // physics selection
1278     //printf("entry->GetPhysicsTriggerSelection() %p \n", entry->GetPhysicsTriggerSelection());
1279     //AliPhysicsSelection *physSel = entry->GetPhysicsTriggerSelection();
1280     //if( physSel ) collPhysSelection->Add(physSel);
1281
1282     //
1283     fTrackPtCorrelationMatrix->Add(entry->fTrackPtCorrelationMatrix);
1284
1285     //
1286     fGenEventMatrix->Add(entry->fGenEventMatrix);
1287
1288     fTriggerEventMatrix->Add(entry->fTriggerEventMatrix);
1289
1290     fRecEventMatrix->Add(entry->fRecEventMatrix);
1291     //
1292     fGenTrackEventMatrix->Add(entry->fGenTrackEventMatrix);
1293
1294     fTriggerTrackEventMatrix->Add(entry->fTriggerTrackEventMatrix);
1295
1296     fRecTrackEventMatrix->Add(entry->fRecTrackEventMatrix);
1297
1298     //
1299     fGenTrackMatrix->Add(entry->fGenTrackMatrix);
1300     fGenPrimTrackMatrix->Add(entry->fGenPrimTrackMatrix);
1301     fRecPrimTrackMatrix->Add(entry->fRecPrimTrackMatrix);
1302     //
1303     fRecTrackMatrix->Add(entry->fRecTrackMatrix);
1304     fRecSecTrackMatrix->Add(entry->fRecSecTrackMatrix);
1305     //
1306     fRecMultTrackMatrix->Add(entry->fRecMultTrackMatrix);
1307
1308     //
1309     // control analysis histograms
1310     //
1311     fMCEventHist1->Add(entry->fMCEventHist1);
1312     fRecEventHist1->Add(entry->fRecEventHist1);
1313     fRecEventHist2->Add(entry->fRecEventHist2);
1314     fRecMCEventHist1->Add(entry->fRecMCEventHist1);
1315     fRecMCEventHist2->Add(entry->fRecMCEventHist2);
1316
1317
1318     for(Int_t i=0; i<AlidNdPtHelper::kCutSteps; i++) {
1319       fMCTrackHist1[i]->Add(entry->fMCTrackHist1[i]);
1320
1321       fMCPrimTrackHist1[i]->Add(entry->fMCPrimTrackHist1[i]);
1322       fMCPrimTrackHist2[i]->Add(entry->fMCPrimTrackHist2[i]);
1323       fMCSecTrackHist1[i]->Add(entry->fMCSecTrackHist1[i]);
1324
1325       fRecTrackHist1[i]->Add(entry->fRecTrackHist1[i]);
1326       fRecTrackHist2[i]->Add(entry->fRecTrackHist2[i]);
1327       fRecTrackMultHist1[i]->Add(entry->fRecTrackMultHist1[i]);
1328     }
1329     fRecMCTrackHist1->Add(entry->fRecMCTrackHist1);
1330     fMCMultRecTrackHist1->Add(entry->fMCMultRecTrackHist1);
1331     fRecTrackHist3->Add(entry->fRecTrackHist3);
1332
1333   count++;
1334   }
1335
1336   //AliPhysicsSelection *trigSelection = GetPhysicsTriggerSelection();
1337   //if( trigSelection ) trigSelection->Merge(collPhysSelection);
1338   //if(collPhysSelection) delete collPhysSelection;
1339
1340 return count;
1341 }
1342
1343
1344
1345 //_____________________________________________________________________________
1346 void AlidNdPtAnalysisPbPb::Analyse()
1347 {
1348
1349   //  init if not done already
1350   if (!fIsInit) { Init(); }
1351
1352   // Analyse histograms
1353   //
1354   TH1::AddDirectory(kFALSE);
1355   TH1 *h=0, *h1=0, *h2=0, *h2c = 0;
1356   THnSparse *hs=0;
1357   TH2 *h2D=0;
1358
1359   char name[256];
1360   TObjArray *aFolderObj = new TObjArray;
1361   if(!aFolderObj) return;
1362
1363   //
1364   // LHC backgraund in all and 0-bins
1365   //
1366   //AliPhysicsSelection *trigSel = GetPhysicsTriggerSelection();
1367   //trigSel->SaveHistograms("physics_selection");
1368
1369   //
1370   // Reconstructed event vertex
1371   //
1372   h = fRecEventHist1->Projection(0);
1373   h->SetName("Xv");
1374   aFolderObj->Add(h);
1375
1376   h = fRecEventHist1->Projection(1);
1377   h->SetName("Yv");
1378   aFolderObj->Add(h);
1379
1380   h = fRecEventHist1->Projection(2);
1381   h->SetName("Zv");
1382   aFolderObj->Add(h);
1383
1384   //
1385   // multiplicity
1386   //
1387   h = fRecEventHist2->Projection(1);
1388   h->SetName("multMB");
1389   aFolderObj->Add(h);
1390
1391   h2D = fRecEventHist2->Projection(0,1);
1392   h2D->SetName("Zv_vs_multiplicity_MB");
1393   aFolderObj->Add(h2D);
1394
1395   //
1396   // reconstructed pt histograms
1397   //
1398   h = fRecTrackHist1[0]->Projection(0);
1399   h->Scale(1.,"width");
1400   h->SetName("pt_all_ch");
1401   aFolderObj->Add(h);
1402
1403   h = fRecTrackHist1[1]->Projection(0);
1404   h->Scale(1.,"width");
1405   h->SetName("pt_acc");
1406   aFolderObj->Add(h);
1407
1408   h = fRecTrackHist1[2]->Projection(0);
1409   h->Scale(1.,"width");
1410   h->SetName("pt_rec");
1411   aFolderObj->Add(h);
1412
1413   //
1414   // reconstructed eta histograms
1415   //
1416   h = fRecTrackHist1[0]->Projection(1);
1417   h->SetName("eta_all_ch");
1418   aFolderObj->Add(h);
1419
1420   h = fRecTrackHist1[1]->Projection(1);
1421   h->SetName("eta_acc");
1422   aFolderObj->Add(h);
1423
1424   h = fRecTrackHist1[2]->Projection(1);
1425   h->SetName("eta_rec");
1426   aFolderObj->Add(h);
1427
1428   //
1429   // reconstructed phi histograms
1430   //
1431   h = fRecTrackHist1[0]->Projection(2);
1432   h->SetName("phi_all_ch");
1433   aFolderObj->Add(h);
1434
1435   h = fRecTrackHist1[1]->Projection(2);
1436   h->SetName("phi_acc");
1437   aFolderObj->Add(h);
1438
1439   h = fRecTrackHist1[2]->Projection(2);
1440   h->SetName("phi_rec");
1441   aFolderObj->Add(h);
1442
1443   //
1444   // reconstructed eta:pt histograms
1445   //
1446   h2D = fRecTrackHist1[0]->Projection(1,0);
1447   h2D->SetName("pt_eta_all_ch");
1448   aFolderObj->Add(h2D);
1449
1450   h2D = fRecTrackHist1[1]->Projection(1,0);
1451   h2D->SetName("pt_eta_acc");
1452   aFolderObj->Add(h2D);
1453
1454   h2D = fRecTrackHist1[2]->Projection(1,0);
1455   h2D->SetName("pt_eta_rec");
1456   aFolderObj->Add(h2D);
1457
1458   //
1459   // reconstructed phi:pt histograms
1460   //
1461   h2D = fRecTrackHist1[0]->Projection(2,0);
1462   h2D->SetName("pt_phi_all_ch");
1463   aFolderObj->Add(h2D);
1464
1465   h2D = fRecTrackHist1[1]->Projection(2,0);
1466   h2D->SetName("pt_phi_acc");
1467   aFolderObj->Add(h2D);
1468
1469   h2D = fRecTrackHist1[2]->Projection(2,0);
1470   h2D->SetName("pt_phi_rec");
1471   aFolderObj->Add(h2D);
1472
1473   //
1474   // reconstructed phi:eta histograms
1475   //
1476   h2D = fRecTrackHist1[0]->Projection(2,1);
1477   h2D->SetName("eta_phi_all_ch");
1478   aFolderObj->Add(h2D);
1479
1480   h2D = fRecTrackHist1[1]->Projection(2,1);
1481   h2D->SetName("eta_phi_acc");
1482   aFolderObj->Add(h2D);
1483
1484   h2D = fRecTrackHist1[2]->Projection(2,1);
1485   h2D->SetName("eta_phi_rec");
1486   aFolderObj->Add(h2D);
1487
1488   //
1489   // reconstructed nClust, chi2 vs pt, eta, phi
1490   //
1491   if(fHistogramsOn) {
1492
1493     h2D = fRecTrackHist3->Projection(0,1);
1494     h2D->SetName("nClust_chi2_rec");
1495     aFolderObj->Add(h2D);
1496
1497     h2D = fRecTrackHist3->Projection(0,2);
1498     h2D->SetName("nClust_pt_rec");
1499     aFolderObj->Add(h2D);
1500
1501     h2D = fRecTrackHist3->Projection(0,3);
1502     h2D->SetName("nClust_eta_rec");
1503     aFolderObj->Add(h2D);
1504
1505     h2D = fRecTrackHist3->Projection(0,4);
1506     h2D->SetName("nClust_phi_rec");
1507     aFolderObj->Add(h2D);
1508
1509     h2D = fRecTrackHist3->Projection(1,2);
1510     h2D->SetName("chi2_pt_rec");
1511     aFolderObj->Add(h2D);
1512
1513     h2D = fRecTrackHist3->Projection(1,3);
1514     h2D->SetName("chi2_eta_rec");
1515     aFolderObj->Add(h2D);
1516
1517     h2D = fRecTrackHist3->Projection(1,4);
1518     h2D->SetName("chi2_phi_rec");
1519     aFolderObj->Add(h2D);
1520
1521   }
1522
1523   //
1524   // calculate corrections for empty events
1525   // with multMB==0
1526   //
1527
1528   //
1529   // normalised zv to generate zv for triggered events
1530   //
1531   h = fRecEventHist2->Projection(0);
1532   if( h->Integral() ) h->Scale(1./h->Integral());
1533   h->SetName("zv_distribution_norm");
1534   aFolderObj->Add(h);
1535
1536   //
1537   // MC available
1538   //
1539   if(IsUseMCInfo()) {
1540
1541   //
1542   // Event vertex resolution
1543   //
1544   h2D = fRecMCEventHist2->Projection(0,2);
1545   h2D->SetName("DeltaXv_vs_mult");
1546   aFolderObj->Add(h2D);
1547
1548   h2D = fRecMCEventHist2->Projection(1,2);
1549   h2D->SetName("DeltaZv_vs_mult");
1550   aFolderObj->Add(h2D);
1551
1552   //
1553   // normalised zv to get trigger/trigger+vertex event differences
1554   // F(zv) = E_trig(zv,0)/Int(E_trig(zv,0) / Sum(E_trigvtx(zv,n))/Sum(Int(E_trigvtx(zv,n))dzv)
1555   //
1556   fTriggerEventMatrix->GetAxis(1)->SetRangeUser(0.,0.);
1557   h = fTriggerEventMatrix->Projection(0);
1558   h2D = fTriggerEventMatrix->Projection(0,1);
1559   if(h2D->Integral()) h->Scale(1./h2D->Integral());
1560
1561   h1 = fRecEventMatrix->Projection(0);
1562   h2D = fRecEventMatrix->Projection(0,1);
1563   if(h2D->Integral()) h1->Scale(1./h2D->Integral());
1564
1565   h->Divide(h1);
1566   h->SetName("zv_empty_events_norm");
1567   aFolderObj->Add(h);
1568
1569   fTriggerEventMatrix->GetAxis(1)->SetRange(1,fTriggerEventMatrix->GetAxis(1)->GetNbins());
1570
1571   //
1572   // rec. vs true track pt correlation matrix
1573   //
1574   hs = (THnSparse*)fTrackPtCorrelationMatrix->Clone("track_pt_correlation_matrix");
1575   aFolderObj->Add(hs);
1576
1577   //
1578   // trigger efficiency for INEL
1579   //
1580   h = AlidNdPtHelper::GenerateCorrMatrix(fTriggerEventMatrix->Projection(0),fGenEventMatrix->Projection(0),"zv_trig_INEL_eff_matrix");
1581   aFolderObj->Add(h);
1582
1583
1584   //
1585   // trigger bias correction (MB to INEL)
1586   //
1587   hs = AlidNdPtHelper::GenerateCorrMatrix(fGenEventMatrix,fTriggerEventMatrix,"zv_mult_trig_MBtoInel_corr_matrix");
1588   aFolderObj->Add(hs);
1589
1590   h = AlidNdPtHelper::GenerateCorrMatrix(fGenEventMatrix->Projection(0),fTriggerEventMatrix->Projection(0),"zv_trig_MBtoInel_corr_matrix");
1591   aFolderObj->Add(h);
1592
1593   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenEventMatrix->Projection(0,1),fTriggerEventMatrix->Projection(0,1),"zv_mult_trig_MBtoInel_corr_matrix_2D");
1594   aFolderObj->Add(h2D);
1595
1596
1597   h = AlidNdPtHelper::GenerateCorrMatrix(fGenEventMatrix->Projection(1),fTriggerEventMatrix->Projection(1),"mult_trig_MBtoInel_corr_matrix");
1598   aFolderObj->Add(h);
1599
1600
1601   //
1602   // event vertex reconstruction correction (MB)
1603   //
1604   hs = AlidNdPtHelper::GenerateCorrMatrix(fTriggerEventMatrix,fRecEventMatrix,"zv_mult_event_corr_matrix");
1605   aFolderObj->Add(hs);
1606
1607   h2D = AlidNdPtHelper::GenerateCorrMatrix(fTriggerEventMatrix->Projection(0,1),fRecEventMatrix->Projection(0,1),"zv_mult_event_corr_matrix_2D");
1608   aFolderObj->Add(h2D);
1609
1610
1611   h = AlidNdPtHelper::GenerateCorrMatrix(fTriggerEventMatrix->Projection(1),fRecEventMatrix->Projection(1),"mult_event_corr_matrix");
1612   aFolderObj->Add(h);
1613
1614
1615   h = AlidNdPtHelper::GenerateCorrMatrix(fTriggerEventMatrix->Projection(0),fRecEventMatrix->Projection(0),"zv_event_corr_matrix");
1616   aFolderObj->Add(h);
1617
1618
1619
1620   //
1621   // track-event trigger bias correction (MB to INEL)
1622   //
1623   hs = AlidNdPtHelper::GenerateCorrMatrix(fGenTrackEventMatrix,fTriggerTrackEventMatrix,"zv_pt_eta_track_trig_MBtoInel_corr_matrix");
1624   aFolderObj->Add(hs);
1625
1626   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenTrackEventMatrix->Projection(1,2),fTriggerTrackEventMatrix->Projection(1,2),"eta_pt_track_trig_MBtoInel_corr_matrix");
1627   aFolderObj->Add(h2D);
1628
1629   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenTrackEventMatrix->Projection(1,0),fTriggerTrackEventMatrix->Projection(1,0),"pt_zv_track_trig_MBtoInel_corr_matrix");
1630   aFolderObj->Add(h2D);
1631
1632   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenTrackEventMatrix->Projection(2,0),fTriggerTrackEventMatrix->Projection(2,0),"zv_eta_track_trig_MBtoInel_corr_matrix");
1633   aFolderObj->Add(h2D);
1634
1635   // efficiency
1636
1637   h = AlidNdPtHelper::GenerateCorrMatrix(fTriggerTrackEventMatrix->Projection(1),fGenTrackEventMatrix->Projection(1),"pt_track_trig_MBtoInel_eff_matrix");
1638   aFolderObj->Add(h);
1639
1640
1641   //
1642   // track-event vertex reconstruction correction (MB)
1643   //
1644   hs = AlidNdPtHelper::GenerateCorrMatrix(fTriggerTrackEventMatrix,fRecTrackEventMatrix,"zv_pt_eta_track_event_corr_matrix");
1645   aFolderObj->Add(hs);
1646
1647   h2D = AlidNdPtHelper::GenerateCorrMatrix(fTriggerTrackEventMatrix->Projection(1,2),fRecTrackEventMatrix->Projection(1,2),"eta_pt_track_event_corr_matrix");
1648   aFolderObj->Add(h2D);
1649
1650   h2D = AlidNdPtHelper::GenerateCorrMatrix(fTriggerTrackEventMatrix->Projection(1,0),fRecTrackEventMatrix->Projection(1,0),"pt_zv_track_event_corr_matrix");
1651   aFolderObj->Add(h2D);
1652
1653   h2D = AlidNdPtHelper::GenerateCorrMatrix(fTriggerTrackEventMatrix->Projection(2,0),fRecTrackEventMatrix->Projection(2,0),"zv_eta_track_event_corr_matrix");
1654   aFolderObj->Add(h2D);
1655
1656   // efficiency
1657
1658   h = AlidNdPtHelper::GenerateCorrMatrix(fRecTrackEventMatrix->Projection(1),fTriggerTrackEventMatrix->Projection(1),"pt_track_event_eff_matrix");
1659   aFolderObj->Add(h);
1660
1661
1662   //
1663   // track rec. efficiency correction
1664   //
1665   hs = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix,fRecPrimTrackMatrix,"zv_pt_eta_track_corr_matrix");
1666   aFolderObj->Add(hs);
1667
1668   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(1,2),fRecPrimTrackMatrix->Projection(1,2),"eta_pt_track_corr_matrix");
1669   aFolderObj->Add(h2D);
1670
1671   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(1,0),fRecPrimTrackMatrix->Projection(1,0),"pt_zv_track_corr_matrix");
1672   aFolderObj->Add(h2D);
1673
1674   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(2,0),fRecPrimTrackMatrix->Projection(2,0),"zv_eta_track_corr_matrix");
1675   aFolderObj->Add(h2D);
1676
1677
1678   h = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(0),fRecPrimTrackMatrix->Projection(0),"zv_track_corr_matrix");
1679   aFolderObj->Add(h);
1680
1681   h = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(1),fRecPrimTrackMatrix->Projection(1),"pt_track_corr_matrix");
1682   aFolderObj->Add(h);
1683
1684   // efficiency
1685
1686   h = AlidNdPtHelper::GenerateCorrMatrix(fRecPrimTrackMatrix->Projection(1), fGenPrimTrackMatrix->Projection(1),"pt_track_eff_matrix");
1687   aFolderObj->Add(h);
1688
1689   h = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(2),fRecPrimTrackMatrix->Projection(2),"eta_track_corr_matrix");
1690   aFolderObj->Add(h);
1691
1692   //
1693   // secondary track contamination correction
1694   //
1695   //hs = AlidNdPtHelper::GenerateContCorrMatrix(fRecSecTrackMatrix,fRecTrackMatrix,"zv_pt_eta_track_cont_matrix");
1696   hs = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix,fRecTrackMatrix,"zv_pt_eta_track_cont_matrix");
1697   aFolderObj->Add(hs);
1698
1699   h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(1,2),fRecTrackMatrix->Projection(1,2),"eta_pt_track_cont_matrix");
1700   aFolderObj->Add(h2D);
1701
1702   h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(1,0),fRecTrackMatrix->Projection(1,0),"pt_zv_track_cont_matrix");
1703   aFolderObj->Add(h2D);
1704
1705   h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(2,0),fRecTrackMatrix->Projection(2,0),"zv_eta_track_cont_matrix");
1706   aFolderObj->Add(h2D);
1707
1708   h = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(0),fRecTrackMatrix->Projection(0),"zv_track_cont_matrix");
1709   aFolderObj->Add(h);
1710
1711
1712   h = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(1),fRecTrackMatrix->Projection(1),"pt_track_cont_matrix");
1713   aFolderObj->Add(h);
1714
1715   h = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(2),fRecTrackMatrix->Projection(2),"eta_track_cont_matrix");
1716   aFolderObj->Add(h);
1717
1718   //
1719   // multiple track reconstruction correction
1720   //
1721   //hs = AlidNdPtHelper::GenerateContCorrMatrix(fRecMultTrackMatrix,fRecTrackMatrix,"zv_pt_eta_mult_track_cont_matrix");
1722   hs = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix,fRecTrackMatrix,"zv_pt_eta_mult_track_cont_matrix");
1723   aFolderObj->Add(hs);
1724
1725   h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(1,2),fRecTrackMatrix->Projection(1,2),"eta_pt_mult_track_cont_matrix");
1726   aFolderObj->Add(h2D);
1727
1728   h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(1,0),fRecTrackMatrix->Projection(1,0),"pt_zv_mult_track_cont_matrix");
1729   aFolderObj->Add(h2D);
1730
1731   h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(2,0),fRecTrackMatrix->Projection(2,0),"zv_eta_mult_track_cont_matrix");
1732   aFolderObj->Add(h2D);
1733
1734   h = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(0),fRecTrackMatrix->Projection(0),"zv_mult_track_cont_matrix");
1735   aFolderObj->Add(h);
1736
1737   h = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(1),fRecTrackMatrix->Projection(1),"pt_mult_track_cont_matrix");
1738   aFolderObj->Add(h);
1739
1740   h = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(2),fRecTrackMatrix->Projection(2),"eta_mult_track_cont_matrix");
1741   aFolderObj->Add(h);
1742
1743   //
1744   // Control histograms
1745   //
1746
1747   if(fHistogramsOn) {
1748
1749   // Efficiency electrons, muons, pions, kaons, protons, all
1750   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(1,1);
1751   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(1,1);
1752   h1 = fMCPrimTrackHist1[1]->Projection(0);
1753   h2 = fMCPrimTrackHist1[2]->Projection(0);
1754   h2c = (TH1D *)h2->Clone();
1755   h2c->Divide(h1);
1756   h2c->SetName("eff_pt_electrons");
1757   aFolderObj->Add(h2c);
1758
1759   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(2,2);
1760   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(2,2);
1761   h1 = fMCPrimTrackHist1[1]->Projection(0);
1762   h2 = fMCPrimTrackHist1[2]->Projection(0);
1763   h2c = (TH1D *)h2->Clone();
1764   h2c->Divide(h1);
1765   h2c->SetName("eff_pt_muons");
1766   aFolderObj->Add(h2c);
1767
1768   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(3,3);
1769   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(3,3);
1770   h1 = fMCPrimTrackHist1[1]->Projection(0);
1771   h2 = fMCPrimTrackHist1[2]->Projection(0);
1772   h2c = (TH1D *)h2->Clone();
1773   h2c->Divide(h1);
1774   h2c->SetName("eff_pt_pions");
1775   aFolderObj->Add(h2c);
1776
1777   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(4,4);
1778   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(4,4);
1779   h1 = fMCPrimTrackHist1[1]->Projection(0);
1780   h2 = fMCPrimTrackHist1[2]->Projection(0);
1781   h2c = (TH1D *)h2->Clone();
1782   h2c->Divide(h1);
1783   h2c->SetName("eff_pt_kaons");
1784   aFolderObj->Add(h2c);
1785
1786   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(5,5);
1787   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(5,5);
1788   h1 = fMCPrimTrackHist1[1]->Projection(0);
1789   h2 = fMCPrimTrackHist1[2]->Projection(0);
1790   h2c = (TH1D *)h2->Clone();
1791   h2c->Divide(h1);
1792   h2c->SetName("eff_pt_protons");
1793   aFolderObj->Add(h2c);
1794
1795   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(1,5);
1796   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(1,5);
1797   h1 = fMCPrimTrackHist1[1]->Projection(0);
1798   h2 = fMCPrimTrackHist1[2]->Projection(0);
1799   h2c = (TH1D *)h2->Clone();
1800   h2c->Divide(h1);
1801   h2c->SetName("eff_pt_selected");
1802   aFolderObj->Add(h2c);
1803
1804   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(1,6);
1805   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(1,6);
1806   h1 = fMCPrimTrackHist1[1]->Projection(0);
1807   h2 = fMCPrimTrackHist1[2]->Projection(0);
1808   h2c = (TH1D *)h2->Clone();
1809   h2c->Divide(h1);
1810   h2c->SetName("eff_pt_all");
1811   aFolderObj->Add(h2c);
1812
1813   fMCPrimTrackHist1[1]->GetAxis(1)->SetRange(1,fMCPrimTrackHist1[1]->GetAxis(1)->GetNbins());
1814   fMCPrimTrackHist1[2]->GetAxis(1)->SetRange(1,fMCPrimTrackHist1[2]->GetAxis(1)->GetNbins());
1815
1816   //  pt spetra
1817   // - rec, primaries, secondaries
1818   // - primaries (pid)
1819   // - secondaries (pid)
1820   // - secondaries (mech)
1821   // - secondaries (mother)
1822   //
1823
1824   TH1D *mcPtAccall = fMCTrackHist1[1]->Projection(0);
1825   mcPtAccall->SetName("mc_pt_acc_all");
1826   aFolderObj->Add(mcPtAccall);
1827
1828   TH1D *mcPtAccprim = fMCPrimTrackHist1[1]->Projection(0);
1829   mcPtAccprim->SetName("mc_pt_acc_prim");
1830   aFolderObj->Add(mcPtAccprim);
1831
1832   TH1D *mcPtRecall = fMCTrackHist1[2]->Projection(0);
1833   mcPtRecall->SetName("mc_pt_rec_all");
1834   aFolderObj->Add(mcPtRecall);
1835
1836   TH1D *mcPtRecprim = fMCPrimTrackHist1[2]->Projection(0);
1837   mcPtRecprim->SetName("mc_pt_rec_prim");
1838   aFolderObj->Add(mcPtRecprim);
1839
1840   TH1D *mcPtRecsec = fMCSecTrackHist1[2]->Projection(0);
1841   mcPtRecsec->SetName("mc_pt_rec_sec");
1842   aFolderObj->Add(mcPtRecsec);
1843
1844   for(Int_t i = 0; i<6; i++)
1845   {
1846     snprintf(name,256,"mc_pt_rec_prim_pid_%d",i);
1847     fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(i+1,i+1);
1848     h = fMCPrimTrackHist1[2]->Projection(0);
1849     h->SetName(name);
1850     aFolderObj->Add(h);
1851
1852     snprintf(name,256,"mc_pt_rec_sec_pid_%d",i);
1853     fMCSecTrackHist1[2]->GetAxis(2)->SetRange(i+1,i+1);
1854     h = fMCSecTrackHist1[2]->Projection(0);
1855     h->SetName(name);
1856     aFolderObj->Add(h);
1857
1858     // production mechanisms for given pid
1859     fMCSecTrackHist1[2]->GetAxis(2)->SetRange(i+1,i+1);
1860
1861     for(Int_t j=0; j<20; j++) {
1862       if(j == 4) {
1863         // decay
1864
1865         snprintf(name,256,"mc_pt_rec_sec_pid_%d_decay",i);
1866         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1867         h = fMCSecTrackHist1[2]->Projection(0);
1868         h->SetName(name);
1869         aFolderObj->Add(h);
1870
1871         snprintf(name,256,"mc_eta_rec_sec_pid_%d_decay",i);
1872         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1873         h = fMCSecTrackHist1[2]->Projection(1);
1874         h->SetName(name);
1875         aFolderObj->Add(h);
1876
1877         snprintf(name,256,"mc_mother_rec_sec_pid_%d_decay",i);
1878         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1879         h = fMCSecTrackHist1[2]->Projection(4);
1880         h->SetName(name);
1881         aFolderObj->Add(h);
1882
1883       } else if (j == 5) {
1884        // conversion
1885
1886         snprintf(name,256,"mc_pt_rec_sec_pid_%d_conv",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         snprintf(name,256,"mc_eta_rec_sec_pid_%d_conv",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         snprintf(name,256,"mc_mother_rec_sec_pid_%d_conv",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 == 13) {
1905        // mat
1906
1907         snprintf(name,256,"mc_pt_rec_sec_pid_%d_mat",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         snprintf(name,256,"mc_eta_rec_sec_pid_%d_mat",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         snprintf(name,256,"mc_eta_mother_rec_sec_pid_%d_mat",i);
1920         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1921         h = fMCSecTrackHist1[2]->Projection(4,1);
1922         h->SetName(name);
1923         aFolderObj->Add(h);
1924
1925         snprintf(name,256,"mc_mother_rec_sec_pid_%d_mat",i);
1926         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1927         h = fMCSecTrackHist1[2]->Projection(4);
1928         h->SetName(name);
1929         aFolderObj->Add(h);
1930
1931         snprintf(name,256,"mc_pt_mother_rec_sec_pid_%d_mat",i);
1932         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1933         h = fMCSecTrackHist1[2]->Projection(4,0);
1934         h->SetName(name);
1935         aFolderObj->Add(h);
1936
1937      } else {
1938        continue;
1939      }
1940    }
1941
1942   }
1943   } // end fHistogramOn
1944
1945   //
1946   //  resolution histograms
1947   //  only for reconstructed tracks
1948   //
1949
1950   TH2F *h2F=0;
1951   TCanvas * c = new TCanvas("resol","resol");
1952   c->cd();
1953
1954   //
1955   fRecMCTrackHist1->GetAxis(1)->SetRangeUser(-0.8,0.79);
1956
1957   h2F = (TH2F*)fRecMCTrackHist1->Projection(2,0);
1958   h = AlidNdPtHelper::MakeResol(h2F,1,0,kTRUE,10);
1959   h->SetXTitle("p_{tmc} (GeV/c)");
1960   h->SetYTitle("(p_{t}-p_{tmc})/p_{tmc} resolution");
1961   h->Draw();
1962   h->SetName("pt_resolution_vs_mcpt");
1963   aFolderObj->Add(h);
1964
1965   h2F = (TH2F*)fRecMCTrackHist1->Projection(2,0);
1966   h = AlidNdPtHelper::MakeResol(h2F,1,1,kTRUE,10);
1967   h->SetXTitle("p_{tmc} (GeV/c)");
1968   h->SetYTitle("(p_{t}-p_{tmc})/p_{tmc} mean");
1969   h->Draw();
1970   h->SetName("dpt_mean_vs_mcpt");
1971   aFolderObj->Add(h);
1972
1973   //
1974   h2F = (TH2F*)fRecMCTrackHist1->Projection(3,0);
1975   h = AlidNdPtHelper::MakeResol(h2F,1,0,kTRUE,10);
1976   h->SetXTitle("p_{tmc} (GeV/c)");
1977   h->SetYTitle("(#eta-#eta_{mc}) resolution");
1978   h->Draw();
1979   h->SetName("eta_resolution_vs_mcpt");
1980   aFolderObj->Add(h);
1981
1982   h2F = (TH2F*)fRecMCTrackHist1->Projection(3,0);
1983   h = AlidNdPtHelper::MakeResol(h2F,1,1,kTRUE,10);
1984   h->SetXTitle("p_{tmc} (GeV/c)");
1985   h->SetYTitle("(#eta-mc#eta) mean");
1986   h->Draw();
1987   h->SetName("deta_mean_vs_mcpt");
1988   aFolderObj->Add(h);
1989
1990   //
1991   fRecMCTrackHist1->GetAxis(1)->SetRange(1,fRecMCTrackHist1->GetAxis(1)->GetNbins());
1992
1993   h2F = (TH2F*)fRecMCTrackHist1->Projection(2,1);
1994   h = AlidNdPtHelper::MakeResol(h2F,1,0,kTRUE,10);
1995   h->SetXTitle("#eta_{mc}");
1996   h->SetYTitle("(p_{t}-p_{tmc})/p_{tmc} resolution");
1997   h->Draw();
1998   h->SetName("pt_resolution_vs_mceta");
1999   aFolderObj->Add(h);
2000
2001   h2F = (TH2F*)fRecMCTrackHist1->Projection(2,1);
2002   h = AlidNdPtHelper::MakeResol(h2F,1,1,kTRUE,10);
2003   h->SetXTitle("#eta_{mc}");
2004   h->SetYTitle("(p_{t}-p_{tmc})/p_{tmc} mean");
2005   h->Draw();
2006   h->SetName("dpt_mean_vs_mceta");
2007   aFolderObj->Add(h);
2008
2009   //
2010   h2F = (TH2F*)fRecMCTrackHist1->Projection(3,1);
2011   h = AlidNdPtHelper::MakeResol(h2F,1,0,kTRUE,10);
2012   h->SetXTitle("#eta_{mc}");
2013   h->SetYTitle("(#eta-#eta_{mc}) resolution");
2014   h->Draw();
2015   h->SetName("eta_resolution_vs_mceta");
2016   aFolderObj->Add(h);
2017
2018   h2F = (TH2F*)fRecMCTrackHist1->Projection(3,1);
2019   h = AlidNdPtHelper::MakeResol(h2F,1,1,kTRUE,10);
2020   h->SetXTitle("#eta_{mc}");
2021   h->SetYTitle("(#eta-mc#eta) mean");
2022   h->Draw();
2023   h->SetName("deta_mean_vs_mceta");
2024   aFolderObj->Add(h);
2025
2026   fRecMCTrackHist1->GetAxis(0)->SetRange(1,fRecMCTrackHist1->GetAxis(0)->GetNbins());
2027
2028   } // end use MC info
2029
2030   // export objects to analysis folder
2031   fAnalysisFolder = ExportToFolder(aFolderObj);
2032
2033   // delete only TObjArray
2034   if(aFolderObj) delete aFolderObj;
2035 }
2036
2037 //_____________________________________________________________________________
2038 TFolder* AlidNdPtAnalysisPbPb::ExportToFolder(TObjArray * const array)
2039 {
2040   // recreate folder avery time and export objects to new one
2041   //
2042   AlidNdPtAnalysisPbPb * comp=this;
2043   TFolder *folder = comp->GetAnalysisFolder();
2044
2045   TString name, title;
2046   TFolder *newFolder = 0;
2047   Int_t i = 0;
2048   Int_t size = array->GetSize();
2049
2050   if(folder) {
2051      // get name and title from old folder
2052      name = folder->GetName();
2053      title = folder->GetTitle();
2054
2055          // delete old one
2056      delete folder;
2057
2058          // create new one
2059      newFolder = CreateFolder(name.Data(),title.Data());
2060      newFolder->SetOwner();
2061
2062          // add objects to folder
2063      while(i < size) {
2064            newFolder->Add(array->At(i));
2065            i++;
2066          }
2067   }
2068
2069 return newFolder;
2070 }
2071
2072 //_____________________________________________________________________________
2073 TFolder* AlidNdPtAnalysisPbPb::CreateFolder(TString name,TString title) {
2074 // create folder for analysed histograms
2075 //
2076 TFolder *folder = 0;
2077   folder = new TFolder(name.Data(),title.Data());
2078
2079   return folder;
2080 }