[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskEmcalDiJetResponse.cxx

//
// Dijet response 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 <TProfile.h>
#include <TChain.h>
#include <TSystem.h>
#include <TFile.h>
#include <TKey.h>

#include "AliVCluster.h"
#include "AliVTrack.h"
#include "AliEmcalJet.h"
#include "AliRhoParameter.h"
#include "AliLog.h"
#include "AliEmcalParticle.h"
#include "AliMCEvent.h"
#include "AliGenPythiaEventHeader.h"
#include "AliAODMCHeader.h"
#include "AliMCEvent.h"
#include "AliAnalysisManager.h"
#include "AliJetContainer.h"

#include "AliAnalysisTaskEmcalDiJetResponse.h"

ClassImp(AliAnalysisTaskEmcalDiJetResponse)

//________________________________________________________________________
AliAnalysisTaskEmcalDiJetResponse::AliAnalysisTaskEmcalDiJetResponse() : 
  AliAnalysisTaskEmcalDiJetBase("AliAnalysisTaskEmcalDiJetResponse"),
  fDoMatchFullCharged(kTRUE),
  fhnDiJetResponseCharged(0),
  fhnDiJetResponseFullCharged(0),
  fh1TriggersLostCharged(0),
  fh1TriggersLostFull(0),
  fh3AssocLostPtDeltaPhiCharged(0),
  fh3AssocLostPtDeltaPhiFull(0),
  fhnMatchingCharged(0),
  fhnMatchingFull(0),
  fnUsedResponseVar(0)
{
  // Default constructor.
  
  SetMakeGeneralHistograms(kTRUE);

  for(Int_t i = 0; i<2; i++) {
    fh1TriggersCharged[i] = 0;
    fh1TriggersFull[i] = 0;
  }

}

//________________________________________________________________________
AliAnalysisTaskEmcalDiJetResponse::AliAnalysisTaskEmcalDiJetResponse(const char *name) : 
  AliAnalysisTaskEmcalDiJetBase(name),
  fDoMatchFullCharged(kTRUE),
  fhnDiJetResponseCharged(0),
  fhnDiJetResponseFullCharged(0),
  fh1TriggersLostCharged(0),
  fh1TriggersLostFull(0),
  fh3AssocLostPtDeltaPhiCharged(0),
  fh3AssocLostPtDeltaPhiFull(0),
  fhnMatchingCharged(0),
  fhnMatchingFull(0),
  fnUsedResponseVar(0)
{
  // Standard constructor.

  SetMakeGeneralHistograms(kTRUE);

  for(Int_t i = 0; i<2; i++) {
    fh1TriggersCharged[i] = 0;
    fh1TriggersFull[i] = 0;
  }

}

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


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

  AliAnalysisTaskEmcalDiJetBase::UserCreateOutputObjects();

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

  //Store dijet vars: pt,trig MC, pt,trig DET, pt,ass MC, pt,ass DET, dPhi MC, dPhi Det, kT MC, kT Det 
  const Int_t nBinsSparse0 = 10;
  const Int_t nBinsPt = 250;
  const Int_t nBinsDPhi     = 36;
  const Int_t nBinsKt       = 25;
  const Int_t nBinsDiJetEta = 40;
  const Int_t nBinsAj       = 50;
  const Int_t nBinsVar[2] = {nBinsKt,nBinsDiJetEta};

  const Int_t nBins0[nBinsSparse0] = {nBinsPt,nBinsPt,nBinsPt,nBinsPt,nBinsDPhi,nBinsDPhi,nBinsVar[fnUsedResponseVar],nBinsVar[fnUsedResponseVar],nBinsAj,nBinsAj};

  const Double_t minPt = 0.;
  const Double_t maxPt = 250.;
  const Double_t minVar[2] = {   0.,-1.};
  const Double_t maxVar[2] = { 100., 1.};

  const Double_t xmin0[nBinsSparse0]  = {  minPt, minPt, minPt, minPt, 0.5*TMath::Pi(), 0.5*TMath::Pi(), minVar[fnUsedResponseVar], minVar[fnUsedResponseVar],0.,0.};
  const Double_t xmax0[nBinsSparse0]  = {  maxPt, maxPt, maxPt, maxPt, 1.5*TMath::Pi(), 1.5*TMath::Pi(), maxVar[fnUsedResponseVar], maxVar[fnUsedResponseVar],1.,1.};

  fhnDiJetResponseCharged = new THnSparseF("fhnDiJetResponseCharged","fhnDiJetResponseCharged;p_{T,trig}^{part};p_{T,trig}^{det};p_{T,ass}^{part};p_{T,ass}^{det};#Delta#varphi_{part};#Delta#varphi_{det};k_{T}^{part},k_{T}^{det};A_{j}^{part}A_{j}^{det}",nBinsSparse0,nBins0,xmin0,xmax0);

  fhnDiJetResponseFullCharged = new THnSparseF("fhnDiJetResponseFullCharged","fhnDiJetResponseFullCharged;p_{T,trig}^{part};p_{T,trig}^{det};p_{T,ass}^{part};p_{T,ass}^{det};#Delta#varphi_{part};#Delta#varphi_{det};k_{T}^{part},k_{T}^{det};A_{j}^{part}A_{j}^{det}",nBinsSparse0,nBins0,xmin0,xmax0);

  if(fnUsedResponseVar==1) {
    fhnDiJetResponseCharged->SetTitle("fhnDiJetResponseCharged DiJetEta"); 
    fhnDiJetResponseCharged->GetAxis(6)->SetTitle("#eta_{dijet}^{part}");
    fhnDiJetResponseCharged->GetAxis(7)->SetTitle("#eta_{dijet}^{det}");

    fhnDiJetResponseFullCharged->SetTitle("fhnDiJetResponseFullCharged DiJetEta"); 
    fhnDiJetResponseFullCharged->GetAxis(6)->SetTitle("#eta_{dijet}^{part}");
    fhnDiJetResponseFullCharged->GetAxis(7)->SetTitle("#eta_{dijet}^{det}");
  }

  fOutput->Add(fhnDiJetResponseCharged);
  fOutput->Add(fhnDiJetResponseFullCharged);

  TString strType = "";
  for(Int_t i = 0; i<2; i++) {
    if(i==0)      strType="Part";
    else if(i==1) strType="Det";
    fh1TriggersCharged[i] = new TH1F(Form("fh1TriggersCharged%s",strType.Data()),Form("fh1TriggersCharged%s;p_{T,trig}^{ch}",strType.Data()),nBinsPt,minPt,maxPt);
    fOutput->Add(fh1TriggersCharged[i]);
  
    fh1TriggersFull[i] = new TH1F(Form("fh1TriggersFull%s",strType.Data()),Form("fh1TriggersFull%s;p_{T,trig}^{ch}",strType.Data()),nBinsPt,minPt,maxPt);
    fOutput->Add(fh1TriggersFull[i]);
  }

  fh1TriggersLostCharged = new TH1F("fh1TriggersLostCharged","fh1TriggersLostCharged;p_{T,trig}^{ch}",nBinsPt,minPt,maxPt);
  fOutput->Add(fh1TriggersLostCharged);

  fh1TriggersLostFull = new TH1F("fh1TriggersLostFull","fh1TriggersLostFull;p_{T,trig}^{ch}",nBinsPt,minPt,maxPt);
  fOutput->Add(fh1TriggersLostFull);

  fh3AssocLostPtDeltaPhiCharged = new TH3F("fh3AssocLostPtDeltaPhiCharged","fh3AssocLostPtDeltaPhiCharged;p_{T,trig}^{ch};p_{T,assoc}^{ch};#Delta#varphi",nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,nBinsDPhi,-0.5*TMath::Pi(),1.5*TMath::Pi());
  fOutput->Add(fh3AssocLostPtDeltaPhiCharged);

  fh3AssocLostPtDeltaPhiFull = new TH3F("fh3AssocLostPtDeltaPhiFull","fh3AssocLostPtDeltaPhiFull;p_{T,trig}^{ch};p_{T,assoc}^{ch};#Delta#varphi",nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,nBinsDPhi,-0.5*TMath::Pi(),1.5*TMath::Pi());
  fOutput->Add(fh3AssocLostPtDeltaPhiFull);

  const Int_t nBinsSparseMatch = 6;
  const Int_t nBinsDPhiMatch = 80;
  const Int_t nBinsDEtaMatch = 80;
  const Int_t nBinsDR = 20;
  const Int_t nBinsType = 3;
  const Int_t nBinsMatch[nBinsSparseMatch] = {nBinsPt,nBinsPt,nBinsDPhiMatch,nBinsDEtaMatch,nBinsDR,nBinsType};
  //pTpart, pTdet, deltaPhi, deltaEta, deltaR, jet type (leading,subleading,other)
  const Double_t xminMatch[nBinsSparseMatch]  = { minPt, minPt, -0.5,-0.5, 0., 0};
  const Double_t xmaxMatch[nBinsSparseMatch]  = { maxPt, maxPt,  0.5, 0.5, 0.5,3};
  fhnMatchingCharged = new THnSparseF("fhnMatchingCharged","fhnMatchingCharged;#it{p}_{T,part} (GeV/#it{c});#it{p}_{T,det} (GeV/#it{c});#Delta#varphi;#Delta#eta;#Delta R;type",
					  nBinsSparseMatch,nBinsMatch,xminMatch,xmaxMatch);
  fOutput->Add(fhnMatchingCharged);

  fhnMatchingFull = new THnSparseF("fhnMatchingFull","fhnMatchingFull;#it{p}_{T,part} (GeV/#it{c});#it{p}_{T,det} (GeV/#it{c});#Delta#varphi;#Delta#eta;#Delta R;type",
					  nBinsSparseMatch,nBinsMatch,xminMatch,xmaxMatch);
  fOutput->Add(fhnMatchingFull);


  // =========== 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;
    }
    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 AliAnalysisTaskEmcalDiJetResponse::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;

  if(fRhoType==0) {
    fRhoFullVal = 0.;
    fRhoChVal = 0.;
  }
  if(fRhoType==1) {
    fRhoFullVal = GetRhoVal(fContainerFull);
    fRhoChVal = GetRhoVal(fContainerCharged);
  }
  

  //Do matching
  MatchJetsGeo(fContainerCharged,fContainerChargedMC,0,0.3,1);
  MatchJetsGeo(fContainerFull,fContainerFullMC,0,0.3,2);

  //Fill particle-detector level matching histos
  
  if(fDoChargedCharged)   CorrelateJets(1);

  if(fDoFullCharged) {
    SetChargedFractionIndexMC();
    CorrelateJets(2);
  }
  
  return kTRUE;  // If return kFALSE FillHistogram() will NOT be executed.
}

//________________________________________________________________________
void AliAnalysisTaskEmcalDiJetResponse::CorrelateJets(const Int_t type) {
  //
  // Correlate jets and fill histos
  //

  if( fJetCorrelationType==kCorrelateAll )
    CorrelateAllJets(type);
  else if( fJetCorrelationType==kCorrelateTwo )
    CorrelateTwoJets(type);
  else if( fJetCorrelationType==kCorrelateLS )
    AliWarning(Form("%s: leading-subleading correlation not implemented for response!",GetName()));

  return;

}


//________________________________________________________________________
void AliAnalysisTaskEmcalDiJetResponse::CorrelateAllJets(const Int_t type) {
  //
  // Correlate jets and fill histos
  //

  Int_t typet = 0;
  Int_t typetMC = 0;
  Int_t typeaMC = 0;
  if(type==0) { //full-full
    typetMC = fContainerFullMC;
    typeaMC = fContainerFullMC;
    typet = fContainerFull;
  }
  else if(type==1) { //charged-charged
    typetMC = fContainerChargedMC;
    typeaMC = fContainerChargedMC;
    typet = fContainerCharged;
  }
  else if(type==2) { //full-charged
    typetMC = fContainerFullMC;
    typeaMC = fContainerChargedMC;
    typet = fContainerFull;
  }
  else {
    AliWarning(Form("%s: type %d of dijet correlation not defined!",GetName(),type));
    return;
  }

  Int_t nJetsTrig  = 0;
  Int_t nJetsAssoc = 0;
  if(type==0) {
    nJetsTrig  = GetNJets(fContainerFullMC);
    nJetsAssoc = nJetsTrig;
  }
  else if(type==1) {
    nJetsTrig  = GetNJets(fContainerChargedMC);
    nJetsAssoc = nJetsTrig;
  }
  else if(type==2) {
    nJetsTrig  = GetNJets(fContainerFullMC);
    nJetsAssoc = GetNJets(fContainerChargedMC);
  }


  for(Int_t ijt=0; ijt<nJetsTrig; ijt++) {

    AliEmcalJet *jetTrigMC = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ijt, typetMC));
    if(!jetTrigMC) continue; //jet not selected

    Double_t jetTrigPtMC = GetJetPt(jetTrigMC,typetMC);

    if(jetTrigPtMC<fPtMinTriggerJet)
      continue;

    if(type==1)
      fh1TriggersCharged[0]->Fill(jetTrigPtMC);
    if(type==2)
      fh1TriggersFull[0]->Fill(jetTrigPtMC);

    AliEmcalJet *jetTrigDet = jetTrigMC->ClosestJet();
    if(!jetTrigDet) {
      //trigger is lost
      if(type==1)
	fh1TriggersLostCharged->Fill(jetTrigPtMC);
      if(type==2)
	fh1TriggersLostFull->Fill(jetTrigPtMC);
      
      continue;
    }

    if(type==1)
      fh1TriggersCharged[1]->Fill(GetJetPt(jetTrigDet,typet));
    if(type==2)
      fh1TriggersFull[1]->Fill(GetJetPt(jetTrigDet,typet));

    for(Int_t ija=0; ija<nJetsAssoc; ija++) {
      if(IsSameJet(ijt,ija,type,kTRUE)) continue;

      AliEmcalJet *jetAssocMC = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ija, typeaMC));
      if(!jetAssocMC) continue;

      Double_t jetAssocPtMC = GetJetPt(jetAssocMC,typeaMC);

      //Now check if jets are also there on detector level
      AliEmcalJet *jetAssocDet = jetAssocMC->ClosestJet();
      if(!jetAssocDet) {
	//dijet is lost
      if(type==1)
	fh3AssocLostPtDeltaPhiCharged->Fill(jetTrigPtMC,jetAssocPtMC,GetDeltaPhi(jetTrigMC,jetAssocMC));
      if(type==2)
	fh3AssocLostPtDeltaPhiFull->Fill(jetTrigPtMC,jetAssocPtMC,GetDeltaPhi(jetTrigMC,jetAssocMC));
	continue;
      }

      FillDiJetResponse(jetTrigMC,jetAssocMC,jetTrigDet,jetAssocDet,type);

    } // associate jet loop
  }//trigger jet loop

}

//________________________________________________________________________
void AliAnalysisTaskEmcalDiJetResponse::CorrelateTwoJets(const Int_t type) {
  //
  // Correlate jets and fill histos
  //

  Int_t typet = 0;
  Int_t typea = 0;
  Int_t typetMC = 0;
  Int_t typeaMC = 0;
  if(type==0) { //full-full
    typetMC = fContainerFullMC;
    typeaMC = fContainerFullMC;
    typet = fContainerFull;
    typea = fContainerFull;
  }
  else if(type==1) { //charged-charged
    typetMC = fContainerChargedMC;
    typeaMC = fContainerChargedMC;
    typet = fContainerCharged;
    typea = fContainerCharged;
  }
  else if(type==2) { //full-charged
    typetMC = fContainerFullMC;
    typeaMC = fContainerChargedMC;
    typet = fContainerFull;
    typea = fContainerCharged;
  }
  else {
    AliWarning(Form("%s: type %d of dijet correlation not defined!",GetName(),type));
    return;
  }

  Int_t nJetsTrig  = 0;
  if(type==0) {
    nJetsTrig  = GetNJets(fContainerFullMC);
  }
  else if(type==1) {
    nJetsTrig  = GetNJets(fContainerChargedMC);
  }
  else if(type==2) {
    nJetsTrig  = GetNJets(fContainerFullMC);
  }

  for(Int_t ijt=0; ijt<nJetsTrig; ijt++) {

    AliEmcalJet *jetTrigMC = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ijt, typetMC));
    if(!jetTrigMC) continue; //jet not selected

    Double_t jetTrigPtMC = GetJetPt(jetTrigMC,typetMC);

    if(jetTrigPtMC<fPtMinTriggerJet)
      continue;

    if(type==1)
      fh1TriggersCharged[0]->Fill(jetTrigPtMC);
    if(type==2)
      fh1TriggersFull[0]->Fill(jetTrigPtMC);

    AliEmcalJet *jetTrigDet = jetTrigMC->ClosestJet();
    if(!jetTrigDet) {
      //trigger is lost
      if(type==1)
	fh1TriggersLostCharged->Fill(jetTrigPtMC);
      if(type==2)
	fh1TriggersLostFull->Fill(jetTrigPtMC);
      continue;
    }

    if(type==1)
      fh1TriggersCharged[1]->Fill(GetJetPt(jetTrigDet,typet));
    if(type==2)
      fh1TriggersFull[1]->Fill(GetJetPt(jetTrigDet,typet));


    AliEmcalJet *jetAssocMC = GetLeadingJetOppositeHemisphere(type,typeaMC,jetTrigMC);
    if(!jetAssocMC) continue;

    Double_t jetAssocPtMC = GetJetPt(jetAssocMC,typeaMC);
      
    //Now check if jets are also there on detector level
    AliEmcalJet *jetAssocDet = jetAssocMC->ClosestJet();
    if(!jetAssocDet) {
      //dijet is lost
      if(type==1)
	fh3AssocLostPtDeltaPhiCharged->Fill(jetTrigPtMC,jetAssocPtMC,GetDeltaPhi(jetTrigMC,jetAssocMC));
      if(type==2)
	fh3AssocLostPtDeltaPhiFull->Fill(jetTrigPtMC,jetAssocPtMC,GetDeltaPhi(jetTrigMC,jetAssocMC));
      continue;
    }

    if(fDoPtBias) {
      if(type==0 || type==1) {
	if(GetJetPt(jetAssocDet,typea)>GetJetPt(jetTrigDet,typet))
	  continue;
      }
    }

    FillDiJetResponse(jetTrigMC,jetAssocMC,jetTrigDet,jetAssocDet,type);


  }//trigger jet loop

}

//________________________________________________________________________
void AliAnalysisTaskEmcalDiJetResponse::FillDiJetResponse(const AliEmcalJet *jetTrigMC, const AliEmcalJet *jetAssocMC, const AliEmcalJet *jetTrigDet, const AliEmcalJet *jetAssocDet, Int_t type) {

  //Fill dijet response

  Int_t typet = 0;
  Int_t typea = 0;
  Int_t typetMC = 0;
  Int_t typeaMC = 0;
  if(type==0) { //full-full
    typetMC = fContainerFullMC;
    typeaMC = fContainerFullMC;
    typet = fContainerFull;
    typea = fContainerFull;
  }
  else if(type==1) { //charged-charged
    typetMC = fContainerChargedMC;
    typeaMC = fContainerChargedMC;
    typet = fContainerCharged;
    typea = fContainerCharged;
  }
  else if(type==2) { //full-charged
    typetMC = fContainerFullMC;
    typeaMC = fContainerChargedMC;
    typet = fContainerFull;
    typea = fContainerCharged;
  }
  else {
    AliWarning(Form("%s: type %d of dijet correlation not defined!",GetName(),type));
    return;
  }

  Double_t jetTrigPtMC  = GetJetPt(jetTrigMC,typetMC);
  Double_t jetAssocPtMC = GetJetPt(jetAssocMC,typeaMC);

  Double_t varDet[2] = {TMath::Abs(GetJetPt(jetTrigDet,typet)*TMath::Sin(GetDeltaPhi(jetTrigDet,jetAssocDet))),(jetTrigDet->Eta()+jetAssocDet->Eta())/2.};
  Double_t varPart[2] = {TMath::Abs(jetTrigPtMC*TMath::Sin(GetDeltaPhi(jetTrigMC,jetAssocMC))),(jetTrigMC->Eta()+jetAssocMC->Eta())/2.};

  Double_t ajDet  = (GetJetPt(jetTrigDet,typet)-GetJetPt(jetAssocDet,typea))/(GetJetPt(jetTrigDet,typet)+GetJetPt(jetAssocDet,typea));
  Double_t ajPart = (jetTrigPtMC-jetAssocPtMC)/(jetTrigPtMC+jetAssocPtMC);

  //Store dijet vars: pt,trig MC; pt,trig DET; pt,ass MC; pt,ass DET; dPhi MC; dPhi Det; kT MC; kT Det;
  Double_t diJetVars[10] = {
    jetTrigPtMC,
    GetJetPt(jetTrigDet,typet),
    jetAssocPtMC,
    GetJetPt(jetAssocDet,typea),
    GetDeltaPhi(jetTrigMC,jetAssocMC),
    GetDeltaPhi(jetTrigDet,jetAssocDet),
    varPart[fnUsedResponseVar],
    varDet[fnUsedResponseVar],
    ajDet,
    ajPart
  }; 
  
  if(type==1)
    fhnDiJetResponseCharged->Fill(diJetVars);
  else if(type==2)
    fhnDiJetResponseFullCharged->Fill(diJetVars);


}

//________________________________________________________________________
void AliAnalysisTaskEmcalDiJetResponse::FillMatchHistos() {
  //
  // Fill Particle-Detector level matching histos
  //

  for(int i = 0; i < GetNJets(fContainerFull);++i) {
    AliEmcalJet *jetDet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(i, fContainerFull));
    if(!jetDet) continue;

    AliEmcalJet *jetPart = jetDet->ClosestJet();
    if(!jetPart) continue;

    Double_t matchVars[6] = {
      jetPart->Pt(),
      jetDet->Pt(),
      GetDeltaPhi(jetPart->Phi(),jetDet->Phi()),
      jetPart->Eta()-jetDet->Eta(),
      GetDeltaR(jetPart,jetDet),
      TMath::Min((Float_t)i+0.5,2.5)
    };
    fhnMatchingFull->Fill(matchVars);

  }//loop over full jets

  for(int i = 0; i < GetNJets(fContainerCharged);++i) {
    AliEmcalJet *jetDet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(i, fContainerCharged));
    if(!jetDet) continue;

    AliEmcalJet *jetPart = jetDet->ClosestJet();
    if(!jetPart) continue;

    Double_t matchVars[6] = {
      jetPart->Pt(),
      jetDet->Pt(),
      GetDeltaPhi(jetPart->Phi(),jetDet->Phi()),
      jetPart->Eta()-jetDet->Eta(),
      GetDeltaR(jetPart,jetDet),
      TMath::Min((Float_t)i+0.5,2.5)
    };
    fhnMatchingCharged->Fill(matchVars);

  }//loop over charged jets

}


//________________________________________________________________________
Bool_t AliAnalysisTaskEmcalDiJetResponse::RetrieveEventObjects() {
  //
  // retrieve event objects
  //

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

  return kTRUE;

}

//_______________________________________________________________________
void AliAnalysisTaskEmcalDiJetResponse::Terminate(Option_t *) 
{
  // Called once at the end of the analysis.
}
Commit	Line	Data
6e8b6371	1	//
	2	// Dijet response analysis task.
	3	//
	4	// Author: M.Verweij
	5
	6	#include <TClonesArray.h>
	7	#include <TH1F.h>
	8	#include <TH2F.h>
	9	#include <TH3F.h>
	10	#include <THnSparse.h>
	11	#include <TList.h>
	12	#include <TLorentzVector.h>
	13	#include <TProfile.h>
	14	#include <TChain.h>
	15	#include <TSystem.h>
	16	#include <TFile.h>
	17	#include <TKey.h>
	18
	19	#include "AliVCluster.h"
	20	#include "AliVTrack.h"
	21	#include "AliEmcalJet.h"
	22	#include "AliRhoParameter.h"
	23	#include "AliLog.h"
6e8b6371	24	#include "AliEmcalParticle.h"
	25	#include "AliMCEvent.h"
	26	#include "AliGenPythiaEventHeader.h"
	27	#include "AliAODMCHeader.h"
	28	#include "AliMCEvent.h"
	29	#include "AliAnalysisManager.h"
	30	#include "AliJetContainer.h"
	31
	32	#include "AliAnalysisTaskEmcalDiJetResponse.h"
	33
	34	ClassImp(AliAnalysisTaskEmcalDiJetResponse)
	35
	36	//________________________________________________________________________
	37	AliAnalysisTaskEmcalDiJetResponse::AliAnalysisTaskEmcalDiJetResponse() :
	38	AliAnalysisTaskEmcalDiJetBase("AliAnalysisTaskEmcalDiJetResponse"),
	39	fDoMatchFullCharged(kTRUE),
	40	fhnDiJetResponseCharged(0),
	41	fhnDiJetResponseFullCharged(0),
	42	fh1TriggersLostCharged(0),
	43	fh1TriggersLostFull(0),
67c36e97	44	fh3AssocLostPtDeltaPhiCharged(0),
67c36e97	45	fh3AssocLostPtDeltaPhiFull(0),
6e8b6371	46	fhnMatchingCharged(0),
6ab30d5f	47	fhnMatchingFull(0),
6ab30d5f	48	fnUsedResponseVar(0)
6e8b6371	49	{
6e8b6371	50	// Default constructor.
6e8b6371	51
	52	SetMakeGeneralHistograms(kTRUE);
	53
	54	for(Int_t i = 0; i<2; i++) {
	55	fh1TriggersCharged[i] = 0;
	56	fh1TriggersFull[i] = 0;
	57	}
	58
	59	}
	60
	61	//________________________________________________________________________
	62	AliAnalysisTaskEmcalDiJetResponse::AliAnalysisTaskEmcalDiJetResponse(const char *name) :
	63	AliAnalysisTaskEmcalDiJetBase(name),
	64	fDoMatchFullCharged(kTRUE),
	65	fhnDiJetResponseCharged(0),
	66	fhnDiJetResponseFullCharged(0),
	67	fh1TriggersLostCharged(0),
	68	fh1TriggersLostFull(0),
67c36e97	69	fh3AssocLostPtDeltaPhiCharged(0),
67c36e97	70	fh3AssocLostPtDeltaPhiFull(0),
6e8b6371	71	fhnMatchingCharged(0),
6ab30d5f	72	fhnMatchingFull(0),
6ab30d5f	73	fnUsedResponseVar(0)
6e8b6371	74	{
	75	// Standard constructor.
	76
	77	SetMakeGeneralHistograms(kTRUE);
	78
	79	for(Int_t i = 0; i<2; i++) {
	80	fh1TriggersCharged[i] = 0;
	81	fh1TriggersFull[i] = 0;
	82	}
	83
	84	}
	85
	86	//________________________________________________________________________
	87	AliAnalysisTaskEmcalDiJetResponse::~AliAnalysisTaskEmcalDiJetResponse()
	88	{
	89	// Destructor.
	90	}
	91
	92
	93	//________________________________________________________________________
	94	void AliAnalysisTaskEmcalDiJetResponse::UserCreateOutputObjects()
	95	{
	96	// Create user output.
	97
6e8b6371	98	AliAnalysisTaskEmcalDiJetBase::UserCreateOutputObjects();
	99
	100	Bool_t oldStatus = TH1::AddDirectoryStatus();
	101	TH1::AddDirectory(kFALSE);
	102
	103	//Store dijet vars: pt,trig MC, pt,trig DET, pt,ass MC, pt,ass DET, dPhi MC, dPhi Det, kT MC, kT Det
e50460de	104	const Int_t nBinsSparse0 = 10;
6e8b6371	105	const Int_t nBinsPt = 250;
e50460de	106	const Int_t nBinsDPhi = 36;
e50460de	107	const Int_t nBinsKt = 25;
6ab30d5f	108	const Int_t nBinsDiJetEta = 40;
e50460de	109	const Int_t nBinsAj = 50;
6ab30d5f	110	const Int_t nBinsVar[2] = {nBinsKt,nBinsDiJetEta};
6ab30d5f	111
e50460de	112	const Int_t nBins0[nBinsSparse0] = {nBinsPt,nBinsPt,nBinsPt,nBinsPt,nBinsDPhi,nBinsDPhi,nBinsVar[fnUsedResponseVar],nBinsVar[fnUsedResponseVar],nBinsAj,nBinsAj};
6e8b6371	113
	114	const Double_t minPt = 0.;
	115	const Double_t maxPt = 250.;
e50460de	116	const Double_t minVar[2] = { 0.,-1.};
6ab30d5f	117	const Double_t maxVar[2] = { 100., 1.};
6e8b6371	118
e50460de	119	const Double_t xmin0[nBinsSparse0] = { minPt, minPt, minPt, minPt, 0.5TMath::Pi(), 0.5TMath::Pi(), minVar[fnUsedResponseVar], minVar[fnUsedResponseVar],0.,0.};
e50460de	120	const Double_t xmax0[nBinsSparse0] = { maxPt, maxPt, maxPt, maxPt, 1.5TMath::Pi(), 1.5TMath::Pi(), maxVar[fnUsedResponseVar], maxVar[fnUsedResponseVar],1.,1.};
6e8b6371	121
e50460de	122	fhnDiJetResponseCharged = new THnSparseF("fhnDiJetResponseCharged","fhnDiJetResponseCharged;p_{T,trig}^{part};p_{T,trig}^{det};p_{T,ass}^{part};p_{T,ass}^{det};#Delta#varphi_{part};#Delta#varphi_{det};k_{T}^{part},k_{T}^{det};A_{j}^{part}A_{j}^{det}",nBinsSparse0,nBins0,xmin0,xmax0);
6e8b6371	123
e50460de	124	fhnDiJetResponseFullCharged = new THnSparseF("fhnDiJetResponseFullCharged","fhnDiJetResponseFullCharged;p_{T,trig}^{part};p_{T,trig}^{det};p_{T,ass}^{part};p_{T,ass}^{det};#Delta#varphi_{part};#Delta#varphi_{det};k_{T}^{part},k_{T}^{det};A_{j}^{part}A_{j}^{det}",nBinsSparse0,nBins0,xmin0,xmax0);
6ab30d5f	125
	126	if(fnUsedResponseVar==1) {
	127	fhnDiJetResponseCharged->SetTitle("fhnDiJetResponseCharged DiJetEta");
	128	fhnDiJetResponseCharged->GetAxis(6)->SetTitle("#eta_{dijet}^{part}");
	129	fhnDiJetResponseCharged->GetAxis(7)->SetTitle("#eta_{dijet}^{det}");
	130
	131	fhnDiJetResponseFullCharged->SetTitle("fhnDiJetResponseFullCharged DiJetEta");
	132	fhnDiJetResponseFullCharged->GetAxis(6)->SetTitle("#eta_{dijet}^{part}");
	133	fhnDiJetResponseFullCharged->GetAxis(7)->SetTitle("#eta_{dijet}^{det}");
	134	}
	135
	136	fOutput->Add(fhnDiJetResponseCharged);
6e8b6371	137	fOutput->Add(fhnDiJetResponseFullCharged);
	138
	139	TString strType = "";
	140	for(Int_t i = 0; i<2; i++) {
	141	if(i==0) strType="Part";
	142	else if(i==1) strType="Det";
	143	fh1TriggersCharged[i] = new TH1F(Form("fh1TriggersCharged%s",strType.Data()),Form("fh1TriggersCharged%s;p_{T,trig}^{ch}",strType.Data()),nBinsPt,minPt,maxPt);
	144	fOutput->Add(fh1TriggersCharged[i]);
	145
	146	fh1TriggersFull[i] = new TH1F(Form("fh1TriggersFull%s",strType.Data()),Form("fh1TriggersFull%s;p_{T,trig}^{ch}",strType.Data()),nBinsPt,minPt,maxPt);
	147	fOutput->Add(fh1TriggersFull[i]);
	148	}
	149
	150	fh1TriggersLostCharged = new TH1F("fh1TriggersLostCharged","fh1TriggersLostCharged;p_{T,trig}^{ch}",nBinsPt,minPt,maxPt);
	151	fOutput->Add(fh1TriggersLostCharged);
	152
	153	fh1TriggersLostFull = new TH1F("fh1TriggersLostFull","fh1TriggersLostFull;p_{T,trig}^{ch}",nBinsPt,minPt,maxPt);
	154	fOutput->Add(fh1TriggersLostFull);
	155
67c36e97	156	fh3AssocLostPtDeltaPhiCharged = new TH3F("fh3AssocLostPtDeltaPhiCharged","fh3AssocLostPtDeltaPhiCharged;p_{T,trig}^{ch};p_{T,assoc}^{ch};#Delta#varphi",nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,nBinsDPhi,-0.5TMath::Pi(),1.5TMath::Pi());
	157	fOutput->Add(fh3AssocLostPtDeltaPhiCharged);
	158
	159	fh3AssocLostPtDeltaPhiFull = new TH3F("fh3AssocLostPtDeltaPhiFull","fh3AssocLostPtDeltaPhiFull;p_{T,trig}^{ch};p_{T,assoc}^{ch};#Delta#varphi",nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,nBinsDPhi,-0.5TMath::Pi(),1.5TMath::Pi());
	160	fOutput->Add(fh3AssocLostPtDeltaPhiFull);
	161
6e8b6371	162	const Int_t nBinsSparseMatch = 6;
	163	const Int_t nBinsDPhiMatch = 80;
	164	const Int_t nBinsDEtaMatch = 80;
	165	const Int_t nBinsDR = 20;
	166	const Int_t nBinsType = 3;
	167	const Int_t nBinsMatch[nBinsSparseMatch] = {nBinsPt,nBinsPt,nBinsDPhiMatch,nBinsDEtaMatch,nBinsDR,nBinsType};
	168	//pTpart, pTdet, deltaPhi, deltaEta, deltaR, jet type (leading,subleading,other)
	169	const Double_t xminMatch[nBinsSparseMatch] = { minPt, minPt, -0.5,-0.5, 0., 0};
	170	const Double_t xmaxMatch[nBinsSparseMatch] = { maxPt, maxPt, 0.5, 0.5, 0.5,3};
	171	fhnMatchingCharged = new THnSparseF("fhnMatchingCharged","fhnMatchingCharged;#it{p}_{T,part} (GeV/#it{c});#it{p}_{T,det} (GeV/#it{c});#Delta#varphi;#Delta#eta;#Delta R;type",
	172	nBinsSparseMatch,nBinsMatch,xminMatch,xmaxMatch);
	173	fOutput->Add(fhnMatchingCharged);
	174
	175	fhnMatchingFull = new THnSparseF("fhnMatchingFull","fhnMatchingFull;#it{p}_{T,part} (GeV/#it{c});#it{p}_{T,det} (GeV/#it{c});#Delta#varphi;#Delta#eta;#Delta R;type",
	176	nBinsSparseMatch,nBinsMatch,xminMatch,xmaxMatch);
	177	fOutput->Add(fhnMatchingFull);
	178
	179
	180	// =========== Switch on Sumw2 for all histos ===========
	181	for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
	182	TH1 h1 = dynamic_cast<TH1>(fOutput->At(i));
	183	if (h1){
	184	h1->Sumw2();
	185	continue;
	186	}
	187	THnSparse hn = dynamic_cast<THnSparse>(fOutput->At(i));
	188	if(hn)hn->Sumw2();
	189	}
	190
	191	TH1::AddDirectory(oldStatus);
	192
	193	PostData(1, fOutput); // Post data for ALL output slots > 0 here.
	194	}
	195
6e8b6371	196	//________________________________________________________________________
	197	Bool_t AliAnalysisTaskEmcalDiJetResponse::Run()
	198	{
	199	// Run analysis code here, if needed. It will be executed before FillHistograms().
	200
	201	//Check if event is selected (vertex & pile-up)
	202	if(!SelectEvent())
	203	return kFALSE;
	204
6e8b6371	205	if(fRhoType==0) {
	206	fRhoFullVal = 0.;
	207	fRhoChVal = 0.;
	208	}
	209	if(fRhoType==1) {
	210	fRhoFullVal = GetRhoVal(fContainerFull);
	211	fRhoChVal = GetRhoVal(fContainerCharged);
	212	}
	213
	214
	215	//Do matching
	216	MatchJetsGeo(fContainerCharged,fContainerChargedMC,0,0.3,1);
	217	MatchJetsGeo(fContainerFull,fContainerFullMC,0,0.3,2);
	218
	219	//Fill particle-detector level matching histos
	220
	221	if(fDoChargedCharged) CorrelateJets(1);
	222
	223	if(fDoFullCharged) {
	224	SetChargedFractionIndexMC();
	225	CorrelateJets(2);
	226	}
	227
	228	return kTRUE; // If return kFALSE FillHistogram() will NOT be executed.
	229	}
	230
	231	//________________________________________________________________________
	232	void AliAnalysisTaskEmcalDiJetResponse::CorrelateJets(const Int_t type) {
	233	//
	234	// Correlate jets and fill histos
	235	//
	236
6ab30d5f	237	if( fJetCorrelationType==kCorrelateAll )
	238	CorrelateAllJets(type);
	239	else if( fJetCorrelationType==kCorrelateTwo )
	240	CorrelateTwoJets(type);
	241	else if( fJetCorrelationType==kCorrelateLS )
	242	AliWarning(Form("%s: leading-subleading correlation not implemented for response!",GetName()));
	243
	244	return;
	245
	246	}
	247
	248
	249	//________________________________________________________________________
	250	void AliAnalysisTaskEmcalDiJetResponse::CorrelateAllJets(const Int_t type) {
	251	//
	252	// Correlate jets and fill histos
	253	//
	254
6e8b6371	255	Int_t typet = 0;
6e8b6371	256	Int_t typetMC = 0;
	257	Int_t typeaMC = 0;
	258	if(type==0) { //full-full
	259	typetMC = fContainerFullMC;
	260	typeaMC = fContainerFullMC;
	261	typet = fContainerFull;
6e8b6371	262	}
	263	else if(type==1) { //charged-charged
	264	typetMC = fContainerChargedMC;
	265	typeaMC = fContainerChargedMC;
	266	typet = fContainerCharged;
6e8b6371	267	}
	268	else if(type==2) { //full-charged
	269	typetMC = fContainerFullMC;
	270	typeaMC = fContainerChargedMC;
	271	typet = fContainerFull;
6e8b6371	272	}
	273	else {
	274	AliWarning(Form("%s: type %d of dijet correlation not defined!",GetName(),type));
	275	return;
	276	}
	277
	278	Int_t nJetsTrig = 0;
	279	Int_t nJetsAssoc = 0;
	280	if(type==0) {
	281	nJetsTrig = GetNJets(fContainerFullMC);
	282	nJetsAssoc = nJetsTrig;
	283	}
	284	else if(type==1) {
	285	nJetsTrig = GetNJets(fContainerChargedMC);
	286	nJetsAssoc = nJetsTrig;
	287	}
	288	else if(type==2) {
	289	nJetsTrig = GetNJets(fContainerFullMC);
	290	nJetsAssoc = GetNJets(fContainerChargedMC);
	291	}
	292
	293
	294	for(Int_t ijt=0; ijt<nJetsTrig; ijt++) {
	295
	296	AliEmcalJet jetTrigMC = static_cast<AliEmcalJet>(GetAcceptJetFromArray(ijt, typetMC));
	297	if(!jetTrigMC) continue; //jet not selected
	298
	299	Double_t jetTrigPtMC = GetJetPt(jetTrigMC,typetMC);
	300
	301	if(jetTrigPtMC<fPtMinTriggerJet)
	302	continue;
	303
	304	if(type==1)
	305	fh1TriggersCharged[0]->Fill(jetTrigPtMC);
	306	if(type==2)
	307	fh1TriggersFull[0]->Fill(jetTrigPtMC);
	308
	309	AliEmcalJet *jetTrigDet = jetTrigMC->ClosestJet();
	310	if(!jetTrigDet) {
	311	//trigger is lost
	312	if(type==1)
	313	fh1TriggersLostCharged->Fill(jetTrigPtMC);
	314	if(type==2)
	315	fh1TriggersLostFull->Fill(jetTrigPtMC);
	316
	317	continue;
	318	}
	319
	320	if(type==1)
	321	fh1TriggersCharged[1]->Fill(GetJetPt(jetTrigDet,typet));
	322	if(type==2)
	323	fh1TriggersFull[1]->Fill(GetJetPt(jetTrigDet,typet));
	324
	325	for(Int_t ija=0; ija<nJetsAssoc; ija++) {
	326	if(IsSameJet(ijt,ija,type,kTRUE)) continue;
	327
	328	AliEmcalJet jetAssocMC = static_cast<AliEmcalJet>(GetAcceptJetFromArray(ija, typeaMC));
	329	if(!jetAssocMC) continue;
	330
	331	Double_t jetAssocPtMC = GetJetPt(jetAssocMC,typeaMC);
	332
	333	//Now check if jets are also there on detector level
	334	AliEmcalJet *jetAssocDet = jetAssocMC->ClosestJet();
67c36e97	335	if(!jetAssocDet) {
6e8b6371	336	//dijet is lost
67c36e97	337	if(type==1)
	338	fh3AssocLostPtDeltaPhiCharged->Fill(jetTrigPtMC,jetAssocPtMC,GetDeltaPhi(jetTrigMC,jetAssocMC));
	339	if(type==2)
	340	fh3AssocLostPtDeltaPhiFull->Fill(jetTrigPtMC,jetAssocPtMC,GetDeltaPhi(jetTrigMC,jetAssocMC));
6e8b6371	341	continue;
	342	}
	343
6ab30d5f	344	FillDiJetResponse(jetTrigMC,jetAssocMC,jetTrigDet,jetAssocDet,type);
	345
	346	} // associate jet loop
	347	}//trigger jet loop
	348
	349	}
	350
	351	//________________________________________________________________________
	352	void AliAnalysisTaskEmcalDiJetResponse::CorrelateTwoJets(const Int_t type) {
	353	//
	354	// Correlate jets and fill histos
	355	//
	356
	357	Int_t typet = 0;
	358	Int_t typea = 0;
	359	Int_t typetMC = 0;
	360	Int_t typeaMC = 0;
	361	if(type==0) { //full-full
	362	typetMC = fContainerFullMC;
	363	typeaMC = fContainerFullMC;
	364	typet = fContainerFull;
	365	typea = fContainerFull;
	366	}
	367	else if(type==1) { //charged-charged
	368	typetMC = fContainerChargedMC;
	369	typeaMC = fContainerChargedMC;
	370	typet = fContainerCharged;
	371	typea = fContainerCharged;
	372	}
	373	else if(type==2) { //full-charged
	374	typetMC = fContainerFullMC;
	375	typeaMC = fContainerChargedMC;
	376	typet = fContainerFull;
	377	typea = fContainerCharged;
	378	}
	379	else {
	380	AliWarning(Form("%s: type %d of dijet correlation not defined!",GetName(),type));
	381	return;
	382	}
	383
	384	Int_t nJetsTrig = 0;
6ab30d5f	385	if(type==0) {
6ab30d5f	386	nJetsTrig = GetNJets(fContainerFullMC);
6ab30d5f	387	}
	388	else if(type==1) {
	389	nJetsTrig = GetNJets(fContainerChargedMC);
6ab30d5f	390	}
	391	else if(type==2) {
	392	nJetsTrig = GetNJets(fContainerFullMC);
6ab30d5f	393	}
6e8b6371	394
6ab30d5f	395	for(Int_t ijt=0; ijt<nJetsTrig; ijt++) {
	396
	397	AliEmcalJet jetTrigMC = static_cast<AliEmcalJet>(GetAcceptJetFromArray(ijt, typetMC));
	398	if(!jetTrigMC) continue; //jet not selected
	399
	400	Double_t jetTrigPtMC = GetJetPt(jetTrigMC,typetMC);
	401
	402	if(jetTrigPtMC<fPtMinTriggerJet)
	403	continue;
	404
	405	if(type==1)
	406	fh1TriggersCharged[0]->Fill(jetTrigPtMC);
	407	if(type==2)
	408	fh1TriggersFull[0]->Fill(jetTrigPtMC);
	409
	410	AliEmcalJet *jetTrigDet = jetTrigMC->ClosestJet();
	411	if(!jetTrigDet) {
	412	//trigger is lost
6e8b6371	413	if(type==1)
6ab30d5f	414	fh1TriggersLostCharged->Fill(jetTrigPtMC);
	415	if(type==2)
	416	fh1TriggersLostFull->Fill(jetTrigPtMC);
6ab30d5f	417	continue;
	418	}
	419
	420	if(type==1)
	421	fh1TriggersCharged[1]->Fill(GetJetPt(jetTrigDet,typet));
	422	if(type==2)
	423	fh1TriggersFull[1]->Fill(GetJetPt(jetTrigDet,typet));
	424
	425
	426	AliEmcalJet *jetAssocMC = GetLeadingJetOppositeHemisphere(type,typeaMC,jetTrigMC);
	427	if(!jetAssocMC) continue;
6e8b6371	428
6ab30d5f	429	Double_t jetAssocPtMC = GetJetPt(jetAssocMC,typeaMC);
6e8b6371	430
6ab30d5f	431	//Now check if jets are also there on detector level
	432	AliEmcalJet *jetAssocDet = jetAssocMC->ClosestJet();
	433	if(!jetAssocDet) {
	434	//dijet is lost
	435	if(type==1)
	436	fh3AssocLostPtDeltaPhiCharged->Fill(jetTrigPtMC,jetAssocPtMC,GetDeltaPhi(jetTrigMC,jetAssocMC));
	437	if(type==2)
	438	fh3AssocLostPtDeltaPhiFull->Fill(jetTrigPtMC,jetAssocPtMC,GetDeltaPhi(jetTrigMC,jetAssocMC));
	439	continue;
	440	}
	441
f9de9b73	442	if(fDoPtBias) {
	443	if(type==0 \|\| type==1) {
	444	if(GetJetPt(jetAssocDet,typea)>GetJetPt(jetTrigDet,typet))
	445	continue;
	446	}
3c709670	447	}
3c709670	448
6ab30d5f	449	FillDiJetResponse(jetTrigMC,jetAssocMC,jetTrigDet,jetAssocDet,type);
	450
	451
6e8b6371	452	}//trigger jet loop
	453
	454	}
	455
6ab30d5f	456	//________________________________________________________________________
	457	void AliAnalysisTaskEmcalDiJetResponse::FillDiJetResponse(const AliEmcalJet jetTrigMC, const AliEmcalJet jetAssocMC, const AliEmcalJet jetTrigDet, const AliEmcalJet jetAssocDet, Int_t type) {
	458
	459	//Fill dijet response
	460
	461	Int_t typet = 0;
	462	Int_t typea = 0;
	463	Int_t typetMC = 0;
	464	Int_t typeaMC = 0;
	465	if(type==0) { //full-full
	466	typetMC = fContainerFullMC;
	467	typeaMC = fContainerFullMC;
	468	typet = fContainerFull;
	469	typea = fContainerFull;
	470	}
	471	else if(type==1) { //charged-charged
	472	typetMC = fContainerChargedMC;
	473	typeaMC = fContainerChargedMC;
	474	typet = fContainerCharged;
	475	typea = fContainerCharged;
	476	}
	477	else if(type==2) { //full-charged
	478	typetMC = fContainerFullMC;
	479	typeaMC = fContainerChargedMC;
	480	typet = fContainerFull;
	481	typea = fContainerCharged;
	482	}
	483	else {
	484	AliWarning(Form("%s: type %d of dijet correlation not defined!",GetName(),type));
	485	return;
	486	}
	487
	488	Double_t jetTrigPtMC = GetJetPt(jetTrigMC,typetMC);
	489	Double_t jetAssocPtMC = GetJetPt(jetAssocMC,typeaMC);
	490
e50460de	491	Double_t varDet[2] = {TMath::Abs(GetJetPt(jetTrigDet,typet)*TMath::Sin(GetDeltaPhi(jetTrigDet,jetAssocDet))),(jetTrigDet->Eta()+jetAssocDet->Eta())/2.};
	492	Double_t varPart[2] = {TMath::Abs(jetTrigPtMC*TMath::Sin(GetDeltaPhi(jetTrigMC,jetAssocMC))),(jetTrigMC->Eta()+jetAssocMC->Eta())/2.};
	493
	494	Double_t ajDet = (GetJetPt(jetTrigDet,typet)-GetJetPt(jetAssocDet,typea))/(GetJetPt(jetTrigDet,typet)+GetJetPt(jetAssocDet,typea));
	495	Double_t ajPart = (jetTrigPtMC-jetAssocPtMC)/(jetTrigPtMC+jetAssocPtMC);
6ab30d5f	496
6ab30d5f	497	//Store dijet vars: pt,trig MC; pt,trig DET; pt,ass MC; pt,ass DET; dPhi MC; dPhi Det; kT MC; kT Det;
e50460de	498	Double_t diJetVars[10] = {
6ab30d5f	499	jetTrigPtMC,
	500	GetJetPt(jetTrigDet,typet),
	501	jetAssocPtMC,
	502	GetJetPt(jetAssocDet,typea),
	503	GetDeltaPhi(jetTrigMC,jetAssocMC),
	504	GetDeltaPhi(jetTrigDet,jetAssocDet),
	505	varPart[fnUsedResponseVar],
e50460de	506	varDet[fnUsedResponseVar],
	507	ajDet,
	508	ajPart
6ab30d5f	509	};
	510
	511	if(type==1)
	512	fhnDiJetResponseCharged->Fill(diJetVars);
	513	else if(type==2)
	514	fhnDiJetResponseFullCharged->Fill(diJetVars);
	515
	516
	517	}
	518
6e8b6371	519	//________________________________________________________________________
	520	void AliAnalysisTaskEmcalDiJetResponse::FillMatchHistos() {
	521	//
	522	// Fill Particle-Detector level matching histos
	523	//
	524
	525	for(int i = 0; i < GetNJets(fContainerFull);++i) {
	526	AliEmcalJet jetDet = static_cast<AliEmcalJet>(GetAcceptJetFromArray(i, fContainerFull));
	527	if(!jetDet) continue;
	528
	529	AliEmcalJet *jetPart = jetDet->ClosestJet();
	530	if(!jetPart) continue;
	531
67c36e97	532	Double_t matchVars[6] = {
6e8b6371	533	jetPart->Pt(),
	534	jetDet->Pt(),
	535	GetDeltaPhi(jetPart->Phi(),jetDet->Phi()),
	536	jetPart->Eta()-jetDet->Eta(),
	537	GetDeltaR(jetPart,jetDet),
	538	TMath::Min((Float_t)i+0.5,2.5)
	539	};
	540	fhnMatchingFull->Fill(matchVars);
	541
	542	}//loop over full jets
	543
	544	for(int i = 0; i < GetNJets(fContainerCharged);++i) {
	545	AliEmcalJet jetDet = static_cast<AliEmcalJet>(GetAcceptJetFromArray(i, fContainerCharged));
	546	if(!jetDet) continue;
	547
	548	AliEmcalJet *jetPart = jetDet->ClosestJet();
	549	if(!jetPart) continue;
	550
67c36e97	551	Double_t matchVars[6] = {
6e8b6371	552	jetPart->Pt(),
	553	jetDet->Pt(),
	554	GetDeltaPhi(jetPart->Phi(),jetDet->Phi()),
	555	jetPart->Eta()-jetDet->Eta(),
	556	GetDeltaR(jetPart,jetDet),
	557	TMath::Min((Float_t)i+0.5,2.5)
	558	};
	559	fhnMatchingCharged->Fill(matchVars);
	560
	561	}//loop over charged jets
	562
	563	}
	564
	565
	566	//________________________________________________________________________
	567	Bool_t AliAnalysisTaskEmcalDiJetResponse::RetrieveEventObjects() {
	568	//
	569	// retrieve event objects
	570	//
	571
	572	if (!AliAnalysisTaskEmcalDiJetBase::RetrieveEventObjects())
	573	return kFALSE;
	574
	575	return kTRUE;
	576
	577	}
	578
	579	//_______________________________________________________________________
	580	void AliAnalysisTaskEmcalDiJetResponse::Terminate(Option_t *)
	581	{
	582	// Called once at the end of the analysis.
	583	}