TICPhysics.cxx

/*****************************************************************************
 * Copyright (C) 2009-2020   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: P. Morfouace  contact address: pierre.morfouace@cea.fr   *
 *                                                                           *
 * Creation Date  : October 2023                                             *
 * Last update    :                                                          *
 *---------------------------------------------------------------------------*
 * Decription:                                                               *
 *  This class hold IC Treated  data                                    *
 *                                                                           *
 *---------------------------------------------------------------------------*
 * Comment:                                                                  *
 *                                                                           *   
 *                                                                           *
 *****************************************************************************/

#include "TICPhysics.h"

//   STL
#include <sstream>
#include <iostream>
#include <cmath>
#include <stdlib.h>
#include <limits>
using namespace std;

//   NPL
#include "RootInput.h"
#include "RootOutput.h"
#include "NPDetectorFactory.h"
#include "NPOptionManager.h"

//   ROOT
#include "TChain.h"

ClassImp(TICPhysics);

///////////////////////////////////////////////////////////////////////////
TICPhysics::TICPhysics()
  : m_EventData(new TICData),
  m_PreTreatedData(new TICData),
  m_EventPhysics(this),
  m_FPMW_Section(-1),
  m_number_of_splines(34),
  m_Eres_Threshold(3000),
  m_Z_SPLINE_CORRECTION(false),
  m_NumberOfDetectors(0){
  }

///////////////////////////////////////////////////////////////////////////
/// A usefull method to bundle all operation to add a detector
void TICPhysics::AddDetector(TVector3 Pos){
  // In That simple case nothing is done
  // Typically for more complex detector one would calculate the relevant 
  // positions (stripped silicon) or angles (gamma array)
} 

///////////////////////////////////////////////////////////////////////////
void TICPhysics::AddDetector(double R, double Theta, double Phi){
  // Compute the TVector3 corresponding
  TVector3 Pos(R*sin(Theta)*cos(Phi),R*sin(Theta)*sin(Phi),R*cos(Theta));
  // Call the cartesian method
  AddDetector(Pos);
} 

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



///////////////////////////////////////////////////////////////////////////
void TICPhysics::BuildPhysicalEvent() {

  if(m_FPMW_Section<0)
    return;

  Clear();
  PreTreat();

  static CalibrationManager* Cal = CalibrationManager::getInstance();

  int size = m_PreTreatedData->GetICMult();
  for(int i=0; i<size; i++){
    int segment = m_PreTreatedData->GetIC_Section(i);
    double gain = Cal->GetValue("IC/SEC"+NPL::itoa(m_FPMW_Section)+"_SEG"+NPL::itoa(segment)+"_ALIGN",0);
    double GainInit = Cal->GetValue("IC/INIT_SEG"+NPL::itoa(segment)+"_ALIGN",0);

    fIC_raw[i] = m_PreTreatedData->GetIC_Charge(i);
    fIC[i] = gain*m_PreTreatedData->GetIC_Charge(i);
    fIC_Init[i] = GainInit * m_PreTreatedData->GetIC_Charge(i);

  }


  if (fIC_Init[1]>0 && fIC_Init[5]>0) {
    EtotInit = 0 ;
    double ScalorInit = Cal->GetValue("IC/INIT_ETOT_SCALING",0);

    for (int Seg = 0; Seg<10; Seg++){

      if (Seg == 0 && m_Data_Year == 2024 ){
        EtotInit += 0.75 * Cal->GetValue("IC/INIT_SEG"+NPL::itoa(Seg+1)+"_ALIGN",0)* fIC_raw[1];
      }

      else { 
        EtotInit += fIC_Init[Seg];  
      }
    } 
    EtotInit = EtotInit * ScalorInit ;
  } //Condition for Init Etot

  else {
    EtotInit = -100 ;
  }

  if(fIC[1]>0 && fIC[5]>0){
    DE = 0.5*(fIC_raw[0] + fIC_raw[1] + fIC_raw[2] + fIC_raw[3]) + fIC_raw[4];
    Eres = fIC_raw[5] + fIC_raw[6] + fIC_raw[7] + fIC_raw[8] + fIC_raw[9];

    double scalor = Cal->GetValue("IC/ETOT_SCALING_SEC"+NPL::itoa(m_FPMW_Section),0);

    for(int i=0; i<10; i++){
      Etot += fIC[i];
    }
    Etot = scalor*Etot;

    //Etot = 0.02411*(0.8686*fIC_raw[0]+0.7199*fIC_raw[1]+0.6233*fIC_raw[2]+0.4697*fIC_raw[3]+0.9787*fIC_raw[4]+0.9892*fIC_raw[5]+2.1038*fIC_raw[6]+1.9429*fIC_raw[7]+1.754*fIC_raw[8]+2.5*fIC_raw[9]);  


    if(m_Z_SPLINE_CORRECTION && Eres>m_Eres_Threshold){
      Chio_Z = Cal->ApplyCalibration("IC/Z_CALIBRATION",ApplyZSpline());
    }
    else Chio_Z = -1000;

  }
  else{
    DE = -100;
    Eres = -100;
    Etot = -100;
    Chio_Z = -1000;
  }

  m_FPMW_Section = -1;
}

///////////////////////////////////////////////////////////////////////////
void TICPhysics::PreTreat() {
  // This method typically applies thresholds and calibrations
  // Might test for disabled channels for more complex detector

  // clear pre-treated object
  ClearPreTreatedData();

  // instantiate CalibrationManager
  //static CalibrationManager* Cal = CalibrationManager::getInstance();

  unsigned int mysize = m_EventData->GetICMult();
  for (unsigned int i = 0; i < mysize ; ++i) {
    int segment = m_EventData->GetIC_Section(i);
    //cout << section << " " << gain << endl;
    //double charge = gain*m_EventData->GetIC_Charge(i);
    double charge = m_EventData->GetIC_Charge(i);

    m_PreTreatedData->SetIC_Charge(charge);
    m_PreTreatedData->SetIC_Section(segment);
  }
}



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

  // path to file
  string FileName = "./configs/ConfigIC.dat";

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

  if (!AnalysisConfigFile.is_open()) {
    cout << " No ConfigIC.dat found: Default parameter loaded for Analayis " << FileName << endl;
    return;
  }
  cout << " Loading user parameter for Analysis from ConfigIC.dat " << endl;

  // Save it in a TAsciiFile
  TAsciiFile* asciiConfig = RootOutput::getInstance()->GetAsciiFileAnalysisConfig();
  asciiConfig->AppendLine("%%% ConfigIC.dat %%%");
  asciiConfig->Append(FileName.c_str());
  asciiConfig->AppendLine("");
  // read analysis config file
  string LineBuffer,DataBuffer,whatToDo;
  while (!AnalysisConfigFile.eof()) {
    // Pick-up next line
    getline(AnalysisConfigFile, LineBuffer);

    // search for "header"
    string name = "ConfigIC";
    if (LineBuffer.compare(0, name.length(), name) == 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=="ERESIDUAL_THRESHOLD") {
        AnalysisConfigFile >> DataBuffer;
        m_Eres_Threshold = atof(DataBuffer.c_str());
        cout << whatToDo << " " << m_Eres_Threshold << endl;
      }

      else if (whatToDo=="DATA_YEAR") {
        AnalysisConfigFile >> DataBuffer;
        m_Data_Year = atof(DataBuffer.c_str());
        cout << whatToDo << " " << m_Data_Year << endl;
      }

      else if(whatToDo=="NUMBER_OF_SPLINES"){
        AnalysisConfigFile >> DataBuffer;
        m_number_of_splines = atoi(DataBuffer.c_str());
        cout << whatToDo << " " << m_number_of_splines << endl;
      }
      else if (whatToDo=="LOAD_Z_SPLINE"){
        AnalysisConfigFile >> DataBuffer;
        m_Z_SPLINE_PATH = DataBuffer;
        cout << "*** Loading Z spline ***" << endl;
        LoadZSpline();
      }
      else {
        ReadingStatus = false;
      }
    }
  }
}


///////////////////////////////////////////////////////////////////////////
void TICPhysics::LoadZSpline(){
  TString filename = m_Z_SPLINE_PATH;
  TFile* ifile = new TFile(filename,"read");

  if(ifile->IsOpen()){
    m_Z_SPLINE_CORRECTION = true;
    for(int i=0; i<m_number_of_splines; i++){
      m_Zspline[i] = (TSpline3*) ifile->FindObjectAny(Form("fspline_%d",i+1));
      cout << "*** " << m_Zspline[i]->GetName() << " is loaded!" << endl;
    }

  }
  else
    cout << "File " << filename << " not found!" << endl;
  ifile->Close();
}

///////////////////////////////////////////////////////////////////////////
double TICPhysics::ApplyZSpline(){
  double DEcorr;
  double FF_DEcorr0[m_number_of_splines];
  for(int i=0; i<m_number_of_splines; i++){
    FF_DEcorr0[i] = m_Zspline[i]->Eval(8500);
  }

  double DEspline0;
  double Eval_DEspline;
  int index=0;
  for(int i=0; i<m_number_of_splines; i++){
    Eval_DEspline = m_Zspline[i]->Eval(Eres);
    if(DE<Eval_DEspline) break;
    index = i;
  }

  Eval_DEspline = m_Zspline[index]->Eval(Eres);
  DEspline0 = FF_DEcorr0[index];
  double dmin, dsup;
  if(index<(m_number_of_splines-1) && DE>m_Zspline[0]->Eval(Eres)){
    dmin = DE - m_Zspline[index]->Eval(Eres);
    dsup = m_Zspline[index+1]->Eval(Eres) - DE;

    Eval_DEspline = dsup*m_Zspline[index]->Eval(Eres)/(dmin+dsup) + dmin*m_Zspline[index+1]->Eval(Eres)/(dmin+dsup);
    DEspline0     = dsup*FF_DEcorr0[index]/(dmin+dsup) + dmin*FF_DEcorr0[index+1]/(dmin+dsup);

    DEcorr        = DEspline0 * DE / Eval_DEspline;
  }
  else if(index==m_number_of_splines-1){
    Eval_DEspline = m_Zspline[index]->Eval(Eres);
    DEspline0     = FF_DEcorr0[index];

    DEcorr        = DEspline0 * DE / Eval_DEspline;
  }

  return DEcorr;
}

///////////////////////////////////////////////////////////////////////////
void TICPhysics::Clear() {
  DE = 0;
  Eres = 0;
  Etot = 0;
  EtotInit = 0;
  Chio_Z = 0;
  for(int i=0; i<11; i++){
    fIC[i] = 0;
    fIC_raw[i] = 0;
    fIC_Init[i] = 0;
  }
}



///////////////////////////////////////////////////////////////////////////
void TICPhysics::ReadConfiguration(NPL::InputParser parser) {
  vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("IC");
  if(NPOptionManager::getInstance()->GetVerboseLevel())
    cout << "//// " << blocks.size() << " detectors found " << endl; 

  vector<string> cart = {"POS"};
  vector<string> sphe = {"R","Theta","Phi"};

  for(unsigned int i = 0 ; i < blocks.size() ; i++){
    if(blocks[i]->HasTokenList(cart)){
      if(NPOptionManager::getInstance()->GetVerboseLevel())
        cout << endl << "////  IC " << i+1 <<  endl;

      TVector3 Pos = blocks[i]->GetTVector3("POS","mm");
      AddDetector(Pos);
    }
    else if(blocks[i]->HasTokenList(sphe)){
      if(NPOptionManager::getInstance()->GetVerboseLevel())
        cout << endl << "////  IC " << i+1 <<  endl;
      double R = blocks[i]->GetDouble("R","mm");
      double Theta = blocks[i]->GetDouble("Theta","deg");
      double Phi = blocks[i]->GetDouble("Phi","deg");
      AddDetector(R,Theta,Phi);
    }
    else{
      cout << "ERROR: check your input file formatting " << endl;
      exit(1);
    }
  }

  ReadAnalysisConfig();
}


///////////////////////////////////////////////////////////////////////////
void TICPhysics::AddParameterToCalibrationManager() {
  CalibrationManager* Cal = CalibrationManager::getInstance();

  //Calibration from online analysis
  Cal->AddParameter("IC","Z_CALIBRATION","IC_Z_CALIBRATION");
  Cal->AddParameter("IC","INIT_ETOT_SCALING","IC_INIT_ETOT_SCALING");
  for (int segment=0; segment<11;segment++){
    Cal->AddParameter("IC","INIT_SEG"+NPL::itoa(segment+1)+"_ALIGN","IC_INIT_SEG"+NPL::itoa(segment+1)+"_ALIGN");
  }


  for(int section = 0; section<20; section++){
    Cal->AddParameter("IC","ETOT_SCALING_SEC"+NPL::itoa(section),"IC_ETOT_SCALING_SEC"+NPL::itoa(section));
    for(int segment = 0; segment<11; segment++){
      Cal->AddParameter("IC","SEC"+NPL::itoa(section)+"_SEG"+NPL::itoa(segment+1)+"_ALIGN","IC_SEC"+NPL::itoa(section)+"_SEG"+NPL::itoa(segment+1)+"_ALIGN");
    }
  }
}

///////////////////////////////////////////////////////////////////////////
void TICPhysics::InitializeRootInputRaw() {
  TChain* inputChain = RootInput::getInstance()->GetChain();
  inputChain->SetBranchStatus("IC",  true );
  inputChain->SetBranchAddress("IC", &m_EventData );
}



///////////////////////////////////////////////////////////////////////////
void TICPhysics::InitializeRootInputPhysics() {
  TChain* inputChain = RootInput::getInstance()->GetChain();
  inputChain->SetBranchAddress("IC", &m_EventPhysics);
}



///////////////////////////////////////////////////////////////////////////
void TICPhysics::InitializeRootOutput() {
  TTree* outputTree = RootOutput::getInstance()->GetTree();
  outputTree->Branch("IC", "TICPhysics", &m_EventPhysics);
}



////////////////////////////////////////////////////////////////////////////////
//            Construct Method to be pass to the DetectorFactory              //
////////////////////////////////////////////////////////////////////////////////
NPL::VDetector* TICPhysics::Construct() {
  return (NPL::VDetector*) new TICPhysics();
}



////////////////////////////////////////////////////////////////////////////////
//            Registering the construct method to the factory                 //
////////////////////////////////////////////////////////////////////////////////
extern "C"{
  class proxy_IC{
    public:
      proxy_IC(){
        NPL::DetectorFactory::getInstance()->AddToken("IC","IC");
        NPL::DetectorFactory::getInstance()->AddDetector("IC",TICPhysics::Construct);
      }
  };

  proxy_IC p_IC;
}