Commit d6217a81 authored by Adrien Matta's avatar Adrien Matta
Browse files

* slight change that seems to suppress segfault in DC reconstruction MT

parent aea3d8ae
Pipeline #116199 failed with stages
in 8 seconds
......@@ -49,7 +49,7 @@ Channel::~Channel(){};
////////////////////////////////////////////////////////////////////////////////
void XmlParser::LoadFile(std::string file){
// First create engine
// First create engine
TXMLEngine* xml = new TXMLEngine;
// Now try to parse xml file
// Only file with restricted xml syntax are supported
......
add_custom_command(OUTPUT TCatanaPhysicsDict.cxx COMMAND ../../scripts/build_dict.sh TCatanaPhysics.h TCatanaPhysicsDict.cxx TCatanaPhysics.rootmap libNPCatana.dylib DEPENDS TCatanaPhysics.h)
add_custom_command(OUTPUT TCatanaDataDict.cxx COMMAND ../../scripts/build_dict.sh TCatanaData.h TCatanaDataDict.cxx TCatanaData.rootmap libNPCatana.dylib DEPENDS TCatanaData.h)
add_library(NPCatana SHARED TCatanaSpectra.cxx TCatanaData.cxx TCatanaPhysics.cxx TCatanaDataDict.cxx TCatanaPhysicsDict.cxx )
target_link_libraries(NPCatana ${ROOT_LIBRARIES} NPCore)
target_link_libraries(NPCatana ${ROOT_LIBRARIES} NPCore NPTrackReconstruction)
install(FILES TCatanaData.h TCatanaPhysics.h TCatanaSpectra.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
#include "MinosUtility.h"
#include "Math/Factory.h"
#include "TROOT.h"
#include <algorithm>
#include <fstream>
#include <vector>
using namespace NPL;
////////////////////////////////////////////////////////////////////////////////
MinosUtility::MinosUtility(){
ROOT::EnableThreadSafety();
// Setting up for signal fitting
m_signal_offset=0;
m_signal_min=ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
m_signal_func=ROOT::Math::Functor(this,&NPL::MinosUtility::signal_chi2,4);
m_signal_func=ROOT::Math::Functor(this,&NPL::MinosUtility::signal_chi2,2);
m_signal_min->SetFunction(m_signal_func);
m_signal_parameter[0]=1000;
m_signal_parameter[1]=250;
m_signal_parameter[2]=15;
m_signal_parameter[3]=0;
m_Dp1=8; // mean distance between t0 and T[0]
m_p2=15;// mean tau
m_Ap0=20;// mean ratio between q max and A
m_signal_min->SetLimitedVariable(0,"a",m_signal_parameter[0],1,0,100000);
m_signal_min->SetLimitedVariable(1,"b",m_signal_parameter[1],1,0,512);
m_signal_min->SetLimitedVariable(2,"c",m_signal_parameter[2],1,1,50);
m_signal_min->SetLimitedVariable(3,"d",m_signal_parameter[3],1,-50,50);
//30 ns sample (33MHz FEM clock) and 333.9 shaping time on preamp, fit with 1/10 of the point
SetParameters(30,333.9,250,10);
m_counter=0;
}
////////////////////////////////////////////////////////////////////////////////
double MinosUtility::signal_shape(const double& x, const double* p){
static double a;
a=(x-p[1])/p[2];
double a = (x-p[1])/m_Tau;
if(x > p[1] && x < 512.)
return (p[0]*exp(-3.*a)*sin(a)*(a*a*a) + p[3]);
return (p[0]*exp(-3.*a)*sin(a)*(a*a*a) + m_Baseline);
else
return p[3];
return m_Baseline;
}
////////////////////////////////////////////////////////////////////////////////
double MinosUtility::signal_chi2(const double* p){
static double chi2, diff,expected;
static unsigned int counter;
counter=0;
chi2=0;
for(unsigned int i = 0 ; i < m_signal_size ; i++ ){
expected = signal_shape((*m_fitSignalT)[i],p);
diff = (*m_fitSignalQ)[i] - expected ;
double chi2=0;
unsigned int counter=0;
double diff,expected;
for(auto it = m_vTQ.begin(); it!=m_vTQ.end() ; it++ ){
expected = signal_shape(it->first,p);
diff = it->second - expected ;
chi2 += diff*diff;
counter++;
//if the track is too long, only the first part of the signal is taken
if((*m_fitSignalT)[i]>t0+50)
break;
}
return chi2/counter;
}
////////////////////////////////////////////////////////////////////////////////
void MinosUtility::sample_signal(){
m_vTQ.clear();
for(unsigned int i = m_minbin; i < m_maxbin ; i+= m_Sampling ){
m_vTQ.push_back(std::make_pair((*m_fitSignalT)[i],(*m_fitSignalQ)[i]));
}
}
////////////////////////////////////////////////////////////////////////////////
void MinosUtility::find_maxQ(){
m_qmax=0; m_qmaxbin=0;
for(unsigned int i = 0 ; i < m_signal_size ; i++ ){
if((*m_fitSignalQ)[i]>m_qmax){
m_qmax=(*m_fitSignalQ)[i];
m_qmaxbin=i;
}
}
}
////////////////////////////////////////////////////////////////////////////////
double MinosUtility::Calibrate(const std::vector<double>& T, const std::vector<double>& Q, double& time, double& charge){
double MinosUtility::Calibrate(const std::vector<unsigned short>* T,const std::vector<unsigned short>* Q, const unsigned int& i , double& time, double& charge){
// setting up for the minisation
m_fitSignalT=&T;
m_fitSignalQ=&Q;
m_signal_size=T.size();
static double delta,qmax,Qmax,Dp1,Ap0;
if(!m_signal_size){
time=-10000;
charge=-10000;
m_fitSignalT=T;
m_fitSignalQ=Q;
m_signal_size=m_fitSignalT->size();
find_maxQ();
if(m_qmaxbin<20){
time=-1;
charge=-1;
return -10000;
}
t0=T[0];
qmax=*(std::max_element(Q.begin(),Q.end()));
Qmax = qmax*m_Ap0;
m_signal_min->SetLimitedVariable(0,"a",Qmax,1,Qmax*0.8,Qmax*1.2);
m_signal_min->SetLimitedVariable(1,"b",t0+m_Dp1,1,t0-10,t0+10);
m_signal_min->SetLimitedVariable(2,"c",m_p2,1,1,50);
sample_signal();
m_guess_t0_bin = m_qmaxbin-20;
m_minbin = m_guess_t0_bin;
m_maxbin = std::min(m_guess_t0_bin+40,m_signal_size);
m_signal_min->Clear();
m_signal_min->SetLimitedVariable(0,"A",m_qmax*10,100,0,m_qmax*20);
m_signal_min->SetLimitedVariable(1,"t0",(*m_fitSignalT)[m_guess_t0_bin],1,0,(*m_fitSignalT)[m_qmaxbin]);
// Perform minimisation
m_signal_min->Minimize();
// access set of parameter that produce the minimum
const double *xs = m_signal_min->X();
const double* xs = m_signal_min->X();
charge= xs[0];
time=xs[1];
// Self optimisation of the starting parameter
m_counter++;
Dp1=time-t0;
m_Dp1 +=(Dp1-m_Dp1)/(m_counter) ;
Ap0=charge/qmax;
m_Ap0 += (Ap0-m_Ap0)/(m_counter);
m_p2 +=(xs[2]-m_p2)/(m_counter);
/*
using namespace std;
ofstream out("signal.txt");
for(unsigned int i = 0 ; i < m_signal_size ; i++){
out << T[i] << " " << Q[i] << endl;
}
out.close();
out.open("fit.txt");
out << xs[0] << " " << xs[1] << " " << xs[2] << " " << xs[3] << endl;
cout << xs[1]-T[0] << " " << xs[1] << " " << xs[2] << " " << xs[3] << " " << delta << " " << endl;
static int first=0;
if(first==40)
exit(0);
first++;
*/
return m_signal_min->MinValue() ;
}
......@@ -22,6 +22,7 @@
* *
*****************************************************************************/
#include<vector>
#include "TMinosData.h"
#include "Math/Minimizer.h"
#include "Math/Functor.h"
......@@ -30,28 +31,39 @@ namespace NPL{
class MinosUtility{
public:
MinosUtility();
~MinosUtility();
~MinosUtility(){};
public:
// take the vector describing the charge, perform a fit and return the charge and time of the trace
double Calibrate(const std::vector<double>& T, const std::vector<double>& Q, double& time, double& charge);
double Calibrate(const std::vector<unsigned short>* T,const std::vector<unsigned short>* Q, const unsigned int& i , double& time, double& charge);
// This function describe the shape of signal
double signal_shape(const double& x, const double* p);
double signal_chi2(const double* p);
void sample_signal(); // take one point every sampling
void find_maxQ();
private:
double m_signal_offset; // offset apply to the signal before fitting
std::vector<std::pair<unsigned short,unsigned short>> m_vTQ;
ROOT::Math::Minimizer* m_signal_min; // minimiser used for the signal fitting
ROOT::Math::Functor m_signal_func; // functor passed to the miniser
const std::vector<double>* m_fitSignalT;
const std::vector<double>* m_fitSignalQ;
unsigned int m_signal_size;
double t0;
double m_signal_parameter[4];
double m_Dp1; // mean distance between t0 and T[0]
double m_p2;// mean tau
double m_Ap0;// mean ratio between q max and A
const std::vector<unsigned short>* m_fitSignalT;
const std::vector<unsigned short>* m_fitSignalQ;
unsigned int m_signal_size,m_minbin,m_maxbin,m_qmax,m_qmaxbin,m_guess_t0_bin;
unsigned int m_counter;
double m_TimeBin; // time between two sample in the waveform (10 to 30 ns depending on the FEM clock)
double m_ShapingTime; // Setting of the preamp on the AGET
unsigned int m_Baseline; // Waveform offset ~250
double m_Tau;// Time constant of the pad signal expressed in bin (depending on TimeBin and ShapingTime)
unsigned int m_Sampling;// Allow to take less point in the wave form to perform the fit to improve speed. 10 gives good results and is twice faster than 1
public:
void SetParameters(double TimeBin,double ShapingTime, unsigned int Baseline, unsigned int Sampling){
m_TimeBin=TimeBin;
m_ShapingTime=ShapingTime;
m_Tau=m_ShapingTime/(log(2)*m_TimeBin);
m_Baseline = Baseline;
m_Sampling = Sampling;};
};
......
/*****************************************************************************
* Copyright (C) 2009-2018 this file is part of the NPTool Project *
* Copyright (C) 2009-2018 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: Elidiano Tronchin contact address: tronchin@lpccaen.in2p3.fr *
* Original Author: Elidiano Tronchin *
* Maintainer : Adrien Matta *
* contact address: matta@lpccaen.in2p3.fr *
* *
* Creation Date : October 2018 *
* Last update : *
* *
* Creation Date : October 2018 *
* Last update : April 2021 *
*---------------------------------------------------------------------------*
* Decription: *
* This class hold Minos Raw data *
* This class hold Minos Raw data *
* *
*---------------------------------------------------------------------------*
* Comment: *
* *
* *
*****************************************************************************/
#include "TMinosData.h"
#include <iostream>
......@@ -37,46 +41,11 @@ TMinosData::TMinosData() {
TMinosData::~TMinosData() {
}
/* TGraph* TMinosData::GetChargeAsGraph(const unsigned int& i){ */
/* TGraph* res = new TGraph (fMinos_Charge[i].size(),&fMinos_Charge[i],&fMinos_Time[i]); */
/* /1* res->Sort(); *1/ */
/* return res; */
/* } */
//////////////////////////////////////////////////////////////////////
void TMinosData::Clear() {
// Minos_Pads
fMinos_PadNumber.clear();
fMinos_PadX.clear();
fMinos_PadY.clear();
/* fMinos_DriftTime.clear(); */
/* fMinos_PadCharge.clear(); */
fMinos_HitCount.clear();
fMinos_Charge.clear();
fMinos_Time.clear();
MINOSx_0.clear();
MINOSy_0.clear();
MINOSz_0.clear();
MINOS_D_min.clear();
MINOS_Radius.clear();
MINOS_NumberTracks.clear();
theta0.clear();
phi0.clear();
energy0.clear();
}
//////////////////////////////////////////////////////////////////////
/* void TMinosData::Dump() const { */
/* // This method is very useful for debuging and worth the dev. */
/* cout << "XXXXXXXXXXXXXXXXXXXXXXXX New Event [TMinosData::Dump()] XXXXXXXXXXXXXXXXX" << endl; */
/* // Energy */
/* size_t mysize = fMinos_E_DetectorNbr.size(); */
/* cout << "Minos_E_Mult: " << mysize << endl; */
/* for (size_t i = 0 ; i < mysize ; i++){ */
/* cout << "DetNbr: " << fMinos_E_DetectorNbr[i] */
/* << " Energy: " << fMinos_Energy[i]; */
/* } */
#ifndef __MinosDATA__
#define __MinosDATA__
#ifndef TMinosDATA_H
#define TMinosDATA_H
/*****************************************************************************
* Copyright (C) 2009-2018 this file is part of the NPTool Project *
* Copyright (C) 2009-2018 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: Elidiano Tronchin contact address: tronchin@lpccaen.in2p3.fr *
* Original Author: Elidiano Tronchin *
* Maintainer : Adrien Matta *
* contact address: matta@lpccaen.in2p3.fr *
* *
* Creation Date : October 2018 *
* Last update : *
* *
* Creation Date : October 2018 *
* Last update : April 2021 *
*---------------------------------------------------------------------------*
* Decription: *
* This class hold Minos Raw data *
* This class hold Minos Raw data *
* *
*---------------------------------------------------------------------------*
* Comment: *
......@@ -27,7 +30,6 @@
// ROOT
#include "TObject.h"
#include "TGraph.h"
using namespace std;
......@@ -38,51 +40,10 @@ class TMinosData : public TObject {
private:
// Pads
vector<UShort_t> fMinos_PadNumber;
vector<Double_t> fMinos_PadX;
vector<Double_t> fMinos_PadY;
/* vector<Double_t> fMinos_PadCharge; */
/* vector<Double_t> fMinos_Charge; */
/* vector<Double_t> fMinos_Time; */
vector< vector<Double_t> > fMinos_Charge;
vector< vector<Double_t> > fMinos_Time;
/* vector<Double_t> fMinos_DriftTime; */
/* vector<UShort_t> fMinos_EventNumber; */
/* //From Santamaria:*/
/* //from simulation*/
/* vector<double> x_tpc,y_tpc,z_tpc,e_tpc;*/
/* vector<double> x_trigger,y_trigger,z_trigger,e_trigger;*/
/* vector<double> x_tar,y_tar,z_tar,e_tar;*/
/* vector<double> x_ch,y_ch,z_ch,e_ch;*/
/* vector<double> x_win,y_win,z_win,e_win;*/
/* vector<double> x_InRoh,y_InRoh,z_InRoh,e_InRoh;*/
/* vector<double> x_OutRoh,y_OutRoh,z_OutRoh,e_OutRoh;*/
/* vector<double> x_Kap,y_Kap,z_Kap,e_Kap;*/
/* double Et_tpc_tot;*/
/* vector<double> Et_tar,Et_ch,Et_tpc,Et_trigger,Et_win,Et_InnerRohacell, Et_OuterRohacell, Et_Kapton;*/
/* vector<int> A, Z;*/
/* vector<int> trackID, parentID;*/
/* */ //unuseful, cause nptool should make already that*/
/* //initial conditions*/
/* //double x0,y0,z0,theta0,phi0,energy0;*/
/* vector<double> x0, y0, z0, theta0, phi0, energy0;*/
/* vector<bool> detection;*/
/* int event;*/
/* */
vector<Double_t> MINOSx_0; //!
vector<Double_t> MINOSy_0; //!
vector<Double_t> MINOSz_0; //!
vector<Double_t> MINOS_D_min; //!
vector<Double_t> MINOS_Radius; //!
vector<Double_t> MINOS_NumberTracks; //!
vector<Double_t> theta0; //!
vector<Double_t> phi0; //!
vector<Double_t> energy0; //!
//For take fitpar values
vector<unsigned short> fMinos_PadNumber;
vector<unsigned short> fMinos_HitCount;
vector< vector<unsigned short> > fMinos_Charge;
vector< vector<unsigned short> > fMinos_Time;
//////////////////////////////////////////////////////////////
// Constructor and destructor
......@@ -96,8 +57,6 @@ class TMinosData : public TObject {
public:
void Clear();
void Clear(const Option_t*) {};
/* void Dump() const; */
//////////////////////////////////////////////////////////////
// Getters and Setters
......@@ -106,93 +65,39 @@ class TMinosData : public TObject {
// add //! to avoid ROOT creating dictionnary for the methods
public:
////////////////////// SETTERS ////////////////////////
// Minos Pads
inline void SetCharge(const UShort_t& Pad,/*const vector<Double_t>& Charge,const vector<Double_t>& Time,*/ const Double_t& X,const Double_t& Y/*,const Double_t& PadCharge*/){
inline void SetPad(const unsigned short& Pad, const unsigned short& HitCount, const vector<int>* time , const vector<int>* charge){
fMinos_HitCount.push_back(HitCount);
fMinos_PadNumber.push_back(Pad);
fMinos_PadX.push_back(X);
fMinos_PadY.push_back(Y);
};//!
inline void AddChargePoint(const vector< Double_t >& Q,const vector<Double_t>& T){
fMinos_Charge.push_back(Q);
fMinos_Time.push_back(T);
};//!
/* inline void AddChargePoint(const double& Q,const double& T){ */
/* fMinos_Charge.push_back(Q); */
/* fMinos_Time.push_back(); */
/* };//! */
//
//Setters for position vertex and obsv in experiment analysis
// Position
inline void SetVertexPos(const Double_t& x,const Double_t& y,const Double_t& z) {
MINOSx_0.push_back(x);
MINOSy_0.push_back(y);
MINOSz_0.push_back(z);
};//!
// Min Distance
inline void SetD_min(const Double_t& dmin) {
MINOS_D_min.push_back(dmin);
};//!
unsigned int size=time->size();
static vector<unsigned short> Charge,Time;
Charge.clear();Time.clear();
for(unsigned int i = 0 ; i < size ; i++){
Time.push_back((*time)[i]);
Charge.push_back((*charge)[i]);
}
fMinos_Time.push_back(Time);
fMinos_Charge.push_back(Charge);
}//!
////////////////////// GETTERS ////////////////////////
inline int GetPadNumberMult()
inline unsigned int GetPadMult() const
{return fMinos_PadNumber.size() ;}//!
/* inline int GetEventNumberMult() */
/* {return fMinos_EventNumber.size() ;}//! */
inline double GetPadX(const unsigned int&i) const
{return fMinos_PadX[i] ;}//!
inline double GetPadY(const unsigned int&i) const
{return fMinos_PadY[i] ;}//!
inline UShort_t GetPadNbr(const unsigned int&i) const
inline unsigned short GetPadNumber(const unsigned int& i) const
{return fMinos_PadNumber[i] ;}//!
/* inline double GetPadCharge(const unsigned int&i) const */
/* {return fMinos_PadCharge[i] ;}//! */
/* inline double GetPadTime(const unsigned int&i) const */
/* {return fMinos_DriftTime[i] ;}//! */
inline vector<double>& GetCharge(const unsigned int&i)
inline vector<unsigned short> GetCharge(const unsigned int& i)const
{return fMinos_Charge[i] ;}//!
inline vector<double>& GetTime(const unsigned int&i)
inline vector<unsigned short> GetTime(const unsigned int& i)const
{return fMinos_Time[i] ;}//!
/* TGraph* GetChargeAsGraph(const unsigned int&i) ; //! */
/* inline vector<double> GetCharge() */
/* {return fMinos_PadChar TGraph* GetEnergyAsGraph();ge() ;}//! */
// Position
inline Double_t GetVertexPos() const
{return MINOSz_0[0] ;}//!
inline Double_t GetVertexPosX() const
{return MINOSx_0[0] ;}//!
inline Double_t GetVertexPosY() const
{return MINOSy_0[0] ;}//!
inline Double_t GetVertexPosZ() const
{return MINOSz_0[0] ;}//!
// Min Distance
inline Double_t GetD_min() const
{return MINOS_D_min[0] ; }//!
// Charge
inline vector<unsigned short>* GetChargePtr(const unsigned int& i)
{return &fMinos_Charge[i];}//!
inline vector<unsigned short>* GetTimePtr(const unsigned int& i)
{return &fMinos_Time[i];}//!
//////////////////////////////////////////////////////////////
// Required for ROOT dictionnary
ClassDef(TMinosData,1) // MinosData structure
ClassDef(TMinosData,2) // MinosData structure
};
#endif
This diff is collapsed.
......@@ -8,7 +8,7 @@
*****************************************************************************/
/*****************************************************************************
* Original Author: Cyril Lenain contact address: lenain@lpccaen.in2p3.fr *
* Original Author: Cyril Lenain contact address: lenain@lpccaen.in2p3.fr *
* *
* Creation Date : 2019 *
* Last update : *
......@@ -30,22 +30,17 @@ using namespace std;
// ROOT headers
#include "TObject.h"
#include "TH1.h"
#include "TVector3.h"
#include "TF1.h"
#include "TH2F.h"
#include "TMinuit.h"
// NPTool headers
#include "TMinosData.h"
#include "MinosUtility.h"
#include "NPLinearRansac3D.h"
#include "NPCalibrationManager.h"
#include "NPVDetector.h"
#include "NPInputParser.h"
#include "Tracking.h"
#include "TMinosResult.h"
#include "TClonesArray.h"
#include "NPXmlParser.h"
using namespace NPL;
// forward declaration
......@@ -67,17 +62,41 @@ class TMinosPhysics : public TObject, public NPL::VDetector {
// data obtained after BuildPhysicalEvent() and stored in
// output ROOT file
public:
// PAD information
vector<unsigned short> Ring_Pad;
vector<double> X_Pad;
vector<double> Y_Pad;
vector<double> Z_Pad;
vector<double> Z_Pad;
vector<double> Q_Pad;
vector<double> T_Pad;
// Tracks information
vector<TVector3> Tracks_P0;
vector<TVector3> Tracks_Dir;