Commit 5abb3b07 authored by Adrien Matta's avatar Adrien Matta
Browse files

Merge branch 'NPTool.2.dev' of gitlab.in2p3.fr:np/nptool into NPTool.2.dev

parents c0bafd0d e397d977
Pipeline #85822 passed with stages
in 55 minutes and 31 seconds
......@@ -78,7 +78,6 @@ void Analysis::TreatEvent(){
TVector3 vX = TVector3(1,0,0);
TVector3 aTrack, vB;
if(TrackMult>1){
vTrack = Actar->GetTracks();
double scalarproduct=0;
......@@ -106,7 +105,8 @@ void Analysis::TreatEvent(){
double ZBeamPoint = vTrack[BeamTrack].GetZh();
TVector3 vBeamPos = TVector3(XBeamPoint,YBeamPoint,ZBeamPoint);
vB = TVector3(XBeam*PadSizeX, YBeam*PadSizeY,ZBeam*DriftVelocity);
//vB = TVector3(XBeam*PadSizeX, YBeam*PadSizeY,ZBeam*DriftVelocity);
vB = TVector3(XBeam*PadSizeX, YBeam*PadSizeY,ZBeam);
BeamAngle = (vX.Angle(vB))*180/TMath::Pi();
for(unsigned int i=0; i<TrackMult; i++){
......@@ -117,9 +117,11 @@ void Analysis::TreatEvent(){
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;
//double vertex_z = vTrack[i].GetVertexPostion(vBeam,vBeamPos).Z()*DriftVelocity;
double vertex_z = vTrack[i].GetVertexPostion(vBeam,vBeamPos).Z();
aTrack = TVector3(Xdir*PadSizeX, Ydir*PadSizeY, Zdir*DriftVelocity);
//aTrack = TVector3(Xdir*PadSizeX, Ydir*PadSizeY, Zdir*DriftVelocity);
aTrack = TVector3(Xdir*PadSizeX, Ydir*PadSizeY, Zdir);
double angle = vX.Angle(aTrack)*180/TMath::Pi();
//double angle = vB.Angle(aTrack)*180/TMath::Pi();
if(angle>90) angle = 180-angle;
......@@ -130,8 +132,10 @@ void Analysis::TreatEvent(){
double y2 = vTrack[i].GetYh()*PadSizeY-0.5*NumberOfPadsY*PadSizeY;
//double z1 = -(vTrack[i].GetZm()-256)*DriftVelocity;
//double z2 = -(vTrack[i].GetZh()-256)*DriftVelocity;
double z1 = vTrack[i].GetZm()*DriftVelocity;
double z2 = vTrack[i].GetZh()*DriftVelocity;
//double z1 = vTrack[i].GetZm()*DriftVelocity;
double z1 = vTrack[i].GetZm();
//double z2 = vTrack[i].GetZh()*DriftVelocity;
double z2 = vTrack[i].GetZh();
if(vertex_x>0 && vertex_x<256){
......
ConfigActar
RecoRansac= 0
RecoCluster= 1
RecoVisu= 0
RecoRansac= 1
RecoCluster= 0
RecoVisu= 1
HIT_THRESHOLD= 2
Q_THRESHOLD= 0
T_THRESHOLD= 0
......@@ -17,4 +17,4 @@ Gas= iC4H10
%Gas= D2
%Pressure in mbar
Pressure= 100
DriftVelocity= 2.66
DriftVelocity= 0.8e-2
......@@ -17,7 +17,7 @@ void LoadCut()
void LoadChain()
{
chain = new TChain("PhysicsTree");
chain->Add("../../Outputs/Analysis/Example4.root");
chain->Add("../../../Outputs/Analysis/Example4.root");
}
////////////////////////////////////////////////
......
......@@ -18,7 +18,7 @@ void LoadCut()
void LoadChain()
{
chain = new TChain("PhysicsTree");
chain->Add("../../Outputs/Analysis/Example4.root");
chain->Add("../../../Outputs/Analysis/Example4.root");
}
////////////////////////////////////////////////
......
......@@ -45,55 +45,55 @@ using namespace std;
ClassImp(TActarPhysics)
///////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////
TActarPhysics::TActarPhysics()
: m_EventData(new TActarData),
m_EventReduced(new MEventReduced),
m_PreTreatedData(new TActarData),
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),
fPadSizeY(2),
fDriftVelocity(40),
fPressure(100),
fGas("iC4H10"),
m_NumberOfPadSilicon(20),
m_NumberOfDetectors(0) {
}
: m_EventData(new TActarData),
m_EventReduced(new MEventReduced),
m_PreTreatedData(new TActarData),
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),
fPadSizeY(2),
fDriftVelocity(40),
fPressure(100),
fGas("iC4H10"),
m_NumberOfPadSilicon(20),
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();
}
......@@ -101,570 +101,570 @@ void TActarPhysics::BuildSimplePhysicalEvent() {
///////////////////////////////////////////////////////////////////////////
void TActarPhysics::BuildPhysicalEvent() {
PreTreat();
PreTreat();
//unsigned int mysize = m_PreTreatedData->GetPadMult();
//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));
/*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();
if(fRecoRansac && PadX.size()>fHitThreshold){
//m_Ransac->Init(PadX, PadY, PadZ, PadCharge);
m_Ransac->Init(PadX, PadY, PadZ, PadCharge);
m_Track = m_Ransac->SimpleRansac();
}
else if(fRecoCluster){
if(PadX.size()>fHitThreshold){
m_Cluster->Init(PadX, PadY, PadZ, PadCharge);
m_Cluster->FindTracks();
}
else if(fRecoCluster){
if(PadX.size()>fHitThreshold){
m_Cluster->Init(PadX, PadY, PadZ, PadCharge);
m_Cluster->FindTracks();
}
if(BeamPadX.size()>fHitThreshold){
m_Cluster->Init(BeamPadX, BeamPadY, BeamPadZ, BeamPadCharge);
m_Cluster->FindTracks();
}
if(BeamPadX.size()>fHitThreshold){
m_Cluster->Init(BeamPadX, BeamPadY, BeamPadZ, BeamPadCharge);
m_Cluster->FindTracks();
}
m_Track = m_Cluster->GetTracks();
m_Track = m_Cluster->GetTracks();
m_Cluster->Clear();
}
m_Cluster->Clear();
}
TrackMult = m_Track.size();
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();
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){
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 if(fRecoRansac){
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]);
}
//}
}
}
// 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();
//cout << vector_size << endl;
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){
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(int(m_EventReduced->CoboAsad[it].peaktime[hh]*fDriftVelocity));
}
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(int(m_EventReduced->CoboAsad[it].peaktime[hh]*fDriftVelocity));
}
}
}
else if(co==31){
unsigned int vector_size = m_EventReduced->CoboAsad[it].peakheight.size();
for(unsigned int hit=0;hit<vector_size;hit++){
if(fInputTreeName=="ACTAR_TTree"){
int vxi_parameter = m_EventReduced->CoboAsad[it].peaktime[hit];
if(Si_map[vxi_parameter]<21 && Si_map[vxi_parameter]>0){
double RawEnergy = m_EventReduced->CoboAsad[it].peakheight[hit];
static string name;
name = "ActarSi/D" ;
name+= NPL::itoa( Si_map[vxi_parameter] - 1) ;
name+= "_E" ;
double CalEnergy = CalibrationManager::getInstance()->ApplyCalibration( name, RawEnergy );
Si_Number.push_back(Si_map[vxi_parameter]);
Si_E.push_back(CalEnergy);
}
}
else{
Si_Number.push_back(m_EventReduced->CoboAsad[it].peaktime[hit]);
Si_E.push_back(m_EventReduced->CoboAsad[it].peakheight[hit]);
}
else if(fRecoRansac){
PadCharge.push_back(m_EventReduced->CoboAsad[it].peakheight[hh]);
PadX.push_back(TABLE[4][where]);
PadY.push_back(TABLE[5][where]);
PadZ.push_back(int(m_EventReduced->CoboAsad[it].peaktime[hh]*fDriftVelocity));
}
//}
}
}
}
}
else if(co==31){
unsigned int vector_size = m_EventReduced->CoboAsad[it].peakheight.size();
for(unsigned int hit=0;hit<vector_size;hit++){
if(fInputTreeName=="ACTAR_TTree"){
int vxi_parameter = m_EventReduced->CoboAsad[it].peaktime[hit];
if(Si_map[vxi_parameter]<21 && Si_map[vxi_parameter]>0){
double RawEnergy = m_EventReduced->CoboAsad[it].peakheight[hit];
static string name;
name = "ActarSi/D" ;
name+= NPL::itoa( Si_map[vxi_parameter] - 1) ;
name+= "_E" ;
double CalEnergy = CalibrationManager::getInstance()->ApplyCalibration( name, RawEnergy );
Si_Number.push_back(Si_map[vxi_parameter]);
Si_E.push_back(CalEnergy);
}
}
else{
Si_Number.push_back(m_EventReduced->CoboAsad[it].peaktime[hit]);
Si_E.push_back(m_EventReduced->CoboAsad[it].peakheight[hit]);
}
}
}
}
}
///////////////////////////////////////////////////////////////////////////
bool TActarPhysics::IsBeamZone(int X, int Y)
{
bool isBeam=false;
if( (X>=fXBeamMin && X<=fXBeamMax) && (Y>=fYBeamMin && Y<=fYBeamMax) ){
isBeam=true;
}
bool isBeam=false;
if( (X>=fXBeamMin && X<=fXBeamMax) && (Y>=fYBeamMin && Y<=fYBeamMax) ){
isBeam=true;
}
return isBeam;
return isBeam;
}
///////////////////////////////////////////////////////////////////////////
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;
}
}
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;
return bHit;
}
///////////////////////////////////////////////////////////////////////////
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]);*/
}
}
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() {
// 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;
}
// 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;
}
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];
}
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;
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());