calorimeter/calo_tower_slope/HCalCosmics.cc

0001 #include "HCalCosmics.h"
0002
0003 #include <calobase/TowerInfo.h>
0004 #include <calobase/TowerInfoContainer.h>
0005
0006 #include <fun4all/Fun4AllReturnCodes.h>
0007 #include <fun4all/SubsysReco.h>
0008
0009 #include <phool/getClass.h>
0010
0011 #include <Math/SpecFuncMathCore.h>
0012 #include <TF1.h>
0013 #include <TFile.h>
0014 #include <TH1.h>
0015 #include <TH2.h>
0016
0017 #include <algorithm>
0018 #include <cmath>
0019 #include <iostream>
0020
0021 HCalCosmics::HCalCosmics(const std::string &name, const std::string &fname)
0022   : SubsysReco(name)
0023   , outfilename(fname)
0024 {
0025 }
0026
0027 int HCalCosmics::Init(PHCompositeNode * /*topNode*/)
0028 {
0029   std::cout << std::endl
0030             << "HCalCosmics::Init" << std::endl;
0031   std::cout << "Node with raw ADC: " << rawprefix + rawdetector << std::endl;
0032   std::cout << "Node with calibrated Energy in eV: " << prefix + detector << std::endl;
0033   std::cout << "Threshold: " << tower_threshold << " " << vert_threshold << " " << veto_threshold << std::endl;
0034   std::cout << "Bin width: " << bin_width << std::endl;
0035
0036   outfile = new TFile(outfilename.c_str(), "RECREATE");
0037
0038   for (int ieta = 0; ieta < n_etabin; ++ieta)
0039   {
0040     for (int iphi = 0; iphi < n_phibin; ++iphi)
0041     {
0042       std::string channel_histname = "h_channel_" + std::to_string(ieta) + "_" + std::to_string(iphi);
0043       h_channel_hist[ieta][iphi] = new TH1F(channel_histname.c_str(), "", 500, 0, 500 * bin_width);
0044
0045       std::string adc_histname = "h_adc_" + std::to_string(ieta) + "_" + std::to_string(iphi);
0046       h_adc_hist[ieta][iphi] = new TH1F(adc_histname.c_str(), "", 500, 0, 16000 * rawbin_width);
0047
0048       std::string adc_ecut_histname = "h_adc_ecut_" + std::to_string(ieta) + "_" + std::to_string(iphi);
0049       h_adc_ecut_hist[ieta][iphi] = new TH1F(adc_ecut_histname.c_str(), "", 500, 0, 16000 * rawbin_width);
0050
0051       std::string time_histname = "h_towertime_" + std::to_string(ieta) + "_" + std::to_string(iphi);
0052       h_towertime_hist[ieta][iphi] = new TH1F(time_histname.c_str(), "", 100, -10, 10);
0053
0054       std::string calibfactor_histname = "h_gain_" + std::to_string(ieta) + "_" + std::to_string(iphi);
0055       h_gain_hist[ieta][iphi] = new TH1F(calibfactor_histname.c_str(), "", 500, 0, 1000);
0056     }
0057   }
0058   h_waveformchi2 = new TH2F("h_waveformchi2", "", 1000, 0, 500 * bin_width, 1000, 0, 1000000);
0059   h_waveformchi2->GetXaxis()->SetTitle("peak (ADC)");
0060   h_waveformchi2->GetYaxis()->SetTitle("chi2");
0061   h_waveformchi2_aftercut = new TH2F("h_waveformchi2_aftercut", "", 1000, 0, 500 * bin_width, 1000, 0, 1000000);
0062   h_waveformchi2_aftercut->GetXaxis()->SetTitle("peak (ADC)");
0063   h_waveformchi2_aftercut->GetYaxis()->SetTitle("chi2");
0064   h_mip = new TH1F("h_mip", "", 500, 0, 500 * bin_width);
0065   h_adc = new TH1F("h_adc", "", 500, 0, 16000 * rawbin_width);
0066   h_adc_ecut = new TH1F("h_adc_ecut", "", 500, 0, 16000 * rawbin_width);
0067   h_gain = new TH1F("h_gain", "", 500, -500, 1000);
0068   h_event = new TH1F("h_event", "", 1, 0, 1);
0069
0070   h_time_energy = new TH2F("h_time_energy", "", 100, -10, 10, 100, -10 * bin_width, 90 * bin_width);
0071
0072   event = 0;
0073   return 0;
0074 }
0075
0076 int HCalCosmics::process_event(PHCompositeNode *topNode)
0077 {
0078   process_towers(topNode);
0079   event++;
0080   h_event->Fill(0);
0081
0082   return Fun4AllReturnCodes::EVENT_OK;
0083 }
0084
0085 int HCalCosmics::process_towers(PHCompositeNode *topNode)
0086 {
0087   //***************** Raw tower ADC (uncalibrated) info *****************
0088
0089   std::string rawnode_name = rawprefix + rawdetector;
0090
0091   TowerInfoContainer *rawtowers = findNode::getClass<TowerInfoContainer>(topNode, rawnode_name);
0092   if (!rawtowers)
0093   {
0094     std::cout << std::endl
0095               << "Didn't find node " << rawnode_name << std::endl;
0096     return Fun4AllReturnCodes::EVENT_OK;
0097   }
0098
0099   //*****************  tower energy in eV (calibrated) info ***************
0100
0101   std::string calibnode_name = prefix + detector;
0102
0103   TowerInfoContainer *towers = findNode::getClass<TowerInfoContainer>(topNode, calibnode_name);
0104   if (!towers)
0105   {
0106     std::cout << std::endl
0107               << "Didn't find node " << calibnode_name << std::endl;
0108     return Fun4AllReturnCodes::EVENT_OK;
0109   }
0110
0111   int size = towers->size();  // 1536 towers
0112
0113   for (int channel = 0; channel < size; channel++)
0114   {
0115     TowerInfo *tower = towers->get_tower_at_channel(channel);
0116     TowerInfo *rawtower = rawtowers->get_tower_at_channel(channel);
0117     float adc = rawtower->get_energy();
0118     float energy = tower->get_energy();
0119     float chi2 = tower->get_chi2();
0120     float time = tower->get_time();
0121     unsigned int towerkey = towers->encode_key(channel);
0122     int ieta = towers->getTowerEtaBin(towerkey);
0123     int iphi = towers->getTowerPhiBin(towerkey);
0124     m_peak[ieta][iphi] = energy;
0125     m_adc[ieta][iphi] = adc;
0126     m_chi2[ieta][iphi] = chi2;
0127     m_time[ieta][iphi] = time;
0128     h_waveformchi2->Fill(m_peak[ieta][iphi], m_chi2[ieta][iphi]);
0129     if (tower->get_isBadChi2())
0130     {
0131       m_peak[ieta][iphi] = 0;
0132     }
0133     h_waveformchi2_aftercut->Fill(m_peak[ieta][iphi], m_chi2[ieta][iphi]);
0134     h_time_energy->Fill(time, energy);
0135   }
0136   // Apply cuts based on calibrated energy
0137   for (int ieta = 0; ieta < n_etabin; ++ieta)
0138   {
0139     for (int iphi = 0; iphi < n_phibin; ++iphi)
0140     {
0141       if (m_peak[ieta][iphi] < tower_threshold)
0142       {
0143         continue;  // target tower cut
0144       }
0145       int up = iphi + 1;
0146       int down = iphi - 1;
0147       if (up > 63)
0148       {
0149         up -= 64;
0150       }
0151       if (down < 0)
0152       {
0153         down += 64;
0154       }
0155       if (m_peak[ieta][up] < vert_threshold || m_peak[ieta][down] < vert_threshold)
0156       {
0157         continue;  // vertical neighbor cut
0158       }
0159       if (ieta != 0 && (m_peak[ieta - 1][up] > veto_threshold ||
0160                         m_peak[ieta - 1][iphi] > veto_threshold ||
0161                         m_peak[ieta - 1][down] > veto_threshold))
0162       {
0163         continue;  // left veto cut
0164       }
0165       if (ieta != 23 && (m_peak[ieta + 1][up] > veto_threshold ||
0166                          m_peak[ieta + 1][iphi] > veto_threshold ||
0167                          m_peak[ieta + 1][down] > veto_threshold))
0168       {
0169         continue;  // right veto cut
0170       }
0171
0172       h_channel_hist[ieta][iphi]->Fill(m_peak[ieta][iphi]);
0173       h_towertime_hist[ieta][iphi]->Fill(m_time[ieta][iphi]);
0174       h_mip->Fill(m_peak[ieta][iphi]);
0175       h_adc_ecut_hist[ieta][iphi]->Fill(m_adc[ieta][iphi]);
0176       h_adc_ecut->Fill(m_adc[ieta][iphi]);
0177
0178       if (m_peak[ieta][iphi] != 0)
0179       {
0180         double gain_fac = m_adc[ieta][iphi] / m_peak[ieta][iphi];
0181         h_gain_hist[ieta][iphi]->Fill(gain_fac);
0182         h_gain->Fill(gain_fac);
0183       }
0184     }
0185   }
0186
0187   for (int ieta = 0; ieta < n_etabin; ++ieta)
0188   {
0189     for (int iphi = 0; iphi < n_phibin; ++iphi)
0190     {
0191       if (m_adc[ieta][iphi] < adc_tower_threshold)
0192       {
0193         continue;  // target tower cut
0194       }
0195       int up = iphi + 1;
0196       int down = iphi - 1;
0197       if (up > 63)
0198       {
0199         up -= 64;
0200       }
0201       if (down < 0)
0202       {
0203         down += 64;
0204       }
0205
0206       if (m_adc[ieta][up] < adc_vert_threshold || m_adc[ieta][down] < adc_vert_threshold)
0207       {
0208         continue;  // vertical neighbor cut
0209       }
0210       if (ieta != 0 && (m_adc[ieta - 1][up] > adc_veto_threshold ||
0211                         m_adc[ieta - 1][iphi] > adc_veto_threshold ||
0212                         m_adc[ieta - 1][down] > adc_veto_threshold))
0213       {
0214         continue;  // left veto cut
0215       }
0216       if (ieta != 23 && (m_adc[ieta + 1][up] > adc_veto_threshold ||
0217                          m_adc[ieta + 1][iphi] > adc_veto_threshold ||
0218                          m_adc[ieta + 1][down] > adc_veto_threshold))
0219       {
0220         continue;  // right veto cut
0221       }
0222       h_adc_hist[ieta][iphi]->Fill(m_adc[ieta][iphi]);
0223       h_adc->Fill(m_adc[ieta][iphi]);
0224     }
0225   }
0226
0227   return Fun4AllReturnCodes::EVENT_OK;
0228 }
0229
0230 int HCalCosmics::End(PHCompositeNode * /*topNode*/)
0231 {
0232   std::cout << "HCalCosmics::End" << std::endl;
0233   outfile->cd();
0234   for (auto &ieta : h_channel_hist)
0235   {
0236     for (auto &iphi : ieta)
0237     {
0238       iphi->Write();
0239       delete iphi;
0240     }
0241   }
0242
0243   for (auto &ieta : h_towertime_hist)
0244   {
0245     for (auto &iphi : ieta)
0246     {
0247       iphi->Write();
0248       delete iphi;
0249     }
0250   }
0251
0252   for (auto &ieta : h_adc_hist)
0253   {
0254     for (auto &iphi : ieta)
0255     {
0256       iphi->Write();
0257       delete iphi;
0258     }
0259   }
0260
0261   for (auto &ieta : h_adc_ecut_hist)
0262   {
0263     for (auto &iphi : ieta)
0264     {
0265       iphi->Write();
0266       delete iphi;
0267     }
0268   }
0269
0270   for (auto &ieta : h_gain_hist)
0271   {
0272     for (auto &iphi : ieta)
0273     {
0274       iphi->Write();
0275       delete iphi;
0276     }
0277   }
0278
0279   h_mip->Write();
0280   h_adc->Write();
0281   h_adc_ecut->Write();
0282   h_gain->Write();
0283   h_waveformchi2->Write();
0284   h_waveformchi2_aftercut->Write();
0285   h_time_energy->Write();
0286   h_event->Write();
0287
0288   outfile->Close();
0289   delete outfile;
0290   return 0;
0291 }
0292
0293 double HCalCosmics::gamma_function(const double *x, const double *par)
0294 {
0295   double peak = par[0];
0296   double shift = par[1];
0297   double scale = par[2];
0298   double N = par[3];
0299
0300   if (scale == 0)
0301   {
0302     return 0;
0303   }
0304
0305   double arg_para = (x[0] - shift) / scale;
0306   arg_para = std::max<double>(arg_para, 0);
0307   double peak_para = (peak - shift) / scale;
0308   //  double numerator = N * pow(arg_para, peak_para) * TMath::Exp(-arg_para);
0309   double numerator = N * pow(arg_para, peak_para) * std::exp(-arg_para);
0310   double denominator = ROOT::Math::tgamma(peak_para + 1) * scale;
0311
0312   if (denominator == 0)
0313   {
0314     return 1e8;
0315   }
0316   double val = numerator / denominator;
0317   if (std::isnan(val))
0318   {
0319     return 0;
0320   }
0321   return val;
0322 }
0323
0324 TF1 *HCalCosmics::fitHist(TH1 *h)
0325 {
0326   TF1 *f_gaus = new TF1("f_gaus", "gaus", 0, 10000);
0327   h->Fit(f_gaus, "QN", "", 0, 10000);
0328
0329   TF1 *f_gamma = new TF1("f_gamma", gamma_function, 0, 5000, 4);
0330   f_gamma->SetParName(0, "Peak(ADC)");
0331   f_gamma->SetParName(1, "Shift");
0332   f_gamma->SetParName(2, "Scale");
0333   f_gamma->SetParName(3, "N");
0334
0335   f_gamma->SetParLimits(0, 1000, 4000);
0336   f_gamma->SetParLimits(1, 500, 2000);
0337   f_gamma->SetParLimits(2, 200, 1000);
0338   f_gamma->SetParLimits(3, 0, 1e9);
0339
0340   f_gamma->SetParameter(0, f_gaus->GetParameter(1));
0341   f_gamma->SetParameter(1, 1300);
0342   f_gamma->SetParameter(2, 1300);
0343   f_gamma->SetParameter(3, 1e6);
0344
0345   h->Fit(f_gamma, "RQN", "", 1000, 5000);
0346
0347   return f_gamma;
0348 }
0349
0350 void HCalCosmics::fitChannels(const std::string &infile, const std::string &outfilename2)
0351 {
0352   TFile *fin = new TFile(infile.c_str(), "READ");
0353   if (!fin)
0354   {
0355     std::cout << "file " << infile << "   not found";
0356     return;
0357   }
0358   for (int ieta = 0; ieta < n_etabin; ++ieta)
0359   {
0360     for (int iphi = 0; iphi < n_phibin; ++iphi)
0361     {
0362       std::string channel_histname = "h_channel_" + std::to_string(ieta) + "_" + std::to_string(iphi);
0363       fin->GetObject(channel_histname.c_str(), h_channel_hist[ieta][iphi]);
0364     }
0365   }
0366
0367   if (!h_channel_hist[0][0])
0368   {
0369     std::cout << "no hists in " << infile << std::endl;
0370     return;
0371   }
0372
0373   TFile *outfileFit = new TFile(outfilename2.c_str(), "recreate");
0374
0375   TH2 *h2_peak = new TH2F("h2_peak", "", n_etabin, 0, n_etabin, n_phibin, 0, n_phibin);
0376
0377   TH1 *h_allTow = (TH1 *) h_channel_hist[0][0]->Clone("h_allTow");
0378   h_allTow->Reset();
0379
0380   for (int ieta = 0; ieta < n_etabin; ++ieta)
0381   {
0382     for (int iphi = 0; iphi < n_phibin; ++iphi)
0383     {
0384       TF1 *res = fitHist(h_channel_hist[ieta][iphi]);
0385       h2_peak->SetBinContent(ieta + 1, iphi + 1, res->GetParameter(0));
0386       h_allTow->Add(h_channel_hist[ieta][iphi]);
0387     }
0388   }
0389
0390   TF1 *res = fitHist(h_allTow);
0391   res->Write("f_allTow");
0392
0393   h_allTow->Write();
0394   h2_peak->Write();
0395
0396   outfileFit->Close();
0397   fin->Close();
0398 }