Commit 32ba5a25 authored by Adrien Matta's avatar Adrien Matta
Browse files

*Adding Samurai BDC

parent 2ff46ad5
Pipeline #118638 passed with stages
in 9 minutes and 38 seconds
......@@ -2,7 +2,9 @@ add_custom_command(OUTPUT TSamuraiHodoscopeDataDict.cxx COMMAND ${CMAKE_BINARY_D
add_custom_command(OUTPUT TSamuraiHodoscopePhysicsDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TSamuraiHodoscopePhysics.h TSamuraiHodoscopePhysicsDict.cxx TSamuraiHodoscopePhysics.rootmap libNPSamurai.dylib DEPENDS TSamuraiHodoscopePhysics.h)
add_custom_command(OUTPUT TSamuraiBDCDataDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TSamuraiBDCData.h TSamuraiBDCDataDict.cxx TSamuraiBDCData.rootmap libNPSamurai.dylib DEPENDS TSamuraiBDCData.h)
add_custom_command(OUTPUT TSamuraiBDCPhysicsDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TSamuraiBDCPhysics.h TSamuraiBDCPhysicsDict.cxx TSamuraiBDCPhysics.rootmap libNPSamurai.dylib DEPENDS TSamuraiBDCPhysics.h)
add_custom_command(OUTPUT TSamuraiFDC2DataDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TSamuraiFDC2Data.h TSamuraiFDC2DataDict.cxx TSamuraiFDC2Data.rootmap libNPSamurai.dylib DEPENDS TSamuraiFDC2Data.h)
......@@ -12,8 +14,8 @@ add_custom_command(OUTPUT TSamuraiFDC0DataDict.cxx COMMAND ${CMAKE_BINARY_DIR}/s
add_custom_command(OUTPUT TSamuraiFDC0PhysicsDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TSamuraiFDC0Physics.h TSamuraiFDC0PhysicsDict.cxx TSamuraiFDC0Physics.rootmap libNPSamurai.dylib DEPENDS TSamuraiFDC0Physics.h)
add_library(NPSamurai SHARED TSamuraiHodoscopeData.cxx TSamuraiFDC2Data.cxx TSamuraiFDC2DataDict.cxx TSamuraiHodoscopePhysics.cxx TSamuraiHodoscopePhysicsDict.cxx TSamuraiFDC2Physics.cxx TSamuraiHodoscopeDataDict.cxx TSamuraiFDC2PhysicsDict.cxx TSamuraiFDC0Data.cxx TSamuraiFDC0DataDict.cxx TSamuraiFDC0Physics.cxx TSamuraiFDC0PhysicsDict.cxx)
add_library(NPSamurai SHARED TSamuraiHodoscopeData.cxx TSamuraiBDCData.cxx TSamuraiBDCPhysics.cxx TSamuraiBDCDataDict.cxx TSamuraiBDCPhysicsDict.cxx TSamuraiFDC2Data.cxx TSamuraiFDC2DataDict.cxx TSamuraiHodoscopePhysics.cxx TSamuraiHodoscopePhysicsDict.cxx TSamuraiFDC2Physics.cxx TSamuraiHodoscopeDataDict.cxx TSamuraiFDC2PhysicsDict.cxx TSamuraiFDC0Data.cxx TSamuraiFDC0DataDict.cxx TSamuraiFDC0Physics.cxx TSamuraiFDC0PhysicsDict.cxx)
target_link_libraries(NPSamurai ${ROOT_LIBRARIES} NPCore NPTrackReconstruction)
install(FILES TSamuraiHodoscopeData.h TSamuraiHodoscopePhysics.h TSamuraiFDC2Data.h TSamuraiFDC2Physics.h TSamuraiFDC0Data.h TSamuraiFDC0Physics.h SamuraiDCIndex.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
install(FILES TSamuraiBDCData.h TSamuraiBDCPhysics.h TSamuraiHodoscopeData.h TSamuraiHodoscopePhysics.h TSamuraiFDC2Data.h TSamuraiFDC2Physics.h TSamuraiFDC0Data.h TSamuraiFDC0Physics.h SamuraiDCIndex.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
#include "TSamuraiBDCData.h"
#include <iostream>
TSamuraiBDCData::TSamuraiBDCData(){};
TSamuraiBDCData::~TSamuraiBDCData(){};
////////////////////////////////////////////////////////////////////////////////
void TSamuraiBDCData::SetData(const int& Det, const int& Layer, const int& Wire, const double& Time, const int& Edge){
fBDC_DetectorNbr.push_back(Det);
fBDC_LayerNbr.push_back(Layer);
fBDC_WireNbr.push_back(Wire);
fBDC_Time.push_back(Time);
fBDC_Edge.push_back(Edge);
}
////////////////////////////////////////////////////////////////////////////////
void TSamuraiBDCData::Clear(){
fBDC_DetectorNbr.clear();
fBDC_LayerNbr.clear();
fBDC_WireNbr.clear();
fBDC_Time.clear();
fBDC_Edge.clear();
}
////////////////////////////////////////////////////////////////////////////////
void TSamuraiBDCData::Print(){
using namespace std;
cout << " -- Event:" << endl;
cout << " - Multiplicity: " << Mult() << endl;
}
////////////////////////////////////////////////////////////////////////////////
unsigned int TSamuraiBDCData::MultLayer(unsigned int det , unsigned int layer, int edge){
unsigned int mult=0;
unsigned int size = fBDC_DetectorNbr.size();
for(unsigned int i = 0 ; i< size ; i++ ){
if(fBDC_DetectorNbr[i]==det)
if(fBDC_LayerNbr[i]==layer)
if(fBDC_Edge[i]==edge && edge!=-1) // edge type is specified (0 or 1)
mult++;
else if(edge==-1)// edge type is not specified
mult++;
}
return mult;
}
////////////////////////////////////////////////////////////////////////////////
std::vector<int> TSamuraiBDCData::GetWire(unsigned int det , unsigned int layer){
std::vector<int> wires;
unsigned int size = fBDC_DetectorNbr.size();
for(unsigned int i = 0 ; i< size ; i++ ){
if(fBDC_DetectorNbr[i]==det)
if(fBDC_LayerNbr[i]==layer)
wires.push_back(fBDC_WireNbr[i]);
}
return wires;
}
ClassImp(TSamuraiBDCData);
#ifndef TSamuraiBDCData_H
#define TSamuraiBDCData_H
#include "TObject.h"
#include <vector>
class TSamuraiBDCData: public TObject{
public:
TSamuraiBDCData();
~TSamuraiBDCData();
private:
std::vector<int> fBDC_DetectorNbr;
std::vector<int> fBDC_LayerNbr;
std::vector<int> fBDC_WireNbr;
std::vector<double> fBDC_Time;
std::vector<int> fBDC_Edge;
public:
void Clear();
void Print();
void Clear(const Option_t*) {};
void Dump() const{};
public:
void SetData(const int& Det, const int& Layer, const int& Wire, const double& Time, const int& Edge);
unsigned int Mult(){return fBDC_DetectorNbr.size();};
unsigned int MultLayer(unsigned int det , unsigned int layer, int edge=-1);
std::vector<int> GetWire(unsigned int det , unsigned int layer);
int const GetDetectorNbr(const unsigned int& i){return fBDC_DetectorNbr[i];};
int const GetLayerNbr(const unsigned int& i){return fBDC_LayerNbr[i];};
int const GetWireNbr(const unsigned int& i){return fBDC_WireNbr[i];};
double const GetTime(const unsigned int& i){return fBDC_Time[i];};
int const GetEdge(const unsigned int& i){return fBDC_Edge[i];};
ClassDef(TSamuraiBDCData,1);
};
#endif
/*****************************************************************************
* Copyright (C) 2009-2020 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: Adrien MATTA contact address: matta@lpccaen.in2p3.fr *
* *
* Creation Date : May 2021 *
* Last update : *
*---------------------------------------------------------------------------*
* Decription: *
* This class hold SamuraiBDC treated data *
* *
*---------------------------------------------------------------------------*
* Comment: *
* *
*****************************************************************************/
#include "TSamuraiBDCPhysics.h"
// STL
#include <sstream>
#include <iostream>
#include <cmath>
#include <stdlib.h>
#include <limits>
using namespace std;
// NPL
#include "RootInput.h"
#include "RootOutput.h"
#include "TAsciiFile.h"
#include "NPOptionManager.h"
#include "NPDetectorFactory.h"
#include "NPSystemOfUnits.h"
// ROOT
using namespace NPUNITS;
///////////////////////////////////////////////////////////////////////////
ClassImp(TSamuraiBDCPhysics)
///////////////////////////////////////////////////////////////////////////
TSamuraiBDCPhysics::TSamuraiBDCPhysics(){
m_EventData = new TSamuraiBDCData ;
m_PreTreatedData = new TSamuraiBDCData ;
m_EventPhysics = this ;
//m_Spectra = NULL;
ToTThreshold_L = 0;
ToTThreshold_H = 180;
DriftLowThreshold=0.1 ;
DriftUpThreshold=2.4;
PowerThreshold=5;
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiBDCPhysics::BuildSimplePhysicalEvent(){
BuildPhysicalEvent();
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiBDCPhysics::BuildPhysicalEvent(){
PreTreat();
// RemoveNoise();
// Map[plane angle, vector of spatial information]
static map<double, vector<double> > X ;
static map<double, vector<double> > Z ;
static map<double, vector<double> > R ;
static int det,layer,wire;
X.clear();Z.clear();R.clear();
unsigned int size = Detector.size();
for(unsigned int i = 0 ; i < size ; i++){
if(DriftLength[i] > DriftLowThreshold && DriftLength[i] < DriftUpThreshold){
det = Detector[i];
layer = Layer[i];
wire = Wire[i];
SamuraiDCIndex idx(det,layer,wire);
X[Wire_Angle[idx]].push_back(Wire_X[idx]);
Z[Wire_Angle[idx]].push_back(Wire_Z[idx]);
R[Wire_Angle[idx]].push_back(DriftLength[i]);
}
}
// Reconstruct the vector for each of the plane of each of the detector
static double X0,X100,a,b; // store the BuildTrack2D results
static map<double, TVector3 > VX0 ;
static map<double, TVector3 > VX100 ;
static map<double, double > D ;// the minimum distance
static unsigned int uid; uid=0;
VX0.clear();VX100.clear(),D.clear();
for(auto it = X.begin();it!=X.end();++it){
#if __cplusplus > 199711L && NPMULTITHREADING
m_reconstruction.AddPlan(uid++,X[it->first],Z[it->first],R[it->first]);
#else
D[it->first]=m_reconstruction.BuildTrack2D(X[it->first],Z[it->first],R[it->first],X0,X100,a,b);
#endif
}
#if __cplusplus > 199711L && NPMULTITHREADING
// do all plan at once in parallele, return when all plan are done
m_reconstruction.BuildTrack2D();
uid=0;
#endif
for(auto it = X.begin();it!=X.end();++it){
#if __cplusplus > 199711L && NPMULTITHREADING
D[it->first]=m_reconstruction.GetResults(uid++,X0,X100,a,b);
#endif
// for Debug, write a file of
/* { std::ofstream f("distance.txt", std::ios::app);
f<< D[it->first] << endl;
f.close();
}
*/
// very large a means track perpendicular to the chamber, what happen when there is pile up
if(abs(a)>5000)
PileUp++;
Mult+=X[it->first].size();
// Position at z=0
TVector3 P(X0,0,0);
P.RotateZ(it->first);
VX0[it->first]=P;
// Position at z=100
TVector3 P100= TVector3(X100,0,0);
P100.RotateZ(it->first);
VX100[it->first]=P100;
}
// Reconstruct the central position (z=0) for each detector
static vector<TVector3> C ;
static vector<double > W ; // weight based on D
C.clear(),W.clear();
TVector3 P;
for(auto it1 = VX0.begin();it1!=VX0.end();++it1){
for(auto it2 = it1;it2!=VX0.end();++it2){
if(it1!=it2){// different plane
//cout << "BDC" << endl;
m_reconstruction.ResolvePlane(it1->second,it1->first,it2->second,it2->first,P);
// cout << "done " << D[it1->first] << " " << D[it2->first] << endl;
if(P.X()!=-10000 && D[it1->first]<PowerThreshold&& D[it2->first]<PowerThreshold){
C.push_back(P);
// Mean pos are weighted based on the the sum of distance from track
// to hit obtained during the minimisation
W.push_back(1./sqrt(D[it1->first]*D[it2->first]));
}
}
}
}
// Reconstruct the position at z=100 for each detector
static vector<TVector3> C100 ;
C100.clear();
for(auto it1 = VX100.begin();it1!=VX100.end();++it1){
for(auto it2 = it1;it2!=VX100.end();++it2){
if(it1!=it2 ){// different plane, same detector
m_reconstruction.ResolvePlane(it1->second,it1->first,it2->second,it2->first,P);
if(P.X()!=-10000&& D[it1->first]<PowerThreshold && D[it2->first]<PowerThreshold)
C100.push_back(P);
}
}
}
// Build the Reference position by averaging all possible pair
size = C.size();
static double PosX100,PosY100,norm;
if(size){
PosX=0;
PosY=0;
PosX100=0;
PosY100=0;
norm=0;
for(unsigned int i = 0 ; i < size ; i++){
PosX+= C[i].X()*W[i];
PosY+= C[i].Y()*W[i];
PosX100+= C100[i].X()*W[i];
PosY100+= C100[i].Y()*W[i];
norm+=W[i];
}
MultMean=size;
// Mean position at Z=0
PosX=PosX/norm;
PosY=PosY/norm;
// Mean position at Z=100
PosX100=PosX100/norm;
PosY100=PosY100/norm;
devX=0;
devY=0;
for(unsigned int i = 0 ; i < size ; i++){
devX+=W[i]*(C[i].X()-PosX)*(C[i].X()-PosX);
devY+=W[i]*(C[i].Y()-PosY)*(C[i].Y()-PosY);
}
devX=sqrt(devX/((size-1)*norm));
devY=sqrt(devY/((size-1)*norm));
// Compute ThetaX, angle between the Direction vector projection in XZ with
// the Z axis
//ThetaX=atan((PosX100-PosX)/100.);
ThetaX = (PosX100-PosX)/100.;
// Compute PhiY, angle between the Direction vector projection in YZ with
// the Z axis
//PhiY=atan((PosY100-PosY)/100.);
PhiY=(PosY100-PosY)/100.;
Dir=TVector3(PosX100-PosX,PosY100-PosY,100).Unit();
}
return;
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiBDCPhysics::PreTreat(){
ClearPreTreatedData();
static CalibrationManager* Cal = CalibrationManager::getInstance();
static string channel;
// one map per detector
map<unsigned int, vector<double> > X ;
map<unsigned int, vector<double> > Z ;
map<unsigned int, vector<double> > R ;
unsigned int size = m_EventData->Mult();
for(unsigned int i = 0 ; i < size ; i++){
// EDGE=1 is the leading edge, IE, the real time.
// EDGE=0 is the trailing edge, so it helps build Tot
if(m_EventData->GetEdge(i)==1){
int det = m_EventData->GetDetectorNbr(i);
int layer = m_EventData->GetLayerNbr(i);
int wire = m_EventData->GetWireNbr(i);
double time = m_EventData->GetTime(i);
double etime = 0;
// look for matching trailing edge
for(unsigned int j = 0 ; j < size ; j++){
if(m_EventData->GetEdge(j)==0){
int edet = m_EventData->GetDetectorNbr(j);
int elayer = m_EventData->GetLayerNbr(j);
int ewire = m_EventData->GetWireNbr(j);
// same wire
if(wire==ewire && layer==elayer && det==edet){
etime = m_EventData->GetTime(j);
}
}
if(etime && etime>time)
break;
else
etime=0;
}
// a valid wire must have an edge
if(etime && time && etime-time>ToTThreshold_L && etime-time<ToTThreshold_H){
Detector.push_back(det);
Layer.push_back(layer);
Wire.push_back(wire);
Time.push_back(time);
ToT.push_back(etime-time);
channel="SamuraiBDC/L" + NPL::itoa(layer);
DriftLength.push_back(2.5-Cal->ApplySigmoid(channel,etime));
}
}
}
return;
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiBDCPhysics::RemoveNoise(){
// Remove the noise by looking if a matching wire exist in the adjacent layer
// this done by looking at the closest plane with the same orientation
unsigned int size = Detector.size();
for(unsigned int i = 0 ; i < size ; i++){
bool match=false;
int det = Detector[i];
int layer = Layer[i];
int wire = Wire[i];
// look for matching adjacent wire
for(unsigned int j = 0 ; j < size ; j++){
int adet = Detector[j];
int alayer = Layer[j];
int awire = Wire[j];
bool blayer = false;
if(layer%2==0){
if(layer+1==alayer)
blayer=true;
}
else{
if(layer-1==alayer)
blayer=true;
}
if(det==adet && blayer && abs(wire-awire)<=1){
match=true;
break;
}
}
if(match)
Matched.push_back(true);
else
Matched.push_back(false);
}
return;
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiBDCPhysics::Clear(){
MultMean=0;
PileUp=0;
Mult=0;
PosX=PosY=-10000;
devX=devY=-10000;
DriftLength.clear();
Detector.clear();
Layer.clear();
Wire.clear();
Time.clear();
ToT.clear();
ParticleDirection.clear();
MiddlePosition.clear();
Matched.clear();
}
///////////////////////////////////////////////////////////////////////////
//// Innherited from VDetector Class ////
///////////////////////////////////////////////////////////////////////////
void TSamuraiBDCPhysics::ReadConfiguration(NPL::InputParser parser){
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("SAMURAIBDC");
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << "//// " << blocks.size() << " detector(s) found " << endl;
vector<string> token= {"XML"};
for(unsigned int i = 0 ; i < blocks.size() ; i++){
cout << endl << "//// Samurai BDC (" << i+1 << ")" << endl;
string xmlpath = blocks[i]->GetString("XML");
NPL::XmlParser xml;
xml.LoadFile(xmlpath);
AddDC("SAMURAIBDC",xml);
}
#if __cplusplus > 199711L && NPMULTITHREADING
if(blocks.size()){ // if a detector is found, init the thread pool
// one thread for each plan X,Y = 2
// ! more than that this will not help !
m_reconstruction.SetNumberOfThread(2);
m_reconstruction.InitThreadPool();
}
#endif
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiBDCPhysics::AddDC(string name, NPL::XmlParser& xml){
std::vector<NPL::XML::block*> b = xml.GetAllBlocksWithName(name);
// BDC case
if(name=="SAMURAIBDC"){
unsigned int det=0;
unsigned int sizeB = b.size();
for(unsigned int i = 0 ; i < sizeB ; i++){
unsigned int layer = b[i]->AsInt("layer");
unsigned int wire = b[i]->AsInt("wireid");
double X = b[i]->AsDouble("wirepos");
double Z = b[i]->AsDouble("wirez");
string sDir = b[i]->AsString("anodedir");
double T=0;
if(sDir=="X")
T= 0*deg;
else if(sDir=="Y")
T= -90*deg;
else if(sDir=="U")
T=-30*deg;
else if(sDir=="V")
T=+30*deg;
else{
cout << "ERROR: Unknown layer orientation for Samurai BDC"<< endl;
exit(1);
}
SamuraiDCIndex idx(det,layer,wire);
Wire_X[idx]=X;
Wire_Z[idx]=Z;
Wire_Angle[idx]=T;
}
}
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiBDCPhysics::InitSpectra(){
//m_Spectra = new TSamuraiBDCSpectra(m_NumberOfDetector);
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiBDCPhysics::FillSpectra(){
// m_Spectra -> FillRawSpectra(m_EventData);
// m_Spectra -> FillPreTreatedSpectra(m_PreTreatedData);
// m_Spectra -> FillPhysicsSpectra(m_EventPhysics);
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiBDCPhysics::CheckSpectra(){
// m_Spectra->CheckSpectra();
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiBDCPhysics::ClearSpectra(){
// To be done
}
///////////////////////////////////////////////////////////////////////////
map< string , TH1*> TSamuraiBDCPhysics::GetSpectra() {
/* if(m_Spectra)
return m_Spectra->GetMapHisto();
else{
map< string , TH1*> empty;
return empty;
}*/
map< string , TH1*> empty;
return empty;
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiBDCPhysics::WriteSpectra(){
// m_Spectra->WriteSpectra();
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiBDCPhysics::AddParameterToCalibrationManager(){
CalibrationManager* Cal = CalibrationManager::getInstance();
// each layer
for( int l = 0 ; l < 8 ; ++l){
Cal->AddParameter("SamuraiBDC", "L"+ NPL::itoa(l),"BDC_L"+ NPL::itoa(l));
}
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiBDCPhysics::InitializeRootInputRaw(){
TChain* inputChain = RootInput::getInstance()->GetChain() ;
inputChain->SetBranchStatus( "SamuraiBDC" , true );
// The following line is necessary only for system were the tree is splitted
// (older root version). The found argument silenced the Branches not found
// warning for non splitted tree.
if(inputChain->FindBranch("fDC_*"))
inputChain->SetBranchStatus( "fDC_*",true);
inputChain->SetBranchAddress( "SamuraiBDC" , &m_EventData );
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiBDCPhysics::InitializeRootInputPhysics(){
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiBDCPhysics::InitializeRootOutput(){
TTree* outputTree = RootOutput::getInstance()->GetTree("SamuraiBDC");
outputTree->Branch( "SamuraiBDC" , "TSamuraiBDCPhysics" , &m_EventPhysics );
}
////////////////////////////////////////////////////////////////////////////////
// Construct Method to be pass to the DetectorFactory //
////////////////////////////////////////////////////////////////////////////////
NPL::VDetector* TSamuraiBDCPhysics::Construct(){
return (NPL::VDetector*) new TSamuraiBDCPhysics();
}
////////////////////////////////////////////////////////////////////////////////
// Registering the construct method to the factory //
////////////////////////////////////////////////////////////////////////////////
extern "C"{
class proxy_samuraiBDC{
public:
proxy_samuraiBDC(){
NPL::DetectorFactory::getInstance()->AddToken("SAMURAIBDC","Samurai");
NPL::DetectorFactory::getInstance()->AddToken("SAMURAIBDC1","Samurai");
NPL::DetectorFactory::getInstance()->