Add Showermax stackplanes and update response lookup

analysis/showermax/response_lookup/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/response_lookup/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/response_lookup/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 → analysis/showermax/response_lookup/scripts/plot_showermax_response.cpp

@@ -5,6 +5,7 @@
 // Purpose: Read the remoll rootfile and use shower-max lookup table to plot the PE response.
 // -----------------------------------------------------------------------------------
 
+#include "RtypesCore.h"
 #include "TChain.h"
 #include "TH1.h"
 #include "TROOT.h"
@@ -17,44 +18,29 @@
 #include "./remollToQsim.hh"
 #include "./shower-max_resposne_lookup.hh"
 
+#ifdef __REMOLLTYPES__ // for clangd LSP support
+  #include "../include/remolltypes.hh"
+#endif // __REMOLLTYPES__
+
+
 // 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;
+void plot_showermax_response(TString rootFile = "~/moller/softwares/remoll-zone/rootfiles/sm_stack/stack_v27/sm_tqStack_3sector_moller_50k_1001.root") {
+
+    // Testing shower-max response 
+    ShowermaxLookup::getPeResponse(73015, 11, 1000, 0.0, 1080.0);
 
     // 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); 
+    TH1D* h_hitRateSmPE = new TH1D("hitRateSmPE", "PE*Rate weighted hit; hit.r [mm]; rate [PE*GHz/65uA]", 100, 1000, 1200); 
 
     // Declare TChain
     TChain* T = new TChain("T");
@@ -74,7 +60,7 @@ void plot_showermax_response(){
 	T->SetBranchAddress("rate", &fRate);
 	// T->SetBranchAddress("ev", &fEv);
 
-	//This loop goes over all the events in the root files
+	// 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;
@@ -92,15 +78,16 @@ void plot_showermax_response(){
 			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
+			bool all_cuts = (det >= 73001 && det <= 73028 &&											    		//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
+							energy>10 &&  										//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);
+                // Double_t pe = ShowermaxLookup::getPeResponseWithoutLP(pid, energy, hitx, hity);
+                Double_t pe = ShowermaxLookup::getPeResponse(det, pid, energy, hitx, hity);
+                std::cout << "pid: " << pid << ", energy: " << energy << ", pe: " << pe << std::endl;
 
                 h_hitRate->Fill(hitr, rate);
                 h_hitRateSmPE->Fill(hitr,rate*pe);
@@ -109,18 +96,19 @@ void plot_showermax_response(){
     }
 
     // Print the rate in shower-max for script validation
-    double rateTotal = h_hitRate->Integral(); // in GHz
+    double rateTotal = h_hitRate->Integral(); // in GH
     cout << "Accepted rate: " << rateTotal << " GHz" << endl;
 
     double cathCurrent = h_hitRateSmPE->Integral()*1e9*1.6e-19*1e9; // in nA
-    cout << "Cathode current: " << cathCurrent << " nA" << endl;
+    cout << "Total 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");
+    // 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 → analysis/showermax/response_lookup/scripts/remollToQsim.hh

@@ -30,7 +30,7 @@ inline std::pair<double, double> ConvertRemollToQsim(double x_remoll, double y_r
     // Adjust 
 
     // Adjust phi_remoll to be in the range [0, 2pi)
-    if (phi_remoll < 0) { phi_remoll += 2 * TMath::Pi(); }
+    if (phi_remoll < 0) phi_remoll += 2 * TMath::Pi();
 
     // Calculate qsim plane angles
     std::vector<double> qsimAngles = CalculateQsimPlaneAngles();
@@ -67,3 +67,4 @@ inline std::pair<double, double> ConvertRemollToQsim(double x_remoll, double y_r
 }
 
 #endif // REMOLL_TO_QSIM_HH
+