Commit 4e57b61e authored by Adrien Matta's avatar Adrien Matta
Browse files

* adding FDC0

parent c52ab052
......@@ -92,7 +92,7 @@ namespace NPL{
void InitThreadPool();
bool IsDone();
#endif
#endif
private: // Target property
double m_TargetThickness;
......
......@@ -2,8 +2,13 @@ add_custom_command(OUTPUT TSamuraiFDC2DataDict.cxx COMMAND ${CMAKE_BINARY_DIR}/s
add_custom_command(OUTPUT TSamuraiFDC2PhysicsDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TSamuraiFDC2Physics.h TSamuraiFDC2PhysicsDict.cxx TSamuraiFDC2Physics.rootmap libNPSamurai.dylib DEPENDS TSamuraiFDC2Physics.h)
add_library(NPSamurai SHARED TSamuraiFDC2Data.cxx TSamuraiFDC2DataDict.cxx TSamuraiFDC2Physics.cxx TSamuraiFDC2PhysicsDict.cxx)
add_custom_command(OUTPUT TSamuraiFDC0DataDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TSamuraiFDC0Data.h TSamuraiFDC0DataDict.cxx TSamuraiFDC0Data.rootmap libNPSamurai.dylib DEPENDS TSamuraiFDC0Data.h)
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 TSamuraiFDC2Data.cxx TSamuraiFDC2DataDict.cxx TSamuraiFDC2Physics.cxx TSamuraiFDC2PhysicsDict.cxx TSamuraiFDC0Data.cxx TSamuraiFDC0DataDict.cxx TSamuraiFDC0Physics.cxx TSamuraiFDC0PhysicsDict.cxx)
target_link_libraries(NPSamurai ${ROOT_LIBRARIES} NPCore NPTrackReconstruction)
install(FILES TSamuraiFDC2Data.h TSamuraiFDC2Physics.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
install(FILES TSamuraiFDC2Data.h TSamuraiFDC2Physics.h TSamuraiFDC0Data.h TSamuraiFDC0Physics.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
#ifndef SAMURAIDCIndex_H
#define SAMURAIDCIndex_H
/*****************************************************************************
* 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 : October 2020 *
* Last update : *
*---------------------------------------------------------------------------*
* Decription: *
* This class hold SamuraiFDC0 treated data *
* *
*---------------------------------------------------------------------------*
* Comment: *
* *
* little class to index each of the DC wire *
* *
*****************************************************************************/
class SamuraiDCIndex{
public:
SamuraiDCIndex(){};
~SamuraiDCIndex(){};
SamuraiDCIndex(unsigned int det, unsigned int layer, unsigned int wire){
m_det=det;
m_layer=layer;
m_wire=wire;
m_norme=Norme();
};
unsigned int m_det;
unsigned int m_layer;
unsigned int m_wire;
unsigned int m_norme;
inline int Norme() const {return (m_det*1000000000+m_layer*1000000+m_wire);} ;
bool operator<(const SamuraiDCIndex i2){
return this->Norme()<i2.Norme();
}
friend bool operator<(const SamuraiDCIndex i1,const SamuraiDCIndex i2){
return i1.Norme()<i2.Norme();
}
friend bool operator==(const SamuraiDCIndex i1,const SamuraiDCIndex i2){
return i1.Norme()==i2.Norme();
}
};
#endif
#include "TSamuraiFDC0Data.h"
#include <iostream>
TSamuraiFDC0Data::TSamuraiFDC0Data(){};
TSamuraiFDC0Data::~TSamuraiFDC0Data(){};
////////////////////////////////////////////////////////////////////////////////
void TSamuraiFDC0Data::SetData(const int& Det, const int& Layer, const int& Wire, const double& Time, const int& Edge){
fFDC0_DetectorNbr.push_back(Det);
fFDC0_LayerNbr.push_back(Layer);
fFDC0_WireNbr.push_back(Wire);
fFDC0_Time.push_back(Time);
fFDC0_Edge.push_back(Edge);
}
////////////////////////////////////////////////////////////////////////////////
void TSamuraiFDC0Data::Clear(){
fFDC0_DetectorNbr.clear();
fFDC0_LayerNbr.clear();
fFDC0_WireNbr.clear();
fFDC0_Time.clear();
fFDC0_Edge.clear();
}
////////////////////////////////////////////////////////////////////////////////
void TSamuraiFDC0Data::Print(){
using namespace std;
cout << " -- Event:" << endl;
cout << " - Multiplicity: " << Mult() << endl;
}
////////////////////////////////////////////////////////////////////////////////
unsigned int TSamuraiFDC0Data::MultLayer(unsigned int det , unsigned int layer, int edge){
unsigned int mult=0;
unsigned int size = fFDC0_DetectorNbr.size();
for(unsigned int i = 0 ; i< size ; i++ ){
if(fFDC0_DetectorNbr[i]==det)
if(fFDC0_LayerNbr[i]==layer)
if(fFDC0_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> TSamuraiFDC0Data::GetWire(unsigned int det , unsigned int layer){
std::vector<int> wires;
unsigned int size = fFDC0_DetectorNbr.size();
for(unsigned int i = 0 ; i< size ; i++ ){
if(fFDC0_DetectorNbr[i]==det)
if(fFDC0_LayerNbr[i]==layer)
wires.push_back(fFDC0_WireNbr[i]);
}
return wires;
}
ClassImp(TSamuraiFDC0Data);
#ifndef TSamuraiFDC0Data_H
#define TSamuraiFDC0Data_H
#include "TObject.h"
#include <vector>
class TSamuraiFDC0Data: public TObject{
public:
TSamuraiFDC0Data();
~TSamuraiFDC0Data();
private:
std::vector<int> fFDC0_DetectorNbr;
std::vector<int> fFDC0_LayerNbr;
std::vector<int> fFDC0_WireNbr;
std::vector<double> fFDC0_Time;
std::vector<int> fFDC0_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 fFDC0_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 fFDC0_DetectorNbr[i];};
int const GetLayerNbr(const unsigned int& i){return fFDC0_LayerNbr[i];};
int const GetWireNbr(const unsigned int& i){return fFDC0_WireNbr[i];};
double const GetTime(const unsigned int& i){return fFDC0_Time[i];};
int const GetEdge(const unsigned int& i){return fFDC0_Edge[i];};
ClassDef(TSamuraiFDC0Data,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 : October 2020 *
* Last update : *
*---------------------------------------------------------------------------*
* Decription: *
* This class hold SamuraiFDC0 treated data *
* *
*---------------------------------------------------------------------------*
* Comment: *
* *
*****************************************************************************/
#include "TSamuraiFDC0Physics.h"
// STL
#include <sstream>
#include <iostream>
#include <cmath>
#include <stdlib.h>
#include <limits>
// NPL
#include "RootInput.h"
#include "RootOutput.h"
#include "TAsciiFile.h"
#include "NPOptionManager.h"
#include "NPDetectorFactory.h"
#include "NPSystemOfUnits.h"
// ROOT
using namespace NPUNITS;
///////////////////////////////////////////////////////////////////////////
ClassImp(TSamuraiFDC0Physics)
///////////////////////////////////////////////////////////////////////////
TSamuraiFDC0Physics::TSamuraiFDC0Physics(){
m_EventData = new TSamuraiFDC0Data ;
m_PreTreatedData = new TSamuraiFDC0Data ;
m_EventPhysics = this ;
//m_Spectra = NULL;
ToTThreshold_L = 0;
ToTThreshold_H = 180;
DriftLowThreshold=0.1 ;
DriftUpThreshold=2.4;
PowerThreshold=5;
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC0Physics::BuildSimplePhysicalEvent(){
BuildPhysicalEvent();
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC0Physics::BuildPhysicalEvent(){
PreTreat();
// RemoveNoise();
// Map[detector & plane angle, vector of spatial information]
static map<std::pair<unsigned int,double>, vector<double> > X ;
static map<std::pair<unsigned int,double>, vector<double> > Z ;
static map<std::pair<unsigned int,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);
std::pair<unsigned int, double> p(det,Wire_Angle[idx]);
//cout << det << " " << layer << " " << wire << endl;
X[p].push_back(Wire_X[idx]);
Z[p].push_back(Wire_Z[idx]);
R[p].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<std::pair<unsigned int,double>, TVector3 > VX0 ;
static map<std::pair<unsigned int,double>, TVector3 > VX100 ;
static map<std::pair<unsigned int,double>, double > D ;// the minimum distance
VX0.clear();VX100.clear(),D.clear();
for(auto it = X.begin();it!=X.end();++it){
D[it->first]=m_reconstruction.BuildTrack2D(X[it->first],Z[it->first],R[it->first],X0,X100,a,b);
// 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)>1000)
PileUp++;
Mult+=X[it->first].size();
// Position at z=0
TVector3 P(X0,0,0);
P.RotateZ(it->first.second);
VX0[it->first]=P;
// Position at z=100
TVector3 P100= TVector3(X100,0,0);
P100.RotateZ(it->first.second);
VX100[it->first]=P100;
}
// Reconstruct the central position (z=0) for each detector
static map<unsigned int,vector<TVector3> > C ;
static map<unsigned int,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 && it1->first.first==it2->first.first){// different plane, same detector
//cout << "FDC0" << endl;
m_reconstruction.ResolvePlane(it1->second,it1->first.second,it2->second,it2->first.second,P);
// cout << "done " << D[it1->first] << " " << D[it2->first] << endl;
if(P.X()!=-10000 && D[it1->first]<PowerThreshold&& D[it2->first]<PowerThreshold){
C[it1->first.first].push_back(P);
// Mean pos are weighted based on the the sum of distance from track
// to hit obtained during the minimisation
W[it1->first.first].push_back(1./sqrt(D[it1->first]*D[it2->first]));
}
}
}
}
// Reconstruct the position at z=100 for each detector
static map<unsigned int,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 && it1->first.first==it2->first.first){// different plane, same detector
m_reconstruction.ResolvePlane(it1->second,it1->first.second,it2->second,it2->first.second,P);
if(P.X()!=-10000&& D[it1->first]<PowerThreshold && D[it2->first]<PowerThreshold)
C100[it1->first.first].push_back(P);
}
}
}
// Build the Reference position by averaging all possible pair
size = C[0].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[0][i].X()*W[0][i];
PosY+= C[0][i].Y()*W[0][i];
PosX100+= C100[0][i].X()*W[0][i];
PosY100+= C100[0][i].Y()*W[0][i];
norm+=W[0][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[0][i]*(C[0][i].X()-PosX)*(C[0][i].X()-PosX);
devY+=W[0][i]*(C[0][i].Y()-PosY)*(C[0][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 TSamuraiFDC0Physics::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="SamuraiFDC0/L" + NPL::itoa(layer);
DriftLength.push_back(2.5-Cal->ApplySigmoid(channel,etime));
}
}
}
return;
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC0Physics::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 TSamuraiFDC0Physics::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 TSamuraiFDC0Physics::ReadConfiguration(NPL::InputParser parser){
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("SAMURAIFDC0");
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 FDC0 (" << i+1 << ")" << endl;
string xmlpath = blocks[i]->GetString("XML");
NPL::XmlParser xml;
xml.LoadFile(xmlpath);
AddDC("SAMURAIFDC0",xml);
}
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC0Physics::AddDC(string name, NPL::XmlParser& xml){
std::vector<NPL::XML::block*> b = xml.GetAllBlocksWithName(name);
// FDC0 case
if(name=="SAMURAIFDC0"){
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 FDC0"<< endl;
exit(1);
}
SamuraiDCIndex idx(det,layer,wire);
Wire_X[idx]=X;
Wire_Z[idx]=Z;
Wire_Angle[idx]=T;
}
}
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC0Physics::InitSpectra(){
//m_Spectra = new TSamuraiFDC0Spectra(m_NumberOfDetector);
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC0Physics::FillSpectra(){
// m_Spectra -> FillRawSpectra(m_EventData);
// m_Spectra -> FillPreTreatedSpectra(m_PreTreatedData);
// m_Spectra -> FillPhysicsSpectra(m_EventPhysics);
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC0Physics::CheckSpectra(){
// m_Spectra->CheckSpectra();
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC0Physics::ClearSpectra(){
// To be done
}
///////////////////////////////////////////////////////////////////////////
map< string , TH1*> TSamuraiFDC0Physics::GetSpectra() {
/* if(m_Spectra)
return m_Spectra->GetMapHisto();
else{
map< string , TH1*> empty;
return empty;
}*/
map< string , TH1*> empty;
return empty;
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC0Physics::WriteSpectra(){
// m_Spectra->WriteSpectra();
}
///////////////////////////////////////////////////////////////////////////
void TSamuraiFDC0Physics::AddParameterToCalibrationManager(){