Back to home page

sPhenix code displayed by LXR

 
 

    


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 // This macro executes the sub-passes for the pass2 calibrations
0031 //pass 2.0: do tt_t0 and tq_t0 offset calibration
0032 //pass 2.1: do the slew correction calibration
0033 //pass 2.2: do the next iteration of the slew correction calibration
0034 //pass 2.3: do the mip fits for charge calibration 
0035 
0036 //runtype 0: au+au 200 GeV
0037 //runtype 1: p+p 200 GeV
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   //== Create output directory (should already exist)
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   //== Load in calib constants
0055   //Float_t tq_t0_offsets[MbdDefs::MBD_N_PMT] = {};
0056 
0057   MbdCalib *mcal = new MbdCalib();
0058   mcal->Verbosity(1);
0059 
0060   // Load local mbd_sampmax and ped if they exist
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     // temp, for final t0
0078     /*
0079     calfile = dir + "/mbd_slewcorr.calib";
0080     mcal->Download_SlewCorr( calfile.Data() );
0081     cout << "Loaded " << calfile << endl;
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   // Load whatever calibrations are available at each subpass
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     //h_tt[ipmt] = new TH1F(name,title,7000,-150,31*17.76);
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   // Event loop, each ientry is one triggered event
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       // print charge from channels 0 and 127
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     // make npmt cut
0198     /*
0199     if ( f_npmt==0 )
0200     {
0201       cout << "f_npmt == 0" << endl;
0202       continue;
0203     }
0204     */
0205 
0206     // Make vertex cut
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       //if ( subpass==0 )  // temp, for final t0
0240       if ( subpass>1 && subpass<3 ) // apply slewcorr if subpass > 1
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       //else if ( subpass>0 && fabs(ttcorr[ipmt])<26. )
0260       else if ( subpass>0 && (fabs(ttcorr[ipmt])<26.||f_q[ipmt]>40.) )  // to get around high threshold
0261       {
0262         h_q[ipmt]->Fill( f_q[ipmt] );
0263 
0264         //if ( f_q[ipmt] > 1000. && fabs(ttcorr[ipmt])<26. )    // for p+p
0265         //if ( f_q[ipmt] > 2000. && fabs(ttcorr[ipmt])<26. )
0266         if ( fabs(ttcorr[ipmt])<26. )
0267         {
0268           nhit[arm] += 1.0;
0269           armtime[arm] += ttcorr[ipmt];
0270         }
0271 
0272         // check charge
0273         // why cuts on f_bn?
0274         /*
0275         if ( arm==0 && f_bn[0]<30 )
0276         {
0277           h_q[ipmt]->Fill( f_q[ipmt] );
0278         }
0279         else if ( arm==1 && f_bn[1]<30 )
0280         {
0281           h_q[ipmt]->Fill( f_q[ipmt] );
0282         }
0283         */
0284       }
0285     } // end loop over PMTs
0286 
0287     // calc meantime for high amplitude
0288     for (int iarm=0; iarm<2; iarm++)
0289     {
0290       if ( nhit[iarm]>1. )
0291       {
0292         armtime[iarm] = armtime[iarm]/nhit[iarm];
0293       }
0294       //cout << "aaa " << iarm << "\t" << nhit[iarm] << "\t" << armtime[iarm] << endl;
0295     }
0296 
0297     for (int ipmt=0; ipmt<MbdDefs::MBD_N_PMT; ipmt++)
0298     {
0299       //int ifeech = (ipmt/8)*16 + 8 + ipmt%8;  // time ifeech only
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       //cout << "filling" << endl;
0307       h2_slew[ipmt]->Fill( f_q[ipmt], dt );
0308     }
0309   }
0310 
0311   TCanvas *ac[100];
0312   int cvindex = 0;
0313 
0314   TString pdfname;  // output pdf
0315 
0316   //== Calculate tq_t0 and tt_t0 on subpass, or apply tt_t0 if subpass>0
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     //name = dir + "h2_tt.png";
0323     title = "t0, Timing Ch's";
0324   }
0325   else if ( subpass>0 )
0326   {
0327     //name = dir + "h2_ttcorr.png";
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   // now tq_t0
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     //name = dir + "h2_tq.png";
0344     title = "t0, Charge Ch's";
0345   }
0346   else if ( subpass>0 )
0347   {
0348     //name = dir + "h2_tqcorr.png";
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   // Here we calculate tt_t0 and tq_t0, starting with tt_t0 first
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     // only look in the middle
0375     if ( ipmt==0 || ipmt==64 )
0376     {
0377       // use wide range for 1st channel in each arm
0378       h_tt[ipmt]->SetAxisRange(-25.,25.);
0379       //h_tt[ipmt]->SetAxisRange(9.,13.);
0380     }
0381     else
0382     {
0383       // for subsequent channels in an arm, use a window
0384       // determined by the 1st channel
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     //Double_t sigma = 3.0;
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     //gPad->SetLogy(1);
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     // normalize th2 histogram
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       //name = dir + "h_tt"; name += ipmt; name += ".png";
0427       title = "h_tt"; title += ipmt;
0428     }
0429     else if ( subpass==1 || subpass==2 )
0430     {
0431       //name = dir + "h_ttcorr"; name += ipmt; name += ".png";
0432       title = "h_ttcorr"; title += ipmt;
0433     }
0434     //cout << title << endl;
0435     ac[cvindex]->Print( pdfname, title );
0436   }
0437   cal_tt_t0_file.close();
0438   make_cdbtree( cal_fname );
0439   ++cvindex;
0440 
0441   // Now calculate tq_t0
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     // only look in the middle
0452     if ( ipmt==0 || ipmt==64 )
0453     {
0454       // use wide range for 1st channel in each arm
0455       h_tq[ipmt]->SetAxisRange(-25.,25.);
0456       //h_tq[ipmt]->SetAxisRange(9.,13.); // kludge to select middle bunch during stochastic cooling
0457     }
0458     else
0459     {
0460       // for subsequent channels in an arm, use a window
0461       // determined by the 1st channel
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     //Double_t sigma = 3.0;
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     //gPad->SetLogy(1);
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     // normalize th2 histogram
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       //name = dir + "h_tq"; name += ipmt; name += ".png";
0503       title = "h_tq"; title += ipmt;
0504     }
0505     else if ( subpass==1 || subpass==2 )
0506     {
0507       //name = dir + "h_tqcorr"; name += ipmt; name += ".png";
0508       title = "h_tqcorr"; title += ipmt;
0509     }
0510     //cout << name << endl;
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     //== Draw the slewcorr histograms
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       //name = dir + "h2_slew"; name += ipmt; name += "_pass"; name += subpass; name += ".png";
0533       name = dir + "h2_slew"; name += ipmt; name += "_pass"; name += subpass;
0534       cout << name << endl;
0535       ac[cvindex]->Print( pdfname, name );
0536     }
0537 
0538     // Here we could perhaps add a summary page
0539  
0540     ac[cvindex]->Print( pdfname + "]" );
0541     ++cvindex;
0542   }
0543 
0544   //== Draw the charge histograms
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       //name = dir + "h_adc"; name += ipmt; name += ".png";
0559       title = "h_adc"; title += ipmt;
0560       //cout << pdfname << " " << title << endl;
0561       ac[cvindex]->Print( pdfname, title );
0562     }
0563     ac[cvindex]->Print( pdfname + "]" );
0564     ++cvindex;
0565   }
0566 
0567   savefile->Write();
0568   //savefile->Close();
0569 
0570   if ( subpass==3 )
0571   {
0572     recal_mbd_mip( tfname, subpass, runtype );
0573   }
0574 
0575 }
0576