added basic postanalysis macros

Projects/e748/Postanalysis/AnalyzeBeam.cpp

+#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");
+}

Projects/e748/Postanalysis/AnalyzeCATS.cpp

+#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");
+  
+}

Projects/e748/Postanalysis/ApplyCuts.cpp

+#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");    
+}

Projects/e748/Postanalysis/CalibrateBeam.cpp

+#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();   
+}

Projects/e748/Postanalysis/GetCuts.cpp

+#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");
+}

Projects/e748/Postanalysis/GetSpectrum.cpp

+#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();
+}

Projects/e748/Postanalysis/ImprovePID.cpp

+#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();
+}

Projects/e748/Postanalysis/PlotTheoreticalKinematics.cpp

+#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");
+}

Projects/e748/Postanalysis/Utils.cpp

+#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);
+}