PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskSAQA.cxx

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