Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:12:18 

0001 #include <calobase/TowerInfoDefs.h>
0002 
0003 // c++ includes --
0004 #include <string>
0005 #include <vector>
0006 #include <iostream>
0007 #include <fstream>
0008 
0009 // -- root includes --
0010 #include <TH2F.h>
0011 #include <TFile.h>
0012 #include <TTree.h>
0013 #include <TLine.h>
0014 #include <TCanvas.h>
0015 
0016 using std::cout;
0017 using std::endl;
0018 using std::string;
0019 using std::vector;
0020 using std::min;
0021 using std::max;
0022 using std::to_string;
0023 using std::ofstream;
0024 
0025 R__LOAD_LIBRARY(libcalo_io.so)
0026 
0027 namespace myAnalysis {
0028 
0029     void init(const string& input_CEMC, const string& outputFile = "event-display/test.json", const string& input_HCALIN = "", const string& input_HCALOUT = "");
0030     void analyze();
0031     void finalize(const string& run, const string& event);
0032 
0033     ofstream output;
0034 
0035     TFile* finput_CEMC = 0;
0036     TFile* finput_HCALIN = 0;
0037     TFile* finput_HCALOUT = 0;
0038 
0039     TTree* tree_CEMC;
0040     TTree* tree_HCALIN;
0041     TTree* tree_HCALOUT;
0042 
0043     TH2F* h2ADC_CEMC;
0044     TH2F* h2ADC_HCALIN;
0045     TH2F* h2ADC_HCALOUT;
0046 
0047     TH1F* hADC_CEMC;
0048     TH1F* hADC_HCALIN;
0049     TH1F* hADC_HCALOUT;
0050 
0051     // CEMC
0052 
0053     UInt_t  eta_bins = 96;
0054     Float_t eta_min  = -0.5;
0055     Float_t eta_max  = 95.5;
0056 
0057     UInt_t  phi_bins = 256;
0058     Float_t phi_min  = -0.5;
0059     Float_t phi_max  = 255.5;
0060 
0061     UInt_t etabin_max = 95;
0062     UInt_t phibin_max = 255;
0063 
0064     UInt_t channels_bins = 24576;
0065     UInt_t channels_min  = 0;
0066     UInt_t channels_max  = 24576;
0067 
0068     // HCAL
0069 
0070     UInt_t  eta_hcal_bins = 24;
0071     Float_t eta_hcal_min  = -0.5;
0072     Float_t eta_hcal_max  = 23.5;
0073 
0074     UInt_t  phi_hcal_bins = 64;
0075     Float_t phi_hcal_min  = -0.5;
0076     Float_t phi_hcal_max  = 63.5;
0077 
0078     UInt_t channels_hcal_bins = 1536;
0079     UInt_t channels_hcal_min  = 0;
0080     UInt_t channels_hcal_max  = 1536;
0081 
0082     const float eta_map[] = {-1.05417, -0.9625,  -0.870833,-0.779167,-0.6875,  -0.595833,-0.504167,-0.4125,
0083                              -0.320833,-0.229167,-0.1375,  -0.0458333,0.0458333,0.1375,   0.229167, 0.320833,
0084                               0.4125,   0.504167, 0.595833, 0.6875,   0.779167, 0.870833, 0.9625,   1.05417};
0085 
0086     const float phi_map[] = {0.0490874, 0.147262, 0.245437, 0.343612, 0.441786, 0.539961, 0.638136, 0.736311,
0087                              0.834486,  0.93266,  1.03084,  1.12901,  1.22718,  1.32536,  1.42353,  1.52171,
0088                              1.61988,   1.71806,  1.81623,  1.91441,  2.01258,  2.11076,  2.20893,  2.30711,
0089                              2.40528,   2.50346,  2.60163,  2.69981,  2.79798,  2.89616,  2.99433,  3.09251,
0090                              3.19068,   3.28885,  3.38703,  3.4852,   3.58338,  3.68155,  3.77973,  3.8779,
0091                              3.97608,   4.07425,  4.17243,  4.2706,   4.36878,  4.46695,  4.56513,  4.6633,
0092                              4.76148,   4.85965,  4.95783,  5.056,    5.15418,  5.25235,  5.35052,  5.4487,
0093                              5.54687,   5.64505,  5.74322,  5.8414,   5.93957,  6.03775,  6.13592,  6.2341};
0094 
0095     UInt_t  ADC_bins = 65;
0096     Float_t ADC_min  = 0;
0097     Float_t ADC_max  = 1.3e3;
0098 }
0099 
0100 void myAnalysis::init(const string& input_CEMC, const string& outputFile, const string& input_HCALIN, const string& input_HCALOUT) {
0101 
0102     finput_CEMC    = TFile::Open(input_CEMC.c_str());
0103     if(!input_HCALIN.empty()) finput_HCALIN  = TFile::Open(input_HCALIN.c_str());
0104     if(!input_HCALOUT.empty()) finput_HCALOUT = TFile::Open(input_HCALOUT.c_str());
0105 
0106     tree_CEMC    = (TTree*)finput_CEMC->Get("W");
0107     tree_HCALIN  = (TTree*)finput_HCALIN->Get("W");
0108     tree_HCALOUT = (TTree*)finput_HCALOUT->Get("W");
0109 
0110     output.open(outputFile.c_str());
0111 
0112     hADC_CEMC  = new TH1F("hADC_CEMC", "ADC CEMC; ADC; Counts", ADC_bins, ADC_min, ADC_max);
0113     h2ADC_CEMC = new TH2F("h2ADC_CEMC", "ADC CEMC; towerid #eta; towerid #phi", eta_bins, eta_min, eta_max,
0114                                                                                 phi_bins, phi_min, phi_max);
0115 
0116     if(finput_HCALIN) {
0117         hADC_HCALIN  = new TH1F("hADC_HCALIN", "ADC HCALIN; ADC; Counts", ADC_bins, ADC_min, ADC_max);
0118         h2ADC_HCALIN = new TH2F("h2ADC_HCALIN", "ADC HCALIN; towerid #eta; towerid #phi", eta_hcal_bins, eta_hcal_min, eta_hcal_max,
0119                                                                                           phi_hcal_bins, phi_hcal_min, phi_hcal_max);
0120     }
0121 
0122     if(finput_HCALOUT) {
0123         hADC_HCALOUT  = new TH1F("hADC_HCALOUT", "ADC HCALOUT; ADC; Counts", ADC_bins, ADC_min, ADC_max);
0124         h2ADC_HCALOUT = new TH2F("h2ADC_HCALOUT", "ADC HCALOUT; towerid #eta; towerid #phi", eta_hcal_bins, eta_hcal_min, eta_hcal_max,
0125                                                                                              phi_hcal_bins, phi_hcal_min, phi_hcal_max);
0126     }
0127 
0128 output << "{\n \
0129     \"EVENT\": {\n \
0130       \"runid\": 1,\n \
0131       \"evtid\": 1,\n \
0132       \"time\": 0,\n \
0133       \"type\": \"Collision\",\n \
0134       \"s_nn\": 0,\n \
0135       \"B\": 3.0,\n \
0136       \"pv\": [0,0,0]\n \
0137     },\n \
0138     \"META\": {\n \
0139       \"HITS\": {\n \
0140         \"INNERTRACKER\": {\n \
0141           \"type\": \"3D\",\n \
0142           \"options\": {\n \
0143             \"size\": 5,\n \
0144             \"color\": 16777215\n \
0145           }\n \
0146         },\n \
0147         \"CEMC\": {\n \
0148           \"type\": \"PROJECTIVE\",\n \
0149             \"options\": {\n \
0150               \"rmin\": 90,\n \
0151               \"rmax\": 116.1,\n \
0152               \"deta\": 0.024,\n \
0153               \"dphi\": 0.024,\n \
0154               \"color\": 16766464,\n \
0155               \"transparent\": 0.6,\n \
0156               \"scaleminmax\": false\n \
0157             }\n \
0158         },\n \
0159         \"HCALIN\": {\n \
0160           \"type\": \"PROJECTIVE\",\n \
0161             \"options\": {\n \
0162               \"rmin\": 117.27,\n \
0163               \"rmax\": 134.42,\n \
0164               \"deta\": 0.025,\n \
0165               \"dphi\": 0.025,\n \
0166               \"color\": 4290445312,\n \
0167               \"transparent\": 0.6,\n \
0168               \"scaleminmax\": false\n \
0169             }\n \
0170         },\n \
0171         \"HCALOUT\": {\n \
0172           \"type\": \"PROJECTIVE\",\n \
0173             \"options\": {\n \
0174               \"rmin\": 183.3,\n \
0175               \"rmax\": 274.317,\n \
0176               \"deta\": 0.025,\n \
0177               \"dphi\": 0.025,\n \
0178               \"color\": 24773,\n \
0179               \"transparent\": 0.6,\n \
0180               \"scaleminmax\": true\n \
0181             }\n \
0182         }\n \
0183       }\n \
0184     },\n \
0185 \"HITS\": {\n";
0186 }
0187 
0188 void myAnalysis::analyze() {
0189 
0190     bool first_entry = true;
0191     vector<Float_t>* ADC = 0;
0192     vector<Int_t>* chan  = 0;
0193     // vector<Float_t>* time              = 0;
0194     // vector<Float_t>* ped               = 0;
0195     // vector<vector<Float_t>>* waveforms = 0; // 2D: channel x time sample
0196     Float_t max_ADC = ADC_max;
0197     if(finput_CEMC && tree_CEMC->GetEntries()) {
0198         tree_CEMC->SetBranchAddress("adc",&ADC);
0199         tree_CEMC->SetBranchAddress("chan",&chan);
0200         tree_CEMC->GetEntry(0);
0201         UInt_t nchannels = ADC->size();
0202 
0203         for(UInt_t j = 0; j < nchannels; ++j) {
0204             // convert from high gain to low gain (factor of 16)
0205             Float_t ADC_val  = ADC->at(j)/16;
0206             // Float_t e        = 0.002*ADC_val/16;
0207             Float_t e        = ADC_val;
0208             // if(ADC_val < 4e3) continue;
0209             Int_t channel    = chan->at(j);
0210             if(channel >= channels_max) {
0211                 cout << "CEMC invalid channel number: " << channel << endl;
0212                 continue;
0213             }
0214 
0215             UInt_t key    = TowerInfoDefs::encode_emcal(channel);
0216             UInt_t etabin = TowerInfoDefs::getCaloTowerEtaBin(key);
0217             UInt_t phibin = TowerInfoDefs::getCaloTowerPhiBin(key);
0218 
0219             h2ADC_CEMC->SetBinContent(etabin+1, phibin+1, ADC_val);
0220             hADC_CEMC->Fill(ADC_val);
0221 
0222             max_ADC = max(max_ADC, ADC_val);
0223 
0224             // mapping: [0,95] -> [-1.1, 1.1]
0225             Float_t eta = 2.2/etabin_max*etabin-1.1;
0226             // mapping: [0,127] -> [0, pi] and [128, 255] -> [-pi, 0]
0227             Float_t phi = (phibin <= phibin_max/2) ? M_PI/(phibin_max/2)*phibin :
0228                                                     M_PI/(phibin_max/2)*(phibin*1.-phibin_max);
0229 
0230             if(first_entry) {
0231                 output << "\"CEMC\": [{ \"eta\": " << eta << ", \"phi\": " << phi << " , \"e\": " << e << "}\n";
0232                 first_entry = false;
0233             }
0234             else output << ",{ \"eta\": " << eta << ", \"phi\": " << phi << " , \"e\": " << e << "}\n";
0235             // if(j == 10) break;
0236         }
0237         cout << "CEMC Max ADC: " << max_ADC << endl;
0238         output << "],\n";
0239     }
0240     else output << "\"CEMC\": [{}],\n";
0241 
0242     first_entry = true;
0243     max_ADC = ADC_max;
0244     if(finput_HCALIN && tree_HCALIN->GetEntries()) {
0245         ADC  = 0;
0246         chan = 0;
0247 
0248         tree_HCALIN->SetBranchAddress("adc",&ADC);
0249         tree_HCALIN->SetBranchAddress("chan",&chan);
0250         tree_HCALIN->GetEntry(0);
0251 
0252         UInt_t nchannels = ADC->size();
0253 
0254         for(UInt_t j = 0; j < nchannels; ++j) {
0255             Float_t ADC_val  = ADC->at(j);
0256             // Float_t e        = 0.002*ADC_val;
0257             Float_t e        = ADC_val;
0258             Int_t channel    = chan->at(j);
0259             if(channel >= channels_max) {
0260                 cout << "HCALIN invalid channel number: " << channel << endl;
0261                 continue;
0262             }
0263 
0264             UInt_t key    = TowerInfoDefs::encode_hcal(channel);
0265             UInt_t etabin = TowerInfoDefs::getCaloTowerEtaBin(key);
0266             UInt_t phibin = TowerInfoDefs::getCaloTowerPhiBin(key);
0267 
0268             h2ADC_HCALIN->SetBinContent(etabin+1, phibin+1, ADC_val);
0269             hADC_HCALIN->Fill(ADC_val);
0270 
0271             max_ADC = max(max_ADC, ADC_val);
0272 
0273             Float_t eta = eta_map[etabin];
0274             Float_t phi = phi_map[phibin];
0275 
0276             if(first_entry){
0277               output << "\"HCALIN\": [{ \"eta\": " << eta << ", \"phi\": " << phi << " , \"e\": " << e << "}\n";
0278               first_entry = false;
0279             }
0280             else output << ",{ \"eta\": " << eta << ", \"phi\": " << phi << " , \"e\": " << e << "}\n";
0281             // if(j == 10) break;
0282         }
0283         cout << "HCALIN Max ADC: " << max_ADC << endl;
0284         output << "],\n";
0285     }
0286     else output << "\"HCALIN\": [{}],\n";
0287 
0288     first_entry = true;
0289     max_ADC = ADC_max;
0290     if(finput_HCALOUT && tree_HCALOUT->GetEntries()) {
0291         ADC  = 0;
0292         chan = 0;
0293 
0294         tree_HCALOUT->SetBranchAddress("adc",&ADC);
0295         tree_HCALOUT->SetBranchAddress("chan",&chan);
0296         tree_HCALOUT->GetEntry(0);
0297 
0298         UInt_t nchannels = ADC->size();
0299 
0300         for(UInt_t j = 0; j < nchannels; ++j) {
0301             Float_t ADC_val  = ADC->at(j);
0302             // Float_t e        = 0.002*ADC_val;
0303             Float_t e        = ADC_val;
0304             // if(ADC_val < 100) continue;
0305             Int_t channel    = chan->at(j);
0306             if(channel >= channels_max) {
0307                 cout << "HCALOUT invalid channel number: " << channel << endl;
0308                 continue;
0309             }
0310 
0311             UInt_t key    = TowerInfoDefs::encode_hcal(channel);
0312             UInt_t etabin = TowerInfoDefs::getCaloTowerEtaBin(key);
0313             UInt_t phibin = TowerInfoDefs::getCaloTowerPhiBin(key);
0314 
0315             h2ADC_HCALOUT->SetBinContent(etabin+1, phibin+1, ADC_val);
0316             hADC_HCALOUT->Fill(ADC_val);
0317 
0318             max_ADC = max(max_ADC, ADC_val);
0319 
0320             Float_t eta = eta_map[etabin];
0321             Float_t phi = phi_map[phibin];
0322 
0323             if(first_entry) {
0324                 output << "\"HCALOUT\": [{ \"eta\": " << eta << ", \"phi\": " << phi << " , \"e\": " << e << "}\n";
0325                 first_entry = false;
0326             }
0327             else output << ",{ \"eta\": " << eta << ", \"phi\": " << phi << " , \"e\": " << e << "}\n";
0328             // if(j == 10) break;
0329             // cout << "etabin: " << etabin << ", phibin: " << phibin << ", ADC: " << ADC_val << endl;
0330             // cout << ",{ \"eta\": " << eta << ", \"phi\": " << phi << " , \"e\": " << ADC_val << "}" << endl;
0331         }
0332         cout << "HCALOUT Max ADC: " << max_ADC << endl;
0333         output << "]\n";
0334     }
0335     else output << "\"HCALOUT\": [{}]\n";
0336 }
0337 
0338 void myAnalysis::finalize(const string& run, const string& event) {
0339 
0340         output << "},\n \
0341     \"TRACKS\": {\n \
0342     \"INNERTRACKER\": []\n \
0343     }}";
0344 
0345 
0346     auto c1 = new TCanvas();
0347 
0348     c1->cd();
0349     // c1->SetTickx();
0350     // c1->SetTicky();
0351 
0352     c1->SetCanvasSize(1000, 1500);
0353 
0354     c1->SetLeftMargin(.12);
0355     c1->SetRightMargin(.15);
0356     // gPad->SetLogz();
0357 
0358     string outputFile = "event-display/plots-" + run + "-" + event + ".pdf";
0359     c1->Print((outputFile + "[").c_str(), "pdf portrait");
0360 
0361     gPad->SetLogz();
0362     h2ADC_CEMC->SetStats(0);
0363     h2ADC_CEMC->SetTitle(("ADC CEMC, Run: " + run + ", event: " + event).c_str());
0364     // h2ADC_CEMC->SetTickLength(0,"xy");
0365     // h2ADC_CEMC->SetTitleOffset(0.5, "z");
0366     // h2ADC_CEMC->SetNdivisions(12,"x");
0367     // h2ADC_CEMC->SetNdivisions(32,"y");
0368     h2ADC_CEMC->Draw("COLZ1");
0369 
0370     auto tline = new TLine();
0371 
0372     for (UInt_t i = 0; i < 32; ++i) {
0373         tline->DrawLine(eta_min, 8*i-0.5, eta_max, 8*i-0.5);
0374     }
0375     for (UInt_t i = 0; i < 12; ++i) {
0376         tline->DrawLine(8*i-0.5, phi_min, 8*i-0.5, phi_max);
0377     }
0378     c1->Print((outputFile).c_str(), "pdf portrait");
0379 
0380     // hADC_CEMC->SetStats(0);
0381     c1->SetCanvasSize(1500, 1000);
0382     gPad->SetLogy();
0383     hADC_CEMC->SetTitle(("ADC CEMC, Run: " + run + ", event: " + event).c_str());
0384     hADC_CEMC->Draw();
0385     c1->Print((outputFile).c_str(), "pdf portrait");
0386     gPad->SetLogy(0);
0387 
0388     // gPad->SetLogz(0);
0389     // h2ADC_CEMC->SetMinimum(4000);
0390     // c1->Print((outputFile).c_str(), "pdf portrait");
0391 
0392     if(finput_HCALIN) {
0393         c1->SetCanvasSize(1000, 1500);
0394         gPad->SetLogz();
0395         h2ADC_HCALIN->SetStats(0);
0396         h2ADC_HCALIN->SetTitle(("ADC HCALIN, Run: " + run + ", event: " + event).c_str());
0397         h2ADC_HCALIN->Draw("COLZ1");
0398         c1->Print((outputFile).c_str(), "pdf portrait");
0399 
0400         c1->SetCanvasSize(1500, 1000);
0401         gPad->SetLogy();
0402         hADC_HCALIN->SetTitle(("ADC HCALIN, Run: " + run + ", event: " + event).c_str());
0403         hADC_HCALIN->Draw();
0404         c1->Print((outputFile).c_str(), "pdf portrait");
0405         gPad->SetLogy(0);
0406 
0407         // gPad->SetLogz(0);
0408         // h2ADC_HCALIN->SetMinimum(100);
0409         // c1->Print((outputFile).c_str(), "pdf portrait");
0410     }
0411 
0412 
0413     if(finput_HCALOUT) {
0414         c1->SetCanvasSize(1000, 1500);
0415         gPad->SetLogz();
0416         h2ADC_HCALOUT->SetStats(0);
0417         h2ADC_HCALOUT->SetTitle(("ADC HCALOUT, Run: " + run + ", event: " + event).c_str());
0418         h2ADC_HCALOUT->Draw("COLZ1");
0419         c1->Print((outputFile).c_str(), "pdf portrait");
0420 
0421         c1->SetCanvasSize(1500, 1000);
0422         gPad->SetLogy();
0423         hADC_HCALOUT->SetTitle(("ADC HCALOUT, Run: " + run + ", event: " + event).c_str());
0424         hADC_HCALOUT->Draw();
0425         c1->Print((outputFile).c_str(), "pdf portrait");
0426         gPad->SetLogy(0);
0427 
0428         // gPad->SetLogz(0);
0429         // h2ADC_HCALOUT->SetMinimum(100);
0430         // c1->Print((outputFile).c_str(), "pdf portrait");
0431     }
0432 
0433     c1->Print((outputFile + "]").c_str(), "pdf portrait");
0434 
0435     // Close files
0436     if(finput_CEMC)    finput_CEMC->Close();
0437     if(finput_HCALIN)  finput_HCALIN->Close();
0438     if(finput_HCALOUT) finput_HCALOUT->Close();
0439     output.close();
0440 }
0441 
0442 void test_hcal() {
0443     UInt_t max_eta = 0;
0444     UInt_t max_phi = 0;
0445     for (UInt_t channel = 0; channel < myAnalysis::channels_hcal_bins; ++channel) {
0446         UInt_t key    = TowerInfoDefs::encode_hcal(channel);
0447         UInt_t etabin = TowerInfoDefs::getCaloTowerEtaBin(key);
0448         UInt_t phibin = TowerInfoDefs::getCaloTowerPhiBin(key);
0449 
0450         cout << "channel: " << channel << ", etabin: " << etabin << ", phibin: " << phibin << endl;
0451 
0452         max_eta = max(max_eta, etabin);
0453         max_phi = max(max_phi, phibin);
0454     }
0455     cout << "max etabin: " << max_eta << ", max phibin: " << max_phi << endl;
0456 }
0457 
0458 void event_display(const string& run,
0459                    const string& event,
0460                    const string& input_CEMC,
0461                    const string& outputFile="event-display/test.json",
0462                    const string& input_HCALIN = "",
0463                    const string& input_HCALOUT = "") {
0464     myAnalysis::init(input_CEMC, outputFile, input_HCALIN, input_HCALOUT);
0465     myAnalysis::analyze();
0466     myAnalysis::finalize(run, event);
0467 }
0468 
0469 # ifndef __CINT__
0470 int main(int argc, char* argv[]) {
0471     if(argc < 1 || argc > 7){
0472         cout << "usage: ./bin/event-display run event input_CEMC outputFile input_HCALIN input_HCALOUT" << endl;
0473         cout << "run: run number." << endl;
0474         cout << "event: event number." << endl;
0475         cout << "input_CEMC: Location of CEMC root file." << endl;
0476         cout << "outputFile: output json file. Default = event-display/test.json." << endl;
0477         cout << "input_HCALIN: Location of HCALIN root file." << endl;
0478         cout << "input_HCALOUT: Location of HCALOUT root file." << endl;
0479         return 1;
0480     }
0481 
0482     string run;
0483     string event;
0484     string input_CEMC;
0485     string output = "event-display/test.json";
0486     string input_HCALIN = "";
0487     string input_HCALOUT = "";
0488 
0489     if(argc >= 2) {
0490         run = argv[1];
0491     }
0492     if(argc >= 3) {
0493         event = argv[2];
0494     }
0495     if(argc >= 4) {
0496         input_CEMC = argv[3];
0497     }
0498     if(argc >= 5) {
0499         output = argv[4];
0500     }
0501     if(argc >= 6) {
0502         input_HCALIN = argv[5];
0503     }
0504     if(argc >= 7) {
0505         input_HCALOUT = argv[6];
0506     }
0507     event_display(run, event, input_CEMC, output, input_HCALIN, input_HCALOUT);
0508 
0509     cout << "done" << endl;
0510     return 0;
0511 }
0512 # endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!