PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskSAQA.cxx

   1 // $Id$
   2 //
   3 // General QA task.
   4 //
   5 // Author: S.Aiola
   6
   7 #include <TClonesArray.h>
   8 #include <TH1F.h>
   9 #include <TH2F.h>
  10 #include <TH3F.h>
  11 #include <THnSparse.h>
  12 #include <TList.h>
  13 #include <TLorentzVector.h>
  14
  15 #include "AliAnalysisManager.h"
  16 #include "AliCentrality.h"
  17 #include "AliVCluster.h"
  18 #include "AliVParticle.h"
  19 #include "AliVTrack.h"
  20 #include "AliEmcalJet.h"
  21 #include "AliVEventHandler.h"
  22 #include "AliAODEvent.h"
  23 #include "AliExternalTrackParam.h"
  24 #include "AliTrackerBase.h"
  25 #include "AliLog.h"
  26 #include "AliEMCALGeometry.h"
  27 #include "AliEMCALGeoParams.h"
  28 #include "AliPicoTrack.h"
  29 #include "AliVVZERO.h"
  30 #include "AliESDUtils.h"
  31
  32 #include "AliAnalysisTaskSAQA.h"
  33
  34 ClassImp(AliAnalysisTaskSAQA)
  35
  36 //________________________________________________________________________
  37 AliAnalysisTaskSAQA::AliAnalysisTaskSAQA() :
  38   AliAnalysisTaskEmcalJet("AliAnalysisTaskSAQA", kTRUE),
  39   fCellEnergyCut(0.05),
  40   fParticleLevel(kFALSE),
  41   fIsMC(kFALSE),
  42   fCentMethod2(""),
  43   fCentMethod3(""),
  44   fDoV0QA(0),
  45   fDoEPQA(0),
  46   fDoLeadingObjectPosition(0),
  47   fMaxCellsInCluster(50),
  48   fCent2(0),
  49   fCent3(0),
  50   fVZERO(0),
  51   fV0ATotMult(0),
  52   fV0CTotMult(0),
  53   fHistEventQA(0),
  54   fHistTrNegativeLabels(0),
  55   fHistTrZeroLabels(0),
  56   fHistTrPhiEtaPt(0),
  57   fHistTrPhiEtaZeroLab(0),
  58   fHistTrPtZeroLab(0),
  59   fHistTrEmcPhiEta(0),
  60   fHistTrEmcPt(0),
  61   fHistTrPhiEtaNonProp(0),
  62   fHistTrPtNonProp(0),
  63   fHistDeltaEtaPt(0),
  64   fHistDeltaPhiPt(0),
  65   fHistDeltaPtvsPt(0),
  66   fHistClusPhiEtaEnergy(0),
  67   fHistClusDeltaPhiEPEnergy(0),
  68   fHistNCellsEnergy(0),
  69   fHistFcrossEnergy(0),
  70   fHistClusTimeEnergy(0),
  71   fHistClusMCEnergyFraction(0),
  72   fHistCellsAbsIdEnergy(0),
  73   fHistJetsPhiEta(0),
  74   fHistJetsPtArea(0)
  75 {
  76   // Default constructor.
  77
  78   SetMakeGeneralHistograms(kTRUE);
  79 }
  80
  81 //________________________________________________________________________
  82 AliAnalysisTaskSAQA::AliAnalysisTaskSAQA(const char *name) :
  83   AliAnalysisTaskEmcalJet(name, kTRUE),
  84   fCellEnergyCut(0.05),
  85   fParticleLevel(kFALSE),
  86   fIsMC(kFALSE),
  87   fCentMethod2(""),
  88   fCentMethod3(""),
  89   fDoV0QA(0),
  90   fDoEPQA(0),
  91   fDoLeadingObjectPosition(0),
  92   fMaxCellsInCluster(50),
  93   fCent2(0),
  94   fCent3(0),
  95   fVZERO(0),
  96   fV0ATotMult(0),
  97   fV0CTotMult(0),
  98   fHistEventQA(0),
  99   fHistTrNegativeLabels(0),
 100   fHistTrZeroLabels(0),
 101   fHistTrPhiEtaPt(0),
 102   fHistTrPhiEtaZeroLab(0),
 103   fHistTrPtZeroLab(0),
 104   fHistTrEmcPhiEta(0),
 105   fHistTrEmcPt(0),
 106   fHistTrPhiEtaNonProp(0),
 107   fHistTrPtNonProp(0),
 108   fHistDeltaEtaPt(0),
 109   fHistDeltaPhiPt(0),
 110   fHistDeltaPtvsPt(0),
 111   fHistClusPhiEtaEnergy(0),
 112   fHistClusDeltaPhiEPEnergy(0),
 113   fHistNCellsEnergy(0),
 114   fHistFcrossEnergy(0),
 115   fHistClusTimeEnergy(0),
 116   fHistClusMCEnergyFraction(0),
 117   fHistCellsAbsIdEnergy(0),
 118   fHistJetsPhiEta(0),
 119   fHistJetsPtArea(0)
 120 {
 121   // Standard
 122
 123   SetMakeGeneralHistograms(kTRUE);
 124 }
 125
 126 //________________________________________________________________________
 127 AliAnalysisTaskSAQA::~AliAnalysisTaskSAQA()
 128 {
 129   // Destructor
 130 }
 131
 132 //________________________________________________________________________
 133 void AliAnalysisTaskSAQA::AllocateHistogramArrays()
 134 {
 135   fHistTrNegativeLabels = new TH1*[fNcentBins];
 136   fHistTrZeroLabels = new TH1*[fNcentBins];
 137   fHistTrPhiEtaZeroLab = new TH2*[fNcentBins];
 138   fHistTrPtZeroLab = new TH1*[fNcentBins];
 139   fHistTrEmcPhiEta = new TH2*[fNcentBins];
 140   fHistTrEmcPt = new TH1*[fNcentBins];
 141   fHistTrPhiEtaNonProp = new TH2*[fNcentBins];
 142   fHistTrPtNonProp = new TH1*[fNcentBins];
 143   fHistDeltaEtaPt = new TH2*[fNcentBins];
 144   fHistDeltaPhiPt = new TH2*[fNcentBins];
 145   fHistDeltaPtvsPt = new TH2*[fNcentBins];
 146
 147   fHistClusPhiEtaEnergy = new TH3*[fNcentBins];
 148   fHistClusDeltaPhiEPEnergy = new TH2*[fNcentBins];
 149   fHistNCellsEnergy = new TH2*[fNcentBins];
 150   fHistFcrossEnergy = new TH2*[fNcentBins];
 151   fHistClusTimeEnergy = new TH2*[fNcentBins];
 152   fHistClusMCEnergyFraction = new TH1*[fNcentBins];
 153
 154   fHistCellsAbsIdEnergy = new TH2*[fNcentBins];
 155
 156   fHistJetsPhiEta = new TH2*[fNcentBins];
 157   fHistJetsPtArea = new TH2*[fNcentBins];
 158
 159   fHistTrPhiEtaPt = new TH3**[fNcentBins];
 160
 161   for (Int_t i = 0; i < fNcentBins; i++) {
 162
 163     fHistTrPhiEtaPt[i] = new TH3*[4];
 164     for (Int_t j = 0; j < 4; j++) fHistTrPhiEtaPt[i][j] = 0;
 165
 166     fHistTrNegativeLabels[i] = 0;
 167     fHistTrZeroLabels[i] = 0;
 168     fHistTrPhiEtaZeroLab[i] = 0;
 169     fHistTrPtZeroLab[i] = 0;
 170     fHistTrEmcPhiEta[i] = 0;
 171     fHistTrEmcPt[i] = 0;
 172     fHistTrPhiEtaNonProp[i] = 0;
 173     fHistTrPtNonProp[i] = 0;
 174     fHistDeltaEtaPt[i] = 0;
 175     fHistDeltaPhiPt[i] = 0;
 176     fHistDeltaPtvsPt[i] = 0;
 177
 178     fHistClusPhiEtaEnergy[i] = 0;
 179     fHistClusDeltaPhiEPEnergy[i] = 0;
 180     fHistNCellsEnergy[i] = 0;
 181     fHistFcrossEnergy[i] = 0;
 182     fHistClusTimeEnergy[i] = 0;
 183     fHistClusMCEnergyFraction[i] = 0;
 184
 185     fHistCellsAbsIdEnergy[i] = 0;
 186
 187     fHistJetsPhiEta[i] = 0;
 188     fHistJetsPtArea[i] = 0;
 189   }
 190 }
 191
 192 //________________________________________________________________________
 193 void AliAnalysisTaskSAQA::UserCreateOutputObjects()
 194 {
 195   // Create histograms
 196
 197   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 198
 199   AllocateHistogramArrays();
 200
 201   TString histname;
 202
 203   if (fParticleCollArray.GetEntriesFast()>0) {
 204     if (!fParticleLevel && fIsMC) {
 205       for (Int_t i = 0; i < fNcentBins; i++) {
 206         histname = Form("fHistTrNegativeLabels_%d",i);
 207         fHistTrNegativeLabels[i] = new TH1F(histname,histname, 500, 0, 1);
 208         fHistTrNegativeLabels[i]->GetXaxis()->SetTitle("% of negative labels");
 209         fHistTrNegativeLabels[i]->GetYaxis()->SetTitle("counts");
 210         fOutput->Add(fHistTrNegativeLabels[i]);
 211
 212         histname = Form("fHistTrZeroLabels_%d",i);
 213         fHistTrZeroLabels[i] = new TH1F(histname,histname, 500, 0, 1);
 214         fHistTrZeroLabels[i]->GetXaxis()->SetTitle("% of negative labels");
 215         fHistTrZeroLabels[i]->GetYaxis()->SetTitle("counts");
 216         fOutput->Add(fHistTrZeroLabels[i]);
 217       }
 218     }
 219
 220     Int_t nlabels = 4;
 221     if (fParticleLevel)
 222       nlabels = 1;
 223
 224     for (Int_t i = 0; i < fNcentBins; i++) {
 225       for (Int_t j = 0; j < nlabels; j++) {
 226         histname = Form("fHistTrPhiEtaPt_%d_%d",i,j);
 227         fHistTrPhiEtaPt[i][j] = new TH3F(histname,histname, 100, -1, 1, 101, 0, TMath::Pi() * 2.02, fNbins, fMinBinPt, fMaxBinPt);
 228         fHistTrPhiEtaPt[i][j]->GetXaxis()->SetTitle("#eta");
 229         fHistTrPhiEtaPt[i][j]->GetYaxis()->SetTitle("#phi");
 230         fHistTrPhiEtaPt[i][j]->GetZaxis()->SetTitle("p_{T} (GeV/c)");
 231         fOutput->Add(fHistTrPhiEtaPt[i][j]);
 232       }
 233
 234       if (!fParticleLevel) {
 235         if (fIsMC) {
 236           histname = Form("fHistTrPhiEtaZeroLab_%d",i);
 237           fHistTrPhiEtaZeroLab[i] = new TH2F(histname,histname, 100, -1, 1, 101, 0, TMath::Pi() * 2.02);
 238           fHistTrPhiEtaZeroLab[i]->GetXaxis()->SetTitle("#eta");
 239           fHistTrPhiEtaZeroLab[i]->GetYaxis()->SetTitle("#phi");
 240           fHistTrPhiEtaZeroLab[i]->GetZaxis()->SetTitle("counts");
 241           fOutput->Add(fHistTrPhiEtaZeroLab[i]);
 242
 243           histname = Form("fHistTrPtZeroLab_%d",i);
 244           fHistTrPtZeroLab[i] = new TH1F(histname,histname, fNbins, fMinBinPt, fMaxBinPt);
 245           fHistTrPtZeroLab[i]->GetZaxis()->SetTitle("p_{T} (GeV/c)");
 246           fHistTrPtZeroLab[i]->GetYaxis()->SetTitle("counts");
 247           fOutput->Add(fHistTrPtZeroLab[i]);
 248         }
 249
 250         histname = Form("fHistTrEmcPhiEta_%d",i);
 251         fHistTrEmcPhiEta[i] = new TH2F(histname,histname, 100, -1, 1, 101, 0, TMath::Pi() * 2.02);
 252         fHistTrEmcPhiEta[i]->GetXaxis()->SetTitle("#eta");
 253         fHistTrEmcPhiEta[i]->GetYaxis()->SetTitle("#phi");
 254         fOutput->Add(fHistTrEmcPhiEta[i]);
 255
 256         histname = Form("fHistTrEmcPt_%d",i);
 257         fHistTrEmcPt[i] = new TH1F(histname,histname, fNbins, fMinBinPt, fMaxBinPt);
 258         fHistTrEmcPt[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 259         fHistTrEmcPt[i]->GetYaxis()->SetTitle("counts");
 260         fOutput->Add(fHistTrEmcPt[i]);
 261
 262         histname = Form("fHistTrPhiEtaNonProp_%d",i);
 263         fHistTrPhiEtaNonProp[i] = new TH2F(histname,histname, 100, -1, 1, 101, 0, TMath::Pi() * 2.02);
 264         fHistTrPhiEtaNonProp[i]->GetXaxis()->SetTitle("#eta");
 265         fHistTrPhiEtaNonProp[i]->GetYaxis()->SetTitle("#phi");
 266         fOutput->Add(fHistTrPhiEtaNonProp[i]);
 267
 268         histname = Form("fHistTrPtNonProp_%d",i);
 269         fHistTrPtNonProp[i] = new TH1F(histname,histname, fNbins, fMinBinPt, fMaxBinPt);
 270         fHistTrPtNonProp[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 271         fHistTrPtNonProp[i]->GetYaxis()->SetTitle("counts");
 272         fOutput->Add(fHistTrPtNonProp[i]);
 273
 274         histname = Form("fHistDeltaEtaPt_%d",i);
 275         fHistDeltaEtaPt[i] = new TH2F(histname,histname, fNbins, fMinBinPt, fMaxBinPt, 50, -0.5, 0.5);
 276         fHistDeltaEtaPt[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 277         fHistDeltaEtaPt[i]->GetYaxis()->SetTitle("#delta#eta");
 278         fOutput->Add(fHistDeltaEtaPt[i]);
 279
 280         histname = Form("fHistDeltaPhiPt_%d",i);
 281         fHistDeltaPhiPt[i] = new TH2F(histname,histname, fNbins, fMinBinPt, fMaxBinPt, 200, -2, 2);
 282         fHistDeltaPhiPt[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 283         fHistDeltaPhiPt[i]->GetYaxis()->SetTitle("#delta#phi");
 284         fOutput->Add(fHistDeltaPhiPt[i]);
 285
 286         histname = Form("fHistDeltaPtvsPt_%d",i);
 287         fHistDeltaPtvsPt[i] = new TH2F(histname,histname, fNbins, fMinBinPt, fMaxBinPt, fNbins, -fMaxBinPt/2, fMaxBinPt/2);
 288         fHistDeltaPtvsPt[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 289         fHistDeltaPtvsPt[i]->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
 290         fHistDeltaPtvsPt[i]->GetZaxis()->SetTitle("counts");
 291         fOutput->Add(fHistDeltaPtvsPt[i]);
 292       }
 293     }
 294   }
 295
 296   if (fClusterCollArray.GetEntriesFast()>0) {
 297     for (Int_t i = 0; i < fNcentBins; i++) {
 298       histname = "fHistClusPhiEtaEnergy_";
 299       histname += i;
 300       fHistClusPhiEtaEnergy[i] = new TH3F(histname, histname, 100, -1, 1, 101, 0, TMath::Pi() * 2.02, fNbins, fMinBinPt, fMaxBinPt);
 301       fHistClusPhiEtaEnergy[i]->GetXaxis()->SetTitle("#eta");
 302       fHistClusPhiEtaEnergy[i]->GetYaxis()->SetTitle("#phi");
 303       fHistClusPhiEtaEnergy[i]->GetZaxis()->SetTitle("E_{cluster} (GeV)");
 304       fOutput->Add(fHistClusPhiEtaEnergy[i]);
 305
 306       histname = "fHistClusDeltaPhiEPEnergy_";
 307       histname += i;
 308       fHistClusDeltaPhiEPEnergy[i] = new TH2F(histname, histname, fNbins, fMinBinPt, fMaxBinPt, 100, 0, TMath::Pi());
 309       fHistClusDeltaPhiEPEnergy[i]->GetXaxis()->SetTitle("E_{cluster} (GeV)");
 310       fHistClusDeltaPhiEPEnergy[i]->GetYaxis()->SetTitle("#phi_{cluster} - #psi_{RP}");
 311       fOutput->Add(fHistClusDeltaPhiEPEnergy[i]);
 312
 313       if (fIsEmbedded) {
 314         histname = "fHistClusMCEnergyFraction_";
 315         histname += i;
 316         fHistClusMCEnergyFraction[i] = new TH1F(histname, histname, fNbins, 0, 1.2);
 317         fHistClusMCEnergyFraction[i]->GetXaxis()->SetTitle("MC fraction");
 318         fHistClusMCEnergyFraction[i]->GetYaxis()->SetTitle("counts");
 319         fOutput->Add(fHistClusMCEnergyFraction[i]);
 320       }
 321
 322       histname = "fHistClusTimeEnergy_";
 323       histname += i;
 324       fHistClusTimeEnergy[i] = new TH2F(histname,histname, fNbins, fMinBinPt, fMaxBinPt, fNbins,  -1e-6, 1e-6);
 325       fHistClusTimeEnergy[i]->GetXaxis()->SetTitle("E_{cluster} (GeV)");
 326       fHistClusTimeEnergy[i]->GetYaxis()->SetTitle("Time");
 327       fOutput->Add(fHistClusTimeEnergy[i]);
 328
 329       Int_t nbins = fMaxCellsInCluster;
 330       while (nbins > fNbins) nbins /= 2;
 331       histname = "fHistNCellsEnergy_";
 332       histname += i;
 333       fHistNCellsEnergy[i] = new TH2F(histname,histname, fNbins, fMinBinPt, fMaxBinPt, nbins, -0.5, fMaxCellsInCluster-0.5);
 334       fHistNCellsEnergy[i]->GetXaxis()->SetTitle("E_{cluster} (GeV)");
 335       fHistNCellsEnergy[i]->GetYaxis()->SetTitle("N_{cells}");
 336       fOutput->Add(fHistNCellsEnergy[i]);
 337
 338       histname = "fHistFcrossEnergy_";
 339       histname += i;
 340       fHistFcrossEnergy[i] = new TH2F(histname,histname, fNbins, fMinBinPt, fMaxBinPt, 200, -3.5, 1.5);
 341       fHistFcrossEnergy[i]->GetXaxis()->SetTitle("E_{cluster} (GeV)");
 342       fHistFcrossEnergy[i]->GetYaxis()->SetTitle("F_{cross}");
 343       fOutput->Add(fHistFcrossEnergy[i]);
 344
 345       histname = "fHistCellsAbsIdEnergy_";
 346       histname += i;
 347       fHistCellsAbsIdEnergy[i] = new TH2F(histname,histname, 11600,0,11599,(Int_t)(fNbins / 2), fMinBinPt, fMaxBinPt / 2);
 348       fHistCellsAbsIdEnergy[i]->GetXaxis()->SetTitle("cell abs. Id");
 349       fHistCellsAbsIdEnergy[i]->GetYaxis()->SetTitle("E_{cluster} (GeV)");
 350       fHistCellsAbsIdEnergy[i]->GetZaxis()->SetTitle("counts");
 351       fOutput->Add(fHistCellsAbsIdEnergy[i]);
 352     }
 353   }
 354
 355   if (fJetCollArray.GetEntriesFast()>0) {
 356     for (Int_t i = 0; i < fNcentBins; i++) {
 357       histname = "fHistJetsPhiEta_";
 358       histname += i;
 359       fHistJetsPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 100, -1, 1, 101, 0, TMath::Pi() * 2.02);
 360       fHistJetsPhiEta[i]->GetXaxis()->SetTitle("#eta");
 361       fHistJetsPhiEta[i]->GetYaxis()->SetTitle("#phi");
 362       fHistJetsPhiEta[i]->GetZaxis()->SetTitle("counts");
 363       fOutput->Add(fHistJetsPhiEta[i]);
 364
 365       histname = "fHistJetsPtArea_";
 366       histname += i;
 367       fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, 50, 0, 1.5);
 368       fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 369       fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
 370       fOutput->Add(fHistJetsPtArea[i]);
 371     }
 372   }
 373
 374   Int_t dim = 0;
 375   TString title[20];
 376   Int_t nbins[20] = {0};
 377   Double_t min[20] = {0};
 378   Double_t max[20] = {0};
 379
 380   if (fForceBeamType != AliAnalysisTaskEmcal::kpp) {
 381     title[dim] = "Centrality %";
 382     nbins[dim] = 101;
 383     min[dim] = 0;
 384     max[dim] = 101;
 385     dim++;
 386
 387     if (!fCentMethod2.IsNull()) {
 388       title[dim] = Form("Centrality %s %%", fCentMethod2.Data());
 389       nbins[dim] = 101;
 390       min[dim] = 0;
 391       max[dim] = 101;
 392       dim++;
 393     }
 394
 395     if (!fCentMethod3.IsNull()) {
 396       title[dim] = Form("Centrality %s %%", fCentMethod3.Data());
 397       nbins[dim] = 101;
 398       min[dim] = 0;
 399       max[dim] = 101;
 400       dim++;
 401     }
 402
 403     if (fDoV0QA==1) {
 404       title[dim] = "V0A total multiplicity";
 405       nbins[dim] = 200;
 406       min[dim] = 0;
 407       max[dim] = 20000;
 408       dim++;
 409
 410       title[dim] = "V0C total multiplicity";
 411       nbins[dim] = 200;
 412       min[dim] = 0;
 413       max[dim] = 20000;
 414       dim++;
 415     }
 416     else if (fDoV0QA==2) {
 417       title[dim] = "V0A+V0C total multiplicity";
 418       nbins[dim] = 300;
 419       min[dim] = 0;
 420       max[dim] = 30000;
 421       dim++;
 422     }
 423
 424     if (!fRhoName.IsNull()) {
 425       title[dim] = "#rho (GeV/c)";
 426       nbins[dim] = fNbins*4;
 427       min[dim] = 0;
 428       max[dim] = 400;
 429       dim++;
 430     }
 431
 432     if (fDoEPQA) {
 433       title[dim] = "#psi_{RP}";
 434       nbins[dim] = 200;
 435       min[dim] = -TMath::Pi();
 436       max[dim] = TMath::Pi();
 437       dim++;
 438     }
 439   }
 440
 441   if (fParticleCollArray.GetEntriesFast()>0) {
 442     title[dim] = "No. of tracks";
 443     nbins[dim] = 3000;
 444     min[dim] = -0.5;
 445     max[dim] = 6000-0.5;
 446     dim++;
 447
 448     title[dim] = "p_{T,track}^{leading} (GeV/c)";
 449     nbins[dim] = fNbins;
 450     min[dim] = fMinBinPt;
 451     max[dim] = fMaxBinPt;
 452     dim++;
 453
 454     if (fDoLeadingObjectPosition) {
 455       title[dim] = "#eta_{track}^{leading}";
 456       nbins[dim] = 100;
 457       min[dim] = -1;
 458       max[dim] = 1;
 459       dim++;
 460
 461       title[dim] = "#phi_{track}^{leading}";
 462       nbins[dim] = 101;
 463       min[dim] = 0;
 464       max[dim] = TMath::Pi() * 2.02;
 465       dim++;
 466     }
 467   }
 468
 469   if (fClusterCollArray.GetEntriesFast()>0) {
 470     title[dim] = "No. of clusters";
 471     nbins[dim] = 2000;
 472     min[dim] = 0;
 473     max[dim] = 4000-0.5;
 474     dim++;
 475
 476     title[dim] = "E_{cluster}^{leading} (GeV)";
 477     nbins[dim] = fNbins;
 478     min[dim] = fMinBinPt;
 479     max[dim] = fMaxBinPt;
 480     dim++;
 481
 482     if (fDoLeadingObjectPosition) {
 483       title[dim] = "#eta_{cluster}^{leading}";
 484       nbins[dim] = 100;
 485       min[dim] = -1;
 486       max[dim] = 1;
 487       dim++;
 488
 489       title[dim] = "#phi_{cluster}^{leading}";
 490       nbins[dim] = 101;
 491       min[dim] = 0;
 492       max[dim] = TMath::Pi() * 2.02;
 493       dim++;
 494     }
 495   }
 496
 497   if (!fCaloCellsName.IsNull()) {
 498     title[dim] = "No. of cells";
 499     nbins[dim] = 3000;
 500     min[dim] = 0;
 501     max[dim] = 6000-0.5;
 502     dim++;
 503   }
 504
 505   if (fJetCollArray.GetEntriesFast()>0) {
 506     title[dim] = "No. of jets";
 507     nbins[dim] = 200;
 508     min[dim] = 0;
 509     max[dim] = 200-0.5;
 510     dim++;
 511
 512     title[dim] = "p_{T,jet}^{leading} (GeV/c)";
 513     nbins[dim] = fNbins;
 514     min[dim] = fMinBinPt;
 515     max[dim] = fMaxBinPt;
 516     dim++;
 517
 518     if (fDoLeadingObjectPosition) {
 519       title[dim] = "#eta_{jet}^{leading}";
 520       nbins[dim] = 100;
 521       min[dim] = -1;
 522       max[dim] = 1;
 523       dim++;
 524
 525       title[dim] = "#phi_{jet}^{leading}";
 526       nbins[dim] = 101;
 527       min[dim] = 0;
 528       max[dim] = TMath::Pi() * 2.02;
 529       dim++;
 530     }
 531   }
 532
 533   fHistEventQA = new THnSparseF("fHistEventQA","fHistEventQA",dim,nbins,min,max);
 534   for (Int_t i = 0; i < dim; i++)
 535     fHistEventQA->GetAxis(i)->SetTitle(title[i]);
 536   fOutput->Add(fHistEventQA);
 537
 538   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
 539 }
 540
 541 //________________________________________________________________________
 542 void AliAnalysisTaskSAQA::ExecOnce()
 543 {
 544   AliAnalysisTaskEmcalJet::ExecOnce();
 545
 546   if (fDoV0QA) {
 547     fVZERO = InputEvent()->GetVZEROData();
 548     if (!fVZERO) {
 549       AliError("AliVVZERO not available");
 550     }
 551   }
 552 }
 553
 554 //________________________________________________________________________
 555 Bool_t AliAnalysisTaskSAQA::RetrieveEventObjects()
 556 {
 557   // Retrieve event objects.
 558
 559   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 560     return kFALSE;
 561
 562   if (!fCentMethod2.IsNull() || !fCentMethod3.IsNull()) {
 563     if (fBeamType == kAA || fBeamType == kpA ) {
 564       AliCentrality *aliCent = InputEvent()->GetCentrality();
 565       if (aliCent) {
 566         if (!fCentMethod2.IsNull())
 567           fCent2 = aliCent->GetCentralityPercentile(fCentMethod2);
 568         if (!fCentMethod3.IsNull())
 569           fCent3 = aliCent->GetCentralityPercentile(fCentMethod3);
 570       }
 571     }
 572   }
 573
 574   if (fVZERO) {
 575     fV0ATotMult = AliESDUtils::GetCorrV0A(fVZERO->GetMTotV0A(),fVertex[2]);
 576     fV0CTotMult = AliESDUtils::GetCorrV0C(fVZERO->GetMTotV0C(),fVertex[2]);
 577   }
 578
 579   return kTRUE;
 580 }
 581
 582
 583 //________________________________________________________________________
 584 Bool_t AliAnalysisTaskSAQA::FillHistograms()
 585 {
 586   // Fill histograms.
 587
 588   Float_t trackSum = 0;
 589   Float_t clusSum = 0;
 590   Float_t cellSum = 0;
 591
 592   Int_t ntracks = 0;
 593   Int_t nclusters = 0;
 594   Int_t ncells = 0;
 595   Int_t njets = 0;
 596
 597   Float_t leadingClusE = 0;
 598   Float_t leadingClusEta = 0;
 599   Float_t leadingClusPhi = 0;
 600
 601   Float_t leadingTrackPt = 0;
 602   Float_t leadingTrackEta = 0;
 603   Float_t leadingTrackPhi = 0;
 604
 605   Float_t leadingJetPt = 0;
 606   Float_t leadingJetEta = 0;
 607   Float_t leadingJetPhi = 0;
 608
 609   if (fTracks) {
 610     AliVParticle *leadingTrack = 0;
 611
 612     ntracks = DoTrackLoop(trackSum, leadingTrack);
 613     AliDebug(2,Form("%d tracks found in the event", ntracks));
 614
 615     if (leadingTrack) {
 616       leadingTrackPt = leadingTrack->Pt();
 617       leadingTrackEta = leadingTrack->Eta();
 618       leadingTrackPhi = leadingTrack->Phi();
 619     }
 620   }
 621
 622   if (fCaloClusters) {
 623     AliVCluster  *leadingClus = 0;
 624
 625     nclusters = DoClusterLoop(clusSum, leadingClus);
 626     AliDebug(2,Form("%d clusters found in the event", nclusters));
 627
 628     if (leadingClus) {
 629       TLorentzVector leadingClusVect;
 630       leadingClus->GetMomentum(leadingClusVect, fVertex);
 631       leadingClusE = leadingClus->E();
 632       leadingClusEta = leadingClusVect.Eta();
 633       leadingClusPhi = leadingClusVect.Phi();
 634     }
 635   }
 636
 637   if (fCaloCells) {
 638     ncells = DoCellLoop(cellSum);
 639     AliDebug(2,Form("%d cells found in the event", ncells));
 640   }
 641
 642   if (fJets) {
 643     AliEmcalJet  *leadingJet = 0;
 644
 645     njets = DoJetLoop(leadingJet);
 646     AliDebug(2,Form("%d jets found in the event", njets));
 647
 648     if (leadingJet) {
 649       leadingJetPt = leadingJet->Pt();
 650       leadingJetEta = leadingJet->Eta();
 651       leadingJetPhi = leadingJet->Phi();
 652     }
 653   }
 654
 655   FillEventQAHisto(fCent, fCent2, fCent3, fV0ATotMult, fV0CTotMult, fEPV0, fRhoVal,
 656                    ntracks, nclusters, ncells, njets,
 657                    leadingTrackPt, leadingTrackEta, leadingTrackPhi,
 658                    leadingClusE, leadingClusEta, leadingClusPhi,
 659                    leadingJetPt, leadingJetEta, leadingJetPhi);
 660
 661   return kTRUE;
 662 }
 663
 664 //________________________________________________________________________
 665 void AliAnalysisTaskSAQA::FillEventQAHisto(Float_t cent, Float_t cent2, Float_t cent3, Float_t v0a, Float_t v0c,
 666                                            Float_t ep, Float_t rho, Int_t ntracks, Int_t nclusters, Int_t ncells, Int_t njets,
 667                                            Float_t maxTrackPt, Float_t maxTrackEta, Float_t maxTrackPhi,
 668                                            Float_t maxClusterE, Float_t maxClusterEta, Float_t maxClusterPhi,
 669                                            Float_t maxJetPt, Float_t maxJetEta, Float_t maxJetPhi)
 670 {
 671   Double_t contents[20]={0};
 672
 673   for (Int_t i = 0; i < fHistEventQA->GetNdimensions(); i++) {
 674     TString title(fHistEventQA->GetAxis(i)->GetTitle());
 675     if (title=="Centrality %")
 676       contents[i] = cent;
 677     else if (title==Form("Centrality %s %%", fCentMethod2.Data()))
 678       contents[i] = cent2;
 679     else if (title==Form("Centrality %s %%", fCentMethod3.Data()))
 680       contents[i] = cent3;
 681     else if (title=="V0A total multiplicity")
 682       contents[i] = v0a;
 683     else if (title=="V0C total multiplicity")
 684       contents[i] = v0c;
 685     else if (title=="V0A+V0C total multiplicity")
 686       contents[i] = v0a+v0c;
 687     else if (title=="#psi_{RP}")
 688       contents[i] = ep;
 689     else if (title=="#rho (GeV/c)")
 690       contents[i] = rho;
 691     else if (title=="No. of tracks")
 692         contents[i] = ntracks;
 693     else if (title=="No. of clusters")
 694       contents[i] = nclusters;
 695     else if (title=="No. of cells")
 696       contents[i] = ncells;
 697     else if (title=="No. of jets")
 698       contents[i] = njets;
 699     else if (title=="p_{T,track}^{leading} (GeV/c)")
 700       contents[i] = maxTrackPt;
 701     else if (title=="#eta_{track}^{leading}")
 702       contents[i] = maxTrackEta;
 703     else if (title=="#phi_{track}^{leading}")
 704       contents[i] = maxTrackPhi;
 705     else if (title=="E_{cluster}^{leading} (GeV)")
 706       contents[i] = maxClusterE;
 707     else if (title=="#eta_{cluster}^{leading}")
 708       contents[i] = maxClusterEta;
 709     else if (title=="#phi_{cluster}^{leading}")
 710       contents[i] = maxClusterPhi;
 711     else if (title=="p_{T,jet}^{leading} (GeV/c)")
 712       contents[i] = maxJetPt;
 713     else if (title=="#eta_{jet}^{leading}")
 714       contents[i] = maxJetEta;
 715     else if (title=="#phi_{jet}^{leading}")
 716       contents[i] = maxJetPhi;
 717     else
 718       AliWarning(Form("Unable to fill dimension %s!",title.Data()));
 719   }
 720
 721   fHistEventQA->Fill(contents);
 722 }
 723
 724 //________________________________________________________________________
 725 Int_t AliAnalysisTaskSAQA::DoCellLoop(Float_t &sum)
 726 {
 727   // Do cell loop.
 728
 729   AliVCaloCells *cells = InputEvent()->GetEMCALCells();
 730
 731   if (!cells)
 732     return 0;
 733
 734   const Int_t ncells = cells->GetNumberOfCells();
 735   Int_t nAccCells = 0;
 736
 737   for (Int_t pos = 0; pos < ncells; pos++) {
 738     Float_t amp   = cells->GetAmplitude(pos);
 739     Int_t   absId = cells->GetCellNumber(pos);
 740
 741     if (amp < fCellEnergyCut)
 742       continue;
 743
 744     fHistCellsAbsIdEnergy[fCentBin]->Fill(absId,amp);
 745     nAccCells++;
 746     sum += amp;
 747   }
 748
 749   return nAccCells;
 750 }
 751
 752 //________________________________________________________________________
 753 Double_t AliAnalysisTaskSAQA::GetFcross(AliVCluster *cluster, AliVCaloCells *cells)
 754 {
 755   Int_t    AbsIdseed  = -1;
 756   Double_t Eseed      = 0;
 757   for (Int_t i = 0; i < cluster->GetNCells(); i++) {
 758     if (cells->GetCellAmplitude(cluster->GetCellAbsId(i)) > AbsIdseed) {
 759       Eseed     = cells->GetCellAmplitude(cluster->GetCellAbsId(i));
 760       AbsIdseed = cluster->GetCellAbsId(i);
 761     }
 762   }
 763
 764   if (Eseed < 1e-9)
 765     return 100;
 766
 767   Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
 768   fGeom->GetCellIndex(AbsIdseed,imod,iTower,iIphi,iIeta);
 769   fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi,iIeta,iphi,ieta);
 770
 771   //Get close cells index and energy, not in corners
 772
 773   Int_t absID1 = -1;
 774   Int_t absID2 = -1;
 775
 776   if (iphi < AliEMCALGeoParams::fgkEMCALRows-1) absID1 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
 777   if (iphi > 0)                                 absID2 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
 778
 779   // In case of cell in eta = 0 border, depending on SM shift the cross cell index
 780
 781   Int_t absID3 = -1;
 782   Int_t absID4 = -1;
 783
 784   if (ieta == AliEMCALGeoParams::fgkEMCALCols-1 && !(imod%2)) {
 785     absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod+1, iphi, 0);
 786     absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod,   iphi, ieta-1);
 787   }
 788   else if (ieta == 0 && imod%2) {
 789     absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod,   iphi, ieta+1);
 790     absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod-1, iphi, AliEMCALGeoParams::fgkEMCALCols-1);
 791   }
 792   else {
 793     if (ieta < AliEMCALGeoParams::fgkEMCALCols-1)
 794       absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
 795     if (ieta > 0)
 796       absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
 797   }
 798
 799   Double_t  ecell1 = cells->GetCellAmplitude(absID1);
 800   Double_t  ecell2 = cells->GetCellAmplitude(absID2);
 801   Double_t  ecell3 = cells->GetCellAmplitude(absID3);
 802   Double_t  ecell4 = cells->GetCellAmplitude(absID4);
 803
 804   Double_t Ecross = ecell1 + ecell2 + ecell3 + ecell4;
 805
 806   Double_t Fcross = 1 - Ecross/Eseed;
 807
 808   return Fcross;
 809 }
 810
 811 //________________________________________________________________________
 812 Int_t AliAnalysisTaskSAQA::DoClusterLoop(Float_t &sum, AliVCluster* &leading)
 813 {
 814   // Do cluster loop.
 815
 816   if (!fCaloClusters)
 817     return 0;
 818
 819   Int_t nAccClusters = 0;
 820
 821   AliVCaloCells *cells = InputEvent()->GetEMCALCells();
 822
 823   sum = 0;
 824   leading = 0;
 825
 826   // Cluster loop
 827   const Int_t nclusters = fCaloClusters->GetEntriesFast();
 828
 829   for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
 830     AliVCluster* cluster = static_cast<AliVCluster*>(fCaloClusters->At(iClusters));
 831     if (!cluster) {
 832       AliError(Form("Could not receive cluster %d", iClusters));
 833       continue;
 834     }
 835
 836     if (!AcceptCluster(cluster))
 837       continue;
 838
 839     sum += cluster->E();
 840
 841     if (!leading || leading->E() < cluster->E()) leading = cluster;
 842
 843     TLorentzVector nPart;
 844     cluster->GetMomentum(nPart, fVertex);
 845
 846     fHistClusPhiEtaEnergy[fCentBin]->Fill(nPart.Eta(), nPart.Phi(), cluster->E());
 847
 848     Double_t ep = nPart.Phi() - fEPV0;
 849     while (ep < 0) ep += TMath::Pi();
 850     while (ep >= TMath::Pi()) ep -= TMath::Pi();
 851     fHistClusDeltaPhiEPEnergy[fCentBin]->Fill(cluster->E(), ep);
 852
 853     fHistNCellsEnergy[fCentBin]->Fill(cluster->E(), cluster->GetNCells());
 854
 855     fHistClusTimeEnergy[fCentBin]->Fill(cluster->E(), cluster->GetTOF());
 856
 857     if (cells) fHistFcrossEnergy[fCentBin]->Fill(cluster->E(), GetFcross(cluster, cells));
 858
 859     if (fHistClusMCEnergyFraction[fCentBin])
 860       fHistClusMCEnergyFraction[fCentBin]->Fill(cluster->GetMCEnergyFraction());
 861
 862     nAccClusters++;
 863   }
 864
 865   return nAccClusters;
 866 }
 867
 868 //________________________________________________________________________
 869 Int_t AliAnalysisTaskSAQA::DoTrackLoop(Float_t &sum, AliVParticle* &leading)
 870 {
 871   // Do track loop.
 872   if (!fTracks)
 873     return 0;
 874
 875   Int_t nAccTracks = 0;
 876
 877   sum = 0;
 878   leading = 0;
 879
 880   const Int_t ntracks = fTracks->GetEntriesFast();
 881   Int_t neg = 0;
 882   Int_t zero = 0;
 883
 884   for (Int_t i = 0; i < ntracks; i++) {
 885     AliVParticle* track = static_cast<AliVParticle*>(fTracks->At(i)); // pointer to reconstructed to track
 886
 887     if (!track) {
 888       AliError(Form("Could not retrieve track %d",i));
 889       continue;
 890     }
 891
 892     if (!AcceptTrack(track))
 893       continue;
 894
 895     nAccTracks++;
 896
 897     sum += track->P();
 898
 899     if (!leading || leading->Pt() < track->Pt()) leading = track;
 900
 901     if (fParticleLevel) {
 902       fHistTrPhiEtaPt[fCentBin][0]->Fill(track->Eta(), track->Phi(), track->Pt());
 903     }
 904     else {
 905       fHistTrPhiEtaPt[fCentBin][3]->Fill(track->Eta(), track->Phi(), track->Pt());
 906       if (track->GetLabel() == 0) {
 907         zero++;
 908         if (fHistTrPhiEtaZeroLab[fCentBin]) {
 909           fHistTrPhiEtaZeroLab[fCentBin]->Fill(track->Eta(), track->Phi());
 910           fHistTrPtZeroLab[fCentBin]->Fill(track->Pt());
 911         }
 912       }
 913
 914       if (track->GetLabel() < 0)
 915         neg++;
 916
 917       Int_t type = 0;
 918
 919       AliPicoTrack* ptrack = dynamic_cast<AliPicoTrack*>(track);
 920       if (ptrack)
 921         type = ptrack->GetTrackType();
 922
 923       if (type >= 0 && type < 3)
 924         fHistTrPhiEtaPt[fCentBin][type]->Fill(track->Eta(), track->Phi(), track->Pt());
 925       else
 926         AliDebug(2,Form("%s: track type %d not recognized!", GetName(), type));
 927     }
 928
 929     AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track);
 930
 931     if (!vtrack)
 932       continue;
 933
 934     if ((vtrack->GetTrackEtaOnEMCal() == -999 || vtrack->GetTrackPhiOnEMCal() == -999) && fHistTrPhiEtaNonProp[fCentBin]) {
 935       fHistTrPhiEtaNonProp[fCentBin]->Fill(vtrack->Eta(), vtrack->Phi());
 936       fHistTrPtNonProp[fCentBin]->Fill(vtrack->Pt());
 937     }
 938
 939     if (fHistTrEmcPhiEta[fCentBin])
 940       fHistTrEmcPhiEta[fCentBin]->Fill(vtrack->GetTrackEtaOnEMCal(), vtrack->GetTrackPhiOnEMCal());
 941     if (fHistTrEmcPt[fCentBin])
 942       fHistTrEmcPt[fCentBin]->Fill(vtrack->GetTrackPtOnEMCal());
 943     if (fHistDeltaEtaPt[fCentBin])
 944       fHistDeltaEtaPt[fCentBin]->Fill(vtrack->Pt(), vtrack->Eta() - vtrack->GetTrackEtaOnEMCal());
 945     if (fHistDeltaPhiPt[fCentBin])
 946       fHistDeltaPhiPt[fCentBin]->Fill(vtrack->Pt(), vtrack->Phi() - vtrack->GetTrackPhiOnEMCal());
 947     if (fHistDeltaPtvsPt[fCentBin])
 948       fHistDeltaPtvsPt[fCentBin]->Fill(vtrack->Pt(), vtrack->Pt() - vtrack->GetTrackPtOnEMCal());
 949   }
 950
 951   if (fHistTrNegativeLabels[fCentBin])
 952     fHistTrNegativeLabels[fCentBin]->Fill(1. * neg / ntracks);
 953
 954   if (fHistTrZeroLabels[fCentBin])
 955     fHistTrZeroLabels[fCentBin]->Fill(1. * zero / ntracks);
 956
 957   return nAccTracks;
 958 }
 959
 960 //________________________________________________________________________
 961 Int_t AliAnalysisTaskSAQA::DoJetLoop(AliEmcalJet* &leading)
 962 {
 963   // Do jet loop.
 964
 965   if (!fJets)
 966     return 0;
 967
 968   Int_t nAccJets = 0;
 969
 970   leading = 0;
 971
 972   Int_t njets = fJets->GetEntriesFast();
 973
 974   for (Int_t ij = 0; ij < njets; ij++) {
 975
 976     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(ij));
 977
 978     if (!jet) {
 979       AliError(Form("Could not receive jet %d", ij));
 980       continue;
 981     }
 982
 983     if (!AcceptJet(jet))
 984       continue;
 985
 986     if (!leading || leading->Pt() < jet->Pt()) leading = jet;
 987
 988     nAccJets++;
 989
 990     fHistJetsPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());
 991     fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
 992   }
 993
 994   return nAccJets;
 995 }