Skip to content

Commit

Permalink
Shower-max PE response lookup table analysis scripts and data files a…
Browse files Browse the repository at this point in the history
…dded.
  • Loading branch information
sudipbhattarai committed Jan 22, 2025
1 parent 6a5212d commit 3fab176
Show file tree
Hide file tree
Showing 8 changed files with 416 additions and 0 deletions.
18 changes: 18 additions & 0 deletions analysis/showermax/data/fit_param_xy_e-_ifarm.csv
Original file line number Diff line number Diff line change
@@ -0,0 +1,18 @@
energy,xp0,xp1,yp0,yp1,yp2
5,0.00186854,-6.27248e-07,0.000318276,3.09224e-07,1.11394e-07
10,0.0102587,-1.50709e-05,0.0182153,2.54559e-05,-1.01119e-06
50,0.617914,0.00109035,0.734966,-0.000662327,-2.13036e-05
100,2.40497,0.00304654,2.52607,-0.000402899,-2.82908e-05
500,17.2653,0.0150206,17.7501,-0.00199242,-0.000123669
1000,35.8243,0.0354252,36.7149,-0.00443617,-0.000235146
2000,71.2615,0.0664176,73.2177,-0.0101103,-0.000481341
1000,35.8243,0.0354252,36.7149,-0.00443617,-0.000235146
3000,105.66,0.102495,107.507,-0.00432239,-0.000539539
4000,138.983,0.129233,141.393,-0.00765265,-0.000714903
5000,170.902,0.156594,174.967,-0.0165793,-0.00102413
6000,202.591,0.195464,206.781,-0.0185614,-0.00107666
7000,233.861,0.214678,238.794,-0.0255901,-0.00132374
8000,265.736,0.253171,270.522,-0.021348,-0.001311
9000,295.394,0.2864,301.713,-0.0221357,-0.00161815
10000,324.714,0.318387,331.45,-0.0336423,-0.00174544
11000,354.866,0.339425,362.352,-0.00987083,-0.0018857
18 changes: 18 additions & 0 deletions analysis/showermax/data/fit_param_xy_gamma_ifarm.csv
Original file line number Diff line number Diff line change
@@ -0,0 +1,18 @@
energy,xp0,xp1,yp0,yp1,yp2
5,0.0627517,-1.0543e-06,0.0581547,-0.000192039,1.40549e-06
10,0.148161,-0.000245255,0.151718,-0.00012809,-9.33095e-07
50,1.42496,0.00174349,1.46062,0.000384533,-1.0486e-05
100,3.20275,0.00199171,3.36658,-0.000972939,-3.72528e-05
500,17.7877,0.0168821,18.5357,-0.00305322,-0.000173396
1000,35.7328,0.0391455,37.0105,-0.00385948,-0.000290996
2000,68.9893,0.0654994,70.7129,-0.00186141,-0.000405991
1000,35.7328,0.0391455,37.0105,-0.00385948,-0.000290996
3000,100.729,0.0891212,103.059,-0.00785097,-0.00059099
4000,131.482,0.121328,134.431,-0.00446825,-0.000761432
5000,161.28,0.161656,165.487,-0.0160822,-0.000957194
6000,189.225,0.163666,193.413,-0.0122421,-0.000973039
7000,218.903,0.222023,223.881,-0.0321237,-0.00125839
8000,247.765,0.239203,252.623,-0.0280649,-0.00124614
9000,273.899,0.262618,279.212,0.00200576,-0.0014183
10000,302.462,0.304135,309.125,-0.0166459,-0.00164194
11000,329.102,0.376249,338.584,-0.0376803,-0.00227799
14 changes: 14 additions & 0 deletions analysis/showermax/data/fit_param_xy_mu-_ifarm.csv
Original file line number Diff line number Diff line change
@@ -0,0 +1,14 @@
energy,xp0,xp1,yp0,yp1,yp2
500,7.61904,0.00705542,7.68358,8.26442e-05,-2.35436e-05
1000,9.84448,0.00476633,9.85678,-0.000641187,-7.83977e-06
2000,10.6694,0.00424422,10.6791,0.000134362,-1.42913e-06
1000,9.84448,0.00476633,9.85678,-0.000641187,-7.83977e-06
3000,10.902,0.00580823,10.7998,-0.000183595,1.279e-05
4000,10.9523,0.00511675,10.9356,-0.000334526,-9.12267e-07
5000,11.1126,0.00453355,11.0124,8.11234e-05,1.29421e-05
6000,11.0563,0.00719545,11.0066,0.000110291,4.91498e-06
7000,11.3055,0.00669282,11.1802,0.000123959,1.57476e-05
8000,11.1125,0.00717593,10.9525,-0.000299558,2.83419e-05
9000,11.2697,0.00621689,11.1181,-0.000846437,2.12847e-05
10000,11.2074,0.00588017,11.126,-0.0010563,1.76701e-05
11000,11.2681,0.00494686,11.1262,0.000572583,2.20093e-05
13 changes: 13 additions & 0 deletions analysis/showermax/data/fit_param_xy_neutron_ifarm.csv
Original file line number Diff line number Diff line change
@@ -0,0 +1,13 @@
energy,xp0,xp1,yp0,yp1,yp2
1000,0.0325032,5.81897e-05,0.0380673,9.18402e-06,-7.12998e-07
2000,0.664703,0.000268984,0.720299,-7.24534e-05,-1.19415e-05
1000,0.0325032,5.81897e-05,0.0380673,9.18402e-06,-7.12998e-07
3000,1.83493,0.000569759,2.03071,-0.000872436,-3.77283e-05
4000,3.25583,-2.87781e-05,3.35537,6.4597e-05,-2.61931e-05
5000,4.34997,0.00252794,4.64361,-0.00141612,-5.65787e-05
6000,5.83171,0.00972787,5.98878,-0.00104332,-4.09446e-05
7000,7.62846,0.00426204,8.04478,-0.000292948,-8.13441e-05
8000,9.8956,0.00675217,10.3159,0.00146602,-9.63828e-05
9000,11.0889,0.00751149,11.8324,-0.00388285,-0.000106202
10000,13.8775,0.00984104,14.9472,0.00143721,-0.000206724
11000,15.6238,0.0224891,16.3091,-0.00900119,-0.000160528
14 changes: 14 additions & 0 deletions analysis/showermax/data/fit_param_xy_pi-_ifarm.csv
Original file line number Diff line number Diff line change
@@ -0,0 +1,14 @@
energy,xp0,xp1,yp0,yp1,yp2
500,5.39623,0.00672054,5.5963,-0.000690162,-4.39051e-05
1000,8.70975,0.00426233,8.69977,0.000156259,-7.83212e-06
2000,10.7281,0.00291869,10.8914,0.000375828,-4.35732e-05
1000,8.70975,0.00426233,8.69977,0.000156259,-7.83212e-06
3000,12.4285,0.00493266,12.6518,0.000701055,-4.96781e-05
4000,14.5609,-0.000216201,15.0292,-0.00106373,-0.000120579
5000,17.187,0.00532813,17.9249,0.000601087,-0.000154178
6000,21.536,0.00222366,21.9014,-0.00511615,-0.000135389
7000,23.5107,0.015075,23.7392,-0.00262259,-6.06112e-05
8000,25.8184,0.015339,26.8405,-0.00330766,-0.000222387
9000,28.3697,0.0249847,28.6818,0.0040931,-0.000111408
10000,29.5291,0.0173609,29.3689,-0.00523963,-4.497e-05
11000,30.8697,0.0487093,31.7276,0.00668593,-0.000143964
126 changes: 126 additions & 0 deletions analysis/showermax/scripts/plot_showermax_response.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,126 @@
// ROOTscript: plot_shower-max_response.C
// -----------------------------------------------------------------------------------
// Author: Sudip Bhattarai
// Date: 01/21/2025
// Purpose: Read the remoll rootfile and use shower-max lookup table to plot the PE response.
// -----------------------------------------------------------------------------------

#include "TChain.h"
#include "TH1.h"
#include "TROOT.h"
#include "TCanvas.h"
#include <iostream>
#include <fstream>
#include <utility>
#include <vector>
#include <string>
#include "./remollToQsim.hh"
#include "./shower-max_resposne_lookup.hh"

// Namespace list to use
using std::cout, std::cerr, std::endl;
using std::vector, std::string;
using std::ifstream, std::ofstream;
using std::pair, std::map;

// Main function
void plot_showermax_response(){
// Genetate random values for x, y, theta, phi

// Get the fit parameters
vector<vector<Double_t>> fit_data_electron = retrieve_fit_data("e-");
vector<vector<Double_t>> fit_data_gamma = retrieve_fit_data("gamma");
vector<vector<Double_t>> fit_data_mu = retrieve_fit_data("mu-");
vector<vector<Double_t>> fit_data_pi = retrieve_fit_data("pi-");
vector<vector<Double_t>> fit_data_neutron = retrieve_fit_data("neutron");

// Make a map of pid and fit_data
map<int, vector<vector<Double_t>>> map_fit_data;
map_fit_data[11] = fit_data_electron;
map_fit_data[-11] = fit_data_electron;
map_fit_data[22] = fit_data_gamma;
map_fit_data[13] = fit_data_mu;
map_fit_data[-13] = fit_data_mu;
map_fit_data[211] = fit_data_pi;
map_fit_data[-211] = fit_data_pi;
map_fit_data[2112] = fit_data_neutron;

// Test pe response
// cout << "PE: " << get_PE_response(map_fit_data[11], 1234, 0, 0) << endl;

// Remoll file path goes here
int goodFileCount = 1; // Number of good (non-corrupted) files
TString rootFile= "~/moller/softwares/remoll-zone/rootfiles/smStack_v23/sm_tqStack_3sector_moller_50k_1001.root";

// Define histograms
TH1D* h_hitRate = new TH1D("hitRate", "Rate weighted hit; hit.r [mm]; rate[GHz/65uA]", 100, 1000, 1200);
TH1D* h_hitRateSmPE = new TH1D("hitRateSmPE", "Rate weighted hit; hit.r [mm]; rate [GHz/65uA/PE]", 100, 1000, 1200);

// Declare TChain
TChain* T = new TChain("T");
T->Add(rootFile);

Long64_t nEntries = T->GetEntries();
cout << "Total number of entries in the chain: " << nEntries << endl;
std::vector<remollGenericDetectorHit_t> *fHit = 0;
// remollEvent_t *fEv = 0;
Double_t fRate = 0;
Float_t energy(-1.0e-12), rate(-1.0e-12),
hitx(-1.0e-12), hity(-1.0e-12), hitr(-1.0e-12), hitpz(-1.0e-12),
pe(0);
Int_t pid(0), det(0);

T->SetBranchAddress("hit", &fHit);
T->SetBranchAddress("rate", &fRate);
// T->SetBranchAddress("ev", &fEv);

//This loop goes over all the events in the root files
for (Long64_t iEvent = 0; iEvent < nEntries; iEvent++){
if (iEvent % (nEntries/10) == 0)
cout << "Analyzed " << iEvent << " events." << endl;
T->GetEntry(iEvent); //Reads all the branches for entry(iEvent) and writes it to the respective variable(fHit or fRate)

//This loop goes over all the hits in the specific event
for (Int_t iHit = 0; iHit < fHit->size(); iHit++){
det = fHit->at(iHit).det;
energy = fHit->at(iHit).e;
pid = fHit->at(iHit).pid;
hitx = fHit->at(iHit).x;
hity = fHit->at(iHit).y;
hitr = fHit->at(iHit).r;
hitpz = fHit->at(iHit).pz;
rate = fRate/1.0e9/goodFileCount; // convert to GHz/uA

//Fill histograms with proper cuts
bool all_cuts = (det == 30 && //det number 30 is SM plane
hitr>1020 && hitr<1180)&&
(pid==11 || pid==-11 || pid==22 || pid==13 || pid==-13 || pid==211 || pid==-211 || pid==2112) && //particle selection
energy>2 && //energy cut
hitpz>0; //particle coming from upstream (pz>0)

if (all_cuts) {
std::pair<double, double> qsimxy = ConvertRemollToQsim(hitx, hity);
pe = get_PE_response(map_fit_data[pid], energy, qsimxy.first, qsimxy.second);

h_hitRate->Fill(hitr, rate);
h_hitRateSmPE->Fill(hitr,rate*pe);
}
}
}

// Print the rate in shower-max for script validation
double rateTotal = h_hitRate->Integral(); // in GHz
cout << "Accepted rate: " << rateTotal << " GHz" << endl;

double cathCurrent = h_hitRateSmPE->Integral()*1e9*1.6e-19*1e9; // in nA
cout << "Cathode current: " << cathCurrent << " nA" << endl;

// Plot the histograms
TCanvas* c1 = new TCanvas("c1", "c1", 1300, 500);
c1->Divide(2, 1);
c1->cd(1);
h_hitRate->Draw("hist E");
c1->cd(2);
h_hitRateSmPE->Draw("hist E");

}
69 changes: 69 additions & 0 deletions analysis/showermax/scripts/remollToQsim.hh
Original file line number Diff line number Diff line change
@@ -0,0 +1,69 @@
#ifndef REMOLL_TO_QSIM_HH
#define REMOLL_TO_QSIM_HH

#include <cmath>
#include <vector>
#include <TMath.h>

// Constants for the qsim system
const double RADIUS_CENTER = 1100.0; // mm, radial distance from remoll center to the center of qsim planes
const double QSIM_PLANE_WIDTH = 160.0; // mm, radial dimension of qsim planes
const double QSIM_PLANE_HEIGHT = 265.0; // mm, azimuthal straight dimension of qsim planes
const int N_PLANES = 28; // Number of qsim planes

// Function to calculate the azimuthal angle for each qsim plane
inline std::vector<double> CalculateQsimPlaneAngles() {
std::vector<double> planeAngles;
double deltaPhi = 2 * TMath::Pi() / N_PLANES;
for (int i = 0; i < N_PLANES; ++i) {
planeAngles.push_back(i * deltaPhi);
}
return planeAngles;
}

// Function to convert remoll (x, y) coordinates to qsim plane coordinates
inline std::pair<double, double> ConvertRemollToQsim(double x_remoll, double y_remoll) {
// Calculate the polar coordinates (r, phi) in the remoll system
double r_remoll = std::sqrt(x_remoll * x_remoll + y_remoll * y_remoll);
double phi_remoll = std::atan2(y_remoll, x_remoll); // Angle in radians

// Adjust

// Adjust phi_remoll to be in the range [0, 2pi)
if (phi_remoll < 0) { phi_remoll += 2 * TMath::Pi(); }

// Calculate qsim plane angles
std::vector<double> qsimAngles = CalculateQsimPlaneAngles();

// Determine the closest qsim plane
int closestPlane = 0;
double minDeltaPhi = TMath::Pi(); // Initialize to the largest possible value

for (size_t i = 0; i < qsimAngles.size(); ++i) {
double deltaPhi = std::fabs(phi_remoll - qsimAngles[i]);
deltaPhi = std::min(deltaPhi, 2 * TMath::Pi() - deltaPhi); // Account for circular nature

if (deltaPhi < minDeltaPhi) {
minDeltaPhi = deltaPhi;
closestPlane = i;
}
}

// Calculate the local coordinates on the closest qsim plane
double planeAngle = qsimAngles[closestPlane];
double x_plane_center = RADIUS_CENTER * std::cos(planeAngle);
double y_plane_center = RADIUS_CENTER * std::sin(planeAngle);

// Transform remoll coordinates to the local coordinate system of the qsim plane
double x_local = (x_remoll - x_plane_center) * std::cos(planeAngle) + (y_remoll - y_plane_center) * std::sin(planeAngle);
double y_local = - (x_remoll - x_plane_center) * std::sin(planeAngle) + (y_remoll - y_plane_center) * std::cos(planeAngle);

// Output the results
// std::cout << "Remoll coordinates: (" << x_remoll << ", " << y_remoll << ")\n";
// std::cout << "Closest qsim plane: " << closestPlane + 1 << "\n";
// std::cout << "Qsim local coordinates on plane " << closestPlane + 1 << ": (" << x_local << ", " << y_local << ")\n";

return {x_local, y_local};
}

#endif // REMOLL_TO_QSIM_HH
Loading

0 comments on commit 3fab176

Please sign in to comment.