PWGLF/SPECTRA/ChargedHadrons/dNdPt/AlidNdPtAnalysis.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 // AlidNdPtAnalysis 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: 2013-06-13 by M.Knichel
  29 //
  30 // meaning of different multiplicities:
  31 // multRec      : number of reconstructed tracks, after AcceptanceCuts and TrackCuts
  32 // multRecMult  : number of reconstructed tracks, after MultAcceptanceCuts and MultTrackCuts
  33 // multMB       : number of contributers to vertex
  34 // multTrueMC   : MC true mult, after MultAcceptanceCuts
  35 // mutlAll      : number of ESD tracks
  36 // mutlAcc      : number of ESD tracks after AcceptanceCuts
  37 //------------------------------------------------------------------------------
  38
  39 #include "TH1.h"
  40 #include "TH2.h"
  41 #include "TCanvas.h"
  42 #include "THnSparse.h"
  43
  44 #include "AliHeader.h"
  45 #include "AliInputEventHandler.h"
  46 #include "AliAnalysisManager.h"
  47 #include "AliGenEventHeader.h"
  48 #include "AliStack.h"
  49 #include "AliESDEvent.h"
  50 #include "AliMCEvent.h"
  51 #include "AliESDtrackCuts.h"
  52 #include "AliLog.h"
  53 #include "AliMultiplicity.h"
  54 #include "AliTracker.h"
  55
  56 #include "AlidNdPtEventCuts.h"
  57 #include "AlidNdPtAcceptanceCuts.h"
  58 #include "AliPhysicsSelection.h"
  59 #include "AliTriggerAnalysis.h"
  60
  61 #include "AliPWG0Helper.h"
  62 #include "AlidNdPtHelper.h"
  63 #include "AlidNdPtAnalysis.h"
  64
  65 using namespace std;
  66
  67 ClassImp(AlidNdPtAnalysis)
  68
  69 //_____________________________________________________________________________
  70   AlidNdPtAnalysis::AlidNdPtAnalysis(): AlidNdPt(),
  71   fAnalysisFolder(0),
  72   fHistogramsOn(kFALSE),
  73
  74   // event multiplicity correlation matrix
  75   fEventMultCorrelationMatrix(0),
  76
  77   // rec. track pt vs true track pt correlation matrix
  78   fTrackPtCorrelationMatrix(0),
  79
  80   // event level correction
  81   fGenEventMatrix(0),
  82   fGenSDEventMatrix(0),
  83   fGenDDEventMatrix(0),
  84   fGenNDEventMatrix(0),
  85   fGenNSDEventMatrix(0),
  86
  87   fTriggerEventMatrix(0),
  88   fTriggerSDEventMatrix(0),
  89   fTriggerDDEventMatrix(0),
  90   fTriggerNDEventMatrix(0),
  91   fTriggerNSDEventMatrix(0),
  92
  93   fRecEventMatrix(0),
  94   fRecSDEventMatrix(0),
  95   fRecDDEventMatrix(0),
  96   fRecNDEventMatrix(0),
  97   fRecNSDEventMatrix(0),
  98
  99   //
 100   // track-event level correction
 101   //
 102   fGenTrackEventMatrix(0),
 103   fGenTrackSDEventMatrix(0),
 104   fGenTrackDDEventMatrix(0),
 105   fGenTrackNDEventMatrix(0),
 106   fGenTrackNSDEventMatrix(0),
 107
 108   fTriggerTrackEventMatrix(0),
 109   fTriggerTrackSDEventMatrix(0),
 110   fTriggerTrackDDEventMatrix(0),
 111   fTriggerTrackNDEventMatrix(0),
 112   fTriggerTrackNSDEventMatrix(0),
 113
 114   fRecTrackEventMatrix(0),
 115   fRecTrackSDEventMatrix(0),
 116   fRecTrackDDEventMatrix(0),
 117   fRecTrackNDEventMatrix(0),
 118   fRecTrackNSDEventMatrix(0),
 119
 120   // track rec. efficiency correction (fRecPrimTrackMatrix / fGenPrimTrackMatrix)
 121   fGenTrackMatrix(0),
 122   fGenPrimTrackMatrix(0),
 123   fRecPrimTrackMatrix(0),
 124
 125   // secondary track contamination correction (fRecSecTrackMatrix / fRecTrackMatrix)
 126   fRecTrackMatrix(0),
 127   fRecSecTrackMatrix(0),
 128
 129   // multiple rec. track contamination corrections (fRecMultTrackMatrix / fRecTrackMatrix)
 130   fRecMultTrackMatrix(0),
 131
 132   // event control histograms
 133   fMCEventHist1(0),
 134   fRecEventHist1(0),
 135   fRecEventHist2(0),
 136   fRecMCEventHist1(0),
 137   fRecMCEventHist2(0),
 138   fRecMCEventHist3(0),
 139
 140   // rec. pt and eta resolution w.r.t MC
 141   fRecMCTrackHist1(0),
 142
 143   //multple reconstructed tracks
 144   fMCMultRecTrackHist1(0),
 145
 146   // rec. track control histograms
 147   fRecTrackHist2(0),
 148
 149   // Generic histograms to be corrected
 150   fRecEventHist(0),
 151   fRecTrackHist(0),
 152   fEventCount(0),
 153   fMCPrimTrackHist(0),
 154
 155   // Candle event histogram
 156   fRecCandleEventMatrix(0),
 157
 158   fMultNbins(0),
 159   fPtNbins(0),
 160   fPtCorrNbins(0),
 161   fEtaNbins(0),
 162   fZvNbins(0),
 163   fBinsMult(0),
 164   fBinsPt(0),
 165   fBinsPtCorr(0),
 166   fBinsEta(0),
 167   fBinsZv(0),
 168
 169   fIsInit(kFALSE)
 170
 171 {
 172   // default constructor
 173   for(Int_t i=0; i<AlidNdPtHelper::kCutSteps; i++) {
 174     fMCTrackHist1[i]=0;
 175     fMCPrimTrackHist1[i]=0;
 176     fMCPrimTrackHist2[i]=0;
 177     fMCSecTrackHist1[i]=0;
 178     fRecTrackHist1[i]=0;
 179     fRecTrackMultHist1[i]=0;
 180   }
 181   //Init();
 182 }
 183
 184 //_____________________________________________________________________________
 185 AlidNdPtAnalysis::AlidNdPtAnalysis(Char_t* name, Char_t* title): AlidNdPt(name,title),
 186   fAnalysisFolder(0),
 187   fHistogramsOn(kFALSE),
 188
 189   // event multiplicity correlation matrix
 190   fEventMultCorrelationMatrix(0),
 191
 192   // rec. track pt vs true track pt correlation matrix
 193   fTrackPtCorrelationMatrix(0),
 194
 195   // event level correction
 196   fGenEventMatrix(0),
 197   fGenSDEventMatrix(0),
 198   fGenDDEventMatrix(0),
 199   fGenNDEventMatrix(0),
 200   fGenNSDEventMatrix(0),
 201
 202   fTriggerEventMatrix(0),
 203   fTriggerSDEventMatrix(0),
 204   fTriggerDDEventMatrix(0),
 205   fTriggerNDEventMatrix(0),
 206   fTriggerNSDEventMatrix(0),
 207
 208   fRecEventMatrix(0),
 209   fRecSDEventMatrix(0),
 210   fRecDDEventMatrix(0),
 211   fRecNDEventMatrix(0),
 212   fRecNSDEventMatrix(0),
 213
 214   //
 215   // track-event level correction
 216   //
 217   fGenTrackEventMatrix(0),
 218   fGenTrackSDEventMatrix(0),
 219   fGenTrackDDEventMatrix(0),
 220   fGenTrackNDEventMatrix(0),
 221   fGenTrackNSDEventMatrix(0),
 222
 223   fTriggerTrackEventMatrix(0),
 224   fTriggerTrackSDEventMatrix(0),
 225   fTriggerTrackDDEventMatrix(0),
 226   fTriggerTrackNDEventMatrix(0),
 227   fTriggerTrackNSDEventMatrix(0),
 228
 229   fRecTrackEventMatrix(0),
 230   fRecTrackSDEventMatrix(0),
 231   fRecTrackDDEventMatrix(0),
 232   fRecTrackNDEventMatrix(0),
 233   fRecTrackNSDEventMatrix(0),
 234
 235   // track rec. efficiency correction (fRecPrimTrackMatrix / fGenPrimTrackMatrix)
 236   fGenTrackMatrix(0),
 237   fGenPrimTrackMatrix(0),
 238   fRecPrimTrackMatrix(0),
 239
 240   // secondary track contamination correction (fRecSecTrackMatrix / fRecTrackMatrix)
 241   fRecTrackMatrix(0),
 242   fRecSecTrackMatrix(0),
 243
 244   // multiple rec. track contamination corrections (fRecMultTrackMatrix / fRecTrackMatrix)
 245   fRecMultTrackMatrix(0),
 246
 247   // event control histograms
 248   fMCEventHist1(0),
 249   fRecEventHist1(0),
 250   fRecEventHist2(0),
 251   fRecMCEventHist1(0),
 252   fRecMCEventHist2(0),
 253   fRecMCEventHist3(0),
 254
 255   // rec. pt and eta resolution w.r.t MC
 256   fRecMCTrackHist1(0),
 257
 258   //multple reconstructed tracks
 259   fMCMultRecTrackHist1(0),
 260
 261   // rec. track control histograms
 262   fRecTrackHist2(0),
 263
 264   // Generic histograms to be corrected
 265   fRecEventHist(0),
 266   fRecTrackHist(0),
 267   fEventCount(0),
 268   fMCPrimTrackHist(0),
 269
 270   // Candle event histogram
 271   fRecCandleEventMatrix(0),
 272
 273   fMultNbins(0),
 274   fPtNbins(0),
 275   fPtCorrNbins(0),
 276   fEtaNbins(0),
 277   fZvNbins(0),
 278   fBinsMult(0),
 279   fBinsPt(0),
 280   fBinsPtCorr(0),
 281   fBinsEta(0),
 282   fBinsZv(0),
 283
 284   fIsInit(kFALSE)
 285
 286 {
 287   //
 288   // constructor
 289   //
 290   for(Int_t i=0; i<AlidNdPtHelper::kCutSteps; i++) {
 291     fMCTrackHist1[i]=0;
 292     fMCPrimTrackHist1[i]=0;
 293     fMCPrimTrackHist2[i]=0;
 294     fMCSecTrackHist1[i]=0;
 295     fRecTrackHist1[i]=0;
 296     fRecTrackMultHist1[i]=0;
 297   }
 298
 299   //Init();
 300 }
 301
 302 //_____________________________________________________________________________
 303 AlidNdPtAnalysis::~AlidNdPtAnalysis() {
 304   //
 305   // destructor
 306   //
 307   if(fEventMultCorrelationMatrix) delete fEventMultCorrelationMatrix; fEventMultCorrelationMatrix=0;
 308   //
 309   if(fTrackPtCorrelationMatrix) delete fTrackPtCorrelationMatrix; fTrackPtCorrelationMatrix=0;
 310   //
 311   if(fGenEventMatrix) delete fGenEventMatrix; fGenEventMatrix=0;
 312   if(fGenSDEventMatrix) delete fGenSDEventMatrix; fGenSDEventMatrix=0;
 313   if(fGenDDEventMatrix) delete fGenDDEventMatrix; fGenDDEventMatrix=0;
 314   if(fGenNDEventMatrix) delete fGenNDEventMatrix; fGenNDEventMatrix=0;
 315   if(fGenNSDEventMatrix) delete fGenNSDEventMatrix; fGenNSDEventMatrix=0;
 316
 317   if(fTriggerEventMatrix) delete fTriggerEventMatrix; fTriggerEventMatrix=0;
 318   if(fTriggerSDEventMatrix) delete fTriggerSDEventMatrix; fTriggerSDEventMatrix=0;
 319   if(fTriggerDDEventMatrix) delete fTriggerDDEventMatrix; fTriggerDDEventMatrix=0;
 320   if(fTriggerNDEventMatrix) delete fTriggerNDEventMatrix; fTriggerNDEventMatrix=0;
 321   if(fTriggerNSDEventMatrix) delete fTriggerNSDEventMatrix; fTriggerNSDEventMatrix=0;
 322
 323   if(fRecEventMatrix) delete fRecEventMatrix; fRecEventMatrix=0;
 324   if(fRecSDEventMatrix) delete fRecSDEventMatrix; fRecSDEventMatrix=0;
 325   if(fRecDDEventMatrix) delete fRecDDEventMatrix; fRecDDEventMatrix=0;
 326   if(fRecNDEventMatrix) delete fRecNDEventMatrix; fRecNDEventMatrix=0;
 327   if(fRecNSDEventMatrix) delete fRecNSDEventMatrix; fRecNSDEventMatrix=0;
 328
 329   if(fRecCandleEventMatrix) delete fRecCandleEventMatrix; fRecCandleEventMatrix=0;
 330   //
 331   if(fGenTrackEventMatrix) delete fGenTrackEventMatrix; fGenTrackEventMatrix=0;
 332   if(fGenTrackSDEventMatrix) delete fGenTrackSDEventMatrix; fGenTrackSDEventMatrix=0;
 333   if(fGenTrackDDEventMatrix) delete fGenTrackDDEventMatrix; fGenTrackDDEventMatrix=0;
 334   if(fGenTrackNDEventMatrix) delete fGenTrackNDEventMatrix; fGenTrackNDEventMatrix=0;
 335   if(fGenTrackNSDEventMatrix) delete fGenTrackNSDEventMatrix; fGenTrackNSDEventMatrix=0;
 336
 337   if(fTriggerTrackEventMatrix) delete fTriggerTrackEventMatrix; fTriggerTrackEventMatrix=0;
 338   if(fTriggerTrackSDEventMatrix) delete fTriggerTrackSDEventMatrix; fTriggerTrackSDEventMatrix=0;
 339   if(fTriggerTrackDDEventMatrix) delete fTriggerTrackDDEventMatrix; fTriggerTrackDDEventMatrix=0;
 340   if(fTriggerTrackNDEventMatrix) delete fTriggerTrackNDEventMatrix; fTriggerTrackNDEventMatrix=0;
 341   if(fTriggerTrackNSDEventMatrix) delete fTriggerTrackNSDEventMatrix; fTriggerTrackNSDEventMatrix=0;
 342
 343   if(fRecTrackEventMatrix) delete fRecTrackEventMatrix; fRecTrackEventMatrix=0;
 344   if(fRecTrackSDEventMatrix) delete fRecTrackSDEventMatrix; fRecTrackSDEventMatrix=0;
 345   if(fRecTrackDDEventMatrix) delete fRecTrackDDEventMatrix; fRecTrackDDEventMatrix=0;
 346   if(fRecTrackNDEventMatrix) delete fRecTrackNDEventMatrix; fRecTrackNDEventMatrix=0;
 347   if(fRecTrackNSDEventMatrix) delete fRecTrackNSDEventMatrix; fRecTrackNSDEventMatrix=0;
 348
 349   //
 350   if(fGenTrackMatrix) delete fGenTrackMatrix; fGenTrackMatrix=0;
 351   if(fGenPrimTrackMatrix) delete fGenPrimTrackMatrix; fGenPrimTrackMatrix=0;
 352   if(fRecPrimTrackMatrix) delete fRecPrimTrackMatrix; fRecPrimTrackMatrix=0;
 353   //
 354   if(fRecTrackMatrix) delete fRecTrackMatrix; fRecTrackMatrix=0;
 355   if(fRecSecTrackMatrix) delete fRecSecTrackMatrix; fRecSecTrackMatrix=0;
 356   //
 357   if(fRecMultTrackMatrix) delete fRecMultTrackMatrix; fRecMultTrackMatrix=0;
 358   //
 359   // Control histograms
 360   //
 361   if(fMCEventHist1) delete fMCEventHist1; fMCEventHist1=0;
 362   if(fRecEventHist1) delete fRecEventHist1; fRecEventHist1=0;
 363   if(fRecEventHist2) delete fRecEventHist2; fRecEventHist2=0;
 364   if(fRecMCEventHist1) delete fRecMCEventHist1; fRecMCEventHist1=0;
 365   if(fRecMCEventHist2) delete fRecMCEventHist2; fRecMCEventHist2=0;
 366   if(fRecMCEventHist3) delete fRecMCEventHist3; fRecMCEventHist3=0;
 367   //
 368   for(Int_t i=0; i<AlidNdPtHelper::kCutSteps; i++) {
 369     if(fMCTrackHist1[i]) delete fMCTrackHist1[i]; fMCTrackHist1[i]=0;
 370     if(fMCPrimTrackHist1[i]) delete fMCPrimTrackHist1[i]; fMCPrimTrackHist1[i]=0;
 371     if(fMCPrimTrackHist2[i]) delete fMCPrimTrackHist2[i]; fMCPrimTrackHist2[i]=0;
 372     if(fMCSecTrackHist1[i]) delete fMCSecTrackHist1[i]; fMCSecTrackHist1[i]=0;
 373     if(fRecTrackHist1[i]) delete fRecTrackHist1[i]; fRecTrackHist1[i]=0;
 374     if(fRecTrackMultHist1[i]) delete fRecTrackMultHist1[i]; fRecTrackMultHist1[i]=0;
 375   }
 376   if(fRecMCTrackHist1) delete fRecMCTrackHist1; fRecMCTrackHist1=0;
 377   if(fMCMultRecTrackHist1) delete fMCMultRecTrackHist1; fMCMultRecTrackHist1=0;
 378   if(fRecTrackHist2) delete fRecTrackHist2; fRecTrackHist2=0;
 379
 380   //
 381   if(fRecEventHist) delete fRecEventHist; fRecEventHist=0;
 382   if(fRecTrackHist) delete fRecTrackHist; fRecTrackHist=0;
 383   if(fEventCount) delete fEventCount; fEventCount=0;
 384   if(fMCPrimTrackHist) delete fMCPrimTrackHist; fMCPrimTrackHist=0;
 385
 386   //
 387   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
 388
 389   if (fBinsMult) delete[] fBinsMult; fBinsMult=0;
 390   if (fBinsPt) delete[] fBinsPt; fBinsPt=0;
 391   if (fBinsPtCorr) delete[] fBinsPtCorr; fBinsPtCorr=0;
 392   if (fBinsEta) delete[] fBinsEta; fBinsEta=0;
 393   if (fBinsZv) delete[] fBinsMult; fBinsZv=0;
 394 }
 395
 396 //_____________________________________________________________________________
 397 void AlidNdPtAnalysis::Init()
 398 {
 399     //define default binning
 400     Double_t binsMultDefault[15] = {-0.5, 0.5 , 1.5 , 2.5 , 3.5 , 4.5 , 5.5 , 6.5 , 7.5 , 8.5, 9.5, 10.5, 20.5, 50.5, 150.5};
 401     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};
 402     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};
 403     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};
 404     Double_t binsZvDefault[13] = {-30.,-25.,-20.,-15.,-10.,-5.,0.,5.,10.,15.,20.,25.,30.};
 405
 406    // if no binning is set, use the default
 407    if (!fBinsMult)   { SetBinsMult(14,binsMultDefault); }
 408    if (!fBinsPt)     { SetBinsPt(68,binsPtDefault); }
 409    if (!fBinsPtCorr) { SetBinsPtCorr(36,binsPtCorrDefault); }
 410    if (!fBinsEta)    { SetBinsEta(30,binsEtaDefault); }
 411    if (!fBinsZv)     { SetBinsZv(12,binsZvDefault); }
 412
 413   //Int_t binsTrackMatrix[3]={zvNbins,ptNbins,etaNbins};
 414   Int_t binsTrackEventCorrMatrix[4]={fZvNbins,fPtCorrNbins,fEtaNbins,fMultNbins};
 415
 416   //
 417   // Generic histograms to be corrected
 418   //
 419   Int_t binsEventHist[2]={fZvNbins,fMultNbins};
 420   //Double_t minEventHist[2]={-fBinsZv[0],fBinsMult[0]};
 421   //Double_t maxEventHist[2]={fBinsZv[fZvNbins],fBinsMult[fMultNbins]};
 422
 423   fRecEventHist = new THnSparseF("fRecEventHist","Zv:multMB",2,binsEventHist); //,minEventHist,maxEventHist);
 424   fRecEventHist->SetBinEdges(0,fBinsZv);
 425   fRecEventHist->SetBinEdges(1,fBinsMult);
 426   fRecEventHist->GetAxis(0)->SetTitle("Zv (cm)");
 427   fRecEventHist->GetAxis(1)->SetTitle("multiplicity MB");
 428   fRecEventHist->Sumw2();
 429
 430   //
 431   Int_t binsTrackHist[4]={fZvNbins,fPtNbins,fEtaNbins,fMultNbins};
 432  // Double_t minTrackHist[4]={-25.,0.,-1.5,-0.5};
 433  // Double_t maxTrackHist[4]={25.,50.,1.5,149.5};
 434
 435   fRecTrackHist = new THnSparseF("fRecTrackHist","Zv:pT:eta:multRecMult",4,binsTrackHist); //,minTrackHist,maxTrackHist);
 436   fRecTrackHist->SetBinEdges(0,fBinsZv);
 437   fRecTrackHist->SetBinEdges(1,fBinsPt);
 438   fRecTrackHist->SetBinEdges(2,fBinsEta);
 439   fRecTrackHist->SetBinEdges(3,fBinsMult);
 440   fRecTrackHist->GetAxis(0)->SetTitle("Zv (cm)");
 441   fRecTrackHist->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
 442   fRecTrackHist->GetAxis(2)->SetTitle("#eta");
 443   fRecTrackHist->GetAxis(3)->SetTitle("multiplicity (multCuts)");
 444   fRecTrackHist->Sumw2();
 445
 446     Int_t binsMCPrimTrackHist[4]={fZvNbins,fPtNbins,fEtaNbins,fMultNbins};
 447  // Double_t minTrackHist[4]={-25.,0.,-1.5,-0.5};
 448  // Double_t maxTrackHist[4]={25.,50.,1.5,149.5};
 449
 450   fMCPrimTrackHist = new THnSparseF("fMCPrimTrackHist","Zv:mcpT:mceta:multTrueMC",4,binsMCPrimTrackHist);
 451   fMCPrimTrackHist->SetBinEdges(0,fBinsZv);
 452   fMCPrimTrackHist->SetBinEdges(1,fBinsPt);
 453   fMCPrimTrackHist->SetBinEdges(2,fBinsEta);
 454   fMCPrimTrackHist->SetBinEdges(3,fBinsMult);
 455   fMCPrimTrackHist->GetAxis(0)->SetTitle("Zv (cm)");
 456   fMCPrimTrackHist->GetAxis(1)->SetTitle("MC p_{T} (GeV/c)");
 457   fMCPrimTrackHist->GetAxis(2)->SetTitle("#eta (MC)");
 458   fMCPrimTrackHist->GetAxis(3)->SetTitle("true multiplicity (MC)");
 459   fMCPrimTrackHist->Sumw2();
 460
 461   //
 462   // rec. vs MC correlation matrices
 463   //
 464   Int_t binsMultTrueEventMatrix[3]={fMultNbins,fMultNbins,fMultNbins};
 465 //   Double_t minMultTrueEventMatrix[3]={-0.5,-0.5,-0.5};
 466 //   Double_t maxMultTrueEventMatrix[3]={149.5,149.5,149.5};
 467   fEventMultCorrelationMatrix = new THnSparseF("fEventMultCorrelationMatrix","multRecMult:multTrueMC:multMB",3,binsMultTrueEventMatrix); //,minMultTrueEventMatrix,maxMultTrueEventMatrix);
 468   fEventMultCorrelationMatrix->SetBinEdges(0,fBinsMult);
 469   fEventMultCorrelationMatrix->SetBinEdges(1,fBinsMult);
 470   fEventMultCorrelationMatrix->SetBinEdges(2,fBinsMult);
 471   fEventMultCorrelationMatrix->GetAxis(0)->SetTitle("multiplicity (multCuts)");
 472   fEventMultCorrelationMatrix->GetAxis(1)->SetTitle("true multiplicity (MC)");
 473   fEventMultCorrelationMatrix->GetAxis(2)->SetTitle("MB multiplicity");
 474   fEventMultCorrelationMatrix->Sumw2();
 475
 476   Int_t binsTrackPtCorrelationMatrix[3]={fPtCorrNbins,fPtCorrNbins,fEtaNbins};
 477   fTrackPtCorrelationMatrix = new THnSparseF("fTrackPtCorrelationMatrix","Pt:mcPt:mcEta",3,binsTrackPtCorrelationMatrix);
 478   fTrackPtCorrelationMatrix->SetBinEdges(0,fBinsPtCorr);
 479   fTrackPtCorrelationMatrix->SetBinEdges(1,fBinsPtCorr);
 480   fTrackPtCorrelationMatrix->SetBinEdges(2,fBinsEta);
 481   fTrackPtCorrelationMatrix->GetAxis(0)->SetTitle("Pt (GeV/c)");
 482   fTrackPtCorrelationMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 483   fTrackPtCorrelationMatrix->GetAxis(2)->SetTitle("mcEta");
 484   fTrackPtCorrelationMatrix->Sumw2();
 485
 486   //
 487   // Efficiency and contamination correction matrices
 488   //
 489   Int_t binsEventMatrix[2]={fZvNbins,fMultNbins};
 490 //   Double_t minEventMatrix[2]={-25.,-0.5};
 491 //   Double_t maxEventMatrix[2]={25.,149.5};
 492
 493   fGenEventMatrix = new THnSparseF("fGenEventMatrix","mcZv:multTrueMC",2,binsEventMatrix); //,minEventMatrix,maxEventMatrix);
 494   fGenEventMatrix->SetBinEdges(0,fBinsZv);
 495   fGenEventMatrix->SetBinEdges(1,fBinsMult);
 496   fGenEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 497   fGenEventMatrix->GetAxis(1)->SetTitle("true multiplicity (MC)");
 498   fGenEventMatrix->Sumw2();
 499
 500   fGenSDEventMatrix = new THnSparseF("fGenSDEventMatrix","mcZv:multTrueMC",2,binsEventMatrix); //,minEventMatrix,maxEventMatrix);
 501   fGenSDEventMatrix->SetBinEdges(0,fBinsZv);
 502   fGenSDEventMatrix->SetBinEdges(1,fBinsMult);
 503   fGenSDEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 504   fGenSDEventMatrix->GetAxis(1)->SetTitle("true multiplicity (MC)");
 505   fGenSDEventMatrix->Sumw2();
 506
 507   fGenDDEventMatrix = new THnSparseF("fGenDDEventMatrix","mcZv:multTrueMC",2,binsEventMatrix); //,minEventMatrix,maxEventMatrix);
 508   fGenDDEventMatrix->SetBinEdges(0,fBinsZv);
 509   fGenDDEventMatrix->SetBinEdges(1,fBinsMult);
 510   fGenDDEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 511   fGenDDEventMatrix->GetAxis(1)->SetTitle("true multiplicity (MC)");
 512   fGenDDEventMatrix->Sumw2();
 513
 514   fGenNDEventMatrix = new THnSparseF("fGenNDEventMatrix","mcZv:multTrueMC",2,binsEventMatrix); //,minEventMatrix,maxEventMatrix);
 515   fGenNDEventMatrix->SetBinEdges(0,fBinsZv);
 516   fGenNDEventMatrix->SetBinEdges(1,fBinsMult);
 517   fGenNDEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 518   fGenNDEventMatrix->GetAxis(1)->SetTitle("true multiplicity (MC)");
 519   fGenNDEventMatrix->Sumw2();
 520
 521   fGenNSDEventMatrix = new THnSparseF("fGenNSDEventMatrix","mcZv:multTrueMC",2,binsEventMatrix); //,minEventMatrix,maxEventMatrix);
 522   fGenNSDEventMatrix->SetBinEdges(0,fBinsZv);
 523   fGenNSDEventMatrix->SetBinEdges(1,fBinsMult);
 524   fGenNSDEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 525   fGenNSDEventMatrix->GetAxis(1)->SetTitle("true multiplicity (MC)");
 526   fGenNSDEventMatrix->Sumw2();
 527
 528   //
 529   fTriggerEventMatrix = new THnSparseF("fTriggerEventMatrix","mcZv:multTrueMC",2,binsEventMatrix); //,minEventMatrix,maxEventMatrix);
 530   fTriggerEventMatrix->SetBinEdges(0,fBinsZv);
 531   fTriggerEventMatrix->SetBinEdges(1,fBinsMult);
 532   fTriggerEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 533   fTriggerEventMatrix->GetAxis(1)->SetTitle("true multiplicity (MC)");
 534   fTriggerEventMatrix->Sumw2();
 535
 536   fTriggerSDEventMatrix = new THnSparseF("fTriggerSDEventMatrix","mcZv:multTrueMC",2,binsEventMatrix); //,minEventMatrix,maxEventMatrix);
 537   fTriggerSDEventMatrix->SetBinEdges(0,fBinsZv);
 538   fTriggerSDEventMatrix->SetBinEdges(1,fBinsMult);
 539   fTriggerSDEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 540   fTriggerSDEventMatrix->GetAxis(1)->SetTitle("true multiplicity (MC)");
 541   fTriggerSDEventMatrix->Sumw2();
 542
 543   fTriggerDDEventMatrix = new THnSparseF("fTriggerDDEventMatrix","mcZv:multTrueMC",2,binsEventMatrix); //,minEventMatrix,maxEventMatrix);
 544   fTriggerDDEventMatrix->SetBinEdges(0,fBinsZv);
 545   fTriggerDDEventMatrix->SetBinEdges(1,fBinsMult);
 546   fTriggerDDEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 547   fTriggerDDEventMatrix->GetAxis(1)->SetTitle("true multiplicity (MC)");
 548   fTriggerDDEventMatrix->Sumw2();
 549
 550   fTriggerNDEventMatrix = new THnSparseF("fTriggerNDEventMatrix","mcZv:multTrueMC",2,binsEventMatrix); //,minEventMatrix,maxEventMatrix);
 551   fTriggerNDEventMatrix->SetBinEdges(0,fBinsZv);
 552   fTriggerNDEventMatrix->SetBinEdges(1,fBinsMult);
 553   fTriggerNDEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 554   fTriggerNDEventMatrix->GetAxis(1)->SetTitle("true multiplicity (MC)");
 555   fTriggerNDEventMatrix->Sumw2();
 556
 557   fTriggerNSDEventMatrix = new THnSparseF("fTriggerNSDEventMatrix","mcZv:multTrueMC",2,binsEventMatrix); //,minEventMatrix,maxEventMatrix);
 558   fTriggerNSDEventMatrix->SetBinEdges(0,fBinsZv);
 559   fTriggerNSDEventMatrix->SetBinEdges(1,fBinsMult);
 560   fTriggerNSDEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 561   fTriggerNSDEventMatrix->GetAxis(1)->SetTitle("true multiplicity (MC)");
 562   fTriggerNSDEventMatrix->Sumw2();
 563
 564   //
 565   fRecEventMatrix = new THnSparseF("fRecEventMatrix","mcZv:multTrueMC",2,binsEventMatrix); //,minEventMatrix,maxEventMatrix);
 566   fRecEventMatrix->SetBinEdges(0,fBinsZv);
 567   fRecEventMatrix->SetBinEdges(1,fBinsMult);
 568   fRecEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 569   fRecEventMatrix->GetAxis(1)->SetTitle("true multiplicity (MC)");
 570   fRecEventMatrix->Sumw2();
 571
 572   fRecSDEventMatrix = new THnSparseF("fRecSDEventMatrix","mcZv:multTrueMC",2,binsEventMatrix); //,minEventMatrix,maxEventMatrix);
 573   fRecSDEventMatrix->SetBinEdges(0,fBinsZv);
 574   fRecSDEventMatrix->SetBinEdges(1,fBinsMult);
 575   fRecSDEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 576   fRecSDEventMatrix->GetAxis(1)->SetTitle("true multiplicity (MC)");
 577   fRecSDEventMatrix->Sumw2();
 578
 579   fRecDDEventMatrix = new THnSparseF("fRecDDEventMatrix","mcZv:multTrueMC",2,binsEventMatrix); //,minEventMatrix,maxEventMatrix);
 580   fRecDDEventMatrix->SetBinEdges(0,fBinsZv);
 581   fRecDDEventMatrix->SetBinEdges(1,fBinsMult);
 582   fRecDDEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 583   fRecDDEventMatrix->GetAxis(1)->SetTitle("true multiplicity (MC)");
 584   fRecDDEventMatrix->Sumw2();
 585
 586   fRecNDEventMatrix = new THnSparseF("fRecNDEventMatrix","mcZv:multTrueMC",2,binsEventMatrix); //,minEventMatrix,maxEventMatrix);
 587   fRecNDEventMatrix->SetBinEdges(0,fBinsZv);
 588   fRecNDEventMatrix->SetBinEdges(1,fBinsMult);
 589   fRecNDEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 590   fRecNDEventMatrix->GetAxis(1)->SetTitle("true multiplicity (MC)");
 591   fRecNDEventMatrix->Sumw2();
 592
 593   fRecNSDEventMatrix = new THnSparseF("fRecNSDEventMatrix","mcZv:multTrueMC",2,binsEventMatrix); //,minEventMatrix,maxEventMatrix);
 594   fRecNSDEventMatrix->SetBinEdges(0,fBinsZv);
 595   fRecNSDEventMatrix->SetBinEdges(1,fBinsMult);
 596   fRecNSDEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 597   fRecNSDEventMatrix->GetAxis(1)->SetTitle("true multiplicity (MC)");
 598   fRecNSDEventMatrix->Sumw2();
 599
 600   fRecCandleEventMatrix = new THnSparseF("fRecCandleEventMatrix","mcZv:multTrueMC",2,binsEventMatrix); //,minEventMatrix,maxEventMatrix);
 601   fRecCandleEventMatrix->SetBinEdges(0,fBinsZv);
 602   fRecCandleEventMatrix->SetBinEdges(1,fBinsMult);
 603   fRecCandleEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 604   fRecCandleEventMatrix->GetAxis(1)->SetTitle("true multiplicity (MC)");
 605   fRecCandleEventMatrix->Sumw2();
 606
 607   //
 608   // track to event corrections
 609   //
 610
 611   fGenTrackEventMatrix = new THnSparseF("fGenTrackEventMatrix","mcZv:mcPt:mcEta:multTrueMC",4,binsTrackEventCorrMatrix);
 612   fGenTrackEventMatrix->SetBinEdges(0,fBinsZv);
 613   fGenTrackEventMatrix->SetBinEdges(1,fBinsPtCorr);
 614   fGenTrackEventMatrix->SetBinEdges(2,fBinsEta);
 615   fGenTrackEventMatrix->SetBinEdges(3,fBinsMult);
 616   fGenTrackEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 617   fGenTrackEventMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 618   fGenTrackEventMatrix->GetAxis(2)->SetTitle("mcEta");
 619   fGenTrackEventMatrix->GetAxis(3)->SetTitle("true multiplicity (MC)");
 620   fGenTrackEventMatrix->Sumw2();
 621
 622   fGenTrackSDEventMatrix = new THnSparseF("fGenTrackSDEventMatrix","mcZv:mcPt:mcEta:multTrueMC",4,binsTrackEventCorrMatrix);
 623   fGenTrackSDEventMatrix->SetBinEdges(0,fBinsZv);
 624   fGenTrackSDEventMatrix->SetBinEdges(1,fBinsPtCorr);
 625   fGenTrackSDEventMatrix->SetBinEdges(2,fBinsEta);
 626   fGenTrackSDEventMatrix->SetBinEdges(3,fBinsMult);
 627   fGenTrackSDEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 628   fGenTrackSDEventMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 629   fGenTrackSDEventMatrix->GetAxis(2)->SetTitle("mcEta");
 630   fGenTrackSDEventMatrix->GetAxis(3)->SetTitle("true multiplicity (MC)");
 631   fGenTrackSDEventMatrix->Sumw2();
 632
 633   fGenTrackDDEventMatrix = new THnSparseF("fGenTrackDDEventMatrix","mcZv:mcPt:mcEta:multTrueMC",4,binsTrackEventCorrMatrix);
 634   fGenTrackDDEventMatrix->SetBinEdges(0,fBinsZv);
 635   fGenTrackDDEventMatrix->SetBinEdges(1,fBinsPtCorr);
 636   fGenTrackDDEventMatrix->SetBinEdges(2,fBinsEta);
 637   fGenTrackDDEventMatrix->SetBinEdges(3,fBinsMult);
 638   fGenTrackDDEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 639   fGenTrackDDEventMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 640   fGenTrackDDEventMatrix->GetAxis(2)->SetTitle("mcEta");
 641   fGenTrackDDEventMatrix->GetAxis(3)->SetTitle("true multiplicity (MC)");
 642   fGenTrackDDEventMatrix->Sumw2();
 643
 644   fGenTrackNDEventMatrix = new THnSparseF("fGenTrackNDEventMatrix","mcZv:mcPt:mcEta:multTrueMC",4,binsTrackEventCorrMatrix);
 645   fGenTrackNDEventMatrix->SetBinEdges(0,fBinsZv);
 646   fGenTrackNDEventMatrix->SetBinEdges(1,fBinsPtCorr);
 647   fGenTrackNDEventMatrix->SetBinEdges(2,fBinsEta);
 648   fGenTrackNDEventMatrix->SetBinEdges(3,fBinsMult);
 649   fGenTrackNDEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 650   fGenTrackNDEventMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 651   fGenTrackNDEventMatrix->GetAxis(2)->SetTitle("mcEta");
 652   fGenTrackNDEventMatrix->GetAxis(3)->SetTitle("true multiplicity (MC)");
 653   fGenTrackNDEventMatrix->Sumw2();
 654
 655   fGenTrackNSDEventMatrix = new THnSparseF("fGenTrackNSDEventMatrix","mcZv:mcPt:mcEta:multTrueMC",4,binsTrackEventCorrMatrix);
 656   fGenTrackNSDEventMatrix->SetBinEdges(0,fBinsZv);
 657   fGenTrackNSDEventMatrix->SetBinEdges(1,fBinsPtCorr);
 658   fGenTrackNSDEventMatrix->SetBinEdges(2,fBinsEta);
 659   fGenTrackNSDEventMatrix->SetBinEdges(3,fBinsMult);
 660   fGenTrackNSDEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 661   fGenTrackNSDEventMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 662   fGenTrackNSDEventMatrix->GetAxis(2)->SetTitle("mcEta");
 663   fGenTrackNSDEventMatrix->GetAxis(3)->SetTitle("true multiplicity (MC)");
 664   fGenTrackNSDEventMatrix->Sumw2();
 665
 666
 667   //
 668   fTriggerTrackEventMatrix = new THnSparseF("fTriggerTrackEventMatrix","mcZv:mcPt:mcEta:multTrueMC",4,binsTrackEventCorrMatrix);
 669   fTriggerTrackEventMatrix->SetBinEdges(0,fBinsZv);
 670   fTriggerTrackEventMatrix->SetBinEdges(1,fBinsPtCorr);
 671   fTriggerTrackEventMatrix->SetBinEdges(2,fBinsEta);
 672   fTriggerTrackEventMatrix->SetBinEdges(3,fBinsMult);
 673   fTriggerTrackEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 674   fTriggerTrackEventMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 675   fTriggerTrackEventMatrix->GetAxis(2)->SetTitle("mcEta");
 676   fTriggerTrackEventMatrix->GetAxis(3)->SetTitle("true multiplicity (MC)");
 677   fTriggerTrackEventMatrix->Sumw2();
 678
 679   fTriggerTrackSDEventMatrix = new THnSparseF("fTriggerTrackSDEventMatrix","mcZv:mcPt:mcEta:multTrueMC",4,binsTrackEventCorrMatrix);
 680   fTriggerTrackSDEventMatrix->SetBinEdges(0,fBinsZv);
 681   fTriggerTrackSDEventMatrix->SetBinEdges(1,fBinsPtCorr);
 682   fTriggerTrackSDEventMatrix->SetBinEdges(2,fBinsEta);
 683   fTriggerTrackSDEventMatrix->SetBinEdges(3,fBinsMult);
 684   fTriggerTrackSDEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 685   fTriggerTrackSDEventMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 686   fTriggerTrackSDEventMatrix->GetAxis(2)->SetTitle("mcEta");
 687   fTriggerTrackSDEventMatrix->GetAxis(3)->SetTitle("true multiplicity (MC)");
 688   fTriggerTrackSDEventMatrix->Sumw2();
 689
 690   fTriggerTrackDDEventMatrix = new THnSparseF("fTriggerTrackDDEventMatrix","mcZv:mcPt:mcEta:multTrueMC",4,binsTrackEventCorrMatrix);
 691   fTriggerTrackDDEventMatrix->SetBinEdges(0,fBinsZv);
 692   fTriggerTrackDDEventMatrix->SetBinEdges(1,fBinsPtCorr);
 693   fTriggerTrackDDEventMatrix->SetBinEdges(2,fBinsEta);
 694   fTriggerTrackDDEventMatrix->SetBinEdges(3,fBinsMult);
 695   fTriggerTrackDDEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 696   fTriggerTrackDDEventMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 697   fTriggerTrackDDEventMatrix->GetAxis(2)->SetTitle("mcEta");
 698   fTriggerTrackDDEventMatrix->GetAxis(3)->SetTitle("true multiplicity (MC)");
 699   fTriggerTrackDDEventMatrix->Sumw2();
 700
 701   fTriggerTrackNDEventMatrix = new THnSparseF("fTriggerTrackNDEventMatrix","mcZv:mcPt:mcEta:multTrueMC",4,binsTrackEventCorrMatrix);
 702   fTriggerTrackNDEventMatrix->SetBinEdges(0,fBinsZv);
 703   fTriggerTrackNDEventMatrix->SetBinEdges(1,fBinsPtCorr);
 704   fTriggerTrackNDEventMatrix->SetBinEdges(2,fBinsEta);
 705   fTriggerTrackNDEventMatrix->SetBinEdges(3,fBinsMult);
 706   fTriggerTrackNDEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 707   fTriggerTrackNDEventMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 708   fTriggerTrackNDEventMatrix->GetAxis(2)->SetTitle("mcEta");
 709   fTriggerTrackNDEventMatrix->GetAxis(3)->SetTitle("true multiplicity (MC)");
 710   fTriggerTrackNDEventMatrix->Sumw2();
 711
 712   fTriggerTrackNSDEventMatrix = new THnSparseF("fTriggerTrackNSDEventMatrix","mcZv:mcPt:mcEta:multTrueMC",4,binsTrackEventCorrMatrix);
 713   fTriggerTrackNSDEventMatrix->SetBinEdges(0,fBinsZv);
 714   fTriggerTrackNSDEventMatrix->SetBinEdges(1,fBinsPtCorr);
 715   fTriggerTrackNSDEventMatrix->SetBinEdges(2,fBinsEta);
 716   fTriggerTrackNSDEventMatrix->SetBinEdges(3,fBinsMult);
 717   fTriggerTrackNSDEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 718   fTriggerTrackNSDEventMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 719   fTriggerTrackNSDEventMatrix->GetAxis(2)->SetTitle("mcEta");
 720   fTriggerTrackNSDEventMatrix->GetAxis(3)->SetTitle("true multiplicity (MC)");
 721   fTriggerTrackNSDEventMatrix->Sumw2();
 722
 723   //
 724   fRecTrackEventMatrix = new THnSparseF("fRecTrackEventMatrix","mcZv:mcPt:mcEta:multTrueMC",4,binsTrackEventCorrMatrix);
 725   fRecTrackEventMatrix->SetBinEdges(0,fBinsZv);
 726   fRecTrackEventMatrix->SetBinEdges(1,fBinsPtCorr);
 727   fRecTrackEventMatrix->SetBinEdges(2,fBinsEta);
 728   fRecTrackEventMatrix->SetBinEdges(3,fBinsMult);
 729   fRecTrackEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 730   fRecTrackEventMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 731   fRecTrackEventMatrix->GetAxis(2)->SetTitle("mcEta");
 732   fRecTrackEventMatrix->GetAxis(3)->SetTitle("true multiplicity (MC)");
 733   fRecTrackEventMatrix->Sumw2();
 734
 735   fRecTrackSDEventMatrix = new THnSparseF("fRecTrackSDEventMatrix","mcZv:mcPt:mcEta:multTrueMC",4,binsTrackEventCorrMatrix);
 736   fRecTrackSDEventMatrix->SetBinEdges(0,fBinsZv);
 737   fRecTrackSDEventMatrix->SetBinEdges(1,fBinsPtCorr);
 738   fRecTrackSDEventMatrix->SetBinEdges(2,fBinsEta);
 739   fRecTrackSDEventMatrix->SetBinEdges(3,fBinsMult);
 740   fRecTrackSDEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 741   fRecTrackSDEventMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 742   fRecTrackSDEventMatrix->GetAxis(2)->SetTitle("mcEta");
 743   fRecTrackSDEventMatrix->GetAxis(3)->SetTitle("true multiplicity (MC)");
 744   fRecTrackSDEventMatrix->Sumw2();
 745
 746   fRecTrackDDEventMatrix = new THnSparseF("fRecTrackDDEventMatrix","mcZv:mcPt:mcEta:multTrueMC",4,binsTrackEventCorrMatrix);
 747   fRecTrackDDEventMatrix->SetBinEdges(0,fBinsZv);
 748   fRecTrackDDEventMatrix->SetBinEdges(1,fBinsPtCorr);
 749   fRecTrackDDEventMatrix->SetBinEdges(2,fBinsEta);
 750   fRecTrackDDEventMatrix->SetBinEdges(3,fBinsMult);
 751   fRecTrackDDEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 752   fRecTrackDDEventMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 753   fRecTrackDDEventMatrix->GetAxis(2)->SetTitle("mcEta");
 754   fRecTrackDDEventMatrix->GetAxis(3)->SetTitle("true multiplicity (MC)");
 755   fRecTrackDDEventMatrix->Sumw2();
 756
 757   fRecTrackNDEventMatrix = new THnSparseF("fRecTrackNDEventMatrix","mcZv:mcPt:mcEta:multTrueMC",4,binsTrackEventCorrMatrix);
 758   fRecTrackNDEventMatrix->SetBinEdges(0,fBinsZv);
 759   fRecTrackNDEventMatrix->SetBinEdges(1,fBinsPtCorr);
 760   fRecTrackNDEventMatrix->SetBinEdges(2,fBinsEta);
 761   fRecTrackNDEventMatrix->SetBinEdges(3,fBinsMult);
 762   fRecTrackNDEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 763   fRecTrackNDEventMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 764   fRecTrackNDEventMatrix->GetAxis(2)->SetTitle("mcEta");
 765   fRecTrackNDEventMatrix->GetAxis(3)->SetTitle("true multiplicity (MC)");
 766   fRecTrackNDEventMatrix->Sumw2();
 767
 768   fRecTrackNSDEventMatrix = new THnSparseF("fRecTrackNSDEventMatrix","mcZv:mcPt:mcEta:multTrueMC",4,binsTrackEventCorrMatrix);
 769   fRecTrackNSDEventMatrix->SetBinEdges(0,fBinsZv);
 770   fRecTrackNSDEventMatrix->SetBinEdges(1,fBinsPtCorr);
 771   fRecTrackNSDEventMatrix->SetBinEdges(2,fBinsEta);
 772   fRecTrackNSDEventMatrix->SetBinEdges(3,fBinsMult);
 773   fRecTrackNSDEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 774   fRecTrackNSDEventMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 775   fRecTrackNSDEventMatrix->GetAxis(2)->SetTitle("mcEta");
 776   fRecTrackNSDEventMatrix->GetAxis(3)->SetTitle("true multiplicity (MC)");
 777   fRecTrackNSDEventMatrix->Sumw2();
 778
 779   //
 780   // tracks correction matrices
 781   //
 782   //fGenTrackMatrix = new THnSparseF("fGenTrackMatrix","mcZv:mcPt:mcEta",3,binsTrackMatrix);
 783   fGenTrackMatrix = new THnSparseF("fGenTrackMatrix","mcZv:mcPt:mcEta",3,binsTrackEventCorrMatrix);
 784   fGenTrackMatrix->SetBinEdges(0,fBinsZv);
 785   //fGenTrackMatrix->SetBinEdges(1,binsPt);
 786   fGenTrackMatrix->SetBinEdges(1,fBinsPtCorr);
 787   fGenTrackMatrix->SetBinEdges(2,fBinsEta);
 788   fGenTrackMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 789   fGenTrackMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 790   fGenTrackMatrix->GetAxis(2)->SetTitle("mcEta");
 791   fGenTrackMatrix->Sumw2();
 792
 793   //fGenPrimTrackMatrix = new THnSparseF("fGenPrimTrackMatrix","mcZv:mcPt:mcEta",3,binsTrackMatrix);
 794   fGenPrimTrackMatrix = new THnSparseF("fGenPrimTrackMatrix","mcZv:mcPt:mcEta",3,binsTrackEventCorrMatrix);
 795   fGenPrimTrackMatrix->SetBinEdges(0,fBinsZv);
 796   //fGenPrimTrackMatrix->SetBinEdges(1,binsPt);
 797   fGenPrimTrackMatrix->SetBinEdges(1,fBinsPtCorr);
 798   fGenPrimTrackMatrix->SetBinEdges(2,fBinsEta);
 799   fGenPrimTrackMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 800   fGenPrimTrackMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 801   fGenPrimTrackMatrix->GetAxis(2)->SetTitle("mcEta");
 802   fGenPrimTrackMatrix->Sumw2();
 803
 804   //fRecPrimTrackMatrix = new THnSparseF("fRecPrimTrackMatrix","mcZv:mcPt:mcEta",3,binsTrackMatrix);
 805   fRecPrimTrackMatrix = new THnSparseF("fRecPrimTrackMatrix","mcZv:mcPt:mcEta",3,binsTrackEventCorrMatrix);
 806   fRecPrimTrackMatrix->SetBinEdges(0,fBinsZv);
 807   //fRecPrimTrackMatrix->SetBinEdges(1,binsPt);
 808   fRecPrimTrackMatrix->SetBinEdges(1,fBinsPtCorr);
 809   fRecPrimTrackMatrix->SetBinEdges(2,fBinsEta);
 810   fRecPrimTrackMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 811   fRecPrimTrackMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 812   fRecPrimTrackMatrix->GetAxis(2)->SetTitle("mcEta");
 813   fRecPrimTrackMatrix->Sumw2();
 814
 815   //
 816   //fRecTrackMatrix = new THnSparseF("fRecTrackMatrix","mcZv:mcPt:mcEta",3,binsTrackMatrix);
 817   fRecTrackMatrix = new THnSparseF("fRecTrackMatrix","mcZv:mcPt:mcEta",3,binsTrackEventCorrMatrix);
 818   fRecTrackMatrix->SetBinEdges(0,fBinsZv);
 819   //fRecTrackMatrix->SetBinEdges(1,binsPt);
 820   fRecTrackMatrix->SetBinEdges(1,fBinsPtCorr);
 821   fRecTrackMatrix->SetBinEdges(2,fBinsEta);
 822   fRecTrackMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 823   fRecTrackMatrix->GetAxis(1)->SetTitle("Pt (GeV/c)");
 824   fRecTrackMatrix->GetAxis(2)->SetTitle("Eta");
 825   fRecTrackMatrix->Sumw2();
 826
 827   //fRecSecTrackMatrix = new THnSparseF("fRecSecTrackMatrix","mcZv:mcPt:mcEta",3,binsTrackMatrix);
 828   fRecSecTrackMatrix = new THnSparseF("fRecSecTrackMatrix","mcZv:mcPt:mcEta",3,binsTrackEventCorrMatrix);
 829   fRecSecTrackMatrix->SetBinEdges(0,fBinsZv);
 830   //fRecSecTrackMatrix->SetBinEdges(1,binsPt);
 831   fRecSecTrackMatrix->SetBinEdges(1,fBinsPtCorr);
 832   fRecSecTrackMatrix->SetBinEdges(2,fBinsEta);
 833   fRecSecTrackMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 834   fRecSecTrackMatrix->GetAxis(1)->SetTitle("Pt (GeV/c)");
 835   fRecSecTrackMatrix->GetAxis(2)->SetTitle("Eta");
 836   fRecSecTrackMatrix->Sumw2();
 837
 838   //
 839   //fRecMultTrackMatrix = new THnSparseF("fRecMultTrackMatrix","mcZv:mcPt:mcEta",3,binsTrackMatrix);
 840   fRecMultTrackMatrix = new THnSparseF("fRecMultTrackMatrix","mcZv:mcPt:mcEta",3,binsTrackEventCorrMatrix);
 841   fRecMultTrackMatrix->SetBinEdges(0,fBinsZv);
 842   //fRecMultTrackMatrix->SetBinEdges(1,binsPt);
 843   fRecMultTrackMatrix->SetBinEdges(1,fBinsPtCorr);
 844   fRecMultTrackMatrix->SetBinEdges(2,fBinsEta);
 845   fRecMultTrackMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
 846   fRecMultTrackMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
 847   fRecMultTrackMatrix->GetAxis(2)->SetTitle("mcEta");
 848   fRecMultTrackMatrix->Sumw2();
 849
 850   //
 851   // Control analysis histograms
 852   //
 853
 854   Int_t binsMCEventHist1[3]={100,100,140};
 855   Double_t minMCEventHist1[3]={-0.1,-0.1,-35.};
 856   Double_t maxMCEventHist1[3]={0.1,0.1,35.};
 857   fMCEventHist1 = new THnSparseF("fMCEventHist1","mcXv:mcYv:mcZv",3,binsMCEventHist1,minMCEventHist1,maxMCEventHist1);
 858   fMCEventHist1->GetAxis(0)->SetTitle("mcXv (cm)");
 859   fMCEventHist1->GetAxis(1)->SetTitle("mcYv (cm)");
 860   fMCEventHist1->GetAxis(2)->SetTitle("mcZv (cm)");
 861   fMCEventHist1->Sumw2();
 862
 863   //
 864   Int_t binsRecEventHist1[3]={100,100,140};
 865   Double_t minRecEventHist1[3]={-3.,-3.,-35.};
 866   Double_t maxRecEventHist1[3]={3.,3.,35.};
 867
 868   fRecEventHist1 = new THnSparseF("fRecEventHist1","Xv:Yv:Zv",3,binsRecEventHist1,minRecEventHist1,maxRecEventHist1);
 869   fRecEventHist1->GetAxis(0)->SetTitle("Xv (cm)");
 870   fRecEventHist1->GetAxis(1)->SetTitle("Yv (cm)");
 871   fRecEventHist1->GetAxis(2)->SetTitle("Zv (cm)");
 872   fRecEventHist1->Sumw2();
 873
 874   //
 875   Int_t binsRecEventHist2[3]={fZvNbins,fMultNbins,fMultNbins};
 876 //   Double_t minRecEventHist2[3]={-25.,-0.5,-0.5};
 877 //   Double_t maxRecEventHist2[3]={25.,149.5,149.5};
 878
 879   fRecEventHist2 = new THnSparseF("fRecEventHist2","Zv:multMB:multRecMult",3,binsRecEventHist2); //,minRecEventHist2,maxRecEventHist2);
 880   fRecEventHist2->SetBinEdges(0,fBinsZv);
 881   fRecEventHist2->SetBinEdges(1,fBinsMult);
 882   fRecEventHist2->SetBinEdges(2,fBinsMult);
 883   fRecEventHist2->GetAxis(0)->SetTitle("Zv (cm)");
 884   fRecEventHist2->GetAxis(1)->SetTitle("multiplicity MB");
 885   fRecEventHist2->GetAxis(2)->SetTitle("multiplicity (multCuts)");
 886   fRecEventHist2->Sumw2();
 887
 888   //
 889   Double_t kFact = 0.1;
 890   Int_t binsRecMCEventHist1[3]={100,100,100};
 891   Double_t minRecMCEventHist1[3]={-10.0*kFact,-10.0*kFact,-10.0*kFact};
 892   Double_t maxRecMCEventHist1[3]={10.0*kFact,10.0*kFact,10.0*kFact};
 893
 894   fRecMCEventHist1 = new THnSparseF("fRecMCEventHist1","Xv-mcXv:Yv-mcYv:Zv-mcZv",3,binsRecMCEventHist1,minRecMCEventHist1,maxRecMCEventHist1);
 895   fRecMCEventHist1->GetAxis(0)->SetTitle("Xv-mcXv (cm)");
 896   fRecMCEventHist1->GetAxis(1)->SetTitle("Yv-mcYv (cm)");
 897   fRecMCEventHist1->GetAxis(2)->SetTitle("Zv-mcZv (cm)");
 898   fRecMCEventHist1->Sumw2();
 899
 900   //
 901   Int_t binsRecMCEventHist2[3]={100,100,fMultNbins};
 902   Double_t minRecMCEventHist2[3]={-10.0*kFact,-10.0*kFact,0.0};
 903   Double_t maxRecMCEventHist2[3]={10.0*kFact,10.0*kFact,149.50};
 904
 905   fRecMCEventHist2 = new THnSparseF("fRecMCEventHist2","Xv-mcXv:Zv-mcZv:multMB",3,binsRecMCEventHist2,minRecMCEventHist2,maxRecMCEventHist2);
 906   fRecMCEventHist2->SetBinEdges(2,fBinsMult);
 907   fRecMCEventHist2->GetAxis(0)->SetTitle("Xv-mcXv (cm)");
 908   fRecMCEventHist2->GetAxis(1)->SetTitle("Zv-mcZv (cm)");
 909   fRecMCEventHist2->GetAxis(2)->SetTitle("multiplicity MB");
 910   fRecMCEventHist2->Sumw2();
 911
 912   Int_t binsRecMCEventHist3[2]={fMultNbins,5};
 913   Double_t minRecMCEventHist3[2]={-0.5,0.0};
 914   Double_t maxRecMCEventHist3[2]={149.50,5.0};
 915   fRecMCEventHist3 = new THnSparseF("fRecMCEventHist3","multRecMult:EventType (ND, DD, SD)",2,binsRecMCEventHist3,minRecMCEventHist3,maxRecMCEventHist3);
 916   fRecMCEventHist3->SetBinEdges(0,fBinsMult);
 917   fRecMCEventHist3->GetAxis(0)->SetTitle("multiplicity (multCuts)");
 918   fRecMCEventHist3->GetAxis(1)->SetTitle("EventType");
 919   fRecMCEventHist3->Sumw2();
 920
 921   //
 922   char name[256];
 923   char title[256];
 924   for(Int_t i=0; i<AlidNdPtHelper::kCutSteps; i++)
 925   {
 926   // THnSparse track histograms
 927
 928   Int_t binsMCTrackHist1[3]=  {fPtCorrNbins, fEtaNbins, 90};
 929   Double_t minMCTrackHist1[3]={0.,-1.,0.};
 930   Double_t maxMCTrackHist1[3]={10.,1.,2.*TMath::Pi()};
 931   snprintf(name,256,"fMCTrackHist1_%d",i);
 932   snprintf(title,256,"mcPt:mcEta:mcPhi");
 933
 934   fMCTrackHist1[i] = new THnSparseF(name,title,3,binsMCTrackHist1,minMCTrackHist1,maxMCTrackHist1);
 935   fMCTrackHist1[i]->SetBinEdges(0,fBinsPtCorr);
 936   fMCTrackHist1[i]->SetBinEdges(1,fBinsEta);
 937   fMCTrackHist1[i]->GetAxis(0)->SetTitle("mcPt (GeV/c)");
 938   fMCTrackHist1[i]->GetAxis(1)->SetTitle("mcEta");
 939   fMCTrackHist1[i]->GetAxis(2)->SetTitle("mcPhi (rad)");
 940   fMCTrackHist1[i]->Sumw2();
 941
 942   Int_t binsMCPrimTrackHist1[5]=  {fPtCorrNbins,fEtaNbins,6,20,4000};
 943   Double_t minMCPrimTrackHist1[5]={0.,-1.,0.,0.,0.};
 944   Double_t maxMCPrimTrackHist1[5]={10.,1.,6.,20.,4000.};
 945   snprintf(name,256,"fMCPrimTrackHist1_%d",i);
 946   snprintf(title,256,"mcPt:mcEta:pid:mech:mother");
 947
 948   fMCPrimTrackHist1[i] = new THnSparseF(name,title,5,binsMCPrimTrackHist1,minMCPrimTrackHist1,maxMCPrimTrackHist1);
 949   fMCPrimTrackHist1[i]->SetBinEdges(0,fBinsPtCorr);
 950   fMCPrimTrackHist1[i]->SetBinEdges(1,fBinsEta);
 951   fMCPrimTrackHist1[i]->GetAxis(0)->SetTitle("mcPt (GeV/c)");
 952   fMCPrimTrackHist1[i]->GetAxis(1)->SetTitle("mcEta");
 953   fMCPrimTrackHist1[i]->GetAxis(2)->SetTitle("pid");
 954   fMCPrimTrackHist1[i]->GetAxis(3)->SetTitle("mech");
 955   fMCPrimTrackHist1[i]->GetAxis(4)->SetTitle("mother");
 956   fMCPrimTrackHist1[i]->Sumw2();
 957
 958   Int_t binsMCPrimTrackHist2[3]=  {4000,20,4000};
 959   Double_t minMCPrimTrackHist2[3]={0.,0.,0.};
 960   Double_t maxMCPrimTrackHist2[3]={4000.,20.,4000.};
 961   snprintf(name,256,"fMCPrimTrackHist2_%d",i);
 962   snprintf(title,256,"pdg:mech:mother");
 963
 964   fMCPrimTrackHist2[i] = new THnSparseF(name,title,3,binsMCPrimTrackHist2,minMCPrimTrackHist2,maxMCPrimTrackHist2);
 965   fMCPrimTrackHist2[i]->GetAxis(0)->SetTitle("pdg");
 966   fMCPrimTrackHist2[i]->GetAxis(1)->SetTitle("mech");
 967   fMCPrimTrackHist2[i]->GetAxis(2)->SetTitle("mother");
 968   fMCPrimTrackHist2[i]->Sumw2();
 969
 970   Int_t binsMCSecTrackHist1[5]=  {fPtCorrNbins,fEtaNbins,6,20,4000};
 971   Double_t minMCSecTrackHist1[5]={0.,-1.,0.,0.,0.};
 972   Double_t maxMCSecTrackHist1[5]={10.,1.,6.,20.,4000.};
 973   snprintf(name,256,"fMCSecTrackHist1_%d",i);
 974   snprintf(title,256,"mcPt:mcEta:pid:mech:mother");
 975
 976   fMCSecTrackHist1[i] = new THnSparseF(name,title,5,binsMCSecTrackHist1,minMCSecTrackHist1,maxMCSecTrackHist1);
 977   fMCSecTrackHist1[i]->SetBinEdges(0,fBinsPtCorr);
 978   fMCSecTrackHist1[i]->SetBinEdges(1,fBinsEta);
 979   fMCSecTrackHist1[i]->GetAxis(0)->SetTitle("mcPt (GeV/c)");
 980   fMCSecTrackHist1[i]->GetAxis(1)->SetTitle("mcEta");
 981   fMCSecTrackHist1[i]->GetAxis(2)->SetTitle("pid");
 982   fMCSecTrackHist1[i]->GetAxis(3)->SetTitle("mech");
 983   fMCSecTrackHist1[i]->GetAxis(4)->SetTitle("mother");
 984   fMCSecTrackHist1[i]->Sumw2();
 985
 986   //
 987
 988   //
 989
 990   Int_t binsRecTrackHist1[3]={fPtNbins,fEtaNbins,90};
 991   Double_t minRecTrackHist1[3]={0.,-1.,0.};
 992   Double_t maxRecTrackHist1[3]={10.,1.,2.*TMath::Pi()};
 993   snprintf(name,256,"fRecTrackHist1_%d",i);
 994   snprintf(title,256,"Pt:Eta:Phi");
 995   fRecTrackHist1[i] = new THnSparseF(name,title,3,binsRecTrackHist1,minRecTrackHist1,maxRecTrackHist1);
 996   fRecTrackHist1[i]->SetBinEdges(0,fBinsPt);
 997   fRecTrackHist1[i]->SetBinEdges(1,fBinsEta);
 998   fRecTrackHist1[i]->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
 999   fRecTrackHist1[i]->GetAxis(1)->SetTitle("#eta");
1000   fRecTrackHist1[i]->GetAxis(2)->SetTitle("#phi (rad)");
1001   fRecTrackHist1[i]->Sumw2();
1002
1003   //
1004   Int_t binsRecTrackMultHist1[2]={fPtNbins,fMultNbins};
1005   snprintf(name,256,"fRecTrackMultHist_%d",i);
1006   snprintf(title,256,"Pt:Mult");
1007   fRecTrackMultHist1[i] = new THnSparseF(name,title,2,binsRecTrackMultHist1); //,minRecTrackMultHist1,maxRecTrackMultHist1);
1008   fRecTrackMultHist1[i]->SetBinEdges(0,fBinsPt);
1009   fRecTrackMultHist1[i]->SetBinEdges(1,fBinsMult);
1010   fRecTrackMultHist1[i]->GetAxis(0)->SetTitle("Pt (GeV/c)");
1011   fRecTrackMultHist1[i]->GetAxis(1)->SetTitle("multiplicity");
1012   fRecTrackMultHist1[i]->Sumw2();
1013   }
1014
1015   Int_t binsRecMCTrackHist1[4] = {fPtCorrNbins,fEtaNbins,100,100};
1016   Double_t minRecMCTrackHist1[4]={0.,-1.,-0.5,-0.5};
1017   Double_t maxRecMCTrackHist1[4]={20.,1.,0.5,0.5};
1018   snprintf(name,256,"fRecMCTrackHist1");
1019   snprintf(title,256,"mcPt:mcEta:(Pt-mcPt)/mcPt:(Eta-mcEta)");
1020
1021   fRecMCTrackHist1 = new THnSparseF(name,title,4,binsRecMCTrackHist1,minRecMCTrackHist1,maxRecMCTrackHist1);
1022   fRecMCTrackHist1->SetBinEdges(0,fBinsPtCorr);
1023   fRecMCTrackHist1->SetBinEdges(1,fBinsEta);
1024   fRecMCTrackHist1->GetAxis(0)->SetTitle("mcPt (GeV/c)");
1025   fRecMCTrackHist1->GetAxis(1)->SetTitle("mcEta");
1026   fRecMCTrackHist1->GetAxis(2)->SetTitle("(Pt-mcPt)/mcPt");
1027   fRecMCTrackHist1->GetAxis(3)->SetTitle("Eta-mcEta");
1028
1029   Int_t binsMCMultRecTrackHist1[3] = {fPtCorrNbins,fEtaNbins,6};
1030   Double_t minMCMultRecTrackHist1[3]={0.,-1.,0.};
1031   Double_t maxMCMultRecTrackHist1[3]={20.,1.,6.};
1032   snprintf(name,256,"fMCMultRecTrackHist1");
1033   snprintf(title,256,"mcPt:mcEta:pid");
1034   fMCMultRecTrackHist1 = new THnSparseF(name,title,3,binsMCMultRecTrackHist1,minMCMultRecTrackHist1,maxMCMultRecTrackHist1);
1035   fMCMultRecTrackHist1->SetBinEdges(0,fBinsPtCorr);
1036   fMCMultRecTrackHist1->SetBinEdges(1,fBinsEta);
1037   fMCMultRecTrackHist1->GetAxis(0)->SetTitle("mcPt (GeV/c)");
1038   fMCMultRecTrackHist1->GetAxis(1)->SetTitle("mcEta");
1039   fMCMultRecTrackHist1->GetAxis(2)->SetTitle("pid");
1040
1041   //nClust:chi2PerClust:pt:eta:phi
1042   Int_t binsRecTrackHist2[5]={160,100,fPtNbins,fEtaNbins,90};
1043   Double_t minRecTrackHist2[5]={0., 0., 0, -1.5, 0.};
1044   Double_t maxRecRecTrackHist2[5]={160.,10., 16, 1.5, 2.*TMath::Pi()};
1045
1046   fRecTrackHist2 = new THnSparseF("fRecTrackHist2","nClust:chi2PerClust:pt:eta:phi",5,binsRecTrackHist2,minRecTrackHist2,maxRecRecTrackHist2);
1047   fRecTrackHist2->SetBinEdges(2,fBinsPt);
1048   fRecTrackHist2->SetBinEdges(3,fBinsEta);
1049   fRecTrackHist2->GetAxis(0)->SetTitle("nClust");
1050   fRecTrackHist2->GetAxis(1)->SetTitle("chi2PerClust");
1051   fRecTrackHist2->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
1052   fRecTrackHist2->GetAxis(3)->SetTitle("#eta");
1053   fRecTrackHist2->GetAxis(4)->SetTitle("#phi (rad)");
1054   fRecTrackHist2->Sumw2();
1055
1056   Int_t binsEventCount[3]={2,2,2};
1057   Double_t minEventCount[3]={0,0,0};
1058   Double_t maxEventCount[3]={2,2,2};
1059   fEventCount = new THnSparseF("fEventCount","trig vs trig+vertex",3,binsEventCount,minEventCount,maxEventCount);
1060   fEventCount->GetAxis(0)->SetTitle("trig");
1061   fEventCount->GetAxis(1)->SetTitle("trig+vert");
1062   fEventCount->GetAxis(2)->SetTitle("selected");
1063   fEventCount->Sumw2();
1064
1065   // init folder
1066   fAnalysisFolder = CreateFolder("folderdNdPt","Analysis dNdPt Folder");
1067
1068   // set init flag
1069   fIsInit = kTRUE;
1070 }
1071
1072 //_____________________________________________________________________________
1073 void AlidNdPtAnalysis::Process(AliESDEvent *const esdEvent, AliMCEvent *const mcEvent)
1074 {
1075   //  init if not done already
1076   if (!fIsInit) { Init(); }
1077
1078   //
1079   // Process real and/or simulated events
1080   //
1081   if(!esdEvent) {
1082     AliDebug(AliLog::kError, "esdEvent not available");
1083     return;
1084   }
1085
1086   // get selection cuts
1087   AlidNdPtEventCuts *evtCuts = GetEventCuts();
1088   AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts();
1089   AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
1090   AlidNdPtAcceptanceCuts *multAccCuts = GetMultAcceptanceCuts();
1091
1092   if(!evtCuts || !accCuts  || !esdTrackCuts) {
1093     AliDebug(AliLog::kError, "cuts not available");
1094     return;
1095   }
1096
1097   // trigger selection
1098   Bool_t isEventTriggered = kTRUE;
1099   AliPhysicsSelection *physicsSelection = NULL;
1100   AliTriggerAnalysis* triggerAnalysis = NULL;
1101
1102   //
1103   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
1104   if (!inputHandler)
1105   {
1106     Printf("ERROR: Could not receive input handler");
1107     return;
1108   }
1109
1110   if(evtCuts->IsTriggerRequired())
1111   {
1112     // always MB
1113     //isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
1114     isEventTriggered = inputHandler->IsEventSelected() & GetTriggerMask();
1115
1116     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
1117     if(!physicsSelection) return;
1118     //SetPhysicsTriggerSelection(physicsSelection);
1119
1120     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
1121       // set trigger (V0AND)
1122       triggerAnalysis = physicsSelection->GetTriggerAnalysis();
1123       if(!triggerAnalysis) return;
1124       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
1125     }
1126
1127     // calculate LHC background
1128     if(!IsUseMCInfo())
1129     {
1130       //
1131       // 0-multiplicity bin for LHC background correction
1132       //
1133       /* bin0 done in the train
1134       if( GetAnalysisMode() == AlidNdPtHelper::kTPCITS ||
1135           GetAnalysisMode() == AlidNdPtHelper::kTPCTrackSPDvtx ||
1136           GetAnalysisMode() == AlidNdPtHelper::kTPCTrackSPDvtxUpdate ||
1137           GetAnalysisMode() == AlidNdPtHelper::kTPCITSHybridTrackSPDvtx ||
1138           GetAnalysisMode() == AlidNdPtHelper::kTPCITSHybridTrackSPDvtxDCArPt )
1139       {
1140          physicsSelection->SetBin0CallbackViaPointer(&AlidNdPtAnalysis::IsBinZeroTrackSPDvtx);
1141       } else {
1142          physicsSelection->SetBin0CallbackViaPointer(&AlidNdPtAnalysis::IsBinZeroSPDvtx);
1143       }
1144       */
1145     }
1146   }
1147
1148
1149   // use MC information
1150   AliHeader* header = 0;
1151   AliGenEventHeader* genHeader = 0;
1152   AliStack* stack = 0;
1153   TArrayF vtxMC(3);
1154   AliPWG0Helper::MCProcessType evtType = AliPWG0Helper::kInvalidProcess;
1155
1156   Int_t multMCTrueTracks = 0;
1157   if(IsUseMCInfo())
1158   {
1159     //
1160     if(!mcEvent) {
1161       AliDebug(AliLog::kError, "mcEvent not available");
1162       return;
1163     }
1164     // get MC event header
1165     header = mcEvent->Header();
1166     if (!header) {
1167       AliDebug(AliLog::kError, "Header not available");
1168       return;
1169     }
1170     // MC particle stack
1171     stack = mcEvent->Stack();
1172     if (!stack) {
1173       AliDebug(AliLog::kError, "Stack not available");
1174       return;
1175     }
1176     // get event type (ND=0x1, DD=0x2, SD=0x4)
1177     evtType = AliPWG0Helper::GetEventProcessType(header);
1178     AliDebug(AliLog::kDebug+1, Form("Found process type %d", evtType));
1179
1180     // get MC vertex
1181     genHeader = header->GenEventHeader();
1182     if (!genHeader) {
1183       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
1184       return;
1185     }
1186     genHeader->PrimaryVertex(vtxMC);
1187
1188     Double_t vMCEventHist1[3]={vtxMC[0],vtxMC[1],vtxMC[2]};
1189     fMCEventHist1->Fill(vMCEventHist1);
1190
1191     // multipliticy of all MC primary tracks
1192     // in Zv, pt and eta ranges)
1193    multMCTrueTracks = AlidNdPtHelper::GetMCTrueTrackMult(mcEvent,evtCuts,multAccCuts);
1194
1195   } // end bUseMC
1196
1197
1198   // get reconstructed vertex
1199   const AliESDVertex* vtxESD = 0;
1200   Bool_t isRecVertex = kFALSE;
1201   if(evtCuts->IsRecVertexRequired())
1202   {
1203     Bool_t bRedoTPCVertex = evtCuts->IsRedoTPCVertex();
1204     Bool_t bUseConstraints = evtCuts->IsUseBeamSpotConstraint();
1205     vtxESD = AlidNdPtHelper::GetVertex(esdEvent,evtCuts,accCuts,esdTrackCuts,GetAnalysisMode(),kFALSE,bRedoTPCVertex,bUseConstraints);
1206     if(!vtxESD) return;
1207     isRecVertex = AlidNdPtHelper::TestRecVertex(vtxESD, esdEvent->GetPrimaryVertexSPD(), GetAnalysisMode(), kFALSE);
1208   }
1209
1210   if( IsUseMCInfo() && !evtCuts->IsRecVertexRequired() ) {
1211     vtxESD = new AliESDVertex(vtxMC[2],10.,genHeader->NProduced(),"smearMC");
1212     if(!vtxESD) return;
1213     isRecVertex = kTRUE;
1214   }
1215   if(!vtxESD) return;
1216
1217   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD) && isRecVertex;
1218   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
1219   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
1220
1221   Bool_t isTrigAndVertex = isEventTriggered && isEventOK;
1222
1223   Double_t vEventCount[3] = { (isEventTriggered && kTRUE) , isTrigAndVertex,  isTrigAndVertex && (TMath::Abs(vtxESD->GetZv()) < 10.)};
1224   fEventCount->Fill(vEventCount);
1225
1226   // vertex contributors
1227   Int_t multMBTracks = 0;
1228   if(GetAnalysisMode() == AlidNdPtHelper::kTPC)
1229   {
1230      if(vtxESD && vtxESD->GetStatus() && isRecVertex)
1231        multMBTracks = AlidNdPtHelper::GetTPCMBTrackMult(esdEvent,evtCuts,accCuts,esdTrackCuts);
1232   }
1233   else if( GetAnalysisMode() == AlidNdPtHelper::kTPCSPDvtx ||
1234            GetAnalysisMode()==AlidNdPtHelper::kTPCSPDvtxUpdate )
1235   {
1236      const AliMultiplicity* mult = esdEvent->GetMultiplicity();
1237      if(mult && vtxESD && vtxESD->GetStatus() && isRecVertex)
1238        multMBTracks = mult->GetNumberOfTracklets();
1239   }
1240   else if( GetAnalysisMode() == AlidNdPtHelper::kTPCITS ||
1241            GetAnalysisMode()==AlidNdPtHelper::kTPCITSHybrid ||
1242            GetAnalysisMode() == AlidNdPtHelper::kTPCTrackSPDvtx ||
1243            GetAnalysisMode() == AlidNdPtHelper::kTPCTrackSPDvtxUpdate ||
1244            GetAnalysisMode() == AlidNdPtHelper::kTPCITSHybridTrackSPDvtx ||
1245            GetAnalysisMode() == AlidNdPtHelper::kTPCITSHybridTrackSPDvtxDCArPt ||
1246            GetAnalysisMode() == AlidNdPtHelper::kITSStandAloneTrackSPDvtx ||
1247            GetAnalysisMode() == AlidNdPtHelper::kITSStandAloneTPCTrackSPDvtx
1248            )
1249   {
1250      if(vtxESD && vtxESD->GetStatus() && isRecVertex)
1251        multMBTracks = vtxESD->GetNContributors();
1252   }
1253   else {
1254     AliDebug(AliLog::kError, Form("Found analysis type %d", GetAnalysisMode()));
1255     return;
1256   }
1257
1258   Bool_t isEventSelected = kTRUE;
1259   if(evtCuts->IsEventSelectedRequired())
1260   {
1261      // select events with at least
1262      // one prompt track in acceptance
1263      // pT>0.5 GeV/c, |eta|<0.8 for the Cross Section studies
1264
1265      isEventSelected = AlidNdPtHelper::SelectEvent(esdEvent,esdTrackCuts);
1266      //printf("isEventSelected %d \n", isEventSelected);
1267   }
1268
1269   TObjArray *allChargedTracks=0;
1270   //Int_t multAll=0, multAcc=0, multRec=0;
1271   Int_t multAll=0, multRec=0;
1272   Int_t *labelsAll=0, *labelsAcc=0, *labelsRec=0;
1273
1274
1275   // check event cuts
1276   Int_t multRecMult=0;
1277   if(isEventOK && isEventTriggered && isEventSelected)
1278   {
1279     // get all charged tracks
1280     allChargedTracks = AlidNdPtHelper::GetAllChargedTracks(esdEvent,GetAnalysisMode());
1281     if(!allChargedTracks) return;
1282
1283     Int_t entries = allChargedTracks->GetEntries();
1284     //printf("entries %d \n",entries);
1285
1286
1287     // calculate mult of reconstructed tracks
1288     for(Int_t i=0; i<entries;++i)
1289     {
1290       AliESDtrack *track = (AliESDtrack*)allChargedTracks->At(i);
1291       if(!track) continue;
1292       if(track->Charge()==0) continue;
1293
1294
1295       // only postive charged
1296       if(GetParticleMode() == AlidNdPtHelper::kPlus && track->Charge() < 0)
1297         continue;
1298
1299       // only negative charged
1300       if(GetParticleMode() == AlidNdPtHelper::kMinus && track->Charge() > 0)
1301         continue;
1302
1303       if(GetMultTrackCuts()->AcceptTrack(track))
1304       {
1305           if(GetMultAcceptanceCuts()->AcceptTrack(track)) multRecMult++;
1306       }
1307     }
1308
1309
1310     labelsAll = new Int_t[entries];
1311     labelsAcc = new Int_t[entries];
1312     labelsRec = new Int_t[entries];
1313     for(Int_t i=0; i<entries;++i)
1314     {
1315       AliESDtrack *track = (AliESDtrack*)allChargedTracks->At(i);
1316       if(!track) continue;
1317       if(track->Charge()==0) continue;
1318
1319       // only postive charged
1320       if(GetParticleMode() == AlidNdPtHelper::kPlus && track->Charge() < 0)
1321         continue;
1322
1323       // only negative charged
1324       if(GetParticleMode() == AlidNdPtHelper::kMinus && track->Charge() > 0)
1325         continue;
1326
1327       FillHistograms(track,stack,vtxESD->GetZv(),AlidNdPtHelper::kAllTracks, multRecMult);
1328       labelsAll[multAll] = TMath::Abs(track->GetLabel());
1329       multAll++;
1330
1331       //if(accCuts->AcceptTrack(track)) {
1332       //FillHistograms(track,stack,AlidNdPtHelper::kAccTracks);
1333       //labelsAcc[multAcc] = TMath::Abs(track->GetLabel());
1334       //multAcc++;
1335       //}
1336
1337       // esd track selection
1338       // ITS stand alone
1339       if(GetAnalysisMode() == AlidNdPtHelper::kITSStandAloneTrackSPDvtx)
1340       {
1341         if(!(track->GetStatus() & AliESDtrack::kITSpureSA)) continue;
1342         if(!(track->GetStatus() & AliESDtrack::kITSrefit)) continue;
1343         if(track->GetNcls(0)<4) continue;
1344         if(!track->HasPointOnITSLayer(0) && !track->HasPointOnITSLayer(1)) continue;
1345       }
1346       else if(GetAnalysisMode() == AlidNdPtHelper::kITSStandAloneTPCTrackSPDvtx)
1347       {
1348         //
1349         // ITS and TPC stand alone tracks
1350         //
1351         if(!(track->GetStatus() & AliESDtrack::kITSpureSA)) continue;
1352         if(!(track->GetStatus() & AliESDtrack::kITSrefit)) continue;
1353         if(track->GetNcls(0)<4) continue;
1354         if(!track->HasPointOnITSLayer(0) && !track->HasPointOnITSLayer(1)) continue;
1355
1356         // Check matching with TPC only track
1357         Bool_t hasMatch = kFALSE;
1358         for(Int_t j=0; j<entries;++j)
1359         {
1360           if(i==j) continue;
1361           AliESDtrack *track2 = (AliESDtrack*)allChargedTracks->At(j);
1362           if(!track2) continue;
1363           if(track2->Charge()==0) continue;
1364           if (!track2->GetTPCInnerParam()) continue;
1365
1366           // check loose cuts for TPC tracks
1367           if(!esdTrackCuts->AcceptTrack(track2)) continue;
1368
1369           AliExternalTrackParam *innerTPC = new AliExternalTrackParam(*(track2->GetTPCInnerParam()));
1370           if(!innerTPC) continue;
1371           Double_t x2[3]; track2->GetXYZ(x2);
1372           Double_t b2[3]; AliTracker::GetBxByBz(x2,b2);
1373           Double_t dz[2],cov[3];
1374           Bool_t isPropOK = innerTPC->PropagateToDCABxByBz(vtxESD,b2,kVeryBig,dz,cov);
1375           if(!isPropOK && innerTPC) {delete innerTPC; continue;}
1376
1377           // check matching
1378           if (TMath::Abs(track->GetY() - innerTPC->GetY()) > 3) { delete innerTPC; continue; }
1379           if (TMath::Abs(track->GetSnp() - innerTPC->GetSnp()) > 0.2) { delete innerTPC; continue; }
1380           if (TMath::Abs(track->GetTgl() - innerTPC->GetTgl()) > 0.2) { delete innerTPC; continue; }
1381
1382           hasMatch = kTRUE;
1383           if(innerTPC) delete innerTPC;
1384         }
1385
1386         if(!hasMatch) continue;
1387       }
1388       else {
1389         // check track cuts
1390
1391         if(!esdTrackCuts->AcceptTrack(track))
1392           continue;
1393       }
1394
1395       //
1396       Bool_t isOK = kFALSE;
1397       Double_t x[3]; track->GetXYZ(x);
1398       Double_t b[3]; AliTracker::GetBxByBz(x,b);
1399
1400       //
1401       // if TPC-ITS hybrid tracking (kTPCITSHybrid)
1402       // replace track parameters with TPC-ony track parameters
1403       //
1404       if( GetAnalysisMode() == AlidNdPtHelper::kTPCITSHybrid ||
1405           GetAnalysisMode() == AlidNdPtHelper::kTPCITSHybridTrackSPDvtx ||
1406           GetAnalysisMode() == AlidNdPtHelper::kTPCITSHybridTrackSPDvtxDCArPt)
1407       {
1408         // Relate TPC-only tracks to Track or SPD vertex
1409         isOK = track->RelateToVertexTPCBxByBz(vtxESD, b, kVeryBig);
1410         if(!isOK) continue;
1411
1412         // replace esd track parameters with TPCinner
1413         AliExternalTrackParam  *tpcTrack  = new AliExternalTrackParam(*(track->GetTPCInnerParam()));
1414         if (tpcTrack) {
1415           track->Set(tpcTrack->GetX(),tpcTrack->GetAlpha(),tpcTrack->GetParameter(),tpcTrack->GetCovariance());
1416         }
1417         if(tpcTrack) delete tpcTrack;
1418       }
1419
1420
1421          // update track parameters using vertex point
1422          if( GetAnalysisMode() == AlidNdPtHelper::kTPCSPDvtxUpdate ||
1423              GetAnalysisMode() == AlidNdPtHelper::kTPCTrackSPDvtxUpdate )
1424          {
1425            // update track parameters
1426              AliExternalTrackParam cParam;
1427              isOK = track->RelateToVertexTPCBxByBz(vtxESD, b, kVeryBig, &cParam);
1428              if(!isOK) continue;
1429              track->Set(cParam.GetX(),cParam.GetAlpha(),cParam.GetParameter(),cParam.GetCovariance());
1430
1431              if(accCuts->AcceptTrack(track)) {
1432                FillHistograms(track,stack,vtxESD->GetZv(),AlidNdPtHelper::kRecTracks,multRecMult);
1433                labelsRec[multRec] = TMath::Abs(track->GetLabel());
1434                multRec++;
1435              }
1436           }
1437           else {
1438              if(accCuts->AcceptTrack(track))
1439              {
1440                FillHistograms(track,stack,vtxESD->GetZv(),AlidNdPtHelper::kRecTracks,multRecMult);
1441                labelsRec[multRec] = TMath::Abs(track->GetLabel());
1442                multRec++;
1443              }
1444           }
1445      }
1446
1447
1448      // fill track multiplicity histograms
1449      // terribly slow
1450      // FillHistograms(allChargedTracks,labelsAll,multAll,labelsAcc,multAcc,labelsRec,multRec);
1451
1452      Double_t vRecEventHist1[3] = {vtxESD->GetXv(),vtxESD->GetYv(),vtxESD->GetZv()};
1453      fRecEventHist1->Fill(vRecEventHist1);
1454
1455      Double_t vRecEventHist2[3] = {vtxESD->GetZv(),multMBTracks,multRecMult};
1456      fRecEventHist2->Fill(vRecEventHist2);
1457
1458      //
1459      Double_t vRecEventHist[2] = {vtxESD->GetZv(),multMBTracks};
1460      fRecEventHist->Fill(vRecEventHist);
1461    }
1462
1463    if(IsUseMCInfo())
1464    {
1465      if(mcEvent) {
1466
1467      if(evtCuts->IsEventSelectedRequired())
1468      {
1469        // select events with at least
1470        // one MC primary track in acceptance
1471        // pT>0.5 GeV/c, |eta|<0.8 for the Cross Section studies
1472
1473        Bool_t isMCEventSelected = AlidNdPtHelper::SelectMCEvent(mcEvent);
1474        //printf("isMCEventSelected %d \n", isMCEventSelected);
1475        if(!isMCEventSelected) {
1476
1477         if(allChargedTracks) delete allChargedTracks; allChargedTracks = 0;
1478         if(labelsAll) delete [] labelsAll; labelsAll = 0;
1479         if(labelsAcc) delete [] labelsAcc; labelsAcc = 0;
1480         if(labelsRec) delete [] labelsRec; labelsRec = 0;
1481
1482         return;
1483        }
1484      }
1485
1486      Double_t vMultTrueEventMatrix[3] = { multRecMult, multMCTrueTracks, multMBTracks};
1487      if(isEventOK && isEventTriggered) {
1488        if(TMath::Abs(vtxMC[2]) < 10.0) // both Rec. and corresponding MC events must be accepted
1489          fEventMultCorrelationMatrix->Fill(vMultTrueEventMatrix);
1490      }
1491
1492      //
1493      // event level corrections (zv,N_MB)
1494      //
1495
1496      // all inelastic
1497      Double_t vEventMatrix[2] = {vtxMC[2],multMCTrueTracks};
1498      fGenEventMatrix->Fill(vEventMatrix);
1499      if(isEventTriggered) fTriggerEventMatrix->Fill(vEventMatrix);
1500      if(isEventOK && isEventTriggered) fRecEventMatrix->Fill(vEventMatrix);
1501
1502      // single diffractive
1503      if(evtType == AliPWG0Helper::kSD) {
1504        fGenSDEventMatrix->Fill(vEventMatrix);
1505        if(isEventTriggered) fTriggerSDEventMatrix->Fill(vEventMatrix);
1506        if(isEventOK && isEventTriggered) fRecSDEventMatrix->Fill(vEventMatrix);
1507      }
1508
1509      // double diffractive
1510      if(evtType == AliPWG0Helper::kDD) {
1511        fGenDDEventMatrix->Fill(vEventMatrix);
1512        if(isEventTriggered) fTriggerDDEventMatrix->Fill(vEventMatrix);
1513        if(isEventOK && isEventTriggered)  fRecDDEventMatrix->Fill(vEventMatrix);
1514      }
1515
1516      // non diffractive
1517      if(evtType == AliPWG0Helper::kND) {
1518        fGenNDEventMatrix->Fill(vEventMatrix);
1519        if(isEventTriggered) fTriggerNDEventMatrix->Fill(vEventMatrix);
1520        if(isEventOK && isEventTriggered) fRecNDEventMatrix->Fill(vEventMatrix);
1521      }
1522
1523      // non single diffractive
1524      if(evtType != AliPWG0Helper::kSD) {
1525        fGenNSDEventMatrix->Fill(vEventMatrix);
1526        if(isEventTriggered) fTriggerNSDEventMatrix->Fill(vEventMatrix);
1527        if(isEventOK && isEventTriggered) fRecNSDEventMatrix->Fill(vEventMatrix);
1528      }
1529
1530      //
1531      // track-event level corrections (zv,pt,eta)
1532      //
1533      for (Int_t iMc = 0; iMc < stack->GetNtrack(); ++iMc)
1534      {
1535        TParticle* particle = stack->Particle(iMc);
1536        if (!particle)
1537        continue;
1538
1539        // only charged particles
1540        if(!particle->GetPDG()) continue;
1541        Double_t charge = particle->GetPDG()->Charge()/3.;
1542        if ( TMath::Abs(charge) < 0.001 )
1543         continue;
1544
1545        // only postive charged
1546        if(GetParticleMode() == AlidNdPtHelper::kPlus && charge < 0.)
1547         continue;
1548
1549        // only negative charged
1550        if(GetParticleMode() == AlidNdPtHelper::kMinus && charge > 0.)
1551        continue;
1552
1553        // physical primary
1554        Bool_t prim = stack->IsPhysicalPrimary(iMc);
1555        if(!prim) continue;
1556
1557        // checked accepted
1558        if(accCuts->AcceptTrack(particle))
1559        {
1560          Double_t vTrackEventMatrix[4] = {vtxMC[2], particle->Pt(), particle->Eta(),multMCTrueTracks};
1561          fGenTrackEventMatrix->Fill(vTrackEventMatrix);
1562
1563          if(evtType == AliPWG0Helper::kSD) {
1564            fGenTrackSDEventMatrix->Fill(vTrackEventMatrix);
1565          }
1566          if(evtType == AliPWG0Helper::kDD) {
1567            fGenTrackDDEventMatrix->Fill(vTrackEventMatrix);
1568          }
1569          if(evtType == AliPWG0Helper::kND) {
1570            fGenTrackNDEventMatrix->Fill(vTrackEventMatrix);
1571          }
1572          if(evtType != AliPWG0Helper::kSD) {
1573            fGenTrackNSDEventMatrix->Fill(vTrackEventMatrix);
1574          }
1575
1576          //
1577          if(!isEventTriggered) continue;
1578
1579          fTriggerTrackEventMatrix->Fill(vTrackEventMatrix);
1580          if(evtType == AliPWG0Helper::kSD) {
1581            fTriggerTrackSDEventMatrix->Fill(vTrackEventMatrix);
1582          }
1583          if(evtType == AliPWG0Helper::kDD) {
1584            fTriggerTrackDDEventMatrix->Fill(vTrackEventMatrix);
1585          }
1586          if(evtType == AliPWG0Helper::kND) {
1587            fTriggerTrackNDEventMatrix->Fill(vTrackEventMatrix);
1588          }
1589          if(evtType != AliPWG0Helper::kSD) {
1590            fTriggerTrackNSDEventMatrix->Fill(vTrackEventMatrix);
1591          }
1592
1593          //
1594          if(!isEventOK) continue;
1595
1596          fRecTrackEventMatrix->Fill(vTrackEventMatrix);
1597          if(evtType == AliPWG0Helper::kSD) {
1598            fRecTrackSDEventMatrix->Fill(vTrackEventMatrix);
1599          }
1600          if(evtType == AliPWG0Helper::kDD) {
1601            fRecTrackDDEventMatrix->Fill(vTrackEventMatrix);
1602          }
1603          if(evtType == AliPWG0Helper::kND) {
1604            fRecTrackNDEventMatrix->Fill(vTrackEventMatrix);
1605          }
1606          if(evtType != AliPWG0Helper::kSD) {
1607            fRecTrackNSDEventMatrix->Fill(vTrackEventMatrix);
1608          }
1609        }
1610      }
1611
1612      //
1613      // track-level corrections (zv,pt,eta)
1614      //
1615      if(isEventOK && isEventTriggered)
1616      {
1617
1618        // fill MC and rec event control histograms
1619        if(fHistogramsOn) {
1620          Double_t vRecMCEventHist1[3] = {vtxESD->GetXv()-vtxMC[0],vtxESD->GetYv()-vtxMC[1],vtxESD->GetZv()-vtxMC[2]};
1621          fRecMCEventHist1->Fill(vRecMCEventHist1);
1622
1623          Double_t vRecMCEventHist2[3] = {vtxESD->GetXv()-vtxMC[0],vtxESD->GetZv()-vtxMC[2],multMBTracks};
1624          fRecMCEventHist2->Fill(vRecMCEventHist2);
1625
1626          Double_t vRecMCEventHist3[2] = {multRecMult,evtType};
1627          fRecMCEventHist3->Fill(vRecMCEventHist3);
1628        }
1629
1630        //
1631        // MC histograms for track efficiency studies
1632        //
1633        Int_t countRecCandle = 0;
1634        for (Int_t iMc = 0; iMc < stack->GetNtrack(); ++iMc)
1635        {
1636          TParticle* particle = stack->Particle(iMc);
1637          if (!particle)
1638          continue;
1639
1640          Double_t vTrackMatrix[3] = {vtxMC[2],particle->Pt(),particle->Eta()};
1641
1642          // all genertated primaries including neutral
1643          //if( iMc < stack->GetNprimary() ) {
1644          //fGenTrackMatrix->Fill(vTrackMatrix);
1645          //}
1646
1647          // only charged particles
1648          if(!particle->GetPDG()) continue;
1649          Double_t charge = particle->GetPDG()->Charge()/3.;
1650          if (TMath::Abs(charge) < 0.001)
1651          continue;
1652
1653          // only postive charged
1654          if(GetParticleMode() == AlidNdPtHelper::kPlus && charge < 0.)
1655          continue;
1656
1657          // only negative charged
1658          if(GetParticleMode() == AlidNdPtHelper::kMinus && charge > 0.)
1659          continue;
1660
1661          // physical primary
1662          Bool_t prim = stack->IsPhysicalPrimary(iMc);
1663
1664          // check accepted
1665          if(accCuts->AcceptTrack(particle))
1666          {
1667            if( AlidNdPtHelper::IsPrimaryParticle(stack, iMc, GetParticleMode()) ) {
1668              fGenPrimTrackMatrix->Fill(vTrackMatrix);
1669              Double_t vMCPrimTrackHist[4] = {vtxMC[2],particle->Pt(),particle->Eta(),multMCTrueTracks};
1670              fMCPrimTrackHist->Fill(vMCPrimTrackHist);
1671            }
1672
1673            // fill control histograms
1674            if(fHistogramsOn)
1675              FillHistograms(stack,iMc,AlidNdPtHelper::kAccTracks);
1676
1677            // check multiple found tracks
1678            Int_t multCount = 0;
1679            for(Int_t iRec=0; iRec<multRec; ++iRec)
1680            {
1681              if(iMc == labelsRec[iRec])
1682              {
1683                multCount++;
1684                if(multCount>1)
1685                {
1686                  fRecMultTrackMatrix->Fill(vTrackMatrix);
1687
1688                  // fill control histogram
1689                  if(fHistogramsOn) {
1690                    Int_t pid = AlidNdPtHelper::ConvertPdgToPid(particle);
1691                    Double_t vMCMultRecTrackHist1[3] = {particle->Pt(), particle->Eta(), pid};
1692                    fMCMultRecTrackHist1->Fill(vMCMultRecTrackHist1);
1693                  }
1694                }
1695              }
1696            }
1697
1698            // check reconstructed
1699            for(Int_t iRec=0; iRec<multRec; ++iRec)
1700            {
1701              if(iMc == labelsRec[iRec])
1702              {
1703                fRecTrackMatrix->Fill(vTrackMatrix);
1704
1705                if( AlidNdPtHelper::IsPrimaryParticle(stack, iMc, GetParticleMode()) ) {
1706                  fRecPrimTrackMatrix->Fill(vTrackMatrix);
1707                  //AliESDtrack *track = esdEvent->GetTrack(iRec);
1708                  //if(track && track->GetKinkIndex(0) < 0)
1709                  //  printf("prim kink \n");
1710                  countRecCandle++;
1711                }
1712
1713                if(!prim) fRecSecTrackMatrix->Fill(vTrackMatrix);
1714
1715                // fill control histograms
1716                if(fHistogramsOn)
1717                  FillHistograms(stack,iMc,AlidNdPtHelper::kRecTracks);
1718
1719                break;
1720              }
1721            }
1722          }
1723        }
1724
1725        if(countRecCandle>0) fRecCandleEventMatrix->Fill(vEventMatrix);
1726       }
1727     }
1728   }// end bUseMC
1729
1730   if(allChargedTracks) delete allChargedTracks; allChargedTracks = 0;
1731   if(labelsAll) delete [] labelsAll; labelsAll = 0;
1732   if(labelsAcc) delete [] labelsAcc; labelsAcc = 0;
1733   if(labelsRec) delete [] labelsRec; labelsRec = 0;
1734
1735   if(!evtCuts->IsRecVertexRequired() && vtxESD != NULL) delete vtxESD;
1736   //if(trigAna) delete trigAna;
1737
1738 }
1739
1740 //_____________________________________________________________________________
1741 void AlidNdPtAnalysis::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) {
1742  // multiplicity  histograms
1743
1744   if(!allChargedTracks) return;
1745   if(!labelsAll) return;
1746   if(!labelsAcc) return;
1747   if(!labelsRec) return;
1748
1749   Int_t entries = allChargedTracks->GetEntries();
1750   for(Int_t i=0; i<entries; ++i)
1751   {
1752      AliESDtrack *track = (AliESDtrack*)allChargedTracks->At(i);
1753      if(!track) continue;
1754      if(track->Charge() == 0) continue;
1755
1756      Int_t label = TMath::Abs(track->GetLabel());
1757      for(Int_t iAll=0; iAll<multAll; ++iAll) {
1758        if(label == labelsAll[iAll]) {
1759          Double_t v1[2] = {track->Pt(), multAll};
1760          fRecTrackMultHist1[AlidNdPtHelper::kAllTracks]->Fill(v1);
1761        }
1762      }
1763      for(Int_t iAcc=0; iAcc<multAcc; ++iAcc) {
1764        if(label == labelsAcc[iAcc]) {
1765          Double_t v2[2] = {track->Pt(), multAcc};
1766          fRecTrackMultHist1[AlidNdPtHelper::kAccTracks]->Fill(v2);
1767        }
1768      }
1769      for(Int_t iRec=0; iRec<multRec; ++iRec) {
1770        if(label == labelsRec[iRec]) {
1771          Double_t v3[2] = {track->Pt(), multRec};
1772          fRecTrackMultHist1[AlidNdPtHelper::kRecTracks]->Fill(v3);
1773        }
1774      }
1775   }
1776 }
1777
1778 //_____________________________________________________________________________
1779 void AlidNdPtAnalysis::FillHistograms(AliESDtrack *const esdTrack, AliStack *const stack, const Double_t zv, AlidNdPtHelper::TrackObject trackObj, Int_t multRecMult)
1780 {
1781   //
1782   // Fill ESD track and MC histograms
1783   //
1784   if(!esdTrack) return;
1785
1786   Float_t q = esdTrack->Charge();
1787   if(TMath::Abs(q) < 0.001) return;
1788
1789   Float_t pt = esdTrack->Pt();
1790   //Float_t qpt = esdTrack->Pt() * q;
1791   Float_t eta = esdTrack->Eta();
1792   Float_t phi = esdTrack->Phi();
1793
1794   Float_t dca[2], bCov[3];
1795   esdTrack->GetImpactParameters(dca,bCov);
1796
1797   Int_t nClust = esdTrack->GetTPCclusters(0);
1798   Float_t chi2PerCluster = 0.;
1799   if(nClust>0.) chi2PerCluster = esdTrack->GetTPCchi2()/Float_t(nClust);
1800
1801
1802   // fill histograms
1803   Double_t values1[3] = {pt,eta,phi};
1804   fRecTrackHist1[trackObj]->Fill(values1);
1805
1806   Double_t values[4] = {zv, pt,eta, multRecMult};
1807   if(trackObj == AlidNdPtHelper::kRecTracks) {
1808     fRecTrackHist->Fill(values);
1809   }
1810
1811   /*
1812   Double_t values2[5] = {nClust,chi2PerCluster,pt,eta,phi};
1813   if(trackObj == AlidNdPtHelper::kRecTracks)
1814   {
1815     if(fHistogramsOn)
1816       fRecTrackHist2->Fill(values2);
1817   }
1818   */
1819
1820   //
1821   // Fill rec vs MC information
1822   //
1823   if(!stack) return;
1824
1825   Int_t label = TMath::Abs(esdTrack->GetLabel());
1826   //if(label == 0) return;
1827
1828   if(label > stack->GetNtrack()) return;
1829   TParticle* particle = stack->Particle(label);
1830   if(!particle) return;
1831
1832   //Bool_t prim = stack->IsPhysicalPrimary(label);
1833   //Int_t pid = AlidNdPtHelper::ConvertPdgToPid(particle);
1834
1835   Int_t motherPdg = -1;
1836   TParticle* mother = 0;
1837
1838   //TParticle* prim_mother = AlidNdPtHelper::FindPrimaryMother(stack,label);
1839   Int_t motherLabel = particle->GetMother(0);
1840   if(motherLabel>0) mother = stack->Particle(motherLabel);
1841   if(mother) motherPdg = TMath::Abs(mother->GetPdgCode()); // take abs for visualisation only
1842   //Int_t mech = particle->GetUniqueID(); // production mechanism
1843
1844   if(!particle->GetPDG()) return;
1845   Double_t gq = particle->GetPDG()->Charge()/3.0; // Charge units |e|/3
1846   if(TMath::Abs(gq)<0.001) return;
1847   Float_t gpt = particle->Pt();
1848   Float_t geta = particle->Eta();
1849   //Float_t qgpt = particle->Pt() * gq;
1850   //Float_t gphi = particle->Phi();
1851
1852   Double_t dpt=0;
1853   //printf("pt %f, gpt %f \n",pt,gpt);
1854   if(gpt) dpt = (pt-gpt)/gpt;
1855   Double_t deta = (eta-geta);
1856
1857   // fill histograms
1858   if(trackObj == AlidNdPtHelper::kRecTracks)
1859   {
1860     Double_t vTrackPtCorrelationMatrix[3]={pt,gpt,geta};
1861     fTrackPtCorrelationMatrix->Fill(vTrackPtCorrelationMatrix);
1862
1863     Double_t vRecMCTrackHist1[4]={gpt,geta,dpt,deta};
1864     fRecMCTrackHist1->Fill(vRecMCTrackHist1);
1865   }
1866 }
1867
1868 //_____________________________________________________________________________
1869 void AlidNdPtAnalysis::FillHistograms(AliStack *const stack, Int_t label, AlidNdPtHelper::TrackObject trackObj)
1870 {
1871   // Fill MC histograms
1872   if(!stack) return;
1873
1874   if(label > stack->GetNtrack()) return;
1875   TParticle* particle = stack->Particle(label);
1876   if(!particle) return;
1877
1878   Int_t motherPdg = -1;
1879   TParticle* mother = 0;
1880
1881   //TParticle* prim_mother = AlidNdPtHelper::FindPrimaryMother(stack,label);
1882   Int_t motherLabel = particle->GetMother(0);
1883   if(motherLabel>0) mother = stack->Particle(motherLabel);
1884   if(mother) motherPdg = TMath::Abs(mother->GetPdgCode()); // take abs for visualisation only
1885   Int_t mech = particle->GetUniqueID(); // production mechanism
1886
1887   if(!particle->GetPDG()) return;
1888   Double_t gq = particle->GetPDG()->Charge()/3.0; // Charge units |e|/3
1889   if(TMath::Abs(gq) < 0.001) return;
1890
1891   Float_t gpt = particle->Pt();
1892   //Float_t qgpt = particle->Pt() * gq;
1893   Float_t geta = particle->Eta();
1894   Float_t gphi = particle->Phi();
1895   //Float_t gpz = particle->Pz();
1896
1897   Bool_t prim = stack->IsPhysicalPrimary(label);
1898   //Float_t vx = particle->Vx(); Float_t vy = particle->Vy(); Float_t vz = particle->Vz();
1899
1900   Int_t pid = AlidNdPtHelper::ConvertPdgToPid(particle);
1901
1902   //if(prim&&pid==5) printf("pdgcode %d, production mech %d \n",particle->GetPdgCode(),mech);
1903   //if(!prim) printf("motherPdg %d, particle %d, production mech %d\n",motherPdg, particle->GetPdgCode(),mech);
1904
1905   //
1906   // fill histogram
1907   //
1908   Double_t vMCTrackHist1[3] = {gpt,geta,gphi};
1909   fMCTrackHist1[trackObj]->Fill(vMCTrackHist1);
1910
1911   Double_t vMCPrimTrackHist1[5] = {gpt,geta,pid,mech,motherPdg};
1912   Double_t vMCPrimTrackHist2[5] = {TMath::Abs(particle->GetPdgCode()),mech,motherPdg};
1913   //if(prim && AliPWG0Helper::IsPrimaryCharged(particle, label)) fMCPrimTrackHist1[trackObj]->Fill(vMCPrimTrackHist1);
1914   if(prim) {
1915     fMCPrimTrackHist1[trackObj]->Fill(vMCPrimTrackHist1);
1916     if(pid == 5) fMCPrimTrackHist2[trackObj]->Fill(vMCPrimTrackHist2);
1917   }
1918   else {
1919     fMCSecTrackHist1[trackObj]->Fill(vMCPrimTrackHist1);
1920   }
1921
1922 }
1923
1924 //_____________________________________________________________________________
1925 Long64_t AlidNdPtAnalysis::Merge(TCollection* const list)
1926 {
1927   //  init if not done already
1928   if (!fIsInit) { Init(); }
1929
1930   // Merge list of objects (needed by PROOF)
1931
1932   if (!list)
1933   return 0;
1934
1935   if (list->IsEmpty())
1936   return 1;
1937
1938   TIterator* iter = list->MakeIterator();
1939   TObject* obj = 0;
1940
1941   //
1942   //TList *collPhysSelection = new TList;
1943
1944   // collection of generated histograms
1945
1946   Int_t count=0;
1947   while((obj = iter->Next()) != 0) {
1948     AlidNdPtAnalysis* entry = dynamic_cast<AlidNdPtAnalysis*>(obj);
1949     if (entry == 0) continue;
1950
1951     // physics selection
1952     //printf("entry->GetPhysicsTriggerSelection() %p \n", entry->GetPhysicsTriggerSelection());
1953     //collPhysSelection->Add(entry->GetPhysicsTriggerSelection());
1954
1955     //
1956     fRecEventHist->Add(entry->fRecEventHist);
1957     fRecTrackHist->Add(entry->fRecTrackHist);
1958
1959     //
1960     fEventMultCorrelationMatrix->Add(entry->fEventMultCorrelationMatrix);
1961     fTrackPtCorrelationMatrix->Add(entry->fTrackPtCorrelationMatrix);
1962
1963     //
1964     fGenEventMatrix->Add(entry->fGenEventMatrix);
1965     fGenSDEventMatrix->Add(entry->fGenSDEventMatrix);
1966     fGenDDEventMatrix->Add(entry->fGenDDEventMatrix);
1967     fGenNDEventMatrix->Add(entry->fGenNDEventMatrix);
1968     fGenNSDEventMatrix->Add(entry->fGenNSDEventMatrix);
1969
1970     fTriggerEventMatrix->Add(entry->fTriggerEventMatrix);
1971     fTriggerSDEventMatrix->Add(entry->fTriggerSDEventMatrix);
1972     fTriggerDDEventMatrix->Add(entry->fTriggerDDEventMatrix);
1973     fTriggerNDEventMatrix->Add(entry->fTriggerNDEventMatrix);
1974     fTriggerNSDEventMatrix->Add(entry->fTriggerNSDEventMatrix);
1975
1976     fRecEventMatrix->Add(entry->fRecEventMatrix);
1977     fRecSDEventMatrix->Add(entry->fRecSDEventMatrix);
1978     fRecDDEventMatrix->Add(entry->fRecDDEventMatrix);
1979     fRecNDEventMatrix->Add(entry->fRecNDEventMatrix);
1980     fRecNSDEventMatrix->Add(entry->fRecNSDEventMatrix);
1981
1982     fRecCandleEventMatrix->Add(entry->fRecCandleEventMatrix);
1983     //
1984     fGenTrackEventMatrix->Add(entry->fGenTrackEventMatrix);
1985     fGenTrackSDEventMatrix->Add(entry->fGenTrackSDEventMatrix);
1986     fGenTrackDDEventMatrix->Add(entry->fGenTrackDDEventMatrix);
1987     fGenTrackNDEventMatrix->Add(entry->fGenTrackNDEventMatrix);
1988     fGenTrackNSDEventMatrix->Add(entry->fGenTrackNSDEventMatrix);
1989
1990     fTriggerTrackEventMatrix->Add(entry->fTriggerTrackEventMatrix);
1991     fTriggerTrackSDEventMatrix->Add(entry->fTriggerTrackSDEventMatrix);
1992     fTriggerTrackDDEventMatrix->Add(entry->fTriggerTrackDDEventMatrix);
1993     fTriggerTrackNDEventMatrix->Add(entry->fTriggerTrackNDEventMatrix);
1994     fTriggerTrackNSDEventMatrix->Add(entry->fTriggerTrackNSDEventMatrix);
1995
1996     fRecTrackEventMatrix->Add(entry->fRecTrackEventMatrix);
1997     fRecTrackSDEventMatrix->Add(entry->fRecTrackSDEventMatrix);
1998     fRecTrackDDEventMatrix->Add(entry->fRecTrackDDEventMatrix);
1999     fRecTrackNDEventMatrix->Add(entry->fRecTrackNDEventMatrix);
2000     fRecTrackNSDEventMatrix->Add(entry->fRecTrackNSDEventMatrix);
2001
2002     //
2003     fGenTrackMatrix->Add(entry->fGenTrackMatrix);
2004     fGenPrimTrackMatrix->Add(entry->fGenPrimTrackMatrix);
2005     fRecPrimTrackMatrix->Add(entry->fRecPrimTrackMatrix);
2006     //
2007     fRecTrackMatrix->Add(entry->fRecTrackMatrix);
2008     fRecSecTrackMatrix->Add(entry->fRecSecTrackMatrix);
2009     //
2010     fRecMultTrackMatrix->Add(entry->fRecMultTrackMatrix);
2011
2012     //
2013     // control analysis histograms
2014     //
2015     fMCEventHist1->Add(entry->fMCEventHist1);
2016     fRecEventHist1->Add(entry->fRecEventHist1);
2017     fRecEventHist2->Add(entry->fRecEventHist2);
2018     fRecMCEventHist1->Add(entry->fRecMCEventHist1);
2019     fRecMCEventHist2->Add(entry->fRecMCEventHist2);
2020     fRecMCEventHist3->Add(entry->fRecMCEventHist3);
2021
2022     for(Int_t i=0; i<AlidNdPtHelper::kCutSteps; i++) {
2023       fMCTrackHist1[i]->Add(entry->fMCTrackHist1[i]);
2024
2025       fMCPrimTrackHist1[i]->Add(entry->fMCPrimTrackHist1[i]);
2026       fMCPrimTrackHist2[i]->Add(entry->fMCPrimTrackHist2[i]);
2027       fMCSecTrackHist1[i]->Add(entry->fMCSecTrackHist1[i]);
2028
2029       fRecTrackHist1[i]->Add(entry->fRecTrackHist1[i]);
2030       fRecTrackMultHist1[i]->Add(entry->fRecTrackMultHist1[i]);
2031     }
2032     fRecMCTrackHist1->Add(entry->fRecMCTrackHist1);
2033     fMCMultRecTrackHist1->Add(entry->fMCMultRecTrackHist1);
2034     fRecTrackHist2->Add(entry->fRecTrackHist2);
2035     fEventCount->Add(entry->fEventCount);
2036     fMCPrimTrackHist->Add(entry->fMCPrimTrackHist);
2037
2038   count++;
2039   }
2040
2041   //AliPhysicsSelection *trigSelection = GetPhysicsTriggerSelection();
2042   //trigSelection->Merge(collPhysSelection);
2043   //if(collPhysSelection) delete collPhysSelection;
2044
2045 return count;
2046 }
2047
2048 //____________________________________________________________________
2049 Bool_t AlidNdPtAnalysis::IsBinZeroTrackSPDvtx(const AliESDEvent* esdEvent) {
2050 //
2051 // check 0-bin
2052 // for LHC background calculation
2053 // return kTRUE if vertex not reconstructed or
2054 // track multiplicity == 0
2055 //
2056 if(!esdEvent) return kFALSE;
2057
2058   // check vertex
2059   const AliESDVertex *vertex = esdEvent->GetPrimaryVertexTracks();
2060   if(!vertex) return kTRUE;
2061
2062   if(vertex->GetNContributors() < 1) {
2063     // SPD vertex
2064     vertex = esdEvent->GetPrimaryVertexSPD();
2065     if(!vertex) return kTRUE;
2066   }
2067   if(TMath::Abs(vertex->GetZv()) > 15.0) return kTRUE;
2068   if( !AlidNdPtHelper::TestRecVertex(vertex, esdEvent->GetPrimaryVertexSPD(), AlidNdPtHelper::kTPCITSHybridTrackSPDvtx) )
2069     return kTRUE;
2070
2071 return kFALSE;
2072 }
2073
2074 //____________________________________________________________________
2075 Bool_t AlidNdPtAnalysis::IsBinZeroSPDvtx(const AliESDEvent* esdEvent) {
2076 //
2077 // check 0-bin
2078 // for LHC background calculation
2079 // return kTRUE if vertex not reconstructed or
2080 // tracklet multiplicity == 0
2081 //
2082 if(!esdEvent) return kFALSE;
2083
2084   // check vertex
2085   const AliESDVertex* vertex = esdEvent->GetPrimaryVertexSPD();
2086   if(!vertex) return kTRUE;
2087   if(TMath::Abs(vertex->GetZv()) > 15.0) return kTRUE;
2088   if( !AlidNdPtHelper::TestRecVertex(vertex, esdEvent->GetPrimaryVertexSPD(), AlidNdPtHelper::kTPCSPDvtx) ) return kTRUE;
2089
2090 return kFALSE;
2091 }
2092
2093 //_____________________________________________________________________________
2094 void AlidNdPtAnalysis::Analyse()
2095 {
2096   //  init if not done already
2097   if (!fIsInit) { Init(); }
2098
2099   // Analyse histograms
2100   //
2101   TH1::AddDirectory(kFALSE);
2102   TH1 *h=0, *h1=0, *h2=0, *h2c = 0;
2103   THnSparse *hs=0;
2104   TH2 *h2D=0;
2105
2106   char name[256];
2107   TObjArray *aFolderObj = new TObjArray;
2108   if(!aFolderObj) return;
2109
2110   //
2111   // LHC backgraund in all and 0-bins
2112   //
2113   //AliPhysicsSelection *trigSel = GetPhysicsTriggerSelection();
2114   //trigSel->SaveHistograms("physics_selection");
2115
2116   //
2117   // Reconstructed event vertex
2118   //
2119   h = fRecEventHist1->Projection(0);
2120   h->SetName("Xv");
2121   aFolderObj->Add(h);
2122
2123   h = fRecEventHist1->Projection(1);
2124   h->SetName("Yv");
2125   aFolderObj->Add(h);
2126
2127   h = fRecEventHist1->Projection(2);
2128   h->SetName("Zv");
2129   aFolderObj->Add(h);
2130
2131   //
2132   // multiplicity
2133   //
2134   h = fRecEventHist2->Projection(1);
2135   h->SetName("multMB");
2136   aFolderObj->Add(h);
2137
2138   h = fRecEventHist2->Projection(2);
2139   h->SetName("multiplicity");
2140   aFolderObj->Add(h);
2141
2142   h2D = fRecEventHist2->Projection(0,1);
2143   h2D->SetName("Zv_vs_multiplicity_MB");
2144   aFolderObj->Add(h2D);
2145
2146   //
2147   // reconstructed pt histograms
2148   //
2149   h = fRecTrackHist1[0]->Projection(0);
2150   h->Scale(1.,"width");
2151   h->SetName("pt_all_ch");
2152   aFolderObj->Add(h);
2153
2154   h = fRecTrackHist1[1]->Projection(0);
2155   h->Scale(1.,"width");
2156   h->SetName("pt_acc");
2157   aFolderObj->Add(h);
2158
2159   h = fRecTrackHist1[2]->Projection(0);
2160   h->Scale(1.,"width");
2161   h->SetName("pt_rec");
2162   aFolderObj->Add(h);
2163
2164   //
2165   // reconstructed eta histograms
2166   //
2167   h = fRecTrackHist1[0]->Projection(1);
2168   h->SetName("eta_all_ch");
2169   aFolderObj->Add(h);
2170
2171   h = fRecTrackHist1[1]->Projection(1);
2172   h->SetName("eta_acc");
2173   aFolderObj->Add(h);
2174
2175   h = fRecTrackHist1[2]->Projection(1);
2176   h->SetName("eta_rec");
2177   aFolderObj->Add(h);
2178
2179   //
2180   // reconstructed phi histograms
2181   //
2182   h = fRecTrackHist1[0]->Projection(2);
2183   h->SetName("phi_all_ch");
2184   aFolderObj->Add(h);
2185
2186   h = fRecTrackHist1[1]->Projection(2);
2187   h->SetName("phi_acc");
2188   aFolderObj->Add(h);
2189
2190   h = fRecTrackHist1[2]->Projection(2);
2191   h->SetName("phi_rec");
2192   aFolderObj->Add(h);
2193
2194   //
2195   // reconstructed eta:pt histograms
2196   //
2197   h2D = fRecTrackHist1[0]->Projection(1,0);
2198   h2D->SetName("pt_eta_all_ch");
2199   aFolderObj->Add(h2D);
2200
2201   h2D = fRecTrackHist1[1]->Projection(1,0);
2202   h2D->SetName("pt_eta_acc");
2203   aFolderObj->Add(h2D);
2204
2205   h2D = fRecTrackHist1[2]->Projection(1,0);
2206   h2D->SetName("pt_eta_rec");
2207   aFolderObj->Add(h2D);
2208
2209   //
2210   // reconstructed phi:pt histograms
2211   //
2212   h2D = fRecTrackHist1[0]->Projection(2,0);
2213   h2D->SetName("pt_phi_all_ch");
2214   aFolderObj->Add(h2D);
2215
2216   h2D = fRecTrackHist1[1]->Projection(2,0);
2217   h2D->SetName("pt_phi_acc");
2218   aFolderObj->Add(h2D);
2219
2220   h2D = fRecTrackHist1[2]->Projection(2,0);
2221   h2D->SetName("pt_phi_rec");
2222   aFolderObj->Add(h2D);
2223
2224   //
2225   // reconstructed phi:eta histograms
2226   //
2227   h2D = fRecTrackHist1[0]->Projection(2,1);
2228   h2D->SetName("eta_phi_all_ch");
2229   aFolderObj->Add(h2D);
2230
2231   h2D = fRecTrackHist1[1]->Projection(2,1);
2232   h2D->SetName("eta_phi_acc");
2233   aFolderObj->Add(h2D);
2234
2235   h2D = fRecTrackHist1[2]->Projection(2,1);
2236   h2D->SetName("eta_phi_rec");
2237   aFolderObj->Add(h2D);
2238
2239   //
2240   // reconstructed nClust, chi2 vs pt, eta, phi
2241   //
2242   if(fHistogramsOn) {
2243
2244     h2D = fRecTrackHist2->Projection(0,1);
2245     h2D->SetName("nClust_chi2_rec");
2246     aFolderObj->Add(h2D);
2247
2248     h2D = fRecTrackHist2->Projection(0,2);
2249     h2D->SetName("nClust_pt_rec");
2250     aFolderObj->Add(h2D);
2251
2252     h2D = fRecTrackHist2->Projection(0,3);
2253     h2D->SetName("nClust_eta_rec");
2254     aFolderObj->Add(h2D);
2255
2256     h2D = fRecTrackHist2->Projection(0,4);
2257     h2D->SetName("nClust_phi_rec");
2258     aFolderObj->Add(h2D);
2259
2260     h2D = fRecTrackHist2->Projection(1,2);
2261     h2D->SetName("chi2_pt_rec");
2262     aFolderObj->Add(h2D);
2263
2264     h2D = fRecTrackHist2->Projection(1,3);
2265     h2D->SetName("chi2_eta_rec");
2266     aFolderObj->Add(h2D);
2267
2268     h2D = fRecTrackHist2->Projection(1,4);
2269     h2D->SetName("chi2_phi_rec");
2270     aFolderObj->Add(h2D);
2271
2272   }
2273
2274   //
2275   // calculate corrections for empty events
2276   // with multMB==0
2277   //
2278
2279   //
2280   // normalised zv to generate zv for triggered events
2281   //
2282   h = fRecEventHist2->Projection(0);
2283   if( h->Integral() ) h->Scale(1./h->Integral());
2284   h->SetName("zv_distribution_norm");
2285   aFolderObj->Add(h);
2286
2287   //
2288   // MC available
2289   //
2290   if(IsUseMCInfo()) {
2291
2292   //
2293   // Event vertex resolution
2294   //
2295   h2D = fRecMCEventHist2->Projection(0,2);
2296   h2D->SetName("DeltaXv_vs_mult");
2297   aFolderObj->Add(h2D);
2298
2299   h2D = fRecMCEventHist2->Projection(1,2);
2300   h2D->SetName("DeltaZv_vs_mult");
2301   aFolderObj->Add(h2D);
2302
2303   //
2304   // normalised zv to get trigger/trigger+vertex event differences
2305   // F(zv) = E_trig(zv,0)/Int(E_trig(zv,0) / Sum(E_trigvtx(zv,n))/Sum(Int(E_trigvtx(zv,n))dzv)
2306   //
2307   fTriggerEventMatrix->GetAxis(1)->SetRangeUser(0.,0.);
2308   h = fTriggerEventMatrix->Projection(0);
2309   h2D = fTriggerEventMatrix->Projection(0,1);
2310   if(h2D->Integral()) h->Scale(1./h2D->Integral());
2311
2312   h1 = fRecEventMatrix->Projection(0);
2313   h2D = fRecEventMatrix->Projection(0,1);
2314   if(h2D->Integral()) h1->Scale(1./h2D->Integral());
2315
2316   h->Divide(h1);
2317   h->SetName("zv_empty_events_norm");
2318   aFolderObj->Add(h);
2319
2320   fTriggerEventMatrix->GetAxis(1)->SetRange(1,fTriggerEventMatrix->GetAxis(1)->GetNbins());
2321
2322   //
2323   // rec. vs true multiplicity correlation matrix
2324   //
2325   hs = (THnSparse*)fEventMultCorrelationMatrix->Clone("event_mult_correlation_matrix");
2326   aFolderObj->Add(hs);
2327
2328   //
2329   // rec. vs true track pt correlation matrix
2330   //
2331   hs = (THnSparse*)fTrackPtCorrelationMatrix->Clone("track_pt_correlation_matrix");
2332   aFolderObj->Add(hs);
2333
2334   //
2335   // trigger efficiency for INEL
2336   //
2337   h = AlidNdPtHelper::GenerateCorrMatrix(fTriggerEventMatrix->Projection(0),fGenEventMatrix->Projection(0),"zv_trig_INEL_eff_matrix");
2338   aFolderObj->Add(h);
2339
2340   //
2341   // trigger efficiency for NSD
2342   //
2343   h = AlidNdPtHelper::GenerateCorrMatrix(fTriggerNSDEventMatrix->Projection(0),fGenNSDEventMatrix->Projection(0),"zv_trig_NSD_eff_matrix");
2344   aFolderObj->Add(h);
2345
2346   //
2347   // trigger bias correction (MB to ND)
2348   //
2349   hs = AlidNdPtHelper::GenerateCorrMatrix(fGenNDEventMatrix,fTriggerEventMatrix,"zv_mult_trig_MBtoND_corr_matrix");
2350   aFolderObj->Add(hs);
2351
2352   h = AlidNdPtHelper::GenerateCorrMatrix(fGenNDEventMatrix->Projection(0),fTriggerEventMatrix->Projection(0),"zv_trig_MBtoND_corr_matrix");
2353   aFolderObj->Add(h);
2354
2355
2356   h = AlidNdPtHelper::GenerateCorrMatrix(fGenNDEventMatrix->Projection(1),fTriggerEventMatrix->Projection(1),"mult_trig_MBtoND_corr_matrix");
2357
2358   aFolderObj->Add(h);
2359
2360   //
2361   // trigger bias correction (MB to NSD)
2362   //
2363   hs = AlidNdPtHelper::GenerateCorrMatrix(fGenNSDEventMatrix,fTriggerEventMatrix,"zv_mult_trig_MBtoNSD_corr_matrix");
2364   aFolderObj->Add(hs);
2365
2366   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenNSDEventMatrix->Projection(0,1),fTriggerEventMatrix->Projection(0,1),"zv_mult_trig_MBtoNSD_corr_matrix_2D");
2367   aFolderObj->Add(h2D);
2368
2369   h = AlidNdPtHelper::GenerateCorrMatrix(fGenNSDEventMatrix->Projection(0),fTriggerEventMatrix->Projection(0),"zv_trig_MBtoNSD_corr_matrix");
2370   aFolderObj->Add(h);
2371
2372
2373   h = AlidNdPtHelper::GenerateCorrMatrix(fGenNSDEventMatrix->Projection(1),fTriggerEventMatrix->Projection(1),"mult_trig_MBtoNSD_corr_matrix");
2374   aFolderObj->Add(h);
2375
2376
2377   //
2378   // trigger bias correction (MB to INEL)
2379   //
2380   hs = AlidNdPtHelper::GenerateCorrMatrix(fGenEventMatrix,fTriggerEventMatrix,"zv_mult_trig_MBtoInel_corr_matrix");
2381   aFolderObj->Add(hs);
2382
2383   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenEventMatrix->Projection(0,1),fTriggerEventMatrix->Projection(0,1),"zv_mult_trig_MBtoInel_corr_matrix_2D");
2384   aFolderObj->Add(h2D);
2385
2386   h = AlidNdPtHelper::GenerateCorrMatrix(fGenEventMatrix->Projection(0),fTriggerEventMatrix->Projection(0),"zv_trig_MBtoInel_corr_matrix");
2387   aFolderObj->Add(h);
2388
2389   h = AlidNdPtHelper::GenerateCorrMatrix(fGenEventMatrix->Projection(1),fTriggerEventMatrix->Projection(1),"mult_trig_MBtoInel_corr_matrix");
2390   aFolderObj->Add(h);
2391
2392
2393   //
2394   // event vertex reconstruction correction (MB)
2395   //
2396   hs = AlidNdPtHelper::GenerateCorrMatrix(fTriggerEventMatrix,fRecEventMatrix,"zv_mult_event_corr_matrix");
2397   aFolderObj->Add(hs);
2398
2399   h2D = AlidNdPtHelper::GenerateCorrMatrix(fTriggerEventMatrix->Projection(0,1),fRecEventMatrix->Projection(0,1),"zv_mult_event_corr_matrix_2D");
2400   aFolderObj->Add(h2D);
2401
2402
2403   h = AlidNdPtHelper::GenerateCorrMatrix(fTriggerEventMatrix->Projection(1),fRecEventMatrix->Projection(1),"mult_event_corr_matrix");
2404   aFolderObj->Add(h);
2405
2406
2407   h = AlidNdPtHelper::GenerateCorrMatrix(fTriggerEventMatrix->Projection(0),fRecEventMatrix->Projection(0),"zv_event_corr_matrix");
2408   aFolderObj->Add(h);
2409
2410   //
2411   // track-event trigger bias correction (MB to ND)
2412   //
2413
2414   hs = AlidNdPtHelper::GenerateCorrMatrix(fGenTrackNDEventMatrix,fTriggerTrackEventMatrix,"zv_pt_eta_track_trig_MBtoND_corr_matrix");
2415   aFolderObj->Add(hs);
2416
2417   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenTrackNDEventMatrix->Projection(1,2),fTriggerTrackEventMatrix->Projection(1,2),"eta_pt_track_trig_MBtoND_corr_matrix");
2418   aFolderObj->Add(h2D);
2419
2420   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenTrackNDEventMatrix->Projection(1,0),fTriggerTrackEventMatrix->Projection(1,0),"pt_zv_track_trig_MBtoND_corr_matrix");
2421   aFolderObj->Add(h2D);
2422
2423   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenTrackNDEventMatrix->Projection(2,0),fTriggerTrackEventMatrix->Projection(2,0),"zv_eta_track_trig_MBtoND_corr_matrix");
2424   aFolderObj->Add(h2D);
2425
2426   //
2427   // track-event trigger bias correction (MB to NSD)
2428   //
2429   hs = AlidNdPtHelper::GenerateCorrMatrix(fGenTrackNSDEventMatrix,fTriggerTrackEventMatrix,"zv_pt_eta_track_trig_MBtoNSD_corr_matrix");
2430   aFolderObj->Add(hs);
2431
2432   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenTrackNSDEventMatrix->Projection(1,2),fTriggerTrackEventMatrix->Projection(1,2),"eta_pt_track_trig_MBtoNSD_corr_matrix");
2433   aFolderObj->Add(h2D);
2434
2435   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenTrackNSDEventMatrix->Projection(1,0),fTriggerTrackEventMatrix->Projection(1,0),"pt_zv_track_trig_MBtoNSD_corr_matrix");
2436   aFolderObj->Add(h2D);
2437
2438   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenTrackNSDEventMatrix->Projection(2,0),fTriggerTrackEventMatrix->Projection(2,0),"zv_eta_track_trig_MBtoNSD_corr_matrix");
2439   aFolderObj->Add(h2D);
2440
2441
2442   //
2443   // track-event trigger bias correction (MB to INEL)
2444   //
2445   hs = AlidNdPtHelper::GenerateCorrMatrix(fGenTrackEventMatrix,fTriggerTrackEventMatrix,"zv_pt_eta_track_trig_MBtoInel_corr_matrix");
2446   aFolderObj->Add(hs);
2447
2448   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenTrackEventMatrix->Projection(1,2),fTriggerTrackEventMatrix->Projection(1,2),"eta_pt_track_trig_MBtoInel_corr_matrix");
2449   aFolderObj->Add(h2D);
2450
2451   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenTrackEventMatrix->Projection(1,0),fTriggerTrackEventMatrix->Projection(1,0),"pt_zv_track_trig_MBtoInel_corr_matrix");
2452   aFolderObj->Add(h2D);
2453
2454   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenTrackEventMatrix->Projection(2,0),fTriggerTrackEventMatrix->Projection(2,0),"zv_eta_track_trig_MBtoInel_corr_matrix");
2455   aFolderObj->Add(h2D);
2456
2457   // efficiency
2458
2459   h = AlidNdPtHelper::GenerateCorrMatrix(fTriggerTrackEventMatrix->Projection(1),fGenTrackEventMatrix->Projection(1),"pt_track_trig_MBtoInel_eff_matrix");
2460   aFolderObj->Add(h);
2461
2462
2463   //
2464   // track-event vertex reconstruction correction (MB)
2465   //
2466   hs = AlidNdPtHelper::GenerateCorrMatrix(fTriggerTrackEventMatrix,fRecTrackEventMatrix,"zv_pt_eta_track_event_corr_matrix");
2467   aFolderObj->Add(hs);
2468
2469   h2D = AlidNdPtHelper::GenerateCorrMatrix(fTriggerTrackEventMatrix->Projection(1,2),fRecTrackEventMatrix->Projection(1,2),"eta_pt_track_event_corr_matrix");
2470   aFolderObj->Add(h2D);
2471
2472   h2D = AlidNdPtHelper::GenerateCorrMatrix(fTriggerTrackEventMatrix->Projection(1,0),fRecTrackEventMatrix->Projection(1,0),"pt_zv_track_event_corr_matrix");
2473   aFolderObj->Add(h2D);
2474
2475   h2D = AlidNdPtHelper::GenerateCorrMatrix(fTriggerTrackEventMatrix->Projection(2,0),fRecTrackEventMatrix->Projection(2,0),"zv_eta_track_event_corr_matrix");
2476   aFolderObj->Add(h2D);
2477
2478   // efficiency
2479
2480   h = AlidNdPtHelper::GenerateCorrMatrix(fRecTrackEventMatrix->Projection(1),fTriggerTrackEventMatrix->Projection(1),"pt_track_event_eff_matrix");
2481   aFolderObj->Add(h);
2482
2483
2484   //
2485   // track rec. efficiency correction
2486   //
2487   hs = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix,fRecPrimTrackMatrix,"zv_pt_eta_track_corr_matrix");
2488   aFolderObj->Add(hs);
2489
2490   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(1,2),fRecPrimTrackMatrix->Projection(1,2),"eta_pt_track_corr_matrix");
2491   aFolderObj->Add(h2D);
2492
2493   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(1,0),fRecPrimTrackMatrix->Projection(1,0),"pt_zv_track_corr_matrix");
2494   aFolderObj->Add(h2D);
2495
2496   h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(2,0),fRecPrimTrackMatrix->Projection(2,0),"zv_eta_track_corr_matrix");
2497   aFolderObj->Add(h2D);
2498
2499
2500   h = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(0),fRecPrimTrackMatrix->Projection(0),"zv_track_corr_matrix");
2501   aFolderObj->Add(h);
2502
2503   h = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(1),fRecPrimTrackMatrix->Projection(1),"pt_track_corr_matrix");
2504   aFolderObj->Add(h);
2505
2506   // efficiency
2507
2508   h = AlidNdPtHelper::GenerateCorrMatrix(fRecPrimTrackMatrix->Projection(1), fGenPrimTrackMatrix->Projection(1),"pt_track_eff_matrix");
2509   aFolderObj->Add(h);
2510
2511   h = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(2),fRecPrimTrackMatrix->Projection(2),"eta_track_corr_matrix");
2512   aFolderObj->Add(h);
2513
2514   //
2515   // secondary track contamination correction
2516   //
2517   //hs = AlidNdPtHelper::GenerateContCorrMatrix(fRecSecTrackMatrix,fRecTrackMatrix,"zv_pt_eta_track_cont_matrix");
2518   hs = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix,fRecTrackMatrix,"zv_pt_eta_track_cont_matrix");
2519   aFolderObj->Add(hs);
2520
2521   h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(1,2),fRecTrackMatrix->Projection(1,2),"eta_pt_track_cont_matrix");
2522   aFolderObj->Add(h2D);
2523
2524   h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(1,0),fRecTrackMatrix->Projection(1,0),"pt_zv_track_cont_matrix");
2525   aFolderObj->Add(h2D);
2526
2527   h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(2,0),fRecTrackMatrix->Projection(2,0),"zv_eta_track_cont_matrix");
2528   aFolderObj->Add(h2D);
2529
2530   h = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(0),fRecTrackMatrix->Projection(0),"zv_track_cont_matrix");
2531   aFolderObj->Add(h);
2532
2533
2534   h = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(1),fRecTrackMatrix->Projection(1),"pt_track_cont_matrix");
2535   aFolderObj->Add(h);
2536
2537   h = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(2),fRecTrackMatrix->Projection(2),"eta_track_cont_matrix");
2538   aFolderObj->Add(h);
2539
2540   //
2541   // multiple track reconstruction correction
2542   //
2543   //hs = AlidNdPtHelper::GenerateContCorrMatrix(fRecMultTrackMatrix,fRecTrackMatrix,"zv_pt_eta_mult_track_cont_matrix");
2544   hs = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix,fRecTrackMatrix,"zv_pt_eta_mult_track_cont_matrix");
2545   aFolderObj->Add(hs);
2546
2547   h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(1,2),fRecTrackMatrix->Projection(1,2),"eta_pt_mult_track_cont_matrix");
2548   aFolderObj->Add(h2D);
2549
2550   h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(1,0),fRecTrackMatrix->Projection(1,0),"pt_zv_mult_track_cont_matrix");
2551   aFolderObj->Add(h2D);
2552
2553   h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(2,0),fRecTrackMatrix->Projection(2,0),"zv_eta_mult_track_cont_matrix");
2554   aFolderObj->Add(h2D);
2555
2556   h = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(0),fRecTrackMatrix->Projection(0),"zv_mult_track_cont_matrix");
2557   aFolderObj->Add(h);
2558
2559   h = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(1),fRecTrackMatrix->Projection(1),"pt_mult_track_cont_matrix");
2560   aFolderObj->Add(h);
2561
2562   h = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(2),fRecTrackMatrix->Projection(2),"eta_mult_track_cont_matrix");
2563   aFolderObj->Add(h);
2564
2565   //
2566   // Control histograms
2567   //
2568
2569   if(fHistogramsOn) {
2570
2571   // Efficiency electrons, muons, pions, kaons, protons, all
2572   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(1,1);
2573   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(1,1);
2574   h1 = fMCPrimTrackHist1[1]->Projection(0);
2575   h2 = fMCPrimTrackHist1[2]->Projection(0);
2576   h2c = (TH1D *)h2->Clone();
2577   h2c->Divide(h1);
2578   h2c->SetName("eff_pt_electrons");
2579   aFolderObj->Add(h2c);
2580
2581   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(2,2);
2582   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(2,2);
2583   h1 = fMCPrimTrackHist1[1]->Projection(0);
2584   h2 = fMCPrimTrackHist1[2]->Projection(0);
2585   h2c = (TH1D *)h2->Clone();
2586   h2c->Divide(h1);
2587   h2c->SetName("eff_pt_muons");
2588   aFolderObj->Add(h2c);
2589
2590   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(3,3);
2591   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(3,3);
2592   h1 = fMCPrimTrackHist1[1]->Projection(0);
2593   h2 = fMCPrimTrackHist1[2]->Projection(0);
2594   h2c = (TH1D *)h2->Clone();
2595   h2c->Divide(h1);
2596   h2c->SetName("eff_pt_pions");
2597   aFolderObj->Add(h2c);
2598
2599   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(4,4);
2600   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(4,4);
2601   h1 = fMCPrimTrackHist1[1]->Projection(0);
2602   h2 = fMCPrimTrackHist1[2]->Projection(0);
2603   h2c = (TH1D *)h2->Clone();
2604   h2c->Divide(h1);
2605   h2c->SetName("eff_pt_kaons");
2606   aFolderObj->Add(h2c);
2607
2608   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(5,5);
2609   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(5,5);
2610   h1 = fMCPrimTrackHist1[1]->Projection(0);
2611   h2 = fMCPrimTrackHist1[2]->Projection(0);
2612   h2c = (TH1D *)h2->Clone();
2613   h2c->Divide(h1);
2614   h2c->SetName("eff_pt_protons");
2615   aFolderObj->Add(h2c);
2616
2617   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(1,5);
2618   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(1,5);
2619   h1 = fMCPrimTrackHist1[1]->Projection(0);
2620   h2 = fMCPrimTrackHist1[2]->Projection(0);
2621   h2c = (TH1D *)h2->Clone();
2622   h2c->Divide(h1);
2623   h2c->SetName("eff_pt_selected");
2624   aFolderObj->Add(h2c);
2625
2626   fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(1,6);
2627   fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(1,6);
2628   h1 = fMCPrimTrackHist1[1]->Projection(0);
2629   h2 = fMCPrimTrackHist1[2]->Projection(0);
2630   h2c = (TH1D *)h2->Clone();
2631   h2c->Divide(h1);
2632   h2c->SetName("eff_pt_all");
2633   aFolderObj->Add(h2c);
2634
2635   fMCPrimTrackHist1[1]->GetAxis(1)->SetRange(1,fMCPrimTrackHist1[1]->GetAxis(1)->GetNbins());
2636   fMCPrimTrackHist1[2]->GetAxis(1)->SetRange(1,fMCPrimTrackHist1[2]->GetAxis(1)->GetNbins());
2637
2638   //  pt spetra
2639   // - rec, primaries, secondaries
2640   // - primaries (pid)
2641   // - secondaries (pid)
2642   // - secondaries (mech)
2643   // - secondaries (mother)
2644   //
2645
2646   TH1D *mcPtAccall = fMCTrackHist1[1]->Projection(0);
2647   mcPtAccall->SetName("mc_pt_acc_all");
2648   aFolderObj->Add(mcPtAccall);
2649
2650   TH1D *mcPtAccprim = fMCPrimTrackHist1[1]->Projection(0);
2651   mcPtAccprim->SetName("mc_pt_acc_prim");
2652   aFolderObj->Add(mcPtAccprim);
2653
2654   TH1D *mcPtRecall = fMCTrackHist1[2]->Projection(0);
2655   mcPtRecall->SetName("mc_pt_rec_all");
2656   aFolderObj->Add(mcPtRecall);
2657
2658   TH1D *mcPtRecprim = fMCPrimTrackHist1[2]->Projection(0);
2659   mcPtRecprim->SetName("mc_pt_rec_prim");
2660   aFolderObj->Add(mcPtRecprim);
2661
2662   TH1D *mcPtRecsec = fMCSecTrackHist1[2]->Projection(0);
2663   mcPtRecsec->SetName("mc_pt_rec_sec");
2664   aFolderObj->Add(mcPtRecsec);
2665
2666   for(Int_t i = 0; i<6; i++)
2667   {
2668     snprintf(name,256,"mc_pt_rec_prim_pid_%d",i);
2669     fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(i+1,i+1);
2670     h = fMCPrimTrackHist1[2]->Projection(0);
2671     h->SetName(name);
2672     aFolderObj->Add(h);
2673
2674     snprintf(name,256,"mc_pt_rec_sec_pid_%d",i);
2675     fMCSecTrackHist1[2]->GetAxis(2)->SetRange(i+1,i+1);
2676     h = fMCSecTrackHist1[2]->Projection(0);
2677     h->SetName(name);
2678     aFolderObj->Add(h);
2679
2680     // production mechanisms for given pid
2681     fMCSecTrackHist1[2]->GetAxis(2)->SetRange(i+1,i+1);
2682
2683     for(Int_t j=0; j<20; j++) {
2684       if(j == 4) {
2685         // decay
2686
2687         snprintf(name,256,"mc_pt_rec_sec_pid_%d_decay",i);
2688         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
2689         h = fMCSecTrackHist1[2]->Projection(0);
2690         h->SetName(name);
2691         aFolderObj->Add(h);
2692
2693         snprintf(name,256,"mc_eta_rec_sec_pid_%d_decay",i);
2694         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
2695         h = fMCSecTrackHist1[2]->Projection(1);
2696         h->SetName(name);
2697         aFolderObj->Add(h);
2698
2699         snprintf(name,256,"mc_mother_rec_sec_pid_%d_decay",i);
2700         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
2701         h = fMCSecTrackHist1[2]->Projection(4);
2702         h->SetName(name);
2703         aFolderObj->Add(h);
2704
2705       } else if (j == 5) {
2706        // conversion
2707
2708         snprintf(name,256,"mc_pt_rec_sec_pid_%d_conv",i);
2709         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
2710         h = fMCSecTrackHist1[2]->Projection(0);
2711         h->SetName(name);
2712         aFolderObj->Add(h);
2713
2714         snprintf(name,256,"mc_eta_rec_sec_pid_%d_conv",i);
2715         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
2716         h = fMCSecTrackHist1[2]->Projection(1);
2717         h->SetName(name);
2718         aFolderObj->Add(h);
2719
2720         snprintf(name,256,"mc_mother_rec_sec_pid_%d_conv",i);
2721         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
2722         h = fMCSecTrackHist1[2]->Projection(4);
2723         h->SetName(name);
2724         aFolderObj->Add(h);
2725
2726      } else if (j == 13) {
2727        // mat
2728
2729         snprintf(name,256,"mc_pt_rec_sec_pid_%d_mat",i);
2730         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
2731         h = fMCSecTrackHist1[2]->Projection(0);
2732         h->SetName(name);
2733         aFolderObj->Add(h);
2734
2735         snprintf(name,256,"mc_eta_rec_sec_pid_%d_mat",i);
2736         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
2737         h = fMCSecTrackHist1[2]->Projection(1);
2738         h->SetName(name);
2739         aFolderObj->Add(h);
2740
2741         snprintf(name,256,"mc_eta_mother_rec_sec_pid_%d_mat",i);
2742         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
2743         h = fMCSecTrackHist1[2]->Projection(4,1);
2744         h->SetName(name);
2745         aFolderObj->Add(h);
2746
2747         snprintf(name,256,"mc_mother_rec_sec_pid_%d_mat",i);
2748         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
2749         h = fMCSecTrackHist1[2]->Projection(4);
2750         h->SetName(name);
2751         aFolderObj->Add(h);
2752
2753         snprintf(name,256,"mc_pt_mother_rec_sec_pid_%d_mat",i);
2754         fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
2755         h = fMCSecTrackHist1[2]->Projection(4,0);
2756         h->SetName(name);
2757         aFolderObj->Add(h);
2758
2759      } else {
2760        continue;
2761      }
2762    }
2763
2764   }
2765   } // end fHistogramOn
2766
2767   //
2768   //  resolution histograms
2769   //  only for reconstructed tracks
2770   //
2771
2772   TH2F *h2F=0;
2773   TCanvas * c = new TCanvas("resol","resol");
2774   c->cd();
2775
2776   //
2777   fRecMCTrackHist1->GetAxis(1)->SetRangeUser(-0.8,0.79);
2778
2779   h2F = (TH2F*)fRecMCTrackHist1->Projection(2,0);
2780   h = AlidNdPtHelper::MakeResol(h2F,1,0,kTRUE,10);
2781   h->SetXTitle("p_{tmc} (GeV/c)");
2782   h->SetYTitle("(p_{t}-p_{tmc})/p_{tmc} resolution");
2783   h->Draw();
2784   h->SetName("pt_resolution_vs_mcpt");
2785   aFolderObj->Add(h);
2786
2787   h2F = (TH2F*)fRecMCTrackHist1->Projection(2,0);
2788   h = AlidNdPtHelper::MakeResol(h2F,1,1,kTRUE,10);
2789   h->SetXTitle("p_{tmc} (GeV/c)");
2790   h->SetYTitle("(p_{t}-p_{tmc})/p_{tmc} mean");
2791   h->Draw();
2792   h->SetName("dpt_mean_vs_mcpt");
2793   aFolderObj->Add(h);
2794
2795   //
2796   h2F = (TH2F*)fRecMCTrackHist1->Projection(3,0);
2797   h = AlidNdPtHelper::MakeResol(h2F,1,0,kTRUE,10);
2798   h->SetXTitle("p_{tmc} (GeV/c)");
2799   h->SetYTitle("(#eta-#eta_{mc}) resolution");
2800   h->Draw();
2801   h->SetName("eta_resolution_vs_mcpt");
2802   aFolderObj->Add(h);
2803
2804   h2F = (TH2F*)fRecMCTrackHist1->Projection(3,0);
2805   h = AlidNdPtHelper::MakeResol(h2F,1,1,kTRUE,10);
2806   h->SetXTitle("p_{tmc} (GeV/c)");
2807   h->SetYTitle("(#eta-mc#eta) mean");
2808   h->Draw();
2809   h->SetName("deta_mean_vs_mcpt");
2810   aFolderObj->Add(h);
2811
2812   //
2813   fRecMCTrackHist1->GetAxis(1)->SetRange(1,fRecMCTrackHist1->GetAxis(1)->GetNbins());
2814
2815   h2F = (TH2F*)fRecMCTrackHist1->Projection(2,1);
2816   h = AlidNdPtHelper::MakeResol(h2F,1,0,kTRUE,10);
2817   h->SetXTitle("#eta_{mc}");
2818   h->SetYTitle("(p_{t}-p_{tmc})/p_{tmc} resolution");
2819   h->Draw();
2820   h->SetName("pt_resolution_vs_mceta");
2821   aFolderObj->Add(h);
2822
2823   h2F = (TH2F*)fRecMCTrackHist1->Projection(2,1);
2824   h = AlidNdPtHelper::MakeResol(h2F,1,1,kTRUE,10);
2825   h->SetXTitle("#eta_{mc}");
2826   h->SetYTitle("(p_{t}-p_{tmc})/p_{tmc} mean");
2827   h->Draw();
2828   h->SetName("dpt_mean_vs_mceta");
2829   aFolderObj->Add(h);
2830
2831   //
2832   h2F = (TH2F*)fRecMCTrackHist1->Projection(3,1);
2833   h = AlidNdPtHelper::MakeResol(h2F,1,0,kTRUE,10);
2834   h->SetXTitle("#eta_{mc}");
2835   h->SetYTitle("(#eta-#eta_{mc}) resolution");
2836   h->Draw();
2837   h->SetName("eta_resolution_vs_mceta");
2838   aFolderObj->Add(h);
2839
2840   h2F = (TH2F*)fRecMCTrackHist1->Projection(3,1);
2841   h = AlidNdPtHelper::MakeResol(h2F,1,1,kTRUE,10);
2842   h->SetXTitle("#eta_{mc}");
2843   h->SetYTitle("(#eta-mc#eta) mean");
2844   h->Draw();
2845   h->SetName("deta_mean_vs_mceta");
2846   aFolderObj->Add(h);
2847
2848   fRecMCTrackHist1->GetAxis(0)->SetRange(1,fRecMCTrackHist1->GetAxis(0)->GetNbins());
2849
2850   } // end use MC info
2851
2852   // export objects to analysis folder
2853   fAnalysisFolder = ExportToFolder(aFolderObj);
2854
2855   // delete only TObjArray
2856   if(aFolderObj) delete aFolderObj;
2857 }
2858
2859 //_____________________________________________________________________________
2860 TFolder* AlidNdPtAnalysis::ExportToFolder(TObjArray * const array)
2861 {
2862   // recreate folder avery time and export objects to new one
2863   //
2864   AlidNdPtAnalysis * comp=this;
2865   TFolder *folder = comp->GetAnalysisFolder();
2866
2867   TString name, title;
2868   TFolder *newFolder = 0;
2869   Int_t i = 0;
2870   Int_t size = array->GetSize();
2871
2872   if(folder) {
2873      // get name and title from old folder
2874      name = folder->GetName();
2875      title = folder->GetTitle();
2876
2877          // delete old one
2878      delete folder;
2879
2880          // create new one
2881      newFolder = CreateFolder(name.Data(),title.Data());
2882      newFolder->SetOwner();
2883
2884          // add objects to folder
2885      while(i < size) {
2886            newFolder->Add(array->At(i));
2887            i++;
2888          }
2889   }
2890
2891 return newFolder;
2892 }
2893
2894 //_____________________________________________________________________________
2895 TFolder* AlidNdPtAnalysis::CreateFolder(TString name,TString title) {
2896 // create folder for analysed histograms
2897 //
2898 TFolder *folder = 0;
2899   folder = new TFolder(name.Data(),title.Data());
2900
2901   return folder;
2902 }