//_____________________________________________________________________________
//    Class for semiforward UPC filter
//    Author: Jaroslav Adam
//
//    Fill structure of AliUPCEvent
//_____________________________________________________________________________

// c++ headers
#include <iostream>
#include <string.h>

// root headers
#include "TList.h"
#include "TH1I.h"
#include "TTree.h"
#include "TClonesArray.h"
#include "TParticle.h"
#include "TObjString.h"
#include "TFile.h"
#include "TArrayF.h"

// aliroot headers
#include "AliAnalysisManager.h"
#include "AliInputEventHandler.h"
#include "AliESDEvent.h"
#include "AliAODEvent.h"
#include "AliMCEvent.h"
#include "AliAODVZERO.h"
#include "AliAODZDC.h"
#include "AliAODVertex.h"
#include "AliESDVZERO.h"
#include "AliESDZDC.h"
#include "AliESDVertex.h"
#include "AliAODTrack.h"
#include "AliMultiplicity.h"
#include "AliESDtrack.h"
#include "AliESDMuonTrack.h"
#include "AliAODMCParticle.h"
#include "AliMCParticle.h"
#include "AliTriggerAnalysis.h"
#include "AliMuonTrackCuts.h"
#include "AliESDtrackCuts.h"
#include "AliAODPid.h"
#include "AliGenEventHeader.h"
#include "AliAODMCHeader.h"

// my headers
#include "AliAnalysisTaskUpcFilterSemiforward.h"
#include "AliUPCTrack.h"
#include "AliUPCMuonTrack.h"
#include "AliUPCEvent.h"

ClassImp(AliAnalysisTaskUpcFilterSemiforward);

// task for upc semiforward filter
// jaroslav.adam@cern.ch

//_____________________________________________________________________________
AliAnalysisTaskUpcFilterSemiforward::AliAnalysisTaskUpcFilterSemiforward(const char *name)
 :AliAnalysisTaskSE(name),
  fIsESD(0), fIsMC(0), fMuonCuts(0x0), fTriggerAna(0x0), fCutsList(0x0),
  fHistList(0x0), fCounter(0x0), fTriggerCounter(0x0),
  fUPCEvent(0x0), fUPCTree(0x0)
{

  // Constructor

  DefineOutput(1, TTree::Class());
  DefineOutput(2, TList::Class());

}//AliAnalysisTaskUpcFilterSemiforward

//_____________________________________________________________________________
AliAnalysisTaskUpcFilterSemiforward::~AliAnalysisTaskUpcFilterSemiforward()
{
  // destructor

  if(fHistList) {delete fHistList; fHistList = 0x0;}
  if(fCounter) {delete fCounter; fCounter = 0x0;}
  if(fTriggerCounter) {delete fTriggerCounter; fTriggerCounter = 0x0;}
  if(fUPCEvent) {delete fUPCEvent; fUPCEvent = 0x0;}
  if(fUPCTree) {delete fUPCTree; fUPCTree = 0x0;}
  if(fMuonCuts) {delete fMuonCuts; fMuonCuts = 0x0;}
  if(fTriggerAna) {delete fTriggerAna; fTriggerAna = 0x0;}
  if(fCutsList) {delete[] fCutsList; fCutsList=0x0;}

}//~AliAnalysisTaskUpcFilterSemiforward

//_____________________________________________________________________________
void AliAnalysisTaskUpcFilterSemiforward::UserCreateOutputObjects()
{
  //muon track cuts
  fMuonCuts = new AliMuonTrackCuts("StdMuonCuts", "StdMuonCuts");
  fMuonCuts->SetFilterMask ( AliMuonTrackCuts::kMuPdca );
  fMuonCuts->Print("mask");
  fMuonCuts->SetAllowDefaultParams(kTRUE);
  fMuonCuts->SetPassName("muon_pass2");

  //trigger analysis for SPD FO fired chips in MC
  fTriggerAna = new AliTriggerAnalysis();
  fTriggerAna->SetAnalyzeMC( fIsMC );

  //ESD track cuts
  if( fIsESD ) {
    fCutsList = new AliESDtrackCuts*[32];
    for(UInt_t i=0; i<32; i++) fCutsList[i] = 0x0;

    // bits 0 - 10 set in $ALICE_ROOT/ANALYSIS/ESDfilter/macros/AddTaskESDFilter.C

    // Cuts on primary tracks
    AliESDtrackCuts* esdTrackCutsL = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
    fCutsList[0] = esdTrackCutsL;

    // standard cuts with very loose DCA
    AliESDtrackCuts* esdTrackCutsH = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE); 
    esdTrackCutsH->SetMaxDCAToVertexXY(2.4);
    esdTrackCutsH->SetMaxDCAToVertexZ(3.2);
    esdTrackCutsH->SetDCAToVertex2D(kTRUE);
    fCutsList[4] = esdTrackCutsH;

    // standard cuts with tight DCA cut
    AliESDtrackCuts *esdTrackCutsH2 = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011();
    fCutsList[5] = esdTrackCutsH2;

    // bits 11 - 31 free

    // selection by Daniele, UPC meeting, 15 July 2014
    // https://indico.cern.ch/event/330483/contribution/3/material/slides/0.pdf

    AliESDtrackCuts *cuts11 = new AliESDtrackCuts();
    cuts11->SetMinNClustersTPC(70);
    cuts11->SetRequireTPCRefit(kTRUE);
    cuts11->SetRequireITSRefit(kTRUE);
    cuts11->SetMaxChi2PerClusterTPC(4);
    cuts11->SetAcceptKinkDaughters(kFALSE);
    cuts11->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
    cuts11->SetMaxDCAToVertexXYPtDep("0.0105+0.0350/pt^1.01");
    cuts11->SetMaxDCAToVertexZ(2);
    cuts11->SetDCAToVertex2D(kFALSE);
    cuts11->SetRequireSigmaToVertex(kFALSE);
    cuts11->SetMaxChi2PerClusterITS(5.);
    fCutsList[11] = cuts11;

    AliESDtrackCuts *cuts12 = new AliESDtrackCuts();
    cuts12->SetMinNClustersTPC(70);
    cuts12->SetRequireTPCRefit(kTRUE);
    cuts12->SetRequireITSRefit(kTRUE);
    cuts12->SetMaxChi2PerClusterTPC(4);
    cuts12->SetAcceptKinkDaughters(kFALSE);
    cuts12->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
    cuts12->SetMaxDCAToVertexXYPtDep("0.0105+0.0350/pt^1.01");
    cuts12->SetMaxDCAToVertexZ(2);
    cuts12->SetDCAToVertex2D(kFALSE);
    cuts12->SetRequireSigmaToVertex(kFALSE);
    fCutsList[12] = cuts12;

    // setting in Evgeny's trees
    AliESDtrackCuts *cuts13 = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011();
    cuts13->SetMaxChi2TPCConstrainedGlobal(1e+10);
    fCutsList[13] = cuts13;

  }

  //output histograms
  fHistList = new TList();
  fHistList->SetOwner();
  fCounter = new TH1I("fCounter", "fCounter", 30, 1, 31);
  fHistList->Add(fCounter);
  fTriggerCounter = new TH2I("fTriggerCounter", "fTriggerCounter", 3000, 195000, 198000, NTRG+1, 0, NTRG+1);
  fHistList->Add(fTriggerCounter);

  //output tree
  fUPCEvent  = new AliUPCEvent();
  fUPCEvent->SetIsMC( fIsMC );
  fUPCEvent->SetIsESD( fIsESD );

  TDirectory *pwd = gDirectory;
  OpenFile(1);
  fUPCTree = new TTree("fUPCTree", "fUPCTree");
  pwd->cd();
  fUPCTree->Branch("fUPCEvent", &fUPCEvent);

  PostData(1, fUPCTree);
  PostData(2, fHistList);

}//UserCreateOutputObjects

//_____________________________________________________________________________
void AliAnalysisTaskUpcFilterSemiforward::NotifyRun()
{

  //cout<<"AliAnalysisTaskUpcFilterSemiforward::NotifyRun()"<<endl;

  fMuonCuts->SetRun(fInputHandler); // use input handler from AliAnalysisTaskSE

  if( fIsESD ) { ((AliESDEvent*) InputEvent())->InitMagneticField(); }

}//NotifyRun

//_____________________________________________________________________________
void AliAnalysisTaskUpcFilterSemiforward::UserExec(Option_t *) 
{

  //cout<<"#################### Next event ##################"<<endl;

  // input event
  AliVEvent *vEvent = InputEvent();
  if(!vEvent) return;

  fUPCEvent->ClearEvent();

  fCounter->Fill( 1 ); // 1 = analyzed events

  // trigger
  TString trigger = vEvent->GetFiredTriggerClasses();
  Bool_t trgClasses[NTRG]; // array of fired trigger classes
  for(Int_t itrg=0; itrg<NTRG; itrg++) trgClasses[itrg] = kFALSE;

  trgClasses[1] = trigger.Contains("CMUP6-B"); // p-Pb FW
  trgClasses[2] = trigger.Contains("CMUP3-B"); // Pb-p FW
  trgClasses[3] = trigger.Contains("CMUP8-B"); // Pb-p FW

  trgClasses[4] = trigger.Contains("CMUP7-B"); // p-Pb SFW
  trgClasses[5] = trigger.Contains("CMUP5-B"); // Pb-p SFW
  trgClasses[6] = trigger.Contains("CMUP9-B"); // Pb-p SFW

  trgClasses[7] = trigger.Contains("CCUP7-B"); // CEN

  trgClasses[8] = trigger.Contains("CMUP7-ACE");
  trgClasses[9] = trigger.Contains("CMUP5-ACE");
  trgClasses[10]= trigger.Contains("CMUP9-ACE");

  Bool_t isTrg = kFALSE;
  for(Int_t itrg=1; itrg<NTRG; itrg++) {
    if(!trgClasses[itrg]) continue;
    //trigger at itrg is fired
    fUPCEvent->SetTriggerClass( itrg , kTRUE );
    fTriggerCounter->Fill( vEvent->GetRunNumber() , itrg );
    isTrg = kTRUE;
  }
  if(!isTrg && !fIsMC) {PostData(2, fHistList); return;}
  //event passed the trigger

  fCounter->Fill( 2 ); // 2 = events after trigger (ESD and AOD)

  // ESD / AOD specific tasks: MC, SPD FO, L0 inputs, SPD tracklets, tracks, ZDC tdc
  Bool_t stat;
  if( fIsESD ) stat = RunESD();
  if(!fIsESD ) stat = RunAOD();
  if( !stat && !fIsMC ) {PostData(2, fHistList); return;}
  //event passed ESD / AOD specific selection

  fCounter->Fill( 3 ); // 3 = events after ESD / AOD specific part

  // input data
  const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
  // reconstruction pass: -1 = unknown, 1 = pass1, 2 = pass2, counter indices: 9 (unknown), 11 (pass1), 12 (pass2)
  fUPCEvent->SetRecoPass( -1 );
  if( strstr(filnam,"/pass1/") ) {fUPCEvent->SetRecoPass( 1 ); fCounter->Fill( 11 );}
  if( strstr(filnam,"/pass2/") ) {fUPCEvent->SetRecoPass( 2 ); fCounter->Fill( 12 );}
  if( fUPCEvent->GetRecoPass() < 0 ) fCounter->Fill( 9 );
  fUPCEvent->SetInputFileName( filnam );
  fUPCEvent->SetEventNumber( ((TTree*) GetInputData(0))->GetTree()->GetReadEntry() );
  fUPCEvent->SetRunNumber( vEvent->GetRunNumber() );

  //VZERO
  AliVVZERO *dataVZERO = vEvent->GetVZEROData();
  if(!dataVZERO) {PostData(2, fHistList); return;}

  fUPCEvent->SetV0ADecision( dataVZERO->GetV0ADecision() );
  fUPCEvent->SetV0CDecision( dataVZERO->GetV0CDecision() );
  for (UInt_t iv=0; iv<32; iv++) {
    if( dataVZERO->BBTriggerV0C((Int_t)iv) ) fUPCEvent->SetBBtriggerV0Cmask(iv);
    if( dataVZERO->GetBBFlag((Int_t)iv) ) fUPCEvent->SetBBFlagV0Cmask(iv);
  }

  //ZDC
  AliVZDC *dataZDC = vEvent->GetZDCData();
  if(!dataZDC) {PostData(2, fHistList); return;}

  //energy in ZDC
  Double_t eZnc=0., eZpc=0., eZna=0., eZpa=0.;
  for(Int_t i=0; i<5; i++) {
    eZnc += dataZDC->GetZNCTowerEnergy()[i];
    eZpc += dataZDC->GetZPCTowerEnergy()[i];
    eZna += dataZDC->GetZNATowerEnergy()[i];
    eZpa += dataZDC->GetZPATowerEnergy()[i];
  }
  fUPCEvent->SetZNCEnergy( eZnc );
  fUPCEvent->SetZPCEnergy( eZpc );
  fUPCEvent->SetZNAEnergy( eZna );
  fUPCEvent->SetZPAEnergy( eZpa );

  //default primary vertex
  const AliVVertex *vtx = vEvent->GetPrimaryVertex();
  if(!vtx) {PostData(2, fHistList); return;}
  Double_t pos[3]; vtx->GetXYZ(pos);
  Double_t cov[6]; vtx->GetCovarianceMatrix(cov);
  fUPCEvent->SetPrimaryVertex(pos, vtx->GetChi2perNDF(), cov, vtx->GetNContributors());
  const char *vtxtitle = vtx->GetTitle();
  fUPCEvent->SetPrimaryVertexTitle( vtxtitle );

  fCounter->Fill( 4 ); // 4 = events written to the tree (ESD and AOD)

  fUPCTree ->Fill();
  PostData(1, fUPCTree);
  PostData(2, fHistList);

}//UserExec

//_____________________________________________________________________________
Bool_t AliAnalysisTaskUpcFilterSemiforward::RunAOD()
{
  //cout<<"#################### AOD event ##################"<<endl;

  //input AOD event
  AliAODEvent *aodEvent = (AliAODEvent*) InputEvent();
  if(!aodEvent) return kFALSE;

  fCounter->Fill( 22 ); // 22 = AOD analyzed events

  if(fIsMC) RunAODMC( (TClonesArray*) aodEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName()),
                      (AliAODMCHeader*) aodEvent->FindListObject(AliAODMCHeader::StdBranchName())
                    );

  //nominal interaction point
  static AliAODVertex *vtx0 = 0x0;
  if( !vtx0 ) {
    // make a vertex at coordinates 0, 0, 0
    vtx0 = new AliAODVertex();
    vtx0->SetX(0.); vtx0->SetY(0.); vtx0->SetZ(0.);
  }

  //array of tracks for DCAs calculation
  //static TClonesArray *trackArray = 0x0;
  //if( !trackArray ) trackArray = new TClonesArray("AliAODTrack", 3); // number of DCAs

  Double_t pxpypz[3];
  UChar_t maskMan;
  Double_t xyzDca[2], cov[3];
  Float_t b[2], covF[3];
  Int_t nmun=0, ncen=0;
  // AOD tracks loop
  for(Int_t itr=0; itr<aodEvent->GetNumberOfTracks(); itr++) {
    AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aodEvent->GetTrack(itr));
    if( !trk ) continue;

    //muon track
    if( trk->IsMuonTrack() ) {

      AliUPCMuonTrack *upcMuon = fUPCEvent->AddMuonTrack();
      upcMuon->SetPtEtaPhi( trk->Pt(), trk->Eta(), trk->Phi() );
      upcMuon->SetCharge( trk->Charge() );
      upcMuon->SetMatchTrigger( trk->GetMatchTrigger() );
      upcMuon->SetRAtAbsorberEnd( trk->GetRAtAbsorberEnd() );
      upcMuon->SetChi2perNDF( trk->Chi2perNDF() );
      upcMuon->SetDCA( trk->DCA() );
      upcMuon->SetPxDCA( fMuonCuts->IsSelected(trk) );

      nmun++;
    }
    //central track
    else {

      maskMan = 0;
      if( trk->HasPointOnITSLayer(0) )                 { maskMan |= 1 << 0; }
      if( trk->HasPointOnITSLayer(1) )                 { maskMan |= 1 << 1; }
      if( trk->GetStatus() & AliESDtrack::kITSrefit )  { maskMan |= 1 << 2; }
      if( trk->GetStatus() & AliESDtrack::kTPCrefit )  { maskMan |= 1 << 3; }

      //selection for at least one point in ITS and TPC refit
      //if( !(maskMan & (1 << 0)) && !(maskMan & (1 << 1)) ) continue;
      //if( !(maskMan & (1 << 3)) ) continue;

      //selection for its refit and tpc refit
      //if( !(maskMan & (1 << 2)) || !(maskMan & (1 << 3)) ) continue;

      //selection for tpc refit
      //if( !(maskMan & (1 << 3)) ) continue;

      //selection for at least one point in SPD and ITS refit and TPC refit
      //if( !(maskMan & (1 << 0)) && !(maskMan & (1 << 1)) ) continue;
      //if( !(maskMan & (1 << 2)) || !(maskMan & (1 << 3)) ) continue;

      //selection for at least one point in SPD and at least one TPC cluster
      if( !(maskMan & (1 << 0)) && !(maskMan & (1 << 1)) ) continue;
      if( trk->GetTPCNcls() == 0 ) continue;

      trk->GetPxPyPz(pxpypz);
      AliAODPid *apid = trk->GetDetPid();
      if(!apid) continue;

      AliUPCTrack *upcTrack = fUPCEvent->AddTrack();
      upcTrack->SetPxPyPz( pxpypz );
      upcTrack->SetMaskMan( maskMan );
      upcTrack->SetFilterMap( trk->GetFilterMap() );
      upcTrack->SetCharge( trk->Charge() );
      upcTrack->SetChi2perNDF( trk->Chi2perNDF() );
      upcTrack->SetTPCmomentum( apid->GetTPCmomentum() );
      upcTrack->SetTPCsignal( apid->GetTPCsignal() );
      upcTrack->SetTPCNcls( trk->GetTPCNcls() );
      upcTrack->SetTPCCrossedRows( (Float_t) trk->GetTPCNCrossedRows() );
      upcTrack->SetTPCNclsF( trk->GetTPCNclsF() );
      upcTrack->SetTPCNclsS( trk->GetTPCnclsS() );
      upcTrack->SetITSClusterMap( trk->GetITSClusterMap() );

      //array for DCA
      //trackArray->Clear("C");
      //new((*trackArray)[0]) AliAODTrack(*trk);
      //new((*trackArray)[1]) AliAODTrack(*trk);
      //new((*trackArray)[2]) AliAODTrack(*trk);

      //DCA to default primary vertex
      AliAODTrack *track1 = (AliAODTrack*) trk->Clone("track1");
      //AliAODTrack *track1 = (AliAODTrack*) trackArray->At( 0 );
      track1->PropagateToDCA(aodEvent->GetPrimaryVertex(), aodEvent->GetMagneticField(), 9999., xyzDca, cov);
      for(Int_t i=0; i<2; i++) {b[i] = (Float_t) xyzDca[i]; covF[i] = (Float_t) cov[i];}
      upcTrack->SetImpactParameters(b, covF);
      delete track1;

      //DCA to SPD vertex
      track1 = (AliAODTrack*) trk->Clone("track1");
      //track1 = (AliAODTrack*) trackArray->At( 1 );
      track1->PropagateToDCA(aodEvent->GetPrimaryVertexSPD(), aodEvent->GetMagneticField(), 9999., xyzDca, cov);
      for(Int_t i=0; i<2; i++) {b[i] = (Float_t) xyzDca[i]; covF[i] = (Float_t) cov[i];}
      upcTrack->SetImpactParametersSPD(b, covF);
      delete track1;

      //DCA to nominal interaction point
      track1 = (AliAODTrack*) trk->Clone("track1");
      //track1 = (AliAODTrack*) trackArray->At( 2 );
      track1->PropagateToDCA(vtx0, aodEvent->GetMagneticField(), 9999., xyzDca, cov);
      for(Int_t i=0; i<2; i++) {b[i] = (Float_t) xyzDca[i]; covF[i] = (Float_t) cov[i];}
      upcTrack->SetImpactParametersIP(b, covF);
      delete track1;

      ncen++;
    }
  }// AOD tracks loop

  //selection for at least one muon or central track
  if( nmun + ncen < 1 ) return kFALSE;

  //selection for at least one muon and at least one central track
  //if( nmun < 1 || ncen < 1 ) return kFALSE;


  // L0 trigger inputs
  fUPCEvent->SetL0inputs( aodEvent->GetHeader()->GetL0TriggerInputs() );

  // Tracklets
  fUPCEvent->SetNumberOfTracklets( aodEvent->GetTracklets()->GetNumberOfTracklets() );

  // AOD ZDC for TDC
  AliAODZDC *dataZDCAOD = aodEvent->GetZDCData();
  if(!dataZDCAOD) {PostData(2, fHistList); return kFALSE;}

  Bool_t znctdc = kFALSE, znatdc = kFALSE;
  if( dataZDCAOD->GetZNCTime() != 0. ) znctdc = kTRUE;
  if( dataZDCAOD->GetZNATime() != 0. ) znatdc = kTRUE;
  fUPCEvent->SetZNCtdc( znctdc );
  fUPCEvent->SetZNAtdc( znatdc );
  fUPCEvent->SetZPCtdc( kFALSE );
  fUPCEvent->SetZPAtdc( kFALSE );

  //SPD primary vertex in AOD
  AliAODVertex *vtx = aodEvent->GetPrimaryVertexSPD();
  if(!vtx) {PostData(2, fHistList); return kFALSE;}
  Double_t posVtx[3]; vtx->GetXYZ( posVtx );
  Double_t covVtx[6]; vtx->GetCovarianceMatrix( covVtx );
  fUPCEvent->SetPrimaryVertexSPD(posVtx, vtx->GetChi2perNDF(), covVtx, vtx->GetNContributors());
  const char *vtxtitle = vtx->GetTitle();
  fUPCEvent->SetPrimaryVertexSPDtitle( vtxtitle );

  return kTRUE;

}//RunAOD

//_____________________________________________________________________________
void AliAnalysisTaskUpcFilterSemiforward::RunAODMC(TClonesArray *arrayMC, AliAODMCHeader *headerMC)
{
  // run over AOD mc particles

  if( !arrayMC || !headerMC ) return;

  //generated vertex in AOD MC
  Double_t vtxD[3]; headerMC->GetVertex( vtxD );
  Float_t vtxF[3]; for(Int_t i=0; i<3; i++) vtxF[i] = (Float_t) vtxD[i];
  fUPCEvent->SetPrimaryVertexMC( vtxF );

  //loop over mc particles
  for(Int_t imc=0; imc<arrayMC->GetEntriesFast(); imc++) {
    AliAODMCParticle *aodmc = (AliAODMCParticle*) arrayMC->At(imc);
    if(!aodmc) continue;

    if(aodmc->GetMother() >= 0) continue;

    TParticle *part = fUPCEvent->AddMCParticle();
    part->SetMomentum(aodmc->Px(), aodmc->Py(), aodmc->Pz(), aodmc->E());
    part->SetProductionVertex(aodmc->Xv(), aodmc->Yv(), aodmc->Zv(), aodmc->T());
    part->SetFirstMother(aodmc->GetMother());
    part->SetLastDaughter(aodmc->GetNDaughters());
    part->SetPdgCode(aodmc->GetPdgCode());
    part->SetUniqueID(imc);
  }//loop over mc particles

}//RunAODMC

//_____________________________________________________________________________
Bool_t AliAnalysisTaskUpcFilterSemiforward::RunESD()
{
  //cout<<"#################### ESD event ##################"<<endl;

  // Input ESD event
  AliESDEvent *esdEvent = (AliESDEvent*) InputEvent();
  if(!esdEvent) return kFALSE;

  fCounter->Fill( 21 ); // 21 = ESD analyzed events

  if(fIsMC) {
    fUPCEvent->SetNSPDfiredInner( fTriggerAna->SPDFiredChips(esdEvent,1,kFALSE,1) );
    fUPCEvent->SetNSPDfiredOuter( fTriggerAna->SPDFiredChips(esdEvent,1,kFALSE,2) );

    RunESDMC();
  }

  static AliESDVertex *vtx0 = 0x0;
  if( !vtx0 ) {
    // make a vertex at coordinates 0, 0, 0
    vtx0 = new AliESDVertex();
    vtx0->SetXv(0.); vtx0->SetYv(0.); vtx0->SetZv(0.);
  }

  // ESD central tracks
  Double_t pxpypz[3];
  Float_t b[2]; Float_t cov[3];
  UChar_t maskMan;
  UInt_t filterMap;
  Int_t nmun=0, ncen=0;
  //ESD central tracks loop
  for(Int_t itr=0; itr<esdEvent->GetNumberOfTracks(); itr++) {
    AliESDtrack *eTrack = esdEvent->GetTrack(itr);
    if( !eTrack ) continue;

    // manual filter
    maskMan = 0;
    if( eTrack->HasPointOnITSLayer(0) )                 { maskMan |= 1 << 0; }
    if( eTrack->HasPointOnITSLayer(1) )                 { maskMan |= 1 << 1; }
    if( eTrack->GetStatus() & AliESDtrack::kITSrefit )  { maskMan |= 1 << 2; }
    if( eTrack->GetStatus() & AliESDtrack::kTPCrefit )  { maskMan |= 1 << 3; }
    if( eTrack->GetKinkIndex(0) > 0 )                   { maskMan |= 1 << 4; } // bit set if track is kink candidate

    //selection for at least one point in ITS and TPC refit
    //if( !(maskMan & (1 << 0)) && !(maskMan & (1 << 1)) ) continue;
    //if( !(maskMan & (1 << 3)) ) continue;

    //selection for its refit and tpc refit
    //if( !(maskMan & (1 << 2)) || !(maskMan & (1 << 3)) ) continue;

    //selection for tpc refit
    //if( !(maskMan & (1 << 3)) ) continue;

    //selection for at least one point in SPD and ITS refit and TPC refit
    //if( !(maskMan & (1 << 0)) && !(maskMan & (1 << 1)) ) continue;
    //if( !(maskMan & (1 << 2)) || !(maskMan & (1 << 3)) ) continue;

    //selection for at least one point in SPD and at least one TPC cluster
    if( !(maskMan & (1 << 0)) && !(maskMan & (1 << 1)) ) continue;
    if( eTrack->GetTPCNcls() == 0 ) continue;

    //central track accepted to write to the UPC event

    // various configurations of AliESDtrackCuts
    filterMap = 0;
    for(UInt_t i=0; i<32; i++) {
      if( !fCutsList[i] ) continue;

      if( fCutsList[i]->AcceptTrack(eTrack) ) { filterMap |= 1 << i; }
    }

    eTrack->GetPxPyPz(pxpypz);

    //fill UPC track
    AliUPCTrack *upcTrack = fUPCEvent->AddTrack();
    upcTrack->SetPxPyPz( pxpypz );
    upcTrack->SetMaskMan( maskMan );
    upcTrack->SetFilterMap( filterMap );
    upcTrack->SetCharge( eTrack->Charge() );
    upcTrack->SetChi2perNDF( eTrack->GetTPCchi2()/((Double_t) eTrack->GetTPCNcls()) ); // TPC chi2 used to fill this data member in esd
    upcTrack->SetTPCmomentum( eTrack->GetTPCmomentum() );
    upcTrack->SetTPCsignal( eTrack->GetTPCsignal() );
    upcTrack->SetTPCNcls( eTrack->GetTPCNcls() );
    upcTrack->SetTPCCrossedRows( eTrack->GetTPCCrossedRows() );
    upcTrack->SetTPCNclsF( eTrack->GetTPCNclsF() );
    upcTrack->SetTPCNclsS( eTrack->GetTPCnclsS() );
    upcTrack->SetITSchi2perNDF( eTrack->GetITSchi2()/((Double_t) eTrack->GetNcls(0)) );
    upcTrack->SetITSClusterMap( eTrack->GetITSClusterMap() );

    //DCA to default primary vertex
    eTrack->GetImpactParameters(b, cov);
    upcTrack->SetImpactParameters(b, cov);

    //DCA to SPD vertex
    AliESDtrack *track1 = (AliESDtrack*) eTrack->Clone("track1");
    track1->RelateToVertex( esdEvent->GetPrimaryVertexSPD(), esdEvent->GetMagneticField(), 9999. );
    track1->GetImpactParameters(b, cov);
    upcTrack->SetImpactParametersSPD(b, cov);
    delete track1;

    //DCA to nominal interaction point
    track1 = (AliESDtrack*) eTrack->Clone("track1");
    track1->RelateToVertex( vtx0, esdEvent->GetMagneticField(), 9999. );
    track1->GetImpactParameters(b, cov);
    upcTrack->SetImpactParametersIP(b, cov);
    delete track1;

    ncen++;
  } //ESD central tracks loop

  // ESD muon tracks
  //muon tracks loop
  for(Int_t itr=0; itr<esdEvent->GetNumberOfMuonTracks(); itr++) {
    AliESDMuonTrack *mTrack = esdEvent->GetMuonTrack(itr);
    if( !mTrack ) continue;

    AliUPCMuonTrack *upcMuon = fUPCEvent->AddMuonTrack();
    upcMuon->SetPtEtaPhi( mTrack->Pt(), mTrack->Eta(), mTrack->Phi() );
    upcMuon->SetCharge( mTrack->Charge() );
    upcMuon->SetMatchTrigger( mTrack->GetMatchTrigger() );
    upcMuon->SetRAtAbsorberEnd( mTrack->GetRAtAbsorberEnd() );
    upcMuon->SetChi2perNDF( mTrack->GetChi2()/((Double_t) mTrack->GetNDF()) );
    upcMuon->SetDCA( mTrack->GetDCA() );
    upcMuon->SetPxDCA( fMuonCuts->IsSelected(mTrack) );

    nmun++;
  } //muon tracks loop

  //selection for at least one muon or central track
  if( nmun + ncen < 1 ) return kFALSE;

  //selection for at least one muon and at least one central track
  //if( nmun < 1 || ncen < 1 ) return kFALSE;

  // L0 trigger inputs
  fUPCEvent->SetL0inputs( esdEvent->GetHeader()->GetL0TriggerInputs() );

  //Tracklets
  fUPCEvent->SetNumberOfTracklets( esdEvent->GetMultiplicity()->GetNumberOfTracklets() );

  // ESD ZDC for TDC
  AliESDZDC *dataZDCESD = esdEvent->GetESDZDC();
  if(!dataZDCESD) {PostData(2, fHistList); return kFALSE;}

  Bool_t znctdc = kFALSE, zpctdc = kFALSE, znatdc = kFALSE, zpatdc = kFALSE;
  for(Int_t iz=0;iz<4;iz++) {
    if( dataZDCESD->GetZDCTDCData(10,iz) ) znctdc = kTRUE;
    if( dataZDCESD->GetZDCTDCData(11,iz) ) zpctdc = kTRUE;
    if( dataZDCESD->GetZDCTDCData(12,iz) ) znatdc = kTRUE;
    if( dataZDCESD->GetZDCTDCData(13,iz) ) zpatdc = kTRUE;
  }
  fUPCEvent->SetZNCtdc( znctdc );
  fUPCEvent->SetZPCtdc( zpctdc );
  fUPCEvent->SetZNAtdc( znatdc );
  fUPCEvent->SetZPAtdc( zpatdc );

  //SPD primary vertex in ESD
  const AliESDVertex *vtx = esdEvent->GetPrimaryVertexSPD();
  if(!vtx) {PostData(2, fHistList); return kFALSE;}
  Double_t posVtx[3]; vtx->GetXYZ( posVtx );
  Double_t covVtx[6]; vtx->GetCovarianceMatrix( covVtx );
  fUPCEvent->SetPrimaryVertexSPD(posVtx, vtx->GetChi2perNDF(), covVtx, vtx->GetNContributors());
  const char *vtxtitle = vtx->GetTitle();
  fUPCEvent->SetPrimaryVertexSPDtitle( vtxtitle );

  return kTRUE;

}//RunESD

//_____________________________________________________________________________
void AliAnalysisTaskUpcFilterSemiforward::RunESDMC()
{
  // ESD MC particles

  AliMCEvent *mcEvent = MCEvent();
  if(!mcEvent) {cout<<"Error: no ESD MC found"<<endl; return;}

  //generated vertex
  TArrayF vtx(3);
  mcEvent->GenEventHeader()->PrimaryVertex(vtx);
  fUPCEvent->SetPrimaryVertexMC( vtx.GetArray() );

  //loop over mc particles
  for(Int_t imc=0; imc<mcEvent->GetNumberOfTracks(); imc++) {
    AliMCParticle *esdmc = (AliMCParticle*) mcEvent->GetTrack(imc);
    if(!esdmc) continue;

    if(esdmc->GetMother() >= 0) continue;

    TParticle *part = fUPCEvent->AddMCParticle();
    part->SetMomentum(esdmc->Px(), esdmc->Py(), esdmc->Pz(), esdmc->E());
    part->SetProductionVertex(esdmc->Xv(), esdmc->Yv(), esdmc->Zv(), 0.);
    part->SetFirstMother(esdmc->GetMother());
    part->SetLastDaughter(esdmc->GetLastDaughter()-esdmc->GetFirstDaughter()+1);
    part->SetPdgCode(esdmc->PdgCode());
    part->SetUniqueID(imc);
  }//loop over mc particles

}//RunESDMC

//_____________________________________________________________________________
void AliAnalysisTaskUpcFilterSemiforward::Terminate(Option_t *) 
{

  cout<<"Analysis complete."<<endl;
}//Terminate