HCalCosmicCalib/offline/fit_channel_hist.C

0001 #include <TH1.h>
0002 #include <TStyle.h>
0003 #include <TCanvas.h>
0004 #include <TMath.h>
0005 #include <TF1.h>
0006 #include "TROOT.h"
0007 #include <iostream>
0008 #include <string>
0009 #include <fstream>
0010 #include <TAxis.h>
0011 #include <TGraph.h>
0012 #include <TSpectrum.h>
0013
0014 double gamma_function(double *x, double *par) {
0015   double peak = par[0];
0016   double shift = par[1];
0017   double scale = par[2];
0018   double N = par[3];
0019
0020   double arg_para = (x[0] - shift) / scale;
0021   double peak_para = (peak - shift) / scale;
0022   double numerator = N * pow(arg_para, peak_para) * TMath::Exp(-arg_para);
0023   double denominator = ROOT::Math::tgamma(peak_para + 1) * scale;
0024
0025   if (denominator != 0) return (double)(numerator / denominator);
0026   else return 0;
0027 }
0028
0029 float get_temp_shift(TH1F* hist, float y) {
0030   float shift = 0.;
0031   for (int i = 1; i <= hist->GetNbinsX(); i++) {
0032     if (hist->GetBinContent(i) > y/10.) {
0033       shift = hist->GetBinCenter(i);
0034       break;
0035     }
0036   }
0037   return shift;
0038 }
0039
0040 float get_temp_scale(TH1F* hist, float x) {
0041   float scale = (hist->GetMean() - x);
0042   return scale;
0043 }
0044
0045 float get_temp_N(float x, float y, float shift, float scale) {
0046   float temp_par0 = (x - shift) / scale;
0047   float temp_par1 = pow(temp_par0, temp_par0) * TMath::Exp(-temp_par0);
0048   float temp_par2 = ROOT::Math::tgamma(temp_par0 + 1) * scale * y;
0049   float temp_N = temp_par2 / temp_par1;
0050   return temp_N;
0051 }
0052
0053 float get_fit_min(TH1F* hist, float y) {
0054   float min = 0.;
0055   for (int i = 1; i <= hist->GetNbinsX(); i++) {
0056     if (hist->GetBinContent(i) > 2*y/3.) {
0057       min = hist->GetBinCenter(i);
0058       if (hist->GetBinContent(i) > 3*y/4.) {
0059         min = hist->GetBinCenter(i-1);
0060       }
0061       break;
0062     }
0063   }
0064   return min;
0065 }
0066
0067 float get_fit_max(TH1F* hist, float y) {
0068   float max = 0.;
0069   for (int i = hist->GetNbinsX(); i >= 1; i--) {
0070     if (hist->GetBinContent(i) > y/5.) {
0071       max = hist->GetBinCenter(i);
0072       break;
0073     }
0074   }
0075   return max;
0076 }
0077
0078 void do_peak_finder(TH1F* hist, float& x, float& y, float& shift, float& scale, float& N, float& min, float& max) {
0079   TH1F *h_peakfinder = (TH1F*)hist->Clone("h_peakfinder");
0080   h_peakfinder->Smooth(2);
0081   TSpectrum *s = new TSpectrum(10);
0082   int nfound = 0;
0083   nfound = s->Search(h_peakfinder, 2, "nobackground new", 0.6);
0084   if (nfound == 0) {
0085     int index = hist->GetMaximumBin();
0086     x = hist->GetBinCenter(index);
0087     y = hist->GetBinContent(index);
0088   } else {
0089     double *xtemp = s->GetPositionX();
0090     double *ytemp = s->GetPositionY();
0091     x = xtemp[0];
0092     y = ytemp[0];
0093   }
0094   shift = get_temp_shift(h_peakfinder, y);
0095   scale = get_temp_scale(hist, x);
0096   N = get_temp_N(x, y, shift, scale);
0097   min = get_fit_min(h_peakfinder, y);
0098   max = get_fit_max(h_peakfinder, y);
0099   delete h_peakfinder;
0100   delete s;
0101 }
0102
0103 void do_peak_fit(TH1F* hist, TF1* func) {
0104   float peak_position_finder, peak_value_finder, shift_finder, scale_finder, N_finder, min_finder, max_finder;
0105   do_peak_finder(hist, peak_position_finder, peak_value_finder, shift_finder, scale_finder, N_finder, min_finder, max_finder);
0106   func->SetParameter(0, peak_position_finder);
0107   func->SetParameter(1, shift_finder);
0108   func->SetParameter(2, scale_finder);
0109   func->SetParameter(3, N_finder);
0110   hist->Fit(func, "", "", min_finder, max_finder);
0111 }
0112
0113 void fit_channel_hist() {
0114   // Constant
0115   static const int n_etabin = 24;
0116   static const int n_phibin = 64;
0117
0118   // File
0119   TFile *ohcal_input = new TFile("../macro/ohcal_hist.root","READ");
0120   TFile *ihcal_input = new TFile("../macro/ihcal_hist.root","READ");
0121   std::ofstream total_mpvinfo("total_mpvinfo.txt",ios::trunc);
0122   TFile *ohcal_output = new TFile("ohcal_fitresult.root","recreate");
0123   std::ofstream ohcal_mpvinfo("ohcal_mpvinfo.txt",ios::trunc);
0124   std::ofstream ohcal_calibinfo("ohcal_calibinfo.txt",ios::trunc);
0125   TFile *ihcal_output = new TFile("ihcal_fitresult.root","recreate");
0126   std::ofstream ihcal_mpvinfo("ihcal_mpvinfo.txt",ios::trunc);
0127   std::ofstream ihcal_calibinfo("ihcal_calibinfo.txt",ios::trunc);
0128
0129   // oHCal setup
0130   TH2F *h_2Dohcal_hit = new TH2F("h_2Dohcal_hit", "", n_etabin, 0, n_etabin, n_phibin, 0, n_phibin);
0131   h_2Dohcal_hit->GetXaxis()->SetTitle("ieta");
0132   h_2Dohcal_hit->GetYaxis()->SetTitle("iphi");
0133   TH1F *h_ohcal_total;
0134   h_ohcal_total = new TH1F("h_ohcal_total","",200,0,10000);
0135   h_ohcal_total->GetXaxis()->SetTitle("ADC");
0136   h_ohcal_total->GetYaxis()->SetTitle("Counts");
0137   TH1F *h_ohcal_channel[n_etabin][n_phibin];
0138   for (int ieta = 0; ieta < n_etabin; ++ieta) {
0139     for (int iphi = 0; iphi < n_phibin; ++iphi) {
0140       std::ostringstream h_ohcal_channelname;
0141       h_ohcal_channelname << "h_ohcal_channel_" << ieta << "_" << iphi;
0142       h_ohcal_channel[ieta][iphi] = new TH1F(h_ohcal_channelname.str().c_str(), "", 200, 0, 10000);
0143       h_ohcal_channel[ieta][iphi]->GetXaxis()->SetTitle("ADC");
0144       h_ohcal_channel[ieta][iphi]->GetYaxis()->SetTitle("Counts");
0145     }
0146   }
0147   TH1F *h_ohcal_chindf = new TH1F("h_ohcal_chindf","",50,0,10);
0148   h_ohcal_chindf->GetXaxis()->SetTitle("Chi^2/NDF");
0149   h_ohcal_chindf->GetYaxis()->SetTitle("Counts");
0150
0151   // iHCal setup
0152   TH2F *h_2Dihcal_hit = new TH2F("h_2Dihcal_hit", "", n_etabin, 0, n_etabin, n_phibin, 0, n_phibin);
0153   h_2Dihcal_hit->GetXaxis()->SetTitle("ieta");
0154   h_2Dihcal_hit->GetYaxis()->SetTitle("iphi");
0155   TH1F *h_ihcal_total;
0156   h_ihcal_total = new TH1F("h_ihcal_total","",200,0,10000);
0157   h_ihcal_total->GetXaxis()->SetTitle("ADC");
0158   h_ihcal_total->GetYaxis()->SetTitle("Counts");
0159   TH1F *h_ihcal_channel[n_etabin][n_phibin];
0160   for (int ieta = 0; ieta < n_etabin; ++ieta) {
0161     for (int iphi = 0; iphi < n_phibin; ++iphi) {
0162       std::ostringstream h_ihcal_channelname;
0163       h_ihcal_channelname << "h_ihcal_channel_" << ieta << "_" << iphi;
0164       h_ihcal_channel[ieta][iphi] = new TH1F(h_ihcal_channelname.str().c_str(), "", 200, 0, 10000);
0165       h_ihcal_channel[ieta][iphi]->GetXaxis()->SetTitle("ADC");
0166       h_ihcal_channel[ieta][iphi]->GetYaxis()->SetTitle("Counts");
0167     }
0168   }
0169   TH1F *h_ihcal_chindf = new TH1F("h_ihcal_chindf","",50,0,10);
0170   h_ihcal_chindf->GetXaxis()->SetTitle("Chi^2/NDF");
0171   h_ihcal_chindf->GetYaxis()->SetTitle("Counts");
0172
0173   // Read hist
0174   TH1F *h_temp_channel;
0175   for (int ieta = 0; ieta < n_etabin; ++ieta) {
0176     for (int iphi = 0; iphi < n_phibin; ++iphi) {
0177       std::ostringstream h_temp_channelname;
0178       h_temp_channelname << "h_channel_" << ieta << "_" << iphi;
0179       h_temp_channel = (TH1F*)ohcal_input->Get(h_temp_channelname.str().c_str());
0180       h_ohcal_channel[ieta][iphi]->Add(h_temp_channel);
0181       h_ohcal_total->Add(h_temp_channel);
0182       h_temp_channel = (TH1F*)ihcal_input->Get(h_temp_channelname.str().c_str());
0183       h_ihcal_channel[ieta][iphi]->Add(h_temp_channel);
0184       h_ihcal_total->Add(h_temp_channel);
0185     }
0186   }
0187   ohcal_input->Close();
0188   ihcal_input->Close();
0189
0190   // Fit hist.
0191   float peak_ohcal_total, peak_ihcal_total, chindf_ohcal_total, chindf_ihcal_total;
0192   float peak_ohcal[n_etabin][n_phibin], peak_ihcal[n_etabin][n_phibin];
0193   float chindf_ohcal[n_etabin][n_phibin], chindf_ihcal[n_etabin][n_phibin];
0194   TF1 *f_gamma = new TF1("f_gamma", gamma_function, 0, 10000, 4);
0195   f_gamma->SetParName(0,"Peak(ADC)");
0196   f_gamma->SetParName(1,"Shift");
0197   f_gamma->SetParName(2,"Scale");
0198   f_gamma->SetParName(3,"N");
0199
0200   do_peak_fit(h_ohcal_total, f_gamma);
0201   peak_ohcal_total = f_gamma->GetParameter(0);
0202   chindf_ohcal_total = f_gamma->GetChisquare() / (float)f_gamma->GetNDF();
0203
0204   do_peak_fit(h_ihcal_total, f_gamma);
0205   peak_ihcal_total = f_gamma->GetParameter(0);
0206   chindf_ihcal_total = f_gamma->GetChisquare() / (float)f_gamma->GetNDF();
0207   for (int ieta = 0; ieta < n_etabin; ++ieta) {
0208     for (int iphi = 0; iphi < n_phibin; ++iphi) {
0209       do_peak_fit(h_ohcal_channel[ieta][iphi], f_gamma);
0210       peak_ohcal[ieta][iphi] = f_gamma->GetParameter(0);
0211       chindf_ohcal[ieta][iphi] = f_gamma->GetChisquare() / (float)f_gamma->GetNDF();
0212       h_ohcal_chindf->Fill(chindf_ohcal[ieta][iphi]);
0213       do_peak_fit(h_ihcal_channel[ieta][iphi], f_gamma);
0214       peak_ihcal[ieta][iphi] = f_gamma->GetParameter(0);
0215       chindf_ihcal[ieta][iphi] = f_gamma->GetChisquare() / (float)f_gamma->GetNDF();
0216       h_ihcal_chindf->Fill(chindf_ihcal[ieta][iphi]);
0217
0218       h_2Dohcal_hit->SetBinContent(ieta+1, iphi+1, h_ohcal_channel[ieta][iphi]->GetEntries());
0219       h_2Dihcal_hit->SetBinContent(ieta+1, iphi+1, h_ihcal_channel[ieta][iphi]->GetEntries());
0220     }
0221   }
0222   delete f_gamma;
0223
0224   // Output.
0225   // Get Calib info.
0226   float ihcal_samplingfraction = 0.162166;
0227   float ohcal_samplingfraction = 3.38021e-02;
0228   float ihcal_mipcalib, ohcal_mipcalib, ihcal_abscalib, ohcal_abscalib;
0229   for (int ieta = 0; ieta < n_etabin; ++ieta) {
0230     for (int iphi = 0; iphi < n_phibin; ++iphi) {
0231       ihcal_mipcalib = peak_ihcal[ieta][iphi] / 32.;
0232       ohcal_mipcalib = peak_ohcal[ieta][iphi] / 32.;
0233       ihcal_abscalib = 1. / ihcal_mipcalib / ihcal_samplingfraction; // abs calib in GeV per ADC
0234       ohcal_abscalib = 1. / ohcal_mipcalib / ohcal_samplingfraction; // abs calib in GeV per ADC
0235       ihcal_calibinfo << ieta << "\t" << iphi << "\t" << ihcal_abscalib << std::endl;
0236       ohcal_calibinfo << ieta << "\t" << iphi << "\t" << ohcal_abscalib << std::endl;
0237     }
0238   }
0239
0240   total_mpvinfo << "oHCal: MPV = " << peak_ohcal_total << ", Chi^2/NDF = " << chindf_ohcal_total << std::endl;
0241   total_mpvinfo << "iHCal: MPV = " << peak_ihcal_total << ", Chi^2/NDF = " << chindf_ihcal_total << std::endl;
0242   total_mpvinfo << std::endl << "Channel-by-channel Chi^2/NDF value: " << std::endl << "    oHCal:" << h_ohcal_chindf->GetMean() << std::endl << "    iHCal:" << h_ihcal_chindf->GetMean() << std::endl;
0243   std::vector<int> badfit_eta_ohcal, badfit_phi_ohcal, badfit_eta_ihcal, badfit_phi_ihcal;
0244   for (int ieta = 0; ieta < n_etabin; ++ieta) {
0245     for (int iphi = 0; iphi < n_phibin; ++iphi) {
0246       if (chindf_ohcal[ieta][iphi] > 2) {
0247         badfit_eta_ohcal.push_back(ieta);
0248         badfit_phi_ohcal.push_back(iphi);
0249       }
0250       if (chindf_ihcal[ieta][iphi] > 2) {
0251         badfit_eta_ihcal.push_back(ieta);
0252         badfit_phi_ihcal.push_back(iphi);
0253       }
0254       ohcal_output->cd();
0255       h_ohcal_channel[ieta][iphi]->Write();
0256       delete h_ohcal_channel[ieta][iphi];
0257       ihcal_output->cd();
0258       h_ihcal_channel[ieta][iphi]->Write();
0259       delete h_ihcal_channel[ieta][iphi];
0260       ohcal_mpvinfo << ieta << "\t" << iphi << "\t" << peak_ohcal[ieta][iphi] << "\t" << chindf_ohcal[ieta][iphi] << std::endl;
0261       ihcal_mpvinfo << ieta << "\t" << iphi << "\t" << peak_ihcal[ieta][iphi] << "\t" << chindf_ihcal[ieta][iphi] << std::endl;
0262     }
0263   }
0264   ohcal_mpvinfo.close();
0265   ihcal_mpvinfo.close();
0266
0267   total_mpvinfo << std::endl << "oHCal bad fit channel:" << std::endl;
0268   for (int n_bad = 0; n_bad < badfit_eta_ohcal.size(); ++n_bad) {
0269     total_mpvinfo << badfit_eta_ohcal[n_bad] << "\t" << badfit_phi_ohcal[n_bad] << "\t" << chindf_ohcal[badfit_eta_ohcal[n_bad]][badfit_phi_ohcal[n_bad]];
0270   }
0271   total_mpvinfo << std::endl << "iHCal bad fit channel:" << std::endl;
0272   for (int n_bad = 0; n_bad < badfit_eta_ihcal.size(); ++n_bad) {
0273     total_mpvinfo << badfit_eta_ihcal[n_bad] << "\t" << badfit_phi_ihcal[n_bad] << "\t" << chindf_ihcal[badfit_eta_ihcal[n_bad]][badfit_phi_ihcal[n_bad]];
0274   }
0275   total_mpvinfo.close();
0276
0277   ohcal_output->cd();
0278   h_ohcal_total->Write();
0279   h_ohcal_chindf->Write();
0280   h_2Dohcal_hit->Write();
0281   delete h_ohcal_total;
0282   delete h_ohcal_chindf;
0283   ohcal_output->Close();
0284
0285   ihcal_output->cd();
0286   h_ihcal_total->Write();
0287   h_ihcal_chindf->Write();
0288   h_2Dihcal_hit->Write();
0289   delete h_ihcal_total;
0290   delete h_ihcal_chindf;
0291   ihcal_output->Close();
0292 }