mbd/THRESHOLDS/mbdfem_thresholds.C

0001 #include <mbd/MbdCalib.h>
0002
0003 #include <TCanvas.h>
0004 #include <TChain.h>
0005 #include <TF1.h>
0006 #include <TFile.h>
0007 #include <TH2.h>
0008 #include <TMath.h>
0009 #include <TStyle.h>
0010
0011 #include <filesystem>
0012 #include <fstream>
0013
0014 R__LOAD_LIBRARY(libmbd.so)
0015
0016 TH2 *h2_qraw{nullptr};
0017 TH2 *h2_qcut{nullptr};
0018 TH2 *h2_qratio{nullptr};
0019 TH1 *h_qraw[128]{nullptr};
0020 TH1 *h_qcut[128]{nullptr};
0021 TH1 *h_qratio[128]{nullptr};
0022 TF1 *turnon[128] = {nullptr};
0023
0024 TCanvas *ac{nullptr};
0025
0026 Double_t TrigTurnOn(Double_t *x, Double_t *par)  // NOLINT(readability-non-const-parameter)
0027 {
0028   // par[0] is the amplitude (turn on max)
0029   // par[1] is the stretch (how long it takes to reach max)
0030   // par[2] is the offset (where it reaches max)
0031   // x[0] is the energy (or pt)
0032   Double_t y = 0.5 * par[0] * TMath::Erfc(-(sqrt(2) / par[1]) * (x[0] - par[2]));
0033
0034   return y;
0035 }
0036
0037 void ana_channels(const std::string &fname = "thresholds-00042123-0001.root")
0038 {
0039   TFile *infile{nullptr};
0040   if (!fname.empty())
0041   {
0042     infile = new TFile(fname.c_str(), "READ");
0043     h2_qraw = (TH2F *) infile->Get("h2_qraw");
0044     h2_qcut = (TH2F *) infile->Get("h2_qcut");
0045     h2_qratio = (TH2F *) infile->Get("h2_qratio");
0046   }
0047
0048   // outfile is used for updating thresholds in 1008
0049   TString name = fname;
0050   name.ReplaceAll(".root", ".out");
0051   std::ofstream outfile(name.Data());
0052
0053   name = fname;
0054   name.ReplaceAll(".root", ".calib");
0055   std::ofstream calibfile(name.Data());
0056
0057   ac = new TCanvas("ac", "threshold turn-on", 800, 600);
0058
0059   TString pdfname = fname;
0060   pdfname.ReplaceAll(".root", ".pdf");
0061   ac->Print(pdfname + "[");
0062
0063   TString title;
0064   for (int ipmt = 0; ipmt < 128; ipmt++)
0065   {
0066     name = "h_qraw";
0067     name += ipmt;
0068     h_qraw[ipmt] = h2_qraw->ProjectionX(name, ipmt + 1, ipmt + 1);
0069     name = "h_qcut";
0070     name += ipmt;
0071     h_qcut[ipmt] = h2_qcut->ProjectionX(name, ipmt + 1, ipmt + 1);
0072     name = "h_ratio";
0073     name += ipmt;
0074     h_qratio[ipmt] = h2_qratio->ProjectionX(name, ipmt + 1, ipmt + 1);
0075     title = "ch";
0076     title += ipmt;
0077     title += "threshold";
0078     h_qratio[ipmt]->SetTitle(title);
0079
0080     h_qraw[ipmt]->Sumw2();
0081     h_qcut[ipmt]->Sumw2();
0082
0083     /*
0084     name = "h_ratio"; name += ipmt;
0085     h_qratio[ipmt] = (TH1*)h_qcut[ipmt]->Clone( name );
0086     h_qratio[ipmt]->Divide( h_qcut[ipmt], h_qraw[ipmt], 1, 1, "B" );
0087     */
0088
0089     // force eff=0 at adc=0
0090     h_qratio[ipmt]->SetBinContent(1, 0.);
0091     h_qratio[ipmt]->SetBinError(1, 0.001);
0092
0093     name = "turnon";
0094     name += ipmt;
0095     turnon[ipmt] = new TF1(name, TrigTurnOn, 0, 200, 3);
0096     turnon[ipmt]->SetLineColor(4);
0097
0098     turnon[ipmt]->SetParameters(1, 20, 40);
0099     h_qratio[ipmt]->Fit(turnon[ipmt], "R");
0100     h_qratio[ipmt]->Draw();
0101     gPad->Modified();
0102     gPad->Update();
0103
0104     ac->Print(pdfname);
0105
0106     double mean = turnon[ipmt]->GetParameter(2);   // midpoint of turn-on
0107     double width = turnon[ipmt]->GetParameter(1);  // width of turn-on
0108     double eff = turnon[ipmt]->GetParameter(0);    // efficiency at max
0109
0110     // For changing thresholds
0111     double threshold_to_set = 5;                                       // if we want threshold at ADC=50
0112     double effic_at = threshold_to_set - round((mean + width) / 10.);  // about 10 mV per ADC
0113
0114     double meanerr = turnon[ipmt]->GetParError(2);   // error, midpoint of turn-on
0115     double widtherr = turnon[ipmt]->GetParError(1);  // error, width of turn-on
0116     double efferr = turnon[ipmt]->GetParError(0);    // error, efficiency
0117
0118     double chi2ndf = turnon[ipmt]->GetChisquare() / turnon[ipmt]->GetNDF();  // error, width of turn-on
0119
0120     outfile << ipmt << "\t" << mean << "\t" << width << "\t" << effic_at << std::endl;
0121     calibfile << ipmt << "\t" << mean << "\t" << meanerr << "\t" << width << "\t" << widtherr << "\t"
0122               << eff << "\t" << efferr << chi2ndf << std::endl;
0123
0124     /*
0125     if ( ipmt==23 )
0126     {
0127       string junk;
0128       std::cout << "? ";
0129       cin >> junk;
0130     }
0131     */
0132   }
0133   ac->Print(pdfname + "]");
0134
0135   outfile.close();
0136   calibfile.close();
0137 }
0138
0139 void mbdfem_thresholds(const std::filesystem::path &fname = "DST_UNCALMBD-00042123-0001.root")
0140 {
0141   gStyle->SetOptStat(0);
0142
0143   TChain *T = new TChain("T");
0144   TString savefname = "thresholds_";
0145   savefname += fname.stem().c_str();
0146
0147   if (fname.extension() == ".list")
0148   {
0149     std::cout << "is list" << std::endl;
0150     std::string dstfname;
0151
0152     std::ifstream inflist(fname.c_str());
0153     while (inflist >> dstfname)
0154     {
0155       std::cout << "Added " << dstfname << std::endl;
0156       T->Add(dstfname.c_str());
0157     }
0158
0159     savefname += "-99999.root";
0160   }
0161   else
0162   {
0163     std::cout << "is root file" << std::endl;
0164     T->Add(fname.c_str());
0165
0166     savefname.ReplaceAll("DST_UNCALMBD_", "");
0167     savefname.ReplaceAll("DST_CALO_", "");
0168     savefname.ReplaceAll("DST_CALOFITTING_", "");
0169     savefname += ".root";
0170   }
0171
0172   TFile *savefile = new TFile(savefname, "RECREATE");
0173
0174   h2_qraw = new TH2F("h2_qraw", "ch vs bq", 200, 0.5, 200.5, 128, 0, 128);
0175   h2_qraw->SetXTitle("ADC");
0176   h2_qraw->SetYTitle("pmt");
0177   h2_qcut = (TH2 *) h2_qraw->Clone("h2_qcut");
0178   h2_qcut->SetTitle("ch vs bq, cut");
0179   h2_qratio = (TH2 *) h2_qraw->Clone("h2_qratio");
0180   h2_qratio->SetTitle("ch vs bq, cut/raw ratio");
0181
0182   //  TCanvas *ac = new TCanvas("ac","ac",550*1.5,425*1.5);
0183
0184   T->Draw("bpmt:bq>>h2_qraw", "bq>0&&bq<200.5", "colz");
0185   T->Draw("bpmt:bq>>h2_qcut", "bq>0&&bq<200.5&&abs(btt)<50", "colz");
0186
0187   h2_qratio->Divide(h2_qcut, h2_qraw, 1, 1, "B");
0188   h2_qratio->Draw("colz");
0189
0190   // TCanvas *bc = new TCanvas("bc","bc",550*1.5,425*1.5);
0191   // ana_channels(savefname.Data());
0192
0193   savefile->cd();
0194
0195   h2_qraw->Write();
0196   h2_qcut->Write();
0197   h2_qratio->Write();
0198
0199   /*
0200      for (int ich=0; ich<128; ich++)
0201      {
0202      h_qraw[ich]->Write();
0203      h_qcut[ich]->Write();
0204      h_qcut[ich]->SetLineColor(2);
0205      h_qcut[ich]->SetMarkerColor(2);
0206      h_qratio[ich]->Write();
0207      }
0208      */
0209
0210   savefile->Write();
0211 }