PWG0/dNdPt/AlidNdPtAnalysisPbPb.cxx

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