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