EMCal-LEDs/macro/read-LEDs.C

0001 #include <calobase/TowerInfoDefs.h>
0002
0003 // c++ includes --
0004 #include <string>
0005 #include <vector>
0006 #include <iostream>
0007
0008 // -- root includes --
0009 #include <TH2F.h>
0010 #include <TFile.h>
0011 #include <TTree.h>
0012
0013 using std::cout;
0014 using std::endl;
0015 using std::string;
0016 using std::vector;
0017 using std::min;
0018 using std::max;
0019 using std::to_string;
0020
0021 R__LOAD_LIBRARY(libcalo_io.so)
0022
0023 namespace myAnalysis {
0024
0025     void init(const string& inputFile = "data/LEDTowerBuilder.root");
0026     void analyze(UInt_t nevents=0);
0027     void finalize(const string &outputFile="output/test.root");
0028
0029     TFile* input;
0030     TTree* led_tree;
0031
0032     // QA
0033     TH1F* hTime;
0034     TH1F* hADC;
0035     TH1F* hadc;
0036     TH1F* hPed;
0037     TH1F* hChannels;
0038
0039     // 2D correlations
0040     TH2F* h2TimeVsChannel;
0041     TH2F* h2ADCVsChannel;
0042     TH2F* h2PedVsChannel;
0043     TH2F* h2ADCVsTime;
0044     TH2F* h2PedVsTime;
0045     TH2F* h2ADCVsPed;
0046
0047     TH2F* h2TimeVsChannel_scat;
0048     TH2F* h2ADCVsChannel_scat;
0049     TH2F* h2PedVsChannel_scat;
0050     TH2F* h2ADCVsTime_scat;
0051     TH2F* h2PedVsTime_scat;
0052     TH2F* h2ADCVsPed_scat;
0053
0054     // one graph per channel
0055     vector<TH2F*> h2adcVsTime;
0056
0057     UInt_t  time_bins    = 32;
0058     Float_t time_min     = 0-0.5;
0059     Float_t time_max     = 32-0.5;
0060
0061     UInt_t  ADC_bins     = 1800;
0062     Float_t ADC_min      = 0;
0063     Float_t ADC_max      = 18000;
0064
0065     UInt_t  adc_bins     = 1800;
0066     Float_t adc_min      = 0;
0067     Float_t adc_max      = 18000;
0068
0069     UInt_t  ped_bins     = 1650;
0070     Float_t ped_min      = 0;
0071     Float_t ped_max      = 16500;
0072
0073     UInt_t channels_bins = 24576;
0074     UInt_t channels_min  = 0;
0075     UInt_t channels_max  = 24576;
0076 }
0077
0078 void myAnalysis::init(const string& inputFile) {
0079     input = TFile::Open(inputFile.c_str());
0080     led_tree = (TTree*)input->Get("W");
0081
0082     // QA
0083
0084     hTime = new TH1F("hTime", "Peak Location; Time Sample; Counts", time_bins, time_min, time_max);
0085     hADC = new TH1F("hADC", "ADC; ADC; Counts", ADC_bins, ADC_min, ADC_max);
0086     hadc = new TH1F("hadc", "adc; adc; Counts", adc_bins, adc_min, adc_max);
0087     hPed = new TH1F("hPed", "Ped; Ped; Counts", ped_bins, ped_min, ped_max);
0088     hChannels = new TH1F("hChannels", "Channels; Channel; Counts", channels_bins, channels_min, channels_max);
0089
0090     // 2D correlations
0091
0092     h2TimeVsChannel = new TH2F("h2TimeVsChannel", "Peak Location vs Channel; Channel; Time Sample", channels_bins, channels_min, channels_max, time_bins, time_min, time_max);
0093     h2ADCVsChannel = new TH2F("h2ADCVsChannel", "ADC vs Channel; Channel; ADC", channels_bins, channels_min, channels_max, ADC_bins, ADC_min, ADC_max);
0094     h2PedVsChannel = new TH2F("h2PedVsChannel", "Ped vs Channel; Channel; Ped", channels_bins, channels_min, channels_max, ped_bins, ped_min, ped_max);
0095
0096     h2ADCVsTime = new TH2F("h2ADCVsTime", "ADC vs Peak Location; Time Sample; ADC", time_bins, time_min, time_max, ADC_bins, ADC_min, ADC_max);
0097     h2PedVsTime = new TH2F("h2PedVsTime", "Ped vs Peak Location; Time Sample; Ped", time_bins, time_min, time_max, ped_bins, ped_min, ped_max);
0098     h2ADCVsPed = new TH2F("h2ADCVsPed", "ADC vs Ped; Ped; ADC", ped_bins, ped_min, ped_max, ADC_bins, ADC_min, ADC_max);
0099
0100     // scatter
0101     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);
0102     h2ADCVsChannel_scat = new TH2F("h2ADCVsChannel_scat", "ADC vs Channel; Channel; ADC", channels_bins, channels_min, channels_max, ADC_bins, ADC_min, ADC_max);
0103     h2PedVsChannel_scat = new TH2F("h2PedVsChannel_scat", "Ped vs Channel; Channel; Ped", channels_bins, channels_min, channels_max, ped_bins, ped_min, ped_max);
0104
0105     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);
0106     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);
0107     h2ADCVsPed_scat = new TH2F("h2ADCVsPed_scat", "ADC vs Ped; Ped; ADC", ped_bins, ped_min, ped_max, ADC_bins, ADC_min, ADC_max);
0108
0109     for (UInt_t i = 0; i < channels_bins; ++i) {
0110         h2adcVsTime.push_back(new TH2F(("h2adcVsTime_" + to_string(i)).c_str(),
0111                                          "adc vs Time Sample; Time Sample; adc",
0112                                          time_bins, time_min, time_max,
0113                                          50, adc_min, adc_max));
0114     }
0115 }
0116
0117 void myAnalysis::analyze(UInt_t nevents) {
0118     vector<Float_t>* time              = 0;
0119     vector<Float_t>* ADC               = 0;
0120     vector<Float_t>* ped               = 0;
0121     vector<Int_t>* chan                = 0;
0122     vector<vector<Float_t>>* waveforms = 0; // 2D: channel x time sample
0123
0124     led_tree->SetBranchAddress("time",&time);
0125     led_tree->SetBranchAddress("adc",&ADC);
0126     led_tree->SetBranchAddress("ped",&ped);
0127     led_tree->SetBranchAddress("chan",&chan);
0128     led_tree->SetBranchAddress("waveforms",&waveforms);
0129
0130     // if nevents is 0 then use all events otherwise use the set number of events
0131     nevents = (nevents) ? nevents : led_tree->GetEntries();
0132
0133     Int_t counter_event[24576] = {0};
0134     // maximum waveforms to overlap per channel
0135     Int_t events_max = 100;
0136     for(UInt_t i = 0; i < nevents; ++i) {
0137         if(i%100 == 0) cout << "Progress: " << i*100./nevents << " %" << endl;
0138
0139         led_tree->GetEntry(i);
0140
0141         UInt_t nchannels = time->size();
0142         channels_max = max(channels_max,nchannels);
0143
0144         for(UInt_t j = 0; j < nchannels; ++j) {
0145             Float_t time_val = time->at(j);
0146             Float_t ADC_val  = ADC->at(j);
0147             Float_t ped_val  = ped->at(j);
0148             Int_t channel    = chan->at(j);
0149             if(channel >= 24576) {
0150                 cout << "invalid channel number: " << channel << ", event: " << i << endl;
0151                 continue;
0152             }
0153
0154             UInt_t key    = TowerInfoDefs::encode_emcal(channel);
0155             UInt_t etabin = TowerInfoDefs::getCaloTowerEtaBin(key);
0156             UInt_t phibin = TowerInfoDefs::getCaloTowerPhiBin(key);
0157
0158             Int_t time_bin    = hTime->FindBin(time_val);
0159             Int_t ADC_bin     = hADC->FindBin(ADC_val);
0160             Int_t ped_bin     = hPed->FindBin(ped_val);
0161             Int_t channel_bin = hChannels->FindBin(channel);
0162
0163             for (UInt_t sample = 0; sample < time_bins; ++sample) {
0164                 Float_t adc_val = waveforms->at(j).at(sample);
0165                 // Int_t adc_bin   = h2adcVsTime[channel]->GetYaxis()->FindBin(adc_val);
0166
0167                 hadc->Fill(adc_val);
0168                 // overlay at most #n events of waveforms for each channel
0169                 if(counter_event[channel] < events_max) {
0170                     // h2adcVsTime[channel]->SetBinContent(sample+1, adc_bin, 1);
0171                     h2adcVsTime[channel]->Fill(sample, adc_val);
0172                 }
0173
0174                 adc_min = min(adc_min, adc_val);
0175                 adc_max = max(adc_max, adc_val);
0176             }
0177             ++counter_event[channel];
0178
0179             hChannels->Fill(channel);
0180
0181             hTime->Fill(time_val);
0182             hADC->Fill(ADC_val);
0183             hPed->Fill(ped_val);
0184
0185             h2TimeVsChannel->Fill(channel, time_val);
0186             h2ADCVsChannel->Fill(channel, ADC_val);
0187             h2PedVsChannel->Fill(channel, ped_val);
0188
0189             h2TimeVsChannel_scat->SetBinContent(channel_bin, time_bin, 1);
0190             h2ADCVsChannel_scat->SetBinContent(channel_bin, ADC_bin, 1);
0191             h2PedVsChannel_scat->SetBinContent(channel_bin, ped_bin, 1);
0192
0193             h2ADCVsTime->Fill(time_val, ADC_val);
0194             h2PedVsTime->Fill(time_val, ped_val);
0195             h2ADCVsPed->Fill(ped_val, ADC_val);
0196
0197             h2ADCVsTime_scat->SetBinContent(time_bin, ADC_bin, 1);
0198             h2PedVsTime_scat->SetBinContent(time_bin, ped_bin, 1);
0199             h2ADCVsPed_scat->SetBinContent(ped_bin, ADC_bin, 1);
0200
0201             time_min = min(time_min, time_val);
0202             time_max = max(time_max, time_val);
0203
0204             ADC_min = min(ADC_min, ADC_val);
0205             ADC_max = max(ADC_max, ADC_val);
0206
0207             ped_min = min(ped_min, ped_val);
0208             ped_max = max(ped_max, ped_val);
0209         }
0210     }
0211
0212     cout << "events processed: " << nevents << endl;
0213     cout << "max channels per event: " << channels_max << endl;
0214     cout << "time_min: " << time_min << " time_max: " << time_max << endl;
0215     cout << "ADC_min: " << ADC_min << " ADC_max: " << ADC_max << endl;
0216     cout << "adc_min: " << adc_min << " adc_max: " << adc_max << endl;
0217     cout << "ped_min: " << ped_min << " ped_max: " << ped_max << endl;
0218 }
0219
0220 void myAnalysis::finalize(const string& outputFile) {
0221     TFile output(outputFile.c_str(),"recreate");
0222     output.cd();
0223
0224     hChannels->Write();
0225     hTime->Write();
0226     hADC->Write();
0227     hadc->Write();
0228     hPed->Write();
0229
0230     h2TimeVsChannel->Write();
0231     h2ADCVsChannel->Write();
0232     h2PedVsChannel->Write();
0233
0234     h2ADCVsTime->Write();
0235     h2PedVsTime->Write();
0236     h2ADCVsPed->Write();
0237
0238     output.mkdir("scat");
0239     output.cd("scat");
0240
0241     h2TimeVsChannel_scat->Write();
0242     h2ADCVsChannel_scat->Write();
0243     h2PedVsChannel_scat->Write();
0244
0245     h2ADCVsTime_scat->Write();
0246     h2PedVsTime_scat->Write();
0247     h2ADCVsPed_scat->Write();
0248
0249     output.cd();
0250
0251     output.mkdir("adcVsTime");
0252     output.cd("adcVsTime");
0253     for(auto h2 : h2adcVsTime) {
0254         if(h2->GetEntries()) h2->Write();
0255     }
0256
0257     // Close root file
0258     input->Close();
0259     output.Close();
0260 }
0261
0262 void convert_channel_to_tower_bin() {
0263     for (UInt_t i = 0; i < 1536; ++i) {
0264         UInt_t key = TowerInfoDefs::encode_emcal(i);
0265         UInt_t etabin = TowerInfoDefs::getCaloTowerEtaBin(key);
0266         UInt_t phibin = TowerInfoDefs::getCaloTowerPhiBin(key);
0267         // cout << "channel: " << i << ", key: " << key << ", etabin: " << etabin << ", phibin: " << phibin << endl;
0268         cout << "channel: " << i << ", etabin: " << etabin << ", phibin: " << phibin << endl;
0269     }
0270 }
0271
0272 void read_LEDs(const string &inputFile="data/LEDTowerBuilder.root",
0273                const string &outputFile="output/test.root",
0274                UInt_t nevents=0) {
0275
0276     myAnalysis::init(inputFile);
0277     myAnalysis::analyze(nevents);
0278     myAnalysis::finalize(outputFile);
0279
0280     // convert_channel_to_tower_bin();
0281 }
0282
0283 # ifndef __CINT__
0284 int main(int argc, char* argv[]) {
0285     if(argc < 1 || argc > 4){
0286         cout << "usage: ./bin/read-LEDs inputFile outputFile events" << endl;
0287         cout << "inputFile: Location of LEDTowerBuilder.root. Default = data/LEDTowerBuilder.root." << endl;
0288         cout << "outputFile: output root file. Default = output/test.root." << endl;
0289         cout << "events: Number of events to analyze. Default = 0 (meaning all events)." << endl;
0290         return 1;
0291     }
0292
0293     string input  = "data/LEDTowerBuilder.root";
0294     string output = "output/test.root";
0295     UInt_t events = 0;
0296
0297     if(argc >= 2) {
0298         input = argv[1];
0299     }
0300     if(argc >= 3) {
0301         output = argv[2];
0302     }
0303     if(argc >= 4) {
0304         events = atoi(argv[3]);
0305     }
0306
0307     read_LEDs(input, output, events);
0308
0309     cout << "done" << endl;
0310     return 0;
0311 }
0312 # endif