From 93296f505d6e670bfbb7cdecd72c0f68def9b2ef Mon Sep 17 00:00:00 2001 From: Miguel Lozano-Gonzalez <miguellozano.gonzalez@usc.es> Date: Tue, 21 Mar 2023 11:16:39 +0100 Subject: [PATCH] added basic postanalysis macros --- Projects/e748/Postanalysis/AnalyzeBeam.cpp | 44 +++++++ Projects/e748/Postanalysis/AnalyzeCATS.cpp | 34 +++++ Projects/e748/Postanalysis/ApplyCuts.cpp | 52 ++++++++ Projects/e748/Postanalysis/CalibrateBeam.cpp | 18 +++ Projects/e748/Postanalysis/GetCuts.cpp | 117 ++++++++++++++++++ Projects/e748/Postanalysis/GetSpectrum.cpp | 62 ++++++++++ Projects/e748/Postanalysis/ImprovePID.cpp | 85 +++++++++++++ .../PlotTheoreticalKinematics.cpp | 47 +++++++ Projects/e748/Postanalysis/Utils.cpp | 35 ++++++ 9 files changed, 494 insertions(+) create mode 100644 Projects/e748/Postanalysis/AnalyzeBeam.cpp create mode 100644 Projects/e748/Postanalysis/AnalyzeCATS.cpp create mode 100644 Projects/e748/Postanalysis/ApplyCuts.cpp create mode 100644 Projects/e748/Postanalysis/CalibrateBeam.cpp create mode 100644 Projects/e748/Postanalysis/GetCuts.cpp create mode 100644 Projects/e748/Postanalysis/GetSpectrum.cpp create mode 100644 Projects/e748/Postanalysis/ImprovePID.cpp create mode 100644 Projects/e748/Postanalysis/PlotTheoreticalKinematics.cpp create mode 100644 Projects/e748/Postanalysis/Utils.cpp diff --git a/Projects/e748/Postanalysis/AnalyzeBeam.cpp b/Projects/e748/Postanalysis/AnalyzeBeam.cpp new file mode 100644 index 000000000..a6d39929d --- /dev/null +++ b/Projects/e748/Postanalysis/AnalyzeBeam.cpp @@ -0,0 +1,44 @@ +#include <ROOT/RDataFrame.hxx> +#include <TCanvas.h> +#include <TChain.h> +#include <TFile.h> +#include <TH2D.h> +#include <TH1D.h> +#include <TString.h> +#include <TTree.h> +#include <TROOT.h> +#include <TStyle.h> + +#include "TMust2Physics.h" +#include "Utils.cpp" + +#include <set> + +void AnalyzeBeam() +{ + ROOT::EnableImplicitMT(); + //gStyle->SetOptStat("nmeruoi"); + + //ROOT::RDataFrame df("BaseCutTree", "./RootFiles/BaseCutTree_12Be.root"); + auto df {ReadAll12BeData()}; + + auto gated {df.Filter( + [](const TMust2Physics& must2) + { + return (must2.TelescopeNumber.size() == 0); + }, + {"MUST2"})}; + auto hCATSchio {gated.Histo2D({"hCATSchio", "Raw TOF vs IC", 700, 0., 300., 900, 0, 900}, "T_CATS1_CAV", "IC_E")}; + auto hCATSplast {gated.Histo2D({"hCATSplast", "Raw TOF vs QPlast", 700, 0., 300., 4000, 0, 16000}, "T_CATS1_CAV", "QPlast")}; + auto hCHIOplast {gated.Histo2D({"hCHIOplast", "CHIO vs QPlast", 900, 0, 900, 4000, 0, 16000}, "IC_E", "QPlast")}; + + //plotting + auto* c1 {new TCanvas("c1", "Beam ID")}; + c1->DivideSquare(4); + c1->cd(1); + hCATSchio->DrawClone("colz"); + c1->cd(2); + hCATSplast->DrawClone("colz"); + c1->cd(3); + hCHIOplast->DrawClone("colz"); +} diff --git a/Projects/e748/Postanalysis/AnalyzeCATS.cpp b/Projects/e748/Postanalysis/AnalyzeCATS.cpp new file mode 100644 index 000000000..16660ad1c --- /dev/null +++ b/Projects/e748/Postanalysis/AnalyzeCATS.cpp @@ -0,0 +1,34 @@ +#include "TROOT.h" +#include "TFile.h" +#include "TTree.h" +#include "TCanvas.h" +#include "TString.h" +#include "TCATSPhysics.h" +#include "TMust2Physics.h" +#include "TModularLeafPhysics.h" + +#include "ROOT/RDataFrame.hxx" + +void AnalyzeCATS() +{ + ROOT::EnableImplicitMT(); + + auto* inFile {new TFile("/home/miguel/nptool/Projects/e748/output/analysis/Test_12Be.root")}; + //inFile->ls(); + auto* tree {inFile->Get<TTree>("PhysicsTree")}; + //tree->Print(); + + ROOT::RDataFrame df(*tree); + auto cond {df.Filter("CATS.PositionOnTargetX != -1000 && CATS.PositionOnTargetY != -1000")}; + auto hPos {cond.Define("x", "CATS.PositionOnTargetX").Define("y", "CATS.PositionOnTargetY") + .Histo2D({"hPos", "Position on CATS;X;Y", 100, -70., 70., 100, -70., 70.}, "x", "y")}; + hPos->DrawClone("colz"); + + auto hCATS1Cal {cond.Histo1D("T_CATS1_CAV_Cal")}; + hCATS1Cal->DrawClone(); + + auto hCATSBeamEnergy {cond.Filter("T_CATS1_CAV > 0") + .Histo2D({"hCATSBeamEnergy", "CATS against BeamEnergy", 100, 0., 300., 200, 0., 500.}, "T_CATS1_CAV", "BeamEnergy")}; + hCATSBeamEnergy->DrawClone("colz"); + +} diff --git a/Projects/e748/Postanalysis/ApplyCuts.cpp b/Projects/e748/Postanalysis/ApplyCuts.cpp new file mode 100644 index 000000000..d0716e879 --- /dev/null +++ b/Projects/e748/Postanalysis/ApplyCuts.cpp @@ -0,0 +1,52 @@ +#include "TROOT.h" +#include "TStyle.h" +#include "TFile.h" +#include "TTree.h" +#include "TChain.h" +#include "TCATSPhysics.h" +#include "TMust2Physics.h" + +#include "ROOT/RDataFrame.hxx" +#include <ROOT/RVec.hxx> + +#include "Utils.cpp" + +#include <Rtypes.h> +#include <TCanvas.h> +#include <algorithm> +#include <iostream> +#include <vector> + +void ApplyCuts() +{ + ROOT::EnableImplicitMT(); + gStyle->SetOptStat("nmeruoi"); + + ROOT::RDataFrame df("BaseCutTree", "./RootFiles/BaseCutTree_12Be.root"); + auto dsc {df.Describe()}; + //dsc.Print(); + + //--> Read TCutGs + auto* cutChio {GetGraphCutFromFile("./Cuts/chio_v0.root")}; + auto* cutPID {GetGraphCutFromFile("./../TOF.root", "TOF")}; + //--> Apply them. Must apply cuts to the vectors stored in classes + auto gated {df.Filter([&](const ROOT::RVecD& vSiE, const ROOT::RVecD& pid, + const unsigned short& chioE, const double& Qplast) + { + std::vector<bool> vCondPID {}; + for(int i = 0; i < vSiE.size(); i++) + { + vCondPID.push_back(cutPID->IsInside(vSiE.at(i), pid.at(i))); + } + //all hits in MUST2 should be inside TCutG + bool condPID {std::all_of(vCondPID.begin(), vCondPID.end(), [](const bool& e){return e;})}; + //same but for heavy ID + bool condChio {static_cast<bool>(cutChio->IsInside(Qplast, chioE))}; + return (condChio && condPID); + }, + {"Must2SiE", "TOFLight", "IC_E", "QPlast"})}; + + + //--> Write to file + gated.Snapshot("FinalTree", "./RootFiles/FinalTree_12Be.root"); +} diff --git a/Projects/e748/Postanalysis/CalibrateBeam.cpp b/Projects/e748/Postanalysis/CalibrateBeam.cpp new file mode 100644 index 000000000..bfe9e65eb --- /dev/null +++ b/Projects/e748/Postanalysis/CalibrateBeam.cpp @@ -0,0 +1,18 @@ +#include <ROOT/RDataFrame.hxx> +#include <TROOT.h> + +#include "TMust2Physics.h" +#include "Utils.cpp" + +void CalibrateBeam() +{ + ROOT::EnableImplicitMT(); + //we can use all the data, because its measured between Caviar and CATS + auto df {ReadAll12BeData()}; + + //plot original values without telescope hits + auto hOri {df.Filter([](const TMust2Physics& must2){return must2.TelescopeNumber.size() == 0;}, {"MUST2"}) + .Histo1D({"hOri", "Raw TOF", 1000, 100., 350.}, "T_CATS1_CAV")}; + + hOri->DrawClone(); +} diff --git a/Projects/e748/Postanalysis/GetCuts.cpp b/Projects/e748/Postanalysis/GetCuts.cpp new file mode 100644 index 000000000..e330c4e8f --- /dev/null +++ b/Projects/e748/Postanalysis/GetCuts.cpp @@ -0,0 +1,117 @@ +#include "TROOT.h" +#include "TStyle.h" +#include "TFile.h" +#include "TTree.h" +#include "TChain.h" +#include "TCATSPhysics.h" +#include "TMust2Physics.h" +#include "Math/Point3D.h" + +#include "ROOT/RDataFrame.hxx" +#include <ROOT/RVec.hxx> + +#include "Utils.cpp" + +#include <Rtypes.h> +#include <TCanvas.h> +#include <algorithm> +#include <iostream> +#include <vector> + +void GetCuts() +{ + ROOT::EnableImplicitMT(); + ///////////////////////////////////////////////////// + ///READ DATA + gStyle->SetOptStat("nmeruoi"); + + auto* inFile {new TFile("/home/miguel/nptool/Projects/e748/output/analysis/Physics_12Be.root")}; + //inFile->ls(); + auto* tree {inFile->Get<TTree>("PhysicsTree")}; + //tree->Print(); + + //load also chio chain + auto* chioChain {new TChain("Numexo2")};//runs for 12Be + std::set<int> runs {315, 316, 317, 318, 320, 321, 323, 325, 326, 327, 328, 329, 330, 331, 339, 341, 342, 346, 347, 348}; + for(const auto& run : runs) + chioChain->Add(TString::Format("/home/miguel/nptool/Projects/e748/Data/Merged/run_%04d.root", run)); + tree->AddFriend(chioChain); + + ////// DATAFRAME + ROOT::RDataFrame df(*tree); + //Apply base cut + auto baseCut = [](const TCATSPhysics& cats, + const std::vector<short>& vtel, const std::vector<double>& vCsIE, + unsigned short chioE, double QPlast, double T_CATS1_CAV) + { + //0--> Only one hit per event! + bool condSize {vtel.size() == 1}; + //1--> Position on target should be inside its dimensions + bool condPosition {cats.PositionOnTargetY > -10 && cats.PositionOnTargetY < 10. + && cats.PositionOnTargetX > -20. && cats.PositionOnTargetX < 20.}; + //2--> Telescope number lower than 5 + bool condTelescopeN {std::all_of(vtel.begin(), vtel.end(), + [](const auto& t){return t < 5;})}; + //3--> CsI energy == 0!! + bool condCsIEnergy {std::all_of(vCsIE.begin(), vCsIE.end(), + [](const auto& e){return e < 0;})}; + //4--> Recoil measured in IC and plastic + bool condRecoil {chioE > 0 && QPlast > 0}; + //5--> TOF for beam, raw cut + bool condBeam {T_CATS1_CAV > 160. && T_CATS1_CAV < 235.}; + return (condSize && condPosition && condTelescopeN && condCsIEnergy && condRecoil && condBeam); + }; + + auto gated {df.Filter(baseCut, + {"CATS", + "Must2Telescopes", "Must2CsIE", + "IC_E", "QPlast", "T_CATS1_CAV"})}; + + //Build ID for light recoil, based on TOF between CATS and MUST2 + gated = gated.Define("TOFLight", + [](const std::vector<short>& vTel, const std::vector<double>& vSiT, + const double& timeCorr) + { + ROOT::RVecD v {};//serialize use of RVec even though we are working with a size=1 vector + for(int i = 0; i < vTel.size(); i++) + { + auto val1 {vSiT.at(i) + timeCorr}; + auto tn {vTel.at(i)}; + double opt {}; + if(tn == 1) + opt = -2.521; + else if(tn == 2) + opt = 0.148; + else if(tn == 3) + opt = -1.922; + else if(tn == 4) + opt = -7.176; + v.push_back(val1 + opt); + } + return v; + }, + {"Must2Telescopes", "Must2SiT", "TimeCorr"}); + + // //HISTOGRAMS + auto hCHIO { gated.Histo2D({"hCHIO", "Heavy ID;QPlast;IC_E", 4000, 0., 16000., 900, 0., 900.}, "QPlast", "IC_E")}; + auto hPID {gated.Histo2D({"hPID", "PID;E_{Si} [MeV];TOF_{Heavy}[au]", 1000, 0., 30., 1000, 460., 580.}, "Must2SiE", "TOFLight")}; + + //--> Read TCutGs + auto* cutChio {GetGraphCutFromFile("./Cuts/chio_v0.root")}; + auto* cutPID {GetGraphCutFromFile("./../TOF.root", "TOF")}; + + //plotting + auto* cID {new TCanvas("cID", "ID canvas", 1)}; + cID->DivideSquare(2); + cID->cd(1); + hCHIO->DrawClone("colz"); + cutChio->SetLineColor(kRed); cutChio->SetLineWidth(2); + cutChio->Draw("same"); + cID->cd(2); + hPID->DrawClone("colz"); + cutPID->SetLineColor(kRed); cutPID->SetLineWidth(2); + cutPID->Draw("same"); + + //save to disk gated TTree + gated.Snapshot("BaseCutTree", "./RootFiles/BaseCutTree_12Be.root"); +} diff --git a/Projects/e748/Postanalysis/GetSpectrum.cpp b/Projects/e748/Postanalysis/GetSpectrum.cpp new file mode 100644 index 000000000..bff5942f2 --- /dev/null +++ b/Projects/e748/Postanalysis/GetSpectrum.cpp @@ -0,0 +1,62 @@ +#include "TROOT.h" +#include "TStyle.h" +#include "TFile.h" +#include "TTree.h" +#include "TChain.h" +#include "TCATSPhysics.h" +#include "TMust2Physics.h" +#include "NPReaction.h" + +#include "ROOT/RDataFrame.hxx" +#include <ROOT/RVec.hxx> + +#include "Utils.cpp" + +#include <Rtypes.h> +#include <TCanvas.h> +#include <algorithm> +#include <iostream> +#include <vector> + +void GetSpectrum() +{ + ROOT::EnableImplicitMT(); + gStyle->SetOptStat("nmeruoi"); + + ROOT::RDataFrame df("FinalTree", "./RootFiles/FinalTree_12Be.root"); + auto dsc {df.Describe()}; + dsc.Print(); + + //--> Read TCutGs + auto* cutChio {GetGraphCutFromFile("./Cuts/chio_v0.root")}; + auto* cutPID {GetGraphCutFromFile("./../TOF.root", "TOF")}; + + //--> Get histograms + // kinematic lines + auto hKin {df.Histo2D({"hKin", "Kinematic lines;#theta_{Lab} [degree];E_{Vertex} [MeV]", 500, 0., 60., 500, 0., 30.}, "ThetaLab", "ELab")}; + auto hBeam {df.Histo1D({"hBem", "Beam raw TOF", 500, 100., 300.}, "T_CATS1_CAV")}; + // excitation energy + auto hEex {df.Histo1D("Ex")}; + // beam energy + auto hBeamE {df.Histo1D("BeamEnergy")}; + + // Reaction + NPL::Reaction reaction("12Be(d,3He)11Li@355"); + auto* g {reaction.GetKinematicLine3()}; + + auto* cRes {new TCanvas("cRes", "Results canvas", 1)}; + cRes->DivideSquare(2); + cRes->cd(1); + hKin->DrawClone("colz"); + g->SetLineColor(kRed); g->SetLineWidth(2); + g->Draw("l same"); + cRes->cd(2); + hEex->DrawClone(); + + auto* cBeam {new TCanvas("cBeam", "Beam canvas", 1)}; + cBeam->DivideSquare(2); + cBeam->cd(1); + hBeam->DrawClone(); + cBeam->cd(2); + hBeamE->DrawClone(); +} diff --git a/Projects/e748/Postanalysis/ImprovePID.cpp b/Projects/e748/Postanalysis/ImprovePID.cpp new file mode 100644 index 000000000..c27f22b61 --- /dev/null +++ b/Projects/e748/Postanalysis/ImprovePID.cpp @@ -0,0 +1,85 @@ +#include <Math/Vector3D.h> +#include <Math/Point3D.h> +#include <ROOT/RDataFrame.hxx> +#include <ROOT/RVec.hxx> +#include <TCanvas.h> +#include <TMath.h> +#include <iostream> + +#include "TCATSPhysics.h" +#include "TMust2Data.h" +#include "TMust2Physics.h" + +using XYZVector = ROOT::Math::XYZVector; +using XYZPoint = ROOT::Math::XYZPoint; + +double LandTtoGamma(double L, double T) +{ + //MUST be in SI units + double v {L / T}; + double beta {v / TMath::C()}; + std::cout<<"Beta = "<<beta<<'\n'; + double gamma {1. / TMath::Sqrt(1. - beta * beta)}; + return gamma; +} + +double ConvertTOFToMass(double Tbeam, + double xM2, double yM2, double zM2, double tM2, + double tCorr, + TCATSPhysics& cats) +{ + double timeWindow {600.};// ns + double lengthFactor {180.};//mm + //Length + XYZPoint vertex {cats.PositionOnTargetX, cats.PositionOnTargetY, 0.}; + //get MUST2 point + XYZPoint must2Point {xM2, yM2, zM2 + lengthFactor}; + //Time (including correction) + double tof {tM2 + tCorr}; + std::cout<<"tCorr = "<<tCorr<<'\n'; + auto LAfterReaction {(must2Point - vertex).R()};//should be in mm! + //let's ignore by now distance cats (or whatever provides STOP signal) to vertex + auto L {(LAfterReaction) * 1.E-3};//mm to m + auto T {TMath::Abs(tof - timeWindow) * 1.E-9}; //ns to s + std::cout<<"L = "<<L<<" T = "<<T<<'\n'; + auto gamma {LandTtoGamma(L, T)}; + std::cout<<"Gamma = "<<gamma<<'\n'; + return Tbeam / (gamma - 1); +} + +void ImprovePID() +{ + ROOT::RDataFrame d("GatedTree", "./RootFiles/GatedTree_12Be.root"); + auto description {d.Describe()}; + description.Print();std::cout<<'\n'; + + auto df = d.Redefine("TOF", "TOF"); + auto hPID {df.Define("x", "MUST2.Si_E") + .Histo2D({"hPID", "PID;E_{Si} [MeV];TOF [au]", 1000, 0., 30., 1000, 460., 580.}, "x", "TOF")}; + df = df.Define("ReconstructedMass", ConvertTOFToMass, {"BeamEnergy", + "X_M2", "Y_M2", "Z_M2", "T_M2", + "TimeCorr", + "CATS"}); + + auto hMass {df.Histo1D( "ReconstructedMass")}; + + auto hTelN {df.Define("TelescopeNumber", [](const TMust2Physics& must2) + { + return must2.TelescopeNumber; + }, + {"MUST2"}) + .Histo1D("TelescopeNumber")}; + + //plotting + auto* c1 {new TCanvas("c1")}; + c1->DivideSquare(2); + c1->cd(1); + hPID->DrawClone("col"); + c1->cd(2); + hMass->DrawClone(); + + auto* c2 {new TCanvas("c2")}; + c2->DivideSquare(2); + c2->cd(1); + hTelN->DrawClone(); +} diff --git a/Projects/e748/Postanalysis/PlotTheoreticalKinematics.cpp b/Projects/e748/Postanalysis/PlotTheoreticalKinematics.cpp new file mode 100644 index 000000000..d2761e5c0 --- /dev/null +++ b/Projects/e748/Postanalysis/PlotTheoreticalKinematics.cpp @@ -0,0 +1,47 @@ +#include "NPReaction.h" +#include <Rtypes.h> +#include <TCanvas.h> +#include <TString.h> + +void PlotTheoreticalKinematics() +{ + double beam10Be {28. * 10.013534}; + double beam12Be {30. * 12.026922}; + //simply compare reaction kinematics + NPL::Reaction r10Be(TString::Format("10Be(d,3He)9Li@%.f", beam10Be).Data()); + NPL::Reaction r12Be (TString::Format("12Be(d,3He)11Li@%.f", beam12Be).Data()); + + //get graphs + auto* g9Li {r10Be.GetKinematicLine3()}; + auto* g11Li {r12Be.GetKinematicLine3()}; + //also Elab vs thetaCM + auto* gcm9Li {r10Be.GetELabVersusThetaCM()}; + auto* gcm11Li {r12Be.GetELabVersusThetaCM()}; + //theta3 vs theta4 in LAB + auto* gtheta9Li {r10Be.GetTheta3VsTheta4()}; + auto* gtheta11Li {r12Be.GetTheta3VsTheta4()}; + + //plotting + auto* c1 {new TCanvas("c1", "Kinematic lines comparaison")}; + c1->DivideSquare(4); + c1->cd(1); + g11Li->SetTitle(";#theta_{3}^{LAB} [degree];T_{3} [MeV]"); + g9Li->SetLineWidth(2); g9Li->SetLineColor(kRed); + g11Li->SetLineWidth(2); g11Li->SetLineColor(kBlue); + g11Li->Draw("apl"); + g9Li->Draw("pl same"); + + c1->cd(2); + gcm11Li->SetTitle(";#theta_{3}^{CM} [degree];T_{3} [MeV]"); + gcm9Li->SetLineWidth(2); gcm9Li->SetLineColor(kRed); + gcm11Li->SetLineWidth(2); gcm11Li->SetLineColor(kBlue); + gcm11Li->Draw("apl"); + gcm9Li->Draw("pl same"); + + c1->cd(3); + gtheta9Li->SetTitle(";#theta_{3}^{LAB} [degree];#theta_{4}^{LAB} [degree]"); + gtheta9Li->SetLineWidth(2); gtheta9Li->SetLineColor(kRed); + gtheta11Li->SetLineWidth(2); gtheta11Li->SetLineColor(kBlue); + gtheta9Li->Draw("apl"); + gtheta11Li->Draw("pl same"); +} diff --git a/Projects/e748/Postanalysis/Utils.cpp b/Projects/e748/Postanalysis/Utils.cpp new file mode 100644 index 000000000..78da8e96c --- /dev/null +++ b/Projects/e748/Postanalysis/Utils.cpp @@ -0,0 +1,35 @@ +#include <ROOT/RDataFrame.hxx> +#include <TCutG.h> +#include <TFile.h> +#include <TTree.h> +#include <TChain.h> + +#include <stdexcept> +#include <string> +#include <set> + +TCutG* GetGraphCutFromFile(const std::string& fileName, const std::string& name = "CUTG") +{ + auto* file {new TFile(fileName.c_str())}; + auto* cut {file->Get<TCutG>(name.c_str())}; + if(!cut) + throw std::runtime_error("Error loading TCutG: Check file existence or name!"); + return cut; + +} + +ROOT::RDataFrame ReadAll12BeData() +{ + auto* inFile {new TFile("/home/miguel/nptool/Projects/e748/output/analysis/Physics_12Be.root")}; + //inFile->ls(); + auto* tree {inFile->Get<TTree>("PhysicsTree")}; + //tree->Print(); + + //load also chio chain + auto* chioChain {new TChain("Numexo2")};//runs for 12Be + std::set<int> runs {315, 316, 317, 318, 320, 321, 323, 325, 326, 327, 328, 329, 330, 331, 339, 341, 342, 346, 347, 348}; + for(const auto& run : runs) + chioChain->Add(TString::Format("/home/miguel/nptool/Projects/e748/Data/Merged/run_%04d.root", run)); + tree->AddFriend(chioChain); + return ROOT::RDataFrame(*tree); +} -- GitLab