diff --git a/NPLib/Detectors/BigRIPS/CMakeLists.txt b/NPLib/Detectors/BigRIPS/CMakeLists.txt
index d93a550dbabee40d7da150026b46f626d999b162..779a57e863920576dd30d481b0023c26656f7484 100644
--- a/NPLib/Detectors/BigRIPS/CMakeLists.txt
+++ b/NPLib/Detectors/BigRIPS/CMakeLists.txt
@@ -10,8 +10,12 @@ add_custom_command(OUTPUT TBigRIPSICDataDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scr
add_custom_command(OUTPUT TBigRIPSICPhysicsDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TBigRIPSICPhysics.h TBigRIPSICPhysicsDict.cxx TBigRIPSICPhysics.rootmap libNPBigRIPS.dylib DEPENDS TBigRIPSICPhysics.h)
-add_library(NPBigRIPS SHARED TBigRIPSPPACData.cxx TBigRIPSPPACDataDict.cxx TBigRIPSPPACPhysics.cxx TBigRIPSPPACPhysicsDict.cxx TBigRIPSPlasticData.cxx TBigRIPSPlasticDataDict.cxx TBigRIPSPlasticPhysics.cxx TBigRIPSPlasticPhysicsDict.cxx TBigRIPSICData.cxx TBigRIPSICDataDict.cxx TBigRIPSICPhysics.cxx TBigRIPSICPhysicsDict.cxx)
+add_custom_command(OUTPUT TBigRIPSFocalDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TBigRIPSFocal.h TBigRIPSFocalDict.cxx TBigRIPSFocal.rootmap libNPBigRIPS.dylib DEPENDS TBigRIPSFocal.h)
+
+add_custom_command(OUTPUT TBigRIPSRecoDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TBigRIPSReco.h TBigRIPSRecoDict.cxx TBigRIPSReco.rootmap libNPBigRIPS.dylib DEPENDS TBigRIPSReco.h)
+
+add_library(NPBigRIPS SHARED TBigRIPSPPACData.cxx TBigRIPSPPACDataDict.cxx TBigRIPSPPACPhysics.cxx TBigRIPSPPACPhysicsDict.cxx TBigRIPSPlasticData.cxx TBigRIPSPlasticDataDict.cxx TBigRIPSPlasticPhysics.cxx TBigRIPSPlasticPhysicsDict.cxx TBigRIPSICData.cxx TBigRIPSICDataDict.cxx TBigRIPSICPhysics.cxx TBigRIPSICPhysicsDict.cxx TBigRIPSFocal.cxx TBigRIPSFocalDict.cxx TBigRIPSReco.cxx TBigRIPSRecoDict.cxx)
target_link_libraries(NPBigRIPS ${ROOT_LIBRARIES} NPCore NPTrackReconstruction)
-install(FILES TBigRIPSPPACData.h TBigRIPSPPACPhysics.h TBigRIPSPlasticData.h TBigRIPSPlasticPhysics.h TBigRIPSICData.h TBigRIPSICPhysics.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
+install(FILES TBigRIPSPPACData.h TBigRIPSPPACPhysics.h TBigRIPSPlasticData.h TBigRIPSPlasticPhysics.h TBigRIPSICData.h TBigRIPSICPhysics.h TBigRIPSFocal.h TBigRIPSReco.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
diff --git a/NPLib/Detectors/BigRIPS/TBigRIPSFocal.cxx b/NPLib/Detectors/BigRIPS/TBigRIPSFocal.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..63df7c5ee26a44e8542fee6949e2a667239afc5a
--- /dev/null
+++ b/NPLib/Detectors/BigRIPS/TBigRIPSFocal.cxx
@@ -0,0 +1,26 @@
+#include "TBigRIPSFocal.h"
+
+TBigRIPSFocal::TBigRIPSFocal(){};
+TBigRIPSFocal::~TBigRIPSFocal(){Clear();};
+
+
+////////////////////////////////////////////////////////////////////////////////
+void TBigRIPSFocal::Clear(){
+ FPTrack.clear();
+}
+////////////////////////////////////////////////////////////////////////////////
+void TBigRIPSFocal::Init(int N_FP){
+ std::vector<double> tmp;
+ for(int j=0;j<4;j++) tmp.push_back(-9999);
+ for(int i=0;i<N_FP;i++) FPTrack.push_back(tmp);
+}
+////////////////////////////////////////////////////////////////////////////////
+void TBigRIPSFocal::Print(int FPNbr){
+ std::cout << "----FP----:" << FPNbr <<std::endl;
+ std::cout << "X: " << FPTrack[FPNbr][0] << std::endl;
+ std::cout << "A: " << FPTrack[FPNbr][1] << std::endl;
+ std::cout << "Y: " << FPTrack[FPNbr][2] << std::endl;
+ std::cout << "B: " << FPTrack[FPNbr][3] << std::endl;
+
+}
+ClassImp(TBigRIPSFocal);
diff --git a/NPLib/Detectors/BigRIPS/TBigRIPSFocal.h b/NPLib/Detectors/BigRIPS/TBigRIPSFocal.h
new file mode 100644
index 0000000000000000000000000000000000000000..33531833b3de2d3b6259261437fb6824b74e2d97
--- /dev/null
+++ b/NPLib/Detectors/BigRIPS/TBigRIPSFocal.h
@@ -0,0 +1,35 @@
+#ifndef TBigRIPSFocal_H
+#define TBigRIPSFocal_H
+#include "TObject.h"
+#include <vector>
+#include <iostream>
+class TBigRIPSFocal: public TObject{
+
+ public:
+ TBigRIPSFocal();
+ ~TBigRIPSFocal();
+
+ private:
+ std::vector<std::vector<double>> FPTrack ;
+
+ public:
+ void Clear();
+ void Init(int);
+ void Print(int);
+
+ public:
+ void SetFPTrack(int FPNbr, std::vector<double> track){
+ FPTrack[FPNbr]=track;
+ };
+ void SetFPTrack(int FPNbr, double x, double a, double y, double b){
+ FPTrack[FPNbr][0]=x;
+ FPTrack[FPNbr][1]=a;
+ FPTrack[FPNbr][2]=y;
+ FPTrack[FPNbr][3]=b;
+ };
+ std::vector<double> GetFPTrack(int FPNbr) {return FPTrack[FPNbr];}
+
+ ClassDef(TBigRIPSFocal,1);
+};
+
+#endif
diff --git a/NPLib/Detectors/BigRIPS/TBigRIPSPPACData.cxx b/NPLib/Detectors/BigRIPS/TBigRIPSPPACData.cxx
index e64b1834de0cc34b6c2939d897efa3ed18466911..62041acff087398f571c4321df3bfa9f7102fcb9 100644
--- a/NPLib/Detectors/BigRIPS/TBigRIPSPPACData.cxx
+++ b/NPLib/Detectors/BigRIPS/TBigRIPSPPACData.cxx
@@ -4,30 +4,8 @@
TBigRIPSPPACData::TBigRIPSPPACData(){};
TBigRIPSPPACData::~TBigRIPSPPACData(){};
-/*
-void TBigRIPSPPACData::SetData(const int& FP, const int& ID, const int& Edge, const int& value,const int& VariableType){
- fPPAC_FP.push_back(FP);
- fPPAC_ID.push_back(ID);
- fPPAC_Edge.push_back(Edge);
- switch(VariableType){
- case 0: fPPAC_TX1.push_back(value); break;
- case 1: fPPAC_TX2.push_back(value); break;
- case 2: fPPAC_TY1.push_back(value); break;
- case 3: fPPAC_TY2.push_back(value); break;
- case 4: fPPAC_TA.push_back(value); break;
- case 5: fPPAC_QX1.push_back(value); break;
- case 6: fPPAC_QX2.push_back(value); break;
- case 7: fPPAC_QY1.push_back(value); break;
- case 8: fPPAC_QY2.push_back(value); break;
- case 9: fPPAC_QA.push_back(value); break;
- }
-}
-*/
////////////////////////////////////////////////////////////////////////////////
void TBigRIPSPPACData::Clear(){
- // fPPAC_FP.clear();
- // fPPAC_ID.clear();
- // fPPAC_Edge.clear();
fPPAC_TX1.clear();
fPPAC_TX2.clear();
fPPAC_TY1.clear();
@@ -59,32 +37,4 @@ void TBigRIPSPPACData::Print(){
//cout << " - Multiplicity: " << Mult() << endl;
}
-////////////////////////////////////////////////////////////////////////////////
-/*
-unsigned int TBigRIPSPPACData::MultLayer(unsigned int det , unsigned int layer, int edge){
- unsigned int mult=0;
- unsigned int size = fDC_DetectorNbr.size();
- for(unsigned int i = 0 ; i< size ; i++ ){
- if(fDC_DetectorNbr[i]==det)
- if(fDC_LayerNbr[i]==layer)
- if(fDC_Edge[i]==edge && edge!=-1) // edge type is specified (0 or 1)
- mult++;
- else if(edge==-1)// edge type is not specified
- mult++;
- }
- return mult;
-
-}
-////////////////////////////////////////////////////////////////////////////////
-std::vector<int> TBigRIPSPPACData::GetWire(unsigned int det , unsigned int layer){
- std::vector<int> wires;
- unsigned int size = fDC_DetectorNbr.size();
- for(unsigned int i = 0 ; i< size ; i++ ){
- if(fDC_DetectorNbr[i]==det)
- if(fDC_LayerNbr[i]==layer)
- wires.push_back(fDC_WireNbr[i]);
- }
- return wires;
-}
-*/
ClassImp(TBigRIPSPPACData);
diff --git a/NPLib/Detectors/BigRIPS/TBigRIPSReco.cxx b/NPLib/Detectors/BigRIPS/TBigRIPSReco.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..8b9c3219e143e896946d5cc53e31918d3342231c
--- /dev/null
+++ b/NPLib/Detectors/BigRIPS/TBigRIPSReco.cxx
@@ -0,0 +1,84 @@
+#include "TBigRIPSReco.h"
+
+TBigRIPSReco::TBigRIPSReco(){Init();}
+TBigRIPSReco::~TBigRIPSReco(){Clear();}
+
+
+////////////////////////////////////////////////////////////////////////////////
+void TBigRIPSReco::Clear(){
+}
+////////////////////////////////////////////////////////////////////////////////
+void TBigRIPSReco::Init(){
+ aoq=-9999;
+ beta=-9999;
+ delta=-9999;
+ brho=-9999;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+
+void TBigRIPSReco::RecBrho(std::vector<double> RecFPUpstream,std::vector<double> RecFPDownstream, std::vector<std::vector<double>> matrix,double BrhoCentral){
+
+ if(matrix.size()!=6) {
+ std::cout << "MATRIX used for Brho Rec HAS NOT 6 rows" << std::endl;
+ return;
+ }else{
+ for(int i=0;i<6;i++){
+ if(matrix[i].size()!=5){
+ std::cout << "MATRIX used for Brho Rec. HAS NOT 5 columns" << std::endl;
+ std::cout << "but" << matrix[i].size()<<std::endl;
+ return;
+ }
+ }
+ }
+
+ if(RecFPUpstream[0]==-9999 ||
+ RecFPUpstream[1]==-9999 ||
+ RecFPDownstream[0]==-9999)
+ return;
+
+ double x_x = matrix[0][0];
+ double a_x = matrix[0][1];
+ double x_a = matrix[1][0];
+ double x_d = matrix[5][0];
+ double a_a = matrix[1][1];
+ double a_d = matrix[5][1];
+
+ double x_up = RecFPUpstream[0]; // downst. X
+ double a_up = RecFPUpstream[1]; // downst. X
+ double x_down = RecFPDownstream[0]; // downst. X
+ double a_down = RecFPDownstream[1]; // downst. A
+
+ delta = (a_a * x_down - x_a * a_down - x_up) / (x_d * a_a - x_a * a_d);
+ angle = (a_a * x_down - x_a * a_down - x_up) / (x_d * a_a - x_a * a_d);
+
+ brho = BrhoCentral*(1.0+delta*0.01);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void TBigRIPSReco::RecAoqOne(double tof, double length){
+ beta = length /(tof * c_mm_ns);
+ double gamma = 1./sqrt(1 - beta*beta);
+ aoq = (brho * c_mm_ns) / (mnucleon * beta * gamma);
+//std::cout << aoq << std::endl;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void TBigRIPSReco::RecZet(double dE,double ionpair, std::vector<double> coeff){
+ if(coeff.size()<2){
+ std::cout << "Zcoeff for Z reco not defined" << std::endl;
+ return;
+ }
+ double de_v = log(ionpair*beta*beta) - log((1-beta*beta)) - beta*beta;
+ z = coeff[0] * sqrt(dE/de_v)*beta + coeff[1];
+}
+////////////////////////////////////////////////////////////////////////////////
+void TBigRIPSReco::Print(){
+ std::cout << "aoq: " << aoq << std::endl;
+ std::cout << "z: " << z << std::endl;
+ std::cout << "beta: " << beta << std::endl;
+ std::cout << "delta: "<< delta << std::endl;
+ std::cout << "brho: " << brho << std::endl;
+}
+
+ClassImp(TBigRIPSReco);
diff --git a/NPLib/Detectors/BigRIPS/TBigRIPSReco.h b/NPLib/Detectors/BigRIPS/TBigRIPSReco.h
new file mode 100644
index 0000000000000000000000000000000000000000..023f9f906ad52cc3e68e78b84cf1826c6e5ec728
--- /dev/null
+++ b/NPLib/Detectors/BigRIPS/TBigRIPSReco.h
@@ -0,0 +1,39 @@
+#ifndef TBigRIPSReco_H
+#define TBigRIPSReco_H
+#include "TObject.h"
+#include <vector>
+#include <iostream>
+class TBigRIPSReco: public TObject{
+
+ public:
+ TBigRIPSReco();
+ ~TBigRIPSReco();
+
+ private:
+ double aoq;
+ double beta;
+ double delta;
+ double brho;
+ double angle;
+ double z;
+ public:
+ void Clear();
+ void Init();
+ void Print();
+ void RecBrho(std::vector<double>,std::vector<double>,std::vector<std::vector<double>>,double);
+ void RecAoqOne(double, double);
+ void RecZet(double,double, std::vector<double>);
+ const double c_mm_ns = 299.7792458; // in mm/ns
+ const double mnucleon = 931.49432; // MeV
+
+ public:
+ void SetAoq(double value){aoq=value;}
+ void SetBeta(double value){beta=value;}
+ void SetDelta(double value){delta=value;}
+ void SetBrho(double value){brho=value;}
+ void SetZ(double value){z=value;}
+
+ ClassDef(TBigRIPSReco,1);
+};
+
+#endif
diff --git a/Projects/RIBF196/Analysis.cxx b/Projects/RIBF196/Analysis.cxx
index 04cb90ffa65d57c9cdf08613eb1f92199ea49283..5da4c57c645963f8af529ecdbc7b724c93117434 100644
--- a/Projects/RIBF196/Analysis.cxx
+++ b/Projects/RIBF196/Analysis.cxx
@@ -40,6 +40,14 @@ Analysis::~Analysis(){
////////////////////////////////////////////////////////////////////////////////
void Analysis::Init() {
+ NmaxFP=11;
+ //FP = new TBigRIPSFocal();
+ //FP->Init(NmaxFP+1);
+ Rec35 = new TBigRIPSReco();
+ Rec57 = new TBigRIPSReco();
+ Rec89 = new TBigRIPSReco();
+ Rec911 = new TBigRIPSReco();
+
InitOutputBranch();
InitInputBranch();
// get PPAC, PL, and IC objects
@@ -48,9 +56,25 @@ void Analysis::Init() {
IC = (TBigRIPSICPhysics*) m_DetectorManager -> GetDetector("BigRIPSIC");
ReadXmls();
-
+/*
matF35.ResizeTo(6,5); matF35 = RecReadTransferMatrix("mat1.mat");
matF57.ResizeTo(6,5); matF57 = RecReadTransferMatrix("mat2.mat");
+ matF89.ResizeTo(6,5); matF89 = RecReadTransferMatrix("F8F9_LargeAccAchr.mat");
+ matF911.ResizeTo(6,5); matF911 = RecReadTransferMatrix("F9F11_LargeAccAchr.mat");
+*/
+ matrixF35 = RecReadTransferMatrix2("mat1.mat");
+ matrixF57 = RecReadTransferMatrix2("mat2.mat");
+ matrixF89 = RecReadTransferMatrix2("F8F9_LargeAccAchr.mat");
+ matrixF911 = RecReadTransferMatrix2("F9F11_LargeAccAchr.mat");
+
+ length37 = (FP_Z[FPtoID[7]]+PL_NametoZ["F7pl"])
+ -(FP_Z[FPtoID[3]]+PL_NametoZ["F3pl"]);
+ length35 = FP_Z[FPtoID[5]] - (FP_Z[FPtoID[3]]+PL_NametoZ["F3pl"]);
+ length57 = (FP_Z[FPtoID[7]]+PL_NametoZ["F7pl"]) - FP_Z[FPtoID[5]];
+ length811 = (FP_Z[FPtoID[11]]+PL_NametoZ["F11pl-1"])
+ -(FP_Z[FPtoID[8]]+PL_NametoZ["F8pl"]);
+ length89 = FP_Z[FPtoID[9]] - (FP_Z[FPtoID[8]]+PL_NametoZ["F8pl"]);
+ length911 = (FP_Z[FPtoID[11]]+PL_NametoZ["F11pl-1"]) - FP_Z[FPtoID[9]];
// initialize various parameters
//Rand = TRandom3();
@@ -111,11 +135,15 @@ void Analysis::ReadXmls() {
IC_IDtoName[ID] = bic[i]->AsString("NAME");
string name = bic[i]->AsString("NAME");
IC_NametoID[name] = ID;
+ IC_Zcoef[name].push_back(bic[i]->AsDouble("zcoef_0"));
+ IC_Zcoef[name].push_back(bic[i]->AsDouble("zcoef_1"));
+ IC_Zcoef[name].push_back(bic[i]->AsDouble("zcoef_2"));
+ IC_Ionpair[name] = bic[i]->AsDouble("ionpair");
}
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::TreatEvent() {
- // Reinitiate calculated variable
+ // Reinitiate and clear variables
ReInitValue();
FP_PPACHitList.clear();
@@ -124,162 +152,106 @@ void Analysis::TreatEvent() {
/////////////////////////////////////////////
// Get all PPAC X,Y hits and store them in a HitList indexed by Focalplane
+ // For a given Focal plane, the HitList size depends on how many PPAC
+ // are installed at this FP
for(int i=0;i<PPAC->ID.size();i++){
FP_PPACHitList[PPAC->FP[i]].AddPPACHit(PPAC->X[i],PPAC->Y[i],PPAC->ID[i]);
}
- // Perform FP reco. and store in map: index FP, content vector(x,a,y,b)
- std::map<int,vector<double>> RecFP;
- double NmaxFP=12;
- for(int i=0;i<NmaxFP;i++){
- for(int j=0;j<4;j++){
- RecFP[i].push_back(-9999);}}
-
+ // For each FP, perform reco. and store in vector<vector>
+ // RecFP[FPNbr][0,1,2,3 for x,a,y,b]
+ for(auto it = FP_PPACHitList.begin();it!=FP_PPACHitList.end();++it){
+ if(it->first < NmaxFP){
+ RecFP[it->first] = RecFPposition(it->second.PPACHit,it->second.PPACID);
+ }
+ }
+ // Same as above but using an object from class TBigRIPSFocal
+ /*
for(auto it = FP_PPACHitList.begin();it!=FP_PPACHitList.end();++it){
- RecFP[it->first] = RecFPposition(it->second.PPACHit,it->second.PPACID);
+ if(it->first < NmaxFP){
+ FP->SetFPTrack(it->first, RecFPposition(it->second.PPACHit,it->second.PPACID) );
+ //FP->Print(it->first);
+ }
}
+*/
- fX=RecFP[8][0];
- fA=RecFP[8][1];
- fY=RecFP[8][2];
- fB=RecFP[8][3];
-
- delta35 = RecDeltaBrho(RecFP[3],RecFP[5],matF35);
- Brho35 = BrhoD3 * (1.0 + delta35*0.01);
+ // Brho/delta reconstruction
+ Rec35->RecBrho(RecFP[3],RecFP[5],matrixF35,BrhoD3);
+ Rec57->RecBrho(RecFP[5],RecFP[7],matrixF57,BrhoD5);
+ Rec89->RecBrho(RecFP[8],RecFP[9],matrixF89,BrhoD7);
+ Rec911->RecBrho(RecFP[9],RecFP[11],matrixF911,BrhoD8);
- delta57 = RecDeltaBrho(RecFP[5],RecFP[7],matF57);
- Brho57 = BrhoD5 * (1.0 + delta57*0.01);
/////////////////////////////////
/// STEP2: TOF, flight length ///
/////////////////////////////////
std::map<string,double> PL_T; //index is plastic name in xml
- for(int i=0;i<PL->ID.size();i++) {PL_T[PL_IDtoName[PL->ID[i]]] = PL->T[i];
- //if(PL->ID[i]==2) std::cout << "PL ID2 T::" << PL->T[i] << std::endl;
- }
+ for(int i=0;i<PL->ID.size();i++) {PL_T[PL_IDtoName[PL->ID[i]]] = PL->T[i];}
tf3 = PL_T["F3pl"];
tf7 = PL_T["F7pl"];
- if(PL_T["F7pl"]!=-9999 && PL_T["F3pl"]!=-9999) tof37 = PL_T["F7pl"] - PL_T["F3pl"];
+ if(PL_T["F7pl"]!=-9999 && PL_T["F3pl"]!=-9999) tof37 = PL_T["F7pl"] - PL_T["F3pl"];
tof37 += tof_offset37;
-
- //double length35 = (FP_Z[5]+PL_NametoZ["F5pl"]) - (FP_Z[3]+PL_NametoZ["F3pl"]) ;
- //double length57 = (FP_Z[7]+PL_NametoZ["F7pl"]) - (FP_Z[5]+PL_NametoZ["F5pl"]) ;
- //double length37 = length35 + length57 ;
-
- double length37 = (FP_Z[FPtoID[7]]+PL_NametoZ["F7pl"]) - (FP_Z[FPtoID[3]]+PL_NametoZ["F3pl"]);
- double length35 = FP_Z[FPtoID[5]] - (FP_Z[FPtoID[3]]+PL_NametoZ["F3pl"]);
- double length57 = (FP_Z[FPtoID[7]]+PL_NametoZ["F7pl"]) - FP_Z[FPtoID[5]];
-
- //std::cout << "length37:" << length37 << std::endl;
- //std::cout << "Z FP7:" << FP_Z[FPtoID[7]] << std::endl;
- //std::cout << "Zoff PlasticF7:" << PL_NametoZ["F7pl"] << std::endl;
- //std::cout << "Z FP3:" << FP_Z[FPtoID[3]] << std::endl;
- //std::cout << "Zoff PlasticF3:" << PL_NametoZ["F3pl"] << std::endl;
- beta35 = length35 /(tof37/2. * clight);
- beta57 = length57 /(tof37/2. * clight);
- beta37 = length37 /(tof37 * clight);
+ tf8 = PL_T["F8pl"];
+ tf11 = PL_T["F11pl-1"];
+ if(PL_T["F11pl-1"]!=-9999 && PL_T["F8pl"]!=-9999) tof811 = PL_T["F11pl-1"] - PL_T["F8pl"];
+ tof811 += tof_offset811;
+
/////////////////////////////////
/// STEP3: Beta, Gamma, AoQ ///
/////////////////////////////////
- // std::cout << "______________________" << std::endl;
- // std::cout << "Brho35:" << Brho35 << std::endl;
- // std::cout << "Tof37:" << tof37 << std::endl;
- // std::cout << "length35:" << length35 << std::endl;
- aoq35 = RecAoqOneFold(Brho35, tof37, length37);
- // std::cout << "aoq35:" << aoq35 << std::endl;
-
- // std::cout << "______________________" << std::endl;
- // std::cout << "Brho57:" << Brho57 << std::endl;
- // std::cout << "Tof37:" << tof37 << std::endl;
- // std::cout << "length57:" << length57 << std::endl;
- aoq57 = RecAoqOneFold(Brho57, tof37, length37);
- // std::cout << "aoq57:" << aoq57 << std::endl;
+ // Calculate Beta from length/tof and combine it with Brho to get AoQ
+ Rec35->RecAoqOne(tof37,length37);
+ Rec57->RecAoqOne(tof37,length37);
+ Rec89->RecAoqOne(tof811,length811);
+ Rec911->RecAoqOne(tof811,length811);
- //////////////////
- /// STEP4: Z ///
- //////////////////
+ ///////////////////////////////////////////////
+ /// STEP4: Z reco. from dE in IC and Beta ///
+ ///////////////////////////////////////////////
std::map<string,double> IC_dE; //index is IC name in xml
- for(int i=0;i<IC->ID.size();i++) {IC_dE[IC_IDtoName[IC->ID[i]]] = IC->CalSqSum[i];
+ for(int i=0;i<IC->ID.size();i++) {
+ IC_dE[IC_IDtoName[IC->ID[i]]] = IC->CalSqSum[i];
//if(IC->ID[i]==1) std::cout << "IC ID1 CalSq:" << IC->CalSqSum[i] << std::endl;
}
dE_ICF7 = IC_dE["F7IC"];
- z_BR = RecZ(dE_ICF7, tof37, length37);
+ Rec35->RecZet(dE_ICF7,IC_Ionpair["F7IC"],IC_Zcoef["F7IC"]);
+ Rec57->RecZet(dE_ICF7,IC_Ionpair["F7IC"],IC_Zcoef["F7IC"]);
+
+ dE_ICF11 = IC_dE["F11IC"];
+ Rec89->RecZet(dE_ICF11,IC_Ionpair["F11IC"],IC_Zcoef["F11IC"]);
+ Rec911->RecZet(dE_ICF11,IC_Ionpair["F11IC"],IC_Zcoef["F11IC"]);
}
////////////////////////////////////////////////////////////////////////////////
-TMatrixD Analysis::RecReadTransferMatrix(string matfile){
+std::vector<std::vector<double>> Analysis::RecReadTransferMatrix2(string matfile){
char buffer[256];
std::ifstream matfin(matfile);
if(matfin.fail()) cout<<"Failed to open optical matrix file: "<<matfile<<endl;
matfin.getline(buffer,256);
- TMatrixD matrix(6,5);
+ std::vector<std::vector<double>> matrix;
+ //matrix.resize(5, vector<double>(6));
+
double val[6];
+ std::vector<double> row;
for(Int_t i=0;i<6;i++){
+ row.clear();
matfin>>val[0]>>val[1]>>val[2]>>val[3]>>val[4]>>val[5];
//cout<<val[0]<<val[1]<<val[2]<<val[3]<<val[4]<<val[5]<<endl;
- for(Int_t j=0;j<5;j++)
- matrix(i,j) = val[j];
+ for(Int_t j=0;j<5;j++) row.push_back(val[j]);
+ matrix.push_back(row);
}
matfin.close();
-
- matrix.Print();
return matrix;
}
-////////////////////////////////////////////////////////////////////////////////
-double Analysis::RecDeltaBrho(std::vector<double> RecFPUpstream,std::vector<double> RecFPDownstream, TMatrixD matrix){
-
- if(matrix.GetNrows()!=6 && matrix.GetNcols()!=5){
- std::cout << "MATRIX used for Brho Rec. is NULL" << std::endl;
- return -9999;
- }
-
- double x_a = matrix(0,0);
-
- TMatrixD mat1(2,1);
- mat1(0,0) = matrix(0,0); mat1(1,0) = matrix(0,1);
- //mat1.Print();
- TMatrixD mat2(2,2);
- mat2(0,0) = matrix(1,0); mat2(0,1) = matrix(5,0);
- mat2(1,0) = matrix(1,1); mat2(1,1) = matrix(5,1);
- //mat2.Print();
- mat2.Invert();
-
- //TVectorD * fvec = (TVectorD *)fUpstreamFplArrayBuffer[i]->GetOptVector();
- //TVectorD * bvec = (TVectorD *)fDownstreamFplArrayBuffer[i]->GetOptVector();
-
- if(RecFPUpstream[0]==-9999 ||
- RecFPUpstream[1]==-9999 ||
- RecFPDownstream[0]==-9999)
- return -9999;
-
- TMatrixD xvec(2,1);
- //xvec(0,0) = (*bvec)(0);
- //xvec(1,0) = (*bvec)(1);
- xvec(0,0) = RecFPDownstream[0]; // downst. X
- xvec(1,0) = RecFPDownstream[1]; // downst. A
-
- TMatrixD rvec = mat2 * (xvec - RecFPUpstream[0] * mat1);
- Double_t angle = rvec(0,0); // change to mrad
- Double_t delta = rvec(1,0);
- // simple_delta = (*bvec)(0) / 33.; // only for fpl-5
-
- //double BrhoRec = BrhoCentral*(1.0+delta*0.01);
- //rips->SetDelta(delta);
- //rips->SetAngle(angle);
-
-
- //fReconstructed = true;
- return delta;
-}
////////////////////////////////////////////////////////////////////////////////
// Function transforming a list of PPAC Hit positions and IDs for a given FP
// into positions and angles at the FP position (FX,FY,FA,FB)
@@ -376,47 +348,29 @@ std::vector<double> Analysis::RecFPposition(std::vector<TVector2> HitList,std::v
}
-////////////////////////////////////////////////////////////////////////////////
-double Analysis::RecAoqOneFold(double Brho, double tof, double length){
- double beta = length /(tof * clight);
- double gamma = 1./sqrt(1 - beta*beta);
- double aoq = (Brho * clight) / (mnucleon * beta * gamma);
- return aoq;
-}
-
-////////////////////////////////////////////////////////////////////////////////
-double Analysis::RecZ(double dE, double tof, double length){
- double ionpair = 4866.;
- double beta = length /(tof * clight);
- Double_t de_v = log(ionpair*beta*beta) - log((1-beta*beta)) - beta*beta;
- double z = 17.9143234533*sqrt(dE/de_v)*beta -13.0990490522;
- return z;
-}
////////////////////////////////////////////////////////////////////////////////
void Analysis::End(){
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitOutputBranch() {
- RootOutput::getInstance()->GetTree()->Branch("EventNumber",&EventNumber,"EventNumber/I");
RootOutput::getInstance()->GetTree()->Branch("RunNumber",&RunNumber,"RunNumber/I");
- RootOutput::getInstance()->GetTree()->Branch("fX",&fX,"fX/D");
- RootOutput::getInstance()->GetTree()->Branch("fY",&fY,"fY/D");
- RootOutput::getInstance()->GetTree()->Branch("delta35",&delta35,"delta35/D");
- RootOutput::getInstance()->GetTree()->Branch("delta57",&delta57,"delta57/D");
- RootOutput::getInstance()->GetTree()->Branch("Brho35",&Brho35,"Brho35/D");
- RootOutput::getInstance()->GetTree()->Branch("Brho57",&Brho57,"Brho57/D");
+ RootOutput::getInstance()->GetTree()->Branch("EventNumber",&EventNumber,"EventNumber/I");
+ RootOutput::getInstance()->GetTree()->Branch("Trigger",&Trigger,"Trigger/I");
+ RootOutput::getInstance()->GetTree()->Branch("TimeStamp",&TimeStamp,"TimeStamp/I");
+ RootOutput::getInstance()->GetTree()->Branch("RecFP",&RecFP);
+ //RootOutput::getInstance()->GetTree()->Branch("FP","TBigRIPSFocal",&FP);
+ RootOutput::getInstance()->GetTree()->Branch("Rec35","TBigRIPSReco",&Rec35);
+ RootOutput::getInstance()->GetTree()->Branch("Rec57","TBigRIPSReco",&Rec57);
+ RootOutput::getInstance()->GetTree()->Branch("Rec89","TBigRIPSReco",&Rec89);
+ RootOutput::getInstance()->GetTree()->Branch("Rec911","TBigRIPSReco",&Rec911);
RootOutput::getInstance()->GetTree()->Branch("tof37",&tof37,"tof37/D");
- RootOutput::getInstance()->GetTree()->Branch("aoq35",&aoq35,"aoq35/D");
- RootOutput::getInstance()->GetTree()->Branch("aoq57",&aoq57,"aoq57/D");
- RootOutput::getInstance()->GetTree()->Branch("aoq37",&aoq57,"aoq37/D");
- RootOutput::getInstance()->GetTree()->Branch("tf3",&tf3,"tf3/D");
+ RootOutput::getInstance()->GetTree()->Branch("tof811",&tof811,"tof811/D");
RootOutput::getInstance()->GetTree()->Branch("tf7",&tf7,"tf7/D");
- RootOutput::getInstance()->GetTree()->Branch("beta35",&beta35,"beta35/D");
- RootOutput::getInstance()->GetTree()->Branch("beta57",&beta57,"beta57/D");
- RootOutput::getInstance()->GetTree()->Branch("beta37",&beta57,"beta37/D");
- RootOutput::getInstance()->GetTree()->Branch("z_BR",&z_BR,"z_BR/D");
+ RootOutput::getInstance()->GetTree()->Branch("tf8",&tf8,"tf8/D");
+ RootOutput::getInstance()->GetTree()->Branch("tf11",&tf11,"tf11/D");
RootOutput::getInstance()->GetTree()->Branch("dE_ICF7",&dE_ICF7,"dE_ICF7/D");
+ RootOutput::getInstance()->GetTree()->Branch("dE_ICF11",&dE_ICF11,"dE_ICF11/D");
}
@@ -424,60 +378,29 @@ void Analysis::InitOutputBranch() {
void Analysis::InitInputBranch(){
RootInput:: getInstance()->GetChain()->SetBranchAddress("RunNumber",&RunNumber);
RootInput:: getInstance()->GetChain()->SetBranchAddress("EventNumber",&EventNumber);
-/*
- *
- if(!simulation){
- // Vamos
- RootInput::getInstance()->GetChain()->SetBranchAddress("LTS",<S);
- // Agata
- RootInput::getInstance()->GetChain()->SetBranchAddress("TStrack",&TStrack);
- RootInput::getInstance()->GetChain()->SetBranchAddress("nbTrack",&nbTrack);
- RootInput::getInstance()->GetChain()->SetBranchAddress("trackE",trackE);
- RootInput::getInstance()->GetChain()->SetBranchAddress("trackX1",trackX1);
- RootInput::getInstance()->GetChain()->SetBranchAddress("trackY1",trackY1);
- RootInput::getInstance()->GetChain()->SetBranchAddress("trackZ1",trackZ1);
- RootInput::getInstance()->GetChain()->SetBranchAddress("trackT",trackT);
- RootInput::getInstance()->GetChain()->SetBranchAddress("trackCrystalID",trackCrystalID);
- RootInput::getInstance()->GetChain()->SetBranchAddress("nbCores",&nbCores);
- RootInput::getInstance()->GetChain()->SetBranchAddress("coreId",coreId);
- RootInput::getInstance()->GetChain()->SetBranchAddress("coreTS",coreTS);
- RootInput::getInstance()->GetChain()->SetBranchAddress("coreE0",coreE0);
-
- }
- else{
-
- RootInput:: getInstance()->GetChain()->SetBranchStatus("InitialConditions",true );
- RootInput:: getInstance()->GetChain()->SetBranchStatus("fIC_*",true );
- RootInput:: getInstance()->GetChain()->SetBranchAddress("InitialConditions",&Initial);
- RootInput:: getInstance()->GetChain()->SetBranchStatus("ReactionConditions",true );
- RootInput:: getInstance()->GetChain()->SetBranchStatus("fRC_*",true );
- RootInput:: getInstance()->GetChain()->SetBranchAddress("ReactionConditions",&ReactionConditions);
-
- }
-*/
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::ReInitValue(){
- fX = -9999 ;
- fY = -9999 ;
- fA = -9999 ;
- fB = -9999 ;
- delta35 = -9999 ;
- delta57 = -9999 ;
- Brho35 = -9999 ;
- Brho57 = -9999 ;
+
+ RecFP.clear();
+ std::vector<double> reset;
+ for(int j=0;j<4;j++) reset.push_back(-9999);
+ for(int i=0;i<NmaxFP+1;i++) RecFP.push_back(reset);
+
+ //FP->Clear();
+ //FP->Init(NmaxFP+1);
+ Rec35->Init();
+ Rec57->Init();
+ Rec89->Init();
+ Rec911->Init();
tf3 =-9999;
tf7 =-9999;
tof37 = -9999 ;
- aoq35 = -9999 ;
- aoq57 = -9999 ;
- aoq37 = -9999 ;
-
- beta35 = -9999 ;
- beta57 = -9999 ;
- beta37 = -9999 ;
dE_ICF7 = -9999;
- z_BR = -9999 ;
+ tf8 =-9999;
+ tf11 =-9999;
+ tof811 = -9999 ;
+ dE_ICF11 = -9999;
}
diff --git a/Projects/RIBF196/Analysis.h b/Projects/RIBF196/Analysis.h
index 6fcbf5762ee967e97b206e3b159203af508adcfb..e5d6e4d1f0d65cfa3dabc28b6b841c53704bb18b 100644
--- a/Projects/RIBF196/Analysis.h
+++ b/Projects/RIBF196/Analysis.h
@@ -31,6 +31,8 @@
#include "TBigRIPSPPACPhysics.h"
#include "TBigRIPSPlasticPhysics.h"
#include "TBigRIPSICPhysics.h"
+#include "TBigRIPSFocal.h"
+#include "TBigRIPSReco.h"
#include <TMatrixD.h>
#include <TRandom3.h>
#include <TVector3.h>
@@ -77,11 +79,12 @@ class Analysis: public NPL::VAnalysis{
void End();
void ReadXmls();
//TMatrixD* RecReadTransferMatrix(char*);
- TMatrixD RecReadTransferMatrix(string);
+ //TMatrixD RecReadTransferMatrix(string);
+ std::vector<std::vector<double>> RecReadTransferMatrix2(string);
std::vector<double> RecFPposition(std::vector<TVector2>,std::vector<int>);
- double RecDeltaBrho(std::vector<double>,std::vector<double>, TMatrixD matrix);
- double RecAoqOneFold(double, double, double);
- double RecZ(double,double,double);
+ //double RecDeltaBrho(std::vector<double>,std::vector<double>, TMatrixD matrix);
+ //double RecAoqOneFold(double, double, double);
+ double RecZ(double,double,double, bool);
void InitOutputBranch();
void InitInputBranch();
@@ -90,28 +93,16 @@ class Analysis: public NPL::VAnalysis{
private:
//calculated variables
- double aoq;
- double zet;
- double beta;
- double delta35;
- double delta57;
- double Brho35;
- double Brho57;
double tof37;
- double aoq35;
- double aoq57;
- double aoq37;
- double beta35;
- double beta57;
- double beta37;
+ double tof811;
double dE_ICF7;
+ double dE_ICF11;
double z_BR;
+ double z_ZD;
double tf3;
double tf7;
- double fX;
- double fY;
- double fA;
- double fB;
+ double tf8;
+ double tf11;
// List of Dipole Brho, manually set for the moment
// But method to read from root file should be implemented
@@ -124,17 +115,24 @@ class Analysis: public NPL::VAnalysis{
double BrhoD7 = 5.9429;
double BrhoD8 = 5.9381;
- //double tof_offset37 = 0.;
- //double tof_offset37 = 300.06;
- double tof_offset37 = 299.968;
+ double length35;
+ double length57;
+ double length37;
+ double length89;
+ double length911;
+ double length811;
+
+ double tof_offset37 = 300.06;
+ double tof_offset811 = -134.503;
const double clight = 299.7792458; // in mm/ns
const double mnucleon = 931.49432; // MeV
-
- // intermediate variables
- TMatrixD matF35;
- TMatrixD matF57;
+ //5 columns/ 6 rows transfer matrices
+ vector<vector<double>> matrixF35;
+ vector<vector<double>> matrixF57;
+ vector<vector<double>> matrixF89;
+ vector<vector<double>> matrixF911;
// Branches and detectors
@@ -144,6 +142,15 @@ class Analysis: public NPL::VAnalysis{
map<int,BigRIPSPPACHitList> FP_PPACHitList ;
int EventNumber;
int RunNumber;
+ int Trigger;
+ unsigned long long int TimeStamp;
+ //Focal* FP;
+ //TBigRIPSFocal* FP;
+ TBigRIPSReco* Rec35;
+ TBigRIPSReco* Rec57;
+ TBigRIPSReco* Rec89;
+ TBigRIPSReco* Rec911;
+ std::vector<std::vector<double>> RecFP;
//maps containings infos from input xml config files
// for PPACs
@@ -159,11 +166,14 @@ class Analysis: public NPL::VAnalysis{
map<int,int> IC_IDtoFP;
map<int, string> IC_IDtoName;
map<string,int> IC_NametoID;
+ map<string,vector<double>> IC_Zcoef;
+ map<string,double> IC_Ionpair;
// for FocalPlanes
map<int,int> FP_IDtoFP;
map<int,int> FPtoID;
map<int,double> FP_Z;
map<int,double> FP_Zoffset;
+ int NmaxFP;
};