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