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