#include "TROOT.h"
#include "TSystem.h"
#include "TChain.h"
#include "TFile.h"
#include "TString.h"
#include <iostream>

const Char_t* defaultInFile = "sample.list";
const Char_t *defaultOutFile = "oFemtoDst_mesons_pau200.root";

class StFemtoDstInclusiveSelector;
StFemtoDstInclusiveSelector *mFemtoDstInclSelector = NULL;

//_________________
void FemtoDstMaker_mesons_pau200(const Char_t* inFileList = defaultInFile,
				 const Char_t* outFileName = defaultOutFile)
{

  std::cout << "**************************" << std::endl
	    << "*       FemtoDstMaker    *" << std::endl
	    << "*           Start        *" << std::endl
	    << "**************************" << std::endl << endl;

  //Constants, cut and initial parameters
  Float_t cVtxZ[2] = {-70., 70.};
  Float_t cVtxR[2] = {0., 2.};
  Float_t cVtxXShift = 0.;
  Float_t cVtxYShift = 0.;
  Float_t cVtxVpdVzDiff[2] = {-7., 7.};
  Float_t cParticleMom[2] = {0.14, 1.6};
  Float_t cTrackDcaGlobal[2] = {-0.01, 5.};
  Int_t   cTrackNHits[2] = {15, 50};
  Int_t   cTrackNHitsFit[2] = {5, 50};
  Float_t cTrackEta[2] = {-1., 1.};
  Int_t   cTrackFlag[2] = {0, 1000};
  Bool_t  cZdcRefmultCorrection = true;   //Should be false for low energies

  //Inclusive cuts
  //Pions
  Bool_t  cSelectPions = true;
  Bool_t  cSelectTpcPions = true;
  Float_t cPionTpcMom[2] = {0.15, 0.75};
  Float_t cPionTpcNSigmaElectron[2] = {-1., 1.};
  Float_t cPionTpcNSigmaPion[2] = {-2., 2.};
  Float_t cPionTpcNSigmaKaon[2] = {-1., 1.};
  Float_t cPionTpcNSigmaProton[2] = {-1., 1.};
  Bool_t  cSelectTofPions = true;
  Float_t cPionTofMom[2] = {0.2, 1.55};
  Float_t cPionMassSqr[2] = {-0.05, 0.09};

  //Kaons
  Bool_t  cSelectKaons = true;
  Bool_t  cSelectTpcKaons = true;
  Float_t cKaonTpcMom[2] = {0.15, 0.75};
  Float_t cKaonTpcNSigmaElectron[2] = {-1., 1.};
  Float_t cKaonTpcNSigmaPion[2] = {-1., 1.};
  Float_t cKaonTpcNSigmaKaon[2] = {-2., 2.};
  Float_t cKaonTpcNSigmaProton[2] = {-1., 1.};
  Bool_t  cSelectTofKaons = true;
  Float_t cKaonTofMom[2] = {0.2, 1.55};
  Float_t cKaonMassSqr[2] = {0.16, 0.35};

  //Protons
  Bool_t  cSelectProtons = false;
  Bool_t  cSelectTpcProtons = false;
  Float_t cProtonTpcMom[2] = {0.15, 1.05};
  Float_t cProtonTpcNSigmaElectron[2] = {-2., 2.};
  Float_t cProtonTpcNSigmaPion[2] = {-2., 2.};
  Float_t cProtonTpcNSigmaKaon[2] = {-2., 2.};
  Float_t cProtonTpcNSigmaProton[2] = {-3., 3.};
  Bool_t  cSelectTofProtons = false;
  Float_t cProtonTofMom[2] = {0.15, 1.7};
  Float_t cProtonMassSqr[2] = {0.75, 1.05};

  //Load libraries
  std::cout << "Loading libraries..." << std::endl;
  gROOT->LoadMacro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
  loadSharedLibraries();
  gSystem->Load("libMinuit");
  gSystem->Load("StChain");
  gSystem->Load("StRefMultCorr");
  gSystem->Load("StFlowMaker");
  gSystem->Load("StHbtMaker");
  gSystem->Load("StFemtoDstMaker");
  gSystem->Load("libStFemtoDstMaker.so");
  std::cout << "Libraries have been successfully loaded" << std::endl;

  //Create chain
  StChain *mChain = new StChain("StChain");
  mChain->SetDebug(0);
  StMuDebug::setLevel(0);

  //Create MuDst reader
  int nMaxFileToRead = 1e9;
  std::cout << Form("Try to read file: %s", inFileList) << std::endl;
  StMuDstMaker *muDstMaker = 
    new StMuDstMaker(0, 0, "", inFileList, "MuDst", nMaxFileToRead);
  muDstMaker->SetStatus("*",0);
  muDstMaker->SetStatus("Event*",1);
  muDstMaker->SetStatus("MuEvent*",1);
  muDstMaker->SetStatus("PrimaryVertices*",1);
  muDstMaker->SetStatus("PrimaryTracks*",1);
  muDstMaker->SetStatus("GlobalTracks*",1);
  muDstMaker->SetStatus("BTofHeader*",1);
  TChain *mCurrentTree = muDstMaker->chain();
  Long64_t mNEvents = mCurrentTree->GetEntries();
  std::cout << Form("Number of events in TTree: %d", mNEvents) << std::endl;


  //StFemtoDstMaker initialization
  mFemtoDstInclSelector = new StFemtoDstInclusiveSelector(muDstMaker,
							  outFileName);
    
  //
  // Trigger selection for pAu 200
  //          P16id
  //                                           // trg             -- number of events
  mFemtoDstInclSelector->AddTriggerId(500003); // VPDMB-5-trgonly -- 24901523
  mFemtoDstInclSelector->AddTriggerId(500004); // VPDMB-novtx     -- 23532384
  mFemtoDstInclSelector->AddTriggerId(500904); // VPDMB-30        -- 58381649
                                               // total events    == 106815556
  
  mFemtoDstInclSelector->SetVtxZCut(cVtxZ[0], cVtxZ[1]);
  mFemtoDstInclSelector->SetVtxRCut(cVtxR[0], cVtxR[1]);
  mFemtoDstInclSelector->SetVtxShift(cVtxXShift, cVtxYShift);
  mFemtoDstInclSelector->SetVtxVpdVzDiffCut(cVtxVpdVzDiff[0], cVtxVpdVzDiff[1]);
  mFemtoDstInclSelector->SetParticleMomentum(cParticleMom[0], cParticleMom[1]);
  mFemtoDstInclSelector->SetTrackDcaGlobal(cTrackDcaGlobal[0], cTrackDcaGlobal[1]);
  mFemtoDstInclSelector->SetTrackNHits(cTrackNHits[0], cTrackNHits[1]);
  mFemtoDstInclSelector->SetTrackNHitsFit(cTrackNHitsFit[0], cTrackNHitsFit[1]);
  mFemtoDstInclSelector->SetTrackEta(cTrackEta[0], cTrackEta[1]);
  mFemtoDstInclSelector->SetTrackFlag(cTrackFlag[0], cTrackFlag[1]);
  mFemtoDstInclSelector->SetCollisionTypeAuAu(false);
  mFemtoDstInclSelector->SetAuAuZdcCoincidenceEnergy(cZdcRefmultCorrection);
  mFemtoDstInclSelector->SetMuDstMaker(muDstMaker);

  //Inclusive selection
  //Pions
  mFemtoDstInclSelector->SetSelectPions(cSelectPions);
  mFemtoDstInclSelector->SetSelectTpcPions(cSelectTpcPions);
  mFemtoDstInclSelector->SetPionTpcMom(cPionTpcMom[0], cPionTpcMom[1]);
  mFemtoDstInclSelector->SetPionTpcNSigmaElectron(cPionTpcNSigmaElectron[0],cPionTpcNSigmaElectron[1]);
  mFemtoDstInclSelector->SetPionTpcNSigmaPion(cPionTpcNSigmaPion[0],cPionTpcNSigmaPion[1]);
  mFemtoDstInclSelector->SetPionTpcNSigmaKaon(cPionTpcNSigmaKaon[0],cPionTpcNSigmaKaon[1]);
  mFemtoDstInclSelector->SetPionTpcNSigmaProton(cPionTpcNSigmaProton[0],cPionTpcNSigmaProton[1]);
  mFemtoDstInclSelector->SetSelectTofPions(cSelectTofPions);
  mFemtoDstInclSelector->SetPionTofMom(cPionTofMom[0],cPionTofMom[1]);
  mFemtoDstInclSelector->SetPionMassSqr(cPionMassSqr[0],cPionMassSqr[1]);

  //Kaons
  mFemtoDstInclSelector->SetSelectKaons(cSelectKaons);
  mFemtoDstInclSelector->SetSelectTpcKaons(cSelectTpcKaons);
  mFemtoDstInclSelector->SetKaonTpcMom(cKaonTpcMom[0], cKaonTpcMom[1]);
  mFemtoDstInclSelector->SetKaonTpcNSigmaElectron(cKaonTpcNSigmaElectron[0],cKaonTpcNSigmaElectron[1]);
  mFemtoDstInclSelector->SetKaonTpcNSigmaPion(cKaonTpcNSigmaPion[0],cKaonTpcNSigmaPion[1]);
  mFemtoDstInclSelector->SetKaonTpcNSigmaKaon(cKaonTpcNSigmaKaon[0],cKaonTpcNSigmaKaon[1]);
  mFemtoDstInclSelector->SetKaonTpcNSigmaProton(cKaonTpcNSigmaProton[0],cKaonTpcNSigmaProton[1]);
  mFemtoDstInclSelector->SetSelectTofKaons(cSelectTofKaons);
  mFemtoDstInclSelector->SetKaonTofMom(cKaonTofMom[0],cKaonTofMom[1]);
  mFemtoDstInclSelector->SetKaonMassSqr(cKaonMassSqr[0],cKaonMassSqr[1]);

  //Protons
  mFemtoDstInclSelector->SetSelectProtons(cSelectProtons);
  mFemtoDstInclSelector->SetSelectTpcProtons(cSelectTpcProtons);
  mFemtoDstInclSelector->SetProtonTpcMom(cProtonTpcMom[0], cProtonTpcMom[1]);
  mFemtoDstInclSelector->SetProtonTpcNSigmaElectron(cProtonTpcNSigmaElectron[0],cProtonTpcNSigmaElectron[1]);
  mFemtoDstInclSelector->SetProtonTpcNSigmaPion(cProtonTpcNSigmaPion[0],cProtonTpcNSigmaPion[1]);
  mFemtoDstInclSelector->SetProtonTpcNSigmaKaon(cProtonTpcNSigmaKaon[0],cProtonTpcNSigmaKaon[1]);
  mFemtoDstInclSelector->SetProtonTpcNSigmaProton(cProtonTpcNSigmaProton[0],cProtonTpcNSigmaProton[1]);
  mFemtoDstInclSelector->SetSelectTofProtons(cSelectTofProtons);
  mFemtoDstInclSelector->SetProtonTofMom(cProtonTofMom[0],cProtonTofMom[1]);
  mFemtoDstInclSelector->SetProtonMassSqr(cProtonMassSqr[0],cProtonMassSqr[1]);

  //Chain initialization and loop over events
  mChain->Init();
  Int_t iReturn = 0;
  Int_t mNEventsProcessed = 0;
  Float_t mPercentCounter = 0.01;
  Float_t mProgress = 0.;

  unsigned int mSubCounter = 1;              // time
  time_t mStartTime, mStopTime, mDiffTime;   // time
  float mFrac;                               // time
  mStartTime = time(0);                      // time

  for(Int_t iEvent=0; iEvent<mNEvents; iEvent++) {

    mProgress = (Float_t)iEvent / (Float_t)mNEvents;
    mNEventsProcessed++;
    
    if(mProgress >= mPercentCounter) {
      mPercentCounter += mPercentCounter;
      mStopTime = time(0);
      mDiffTime = difftime(mStopTime, mStartTime);
      mFrac = (float)mDiffTime*(float)(mNEvents - iEvent)/(float)iEvent;
      mSubCounter++;
      std::cout << Form("Processing progress: %4.2f%%. Time left: %.1f sec", mProgress*100., mFrac)
		<< std::endl;
    }

    mChain->Clear();
    iReturn = mChain->Make(iEvent);
    if(iReturn != 0) {
      std::cout << "Error has been occured. Event processing has been stopped"
		<< std::endl;
      break;
    }
  } //for(Int_t iEvent=0; iEvent<mNEvents; iEvent++)

  mChain->Finish();
  delete mFemtoDstInclSelector;
  delete mChain;

  std::cout << "**************************" << std::endl
	    << "*       FemtoDstMaker    *" << std::endl
	    << "*          Finish        *" << std::endl
	    << "**************************" << std::endl << endl;
}