calibrations/mbd/cal_mbd.C

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