diff --git a/NPLib/Core/NPXmlParser.cxx b/NPLib/Core/NPXmlParser.cxx
index d44367a45f652b1fc5dc36ca2c2113ea1980d4e6..83c756c8267307e2901c75804e76b77dde246d15 100644
--- a/NPLib/Core/NPXmlParser.cxx
+++ b/NPLib/Core/NPXmlParser.cxx
@@ -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
diff --git a/NPLib/Detectors/Catana/CMakeLists.txt b/NPLib/Detectors/Catana/CMakeLists.txt
index a036990ec04e4838f83d40f17553246d0fce780b..9741fddb0827cbc895f4cf86e3d08a76457e1bbe 100644
--- a/NPLib/Detectors/Catana/CMakeLists.txt
+++ b/NPLib/Detectors/Catana/CMakeLists.txt
@@ -1,6 +1,6 @@
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})
diff --git a/NPLib/Detectors/Minos/MinosUtility.cxx b/NPLib/Detectors/Minos/MinosUtility.cxx
index f15d8ac70e65f9deb4587a53bdeee62cabc290b1..52a6e32d33b9231263e83f426ebb73fb2948a8a4 100644
--- a/NPLib/Detectors/Minos/MinosUtility.cxx
+++ b/NPLib/Detectors/Minos/MinosUtility.cxx
@@ -1,121 +1,94 @@
#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() ;
}
diff --git a/NPLib/Detectors/Minos/MinosUtility.h b/NPLib/Detectors/Minos/MinosUtility.h
index af75c2058aea5ba58501564ef5e973aee9dfab35..d6a5682833fe8c5ec781321e94998e1230ddba0a 100644
--- a/NPLib/Detectors/Minos/MinosUtility.h
+++ b/NPLib/Detectors/Minos/MinosUtility.h
@@ -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;};
};
diff --git a/NPLib/Detectors/Minos/TMinosData.cxx b/NPLib/Detectors/Minos/TMinosData.cxx
index e042bdc026e8d3d77ac771747bebbcaeabed10cd..3f6eeb33fb88883302cefff82584739dd9f05e49 100644
--- a/NPLib/Detectors/Minos/TMinosData.cxx
+++ b/NPLib/Detectors/Minos/TMinosData.cxx
@@ -1,24 +1,28 @@
/*****************************************************************************
- * 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]; */
-/* } */
-
diff --git a/NPLib/Detectors/Minos/TMinosData.h b/NPLib/Detectors/Minos/TMinosData.h
index cd6458c80a7b4d1163f9beb2325b89db2f9c7415..01fcb21b49531baaead3d7b02ff5f78dd6e4ddbf 100644
--- a/NPLib/Detectors/Minos/TMinosData.h
+++ b/NPLib/Detectors/Minos/TMinosData.h
@@ -1,20 +1,23 @@
-#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
diff --git a/NPLib/Detectors/Minos/TMinosPhysics.cxx b/NPLib/Detectors/Minos/TMinosPhysics.cxx
index d92b84b55802c209823f54194a40477a1f122d09..d23cedc5d7225810c8a514001da404948095153d 100644
--- a/NPLib/Detectors/Minos/TMinosPhysics.cxx
+++ b/NPLib/Detectors/Minos/TMinosPhysics.cxx
@@ -36,739 +36,298 @@ using namespace std;
#include "RootOutput.h"
#include "NPDetectorFactory.h"
#include "NPOptionManager.h"
-#include "TMinosClust.h"
-#include "TMinosResult.h"
+#include "NPSystemOfUnits.h"
// ROOT
#include "TChain.h"
-#include "TROOT.h"
-#include "TMath.h"
-/* #include "TMinuit2.h" */
-#include "Math/Vector3D.h"
-#include "TStyle.h"
-#include "TClonesArray.h"
-#include "Math/Minimizer.h"
-#include "Math/Factory.h"
TMinosPhysics* current_phy = 0;
///////////////////////////////////////////////////////////////////////////
-TMinosPhysics::TMinosPhysics()
- : m_EventData(new TMinosData),
- m_PreTreatedData(new TMinosData),
- m_EventPhysics(this),
- Tracking_functions(new Tracking),
- m_E_RAW_Threshold(0), // adc channels
- m_E_Threshold(0), // MeV
- hfit(new TH1F("hfit","hfit",512,0,512)),
- grxztmp(new TGraph()),
- gryztmp(new TGraph()),
- npoint_temp(0),
- allevt_2pfiltered(0),
- m_NumberOfDetectors(0) {
- current_phy = this;
- data_result.SetClass("TMinosResult");
- fitdata.SetClass("TMinosClust");
- }
-
-int NclusterFit;
-
-///////////////////////////////////////////////////////////////////////////
-/// A usefull method to bundle all operation to add a detector
-void TMinosPhysics::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 TMinosPhysics::AddDetector(double R, double Theta, TVector3 POS, double Phi, double TargetLenght){
- // Compute the TVector3 corresponding
- // Call the cartesian method
-}
+TMinosPhysics::TMinosPhysics(){
+ m_EventData=new TMinosData;
+ m_EventPhysics=this;
+ m_E_RAW_Threshold=0; // adc channels
+ m_E_Threshold=0; // MeV
+ m_ransac.SetParameters(100,// max 100 iteration
+ 5, // a point belong to a track if it is within 5 mm
+ 40, // max distance allowed between a pair of point to consider a track
+ 10);// minimum point to form a valid cluster
+}
///////////////////////////////////////////////////////////////////////////
void TMinosPhysics::BuildSimplePhysicalEvent() {
BuildPhysicalEvent();
}
-
///////////////////////////////////////////////////////////////////////////
void TMinosPhysics::BuildPhysicalEvent() {
- // apply thresholds and calibration
PreTreat();
-}
+ vector<NPL::LinearCluster3D> clusters = m_ransac.TreatEvent(X_Pad,Y_Pad,Z_Pad);
+ unsigned int sizeC = clusters.size();
+ for(unsigned int i = 0 ; i < sizeC ; i++){
+ // Extract the Direction of the track
+ clusters[i].LinearFit();
+ Tracks_P0.push_back(clusters[i].GetP0());
+ Tracks_Dir.push_back(clusters[i].GetDir());
+ }
-// definition of the fit function for the signal Q(t)
-double TMinosPhysics::conv_fit(double *x, double *p){
- static double p0, p1, p2, x0;
- p0=p[0]; p1=p[1]; p2=p[2];
- x0 = x[0];
- if(x0 > p1 && x0<512.)
- return (p0 * exp(-3.*(x0-p1)/p2) * sin((x0-p1)/p2) * ((x0-p1)/p2)*((x0-p1)/p2)*((x0-p1)/p2) + 250);
- else return (250);
+ if(sizeC==2){
+ static TVector3 Vertex,delta;
+ Delta_Vertex =
+ MinimumDistanceTwoLines(Tracks_P0[0],Tracks_P0[0]+Tracks_Dir[0],
+ Tracks_P0[1],Tracks_P0[1]+Tracks_Dir[1],
+ Vertex, delta);
+ X_Vertex= Vertex.X();
+ Y_Vertex= Vertex.Y();
+ Z_Vertex= Vertex.Z();
+ Theta_12=180.*Tracks_Dir[0].Angle(Tracks_Dir[1])/(M_PI);
+ }
}
///////////////////////////////////////////////////////////////////////////
void TMinosPhysics::PreTreat() {
+ // apply thresholds and calibration
+ static unsigned int sizePad,sizeQ ;
+ sizePad = m_EventData->GetPadMult();
+ static unsigned short PadNumber;
+ static double Q,T;
+ static auto cal= CalibrationManager::getInstance();
+ static string cal_v, cal_o;
+ if(sizePad>20){
+ for(unsigned int i = 0 ; i < sizePad ; i++){
+ vector<unsigned short>* Charge = m_EventData->GetChargePtr(i);
+ if(Charge->size()>40 ){
+ vector<unsigned short>* Time = m_EventData->GetTimePtr(i);
+ m_utility.Calibrate(Time,Charge,i,T,Q);
+
+ if(T>0){
+ PadNumber = m_EventData->GetPadNumber(i);
+ double x_mm = m_X[PadNumber];
+ double y_mm = m_Y[PadNumber];
+ unsigned int ring = round((sqrt(x_mm*x_mm + y_mm*y_mm)-44.15)/2.1);
+ cal_v="Minos/R"+NPL::itoa(ring)+"_VDRIFT";
+ cal_o="Minos/R"+NPL::itoa(ring)+"_OFFSET";
+
+ double z_mm = (T*m_TimeBin+cal->GetValue(cal_o,0))*cal->GetValue(cal_v,0);
+
+ z_mm += m_ZOffset;
+ // Calibrate the Pad:
+ X_Pad.push_back(x_mm);
+ Y_Pad.push_back(y_mm);
+ Z_Pad.push_back(z_mm);
+ Ring_Pad.push_back(ring);
+ Q_Pad.push_back(Q);
+ T_Pad.push_back(T);
+ }
+ }
- // This method applies thresholds and calibrations
- // Isolate 2D tracks with the modified HoughTransformation
- // Calculate z for each pad
-
- data_result.Clear();
- fitdata.Clear();
+ else{
+ /* X_Pad.push_back(-1);
+ Y_Pad.push_back(-1);
+ Z_Pad.push_back((-1);
+ Ring_Pad.push_back(-1);
+ Z_Pad.push_back(-1);
+ Q_Pad.push_back(-1);
+ T_Pad.push_back(-1);*/
+ }
- for(int i = 0; i <2 ; i++){
- parFit1.push_back(-1000);
- parFit2.push_back(-1000);
- parFit3.push_back(-1000);
- parFit4.push_back(-1000);
- }
+ }
- Xvertex = -1000;
- Yvertex = -1000;
- Zvertex = -1000;
- Theta_12= -1000;
-
- Theta_1= -1000;
- Theta_2= -1000;
-
- Phi1 = -1000;
- Phi2 = -1000;
-
- Dmin = -1000;
-
- TMinosClust *minosfitdata;
- ////////////////////fit variables//////////////////////
-
- /* fit_function = new TF1("fit_function",conv_fit, 0, 511, 3); */
- fit_function->Clear();
- int fit2DStatus = 0.;
- double fit_function_max = 0., fit_function_Tpad = 0.;
- double Chi2 = 0.;
- double Tshaping = 333.9;
- const double TimeBinElec = 30.; // in ns
- if(SimulationBool){ // Simulated data
- UperFitLim = 1000000000.; //
- MINOSthresh = 899.;
- VDrift = 0.03475;
- Z_Shift = 0;
}
- else{ // Experiment data
- MINOSthresh = 100;
- UperFitLim = 100000.;
- coef_Vdrift = 1.125;
- ZRot_Minos = 35; // Rotation of Minos along Z axis s034
- //to fixe VDrift
- Z_Shift = -18.95 + 4;
- }
-
- ///////////////////////////////////////////////////////
-
- ClearPreTreatedData();
+}
- /// Read Q(t) signal to fit and extract DriftTime
+///////////////////////////////////////////////////////////////////////////
+void TMinosPhysics::ReadAnalysisConfig() {
+ bool ReadingStatus = false;
+ // path to file
+ string FileName = "./configs/ConfigMinos.dat";
- unsigned int NbOfPad = m_EventData->GetPadNumberMult();
- filled =0;trackNbr=0;
- trackNbr_FINAL=0;
- ringsum=0;
- zmax=0.;Iteration=0;
- for (unsigned int e = 0; e < NbOfPad; e++){
+ // open analysis config file
+ ifstream AnalysisConfigFile;
+ AnalysisConfigFile.open(FileName.c_str());
- x_mm = 0.; y_mm = 0.; z_mm = 0.;maxCharge = 0.; t_pad =0; q_pad=0;
- x_mm = m_EventData->GetPadX(e);
- y_mm = m_EventData->GetPadY(e);
+ if (!AnalysisConfigFile.is_open()) {
+ cout << " No ConfigMinos.dat found: Default parameter loaded for Analayis " << FileName << endl;
+ return;
+ }
+ cout << " Loading user parameter for Analysis from ConfigMinos.dat " << endl;
+
+ // Save it in a TAsciiFile
+ TAsciiFile* asciiConfig = RootOutput::getInstance()->GetAsciiFileAnalysisConfig();
+ asciiConfig->AppendLine("%%% ConfigMinos.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 = "ConfigMinos";
+ 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' );
+ }
- if ( !(abs(x_mm)<0.0001 && abs(y_mm)<0.0001) ) {
- for (UShort_t o = 0; o < m_EventData->GetTime(e).size() ; ++o){
- if(m_EventData->GetCharge(e)[o]> maxCharge){
- maxCharge = m_EventData->GetCharge(e)[o];
- }
+ else if (whatToDo=="E_RAW_THRESHOLD") {
+ AnalysisConfigFile >> DataBuffer;
+ m_E_RAW_Threshold = atof(DataBuffer.c_str());
+ cout << whatToDo << " " << m_E_RAW_Threshold << endl;
}
- if(maxCharge >= MINOSthresh ) { // to remove rings 1,2,3,4 and 16,17,18
- /* if((maxCharge >= MINOSthresh && round((sqrt(x_mm*x_mm+y_mm*y_mm)-44.15)/2.1) < 16 && round((sqrt(x_mm*x_mm+y_mm*y_mm)-44.15)/2.1) >4 )) { // to remove rings 1,2,3,4 and 16,17,18 */
- Xpad.push_back(x_mm);
- Ypad.push_back(y_mm);
- Qpad.push_back(maxCharge);
- filled++; //
+
+ else if (whatToDo=="E_THRESHOLD") {
+ AnalysisConfigFile >> DataBuffer;
+ m_E_Threshold = atof(DataBuffer.c_str());
+ cout << whatToDo << " " << m_E_Threshold << endl;
}
- }
- } //end of NbOfPad
-
- if(filled>0){ // Nbr of activated pads
- while(Xpad.size()>=7 && Iteration< 20 ) { // !!!!!!!!
- /* while(Xpad.size()>=10 && Iteration< 20 ) { // !!!!!!!! */
- filter_result = 0; // Nbr of pads in the Track
- Iteration++;
-
- static vector<double> XpadTemp, YpadTemp, QpadTemp;
- static vector <int> clusterringboolTemp;
- XpadTemp.clear();
- YpadTemp.clear();
- QpadTemp.clear();
- clusterringboolTemp.clear();
-
- filter_result = Tracking_functions->Hough_modified(&Xpad, &Ypad, &Qpad, &XpadTemp, &YpadTemp, &QpadTemp, &clusterringboolTemp);
-
- if(filter_result<0) break;
- if(filter_result>7 ){
- trackNbr++;
- for(int ik=0; ik<filter_result; ik++) { // selected pads
- XpadNew.push_back(XpadTemp[ik]);
- YpadNew.push_back(YpadTemp[ik]);
- ZpadNew.push_back(-10000);
- QpadNew.push_back(QpadTemp[ik]);
- clusterringbool.push_back(clusterringboolTemp[ik]);
- clusternbr.push_back(Iteration);
- clusterpads.push_back(filter_result);
- }
+
+ else {
+ ReadingStatus = false;
}
}
}
- for(unsigned int il=0; il<XpadNew.size(); il++){
- minosfitdata = (TMinosClust*)fitdata.ConstructedAt(il);
- minosfitdata->Set(XpadNew[il], YpadNew[il],
- -10000, -10000, QpadNew[il],
- clusternbr[il],clusterpads[il], 0.);
- /* ZpadNew.push_back(-10000.); */
- }
-
- if(trackNbr == 2 || trackNbr == 1){ ///////////
-
- /* if(filled==0) cerr << "Error !!!" << endl; */
- //-------------------------------------------------------
- // STEP 4.1: Fitting taken pads for Qmax and Ttrig info
- //-------------------------------------------------------
- for(unsigned int i=0 ; i< NbOfPad; i++ ) {
- hfit->Reset();
- bool fitbool = false;
- x_mm = m_EventData->GetPadX(i);
- y_mm = m_EventData->GetPadY(i);
-
- for(unsigned int jj=0; jj<XpadNew.size(); jj++) {
- if( abs(XpadNew[jj]-x_mm)<0.0001 && abs(YpadNew[jj]-y_mm)<0.0001) {
- fitbool = true;
- indexfill=jj;
- break;
- }
- }
- if( fitbool==true ) {
-
- /* TGraph gfit( m_EventData->GetTime(i).size() ,&(m_EventData->GetTime(i)[0]),&(m_EventData->GetCharge(i)[0])); */
- /* for (UShort_t o = 0; o < m_EventData->GetTime(i).size() ; ++o){ */
- /* gfit.SetPoint(o,m_EventData->GetTime(i)[o],m_EventData->GetCharge(i)[o]+250); */
- /* if(m_EventData->GetCharge(i)[o]+250> hfit_max){ */
- /* hfit_max_T = m_EventData->GetTime(i)[o]; */
- /* hfit_max=m_EventData->GetCharge(i)[o]+250; */
- /* } */
- /* } */
-
- for(Int_t j=0; j< m_EventData->GetTime(i).size(); j++) {
- if(m_EventData->GetCharge(i)[j]>=0){
- hfit->SetBinContent(hfit->FindBin(m_EventData->GetTime(i)[j]), m_EventData->GetCharge(i)[j]+250);
- }
- }
-
- // Fitting the hfit histogram of last channel if not empty
- if(hfit->GetSumOfWeights()>3000) {
- hfit->GetXaxis()->SetRange(0,510);
- hfit_max = hfit->GetMaximum();
- hfit_max_T = hfit->GetMaximumBin();
- T_min=-1;
- T_max=-1;
-
- // Find the T_min & T_max limits of the signal non zero
- for(int h=hfit_max_T;h>0;h--) {
- if(T_min == -1 && (hfit->GetBinContent(h))<=250 ) {
- T_min = h;
- break;
- }
- }
- for(int h=hfit_max_T;h<510;h++) {
- if(T_max == -1 && (hfit->GetBinContent(h)) < 1 ) {
- T_max = h;
- break;
- }
- }
-
- if((hfit_max_T-3.5*(Tshaping/TimeBinElec)) > T_min) T_min = hfit_max_T-2*Tshaping/TimeBinElec;
- if((hfit_max_T+10) < T_max || T_max==-1) T_max = hfit_max_T+10.;
- T_min = max(T_min,0.);
- if(T_max>510) T_max = 510;
-
- // Set fit parameters
- fit_function->SetParameter(0, hfit_max-250.);
- fit_function->SetParameter(1,hfit_max_T - Tshaping/TimeBinElec);
- fit_function->SetParameter(2, Tshaping/TimeBinElec);
- fit_function->SetParLimits(0,0,UperFitLim);
- fit_function->SetParLimits(1,-20,512);
- fit_function->SetParLimits(2,0,512);
-
- fit2DStatus = hfit->Fit(fit_function,"QN","",T_min,T_max);
- /* fit2DStatus = gfit.Fit(fit_function,"QN","",hfit_max-2*Tshaping,hfit_max+2*Tshaping); */
-
- double fit_function_max = 0., fit_function_Tpad = 0.;
- if(fit2DStatus==0) {
- Chi2 = fit_function->GetChisquare();
- fit_function_max = fit_function->GetMaximum();
- fit_function_Tpad = fit_function->GetParameter(1);
- }
- //attribute q_pad and z_mm value
- if(fit2DStatus!=0 || fit_function_max<=20. ||
- fit_function_max >= UperFitLim || fit_function_Tpad<=0.15 ||
- fit_function_Tpad>=513. || fit_function->GetParameter(2)<=0.15 ||
- fit_function->GetParameter(2)>=513.) {
- q_pad = hfit_max-250.;
- z_mm = -10000;
- }
- else {
- t_pad = fit_function_Tpad;
- if(SimulationBool)z_mm = t_pad*TimeBinElec*VDrift; // for simu
- else{
- int ring = round((sqrt(x_mm*x_mm + y_mm*y_mm)-44.15)/2.1);
- z_mm =(t_pad*TimeBinElec+DelayTrig[ring])*(VdriftperRing[ring]*coef_Vdrift);
- }
- Q_Pad.push_back(fit_function_max-250);
- T_Pad.push_back(t_pad*TimeBinElec);
- X_Pad.push_back(x_mm);
- Y_Pad.push_back(y_mm);
- Z_Pad.push_back(z_mm);
- q_pad = fit_function_max-250.;
-
- }
- ZpadNew[indexfill] = z_mm + Z_Shift;
- QpadNew[indexfill] = q_pad;
- minosfitdata = (TMinosClust*)fitdata.ConstructedAt(indexfill);
- minosfitdata->Set(XpadNew[indexfill], YpadNew[indexfill], t_pad*TimeBinElec, z_mm, q_pad, clusternbr[indexfill], clusterpads[indexfill], Chi2);
- } // end if SumOfWeights > 0
- } // end if fitbool==True
- else continue;
- } // end of NbOfPad
-
- //-------------------------------------------------------
- // STEP 3.2: Filtering the tracks off possible noise
- // with Hough3D (3*2D planes)
- //-------------------------------------------------------
-
- cluster_temp = 0;
- int ringtouch[19]={0};
- static vector<double> xin, yin, zin, qin;
- static vector<double> xout, yout, zout, qout;
- static vector<double> xoutprime, youtprime, zoutprime, qoutprime;
-
- for(unsigned int i=0;i<(XpadNew.size());i++){
- if(xin.size()>0 && ((cluster_temp != int(clusternbr[i]) && i!=0) || i==(XpadNew.size()-1))){ // We fill xin until the next cluster
- Tracking_functions->Hough_3D(&xin, &yin, &zin, &qin, &xout, &yout, &zout, &qout);
- for(unsigned int ij=0; ij<xout.size();ij++){
- if(zout[ij]>zmax) {zmax = zout[ij];}
- ringtouch[int(round((sqrt(xout[ij]*xout[ij]+yout[ij]*yout[ij])-44.15)/2.1))]++; // Corr by Cyril
- }
- for(int ko=0; ko<19; ko++){
- if(ringtouch[ko]>0) ringsum++;
- }
- // Tracks of interest: >10 pads and >=12 rings hit
- if(xout.size()>10 && ringsum>=8){
- npoint=0;
- trackNbr_FINAL++;
- double charge_temp=0.;
- double lenght_temp=0;
- double zintarget=300;
- double zouttarget=0;
- int indexin=0; int indexout=0;
- for(unsigned int ij=0; ij<xout.size(); ij++){
- point.SetXYZ(xout[ij],yout[ij],zout[ij]);
- if(!SimulationBool)point.RotateZ(ZRot_Minos*TMath::DegToRad());//15.6*TMath::DegToRad() CHANGE####
- xoutprime.push_back(point.X());youtprime.push_back(point.Y());zoutprime.push_back(point.Z());
- TOTxoutprime.push_back(point.X());
- // Added line to see events used for fit of lines in tree
- TOTyoutprime.push_back(point.Y());
- TOTzoutprime.push_back(point.Z());
- TOTqout.push_back(qout[ij]);
- // save cluster_temp
- // save numbtrack
- trackclusternbr.push_back(clusternbr[i]);
- tracknbr.push_back( trackNbr_FINAL ) ;
-
- grxztmp->SetPoint(npoint,zoutprime[ij],xoutprime[ij]);
- gryztmp->SetPoint(npoint,zoutprime[ij],youtprime[ij]);
- charge_temp += qout[ij];
- minosdata_result = (TMinosResult*)data_result.ConstructedAt(array_final);
- minosdata_result->Set(xoutprime[ij], youtprime[ij], zoutprime[ij], qout[ij], trackNbr_FINAL, xout.size(), zmax);
- array_final++;
- npoint++;
-
- if(zoutprime[ij]<zintarget) {zintarget=zoutprime[ij];indexin=ij;}
- if(zoutprime[ij]>zouttarget) {zouttarget=zoutprime[ij];indexout=ij;}
- }
-
- if(xout.size()>0)lenght_temp=sqrt((zoutprime[indexin]-zoutprime[indexout])*(zoutprime[indexin]-zoutprime[indexout]) +(youtprime[indexin]-youtprime[indexout])*(youtprime[indexin]-youtprime[indexout]) +(xoutprime[indexin]-xoutprime[indexout])*(xoutprime[indexin]-xoutprime[indexout]));
-
- grxz.push_back(*grxztmp);
- gryz.push_back(*gryztmp);
- grxztmp->Set(0);
- /* if(maxCharge >= MINOSthresh ) { */
- gryztmp->Set(0);
-
- chargeTot.push_back(charge_temp);
- lenght.push_back(lenght_temp);
- } // end of if(xout.size()>10 &&)
-
- xin.clear(); yin.clear(); zin.clear(); qin.clear();
- xout.clear(); yout.clear(); zout.clear(); qout.clear();
- xoutprime.clear(); youtprime.clear(); zoutprime.clear();
-
- npoint_temp=0; ringsum=0; zmax=0.;
- for(int ko=0; ko<18; ko++)
- ringtouch[ko] = 0;
-
- } // if(xin.size()>0 && (( cluster_temp .... )
-
- cluster_temp = clusternbr[i]; // Number of the track/cluster
- if(!(clusterpads[i]>=10 && clusterringbool[i]==1 && ZpadNew[i]>-100 && ZpadNew[i]<=310)) continue;
- else
- {
- xin.push_back(XpadNew[i]);
- yin.push_back(YpadNew[i]);
- zin.push_back(ZpadNew[i]);
- qin.push_back(QpadNew[i]);
- npoint_temp++;
- }
-
- }//end of PadNews
-
- /* //------------------------------------------------------- */
- /* // STEP 3.2: Fitting the filtered tracks in 3D */
- /* // (weight by charge, TMinuit) */
- /* //------------------------------------------------------- */
-
-
- if(trackNbr_FINAL== 2 || trackNbr_FINAL == 1){
-
- //////////Minimization in 2D to reconstruct track lines
- allevt_2pfiltered++;
- for(int itr= 0 ; itr < trackNbr_FINAL; itr++) {
- pStart[0]=0; pStart[2]=0; pStart[1]=1; pStart[3]=3;
-
- min = new TMinuit(4);
- min->SetPrintLevel(-1);
- arglist[0] = 3;
-
- Tracking_functions->FindStart(pStart,chi,fitStatus, &grxz.at(itr), &gryz.at(itr));
-
- NclusterFit = itr+1;
- current_phy=this;
-
- min->SetFCN(SumDistance);
-
- // Set starting values and step sizes for parameters
- min->mnparm(0,"x0",pStart[0],0.1,-500,500,iflag);
- min->mnparm(1,"Ax",pStart[1],0.1,-10,10,iflag);
- min->mnparm(2,"y0",pStart[2],0.1,-500,500,iflag);
- min->mnparm(3,"Ay",pStart[3],0.1,-10,10,iflag);
-
- arglist[0] = 200; // number of function calls
- arglist[1] = 0.000001; // tolerance
-
- min->mnexcm("MIGRAD",arglist,2,iflag); // minimization with MIGRAD
- min->mnstat(amin,edm,errdef,nvpar,nparx,iflag); //returns current status of the minimization
-
- // get fit parameters
- for(int i = 0; i <4; i++) min->GetParameter(i,parFit_temp[i],err_temp[i]);
-
- /* if( (parFit_temp[0] >-499 && parFit_temp[0]<499) && (parFit_temp[2] >-499 && parFit_temp[2]<499)) { */
- /* parFit1[itr] = push_back(parFit_temp[0]); */
- /* parFit2.push_back(parFit_temp[1]); */
- /* parFit3.push_back(parFit_temp[2]); */
- /* parFit4.push_back(parFit_temp[3]); */
- /* } */
- parFit1[itr]=parFit_temp[0];
- parFit2[itr]=parFit_temp[1];
- parFit3[itr]=parFit_temp[2];
- parFit4[itr]=parFit_temp[3];
-
- delete min;
- }
-
- static double ParTrack1[4], ParTrack2[4];
- static double VectorTrack11[3], VectorTrack22[3];
-
- ParTrack1[0] = parFit1[0];
- ParTrack1[1] = parFit2[0];
- ParTrack1[2] = parFit3[0];
- ParTrack1[3] = parFit4[0];
-
- if(trackNbr_FINAL==2){
- ParTrack2[0] = parFit1[1];
- ParTrack2[1] = parFit2[1];
- ParTrack2[2] = parFit3[1];
- ParTrack2[3] = parFit4[1];
- }
- else if(trackNbr_FINAL==1){// The vertex is still calculated with Zaxis
- ParTrack2[0] = 0;
- ParTrack2[1] = 0;
- ParTrack2[2] = 0;
- ParTrack2[3] = 0;
- }
- Dmin=-100, Theta_1 = -1, Theta_2 = -1;
-
- Tracking_functions->vertex(ParTrack1, ParTrack2, Xvertex, Yvertex, Zvertex, Dmin, Theta_1, Theta_2, Phi1, Phi2, VectorTrack11, VectorTrack22);
+}
- VectorTrack1.SetXYZ(VectorTrack11[0],VectorTrack11[1],VectorTrack11[2]);
- VectorTrack2.SetXYZ(VectorTrack22[0],VectorTrack22[1],VectorTrack22[2]);
- VectorTrack1 = VectorTrack1.Unit();
- VectorTrack2 = VectorTrack2.Unit();
- Theta_12 = VectorTrack1.Angle(VectorTrack2)*180/TMath::Pi();
- }// end if trackNbr_FINAL>=1
+///////////////////////////////////////////////////////////////////////////
+void TMinosPhysics::Clear() {
+ Ring_Pad.clear();
+ X_Pad.clear();
+ Y_Pad.clear();
+ Z_Pad.clear();
+ Q_Pad.clear();
+ T_Pad.clear();
+
+ // Tracks information
+ Tracks_P0.clear();
+ Tracks_Dir.clear();
+
+ // Vertex information
+ X_Vertex=-1000;
+ Y_Vertex=-1000;
+ Z_Vertex=-1000;
+ Theta_12=-1000;
+ Delta_Vertex=-1000;
+}
- } // end loop if 0<trackNbr < 5
- // instantiate CalibrationManager
- static CalibrationManager* Cal = CalibrationManager::getInstance();
- }
+///////////////////////////////////////////////////////////////////////////
+void TMinosPhysics::ReadConfiguration(NPL::InputParser parser) {
+ vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Minos");
+ if(NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << "//// " << blocks.size() << " detector(s) found " << endl;
- ///////////////////////////////////////////////////////////////////////////
- void TMinosPhysics::ReadAnalysisConfig() {
- bool ReadingStatus = false;
- // path to file
- string FileName = "./configs/ConfigMinos.dat";
+ vector<string> token= {"XML","TimeBin","ShapingTime","Baseline","Sampling","ZOffset"};
- // open analysis config file
- ifstream AnalysisConfigFile;
- AnalysisConfigFile.open(FileName.c_str());
+ for(unsigned int i = 0 ; i < blocks.size() ; i++){
+
+ if (blocks[i]->HasTokenList(token)) {
+ cout << endl << "//// MINOS (" << i+1 << ")" << endl;
+ m_ShapingTime = blocks[i]->GetDouble("ShapingTime","ns")/NPUNITS::ns;
+ m_TimeBin = blocks[i]->GetDouble("TimeBin","ns")/NPUNITS::ns;
+ m_Sampling= blocks[i]->GetInt("Sampling");
+ m_Baseline= blocks[i]->GetInt("BaseLine");
+ m_utility.SetParameters(m_TimeBin,m_ShapingTime,m_Baseline,m_Sampling);
+ m_ZOffset = blocks[i]->GetDouble("ZOffset","mm");
+ string xmlpath = blocks[i]->GetString("XML");
+ NPL::XmlParser xml;
+ xml.LoadFile(xmlpath);
+ ReadXML(xml);
- if (!AnalysisConfigFile.is_open()) {
- cout << " No ConfigMinos.dat found: Default parameter loaded for Analayis " << FileName << endl;
- return;
}
- cout << " Loading user parameter for Analysis from ConfigMinos.dat " << endl;
-
- // Save it in a TAsciiFile
- TAsciiFile* asciiConfig = RootOutput::getInstance()->GetAsciiFileAnalysisConfig();
- asciiConfig->AppendLine("%%% ConfigMinos.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 = "ConfigMinos";
- 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;
- }
- }
+ else {
+ cout << "ERROR: Missing token for Minos, check your input file"<< endl;
+ cout << "Required token: ";
+ for(unsigned int i = 0 ; i < token.size() ; i++)
+ cout << token[i] << " " ;
+ cout << endl;
+ exit(1);
}
}
-
- double TMinosPhysics::distance2(double x,double y,double z, double *p) {
- // distance line point is D= | (xp-x0) cross ux |
- // where ux is direction of line and x0 is a point in the line (like t = 0)
- ROOT::Math::XYZVector xp(x,y,z); //point of the track
- ROOT::Math:: XYZVector x0(p[0], p[2], 0. );
- ROOT::Math::XYZVector x1(p[0] + p[1], p[2] + p[3], 1. ); //line
- ROOT::Math::XYZVector u = (x1-x0).Unit();
- double d2 = ((xp-x0).Cross(u)) .Mag2();
- return d2;
- }
-
- void TMinosPhysics::SumDistance(int &, double *, double & sum, double * par, int) {
- TMinosPhysics* phy = current_phy;
- int nused=0;
- double qtot=0;
- sum = 0;
-
- for(int i=0; i < phy->data_result.GetEntriesFast(); i++)
- {
- phy->minosdata_result = (TMinosResult*)phy->data_result.At(i);
- if(phy->minosdata_result->n_Cluster==NclusterFit)
- {
- float x=phy->minosdata_result->x_mm;
- float y=phy->minosdata_result->y_mm;
- float z=phy->minosdata_result->z_mm;
- float q=phy->minosdata_result->Chargemax;
- //if(nused<2)cout<<minosdata_result->n_Cluster<<" "<<x<<" "<<y<<" "<<z<<" "<<q<<endl;
- double d = TMinosPhysics::distance2(x, y, z, par);
- sum += d*q;
- nused++;
- qtot+=q;
- }
- }
- //sum/=nused;
- sum/=qtot;
- return;
+}
+///////////////////////////////////////////////////////////////////////////
+void TMinosPhysics::ReadXML(NPL::XmlParser& xml){
+ std::vector<NPL::XML::block*> b = xml.GetAllBlocksWithName("MINOS");
+ unsigned int size = b.size();
+ unsigned int valid_pad=0;
+ for(unsigned int i = 0 ; i < size ; i++){
+
+ unsigned short ID = b[i]->AsInt("ID");
+ m_X[ID] = b[i]->AsDouble("X");
+ m_Y[ID] = b[i]->AsDouble("Y");
+ valid_pad++;
+ if((m_X[ID]==0&&m_Y[ID]==0)||(m_X[ID]==-1&&m_Y[ID]==-1))
+ valid_pad--;
+ m_Period[ID] = b[i]->AsDouble("TPERIOD");
+ m_Offset[ID] = b[i]->AsDouble("TOFFSET");
+ m_Gain[ID] = b[i]->AsDouble("GAIN");
}
- ///////////////////////////////////////////////////////////////////////////
- void TMinosPhysics::Clear() {
-
- Xpad.clear();
- Ypad.clear();
- Qpad.clear();
- XpadNew.clear();
- YpadNew.clear();
- ZpadNew.clear();
- QpadNew.clear();
-
- clusterringboolTemp.clear();
- clusterringbool.clear();
- clusternbr.clear();
- clusterpads.clear();
- hfit->Reset();
-
- /* xoutprime.clear(); */
- /* youtprime.clear(); */
- /* zoutprime.clear(); */
-
- trackclusternbr.clear();
- tracknbr.clear();
- TOTxoutprime.clear();
- TOTyoutprime.clear();
- TOTzoutprime.clear();
- TOTqout.clear();
-
- lenght.clear();
- chargeTot.clear();
- parFit1.clear();
- parFit2.clear();
- parFit3.clear();
- parFit4.clear();
-
- grxz.clear();
- gryz.clear();
-
- hfit->Reset();
-
- filled=0;
- indexfill = 0;
- ChargeBin = 0.;
- maxCharge = 0.;
- Iteration=0;
- filter_result=0;
- fit2DStatus=0;
- trackNbr=0;
- trackNbr_FINAL=0;
- x_mm = 0.; y_mm = 0.; z_mm = 0.; q_pad = 0.; t_pad = 0.;
- array_final=0;
- ringsum=0;
- zmax=0.;
+ if(NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << "//// " << m_X.size() << " pads found with " << valid_pad << " valid pads" << endl;
- }
+}
- ///////////////////////////////////////////////////////////////////////////
- void TMinosPhysics::ReadConfiguration(NPL::InputParser parser) {
- }
- ///////////////////////////////////////////////////////////////////////////
- void TMinosPhysics::AddParameterToCalibrationManager() {
- CalibrationManager* Cal = CalibrationManager::getInstance();
- for (int i = 0; i < m_NumberOfDetectors; ++i) {
- Cal->AddParameter("Minos", "D"+ NPL::itoa(i+1)+"_ENERGY","Minos_D"+ NPL::itoa(i+1)+"_ENERGY");
- Cal->AddParameter("Minos", "D"+ NPL::itoa(i+1)+"_TIME","Minos_D"+ NPL::itoa(i+1)+"_TIME");
- }
+///////////////////////////////////////////////////////////////////////////
+void TMinosPhysics::AddParameterToCalibrationManager() {
+ CalibrationManager* Cal = CalibrationManager::getInstance();
+ unsigned int NbrRing = 18;
+ vector<double> standardV={0.03475};
+ vector<double> standardO={0};
+ for(unsigned int i = 0 ; i < NbrRing ; i++){
+ Cal->AddParameter("Minos", "R"+NPL::itoa(i+1)+"_VDRIFT","Minos_R"+NPL::itoa(i+1)+"_VDRIFT",standardV);
+ Cal->AddParameter("Minos", "R"+NPL::itoa(i+1)+"_OFFSET","Minos_R"+NPL::itoa(i+1)+"_OFFSET",standardO);
}
+}
- ///////////////////////////////////////////////////////////////////////////
- void TMinosPhysics::InitializeRootInputRaw() {
+///////////////////////////////////////////////////////////////////////////
+void TMinosPhysics::InitializeRootInputRaw() {
+ TChain* inputChain = RootInput::getInstance()->GetChain();
+ inputChain->SetBranchStatus("Minos", true );
+ inputChain->SetBranchAddress("Minos", &m_EventData );
+}
- TChain* inputChain = RootInput::getInstance()->GetChain();
- inputChain->SetBranchStatus("Minos", true );
- inputChain->SetBranchAddress("Minos", &m_EventData );
+///////////////////////////////////////////////////////////////////////////
+void TMinosPhysics::InitializeRootInputPhysics() {
+ TChain* inputChain = RootInput::getInstance()->GetChain();
+ inputChain->SetBranchAddress("Minos", &m_EventPhysics);
+}
- fit_function = new TF1("fit_function",conv_fit, 0, 511, 3);
+///////////////////////////////////////////////////////////////////////////
+void TMinosPhysics::InitializeRootOutput() {
+ TTree* outputTree = RootOutput::getInstance()->GetTree();
+ outputTree->Branch("Minos", "TMinosPhysics", &m_EventPhysics);
+}
- if(NPOptionManager::getInstance()->HasDefinition("simulation")){
- cout << "Considering input data as simulation"<< endl;
- SimulationBool = true;
- }
- else{
- cout << "Considering input data as real" << endl;
-
- SimulationBool = false;
-
- ifstream calibFile2("Vdrift.txt");
- string buffer2;
- getline(calibFile2, buffer2);
- double vdriftR;
- int i = 0;
- while(calibFile2 >> vdriftR){
- VdriftperRing[i] = vdriftR; // ns, s034 par.
- i++;
- }
+////////////////////////////////////////////////////////////////////////////////
+// Construct Method to be pass to the DetectorFactory //
+////////////////////////////////////////////////////////////////////////////////
+NPL::VDetector* TMinosPhysics::Construct() {
+ return (NPL::VDetector*) new TMinosPhysics();
+}
- ifstream calibFile("Time_Offset.txt");
- string buffer;
- getline(calibFile, buffer);
- double offset;
- i = 0;
- while(calibFile >> offset){
- DelayTrig[i] = offset; // ns, s034 par.
- i++;
+////////////////////////////////////////////////////////////////////////////////
+// Registering the construct method to the factory //
+////////////////////////////////////////////////////////////////////////////////
+extern "C"{
+ class proxy_Minos{
+ public:
+ proxy_Minos(){
+ NPL::DetectorFactory::getInstance()->AddToken("Minos","Minos");
+ NPL::DetectorFactory::getInstance()->AddDetector("Minos",TMinosPhysics::Construct);
}
- }
- }
-
- ///////////////////////////////////////////////////////////////////////////
- void TMinosPhysics::InitializeRootInputPhysics() {
- TChain* inputChain = RootInput::getInstance()->GetChain();
- inputChain->SetBranchAddress("Minos", &m_EventPhysics);
- }
-
- ///////////////////////////////////////////////////////////////////////////
- void TMinosPhysics::InitializeRootOutput() {
- TTree* outputTree = RootOutput::getInstance()->GetTree();
- outputTree->Branch("Minos", "TMinosPhysics", &m_EventPhysics);
- }
-
- ////////////////////////////////////////////////////////////////////////////////
- // Construct Method to be pass to the DetectorFactory //
- ////////////////////////////////////////////////////////////////////////////////
- NPL::VDetector* TMinosPhysics::Construct() {
- return (NPL::VDetector*) new TMinosPhysics();
- }
-
- ////////////////////////////////////////////////////////////////////////////////
- // Registering the construct method to the factory //
- ////////////////////////////////////////////////////////////////////////////////
- extern "C"{
- class proxy_Minos{
- public:
- proxy_Minos(){
- NPL::DetectorFactory::getInstance()->AddToken("Minos","Minos");
- NPL::DetectorFactory::getInstance()->AddDetector("Minos",TMinosPhysics::Construct);
- }
- };
+ };
proxy_Minos p_Minos;
- }
-
+}
diff --git a/NPLib/Detectors/Minos/TMinosPhysics.h b/NPLib/Detectors/Minos/TMinosPhysics.h
index 3ab14ef5a44d7ef68c65d37463c5316ca5a455be..4e7f9aac9b49ee178683a1a822855e39595b6b63 100644
--- a/NPLib/Detectors/Minos/TMinosPhysics.h
+++ b/NPLib/Detectors/Minos/TMinosPhysics.h
@@ -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;
+
+ // Vertex information
+ double X_Vertex,Y_Vertex,Z_Vertex,Theta_12,Delta_Vertex;
+
+ private: // used for analysis
+ double m_TimeBin;//!
+ double m_ShapingTime;//!
+ double m_Baseline;//!
+ unsigned int m_Sampling;//!
+ double m_ZOffset;//!
+
+
+ NPL::MinosUtility m_utility;//! // an utility to fit the pad signal
+ NPL::LinearRansac3D m_ransac;//! // a linear ransac to build the 3D tracks
+
+ double m_E_RAW_Threshold; //! adc channels
+ double m_E_Threshold; //! MeV
+ std::map<unsigned short,double> m_X;//!
+ std::map<unsigned short,double> m_Y;//!
+ std::map<unsigned short,double> m_Period;//!
+ std::map<unsigned short,double> m_Offset;//!
+ std::map<unsigned short,double> m_Gain;//!
+
- /// A usefull method to bundle all operation to add a detector
- void AddDetector(TVector3 POS, string shape);
- void AddDetector(double R, double Theta,TVector3 POS, double Phi, double TargetLenght);
-
//////////////////////////////////////////////////////////////
// methods inherited from the VDetector ABC class
public:
@@ -112,26 +131,9 @@ class TMinosPhysics : public TObject, public NPL::VDetector {
// specific methods to Minos array
public:
- NPL::Tracking* Tracking_functions; //!
-
- // fit function for the Q(t) signals at each pad
- static double conv_fit(double *x, double *p);
- static double distance2(double x,double y,double z, double *p);
- static void SumDistance(int &, double *, double & sum, double * par, int);
-
// remove bad channels, calibrate the data and apply thresholds
void PreTreat();
- // clear the pre-treated object
- void ClearPreTreatedData() {
-
- m_PreTreatedData->Clear();
- Q_Pad.clear();
- X_Pad.clear();
- Y_Pad.clear();
- Z_Pad.clear();
- T_Pad.clear();
- }
// read the user configuration file. If no file is found, load standard one
void ReadAnalysisConfig();
@@ -140,173 +142,29 @@ class TMinosPhysics : public TObject, public NPL::VDetector {
// needed for online analysis for example
void SetRawDataPointer(TMinosData* rawDataPointer) {m_EventData = rawDataPointer;}
+ // Read and load the XML info
+ void ReadXML(NPL::XmlParser&);
// objects are not written in the TTree
private:
TMinosData* m_EventData; //!
- TMinosData* m_PreTreatedData; //!
TMinosPhysics* m_EventPhysics; //!
- TMinuit* min;//!
-TClonesArray data_result;//!
-TMinosResult *minosdata_result;//!
-
- TClonesArray fitdata;//!
-
- // getters for raw and pre-treated data object
+ // getters for raw and pre-treated data object
public:
- TMinosData* GetRawData() const {return m_EventData;}
- TMinosData* GetPreTreatedData() const {return m_PreTreatedData;}
-
- vector<double> GetPad_X() {return X_Pad;}
- vector<double> GetPad_Y() {return Y_Pad;}
- vector<double> GetPad_Z() {return Z_Pad;}
- vector<double> GetPad_T() {return T_Pad;}
- vector<double> GetParFit1() {return parFit1;}
- vector<double> GetParFit2() {return parFit2;}
- vector<double> GetParFit3() {return parFit3;}
- vector<double> GetParFit4() {return parFit4;}
- double GetVertexZ() {return Zvertex;}
- double GetVertexX() {return Xvertex;}
- double GetVertexY() {return Yvertex;}
-
- double GetTheta1() {return Theta_1 ;}
- double GetTheta2() {return Theta_2 ;}
-
- double GetPhi1(){return Phi1;}
- double GetPhi2(){return Phi2;}
-
- double GetDmin(){return Dmin;}
- double GetTheta_12(){return Theta_12;}
-
- int GetNbrOfTracks(){return trackNbr_FINAL;}
-
- TVector3 GetVectorTrack1(){return VectorTrack1;}
- TVector3 GetVectorTrack2(){return VectorTrack2;}
- // parameters used in the analysis
- private:
-
-
- bool SimulationBool;//!
-
- // thresholds
- double m_E_RAW_Threshold; //!
- double m_E_Threshold; //!
-
- // used in PreTreatedData()
- vector<int> clusterringboolTemp;//!
- vector<int> clusterringbool;//!
- vector<int> clusternbr;//!
- vector<int> clusterpads;//!
-
- vector<double> Xpad;//!
- vector<double> Ypad;//!
- vector<double> Qpad;//!
- vector<double> XpadNew;
- vector<double> YpadNew;
- vector<double> QpadNew;//!
- vector<double> ZpadNew;//!
-
- double MINOSthresh, UperFitLim, VDrift, Z_Shift, DelayTrig[18]={0},VdriftperRing[18]={0}, ZRot_Minos, coef_Vdrift; //!
-
- double x_mm;//!
- double y_mm;//!
- double z_mm;//!
- double q_pad;//!
- double t_pad;//!
- double Chi2;//!
- double maxCharge;//!
- double ChargeBin;//!
-
- double hfit_max; //!
- double hfit_max_T; //!
- double T_min; //!
- double T_max; //!
-
-
- int Iteration;///
- int filter_result;//!
- int filled;//!
- int fit2DStatus;//!
- int indexfill;//!
- int trackNbr;//!
-
- TF1 *fit_function;//!
- /* TH1F *hfit = new TH1F("hfit","hfit",512,0,512);//! */
- TH1F *hfit ;//!
-
- int npoint_temp=0;//!
- int cluster_temp;//!
- int ringsum;//!
- int ringtouch[19];//!
-
- TGraph *gryztmp;//!
- TGraph* grxztmp;//!
-
- double zmax;//!
-
- vector<double> TOTxoutprime;
- vector<double> TOTyoutprime;
- vector<double> TOTzoutprime;
- vector<double> TOTqout;
-
- double Xvertex;
- double Yvertex;
- double Zvertex;
- double Dmin;
-
- double Theta_12;
-
- double Theta_1;
- double Theta_2;
-
- double Phi1;
- double Phi2;
-
- TVector3 VectorTrack1, VectorTrack2;
-
- int trackNbr_FINAL;
-
- vector<int> trackclusternbr;//!
- vector<int> tracknbr;//!
-
- double Tot_tracks=0; //!
- vector<TGraph> gryz;//!
- vector<TGraph> grxz;//!
-
- int npoint;//!
- int array_final;//!
- TVector3 point;//!
-
- // Variables only filled when trackNbr_FINAL>=1
- int allevt_2pfiltered;//!
- Double_t pStart[4];//!
- double parFit_temp[4];//!
- double err_temp[4];//!
- Double_t chi[2];//!
- Int_t fitStatus[2];//!
- Double_t arglist[10];//!
- Int_t iflag;//!
- int nvpar;//!
- int nparx;//!
- double amin;//!
- double edm;//!
- double errdef;//!
-
- vector<double> lenght;
- vector<double> chargeTot;
- vector<double> parFit1;
- vector<double> parFit2;
- vector<double> parFit3;
- vector<double> parFit4;
- vector<double> errFit1_local;
- vector<double> errFit2_local;
- vector<double> errFit3_local;
- vector<double> errFit4_local;
+ TMinosData* GetRawData() const {return m_EventData;}//!
+ vector<double> GetPad_X() {return X_Pad;} //!
+ vector<double> GetPad_Y() {return Y_Pad;} //!
+ vector<double> GetPad_Z() {return Z_Pad;} //!
+ vector<double> GetPad_T() {return T_Pad;} //!
+ double GetVertexX() {return X_Vertex;} //!
+ double GetVertexY() {return Y_Vertex;} //!
+ double GetVertexZ() {return Z_Vertex;} //!
- // number of detectors
- private:
- int m_NumberOfDetectors; //!
-
+ int GetNbrOfTracks(){return Tracks_P0.size();}
+
+ TVector3 GetTracksP0(unsigned int i){return Tracks_P0[i];}//!
+ TVector3 GetTracksDir(unsigned int i){return Tracks_Dir[i];}//!
+ double Angle(unsigned int i, unsigned j){return Tracks_Dir[i].Angle(Tracks_Dir[j]);};
// Static constructor to be passed to the Detector Factory
public:
static NPL::VDetector* Construct();
diff --git a/NPLib/TrackReconstruction/CMakeLists.txt b/NPLib/TrackReconstruction/CMakeLists.txt
index 266e1c3a8faf0b89640707926528a2cb9cff1814..341a03e345feb5010c5bc4220dab7d4df0449fe3 100644
--- a/NPLib/TrackReconstruction/CMakeLists.txt
+++ b/NPLib/TrackReconstruction/CMakeLists.txt
@@ -1,3 +1,5 @@
+add_custom_command(OUTPUT NPLinearRansac3DDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh NPLinearRansac3D.h NPLinearRansac3DDict.cxx NPLinearRansac3D.rootmap libNPTrackReconstruction.so NPTrackReconstructionLinkDef.h DEPENDS NPLinearRansac3D.h)
+
add_custom_command(OUTPUT NPRansacDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh NPRansac.h NPRansacDict.cxx NPRansac.rootmap libNPTrackReconstruction.so NPTrackReconstructionLinkDef.h DEPENDS NPRansac.h)
add_custom_command(OUTPUT NPClusterDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh NPCluster.h NPClusterDict.cxx NPCluster.rootmap libNPTrackReconstruction.so NPTrackReconstructionLinkDef.h DEPENDS NPCluster.h)
@@ -14,7 +16,7 @@ else()
endif()
-add_library(NPTrackReconstruction SHARED NPRansac.cxx NPCluster.cxx NPTrack.cxx Tracking.cxx NPRansacDict.cxx NPClusterDict.cxx TrackingDict.cxx NPDCReconstruction.cxx NPDCReconstructionMT.cxx)
+add_library(NPTrackReconstruction SHARED NPTrackingUtility.cxx NPLinearRansac3D.cxx NPRansac.cxx NPCluster.cxx NPTrack.cxx Tracking.cxx NPLinearRansac3DDict.cxx NPRansacDict.cxx NPClusterDict.cxx TrackingDict.cxx NPDCReconstruction.cxx NPDCReconstructionMT.cxx)
if(Minuit2_FOUND)
target_link_libraries(NPTrackReconstruction ${ROOT_LIBRARIES} -lMinuit2 NPCore)
@@ -22,4 +24,4 @@ else()
target_link_libraries(NPTrackReconstruction ${ROOT_LIBRARIES} NPCore)
endif()
-install(FILES NPRansac.h NPCluster.h NPTrack.h Tracking.h NPTrackingUtility.h NPDCReconstruction.h NPDCReconstructionMT.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
+install(FILES NPLinearRansac3D.h NPRansac.h NPCluster.h NPTrack.h Tracking.h NPTrackingUtility.h NPDCReconstruction.h NPDCReconstructionMT.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
diff --git a/NPLib/TrackReconstruction/NPDCReconstructionMT.cxx b/NPLib/TrackReconstruction/NPDCReconstructionMT.cxx
index 2d466e04b4630f202d570a9e699948c625c4d22e..176f0543346fe956c7f0f6e78d1702c20e3ad3d9 100644
--- a/NPLib/TrackReconstruction/NPDCReconstructionMT.cxx
+++ b/NPLib/TrackReconstruction/NPDCReconstructionMT.cxx
@@ -149,11 +149,15 @@ double DCReconstructionMT::SumD(const double* parameter ){
double a2=a*a;
unsigned int id = parameter[2];
unsigned int size = sizeX[id];
+ const std::vector<double>* X=fitX[id];
+ const std::vector<double>* Z=fitZ[id];
+ const std::vector<double>* R=fitR[id];
+
double c,d,r,x,z,p;
for(unsigned int i = 0 ; i < size ; i++){
- c = (*fitX[id])[i];
- d = (*fitZ[id])[i];
- r = (*fitR[id])[i];
+ c = (*X)[i];
+ d = (*Z)[i];
+ r = (*R)[i];
x = (a*d-ab+c)/(1+a2);
z = a*x+b;
p= (x-c)*(x-c)+(z-d)*(z-d)-r*r;
diff --git a/NPLib/TrackReconstruction/NPTrackReconstructionLinkDef.h b/NPLib/TrackReconstruction/NPTrackReconstructionLinkDef.h
index 8fc0faa4b11980432297877cd4cd6fcb24012501..8c516b628e0f0fad4be5b1b60399bddd1180d351 100644
--- a/NPLib/TrackReconstruction/NPTrackReconstructionLinkDef.h
+++ b/NPLib/TrackReconstruction/NPTrackReconstructionLinkDef.h
@@ -1,5 +1,6 @@
#ifdef __CINT__
#pragma link C++ defined_in "./NPRansac.h";
+#pragma link C++ defined_in "./NPLinearRansac3D.h";
#pragma link C++ defined_in "./NPCluster.h";
#pragma link C++ defined_in "./Tracking.h";
#endif
diff --git a/NPLib/TrackReconstruction/NPTrackingUtility.h b/NPLib/TrackReconstruction/NPTrackingUtility.h
index ee7c30502399c3bbae7ac174869277243d923457..5fc82d647d80e03065b96228165a947d741e8842 100644
--- a/NPLib/TrackReconstruction/NPTrackingUtility.h
+++ b/NPLib/TrackReconstruction/NPTrackingUtility.h
@@ -1,5 +1,5 @@
-#ifndef NPUTILITY_H
-#define NPUTILITY_H
+#ifndef NPTrackingUtility_H
+#define NPTrackingUtility_H
/*****************************************************************************
* Copyright (C) 2009-2016 this file is part of the NPTool Project *
* *
@@ -16,38 +16,20 @@
* Decription: *
* This class deal with finding all the track event by event *
*****************************************************************************/
-
+#include "TVector3.h"
namespace NPL{
//////////////////////////////////////////////////////////////////////////////
// return the minimum distance between v and w defined respectively by points
// v1,v2 and w1 w1
// Also compute the best crossing position BestPosition, i.e. average position
// at the minimum distance.
- double MinimumDistanceTwoLines(const TVector3& v1,const TVector3& v2, const TVector3& w1, const TVector3& w2, TVector3& BestPosition, TVector3& delta){
- TVector3 v = v2-v1;
- TVector3 w = w2-w1;
- // Finding best position
- TVector3 e = v1-w1;
- double A = -(v.Mag2()*w.Mag2()-(v.Dot(w)*v.Dot(w)));
- double s = (-v.Mag2()*(w.Dot(e))+(v.Dot(e))*(w.Dot(v)))/A;
- double t = (w.Mag2()*(v.Dot(e))-(w.Dot(e)*w.Dot(v)))/A;
- double d = sqrt((e+v*t-w*s).Mag2());
-
- BestPosition = 0.5*(v1+t*v+w1+s*w);
- delta = (v1+t*v-w1-s*w);
- return d;
- }
+ double MinimumDistanceTwoLines(const TVector3& v1,const TVector3& v2, const TVector3& w1, const TVector3& w2, TVector3& BestPosition, TVector3& delta);
//////////////////////////////////////////////////////////////////////////////
// return the minimum distance between the line defines by v1,v2 and the point
// in space x
// demo is here: https://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html
- double MinimumDistancePointLine(const TVector3& v1, const TVector3& v2, const TVector3& x){
- TVector3 w1 = x-v1;
- TVector3 w2 = x-v2;
- TVector3 w = w1.Cross(w2);
- return w.Mag()/(v2-v1).Mag();
- }
+ double MinimumDistancePointLine(const TVector3& v1, const TVector3& v2, const TVector3& x);
}
diff --git a/NPSimulation/Detectors/Minos/Minos.cc b/NPSimulation/Detectors/Minos/Minos.cc
index 38a58decd0fe6476fa5193fcdeb8bb7d6d6a6c6d..b3589f98c4833419851fcbfe1657fc38ea8c42c5 100644
--- a/NPSimulation/Detectors/Minos/Minos.cc
+++ b/NPSimulation/Detectors/Minos/Minos.cc
@@ -7,6 +7,7 @@
/*****************************************************************************
* Original Author: E. Tronchin contact address: tronchin@lpccaen.in2p3.fr *
+ * Updated by C. Lenain contact address: lenain@lpccaen.in2p3.fr *
* *
* Creation Date : October 2018 *
* Last update : *
@@ -80,19 +81,19 @@ using namespace CLHEP;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
namespace Minos_NS{
-// TPC
-const G4double ChamberInnerRadius = 37.*mm;
-/* const G4double ChamberInnerRadius = 29*mm; */ //big TPC
-//const G4double ChamberThickness = 2.*mm;
-const G4double ChamberLength = 300.*mm;
-const G4double KaptonThickness = 0.125*mm;
-const G4double RohacellThickness = 2.*mm;
-const G4double TPCRadiusExt = 91.525*mm;
-/* const G4double TPCRadiusExt = 150*mm; //big TPC */
+ // TPC
+ const G4double ChamberInnerRadius = 37.*mm;
+ /* const G4double ChamberInnerRadius = 29*mm; */ //big TPC
+ //const G4double ChamberThickness = 2.*mm;
+ const G4double ChamberLength = 300.*mm;
+ const G4double KaptonThickness = 0.125*mm;
+ const G4double RohacellThickness = 2.*mm;
+ const G4double TPCRadiusExt = 91.525*mm;
+ /* const G4double TPCRadiusExt = 150*mm; //big TPC */
-// MINOS
-const G4double TargetRadius = 28.*mm;
-const G4double WindowThickness = 0.150*mm;
+ // MINOS
+ const G4double TargetRadius = 28.*mm;
+ const G4double WindowThickness = 0.150*mm;
}
@@ -102,41 +103,42 @@ using namespace Minos_NS;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Minos Specific Method
Minos::Minos(){
-m_Event = new TMinosData() ;
-m_MinosPadScorer = 0;
-m_ReactionRegion=NULL;
-
-// RGB Color + Transparency
-m_VisTarget= new G4VisAttributes(G4Colour(0.6,1.,1., .4));
-m_VissimpleBox= new G4VisAttributes(G4Colour(0,1,0,.6));
-m_VisTPC= new G4VisAttributes(G4Colour(1.,0.5,0.6,0.3));
-m_VisRohacell= new G4VisAttributes(G4Colour(1.,1.,1., .8));
-m_VisKapton = new G4VisAttributes(G4Colour(1.,1.,0.6,0.4));
-m_VisTargetCell = new G4VisAttributes(G4Colour(0,0,1, .4));
-m_VisTargetCell->SetForceSolid(true);
-m_VisOuterKapton = new G4VisAttributes(G4Colour(1.,1.,0.6,0.8));
-
-Raw_Signal = new TH1F("raw_Signal","raw_Signal",350,0,350); // 30ns per bin
-Elec_Signal = new TH1F("Elec_Signal","Elec_Signal",350,0,350);
-
-/* Raw_Signal = new TH1F;*/
-/* Elec_Signal = new TH1F;*/
-
-
-fa1 = new TF1("fa1","[0]*exp(-3.*(x-[1])/[2]) * sin((x-[1])/[2]) * pow((x-[1])/[2], 3)",0,1000);
-
-solidTarget=0;
-logicTarget=0;
-solidChamber=0;
-logicChamber=0;
-solidTPC=0;
-logicTPC=0;
-solidWindow0=0;
-logicWindow0=0;
-solidRohacell=0;
-logicRohacell=0;
-solidKapton=0;
-logicKapton=0;
+ m_Event = new TMinosData() ;
+ m_MinosPadScorer = 0;
+ m_ReactionRegion=NULL;
+
+ // RGB Color + Transparency
+ m_VisTarget= new G4VisAttributes(G4Colour(0.6,1.,1., .4));
+ m_VissimpleBox= new G4VisAttributes(G4Colour(0,1,0,.6));
+ m_VisTPC= new G4VisAttributes(G4Colour(1.,0.5,0.6,0.3));
+ m_VisRohacell= new G4VisAttributes(G4Colour(1.,1.,1., .8));
+ m_VisKapton = new G4VisAttributes(G4Colour(1.,1.,0.6,0.4));
+ m_VisTargetCell = new G4VisAttributes(G4Colour(0,0,1, .4));
+ m_VisTargetCell->SetForceSolid(true);
+ m_VisOuterKapton = new G4VisAttributes(G4Colour(1.,1.,0.6,0.8));
+
+ Raw_Signal = new TH1F("raw_Signal","raw_Signal",512,0,512);
+ Elec_Signal = new TH1F("Elec_Signal","Elec_Signal",512,0,512);
+
+ /* Raw_Signal = new TH1F;*/
+ /* Elec_Signal = new TH1F;*/
+
+
+ fa1 = new TF1("fa1","abs((x>[1]&&x<512)*([0]*exp(-3.*(x-[1])/[2]) * sin((x-[1])/[2]) * pow((x-[1])/[2],3))+[3])",0,1000);
+ fa1->SetNpx(512);
+
+ solidTarget=0;
+ logicTarget=0;
+ solidChamber=0;
+ logicChamber=0;
+ solidTPC=0;
+ logicTPC=0;
+ solidWindow0=0;
+ logicWindow0=0;
+ solidRohacell=0;
+ logicRohacell=0;
+ solidKapton=0;
+ logicKapton=0;
}
Minos::~Minos(){
@@ -147,11 +149,11 @@ Minos::~Minos(){
/* void Minos::AddDetector(G4ThreeVector POS, double LengthOfTarget, int PresenceOfMinos){ */
void Minos::AddDetector(G4ThreeVector POS, double LengthOfTarget, G4String MaterialOfTarget,G4String MaterialOfCell, int PresenceOfMinos){
-m_POS.push_back(POS);
-m_TargetLength.push_back(LengthOfTarget);
-m_TargetMaterial.push_back(MaterialOfTarget);
-m_CellMaterial.push_back(MaterialOfCell);
-m_TPCOnly.push_back(PresenceOfMinos);
+ m_POS.push_back(POS);
+ m_TargetLength.push_back(LengthOfTarget);
+ m_TargetMaterial.push_back(MaterialOfTarget);
+ m_CellMaterial.push_back(MaterialOfCell);
+ m_TPCOnly.push_back(PresenceOfMinos);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -294,21 +296,29 @@ void Minos::ReadConfiguration(NPL::InputParser parser){
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << "//// " << blocks.size() << " detectors found " << endl;
- vector<string> cart = {"POS","TargetLength"};
+ vector<string> optional = {"TPCOnly"};
+
+ vector<string> token= {"Position","TargetMaterial","TargetLength","CellMaterial","TimeBin","ShapingTime","Baseline","Sampling","ZOffset"};
for(unsigned int i = 0 ; i < blocks.size() ; i++){
- if(blocks[i]->HasTokenList(cart)){
+ if(blocks[i]->HasTokenList(token)){
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << endl << "//// Minos " << i+1 << endl;
- G4ThreeVector Pos = NPS::ConvertVector(blocks[i]->GetTVector3("POS","mm"));
+ G4ThreeVector Pos = NPS::ConvertVector(blocks[i]->GetTVector3("Position","mm"));
TargetLength = blocks[i]->GetDouble("TargetLength","mm");
G4String TargetMaterialname = blocks[i]->GetString("TargetMaterial");
G4String CellMaterial = blocks[i]->GetString("CellMaterial");
- TPCOnly = blocks[i]->GetInt("TPCOnly");
+ m_ShapingTime = blocks[i]->GetDouble("ShapingTime","ns")/ns;
+ m_TimeBin = blocks[i]->GetDouble("TimeBin","ns")/ns;
+ m_Sampling= blocks[i]->GetInt("Sampling");
+ m_Baseline= blocks[i]->GetInt("BaseLine");
+ m_ZOffset = blocks[i]->GetDouble("ZOffset","mm");
+ TPCOnly=1;
+ if(blocks[i]->HasTokenList(optional))
+ TPCOnly = blocks[i]->GetInt("TPCOnly");
AddDetector(Pos,TargetLength,TargetMaterialname, CellMaterial, TPCOnly);
/* AddDetector(Pos,TargetLength, TPCOnly); */
-
}
else{
cout << "ERROR: check your input file formatting " << endl;
@@ -426,7 +436,7 @@ void Minos::ConstructDetector(G4LogicalVolume* world){
Cu,"div_ringL",0,0,0);
{G4VisAttributes* atb= new G4VisAttributes(G4Colour(0.8, 0.4, 0.,0.4));
- Pad_log->SetVisAttributes(atb);}
+ Pad_log->SetVisAttributes(atb);}
Pad_log->SetSensitiveDetector(m_MinosPadScorer);
new G4PVReplica("div_ring_phys", Pad_log,
@@ -531,28 +541,18 @@ void Minos::ReadSensitive(const G4Event* ){
///////////
// DriftElectron scorer
CylinderTPCScorers::PS_TPCAnode* Scorer2= (CylinderTPCScorers::PS_TPCAnode*) m_MinosPadScorer->GetPrimitive(0);
-
unsigned int size2 = Scorer2->GetMult();
for(unsigned int i = 0 ; i < size2 ; i++){
-
int Pad = Scorer2->GetPad(i);
- double X = Scorer2->GetX(i);
- double Y = Scorer2->GetY(i);
-
- m_Event->SetCharge(Pad,X,Y);
- //m_Event->AddChargePoint(Scorer2->GetQ(i), Scorer2->GetT(i));
-
- SimulateGainAndDigitizer(Scorer2->GetQ(i), Scorer2->GetT(i));
-
- /* m_Event->AddChargePoint(Scorer2->GetQ(i),Scorer2->GetT(i)); */
+ SimulateGainAndDigitizer(Scorer2->GetT(i), Scorer2->GetQ(i),T,Q);
+ m_Event->SetPad(Pad, Scorer2->GetQ(i)->size(),&T,&Q);
}
}
-void Minos::SimulateGainAndDigitizer(vector<double> rawQ, vector<double> rawT){
-
- Charge2.clear();
- Time.clear();
+void Minos::SimulateGainAndDigitizer(vector<double>* rawT, vector<double>* rawQ,vector<int>& T, vector<int>& Q){
+ T.clear();
+ Q.clear();
Raw_Signal->Reset();
Elec_Signal->Reset();
@@ -576,38 +576,35 @@ void Minos::SimulateGainAndDigitizer(vector<double> rawQ, vector<double> rawT){
/* } */
/* } */
- for ( unsigned int j = 0; j < rawQ.size(); j++){
- for ( int i = 0 ; i < rawQ[j]; i++ ){
- Raw_Signal->Fill((rawT[j])/30.);
+ for ( unsigned int j = 0; j < rawQ->size(); j++){
+ for ( int i = 0 ; i < (*rawQ)[j]; i++ ){
+ Raw_Signal->Fill(((*rawT)[j])/m_TimeBin);
}
}
// Application of the electronic reponse function
for ( int x=0; x < Raw_Signal->GetNbinsX(); x++){
-
- if(Raw_Signal->GetBinContent(x) < 0.5)
- continue;
- else{
+ if(Raw_Signal->GetBinContent(x) > 0.5){
start = Raw_Signal->GetBinCenter(x);
- end = Raw_Signal->GetBinCenter(x)+1200/30.;
+ end = Raw_Signal->GetBinCenter(x)+512;
// DriftTime = Raw_Signal->GetBinCenter(x);
fa1->SetRange(start, end);
fa1->SetParameter(0,1500);
fa1->SetParameter(1,start);
- fa1->SetParameter(2,450/30.);
- // Elec_Signal->SetBinContent(bin,20+(10-gRandom->Uniform(20)*20))
- // for (int p=0; p< 0; p++)
- for (int p=0; p< Raw_Signal->GetBinContent(x)*1500; p++)
- Elec_Signal->Fill(fa1->GetRandom());
- }
- }
+ fa1->SetParameter(2,m_ShapingTime/(log(2)*m_TimeBin));
+ fa1->SetParameter(3,0);
- for ( int bin=0; bin < Elec_Signal->GetNbinsX(); bin++){
- Charge2.push_back(Elec_Signal->GetBinContent(bin));
- Time.push_back(Elec_Signal->GetBinCenter(bin));
+ for (int p=0; p < Raw_Signal->GetBinContent(x)*1500; p++)
+ Elec_Signal->Fill(fa1->GetRandom(start,end));
+ }
}
- m_Event->AddChargePoint(Charge2,Time);
+ for ( int bin=0; bin < Elec_Signal->GetNbinsX(); bin++){
+ if(Elec_Signal->GetBinContent(bin)){
+ Q.push_back(Elec_Signal->GetBinContent(bin)+m_Baseline);
+ T.push_back(Elec_Signal->GetBinCenter(bin));
+ }
+ }
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/NPSimulation/Detectors/Minos/Minos.hh b/NPSimulation/Detectors/Minos/Minos.hh
index 71ca4078b5e99e9a2fd74708da22646b7250eb62..c32ba9c7caff33326f13a3288548e06c574ce552 100644
--- a/NPSimulation/Detectors/Minos/Minos.hh
+++ b/NPSimulation/Detectors/Minos/Minos.hh
@@ -134,7 +134,7 @@ class Minos : public NPS::VDetector{
TH1F* Raw_Signal ;
TH1F* Elec_Signal;
TF1* fa1;
- vector<double> Charge2, Time;
+ vector<int> Q, T;
////////////////////////////////////////////////////
////// Inherite from NPS::VDetector class /////////
@@ -159,8 +159,16 @@ class Minos : public NPS::VDetector{
public: // Scorer
// Initialize all Scorer used by the MUST2Array
void InitializeScorers() ;
- void SimulateGainAndDigitizer(vector<double> Q, vector<double> T);
-
+ void SimulateGainAndDigitizer(vector<double>* rawT, vector<double>* rawQ,vector<int>& Q,vector<int>& T);
+
+ private: // parameter for the digitization
+ double m_TimeBin;
+ double m_ShapingTime;
+ double m_Baseline;
+ unsigned int m_Sampling;
+ double m_ZOffset;
+
+
// Associated Scorer
G4MultiFunctionalDetector* m_MinosPadScorer ;
diff --git a/NPSimulation/Scorers/CylinderTPCScorers.hh b/NPSimulation/Scorers/CylinderTPCScorers.hh
index 1f9434d196fdc7ab05f265c9f38eb07ddb0f24ae..f04c9c1ff6ad897db4b4b27d08a7d63eaeb461da 100644
--- a/NPSimulation/Scorers/CylinderTPCScorers.hh
+++ b/NPSimulation/Scorers/CylinderTPCScorers.hh
@@ -62,8 +62,8 @@ namespace CylinderTPCScorers{
inline double GetY() const {return m_Y;};
inline unsigned int GetPad() const {return m_Pad;};
- inline vector<double> GetQ() const {return m_Q;};
- inline vector<double> GetT() const {return m_T;};
+ inline vector<double>* GetQ() {return &m_Q;};
+ inline vector<double>* GetT() {return &m_T;};
void Add(const double& Charge){m_Charge+=Charge;};
@@ -139,8 +139,8 @@ namespace CylinderTPCScorers{
double GetY(const double& i) {return m_Data[i]->GetY();};
unsigned int GetPad(const unsigned int& i) {return m_Data[i]->GetPad();};
- vector<double> GetQ(const double& i) {return m_Data[i]->GetQ();}
- vector<double> GetT(const double& i) {return m_Data[i]->GetT();}
+ vector<double>* GetQ(const double& i) {return m_Data[i]->GetQ();}
+ vector<double>* GetT(const double& i) {return m_Data[i]->GetT();}
};
diff --git a/Projects/Dali/DaliMinos.detector b/Projects/Dali/DaliMinos.detector
index ccb30275e9f2d0f8d62194602ee9f4562c2e4eab..e25d6ed3a4b9892537f35e6f96737eb39554a619 100755
--- a/Projects/Dali/DaliMinos.detector
+++ b/Projects/Dali/DaliMinos.detector
@@ -4,6 +4,7 @@ Minos
TargetMaterial= LH2
CellMaterial= Mylar
TPCOnly= 0
+ XML= /Users/matta/mrdc/xml/s034/MINOS.xml
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Dali
R = 212.4 mm
diff --git a/Projects/Minos/Minos.detector b/Projects/Minos/Minos.detector
index 3b4256f5c4243c0be04d643017b34ab736dba346..56ee4a30224a0cca43d4db825f4aa35abd17e31a 100644
--- a/Projects/Minos/Minos.detector
+++ b/Projects/Minos/Minos.detector
@@ -9,8 +9,14 @@
% Z= 75 mm
Minos
- POS= 0 0 0 mm
+ Position= 0 0 0 mm
TargetLength= 152.76 mm
TargetMaterial= LH2
CellMaterial= Mylar
TPCOnly= 0
+ TimeBin= 30 ns
+ ShapingTime= 333.9 ns
+ BaseLine= 250
+ Sampling= 10
+ ZOffset= 0 mm
+ XML= /Users/matta/mrdc/xml/s034/MINOS.xml
diff --git a/Projects/Minos/PhysicsListOption.txt b/Projects/Minos/PhysicsListOption.txt
index 68f97e860e04654d462d3f6e1888f1ab694af8de..5e948e6061566e703143f8c5a86e4ee5ff6d7f76 100644
--- a/Projects/Minos/PhysicsListOption.txt
+++ b/Projects/Minos/PhysicsListOption.txt
@@ -1,5 +1,5 @@
EmPhysicsList Option4
-DriftElectronPhysics 0
+DriftElectronPhysics 1
DefaultCutOff 1e9
IonBinaryCascadePhysics 0
NPIonInelasticPhysics 0