Commit fd9d7551 authored by Morfouace's avatar Morfouace
Browse files

* Adding Cluster method to find and fit tracks

parent 6e0ecf9e
......@@ -40,22 +40,25 @@ Analysis::~Analysis(){
void Analysis::Init(){
Actar= (TActarPhysics*) m_DetectorManager->GetDetector("Actar");
ReactionConditions = new TReactionConditions();
Actar->ReadAnalysisConfig();
if(Actar->GetRansacStatus()){
Actar->SetRansacParameter("./configs/RansacConfig.dat");
}
else if(Actar->GetClusterStatus()){
Actar->SetClusterParameter("./configs/ClusterConfig.dat");
}
DriftVelocity = Actar->GetDriftVelocity();
PadSizeX = Actar->GetPadSizeX();
PadSizeY = Actar->GetPadSizeY();
NumberOfPadsX = Actar->GetNumberOfPadsX();
NumberOfPadsY = Actar->GetNumberOfPadsY();
EnergyLoss_1H = NPL::EnergyLoss("./EnergyLossTable/proton_iC4H10_6.24151e+07_295.G4table","G4Table",100);
EnergyLoss_18O = NPL::EnergyLoss("./EnergyLossTable/O18_iC4H10_6.24151e+07_295.G4table","G4Table",100);
TheReaction = new NPL::Reaction("18O(p,p)18O@59");
InitInputBranch();
InitOutputBranch();
}
......@@ -63,20 +66,20 @@ void Analysis::Init(){
////////////////////////////////////////////////////////////////////////////////
void Analysis::TreatEvent(){
ReInitValue();
LightName="";
if(ReactionConditions->GetParticleMultiplicity()>0){
InitXVertex = ReactionConditions->GetVertexPositionZ();
InitTheta3 = ReactionConditions->GetTheta(0);
InitE3 = ReactionConditions->GetKineticEnergy(0);
LightName = ReactionConditions->GetParticleName(0);
InitXVertex = ReactionConditions->GetVertexPositionZ();
InitTheta3 = ReactionConditions->GetTheta(0);
InitE3 = ReactionConditions->GetKineticEnergy(0);
LightName = ReactionConditions->GetParticleName(0);
}
int TrackMult = Actar->GetTrackMult();
TVector3 vX = TVector3(1,0,0);
TVector3 aTrack, vB;
if(TrackMult>1){
if(TrackMult==2){
vTrack = Actar->GetTracks();
double scalarproduct=0;
int BeamTrack=0;
......@@ -92,35 +95,35 @@ void Analysis::TreatEvent(){
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();
double vertex_x = vTrack[i].GetVertexPostion(vBeam,vBeamPos).X()*PadSizeX;
double vertex_y = vTrack[i].GetVertexPostion(vBeam,vBeamPos).Y()*PadSizeY;
double vertex_z = 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 = vB.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;
......@@ -129,37 +132,36 @@ void Analysis::TreatEvent(){
//double z2 = -(vTrack[i].GetZh()-256)*DriftVelocity;
double z1 = vTrack[i].GetZm()*DriftVelocity;
double z2 = vTrack[i].GetZh()*DriftVelocity;
if(vertex_x>0 && vertex_x<256){
double LengthInGas = fSiDistanceX - vertex_x;
for(unsigned int k=0; k<Actar->Si_E.size(); k++){
XVertex.push_back(vertex_x);
YVertex.push_back(vertex_y);
ZVertex.push_back(vertex_z);
GetMayaSiHitPosition(x1,x2,y1,y2,z1,z2);
ESi.push_back(Actar->Si_E[k]);
SiNumber.push_back(Actar->Si_Number[k]);
DE.push_back(vTrack[i].GetPartialCharge(108,128)/(20./cos(angle*TMath::Pi()/180)));
double E3;
if(LightName=="proton")E3 = EnergyLoss_1H.EvaluateInitialEnergy(Actar->Si_E[k]*MeV,LengthInGas*mm,angle*TMath::Pi()/180);
if(LightName=="deuteron")E3 = EnergyLoss_2H.EvaluateInitialEnergy(Actar->Si_E[k]*MeV,LengthInGas*mm,angle*TMath::Pi()/180);
if(LightName=="triton")E3 = EnergyLoss_3H.EvaluateInitialEnergy(Actar->Si_E[k]*MeV,LengthInGas*mm,angle*TMath::Pi()/180);
if(LightName=="He3")E3 = EnergyLoss_3He.EvaluateInitialEnergy(Actar->Si_E[k]*MeV,LengthInGas*mm,angle*TMath::Pi()/180);
if(LightName=="alpha")E3 = EnergyLoss_4He.EvaluateInitialEnergy(Actar->Si_E[k]*MeV,LengthInGas*mm,angle*TMath::Pi()/180);
//double BeamEnergy = EnergyLoss_18O.Slow(59.4*MeV,(vertex_x[i]+60)*mm, BeamAngle*TMath::Pi()/180);
double Energy_beam = EnergyLoss_18O.Slow(59.4*MeV,(vertex_x+60)*mm, 0);
BeamEnergy.push_back(Energy_beam);
TheReaction->SetBeamEnergy(Energy_beam);
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());
XVertex.push_back(vertex_x);
YVertex.push_back(vertex_y);
ZVertex.push_back(vertex_z);
GetMayaSiHitPosition(x1,x2,y1,y2,z1,z2);
ESi.push_back(Actar->Si_E[k]);
SiNumber.push_back(Actar->Si_Number[k]);
DE.push_back(vTrack[i].GetPartialCharge(108,128)/(20./cos(angle*TMath::Pi()/180)));
double E3;
if(LightName=="proton")E3 = EnergyLoss_1H.EvaluateInitialEnergy(Actar->Si_E[k]*MeV,LengthInGas*mm,angle*TMath::Pi()/180);
if(LightName=="deuteron")E3 = EnergyLoss_2H.EvaluateInitialEnergy(Actar->Si_E[k]*MeV,LengthInGas*mm,angle*TMath::Pi()/180);
if(LightName=="triton")E3 = EnergyLoss_3H.EvaluateInitialEnergy(Actar->Si_E[k]*MeV,LengthInGas*mm,angle*TMath::Pi()/180);
if(LightName=="He3")E3 = EnergyLoss_3He.EvaluateInitialEnergy(Actar->Si_E[k]*MeV,LengthInGas*mm,angle*TMath::Pi()/180);
if(LightName=="alpha")E3 = EnergyLoss_4He.EvaluateInitialEnergy(Actar->Si_E[k]*MeV,LengthInGas*mm,angle*TMath::Pi()/180);
//double BeamEnergy = EnergyLoss_18O.Slow(59.4*MeV,(vertex_x[i]+60)*mm, BeamAngle*TMath::Pi()/180);
double Energy_beam = EnergyLoss_18O.Slow(59.4*MeV,(vertex_x+60)*mm, 0);
BeamEnergy.push_back(Energy_beam);
TheReaction->SetBeamEnergy(Energy_beam);
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());
}
}
}
......@@ -171,46 +173,46 @@ void Analysis::TreatEvent(){
void Analysis::GetMayaSiHitPosition(double xm, double xh, double ym, double yh, double zm, double zh)
{
double X1, X2, Y1, Y2, Z1, Z2;
if(xm>xh){
X1 = xh;
Y1 = yh;
Z1 = zh;
X2 = xm;
Y2 = ym;
Z2 = zm;
X1 = xh;
Y1 = yh;
Z1 = zh;
X2 = xm;
Y2 = ym;
Z2 = zm;
}
else if(xh>xm){
X1 = xm;
Y1 = ym;
Z1 = zm;
X2 = xh;
Y2 = yh;
Z2 = zh;
X1 = xm;
Y1 = ym;
Z1 = zm;
X2 = xh;
Y2 = yh;
Z2 = zh;
}
double l, L, t;
double zf, yf;
if(fSiDistanceX>X2){
L = fSiDistanceX-X2;
l = X2 - X1;
t = (l+L)/l;
zf = Z1 + (Z2-Z1)*t;
yf = Y1 + (Y2-Y1)*t;
L = fSiDistanceX-X2;
l = X2 - X1;
t = (l+L)/l;
zf = Z1 + (Z2-Z1)*t;
yf = Y1 + (Y2-Y1)*t;
}
else if(fSiDistanceX<X2){
L = X2 - fSiDistanceX;
l = fSiDistanceX - X1;
t = (l+L)/l;
zf = Z1 + (Z2-Z1)/t;
yf = Y1 + (Y2-Y1)/t;
L = X2 - fSiDistanceX;
l = fSiDistanceX - X1;
t = (l+L)/l;
zf = Z1 + (Z2-Z1)/t;
yf = Y1 + (Y2-Y1)/t;
}
SiPosY.push_back(yf);
SiPosZ.push_back(zf);
}
......@@ -221,9 +223,9 @@ void Analysis::End(){
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitInputBranch() {
RootInput::getInstance()->GetChain()->SetBranchStatus("ReactionConditions",true);
RootInput::getInstance()->GetChain()->SetBranchStatus("fRC_*",true);
RootInput::getInstance()->GetChain()->SetBranchAddress("ReactionConditions",&ReactionConditions);
RootInput::getInstance()->GetChain()->SetBranchStatus("ReactionConditions",true);
RootInput::getInstance()->GetChain()->SetBranchStatus("fRC_*",true);
RootInput::getInstance()->GetChain()->SetBranchAddress("ReactionConditions",&ReactionConditions);
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitOutputBranch() {
......@@ -245,7 +247,7 @@ void Analysis::InitOutputBranch() {
RootOutput::getInstance()->GetTree()->Branch("InitXVertex",&InitXVertex,"InitXVertex/D");
RootOutput::getInstance()->GetTree()->Branch("InitE3",&InitE3,"InitE3/D");
RootOutput::getInstance()->GetTree()->Branch("InitTheta3",&InitTheta3,"InitTheta3/D");
}
////////////////////////////////////////////////////////////////////////////////
......@@ -264,7 +266,7 @@ void Analysis::ReInitValue(){
SiPosY.clear();
SiPosZ.clear();
BeamEnergy.clear();
BeamAngle=-1000;
InitE3=-1000;
InitTheta3=-1000;
......@@ -288,6 +290,6 @@ extern "C"{
NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct);
}
};
proxy p;
}
ConfigCluster
NumberOfTracksMax= 20
ClusterDistance= 5
ClusterThreshold= 15
ConfigActar
RecoRansac= 1
RecoVisu= 1
RecoCluster= 0
RecoVisu= 0
HIT_THRESHOLD= 2
Q_THRESHOLD= 0
T_THRESHOLD= 0
......
......@@ -148,7 +148,9 @@ void Sync()
vhxz.push_back(hXZ);
}
/////////////////
// Track in 3D //
/////////////////
c1->cd(1);
hXYZ->Draw();
for(int i=0; i<vh3.size(); i++){
......@@ -197,6 +199,10 @@ void Sync()
vfxz.push_back(f1xz);
}
///////////////////////////
// Track in the XY plan //
//////////////////////////
c1->cd(2);
h2Dxy->Draw();
for(int i=0; i<vhxy.size(); i++){
......@@ -209,6 +215,9 @@ void Sync()
///////////////////////////
// Track in the XZ plan //
//////////////////////////
c1->cd(3);
h2Dxz->Draw();
for(int i=0; i<vhxz.size(); i++){
......
......@@ -54,10 +54,15 @@ m_EventPhysics(this),
m_Spectra(0),
fRecoRansac(1),
fRecoCluster(0),
fRecoVisu(0),
fHitThreshold(20),
fQ_Threshold(0),
fT_Threshold(0),
fXBeamMin(0),
fXBeamMax(128),
fYBeamMin(60),
fYBeamMax(67),
fNumberOfPadsX(128),
fNumberOfPadsY(128),
fPadSizeX(2),
......@@ -111,6 +116,18 @@ void TActarPhysics::BuildPhysicalEvent() {
m_Track = m_Ransac->SimpleRansac();
}
else if(fRecoCluster && PadX.size()>fHitThreshold){
m_Cluster->Init(PadX, PadY, PadZ, PadCharge);
m_Cluster->FindTracks();
m_Cluster->Init(BeamPadX, BeamPadY, BeamPadZ, BeamPadCharge);
m_Cluster->FindTracks();
m_Track = m_Cluster->GetTracks();
m_Cluster->Clear();
}
TrackMult = m_Track.size();
}
......@@ -143,11 +160,27 @@ void TActarPhysics::PreTreat() {
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){
if(fRecoCluster){
if(IsBeamZone(TABLE[4][where],TABLE[5][where])){
BeamPadCharge.push_back(m_EventReduced->CoboAsad[it].peakheight[hh]);
BeamPadX.push_back(TABLE[4][where]);
BeamPadY.push_back(TABLE[5][where]);
BeamPadZ.push_back(m_EventReduced->CoboAsad[it].peaktime[hh]);
}
else{
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{
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]);
}
//}
}
}
......@@ -169,6 +202,18 @@ void TActarPhysics::PreTreat() {
}
}
}
///////////////////////////////////////////////////////////////////////////
bool TActarPhysics::IsBeamZone(int X, int Y)
{
bool isBeam=false;
if(X>fXBeamMin && X<fXBeamMax)
if(Y>fYBeamMin && fYBeamMax)
isBeam=true;
return isBeam;
}
///////////////////////////////////////////////////////////////////////////
bool TActarPhysics::GoodHit(int iX, int iY)
{
......@@ -303,6 +348,12 @@ void TActarPhysics::ReadAnalysisConfig() {
cout << "/// Reco using Ransac= " << " " << fRecoRansac << " ///" << endl;
}
else if (whatToDo.compare(0,12,"RecoCluster=") == 0) {
AnalysisConfigFile >> DataBuffer;
fRecoCluster = atoi(DataBuffer.c_str());
cout << "/// Reco using Cluster= " << " " << fRecoCluster << " ///" << endl;
}
else if (whatToDo.compare(0,9,"RecoVisu=") == 0) {
AnalysisConfigFile >> DataBuffer;
fRecoVisu = atoi(DataBuffer.c_str());
......@@ -354,6 +405,30 @@ void TActarPhysics::ReadAnalysisConfig() {
cout << "/// Pad Size Y= " << fPadSizeY << " ///" << endl;
}
else if(whatToDo.compare(0,9,"XBeamMin=")==0){
AnalysisConfigFile >> DataBuffer;
fXBeamMin = atof(DataBuffer.c_str());
cout << "/// X Beam Min= " << fXBeamMin << " ///" << endl;
}
else if(whatToDo.compare(0,9,"XBeamMax=")==0){
AnalysisConfigFile >> DataBuffer;
fXBeamMax = atof(DataBuffer.c_str());
cout << "/// X Beam Max= " << fXBeamMax << " ///" << endl;
}
else if(whatToDo.compare(0,9,"YBeamMin=")==0){
AnalysisConfigFile >> DataBuffer;
fYBeamMin = atof(DataBuffer.c_str());
cout << "/// Y Beam Min= " << fYBeamMin << " ///" << endl;
}
else if(whatToDo.compare(0,9,"YBeamMax=")==0){
AnalysisConfigFile >> DataBuffer;
fYBeamMax = atof(DataBuffer.c_str());
cout << "/// Y Beam Max= " << fYBeamMax << " ///" << endl;
}
else if(whatToDo.compare(0,9,"Pressure=")==0){
AnalysisConfigFile >> DataBuffer;
fPressure = atof(DataBuffer.c_str());
......@@ -384,6 +459,9 @@ void TActarPhysics::ReadAnalysisConfig() {
if(fRecoRansac){
m_Ransac = new NPL::Ransac(fNumberOfPadsX,fNumberOfPadsY,fRecoVisu);
}
else if(fRecoCluster){
m_Cluster = new NPL::Cluster(fNumberOfPadsX,fNumberOfPadsY,fRecoVisu);
}
}
///////////////////////////////////////////////////////////////////////////
......@@ -393,8 +471,20 @@ void TActarPhysics::SetRansacParameter(string filename){
}
}
///////////////////////////////////////////////////////////////////////////
void TActarPhysics::SetClusterParameter(string filename){
if(fRecoCluster){
m_Cluster->ReadParameterValue(filename);
}
}
///////////////////////////////////////////////////////////////////////////
void TActarPhysics::Clear() {
BeamPadX.clear();
BeamPadY.clear();
BeamPadZ.clear();
BeamPadCharge.clear();
PadX.clear();
PadY.clear();
PadZ.clear();
......
......@@ -43,6 +43,7 @@ using namespace std;
#include "NPInputParser.h"
#include "NPTrack.h"
#include "NPRansac.h"
#include "NPCluster.h"
#define NumberOfCobo 16
#define NumberOfASAD 4
......@@ -74,6 +75,10 @@ public:
vector<int> PadY;
vector<double> PadZ;
vector<double> PadCharge;
vector<int> BeamPadX;
vector<int> BeamPadY;
vector<double> BeamPadZ;
vector<double> BeamPadCharge;
vector<double> Si_E;
vector<int> Si_Number;
int TrackMult;
......@@ -159,6 +164,8 @@ public:
bool GoodHit(int iX, int iY);
bool IsBeamZone(int X, int Y);
// give and external TActarData object to TActarPhysics.
// needed for online analysis for example
void SetRawDataPointer(TActarData* rawDataPointer) {m_EventData = rawDataPointer;}
......@@ -170,6 +177,7 @@ private:
TActarData* m_PreTreatedData; //!
TActarPhysics* m_EventPhysics; //!
NPL::Ransac* m_Ransac; //!
NPL::Cluster* m_Cluster; //!
vector<NPL::Track> m_Track; //!
// getters for raw and pre-treated data object
......@@ -193,12 +201,17 @@ private:
int fT_Threshold; //!
int fNumberOfPadsX; //!
int fNumberOfPadsY; //!
int fXBeamMax; //!
int fXBeamMin; //!
int fYBeamMax; //!
int fYBeamMin; //!
double fPadSizeX; //!
double fPadSizeY; //!
double fDriftVelocity; //!
double fPressure; //!
string fGas; //!
bool fRecoRansac; //!
bool fRecoCluster; //!
bool fRecoVisu; //!
map<int, int> Si_map; //!
......@@ -217,8 +230,11 @@ private:
//spme getters and setters
public:
void SetRansacParameter(string filename);
void SetClusterParameter(string filename);
NPL::Ransac* GetRansacObject() {return m_Ransac;}
bool GetRansacStatus() {return fRecoRansac;}
NPL::Cluster* GetClusterObject() {return m_Cluster;}
bool GetClusterStatus() {return fRecoCluster;}
// spectra getter
......
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)
add_custom_command(OUTPUT NPClusterDict.cxx COMMAND ../scripts/build_dict.sh NPCluster.h NPClusterDict.cxx NPCluster.rootmap libNPTrackReconstruction.so NPTrackReconstructionLinkDef.h DEPENDS NPCluster.h)
add_library(NPTrackReconstruction SHARED NPRansac.cxx NPCluster.cxx NPTrack.cxx NPRansacDict.cxx NPClusterDict.cxx)
target_link_libraries(NPTrackReconstruction ${ROOT_LIBRARIES} NPCore)
install(FILES NPRansac.h NPTrack.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
install(FILES NPRansac.h NPCluster.h NPTrack.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
This diff is collapsed.
#ifndef __CLUSTER__
#define __CLUSTER__
/*****************************************************************************
* 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 : September 2018 *
* Last update : *
*---------------------------------------------------------------------------*
* 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>