Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
read-LEDs.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file read-LEDs.C
1 #include <calobase/TowerInfoDefs.h>
2 
3 // c++ includes --
4 #include <string>
5 #include <vector>
6 #include <iostream>
7 
8 // -- root includes --
9 #include <TH2F.h>
10 #include <TFile.h>
11 #include <TTree.h>
12 
13 using std::cout;
14 using std::endl;
15 using std::string;
16 using std::vector;
17 using std::min;
18 using std::max;
19 using std::to_string;
20 
21 R__LOAD_LIBRARY(libcalo_io.so)
22 
23 namespace myAnalysis {
24 
25  void init(const string& inputFile = "data/LEDTowerBuilder.root");
26  void analyze(UInt_t nevents=0);
27  void finalize(const string &outputFile="output/test.root");
28 
29  TFile* input;
30  TTree* led_tree;
31 
32  // QA
33  TH1F* hTime;
34  TH1F* hADC;
35  TH1F* hadc;
36  TH1F* hPed;
37  TH1F* hChannels;
38 
39  // 2D correlations
40  TH2F* h2TimeVsChannel;
41  TH2F* h2ADCVsChannel;
42  TH2F* h2PedVsChannel;
43  TH2F* h2ADCVsTime;
44  TH2F* h2PedVsTime;
45  TH2F* h2ADCVsPed;
46 
47  TH2F* h2TimeVsChannel_scat;
48  TH2F* h2ADCVsChannel_scat;
49  TH2F* h2PedVsChannel_scat;
50  TH2F* h2ADCVsTime_scat;
51  TH2F* h2PedVsTime_scat;
52  TH2F* h2ADCVsPed_scat;
53 
54  // one graph per channel
55  vector<TH2F*> h2adcVsTime;
56 
57  UInt_t time_bins = 32;
58  Float_t time_min = 0-0.5;
59  Float_t time_max = 32-0.5;
60 
61  UInt_t ADC_bins = 1800;
62  Float_t ADC_min = 0;
63  Float_t ADC_max = 18000;
64 
65  UInt_t adc_bins = 1800;
66  Float_t adc_min = 0;
67  Float_t adc_max = 18000;
68 
69  UInt_t ped_bins = 1650;
70  Float_t ped_min = 0;
71  Float_t ped_max = 16500;
72 
73  UInt_t channels_bins = 24576;
74  UInt_t channels_min = 0;
75  UInt_t channels_max = 24576;
76 }
77 
78 void myAnalysis::init(const string& inputFile) {
79  input = TFile::Open(inputFile.c_str());
80  led_tree = (TTree*)input->Get("W");
81 
82  // QA
83 
84  hTime = new TH1F("hTime", "Peak Location; Time Sample; Counts", time_bins, time_min, time_max);
85  hADC = new TH1F("hADC", "ADC; ADC; Counts", ADC_bins, ADC_min, ADC_max);
86  hadc = new TH1F("hadc", "adc; adc; Counts", adc_bins, adc_min, adc_max);
87  hPed = new TH1F("hPed", "Ped; Ped; Counts", ped_bins, ped_min, ped_max);
88  hChannels = new TH1F("hChannels", "Channels; Channel; Counts", channels_bins, channels_min, channels_max);
89 
90  // 2D correlations
91 
92  h2TimeVsChannel = new TH2F("h2TimeVsChannel", "Peak Location vs Channel; Channel; Time Sample", channels_bins, channels_min, channels_max, time_bins, time_min, time_max);
93  h2ADCVsChannel = new TH2F("h2ADCVsChannel", "ADC vs Channel; Channel; ADC", channels_bins, channels_min, channels_max, ADC_bins, ADC_min, ADC_max);
94  h2PedVsChannel = new TH2F("h2PedVsChannel", "Ped vs Channel; Channel; Ped", channels_bins, channels_min, channels_max, ped_bins, ped_min, ped_max);
95 
96  h2ADCVsTime = new TH2F("h2ADCVsTime", "ADC vs Peak Location; Time Sample; ADC", time_bins, time_min, time_max, ADC_bins, ADC_min, ADC_max);
97  h2PedVsTime = new TH2F("h2PedVsTime", "Ped vs Peak Location; Time Sample; Ped", time_bins, time_min, time_max, ped_bins, ped_min, ped_max);
98  h2ADCVsPed = new TH2F("h2ADCVsPed", "ADC vs Ped; Ped; ADC", ped_bins, ped_min, ped_max, ADC_bins, ADC_min, ADC_max);
99 
100  // scatter
101  h2TimeVsChannel_scat = new TH2F("h2TimeVsChannel_scat", "Peak Location vs Channel; Channel; Time Sample", channels_bins, channels_min, channels_max, time_bins, time_min, time_max);
102  h2ADCVsChannel_scat = new TH2F("h2ADCVsChannel_scat", "ADC vs Channel; Channel; ADC", channels_bins, channels_min, channels_max, ADC_bins, ADC_min, ADC_max);
103  h2PedVsChannel_scat = new TH2F("h2PedVsChannel_scat", "Ped vs Channel; Channel; Ped", channels_bins, channels_min, channels_max, ped_bins, ped_min, ped_max);
104 
105  h2ADCVsTime_scat = new TH2F("h2ADCVsTime_scat", "ADC vs Peak Location; Time Sample; ADC", time_bins, time_min, time_max, ADC_bins, ADC_min, ADC_max);
106  h2PedVsTime_scat = new TH2F("h2PedVsTime_scat", "Ped vs Peak Location; Time Sample; Ped", time_bins, time_min, time_max, ped_bins, ped_min, ped_max);
107  h2ADCVsPed_scat = new TH2F("h2ADCVsPed_scat", "ADC vs Ped; Ped; ADC", ped_bins, ped_min, ped_max, ADC_bins, ADC_min, ADC_max);
108 
109  for (UInt_t i = 0; i < channels_bins; ++i) {
110  h2adcVsTime.push_back(new TH2F(("h2adcVsTime_" + to_string(i)).c_str(),
111  "adc vs Time Sample; Time Sample; adc",
112  time_bins, time_min, time_max,
113  50, adc_min, adc_max));
114  }
115 }
116 
117 void myAnalysis::analyze(UInt_t nevents) {
118  vector<Float_t>* time = 0;
119  vector<Float_t>* ADC = 0;
120  vector<Float_t>* ped = 0;
121  vector<Int_t>* chan = 0;
122  vector<vector<Float_t>>* waveforms = 0; // 2D: channel x time sample
123 
124  led_tree->SetBranchAddress("time",&time);
125  led_tree->SetBranchAddress("adc",&ADC);
126  led_tree->SetBranchAddress("ped",&ped);
127  led_tree->SetBranchAddress("chan",&chan);
128  led_tree->SetBranchAddress("waveforms",&waveforms);
129 
130  // if nevents is 0 then use all events otherwise use the set number of events
131  nevents = (nevents) ? nevents : led_tree->GetEntries();
132 
133  Int_t counter_event[24576] = {0};
134  // maximum waveforms to overlap per channel
135  Int_t events_max = 100;
136  for(UInt_t i = 0; i < nevents; ++i) {
137  if(i%100 == 0) cout << "Progress: " << i*100./nevents << " %" << endl;
138 
139  led_tree->GetEntry(i);
140 
141  UInt_t nchannels = time->size();
142  channels_max = max(channels_max,nchannels);
143 
144  for(UInt_t j = 0; j < nchannels; ++j) {
145  Float_t time_val = time->at(j);
146  Float_t ADC_val = ADC->at(j);
147  Float_t ped_val = ped->at(j);
148  Int_t channel = chan->at(j);
149  if(channel >= 24576) {
150  cout << "invalid channel number: " << channel << ", event: " << i << endl;
151  continue;
152  }
153 
154  UInt_t key = TowerInfoDefs::encode_emcal(channel);
155  UInt_t etabin = TowerInfoDefs::getCaloTowerEtaBin(key);
156  UInt_t phibin = TowerInfoDefs::getCaloTowerPhiBin(key);
157 
158  Int_t time_bin = hTime->FindBin(time_val);
159  Int_t ADC_bin = hADC->FindBin(ADC_val);
160  Int_t ped_bin = hPed->FindBin(ped_val);
161  Int_t channel_bin = hChannels->FindBin(channel);
162 
163  for (UInt_t sample = 0; sample < time_bins; ++sample) {
164  Float_t adc_val = waveforms->at(j).at(sample);
165  // Int_t adc_bin = h2adcVsTime[channel]->GetYaxis()->FindBin(adc_val);
166 
167  hadc->Fill(adc_val);
168  // overlay at most #n events of waveforms for each channel
169  if(counter_event[channel] < events_max) {
170  // h2adcVsTime[channel]->SetBinContent(sample+1, adc_bin, 1);
171  h2adcVsTime[channel]->Fill(sample, adc_val);
172  }
173 
174  adc_min = min(adc_min, adc_val);
175  adc_max = max(adc_max, adc_val);
176  }
177  ++counter_event[channel];
178 
179  hChannels->Fill(channel);
180 
181  hTime->Fill(time_val);
182  hADC->Fill(ADC_val);
183  hPed->Fill(ped_val);
184 
185  h2TimeVsChannel->Fill(channel, time_val);
186  h2ADCVsChannel->Fill(channel, ADC_val);
187  h2PedVsChannel->Fill(channel, ped_val);
188 
189  h2TimeVsChannel_scat->SetBinContent(channel_bin, time_bin, 1);
190  h2ADCVsChannel_scat->SetBinContent(channel_bin, ADC_bin, 1);
191  h2PedVsChannel_scat->SetBinContent(channel_bin, ped_bin, 1);
192 
193  h2ADCVsTime->Fill(time_val, ADC_val);
194  h2PedVsTime->Fill(time_val, ped_val);
195  h2ADCVsPed->Fill(ped_val, ADC_val);
196 
197  h2ADCVsTime_scat->SetBinContent(time_bin, ADC_bin, 1);
198  h2PedVsTime_scat->SetBinContent(time_bin, ped_bin, 1);
199  h2ADCVsPed_scat->SetBinContent(ped_bin, ADC_bin, 1);
200 
201  time_min = min(time_min, time_val);
202  time_max = max(time_max, time_val);
203 
204  ADC_min = min(ADC_min, ADC_val);
205  ADC_max = max(ADC_max, ADC_val);
206 
207  ped_min = min(ped_min, ped_val);
208  ped_max = max(ped_max, ped_val);
209  }
210  }
211 
212  cout << "events processed: " << nevents << endl;
213  cout << "max channels per event: " << channels_max << endl;
214  cout << "time_min: " << time_min << " time_max: " << time_max << endl;
215  cout << "ADC_min: " << ADC_min << " ADC_max: " << ADC_max << endl;
216  cout << "adc_min: " << adc_min << " adc_max: " << adc_max << endl;
217  cout << "ped_min: " << ped_min << " ped_max: " << ped_max << endl;
218 }
219 
220 void myAnalysis::finalize(const string& outputFile) {
221  TFile output(outputFile.c_str(),"recreate");
222  output.cd();
223 
224  hChannels->Write();
225  hTime->Write();
226  hADC->Write();
227  hadc->Write();
228  hPed->Write();
229 
230  h2TimeVsChannel->Write();
231  h2ADCVsChannel->Write();
232  h2PedVsChannel->Write();
233 
234  h2ADCVsTime->Write();
235  h2PedVsTime->Write();
236  h2ADCVsPed->Write();
237 
238  output.mkdir("scat");
239  output.cd("scat");
240 
241  h2TimeVsChannel_scat->Write();
242  h2ADCVsChannel_scat->Write();
243  h2PedVsChannel_scat->Write();
244 
245  h2ADCVsTime_scat->Write();
246  h2PedVsTime_scat->Write();
247  h2ADCVsPed_scat->Write();
248 
249  output.cd();
250 
251  output.mkdir("adcVsTime");
252  output.cd("adcVsTime");
253  for(auto h2 : h2adcVsTime) {
254  if(h2->GetEntries()) h2->Write();
255  }
256 
257  // Close root file
258  input->Close();
259  output.Close();
260 }
261 
262 void convert_channel_to_tower_bin() {
263  for (UInt_t i = 0; i < 1536; ++i) {
264  UInt_t key = TowerInfoDefs::encode_emcal(i);
265  UInt_t etabin = TowerInfoDefs::getCaloTowerEtaBin(key);
266  UInt_t phibin = TowerInfoDefs::getCaloTowerPhiBin(key);
267  // cout << "channel: " << i << ", key: " << key << ", etabin: " << etabin << ", phibin: " << phibin << endl;
268  cout << "channel: " << i << ", etabin: " << etabin << ", phibin: " << phibin << endl;
269  }
270 }
271 
272 void read_LEDs(const string &inputFile="data/LEDTowerBuilder.root",
273  const string &outputFile="output/test.root",
274  UInt_t nevents=0) {
275 
276  myAnalysis::init(inputFile);
277  myAnalysis::analyze(nevents);
278  myAnalysis::finalize(outputFile);
279 
280  // convert_channel_to_tower_bin();
281 }
282 
283 # ifndef __CINT__
284 int main(int argc, char* argv[]) {
285  if(argc < 1 || argc > 4){
286  cout << "usage: ./bin/read-LEDs inputFile outputFile events" << endl;
287  cout << "inputFile: Location of LEDTowerBuilder.root. Default = data/LEDTowerBuilder.root." << endl;
288  cout << "outputFile: output root file. Default = output/test.root." << endl;
289  cout << "events: Number of events to analyze. Default = 0 (meaning all events)." << endl;
290  return 1;
291  }
292 
293  string input = "data/LEDTowerBuilder.root";
294  string output = "output/test.root";
295  UInt_t events = 0;
296 
297  if(argc >= 2) {
298  input = argv[1];
299  }
300  if(argc >= 3) {
301  output = argv[2];
302  }
303  if(argc >= 4) {
304  events = atoi(argv[3]);
305  }
306 
307  read_LEDs(input, output, events);
308 
309  cout << "done" << endl;
310  return 0;
311 }
312 # endif