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
  27 #include "AliAnalysisTaskSAQA.h"
  28
  29 ClassImp(AliAnalysisTaskSAQA)
  30
  31 //________________________________________________________________________
  32 AliAnalysisTaskSAQA::AliAnalysisTaskSAQA() :
  33   AliAnalysisTaskEmcalJet("AliAnalysisTaskSAQA", kTRUE),
  34   fCellEnergyCut(0.1),
  35   fDoTrigger(kFALSE),
  36   fRepropagateTracks(kFALSE),
  37   fTrgClusName("ClustersL1GAMMAFEE"),
  38   fTrgClusters(0),
  39   fHistCentrality(0),
  40   fHistZVertex(0),
  41   fHistTracksCent(0),
  42   fHistClusCent(0),
  43   fHistClusTracks(0),
  44   fHistCellsCent(0),
  45   fHistCellsTracks(0),
  46   fHistMaxL1FastORCent(0),
  47   fHistMaxL1ClusCent(0),
  48   fHistMaxL1ThrCent(0),
  49   fHistTracksPt(0),
  50   fHistTrPhiEta(0),
  51   fHistTrEmcPhiEta(0),
  52   fHistTrPhiEtaNonProp(0),
  53   fHistDeltaEtaPt(0),
  54   fHistDeltaPhiPt(0),
  55   fHistDeltaEtaNewProp(0),
  56   fHistDeltaPhiNewProp(0),
  57   fHistClusPhiEtaEnergy(0),
  58   fHistNCellsEnergy(0),
  59   fHistClusTimeEnergy(0),
  60   fHistCellsEnergy(0),
  61   fHistChVSneCells(0),
  62   fHistChVSneClus(0),
  63   fHistChVSneCorrCells(0)
  64 {
  65   // Default constructor.
  66
  67   for (Int_t i = 0; i < 5; i++) {
  68     fHistTrackPhi[i] = 0;
  69     fHistTrackEta[i] = 0;
  70   }
  71
  72   for (Int_t i = 0; i < 4; i++) {
  73     fHistJetsPhiEta[i] = 0;
  74     fHistJetsPtNonBias[i] = 0;
  75     fHistJetsPtTrack[i] = 0;
  76     fHistJetsPtClus[i] = 0;
  77     fHistJetsPt[i] = 0;
  78     fHistJetsPtArea[i] = 0;
  79     fHistJetsPtAreaNonBias[i] = 0;
  80   }
  81 }
  82
  83 //________________________________________________________________________
  84 AliAnalysisTaskSAQA::AliAnalysisTaskSAQA(const char *name) :
  85   AliAnalysisTaskEmcalJet(name, kTRUE),
  86   fCellEnergyCut(0.1),
  87   fDoTrigger(kFALSE),
  88   fRepropagateTracks(kFALSE),
  89   fTrgClusName("ClustersL1GAMMAFEE"),
  90   fTrgClusters(0),
  91   fHistCentrality(0),
  92   fHistZVertex(0),
  93   fHistTracksCent(0),
  94   fHistClusCent(0),
  95   fHistClusTracks(0),
  96   fHistCellsCent(0),
  97   fHistCellsTracks(0),
  98   fHistMaxL1FastORCent(0),
  99   fHistMaxL1ClusCent(0),
 100   fHistMaxL1ThrCent(0),
 101   fHistTracksPt(0),
 102   fHistTrPhiEta(0),
 103   fHistTrEmcPhiEta(0),
 104   fHistTrPhiEtaNonProp(0),
 105   fHistDeltaEtaPt(0),
 106   fHistDeltaPhiPt(0),
 107   fHistDeltaEtaNewProp(0),
 108   fHistDeltaPhiNewProp(0),
 109   fHistClusPhiEtaEnergy(0),
 110   fHistNCellsEnergy(0),
 111   fHistClusTimeEnergy(0),
 112   fHistCellsEnergy(0),
 113   fHistChVSneCells(0),
 114   fHistChVSneClus(0),
 115   fHistChVSneCorrCells(0)
 116 {
 117   // Standard constructor.
 118
 119   for (Int_t i = 0; i < 5; i++) {
 120     fHistTrackPhi[i] = 0;
 121     fHistTrackEta[i] = 0;
 122   }
 123
 124   for (Int_t i = 0; i < 4; i++) {
 125     fHistJetsPhiEta[i] = 0;
 126     fHistJetsPtNonBias[i] = 0;
 127     fHistJetsPtTrack[i] = 0;
 128     fHistJetsPtClus[i] = 0;
 129     fHistJetsPt[i] = 0;
 130     fHistJetsPtArea[i] = 0;
 131     fHistJetsPtAreaNonBias[i] = 0;
 132   }
 133 }
 134
 135 //________________________________________________________________________
 136 AliAnalysisTaskSAQA::~AliAnalysisTaskSAQA()
 137 {
 138   // Destructor
 139 }
 140
 141 //________________________________________________________________________
 142 void AliAnalysisTaskSAQA::UserCreateOutputObjects()
 143 {
 144   // Create histograms
 145
 146   OpenFile(1);
 147   fOutput = new TList();
 148   fOutput->SetOwner();
 149
 150   fHistCentrality = new TH1F("fHistCentrality","Event centrality distribution", 100, 0, 100);
 151   fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
 152   fHistCentrality->GetYaxis()->SetTitle("counts");
 153   fOutput->Add(fHistCentrality);
 154
 155   fHistZVertex = new TH1F("fHistZVertex","Z vertex position", 60, -30, 30);
 156   fHistZVertex->GetXaxis()->SetTitle("z");
 157   fHistZVertex->GetYaxis()->SetTitle("counts");
 158   fOutput->Add(fHistZVertex);
 159
 160   fHistTracksCent = new TH2F("fHistTracksCent","Tracks vs. centrality", 100, 0, 100, fNbins, 0, 4000);
 161   fHistTracksCent->GetXaxis()->SetTitle("Centrality (%)");
 162   fHistTracksCent->GetYaxis()->SetTitle("No. of tracks");
 163   fOutput->Add(fHistTracksCent);
 164
 165   if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {
 166     fHistClusCent = new TH2F("fHistClusCent","Clusters vs. centrality", 100, 0, 100, fNbins, 0, 2000);
 167     fHistClusCent->GetXaxis()->SetTitle("Centrality (%)");
 168     fHistClusCent->GetYaxis()->SetTitle("No. of clusters");
 169     fOutput->Add(fHistClusCent);
 170
 171     fHistClusTracks = new TH2F("fHistClusTracks","Clusters vs. tracks", fNbins, 0, 4000, fNbins, 0, 2000);
 172     fHistClusTracks->GetXaxis()->SetTitle("No. of tracks");
 173     fHistClusTracks->GetYaxis()->SetTitle("No. of clusters");
 174     fOutput->Add(fHistClusTracks);
 175
 176     fHistCellsCent = new TH2F("fHistCellsCent","Cells vs. centrality", 100, 0, 100, fNbins, 0, 6000);
 177     fHistCellsCent->GetXaxis()->SetTitle("Centrality (%)");
 178     fHistCellsCent->GetYaxis()->SetTitle("No. of EMCal cells");
 179     fOutput->Add(fHistCellsCent);
 180
 181     fHistCellsTracks = new TH2F("fHistCellsTracks","Cells vs. tracks", fNbins, 0, 4000, fNbins, 0, 6000);
 182     fHistCellsTracks->GetXaxis()->SetTitle("No. of tracks");
 183     fHistCellsTracks->GetYaxis()->SetTitle("No. of EMCal cells");
 184     fOutput->Add(fHistCellsTracks);
 185
 186     fHistClusTimeEnergy = new TH2F("fHistClusTimeEnergy","Time vs. energy of clusters", fNbins, fMinBinPt, fMaxBinPt, fNbins,  0, 1e-6);
 187     fHistClusTimeEnergy->GetXaxis()->SetTitle("Energy (GeV)");
 188     fHistClusTimeEnergy->GetYaxis()->SetTitle("Time");
 189     fOutput->Add(fHistClusTimeEnergy);
 190
 191     if (fDoTrigger) {
 192       fHistMaxL1FastORCent = new TH2F("fHistMaxL1FastORCent","fHistMaxL1ClusCent", 100, 0, 100, 250, 0, 250);
 193       fHistMaxL1FastORCent->GetXaxis()->SetTitle("Centrality [%]");
 194       fHistMaxL1FastORCent->GetYaxis()->SetTitle("Maximum L1 FastOR");
 195       fOutput->Add(fHistMaxL1FastORCent);
 196
 197       fHistMaxL1ClusCent = new TH2F("fHistMaxL1ClusCent","fHistMaxL1ClusCent", 100, 0, 100, 250, 0, 250);
 198       fHistMaxL1ClusCent->GetXaxis()->SetTitle("Centrality [%]");
 199       fHistMaxL1ClusCent->GetYaxis()->SetTitle("Maximum L1 trigger cluster");
 200       fOutput->Add(fHistMaxL1ClusCent);
 201
 202       fHistMaxL1ThrCent = new TH2F("fHistMaxL1ThrCent","fHistMaxL1ThrCent", 100, 0, 100, 250, 0, 250);
 203       fHistMaxL1ThrCent->GetXaxis()->SetTitle("Centrality [%]");
 204       fHistMaxL1ThrCent->GetYaxis()->SetTitle("Maximum L1 threshold");
 205       fOutput->Add(fHistMaxL1ThrCent);
 206     }
 207   }
 208
 209   fHistTracksPt = new TH1F("fHistTracksPt","p_{T} spectrum of reconstructed tracks", fNbins, fMinBinPt, fMaxBinPt);
 210   fHistTracksPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 211   fHistTracksPt->GetYaxis()->SetTitle("counts");
 212   fOutput->Add(fHistTracksPt);
 213
 214   fHistTrPhiEta = new TH2F("fHistTrPhiEta","Phi-Eta distribution of tracks", 80, -2, 2, 128, 0, 6.4);
 215   fHistTrPhiEta->GetXaxis()->SetTitle("#eta");
 216   fHistTrPhiEta->GetYaxis()->SetTitle("#phi");
 217   fOutput->Add(fHistTrPhiEta);
 218
 219   fHistTrEmcPhiEta = new TH2F("fHistTrEmcPhiEta","Phi-Eta emcal distribution of tracks", 80, -2, 2, 128, 0, 6.4);
 220   fHistTrEmcPhiEta->GetXaxis()->SetTitle("#eta");
 221   fHistTrEmcPhiEta->GetYaxis()->SetTitle("#phi");
 222   fOutput->Add(fHistTrEmcPhiEta);
 223
 224   fHistTrPhiEtaNonProp = new TH2F("fHistTrPhiEtaNonProp","fHistTrPhiEtaNonProp", 80, -2, 2, 128, 0, 6.4);
 225   fHistTrPhiEtaNonProp->GetXaxis()->SetTitle("#eta");
 226   fHistTrPhiEtaNonProp->GetYaxis()->SetTitle("#phi");
 227   fOutput->Add(fHistTrPhiEtaNonProp);
 228
 229   fHistDeltaEtaPt = new TH2F("fHistDeltaEtaPt","fHistDeltaEtaPt", fNbins, fMinBinPt, fMaxBinPt, 80, -0.5, 0.5);
 230   fHistDeltaEtaPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 231   fHistDeltaEtaPt->GetYaxis()->SetTitle("#delta#eta");
 232   fOutput->Add(fHistDeltaEtaPt);
 233
 234   fHistDeltaPhiPt = new TH2F("fHistDeltaPhiPt","fHistDeltaPhiPt", fNbins, fMinBinPt, fMaxBinPt, 256, -1.6, 4.8);
 235   fHistDeltaPhiPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 236   fHistDeltaPhiPt->GetYaxis()->SetTitle("#delta#phi");
 237   fOutput->Add(fHistDeltaPhiPt);
 238
 239   if (fRepropagateTracks) {
 240     fHistDeltaEtaNewProp = new TH1F("fHistDeltaEtaNewProp","fHistDeltaEtaNewProp", 800, -0.4, 0.4);
 241     fHistDeltaEtaNewProp->GetXaxis()->SetTitle("#delta#eta");
 242     fHistDeltaEtaNewProp->GetYaxis()->SetTitle("counts");
 243     fOutput->Add(fHistDeltaEtaNewProp);
 244
 245     fHistDeltaPhiNewProp = new TH1F("fHistDeltaPhiNewProp","fHistDeltaPhiNewProp", 2560, -1.6, 3.2);
 246     fHistDeltaPhiNewProp->GetXaxis()->SetTitle("#delta#phi");
 247     fHistDeltaPhiNewProp->GetYaxis()->SetTitle("counts");
 248     fOutput->Add(fHistDeltaPhiNewProp);
 249   }
 250
 251   if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {
 252     fHistClusPhiEtaEnergy = new TH3F("fHistClusPhiEtaEnergy","Phi-Eta-Energy distribution of clusters", fNbins, fMinBinPt, fMaxBinPt, 80, -2, 2, 128, 0, 6.4);
 253     fHistClusPhiEtaEnergy->GetXaxis()->SetTitle("E [GeV]");
 254     fHistClusPhiEtaEnergy->GetYaxis()->SetTitle("#eta");
 255     fHistClusPhiEtaEnergy->GetZaxis()->SetTitle("#phi");
 256     fOutput->Add(fHistClusPhiEtaEnergy);
 257
 258     fHistNCellsEnergy = new TH2F("fHistNCellsEnergy","Number of cells vs. energy of clusters", fNbins, fMinBinPt, fMaxBinPt, 30, 0, 30);
 259     fHistNCellsEnergy->GetXaxis()->SetTitle("E [GeV]");
 260     fHistNCellsEnergy->GetYaxis()->SetTitle("N_{cells}");
 261     fOutput->Add(fHistNCellsEnergy);
 262   }
 263
 264   if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {
 265
 266     fHistCellsEnergy = new TH1F("fHistCellsEnergy","Energy spectrum of cells", fNbins, fMinBinPt, fMaxBinPt);
 267     fHistCellsEnergy->GetXaxis()->SetTitle("E [GeV]");
 268     fHistCellsEnergy->GetYaxis()->SetTitle("counts");
 269     fOutput->Add(fHistCellsEnergy);
 270
 271     fHistChVSneCells = new TH2F("fHistChVSneCells","Charged energy vs. neutral (cells) energy",
 272                                 (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
 273     fHistChVSneCells->GetXaxis()->SetTitle("E [GeV]");
 274     fHistChVSneCells->GetYaxis()->SetTitle("P [GeV/c]");
 275     fOutput->Add(fHistChVSneCells);
 276
 277     fHistChVSneClus = new TH2F("fHistChVSneClus","Charged energy vs. neutral (clusters) energy",
 278                                (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
 279     fHistChVSneClus->GetXaxis()->SetTitle("E [GeV]");
 280     fHistChVSneClus->GetYaxis()->SetTitle("P [GeV/c]");
 281     fOutput->Add(fHistChVSneClus);
 282
 283     fHistChVSneCorrCells = new TH2F("fHistChVSneCorrCells","Charged energy vs. neutral (corrected cells) energy",
 284                                     (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, (Int_t)(fNbins * 2.5), fMinBinPt , fMaxBinPt * 2.5);
 285     fHistChVSneCorrCells->GetXaxis()->SetTitle("E [GeV]");
 286     fHistChVSneCorrCells->GetYaxis()->SetTitle("P [GeV/c]");
 287     fOutput->Add(fHistChVSneCorrCells);
 288   }
 289
 290   for (Int_t i = 0; i < 5; i++) {
 291     TString histnamephi("fHistTrackPhi_");
 292     histnamephi += i;
 293     fHistTrackPhi[i] = new TH1F(histnamephi.Data(),histnamephi.Data(), 128, 0, 6.4);
 294     fHistTrackPhi[i]->GetXaxis()->SetTitle("Phi");
 295     fOutput->Add(fHistTrackPhi[i]);
 296
 297     TString histnameeta("fHistTrackEta_");
 298     histnameeta += i;
 299     fHistTrackEta[i] = new TH1F(histnameeta.Data(),histnameeta.Data(), 100, -2, 2);
 300     fHistTrackEta[i]->GetXaxis()->SetTitle("Eta");
 301     fOutput->Add(fHistTrackEta[i]);
 302   }
 303
 304   fHistTrackPhi[0]->SetLineColor(kRed);
 305   fHistTrackEta[0]->SetLineColor(kRed);
 306   fHistTrackPhi[1]->SetLineColor(kBlue);
 307   fHistTrackEta[1]->SetLineColor(kBlue);
 308   fHistTrackPhi[2]->SetLineColor(kGreen);
 309   fHistTrackEta[2]->SetLineColor(kGreen);
 310   fHistTrackPhi[3]->SetLineColor(kOrange);
 311   fHistTrackEta[3]->SetLineColor(kOrange);
 312   fHistTrackPhi[4]->SetLineColor(kBlack);
 313   fHistTrackEta[4]->SetLineColor(kBlack);
 314
 315   if (!fJetsName.IsNull()) {
 316
 317     TString histname;
 318
 319     for (Int_t i = 0; i < 4; i++) {
 320       histname = "fHistJetsPhiEta_";
 321       histname += i;
 322       fHistJetsPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 80, -2, 2, 128, 0, 6.4);
 323       fHistJetsPhiEta[i]->GetXaxis()->SetTitle("#eta");
 324       fHistJetsPhiEta[i]->GetYaxis()->SetTitle("#phi");
 325       fHistJetsPhiEta[i]->GetZaxis()->SetTitle("p_{T} [GeV/c]");
 326       fOutput->Add(fHistJetsPhiEta[i]);
 327
 328       histname = "fHistJetsPtNonBias_";
 329       histname += i;
 330       fHistJetsPtNonBias[i] = new TH1F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
 331       fHistJetsPtNonBias[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 332       fHistJetsPtNonBias[i]->GetYaxis()->SetTitle("counts");
 333       fOutput->Add(fHistJetsPtNonBias[i]);
 334
 335       histname = "fHistJetsPtTrack_";
 336       histname += i;
 337       fHistJetsPtTrack[i] = new TH1F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
 338       fHistJetsPtTrack[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 339       fHistJetsPtTrack[i]->GetYaxis()->SetTitle("counts");
 340       fOutput->Add(fHistJetsPtTrack[i]);
 341
 342       if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {
 343         histname = "fHistJetsPtClus_";
 344         histname += i;
 345         fHistJetsPtClus[i] = new TH1F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
 346         fHistJetsPtClus[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 347         fHistJetsPtClus[i]->GetYaxis()->SetTitle("counts");
 348         fOutput->Add(fHistJetsPtClus[i]);
 349       }
 350
 351       histname = "fHistJetsPt_";
 352       histname += i;
 353       fHistJetsPt[i] = new TH1F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
 354       fHistJetsPt[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 355       fHistJetsPt[i]->GetYaxis()->SetTitle("counts");
 356       fOutput->Add(fHistJetsPt[i]);
 357
 358       histname = "fHistJetsPtAreaNonBias_";
 359       histname += i;
 360       fHistJetsPtAreaNonBias[i] = new TH2F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, 20, 0, fJetRadius * fJetRadius * TMath::Pi() * 1.5);
 361       fHistJetsPtAreaNonBias[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 362       fHistJetsPtAreaNonBias[i]->GetYaxis()->SetTitle("area");
 363       fOutput->Add(fHistJetsPtAreaNonBias[i]);
 364
 365       histname = "fHistJetsPtArea_";
 366       histname += i;
 367       fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, 20, 0, fJetRadius * fJetRadius * TMath::Pi() * 1.5);
 368       fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 369       fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
 370       fOutput->Add(fHistJetsPtArea[i]);
 371     }
 372   }
 373
 374   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
 375 }
 376
 377 //________________________________________________________________________
 378 Bool_t AliAnalysisTaskSAQA::RetrieveEventObjects()
 379 {
 380   // Retrieve event objects.
 381
 382   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 383     return kFALSE;
 384
 385   if (!fTrgClusName.IsNull() && fDoTrigger && !fTrgClusters) {
 386     fTrgClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrgClusName));
 387     if (!fTrgClusters) {
 388       AliError(Form("%s: Could not retrieve trigger clusters %s!", GetName(), fTrgClusName.Data()));
 389       return kFALSE;
 390     }
 391     else {
 392       TClass *cl = fTrgClusters->GetClass();
 393       if (!cl->GetBaseClass("AliVCluster") && !cl->GetBaseClass("AliEmcalParticle")) {
 394         AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fTrgClusName.Data()));
 395         fTrgClusters = 0;
 396         return kFALSE;
 397       }
 398     }
 399   }
 400
 401   return kTRUE;
 402 }
 403
 404 //________________________________________________________________________
 405 Bool_t AliAnalysisTaskSAQA::FillHistograms()
 406 {
 407   // Fill histograms.
 408
 409   fHistCentrality->Fill(fCent);
 410   if (fTracks)
 411     fHistTracksCent->Fill(fCent, fTracks->GetEntriesFast());
 412   if (fCaloClusters)
 413     fHistClusCent->Fill(fCent, fCaloClusters->GetEntriesFast());
 414
 415   fHistZVertex->Fill(fVertex[2]);
 416
 417   Float_t trackSum = DoTrackLoop();
 418
 419   DoJetLoop();
 420
 421   if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {
 422
 423     if (fTracks && fCaloClusters)
 424       fHistClusTracks->Fill(fTracks->GetEntriesFast(), fCaloClusters->GetEntriesFast());
 425
 426     Float_t clusSum = DoClusterLoop();
 427
 428     Float_t cellSum = 0, cellCutSum = 0;
 429
 430     Int_t ncells = DoCellLoop(cellSum, cellCutSum);
 431
 432     if (fTracks)
 433       fHistCellsTracks->Fill(fTracks->GetEntriesFast(), ncells);
 434
 435     fHistCellsCent->Fill(fCent, ncells);
 436
 437     fHistChVSneCells->Fill(cellSum, trackSum);
 438     fHistChVSneClus->Fill(clusSum, trackSum);
 439     fHistChVSneCorrCells->Fill(cellCutSum, trackSum);
 440
 441     if (fDoTrigger) {
 442       Float_t maxTrgClus = DoTriggerClusLoop();
 443       fHistMaxL1ClusCent->Fill(fCent, maxTrgClus);
 444
 445       Int_t maxL1amp = -1;
 446       Int_t maxL1thr = -1;
 447
 448       DoTriggerPrimitives(maxL1amp, maxL1thr);
 449
 450       if (maxL1amp > -1)
 451         fHistMaxL1FastORCent->Fill(fCent, maxL1amp);
 452
 453       if (maxL1thr > -1)
 454         fHistMaxL1ThrCent->Fill(fCent, maxL1thr);
 455     }
 456   }
 457
 458   return kTRUE;
 459 }
 460
 461 //________________________________________________________________________
 462 Int_t AliAnalysisTaskSAQA::DoCellLoop(Float_t &sum, Float_t &sum_cut)
 463 {
 464   // Do cell loop.
 465
 466   AliVCaloCells *cells = InputEvent()->GetEMCALCells();
 467
 468   if (!cells)
 469     return 0;
 470
 471   const Int_t ncells = cells->GetNumberOfCells();
 472
 473   for (Int_t pos = 0; pos < ncells; pos++) {
 474     Float_t amp = cells->GetAmplitude(pos);
 475     fHistCellsEnergy->Fill(amp);
 476     sum += amp;
 477     if (amp < fCellEnergyCut)
 478       continue;
 479     sum_cut += amp;
 480   }
 481
 482   return ncells;
 483 }
 484
 485 //________________________________________________________________________
 486 Float_t AliAnalysisTaskSAQA::DoClusterLoop()
 487 {
 488   // Do cluster loop.
 489
 490   if (!fCaloClusters)
 491     return 0;
 492
 493   Float_t sum = 0;
 494
 495   // Cluster loop
 496   Int_t nclusters =  fCaloClusters->GetEntriesFast();
 497
 498   for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
 499     AliVCluster* cluster = static_cast<AliVCluster*>(fCaloClusters->At(iClusters));
 500     if (!cluster) {
 501       AliError(Form("Could not receive cluster %d", iClusters));
 502       continue;
 503     }
 504
 505     if (!AcceptCluster(cluster, kTRUE))
 506       continue;
 507
 508     sum += cluster->E();
 509
 510     TLorentzVector nPart;
 511     cluster->GetMomentum(nPart, fVertex);
 512
 513     fHistClusPhiEtaEnergy->Fill(cluster->E(), nPart.Eta(), nPart.Phi());
 514     fHistNCellsEnergy->Fill(cluster->E(), cluster->GetNCells());
 515
 516     fHistClusTimeEnergy->Fill(cluster->E(), cluster->GetTOF());
 517   }
 518
 519   return sum;
 520 }
 521
 522 //________________________________________________________________________
 523 Float_t AliAnalysisTaskSAQA::DoTrackLoop()
 524 {
 525   // Do track loop.
 526
 527   if (!fTracks)
 528     return 0;
 529
 530   Float_t sum = 0;
 531
 532   Int_t ntracks = fTracks->GetEntriesFast();
 533   Int_t nclusters = 0;
 534   if (fCaloClusters)
 535     nclusters = fCaloClusters->GetEntriesFast();
 536
 537   for (Int_t i = 0; i < ntracks; i++) {
 538
 539     AliVParticle* track = static_cast<AliVParticle*>(fTracks->At(i)); // pointer to reconstructed to track
 540
 541     if (!track) {
 542       AliError(Form("Could not retrieve track %d",i));
 543       continue;
 544     }
 545
 546     AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track);
 547
 548     if (vtrack && !AcceptTrack(vtrack, kTRUE))
 549       continue;
 550
 551     fHistTracksPt->Fill(track->Pt());
 552
 553     sum += track->P();
 554
 555     Int_t label = track->GetLabel();
 556
 557     fHistTrPhiEta->Fill(track->Eta(), track->Phi());
 558
 559     fHistTrackEta[4]->Fill(track->Eta());
 560     fHistTrackPhi[4]->Fill(track->Phi());
 561
 562     if (label >= 0 && label < 4) {
 563       fHistTrackEta[label]->Fill(track->Eta());
 564       fHistTrackPhi[label]->Fill(track->Phi());
 565     }
 566
 567     if (!vtrack)
 568       continue;
 569
 570     if (vtrack->GetTrackEtaOnEMCal() == -999 || vtrack->GetTrackPhiOnEMCal() == -999)
 571       fHistTrPhiEtaNonProp->Fill(vtrack->Eta(), vtrack->Phi());
 572
 573     fHistTrEmcPhiEta->Fill(vtrack->GetTrackEtaOnEMCal(), vtrack->GetTrackPhiOnEMCal());
 574     fHistDeltaEtaPt->Fill(vtrack->Pt(), vtrack->Eta() - vtrack->GetTrackEtaOnEMCal());
 575     fHistDeltaPhiPt->Fill(vtrack->Pt(), vtrack->Phi() - vtrack->GetTrackPhiOnEMCal());
 576
 577     if (fRepropagateTracks && vtrack->GetTrackEtaOnEMCal() > -2) {
 578       Float_t propeta = -999, propphi = -999;
 579       PropagateTrack(vtrack, propeta, propphi);
 580       fHistDeltaEtaNewProp->Fill(propeta - vtrack->GetTrackEtaOnEMCal());
 581       fHistDeltaPhiNewProp->Fill(propphi - vtrack->GetTrackPhiOnEMCal());
 582     }
 583   }
 584
 585   return sum;
 586 }
 587
 588 //____________________________________________________________________________
 589 void AliAnalysisTaskSAQA::PropagateTrack(AliVTrack *track, Float_t &eta, Float_t &phi)
 590 {
 591   eta = -999;
 592   phi = -999;
 593
 594   if (!track)
 595     return;
 596
 597   if (track->Pt() == 0)
 598     return;
 599
 600   // init the magnetic field if not already on
 601   if(!TGeoGlobalMagField::Instance()->GetField()) {
 602     AliInfo("Init the magnetic field\n");
 603     AliAODEvent* aodevent = dynamic_cast<AliAODEvent*>(InputEvent());
 604     if (aodevent) {
 605       Double_t curSol = 30000*aodevent->GetMagneticField()/5.00668;
 606       Double_t curDip = 6000 *aodevent->GetMuonMagFieldScale();
 607       AliMagF *field  = AliMagF::CreateFieldMap(curSol,curDip);
 608       TGeoGlobalMagField::Instance()->SetField(field);
 609     }
 610   }
 611
 612   Double_t cv[21];
 613   for (Int_t i = 0; i < 21; i++) cv[i] = 0;
 614
 615   Double_t pos[3], mom[3];
 616   track->GetXYZ(pos);
 617   track->GetPxPyPz(mom);
 618   AliExternalTrackParam *trackParam = new AliExternalTrackParam(pos, mom, cv, track->Charge());
 619
 620   if(!AliTrackerBase::PropagateTrackToBxByBz(trackParam, 430., 0.139, 20, kTRUE, 0.8, -1)) return;
 621   Double_t trkPos[3] = {0., 0., 0.};
 622   if(!trackParam->GetXYZ(trkPos)) return;
 623   TVector3 trkPosVec(trkPos[0], trkPos[1], trkPos[2]);
 624   eta = trkPosVec.Eta();
 625   phi = trkPosVec.Phi();
 626   if(phi < 0)
 627     phi += 2 * TMath::Pi();
 628 }
 629
 630 //________________________________________________________________________
 631 void AliAnalysisTaskSAQA::DoJetLoop()
 632 {
 633   // Do jet loop.
 634
 635   if (!fJets)
 636     return;
 637
 638   Int_t njets = fJets->GetEntriesFast();
 639
 640   for (Int_t ij = 0; ij < njets; ij++) {
 641
 642     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(ij));
 643
 644     if (!jet) {
 645       AliError(Form("Could not receive jet %d", ij));
 646       continue;
 647     }
 648
 649     if (!AcceptJet(jet, kFALSE))
 650       continue;
 651
 652     fHistJetsPtNonBias[fCentBin]->Fill(jet->Pt());
 653     fHistJetsPtAreaNonBias[fCentBin]->Fill(jet->Pt(), jet->Area());
 654
 655     if (jet->MaxTrackPt() > fPtBiasJetTrack)
 656       fHistJetsPtTrack[fCentBin]->Fill(jet->Pt());
 657
 658     if (fAnaType == kEMCAL && jet->MaxClusterPt() > fPtBiasJetClus)
 659       fHistJetsPtClus[fCentBin]->Fill(jet->Pt());
 660
 661     if (!AcceptBiasJet(jet))
 662       continue;
 663
 664     fHistJetsPt[fCentBin]->Fill(jet->Pt());
 665     fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
 666
 667     fHistJetsPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());
 668   }
 669 }
 670
 671 //________________________________________________________________________
 672 Float_t AliAnalysisTaskSAQA::DoTriggerClusLoop()
 673 {
 674   // Do trigger cluster loop.
 675
 676   if (!fTrgClusters)
 677     return 0;
 678
 679   Int_t ntrgclusters = fTrgClusters->GetEntriesFast();
 680   Float_t maxTrgClus = 0;
 681
 682   for (Int_t iClusters = 0; iClusters < ntrgclusters; iClusters++) {
 683     AliVCluster* cluster = static_cast<AliVCluster*>(fTrgClusters->At(iClusters));
 684     if (!cluster) {
 685       AliError(Form("Could not receive cluster %d", iClusters));
 686       continue;
 687     }
 688
 689     if (!(cluster->IsEMCAL())) continue;
 690
 691     if (cluster->E() > maxTrgClus)
 692       maxTrgClus = cluster->E();
 693
 694   }
 695   return maxTrgClus;
 696 }
 697
 698 //________________________________________________________________________
 699 void AliAnalysisTaskSAQA::DoTriggerPrimitives(Int_t &maxL1amp, Int_t &maxL1thr)
 700 {
 701   // Do trigger primitives loop.
 702
 703   AliVCaloTrigger *triggers = InputEvent()->GetCaloTrigger("EMCAL");
 704
 705   if (!triggers || triggers->GetEntries() == 0)
 706     return;
 707
 708   triggers->Reset();
 709   Int_t L1amp = 0;
 710   Int_t L1thr = 0;
 711   maxL1amp = -1;
 712   maxL1thr = -1;
 713
 714   while (triggers->Next()) {
 715     triggers->GetL1TimeSum(L1amp);
 716     if (maxL1amp < L1amp)
 717       maxL1amp = L1amp;
 718
 719     triggers->GetL1Threshold(L1thr);
 720     if (maxL1thr < L1thr)
 721       maxL1thr = L1thr;
 722   }
 723 }
 724
 725 //________________________________________________________________________
 726 void AliAnalysisTaskSAQA::Terminate(Option_t *)
 727 {
 728   // Called once at the end of the analysis.
 729 }