Skip to content

GitLab

Commit 582f6e2d authored by flavigny's avatar flavigny
Browse files

* Work on AoQ reconstrction in BigRIPS

parent 5a8739f2
Pipeline #123275 passed with stages
in 8 minutes and 17 seconds
Hide whitespace changes
Inline Side-by-side
Showing with 437 additions and 33 deletions
NPLib/Detectors/BigRIPS/TBigRIPSReco.cxx
View file @ 582f6e2d
......@@ -12,6 +12,8 @@ void TBigRIPSReco::Init(){
aoq=-9999;
beta=-9999;
delta=-9999;
brho1=-9999;
brho2=-9999;
brho=-9999;
}
......@@ -55,12 +57,58 @@ void TBigRIPSReco::RecBrho(std::vector<double> RecFPUpstream,std::vector<double>
brho = BrhoCentral*(1.0+delta*0.01);
}
////////////////////////////////////////////////////////////////////////////////
void TBigRIPSReco::RecAoqOne(double tof, double length){
beta = length /(tof * c_mm_ns);
void TBigRIPSReco::RecAoqOneFold(double tof, double length){
beta = length /(tof * c_light);
double gamma = 1./sqrt(1 - beta*beta);
aoq = (brho * c_mm_ns) / (mnucleon * beta * gamma);
//std::cout << aoq << std::endl;
aoq = (brho * c_light) / (amu_c2 * beta * gamma);
}
////////////////////////////////////////////////////////////////////////////////
void TBigRIPSReco::RecAoqTwoFold(double tof, double length1, double length2, int useBeta){
if(!(length1>0 && length2>0)){
std::cout << "Length1 or Length2 is not positive (or both)" << std::endl;
return;
}
if(!(brho1>0 && brho2>0)){
//std::cout << "brho1 or brho2 was not properly initialized/setup" << std::endl;
//std::cout << "brho1:" <<brho1<< std::endl;
//std::cout << "brho2:" <<brho2<< std::endl;
return;
}
double alpha = brho2 / brho1;
double a1 = sqrt(alpha * alpha * c_squared * tof * tof
+ (pow(alpha,4) - alpha * alpha) * length1 * length1
+ (1 - alpha*alpha) * length2 * length2);
double rbeta1 = ( a1 * length1 + length2 * c_light * tof ) /
( a1 * c_light * tof + (1 - alpha * alpha) * length1 * length2);
double gamma1 = 1 / sqrt(1 - pow(rbeta1,2));
double rbeta2 = ( a1 * length1 + length2 * c_light * tof ) /
( c_squared * tof * tof + (alpha * alpha - 1) * length1 * length1);
double gammab = 1 / sqrt(1 - pow(rbeta2,2));
aoq = brho1 * c_light / amu_c2 / rbeta1 / gamma1; // should be same as brho2/beta2/gamma2
if(useBeta == 2 ) beta = rbeta2;
else if(useBeta == 1) beta = rbeta1;
brho = brho2;
/*
std::cout << "-------------" << std::endl;
std::cout << "brho1:" <<brho1<< std::endl;
std::cout << "brho2:" <<brho2<< std::endl;
std::cout << "brho:" <<brho<< std::endl;
std::cout << "alpha:" <<alpha<< std::endl;
std::cout << "rbeta1:" <<rbeta1<< std::endl;
std::cout << "rbeta2:" <<rbeta2<< std::endl;
std::cout << "aoq:" <<aoq<< std::endl;
*/
}
////////////////////////////////////////////////////////////////////////////////
......
NPLib/Detectors/BigRIPS/TBigRIPSReco.h
View file @ 582f6e2d
......@@ -3,6 +3,10 @@
#include "TObject.h"
#include <vector>
#include <iostream>
#include "NPPhysicalConstants.h"
using namespace NPUNITS;
class TBigRIPSReco: public TObject{
public:
......@@ -13,6 +17,8 @@ class TBigRIPSReco: public TObject{
double aoq;
double beta;
double delta;
double brho1; // upstream brho if twofold
double brho2; // downstream brho if twofold
double brho;
double angle;
double z;
......@@ -21,17 +27,29 @@ class TBigRIPSReco: public TObject{
void Init();
void Print();
void RecBrho(std::vector<double>,std::vector<double>,std::vector<std::vector<double>>,double);
void RecAoqOne(double, double);
void RecAoqOneFold(double, double);
void RecAoqTwoFold(double, double, double, int);
void RecZet(double,double, std::vector<double>);
const double c_mm_ns = 299.7792458; // in mm/ns
const double mnucleon = 931.49432; // MeV
//const double c_mm_ns = NPUNITS::c_light; // 299.779 mm/ns
//const double amu_c2 = NPUNITS::amu_c2; // 931.494 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 SetBrho1(double value){brho1=value;}
void SetBrho2(double value){brho2=value;}
void SetZ(double value){z=value;}
double GetAoq(){return aoq;}
double GetBeta(){return beta;}
double GetDelta(){return delta;}
double GetBrho1(){return brho1;}
double GetBrho2(){return brho2;}
double GetBrho(){return brho;}
double GetAngle(){return angle;}
double GetZ(){return z;}
ClassDef(TBigRIPSReco,1);
};
......
Projects/RIBF196/Analysis.cxx
View file @ 582f6e2d
......@@ -40,13 +40,14 @@ Analysis::~Analysis(){
////////////////////////////////////////////////////////////////////////////////
void Analysis::Init() {
NmaxFP=11;
//FP = new TBigRIPSFocal();
//FP->Init(NmaxFP+1);
NmaxFP=12;
Rec35 = new TBigRIPSReco();
Rec57 = new TBigRIPSReco();
Rec37 = new TBigRIPSReco();
Rec89 = new TBigRIPSReco();
Rec911 = new TBigRIPSReco();
Rec811 = new TBigRIPSReco();
InitOutputBranch();
InitInputBranch();
......@@ -56,12 +57,7 @@ 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");
......@@ -165,15 +161,6 @@ void Analysis::TreatEvent() {
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){
if(it->first < NmaxFP){
FP->SetFPTrack(it->first, RecFPposition(it->second.PPACHit,it->second.PPACID) );
//FP->Print(it->first);
}
}
*/
// Brho/delta reconstruction
Rec35->RecBrho(RecFP[3],RecFP[5],matrixF35,BrhoD3);
......@@ -181,7 +168,6 @@ void Analysis::TreatEvent() {
Rec89->RecBrho(RecFP[8],RecFP[9],matrixF89,BrhoD7);
Rec911->RecBrho(RecFP[9],RecFP[11],matrixF911,BrhoD8);
/////////////////////////////////
/// STEP2: TOF, flight length ///
/////////////////////////////////
......@@ -204,11 +190,19 @@ void Analysis::TreatEvent() {
/////////////////////////////////
// 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);
Rec35->RecAoqOneFold(tof37,length37);
Rec57->RecAoqOneFold(tof37,length37);
Rec37->SetBrho1(Rec35->GetBrho());
Rec37->SetBrho2(Rec57->GetBrho());
Rec37->RecAoqTwoFold(tof37,length35,length57,2);
//
Rec89->RecAoqOneFold(tof811,length811);
Rec911->RecAoqOneFold(tof811,length811);
Rec811->SetBrho1(Rec89->GetBrho());
Rec811->SetBrho2(Rec911->GetBrho());
Rec811->RecAoqTwoFold(tof811,length89,length911,2);
///////////////////////////////////////////////
/// STEP4: Z reco. from dE in IC and Beta ///
///////////////////////////////////////////////
......@@ -221,11 +215,37 @@ void Analysis::TreatEvent() {
dE_ICF7 = IC_dE["F7IC"];
Rec35->RecZet(dE_ICF7,IC_Ionpair["F7IC"],IC_Zcoef["F7IC"]);
Rec57->RecZet(dE_ICF7,IC_Ionpair["F7IC"],IC_Zcoef["F7IC"]);
Rec37->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"]);
Rec811->RecZet(dE_ICF11,IC_Ionpair["F11IC"],IC_Zcoef["F11IC"]);
///////////////////////////////////////////////
/// STEP5: Higher order AoQ optical corr. ///
///////////////////////////////////////////////
double x[8];
if(Rec37->GetAoq()>0){
x[0] = RecFP[5][0]; x[1] = RecFP[5][2]; x[2] = RecFP[5][1]; x[3] = RecFP[5][3];
x[4] = RecFP[7][0]; x[5] = RecFP[7][2]; x[6] = RecFP[7][1]; x[7] = RecFP[7][3];
aoqc37 = Rec37->GetAoq() - BR_OptCorr_Ga(x);
}
if(Rec811->GetAoq()>0){
x[0] = RecFP[9][0]; x[1] = RecFP[9][2]; x[2] = RecFP[9][1]; x[3] = RecFP[9][3];
x[4] = RecFP[11][0]; x[5] = RecFP[11][2]; x[6] = RecFP[11][1]; x[7] = RecFP[11][3];
aoqc811 = Rec811->GetAoq() - ZD_OptCorr_Ga(x);
}
/*
for(int i=0; i<8; i++){
std::cout << "x["<<i<<"]:"<<x[i]<< std::endl;
}
std::cout << "aoq:"<<Rec811->GetAoq() << std::endl;
std::cout << "corr:"<<ZD_OptCorr_Ga(x) << std::endl;
std::cout << "aoqc:"<<aoqc811 << std::endl;
std::cout << "________________"<< std::endl;
*/
}
////////////////////////////////////////////////////////////////////////////////
......@@ -348,6 +368,75 @@ std::vector<double> Analysis::RecFPposition(std::vector<TVector2> HitList,std::v
}
////////////////////////////////////////////////////////////////////////////////
double Analysis::ZD_OptCorr_Ga(double *x) {
//double returnValue = gGaZDDMean;
double returnValue = 0;
int i = 0, j = 0, k = 0;
for (i = 0; i < gGaZDNCoefficients ; i++) {
// Evaluate the ith term in the expansion
double term = gGaZDCoefficient[i];
for (j = 0; j < gGaZDNVariables; j++) {
// Evaluate the polynomial in the jth variable.
int power = gGaZDPower[gGaZDNVariables * i + j];
double p1 = 1, p2 = 0, p3 = 0, r = 0;
double v = 1 + 2. / (gGaZDXMax[j] - gGaZDXMin[j]) * (x[j] - gGaZDXMax[j]);
// what is the power to use!
switch(power) {
case 1: r = 1; break;
case 2: r = v; break;
default:
p2 = v;
for (k = 3; k <= power; k++) {
p3 = p2 * v;
p1 = p2; p2 = p3;
}
r = p3;
}
// multiply this term by the poly in the jth var
term *= r;
}
// Add this term to the final result
returnValue += term;
}
return returnValue;
}
////////////////////////////////////////////////////////////////////////////////
double Analysis::BR_OptCorr_Ga(double *x) {
//double returnValue = gGaBRDMean;
double returnValue = 0;
int i = 0, j = 0, k = 0;
for (i = 0; i < gGaBRNCoefficients ; i++) {
// Evaluate the ith term in the expansion
double term = gGaBRCoefficient[i];
for (j = 0; j < gGaBRNVariables; j++) {
// Evaluate the polynomial in the jth variable.
int power = gGaBRPower[gGaBRNVariables * i + j];
double p1 = 1, p2 = 0, p3 = 0, r = 0;
double v = 1 + 2. / (gGaBRXMax[j] - gGaBRXMin[j]) * (x[j] - gGaBRXMax[j]);
// what is the power to use!
switch(power) {
case 1: r = 1; break;
case 2: r = v; break;
default:
p2 = v;
for (k = 3; k <= power; k++) {
p3 = p2 * v;
p1 = p2; p2 = p3;
}
r = p3;
}
// multiply this term by the poly in the jth var
term *= r;
}
// Add this term to the final result
returnValue += term;
}
return returnValue;
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::End(){
}
......@@ -359,11 +448,12 @@ void Analysis::InitOutputBranch() {
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("Rec37","TBigRIPSReco",&Rec37);
RootOutput::getInstance()->GetTree()->Branch("Rec89","TBigRIPSReco",&Rec89);
RootOutput::getInstance()->GetTree()->Branch("Rec911","TBigRIPSReco",&Rec911);
RootOutput::getInstance()->GetTree()->Branch("Rec811","TBigRIPSReco",&Rec811);
RootOutput::getInstance()->GetTree()->Branch("tof37",&tof37,"tof37/D");
RootOutput::getInstance()->GetTree()->Branch("tof811",&tof811,"tof811/D");
RootOutput::getInstance()->GetTree()->Branch("tf7",&tf7,"tf7/D");
......@@ -371,6 +461,8 @@ void Analysis::InitOutputBranch() {
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");
RootOutput::getInstance()->GetTree()->Branch("aoqc37",&aoqc37,"aoqc37/D");
RootOutput::getInstance()->GetTree()->Branch("aoqc811",&aoqc811,"aoqc811/D");
}
......@@ -391,8 +483,12 @@ void Analysis::ReInitValue(){
//FP->Init(NmaxFP+1);
Rec35->Init();
Rec57->Init();
Rec37->Init();
Rec89->Init();
Rec911->Init();
Rec811->Init();
aoqc37 =-9999;
aoqc811 =-9999;
tf3 =-9999;
tf7 =-9999;
tof37 = -9999 ;
......
Projects/RIBF196/Analysis.h
View file @ 582f6e2d
......@@ -14,7 +14,7 @@
* Last update : *
*---------------------------------------------------------------------------*
* Decription: *
* Basic AoQ,Z reconstruction for RIBF196 experiment *
* AoQ,Z reconstruction for RIBF196 experiment *
* *
*---------------------------------------------------------------------------*
* Comment: *
......@@ -93,6 +93,8 @@ class Analysis: public NPL::VAnalysis{
private:
//calculated variables
double aoqc37;
double aoqc811;
double tof37;
double tof811;
double dE_ICF7;
......@@ -148,8 +150,10 @@ class Analysis: public NPL::VAnalysis{
//TBigRIPSFocal* FP;
TBigRIPSReco* Rec35;
TBigRIPSReco* Rec57;
TBigRIPSReco* Rec37;
TBigRIPSReco* Rec89;
TBigRIPSReco* Rec911;
TBigRIPSReco* Rec811;
std::vector<std::vector<double>> RecFP;
//maps containings infos from input xml config files
......@@ -175,6 +179,173 @@ class Analysis: public NPL::VAnalysis{
map<int,double> FP_Zoffset;
int NmaxFP;
// Higher order Optical corrections
double ZD_OptCorr_Ga(double *);
constexpr static int gGaZDNVariables = 8;
constexpr static int gGaZDNCoefficients = 10;
constexpr static double gGaZDDMean = 2.67982;
// Assignment to mean vector.
constexpr static double gGaZDXMean[] = {
-0.151618, -0.97249, -0.226697, -0.0028457, -0.0250826, 0.00288959, -0.0263575, 0.180498 };
// Assignment to minimum vector.
constexpr static double gGaZDXMin[] = {
-112.762, -67.3691, -21.6583, -1461.26, -19.009, -25.9623, -63.8841, -82.8509 };
// Assignment to maximum vector.
constexpr static double gGaZDXMax[] = {
110.757, 6247.9, 37.3008, 1461.03, 19.6172, 25.5737, 55.9905, 60.4835 };
// Assignment to coefficients vector.
constexpr static double gGaZDCoefficient[] = {
0.0045763,
0.0260008,
0.0202036,
-0.00373027,
0.0117822,
-0.00717916,
-0.010714,
-6.11157e-06,
0.0045109,
-0.00397418
};
// Assignment to error coefficients vector.
constexpr static double gGaZDCoefficientRMS[] = {
0.0263103,
0.0627114,
0.0645658,
0.0289265,
0.0959883,
0.0450067,
0.125588,
0.0482148,
0.0962523,
0.0798097
};
// Assignment to powers vector.
// The powers are stored row-wise, that is
// p_ij = gPower[i * NVariables + j];
constexpr static int gGaZDPower[] = {
1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 2, 1,
1, 1, 1, 1, 2, 1, 1, 1,
1, 1, 1, 1, 1, 2, 1, 1,
1, 1, 2, 1, 1, 1, 1, 1,
3, 1, 1, 1, 1, 1, 1, 1,
2, 1, 1, 1, 1, 1, 3, 1,
1, 1, 1, 1, 1, 1, 1, 2,
1, 1, 1, 1, 1, 1, 1, 3,
1, 1, 1, 1, 1, 3, 1, 1
};
double BR_OptCorr_Ga(double *);
constexpr static int gGaBRNVariables = 8;
constexpr static int gGaBRNCoefficients = 25;
constexpr static double gGaBRDMean = 2.67489;
// Assignment to mean vector.
constexpr static double gGaBRXMean[] = {
-0.0642977, -0.00307703, -0.00100196, -0.0107261, -0.0131277, -0.00317987, -0.158857, -0.0435761 };
// Assignment to minimum vector.
constexpr static double gGaBRXMin[] = {
-107.29, -25.9837, -31.9689, -42.8947, -29.3661, -31.5436, -20.0646, -37.375 };
// Assignment to maximum vector.
constexpr static double gGaBRXMax[] = {
64.5168, 21.8811, 30.0843, 43.7818, 35.0687, 23.5787, 29.204, 35.7001 };
// Assignment to coefficients vector.
constexpr static double gGaBRCoefficient[] = {
0.000721757,
0.00503421,
-0.00479955,
0.00219026,
0.00552134,
0.00247818,
0.0026914,
-0.000758646,
0.00743853,
0.0162762,
0.0157268,
0.00235257,
-0.0029327,
-0.0115376,
-0.0139695,
-0.0526337,
0.000312724,
-0.00270404,
0.00261865,
-0.00281391,
0.00284842,
-0.00417684,
0.029841,
-0.0494852,
0.00822339
};
// Assignment to error coefficients vector.
constexpr static double gGaBRCoefficientRMS[] = {
0.0221415,
0.0946943,
0.195954,
0.0827622,
0.10566,
0.0550079,
0.091793,
0.0767151,
0.172153,
0.94549,
0.578329,
0.49222,
0.16515,
0.50208,
0.489808,
1.89549,
0.0284801,
0.100958,
0.121855,
0.161108,
0.277453,
0.292148,
4.02157,
3.97904,
0.647432
};
// Assignment to powers vector.
// The powers are stored row-wise, that is
// p_ij = gPower[i * NVariables + j];
constexpr static int gGaBRPower[] = {
1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 2, 1,
1, 1, 1, 1, 1, 2, 1, 2,
1, 1, 1, 1, 1, 1, 3, 1,
2, 1, 1, 1, 1, 1, 2, 1,
2, 1, 1, 1, 1, 1, 1, 1,
1, 1, 2, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 2, 2, 1,
1, 1, 2, 1, 2, 1, 1, 1,
1, 1, 1, 1, 2, 1, 1, 3,
2, 1, 1, 1, 2, 1, 2, 1,
1, 2, 1, 2, 1, 1, 2, 1,
4, 1, 1, 1, 1, 1, 1, 1,
2, 1, 1, 1, 3, 1, 1, 1,
1, 1, 1, 1, 3, 1, 3, 1,
1, 2, 4, 1, 3, 2, 1, 1,
1, 1, 1, 1, 1, 2, 1, 1,
3, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 3, 1, 1,
2, 1, 1, 1, 2, 1, 1, 1,
1, 2, 1, 1, 1, 2, 2, 1,
3, 1, 1, 1, 2, 1, 1, 1,
1, 1, 1, 1, 2, 1, 2, 3,
1, 1, 2, 1, 2, 1, 1, 3,
3, 1, 1, 1, 1, 1, 3, 1
};
};
......
Projects/RIBF196/CheckFP.C 0 → 100644
View file @ 582f6e2d
void CheckFP(int id=3){
auto tree = new TChain("PhysicsTree");
//tree->Add("../../Outputs/Analysis/PhysicsTree.root");
tree->Add("../../Outputs/Analysis/ribf196_0030_phy.root");
TString Cond="";
auto c = new TCanvas();
c->Divide(3,2);
string buffer, command;
buffer = Form("RecFP[%i]",id);
c->cd(1);
command = buffer + "[0]>>hfX(300,-150,150)";
cout<< command<<endl;
tree->Draw(command.c_str(),Cond,"colz");
c->cd(2);
command = buffer + "[2]>>hfY(300,-150,150)";
cout<< command<<endl;
tree->Draw(command.c_str(),Cond,"colz");
c->cd(3);
command = buffer + "[2]:"+ buffer +"[0]>>hfXY(300,-150,150,300,-150,150)";
cout<< command<<endl;
tree->Draw(command.c_str(),Cond,"colz");
c->cd(4);
command = buffer + "[1]>>hfA(200,-50,50)";
cout<< command<<endl;
tree->Draw(command.c_str(),Cond,"colz");
c->cd(5);
command = buffer + "[3]>>hfB(200,-50,50)";
cout<< command<<endl;
tree->Draw(command.c_str(),Cond,"colz");
c->cd(5);
command = buffer + "[3]:"+ buffer +"[1]>>hfAB(200,-50,50,200,-50,50)";
c->cd(6);
}
Projects/RIBF196/CheckReco.C 0 → 100644
View file @ 582f6e2d
void CheckReco(){
auto tree = new TChain("PhysicsTree");
//tree->Add("../../Outputs/Analysis/PhysicsTree.root");
tree->Add("../../Outputs/Analysis/ribf196_0030_phy.root");
//TString Cond="";
auto c = new TCanvas();
c->Divide(3,2);
c->cd(1);