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 #include <TH1F.h>
0017 #include <TH2F.h>
0018
0019 #include <algorithm>
0020 #include <cmath>
0021 #include <iostream>
0022
0023 using namespace std;
0024
0025 HCalCosmics::HCalCosmics(const std::string &name, const std::string &fname)
0026   : SubsysReco(name)
0027   , outfilename(fname)
0028 {
0029 }
0030
0031 int HCalCosmics::Init(PHCompositeNode * /*topNode*/)
0032 {
0033   std::cout << std::endl
0034             << "HCalCosmics::Init" << std::endl;
0035   std::cout << "Node with raw ADC: " << rawprefix + rawdetector << std::endl;
0036   std::cout << "Node with calibrated Energy in eV: " << prefix + detector << std::endl;
0037   std::cout << "Threshold: " << tower_threshold << " " << vert_threshold << " " << veto_threshold << std::endl;
0038   std::cout << "Bin width: " << bin_width << std::endl;
0039
0040   outfile = new TFile(outfilename.c_str(), "RECREATE");
0041
0042   for (int ieta = 0; ieta < n_etabin; ++ieta)
0043   {
0044     for (int iphi = 0; iphi < n_phibin; ++iphi)
0045     {
0046       std::string channel_histname = "h_channel_" + std::to_string(ieta) + "_" + std::to_string(iphi);
0047       h_channel_hist[ieta][iphi] = new TH1F(channel_histname.c_str(), "", 500, 0, 500 * bin_width);
0048
0049       std::string adc_histname = "h_adc_" + std::to_string(ieta) + "_" + std::to_string(iphi);
0050       h_adc_hist[ieta][iphi] = new TH1F(adc_histname.c_str(), "", 500, 0, 16000 * rawbin_width);
0051
0052       std::string adc_ecut_histname = "h_adc_ecut_" + std::to_string(ieta) + "_" + std::to_string(iphi);
0053       h_adc_ecut_hist[ieta][iphi] = new TH1F(adc_ecut_histname.c_str(), "", 500, 0, 16000 * rawbin_width);
0054
0055
0056       std::string time_histname = "h_towertime_" + std::to_string(ieta) + "_" + std::to_string(iphi);
0057       h_towertime_hist[ieta][iphi] = new TH1F(time_histname.c_str(), "", 100, -10, 10);
0058
0059       std::string calibfactor_histname = "h_gain_" + std::to_string(ieta) + "_" + std::to_string(iphi);
0060       h_gain_hist[ieta][iphi] = new TH1F(calibfactor_histname.c_str(), "", 500, 0, 1000);
0061
0062     }
0063   }
0064   h_waveformchi2 = new TH2F("h_waveformchi2", "", 1000, 0, 500 * bin_width, 1000, 0, 1000000);
0065   h_waveformchi2->GetXaxis()->SetTitle("peak (ADC)");
0066   h_waveformchi2->GetYaxis()->SetTitle("chi2");
0067   h_waveformchi2_aftercut = new TH2F("h_waveformchi2_aftercut", "", 1000, 0, 500 * bin_width, 1000, 0, 1000000);
0068   h_waveformchi2_aftercut->GetXaxis()->SetTitle("peak (ADC)");
0069   h_waveformchi2_aftercut->GetYaxis()->SetTitle("chi2");
0070   h_mip = new TH1F("h_mip", "", 500, 0, 500 * bin_width);
0071   h_adc = new TH1F("h_adc", "", 500, 0, 16000 * rawbin_width);
0072   h_adc_ecut = new TH1F("h_adc_ecut", "", 500, 0, 16000 * rawbin_width);
0073   h_gain = new TH1F("h_gain", "", 500, -500, 1000);
0074   h_event = new TH1F("h_event", "", 1, 0, 1);
0075
0076   h_time_energy = new TH2F("h_time_energy", "", 100, -10, 10, 100, -10 * bin_width, 90 * bin_width);
0077
0078   event = 0;
0079   return 0;
0080 }
0081
0082 int HCalCosmics::process_event(PHCompositeNode *topNode)
0083 {
0084   process_towers(topNode);
0085   event++;
0086   h_event->Fill(0);
0087
0088   return Fun4AllReturnCodes::EVENT_OK;
0089 }
0090
0091 int HCalCosmics::process_towers(PHCompositeNode *topNode)
0092 {
0093
0094   //***************** Raw tower ADC (uncalibrated) info *****************
0095
0096   std::string rawnode_name = rawprefix + rawdetector;
0097
0098   TowerInfoContainer *rawtowers = findNode::getClass<TowerInfoContainer>(topNode, rawnode_name);
0099   if (!rawtowers)
0100   {
0101     std::cout << std::endl
0102               << "Didn't find node " << rawnode_name << std::endl;
0103     return Fun4AllReturnCodes::EVENT_OK;
0104   }
0105
0106   //*****************  tower energy in eV (calibrated) info ***************
0107
0108   std::string calibnode_name = prefix + detector;
0109
0110   TowerInfoContainer *towers = findNode::getClass<TowerInfoContainer>(topNode, calibnode_name);
0111   if (!towers)
0112   {
0113     std::cout << std::endl
0114               << "Didn't find node " << calibnode_name << std::endl;
0115     return Fun4AllReturnCodes::EVENT_OK;
0116   }
0117
0118   int size = towers->size(); // 1536 towers
0119
0120   for (int channel = 0; channel < size; channel++)
0121   {
0122     TowerInfo *tower = towers->get_tower_at_channel(channel);
0123     TowerInfo *rawtower = rawtowers->get_tower_at_channel(channel);
0124     float adc = rawtower->get_energy();
0125     float energy = tower->get_energy();
0126     float chi2 = tower->get_chi2();
0127     float time = tower->get_time_float();
0128     unsigned int towerkey = towers->encode_key(channel);
0129     int ieta = towers->getTowerEtaBin(towerkey);
0130     int iphi = towers->getTowerPhiBin(towerkey);
0131     m_peak[ieta][iphi] = energy;
0132     m_adc[ieta][iphi] = adc;
0133     m_chi2[ieta][iphi] = chi2;
0134     m_time[ieta][iphi] = time;
0135     h_waveformchi2->Fill(m_peak[ieta][iphi], m_chi2[ieta][iphi]);
0136     if (tower->get_isBadChi2())
0137     {
0138       m_peak[ieta][iphi] = 0;
0139     }
0140     h_waveformchi2_aftercut->Fill(m_peak[ieta][iphi], m_chi2[ieta][iphi]);
0141     h_time_energy->Fill(time, energy);
0142   }
0143   // Apply cuts based on calibrated energy
0144   for (int ieta = 0; ieta < n_etabin; ++ieta)
0145   {
0146     for (int iphi = 0; iphi < n_phibin; ++iphi)
0147     {
0148       if (m_peak[ieta][iphi] < tower_threshold)
0149       {
0150         continue;  // target tower cut
0151       }
0152       int up = iphi + 1;
0153       int down = iphi - 1;
0154       if (up > 63) { up -= 64; }
0155       if (down < 0) { down += 64; }
0156       if (m_peak[ieta][up] < vert_threshold || m_peak[ieta][down] < vert_threshold)
0157       {
0158         continue;  // vertical neighbor cut
0159       }
0160       if (ieta != 0 && (m_peak[ieta - 1][up] > veto_threshold ||
0161                         m_peak[ieta - 1][iphi] > veto_threshold ||
0162                         m_peak[ieta - 1][down] > veto_threshold))
0163       {
0164         continue;  // left veto cut
0165       }
0166       if (ieta != 23 && (m_peak[ieta + 1][up] > veto_threshold ||
0167                          m_peak[ieta + 1][iphi] > veto_threshold ||
0168                          m_peak[ieta + 1][down] > veto_threshold))
0169       {
0170         continue;  // right veto cut
0171       }
0172
0173       h_channel_hist[ieta][iphi]->Fill(m_peak[ieta][iphi]);
0174       h_towertime_hist[ieta][iphi]->Fill(m_time[ieta][iphi]);
0175       h_mip->Fill(m_peak[ieta][iphi]);
0176       h_adc_ecut_hist[ieta][iphi]->Fill(m_adc[ieta][iphi]);
0177       h_adc_ecut->Fill(m_adc[ieta][iphi]);
0178
0179       if (m_peak[ieta][iphi] != 0)
0180       {
0181          double gain_fac = m_adc[ieta][iphi] / m_peak[ieta][iphi];
0182          h_gain_hist[ieta][iphi]->Fill(gain_fac);
0183          h_gain->Fill(gain_fac);
0184       }
0185
0186
0187     }
0188   }
0189
0190
0191   for (int ieta = 0; ieta < n_etabin; ++ieta)
0192   {
0193     for (int iphi = 0; iphi < n_phibin; ++iphi)
0194     {
0195       if (m_adc[ieta][iphi] < adc_tower_threshold)
0196       {
0197         continue;   // target tower cut
0198       }
0199       int up = iphi + 1;
0200       int down = iphi - 1;
0201       if (up > 63) { up -= 64; }
0202       if (down < 0) { down += 64; }
0203
0204       if (m_adc[ieta][up] < adc_vert_threshold || m_adc[ieta][down] < adc_vert_threshold)
0205         {
0206           continue;  // vertical neighbor cut
0207         }
0208       if (ieta != 0 && (m_adc[ieta - 1][up] > adc_veto_threshold ||
0209                         m_adc[ieta - 1][iphi] > adc_veto_threshold ||
0210                         m_adc[ieta - 1][down] > adc_veto_threshold))
0211         {
0212           continue;  // left veto cut
0213         }
0214       if (ieta != 23 && (m_adc[ieta + 1][up] > adc_veto_threshold ||
0215                          m_adc[ieta + 1][iphi] > adc_veto_threshold ||
0216                          m_adc[ieta + 1][down] > adc_veto_threshold))
0217         {
0218           continue;  // right veto cut
0219         }
0220         h_adc_hist[ieta][iphi]->Fill(m_adc[ieta][iphi]);
0221         h_adc->Fill(m_adc[ieta][iphi]);
0222
0223
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 }
0399
0400