[u/mrichter/AliRoot.git] / PWGGA / EMCALJetTasks / UserTasks / AliAnalysisTaskSAQA.cxx

// $Id$
//
// General QA task.
//
// Author: S.Aiola

#include <TChain.h>
#include <TClonesArray.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TH3F.h>
#include <TList.h>
#include <TLorentzVector.h>

#include "AliAnalysisManager.h"
#include "AliCentrality.h"
#include "AliVCluster.h"
#include "AliVParticle.h"
#include "AliVTrack.h"
#include "AliEmcalJet.h"
#include "AliVEventHandler.h"
#include "AliLog.h"

#include "AliAnalysisTaskSAQA.h"

ClassImp(AliAnalysisTaskSAQA)

//________________________________________________________________________
AliAnalysisTaskSAQA::AliAnalysisTaskSAQA() : 
  AliAnalysisTaskEmcalJet("AliAnalysisTaskSAQA", kTRUE),
  fCellEnergyCut(0.1),
  fDoTrigger(kFALSE),
  fTrgClusName("ClustersL1GAMMAFEE"),
  fTrgClusters(0),
  fHistCentrality(0),
  fHistTracksCent(0),
  fHistClusCent(0),
  fHistMaxL1FastORCent(0),
  fHistMaxL1ClusCent(0),
  fHistMaxL1ThrCent(0),
  fHistTracksPt(0),
  fHistTrPhiEta(0),
  fHistTrEmcPhiEta(0),
  fHistClusPhiEtaEnergy(0),
  fHistNCellsEnergy(0),
  fHistCellsEnergy(0),
  fHistChVSneCells(0),
  fHistChVSneClus(0),
  fHistChVSneCorrCells(0)
{
  // Default constructor.

  for (Int_t i = 0; i < 5; i++) {
    fHistTrackPhi[i] = 0;
    fHistTrackEta[i] = 0;
  }

  for (Int_t i = 0; i < 4; i++) {
    fHistJetsPhiEtaPt[i] = 0;
    fHistJetsPtNonBias[i] = 0;
    fHistJetsPtTrack[i] = 0;
    fHistJetsPtClus[i] = 0;
    fHistJetsPt[i] = 0;
    fHistJetsPtArea[i] = 0;
    fHistJetsPtAreaNonBias[i] = 0;
  }
}

//________________________________________________________________________
AliAnalysisTaskSAQA::AliAnalysisTaskSAQA(const char *name) : 
  AliAnalysisTaskEmcalJet(name, kTRUE),
  fCellEnergyCut(0.1),
  fDoTrigger(kFALSE),
  fTrgClusName("ClustersL1GAMMAFEE"),
  fTrgClusters(0),
  fHistCentrality(0),
  fHistTracksCent(0),
  fHistClusCent(0),
  fHistMaxL1FastORCent(0),
  fHistMaxL1ClusCent(0),
  fHistMaxL1ThrCent(0),
  fHistTracksPt(0),
  fHistTrPhiEta(0),
  fHistTrEmcPhiEta(0),
  fHistClusPhiEtaEnergy(0),
  fHistNCellsEnergy(0),
  fHistCellsEnergy(0),
  fHistChVSneCells(0),
  fHistChVSneClus(0),
  fHistChVSneCorrCells(0)
{
  // Standard constructor.

  for (Int_t i = 0; i < 5; i++) {
    fHistTrackPhi[i] = 0;
    fHistTrackEta[i] = 0;
  }

  for (Int_t i = 0; i < 4; i++) {
    fHistJetsPhiEtaPt[i] = 0;
    fHistJetsPtNonBias[i] = 0;
    fHistJetsPtTrack[i] = 0;
    fHistJetsPtClus[i] = 0;
    fHistJetsPt[i] = 0;
    fHistJetsPtArea[i] = 0;
    fHistJetsPtAreaNonBias[i] = 0;
  }
}

//________________________________________________________________________
AliAnalysisTaskSAQA::~AliAnalysisTaskSAQA()
{
  // Destructor
}

//________________________________________________________________________
void AliAnalysisTaskSAQA::UserCreateOutputObjects()
{
  // Create histograms

  OpenFile(1);
  fOutput = new TList();
  fOutput->SetOwner();

  fHistCentrality = new TH1F("fHistCentrality","Event centrality distribution", fNbins, 0, 100);
  fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
  fHistCentrality->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistCentrality);

  fHistTracksCent = new TH2F("fHistTracksCent","Tracks vs. centrality", fNbins, 0, 100, fNbins, 0, 4000);
  fHistTracksCent->GetXaxis()->SetTitle("Centrality (%)");
  fHistTracksCent->GetYaxis()->SetTitle("No. of tracks");
  fOutput->Add(fHistTracksCent);

  if (fAnaType == kEMCAL) {
    fHistClusCent = new TH2F("fHistClusCent","Clusters vs. centrality", fNbins, 0, 100, fNbins, 0, 2000);
    fHistClusCent->GetXaxis()->SetTitle("Centrality (%)");
    fHistClusCent->GetYaxis()->SetTitle("No. of clusters");
    fOutput->Add(fHistClusCent);
    
    if (fDoTrigger) {
      fHistMaxL1FastORCent = new TH2F("fHistMaxL1FastORCent","fHistMaxL1ClusCent", 100, 0, 100, 250, 0, 250);
      fHistMaxL1FastORCent->GetXaxis()->SetTitle("Centrality [%]");
      fHistMaxL1FastORCent->GetYaxis()->SetTitle("Maximum L1 FastOR");
      fOutput->Add(fHistMaxL1FastORCent);
      
      fHistMaxL1ClusCent = new TH2F("fHistMaxL1ClusCent","fHistMaxL1ClusCent", 100, 0, 100, 250, 0, 250);
      fHistMaxL1ClusCent->GetXaxis()->SetTitle("Centrality [%]");
      fHistMaxL1ClusCent->GetYaxis()->SetTitle("Maximum L1 trigger cluster");
      fOutput->Add(fHistMaxL1ClusCent);
      
      fHistMaxL1ThrCent = new TH2F("fHistMaxL1ThrCent","fHistMaxL1ThrCent", 100, 0, 100, 250, 0, 250);
      fHistMaxL1ThrCent->GetXaxis()->SetTitle("Centrality [%]");
      fHistMaxL1ThrCent->GetYaxis()->SetTitle("Maximum L1 threshold");
      fOutput->Add(fHistMaxL1ThrCent);
    }
  }

  fHistTracksPt = new TH1F("fHistTracksPt","p_{T} spectrum of reconstructed tracks", fNbins, fMinBinPt, fMaxBinPt);
  fHistTracksPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
  fHistTracksPt->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistTracksPt);
       
  fHistTrPhiEta = new TH2F("fHistTrPhiEta","Phi-Eta distribution of tracks", 80, -2, 2, 128, 0, 6.4);
  fHistTrPhiEta->GetXaxis()->SetTitle("#eta");
  fHistTrPhiEta->GetYaxis()->SetTitle("#phi");
  fOutput->Add(fHistTrPhiEta);

  fHistTrEmcPhiEta = new TH2F("fHistTrEmcPhiEta","Phi-Eta emcal distribution of tracks", 80, -2, 2, 128, 0, 6.4);
  fHistTrEmcPhiEta->GetXaxis()->SetTitle("#eta");
  fHistTrEmcPhiEta->GetYaxis()->SetTitle("#phi");
  fOutput->Add(fHistTrEmcPhiEta);

  if (fAnaType == kEMCAL) {
    fHistClusPhiEtaEnergy = new TH3F("fHistClusPhiEtaEnergy","Phi-Eta-Energy distribution of clusters", fNbins, fMinBinPt, fMaxBinPt, 80, -2, 2, 128, 0, 6.4);
    fHistClusPhiEtaEnergy->GetXaxis()->SetTitle("E [GeV]");
    fHistClusPhiEtaEnergy->GetYaxis()->SetTitle("#eta");
    fHistClusPhiEtaEnergy->GetZaxis()->SetTitle("#phi");
    fOutput->Add(fHistClusPhiEtaEnergy);

    fHistNCellsEnergy = new TH2F("fHistNCellsEnergy","Number of cells vs. energy of clusters", fNbins, fMinBinPt, fMaxBinPt, 30, 0, 30);
    fHistNCellsEnergy->GetXaxis()->SetTitle("E [GeV]");
    fHistNCellsEnergy->GetYaxis()->SetTitle("N_{cells}");
    fOutput->Add(fHistNCellsEnergy);
  }

  if (fAnaType == kEMCAL) {
   
    fHistCellsEnergy = new TH1F("fHistCellsEnergy","Energy spectrum of cells", fNbins, fMinBinPt, fMaxBinPt);
    fHistCellsEnergy->GetXaxis()->SetTitle("E [GeV]");
    fHistCellsEnergy->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistCellsEnergy);
    
    fHistChVSneCells = new TH2F("fHistChVSneCells","Charged energy vs. neutral (cells) energy", 
				fNbins * 2.5, fMinBinPt, fMaxBinPt * 2.5, fNbins * 2.5, fMinBinPt, fMaxBinPt * 2.5);
    fHistChVSneCells->GetXaxis()->SetTitle("E [GeV]");
    fHistChVSneCells->GetYaxis()->SetTitle("P [GeV/c]");
    fOutput->Add(fHistChVSneCells);
    
    fHistChVSneClus = new TH2F("fHistChVSneClus","Charged energy vs. neutral (clusters) energy", 
			       fNbins * 2.5, fMinBinPt, fMaxBinPt * 2.5, fNbins * 2.5, fMinBinPt, fMaxBinPt * 2.5);
    fHistChVSneClus->GetXaxis()->SetTitle("E [GeV]");
    fHistChVSneClus->GetYaxis()->SetTitle("P [GeV/c]");
    fOutput->Add(fHistChVSneClus);
    
    fHistChVSneCorrCells = new TH2F("fHistChVSneCorrCells","Charged energy vs. neutral (corrected cells) energy", 
				    fNbins * 2.5, fMinBinPt, fMaxBinPt * 2.5, fNbins * 2.5, fMinBinPt , fMaxBinPt * 2.5);
    fHistChVSneCorrCells->GetXaxis()->SetTitle("E [GeV]");
    fHistChVSneCorrCells->GetYaxis()->SetTitle("P [GeV/c]");
    fOutput->Add(fHistChVSneCorrCells);
  }
 
  for (Int_t i = 0; i < 5; i++) {
    TString histnamephi("fHistTrackPhi_");
    histnamephi += i;
    fHistTrackPhi[i] = new TH1F(histnamephi.Data(),histnamephi.Data(), 128, 0, 6.4);
    fHistTrackPhi[i]->GetXaxis()->SetTitle("Phi");
    fOutput->Add(fHistTrackPhi[i]);

    TString histnameeta("fHistTrackEta_");
    histnameeta += i;
    fHistTrackEta[i] = new TH1F(histnameeta.Data(),histnameeta.Data(), 100, -2, 2);
    fHistTrackEta[i]->GetXaxis()->SetTitle("Eta");
    fOutput->Add(fHistTrackEta[i]);
  }
  
  fHistTrackPhi[0]->SetLineColor(kRed);
  fHistTrackEta[0]->SetLineColor(kRed);
  fHistTrackPhi[1]->SetLineColor(kBlue);
  fHistTrackEta[1]->SetLineColor(kBlue);
  fHistTrackPhi[2]->SetLineColor(kGreen);
  fHistTrackEta[2]->SetLineColor(kGreen);
  fHistTrackPhi[3]->SetLineColor(kOrange);
  fHistTrackEta[3]->SetLineColor(kOrange);
  fHistTrackPhi[4]->SetLineColor(kBlack);
  fHistTrackEta[4]->SetLineColor(kBlack);

  if (!fJetsName.IsNull()) {

    TString histname;

    for (Int_t i = 0; i < 4; i++) {
      histname = "fHistJetsPhiEtaPt_";
      histname += i;
      fHistJetsPhiEtaPt[i] = new TH3F(histname.Data(), histname.Data(), 80, -2, 2, 128, 0, 6.4, fNbins * 2.5, fMinBinPt, fMaxBinPt * 2.5);
      fHistJetsPhiEtaPt[i]->GetXaxis()->SetTitle("#eta");
      fHistJetsPhiEtaPt[i]->GetYaxis()->SetTitle("#phi");
      fHistJetsPhiEtaPt[i]->GetZaxis()->SetTitle("p_{T} [GeV/c]");
      fOutput->Add(fHistJetsPhiEtaPt[i]);

      histname = "fHistJetsPtNonBias_";
      histname += i;
      fHistJetsPtNonBias[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2.5, fMinBinPt, fMaxBinPt * 2.5);
      fHistJetsPtNonBias[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
      fHistJetsPtNonBias[i]->GetYaxis()->SetTitle("counts");
      fOutput->Add(fHistJetsPtNonBias[i]);

      histname = "fHistJetsPtTrack_";
      histname += i;
      fHistJetsPtTrack[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2.5, fMinBinPt, fMaxBinPt * 2.5);
      fHistJetsPtTrack[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
      fHistJetsPtTrack[i]->GetYaxis()->SetTitle("counts");
      fOutput->Add(fHistJetsPtTrack[i]);

      if (fAnaType == kEMCAL) {
	histname = "fHistJetsPtClus_";
	histname += i;
	fHistJetsPtClus[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2.5, fMinBinPt, fMaxBinPt * 2.5);
	fHistJetsPtClus[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
	fHistJetsPtClus[i]->GetYaxis()->SetTitle("counts");
	fOutput->Add(fHistJetsPtClus[i]);
      }

      histname = "fHistJetsPt_";
      histname += i;
      fHistJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2.5, fMinBinPt, fMaxBinPt * 2.5);
      fHistJetsPt[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
      fHistJetsPt[i]->GetYaxis()->SetTitle("counts");
      fOutput->Add(fHistJetsPt[i]);

      histname = "fHistJetsPtAreaNonBias_";
      histname += i;
      fHistJetsPtAreaNonBias[i] = new TH2F(histname.Data(), histname.Data(), fNbins * 2.5, fMinBinPt, fMaxBinPt * 2.5, 20, 0, fJetRadius * fJetRadius * TMath::Pi() * 1.5);
      fHistJetsPtAreaNonBias[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
      fHistJetsPtAreaNonBias[i]->GetYaxis()->SetTitle("area");
      fOutput->Add(fHistJetsPtAreaNonBias[i]);

      histname = "fHistJetsPtArea_";
      histname += i;
      fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins * 2.5, fMinBinPt, fMaxBinPt * 2.5, 20, 0, fJetRadius * fJetRadius * TMath::Pi() * 1.5);
      fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
      fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
      fOutput->Add(fHistJetsPtArea[i]);
    }
  }

  PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
}

//________________________________________________________________________
Bool_t AliAnalysisTaskSAQA::RetrieveEventObjects()
{
  // Retrieve event objects.

  if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
    return kFALSE;

  if (!fTrgClusName.IsNull() && fDoTrigger && !fTrgClusters) {
    fTrgClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrgClusName));
    if (!fTrgClusters) {
      AliError(Form("%s: Could not retrieve trigger clusters %s!", GetName(), fTrgClusName.Data())); 
      return kFALSE;
    }
    else {
      TClass *cl = fTrgClusters->GetClass();
      if (!cl->GetBaseClass("AliVCluster") && !cl->GetBaseClass("AliEmcalParticle")) {
	AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fTrgClusName.Data())); 
	fTrgClusters = 0;
	return kFALSE;
      }
    }
  }

  return kTRUE;
}

//________________________________________________________________________
Bool_t AliAnalysisTaskSAQA::FillHistograms()
{
  // Fill histograms.

  fHistCentrality->Fill(fCent);
  if (fTracks)
    fHistTracksCent->Fill(fCent, fTracks->GetEntriesFast());
  if (fCaloClusters)
    fHistClusCent->Fill(fCent, fCaloClusters->GetEntriesFast());

  Float_t trackSum = DoTrackLoop();

  DoJetLoop();

  if (fAnaType == kEMCAL) {
    Float_t clusSum = DoClusterLoop();

    Float_t cellSum = 0, cellCutSum = 0;
    
    DoCellLoop(cellSum, cellCutSum);
    
    fHistChVSneCells->Fill(cellSum, trackSum);
    fHistChVSneClus->Fill(clusSum, trackSum);
    fHistChVSneCorrCells->Fill(cellCutSum, trackSum);

    if (fDoTrigger) {
      Float_t maxTrgClus = DoTriggerClusLoop();
      fHistMaxL1ClusCent->Fill(fCent, maxTrgClus);
    
      Int_t maxL1amp = -1;
      Int_t maxL1thr = -1;
    
      DoTriggerPrimitives(maxL1amp, maxL1thr);
      
      if (maxL1amp > -1) 
	fHistMaxL1FastORCent->Fill(fCent, maxL1amp);
      
      if (maxL1thr > -1) 
	fHistMaxL1ThrCent->Fill(fCent, maxL1thr);
    }
  }

  return kTRUE;
}

//________________________________________________________________________
void AliAnalysisTaskSAQA::DoCellLoop(Float_t &sum, Float_t &sum_cut)
{
  // Do cell loop.

  AliVCaloCells *cells = InputEvent()->GetEMCALCells();

  if (!cells)
    return;

  const Int_t ncells = cells->GetNumberOfCells();

  for (Int_t pos = 0; pos < ncells; pos++) {
    Float_t amp = cells->GetAmplitude(pos);
    fHistCellsEnergy->Fill(amp);
    sum += amp;
    if (amp < fCellEnergyCut)
      continue;
    sum_cut += amp;
  } 
}

//________________________________________________________________________
Float_t AliAnalysisTaskSAQA::DoClusterLoop()
{
  // Do cluster loop.

  if (!fCaloClusters)
    return 0;

  Float_t sum = 0;

  // Cluster loop
  Int_t nclusters =  fCaloClusters->GetEntriesFast();

  for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
    AliVCluster* cluster = static_cast<AliVCluster*>(fCaloClusters->At(iClusters));
    if (!cluster) {
      AliError(Form("Could not receive cluster %d", iClusters));
      continue;
    }  

    if (!AcceptCluster(cluster, kTRUE))
      continue;

    sum += cluster->E();

    TLorentzVector nPart;
    cluster->GetMomentum(nPart, fVertex);

    fHistClusPhiEtaEnergy->Fill(cluster->E(), nPart.Eta(), nPart.Phi());
    fHistNCellsEnergy->Fill(cluster->E(), cluster->GetNCells());
  }

  return sum;
}

//________________________________________________________________________
Float_t AliAnalysisTaskSAQA::DoTrackLoop()
{
  // Do track loop.

  if (!fTracks)
    return 0;

  Float_t sum = 0;

  Int_t ntracks = fTracks->GetEntriesFast();
  Int_t nclusters = 0;
  if (fCaloClusters)
    nclusters = fCaloClusters->GetEntriesFast();

  for (Int_t i = 0; i < ntracks; i++) {

    AliVParticle* track = static_cast<AliVParticle*>(fTracks->At(i)); // pointer to reconstructed to track  

    if (!track) {
      AliError(Form("Could not retrieve track %d",i)); 
      continue; 
    }

    AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track); 
    
    if (vtrack && !AcceptTrack(vtrack, kTRUE)) 
      continue;
    
    fHistTracksPt->Fill(track->Pt());

    sum += track->P();

    Int_t label = track->GetLabel();
      
    fHistTrPhiEta->Fill(track->Eta(), track->Phi());
    
    fHistTrackEta[4]->Fill(track->Eta());
    fHistTrackPhi[4]->Fill(track->Phi());

    if (label >= 0 && label < 4) {
      fHistTrackEta[label]->Fill(track->Eta());
      fHistTrackPhi[label]->Fill(track->Phi());
    }

    if (!vtrack)
      continue;

    fHistTrEmcPhiEta->Fill(vtrack->GetTrackEtaOnEMCal(), vtrack->GetTrackPhiOnEMCal());
  }
  
  return sum;
}

//________________________________________________________________________
void AliAnalysisTaskSAQA::DoJetLoop()
{
  // Do jet loop.

  if (!fJets)
    return;

  Int_t njets = fJets->GetEntriesFast();

  for (Int_t ij = 0; ij < njets; ij++) {

    AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(ij));

    if (!jet) {
      AliError(Form("Could not receive jet %d", ij));
      continue;
    }  

    if (!AcceptJet(jet, kFALSE))
      continue;

    fHistJetsPtNonBias[fCentBin]->Fill(jet->Pt());
    fHistJetsPtAreaNonBias[fCentBin]->Fill(jet->Pt(), jet->Area());

    if (jet->MaxTrackPt() > fPtBiasJetTrack)
      fHistJetsPtTrack[fCentBin]->Fill(jet->Pt());
    
    if (fAnaType == kEMCAL && jet->MaxClusterPt() > fPtBiasJetClus)
      fHistJetsPtClus[fCentBin]->Fill(jet->Pt());
    
    if (!AcceptBiasJet(jet))
      continue;

    fHistJetsPt[fCentBin]->Fill(jet->Pt());
    fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());

    fHistJetsPhiEtaPt[fCentBin]->Fill(jet->Eta(), jet->Phi(), jet->Pt());
  }
}

//________________________________________________________________________
Float_t AliAnalysisTaskSAQA::DoTriggerClusLoop()
{
  // Do trigger cluster loop.

  if (!fTrgClusters)
    return 0;

  Int_t ntrgclusters = fTrgClusters->GetEntriesFast();
  Float_t maxTrgClus = 0;

  for (Int_t iClusters = 0; iClusters < ntrgclusters; iClusters++) {
    AliVCluster* cluster = static_cast<AliVCluster*>(fTrgClusters->At(iClusters));
    if (!cluster) {
      AliError(Form("Could not receive cluster %d", iClusters));
      continue;
    }  
    
    if (!(cluster->IsEMCAL())) continue;

    if (cluster->E() > maxTrgClus)
      maxTrgClus = cluster->E();

  }
  return maxTrgClus;
}

//________________________________________________________________________
void AliAnalysisTaskSAQA::DoTriggerPrimitives(Int_t &maxL1amp, Int_t &maxL1thr)
{
  // Do trigger primitives loop.

  AliVCaloTrigger *triggers = InputEvent()->GetCaloTrigger("EMCAL");

  if (!triggers || triggers->GetEntries() == 0)
    return;
    
  triggers->Reset();
  Int_t L1amp = 0;
  Int_t L1thr = 0;
  maxL1amp = -1;
  maxL1thr = -1;

  while (triggers->Next()) {
    triggers->GetL1TimeSum(L1amp);
    if (maxL1amp < L1amp) 
      maxL1amp = L1amp;

    triggers->GetL1Threshold(L1thr);
    if (maxL1thr < L1thr) 
      maxL1thr = L1thr;
  }
}

//________________________________________________________________________
void AliAnalysisTaskSAQA::Terminate(Option_t *) 
{
  // Called once at the end of the analysis.
}
Commit	Line	Data
6e8d91c9	1	// $Id$
c3ba2d3d	2	//
297edd60	3	// General QA task.
c3ba2d3d	4	//
297edd60	5	// Author: S.Aiola
c3ba2d3d	6
	7	#include <TChain.h>
	8	#include <TClonesArray.h>
	9	#include <TH1F.h>
	10	#include <TH2F.h>
3000c095	11	#include <TH3F.h>
c3ba2d3d	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 "AliLog.h"
	23
	24	#include "AliAnalysisTaskSAQA.h"
	25
	26	ClassImp(AliAnalysisTaskSAQA)
	27
	28	//________________________________________________________________________
	29	AliAnalysisTaskSAQA::AliAnalysisTaskSAQA() :
99cfd012	30	AliAnalysisTaskEmcalJet("AliAnalysisTaskSAQA", kTRUE),
c3ba2d3d	31	fCellEnergyCut(0.1),
1f6fff78	32	fDoTrigger(kFALSE),
c3ba2d3d	33	fTrgClusName("ClustersL1GAMMAFEE"),
	34	fTrgClusters(0),
	35	fHistCentrality(0),
	36	fHistTracksCent(0),
	37	fHistClusCent(0),
	38	fHistMaxL1FastORCent(0),
	39	fHistMaxL1ClusCent(0),
	40	fHistMaxL1ThrCent(0),
	41	fHistTracksPt(0),
	42	fHistTrPhiEta(0),
6fd5039f	43	fHistTrEmcPhiEta(0),
7cf4626b	44	fHistClusPhiEtaEnergy(0),
7cf4626b	45	fHistNCellsEnergy(0),
c3ba2d3d	46	fHistCellsEnergy(0),
	47	fHistChVSneCells(0),
	48	fHistChVSneClus(0),
	49	fHistChVSneCorrCells(0)
	50	{
	51	// Default constructor.
	52
	53	for (Int_t i = 0; i < 5; i++) {
	54	fHistTrackPhi[i] = 0;
	55	fHistTrackEta[i] = 0;
	56	}
6fd5039f	57
6fd5039f	58	for (Int_t i = 0; i < 4; i++) {
3000c095	59	fHistJetsPhiEtaPt[i] = 0;
99cfd012	60	fHistJetsPtNonBias[i] = 0;
6fd5039f	61	fHistJetsPtTrack[i] = 0;
	62	fHistJetsPtClus[i] = 0;
	63	fHistJetsPt[i] = 0;
98750b70	64	fHistJetsPtArea[i] = 0;
99cfd012	65	fHistJetsPtAreaNonBias[i] = 0;
6fd5039f	66	}
c3ba2d3d	67	}
	68
	69	//________________________________________________________________________
	70	AliAnalysisTaskSAQA::AliAnalysisTaskSAQA(const char *name) :
99cfd012	71	AliAnalysisTaskEmcalJet(name, kTRUE),
c3ba2d3d	72	fCellEnergyCut(0.1),
1f6fff78	73	fDoTrigger(kFALSE),
c3ba2d3d	74	fTrgClusName("ClustersL1GAMMAFEE"),
	75	fTrgClusters(0),
	76	fHistCentrality(0),
	77	fHistTracksCent(0),
	78	fHistClusCent(0),
	79	fHistMaxL1FastORCent(0),
	80	fHistMaxL1ClusCent(0),
	81	fHistMaxL1ThrCent(0),
	82	fHistTracksPt(0),
	83	fHistTrPhiEta(0),
6fd5039f	84	fHistTrEmcPhiEta(0),
7cf4626b	85	fHistClusPhiEtaEnergy(0),
7cf4626b	86	fHistNCellsEnergy(0),
c3ba2d3d	87	fHistCellsEnergy(0),
	88	fHistChVSneCells(0),
	89	fHistChVSneClus(0),
	90	fHistChVSneCorrCells(0)
	91	{
	92	// Standard constructor.
	93
	94	for (Int_t i = 0; i < 5; i++) {
	95	fHistTrackPhi[i] = 0;
	96	fHistTrackEta[i] = 0;
	97	}
6fd5039f	98
6fd5039f	99	for (Int_t i = 0; i < 4; i++) {
3000c095	100	fHistJetsPhiEtaPt[i] = 0;
99cfd012	101	fHistJetsPtNonBias[i] = 0;
6fd5039f	102	fHistJetsPtTrack[i] = 0;
	103	fHistJetsPtClus[i] = 0;
	104	fHistJetsPt[i] = 0;
98750b70	105	fHistJetsPtArea[i] = 0;
99cfd012	106	fHistJetsPtAreaNonBias[i] = 0;
6fd5039f	107	}
c3ba2d3d	108	}
	109
	110	//________________________________________________________________________
	111	AliAnalysisTaskSAQA::~AliAnalysisTaskSAQA()
	112	{
	113	// Destructor
	114	}
	115
	116	//________________________________________________________________________
	117	void AliAnalysisTaskSAQA::UserCreateOutputObjects()
	118	{
	119	// Create histograms
	120
	121	OpenFile(1);
	122	fOutput = new TList();
6fd5039f	123	fOutput->SetOwner();
c3ba2d3d	124
	125	fHistCentrality = new TH1F("fHistCentrality","Event centrality distribution", fNbins, 0, 100);
	126	fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
	127	fHistCentrality->GetYaxis()->SetTitle("counts");
	128	fOutput->Add(fHistCentrality);
	129
	130	fHistTracksCent = new TH2F("fHistTracksCent","Tracks vs. centrality", fNbins, 0, 100, fNbins, 0, 4000);
	131	fHistTracksCent->GetXaxis()->SetTitle("Centrality (%)");
	132	fHistTracksCent->GetYaxis()->SetTitle("No. of tracks");
	133	fOutput->Add(fHistTracksCent);
	134
6fd5039f	135	if (fAnaType == kEMCAL) {
	136	fHistClusCent = new TH2F("fHistClusCent","Clusters vs. centrality", fNbins, 0, 100, fNbins, 0, 2000);
	137	fHistClusCent->GetXaxis()->SetTitle("Centrality (%)");
	138	fHistClusCent->GetYaxis()->SetTitle("No. of clusters");
	139	fOutput->Add(fHistClusCent);
	140
1f6fff78	141	if (fDoTrigger) {
	142	fHistMaxL1FastORCent = new TH2F("fHistMaxL1FastORCent","fHistMaxL1ClusCent", 100, 0, 100, 250, 0, 250);
	143	fHistMaxL1FastORCent->GetXaxis()->SetTitle("Centrality [%]");
	144	fHistMaxL1FastORCent->GetYaxis()->SetTitle("Maximum L1 FastOR");
	145	fOutput->Add(fHistMaxL1FastORCent);
	146
	147	fHistMaxL1ClusCent = new TH2F("fHistMaxL1ClusCent","fHistMaxL1ClusCent", 100, 0, 100, 250, 0, 250);
	148	fHistMaxL1ClusCent->GetXaxis()->SetTitle("Centrality [%]");
	149	fHistMaxL1ClusCent->GetYaxis()->SetTitle("Maximum L1 trigger cluster");
	150	fOutput->Add(fHistMaxL1ClusCent);
	151
	152	fHistMaxL1ThrCent = new TH2F("fHistMaxL1ThrCent","fHistMaxL1ThrCent", 100, 0, 100, 250, 0, 250);
	153	fHistMaxL1ThrCent->GetXaxis()->SetTitle("Centrality [%]");
	154	fHistMaxL1ThrCent->GetYaxis()->SetTitle("Maximum L1 threshold");
	155	fOutput->Add(fHistMaxL1ThrCent);
	156	}
6fd5039f	157	}
	158
	159	fHistTracksPt = new TH1F("fHistTracksPt","p_{T} spectrum of reconstructed tracks", fNbins, fMinBinPt, fMaxBinPt);
	160	fHistTracksPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
c3ba2d3d	161	fHistTracksPt->GetYaxis()->SetTitle("counts");
	162	fOutput->Add(fHistTracksPt);
	163
6fd5039f	164	fHistTrPhiEta = new TH2F("fHistTrPhiEta","Phi-Eta distribution of tracks", 80, -2, 2, 128, 0, 6.4);
c3ba2d3d	165	fHistTrPhiEta->GetXaxis()->SetTitle("#eta");
	166	fHistTrPhiEta->GetYaxis()->SetTitle("#phi");
	167	fOutput->Add(fHistTrPhiEta);
	168
6fd5039f	169	fHistTrEmcPhiEta = new TH2F("fHistTrEmcPhiEta","Phi-Eta emcal distribution of tracks", 80, -2, 2, 128, 0, 6.4);
	170	fHistTrEmcPhiEta->GetXaxis()->SetTitle("#eta");
	171	fHistTrEmcPhiEta->GetYaxis()->SetTitle("#phi");
	172	fOutput->Add(fHistTrEmcPhiEta);
	173
	174	if (fAnaType == kEMCAL) {
7cf4626b	175	fHistClusPhiEtaEnergy = new TH3F("fHistClusPhiEtaEnergy","Phi-Eta-Energy distribution of clusters", fNbins, fMinBinPt, fMaxBinPt, 80, -2, 2, 128, 0, 6.4);
	176	fHistClusPhiEtaEnergy->GetXaxis()->SetTitle("E [GeV]");
	177	fHistClusPhiEtaEnergy->GetYaxis()->SetTitle("#eta");
	178	fHistClusPhiEtaEnergy->GetZaxis()->SetTitle("#phi");
	179	fOutput->Add(fHistClusPhiEtaEnergy);
	180
	181	fHistNCellsEnergy = new TH2F("fHistNCellsEnergy","Number of cells vs. energy of clusters", fNbins, fMinBinPt, fMaxBinPt, 30, 0, 30);
	182	fHistNCellsEnergy->GetXaxis()->SetTitle("E [GeV]");
	183	fHistNCellsEnergy->GetYaxis()->SetTitle("N_{cells}");
	184	fOutput->Add(fHistNCellsEnergy);
6fd5039f	185	}
c3ba2d3d	186
6fd5039f	187	if (fAnaType == kEMCAL) {
6fd5039f	188
	189	fHistCellsEnergy = new TH1F("fHistCellsEnergy","Energy spectrum of cells", fNbins, fMinBinPt, fMaxBinPt);
	190	fHistCellsEnergy->GetXaxis()->SetTitle("E [GeV]");
	191	fHistCellsEnergy->GetYaxis()->SetTitle("counts");
	192	fOutput->Add(fHistCellsEnergy);
	193
98750b70	194	fHistChVSneCells = new TH2F("fHistChVSneCells","Charged energy vs. neutral (cells) energy",
98750b70	195	fNbins * 2.5, fMinBinPt, fMaxBinPt * 2.5, fNbins * 2.5, fMinBinPt, fMaxBinPt * 2.5);
6fd5039f	196	fHistChVSneCells->GetXaxis()->SetTitle("E [GeV]");
	197	fHistChVSneCells->GetYaxis()->SetTitle("P [GeV/c]");
	198	fOutput->Add(fHistChVSneCells);
	199
98750b70	200	fHistChVSneClus = new TH2F("fHistChVSneClus","Charged energy vs. neutral (clusters) energy",
98750b70	201	fNbins * 2.5, fMinBinPt, fMaxBinPt * 2.5, fNbins * 2.5, fMinBinPt, fMaxBinPt * 2.5);
6fd5039f	202	fHistChVSneClus->GetXaxis()->SetTitle("E [GeV]");
	203	fHistChVSneClus->GetYaxis()->SetTitle("P [GeV/c]");
	204	fOutput->Add(fHistChVSneClus);
	205
98750b70	206	fHistChVSneCorrCells = new TH2F("fHistChVSneCorrCells","Charged energy vs. neutral (corrected cells) energy",
98750b70	207	fNbins * 2.5, fMinBinPt, fMaxBinPt * 2.5, fNbins * 2.5, fMinBinPt , fMaxBinPt * 2.5);
6fd5039f	208	fHistChVSneCorrCells->GetXaxis()->SetTitle("E [GeV]");
	209	fHistChVSneCorrCells->GetYaxis()->SetTitle("P [GeV/c]");
	210	fOutput->Add(fHistChVSneCorrCells);
c3ba2d3d	211	}
c3ba2d3d	212
	213	for (Int_t i = 0; i < 5; i++) {
	214	TString histnamephi("fHistTrackPhi_");
	215	histnamephi += i;
	216	fHistTrackPhi[i] = new TH1F(histnamephi.Data(),histnamephi.Data(), 128, 0, 6.4);
	217	fHistTrackPhi[i]->GetXaxis()->SetTitle("Phi");
	218	fOutput->Add(fHistTrackPhi[i]);
	219
	220	TString histnameeta("fHistTrackEta_");
	221	histnameeta += i;
	222	fHistTrackEta[i] = new TH1F(histnameeta.Data(),histnameeta.Data(), 100, -2, 2);
	223	fHistTrackEta[i]->GetXaxis()->SetTitle("Eta");
	224	fOutput->Add(fHistTrackEta[i]);
	225	}
	226
	227	fHistTrackPhi[0]->SetLineColor(kRed);
	228	fHistTrackEta[0]->SetLineColor(kRed);
	229	fHistTrackPhi[1]->SetLineColor(kBlue);
	230	fHistTrackEta[1]->SetLineColor(kBlue);
	231	fHistTrackPhi[2]->SetLineColor(kGreen);
	232	fHistTrackEta[2]->SetLineColor(kGreen);
	233	fHistTrackPhi[3]->SetLineColor(kOrange);
	234	fHistTrackEta[3]->SetLineColor(kOrange);
	235	fHistTrackPhi[4]->SetLineColor(kBlack);
	236	fHistTrackEta[4]->SetLineColor(kBlack);
	237
e44e8726	238	if (!fJetsName.IsNull()) {
6fd5039f	239
e44e8726	240	TString histname;
	241
	242	for (Int_t i = 0; i < 4; i++) {
3000c095	243	histname = "fHistJetsPhiEtaPt_";
	244	histname += i;
	245	fHistJetsPhiEtaPt[i] = new TH3F(histname.Data(), histname.Data(), 80, -2, 2, 128, 0, 6.4, fNbins * 2.5, fMinBinPt, fMaxBinPt * 2.5);
	246	fHistJetsPhiEtaPt[i]->GetXaxis()->SetTitle("#eta");
	247	fHistJetsPhiEtaPt[i]->GetYaxis()->SetTitle("#phi");
	248	fHistJetsPhiEtaPt[i]->GetZaxis()->SetTitle("p_{T} [GeV/c]");
	249	fOutput->Add(fHistJetsPhiEtaPt[i]);
	250
e44e8726	251	histname = "fHistJetsPtNonBias_";
98750b70	252	histname += i;
e44e8726	253	fHistJetsPtNonBias[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2.5, fMinBinPt, fMaxBinPt * 2.5);
	254	fHistJetsPtNonBias[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
	255	fHistJetsPtNonBias[i]->GetYaxis()->SetTitle("counts");
	256	fOutput->Add(fHistJetsPtNonBias[i]);
	257
	258	histname = "fHistJetsPtTrack_";
	259	histname += i;
	260	fHistJetsPtTrack[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2.5, fMinBinPt, fMaxBinPt * 2.5);
	261	fHistJetsPtTrack[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
	262	fHistJetsPtTrack[i]->GetYaxis()->SetTitle("counts");
	263	fOutput->Add(fHistJetsPtTrack[i]);
	264
	265	if (fAnaType == kEMCAL) {
	266	histname = "fHistJetsPtClus_";
	267	histname += i;
	268	fHistJetsPtClus[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2.5, fMinBinPt, fMaxBinPt * 2.5);
	269	fHistJetsPtClus[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
	270	fHistJetsPtClus[i]->GetYaxis()->SetTitle("counts");
	271	fOutput->Add(fHistJetsPtClus[i]);
	272	}
	273
	274	histname = "fHistJetsPt_";
	275	histname += i;
	276	fHistJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2.5, fMinBinPt, fMaxBinPt * 2.5);
	277	fHistJetsPt[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
	278	fHistJetsPt[i]->GetYaxis()->SetTitle("counts");
	279	fOutput->Add(fHistJetsPt[i]);
99cfd012	280
e44e8726	281	histname = "fHistJetsPtAreaNonBias_";
	282	histname += i;
	283	fHistJetsPtAreaNonBias[i] = new TH2F(histname.Data(), histname.Data(), fNbins * 2.5, fMinBinPt, fMaxBinPt * 2.5, 20, 0, fJetRadius * fJetRadius * TMath::Pi() * 1.5);
	284	fHistJetsPtAreaNonBias[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
	285	fHistJetsPtAreaNonBias[i]->GetYaxis()->SetTitle("area");
	286	fOutput->Add(fHistJetsPtAreaNonBias[i]);
	287
	288	histname = "fHistJetsPtArea_";
	289	histname += i;
	290	fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins * 2.5, fMinBinPt, fMaxBinPt * 2.5, 20, 0, fJetRadius * fJetRadius * TMath::Pi() * 1.5);
	291	fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
	292	fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
	293	fOutput->Add(fHistJetsPtArea[i]);
	294	}
6fd5039f	295	}
6fd5039f	296
c3ba2d3d	297	PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
	298	}
	299
	300	//________________________________________________________________________
6fd5039f	301	Bool_t AliAnalysisTaskSAQA::RetrieveEventObjects()
c3ba2d3d	302	{
c3604199	303	// Retrieve event objects.
	304
	305	if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
6fd5039f	306	return kFALSE;
c3ba2d3d	307
e44e8726	308	if (!fTrgClusName.IsNull() && fDoTrigger && !fTrgClusters) {
c3ba2d3d	309	fTrgClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrgClusName));
c3ba2d3d	310	if (!fTrgClusters) {
e44e8726	311	AliError(Form("%s: Could not retrieve trigger clusters %s!", GetName(), fTrgClusName.Data()));
	312	return kFALSE;
	313	}
	314	else {
	315	TClass *cl = fTrgClusters->GetClass();
	316	if (!cl->GetBaseClass("AliVCluster") && !cl->GetBaseClass("AliEmcalParticle")) {
	317	AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fTrgClusName.Data()));
	318	fTrgClusters = 0;
	319	return kFALSE;
	320	}
c3ba2d3d	321	}
c3ba2d3d	322	}
6fd5039f	323
6fd5039f	324	return kTRUE;
c3ba2d3d	325	}
	326
	327	//________________________________________________________________________
6fd5039f	328	Bool_t AliAnalysisTaskSAQA::FillHistograms()
c3ba2d3d	329	{
c3604199	330	// Fill histograms.
c3604199	331
c3ba2d3d	332	fHistCentrality->Fill(fCent);
	333	if (fTracks)
	334	fHistTracksCent->Fill(fCent, fTracks->GetEntriesFast());
	335	if (fCaloClusters)
	336	fHistClusCent->Fill(fCent, fCaloClusters->GetEntriesFast());
	337
c3ba2d3d	338	Float_t trackSum = DoTrackLoop();
c3ba2d3d	339
c3ba2d3d	340	DoJetLoop();
c3ba2d3d	341
6fd5039f	342	if (fAnaType == kEMCAL) {
6fd5039f	343	Float_t clusSum = DoClusterLoop();
c3ba2d3d	344
6fd5039f	345	Float_t cellSum = 0, cellCutSum = 0;
	346
	347	DoCellLoop(cellSum, cellCutSum);
	348
	349	fHistChVSneCells->Fill(cellSum, trackSum);
	350	fHistChVSneClus->Fill(clusSum, trackSum);
	351	fHistChVSneCorrCells->Fill(cellCutSum, trackSum);
c3ba2d3d	352
1f6fff78	353	if (fDoTrigger) {
	354	Float_t maxTrgClus = DoTriggerClusLoop();
	355	fHistMaxL1ClusCent->Fill(fCent, maxTrgClus);
6fd5039f	356
1f6fff78	357	Int_t maxL1amp = -1;
1f6fff78	358	Int_t maxL1thr = -1;
6fd5039f	359
1f6fff78	360	DoTriggerPrimitives(maxL1amp, maxL1thr);
	361
	362	if (maxL1amp > -1)
	363	fHistMaxL1FastORCent->Fill(fCent, maxL1amp);
	364
	365	if (maxL1thr > -1)
	366	fHistMaxL1ThrCent->Fill(fCent, maxL1thr);
	367	}
6fd5039f	368	}
c3ba2d3d	369
6fd5039f	370	return kTRUE;
c3ba2d3d	371	}
	372
	373	//________________________________________________________________________
	374	void AliAnalysisTaskSAQA::DoCellLoop(Float_t &sum, Float_t &sum_cut)
	375	{
c3604199	376	// Do cell loop.
c3604199	377
c3ba2d3d	378	AliVCaloCells *cells = InputEvent()->GetEMCALCells();
	379
	380	if (!cells)
	381	return;
	382
c3604199	383	const Int_t ncells = cells->GetNumberOfCells();
c3ba2d3d	384
c3ba2d3d	385	for (Int_t pos = 0; pos < ncells; pos++) {
c3ba2d3d	386	Float_t amp = cells->GetAmplitude(pos);
c3ba2d3d	387	fHistCellsEnergy->Fill(amp);
c3ba2d3d	388	sum += amp;
c3ba2d3d	389	if (amp < fCellEnergyCut)
c3ba2d3d	390	continue;
c3ba2d3d	391	sum_cut += amp;
c3ba2d3d	392	}
	393	}
	394
	395	//________________________________________________________________________
	396	Float_t AliAnalysisTaskSAQA::DoClusterLoop()
	397	{
c3604199	398	// Do cluster loop.
c3604199	399
c3ba2d3d	400	if (!fCaloClusters)
	401	return 0;
	402
	403	Float_t sum = 0;
	404
	405	// Cluster loop
	406	Int_t nclusters = fCaloClusters->GetEntriesFast();
	407
	408	for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
e44e8726	409	AliVCluster* cluster = static_cast<AliVCluster*>(fCaloClusters->At(iClusters));
c3ba2d3d	410	if (!cluster) {
	411	AliError(Form("Could not receive cluster %d", iClusters));
	412	continue;
	413	}
c3ba2d3d	414
e44e8726	415	if (!AcceptCluster(cluster, kTRUE))
e44e8726	416	continue;
c3ba2d3d	417
c3ba2d3d	418	sum += cluster->E();
c3ba2d3d	419
7cf4626b	420	TLorentzVector nPart;
7cf4626b	421	cluster->GetMomentum(nPart, fVertex);
c3ba2d3d	422
7cf4626b	423	fHistClusPhiEtaEnergy->Fill(cluster->E(), nPart.Eta(), nPart.Phi());
7cf4626b	424	fHistNCellsEnergy->Fill(cluster->E(), cluster->GetNCells());
c3604199	425	}
c3ba2d3d	426
	427	return sum;
	428	}
	429
	430	//________________________________________________________________________
	431	Float_t AliAnalysisTaskSAQA::DoTrackLoop()
	432	{
c3604199	433	// Do track loop.
c3604199	434
c3ba2d3d	435	if (!fTracks)
	436	return 0;
	437
	438	Float_t sum = 0;
	439
c3ba2d3d	440	Int_t ntracks = fTracks->GetEntriesFast();
	441	Int_t nclusters = 0;
	442	if (fCaloClusters)
	443	nclusters = fCaloClusters->GetEntriesFast();
	444
c3604199	445	for (Int_t i = 0; i < ntracks; i++) {
c3ba2d3d	446
e44e8726	447	AliVParticle* track = static_cast<AliVParticle*>(fTracks->At(i)); // pointer to reconstructed to track
1f6fff78	448
c3604199	449	if (!track) {
c3ba2d3d	450	AliError(Form("Could not retrieve track %d",i));
	451	continue;
	452	}
c3ba2d3d	453
1f6fff78	454	AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track);
	455
	456	if (vtrack && !AcceptTrack(vtrack, kTRUE))
	457	continue;
	458
c3ba2d3d	459	fHistTracksPt->Fill(track->Pt());
	460
	461	sum += track->P();
	462
	463	Int_t label = track->GetLabel();
	464
	465	fHistTrPhiEta->Fill(track->Eta(), track->Phi());
	466
	467	fHistTrackEta[4]->Fill(track->Eta());
	468	fHistTrackPhi[4]->Fill(track->Phi());
	469
	470	if (label >= 0 && label < 4) {
	471	fHistTrackEta[label]->Fill(track->Eta());
	472	fHistTrackPhi[label]->Fill(track->Phi());
	473	}
	474
1f6fff78	475	if (!vtrack)
c3ba2d3d	476	continue;
c3ba2d3d	477
1f6fff78	478	fHistTrEmcPhiEta->Fill(vtrack->GetTrackEtaOnEMCal(), vtrack->GetTrackPhiOnEMCal());
c3ba2d3d	479	}
	480
	481	return sum;
	482	}
	483
	484	//________________________________________________________________________
	485	void AliAnalysisTaskSAQA::DoJetLoop()
	486	{
c3604199	487	// Do jet loop.
c3604199	488
c3ba2d3d	489	if (!fJets)
	490	return;
	491
	492	Int_t njets = fJets->GetEntriesFast();
	493
	494	for (Int_t ij = 0; ij < njets; ij++) {
	495
e44e8726	496	AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(ij));
c3ba2d3d	497
	498	if (!jet) {
	499	AliError(Form("Could not receive jet %d", ij));
	500	continue;
	501	}
	502
6fd5039f	503	if (!AcceptJet(jet, kFALSE))
c3ba2d3d	504	continue;
6fd5039f	505
99cfd012	506	fHistJetsPtNonBias[fCentBin]->Fill(jet->Pt());
99cfd012	507	fHistJetsPtAreaNonBias[fCentBin]->Fill(jet->Pt(), jet->Area());
6fd5039f	508
	509	if (jet->MaxTrackPt() > fPtBiasJetTrack)
	510	fHistJetsPtTrack[fCentBin]->Fill(jet->Pt());
	511
	512	if (fAnaType == kEMCAL && jet->MaxClusterPt() > fPtBiasJetClus)
	513	fHistJetsPtClus[fCentBin]->Fill(jet->Pt());
c3ba2d3d	514
99cfd012	515	if (!AcceptBiasJet(jet))
	516	continue;
	517
	518	fHistJetsPt[fCentBin]->Fill(jet->Pt());
	519	fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
	520
3000c095	521	fHistJetsPhiEtaPt[fCentBin]->Fill(jet->Eta(), jet->Phi(), jet->Pt());
c3ba2d3d	522	}
	523	}
	524
	525	//________________________________________________________________________
	526	Float_t AliAnalysisTaskSAQA::DoTriggerClusLoop()
	527	{
c3604199	528	// Do trigger cluster loop.
c3604199	529
c3ba2d3d	530	if (!fTrgClusters)
	531	return 0;
	532
	533	Int_t ntrgclusters = fTrgClusters->GetEntriesFast();
	534	Float_t maxTrgClus = 0;
	535
	536	for (Int_t iClusters = 0; iClusters < ntrgclusters; iClusters++) {
e44e8726	537	AliVCluster* cluster = static_cast<AliVCluster*>(fTrgClusters->At(iClusters));
c3ba2d3d	538	if (!cluster) {
	539	AliError(Form("Could not receive cluster %d", iClusters));
	540	continue;
	541	}
	542
	543	if (!(cluster->IsEMCAL())) continue;
	544
	545	if (cluster->E() > maxTrgClus)
	546	maxTrgClus = cluster->E();
	547
	548	}
	549	return maxTrgClus;
	550	}
	551
	552	//________________________________________________________________________
	553	void AliAnalysisTaskSAQA::DoTriggerPrimitives(Int_t &maxL1amp, Int_t &maxL1thr)
	554	{
c3604199	555	// Do trigger primitives loop.
c3604199	556
c3ba2d3d	557	AliVCaloTrigger *triggers = InputEvent()->GetCaloTrigger("EMCAL");
	558
	559	if (!triggers \|\| triggers->GetEntries() == 0)
	560	return;
	561
	562	triggers->Reset();
	563	Int_t L1amp = 0;
	564	Int_t L1thr = 0;
	565	maxL1amp = -1;
	566	maxL1thr = -1;
	567
	568	while (triggers->Next()) {
	569	triggers->GetL1TimeSum(L1amp);
	570	if (maxL1amp < L1amp)
	571	maxL1amp = L1amp;
	572
	573	triggers->GetL1Threshold(L1thr);
	574	if (maxL1thr < L1thr)
	575	maxL1thr = L1thr;
	576	}
	577	}
	578
	579	//________________________________________________________________________
	580	void AliAnalysisTaskSAQA::Terminate(Option_t *)
	581	{
	582	// Called once at the end of the analysis.
	583	}