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