From 13436990731a370957f9cae1f34bb0fde5421d21 Mon Sep 17 00:00:00 2001
From: Hugo Jacob <hugojacob57@gmail.com>
Date: Mon, 13 Nov 2023 14:44:28 +0100
Subject: [PATCH] Modified CATS mask calibration code
---
NPLib/Detectors/CATS/TCATSPhysics.cxx | 7 +-
Projects/E805/Analysis.cxx | 91 ++++++++++--
Projects/E805/Analysis.h | 19 ++-
Projects/E805/Analysis_E805.cxx | 7 +-
Projects/E805/Calibration1.txt | 41 ++++++
Projects/E805/Calibration2.txt | 41 ++++++
.../DetectorConfiguration/CATS_MCR.detector | 16 +--
.../DetectorConfiguration/MUST2_E805.detector | 18 +--
Projects/E805/ProjetInfo/MCRMethod.cpp | 132 +++++++++++++++++-
Projects/E805/ProjetInfo/MCRMethod.h | 17 ++-
10 files changed, 340 insertions(+), 49 deletions(-)
create mode 100644 Projects/E805/Calibration1.txt
create mode 100644 Projects/E805/Calibration2.txt
diff --git a/NPLib/Detectors/CATS/TCATSPhysics.cxx b/NPLib/Detectors/CATS/TCATSPhysics.cxx
index 78ad6c461..adc6fc571 100644
--- a/NPLib/Detectors/CATS/TCATSPhysics.cxx
+++ b/NPLib/Detectors/CATS/TCATSPhysics.cxx
@@ -223,9 +223,10 @@ void TCATSPhysics::BuildPhysicalEvent(){
// });
// At least two CATS need to gave back position in order to reconstruct on Target
if(Positions.size()>1){
- double t = (m_Zproj-Positions[2].first)/(Positions[2].first-Positions[1].first);
- PositionOnTargetX= Positions[2].second.first + (Positions[2].second.first-Positions[1].second.first)*t;
- PositionOnTargetY= Positions[2].second.second + (Positions[2].second.second-Positions[1].second.second)*t;
+ double t = (m_Zproj-Positions[2].first)/(m_Zproj-Positions[1].first);
+ PositionOnTargetX= (Positions[2].second.first - Positions[1].second.first*t)/(1 - t);
+ PositionOnTargetY= (Positions[2].second.second - Positions[1].second.second*t)/(1 - t);
+ // PositionOnTargetY= Positions[2].second.second + (Positions[2].second.second-3-Positions[1].second.second)*t;
if(Mask1_Z != 0 && Mask2_Z != 0)
{
double tmask1 = (Mask1_Z-Positions[2].first)/(Positions[2].first-Positions[1].first);
diff --git a/Projects/E805/Analysis.cxx b/Projects/E805/Analysis.cxx
index 1ca337ef6..23ffc7b24 100755
--- a/Projects/E805/Analysis.cxx
+++ b/Projects/E805/Analysis.cxx
@@ -75,13 +75,14 @@ void Analysis::Init(){
bool Analysis::UnallocateBeforeBuild(){
// std::cout << "test unallocate" << std::endl;
+ return true;
GATCONFMASTER = **GATCONFMASTER_;
return (GATCONFMASTER > 0);
//return true;
}
bool Analysis::UnallocateBeforeTreat(){
- for(int i = 0; i < 10; i++){
+/* for(int i = 0; i < 10; i++){
PlasticRaw[i] = (*PlasticRaw_ )[i];
PlasticRawTS[i] = (*PlasticRaw_TS_)[i];
}
@@ -125,6 +126,7 @@ bool Analysis::UnallocateBeforeTreat(){
TAC_PL_4TS = **TAC_PL_4_TS_;
TAC_PL_5 = **TAC_PL_5_;
TAC_PL_5TS = **TAC_PL_5_TS_;
+*/
return true;
}
@@ -137,8 +139,7 @@ bool Analysis::FillOutputCondition(){
void Analysis::TreatEvent(){
// if(M2->CsI_E.size() > 0)
- // std::cout << "Analysis test " << M2->CsI_E[0] << " " << M2->CsI_E.size() << " " << "\n \n";
-
+ // std::cout << "Analysis test " << M2->CsI_E[0] << " " << M2->CsI_E.size() << " " << "\n \n";*
ReInit();
// std::cout << CATS->PositionX.size() << std::endl;
//////////////////// MUST2 Part ////////////////////
@@ -158,7 +159,7 @@ void Analysis::TreatEvent(){
//BeamImpact = TVector3(0,0,0);
- BeamDirection = TVector3(CATS->PositionX[1] - CATS->PositionX[0],CATS->PositionY[1] - CATS->PositionY[0],CATS->PositionZ[1] - CATS->PositionZ[0]);
+ BeamDirection = TVector3(CATS->PositionX[0] - CATS->PositionX[1],CATS->PositionY[0] - CATS->PositionY[1],CATS->PositionZ[0] - CATS->PositionZ[1]);
// std::cout << "Position XY " << CATS->PositionX[1] - CATS->PositionX[0] << " " << CATS->PositionY[1] - CATS->PositionY[0] << " " << CATS->PositionZ[1] - CATS->PositionZ[0] << std::endl;
// BeamDirection = TVector3(0,0,1);
@@ -200,7 +201,6 @@ void Analysis::TreatEvent(){
Energy[ParticleType[i]] = LightTarget[ParticleType[i]].EvaluateInitialEnergy(Energy[ParticleType[i]] ,TargetThickness*0.5, ThetaNormalTarget);
else
Energy[ParticleType[i]] = -1000;
- // std::cout << "DILO tata " << Energy[ParticleType[i]] << std::endl;
}
@@ -214,6 +214,7 @@ void Analysis::TreatEvent(){
// Part 3 : Excitation Energy Calculation
M2_Ex_p.push_back(reaction->ReconstructRelativistic( Energy["proton"] , M2_ThetaLab[countMust2] ));
M2_Ex_d.push_back(Reaction_pd->ReconstructRelativistic( Energy["deuteron"] , M2_ThetaLab[countMust2] ));
+ std::cout << "oui " << M2_Ex_d[countMust2] << std::endl;
M2_Ex_t.push_back(Reaction_pt->ReconstructRelativistic( Energy["triton"] , M2_ThetaLab[countMust2] ));
M2_Ex_a.push_back(Reaction_p3He->ReconstructRelativistic( Energy["alpha"] , M2_ThetaLab[countMust2] ));
@@ -232,6 +233,47 @@ void Analysis::TreatEvent(){
//
}//end loop MUST2
}
+/*
+ }
+ if(M2->Si_E.size() > 0){
+ if(Inner6MVM > 0){
+
+ int ExoMultMax = 2;
+ if(OutersVM < 2){
+ ExoMultMax == OutersVM;
+ }
+ for(int ExoMult = 0; ExoMult < ExoMultMax;ExoMult++){
+ ExogamDetNb[ExoMult] = (OutersVN[ExoMult] - 32)/16;
+ CristalNb[ExoMult] = (OutersVN[ExoMult] - (2+ ExogamDetNb[ExoMult])*16)/4;
+ SegmentNb[ExoMult] = (OutersVN[ExoMult] - 16*(ExogamDetNb[ExoMult] + 2) - 4*(CristalNb[ExoMult]));
+ }
+
+ if(Inner6MVM == 1){
+ if(OutersVM > 0 && BGOV[0] < 20){
+ Theta_seg = Exogam_Clovers_struc[ExogamDetNb[0]].Theta_Crystal_Seg[CristalNb[0]][SegmentNb[0]];
+ Phi_seg = Exogam_Clovers_struc[ExogamDetNb[0]].Phi_Crystal_Seg[CristalNb[0]][SegmentNb[0]];
+ EnergyDoppler = Doppler_Correction(Theta_seg,Phi_seg,0,0,Beta,Inner6MV[0]);
+ EnergyAddBackDoppler = EnergyDoppler;
+ }
+ }
+ if(Inner6MVM == 2 && OutersVM > 1 && BGOV[0] < 20 && BGOV[1] < 20){
+ if(Inner6MVN[0]/4 == Inner6MVN[1]/4){
+ EnergyAddBack = Inner6MV[0] + Inner6MV[1];
+ }
+ }
+ if(EnergyAddBack > -1000){
+ for(int i = 0;i < ExoMultMax; i++){
+ if(OutersV[i] > OutersV[highest_E]){
+ highest_E = i;
+ }
+ }
+ Theta_seg = Exogam_Clovers_struc[ExogamDetNb[highest_E]].Theta_Crystal_Seg[CristalNb[highest_E]][SegmentNb[highest_E]];
+ Phi_seg = Exogam_Clovers_struc[ExogamDetNb[highest_E]].Phi_Crystal_Seg[CristalNb[highest_E]][SegmentNb[highest_E]];
+ EnergyAddBackDoppler = Doppler_Correction(Theta_seg,Phi_seg,0,0,Beta,EnergyAddBack);
+ }
+ }
+ }
+
//for(unsigned int countMust2 = 0 ; countMust2 < M2->Si_E.size() ; countMust2++){
@@ -254,11 +296,13 @@ void Analysis::TreatEvent(){
// }
// }
//}
+*/
}
void Analysis::InitOutputBranch() {
+ /*
RootOutput::getInstance()->GetTree()->Branch("M2_TelescopeM",&M2_TelescopeM,"M2_TelescopeM/s");
RootOutput::getInstance()->GetTree()->Branch("M2_CsI_E_p",&M2_CsI_E_p);
RootOutput::getInstance()->GetTree()->Branch("M2_CsI_E_d",&M2_CsI_E_d);
@@ -278,7 +322,7 @@ void Analysis::InitOutputBranch() {
RootOutput::getInstance()->GetTree()->Branch("M2_Z",&M2_Z);
RootOutput::getInstance()->GetTree()->Branch("M2_dE",&M2_dE);
RootOutput::getInstance()->GetTree()->Branch("CsI_E_M2",&CsI_E_M2);
- RootOutput::getInstance()->GetTree()->Branch("M2_ECsI_from_deltaE",&M2_ECsI_from_deltaE);
+ // RootOutput::getInstance()->GetTree()->Branch("M2_ECsI_from_deltaE",&M2_ECsI_from_deltaE);
RootOutput::getInstance()->GetTree()->Branch("GATCONF",&GATCONFMASTER);
RootOutput:: getInstance()->GetTree()->Branch("TAC_CATS_PL",&TAC_CATS_PL,"TAC_CATS_PL/s");
@@ -321,15 +365,32 @@ void Analysis::InitOutputBranch() {
RootOutput:: getInstance()->GetTree()->Branch("IC_ZDDRaw",IC_ZDDRaw,"IC_ZDDRaw[6]/s");
RootOutput:: getInstance()->GetTree()->Branch("IC_ZDDRawTS",IC_ZDDRawTS,"IC_ZDDRawTS[6]/l");
+ RootOutput:: getInstance()->GetTree()->Branch("EnergyDoppler",&EnergyDoppler,"EnergyDoppler/F");
+ RootOutput:: getInstance()->GetTree()->Branch("EnergyAddBack",&EnergyAddBack,"EnergyAddBack/F");
+ RootOutput:: getInstance()->GetTree()->Branch("EnergyAddBackDoppler",&EnergyAddBackDoppler,"EnergyAddBackDoppler/F");
+
+ RootOutput:: getInstance()->GetTree()->Branch("Inner6MVM",&Inner6MVM,"Inner6MVM/I");
+ RootOutput:: getInstance()->GetTree()->Branch("Inner6MV",Inner6MV,"Inner6MV[Inner6MVM]/F");
+ RootOutput:: getInstance()->GetTree()->Branch("Inner6MVN",Inner6MVN,"Inner6MVN[Inner6MVM]/s");
+ RootOutput:: getInstance()->GetTree()->Branch("Inner6MVTS",Inner6MVTS,"Inner6MVTS[Inner6MVM]/l");
+ RootOutput:: getInstance()->GetTree()->Branch("BGOVM",&BGOVM,"BGOVM/I");
+ RootOutput:: getInstance()->GetTree()->Branch("BGOV",BGOV,"BGOV[BGOVM]/F");
+ RootOutput:: getInstance()->GetTree()->Branch("BGOVN",BGOVN,"BGOVN[BGOVM]/s");
+ RootOutput:: getInstance()->GetTree()->Branch("DeltaTVM",&DeltaTVM,"DeltaTVM/I");
+ RootOutput:: getInstance()->GetTree()->Branch("DeltaTV",DeltaTV,"DeltaTV[DeltaTVM]/F");
+ RootOutput:: getInstance()->GetTree()->Branch("DeltaTVN",DeltaTVN,"DeltaTVN[DeltaTVM]/s");
+ RootOutput:: getInstance()->GetTree()->Branch("DeltaTVTS",DeltaTVTS,"DeltaTVTS[DeltaTVM]/l");
+*/
}
void Analysis::UnallocateVariables(){
}
void Analysis::InitInputBranch(){
+
TTreeReader* inputTreeReader = RootInput::getInstance()->GetTreeReader();
GATCONFMASTER_ = new TTreeReaderValue<unsigned short>(*inputTreeReader,"GATCONFMASTER");
- //DATATRIG_CATS_ = new TTreeReaderValue<unsigned short>(*inputTreeReader,"DATATRIG_CATS");
+/* //DATATRIG_CATS_ = new TTreeReaderValue<unsigned short>(*inputTreeReader,"DATATRIG_CATS");
PlasticRaw_ = new TTreeReaderArray<UShort_t>(*inputTreeReader,"PlasticRaw");
PlasticRaw_TS_ = new TTreeReaderArray<ULong64_t>(*inputTreeReader,"PlasticRawTS");
@@ -370,6 +431,20 @@ void Analysis::InitInputBranch(){
TAC_PL_4_TS_= new TTreeReaderValue<ULong64_t>(*inputTreeReader,"TAC_PL_4TS");
TAC_PL_5_= new TTreeReaderValue<UShort_t>(*inputTreeReader,"TAC_PL_5");
TAC_PL_5_TS_= new TTreeReaderValue<ULong64_t>(*inputTreeReader,"TAC_PL_5TS");
+
+ Inner6MVM_ = new TTreeReaderValue<int>(*inputTreeReader,"Inner6MRawM");
+ Inner6MV_ = new TTreeReaderArray<float>(*inputTreeReader,"Inner6MRaw");
+ Inner6MVN_ = new TTreeReaderArray<unsigned short>(*inputTreeReader,"Inner6MRawNr");
+ //RootInput:: getInstance()->GetChain()->SetBranchAddress("Inner6MVTS",Inner6MVTS);
+
+ OutersVM_ = new TTreeReaderValue<int>(*inputTreeReader,"OutersRawM");
+ OutersV_ = new TTreeReaderArray<float>(*inputTreeReader,"OutersRaw");
+ OutersVN_ = new TTreeReaderArray<unsigned short>(*inputTreeReader,"OutersRawNr");
+
+ BGOVM_ = new TTreeReaderValue<int>(*inputTreeReader,"BGORawM");
+ BGOV_ = new TTreeReaderArray<float>(*inputTreeReader,"BGORaw");
+ BGOVN_ = new TTreeReaderArray<unsigned short>(*inputTreeReader,"BGORawNr");
+*/
}
////////////////////////////////////////////////////////////////////////////////
@@ -392,7 +467,7 @@ void Analysis::ReInit(){
M2_CsI_E_t.clear();
M2_CsI_E_a.clear();
- M2_ECsI_from_deltaE.clear();
+ // M2_ECsI_from_deltaE.clear();
//ExNoBeam=ExNoProto.clear();
//EDC.clear();
M2_ELab.clear();
diff --git a/Projects/E805/Analysis.h b/Projects/E805/Analysis.h
index 903c00996..8681eca29 100755
--- a/Projects/E805/Analysis.h
+++ b/Projects/E805/Analysis.h
@@ -105,9 +105,7 @@ class Analysis: public NPL::VAnalysis{
std::vector<double> M2_Y;
std::vector<double> M2_Z;
std::vector<double> M2_dE;
- std::vector<double> M2_ECsI_from_deltaE ;
- std::vector<double> Beta_from_deltaE;
- std::vector<double> Beth_from_deltaE;
+
double OriginalBeamEnergy ; // AMEV
double FinalBeamEnergy;
@@ -133,7 +131,14 @@ class Analysis: public NPL::VAnalysis{
unsigned short Inner6MVN[48];
TTreeReaderArray<unsigned short>* Inner6MVN_;
unsigned long long Inner6MVTS[48];
+ TTreeReaderArray<unsigned long long>* Inner6MVTS_;
+ int OutersVM;
+ TTreeReaderValue<int>* OutersVM_;
+ float OutersV[192];
+ TTreeReaderArray<float>* OutersV_;
+ unsigned short OutersVN[192];
+ TTreeReaderArray<unsigned short>* OutersVN_;
int BGOVM;
TTreeReaderValue<int>* BGOVM_;
@@ -169,19 +174,13 @@ class Analysis: public NPL::VAnalysis{
float EnergyAddBackDoppler;
float EnergyAddBack;
int ExogamDetNb[3];
- // int CristalNb[3];
+ int CristalNb[3];
int SegmentNb[3];
std::vector<int> event1;
std::vector<int> event2;
int highest_E;
- int OutersVM;
- TTreeReaderValue<int>* OutersVM_;
- float OutersV[192];
- TTreeReaderArray<float>* OutersV_;
- unsigned short OutersVN[192];
- TTreeReaderArray<unsigned short>* OutersVN_;
int DCRawM;
unsigned short DCRaw[4];
diff --git a/Projects/E805/Analysis_E805.cxx b/Projects/E805/Analysis_E805.cxx
index c03741e94..2e936e70f 100755
--- a/Projects/E805/Analysis_E805.cxx
+++ b/Projects/E805/Analysis_E805.cxx
@@ -21,7 +21,7 @@
#include<iostream>
using namespace std;
-#include"Analysis.h"
+#include"Analysis_E805.h"
#include"NPAnalysisFactory.h"
#include"NPDetectorManager.h"
////////////////////////////////////////////////////////////////////////////////
@@ -218,8 +218,6 @@ bool Analysis::UnallocateBeforeBuild(){
////////////////////////////////////////////////////////////////////////////////
void Analysis::TreatEvent(){
-}
-/*
//std::cout << "Bonjour je suis Nicolas" << std::endl;
ReInit();
@@ -303,7 +301,7 @@ void Analysis::TreatEvent(){
/////////////////////////////// Drift Chambers ///////////////////
- /*if(DCRawM <=2 && DATATRIG_CATSTS[0] > 0){
+ if(DCRawM <=2 && DATATRIG_CATSTS[0] > 0){
for(int DCmult = 0; DCmult < DCRawM; DCmult++){
DC_Nr = DCRawNr[DCmult];
@@ -433,7 +431,6 @@ void Analysis::TreatEvent(){
//std::cout << "Bonjour Nicolas est parti " << std::endl;
}
-*/
diff --git a/Projects/E805/Calibration1.txt b/Projects/E805/Calibration1.txt
new file mode 100644
index 000000000..165ab033d
--- /dev/null
+++ b/Projects/E805/Calibration1.txt
@@ -0,0 +1,41 @@
+CalibrationFilePath
+ ./data/NPRoot/Calibration/CalAlpha_r0388_06+/Cal_Str_X_E_MM1.cal
+ ./data/NPRoot/Calibration/CalAlpha_r0388_06+/Cal_Str_X_E_MM2.cal
+ ./data/NPRoot/Calibration/CalAlpha_r0388_06+/Cal_Str_X_E_MM3.cal
+ ./data/NPRoot/Calibration/CalAlpha_r0388_06+/Cal_Str_X_E_MM4.cal
+
+ ./data/NPRoot/Calibration/CalAlpha_r0388_06+/Cal_Str_Y_E_MM1.cal
+ ./data/NPRoot/Calibration/CalAlpha_r0388_06+/Cal_Str_Y_E_MM2.cal
+ ./data/NPRoot/Calibration/CalAlpha_r0388_06+/Cal_Str_Y_E_MM3.cal
+ ./data/NPRoot/Calibration/CalAlpha_r0388_06+/Cal_Str_Y_E_MM4.cal
+
+ ./data/NPRoot/Calibration/Time/Coeff/Cal_Str_X_T_MM1_r0183.cal
+ ./data/NPRoot/Calibration/Time/Coeff/Cal_Str_Y_T_MM1_r0184.cal
+ ./data/NPRoot/Calibration/Time/Coeff/Cal_Str_X_T_MM2_r0183.cal
+ ./data/NPRoot/Calibration/Time/Coeff/Cal_Str_Y_T_MM2_r0184.cal
+ ./data/NPRoot/Calibration/Time/Coeff/Cal_Str_X_T_MM3_r0189.cal
+ ./data/NPRoot/Calibration/Time/Coeff/Cal_Str_Y_T_MM3_r0190.cal
+ ./data/NPRoot/Calibration/Time/Coeff/Cal_Str_X_T_MM4_r0187.cal
+ ./data/NPRoot/Calibration/Time/Coeff/Cal_Str_Y_T_MM4_r0188_manual.cal
+
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM1_proton.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM2_proton.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM3_proton.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM4_proton.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM1_deuteron.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM2_deuteron.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM3_deuteron.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM4_deuteron.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM1_triton.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM2_triton.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM3_triton.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM4_triton.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM1_alpha.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM2_alpha.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM3_alpha.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM4_alpha.txt
+
+ ./calibration/CATS/CATS1X/CATSLin1X_r0164.txt
+ ./calibration/CATS/CATS1Y/CATSLin1Y_r0164.txt
+ ./calibration/CATS/CATS2X/CATSLin2X_r0166.txt
+ ./calibration/CATS/CATS2Y/CATSLin2Y_r0166.txt
diff --git a/Projects/E805/Calibration2.txt b/Projects/E805/Calibration2.txt
new file mode 100644
index 000000000..f2ff4cefc
--- /dev/null
+++ b/Projects/E805/Calibration2.txt
@@ -0,0 +1,41 @@
+CalibrationFilePath
+ ./data/NPRoot/Calibration/CalAlpha_r0388_06+/Cal_Str_X_E_MM1.cal
+ ./data/NPRoot/Calibration/CalAlpha_r0388_06+/Cal_Str_X_E_MM2.cal
+ ./data/NPRoot/Calibration/CalAlpha_r0388_06+/Cal_Str_X_E_MM3.cal
+ ./data/NPRoot/Calibration/CalAlpha_r0388_06+/Cal_Str_X_E_MM4.cal
+
+ ./data/NPRoot/Calibration/CalAlpha_r0388_06+/Cal_Str_Y_E_MM1.cal
+ ./data/NPRoot/Calibration/CalAlpha_r0388_06+/Cal_Str_Y_E_MM2.cal
+ ./data/NPRoot/Calibration/CalAlpha_r0388_06+/Cal_Str_Y_E_MM3.cal
+ ./data/NPRoot/Calibration/CalAlpha_r0388_06+/Cal_Str_Y_E_MM4.cal
+
+ ./data/NPRoot/Calibration/Time/Coeff/Cal_Str_X_T_MM1_r0183.cal
+ ./data/NPRoot/Calibration/Time/Coeff/Cal_Str_Y_T_MM1_r0184.cal
+ ./data/NPRoot/Calibration/Time/Coeff/Cal_Str_X_T_MM2_r0183.cal
+ ./data/NPRoot/Calibration/Time/Coeff/Cal_Str_Y_T_MM2_r0184.cal
+ ./data/NPRoot/Calibration/Time/Coeff/Cal_Str_X_T_MM3_r0189.cal
+ ./data/NPRoot/Calibration/Time/Coeff/Cal_Str_Y_T_MM3_r0190.cal
+ ./data/NPRoot/Calibration/Time/Coeff/Cal_Str_X_T_MM4_r0187.cal
+ ./data/NPRoot/Calibration/Time/Coeff/Cal_Str_Y_T_MM4_r0188_manual.cal
+
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM1_proton.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM2_proton.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM3_proton.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM4_proton.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM1_deuteron.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM2_deuteron.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM3_deuteron.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM4_deuteron.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM1_triton.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM2_triton.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM3_triton.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM4_triton.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM1_alpha.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM2_alpha.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM3_alpha.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM4_alpha.txt
+
+ ./calibration/CATS/CATS1X/CATSLin1X_r0366.txt
+ ./calibration/CATS/CATS1Y/CATSLin1Y_r0366.txt
+ ./calibration/CATS/CATS2X/CATSLin2X_r0166.txt
+ ./calibration/CATS/CATS2Y/CATSLin2Y_r0166.txt
diff --git a/Projects/E805/DetectorConfiguration/CATS_MCR.detector b/Projects/E805/DetectorConfiguration/CATS_MCR.detector
index 1b6451b0e..a3d6506d3 100644
--- a/Projects/E805/DetectorConfiguration/CATS_MCR.detector
+++ b/Projects/E805/DetectorConfiguration/CATS_MCR.detector
@@ -11,17 +11,17 @@ Target
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CATSDetector
CATSNumber = 1
- X1_Y1= 43.16 -42.26 -1587.1 mm
- X28_Y1= -27.96 -42.26 -1587.1 mm
- X1_Y28= 43.16 28.86 -1587.1 mm
- X28_Y28= -27.96 28.86 -1587.1 mm
+ X1_Y1= 43.71 -41.76 -1587.1 mm
+ X28_Y1= -27.41 -41.76 -1587.1 mm
+ X1_Y28= 43.71 29.36 -1587.1 mm
+ X28_Y28= -27.41 29.36 -1587.1 mm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CATSDetector
CATSNumber = 2
- X1_Y1= 33.66 -37.06 -1090.1 mm
- X28_Y1= -37.46 -37.06 -1090.1 mm
- X1_Y28= 33.66 34.06 -1090.1 mm
- X28_Y28= -37.46 34.06 -1090.1 mm
+ X1_Y1= 33.94 -37.06 -1090.1 mm
+ X28_Y1= -37.18 -37.06 -1090.1 mm
+ X1_Y28= 33.94 34.06 -1090.1 mm
+ X28_Y28= -37.18 34.06 -1090.1 mm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MASK
MaskNumber = 1
diff --git a/Projects/E805/DetectorConfiguration/MUST2_E805.detector b/Projects/E805/DetectorConfiguration/MUST2_E805.detector
index a1ab2b27b..98663a5d7 100644
--- a/Projects/E805/DetectorConfiguration/MUST2_E805.detector
+++ b/Projects/E805/DetectorConfiguration/MUST2_E805.detector
@@ -7,7 +7,7 @@ Target
ANGLE= 0 deg
X= 0 mm
Y= 0 mm
- Z= 10 mm
+ Z= 0 mm
%%%%%%% Telescope 1 %%%%%%%
M2Telescope
X1_Y1= -13.57 -104.78 299.83 mm
@@ -55,17 +55,17 @@ M2Telescope
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CATSDetector
CATSNumber = 1
- X1_Y1= 43.16 -42.26 -1587.1 mm
- X28_Y1= -27.96 -42.26 -1587.1 mm
- X1_Y28= 43.16 28.86 -1587.1 mm
- X28_Y28= -27.96 28.86 -1587.1 mm
+ X1_Y1= 43.71 -41.76 -1587.1 mm
+ X28_Y1= -27.41 -41.76 -1587.1 mm
+ X1_Y28= 43.71 29.36 -1587.1 mm
+ X28_Y28= -27.41 29.36 -1587.1 mm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CATSDetector
CATSNumber = 2
- X1_Y1= 33.66 -37.06 -1090.1 mm
- X28_Y1= -37.46 -37.06 -1090.1 mm
- X1_Y28= 33.66 34.06 -1090.1 mm
- X28_Y28= -37.46 34.06 -1090.1 mm
+ X1_Y1= 33.94 -37.06 -1090.1 mm
+ X28_Y1= -37.18 -37.06 -1090.1 mm
+ X1_Y28= 33.94 34.06 -1090.1 mm
+ X28_Y28= -37.18 34.06 -1090.1 mm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MASK
MaskNumber = 1
diff --git a/Projects/E805/ProjetInfo/MCRMethod.cpp b/Projects/E805/ProjetInfo/MCRMethod.cpp
index 39a6800f7..5872869b7 100644
--- a/Projects/E805/ProjetInfo/MCRMethod.cpp
+++ b/Projects/E805/ProjetInfo/MCRMethod.cpp
@@ -17,8 +17,9 @@ using namespace std;
MCRMethod::MCRMethod()
{
- Init();
- MetroHast();
+ //Init();
+ // MetroHast();
+ Nappe();
//MinimizationFunction();
}
@@ -185,8 +186,8 @@ void MCRMethod::Init(){
Temperature = TemperatureInitiale;
TemperatureFinale = 0.1;
Step = 0.1;
- runmask[0] = "r0314_algo";
- runmask[1] = "r0315_algo";
+ runmask[0] = "r0367_mask1";
+ runmask[1] = "r0368_mask1";
for(unsigned int i = 0; i < 2; i++){
Chain[i] = new TChain("PhysicsTree");
Chain[i]->Add(path+runmask[i]+".root"); //CATS1
@@ -281,8 +282,131 @@ double MCRMethod::MinimizationFunction(){
return distdiff;
}
+void MCRMethod::Nappe(){
+ int X1 = 40;
+ int X2 = 40;
+ double stepx1 = 0.1;
+ double stepx2 = 0.1;
+ string Side = "Y";
+ TH1Map = new std::map<int,std::map<string,TH1F*>>;
+ double norm[2*X1+1][2*X2+1];
+ for(unsigned int i = 0; i < 2; i++){
+ for(int i1 = -X1; i1 <= X1; i1++){
+ for(int i2 = -X2; i2 <= X2; i2++){
+ (*TH1Map)[i][Form("h_m%i_x1%i_x2%i",i, i1, i2)] = new TH1F(Form("h_m%i_x1%i_x2%i",i, i1, i2),Form("h_m%i_x1%i_x2%i",i, i1, i2),200,-10,10);
+ }
+ }
+ }
+ // std::cout << "test2" << std::endl;
+ runmask[0] = "r0314_mask1";
+ runmask[1] = "r0315_mask1";
+ for(unsigned int i = 0; i < 2; i++){
+ Chain[i] = new TChain("PhysicsTree");
+ Chain[i]->Add(path+runmask[i]+".root"); //CATS1
+ if (!(Chain[i]->GetEntries())){
+ printf("ERROR in the Chain !! \n");
+ return;
+ }
+ // std::cout << "test3" << std::endl;
+ TreeReader[i] = new TTreeReader(Chain[i]);
+ CATSPhysics_[i]= new TTreeReaderValue<TCATSPhysics>(*TreeReader[i],"CATS");
+
+ }
+ TString CName = "CUT_314_mask1";
+ CFile[0] = new TFile(CName+".root");
+ CUT[0] = (TCutG*)CFile[0]->FindObjectAny(CName);
+ std::cout << CUT[0] << std::endl;
+ CName = "CUT_315_mask1";
+ CFile[1] = new TFile(CName+".root");
+ CUT[1] = (TCutG*)CFile[1]->FindObjectAny(CName);
+ std::cout << CUT[1] << std::endl;
+ for(unsigned int i = 0; i < 2; i++){
+ while (TreeReader[i]->Next()){
+ CATSPhysics = &**(CATSPhysics_[i]);
+ if(CATSPhysics->PositionX.size() == 2){
+ // std::cout << "test5 " << CATSPhysics->GetCATSMult() << " " << CATSPhysics->DetNumber[0] << " " << (CATSPhysics->PositionX).size() << " " << CATSPhysics->PositionY[i] << " " << CUT[i]->IsInside(CATSPhysics->PositionX[i],CATSPhysics->PositionY[i]) << std::endl;
+ if(CATSPhysics->DetNumber[0] == 1 && CUT[i]->IsInside(CATSPhysics->PositionX[i],CATSPhysics->PositionY[i])){
+ for(int i1 = -X1; i1 <= X1; i1++){
+ for(int i2 = -X2; i2 <= X2; i2++){
+ // std::cout << "test " << i1 << " " << i2 << std::endl;
+ // std::cout << "test4" << std::endl;
+ double x1, x2;
+ if(Side == "X"){
+ x1 = CATSPhysics->PositionX[0] +i1*stepx1;
+ x2 = CATSPhysics->PositionX[1] +i2*stepx2;
+ }
+ else if(Side == "Y"){
+ x1 = CATSPhysics->PositionY[0] +i1*stepx1;
+ x2 = CATSPhysics->PositionY[1] +i2*stepx2;
+ }
+ (*TH1Map)[i][Form("h_m%i_x1%i_x2%i",i, i1, i2)]->Fill(ProjectOnCats(i,x1,x2));
+ norm[X1+i1][X2+i2] = 0;
+
+ }
+ }
+ }
+ }
+ }
+ }
+
+ (*TH1Map)[0][Form("h_m%i_x1%i_x2%i",0, 0, 0)]->Draw();
+ (*TH1Map)[0][Form("h_m%i_x1%i_x2%i",0, 1, 1)]->Draw("same");
+ (*TH1Map)[0][Form("h_m%i_x1%i_x2%i",0, 2, 2)]->Draw("same");
+ auto c1 = new TCanvas;
+ TH2F* HeatMap[2];
+ auto fitfunc = new TF1("fitfunc","gausn",-10,10);
+ auto SumHeatMap = new TH2F("Sumheatmap","Sumheatmap",2*X1+1, -X1*stepx1, (X1+1)*stepx1, 2*X2+1, -X2*stepx2, (X2+1)*stepx2);
+ HeatMap[0] = new TH2F("heatmap_0","heatmap_0",2*X1+1, -X1*stepx1, (X1+1)*stepx1, 2*X2+1, -X2*stepx2, (X2+1)*stepx2);
+ HeatMap[1] = new TH2F("heatmap_1","heatmap_1",2*X1+1, -X1*stepx1, (X1+1)*stepx1, 2*X2+1, -X2*stepx2, (X2+1)*stepx2);
+ TSpectrum* Spec = new TSpectrum(1,1.);
+ for(unsigned int i = 0; i < 2; i++){
+ for(int i1 = -X1; i1 <= X1; i1++){
+ for(int i2 = -X2; i2 <= X2; i2++){
+ Int_t nfound = Spec->Search((*TH1Map)[i][Form("h_m%i_x1%i_x2%i",i, i1, i2)],2,"",0.15);
+ Double_t* PeaksPosX = Spec->GetPositionX();
+ fitfunc->SetParameters(200,PeaksPosX[0],1);
+ fitfunc->SetParLimits(1,-10,10);
+ (*TH1Map)[i][Form("h_m%i_x1%i_x2%i",i, i1, i2)]->Fit(fitfunc,"QL","",-10,10);
+ // std::cout << i1 << " " << i2 << " " << PeaksPosX[0] << std::endl;
+ double value;
+ if(Side == "X")
+ value = abs(fitfunc->GetParameter(1) - MaskPosX[i]);
+ else if(Side == "Y")
+ value = abs(fitfunc->GetParameter(1) - MaskPosY[i]);
+ HeatMap[i]->SetBinContent(X1+i1+1,X2+i2+1,value);
+ //HeatMap[i]->SetBinContent(i1+1,i2+1,value);
+ std::cout << norm[X1+i1][X2+i2] << std::endl;
+ norm[X1+i1][X2+i2] += value*value;
+ if(i ==1){
+ norm[X1+i1][X2+i2] = sqrt(norm[X1+i1][X2+i2]);
+
+ }
+ }
+ }
+ }
+ for(int i1 = -X1; i1 <= X1; i1++){
+ for(int i2 = -X2; i2 <= X2; i2++){
+ SumHeatMap->SetBinContent(X1+i1+1,X2+i2+1,norm[X1+i1][X2+i2]);
+
+ }
+ }
+ // SumHeatMap->Add(HeatMap[0],HeatMap[1],1,1);
+ for(unsigned int i = 0; i < 2; i++){
+ auto c = new TCanvas();
+ HeatMap[i]->Draw("colz");
+ }
+ auto c = new TCanvas();
+ SumHeatMap->Draw("colz");
+}
+
+
+float MCRMethod::ProjectOnCats(unsigned int i, double x1, double x2){
+ double tmask = (CATSPosZ[1]-MASKPosZ[i])/(CATSPosZ[0] - MASKPosZ[i]);
+ return x2 -x1*tmask;
+};
+
diff --git a/Projects/E805/ProjetInfo/MCRMethod.h b/Projects/E805/ProjetInfo/MCRMethod.h
index e6775df4c..3271d6cec 100644
--- a/Projects/E805/ProjetInfo/MCRMethod.h
+++ b/Projects/E805/ProjetInfo/MCRMethod.h
@@ -12,6 +12,9 @@
#include "TTreeReaderValue.h"
#include "TH2F.h"
#include "TSpectrum2.h"
+#include "TSpectrum.h"
+#include "TCutG.h"
+#include "TF1.h"
// NPTOOL
#include "TCATSPhysics.h"
@@ -32,6 +35,10 @@ class MCRMethod
void RandomStep();
double MinimizationFunction();
void Init();
+ float ProjectOnCats(unsigned int i, double x1, double x2);
+
+
+ void Nappe();
private:
@@ -55,9 +62,11 @@ class MCRMethod
private:
TChain *Chain[2];
TTreeReader *TreeReader[2];
- TString path = "./NPRootA/";
+ TString path = "./RootR/";
TTreeReaderValue<TCATSPhysics> *CATSPhysics_[2];
- TCATSPhysics *CATSPhysics;
+ TCATSPhysics *CATSPhysics;
+ TCutG *CUT[2];
+ TFile *CFile[2];
Double_t TopLeft[2][2];
Double_t BotLeft[2][2];
Double_t BotRight[2][2];
@@ -66,8 +75,12 @@ class MCRMethod
Double_t BotLeftGeo[4] = {-12.5,-12.3,-12.4,-13.6};
Double_t BotRightGeo[4] = {12.5,-12.3,12.6,-13.6};
+ std::map<int,std::map<string,TH1F*>>* TH1Map;
TGraph* Graph;
+ double MaskPosX[2] = {0,0.1};
+ double MaskPosY[2] = {0.2,-1.1};
+
double tini;
double tfinale;
double pas;
--
GitLab