File indexing completed on 2025-08-03 08:20:22
0001 #include <iostream>
0002 #include <fstream>
0003
0004 #include <TString.h>
0005 #include <TFile.h>
0006 #include <TTree.h>
0007 #include <TH1F.h>
0008 #include <TH2F.h>
0009 #include <TF1.h>
0010 #include <TCanvas.h>
0011 #include <TPad.h>
0012 #include <TSystem.h>
0013
0014 #include <ffamodules/CDBInterface.h>
0015 #include <phool/recoConsts.h>
0016
0017 #include "get_runstr.h"
0018 #include "make_cdbtree.C"
0019 #include <mbd/MbdCalib.h>
0020 #include <mbd/MbdDefs.h>
0021 #include "read_dstmbd.C"
0022 #include "recal_mbd_mip.C"
0023
0024 #if defined(__CLING__)
0025 R__LOAD_LIBRARY(libffamodules.so)
0026 R__LOAD_LIBRARY(libphool.so)
0027 R__LOAD_LIBRARY(libmbd.so)
0028 #endif
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038 void cal_mbd(const char *tfname = "DST_MBDUNCAL-00020869-0000.root", const int subpass = 0, const int nevt = 0, const int runtype = 0, const std::string_view dbtag = "")
0039 {
0040 cout << "cal_mbd(), tfname " << tfname << endl;
0041 cout << "cal_mbd(), runtype " << runtype << endl;
0042 read_dstmbd( tfname );
0043
0044
0045 int runnumber = get_runnumber(tfname);
0046 TString dir = "results/";
0047 dir += runnumber;
0048 dir += "/";
0049 TString name = "mkdir -p " + dir;
0050 gSystem->Exec( name );
0051
0052 cout << name << endl;
0053
0054
0055
0056
0057 MbdCalib *mcal = new MbdCalib();
0058 mcal->Verbosity(1);
0059
0060
0061 TString calfile = dir + "/mbd_sampmax.calib";
0062 if ( gSystem->AccessPathName(calfile)==0 )
0063 {
0064 mcal->Download_SampMax( calfile.Data() );
0065 }
0066 calfile = dir + "/mbd_ped.calib";
0067 if ( gSystem->AccessPathName(calfile)==0 )
0068 {
0069 mcal->Download_Ped( calfile.Data() );
0070 }
0071
0072 TString savefname = dir;
0073 if ( subpass==0 )
0074 {
0075 savefname += "calmbdpass2."; savefname += subpass; savefname += "_time-"; savefname += runnumber; savefname += ".root";
0076
0077
0078
0079
0080
0081
0082
0083 }
0084 else if ( subpass==1 || subpass==2 )
0085 {
0086 savefname += "calmbdpass2."; savefname += subpass; savefname += "_slew-"; savefname += runnumber; savefname += ".root";
0087 }
0088 else if ( subpass==3 )
0089 {
0090 savefname += "calmbdpass2."; savefname += subpass; savefname += "_q-"; savefname += runnumber; savefname += ".root";
0091 }
0092 cout << "saving to " << savefname << endl;
0093
0094
0095 if ( subpass==3 )
0096 {
0097 calfile = dir + "/pass0_mbd_tq_t0.calib";
0098 mcal->Download_TQT0( calfile.Data() );
0099 cout << "Loaded " << calfile << endl;
0100
0101 calfile = dir + "/pass0_mbd_tt_t0.calib";
0102 mcal->Download_TTT0( calfile.Data() );
0103 cout << "Loaded " << calfile << endl;
0104 }
0105 if ( subpass>1 )
0106 {
0107 if ( dbtag.empty() )
0108 {
0109 calfile = dir + "/mbd_slewcorr.calib";
0110 mcal->Download_SlewCorr( calfile.Data() );
0111 cout << "Loaded " << calfile << endl;
0112 }
0113 else
0114 {
0115 recoConsts *rc = recoConsts::instance();
0116 rc->set_StringFlag("CDB_GLOBALTAG","ProdA_2024");
0117 rc->set_uint64Flag("TIMESTAMP",runnumber);
0118 CDBInterface *cdb = CDBInterface::instance();
0119 std::string slew_url = cdb->getUrl("MBD_SLEWCORR");
0120 mcal->Download_SlewCorr(slew_url);
0121 cout << "Loaded " << slew_url << endl;
0122 }
0123 }
0124 if ( subpass>3 )
0125 {
0126 calfile = dir + "/mbd_qfit.calib";
0127 mcal->Download_Gains( calfile.Data() );
0128 cout << "Loaded " << calfile << endl;
0129 }
0130
0131 TFile *savefile = new TFile(savefname,"RECREATE");
0132
0133 TH1 *h_q[MbdDefs::MBD_N_PMT];
0134 TH1 *h_tt[MbdDefs::MBD_N_PMT];
0135 TH1 *h_tq[MbdDefs::MBD_N_PMT];
0136
0137 TH2 *h2_slew[MbdDefs::MBD_N_PMT];
0138
0139 TString title;
0140 for (int ipmt=0; ipmt<MbdDefs::MBD_N_PMT; ipmt++)
0141 {
0142 name = "h_q"; name += ipmt;
0143 title = "q"; title += ipmt;
0144 h_q[ipmt] = new TH1F(name,title,3000,-100,15000-100);
0145 h_q[ipmt]->SetXTitle("ADC");
0146
0147 name = "h_tt"; name += ipmt;
0148 title = "tt"; title += ipmt;
0149
0150 h_tt[ipmt] = new TH1F(name,title,7000,-30.,30.);
0151 h_tt[ipmt]->SetXTitle("ns");
0152
0153 name = "h_tq"; name += ipmt;
0154 title = "tq"; title += ipmt;
0155 h_tq[ipmt] = new TH1F(name,title,7000,-150,31*17.76);
0156 h_tq[ipmt]->SetXTitle("ns");
0157
0158 if ( subpass>0 )
0159 {
0160 name = "h2_slew"; name += ipmt;
0161 title = "slew curve, ch "; title += ipmt;
0162 h2_slew[ipmt] = new TH2F(name,title,4000,-0.5,16000-0.5,1100,-5,6);
0163 h2_slew[ipmt]->SetXTitle("ADC");
0164 h2_slew[ipmt]->SetYTitle("#Delta T (ns)");
0165 }
0166 }
0167 TH2 *h2_tq = new TH2F("h2_tq","ch vs tq",900,-150,150,MbdDefs::MBD_N_PMT,-0.5,MbdDefs::MBD_N_PMT-0.5);
0168 h2_tq->SetXTitle("tq [ns]");
0169 h2_tq->SetYTitle("pmt ch");
0170 TH2 *h2_tt = new TH2F("h2_tt","ch vs tt",900,-150,150,MbdDefs::MBD_N_PMT,-0.5,MbdDefs::MBD_N_PMT-0.5);
0171 h2_tt->SetXTitle("tt [ns]");
0172 h2_tt->SetYTitle("pmt ch");
0173
0174
0175 int nentries = tree->GetEntries();
0176 if ( nevt!=0 && nevt<nentries )
0177 {
0178 nentries = nevt;
0179 }
0180
0181 double armtime[2]{0.,0.};
0182 double nhit[2]{0.,0.};
0183 float ttcorr[128] = {0.};
0184
0185 std::cout << "Processing " << nentries << std::endl;
0186 for (int ientry=0; ientry<nentries; ientry++)
0187 {
0188 dstmbd_GetEntry(ientry);
0189
0190 if (ientry<4)
0191 {
0192
0193 cout << f_evt << "\tch0\t" << f_tt[0] << "\t" << f_tq[0] << "\t" << f_q[0] << endl;
0194 cout << "ch127\t" << f_tt[127] << "\t" << f_tq[127] << "\t" << f_q[127] << endl;
0195 }
0196
0197
0198
0199
0200
0201
0202
0203
0204
0205
0206
0207 if ( subpass!=0 && fabs(f_bz)>60. )
0208 {
0209 continue;
0210 }
0211
0212 if ( subpass>0 )
0213 {
0214 armtime[0] = 0.;
0215 armtime[1] = 0.;
0216 nhit[0] = 0.;
0217 nhit[1] = 0.;
0218 }
0219
0220 for (int ipmt=0; ipmt<MbdDefs::MBD_N_PMT; ipmt++)
0221 {
0222 int arm = ipmt/64;
0223
0224 ttcorr[ipmt] = f_tt[ipmt];
0225 float tq = f_tq[ipmt];
0226
0227 if ( subpass>0 )
0228 {
0229 if ( !isnan(mcal->get_tt0(ipmt)) && mcal->get_tt0(ipmt)>-100. && f_q[ipmt]>0. && f_q[ipmt]<16000. )
0230 {
0231 ttcorr[ipmt] -= mcal->get_tt0(ipmt);
0232 }
0233 if ( !isnan(mcal->get_tq0(ipmt)) && mcal->get_tq0(ipmt)>-100. )
0234 {
0235 tq -= mcal->get_tq0(ipmt);
0236 }
0237 }
0238
0239
0240 if ( subpass>1 && subpass<3 )
0241 {
0242 int feech = (ipmt / 8) * 16 + ipmt % 8;
0243 if ( !isnan(mcal->get_tt0(ipmt)) && mcal->get_tt0(ipmt)>-100. && f_q[ipmt]>0. && f_q[ipmt]<16000. )
0244 {
0245 ttcorr[ipmt] -= mcal->get_scorr(feech,f_q[ipmt]);
0246 }
0247 }
0248
0249 h_tt[ipmt]->Fill( ttcorr[ipmt] );
0250 h2_tt->Fill( ttcorr[ipmt], ipmt );
0251
0252 h_tq[ipmt]->Fill( tq );
0253 h2_tq->Fill( tq, ipmt );
0254
0255 if ( subpass==0 && fabs(f_tt[ipmt])<26. )
0256 {
0257 h_q[ipmt]->Fill( f_q[ipmt] );
0258 }
0259
0260 else if ( subpass>0 && (fabs(ttcorr[ipmt])<26.||f_q[ipmt]>40.) )
0261 {
0262 h_q[ipmt]->Fill( f_q[ipmt] );
0263
0264
0265
0266 if ( fabs(ttcorr[ipmt])<26. )
0267 {
0268 nhit[arm] += 1.0;
0269 armtime[arm] += ttcorr[ipmt];
0270 }
0271
0272
0273
0274
0275
0276
0277
0278
0279
0280
0281
0282
0283
0284 }
0285 }
0286
0287
0288 for (int iarm=0; iarm<2; iarm++)
0289 {
0290 if ( nhit[iarm]>1. )
0291 {
0292 armtime[iarm] = armtime[iarm]/nhit[iarm];
0293 }
0294
0295 }
0296
0297 for (int ipmt=0; ipmt<MbdDefs::MBD_N_PMT; ipmt++)
0298 {
0299
0300 int arm = ipmt/64;
0301
0302 if ( nhit[arm]<2 || f_q[ipmt]<=0. ) continue;
0303
0304 double dt = ttcorr[ipmt] - armtime[arm];
0305
0306
0307 h2_slew[ipmt]->Fill( f_q[ipmt], dt );
0308 }
0309 }
0310
0311 TCanvas *ac[100];
0312 int cvindex = 0;
0313
0314 TString pdfname;
0315
0316
0317 ac[cvindex] = new TCanvas("cal_tt","ch vs tt",425*1.5,550*1.5);
0318 h2_tt->Draw("colz");
0319 pdfname = dir; pdfname += "calmbdpass2."; pdfname += subpass; pdfname += "_time-"; pdfname += runnumber; pdfname += ".pdf";
0320 if ( subpass==0 )
0321 {
0322
0323 title = "t0, Timing Ch's";
0324 }
0325 else if ( subpass>0 )
0326 {
0327
0328 title = "t0 corrected, Timing Ch's";
0329 h2_tt->GetXaxis()->SetRangeUser( -12.5, 12.5 );
0330 gPad->Modified();
0331 gPad->Update();
0332 }
0333 cout << pdfname << endl;
0334 ac[cvindex]->Print( pdfname + "[");
0335 ac[cvindex]->Print( pdfname, title );
0336 ++cvindex;
0337
0338
0339 ac[cvindex] = new TCanvas("cal_tq","ch vs tq",425*1.5,550*1.5);
0340 h2_tq->Draw("colz");
0341 if ( subpass==0 )
0342 {
0343
0344 title = "t0, Charge Ch's";
0345 }
0346 else if ( subpass>0 )
0347 {
0348
0349 title = "t0 corrected, Charge Ch's";
0350 h2_tq->GetXaxis()->SetRangeUser( -20,20 );
0351 gPad->Modified();
0352 gPad->Update();
0353 }
0354 cout << pdfname << endl;
0355 ac[cvindex]->Print( pdfname, title );
0356 ++cvindex;
0357
0358
0359 ac[cvindex] = new TCanvas("cal_tt_ch","tt",550*1.5,425*1.5);
0360 gPad->SetLogy(1);
0361
0362 ofstream cal_tt_t0_file;
0363 TString cal_fname = dir; cal_fname += "pass"; cal_fname += subpass; cal_fname += "_mbd_tt_t0.calib";
0364 cal_tt_t0_file.open( cal_fname );
0365 cout << "Creating " << cal_fname << endl;
0366
0367 TF1 *gaussian = new TF1("gaussian","gaus",-25,25);
0368 gaussian->SetLineColor(2);
0369 double min_twindow = -25.;
0370 double max_twindow = 25.;
0371
0372 for (int ipmt=0; ipmt<MbdDefs::MBD_N_PMT; ipmt++)
0373 {
0374
0375 if ( ipmt==0 || ipmt==64 )
0376 {
0377
0378 h_tt[ipmt]->SetAxisRange(-25.,25.);
0379
0380 }
0381 else
0382 {
0383
0384
0385 h_tt[ipmt]->SetAxisRange(min_twindow,max_twindow);
0386 }
0387 Double_t peak = h_tt[ipmt]->GetMaximum();
0388 int peakbin = h_tt[ipmt]->GetMaximumBin();
0389 Double_t mean = h_tt[ipmt]->GetBinCenter( peakbin );
0390 Double_t sigma = 1.0;
0391
0392 gaussian->SetParameters( peak, mean, 5 );
0393 gaussian->SetRange( mean-3*sigma, mean+3*sigma );
0394
0395 if ( ipmt==0 || ipmt==64 )
0396 {
0397 min_twindow = mean - 3*sigma;
0398 max_twindow = mean + 3*sigma;
0399 }
0400
0401 h_tt[ipmt]->Fit(gaussian,"R");
0402 h_tt[ipmt]->SetAxisRange(mean-15.*sigma,mean+15.*sigma);
0403
0404 mean = gaussian->GetParameter(1);
0405 Double_t meanerr = gaussian->GetParError(1);
0406 sigma = gaussian->GetParameter(2);
0407 Double_t sigmaerr = gaussian->GetParError(2);
0408
0409
0410 int nbinsx = h2_tt->GetNbinsX();
0411 for (int ibinx=1; ibinx<=nbinsx; ibinx++)
0412 {
0413 Double_t fitpeak = gaussian->GetParameter(0);
0414 if ( fitpeak!=0 )
0415 {
0416 Float_t bincontent = h2_tt->GetBinContent(ibinx,ipmt+1);
0417 h2_tt->SetBinContent(ibinx,ipmt+1,bincontent/fitpeak);
0418 }
0419 }
0420
0421 cal_tt_t0_file << ipmt << "\t" << mean << "\t" << meanerr << "\t" << sigma << "\t" << sigmaerr << endl;
0422 cout << ipmt << "\t" << mean << "\t" << meanerr << "\t" << sigma << "\t" << sigmaerr << endl;
0423
0424 if ( subpass==0 )
0425 {
0426
0427 title = "h_tt"; title += ipmt;
0428 }
0429 else if ( subpass==1 || subpass==2 )
0430 {
0431
0432 title = "h_ttcorr"; title += ipmt;
0433 }
0434
0435 ac[cvindex]->Print( pdfname, title );
0436 }
0437 cal_tt_t0_file.close();
0438 make_cdbtree( cal_fname );
0439 ++cvindex;
0440
0441
0442 ac[cvindex] = new TCanvas("cal_tq_ch","tq",550*1.5,425*1.5);
0443 gPad->SetLogy(1);
0444
0445 ofstream cal_tq_t0_file;
0446 cal_fname = dir; cal_fname += "pass"; cal_fname += subpass; cal_fname += "_mbd_tq_t0.calib";
0447 cal_tq_t0_file.open( cal_fname );
0448
0449 for (int ipmt=0; ipmt<MbdDefs::MBD_N_PMT; ipmt++)
0450 {
0451
0452 if ( ipmt==0 || ipmt==64 )
0453 {
0454
0455 h_tq[ipmt]->SetAxisRange(-25.,25.);
0456
0457 }
0458 else
0459 {
0460
0461
0462 h_tq[ipmt]->SetAxisRange(min_twindow,max_twindow);
0463 }
0464 Double_t peak = h_tq[ipmt]->GetMaximum();
0465 int peakbin = h_tq[ipmt]->GetMaximumBin();
0466 Double_t mean = h_tq[ipmt]->GetBinCenter( peakbin );
0467 Double_t sigma = 1.0;
0468
0469 gaussian->SetParameters( peak, mean, 5 );
0470 gaussian->SetRange( mean-3*sigma, mean+3*sigma );
0471
0472 if ( ipmt==0 || ipmt==64 )
0473 {
0474 min_twindow = mean - 3*sigma;
0475 max_twindow = mean + 3*sigma;
0476 }
0477
0478 h_tq[ipmt]->Fit(gaussian,"R");
0479 h_tq[ipmt]->SetAxisRange(mean-12.*sigma,mean+12.*sigma);
0480
0481 mean = gaussian->GetParameter(1);
0482 Double_t meanerr = gaussian->GetParError(1);
0483 sigma = gaussian->GetParameter(2);
0484 Double_t sigmaerr = gaussian->GetParError(2);
0485
0486
0487 int nbinsx = h2_tq->GetNbinsX();
0488 for (int ibinx=1; ibinx<=nbinsx; ibinx++)
0489 {
0490 Double_t fitpeak = gaussian->GetParameter(0);
0491 if ( fitpeak!=0 )
0492 {
0493 Float_t bincontent = h2_tq->GetBinContent(ibinx,ipmt+1);
0494 h2_tq->SetBinContent(ibinx,ipmt+1,bincontent/fitpeak);
0495 }
0496 }
0497
0498 cal_tq_t0_file << ipmt << "\t" << mean << "\t" << meanerr << "\t" << sigma << "\t" << sigmaerr << endl;
0499
0500 if ( subpass==0 )
0501 {
0502
0503 title = "h_tq"; title += ipmt;
0504 }
0505 else if ( subpass==1 || subpass==2 )
0506 {
0507
0508 title = "h_tqcorr"; title += ipmt;
0509 }
0510
0511 ac[cvindex]->Print( pdfname, title );
0512 }
0513
0514 cal_tq_t0_file.close();
0515 make_cdbtree( cal_fname );
0516
0517 ac[0]->Print( pdfname + "]" );
0518 ++cvindex;
0519
0520 if ( subpass==1 || subpass==2 )
0521 {
0522
0523 ac[cvindex] = new TCanvas("cal_slew","slew",425*1.5,550*1.5);
0524 pdfname = dir; pdfname += "calmbdpass2."; pdfname += subpass; pdfname += "_slew-"; pdfname += runnumber; pdfname += ".pdf";
0525 ac[cvindex]->Print( pdfname + "[" );
0526
0527 for (int ipmt=0; ipmt<MbdDefs::MBD_N_PMT; ipmt++)
0528 {
0529 h2_slew[ipmt]->Draw("colz");
0530 gPad->SetLogz(1);
0531
0532
0533 name = dir + "h2_slew"; name += ipmt; name += "_pass"; name += subpass;
0534 cout << name << endl;
0535 ac[cvindex]->Print( pdfname, name );
0536 }
0537
0538
0539
0540 ac[cvindex]->Print( pdfname + "]" );
0541 ++cvindex;
0542 }
0543
0544
0545 if ( subpass==0 )
0546 {
0547 ac[cvindex] = new TCanvas("cal_q","q",425*1.5,550*1.5);
0548
0549 pdfname = dir; pdfname += "calmbdpass2."; pdfname += subpass; pdfname += "_adc-"; pdfname += runnumber; pdfname += ".pdf";
0550 ac[cvindex]->Print( pdfname + "[" );
0551
0552 gPad->SetLogy(1);
0553
0554 for (int ipmt=0; ipmt<MbdDefs::MBD_N_PMT && subpass==0; ipmt++)
0555 {
0556 h_q[ipmt]->Draw();
0557
0558
0559 title = "h_adc"; title += ipmt;
0560
0561 ac[cvindex]->Print( pdfname, title );
0562 }
0563 ac[cvindex]->Print( pdfname + "]" );
0564 ++cvindex;
0565 }
0566
0567 savefile->Write();
0568
0569
0570 if ( subpass==3 )
0571 {
0572 recal_mbd_mip( tfname, subpass, runtype );
0573 }
0574
0575 }
0576