// C++ headers
#include <iostream>
#include <string>
#include <utility>
#include <set>
#include <map>
#include <functional>

// ROOT headers
#include <TROOT.h>
#include <TFile.h>
#include <TChain.h>
#include <TTree.h>
#include <TSystem.h>
#include <TH1.h>
#include <TH2.h>
#include <TMath.h>
#include <TString.h>
#include <TStopwatch.h>

// FemtoDst headers
#include "StFemtoDstReader.h"
#include "StFemtoDst.h"
#include "StFemtoEvent.h"
#include "StFemtoTrack.h"

// AnalysisTree headers
#include "AnalysisTree/Configuration.hpp"
#include "AnalysisTree/DataHeader.hpp"
#include "AnalysisTree/EventHeader.hpp"
#include "AnalysisTree/Detector.hpp"
#include "AnalysisTree/Matching.hpp"

float get_beamP(float sqrtSnn, float m_target = 0.938, float m_beam = 0.938);
float GetCentFromCent9(int cent9);
float GetCentFromCent16(int cent16);
TVector3 GetBbcPositionEast(int imod);
TVector3 GetBbcPositionWest(int imod);
TVector3 GetZdcPositionXEast(int strip);
TVector3 GetZdcPositionYEast(int strip);
TVector3 GetZdcPositionXWest(int strip);
TVector3 GetZdcPositionYWest(int strip);

int main(int argc, char **argv)
{
  TString iFileName, oFileName;
  std::string system="";
  float sqrtSnn = -1.;
  bool isDataHeaderDefined = false;

  if (argc < 5)
  {
    std::cerr << "./MpdDst2AnalysisTree -i inputfile/inputlist -o outputfile [OPTIONAL: --sqrtSnn sqrtSnn --system colliding_system_name]" << std::endl;
    return 1;
  }
  for (int i = 1; i < argc; i++)
  {
    if (std::string(argv[i]) != "-i" &&
      std::string(argv[i]) != "-o" &&
      std::string(argv[i]) != "--sqrtSnn" &&
      std::string(argv[i]) != "--system")
    {
      std::cerr << "\n[ERROR]: Unknown parameter " << i << ": " << argv[i] << std::endl;
      return 1;
    }
    else
    {
      if (std::string(argv[i]) == "-i" && i != argc - 1)
      {
        iFileName = argv[++i];
        continue;
      }
      if (std::string(argv[i]) == "-i" && i == argc - 1)
      {
        std::cerr << "\n[ERROR]: Input file name was not specified " << std::endl;
        return 1;
      }
      if (std::string(argv[i]) == "-o" && i != argc - 1)
      {
        oFileName = argv[++i];
        continue;
       }
      if (std::string(argv[i]) == "-o" && i == argc - 1)
      {
        std::cerr << "\n[ERROR]: Output file name was not specified " << std::endl;
        return 1;
      }
      if (std::string(argv[i]) == "--system" && i != argc - 1)
      {
        system = std::string(argv[++i]);
        isDataHeaderDefined = true;
        continue;
       }
      if (std::string(argv[i]) == "--system" && i == argc - 1)
      {
        std::cerr << "\n[ERROR]: Output file name was not specified " << std::endl;
        return 1;
      }
      if (std::string(argv[i]) == "--sqrtSnn" && i != argc - 1)
      {
        sqrtSnn = std::stof(std::string(argv[++i]));
        isDataHeaderDefined = true;
        continue;
       }
      if (std::string(argv[i]) == "--sqrtSnn" && i == argc - 1)
      {
        std::cerr << "\n[ERROR]: Output file name was not specified " << std::endl;
        return 1;
      }
    }
  }

  const int NbbcModules = 32; // 0-15 - east, 16-31 - west
  const int NbbcModulesHalf = 16; 
  const int NzdcStipsX = 7;
  const int NzdcStipsY = 8;

  TStopwatch timer;
  timer.Start();

  StFemtoDstReader* femtoReader = new StFemtoDstReader(iFileName.Data());
  femtoReader->Init();

  // This is a way if you want to spead up IO
  std::cout << "Explicit read status for some branches" << std::endl;
  femtoReader->SetStatus("*",0);
  femtoReader->SetStatus("Event",1);
  femtoReader->SetStatus("Track",1);
  std::cout << "Status has been set" << std::endl;

  if( !femtoReader->chain() ) {
    std::cout << "No chain has been found." << std::endl;
  }
  Long64_t eventsInTree = femtoReader->tree()->GetEntries();
  std::cout << "eventsInTree: "  << eventsInTree << std::endl;
  Long64_t events2read = femtoReader->chain()->GetEntries();

  std::cout << "Number of events to read: " << events2read << std::endl;

  // Set up output dst
  TFile *outFile = new TFile(oFileName.Data(), "RECREATE");
  TTree *outTree = new TTree("aTree","AnalysisTree Dst at STAR");

  // Set up AnalysisTree data header
  AnalysisTree::DataHeader *dataHeader = new AnalysisTree::DataHeader;
  if (system != "") dataHeader->SetSystem(system);
  if (sqrtSnn > 0) dataHeader->SetBeamMomentum(get_beamP(sqrtSnn));
  auto &bbc_mod_pos = dataHeader->AddDetector();
  TVector3 modulePos;
  for (int imodule=0; imodule<NbbcModules; imodule++)
  {
    auto *module = bbc_mod_pos.AddChannel();
    if (imodule < NbbcModulesHalf) modulePos = GetBbcPositionEast(imodule);
    if (imodule >= NbbcModulesHalf) modulePos = GetBbcPositionWest(imodule-NbbcModulesHalf);
    module->SetPosition(modulePos);
  }

  // Set up AnalysisTree configureation
  AnalysisTree::Configuration *out_config = new AnalysisTree::Configuration;

  // Set up AnalysisTree configuration
  std::string str_reco_event_branch = "RecoEvent";
  std::string str_tpc_tracks_branch = "TpcTracks";
  std::string str_bbc_branch = "BbcModules";
  std::string str_zdc_horizontal_branch = "ZdcHitsHorizontal";
  std::string str_zdc_vertical_branch = "ZdcHitsVertical";

  AnalysisTree::BranchConfig reco_event_branch(str_reco_event_branch.c_str(), AnalysisTree::DetType::kEventHeader); 
  AnalysisTree::BranchConfig tpc_tracks_branch(str_tpc_tracks_branch.c_str(), AnalysisTree::DetType::kTrack);
  AnalysisTree::BranchConfig bbc_branch(str_bbc_branch.c_str(), AnalysisTree::DetType::kModule); 
  AnalysisTree::BranchConfig zdc_horizontal_branch(str_zdc_horizontal_branch.c_str(), AnalysisTree::DetType::kHit); 
  AnalysisTree::BranchConfig zdc_vertical_branch(str_zdc_vertical_branch.c_str(), AnalysisTree::DetType::kHit); 

  reco_event_branch.AddField<int>("runId");
  reco_event_branch.AddField<int>("refMult");
  reco_event_branch.AddField<int>("refMult2");
  reco_event_branch.AddField<float>("cent9");
  reco_event_branch.AddField<float>("cent16");
  reco_event_branch.AddField<int>("numberOfBTofHit");
  reco_event_branch.AddField<int>("numberOfTofMatched");
  reco_event_branch.AddField<float>("vpdVz");

  tpc_tracks_branch.AddField<int>("nHits");
  tpc_tracks_branch.AddField<int>("nHitsFit");
  tpc_tracks_branch.AddField<int>("nHitsPoss");
  tpc_tracks_branch.AddField<float>("nHits_Fit2Poss");
  tpc_tracks_branch.AddField<float>("nSigmaPion");
  tpc_tracks_branch.AddField<float>("nSigmaKaon");
  tpc_tracks_branch.AddField<float>("nSigmaProton");
  tpc_tracks_branch.AddField<float>("dEdx");
  tpc_tracks_branch.AddField<float>("tofMass2");
  tpc_tracks_branch.AddFields<float>({"dca_x","dca_y","dca_z"});
  tpc_tracks_branch.AddField<float>("dca");
  tpc_tracks_branch.AddField<bool>("isPrimary");
  tpc_tracks_branch.AddField<bool>("isTofTrack");

  auto hasher = std::hash<std::string>();

  // Initialize AnalysisTree Dst components
  out_config->AddBranchConfig(reco_event_branch);
  AnalysisTree::EventHeader *reco_event = new AnalysisTree::EventHeader( Short_t(hasher(reco_event_branch.GetName())) );
  out_config->AddBranchConfig(tpc_tracks_branch);
  AnalysisTree::TrackDetector *tpc_tracks = new AnalysisTree::TrackDetector( Short_t(hasher(tpc_tracks_branch.GetName())) ); 
  out_config->AddBranchConfig(bbc_branch);
  AnalysisTree::ModuleDetector *bbc_modules = new AnalysisTree::ModuleDetector( Short_t(hasher(bbc_branch.GetName())) );
  out_config->AddBranchConfig(zdc_horizontal_branch);
  AnalysisTree::HitDetector *zdc_horizontal_hits = new AnalysisTree::HitDetector( Short_t(hasher(zdc_horizontal_branch.GetName())) );
  out_config->AddBranchConfig(zdc_vertical_branch);
  AnalysisTree::HitDetector *zdc_vertical_hits = new AnalysisTree::HitDetector( Short_t(hasher(zdc_vertical_branch.GetName())) );

  reco_event->Init(reco_event_branch);

  // reco_event's additional field ids
  const int irunId = out_config->GetBranchConfig(str_reco_event_branch).GetFieldId("runId");
  const int irefMult = out_config->GetBranchConfig(str_reco_event_branch).GetFieldId("refMult");
  const int irefMult2 = out_config->GetBranchConfig(str_reco_event_branch).GetFieldId("refMult2");
  const int icent9 = out_config->GetBranchConfig(str_reco_event_branch).GetFieldId("cent9");
  const int icent16 = out_config->GetBranchConfig(str_reco_event_branch).GetFieldId("cent16");
  const int inumberOfBTofHit = out_config->GetBranchConfig(str_reco_event_branch).GetFieldId("numberOfBTofHit");
  const int inumberOfTofMatched = out_config->GetBranchConfig(str_reco_event_branch).GetFieldId("numberOfTofMatched");
  const int ivpdVz = out_config->GetBranchConfig(str_reco_event_branch).GetFieldId("vpdVz");

  // tpc_tracks' additional field ids
  const int inHits = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("nHits");
  const int inHitsFit = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("nHitsFit");
  const int inHitsPoss = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("nHitsPoss");
  const int inHits_Fit2Poss = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("nHits_Fit2Poss");
  const int inSigmaPion = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("nSigmaPion");
  const int inSigmaKaon = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("nSigmaKaon");
  const int inSigmaProton = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("nSigmaProton");
  const int idEdx = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("dEdx");
  const int itofMass2 = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("tofMass2");
  const int idca_x = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("dca_x");
  const int idca_y = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("dca_y");
  const int idca_z = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("dca_z");
  const int idca = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("dca");
  const int iisPrimary = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("isPrimary");
  const int iisTofTrack = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("isTofTrack");

  // Create branches in the output tree
  outTree->Branch(str_reco_event_branch.c_str(), "AnalysisTree::EventHeader", &reco_event, 32000, 99);
  outTree->Branch(str_tpc_tracks_branch.c_str(), "AnalysisTree::TrackDetector", &tpc_tracks, 256000, 99);
  outTree->Branch(str_bbc_branch.c_str(), "AnalysisTree::ModuleDetector",  &bbc_modules, 128000, 99);
  outTree->Branch(str_zdc_horizontal_branch.c_str(), "AnalysisTree::HitDetector",  &zdc_horizontal_hits, 128000, 99);
  outTree->Branch(str_zdc_vertical_branch.c_str(), "AnalysisTree::HitDetector",  &zdc_vertical_hits, 128000, 99);

  // Printout basic configuration info
  out_config->Print();

  // Starting event loop
  TVector3 pVtx;
  // Loop over events
  for(Long64_t iEvent=0; iEvent<events2read; iEvent++)
  {
    if (iEvent % 1000 == 0)
    {
      std::cout << "Working on event #[" << (iEvent+1)
                << "/" << events2read << "]" << std::endl;
    }

    Bool_t readEvent = femtoReader->readFemtoEvent(iEvent);
    if( !readEvent ) {
      std::cout << "Something went wrong! Nothing to analyze..." << std::endl;
      break;
    }

    // Retrieve femtoDst
    StFemtoDst *dst = femtoReader->femtoDst();

    // Retrieve event information
    StFemtoEvent *event = dst->event();
    if( !event ) {
      std::cout << "Something went wrong! Event is hiding from me..." << std::endl;
      break;
    }

    // Return primary vertex position
    pVtx = event->primaryVertex();

    // Reject vertices that are far from the central membrane along the beam
    if( TMath::Abs( pVtx.Z() ) > 200. ) continue;

    tpc_tracks->ClearChannels();
    bbc_modules->ClearChannels();
    zdc_horizontal_hits->ClearChannels();
    zdc_vertical_hits->ClearChannels();

    // Reading Reco Event
    reco_event->SetVertexPosition3(pVtx);
    reco_event->SetField(int(event->runId()), irunId);
    reco_event->SetField(int(event->refMult()), irefMult);
    reco_event->SetField(int(event->refMult2()), irefMult2);
    reco_event->SetField(float(GetCentFromCent9(event->cent9())), icent9);
    reco_event->SetField(float(GetCentFromCent16(event->cent16())), icent16);
    reco_event->SetField(int(event->numberOfBTofHit()), inumberOfBTofHit);
    reco_event->SetField(int(event->numberOfTofMatched()), inumberOfTofMatched);
    reco_event->SetField(float(event->vpdVz()), ivpdVz);

    // Track analysis
    Int_t nTracks = dst->numberOfTracks();
    tpc_tracks->Reserve(nTracks);

    // Track loop
    for(Int_t iTrk=0; iTrk<nTracks; iTrk++)
    {
      // Retrieve i-th femto track
      StFemtoTrack *femtoTrack = dst->track(iTrk);

      if (!femtoTrack) continue;

      auto *track = tpc_tracks->AddChannel();
      track->Init(out_config->GetBranchConfig(str_tpc_tracks_branch));

      // Reading Reco Track
      track->SetMomentum(femtoTrack->pMom().X(),femtoTrack->pMom().Y(),femtoTrack->pMom().Z());
      track->SetField(int(femtoTrack->nHits()), inHits);
      track->SetField(int(femtoTrack->nHitsFit()), inHitsFit);
      track->SetField(int(femtoTrack->nHitsPoss()), inHitsPoss);
      track->SetField(float(femtoTrack->nHitsFit()/femtoTrack->nHitsPoss()), inHits_Fit2Poss);
      track->SetField(float(femtoTrack->nSigmaPion()), inSigmaPion);
      track->SetField(float(femtoTrack->nSigmaKaon()), inSigmaKaon);
      track->SetField(float(femtoTrack->nSigmaProton()), inSigmaProton);
      track->SetField(float(femtoTrack->dEdx()), idEdx);
      track->SetField(float(femtoTrack->massSqr()), itofMass2);
      track->SetField(float(femtoTrack->gDCA(pVtx).Mag()), idca);
      track->SetField(float(femtoTrack->gDCA(pVtx).X()), idca_x);
      track->SetField(float(femtoTrack->gDCA(pVtx).Y()), idca_y);
      track->SetField(float(femtoTrack->gDCA(pVtx).Z()), idca_z);
      track->SetField(bool(femtoTrack->isPrimary()), iisPrimary);
      track->SetField(bool(femtoTrack->isTofTrack()), iisTofTrack);
    } //for(Int_t iTrk=0; iTrk<nTracks; iTrk++)

    // Fill BBC information
    bbc_modules->Reserve(NbbcModules);
    for (int imodule = 0; imodule<NbbcModules; imodule++)
    {
      auto *module = bbc_modules->AddChannel();
      module->Init(out_config->GetBranchConfig(str_bbc_branch));
      module->SetNumber(imodule);
      if (imodule < NbbcModulesHalf) module->SetSignal(event->bbcAdcEast(imodule));
      if (imodule >= NbbcModulesHalf) module->SetSignal(event->bbcAdcWest(imodule-NbbcModulesHalf));
    }

    // Fill ZDC-SMD information
    // There are 7 vertical and 8 horizontal strips of ZDC-SMD. To calculate Q-vectors from ZDC
    // one has to sum up signal from those strips separately:
    // Qx = sum zdcSmdVertical   * position_of_the_vertical_strip,   (sum over 7 strips)
    // Qy = sum zdcSmdHorizontal * position_of_the_Horizontal_strip, (sum over 8 strips)
    // To translate this information into AnalysisTree, HitDetector was chosen for horizontal
    // and vertical strips separately.
    zdc_horizontal_hits->Reserve(2*NzdcStipsY);
    zdc_vertical_hits->Reserve(2*NzdcStipsX);
    TVector3 geomZdcPos;

    // Fill ZDC-SMD East - horizontal
    for (int ihity=0; ihity<NzdcStipsY; ihity++)
    {
      auto *hit = zdc_horizontal_hits->AddChannel();
      hit->Init(out_config->GetBranchConfig(str_zdc_horizontal_branch));
      geomZdcPos = GetZdcPositionYEast(ihity);
      hit->SetPosition(geomZdcPos);
      hit->SetSignal(event->zdcSmdEastHorizontal(ihity));
    }
    // Fill ZDC-SMD East - vertical
    for (int ihitx=0; ihitx<NzdcStipsX; ihitx++)
    {
      auto *hit = zdc_vertical_hits->AddChannel();
      hit->Init(out_config->GetBranchConfig(str_zdc_vertical_branch));
      geomZdcPos = GetZdcPositionXEast(ihitx);
      hit->SetPosition(geomZdcPos);
      hit->SetSignal(event->zdcSmdEastVertical(ihitx));
    }
    // Fill ZDC-SMD West - horizontal
    for (int ihity=0; ihity<NzdcStipsY; ihity++)
    {
      auto *hit = zdc_horizontal_hits->AddChannel();
      hit->Init(out_config->GetBranchConfig(str_zdc_horizontal_branch));
      geomZdcPos = GetZdcPositionYWest(ihity);
      hit->SetPosition(geomZdcPos);
      hit->SetSignal(event->zdcSmdWestHorizontal(ihity));
    }
    // Fill ZDC-SMD West - vertical
    for (int ihitx=0; ihitx<NzdcStipsX; ihitx++)
    {
      auto *hit = zdc_vertical_hits->AddChannel();
      hit->Init(out_config->GetBranchConfig(str_zdc_vertical_branch));
      geomZdcPos = GetZdcPositionXWest(ihitx);
      hit->SetPosition(geomZdcPos);
      hit->SetSignal(event->zdcSmdWestVertical(ihitx));
    }

    outTree->Fill();
  } //for(Long64_t iEvent=0; iEvent<events2read; iEvent++)

  outFile->cd();
  outTree->Print();
  outTree->Write();
  if (isDataHeaderDefined) dataHeader->Write("DataHeader");
  out_config->Write("Configuration");
  outFile->Close();

  femtoReader->Finish();

  timer.Stop();
  timer.Print();

  return 0;
}

float get_beamP(float sqrtSnn, float m_target, float m_beam)
{
  return sqrt( pow((pow(sqrtSnn,2) - pow(m_target,2) - pow(m_beam,2))/(2*m_target), 2) - pow(m_beam,2) );
}

float GetCentFromCent9(int cent9)
{
  if (cent9 == 0) return 75.;
  if (cent9 == 1) return 65.;
  if (cent9 == 2) return 55.;
  if (cent9 == 3) return 45.;
  if (cent9 == 4) return 35.;
  if (cent9 == 5) return 25.;
  if (cent9 == 6) return 15.;
  if (cent9 == 7) return 7.5;
  if (cent9 == 8) return 2.5;

  return -1.;
}

float GetCentFromCent16(int cent16)
{
  if (cent16 == 0)  return 77.5;
  if (cent16 == 1)  return 72.5;
  if (cent16 == 2)  return 67.5;
  if (cent16 == 3)  return 62.5;
  if (cent16 == 4)  return 57.5;
  if (cent16 == 5)  return 52.5;
  if (cent16 == 6)  return 47.5;
  if (cent16 == 7)  return 42.5;
  if (cent16 == 8)  return 37.5;
  if (cent16 == 9)  return 32.5;
  if (cent16 == 10) return 27.5;
  if (cent16 == 11) return 22.5;
  if (cent16 == 12) return 17.5;
  if (cent16 == 13) return 12.5;
  if (cent16 == 14) return 7.5;
  if (cent16 == 15) return 2.5;

  return -1.;
}

//--------------------------------------------------------------------------------------------------------------------------//
//this function simply connects the gain values read in to the BBC azimuthal distribution
//since tiles 7 and 9 (+ 13 and 15) share a gain value it is ambiguous how to assign the geometry here
//It is preferred sometimes assigning the angle between the tiles thus "greying out" the adcs. 
//Others have assigned all of the adc to one (exclusive) or the the other. 
TVector3 GetBbcPositionEast(int imod)
{
  double z_pos = -375.; // in cm
  double radius=0.;;

  // Radius here is an artificial value just to distinguish between different rings of BBC modules.
  // It does not represent the true geometrical positions of the modules themself.
  if (imod >=0 && imod < 6) radius = 1.;
  if (imod >=6 && imod < 16) radius = 2.;

  //float GetPhiInBBC(int eastWest, int bbcN) { //imod=0 to 23
  const float Pi = TMath::Pi();
  const float Phi_div = Pi / 6;
  float bbc_phi = Phi_div;
  switch(imod) {
    case 0: bbc_phi = 3.*Phi_div;
  break;
    case 1: bbc_phi = Phi_div;
  break;
    case 2: bbc_phi = -1.*Phi_div;
  break;
    case 3: bbc_phi = -3.*Phi_div;
  break;
    case 4: bbc_phi = -5.*Phi_div;
  break;
    case 5: bbc_phi = 5.*Phi_div;
  break;
    //case 6: bbc_phi= (mRndm.Rndm() > 0.5) ? 2.*Phi_div:4.*Phi_div;	//tiles 7 and 9 are gained together we randomly assign the gain to one XOR the other
    case 6: bbc_phi = 3.*Phi_div;
  break;
    case 7: bbc_phi = 3.*Phi_div;
  break;
    case 8: bbc_phi = Phi_div;
  break;
    case 9: bbc_phi = 0.;
  break;
    case 10: bbc_phi = -1.*Phi_div;
  break;
    //case 11: bbc_phi = (mRndm.Rndm() > 0.5) ? -2.*Phi_div:-4.*Phi_div;	//tiles 13 and 15 are gained together
    case 11: bbc_phi = -3.*Phi_div;
  break;
    case 12: bbc_phi = -3.*Phi_div;
  break;
    case 13: bbc_phi = -5.*Phi_div;
  break;
    case 14: bbc_phi = Pi;
  break;
    case 15: bbc_phi = 5.*Phi_div;
  break;
  }

  //if we're looking at the east BBC we need to flip around x in the STAR coordinates, 
  //a line parallel to the beam would go through tile 1 on the W BBC and tile 3 on the 
  // E BBC
  if (bbc_phi > -0.001){ //this is not a >= since we are talking about finite adcs -- not to important
    bbc_phi = Pi - bbc_phi;
  }
  else {
    bbc_phi= -Pi - bbc_phi;
  }

  if(bbc_phi < 0.0) bbc_phi += 2.*Pi;
  if(bbc_phi > 2.*Pi) bbc_phi -= 2.*Pi;

  TVector3 vresult;
  vresult.SetXYZ(radius*TMath::Cos(bbc_phi), radius*TMath::Sin(bbc_phi), z_pos);

  return vresult;
}

TVector3 GetBbcPositionWest(int imod)
{
  double z_pos = 375.; // in cm
  double radius=0.;;

  // Radius here is an artificial value just to distinguish between different rings of BBC modules.
  // It does not represent the true geometrical positions of the modules themself.
  if (imod >=0 && imod < 6) radius = 1.;
  if (imod >=6 && imod < 16) radius = 2.;

  //float GetPhiInBBC(int eastWest, int bbcN) { //imod=0 to 23
  const float Pi = TMath::Pi();
  const float Phi_div = Pi / 6;
  float bbc_phi = Phi_div;
  switch(imod) {
    case 0: bbc_phi = 3.*Phi_div;
  break;
    case 1: bbc_phi = Phi_div;
  break;
    case 2: bbc_phi = -1.*Phi_div;
  break;
    case 3: bbc_phi = -3.*Phi_div;
  break;
    case 4: bbc_phi = -5.*Phi_div;
  break;
    case 5: bbc_phi = 5.*Phi_div;
  break;
    //case 6: bbc_phi= (mRndm.Rndm() > 0.5) ? 2.*Phi_div:4.*Phi_div;	//tiles 7 and 9 are gained together we randomly assign the gain to one XOR the other
    case 6: bbc_phi = 3.*Phi_div;
  break;
    case 7: bbc_phi = 3.*Phi_div;
  break;
    case 8: bbc_phi = Phi_div;
  break;
    case 9: bbc_phi = 0.;
  break;
    case 10: bbc_phi = -1.*Phi_div;
  break;
    //case 11: bbc_phi = (mRndm.Rndm() > 0.5) ? -2.*Phi_div:-4.*Phi_div;	//tiles 13 and 15 are gained together
    case 11: bbc_phi = -3.*Phi_div;
  break;
    case 12: bbc_phi = -3.*Phi_div;
  break;
    case 13: bbc_phi = -5.*Phi_div;
  break;
    case 14: bbc_phi = Pi;
  break;
    case 15: bbc_phi = 5.*Phi_div;
  break;
  }

  if(bbc_phi < 0.0) bbc_phi += 2.*Pi;
  if(bbc_phi > 2.*Pi) bbc_phi -= 2.*Pi;

  TVector3 vresult;
  vresult.SetXYZ(radius*TMath::Cos(bbc_phi), radius*TMath::Sin(bbc_phi), z_pos);

  return vresult;
}

// Get position of each slat;strip starts from 0
TVector3 GetZdcPositionXEast(int strip)
{
  TVector3 vresult;
  double x_pos, y_pos;
  double z_pos = -1800.; // in cm

  Double_t zdcsmd_x[7] = {0.5,2,3.5,5,6.5,8,9.5};

  Double_t mZDCSMDCenterex = 4.72466;
  Double_t mZDCSMDCenterey = 5.53629;

  y_pos = mZDCSMDCenterey;

  if (strip>=0 && strip<7)
  {
    x_pos = zdcsmd_x[strip]-mZDCSMDCenterex;
  }
  else
  {
    x_pos = 0.;
  }

  vresult.SetXYZ(x_pos, y_pos, z_pos);
  return vresult;
}

TVector3 GetZdcPositionYEast(int strip)
{
  TVector3 vresult;
  double x_pos, y_pos;
  double z_pos = -1800.; // in cm

  Double_t zdcsmd_y[8] = {1.25,3.25,5.25,7.25,9.25,11.25,13.25,15.25};

  Double_t mZDCSMDCenterex = 4.72466;
  Double_t mZDCSMDCenterey = 5.53629;

  x_pos = mZDCSMDCenterex;

  if (strip>=0 && strip<8)
  {
    y_pos = zdcsmd_y[strip]/TMath::Sqrt(2.)-mZDCSMDCenterey;
  }
  else
  {
    x_pos = 0.;
  }

  vresult.SetXYZ(x_pos, y_pos, z_pos);
  return vresult;
}

// Get position of each slat;strip starts from 0
TVector3 GetZdcPositionXWest(int strip)
{
  TVector3 vresult;
  double x_pos, y_pos;
  double z_pos = 1800.; // in cm

  Double_t zdcsmd_x[7] = {0.5,2,3.5,5,6.5,8,9.5};

  Double_t mZDCSMDCenterwx = 4.39604;
  Double_t mZDCSMDCenterwy = 5.19968;

  y_pos = mZDCSMDCenterwy;

  if (strip>=0 && strip<7)
  {
    x_pos = mZDCSMDCenterwx-zdcsmd_x[strip];
  }
  else
  {
    x_pos = 0.;
  }

  vresult.SetXYZ(x_pos, y_pos, z_pos);
  return vresult;
}

// Get position of each slat;strip starts from 0
TVector3 GetZdcPositionYWest(int strip)
{
  TVector3 vresult;
  double x_pos, y_pos;
  double z_pos = 1800.; // in cm

  Double_t zdcsmd_y[8] = {1.25,3.25,5.25,7.25,9.25,11.25,13.25,15.25};

  Double_t mZDCSMDCenterwx = 4.39604;
  Double_t mZDCSMDCenterwy = 5.19968;

  x_pos = mZDCSMDCenterwx;

  if (strip>=0 && strip<8)
  {
    y_pos = zdcsmd_y[strip]/sqrt(2.)-mZDCSMDCenterwy;
  }
  else
  {
    y_pos = 0.;
  }

  vresult.SetXYZ(x_pos, y_pos, z_pos);
  return vresult;
}