From c142f4b567f35e6db9227b9595eac8b4920fa2ab Mon Sep 17 00:00:00 2001 From: vsitaraman Date: Fri, 15 May 2026 14:02:00 -0400 Subject: [PATCH] deleted: Analyzer1.C 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. --- Analyzer1.C | 402 -------------------------------------- Analyzer1.h | 114 ----------- GainMatchSX3.C | 432 ----------------------------------------- GainMatchSX3.h | 114 ----------- GainMatchSX3Front.C | 420 --------------------------------------- GainMatchSX3Front.h | 104 ---------- GainMatchSX3Front1.C | 245 ----------------------- MakeVertex.C | 170 +++++++++++++++- run_27Al.sh | 2 +- sx3_BackGains.txt | 12 -- sx3_GainMatchfront.txt | 8 - 11 files changed, 169 insertions(+), 1854 deletions(-) delete mode 100644 Analyzer1.C delete mode 100644 Analyzer1.h delete mode 100644 GainMatchSX3.C delete mode 100644 GainMatchSX3.h delete mode 100644 GainMatchSX3Front.C delete mode 100644 GainMatchSX3Front.h delete mode 100644 GainMatchSX3Front1.C delete mode 100644 sx3_BackGains.txt delete mode 100644 sx3_GainMatchfront.txt diff --git a/Analyzer1.C b/Analyzer1.C deleted file mode 100644 index 62c47d1..0000000 --- a/Analyzer1.C +++ /dev/null @@ -1,402 +0,0 @@ -#define Analyzer1_cxx - -#include "Analyzer1.h" -#include -#include -#include -#include - -#include -#include - -#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> ID; // first = id, 2nd = index - for( int i = 0; i < sx3.multi; i ++){ - ID.push_back(std::pair(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 & a, const std::pair & 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> 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> E; - E.clear(); - for( int i = 0; i < pc.multi; i ++){ - - if( pc.e[i] > 100 ) ID.push_back(std::pair(pc.id[i], i)); - if( pc.e[i] > 100 ) E.push_back(std::pair(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=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"); - -} diff --git a/Analyzer1.h b/Analyzer1.h deleted file mode 100644 index 9d9d3a7..0000000 --- a/Analyzer1.h +++ /dev/null @@ -1,114 +0,0 @@ -#ifndef Analyzer1_h -#define Analyzer1_h - -#include -#include -#include -#include - -#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 diff --git a/GainMatchSX3.C b/GainMatchSX3.C deleted file mode 100644 index dd385af..0000000 --- a/GainMatchSX3.C +++ /dev/null @@ -1,432 +0,0 @@ -#define GainMatchSX3_cxx - -#include "GainMatchSX3.h" -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include "Armory/ClassSX3.h" -#include "Armory/HistPlotter.h" -#include -#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::vector>> dataPoints; -std::map, 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(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(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> 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(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 &a, const std::pair &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> 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 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 exVals(udE.size(), 0.0); // no x error - std::vector 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(); -} diff --git a/GainMatchSX3.h b/GainMatchSX3.h deleted file mode 100644 index 5088a32..0000000 --- a/GainMatchSX3.h +++ /dev/null @@ -1,114 +0,0 @@ -#ifndef GainMatchSX3_h -#define GainMatchSX3_h - -#include -#include -#include -#include - -#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 \ No newline at end of file diff --git a/GainMatchSX3Front.C b/GainMatchSX3Front.C deleted file mode 100644 index 9166c2c..0000000 --- a/GainMatchSX3Front.C +++ /dev/null @@ -1,420 +0,0 @@ -#define GainMatchSX3Front_cxx - -#include "GainMatchSX3Front.h" -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#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::vector>> 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(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(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> 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(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 &a, const std::pair &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> 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, 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 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 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; -} \ No newline at end of file diff --git a/GainMatchSX3Front.h b/GainMatchSX3Front.h deleted file mode 100644 index 47aef61..0000000 --- a/GainMatchSX3Front.h +++ /dev/null @@ -1,104 +0,0 @@ -#ifndef GainMatchSX3Front_h -#define GainMatchSX3Front_h - -#include -#include -#include -#include - -#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 \ No newline at end of file diff --git a/GainMatchSX3Front1.C b/GainMatchSX3Front1.C deleted file mode 100644 index 28ac0e2..0000000 --- a/GainMatchSX3Front1.C +++ /dev/null @@ -1,245 +0,0 @@ -#define GainMatchSX3_cxx - -#include "GainMatchSX3.h" -#include "Armory/ClassSX3.h" -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -// 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::vector>> 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> 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> 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 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; -} diff --git a/MakeVertex.C b/MakeVertex.C index 7a0d97a..d354107 100755 --- a/MakeVertex.C +++ b/MakeVertex.C @@ -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= 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> 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) @@ -2209,4 +2375,4 @@ void protonAlphaHistograms(HistPlotter *plotter, std::vector QQQ_Events, } // end QQQ_Events for loop, end sidetrack a(p,p) return; -} +} \ No newline at end of file diff --git a/run_27Al.sh b/run_27Al.sh index eefded6..bacd00e 100644 --- a/run_27Al.sh +++ b/run_27Al.sh @@ -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) diff --git a/sx3_BackGains.txt b/sx3_BackGains.txt deleted file mode 100644 index b4fed9d..0000000 --- a/sx3_BackGains.txt +++ /dev/null @@ -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 diff --git a/sx3_GainMatchfront.txt b/sx3_GainMatchfront.txt deleted file mode 100644 index a2fd28e..0000000 --- a/sx3_GainMatchfront.txt +++ /dev/null @@ -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