From 38aaf777c00c670198eaf6325e29b378001e8a42 Mon Sep 17 00:00:00 2001
From: Morfouace <morfouac@ipno.in2p3.fr>
Date: Tue, 28 Aug 2018 10:26:06 +0200
Subject: [PATCH] ! Adding Actar files for analysis and simulation ! Adding a
new folder TrackReconstruction for Tracking algorithm
---
Inputs/EventGenerator/18Opp.reaction | 2 +-
Inputs/EventGenerator/proton.source | 6 +-
NPLib/CMakeLists.txt | 21 +-
NPLib/Detectors/Actar/CMakeLists.txt | 7 +-
NPLib/Detectors/Actar/MEventReduced.h | 40 +
NPLib/Detectors/Actar/MTreeStructureLinkDef.h | 12 +
NPLib/Detectors/Actar/TActarData.cxx | 23 +-
NPLib/Detectors/Actar/TActarData.h | 60 +-
NPLib/Detectors/Actar/TActarPhysics.cxx | 601 +++++++++-----
NPLib/Detectors/Actar/TActarPhysics.h | 206 +++--
NPLib/Detectors/Actar/TActarSpectra.cxx | 24 +-
NPLib/TrackReconstruction/CMakeLists.txt | 6 +
NPLib/TrackReconstruction/NPRansac.cxx | 774 ++++++++++++++++++
NPLib/TrackReconstruction/NPRansac.h | 117 +++
NPLib/TrackReconstruction/NPTrack.cxx | 210 +++++
NPLib/TrackReconstruction/NPTrack.h | 96 +++
.../NPTrackReconstructionLinkDef.h | 4 +
NPLib/scripts/build_dict.sh | 4 +-
NPSimulation/Core/MaterialManager.cc | 16 +-
NPSimulation/Core/Target.cc | 4 +-
NPSimulation/Detectors/Actar/Actar.cc | 703 ++++++++++------
NPSimulation/Detectors/Actar/Actar.hh | 37 +-
.../EventGenerator/EventGeneratorIsotropic.cc | 7 +-
.../EventGenerator/EventGeneratorIsotropic.hh | 1 +
NPSimulation/Process/BeamReaction.cc | 4 +-
NPSimulation/Scorers/CMakeLists.txt | 2 +-
NPSimulation/Scorers/TPCScorers.cc | 119 +++
NPSimulation/Scorers/TPCScorers.hh | 94 +++
Projects/Actar/Analysis.cxx | 168 +++-
Projects/Actar/Analysis.h | 60 +-
30 files changed, 2745 insertions(+), 683 deletions(-)
create mode 100644 NPLib/Detectors/Actar/MEventReduced.h
create mode 100644 NPLib/Detectors/Actar/MTreeStructureLinkDef.h
create mode 100644 NPLib/TrackReconstruction/CMakeLists.txt
create mode 100644 NPLib/TrackReconstruction/NPRansac.cxx
create mode 100644 NPLib/TrackReconstruction/NPRansac.h
create mode 100644 NPLib/TrackReconstruction/NPTrack.cxx
create mode 100644 NPLib/TrackReconstruction/NPTrack.h
create mode 100644 NPLib/TrackReconstruction/NPTrackReconstructionLinkDef.h
create mode 100644 NPSimulation/Scorers/TPCScorers.cc
create mode 100644 NPSimulation/Scorers/TPCScorers.hh
diff --git a/Inputs/EventGenerator/18Opp.reaction b/Inputs/EventGenerator/18Opp.reaction
index 488eed07e..13be99a5f 100644
--- a/Inputs/EventGenerator/18Opp.reaction
+++ b/Inputs/EventGenerator/18Opp.reaction
@@ -4,7 +4,7 @@
%%Beam energy given in MeV ; Excitation in MeV
Beam
Particle= 18O
- Energy= 137.88
+ Energy= 59.4
SigmaEnergy= 0
SigmaThetaX= 0
SigmaPhiY= 0
diff --git a/Inputs/EventGenerator/proton.source b/Inputs/EventGenerator/proton.source
index 7b62f807a..7a4e1e53e 100644
--- a/Inputs/EventGenerator/proton.source
+++ b/Inputs/EventGenerator/proton.source
@@ -4,10 +4,10 @@
% Energy are given in MeV , Position in mm %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Isotropic
- EnergyLow= 30
- EnergyHigh= 80
+ EnergyLow= 20
+ EnergyHigh= 25
HalfOpenAngleMin= 0
- HalfOpenAngleMax= 1
+ HalfOpenAngleMax= 180
x0= 0
y0= 0
z0= 0
diff --git a/NPLib/CMakeLists.txt b/NPLib/CMakeLists.txt
index fd1349561..0f10c8e49 100644
--- a/NPLib/CMakeLists.txt
+++ b/NPLib/CMakeLists.txt
@@ -2,6 +2,7 @@ cmake_minimum_required (VERSION 2.8)
include(CheckCXXCompilerFlag)
project(NPLib CXX)
set(CMAKE_BUILD_TYPE Release)
+
# Setting the policy to match Cmake version
cmake_policy(VERSION ${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION})
@@ -18,23 +19,6 @@ set(NPLIB_VERSION_DETA 45)
configure_file(Core/NPLibVersion.h.in Core/NPLibVersion.h @ONLY)
set(DETLIST ${ETLIST})
-#activate Multithreading (on by default)
-if(NOT DEFINED NPMULTITHREADING)
- set(NPMULTITHREADING 1)
- else()
- set(NPMULTITHREADING ${NPMULTITHREADING})
-endif()
-
-if(NPMULTITHREADING)
- set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DNPMULTITHREADING=1")
- message("Building application with MultiThreading Support")
- else()
- set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DNPMULTITHREADING=0")
- message("Building application with no MutilThreading Support")
-endif()
-
-
-
string(COMPARE EQUAL "${DETLIST}" "" rdet)
if(rdet)
message("Building all detectors")
@@ -73,11 +57,12 @@ ENDMACRO()
# Add include directories to all target
include_directories("Core/")
include_directories("Physics/")
+include_directories("TrackReconstruction/")
include_directories("Online/")
# Call the macro
subdirlist(SUB_DIRECTORY ${CMAKE_BINARY_DIR})
-set(SUB_DIRECTORY ${SUB_DIRECTORY} Core Physics Calibration Online Utility)
+set(SUB_DIRECTORY ${SUB_DIRECTORY} Core Physics TrackReconstruction Calibration Online Utility)
# Add each sub folder to the project
set(TARGET_LIST "")
foreach(subdir ${SUB_DIRECTORY})
diff --git a/NPLib/Detectors/Actar/CMakeLists.txt b/NPLib/Detectors/Actar/CMakeLists.txt
index a1b027f69..3a0cfd553 100644
--- a/NPLib/Detectors/Actar/CMakeLists.txt
+++ b/NPLib/Detectors/Actar/CMakeLists.txt
@@ -1,6 +1,7 @@
add_custom_command(OUTPUT TActarPhysicsDict.cxx COMMAND ../../scripts/build_dict.sh TActarPhysics.h TActarPhysicsDict.cxx TActarPhysics.rootmap libNPActar.dylib DEPENDS TActarPhysics.h)
add_custom_command(OUTPUT TActarDataDict.cxx COMMAND ../../scripts/build_dict.sh TActarData.h TActarDataDict.cxx TActarData.rootmap libNPActar.dylib DEPENDS TActarData.h)
-add_library(NPActar SHARED TActarSpectra.cxx TActarData.cxx TActarPhysics.cxx TActarDataDict.cxx TActarPhysicsDict.cxx )
-target_link_libraries(NPActar ${ROOT_LIBRARIES} NPCore)
-install(FILES TActarData.h TActarPhysics.h TActarSpectra.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
+add_custom_command(OUTPUT MEventReducedDict.cxx COMMAND ../../scripts/build_dict.sh MEventReduced.h MEventReducedDict.cxx MEventReduced.rootmap libNPActar.dylib MTreeStructureLinkDef.h DEPENDS MEventReduced.h)
+add_library(NPActar SHARED TActarSpectra.cxx TActarData.cxx TActarPhysics.cxx TActarDataDict.cxx TActarPhysicsDict.cxx MEventReducedDict.cxx)
+target_link_libraries(NPActar ${ROOT_LIBRARIES} NPCore NPTrackReconstruction)
+install(FILES TActarData.h TActarPhysics.h TActarSpectra.h MEventReduced.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
diff --git a/NPLib/Detectors/Actar/MEventReduced.h b/NPLib/Detectors/Actar/MEventReduced.h
new file mode 100644
index 000000000..0d642c10a
--- /dev/null
+++ b/NPLib/Detectors/Actar/MEventReduced.h
@@ -0,0 +1,40 @@
+#ifndef MEVENTREDUCED_H
+#define MEVENTREDUCED_H
+
+#include <stdio.h>
+#include <iostream>
+#include <cstdlib>
+#include <vector>
+#include <TObject.h>
+
+using namespace std;
+
+class ReducedData: public TObject
+{
+ public:
+ ReducedData(){};
+ virtual ~ReducedData(){};
+
+ unsigned short globalchannelid;
+ bool hasSaturation;
+
+ std::vector<float> peakheight;
+ std::vector<float> peaktime;
+
+ ClassDef(ReducedData,1);
+};
+
+class MEventReduced: public TObject
+{
+ public:
+ MEventReduced(){};
+ virtual ~MEventReduced(){};
+ unsigned long int event;
+ unsigned long int timestamp;
+ std::vector<ReducedData> CoboAsad;
+
+ ClassDef(MEventReduced,1);
+};
+
+
+#endif
diff --git a/NPLib/Detectors/Actar/MTreeStructureLinkDef.h b/NPLib/Detectors/Actar/MTreeStructureLinkDef.h
new file mode 100644
index 000000000..c9ec8b317
--- /dev/null
+++ b/NPLib/Detectors/Actar/MTreeStructureLinkDef.h
@@ -0,0 +1,12 @@
+#ifdef __CINT__
+
+#pragma link off all globals;
+#pragma link off all classes;
+#pragma link off all functions;
+#pragma link C++ nestedclasses;
+
+#pragma link C++ class ReducedData+;
+#pragma link C++ class MEventReduced+;
+
+
+#endif
diff --git a/NPLib/Detectors/Actar/TActarData.cxx b/NPLib/Detectors/Actar/TActarData.cxx
index a9462d0ee..9c569553b 100644
--- a/NPLib/Detectors/Actar/TActarData.cxx
+++ b/NPLib/Detectors/Actar/TActarData.cxx
@@ -6,7 +6,7 @@
*****************************************************************************/
/*****************************************************************************
- * Original Author: Pierre Morfouace contact address: morfouac@nscl.msu.edu *
+ * Original Author: Pierre Morfouace contact address: morfouace@ganil.fr *
* *
* Creation Date : September 2017 *
* Last update : *
@@ -46,12 +46,10 @@ TActarData::~TActarData() {
void TActarData::Clear() {
// Charge
fActar_PadNumber.clear();
- fActar_PadRow.clear();
- fActar_PadColumn.clear();
+ fActar_PadX.clear();
+ fActar_PadY.clear();
+ fActar_PadZ.clear();
fActar_PadCharge.clear();
- fActar_PadTime.clear();
- fActar_dig_Charge.clear();
- fActar_dig_Time.clear();
fSilicon_Energy.clear();
fSilicon_Time.clear();
@@ -61,13 +59,6 @@ void TActarData::Clear() {
fCsI_CrystalNumber.clear();
}
-////////////////////////////////////////////////////////////////////////////////
-TGraph* TActarData::GetChargeAsGraph(){
- TGraph* res = new TGraph (fActar_dig_Charge.size(),&fActar_dig_Time[0],&fActar_dig_Charge[0]);
- res->Sort();
- return res;
-
-}
//////////////////////////////////////////////////////////////////////
@@ -81,8 +72,10 @@ void TActarData::Dump() const {
for (size_t i = 0 ; i < mysize ; i++){
cout << "Pad Number: " << fActar_PadNumber[i]
- << " Charge: " << fActar_PadCharge[i]
- << " Time: " << fActar_PadTime[i];
+ << "Charge: " << fActar_PadCharge[i]
+ << " X: " << fActar_PadX[i]
+ << " Y: " << fActar_PadY[i]
+ << " Z: " << fActar_PadZ[i];
}
}
diff --git a/NPLib/Detectors/Actar/TActarData.h b/NPLib/Detectors/Actar/TActarData.h
index 385512cd1..5674c6aa4 100644
--- a/NPLib/Detectors/Actar/TActarData.h
+++ b/NPLib/Detectors/Actar/TActarData.h
@@ -8,7 +8,7 @@
*****************************************************************************/
/*****************************************************************************
- * Original Author: Pierre Morfouace contact address: morfouac@nscl.msu.edu *
+ * Original Author: Pierre Morfouace contact address: morfouace@ganil.fr *
* *
* Creation Date : September 2017 *
* Last update : *
@@ -37,14 +37,11 @@ class TActarData : public TObject {
private:
// Charge
vector<int> fActar_PadNumber;
- vector<int> fActar_PadTime;
- vector<int> fActar_PadRow;
- vector<int> fActar_PadColumn;
+ vector<int> fActar_PadX;
+ vector<int> fActar_PadY;
+ vector<int> fActar_PadZ;
vector<double> fActar_PadCharge;
- vector<double> fActar_dig_Time;
- vector<double> fActar_dig_Charge;
-
vector<double> fSilicon_Energy;
vector<int> fSilicon_DetectorNumber;
vector<double> fSilicon_Time;
@@ -83,25 +80,19 @@ class TActarData : public TObject {
inline void SetPadNumber(const UShort_t& PadNbr){
fActar_PadNumber.push_back(PadNbr);
};//!
- // Row
- inline void SetRowNumber(const UShort_t& RowNbr){
- fActar_PadRow.push_back(RowNbr);
+ // Row Y
+ inline void SetPadY(const UShort_t& RowNbr){
+ fActar_PadY.push_back(RowNbr);
}
- // Column
- inline void SetColumnNumber(const UShort_t& ColumnNbr){
- fActar_PadColumn.push_back(ColumnNbr);
+ // Column X
+ inline void SetPadX(const UShort_t& ColumnNbr){
+ fActar_PadX.push_back(ColumnNbr);
}
- // Time
- inline void SetTime(const Double_t& Time) {
- fActar_PadTime.push_back(Time);
+ // Time Z
+ inline void SetPadZ(const Double_t& Time) {
+ fActar_PadZ.push_back(Time);
};//!
- // Digitizer
- void AddEnergyPoint(double& E,double& T){
- fActar_dig_Charge.push_back(E);
- fActar_dig_Time.push_back(T);
- }
-
//Silicon
inline void SetSiliconEnergy(const Double_t& Energy){
fSilicon_Energy.push_back(Energy);
@@ -124,23 +115,20 @@ class TActarData : public TObject {
////////////////////// GETTERS ////////////////////////
// Charge
- inline UShort_t GetMultCharge() const
- {return fActar_PadNumber.size();}
- inline UShort_t Get_PadNumber(const unsigned int &i) const
+ inline UShort_t GetPadMult() const
+ {return fActar_PadX.size();}
+ inline UShort_t GetPadNumber(const unsigned int &i) const
{return fActar_PadNumber[i];}//!
- inline UShort_t Get_PadRow(const unsigned int &i) const
- {return fActar_PadRow[i];}//!
- inline UShort_t Get_PadColumn(const unsigned int &i) const
- {return fActar_PadColumn[i];}//!
- inline Double_t Get_Charge(const unsigned int &i) const
+ inline UShort_t GetPadX(const unsigned int &i) const
+ {return fActar_PadX[i];}//!
+ inline UShort_t GetPadY(const unsigned int &i) const
+ {return fActar_PadY[i];}//!
+ inline Double_t GetPadZ(const unsigned int &i) const
+ {return fActar_PadZ[i];}//!
+ inline Double_t GetPadCharge(const unsigned int &i) const
{return fActar_PadCharge[i];}//!
+
- // Time
- inline Double_t Get_Time(const unsigned int &i) const
- {return fActar_PadTime[i];}//!
-
- inline vector<double> Get_dig_Charge() {return fActar_dig_Charge;}
- TGraph* GetChargeAsGraph();
// Silicon
inline Double_t Get_SiliconEnergy(const unsigned int &i) const
diff --git a/NPLib/Detectors/Actar/TActarPhysics.cxx b/NPLib/Detectors/Actar/TActarPhysics.cxx
index 644cbd4d1..4fc143b32 100644
--- a/NPLib/Detectors/Actar/TActarPhysics.cxx
+++ b/NPLib/Detectors/Actar/TActarPhysics.cxx
@@ -6,7 +6,7 @@
*****************************************************************************/
/*****************************************************************************
- * Original Author: Pierre Morfouace contact address: morfouac@nscl.msu.edu *
+ * Original Author: Pierre Morfouace contact address: morfouace@ganil.fr *
* *
* Creation Date : September 2017 *
* Last update : *
@@ -47,305 +47,506 @@ ClassImp(TActarPhysics)
///////////////////////////////////////////////////////////////////////////
TActarPhysics::TActarPhysics()
- : m_EventData(new TActarData),
- m_PreTreatedData(new TActarData),
- m_EventPhysics(this),
- m_Spectra(0),
- m_E_RAW_Threshold(0), // adc channels
- m_E_Threshold(0), // MeV
- m_NumberOfDetectors(0) {
+: m_EventData(new TActarData),
+m_EventReduced(new MEventReduced),
+m_PreTreatedData(new TActarData),
+m_EventPhysics(this),
+
+m_Spectra(0),
+fRecoRansac(1),
+fRecoVisu(0),
+fHitThreshold(20),
+fQ_Threshold(0),
+fT_Threshold(0),
+fNumberOfPadsX(128),
+fNumberOfPadsY(128),
+fPadSizeX(2),
+fPadSizeY(2),
+fDriftVelocity(40),
+fPressure(100),
+fGas("iC4H10"),
+m_NumberOfDetectors(0) {
}
///////////////////////////////////////////////////////////////////////////
/// A usefull method to bundle all operation to add a detector
void TActarPhysics::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++;
+ // 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 TActarPhysics::AddDetector(double R, double Theta, double Phi, string shape){
- // Compute the TVector3 corresponding
- TVector3 Pos(R*sin(Theta)*cos(Phi),R*sin(Theta)*sin(Phi),R*cos(Theta));
- // Call the cartesian method
- AddDetector(Pos,shape);
+ // Compute the TVector3 corresponding
+ TVector3 Pos(R*sin(Theta)*cos(Phi),R*sin(Theta)*sin(Phi),R*cos(Theta));
+ // Call the cartesian method
+ AddDetector(Pos,shape);
}
///////////////////////////////////////////////////////////////////////////
void TActarPhysics::BuildSimplePhysicalEvent() {
- BuildPhysicalEvent();
+ BuildPhysicalEvent();
}
///////////////////////////////////////////////////////////////////////////
void TActarPhysics::BuildPhysicalEvent() {
- // apply thresholds and calibration
-
- PreTreat();
-
- // match Charge and time together
- unsigned int mysizeE = m_PreTreatedData->GetMultCharge();
-
- for (UShort_t e = 0; e < mysizeE ; e++) {
- PadNumber.push_back(m_PreTreatedData->Get_PadNumber(e));
- PadRow.push_back(m_PreTreatedData->Get_PadRow(e));
- PadColumn.push_back(m_PreTreatedData->Get_PadColumn(e));
- PadCharge.push_back(m_PreTreatedData->Get_Charge(e));
- PadTime.push_back(m_PreTreatedData->Get_Time(e));
- }
-
- HoughTransform(PadRow, PadColumn, PadTime);
+
+ PreTreat();
+
+ //unsigned int mysize = m_PreTreatedData->GetPadMult();
+
+ /*for (unsigned int e = 0; e < mysize; e++) {
+ PadX.push_back(m_PreTreatedData->GetPadX(e));
+ PadY.push_back(m_PreTreatedData->GetPadY(e));
+ PadZ.push_back(m_PreTreatedData->GetPadZ(e));
+ PadCharge.push_back(m_PreTreatedData->GetPadCharge(e));
+ }*/
+
+ if(fRecoRansac && PadX.size()>fHitThreshold){
+ m_Ransac->Init(PadX, PadY, PadZ, PadCharge);
+ m_Track = m_Ransac->SimpleRansac();
+ }
+
+ TrackMult = m_Track.size();
}
///////////////////////////////////////////////////////////////////////////
void TActarPhysics::PreTreat() {
- // This method typically applies thresholds and calibrations
- // Might test for disabled channels for more complex detector
-
- // clear pre-treated object
- ClearPreTreatedData();
-
- // instantiate CalibrationManager
- static CalibrationManager* Cal = CalibrationManager::getInstance();
-
- // Charge
- unsigned int mysize = m_EventData->GetMultCharge();
-
- for (UShort_t i = 0; i < mysize ; ++i) {
- if (m_EventData->Get_Charge(i) > m_E_RAW_Threshold) {
- //Double_t Charge = Cal->ApplyCalibration("Actar/CHARGE"+NPL::itoa(m_EventData->Get_PadNumber(i)),m_EventData->Get_Charge(i));
- //Double_t Time= Cal->ApplyCalibration("Actar/TIME"+NPL::itoa(m_EventData->Get_PadNumber(i)),m_EventData->Get_Time(i));
- if (m_EventData->Get_Charge(i) > m_E_Threshold) {
- m_PreTreatedData->SetCharge(m_EventData->Get_Charge(i));
- m_PreTreatedData->SetPadNumber(m_EventData->Get_PadNumber(i));
- m_PreTreatedData->SetRowNumber(m_EventData->Get_PadRow(i));
- m_PreTreatedData->SetColumnNumber(m_EventData->Get_PadColumn(i));
- m_PreTreatedData->SetTime(m_EventData->Get_Time(i));
- }
+ // This method typically applies thresholds and calibrations
+ // Might test for disabled channels for more complex detector
+
+ // clear pre-treated object
+ ClearPreTreatedData();
+
+ CleanPads();
+
+ // instantiate CalibrationManager
+ static CalibrationManager* Cal = CalibrationManager::getInstance();
+
+
+ unsigned int mysize = m_EventReduced->CoboAsad.size();
+
+ for (unsigned int it = 0; it < mysize ; ++it) {
+ int co=m_EventReduced->CoboAsad[it].globalchannelid>>11;
+ int as=(m_EventReduced->CoboAsad[it].globalchannelid - (co<<11))>>9;
+ int ag=(m_EventReduced->CoboAsad[it].globalchannelid - (co<<11)-(as<<9))>>7;
+ int ch=m_EventReduced->CoboAsad[it].globalchannelid - (co<<11)-(as<<9)-(ag<<7);
+ int where=co*NumberOfASAD*NumberOfAGET*NumberOfChannel + as*NumberOfAGET*NumberOfChannel + ag*NumberOfChannel + ch;
+
+ if(co!=31){
+ unsigned int vector_size = m_EventReduced->CoboAsad[it].peakheight.size();
+ for(unsigned int hh=0; hh<vector_size; hh++){
+ if(m_EventReduced->CoboAsad[it].peakheight[hh]>fQ_Threshold && m_EventReduced->CoboAsad[it].peaktime[hh]>fT_Threshold){
+ if(GoodHit(TABLE[4][where],TABLE[5][where])){
+ //if(Hit[TABLE[4][where]][TABLE[5][where]]<2){
+ PadCharge.push_back(m_EventReduced->CoboAsad[it].peakheight[hh]);
+ PadX.push_back(TABLE[4][where]);
+ PadY.push_back(TABLE[5][where]);
+ PadZ.push_back(m_EventReduced->CoboAsad[it].peaktime[hh]);
+
+ //}
+ }
+ }
+ }
+ }
+ else if(co==31){
+ unsigned int vector_size = m_EventReduced->CoboAsad[it].peakheight.size();
+ for(unsigned int hit=0;hit<vector_size;hit++){
+
+ Si_Number.push_back(m_EventReduced->CoboAsad[it].peaktime[hit]);
+ Si_E.push_back(m_EventReduced->CoboAsad[it].peakheight[hit]);
+
+ int vxi_parameter = m_EventReduced->CoboAsad[it].peaktime[hit];
+ if(Si_map[vxi_parameter]<21 && Si_map[vxi_parameter]>0){
+ m_PreTreatedData->SetSiliconEnergy(m_EventReduced->CoboAsad[it].peakheight[hit]);
+ m_PreTreatedData->SetSiliconDetectorNumber(Si_map[vxi_parameter]);
+ }
+ }
+ }
}
- }
+}
+///////////////////////////////////////////////////////////////////////////
+bool TActarPhysics::GoodHit(int iX, int iY)
+{
+ bool bHit = true;
+ if(Hit[iX][iY]<2){
+ if(Hit[iX+1][iY]>0 || Hit[iX+1][iY-1]>0 || Hit[iX+1][iY+1]>0){
+ if(Hit[iX+2][iY]>0 || Hit[iX+2][iY-1]>0 || Hit[iX+2][iY+1]>0){
+ bHit = true;
+ }
+ }
+ }
+
+ return bHit;
}
+
///////////////////////////////////////////////////////////////////////////
-void TActarPhysics::HoughTransform(vector<int> v1, vector<int> v2, vector<int> v3){
- for(unsigned int i=0; i<v1.size(); i++){
- for(int itheta=0; itheta<180; itheta++){
- HoughAngle.push_back(itheta);
- double tmp= v1[i]*cos(itheta*NPUNITS::deg)+v2[i]*sin(itheta*NPUNITS::deg);
- HoughRadius.push_back(tmp);
+void TActarPhysics::CleanPads()
+{
+ unsigned int mysize = m_EventReduced->CoboAsad.size();
+
+ for(unsigned int it=0; it<mysize; it++){
+ int co=m_EventReduced->CoboAsad[it].globalchannelid>>11;
+ int as=(m_EventReduced->CoboAsad[it].globalchannelid - (co<<11))>>9;
+ int ag=(m_EventReduced->CoboAsad[it].globalchannelid - (co<<11)-(as<<9))>>7;
+ int ch=m_EventReduced->CoboAsad[it].globalchannelid - (co<<11)-(as<<9)-(ag<<7);
+ int where=co*NumberOfASAD*NumberOfAGET*NumberOfChannel + as*NumberOfAGET*NumberOfChannel + ag*NumberOfChannel + ch;
+
+
+ if(co!=31){
+ unsigned int vector_size = m_EventReduced->CoboAsad[it].peakheight.size();
+ for(unsigned int hh=0; hh < vector_size; hh++){
+ if(m_EventReduced->CoboAsad[it].peakheight[hh]>fQ_Threshold && m_EventReduced->CoboAsad[it].peaktime[hh]>fT_Threshold){
+ Hit[TABLE[4][where]][TABLE[5][where]] += 1;
+ /*
+ m_PreTreatedData->SetCharge(m_EventReduced->CoboAsad[it].peakheight[hh]);
+ m_PreTreatedData->SetPadX(TABLE[4][where]);
+ m_PreTreatedData->SetPadY(TABLE[5][where]);
+ m_PreTreatedData->SetPadZ(m_EventReduced->CoboAsad[it].peaktime[hh]);*/
+ }
+ }
+ }
}
- }
}
///////////////////////////////////////////////////////////////////////////
void TActarPhysics::ReadAnalysisConfig() {
- bool ReadingStatus = false;
-
- // path to file
- string FileName = "./configs/ConfigActar.dat";
-
- // open analysis config file
- ifstream AnalysisConfigFile;
- AnalysisConfigFile.open(FileName.c_str());
-
- if (!AnalysisConfigFile.is_open()) {
- cout << " No ConfigActar.dat found: Default parameter loaded for Analayis " << FileName << endl;
- return;
- }
- cout << " Loading user parameter for Analysis from ConfigActar.dat " << endl;
-
- // Save it in a TAsciiFile
- TAsciiFile* asciiConfig = RootOutput::getInstance()->GetAsciiFileAnalysisConfig();
- asciiConfig->AppendLine("%%% ConfigActar.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 = "ConfigActar";
- 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;
- }
+ // VXI ACTION FILE //
+ string VXI_FileName = "./configs/ACTION_Si_config.dat";
+ ifstream VXIConfigFile;
+ VXIConfigFile.open(VXI_FileName.c_str());
+ if(!VXIConfigFile.is_open()){
+ cout << "No VXI ACTION FILE ./configs/ACTION_Si_config.dat found!" << endl;
+ return;
+ }
+ else{
+ cout << "/// Using VXI ACTION FILE: " << VXI_FileName << " ///" << endl;
+ string token;
+ int vxi_param, si_nbr;
+ for(int i=0; i<20; i++){
+ VXIConfigFile >> token >> vxi_param >> si_nbr;
+ Si_map[vxi_param] = si_nbr+1;
+ }
+ }
+ VXIConfigFile.close();
+
+ // Lookup table //
+ bool ReadingLookupTable = false;
+ string LT_FileName = "./configs/LT.dat";
+ ifstream LTConfigFile;
+ LTConfigFile.open(LT_FileName.c_str());
+ if(!LTConfigFile.is_open()){
+ cout << "No Lookup Table in ./configs/LT.dat found!" << endl;
+ return;
+ }
+ else{
+ cout << "/// Using LookupTable from: " << LT_FileName << " ///" << endl;
+ for(int i=0;i<NumberOfCobo*NumberOfASAD*NumberOfAGET*NumberOfChannel;i++){
+ LTConfigFile >> TABLE[0][i] >> TABLE[1][i] >> TABLE[2][i] >> TABLE[3][i] >> TABLE[4][i] >> TABLE[5][i];
+ }
+ ReadingLookupTable = true;
+ }
+ LTConfigFile.close();
+
+
+ bool ReadingStatus = false;
+
+ // ACTAR config file //
+ string FileName = "./configs/ConfigActar.dat";
+
+ // open analysis config file
+ ifstream AnalysisConfigFile;
+ AnalysisConfigFile.open(FileName.c_str());
+
+ if (!AnalysisConfigFile.is_open()) {
+ cout << " No ConfigActar.dat found: Default parameter loaded for Analayis " << FileName << endl;
+ return;
+ }
+ cout << "/// Loading user parameter for Analysis from ConfigActar.dat ///" << endl;
+
+ // Save it in a TAsciiFile
+ TAsciiFile* asciiConfig = RootOutput::getInstance()->GetAsciiFileAnalysisConfig();
+ asciiConfig->AppendLine("%%% ConfigActar.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 = "ConfigActar";
+ if (LineBuffer.compare(0, name.length(), name) == 0){
+ ReadingStatus = true;
+ cout << "**** ConfigActar found" << endl;
+ }
+
+ // 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.compare(0,11,"RecoRansac=") == 0) {
+ AnalysisConfigFile >> DataBuffer;
+ fRecoRansac = atoi(DataBuffer.c_str());
+ cout << "/// Reco using Ransac= " << " " << fRecoRansac << " ///" << endl;
+ }
+
+ else if (whatToDo.compare(0,9,"RecoVisu=") == 0) {
+ AnalysisConfigFile >> DataBuffer;
+ fRecoVisu = atoi(DataBuffer.c_str());
+ cout << "/// Visualisation= " << " " << fRecoVisu << " ///" << endl;
+ }
+
+ else if (whatToDo.compare(0,14,"HIT_THRESHOLD=") == 0) {
+ AnalysisConfigFile >> DataBuffer;
+ fHitThreshold = atoi(DataBuffer.c_str());
+ cout << "/// Hit Threshold= " << " " << fHitThreshold << " ///" << endl;
+ }
+
+ else if (whatToDo.compare(0,12,"Q_THRESHOLD=") == 0) {
+ AnalysisConfigFile >> DataBuffer;
+ fQ_Threshold = atof(DataBuffer.c_str());
+ cout << "/// Q Threshold= " << " " << fQ_Threshold << " ///" << endl;
+ }
+
+ else if (whatToDo.compare(0,12,"T_THRESHOLD=") == 0) {
+ AnalysisConfigFile >> DataBuffer;
+ fT_Threshold = atof(DataBuffer.c_str());
+ cout << "/// T Threshold= " << " " << fT_Threshold << " ///" << endl;
+ }
+ else if(whatToDo.compare(0,14,"NumberOfPadsX=")==0){
+ AnalysisConfigFile >> DataBuffer;
+ fNumberOfPadsX = atoi(DataBuffer.c_str());
+ //check_padsX=true;
+ cout << "/// Number Of Pads X= " << fNumberOfPadsX << " ///" << endl;
+ }
+
+ else if(whatToDo.compare(0,14,"NumberOfPadsY=")==0){
+ AnalysisConfigFile >> DataBuffer;
+ fNumberOfPadsY = atoi(DataBuffer.c_str());
+ //check_padsY=true;
+ cout << "/// Number Of Pads Y= " << fNumberOfPadsY << " ///" << endl;
+ }
+
+ else if(whatToDo.compare(0,9,"PadSizeX=")==0){
+ AnalysisConfigFile >> DataBuffer;
+ fPadSizeX = atof(DataBuffer.c_str());
+ //check_sizeX=true;
+ cout << "/// Pad Size X= " << fPadSizeX << " ///" << endl;
+ }
+
+ else if(whatToDo.compare(0,9,"PadSizeY=")==0){
+ AnalysisConfigFile >> DataBuffer;
+ fPadSizeY = atof(DataBuffer.c_str());
+ //check_sizeY=true;
+ cout << "/// Pad Size Y= " << fPadSizeY << " ///" << endl;
+ }
+
+ else if(whatToDo.compare(0,9,"Pressure=")==0){
+ AnalysisConfigFile >> DataBuffer;
+ fPressure = atof(DataBuffer.c_str());
+ //check_pressure=true;
+ cout << "/// Pressure= " << fPressure << " ///" << endl;
+ }
+
+ else if(whatToDo.compare(0,14,"DriftVelocity=")==0){
+ AnalysisConfigFile >> DataBuffer;
+ fDriftVelocity = atof(DataBuffer.c_str());
+ //check_driftvelocity=true;
+ cout << "/// Drift Velocity= " << fDriftVelocity << " ///" << endl;
+ }
+
+ else if(whatToDo.compare(0,4,"Gas=")==0){
+ AnalysisConfigFile >> DataBuffer;
+ fGas = DataBuffer.c_str();
+ //check_gas=true;
+ cout << "/// Gas Type= " << fGas << " ///" << endl;
+ }
+
+ else {
+ ReadingStatus = false;
+ }
+ }
+ }
+
+ if(fRecoRansac){
+ m_Ransac = new NPL::Ransac(fNumberOfPadsX,fNumberOfPadsY,fRecoVisu);
}
- }
}
-
+///////////////////////////////////////////////////////////////////////////
+void TActarPhysics::SetRansacParameter(string filename){
+ if(fRecoRansac){
+ m_Ransac->ReadParameterValue(filename);
+ }
+}
///////////////////////////////////////////////////////////////////////////
void TActarPhysics::Clear() {
- PadNumber.clear();
- PadRow.clear();
- PadColumn.clear();
- PadCharge.clear();
- PadTime.clear();
-
- HoughRadius.clear();
- HoughAngle.clear();
+ PadX.clear();
+ PadY.clear();
+ PadZ.clear();
+ PadCharge.clear();
+ m_Track.clear();
+ Si_Number.clear();
+ Si_E.clear();
+
+ for(int i=0; i<fNumberOfPadsX; i++){
+ for(int j=0; j<fNumberOfPadsY; j++){
+ Hit[i][j] = 0;
+ }
+ }
}
///////////////////////////////////////////////////////////////////////////
void TActarPhysics::ReadConfiguration(NPL::InputParser parser) {
- vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Actar");
- if(NPOptionManager::getInstance()->GetVerboseLevel())
- cout << "//// " << blocks.size() << " detectors found " << endl;
-
- vector<string> cart = {"POS","Shape"};
- vector<string> sphe = {"R","Theta","Phi","Shape"};
-
- for(unsigned int i = 0 ; i < blocks.size() ; i++){
- if(blocks[i]->HasTokenList(cart)){
- if(NPOptionManager::getInstance()->GetVerboseLevel())
- cout << endl << "//// Actar " << i+1 << endl;
-
- TVector3 Pos = blocks[i]->GetTVector3("POS","mm");
- string Shape = blocks[i]->GetString("Shape");
- AddDetector(Pos,Shape);
- }
- else if(blocks[i]->HasTokenList(sphe)){
- if(NPOptionManager::getInstance()->GetVerboseLevel())
- cout << endl << "//// Actar " << i+1 << endl;
- double R = blocks[i]->GetDouble("R","mm");
- double Theta = blocks[i]->GetDouble("Theta","deg");
- double Phi = blocks[i]->GetDouble("Phi","deg");
- string Shape = blocks[i]->GetString("Shape");
- AddDetector(R,Theta,Phi,Shape);
- }
- else{
- cout << "ERROR: check your input file formatting " << endl;
- exit(1);
+ vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Actar");
+ if(NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << "//// " << blocks.size() << " detectors found " << endl;
+
+ vector<string> cart = {"POS","Shape"};
+ vector<string> sphe = {"R","Theta","Phi","Shape"};
+
+ for(unsigned int i = 0 ; i < blocks.size() ; i++){
+ if(blocks[i]->HasTokenList(cart)){
+ if(NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << endl << "//// Actar " << i+1 << endl;
+
+ TVector3 Pos = blocks[i]->GetTVector3("POS","mm");
+ string Shape = blocks[i]->GetString("Shape");
+ AddDetector(Pos,Shape);
+ }
+ else if(blocks[i]->HasTokenList(sphe)){
+ if(NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << endl << "//// Actar " << i+1 << endl;
+ double R = blocks[i]->GetDouble("R","mm");
+ double Theta = blocks[i]->GetDouble("Theta","deg");
+ double Phi = blocks[i]->GetDouble("Phi","deg");
+ string Shape = blocks[i]->GetString("Shape");
+ AddDetector(R,Theta,Phi,Shape);
+ }
+ else{
+ cout << "ERROR: check your input file formatting " << endl;
+ exit(1);
+ }
}
- }
}
///////////////////////////////////////////////////////////////////////////
void TActarPhysics::InitSpectra() {
- m_Spectra = new TActarSpectra(m_NumberOfDetectors);
+ m_Spectra = new TActarSpectra(m_NumberOfDetectors);
}
///////////////////////////////////////////////////////////////////////////
void TActarPhysics::FillSpectra() {
- m_Spectra -> FillRawSpectra(m_EventData);
- m_Spectra -> FillPreTreatedSpectra(m_PreTreatedData);
- m_Spectra -> FillPhysicsSpectra(m_EventPhysics);
+ m_Spectra -> FillRawSpectra(m_EventData);
+ m_Spectra -> FillPreTreatedSpectra(m_PreTreatedData);
+ m_Spectra -> FillPhysicsSpectra(m_EventPhysics);
}
///////////////////////////////////////////////////////////////////////////
void TActarPhysics::CheckSpectra() {
- m_Spectra->CheckSpectra();
+ m_Spectra->CheckSpectra();
}
///////////////////////////////////////////////////////////////////////////
void TActarPhysics::ClearSpectra() {
- // To be done
+ // To be done
}
///////////////////////////////////////////////////////////////////////////
map< string , TH1*> TActarPhysics::GetSpectra() {
- if(m_Spectra)
- return m_Spectra->GetMapHisto();
- else{
- map< string , TH1*> empty;
- return empty;
- }
+ if(m_Spectra)
+ return m_Spectra->GetMapHisto();
+ else{
+ map< string , TH1*> empty;
+ return empty;
+ }
}
////////////////////////////////////////////////////////////////////////////////
vector<TCanvas*> TActarPhysics::GetCanvas() {
- //if(m_Spectra)
+ //if(m_Spectra)
//return m_Spectra->GetCanvas();
- //else{
- //vector<TCanvas*> empty;
- //return empty;
- //}
- return vector<TCanvas*>();
+ //else{
+ vector<TCanvas*> empty;
+ return empty;
+ //}
}
///////////////////////////////////////////////////////////////////////////
void TActarPhysics::WriteSpectra() {
- m_Spectra->WriteSpectra();
+ m_Spectra->WriteSpectra();
}
///////////////////////////////////////////////////////////////////////////
void TActarPhysics::AddParameterToCalibrationManager() {
- CalibrationManager* Cal = CalibrationManager::getInstance();
- for (int i = 0; i < m_NumberOfDetectors; ++i) {
- Cal->AddParameter("Actar", "D"+ NPL::itoa(i+1)+"_CHARGE","Actar_D"+ NPL::itoa(i+1)+"_CHARGE");
- Cal->AddParameter("Actar", "D"+ NPL::itoa(i+1)+"_TIME","Actar_D"+ NPL::itoa(i+1)+"_TIME");
- }
+ CalibrationManager* Cal = CalibrationManager::getInstance();
+ for (int i = 0; i < m_NumberOfDetectors; ++i) {
+ Cal->AddParameter("Actar", "D"+ NPL::itoa(i+1)+"_CHARGE","Actar_D"+ NPL::itoa(i+1)+"_CHARGE");
+ Cal->AddParameter("Actar", "D"+ NPL::itoa(i+1)+"_TIME","Actar_D"+ NPL::itoa(i+1)+"_TIME");
+ }
}
+///////////////////////////////////////////////////////////////////////////
+/*void TActarPhysics::InitializeRootInputRaw() {
+ TChain* inputChain = RootInput::getInstance()->GetChain();
+ inputChain->SetBranchStatus("Actar", true );
+ inputChain->SetBranchAddress("Actar", &m_EventData );
+}*/
+
///////////////////////////////////////////////////////////////////////////
void TActarPhysics::InitializeRootInputRaw() {
- TChain* inputChain = RootInput::getInstance()->GetChain();
- inputChain->SetBranchStatus("Actar", true );
- inputChain->SetBranchAddress("Actar", &m_EventData );
+ TChain* inputChain = RootInput::getInstance()->GetChain();
+ inputChain->SetBranchStatus("data", true );
+ inputChain->SetBranchStatus("ACTAR_TTree", true );
+ inputChain->SetBranchAddress("data", &m_EventReduced );
}
///////////////////////////////////////////////////////////////////////////
void TActarPhysics::InitializeRootInputPhysics() {
- TChain* inputChain = RootInput::getInstance()->GetChain();
- inputChain->SetBranchAddress("Actar", &m_EventPhysics);
+ TChain* inputChain = RootInput::getInstance()->GetChain();
+ inputChain->SetBranchAddress("Actar", &m_EventPhysics);
}
///////////////////////////////////////////////////////////////////////////
void TActarPhysics::InitializeRootOutput() {
- TTree* outputTree = RootOutput::getInstance()->GetTree();
- outputTree->Branch("Actar", "TActarPhysics", &m_EventPhysics);
+ TTree* outputTree = RootOutput::getInstance()->GetTree();
+ outputTree->Branch("Actar", "TActarPhysics", &m_EventPhysics);
}
@@ -354,7 +555,7 @@ void TActarPhysics::InitializeRootOutput() {
// Construct Method to be pass to the DetectorFactory //
////////////////////////////////////////////////////////////////////////////////
NPL::VDetector* TActarPhysics::Construct() {
- return (NPL::VDetector*) new TActarPhysics();
+ return (NPL::VDetector*) new TActarPhysics();
}
@@ -363,13 +564,13 @@ NPL::VDetector* TActarPhysics::Construct() {
// Registering the construct method to the factory //
////////////////////////////////////////////////////////////////////////////////
extern "C"{
-class proxy_Actar{
- public:
- proxy_Actar(){
- NPL::DetectorFactory::getInstance()->AddToken("Actar","Actar");
- NPL::DetectorFactory::getInstance()->AddDetector("Actar",TActarPhysics::Construct);
- }
-};
-
-proxy_Actar p_Actar;
+ class proxy_Actar{
+ public:
+ proxy_Actar(){
+ NPL::DetectorFactory::getInstance()->AddToken("Actar","Actar");
+ NPL::DetectorFactory::getInstance()->AddDetector("Actar",TActarPhysics::Construct);
+ }
+ };
+
+ proxy_Actar p_Actar;
}
diff --git a/NPLib/Detectors/Actar/TActarPhysics.h b/NPLib/Detectors/Actar/TActarPhysics.h
index 6240921b0..ab740fa95 100644
--- a/NPLib/Detectors/Actar/TActarPhysics.h
+++ b/NPLib/Detectors/Actar/TActarPhysics.h
@@ -8,7 +8,7 @@
*****************************************************************************/
/*****************************************************************************
- * Original Author: Pierre Morfouace contact address: morfouac@nscl.msu.edu *
+ * Original Author: Pierre Morfouace contact address: morfouace@ganil.fr *
* *
* Creation Date : September 2017 *
* Last update : *
@@ -33,157 +33,203 @@ using namespace std;
#include "TH1.h"
#include "TCanvas.h"
#include "TVector3.h"
+
// NPTool headers
+#include "MEventReduced.h"
#include "TActarData.h"
#include "TActarSpectra.h"
#include "NPCalibrationManager.h"
#include "NPVDetector.h"
#include "NPInputParser.h"
-// forward declaration
-class TActarSpectra;
+#include "NPTrack.h"
+#include "NPRansac.h"
+#define NumberOfCobo 16
+#define NumberOfASAD 4
+#define NumberOfAGET 4
+#define NumberOfChannel 68
+class TActarSpectra;
class TActarPhysics : public TObject, public NPL::VDetector {
- //////////////////////////////////////////////////////////////
- // constructor and destructor
- public:
+ //////////////////////////////////////////////////////////////
+ // constructor and destructor
+public:
TActarPhysics();
~TActarPhysics() {};
-
-
- //////////////////////////////////////////////////////////////
- // Inherited from TObject and overriden to avoid warnings
- public:
+
+
+ //////////////////////////////////////////////////////////////
+ // Inherited from TObject and overriden to avoid warnings
+public:
void Clear();
void Clear(const Option_t*) {};
-
-
- //////////////////////////////////////////////////////////////
- // data obtained after BuildPhysicalEvent() and stored in
- // output ROOT file
- public:
- vector<int> PadNumber;
- vector<int> PadRow;
- vector<int> PadColumn;
- vector<int> PadTime;
+
+
+ //////////////////////////////////////////////////////////////
+ // data obtained after BuildPhysicalEvent() and stored in
+ // output ROOT file
+public:
+ vector<int> PadX;
+ vector<int> PadY;
+ vector<double> PadZ;
vector<double> PadCharge;
-
- vector<double> HoughRadius;
- vector<double> HoughAngle;
-
- /// A usefull method to bundle all operation to add a detector
- void AddDetector(TVector3 POS, string shape);
- void AddDetector(double R, double Theta, double Phi, string shape);
-
- //////////////////////////////////////////////////////////////
- // methods inherited from the VDetector ABC class
- public:
+ vector<double> Si_E;
+ vector<int> Si_Number;
+ int TrackMult;
+
+ /// A usefull method to bundle all operation to add a detector
+ void AddDetector(TVector3 POS, string shape);
+ void AddDetector(double R, double Theta, double Phi, string shape);
+
+public:
+ int GetTrackMult() {return m_Track.size();}
+ vector<NPL::Track> GetTracks() {return m_Track;}
+
+
+ //////////////////////////////////////////////////////////////
+ // methods inherited from the VDetector ABC class
+public:
// read stream from ConfigFile to pick-up detector parameters
void ReadConfiguration(NPL::InputParser);
-
+
// add parameters to the CalibrationManger
void AddParameterToCalibrationManager();
-
+
// method called event by event, aiming at extracting the
// physical information from detector
void BuildPhysicalEvent();
-
+
// same as BuildPhysicalEvent() method but with a simpler
// treatment
void BuildSimplePhysicalEvent();
-
+
// same as above but for online analysis
void BuildOnlinePhysicalEvent() {BuildPhysicalEvent();};
-
+
// activate raw data object and branches from input TChain
// in this method mother branches (Detector) AND daughter leaves
// (fDetector_parameter) have to be activated
void InitializeRootInputRaw();
-
+
// activate physics data object and branches from input TChain
// in this method mother branches (Detector) AND daughter leaves
// (fDetector_parameter) have to be activated
void InitializeRootInputPhysics();
-
+
// create branches of output ROOT file
void InitializeRootOutput();
-
+
// clear the raw and physical data objects event by event
void ClearEventPhysics() {Clear();}
void ClearEventData() {m_EventData->Clear();}
-
+
// methods related to the TActarSpectra class
// instantiate the TActarSpectra class and
// declare list of histograms
void InitSpectra();
-
+
// fill the spectra
void FillSpectra();
-
+
// used for Online mainly, sanity check for histograms and
// change their color if issues are found, for example
void CheckSpectra();
-
+
// used for Online only, clear all the spectra
void ClearSpectra();
-
+
// write spectra to ROOT output file
void WriteSpectra();
-
-
- //////////////////////////////////////////////////////////////
- // specific methods to Actar array
- public:
+
+
+ //////////////////////////////////////////////////////////////
+ // specific methods to Actar array
+public:
// remove bad channels, calibrate the data and apply thresholds
void PreTreat();
-
+
// clear the pre-treated object
void ClearPreTreatedData() {m_PreTreatedData->Clear();}
-
+
// read the user configuration file. If no file is found, load standard one
void ReadAnalysisConfig();
-
+
+ void CleanPads();
+
+ bool GoodHit(int iX, int iY);
+
// give and external TActarData object to TActarPhysics.
// needed for online analysis for example
void SetRawDataPointer(TActarData* rawDataPointer) {m_EventData = rawDataPointer;}
-
- void HoughTransform(vector<int> v1, vector<int> v2, vector<int> v3);
-
- // objects are not written in the TTree
- private:
+
+ // objects are not written in the TTree
+private:
TActarData* m_EventData; //!
+ MEventReduced* m_EventReduced; //!
TActarData* m_PreTreatedData; //!
TActarPhysics* m_EventPhysics; //!
-
- // getters for raw and pre-treated data object
- public:
+ NPL::Ransac* m_Ransac; //!
+ vector<NPL::Track> m_Track; //!
+
+ // getters for raw and pre-treated data object
+public:
TActarData* GetRawData() const {return m_EventData;}
TActarData* GetPreTreatedData() const {return m_PreTreatedData;}
-
- // parameters used in the analysis
- private:
+
+ double GetDriftVelocity() {return fDriftVelocity;}
+ double GetPadSizeX() {return fPadSizeX;}
+ double GetPadSizeY() {return fPadSizeY;}
+ int GetNumberOfPadsX() {return fNumberOfPadsX;}
+ int GetNumberOfPadsY() {return fNumberOfPadsY;}
+ double GetPRessure() {return fPressure;}
+ string GetGasName() {return fGas;}
+
+ // parameters used in the analysis
+private:
// thresholds
- double m_E_RAW_Threshold; //!
- double m_E_Threshold; //!
-
- // number of detectors
- private:
+ int fHitThreshold; //!
+ int fQ_Threshold; //!
+ int fT_Threshold; //!
+ int fNumberOfPadsX; //!
+ int fNumberOfPadsY; //!
+ double fPadSizeX; //!
+ double fPadSizeY; //!
+ double fDriftVelocity; //!
+ double fPressure; //!
+ string fGas; //!
+ bool fRecoRansac; //!
+ bool fRecoVisu; //!
+ map<int, int> Si_map; //!
+
+ int TABLE[6][NumberOfCobo*NumberOfASAD*NumberOfAGET*NumberOfChannel]; //!
+ int Hit[128][128]; //!
+
+
+ // number of detectors
+private:
int m_NumberOfDetectors; //!
-
- // spectra class
- private:
- TActarSpectra* m_Spectra; // !
-
- // spectra getter
- public:
+
+ // spectra class
+private:
+ TActarSpectra* m_Spectra; //!
+
+ //spme getters and setters
+public:
+ void SetRansacParameter(string filename);
+ NPL::Ransac* GetRansacObject() {return m_Ransac;}
+ bool GetRansacStatus() {return fRecoRansac;}
+
+
+ // spectra getter
+public:
map<string, TH1*> GetSpectra();
vector<TCanvas*> GetCanvas();
-
- // Static constructor to be passed to the Detector Factory
- public:
+
+ // Static constructor to be passed to the Detector Factory
+public:
static NPL::VDetector* Construct();
-
+
ClassDef(TActarPhysics,1) // ActarPhysics structure
};
#endif
diff --git a/NPLib/Detectors/Actar/TActarSpectra.cxx b/NPLib/Detectors/Actar/TActarSpectra.cxx
index f37527d1a..31e3cd0b0 100644
--- a/NPLib/Detectors/Actar/TActarSpectra.cxx
+++ b/NPLib/Detectors/Actar/TActarSpectra.cxx
@@ -112,21 +112,21 @@ void TActarSpectra::FillRawSpectra(TActarData* RawData) {
static string family;
// Charge
- unsigned int sizeE = RawData->GetMultCharge();
+ unsigned int sizeE = RawData->GetPadMult();
for (unsigned int i = 0; i < sizeE; i++) {
- name = "Actar"+NPL::itoa(RawData->Get_PadNumber(i))+"_CHARGE_RAW";
+ name = "Actar"+NPL::itoa(RawData->GetPadNumber(i))+"_CHARGE_RAW";
family = "Actar/RAW";
- //GetHisto(family,name) -> Fill(RawData->Get_Charge(i));
+ //GetHisto(family,name) -> Fill(RawData->GetPadCharge(i));
}
// Time
- unsigned int sizeT = RawData->GetMultCharge();
+ unsigned int sizeT = RawData->GetPadMult();
for (unsigned int i = 0; i < sizeT; i++) {
- name = "Actar"+NPL::itoa(RawData->Get_PadNumber(i))+"_TIME_RAW";
+ name = "Actar"+NPL::itoa(RawData->GetPadNumber(i))+"_TIME_RAW";
family = "Actar/RAW";
- //GetHisto(family,name) -> Fill(RawData->Get_Time(i));
+ //GetHisto(family,name) -> Fill(RawData->GetPadZ(i));
}
}
@@ -138,21 +138,21 @@ void TActarSpectra::FillPreTreatedSpectra(TActarData* PreTreatedData) {
static string family;
// Energy
- unsigned int sizeE = PreTreatedData->GetMultCharge();
+ unsigned int sizeE = PreTreatedData->GetPadMult();
for (unsigned int i = 0; i < sizeE; i++) {
- name = "Actar"+NPL::itoa(PreTreatedData->Get_PadNumber(i))+"_ENERGY_CAL";
+ name = "Actar"+NPL::itoa(PreTreatedData->GetPadNumber(i))+"_ENERGY_CAL";
family = "Actar/CAL";
- //GetHisto(family,name) -> Fill(PreTreatedData->Get_Charge(i));
+ //GetHisto(family,name) -> Fill(PreTreatedData->GetPadCharge(i));
}
// Time
- unsigned int sizeT = PreTreatedData->GetMultCharge();
+ unsigned int sizeT = PreTreatedData->GetPadMult();
for (unsigned int i = 0; i < sizeT; i++) {
- name = "Actar"+NPL::itoa(PreTreatedData->Get_PadNumber(i))+"_TIME_CAL";
+ name = "Actar"+NPL::itoa(PreTreatedData->GetPadNumber(i))+"_TIME_CAL";
family = "Actar/CAL";
- //GetHisto(family,name) -> Fill(PreTreatedData->Get_Time(i));
+ //GetHisto(family,name) -> Fill(PreTreatedData->GetPadZ(i));
}
}
diff --git a/NPLib/TrackReconstruction/CMakeLists.txt b/NPLib/TrackReconstruction/CMakeLists.txt
new file mode 100644
index 000000000..4d17e167c
--- /dev/null
+++ b/NPLib/TrackReconstruction/CMakeLists.txt
@@ -0,0 +1,6 @@
+add_custom_command(OUTPUT NPRansacDict.cxx COMMAND ../scripts/build_dict.sh NPRansac.h NPRansacDict.cxx NPRansac.rootmap libNPTrackReconstruction.so NPTrackReconstructionLinkDef.h DEPENDS NPRansac.h)
+
+add_library(NPTrackReconstruction SHARED NPRansac.cxx NPTrack.cxx NPRansacDict.cxx)
+target_link_libraries(NPTrackReconstruction ${ROOT_LIBRARIES} NPCore)
+
+install(FILES NPRansac.h NPTrack.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
diff --git a/NPLib/TrackReconstruction/NPRansac.cxx b/NPLib/TrackReconstruction/NPRansac.cxx
new file mode 100644
index 000000000..ec0d36bfb
--- /dev/null
+++ b/NPLib/TrackReconstruction/NPRansac.cxx
@@ -0,0 +1,774 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2016 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * *
+ * Original Author : Pierre MORFOUACE contact address: morfouace@ganil.fr *
+ * *
+ * Creation Date : April 2018 *
+ * Last update : April 2018 *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class deal with finding all the track event by event *
+ *****************************************************************************/
+
+#include "NPRansac.h"
+using namespace NPL;
+
+using namespace std;
+
+ClassImp(Ransac)
+
+/////////////////////////////////////////////////
+Ransac::Ransac()
+{
+ fRANSACMaxIteration = 1000;
+ fRANSACThreshold = 30;//100;
+ fRANSACPointThreshold = 0.04;//0.07;
+ fRANSACChargeThreshold = 0.04;//0.07;
+ fRANSACDistance = 11;//7;
+ fMAXBarycenterDistance = 10;
+ fAngleMax = 2;//degree
+ fMergeWithBarycenter = true;
+ fMergeWithScalarProduct = false;
+
+ fNumberOfPadsX = 128;
+ fNumberOfPadsY = 128;
+}
+
+/////////////////////////////////////////////////
+Ransac::Ransac(int NumberOfPadsX, int NumberOfPadsY, bool Visu)
+{
+ fRANSACMaxIteration = 1000;
+ fRANSACThreshold = 60;//100;
+ fRANSACPointThreshold = 0.04;//0.07;
+ fRANSACChargeThreshold = 0.04;//0.07;
+ fRANSACDistance = 11;//7;
+
+ fNumberOfPadsX = NumberOfPadsX;
+ fNumberOfPadsY = NumberOfPadsY;
+ fVisu = Visu;
+ //fNumberOfTracksMax = 10;
+
+ if(fVisu==1){
+ m_ServerSocket = new TServerSocket(9092,true,100);
+ //pl = new TGraph2D();
+ if(!m_ServerSocket->IsValid()){
+ cout << endl;
+ cout << "/// ServerSocket invalid! ///" << endl;
+ exit(1);
+ }
+ cout << endl;
+ cout << "/// ServerSocket valid ///" << endl;
+ m_ServerSocket->SetCompressionSettings(1);
+
+ // Add server socket to monitor so we are notified when a client needs to be
+ // accepted
+ m_Monitor = new TMonitor();
+ m_Monitor->Add(m_ServerSocket,TMonitor::kRead|TMonitor::kWrite);
+ // Create a list to contain all client connections
+ m_Sockets = new TList;
+
+ // Create a list to contain all client connections
+ m_Histo = new TList;
+
+ }
+}
+
+
+//////////////////////////////////////////////////
+Ransac::~Ransac()
+{
+ for(unsigned int i=0; i<vline.size(); i++){
+ delete vline[i];
+ }
+}
+
+//////////////////////////////////////////////////
+void Ransac::ReadParameterValue(string filename)
+{
+ bool ReadingStatus = false;
+
+ ifstream RansacParamFile;
+ RansacParamFile.open(filename.c_str());
+
+ if(!RansacParamFile.is_open()){
+ cout << "No Paramter File found for Ransac parameters" << endl;
+ cout << "Using the default parameter:" << endl;
+ cout << "RANSAC MaxIteration= " << fRANSACMaxIteration << endl;
+ cout << "RANSAC Threshold=" << fRANSACThreshold << endl;
+ cout << "RANSAC ChargeThreshold= " << fRANSACChargeThreshold << endl;
+ cout << "RANSAC Distance= " << fRANSACDistance << endl;
+ cout << "RANSAC PointThreshold= " << fRANSACPointThreshold << endl;
+ cout << "MAX Barycenter Distance= " << fMAXBarycenterDistance << endl;
+ cout << "Angle Max to merge tracks= " << fAngleMax << endl;
+ cout << "Merging Tracks using the barycenter= " << fMergeWithBarycenter << endl;
+ cout << "Merging Tracks using the scalar product= " << fMergeWithScalarProduct << endl;
+ }
+ else{
+ string LineBuffer, whatToDo, DataBuffer;
+ while(!RansacParamFile.eof()){
+ getline(RansacParamFile,LineBuffer);
+ string name = "ConfigRansac";
+ if(LineBuffer.compare(0, name.length(), name) == 0){
+ cout << endl;
+ cout << "**** ConfigRansac found" << endl;
+ ReadingStatus = true;
+ }
+ while(ReadingStatus){
+ whatToDo="";
+ RansacParamFile >> whatToDo;
+ // Search for comment symbol (%)
+ if (whatToDo.compare(0, 1, "%") == 0) {
+ RansacParamFile.ignore(numeric_limits<streamsize>::max(), '\n' );
+ }
+ else if (whatToDo.compare(0,19,"RANSACMaxIteration=") == 0) {
+ RansacParamFile >> DataBuffer;
+ fRANSACMaxIteration = atoi(DataBuffer.c_str());
+ cout << "/// RANSAC MaxIteration= " << " " << fRANSACMaxIteration << " ///" << endl;
+ }
+ else if (whatToDo.compare(0,15,"RANSACDistance=") == 0) {
+ RansacParamFile >> DataBuffer;
+ fRANSACDistance = atoi(DataBuffer.c_str());
+ cout << "/// RANSAC Distance= " << " " << fRANSACDistance << " ///" << endl;
+ }
+ else if (whatToDo.compare(0,22,"RANSACChargeThreshold=") == 0) {
+ RansacParamFile >> DataBuffer;
+ fRANSACChargeThreshold = atof(DataBuffer.c_str());
+ cout << "/// RANSAC ChargeThreshold= " << " " << fRANSACChargeThreshold << " ///" << endl;
+ }
+ else if (whatToDo.compare(0,21,"RANSACPointThreshold=") == 0) {
+ RansacParamFile >> DataBuffer;
+ fRANSACPointThreshold = atof(DataBuffer.c_str());
+ cout << "/// RANSAC PointThreshold= " << " " << fRANSACPointThreshold << " ///" << endl;
+ }
+ else if (whatToDo.compare(0,16,"RANSACThreshold=") == 0) {
+ RansacParamFile >> DataBuffer;
+ fRANSACThreshold = atoi(DataBuffer.c_str());
+ cout << "/// RANSAC Threshold= " << " " << fRANSACThreshold << " ///" << endl;
+ }
+ else if (whatToDo.compare(0,22,"MAXBarycenterDistance=") == 0) {
+ RansacParamFile >> DataBuffer;
+ fMAXBarycenterDistance = atoi(DataBuffer.c_str());
+ cout << "/// Max Barycenter Distance= " << " " << fMAXBarycenterDistance << " ///" << endl;
+ }
+ else if (whatToDo.compare(0,16,"AngleMaxToMerge=") == 0) {
+ RansacParamFile >> DataBuffer;
+ fAngleMax = atoi(DataBuffer.c_str());
+ cout << "/// Angle Max to merge tracks= " << " " << fAngleMax << " ///" << endl;
+ }
+ else if (whatToDo.compare(0,20,"MergeWithBarycenter=") == 0) {
+ RansacParamFile >> DataBuffer;
+ fMergeWithBarycenter = atoi(DataBuffer.c_str());
+ cout << "/// Merging Tracks using the barycenter= " << " " << fMergeWithBarycenter << " ///" << endl;
+ }
+ else if (whatToDo.compare(0,23,"MergeWithScalarProduct=") == 0) {
+ RansacParamFile >> DataBuffer;
+ fMergeWithScalarProduct = atoi(DataBuffer.c_str());
+ cout << "/// Merging Tracks using the scalar product= " << " " << fMergeWithScalarProduct << " ///" << endl;
+ }
+ else{
+ ReadingStatus = false;
+ }
+ }
+ }
+ }
+}
+
+//////////////////////////////////////////////////
+void Ransac::Init(vector<int> v1, vector<int> v2, vector<double> v3, vector<double> v4)
+{
+ Reset();
+
+ vX = v1;
+ vY = v2;
+ vZ = v3;
+ vQ = v4;
+
+ fOriginalCloudSize = vX.size();
+ double TotalCharge=0;
+ for(unsigned int i=0; i< vQ.size(); i++){
+ TotalCharge += vQ[i];
+ }
+ fTotalCharge = TotalCharge;
+
+}
+
+//////////////////////////////////////////////////
+vector<NPL::Track> Ransac::SimpleRansac()
+{
+ if(fVisu==1){
+ m_Histo->Clear();
+ }
+
+ TRandom* Rand=new TRandom();
+ double RemainingCharge = fTotalCharge;
+ /*cout << "/// Original total charge= " << fTotalCharge << endl;
+ cout << "/// Charge Threshold= " << fTotalCharge*fRANSACChargeThreshold << endl;
+ cout << "/// vX.size()= " << vX.size() << endl;
+ cout << "/// Point Threshold= " << fOriginalCloudSize*fRANSACPointThreshold << endl;*/
+
+ int aa=0;
+ //while(vX.size() > fOriginalCloudSize*fRANSACPointThreshold){
+ //while(RemainingCharge > fTotalCharge*fRANSACChargeThreshold){
+ while(RemainingCharge > fTotalCharge*fRANSACChargeThreshold && vX.size() > fOriginalCloudSize*fRANSACPointThreshold){
+ aa++;
+
+ double minErr=1E20;
+ trackX.clear();
+ trackY.clear();
+ trackZ.clear();
+ trackQ.clear();
+
+ TVector3 V1, V2;
+ std::vector<int> inliners;
+ inliners.clear();
+
+ fRANSACMaxIteration = 2.0*vX.size();
+ for(int i=0;i<fRANSACMaxIteration;i++){
+ int vXsize = vX.size();
+ int p1=(int)(Rand->Uniform(0,vXsize));
+ int p2;
+ do p2=(int)(Rand->Uniform(0,vXsize));
+ while(p2==p1);
+
+ TVector3 Vu;
+ Vu.SetXYZ(vX[p2]-vX[p1],vY[p2]-vY[p1],vZ[p2]-vZ[p1]);
+
+ std::vector<int> alsoin;
+ alsoin.clear();
+ for(int p=0;p<vXsize;p++){
+ if(p!=p1 && p!=p2 && vQ[p]){
+ TVector3 Va;
+ Va.SetXYZ(vX[p2]-vX[p],vY[p2]-vY[p],vZ[p2]-vZ[p]);
+ TVector3 VaCVu=Va.Cross(Vu);
+
+ if(VaCVu.Mag()/Vu.Mag()<fRANSACDistance){
+ alsoin.push_back(p);
+ }
+ }
+ }
+ //cout << "****** " << vX.size() << " / " << alsoin.size() << endl;
+
+ if(alsoin.size()>fRANSACThreshold){
+ TVector3 v1, v2;
+ double chi2=Fit3D(vX,vY,vZ,vQ,alsoin,v1,v2);
+ if(chi2<minErr){
+ minErr=chi2;
+ V1=v1;
+ V2=v2;
+ inliners=alsoin;
+ }
+ }
+ }
+
+ if(inliners.size()==0) break;
+ //cout << "iteration " << aa << ": with also inliners: " << inliners.size() << endl;
+
+ TVector3 Vdir = TVector3(V2.x()-V1.x(),V2.y()-V1.y(),V2.z()-V1.z());
+ Vdir = Vdir.Unit();
+
+ NPL::Track myTrack = NPL::Track();
+
+ myTrack.Xm=V1.x();
+ myTrack.Ym=V1.y();
+ myTrack.Zm=V1.z();
+
+ myTrack.Xh=V2.x();
+ myTrack.Yh=V2.y();
+ myTrack.Zh=V2.z();
+
+
+
+ int TrackNumber = vTrack.size();
+ TString TrackName = Form("Track%d",TrackNumber);
+
+ int inliner_size = inliners.size();
+ for(int p=0; p<inliner_size; p++){
+ trackX.push_back(vX[inliners[p]]);
+ trackY.push_back(vY[inliners[p]]);
+ trackZ.push_back(vZ[inliners[p]]);
+ trackQ.push_back(vQ[inliners[p]]);
+ }
+
+ myTrack.SetXPoints(trackX);
+ myTrack.SetYPoints(trackY);
+ myTrack.SetZPoints(trackZ);
+ myTrack.SetQPoints(trackQ);
+
+ vTrack.push_back(myTrack);
+
+ for(int i=inliner_size-1 ; i>=0 ; i--){
+ vQ.erase(vQ.begin()+inliners[i]);
+ vX.erase(vX.begin()+inliners[i]);
+ vY.erase(vY.begin()+inliners[i]);
+ vZ.erase(vZ.begin()+inliners[i]);
+ }
+
+ RemainingCharge = 0;
+ for(unsigned int i =0; i<vQ.size(); i++){
+ RemainingCharge += vQ[i];
+ }
+ //cout << "/// RemainingCharge= " << RemainingCharge << endl;
+ //cout << "////****** " << vX.size() << " / " << fOriginalCloudSize*fRANSACPointThreshold << endl;
+ //cout << "////****** " << RemainingCharge << " / " << fTotalCharge*fRANSACChargeThreshold << endl;
+ }
+ //cout << " //// Number of Tracks found " << vTrack.size() << " ////" << endl;
+
+ if(fMergeWithBarycenter) DoMerge("barycenter");
+ if(fMergeWithScalarProduct) DoMerge("scalar");
+
+ if(fVisu==1){
+ int track_size = vTrack.size();
+ for(int t=0; t<track_size; t++){
+ h3D = new TGraph2D();
+ for(int i=0; i<vTrack[t].GetXPoints().size(); i++){
+ h3D->SetPoint(i, vTrack[t].GetXPoints()[i], vTrack[t].GetYPoints()[i], vTrack[t].GetZPoints()[i]);
+ }
+ m_Histo->Add(h3D);
+
+ pl = new TGraph2D();
+
+ TVector3 Vdir = TVector3(vTrack[t].Xh-vTrack[t].Xm, vTrack[t].Yh-vTrack[t].Ym, vTrack[t].Zh-vTrack[t].Zm);
+ Vdir = Vdir.Unit();
+
+ //pl->SetName(TrackName);
+ pl->SetPoint(0,vTrack[t].Xm, vTrack[t].Ym, vTrack[t].Zm);
+ pl->SetPoint(1,vTrack[t].Xh, vTrack[t].Yh, vTrack[t].Zh);
+ pl->SetPoint(2,vTrack[t].Xm + 10*Vdir.x(), vTrack[t].Ym + 10*Vdir.y(), (vTrack[t].Zm + 10*Vdir.z()));
+ pl->SetPoint(3,vTrack[t].Xm + 100*Vdir.x(), vTrack[t].Ym + 100*Vdir.y(), (vTrack[t].Zm + 100*Vdir.z()));
+ pl->SetPoint(4,vTrack[t].Xm - 100*Vdir.x(), vTrack[t].Ym - 100*Vdir.y(), (vTrack[t].Zm - 100*Vdir.z()));
+ pl->SetLineColor(2);
+ pl->SetMarkerColor(2);
+ pl->SetMarkerStyle(8);
+ pl->SetMarkerSize(1);
+ pl->SetLineWidth(2);
+ m_Histo->Add(pl);
+ }
+ }
+
+
+
+ if(fVisu && vTrack.size()>1 && vTrack.size()<4){
+ TSocket* s = NULL;
+ char request[64];
+ if(m_ServerSocket && m_Monitor){
+
+ m_Monitor->ResetInterrupt();
+ TList readClients;
+ if((m_Monitor->Select(&readClients, NULL, 10)) > 0){
+ s = (TSocket*) readClients.Last();
+ readClients.RemoveLast();
+ }
+ }
+
+ if(s){
+ if (s->IsA() == TServerSocket::Class()) {
+ // accept new connection from spy
+ TSocket* socket = ((TServerSocket*)s)->Accept();
+ m_Monitor->Add(socket,TMonitor::kRead|TMonitor::kWrite);
+ m_Sockets->Add(socket);
+ cout << "Client Connected" << endl;
+ //NPL::SendInformation("NPL::SpectraServer","Accepted new client connection");
+ }
+ else{
+ // we only get string based requests from the spy
+ char request[64];
+ if (s->Recv(request, sizeof(request)) <= 0) {
+ m_Monitor->Remove(s);
+ m_Sockets->Remove(s);
+ delete s;
+ //return;
+ }
+ // send requested object back
+ static TMessage answer(kMESS_OBJECT);
+ answer.Reset();
+
+ if (!strcmp(request, "RequestHisto")){
+ answer.WriteObject(m_Histo);
+ s->Send(answer);
+ }
+
+ else // answer nothing
+ s->Send(answer);
+ }
+ }
+ }
+
+ return vTrack;
+}
+
+//////////////////////////////////////////////////
+bool Ransac::CompareTracksWithBarycenter(NPL::Track track1, NPL::Track track2)
+{
+ bool ToBeMerged = false;
+
+ TVector3 bar1 = track1.GetPointBarycenter();
+ TVector3 bar2 = track2.GetPointBarycenter();
+ double Y1 = bar1.Y();
+ double Z1 = bar1.Z();
+ double Y2 = bar2.Y();
+ double Z2 = bar2.Z();
+
+ double Distance = sqrt(pow(Y2-Y1,2) + pow(Z2-Z1,2));
+
+ if(Distance<fMAXBarycenterDistance) ToBeMerged=true;
+
+ return ToBeMerged;
+}
+
+//////////////////////////////////////////////////
+bool Ransac::CompareTracksWithScalarProduct(NPL::Track track1, NPL::Track track2)
+{
+ bool ToBeMerged = false;
+
+ TVector3 v1 = TVector3(track1.Xh-track1.Xm, track1.Zh-track1.Zm, track1.Zh-track1.Zm);
+ TVector3 v2 = TVector3(track2.Xh-track2.Xm, track2.Zh-track2.Zm, track2.Zh-track2.Zm);
+
+ TVector3 v1_unit = v1.Unit();
+ TVector3 v2_unit = v2.Unit();
+
+ double ScalarProduct = v1_unit.Dot(v2_unit);
+
+ if(abs(ScalarProduct)>cos(fAngleMax*TMath::Pi()/180)) ToBeMerged=true;
+
+ return ToBeMerged;
+}
+
+//////////////////////////////////////////////////
+NPL::Track Ransac::MergeTracks(vector<NPL::Track> tracks)
+{
+ NPL::Track newTrack = NPL::Track();
+ vector<int> newTrackX, newTrackY, newTrackZ, newTrackQ;
+
+ double X1=0; double Y1=0; double Z1=0;
+ double X2=0; double Y2=0; double Z2=0;
+
+ for(int j=0; j<tracks.size(); j++){
+ X1 += tracks[j].Xm;
+ Y1 += tracks[j].Ym;
+ Z1 += tracks[j].Zm;
+
+ X2 += tracks[j].Xh;
+ Y2 += tracks[j].Yh;
+ Z2 += tracks[j].Zh;
+
+ for(int k=0; k<tracks[j].GetXPoints().size(); k++){
+ newTrackX.push_back(tracks[j].GetXPoints()[k]);
+ newTrackY.push_back(tracks[j].GetYPoints()[k]);
+ newTrackZ.push_back(tracks[j].GetZPoints()[k]);
+ }
+
+ }
+ X1 = X1/tracks.size();
+ Y1 = Y1/tracks.size();
+ Z1 = Z1/tracks.size();
+ X2 = X2/tracks.size();
+ Y2 = Y2/tracks.size();
+ Z2 = Z2/tracks.size();
+
+ newTrack.Xm=X1;
+ newTrack.Ym=Y1;
+ newTrack.Zm=Z1;
+
+ newTrack.Xh=X2;
+ newTrack.Yh=Y2;
+ newTrack.Zh=Z2;
+
+ newTrack.SetXPoints(newTrackX);
+ newTrack.SetYPoints(newTrackY);
+ newTrack.SetZPoints(newTrackZ);
+
+ return newTrack;
+}
+
+//////////////////////////////////////////////////
+void Ransac::DoMerge(string method)
+{
+ vector<NPL::Track> vTrack_new;
+ vector<NPL::Track> vTrackToMerge;
+ vector<NPL::Track> vTrack_temp;
+ vTrack_temp = vTrack;
+ NPL::Track newTrack = NPL::Track();
+
+ int size = vTrack.size();
+ vector<int> id_track;
+ for(int i=0; i<size-1; i++){
+ int NumberOfTrackToBeMergedWith=0;
+ vTrackToMerge.clear();
+
+ for(int j=i+1; j<size; j++){
+ if(find(id_track.begin(),id_track.end(), j) == id_track.end()){
+ if(method=="scalar") {
+ if( CompareTracksWithScalarProduct(vTrack[i],vTrack[j]) ){
+ NumberOfTrackToBeMergedWith++;
+ id_track.push_back(j);
+ vTrackToMerge.push_back(vTrack[j]);
+ }
+ }
+ else if(method=="barycenter") {
+ if( CompareTracksWithBarycenter(vTrack[i],vTrack[j]) ){
+ NumberOfTrackToBeMergedWith++;
+ id_track.push_back(j);
+ vTrackToMerge.push_back(vTrack[j]);
+ }
+ }
+ else cout << "Method for merging not known: Problem with token!" << endl;
+ }
+ }
+ if(NumberOfTrackToBeMergedWith>0){
+ if(find(id_track.begin(),id_track.end(), i) == id_track.end()) id_track.push_back(i);
+ vTrackToMerge.push_back(vTrack[i]);
+
+ newTrack = MergeTracks(vTrackToMerge);
+
+ vTrack_new.push_back(newTrack);
+ }
+ }
+
+ std::sort(id_track.begin(), id_track.end());
+ for(int p=id_track.size()-1; p>=0; p--){
+ vTrack_temp.erase(vTrack_temp.begin()+id_track[p]);
+ }
+
+ for(int p=0; p<vTrack_temp.size(); p++){
+ vTrack_new.push_back(vTrack_temp[p]);
+ }
+
+ vTrack.clear();
+ vTrack = vTrack_new;
+}
+
+
+//////////////////////////////////////////////////
+/*void Ransac::CheckTracks()
+{
+ vector<NPL::Track> vTrack_temp;
+ vTrack_temp = vTrack;
+
+ int size = vTrack.size();
+ vector<double> Ybar;
+ vector<double> Zbar;
+
+ Ybar.clear();
+ Zbar.clear();
+ for(int i=0; i<size; i++){
+ TVector3 bar = vTrack[i].GetPointBarycenter();
+ Ybar.push_back(bar.Y());
+ Zbar.push_back(bar.Z());
+ }
+
+ vector<double> Distance;
+ vector<int> id_i, id_j;
+ for(int i=0; i<size-1; i++){
+ for(int j=i+1; j<size; j++){
+ id_i.push_back(i);
+ id_j.push_back(j);
+ Distance.push_back(sqrt(pow(Ybar[i]-Ybar[j],2) + pow(Zbar[i]-Zbar[j],2)));
+ }
+ }
+
+ vector<int> removed_id;
+ for(int c=0; c<Distance.size(); c++){
+ if(Distance[c]<fMAXBarycenterDistance){
+ if(find(removed_id.begin(),removed_id.end(), id_i[c]) == removed_id.end()) removed_id.push_back(id_i[c]);
+ if(find(removed_id.begin(),removed_id.end(), id_j[c]) == removed_id.end()) removed_id.push_back(id_j[c]);
+ }
+ }
+
+ std::sort (removed_id.begin(), removed_id.end());
+ for(int j=removed_id.size()-1; j>=0; j--){
+ vTrack.erase(vTrack.begin()+removed_id[j]);
+ }
+
+
+ double X1=0; double Y1=0; double Z1=0;
+ double X2=0; double Y2=0; double Z2=0;
+
+ NPL::Track newTrack = NPL::Track();
+ vector<int> newTrackX, newTrackY, newTrackZ, newTrackQ;
+
+ if(removed_id.size()>0){
+ for(int j=0; j<removed_id.size(); j++){
+ X1 += vTrack_temp[removed_id[j]].Xm;
+ Y1 += vTrack_temp[removed_id[j]].Ym;
+ Z1 += vTrack_temp[removed_id[j]].Zm;
+
+ X2 += vTrack_temp[removed_id[j]].Xh;
+ Y2 += vTrack_temp[removed_id[j]].Yh;
+ Z2 += vTrack_temp[removed_id[j]].Zh;
+
+ for(int k=0; k<vTrack_temp[removed_id[j]].GetXPoints().size(); k++){
+ newTrackX.push_back(vTrack_temp[removed_id[j]].GetXPoints()[k]);
+ newTrackY.push_back(vTrack_temp[removed_id[j]].GetYPoints()[k]);
+ newTrackZ.push_back(vTrack_temp[removed_id[j]].GetZPoints()[k]);
+ }
+
+ }
+ X1 = X1/removed_id.size();
+ Y1 = Y1/removed_id.size();
+ Z1 = Z1/removed_id.size();
+ X2 = X2/removed_id.size();
+ Y2 = Y2/removed_id.size();
+ Z2 = Z2/removed_id.size();
+
+ newTrack.Xm=X1;
+ newTrack.Ym=Y1;
+ newTrack.Zm=Z1;
+
+ newTrack.Xh=X2;
+ newTrack.Yh=Y2;
+ newTrack.Zh=Z2;
+
+ newTrack.SetXPoints(newTrackX);
+ newTrack.SetYPoints(newTrackY);
+ newTrack.SetZPoints(newTrackZ);
+
+ vTrack.push_back(newTrack);
+ }
+}*/
+
+//////////////////////////////////////////////////
+vector<double> Ransac::GetChargeOfTracks()
+{
+ return vTrackCharge;
+}
+
+//////////////////////////////////////////////////
+vector<double> Ransac::GetTrackLength(double PadSizeX, double PadSizeY, double DriftVelocity)
+{
+ vector<double> length;
+ double Ymin, Ymax;
+ double Zmin, Zmax;
+ double Xmin = 240;
+ double Xmax = 256;
+
+ for(unsigned int i=0; i<vTrack.size(); i++){
+ double xh = vTrack[i].GetXh()*PadSizeX;
+ double xm = vTrack[i].GetXm()*PadSizeX;
+ double yh = vTrack[i].GetYh()*PadSizeY;
+ double ym = vTrack[i].GetYm()*PadSizeY;
+ double zh = vTrack[i].GetZh()*DriftVelocity;
+ double zm = vTrack[i].GetZm()*DriftVelocity;
+
+ Ymin = yh + (ym - yh)*(Xmin-xh)/(xm-xh);
+ Ymax = yh + (ym - yh)*(Xmax-xh)/(xm-xh);
+
+ Zmin = zh + (zm - zh)*(Xmin-xh)/(xm-xh);
+ Zmax = zh + (zm - zh)*(Xmax-xh)/(xm-xh);
+
+ length.push_back(sqrt( pow(Xmax-Xmin,2) + pow(Ymax-Ymin,2) + pow(Zmax-Zmin,2) ));
+ }
+
+ return length;
+}
+
+
+//////////////////////////////////////////////////
+void Ransac::Reset()
+{
+ vline.clear();
+ vX.clear();
+ vY.clear();
+ vZ.clear();
+ vQ.clear();
+ vTrack.clear();
+ vTrackCharge.clear();
+}
+
+//////////////////////////////////////////////////
+double Ransac::Fit3D(vector<int> X, vector<int> Y, vector<double> Z, vector<double> Charge, vector<int> inliners, TVector3& V1, TVector3& V2)
+{
+ // adapted from: http://fr.scribd.com/doc/31477970/Regressions-et-trajectoires-3D
+
+ int R, C;
+ double Q;
+ double Xm,Ym,Zm;
+ double Xh,Yh,Zh;
+ double a,b;
+ double Sxx,Sxy,Syy,Sxz,Szz,Syz;
+ double theta;
+ double K11,K22,K12,K10,K01,K00;
+ double c0,c1,c2;
+ double p,q,r,dm2;
+ double rho,phi;
+
+ Q=Xm=Ym=Zm=0.;
+ double total_charge=0;
+ Sxx=Syy=Szz=Sxy=Sxz=Syz=0.;
+
+ for (auto i : inliners)
+ {
+ if(X[i]>119)total_charge+=Charge[i];
+ Q+=Charge[i]/10.;
+ Xm+=X[i]*Charge[i]/10.;
+ Ym+=Y[i]*Charge[i]/10.;
+ Zm+=Z[i]*Charge[i]/10.;
+ Sxx+=X[i]*X[i]*Charge[i]/10.;
+ Syy+=Y[i]*Y[i]*Charge[i]/10.;
+ Szz+=Z[i]*Z[i]*Charge[i]/10.;
+ Sxy+=X[i]*Y[i]*Charge[i]/10.;
+ Sxz+=X[i]*Z[i]*Charge[i]/10.;
+ Syz+=Y[i]*Z[i]*Charge[i]/10.;
+ }
+ vTrackCharge.push_back(total_charge);
+
+ Xm/=Q;
+ Ym/=Q;
+ Zm/=Q;
+ Sxx/=Q;
+ Syy/=Q;
+ Szz/=Q;
+ Sxy/=Q;
+ Sxz/=Q;
+ Syz/=Q;
+ Sxx-=(Xm*Xm);
+ Syy-=(Ym*Ym);
+ Szz-=(Zm*Zm);
+ Sxy-=(Xm*Ym);
+ Sxz-=(Xm*Zm);
+ Syz-=(Ym*Zm);
+
+ theta=0.5*atan((2.*Sxy)/(Sxx-Syy));
+
+ K11=(Syy+Szz)*pow(cos(theta),2)+(Sxx+Szz)*pow(sin(theta),2)-2.*Sxy*cos(theta)*sin(theta);
+ K22=(Syy+Szz)*pow(sin(theta),2)+(Sxx+Szz)*pow(cos(theta),2)+2.*Sxy*cos(theta)*sin(theta);
+ K12=-Sxy*(pow(cos(theta),2)-pow(sin(theta),2))+(Sxx-Syy)*cos(theta)*sin(theta);
+ K10=Sxz*cos(theta)+Syz*sin(theta);
+ K01=-Sxz*sin(theta)+Syz*cos(theta);
+ K00=Sxx+Syy;
+
+ c2=-K00-K11-K22;
+ c1=K00*K11+K00*K22+K11*K22-K01*K01-K10*K10;
+ c0=K01*K01*K11+K10*K10*K22-K00*K11*K22;
+
+
+ p=c1-pow(c2,2)/3.;
+ q=2.*pow(c2,3)/27.-c1*c2/3.+c0;
+ r=pow(q/2.,2)+pow(p,3)/27.;
+
+
+ if(r>0) dm2=-c2/3.+pow(-q/2.+sqrt(r),1./3.)+pow(-q/2.-sqrt(r),1./3.);
+ if(r<0)
+ {
+ rho=sqrt(-pow(p,3)/27.);
+ phi=acos(-q/(2.*rho));
+ dm2=min(-c2/3.+2.*pow(rho,1./3.)*cos(phi/3.),min(-c2/3.+2.*pow(rho,1./3.)*cos((phi+2.*TMath::Pi())/3.),-c2/3.+2.*pow(rho,1./3.)*cos((phi+4.*TMath::Pi())/3.)));
+ }
+
+ a=-K10*cos(theta)/(K11-dm2)+K01*sin(theta)/(K22-dm2);
+ b=-K10*sin(theta)/(K11-dm2)-K01*cos(theta)/(K22-dm2);
+
+ Xh=((1.+b*b)*Xm-a*b*Ym+a*Zm)/(1.+a*a+b*b);
+ Yh=((1.+a*a)*Ym-a*b*Xm+b*Zm)/(1.+a*a+b*b);
+ Zh=((a*a+b*b)*Zm+a*Xm+b*Ym)/(1.+a*a+b*b);
+
+ V1.SetXYZ(Xm,Ym,Zm);
+ V2.SetXYZ(Xh,Yh,Zh);
+
+ return(fabs(dm2/Q));
+}
diff --git a/NPLib/TrackReconstruction/NPRansac.h b/NPLib/TrackReconstruction/NPRansac.h
new file mode 100644
index 000000000..2428c5571
--- /dev/null
+++ b/NPLib/TrackReconstruction/NPRansac.h
@@ -0,0 +1,117 @@
+#ifndef __RANSAC__
+#define __RANSAC__
+/*****************************************************************************
+ * Copyright (C) 2009-2016 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * *
+ * Original Author : Pierre MORFOUACE contact address: morfouace@ganil.fr *
+ * *
+ * Creation Date : April 2018 *
+ * Last update : April 2018 *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class deal with finding all the track event by event *
+ *****************************************************************************/
+//C++ Header
+#include <stdio.h>
+#include <vector>
+#include <sstream>
+#include <iostream>
+#include <fstream>
+#include <cmath>
+#include <stdlib.h>
+#include <limits>
+#include <string>
+#include <algorithm>
+
+//NPL
+#include "NPTrack.h"
+using namespace NPL;
+
+// ROOT Headers
+#include <TLine.h>
+#include <TH2F.h>
+#include <TH3F.h>
+#include <TGraph2D.h>
+#include <TMath.h>
+#include <TCanvas.h>
+#include <TVector3.h>
+#include <TRandom.h>
+#include <TServerSocket.h>
+#include <TSocket.h>
+#include <TMonitor.h>
+#include <TMessage.h>
+#include <TList.h>
+#include <TVector.h>
+
+using namespace std;
+
+namespace NPL{
+
+ class Ransac
+ {
+ public:
+ Ransac();
+ Ransac(int NumberOfPadsX, int NumberOfPadsY, bool Visu);
+ ~Ransac();
+
+ public:
+ void Reset();
+ void ReadParameterValue(string filename);
+ void Init(vector<int> v1, vector<int> v2, vector<double> v3, vector<double> v4);
+ vector<NPL::Track> SimpleRansac();
+ bool CompareTracksWithBarycenter(NPL::Track track1, NPL::Track track2);
+ bool CompareTracksWithScalarProduct(NPL::Track track1, NPL::Track track2);
+ NPL::Track MergeTracks(vector<NPL::Track> tracks);
+ void DoMerge(string method);
+ //void CheckTracks();
+
+ vector<double> GetChargeOfTracks();
+ vector<double> GetTrackLength(double PadSizeX, double PadSizeY, double DriftVelocity);
+ double Fit3D(vector<int> X, vector<int> Y, vector<double> Z, vector<double> Charge, vector<int> inliners, TVector3& V1, TVector3& V2);
+
+ private:
+ TCanvas* c1;
+ vector<TLine*> vline;
+ vector<int> vX, vY;
+ vector<int> trackX, trackY, trackZ, trackQ;
+ vector<double> vZ, vQ;
+ TLine* L;
+ TGraph2D* h3D;
+ TGraph2D* pl;
+ vector<NPL::Track> vTrack;
+ vector<double> vTrackCharge;
+
+ private:
+ float fRANSACThreshold;
+ float fRANSACPointThreshold;
+ float fRANSACChargeThreshold;
+ float fRANSACDistance;
+ float fRANSACMaxIteration;
+ float fMAXBarycenterDistance;
+ float fAngleMax;
+ int fNumberOfTracksMax;
+ int fOriginalCloudSize;
+ double fTotalCharge;
+ int fNumberOfPadsX;
+ int fNumberOfPadsY;
+ int fVisu;
+ bool fMergeWithBarycenter;
+ bool fMergeWithScalarProduct;
+
+ private:
+ TServerSocket* m_ServerSocket;
+ TSocket* m_Socket;
+ TMonitor* m_Monitor;
+ TList* m_Histo;
+ TList* m_Sockets;
+
+ ClassDef(Ransac, 1)
+ };
+}
+#endif
diff --git a/NPLib/TrackReconstruction/NPTrack.cxx b/NPLib/TrackReconstruction/NPTrack.cxx
new file mode 100644
index 000000000..93b00e8e2
--- /dev/null
+++ b/NPLib/TrackReconstruction/NPTrack.cxx
@@ -0,0 +1,210 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2016 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * *
+ * Original Author : Pierre MORFOUACE contact address: morfouace@ganil.fr *
+ * *
+ * Creation Date : April 2018 *
+ * Last update : Apro; 2018 *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class deal with track charactieristic from TPC analysis *
+ *****************************************************************************/
+
+
+#include "NPTrack.h"
+using namespace NPL;
+
+using namespace std;
+
+//ClassImp(NPTrack)
+
+//////////////////////////////////////////////////////
+Track::Track()
+{
+ /*L2DXY=new TLine();
+ L2DXZ=new TLine();
+ L2DYZ=new TLine();
+ L3D=new TLine();*/
+}
+
+//////////////////////////////////////////////////////
+Track::~Track()
+{
+ /*delete L2DXY;
+ delete L2DXZ;
+ delete L2DYZ;
+ delete L3D;*/
+}
+
+//////////////////////////////////////////////////////
+TVector3 Track::GetDirectionVector()
+{
+ TVector3 vTrack = TVector3(Xh-Xm, Yh-Ym, Zh-Zm);
+ return vTrack;
+}
+
+//////////////////////////////////////////////////////
+TVector3 Track::GetVertexPostion(TVector3 vBeamDir, TVector3 vBeamPoint)
+{
+ //y_beam = Ay_beam*x_beam + By_beam
+ double Ay_beam = vBeamDir.Y()/vBeamDir.X();
+ double By_beam = vBeamPoint.Y() - (vBeamDir.Y()/vBeamDir.X())*vBeamPoint.X();
+ //z_beam = Az_beam*x_beam + Bz_beam
+ double Az_beam = vBeamDir.Z()/vBeamDir.X();
+ double Bz_beam = vBeamPoint.Z() - (vBeamDir.Z()/vBeamDir.X())*vBeamPoint.X();
+
+ //y_track = Ay_track*x_track + By_track
+ double Ay_track = GetDirectionVector().Y()/GetDirectionVector().X();
+ double By_track = Ym - (GetDirectionVector().Y()/GetDirectionVector().X())*Xm;
+ //z_track = Az_track*x_track + Bz_track
+ double Az_track = GetDirectionVector().Z()/GetDirectionVector().X();
+ double Bz_track = Zm - (GetDirectionVector().Z()/GetDirectionVector().X())*Xm;
+
+ double DBy = By_beam - By_track;
+ double DBz = Bz_beam - Bz_track;
+ double alpha = (-Ay_beam*DBy - Az_beam*DBz)/(Ay_beam*Ay_beam + Az_beam*Az_beam + 1);
+ double beta = (Ay_beam*Ay_track + Az_beam*Az_track + 1)/(Ay_beam*Ay_beam + Az_beam*Az_beam + 1);
+
+ double A = beta*(Ay_beam*Ay_beam + Az_beam*Az_beam + 1) - (Ay_beam*Ay_track + Az_beam*Az_track + 1);
+ double B = (Ay_track*Ay_track + Az_track*Az_track + 1) - beta*(Ay_beam*Ay_track + Az_beam*Az_track + 1);
+ double C = beta*(Ay_beam*DBy + Az_beam*DBz) - Ay_track*DBy - Az_track*DBz;
+
+ double xt = -(A*alpha+C)/(A*beta+B);
+ double xb = beta*xt+alpha;
+
+ double yb = Ay_beam*xb + By_beam;
+ double zb = Az_beam*xb + Bz_beam;
+
+ double yt = Ay_track*xt + By_track;
+ double zt = Az_track*xt + Bz_track;
+
+ double xvertex = (xb+xt)/2;
+ double yvertex = (yb+yt)/2;
+ double zvertex = (zb+zt)/2;
+
+ TVector3 vVertex = TVector3(xvertex,yvertex,zvertex);
+
+ return vVertex;
+}
+
+//////////////////////////////////////////////////////
+TVector3 Track::GetChargeBarycenter()
+{
+ double bX=0;
+ double bY=0;
+ double bZ=0;
+ double Qtot=0;
+
+ int size = vX.size();
+
+ for(int i=0; i<size; i++){
+ Qtot += (double)vQ[i];
+ bX += (double)vX[i]*vQ[i];
+ bY += (double)vY[i]*vQ[i];
+ bZ += (double)vZ[i]*vQ[i];
+ }
+
+ bX=bX/Qtot;
+ bY=bY/Qtot;
+ bZ=bZ/Qtot;
+
+ TVector3 vBarycentre = TVector3(bX,bY,bZ);
+
+ return vBarycentre;
+}
+
+//////////////////////////////////////////////////////
+TVector3 Track::GetPointBarycenter()
+{
+ double bX=0;
+ double bY=0;
+ double bZ=0;
+
+ int size = vX.size();
+
+ for(int i=0; i<size; i++){
+ bX += (double)vX[i];
+ bY += (double)vY[i];
+ bZ += (double)vZ[i];
+ }
+
+ bX=bX/size;
+ bY=bY/size;
+ bZ=bZ/size;
+
+ TVector3 vBarycentre = TVector3(bX,bY,bZ);
+
+ return vBarycentre;
+}
+
+//////////////////////////////////////////////////////
+double Track::GetTrackLength(double PadSizeX, double PadSizeY, double DriftVelocity)
+{
+ double length;
+
+ double Xmax = *max_element(vX.begin(), vX.end())*PadSizeX;
+ double Xmin = *min_element(vX.begin(), vX.end())*PadSizeX;
+
+ double Ymax = *max_element(vY.begin(), vY.end())*PadSizeY;
+ double Ymin = *min_element(vY.begin(), vY.end())*PadSizeY;
+
+ double Zmax = *max_element(vZ.begin(), vZ.end())*DriftVelocity;
+ double Zmin = *min_element(vZ.begin(), vZ.end())*DriftVelocity;
+
+ length = sqrt(pow(Xmax-Xmin,2) + pow(Ymax-Ymin,2) + pow(Zmax-Zmin,2));
+
+ return length;
+}
+
+//////////////////////////////////////////////////////
+double Track::GetTotalCharge()
+{
+ double Qtot=0;
+
+ int size = vQ.size();
+ for(int i=0; i<size; i++){
+ Qtot += vQ[i];
+ }
+ Qtot = Qtot/size;
+
+ return Qtot;
+}
+
+//////////////////////////////////////////////////////
+void Track::Clear()
+{
+ vX.clear();
+ vY.clear();
+ vZ.clear();
+ vQ.clear();
+}
+
+//////////////////////////////////////////////////////
+/*void Track::ResetLines()
+{
+ L2DXY->SetX1(-1);
+ L2DXY->SetY1(-1);
+ L2DXY->SetX2(-1);
+ L2DXY->SetY2(-1);
+
+ L2DXZ->SetX1(-1);
+ L2DXZ->SetY1(-1);
+ L2DXZ->SetX2(-1);
+ L2DXZ->SetY2(-1);
+
+ L2DYZ->SetX1(-1);
+ L2DYZ->SetY1(-1);
+ L2DYZ->SetX2(-1);
+ L2DYZ->SetY2(-1);
+
+ L3D->SetX1(-1);
+ L3D->SetY1(-1);
+ L3D->SetX2(-1);
+ L3D->SetY2(-1);
+}*/
diff --git a/NPLib/TrackReconstruction/NPTrack.h b/NPLib/TrackReconstruction/NPTrack.h
new file mode 100644
index 000000000..1cb7226bd
--- /dev/null
+++ b/NPLib/TrackReconstruction/NPTrack.h
@@ -0,0 +1,96 @@
+#ifndef __TRACK__
+#define __TRACK__
+/*****************************************************************************
+ * Copyright (C) 2009-2016 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * *
+ * Original Author : Pierre MORFOUACE contact address: morfouace@ganil.fr *
+ * *
+ * Creation Date : April 2018 *
+ * Last update : Apro; 2018 *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class deal with track charactieristic from TPC analysis *
+ *****************************************************************************/
+
+#include <TLine.h>
+#include <TVector3.h>
+#include <stdio.h>
+#include <iostream>
+#include <vector>
+#include <algorithm>
+#include <cmath>
+
+using namespace std;
+
+namespace NPL{
+ class Track{
+
+ public:
+ Track();
+ ~Track();
+ public:
+ double Xm;
+ double Ym;
+ double Zm;
+ double Xh;
+ double Yh;
+ double Zh;
+
+ TLine* L2DXY;
+ TLine* L2DXZ;
+ TLine* L2DYZ;
+ TLine* L3D;
+
+ private:
+ vector<int> vX;
+ vector<int> vY;
+ vector<int> vZ;
+ vector<int> vQ;
+
+ public:
+ ///////////////
+ /// SETTERS ///
+ ///////////////
+ void SetXPoints(vector<int> vec) {vX=vec;}
+ void SetYPoints(vector<int> vec) {vY=vec;}
+ void SetZPoints(vector<int> vec) {vZ=vec;}
+ void SetQPoints(vector<int> vec) {vQ=vec;}
+
+ ///////////////
+ /// GETTERS ///
+ ///////////////
+ vector<int> GetXPoints() {return vX;}
+ vector<int> GetYPoints() {return vY;}
+ vector<int> GetZPoints() {return vZ;}
+ vector<int> GetQPoints() {return vQ;}
+
+ double GetXm() {return Xm;}
+ double GetXh() {return Xh;}
+ double GetYm() {return Ym;}
+ double GetYh() {return Yh;}
+ double GetZm() {return Zm;}
+ double GetZh() {return Zh;}
+
+ public:
+ TVector3 GetDirectionVector();
+
+ TVector3 GetVertexPostion(TVector3 vBeamDir, TVector3 vBeamPoint);
+
+ TVector3 GetChargeBarycenter();
+
+ TVector3 GetPointBarycenter();
+
+ double GetTrackLength(double PadSizeX, double PadSizeY, double DriftVelocity);
+ double GetTotalCharge();
+
+ //void ResetLines();
+ void Clear();
+ };
+}
+#endif
diff --git a/NPLib/TrackReconstruction/NPTrackReconstructionLinkDef.h b/NPLib/TrackReconstruction/NPTrackReconstructionLinkDef.h
new file mode 100644
index 000000000..058353da8
--- /dev/null
+++ b/NPLib/TrackReconstruction/NPTrackReconstructionLinkDef.h
@@ -0,0 +1,4 @@
+#ifdef __CINT__
+#pragma link C++ defined_in "./NPRansac.h";
+#pragma link C++ defined_in "./NPTrack.h";
+#endif
diff --git a/NPLib/scripts/build_dict.sh b/NPLib/scripts/build_dict.sh
index facd09f59..c8ee5388c 100755
--- a/NPLib/scripts/build_dict.sh
+++ b/NPLib/scripts/build_dict.sh
@@ -53,13 +53,13 @@ fi
# Version 5 : generate the dictionnary then the libmap
if [ $version_major -eq 5 ]
then
- rootcint -f $2 -c -I../Core -I../Physics -I../../Core -I../../Physics $1 $5
+ rootcint -f $2 -c -I../Core -I../Physics -I../../Core -I../../Physics -I../TrackReconstruction -I../../TrackReconstruction $1 $5
fi
# Version 6 or more : generate both at once
if [ $version_major -gt 5 ]
then
- rootcint -f $2 -rmf $3 -rml $lib_name -I../Core -I../Physics -I../../Core -I../../Physics $1 $5
+ rootcint -f $2 -rmf $3 -rml $lib_name -I../Core -I../Physics -I../../Core -I../../Physics -I../TrackReconstruction -I../../TrackReconstruction $1 $5
fi
diff --git a/NPSimulation/Core/MaterialManager.cc b/NPSimulation/Core/MaterialManager.cc
index 703c3dff5..7a3012331 100644
--- a/NPSimulation/Core/MaterialManager.cc
+++ b/NPSimulation/Core/MaterialManager.cc
@@ -821,18 +821,18 @@ G4Material* MaterialManager::GetGasFromLibrary(string Name, double Pressure, dou
if(Name == "H2"){
density = (2*1.00794/Vm)*mg/cm3;
- G4Material* material = new G4Material("NPS_"+newName,density,2,kStateGas,Temperature,Pressure);
- material->AddElement(GetElementFromLibrary("H"), 1);
- material->AddElement(GetElementFromLibrary("H"), 1);
+ G4Material* material = new G4Material("NPS_"+newName,density,1,kStateGas,Temperature,Pressure);
+ material->AddElement(GetElementFromLibrary("H"), 2);
+ //material->AddElement(GetElementFromLibrary("H"), 1);
m_Material[Name]=material;
return material;
}
if(Name == "D2"){
density = (2*2.0140/Vm)*mg/cm3;
- G4Material* material = new G4Material("NPS_"+newName,density,2,kStateGas,Temperature,Pressure);
- material->AddElement(GetElementFromLibrary("D"), 1);
- material->AddElement(GetElementFromLibrary("D"), 1);
+ G4Material* material = new G4Material("NPS_"+newName,density,1,kStateGas,Temperature,Pressure);
+ material->AddElement(GetElementFromLibrary("D"), 2);
+ //material->AddElement(GetElementFromLibrary("D"), 1);
m_Material[Name]=material;
return material;
}
@@ -850,7 +850,7 @@ G4Material* MaterialManager::GetGasFromLibrary(string Name, double Pressure, dou
void MaterialManager::WriteDEDXTable(G4ParticleDefinition* Particle ,G4double Emin,G4double Emax){
map<string,G4Material*>::iterator it;
if(Particle->GetPDGCharge()==0)
- return;
+ return;
for(it = m_Material.begin() ; it != m_Material.end() ; it++){
// Opening hte output file
G4String GlobalPath = getenv("NPTOOL");
@@ -894,7 +894,7 @@ void MaterialManager::WriteDEDXTable(std::set<string> Particle ,G4double Emin,G4
std::set<string>::iterator it;
for(it=Particle.begin(); it!=Particle.end() ; it++){
G4ParticleDefinition* p = G4ParticleTable::GetParticleTable()->FindParticle((*it));
- WriteDEDXTable(p,Emin,Emax);
+ WriteDEDXTable(p,Emin,Emax);
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/NPSimulation/Core/Target.cc b/NPSimulation/Core/Target.cc
index 35b795c4c..166407999 100644
--- a/NPSimulation/Core/Target.cc
+++ b/NPSimulation/Core/Target.cc
@@ -430,7 +430,7 @@ void Target::SetReactionRegion(){
if(!m_ReactionRegion){
m_ReactionRegion= new G4Region("NPSimulationProcess");
m_ReactionRegion->AddRootLogicalVolume(m_TargetLogic);
- m_ReactionRegion->SetUserLimits(new G4UserLimits(m_TargetThickness/10.));
+ m_ReactionRegion->SetUserLimits(new G4UserLimits(m_TargetThickness/10.));
}
G4FastSimulationManager* mng = m_ReactionRegion->GetFastSimulationManager();
@@ -447,7 +447,7 @@ void Target::SetReactionRegion(){
fsm = new NPS::Decay("Decay",m_ReactionRegion);
m_ReactionModel.push_back(fsm);
-}
+}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Add Detector branch to the EventTree.
diff --git a/NPSimulation/Detectors/Actar/Actar.cc b/NPSimulation/Detectors/Actar/Actar.cc
index 090dcc5b5..bbe53861b 100644
--- a/NPSimulation/Detectors/Actar/Actar.cc
+++ b/NPSimulation/Detectors/Actar/Actar.cc
@@ -1,18 +1,18 @@
/*****************************************************************************
- * Copyright (C) 2009-2017 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2017 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: Pierre Morfouace contact address: morfouac@nscl.msu.edu *
+ * Original Author: Pierre Morfouace contact address: morfouac@nscl.msu.edu *
* *
* Creation Date : September 2017 *
* Last update : *
*---------------------------------------------------------------------------*
* Decription: *
- * This class describe Actar simulation *
+ * This class describe Actar simulation *
* *
*---------------------------------------------------------------------------*
* Comment: *
@@ -53,7 +53,7 @@
// NPTool header
#include "Actar.hh"
#include "SiliconScorers.hh"
-#include "DriftElectronScorers.hh"
+#include "TPCScorers.hh"
#include "RootOutput.h"
#include "MaterialManager.hh"
#include "NPSDetectorFactory.hh"
@@ -79,39 +79,50 @@ namespace Actar_NS{
// Energy and time Resolution
const double ChargeThreshold = 0;
const double ResoTime = 0.1*ns ;
- //const double ResoEnergy = 1.0*MeV ;
- const double ChamberThickness = 55*cm ;
- const double ChamberWidth = 50*cm ;
+ const double ResoCharge = 5./100 ;
+ const double ChamberThickness = 376*mm ;
+ const double ChamberWidth = 376*mm ;
const double ChamberHeight = 40*cm ;
//const int NumberOfPads = 16384;
- const int PadX = 128;
- const int PadZ = 128;
-
+ //const int PadX = 128;
+ //const int PadZ = 128;
+
const double Nose_Rmin = 2.5*cm;
const double Nose_Rmax = 3.5*cm;
const double Nose_Length = 12*cm;
-
+
const double Mylar_Rmax = 3.5*cm;
const double Mylar_Thickness = 7*micrometer;
-
+
const double XGazVolume = 256.*mm;
const double YGazVolume = 256.*mm;
const double ZGazVolume = 256.*mm;
-
+
const double SiliconHeight = 53.*mm;
const double SiliconWidth = 53.*mm;
const double SiliconThickness = 0.7*mm;
const double DistInterSi = 1.*mm;
- const double Si_PosZ=15.*cm;
+ const double Si_PosZ=175.*mm;
const double ResoSilicon = 0.80/2.35;
const double EnergyThreshold = 0.1;
-
+
+ const double VamosSiliconHeight = 70*mm;
+ const double VamosSiliconWidth = 50*mm;
+ const double VamosSiliconThickness = 0.5*mm;
+ const double VamosSiliconDistanInterSi = 1*mm;
+ const double VamosSilicon_PosZ = -160*mm;
+
const double CsIThickness = 1.*cm;
const double CsIHeight = 2.5*cm;
const double CsIWidth = 2.5*cm;
const double DistInterCsI = 1.*mm;
const double CsI_PosZ = 16.*cm;
const double ResoCsI = 0.200/2.35;
+
+ const double BeamDumpRadius = 30*mm;
+ const double BeamDumpThickness = 5*mm;
+ const double BeamDump_PosZ = 160*mm;
+
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -119,20 +130,52 @@ namespace Actar_NS{
// Actar Specific Method
Actar::Actar(){
m_Event = new TActarData() ;
+ m_EventReduced = new MEventReduced();
m_ActarScorer = 0;
m_SquareDetector = 0;
-
+
// RGB Color + Transparency
- m_VisChamber = new G4VisAttributes(G4Colour(0.7, 0.7, 0.7, 0.3));
- m_VisWindows = new G4VisAttributes(G4Colour(1, 0, 0, 0.25));
- m_VisGas = new G4VisAttributes(G4Colour(0, 0.5, 0.5, 0.3));
- m_VisPads = new G4VisAttributes(G4Colour(255, 223, 50, 0.8));
- m_SiliconVisAtt = new G4VisAttributes(G4Colour(0.529412, 0.807843, 0.980392, 0.95)) ;
- m_CsIVisAtt = new G4VisAttributes(G4Colour(0.429412, 0.607843, 0.780392, 0.95));
+ m_VisChamber = new G4VisAttributes(G4Colour(0.7, 0.7, 0.7, 0.3));
+ m_VisWindows = new G4VisAttributes(G4Colour(1, 0, 0, 0.25));
+ m_VisGas = new G4VisAttributes(G4Colour(0, 0.5, 0.5, 0.3));
+ m_VisPads = new G4VisAttributes(G4Colour(255, 223, 50, 0.8));
+ m_VisMicromegas = new G4VisAttributes(G4Colour(100, 100, 100, 0.4));
+ m_SiliconVisAtt = new G4VisAttributes(G4Colour(0.529412, 0.807843, 0.980392, 0.95)) ;
+ m_CsIVisAtt = new G4VisAttributes(G4Colour(0.429412, 0.607843, 0.780392, 0.95));
+ m_BeamDumpVisAtt = new G4VisAttributes(G4Colour(0.9, 0.5, 0.5));
m_VisPads->SetForceWireframe(true);
-
+
+ m_build_BeamDump= 0;
m_build_Silicon=1;
+ m_build_Vamos_Silicon=0;
m_build_CsI=1;
+
+ // Lookup table //
+ bool ReadingLookupTable = false;
+ string LT_FileName = "./Detectors/Actar/LT.dat";
+ ifstream LTConfigFile;
+ LTConfigFile.open(LT_FileName.c_str());
+ if(!LTConfigFile.is_open()){
+ cout << "No Lookup Table in " << LT_FileName << " found!" << endl;
+ return;
+ }
+ else{
+ cout << "/// Using LookupTable from: " << LT_FileName << " ///" << endl;
+ int co, as, ag, ch;
+ int pX, pY;
+ for(int i=0;i<NumberOfCobo*NumberOfASAD*NumberOfAGET*NumberOfChannel;i++){
+ LTConfigFile >> co >> as >> ag >> ch >> pX >> pY;
+ if(pX!=-1 && pY !=-1){
+ m_PadToCobo[pX][pY] = co;
+ m_PadToAsad[pX][pY] = as;
+ m_PadToAGET[pX][pY] = ag;
+ m_PadToChannel[pX][pY] = ch;
+ }
+ }
+ ReadingLookupTable = true;
+ }
+ LTConfigFile.close();
+
}
Actar::~Actar(){
@@ -157,47 +200,51 @@ void Actar::AddDetector(double R, double Theta, double Phi, string Shape){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4LogicalVolume* Actar::BuildDetector(){
-
- G4Material* Cu= MaterialManager::getInstance()->GetMaterialFromLibrary("Cu");
- G4Material* Si= MaterialManager::getInstance()->GetMaterialFromLibrary("Si");
- G4Material* Al= MaterialManager::getInstance()->GetMaterialFromLibrary("Al");
- G4Material* Mylar= MaterialManager::getInstance()->GetMaterialFromLibrary("Mylar");
- G4Material* MaterialCsI = MaterialManager::getInstance()->GetMaterialFromLibrary("CsI");
-
- if(!m_SquareDetector){
+
+ G4Material* Cu= MaterialManager::getInstance()->GetMaterialFromLibrary("Cu");
+ G4Material* Si= MaterialManager::getInstance()->GetMaterialFromLibrary("Si");
+ G4Material* Al= MaterialManager::getInstance()->GetMaterialFromLibrary("Al");
+ G4Material* Mylar= MaterialManager::getInstance()->GetMaterialFromLibrary("Mylar");
+ G4Material* MaterialCsI = MaterialManager::getInstance()->GetMaterialFromLibrary("CsI");
+
+ if(!m_SquareDetector){
// Main volume
G4Box* sChamber = new G4Box("Actar_Box",Actar_NS::ChamberWidth*0.5,
Actar_NS::ChamberHeight*0.5,Actar_NS::ChamberThickness*0.5);
-
+
//Nose volume
G4Tubs* sNose = new G4Tubs("Actar_Nose",Actar_NS::Nose_Rmin, Actar_NS::Nose_Rmax,Actar_NS::Nose_Length*0.5,
- 0*deg, 360*deg);
-
+ 0*deg, 360*deg);
+
//Mylar volume
G4Tubs* sWindows = new G4Tubs("Actar_Windows",0, Actar_NS::Mylar_Rmax,Actar_NS::Mylar_Thickness*0.5,
- 0*deg, 360*deg);
-
+ 0*deg, 360*deg);
+
// Cage volume
G4Box* sCage = new G4Box("Actar_Gas",Actar_NS::XGazVolume*0.5,
- Actar_NS::YGazVolume*0.5,Actar_NS::ZGazVolume*0.5);
-
+ Actar_NS::YGazVolume*0.5,Actar_NS::ZGazVolume*0.5);
+
// Pad
- G4Box* sPad = new G4Box("Actar_Pad",2*mm*0.5,
- 1*um*0.5,2*mm*0.5);
-
+ G4Box* sPad = new G4Box("Actar_Pad",256*mm*0.5,
+ 2*mm*0.5,256*mm*0.5);
+
+ // Micromegas
+ G4Box* sMicromegas = new G4Box("Actar_Micromegas",256*mm*0.5,
+ 220*micrometer*0.5,256*mm*0.5);
+
// Cathode
- //G4Box* sCathode = new G4Box("Actar_Cathode",26.5*cm*0.5,
- // 1*um*0.5,25.6*cm*0.5);
-
-
-
+ G4Box* sCathode = new G4Box("Actar_Cathode",26.5*cm*0.5,
+ 1*um*0.5,25.6*cm*0.5);
+
+
+
unsigned const int NumberOfGasMix = m_GasMaterial.size();
-
+
double density=0;
double density_sum=0;
vector<G4Material*> GasComponent;
vector<double> FractionMass;
-
+
for(unsigned int i=0; i<NumberOfGasMix; i++){
GasComponent.push_back(MaterialManager::getInstance()->GetGasFromLibrary(m_GasMaterial[i],m_Pressure,m_Temperature) );
}
@@ -206,188 +253,252 @@ G4LogicalVolume* Actar::BuildDetector(){
density_sum += GasComponent[i]->GetDensity();
}
//cout << "density = " << density*cm3/g << endl;
-
+
for(unsigned int i=0; i<NumberOfGasMix; i++){
FractionMass.push_back(GasComponent[i]->GetDensity()/density_sum);
}
-
+
G4Material* GasMaterial = new G4Material("GasMix", density, NumberOfGasMix, kStateGas, m_Temperature, m_Pressure);
- G4Material* DriftGasMaterial = new G4Material("GasMix", density, NumberOfGasMix, kStateGas, m_Temperature, m_Pressure);
-
+ G4Material* DriftGasMaterial = new G4Material("DriftGasMix", density, NumberOfGasMix, kStateGas, m_Temperature, m_Pressure);
+
for(unsigned int i=0; i<NumberOfGasMix; i++){
GasMaterial->AddMaterial(GasComponent[i], FractionMass[i]);
DriftGasMaterial->AddMaterial(GasComponent[i], FractionMass[i]);
+ cout << GasComponent[i] << endl;
}
-
+
G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable();
MPT->AddConstProperty("DE_PAIRENERGY",20*eV);
- MPT->AddConstProperty("DE_YIELD",1e-2);
- //MPT->AddConstProperty("DE_YIELD",1e-1);
- //MPT->AddConstProperty("DE_AMPLIFICATION",5);
+ MPT->AddConstProperty("DE_YIELD",3e-1);
+ //MPT->AddConstProperty("DE_AMPLIFICATION",2);
MPT->AddConstProperty("DE_ABSLENGTH",1*pc);
- MPT->AddConstProperty("DE_DRIFTSPEED",5.*cm/microsecond);
- //MPT->AddConstProperty("DE_TRANSVERSALSPREAD",5e-5*mm2/ns);
- //MPT->AddConstProperty("DE_LONGITUDINALSPREAD",5e-5*mm2/ns);
- MPT->AddConstProperty("DE_TRANSVERSALSPREAD",7e-6*mm2/ns);
- MPT->AddConstProperty("DE_LONGITUDINALSPREAD",7e-6*mm2/ns);
-
+ MPT->AddConstProperty("DE_DRIFTSPEED",0.8*cm/microsecond);
+ MPT->AddConstProperty("DE_TRANSVERSALSPREAD",1e-5*mm2/ns);
+ MPT->AddConstProperty("DE_LONGITUDINALSPREAD",7e-5*mm2/ns);
+
DriftGasMaterial->SetMaterialPropertiesTable(MPT);
-
+
G4MaterialPropertiesTable* MPT2 = new G4MaterialPropertiesTable();
- MPT2->AddConstProperty("DE_YIELD",1);
- MPT2->AddConstProperty("DE_AMPLIFICATION",2);
+ MPT2->AddConstProperty("DE_AMPLIFICATION",1000);
MPT2->AddConstProperty("DE_ABSLENGTH",1*pc);
-
- Al->SetMaterialPropertiesTable(MPT2);
-
- m_SquareDetector = new G4LogicalVolume(sChamber,GasMaterial,"logic_Actar_Box",0,0,0);
- G4LogicalVolume* logicGas = new G4LogicalVolume(sCage,DriftGasMaterial,"logic_Gas",0,0,0);
- G4LogicalVolume* logicNose = new G4LogicalVolume(sNose,Al,"logic_Nose",0,0,0);
+
+ //Al->SetMaterialPropertiesTable(MPT2);
+
+ m_SquareDetector = new G4LogicalVolume(sChamber,GasMaterial,"logic_Actar_Box",0,0,0);
+ m_logicGas = new G4LogicalVolume(sCage,DriftGasMaterial,"logic_Gas",0,0,0);
G4LogicalVolume* logicPad = new G4LogicalVolume(sPad,Cu,"logic_Pad",0,0,0);
- //G4LogicalVolume* logicCathode = new G4LogicalVolume(sCathode,Cu,"logic_Cathode",0,0,0);
+ G4LogicalVolume* logicMicromegas = new G4LogicalVolume(sMicromegas,Al,"logic_Micromegas",0,0,0);
+
+ G4LogicalVolume* logicNose = new G4LogicalVolume(sNose,Cu,"logic_Nose",0,0,0);
+ G4LogicalVolume* logicCathode = new G4LogicalVolume(sCathode,Cu,"logic_Cathode",0,0,0);
G4LogicalVolume* logicWindows = new G4LogicalVolume(sWindows,Mylar,"logic_Windows",0,0,0);
-
+
G4RotationMatrix* Rot = new G4RotationMatrix();
- new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,0,-Actar_NS::ChamberThickness*0.5+Actar_NS::Nose_Length*0.5)),
- logicNose,
- "ActarNose",m_SquareDetector,false,0);
-
- new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,0,-Actar_NS::ChamberThickness*0.5+Actar_NS::Mylar_Thickness*0.5+Actar_NS::Nose_Length)),
- logicWindows,
- "ActarEntranceWindows",m_SquareDetector,false,0);
-
+ //new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,0,-Actar_NS::ChamberThickness*0.5+Actar_NS::Nose_Length*0.5)),
+ // logicNose,
+ // "ActarNose",m_SquareDetector,false,0);
+
+ //new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,0,-Actar_NS::ChamberThickness*0.5+Actar_NS::Mylar_Thickness*0.5+Actar_NS::Nose_Length)),
+ // logicWindows,
+ // "ActarEntranceWindows",m_SquareDetector,false,0);
+
new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,0,0)),
- logicGas,
+ m_logicGas,
"ActarGas",m_SquareDetector,false,0);
-
- int pad=0;
- m_PadToXRow.clear();
- m_PadToZColumn.clear();
- for(int i=0; i<Actar_NS::PadX; i++){
- for(int j=0; j<Actar_NS::PadZ; j++){
- m_PadToXRow[pad] = i;
- m_PadToZColumn[pad] = j;
- double X=(i-64)*2*mm;
- double Z=(j-64)*2*mm;
- new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(X,Actar_NS::YGazVolume*0.5,Z)),
- logicPad,
- "ActarPad",logicGas,false,pad+1);
-
- pad++;
- }
- }
-
+
+ //new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,Actar_NS::YGazVolume*0.5-0.3*cm,0)),
+ // logicMicromegas,
+ // "ActarMicromegas",m_logicGas,false,0);
+
+ new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,Actar_NS::YGazVolume*0.5,0)),
+ logicPad,
+ "ActarPad",m_logicGas,false,0);
+
+ /*int pad=0;
+ m_PadToXRow.clear();
+ m_PadToZColumn.clear();
+ for(int i=0; i<Actar_NS::PadX; i++){
+ for(int j=0; j<Actar_NS::PadZ; j++){
+ m_PadToXRow[pad] = i;
+ m_PadToZColumn[pad] = j;
+ double X=(i-64)*2*mm;
+ double Z=(j-64)*2*mm;
+ new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(X,Actar_NS::YGazVolume*0.5,Z)),
+ logicPad,
+ "ActarPad",m_logicGas,false,pad+1);
+
+ pad++;
+ }
+ }*/
+
/*new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(3*cm-0.5*1*um,0,0)),
- logicCathode,
- "ActarCathode",logicGas,false,0);
-
-
-
- new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,0,-6*cm+6*micrometer)),
- logicWindows,
- "ActarEntranceWindows",m_SquareDetector,false,0);*/
-
- G4ElectricField* field = new G4UniformElectricField(G4ThreeVector(0.0,-100*volt/cm,0.0));
+ logicCathode,
+ "ActarCathode",m_logicGas,false,0);
+
+
+
+ new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,0,-6*cm+6*micrometer)),
+ logicWindows,
+ "ActarEntranceWindows",m_SquareDetector,false,0);*/
+
+ G4ElectricField* field = new G4UniformElectricField(G4ThreeVector(0.0,-400*volt/cm,0.0));
// Create an equation of motion for this field
G4EqMagElectricField* Equation = new G4EqMagElectricField(field);
G4MagIntegratorStepper* Stepper = new G4ClassicalRK4( Equation, 8 );
-
+
// Get the global field manager
G4FieldManager* FieldManager= new G4FieldManager();
// Set this field to the global field manager
FieldManager->SetDetectorField(field);
- logicGas->SetFieldManager(FieldManager,true);
-
+ m_logicGas->SetFieldManager(FieldManager,true);
+
G4MagInt_Driver* IntgrDriver = new G4MagInt_Driver(0.1*mm,
Stepper,
Stepper->GetNumberOfVariables() );
-
+
G4ChordFinder* ChordFinder = new G4ChordFinder(IntgrDriver);
FieldManager->SetChordFinder( ChordFinder );
-
-
+
+
logicPad->SetSensitiveDetector(m_ActarScorer);
logicNose->SetVisAttributes(m_VisChamber);
m_SquareDetector->SetVisAttributes(m_VisChamber);
- logicGas->SetVisAttributes(m_VisGas);
+ m_logicGas->SetVisAttributes(m_VisGas);
logicWindows->SetVisAttributes(m_VisWindows);
logicPad->SetVisAttributes(m_VisPads);
+ logicMicromegas->SetVisAttributes(m_VisMicromegas);
//m_SquareDetector->SetSensitiveDetector(m_ActarScorer);
}
-
-
+
+
+ ////////////////////////////////////////////////////
+ ///////////////////// Beam Dump ////////////////////
+ ////////////////////////////////////////////////////
+ if(m_build_BeamDump){
+ G4Tubs* sBeamDump = new G4Tubs("Actar_BeamDump",0*mm, Actar_NS::BeamDumpRadius,Actar_NS::BeamDumpThickness*0.5,
+ 0*deg, 360*deg);
+ m_LogicBeamDump = new G4LogicalVolume(sBeamDump, Cu, "logicBeamDump",0,0,0);
+
+ G4ThreeVector positionBeamDump = G4ThreeVector(0, 0, Actar_NS::BeamDump_PosZ);
+ new G4PVPlacement(new G4RotationMatrix(0,0,0),
+ positionBeamDump,
+ m_LogicBeamDump,"BeamDump",
+ m_SquareDetector,false,0);
+
+ m_LogicBeamDump->SetVisAttributes(m_BeamDumpVisAtt) ;
+
+ }
///////////////////////////////////////////////////
///////////////////// Thin Si /////////////////////
///////////////////////////////////////////////////
+ int SiliconNumber=0;
if(m_build_Silicon){
G4Box* solidSi = new G4Box("Si", 0.5*Actar_NS::SiliconWidth, 0.5*Actar_NS::SiliconHeight, 0.5*Actar_NS::SiliconThickness); ;
m_LogicSilicon = new G4LogicalVolume(solidSi, Si, "logicSi", 0, 0, 0);
-
- int SiliconNumber=0;
+
for(int k=0;k<4; k++){
- for(int p=0; p<5; p++){
- double PosX;
- double PosY;
- if(k==0) PosY= -1.5*Actar_NS::SiliconHeight-2*Actar_NS::DistInterSi;
- if(k==1) PosY= -0.5*Actar_NS::SiliconHeight-1*Actar_NS::DistInterSi;
- if(k==2) PosY= 0.5*Actar_NS::SiliconHeight+1*Actar_NS::DistInterSi;
- if(k==3) PosY= 1.5*Actar_NS::SiliconHeight+2*Actar_NS::DistInterSi;
- if(p==0) PosX= -2*Actar_NS::SiliconWidth-2*Actar_NS::DistInterSi;
- if(p==1) PosX= -1*Actar_NS::SiliconWidth-1*Actar_NS::DistInterSi;
- if(p==2) PosX= 0;
- if(p==3) PosX= 1*Actar_NS::SiliconWidth+2*Actar_NS::DistInterSi;
- if(p==4) PosX= 2*Actar_NS::SiliconWidth+1*Actar_NS::DistInterSi;
-
- G4ThreeVector positionSi = G4ThreeVector(PosX, PosY, Actar_NS::Si_PosZ);
- new G4PVPlacement(new G4RotationMatrix(0,0,0),
- positionSi,
- m_LogicSilicon,"Si",
- m_SquareDetector,false,SiliconNumber);
- SiliconNumber++;
+ for(int p=0; p<5; p++){
+ double PosX;
+ double PosY;
+ if(k==0) PosY= -1.5*Actar_NS::SiliconHeight-2*Actar_NS::DistInterSi;
+ if(k==1) PosY= -0.5*Actar_NS::SiliconHeight-1*Actar_NS::DistInterSi;
+ if(k==2) PosY= 0.5*Actar_NS::SiliconHeight+1*Actar_NS::DistInterSi;
+ if(k==3) PosY= 1.5*Actar_NS::SiliconHeight+2*Actar_NS::DistInterSi;
+ if(p==0) PosX= -2*Actar_NS::SiliconWidth-2*Actar_NS::DistInterSi;
+ if(p==1) PosX= -1*Actar_NS::SiliconWidth-1*Actar_NS::DistInterSi;
+ if(p==2) PosX= 0;
+ if(p==3) PosX= 1*Actar_NS::SiliconWidth+2*Actar_NS::DistInterSi;
+ if(p==4) PosX= 2*Actar_NS::SiliconWidth+1*Actar_NS::DistInterSi;
+
+ G4ThreeVector positionSi = G4ThreeVector(PosX, PosY, Actar_NS::Si_PosZ);
+ new G4PVPlacement(new G4RotationMatrix(0,0,0),
+ positionSi,
+ m_LogicSilicon,"Si",
+ m_SquareDetector,false,SiliconNumber);
+ SiliconNumber++;
}
- }
-
+ }
+
// Set Si sensible
m_LogicSilicon->SetSensitiveDetector(m_SiliconScorer);
-
+
// Visualisation of ThinSi
- m_LogicSilicon->SetVisAttributes(m_SiliconVisAtt) ;
+ m_LogicSilicon->SetVisAttributes(m_SiliconVisAtt);
}
-
+
+ ////////////////////////////////////////////////////
+ ///////////////////// Vamos Si /////////////////////
+ ////////////////////////////////////////////////////
+ if(m_build_Vamos_Silicon){
+ G4Box* solidSi = new G4Box("Si", 0.5*Actar_NS::VamosSiliconWidth, 0.5*Actar_NS::VamosSiliconHeight, 0.5*Actar_NS::VamosSiliconThickness); ;
+ m_LogicVamosSilicon = new G4LogicalVolume(solidSi, Si, "logicVamosSi", 0, 0, 0);
+
+
+ int VamosSiliconNumber=0;
+ for(int k=0;k<4; k++){
+ for(int p=0; p<3; p++){
+ double PosX;
+ double PosY;
+ if(k==0) PosX= -35*mm -0.5*Actar_NS::VamosSiliconWidth - Actar_NS::VamosSiliconDistanInterSi;
+ if(k==1) PosX= -35*mm -1.5*Actar_NS::VamosSiliconWidth - Actar_NS::VamosSiliconDistanInterSi;
+ if(k==2) PosX= 35*mm + 0.5*Actar_NS::VamosSiliconWidth + Actar_NS::VamosSiliconDistanInterSi;
+ if(k==3) PosX= 35*mm + 1.5*Actar_NS::VamosSiliconWidth + Actar_NS::VamosSiliconDistanInterSi;
+
+ if(p==0) PosY= -Actar_NS::VamosSiliconHeight-Actar_NS::VamosSiliconDistanInterSi;
+ if(p==1) PosY= 0;
+ if(p==2) PosY= Actar_NS::VamosSiliconHeight+Actar_NS::VamosSiliconDistanInterSi;
+
+ G4ThreeVector positionSi = G4ThreeVector(PosX, PosY, Actar_NS::VamosSilicon_PosZ);
+ new G4PVPlacement(new G4RotationMatrix(0,0,0),
+ positionSi,
+ m_LogicVamosSilicon,"Si",
+ m_SquareDetector,false,SiliconNumber);
+ VamosSiliconNumber++;
+ SiliconNumber++;
+ }
+ }
+
+ // Set Si sensible
+ m_LogicVamosSilicon->SetSensitiveDetector(m_SiliconScorer);
+
+ // Visualisation of ThinSi
+ m_LogicVamosSilicon->SetVisAttributes(m_SiliconVisAtt);
+ }
+
///////////////////////////////////////////////
///////////////////// CsI /////////////////////
///////////////////////////////////////////////
if(m_build_CsI){
G4Box* solidCsI = new G4Box("Si", 0.5*Actar_NS::CsIWidth, 0.5*Actar_NS::CsIHeight, 0.5*Actar_NS::CsIThickness); ;
m_LogicCsICrystal = new G4LogicalVolume(solidCsI, MaterialCsI, "logicCsI", 0, 0, 0);
-
+
int CsINumber=0;
for(int k=0;k<8; k++){
- for(int p=0; p<10; p++){
- double PosX;
- double PosY;
- if(k<4) PosY= -0.5*Actar_NS::CsIHeight-0.5*Actar_NS::DistInterCsI+(k-3)*(Actar_NS::CsIHeight+Actar_NS::DistInterCsI);
- if(k>3) PosY= 0.5*Actar_NS::CsIHeight+0.5*Actar_NS::DistInterCsI+(k-4)*(Actar_NS::CsIHeight+Actar_NS::DistInterCsI);
-
- if(p<5) PosX= 0.5*Actar_NS::CsIWidth+0.5*Actar_NS::DistInterCsI+(4-p)*(Actar_NS::CsIWidth+Actar_NS::DistInterCsI);
- if(p>4) PosX= -0.5*Actar_NS::CsIWidth-0.5*Actar_NS::DistInterCsI+(5-p)*(Actar_NS::CsIWidth+Actar_NS::DistInterCsI);
-
- G4ThreeVector positionCsI = G4ThreeVector(PosX, PosY, Actar_NS::CsI_PosZ);
- new G4PVPlacement(new G4RotationMatrix(0,0,0),
- positionCsI,
- m_LogicCsICrystal,"CsI",
- m_SquareDetector,false,CsINumber);
- CsINumber++;
+ for(int p=0; p<10; p++){
+ double PosX;
+ double PosY;
+ if(k<4) PosY= -0.5*Actar_NS::CsIHeight-0.5*Actar_NS::DistInterCsI+(k-3)*(Actar_NS::CsIHeight+Actar_NS::DistInterCsI);
+ if(k>3) PosY= 0.5*Actar_NS::CsIHeight+0.5*Actar_NS::DistInterCsI+(k-4)*(Actar_NS::CsIHeight+Actar_NS::DistInterCsI);
+
+ if(p<5) PosX= 0.5*Actar_NS::CsIWidth+0.5*Actar_NS::DistInterCsI+(4-p)*(Actar_NS::CsIWidth+Actar_NS::DistInterCsI);
+ if(p>4) PosX= -0.5*Actar_NS::CsIWidth-0.5*Actar_NS::DistInterCsI+(5-p)*(Actar_NS::CsIWidth+Actar_NS::DistInterCsI);
+
+ G4ThreeVector positionCsI = G4ThreeVector(PosX, PosY, Actar_NS::CsI_PosZ);
+ new G4PVPlacement(new G4RotationMatrix(0,0,0),
+ positionCsI,
+ m_LogicCsICrystal,"CsI",
+ m_SquareDetector,false,CsINumber);
+ CsINumber++;
}
- }
-
+ }
+
// Set Si sensible
m_LogicCsICrystal->SetSensitiveDetector(m_CsIScorer);
-
+
// Visualisation of ThinSi
m_LogicCsICrystal->SetVisAttributes(m_CsIVisAtt) ;
}
-
+
return m_SquareDetector;
}
@@ -403,15 +514,15 @@ void Actar::ReadConfiguration(NPL::InputParser parser){
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Actar");
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << "//// " << blocks.size() << " detectors found " << endl;
-
- vector<string> cart = {"POS","Shape","GasMaterial","GasFraction","Temperature","Pressure","Si","CsI"};
- vector<string> sphe = {"R","Theta","Phi","Shape","GasMaterial","GasFraction","Temperature","Pressure","Si","CsI"};
-
+
+ vector<string> cart = {"POS","Shape","GasMaterial","GasFraction","Temperature","Pressure","Si","VamosSi","CsI","BeamDump"};
+ //vector<string> sphe = {"R","Theta","Phi","Shape","GasMaterial","GasFraction","Temperature","Pressure","Si","CsI","BeamDump"};
+
for(unsigned int i = 0 ; i < blocks.size() ; i++){
if(blocks[i]->HasTokenList(cart)){
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << endl << "//// Actar " << i+1 << endl;
-
+
G4ThreeVector Pos = NPS::ConvertVector(blocks[i]->GetTVector3("POS","mm"));
string Shape = blocks[i]->GetString("Shape");
vector<string> GasName = blocks[i]->GetVectorString("GasMaterial");
@@ -423,11 +534,13 @@ void Actar::ReadConfiguration(NPL::InputParser parser){
m_Temperature = blocks[i]->GetDouble("Temperature","kelvin");
m_Pressure = blocks[i]->GetDouble("Pressure","bar");
m_build_Silicon = blocks[i]->GetInt("Si");
+ m_build_Vamos_Silicon = blocks[i]->GetInt("VamosSi");
m_build_CsI = blocks[i]->GetInt("CsI");
-
+ m_build_BeamDump = blocks[i]->GetInt("BeamDump");
+
AddDetector(Pos,Shape);
}
- else if(blocks[i]->HasTokenList(sphe)){
+ /*else if(blocks[i]->HasTokenList(sphe)){
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << endl << "//// Actar " << i+1 << endl;
double R = blocks[i]->GetDouble("R","mm");
@@ -444,9 +557,9 @@ void Actar::ReadConfiguration(NPL::InputParser parser){
m_Pressure = blocks[i]->GetDouble("Pressure","bar");
m_build_Silicon = blocks[i]->GetInt("Si");
m_build_CsI = blocks[i]->GetInt("CsI");
-
+
AddDetector(R,Theta,Phi,Shape);
- }
+ }*/
else{
cout << "ERROR: check your input file formatting " << endl;
exit(1);
@@ -461,7 +574,7 @@ void Actar::ReadConfiguration(NPL::InputParser parser){
// Called After DetecorConstruction::AddDetector Method
void Actar::ConstructDetector(G4LogicalVolume* world){
for (unsigned short i = 0 ; i < m_R.size() ; i++) {
-
+
G4double wX = m_R[i] * sin(m_Theta[i] ) * cos(m_Phi[i] ) ;
G4double wY = m_R[i] * sin(m_Theta[i] ) * sin(m_Phi[i] ) ;
G4double wZ = m_R[i] * cos(m_Theta[i] ) ;
@@ -478,13 +591,42 @@ void Actar::ConstructDetector(G4LogicalVolume* world){
G4ThreeVector v = w.cross(u);
v = v.unit();
u = u.unit();
-
+
G4RotationMatrix* Rot = new G4RotationMatrix(u,v,w);
-
+
new G4PVPlacement(G4Transform3D(*Rot,Det_pos),
- BuildDetector(),
- "Actar",world,false,i+1);
+ BuildDetector(),
+ "Actar",world,false,i+1);
+ }
+
+ if(!m_ReactionRegion){
+ G4ProductionCuts* ecut = new G4ProductionCuts();
+ ecut->SetProductionCut(1000,"e-");
+
+ m_ReactionRegion= new G4Region("NPSimulationProcess");
+ m_ReactionRegion->SetProductionCuts(ecut);
+ m_ReactionRegion->AddRootLogicalVolume(m_logicGas);
+ m_ReactionRegion->SetUserLimits(new G4UserLimits(1.2*mm));
+
+ G4Region* Region_cut = new G4Region("RegionCut");
+ Region_cut->SetProductionCuts(ecut);
+ Region_cut->AddRootLogicalVolume(m_SquareDetector);
+ }
+
+ G4FastSimulationManager* mng = m_ReactionRegion->GetFastSimulationManager();
+
+ unsigned int size = m_ReactionModel.size();
+ for(unsigned int i = 0 ; i < size ; i++){
+ mng->RemoveFastSimulationModel(m_ReactionModel[i]);
}
+
+ m_ReactionModel.clear();
+ G4VFastSimulationModel* fsm;
+ fsm = new NPS::BeamReaction("BeamReaction",m_ReactionRegion);
+ m_ReactionModel.push_back(fsm);
+ fsm = new NPS::Decay("Decay",m_ReactionRegion);
+ m_ReactionModel.push_back(fsm);
+
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Add Detector branch to the EventTree.
@@ -492,29 +634,32 @@ void Actar::ConstructDetector(G4LogicalVolume* world){
void Actar::InitializeRootOutput(){
RootOutput *pAnalysis = RootOutput::getInstance();
TTree *pTree = pAnalysis->GetTree();
- if(!pTree->FindBranch("Actar")){
- pTree->Branch("Actar", "TActarData", &m_Event) ;
+ if(!pTree->FindBranch("data")){
+ //pTree->Branch("Actar", "TActarData", &m_Event) ;
+ pTree->Branch("data", "MEventReduced", &m_EventReduced) ;
}
- pTree->SetBranchAddress("Actar", &m_Event) ;
+ //pTree->SetBranchAddress("Actar", &m_Event) ;
+ pTree->SetBranchAddress("data", &m_EventReduced) ;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Read sensitive part and fill the Root tree.
// Called at in the EventAction::EndOfEventAvtion
void Actar::ReadSensitive(const G4Event* event){
- m_Event->Clear();
-
+ m_EventReduced->CoboAsad.clear();
+
+
//////////////
// Pad scorer
NPS::HitsMap<G4double*>* PadHitMap;
std::map<G4int, G4double**>::iterator Pad_itr;
-
+
G4int PadCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("ActarScorer/Actar_dig");
PadHitMap = (NPS::HitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(PadCollectionID));
-
// Loop on the Pad map
- TH1D* h = new TH1D("h","h",25000,0,25000);
+ //TH1D* h = new TH1D("h","h",25000,0,25000);
for (Pad_itr = PadHitMap->GetMap()->begin() ; Pad_itr != PadHitMap->GetMap()->end() ; Pad_itr++){
+ ReducedData DataReduced;
G4double* Info = *(Pad_itr->second);
// Interraction Coordinates
ms_InterCoord->SetDetectedPositionX(Info[2]) ;
@@ -522,145 +667,167 @@ void Actar::ReadSensitive(const G4Event* event){
ms_InterCoord->SetDetectedPositionZ(Info[4]) ;
ms_InterCoord->SetDetectedAngleTheta(Info[5]/deg) ;
ms_InterCoord->SetDetectedAnglePhi(Info[6]/deg) ;
-
- double Count = Info[0];
- double Time = RandGauss::shoot(Info[1],Actar_NS::ResoTime);
+
+ double Count = RandGauss::shoot(Info[0],Actar_NS::ResoCharge*Info[0]);
+ double Time = Info[1];//RandGauss::shoot(Info[1],Actar_NS::ResoTime);
//int iTime = ((int) Time*20/512)+1;
int PadNbr = Info[7];
+ int Pad_X = Info[8];//m_PadToXRow[PadNbr];
+ int Pad_Y = Info[9];//m_PadToZColumn[PadNbr];
+
+ int co = m_PadToCobo[Pad_X][Pad_Y];
+ int as = m_PadToAsad[Pad_X][Pad_Y];
+ int ag = m_PadToAGET[Pad_X][Pad_Y];
+ int ch = m_PadToChannel[Pad_X][Pad_Y];
+
if(Count>Actar_NS::ChargeThreshold){
- m_Event->SetTime(Time);
- m_Event->SetPadNumber(PadNbr);
- m_Event->SetCharge(Count);
- m_Event->SetRowNumber(m_PadToXRow[PadNbr]);
- m_Event->SetColumnNumber(m_PadToZColumn[PadNbr]);
- }
-
- if(Count){
- h->Fill(Info[1],Info[0]);
+ DataReduced.globalchannelid = ch+(ag<<7)+(as<<9)+(co<<11);
+ DataReduced.peakheight.push_back(Count);
+ DataReduced.peaktime.push_back(Time);
}
+ m_EventReduced->CoboAsad.push_back(DataReduced);
+ /*if(Count){
+ h->Fill(Info[1],Info[0]);
+ }*/
}
-
- vector<double> Q, T;
- for(int i=0; i<h->GetNbinsX(); i++){
- double count = h->GetBinContent(i);
- double time = h->GetBinCenter(i);
- if(count){
- Q.push_back(count);
- T.push_back(time+500);
-
- }
- }
- // clear map for next event
- //SimulateDigitizer(Q,T,1.40*microsecond,0,8750,25,5);
- delete h;
+
+ /*vector<double> Q, T;
+ for(int i=0; i<h->GetNbinsX(); i++){
+ double count = h->GetBinContent(i);
+ double time = h->GetBinCenter(i);
+ if(count){
+ Q.push_back(count);
+ T.push_back(time+500);
+
+ }
+ }
+ // clear map for next event
+ SimulateDigitizer(Q,T,1.40*microsecond,0,8750,25,5);
+ delete h;*/
+
PadHitMap->clear();
-
+
// Silicon //
if(m_build_Silicon){
+ ReducedData DataReduced;
+
NPS::HitsMap<G4double*>* SiHitMap;
std::map<G4int, G4double**>::iterator Si_itr;
-
+
G4int SiCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("SiliconScorer/SiliconScorer");
SiHitMap = (NPS::HitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(SiCollectionID));
-
+
// Loop on the ThinSi map
for (Si_itr = SiHitMap->GetMap()->begin() ; Si_itr != SiHitMap->GetMap()->end() ; Si_itr++){
G4double* Info = *(Si_itr->second);
double E_Si = RandGauss::shoot(Info[0],Actar_NS::ResoSilicon);
-
+
+ int co=31;
+ int as=0;
+ int ag=0;
+ int ch=Info[7];
+
if(E_Si>Actar_NS::EnergyThreshold){
- m_Event->SetSiliconEnergy(E_Si);
- m_Event->SetSiliconDetectorNumber(Info[7]);
- m_Event->SetSiliconTime(Info[1]);
+ DataReduced.globalchannelid = ch+(ag<<7)+(as<<9)+(co<<11);
+ DataReduced.peaktime.push_back(ch);
+ DataReduced.peakheight.push_back(E_Si);
}
+ m_EventReduced->CoboAsad.push_back(DataReduced);
}
-
+
// Clear Map for next event
SiHitMap->clear();
}
-
+
// CsI //
if(m_build_CsI){
NPS::HitsMap<G4double*>* CsIHitMap;
std::map<G4int, G4double**>::iterator CsI_itr;
-
+
G4int CsICollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("CsIScorer/CsI");
CsIHitMap = (NPS::HitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(CsICollectionID));
-
+
// Loop on the CsI map
for (CsI_itr = CsIHitMap->GetMap()->begin() ; CsI_itr !=CsIHitMap->GetMap()->end() ; CsI_itr++){
G4double* Info = *(CsI_itr->second);
double E_CsI = RandGauss::shoot(Info[0],Actar_NS::ResoCsI);
-
- if(E_CsI>Actar_NS::EnergyThreshold){
- m_Event->SetCsIEnergy(E_CsI);
- m_Event->SetCsICrystalNumber(Info[2]);
- }
+
+ /*if(E_CsI>Actar_NS::EnergyThreshold){
+ m_Event->SetCsIEnergy(E_CsI);
+ m_Event->SetCsICrystalNumber(Info[2]);
+ }
+ m_EventReduced->CoboAsad.push_back(DataReduced);*/
}
// Clear Map for next event
CsIHitMap->clear();
}
+
+
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void Actar::SimulateDigitizer(vector<double> E, vector<double> T, double fallTime,double start,double stop, double step,double noise){
-
- static string formula;
- formula= "";
- static string Es,Ts,var,cond;
- static string fall;
- fall=std::to_string(fallTime);
-
- for(unsigned int i = 0 ; i < E.size() ; i++){
- if(E[i]!=0 && T[i]!=0){
- Es = std::to_string(E[i]);
- Ts = std::to_string(T[i]);
- cond = ")*(x>"+Ts+")+";
- var = "(x-"+Ts+")";
- formula += Es+"*-1*exp(-"+var+"/"+fall+cond;
- }
- }
- formula+="0";
- //cout << formula << endl;
- TF1* f = new TF1("f",formula.c_str(),start,stop);
- unsigned int size = (stop-start)/step;
- for(unsigned int i = 0 ; i < size ; i++){
- double time = start+i*step;
- double energy = f->Eval(time)+noise*(1-2*G4UniformRand());
- m_Event->AddEnergyPoint(energy,time);
- }
-
- delete f;
-}
+/*void Actar::SimulateDigitizer(vector<double> E, vector<double> T, double fallTime,double start,double stop, double step,double noise){
+
+ static string formula;
+ formula= "";
+ static string Es,Ts,var,cond;
+ static string fall;
+ fall=std::to_string(fallTime);
+
+ for(unsigned int i = 0 ; i < E.size() ; i++){
+ if(E[i]!=0 && T[i]!=0){
+ Es = std::to_string(E[i]);
+ Ts = std::to_string(T[i]);
+ cond = ")*(x>"+Ts+")+";
+ var = "(x-"+Ts+")";
+ formula += Es+"*-1*exp(-"+var+"/"+fall+cond;
+ }
+ }
+ formula+="0";
+ //cout << formula << endl;
+ TF1* f = new TF1("f",formula.c_str(),start,stop);
+ unsigned int size = (stop-start)/step;
+ for(unsigned int i = 0 ; i < size ; i++){
+ double time = start+i*step;
+ double energy = f->Eval(time)+noise*(1-2*G4UniformRand());
+ }
+
+ delete f;
+ }*/
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
////////////////////////////////////////////////////////////////
void Actar::InitializeScorers() {
// This check is necessary in case the geometry is reloaded
-
+
bool already_exist = false;
vector<G4int> NestingLevel;
NestingLevel.push_back(0);
NestingLevel.push_back(2);
-
+
m_ActarScorer = CheckScorer("ActarScorer",already_exist) ;
m_SiliconScorer = CheckScorer("SiliconScorer",already_exist);
+ //m_VamosSiliconScorer = CheckScorer("VamosSiliconScorer",already_exist);
m_CsIScorer = CheckScorer("CsIScorer",already_exist);
-
+
if(already_exist) return;
-
- G4VPrimitiveScorer *SiScorer = new SILICONSCORERS::PS_Silicon_Rectangle("SiliconScorer",0,Actar_NS::SiliconHeight,Actar_NS::SiliconWidth,1,1);
+
+ G4VPrimitiveScorer* SiScorer = new SILICONSCORERS::PS_Silicon_Rectangle("SiliconScorer",0,Actar_NS::SiliconHeight,Actar_NS::SiliconWidth,1,1);
m_SiliconScorer->RegisterPrimitive(SiScorer);
-
+
+ /*G4VPrimitiveScorer* VamosSiScorer = new SILICONSCORERS::PS_Silicon_Rectangle("VamosSiliconScorer",0,Actar_NS::VamosSiliconHeight,Actar_NS::VamosSiliconWidth,1,1);
+ m_VamosSiliconScorer->RegisterPrimitive(VamosSiScorer);*/
+
G4VPrimitiveScorer* CsIScorer= new CALORIMETERSCORERS::PS_Calorimeter("CsI",NestingLevel);
m_CsIScorer->RegisterPrimitive(CsIScorer);
-
+
vector<int> level; level.push_back(0);
- G4VPrimitiveScorer* Actar_dig= new DRIFTELECTRONSCORERS::PS_DECathode("Actar_dig",0) ;
+ G4VPrimitiveScorer* Actar_dig= new TPCSCORERS::PS_TPCCathode("Actar_dig",0) ;
m_ActarScorer->RegisterPrimitive(Actar_dig);
-
+
G4SDManager::GetSDMpointer()->AddNewDetector(m_ActarScorer);
G4SDManager::GetSDMpointer()->AddNewDetector(m_SiliconScorer);
+ //G4SDManager::GetSDMpointer()->AddNewDetector(m_VamosSiliconScorer);
G4SDManager::GetSDMpointer()->AddNewDetector(m_CsIScorer) ;
}
@@ -686,6 +853,6 @@ extern"C" {
NPS::DetectorFactory::getInstance()->AddDetector("Actar",Actar::Construct);
}
};
-
+
proxy_nps_Actar p_nps_Actar;
}
diff --git a/NPSimulation/Detectors/Actar/Actar.hh b/NPSimulation/Detectors/Actar/Actar.hh
index 43aecfd54..8b6115c11 100644
--- a/NPSimulation/Detectors/Actar/Actar.hh
+++ b/NPSimulation/Detectors/Actar/Actar.hh
@@ -32,11 +32,23 @@ using namespace std;
#include "G4RotationMatrix.hh"
#include "G4LogicalVolume.hh"
#include "G4MultiFunctionalDetector.hh"
+#include "G4UserLimits.hh"
+#include "G4VFastSimulationModel.hh"
+#include "G4FastSimulationManager.hh"
// NPTool header
#include "NPSVDetector.hh"
#include "TActarData.h"
+#include "MEventReduced.h"
#include "NPInputParser.h"
+#include "Decay.hh"
+#include "BeamReaction.hh"
+
+#define NumberOfCobo 16
+#define NumberOfASAD 4
+#define NumberOfAGET 4
+#define NumberOfChannel 68
+
class Actar : public NPS::VDetector{
////////////////////////////////////////////////////
@@ -60,10 +72,15 @@ public:
private:
G4LogicalVolume* m_SquareDetector;
+ G4LogicalVolume* m_logicGas;
bool m_build_Silicon;
+ bool m_build_Vamos_Silicon;
bool m_build_CsI;
+ bool m_build_BeamDump;
G4LogicalVolume* m_LogicSilicon;
+ G4LogicalVolume* m_LogicVamosSilicon;
G4LogicalVolume* m_LogicCsICrystal;
+ G4LogicalVolume* m_LogicBeamDump;
////////////////////////////////////////////////////
////// Inherite from NPS::VDetector class /////////
@@ -94,11 +111,13 @@ public: // Scorer
G4MultiFunctionalDetector* m_ActarScorer ;
G4MultiFunctionalDetector* m_CsIScorer ;
G4MultiFunctionalDetector* m_SiliconScorer ;
+ //G4MultiFunctionalDetector* m_VamosSiliconScorer ;
////////////////////////////////////////////////////
///////////Event class to store Data////////////////
////////////////////////////////////////////////////
private:
TActarData* m_Event ;
+ MEventReduced* m_EventReduced;
////////////////////////////////////////////////////
///////////////Private intern Data//////////////////
@@ -110,12 +129,16 @@ private: // Geometry
vector<double> m_Phi;
// Shape type
- vector<string> m_Shape ;
+ vector<string> m_Shape;
// map
map<int, int> m_PadToXRow;
map<int, int> m_PadToZColumn;
-
+ int m_PadToCobo[128][128]; //!
+ int m_PadToAsad[128][128]; //!
+ int m_PadToAGET[128][128]; //!
+ int m_PadToChannel[128][128]; //!
+
// token
vector<string> m_GasMaterial;
vector<int> m_GasFraction;
@@ -125,11 +148,21 @@ private: // Geometry
// Visualisation Attribute
G4VisAttributes* m_VisChamber;
G4VisAttributes* m_VisWindows;
+ G4VisAttributes* m_VisMicromegas;
G4VisAttributes* m_VisGas;
G4VisAttributes* m_VisPads;
G4VisAttributes* m_SiliconVisAtt;
G4VisAttributes* m_CsIVisAtt;
+ G4VisAttributes* m_BeamDumpVisAtt;
+
// Needed for dynamic loading of the library
+
+private:
+ // Region were reaction can occure:
+ G4Region* m_ReactionRegion;
+ vector<G4VFastSimulationModel*> m_ReactionModel;
+
+
public:
static NPS::VDetector* Construct();
};
diff --git a/NPSimulation/EventGenerator/EventGeneratorIsotropic.cc b/NPSimulation/EventGenerator/EventGeneratorIsotropic.cc
index 96397b44b..f56243e19 100644
--- a/NPSimulation/EventGenerator/EventGeneratorIsotropic.cc
+++ b/NPSimulation/EventGenerator/EventGeneratorIsotropic.cc
@@ -41,6 +41,7 @@ using namespace CLHEP;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
EventGeneratorIsotropic::EventGeneratorIsotropic(){
m_ParticleStack = ParticleStack::getInstance();
+ event_ID=0;
}
@@ -59,7 +60,7 @@ EventGeneratorIsotropic::SourceParameters::SourceParameters(){
m_z0 = 0 ;
m_SigmaX = 0 ;
m_SigmaY = 0 ;
- m_particle = NULL;
+ m_particle = NULL;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -144,6 +145,10 @@ void EventGeneratorIsotropic::GenerateEvent(G4Event*){
G4double theta = acos(cos_theta) ;
G4double phi = RandFlat::shoot() * 2 * pi ;
G4double particle_energy = par.m_EnergyLow + RandFlat::shoot() * (par.m_EnergyHigh - par.m_EnergyLow) ;
+
+ cout << "Event ID= " << event_ID << " / theta= " << theta*180/3.1415 << " / energy= " << particle_energy << endl;
+ event_ID++;
+
// Direction of particle, energy and laboratory angle
G4double momentum_x = sin(theta) * cos(phi) ;
diff --git a/NPSimulation/EventGenerator/EventGeneratorIsotropic.hh b/NPSimulation/EventGenerator/EventGeneratorIsotropic.hh
index d509f8409..c3e18be43 100644
--- a/NPSimulation/EventGenerator/EventGeneratorIsotropic.hh
+++ b/NPSimulation/EventGenerator/EventGeneratorIsotropic.hh
@@ -46,6 +46,7 @@ public: // Inherit from VEventGenerator Class
void InitializeRootOutput() ;
private: // Source parameter from input file
+ G4int event_ID;
struct SourceParameters {
SourceParameters() ;
G4double m_EnergyLow ; // Lower limit of energy range
diff --git a/NPSimulation/Process/BeamReaction.cc b/NPSimulation/Process/BeamReaction.cc
index 5798ec4da..2a8921d10 100644
--- a/NPSimulation/Process/BeamReaction.cc
+++ b/NPSimulation/Process/BeamReaction.cc
@@ -101,7 +101,7 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
fastTrack.GetPrimaryTrack()->GetVolume()->GetLogicalVolume()->GetSolid();
double in = solid->DistanceToOut(P,V);
double out = solid->DistanceToOut(P,-V);
- double ratio = in / (out+in) ;
+ double ratio = in / (out+in) ;
if(out == 0){// first step of current event
rand = G4RandFlat::shoot();
@@ -119,7 +119,7 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
// shoot = false;
if(m_Reaction.IsAllowed(PrimaryTrack->GetKineticEnergy())){
return true;
-
+
}
else{
return false;
diff --git a/NPSimulation/Scorers/CMakeLists.txt b/NPSimulation/Scorers/CMakeLists.txt
index de493f8a6..a0fb63975 100644
--- a/NPSimulation/Scorers/CMakeLists.txt
+++ b/NPSimulation/Scorers/CMakeLists.txt
@@ -1,2 +1,2 @@
-add_library(NPSScorers SHARED NPSHitsMap.hh CalorimeterScorers.cc SiliconScorers.cc PhotoDiodeScorers.cc ObsoleteGeneralScorers.cc DriftElectronScorers.cc MDMScorer.cc )
+add_library(NPSScorers SHARED NPSHitsMap.hh CalorimeterScorers.cc SiliconScorers.cc PhotoDiodeScorers.cc ObsoleteGeneralScorers.cc DriftElectronScorers.cc TPCScorers.cc MDMScorer.cc )
target_link_libraries(NPSScorers ${ROOT_LIBRARIES} ${Geant4_LIBRARIES} ${NPLib_LIBRARIES} -lNPInitialConditions -lNPInteractionCoordinates)
diff --git a/NPSimulation/Scorers/TPCScorers.cc b/NPSimulation/Scorers/TPCScorers.cc
new file mode 100644
index 000000000..ab85f926d
--- /dev/null
+++ b/NPSimulation/Scorers/TPCScorers.cc
@@ -0,0 +1,119 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2016 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Pierre MORFOUACE contact address: morfouace@ganil.fr *
+ * *
+ * Creation Date : June 2018 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * Scorer for TPC Drift Electron collector *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * This new type of scorer to count drift electron in the cathode *
+ * *
+ *****************************************************************************/
+#include "TPCScorers.hh"
+#include "G4UnitsTable.hh"
+using namespace TPCSCORERS ;
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+PS_TPCCathode::PS_TPCCathode(G4String name,G4int depth)
+:G4VPrimitiveScorer(name, depth),HCID(-1){
+ m_Index = -1 ;
+ m_Level = depth;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+PS_TPCCathode::~PS_TPCCathode(){
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+G4bool PS_TPCCathode::ProcessHits(G4Step* aStep, G4TouchableHistory*){
+
+ // contain Energy Time, DetNbr, StripFront and StripBack
+ G4double* Infos = new G4double[10];
+ Infos[0] = 0;
+ Infos[1] = aStep->GetPreStepPoint()->GetProperTime();
+
+ m_DetectorNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(m_Level);
+ m_Position = aStep->GetPreStepPoint()->GetPosition();
+
+ // Interaction coordinates (used to fill the InteractionCoordinates branch)
+ Infos[2] = m_Position.x();
+ Infos[3] = m_Position.y();
+ Infos[4] = m_Position.z();
+ Infos[5] = m_Position.theta();
+ Infos[6] = m_Position.phi();
+ Infos[7] = m_DetectorNumber;
+
+
+ G4String PID = aStep->GetTrack()->GetDefinition()->GetParticleName();
+
+ m_Position = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(m_Position);
+
+
+ Infos[8] = (int)((m_Position.z()+128)/2);
+ Infos[9] = (int)((m_Position.x()+128)/2);
+
+ m_Index = Infos[8] + 128*Infos[9] ;
+
+
+ if(PID=="driftelectron"){
+ Infos[0] = aStep->GetTrack()->GetWeight();
+ }
+
+ // Check if the particle has interact before, if yes, add up the number of electron.
+ map<G4int, G4double**>::iterator it;
+ it= EvtMap->GetMap()->find(m_Index);
+ if(it!=EvtMap->GetMap()->end()){
+ G4double* dummy = *(it->second);
+ Infos[0]+=dummy[0];
+ Infos[1]=dummy[1];
+ }
+
+ EvtMap->set(m_Index, Infos);
+ return TRUE;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_TPCCathode::Initialize(G4HCofThisEvent* HCE){
+ EvtMap = new NPS::HitsMap<G4double*>(GetMultiFunctionalDetector()->GetName(), GetName());
+ if (HCID < 0) {
+ HCID = GetCollectionID(0);
+ }
+ HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_TPCCathode::EndOfEvent(G4HCofThisEvent*){
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_TPCCathode::clear(){
+ std::map<G4int, G4double**>::iterator MapIterator;
+ for (MapIterator = EvtMap->GetMap()->begin() ; MapIterator != EvtMap->GetMap()->end() ; MapIterator++){
+ delete *(MapIterator->second);
+ }
+
+ EvtMap->clear();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_TPCCathode::DrawAll(){
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_TPCCathode::PrintAll(){
+ G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl ;
+ G4cout << " PrimitiveScorer " << GetName() << G4endl ;
+ G4cout << " Number of entries " << EvtMap->entries() << G4endl ;
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/NPSimulation/Scorers/TPCScorers.hh b/NPSimulation/Scorers/TPCScorers.hh
new file mode 100644
index 000000000..6a3bc6db4
--- /dev/null
+++ b/NPSimulation/Scorers/TPCScorers.hh
@@ -0,0 +1,94 @@
+#ifndef DECathodeScorers_h
+/*****************************************************************************
+ * Copyright (C) 2009-2016 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Pierre MORFOUACE contact address: morfouace@ganil.fr *
+ * *
+ * Creation Date : July 2018 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * File old the scorer specific to the Silicon Detector *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * This new style of scorer is aim to become the standard way of doing scorer*
+ * in NPTool. *
+ *The index is build using the TrackID, Detector Number and Strip Number. *
+ *The scorer Hold Energy and time together *
+ *Only one scorer is needed for a detector *
+ *****************************************************************************/
+#include "G4VPrimitiveScorer.hh"
+#include "NPSHitsMap.hh"
+
+#include <map>
+using namespace std;
+using namespace CLHEP;
+
+namespace TPCSCORERS {
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+ class PS_TPCCathode : public G4VPrimitiveScorer{
+
+ public: // with description
+ PS_TPCCathode(G4String name,G4int depth=0);
+ ~PS_TPCCathode();
+
+ protected: // with description
+ G4bool ProcessHits(G4Step*, G4TouchableHistory*);
+
+ public:
+ void Initialize(G4HCofThisEvent*);
+ void EndOfEvent(G4HCofThisEvent*);
+ void clear();
+ void DrawAll();
+ void PrintAll();
+
+private: // inherited from G4VPrimitiveScorer
+ G4int HCID;
+ NPS::HitsMap<G4double*>* EvtMap;
+
+ private: // Needed for intermediate calculation (avoid multiple instantiation in Processing Hit)
+ G4int m_Index;
+ G4int m_Level;
+ G4int m_DetectorNumber;
+ G4ThreeVector m_Position;
+
+ };
+////////////////////////////////////////////////////////////////////////////////
+ class PS_DEDigitizer : public G4VPrimitiveScorer{
+
+ public: // with description
+ PS_DEDigitizer(G4String name,G4int depth=0);
+ ~PS_DEDigitizer();
+
+ protected: // with description
+ G4bool ProcessHits(G4Step*, G4TouchableHistory*);
+
+ public:
+ void Initialize(G4HCofThisEvent*);
+ void EndOfEvent(G4HCofThisEvent*);
+ void clear();
+ void DrawAll();
+ void PrintAll();
+
+private: // inherited from G4VPrimitiveScorer
+ G4int HCID;
+ NPS::HitsMap<G4double*>* EvtMap;
+
+ private: // Needed for intermediate calculation (avoid multiple instantiation in Processing Hit)
+ G4int m_Index;
+ G4int m_Level;
+ G4int m_DetectorNumber;
+ G4ThreeVector m_Position;
+
+ };
+
+
+}
+
+#endif
diff --git a/Projects/Actar/Analysis.cxx b/Projects/Actar/Analysis.cxx
index 23c45664f..8effa663b 100644
--- a/Projects/Actar/Analysis.cxx
+++ b/Projects/Actar/Analysis.cxx
@@ -6,13 +6,13 @@
*****************************************************************************/
/*****************************************************************************
- * Original Author: XAUTHORX contact address: XMAILX *
+ * Original Author: Morfouace Pierre contact address: morfouace@ganil.fr *
* *
- * Creation Date : XMONTHX XYEARX *
+ * Creation Date : April 2018 *
* Last update : *
*---------------------------------------------------------------------------*
* Decription: *
- * This class describe Actar analysis project *
+ * This class describe Actar analysis project *
* *
*---------------------------------------------------------------------------*
* Comment: *
@@ -20,10 +20,15 @@
*****************************************************************************/
#include<iostream>
+
using namespace std;
#include"Analysis.h"
#include"NPAnalysisFactory.h"
#include"NPDetectorManager.h"
+#include"RootOutput.h"
+#include"RootInput.h"
+
+
////////////////////////////////////////////////////////////////////////////////
Analysis::Analysis(){
}
@@ -33,39 +38,168 @@ Analysis::~Analysis(){
////////////////////////////////////////////////////////////////////////////////
void Analysis::Init(){
- Actar= (TActarPhysics*) m_DetectorManager->GetDetector("Actar");
+ Actar= (TActarPhysics*) m_DetectorManager->GetDetector("Actar");
+
+ Actar->ReadAnalysisConfig();
+ if(Actar->GetRansacStatus()){
+ Actar->SetRansacParameter("./configs/RansacConfig.dat");
+ }
+
+ DriftVelocity = Actar->GetDriftVelocity();
+ PadSizeX = Actar->GetPadSizeX();
+ PadSizeY = Actar->GetPadSizeY();
+ NumberOfPadsX = Actar->GetNumberOfPadsX();
+ NumberOfPadsY = Actar->GetNumberOfPadsY();
+
+ EnergyLoss_3He = NPL::EnergyLoss("EnergyLossTable/He3_D2_6.24151e+08_295.G4table","G4Table",100);
+ EnergyLoss_17C = NPL::EnergyLoss("EnergyLossTable/C17_D2_6.24151e+08_295.G4table","G4Table",100);
+ TheReaction = new NPL::Reaction("17C(d,3He)16B@510");
+
+ InitOutputBranch();
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::TreatEvent(){
- //cout << "/// Size=" << Actar->PadRow.size() << endl;
- for(unsigned int i=0; i<Actar->PadRow.size(); i++){
- //cout << "Row= " << Actar->PadRow[i] << endl;
- }
+ ReInitValue();
+
+ int TrackMult = Actar->GetTrackMult();
+
+ TVector3 vX = TVector3(1,0,0);
+ TVector3 aTrack, vB;
+
+ if(TrackMult>1){
+ vTrack = Actar->GetTracks();
+ double scalarproduct=0;
+ int BeamTrack=0;
+ for(unsigned int i=0; i<TrackMult; i++){
+ TVector3 vtest = TVector3(vTrack[i].GetDirectionVector().X(),vTrack[i].GetDirectionVector().Y(),vTrack[i].GetDirectionVector().Z());
+ TVector3 vunit = vtest.Unit();
+ double scalar = abs(vunit.Dot(vX));
+ vScalar.push_back(scalar);
+ //cout << scalar << endl;
+ //cout << scalarproduct << endl;
+ if(scalar>scalarproduct){
+ BeamTrack=i;
+ scalarproduct=scalar;
+ }
+ }
+
+ double XBeam = vTrack[BeamTrack].GetDirectionVector().X();
+ double YBeam = vTrack[BeamTrack].GetDirectionVector().Y();
+ double ZBeam = vTrack[BeamTrack].GetDirectionVector().Z();
+ TVector3 vBeam = TVector3(XBeam,YBeam,ZBeam);
+
+ double XBeamPoint = vTrack[BeamTrack].GetXh();
+ double YBeamPoint = vTrack[BeamTrack].GetYh();
+ double ZBeamPoint = vTrack[BeamTrack].GetZh();
+ TVector3 vBeamPos = TVector3(XBeamPoint,YBeamPoint,ZBeamPoint);
+
+ vB = TVector3(XBeam*PadSizeX, YBeam*PadSizeY,ZBeam*DriftVelocity);
+ BeamAngle = (vX.Angle(vB))*180/TMath::Pi();
+
+ for(unsigned int i=0; i<TrackMult; i++){
+ if(i!=BeamTrack){
+ double Xdir = vTrack[i].GetDirectionVector().X();
+ double Ydir = vTrack[i].GetDirectionVector().Y();
+ double Zdir = vTrack[i].GetDirectionVector().Z();
+
+ XVertex.push_back(vTrack[i].GetVertexPostion(vBeam,vBeamPos).X()*PadSizeX);
+ YVertex.push_back(vTrack[i].GetVertexPostion(vBeam,vBeamPos).Y()*PadSizeY);
+ ZVertex.push_back(vTrack[i].GetVertexPostion(vBeam,vBeamPos).Z()*DriftVelocity);
+
+ aTrack = TVector3(Xdir*PadSizeX, Ydir*PadSizeY, Zdir*DriftVelocity);
+
+
+ double angle = vX.Angle(aTrack)*180/TMath::Pi();
+ //double angle = vX.Angle(aTrack)*180/TMath::Pi();
+ if(angle>90) angle = 180-angle;
+
+ double x1 = vTrack[i].GetXm()*PadSizeX;
+ double x2 = vTrack[i].GetXh()*PadSizeX;
+ double y1 = vTrack[i].GetYm()*PadSizeY-0.5*NumberOfPadsY*PadSizeY;
+ double y2 = vTrack[i].GetYh()*PadSizeY-0.5*NumberOfPadsY*PadSizeY;
+ double z1 = -(vTrack[i].GetZm()-256)*DriftVelocity;
+ double z2 = -(vTrack[i].GetZh()-256)*DriftVelocity;
+ //if(vScalar[i]<0.999)GetMayaSiHitPosition(x1,x2,y1,y2,z1,z2);
+
+ if(XVertex[i]>0 && XVertex[i]<256){
+ double SiDistanceToPadPlane = 47*mm;
+ double LengthInGas = 256-XVertex[i] + SiDistanceToPadPlane;
+ for(unsigned int k=0; k<Actar->Si_E.size(); k++){
+ ESi.push_back(Actar->Si_E[k]);
+ DE.push_back(vTrack[i].GetTotalCharge());
+ double E3 = EnergyLoss_3He.EvaluateInitialEnergy(Actar->Si_E[k]*MeV,LengthInGas*mm,angle*TMath::Pi()/180);
+ double BeamEnergy = EnergyLoss_17C.Slow(510*MeV,(XVertex[i]+60)*mm, BeamAngle*TMath::Pi()/180);
+ TheReaction->SetBeamEnergy(BeamEnergy);
+ ELab.push_back(E3);
+ ThetaLab.push_back(angle);
+ TheReaction->SetNuclei3(E3,angle*TMath::Pi()/180);
+ Ex.push_back(TheReaction->GetExcitation4());
+ ThetaCM.push_back(TheReaction->GetThetaCM()*180./TMath::Pi());
+ }
+ }
+ }
+ }
+ }
+
+ //for(unsigned int i=0; i<Actar->PadX.size(); i++){
+ //cout << "X= " << Actar->PadX[i] << endl;
+ //}
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::End(){
}
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::InitOutputBranch() {
+ RootOutput::getInstance()->GetTree()->Branch("DE",&DE);
+ RootOutput::getInstance()->GetTree()->Branch("ESi",&ESi);
+ RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab);
+ RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab);
+ RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex);
+ RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM);
+ RootOutput::getInstance()->GetTree()->Branch("vScalar",&vScalar);
+ RootOutput::getInstance()->GetTree()->Branch("XVertex",&XVertex);
+ RootOutput::getInstance()->GetTree()->Branch("YVertex",&YVertex);
+ RootOutput::getInstance()->GetTree()->Branch("ZVertex",&ZVertex);
+ RootOutput::getInstance()->GetTree()->Branch("BeamAngle",&BeamAngle,"BeamAngle/D");
+
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::ReInitValue(){
+ DE.clear();
+ ESi.clear();
+ ELab.clear();
+ ThetaLab.clear();
+ Ex.clear();
+ ThetaCM.clear();
+ vScalar.clear();
+ XVertex.clear();
+ YVertex.clear();
+ ZVertex.clear();
+ BeamAngle=-1000;
+
+}
////////////////////////////////////////////////////////////////////////////////
// Construct Method to be pass to the DetectorFactory //
////////////////////////////////////////////////////////////////////////////////
NPL::VAnalysis* Analysis::Construct(){
- return (NPL::VAnalysis*) new Analysis();
+ return (NPL::VAnalysis*) new Analysis();
}
////////////////////////////////////////////////////////////////////////////////
// Registering the construct method to the factory //
////////////////////////////////////////////////////////////////////////////////
extern "C"{
-class proxy{
- public:
- proxy(){
- NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct);
- }
-};
-
-proxy p;
+ class proxy{
+ public:
+ proxy(){
+ NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct);
+ }
+ };
+
+ proxy p;
}
diff --git a/Projects/Actar/Analysis.h b/Projects/Actar/Analysis.h
index 075b9abe6..b49b5f0a6 100644
--- a/Projects/Actar/Analysis.h
+++ b/Projects/Actar/Analysis.h
@@ -1,4 +1,4 @@
-#ifndef Analysis_h
+#ifndef Analysis_h
#define Analysis_h
/*****************************************************************************
* Copyright (C) 2009-2016 this file is part of the NPTool Project *
@@ -21,22 +21,58 @@
* *
*****************************************************************************/
-#include"NPVAnalysis.h"
-#include"TActarPhysics.h"
+#include "NPVAnalysis.h"
+#include "TActarPhysics.h"
+#include "NPEnergyLoss.h"
+#include "NPReaction.h"
+#include "NPTrack.h"
+
+
+
class Analysis: public NPL::VAnalysis{
- public:
+public:
Analysis();
~Analysis();
-
- public:
+
+public:
void Init();
void TreatEvent();
void End();
-
- static NPL::VAnalysis* Construct();
-
- private:
- TActarPhysics* Actar;
-
+ void InitOutputBranch();
+ void ReInitValue();
+
+
+ static NPL::VAnalysis* Construct();
+
+public:
+ double DriftVelocity;
+ double PadSizeX;
+ double PadSizeY;
+ int NumberOfPadsX;
+ int NumberOfPadsY;
+
+
+private:
+ TActarPhysics* Actar;
+
+ vector<NPL::Track> vTrack;
+
+ double BeamAngle;
+ vector<double> vScalar;
+ vector<double> ThetaLab;
+ vector<double> ELab;
+ vector<double> ESi;
+ vector<double> DE;
+ vector<double> Ex;
+ vector<double> ThetaCM;
+ vector<double> XVertex;
+ vector<double> YVertex;
+ vector<double> ZVertex;
+
+ NPL::EnergyLoss EnergyLoss_3He;
+ NPL::EnergyLoss EnergyLoss_17C;
+ NPL::Reaction* TheReaction;
+
+
};
#endif
--
GitLab