new file: MakeVertex.C
new file: MakeVertex.C:Zone.Identifier new file: MakeVertex.h new file: MakeVertex.h:Zone.Identifier modified: TrackRecon.C
This commit is contained in:
parent
282aa5ecea
commit
4401ae8eb2
834
MakeVertex.C
Normal file
834
MakeVertex.C
Normal file
|
|
@ -0,0 +1,834 @@
|
|||
#define MakeVertex_cxx
|
||||
|
||||
#include "MakeVertex.h"
|
||||
#include "Armory/ClassPW.h"
|
||||
#include "Armory/HistPlotter.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 <map>
|
||||
#include <utility>
|
||||
#include <algorithm>
|
||||
|
||||
// Global instances
|
||||
PW pw_contr;
|
||||
PW pwinstance;
|
||||
TVector3 hitPos;
|
||||
double qqqenergy, qqqtimestamp;
|
||||
|
||||
struct Event {
|
||||
TVector3 pos;
|
||||
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;
|
||||
|
||||
// PC Arrays
|
||||
double pcSlope[48];
|
||||
double pcIntercept[48];
|
||||
|
||||
HistPlotter *plotter;
|
||||
|
||||
bool HitNonZero;
|
||||
bool sx3ecut;
|
||||
bool qqqEcut;
|
||||
|
||||
void MakeVertex::Begin(TTree * /*tree*/)
|
||||
{
|
||||
TString option = GetOption();
|
||||
plotter = new HistPlotter("Analyzer_QQQ.root", "TFILE");
|
||||
|
||||
pw_contr.ConstructGeo();
|
||||
pwinstance.ConstructGeo();
|
||||
|
||||
// ---------------------------------------------------------
|
||||
// 1. CRITICAL FIX: Initialize PC Arrays to Default (Raw)
|
||||
// ---------------------------------------------------------
|
||||
for (int i = 0; i < 48; i++)
|
||||
{
|
||||
pcSlope[i] = 1.0; // Default slope = 1 (preserves Raw energy)
|
||||
pcIntercept[i] = 0.0; // Default intercept = 0
|
||||
}
|
||||
|
||||
// Calculate Crossover Geometry ONCE
|
||||
TVector3 a, c, diff;
|
||||
double a2, ac, c2, adiff, cdiff, denom, alpha;
|
||||
|
||||
for (size_t i = 0; i < pwinstance.An.size(); i++)
|
||||
{
|
||||
a = pwinstance.An[i].first - pwinstance.An[i].second;
|
||||
|
||||
for (size_t j = 0; j < pwinstance.Ca.size(); j++)
|
||||
{
|
||||
c = pwinstance.Ca[j].first - pwinstance.Ca[j].second;
|
||||
diff = pwinstance.An[i].first - pwinstance.Ca[j].first;
|
||||
a2 = a.Dot(a);
|
||||
c2 = c.Dot(c);
|
||||
ac = a.Dot(c);
|
||||
adiff = a.Dot(diff);
|
||||
cdiff = c.Dot(diff);
|
||||
denom = a2 * c2 - ac * ac;
|
||||
alpha = (ac * cdiff - c2 * adiff) / denom;
|
||||
|
||||
Crossover[i][j][0].x = pwinstance.An[i].first.X() + alpha * a.X();
|
||||
Crossover[i][j][0].y = pwinstance.An[i].first.Y() + alpha * a.Y();
|
||||
Crossover[i][j][0].z = pwinstance.An[i].first.Z() + alpha * a.Z();
|
||||
|
||||
if (Crossover[i][j][0].z < -190 || Crossover[i][j][0].z > 190 || (i+j)%24 == 12)
|
||||
{
|
||||
Crossover[i][j][0].z = 9999999;
|
||||
}
|
||||
|
||||
Crossover[i][j][1].x = alpha;
|
||||
Crossover[i][j][1].y = 0;
|
||||
}
|
||||
}
|
||||
|
||||
// Load PC Calibrations
|
||||
std::ifstream inputFile("slope_intercept_results.txt");
|
||||
if (inputFile.is_open())
|
||||
{
|
||||
std::string line;
|
||||
int index;
|
||||
double slope, intercept;
|
||||
while (std::getline(inputFile, line))
|
||||
{
|
||||
std::stringstream ss(line);
|
||||
ss >> index >> slope >> intercept;
|
||||
if (index >= 0 && index <= 47)
|
||||
{
|
||||
pcSlope[index] = slope;
|
||||
pcIntercept[index] = intercept;
|
||||
}
|
||||
}
|
||||
inputFile.close();
|
||||
}
|
||||
else
|
||||
{
|
||||
std::cerr << "Error opening slope_intercept.txt" << std::endl;
|
||||
}
|
||||
|
||||
// Load QQQ Cuts from file
|
||||
// {
|
||||
// std::string filename = "QQQ_PCCut.root";
|
||||
// TFile *cutFile = TFile::Open(filename.c_str(), "READ");
|
||||
// if (cutFile && !cutFile->IsZombie())
|
||||
// {
|
||||
// cutQQQ = (TCutg *)cutFile->Get("cutQQQPC");
|
||||
// if (cutQQQ)
|
||||
// {
|
||||
// std::cout << "Loaded QQQ PC cut from " << filename << std::endl;
|
||||
// }
|
||||
// else
|
||||
// {
|
||||
// std::cerr << "Error: cutQQQPC not found in " << filename << std::endl;
|
||||
// }
|
||||
// cutFile->Close();
|
||||
// }
|
||||
// }
|
||||
|
||||
// ... (Load QQQ Gains and Calibs - same as before) ...
|
||||
{
|
||||
std::string filename = "qqq_GainMatch.dat";
|
||||
std::ifstream infile(filename);
|
||||
if (infile.is_open())
|
||||
{
|
||||
int det, ring, wedge;
|
||||
double gainw, gainr;
|
||||
while (infile >> det >> wedge >> ring >> gainw >> gainr)
|
||||
{
|
||||
qqqGain[det][wedge][ring] = gainw;
|
||||
qqqGainValid[det][wedge][ring] = (gainw > 0);
|
||||
// std::cout << "QQQ Gain Loaded: Det " << det << " Ring " << ring << " Wedge " << wedge << " GainW " << gainw << " GainR " << gainr << std::endl;
|
||||
}
|
||||
infile.close();
|
||||
}
|
||||
}
|
||||
{
|
||||
std::string filename = "qqq_Calib.dat";
|
||||
std::ifstream infile(filename);
|
||||
if (infile.is_open())
|
||||
{
|
||||
int det, ring, wedge;
|
||||
double slope;
|
||||
while (infile >> det >> wedge >> ring >> slope)
|
||||
{
|
||||
qqqCalib[det][wedge][ring] = slope;
|
||||
qqqCalibValid[det][wedge][ring] = (slope > 0);
|
||||
// std::cout << "QQQ Calib Loaded: Det " << det << " Ring " << ring << " Wedge " << wedge << " Slope " << slope << std::endl;
|
||||
}
|
||||
infile.close();
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
Bool_t MakeVertex::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();
|
||||
|
||||
// QQQ Processing
|
||||
|
||||
int qqqCount = 0;
|
||||
int qqqAdjCh = 0;
|
||||
// REMOVE WHEN RERUNNING USING THE NEW CALIBRATION FILE
|
||||
for (int i = 0; i < qqq.multi; i++)
|
||||
{
|
||||
//if ((qqq.id[i] == 3 || qqq.id[i] == 1) && qqq.ch[i] < 16)
|
||||
if (qqq.id[i] == 1 && qqq.ch[i] < 16) //for run 12, 26Al
|
||||
{
|
||||
qqq.ch[i] = 16 - qqq.ch[i];
|
||||
}
|
||||
}
|
||||
for (int i = 0; i < qqq.multi; i++)
|
||||
{
|
||||
if (qqq.id[i] == 0 && qqq.ch[i] >= 16)
|
||||
{
|
||||
qqq.ch[i] = 31 - qqq.ch[i] + 16;
|
||||
}
|
||||
}
|
||||
|
||||
std::vector<Event> QQQ_Events, PC_Events;
|
||||
std::vector<Event> QQQ_Events_Raw, PC_Events_Raw;
|
||||
|
||||
|
||||
bool PCQQQTimeCut = false;
|
||||
for (int i = 0; i < qqq.multi; i++)
|
||||
{
|
||||
|
||||
plotter->Fill2D("QQQ_Index_Vs_Energy", 16 * 8, 0, 16 * 8, 2000, 0, 16000, qqq.index[i], qqq.e[i], "hRawQQQ");
|
||||
|
||||
for (int j = 0; j < qqq.multi; j++)
|
||||
{
|
||||
if (j == i)
|
||||
continue;
|
||||
plotter->Fill2D("QQQ_Coincidence_Matrix", 16 * 8, 0, 16 * 8, 16 * 8, 0, 16 * 8, qqq.index[i], qqq.index[j], "hRawQQQ");
|
||||
}
|
||||
|
||||
for (int k = 0; k < pc.multi; k++)
|
||||
{
|
||||
if (pc.index[k] < 24 && pc.e[k] > 50)
|
||||
{
|
||||
plotter->Fill2D("QQQ_Vs_Anode_Energy", 400, 0, 4000, 1000, 0, 16000, qqq.e[i], pc.e[k], "hRawQQQ");
|
||||
plotter->Fill2D("QQQ_Vs_PC_Index", 16 * 8, 0, 16 * 8, 24, 0, 24, qqq.index[i], pc.index[k], "hRawQQQ");
|
||||
}
|
||||
else if (pc.index[k] >= 24 && pc.e[k] > 50)
|
||||
{
|
||||
plotter->Fill2D("QQQ_Vs_Cathode_Energy", 400, 0, 4000, 1000, 0, 16000, qqq.e[i], pc.e[k], "hRawQQQ");
|
||||
}
|
||||
}
|
||||
|
||||
for (int j = i + 1; j < qqq.multi; j++)
|
||||
{
|
||||
if (qqq.id[i] == qqq.id[j])
|
||||
{
|
||||
qqqCount++;
|
||||
|
||||
int chWedge = -1;
|
||||
int chRing = -1;
|
||||
double eWedge = 0.0;
|
||||
double eWedgeMeV = 0.0;
|
||||
double eRing = 0.0;
|
||||
double eRingMeV = 0.0;
|
||||
double tRing = 0.0;
|
||||
double tWedge = 0.0;
|
||||
|
||||
if (qqq.ch[i] < 16 && qqq.ch[j] >= 16 && qqqGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16])
|
||||
{
|
||||
chWedge = qqq.ch[i];
|
||||
eWedge = qqq.e[i] * qqqGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16];
|
||||
chRing = qqq.ch[j] - 16;
|
||||
eRing = qqq.e[j];
|
||||
tRing = static_cast<double>(qqq.t[j]);
|
||||
tWedge = static_cast<double>(qqq.t[i]);
|
||||
}
|
||||
else if (qqq.ch[j] < 16 && qqq.ch[i] >= 16 && qqqGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16])
|
||||
{
|
||||
chWedge = qqq.ch[j];
|
||||
eWedge = qqq.e[j] * qqqGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16];
|
||||
chRing = qqq.ch[i] - 16;
|
||||
eRing = qqq.e[i];
|
||||
tRing = static_cast<double>(qqq.t[i]);
|
||||
tWedge = static_cast<double>(qqq.t[j]);
|
||||
}
|
||||
else
|
||||
continue;
|
||||
|
||||
plotter->Fill1D("Wedgetime_Vs_Ringtime", 100, -1000, 1000, tWedge - tRing, "hTiming");
|
||||
plotter->Fill2D("RingE_vs_Index", 16 * 4, 0, 16 * 4, 1000, 0, 16000, chRing + qqq.id[i] * 16, eRing, "hRawQQQ");
|
||||
plotter->Fill2D("WedgeE_vs_Index", 16 * 4, 0, 16 * 4, 1000, 0, 16000, chWedge + qqq.id[i] * 16, eWedge, "hRawQQQ");
|
||||
|
||||
if (qqqCalibValid[qqq.id[i]][chWedge][chRing])
|
||||
{
|
||||
eWedgeMeV = eWedge * qqqCalib[qqq.id[i]][chWedge][chRing] / 1000;
|
||||
eRingMeV = eRing * qqqCalib[qqq.id[i]][chWedge][chRing] / 1000;
|
||||
}
|
||||
else
|
||||
continue;
|
||||
plotter->Fill2D("WedgeE_Vs_RingECal", 1000, 0, 10, 1000, 0, 10, eWedgeMeV, eRingMeV, "hCalQQQ");
|
||||
if(eRingMeV/eWedgeMeV > 3.0 || eRingMeV/eWedgeMeV<1.0/3.0) continue;
|
||||
plotter->Fill2D("WedgeE_Vs_RingECal_selected", 1000, 0, 10, 1000, 0, 10, eWedgeMeV, eRingMeV, "hCalQQQ");
|
||||
|
||||
for (int k = 0; k < pc.multi; k++)
|
||||
{
|
||||
plotter->Fill2D("RingCh_vs_Anode_Index", 16 * 4, 0, 16 * 4, 24, 0, 24, chRing + qqq.id[i] * 16, pc.index[k], "hRawQQQ");
|
||||
plotter->Fill2D("WedgeCh_vs_Anode_Index", 16 * 4, 0, 16 * 4, 24, 0, 24, chWedge + qqq.id[i] * 16, pc.index[k], "hRawQQQ");
|
||||
plotter->Fill2D("WedgeCh_vs_Anode_Index" + std::to_string(qqq.id[i]), 16 * 4, 0, 16 * 4, 24, 0, 24, chWedge + qqq.id[i] * 16, pc.index[k]);
|
||||
plotter->Fill2D("RingCh_vs_Cathode_Index", 16 * 4, 0, 16 * 4, 24, 24, 48, chRing + qqq.id[i] * 16, pc.index[k], "hRawQQQ");
|
||||
plotter->Fill2D("WedgeCh_vs_Cathode_Index", 16 * 4, 0, 16 * 4, 24, 24, 48, chWedge + qqq.id[i] * 16, pc.index[k], "hRawQQQ");
|
||||
|
||||
if (pc.index[k] < 24 && pc.e[k] > 50)
|
||||
{
|
||||
// plotter->Fill2D("QQQ_CalibW_Vs_PC_Energy", 1000, 0, 16, 2000, 0, 30000, eWedgeMeV, pc.e[k], "hCalQQQ");
|
||||
// plotter->Fill2D("QQQ_CalibR_Vs_PC_Energy", 1000, 0, 16, 2000, 0, 30000, eRingMeV, pc.e[k], "hCalQQQ");
|
||||
|
||||
// if (tRing - static_cast<double>(pc.t[k]) < 0 && tRing - static_cast<double>(pc.t[k]) > -600)
|
||||
// // {
|
||||
// // plotter->Fill2D("QQQ_CalibW_Vs_PC_Energy_Tight", 1000, 0, 16, 2000, 0, 30000, eWedgeMeV, pc.e[k], "hCalQQQ");
|
||||
// // plotter->Fill2D("QQQ_CalibR_Vs_PC_Energy_Tight", 1000, 0, 16, 2000, 0, 30000, eRingMeV, pc.e[k], "hCalQQQ");
|
||||
// // }
|
||||
// // else
|
||||
// // {
|
||||
// // plotter->Fill2D("QQQ_CalibW_Vs_PC_Energy_OffTime", 1000, 0, 16, 2000, 0, 30000, eWedgeMeV, pc.e[k], "hCalQQQ");
|
||||
// // plotter->Fill2D("QQQ_CalibR_Vs_PC_Energy_OffTime", 1000, 0, 16, 2000, 0, 30000, eRingMeV, pc.e[k], "hCalQQQ");
|
||||
// // }
|
||||
plotter->Fill2D("Timing_Difference_QQQ_PC", 500, -2000, 2000, 16, 0, 16, tRing - static_cast<double>(pc.t[k]), chRing, "hTiming");
|
||||
plotter->Fill2D("DelT_Vs_QQQRingECal", 500, -2000, 2000, 1000, 0, 10, tRing - static_cast<double>(pc.t[k]), eRingMeV, "hTiming");
|
||||
plotter->Fill2D("CalibratedQQQEvsPCE_R", 1000, 0, 10, 2000, 0, 30000, eRingMeV, pc.e[k], "hPCQQQ");
|
||||
plotter->Fill2D("CalibratedQQQEvsPCE_W", 1000, 0, 10, 2000, 0, 30000, eWedgeMeV, pc.e[k], "hPCQQQ");
|
||||
if (tRing - static_cast<double>(pc.t[k]) < -150) // proton tests, 27Al
|
||||
|
||||
//if (tRing - static_cast<double>(pc.t[k]) < -150 && tRing - static_cast<double>(pc.t[k]) > -450) // 27Al
|
||||
//if (tRing - static_cast<double>(pc.t[k]) < -70 && tRing - static_cast<double>(pc.t[k]) > -150) // 17F
|
||||
{
|
||||
PCQQQTimeCut = true;
|
||||
}
|
||||
}
|
||||
|
||||
if (pc.index[k] >= 24 && pc.e[k] > 50) {
|
||||
plotter->Fill2D("Timing_Difference_QQQ_PC_Cathode", 500, -2000, 2000, 16, 0, 16, tRing - static_cast<double>(pc.t[k]), chRing, "hTiming");
|
||||
}
|
||||
} //end of pc loop
|
||||
|
||||
double theta = -TMath::Pi() / 2 + 2 * TMath::Pi() / 16 / 4. * (qqq.id[i] * 16 + chWedge + 0.5);
|
||||
double rho = 50. + (50. / 16.) * (chRing + 0.5); //"?"
|
||||
|
||||
//Event qqqevent(TVector3(rho*TMath::Cos(theta),rho*TMath::Sin(theta),23+75+30), eRingMeV, eWedgeMeV, tRing, tWedge);
|
||||
//Event qqqeventr(TVector3(rho*TMath::Cos(theta),rho*TMath::Sin(theta),23+75+30), eRing, eWedge, tRing, tWedge);
|
||||
//QQQ_Events.push_back(qqqevent);
|
||||
//QQQ_Events_Raw.push_back(qqqeventr);
|
||||
plotter->Fill2D("QQQPolarPlot", 16 * 4, -TMath::Pi(), TMath::Pi(), 32, 40, 100, theta, rho, "hCalQQQ");
|
||||
plotter->Fill2D("QQQCartesianPlot", 200, -100, 100, 200, -100, 100, rho * TMath::Cos(theta), rho * TMath::Sin(theta), "hCalQQQ");
|
||||
plotter->Fill2D("QQQCartesianPlot" + std::to_string(qqq.id[i]), 200, -100, 100, 200, -100, 100, rho * TMath::Cos(theta), rho * TMath::Sin(theta), "hCalQQQ");
|
||||
if (PCQQQTimeCut)
|
||||
{
|
||||
plotter->Fill2D("PC_XY_Projection_QQQ_TimeCut" + std::to_string(qqq.id[i]), 400, -100, 100, 400, -100, 100, rho * TMath::Cos(theta), rho * TMath::Sin(theta), "hPCQQQ");
|
||||
}
|
||||
plotter->Fill2D("PC_XY_Projection_QQQ" + std::to_string(qqq.id[i]), 400, -100, 100, 400, -100, 100, rho * TMath::Cos(theta), rho * TMath::Sin(theta), "hPCQQQ");
|
||||
|
||||
if (!HitNonZero)
|
||||
{
|
||||
double x = rho * TMath::Cos(theta);
|
||||
double y = rho * TMath::Sin(theta);
|
||||
hitPos.SetXYZ(x, y, (23 + 75 + 30));
|
||||
qqqenergy = eRingMeV;
|
||||
qqqtimestamp = tRing;
|
||||
HitNonZero = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
plotter->Fill1D("QQQ_Multiplicity", 10, 0, 10, qqqCount, "hRawQQQ");
|
||||
|
||||
// PC Gain Matching and Filling
|
||||
double anodeT = -99999;
|
||||
double cathodeT = 99999;
|
||||
int anodeIndex = -1;
|
||||
int cathodeIndex = -1;
|
||||
for (int i = 0; i < pc.multi; i++)
|
||||
{
|
||||
if (pc.e[i] > 50)
|
||||
{
|
||||
plotter->Fill2D("PC_Index_Vs_Energy", 48, 0, 48, 2000, 0, 30000, pc.index[i], static_cast<double>(pc.e[i]), "hRawPC");
|
||||
} else {
|
||||
continue;
|
||||
}
|
||||
|
||||
if (pc.index[i] < 48)
|
||||
{
|
||||
pc.e[i] = pcSlope[pc.index[i]] * pc.e[i] + pcIntercept[pc.index[i]];
|
||||
plotter->Fill2D("PC_Index_VS_GainMatched_Energy", 48, 0, 48, 2000, 0, 30000, pc.index[i], pc.e[i], "hGMPC");
|
||||
}
|
||||
|
||||
if (pc.index[i] < 24)
|
||||
{
|
||||
anodeT = static_cast<double>(pc.t[i]);
|
||||
anodeIndex = pc.index[i];
|
||||
}
|
||||
else
|
||||
{
|
||||
cathodeT = static_cast<double>(pc.t[i]);
|
||||
cathodeIndex = pc.index[i] - 24;
|
||||
}
|
||||
|
||||
if (anodeT != -99999 && cathodeT != 99999)
|
||||
{
|
||||
for (int j = 0; j < qqq.multi; j++)
|
||||
{
|
||||
plotter->Fill1D("PC_Time_qqq", 200, -2000, 2000, anodeT - cathodeT, "hTiming");
|
||||
plotter->Fill2D("PC_Time_Vs_QQQ_ch", 200, -2000, 2000, 16 * 8, 0, 16 * 8, anodeT - cathodeT, qqq.ch[j], "hTiming");
|
||||
plotter->Fill2D("PC_Time_vs_AIndex", 200, -2000, 2000, 24, 0, 24, anodeT - cathodeT, anodeIndex, "hTiming");
|
||||
plotter->Fill2D("PC_Time_vs_CIndex", 200, -2000, 2000, 24, 0, 24, anodeT - cathodeT, cathodeIndex, "hTiming");
|
||||
// plotter->Fill1D("PC_Time_A" + std::to_string(anodeIndex) + "_C" + std::to_string(cathodeIndex), 200, -1000, 1000, anodeT - cathodeT, "TimingPC");
|
||||
}
|
||||
|
||||
for (int j = 0; j < sx3.multi; j++)
|
||||
{
|
||||
plotter->Fill1D("PC_Time_sx3", 200, -2000, 2000, anodeT - cathodeT, "hTiming");
|
||||
}
|
||||
|
||||
plotter->Fill1D("PC_Time", 200, -2000, 2000, anodeT - cathodeT, "hTiming");
|
||||
}
|
||||
|
||||
for (int j = i + 1; j < pc.multi; j++)
|
||||
{
|
||||
plotter->Fill2D("PC_Coincidence_Matrix", 48, 0, 48, 48, 0, 48, pc.index[i], pc.index[j], "hRawPC");
|
||||
plotter->Fill2D("PC_Coincidence_Matrix_anodeMinusCathode_lt_-200_" + std::to_string(anodeT - cathodeT < -200), 48, 0, 48, 48, 0, 48, pc.index[i], pc.index[j], "hRawPC");
|
||||
plotter->Fill2D("Anode_V_Anode", 24, 0, 24, 24, 0, 24, pc.index[i], pc.index[j], "hGMPC");
|
||||
}
|
||||
}
|
||||
|
||||
anodeHits.clear();
|
||||
cathodeHits.clear();
|
||||
corrcatMax.clear();
|
||||
|
||||
int aID = 0;
|
||||
int cID = 0;
|
||||
double aE = 0;
|
||||
double cE = 0;
|
||||
double aESum = 0;
|
||||
double cESum = 0;
|
||||
double aEMax = 0;
|
||||
int aIDMax = 0;
|
||||
|
||||
typedef std::unordered_map<int,std::tuple<int,double,double>> WireEvent; //this stores nearest neighbour wire events, or a 'cluster'
|
||||
WireEvent aWireEvents, cWireEvents, aWireEvent, cWireEvent; //naming for book keeping
|
||||
std::vector<WireEvent> aClusters, cClusters; //all clusters that fire toghther. aClusters.at(0).size() gives 'how many wires in cluster-0', aClusters.size() gives 'how many clusters'
|
||||
|
||||
std::set<int> awirelist, cwirelist;
|
||||
for (int i = 0; i < pc.multi; i++)
|
||||
{
|
||||
// if (pc.e[i] > 100)
|
||||
{
|
||||
if (pc.index[i] < 24) {
|
||||
anodeHits.push_back(std::pair<int, double>(pc.index[i], pc.e[i]));
|
||||
//aWireEvents.push_back(std::tuple(pc.index[i],pc.t[i],pc.e[i]));
|
||||
//awirelist.insert(pc.index[i]);
|
||||
}
|
||||
else if (pc.index[i] >= 24) {
|
||||
cathodeHits.push_back(std::pair<int, double>(pc.index[i] - 24, pc.e[i]));
|
||||
//cxWireEvents.push_back(std::tuple(pc.index[i]-24,pc.t[i],pc.e[i]));
|
||||
//cwirelist.insert(pc.index[i]-24);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
std::sort(anodeHits.begin(),anodeHits.end(),[](std::pair<int,double> a, std::pair<int,double> b){ return a.first < b.first;});
|
||||
std::sort(cathodeHits.begin(),cathodeHits.end(),[](std::pair<int,double> a, std::pair<int,double> b){ return a.first < b.first;});
|
||||
|
||||
//std::sort(cWireEvents.begin(),cWireEvents.end(),[](std::tuple<int,double,double> a, std::tuple<int,double,double> b){return std::get<0>(a) < std::get<0>(b);});
|
||||
//std::sort(aWireEvents.begin(),aWireEvents.end(),[](std::tuple<int,double,double> a, std::tuple<int,double,double> b){return std::get<0>(a) < std::get<0>(b);});
|
||||
/*
|
||||
|
||||
for(int i=0; i<24; i++) {
|
||||
if(std::find_if(awireEvents.begin(),aWireEvents.end(),[&](std::tuple<int,double,double> a){return std::get<0>(a)==i;})) {
|
||||
aClusters.
|
||||
}
|
||||
}
|
||||
|
||||
std::vector<std::vector<std::tuple<int,double,double>>> aWires, cWires;
|
||||
std::set<int> acs, bcs;
|
||||
|
||||
for(size_t i=0; i<anodeHits.size() ; i++) {
|
||||
if(i==0) {
|
||||
aWires.push_back(anodeHits.at(i)); //store the first one
|
||||
acs.push_back(anodeHits.at(i).first); //store channel for easy lookup
|
||||
} else {
|
||||
if(acs.find(anodeHits.at(i).first-1)!=acs.end()) {
|
||||
aWires.push_back(anodeHits.at(i)); //store the first one
|
||||
acs.push_back(anodeHits.at(i).first); //store channel for easy lookup
|
||||
} else {
|
||||
|
||||
}
|
||||
}
|
||||
}*/
|
||||
|
||||
|
||||
if (anodeHits.size() >= 1 && cathodeHits.size() >= 1)
|
||||
{
|
||||
// 2. CRITICAL FIX: Define reference vector 'a'
|
||||
// In Analyzer.cxx, 'a' was left over from the loop. We use the first anode wire as reference here.
|
||||
// (Assuming pwinstance.An is populated and wires are generally parallel).
|
||||
TVector3 refAnode = pwinstance.An[0].first - pwinstance.An[0].second;
|
||||
|
||||
{
|
||||
for (const auto &anode : anodeHits)
|
||||
{
|
||||
aID = anode.first;
|
||||
aE = anode.second;
|
||||
aESum += aE;
|
||||
if (aE > aEMax)
|
||||
{
|
||||
aEMax = aE;
|
||||
aIDMax = aID;
|
||||
}
|
||||
}
|
||||
|
||||
for (const auto &cathode : cathodeHits)
|
||||
{
|
||||
cID = cathode.first;
|
||||
cE = cathode.second;
|
||||
plotter->Fill2D("AnodeMax_Vs_Cathode_Coincidence_Matrix", 24, 0, 24, 24, 0, 24, aIDMax, cID, "hRawPC");
|
||||
plotter->Fill2D("Anode_Vs_Cathode_Coincidence_Matrix", 24, 0, 24, 24, 0, 24, aID, cID, "hRawPC");
|
||||
plotter->Fill2D("Anode_vs_CathodeE", 2000, 0, 30000, 2000, 0, 30000, aE, cE, "hGMPC");
|
||||
plotter->Fill2D("CathodeMult_V_CathodeE", 6, 0, 6, 2000, 0, 30000, cathodeHits.size(), cE, "hGMPC");
|
||||
for (int j = -4; j < 3; j++)
|
||||
{
|
||||
if ((aIDMax + 24 + j) % 24 == 23 - cID)
|
||||
{
|
||||
corrcatMax.push_back(std::pair<int, double>(cID, cE));
|
||||
cESum += cE;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
TVector3 anodeIntersection,vector_closest_to_z;
|
||||
anodeIntersection.Clear();
|
||||
vector_closest_to_z.Clear();
|
||||
if (corrcatMax.size() > 0)
|
||||
{
|
||||
double x = 0, y = 0, z = 0;
|
||||
for (const auto &corr : corrcatMax)
|
||||
{
|
||||
if (Crossover[aIDMax][corr.first][0].z > 9000000)
|
||||
continue;
|
||||
if (cESum > 0)
|
||||
{
|
||||
x += (corr.second) / cESum * Crossover[aIDMax][corr.first][0].x;
|
||||
y += (corr.second) / cESum * Crossover[aIDMax][corr.first][0].y;
|
||||
z += (corr.second) / cESum * Crossover[aIDMax][corr.first][0].z;
|
||||
}
|
||||
}
|
||||
if (x == 0 && y == 0 && z == 0)
|
||||
;
|
||||
// to ignore events with no valid crossover points
|
||||
else
|
||||
anodeIntersection = TVector3(x, y, z);
|
||||
// << "Anode Intersection: " << anodeIntersection.X() << ", " << anodeIntersection.Y() << ", " << anodeIntersection.Z() << std::endl;
|
||||
}
|
||||
bool PCQQQPhiCut = false;
|
||||
// flip the algorithm for cathode 1 multi anode events
|
||||
if ((hitPos.Phi() > (anodeIntersection.Phi() - TMath::PiOver4())) && (hitPos.Phi() < (anodeIntersection.Phi() + TMath::PiOver4()))) {
|
||||
PCQQQPhiCut = true;
|
||||
}
|
||||
|
||||
for (double Tz = 60; Tz <= 100; Tz += 1.0)
|
||||
{
|
||||
TVector3 TargetPos(0, 0, Tz);
|
||||
if(PCQQQPhiCut && anodeIntersection.Perp()>0 && anodeIntersection.Z()!=0 && cathodeHits.size()>=2) {
|
||||
plotter->Fill2D("Inttheta_vs_QQQtheta_TC" + std::to_string(PCQQQTimeCut) + "_TZ" + std::to_string(Tz), 400, 0, 180, 90, 0, 90, (anodeIntersection - TargetPos).Theta() * 180. / TMath::Pi(), (hitPos - TargetPos).Theta() * 180. / TMath::Pi(), "TPosVariation");
|
||||
//plotter->Fill2D("R_ratio_to_Z_ratio" + std::to_string(PCQQQTimeCut) + "_TZ" + std::to_string(Tz), 100, -2, 2, 100, -2, 2, (anodeIntersection - TargetPos).Z()/(hitPos-TargetPos).Z(), ((anodeIntersection - TargetPos).Perp()+2.5)/(hitPos-TargetPos).Perp(), "TPosVariation");
|
||||
}
|
||||
}
|
||||
|
||||
if (anodeIntersection.Z() != 0 && anodeIntersection.Perp()>0 && HitNonZero)
|
||||
{
|
||||
plotter->Fill1D("PC_Z_Projection", 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
|
||||
plotter->Fill2D("Z_Proj_VsDelTime", 600, -300, 300, 200, -2000, 2000, anodeIntersection.Z(), anodeT - cathodeT, "hPCzQQQ");
|
||||
plotter->Fill2D("IntPhi_vs_QQQphi", 100, -200, 200, 80, -200, 200, anodeIntersection.Phi() * 180. / TMath::Pi(), hitPos.Phi() * 180. / TMath::Pi(), "hPCQQQ");
|
||||
//plotter->Fill2D("Inttheta_vs_QQQtheta", 90, 0, 180, 20, 0, 45, anodeIntersection.Theta() * 180. / TMath::Pi(), hitPos.Theta() * 180. / TMath::Pi(), "hPCQQQ");
|
||||
//plotter->Fill2D("Inttheta_vs_QQQtheta_TC" + std::to_string(PCQQQTimeCut)+ "_PC"+std::to_string(PCQQQPhiCut), 90, 0, 180, 20, 0, 45, anodeIntersection.Theta() * 180. / TMath::Pi(), hitPos.Theta() * 180. / TMath::Pi(), "hPCQQQ");
|
||||
plotter->Fill2D("IntPhi_vs_QQQphi_TC" + std::to_string(PCQQQTimeCut) + "PhiC" + std::to_string(PCQQQPhiCut), 100, -200, 200, 80, -200, 200, anodeIntersection.Phi() * 180. / TMath::Pi(), hitPos.Phi() * 180. / TMath::Pi(), "hPCQQQ");
|
||||
}
|
||||
if (anodeIntersection.Z() != 0 && cathodeHits.size() >= 2)
|
||||
plotter->Fill1D("PC_Z_Projection_TC" + std::to_string(PCQQQTimeCut) + "PhiC" + std::to_string(PCQQQPhiCut), 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
|
||||
|
||||
if (anodeIntersection.Z() != 0 && cathodeHits.size() == 1)
|
||||
{
|
||||
plotter->Fill1D("PC_Z_proj_1C", 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
|
||||
plotter->Fill2D("IntersectionPhi_vs_AnodeZ_1C", 400, -200, 200, 600, -300, 300, anodeIntersection.Phi() * 180. / TMath::Pi(), anodeIntersection.Z(), "hPCzQQQ");
|
||||
}
|
||||
|
||||
if (anodeIntersection.Z() != 0 && cathodeHits.size() == 2)
|
||||
{
|
||||
plotter->Fill1D("PC_Z_proj_2C", 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
|
||||
plotter->Fill2D("IntersectionPhi_vs_AnodeZ_2C", 400, -200, 200, 600, -300, 300, anodeIntersection.Phi() * 180. / TMath::Pi(), anodeIntersection.Z(), "hGMPC");
|
||||
}
|
||||
if (anodeIntersection.Z() != 0 && cathodeHits.size() > 2)
|
||||
{
|
||||
plotter->Fill1D("PC_Z_proj_nC", 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
|
||||
plotter->Fill2D("IntersectionPhi_vs_AnodeZ_nC", 400, -200, 200, 600, -300, 300, anodeIntersection.Phi() * 180. / TMath::Pi(), anodeIntersection.Z(), "hGMPC");
|
||||
}
|
||||
if (anodeHits.size() > 0 && cathodeHits.size() > 0)
|
||||
plotter->Fill2D("AHits_vs_CHits", 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
|
||||
|
||||
// make another plot with nearest neighbour constraint
|
||||
bool hasNeighbourAnodes = false;
|
||||
bool hasNeighbourCathodes = false;
|
||||
|
||||
// 1. Check Anodes for neighbours (including wrap-around 0-23)
|
||||
for (size_t i = 0; i < anodeHits.size(); i++)
|
||||
{
|
||||
for (size_t j = i + 1; j < anodeHits.size(); j++)
|
||||
{
|
||||
int diff = std::abs(anodeHits[i].first - anodeHits[j].first);
|
||||
if (diff == 1 || diff == 23)
|
||||
{ // 23 handles the cylindrical wrap
|
||||
hasNeighbourAnodes = true;
|
||||
break;
|
||||
}
|
||||
}
|
||||
if (hasNeighbourAnodes)
|
||||
break;
|
||||
}
|
||||
|
||||
// 2. Check Cathodes for neighbours (including wrap-around 0-23)
|
||||
for (size_t i = 0; i < cathodeHits.size(); i++)
|
||||
{
|
||||
for (size_t j = i + 1; j < cathodeHits.size(); j++)
|
||||
{
|
||||
int diff = std::abs(cathodeHits[i].first - cathodeHits[j].first);
|
||||
if (diff == 1 || diff == 23)
|
||||
{
|
||||
hasNeighbourCathodes = true;
|
||||
break;
|
||||
}
|
||||
}
|
||||
if (hasNeighbourCathodes)
|
||||
break;
|
||||
}
|
||||
|
||||
// ---------------------------------------------------------
|
||||
// FILL PLOTS
|
||||
// ---------------------------------------------------------
|
||||
if (anodeHits.size() > 0 && cathodeHits.size() > 0)
|
||||
{
|
||||
plotter->Fill2D("AHits_vs_CHits_NA" + std::to_string(hasNeighbourAnodes), 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
|
||||
plotter->Fill2D("AHits_vs_CHits_NC" + std::to_string(hasNeighbourCathodes), 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
|
||||
|
||||
// Constraint Plot: Only fill if BOTH planes have adjacent hits
|
||||
// This effectively removes events with only isolated single-wire hits (noise)
|
||||
if (hasNeighbourAnodes && hasNeighbourCathodes)
|
||||
{
|
||||
plotter->Fill2D("AHits_vs_CHits_NN", 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
|
||||
}
|
||||
}
|
||||
|
||||
if (HitNonZero && anodeIntersection.Z() != 0)
|
||||
{
|
||||
pw_contr.CalTrack2(hitPos, anodeIntersection);
|
||||
plotter->Fill1D("VertexRecon", 600, -1300, 1300, pw_contr.GetZ0());
|
||||
plotter->Fill1D("VertexRecon_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 600, -1300, 1300, pw_contr.GetZ0());
|
||||
|
||||
if (cathodeHits.size() == 2)
|
||||
plotter->Fill1D("VertexRecon_2c_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 600, -1300, 1300, pw_contr.GetZ0());
|
||||
|
||||
TVector3 x2(anodeIntersection), x1(hitPos);
|
||||
|
||||
TVector3 v = x2-x1;
|
||||
double t_minimum = -1.0*(x1.X()*v.X()+x1.Y()*v.Y())/(v.X()*v.X()+v.Y()*v.Y());
|
||||
vector_closest_to_z = x1 + t_minimum*v;
|
||||
|
||||
plotter->Fill1D("VertexRecon_Z_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 600, -1300, 1300, vector_closest_to_z.Z() ,"customVertex");
|
||||
if(vector_closest_to_z.Perp() < 20) {
|
||||
plotter->Fill1D("VertexRecon_RadialCut_Z_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 600, -1300, 1300, vector_closest_to_z.Z() ,"customVertex");
|
||||
}
|
||||
|
||||
plotter->Fill2D("VertexRecon_XY_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 100, -100, 100, 100,-100,100, vector_closest_to_z.X(), vector_closest_to_z.Y() ,"customVertex");
|
||||
if(cathodeHits.size()==2) {
|
||||
plotter->Fill1D("VertexRecon2C_Z_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 600, -1300, 1300, vector_closest_to_z.Z() ,"customVertex");
|
||||
if(vector_closest_to_z.Perp() < 20) {
|
||||
plotter->Fill1D("VertexRecon2C_RadialCut_Z_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 600, -1300, 1300, vector_closest_to_z.Z() ,"customVertex");
|
||||
}
|
||||
plotter->Fill2D("VertexRecon2C_XY_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 100, -100, 100, 100,-100,100, vector_closest_to_z.X(), vector_closest_to_z.Y() ,"customVertex");
|
||||
plotter->Fill2D("VertexRecon2C_RhoZ_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 100, -100, 100, 600,-1300,1300, vector_closest_to_z.Perp(), vector_closest_to_z.Z() ,"customVertex");
|
||||
plotter->Fill2D("VertexRecon2C_Z_vs_QQQE_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 600, -1300, 1300, 800,0,20000, vector_closest_to_z.Z(), qqqenergy ,"customVertex");
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
for (int i = 0; i < qqq.multi; i++)
|
||||
{
|
||||
if(anodeIntersection.Perp() > 0) { //suppress x,y=0,0 events
|
||||
if (PCQQQTimeCut) {
|
||||
plotter->Fill2D("PC_XY_Projection_QQQ_TimeCut" + std::to_string(qqq.id[i]), 400, -100, 100, 400, -100, 100, anodeIntersection.X(), anodeIntersection.Y(), "hPCQQQ");
|
||||
}
|
||||
plotter->Fill2D("PC_XY_Projection_QQQ" + std::to_string(qqq.id[i]), 400, -100, 100, 400, -100, 100, anodeIntersection.X(), anodeIntersection.Y(), "hPCQQQ");
|
||||
}
|
||||
for (int j = i + 1; j < qqq.multi; j++)
|
||||
{
|
||||
if (qqq.id[i] == qqq.id[j])
|
||||
{
|
||||
int chWedge = -1;
|
||||
int chRing = -1;
|
||||
double eWedge = 0.0;
|
||||
double eWedgeMeV = 0.0;
|
||||
double eRing = 0.0;
|
||||
double eRingMeV = 0.0;
|
||||
double tRing = 0.0;
|
||||
int qqqID = -1;
|
||||
if (qqq.ch[i] < 16 && qqq.ch[j] >= 16 && qqqGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16])
|
||||
{
|
||||
chWedge = qqq.ch[i];
|
||||
eWedge = qqq.e[i] * qqqGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16];
|
||||
chRing = qqq.ch[j] - 16;
|
||||
eRing = qqq.e[j];
|
||||
tRing = static_cast<double>(qqq.t[j]);
|
||||
qqqID = qqq.id[i];
|
||||
}
|
||||
else if (qqq.ch[j] < 16 && qqq.ch[i] >= 16 && qqqGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16])
|
||||
{
|
||||
chWedge = qqq.ch[j];
|
||||
eWedge = qqq.e[j] * qqqGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16];
|
||||
chRing = qqq.ch[i] - 16;
|
||||
tRing = static_cast<double>(qqq.t[i]);
|
||||
eRing = qqq.e[i];
|
||||
qqqID = qqq.id[i];
|
||||
}
|
||||
else
|
||||
continue;
|
||||
|
||||
if (qqqCalibValid[qqq.id[i]][chRing][chWedge])
|
||||
{
|
||||
eWedgeMeV = eWedge * qqqCalib[qqq.id[i]][chRing][chWedge] / 1000;
|
||||
eRingMeV = eRing * qqqCalib[qqq.id[i]][chRing][chWedge] / 1000;
|
||||
}
|
||||
else
|
||||
continue;
|
||||
|
||||
// if (anodeIntersection.Z() != 0)
|
||||
{
|
||||
plotter->Fill2D("PC_Z_vs_QQQRing", 600, -300, 300, 16, 0, 16, anodeIntersection.Z(), chRing, "hPCzQQQ");
|
||||
}
|
||||
|
||||
if (anodeIntersection.Z() != 0 && cathodeHits.size() == 2)
|
||||
{
|
||||
plotter->Fill2D("PC_Z_vs_QQQRing_2C", 600, -300, 300, 16, 0, 16, anodeIntersection.Z(), chRing, "hPCzQQQ");
|
||||
plotter->Fill2D("PC_Z_vs_QQQRing_2C" + std::to_string(qqq.id[i]), 600, -300, 300, 16, 0, 16, anodeIntersection.Z(), chRing, "hPCzQQQ");
|
||||
plotter->Fill2D("PC_Z_vs_QQQWedge_2C", 600, -300, 300, 16, 0, 16, anodeIntersection.Z(), chWedge, "hPCzQQQ");
|
||||
}
|
||||
plotter->Fill2D("VertexRecon_QQQRingTC" + std::to_string(PCQQQTimeCut) + "PhiC" + std::to_string(PCQQQPhiCut), 600, -1300, 1300, 16, 0, 16, vector_closest_to_z.Z(), chRing, "hPCQQQ");
|
||||
double phi = TMath::ATan2(anodeIntersection.Y(), anodeIntersection.X()) * 180. / TMath::Pi();
|
||||
plotter->Fill2D("PolarAngle_Vs_QQQWedge" + std::to_string(qqqID), 360, -360, 360, 16, 0, 16, phi, chWedge, "hPCQQQ");
|
||||
// plotter->Fill2D("EdE_PC_vs_QQQ_timegate_ls1000"+std::to_string())
|
||||
|
||||
plotter->Fill2D("PC_Z_vs_QQQRing_Det" + std::to_string(qqqID), 600, -300, 300, 16, 0, 16, anodeIntersection.Z(), chRing, "hPCQQQ");
|
||||
//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);
|
||||
|
||||
for (int k = 0; k < pc.multi; k++)
|
||||
{
|
||||
if(pc.index[k] >= 24)
|
||||
continue;
|
||||
|
||||
// double sinTheta = TMath::Sin((hitPos-vector_closest_to_z).Theta());
|
||||
double sinTheta = TMath::Sin((anodeIntersection-TVector3(0,0,90.0)).Theta());
|
||||
// double sinTheta = TMath::Sin((hitPos-TVector3(0,0,30.0)).Theta());
|
||||
// double sinTheta = TMath::Sin(hitPos.Theta());
|
||||
|
||||
if(cathodeHits.size()==2 && PCQQQPhiCut) {
|
||||
plotter->Fill2D("CalibratedQQQE_RvsCPCE_TC" + std::to_string(PCQQQTimeCut) , 400, 0, 10, 400, 0, 30000, eRingMeV, pc.e[k]*sinTheta, "hPCQQQ");
|
||||
plotter->Fill2D("CalibratedQQQE_WvsCPCE_TC" + std::to_string(PCQQQTimeCut) , 400, 0, 10, 400, 0, 30000, eWedgeMeV, pc.e[k]*sinTheta, "hPCQQQ");
|
||||
plotter->Fill2D("CalibratedQQQE_RvsPCE_TC" + std::to_string(PCQQQTimeCut) , 400, 0, 10, 400, 0, 30000, eRingMeV, pc.e[k], "hPCQQQ");
|
||||
plotter->Fill2D("CalibratedQQQE_WvsPCE_TC" + std::to_string(PCQQQTimeCut) , 400, 0, 10, 400, 0, 30000, eWedgeMeV, pc.e[k], "hPCQQQ");
|
||||
plotter->Fill2D("PCQQQ_dTimevsdPhi", 200, -2000, 2000, 80, -200, 200, tRing - static_cast<double>(pc.t[k]), (hitPos.Phi()-anodeIntersection.Phi()) * 180. / TMath::Pi(), "hTiming");
|
||||
}
|
||||
|
||||
}
|
||||
}///qqq i==j case end
|
||||
} //j loop end
|
||||
} // qqq i loop end
|
||||
TVector3 guessVertex(0,0,90.); //for run12, subtract anodeIntersection.Z() by ~74.0 seems to work
|
||||
//rho=40.0 mm is halfway between the cathodes(rho=42) and anodes(rho=37)
|
||||
double pcz_guess = 42.0/TMath::Tan((hitPos-guessVertex).Theta()) + guessVertex.Z(); //this is ideally kept to be all QQQ+userinput for calibration of pcz
|
||||
if(PCQQQTimeCut && PCQQQPhiCut && hitPos.Perp()>0 && anodeIntersection.Perp()>0 && cathodeHits.size()>=2) {
|
||||
plotter->Fill2D("pczguess_vs_qqqE_rad="+std::to_string(hitPos.Perp()),100,0,200,800,0,20,pcz_guess,qqqenergy,"pczguess");
|
||||
plotter->Fill2D("pczguess_vs_pcz_rad="+std::to_string(hitPos.Perp()),100,0,200,150,0,200,pcz_guess,anodeIntersection.Z(),"pczguess"); //entirely qqq-derived position vs entirely PC derived position
|
||||
plotter->Fill2D("pczguess_vs_pcz_phi="+std::to_string(hitPos.Phi()*180./M_PI),100,0,200,150,0,200,pcz_guess,anodeIntersection.Z()/0.8,"pczguess"); //entirely qqq-derived position vs entirely PC derived position
|
||||
plotter->Fill2D("pczguess_vs_pcz",100,0,200,150,0,200,pcz_guess,anodeIntersection.Z()/0.8);
|
||||
plotter->Fill2D("pcz_vs_pcPhi_rad="+std::to_string(hitPos.Perp()),360,0,360,150,0,200,anodeIntersection.Phi()*180./M_PI,anodeIntersection.Z(),"pczguess");
|
||||
}
|
||||
for (int i = 0; i < sx3.multi; i++)
|
||||
{
|
||||
// plotting sx3 strip hits vs anode phi
|
||||
if (sx3.ch[i] < 8 && anodeIntersection.Perp()>0)
|
||||
plotter->Fill2D("PCPhi_vs_SX3Strip", 100, -200, 200, 8 * 24, 0, 8 * 24, anodeIntersection.Phi() * 180. / TMath::Pi(), sx3.id[i] * 8 + sx3.ch[i]);
|
||||
}
|
||||
|
||||
if (anodeIntersection.Z() != 0 && cathodeHits.size() == 3)
|
||||
{
|
||||
plotter->Fill1D("PC_Z_proj_3C", 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
|
||||
}
|
||||
|
||||
plotter->Fill2D("AnodeMaxE_Vs_Cathode_Sum_Energy", 2000, 0, 30000, 2000, 0, 30000, aEMax, cESum, "hGMPC");
|
||||
plotter->Fill1D("Correlated_Cathode_MaxAnode", 6, 0, 5, corrcatMax.size(), "hGMPC");
|
||||
plotter->Fill2D("Correlated_Cathode_VS_MaxAnodeEnergy", 6, 0, 5, 2000, 0, 30000, corrcatMax.size(), aEMax, "hGMPC");
|
||||
plotter->Fill1D("AnodeHits", 12, 0, 11, anodeHits.size(), "hGMPC");
|
||||
plotter->Fill2D("AnodeMaxE_vs_AnodeHits", 12, 0, 11, 2000, 0, 30000, anodeHits.size(), aEMax, "hGMPC");
|
||||
|
||||
if (anodeHits.size() < 1)
|
||||
{
|
||||
plotter->Fill1D("NoAnodeHits_CathodeHits", 6, 0, 5, cathodeHits.size(), "hGMPC");
|
||||
}
|
||||
|
||||
return kTRUE;
|
||||
}
|
||||
|
||||
void MakeVertex::Terminate()
|
||||
{
|
||||
plotter->FlushToDisk();
|
||||
}
|
||||
BIN
MakeVertex.C:Zone.Identifier
Normal file
BIN
MakeVertex.C:Zone.Identifier
Normal file
Binary file not shown.
132
MakeVertex.h
Normal file
132
MakeVertex.h
Normal file
|
|
@ -0,0 +1,132 @@
|
|||
#ifndef MakeVertex_h
|
||||
#define MakeVertex_h
|
||||
|
||||
#include <TROOT.h>
|
||||
#include <TChain.h>
|
||||
#include <TApplication.h>
|
||||
#include <TFile.h>
|
||||
#include <TSelector.h>
|
||||
#include <iomanip>
|
||||
#include <vector> // Required for vectors
|
||||
#include <utility> // Required for std::pair
|
||||
|
||||
#include "Armory/ClassDet.h"
|
||||
#include "Armory/ClassPW.h" // YOU ADDED THIS (Correct! Defines Coord)
|
||||
|
||||
class MakeVertex : public TSelector {
|
||||
public :
|
||||
TTree *fChain; //!pointer to the analyzed TTree or TChain
|
||||
|
||||
// Declaration of leaf types
|
||||
Det sx3;
|
||||
Det qqq;
|
||||
Det pc ;
|
||||
Det misc;
|
||||
|
||||
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; //!
|
||||
TBranch *b_miscMulti; //!
|
||||
TBranch *b_miscID; //!
|
||||
TBranch *b_miscCh; //!
|
||||
TBranch *b_miscE; //!
|
||||
TBranch *b_miscT; //!
|
||||
TBranch *b_miscTf; //!
|
||||
|
||||
// 1. Geometry Cache
|
||||
Coord Crossover[24][24][2];
|
||||
|
||||
// 2. Persistent Vectors (REQUIRED for the optimized .cxx to work)
|
||||
std::vector<std::pair<int, double>> anodeHits;
|
||||
std::vector<std::pair<int, double>> cathodeHits;
|
||||
std::vector<std::pair<int, double>> corrcatMax;
|
||||
std::vector<std::pair<int, double>> corranoMax;
|
||||
std::vector<double> cathodeTimes;
|
||||
std::vector<double> anodeTimes;
|
||||
|
||||
MakeVertex(TTree * /*tree*/ =0) : fChain(0) { }
|
||||
virtual ~MakeVertex() { }
|
||||
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(MakeVertex,0);
|
||||
};
|
||||
|
||||
#endif
|
||||
|
||||
#ifdef MakeVertex_cxx
|
||||
void MakeVertex::Init(TTree *tree){
|
||||
|
||||
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);
|
||||
fChain->SetBranchAddress("miscMulti", &misc.multi, &b_miscMulti);
|
||||
fChain->SetBranchAddress("miscID", &misc.id, &b_miscID);
|
||||
fChain->SetBranchAddress("miscCh", &misc.ch, &b_miscCh);
|
||||
fChain->SetBranchAddress("miscE", &misc.e, &b_miscE);
|
||||
fChain->SetBranchAddress("miscT", &misc.t, &b_miscT);
|
||||
}
|
||||
|
||||
Bool_t MakeVertex::Notify(){
|
||||
return kTRUE;
|
||||
}
|
||||
|
||||
void MakeVertex::SlaveBegin(TTree * /*tree*/){
|
||||
// TString option = GetOption();
|
||||
}
|
||||
|
||||
void MakeVertex::SlaveTerminate(){
|
||||
}
|
||||
#endif // #ifdef MakeVertex_cxx
|
||||
BIN
MakeVertex.h:Zone.Identifier
Normal file
BIN
MakeVertex.h:Zone.Identifier
Normal file
Binary file not shown.
34
TrackRecon.C
34
TrackRecon.C
|
|
@ -21,6 +21,7 @@
|
|||
PW pw_contr;
|
||||
PW pwinstance;
|
||||
TVector3 hitPos;
|
||||
double qqqenergy,qqqtime;
|
||||
|
||||
// Calibration globals
|
||||
const int MAX_QQQ = 4;
|
||||
|
|
@ -206,6 +207,9 @@ Bool_t TrackRecon::Process(Long64_t entry)
|
|||
|
||||
int qqqCount = 0;
|
||||
int qqqAdjCh = 0;
|
||||
qqqenergy=-1;
|
||||
qqqtime=-1;
|
||||
|
||||
// REMOVE WHEN RERUNNING USING THE NEW CALIBRATION FILE
|
||||
// for (int i = 0; i < qqq.multi; i++)
|
||||
// {
|
||||
|
|
@ -355,6 +359,8 @@ Bool_t TrackRecon::Process(Long64_t entry)
|
|||
double x = rho * TMath::Cos(theta);
|
||||
double y = rho * TMath::Sin(theta);
|
||||
hitPos.SetXYZ(x, y, 23 + 75 + 30);
|
||||
qqqenergy = eRingMeV;
|
||||
qqqtime = tRing;
|
||||
HitNonZero = true;
|
||||
}
|
||||
}
|
||||
|
|
@ -397,7 +403,7 @@ Bool_t TrackRecon::Process(Long64_t entry)
|
|||
for (int j = 0; j < qqq.multi; j++)
|
||||
{
|
||||
plotter->Fill1D("AC_Time_qqq_coin", 500, -2500, 2500, anodeT - cathodeT, "hTiming");
|
||||
plotter->Fill2D("AC_Time_Vs_QQQ_ch", 500, -2500, 2500, 16 * 8, 0, 16 * 8, anodeT - cathodeT, qqq.ch[j], "hTiming");
|
||||
// plotter->Fill2D("AC_Time_Vs_QQQ_ch", 500, -2500, 2500, 16 * 8, 0, 16 * 8, anodeT - cathodeT, qqq.ch[j], "hTiming");
|
||||
plotter->Fill2D("AC_Time_vs_AIndex", 500, -2500, 2500, 24, 0, 24, anodeT - cathodeT, anodeIndex, "hTiming");
|
||||
plotter->Fill2D("AC_Time_vs_CIndex", 500, -2500, 2500, 24, 0, 24, anodeT - cathodeT, cathodeIndex, "hTiming");
|
||||
// plotter->Fill1D("AC_Time_A" + std::to_string(anodeIndex) + "_C" + std::to_string(cathodeIndex), 200, -1000, 1000, anodeT - cathodeT, "TimingPC");
|
||||
|
|
@ -413,8 +419,7 @@ Bool_t TrackRecon::Process(Long64_t entry)
|
|||
|
||||
for (int j = i + 1; j < pc.multi; j++)
|
||||
{
|
||||
plotter->Fill2D("PC_Coincidence_Matrix", 48, 0, 48, 48, 0, 48, pc.index[i], pc.index[j], "hRawPC");
|
||||
plotter->Fill2D("PC_Coincidence_Matrix_anodeMinusCathode_lt_-200_" + std::to_string(anodeT - cathodeT < -200), 48, 0, 48, 48, 0, 48, pc.index[i], pc.index[j], "hRawPC");
|
||||
// plotter->Fill2D("PC_Coincidence_Matrix_anodeMinusCathode_lt_-200_" + std::to_string(anodeT - cathodeT < -200), 48, 0, 48, 48, 0, 48, pc.index[i], pc.index[j], "hRawPC");
|
||||
plotter->Fill2D("Anode_V_Anode", 24, 0, 24, 24, 0, 24, pc.index[i], pc.index[j], "hGMPC");
|
||||
}
|
||||
}
|
||||
|
|
@ -475,7 +480,7 @@ Bool_t TrackRecon::Process(Long64_t entry)
|
|||
plotter->Fill2D("AnodeMax_Vs_Cathode_Coincidence_Matrix", 24, 0, 24, 24, 0, 24, aIDMax, cID, "hRawPC");
|
||||
plotter->Fill2D("Anode_Vs_Cathode_Coincidence_Matrix", 24, 0, 24, 24, 0, 24, aID, cID, "hRawPC");
|
||||
plotter->Fill2D("Anode_vs_CathodeE", 2000, 0, 30000, 2000, 0, 30000, aE, cE, "hGMPC");
|
||||
plotter->Fill2D("CathodeMult_V_CathodeE", 6, 0, 6, 2000, 0, 30000, cathodeHits.size(), cE, "hGMPC");
|
||||
// plotter->Fill2D("CathodeMult_V_CathodeE", 6, 0, 6, 2000, 0, 30000, cathodeHits.size(), cE, "hGMPC");
|
||||
for (int j = -4; j < 3; j++)
|
||||
{
|
||||
if ((aIDMax + 24 + j) % 24 == 23 - cID)
|
||||
|
|
@ -550,12 +555,12 @@ Bool_t TrackRecon::Process(Long64_t entry)
|
|||
if (anodeIntersection.Perp() != 0 && cathodeHits.size() == 2)
|
||||
{
|
||||
plotter->Fill1D("PC_Z_proj_2C", 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
|
||||
plotter->Fill2D("IntersectionPhi_vs_AnodeZ_2C", 400, -200, 200, 600, -300, 300, anodeIntersection.Phi() * 180. / TMath::Pi(), anodeIntersection.Z(), "hGMPC");
|
||||
plotter->Fill2D("IntersectionPhi_vs_AnodeZ_2C", 400, -200, 200, 600, -300, 300, anodeIntersection.Phi() * 180. / TMath::Pi(), anodeIntersection.Z(), "hPCzQQQ");
|
||||
}
|
||||
if (anodeIntersection.Perp() != 0 && cathodeHits.size() > 2)
|
||||
{
|
||||
plotter->Fill1D("PC_Z_proj_nC", 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
|
||||
plotter->Fill2D("IntersectionPhi_vs_AnodeZ_nC", 400, -200, 200, 600, -300, 300, anodeIntersection.Phi() * 180. / TMath::Pi(), anodeIntersection.Z(), "hGMPC");
|
||||
plotter->Fill2D("IntersectionPhi_vs_AnodeZ_nC", 400, -200, 200, 600, -300, 300, anodeIntersection.Phi() * 180. / TMath::Pi(), anodeIntersection.Z(), "hPCzQQQ");
|
||||
}
|
||||
if (anodeHits.size() > 0 && cathodeHits.size() > 0)
|
||||
plotter->Fill2D("AHits_vs_CHits", 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
|
||||
|
|
@ -601,8 +606,8 @@ Bool_t TrackRecon::Process(Long64_t entry)
|
|||
// ---------------------------------------------------------
|
||||
if (anodeHits.size() > 0 && cathodeHits.size() > 0)
|
||||
{
|
||||
plotter->Fill2D("AHits_vs_CHits_NA" + std::to_string(hasNeighbourAnodes), 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
|
||||
plotter->Fill2D("AHits_vs_CHits_NC" + std::to_string(hasNeighbourCathodes), 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
|
||||
// plotter->Fill2D("AHits_vs_CHits_NA" + std::to_string(hasNeighbourAnodes), 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
|
||||
// plotter->Fill2D("AHits_vs_CHits_NC" + std::to_string(hasNeighbourCathodes), 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
|
||||
|
||||
// Constraint Plot: Only fill if BOTH planes have adjacent hits
|
||||
// This effectively removes events with only isolated single-wire hits (noise)
|
||||
|
|
@ -707,6 +712,19 @@ Bool_t TrackRecon::Process(Long64_t entry)
|
|||
}
|
||||
}
|
||||
}
|
||||
|
||||
TVector3 guessVertex(0, 0, 90.);
|
||||
// rho=40.0 mm is halfway between the cathodes(rho=42) and anodes(rho=37)
|
||||
double pcz_guess = 42.0 / TMath::Tan((hitPos - guessVertex).Theta()) + guessVertex.Z(); // this is ideally kept to be all QQQ+userinput for calibration of pcz
|
||||
if (PCQQQTimeCut && PCQQQPhiCut && hitPos.Perp() > 0 && anodeIntersection.Perp() > 0 && cathodeHits.size() >= 2)
|
||||
{
|
||||
plotter->Fill2D("pczguess_vs_qqqE", 100, 0, 200, 800, 0, 20, pcz_guess, qqqenergy, "pczguess");
|
||||
// plotter->Fill2D("pczguess_vs_pcz_rad="+std::to_string(hitPos.Perp()),100,0,200,150,0,200,pcz_guess,anodeIntersection.Z(),"pczguess"); //entirely qqq-derived position vs entirely PC derived position
|
||||
plotter->Fill2D("pczguess_vs_pcz_phi=" + std::to_string(hitPos.Phi() * 180. / M_PI), 100, 0, 200, 150, 0, 200, pcz_guess, anodeIntersection.Z() / 0.8, "pczguess"); // entirely qqq-derived position vs entirely PC derived position
|
||||
plotter->Fill2D("pczguess_vs_pcz", 100, 0, 200, 150, 0, 200, pcz_guess, anodeIntersection.Z() / 0.8);
|
||||
plotter->Fill2D("pcz_vs_pcPhi_rad=" + std::to_string(hitPos.Perp()), 360, 0, 360, 150, 0, 200, anodeIntersection.Phi() * 180. / M_PI, anodeIntersection.Z() / 0.8, "pczguess");
|
||||
}
|
||||
|
||||
for (int i = 0; i < sx3.multi; i++)
|
||||
{
|
||||
// plotting sx3 strip hits vs anode phi
|
||||
|
|
|
|||
Loading…
Reference in New Issue
Block a user