Skip to content
Snippets Groups Projects
Commit e5446d73 authored by Adrien Matta's avatar Adrien Matta :skull_crossbones:
Browse files

* Fixing SuperX3 placement by Spherical coordinate

parent a950a988
No related branches found
No related tags found
No related merge requests found
......@@ -19,22 +19,22 @@
* *
*****************************************************************************/
#include <iostream>
#include <sstream>
#include <dlfcn.h>
#include <fstream>
#include <iostream>
#include <limits>
#include <sys/stat.h>
#include <sstream>
#include <stdlib.h>
#include <dlfcn.h>
#include <sys/stat.h>
#include "NPOptionManager.h"
#include "RootInput.h"
#include "TAsciiFile.h"
#include "NPOptionManager.h"
using namespace std;
RootInput* RootInput::instance = 0;
////////////////////////////////////////////////////////////////////////////////
RootInput* RootInput::getInstance(std::string configFileName){
RootInput* RootInput::getInstance(std::string configFileName) {
// A new instance of RootInput is created if it does not exist:
if (instance == 0) {
instance = new RootInput(configFileName);
......@@ -45,7 +45,7 @@ RootInput* RootInput::getInstance(std::string configFileName){
}
////////////////////////////////////////////////////////////////////////////////
void RootInput::Destroy(){
void RootInput::Destroy() {
if (instance != 0) {
delete instance;
instance = 0;
......@@ -53,64 +53,64 @@ void RootInput::Destroy(){
}
////////////////////////////////////////////////////////////////////////////////
void RootInput::ReadOldStyleInputFile(NPL::InputParser& parser){
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("TTreeName");
pTreeName = blocks[0]->GetLines()[0];
vector<NPL::InputBlock*> files = parser.GetAllBlocksWithToken("RootFileName");
if(files.size()>0){
vector<string> lines=files[0]->GetLines();
unsigned int size = lines.size();
for(unsigned int i = 0 ; i < size ; i++){
pTreePath.push_back(lines[i]);
}
void RootInput::ReadOldStyleInputFile(NPL::InputParser& parser) {
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("TTreeName");
pTreeName = blocks[0]->GetLines()[0];
vector<NPL::InputBlock*> files = parser.GetAllBlocksWithToken("RootFileName");
if (files.size() > 0) {
vector<string> lines = files[0]->GetLines();
unsigned int size = lines.size();
for (unsigned int i = 0; i < size; i++) {
pTreePath.push_back(lines[i]);
}
}
}
////////////////////////////////////////////////////////////////////////////////
void RootInput::ReadInputFile(NPL::InputParser& parser){
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Tree");
pTreeName = blocks[0]->GetMainValue();
std::vector<std::string> lines=blocks[0]->GetLines();
unsigned int size = lines.size();
for(unsigned int i = 0 ; i < size ; i++){
if(lines[i].find(".root")!=string::npos)
pTreePath.push_back(lines[i]);
else if(lines[i].find(".tree")!=string::npos)
ReadTreeFile(lines[i]);
}
void RootInput::ReadInputFile(NPL::InputParser& parser) {
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Tree");
pTreeName = blocks[0]->GetMainValue();
std::vector<std::string> lines = blocks[0]->GetLines();
unsigned int size = lines.size();
for (unsigned int i = 0; i < size; i++) {
if (lines[i].find(".root") != string::npos)
pTreePath.push_back(lines[i]);
else if (lines[i].find(".tree") != string::npos)
ReadTreeFile(lines[i]);
}
vector<NPL::InputBlock*> friends = parser.GetAllBlocksWithToken("Friend");
unsigned int sizeF = friends.size();
for(unsigned int i = 0 ; i < sizeF ; i++){
pFriendsPath.insert(pair< string,vector<string> > (friends[i]->GetMainValue(),friends[i]->GetLines()));
}
vector<NPL::InputBlock*> friends = parser.GetAllBlocksWithToken("Friend");
unsigned int sizeF = friends.size();
for (unsigned int i = 0; i < sizeF; i++) {
pFriendsPath.insert(pair<string, vector<string>>(friends[i]->GetMainValue(), friends[i]->GetLines()));
}
}
////////////////////////////////////////////////////////////////////////////////
void RootInput::ReadTreeFile(std::string path){
void RootInput::ReadTreeFile(std::string path) {
ifstream tree(path.c_str());
path=path.substr(0,path.rfind("/")+1);
path = path.substr(0, path.rfind("/") + 1);
std::string buffer;
bool first=true;
unsigned int count = 0 ;
while(tree>>buffer){
if(first){
pTreePath.push_back(path+buffer);
bool first = true;
unsigned int count = 0;
while (tree >> buffer) {
if (first) {
pTreePath.push_back(path + buffer);
count++;
first=false;
first = false;
}
else{
pFriendsTreePath[count++].push_back(path+buffer);
else {
pFriendsTreePath[count++].push_back(path + buffer);
}
}
}
////////////////////////////////////////////////////////////////////////////////
RootInput::RootInput(std::string configFileName){
RootInput::RootInput(std::string configFileName) {
NumberOfFriend = 0;
pRootFile = NULL;
std::string lastfile= NPOptionManager::getInstance()->GetLastFile();
if(lastfile!="VOID"){
pRootFile = NULL;
std::string lastfile = NPOptionManager::getInstance()->GetLastFile();
if (lastfile != "VOID") {
configFileName = lastfile;
}
......@@ -120,25 +120,25 @@ RootInput::RootInput(std::string configFileName){
// Old style file
vector<NPL::InputBlock*> old_blocks = parser.GetAllBlocksWithToken("TTreeName");
if(old_blocks.size()==1){
if (old_blocks.size() == 1) {
ReadOldStyleInputFile(parser);
}
// New style file
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Tree");
if(blocks.size()==1){
if (blocks.size() == 1) {
ReadInputFile(parser);
}
// If the tree come from a simulation, the InteractionCoordinates
// and InitialConditions lib are loaded
if(pTreeName=="SimulatedTree"){
std::string path = getenv("NPTOOL");
path+="/NPLib/lib/";
std::string libName="libNPInteractionCoordinates"+NPOptionManager::getInstance()->GetSharedLibExtension();
libName=path+libName;
dlopen(libName.c_str(),RTLD_NOW);
libName="libNPInitialConditions"+NPOptionManager::getInstance()->GetSharedLibExtension();
libName=path+libName;
dlopen(libName.c_str(),RTLD_NOW);
if (pTreeName == "SimulatedTree") {
std::string path = getenv("NPTOOL");
path += "/NPLib/lib/";
std::string libName = "libNPInteractionCoordinates" + NPOptionManager::getInstance()->GetSharedLibExtension();
libName = path + libName;
dlopen(libName.c_str(), RTLD_NOW);
libName = "libNPInitialConditions" + NPOptionManager::getInstance()->GetSharedLibExtension();
libName = path + libName;
dlopen(libName.c_str(), RTLD_NOW);
}
// Initialise the chain
......@@ -147,106 +147,111 @@ RootInput::RootInput(std::string configFileName){
// Add all the files
unsigned int size = pTreePath.size();
std::string firstfile;
for(unsigned int i = 0 ; i < size ; i++){
for (unsigned int i = 0; i < size; i++) {
cout << " - Adding file : " << pTreePath[i].c_str() << endl;
pRootChain->Add(pTreePath[i].c_str());
// Test if the file is a regex or a single file
double counts;
std::string command = "ls " + pTreePath[i] + " > .ls_return";
counts= system(command.c_str());
counts = system(command.c_str());
std::ifstream return_ls(".ls_return");
std::string files;
while(return_ls >> files){
if(counts == 0)
while (return_ls >> files) {
if (counts == 0)
firstfile = files;
counts++;
}
}
// Case of user made Friends
for(auto it = pFriendsPath.begin(); it!=pFriendsPath.end() ; it++){
for (auto it = pFriendsPath.begin(); it != pFriendsPath.end(); it++) {
TChain* chain = new TChain(it->first.c_str());
cout << " - Adding friend : " << endl;
for(auto itp = it->second.begin() ; itp!=it->second.end() ; itp++){
cout << " - " << (*itp).c_str() << endl;
chain->Add((*itp).c_str());
for (auto itp = it->second.begin(); itp != it->second.end(); itp++) {
cout << " - " << (*itp).c_str() << endl;
chain->Add((*itp).c_str());
}
pRootChain->AddFriend(chain);
}
// Case of tree file
for(auto it = pFriendsTreePath.begin(); it!=pFriendsTreePath.end() ; it++){
for (auto it = pFriendsTreePath.begin(); it != pFriendsTreePath.end(); it++) {
TChain* chain = new TChain(pTreeName.c_str());
cout << " - Adding friend : " << endl;
for(auto itp = it->second.begin() ; itp!=it->second.end() ; itp++){
for (auto itp = it->second.begin(); itp != it->second.end(); itp++) {
cout << " - " << (*itp).c_str() << endl;
chain->Add((*itp).c_str());
}
}
pRootChain->AddFriend(chain);
}
if (!pRootFile)
if (!pRootFile)
pRootFile = new TFile(firstfile.c_str());
if( pRootChain->GetEntries() ==0){
std::cout << "\033[1;31m**** ERROR: No entries to analyse ****\033[0m" << std::endl;
if (pRootChain->GetEntries() == 0) {
std::cout << "\033[1;31m**** ERROR: No entries to analyse ****\033[0m" << std::endl;
exit(1);
}
else
std::cout << "\033[1;32mROOTInput: " << pRootChain->GetEntries() << " entries loaded in the input chain\033[0m" << std::endl ;
else
std::cout << "\033[1;32mROOTInput: " << pRootChain->GetEntries() << " entries loaded in the input chain\033[0m"
<< std::endl;
}
////////////////////////////////////////////////////////////////////////////////
std::string RootInput::DumpAsciiFile(const char* type, const char* folder){
std::string name="fail";
std::string RootInput::DumpAsciiFile(const char* type, const char* folder) {
std::string name = "fail";
std::string sfolder = folder;
// create folder if not existing
// first test if folder exist
struct stat dirInfo;
int res = stat(folder, &dirInfo);
if (res != 0) { // if folder does not exist, create it
if (res != 0) { // if folder does not exist, create it
std::string cmd = "mkdir " + sfolder;
if (system(cmd.c_str()) != 0) std::cout << "RootInput::DumpAsciiFile() problem creating directory" << std::endl;
if (system(cmd.c_str()) != 0)
std::cout << "RootInput::DumpAsciiFile() problem creating directory" << std::endl;
}
std::string stype = type;
if (stype == "EventGenerator") {
TAsciiFile *aFile = (TAsciiFile*)pRootFile->Get(stype.c_str());
TAsciiFile* aFile = (TAsciiFile*)pRootFile->Get(stype.c_str());
if(aFile)
{
if (aFile) {
// build file name
std::string title = aFile->GetTitle();
size_t pos = title.rfind("/");
if (pos != std::string::npos) name = sfolder + title.substr(pos);
else name = sfolder + "/" + title;
if (pos != std::string::npos)
name = sfolder + title.substr(pos);
else
name = sfolder + "/" + title;
aFile->WriteToFile(name.c_str());
}
}
else if (stype == "DetectorConfiguration") {
TAsciiFile *aFile = (TAsciiFile*)pRootFile->Get(stype.c_str());
if(aFile)
{
TAsciiFile* aFile = (TAsciiFile*)pRootFile->Get(stype.c_str());
if (aFile) {
// build file name
std::string title = aFile->GetTitle();
size_t pos = title.rfind("/");
if (pos != std::string::npos) name = sfolder + title.substr(pos);
else name = sfolder + "/" + title;
if (pos != std::string::npos)
name = sfolder + title.substr(pos);
else
name = sfolder + "/" + title;
aFile->WriteToFile(name.c_str());
}
}
else if (stype == "Calibration") {
TAsciiFile *aFile = (TAsciiFile*)pRootFile->Get(stype.c_str());
if(aFile)
{
TAsciiFile* aFile = (TAsciiFile*)pRootFile->Get(stype.c_str());
if (aFile) {
// build file name
std::string title = aFile->GetTitle();
size_t pos = title.rfind("/");
if (pos != std::string::npos) name = sfolder + title.substr(pos);
else name = sfolder + "/" + title;
if (pos != std::string::npos)
name = sfolder + title.substr(pos);
else
name = sfolder + "/" + title;
aFile->WriteToFile(name.c_str());
}
}
......@@ -261,12 +266,13 @@ std::string RootInput::DumpAsciiFile(const char* type, const char* folder){
}
////////////////////////////////////////////////////////////////////////////////
RootInput::~RootInput(){
RootInput::~RootInput() {
// delete default directory ./.tmp
struct stat dirInfo;
int res = stat("./.tmp", &dirInfo);
if (res == 0) { // if does exist, delete it
if (system("rm -rf ./.tmp") != 0) std::cout << "RootInput::~RootInput() problem deleting ./.tmp directory" << std::endl;
if (res == 0) { // if does exist, delete it
if (system("rm -rf ./.tmp") != 0)
std::cout << "RootInput::~RootInput() problem deleting ./.tmp directory" << std::endl;
}
std::cout << std::endl << "Root Input summary" << std::endl;
std::cout << " - Number of bites read: " << pRootFile->GetBytesRead() << std::endl;
......
......@@ -6,7 +6,7 @@ add_custom_command(OUTPUT TMDMDataDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/b
##
## Check if MINUIT2 is installed along with ROOT
find_library(libMinuit2_FOUND NAMES Minuit2 HINTS "${ROOTSYS}/lib")
find_library(libMinuit2_FOUND NAMES Minuit2 HINTS "${ROOTSYS}/lib" "/usr/lib64/root")
if(libMinuit2_FOUND)
message(STATUS "Minuit2 support enabled for MDM.")
add_definitions(-DHAVE_MINUIT2)
......
......@@ -10,7 +10,7 @@ add_custom_command(OUTPUT TMinosClustDict.cxx COMMAND ../../scripts/build_dict.s
add_custom_command(OUTPUT TMinosResultDict.cxx COMMAND ../../scripts/build_dict.sh TMinosResult.h TMinosResultDict.cxx TMinosResult.rootmap libNPMinos.dylib DEPENDS TMinosResult.h)
## Check if MINUIT is installed along with ROOT
find_library(libMinuit_FOUND NAMES Minuit2 HINTS "$ENV{ROOTSYS}/lib")
find_library(libMinuit_FOUND NAMES Minuit2 HINTS "$ENV{ROOTSYS}/lib" "/usr/lib64/root")
if(libMinuit_FOUND)
message(STATUS "Minuit support enabled for Minos.")
add_definitions(-DHAVE_MINUIT)
......
......@@ -10,33 +10,33 @@
* *
* Creation Date : 2014/11/24 *
* Last update : 2019/09 implemeted in NPTool by Cyril Lenain *
* lenain@lpccaen.in2p3.fr *
* lenain@lpccaen.in2p3.fr *
*---------------------------------------------------------------------------*
* Decription: *
* This class hold Minos Treated data *
* *
*---------------------------------------------------------------------------*
* Comment: *
* *
* *
* *
*****************************************************************************/
#include "TMinosPhysics.h"
// STL
#include <sstream>
#include <iostream>
#include <cmath>
#include <stdlib.h>
#include <iostream>
#include <limits>
#include <sstream>
#include <stdlib.h>
using namespace std;
// NPL
#include "RootInput.h"
#include "RootOutput.h"
#include "NPDetectorFactory.h"
#include "NPOptionManager.h"
#include "NPSystemOfUnits.h"
#include "RootInput.h"
#include "RootOutput.h"
// ROOT
#include "TChain.h"
......@@ -44,100 +44,94 @@ using namespace std;
TMinosPhysics* current_phy = 0;
///////////////////////////////////////////////////////////////////////////
TMinosPhysics::TMinosPhysics(){
m_EventData=new TMinosData;
m_EventPhysics=this;
m_E_RAW_Threshold=0; // adc channels
m_E_Threshold=0; // MeV
m_ransac.SetParameters(100,// max 100 iteration
5, // a point belong to a track if it is within 5 mm
40, // max distance allowed between a pair of point to consider a track
10);// minimum point to form a valid cluster
TMinosPhysics::TMinosPhysics() {
m_EventData = new TMinosData;
m_EventPhysics = this;
m_E_RAW_Threshold = 0; // adc channels
m_E_Threshold = 0; // MeV
m_ransac.SetParameters(100, // max 100 iteration
5, // a point belong to a track if it is within 5 mm
40, // max distance allowed between a pair of point to consider a track
10); // minimum point to form a valid cluster
}
///////////////////////////////////////////////////////////////////////////
void TMinosPhysics::BuildSimplePhysicalEvent() {
BuildPhysicalEvent();
}
void TMinosPhysics::BuildSimplePhysicalEvent() { BuildPhysicalEvent(); }
///////////////////////////////////////////////////////////////////////////
void TMinosPhysics::BuildPhysicalEvent() {
PreTreat();
vector<NPL::LinearCluster3D> clusters = m_ransac.TreatEvent(X_Pad,Y_Pad,Z_Pad);
vector<NPL::LinearCluster3D> clusters = m_ransac.TreatEvent(X_Pad, Y_Pad, Z_Pad);
unsigned int sizeC = clusters.size();
for(unsigned int i = 0 ; i < sizeC ; i++){
// Extract the Direction of the track
for (unsigned int i = 0; i < sizeC; i++) {
// Extract the Direction of the track
clusters[i].LinearFit();
Tracks_P0.push_back(clusters[i].GetP0());
Tracks_Dir.push_back(clusters[i].GetDir());
}
if(sizeC==2){
static TVector3 Vertex,delta;
Delta_Vertex =
MinimumDistanceTwoLines(Tracks_P0[0],Tracks_P0[0]+Tracks_Dir[0],
Tracks_P0[1],Tracks_P0[1]+Tracks_Dir[1],
Vertex, delta);
X_Vertex= Vertex.X();
Y_Vertex= Vertex.Y();
Z_Vertex= Vertex.Z();
Theta_12=180.*Tracks_Dir[0].Angle(Tracks_Dir[1])/(M_PI);
if (sizeC == 2) {
static TVector3 Vertex, delta;
Delta_Vertex = MinimumDistanceTwoLines(Tracks_P0[0], Tracks_P0[0] + Tracks_Dir[0], Tracks_P0[1],
Tracks_P0[1] + Tracks_Dir[1], Vertex, delta);
X_Vertex = Vertex.X();
Y_Vertex = Vertex.Y();
Z_Vertex = Vertex.Z();
Theta_12 = 180. * Tracks_Dir[0].Angle(Tracks_Dir[1]) / (M_PI);
}
}
///////////////////////////////////////////////////////////////////////////
void TMinosPhysics::PreTreat() {
// apply thresholds and calibration
static unsigned int sizePad,sizeQ ;
static unsigned int sizePad, sizeQ;
static unsigned short PadNumber;
static double Q,T;
static auto cal= CalibrationManager::getInstance();
static double Q, T;
static auto cal = CalibrationManager::getInstance();
static string cal_v, cal_o;
sizePad = m_EventData->GetPadMult();
if(sizePad>20){
for(unsigned int i = 0 ; i < sizePad ; i++){
if (sizePad > 20) {
for (unsigned int i = 0; i < sizePad; i++) {
vector<unsigned short>* Charge = m_EventData->GetChargePtr(i);
if(Charge->size()>40 ){
vector<unsigned short>* Time = m_EventData->GetTimePtr(i);
m_utility.Calibrate(Time,Charge,i,T,Q);
if (Charge->size() > 40) {
vector<unsigned short>* Time = m_EventData->GetTimePtr(i);
m_utility.Calibrate(Time, Charge, i, T, Q);
if(T>0&&Q<77000){
if (T > 0 && Q < 77000) {
PadNumber = m_EventData->GetPadNumber(i);
double x_mm = m_X[PadNumber];
double y_mm = m_Y[PadNumber];
unsigned int ring = round((sqrt(x_mm*x_mm + y_mm*y_mm)-44.15)/2.1);
cal_v="Minos/R"+NPL::itoa(ring)+"_VDRIFT";
cal_o="Minos/R"+NPL::itoa(ring)+"_OFFSET";
double calV= cal->GetValue(cal_v,0);
double calO=cal->GetValue(cal_o,0);
double z_mm = (T*m_TimeBin+calO/*-0.7-0.04*ring*/)*calV;
//cout << T*m_TimeBin+calO << endl;
TVector3 Pos=TVector3(x_mm+m_Position.X(),y_mm+m_Position.Y(),z_mm+m_Position.Z());
Pos.RotateZ(m_ZRotation);
unsigned int ring = round((sqrt(x_mm * x_mm + y_mm * y_mm) - 44.15) / 2.1);
cal_v = "Minos/R" + NPL::itoa(ring) + "_VDRIFT";
cal_o = "Minos/R" + NPL::itoa(ring) + "_OFFSET";
double calV = cal->GetValue(cal_v, 0);
double calO = cal->GetValue(cal_o, 0);
double z_mm = (T * m_TimeBin + calO /*-0.7-0.04*ring*/) * calV;
// cout << T*m_TimeBin+calO << endl;
TVector3 Pos = TVector3(x_mm + m_Position.X(), y_mm + m_Position.Y(), z_mm + m_Position.Z());
Pos.RotateZ(m_ZRotation);
// Calibrate the Pad:
X_Pad.push_back(Pos.X());
Y_Pad.push_back(Pos.Y());
Z_Pad.push_back(Pos.Z());
Z_Pad.push_back(Pos.Z());
Ring_Pad.push_back(ring);
Q_Pad.push_back(Q);
T_Pad.push_back(T);
Q_Pad.push_back(Q);
T_Pad.push_back(T);
}
}
else{
else {
/* X_Pad.push_back(-1);
Y_Pad.push_back(-1);
Z_Pad.push_back((-1);
Ring_Pad.push_back(-1);
Z_Pad.push_back(-1);
Q_Pad.push_back(-1);
T_Pad.push_back(-1);*/
Z_Pad.push_back(-1);
Q_Pad.push_back(-1);
T_Pad.push_back(-1);*/
}
}
}
}
///////////////////////////////////////////////////////////////////////////
......@@ -162,33 +156,33 @@ void TMinosPhysics::ReadAnalysisConfig() {
asciiConfig->Append(FileName.c_str());
asciiConfig->AppendLine("");
// read analysis config file
string LineBuffer,DataBuffer,whatToDo;
string LineBuffer, DataBuffer, whatToDo;
while (!AnalysisConfigFile.eof()) {
// Pick-up next line
getline(AnalysisConfigFile, LineBuffer);
// search for "header"
string name = "ConfigMinos";
if (LineBuffer.compare(0, name.length(), name) == 0)
if (LineBuffer.compare(0, name.length(), name) == 0)
ReadingStatus = true;
// loop on tokens and data
while (ReadingStatus ) {
whatToDo="";
while (ReadingStatus) {
whatToDo = "";
AnalysisConfigFile >> whatToDo;
// Search for comment symbol (%)
if (whatToDo.compare(0, 1, "%") == 0) {
AnalysisConfigFile.ignore(numeric_limits<streamsize>::max(), '\n' );
AnalysisConfigFile.ignore(numeric_limits<streamsize>::max(), '\n');
}
else if (whatToDo=="E_RAW_THRESHOLD") {
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") {
else if (whatToDo == "E_THRESHOLD") {
AnalysisConfigFile >> DataBuffer;
m_E_Threshold = atof(DataBuffer.c_str());
cout << whatToDo << " " << m_E_Threshold << endl;
......@@ -199,14 +193,13 @@ void TMinosPhysics::ReadAnalysisConfig() {
}
}
}
}
///////////////////////////////////////////////////////////////////////////
void TMinosPhysics::Clear() {
Ring_Pad.clear();
X_Pad.clear();
Y_Pad.clear();
Ring_Pad.clear();
X_Pad.clear();
Y_Pad.clear();
Z_Pad.clear();
Q_Pad.clear();
T_Pad.clear();
......@@ -216,91 +209,86 @@ void TMinosPhysics::Clear() {
Tracks_Dir.clear();
// Vertex information
X_Vertex=-1000;
Y_Vertex=-1000;
Z_Vertex=-1000;
Theta_12=-1000;
Delta_Vertex=-1000;
X_Vertex = -1000;
Y_Vertex = -1000;
Z_Vertex = -1000;
Theta_12 = -1000;
Delta_Vertex = -1000;
}
///////////////////////////////////////////////////////////////////////////
void TMinosPhysics::ReadConfiguration(NPL::InputParser parser) {
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Minos");
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << "//// " << blocks.size() << " detector(s) found " << endl;
if (NPOptionManager::getInstance()->GetVerboseLevel())
cout << "//// " << blocks.size() << " detector(s) found " << endl;
vector<string> token= {"XML","TimeBin","ShapingTime","Baseline","Sampling","Position","ZRotation"};
vector<string> token = {"XML", "TimeBin", "ShapingTime", "Baseline", "Sampling", "Position", "ZRotation"};
for (unsigned int i = 0; i < blocks.size(); i++) {
for(unsigned int i = 0 ; i < blocks.size() ; i++){
if (blocks[i]->HasTokenList(token)) {
cout << endl << "//// MINOS (" << i+1 << ")" << endl;
m_ShapingTime = blocks[i]->GetDouble("ShapingTime","ns")/NPUNITS::ns;
m_TimeBin = blocks[i]->GetDouble("TimeBin","ns")/NPUNITS::ns;
m_Sampling= blocks[i]->GetInt("Sampling");
m_Baseline= blocks[i]->GetInt("BaseLine");
m_utility.SetParameters(m_TimeBin,m_ShapingTime,m_Baseline,m_Sampling);
m_Position = blocks[i]->GetTVector3("Position","mm");
m_ZRotation= blocks[i]->GetDouble("ZRotation","deg");
cout << endl << "//// MINOS (" << i + 1 << ")" << endl;
m_ShapingTime = blocks[i]->GetDouble("ShapingTime", "ns") / NPUNITS::ns;
m_TimeBin = blocks[i]->GetDouble("TimeBin", "ns") / NPUNITS::ns;
m_Sampling = blocks[i]->GetInt("Sampling");
m_Baseline = blocks[i]->GetInt("BaseLine");
m_utility.SetParameters(m_TimeBin, m_ShapingTime, m_Baseline, m_Sampling);
m_Position = blocks[i]->GetTVector3("Position", "mm");
m_ZRotation = blocks[i]->GetDouble("ZRotation", "deg");
string xmlpath = blocks[i]->GetString("XML");
NPL::XmlParser xml;
xml.LoadFile(xmlpath);
ReadXML(xml);
}
else {
cout << "ERROR: Missing token for Minos, check your input file"<< endl;
cout << "Required token: ";
for(unsigned int i = 0 ; i < token.size() ; i++)
cout << token[i] << " " ;
cout << "ERROR: Missing token for Minos, check your input file" << endl;
cout << "Required token: ";
for (unsigned int i = 0; i < token.size(); i++)
cout << token[i] << " ";
cout << endl;
exit(1);
}
}
}
///////////////////////////////////////////////////////////////////////////
void TMinosPhysics::ReadXML(NPL::XmlParser& xml){
std::vector<NPL::XML::block*> b = xml.GetAllBlocksWithName("MINOS");
void TMinosPhysics::ReadXML(NPL::XmlParser& xml) {
std::vector<NPL::XML::block*> b = xml.GetAllBlocksWithName("MINOS");
unsigned int size = b.size();
unsigned int valid_pad=0;
for(unsigned int i = 0 ; i < size ; i++){
unsigned int valid_pad = 0;
for (unsigned int i = 0; i < size; i++) {
unsigned short ID = b[i]->AsInt("ID");
m_X[ID] = b[i]->AsDouble("X");
m_Y[ID] = b[i]->AsDouble("Y");
unsigned short ID = b[i]->AsInt("ID");
m_X[ID] = b[i]->AsDouble("X");
m_Y[ID] = b[i]->AsDouble("Y");
valid_pad++;
if((m_X[ID]==0&&m_Y[ID]==0)||(m_X[ID]==-1&&m_Y[ID]==-1))
if ((m_X[ID] == 0 && m_Y[ID] == 0) || (m_X[ID] == -1 && m_Y[ID] == -1))
valid_pad--;
m_Period[ID] = b[i]->AsDouble("TPERIOD");
m_Offset[ID] = b[i]->AsDouble("TOFFSET");
m_Gain[ID] = b[i]->AsDouble("GAIN");
m_Period[ID] = b[i]->AsDouble("TPERIOD");
m_Offset[ID] = b[i]->AsDouble("TOFFSET");
m_Gain[ID] = b[i]->AsDouble("GAIN");
}
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << "//// " << m_X.size() << " pads found with " << valid_pad << " valid pads" << endl;
if (NPOptionManager::getInstance()->GetVerboseLevel())
cout << "//// " << m_X.size() << " pads found with " << valid_pad << " valid pads" << endl;
}
///////////////////////////////////////////////////////////////////////////
void TMinosPhysics::AddParameterToCalibrationManager() {
CalibrationManager* Cal = CalibrationManager::getInstance();
unsigned int NbrRing = 18;
vector<double> standardV={0.03475};
vector<double> standardO={0};
for(unsigned int i = 0 ; i < NbrRing ; i++){
Cal->AddParameter("Minos", "R"+NPL::itoa(i+1)+"_VDRIFT","Minos_R"+NPL::itoa(i+1)+"_VDRIFT",standardV);
Cal->AddParameter("Minos", "R"+NPL::itoa(i+1)+"_OFFSET","Minos_R"+NPL::itoa(i+1)+"_OFFSET",standardO);
vector<double> standardV = {0.03475};
vector<double> standardO = {0};
for (unsigned int i = 0; i < NbrRing; i++) {
Cal->AddParameter("Minos", "R" + NPL::itoa(i + 1) + "_VDRIFT", "Minos_R" + NPL::itoa(i + 1) + "_VDRIFT", standardV);
Cal->AddParameter("Minos", "R" + NPL::itoa(i + 1) + "_OFFSET", "Minos_R" + NPL::itoa(i + 1) + "_OFFSET", standardO);
}
}
///////////////////////////////////////////////////////////////////////////
void TMinosPhysics::InitializeRootInputRaw() {
TChain* inputChain = RootInput::getInstance()->GetChain();
inputChain->SetBranchStatus("Minos", true );
inputChain->SetBranchAddress("Minos", &m_EventData );
inputChain->SetBranchStatus("Minos", true);
inputChain->SetBranchAddress("Minos", &m_EventData);
}
///////////////////////////////////////////////////////////////////////////
......@@ -318,21 +306,19 @@ void TMinosPhysics::InitializeRootOutput() {
////////////////////////////////////////////////////////////////////////////////
// Construct Method to be pass to the DetectorFactory //
////////////////////////////////////////////////////////////////////////////////
NPL::VDetector* TMinosPhysics::Construct() {
return (NPL::VDetector*) new TMinosPhysics();
}
NPL::VDetector* TMinosPhysics::Construct() { return (NPL::VDetector*)new TMinosPhysics(); }
////////////////////////////////////////////////////////////////////////////////
// Registering the construct method to the factory //
////////////////////////////////////////////////////////////////////////////////
extern "C"{
class proxy_Minos{
public:
proxy_Minos(){
NPL::DetectorFactory::getInstance()->AddToken("Minos","Minos");
NPL::DetectorFactory::getInstance()->AddDetector("Minos",TMinosPhysics::Construct);
}
};
extern "C" {
class proxy_Minos {
public:
proxy_Minos() {
NPL::DetectorFactory::getInstance()->AddToken("Minos", "Minos");
NPL::DetectorFactory::getInstance()->AddDetector("Minos", TMinosPhysics::Construct);
}
};
proxy_Minos p_Minos;
proxy_Minos p_Minos;
}
......@@ -7,7 +7,7 @@ add_custom_command(OUTPUT NPClusterDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/
add_custom_command(OUTPUT TrackingDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh Tracking.h TrackingDict.cxx Tracking.rootmap libNPTrackReconstruction.so NPTrackReconstructionLinkDef.h DEPENDS Tracking.h)
## Check if MINUIT2 is installed along with ROOT
find_library(libMinuit2_FOUND NAMES Minuit2 HINTS "${ROOTSYS}/lib")
find_library(libMinuit2_FOUND NAMES Minuit2 HINTS "${ROOTSYS}/lib" "/usr/lib64/root")
if(libMinuit2_FOUND)
message(STATUS "Minuit2 support enabled for TrackReconstruction.")
add_definitions(-DHAVE_MINUIT2)
......
......@@ -16,139 +16,143 @@
* Find linear tracks using the ransac method *
*****************************************************************************/
#include"NPLinearRansac3D.h"
#include "NPLinearRansac3D.h"
using namespace std;
using namespace NPL;
////////////////////////////////////////////////////////////////////////////////
void LinearCluster3D::LinearFit(){
void LinearCluster3D::LinearFit() {
unsigned int sizeI = m_index.size();
unsigned int sizeI = m_index.size();
double min_distance = 1e20;
unsigned int min_i = 0 ; unsigned int min_j = 0;
double meanDX=0;
double meanDY=0;
double meanDZ=0;
double totalW=0;
double meanX=0;
double meanY=0;
double meanZ=0;
unsigned int min_i = 0;
unsigned int min_j = 0;
double meanDX = 0;
double meanDY = 0;
double meanDZ = 0;
double totalW = 0;
double meanX = 0;
double meanY = 0;
double meanZ = 0;
double w;
for(unsigned int i = 0 ; i < sizeI ; i++){
TVector3 Pi((*m_X)[m_index[i]],(*m_Y)[m_index[i]],(*m_Z)[m_index[i]]);
for(unsigned int j = i+1 ; j < sizeI ; j++){
double dij=0;
TVector3 Pj((*m_X)[m_index[j]],(*m_Y)[m_index[j]],(*m_Z)[m_index[j]]);
for (unsigned int i = 0; i < sizeI; i++) {
TVector3 Pi((*m_X)[m_index[i]], (*m_Y)[m_index[i]], (*m_Z)[m_index[i]]);
for (unsigned int j = i + 1; j < sizeI; j++) {
double dij = 0;
TVector3 Pj((*m_X)[m_index[j]], (*m_Y)[m_index[j]], (*m_Z)[m_index[j]]);
// compute the average distance of all point to this line
for(unsigned int p = 0 ; p < sizeI ; p++){
TVector3 Pp((*m_X)[m_index[p]],(*m_Y)[m_index[p]],(*m_Z)[m_index[p]]);
dij+=MinimumDistancePointLine(Pi,Pj,Pp);
for (unsigned int p = 0; p < sizeI; p++) {
TVector3 Pp((*m_X)[m_index[p]], (*m_Y)[m_index[p]], (*m_Z)[m_index[p]]);
dij += MinimumDistancePointLine(Pi, Pj, Pp);
}
w = 1./dij;
totalW+=w;
meanX+=w*((*m_X)[m_index[i]]);
meanY+=w*((*m_Y)[m_index[i]]);
meanZ+=w*((*m_Z)[m_index[i]]);
meanDX+=w*((*m_X)[m_index[i]]-(*m_X)[m_index[j]]);
meanDY+=w*((*m_Y)[m_index[i]]-(*m_Y)[m_index[j]]);
meanDZ+=w*((*m_Z)[m_index[i]]-(*m_Z)[m_index[j]]);
}
}
meanDX/=totalW;
meanDY/=totalW;
meanDZ/=totalW;
meanX/=totalW;
meanY/=totalW;
meanZ/=totalW;
w = 1. / dij;
totalW += w;
m_P0 = TVector3(meanX,meanY,meanZ);
m_Dir = (TVector3(meanDX,meanDY,meanDZ)).Unit();
meanX += w * ((*m_X)[m_index[i]]);
meanY += w * ((*m_Y)[m_index[i]]);
meanZ += w * ((*m_Z)[m_index[i]]);
meanDX += w * ((*m_X)[m_index[i]] - (*m_X)[m_index[j]]);
meanDY += w * ((*m_Y)[m_index[i]] - (*m_Y)[m_index[j]]);
meanDZ += w * ((*m_Z)[m_index[i]] - (*m_Z)[m_index[j]]);
}
}
meanDX /= totalW;
meanDY /= totalW;
meanDZ /= totalW;
meanX /= totalW;
meanY /= totalW;
meanZ /= totalW;
m_P0 = TVector3(meanX, meanY, meanZ);
m_Dir = (TVector3(meanDX, meanDY, meanDZ)).Unit();
// Cluster dir is arbitrary. we choose to always have positive Z dir
if (m_Dir.Z() < 0)
m_Dir = -m_Dir;
};
////////////////////////////////////////////////////////////////////////////////
vector<LinearCluster3D> LinearRansac3D::TreatEvent(vector<double>& X, vector<double>&Y, vector<double>&Z){
cluster_id.clear();
clusters.clear();
unsigned int sizeX = X.size();
cluster_id.resize(sizeX,0);
m_cluster.clear();
m_assigned.clear();
if(sizeX<m_min_count)
return clusters;
unsigned int p1,p2;
double d;
TVector3 Vp1,Vp2,D,P;
m_iteration_d=0;
m_iteration=0;
while(m_iteration++ < m_max_iteration && m_assigned.size()<sizeX){
LinearCluster3D cluster(&X,&Y,&Z);
// take 2 distant point randomly that has not been match before
d=0 ; m_iteration_d=0;
while(d<3*m_match_distance && m_iteration_d++<m_max_iteration ){
p1 = rand.Integer(sizeX);
p2 = rand.Integer(sizeX);
Vp1.SetXYZ(X[p1],Y[p1],Z[p1]);
Vp2.SetXYZ(X[p2],Y[p2],Z[p2]);
D=Vp1-Vp2;
d = D.Mag();
if(d>m_max_distance)
d=0;
}
vector<LinearCluster3D> LinearRansac3D::TreatEvent(vector<double>& X, vector<double>& Y, vector<double>& Z) {
cluster_id.clear();
clusters.clear();
unsigned int sizeX = X.size();
cluster_id.resize(sizeX, 0);
m_cluster.clear();
m_assigned.clear();
if (sizeX < m_min_count)
return clusters;
unsigned int p1, p2;
double d;
TVector3 Vp1, Vp2, D, P;
m_iteration_d = 0;
m_iteration = 0;
while (m_iteration++ < m_max_iteration && m_assigned.size() < sizeX) {
LinearCluster3D cluster(&X, &Y, &Z);
// take 2 distant point randomly that has not been match before
d = 0;
m_iteration_d = 0;
while (d < 3 * m_match_distance && m_iteration_d++ < m_max_iteration) {
p1 = rand.Integer(sizeX);
p2 = rand.Integer(sizeX);
Vp1.SetXYZ(X[p1], Y[p1], Z[p1]);
Vp2.SetXYZ(X[p2], Y[p2], Z[p2]);
D = Vp1 - Vp2;
d = D.Mag();
if (d > m_max_distance)
d = 0;
}
// loop over all points
for(unsigned int i = 0 ; i < sizeX ; i++){
P.SetXYZ(X[i],Y[i],Z[i]);
if(MinimumDistancePointLine(Vp1,Vp2,P) < m_match_distance){
cluster.AddIndex(i);
m_assigned.insert(i);
}
}
// insert the newly formed cluster
if(cluster.size()>m_min_count){
m_cluster.insert(cluster);
}
// loop over all points
for (unsigned int i = 0; i < sizeX; i++) {
P.SetXYZ(X[i], Y[i], Z[i]);
if (MinimumDistancePointLine(Vp1, Vp2, P) < m_match_distance) {
cluster.AddIndex(i);
m_assigned.insert(i);
}
}
// insert the newly formed cluster
if (cluster.size() > m_min_count) {
m_cluster.insert(cluster);
}
}//while
// loop over the cluster starting with the biggest
unsigned int current_cluster=0;
unsigned int index;
for(auto it = m_cluster.begin() ; it!=m_cluster.end() ; ++it){
current_cluster++;
// cout << current_cluster << endl;
unsigned int sizeC = (*it).size();
unsigned int cluster_size = 0;
for(unsigned int i = 0 ; i < sizeC ; i++){
// Assigned cluster id to identified points
unsigned int index = (*it).GetIndex(i);
if(!cluster_id[index]){
cluster_id[index]=current_cluster;
cluster_size++;
}
}
if(cluster_size<m_min_count){
for(unsigned int i = 0 ; i < sizeC ; i++){
unsigned int index = (*it).GetIndex(i);
// remove the assigned point
if(cluster_id[index]==current_cluster){
cluster_id[index]=0;
}
}
}
else{
clusters.push_back(*it);
} // while
// loop over the cluster starting with the biggest
unsigned int current_cluster = 0;
unsigned int index;
for (auto it = m_cluster.begin(); it != m_cluster.end(); ++it) {
current_cluster++;
// cout << current_cluster << endl;
unsigned int sizeC = (*it).size();
unsigned int cluster_size = 0;
for (unsigned int i = 0; i < sizeC; i++) {
// Assigned cluster id to identified points
unsigned int index = (*it).GetIndex(i);
if (!cluster_id[index]) {
cluster_id[index] = current_cluster;
cluster_size++;
}
}
if (cluster_size < m_min_count) {
for (unsigned int i = 0; i < sizeC; i++) {
unsigned int index = (*it).GetIndex(i);
// remove the assigned point
if (cluster_id[index] == current_cluster) {
cluster_id[index] = 0;
}
}
return clusters;
}
else {
clusters.push_back(*it);
}
}
return clusters;
}
......@@ -183,7 +183,7 @@ void SuperX3::ReadConfiguration(NPL::InputParser parser) {
if (blocks[i]->GetInt("VIS"))
m_non_sensitive_part_visiualisation = true;
AddDetector(Theta, Phi, R, beta[0], beta[1], beta[2]);
AddDetector(R, Theta, Phi, beta[0], beta[1], beta[2]);
}
else {
......@@ -237,18 +237,18 @@ void SuperX3::ConstructDetector(G4LogicalVolume* world) {
// w perpendicular to (u,v) plan and pointing ThirdStage
// Phi is angle between X axis and projection in (X,Y) plan
// Theta is angle between position vector and z axis
G4double wX = m_R[i] * sin(Theta / rad) * cos(Phi / rad);
G4double wY = m_R[i] * sin(Theta / rad) * sin(Phi / rad);
G4double wZ = m_R[i] * cos(Theta / rad);
G4double wX = m_R[i] * sin(Theta) * cos(Phi);
G4double wY = m_R[i] * sin(Theta) * sin(Phi);
G4double wZ = m_R[i] * cos(Theta);
SuperX3w = G4ThreeVector(wX, wY, wZ);
// vector corresponding to the center of the module
SuperX3Center = SuperX3w;
// vector parallel to one axis of silicon plane
G4double ii = cos(Theta / rad) * cos(Phi / rad);
G4double jj = cos(Theta / rad) * sin(Phi / rad);
G4double kk = -sin(Theta / rad);
G4double ii = cos(Theta) * cos(Phi);
G4double jj = cos(Theta) * sin(Phi);
G4double kk = -sin(Theta);
G4ThreeVector Y = G4ThreeVector(ii, jj, kk);
SuperX3w = SuperX3w.unit();
......@@ -265,7 +265,7 @@ void SuperX3::ConstructDetector(G4LogicalVolume* world) {
SuperX3rot->rotate(m_beta_v[i], SuperX3v);
SuperX3rot->rotate(m_beta_w[i], SuperX3w);
// translation to place Telescope
SuperX3pos = SuperX3w * Length * 0.5 + SuperX3Center;
SuperX3pos = SuperX3w * (SiliconThickness + AluStripThickness) * 0.5 + SuperX3Center;
}
VolumeMaker(i + 1, SuperX3pos, SuperX3rot, world);
......
%%%%%%%
Target
THICKNESS= 0.1 micrometer
ANGLE= 0 deg
RADIUS= 15 mm
MATERIAL= CD2
X= 0 mm
Y= 0 mm
Z= -20 cm
NbSlices= 10
%%%%%%% Alias Phi
% will create a new block for each value of Phi
Alias Phi
Action= Copy
Value= 0 15 30 45 60 75 90 105 120 135 150 165 180 195 210 225 240 255 270 285 300 315 330 345
%%%%%%% Detector 1 %%%%%%%
SuperX3
THETA= 90 deg
PHI= @Phi deg
R= 198.4 mm
BETA= 0 0 0 deg
VIS= all
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment