FitNdNonlinQuad.cxx

Go to the documentation of this file.

//  FitNdNonlinQuad.cxx
//
//  Nathaniel Tagg, Feb 2006
//  tagg@minos.phy.tufts.edu
//
//  Code to read the PulserGain and PulserGainPin tables,
//  perform a quadratic fit to the function, then convert it to
//  an ADC-centric model, create an ntuple for analysis, make
//  plots of each channel for validation, and store the 
//  output in CalPulserFits tables.
//
//  Creates:
//   A root file called ndlin_fits.root, with histograms and trees
//   Fills the ndfits/ directory with plots
//       - summary plots for detector and for each PIN
//       - plot page for each channel if 'drawChannelPlots' is turned on.
//   Fills the CalPulserFits tables with task 44.
//
//  Code is crude and badly documented, but it works.

#include "FitNdNonlinQuad.h"
#include "TGraphErrors.h"
#include "TFile.h"
#include "TF1.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "TPaveText.h"
#include "TROOT.h"
#include "TSystem.h"
#include "TStyle.h"
#include "TTree.h"
#include "TClonesArray.h"
#include "TLine.h"
#include "Record/RecArrayAllocator.h"
#include "Plex/PlexHandle.h"
#include "DatabaseInterface/DbiResultPtr.h"
#include "DatabaseInterface/DbiWriter.h"
#include "PulserCalibration/PulserGain.h"
#include "PulserCalibration/PulserGainPin.h"
#include "PulserCalibration/PulserConventions.h"
#include "Calibrator/CalPulserFits.h"
#include <iostream>
#include <fstream>
#include <vector>
#include <map>
#include <cmath>

using namespace std;

const int db_version = 1; // Change this if the code changes!
// This increments the creation date slightly so newer code makes tables that 
// get treated with priority.

ClassImp(NDQStripfit)
ClassImp(NDQStripPulsePoint)


void FitNdNonlinQuad(
                     VldContext inputContext,
                     const char* outputDirectory,
                     bool writeToDb,
                     bool useHighGainPins,
                     bool drawChannelPlots
                     )
{
  PlexHandle plex(inputContext);
  DbiResultPtr<PulserGain>    fPmtTable(inputContext,1);//Pulser::kNearEnd);
  DbiResultPtr<PulserGainPin> fPinTable(inputContext,1);    

  if(fPmtTable.GetNumRows()==0) {
    cout << "No PMT data available for cx " << inputContext.AsString() << endl;
    return;
  }
  if(fPinTable.GetNumRows()==0) {
    cout << "No PIN data available for cx " << inputContext.AsString() << endl;
    return;
  }


  // Create output directory if it doesn't exist.
  gSystem->MakeDirectory(outputDirectory);
  

  DbiWriter<CalPulserFits>* writer = 0;
  if(writeToDb) {
    // Do this for the first prototype table.
    //VldTimeStamp start(1970,1,1,0,0,0,0);
    //VldTimeStamp stop(2038,1,1,0,0,0,0); // 1 year later.
    //VldTimeStamp creation = start + VldTimeStamp((time_t)db_version,0);

    // Otherwise, use the validity time from the input data.
    VldTimeStamp start = fPmtTable.GetValidityRec(fPmtTable.GetRow(0))->GetVldRange().GetTimeStart();
    VldTimeStamp stop(start.GetSec()+(3*365*24*3600),0); // 3 years later.
    VldTimeStamp creation(start.GetSec() + (int)db_version,0);


    VldRange range((int)inputContext.GetDetector(),
                   SimFlag::kData + SimFlag::kMC,
                   start,
                   stop,
                   "");
    writer = new DbiWriter<CalPulserFits>(range,
                                          -1, // aggragate
                                          44, // task
                                          creation,
                                          "offline", //db
                                          "ND Quadratic fit");                                    
    cout << "Starting DB write with range: " << endl;
    range.Print();
  }

  gStyle->SetPadColor(kWhite);
  
  TCanvas* c1;
  if((c1 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("c1"))) {
  }
  else   c1 = new TCanvas("c1","c1",800,850);
  c1->cd();

  TFile* f = new TFile("ndlin_fits.root","RECREATE");
  f->cd();

  NDQStripfit* fitptr = new NDQStripfit(40);
  NDQStripfit& fit = *fitptr;
  TTree* tree = new TTree("fits","fits");
  tree->Branch("fitbranch","NDQStripfit",&fitptr);

  TH1F* hchi = new TH1F("hchi","ChiSq/dof",200,0,50);
  TH1F* hres = new TH1F("hres","Fractional Residual",500,-1,1);
  TH2F* hres2= new TH2F("hres2","Fractional Residual",
                        500,0,35000,
                        500,-1,1);

  TH2F* hnon2= new TH2F("hnon2","Approx. Nonlinearity",
                        500,0,35000,
                        500,-1,0.3);

  TH2F* htotres2= new TH2F("htotres2","TotalResidual",
                           500,0,35000,
                           500,-1000,1000);


  TH2F* hres2i= new TH2F("hres2i","TotalResidual",
                           40,0,40,
                           100,-1,1);

  std::map<PlexPinDiodeId,TH2*> mres2pin;

  TH1F* ha  = new TH1F("ha","Fit Parameter a (offset)",250,-500,500);
  ha->SetXTitle("ADC");
  ha->SetYTitle("stripends");
  TH1F* hb  = new TH1F("hb","Fit Parameter b (linear)",180,-10,80);
  hb->SetXTitle("ADC/PIN");
  hb->SetYTitle("stripends");
  TH1F* hc  = new TH1F("hc","Fit Parameter c (quad)",200,-0.01,0.01);
  hc->SetXTitle("ADC/PIN^{2}");
  hc->SetYTitle("stripends");
    
  TH1F* hd  = new TH1F("hd","Constant d: PIN offset",300,-300,300);
  hd->SetXTitle("PIN");
  hd->SetYTitle("stripends");
  TH1F* he  = new TH1F("he","Constant e: PIN/ADC conversion",240,-20,100);
  he->SetXTitle("PIN/ADC");
  he->SetYTitle("stripends");
  TH1F* hbeta  = new TH1F("hbeta","Nonlinearity Term #beta",200,-2e-5,1e-5);
  he->SetXTitle("ADC^{-1}");
  he->SetYTitle("stripends");

  
  TH1F* hbadres   = new TH1F("hbadres","Worst Residual",200,0,10000);
  TH1F* hbadresf  = new TH1F("hbadresf","Worst Fractional Residual",200,0,2);
  TH1F* hbadres2500   = new TH1F("hbadres2500","Worst Residual above 2500 ADC",200,0,10000);
  TH1F* hbadresf2500  = new TH1F("hbadresf2500","Worst Fractional Residual above 2500 ADC",200,0,2);
 
  ofstream ofit("fitdump.dat");

  //  for(int stp=10;stp<20;stp++) {
  //  PlexStripEndId seid(Detector::kNear, 10, stp, StripEnd::kWest);

  TGraphErrors* gr =0;
  TGraphErrors* grnon =0;
  TGraphErrors* grres =0;
  TGraphErrors* grgain =0;
  TGraphErrors* grcal  =0;
  TGraphErrors* grzoom =0;
  
  const std::vector<PlexStripEndId>& list = plex.GetAllStripEnds();
  for(UInt_t stp=0;
      stp<list.size();
      //stp<100;
      stp++) {
    PlexStripEndId seid = list[stp];
  
    //for(int stp=0;stp<1;stp++) {
    //  PlexStripEndId seid(Detector::kNear,1,15,StripEnd::kWest);

    cout << "PlexStripEndId: " << seid.AsString() << endl;
    PlexLedId        led = plex.GetLedId(seid);
    PlexPinDiodeId pinid = plex.GetPinDiodeIds(led).first;
    if(useHighGainPins) pinid = plex.GetPinDiodeIds(led).second;

    // Select only the good HG pin diode.
    //if(led.GetPulserBox()!=1) continue;
    //if(led.GetLedInBox() <9 ) continue;
    //if(led.GetLedInBox() >11) continue;
    //if(led.GetLedInBox() != 4) continue;
    

    // Fill out id.
    fit.Reset(40);
    fit.seid = seid;
    fit.plane = seid.GetPlane();
    fit.strip = seid.GetStrip();
    fit.end = seid.GetEnd();
    fit.pb  = led.GetPulserBox();
    fit.led = led.GetLedInBox();
    fit.pin = pinid.GetEncoded();
    //fit.pmt = pmtid.GetEncoded();

    fit.ok = true;
    
    bool isBad = false;

    if(gr)     { delete gr;     gr     =0; }
    if(grnon)  { delete grnon;  grnon  =0; }
    if(grres)  { delete grres;  grres  =0; }
    if(grgain) { delete grgain; grgain =0; }
    if(grcal)  { delete grcal; grcal =0; }
    if(grzoom) { delete grzoom; grzoom =0; }

 
    const PulserGain*    pmt = fPmtTable.GetRowByIndex(seid.BuildPlnStripEndKey());
    const PulserGainPin* pin = fPinTable.GetRowByIndex(pinid.GetEncoded());

    if(!pmt){ cout << "No Pmt data. " << endl; fit.ok = false; tree->Fill(); continue; }
    if(!pin){ cout << "No pin data. " << endl; fit.ok = false; tree->Fill(); continue; }
    // check the pin, which may be bad:
    double pintot = 0;
    for(int i=0;i<pin->GetNumPoints();i++) pintot += (pin->GetMean())[i];
    if(pintot<=0) {
      cout << "Pin data is bad." << endl;
      fit.ok = false;
      tree->Fill();
      continue;
    }

    int n = pin->GetNumPoints();
    fit.n = n;

    // Stripend
    const Float_t* pmtmean = pmt->GetMean(); 
    const Float_t* pmterr  = pmt->GetError(); 
    const Float_t* pmtent  = pmt->GetNumEntries(); 
    const Float_t* pmttrg  = pmt->GetNumTriggers();
    // Pin
    const Float_t* pinmean = pin->GetMean(); 
    const Float_t* pinerr  = pin->GetError(); 
    const Float_t* pinent  = pin->GetNumEntries(); 
    const Float_t* pintrg  = pin->GetNumTriggers();

    Int_t dof = 0;
    Float_t lowerbound = 0;
    Float_t upperbound = 9e9;

    gr = new TGraphErrors(n);    
    for(int i=0;i<n;i++) {
      float x = pinmean[i];
      float dx = pinerr[i] + 4.; // Increase errors on x-coord to improve Chisq.
      float z = pmtmean[i];//*pmtent[i]/pmttrg[i]; // Zero correction.
      float dz = pmterr[i];//*pmtent[i]/pmttrg[i]; // Zero correction.

      if(z<300){
        // lowerbound = x+1.; // exclude point from fit
        // Correct for zeroes, and add a big error to these points:
        float zerocorr = pmtent[i]/pmttrg[i];
        z  = z*zerocorr;
        dz = dz*zerocorr + z*0.5*(zerocorr-1);
      }
      dof++;             // include point

      gr->SetPoint     (i,x, z );
      gr->SetPointError(i,dx,dz); // Increase the error on the PIN, which greatly improves the chisq.      
    }
    
    const char* fcnname = "pol2";
    
    //gr->Draw("a p e0");
    gr->Fit(fcnname,"eq","goff",lowerbound,upperbound);
    //gr->Draw("a p e0");
    TF1* fcn = gr->GetFunction(fcnname);  
    

    double chisq = fcn->GetChisquare();
    double a = fcn->GetParameter(0);
    double b = fcn->GetParameter(1);
    double c = fcn->GetParameter(2);
    double da = fcn->GetParError(0);
    double db = fcn->GetParError(1);
    double dc = fcn->GetParError(2);

    cout << chisq << "/" << dof << endl;

    // Badness based on initial fit:
    if( (1.0 - 4*a*c/(b*b)) < 0 ) {
      // Failed: no intercept!
      isBad=true;
    };

    // d is x-intercept of fit. (PIN count offset)
    double d = (b/(2.*c)) * (-1.0 + sqrt(1.0 - 4*a*c/(b*b)) );
    // e is conversion factor in units of PIN/ADC
    double e = 1.0/(b+2.*d*c);
    double de = (db/b)*e; // roughly correct.
    // so, z = y + y*y*beta
    double beta = e*e*c;
    double dbeta = (dc/c)*beta;

    fit.dof = dof;
    fit.a = a;
    fit.b = b;
    fit.c = c;
    fit.d = d;
    fit.e = e;
    fit.beta = beta;
    fit.da = da;
    fit.db = db;
    fit.dc = dc;
    fit.dd = 0;
    fit.de = de;
    fit.dbeta = dbeta;

    // Badness based on parameter results:
    if(beta > -0) isBad = true;

    hchi->Fill(chisq/(float)n);
    ha->Fill(a);
    hb->Fill(b);
    hc->Fill(c);
    hd->Fill(d);
    he->Fill(e);
    hbeta->Fill(beta);


    TH2* hres2pin = 0;
    std::map<PlexPinDiodeId,TH2*>::iterator hit = mres2pin.find(pinid);
    if(hit == mres2pin.end()) {
      hres2pin = 
        new TH2F(Form("hpin_PB%02d_LED%02d_%c",led.GetPulserBox(),led.GetLedInBox(),(pinid.IsLowGain()?'L':'H')),
                 Form("hpin %s (%c)",led.AsString(),(pinid.IsLowGain()?'L':'H')),
                 //100,0,35000,
                 40,0,40,
                 100,-1,1);
      mres2pin[pinid] = hres2pin;
    } else {             
      hres2pin = hit->second;
    }


    float maxresidual = 0;
    float maxfracresidual = 0;
    float maxresidual2500 = 0;
    float maxfracresidual2500 = 0;

    float meanres2500 = 0;
    float meanres = 0;   
    float nres2500 = 0;
    float adcmax = 0;
    
    if(drawChannelPlots) {
      grres = new TGraphErrors(n);
      grnon = new TGraphErrors(n);
      grgain = new TGraphErrors(n);
      grcal = new TGraphErrors(n);
      grzoom = new TGraphErrors(n);
    }

    for(int i=0;i<n;i++) {
      // x is pin counts
      double x = pinmean[i];
      double dx= pinerr[i] + 4.;
      
      // y is linearized ADC counts.
      double y = (x-d)/e;
      double dy= dx/e;

      double zerocorr = pmtent[i]/pmttrg[i]; // Zero correction.

      // z is nonlinear ADC counts:
      double z = pmtmean[i]; 
      double dz = pmterr[i];
      if(z<300){
        // Correct for zeroes, and add a big error to these points:
        z  = z*zerocorr;
        dz = dz*zerocorr + 0.5*z*(1-zerocorr);
      }

      // The fit at this point (to be compared with z, the data)
      double f = fcn->Eval(x);

      // Do a trial calibration:
      double s = 1+4.*z*beta;
      if(s<0) s=0;
      double  caly =  1.0/(2.*beta) * ( sqrt(s) - 1. );

      double dcaly =  0;
      if(s>0) dcaly = 2.0 /  sqrt(s) *  dz;

      double res = z-f;
      double fres = (z-f)/f;
      double gain = (dz*dz*pmtent[i])/z * 0.8;
      double pe   = z/gain;
      
      // Horizontal residual:
      double xres = 9999;
      double temp = b*b - 4.*c*(a-z);
      if(temp>0) {
        temp = sqrt(temp);
        double x1 = (-b + temp)/(2.*c);
        double x2 = (-b - temp)/(2.*c);
        if(fabs(x-x1) < fabs(x-x2))
          xres = x-x1;
        else
          xres = x-x2;
      }
      double fxres = xres/x;
      
      NDQStripPulsePoint& datum = *((NDQStripPulsePoint*) fit.data->operator[](i));
      datum.i  =i;
      datum.x  =x;
      datum.dx =dx;
      datum.y  =y;
      datum.dy =dy;
      datum.z  =z;
      datum.dz =dz;
      datum.f  =f;
      datum.caly  = caly;
      datum.dcaly =dcaly;
      datum.res = res;
      datum.fres= fres;
      datum.xres = xres;
      datum.fxres = fxres;
      datum.gain = gain;
      datum.pe   = pe;
      datum.zerocorr = zerocorr;
      
      hres->Fill(fres);
      hres2->Fill(y,fres);
      hnon2->Fill(y,(z-y)/y);
      htotres2->Fill(y,res);

      //hres2i->  Fill(i,fres);
      //hres2pin->Fill(i,fres);
      hres2i->  Fill(i,fxres);
      hres2pin->Fill(i,fxres);


      if( fabs(res) > maxresidual) maxresidual = fabs(res);
      if( fabs(fres) > maxfracresidual) maxfracresidual = fabs(fres);
      meanres += fabs(res);
      if(y>2500) {
        if( fabs(res)  > maxresidual2500)       maxresidual2500 = fabs(res);
        if( fabs(fres) > maxfracresidual2500) maxfracresidual2500 = fabs(fres);
        meanres2500 += fabs(res);
        nres2500 += 1;
      }
      if(y>adcmax) adcmax = y;
 
      if(drawChannelPlots) {
        grres->SetPoint     (i,y,  (z - f)/f );
        grres->SetPointError(i,dy, dz/f);

        grnon->SetPoint     (i,y,  (z-y)/y);
        grnon->SetPointError(i,dy, dz/y);      
        
        grgain->SetPoint(i,y, (dz*dz*pmtent[i])/pmtmean[i] * 0.8 );

        grcal->SetPoint     (i,y, caly/y);
        grcal->SetPointError(i,dy,dcaly/y);

        grzoom->SetPoint     (i,y, z);
        grzoom->SetPointError(i,fabs(1-pinent[i]/pintrg[i])*y,dz);
      }
      
      //cout << x << "\t" << y << "\t" << z << "\t" << caly
      //           << "\t" << y+y*y*beta << endl;    
    }

    
    hbadres->Fill(maxresidual);
    hbadresf->Fill(maxfracresidual);
    hbadres2500->Fill(maxresidual2500);
    hbadresf2500->Fill(maxfracresidual2500);
    if(nres2500>0) meanres2500 /= nres2500;
    if(n>0) meanres/=n;
    if(chisq/(float)(dof)>75) isBad = true;

    fit.chisq = chisq;
    fit.maxres = maxresidual;
    fit.maxresf= maxfracresidual;
    fit.maxres2500  = maxresidual2500;
    fit.maxresf2500 = maxfracresidual2500;
    fit.meanres2500 = meanres2500;
    fit.adcmax = adcmax;
    fit.bad = isBad;

    if(isBad) {      
      ofit << Form("%s/p%03ds%03d.gif",outputDirectory,seid.GetPlane(),seid.GetStrip()) << endl;
    /*
      ofit << Form("%3d %3d %5.1f %2d  ",
      seid.GetPlane(),
      seid.GetStrip(),
      fcn->GetChisquare(),
      n);
      ofit << Form("  %.1f %.1f  %.1f %.1f  %.4f %.4f",
                   fcn->GetParameter(0), fcn->GetParError(0),
                   fcn->GetParameter(1), fcn->GetParError(1),
                   fcn->GetParameter(2), fcn->GetParError(2) );
      ofit << Form("  %.1f %.2f",maxresidual, maxfracresidual);
      ofit << endl;
    */
    }

    if(drawChannelPlots) {
      c1->cd();
      gPad->Clear();
      TVirtualPad* pad=gPad;
      TVirtualPad* spad1 = new TPad("spad1","spad1",0,  0,0.8,1);
      TVirtualPad* spad2 = new TPad("spad2","spad2",0.75,0,1,  1);
      spad1->Draw();
      spad2->Draw();
      
      // Text.
      spad2->cd();
      TPaveText* txt = new TPaveText(0,0.1,0.95,0.9,"NDC T");
      txt->SetFillColor(kWhite);
      txt->SetTextAlign(12);
      txt->SetTextSize(0.05);
      txt->AddText(seid.AsString());
      txt->AddText(Form("Plane: %d  Strip: %d",seid.GetPlane(),seid.GetStrip()));
      txt->AddText(plex.GetPixelSpotId(seid).AsString("p"));
      txt->AddText(plex.GetRawChannelId(seid).AsString("e"));
      txt->AddText(led.AsString("e"));
      txt->AddText(Form("%s (%02d%02d%01d)",pinid.AsString("e"),pinid.GetInRack(),pinid.GetInBox(),pinid.GetGain()));
      txt->AddText(Form("#{Chi}^{2}/DoF = %.1f/%d",chisq,dof));
      txt->AddText(Form("a: %.1f +/- %.1f ADC",a,da));
      txt->AddText(Form("b: %.1f +/- %.1f ADC/PIN",b,db));
      txt->AddText(Form("c: %.1e +/- %.1e ADC/PIN^{2}",c,dc));
      txt->AddText(" ");
      txt->AddText(Form("d: %.2f PIN",d));
      txt->AddText(Form("e: %.3f PIN/ADC",e));
      txt->AddText(Form("#beta: %.2e ADC^{-1}",beta));
      txt->Draw();      
      
      spad1->Divide(1,4,0.001,0.001,kWhite);
      spad1->cd(1); 
      gr->Draw("a p e0");
      gr->GetHistogram()->SetTitle("Full Curve with fit");
      gr->GetHistogram()->SetXTitle("PIN counts");
      gr->GetHistogram()->SetYTitle("PMT counts");
      
      spad1->cd(2);
      grnon->Draw("a p e0");
      grnon->GetHistogram()->SetTitle("Nonlinearity");
      grnon->GetHistogram()->SetXTitle("Linear PMT counts");
      grnon->GetHistogram()->SetYTitle("(Nonlinear - Linear)/Linear");
      
      spad1->cd(3);
      grres->Draw("a p e0");
      grres->GetHistogram()->SetTitle("Fractional Fit Residual");
      grres->GetHistogram()->SetXTitle("Linear PMT counts");
      grres->GetHistogram()->SetYTitle("(Data - Fit)/Fit");
      
      spad1->cd(4);
      //grgain->SetMarkerStyle(20);
      //grgain->Draw("a p");
      //grgain->GetHistogram()->SetTitle("Apparent gain");
      //grgain->GetHistogram()->SetXTitle("Linear PMT counts");
      //grgain->GetHistogram()->SetYTitle("PEs");
      //grcal->Draw("a p e0");
      //grcal->GetHistogram()->SetTitle("Calibrated result");
      //grcal->GetHistogram()->SetXTitle("Approximate linear ADC");
      //grcal->GetHistogram()->SetYTitle("Calibrated ADC / Pin-esimated ADC");
      static TH2* hframezoom = 0;
      static TLine* lzoom = 0;
      if(hframezoom == 0){
        hframezoom = new TH2F("hframezoom","Zoom",20,0,6000,20,0,6000);
        hframezoom->SetXTitle("Linear PMT counts");
        hframezoom->SetYTitle("Nonlinear PMT counts");
        lzoom = new TLine(0,0,6000,6000);
      }
      hframezoom->Draw();
      lzoom->Draw();
      grzoom->Draw("e0");
      
    
      pad->Update();
      pad->cd();
      
      c1->Print(Form("%s/p%03ds%03d.gif",outputDirectory,seid.GetPlane(),seid.GetStrip()));
    }
    
    tree->Fill();
    
    
    if(writeToDb) {     
      // Adjust for bad fits:
      if(beta>0) {
        dbeta +=fabs(beta);
        beta = 0;
      }
      if( (1.0 - 4*a*c/(b*b)) < 0 ) {
        d = 0;
      };

      CalPulserFits row(-1, //aggno
                        seid.BuildPlnStripEndKey(),
                        0, // fit type
                        (TF1*)0, 0, 0, 0);
      row.SetFitType(44);
      row.SetNPFit(dof);
      row.SetSlope(e);
      row.SetSlopeErr(de);
      row.SetXOffset(d);
      row.SetPar1(beta);
      row.SetPar2(dbeta);
      row.SetChisq(chisq);
      if(nres2500>0) {
        row.SetMeanRes(meanres2500);
        row.SetMaxRes(maxresidual2500);
      } else {
        row.SetMeanRes(meanres);
        row.SetMaxRes(maxresidual);
      }
      row.SetAdcMax(adcmax);
      //row.Dump();
      (*writer) << row;
    }
    
  }
  
  if(writeToDb) {
    writer->Close();
    delete writer;
  }

  tree->Write();
  hchi->Write();
  hres->Write();
  htotres2->Write();

  hnon2->Write();

  c1->cd(); c1->Clear();
  hres2->Draw("colz");
  c1->Update();
  c1->Print(Form("%s/_res2.gif",outputDirectory));
  hres2->Write();

  hres2i->Draw("colz");
  c1->Print(Form("%s/_res2i.gif",outputDirectory));
  c1->Update();
  hres2i->Write();

  ha->Write();
  hb->Write();
  hc->Write();
  hd->Write();
  he->Write();
  hbeta->Write();
  hbadres->Write();
  hbadresf->Write();
  hbadres2500->Write();
  hbadresf2500->Write();
  for(std::map<PlexPinDiodeId,TH2*>::iterator hit = mres2pin.begin();
      hit!=mres2pin.end();
      hit++) {
    TH2* h = hit->second;
    c1->cd(); c1->Clear();
    h->Draw("colz");
    c1->Update();
    c1->Print(Form("%s/_%s.gif",outputDirectory,h->GetName()));
    h->Write();
  }
  f->Write();
}

---