Commit 63feceff authored by Adrien Matta's avatar Adrien Matta
Browse files

* progress on Hodo and Samurai analysis

parent 6b4f8fd9
......@@ -248,15 +248,17 @@ void TSamuraiFDC2Physics::PreTreat(){
}
// a valid wire must have an edge
if(etime && time && etime-time>ToTThreshold_L && etime-time<ToTThreshold_H){
Detector.push_back(det);
Layer.push_back(layer);
Wire.push_back(wire);
Time.push_back(time);
ToT.push_back(etime-time);
channel="SamuraiFDC2/L" + NPL::itoa(layer);
// rescalling is needed because calib are bad.
// to be fixed
DriftLength.push_back(10-Cal->ApplySigmoid(channel,etime));
if(!(wire==93 && layer ==7)){// remove noisy wire
Detector.push_back(det);
Layer.push_back(layer);
Wire.push_back(wire);
Time.push_back(time);
ToT.push_back(etime-time);
channel="SamuraiFDC2/L" + NPL::itoa(layer);
// rescalling is needed because calib are bad.
// to be fixed
DriftLength.push_back(10-Cal->ApplySigmoid(channel,etime));
}
}
}
}
......@@ -302,6 +304,24 @@ void TSamuraiFDC2Physics::RemoveNoise(){
}
return;
}
////////////////////////////////////////////////////////////////////////////////
TVector3 TSamuraiFDC2Physics::ProjectedPosition(double Z){
TVector3 pos(-10000,-10000,-10000);
if(PosX!=-10000){
pos = TVector3(PosX,PosY,0)+(Z/Dir.Z())*Dir;
//cout << pos.X() << " " << pos.Y() << " " << pos.Z() << endl;
}
return pos;
}
////////////////////////////////////////////////////////////////////////////////
double TSamuraiFDC2Physics::ProjectedPositionX(double Z){
return ProjectedPosition(Z).X();
}
////////////////////////////////////////////////////////////////////////////////
double TSamuraiFDC2Physics::ProjectedPositionY(double Z){
return ProjectedPosition(Z).Y();
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC2Physics::Clear(){
MultMean=0;
......
......@@ -83,7 +83,8 @@ class TSamuraiFDC2Physics : public TObject, public NPL::VDetector{
public:
// Projected position at given Z plan
TVector3 ProjectedPosition(double Z);
double ProjectedPositionX(double Z);
double ProjectedPositionY(double Z);
private: // Charateristic of the DC
void AddDC(std::string name, NPL::XmlParser&);//! take the XML file and fill in Wire_X and Layer_Angle
std::map<SamuraiDCIndex,double> Wire_X;//! X position of the wires
......
......@@ -77,47 +77,61 @@ class TSamuraiHodoscopeData : public TObject {
////////////////////// SETTERS ////////////////////////
// UP //
// Charge
inline void SetChargeUp(const Double_t& ID, const Double_t& Charge){
inline void SetChargeUp(const UShort_t& ID, const Double_t& Charge){
fSamuraiHodoscope_Qu_ID.push_back(ID);
fSamuraiHodoscope_Qu_Charge.push_back(Charge);
};//!
// Time
inline void SetTimeUp(const Double_t& ID, const Double_t& Time){
inline void SetTimeUp(const UShort_t& ID, const Double_t& Time){
fSamuraiHodoscope_Tu_ID.push_back(ID);
fSamuraiHodoscope_Tu_Time.push_back(Time);
};//!
// DOWN //
// Charge
inline void SetChargeDown(const Double_t& ID, const Double_t& Charge){
inline void SetChargeDown(const UShort_t& ID, const Double_t& Charge){
fSamuraiHodoscope_Qd_ID.push_back(ID);
fSamuraiHodoscope_Qd_Charge.push_back(Charge);
};//!
// Time
inline void SetTimeDown(const Double_t& ID, const Double_t& Time){
inline void SetTimeDown(const UShort_t& ID, const Double_t& Time){
fSamuraiHodoscope_Td_ID.push_back(ID);
fSamuraiHodoscope_Td_Time.push_back(Time);
};//!
/*
////////////////////// GETTERS ////////////////////////
// Energy
inline UShort_t GetMultEnergy() const
{return fSamuraiHodoscope_E_DetectorNbr.size();}
inline UShort_t GetE_DetectorNbr(const unsigned int &i) const
{return fSamuraiHodoscope_E_DetectorNbr[i];}//!
inline Double_t Get_Energy(const unsigned int &i) const
{return fSamuraiHodoscope_Energy[i];}//!
inline UShort_t GetMultQUp() const
{return fSamuraiHodoscope_Qu_ID.size();}
inline UShort_t GetQUp_ID(const unsigned int &i) const
{return fSamuraiHodoscope_Qu_ID[i];}//!
inline UShort_t GetQUp_Charge(const unsigned int &i) const
{return fSamuraiHodoscope_Qu_Charge[i];}//!
// Time
inline UShort_t GetMultTUp() const
{return fSamuraiHodoscope_Tu_ID.size();}
inline UShort_t GetTUp_ID(const unsigned int &i) const
{return fSamuraiHodoscope_Tu_ID[i];}//!
inline Double_t GetTUp_Time(const unsigned int &i) const
{return fSamuraiHodoscope_Tu_Time[i];}//!
// Energy
inline UShort_t GetMultQDown() const
{return fSamuraiHodoscope_Qd_ID.size();}
inline UShort_t GetQDown_ID(const unsigned int &i) const
{return fSamuraiHodoscope_Qd_ID[i];}//!
inline Double_t GetQDown_Charge(const unsigned int &i) const
{return fSamuraiHodoscope_Qd_Charge[i];}//!
// Time
inline UShort_t GetMultTime() const
{return fSamuraiHodoscope_T_DetectorNbr.size();}
inline UShort_t GetT_DetectorNbr(const unsigned int &i) const
{return fSamuraiHodoscope_T_DetectorNbr[i];}//!
inline Double_t Get_Time(const unsigned int &i) const
{return fSamuraiHodoscope_Time[i];}//!
*/
inline UShort_t GetMultTDown() const
{return fSamuraiHodoscope_Td_ID.size();}
inline UShort_t GetTDown_ID(const unsigned int &i) const
{return fSamuraiHodoscope_Td_ID[i];}//!
inline Double_t GetTDown_Time(const unsigned int &i) const
{return fSamuraiHodoscope_Td_Time[i];}//!
//////////////////////////////////////////////////////////////
// Required for ROOT dictionnary
......
......@@ -42,169 +42,152 @@ using namespace std;
ClassImp(TSamuraiHodoscopePhysics)
///////////////////////////////////////////////////////////////////////////
TSamuraiHodoscopePhysics::TSamuraiHodoscopePhysics()
: m_EventData(new TSamuraiHodoscopeData),
m_PreTreatedData(new TSamuraiHodoscopeData),
m_EventPhysics(this),
//m_Spectra(0),
m_E_RAW_Threshold(0), // adc channels
m_E_Threshold(0), // MeV
m_NumberOfDetectors(0) {
}
///////////////////////////////////////////////////////////////////////////
TSamuraiHodoscopePhysics::TSamuraiHodoscopePhysics(){
m_EventData = new TSamuraiHodoscopeData;
m_PreTreatedData = new TSamuraiHodoscopeData;
m_EventPhysics = this;
rawQ_low_threshold=0;
Q_low_threshold=0;
rawQ_high_threshold=4095;
Q_high_threshold=300;
rawT_low_threshold=0;
T_low_threshold=0;
rawT_high_threshold=4095;
T_high_threshold=300;
///////////////////////////////////////////////////////////////////////////
/// A usefull method to bundle all operation to add a detector
void TSamuraiHodoscopePhysics::AddDetector(TVector3 , string ){
// In That simple case nothing is done
// Typically for more complex detector one would calculate the relevant
// positions (stripped silicon) or angles (gamma array)
m_NumberOfDetectors++;
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiHodoscopePhysics::AddDetector(double R, double Theta, double Phi, string shape){
// 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,shape);
}
//m_Spectra = new TSamuraiHodoscopeSpectra;
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiHodoscopePhysics::BuildSimplePhysicalEvent() {
BuildPhysicalEvent();
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiHodoscopePhysics::BuildPhysicalEvent() {
// apply thresholds and calibration
PreTreat();
/*
// match energy and time together
unsigned int mysizeE = m_PreTreatedData->GetMultEnergy();
unsigned int mysizeT = m_PreTreatedData->GetMultTime();
for (UShort_t e = 0; e < mysizeE ; e++) {
for (UShort_t t = 0; t < mysizeT ; t++) {
if (m_PreTreatedData->GetE_DetectorNbr(e) == m_PreTreatedData->GetT_DetectorNbr(t)) {
DetectorNumber.push_back(m_PreTreatedData->GetE_DetectorNbr(e));
Energy.push_back(m_PreTreatedData->Get_Energy(e));
Time.push_back(m_PreTreatedData->Get_Time(t));
unsigned int sizeQUp = m_PreTreatedData->GetMultQUp();
unsigned int sizeTUp = m_PreTreatedData->GetMultTUp();
unsigned int sizeQDown = m_PreTreatedData->GetMultQDown();
unsigned int sizeTDown = m_PreTreatedData->GetMultTDown();
unsigned int IDqu,IDqd,IDtu,IDtd;
double Qu,Qd,Tu,Td;
for(unsigned int qu = 0 ; qu < sizeQUp ; qu++){
Qu = m_PreTreatedData->GetQUp_Charge(qu);
Qd = Tu = Td = -1000;
IDqu = m_PreTreatedData->GetQUp_ID(qu);
// looking for match q down
for(unsigned int qd = 0 ; qd < sizeQDown ; qd++){
IDqd = m_PreTreatedData->GetQDown_ID(qd);
if(IDqd==IDqu){
Qd=m_PreTreatedData->GetQDown_Charge(qd);
break;
}
}
}*/
// valid charge
if(Qd>0){
// Looking for associated time
for(unsigned int tu = 0 ; tu < sizeTUp ; tu++){
IDtu = m_PreTreatedData->GetTUp_ID(tu);
if(IDtu==IDqu){
Tu=m_PreTreatedData->GetTUp_Time(tu);
break;
}
}
// Valide Time Up
if(Tu>0){
for(unsigned int td = 0 ; td < sizeTDown ; td++){
IDtd = m_PreTreatedData->GetTDown_ID(td);
if(IDtd==IDqu){
Td=m_PreTreatedData->GetTDown_Time(td);
break;
}
}
}
// Both Time are valid
if(Td>0){
Charge.push_back(sqrt(Qu*Qd));
Time.push_back(0.5*(Tu+Td));
ID.push_back(IDqu);
}
}
}
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiHodoscopePhysics::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();
/*
// Energy
unsigned int mysize = m_EventData->GetMultEnergy();
for (UShort_t i = 0; i < mysize ; ++i) {
if (m_EventData->Get_Energy(i) > m_E_RAW_Threshold) {
Double_t Energy = Cal->ApplyCalibration("SamuraiHodoscope/ENERGY"+NPL::itoa(m_EventData->GetE_DetectorNbr(i)),m_EventData->Get_Energy(i));
if (Energy > m_E_Threshold) {
m_PreTreatedData->SetEnergy(m_EventData->GetE_DetectorNbr(i), Energy);
}
m_PreTreatedData->Clear();
// UP //
// Q //
unsigned int sizeQUp = m_EventData->GetMultQUp();
for(unsigned int i = 0 ; i < sizeQUp ; i++){
double rawQ = m_EventData->GetQUp_Charge(i)+rand.Uniform(0,1);
if(rawQ>rawQ_low_threshold && rawQ<rawQ_high_threshold){
unsigned int oldID=m_EventData->GetQUp_ID(i);
double Q = (rawQ-m_qup_offset[oldID])*m_qup_gain[oldID];
if(Q>Q_low_threshold)
m_PreTreatedData->SetChargeUp(m_ID[oldID],Q);
}
}
// Time
mysize = m_EventData->GetMultTime();
for (UShort_t i = 0; i < mysize; ++i) {
Double_t Time= Cal->ApplyCalibration("SamuraiHodoscope/TIME"+NPL::itoa(m_EventData->GetT_DetectorNbr(i)),m_EventData->Get_Time(i));
m_PreTreatedData->SetTime(m_EventData->GetT_DetectorNbr(i), Time);
// T //
unsigned int sizeTUp = m_EventData->GetMultTUp();
for(unsigned int i = 0 ; i < sizeTUp ; i++){
unsigned int oldID=m_EventData->GetTUp_ID(i);
double rawT = m_EventData->GetTUp_Time(i)+rand.Uniform(0,1);
if(rawT>rawT_low_threshold && rawT<rawT_high_threshold){
double T = rawT*m_tup_gain[oldID]+m_tup_offset[oldID];
if(T>T_low_threshold && T<T_high_threshold)
m_PreTreatedData->SetTimeUp(m_ID[oldID],T);
}
}
*/
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiHodoscopePhysics::ReadAnalysisConfig() {
bool ReadingStatus = false;
// path to file
string FileName = "./configs/ConfigSamuraiHodoscope.dat";
// open analysis config file
ifstream AnalysisConfigFile;
AnalysisConfigFile.open(FileName.c_str());
if (!AnalysisConfigFile.is_open()) {
cout << " No ConfigSamuraiHodoscope.dat found: Default parameter loaded for Analayis " << FileName << endl;
return;
// Down //
// Q //
unsigned int sizeQDown = m_EventData->GetMultQDown();
for(unsigned int i = 0 ; i < sizeQDown ; i++){
double rawQ = m_EventData->GetQDown_Charge(i)+rand.Uniform(0,1);
if(rawQ>rawQ_low_threshold && rawQ<rawQ_high_threshold){
unsigned int oldID=m_EventData->GetQDown_ID(i);
double Q = (rawQ-m_qdw_offset[oldID])*m_qdw_gain[oldID];
if(Q>Q_low_threshold)
m_PreTreatedData->SetChargeDown(m_ID[oldID],Q);
}
}
cout << " Loading user parameter for Analysis from ConfigSamuraiHodoscope.dat " << endl;
// Save it in a TAsciiFile
TAsciiFile* asciiConfig = RootOutput::getInstance()->GetAsciiFileAnalysisConfig();
asciiConfig->AppendLine("%%% ConfigSamuraiHodoscope.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 = "ConfigSamuraiHodoscope";
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=="E_RAW_THRESHOLD") {
AnalysisConfigFile >> DataBuffer;
m_E_RAW_Threshold = atof(DataBuffer.c_str());
cout << whatToDo << " " << m_E_RAW_Threshold << endl;
}
else if (whatToDo=="E_THRESHOLD") {
AnalysisConfigFile >> DataBuffer;
m_E_Threshold = atof(DataBuffer.c_str());
cout << whatToDo << " " << m_E_Threshold << endl;
}
else {
ReadingStatus = false;
}
// T //
unsigned int sizeTDown = m_EventData->GetMultTDown();
for(unsigned int i = 0 ; i < sizeTDown ; i++){
unsigned int oldID=m_EventData->GetTDown_ID(i);
double rawT=m_EventData->GetTDown_Time(i)+rand.Uniform(0,1);
if(rawT>rawT_low_threshold && rawT<rawT_high_threshold){
double T = rawT*m_tdw_gain[oldID]+m_tdw_offset[oldID];
if(T>T_low_threshold && T<T_high_threshold)
m_PreTreatedData->SetTimeDown(m_ID[oldID],T);
}
}
}
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiHodoscopePhysics::Clear() {
DetectorNumber.clear();
Energy.clear();
ID.clear();
Charge.clear();
Time.clear();
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiHodoscopePhysics::ReadConfiguration(NPL::InputParser parser) {
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("SamuraiHodoscope");
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("SAMURAIHOD");
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << "//// " << blocks.size() << " detectors found " << endl;
......@@ -214,7 +197,7 @@ void TSamuraiHodoscopePhysics::ReadConfiguration(NPL::InputParser parser) {
if(blocks[i]->HasTokenList(mandatory)){
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << endl << "//// SamuraiHodoscope " << i+1 << endl;
string xmlpath= blocks[i]->GetString("XML");
NPL::XmlParser xml;
xml.LoadFile(xmlpath);
......@@ -225,10 +208,31 @@ void TSamuraiHodoscopePhysics::ReadConfiguration(NPL::InputParser parser) {
exit(1);
}
}
// fill the ID to ID map to reorder (probably specific to s034, need to be more modular)
for(int i = 1 ; i<=24; i++)
m_ID[i] = 24-i+1;
for(int i = 27; i <= 40; i++)
m_ID[i] = 40-i+27;
m_ID[25]=26;
m_ID[26]=25;
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiHodoscopePhysics::ReadXML(NPL::XmlParser& xml){
std::vector<NPL::XML::block*> b = xml.GetAllBlocksWithName("SAMURAIHOD");
unsigned int size = b.size();
for(unsigned int i = 0 ; i < size ; i++){
unsigned short ID = b[i]->AsInt("ID");
m_qup_gain [ID] = b[i]->AsDouble("qcal_up");
m_qup_offset [ID] = b[i]->AsDouble("qped_up");
m_qdw_gain [ID] = b[i]->AsDouble("qcal_down");
m_qdw_offset [ID] = b[i]->AsDouble("qped_down");
m_tup_gain [ID] = b[i]->AsDouble("tup_ch2ns");
m_tup_offset [ID] = b[i]->AsDouble("tup_offset");
m_tdw_gain [ID] = b[i]->AsDouble("tdown_ch2ns");
m_tdw_offset [ID] = b[i]->AsDouble("tdown_offset");
}
}
......@@ -264,28 +268,28 @@ void TSamuraiHodoscopePhysics::ClearSpectra() {
///////////////////////////////////////////////////////////////////////////
map< string , TH1*> TSamuraiHodoscopePhysics::GetSpectra() {
// if(m_Spectra)
//return m_Spectra->GetMapHisto();
// else{
map< string , TH1*> empty;
return empty;
// }
// if(m_Spectra)
//return m_Spectra->GetMapHisto();
// else{
map< string , TH1*> empty;
return empty;
// }
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiHodoscopePhysics::WriteSpectra() {
// m_Spectra->WriteSpectra();
// m_Spectra->WriteSpectra();
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiHodoscopePhysics::AddParameterToCalibrationManager() {
CalibrationManager* Cal = CalibrationManager::getInstance();
for (int i = 0; i < m_NumberOfDetectors; ++i) {
Cal->AddParameter("SamuraiHodoscope", "D"+ NPL::itoa(i+1)+"_ENERGY","SamuraiHodoscope_D"+ NPL::itoa(i+1)+"_ENERGY");
Cal->AddParameter("SamuraiHodoscope", "D"+ NPL::itoa(i+1)+"_TIME","SamuraiHodoscope_D"+ NPL::itoa(i+1)+"_TIME");
}
/* CalibrationManager* Cal = CalibrationManager::getInstance();
for (int i = 0; i < m_NumberOfDetectors; ++i) {
Cal->AddParameter("SamuraiHodoscope", "D"+ NPL::itoa(i+1)+"_ENERGY","SamuraiHodoscope_D"+ NPL::itoa(i+1)+"_ENERGY");
Cal->AddParameter("SamuraiHodoscope", "D"+ NPL::itoa(i+1)+"_TIME","SamuraiHodoscope_D"+ NPL::itoa(i+1)+"_TIME");
}*/
}
......@@ -328,14 +332,14 @@ NPL::VDetector* TSamuraiHodoscopePhysics::Construct() {
// Registering the construct method to the factory //
////////////////////////////////////////////////////////////////////////////////
extern "C"{
class proxy_SamuraiHodoscope{
public:
proxy_SamuraiHodoscope(){
NPL::DetectorFactory::getInstance()->AddToken("SAMURAIHOD","SamuraiHodoscope");
NPL::DetectorFactory::getInstance()->AddDetector("SAMURAIHOD",TSamuraiHodoscopePhysics::Construct);
}
};
class proxy_SamuraiHodoscope{
public:
proxy_SamuraiHodoscope(){
NPL::DetectorFactory::getInstance()->AddToken("SAMURAIHOD","Samurai");
NPL::DetectorFactory::getInstance()->AddDetector("SAMURAIHOD",TSamuraiHodoscopePhysics::Construct);
}
};
proxy_SamuraiHodoscope p_SamuraiHodoscope;
proxy_SamuraiHodoscope p_SamuraiHodoscope;
}
......@@ -32,6 +32,7 @@ using namespace std;
#include "TObject.h"
#include "TH1.h"
#include "TVector3.h"
#include "TRandom3.h"
// NPTool headers
#include "TSamuraiHodoscopeData.h"
//#include "TSamuraiHodoscopeSpectra.h"
......@@ -43,48 +44,52 @@ using namespace std;
// forward declaration
//class TSamuraiHodoscopeSpectra;
class TSamuraiHodoscopePhysics : public TObject, public NPL::VDetector {
//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////
// constructor and destructor
public:
TSamuraiHodoscopePhysics();
~TSamuraiHodoscopePhysics() {};
//////////////////////////////////////////////////////////////
// Inherited from TObject and overriden to avoid warnings
public:
void Clear();
void Clear(const Option_t*) {};
//////////////////////////////////////////////////////////////
// data obtained after BuildPhysicalEvent() and stored in
// output ROOT file
public:
vector<int> DetectorNumber;
vector<double> Energy;
vector<int> ID;
vector<double> Charge;
vector<double> Time;
/// A usefull method to bundle all operation to add a detector
void AddDetector(TVector3 POS, string shape);
void AddDetector(double R, double Theta, double Phi, string shape);
private: // XML reading of calibration and position
void ReadXML(NPL::XmlParser&);
// calibration parameter original ID vs coeff
map<unsigned int , double > qup_gain;
map<unsigned int , double > qup_offset;
map<unsigned int , double > tup_gain;
map<unsigned int , double > tup_offset;
map<unsigned int , double > qdw_gain;
map<unsigned int , double > qdw_offset;
map<unsigned int , double > tdw_gain;
map<unsigned int , double > tdw_offset;
map<unsigned int , double > m_qup_gain;//!
map<unsigned int , double > m_qup_offset;//!
map<unsigned int , double > m_tup_gain;//!
map<unsigned int , double > m_tup_offset;//!
map<unsigned int , double > m_qdw_gain;//!