//
// Jet trigger QA analysis task.
//
// Author: M.Verweij

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

#include "AliVCluster.h"
#include "AliVTrack.h"
#include "AliEmcalJet.h"
#include "AliRhoParameter.h"
#include "AliLog.h"
#include "AliAnalysisUtils.h"
#include "AliEmcalParticle.h"
#include "AliAODCaloTrigger.h"
#include "AliEMCALGeometry.h"
#include "AliVCaloCells.h"


#include "AliAnalysisTaskEmcalJetTriggerQA.h"

ClassImp(AliAnalysisTaskEmcalJetTriggerQA)

//________________________________________________________________________
AliAnalysisTaskEmcalJetTriggerQA::AliAnalysisTaskEmcalJetTriggerQA() : 
  AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetTriggerQA", kTRUE),
  fDebug(kFALSE),
  fUseAnaUtils(kTRUE),
  fAnalysisUtils(0),
  fJetsNameCh("Jet_AKTChargedR040_PicoTracks_pT0150_caloClustersCorr_ET0300"),
  fJetsCh(0),
  fEtaMaxChJet(0.5),
  fPhiMinChJet(0.),
  fPhiMaxChJet(TMath::TwoPi()),
  fMaxTrackPtChJet(100.),
  fRhoChName(""),
  fRhoCh(0),
  fRhoChVal(0),
  fTriggerClass(""),
  fBitJ1((1<<8)),
  fBitJ2((1<<11)),
  fMaxPatchEnergy(0),
  fTriggerType(-1),
  fNFastOR(16),
  fhNEvents(0),
  fh3PtEtaPhiJetFull(0),
  fh3PtEtaPhiJetCharged(0),
  fh2NJetsPtFull(0),
  fh2NJetsPtCharged(0),
  fh3PtEtaAreaJetFull(0),
  fh3PtEtaAreaJetCharged(0),
  fh2PtNConstituentsCharged(0),
  fh2PtNConstituents(0),
  fh2PtMeanPtConstituentsCharged(0),
  fh2PtMeanPtConstituentsNeutral(0),
  fh2PtNEF(0),
  fh2Ptz(0),
  fh2PtLeadJet1VsLeadJet2(0),
  fh3EEtaPhiCluster(0),
  fh3PtLeadJet1VsPatchEnergy(0),
  fh3PtLeadJet2VsPatchEnergy(0),
  fh3PatchEnergyEtaPhiCenter(0)
{
  // Default constructor.

  
  SetMakeGeneralHistograms(kTRUE);
}

//________________________________________________________________________
AliAnalysisTaskEmcalJetTriggerQA::AliAnalysisTaskEmcalJetTriggerQA(const char *name) : 
  AliAnalysisTaskEmcalJet(name, kTRUE),
  fDebug(kFALSE),
  fUseAnaUtils(kTRUE),
  fAnalysisUtils(0),
  fJetsNameCh("Jet_AKTChargedR040_PicoTracks_pT0150_caloClustersCorr_ET0300"),
  fJetsCh(0),
  fEtaMaxChJet(0.5),
  fPhiMinChJet(0.),
  fPhiMaxChJet(TMath::TwoPi()),
  fMaxTrackPtChJet(100.),
  fRhoChName(""),
  fRhoCh(0),
  fRhoChVal(0),
  fTriggerClass(""),
  fBitJ1((1<<8)),
  fBitJ2((1<<11)),
  fMaxPatchEnergy(0),
  fTriggerType(-1),
  fNFastOR(16),
  fhNEvents(0),
  fh3PtEtaPhiJetFull(0),
  fh3PtEtaPhiJetCharged(0),
  fh2NJetsPtFull(0),
  fh2NJetsPtCharged(0),
  fh3PtEtaAreaJetFull(0),
  fh3PtEtaAreaJetCharged(0),
  fh2PtNConstituentsCharged(0),
  fh2PtNConstituents(0),
  fh2PtMeanPtConstituentsCharged(0),
  fh2PtMeanPtConstituentsNeutral(0),
  fh2PtNEF(0),
  fh2Ptz(0),
  fh2PtLeadJet1VsLeadJet2(0),
  fh3EEtaPhiCluster(0),
  fh3PtLeadJet1VsPatchEnergy(0),
  fh3PtLeadJet2VsPatchEnergy(0),
  fh3PatchEnergyEtaPhiCenter(0)
{
  // Standard constructor.

  SetMakeGeneralHistograms(kTRUE);
}

//________________________________________________________________________
AliAnalysisTaskEmcalJetTriggerQA::~AliAnalysisTaskEmcalJetTriggerQA()
{
  // Destructor.
}

//----------------------------------------------------------------------------------------------
void AliAnalysisTaskEmcalJetTriggerQA::InitOnce() {
  //
  // only initialize once
  //

  // Initialize analysis util class for vertex selection
  if(fUseAnaUtils) {
    fAnalysisUtils = new AliAnalysisUtils();
    fAnalysisUtils->SetMinVtxContr(2);
    fAnalysisUtils->SetMaxVtxZ(10.);
  }
}

//________________________________________________________________________
Bool_t AliAnalysisTaskEmcalJetTriggerQA::SelectEvent() {
  //
  // Decide if event should be selected for analysis
  //

  fhNEvents->Fill(0.5);
  
  if(fAnalysisUtils) {
    if(!fAnalysisUtils->IsVertexSelected2013pA(InputEvent()))
      return kFALSE;

    fhNEvents->Fill(2.5);

    if(fAnalysisUtils->IsPileUpEvent(InputEvent()))
      return kFALSE;
  }
  else{
    if(fUseAnaUtils)
      AliError("fAnalysisUtils not initialized. Call AliAnalysisTaskEmcalJetTriggerQA::InitOnce()");
  }

  fhNEvents->Fill(3.5);

  if(!fTriggerClass.IsNull()) {
    //Check if requested trigger was fired
    TString firedTrigClass = InputEvent()->GetFiredTriggerClasses();
    if(!firedTrigClass.Contains(fTriggerClass))
      return kFALSE;
  }

  fhNEvents->Fill(1.5);

  return kTRUE;

}

//________________________________________________________________________
void AliAnalysisTaskEmcalJetTriggerQA::FindTriggerPatch() {

  //Code to get position of trigger

  if(!fGeom)
    fGeom = AliEMCALGeometry::GetInstance();


  TString firedTrigClass = InputEvent()->GetFiredTriggerClasses();

  AliAODCaloTrigger *trg = dynamic_cast<AliAODCaloTrigger*>(InputEvent()->GetCaloTrigger("EMCAL"));
  trg->Reset();

  int col, row; //FASTOR position

  Int_t nPatchNotEmpty = 0; //counter number of patches which are not empty
  fMaxPatchEnergy = 0.;
  while (trg->Next())
    {

      trg->GetPosition(col, row); //col (0 to 63), row (0 to 47)
      
      if (col > -1 && row > -1)
	{
	  Int_t id = -1; //FASTOR index
	  Int_t cellIndex[4] = {-1};
	  fGeom->GetAbsFastORIndexFromPositionInEMCAL(col,row,id); //phi is column, eta is row
	  if(id<0) continue;

	  fGeom->GetCellIndexFromFastORIndex(id,cellIndex);

	  Int_t ts = 0;
	  trg->GetL1TimeSum(ts);

	  Float_t ampTrg = 0.;
	  trg->GetAmplitude(ampTrg);

	  //L1 analysis
	  Int_t bit = 0;
	  trg->GetTriggerBits(bit);

	  Bool_t bTrigJ = TestFilterBit(bit,fBitJ1|fBitJ2);
	  if(!bTrigJ)
	    continue;

	  Bool_t bTrigJ1 = TestFilterBit(bit,fBitJ1);	   
	  Bool_t bTrigJ2 = TestFilterBit(bit,fBitJ2);	   

	  if(bTrigJ1) fTriggerType = 0;
	  if(bTrigJ2) fTriggerType = 1;
	  if(!bTrigJ1 && !bTrigJ2) fTriggerType = -1;

	  Double_t minPhiPatch =  10.;
	  Double_t maxPhiPatch = -10.;
	  Double_t minEtaPatch =  10.;
	  Double_t maxEtaPatch = -10.;

	  //Get energy in trigger patch 8x8 FASTOR
	  Double_t patchEnergy = 0.;
	  Double_t sumAmp = 0.;
	  //	  const Int_t nFastOR = 8;//16;
	  for(Int_t fastrow = 0; fastrow<fNFastOR; fastrow++) {
	    for(Int_t fastcol = 0; fastcol<fNFastOR; fastcol++) {
	      Int_t nrow = row+fastrow;
	      Int_t ncol = col+fastcol;
	      fGeom->GetAbsFastORIndexFromPositionInEMCAL(ncol,nrow,id);

	      if(id<0) {
		AliWarning(Form("%s: id smaller than 0 %d",GetName(),id));
		continue;
	      }
	      
	      fGeom->GetCellIndexFromFastORIndex(id,cellIndex);
	      for(Int_t icell=0; icell<4; icell++) {

		if(!fGeom->CheckAbsCellId(cellIndex[icell])) continue;

		Double_t amp =0., time = 0., efrac = 0;
		Int_t mclabel = -1;
		Short_t absId = -1;
		Int_t nSupMod = -1, nModule = -1, nIphi = -1, nIeta = -1;
 
		fGeom->GetCellIndex(cellIndex[icell], nSupMod, nModule, nIphi, nIeta );
		
		fCaloCells->GetCell(cellIndex[icell], absId, amp, time,mclabel,efrac);

		Double_t eta, phi;
		fGeom->EtaPhiFromIndex(cellIndex[icell], eta, phi);

		if(phi<minPhiPatch) minPhiPatch = phi;
		if(phi>maxPhiPatch) maxPhiPatch = phi;
		if(eta<minEtaPatch) minEtaPatch = eta;
		if(eta>maxEtaPatch) maxEtaPatch = eta;

		sumAmp+=fCaloCells->GetAmplitude(cellIndex[icell]);
		patchEnergy+=amp;

	      }//cells in fastor loop
	    }//fastor col loop
	  }//fastor row loop

	  Double_t etaCent = minEtaPatch + 0.5*(maxEtaPatch-minEtaPatch);
	  Double_t phiCent = minPhiPatch + 0.5*(maxPhiPatch-minPhiPatch);
	  fh3PatchEnergyEtaPhiCenter->Fill(patchEnergy,etaCent,phiCent);

	  if(patchEnergy>0)
	    if(fDebug>2) AliInfo(Form("%s: patch edges eta: %f - %f   phi: %f - %f ampTrg: %f patchEnergy: %f sumAmp: %f",GetName(),minEtaPatch,maxEtaPatch,minPhiPatch,maxPhiPatch,ampTrg,patchEnergy,sumAmp));

	  if(patchEnergy>0) nPatchNotEmpty++;
	  
	  if(patchEnergy>fMaxPatchEnergy) fMaxPatchEnergy=patchEnergy;
 
	}
    }

}

//________________________________________________________________________
void AliAnalysisTaskEmcalJetTriggerQA::LoadExtraBranches() {
  //
  // get charged jets
  //

  if (!fJetsNameCh.IsNull() && !fJetsCh) {
    fJetsCh = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsNameCh.Data()));
    if (!fJetsCh) {
      AliError(Form("%s: Could not retrieve charged jets %s!", GetName(), fJetsNameCh.Data()));
      return;
    }
    else if (!fJetsCh->GetClass()->GetBaseClass("AliEmcalJet")) {
      AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fJetsNameCh.Data()));
      fJetsCh = 0;
      return;
    }
  }

  //Get Charged Rho
  if (!fRhoChName.IsNull() && !fRhoCh) {
    fRhoCh = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhoChName.Data()));
    if (!fRhoCh) {
      AliError(Form("%s: Could not retrieve charged rho %s!", GetName(), fRhoChName.Data()));
      return;
    }
    else
      fRhoChVal = 0.;
  }
  if (fRhoChName.IsNull() ) {
    if(fDebug>10) AliInfo(Form("%s: fRhoChName empty %s",GetName(),fRhoChName.Data()));
    fRhoChVal = 0.;
  }
  else if(!fRhoChName.IsNull() && fRhoCh) //do not use rho if rhoType==0
    fRhoChVal = fRhoCh->GetVal();

}

//________________________________________________________________________
Bool_t AliAnalysisTaskEmcalJetTriggerQA::AcceptChargedJet(const AliEmcalJet *jet) const {

  // Accept charged jet

  Bool_t accept = kFALSE;

  if(jet->Pt()<0.15) //reject ghost jets
    return accept;

  //do area cut
  if(jet->Area()<fJetAreaCut)
    return accept;

  //Check if jet is within EMCAL acceptance
  Double_t eta = jet->Eta();
  Double_t phi = jet->Phi();
  if(TMath::Abs(eta)>fEtaMaxChJet) 
    return accept;

  if(phi<fPhiMinChJet || phi>fPhiMaxChJet)
    return accept;

  if (jet->MaxTrackPt() > fMaxTrackPt)
    return kFALSE;

  return kTRUE;

}

//________________________________________________________________________
void AliAnalysisTaskEmcalJetTriggerQA::UserCreateOutputObjects()
{
  // Create user output.

  InitOnce();

  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();

  Bool_t oldStatus = TH1::AddDirectoryStatus();
  TH1::AddDirectory(kFALSE);

  fhNEvents = new TH1F("fhNEvents","fhNEvents;selection;N_{evt}",5,0,5);
  fOutput->Add(fhNEvents);

  Int_t nBinsPt = 120;
  Double_t minPt = -20.;
  Double_t maxPt = 100.;
  Int_t nBinsEta = 40;
  Double_t minEta = -1.;
  Double_t maxEta = 1.;
  Int_t nBinsPhi = 18*6;
  Double_t minPhi = 0.;
  Double_t maxPhi = TMath::TwoPi();
  Int_t nBinsArea = 100;
  Double_t minArea = 0.;
  Double_t maxArea = 1.;

  fh3PtEtaPhiJetFull = new TH3F("fh3PtEtaPhiJetFull","fh3PtEtaPhiJetFull;#it{p}_{T}^{jet};#eta;#varphi",nBinsPt,minPt,maxPt,nBinsEta,minEta,maxEta,nBinsPhi,minPhi,maxPhi);
  fOutput->Add(fh3PtEtaPhiJetFull);

  fh3PtEtaPhiJetCharged = new TH3F("fh3PtEtaPhiJetCharged","fh3PtEtaPhiJetCharged;#it{p}_{T}^{jet};#eta;#varphi",nBinsPt,minPt,maxPt,nBinsEta,minEta,maxEta,nBinsPhi,minPhi,maxPhi);
  fOutput->Add(fh3PtEtaPhiJetCharged);

  fh2NJetsPtFull = new TH2F("fh2NJetsPtFull","fh2NJetsPtFull;N_{jets};#it{p}_{T}^{jet}",20,-0.5,19.5,nBinsPt,minPt,maxPt);
  fOutput->Add(fh2NJetsPtFull);

  fh2NJetsPtCharged = new TH2F("fh2NJetsPtCharged","fh2NJetsPtCharged;N_{jets};#it{p}_{T}^{jet}",20,-0.5,19.5,nBinsPt,minPt,maxPt);
  fOutput->Add(fh2NJetsPtCharged);

  fh3PtEtaAreaJetFull = new TH3F("fh3PtEtaAreaJetFull","fh3PtEtaAreaJetFull;#it{p}_{T}^{jet};#eta;A",nBinsPt,minPt,maxPt,nBinsEta,minEta,maxEta,nBinsArea,minArea,maxArea);
  fOutput->Add(fh3PtEtaAreaJetFull);

  fh3PtEtaAreaJetCharged = new TH3F("fh3PtEtaAreaJetCharged","fh3PtEtaAreaJetCharged;#it{p}_{T}^{jet};#eta;A",nBinsPt,minPt,maxPt,nBinsEta,minEta,maxEta,nBinsArea,minArea,maxArea);
  fOutput->Add(fh3PtEtaAreaJetCharged);

  Int_t nBinsConst =100;
  Double_t minConst = 0.;
  Double_t maxConst = 100.;

  Int_t nBinsMeanPt = 200;
  Double_t minMeanPt = 0.;
  Double_t maxMeanPt = 20.;

  Int_t nBinsNEF = 101;
  Double_t minNEF = 0.;
  Double_t maxNEF = 1.01;

  Int_t nBinsz = 101;
  Double_t minz = 0.;
  Double_t maxz = 1.01;

  Int_t nBinsECluster =100;
  Double_t minECluster = 0.;
  Double_t maxECluster = 100.;


  fh2PtNConstituentsCharged = new TH2F("fh2PtNConstituentsCharged","fh2PtNConstituentsCharged;#it{p}_{T}^{jet};N_{charged constituents}",nBinsPt,minPt,maxPt,nBinsConst,minConst,maxConst);
  fOutput->Add(fh2PtNConstituentsCharged);

  fh2PtNConstituents = new TH2F("fh2PtNConstituents","fh2PtNConstituents;#it{p}_{T}^{jet};N_{constituents}",nBinsPt,minPt,maxPt,nBinsConst,minConst,maxConst);
  fOutput->Add(fh2PtNConstituents);

  fh2PtMeanPtConstituentsCharged = new TH2F("fh2PtMeanPtConstituentsCharged","fh2PtMeanPtConstituentsCharged;#it{p}_{T}^{jet};charged #langle #it{p}_{T} #rangle",nBinsPt,minPt,maxPt,nBinsMeanPt,minMeanPt,maxMeanPt);
  fOutput->Add(fh2PtMeanPtConstituentsCharged);

  fh2PtMeanPtConstituentsNeutral = new TH2F("fh2PtMeanPtConstituentsNeutral","fh2PtMeanPtConstituentsNeutral;#it{p}_{T}^{jet};neutral langle #it{p}_{T} #rangle",nBinsPt,minPt,maxPt,nBinsMeanPt,minMeanPt,maxMeanPt);
  fOutput->Add(fh2PtMeanPtConstituentsNeutral);

  fh2PtNEF = new TH2F("fh2PtNEF","fh2PtNEF;#it{p}_{T}^{jet};NEF",nBinsPt,minPt,maxPt,nBinsNEF,minNEF,maxNEF);
  fOutput->Add(fh2PtNEF);

  fh2Ptz = new TH2F("fh2Ptz","fh2Ptz;#it{p}_{T}^{jet};z=p_{t,trk}^{proj}/p_{jet}",nBinsPt,minPt,maxPt,nBinsz,minz,maxz);
  fOutput->Add(fh2Ptz);

  fh2PtLeadJet1VsLeadJet2 = new TH2F("fh2PtLeadJet1VsLeadJet2","fh2PtLeadJet1VsLeadJet2;#it{p}_{T}^{jet 1};#it{p}_{T}^{jet 2}",nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
  fOutput->Add(fh2PtLeadJet1VsLeadJet2);

  fh3EEtaPhiCluster = new TH3F("fh3EEtaPhiCluster","fh3EEtaPhiCluster;E_{clus};#eta,#phi",nBinsECluster,minECluster,maxECluster,nBinsEta,minEta,maxEta,nBinsPhi,minPhi,maxPhi);
  fOutput->Add(fh3EEtaPhiCluster);

  fh3PtLeadJet1VsPatchEnergy = new TH3F("fh3PtLeadJet1VsPatchEnergy","fh3PtLeadJet1VsPatchEnergy;#it{p}_{T}^{jet 1};Amplitude_{patch};trig type",nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,2,-0.5,1.5);
  fOutput->Add(fh3PtLeadJet1VsPatchEnergy);
  fh3PtLeadJet2VsPatchEnergy = new TH3F("fh3PtLeadJet2VsPatchEnergy","fh3PtLeadJet2VsPatchEnergy;#it{p}_{T}^{jet 1};Amplitude_{patch};trig type",nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,2,-0.5,1.5);
  fOutput->Add(fh3PtLeadJet2VsPatchEnergy);

  fh3PatchEnergyEtaPhiCenter = new TH3F("fh3PatchEnergyEtaPhiCenter","fh3PatchEnergyEtaPhiCenter;E_{patch};#eta;#phi",nBinsPt,minPt,maxPt,nBinsEta,minEta,maxEta,nBinsPhi,minPhi,maxPhi);
  fOutput->Add(fh3PatchEnergyEtaPhiCenter);


  // =========== Switch on Sumw2 for all histos ===========
  for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
    TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
    if (h1){
      h1->Sumw2();
      continue;
    }
    TH2 *h2 = dynamic_cast<TH2*>(fOutput->At(i));
    if (h2){
      h2->Sumw2();
      continue;
    }
    TH3 *h3 = dynamic_cast<TH3*>(fOutput->At(i));
    if (h3){
      h3->Sumw2();
      continue;
    }
    THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
    if(hn)hn->Sumw2();
  }

  TH1::AddDirectory(oldStatus);

  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
}

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

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

      TLorentzVector lp;
      cluster->GetMomentum(lp, const_cast<Double_t*>(fVertex));

      //Fill eta,phi,E of clusters here
      fh3EEtaPhiCluster->Fill(lp.E(),lp.Eta(),lp.Phi());
    }
  }

  Double_t ptLeadJet1 = 0.;
  Double_t ptLeadJet2 = 0.;

  TArrayI *nJetsArr = new TArrayI(fh2NJetsPtFull->GetNbinsY()+1);
  nJetsArr->Reset(0);
  nJetsArr->Set(fh2NJetsPtFull->GetNbinsY()+1);

  if (fJets) {
    const 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))
	continue;

      Double_t jetPt = jet->Pt();
      if(jetPt>ptLeadJet1) ptLeadJet1=jetPt;
      fh3PtEtaPhiJetFull->Fill(jetPt,jet->Eta(),jet->Phi());
      fh3PtEtaAreaJetFull->Fill(jetPt,jet->Eta(),jet->Area());
      
      //count jets above certain pT threshold
      Int_t ptbin = fh2NJetsPtFull->GetYaxis()->FindBin(jetPt);
      for(Int_t iptbin = ptbin; iptbin<=fh2NJetsPtFull->GetNbinsY(); iptbin++)
	nJetsArr->AddAt(nJetsArr->At(iptbin)+1,iptbin);
      
      fh2PtNConstituentsCharged->Fill(jetPt,jet->GetNumberOfTracks());
      fh2PtNConstituents->Fill(jetPt,jet->GetNumberOfConstituents());

      fh2PtNEF->Fill(jetPt,jet->NEF());

      AliVParticle *vp;
      Double_t sumPtCh = 0.;
      for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
	vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
	if(!vp) continue;
	fh2Ptz->Fill(jetPt,GetZ(vp,jet));
	
	sumPtCh+=vp->Pt();
	
      }
      
      if(jet->GetNumberOfTracks()>0)
	fh2PtMeanPtConstituentsCharged->Fill(jetPt,sumPtCh/(double)(jet->GetNumberOfTracks()) );


      AliVCluster *vc = 0x0;
      Double_t sumPtNe = 0.;
      for(Int_t icc=0; icc<jet->GetNumberOfClusters(); icc++) {
	vc = static_cast<AliVCluster*>(jet->ClusterAt(icc, fCaloClusters));
	if(!vc) continue;

	TLorentzVector lp;
	vc->GetMomentum(lp, const_cast<Double_t*>(fVertex));
	
	sumPtNe+=lp.Pt();

      }

      if(jet->GetNumberOfClusters()>0)
	fh2PtMeanPtConstituentsNeutral->Fill(jetPt,sumPtNe/(double)(jet->GetNumberOfClusters()) );

    } //full jet loop

    for(Int_t i=1; i<=fh2NJetsPtFull->GetNbinsY(); i++) {
      Int_t nJetsInEvent = nJetsArr->At(i);
      fh2NJetsPtFull->Fill(nJetsInEvent,fh2NJetsPtFull->GetYaxis()->GetBinCenter(i));
    }

  }

  //Reset array to zero to also count charged jets
  nJetsArr->Reset(0);
  
  //Loop over charged jets
  if (fJetsCh) {
    const Int_t njetsCh = fJetsCh->GetEntriesFast();
    for (Int_t ij = 0; ij < njetsCh; ij++) {

      AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJetsCh->At(ij));
      if (!jet) {
	AliError(Form("Could not receive charged jet %d", ij));
	continue;
      }

      if(!AcceptChargedJet(jet)) 
	continue;
      
      Double_t jetPt = jet->Pt();
      if(jetPt>ptLeadJet2) ptLeadJet2=jetPt;
      fh3PtEtaPhiJetCharged->Fill(jetPt,jet->Eta(),jet->Phi());
      fh3PtEtaAreaJetCharged->Fill(jetPt,jet->Eta(),jet->Area());
      
      //count jets above certain pT threshold
      Int_t ptbin = fh2NJetsPtCharged->GetYaxis()->FindBin(jetPt);
      for(Int_t iptbin = ptbin; iptbin<=fh2NJetsPtCharged->GetNbinsY(); iptbin++)
	nJetsArr->AddAt(nJetsArr->At(iptbin)+1,iptbin);
      
    }//ch jet loop
    for(Int_t i=1; i<=fh2NJetsPtCharged->GetNbinsY(); i++) {
      Int_t nJetsInEvent = nJetsArr->At(i);
      fh2NJetsPtCharged->Fill(nJetsInEvent,fh2NJetsPtCharged->GetYaxis()->GetBinCenter(i));
    }
  }

  if(fJets && fJetsCh) {
    fh2PtLeadJet1VsLeadJet2->Fill(ptLeadJet1,ptLeadJet2);
  }

  fh3PtLeadJet1VsPatchEnergy->Fill(ptLeadJet1,fMaxPatchEnergy,fTriggerType);
  fh3PtLeadJet2VsPatchEnergy->Fill(ptLeadJet2,fMaxPatchEnergy,fTriggerType);

  if(nJetsArr) delete nJetsArr;

  return kTRUE;
}

//________________________________________________________________________
Bool_t AliAnalysisTaskEmcalJetTriggerQA::Run()
{
  // Run analysis code here, if needed. It will be executed before FillHistograms().

  //Check if event is selected (vertex & pile-up)
  if(!SelectEvent())
    return kFALSE;
  
  LoadExtraBranches();

  if(fRhoChName.IsNull())
    fRhoChVal = 0.;
  else
    fRhoChVal = fRhoCh->GetVal();

  FindTriggerPatch();

  return kTRUE;  // If return kFALSE FillHistogram() will NOT be executed.
}

//_______________________________________________________________________
void AliAnalysisTaskEmcalJetTriggerQA::Terminate(Option_t *) 
{
  // Called once at the end of the analysis.
}
//________________________________________________________________________
Double_t AliAnalysisTaskEmcalJetTriggerQA::GetZ(const AliVParticle *trk, const AliEmcalJet *jet)          const
{  
  // Get Z of constituent trk

  return GetZ(trk->Px(),trk->Py(),trk->Pz(),jet->Px(),jet->Py(),jet->Pz());
}

//________________________________________________________________________
Double_t AliAnalysisTaskEmcalJetTriggerQA::GetZ(const Double_t trkPx, const Double_t trkPy, const Double_t trkPz, const Double_t jetPx, const Double_t jetPy, const Double_t jetPz) const
{
  // 
  // Get the z of a constituent inside of a jet
  //

  Double_t pJetSq = jetPx*jetPx+jetPy*jetPy+jetPz*jetPz;

  if(pJetSq>0.)
    return (trkPx*jetPx+trkPy*jetPy+trkPz*jetPz)/pJetSq;
  else {
    AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f",GetName(), pJetSq)); 
    return 0;
  }
}