Commit ce6a6ada authored by flavigny's avatar flavigny
Browse files

* Added AoQ, Brho, Z reconstruction to BigRIPS

parent d32dd260
......@@ -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})
#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);
#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
......@@ -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);
#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);
#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
......@@ -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",&beta3