calorimeter/calo_cdb/myUtils.cc

0001 #include "myUtils.h"
0002 #include "geometry_constants.h"
0003
0004 // root includes --
0005 #include <TF1.h>
0006 #include <TFitResult.h>
0007
0008 // c++ includes --
0009 #include <memory>
0010
0011 TFitResultPtr myUtils::doGausFit(TH1 *hist, Double_t start, Double_t end, const std::string &name)
0012 {
0013   // fit calib hist
0014   std::unique_ptr<TF1> fitFunc = std::make_unique<TF1>(name.c_str(), "gaus", start, end);
0015
0016   Double_t initialAmplitude = hist->GetMaximum();
0017   Double_t initialMean = hist->GetMean();
0018   Double_t initialSigma = hist->GetRMS();
0019
0020   fitFunc->SetParameter(0, initialAmplitude);
0021   fitFunc->SetParameter(1, initialMean);
0022   fitFunc->SetParameter(2, initialSigma);
0023
0024   // You can also set parameter names for better readability in the stats box
0025   fitFunc->SetParName(0, "Amplitude");
0026   fitFunc->SetParName(1, "Mean");
0027   fitFunc->SetParName(2, "Sigma");
0028
0029   // Set some visual properties for the fit line
0030   fitFunc->SetLineColor(kRed);
0031   fitFunc->SetLineWidth(2);
0032   fitFunc->SetLineStyle(kDashed);  // Optional: make it dashed
0033
0034   TFitResultPtr fitResult = hist->Fit(fitFunc.get(), "RS");  // Fit within range, store result, quiet
0035
0036   if (fitResult.Get())
0037   {  // Check if TFitResultPtr is valid
0038     std::cout << "\n----------------------------------------------------" << std::endl;
0039     std::cout << "Fit Results for function: " << fitFunc->GetName() << std::endl;
0040     std::cout << "----------------------------------------------------" << std::endl;
0041     std::cout << "Fit Status: " << fitResult->Status() << " (0 means successful)" << std::endl;
0042     if (fitResult->IsValid())
0043     {  // Check if the fit is valid (e.g., covariance matrix is good)
0044       std::cout << "Fit is Valid." << std::endl;
0045       for (UInt_t i = 0; i < static_cast<UInt_t>(fitFunc->GetNpar()); ++i)
0046       {
0047         std::cout << "Parameter " << fitFunc->GetParName(static_cast<Int_t>(i)) << " (" << i << "): "
0048                   << fitResult->Parameter(i) << " +/- " << fitResult->ParError(i) << std::endl;
0049       }
0050       std::cout << "Resolution: Sigma/Mean = " << fitResult->Parameter(2) / fitResult->Parameter(1) << std::endl;
0051       std::cout << "Chi^2 / NDF: " << fitResult->Chi2() << " / " << fitResult->Ndf()
0052                 << " = " << (fitResult->Ndf() > 0 ? fitResult->Chi2() / fitResult->Ndf() : 0) << std::endl;
0053       std::cout << "Probability: " << TMath::Prob(fitResult->Chi2(), static_cast<Int_t>(fitResult->Ndf())) << std::endl;
0054     }
0055     else
0056     {
0057       std::cout << "Fit is NOT Valid." << std::endl;
0058     }
0059     std::cout << "----------------------------------------------------" << std::endl;
0060   }
0061   else
0062   {
0063     std::cout << "Fit did not return a valid TFitResultPtr." << std::endl;
0064   }
0065
0066   return fitResult;
0067 }
0068
0069 std::vector<std::string> myUtils::split(const std::string &s, const char delimiter)
0070 {
0071   std::vector<std::string> result;
0072
0073   std::stringstream ss(s);
0074   std::string temp;
0075
0076   while (getline(ss, temp, delimiter))
0077   {
0078     if (!temp.empty())
0079     {
0080       result.push_back(temp);
0081     }
0082   }
0083
0084   return result;
0085 }
0086
0087 void myUtils::setEMCalDim(TH1 *hist)
0088 {
0089   hist->GetXaxis()->SetLimits(0, 256);
0090   hist->GetXaxis()->SetNdivisions(32, false);
0091   hist->GetXaxis()->SetLabelSize(0.04F);
0092   hist->GetXaxis()->SetTickSize(0.01F);
0093   hist->GetYaxis()->SetTickSize(0.01F);
0094   hist->GetYaxis()->SetLabelSize(0.04F);
0095   hist->GetYaxis()->SetLimits(0, 96);
0096   hist->GetYaxis()->SetNdivisions(12, false);
0097   hist->GetYaxis()->SetTitleOffset(0.5);
0098   hist->GetXaxis()->SetTitleOffset(1);
0099 }
0100
0101 std::pair<int, int> myUtils::getSectorIB(int iphi, int ieta)
0102 {
0103   int k = iphi / CaloGeometry::CEMC_NTOW_IB_SIDE;
0104
0105   int sector = (ieta < 48) ? k + 32 : k;
0106   int ib = (ieta < 48) ? CaloGeometry::CEMC_NIB_PER_SECTOR - (ieta / CaloGeometry::CEMC_NTOW_IB_SIDE) - 1 : (ieta - 48) / CaloGeometry::CEMC_NTOW_IB_SIDE;
0107
0108   return std::make_pair(sector, ib);
0109 }
0110
0111 std::pair<int, int> myUtils::getSectorIB(int towerIndex)
0112 {
0113   int k = towerIndex / CaloGeometry::CEMC_NCHANNEL_PER_SECTOR;
0114
0115   int sector = (k % 2) ? (k - 1) / 2 : (k / 2) + 32;
0116   int ib = (k % CaloGeometry::CEMC_NCHANNEL_PER_SECTOR) / CaloGeometry::CEMC_NCHANNEL_PER_IB;
0117
0118   return std::make_pair(sector, ib);
0119 }