/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *																          *
 * Authors: Svein Lindal                                                  *
 * Version 1.0			                                                  *
 *																          *
 *																          *
 * 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.				  *
 **************************************************************************/

////////////////////////////////////////////////
//--------------------------------------------- 
// Class doing conversion gamma dPhi correlations
// Gamma Conversion analysis
//---------------------------------------------
////////////////////////////////////////////////

#include "AliAnalysisTaskdPhi.h"

#include <TH2I.h>
#include <TList.h>
#include <TChain.h>
#include <TFile.h>
#include <AliAnalysisManager.h>
#include <AliInputEventHandler.h>
#include <AliESDInputHandler.h>
#include <AliAODInputHandler.h>
#include <TGeoGlobalMagField.h>

#include "AliConversionTrackCuts.h"
#include "AliConversionCuts.h"
#include "AliConversionMesonCuts.h"
#include "AliAODConversionPhoton.h"
#include "AliAODConversionMother.h"
#include "TGrid.h"
// #include "AliAnaConvCorrPhoton.h"
// #include "AliAnaConvCorrPion.h"
// #include "AliAnaConvIsolation.h"
#include "AliV0ReaderV1.h"
// Author Svein Lindal <slindal@fys.uio.no>
using namespace std;

ClassImp(AliAnalysisTaskdPhi)


//________________________________________________________________________
AliAnalysisTaskdPhi::AliAnalysisTaskdPhi(const char *name) : AliAnalysisTaskSE(name),
  fHistograms(NULL),
  fCorrSparse(NULL),
  fTrigSparse(NULL),
  fTrackSparse(NULL),
  fMassSparse(NULL),
  fV0Reader(NULL),
  fSaveReaderHists(kFALSE),
  fV0Filter(NULL),
  fV0Filters(),
  fPhotonFilter(NULL),
  fMesonFilter(NULL),
  fMesonFilters(),
  fTrackFilter(NULL),
  fTrackFilters(),
  fGammas(),
  fTracks(),
  hMEvents(NULL),
  hTrackCent(NULL),
  hTrigPt(NULL),
  hTrackPt(NULL), 
  hTrigPhi(NULL),
  fDeltaAODBranchName("AliAODGammaConversion_gamma"), 
  fAxistPt(),
  fAxiscPt(),
  fAxisdEta(),
  fAxisTrigEta(),
  fAxisAssEta(),
  fAxisdPhi(),
  fAxisCent(),
  fAxisZ(), 
  fAxisPiM(), 
  fAxisTrackFilters(),							     
  fAxisV0Filters(),							     
  fAxisMesonFilters(),							     
  fkTrackAxis(kFALSE),
  fkV0Axis(kFALSE),
  fkPionAxis(kFALSE),
  fAxesList(), 
  fTrigAxesList(), 
  fTrackAxesList(), 
  fMassAxesList(),
  fDoPhoton(kFALSE), 
  fCorrectionMap(NULL)
 {
   //constructor
   SetUpBins();

   DefineInput(0, TChain::Class());
   DefineOutput(1, TList::Class());

   fGammas.SetOwner(kTRUE);
   fTracks.SetOwner(kTRUE);



   
 }



 //________________________________________________________________________
 AliAnalysisTaskdPhi::~AliAnalysisTaskdPhi(){
   //destructor

   if(fV0Filter)
	 delete fV0Filter;
   fV0Filter = NULL;

   if(fMesonFilter)
	 delete fMesonFilter;
   fMesonFilter = NULL;

   if(fPhotonFilter)
	 delete fPhotonFilter;
   fPhotonFilter = NULL;

   if(fHistograms)
	 delete fHistograms;
   fHistograms = NULL;

   if(fTrackFilter)
     delete fTrackFilter;
   fTrackFilter = NULL;

   fGammas.Delete();
   fTracks.Delete();

 }

 ///________________________________________________________________________
void AliAnalysisTaskdPhi::SetUpBins() {

   fAxisTrigEta.SetNameTitle("tEta", "Eta");
   fAxisTrigEta.Set(320, -0.8, 0.8);

   fAxisAssEta.SetNameTitle("aEta", "Eta");
   fAxisAssEta.Set(360, -0.9, 0.9);

   fAxisdEta.SetNameTitle("dEta", "Eta");
   fAxisdEta.Set(34, -1.7, 1.7);

   fAxisdPhi.SetNameTitle("dPhi", "delta Phi");
   fAxisdPhi.Set(32, -TMath::PiOver2(), 3*TMath::PiOver2());

   Double_t tptbins[10] = {0.1, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0, 15, 50, 100};                                    fAxistPt.SetNameTitle("tPt", "tPt");
   fAxistPt.Set(9, tptbins);

   Double_t cptbins[13] = {0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 8.0, 10.0, 25, 50, 100};                     fAxiscPt.SetNameTitle("cPt", "cPt");
   fAxiscPt.Set(12, cptbins);

   fAxisZ.SetNameTitle("vertexz", "Z");
   fAxisZ.Set(4, -10, 10);

   Double_t centbins[5] = {0, 10, 30, 60, 100.1};
   fAxisCent.SetNameTitle("centrality", "Cent");
   fAxisCent.Set(4, centbins);

   fAxisPiM.SetNameTitle("InvMassPi0", "Invariant mass");
   Double_t mbins[7] = {0.1, 0.11, 0.12, 0.15, 0.16, 0.18, 0.2};
   fAxisPiM.Set(6, mbins);

}


 ///________________________________________________________________________
// void AliAnalysisTaskdPhi::SetUpCorrObjects() {
   // //Set up corr objects
   // AliDebug(AliLog::kDebug + 5, "Set Up corr objects");

   // if(fDoPhoton) {
   //   fPhotonCorr = new AliAnaConvCorrPhoton("PhotonCorr","photon %s");
   //   SetUpCorrAxes(fPhotonCorr);
   //   fPhotonCorr->CreateHistograms();
   //   fHistograms->Add(fPhotonCorr->GetHistograms());
   // }

   // fPionCorr = new AliAnaConvCorrPion("PionCorr", "pion");
   // SetUpCorrAxes(fPionCorr);
   // fPionCorr->GetAxisM().Set(fAxisPiM.GetNbins(), fAxisPiM.GetXbins()->GetArray());
   // fPionCorr->CreateHistograms();
   // fHistograms->Add(fPionCorr->GetHistograms());
   //}

 ///________________________________________________________________________    
// void AliAnalysisTaskdPhi::SetUpCorrAxes(AliAnaConvCorrBase * corr) {
   ///Set up axes in corr object
   // corr->GetAxisCent().Set(fAxisCent.GetNbins(), fAxisCent.GetXbins()->GetArray());
   // const Double_t * zbins = fAxisZ.GetXbins()->GetArray();
   // if(zbins) {
   //   corr->GetAxisZ().Set(fAxisZ.GetNbins(), fAxisZ.GetXbins()->GetArray());
   // } else {
   //   corr->GetAxisZ().Set(fAxisZ.GetNbins(), fAxisZ.GetBinLowEdge(1), fAxisZ.GetBinUpEdge(fAxisZ.GetNbins()));
   // }

   // corr->GetAxistPt().Set(fAxistPt.GetNbins(), fAxistPt.GetXbins()->GetArray());
   // corr->GetAxiscPt().Set(fAxiscPt.GetNbins(), fAxiscPt.GetXbins()->GetArray());
   //   corr->GetAxisdEta().Set(fAxisdEta.GetNbins(), fAxisdEta.GetBinLowEdge(1), fAxisdEta.GetBinUpEdge(fAxisdEta.GetNbins()));
   // corr->GetAxisTrigEta().Set(fAxisTrigEta.GetNbins(), fAxisTrigEta.GetBinLowEdge(1), fAxisTrigEta.GetBinUpEdge(fAxisTrigEta.GetNbins()));
   // corr->GetAxisAssEta().Set(fAxisAssEta.GetNbins(), fAxisAssEta.GetBinLowEdge(1), fAxisAssEta.GetBinUpEdge(fAxisAssEta.GetNbins()));
// }    


//________________________________________________________________________
void AliAnalysisTaskdPhi::UserCreateOutputObjects() {
  // Create histograms
  // TGrid::Connect("alien://",0,0,"t");
  // if(!gGrid) AliWarning("no GGrid");
  // TFile *tfile = TFile::Open("alien:///alice/cern.ch/user/s/slindal/trackMap.root", "READ");
  // if(tfile) {
  //   THnF * corrmap = dynamic_cast<THnF*>(tfile->Get("hTrackCorr"));
  //   if (corrmap) {
  //     fCorrectionMap = dynamic_cast<THnF*>(THn::CreateHn("corr", "corr", corrmap));
  //     for(Int_t i = 0; i < fCorrectionMap->GetNdimensions(); i++) {
  // 	TAxis * axis = fCorrectionMap->GetAxis(i);
  // 	axis->SetRange(1, axis->GetNbins());
  //     }
    
  //     cout << "yessssssssssssssssssssssssssssssssssssssssssssssssss"<<endl;
  //   } else {
  //     cout << "xxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxx xxxxxxxxxxxxxxxx"<<endl;
  //   }
  //   tfile->Close();
  // } else {
  //   cout << "no tfile shit shit shit "<<endl;
  //   AliFatal("file not ther!!!");
  // }

  
  fHistograms = new TList();
  fHistograms->SetName("dPhi_histograms");
  fHistograms->SetOwner(kTRUE);
  
  
  if(!fV0Reader){
    fV0Reader=(AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1");
  }
  
  if(!fV0Reader){
    printf("Error: No V0 Reader");
  } // GetV0Reader
  
  
  if(fSaveReaderHists) {
    AliConversionCuts * v0cuts = fV0Reader->GetConversionCuts();
    if(v0cuts) {
      TList * histograms = v0cuts->GetCutHistograms();
      if(!histograms) {
	AliWarning("initializing v0 reader hists");
	v0cuts->InitCutHistograms("V0Reader", kTRUE);
      }
      histograms = v0cuts->GetCutHistograms();
      if(histograms) {
	fHistograms->Add(histograms);
      }
    }
  }
  
  for(Int_t igf = 0; igf < fV0Filters[0].GetEntriesFast(); igf ++){
    AliConversionCuts * f = dynamic_cast<AliConversionCuts*>(fV0Filters[0].At(igf));
    if(f) {
      TList * histograms = f->GetCutHistograms();
      if(histograms) fHistograms->Add(f->GetCutHistograms());
    }
  }

  for(Int_t igf = 0; igf < fV0Filters[1].GetEntriesFast(); igf ++){
    AliConversionCuts * f = dynamic_cast<AliConversionCuts*>(fV0Filters[1].At(igf));
    if(f) {
      TList * histograms = f->GetCutHistograms();
      if(histograms) fHistograms->Add(f->GetCutHistograms());
    }
  }

  for(Int_t igf = 0; igf < fMesonFilters[0].GetEntriesFast(); igf ++){
    AliConversionMesonCuts * f = dynamic_cast<AliConversionMesonCuts*>(fMesonFilters[0].At(igf));
    if(f) {
      TList * histograms = f->GetCutHistograms();
      if(histograms) fHistograms->Add(f->GetCutHistograms());
    }
  }

  for(Int_t igf = 0; igf < fMesonFilters[1].GetEntriesFast(); igf ++){
    AliConversionMesonCuts * f = dynamic_cast<AliConversionMesonCuts*>(fMesonFilters[1].At(igf));
    if(f) {
      TList * histograms = f->GetCutHistograms();
      if(histograms) fHistograms->Add(f->GetCutHistograms());
    }
  }

  if(fV0Filter) {
    fV0Filter->InitCutHistograms("V0Filter", kFALSE);
    fHistograms->Add(fV0Filter->GetCutHistograms());
  }
  if(fMesonFilter) {
    fMesonFilter->InitCutHistograms("PionFilter", kFALSE);
    fHistograms->Add(fMesonFilter->GetCutHistograms());
  }
  if(fPhotonFilter) {
    fPhotonFilter->InitCutHistograms("PhotonFilter", kFALSE);
    fHistograms->Add(fPhotonFilter->GetCutHistograms());
  }


  AliConversionTrackCuts * tc = dynamic_cast<AliConversionTrackCuts*>(fTrackFilter);
  if(tc) {
    fHistograms->Add(tc->CreateHistograms());
  }


   // for(Int_t i = 0; i < fTrackFilters.GetEntriesFast(); i++) {
   //   AliConversionTrackCuts * tc = dynamic_cast<AliConversionTrackCuts*>(fTrackFilters.At(i));
   //   if(tc) fHistograms->Add(tc->CreateHistograms());
   // }

  //SetUpCorrObjects();

   ///Set up ME histograms
   TList * MEHistograms = new TList();
   MEHistograms->SetName("MEHistograms");
   MEHistograms->SetOwner(kTRUE);
   fHistograms->Add(MEHistograms);
   
   hMEvents = new TH2I("hMEvents", "Nevents vs centrality vertexz",
		       fAxisZ.GetNbins(), fAxisZ.GetXbins()->GetArray(),
		       fAxisCent.GetNbins(), fAxisCent.GetXbins()->GetArray());
   MEHistograms->Add(hMEvents);

   hTrackCent = new TH2I("hTrackCent", "N accepted tracks vs centrality",
			 fAxisCent.GetNbins() > 1 ? (int) (10*(fAxisCent.GetXmax() - fAxisCent.GetXmin()))  : 1, 
			 fAxisCent.GetXmin(), fAxisCent.GetXmax(),
			 fAxisCent.GetNbins() > 1 ? 900 : 50, 
			 0,
			 fAxisCent.GetNbins() > 1 ? 1800 : 50);
   MEHistograms->Add(hTrackCent);

   hTrigPt = new TH3F("hTrigPt", "trigger pt", 100, 0., 10., 
		      10, 0., 50., 
		      5,  0.05, 0.2);
   MEHistograms->Add(hTrigPt);
   hTrackPt = new TH2F("hTrackPt", "track pt", 100, 0, 10, 10, 0, 50);//fAxisCent.GetNbins(), fAxisCent.GetXbins()->GetArray()); 
   MEHistograms->Add(hTrackPt);
   hTrigPhi = new TH1F("hTrigPhi", "trigger pt", 32, 0, 2*TMath::Pi());
   MEHistograms->Add(hTrigPhi);



   Int_t ntrackfilters[2] = {fTrackFilters[0].GetEntriesFast(), fTrackFilters[1].GetEntriesFast()};
   fkTrackAxis = kTRUE;
   fAxisTrackFilters.SetNameTitle("trackCuts", "trackCuts");
   fAxisTrackFilters.Set(ntrackfilters[0] + ntrackfilters[1] + 1, -ntrackfilters[0] -0.5, ntrackfilters[1]  + 0.5);

   Int_t nV0filters[2] = {fV0Filters[0].GetEntriesFast(), fV0Filters[1].GetEntriesFast()};
   fkV0Axis = kTRUE;
   fAxisV0Filters.SetNameTitle("V0Cuts", "V0Cuts");
   fAxisV0Filters.Set(nV0filters[0] + nV0filters[1] + 1, -nV0filters[0] -0.5, nV0filters[1]  + 0.5);
   
   Int_t nmesonfilters[2] = {fMesonFilters[0].GetEntriesFast(), fMesonFilters[1].GetEntriesFast()};
   fkPionAxis = kTRUE;
   fAxisMesonFilters.SetNameTitle("mesonCuts", "mesonCuts");
   fAxisMesonFilters.Set(nmesonfilters[0] + nmesonfilters[1] + 1, -nmesonfilters[0] -0.5, nmesonfilters[1]  + 0.5);
 
  fAxesList.AddAt(&fAxisdEta, 0);
  fAxesList.AddAt(&fAxisdPhi, 1);
  fAxesList.AddAt(&fAxistPt, 2);
  fAxesList.AddAt(&fAxiscPt, 3);
  fAxesList.AddAt(&fAxisCent, 4);
  fAxesList.AddAt(&fAxisZ, 5);
  fAxesList.AddAt(&fAxisPiM, 6);
  fAxesList.AddAt(&fAxisTrackFilters, 7);
  fAxesList.AddAt(&fAxisV0Filters, 8);
  fAxesList.AddAt(&fAxisMesonFilters, 9);
  
  fTrackAxesList.AddAt(&fAxisAssEta, 0);
  fTrackAxesList.AddAt(&fAxistPt, 1);
  fTrackAxesList.AddAt(&fAxiscPt, 2);
  fTrackAxesList.AddAt(&fAxisCent, 3);
  fTrackAxesList.AddAt(&fAxisZ, 4);
  //fTrackAxesList.AddAt(&fAxisPiM, 5);
  fTrackAxesList.AddAt(&fAxisTrackFilters, 5);
  fTrackAxesList.AddAt(&fAxisV0Filters, 6);
  fTrackAxesList.AddAt(&fAxisMesonFilters, 7);
   
  fTrigAxesList.AddAt(&fAxisTrigEta, 0);
  fTrigAxesList.AddAt(&fAxistPt, 1);
  fTrigAxesList.AddAt(&fAxisCent, 2);
  fTrigAxesList.AddAt(&fAxisZ, 3);
  fTrigAxesList.AddAt(&fAxisPiM, 4);
  fTrigAxesList.AddAt(&fAxisV0Filters, 5);
  fTrigAxesList.AddAt(&fAxisMesonFilters, 6);
   
  TList masslist;
  TAxis massax;
  massax.SetNameTitle("mass", "mass");
  massax.Set(360, 0.04, 0.4); //hardcoded! change also in filling!
  
  masslist.AddAt(&massax, 0);
  masslist.AddAt(&fAxistPt, 1);
  masslist.AddAt(&fAxisCent, 2);
  masslist.AddAt(&fAxisV0Filters, 3);
  masslist.AddAt(&fAxisMesonFilters, 4);
  
  fCorrSparse = CreateSparse(TString("pionSparse"), TString("pionSparse"), &fAxesList);
  fTrackSparse = CreateSparse(TString("trackSparse"), TString("trackSparse"), &fTrackAxesList);
  fTrigSparse = CreateSparse(TString("trigSparse"), TString("trigSparse"), &fTrigAxesList);
  fMassSparse = CreateSparse("massSparse", "massSparse", &masslist);
  
  fHistograms->Add(fCorrSparse);
  fHistograms->Add(fTrackSparse);
  fHistograms->Add(fTrigSparse);
  fHistograms->Add(fMassSparse);


  ///Add gamma and track containers:
  for(Int_t i = 0; i < fV0Filters[1].GetEntriesFast() + 1; i++) {
    fGammas.Add(new TObjArray());
  }

  for(Int_t i = 0; i < fTrackFilters[1].GetEntriesFast() + 1; i++) {
    fTracks.Add(new TObjArray());
  }



  PostData(1, fHistograms);
}

///________________________________________________________________________
THnSparseF * AliAnalysisTaskdPhi::CreateSparse(TString nameString, TString titleString, TList * axesList) {
  //Create sparse
  const Int_t dim = axesList->GetSize();
  
  TAxis * axes[dim];
  Int_t   bins[dim];
  Double_t min[dim];
  Double_t max[dim];

  for(Int_t i = 0; i<dim; i++) {
    TAxis * axis = dynamic_cast<TAxis*>(axesList->At(i));
    if(axis) axes[i] = axis;
    else {
      cout << "AliAnalysisTaskdPhi::CreateSparse: Error error, all the axes are not present in axis list" << endl;
      return NULL;
    }
  }
  
  for(Int_t i = 0; i<dim; i++) {
    //cout << axes[i]->GetTitle() << endl;
    bins[i] = axes[i]->GetNbins(); 
    min[i] = axes[i]->GetBinLowEdge(1);
    max[i] = axes[i]->GetBinUpEdge(axes[i]->GetNbins());
  }

  THnSparseF * sparse = new THnSparseF(Form("%s", nameString.Data()), 
				       Form("%s", titleString.Data()), 
				       dim, bins, min, max);
  
  for(Int_t i = 0; i<dim; i++) {
    sparse->GetAxis(i)->SetNameTitle(axes[i]->GetName(), axes[i]->GetTitle() );
    if(axes[i]->GetXbins()->GetSize() > 0) {
      sparse->SetBinEdges(i, axes[i]->GetXbins()->GetArray() );
    }
  }
  return sparse;
}

//________________________________________________________________________
void AliAnalysisTaskdPhi::UserExec(Option_t *) {
  ///User exec. 
  ///This is a very ugly function, cut the complexity of the logic demands it.
  

  AliESDEvent * esdEvent = dynamic_cast<AliESDEvent*>(fInputEvent);
  if(esdEvent) {
    if (!TGeoGlobalMagField::Instance()->GetField()) esdEvent->InitMagneticField(); 
  }


  //if(! fV0Filter->EventIsSelected(fInputEvent)) return;
  if(!fV0Reader){
    AliError("Error: No V0 Reader");
    return;
  } // GetV0Reader
  
  if(!fV0Reader->IsEventSelected()) {
    return;
  }

  AliDebug(5, "Processing event");
  
 
  for(Int_t i = 0; i < fGammas.GetEntriesFast(); i++) {
    static_cast<TObjArray*>(fGammas.At(i))->Clear();
  }

  for(Int_t i = 0; i < fTracks.GetEntriesFast(); i++) {
    static_cast<TObjArray*>(fTracks.At(i))->Clear();
  }


  
  AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
  Bool_t isAOD=man->GetInputEventHandler()->IsA()==AliAODInputHandler::Class();
  
  AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
  if (!inputHandler) {
    cout << "cout no input event handler"<<endl;
    return;
  }

  


  for(Int_t igf = 0; igf < fV0Filters[0].GetEntriesFast(); igf++) {
    AliConversionCuts * f = dynamic_cast<AliConversionCuts*>(fV0Filters[0].At(igf));
    if ( f && !f->GetPIDResponse() ) {
      if ( inputHandler->GetPIDResponse() ){
	f->SetPIDResponse( inputHandler->GetPIDResponse() );
      } else {
	if (isAOD){
	  if (!f->GetPIDResponse()){
	    f->InitAODpidUtil(1);
	  }
	}
      }
    } else {
      break;
    }
  }

  for(Int_t igf = 0; igf < fV0Filters[1].GetEntriesFast(); igf++) {
    AliConversionCuts * f = dynamic_cast<AliConversionCuts*>(fV0Filters[1].At(igf));
    if ( f && !f->GetPIDResponse() ) {
      if ( inputHandler->GetPIDResponse() ){
	f->SetPIDResponse( inputHandler->GetPIDResponse() );
      } else {
	if (isAOD){
	  if (!f->GetPIDResponse()){
	    f->InitAODpidUtil(1);
	  }
	}
      }
    } else {
      break;
    }
  }

  
  if ( fV0Filter && !fV0Filter->GetPIDResponse() ) {
    if ( inputHandler->GetPIDResponse() ){
      fV0Filter->SetPIDResponse( inputHandler->GetPIDResponse() );
    } else {
      
      //AOD case
      if (isAOD){
	if (!fV0Filter->GetPIDResponse()){
	  fV0Filter->InitAODpidUtil(1);
	}
      }
    }
  }


  ///Initialize track cuts. Delete tracks that have been constrained to vertex (copies)
  AliConversionTrackCuts * tc = dynamic_cast<AliConversionTrackCuts*>(fTrackFilter);
  if(tc) {
    tc->SetEvent(fInputEvent);
    tc->DeleteTracks();
  }
  
  for(Int_t i = 0; i < fTrackFilters[0].GetEntriesFast(); i++){
    AliConversionTrackCuts * tct = dynamic_cast<AliConversionTrackCuts*>(fTrackFilters[0].At(i));
    if(tct) {
      tct->SetEvent(fInputEvent);
      tct->DeleteTracks();
    }
  }
  for(Int_t i = 0; i < fTrackFilters[1].GetEntriesFast(); i++){
    AliConversionTrackCuts * tct = dynamic_cast<AliConversionTrackCuts*>(fTrackFilters[1].At(i));
    if(tct) {
      tct->SetEvent(fInputEvent);
      tct->DeleteTracks();
    }
  }
  

  Double_t centrality = 0.0;
  Double_t vertexz = fInputEvent->GetPrimaryVertex()->GetZ();
  if(isAOD) {
    AliAODHeader * header = static_cast<AliAODHeader*>(fInputEvent->GetHeader());
    centrality = header->GetCentrality();
  } else {
    centrality = static_cast<AliESDEvent*>(fInputEvent)->GetCentrality()->GetCentralityPercentile("V0M");
  }
  
  
  const Int_t centBin = GetBin(fAxisCent, centrality);
  const Int_t vertexBin = GetBin(fAxisZ, vertexz);
  
  
  if(DebugLevel () > 4) {
    cout << "centrality: " << centrality <<  " " << GetBin(fAxisCent, centrality) << endl;
    cout << "vertexz: " << vertexz <<  " " << GetBin(fAxisZ, vertexz) << endl;
  }
  
  
  
  if(centBin < 0 || vertexBin < 0) {
    //	AliError("bin out of range");
    return;
  }
  
  
  
  //TClonesArray * aodGammas = GetConversionGammas(isAOD);
  TClonesArray * aodGammas = fV0Reader->GetReconstructedGammas();
  if(!aodGammas) {
    AliError("no aod gammas found!");
    return;
  }
  
  TObjArray * ggammas = static_cast<TObjArray*>(fGammas.At(0));

  ///Fill arrays of accepted gammas
  if(DebugLevel() > 1) printf("Number of conversion gammas %d \n", aodGammas->GetEntriesFast());
  for(Int_t ig = 0; ig < aodGammas->GetEntriesFast(); ig++) {
    AliAODConversionPhoton * photon = dynamic_cast<AliAODConversionPhoton*>(aodGammas->At(ig));
    if(!photon) continue;
    if(!fV0Filter || fV0Filter->PhotonIsSelected(photon, fInputEvent)) {
      ggammas->Add(photon);
    } else {
      for(Int_t igf = 0; igf < fV0Filters[1].GetEntriesFast(); igf++) {
	AliConversionCuts * gfilter = dynamic_cast<AliConversionCuts*>(fV0Filters[1].At(igf));
	if(gfilter && gfilter->PhotonIsSelected(photon, fInputEvent)) {
	  static_cast<TObjArray*>(fGammas.At(igf+1))->Add(photon);
	}
      }
    }
  }

  Bool_t lowgmap[fV0Filters[0].GetEntriesFast()][ggammas->GetEntriesFast()];

  for(Int_t ig = 0; ig < ggammas->GetEntriesFast(); ig++ ) {
    AliAODConversionPhoton * gamma = static_cast<AliAODConversionPhoton*>(ggammas->At(ig));
    for(Int_t igf = 0; igf < fV0Filters[0].GetEntriesFast(); igf++) {
      
      AliConversionCuts * v0cuts = static_cast<AliConversionCuts*>(fV0Filters[0].At(igf));
      if(!(v0cuts->PhotonIsSelected(gamma, fInputEvent))) {
	lowgmap[igf][ig] = kTRUE;
      } else {
	lowgmap[igf][ig] = kFALSE;
      }
    }
  }
  

  if(DebugLevel() > 4) printf("Number of accepted gammas %d \n", ggammas->GetEntriesFast());
  hMEvents->Fill(vertexz, centrality);
  
  ///create track array
  const Int_t ntrackfilters[2] = { fTrackFilters[0].GetEntriesFast(), fTrackFilters[1].GetEntriesFast()};
  
  TObjArray * ttracks = static_cast<TObjArray*>(fTracks.At(0));
  const Double_t aetalim[2] = { fAxisAssEta.GetXmin(), fAxisAssEta.GetXmax()};
  const Double_t aptlim[2] = { fAxiscPt.GetXmin(), fAxiscPt.GetXmax()};
  for(Int_t iTrack = 0; iTrack < fInputEvent->GetNumberOfTracks(); iTrack++) {
    AliVTrack * track = static_cast<AliVTrack*>(fInputEvent->GetTrack(iTrack));
    if(track->Pt() < aptlim[0] || track->Pt() > aptlim[1]) continue;
    if(track->Eta() < aetalim[0] || track->Eta() > aetalim[1]) continue;
    if(fTrackFilter->IsSelected(track)) {
      hTrackPt->Fill(track->Pt(), centrality);
      ttracks->Add(track);
    } else {
      ///upside cuts
      for(Int_t itf = 1; itf < ntrackfilters[1] + 1; itf++) {
	AliAnalysisCuts * trackCuts = static_cast<AliAnalysisCuts*>(fTrackFilters[1].At(itf -1));
	if(trackCuts->IsSelected(track)) {
	  static_cast<TObjArray*>(fTracks.At(itf))->Add(track);
	}
      }
    }
  }


  Bool_t lowtrackmap[ntrackfilters[0]][ttracks->GetEntriesFast()];
  ///Check lowside cuts
  for(Int_t iTrack = 0; iTrack < ttracks->GetEntriesFast(); iTrack++ ) {
    AliVTrack * track = static_cast<AliVTrack*>(ttracks->At(iTrack));
    for(Int_t itf = 0; itf < ntrackfilters[0]; itf++) {
      AliAnalysisCuts * trackCuts = static_cast<AliAnalysisCuts*>(fTrackFilters[0].At(itf));
      if(!trackCuts->IsSelected(track)) {
	lowtrackmap[itf][iTrack] = kTRUE;
      } else {
	lowtrackmap[itf][iTrack] = kFALSE;
      }
    }
  }

  hTrackCent->Fill(centrality, ttracks->GetEntriesFast());
  
  const Double_t etalim[2] = { fAxisTrigEta.GetXmin(), fAxisTrigEta.GetXmax()};
  if(DebugLevel() > 4) printf("Number of accepted gammas, tracks %d  %d \n", ggammas->GetEntriesFast(), ttracks->GetEntriesFast());
 
  //AliAnaConvCorrBase * gCorr = fPhotonCorr; //GetCorrObject(vertexBin, centBin, fPhotonCorr);
  //  AliAnaConvCorrPion * piCorr = fPionCorr; //static_cast<AliAnaConvCorrPion*>(GetCorrObject(vertexBin, centBin, fPionCorr));
  // if(!piCorr) {
  //   AliError("corr object missing");
  //   return;
  // }
  
  TObjArray pions;

  ///corr values
  Double_t dphivalues[fAxesList.GetSize()];
  dphivalues[4] = centrality;
  dphivalues[5] = vertexz;

  ///Trigger me counters
  Double_t trigValues[7];
  trigValues[2] = centrality;
  trigValues[3] = vertexz;

  ///Mass histogram
  Double_t massval[5];
  massval[2] = centrality;
  
  ///Set up track me counters and initialize
  const Int_t nbins = fAxistPt.GetNbins();
  Bool_t tmap[fAxistPt.GetNbins()];
  Bool_t lv0tmap[fAxistPt.GetNbins()][fV0Filters[0].GetEntriesFast()];
  Bool_t uv0tmap[fAxistPt.GetNbins()][fV0Filters[1].GetEntriesFast()];

  Bool_t lpitmap[fAxistPt.GetNbins()][fMesonFilters[0].GetEntriesFast()];
  Bool_t upitmap[fAxistPt.GetNbins()][fMesonFilters[1].GetEntriesFast()];
  
  for(Int_t igf = 0; igf < fV0Filters[0].GetEntriesFast(); igf++) {
    for(Int_t ptbin = 0; ptbin < fAxistPt.GetNbins(); ptbin++) {
      lv0tmap[ptbin][igf] = kFALSE;
    }
  }

  for(Int_t igf = 0; igf < fV0Filters[1].GetEntriesFast(); igf++) {
    for(Int_t ptbin = 0; ptbin < fAxistPt.GetNbins(); ptbin++) {
      uv0tmap[ptbin][igf] = kFALSE;
    }
  }

  for(Int_t igf = 0; igf < fMesonFilters[0].GetEntriesFast(); igf++) {
    for(Int_t ptbin = 0; ptbin < fAxistPt.GetNbins(); ptbin++) {
      lpitmap[ptbin][igf] = kFALSE;
    }
  }

  for(Int_t igf = 0; igf < fMesonFilters[1].GetEntriesFast(); igf++) {
    for(Int_t ptbin = 0; ptbin < fAxistPt.GetNbins(); ptbin++) {
      upitmap[ptbin][igf] = kFALSE;
    }
  }

  for(Int_t ptbin = 0; ptbin < nbins; ptbin++){
    tmap[ptbin] = kFALSE;
  }

  /////////////////////////////////////////


  for(Int_t igf1 = 0; igf1 < fV0Filters[1].GetEntriesFast() + 1; igf1++) {
    TObjArray * gamm1 = static_cast<TObjArray*>(fGammas.At(igf1));
    for(Int_t i1 = 0; i1 < gamm1->GetEntriesFast(); i1++) {
      AliAODConversionPhoton * ph1 = static_cast<AliAODConversionPhoton*>(gamm1->UncheckedAt(i1));
      Int_t tIDs[4] = {ph1->GetLabel(0), ph1->GetLabel(1), -1, -1};
      
      ///Combine gamma into pions
      Int_t igmax = 0;
      for(Int_t igf2 = 0; igf2 <= igf1; igf2++) {
	TObjArray * gamm2 = NULL;
	if(igf2 == igf1) {
	  gamm2 = gamm1;
	  igmax = i1;
	} else {
	  gamm2 = static_cast<TObjArray*>(fGammas.At(igf2));
	  igmax = gamm2->GetEntriesFast();
	}

	for(Int_t i2 = 0; i2 < igmax; i2++) {
	  AliAODConversionPhoton * ph2 = static_cast<AliAODConversionPhoton*>(gamm2->UncheckedAt(i2));
	  
	  if( ph2->GetTrackLabelPositive()==ph1->GetTrackLabelPositive() 
	      || ph2->GetTrackLabelNegative()==ph1->GetTrackLabelNegative()
	      || ph2->GetTrackLabelNegative()==ph1->GetTrackLabelPositive()
	      || ph2->GetTrackLabelPositive()==ph1->GetTrackLabelNegative()) {
	    continue;
	  }
	  
	  AliAODConversionMother * pion = new AliAODConversionMother(ph1, ph2);
	  if(pion->Eta() < etalim[0] || pion->Eta() > etalim[1]) continue;
	  pion->CalculateDistanceOfClossetApproachToPrimVtx(fInputEvent->GetPrimaryVertex());	  


	  tIDs[2] = ph2->GetLabel(0);
	  tIDs[3] = ph2->GetLabel(1);

	  massval[0] = pion->M();
	  massval[1] = pion->Pt();
	  massval[3] = igf1;
	  massval[4] = 0;

	  dphivalues[2] = pion->Pt();
	  dphivalues[6] = pion->M();
	  dphivalues[8] = igf1;
	  dphivalues[9] = 0;
	  
	  trigValues[0] = pion->Eta();
	  trigValues[1] = pion->Pt();
	  trigValues[4] = pion->M();
	  trigValues[5] = igf1; 
	  trigValues[6] = 0;
	  

	  ///Check that particle is in histo phase space
	  if(pion->Pt() > fAxistPt.GetBinLowEdge(1) && pion->Pt() < fAxistPt.GetXmax() &&
	     pion->Eta() > etalim[0] && pion->Eta() < etalim[1] &&
	     pion->M() > 0.04 && pion->M() < 0.4
	     ) {
	    
	    
	    if(fMesonFilter->MesonIsSelected(pion, kTRUE)) {

	      fMassSparse->Fill(massval);

	      Bool_t lpimap[fMesonFilters[0].GetEntriesFast()];
	      ///See if it passes lowside cuts
	      if(igf1 == 0 && igf2 == 0) {

		///Low side pion
		massval[3] = 0;
		for(Int_t ilpf = 0; ilpf < fMesonFilters[0].GetEntriesFast(); ilpf++) {
		  if(!(static_cast<AliConversionMesonCuts*>(fMesonFilters[0].At(ilpf))->MesonIsSelected(pion, kTRUE))) {
		    lpimap[ilpf] = kTRUE;
		    massval[4] = -(ilpf + 1);
		    fMassSparse->Fill(massval);
		  } else {
		    lpimap[ilpf] = kFALSE;
		  }
		  massval[4] = 0;
		}
		///Lowside v0
		for(Int_t iglf = 0; iglf < fV0Filters[0].GetEntriesFast(); iglf++) {
		  if(lowgmap[iglf][i1] || lowgmap[iglf][i2]){
		    massval[3] = -(iglf+1);
		    fMassSparse->Fill(massval);
		  }
		}
		massval[3] = 0;
	      }  /// End lowside mass histo fillers
	    	      
	      
	      ///Check that particle is in histo phase space
	      if(pion->Pt() > fAxistPt.GetBinLowEdge(1) && pion->Pt() < fAxistPt.GetXmax() &&
		 pion->M() > fAxisPiM.GetXmin() && pion->M() < fAxisPiM.GetXmax() &&
		 pion->Eta() > etalim[0] && pion->Eta() < etalim[1]) {
		
		
		const Int_t tbin = fAxistPt.FindFixBin(pion->Pt());

		///Fill standard triggers including upside v0 filters
		fTrigSparse->Fill(trigValues);

		
		if(igf1 == 0 && igf2 == 0) {
		  // piCorr->FillTriggerCounters(pion);
		  // piCorr->CorrelateWithTracks(pion, &tracks[0], tIDs, centrality, vertexz);
		  
		  
		  ////Only mix events with pion in signal region
		  hTrigPt->Fill(pion->Pt(), centrality, pion->M());
		  if(pion->M() > 0.1 && pion->M() < 0.15) {
		    hTrigPhi->Fill(pion->Phi());

		    ///Check trigger bin
		    if (tbin > 0 && tbin < (nbins + 1)) {
		      tmap[tbin-1] = kTRUE;
		    }
		  
		    ///Check if trigger also in low side (both gamma present in low side!)
		    for(Int_t ilgf = 0; ilgf < fV0Filters[0].GetEntriesFast(); ilgf++) {
		      if(!lowgmap[ilgf][i1] || !lowgmap[ilgf][i2]) {
			lv0tmap[tbin-1][ilgf] = kTRUE;
		      }
		    }
		  
		    ///See if the lowside pion filter also passes this, if not
		    for(Int_t ilpf = 0; ilpf < fMesonFilters[0].GetEntriesFast(); ilpf++) {
		      if(!lpimap[ilpf]) {
			lpitmap[tbin-1][ilpf] = kTRUE;
		      }
		    }
		  }

		} else {
		  if(pion->M() > 0.1 && pion->M() < 0.15) {
		    uv0tmap[tbin-1][igf1 - 1] = kTRUE;
		  }
		}
	      	    
	      
		///Fill the triggers not selected in lowside filters only if passsing standard v0 filter
		if(igf1 == 0 && igf2 == 0) {
		  
		  ///Lowside v0 filters
		  for(Int_t iglf = 0; iglf < fV0Filters[0].GetEntriesFast(); iglf++) {
		    if(lowgmap[iglf][i1] || lowgmap[iglf][i2]){
		      trigValues[5] = -(iglf+1);
		      fTrigSparse->Fill(trigValues);
		    }
		  }
		  
		  ////Low side pion filters
		  trigValues[5] = 0;
		  for(Int_t iplf = 0; iplf < fMesonFilters[0].GetEntriesFast(); iplf ++) {
		    if(lpimap[iplf]) {
		      trigValues[6] = -(iplf + 1);
		      fTrigSparse->Fill(trigValues);
		    }
		  }
		} // ifg1 == 0 
	      
		trigValues[5] = igf1; 
		trigValues[6] = 0;

		///////////////////////////////////////////////
		/// Correlate with tracks
		///////////////////////////////////////////////

		Int_t ntf = 1;
		if(igf1 == 0 && igf2 == 0) { 
		  ntf = fTrackFilters[1].GetEntriesFast() + 1;
		}

		for(Int_t itf = 0; itf < ntf; itf++) {
		  TObjArray * tracks = static_cast<TObjArray*>(fTracks.At(itf));
		  for(int ij = 0; ij < tracks->GetEntriesFast(); ij++) {
		    AliVTrack * track = static_cast<AliVTrack*>(tracks->At(ij));
		    Int_t tid = track->GetID();


		    if(tid == tIDs[0] || tid == tIDs[1] || tid == tIDs[2] || tid == tIDs[3])  {
		      continue;
		    }
		    
		    if(tid < 0) {
		      if(-tid-1 == tIDs[0]+1 || -tid-1 == tIDs[1]+1 || -tid-1 == tIDs[2]+1 || -tid-1 == tIDs[3]+1)  {
			continue;
		      }
		    }
		  
		  
		    dphivalues[0] = pion->Eta() - track->Eta();
		    dphivalues[1] = GetDPhi(pion->Phi() - track->Phi());
		    dphivalues[3] = track->Pt();
		    dphivalues[7] = itf;
		    dphivalues[8] = igf1;
		    dphivalues[9] = 0;
		    fCorrSparse->Fill(dphivalues, GetTrackCorrection(vertexz, track));
	
		    if(itf == 0 && igf1 == 0 && igf2 == 0) {
		      ///Fill the low side track filters
		      for(Int_t itlf = 0; itlf < fTrackFilters[0].GetEntriesFast(); itlf++) {
			if(lowtrackmap[itlf][ij]){
			  dphivalues[7] = -(itlf+1);
			  fCorrSparse->Fill(dphivalues, GetTrackCorrection(vertexz, track));
			}
		      }
		      ///Fill the low side v0 filters
		      dphivalues[7] = 0;
		      for(Int_t iglf = 0; iglf < fV0Filters[0].GetEntriesFast(); iglf++) {
			if(lowgmap[iglf][i1] || lowgmap[iglf][i2]){
			  dphivalues[8] = -(iglf+1);
			  fCorrSparse->Fill(dphivalues, GetTrackCorrection(vertexz, track));
			}
		      }

		      ///Fill the low side pi filter
		      dphivalues[7] = 0;
		      dphivalues[8] = 0;
		      for(Int_t iplf = 0; iplf < fMesonFilters[0].GetEntriesFast(); iplf ++) {
			if(lpimap[iplf]) {
			  dphivalues[9] = -(iplf + 1);
			  fCorrSparse->Fill(dphivalues, GetTrackCorrection(vertexz, track));
			}
		      }
		    }  /// end non standard filters track corr
		  
		  } // end for tracks
		} // end trackfilters loop
	      } //end check pion in histogram range to prevent overflow
	    } else {
	      /////////////////////////////
	      //// Not passing standard meson cuts, check upside filters	
	      ////////////////////////////
	     
	      ///Only check the pions from standard v0 filter 
	      if(igf1 == 0 && igf2 == 0) {
		for(Int_t ipuf = 0; ipuf < fMesonFilters[1].GetEntriesFast(); ipuf++) {
		  if(static_cast<AliConversionMesonCuts*>(fMesonFilters[1].At(ipuf))->MesonIsSelected(pion, kTRUE)) {
		    ///Fill invariant mass hist
		    massval[4] = (ipuf + 1);
		    fMassSparse->Fill(massval);
		    ///Check that particle is in histo phase space --- redundant!
		    
		    

		    if(pion->Pt() > fAxistPt.GetBinLowEdge(1) && pion->Pt() < fAxistPt.GetXmax() &&
		       pion->M() > fAxisPiM.GetXmin() && pion->M() < fAxisPiM.GetXmax() &&
		       pion->Eta() > etalim[0] && pion->Eta() < etalim[1]) {
		      
		    
		      ////Only mix events with pion in signal region
		      if(pion->M() > 0.1 && pion->M() < 0.15) {
			const Int_t tbin = fAxistPt.FindFixBin(pion->Pt());
			upitmap[tbin-1][ipuf] = kTRUE;
		      }
		    
		      ///Fill trigger counters
		      trigValues[6] = (ipuf + 1);
		      fTrigSparse->Fill(trigValues);
		      
		      ///Correlate with standard tracks
		      for(int ij = 0; ij < ttracks->GetEntriesFast(); ij++) {
			AliVTrack * track = static_cast<AliVTrack*>(ttracks->At(ij));
			Int_t tid = track->GetID();
			
			if(tid == tIDs[0] || tid == tIDs[1] || tid == tIDs[2] || tid == tIDs[3] ) {
			  continue;
			}

			if(tid < 0) {
			  if(-tid-1 == tIDs[0]+1 || -tid-1 == tIDs[1]+1 || -tid-1 == tIDs[2]+1 || -tid-1 == tIDs[3]+1)  {
			    continue;
			  }
			}
		  
			dphivalues[0] = pion->Eta() - track->Eta();
			dphivalues[1] = GetDPhi(pion->Phi() - track->Phi());
			dphivalues[3] = track->Pt();
			dphivalues[7] = 0; // track filter
			dphivalues[8] = 0; // v0 filter
			dphivalues[9] = ipuf + 1; // pion filter
			fCorrSparse->Fill(dphivalues, GetTrackCorrection(vertexz, track));
		      } /// end track corr
		    }
		  } // MesonIsSelected
		}
	      }
	    } /// end else ..  end upside meson filters
	    /////////////////////////////////////////////
	  } ///Etalim && pt limits


	} // i2 second gamma
      }
    } // i1 first gamma
  }
  //FillCounters(&pions, tracks, ntrackfilters, centrality, vertexz);
  
  ///////Fill track counters after entire event has been passed through
  ////
  
  Double_t trackValues[fTrackAxesList.GetSize()];
  trackValues[3] = centrality;
  trackValues[4] = vertexz;
  //trackValues[5] = particle->M(); remove !!!
  
  for(Int_t tbin = 0; tbin < fAxistPt.GetNbins(); tbin++) {
    trackValues[1] = fAxistPt.GetBinCenter(tbin+1);
    
    if(tmap[tbin]) {
      
      for(Int_t itf = 0; itf < fTrackFilters[1].GetEntriesFast() + 1; itf++) {
	TObjArray * tracks = static_cast<TObjArray*>(fTracks.At(itf));
	for(Int_t iTrack = 0; iTrack < tracks->GetEntriesFast(); iTrack++) {
	  AliVTrack * track = static_cast<AliVTrack*>(tracks->At(iTrack));
	  trackValues[0] = track->Eta();
	  trackValues[2] = track->Pt();
	  trackValues[5] = itf;
	  trackValues[6] = 0;  ///v0 filter
	  trackValues[7] = 0; ////Pi filter
	  fTrackSparse->Fill(trackValues, GetTrackCorrection(vertexz, track));

	  if(itf == 0) {
	    for(Int_t itlf = 0; itlf < fTrackFilters[0].GetEntriesFast(); itlf++) {
	      if(lowtrackmap[itlf][iTrack]) {
		trackValues[5] = -(itlf + 1);
		fTrackSparse->Fill(trackValues, GetTrackCorrection(vertexz, track) );
	      }
	    }
	    trackValues[5] = 0;
	    
	    ///Check lowside gamma
	    for(Int_t iglf = 0; iglf < fV0Filters[0].GetEntriesFast(); iglf++) {
	      if(!lv0tmap[tbin][iglf]) {
		trackValues[6] = -(iglf + 1);
		fTrackSparse->Fill(trackValues, GetTrackCorrection(vertexz, track));
	      }
	    }
	    trackValues[6] = 0;

	    ////Check lowside pion
	    for(Int_t iplf = 0; iplf < fMesonFilters[0].GetEntriesFast(); iplf++) {
	      if(!lpitmap[tbin][iplf]) {
		trackValues[7] = -(iplf + 1);
		fTrackSparse->Fill(trackValues, GetTrackCorrection(vertexz, track));
	      }
	    }



	  } // itf == 00
	}
      }
    } else {
      ///If not in main, see if in upside filters
      ///Do upside v0 filters
      for(Int_t iguf = 0; iguf < fV0Filters[1].GetEntriesFast(); iguf++) {
	if (uv0tmap[tbin][iguf] ) {

	  //cout << "c vtx " <<  centrality << vertexz << endl;

	  for(Int_t iTrack = 0; iTrack < ttracks->GetEntriesFast(); iTrack++) {
	    AliVTrack * track = static_cast<AliVTrack*>(ttracks->At(iTrack));
	    trackValues[0] = track->Eta();
	    trackValues[2] = track->Pt();
	    trackValues[5] = 0;
	    trackValues[6] = iguf+1;  ///v0 filter
	    trackValues[7] = 0; ////Pi filter
	    fTrackSparse->Fill(trackValues, GetTrackCorrection(vertexz, track));
	  }
	}
      }
      ///Do upside pi filter
      for(Int_t ipuf = 0; ipuf < fMesonFilters[1].GetEntriesFast(); ipuf++) {
	if (upitmap[tbin][ipuf] ) {

	  //cout << "c2 vtx2 " <<  centrality << vertexz << endl;


	  for(Int_t iTrack = 0; iTrack < ttracks->GetEntriesFast(); iTrack++) {
	    AliVTrack * track = static_cast<AliVTrack*>(ttracks->At(iTrack));
	    trackValues[0] = track->Eta();
	    trackValues[2] = track->Pt();
	    trackValues[5] = 0;
	    trackValues[6] = 0;  ///v0 filter
	    trackValues[7] = ipuf+1; ////Pi filter
	    fTrackSparse->Fill(trackValues, GetTrackCorrection(vertexz, track));
	  }
	}
      }
    }
  }



  //////
  //////


  
  //piCorr->FillCounters(&pions, tracks, centrality, vertexz);
  
  PostData(1, fHistograms);

}

//_______________________________________________________________________________
// void AliAnalysisTaskdPhi::FillCounters(TObjArray * particles, TObjArray tracks[], Int_t ntrackfilters, Float_t cent, Float_t vtxz) {
  
//   return;
  
  
//  //Fill ME Counters
//   const Int_t nbins = fAxistPt.GetNbins();
//   Bool_t tmap[nbins];
//   for(Int_t ptbin = 0; ptbin < nbins; ptbin++){
//       tmap[ptbin] = kFALSE;
//     }
//   }


//   Double_t trackValues[fTrackAxesList.GetSize()];
//   trackValues[3] = cent;
//   trackValues[4] = vtxz;

//   for(Int_t ip = 0; ip < particles->GetEntriesFast(); ip++){
//     AliAODConversionParticle * particle = static_cast<AliAODConversionParticle*>(particles->At(ip));

//     Int_t tbin = fAxistPt.FindFixBin(particle->Pt());
//     if (tbin > 0 && tbin < nbins + 1) {
//       if(tmap[tbin - 1] == kTRUE) {
//  	continue;
//       } else {
//  	tmap[tbin -1 ] = kTRUE;

//  	trackValues[5] = particle->M();
//  	trackValues[1] = particle->Pt();
	
//  	for(Int_t itf = 0; itf < ntrackfilters; itf++) {
//  	  if(fkTrackAxis) trackValues[6] = itf;
//  	  for(int ij = 0; ij < tracks[itf].GetEntriesFast(); ij++) {
//  	    AliVTrack * track = static_cast<AliVTrack*>(tracks->UncheckedAt(ij));
//  	    trackValues[0] = track->Eta();
//  	    trackValues[2] = track->Pt();
//  	    fTrackSparse->Fill(trackValues);	
//  	  }
//  	}
//       }
//     }
//   }
// }




// //________________________________________________________________
// void AliAnalysisTaskdPhi::CorrelateWithTracks(AliAODConversionParticle * particle, TObjArray tracks[], Int_t ntrackfilters, Bool_t ** lowtrackmap, Int_t nltf, Int_t const tIDs[4], Double_t dphivalues[]) {
//   //Correlate particle with tracks
//   ///Correlate with tracks

// }

//________________________________________________________________________
void AliAnalysisTaskdPhi::Terminate(Option_t *) {
 
  // Draw result to the screen
  // Called once at the end of the query
}

//________________________________________________________________________
TClonesArray * AliAnalysisTaskdPhi::GetConversionGammas(Bool_t isAOD) {
  if(isAOD) {

    TClonesArray * gammas = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(fDeltaAODBranchName.Data()));
    if(gammas) {
      return gammas;
    }
    
    FindDeltaAODBranchName(fInputEvent);
    gammas = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(fDeltaAODBranchName.Data()));
    return gammas;
    
  } else {
    TClonesArray * gammas = dynamic_cast<TClonesArray*>(GetInputData(1));
    return gammas;
  }
  
}

//________________________________________________________________________
void AliAnalysisTaskdPhi::FindDeltaAODBranchName(AliVEvent * event){
  ///Find aod branch
  TList *list=event->GetList();
  for(Int_t ii=0;ii<list->GetEntries();ii++){
    TString name((list->At(ii))->GetName());
    if(name.BeginsWith("GammaConv")&&name.EndsWith("gamma")){
      fDeltaAODBranchName=name;
      AliDebug(AliLog::kDebug + 5, Form("Set DeltaAOD BranchName to: %s",fDeltaAODBranchName.Data()));
      return;
    }
  }
}
  

//________________________________________________________________________
Double_t AliAnalysisTaskdPhi::GetTrackCorrection(Double_t vtxz, AliVTrack * track) {
  ////Get track correction from map
  Int_t coord[4] = {-1, -1, -1, -1};
  if(fCorrectionMap) {
    Double_t values[4] = { vtxz, track->Pt(), track->Eta(), track->Phi() };
    Double_t correction = fCorrectionMap->GetBinContent(fCorrectionMap->GetBin(values, kFALSE), coord);
    if (fCorrectionMap->IsInRange(coord)) {
      return correction;
    } 
  }
  return 1.0;
}