mbd/THRESHOLDS/mbdfem_thresholds.C

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