|
| 1 | +// ROOTscript: plot_shower-max_response.C |
| 2 | +// ----------------------------------------------------------------------------------- |
| 3 | +// Author: Sudip Bhattarai |
| 4 | +// Date: 01/21/2025 |
| 5 | +// Purpose: Read the remoll rootfile and use shower-max lookup table to plot the PE response. |
| 6 | +// ----------------------------------------------------------------------------------- |
| 7 | + |
| 8 | +#include "TChain.h" |
| 9 | +#include "TH1.h" |
| 10 | +#include "TROOT.h" |
| 11 | +#include "TCanvas.h" |
| 12 | +#include <iostream> |
| 13 | +#include <fstream> |
| 14 | +#include <utility> |
| 15 | +#include <vector> |
| 16 | +#include <string> |
| 17 | +#include "./remollToQsim.hh" |
| 18 | +#include "./shower-max_resposne_lookup.hh" |
| 19 | + |
| 20 | +// Namespace list to use |
| 21 | +using std::cout, std::cerr, std::endl; |
| 22 | +using std::vector, std::string; |
| 23 | +using std::ifstream, std::ofstream; |
| 24 | +using std::pair, std::map; |
| 25 | + |
| 26 | +// Main function |
| 27 | +void plot_showermax_response(){ |
| 28 | + // Genetate random values for x, y, theta, phi |
| 29 | + |
| 30 | + // Get the fit parameters |
| 31 | + vector<vector<Double_t>> fit_data_electron = retrieve_fit_data("e-"); |
| 32 | + vector<vector<Double_t>> fit_data_gamma = retrieve_fit_data("gamma"); |
| 33 | + vector<vector<Double_t>> fit_data_mu = retrieve_fit_data("mu-"); |
| 34 | + vector<vector<Double_t>> fit_data_pi = retrieve_fit_data("pi-"); |
| 35 | + vector<vector<Double_t>> fit_data_neutron = retrieve_fit_data("neutron"); |
| 36 | + |
| 37 | + // Make a map of pid and fit_data |
| 38 | + map<int, vector<vector<Double_t>>> map_fit_data; |
| 39 | + map_fit_data[11] = fit_data_electron; |
| 40 | + map_fit_data[-11] = fit_data_electron; |
| 41 | + map_fit_data[22] = fit_data_gamma; |
| 42 | + map_fit_data[13] = fit_data_mu; |
| 43 | + map_fit_data[-13] = fit_data_mu; |
| 44 | + map_fit_data[211] = fit_data_pi; |
| 45 | + map_fit_data[-211] = fit_data_pi; |
| 46 | + map_fit_data[2112] = fit_data_neutron; |
| 47 | + |
| 48 | + // Test pe response |
| 49 | + // cout << "PE: " << get_PE_response(map_fit_data[11], 1234, 0, 0) << endl; |
| 50 | + |
| 51 | + // Remoll file path goes here |
| 52 | + int goodFileCount = 1; // Number of good (non-corrupted) files |
| 53 | + TString rootFile= "~/moller/softwares/remoll-zone/rootfiles/smStack_v23/sm_tqStack_3sector_moller_50k_1001.root"; |
| 54 | + |
| 55 | + // Define histograms |
| 56 | + TH1D* h_hitRate = new TH1D("hitRate", "Rate weighted hit; hit.r [mm]; rate[GHz/65uA]", 100, 1000, 1200); |
| 57 | + TH1D* h_hitRateSmPE = new TH1D("hitRateSmPE", "Rate weighted hit; hit.r [mm]; rate [GHz/65uA/PE]", 100, 1000, 1200); |
| 58 | + |
| 59 | + // Declare TChain |
| 60 | + TChain* T = new TChain("T"); |
| 61 | + T->Add(rootFile); |
| 62 | + |
| 63 | + Long64_t nEntries = T->GetEntries(); |
| 64 | + cout << "Total number of entries in the chain: " << nEntries << endl; |
| 65 | + std::vector<remollGenericDetectorHit_t> *fHit = 0; |
| 66 | + // remollEvent_t *fEv = 0; |
| 67 | + Double_t fRate = 0; |
| 68 | + Float_t energy(-1.0e-12), rate(-1.0e-12), |
| 69 | + hitx(-1.0e-12), hity(-1.0e-12), hitr(-1.0e-12), hitpz(-1.0e-12), |
| 70 | + pe(0); |
| 71 | + Int_t pid(0), det(0); |
| 72 | + |
| 73 | + T->SetBranchAddress("hit", &fHit); |
| 74 | + T->SetBranchAddress("rate", &fRate); |
| 75 | + // T->SetBranchAddress("ev", &fEv); |
| 76 | + |
| 77 | + //This loop goes over all the events in the root files |
| 78 | + for (Long64_t iEvent = 0; iEvent < nEntries; iEvent++){ |
| 79 | + if (iEvent % (nEntries/10) == 0) |
| 80 | + cout << "Analyzed " << iEvent << " events." << endl; |
| 81 | + T->GetEntry(iEvent); //Reads all the branches for entry(iEvent) and writes it to the respective variable(fHit or fRate) |
| 82 | + |
| 83 | + //This loop goes over all the hits in the specific event |
| 84 | + for (Int_t iHit = 0; iHit < fHit->size(); iHit++){ |
| 85 | + det = fHit->at(iHit).det; |
| 86 | + energy = fHit->at(iHit).e; |
| 87 | + pid = fHit->at(iHit).pid; |
| 88 | + hitx = fHit->at(iHit).x; |
| 89 | + hity = fHit->at(iHit).y; |
| 90 | + hitr = fHit->at(iHit).r; |
| 91 | + hitpz = fHit->at(iHit).pz; |
| 92 | + rate = fRate/1.0e9/goodFileCount; // convert to GHz/uA |
| 93 | + |
| 94 | + //Fill histograms with proper cuts |
| 95 | + bool all_cuts = (det == 30 && //det number 30 is SM plane |
| 96 | + hitr>1020 && hitr<1180)&& |
| 97 | + (pid==11 || pid==-11 || pid==22 || pid==13 || pid==-13 || pid==211 || pid==-211 || pid==2112) && //particle selection |
| 98 | + energy>2 && //energy cut |
| 99 | + hitpz>0; //particle coming from upstream (pz>0) |
| 100 | + |
| 101 | + if (all_cuts) { |
| 102 | + std::pair<double, double> qsimxy = ConvertRemollToQsim(hitx, hity); |
| 103 | + pe = get_PE_response(map_fit_data[pid], energy, qsimxy.first, qsimxy.second); |
| 104 | + |
| 105 | + h_hitRate->Fill(hitr, rate); |
| 106 | + h_hitRateSmPE->Fill(hitr,rate*pe); |
| 107 | + } |
| 108 | + } |
| 109 | + } |
| 110 | + |
| 111 | + // Print the rate in shower-max for script validation |
| 112 | + double rateTotal = h_hitRate->Integral(); // in GHz |
| 113 | + cout << "Accepted rate: " << rateTotal << " GHz" << endl; |
| 114 | + |
| 115 | + double cathCurrent = h_hitRateSmPE->Integral()*1e9*1.6e-19*1e9; // in nA |
| 116 | + cout << "Cathode current: " << cathCurrent << " nA" << endl; |
| 117 | + |
| 118 | + // Plot the histograms |
| 119 | + TCanvas* c1 = new TCanvas("c1", "c1", 1300, 500); |
| 120 | + c1->Divide(2, 1); |
| 121 | + c1->cd(1); |
| 122 | + h_hitRate->Draw("hist E"); |
| 123 | + c1->cd(2); |
| 124 | + h_hitRateSmPE->Draw("hist E"); |
| 125 | + |
| 126 | +} |
0 commit comments