TExogamPhysics.cxx

/*****************************************************************************
 * Copyright (C) 2009-2016   this file is part of the NPTool Project         *
 *                                                                           *
 * For the licensing terms see $NPTOOL/Licence/NPTool_Licence                *
 * For the list of contributors see $NPTOOL/Licence/Contributors             *
 *****************************************************************************/

/*****************************************************************************
 * Original Author: Sandra GIRON  contact address: giron@ipno.in2p3.fr       *
 *                  Benjamin LE CROM		   lecrom@ipno.in2p3.fr              *
 * Creation Date  : march 2014                                               *
 * Last update    : updated in 2023-2024 by H. Jacob hjacob@ijclab.in2p3.fr  *
 *---------------------------------------------------------------------------*
 * Decription:                                                               *
 *  This class hold exogam treated data                                      *
 *                                                                           *
 *---------------------------------------------------------------------------*
 * Comment:                                                                  *
 *                                                                           *
 *****************************************************************************/
#include "TExogamPhysics.h"
using namespace EXOGAM_LOCAL;

//	STL
#include <cmath>
#include <iostream>
#include <sstream>
#include <stdlib.h>
#include <functional>

//	NPL
#include "NPDetectorFactory.h"
#include "NPOptionManager.h"
#include "NPVDetector.h"
#include "RootInput.h"
#include "RootOutput.h"
//	ROOT
#include "TChain.h"

///////////////////////////////////////////////////////////////////////////
ClassImp(TExogamPhysics)
  ///////////////////////////////////////////////////////////////////////////
  TExogamPhysics::TExogamPhysics()
  : m_PreTreatedData(new TExogamCalData),
  m_EventData(new TExogamData),
  TSEvent(new TimeStamp),
  m_EventPhysics(this)
  {
    // m_Spectra = NULL;
    m_EXO_E_RAW_Threshold = 0;
    m_EXO_E_Threshold = 0;
    m_EXO_EHG_RAW_Threshold = 0;
    m_EXO_TDC_RAW_Threshold = 0;
    m_ExoTDC_HighThreshold = 1e6;
    m_ExoTDC_LowThreshold = 0;
    m_EXO_OuterUp_RAW_Threshold = 1e5;

  DataIsCal = false;
}
  
TExogamPhysics::~TExogamPhysics() {
  delete m_PreTreatedData;
  delete m_EventData;
  // delete TSEvent; 
}

///////////////////////////////////////////////////////////////////////////
void TExogamPhysics::BuildSimplePhysicalEvent() { BuildPhysicalEvent(); }

///////////////////////////////////////////////////////////////////////////
void TExogamPhysics::PreTreat() {
  // Clearing PreTreat TExogamData
  ClearPreTreatedData();

  //E
  m_EXO_Mult = m_EventData->GetExoMult();

  for (unsigned int i = 0; i < m_EXO_Mult; ++i) {
      bool DoPreTreat = false;
      if(!RefTS_Name.empty()){
        
        std::string TSName = "EXO_"+std::to_string(m_EventData->GetExoCrystal(i));
        TSEvent->AddTimeStamp(TSName,m_EventData->GetExoTS(i));
        TSEvent->AddTimeStamp(RefTS_Name,RefTS);
        
        if(TSEvent->MatchTS(TSName)){
          DoPreTreat = true;
        }
        TSEvent->ClearTimeStamps();
      }
      
      // else, all datas are filled
      else{
        DoPreTreat = true;
      }
    if(DoPreTreat){ 
      ResetPreTreatVariable();
    
      if (m_EventData->GetExoE(i) > m_EXO_E_RAW_Threshold)
        EXO_E = fEXO_E(m_EventData, i);

      if (m_EventData->GetExoEHG(i) > m_EXO_EHG_RAW_Threshold)
        EXO_EHG = fEXO_EHG(m_EventData, i);
    
      if (m_EventData->GetExoTDC(i) > m_EXO_TDC_RAW_Threshold)
        EXO_TDC = fEXO_T(m_EventData, i);
 
      if (m_EventData->GetExoOuter1(i) < m_EXO_OuterUp_RAW_Threshold)
        EXO_Outer1 = fEXO_Outer(m_EventData, i, 0);
      else
        EXO_Outer1 = 0;
    
      if (m_EventData->GetExoOuter2(i) < m_EXO_OuterUp_RAW_Threshold)
        EXO_Outer2 = fEXO_Outer(m_EventData, i, 1);
      else
        EXO_Outer2 = 0;
    
      if (m_EventData->GetExoOuter3(i) < m_EXO_OuterUp_RAW_Threshold)
        EXO_Outer3 = fEXO_Outer(m_EventData, i, 2);
      else
        EXO_Outer3 = 0;
    
      if (m_EventData->GetExoOuter4(i) < m_EXO_OuterUp_RAW_Threshold)
        EXO_Outer4 = fEXO_Outer(m_EventData, i, 3);
      else
        EXO_Outer4 = 0;
    
      // *1000 to convert MeV into keV
      if(EXO_E > m_EXO_E_Threshold){
        m_PreTreatedData->SetExo(m_EventData->GetExoCrystal(i), EXO_E,
        EXO_EHG, m_EventData->GetExoTS(i), EXO_TDC, 
        m_EventData->GetExoBGO(i), m_EventData->GetExoCsI(i), EXO_Outer1,
        EXO_Outer2, EXO_Outer3, EXO_Outer4);
      }
    } 
  }
}

///////////////////////////////////////////////////////////////////////////
void TExogamPhysics::ResetPreTreatVariable(){
  EXO_E = -1000;
  EXO_EHG = -1000;
  EXO_TDC = -1000;
  EXO_Outer1 = -1000;
  EXO_Outer2 = -1000;
  EXO_Outer3 = -1000;
  EXO_Outer4 = -1000;
}

///////////////////////////////////////////////////////////////////////////
void TExogamPhysics::BuildPhysicalEvent() {
  ClaimReaderData();
  
  // std::cout << m_EventData << std::endl;
  if(!DataIsCal)
    PreTreat();

  // This maps stores ID of events sorted by flange number. Map key is flange nbr, vector should contain ID of events
  std::map<unsigned int,std::vector<unsigned int>> HitsID;

  for(unsigned int i = 0; i < m_PreTreatedData->GetExoMult(); i++){
    // Asking good TDC prompt, which can be ignored in EXOGAM config files (TDC not working for some good crystals in E805) 
    if(TDCMatch(i) ||  find(IgnoreTDC.begin(), IgnoreTDC.end(), m_PreTreatedData->GetExoCrystal(i)) != IgnoreTDC.end()){
      // Doing flange and crystal matching
      flange_nbr = MapCrystalFlangeCLover[m_PreTreatedData->GetExoCrystal(i)].first;
      crystal_nbr = MapCrystalFlangeCLover[m_PreTreatedData->GetExoCrystal(i)].second;
      E.push_back(m_PreTreatedData->GetExoE(i));
      EHG.push_back(m_PreTreatedData->GetExoEHG(i));
      Outer1.push_back(m_PreTreatedData->GetExoOuter1(i));
      Outer2.push_back(m_PreTreatedData->GetExoOuter2(i));
      Outer3.push_back(m_PreTreatedData->GetExoOuter3(i));
      Outer4.push_back(m_PreTreatedData->GetExoOuter4(i));
      TDC.push_back(m_PreTreatedData->GetExoTDC(i));
      TS.push_back(m_PreTreatedData->GetExoTS(i));
      Flange.push_back(flange_nbr);
      Crystal.push_back(crystal_nbr);

      HitsID[flange_nbr].push_back(i);
    }
  }

  // Now that HitsID is full, we use it to process simple AddBack of events in the same flange
  // Basically looping on all flanges, then on al events ID in each flange
  for(auto it = HitsID.begin(); it != HitsID.end(); it++){
    double E_AddBack = 0;
    double E_Max = 0;
    unsigned int Id_Max = 0;
    for(auto itvec = (*it).second.begin(); itvec !=(*it).second.end(); itvec++){
      E_AddBack+= m_PreTreatedData->GetExoE(*itvec);
      if(E_Max < m_PreTreatedData->GetExoE(*itvec)){
        E_Max = m_PreTreatedData->GetExoE(*itvec);
        Id_Max = *itvec;
      }
    }
    // Doing it again for this loop, it's a bit unhappy but didnt find a better way to do it yet
    flange_nbr = (*it).first;
    crystal_nbr = MapCrystalFlangeCLover[m_PreTreatedData->GetExoCrystal(Id_Max)].second;

    // Adding all AddBack (AB) related stuff
    E_AB.push_back(E_AddBack);
    Flange_AB.push_back(flange_nbr);
    Size_AB.push_back((*it).second.size());
    TDC_AB.push_back(m_PreTreatedData->GetExoTDC(Id_Max));
    TS_AB.push_back(m_PreTreatedData->GetExoTS(Id_Max));

    // Adding these parameters for Doppler correction purposes (D)
    Crystal_AB.push_back(crystal_nbr);
    int MaxOuterId = GetMaxOuter(Id_Max);
    Outer_AB.push_back(GetMaxOuter(Id_Max));

    if(MaxOuterId > -1){
      Exogam_struc = Ask_For_Angles(flange_nbr, ComputeMeanFreePath(E_AddBack),147,0.0); //147 default value of Emmanuel's code
      double Theta_seg = Exogam_struc.Theta_Crystal_Seg[crystal_nbr][MaxOuterId];
      double Phi_seg = Exogam_struc.Phi_Crystal_Seg[crystal_nbr][MaxOuterId];

      Theta.push_back(Theta_seg);
      Phi.push_back(Phi_seg);
    }
    else{
      Theta.push_back(-1000);
      Phi.push_back(-1000);
    }
  }
}

void TExogamPhysics::ClaimReaderData() {
  if (NPOptionManager::getInstance()->IsReader() == true) {
    m_EventData = &(**r_ReaderEventData);
  }
}

///////////////////////////////////////////////////////////////////////////
bool TExogamPhysics::TDCMatch(unsigned int event){
  return m_PreTreatedData->GetExoTDC(event) > m_ExoTDC_LowThreshold && m_PreTreatedData->GetExoTDC(event) < m_ExoTDC_HighThreshold;
}

///////////////////////////////////////////////////////////////////////////
int TExogamPhysics::GetMaxOuter(unsigned int EventId){
  // somehow starting at 50 to get something equivalent to a 50keV threshold
  double OuterMax = 50;
  int OuterId = -1;
  if(m_PreTreatedData->GetExoOuter1(EventId) > OuterMax){
    OuterMax =  m_PreTreatedData->GetExoOuter1(EventId);
    OuterId = 0;
  }
  if(m_PreTreatedData->GetExoOuter2(EventId) > OuterMax){
    OuterMax =  m_PreTreatedData->GetExoOuter2(EventId);
    OuterId = 1;
  }
  if(m_PreTreatedData->GetExoOuter3(EventId) > OuterMax){
    OuterMax =  m_PreTreatedData->GetExoOuter3(EventId);
    OuterId = 2;
  }
  if(m_PreTreatedData->GetExoOuter4(EventId) > OuterMax){
    OuterMax =  m_PreTreatedData->GetExoOuter4(EventId);
    OuterId = 3;
  }
  return OuterId;
}

///////////////////////////////////////////////////////////////////////////
double TExogamPhysics::GetDoppler(double Energy, unsigned int Flange, unsigned int Crystal, unsigned int Outer){
  Exogam_struc = Ask_For_Angles(Flange, ComputeMeanFreePath(Energy),147.0,0.0); // same
  double Theta_seg = Exogam_struc.Theta_Crystal_Seg[Crystal][Outer];
  double Phi_seg = Exogam_struc.Phi_Crystal_Seg[Crystal][Outer];
  return Doppler_Correction(Theta_seg,Phi_seg,0,0,Beta,Energy);
}

///////////////////////////////////////////////////////////////////////////
double TExogamPhysics::ComputeMeanFreePath(double Energy){
  auto b = Map_PhotonCS.lower_bound(Energy);
  auto a = prev(b);
  if(b == Map_PhotonCS.begin()){
    a = b;
    b++;
  }
  else if(b == Map_PhotonCS.end()){
    b--;
    a = prev(b);
  }
  double coeff = (Energy - a->first)/(b->first - a->first);

  double PhotonCrossSection = a->second + coeff*(b->second - a->second); // mm2/g
  return 1./(GeDensity*PhotonCrossSection);
}

// unsigned int TExogamPhysics::GetFlangeNbr(unsigned int crystal_nbr){

// }

///////////////////////////////////////////////////////////////////////////
double TExogamPhysics::DopplerCorrection(double E, double Theta) {
  double Pi = 3.141592654;
  TString filename = "configs/beta.txt";
  ifstream file;
  // cout << filename << endl;
  file.open(filename);
  if (!file)
    cout << filename << " was not opened" << endl;

  double E_corr = 0;
  double beta = 0.;
  file >> beta;
  double gamma = 1. / sqrt(1 - beta * beta);

  E_corr = gamma * E * (1. - beta * cos(Theta * Pi / 180.));

  return (E_corr);
}

///////////////////////////////////////////////////////////////////////////
// Routine of doppler correction
///////////////////////////////////////////////////////////////////////////
double TExogamPhysics::CorrectionDoppler(double theta_gamma, double phi_gamma, double theta_part, double phi_part, double beta_part, double E){  //rad, v/c
  double Ecorr,cosinusPSI;

  cosinusPSI =TMath::Sin(theta_part)*TMath::Cos(phi_part)*TMath::Sin(theta_gamma)*TMath::Cos(phi_gamma)+
    TMath::Sin(theta_part)*TMath::Sin(phi_part)*TMath::Sin(theta_gamma)*TMath::Sin(phi_gamma)+
    TMath::Cos(theta_part)*TMath::Cos(theta_gamma);

  Ecorr = E*(1.-beta_part*cosinusPSI)/sqrt(1.-beta_part*beta_part);

  return Ecorr;
}

///////////////////////////////////////////////////////////////////////////
void TExogamPhysics::Clear() {
  // Exogam_struc = {};

  E.clear();
  EHG.clear();
  Outer1.clear();
  Outer2.clear();
  Outer3.clear();
  Outer4.clear();
  Flange.clear();
  Crystal.clear();
  TDC.clear();
  TS.clear();

  E_AB.clear();
  Flange_AB.clear();
  Size_AB.clear();
  Crystal_AB.clear();
  Outer_AB.clear();
  Theta.clear();
  Phi.clear();
  TDC_AB.clear();
  TS_AB.clear();
}

///////////////////////////////////////////////////////////////////////////
//	Read stream at ConfigFile to pick-up parameters of detector (Position,...) using Token
///////////////////////////////////////////////////////////////////////////
void TExogamPhysics::ReadConfiguration(NPL::InputParser parser) {
  vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Exogam");
  if (NPOptionManager::getInstance()->GetVerboseLevel())
    cout << "//// " << blocks.size() << " Exogam clover found " << endl;

  for(unsigned int i=0; i<blocks.size(); i++){
    if(NPOptionManager::getInstance()->GetVerboseLevel())
      cout << endl << "//// EXOGAM Clover " << i+1 << endl;

    m_NumberOfClovers++;
    int flange = blocks[i]->GetInt("FLANGE");
    m_flange.push_back(flange);
    m_pos_segment[0][0].push_back(blocks[i]->GetTVector3("POS_CRYSTALA_SEG1","mm"));
    m_pos_segment[0][1].push_back(blocks[i]->GetTVector3("POS_CRYSTALA_SEG2","mm"));
    m_pos_segment[0][2].push_back(blocks[i]->GetTVector3("POS_CRYSTALA_SEG3","mm"));
    m_pos_segment[0][3].push_back(blocks[i]->GetTVector3("POS_CRYSTALA_SEG4","mm"));
    m_pos_segment[1][0].push_back(blocks[i]->GetTVector3("POS_CRYSTALB_SEG1","mm"));
    m_pos_segment[1][1].push_back(blocks[i]->GetTVector3("POS_CRYSTALB_SEG2","mm"));
    m_pos_segment[1][2].push_back(blocks[i]->GetTVector3("POS_CRYSTALB_SEG3","mm"));
    m_pos_segment[1][3].push_back(blocks[i]->GetTVector3("POS_CRYSTALB_SEG4","mm"));
    m_pos_segment[2][0].push_back(blocks[i]->GetTVector3("POS_CRYSTALC_SEG1","mm"));
    m_pos_segment[2][1].push_back(blocks[i]->GetTVector3("POS_CRYSTALC_SEG2","mm"));
    m_pos_segment[2][2].push_back(blocks[i]->GetTVector3("POS_CRYSTALC_SEG3","mm"));
    m_pos_segment[2][3].push_back(blocks[i]->GetTVector3("POS_CRYSTALC_SEG4","mm"));
    m_pos_segment[3][0].push_back(blocks[i]->GetTVector3("POS_CRYSTALD_SEG1","mm"));
    m_pos_segment[3][1].push_back(blocks[i]->GetTVector3("POS_CRYSTALD_SEG2","mm"));
    m_pos_segment[3][2].push_back(blocks[i]->GetTVector3("POS_CRYSTALD_SEG3","mm"));
    m_pos_segment[3][3].push_back(blocks[i]->GetTVector3("POS_CRYSTALD_SEG4","mm"));

    MapFlangeToCloverNumber[flange] = m_NumberOfClovers; 
  }

  ReadAnalysisConfig();
}

///////////////////////////////////////////////////////////////////////////
void TExogamPhysics::ReadAnalysisConfig() {
  bool ReadingStatus = false;


  // path to photon cross section
  string CSFilename = "../../Inputs/PhotonCrossSection/CoherentGe.xcom";
  string LineBuffer;

  ifstream CSFile;
  CSFile.open(CSFilename.c_str());

  if (!CSFile.is_open()) {
    cout << " No CS file found "
      << CSFilename << endl;
    return;
  }
  while(CSFile.good()){
    double gammaE, CrossSection;
    getline(CSFile, LineBuffer);
    istringstream ss(LineBuffer);
    ss >> gammaE >> CrossSection; // E in MeV, converted to keV, CrossSection in cm2/g  
    gammaE *= 1000.; // Convertion to keV
    CrossSection *= 100.;
    Map_PhotonCS[gammaE] = CrossSection;
  }
  // path to file
  string FileName = "./configs/ConfigExogam.dat";

  // open analysis config file
  ifstream AnalysisConfigFile;
  AnalysisConfigFile.open(FileName.c_str());

  if (!AnalysisConfigFile.is_open()) {
    cout << " No ConfigExogam.dat found: Default parameters loaded for "
      "Analysis "
      << FileName << endl;
    return;
  }


  string DataBuffer, whatToDo;
  while (!AnalysisConfigFile.eof()) {
    // Pick-up next line
    getline(AnalysisConfigFile, LineBuffer);

    // search for "header"
    if (LineBuffer.compare(0, 12, "ConfigExogam") == 0)
      ReadingStatus = true;

    // loop on tokens and data
    while (ReadingStatus) {

      whatToDo = "";
      AnalysisConfigFile >> whatToDo;
      // Search for comment symbol (%)
      if (whatToDo.compare(0, 1, "%") == 0) {
        AnalysisConfigFile.ignore(numeric_limits<streamsize>::max(), '\n');
      }

      else if (whatToDo == "EXO_Threshold") {
        //AnalysisConfigFile >> DataBuffer;
        //m_MaximumStripMultiplicityAllowed = atoi(DataBuffer.c_str());
        //cout << "MAXIMUN STRIP MULTIPLICITY " << m_MaximumStripMultiplicityAllowed << endl;
      }
      else if (whatToDo=="MAP_EXO") {
        unsigned int CrystalNb;
        unsigned int Flange;
        unsigned int CrystalNb2;

        AnalysisConfigFile >> DataBuffer;
        CrystalNb = stoi(DataBuffer);

        AnalysisConfigFile >> DataBuffer;
        Flange = stoi(DataBuffer);

        AnalysisConfigFile >> DataBuffer;
        CrystalNb2 = stoi(DataBuffer);
        MapCrystalFlangeCLover[CrystalNb] = std::make_pair(Flange,CrystalNb2);
        cout << whatToDo << " / " << "crystal= " << CrystalNb << " && flange= " << Flange << endl;
      }
      else if (whatToDo == "TDC_THRESHOLDS") {
        AnalysisConfigFile >> DataBuffer;
        m_ExoTDC_LowThreshold = stoi(DataBuffer);
        AnalysisConfigFile >> DataBuffer;
        m_ExoTDC_HighThreshold = stoi(DataBuffer);
        cout << "TDC Thresholds " << m_ExoTDC_LowThreshold << " " <<m_ExoTDC_HighThreshold << endl;
      }
      else if (whatToDo == "IGNORE_TDC") {
        AnalysisConfigFile >> DataBuffer;
        IgnoreTDC.push_back(stoi(DataBuffer));
        cout << "TDC Ignored for Crystals : " << DataBuffer << endl;
      }
      else if (whatToDo == "DATA_IS_CAL") {
        AnalysisConfigFile >> DataBuffer;
        DataIsCal = (stoi(DataBuffer) == 1);
        if(DataIsCal)
          cout << "Using Calibrated Data for Exogam" << endl;
      }
      else{
        ReadingStatus = false;
      }
    }
  }
}
void TExogamPhysics::InitSpectra() {
}
///////////////////////////////////////////////////////////////////////////
void TExogamPhysics::FillSpectra() {
}
///////////////////////////////////////////////////////////////////////////
void TExogamPhysics::CheckSpectra() { m_Spectra->CheckSpectra(); }
///////////////////////////////////////////////////////////////////////////
void TExogamPhysics::ClearSpectra() {
  // To be done
}
///////////////////////////////////////////////////////////////////////////
map<string, TH1*> TExogamPhysics::GetSpectra() {
  if (m_Spectra)
    return m_Spectra->GetMapHisto();
  else {
    map<string, TH1*> empty;
    return empty;
  }
}


//////////////////////////////////////////////////////////////////////////
void TExogamPhysics::AddClover(int Board, int Flange, int Channel0, int Channel1) {

}

// FIXME Legacy thing... Might delete later
//////////////////////////////////////////////////////////////////////////
// void TExogamPhysics::AddClover(string AngleFile) {
//   ifstream file;
//   //  TString filename = Form("posBaptiste/angles_exogam_clover%d.txt",NumberOfClover);
//   //  TString filename = Form("posz42_simu50mm/angles_exogam_clover%d.txt",NumberOfClover);
//   //  TString filename = Form("posz42_exp_stat_demiring/angles_exogam_clover%d.txt",NumberOfClover);

//   string path = "configs/";
//   TString filename = path + AngleFile;

//   cout << filename << endl;
//   file.open(filename);
//   if (!file)
//     cout << filename << " was not opened" << endl;

//   vector<double> Angles;
//   vector<vector<double>> Segment_angles;
//   vector<vector<vector<double>>> Cristal_angles;

//   Cristal_angles.clear();

//   double angle;
//   string buffer;

//   for (int i = 0; i < 4; i++) {
//     Segment_angles.clear();

//     for (int j = 0; j < 4; j++) {
//       Angles.clear();

//       for (int k = 0; k < 2; k++) {
//         file >> buffer >> angle;

//         Angles.push_back(angle); // Theta (k = 0)   Phi (k = 1)

//         // cout << angle << endl;
//         if (Angles.size() == 2)
//           cout << "Clover " << NumberOfClover << ": Theta=" << Angles[0] << " Phi=" << Angles[1] << endl;
//       }

//       Segment_angles.push_back(Angles);
//     }

//     Cristal_angles.push_back(Segment_angles);
//   }

//   Clover_Angles_Theta_Phi.push_back(Cristal_angles);

//   file.close();

//   NumberOfClover++;
// }

//////////////////////////////////////////////////////////////////////////
//	Add Parameter to the CalibrationManger
//////////////////////////////////////////////////////////////////////////
void TExogamPhysics::AddParameterToCalibrationManager() {

  CalibrationManager* Cal = CalibrationManager::getInstance();

  for (auto it = MapCrystalFlangeCLover.begin(); it != MapCrystalFlangeCLover.end(); it++)
  {  unsigned int i = it->first;
    Cal->AddParameter("EXO", "E" + NPL::itoa(i),
        "EXO_E" + NPL::itoa(i));
    Cal->AddParameter("EXO", "EHG" + NPL::itoa(i),
        "EXO_EHG" + NPL::itoa(i));
    // Cal->AddParameter("EXOGAM", "Cl" + NPL::itoa(i) + "_Cr" + NPL::itoa(j) + "_T",
    // "EXOGAM_Cl" + NPL::itoa(i) + "_Cr" + NPL::itoa(j) + "_T");

    for (int j = 0; j < 4; j++) {
      Cal->AddParameter("EXO", "Outer" + NPL::itoa(i) + "_" + NPL::itoa(j),
          "EXO_Outer" + NPL::itoa(i) + "_" + NPL::itoa(j));
    }
  }
}


//////////////////////////////////////////////////////////////////////////
//	Activated associated Branches and link it to the private member DetectorData address
//	In this method mother Branches (Detector) AND daughter leaf (fDetector_parameter) have to be activated
//////////////////////////////////////////////////////////////////////////
void TExogamPhysics::InitializeRootInputRaw() {
  TChain* inputChain = RootInput::getInstance()->GetChain();
  // Option to use the nptreereader anaysis
  if (NPOptionManager::getInstance()->IsReader() == true) {
    TTreeReader* inputTreeReader = RootInput::getInstance()->GetTreeReader();
    inputTreeReader->SetTree(inputChain);
  }
  // Option to use pre calibrated NPTOOL data
  else if (DataIsCal) {
    inputChain->SetBranchStatus("Exogam", true);
    inputChain->SetBranchStatus("cEXO_*", true);
    inputChain->SetBranchAddress("Exogam", &m_PreTreatedData); 
  }
  // Option to use the standard npanalysis
  else{
    inputChain->SetBranchStatus("Exogam", true);
    inputChain->SetBranchStatus("fEXO_*", true);
    inputChain->SetBranchAddress("Exogam", &m_EventData);

    /*
       TList* outputList = RootOutput::getInstance()->GetList();
       clover_mult = new TH1F("clover_mult","clover_mult",20,0,20);
       outputList->Add(clover_mult);
       cristal_mult = new TH1F("cristal_mult","cristal_mult",20,0,20);
       outputList->Add(cristal_mult);
       */
  }
}

///////////////////////////////////////////////////////////////////////////
void TExogamPhysics::ReadConfigurationTS(){
  TSEvent->ReadConfigurationFile();
}

///////////////////////////////////////////////////////////////////////////
void TExogamPhysics::SetRefTS(std::string TSRef_Name, unsigned long long TSRef){
  RefTS = TSRef;
  RefTS_Name = TSRef_Name; 
}


/////////////////////////////////////////////////////////////////////
//   Activated associated Branches and link it to the private member DetectorPhysics address
//   In this method mother Branches (Detector) AND daughter leaf (parameter) have to be activated
void TExogamPhysics::InitializeRootInputPhysics() {
  TChain* inputChain = RootInput::getInstance()->GetChain();
  // Option to use the nptreereader anaysis
  if (NPOptionManager::getInstance()->IsReader() == true) {
    TTreeReader* inputTreeReader = RootInput::getInstance()->GetTreeReader();
    inputTreeReader->SetTree(inputChain);
  }
  // Option to use the standard npanalysis
  else{
    TChain* inputChain = RootInput::getInstance()->GetChain();
    inputChain->SetBranchStatus("Exogam" , true );
    inputChain->SetBranchStatus("EventMultiplicty", true);
    inputChain->SetBranchStatus("ECC_Multiplicity", true);
    inputChain->SetBranchStatus("GOCCE_Multiplicity", true);
    inputChain->SetBranchStatus("ECC_CloverNumber", true);
    inputChain->SetBranchStatus("ECC_CristalNumber", true);
    inputChain->SetBranchStatus("GOCCE_CloverNumber", true);
    inputChain->SetBranchStatus("GOCCE_CristalNumber", true);
    inputChain->SetBranchStatus("GOCCE_SegmentNumber", true);
    inputChain->SetBranchStatus("ECC_E", true);
    inputChain->SetBranchStatus("ECC_T", true);
    inputChain->SetBranchStatus("GOCCE_E", true);
    inputChain->SetBranchStatus("CristalNumber", true);
    inputChain->SetBranchStatus("SegmentNumber", true);
    inputChain->SetBranchStatus("CloverNumber", true);
    inputChain->SetBranchStatus("CloverMult", true);
    inputChain->SetBranchStatus("TotalEnergy_lab", true);
    inputChain->SetBranchStatus("Time", true);
    inputChain->SetBranchStatus("DopplerCorrectedEnergy", true);
    inputChain->SetBranchStatus("Position", true);
    inputChain->SetBranchStatus("Theta", true);
    inputChain->SetBranchAddress("Exogam", &m_EventPhysics);
  }
}

/////////////////////////////////////////////////////////////////////
//	Create associated branches and associated private member DetectorPhysics address
/////////////////////////////////////////////////////////////////////
void TExogamPhysics::InitializeRootOutput() {
  TTree* outputTree = RootOutput::getInstance()->GetTree();
  outputTree->Branch("Exogam", "TExogamPhysics", &m_EventPhysics);

  // control histograms if needed
  /*
     TList* outputList = RootOutput::getInstance()->GetList();
     controle = new TH1F("controle","histo de controle",20,0,20);
     outputList->Add(controle);
     */
}

/////////////////////////////////////////////////////////////////////
void TExogamPhysics::SetTreeReader(TTreeReader* TreeReader) {
  TExogamPhysicsReader::r_SetTreeReader(TreeReader);
}

/////////////////////////////// DoCalibration Part //////////////////////////:
void TExogamPhysics::InitializeRootHistogramsCalib() {
  std::cout << "Initialize Exogam Histograms" << std::endl;
  map<int, bool>::iterator it;
  map<int, map<int,bool>>::iterator it2;
  for (it = DoCalibrationE.begin(); it != DoCalibrationE.end(); it++) {
    if (it->second) {
      InitializeRootHistogramsE_F(it->first);
    }
  }
  for (it = DoCalibrationEHG.begin(); it != DoCalibrationEHG.end(); it++) {
    if (it->second) {
      InitializeRootHistogramsEHG_F(it->first);
    }
  }
  //for (it = DoCalibrationT.begin(); it != DoCalibrationT.end(); it++) {
  //  if (it->second) {
  //    InitializeRootHistogramsT_F(it->first);
  //  }
  //}
  for (it2 = DoCalibrationOuter.begin(); it2 != DoCalibrationOuter.end(); it2++) {
    for (it = (it2->second).begin(); it != (it2->second).end(); it++) {
      if (it->second) {
        InitializeRootHistogramsOuter_F(it2->first,it->first);
      }
    }
  }
}

/////////////////////////////////////////////////////////////////////
void TExogamPhysics::FillHistogramsCalib() {
  if (NPOptionManager::getInstance()->IsReader())
    m_EventData = &(**r_ReaderEventData);

  FillRootHistogramsCalib_F();
}

/////////////////////////////////////////////////////////////////////
void TExogamPhysics::InitializeRootHistogramsE_F(unsigned int DetectorNumber) {
  auto TH1Map = RootHistogramsCalib::getInstance()->GetTH1Map();

  TString hnameEXOE = Form("EXO_E%d", DetectorNumber);
  TString htitleEXOE = Form("EXO_E%d", DetectorNumber);
  (*TH1Map)["Exogam"][hnameEXOE] = new TH1F(hnameEXOE, htitleEXOE, 65536, 0, 65536);
}

/////////////////////////////////////////////////////////////////////
void TExogamPhysics::InitializeRootHistogramsOuter_F(unsigned int DetectorNumber, unsigned int OuterNumber) {
  auto TH1Map = RootHistogramsCalib::getInstance()->GetTH1Map();

  TString hnameEXOOuter = Form("EXO_Outer%d_%d", DetectorNumber, OuterNumber);
  TString htitleEXOOuter = Form("EXO_Outer%d_%d", DetectorNumber, OuterNumber);
  (*TH1Map)["Exogam"][hnameEXOOuter] = new TH1F(hnameEXOOuter, htitleEXOOuter, 65536, 0, 65536);
}

/////////////////////////////////////////////////////////////////////
void TExogamPhysics::InitializeRootHistogramsEHG_F(unsigned int DetectorNumber) {
  auto TH1Map = RootHistogramsCalib::getInstance()->GetTH1Map();

  TString hnameEXOEHG = Form("EXO_EHG%d", DetectorNumber);
  TString htitleEXOEHG = Form("EXO_EHG%d", DetectorNumber);
  (*TH1Map)["Exogam"][hnameEXOEHG] = new TH1F(hnameEXOEHG, htitleEXOEHG, 65536, 0, 65536);

}

/////////////////////////////////////////////////////////////////////
void TExogamPhysics::FillRootHistogramsCalib_F(){
  auto TH1Map = RootHistogramsCalib::getInstance()->GetTH1Map();
  TString hname;

  for (UShort_t i = 0; i < m_EventData->GetExoMult(); i++) {
    unsigned int DetectorNbr = m_EventData->GetExoCrystal(i);

    if(DoCalibrationE[DetectorNbr] && m_EventData->GetExoE(i) >0){
      hname = Form("EXO_E%d", DetectorNbr);
      (*TH1Map)["Exogam"][hname]->Fill(m_EventData->GetExoE(i));
    }
    if(DoCalibrationEHG[DetectorNbr] && m_EventData->GetExoEHG(i) >0){
      hname = Form("EXO_EHG%d", DetectorNbr);
      (*TH1Map)["Exogam"][hname]->Fill(m_EventData->GetExoEHG(i));
    }
    if(DoCalibrationOuter[DetectorNbr][0] && m_EventData->GetExoOuter1(i) >0){
      hname = Form("EXO_Outer%d_0", DetectorNbr);
      (*TH1Map)["Exogam"][hname]->Fill(m_EventData->GetExoOuter1(i));
    }
    if(DoCalibrationOuter[DetectorNbr][1] && m_EventData->GetExoOuter2(i) >0){
      hname = Form("EXO_Outer%d_1", DetectorNbr);
      (*TH1Map)["Exogam"][hname]->Fill(m_EventData->GetExoOuter2(i));
    }
    if(DoCalibrationOuter[DetectorNbr][2] && m_EventData->GetExoOuter3(i) >0){
      hname = Form("EXO_Outer%d_2", DetectorNbr);
      (*TH1Map)["Exogam"][hname]->Fill(m_EventData->GetExoOuter3(i));
    }
    if(DoCalibrationOuter[DetectorNbr][3] && m_EventData->GetExoOuter4(i) >0){
      hname = Form("EXO_Outer%d_3", DetectorNbr);
      (*TH1Map)["Exogam"][hname]->Fill(m_EventData->GetExoOuter4(i));
    }
  }
}

/////////////////////////////////////////////////////////////////////
void TExogamPhysics::DoCalibration() {
  std::cout << "Do Calibration Exogam" << std::endl;
  DefineCalibrationSource(Source_name);
  map<int, bool>::iterator it;
  map<int, map<int,bool>>::iterator it2;

  std::string Path = NPOptionManager::getInstance()->GetCalibrationOutputPath();
  std::string OutputName = NPOptionManager::getInstance()->GetOutputFile();
  if (OutputName.size() > 5) {
    if (OutputName.substr(OutputName.size() - 5, OutputName.size()) == ".root") {
      OutputName = OutputName.substr(0, OutputName.size() - 5);
    }
  }
  std::string make_folder = "mkdir " + Path + OutputName;

  MakeInitialCalibFolder(make_folder);

  ofstream* calib_file = new ofstream;
  ofstream* dispersion_file = new ofstream;

  if(!DoCalibrationE.empty()){
    MakeECalibFolders(make_folder);
    CreateCalibrationEFiles(calib_file, dispersion_file);
  }
  for (it = DoCalibrationE.begin(); it != DoCalibrationE.end(); it++) {
    if (it->second) {
      DoCalibrationE_F(it->first,"E", calib_file, dispersion_file, Threshold_E_Cal);
    }
  }
  calib_file->close();
  dispersion_file->close();

  if(!DoCalibrationEHG.empty()){
    MakeEHGCalibFolders(make_folder);
    CreateCalibrationEHGFiles(calib_file, dispersion_file);
  }
  for (it = DoCalibrationEHG.begin(); it != DoCalibrationEHG.end(); it++) {
    if (it->second) {
      DoCalibrationE_F(it->first,"EHG", calib_file, dispersion_file, Threshold_EHG_Cal);
    }
  }
  calib_file->close();
  dispersion_file->close();

  if(!DoCalibrationOuter.empty()){
    MakeOuterCalibFolders(make_folder);
    CreateCalibrationOuterFiles(calib_file, dispersion_file);
  }
  for (it2 = DoCalibrationOuter.begin(); it2 != DoCalibrationOuter.end(); it2++) {
    for (it = (it2->second).begin(); it != (it2->second).end(); it++) {
      if (it->second) {
        DoCalibrationE_F(it->first,Form("Outer%d_",it2->first), calib_file, dispersion_file, Threshold_Outers_Cal);
      }
    }
  }
  calib_file->close();
  dispersion_file->close();
}

/////////////////////////////////////////////////////////////////////
void TExogamPhysics::MakeInitialCalibFolder(std::string make_folder) {
  int sys = system(make_folder.c_str());
}

/////////////////////////////////////////////////////////////////////
void TExogamPhysics::MakeECalibFolders(std::string make_folder) {
  int sys =system((make_folder+"/Exogam_E").c_str()); 
}

/////////////////////////////////////////////////////////////////////
void TExogamPhysics::MakeEHGCalibFolders(std::string make_folder) {
  int sys =system((make_folder+"/Exogam_EHG").c_str());
}

/////////////////////////////////////////////////////////////////////
void TExogamPhysics::MakeOuterCalibFolders(std::string make_folder) {
  int sys =system((make_folder+"/Exogam_Outer").c_str());
}

/////////////////////////////////////////////////////////////////////
void TExogamPhysics::DoCalibrationE_F(unsigned int DetectorNumber,std::string CalibType, ofstream* calib_file, ofstream* dispersion_file, unsigned int Threshold) {
  auto TH1Map = RootHistogramsCalib::getInstance()->GetTH1Map();
  auto TGraphMap = RootHistogramsCalib::getInstance()->GetTGraphMap();

#if CUBIX
  CubixEnergyCal->Reset();
  std::string hnameEXOE = Form("EXO_%s%d",CalibType.c_str(), DetectorNumber);
  std::string htitleEXOE = Form("EXO_%s%d",CalibType.c_str(), DetectorNumber);

  auto hist = ((*TH1Map)["Exogam"][hnameEXOE]);

  CubixEnergyCal->SetDataFromHistTH1(hist,0);

  for (auto ie : Source_E)
    CubixEnergyCal->AddPeak(ie);

  CubixEnergyCal->SetGain(1.);
  CubixEnergyCal->SetVerbosityLevel(1);

  CubixEnergyCal->SetFitPlynomialOrder(FitPolOrder);
  CubixEnergyCal->SetNoOffset(false);
  CubixEnergyCal->UseLeftTail(true);
  CubixEnergyCal->UseRightTail(true);

  CubixEnergyCal->UseFirstDerivativeSearch();

  CubixEnergyCal->SetGlobalChannelLimits(hist->GetXaxis()->GetBinLowEdge(1)+Threshold,hist->GetXaxis()->GetBinLowEdge(hist->GetXaxis()->GetNbins()));      // limit the search to this range in channels
  CubixEnergyCal->SetGlobalPeaksLimits(15,5);   // default fwhm and minmum amplitude for the peaksearch [15 5]

  CubixEnergyCal->StartCalib();
  vector < Fitted > FitResults = CubixEnergyCal->GetFitResults();



  std:: cout << calib_file << " " << (*calib_file).is_open() << std::endl;
  std:: cout << hnameEXOE << " ";
  (*calib_file) << hnameEXOE << " ";
  if(FitResults.size() > 1)
  {
    for(unsigned int i = 0; i <= FitPolOrder; i++){
      (*calib_file) << scientific << setprecision(6) << setw(14) << CubixEnergyCal->fCalibFunction->GetParameter(i) << " ";
      std::cout << scientific << setprecision(6) << setw(14) << CubixEnergyCal->fCalibFunction->GetParameter(i) << " ";
    }
  }
  else
  {
    for(unsigned int i = 0; i <= FitPolOrder; i++){
      (*calib_file) << scientific << setprecision(6) << setw(14) << 0. << " ";
      std::cout << scientific << setprecision(6) << setw(14) << 0. << " ";
    }
  }
  (*calib_file) << "\n";
  std::cout << "\n";


  if(FitResults.size()>1 && CubixEnergyCal->fCalibFunction) {
    auto c = new TCanvas;
    c->SetName("CalibrationResults");
    c->SetTitle("Calibration Results");
    c->Divide(1,2,0.0001,0.0001);
    c->cd(1);
    CubixEnergyCal->fCalibGraph->Draw("ap");
    CubixEnergyCal->fCalibFunction->Draw("same");
    c->cd(2);
    CubixEnergyCal->fResidueGraph->Draw("ape");
    c->Update();
    c->Modified();
  }
  if(FitResults.size() > 1)
  {
    (*TGraphMap)["Exogam"][Form("Calib_Graph_%s",hnameEXOE.c_str())] = (TGraphErrors*)(CubixEnergyCal->fCalibGraph->Clone());
    (*TGraphMap)["Exogam"][Form("Calib_Graph_%s",hnameEXOE.c_str())]->GetYaxis()->SetTitle("Energy (MeV)");
    (*TGraphMap)["Exogam"][Form("Calib_Graph_%s",hnameEXOE.c_str())]->SetTitle(Form("Calibration_Graph_%s",hnameEXOE.c_str()));

    (*TGraphMap)["Exogam"][Form("Residue_Graph_%s",hnameEXOE.c_str())] = (TGraphErrors*)(CubixEnergyCal->fResidueGraph->Clone());
    (*TGraphMap)["Exogam"][Form("Residue_Graph_%s",hnameEXOE.c_str())]->GetXaxis()->SetTitle("Energy (MeV)");
    (*TGraphMap)["Exogam"][Form("Residue_Graph_%s",hnameEXOE.c_str())]->GetYaxis()->SetTitle("Residue (MeV)");
    (*TGraphMap)["Exogam"][Form("Residue_Graph_%s",hnameEXOE.c_str())]->SetTitle(Form("Residue_Graph_%s",hnameEXOE.c_str()));
  }
#else
  std::cout << "Exogam calibration currently not supported without CUBIX. Download CUBIX and set -DCUBIX=1 to use EXOGAM calibration\n";
  exit(1);

#endif

}

/////////////////////////////////////////////////////////////////////
void TExogamPhysics::DefineCalibrationSource(std::string source) {
  // 239Pu
  if(source == "60Co"){
    Source_isotope.push_back("$^{60}$Co");
    Source_E.push_back(1.17322);
    Source_Sig.push_back(0.0001);
    Source_branching_ratio.push_back(99.85);

    Source_isotope.push_back("$^{60}$Co");
    Source_E.push_back(1.33249);
    Source_Sig.push_back(0.0001);
    Source_branching_ratio.push_back(99.98);
  }
  else if(source == "152Eu"){
    Source_isotope.push_back("$^{152}$Eu");
    Source_E.push_back(0.121782);
    Source_Sig.push_back(0.0001);
    Source_branching_ratio.push_back(28.58);

    Source_isotope.push_back("$^{152}$Eu");
    Source_E.push_back(0.344279);
    Source_Sig.push_back(0.0001);
    Source_branching_ratio.push_back(26.5);

    Source_isotope.push_back("$^{152}$Eu");
    Source_E.push_back(1.40801);
    Source_Sig.push_back(0.0001);
    Source_branching_ratio.push_back(21.0);