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