// Di-Jet angular correlations class
// Author: Greeshma Koyithatta Meethaleveedu and Ragahva Varma

#include "TChain.h"
#include "TTree.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TH3F.h"
#include "TCanvas.h"
#include "TList.h"
#include "TObjArray.h"
#include "THnSparse.h"

#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"
#include "AliAnalysisTaskSE.h"
#include "AliCentrality.h"
#include "AliMultiplicity.h"
#include "AliESDEvent.h"
#include "AliESDtrackCuts.h"
#include "AliESDInputHandler.h"
#include "AliVEvent.h"
#include "AliAODEvent.h"
#include "AliVTrack.h"
#include "AliAODTrack.h"
#include "AliAODHandler.h"
#include "AliAODInputHandler.h"
#include "AliAODMCParticle.h"
#include "AliAODMCHeader.h"
#include "AliMCParticle.h"
#include "AliMCEvent.h"
#include "AliStack.h"
#include "AliInputEventHandler.h"
#include "AliEventPoolManager.h"

#include "AliAnalysisTaskDiJetCorrelations.h"

using std::cout;
using std::endl;

ClassImp(AliAnalysisTaskDiJetCorrelations)

//_____________________| Constructor
AliAnalysisTaskDiJetCorrelations::AliAnalysisTaskDiJetCorrelations():
  AliAnalysisTaskSE(),
  fSetSystemValue(kTRUE),
  fRecoOrMontecarlo(kTRUE),
  fReadMC(kFALSE),
  fSetFilterBit(kTRUE),
  fbit(272),
  farrayMC(0),
  fCentrOrMult(1),
  fMinCentrality(0),
  fMaxCentrality(100),
  fTrigger1pTLowThr(0),
  fTrigger1pTHighThr(0),
  fTrigger2pTLowThr(0),
  fTrigger2pTHighThr(0),
  fCutResonances(kTRUE),
  fCutConversions(kTRUE),
  ftwoTrackEfficiencyCut(kTRUE),
  fuseVarCentBins(kFALSE),
  fuseVarPtBins(kFALSE),
  fHistNEvents(0),
  fHistT1CorrTrack(0),
  fHistT2CorrTrack(0),
  fOutputQA(0),
  fOutputCorr(0),
  f3DEffCor(0),
  fPool(0x0),
  fPoolMgr(0x0),
  fMixedEvent(kTRUE),
  fMEMaxPoolEvent(1000),
  fMEMinTracks(20000),
  fMEMinEventToMix(6),
  fHistTrigDPhi(0x0),
  fControlConvResT1(0x0),
  fControlConvResT2(0x0),
  fControlConvResMT1(0x0),
  fControlConvResMT2(0x0)
{
  for ( Int_t i = 0; i < 9; i++)fHistQA[i] = NULL;
}

//_____________________| Specific Constructor
AliAnalysisTaskDiJetCorrelations::AliAnalysisTaskDiJetCorrelations(const char *name):
  AliAnalysisTaskSE(name),
  fSetSystemValue(kTRUE),
  fRecoOrMontecarlo(kTRUE),
  fReadMC(kFALSE),
  fSetFilterBit(kTRUE),
  fbit(272),
  farrayMC(0),
  fCentrOrMult(1),
  fMinCentrality(0),
  fMaxCentrality(100),
  fTrigger1pTLowThr(0),
  fTrigger1pTHighThr(0),
  fTrigger2pTLowThr(0),
  fTrigger2pTHighThr(0),
  fCutResonances(kTRUE),
  fCutConversions(kTRUE),
  ftwoTrackEfficiencyCut(kTRUE),
  fuseVarCentBins(kFALSE),
  fuseVarPtBins(kFALSE),
  fHistNEvents(0),
  fHistT1CorrTrack(0),
  fHistT2CorrTrack(0),
  fOutputQA(0),
  fOutputCorr(0),
  f3DEffCor(0),
  fPool(0x0),
  fPoolMgr(0x0),
  fMixedEvent(kTRUE),
  fMEMaxPoolEvent(1000),
  fMEMinTracks(2000),
  fMEMinEventToMix(6),
  fHistTrigDPhi(0x0),
  fControlConvResT1(0x0),
  fControlConvResT2(0x0),
  fControlConvResMT1(0x0),
  fControlConvResMT2(0x0)
{
  Info("AliAnalysisTaskDiJetCorrelations","Calling Constructor");
  for ( Int_t i = 0; i < 9; i++)fHistQA[i] = NULL;
  DefineInput(0, TChain::Class());
  DefineOutput(1, TList::Class()); // Basic output slot (..)
  DefineOutput(2, TList::Class()); // Correlations form Data and MC    
}

//___________________________________| Copy Constructor
AliAnalysisTaskDiJetCorrelations::AliAnalysisTaskDiJetCorrelations(const AliAnalysisTaskDiJetCorrelations &source):
  AliAnalysisTaskSE(source),
  fSetSystemValue(source.fSetSystemValue),
  fRecoOrMontecarlo(source.fSetSystemValue),
  fReadMC(source.fReadMC),
  fSetFilterBit(source.fSetFilterBit),
  fbit(source.fbit),
  farrayMC(source.farrayMC),
  fCentrOrMult(source.fCentrOrMult),
  fMinCentrality(source.fMinCentrality),
  fMaxCentrality(source.fMaxCentrality),
  fTrigger1pTLowThr(source.fTrigger1pTLowThr),
  fTrigger1pTHighThr(source.fTrigger1pTHighThr),
  fTrigger2pTLowThr(source.fTrigger2pTLowThr),
  fTrigger2pTHighThr(source.fTrigger2pTHighThr),
  fCutResonances(source.fCutResonances),
  fCutConversions(source.fCutConversions),
  ftwoTrackEfficiencyCut(source.ftwoTrackEfficiencyCut),
  fuseVarCentBins(source.fuseVarCentBins),
  fuseVarPtBins(source.fuseVarPtBins),
  fHistNEvents(source.fHistNEvents),
  fHistT1CorrTrack(source.fHistT1CorrTrack),
  fHistT2CorrTrack(source.fHistT2CorrTrack),
  fOutputQA(source.fOutputQA),
  fOutputCorr(source.fOutputCorr),
  f3DEffCor(source.f3DEffCor),
  fPool(source.fPool),
  fPoolMgr(source.fPoolMgr),
  fMixedEvent(source.fMixedEvent),
  fMEMaxPoolEvent(source.fMEMaxPoolEvent),
  fMEMinTracks(source.fMEMinTracks),
  fMEMinEventToMix(source.fMEMinEventToMix),
  fHistTrigDPhi(source.fHistTrigDPhi),
  fControlConvResT1(source.fControlConvResT1),
  fControlConvResT2(source.fControlConvResT2),
  fControlConvResMT1(source.fControlConvResMT1),
  fControlConvResMT2(source.fControlConvResMT2)
{
  
}

//_____________________| Destructor
AliAnalysisTaskDiJetCorrelations::~AliAnalysisTaskDiJetCorrelations()
{
  if(fOutputQA) {delete fOutputQA; fOutputQA = 0;}
  if(fOutputCorr) {delete fOutputCorr; fOutputCorr = 0;}
  if(fHistNEvents) {delete fHistNEvents; fHistNEvents = 0;}
}

//________________________________________|  Assignment Constructor
AliAnalysisTaskDiJetCorrelations& AliAnalysisTaskDiJetCorrelations::operator=(const AliAnalysisTaskDiJetCorrelations& orig)
{
  if (&orig == this) return *this; 
  AliAnalysisTaskSE::operator=(orig); 

  fSetSystemValue= orig.fSetSystemValue;
  fRecoOrMontecarlo = orig.fRecoOrMontecarlo;
  fReadMC = orig.fReadMC;
  fSetFilterBit = orig.fSetFilterBit;
  fbit  = orig.fbit;
  farrayMC = orig.farrayMC;
  fCentrOrMult = orig.fCentrOrMult;
  fMinCentrality = orig.fMinCentrality;
  fMaxCentrality = orig.fMaxCentrality;
  fTrigger1pTLowThr = orig.fTrigger1pTLowThr;
  fTrigger1pTHighThr = orig.fTrigger1pTHighThr;
  fTrigger2pTLowThr = orig.fTrigger2pTLowThr;
  fTrigger2pTHighThr = orig.fTrigger2pTHighThr;
  fCutResonances = orig.fCutResonances;
  fCutConversions = orig.fCutConversions;
  ftwoTrackEfficiencyCut = orig.ftwoTrackEfficiencyCut;
  fuseVarCentBins = orig.fuseVarCentBins;
  fuseVarPtBins = orig.fuseVarPtBins;
  fHistNEvents = orig.fHistNEvents;
  fHistT1CorrTrack = orig.fHistT1CorrTrack;
  fHistT2CorrTrack = orig.fHistT2CorrTrack;
  fOutputQA = orig.fOutputQA;
  fOutputCorr = orig.fOutputCorr;
  f3DEffCor = orig.f3DEffCor;
  fPool = orig.fPool;
  fPoolMgr = orig.fPoolMgr;
  fMixedEvent = orig.fMixedEvent;
  fMEMaxPoolEvent = orig.fMEMinTracks;
  fMEMinTracks = orig.fMEMinTracks;
  fMEMinEventToMix = orig.fMEMinEventToMix;
  return *this; 
}

//_____________________| UserCreate Output
void AliAnalysisTaskDiJetCorrelations::UserCreateOutputObjects()
{
  fOutputQA = new TList();
  fOutputQA->SetOwner();
  fOutputQA->SetName("BasicQAHistograms");
  
  fOutputCorr = new TList();
  fOutputCorr->SetOwner();
  fOutputCorr->SetName("CorrelationsHistograms");

  fHistNEvents = new TH1F("fHistNEvents", "number of events ", 10, -0.5, 9.5);
  fHistNEvents->GetXaxis()->SetBinLabel(1,"nEvents analyzed");
  fHistNEvents->GetXaxis()->SetBinLabel(2,"Rejected due Null Vetex and B");
  fHistNEvents->GetXaxis()->SetBinLabel(3,"Within choosen centrality");
  fHistNEvents->GetXaxis()->SetBinLabel(4,"With Good Vertex");
  fHistNEvents->GetXaxis()->SetBinLabel(5,"With Good SPDVertex");
  fHistNEvents->GetXaxis()->SetBinLabel(6,"nTrack all choosen");
  fHistNEvents->GetXaxis()->SetBinLabel(7,"nTrack Trigger1");
  fHistNEvents->GetXaxis()->SetBinLabel(8,"nTrack Trigger2");
  fHistNEvents->GetXaxis()->SetBinLabel(9,"nEvents with T1");
  fHistNEvents->GetXaxis()->SetBinLabel(10,"number of DiJet");
  fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
  fHistNEvents->Sumw2();
  fHistNEvents->SetMinimum(0);
  fOutputQA->Add(fHistNEvents);
  
  fHistT1CorrTrack = new TH1F("fHistT1CorrTrack", "T1 nCorr tracks ", 5, -0.5, 4.5);
  fHistT1CorrTrack->GetXaxis()->SetBinLabel(1,"nTrack Total");
  fHistT1CorrTrack->GetXaxis()->SetBinLabel(2,"after conversion");
  fHistT1CorrTrack->GetXaxis()->SetBinLabel(3,"after K0 resonance");
  fHistT1CorrTrack->GetXaxis()->SetBinLabel(4,"after Lambda resonance");
  fHistT1CorrTrack->GetXaxis()->SetBinLabel(5,"after 2track eff");
  fHistT1CorrTrack->GetXaxis()->SetNdivisions(1,kFALSE);
  fHistT1CorrTrack->Sumw2();
  fHistT1CorrTrack->SetMinimum(0);
  fOutputQA->Add(fHistT1CorrTrack);
  
  fHistT2CorrTrack = new TH1F("fHistT2CorrTrack", "T2 nCorr tracks ", 5, -0.5, 4.5);
  fHistT2CorrTrack->GetXaxis()->SetBinLabel(1,"nTrack Total");
  fHistT2CorrTrack->GetXaxis()->SetBinLabel(2,"after conversion");
  fHistT2CorrTrack->GetXaxis()->SetBinLabel(3,"after K0 resonance");
  fHistT2CorrTrack->GetXaxis()->SetBinLabel(4,"after Lambda resonance");
  fHistT2CorrTrack->GetXaxis()->SetBinLabel(5,"after 2track eff");
  fHistT2CorrTrack->GetXaxis()->SetNdivisions(1,kFALSE);
  fHistT2CorrTrack->Sumw2();
  fHistT2CorrTrack->SetMinimum(0);
  fOutputQA->Add(fHistT2CorrTrack);
  
  DefineHistoNames();
  
  Bool_t DefPool = DefineMixedEventPool();
  if(!DefPool){
    AliInfo("UserCreateOutput: Pool is not define properly");
    return;
  }
  
  PostData(1,fOutputQA);
  PostData(2,fOutputCorr);
}


//_____________________| User Exec Part
void  AliAnalysisTaskDiJetCorrelations::UserExec(Option_t *) 
{

    AliAODEvent *aod = dynamic_cast<AliAODEvent*>(InputEvent());
    if(!aod && AODEvent() && IsStandardAOD()) {
        aod = dynamic_cast<AliAODEvent*> (AODEvent());
    } else if(!aod)  {
        printf("AliAnalysisTaskDiJetCorrelations::UserExec: AOD not found!\n");
        return;
    }
  
  fHistNEvents->Fill(0);
  
  if(!fRecoOrMontecarlo){ // MC Kine
    farrayMC = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
    if(!farrayMC){
      AliError("Array of MC particles not found");
      return;
    }
    AliAODMCHeader *mcHeader = NULL;
    mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
    if(!mcHeader) {
      printf("AliAnalysisTaskDiJetCorrelations::UserExec: MC header branch not found!\n");
      return;
    }
  }
  
  if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
  Float_t bSign = 0;
  bSign = (aod->GetMagneticField() > 0) ? 1 : -1;
  fHistNEvents->Fill(1);
  
  AliCentrality *centralityObj = 0x0;
  if(fSetSystemValue){ // pPb, PbPb
    centralityObj = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP();
    fCentrOrMult = centralityObj->GetCentralityPercentileUnchecked("V0M");
    if(centralityObj->GetQuality()!=0) return;
    if((abs(fCentrOrMult)) < 0. || (abs(fCentrOrMult)) > 100.1)return;
  }
  else if(!fSetSystemValue){ // pp, pPb
    Double_t count = -1, mineta = -1.0, maxeta = 1.0;
    AliAODTracklets* tracklets = aod->GetTracklets();
    Int_t nTr=tracklets->GetNumberOfTracklets();
    for(Int_t iTr=0; iTr<nTr; iTr++){
      Double_t theta=tracklets->GetTheta(iTr);
      Double_t eta=-TMath::Log(TMath::Tan(theta/2.));
      if(eta>mineta && eta<maxeta) count++;
    }
    fCentrOrMult = count;
  }
  fHistNEvents->Fill(2); //
  
  AliAODVertex *vtxPrm = (AliAODVertex*)aod->GetPrimaryVertex();
  Bool_t isGoodVtx=kFALSE;
  TString primTitle = vtxPrm->GetTitle();
  if(primTitle.Contains("VertexerTracks") && vtxPrm->GetNContributors()>0) {
    isGoodVtx=kTRUE;
    fHistNEvents->Fill(3);
  }
  if(!isGoodVtx) return;
  
  Float_t zVertex = vtxPrm->GetZ();
  if (TMath::Abs(zVertex) > 10) return;
  fHistQA[0]->Fill(zVertex);
  
  const AliAODVertex* vtxSPD = aod->GetPrimaryVertexSPD();
  if(vtxSPD->GetNContributors()<=0) return;
  TString vtxTyp = vtxSPD->GetTitle();
  Double_t cov[6] = {0};
  vtxSPD->GetCovarianceMatrix(cov);
  Double_t zRes = TMath::Sqrt(cov[5]);
  if(vtxTyp.Contains("vertexer:Z") && (zRes > 0.25)) return;
  if(TMath::Abs(vtxSPD->GetZ() - vtxPrm->GetZ()) > 0.5) return;
  fHistNEvents->Fill(4);
   
  TObjArray* fTrackArray = new TObjArray;
  fTrackArray->SetOwner(kTRUE);
  
  Double_t refmaxpT1 = -999, refmaxpT2 = -999;
  Double_t etaMaxpT1 = -999, etaMaxpT2 = -999;
  Double_t phiMaxpT1 = -999, phiMaxpT2 = -999;
  Short_t Charge1, Charge2;
  
  TString typeData = "";
  TString SEorME = "";
  
  if(fRecoOrMontecarlo){
    if(!fReadMC){
      typeData += "Data";//data
    }else
      typeData += "MCrc"; // MC reco
  }else{
    typeData += "MCKn"; // MC Generations  
  }
  
  if(!fMixedEvent){
    SEorME += "SE";
  }else if(fMixedEvent){
    SEorME += "ME";
  }
  
  //Mixed Events..
  if(TMath::Abs(zVertex)>=10)AliInfo(Form("Event with Zvertex = %.2f cm out of pool bounds, SKIPPING",zVertex));
  if(fMixedEvent){
    fPool= fPoolMgr->GetEventPool(fCentrOrMult, zVertex);
    if(!fPool){
      AliInfo(Form("No pool found for Event: multiplicity = %f, zVtx = %f cm", fCentrOrMult, zVertex));
      return;
    }
  }
  
  
  for (Int_t iTracks = 0; iTracks < aod->GetNumberOfTracks(); iTracks++){
    
    AliAODTrack* fAodTracks = (AliAODTrack*)aod->GetTrack(iTracks);
    if (!fAodTracks)continue;
    if(fSetFilterBit) if (!fAodTracks->TestFilterBit(fbit)) continue;
    
    if(fAodTracks->Eta() < -0.8 || fAodTracks->Eta() > 0.8)continue;
    if (fAodTracks->Pt() < 0.5 || fAodTracks->Pt() > 30.)continue;
    
    fHistNEvents->Fill(5); 
    fHistQA[1]->Fill(fAodTracks->GetTPCClusterInfo(2,1));
    fHistQA[3]->Fill(fAodTracks->DCA());
    fHistQA[4]->Fill(fAodTracks->ZAtDCA());
    fHistQA[5]->Fill(fAodTracks->Chi2perNDF());
    fHistQA[6]->Fill(fAodTracks->Pt());
    fHistQA[7]->Fill(fAodTracks->Phi());
    fHistQA[8]->Fill(fAodTracks->Eta());
    
    if(fRecoOrMontecarlo){ // reconstruction of data and MC
      if(fReadMC){
	// is track associated to particle ? if yes + implimenting the physical primary..
	Int_t label = TMath::Abs(fAodTracks->GetLabel());
	if (label<=0){
	  AliDebug(3,"Particle not matching MC label \n");
	  continue;
	}
	
	AliAODMCParticle *mcPart  = (AliAODMCParticle*)fMCEvent->GetTrack(label);
	if (!mcPart->IsPhysicalPrimary()) continue;
	fTrackArray->Add(fAodTracks);
      }else
	fTrackArray->Add(fAodTracks); //Storing all tracks for Data
    }
    
    if(fAodTracks->Pt() >= fTrigger1pTLowThr && fAodTracks->Pt()  <= fTrigger1pTHighThr){  
      if(fAodTracks->Pt() > refmaxpT1){
	refmaxpT1 = fAodTracks->Pt();
	etaMaxpT1 = fAodTracks->Eta();
	phiMaxpT1 = fAodTracks->Phi();
	Charge1 = fAodTracks->Charge();
	fHistNEvents->Fill(6);
      }
    }
  }
  
  if(!fMixedEvent)if(refmaxpT1 < fTrigger1pTLowThr || refmaxpT1 > fTrigger1pTHighThr)return;
  fHistNEvents->Fill(8);
  
  for(Int_t entryT2=0; entryT2<fTrackArray->GetEntries(); entryT2++){
    
    TObject* obj = fTrackArray->At(entryT2);
    AliAODTrack* fAodTracksT2 = (AliAODTrack*)obj;
    
    if(fAodTracksT2->Pt() >= fTrigger2pTLowThr && fAodTracksT2->Pt()  <= fTrigger2pTHighThr){
      
      Double_t phiT2 = fAodTracksT2->Phi();
      Double_t TrigDPhi = phiMaxpT1-phiT2;
      
      Double_t TrigDPhi12 =  AssignCorrectPhiRange(TrigDPhi);
      Double_t B2BTrigDPhi = TMath::Abs(TrigDPhi12-TMath::Pi());
      
      if(B2BTrigDPhi > (TMath::Pi())/8)continue;
      fHistTrigDPhi->Fill(B2BTrigDPhi);
      
      if(fAodTracksT2->Pt() > refmaxpT2 ){
	fHistNEvents->Fill(7);
	refmaxpT2 = fAodTracksT2->Pt();
	etaMaxpT2 = fAodTracksT2->Eta();
	phiMaxpT2 = fAodTracksT2->Phi();
	Charge2 = fAodTracksT2->Charge();      
      }
    }
  }
  
  if(!fMixedEvent)if(refmaxpT2 < fTrigger2pTLowThr || refmaxpT2 > fTrigger2pTHighThr)return;
  fHistNEvents->Fill(9);
  
  if (refmaxpT1 != -999 && refmaxpT2 != -999){
    
    Double_t fCentZvtxpT1[3] = {refmaxpT1, fCentrOrMult, zVertex};
    if(!fMixedEvent)((THnSparseD*)fOutputCorr->FindObject(Form("ThnTrg1CentZvtxpT_%s_%s",typeData.Data(), SEorME.Data())))->Fill(fCentZvtxpT1);
    
    Double_t fCentZvtxpT2[3] = {refmaxpT2, fCentrOrMult, zVertex};
    if(!fMixedEvent)((THnSparseD*)fOutputCorr->FindObject(Form("ThnTrg2CentZvtxpT_%s_%s",typeData.Data(), SEorME.Data())))->Fill(fCentZvtxpT2);
  }
  
  
  Int_t NofEventsinPool = 1; // SE
  if(TMath::Abs(zVertex)>=10)AliInfo(Form("Event with Zvertex = %.2f cm out of pool bounds, SKIPPING",zVertex));
  
  if(fMixedEvent){
    //Bool_t execPool = ProcessMixedEventPool();
    NofEventsinPool = fPool->GetCurrentNEvents();
  }
  
  for (Int_t jMix =0; jMix < NofEventsinPool; jMix++){
    
    if(!fTrackArray) return;
    if(fMixedEvent)fTrackArray = fPool->GetEvent(jMix); //check twice this coding wise..
    Int_t NumberOfTracksStore = fTrackArray->GetEntriesFast();
    
    for(int k=0; k < NumberOfTracksStore; k++){
      
      TObject* obj2 = fTrackArray->At(k);
      AliAODTrack* fAodTracksAS = (AliAODTrack*) obj2;
      if(!fAodTracksAS) continue;
      if (fAodTracksAS->Pt() > refmaxpT1)continue;
      
      if(fMixedEvent){
	if(fAodTracksAS->Eta() < -0.8 || fAodTracksAS->Eta() > 0.8)continue;
	if(fAodTracksAS->Pt() < 0.5 || fAodTracksAS->Pt() > 20)continue;
	if(fSetFilterBit) if (!fAodTracksAS->TestFilterBit(fbit)) continue;
      }
      
      if(fAodTracksAS->Pt() < refmaxpT1){
	
	Bool_t CutForConversionResonanceTrg1 = ConversionResonanceCut(refmaxpT1, phiMaxpT1, etaMaxpT1, Charge1, fAodTracksAS,fControlConvResT1, fHistT1CorrTrack);
	if(fCutConversions || fCutResonances)if(!CutForConversionResonanceTrg1)continue;
	
	Bool_t CutForTwoTrackEffiTrg1 = TwoTrackEfficiencyCut(refmaxpT1, phiMaxpT1, etaMaxpT1, Charge1, fAodTracksAS, bSign);
	if(ftwoTrackEfficiencyCut)if(!CutForTwoTrackEffiTrg1)continue;
	fHistT1CorrTrack->Fill(4); //
	
	Double_t deltaPhi1  = phiMaxpT1 - fAodTracksAS->Phi();
	Double_t deltaPhi1C = AssignCorrectPhiRange(deltaPhi1);
	Double_t deltaEta1  = etaMaxpT1 - fAodTracksAS->Eta();
        
	Double_t ptLim_Sparse1 = 0.0 ;
	Double_t ptTrack1 = fAodTracksAS->Pt();
	ptLim_Sparse1 =((THnSparseD*)fOutputCorr->FindObject(Form("ThnTrg1CentZvtxDEtaDPhi_%s_%s",typeData.Data(), SEorME.Data())))->GetAxis(4)->GetXmax();
	if(ptTrack1 > ptLim_Sparse1) ptTrack1 = ptLim_Sparse1-0.01; //filling all the pT in last bin
        
	Double_t CentzVtxDEtaDPhiTrg1[5] = {fCentrOrMult, zVertex, deltaEta1, deltaPhi1C, ptTrack1};
	Double_t effvalueT1 = 1.0;//GetTrackbyTrackEffValue(fAodTracksAS, fCentrOrMult);
	
	((THnSparseD*)fOutputCorr->FindObject(Form("ThnTrg1CentZvtxDEtaDPhi_%s_%s",typeData.Data(), SEorME.Data())))->Fill(CentzVtxDEtaDPhiTrg1, 1/effvalueT1);
	
	//Double_t fBasicsTrg1[5] = {fAodTracksAS->Pt(), fAodTracksAS->Eta(), fCentrOrMult, zVertex, fAodTracksAS->Phi()};
	//if(!fMixedEvent)((THnSparseD*)fOutputQA->FindObject(Form("ThnTrg1BasicsPlots_%s_%s",typeData.Data(), SEorME.Data())))->Fill(fBasicsTrg1);
        
      }
      
      
      if (fAodTracksAS->Pt() < refmaxpT2){
	
	Bool_t CutForConversionResonanceTrg2 = ConversionResonanceCut(refmaxpT2, phiMaxpT2, etaMaxpT2, Charge2, fAodTracksAS,fControlConvResT2, fHistT2CorrTrack);
	if(fCutConversions || fCutResonances)if(!CutForConversionResonanceTrg2)continue;
	
	Bool_t CutForTwoTrackEffiTrg2 = TwoTrackEfficiencyCut(refmaxpT2, phiMaxpT2, etaMaxpT2, Charge2, fAodTracksAS, bSign);
	if(ftwoTrackEfficiencyCut)if(!CutForTwoTrackEffiTrg2)continue;	
	fHistT2CorrTrack->Fill(4); //
	
	Double_t deltaPhi2 = phiMaxpT2 - fAodTracksAS->Phi();
	Double_t deltaPhi2C =  AssignCorrectPhiRange(deltaPhi2);
	Double_t deltaEta2 = etaMaxpT2 - fAodTracksAS->Eta();
	
	Double_t ptLim_Sparse2 = 0.0;
	Double_t ptTrack2 = fAodTracksAS->Pt();
	ptLim_Sparse2 =((THnSparseD*)fOutputCorr->FindObject(Form("ThnTrg2CentZvtxDEtaDPhi_%s_%s",typeData.Data(), SEorME.Data())))->GetAxis(4)->GetXmax();
	if(ptTrack2 > ptLim_Sparse2) ptTrack2 = ptLim_Sparse2-0.01; //filling all the pT in last bin
	
	Double_t CentzVtxDEtaDPhiTrg2[5] = {fCentrOrMult, zVertex, deltaEta2, deltaPhi2C,ptTrack2};
	Double_t effvalueT2 = 1.0;//GetTrackbyTrackEffValue(fAodTracksAS, fCentrOrMult);
	
	((THnSparseD*)fOutputCorr->FindObject(Form("ThnTrg2CentZvtxDEtaDPhi_%s_%s",typeData.Data(), SEorME.Data())))->Fill(CentzVtxDEtaDPhiTrg2,1/effvalueT2);
	
	// Double_t fBasicsTrg2[5] = {fAodTracksAS->Pt(), fAodTracksAS->Eta(), fCentrOrMult, zVertex, fAodTracksAS->Phi()};
	// if(!fMixedEvent)((THnSparseD*)fOutputQA->FindObject(Form("ThnTrg2BasicsPlots_%s_%s",typeData.Data(), SEorME.Data())))->Fill(fBasicsTrg2);
	
      }
    }
  }// event looop
  
  // clone and give ownership to event pool
  TObjArray* tracksClone = (TObjArray*) fTrackArray->Clone();
  tracksClone->SetOwner(kTRUE);
  if(fMixedEvent && tracksClone->GetEntriesFast()>0)fPool->UpdatePool(tracksClone);
  
  PostData(1, fOutputQA);
  PostData(2, fOutputCorr);  
  //delete f3DEffCor;
  
}

//_____________________|Terminate
void AliAnalysisTaskDiJetCorrelations::Terminate(Option_t *){
  
  fOutputQA = dynamic_cast<TList*> (GetOutputData(1));
  if (!fOutputQA) {
    printf("ERROR: Output list not available\n");
    return;
  }
  
  fOutputCorr = dynamic_cast<TList*> (GetOutputData(2));
  if (!fOutputCorr) {
    printf("ERROR: fOutputCorr not available\n");
    return;
  }
  return;
}

//______________________________|  Cuts of Resonance and Conversions..
Bool_t AliAnalysisTaskDiJetCorrelations::ConversionResonanceCut(Double_t refmaxpT, Double_t phiMaxpT, Double_t etaMaxpT, Double_t Charge, AliAODTrack* AodTracks, TH2F*fControlConvResT, TH1F* fHistTCorrTrack){
  
  fHistTCorrTrack->Fill(0); //
  //Conversions
  if(fCutConversions && AodTracks->Charge() * Charge < 0){
    Float_t mass = GetInvMassSquaredCheap(refmaxpT, etaMaxpT, phiMaxpT, AodTracks->Pt(), AodTracks->Eta(), AodTracks->Phi(), 0.510e-3, 0.510e-3); 
    if (mass < 0.1){
      mass = GetInvMassSquared(refmaxpT, etaMaxpT, phiMaxpT, AodTracks->Pt(), AodTracks->Eta(), AodTracks->Phi(), 0.510e-3, 0.510e-3);
      fControlConvResT->Fill(0.0, mass);
      if(mass < 0.04 * 0.04)	return kFALSE;
    }
  }
  fHistTCorrTrack->Fill(1); //
  
  //K0s
  if (fCutResonances && AodTracks->Charge() * Charge < 0){   
    Float_t mass = GetInvMassSquaredCheap(refmaxpT, etaMaxpT, phiMaxpT, AodTracks->Pt(), AodTracks->Eta(), AodTracks->Phi(), 0.1396, 0.1396);
    const Float_t kK0smass = 0.4976;
    if (TMath::Abs(mass -kK0smass * kK0smass)  < 0.1){
      mass = GetInvMassSquared(refmaxpT, etaMaxpT, phiMaxpT, AodTracks->Pt(), AodTracks->Eta(), AodTracks->Phi(), 0.1396, 0.1396);
      fControlConvResT->Fill(1, mass -kK0smass * kK0smass);
      if(mass > (kK0smass-0.02)*(kK0smass-0.02) && mass < (kK0smass+0.02)*(kK0smass+0.02))return kFALSE;
    }
  }
  fHistTCorrTrack->Fill(2); //
  
  //lambda
  if (fCutResonances && AodTracks->Charge() * Charge < 0){
    Float_t mass1 = GetInvMassSquaredCheap(refmaxpT, etaMaxpT, phiMaxpT, AodTracks->Pt(), AodTracks->Eta(), AodTracks->Phi(), 0.1396, 0.9383);
    Float_t mass2 = GetInvMassSquaredCheap(refmaxpT, etaMaxpT, phiMaxpT, AodTracks->Pt(), AodTracks->Eta(), AodTracks->Phi(), 0.9383, 0.1396);
    const Float_t kLambdaMass = 1.115;
    if (TMath::Abs(mass1 -kLambdaMass * kLambdaMass)  < 0.1){
      mass1 = GetInvMassSquared(refmaxpT, etaMaxpT, phiMaxpT, AodTracks->Pt(), AodTracks->Eta(), AodTracks->Phi(), 0.1396, 0.9383);
      fControlConvResT->Fill(2, mass1 - kLambdaMass * kLambdaMass);
      if(mass1 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass1 < (kLambdaMass+0.02)*(kLambdaMass+0.02))return kFALSE;  
    }
    
    if (TMath::Abs(mass2 -kLambdaMass * kLambdaMass)  < 0.1){
      mass2 = GetInvMassSquared(refmaxpT, etaMaxpT, phiMaxpT, AodTracks->Pt(), AodTracks->Eta(), AodTracks->Phi(), 0.1396, 0.9383);
      fControlConvResT->Fill(2, mass2 - kLambdaMass * kLambdaMass);
      if(mass2 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass2 < (kLambdaMass+0.02)*(kLambdaMass+0.02))return kFALSE;
    }
  }
  
  fHistTCorrTrack->Fill(3); //
  return kTRUE; 
}

//_____________________|  Two track Efficiency
Bool_t AliAnalysisTaskDiJetCorrelations::TwoTrackEfficiencyCut(Double_t refmaxpT, Double_t phiMaxpT, Double_t etaMaxpT, Double_t Charge, AliAODTrack* AodTracks, Float_t bSigntmp){
  Float_t pt1  = refmaxpT;
  Float_t phi1 = phiMaxpT;
  Float_t eta1 = etaMaxpT;
  Float_t charge1 = Charge;
  
  Float_t phi2 = AodTracks->Phi();
  Float_t eta2 = AodTracks->Eta();
  Float_t pt2  = AodTracks->Pt();
  Float_t charge2 = AodTracks->Charge();
  
  Float_t deta = eta1 - eta2;
  
  if (TMath::Abs(deta) < 0.02 * 2.5 * 3){
    Float_t dphistar1 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 0.8, bSigntmp);
    Float_t dphistar2 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 2.5, bSigntmp);
    
    const Float_t kLimit = 0.02*3;
    Float_t dphistarminabs = 1e5;
    //Float_t dphistarmin = 1e5;
  
    if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0){
      for (Double_t rad=0.8; rad<2.51; rad+=0.01){
	Float_t dphistar = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, rad, bSigntmp);
	Float_t dphistarabs = TMath::Abs(dphistar);
	if (dphistarabs < dphistarminabs){
	  //dphistarmin = dphistar;
	  dphistarminabs = dphistarabs;
	}
      }
       
      if (dphistarminabs < 0.02 && TMath::Abs(deta) < 0.02){
	return kFALSE;
      }
    }
  }
  return kTRUE;
}


//______________________________|  Nomenclature of Histograms
void AliAnalysisTaskDiJetCorrelations::DefineHistoNames(){
  
  Double_t Pi = TMath::Pi();
  //QA histograms
  fHistQA[0] = new TH1F("fHistZVtx", "Z vertex distribution", 100, -15., 15.);
  fHistQA[1] = new TH1F("fHistnTPCCluster", "n TPC Cluster", 100, 0., 200.);
  fHistQA[2] = new TH1F("fHistnTPCClusterF", "n TPC Cluster findable", 100, 0., 200.);
  fHistQA[3] = new TH1F("fHistDCAXY", "dca-XY", 100, -5., 5.);
  fHistQA[4] = new TH1F("fHistDCAZ", "dca-Z", 100, -5., 5.);
  fHistQA[5] = new TH1F("fHistChi2TPC", "Chi2 TPC", 100, 0., 10.);
  fHistQA[6] = new TH1F("fHistpT", "pT distribution",1000,0.,20.);
  fHistQA[7] = new TH1F("fHistPhi", "Phi distribution" , 100, -0.5, 2*Pi+0.5);
  fHistQA[8] = new TH1F("fHistEta", "Eta distribution" , 100, -2, 2);
  
  for( Int_t i = 0; i < 9; i++)
    {
      fHistQA[i]->Sumw2();
      fOutputQA->Add(fHistQA[i]);
    }
  
  //1 SE Distributuons Phi, Eta for Trg1 and Trg2
  fHistTrigDPhi = new TH1F("fHistTrigDPhi", " Trig Phi Difference Same",100, 0, 0.5);
  fHistTrigDPhi->Sumw2();
  fOutputQA->Add(fHistTrigDPhi);
  
  fControlConvResT1 = new TH2F("fControlConvResT1", ";id;delta mass;T1", 3, -0.5, 2.5, 100, -0.1, 0.1);
  fControlConvResT2 = new TH2F("fControlConvResT2", ";id;delta mass;T2", 3, -0.5, 2.5, 100, -0.1, 0.1);
  fControlConvResT1->Sumw2();
  fControlConvResT2->Sumw2();
  fOutputQA->Add(fControlConvResT1);
  fOutputQA->Add(fControlConvResT2);
  
  //Thnsprase Distributuons Phi, Eta for Trg1 and Trg2
  TString nameThnTrg1CentZvtxDEtaDPhi   = "ThnTrg1CentZvtxDEtaDPhi";
  TString nameThnTrg2CentZvtxDEtaDPhi   = "ThnTrg2CentZvtxDEtaDPhi";
  TString nameThnTrg1CentZvtxpT     = "ThnTrg1CentZvtxpT" ;
  TString nameThnTrg2CentZvtxpT    = "ThnTrg2CentZvtxpT" ;
  TString nameThnAllBasicsPlots    = "ThnAllBasicsPlots" ;
  TString nameThnTrg1BasicsPlots    = "ThnTrg1BasicsPlots";
  TString nameThnTrg2BasicsPlots    = "ThnTrg2BasicsPlots";
  
  if(fRecoOrMontecarlo){
    if(!fReadMC){ //data
      nameThnTrg1CentZvtxDEtaDPhi  += "_Data";
      nameThnTrg2CentZvtxDEtaDPhi  += "_Data";
      nameThnTrg1CentZvtxpT  += "_Data";
      nameThnTrg2CentZvtxpT  += "_Data";
      nameThnAllBasicsPlots  += "_Data";
      nameThnTrg1BasicsPlots  += "_Data";
      nameThnTrg2BasicsPlots  += "_Data";  
    }
    else {// MC reco
      nameThnTrg1CentZvtxDEtaDPhi  += "_MCrc";
      nameThnTrg2CentZvtxDEtaDPhi  += "_MCrc";
      nameThnTrg1CentZvtxpT  += "_MCrc";
      nameThnTrg2CentZvtxpT  += "_MCrc";
      nameThnAllBasicsPlots  += "_MCrc";
      nameThnTrg1BasicsPlots  += "_MCrc";
      nameThnTrg2BasicsPlots  += "_MCrc";
    }
  }
  else{ // MC Generations
    nameThnTrg1CentZvtxDEtaDPhi  += "_MCKn";
    nameThnTrg2CentZvtxDEtaDPhi  += "_MCKn";
    nameThnTrg1CentZvtxpT  += "_MCKn";
    nameThnTrg2CentZvtxpT  += "_MCKn";
    nameThnAllBasicsPlots  += "_MCKn";
    nameThnTrg1BasicsPlots  += "_MCKn";
    nameThnTrg2BasicsPlots  += "_MCKn";
  }
  
  if(!fMixedEvent){
    nameThnTrg1CentZvtxDEtaDPhi  += "_SE";
    nameThnTrg2CentZvtxDEtaDPhi  += "_SE";
    nameThnTrg1CentZvtxpT  += "_SE";
    nameThnTrg2CentZvtxpT  += "_SE";
    nameThnAllBasicsPlots  += "_SE";
    nameThnTrg1BasicsPlots  += "_SE";
    nameThnTrg2BasicsPlots  += "_SE";
  }else if(fMixedEvent){
    nameThnTrg1CentZvtxDEtaDPhi  += "_ME";
    nameThnTrg2CentZvtxDEtaDPhi  += "_ME";
    nameThnTrg1CentZvtxpT  += "_ME";
    nameThnTrg2CentZvtxpT  += "_ME";
    }

    //For Trigger One --> cent, zvrtx, pT
  Int_t    fBins12[5] = {12,      10,   16,               36, 11};
  Double_t  fMin12[5] = {0.0,   -10., -1.6, -0.5*TMath::Pi(), 0.0};
  Double_t  fMax12[5] = {100.1,  10.,  1.6,  1.5*TMath::Pi(), 6.00};
  
  Int_t pTbinTrigger1 = Int_t(fTrigger1pTHighThr - fTrigger1pTLowThr);
  Int_t pTbinTrigger2 = Int_t(fTrigger2pTHighThr - fTrigger2pTLowThr);
    
  Int_t   fBinsTrg1[3]   = {   12,      10,          pTbinTrigger1};
  Double_t fMinTrg1[3]   = {  0.0,   -10.0,      fTrigger1pTLowThr};
  Double_t fMaxTrg1[3]   = {  100.1,  10.0,   fTrigger1pTHighThr};
  
  Int_t   fBinsTrg2[3]   = {   12,      10,        pTbinTrigger2};
  Double_t fMinTrg2[3]   = {  0.0,   -10.0,    fTrigger2pTLowThr};
  Double_t fMaxTrg2[3]   = {  100.1,    10.0,   fTrigger2pTHighThr};
  
  Int_t   fBinsTrk[5]   = {  20,   16,    10,    10,  36};
  Double_t fMinTrk[5]   = { 0.0, -0.8,   0.0,  -10.,  0};
  Double_t fMaxTrk[5]   = {20.0,  0.8, 100.1,   10.,  2*TMath::Pi()};
    // why the limits of DeltaPhi have been changed
  
  THnSparseD *THnTrig1CentZvtxDEtaDPhi   = new THnSparseD(nameThnTrg1CentZvtxDEtaDPhi.Data(),"Cent-zVtx-DEta1-DPhi1",5, fBins12, fMin12, fMax12);
  THnSparseD *THnTrig2CentZvtxDEtaDPhi   = new THnSparseD(nameThnTrg2CentZvtxDEtaDPhi.Data(),"Cent-zVtx-DEta2-DPhi2",5, fBins12, fMin12, fMax12);
  THnSparseD *THnTrig1CentZvtxpT = new THnSparseD(nameThnTrg1CentZvtxpT.Data(),"Cent-Zvtx-pTtr1",3, fBinsTrg1, fMinTrg1, fMaxTrg1);
  THnSparseD *THnTrig2CentZvtxpT = new THnSparseD(nameThnTrg2CentZvtxpT.Data(),"Cent-Zvtx-pTtr2",3, fBinsTrg2, fMinTrg2, fMaxTrg2);

  if(fuseVarCentBins){
    const Int_t nvarBinsCent = 12;
    Double_t varBinsCent[nvarBinsCent+1] = {0., 1., 2., 3., 4., 5., 7.5, 10., 20., 30., 40., 50., 100.1};
    THnTrig1CentZvtxDEtaDPhi->GetAxis(0)->Set(nvarBinsCent, varBinsCent);
    THnTrig2CentZvtxDEtaDPhi->GetAxis(0)->Set(nvarBinsCent, varBinsCent);  
    THnTrig1CentZvtxpT->GetAxis(0)->Set(nvarBinsCent, varBinsCent);
    THnTrig2CentZvtxpT->GetAxis(0)->Set(nvarBinsCent, varBinsCent);
  }
    if(fuseVarPtBins){
        const Int_t nvarBinspT = 11;
        Double_t varBinspT[nvarBinspT+1] = {0.5, 0.75, 1.0, 1.25, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 5.0, 6.0};
        THnTrig1CentZvtxDEtaDPhi->GetAxis(4)->Set(nvarBinspT, varBinspT);
        THnTrig2CentZvtxDEtaDPhi->GetAxis(4)->Set(nvarBinspT, varBinspT);
        
    }
    
  THnTrig1CentZvtxDEtaDPhi->Sumw2();
  THnTrig2CentZvtxDEtaDPhi->Sumw2();
  THnTrig1CentZvtxpT->Sumw2();
  THnTrig2CentZvtxpT->Sumw2();
  fOutputCorr->Add(THnTrig1CentZvtxDEtaDPhi);
  fOutputCorr->Add(THnTrig2CentZvtxDEtaDPhi);
  fOutputCorr->Add(THnTrig1CentZvtxpT);
  fOutputCorr->Add(THnTrig2CentZvtxpT);
  
  THnSparseD *THnTrg1BasicsPlots = new THnSparseD(nameThnTrg1BasicsPlots.Data(),"pTtr1-Eta-Cent-zvtx-Phi",5, fBinsTrk, fMinTrk, fMaxTrk);
  THnSparseD *THnTrg2BasicsPlots = new THnSparseD(nameThnTrg2BasicsPlots.Data(),"pTtr2-Eta-Cent-zvtx-Phi",5, fBinsTrk, fMinTrk, fMaxTrk);
  THnSparseD *THnAllBasicsPlots = new THnSparseD(nameThnAllBasicsPlots.Data(),"pTtrk-Eta-Cent-zvtx-Phi",5, fBinsTrk, fMinTrk, fMaxTrk);  
  THnTrg1BasicsPlots->Sumw2();
  THnTrg2BasicsPlots->Sumw2();
  THnAllBasicsPlots->Sumw2();
  fOutputQA->Add(THnTrg1BasicsPlots);
  fOutputQA->Add(THnTrg2BasicsPlots);
  fOutputQA->Add(THnAllBasicsPlots);
}
//end of file : Greeshma