Shower-max PE response lookup table analysis script

analysis/showermax/data/fit_param_xy_e-_ifarm.csv

@@ -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

analysis/showermax/data/fit_param_xy_gamma_ifarm.csv

@@ -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

analysis/showermax/data/fit_param_xy_mu-_ifarm.csv

@@ -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

analysis/showermax/data/fit_param_xy_neutron_ifarm.csv

@@ -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

analysis/showermax/data/fit_param_xy_pi-_ifarm.csv

@@ -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

analysis/showermax/scripts/plot_showermax_response.cpp

@@ -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");
+
+}

analysis/showermax/scripts/remollToQsim.hh

@@ -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