PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskCLQA.cxx

   1 // $Id$
   2 //
   3 // Constantin's Task
   4 //
   5 // Author: C.Loizides
   6
   7 #include <complex>
   8
   9 #include <TChain.h>
  10 #include <TClonesArray.h>
  11 #include <TDirectory.h>
  12 #include <TFile.h>
  13 #include <TH1F.h>
  14 #include <TH2F.h>
  15 #include <TH3F.h>
  16 #include <TList.h>
  17 #include <TLorentzVector.h>
  18 #include <TNtuple.h>
  19 #include <TNtupleD.h>
  20 #include <TProfile.h>
  21 #include <TTree.h>
  22
  23 #include "AliESDMuonTrack.h"
  24 #include "AliAODEvent.h"
  25 #include "AliAnalysisManager.h"
  26 #include "AliAnalysisUtils.h"
  27 #include "AliCentrality.h"
  28 #include "AliEMCALGeoParams.h"
  29 #include "AliEMCALGeometry.h"
  30 #include "AliESDEvent.h"
  31 #include "AliEmcalJet.h"
  32 #include "AliExternalTrackParam.h"
  33 #include "AliInputEventHandler.h"
  34 #include "AliLog.h"
  35 #include "AliPicoTrack.h"
  36 #include "AliTrackerBase.h"
  37 #include "AliVCluster.h"
  38 #include "AliVEventHandler.h"
  39 #include "AliVParticle.h"
  40 #include "AliVTrack.h"
  41 #include "AliAnalysisTaskCLQA.h"
  42
  43 ClassImp(AliAnalysisTaskCLQA)
  44
  45 //________________________________________________________________________
  46 AliAnalysisTaskCLQA::AliAnalysisTaskCLQA() :
  47   AliAnalysisTaskEmcal("AliAnalysisTaskCLQA", kTRUE),
  48   fDoTracking(1), fDoMuonTracking(0), fDoCumulants(0), fDoCumNtuple(0),
  49   fCumPtMin(0.3), fCumPtMax(5.0), fCumEtaMin(-1.0), fCumEtaMax(1.0), fCumMmin(15),
  50   fCentCL1In(0), fCentV0AIn(0),
  51   fNtupCum(0), fNtupCumInfo(0), fNtupZdcInfo(0)
  52 {
  53   // Default constructor.
  54
  55   for (Int_t i=0;i<1000;++i)
  56     fHists[i] = 0;
  57 }
  58
  59 //________________________________________________________________________
  60 AliAnalysisTaskCLQA::AliAnalysisTaskCLQA(const char *name) :
  61   AliAnalysisTaskEmcal(name, kTRUE),
  62   fDoTracking(1), fDoMuonTracking(0), fDoCumulants(0), fDoCumNtuple(0),
  63   fCumPtMin(0.3), fCumPtMax(5.0), fCumEtaMin(-1.0), fCumEtaMax(1.0), fCumMmin(15),
  64   fCentCL1In(0), fCentV0AIn(0),
  65   fNtupCum(0), fNtupCumInfo(0), fNtupZdcInfo(0)
  66 {
  67   // Standard constructor.
  68
  69   for (Int_t i=0;i<1000;++i)
  70     fHists[i] = 0;
  71 }
  72
  73 //________________________________________________________________________
  74 AliAnalysisTaskCLQA::~AliAnalysisTaskCLQA()
  75 {
  76   // Destructor
  77 }
  78
  79 //________________________________________________________________________
  80 Bool_t AliAnalysisTaskCLQA::FillHistograms()
  81 {
  82   // Fill histograms.
  83
  84   AliVEvent *event        = InputEvent();
  85   AliAnalysisManager *am  = AliAnalysisManager::GetAnalysisManager();
  86
  87   UInt_t trig = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
  88   for (Int_t i=0;i<31;++i) {
  89     if (trig & (1<<i))
  90       fHists[0]->Fill(trig);
  91   }
  92
  93   Double_t vz  = event->GetPrimaryVertex()->GetZ();
  94   fHists[1]->Fill(vz);
  95
  96   Int_t  run = event->GetRunNumber();
  97   Int_t  vzn = event->GetPrimaryVertex()->GetNContributors();
  98   if ((vzn<1)&&(run>0))
  99     return kTRUE;
 100   fHists[2]->Fill(vz);
 101
 102   if (TMath::Abs(vz)>10)
 103     return kTRUE;
 104
 105   if ((run>=188356&&run<=188366) || (run>=195344&&run<=197388)) {
 106     AliAnalysisUtils anau;
 107     if (anau.IsFirstEventInChunk(event))
 108       return kFALSE;
 109     if (!anau.IsVertexSelected2013pA(event))
 110       return kFALSE;
 111   }
 112
 113   // accepted events
 114   fHists[9]->Fill(1);
 115
 116   AliCentrality *cent = InputEvent()->GetCentrality();
 117   Double_t v0acent = cent->GetCentralityPercentile("V0A");
 118   fHists[10]->Fill(v0acent);
 119   Double_t znacent = cent->GetCentralityPercentile("ZNA");
 120   fHists[11]->Fill(znacent);
 121
 122   if (fDoTracking) {
 123     const Int_t ntracks = fTracks->GetEntries();
 124     if (fTracks) {
 125       for (Int_t i =0; i<ntracks; ++i) {
 126         AliVParticle *track = dynamic_cast<AliVParticle*>(fTracks->At(i));
 127         if (!track)
 128           continue;
 129         if (track->Charge()==0)
 130           continue;
 131         Double_t phi = track->Phi();
 132         Double_t eta = track->Eta();
 133         Double_t pt  = track->Pt();
 134         fHists[20]->Fill(phi,eta);
 135         fHists[21]->Fill(phi,pt);
 136         fHists[22]->Fill(eta,pt);
 137         if (TMath::Abs(eta)<0.8) {
 138           fHists[23]->Fill(pt,v0acent);
 139           fHists[24]->Fill(pt,znacent);
 140         }
 141       }
 142     }
 143   }
 144
 145   if (fDoMuonTracking) {
 146     AliAODEvent *aod = dynamic_cast<AliAODEvent*>(InputEvent());
 147     if (aod) {
 148       for (Int_t iMu = 0; iMu<aod->GetNumberOfTracks(); ++iMu) {
 149         AliAODTrack* muonTrack = aod->GetTrack(iMu);
 150         if (!muonTrack)
 151           continue;
 152         if (!muonTrack->IsMuonTrack())
 153           continue;
 154         Double_t dThetaAbs = TMath::ATan(muonTrack->GetRAtAbsorberEnd()/505.)* TMath::RadToDeg();
 155         if ((dThetaAbs<2.) || (dThetaAbs>10.))
 156           continue;
 157         Double_t dEta = muonTrack->Eta();
 158         if ((dEta<-4.) || (dEta>-2.5))
 159           continue;
 160         if (0) {
 161           if (muonTrack->GetMatchTrigger()<0.5)
 162             continue;
 163         }
 164         Double_t ptMu  = muonTrack->Pt();
 165         Double_t etaMu = muonTrack->Eta();
 166         Double_t phiMu = muonTrack->Phi();
 167         fHists[50]->Fill(phiMu,etaMu);
 168         fHists[51]->Fill(phiMu,ptMu);
 169         fHists[52]->Fill(etaMu,ptMu);
 170         fHists[53]->Fill(ptMu,v0acent);
 171         fHists[54]->Fill(ptMu,znacent);
 172       }
 173     } else {
 174       AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
 175       if (esd) {
 176         for (Int_t iMu = 0; iMu<esd->GetNumberOfMuonTracks(); ++iMu) {
 177           AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iMu);
 178           if (!muonTrack)
 179             continue;
 180           if (!muonTrack->ContainTrackerData())
 181             continue;
 182           Double_t thetaTrackAbsEnd = TMath::ATan(muonTrack->GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
 183           if ((thetaTrackAbsEnd < 2.) || (thetaTrackAbsEnd > 10.))
 184             continue;
 185           Double_t eta = muonTrack->Eta();
 186           if ((eta < -4.) || (eta > -2.5))
 187             return kFALSE;
 188           if (0) {
 189             if (!muonTrack->ContainTriggerData())
 190               continue;
 191             if (muonTrack->GetMatchTrigger() < 0.5)
 192               continue;
 193           }
 194
 195         }
 196       }
 197     }
 198   }
 199
 200   return kTRUE;
 201 }
 202
 203 //________________________________________________________________________
 204 Bool_t AliAnalysisTaskCLQA::RetrieveEventObjects()
 205 {
 206   // Retrieve event objects.
 207
 208   if (!AliAnalysisTaskEmcal::RetrieveEventObjects())
 209     return kFALSE;
 210
 211   return kTRUE;
 212 }
 213
 214 //________________________________________________________________________
 215 Bool_t AliAnalysisTaskCLQA::Run()
 216 {
 217   // Run various functions.
 218
 219   RunCumulants(fCumMmin,fCumPtMin,fCumPtMax,fCumEtaMin,fCumEtaMax);
 220
 221   return kTRUE;
 222 }
 223
 224 //________________________________________________________________________
 225 void AliAnalysisTaskCLQA::RunCumulants(Double_t Mmin, Double_t ptmin, Double_t ptmax, Double_t etamin, Double_t etamax)
 226 {
 227   // Run cumulant analysis.
 228
 229   if (!fDoCumulants)
 230     return;
 231
 232   if (!fTracks)
 233     return;
 234
 235   Bool_t isMC = 0;
 236   TString tname(fTracks->GetName());
 237   if (tname.Contains("mc"))
 238     isMC = 1;
 239
 240   const Int_t ntracks = fTracks->GetEntries();
 241   Int_t Mall=0,M=0,Mall2=0;
 242   Double_t ptmaxall=0,ptsumall=0,pt2sumall=0,ptsumall2=0;
 243   Double_t tsa00=0,tsa10=0,tsa11=0;
 244   Double_t Q2r=0,Q2i=0;
 245   Double_t Q4r=0,Q4i=0;
 246   Double_t Q6r=0,Q6i=0;
 247   Double_t mpt=0,mpt2=0,ptmaxq=0;
 248   Double_t ts00=0,ts10=0,ts11=0;
 249   Double_t v0ach=0, v0cch=0;
 250   Double_t cl1ch=0;
 251
 252   for (Int_t i =0; i<ntracks; ++i) {
 253     AliVParticle *track = dynamic_cast<AliVParticle*>(fTracks->At(i));
 254     if (!track)
 255       continue;
 256     if (track->Charge()==0)
 257       continue;
 258     Double_t eta = track->Eta();
 259     if ((eta<5.1)&&(eta>2.8))
 260       ++v0ach;
 261     else if ((eta>-3.7)&&(eta<-1.7))
 262       ++v0cch;
 263     if (TMath::Abs(eta)<1.4) {
 264       ++cl1ch;
 265     }
 266     if ((eta<etamin) || (eta>etamax))
 267       continue;
 268     Double_t pt = track->Pt();
 269     if (pt>ptmaxall)
 270       ptmaxall = pt;
 271     if (pt>1) {
 272       ptsumall2 += pt;
 273       ++Mall2;
 274     }
 275     ptsumall  +=pt;
 276     pt2sumall +=pt*pt;
 277     Double_t px = track->Px();
 278     Double_t py = track->Py();
 279     tsa00 += px*px/pt;
 280     tsa10 += px*py/pt;
 281     tsa11 += py*py/pt;
 282     ++Mall;
 283     if ((pt<ptmin) || (pt>ptmax))
 284       continue;
 285     if (pt>ptmaxq)
 286       ptmaxq = pt;
 287     Double_t phi = track->Phi();
 288     ++M;
 289     mpt  += pt;
 290     mpt2 += pt*pt;
 291     ts00 += px*px/pt;
 292     ts10 += px*py/pt;
 293     ts11 += py*py/pt;
 294     Q2r  += TMath::Cos(2*phi);
 295     Q2i  += TMath::Sin(2*phi);
 296     Q4r  += TMath::Cos(4*phi);
 297     Q4i  += TMath::Sin(4*phi);
 298     Q6r  += TMath::Cos(6*phi);
 299     Q6i  += TMath::Sin(6*phi);
 300   }
 301
 302   if (M<=1)
 303     return;
 304
 305   std::complex<double> q2(Q2r,Q2i);
 306   std::complex<double> q4(Q4r,Q4i);
 307   std::complex<double> q6(Q6r,Q6i);
 308   Double_t Q22   = std::abs(q2)*std::abs(q2);
 309   Double_t Q42   = std::abs(q4)*std::abs(q4);
 310   Double_t Q62   = std::abs(q6)*std::abs(q6);
 311   Double_t Q42re = std::real(q4*std::conj(q2)*std::conj(q2));
 312   Double_t Q6are = std::real(q4*q2*std::conj(q2)*std::conj(q2)*std::conj(q2));
 313   Double_t Q6bre = std::real(q6*std::conj(q2)*std::conj(q2)*std::conj(q2));
 314   Double_t Q6cre = std::real(q6*std::conj(q4)*std::conj(q2));
 315
 316   Double_t tsall = -1;
 317   Double_t tsax = (tsa00+tsa11)*(tsa00+tsa11)-4*(tsa00*tsa11-tsa10*tsa10);
 318   if (tsax>=0) {
 319     Double_t l1 = 0.5*(tsa00+tsa11+TMath::Sqrt(tsax))/ptsumall;
 320     Double_t l2 = 0.5*(tsa00+tsa11-TMath::Sqrt(tsax))/ptsumall;
 321     tsall = 2*l2/(l1+l2);
 322   }
 323
 324   Double_t ts = -1;
 325   Double_t tsx = (ts00+ts11)*(ts00+ts11)-4*(ts00*ts11-ts10*ts10);
 326   if (tsx>=0) {
 327     Double_t l1 = 0.5*(ts00+ts11+TMath::Sqrt(tsx))/ptsumall;
 328     Double_t l2 = 0.5*(ts00+ts11-TMath::Sqrt(tsx))/ptsumall;
 329     ts = 2*l2/(l1+l2);
 330   }
 331
 332   if (isMC) {
 333     fHists[106]->Fill(cl1ch);
 334     fHists[107]->Fill(v0ach);
 335     fHists[108]->Fill(v0cch);
 336     AliCentrality *cent = InputEvent()->GetCentrality();
 337     if (fCentCL1In) {
 338       cent->SetCentralityCL1(100*fCentCL1In->GetBinContent(fCentCL1In->FindBin(cl1ch)));
 339       cent->SetQuality(0);
 340     }
 341     if (fCentV0AIn) {
 342       cent->SetCentralityV0A(100*fCentV0AIn->GetBinContent(fCentV0AIn->FindBin(v0ach)));
 343       cent->SetQuality(0);
 344     }
 345   }
 346
 347   AliAnalysisUtils anau;
 348   AliVEvent *event        = InputEvent();
 349   AliAnalysisManager *am  = AliAnalysisManager::GetAnalysisManager();
 350
 351   fNtupCumInfo->fTrig     = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
 352   fNtupCumInfo->fRun      = event->GetRunNumber();
 353   fNtupCumInfo->fVz       = event->GetPrimaryVertex()->GetZ();
 354   fNtupCumInfo->fIsFEC    = anau.IsFirstEventInChunk(event);
 355   fNtupCumInfo->fIsVSel   = anau.IsVertexSelected2013pA(event);
 356   fNtupCumInfo->fIsP      = event->IsPileupFromSPD(3/*minContributors*/,
 357                                                    0.8/*minZdist*/,
 358                                                    3./*nSigmaZdist*/,
 359                                                    2./*nSigmaDiamXY*/,
 360                                                    5./*nSigmaDiamZ*/);
 361
 362   fNtupCumInfo->fMall     = Mall;
 363   fNtupCumInfo->fMall2    = Mall2;
 364   fNtupCumInfo->fPtMaxall = ptmaxall;
 365   fNtupCumInfo->fMPtall   = ptsumall/Mall;
 366   fNtupCumInfo->fMPt2all  = pt2sumall/Mall;
 367   if (Mall2>0)
 368     fNtupCumInfo->fMPtall2  = ptsumall2/Mall2;
 369   else
 370     fNtupCumInfo->fMPtall2  = -1;
 371   fNtupCumInfo->fTSall    = tsall;
 372   fNtupCumInfo->fM        = M;
 373   fNtupCumInfo->fQ2abs    = Q22;
 374   fNtupCumInfo->fQ4abs    = Q42;
 375   fNtupCumInfo->fQ42re    = Q42re;
 376   fNtupCumInfo->fCos2phi  = Q2r;
 377   fNtupCumInfo->fSin2phi  = Q2i;
 378   fNtupCumInfo->fPtMax    = ptmaxq;
 379   fNtupCumInfo->fMPt      = mpt/M;
 380   fNtupCumInfo->fMPt2     = mpt2/M;
 381   fNtupCumInfo->fTS       = ts;
 382
 383   if (isMC) {
 384     fNtupCumInfo->fMV0M     = v0ach + v0cch;
 385   } else {
 386     AliVVZERO *vzero = InputEvent()->GetVZEROData();
 387     fNtupCumInfo->fMV0M     = vzero->GetMTotV0A()+vzero->GetMTotV0C();
 388   }
 389
 390   AliCentrality *cent = InputEvent()->GetCentrality();
 391   fNtupCumInfo->fCl1      = cent->GetCentralityPercentile("CL1");
 392   fNtupCumInfo->fV0M      = cent->GetCentralityPercentile("V0M");
 393   fNtupCumInfo->fV0MEq    = cent->GetCentralityPercentile("V0MEq");
 394   fNtupCumInfo->fV0A      = cent->GetCentralityPercentile("V0A");
 395   fNtupCumInfo->fV0AEq    = cent->GetCentralityPercentile("V0AEq");
 396   fNtupCumInfo->fZNA      = cent->GetCentralityPercentile("ZNA");
 397
 398   AliVZDC *vZDC = InputEvent()->GetZDCData();
 399   const Double_t *znaTowers = vZDC->GetZNATowerEnergy();
 400   fNtupZdcInfo->fZna0 = znaTowers[0];
 401   fNtupZdcInfo->fZna1 = znaTowers[1];
 402   fNtupZdcInfo->fZna2 = znaTowers[2];
 403   fNtupZdcInfo->fZna3 = znaTowers[3];
 404   fNtupZdcInfo->fZna4 = znaTowers[4];
 405
 406   if (fDoCumNtuple && (M>=Mmin)) {
 407     fNtupCum->Fill();
 408   }
 409
 410   fHists[109]->Fill(fNtupCumInfo->fCl1);
 411   fHists[110]->Fill(fNtupCumInfo->fV0A);
 412   fHists[111]->Fill(fNtupCumInfo->fZNA);
 413
 414   Bool_t fillCumHist = kTRUE;
 415   if (fillCumHist) {
 416     Int_t  run = InputEvent()->GetRunNumber();
 417     if ((run>=188356&&run<=188366) || (run>=195344&&run<=197388)) {
 418       if (anau.IsFirstEventInChunk(event))
 419         fillCumHist = kFALSE;
 420       if (!anau.IsVertexSelected2013pA(event))
 421         fillCumHist = kFALSE;
 422     }
 423     Double_t vz = InputEvent()->GetPrimaryVertex()->GetZ();
 424     if (TMath::Abs(vz)>10)
 425       fillCumHist = kFALSE;
 426   }
 427   if (fillCumHist) {
 428     fHists[112]->Fill(Mall);
 429     fHists[113]->Fill(M);
 430     if (M>1)
 431       fHists[114]->Fill(M,(Q22-M)/(M-1));
 432     if (M>3) {
 433       Double_t qc4tmp = (Q22*Q22+Q42-2*Q42re-4*(M-2)*Q22+2*M*(M-3));
 434       fHists[115]->Fill(M,qc4tmp/(M-1)/(M-2)/(M-3));
 435     }
 436     if (M>5) {
 437       Double_t qc6tmp = Q22*Q22*Q22 + 9*Q42*Q22 - 6*Q6are
 438                       + 4*Q6bre - 12*Q6cre
 439                       + 18*(M-4)*Q42re + 4*Q62*
 440                       - 9*(M-4)*Q22*Q22 - 9*(M-4)*Q42
 441                       + 18*(M-2)*(M-5)*Q22
 442                       -6*(M-4)*(M-5);
 443       fHists[116]->Fill(M,qc6tmp/(M-1)/(M-2)/(M-3)/(M-4)/(M-5));
 444     }
 445   }
 446
 447   if ((isMC) ||
 448       ((TMath::Abs(fNtupCumInfo->fVz)<10) && !fNtupCumInfo->fIsFEC && fNtupCumInfo->fIsVSel)) {
 449     for (Int_t i =0; i<ntracks; ++i) {
 450       AliVParticle *track1 = dynamic_cast<AliVParticle*>(fTracks->At(i));
 451       if (!track1)
 452         continue;
 453       Double_t phi1 = track1->Phi();
 454       Double_t eta1 = track1->Eta();
 455       Double_t pt1  = track1->Pt();
 456       ((TH3*)fHists[103])->Fill(pt1,eta1,fNtupCumInfo->fCl1);
 457       ((TH3*)fHists[104])->Fill(pt1,eta1,fNtupCumInfo->fV0A);
 458       ((TH3*)fHists[105])->Fill(pt1,eta1,fNtupCumInfo->fZNA);
 459       if ((eta1<etamin) || (eta1>etamax))
 460         continue;
 461       if ((pt1<ptmin) || (pt1>ptmax))
 462         continue;
 463       for (Int_t j =0; j<ntracks; ++j) {
 464         AliVParticle *track2 = dynamic_cast<AliVParticle*>(fTracks->At(j));
 465         if (!track2)
 466           continue;
 467         Double_t eta2 = track2->Eta();
 468         if ((eta2<etamin) || (eta2>etamax))
 469           continue;
 470         Double_t pt2 = track2->Pt();
 471         if ((pt2<ptmin) || (pt2>ptmax))
 472           continue;
 473         Double_t phi2 = track2->Phi();
 474         Double_t deta = eta1-eta2;
 475         Double_t dphi = phi1-phi2;
 476         while (dphi<-TMath::Pi())
 477           dphi+=TMath::TwoPi();
 478         while (dphi>3*TMath::Pi()/2)
 479           dphi-=TMath::TwoPi();
 480         ((TH3*)fHists[100])->Fill(dphi,deta,fNtupCumInfo->fCl1);
 481         ((TH3*)fHists[101])->Fill(dphi,deta,fNtupCumInfo->fV0A);
 482         ((TH3*)fHists[102])->Fill(dphi,deta,fNtupCumInfo->fZNA);
 483       }
 484     }
 485   }
 486 }
 487
 488 //________________________________________________________________________
 489 void AliAnalysisTaskCLQA::SetCumParams(Double_t Mmin, Double_t ptmin, Double_t ptmax, Double_t etamin, Double_t etamax)
 490 {
 491   // Set parameters for cumulants.
 492
 493   fCumMmin   = Mmin;
 494   fCumPtMin  = ptmin;
 495   fCumPtMax  = ptmax;
 496   fCumEtaMin = etamin;
 497   fCumEtaMax = etamax;
 498 }
 499
 500 //________________________________________________________________________
 501 void AliAnalysisTaskCLQA::UserCreateOutputObjects()
 502 {
 503   // Create histograms
 504
 505   AliAnalysisTaskEmcal::UserCreateOutputObjects();
 506
 507   fHists[0] = new TH1D("fTrigBits",";bit",32,-0.5,31.5);
 508   fOutput->Add(fHists[0]);
 509   fHists[1] = new TH1D("fVertexZ",";vertex z (cm)",51,-25.5,25.5);
 510   fOutput->Add(fHists[1]);
 511   fHists[2] = new TH1D("fVertexZnc",";vertex z (cm)",51,-25.5,25.5);
 512   fOutput->Add(fHists[2]);
 513   fHists[9] = new TH1D("fAccepted","",1,0.5,1.5);
 514   fOutput->Add(fHists[9]);
 515   fHists[10] = new TH1D("fV0ACent",";percentile",20,0,100);
 516   fOutput->Add(fHists[10]);
 517   fHists[11] = new TH1D("fZNACent",";percentile",20,0,100);
 518   fOutput->Add(fHists[11]);
 519
 520   if (fDoTracking) {
 521     fHists[20] = new TH2D("fPhiEtaTracks",";#phi;#eta",60,0,TMath::TwoPi(),20,-2,2);
 522     fHists[20]->Sumw2();
 523     fOutput->Add(fHists[20]);
 524     fHists[21] = new TH2D("fPhiPtTracks",";#phi;p_{T} (GeV/c)",60,0,TMath::TwoPi(),40,0,20);
 525     fHists[21]->Sumw2();
 526     fOutput->Add(fHists[21]);
 527     fHists[22] = new TH2D("fEtaPtTracks",";#eta;p_{T} (GeV/c)",20,-2,2,40,0,20);
 528     fHists[22]->Sumw2();
 529     fOutput->Add(fHists[22]);
 530     fHists[23] = new TH2D("fPtV0ATracks",";#p_{T} (GeV/c);percentile",100,0,20,20,0,100);
 531     fHists[23]->Sumw2();
 532     fOutput->Add(fHists[23]);
 533     fHists[24] = new TH2D("fPtZNATracks",";#p_{T} (GeV/c);percentile",100,0,20,20,0,100);
 534     fHists[24]->Sumw2();
 535     fOutput->Add(fHists[24]);
 536   }
 537
 538   if (fDoMuonTracking) {
 539     fHists[50] = new TH2D("fPhiEtaMuonTracks",";#phi;#eta",60,0,TMath::TwoPi(),15,-4,-2.5);
 540     fHists[50]->Sumw2();
 541     fOutput->Add(fHists[50]);
 542     fHists[51] = new TH2D("fPhiPtMuonTracks",";#phi;p_{T} (GeV/c)",60,0,TMath::TwoPi(),200,0,20);
 543     fHists[51]->Sumw2();
 544     fOutput->Add(fHists[51]);
 545     fHists[52] = new TH2D("fEtaPtMuonTracks",";#eta;p_{T} (GeV/c)",15,-4,-2.5,200,0,20);
 546     fHists[52]->Sumw2();
 547     fOutput->Add(fHists[52]);
 548     fHists[53] = new TH2D("fPtV0AMuonTracks",";#p_{T} (GeV/c);percentile",100,0,10,20,0,100);
 549     fHists[53]->Sumw2();
 550     fOutput->Add(fHists[53]);
 551     fHists[54] = new TH2D("fPtZNAMuonTracks",";#p_{T} (GeV/c);percentile",100,0,10,20,0,100);
 552     fHists[54]->Sumw2();
 553     fOutput->Add(fHists[54]);
 554   }
 555
 556   if (fDoCumulants) {
 557     fHists[100] = new TH3D("fCumPhiEtaCl1",";#Delta#phi;#Delta#eta",32,-TMath::Pi()/2,3*TMath::Pi()/2,60,-3,3,10,0,100);
 558     fOutput->Add(fHists[100]);
 559     fHists[101] = new TH3D("fCumPhiEtaV0A",";#Delta#phi;#Delta#eta",32,-TMath::Pi()/2,3*TMath::Pi()/2,60,-3,3,10,0,100);
 560     fOutput->Add(fHists[101]);
 561     fHists[102] = new TH3D("fCumPhiEtaZNA",";#Delta#phi;#Delta#eta",32,-TMath::Pi()/2,3*TMath::Pi()/2,60,-3,3,10,0,100);
 562     fOutput->Add(fHists[102]);
 563     fHists[103] = new TH3D("fCumPtEtaCl1",";p_{T} (GeV/c);#eta",100,0,25,20,-2,2,10,0,100);
 564     fOutput->Add(fHists[103]);
 565     fHists[104] = new TH3D("fCumPtEtaV0A",";p_{T} (GeV/c);#eta",100,0,25,20,-2,2,10,0,100);
 566     fOutput->Add(fHists[104]);
 567     fHists[105] = new TH3D("fCumPtEtaZNA",";p_{T} (GeV/c);#eta",100,0,25,20,-2,2,10,0,100);
 568     fOutput->Add(fHists[105]);
 569     fHists[106] = new TH1D("fCumCL1MC",";#tracks",2500,0,2500);
 570     fOutput->Add(fHists[106]);
 571     fHists[107] = new TH1D("fCumV0AMC",";#tracks",2500,0,2500);
 572     fOutput->Add(fHists[107]);
 573     fHists[108] = new TH1D("fCumV0CMC",";#tracks",2500,0,2500);
 574     fOutput->Add(fHists[108]);
 575     fHists[109] = new TH1D("fCumCl1Cent",";percentile",10,0,100);
 576     fOutput->Add(fHists[109]);
 577     fHists[110] = new TH1D("fCumV0ACent",";percentile",10,0,100);
 578     fOutput->Add(fHists[110]);
 579     fHists[111] = new TH1D("fCumZNACent",";percentile",10,0,100);
 580     fOutput->Add(fHists[111]);
 581     fHists[112] = new TH1D("fCumMall",";mult",2500,0,2500);
 582     fOutput->Add(fHists[112]);
 583     fHists[113] = new TH1D("fCumM",";mult",2500,0,2500);
 584     fOutput->Add(fHists[113]);
 585     fHists[114] = new TProfile("fCumQC2",";qc2",2500,0,2500);
 586     fOutput->Add(fHists[114]);
 587     fHists[115] = new TProfile("fCumQC4",";qc4",2500,0,2500);
 588     fOutput->Add(fHists[115]);
 589     fHists[116] = new TProfile("fCumQC6",";qc6",2500,0,2500);
 590     fOutput->Add(fHists[116]);
 591
 592     fNtupCum = new TTree("NtupCum", "Ntuple for cumulant analysis");
 593     if (1) {
 594       fNtupCum->SetDirectory(0);
 595     } else {
 596       TFile *f = OpenFile(1);
 597       if (f) {
 598         f->SetCompressionLevel(2);
 599         fNtupCum->SetDirectory(f);
 600         fNtupCum->SetAutoFlush(-4*1024*1024);
 601         fNtupCum->SetAutoSave(0);
 602       }
 603     }
 604     fNtupCumInfo = new AliNtupCumInfo;
 605     fNtupCum->Branch("cumulants", &fNtupCumInfo, 32*1024, 99);
 606     fNtupZdcInfo = new AliNtupZdcInfo;
 607     fNtupCum->Branch("zdc", &fNtupZdcInfo, 32*1024, 99);
 608     if (fDoCumNtuple)
 609       fOutput->Add(fNtupCum);
 610   }
 611
 612   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
 613 }