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 <TList.h>
  12 #include <TLorentzVector.h>
  13
  14 #include "AliAnalysisManager.h"
  15 #include "AliCentrality.h"
  16 #include "AliVCluster.h"
  17 #include "AliVParticle.h"
  18 #include "AliVTrack.h"
  19 #include "AliEmcalJet.h"
  20 #include "AliVEventHandler.h"
  21 #include "AliAODEvent.h"
  22 #include "AliExternalTrackParam.h"
  23 #include "AliTrackerBase.h"
  24 #include "AliLog.h"
  25 #include "AliEMCALGeometry.h"
  26 #include "AliEMCALGeoParams.h"
  27 #include "AliPicoTrack.h"
  28
  29 #include "AliAnalysisTaskSAQA.h"
  30
  31 ClassImp(AliAnalysisTaskSAQA)
  32
  33 //________________________________________________________________________
  34 AliAnalysisTaskSAQA::AliAnalysisTaskSAQA() :
  35   AliAnalysisTaskEmcalJet("AliAnalysisTaskSAQA", kTRUE),
  36   fCellEnergyCut(0.1),
  37   fParticleLevel(kFALSE),
  38   fIsMC(kFALSE),
  39   fNclusters(0),
  40   fNtracks(0),
  41   fNjets(0),
  42   fHistTracksCent(0),
  43   fHistClusCent(0),
  44   fHistJetsCent(0),
  45   fHistClusTracks(0),
  46   fHistJetsParts(0),
  47   fHistCellsCent(0),
  48   fHistCellsTracks(0),
  49   fHistTrNegativeLabels(0),
  50   fHistTrZeroLabels(0),
  51   fHistTrEmcPhiEta(0),
  52   fHistTrPhiEtaNonProp(0),
  53   fHistDeltaEtaPt(0),
  54   fHistDeltaPhiPt(0),
  55   fHistNCellsEnergy(0),
  56   fHistFcrossEnergy(0),
  57   fHistClusTimeEnergy(0),
  58   fHistCellsAbsIdEnergy(0),
  59   fHistChVSneCells(0),
  60   fHistChVSneClus(0),
  61   fHistChVSneCorrCells(0)
  62 {
  63   // Default constructor.
  64
  65   for (Int_t i = 0; i < 4; i++) {
  66     for (Int_t j = 0; j < 4; j++) fHistTrPhiEtaPt[i][j] = 0;
  67     fHistTrPhiEtaPtZeroLab[i] = 0;
  68     fHistClusPhiEtaEnergy[i] = 0;
  69     fHistJetsPhiEtaPt[i] = 0;
  70     fHistJetsPtArea[i] = 0;
  71   }
  72
  73   SetMakeGeneralHistograms(kTRUE);
  74 }
  75
  76 //________________________________________________________________________
  77 AliAnalysisTaskSAQA::AliAnalysisTaskSAQA(const char *name) :
  78   AliAnalysisTaskEmcalJet(name, kTRUE),
  79   fCellEnergyCut(0.1),
  80   fParticleLevel(kFALSE),
  81   fIsMC(kFALSE),
  82   fNclusters(0),
  83   fNtracks(0),
  84   fNjets(0),
  85   fHistTracksCent(0),
  86   fHistClusCent(0),
  87   fHistJetsCent(0),
  88   fHistClusTracks(0),
  89   fHistJetsParts(0),
  90   fHistCellsCent(0),
  91   fHistCellsTracks(0),
  92   fHistTrNegativeLabels(0),
  93   fHistTrZeroLabels(0),
  94   fHistTrEmcPhiEta(0),
  95   fHistTrPhiEtaNonProp(0),
  96   fHistDeltaEtaPt(0),
  97   fHistDeltaPhiPt(0),
  98   fHistNCellsEnergy(0),
  99   fHistFcrossEnergy(0),
 100   fHistClusTimeEnergy(0),
 101   fHistCellsAbsIdEnergy(0),
 102   fHistChVSneCells(0),
 103   fHistChVSneClus(0),
 104   fHistChVSneCorrCells(0)
 105 {
 106   // Standard constructor.
 107
 108   for (Int_t i = 0; i < 4; i++) {
 109     for (Int_t j = 0; j < 4; j++) fHistTrPhiEtaPt[i][j] = 0;
 110     fHistTrPhiEtaPtZeroLab[i] = 0;
 111     fHistClusPhiEtaEnergy[i] = 0;
 112     fHistJetsPhiEtaPt[i] = 0;
 113     fHistJetsPtArea[i] = 0;
 114   }
 115
 116   SetMakeGeneralHistograms(kTRUE);
 117 }
 118
 119 //________________________________________________________________________
 120 AliAnalysisTaskSAQA::~AliAnalysisTaskSAQA()
 121 {
 122   // Destructor
 123 }
 124
 125 //________________________________________________________________________
 126 void AliAnalysisTaskSAQA::UserCreateOutputObjects()
 127 {
 128   // Create histograms
 129
 130   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 131
 132   if (!fTracksName.IsNull()) {
 133     if (!fParticleLevel && fIsMC) {
 134       fHistTrNegativeLabels = new TH1F("fHistTrNegativeLabels","fHistTrNegativeLabels", 500, 0, 1);
 135       fHistTrNegativeLabels->GetXaxis()->SetTitle("% of negative labels");
 136       fHistTrNegativeLabels->GetYaxis()->SetTitle("counts");
 137       fOutput->Add(fHistTrNegativeLabels);
 138
 139       fHistTrZeroLabels = new TH1F("fHistTrZeroLabels","fHistTrZeroLabels", 500, 0, 1);
 140       fHistTrZeroLabels->GetXaxis()->SetTitle("% of negative labels");
 141       fHistTrZeroLabels->GetYaxis()->SetTitle("counts");
 142       fOutput->Add(fHistTrZeroLabels);
 143     }
 144
 145     fHistTracksCent = new TH2F("fHistTracksCent","Tracks vs. centrality", 100, 0, 100, fNbins, 0, 4000);
 146     fHistTracksCent->GetXaxis()->SetTitle("Centrality (%)");
 147     fHistTracksCent->GetYaxis()->SetTitle("No. of tracks");
 148     fOutput->Add(fHistTracksCent);
 149
 150     TString histname;
 151
 152     Int_t nlabels = 4;
 153     if (fParticleLevel)
 154       nlabels = 1;
 155
 156     for (Int_t i = 0; i < fNcentBins; i++) {
 157       for (Int_t j = 0; j < nlabels; j++) {
 158         histname = Form("fHistTrPhiEtaPt_%d_%d",i,j);
 159         fHistTrPhiEtaPt[i][j] = new TH3F(histname,histname, 100, -1, 1, 201, 0, TMath::Pi() * 2.01, fNbins, fMinBinPt, fMaxBinPt);
 160         fHistTrPhiEtaPt[i][j]->GetXaxis()->SetTitle("#eta");
 161         fHistTrPhiEtaPt[i][j]->GetYaxis()->SetTitle("#phi");
 162         fHistTrPhiEtaPt[i][j]->GetZaxis()->SetTitle("p_{T} (GeV/c)");
 163         fOutput->Add(fHistTrPhiEtaPt[i][j]);
 164       }
 165       if (!fParticleLevel && fIsMC) {
 166         histname = Form("fHistTrPhiEtaPtZeroLab_%d",i);
 167         fHistTrPhiEtaPtZeroLab[i] = new TH3F(histname,histname, 100, -1, 1, 201, 0, TMath::Pi() * 2.01, fNbins, fMinBinPt, fMaxBinPt);
 168         fHistTrPhiEtaPtZeroLab[i]->GetXaxis()->SetTitle("#eta");
 169         fHistTrPhiEtaPtZeroLab[i]->GetYaxis()->SetTitle("#phi");
 170         fHistTrPhiEtaPtZeroLab[i]->GetZaxis()->SetTitle("p_{T} (GeV/c)");
 171         fOutput->Add(fHistTrPhiEtaPtZeroLab[i]);
 172       }
 173     }
 174
 175     if (!fParticleLevel) {
 176       fHistTrEmcPhiEta = new TH2F("fHistTrEmcPhiEta","Phi-Eta emcal distribution of tracks", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
 177       fHistTrEmcPhiEta->GetXaxis()->SetTitle("#eta");
 178       fHistTrEmcPhiEta->GetYaxis()->SetTitle("#phi");
 179       fOutput->Add(fHistTrEmcPhiEta);
 180
 181       fHistTrPhiEtaNonProp = new TH2F("fHistTrPhiEtaNonProp","fHistTrPhiEtaNonProp", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
 182       fHistTrPhiEtaNonProp->GetXaxis()->SetTitle("#eta");
 183       fHistTrPhiEtaNonProp->GetYaxis()->SetTitle("#phi");
 184       fOutput->Add(fHistTrPhiEtaNonProp);
 185
 186       fHistDeltaEtaPt = new TH2F("fHistDeltaEtaPt","fHistDeltaEtaPt", fNbins, fMinBinPt, fMaxBinPt, 80, -0.5, 0.5);
 187       fHistDeltaEtaPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 188       fHistDeltaEtaPt->GetYaxis()->SetTitle("#delta#eta");
 189       fOutput->Add(fHistDeltaEtaPt);
 190
 191       fHistDeltaPhiPt = new TH2F("fHistDeltaPhiPt","fHistDeltaPhiPt", fNbins, fMinBinPt, fMaxBinPt, 256, -1.6, 4.8);
 192       fHistDeltaPhiPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 193       fHistDeltaPhiPt->GetYaxis()->SetTitle("#delta#phi");
 194       fOutput->Add(fHistDeltaPhiPt);
 195     }
 196   }
 197
 198   if (!fCaloName.IsNull()) {
 199     fHistClusCent = new TH2F("fHistClusCent","Clusters vs. centrality", 100, 0, 100, fNbins, 0, 2000);
 200     fHistClusCent->GetXaxis()->SetTitle("Centrality (%)");
 201     fHistClusCent->GetYaxis()->SetTitle("No. of clusters");
 202     fOutput->Add(fHistClusCent);
 203
 204     fHistClusTracks = new TH2F("fHistClusTracks","Clusters vs. tracks", fNbins, 0, 4000, fNbins, 0, 2000);
 205     fHistClusTracks->GetXaxis()->SetTitle("No. of tracks");
 206     fHistClusTracks->GetYaxis()->SetTitle("No. of clusters");
 207     fOutput->Add(fHistClusTracks);
 208
 209     fHistCellsCent = new TH2F("fHistCellsCent","Cells vs. centrality", 100, 0, 100, fNbins, 0, 6000);
 210     fHistCellsCent->GetXaxis()->SetTitle("Centrality (%)");
 211     fHistCellsCent->GetYaxis()->SetTitle("No. of EMCal cells");
 212     fOutput->Add(fHistCellsCent);
 213
 214     fHistCellsTracks = new TH2F("fHistCellsTracks","Cells vs. tracks", fNbins, 0, 4000, fNbins, 0, 6000);
 215     fHistCellsTracks->GetXaxis()->SetTitle("No. of tracks");
 216     fHistCellsTracks->GetYaxis()->SetTitle("No. of EMCal cells");
 217     fOutput->Add(fHistCellsTracks);
 218
 219     TString histname;
 220
 221     for (Int_t i = 0; i < fNcentBins; i++) {
 222       histname = "fHistClusPhiEtaEnergy_";
 223       histname += i;
 224       fHistClusPhiEtaEnergy[i] = new TH3F(histname, histname, 100, -1.2, 1.2, 201, 0, TMath::Pi() * 2.01, fNbins, fMinBinPt, fMaxBinPt);
 225       fHistClusPhiEtaEnergy[i]->GetXaxis()->SetTitle("#eta");
 226       fHistClusPhiEtaEnergy[i]->GetYaxis()->SetTitle("#phi");
 227       fHistClusPhiEtaEnergy[i]->GetZaxis()->SetTitle("E_{cluster} (GeV)");
 228       fOutput->Add(fHistClusPhiEtaEnergy[i]);
 229     }
 230
 231     fHistClusTimeEnergy = new TH2F("fHistClusTimeEnergy","Time vs. energy of clusters", fNbins, fMinBinPt, fMaxBinPt, fNbins,  -1e-6, 1e-6);
 232     fHistClusTimeEnergy->GetXaxis()->SetTitle("E_{cluster} (GeV)");
 233     fHistClusTimeEnergy->GetYaxis()->SetTitle("Time");
 234     fOutput->Add(fHistClusTimeEnergy);
 235
 236     fHistNCellsEnergy = new TH2F("fHistNCellsEnergy","Number of cells vs. energy of clusters", fNbins, fMinBinPt, fMaxBinPt, 30, 0, 30);
 237     fHistNCellsEnergy->GetXaxis()->SetTitle("E_{cluster} (GeV)");
 238     fHistNCellsEnergy->GetYaxis()->SetTitle("N_{cells}");
 239     fOutput->Add(fHistNCellsEnergy);
 240
 241     fHistFcrossEnergy = new TH2F("fHistFcrossEnergy","fHistFcrossEnergy", fNbins, fMinBinPt, fMaxBinPt, 200, -3.5, 1.5);
 242     fHistFcrossEnergy->GetXaxis()->SetTitle("E_{cluster} (GeV)");
 243     fHistFcrossEnergy->GetYaxis()->SetTitle("F_{cross}");
 244     fOutput->Add(fHistFcrossEnergy);
 245
 246     fHistCellsAbsIdEnergy = new TH2F("fHistCellsAbsIdEnergy","fHistCellsAbsIdEnergy", 11600,0,11599,(Int_t)(fNbins / 2), fMinBinPt, fMaxBinPt / 2);
 247     fHistCellsAbsIdEnergy->GetXaxis()->SetTitle("cell abs. Id");
 248     fHistCellsAbsIdEnergy->GetYaxis()->SetTitle("E_{cluster} (GeV)");
 249     fHistCellsAbsIdEnergy->GetZaxis()->SetTitle("counts");
 250     fOutput->Add(fHistCellsAbsIdEnergy);
 251
 252     fHistChVSneCells = new TH2F("fHistChVSneCells","Charged energy vs. neutral (cells) energy",
 253                                 (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
 254     fHistChVSneCells->GetXaxis()->SetTitle("Energy (GeV)");
 255     fHistChVSneCells->GetYaxis()->SetTitle("Momentum (GeV/c)");
 256     fOutput->Add(fHistChVSneCells);
 257
 258     fHistChVSneClus = new TH2F("fHistChVSneClus","Charged energy vs. neutral (clusters) energy",
 259                                (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
 260     fHistChVSneClus->GetXaxis()->SetTitle("Energy (GeV)");
 261     fHistChVSneClus->GetYaxis()->SetTitle("Momentum (GeV/c)");
 262     fOutput->Add(fHistChVSneClus);
 263
 264     fHistChVSneCorrCells = new TH2F("fHistChVSneCorrCells","Charged energy vs. neutral (corrected cells) energy",
 265                                     (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, (Int_t)(fNbins * 2.5), fMinBinPt , fMaxBinPt * 2.5);
 266     fHistChVSneCorrCells->GetXaxis()->SetTitle("Energy (GeV)");
 267     fHistChVSneCorrCells->GetYaxis()->SetTitle("Momentum (GeV/c)");
 268     fOutput->Add(fHistChVSneCorrCells);
 269   }
 270
 271   if (!fJetsName.IsNull()) {
 272     fHistJetsCent = new TH2F("fHistJetsCent","Jets vs. centrality", 100, 0, 100, 150, -0.5, 149.5);
 273     fHistJetsCent->GetXaxis()->SetTitle("Centrality (%)");
 274     fHistJetsCent->GetYaxis()->SetTitle("No. of jets");
 275     fOutput->Add(fHistJetsCent);
 276
 277     fHistJetsParts = new TH2F("fHistJetsParts","Jets vs. centrality", fNbins, 0, 6000, 150, -0.5, 149.5);
 278     fHistJetsParts->GetXaxis()->SetTitle("No. of particles");
 279     fHistJetsParts->GetYaxis()->SetTitle("No. of jets");
 280     fOutput->Add(fHistJetsParts);
 281
 282     TString histname;
 283
 284     for (Int_t i = 0; i < fNcentBins; i++) {
 285       histname = "fHistJetsPhiEtaPt_";
 286       histname += i;
 287       fHistJetsPhiEtaPt[i] = new TH3F(histname.Data(), histname.Data(), 100, -1.2, 1.2, 201, 0, TMath::Pi() * 2.01,  (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
 288       fHistJetsPhiEtaPt[i]->GetXaxis()->SetTitle("#eta");
 289       fHistJetsPhiEtaPt[i]->GetYaxis()->SetTitle("#phi");
 290       fHistJetsPhiEtaPt[i]->GetZaxis()->SetTitle("p_{T} (GeV/c)");
 291       fOutput->Add(fHistJetsPhiEtaPt[i]);
 292
 293       histname = "fHistJetsPtArea_";
 294       histname += i;
 295       fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, 30, 0, fJetRadius * fJetRadius * TMath::Pi() * 3);
 296       fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 297       fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
 298       fOutput->Add(fHistJetsPtArea[i]);
 299     }
 300   }
 301
 302   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
 303 }
 304
 305 //________________________________________________________________________
 306 Bool_t AliAnalysisTaskSAQA::RetrieveEventObjects()
 307 {
 308   // Retrieve event objects.
 309
 310   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 311     return kFALSE;
 312
 313   fNclusters = 0;
 314   fNtracks = 0;
 315   fNjets = 0;
 316
 317   return kTRUE;
 318 }
 319
 320
 321 //________________________________________________________________________
 322 Bool_t AliAnalysisTaskSAQA::FillHistograms()
 323 {
 324   // Fill histograms.
 325
 326   Float_t clusSum = 0;
 327   Float_t trackSum = 0;
 328
 329   if (fTracks) {
 330     trackSum = DoTrackLoop();
 331     AliDebug(2,Form("%d tracks found in the event", fTracks->GetEntriesFast()));
 332     fHistTracksCent->Fill(fCent, fNtracks);
 333   }
 334
 335   if (fCaloClusters) {
 336     clusSum = DoClusterLoop();
 337
 338     fHistClusCent->Fill(fCent, fNclusters);
 339     fHistClusTracks->Fill(fNtracks, fNclusters);
 340
 341     Float_t cellSum = 0, cellCutSum = 0;
 342
 343     Int_t ncells = DoCellLoop(cellSum, cellCutSum);
 344
 345     if (fTracks)
 346       fHistCellsTracks->Fill(fNtracks, ncells);
 347
 348     fHistCellsCent->Fill(fCent, ncells);
 349
 350     fHistChVSneCells->Fill(cellSum, trackSum);
 351     fHistChVSneClus->Fill(clusSum, trackSum);
 352     fHistChVSneCorrCells->Fill(cellCutSum, trackSum);
 353   }
 354
 355   if (fJets) {
 356     DoJetLoop();
 357
 358     fHistJetsCent->Fill(fCent, fNjets);
 359     fHistJetsParts->Fill(fNtracks + fNclusters, fNjets);
 360   }
 361
 362   return kTRUE;
 363 }
 364
 365 //________________________________________________________________________
 366 Int_t AliAnalysisTaskSAQA::DoCellLoop(Float_t &sum, Float_t &sum_cut)
 367 {
 368   // Do cell loop.
 369
 370   AliVCaloCells *cells = InputEvent()->GetEMCALCells();
 371
 372   if (!cells)
 373     return 0;
 374
 375   const Int_t ncells = cells->GetNumberOfCells();
 376
 377   for (Int_t pos = 0; pos < ncells; pos++) {
 378     Float_t amp   = cells->GetAmplitude(pos);
 379     Int_t   absId = cells->GetCellNumber(pos);
 380     fHistCellsAbsIdEnergy->Fill(absId,amp);
 381     sum += amp;
 382     if (amp < fCellEnergyCut)
 383       continue;
 384     sum_cut += amp;
 385   }
 386
 387   return ncells;
 388 }
 389
 390 //________________________________________________________________________
 391 Double_t AliAnalysisTaskSAQA::GetFcross(AliVCluster *cluster, AliVCaloCells *cells)
 392 {
 393   Int_t    AbsIdseed  = -1;
 394   Double_t Eseed      = 0;
 395   for (Int_t i = 0; i < cluster->GetNCells(); i++) {
 396     if (cells->GetCellAmplitude(cluster->GetCellAbsId(i)) > AbsIdseed) {
 397       Eseed     = cells->GetCellAmplitude(cluster->GetCellAbsId(i));
 398       AbsIdseed = cluster->GetCellAbsId(i);
 399     }
 400   }
 401
 402   if (Eseed < 1e-9)
 403     return 100;
 404
 405   Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
 406   fGeom->GetCellIndex(AbsIdseed,imod,iTower,iIphi,iIeta);
 407   fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi,iIeta,iphi,ieta);
 408
 409   //Get close cells index and energy, not in corners
 410
 411   Int_t absID1 = -1;
 412   Int_t absID2 = -1;
 413
 414   if (iphi < AliEMCALGeoParams::fgkEMCALRows-1) absID1 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
 415   if (iphi > 0)                                 absID2 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
 416
 417   // In case of cell in eta = 0 border, depending on SM shift the cross cell index
 418
 419   Int_t absID3 = -1;
 420   Int_t absID4 = -1;
 421
 422   if (ieta == AliEMCALGeoParams::fgkEMCALCols-1 && !(imod%2)) {
 423     absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod+1, iphi, 0);
 424     absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod,   iphi, ieta-1);
 425   }
 426   else if (ieta == 0 && imod%2) {
 427     absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod,   iphi, ieta+1);
 428     absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod-1, iphi, AliEMCALGeoParams::fgkEMCALCols-1);
 429   }
 430   else {
 431     if (ieta < AliEMCALGeoParams::fgkEMCALCols-1)
 432       absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
 433     if (ieta > 0)
 434       absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
 435   }
 436
 437   Double_t  ecell1 = cells->GetCellAmplitude(absID1);
 438   Double_t  ecell2 = cells->GetCellAmplitude(absID2);
 439   Double_t  ecell3 = cells->GetCellAmplitude(absID3);
 440   Double_t  ecell4 = cells->GetCellAmplitude(absID4);
 441
 442   Double_t Ecross = ecell1 + ecell2 + ecell3 + ecell4;
 443
 444   Double_t Fcross = 1 - Ecross/Eseed;
 445
 446   return Fcross;
 447 }
 448
 449 //________________________________________________________________________
 450 Float_t AliAnalysisTaskSAQA::DoClusterLoop()
 451 {
 452   // Do cluster loop.
 453
 454   if (!fCaloClusters)
 455     return 0;
 456
 457   AliVCaloCells *cells = InputEvent()->GetEMCALCells();
 458
 459   Float_t sum = 0;
 460
 461   // Cluster loop
 462   Int_t nclusters =  fCaloClusters->GetEntriesFast();
 463
 464   for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
 465     AliVCluster* cluster = static_cast<AliVCluster*>(fCaloClusters->At(iClusters));
 466     if (!cluster) {
 467       AliError(Form("Could not receive cluster %d", iClusters));
 468       continue;
 469     }
 470
 471     if (!AcceptCluster(cluster))
 472       continue;
 473
 474     sum += cluster->E();
 475
 476     TLorentzVector nPart;
 477     cluster->GetMomentum(nPart, fVertex);
 478
 479     fHistClusPhiEtaEnergy[fCentBin]->Fill(nPart.Eta(), nPart.Phi(), cluster->E());
 480     fHistNCellsEnergy->Fill(cluster->E(), cluster->GetNCells());
 481
 482     fHistClusTimeEnergy->Fill(cluster->E(), cluster->GetTOF());
 483
 484     if (cells)
 485       fHistFcrossEnergy->Fill(cluster->E(), GetFcross(cluster, cells));
 486
 487     fNclusters++;
 488   }
 489
 490   return sum;
 491 }
 492
 493 //________________________________________________________________________
 494 Float_t AliAnalysisTaskSAQA::DoTrackLoop()
 495 {
 496   // Do track loop.
 497
 498   if (!fTracks)
 499     return 0;
 500
 501   Float_t sum = 0;
 502
 503   Int_t ntracks = fTracks->GetEntriesFast();
 504   Int_t neg = 0;
 505   Int_t zero = 0;
 506
 507   for (Int_t i = 0; i < ntracks; i++) {
 508
 509     AliVParticle* track = static_cast<AliVParticle*>(fTracks->At(i)); // pointer to reconstructed to track
 510
 511     if (!track) {
 512       AliError(Form("Could not retrieve track %d",i));
 513       continue;
 514     }
 515
 516     if (!AcceptTrack(track))
 517       continue;
 518
 519     fNtracks++;
 520
 521     sum += track->P();
 522
 523     if (fParticleLevel) {
 524       fHistTrPhiEtaPt[fCentBin][0]->Fill(track->Eta(), track->Phi(), track->Pt());
 525     }
 526     else {
 527       fHistTrPhiEtaPt[fCentBin][3]->Fill(track->Eta(), track->Phi(), track->Pt());
 528       if (track->GetLabel() == 0) {
 529         zero++;
 530         if (fHistTrPhiEtaPtZeroLab[fCentBin])
 531           fHistTrPhiEtaPtZeroLab[fCentBin]->Fill(track->Eta(), track->Phi(), track->Pt());
 532       }
 533
 534       if (track->GetLabel() < 0)
 535         neg++;
 536
 537       Int_t type = 0;
 538
 539       AliPicoTrack* ptrack = dynamic_cast<AliPicoTrack*>(track);
 540       if (ptrack)
 541         type = ptrack->GetTrackType();
 542
 543       if (type >= 0 && type < 3)
 544         fHistTrPhiEtaPt[fCentBin][type]->Fill(track->Eta(), track->Phi(), track->Pt());
 545       else
 546         AliDebug(2,Form("%s: track type %d not recognized!", GetName(), type));
 547     }
 548
 549     AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track);
 550
 551     if (!vtrack)
 552       continue;
 553
 554     if ((vtrack->GetTrackEtaOnEMCal() == -999 || vtrack->GetTrackPhiOnEMCal() == -999) && fHistTrPhiEtaNonProp)
 555       fHistTrPhiEtaNonProp->Fill(vtrack->Eta(), vtrack->Phi());
 556
 557     if (fHistTrEmcPhiEta)
 558       fHistTrEmcPhiEta->Fill(vtrack->GetTrackEtaOnEMCal(), vtrack->GetTrackPhiOnEMCal());
 559     if (fHistDeltaEtaPt)
 560       fHistDeltaEtaPt->Fill(vtrack->Pt(), vtrack->Eta() - vtrack->GetTrackEtaOnEMCal());
 561     if (fHistDeltaPhiPt)
 562       fHistDeltaPhiPt->Fill(vtrack->Pt(), vtrack->Phi() - vtrack->GetTrackPhiOnEMCal());
 563   }
 564
 565   if (fHistTrNegativeLabels)
 566     fHistTrNegativeLabels->Fill(1. * neg / ntracks);
 567
 568   if (fHistTrZeroLabels)
 569     fHistTrZeroLabels->Fill(1. * zero / ntracks);
 570
 571   return sum;
 572 }
 573
 574 //________________________________________________________________________
 575 void AliAnalysisTaskSAQA::DoJetLoop()
 576 {
 577   // Do jet loop.
 578
 579   if (!fJets)
 580     return;
 581
 582   Int_t njets = fJets->GetEntriesFast();
 583
 584   for (Int_t ij = 0; ij < njets; ij++) {
 585
 586     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(ij));
 587
 588     if (!jet) {
 589       AliError(Form("Could not receive jet %d", ij));
 590       continue;
 591     }
 592
 593     if (!AcceptJet(jet))
 594       continue;
 595
 596     fNjets++;
 597
 598     fHistJetsPhiEtaPt[fCentBin]->Fill(jet->Eta(), jet->Phi(), jet->Pt());
 599     fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
 600   }
 601 }
 602
 603
 604 //________________________________________________________________________
 605 void AliAnalysisTaskSAQA::Terminate(Option_t *)
 606 {
 607   // Called once at the end of the analysis.
 608 }