calo/emcal_calib_year1/doTscFit.C

0001 // -- c++ includes --
0002 #include <format>
0003 #include <memory>
0004 #include <string>
0005
0006 // -- root includes --
0007 #include <TFile.h>
0008 #include <TH2F.h>
0009 #include <TSystem.h>
0010
0011 #include <cdbobjects/CDBTTree.h>  // for CDBTTree
0012 #include <litecaloeval/LiteCaloEval.h>
0013 // TowerInfo
0014 #include <calobase/TowerInfoDefs.h>
0015
0016 #include <sPhenixStyle.C>
0017
0018 R__LOAD_LIBRARY(libLiteCaloEvalTowSlope.so)
0019 R__LOAD_LIBRARY(libcdbobjects)
0020
0021 void TSCtoCDBTTree(const char *infile, const char *outputfile, const std::string &m_fieldname);
0022 void mergeCDBTTrees(const char *infile1, const char *infile2, const char *outputfile, const std::string &m_fieldname);
0023
0024 void doTscFit(const std::string &hist_fname = "base/combine_out/out1.root",
0025               const std::string &calib_fname = "base/local_calib_copy.root",
0026               int iter = 1,
0027               const std::string &m_fieldname = "CEMC_calib_ADC_to_ETower")
0028 {
0029   std::string fitoutfile = std::format("tsc_fitout_it{}.root", iter);
0030
0031   if (iter <= 2)
0032   {
0033     LiteCaloEval modlce;
0034     modlce.CaloType(LiteCaloEval::CEMC);
0035     modlce.set_spectra_binWidth(0.02);
0036     modlce.Get_Histos(hist_fname.c_str(), fitoutfile.c_str());
0037     modlce.m_myminbin = 0;
0038     modlce.m_mymaxbin = 96;
0039     modlce.setFitMin(0.40F);
0040     modlce.setFitMax(1.2F);
0041     modlce.set_doQA();
0042     // if (iter==2) modlce.set_doQA(true);
0043     modlce.FitRelativeShifts(&modlce, 110);
0044   }
0045
0046   if (iter == 3)
0047   {
0048     SetsPhenixStyle();
0049     LiteCaloEval modlce;
0050     modlce.CaloType(LiteCaloEval::CEMC);
0051     modlce.set_spectra_binWidth(0.005);
0052     modlce.Get_Histos(hist_fname.c_str(), fitoutfile.c_str());
0053     modlce.plot_cemc("figures");
0054   }
0055
0056   // create the cdbttree from tsc output andd multiply the corrections
0057   // into the base calibration to pickup for pi0 first iteration
0058
0059   TSCtoCDBTTree(fitoutfile.c_str(), std::format("tsc_output_cdb_it{}.root", iter).c_str(), m_fieldname);
0060   mergeCDBTTrees(std::format("tsc_output_cdb_it{}.root", iter).c_str(), calib_fname.c_str(), calib_fname.c_str(), m_fieldname);
0061
0062   size_t pos = calib_fname.find_last_of('.');
0063   std::string f_calib_save_name = calib_fname;
0064   f_calib_save_name.insert(pos, std::format("_postTSC_it{}", iter));
0065
0066   std::unique_ptr<TFile> f_calib_mod = std::make_unique<TFile>(calib_fname.c_str());
0067   f_calib_mod->Cp(f_calib_save_name.c_str());
0068
0069   gSystem->Exit(0);
0070 }
0071
0072 void mergeCDBTTrees(const char *infile1, const char *infile2, const char *outputfile, const std::string &m_fieldname)
0073 {
0074   std::unique_ptr<CDBTTree> cdbttree1 = std::make_unique<CDBTTree>(infile1);
0075   std::unique_ptr<CDBTTree> cdbttree2 = std::make_unique<CDBTTree>(infile2);
0076   std::unique_ptr<CDBTTree> cdbttreeOut = std::make_unique<CDBTTree>(outputfile);
0077
0078   for (unsigned int i = 0; i < 96; i++)
0079   {
0080     for (unsigned int j = 0; j < 256; j++)
0081     {
0082       unsigned int key = TowerInfoDefs::encode_emcal(i, j);
0083       float val1 = cdbttree1->GetFloatValue(static_cast<int>(key), m_fieldname);
0084       float val2 = cdbttree2->GetFloatValue(static_cast<int>(key), m_fieldname);
0085       cdbttreeOut->SetFloatValue(static_cast<int>(key), m_fieldname, val1 * val2);
0086     }
0087   }
0088
0089   cdbttreeOut->Commit();
0090   cdbttreeOut->WriteCDBTTree();
0091
0092 }  // end macro
0093
0094 void TSCtoCDBTTree(const char *infile, const char *outputfile, const std::string &m_fieldname)
0095 {
0096   bool chk4file = gSystem->AccessPathName(infile);
0097   std::unique_ptr<TFile> f = nullptr;
0098
0099   if (!chk4file)
0100   {
0101     f = std::make_unique<TFile>(infile, "READ");
0102   }
0103   else
0104   {
0105     std::cout << "File " << infile << " cant be found in current directory." << std::endl;
0106     exit(0);
0107   }
0108
0109   // write to cdb tree
0110   std::unique_ptr<CDBTTree> cdbttree = std::make_unique<CDBTTree>(outputfile);
0111
0112   // gain values lie in the 2d histogram called corrPat
0113   TH2F *cp = static_cast<TH2F *>(f->Get("corrPat"));
0114
0115   for (int i = 0; i < 96; i++)
0116   {
0117     for (int j = 0; j < 256; j++)
0118     {
0119       unsigned int key = TowerInfoDefs::encode_emcal(static_cast<unsigned int>(i), static_cast<unsigned int>(j));
0120       float gain = static_cast<float>(1.0 / cp->GetBinContent(i + 1, j + 1));
0121       if (cp->GetBinContent(i + 1, j + 1) <= 0)
0122       {
0123         gain = 0;
0124       }
0125       if (std::isnan(cp->GetBinContent(i + 1, j + 1)))
0126       {
0127         gain = 0;
0128         std::cout << "nan calib from tsc " << i << "," << j << std::endl;
0129       }
0130       cdbttree->SetFloatValue(static_cast<int>(key), m_fieldname.c_str(), gain);
0131     }
0132   }
0133
0134   cdbttree->Commit();
0135   cdbttree->WriteCDBTTree();
0136   // cdbttree->Print();
0137   f->Close();
0138 }
0139
0140 #ifndef __CINT__
0141 int main(int argc, const char *const argv[])
0142 {
0143   const std::vector<std::string> args(argv, argv + argc);
0144
0145   if (args.size() < 3 || args.size() > 5)
0146   {
0147     std::cout << "usage: " << args[0] << " hist_fname calib_fname [iter] [m_fieldname]" << std::endl;
0148     return 1;
0149   }
0150
0151   const std::string& hist_fname = args[1];
0152   const std::string& calib_fname = args[2];
0153   int iter = 1;
0154   std::string m_fieldname = "CEMC_calib_ADC_to_ETower";
0155
0156   if (args.size() >= 4)
0157   {
0158     iter = std::stoi(args[3]);
0159   }
0160   if (args.size() >= 5)
0161   {
0162     m_fieldname = args[4];
0163   }
0164
0165   doTscFit(hist_fname, calib_fname, iter, m_fieldname);
0166
0167   std::cout << "======================================" << std::endl;
0168   std::cout << "done" << std::endl;
0169   return 0;
0170 }
0171 #endif