// **************************************
// Task used for the systematic study of jet finders
//
// Compares input (gen) and output (rec) jets   
// gen can also be another jet finder on the rec level, matching is done in eta phi
//
// *******************************************


/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

 
#include <TROOT.h>
#include <TRandom.h>
#include <TSystem.h>
#include <TInterpreter.h>
#include <TChain.h>
#include <TFile.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TH3F.h>
#include <TProfile.h>
#include <TList.h>
#include <TKey.h>
#include <TLorentzVector.h>
#include <TClonesArray.h>
#include  "TDatabasePDG.h"

#include "AliAnalysisTaskJFSystematics.h"
#include "AliAnalysisManager.h"
#include "AliJetFinder.h"
#include "AliJetHeader.h"
#include "AliJetReader.h"
#include "AliJetReaderHeader.h"
#include "AliUA1JetHeaderV1.h"
#include "AliESDEvent.h"
#include "AliAODEvent.h"
#include "AliAODHandler.h"
#include "AliAODTrack.h"
#include "AliAODJet.h"
#include "AliMCEventHandler.h"
#include "AliMCEvent.h"
#include "AliStack.h"
#include "AliGenPythiaEventHeader.h"
#include "AliJetKineReaderHeader.h"
#include "AliGenCocktailEventHeader.h"
#include "AliInputEventHandler.h"


#include "AliAnalysisHelperJetTasks.h"

ClassImp(AliAnalysisTaskJFSystematics)

  const Int_t AliAnalysisTaskJFSystematics::fgkSysBins[AliAnalysisTaskJFSystematics::kSysTypes] = {0,AliAnalysisTaskJFSystematics::kMaxJets};
const char* AliAnalysisTaskJFSystematics::fgkSysName[AliAnalysisTaskJFSystematics::kSysTypes] = {"","j"};

AliAnalysisTaskJFSystematics::AliAnalysisTaskJFSystematics(): AliAnalysisTaskSE(),
  fJetHeaderRec(0x0),
  fJetHeaderGen(0x0),
  fAOD(0x0),
  fBranchRec("jets"),
  fBranchGen(""),
  fUseAODInput(kFALSE),
  fUseExternalWeightOnly(kFALSE),
  fLimitGenJetEta(kFALSE),
  fAnalysisType(0),
  fExternalWeight(1),    
  fRecEtaWindow(0.5),
  fAvgTrials(1),
  fh1Xsec(0x0),
  fh1Trials(0x0),
  fh1PtHard(0x0),
  fh1PtHardNoW(0x0),  
  fh1PtHardTrials(0x0),
  fh1NGenJets(0x0),
  fh1NRecJets(0x0),
  fh1PtRecIn(0x0),						     
  fh1PtRecOut(0x0),     
  fh1PtGenIn(0x0),      
  fh1PtGenOut(0x0),     
  fh2PtFGen(0x0),       
  fh2PhiFGen(0x0),      
  fh2EtaFGen(0x0),      
  fh2PtGenDeltaPhi(0x0),
  fh2PtGenDeltaEta(0x0),
  fh3RecInEtaPhiPt(0x0),   
  fh3RecOutEtaPhiPt(0x0),  
  fh3GenInEtaPhiPt(0x0),   
  fh3GenOutEtaPhiPt(0x0),  
  fhnCorrelation(0x0),     
  fHistList(0x0)
{
  // Default constructor
  /*
    for(int ij  = 0;ij<kMaxJets;++ij){
      fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
      fh2PtFGen[ij] = fh2PhiFGen[ij] = fh2EtaFGen[ij] = fh2PtGenDeltaPhi[ij] =  fh2PtGenDeltaEta[ij] = 0;
      fh3RecInEtaPhiPt[ij] = fh3RecOutEtaPhiPt[ij] =fh3GenInEtaPhiPt[ij] =  fh3GenOutEtaPhiPt[ij] = 0;
      fhnCorrelation[ij] = 0;
    }
  */  
}

AliAnalysisTaskJFSystematics::AliAnalysisTaskJFSystematics(const char* name):
  AliAnalysisTaskSE(name),
  fJetHeaderRec(0x0),
  fJetHeaderGen(0x0),
  fAOD(0x0),
  fBranchRec("jets"),
  fBranchGen(""),
  fUseAODInput(kFALSE),
  fUseExternalWeightOnly(kFALSE),
  fLimitGenJetEta(kFALSE),
  fAnalysisType(0),
  fExternalWeight(1),    
  fRecEtaWindow(0.5),
  fAvgTrials(1),
  fh1Xsec(0x0),
  fh1Trials(0x0),
  fh1PtHard(0x0),
  fh1PtHardNoW(0x0),  
  fh1PtHardTrials(0x0),
  fh1NGenJets(0x0),
  fh1NRecJets(0x0),
  fh1PtRecIn(0x0),						     
  fh1PtRecOut(0x0),     
  fh1PtGenIn(0x0),      
  fh1PtGenOut(0x0),     
  fh2PtFGen(0x0),       
  fh2PhiFGen(0x0),      
  fh2EtaFGen(0x0),      
  fh2PtGenDeltaPhi(0x0),
  fh2PtGenDeltaEta(0x0),
  fh3RecInEtaPhiPt(0x0),   
  fh3RecOutEtaPhiPt(0x0),  
  fh3GenInEtaPhiPt(0x0),   
  fh3GenOutEtaPhiPt(0x0),  
  fhnCorrelation(0x0),     
  fHistList(0x0) 
{
  // Default constructor

  // Default constructor
  /*
  for(int ij  = 0;ij<kMaxJets;++ij){
    fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
    fh2PtFGen[ij] = fh2PhiFGen[ij] = fh2EtaFGen[ij] = fh2PtGenDeltaPhi[ij] =  fh2PtGenDeltaEta[ij] = 0;
    fh3RecInEtaPhiPt[ij] = fh3RecOutEtaPhiPt[ij] =fh3GenInEtaPhiPt[ij] =  fh3GenOutEtaPhiPt[ij] = 0;
    fhnCorrelation[ij] = 0;
  } 
  */ 
  DefineOutput(1, TList::Class());  
}



Bool_t AliAnalysisTaskJFSystematics::Notify()
{
//
  // Implemented Notify() to read the cross sections
  // and number of trials from pyxsec.root
  // 

  fAvgTrials = 1; // reset for each file

  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
  Float_t xsection = 0;
  Float_t ftrials  = 1;
  if(tree){
    TFile *curfile = tree->GetCurrentFile();
    if (!curfile) {
      Error("Notify","No current file");
      return kFALSE;
    }
    if(!fh1Xsec||!fh1Trials){
      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
      return kFALSE;
    }
    AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
    fh1Xsec->Fill("<#sigma>",xsection);
    // construct a poor man average trials 
    Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
    if(ftrials>=nEntries)fAvgTrials = ftrials/nEntries; 
  }
  return kTRUE;
}

void AliAnalysisTaskJFSystematics::UserCreateOutputObjects()
{

  //
  // Create the output container
  //


  if (fDebug > 1) printf("AnalysisTaskJFSystematics::UserCreateOutputObjects() \n");

  OpenFile(1);
  if(!fHistList)fHistList = new TList();

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

  //
  //  Histograms
  // 

  const Int_t nBinPt = 100;
  Double_t binLimitsPt[nBinPt+1];
  for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
    if(iPt == 0){
      binLimitsPt[iPt] = 0.0;
    }
    else {// 1.0
      binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 1.0;
    }
  }
  
  const Int_t nBinEta = 26;
  Double_t binLimitsEta[nBinEta+1] = {
    -1.6,-1.4,-1.2,-1.0,
    -0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0.0,
    0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 
    1.0, 1.2, 1.4, 1.6
  };


  const Int_t nBinPhi = 18;
  Double_t binLimitsPhi[nBinPhi+1];
  for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
    if(iPhi==0){
      binLimitsPhi[iPhi] = 0;
    }
    else{
      binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
    }
  }


  fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
  fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");

  fh1Trials = new TH1F("fh1Trials","trials event header or pyxsec file",1,0,1);
  fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");

  fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
  fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
  fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
  fh1NGenJets  = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
  fh1NRecJets = new TH1F("fh1NRecJets","N reconstructed jets",20,-0.5,19.5);

  // book the single histograms these we clone for various systematics
  // 
   fh1PtRecIn =  new TH1F("fh1PtRecIn","rec p_T input ;p_{T,rec}",nBinPt,binLimitsPt);
  fh1PtRecOut = new TH1F("fh1PtRecOut","rec p_T output jets;p_{T,rec}",nBinPt,binLimitsPt);
  fh1PtGenIn = new TH1F("fh1PtGenIn","found p_T input ;p_{T,gen}",nBinPt,binLimitsPt);
  fh1PtGenOut = new TH1F("fh1PtGenOut","found p_T output jets;p_{T,gen}",nBinPt,binLimitsPt);



  fh2PtFGen = new TH2F("fh2PtFGen","Pt Found vs. gen;p_{T,rec} (GeV/c);p_{T,gen} (GeV/c)",
			     nBinPt,binLimitsPt,nBinPt,binLimitsPt);
  fh2PhiFGen = new TH2F("fh2PhiFGen","#phi Found vs. gen;#phi_{rec};#phi_{gen}",
			nBinPhi,binLimitsPhi,nBinPhi,binLimitsPhi);
  fh2EtaFGen = new TH2F("fh2EtaFGen","#eta Found vs. gen;#eta_{rec};#eta_{gen}",
			nBinEta,binLimitsEta,nBinEta,binLimitsEta);
    
  fh2PtGenDeltaPhi = new TH2F("fh2PtGenDeltaPhi","delta phi vs. P_{T,gen};p_{T,gen} (GeV/c);#phi_{gen}-#phi_{rec}",
			      nBinPt,binLimitsPt,100,-1.0,1.0);
  fh2PtGenDeltaEta = new TH2F("fh2PtGenDeltaEta","delta eta vs. p_{T,gen};p_{T,gen} (GeV/c);#eta_{gen}-#eta_{rec}",
			      nBinPt,binLimitsPt,100,-1.0,1.0);


  fh3RecInEtaPhiPt = new TH3F("fh3RecInEtaPhiPt","Rec eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
			      nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
  fh3RecOutEtaPhiPt = new TH3F("fh3RecOutEtaPhiPt","generated found jet Rec eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
				   nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
  fh3GenInEtaPhiPt = new TH3F("fh3GenInEtaPhiPt","generated jet eta, phi, pt; #eta; #phi; p_{T,gen} (GeV/c)",
				  nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
  fh3GenOutEtaPhiPt = new TH3F("fh3GenOutEtaPhiPt","reconstructed found for Gen eta, phi, pt; #eta; #phi; p_{T,} (GeV/c)",
				   nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
  
  const int nbin[4] = {nBinPt,nBinPt,24,24};
  Double_t vLowEdge[4] = {0,0,-1.2,-1.2};
  Double_t vUpEdge[4] = {250,250,1.2,1.2};
  
  fhnCorrelation  = new THnSparseF("fhnCorrelation","Response Map", 4, nbin, vLowEdge, vUpEdge);



  fHistList->Add(fh1Xsec);
  fHistList->Add(fh1Trials);
  fHistList->Add(fh1PtHard);
  fHistList->Add(fh1PtHardNoW);
  fHistList->Add(fh1PtHardTrials);
  fHistList->Add(fh1NGenJets);
  fHistList->Add(fh1NRecJets);
  fHistList->Add(fh1PtRecIn);
  fHistList->Add(fh1PtRecOut);

  if(fBranchGen.Length()>0){
    fHistList->Add(fh1PtGenIn);
    fHistList->Add(fh1PtGenOut);
    fHistList->Add(fh2PtFGen);
    fHistList->Add(fh2PhiFGen);
    fHistList->Add(fh2EtaFGen);
    fHistList->Add(fh2PtGenDeltaEta);
    fHistList->Add(fh2PtGenDeltaPhi);
    fHistList->Add(fh3RecOutEtaPhiPt);
    fHistList->Add(fh3GenOutEtaPhiPt);      
    fHistList->Add(fh3RecInEtaPhiPt);
    fHistList->Add(fh3GenInEtaPhiPt);
    fHistList->Add(fhnCorrelation);
  }

  if(fAnalysisType==kSysJetOrder){
    // 
    for(int i = 0; i< fgkSysBins[kSysJetOrder];++i){
      TH1F *hTmp = (TH1F*)fh1PtRecIn->Clone(Form("%s_%s%d",fh1PtRecIn->GetName(),fgkSysName[kSysJetOrder],i));
      fHistList->Add(hTmp);
      hTmp = (TH1F*)fh1PtRecOut->Clone(Form("%s_%s%d",fh1PtRecOut->GetName(),fgkSysName[kSysJetOrder],i));
      fHistList->Add(hTmp);

      if(fBranchGen.Length()>0){
	hTmp = (TH1F*)fh1PtGenIn->Clone(Form("%s_%s%d",fh1PtGenIn->GetName(),fgkSysName[kSysJetOrder],i));
	fHistList->Add(hTmp);
	hTmp = (TH1F*)fh1PtGenOut->Clone(Form("%s_%s%d",fh1PtGenOut->GetName(),fgkSysName[kSysJetOrder],i));
	fHistList->Add(hTmp);
	THnSparseF *hnTmp = (THnSparseF*)fhnCorrelation->Clone(Form("%s_%s%d",fhnCorrelation->GetName(),fgkSysName[kSysJetOrder],i));
	fHistList->Add(hnTmp);
      }
    }
  }

  // =========== Switch on Sumw2 for all histos ===========
  for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
    TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
    if (h1){
      // Printf("%s ",h1->GetName());
      h1->Sumw2();
      continue;
    }
    THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
    if(hn)hn->Sumw2();
  }

  TH1::AddDirectory(oldStatus);

  }

void AliAnalysisTaskJFSystematics::Init()
{
  //
  // Initialization
  //

  Printf(">>> AnalysisTaskJFSystematics::Init() debug level %d\n",fDebug);
  if (fDebug > 1) printf("AnalysisTaskJFSystematics::Init() \n");

}

void AliAnalysisTaskJFSystematics::UserExec(Option_t */*option*/)
{
  //
  // Execute analysis for current event
  //
 
  if(fUseAODInput){    
    fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
    if(!fAOD){
      Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
      return;
    }
    // fethc the header
  }
  else{
    //  assume that the AOD is in the general output...
    fAOD  = AODEvent();
    if(!fAOD){
      Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
      return;
    }
  }
  
  if (fDebug > 1)printf("AliAnalysisTaskJFSystematics::Analysing event # %5d\n", (Int_t) fEntry);

  // ========= These pointers need to be valid in any case ======= 

  TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
  if(!aodRecJets){
    Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
    return;
  }

  // We use static arrays, not to fragment the memory
  // 
  AliAODJet genJets[kMaxJets];
  Int_t nGenJets = 0;
  AliAODJet recJets[kMaxJets];
  Int_t nRecJets = 0;

  Double_t eventW = 1;
  Double_t ptHard = 0; 
  
  Double_t nTrials = 1; // Trials for MC trigger weigth for real data
  Int_t     iProcessType = 0;
  if(fUseExternalWeightOnly){
    eventW = fExternalWeight;
  }

  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
  // this is the part where when we need to have MC information
  // we can also work on Reconstructed only when just comparing two JF
  AliMCEvent* mcEvent = MCEvent();
  if(!mcEvent){
    Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
  }
  else {
    AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
    if(pythiaGenHeader){
      nTrials = pythiaGenHeader->Trials();
      ptHard  = pythiaGenHeader->GetPtHard();
      iProcessType = pythiaGenHeader->ProcessType();
      // 11 f+f -> f+f
      // 12 f+barf -> f+barf
      // 13 f+barf -> g+g
      // 28 f+g -> f+g
      // 53 g+g -> f+barf
      // 68 g+g -> g+g
      if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
      if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
      // fetch the pythia generated jets only to be used here
      Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
      AliAODJet pythiaGenJets[kMaxJets];
      Int_t iCount = 0;
      for(int ip = 0;ip < nPythiaGenJets;++ip){
	if(iCount>=kMaxJets)continue;
	Float_t p[4];
	pythiaGenHeader->TriggerJet(ip,p);
	pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
	if(fLimitGenJetEta){
	  if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
	     pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
      }
	if(fBranchGen.Length()==0){
	  // if we have MC particles and we do not read from the aod branch
	  // use the pythia jets
	  genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
	}
      iCount++;
      }
      if(fBranchGen.Length()==0)nGenJets = iCount;    
    }
  }// if we had the MCEvent

  
  fh1Trials->Fill("#sum{ntrials}",fAvgTrials); 
  
  fh1PtHard->Fill(ptHard,eventW);
  fh1PtHardNoW->Fill(ptHard,1);
  fh1PtHardTrials->Fill(ptHard,nTrials);

  // If we set a second branch for the input jets fetch this 
  if(fBranchGen.Length()>0){
    TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
    if(aodGenJets){
      Int_t iCount = 0;
      for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
	if(iCount>=kMaxJets)continue;
	AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
	if(!tmp)continue;
	if(fLimitGenJetEta){
	  if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
	     tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
	}
	genJets[iCount] = *tmp;
	iCount++;
      }
      nGenJets = iCount;
    }
    else{
      Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
    }
  }

  fh1NGenJets->Fill(nGenJets);
  // We do not want to exceed the maximum jet number
  nGenJets = TMath::Min(nGenJets,kMaxJets);


  //
  // Fetch the reconstructed jets...
  //


  nRecJets = TMath::Min(nRecJets,kMaxJets);

  for(int ir = 0;ir < nRecJets;++ir){
    AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
    if(!tmp)continue;
    recJets[ir] = *tmp;
  }


  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);

  //
  // Relate the jets
  //
  Int_t iGenIndex[kMaxJets];    // Index of the generated jet for i-th rec -1 if none
  Int_t iRecIndex[kMaxJets];    // Index of the rec jet for i-th gen -1 if none
  
  for(int i = 0;i<kMaxJets;++i){
    iGenIndex[i] = iRecIndex[i] = -1;
  }


  AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
					    iGenIndex,iRecIndex,fDebug);
  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);

  if(fDebug){
    for(int i = 0;i<kMaxJets;++i){
      if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]); 
      if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]); 
    }
  }

  //
  // Now the premliminaries are over, lets do the jet analysis
  //


  Double_t value[4]; // for the thnsparse 
  // loop over reconstructed jets
  for(int ir = 0;ir < nRecJets;++ir){
    Double_t ptRec = recJets[ir].Pt();
    Double_t phiRec = recJets[ir].Phi();
    if(phiRec<0)phiRec+=TMath::Pi()*2.;    
    Double_t etaRec = recJets[ir].Eta();
    if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
    fh1PtRecIn->Fill(ptRec,eventW);
    if(fAnalysisType==kSysJetOrder)((TH1F*)fHistList->FindObject(Form("fh1PtRecIn_%s%d",fgkSysName[kSysJetOrder],ir)))->Fill(ptRec,eventW);
    fh3RecInEtaPhiPt->Fill(etaRec,phiRec,ptRec,eventW);
    // Fill Correlation
    Int_t ig = iGenIndex[ir];
    if(ig>=0 && ig<nGenJets){
      if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
      fh1PtRecOut->Fill(ptRec,eventW);
      if(fAnalysisType==kSysJetOrder)((TH1F*)fHistList->FindObject(Form("fh1PtRecOut_%s%d",fgkSysName[kSysJetOrder],ir)))->Fill(ptRec,eventW);
      Double_t ptGen  = genJets[ig].Pt();
      Double_t phiGen = genJets[ig].Phi();
      if(phiGen<0)phiGen+=TMath::Pi()*2.; 
      Double_t etaGen = genJets[ig].Eta();

      fh3RecOutEtaPhiPt->Fill(etaRec,phiRec,ptRec,eventW);

      value[0] = ptRec;
      value[1] = ptGen;
      value[2] = etaRec;
      value[3] = etaGen;
      
      fhnCorrelation->Fill(value,eventW);
      if(fAnalysisType==kSysJetOrder)((THnSparseF*)fHistList->FindObject(Form("fhnCorrelation_%s%d",fgkSysName[kSysJetOrder],ir)))->Fill(value,eventW);
      // 
      // we accept only jets which are detected within a smaller window, to 
      // avoid ambigious pair association at the edges of the acceptance
      // 

      if(TMath::Abs(etaRec)<fRecEtaWindow){
	fh2PtFGen->Fill(ptRec,ptGen,eventW);
	fh2PhiFGen->Fill(phiRec,phiGen,eventW);
	fh2EtaFGen->Fill(etaRec,etaGen,eventW);
	fh2PtGenDeltaEta->Fill(ptGen,etaGen-etaRec,eventW);
	fh2PtGenDeltaPhi->Fill(ptGen,phiGen-phiRec,eventW);
      }// if etarec in window
    }//if ig valid
  }// loop over reconstructed jets
  

  // Now llop over generated jets

  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
  for(int ig = 0;ig < nGenJets;++ig){
    Double_t ptGen = genJets[ig].Pt();
    // Fill Correlation
    Double_t phiGen = genJets[ig].Phi();
    if(phiGen<0)phiGen+=TMath::Pi()*2.;    
    Double_t etaGen = genJets[ig].Eta();
    fh3GenInEtaPhiPt->Fill(etaGen,phiGen,ptGen,eventW);
    fh1PtGenIn->Fill(ptGen,eventW);
    if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
    if(fAnalysisType==kSysJetOrder)((TH1F*)fHistList->FindObject(Form("fh1PtGenIn_%s%d",fgkSysName[kSysJetOrder],ig)))->Fill(ptGen,eventW);
    Int_t ir = iRecIndex[ig];
    if(ir>=0&&ir<nRecJets){
      fh1PtGenOut->Fill(ptGen,eventW);
      fh3GenOutEtaPhiPt->Fill(etaGen,phiGen,ptGen,eventW);
      if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
      if(fAnalysisType==kSysJetOrder)((TH1F*)fHistList->FindObject(Form("fh1PtGenOut_%s%d",fgkSysName[kSysJetOrder],ig)))->Fill(ptGen,eventW);
    }
  }// loop over reconstructed jets

  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
  PostData(1, fHistList);
}

void AliAnalysisTaskJFSystematics::Terminate(Option_t */*option*/)
{
// Terminate analysis
//
    if (fDebug > 1) printf("AnalysisTaskJFSystematics: Terminate() \n");
    // Plot 


}