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);
}