rtang/ANASEN_analysis
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
ANASEN_analysis/MakeVertexSX3.C
Sudarsan Balakrishnan 7be45f35df Partial commit with key changes as follows: 
- wireOffset is now 4 not 3. Phew. Gosh Golly. Change in Armory/ClassPW.h
	- Concomitant changes made to Armory/PC_Stepladder_Correction.h. In this particular version, a shift of 2 anodes in run*.sh has been replaced by an offset in z
	- The above needs to be worked out carefully.
	- The rho_vs_z fixes have NOT been corrected for. This version will likely yield nonsensical results, hence. To be fixed by a subsequent push that also cleans up the histogram plotting somewhat.

*Read again*

The sort uses rho from the 3 wireOffset cases in MakeVertex.C
The code assumes wireOffset=4.
Anywhere this matters, this version will break - caveat Emptor. Pushing here to plant a flag, but to also make sure the key changes are communicated across branches.
2026-05-14 18:18:17 -04:00

222 lines
8.1 KiB
C
Raw Blame History

#define MakeVertexSX3_cxx
#include "MakeVertexSX3.h"
#include "Armory/ClassPW.h"
#include "Armory/HistPlotter.h"
#include "Armory/SX3Geom.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TBranch.h>
#include "TVector3.h"
#include <fstream>
#include <iostream>
#include <sstream>
#include <vector>
#include <array>
#include <map>
#include <utility>
#include <algorithm>
// Global instances
PW pw_contr;
PW pwinstance;
TVector3 hitPos;
double qqqenergy, qqqtimestamp;
class Event {
public:
Event(TVector3 p, double e1, double e2, double t1, double t2) : pos(p), Energy1(e1), Energy2(e2), Time1(t1), Time2(t2) {}
Event(TVector3 p, double e1, double e2, double t1, double t2, int c1, int c2) : pos(p), Energy1(e1), Energy2(e2), Time1(t1), Time2(t2), ch1(c1), ch2(c2) {}
TVector3 pos;
int ch1=-1; //int(ch1/16) gives qqq id, ch1%16 gives ring#
int ch2=-1; //int(ch2/16) gives qqq id, ch2%16 gives wedge#
double Energy1=-1; //Front for QQQ, Anode for PC
double Energy2=-1; //Back for QQQ, Cathode for PC
double Time1=-1;
double Time2=-1;
};
// Calibration globals
const int MAX_QQQ = 4;
const int MAX_RING = 16;
const int MAX_WEDGE = 16;
double qqqGain[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}};
bool qqqGainValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}};
double qqqCalib[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}};
bool qqqCalibValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}};
// TCutg *cutQQQ;
double sx3BackGain[24][4][4] = {{{1.}}};
double sx3FrontGain[24][4] = {{1.}};
double sx3FrontOffset[24][4] = {{0.}};
double sx3RightGain[24][4] = {{1.}};
// PC Arrays
double pcSlope[48];
double pcIntercept[48];
HistPlotter *plotter;
bool HitNonZero;
bool sx3ecut;
bool qqqEcut;
void MakeVertexSX3::Begin(TTree * /*tree*/)
{
TString option = GetOption();
plotter = new HistPlotter("Analyzer_SX3.root", "TFILE");
pw_contr.ConstructGeo();
pwinstance.ConstructGeo();
for (int i = 0; i < 48; i++)
{
pcSlope[i] = 1.0; // Default slope = 1 (preserves Raw energy)
pcIntercept[i] = 0.0; // Default intercept = 0
}
std::string dataset = "27Al";
{
std::ifstream infile("sx3cal/"+dataset+"/backgains.dat");
std::string temp;
int backpos, frontpos, clkpos;
if (infile.is_open())
while(infile>>clkpos>>temp>>frontpos>>temp>>backpos>>sx3BackGain[clkpos][frontpos][backpos])
;//std::cout << sx3BackGain[clkpos][frontpos][backpos] << std::endl;
infile.close();
infile.open("sx3cal/"+dataset+"/frontgains.dat");
if (infile.is_open())
while(infile>>clkpos>>temp>>temp>>frontpos>>sx3FrontOffset[clkpos][frontpos]>>sx3FrontGain[clkpos][frontpos])
;//std::cout << sx3FrontOffset[clkpos][frontpos] << " " << sx3FrontGain[clkpos][frontpos] << std::endl;
infile.close();
infile.open("sx3cal/"+dataset+"/rightgains.dat");
if (infile.is_open())
while(infile>>clkpos>>frontpos>>temp>>sx3RightGain[clkpos][frontpos]) {
sx3RightGain[clkpos][frontpos]=TMath::Abs(sx3RightGain[clkpos][frontpos]);
}
infile.close();
}
}
Bool_t MakeVertexSX3::Process(Long64_t entry)
{
hitPos.Clear();
qqqenergy = -1;
qqqtimestamp=-1;
HitNonZero = false;
bool qqq1000cut = false;
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<Event> sx3Events;
if(sx3.multi>1) {
std::array<sx3det,12> Fsx3;
//std::cout << "-----" << std::endl;
for(int i=0; i<sx3.multi; i++) {
if(sx3.id[i]>=12) continue;
int id = sx3.id[i];
if(sx3.ch[i]>=8) {
int sx3ch=sx3.ch[i]-8;
sx3ch=(sx3ch+3)%4;
//if(sx3ch==0 || sx3ch==3) continue;
float value=sx3.e[i];
int gch = sx3.id[i]*4+(sx3.ch[i]-8);
Fsx3.at(id).fillevent("BACK",sx3ch,value);
Fsx3.at(id).ts = static_cast<double>(sx3.t[i]);
plotter->Fill2D("sx3backs_raw",100,0,100,800,0,4096,gch,sx3.e[i]);
} else {
int sx3ch=sx3.ch[i]/2;
double value=sx3.e[i];
if(sx3.ch[i]%2==0) {
Fsx3.at(id).fillevent("FRONT_L",sx3ch,value*sx3RightGain[id][sx3ch]);
} else {
Fsx3.at(id).fillevent("FRONT_R",sx3ch,value);
}
}
}
for(int id=0; id<12; id++) {
Fsx3.at(id).validate();
auto det = Fsx3.at(id);
bool no_charge_sharing_strict = det.valid_front_chans.size()==1 && det.valid_back_chans.size()==1;
if(det.valid) {
//std::cout << det.frontEL << " " << det.frontEL*sx3RightGain[id][det.stripF] << std::endl;
plotter->Fill2D("be_vs_x_sx3_id_"+std::to_string(id)+"_f"+std::to_string(det.stripF)+"_b"+std::to_string(det.stripB),200,-1,1,800,0,8192,
det.frontX,det.backE,"evsx");
//std::cout << sx3BackGain[id][det.stripF][det.stripB] << " " << sx3FrontGain[id][det.stripF] << std::endl;
plotter->Fill2D("matched_be_vs_x_sx3_id_"+std::to_string(id)+"_f"+std::to_string(det.stripF),200,-30,30,800,0,8192,
det.frontX*sx3FrontGain[id][det.stripF]+sx3FrontOffset[id][det.stripF],det.backE*sx3BackGain[id][det.stripF][det.stripB],"evsx_matched");
//plotter->Fill2D("fe_vs_x_sx3_id_"+std::to_string(id)+"_f"+std::to_string(det.stripF)+"_"+std::to_string(det.stripB),200,-1,1,800,0,4096,det.frontX,det.backE,"evsx");
plotter->Fill2D("l_vs_r_sx3_id_"+std::to_string(id)+"_f"+std::to_string(det.stripF),800,0,4096,800,0,4096,det.frontEL,det.frontER,"l_vs_r");
}
if(det.valid && (id ==9 || id==7 || id == 1 || id==3) && det.stripF!=DEFAULT_NULL && det.stripB!=DEFAULT_NULL) {
double z = det.frontX*sx3FrontGain[id][det.stripF]+sx3FrontOffset[id][det.stripF];
double backE = det.backE*sx3BackGain[id][det.stripF][det.stripB];
Event sx3ev(TVector3(0,0,z),backE,-1,det.ts,-1,det.stripB+4*id,det.stripF+4*id);
sx3Events.push_back(sx3ev);
}
if(!det.valid && det.valid_back_chans.size()) {
for(auto bch : det.valid_back_chans) {
plotter->Fill1D("good_backs_bad_fronts_id"+std::to_string(id)+"_b"+std::to_string(bch),800,0,4096,det.back[bch].at(0),"good_backs_no_front");
}
}
if(!det.valid && det.unmatched_front_chans.size()) {
for(auto bch : det.valid_back_chans)
for(auto fch : det.unmatched_front_chans) {
if(det.frontR[fch].size() != 0) {
double trial_left = det.back[bch].at(0) - det.frontR[fch].at(0);
double frontR = det.frontR[fch].at(0);
double backE = det.back[bch].at(0);
double frontX = (trial_left - frontR)/(trial_left + frontR);
plotter->Fill2D("bad_fronts_id"+std::to_string(id)+"_fR"+std::to_string(fch)+"_b"+std::to_string(bch),800,0,4096,800,0,4096,det.frontR[fch].at(0),det.back[bch].at(0),"bad_fronts_vb");
plotter->Fill2D("be_vs_x_sx3_id_"+std::to_string(id)+"_f"+std::to_string(fch)+"R_b"+std::to_string(bch),200,-1,1,800,0,8192,
frontX,backE,"evsx_bad_front");
//double backE = det.back[bch].at(0)*sx3BackGain[id][fch][bch];
//double frontE = det.frontR[fch]
}
if(det.frontL[fch].size() != 0) {
double trial_right = det.back[bch].at(0) - det.frontL[fch].at(0);
double frontL = det.frontL[fch].at(0);
double backE = det.back[bch].at(0);
double frontX = (frontL - trial_right)/(frontL + trial_right);
plotter->Fill2D("bad_fronts_id"+std::to_string(id)+"_fL"+std::to_string(fch)+"_b"+std::to_string(bch),800,0,4096,800,0,4096,det.frontL[fch].at(0),det.back[bch].at(0),"bad_fronts_vb");
plotter->Fill2D("be_vs_x_sx3_id_"+std::to_string(id)+"_f"+std::to_string(fch)+"L_b"+std::to_string(bch),200,-1,1,800,0,8192,
frontX,backE,"evsx_bad_front");
}
}//end for
}//end if(!det.valid)
}//end for id
}//end sx3.multi>1
return kTRUE;
}
void MakeVertexSX3::Terminate()
{
plotter->FlushToDisk(10);
}
Reference in New Issue View Git Blame Copy Permalink