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