rtang/ANASEN_analysis
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

deleted: Analyzer1.C

Browse Source 
deleted:    Analyzer1.h
	deleted:    GainMatchSX3.C
	deleted:    GainMatchSX3.h
	deleted:    GainMatchSX3Front.C
	deleted:    GainMatchSX3Front.h
	deleted:    GainMatchSX3Front1.C
	deleted:    sx3_BackGains.txt
	deleted:    sx3_GainMatchfront.txt
    Book keeping deleting older unused scripts.
This commit is contained in:
Vignesh Sitaraman 2026-05-15 14:02:00 -04:00
parent 574506acb7
commit c142f4b567
11 changed files with 169 additions and 1854 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

402
Analyzer1.C
View File

@ -1,402 +0,0 @@
#define Analyzer1_cxx
#include "Analyzer1.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include <utility>
#include <algorithm>
#include "Armory/ClassSX3.h"
#include "Armory/ClassPW.h"
#include "TVector3.h"
TH2F * hsx3IndexVE;
TH2F * hqqqIndexVE;
TH2F * hpcIndexVE;
TH2F * hsx3Coin;
TH2F * hqqqCoin;
TH2F * hpcCoin;
TH2F * hqqqPolar;
TH2F * hsx3VpcIndex;
TH2F * hqqqVpcIndex;
TH2F * hqqqVpcE;
TH2F * hsx3VpcE;
TH2F * hanVScatsum;
int padID = 0;
SX3 sx3_contr;
PW pw_contr;
TVector3 hitPos;
bool HitNonZero;
TH1F * hZProj;
void Analyzer1::Begin(TTree * /*tree*/){
TString option = GetOption();
hsx3IndexVE = new TH2F("hsx3IndexVE", "SX3 index vs Energy; sx3 index ; Energy", 24*12, 0, 24*12, 400, 0, 5000); hsx3IndexVE->SetNdivisions( -612, "x");
hqqqIndexVE = new TH2F("hqqqIndexVE", "QQQ index vs Energy; QQQ index ; Energy", 4*2*16, 0, 4*2*16, 400, 0, 5000); hqqqIndexVE->SetNdivisions( -1204, "x");
hpcIndexVE = new TH2F("hpcIndexVE", "PC index vs Energy; PC index ; Energy", 2*24, 0, 2*24, 400, 0, 4000); hpcIndexVE->SetNdivisions( -1204, "x");
hsx3Coin = new TH2F("hsx3Coin", "SX3 Coincident", 24*12, 0, 24*12, 24*12, 0, 24*12);
hqqqCoin = new TH2F("hqqqCoin", "QQQ Coincident", 4*2*16, 0, 4*2*16, 4*2*16, 0, 4*2*16);
hpcCoin = new TH2F("hpcCoin", "PC Coincident", 2*24, 0, 2*24, 2*24, 0, 2*24);
hqqqPolar = new TH2F("hqqqPolar", "QQQ Polar ID", 16*4, -TMath::Pi(), TMath::Pi(),16, 10, 50);
hsx3VpcIndex = new TH2F("hsx3Vpcindex", "sx3 vs pc; sx3 index; pc index", 24*12, 0, 24*12, 48, 0, 48);
hsx3VpcIndex->SetNdivisions( -612, "x");
hsx3VpcIndex->SetNdivisions( -12, "y");
hqqqVpcIndex = new TH2F("hqqqVpcindex", "qqq vs pc; qqq index; pc index", 4*2*16, 0, 4*2*16, 48, 0, 48);
hqqqVpcIndex->SetNdivisions( -612, "x");
hqqqVpcIndex->SetNdivisions( -12, "y");
hqqqVpcE = new TH2F("hqqqVpcEnergy", "qqq vs pc; qqq energy; pc energy", 400, 0, 5000, 400, 0, 5000);
hqqqVpcE->SetNdivisions( -612, "x");
hqqqVpcE->SetNdivisions( -12, "y");
hsx3VpcE = new TH2F("hsx3VpcEnergy", "sx3 vs pc; sx3 energy; pc energy", 400, 0, 5000, 400, 0, 5000);
hsx3VpcE->SetNdivisions( -612, "x");
hsx3VpcE->SetNdivisions( -12, "y");
hZProj = new TH1F("hZProj", "Z Projection", 1200, -600, 600);
hanVScatsum = new TH2F("hanVScatsum", "Anode vs Cathode Sum; Anode E; Cathode E", 400,0 , 10000, 400, 0 , 16000);
sx3_contr.ConstructGeo();
pw_contr.ConstructGeo();
}
Bool_t Analyzer1::Process(Long64_t entry){
// if ( entry > 100 ) return kTRUE;
hitPos.Clear();
HitNonZero = false;
// if( entry > 1) return kTRUE;
// printf("################### ev : %llu \n", entry);
b_sx3Multi->GetEntry(entry);
b_sx3ID->GetEntry(entry);
b_sx3Ch->GetEntry(entry);
b_sx3E->GetEntry(entry);
b_sx3T->GetEntry(entry);
b_qqqMulti->GetEntry(entry);
b_qqqID->GetEntry(entry);
b_qqqCh->GetEntry(entry);
b_qqqE->GetEntry(entry);
b_qqqT->GetEntry(entry);
b_pcMulti->GetEntry(entry);
b_pcID->GetEntry(entry);
b_pcCh->GetEntry(entry);
b_pcE->GetEntry(entry);
b_pcT->GetEntry(entry);
sx3.CalIndex();
qqq.CalIndex();
pc.CalIndex();
// sx3.Print();
//########################################################### Raw data
// //======================= SX3
std::vector<std::pair<int, int>> ID; // first = id, 2nd = index
for( int i = 0; i < sx3.multi; i ++){
ID.push_back(std::pair<int, int>(sx3.id[i], i));
hsx3IndexVE->Fill( sx3.index[i], sx3.e[i] );
for( int j = i+1; j < sx3.multi; j++){
hsx3Coin->Fill( sx3.index[i], sx3.index[j]);
}
for( int j = 0; j < pc.multi; j++){
hsx3VpcIndex->Fill( sx3.index[i], pc.index[j] );
// if( sx3.ch[index] > 8 ){
// hsx3VpcE->Fill( sx3.e[i], pc.e[j] );
// }
}
}
if( ID.size() > 0 ){
std::sort(ID.begin(), ID.end(), [](const std::pair<int, int> & a, const std::pair<int, int> & b) {
return a.first < b.first;
} );
// printf("##############################\n");
// for( size_t i = 0; i < ID.size(); i++) printf("%zu | %d %d \n", i, ID[i].first, ID[i].second );
std::vector<std::pair<int, int>> sx3ID;
sx3ID.push_back(ID[0]);
bool found = false;
for( size_t i = 1; i < ID.size(); i++){
if( ID[i].first == sx3ID.back().first) {
sx3ID.push_back(ID[i]);
if( sx3ID.size() >= 3) {
found = true;
}
}else{
if( !found ){
sx3ID.clear();
sx3ID.push_back(ID[i]);
}
}
}
// printf("---------- sx3ID Multi : %zu \n", sx3ID.size());
if( found ){
int sx3ChUp, sx3ChDn, sx3ChBk;
float sx3EUp, sx3EDn;
// printf("------ sx3 ID : %d, multi: %zu\n", sx3ID[0].first, sx3ID.size());
for( size_t i = 0; i < sx3ID.size(); i++ ){
int index = sx3ID[i].second;
// printf(" %zu | index %d | ch : %d, energy : %d \n", i, index, sx3.ch[index], sx3.e[index]);
if( sx3.ch[index] < 8 ){
if( sx3.ch[index] % 2 == 0) {
sx3ChDn = sx3.ch[index];
sx3EDn = sx3.e[index];
}else{
sx3ChUp = sx3.ch[index];
sx3EUp = sx3.e[index];
}
}else{
sx3ChBk = sx3.ch[index];
}
for( int j = 0; j < pc.multi; j++){
// hsx3VpcIndex->Fill( sx3.index[i], pc.index[j] );
if( sx3.ch[index] > 8 ){
hsx3VpcE->Fill( sx3.e[i], pc.e[j] );
// hpcIndexVE->Fill( pc.index[i], pc.e[i] );
}
}
}
sx3_contr.CalSX3Pos(sx3ID[0].first, sx3ChUp, sx3ChDn, sx3ChBk, sx3EUp, sx3EDn);
hitPos = sx3_contr.GetHitPos();
HitNonZero = true;
// hitPos.Print();
}
}
// //======================= QQQ
for( int i = 0; i < qqq.multi; i ++){
// for( int j = 0; j < pc.multi; j++){
// if(pc.index[j]==4){
hqqqIndexVE->Fill( qqq.index[i], qqq.e[i] );
// }
// }
for( int j = 0; j < qqq.multi; j++){
if ( j == i ) continue;
hqqqCoin->Fill( qqq.index[i], qqq.index[j]);
}
for( int j = i + 1; j < qqq.multi; j++){
for( int k = 0; k < pc.multi; k++){
if(pc.index[k]<24 && pc.e[k]>50 ){
hqqqVpcE->Fill( qqq.e[i], pc.e[k] );
// hpcIndexVE->Fill( pc.index[i], pc.e[i] );
hqqqVpcIndex->Fill( qqq.index[i], pc.index[j] );
}
// }
}
// if( qqq.used[i] == true ) continue;
//if( qqq.id[i] == qqq.id[j] && (16 - qqq.ch[i]) * (16 - qqq.ch[j]) < 0 ){ // must be same detector and wedge and ring
if( qqq.id[i] == qqq.id[j] ){ // must be same detector
int chWedge = -1;
int chRing = -1;
if( qqq.ch[i] < qqq.ch[j]){
chRing = qqq.ch[j] - 16;
chWedge = qqq.ch[i];
}else{
chRing = qqq.ch[i];
chWedge = qqq.ch[j] - 16;
}
// printf(" ID : %d , chWedge : %d, chRing : %d \n", qqq.id[i], chWedge, chRing);
double theta = -TMath::Pi()/2 + 2*TMath::Pi()/16/4.*(qqq.id[i]*16 + chWedge +0.5);
double rho = 10.+40./16.*(chRing+0.5);
// if(qqq.e[i]>50){
hqqqPolar->Fill( theta, rho);
// }
// qqq.used[i] = true;
// qqq.used[j] = true;
if( !HitNonZero ){
double x = rho * TMath::Cos(theta);
double y = rho * TMath::Sin(theta);
hitPos.SetXYZ(x, y, 23 + 75 + 30);
HitNonZero = true;
}
}
}
}
// //======================= PC
ID.clear();
int counter=0;
std::vector<std::pair<int, double>> E;
E.clear();
for( int i = 0; i < pc.multi; i ++){
if( pc.e[i] > 100 ) ID.push_back(std::pair<int, int>(pc.id[i], i));
if( pc.e[i] > 100 ) E.push_back(std::pair<int, double>(pc.index[i], pc.e[i]));
hpcIndexVE->Fill( pc.index[i], pc.e[i] );
for( int j = i+1; j < pc.multi; j++){
hpcCoin->Fill( pc.index[i], pc.index[j]);
}
}
// for( size_t i = 0; i < E.size(); i++) printf("%zu | %d %d \n", i, E[i].first, E[i].second );
if( E.size()>=3 ){
int aID = 0;
int cID = 0;
float aE = 0;
float cE = 0;
bool multi_an =false;
// if( ID[0].first < 1 ) {
// aID = pc.ch[ID[0].second];
// cID = pc.ch[ID[1].second];
// }else{
// cID = pc.ch[ID[0].second];
// aID = pc.ch[ID[1].second];
// }
// printf("anode= %d, cathode = %d\n", aID, cID);
// for( int k = 0; k < qqq.multi; k++){
// if(qqq.index[k]==75 && pc.index[k]==2 && pc.e[k]>100){
for(int l=0;l<E.size();l++){
if(E[l].first<24 ){
if(!multi_an){
aE = E[l].second;
}
multi_an=true;
}
else if (E[l].first>=24){
cE = E[l].second + cE;
}
}
// }
// }
hanVScatsum->Fill(aE,cE);
if( ID[0].first < 1 ) {
aID = pc.ch[ID[0].second];
cID = pc.ch[ID[1].second];
}else{
cID = pc.ch[ID[0].second];
aID = pc.ch[ID[1].second];
}
if( HitNonZero){
pw_contr.CalTrack( hitPos, aID, cID);
hZProj->Fill(pw_contr.GetZ0());
}
}
//########################################################### Track constrcution
//############################## DO THE KINEMATICS
return kTRUE;
}
void Analyzer1::Terminate(){
gStyle->SetOptStat("neiou");
TCanvas * canvas = new TCanvas("cANASEN", "ANASEN", 2000, 2000);
canvas->Divide(3,3);
//hsx3VpcIndex->Draw("colz");
//=============================================== pad-1
padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
hsx3IndexVE->Draw("colz");
//=============================================== pad-2
padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
hqqqIndexVE->Draw("colz");
//=============================================== pad-3
padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
hpcIndexVE->Draw("colz");
//=============================================== pad-4
padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
hsx3Coin->Draw("colz");
//=============================================== pad-5
padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
canvas->cd(padID)->SetLogz(true);
hqqqCoin->Draw("colz");
//=============================================== pad-6
padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
hpcCoin->Draw("colz");
//=============================================== pad-7
padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
// hsx3VpcIndex ->Draw("colz");
hsx3VpcE->Draw("colz") ;
//=============================================== pad-8
padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
// hqqqVpcIndex ->Draw("colz");
hqqqVpcE ->Draw("colz");
//=============================================== pad-9
padID ++;
// canvas->cd(padID)->DrawFrame(-50, -50, 50, 50);
// hqqqPolar->Draw("same colz pol");
canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
// hZProj->Draw();
hanVScatsum->Draw("colz");
}

114
Analyzer1.h
View File

@ -1,114 +0,0 @@
#ifndef Analyzer1_h
#define Analyzer1_h
#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
#include <TSelector.h>
#include "Armory/ClassDet.h"
class Analyzer1 : public TSelector {
public :
TTree *fChain; //!pointer to the analyzed TTree or TChain
// Fixed size dimensions of array or collections stored in the TTree if any.
// Declaration of leaf types
Det sx3;
Det qqq;
Det pc ;
ULong64_t evID;
UInt_t run;
// List of branches
TBranch *b_eventID; //!
TBranch *b_run; //!
TBranch *b_sx3Multi; //!
TBranch *b_sx3ID; //!
TBranch *b_sx3Ch; //!
TBranch *b_sx3E; //!
TBranch *b_sx3T; //!
TBranch *b_qqqMulti; //!
TBranch *b_qqqID; //!
TBranch *b_qqqCh; //!
TBranch *b_qqqE; //!
TBranch *b_qqqT; //!
TBranch *b_pcMulti; //!
TBranch *b_pcID; //!
TBranch *b_pcCh; //!
TBranch *b_pcE; //!
TBranch *b_pcT; //!
Analyzer1(TTree * /*tree*/ =0) : fChain(0) { }
virtual ~Analyzer1() { }
virtual Int_t Version() const { return 2; }
virtual void Begin(TTree *tree);
virtual void SlaveBegin(TTree *tree);
virtual void Init(TTree *tree);
virtual Bool_t Notify();
virtual Bool_t Process(Long64_t entry);
virtual Int_t GetEntry(Long64_t entry, Int_t getall = 0) { return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0; }
virtual void SetOption(const char *option) { fOption = option; }
virtual void SetObject(TObject *obj) { fObject = obj; }
virtual void SetInputList(TList *input) { fInput = input; }
virtual TList *GetOutputList() const { return fOutput; }
virtual void SlaveTerminate();
virtual void Terminate();
ClassDef(Analyzer1,0);
};
#endif
#ifdef Analyzer1_cxx
void Analyzer1::Init(TTree *tree){
// Set branch addresses and branch pointers
if (!tree) return;
fChain = tree;
fChain->SetMakeClass(1);
fChain->SetBranchAddress("evID", &evID, &b_eventID);
fChain->SetBranchAddress("run", &run, &b_run);
sx3.SetDetDimension(24,12);
qqq.SetDetDimension(4,32);
pc.SetDetDimension(2,24);
fChain->SetBranchAddress("sx3Multi", &sx3.multi, &b_sx3Multi);
fChain->SetBranchAddress("sx3ID", &sx3.id, &b_sx3ID);
fChain->SetBranchAddress("sx3Ch", &sx3.ch, &b_sx3Ch);
fChain->SetBranchAddress("sx3E", &sx3.e, &b_sx3E);
fChain->SetBranchAddress("sx3T", &sx3.t, &b_sx3T);
fChain->SetBranchAddress("qqqMulti", &qqq.multi, &b_qqqMulti);
fChain->SetBranchAddress("qqqID", &qqq.id, &b_qqqID);
fChain->SetBranchAddress("qqqCh", &qqq.ch, &b_qqqCh);
fChain->SetBranchAddress("qqqE", &qqq.e, &b_qqqE);
fChain->SetBranchAddress("qqqT", &qqq.t, &b_qqqT);
fChain->SetBranchAddress("pcMulti", &pc.multi, &b_pcMulti);
fChain->SetBranchAddress("pcID", &pc.id, &b_pcID);
fChain->SetBranchAddress("pcCh", &pc.ch, &b_pcCh);
fChain->SetBranchAddress("pcE", &pc.e, &b_pcE);
fChain->SetBranchAddress("pcT", &pc.t, &b_pcT);
}
Bool_t Analyzer1::Notify(){
return kTRUE;
}
void Analyzer1::SlaveBegin(TTree * /*tree*/){
TString option = GetOption();
}
void Analyzer1::SlaveTerminate(){
}
#endif // #ifdef Analyzer_cxx

432
GainMatchSX3.C
View File

@ -1,432 +0,0 @@
#define GainMatchSX3_cxx
#include "GainMatchSX3.h"
#include <TH2.h>
#include <TF1.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TCutG.h>
#include <fstream>
#include <utility>
#include <algorithm>
#include <TProfile.h>
#include "Armory/ClassSX3.h"
#include "Armory/HistPlotter.h"
#include <TGraphErrors.h>
#include "TVector3.h"
TH2F *hSX3FvsB;
TH2F *hSX3FvsB_g;
TH2F *hsx3IndexVE;
TH2F *hsx3IndexVE_g;
TH2F *hSX3;
TH2F *hsx3Coin;
int padID = 0;
SX3 sx3_contr;
TCutG *cut;
TCutG *cut1;
std::map<std::tuple<int, int, int, int>, std::vector<std::tuple<double, double, double>>> dataPoints;
std::map<std::tuple<int, int, int, int>, int> comboCounts;
const int MAX_DET = 24;
const int MAX_UP = 4;
const int MAX_DOWN = 4;
const int MAX_BK = 4;
double frontGainUp[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
double frontGainDown[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
bool frontGainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}};
TCanvas c("canvas", "canvas", 800, 600);
// ==== Configuration Flags ====
const bool interactiveMode = true; // If true: show canvas + wait for user
const bool verboseFit = true; // If true: print fit summary and chi²
const bool drawCanvases = true; // If false: canvases won't be drawn at all
// HistPlotter plotter("SX3GainMatchBack.root");
void GainMatchSX3::Begin(TTree * /*tree*/)
{
TString option = GetOption();
hSX3FvsB = new TH2F("hSX3FvsB", "SX3 Front vs Back; Front E; Back E", 400, 0, 16000, 400, 0, 16000);
hSX3FvsB_g = new TH2F("hSX3FvsB_g", "SX3 Front vs Back; Front E; Back E", 400, 0, 16000, 400, 0, 16000);
hsx3IndexVE = new TH2F("hsx3IndexVE", "SX3 index vs Energy; sx3 index ; Energy", 24 * 12, 0, 24 * 12, 400, 0, 5000);
hsx3IndexVE_g = new TH2F("hsx3IndexVE_g", "SX3 index vs Energy; sx3 index ; Energy", 24 * 12, 0, 24 * 12, 400, 0, 5000);
hSX3 = new TH2F("hSX3", "SX3 Front v Back; Fronts; Backs", 8, 0, 8, 4, 0, 4);
hsx3Coin = new TH2F("hsx3Coin", "SX3 Coincident", 24 * 12, 0, 24 * 12, 24 * 12, 0, 24 * 12);
sx3_contr.ConstructGeo();
// Load the TCutG object
TFile *cutFile = TFile::Open("sx3cut.root");
if (!cutFile || cutFile->IsZombie())
{
std::cerr << "Error: Could not open sx3cut.root" << std::endl;
return;
}
cut = dynamic_cast<TCutG *>(cutFile->Get("sx3cut"));
if (!cut)
{
std::cerr << "Error: Could not find TCutG named 'sx3cut' in sx3cut.root" << std::endl;
return;
}
cut->SetName("sx3cut"); // Ensure the cut has the correct name
// Load the TCutG object
TFile *cutFile1 = TFile::Open("UvD.root");
bool cut1Loaded = (cut1 != nullptr);
cut1 = dynamic_cast<TCutG *>(cutFile1->Get("UvD"));
if (!cut1)
{
std::cerr << "Error: Could not find TCutG named 'UvD' in UvD.root" << std::endl;
return;
}
cut1->SetName("UvD");
// plotter.ReadCuts("cuts.txt");
std::string filename = "sx3_GainMatchfront.txt";
// std::string filename = "sx3_GainMatchfront.txt";
std::ifstream infile(filename);
if (!infile.is_open())
{
std::cerr << "Error opening " << filename << "!" << std::endl;
return;
}
int id, bk, u, d;
double gainup, gaindown;
while (infile >> id >> bk >> u >> d >> gainup >> gaindown)
{
frontGainUp[id][bk][u][d] = gainup;
frontGainDown[id][bk][u][d] = gaindown;
frontGainValid[id][bk][u][d] = true;
if(frontGainValid[id][bk][u][d]) {
// std::cout << "Loaded front gain for Det" << id << " Bk" << bk << " U" << u << " D" << d
// << ": Up=" << gainup << ", Down=" << gaindown << std::endl;
}
else {
std::cout << "No valid front gain for Det" << id << " Bk" << bk << " U" << u << " D" << d << std::endl;
}
}
}
Bool_t GainMatchSX3::Process(Long64_t entry)
{
b_sx3Multi->GetEntry(entry);
b_sx3ID->GetEntry(entry);
b_sx3Ch->GetEntry(entry);
b_sx3E->GetEntry(entry);
b_sx3T->GetEntry(entry);
b_qqqMulti->GetEntry(entry);
b_qqqID->GetEntry(entry);
b_qqqCh->GetEntry(entry);
b_qqqE->GetEntry(entry);
b_qqqT->GetEntry(entry);
b_pcMulti->GetEntry(entry);
b_pcID->GetEntry(entry);
b_pcCh->GetEntry(entry);
b_pcE->GetEntry(entry);
b_pcT->GetEntry(entry);
sx3.CalIndex();
qqq.CalIndex();
pc.CalIndex();
std::vector<std::pair<int, int>> ID;
for (int i = 0; i < sx3.multi; i++)
{
// for (int j = i + 1; j < sx3.multi; j++)
// {
// if (sx3.id[i] == 3)
// hsx3Coin->Fill(sx3.index[i], sx3.index[j]);
// }
if (sx3.e[i] > 100)
{
ID.push_back(std::pair<int, int>(sx3.id[i], i));
hsx3IndexVE->Fill(sx3.index[i], sx3.e[i]);
}
}
if (ID.size() > 0)
{
std::sort(ID.begin(), ID.end(), [](const std::pair<int, int> &a, const std::pair<int, int> &b)
{ return a.first < b.first; });
// start with the first entry in the sorted array: channels that belong to the same detector are together in sequenmce
std::vector<std::pair<int, int>> sx3ID;
sx3ID.push_back(ID[0]);
bool found = false;
for (size_t i = 1; i < ID.size(); i++)
{ // Check if id of i belongs to the same detector and then add it to the detector ID vector
if (ID[i].first == sx3ID.back().first)
{ // count the nunmber of hits that belong to the same detector
sx3ID.push_back(ID[i]);
if (sx3ID.size() >= 3)
{
found = true;
}
}
else
{ // the next event does not belong to the same detector, abandon the first event and continue with the next one
if (!found)
{
sx3ID.clear();
sx3ID.push_back(ID[i]);
}
}
}
if (found)
{
int sx3ChUp = -1, sx3ChDn = -1, sx3ChBk = -1;
float sx3EUp = 0.0, sx3EDn = 0.0, sx3EBk = 0.0;
// Build the correlated set once
for (size_t i = 0; i < sx3ID.size(); i++)
{
if (sx3.e[i] > 100)
{
int index = sx3ID[i].second;
if (sx3.ch[index] < 8)
{
if (sx3.ch[index] % 2 == 0)
{
sx3ChDn = sx3.ch[index];
sx3EDn = sx3.e[index];
//
}
else
{
sx3ChUp = sx3.ch[index];
sx3EUp = sx3.e[index];
}
}
else
{
sx3ChBk = sx3.ch[index] - 8;
sx3EBk = sx3.e[index];
}
}
}
// Only if we found all three channels do we proceed
if (sx3ChUp >= 0 && sx3ChDn >= 0 && sx3ChBk >= 0)
{
// Fill once per correlated set
hSX3->Fill(sx3ChDn + 4, sx3ChBk);
hSX3->Fill(sx3ChUp, sx3ChBk);
hSX3FvsB->Fill(sx3EUp + sx3EDn, sx3EBk);
if (frontGainValid[sx3ID[0].first][sx3ChBk][sx3ChUp / 2][sx3ChDn / 2])
{
sx3EUp *= frontGainUp[sx3ID[0].first][sx3ChBk][sx3ChUp / 2][sx3ChDn / 2];
sx3EDn *= frontGainDown[sx3ID[0].first][sx3ChBk][sx3ChUp / 2][sx3ChDn / 2];
}
else
{
// printf("No front gain for Det%d Bk%d U%d D%d\n", sx3ID[0].first, sx3ChBk, sx3ChUp / 2, sx3ChDn / 2);
sx3EUp = sx3EDn = 0.;
}
// plotter.Fill2D("hSX3F", 400, 0, 16000, 400, 0, 16000, sx3EUp + sx3EDn, sx3EBk);
// Pick detector ID from one of the correlated hits (all same detector)
int detID = sx3ID[0].first;
TString histName = Form("hSX3FVB_id%d_U%d_D%d_B%d", detID, sx3ChUp, sx3ChDn, sx3ChBk);
TString histName1 = Form("UnCorr_id%d_U%d-D%dvsB%d", detID, sx3ChUp, sx3ChDn, sx3ChBk);
TH2F *hist2d = (TH2F *)gDirectory->Get(histName);
TH2F *hist2d1 = (TH2F *)gDirectory->Get(histName1);
if (!hist2d)
{
hist2d = new TH2F(histName, histName,
400, 0, 16000, 400, 0, 16000);
}
if (!hist2d1)
{
hist2d1 = new TH2F(histName1, histName1,
800, -1, 1, 800, 0, 4000);
}
if (sx3EBk > 100 || sx3EUp > 100 || sx3EDn > 100)
{
hSX3FvsB_g->Fill(sx3EUp + sx3EDn, sx3EBk);
// Use the correlated triplet directly
dataPoints[{detID, sx3ChBk, sx3ChUp, sx3ChDn}]
.emplace_back(sx3EBk, sx3EUp, sx3EDn);
}
hist2d->Fill(sx3EUp + sx3EDn, sx3EBk);
hist2d1->Fill((sx3EUp - sx3EDn) / (sx3EUp + sx3EDn), sx3EBk);
}
}
}
return kTRUE;
}
const double GAIN_ACCEPTANCE_THRESHOLD = 0.3;
void GainMatchSX3::Terminate()
{
double backSlope[MAX_DET][MAX_BK] = {{0}};
bool backSlopeValid[MAX_DET][MAX_BK] = {{false}};
std::ofstream outFile("sx3_BackGains0.txt");
if (!outFile.is_open())
{
std::cerr << "Error opening sx3_BackGains.txt for writing!" << std::endl;
return;
}
// === Gain fit: (Up+Dn) vs Back, grouped by [id][bk] ===
for (int id = 0; id < MAX_DET; id++)
{
for (int bk = 0; bk < MAX_BK; bk++)
{
std::vector<double> bkE, udE;
// Collect all (Up+Dn, Back) for this id,bk
for (const auto &kv : dataPoints)
{
auto [cid, cbk, u, d] = kv.first;
if (cid != id || cbk != bk)
continue;
for (const auto &pr : kv.second)
{
double eBk, eUp, eDn;
std::tie(eBk, eUp, eDn) = pr;
if ((eBk < 100) || (eUp < 100) || (eDn < 100))
continue;
bkE.push_back(eBk);
udE.push_back(eUp + eDn);
}
}
if (bkE.size() < 5)
continue; // not enough statistics
// Build graph with errors
const double fixedError = 0.0; // ADC channels
std::vector<double> exVals(udE.size(), 0.0); // no x error
std::vector<double> eyVals(udE.size(), fixedError); // constant y error
TGraphErrors g(udE.size(), udE.data(), bkE.data(),
exVals.data(), eyVals.data());
TF1 f("f", "[0]*x", 0, 16000);
// f.SetParameter(0, 1.0); // initial slope
if (drawCanvases)
{
g.SetTitle(Form("Detector %d Back %d: (Up+Dn) vs Back", id, bk));
g.SetMarkerStyle(20);
g.SetMarkerColor(kBlue);
g.Draw("AP");
g.Fit(&f, interactiveMode ? "Q" : "QNR");
if (verboseFit)
{
double chi2 = f.GetChisquare();
int ndf = f.GetNDF();
double reducedChi2 = (ndf != 0) ? chi2 / ndf : -1;
std::cout << Form("Det%d Back%d → Slope: %.4f | χ²/ndf = %.2f/%d = %.2f",
id, bk, f.GetParameter(0), chi2, ndf, reducedChi2)
<< std::endl;
}
if (interactiveMode)
{
c.Update();
gPad->WaitPrimitive();
}
else
{
c.Close();
}
}
else
{
g.Fit(&f, "QNR");
}
double slope = 1 / f.GetParameter(0);
if (std::abs(slope - 1.0) < 0.3) // sanity check
{
backSlope[id][bk] = slope;
backSlopeValid[id][bk] = true;
outFile << id << " " << bk << " " << slope << "\n";
printf("Back slope Det%d Bk%d → %.4f\n", id, bk, slope);
}
else
{
std::cerr << "Warning: Bad slope for Det" << id << " Bk" << bk
<< " slope=" << slope << std::endl;
}
}
}
outFile.close();
std::cout << "Back gain matching complete." << std::endl;
// === Create histograms ===
TH2F *hFVB = new TH2F("hFVB", "Corrected Up+Dn vs Corrected Back;Up+Dn E;Corrected Back E",
600, 0, 16000, 600, 0, 16000);
TH2F *hAsym = new TH2F("hAsym", "Up vs Dn divide corrected back;Up/Back E;Dn/Back E",
400, 0.0, 1.0, 400, 0.0, 1.0);
TH2F *hAsymUnorm = new TH2F("hAsymUnorm", "Up vs Dn;Up E;Dn E",
800, 0.0, 4000.0, 800, 0.0, 4000.0);
// Fill histograms using corrected back energies
for (const auto &kv : dataPoints)
{
auto [id, bk, u, d] = kv.first;
if (!backSlopeValid[id][bk])
continue;
double slope = backSlope[id][bk];
for (const auto &pr : kv.second)
{
double eBk, eUp, eDn;
std::tie(eBk, eUp, eDn) = pr;
double updn = eUp + eDn;
if (updn == 0 || eBk == 0)
continue;
double correctedBack = eBk * slope;
double asym = (eUp - eDn) / updn;
hFVB->Fill(updn, correctedBack);
hAsym->Fill(eUp / correctedBack, eDn / correctedBack);
hAsymUnorm->Fill(eUp, eDn);
TString histNamex = Form("CorrBack_id%d_U%d-D%dvsB%d", id, u, d, bk);
TH2F *hist2dx = (TH2F *)gDirectory->Get(histNamex);
if (!hist2dx)
{
hist2dx = new TH2F(histNamex, histNamex,
800, -1, 1, 800, 0, 4000);
}
hist2dx->Fill((eUp - eDn) / (eUp + eDn), correctedBack);
}
}
// plotter.FlushToDisk();
}

114
GainMatchSX3.h
View File

@ -1,114 +0,0 @@
#ifndef GainMatchSX3_h
#define GainMatchSX3_h
#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
#include <TSelector.h>
#include "Armory/ClassDet.h"
class GainMatchSX3 : public TSelector {
public :
TTree *fChain; //!pointer to the analyzed TTree or TChain
// Fixed size dimensions of array or collections stored in the TTree if any.
// Declaration of leaf types
Det sx3;
Det qqq;
Det pc ;
ULong64_t evID;
UInt_t run;
// List of branches
TBranch *b_eventID; //!
TBranch *b_run; //!
TBranch *b_sx3Multi; //!
TBranch *b_sx3ID; //!
TBranch *b_sx3Ch; //!
TBranch *b_sx3E; //!
TBranch *b_sx3T; //!
TBranch *b_qqqMulti; //!
TBranch *b_qqqID; //!
TBranch *b_qqqCh; //!
TBranch *b_qqqE; //!
TBranch *b_qqqT; //!
TBranch *b_pcMulti; //!
TBranch *b_pcID; //!
TBranch *b_pcCh; //!
TBranch *b_pcE; //!
TBranch *b_pcT; //!
GainMatchSX3(TTree * /*tree*/ =0) : fChain(0) { }
virtual ~GainMatchSX3() { }
virtual Int_t Version() const { return 2; }
virtual void Begin(TTree *tree);
virtual void SlaveBegin(TTree *tree);
virtual void Init(TTree *tree);
virtual Bool_t Notify();
virtual Bool_t Process(Long64_t entry);
virtual Int_t GetEntry(Long64_t entry, Int_t getall = 0) { return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0; }
virtual void SetOption(const char *option) { fOption = option; }
virtual void SetObject(TObject *obj) { fObject = obj; }
virtual void SetInputList(TList *input) { fInput = input; }
virtual TList *GetOutputList() const { return fOutput; }
virtual void SlaveTerminate();
virtual void Terminate();
ClassDef(GainMatchSX3,0);
};
#endif
#ifdef GainMatchSX3_cxx
void GainMatchSX3::Init(TTree *tree){
// Set branch addresses and branch pointers
if (!tree) return;
fChain = tree;
fChain->SetMakeClass(1);
fChain->SetBranchAddress("evID", &evID, &b_eventID);
fChain->SetBranchAddress("run", &run, &b_run);
sx3.SetDetDimension(24,12);
qqq.SetDetDimension(4,32);
pc.SetDetDimension(2,24);
fChain->SetBranchAddress("sx3Multi", &sx3.multi, &b_sx3Multi);
fChain->SetBranchAddress("sx3ID", &sx3.id, &b_sx3ID);
fChain->SetBranchAddress("sx3Ch", &sx3.ch, &b_sx3Ch);
fChain->SetBranchAddress("sx3E", &sx3.e, &b_sx3E);
fChain->SetBranchAddress("sx3T", &sx3.t, &b_sx3T);
fChain->SetBranchAddress("qqqMulti", &qqq.multi, &b_qqqMulti);
fChain->SetBranchAddress("qqqID", &qqq.id, &b_qqqID);
fChain->SetBranchAddress("qqqCh", &qqq.ch, &b_qqqCh);
fChain->SetBranchAddress("qqqE", &qqq.e, &b_qqqE);
fChain->SetBranchAddress("qqqT", &qqq.t, &b_qqqT);
fChain->SetBranchAddress("pcMulti", &pc.multi, &b_pcMulti);
fChain->SetBranchAddress("pcID", &pc.id, &b_pcID);
fChain->SetBranchAddress("pcCh", &pc.ch, &b_pcCh);
fChain->SetBranchAddress("pcE", &pc.e, &b_pcE);
fChain->SetBranchAddress("pcT", &pc.t, &b_pcT);
}
Bool_t GainMatchSX3::Notify(){
return kTRUE;
}
void GainMatchSX3::SlaveBegin(TTree * /*tree*/){
TString option = GetOption();
}
void GainMatchSX3::SlaveTerminate(){
}
#endif // #ifdef GainMatchSX3_cxx

420
GainMatchSX3Front.C
View File

@ -1,420 +0,0 @@
#define GainMatchSX3Front_cxx
#include "GainMatchSX3Front.h"
#include <TH2.h>
#include <TF1.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TCutG.h>
#include <fstream>
#include <utility>
#include <algorithm>
#include <TProfile.h>
#include "Armory/ClassSX3.h"
#include "TGraphErrors.h"
#include "TMultiDimFit.h"
#include "TVector3.h"
TH2F *hSX3FvsB;
TH2F *hSX3FvsB_g;
TH2F *hsx3IndexVE;
TH2F *hsx3IndexVE_g;
TH2F *hSX3;
TH2F *hsx3Coin;
int padID = 0;
SX3 sx3_contr;
TCutG *cut;
TCutG *cut1;
std::map<std::tuple<int, int, int, int>, std::vector<std::tuple<double, double, double>>> dataPoints;
TCanvas c(Form("canvas"), "Fit", 800, 600);
// Gain arrays
const int MAX_DET = 24;
const int MAX_UP = 4;
const int MAX_DOWN = 4;
const int MAX_BK = 4;
double backGain[MAX_DET][MAX_BK] = {{0}};
bool backGainValid[MAX_DET][MAX_BK] = {{false}};
double frontGain[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
bool frontGainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}};
double uvdslope[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
// ==== Configuration Flags ====
const bool interactiveMode = true; // If true: show canvas + wait for user
const bool verboseFit = true; // If true: print fit summary and chi²
const bool drawCanvases = true; // If false: canvases won't be drawn at all
void GainMatchSX3Front::Begin(TTree * /*tree*/)
{
TString option = GetOption();
hSX3FvsB = new TH2F("hSX3FvsB", "SX3 Front vs Back; Front E; Back E", 800, 0, 16000, 800, 0, 16000);
hSX3FvsB_g = new TH2F("hSX3FvsB_g", "SX3 Front vs Back; Front E; Back E", 800, 0, 16000, 800, 0, 16000);
hsx3IndexVE = new TH2F("hsx3IndexVE", "SX3 index vs Energy; sx3 index ; Energy", 24 * 12, 0, 24 * 12, 400, 0, 5000);
hsx3IndexVE_g = new TH2F("hsx3IndexVE_g", "SX3 index vs Energy; sx3 index ; Energy", 24 * 12, 0, 24 * 12, 400, 0, 5000);
hSX3 = new TH2F("hSX3", "SX3 Front v Back; Fronts; Backs", 8, 0, 8, 4, 0, 4);
hsx3Coin = new TH2F("hsx3Coin", "SX3 Coincident", 24 * 12, 0, 24 * 12, 24 * 12, 0, 24 * 12);
sx3_contr.ConstructGeo();
// Load the TCutG object
TFile *cutFile = TFile::Open("sx3cut.root");
bool cutLoaded = (cut != nullptr);
cut = dynamic_cast<TCutG *>(cutFile->Get("sx3cut"));
if (!cut)
{
std::cerr << "Error: Could not find TCutG named 'sx3cut' in sx3cut.root" << std::endl;
return;
}
cut->SetName("sx3cut"); // Ensure the cut has the correct name
// Load the TCutG object
TFile *cutFile1 = TFile::Open("UvD.root");
bool cut1Loaded = (cut1 != nullptr);
cut1 = dynamic_cast<TCutG *>(cutFile1->Get("UvD"));
if (!cut1)
{
std::cerr << "Error: Could not find TCutG named 'UvD' in UvD.root" << std::endl;
return;
}
cut1->SetName("UvD");
std::string filename = "sx3_BackGains.txt";
std::ifstream infile(filename);
if (!infile.is_open())
{
std::cerr << "Error opening " << filename << "!" << std::endl;
return;
}
int id, bk;
double gain;
while (infile >> id >> bk >> gain)
{
backGain[id][bk] = gain;
if (backGain[id][bk] > 0)
backGainValid[id][bk] = true;
else
backGainValid[id][bk] = false;
}
SX3 sx3_contr;
}
Bool_t GainMatchSX3Front::Process(Long64_t entry)
{
b_sx3Multi->GetEntry(entry);
b_sx3ID->GetEntry(entry);
b_sx3Ch->GetEntry(entry);
b_sx3E->GetEntry(entry);
b_sx3T->GetEntry(entry);
sx3.CalIndex();
std::vector<std::pair<int, int>> ID;
for (int i = 0; i < sx3.multi; i++)
{
for (int j = i + 1; j < sx3.multi; j++)
{
// if (sx3.id[i] == 3)
hsx3Coin->Fill(sx3.index[i], sx3.index[j]);
}
if (sx3.e[i] > 100)
{
ID.push_back(std::pair<int, int>(sx3.id[i], i));
hsx3IndexVE->Fill(sx3.index[i], sx3.e[i]);
}
}
if (ID.size() > 0)
{
std::sort(ID.begin(), ID.end(), [](const std::pair<int, int> &a, const std::pair<int, int> &b)
{ return a.first < b.first; });
// start with the first entry in the sorted array: channels that belong to the same detector are together in sequenmce
std::vector<std::pair<int, int>> sx3ID;
sx3ID.push_back(ID[0]);
bool found = false;
for (size_t i = 1; i < ID.size(); i++)
{ // Check if id of i belongs to the same detector and then add it to the detector ID vector
if (ID[i].first == sx3ID.back().first)
{ // count the nunmber of hits that belong to the same detector
sx3ID.push_back(ID[i]);
if (sx3ID.size() >= 3)
{
found = true;
}
}
else
{ // the next event does not belong to the same detector, abandon the first event and continue with the next one
if (!found)
{
sx3ID.clear();
sx3ID.push_back(ID[i]);
}
}
}
if (found)
{
int sx3ChUp = -1, sx3ChDn = -1, sx3ChBk = -1;
float sx3EUp = 0.0, sx3EDn = 0.0, sx3EBk = 0.0;
for (size_t i = 0; i < sx3ID.size(); i++)
{
int index = sx3ID[i].second;
// Check the channel number and assign it to the appropriate channel type
if (sx3.ch[index] < 8)
{
if (sx3.ch[index] % 2 == 0)
{
sx3ChDn = sx3.ch[index] / 2;
sx3EDn = sx3.e[index];
}
else
{
sx3ChUp = sx3.ch[index] / 2;
sx3EUp = sx3.e[index];
}
}
else
{
sx3ChBk = sx3.ch[index] - 8;
// if (sx3ChBk == 2)
// printf("Found back channel Det %d Back %d \n", sx3.id[index], sx3ChBk);
sx3EBk = sx3.e[index];
}
}
for (int i = 0; i < sx3.multi; i++)
{
// If we have a valid front and back channel, fill the histograms
hSX3->Fill(sx3ChDn + 4, sx3ChBk);
hSX3->Fill(sx3ChUp, sx3ChBk);
// Fill the histogram for the front vs back
hSX3FvsB->Fill(sx3EUp + sx3EDn, sx3EBk);
if (sx3.e[i] > 100 && sx3.id[i] == 3)
{
// back gain correction
// Fill the histogram for the front vs back with gain correction
// hSX3FvsB_g->Fill(sx3EUp + sx3EDn, sx3EBk);
// // Fill the index vs energy histogram
// hsx3IndexVE_g->Fill(sx3.index[i], sx3.e[i]);
// }
// {
TString histName = Form("hSX3FVB_id%d_U%d_D%d_B%d", sx3.id[i], sx3ChUp, sx3ChDn, sx3ChBk);
TH2F *hist2d = (TH2F *)gDirectory->Get(histName);
if (!hist2d)
{
hist2d = new TH2F(histName, Form("hSX3FVB_id%d_U%d_D%d_B%d", sx3.id[i], sx3ChUp, sx3ChDn, sx3ChBk), 400, 0, 16000, 400, 0, 16000);
}
// if (sx3ChBk == 2)
// printf("Found back channel Det %d Back %d \n", sx3.id[i], sx3ChBk);
hsx3IndexVE_g->Fill(sx3.index[i], sx3.e[i]);
hSX3FvsB_g->Fill(sx3EUp + sx3EDn, sx3EBk);
hist2d->Fill(sx3EUp + sx3EDn, sx3EBk);
if (cut && cut->IsInside(sx3EUp + sx3EDn, sx3EBk) && cut1 && cut1->IsInside(sx3EUp / sx3EBk, sx3EDn / sx3EBk))
{
if (backGainValid[sx3.id[i]][sx3ChBk])
{
sx3EBk *= backGain[sx3.id[i]][sx3ChBk];
}
// Accumulate data for gain matching
dataPoints[{sx3.id[i], sx3ChBk, sx3ChUp, sx3ChDn}].emplace_back(sx3EBk, sx3EUp, sx3EDn);
}
}
}
}
}
return kTRUE;
}
void GainMatchSX3Front::Terminate()
{
std::map<std::tuple<int, int, int, int>, TVectorD> fitCoefficients;
// === Gain matching ===
std::ofstream outFile("sx3_GainMatchfront.txt");
if (!outFile.is_open())
{
std::cerr << "Error opening output file!" << std::endl;
return;
}
TH2F *hUvD = new TH2F("hUvD", " UvD; Up/CorrBack; Down/CorrBack", 600, 0, 1, 600, 0, 1);
for (const auto &kv : dataPoints)
{
auto [id, bk, u, d] = kv.first;
const auto &pts = kv.second;
if (pts.size() < 50)
continue;
std::vector<double> uE, dE, udE, corrBkE;
for (const auto &pr : pts)
{
double eBkCorr, eUp, eDn;
std::tie(eBkCorr, eUp, eDn) = pr;
if ((eBkCorr < 100) || (eUp < 100) || (eDn < 100))
continue; // Skip if any energy is zero
uE.push_back(eUp / eBkCorr);
dE.push_back(eDn / eBkCorr);
udE.push_back(eUp + eDn);
corrBkE.push_back(eBkCorr);
hUvD->Fill(eUp / eBkCorr, eDn / eBkCorr);
}
if (uE.size() < 5 || dE.size() < 5 || corrBkE.size() < 5)
continue; // Ensure we have enough points for fitting
// TGraph g(udE.size(), udE.data(), corrBkE.data());
// TF1 f("f", "[0]*x", 0, 20000);
// f.SetParameter(0, 1.0); // Initial guess for the gain
// g.Fit(&f, "R");
const double fixedError = 0.0; // in ADC channels
std::vector<double> xVals, yVals, exVals, eyVals;
// Build data with fixed error
for (size_t i = 0; i < udE.size(); ++i)
{
double x = uE[i]; // front energy
double y = dE[i]; // back energy
xVals.push_back(x);
yVals.push_back(y);
exVals.push_back(fixedError); // error in up energy
eyVals.push_back(0.); // error in down energy
}
// Build TGraphErrors with errors
TGraphErrors g(xVals.size(), xVals.data(), yVals.data(), exVals.data(), eyVals.data());
TF1 f("f", "[0]*x+[1]", 0, 16000);
f.SetParameter(0, -1.0); // Initial guess
if (drawCanvases)
{
g.SetTitle(Form("Detector %d: U%d D%d B%d", id, u, d, bk));
g.SetMarkerStyle(20);
g.SetMarkerColor(kBlue);
g.Draw("AP");
g.Fit(&f, interactiveMode ? "Q" : "QNR"); // 'R' avoids refit, 'N' skips drawing
if (verboseFit)
{
double chi2 = f.GetChisquare();
int ndf = f.GetNDF();
double reducedChi2 = (ndf != 0) ? chi2 / ndf : -1;
std::cout << Form("Det%d U%d D%d B%d → Gain: %.4f | χ²/ndf = %.2f/%d = %.2f",
id, u, d, bk, f.GetParameter(0), chi2, ndf, reducedChi2)
<< std::endl;
}
if (interactiveMode)
{
c.Update();
gPad->WaitPrimitive();
}
else
{
c.Close(); // Optionally avoid clutter in batch
}
}
else
{
g.Fit(&f, "QNR");
}
double slope = f.GetParameter(0);
double intercept = f.GetParameter(1);
// printf("Front gain Det%d Back%d Up%dDn%d → %.4f\n", id, bk, u, d, frontGain[id][bk][u][d]);
if (std::abs(slope + 1.0) < 0.3) // sanity check
{
frontGain[id][bk][u][d] = slope;
frontGainValid[id][bk][u][d] = true;
outFile << id << " " << bk << " " << u << " " << d << " " << TMath::Abs(slope)/intercept << " " << 1.0/intercept << std::endl;
printf("Back slope Det%d Bk%d → %.4f\n", id, bk, slope);
}
else
{
std::cerr << "Warning: Bad slope for Det" << id << " Bk" << bk
<< " slope=" << f.GetParameter(0) << std::endl;
}
}
outFile.close();
std::cout << "Gain matching complete." << std::endl;
// === Stage 3: Create corrected histogram ===
TH2F *hCorrectedFvB = new TH2F("hCorrectedFvB", "Corrected;Corrected Front Sum;Corrected Back", 800, 0, 8000, 800, 0, 8000);
TH2F *hCorrectedUvD = new TH2F("hCorrectedUvD", "Corrected UvD; UvD Up; UvD Down", 600, 0, 1, 600, 0, 1);
for (const auto &kv : dataPoints)
{
auto [id, bk, u, d] = kv.first;
double front;
if (frontGainValid[id][bk][u][d])
front = frontGain[id][bk][u][d];
else
continue;
for (const auto &pr : kv.second)
{
double eBk, eUp, eDn;
std::tie(eBk, eUp, eDn) = pr;
double corrUp = eUp * front;
// double corrDn = eDn * front;
hCorrectedFvB->Fill(corrUp + eDn, eBk);
hCorrectedUvD->Fill(corrUp / eBk, eDn / eBk);
}
}
// // === Final canvas ===
// gStyle->SetOptStat(1110);
// TCanvas *c1 = new TCanvas("c1", "Gain Correction Results", 1200, 600);
// c1->Divide(2, 1);
// c1->cd(1);
// hSX3FvsB_g->SetTitle("Before Correction (Gated)");
// hSX3FvsB_g->GetXaxis()->SetTitle("Measured Front Sum (E_Up + E_Dn)");
// hSX3FvsB_g->GetYaxis()->SetTitle("Measured Back E");
// hSX3FvsB_g->Draw("colz");
// c1->cd(2);
// hCorrectedFvB->SetTitle("After Correction");
// hCorrectedFvB->Draw("colz");
// TF1 *diag = new TF1("diag", "x", 0, 40000);
// diag->SetLineColor(kRed);
// diag->SetLineWidth(2);
// diag->Draw("same");
std::cout << "Terminate() completed successfully." << std::endl;
}

104
GainMatchSX3Front.h
View File

@ -1,104 +0,0 @@
#ifndef GainMatchSX3Front_h
#define GainMatchSX3Front_h
#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
#include <TSelector.h>
#include "Armory/ClassDet.h"
class GainMatchSX3Front : public TSelector {
public :
TTree *fChain; //!pointer to the analyzed TTree or TChain
// Fixed size dimensions of array or collections stored in the TTree if any.
// Declaration of leaf types
Det sx3;
Det qqq;
Det pc ;
ULong64_t evID;
UInt_t run;
// List of branches
TBranch *b_eventID; //!
TBranch *b_run; //!
TBranch *b_sx3Multi; //!
TBranch *b_sx3ID; //!
TBranch *b_sx3Ch; //!
TBranch *b_sx3E; //!
TBranch *b_sx3T; //!
TBranch *b_qqqMulti; //!
TBranch *b_qqqID; //!
TBranch *b_qqqCh; //!
TBranch *b_qqqE; //!
TBranch *b_qqqT; //!
TBranch *b_pcMulti; //!
TBranch *b_pcID; //!
TBranch *b_pcCh; //!
TBranch *b_pcE; //!
TBranch *b_pcT; //!
GainMatchSX3Front(TTree * /*tree*/ =0) : fChain(0) { }
virtual ~GainMatchSX3Front() { }
virtual Int_t Version() const { return 2; }
virtual void Begin(TTree *tree);
virtual void SlaveBegin(TTree *tree);
virtual void Init(TTree *tree);
virtual Bool_t Notify();
virtual Bool_t Process(Long64_t entry);
virtual Int_t GetEntry(Long64_t entry, Int_t getall = 0) { return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0; }
virtual void SetOption(const char *option) { fOption = option; }
virtual void SetObject(TObject *obj) { fObject = obj; }
virtual void SetInputList(TList *input) { fInput = input; }
virtual TList *GetOutputList() const { return fOutput; }
virtual void SlaveTerminate();
virtual void Terminate();
ClassDef(GainMatchSX3Front,0);
};
#endif
#ifdef GainMatchSX3Front_cxx
void GainMatchSX3Front::Init(TTree *tree){
// Set branch addresses and branch pointers
if (!tree) return;
fChain = tree;
fChain->SetMakeClass(1);
fChain->SetBranchAddress("evID", &evID, &b_eventID);
fChain->SetBranchAddress("run", &run, &b_run);
sx3.SetDetDimension(24,12);
qqq.SetDetDimension(4,32);
pc.SetDetDimension(2,24);
fChain->SetBranchAddress("sx3Multi", &sx3.multi, &b_sx3Multi);
fChain->SetBranchAddress("sx3ID", &sx3.id, &b_sx3ID);
fChain->SetBranchAddress("sx3Ch", &sx3.ch, &b_sx3Ch);
fChain->SetBranchAddress("sx3E", &sx3.e, &b_sx3E);
fChain->SetBranchAddress("sx3T", &sx3.t, &b_sx3T);
}
Bool_t GainMatchSX3Front::Notify(){
return kTRUE;
}
void GainMatchSX3Front::SlaveBegin(TTree * /*tree*/){
TString option = GetOption();
}
void GainMatchSX3Front::SlaveTerminate(){
}
#endif // #ifdef GainMatchSX3Front_cxx

245
GainMatchSX3Front1.C
View File

@ -1,245 +0,0 @@
#define GainMatchSX3_cxx
#include "GainMatchSX3.h"
#include "Armory/ClassSX3.h"
#include <TFile.h>
#include <TTree.h>
#include <TGraph.h>
#include <TF1.h>
#include <TH2F.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <TApplication.h>
#include <map>
#include <vector>
#include <tuple>
#include <fstream>
#include <iostream>
#include <algorithm>
// Constants
const int MAX_DET = 24;
const int MAX_BK = 4;
const int MAX_UP = 4;
const int MAX_DOWN = 4;
// Gain arrays
double backGain[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
bool backGainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}};
double frontGain[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
bool frontGainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}};
// Data container
std::map<std::tuple<int, int, int, int>, std::vector<std::tuple<double, double, double>>> dataPoints;
// Load back gains
void LoadBackGains(const std::string &filename)
{
std::ifstream infile(filename);
if (!infile.is_open())
{
std::cerr << "Error opening " << filename << "!" << std::endl;
return;
}
int id, bk, u, d;
double gain;
while (infile >> id >> bk >> u >> d >> gain)
{
backGain[id][bk][u][d] = gain;
backGainValid[id][bk][u][d] = true;
}
infile.close();
std::cout << "Loaded back gains from " << filename << std::endl;
SX3 sx3_contr;
}
// Front gain matching function
Bool_t GainMatchSX3::Process(Long64_t entry)
{
// Link SX3 branches
b_sx3Multi->GetEntry(entry);
b_sx3ID->GetEntry(entry);
b_sx3Ch->GetEntry(entry);
b_sx3E->GetEntry(entry);
b_sx3T->GetEntry(entry);
sx3.CalIndex();
Long64_t nentries = tree->GetEntries(Long64_t entry);
std::cout << "Total entries: " << nentries << std::endl;
TH2F *hBefore = new TH2F("hBefore", "Before Correction;E_Up+E_Dn;Back Energy", 400, 0, 40000, 400, 0, 40000);
TH2F *hAfter = new TH2F("hAfter", "After Correction;E_Up+E_Dn;Corrected Back Energy", 400, 0, 40000, 400, 0, 40000);
for (Long64_t entry = 0; entry < nentries; ++entry)
{
tree->GetEntry(entry);
sx3.CalIndex();
std::vector<std::pair<int, int>> ID;
for (int i = 0; i < sx3.multi; i++)
{
if (sx3.e[i] > 100)
{
ID.push_back({sx3.id[i], i});
}
}
if (ID.empty())
continue;
// Sort by id
std::sort(ID.begin(), ID.end(), [](auto &a, auto &b) { return a.first < b.first; });
std::vector<std::pair<int, int>> sx3ID;
sx3ID.push_back(ID[0]);
bool found = false;
for (size_t i = 1; i < ID.size(); i++)
{
if (ID[i].first == sx3ID.back().first)
{
sx3ID.push_back(ID[i]);
if (sx3ID.size() >= 3)
found = true;
}
else if (!found)
{
sx3ID.clear();
sx3ID.push_back(ID[i]);
}
}
if (!found)
continue;
int sx3ChUp = -1, sx3ChDn = -1, sx3ChBk = -1;
float sx3EUp = 0.0, sx3EDn = 0.0, sx3EBk = 0.0;
int sx3id = sx3ID[0].first;
for (auto &[id, idx] : sx3ID)
{
if (sx3.ch[idx] < 8)
{
if (sx3.ch[idx] % 2 == 0)
{
sx3ChDn = sx3.ch[idx] / 2;
sx3EDn = sx3.e[idx];
}
else
{
sx3ChUp = sx3.ch[idx] / 2;
sx3EUp = sx3.e[idx];
}
}
else
{
sx3ChBk = sx3.ch[idx] - 8;
sx3EBk = sx3.e[idx];
}
}
if (sx3ChUp < 0 || sx3ChDn < 0 || sx3ChBk < 0)
continue;
if (!backGainValid[sx3id][sx3ChBk][sx3ChUp][sx3ChDn])
continue;
double corrBk = sx3EBk * backGain[sx3id][sx3ChBk][sx3ChUp][sx3ChDn];
hBefore->Fill(sx3EUp + sx3EDn, sx3EBk);
hAfter->Fill(sx3EUp + sx3EDn, corrBk);
dataPoints[{sx3id, sx3ChBk, sx3ChUp, sx3ChDn}].emplace_back(corrBk, sx3EUp, sx3EDn);
}
// === Fit front gains ===
std::ofstream outFile("sx3_GainMatchfront.txt");
if (!outFile.is_open())
{
std::cerr << "Error opening sx3_GainMatchfront.txt!" << std::endl;
return;
}
for (const auto &kv : dataPoints)
{
auto [id, bk, u, d] = kv.first;
const auto &pts = kv.second;
if (pts.size() < 5)
continue;
std::vector<double> udE, corrBkE;
for (const auto &pr : pts)
{
double eBkCorr, eUp, eDn;
std::tie(eBkCorr, eUp, eDn) = pr;
udE.push_back(eUp + eDn);
corrBkE.push_back(eBkCorr);
}
TGraph g(udE.size(), udE.data(), corrBkE.data());
TF1 f("f", "[0]*x", 0, 40000);
g.Fit(&f, "QNR");
frontGain[id][bk][u][d] = f.GetParameter(0);
frontGainValid[id][bk][u][d] = true;
outFile << id << " " << bk << " " << u << " " << d << " " << frontGain[id][bk][u][d] << std::endl;
printf("Front gain Det%d Back%d Up%dDn%d → %.4f\n", id, bk, u, d, frontGain[id][bk][u][d]);
}
outFile.close();
std::cout << "Front gain matching complete." << std::endl;
// === Draw diagnostic plots ===
gStyle->SetOptStat(1110);
TCanvas *c = new TCanvas("c", "Gain Matching Diagnostics", 1200, 600);
c->Divide(2, 1);
c->cd(1);
hBefore->Draw("colz");
TF1 *diag1 = new TF1("diag1", "x", 0, 40000);
diag1->SetLineColor(kRed);
diag1->Draw("same");
c->cd(2);
hAfter->Draw("colz");
TF1 *diag2 = new TF1("diag2", "x", 0, 40000);
diag2->SetLineColor(kRed);
diag2->Draw("same");
}
int main(int argc, char **argv)
{
TApplication app("app", &argc, argv);
// Load back gains
LoadBackGains("sx3_GainMatchback.txt");
// Open tree
TFile *f = TFile::Open("input_tree.root"); // <<< Change file name
if (!f || f->IsZombie())
{
std::cerr << "Cannot open input_tree.root!" << std::endl;
return 1;
}
TTree *tree = (TTree *)f->Get("tree");
if (!tree)
{
std::cerr << "Tree not found!" << std::endl;
return 1;
}
// Run front gain matching
GainMatchSX3(tree);
app.Run();
return 0;
}

168
MakeVertex.C
View File

@ -1,5 +1,15 @@
#define MakeVertex_cxx
// Comment out any line below to disable that diagnostic section
// #define DIAG_WIREMULT // anode/cathode cluster multiplicity plots
#define DIAG_1WIRE // raw per-wire dE vs Si (no PC required)
#define DIAG_PC_SX3 // PC-SX3 coincidence analysis
// #define DIAG_1A0C_SX3 // 1A0C single-wire vertex with SX3
// #define DIAG_1A0C_QQQ // 1A0C single-wire vertex with QQQ
#define DIAG_NA0C_SX3 // nA0C (n>=2) pseudo-wire vertex with SX3
#define DIAG_NA0C_QQQ // nA0C (n>=2) pseudo-wire vertex with QQQ
#define DIAG_PC_QQQ // PC-QQQ coincidence analysis
Int_t colors[40] = {
kBlack, kRed, kGreen, kBlue, kYellow, kMagenta, kCyan, kOrange,
kSpring, kTeal, kAzure, kViolet, kPink, kGray, kWhite,
@ -851,6 +861,7 @@ Bool_t MakeVertex::Process(Long64_t entry)
if (QQQ_Events.size() && PC_Events.size())
plotter->Fill2D("PCEv_vs_QQQEv", 20, 0, 20, 20, 0, 20, QQQ_Events.size(), PC_Events.size());
#ifdef DIAG_WIREMULT
plotter->Fill2D("ac_vs_cc", 20, 0, 20, 20, 0, 20, aClusters.size(), cClusters.size(), "wiremult");
for (auto cluster : aClusters)
{
@ -865,7 +876,9 @@ Bool_t MakeVertex::Process(Long64_t entry)
{
plotter->Fill2D("ac_vs_cc_ign0", 20, 0, 20, 20, 0, 20, aClusters.size(), cClusters.size(), "wiremult");
}
#endif // DIAG_WIREMULT
#ifdef DIAG_1WIRE
for (auto sx3event : SX3_Events)
{
for (int i = 0; i < 24; i++)
@ -940,7 +953,9 @@ Bool_t MakeVertex::Process(Long64_t entry)
}
} // for 'i' loop
}
#endif // DIAG_1WIRE
#ifdef DIAG_PC_SX3
bool PCSX3PhiCut = false;
for (auto pcevent : PC_Events)
{
@ -1225,6 +1240,7 @@ Bool_t MakeVertex::Process(Long64_t entry)
}*/
}
} // end PC-SX3 coincidence
#endif // DIAG_PC_SX3
/*for(size_t ii=0; ii<QQQ_Events.size(); ii++) {
for(size_t jj=ii+1; jj<QQQ_Events.size(); jj++) {
@ -1260,6 +1276,7 @@ Bool_t MakeVertex::Process(Long64_t entry)
///////////////////Single wire analysis for the anodes///////////////////
#if defined(DIAG_1A0C_SX3) || defined(DIAG_1A0C_QQQ)
if (aClusters.size() == 1 && cClusters.size() == 0)
{
// Extract the primary anode hit properties
@ -1273,6 +1290,7 @@ Bool_t MakeVertex::Process(Long64_t entry)
TVector3 wireVec = pwinstance.An[aWireID].first - pwinstance.An[aWireID].second;
// Loop over SX3_Events directly
#ifdef DIAG_1A0C_SX3
for (auto sx3event : SX3_Events)
{
@ -1369,9 +1387,10 @@ Bool_t MakeVertex::Process(Long64_t entry)
}
}
}
#endif // DIAG_1A0C_SX3
// Loop over QQQ_Events directly
#ifdef DIAG_1A0C_QQQ
for (auto qqqevent : QQQ_Events)
{
if (qqqevent.Time1 - aTime < -150) // Time cut for protons
@ -1468,8 +1487,154 @@ Bool_t MakeVertex::Process(Long64_t entry)
}
}
}
#endif // DIAG_1A0C_QQQ
}
#endif // DIAG_1A0C_SX3 || DIAG_1A0C_QQQ
///////////////////nA0C analysis using pseudo-wire (GetPseudoWire + getClosestWirePosAtWirePhi)///////////////////
#if defined(DIAG_NA0C_SX3) || defined(DIAG_NA0C_QQQ)
if (cClusters.size() == 0 && aClusters.size() >= 2)
{
std::string nA0C_label = std::to_string(aClusters.size()) + "A0C";
// Flatten all anode clusters into a combined hit list for the pseudo-wire
std::vector<std::tuple<int, double, double>> allAnodeHits;
for (const auto &acluster : aClusters)
for (const auto &hit : acluster)
allAnodeHits.push_back(hit);
auto [apwire, apSumE, apMaxE, apTSMaxE] = pwinstance.GetPseudoWire(allAnodeHits, "ANODE");
#ifdef DIAG_NA0C_SX3
for (auto sx3event : SX3_Events)
{
if (sx3event.Time1 - apTSMaxE < -150)
{
TVector3 pcz_intersect = pwinstance.getClosestWirePosAtWirePhi(apwire, sx3event.pos.Phi());
double deltaRho = sx3event.pos.Perp() - pcz_intersect.Perp();
double deltaZ = sx3event.pos.Z() - pcz_intersect.Z();
double vertex_recon = sx3event.pos.Z() - sx3event.pos.Perp() * (deltaZ / deltaRho);
double z_entrance = -274.3;
double beam_path_length = TMath::Abs(vertex_recon - z_entrance) * 0.1;
double initial_beam_energy = 72.0;
double beam_energy_at_vertex = cm_to_MeV_beam->Eval(
MeV_to_cm_beam->Eval(initial_beam_energy) - beam_path_length);
Kinematics apkin_p(26.981538, 4.002603, 1.007825, 29.973770, beam_energy_at_vertex);
Kinematics apkin_a(26.981538, 4.002603, 4.002603, 26.981538, beam_energy_at_vertex);
std::string vtx_gate = "";
if (vertex_recon >= -176.0 && vertex_recon < -100.0)
vtx_gate = "_Z[-176_to_-100]";
else if (vertex_recon >= -100.0 && vertex_recon < -50.0)
vtx_gate = "_Z[-100_to_-50]";
else if (vertex_recon >= -50.0 && vertex_recon < 0.0)
vtx_gate = "_Z[-50_to_0]";
else if (vertex_recon >= 0.0 && vertex_recon < 50.0)
vtx_gate = "_Z[0_to_50]";
else if (vertex_recon >= 50.0 && vertex_recon < 100.0)
vtx_gate = "_Z[50_to_100]";
else if (vertex_recon >= 100.0 && vertex_recon < 176.0)
vtx_gate = "_Z[100_to_176]";
double path_length = (sx3event.pos - TVector3(0, 0, vertex_recon)).Mag() * 0.1;
double sx3Efix = cm_to_MeVp->Eval(MeV_to_cm_p->Eval(sx3event.Energy1) - path_length);
double sx3Efixalpha = cm_to_MeV->Eval(MeV_to_cm->Eval(sx3event.Energy1) - path_length);
double theta_recon = (sx3event.pos - TVector3(0, 0, vertex_recon)).Theta();
double sinTheta = TMath::Sin(theta_recon);
double Ex_from_proton = apkin_p.getExc(sx3Efix, theta_recon * 180. / M_PI);
double Ex_from_alpha = apkin_a.getExc(sx3Efixalpha, theta_recon * 180. / M_PI);
plotter->Fill2D(nA0C_label + "_dE_Ecorr_Anode_SX3", 400, 0, 30, 800, 0, 40000, sx3Efix, apSumE * sinTheta, nA0C_label);
plotter->Fill1D(nA0C_label + "_Ex_from_alphas_SX3" + vtx_gate, 200, -10, 10, Ex_from_alpha, nA0C_label);
plotter->Fill1D(nA0C_label + "_Ex_from_protons_SX3" + vtx_gate, 200, -10, 10, Ex_from_proton, nA0C_label);
plotter->Fill2D(nA0C_label + "_sx3_E_vs_theta_raw_SX3", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, sx3event.Energy1, nA0C_label);
plotter->Fill2D(nA0C_label + "_sx3_E_vs_theta_corr_SX3", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, sx3Efix, nA0C_label);
if (vtx_gate != "")
{
plotter->Fill1D(nA0C_label + "_twisted_pcz_recon_SX3" + vtx_gate, 600, -300, 300, pcz_intersect.Z(), nA0C_label);
plotter->Fill1D(nA0C_label + "_twisted_vertex_recon_SX3" + vtx_gate, 600, -300, 300, vertex_recon, nA0C_label);
plotter->Fill2D(nA0C_label + "_dE_Ecorr_Anode_SX3" + vtx_gate, 400, 0, 30, 800, 0, 40000, sx3Efix, apSumE * sinTheta, nA0C_label);
plotter->Fill2D(nA0C_label + "_dE_Ecorr_Anode_SX3_alpha" + vtx_gate, 400, 0, 30, 800, 0, 40000, sx3Efixalpha, apSumE * sinTheta, nA0C_label);
plotter->Fill1D(nA0C_label + "_Ex_from_alphas_SX3" + vtx_gate, 200, -10, 10, Ex_from_alpha, nA0C_label);
plotter->Fill1D(nA0C_label + "_Ex_from_protons_SX3" + vtx_gate, 200, -10, 10, Ex_from_proton, nA0C_label);
}
}
}
#endif // DIAG_NA0C_SX3
#ifdef DIAG_NA0C_QQQ
for (auto qqqevent : QQQ_Events)
{
if (qqqevent.Time1 - apTSMaxE < -150)
{
TVector3 pcz_intersect = pwinstance.getClosestWirePosAtWirePhi(apwire, qqqevent.pos.Phi());
double deltaRho = qqqevent.pos.Perp() - pcz_intersect.Perp();
double deltaZ = qqqevent.pos.Z() - pcz_intersect.Z();
double vertex_recon = qqqevent.pos.Z() - qqqevent.pos.Perp() * (deltaZ / deltaRho);
double z_entrance = -274.3;
double beam_path_length = TMath::Abs(vertex_recon - z_entrance) * 0.1;
double initial_beam_energy = 72.0;
double beam_energy_at_vertex = cm_to_MeV_beam->Eval(
MeV_to_cm_beam->Eval(initial_beam_energy) - beam_path_length);
Kinematics apkin_p(26.981538, 4.002603, 1.007825, 29.973770, beam_energy_at_vertex);
Kinematics apkin_a(26.981538, 4.002603, 4.002603, 26.981538, beam_energy_at_vertex);
std::string vtx_gate = "";
if (vertex_recon >= -176.0 && vertex_recon < -100.0)
vtx_gate = "_Z[-176_to_-100]";
else if (vertex_recon >= -100.0 && vertex_recon < -50.0)
vtx_gate = "_Z[-100_to_-50]";
else if (vertex_recon >= -50.0 && vertex_recon < 0.0)
vtx_gate = "_Z[-50_to_0]";
else if (vertex_recon >= 0.0 && vertex_recon < 50.0)
vtx_gate = "_Z[0_to_50]";
else if (vertex_recon >= 50.0 && vertex_recon < 100.0)
vtx_gate = "_Z[50_to_100]";
else if (vertex_recon >= 100.0 && vertex_recon < 176.0)
vtx_gate = "_Z[100_to_176]";
double path_length = (qqqevent.pos - TVector3(0, 0, vertex_recon)).Mag() * 0.1;
double qqqEfix = cm_to_MeVp->Eval(MeV_to_cm_p->Eval(qqqevent.Energy1) - path_length);
double qqqEfixalpha = cm_to_MeV->Eval(MeV_to_cm->Eval(qqqevent.Energy2) - path_length);
double theta_recon = (qqqevent.pos - TVector3(0, 0, vertex_recon)).Theta();
double sinTheta = TMath::Sin(theta_recon);
double Ex_from_proton = apkin_p.getExc(qqqEfix, theta_recon * 180. / M_PI);
double Ex_from_alpha = apkin_a.getExc(qqqEfixalpha, theta_recon * 180. / M_PI);
plotter->Fill2D(nA0C_label + "_dE_Ecorr_Anode_QQQ", 400, 0, 30, 800, 0, 40000, qqqEfix, apSumE * sinTheta, nA0C_label);
plotter->Fill1D(nA0C_label + "_Ex_from_alphas_QQQ" + vtx_gate, 200, -10, 10, Ex_from_alpha, nA0C_label);
plotter->Fill1D(nA0C_label + "_Ex_from_protons_QQQ" + vtx_gate, 200, -10, 10, Ex_from_proton, nA0C_label);
plotter->Fill2D(nA0C_label + "_qqq_E_vs_theta_raw_QQQ", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, qqqevent.Energy1, nA0C_label);
plotter->Fill2D(nA0C_label + "_qqq_E_vs_theta_corr_QQQ", 180, 0, 180, 400, 0, 30, theta_recon * 180. / M_PI, qqqEfix, nA0C_label);
if (vtx_gate != "")
{
plotter->Fill1D(nA0C_label + "_twisted_pcz_recon_QQQ" + vtx_gate, 600, -300, 300, pcz_intersect.Z(), nA0C_label);
plotter->Fill1D(nA0C_label + "_twisted_vertex_recon_QQQ" + vtx_gate, 600, -300, 300, vertex_recon, nA0C_label);
plotter->Fill2D(nA0C_label + "_dE_Ecorr_Anode_QQQ" + vtx_gate, 400, 0, 30, 800, 0, 40000, qqqEfix, apSumE * sinTheta, nA0C_label);
plotter->Fill2D(nA0C_label + "_dE_Ecorr_Anode_QQQ_alpha" + vtx_gate, 400, 0, 30, 800, 0, 40000, qqqEfixalpha, apSumE * sinTheta, nA0C_label);
plotter->Fill1D(nA0C_label + "_Ex_from_alphas_QQQ" + vtx_gate, 200, -10, 10, Ex_from_alpha, nA0C_label);
plotter->Fill1D(nA0C_label + "_Ex_from_protons_QQQ" + vtx_gate, 200, -10, 10, Ex_from_proton, nA0C_label);
}
}
}
#endif // DIAG_NA0C_QQQ
}
#endif // DIAG_NA0C_SX3 || DIAG_NA0C_QQQ
#ifdef DIAG_PC_QQQ
for (auto pcevent : PC_Events)
{
for (auto qqqevent : QQQ_Events)
@ -1664,6 +1829,7 @@ Bool_t MakeVertex::Process(Long64_t entry)
}
}
} // end PC QQQ coincidence
#endif // DIAG_PC_QQQ
// HALFTIME! Can stop here in future versions
// return kTRUE;
if (anodeHits.size() >= 1 && cathodeHits.size() >= 1)

2
run_27Al.sh
View File

@ -2,7 +2,7 @@ rm results_run*.root
export DATASET="27Al"
export flip180="0"
export flipa=0
export anode_offset=1
export anode_offset=2
#declare -i run=28
#while [[ $run -lt 34 ]]; do #runs 1 to 84
# wrun=$(printf "%03d" $run)

12
sx3_BackGains.txt
View File

@ -1,12 +0,0 @@
7 0 1.20298
7 1 0.995493
7 2 0.993613
7 3 1.2514
9 0 1.01574
9 1 0.961032
9 2 0.988379
9 3 1.05832
19 0 1.07936
19 1 0.97626
19 2 1.00078
19 3 1.03335

8
sx3_GainMatchfront.txt
View File

@ -1,8 +0,0 @@
3 0 1 1 0.852399 0.881228
3 0 2 2 0.813845 0.975967
3 0 3 3 0.859643 0.863715
3 1 1 1 0.76728 0.942438
3 1 2 2 0.780302 0.929008
3 2 2 2 0.729082 1.02005
3 3 2 2 0.759098 1.05376
3 3 3 3 0.821183 0.952335