X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PWGJE%2FEMCALJetTasks%2FUserTasks%2FAliAnalysisTaskCLQA.cxx;h=9026291eeebb268551a79b635c68f36173e061ae;hb=7548348774b23d2dbcc574310219fb3f3a16fc29;hp=e2e362ddf9806d012d2cc93da45032628b4ac934;hpb=5180389c9a53c9072c8e2e5d54cc5c3cc44bb157;p=u%2Fmrichter%2FAliRoot.git

diff --git a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskCLQA.cxx b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskCLQA.cxx
index e2e362ddf98..9026291eeeb 100644
--- a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskCLQA.cxx
+++ b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskCLQA.cxx
@@ -1,9 +1,11 @@
-// $Id: AliAnalysisTaskCLQA.cxx 60694 2013-02-04 15:35:56Z morsch $
+// $Id: $
 //
 // Constantin's Task
 //
 // Author: C.Loizides
 
+#include <complex>
+
 #include <TChain.h>
 #include <TClonesArray.h>
 #include <TDirectory.h>
@@ -15,11 +17,12 @@
 #include <TLorentzVector.h>
 #include <TNtuple.h>
 #include <TNtupleD.h>
+#include <TProfile.h>
 #include <TTree.h>
 
+#include "AliESDMuonTrack.h"
 #include "AliAODEvent.h"
 #include "AliAnalysisManager.h"
-#include "AliAnalysisTaskCLQA.h"
 #include "AliAnalysisUtils.h"
 #include "AliCentrality.h"
 #include "AliEMCALGeoParams.h"
@@ -35,27 +38,38 @@
 #include "AliVEventHandler.h"
 #include "AliVParticle.h"
 #include "AliVTrack.h"
+#include "AliAnalysisTaskCLQA.h"
 
 ClassImp(AliAnalysisTaskCLQA)
 
 //________________________________________________________________________
 AliAnalysisTaskCLQA::AliAnalysisTaskCLQA() : 
-  AliAnalysisTaskEmcalJet("AliAnalysisTaskCLQA", kTRUE),
-  fDoCumulants(0), 
-  fCumPtMin(0.3), fCumPtMax(5.0), fCumEtaMin(-1.0), fCumEtaMax(1.0), fCumMmin(15),
-  fNtupCum(0), fNtupCumInfo(0)
+  AliAnalysisTaskEmcal("AliAnalysisTaskCLQA", kTRUE),
+  fDo2013VertexCut(1),
+  fDoTracking(1), fDoMuonTracking(0), fDoCumulants(0), fDoCumNtuple(0),
+  fCumPtMin(0.3), fCumPtMax(5.0), fCumEtaMin(-1.0), fCumEtaMax(1.0), fCumMmin(15), 
+  fCumMbins(250), fCentCL1In(0), fCentV0AIn(0),
+  fNtupCum(0), fNtupCumInfo(0), fNtupZdcInfo(0)
 {
   // Default constructor.
+
+  for (Int_t i=0;i<1000;++i)
+    fHists[i] = 0;
 }
 
 //________________________________________________________________________
 AliAnalysisTaskCLQA::AliAnalysisTaskCLQA(const char *name) : 
-  AliAnalysisTaskEmcalJet(name, kTRUE),
-  fDoCumulants(0), 
+  AliAnalysisTaskEmcal(name, kTRUE),
+  fDo2013VertexCut(1),
+  fDoTracking(1), fDoMuonTracking(0), fDoCumulants(0), fDoCumNtuple(0),
   fCumPtMin(0.3), fCumPtMax(5.0), fCumEtaMin(-1.0), fCumEtaMax(1.0), fCumMmin(15),
-  fNtupCum(0), fNtupCumInfo(0)
+  fCumMbins(250), fCentCL1In(0), fCentV0AIn(0),
+  fNtupCum(0), fNtupCumInfo(0), fNtupZdcInfo(0)
 {
   // Standard constructor.
+
+  for (Int_t i=0;i<1000;++i)
+    fHists[i] = 0;
 }
 
 //________________________________________________________________________
@@ -69,6 +83,132 @@ Bool_t AliAnalysisTaskCLQA::FillHistograms()
 {
   // Fill histograms.
 
+  AliVEvent *event        = InputEvent();
+  AliAnalysisManager *am  = AliAnalysisManager::GetAnalysisManager();
+
+  UInt_t trig = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
+  for (Int_t i=0;i<31;++i) {
+    if (trig & (1<<i))
+      fHists[0]->Fill(trig);
+  }
+
+  Double_t vz  = event->GetPrimaryVertex()->GetZ();
+  fHists[1]->Fill(vz);
+
+  Int_t  run = event->GetRunNumber();
+  Int_t  vzn = event->GetPrimaryVertex()->GetNContributors();
+  if ((vzn<1)&&(run>0))
+    return kTRUE;
+  fHists[2]->Fill(vz);
+
+  if (TMath::Abs(vz)>10)
+    return kTRUE;
+
+  if (fDo2013VertexCut) {
+    if ((run>=188356&&run<=188366) || (run>=195344&&run<=197388)) {
+      AliAnalysisUtils anau;
+      if (anau.IsFirstEventInChunk(event))
+	return kFALSE;
+      if (!anau.IsVertexSelected2013pA(event))
+	return kFALSE;
+    }
+  }
+
+  // accepted events
+  fHists[9]->Fill(1,run);
+
+  AliCentrality *cent = InputEvent()->GetCentrality();
+  Double_t v0acent = cent->GetCentralityPercentile("V0A");
+  fHists[10]->Fill(v0acent);
+  Double_t znacent = cent->GetCentralityPercentile("ZNA");
+  fHists[11]->Fill(znacent);
+
+  if (fDoTracking) {
+    const Int_t ntracks = fTracks->GetEntries();
+    if (fTracks) {
+      for (Int_t i =0; i<ntracks; ++i) {
+        AliVParticle *track = dynamic_cast<AliVParticle*>(fTracks->At(i));
+        if (!track)
+          continue;
+        if (track->Charge()==0)
+          continue;
+        Double_t phi = track->Phi();
+        Double_t eta = track->Eta();
+        Double_t pt  = track->Pt();
+        fHists[20]->Fill(phi,eta);
+        fHists[21]->Fill(phi,pt);
+        fHists[22]->Fill(eta,pt);
+        if (TMath::Abs(eta)<0.8) {
+          fHists[23]->Fill(pt,v0acent);
+          fHists[24]->Fill(pt,znacent);
+        }
+	AliPicoTrack *picot = dynamic_cast<AliPicoTrack*>(track);
+	if (picot) {
+	  Int_t ttype = picot->GetTrackType();
+	  fHists[25+ttype]->Fill(phi,pt);
+	  if (picot->IsEMCAL()) {
+	    Double_t dphi = TVector2::Phi_mpi_pi(phi-picot->GetTrackPhiOnEMCal());
+	    fHists[28]->Fill(dphi,pt);
+	  }
+	}
+      }
+    }
+  }
+
+  if (fDoMuonTracking) {
+    AliAODEvent *aod = dynamic_cast<AliAODEvent*>(InputEvent());
+    if (aod) {
+      for (Int_t iMu = 0; iMu<aod->GetNumberOfTracks(); ++iMu) {
+        AliAODTrack* muonTrack = aod->GetTrack(iMu);
+        if (!muonTrack)
+          continue;
+        if (!muonTrack->IsMuonTrack()) 
+          continue;
+        Double_t dThetaAbs = TMath::ATan(muonTrack->GetRAtAbsorberEnd()/505.)* TMath::RadToDeg();
+        if ((dThetaAbs<2.) || (dThetaAbs>10.)) 
+          continue;
+        Double_t dEta = muonTrack->Eta();
+        if ((dEta<-4.) || (dEta>-2.5)) 
+          continue;
+        if (0) {
+          if (muonTrack->GetMatchTrigger()<0.5) 
+            continue;
+        }
+        Double_t ptMu  = muonTrack->Pt();
+        Double_t etaMu = muonTrack->Eta();
+        Double_t phiMu = muonTrack->Phi();
+        fHists[50]->Fill(phiMu,etaMu);
+        fHists[51]->Fill(phiMu,ptMu);
+        fHists[52]->Fill(etaMu,ptMu);
+        fHists[53]->Fill(ptMu,v0acent);
+        fHists[54]->Fill(ptMu,znacent);
+      }
+    } else {
+      AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
+      if (esd) {
+        for (Int_t iMu = 0; iMu<esd->GetNumberOfMuonTracks(); ++iMu) {
+          AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iMu);
+          if (!muonTrack)
+            continue;
+          if (!muonTrack->ContainTrackerData()) 
+            continue;
+          Double_t thetaTrackAbsEnd = TMath::ATan(muonTrack->GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
+          if ((thetaTrackAbsEnd < 2.) || (thetaTrackAbsEnd > 10.)) 
+            continue;
+          Double_t eta = muonTrack->Eta();
+          if ((eta < -4.) || (eta > -2.5))
+            return kFALSE;
+          if (0) {
+            if (!muonTrack->ContainTriggerData()) 
+              continue;
+            if (muonTrack->GetMatchTrigger() < 0.5) 
+              continue;
+          }
+        }
+      }
+    }
+  }
+
   return kTRUE;
 }
 
@@ -77,7 +217,7 @@ Bool_t AliAnalysisTaskCLQA::RetrieveEventObjects()
 {
   // Retrieve event objects.
 
-  if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
+  if (!AliAnalysisTaskEmcal::RetrieveEventObjects())
     return kFALSE;
 
   return kTRUE;
@@ -104,23 +244,44 @@ void AliAnalysisTaskCLQA::RunCumulants(Double_t Mmin, Double_t ptmin, Double_t p
   if (!fTracks) 
     return;
 
+  Bool_t isMC = 0;
+  TString tname(fTracks->GetName());
+  if (tname.Contains("mc"))
+    isMC = 1;
+
   const Int_t ntracks = fTracks->GetEntries();
   Int_t Mall=0,M=0,Mall2=0;
   Double_t ptmaxall=0,ptsumall=0,pt2sumall=0,ptsumall2=0;
   Double_t tsa00=0,tsa10=0,tsa11=0;
   Double_t Q2r=0,Q2i=0;
+  Double_t Q3r=0,Q3i=0;
   Double_t Q4r=0,Q4i=0;
+  Double_t Q6r=0,Q6i=0;
   Double_t mpt=0,mpt2=0,ptmaxq=0;
   Double_t ts00=0,ts10=0,ts11=0;
+  Double_t v0ach=0, v0cch=0;
+  Double_t cl1ch=0;
+ 
   for (Int_t i =0; i<ntracks; ++i) {
-    AliVTrack *track = dynamic_cast<AliVTrack*>(fTracks->At(i));
+    AliVParticle *track = dynamic_cast<AliVParticle*>(fTracks->At(i));
+    if (!track)
+      continue;
+    if (track->Charge()==0)
+      continue;
     Double_t eta = track->Eta();
+    if ((eta<5.1)&&(eta>2.8))
+      ++v0ach;
+    else if ((eta>-3.7)&&(eta<-1.7))
+      ++v0cch;
+    if (TMath::Abs(eta)<1.4) {
+      ++cl1ch;
+    }
     if ((eta<etamin) || (eta>etamax))
       continue;
     Double_t pt = track->Pt();
     if (pt>ptmaxall)
       ptmaxall = pt;
-    if (pt>2) {
+    if (pt>1) {
       ptsumall2 += pt;
       ++Mall2;
     }
@@ -145,16 +306,73 @@ void AliAnalysisTaskCLQA::RunCumulants(Double_t Mmin, Double_t ptmin, Double_t p
     ts11 += py*py/pt;
     Q2r  += TMath::Cos(2*phi);
     Q2i  += TMath::Sin(2*phi);
+    Q3r  += TMath::Cos(3*phi);
+    Q3i  += TMath::Sin(3*phi);
     Q4r  += TMath::Cos(4*phi);
     Q4i  += TMath::Sin(4*phi);
+    Q6r  += TMath::Cos(6*phi);
+    Q6i  += TMath::Sin(6*phi);
   }
 
-  if (M<Mmin)
+  if (M<=1)
     return;
 
-  Double_t Q2abs = Q2r*Q2r+Q2i*Q2i;
-  Double_t Q4abs = Q4r*Q4r+Q4i*Q4i;
-  Double_t Q42re = Q4r*Q2r*Q2r-Q4r*Q2i*Q2i+2*Q4i*Q2r*Q2i;
+  Int_t pmult=0;
+  Double_t v2g=0;
+  Double_t v3g=0;
+  for (Int_t i=0; i<ntracks; ++i) {
+    AliVParticle *track1 = dynamic_cast<AliVParticle*>(fTracks->At(i));
+    if (!track1)
+      continue;
+    if (track1->Charge()==0)
+      continue;
+    Double_t eta1 = track1->Eta();
+    if ((eta1<etamin) || (eta1>etamax))
+      continue;
+    Double_t pt1 = track1->Pt();
+    if ((pt1<ptmin) || (pt1>ptmax))
+      continue;
+    Double_t phi1 = track1->Phi();
+    for (Int_t j = i+1; j<ntracks; ++j) {
+      AliVParticle *track2 = dynamic_cast<AliVParticle*>(fTracks->At(j));
+      if (!track2)
+	continue;
+      if (track2->Charge()==0)
+	continue;
+      Double_t eta2 = track2->Eta();
+      if ((eta2<etamin) || (eta2>etamax))
+	continue;
+      Double_t pt2 = track2->Pt();
+      if ((pt2<ptmin) || (pt2>ptmax))
+	continue;
+      Double_t deta=TMath::Abs(eta2-eta1);
+      if(deta<1)
+	continue;
+      Double_t dphi=TVector2::Phi_0_2pi(phi1-track2->Phi());
+      pmult++;
+      v2g+=TMath::Cos(2*dphi);
+      v3g+=TMath::Cos(3*dphi);
+    }      
+  }
+
+  if (pmult>0) {
+    v2g/=pmult;
+    v3g/=pmult;
+  }
+
+  std::complex<double> q2(Q2r,Q2i);
+  std::complex<double> q3(Q3r,Q3i);
+  std::complex<double> q4(Q4r,Q4i);
+  std::complex<double> q6(Q6r,Q6i);
+  Double_t Q22   = std::abs(q2)*std::abs(q2);
+  Double_t Q32   = std::abs(q3)*std::abs(q3);
+  Double_t Q42   = std::abs(q4)*std::abs(q4);
+  Double_t Q62   = std::abs(q6)*std::abs(q6);
+  Double_t Q32re = std::real(q6*std::conj(q3)*std::conj(q3));
+  Double_t Q42re = std::real(q4*std::conj(q2)*std::conj(q2));
+  Double_t Q6are = std::real(q4*q2*std::conj(q2)*std::conj(q2)*std::conj(q2));
+  Double_t Q6bre = std::real(q6*std::conj(q2)*std::conj(q2)*std::conj(q2));
+  Double_t Q6cre = std::real(q6*std::conj(q4)*std::conj(q2));
 
   Double_t tsall = -1;
   Double_t tsax = (tsa00+tsa11)*(tsa00+tsa11)-4*(tsa00*tsa11-tsa10*tsa10);
@@ -172,6 +390,21 @@ void AliAnalysisTaskCLQA::RunCumulants(Double_t Mmin, Double_t ptmin, Double_t p
     ts = 2*l2/(l1+l2);
   }
 
+  if (isMC) {
+    fHists[106]->Fill(cl1ch);
+    fHists[107]->Fill(v0ach);
+    fHists[108]->Fill(v0cch);
+    AliCentrality *cent = InputEvent()->GetCentrality();
+    if (fCentCL1In) {
+      cent->SetCentralityCL1(100*fCentCL1In->GetBinContent(fCentCL1In->FindBin(cl1ch)));
+      cent->SetQuality(0);
+    }
+    if (fCentV0AIn) {
+      cent->SetCentralityV0A(100*fCentV0AIn->GetBinContent(fCentV0AIn->FindBin(v0ach)));
+      cent->SetQuality(0);
+    }
+  }
+
   AliAnalysisUtils anau;
   AliVEvent *event        = InputEvent();
   AliAnalysisManager *am  = AliAnalysisManager::GetAnalysisManager();
@@ -192,18 +425,29 @@ void AliAnalysisTaskCLQA::RunCumulants(Double_t Mmin, Double_t ptmin, Double_t p
   fNtupCumInfo->fPtMaxall = ptmaxall;
   fNtupCumInfo->fMPtall   = ptsumall/Mall;
   fNtupCumInfo->fMPt2all  = pt2sumall/Mall;
-  fNtupCumInfo->fMPtall2  = ptsumall2/Mall2;
+  if (Mall2>0)
+    fNtupCumInfo->fMPtall2  = ptsumall2/Mall2;
+  else
+    fNtupCumInfo->fMPtall2  = -1;
   fNtupCumInfo->fTSall    = tsall;
   fNtupCumInfo->fM        = M;
-  fNtupCumInfo->fQ2abs    = Q2abs;
-  fNtupCumInfo->fQ4abs    = Q4abs;
+  fNtupCumInfo->fQ2abs    = Q22;
+  fNtupCumInfo->fQ4abs    = Q42;
   fNtupCumInfo->fQ42re    = Q42re;
+  fNtupCumInfo->fCos2phi  = Q2r;
+  fNtupCumInfo->fSin2phi  = Q2i;
   fNtupCumInfo->fPtMax    = ptmaxq;
   fNtupCumInfo->fMPt      = mpt/M;
   fNtupCumInfo->fMPt2     = mpt2/M;
   fNtupCumInfo->fTS       = ts;
-  AliVVZERO *vzero = InputEvent()->GetVZEROData();
-  fNtupCumInfo->fMV0M     = vzero->GetMTotV0A()+vzero->GetMTotV0C();
+
+  if (isMC) {
+    fNtupCumInfo->fMV0M     = v0ach + v0cch;
+  } else {
+    AliVVZERO *vzero = InputEvent()->GetVZEROData();
+    fNtupCumInfo->fMV0M     = vzero->GetMTotV0A()+vzero->GetMTotV0C();
+  }
+
   AliCentrality *cent = InputEvent()->GetCentrality();
   fNtupCumInfo->fCl1      = cent->GetCentralityPercentile("CL1");
   fNtupCumInfo->fV0M      = cent->GetCentralityPercentile("V0M");
@@ -211,7 +455,110 @@ void AliAnalysisTaskCLQA::RunCumulants(Double_t Mmin, Double_t ptmin, Double_t p
   fNtupCumInfo->fV0A      = cent->GetCentralityPercentile("V0A");
   fNtupCumInfo->fV0AEq    = cent->GetCentralityPercentile("V0AEq");
   fNtupCumInfo->fZNA      = cent->GetCentralityPercentile("ZNA");
-  fNtupCum->Fill();
+
+  AliVZDC *vZDC = InputEvent()->GetZDCData();
+  const Double_t *znaTowers = vZDC->GetZNATowerEnergy(); 
+  fNtupZdcInfo->fZna0 = znaTowers[0];
+  fNtupZdcInfo->fZna1 = znaTowers[1];
+  fNtupZdcInfo->fZna2 = znaTowers[2];
+  fNtupZdcInfo->fZna3 = znaTowers[3];
+  fNtupZdcInfo->fZna4 = znaTowers[4];
+
+  if (fDoCumNtuple && (M>=Mmin)) {
+    fNtupCum->Fill();
+  }
+
+  fHists[109]->Fill(fNtupCumInfo->fCl1);
+  fHists[110]->Fill(fNtupCumInfo->fV0A);
+  fHists[111]->Fill(fNtupCumInfo->fZNA);
+
+  Bool_t fillCumHist = kTRUE;
+  if (fillCumHist) {
+    Int_t  run = InputEvent()->GetRunNumber();
+    if (fDo2013VertexCut) {
+      if ((run>=188356&&run<=188366) || (run>=195344&&run<=197388)) {
+	if (anau.IsFirstEventInChunk(event))
+	  fillCumHist = kFALSE;
+	if (!anau.IsVertexSelected2013pA(event))
+	  fillCumHist = kFALSE;
+      }
+    }
+    Double_t vz = InputEvent()->GetPrimaryVertex()->GetZ();
+    if (TMath::Abs(vz)>10)
+      fillCumHist = kFALSE;
+  }
+  if (fillCumHist) {
+    AliVVZERO *vzero = InputEvent()->GetVZEROData();
+    Double_t v0a = vzero->GetMTotV0A();
+    Double_t v0c = vzero->GetMTotV0C();
+    Double_t v0m = vzero->GetMTotV0A()+vzero->GetMTotV0C();
+    fHists[112]->Fill(Mall);
+    fHists[113]->Fill(M);
+    fHists[117]->Fill(v0a,M);
+    fHists[118]->Fill(v0c,M);
+    fHists[119]->Fill(v0a+v0c,M);
+    if (M>1) {
+      fHists[114]->Fill(M,(Q22-M)/M/(M-1));
+      fHists[120]->Fill(M,(Q32-M)/M/(M-1));
+      fHists[122]->Fill(M,v2g);
+      fHists[123]->Fill(M,v3g);
+    }
+    if (M>3) {
+      Double_t qc4tmp = (Q22*Q22+Q42-2*Q42re-4*(M-2)*Q22+2*M*(M-3));
+      fHists[115]->Fill(M,qc4tmp/M/(M-1)/(M-2)/(M-3));
+      qc4tmp = (Q32*Q32+Q62-2*Q32re-4*(M-2)*Q32+2*M*(M-3));
+      fHists[121]->Fill(M,qc4tmp/M/(M-1)/(M-2)/(M-3));
+    }
+    if (M>5) {
+      Double_t qc6tmp = Q22*Q22*Q22 + 9*Q42*Q22 - 6*Q6are 
+                      + 4*Q6bre - 12*Q6cre 
+	              + 18*(M-4)*Q42re + 4*Q62
+                      - 9*(M-4)*Q22*Q22 - 9*(M-4)*Q42
+                      + 18*(M-2)*(M-5)*Q22
+                      - 6*M*(M-4)*(M-5);
+      fHists[116]->Fill(M,qc6tmp/M/(M-1)/(M-2)/(M-3)/(M-4)/(M-5));
+    }
+  }
+
+  if ((isMC) || 
+      ((TMath::Abs(fNtupCumInfo->fVz)<10) && !fNtupCumInfo->fIsFEC && fNtupCumInfo->fIsVSel)) {
+    for (Int_t i =0; i<ntracks; ++i) {
+      AliVParticle *track1 = dynamic_cast<AliVParticle*>(fTracks->At(i));
+      if (!track1)
+        continue;
+      Double_t phi1 = track1->Phi();
+      Double_t eta1 = track1->Eta();
+      Double_t pt1  = track1->Pt();
+      ((TH3*)fHists[103])->Fill(pt1,eta1,fNtupCumInfo->fCl1);
+      ((TH3*)fHists[104])->Fill(pt1,eta1,fNtupCumInfo->fV0A);
+      ((TH3*)fHists[105])->Fill(pt1,eta1,fNtupCumInfo->fZNA);
+      if ((eta1<etamin) || (eta1>etamax))
+        continue;
+      if ((pt1<ptmin) || (pt1>ptmax))
+        continue;
+      for (Int_t j =0; j<ntracks; ++j) {
+        AliVParticle *track2 = dynamic_cast<AliVParticle*>(fTracks->At(j));
+        if (!track2)
+          continue;
+        Double_t eta2 = track2->Eta();
+        if ((eta2<etamin) || (eta2>etamax))
+          continue;
+        Double_t pt2 = track2->Pt();
+        if ((pt2<ptmin) || (pt2>ptmax))
+          continue;
+        Double_t phi2 = track2->Phi();
+        Double_t deta = eta1-eta2;
+        Double_t dphi = phi1-phi2;
+        while (dphi<-TMath::Pi())
+          dphi+=TMath::TwoPi();
+        while (dphi>3*TMath::Pi()/2)
+          dphi-=TMath::TwoPi();
+        ((TH3*)fHists[100])->Fill(dphi,deta,fNtupCumInfo->fCl1);
+        ((TH3*)fHists[101])->Fill(dphi,deta,fNtupCumInfo->fV0A);
+        ((TH3*)fHists[102])->Fill(dphi,deta,fNtupCumInfo->fZNA);
+      }
+    }
+  }
 }
 
 //________________________________________________________________________
@@ -231,24 +578,140 @@ void AliAnalysisTaskCLQA::UserCreateOutputObjects()
 {
   // Create histograms
 
-  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
+  AliAnalysisTaskEmcal::UserCreateOutputObjects();
+
+  fHists[0] = new TH1D("fTrigBits",";bit",32,-0.5,31.5);
+  fOutput->Add(fHists[0]);
+  fHists[1] = new TH1D("fVertexZ",";vertex z (cm)",51,-25.5,25.5);
+  fOutput->Add(fHists[1]);
+  fHists[2] = new TH1D("fVertexZnc",";vertex z (cm)",51,-25.5,25.5);
+  fOutput->Add(fHists[2]);
+  fHists[9] = new TProfile("fAccepted","",1,0.5,1.5);
+  fOutput->Add(fHists[9]);
+  fHists[10] = new TH1D("fV0ACent",";percentile",20,0,100);
+  fOutput->Add(fHists[10]);
+  fHists[11] = new TH1D("fZNACent",";percentile",20,0,100);
+  fOutput->Add(fHists[11]);
+
+  if (fDoTracking) {
+    fHists[20] = new TH2D("fPhiEtaTracks",";#phi;#eta",60,0,TMath::TwoPi(),20,-2,2);
+    fHists[20]->Sumw2();
+    fOutput->Add(fHists[20]);
+    fHists[21] = new TH2D("fPhiPtTracks",";#phi;p_{T} (GeV/c)",60,0,TMath::TwoPi(),40,0,20);
+    fHists[21]->Sumw2();
+    fOutput->Add(fHists[21]);
+    fHists[22] = new TH2D("fEtaPtTracks",";#eta;p_{T} (GeV/c)",20,-2,2,40,0,20);
+    fHists[22]->Sumw2();
+    fOutput->Add(fHists[22]);
+    fHists[23] = new TH2D("fPtV0ATracks",";#p_{T} (GeV/c);percentile",100,0,20,20,0,100);
+    fHists[23]->Sumw2();
+    fOutput->Add(fHists[23]);
+    fHists[24] = new TH2D("fPtZNATracks",";#p_{T} (GeV/c);percentile",100,0,20,20,0,100);
+    fHists[24]->Sumw2();
+    fOutput->Add(fHists[24]);
+    fHists[25] = new TH2D("fPhiPtTracks_type0",";#phi;p_{T} (GeV/c)",60,0,TMath::TwoPi(),40,0,20);
+    fHists[25]->Sumw2();
+    fOutput->Add(fHists[25]);
+    fHists[26] = new TH2D("fPhiPtTracks_type1",";#phi;p_{T} (GeV/c)",60,0,TMath::TwoPi(),40,0,20);
+    fHists[26]->Sumw2();
+    fOutput->Add(fHists[26]);
+    fHists[27] = new TH2D("fPhiPtTracks_type2",";#phi;p_{T} (GeV/c)",60,0,TMath::TwoPi(),40,0,20);
+    fHists[27]->Sumw2();
+    fOutput->Add(fHists[27]);
+    fHists[28] = new TH2D("fDPhiPtTracks",";#Delta#phi;p_{T} (GeV/c)",60,-TMath::Pi(),TMath::Pi(),40,0,20);
+    fHists[28]->Sumw2();
+    fOutput->Add(fHists[28]);
+  }
+
+  if (fDoMuonTracking) {
+    fHists[50] = new TH2D("fPhiEtaMuonTracks",";#phi;#eta",60,0,TMath::TwoPi(),15,-4,-2.5);
+    fHists[50]->Sumw2();
+    fOutput->Add(fHists[50]);
+    fHists[51] = new TH2D("fPhiPtMuonTracks",";#phi;p_{T} (GeV/c)",60,0,TMath::TwoPi(),200,0,20);
+    fHists[51]->Sumw2();
+    fOutput->Add(fHists[51]);
+    fHists[52] = new TH2D("fEtaPtMuonTracks",";#eta;p_{T} (GeV/c)",15,-4,-2.5,200,0,20);
+    fHists[52]->Sumw2();
+    fOutput->Add(fHists[52]);
+    fHists[53] = new TH2D("fPtV0AMuonTracks",";#p_{T} (GeV/c);percentile",100,0,10,20,0,100);
+    fHists[53]->Sumw2();
+    fOutput->Add(fHists[53]);
+    fHists[54] = new TH2D("fPtZNAMuonTracks",";#p_{T} (GeV/c);percentile",100,0,10,20,0,100);
+    fHists[54]->Sumw2();
+    fOutput->Add(fHists[54]);
+  }
 
   if (fDoCumulants) {
+    fHists[100] = new TH3D("fCumPhiEtaCl1",";#Delta#phi;#Delta#eta",32,-TMath::Pi()/2,3*TMath::Pi()/2,60,-3,3,10,0,100);
+    fOutput->Add(fHists[100]);
+    fHists[101] = new TH3D("fCumPhiEtaV0A",";#Delta#phi;#Delta#eta",32,-TMath::Pi()/2,3*TMath::Pi()/2,60,-3,3,10,0,100);
+    fOutput->Add(fHists[101]);
+    fHists[102] = new TH3D("fCumPhiEtaZNA",";#Delta#phi;#Delta#eta",32,-TMath::Pi()/2,3*TMath::Pi()/2,60,-3,3,10,0,100);
+    fOutput->Add(fHists[102]);
+    fHists[103] = new TH3D("fCumPtEtaCl1",";p_{T} (GeV/c);#eta",100,0,25,20,-2,2,10,0,100);
+    fOutput->Add(fHists[103]);
+    fHists[104] = new TH3D("fCumPtEtaV0A",";p_{T} (GeV/c);#eta",100,0,25,20,-2,2,10,0,100);
+    fOutput->Add(fHists[104]);
+    fHists[105] = new TH3D("fCumPtEtaZNA",";p_{T} (GeV/c);#eta",100,0,25,20,-2,2,10,0,100);
+    fOutput->Add(fHists[105]);
+    fHists[106] = new TH1D("fCumCL1MC",";#tracks",fCumMbins,0,fCumMbins);
+    fOutput->Add(fHists[106]);
+    fHists[107] = new TH1D("fCumV0AMC",";#tracks",fCumMbins,0,fCumMbins);
+    fOutput->Add(fHists[107]);
+    fHists[108] = new TH1D("fCumV0CMC",";#tracks",fCumMbins,0,fCumMbins);
+    fOutput->Add(fHists[108]);
+    fHists[109] = new TH1D("fCumCl1Cent",";percentile",10,0,100);
+    fOutput->Add(fHists[109]);
+    fHists[110] = new TH1D("fCumV0ACent",";percentile",10,0,100);
+    fOutput->Add(fHists[110]);
+    fHists[111] = new TH1D("fCumZNACent",";percentile",10,0,100);
+    fOutput->Add(fHists[111]);
+    fHists[112] = new TH1D("fCumMall",";mult",fCumMbins,0,fCumMbins);
+    fOutput->Add(fHists[112]);
+    fHists[113] = new TH1D("fCumM",";mult",fCumMbins,0,fCumMbins);
+    fOutput->Add(fHists[113]);
+    fHists[114] = new TProfile("fCumQC2",";qc2",fCumMbins,0,fCumMbins);
+    fOutput->Add(fHists[114]);
+    fHists[115] = new TProfile("fCumQC4",";qc4",fCumMbins,0,fCumMbins);
+    fOutput->Add(fHists[115]);
+    fHists[116] = new TProfile("fCumQC6",";qc6",fCumMbins,0,fCumMbins);
+    fOutput->Add(fHists[116]);
+    Int_t v0bins=1000;
+    if (fCumMbins>1000)
+      v0bins=25000;
+    fHists[117] = new TH2D("fCumV0ACentVsM",";v0a;M",v0bins,0,v0bins,fCumMbins,0,fCumMbins);
+    fOutput->Add(fHists[117]);
+    fHists[118] = new TH2D("fCumV0CCentVsM",";v0c;M",v0bins,0,v0bins,fCumMbins,0,fCumMbins);
+    fOutput->Add(fHists[118]);
+    fHists[119] = new TH2D("fCumV0MCentVsM",";v0m;M",v0bins,0,v0bins,fCumMbins,0,fCumMbins);
+    fOutput->Add(fHists[119]);
+    fHists[120] = new TProfile("fCum3QC2",";qc2",fCumMbins,0,fCumMbins);
+    fOutput->Add(fHists[120]);
+    fHists[121] = new TProfile("fCum3QC4",";qc4",fCumMbins,0,fCumMbins);
+    fOutput->Add(fHists[121]);
+    fHists[122] = new TProfile("fEtaGapV2",";M",fCumMbins,0,fCumMbins);
+    fOutput->Add(fHists[122]);
+    fHists[123] = new TProfile("fEtaGapV3",";M",fCumMbins,0,fCumMbins);
+    fOutput->Add(fHists[123]);
+
     fNtupCum = new TTree("NtupCum", "Ntuple for cumulant analysis");
     if (1) {
       fNtupCum->SetDirectory(0);
     } else {
       TFile *f = OpenFile(1); 
       if (f) {
-        f->SetCompressionLevel(2);
-        fNtupCum->SetDirectory(f);
-        fNtupCum->SetAutoFlush(-4*1024*1024);
-        fNtupCum->SetAutoSave(0);
+	f->SetCompressionLevel(2);
+	fNtupCum->SetDirectory(f);
+	fNtupCum->SetAutoFlush(-4*1024*1024);
+	fNtupCum->SetAutoSave(0);
       }
     }
     fNtupCumInfo = new AliNtupCumInfo;
     fNtupCum->Branch("cumulants", &fNtupCumInfo, 32*1024, 99);
-    fOutput->Add(fNtupCum);
+    fNtupZdcInfo = new AliNtupZdcInfo;
+    fNtupCum->Branch("zdc", &fNtupZdcInfo, 32*1024, 99);
+    if (fDoCumNtuple)
+      fOutput->Add(fNtupCum);
   }
 
   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram