Skip to content

GitLab

Commit 6c45c0e8 authored by Valérian Alcindor's avatar Valérian Alcindor
Browse files

Adding the latest versions of EXL analysis (from experiment e744)

parent 3c4f7a20
Hide whitespace changes
Inline Side-by-side
Showing with 812 additions and 482 deletions
NPLib/Detectors/EXL/TExlData.cxx
View file @ 6c45c0e8
......@@ -22,39 +22,71 @@
#include <iostream>
#include "TExlData.h"
ClassImp(TExlData)
TExlData::TExlData()
{
}
TExlData::~TExlData()
{
}
void TExlData::Clear()
{
fExl_Energy.clear();
fExl_Time.clear();
fExl_Number.clear();
fExl_Crystal.clear();
}
fExlE_Energy.clear();
fExlE_Number.clear();
fExlE_Crystal.clear();
fExlT_Time.clear();
fExlT_Number.clear();
fExlT_Crystal.clear();
}
void TExlData::Dump() const
{
cout << "XXXXXXXXXXXXXXXXXXXXXXXX New Event XXXXXXXXXXXXXXXXX" << endl;
cout << "XXXXXXXXXXXXXXXXXXXXXXXX New Event XXXXXXXXXXXXXXXXX" << endl;
cout << "Energy:" << endl;
for(unsigned short i = 0 ; i<fExlE_Energy.size() ; i ++)
{
cout << "Exl Number " << fExlE_Number[i]
<< " Crystal Number " << fExlE_Crystal[i]
<< " Energy: " << fExlE_Energy[i] << endl;
}
cout << "Time:" << endl;
for(unsigned short i = 0 ; i<fExlT_Time.size() ; i ++)
{
cout << "Exl Number " << fExlT_Number[i]
<< " Crystal Number " << fExlT_Crystal[i]
<< " Time: " << fExlT_Time[i] << endl;
}
for(unsigned short i = 0 ; i<fExl_Energy.size() ; i ++)
{
cout << "Exl Number " << fExl_Number[i] << " Crystal Number " << fExl_Crystal[i] << " Energy: " << fExl_Energy[i] << " Time: " << fExl_Time[i] << endl;
}
}
double TExlData::GetEnergyCrystalNumber(const int crystal){
double energy = -1000;
unsigned int sizeE = GetEMult();
for (unsigned int i = 0; i < sizeE; i++){
if (GetE_CrystalNumber(i) == crystal){
energy = GetEnergy(i);
// return GetEnergy(i);
}
}
return energy;
}
double TExlData::GetTimeCrystalNumber(const int crystal){
double time = -1000;
unsigned int sizeT = GetTMult();
for (unsigned int i = 0; i < sizeT; i++){
if (GetT_CrystalNumber(i)==crystal){
time = GetTime(i);
// return GetTime(i);
}
}
return time;
}
}
NPLib/Detectors/EXL/TExlData.h
View file @ 6c45c0e8
......@@ -29,38 +29,58 @@ using namespace std ;
class TExlData : public TObject {
private:
// ADC
vector<double> fExl_Energy;
vector<double> fExl_Time;
vector<short> fExl_Number ;
vector<short> fExl_Crystal;
private:
// ADC
vector<double> fExlE_Energy;
vector<short> fExlE_Number ;
vector<short> fExlE_Crystal;
vector<double> fExlT_Time;
vector<short> fExlT_Number ;
vector<short> fExlT_Crystal;
public:
TExlData();
virtual ~TExlData();
void Clear();
void Clear(const Option_t*) {};
void Dump() const;
public:
TExlData();
virtual ~TExlData();
///////////////////// GETTERS ////////////////////////
inline double GetEnergy(const int& i) const { return fExl_Energy[i] ;}
inline double GetTime(const int& i) const { return fExl_Time[i] ;}
inline int GetExlNumber(const int& i) const { return fExl_Number[i] ;}
inline int GetCrystalNumber(const int& i) const { return fExl_Crystal[i] ;}
//Mult
inline double GetMult() const { return fExl_Energy.size() ;}
void Clear();
void Clear(const Option_t*) {};
void Dump() const;
///////////////////// SETTERS ////////////////////////
inline void SetEandTime(const int& N,const int& C,const double& E,const double& T){
fExl_Energy.push_back(E);
fExl_Time.push_back(T);
fExl_Number.push_back(N) ;
fExl_Crystal.push_back(C);
}
//////////////////////////////////////////////////////////////
// Getters and Setters
// Prefer inline declaration to avoid unnecessary called of
// frequently used methods
// add //! to avoid ROOT creating dictionnary for the methods
public:
///////////////////// GETTERS ////////////////////////
inline double GetEnergy(const int& i) const { return fExlE_Energy[i] ;}
double GetEnergyCrystalNumber(const int crystal);
double GetTimeCrystalNumber(const int crystal);
inline double GetTime(const int& i) const { return fExlT_Time[i] ;}
inline int GetE_ExlNumber(const int& i) const { return fExlE_Number[i] ;}
inline int GetE_CrystalNumber(const int& i) const { return fExlE_Crystal[i] ;}
inline int GetT_ExlNumber(const int& i) const { return fExlT_Number[i] ;}
inline int GetT_CrystalNumber(const int& i) const { return fExlT_Crystal[i] ;}
ClassDef(TExlData,1) // ExlData structure
//Mult
inline unsigned int GetEMult() const {return fExlE_Energy.size();}
inline unsigned int GetTMult() const {return fExlT_Time.size();}
///////////////////// SETTERS ////////////////////////
inline void SetEnergy(const int& N,const int& C,const double& E){
fExlE_Number.push_back(N) ;
fExlE_Crystal.push_back(C);
fExlE_Energy.push_back(E);
}
inline void SetTime(const int& N,const int& C,const double& T){
fExlT_Number.push_back(N) ;
fExlT_Crystal.push_back(C);
fExlT_Time.push_back(T);
}
ClassDef(TExlData,1) // ExlData structure
};
#endif
NPLib/Detectors/EXL/TExlPhysics.cxx
View file @ 6c45c0e8
......@@ -7,21 +7,24 @@
/*****************************************************************************
* Original Author: L. Lefebvre contact address: lefebvrl@ipno.in2p3.fr *
* Update : V. Girard-Alcindor contact address: girardalcindor@ganil.fr *
* *
* Creation Date : October 2011 *
* Last update : *
* Last update : 2018 *
*---------------------------------------------------------------------------*
* Decription: *
* This class hold the Exl Detector Physics *
* *
*---------------------------------------------------------------------------*
* Comment: *
* Comment: Due to the way the detector is build, one has to take into *
* account the cross-talk between two neighboring crystals. *
* *
* *
*****************************************************************************/
// NPL
#include "TExlPhysics.h"
#include "TExlData.h"
#include "RootOutput.h"
#include "RootInput.h"
#include "NPDetectorFactory.h"
......@@ -38,199 +41,320 @@ using namespace std;
#include "TChain.h"
// tranform an integer to a string
string itoa(int value)
{
char buffer [33];
sprintf(buffer,"%d",value);
return buffer;
string itoa(int value){
char buffer [33];
sprintf(buffer,"%d",value);
return buffer;
}
ClassImp(TExlPhysics)
///////////////////////////////////////////////////////////////////////////
TExlPhysics::TExlPhysics()
{
NumberOfCrystal = 0;
EventData = new TExlData ;
EventPhysics = this ;
// Raw Threshold
m_E_RAW_Threshold = 0;
}
ClassImp(TExlPhysics);
///////////////////////////////////////////////////////////////////////////
TExlPhysics::~TExlPhysics()
{}
TExlPhysics::TExlPhysics(){
m_NumberOfDetector=0;
m_NumberOfCrystal = 18;
m_EventData = new TExlData;
m_EventPhysics = this;
m_E_RAW_Threshold = 0;
///////////////////////////////////////////////////////////////////////////
void TExlPhysics::Clear()
{
CrystalNumber.clear() ;
EXL_Energy.clear() ;
EXL_Time=-1000 ;
}
A = 0.;
B = 0.;
a = 0.;
b = 0.;
v_A.clear();
v_B.clear();
v_a.clear();
v_b.clear();
m_Crystal_i = 0;
m_NeighborCrystal_i = 0;
m_Channel_i = 0.;
m_NeighborChannel_i = 0.;
m_Energy_i = 0.;
///////////////////////////////////////////////////////////////////////////
void TExlPhysics::ReadConfiguration(NPL::InputParser parser) {
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("EXL");
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << "//// " << blocks.size() << " detectors found " << endl;
vector<string> token= {"POS"};
vector<TVector3> Center_CsI_Crystals;
for(unsigned int i = 0 ; i < blocks.size() ; i++){
if(blocks[i]->HasTokenList(token)){
TVector3 pos = blocks[i]->GetTVector3("POS","mm");
Center_CsI_Crystals.push_back(pos);
}
// std::ifstream pre_calibration_coeff;
// pre_calibration_coeff.open("/home/alcindor/Documents/Projects/EXL/Analysis/Calibration/pre_calibration_coeff2.txt");
else{
cout << "ERROR: check your input file formatting " << endl;
exit(1);
// std::ifstream calibration_coeff;
// calibration_coeff.open("/home/alcindor/Documents/Projects/EXL/Analysis/Calibration/Co60_EXL_Calibration.txt");
// if(v_A.size()==0){
// while(pre_calibration_coeff >> B >> A){
// v_A.push_back(A);
// v_B.push_back(B);
// }
// }
// if(v_a.size()==0){
// while(calibration_coeff >> Channel_name >> b >> a){
// v_a.push_back(a);
// v_b.push_back(b);
// }
// }
for (int i = 0; i < 30 ; i++){
v_A.push_back(1);
v_B.push_back(1);
v_a.push_back(1);
v_b.push_back(1);
}
}
AddEXL(Center_CsI_Crystals) ;
}
///////////////////////////////////////////////////////////////////////////
void TExlPhysics::AddEXL(vector <TVector3> Center_CsI_Crystals)
{
// In progress
// Needs the detectors dimensions
// Crystals in this detector can be mapped differently for every module
// Rows_y (0,1,2) Cols_x (0,1,2,3,4,5)
//local zero is at row 1, col 2.5
double arbit_length = 100; // detector arbit. face length and width, 100 mm
double col_pitch = arbit_length/6;
double row_pitch = arbit_length/3;
for(unsigned i=0; i<Center_CsI_Crystals.size(); i++)
{
TVector3 direction = Center_CsI_Crystals.at(i).Unit();
//place all 18 crystals in place
for(unsigned iCry=0; iCry<18; iCry++)
{
int Row_y = iCry/6;
int Col_x = iCry%6;
// Define Detector position localy
TVector3 localPos( (Row_y-1)*row_pitch,(Col_x-2.5)*row_pitch,0 );
// Rotate
localPos.RotateUz(direction);
TVector3 globalPos = Center_CsI_Crystals.at(i)+localPos;
cout << globalPos.X() << " " << globalPos.Y() << " " << globalPos.Z() << endl;
CsIPosition.push_back(globalPos);
}
}
TExlPhysics::~TExlPhysics(){
}
///////////////////////////////////////////////////////////////////////////
TVector3 TExlPhysics::GetPositionOfInteraction(int det, int cry)
{
int N = det*18+cry; // hyper crystal number
TVector3 Position = TVector3 (CsIPosition.at(N).X(),CsIPosition.at(N).Y(),CsIPosition.at(N).Z()) ;
void TExlPhysics::Clear(){
CrystalNumber.clear() ;
Energy.clear() ;
RawEnergy.clear();
Time.clear() ;
AddBack.clear();
}
return(Position) ;
///////////////////////////////////////////////////////////////////////////
void TExlPhysics::ReadConfiguration(NPL::InputParser parser){
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("EXL");
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << "//// " << blocks.size() << " detectors found " << endl;
}
vector<string> token = {"POS"};
vector<TVector3> Center_CsI_Crystals;
for(unsigned int i = 0; i < blocks.size(); i++){
if(blocks[i]->HasTokenList(token)){
TVector3 pos = blocks[i]->GetTVector3("POS","mm");
Center_CsI_Crystals.push_back(pos);
}
else{
cout << "ERROR: check your input file formatting " << endl;
exit(1);
}
}
AddEXL(Center_CsI_Crystals) ;
}
///////////////////////////////////////////////////////////////////////////
void TExlPhysics::AddParameterToCalibrationManager()
{
CalibrationManager* Cal = CalibrationManager::getInstance();
void TExlPhysics::AddEXL(vector <TVector3> Center_CsI_Crystals){
m_NumberOfDetector++;
// In progress
// Needs the detectors dimensions
// Crystals in this detector can be mapped differently for every module
// Rows_y (0,1,2) Cols_x (0,1,2,3,4,5)
// local zero is at row 1, col 2.5
double length = 110; // detector arbit. face length and width, 100 mm
double width = 7.2;
double height = 9.8;
for(unsigned int i = 0 ; i < NumberOfCrystal ; i++)
{
Cal->AddParameter("EXL", "_E_"+ NPL::itoa(i+1),"EXL_E_"+ NPL::itoa(i+1)) ;
Cal->AddParameter("EXL", "_T_"+ NPL::itoa(i+1),"EXL_T_"+ NPL::itoa(i+1)) ;
}
}
for(unsigned i=0; i<Center_CsI_Crystals.size(); i++){
TVector3 direction = Center_CsI_Crystals.at(i).Unit();
//place all 18 crystals in place
for(unsigned iCry=0; iCry<18; iCry++){
int Row_y = iCry/6;
int Col_x = iCry%6;
// Define Detector position localy
TVector3 localPos((Row_y-1)*length,(Col_x-2.5)*length,0);
// Rotate
localPos.RotateUz(direction);
TVector3 globalPos = Center_CsI_Crystals.at(i)+localPos;
CsIPosition.push_back(globalPos);
}
}
}
///////////////////////////////////////////////////////////////////////////
void TExlPhysics::InitializeRootInputRaw()
{
TChain* inputChain = RootInput::getInstance()->GetChain() ;
inputChain->SetBranchStatus ( "Exl" , true ) ;
inputChain->SetBranchStatus ( "fExl_*" , true ) ;
inputChain->SetBranchAddress( "Exl" , &EventData ) ;
}
///////////////////////////////////////////////////////////////////////////
void TExlPhysics::InitializeRootInputPhysics()
{
TChain* inputChain = RootInput::getInstance()->GetChain();
inputChain->SetBranchStatus ( "Exl", true );
inputChain->SetBranchStatus ( "CrystalNumber", true );
inputChain->SetBranchStatus ( "EXL_Energy", true );
inputChain->SetBranchStatus ( "EXL_Time", true );
inputChain->SetBranchAddress( "Exl", &EventPhysics );
}
TVector3 TExlPhysics::GetPositionOfInteraction(int det, int cry){
int N = det*18+cry; // hyper crystal number
TVector3 Position = TVector3(CsIPosition.at(N).X(),
CsIPosition.at(N).Y(),
CsIPosition.at(N).Z());
return(Position) ;
}
///////////////////////////////////////////////////////////////////////////
void TExlPhysics::InitializeRootOutput()
{
TTree* outputTree = RootOutput::getInstance()->GetTree() ;
outputTree->Branch( "Exl" , "TExlPhysics" , &EventPhysics ) ;
}
void TExlPhysics::AddParameterToCalibrationManager(){
CalibrationManager* Cal = CalibrationManager::getInstance();
for(unsigned int i = 0 ; i < m_NumberOfCrystal ; i++){
Cal->AddParameter("EXL", "E_"+ NPL::itoa(i+1), "EXL_E_" + NPL::itoa(i+1));
Cal->AddParameter("EXL", "T_"+ NPL::itoa(i+1), "EXL_T_" + NPL::itoa(i+1));
}
}
///////////////////////////////////////////////////////////////////////////
void TExlPhysics::BuildPhysicalEvent()
{
BuildSimplePhysicalEvent() ;
}
void TExlPhysics::InitializeRootInputRaw(){
TChain* inputChain = RootInput::getInstance()->GetChain();
inputChain->SetBranchStatus("EXL", true);
inputChain->SetBranchStatus("fExl_*", true);
inputChain->SetBranchAddress("EXL", &m_EventData );
}
///////////////////////////////////////////////////////////////////////////
void TExlPhysics::BuildSimplePhysicalEvent(){
if(EventData->GetMult()>0){
for(unsigned int i = 0 ; i < EventData->GetMult() ; i++){
if(EventData->GetEnergy(i)>m_E_RAW_Threshold){
CrystalNumber.push_back( EventData->GetCrystalNumber(i)) ;
EXL_Energy.push_back( CalibrationManager::getInstance()->ApplyCalibration("EXL/_E_" + NPL::itoa( EventData->GetExlNumber(i) ),EventData->GetEnergy(i) ) );
}
}
}
else{
CrystalNumber.push_back(-1000) ;
EXL_Energy.push_back(-1000);
}
if(EventData->GetMult()>0){
for(unsigned int i = 0 ; i < EventData->GetMult() ; i++){
//EXL_Time.push_back( CalibrationManager::getInstance()->ApplyCalibration("EXL/_T_" + NPL::itoa( EventData->GetExlNumber(i) ),EventData->GetTime(i) ) );
EXL_Time=EventData->GetTime(i);
}
}
else{
EXL_Time=-1000;
}
void TExlPhysics::InitializeRootInputPhysics(){
TChain* inputChain = RootInput::getInstance()->GetChain();
inputChain->SetBranchStatus("EXL", true);
inputChain->SetBranchStatus("CrystalNumber", true);
inputChain->SetBranchStatus("EXL_Energy", true);
inputChain->SetBranchStatus("EXL_AddBack", true );
inputChain->SetBranchStatus("EXL_Time", true );
inputChain->SetBranchAddress("EXL", &m_EventPhysics);
inputChain->SetBranchAddress("EXL", &m_EventData );
}
///////////////////////////////////////////////////////////////////////////
void TExlPhysics::InitializeRootOutput(){
TTree* outputTree = RootOutput::getInstance()->GetTree();
outputTree->Branch("EXL", "TExlPhysics", &m_EventPhysics);
}
///////////////////////////////////////////////////////////////////////////
double TExlPhysics::DopplerCorrection(double E, double Theta, double beta)
{
double Pi = 3.141592654 ;
void TExlPhysics::BuildPhysicalEvent(){
double E_corr = 0;
double gamma = 1./sqrt(1-beta*beta);
unsigned int sizeE = m_EventData->GetEMult();
for(unsigned int i = 0 ; i < sizeE ; i++){
E_corr = gamma*E*(1.-beta*cos(Theta*Pi/180.));
CrystalNumber.push_back(m_EventData->GetE_CrystalNumber(i));
return(E_corr);
double energy = CalibrationManager::getInstance()->ApplyCalibration(
"EXL/E_" + NPL::itoa(m_EventData->GetE_CrystalNumber(i)),
m_EventData->GetEnergy(i));
}
Energy.push_back(energy);
////////////////////////////////////////////////////////////////////////////////
// Construct Method to be pass to the DetectorFactory //
////////////////////////////////////////////////////////////////////////////////
NPL::VDetector* TExlPhysics::Construct(){
return (NPL::VDetector*) new TExlPhysics();
}
// look for associated time
bool check = false;
unsigned int sizeT = m_EventData->GetTMult();
for(unsigned int j = 0 ; j < sizeT ; j++){
if(m_EventData->GetT_CrystalNumber(j)==m_EventData->GetE_CrystalNumber(i)){
double time = CalibrationManager::getInstance()->ApplyCalibration(
"EXL/T_" + NPL::itoa(m_EventData->GetT_CrystalNumber(j)),
m_EventData->GetTime(j));
// cout << time << endl;
////////////////////////////////////////////////////////////////////////////////
// Registering the construct method to the factory //
////////////////////////////////////////////////////////////////////////////////
extern "C"{
class proxy_exl{
public:
proxy_exl(){
NPL::DetectorFactory::getInstance()->AddToken("EXL","Exl");
NPL::DetectorFactory::getInstance()->AddDetector("EXL",TExlPhysics::Construct);
Time.push_back(time);
check = true;
break;
}
}
};
if(!check)
Time.push_back(-1000);
proxy_exl p;
}
}
// for(int m_Crystal_i = 1; m_Crystal_i < 19; m_Crystal_i++){
// // Getting NeighborCrystalNumber
// if(m_Crystal_i % 2 == 0){
// m_NeighborCrystal_i = m_Crystal_i - 1;
// }
// else{
// m_NeighborCrystal_i = m_Crystal_i + 1;
// }
// // Getting channel of both crystals
// m_Channel_i = m_EventData->GetEnergyCrystalNumber(m_Crystal_i);
// m_NeighborChannel_i = m_EventData->GetEnergyCrystalNumber(m_NeighborCrystal_i);
// // Getting time
// Time.push_back(m_EventData->GetTimeCrystalNumber(m_Crystal_i));
// if (m_Channel_i != -10){
// if(m_NeighborChannel_i > 1){
// m_Channel_i = m_Channel_i
// - v_A.at(m_NeighborCrystal_i - 1) * m_NeighborChannel_i
// - v_B.at(m_NeighborCrystal_i - 1);
// }
// if (m_Channel_i < 0){
// m_Channel_i = 0.;
// }
// m_Energy_i = v_a.at(m_Crystal_i - 1) * m_Channel_i
// + v_b.at(m_Crystal_i - 1);
// CrystalNumber.push_back(m_Crystal_i);
// Energy.push_back(m_Energy_i);
// }
// }
// // AddBack