diff --git a/.vscode/settings.json b/.vscode/settings.json index bbc9126..108b698 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -113,7 +113,9 @@ "charconv": "cpp", "format": "cpp", "GainMatchSX3Front.C": "cpp", - "GainMatchSX3Front1.C": "cpp" + "GainMatchSX3Front1.C": "cpp", + "Calibration.C": "cpp", + "GainMatchQQQ.C": "cpp" }, "github-enterprise.uri": "https://fsunuc.physics.fsu.edu" } \ No newline at end of file diff --git a/Calibration.C b/Calibration.C new file mode 100644 index 0000000..46a8906 --- /dev/null +++ b/Calibration.C @@ -0,0 +1,445 @@ +#define Calibration_cxx + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "Armory/ClassSX3.h" +#include "Armory/ClassPW.h" +#include "TGraphErrors.h" +#include "Calibration.h" + +int padID = 0; + +SX3 sx3_contr; +PW pw_contr; +PW pwinstance; +TVector3 hitPos; +// TVector3 anodeIntersection; +std::map> slopeInterceptMap; + +bool HitNonZero; +bool sx3ecut; +bool qqqEcut; + +TH2F *hSX3FvsB; +TH2F *hSX3FvsB_g; +TH2F *hSX3; +TH1F *hZProj; +TH2F *hsx3IndexVE; +TH2F *hqqqIndexVE; +TH2F *hqqqIndexVE_cal; +TH2F *hsx3Coin; +TH2F *hqqqCoin; +TH2F *hqqqPolar; + +TCutG *cut; +TCutG *cut1; + +// Gain arrays + +const int MAX_SX3 = 24; +const int MAX_UP = 4; +const int MAX_DOWN = 4; +const int MAX_BK = 4; +const int MAX_QQQ = 4; +const int MAX_RING = 16; +const int MAX_WEDGE = 16; +double backGain[MAX_SX3][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}}; +bool backGainValid[MAX_SX3][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}}; +double frontGain[MAX_SX3][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}}; +bool frontGainValid[MAX_SX3][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}}; +double uvdslope[MAX_SX3][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}}; +double qqqGain[MAX_QQQ][MAX_BK][MAX_UP] = {{{0}}}; +bool qqqGainValid[MAX_QQQ][MAX_BK][MAX_UP] = {{{false}}}; + +void Calibration::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); + 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); + hsx3IndexVE = new TH2F("hsx3IndexVE", "SX3 index vs Energy; sx3 index ; Energy", 24 * 12, 0, 24 * 12, 400, 0, 5000); + hqqqIndexVE = new TH2F("hqqqIndexVE", "QQQ index vs Energy; QQQ index ; Energy", 4 * 2 * 16, 0, 4 * 2 * 16, 400, 0, 5000); + hqqqIndexVE_cal = new TH2F("hqqqIndexVE_cal", "QQQ index vs Energy (calibrated); QQQ index ; Energy", 4 * 2 * 16, 0, 4 * 2 * 16, 400, 0, 5000); + 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); + + hqqqPolar = new TH2F("hqqqPolar", "QQQ Polar ID", 16 * 4, -TMath::Pi(), TMath::Pi(), 16, 10, 50); + + sx3_contr.ConstructGeo(); + pw_contr.ConstructGeo(); + // ----------------------- Load Back Gains + { + std::string filename = "sx3_GainMatchback.txt"; + std::ifstream infile(filename); + if (!infile.is_open()) + { + std::cerr << "Error opening " << filename << "!" << std::endl; + } + else + { + 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] = (gain > 0); + } + infile.close(); + std::cout << "Loaded back gains from " << filename << std::endl; + } + } + + // ----------------------- Load Front Gains + { + std::string filename = "sx3_GainMatchfront.txt"; + std::ifstream infile(filename); + if (!infile.is_open()) + { + std::cerr << "Error opening " << filename << "!" << std::endl; + } + else + { + int id, bk, u, d; + double gain; + while (infile >> id >> bk >> u >> d >> gain) + { + frontGain[id][bk][u][d] = gain; + frontGainValid[id][bk][u][d] = (gain > 0); + } + infile.close(); + std::cout << "Loaded front gains from " << filename << std::endl; + } + } + + // ----------------------- Load QQQ Gains + { + std::string filename = "qqq_GainMatch.txt"; + std::ifstream infile(filename); + if (!infile.is_open()) + { + std::cerr << "Error opening " << filename << "!" << std::endl; + } + else + { + int det, ring, wedge; + double gain; + while (infile >> det >> ring >> wedge >> gain) + { + qqqGain[det][ring][wedge] = gain; + qqqGainValid[det][ring][wedge] = (gain > 0); + } + infile.close(); + std::cout << "Loaded QQQ gains from " << filename << std::endl; + } + } + + SX3 sx3_contr; +} +Bool_t Calibration::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); + + sx3.CalIndex(); + qqq.CalIndex(); + pc.CalIndex(); + + // sx3.Print(); + + // ########################################################### Raw data + // //======================= SX3 + sx3ecut = false; + 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]); + + if (sx3.e[i] > 100) + { + sx3ecut = true; + } + + for (int j = i + 1; j < sx3.multi; j++) + { + hsx3Coin->Fill(sx3.index[i], sx3.index[j]); + } + } + + // --- safe SX3 handling (replace your existing block that builds sx3ID) --- + if (ID.size() > 0) + { + std::sort(ID.begin(), ID.end(), [](const std::pair &a, const std::pair &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]); + found = false; + } + } + } + + if (found) + { + // initialize to sentinel values + int sx3ChUp = -1; + int sx3ChDn = -1; + int sx3ChBk = -1; + float sx3EUp = 0.0f; + float sx3EDn = 0.0f; + + // collect channels/energies + for (size_t i = 0; i < sx3ID.size(); i++) + { + int index = sx3ID[i].second; + int ch = sx3.ch[index]; + float e = sx3.e[index]; + + if (ch < 8) // front channels + { + // you used even/odd to denote down/up — keep that convention + if ((ch % 2) == 0) // down + { + sx3ChDn = ch; + sx3EDn = e; + } + else // up + { + sx3ChUp = ch; + sx3EUp = e; + } + } + else // back channels (assuming back channels are 8..11 or so) + { + sx3ChBk = ch; // store as raw channel number; adapt if you index bk differently + } + } + + // Basic sanity checks before using indices: + bool haveFrontPair = (sx3ChUp >= 0 && sx3ChDn >= 0); + bool haveBack = (sx3ChBk >= 0); + + // convert raw channel numbers to array indices if needed: + // example: if back channels are stored with ch>=8 and bk index is ch-8: + int bk_index = (haveBack ? (sx3ChBk - 8) : -1); + int up_index = (haveFrontPair ? sx3ChUp : -1); + int dn_index = (haveFrontPair ? sx3ChDn : -1); + + // bounds checks against your array sizes + auto sx3Id = sx3ID[0].first; + auto inRange = [&](int id, int bk, int up, int dn) + { + if (id < 0 || id >= MAX_SX3) + return false; + if (bk < 0 || bk >= MAX_BK) + return false; + if (up < 0 || up >= MAX_UP) + return false; + if (dn < 0 || dn >= MAX_DOWN) + return false; + return true; + }; + + // apply calibration safely, fallback to raw ADC if no valid gain + double calibEUp = sx3EUp; + double calibEDn = sx3EDn; + double calibEBack = 0.0; + + if (haveFrontPair && haveBack && inRange(sx3Id, bk_index, up_index, dn_index)) + { + // If you stored front gains indexed by [id][bk][up][down] + if (frontGainValid[sx3Id][bk_index][up_index][dn_index]) + { + calibEUp = frontGain[sx3Id][bk_index][up_index][dn_index] * sx3EUp; + // calibEDn = frontGain[sx3Id][bk_index][up_index][dn_index] * sx3EDn; + } + if (backGainValid[sx3Id][bk_index][up_index][dn_index]) + { + calibEBack = backGain[sx3Id][bk_index][up_index][dn_index] * sx3.e[sx3ID[0].second /* if you have back index */]; + // Note: if back is from different index in sx3ID vector, you should find the exact index for the back hit + } + } + else + { + // partial information: try best-effort per-channel checks + if (haveFrontPair && up_index >= 0 && dn_index >= 0 && sx3Id >= 0 && sx3Id < MAX_SX3) + { + if (frontGainValid[sx3Id][0][(up_index % MAX_UP)][(dn_index % MAX_DOWN)]) + { + // attempt with default bk=0 if that makes sense in your geometry + calibEUp = frontGain[sx3Id][0][(up_index % MAX_UP)][(dn_index % MAX_DOWN)] * sx3EUp; + // calibEDn = frontGain[sx3Id][0][(up_index % MAX_UP)][(dn_index % MAX_DOWN)] * sx3EDn; + } + } + // keep calibEBack==0 if unavailable + } + + // Only call CalSX3Pos if we have reasonable energies (avoid calling with zeros/uninitialized) + if (haveFrontPair && (calibEUp > 0.0) && haveBack) + { + // find exact back energy value from sx3 entries if you tracked it above + float backEnergyRaw = 0.0f; + // locate the back index in sx3ID if needed + for (size_t k = 0; k < sx3ID.size(); ++k) + { + int idx = sx3ID[k].second; + if (sx3.ch[idx] >= 8) + { + backEnergyRaw = sx3.e[idx]; + break; + } + } + // use calibrated back if available else raw + double backEnergyToUse = (calibEBack > 0.0 ? calibEBack : backEnergyRaw); + hSX3->Fill(sx3ChDn + 4, sx3ChBk); + hSX3->Fill(sx3ChUp, sx3ChBk); + + // Fill the histogram for the front vs back + hSX3FvsB->Fill(sx3EUp + sx3EDn, calibEBack); + + sx3_contr.CalSX3Pos(sx3Id, sx3ChUp, sx3ChDn, sx3ChBk, static_cast(calibEUp), static_cast(calibEDn)); + hitPos = sx3_contr.GetHitPos(); + HitNonZero = true; + } + } // found + } + + // //======================= QQQ + for (int i = 0; i < qqq.multi; i++) + { + + int det = qqq.id[i]; // detector ID (0–3) + int ch = qqq.ch[i]; // raw channel (0–31) + + // Separate ring vs wedge channel + int ring = -1; + int wedge = -1; + if (ch < 16) + { // wedge + wedge = ch; + } + else + { // ring + ring = ch - 16; + } + + double Ecal = qqq.e[i]; // default = raw + if (ring >= 0 && wedge >= 0 && qqqGainValid[det][ring][wedge]) + { + Ecal *= qqqGain[det][ring][wedge]; + } + // for( int j = 0; j < pc.multi; j++){ + // if(pc.index[j]==4){ + hqqqIndexVE_cal->Fill(qqq.index[i], Ecal); + hqqqIndexVE->Fill(qqq.index[i], qqq.e[i]); + + // } + // printf("QQQ ID : %d, ch : %d, e : %d \n", qqq.id[i], qqq.ch[i], qqq.e[i]); + if (qqq.e[i] > 100) + { + qqqEcut = true; + } + // } + 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++) + { + // 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 = 50. + 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; + } + } + } + } + + return kTRUE; +} +void Calibration::Terminate() +{ +} \ No newline at end of file diff --git a/GainMatch.h b/Calibration.h similarity index 88% rename from GainMatch.h rename to Calibration.h index 2611669..202d1dd 100644 --- a/GainMatch.h +++ b/Calibration.h @@ -1,5 +1,5 @@ -#ifndef GainMatch_h -#define GainMatch_h +#ifndef Calibration_h +#define Calibration_h #include #include @@ -8,7 +8,7 @@ #include "Armory/ClassDet.h" -class GainMatch : public TSelector { +class Calibration : public TSelector { public : TTree *fChain; //!pointer to the analyzed TTree or TChain @@ -41,8 +41,8 @@ public : TBranch *b_pcE; //! TBranch *b_pcT; //! - GainMatch(TTree * /*tree*/ =0) : fChain(0) { } - virtual ~GainMatch() { } + Calibration(TTree * /*tree*/ =0) : fChain(0) { } + virtual ~Calibration() { } virtual Int_t Version() const { return 2; } virtual void Begin(TTree *tree); virtual void SlaveBegin(TTree *tree); @@ -57,13 +57,13 @@ public : virtual void SlaveTerminate(); virtual void Terminate(); - ClassDef(GainMatch,0); + ClassDef(Calibration,0); }; #endif -#ifdef GainMatch_cxx -void GainMatch::Init(TTree *tree){ +#ifdef Calibration_cxx +void Calibration::Init(TTree *tree){ // Set branch addresses and branch pointers if (!tree) return; @@ -95,20 +95,20 @@ void GainMatch::Init(TTree *tree){ } -Bool_t GainMatch::Notify(){ +Bool_t Calibration::Notify(){ return kTRUE; } -void GainMatch::SlaveBegin(TTree * /*tree*/){ +void Calibration::SlaveBegin(TTree * /*tree*/){ TString option = GetOption(); } -void GainMatch::SlaveTerminate(){ +void Calibration::SlaveTerminate(){ } -#endif // #ifdef GainMatch_cxx \ No newline at end of file +#endif // #ifdef Calibration_cxx \ No newline at end of file diff --git a/GainMatch.C b/GainMatchQQQ.C similarity index 97% rename from GainMatch.C rename to GainMatchQQQ.C index a5df348..814b25c 100644 --- a/GainMatch.C +++ b/GainMatchQQQ.C @@ -1,6 +1,6 @@ -#define GainMatch_cxx +#define GainMatchQQQ_cxx -#include "GainMatch.h" +#include "GainMatchQQQ.h" #include #include #include @@ -24,7 +24,7 @@ SX3 sx3_contr; TCutG *cut; std::map, std::vector>> dataPoints; -void GainMatch::Begin(TTree * /*tree*/) +void GainMatchQQQ::Begin(TTree * /*tree*/) { TString option = GetOption(); @@ -49,7 +49,7 @@ void GainMatch::Begin(TTree * /*tree*/) cut->SetName("qqqcorr"); // Ensure the cut has the correct name } -Bool_t GainMatch::Process(Long64_t entry) +Bool_t GainMatchQQQ::Process(Long64_t entry) { b_sx3Multi->GetEntry(entry); @@ -187,7 +187,7 @@ Bool_t GainMatch::Process(Long64_t entry) return kTRUE; } -void GainMatch::Terminate() +void GainMatchQQQ::Terminate() { const int MAX_DET = 4; const int MAX_RING = 16; @@ -196,7 +196,7 @@ void GainMatch::Terminate() double gainArray[MAX_DET][MAX_RING][MAX_WEDGE] = {{{0}}}; bool gainValid[MAX_DET][MAX_RING][MAX_WEDGE] = {{{false}}}; - std::ofstream outFile("qqq_gainmatch.txt"); + std::ofstream outFile("qqq_GainMatch.txt"); if (!outFile.is_open()) { std::cerr << "Error opening output file!" << std::endl; diff --git a/GainMatchQQQ.h b/GainMatchQQQ.h new file mode 100644 index 0000000..ad06a0d --- /dev/null +++ b/GainMatchQQQ.h @@ -0,0 +1,114 @@ +#ifndef GainMatchQQQ_h +#define GainMatchQQQ_h + +#include +#include +#include +#include + +#include "Armory/ClassDet.h" + +class GainMatchQQQ : 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; //! + + GainMatchQQQ(TTree * /*tree*/ =0) : fChain(0) { } + virtual ~GainMatchQQQ() { } + 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(GainMatchQQQ,0); +}; + +#endif + +#ifdef GainMatchQQQ_cxx +void GainMatchQQQ::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 GainMatchQQQ::Notify(){ + + return kTRUE; +} + +void GainMatchQQQ::SlaveBegin(TTree * /*tree*/){ + + TString option = GetOption(); + +} + +void GainMatchQQQ::SlaveTerminate(){ + +} + + +#endif // #ifdef GainMatchQQQ_cxx \ No newline at end of file diff --git a/GainMatchSX3.C b/GainMatchSX3.C index 388c42f..daf221d 100644 --- a/GainMatchSX3.C +++ b/GainMatchSX3.C @@ -208,14 +208,13 @@ Bool_t GainMatchSX3::Process(Long64_t entry) auto key = std::make_tuple(sx3.id[i], sx3ChBk, sx3ChUp, sx3ChDn); comboCounts[key]++; // If we have a valid front and back channel, fill the histograms - hSX3->Fill(sx3ChDn+4, sx3ChBk); + hSX3->Fill(sx3ChDn + 4, sx3ChBk); hSX3->Fill(sx3ChUp, sx3ChBk); - + // Fill the histogram for the front vs back hSX3FvsB->Fill(sx3EUp + sx3EDn, sx3EBk); } - for (int i = 0; i < sx3.multi; i++) { // if (sx3.id[i] == 4) @@ -262,13 +261,10 @@ Bool_t GainMatchSX3::Process(Long64_t entry) return kTRUE; } +const double GAIN_ACCEPTANCE_THRESHOLD = 0.3; + void GainMatchSX3::Terminate() { - const int MAX_DET = 24; - const int MAX_UP = 4; - const int MAX_DOWN = 4; - const int MAX_BK = 4; - double gainArray[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}}; bool gainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}}; double fbgain[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}}; @@ -334,7 +330,7 @@ void GainMatchSX3::Terminate() xVals.push_back(x); yVals.push_back(y); - exVals.push_back(fixedError); // error in front energy + // exVals.push_back(fixedError); // error in front energy eyVals.push_back(fixedError); // error in back energy } @@ -352,7 +348,7 @@ void GainMatchSX3::Terminate() g.SetMarkerColor(kBlue); g.Draw("AP"); - g.Fit(&f, interactiveMode ? "Q" : "QNR"); // 'R' avoids refit, 'N' skips drawing + g.Fit(&f, interactiveMode ? "Q" : "QNR"); // 'Q': suppress output, 'N': no fit stats box, 'R': avoid refit if (verboseFit) { @@ -394,17 +390,18 @@ void GainMatchSX3::Terminate() // for (int d = 0; d < MAX_DOWN; ++d) // { // // Check if the gain is valid for this detector, back, up, and down - // if (gainValid[id][bk][u][d]) - // { - if (TMath::Abs(gainArray[id][u][d][bk] - 1) < 0.3) + // Only accept gain values within 30% of unity (i.e., 0.7 < gain < 1.3) to filter out unphysical or poorly fitted results. + if (abs(gainArray[id][bk][u][d] - 1) < 0.3) { - printf("Gain match Det%d Up%dDn%d Backs%d → %.4f \n", id, u, d, bk, gainArray[id][u][d][bk]); - outFile << id << " " << bk << " " << u << " " << d << " " << gainArray[id][u][d][bk] << std::endl; + printf("Gain match Det%d Up%dDn%d Backs%d → %.4f \n", id, u, d, bk, gainArray[id][bk][u][d]); + outFile << id << " " << bk << " " << u << " " << d << " " << gainArray[id][bk][u][d] << std::endl; } - else if (gainArray[id][u][d][bk] != 0) + // outFile << id << " " << bk << " " << u << " " << d << " " << gainArray[id][bk][u][d] << std::endl; + // } + else if (gainArray[id][bk][u][d] != 0) { std::cerr << "Warning: Gain value out of range for Det " << id << " Up " << u << " Dn " << d << " Back " << bk << ": " - << gainArray[id][u][d][bk] << std::endl; + << gainArray[id][bk][u][d] << std::endl; } } // }