////////////////////////////////////////////////////////////////////////////////
///                                                                          ///
/// AliFemtoEventReaderESDChain - the reader class for the Alice ESD         ///
/// tailored for the Task framework                                 ///
/// Reads in AliESDfriend to create shared hit/quality information           ///
/// Authors: Adam Kisiel kisiel@mps.ohio-state.edu                           ///
///                                                                          ///
////////////////////////////////////////////////////////////////////////////////
#include "AliFemtoEventReaderESDChain.h"

#include "TFile.h"
#include "TTree.h"
#include "AliESDEvent.h"
#include "AliESDtrack.h"
#include "AliESDVertex.h"
#include "AliMultiplicity.h"
#include "AliCentrality.h"
#include "AliEventplane.h"
#include "AliESDVZERO.h"
#include "AliFmPhysicalHelixD.h"
#include "AliFmThreeVectorF.h"
#include "SystemOfUnits.h"
#include "AliFemtoEvent.h"
#include "AliFemtoModelHiddenInfo.h"
#include "AliPID.h"
#include "AliAnalysisUtils.h"

ClassImp(AliFemtoEventReaderESDChain)

#if !(ST_NO_NAMESPACES)
  using namespace units;
#endif

using namespace std;

//____________________________
AliFemtoEventReaderESDChain::AliFemtoEventReaderESDChain():
  fFileName(" "),
  fConstrained(true),
  fReadInner(false),
  fUseTPCOnly(false),
  fNumberofEvent(0),
  fCurEvent(0),
  fCurFile(0),
  fEvent(0x0),
  fTrackType(kGlobal),
  fEstEventMult(kReferenceITSTPC), 
  fEventTrig(AliVEvent::kMB), //trigger
  fESDpid(0),
  fIsPidOwner(0),
  fReadV0(0),
  fMagFieldSign(0),
  fpA2013(kFALSE),
  fisPileUp(kFALSE),
  fMVPlp(kFALSE),
  fMinVtxContr(0),
  fMinPlpContribMV(0),
  fMinPlpContribSPD(0)
{
  //constructor with 0 parameters , look at default settings 
  //   fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  //     fClusterPerPadrow[tPad] = new list<Int_t>();
  //   }
  //   fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  //     fSharedList[tPad] = new list<Int_t>();
  //   }
}

//__________________
AliFemtoEventReaderESDChain::AliFemtoEventReaderESDChain(const AliFemtoEventReaderESDChain& aReader):
  AliFemtoEventReader(aReader),
  fFileName(" "),
  fConstrained(true),
  fReadInner(false),
  fUseTPCOnly(false),
  fNumberofEvent(0),
  fCurEvent(0),
  fCurFile(0),
  fEvent(0x0),
  fTrackType(kGlobal),
  fEstEventMult(kReferenceITSTPC),
  fEventTrig(AliVEvent::kMB), //trigger
  fESDpid(0),
  fIsPidOwner(0),
  fReadV0(0),
  fMagFieldSign(0),
  fpA2013(kFALSE),
  fisPileUp(kFALSE),
  fMVPlp(kFALSE),
  fMinVtxContr(0),
  fMinPlpContribMV(0),
  fMinPlpContribSPD(0)
{
  // Copy constructor
  fConstrained = aReader.fConstrained;
  fReadInner = aReader.fReadInner;
  fUseTPCOnly = aReader.fUseTPCOnly;
  fNumberofEvent = aReader.fNumberofEvent;
  fCurEvent = aReader.fCurEvent;
  fCurFile = aReader.fCurFile;
  //  fEvent = new AliESD(*aReader.fEvent);
  fEvent = new AliESDEvent();
  fTrackType = aReader.fTrackType;
  fEstEventMult = aReader.fEstEventMult;
  fEventTrig = aReader.fEventTrig; //trigger
  fReadV0 = aReader.fReadV0;
  fMagFieldSign = aReader.fMagFieldSign;
  fpA2013 = aReader.fpA2013;
  fisPileUp = aReader.fisPileUp;
  fMVPlp = aReader.fMVPlp;
  fMinVtxContr =  aReader.fMinVtxContr;
  fMinPlpContribMV  =  aReader.fMinPlpContribMV;
  fMinPlpContribSPD  =  aReader.fMinPlpContribSPD;

  //   fEventFriend = aReader.fEventFriend;
  //   fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  //     fClusterPerPadrow[tPad] = new list<Int_t>();
  //     list<Int_t>::iterator iter;
  //     for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) {
  //       fClusterPerPadrow[tPad]->push_back(*iter);
  //     }
  //   }
  //   fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  //     fSharedList[tPad] = new list<Int_t>();
  //     list<Int_t>::iterator iter;
  //     for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) {
  //       fSharedList[tPad]->push_back(*iter);
  //     }
  //   }
}
//__________________
AliFemtoEventReaderESDChain::~AliFemtoEventReaderESDChain()
{
  //Destructor
  delete fEvent;

  //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  //     fClusterPerPadrow[tPad]->clear();
  //     delete fClusterPerPadrow[tPad];
  //   }
  //   delete [] fClusterPerPadrow;
  //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  //     fSharedList[tPad]->clear();
  //     delete fSharedList[tPad];
  //   }
  //   delete [] fSharedList;

}

//__________________
AliFemtoEventReaderESDChain& AliFemtoEventReaderESDChain::operator=(const AliFemtoEventReaderESDChain& aReader)
{
  // Assignment operator
  if (this == &aReader)
    return *this;

  fConstrained = aReader.fConstrained;
  fReadInner = aReader.fReadInner;
  fUseTPCOnly = aReader.fUseTPCOnly;
  fNumberofEvent = aReader.fNumberofEvent;
  fCurEvent = aReader.fCurEvent;
  fCurFile = aReader.fCurFile;
  if (fEvent) delete fEvent;
  fEvent = new AliESDEvent();
  fTrackType = aReader.fTrackType;
  fEstEventMult = aReader.fEstEventMult;

  fReadV0 = aReader.fReadV0;
  fMagFieldSign = aReader.fMagFieldSign;
  fpA2013 = aReader.fpA2013;
  fisPileUp = aReader.fisPileUp;
  fMVPlp = aReader.fMVPlp;
  //  fEventFriend = aReader.fEventFriend;
  
  //   if (fClusterPerPadrow) {
  //     for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  //       fClusterPerPadrow[tPad]->clear();
  //       delete fClusterPerPadrow[tPad];
  //     }
  //     delete [] fClusterPerPadrow;
  //   }
  
  //   if (fSharedList) {
  //     for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  //       fSharedList[tPad]->clear();
  //       delete fSharedList[tPad];
  //     }
  //     delete [] fSharedList;
  //   }

  //   fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  //     fClusterPerPadrow[tPad] = new list<Int_t>();
  //     list<Int_t>::iterator iter;
  //     for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) {
  //       fClusterPerPadrow[tPad]->push_back(*iter);
  //     }
  //   }
  //   fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  //     fSharedList[tPad] = new list<Int_t>();
  //     list<Int_t>::iterator iter;
  //     for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) {
  //       fSharedList[tPad]->push_back(*iter);
  //     }
  //   }
  
  return *this;
}
//__________________
// Simple report
AliFemtoString AliFemtoEventReaderESDChain::Report()
{
  AliFemtoString temp = "\n This is the AliFemtoEventReaderESDChain\n";
  return temp;
}

//__________________
void AliFemtoEventReaderESDChain::SetConstrained(const bool constrained)
{
  // Select whether to read constrained or not constrained momentum
  fConstrained=constrained;
}
//__________________
bool AliFemtoEventReaderESDChain::GetConstrained() const
{
  // Check whether we read constrained or not constrained momentum
  return fConstrained;
}
//__________________
void AliFemtoEventReaderESDChain::SetReadTPCInner(const bool readinner)
{
  fReadInner=readinner;
}

bool AliFemtoEventReaderESDChain::GetReadTPCInner() const
{
  return fReadInner;
}

//__________________
void AliFemtoEventReaderESDChain::SetUseTPCOnly(const bool usetpconly)
{
  fUseTPCOnly=usetpconly;
}

bool AliFemtoEventReaderESDChain::GetUseTPCOnly() const
{
  return fUseTPCOnly;
}

void AliFemtoEventReaderESDChain::SetUseMultiplicity(EstEventMult aType)
{
  fEstEventMult = aType;
}

AliFemtoEvent* AliFemtoEventReaderESDChain::ReturnHbtEvent()
{
  // Get the event, read all the relevant information
  // and fill the AliFemtoEvent class
  // Returns a valid AliFemtoEvent
  AliFemtoEvent *hbtEvent = 0;
 

  hbtEvent = new AliFemtoEvent;

  CopyESDtoFemtoEvent(hbtEvent);

  fCurEvent++;	


  return hbtEvent; 
}

void AliFemtoEventReaderESDChain::CopyESDtoFemtoEvent(AliFemtoEvent *hbtEvent)
{
  

 //string tFriendFileName;
  // Get the friend information
  if (Debug()>1) cout<<"starting to read event "<<fCurEvent<<endl;
  //  fEvent->SetESDfriend(fEventFriend);
  if(fEvent->GetAliESDOld())fEvent->CopyFromOldESD();
  


   //setting basic things
  //  hbtEvent->SetEventNumber(fEvent->GetEventNumber());
  hbtEvent->SetRunNumber(fEvent->GetRunNumber());
  //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
  hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
  hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
  hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
  hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
  hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
  hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
  hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
  hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
  //  hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());

  if ((fEvent->IsTriggerClassFired("CINT1WU-B-NOPF-ALL")) ||
      (fEvent->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")) ||
      (fEvent->IsTriggerClassFired("CINT1-B-NOPF-ALLNOTRD")) ||
      (fEvent->IsTriggerClassFired("CINT1-B-NOPF-FASTNOTRD")))
    hbtEvent->SetTriggerCluster(1);
  else if ((fEvent->IsTriggerClassFired("CSH1WU-B-NOPF-ALL")) ||
	   (fEvent->IsTriggerClassFired("CSH1-B-NOPF-ALLNOTRD")))
    hbtEvent->SetTriggerCluster(2);
  else 
    hbtEvent->SetTriggerCluster(0);
	
  //Vertex
  double fV1[3];
  double fVCov[6];

  //AliAnalysisUtils
  if(fisPileUp||fpA2013)
    {
      AliAnalysisUtils *anaUtil=new AliAnalysisUtils();
      if(fMinVtxContr)
	anaUtil->SetMinVtxContr(fMinVtxContr);
      if(fpA2013)
	if(anaUtil->IsVertexSelected2013pA(fEvent)==kFALSE) return; //Vertex rejection for pA analysis.
      if(fMVPlp) anaUtil->SetUseMVPlpSelection(kTRUE);
      else anaUtil->SetUseMVPlpSelection(kFALSE);
      if(fMinPlpContribMV) anaUtil->SetMinPlpContribMV(fMinPlpContribMV);
      if(fMinPlpContribSPD) anaUtil->SetMinPlpContribSPD(fMinPlpContribSPD);
      if(fisPileUp)
	if(anaUtil->IsPileUpEvent(fEvent)) return; //Pile-up rejection.
      delete anaUtil;   
    }


  if (fUseTPCOnly) {
    const AliESDVertex* esdvertex = (AliESDVertex*) fEvent->GetPrimaryVertexTPC();
    if(!esdvertex || esdvertex->GetNContributors() < 1) return; //Bad vertex, skip event.

    esdvertex->GetXYZ(fV1);
    esdvertex->GetCovMatrix(fVCov);
    if (!esdvertex->GetStatus())
      fVCov[4] = -1001.0;
  }
  else {
    const AliESDVertex* esdvertex = (AliESDVertex*) fEvent->GetPrimaryVertex();
    if(!esdvertex || esdvertex->GetNContributors() < 1) return; //Bad vertex, skip event.
   
    esdvertex->GetXYZ(fV1);
    esdvertex->GetCovMatrix(fVCov);
    if (!esdvertex->GetStatus())
      fVCov[4] = -1001.0;
  }
    
  AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
  hbtEvent->SetPrimVertPos(vertex);
  hbtEvent->SetPrimVertCov(fVCov);
  
  Int_t spdetaonecount = 0;
  
  for (int iter=0; iter<fEvent->GetMultiplicity()->GetNumberOfTracklets(); iter++) 
    if (fabs(fEvent->GetMultiplicity()->GetEta(iter)) < 1.0)
      spdetaonecount++;

  //  hbtEvent->SetSPDMult(fEvent->GetMultiplicity()->GetNumberOfTracklets());
  hbtEvent->SetSPDMult(spdetaonecount);

  //starting to reading tracks
  int nofTracks=0;  //number of reconstructed tracks in event
  nofTracks=fEvent->GetNumberOfTracks();
  int realnofTracks=0;//number of track which we use ina analysis

  //   // Clear the shared cluster list
  //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  //     fClusterPerPadrow[tPad]->clear();
  //   }
  //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  //     fSharedList[tPad]->clear();
  //   }


  //   for (int i=0;i<nofTracks;i++) {
  //     const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track

  //     list<Int_t>::iterator tClustIter;

  //     Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
  //     Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
  //     for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
  //       if (tTrackIndices[tNcl] >= 0) {
  // 	tClustIter = find(fClusterPerPadrow[tNcl]->begin(), fClusterPerPadrow[tNcl]->end(), tTrackIndices[tNcl]);
  // 	if (tClustIter == fClusterPerPadrow[tNcl]->end()) {
  // 	  fClusterPerPadrow[tNcl]->push_back(tTrackIndices[tNcl]);
  // 	}
  // 	else {
  // 	  fSharedList[tNcl]->push_back(tTrackIndices[tNcl]);
  // 	}
  //       }
  //     }
      
  //   }

  int tNormMult = 0;
  int tNormMultPos = 0;
  int tNormMultNeg = 0;

  Float_t tTotalPt = 0.0;

  Float_t b[2];
  Float_t bCov[3];

  Int_t tTracklet=0, tITSTPC=0;
  
  //W-AliESDEvent::EstimateMultiplicity: This obsolete method will be eliminated soon. Use AliESDtrackCuts::GetReferenceMultiplicity
  //fEvent->EstimateMultiplicity(tTracklet, tITSTPC, tITSPure, 1.2);
  
  hbtEvent->SetMultiplicityEstimateITSTPC(tITSTPC);
  hbtEvent->SetMultiplicityEstimateTracklets(tTracklet);
  //  hbtEvent->SetMultiplicityEstimateITSPure(tITSPure);
  hbtEvent->SetMultiplicityEstimateITSPure(fEvent->GetMultiplicity()->GetNumberOfITSClusters(1));
  
  for (int i=0;i<nofTracks;i++)
    {
      bool  tGoodMomentum=true; //flaga to chcek if we can read momentum of this track

      const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
      //      const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();

      if ((esdtrack->GetStatus() & AliESDtrack::kTPCrefit) &&
	  (esdtrack->GetStatus() & AliESDtrack::kITSrefit)) {
	if (esdtrack->GetTPCNcls() > 70) 
	  if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 4.0) {
	    if (esdtrack->Pt() > 0.15 && esdtrack->Pt() < 20) 
	      if (TMath::Abs(esdtrack->Eta()) < 0.8) {
		esdtrack->GetImpactParameters(b,bCov);
		if ((b[0]<0.2) && (b[1] < 2.0)) {
		  tNormMult++;
		  tTotalPt += esdtrack->Pt();
		}
	      }
	  }
      }

      hbtEvent->SetZDCEMEnergy(tTotalPt);
      //       if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit)
      // 	if (esdtrack->GetTPCNcls() > 80) 
      // 	  if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 6.0) 
      // 	    if (esdtrack->GetConstrainedParam())
      // 	      if (fabs(esdtrack->GetConstrainedParam()->Eta()) < 0.5)
      // 		if (esdtrack->GetConstrainedParam()->Pt() < 1.0) {
      // 		  if (esdtrack->GetSign() > 0)
      // 		    tNormMultPos++;
      // 		  else if (esdtrack->GetSign() < 0)
      // 		    tNormMultNeg--;
      // 		}

      // If reading ITS-only tracks, reject all with TPC
      if (fTrackType == kITSOnly) {
	if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit) continue;
	if (!(esdtrack->GetStatus() & AliESDtrack::kITSrefit)) continue;
	if (esdtrack->GetStatus() & AliESDtrack::kTPCin) continue;
	UChar_t iclm = esdtrack->GetITSClusterMap();
	Int_t incls = 0;
	for (int iter=0; iter<6; iter++) if (iclm&(1<<iter)) incls++;
	if (incls<=3) {
	  if (Debug()>1) cout << "Rejecting track with " << incls << " clusters" << endl;
	  continue;
	}
      }

      AliFemtoTrack* trackCopy = new AliFemtoTrack();	
      trackCopy->SetCharge((short)esdtrack->GetSign());

      //in aliroot we have AliPID 
      //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon   
      //we use only 5 first
      double esdpid[5];
      //       esdtrack->GetESDpid(esdpid);
      esdtrack->GetTPCpid(esdpid);
      trackCopy->SetPidProbElectron(esdpid[0]);
      trackCopy->SetPidProbMuon(esdpid[1]);
      trackCopy->SetPidProbPion(esdpid[2]);
      trackCopy->SetPidProbKaon(esdpid[3]);
      trackCopy->SetPidProbProton(esdpid[4]);

      esdpid[0] = -100000.0;
      esdpid[1] = -100000.0;
      esdpid[2] = -100000.0;
      esdpid[3] = -100000.0;
      esdpid[4] = -100000.0;

      double tTOF = 0.0;

      if (esdtrack->GetStatus()&AliESDtrack::kTOFout&AliESDtrack::kTIME) {
	tTOF = esdtrack->GetTOFsignal();
	esdtrack->GetIntegratedTimes(esdpid);
      }

      trackCopy->SetTofExpectedTimes(tTOF-esdpid[2], tTOF-esdpid[3], tTOF-esdpid[4]);

      //////  TPC ////////////////////////////////////////////

      float nsigmaTPCK=-1000.;
      float nsigmaTPCPi=-1000.;
      float nsigmaTPCP=-1000.;


      if ((fESDpid) && (esdtrack->IsOn(AliESDtrack::kTPCpid))){
        nsigmaTPCK = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kKaon);
        nsigmaTPCPi = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kPion);
        nsigmaTPCP = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kProton);

      }
      trackCopy->SetNSigmaTPCPi(nsigmaTPCPi);
      trackCopy->SetNSigmaTPCK(nsigmaTPCK);
      trackCopy->SetNSigmaTPCP(nsigmaTPCP);

      ///// TOF ///////////////////////////////////////////////

	float vp=-1000.;
	float nsigmaTOFPi=-1000.;
	float nsigmaTOFK=-1000.;
	float nsigmaTOFP=-1000.;

	if (// (esdtrack->GetStatus()&AliESDtrack::kTOFpid) &&
	    (esdtrack->GetStatus()&AliESDtrack::kTOFout) &&
	    (esdtrack->GetStatus()&AliESDtrack::kTIME))
	  {

	    //if ((esdtrack->GetStatus()&AliESDtrack::kTOFpid) &&
	    //(esdtrack->GetStatus()&AliESDtrack::kTOFout) &&
	    //(esdtrack->GetStatus()&AliESDtrack::kTIME)){
	    // collect info from ESDpid class

	    if ((fESDpid) && (esdtrack->IsOn(AliESDtrack::kTOFout & AliESDtrack::kTIME))) {


	      double tZero = fESDpid->GetTOFResponse().GetStartTime(esdtrack->P());

	      nsigmaTOFPi = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kPion,tZero);
	      nsigmaTOFK = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kKaon,tZero);
	      nsigmaTOFP = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kProton,tZero);

	      Double_t len=esdtrack->GetIntegratedLength();
	      Double_t tof=esdtrack->GetTOFsignal();
	      if(tof > 0.) vp=len/tof/0.03;
	    }
	  }

	trackCopy->SetVTOF(vp);
	trackCopy->SetNSigmaTOFPi(nsigmaTOFPi);
	trackCopy->SetNSigmaTOFK(nsigmaTOFK);
	trackCopy->SetNSigmaTOFP(nsigmaTOFP);

	double pxyz[3];
	double rxyz[3];
	double impact[2];
	double covimpact[3];

	if (fUseTPCOnly) {
	  if (!esdtrack->GetTPCInnerParam()) {
	    delete trackCopy;
	    continue;
	  }


	  AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
	  param->GetXYZ(rxyz);
	  param->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impact, covimpact);
	  param->GetPxPyPz(pxyz);//reading noconstarined momentum

	  if (fReadInner == true) {
	    AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
	    tInfo->SetPDGPid(211);
	    tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
	    tInfo->SetMass(0.13957);
	    //	  tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0);
	    //	  tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0);
	    tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
	    trackCopy->SetHiddenInfo(tInfo);
	  }
	
	  AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
	  if (v.Mag() < 0.0001) {
	    //	cout << "Found 0 momentum ???? " <<endl;
	    delete trackCopy;
	    continue;
	  }
	  trackCopy->SetP(v);//setting momentum
	  trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));

	  const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
	  const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
	  //setting helix I do not if it is ok
	  AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge())); 
	  trackCopy->SetHelix(helix);

	  //some stuff which could be useful 
	  trackCopy->SetImpactD(impact[0]);
	  trackCopy->SetImpactZ(impact[1]);
	  trackCopy->SetCdd(covimpact[0]);
	  trackCopy->SetCdz(covimpact[1]);
	  trackCopy->SetCzz(covimpact[2]);
	  trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact, covimpact));	

	  delete param;
	}
	else {
	  if (fReadInner == true) {
	  
	    if (esdtrack->GetTPCInnerParam()) {
	      AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetInnerParam());
	      //trackCopy->SetInnerMomentum(param->P());
	      trackCopy->SetInnerMomentum(esdtrack->GetTPCmomentum());
	      param->GetXYZ(rxyz);
	      //	    param->PropagateToDCA(fEvent->GetPrimaryVertex(), (fEvent->GetMagneticField()), 10000);
	      param->GetPxPyPz(pxyz);//reading noconstarined momentum
	      delete param;
	    
	      AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
	      tInfo->SetPDGPid(211);
	      tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
	      tInfo->SetMass(0.13957);
	      //	    tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0);
	      //tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0);
	      tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
	      trackCopy->SetHiddenInfo(tInfo);
	    }
	  }

	  if (fTrackType == kGlobal) {
	    if (fConstrained==true)		    
	      tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
	    else
	      tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
	  }
	  else if (fTrackType == kTPCOnly) {
	    if (esdtrack->GetTPCInnerParam())
	      esdtrack->GetTPCInnerParam()->GetPxPyPz(pxyz);
	    else {
	      delete trackCopy;
	      continue;
	    }
	  }
	  else if (fTrackType == kITSOnly) {
	    if (fConstrained==true)		    
	      tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
	    else
	      tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
	  }


	  AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
	  if (v.Mag() < 0.0001) {
	    //	cout << "Found 0 momentum ???? " <<endl;
	    delete trackCopy;
	    continue;
	  }
	  trackCopy->SetP(v);//setting momentum
	  trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
	  const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
	  const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
	  //setting helix I do not if it is ok
	  AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge())); 
	  trackCopy->SetHelix(helix);

	  //some stuff which could be useful 
	  float imp[2];
	  float cim[3];
	  // if (fTrackType == kTPCOnly) {
	  //   esdtrack->GetTPCInnerParam()->GetImpactParameters(imp,cim);
	  // }
	  // else {
	  esdtrack->GetImpactParameters(imp,cim);
	  // }

	  impact[0] = imp[0];
	  impact[1] = imp[1];
	  covimpact[0] = cim[0];
	  covimpact[1] = cim[1];
	  covimpact[2] = cim[2];

	  trackCopy->SetImpactD(impact[0]);
	  trackCopy->SetImpactZ(impact[1]);
	  trackCopy->SetCdd(covimpact[0]);
	  trackCopy->SetCdz(covimpact[1]);
	  trackCopy->SetCzz(covimpact[2]);
	  trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact,covimpact));
	}

	trackCopy->SetTrackId(esdtrack->GetID());
	trackCopy->SetFlags(esdtrack->GetStatus());
	trackCopy->SetLabel(esdtrack->GetLabel());
		
	trackCopy->SetITSchi2(esdtrack->GetITSchi2());    
	if (esdtrack->GetITSFakeFlag())
	  trackCopy->SetITSncls(-esdtrack->GetNcls(0));     
	else
	  trackCopy->SetITSncls(esdtrack->GetNcls(0));     
	trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());       
	trackCopy->SetTPCncls(esdtrack->GetTPCNcls());       
	trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());      
	trackCopy->SetTPCsignal(esdtrack->GetTPCsignal());
	trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class   
	trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma()); 

	trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
	trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());

	double xtpc[3];
	esdtrack->GetInnerXYZ(xtpc);
	xtpc[2] -= fV1[2];
	trackCopy->SetNominalTPCEntrancePoint(xtpc);

	esdtrack->GetOuterXYZ(xtpc);
	xtpc[2] -= fV1[2];
	trackCopy->SetNominalTPCExitPoint(xtpc);

	int indexes[3];
	for (int ik=0; ik<3; ik++) {
	  indexes[ik] = esdtrack->GetKinkIndex(ik);
	}
	trackCopy->SetKinkIndexes(indexes);

	for (int ii=0; ii<6; ii++){
	  trackCopy->SetITSHitOnLayer(ii,esdtrack->HasPointOnITSLayer(ii));
	}

	//decision if we want this track
	//if we using diffrent labels we want that this label was use for first time 
	//if we use hidden info we want to have match between sim data and ESD
	if (tGoodMomentum==true)
	  {
	    hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
	    realnofTracks++;//real number of tracks
	    //	  delete trackCopy;
	  }
	else
	  {
	    delete  trackCopy;
	  }
		
    }

  hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event	

  AliCentrality *cent = fEvent->GetCentrality();
  if (cent) {
    hbtEvent->SetCentralityV0(cent->GetCentralityPercentile("V0M"));
    hbtEvent->SetCentralityV0A(cent->GetCentralityPercentile("V0A"));
    hbtEvent->SetCentralityV0C(cent->GetCentralityPercentile("V0C"));
    hbtEvent->SetCentralityZNA(cent->GetCentralityPercentile("ZNA"));
    hbtEvent->SetCentralityZNC(cent->GetCentralityPercentile("ZNC"));
    hbtEvent->SetCentralityCL1(cent->GetCentralityPercentile("CL1"));
    hbtEvent->SetCentralityCL0(cent->GetCentralityPercentile("CL0"));
    hbtEvent->SetCentralityTKL(cent->GetCentralityPercentile("TKL"));
    hbtEvent->SetCentralityFMD(cent->GetCentralityPercentile("FMD"));
    hbtEvent->SetCentralityFMD(cent->GetCentralityPercentile("NPA"));
    //    hbtEvent->SetCentralitySPD1(cent->GetCentralityPercentile("CL1"));
    hbtEvent->SetCentralityTrk(cent->GetCentralityPercentile("TRK"));
    hbtEvent->SetCentralityCND(cent->GetCentralityPercentile("CND"));

    if (Debug()>1) printf("  FemtoReader Got Event with %f %f %f %f\n", cent->GetCentralityPercentile("V0M"), 0.0, cent->GetCentralityPercentile("CL1"), 0.0);
  }

  if (fEstEventMult == kGlobalCount) 
    hbtEvent->SetNormalizedMult(tNormMult);
  else if (fEstEventMult == kReferenceITSTPC)
    hbtEvent->SetNormalizedMult(AliESDtrackCuts::GetReferenceMultiplicity(fEvent,AliESDtrackCuts::kTrackletsITSTPC,1.0));
  else if(fEstEventMult == kReferenceITSSA)
    hbtEvent->SetNormalizedMult(AliESDtrackCuts::GetReferenceMultiplicity(fEvent,AliESDtrackCuts::kTrackletsITSSA,1.0));
  else if(fEstEventMult == kReferenceTracklets)
    hbtEvent->SetNormalizedMult(AliESDtrackCuts::GetReferenceMultiplicity(fEvent,AliESDtrackCuts::kTracklets,1.0));
  else if (fEstEventMult == kSPDLayer1)
    hbtEvent->SetNormalizedMult(fEvent->GetMultiplicity()->GetNumberOfITSClusters(1));
  else if (fEstEventMult == kVZERO)
    {
      Float_t multV0 = 0;
      for (Int_t i=0; i<64; i++)
	multV0 += fEvent->GetVZEROData()->GetMultiplicity(i);
      hbtEvent->SetNormalizedMult(multV0);
    }
  else if (fEstEventMult == kCentrality) {
    // centrality between 0 (central) and 1 (very peripheral)

    if (cent) {
      if (cent->GetCentralityPercentile("V0M") < 0.00001)
	hbtEvent->SetNormalizedMult(-1);
      else
	hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("V0M")));
      if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(), 
			     10.0*cent->GetCentralityPercentile("V0M"), lrint(10.0*cent->GetCentralityPercentile("V0M")));
    }
  }
  else if (fEstEventMult == kCentralityV0A) {
    // centrality between 0 (central) and 1 (very peripheral)

    if (cent) {
      if (cent->GetCentralityPercentile("V0A") < 0.00001)
	hbtEvent->SetNormalizedMult(-1);
      else
	hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("V0A")));
      if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(), 
			     10.0*cent->GetCentralityPercentile("V0A"), lrint(10.0*cent->GetCentralityPercentile("V0A")));
    }
  }
  else if (fEstEventMult == kCentralityV0C) {
    // centrality between 0 (central) and 1 (very peripheral)

    if (cent) {
      if (cent->GetCentralityPercentile("V0C") < 0.00001)
	hbtEvent->SetNormalizedMult(-1);
      else
	hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("V0C")));
      if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(), 
			     10.0*cent->GetCentralityPercentile("V0C"), lrint(10.0*cent->GetCentralityPercentile("V0C")));
    }
  }
  else if (fEstEventMult == kCentralityZNA) {
    // centrality between 0 (central) and 1 (very peripheral)

    if (cent) {
      if (cent->GetCentralityPercentile("ZNA") < 0.00001)
	hbtEvent->SetNormalizedMult(-1);
      else
	hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("ZNA")));
      if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(), 
			     10.0*cent->GetCentralityPercentile("ZNA"), lrint(10.0*cent->GetCentralityPercentile("ZNA")));
    }
  }
  else if (fEstEventMult == kCentralityZNC) {
    // centrality between 0 (central) and 1 (very peripheral)

    if (cent) {
      if (cent->GetCentralityPercentile("ZNC") < 0.00001)
	hbtEvent->SetNormalizedMult(-1);
      else
	hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("ZNC")));
      if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(), 
			     10.0*cent->GetCentralityPercentile("ZNC"), lrint(10.0*cent->GetCentralityPercentile("ZNC")));
    }
  }
  else if (fEstEventMult == kCentralityCL1) {
    // centrality between 0 (central) and 1 (very peripheral)

    if (cent) {
      if (cent->GetCentralityPercentile("CL1") < 0.00001)
	hbtEvent->SetNormalizedMult(-1);
      else
	hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("CL1")));
      if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(), 
			     10.0*cent->GetCentralityPercentile("CL1"), lrint(10.0*cent->GetCentralityPercentile("CL1")));
    }
  }
  else if (fEstEventMult == kCentralityCL0) {
    // centrality between 0 (central) and 1 (very peripheral)

    if (cent) {
      if (cent->GetCentralityPercentile("CL0") < 0.00001)
	hbtEvent->SetNormalizedMult(-1);
      else
	hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("CL0")));
      if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(), 
			     10.0*cent->GetCentralityPercentile("CL0"), lrint(10.0*cent->GetCentralityPercentile("CL0")));
    }
  }
  else if (fEstEventMult == kCentralityTRK) {
    // centrality between 0 (central) and 1 (very peripheral)

    if (cent) {
      if (cent->GetCentralityPercentile("TRK") < 0.00001)
	hbtEvent->SetNormalizedMult(-1);
      else
	hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("TRK")));
      if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(), 
			     10.0*cent->GetCentralityPercentile("TRK"), lrint(10.0*cent->GetCentralityPercentile("TRK")));
    }
  }
  else if (fEstEventMult == kCentralityTKL) {
    // centrality between 0 (central) and 1 (very peripheral)

    if (cent) {
      if (cent->GetCentralityPercentile("TKL") < 0.00001)
	hbtEvent->SetNormalizedMult(-1);
      else
	hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("TKL")));
      if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(), 
			     10.0*cent->GetCentralityPercentile("TKL"), lrint(10.0*cent->GetCentralityPercentile("TKL")));
    }
  }
  else if (fEstEventMult == kCentralityNPA) {
    // centrality between 0 (central) and 1 (very peripheral)

    if (cent) {
      if (cent->GetCentralityPercentile("NPA") < 0.00001)
	hbtEvent->SetNormalizedMult(-1);
      else
	hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("NPA")));
      if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(), 
			     10.0*cent->GetCentralityPercentile("NPA"), lrint(10.0*cent->GetCentralityPercentile("NPA")));
    }
  }
  else if (fEstEventMult == kCentralityCND) {
    // centrality between 0 (central) and 1 (very peripheral)

    if (cent) {
      if (cent->GetCentralityPercentile("CND") < 0.00001)
	hbtEvent->SetNormalizedMult(-1);
      else
	hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("CND")));
      if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(), 
			     10.0*cent->GetCentralityPercentile("CND"), lrint(10.0*cent->GetCentralityPercentile("CND")));
    }
  }
  else if (fEstEventMult == kCentralityFMD) {
    // centrality between 0 (central) and 1 (very peripheral)

    if (cent) {
      if (cent->GetCentralityPercentile("FMD") < 0.00001)
	hbtEvent->SetNormalizedMult(-1);
      else
	hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("FMD")));
      if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(), 
			     10.0*cent->GetCentralityPercentile("FMD"), lrint(10.0*cent->GetCentralityPercentile("FMD")));
    }
  }



  if (tNormMultPos > tNormMultNeg)
    hbtEvent->SetZDCParticipants(tNormMultPos);
  else
    hbtEvent->SetZDCParticipants(tNormMultNeg);
  
  AliEventplane* ep = fEvent->GetEventplane();
  if (ep) {
    hbtEvent->SetEP(ep);
    hbtEvent->SetReactionPlaneAngle(ep->GetEventplane("Q"));
  }

  //V0
  if(fReadV0)
    {
      for (Int_t i = 0; i < fEvent->GetNumberOfV0s(); i++) {
	AliESDv0* esdv0 = fEvent->GetV0(i);
	if (!esdv0) continue;
	//if(esdv0->GetNDaughters()>2) continue;
	//if(esdv0->GetNProngs()>2) continue;
	if(esdv0->Charge()!=0) continue;
	AliESDtrack *trackPos = fEvent->GetTrack(esdv0->GetPindex());
	if(!trackPos) continue;
	AliESDtrack *trackNeg = fEvent->GetTrack(esdv0->GetNindex());
	if(!trackNeg) continue;
	if(trackPos->Charge()==trackNeg->Charge()) continue;
	AliFemtoV0* trackCopyV0 = new AliFemtoV0();
	CopyESDtoFemtoV0(esdv0, trackCopyV0, fEvent);
	hbtEvent->V0Collection()->push_back(trackCopyV0);
	//cout<<"Pushback v0 to v0collection"<<endl;
      }
    }


  if (Debug()>1) cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
}

//___________________
void AliFemtoEventReaderESDChain::SetESDSource(AliESDEvent *aESD)
{
  // The chain loads the ESD for us
  // You must provide the address where it can be found
  fEvent = aESD;
}
//___________________
// void AliFemtoEventReaderESDChain::SetESDfriendSource(AliESDfriend *aFriend)
// {
//   // We need the ESD tree to obtain
//   // information about the friend file location
//   fEventFriend = aFriend;
// }

//____________________________________________________________________
Float_t AliFemtoEventReaderESDChain::GetSigmaToVertex(double *impact, double *covar)
{
  // Calculates the number of sigma to the vertex.

  Float_t b[2];
  Float_t bRes[2];
  Float_t bCov[3];

  b[0] = impact[0];
  b[1] = impact[1];
  bCov[0] = covar[0];
  bCov[1] = covar[1];
  bCov[2] = covar[2];

  bRes[0] = TMath::Sqrt(bCov[0]);
  bRes[1] = TMath::Sqrt(bCov[2]);

  // -----------------------------------
  // How to get to a n-sigma cut?
  //
  // The accumulated statistics from 0 to d is
  //
  // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
  // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
  //
  // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
  // Can this be expressed in a different way?

  if (bRes[0] == 0 || bRes[1] ==0)
    return -1;

  Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));

  // stupid rounding problem screws up everything:
  // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
  if (TMath::Exp(-d * d / 2) < 1e-10)
    return 1000;

  d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
  return d;
}

void AliFemtoEventReaderESDChain::CopyESDtoFemtoV0(AliESDv0 *tESDv0, AliFemtoV0 *tFemtoV0, AliESDEvent *tESDevent)
{
  double fPrimaryVtxPosition[3];
  tESDevent->GetPrimaryVertex()->GetXYZ(fPrimaryVtxPosition);
  tFemtoV0->SetdcaV0ToPrimVertex(tESDv0->GetD(fPrimaryVtxPosition[0],fPrimaryVtxPosition[1],fPrimaryVtxPosition[2]));

  //tFemtoV0->SetdecayLengthV0(tESDv0->DecayLengthV0()); //wrocic do tego
  //tFemtoV0->SetdecayVertexV0X(tESDv0->DecayVertexV0X());
  //tFemtoV0->SetdecayVertexV0Y(tESDv0->DecayVertexV0Y());
  //tFemtoV0->SetdecayVertexV0Z(tESDv0->DecayVertexV0Z()); //nie ma w AliESDv0
  //AliFemtoThreeVector decayvertex(tESDv0->DecayVertexV0X(),tESDv0->DecayVertexV0Y(),tESDv0->DecayVertexV0Z());
  //tFemtoV0->SetdecayVertexV0(decayvertex);
  tFemtoV0->SetdcaV0Daughters(tESDv0->GetDcaV0Daughters());
  tFemtoV0->SetmomV0X(tESDv0->Px());
  tFemtoV0->SetmomV0Y(tESDv0->Py());
  tFemtoV0->SetmomV0Z(tESDv0->Pz());
  AliFemtoThreeVector momv0(tESDv0->Px(),tESDv0->Py(),tESDv0->Pz());
  tFemtoV0->SetmomV0(momv0);
  tFemtoV0->SetalphaV0(tESDv0->AlphaV0());
  tFemtoV0->SetptArmV0(tESDv0->PtArmV0());
  //tFemtoV0->SeteLambda(tESDv0->ELambda());
  //tFemtoV0->SeteK0Short(tESDv0->EK0Short());
  //tFemtoV0->SetePosProton(tESDv0->EPosProton());
  //tFemtoV0->SeteNegProton(tESDv0->ENegProton());
  tFemtoV0->SetmassLambda(tESDv0->GetEffMass(4,2));
  tFemtoV0->SetmassAntiLambda(tESDv0->GetEffMass(2,4));
  tFemtoV0->SetmassK0Short(tESDv0->GetEffMass(2,2));
  //tFemtoV0->SetrapLambda(tESDv0->RapLambda());
  //tFemtoV0->SetrapK0Short(tESDv0->RapK0Short());
  tFemtoV0->SetptV0(tESDv0->Pt());
  tFemtoV0->SetptotV0(tESDv0->P());
  tFemtoV0->SetEtaV0(tESDv0->Eta());
  tFemtoV0->SetPhiV0(tESDv0->Phi());
  tFemtoV0->SetCosPointingAngle(tESDv0->GetV0CosineOfPointingAngle(fPrimaryVtxPosition[0],fPrimaryVtxPosition[1], fPrimaryVtxPosition[2]));
  tFemtoV0->SetYV0(tESDv0->Y());

  AliESDtrack *trackpos = tESDevent->GetTrack(tESDv0->GetPindex()); //AliAODTrack *trackpos = (AliAODTrack*)tESDv0->GetDaughter(0);
  AliESDtrack *trackneg = tESDevent->GetTrack(tESDv0->GetNindex()); //AliAODTrack *trackneg = (AliAODTrack*)tESDv0->GetDaughter(1);

  if(trackpos && trackneg)
    {
      tFemtoV0->SetdcaPosToPrimVertex(TMath::Abs(trackpos->GetD(fPrimaryVtxPosition[0],fPrimaryVtxPosition[1],tESDevent->GetMagneticField())));
      tFemtoV0->SetdcaNegToPrimVertex(TMath::Abs(trackneg->GetD(fPrimaryVtxPosition[0],fPrimaryVtxPosition[1],tESDevent->GetMagneticField())));
      double MomPos[3];
      trackpos->PxPyPz(MomPos);
      tFemtoV0->SetmomPosX(MomPos[0]);
      tFemtoV0->SetmomPosY(MomPos[1]);
      tFemtoV0->SetmomPosZ(MomPos[2]);
      AliFemtoThreeVector mompos(MomPos[0],MomPos[1],MomPos[2]);
      tFemtoV0->SetmomPos(mompos);

      double MomNeg[3];
      trackneg->PxPyPz(MomNeg);
      tFemtoV0->SetmomNegX(MomNeg[0]);
      tFemtoV0->SetmomNegY(MomNeg[1]);
      tFemtoV0->SetmomNegZ(MomNeg[2]);
      AliFemtoThreeVector momneg(MomNeg[0],MomNeg[1],MomNeg[2]);
      tFemtoV0->SetmomNeg(momneg);

      tFemtoV0->SetptPos(trackpos->Pt());
      tFemtoV0->SetptotPos(trackpos->P());
      tFemtoV0->SetptNeg(trackneg->Pt());
      tFemtoV0->SetptotNeg(trackneg->P());
      
      tFemtoV0->SetidNeg(trackneg->GetID());
      //cout<<"tESDv0->GetNegID(): "<<tESDv0->GetNegID()<<endl;
      //cout<<"tFemtoV0->IdNeg(): "<<tFemtoV0->IdNeg()<<endl;
      tFemtoV0->SetidPos(trackpos->GetID());
      
      tFemtoV0->SetEtaPos(trackpos->Eta());
      tFemtoV0->SetEtaNeg(trackneg->Eta());

      tFemtoV0->SetEtaPos(trackpos->Eta()); //tESDv0->PseudoRapPos()
      tFemtoV0->SetEtaNeg(trackneg->Eta()); //tESDv0->PseudoRapNeg()
      tFemtoV0->SetTPCNclsPos(trackpos->GetTPCNcls());
      tFemtoV0->SetTPCNclsNeg(trackneg->GetTPCNcls());
      tFemtoV0->SetTPCclustersPos(trackpos->GetTPCClusterMap());
      tFemtoV0->SetTPCclustersNeg(trackneg->GetTPCClusterMap());
      tFemtoV0->SetTPCsharingPos(trackpos->GetTPCSharedMap());
      tFemtoV0->SetTPCsharingNeg(trackneg->GetTPCSharedMap());
      tFemtoV0->SetNdofPos(trackpos->GetTPCchi2()/trackpos->GetTPCNcls());
      tFemtoV0->SetNdofNeg(trackneg->GetTPCchi2()/trackneg->GetTPCNcls());
      tFemtoV0->SetStatusPos(trackpos->GetStatus());
      tFemtoV0->SetStatusNeg(trackneg->GetStatus());

      float bfield = 5*fMagFieldSign;
      float globalPositionsAtRadiiPos[9][3];
      GetGlobalPositionAtGlobalRadiiThroughTPC(trackpos,bfield,globalPositionsAtRadiiPos);
      double tpcEntrancePos[3]={globalPositionsAtRadiiPos[0][0],globalPositionsAtRadiiPos[0][1],globalPositionsAtRadiiPos[0][2]};
      double tpcExitPos[3]={globalPositionsAtRadiiPos[7][0],globalPositionsAtRadiiPos[7][1],globalPositionsAtRadiiPos[7][2]};

      float globalPositionsAtRadiiNeg[9][3];
      GetGlobalPositionAtGlobalRadiiThroughTPC(trackneg,bfield,globalPositionsAtRadiiNeg);
      double tpcEntranceNeg[3]={globalPositionsAtRadiiNeg[0][0],globalPositionsAtRadiiNeg[0][1],globalPositionsAtRadiiNeg[0][2]};
      double tpcExitNeg[3]={globalPositionsAtRadiiNeg[7][0],globalPositionsAtRadiiNeg[7][1],globalPositionsAtRadiiNeg[7][2]};

      AliFemtoThreeVector tmpVec;
      tmpVec.SetX(tpcEntrancePos[0]); tmpVec.SetX(tpcEntrancePos[1]); tmpVec.SetX(tpcEntrancePos[2]);
      tFemtoV0->SetNominalTpcEntrancePointPos(tmpVec);

      tmpVec.SetX(tpcExitPos[0]); tmpVec.SetX(tpcExitPos[1]); tmpVec.SetX(tpcExitPos[2]);
      tFemtoV0->SetNominalTpcExitPointPos(tmpVec);

      tmpVec.SetX(tpcEntranceNeg[0]); tmpVec.SetX(tpcEntranceNeg[1]); tmpVec.SetX(tpcEntranceNeg[2]);
      tFemtoV0->SetNominalTpcEntrancePointNeg(tmpVec);

      tmpVec.SetX(tpcExitNeg[0]); tmpVec.SetX(tpcExitNeg[1]); tmpVec.SetX(tpcExitNeg[2]);
      tFemtoV0->SetNominalTpcExitPointNeg(tmpVec);

      AliFemtoThreeVector vecTpcPos[9];
      AliFemtoThreeVector vecTpcNeg[9];
      for(int i=0;i<9;i++)
	{
	  vecTpcPos[i].SetX(globalPositionsAtRadiiPos[i][0]); vecTpcPos[i].SetY(globalPositionsAtRadiiPos[i][1]); vecTpcPos[i].SetZ(globalPositionsAtRadiiPos[i][2]);
	  vecTpcNeg[i].SetX(globalPositionsAtRadiiNeg[i][0]); vecTpcNeg[i].SetY(globalPositionsAtRadiiNeg[i][1]); vecTpcNeg[i].SetZ(globalPositionsAtRadiiNeg[i][2]);
	}
      tFemtoV0->SetNominalTpcPointPos(vecTpcPos);
      tFemtoV0->SetNominalTpcPointNeg(vecTpcNeg);

      tFemtoV0->SetTPCMomentumPos(trackpos->GetTPCInnerParam()->P()); //trackpos->GetTPCmomentum();
      tFemtoV0->SetTPCMomentumNeg(trackneg->GetTPCInnerParam()->P()); //trackneg->GetTPCmomentum();

      tFemtoV0->SetdedxPos(trackpos->GetTPCsignal());
      tFemtoV0->SetdedxNeg(trackneg->GetTPCsignal());


      if (fESDpid) {
	tFemtoV0->SetPosNSigmaTPCK(fESDpid->NumberOfSigmasTPC(trackpos,AliPID::kKaon));
	tFemtoV0->SetNegNSigmaTPCK(fESDpid->NumberOfSigmasTPC(trackneg,AliPID::kKaon));
	tFemtoV0->SetPosNSigmaTPCP(fESDpid->NumberOfSigmasTPC(trackpos,AliPID::kProton));
	tFemtoV0->SetNegNSigmaTPCP(fESDpid->NumberOfSigmasTPC(trackneg,AliPID::kProton));
	tFemtoV0->SetPosNSigmaTPCPi(fESDpid->NumberOfSigmasTPC(trackpos,AliPID::kPion));
	tFemtoV0->SetNegNSigmaTPCPi(fESDpid->NumberOfSigmasTPC(trackneg,AliPID::kPion));
      }
      else {
	tFemtoV0->SetPosNSigmaTPCK(-1000);
	tFemtoV0->SetNegNSigmaTPCK(-1000);
	tFemtoV0->SetPosNSigmaTPCP(-1000);
	tFemtoV0->SetNegNSigmaTPCP(-1000);
	tFemtoV0->SetPosNSigmaTPCPi(-1000);
	tFemtoV0->SetNegNSigmaTPCPi(-1000);
      }

      if(// (tFemtoV0->StatusPos()&AliESDtrack::kTOFpid)==0 ||
         (tFemtoV0->StatusPos()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTOFout)==0)
	{
	  if(// (tFemtoV0->StatusNeg()&AliESDtrack::kTOFpid)==0 ||
       (tFemtoV0->StatusNeg()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTOFout)==0)
	    {
	      tFemtoV0->SetPosNSigmaTOFK(-1000);
	      tFemtoV0->SetNegNSigmaTOFK(-1000);
	      tFemtoV0->SetPosNSigmaTOFP(-1000);
	      tFemtoV0->SetNegNSigmaTOFP(-1000);
	      tFemtoV0->SetPosNSigmaTOFPi(-1000);
	      tFemtoV0->SetNegNSigmaTOFPi(-1000);
	    }
	}
      else
	{
	  if (fESDpid) {
	    tFemtoV0->SetPosNSigmaTOFK(fESDpid->NumberOfSigmasTOF(trackpos,AliPID::kKaon));
	    tFemtoV0->SetNegNSigmaTOFK(fESDpid->NumberOfSigmasTOF(trackneg,AliPID::kKaon));
	    tFemtoV0->SetPosNSigmaTOFP(fESDpid->NumberOfSigmasTOF(trackpos,AliPID::kProton));
	    tFemtoV0->SetNegNSigmaTOFP(fESDpid->NumberOfSigmasTOF(trackneg,AliPID::kProton));
	    tFemtoV0->SetPosNSigmaTOFPi(fESDpid->NumberOfSigmasTOF(trackpos,AliPID::kPion));
	    tFemtoV0->SetNegNSigmaTOFPi(fESDpid->NumberOfSigmasTOF(trackneg,AliPID::kPion));
	  }
	  else {
	    tFemtoV0->SetPosNSigmaTOFK(-1000);
	    tFemtoV0->SetNegNSigmaTOFK(-1000);
	    tFemtoV0->SetPosNSigmaTOFP(-1000);
	    tFemtoV0->SetNegNSigmaTOFP(-1000);
	    tFemtoV0->SetPosNSigmaTOFPi(-1000);
	    tFemtoV0->SetNegNSigmaTOFPi(-1000);
	  }
	}
    }
  else
    {
      tFemtoV0->SetStatusPos(999);
      tFemtoV0->SetStatusNeg(999);
    }
  tFemtoV0->SetOnFlyStatusV0(tESDv0->GetOnFlyStatus());
}

void AliFemtoEventReaderESDChain::SetReadTrackType(ReadTrackType aType)
{
  fTrackType = aType;
}

//trigger
void AliFemtoEventReaderESDChain::SetEventTrigger(UInt_t eventtrig)
{
  fEventTrig = eventtrig;
}

//V0 reading
void AliFemtoEventReaderESDChain::SetReadV0(bool a)
{
  fReadV0 = a;
}

void AliFemtoEventReaderESDChain::SetMagneticFieldSign(int s)
{
  if(s>0)
    fMagFieldSign = 1;
  else if(s<0)
    fMagFieldSign = -1;
  else
    fMagFieldSign = 0;
}

void AliFemtoEventReaderESDChain::GetGlobalPositionAtGlobalRadiiThroughTPC(AliESDtrack *track, Float_t bfield, Float_t globalPositionsAtRadii[9][3])
{
  // Gets the global position of the track at nine different radii in the TPC
  // track is the track you want to propagate
  // bfield is the magnetic field of your event
  // globalPositionsAtRadii is the array of global positions in the radii and xyz

  // Initialize the array to something indicating there was no propagation
  for(Int_t i=0;i<9;i++){
    for(Int_t j=0;j<3;j++){
      globalPositionsAtRadii[i][j]=-9999.;
    }
  }

  // Make a copy of the track to not change parameters of the track
  AliExternalTrackParam etp; etp.CopyFromVTrack(track);
  //printf("\nAfter CopyFromVTrack\n");
  //etp.Print();

  // The global position of the the track
  Double_t xyz[3]={-9999.,-9999.,-9999.};

  // Counter for which radius we want
  Int_t iR=0;
  // The radii at which we get the global positions
  // IROC (OROC) from 84.1 cm to 132.1 cm (134.6 cm to 246.6 cm)
  Float_t Rwanted[9]={85.,105.,125.,145.,165.,185.,205.,225.,245.};
  // The global radius we are at
  Float_t globalRadius=0;

  // Propagation is done in local x of the track
  for (Float_t x = etp.GetX();x<247.;x+=1.){ // GetX returns local coordinates
    // Starts at the tracks fX and goes outwards. x = 245 is the outer radial limit
    // of the TPC when the track is straight, i.e. has inifinite pt and doesn't get bent.
    // If the track's momentum is smaller than infinite, it will develop a y-component, which
    // adds to the global radius

    // Stop if the propagation was not succesful. This can happen for low pt tracks
    // that don't reach outer radii
    if(!etp.PropagateTo(x,bfield))break;
    etp.GetXYZ(xyz); // GetXYZ returns global coordinates
    globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii

    // Roughly reached the radius we want
    if(globalRadius > Rwanted[iR]){

      // Bigger loop has bad precision, we're nearly one centimeter too far, go back in small steps.
      while (globalRadius>Rwanted[iR]){
	x-=.1;
	//      printf("propagating to x %5.2f\n",x);
	if(!etp.PropagateTo(x,bfield))break;
	etp.GetXYZ(xyz); // GetXYZ returns global coordinates
	globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
      }
      //printf("At Radius:%05.2f (local x %5.2f). Setting position to x %4.1f y %4.1f z %4.1f\n",globalRadius,x,xyz[0],xyz[1],xyz[2]);
      globalPositionsAtRadii[iR][0]=xyz[0];
      globalPositionsAtRadii[iR][1]=xyz[1];
      globalPositionsAtRadii[iR][2]=xyz[2];
      // Indicate we want the next radius
      iR+=1;
    }
    if(iR>=8){
      // TPC edge reached
      return;
    }
  }
}


void AliFemtoEventReaderESDChain::SetpA2013(Bool_t pA2013)
{
  fpA2013 = pA2013;
}

void AliFemtoEventReaderESDChain::SetUseMVPlpSelection(Bool_t mvplp)
{
  fMVPlp = mvplp;
}

void AliFemtoEventReaderESDChain::SetIsPileUpEvent(Bool_t ispileup)
{
  fisPileUp = ispileup;
}