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