PulserGainFit.cxx

Go to the documentation of this file.

// $Id: PulserGainFit.cxx,v 1.11 2007/01/23 17:01:56 rhatcher Exp $
//
// PulserGainFit
//
// Perform linear fits to pulser gain data and store results
// in CalPulserFits table in database.
//
// Author: Pat Ward 13/09/2005
//
// $Log: PulserGainFit.cxx,v $
// Revision 1.11  2007/01/23 17:01:56  rhatcher
// Ongoing changes to ROOT headers mean we need more explicit #include's
// in various user classes (ie. can't rely on implicit includes).
//
// Revision 1.10  2006/09/15 13:18:20  cpw1
// Make min/max aggregate configurable
//
// Revision 1.9  2006/08/23 13:48:22  cpw1
// Add setable DBSimMask and DBCommment
//
// Revision 1.8  2006/04/19 20:00:55  rhatcher
// convert DetectorType:: to Detector:: everywhere.
//
// Revision 1.7  2006/04/07 22:02:06  tagg
// Add <math.h> to make gcc 4.0 happy
//
// Revision 1.6  2006/03/29 11:21:50  cpw1
// Updates for 60-point gain curve
//
// Revision 1.5  2006/01/27 09:19:20  cpw1
// Fix bug in loop index in FarMeans
//
// Revision 1.4  2006/01/25 14:07:48  cpw1
// Mods for new definition of CalPulserFits; split Run method into RunPinFits and RunNearFarFits and add code for near-far fits
//
// Revision 1.3  2006/01/12 20:01:43  murgia
// Abot DbiWriter only if last pass.
// Use SetRange to define fit interval.
//
// Revision 1.2  2005/10/27 09:35:21  cpw1
// Use zero-corrected values in fits and add number of fit points to CalPulserFits
//
// Revision 1.1  2005/09/28 09:39:59  cpw1
// Fill CalPulserFits table
//
//

#include "TH1.h"

#include "MessageService/MsgService.h"
#include "Calibrator/CalVaLinearity.h"
#include "PulserCalibration/PulserGainFit.h"
#include "DatabaseInterface/DbiResultSet.h"
#include "DatabaseInterface/DbiValidityRec.h"
#include "Plex/PlexPinDiodeId.h"
#include "Plex/PlexStripEndId.h"
#include <math.h>

ClassImp(PulserGainFit)

CVSID("$Id: PulserGainFit.cxx,v 1.11 2007/01/23 17:01:56 rhatcher Exp $\n  \
      CVSID_DBIRESULTPTR ");

// Constructors and destructor
PulserGainFit::PulserGainFit(const Detector::Detector_t &det, VldTimeStamp &tsin, VldTimeStamp &tein)
{
  // With this constructor we set the validity range for the output to be
  // the same as that in the input database entry for each aggregate
  fUseSameValidity = 1;

  // Set detector and timestamps
  fDetector = det;
  fStartin = tsin;
  fEndin = tein;
  MSG("Pulser",Msg::kInfo)<<"PulserGainFit created for detector "<<Detector::AsString(fDetector)<<endl;
  MSG("Pulser",Msg::kInfo)<<"Input timestart = "<<fStartin.AsString()<<endl;
  MSG("Pulser",Msg::kInfo)<<"Input timeend = "<<fEndin.AsString()<<endl;
  MSG("Pulser",Msg::kInfo)<<"Output timestart and timend will be taken from validity range of input data"<<endl;

  // Initialize output times so can create a range, but these values are not used
  fStartout = fStartin;
  fEndout = fEndin;

  this->Init();

}

PulserGainFit::PulserGainFit(const Detector::Detector_t &det, VldTimeStamp &tsin, VldTimeStamp &tein, VldTimeStamp &tsout, VldTimeStamp &teout)
{
  // With this constructor we set the validity range for the output is set
  // by the user
  fUseSameValidity = 0;

  // Set detector and timestamps
  fDetector = det;
  fStartin = tsin;
  fEndin = tein;
  fStartout = tsout;
  fEndout = teout;
  MSG("Pulser",Msg::kInfo)<<"PulserGainFit created for detector "<<Detector::AsString(fDetector)<<endl;
  MSG("Pulser",Msg::kInfo)<<"Input timestart = "<<fStartin.AsString()<<endl;
  MSG("Pulser",Msg::kInfo)<<"Input timeend = "<<fEndin.AsString()<<endl;
  MSG("Pulser",Msg::kInfo)<<"Output timestart = "<<fStartout.AsString()<<endl;
  MSG("Pulser",Msg::kInfo)<<"Output timeend = "<<fEndout.AsString()<<endl;

  this->Init();

}

PulserGainFit::~PulserGainFit() {
  if(fFile) {
    fFile->Write();
    fFile->Close();
  }
}

// Initialization
void PulserGainFit::Init() {
  // Set default configuration variables
  Registry r;
  r.Set("WriteDB",0);
  r.Set("DBSimMask",5);
  r.Set("DBName","temp");
  r.Set("DBComment","");
  r.Set("WritePlots",0);
  r.Set("PlotsFile","plots.root");
  r.Set("Adcmin",300.);
  r.Set("Adcmax",8000.);
  r.Set("Npass",1);
  r.Set("MinAggNo",-1);
  r.Set("MaxAggNo",-1);
  InitializeConfig(r); 

  // Create fit functions
  // Line
  fFitLine = new TF1("line","[0]*(x-[1])");
  fFitLine->FixParameter(1,0.);

  // Pole only 
  fFitPole = new TF1("pole",CalVaLinearity::FitFunction,200.,15000.,8);
  // First parameter is fit function type
  fFitPole->FixParameter(0,5.);
  // Parameter 5 is not used
  fFitPole->FixParameter(5,0.);
  // Fix pedestal to zero 
  fFitPole->FixParameter(6,0.);
  // k3, k4 not used for pole only
  fFitPole->FixParameter(3,0.);
  fFitPole->FixParameter(4,0.);

  // Pole+kicker with fixed pedestal
  fFitPK1 = new TF1("pk1",CalVaLinearity::FitFunction,200.,15000.,8);
  // First parameter is fit function type
  fFitPK1->FixParameter(0,25.);
  // Parameter 5 is not used
  fFitPK1->FixParameter(5,0.);
  // Fix pedestal to zero 
  fFitPK1->FixParameter(6,0.);

  // Pole+kicker with free pedestal
  fFitPK2 = new TF1("pk2",CalVaLinearity::FitFunction,200.,15000.,8);
  // First parameter is fit function type
  fFitPK2->FixParameter(0,25.);
  // Parameter 5 is not used
  fFitPK2->FixParameter(5,0.);

  fFile = 0;
  fNextra = 0;
}

// Set configuration
void PulserGainFit::ConfigModified()
{
  fWriteDB = GetConfig().GetInt("WriteDB");
  fDBSimMask = GetConfig().GetInt("DBSimMask");
  fDBName = GetConfig().GetCharString("DBName");
  fDBComment = GetConfig().GetCharString("DBComment");
  fWritePlots = GetConfig().GetInt("WritePlots");
  fPlotsFile = GetConfig().GetCharString("PlotsFile");
  GetConfig().Get("Adcmin",fAdcmin);
  GetConfig().Get("Adcmax",fAdcmax);
  GetConfig().Get("Npass",fNpass);
  GetConfig().Get("MinAggNo",fMinAggNo);
  GetConfig().Get("MaxAggNo",fMaxAggNo);

  fWriterl.SetDbName(fDBName);
  fWriterh.SetDbName(fDBName);
  fWriternf.SetDbName(fDBName);

  MSG("Pulser",Msg::kInfo)<<"Configuration: WriteDB = "<<fWriteDB<<" WritePlots = " <<fWritePlots<<" Adcmin = "<<fAdcmin<<" Adcmax = "<<fAdcmax<<" Npass = "<<fNpass<<endl;
  if(fWriteDB > 0) {MSG("Pulser",Msg::kInfo)<<"Writing to database "<<fDBName<<" with SimMask "<<fDBSimMask<<" Comment "<<fDBComment<<endl;}
  if(fWritePlots > 0) {MSG("Pulser",Msg::kInfo)<<"Writing plots to "<<fPlotsFile<<endl;}
}

// Set timestamps for including extra runs for split gain curves
void PulserGainFit::AddRun(VldTimeStamp &ts, VldTimeStamp &te) {

  // Check max runs not exceeded
  if(fNextra == 3) {
  MSG("Pulser",Msg::kError)<<"Maximum runs exceeded, ignoring run"<<endl;
  return;
  }

  // Warning if trying to use same validity context for output
  if(fUseSameValidity) {
  MSG("Pulser",Msg::kWarning)<<"Output validity will not be changed by call to AddRun"<<endl;
  }

  // Set time stamps
  fStartExtra[fNextra] = ts;
  fEndExtra[fNextra] = te;
  fNextra++;
}

// Run stripend v Pin fits on all data and store values
void PulserGainFit::RunPinFits() {

  // Create file for plots and book some histograms
  if(fWritePlots) {
    if(fFile == 0) {
      fFile = new TFile(GetConfig().GetCharString("PlotsFile"),"RECREATE");
    }
    hSlopeNL = new TH1F("hSlopeNL","hSlopeNL",100,0.,20.);
    hSlopeNH = new TH1F("hSlopeNH","hSlopeNH",100,0.,20.);
    hSlopeFL = new TH1F("hSlopeFL","hSlopeFL",100,0.,20.);
    hSlopeFH = new TH1F("hSlopeFH","hSlopeFH",100,0.,20.);
    hChisqNL = new TH1F("hChisqNL","hChisqNF",100,0.,20.);
    hChisqNH = new TH1F("hChisqNH","hChisqNF",100,0.,20.);
    hChisqFL = new TH1F("hChisqFL","hChisqNF",100,0.,20.);
    hChisqFH = new TH1F("hChisqFH","hChisqNF",100,0.,20.);
    hNumPointsNL = new TH1F("hNumPointsNL","hNumPointsNL",20,0.,20.);
    hNumPointsNH = new TH1F("hNumPointsNH","hNumPointsNH",20,0.,20.);
    hNumPointsFL = new TH1F("hNumPointsFL","hNumPointsFL",20,0.,20.);
    hNumPointsFH = new TH1F("hNumPointsFH","hNumPointsFH",20,0.,20.);
  }

  // Set validity range for output data
  // Note that this may be changed later
  VldRange range(fDetector,SimFlag::kData,fStartout,fEndout,fDBName);

  // Get plex
  VldContext vldc(fDetector,SimFlag::kData,fStartin);
  PlexHandle plex(vldc);

  // Retrieve PIN data and store pointers in maps
  GetPinData();

  // Instantiate DbiResultPtr for stripend data
  DbiResultPtr<PulserGain> rsStripEnd; 

  // Set min and max aggregate number according to detector type,
  // unless already set by user
  if (fMinAggNo < 0) {
    if (fDetector==Detector::kNear) fMinAggNo = 1;
    if (fDetector==Detector::kFar) fMinAggNo = 1;
  }
  MSG("Pulser",Msg::kInfo)<<"Set min aggregate to " << fMinAggNo << " for "<<Detector::AsString(fDetector)<< " detector." <<endl;
  if (fMaxAggNo < 0) {
    if (fDetector==Detector::kNear) fMaxAggNo = 139;
    if (fDetector==Detector::kFar) fMaxAggNo = 980;
  }
  MSG("Pulser",Msg::kInfo)<<"Set max aggregate to " << fMaxAggNo << " for "<<Detector::AsString(fDetector)<< " detector." <<endl;

  // Loop over near/far ends and fit PMT v PIN
  for (int task = 1; task <=2; task++) {

    // Loop over aggregates
    for (int aggregateno = fMinAggNo; aggregateno <=fMaxAggNo; aggregateno++ ) {
      // Task = 1: Near end v low gain pin
      //        2: Far end v low gain pin
      //        3: Near end v high gain pin
      //        4: Far end v high gain pin

      // Set extended validity context for input data and get data
      DbiSqlContext inContext(DbiSqlContext::kStarts,fStartin,fEndin,fDetector,SimFlag::kData);
      inContext << " and AGGREGATENO=" << aggregateno;
      MSG("Pulser",Msg::kInfo)<<"Aggregate " << aggregateno << endl;
      rsStripEnd.NewQuery(inContext,task); 
      const int numStripEnd = rsStripEnd.GetNumRows();
      MSG("Pulser",Msg::kInfo)<<"Retrieved " << numStripEnd << " stripend rows" << endl;
      if (numStripEnd == 0 ) {continue;}

      // Change validity range of output to be that of input 
      if(fUseSameValidity) {
        const DbiValidityRec* vrec = rsStripEnd.GetValidityRec();
        const VldRange& vrng = vrec->GetVldRange();
        //vrng.Print();
        fStartout = vrng.GetTimeStart();
        fEndout = fStartout;
        fEndout.Add(VldTimeStamp(3600*24*62,0));
        range.SetTimeStart(fStartout);
        range.SetTimeEnd(fEndout);
      }

      // Allow for two passes
      Float_t lxoff = 0.;
      Float_t hxoff = 0.;
      for (int ipass = 1; ipass <= fNpass; ipass++) {
        // Free intercept if this is not last pass
        if(ipass < fNpass) {fFitLine->ReleaseParameter(1);}

        // Open DbiWriters 
        if(fWriteDB > 0 && ipass == fNpass) {
          fWriterl.Open(range,aggregateno,task,0,fDBName);
          fWriterh.Open(range,aggregateno,task+2,0,fDBName);
        }

        // Loop over rows
        for (int ir=0; ir<numStripEnd; ir++) {
         //Pointer to row
         const PulserGain* rs = rsStripEnd.GetRow(ir);
         int skey = rs->GetStripEnd();
         PlexStripEndId seid = PlexStripEndId::UnbuildPlnStripEndKey(fDetector, skey);

         // Fit with low gain pin; fix intercept if last pass
         if(ipass == fNpass) {fFitLine->FixParameter(1,lxoff);}
         Fit(*(lgpinmap.find(aggregateno)->second),*(const_cast<PulserGain*>(rs)));
         float lslope = fFitLine->GetParameter(0);
         float lslopeerr = fFitLine->GetParError(0);
         float lxoffset = fFitLine->GetParameter(1);
         float lchisq = fFitLine->GetChisquare();
         int lnumpoints = fFitLine->GetNumberFitPoints();
         MSG("Pulser",Msg::kDebug)<<"Stripend "<<seid<<" Pin "<<lgpinmap.find(aggregateno)->second->GetPinDiodeId()<<" fit slope= "<<lslope<< " xoffset= "<<lxoffset<<"  numpoints = "<<lnumpoints<<endl;

         if(fWritePlots) {
           if(task==1) {
             hSlopeNL->Fill(lslope);
             hChisqNL->Fill(lchisq);
             hNumPointsNL->Fill(lnumpoints);
           }
           if(task==2) {
             hSlopeFL->Fill(lslope);
             hChisqFL->Fill(lchisq);
             hNumPointsFL->Fill(lnumpoints);
           }
         }

         if (ipass < fNpass) {
           // For mean x offset
           lxoff += lxoffset;
         }
         else {
           // Add row to database
           if(fWriteDB) {
             CalPulserFits lrow(aggregateno,skey,lnumpoints,lslope,lslopeerr,lxoffset,lchisq);
             fWriterl << lrow;
           }
         }

         // Fit with high gain pin; fix intercept if last pass
         if(ipass == fNpass) {fFitLine->FixParameter(1,hxoff);}
         Fit(*(hgpinmap.find(aggregateno)->second),*(const_cast<PulserGain*>(rs)));
         float hslope = fFitLine->GetParameter(0);
         float hslopeerr = fFitLine->GetParError(0);
         float hxoffset = fFitLine->GetParameter(1);
         float hchisq = fFitLine->GetChisquare();
         int hnumpoints = fFitLine->GetNumberFitPoints();
         MSG("Pulser",Msg::kDebug)<<"Stripend "<<seid<<" Pin "<<hgpinmap.find(aggregateno)->second->GetPinDiodeId()<<" fit slope= "<<hslope<< " xoffset= "<<hxoffset<<"  numpoints = "<<hnumpoints<<endl;

         if(fWritePlots) {
           if(task==1) {
             hSlopeNH->Fill(hslope);
             hChisqNH->Fill(hchisq);
             hNumPointsNH->Fill(hnumpoints);
           }
           if(task==2) {
             hSlopeFH->Fill(hslope);
             hChisqFH->Fill(hchisq);
             hNumPointsFH->Fill(hnumpoints);
           }
         }
         if (ipass < fNpass) {
           // For mean x offset
           hxoff += hxoffset;
         }
         else {
           // Add row to database
           if(fWriteDB) {
             CalPulserFits hrow(aggregateno,skey,hnumpoints,hslope,hslopeerr,hxoffset,hchisq);
             fWriterh << hrow;
           }
         }
        }
        lxoff = lxoff/(float)numStripEnd;
        hxoff = hxoff/(float)numStripEnd;
      }
    }
  }

  if(fWriteDB) {
    if(fWriterl.IsOpen()) {fWriterl.Close();}
    if(fWriterh.IsOpen()) {fWriterh.Close();}
  }

}

// Retrieve pin data and store pointers
void PulserGainFit::GetPinData() {

  // Set extended validity context for input data
  DbiSqlContext inContext(DbiSqlContext::kStarts,fStartin,fEndin,fDetector,SimFlag::kData);

  // Retrieve pin data
  DbiResultPtr<PulserGainPin> rsPin("PULSERGAINPIN",inContext,Dbi::kAnyTask); 
  const int numPin = rsPin.GetNumRows();
  MSG("Pulser",Msg::kInfo)<<"Retrieved " << numPin << " pin rows" << endl;

  // Loop over pin rows and store pointers in maps
  for (int ir=0; ir<numPin; ir++) {
    //Pointer to row
    const PulserGainPin* rs = rsPin.GetRow(ir);
    int pingain = rs->GetPinDiodeId().GetGain();
    int aggregateno = rs->GetAggregateNo();
    //If low gain (gain = 1) 
    if (pingain == 1) {
      lgpinmap.insert(pair<int,PulserGainPin*>(aggregateno,const_cast<PulserGainPin*>(rs) ));
    }
    //If high gain (gain = 0)
    if (pingain == 0) { 
      hgpinmap.insert(pair<int,PulserGainPin*>(aggregateno,const_cast<PulserGainPin*>(rs) ));
    }
  } // end loop over pin rows

}

// Perform fit to stripend v pin
void PulserGainFit::Fit(PulserGainPin &pin, PulserGain &stripend){
  // Fit stripend v pin. N.B. arguments reversed from PulserXScale
  // as (x,y) corresponds to convention in e.g. TH2, TGraph etc.

  assert(stripend.GetNumPoints()==pin.GetNumPoints());
  assert(stripend.GetAggregateNo()==pin.GetAggregateNo());

  // Points for fit
  //  const Float_t *x = pin.GetMean();
  //  const Float_t *y = stripend.GetMean();
  //  const Float_t *er_x = pin.GetError();
  //  const Float_t *er_y = stripend.GetError();
  Float_t x[40];
  Float_t y[40];
  Float_t er_x[40];
  Float_t er_y[40];
  for(int i=0; i<pin.GetNumPoints(); i++) {
    x[i] = pin.ZCMean(i);
    y[i] = stripend.ZCMean(i);
    er_x[i] = pin.ZCError(i);
    er_y[i] = stripend.ZCError(i);
  }
  TGraphErrors g(pin.GetNumPoints(),x,y,er_x,er_y);

  // Set range of fit
  float xmin = fAdcmin;
  float xmax = fAdcmax;
  for (int i=0;i<pin.GetNumPoints();i++) {
    if (y[i] < fAdcmin && x[i] > fAdcmin) xmin =x[i]+1;
    if (y[i] > fAdcmax && x[i] < fAdcmax) {
      xmax = x[i]-1;
      break;
    }
  }

  // Set range (verifies that Adcmin<Adcmax) (SM)
  fFitLine->SetRange(fAdcmin,fAdcmax);
  // Initialize slope (to prevent a silly value causing problems for
  // subsequent fits
  fFitLine->SetParameter(0,1.);

  // Perform fit
  g.Fit(fFitLine,"Q"," ",xmin,xmax);

}

// Run Near v Far fits on all data and store values
void PulserGainFit::RunNearFarFits() {

  // Open file and create ntuple
  if(fWritePlots) {
    if(fFile==0) {
      fFile = new TFile(GetConfig().GetCharString("PlotsFile"),"RECREATE");
    }
    fNearFarTree = new TTree("nftree","Near-Far fits");
    fNearFarTree->Branch("plane",&fPlane,"plane/I");
    fNearFarTree->Branch("strip",&fStrip,"strip/I");
    fNearFarTree->Branch("bay",&fBay,"bay/I");
    fNearFarTree->Branch("rack",&fRack,"rack/I");
    fNearFarTree->Branch("tube",&fTube,"tube/I");
    fNearFarTree->Branch("pixel",&fPixel,"pixel/I");
    fNearFarTree->Branch("spot",&fSpot,"spot/I");
    fNearFarTree->Branch("pb",&fPb,"pb/I");
    fNearFarTree->Branch("led",&fLed,"led/I");
    fNearFarTree->Branch("adcmin",&fMinadc,"adcmin/F");
    fNearFarTree->Branch("adcmax",&fMaxadc,"adcmin/F");
    fNearFarTree->Branch("ftype",&fFtype,"ftype/I");
    fNearFarTree->Branch("npar",&fNpar,"npar/I");
    fNearFarTree->Branch("npfit",&fNpfit,"npfit/I");
    fNearFarTree->Branch("gain",&fGain,"gain/F");
    fNearFarTree->Branch("ped",&fPed,"ped/F");
    fNearFarTree->Branch("dped",&fDped,"dped/F");
    fNearFarTree->Branch("k1",&fK1,"k1/F");
    fNearFarTree->Branch("k2",&fK2,"k2/F");
    fNearFarTree->Branch("k3",&fK3,"k3/F");
    fNearFarTree->Branch("k4",&fK4,"k4/F");
    fNearFarTree->Branch("chisq",&fChisq,"chisq/F");
    fNearFarTree->Branch("meanres",&fMeanres,"meanres/F");
    fNearFarTree->Branch("maxres",&fMaxres,"maxres/F");
    fNearFarTree->Branch("delta200",&fDelta200,"delta200/F");
    fNearFarTree->Branch("delta7500",&fDelta7500,"delta7500/F");
  }

  // Set validity range for output data
  // Note that this may be changed later
  VldRange range(fDetector,fDBSimMask,fStartout,fEndout,fDBName);

  // Get plex
  VldContext vldc(fDetector,SimFlag::kData,fStartin);
  PlexHandle plex(vldc);

  // Instantiate DbiResultPtr for stripend data (near and far)
  DbiResultPtr<PulserGain> rsStripEndf[3]; 
  DbiResultPtr<PulserGain> rsStripEndn[3];

  // Set min and max aggregate number according to detector type,
  // unless already set by user
  if (fMinAggNo < 0) {
    if (fDetector==Detector::kNear) fMinAggNo = 1;
    if (fDetector==Detector::kFar) fMinAggNo = 1;
  }
  MSG("Pulser",Msg::kInfo)<<"Set min aggregate to " << fMinAggNo << " for "<<Detector::AsString(fDetector)<< " detector." <<endl;
  if (fMaxAggNo < 0) {
    if (fDetector==Detector::kNear) fMaxAggNo = 139;
    if (fDetector==Detector::kFar) fMaxAggNo = 980;
  }
  MSG("Pulser",Msg::kInfo)<<"Set max aggregate to " << fMaxAggNo << " for "<<Detector::AsString(fDetector)<< " detector." <<endl;

  // Loop over aggregates
  for (int aggregateno = fMinAggNo; aggregateno <=fMaxAggNo; aggregateno++ ) {
    DbiSqlContext inContext(DbiSqlContext::kStarts,fStartin,fEndin,fDetector,SimFlag::kData);
    inContext << " and AGGREGATENO=" << aggregateno;
    MSG("Pulser",Msg::kInfo)<<"Aggregate " << aggregateno << endl;
    rsStripEndn[0].NewQuery(inContext,Pulser::kNearEnd); 
    const int numStripEndn = rsStripEndn[0].GetNumRows();
    rsStripEndf[0].NewQuery(inContext,Pulser::kFarEnd); 
    const int numStripEndf = rsStripEndf[0].GetNumRows();
    MSG("Pulser",Msg::kInfo)<<"Retrieved "<<numStripEndn<<" Near stripend rows and "<<numStripEndf<< " Far stripend rows"<< endl;
    if (numStripEndn == 0 || numStripEndf == 0 ) {continue;}

    // Change validity range of output to be that of input 
    if(fUseSameValidity) {
      const DbiValidityRec* vrec = rsStripEndn[0].GetValidityRec();
      const VldRange& vrng = vrec->GetVldRange();
      fStartout = vrng.GetTimeStart();
      fEndout = fStartout;
      fEndout.Add(VldTimeStamp(3600*24*62,0));
      range.SetTimeStart(fStartout);
      range.SetTimeEnd(fEndout);
    }

    // Open DbiWriter
    if(fWriteDB) {
      Bool_t iok = fWriternf.Open(range,aggregateno,5,0,fDBName,fDBComment);
      if(!iok) {MSG("Pulser",Msg::kError)<<"Dbiwriter error"<< endl;}
    }

    // Add data from extra runs
    if(fNextra>0) {
      for(int j=0; j<fNextra; j++) {
      DbiSqlContext context(DbiSqlContext::kStarts,fStartExtra[j],fEndExtra[j],fDetector,SimFlag::kData);
      context << " and AGGREGATENO=" << aggregateno;
      rsStripEndn[j+1].NewQuery(context,Pulser::kNearEnd); 
      const int numn = rsStripEndn[j+1].GetNumRows();
      rsStripEndf[j+1].NewQuery(context,Pulser::kFarEnd); 
      const int numf = rsStripEndf[j+1].GetNumRows();
      MSG("Pulser",Msg::kInfo)<<"Retrieved "<<numn<<" Near stripend rows and "<<numf<< " Far stripend rows"<< endl;
      }
    }

    // Arrays for means of far data
    Double_t x[60];
    Double_t er_x[60];

    // Calculate means
    Int_t test = FarMeans(rsStripEndf,x,er_x);
    if(test != 0) {
      MSG("Pulser",Msg::kInfo)<<"Aggregate "<<aggregateno<<" Far Fits failed, ngood = "<<test<<endl;
      continue;
    }

    const PulserGain* rsn[3];

    // Loop over near rows
    for (int ir=0; ir<numStripEndn; ir++) {
      // Pointer to row
      rsn[0] = rsStripEndn[0].GetRow(ir);
      int skey = rsn[0]->GetStripEnd();

      // Pointers to extra rows
      for (int j=1; j<3; j++) {rsn[j] = 0;}
      if(fNextra>0) {
        for (int j=1; j<=fNextra; j++) {
          rsn[j] = rsStripEndn[j].GetRowByIndex(skey);
        }
      }
      // Get points for fit; we plot Near(y) v Far(x)
      Double_t y[60];
      Double_t er_y[60];
      Int_t ntot = 0;
      for (int j=0; j<fNextra+1; j++) {
        for(int i=0; i<rsn[j]->GetNumPoints(); i++) {
          y[ntot+i] = rsn[j]->ZCMean(i);
          er_y[ntot+i] = rsn[j]->ZCError(i);
          // Put floor on errors
          if(er_y[ntot+i]<0.003*y[ntot+i]) {er_y[ntot+i] = 0.003*y[ntot+i];}
        }
        ntot += rsn[j]->GetNumPoints();
      }
      fMinadc = y[0];
      fMaxadc = y[ntot-1];
      // Create and fill graph
      TGraphErrors* g = new TGraphErrors(ntot,x,y,er_x,er_y);
      // Fit Graph
      fFtype = Fit(g);
      // Add row to database
      if(fWriteDB) {
        CalPulserFits row(aggregateno,skey,fFtype,fFit,fMeanRes,fMaxRes,fMaxadc);
        fWriternf << row;
      }
      // Fill ntuple
      if(fNearFarTree!=0) {
        PlexStripEndId seid = PlexStripEndId::UnbuildPlnStripEndKey(fDetector, skey);
        fPlane = seid.GetPlane();
        fStrip = seid.GetStrip();
        PlexPixelSpotId spotid = plex.GetPixelSpotId(seid);
        fBay = spotid.GetRackBay();
        fRack = spotid.GetInRack();
        fTube = spotid.GetTube();
        fPixel = spotid.GetPixel();
        fSpot = spotid.GetSpot();
        PlexLedId ledid = plex.GetLedId(seid);
        fPb = ledid.GetPulserBox();
        fLed = ledid.GetLedInBox();
        // No fit
        if(fFtype < -25 || fFit == 0) {
          fNpar = 0;
          fNpfit = 0;
          fGain = 0;
          fPed = 0;
          fDped = 0;
          fK1 = 0;
          fK2 = 0;
          fK3 = 0;
          fK4 = 0;
          fChisq = 0;
          fDelta200 = 0;
          fDelta7500 = 0;
        }
        // Line fit only
        if(abs(fFtype) == 1) {
          fNpar = 1;
          fNpfit = fFit->GetNumberFitPoints();
          fGain = 1.0/fFit->GetParameter(0);
          fPed = fFit->GetParError(0);
          fDped = fFit->GetParameter(1);
          fK1 = 0;
          fK2 = 0;
          fK3 = 0;
          fK4 = 0;
          fChisq = fFit->GetChisquare();
          fDelta200 = 0;
          fDelta7500 = 0;
        }
        if(abs(fFtype) == 5 || abs(fFtype) == 25) {
          fNpar = fFit->GetNumberFreeParameters()+1;
          fNpfit = fFit->GetNumberFitPoints();
          fGain = fFit->GetParameter(7);
          fPed = fFit->GetParameter(6);
          fDped = fFit->GetParError(6);
          fK1 = fFit->GetParameter(1);
          fK2 = fFit->GetParameter(2);
          if(abs(fFtype) == 25) {
            fK3 = fFit->GetParameter(3);
            fK4 = fFit->GetParameter(4);
          }
          else {
            fK3 = 0;
            fK4 = 0;
          }
          fChisq = fFit->GetChisquare();
          Double_t par[8];
          for (int j=0; j<8; j++) {par[j] = fFit->GetParameter(j);}
          par[0] -= 1.;
          Double_t xtest = 200.;
          Double_t ytest = (CalVaLinearity::FitFunction(&xtest,par)-par[6])/par[7];
          fDelta200 = 100.*(ytest-xtest)/xtest;
          xtest = 7500.;
          ytest = (CalVaLinearity::FitFunction(&xtest,par)-par[6])/par[7];
          fDelta7500 = 100.*(ytest-xtest)/xtest;
        }
        fMeanres = fMeanRes;
        fMaxres = fMaxRes;
        fNearFarTree->Fill();
      }
      delete g;
    }
  }

  if(fWriteDB) {
    if(fWriternf.IsOpen()) {fWriternf.Close();}
  }
}

Int_t PulserGainFit::FarMeans(DbiResultPtr<PulserGain> rsPtrf[3], Double_t* x, Double_t* er_x) {

  // Number of far rows
  Int_t numRowsf = rsPtrf[0].GetNumRows();

  // Pointers to database rows
  const PulserGain* rs[3];

  // Array for strip-end status
  Int_t sgood[1000];
  for(int i=0; i<1000; i++) {sgood[i]=0;}
  Int_t ngood = 0;

  // Use linear fit function with fixed pedestal
  fFitLine->FixParameter(1,0.);

  // Iterate to find good strips
  Int_t nitmax = 3;
  Int_t ngood_last = numRowsf;
  for(int it=0; it<nitmax; it++) {

    // Loop over far rows and calculate means
    ngood = 0;
    for(int i=0; i<60; i++) {
      x[i] = 0.;
      er_x[i] = 0.;
    }

    for (int ir = 0; ir < numRowsf; ++ir) {
      // Get pointer to first row for first run
      rs[0] = rsPtrf[0].GetRow(ir);
      Int_t skey = rs[0]->GetStripEnd();

      // Check this is good stripend
      if(sgood[ir]!=0) {continue;}
      Int_t nmax = rs[0]->GetNumPoints() - 1;
      Float_t adcmax = rs[0]->ZCMean(nmax);
      // Get pointers to extra runs
      for (int j=1; j<3; j++) {rs[j] = 0;}
      if(fNextra > 0) {
        for (int j=1; j<=fNextra; j++) {
          rs[j] = rsPtrf[j].GetRowByIndex(skey);
          nmax = rs[j]->GetNumPoints() - 1;
          if(rs[j]->ZCMean(nmax) > adcmax) {adcmax = rs[j]->ZCMean(nmax);}
        }
      }
      if(adcmax<1000 || adcmax>6000) {continue;}
      ngood++;
      // Now calculate sums for means
      Int_t ntot = 0;
      for (int j=0; j<fNextra+1; j++) {
        for(int i=0; i<rs[j]->GetNumPoints(); i++) {
          Double_t error = rs[j]->ZCError(i);
          x[ntot+i] += rs[j]->ZCMean(i)/(error*error);
          er_x[ntot+i] += 1./(error*error);
        }
        ntot += rs[j]->GetNumPoints();
      }
    }

    // Calculate means
    for(int i=0; i<60; i++) {
      if(er_x[i] > 0.) {
        x[i] = x[i] / er_x[i];
        er_x[i] = 1./sqrt(er_x[i]);
      }
    }

    // Stop here if number of good strips not reducing
    if(ngood_last-ngood < 0.01*ngood_last) {break;}
    ngood_last = ngood;

    // Loop over rows, make plots and fit
    for (int ir = 0; ir < numRowsf; ++ir) {
      // Get pointer to row
      rs[0] = rsPtrf[0].GetRow(ir);
      Int_t skey = rs[0]->GetStripEnd();

      // Get points for fit
      Double_t y[60];
      Double_t er_y[60];
      for(int i=0; i<rs[0]->GetNumPoints(); i++) {
        y[i] = rs[0]->ZCMean(i);
        er_y[i] = rs[0]->ZCError(i);
      }
      Int_t ntot = rs[0]->GetNumPoints();
      if(fNextra > 0) {
        for (int j=1; j<=fNextra; j++) {
          rs[j] = rsPtrf[j].GetRowByIndex(skey);
          for(int i=0; i<rs[j]->GetNumPoints(); i++) {
            y[ntot+i] = rs[j]->ZCMean(i);
            er_y[ntot+i] = rs[j]->ZCError(i);
          }
          ntot += rs[j]->GetNumPoints();
        }
      }

      // Create and fill graph
      TGraphErrors* g = new TGraphErrors(ntot,x,y,er_x,er_y);

      // Fit graph
      Int_t ifit = LineFit(g,fFitLine);

      // Check fit
      if(ifit==0) {
        if(fFitLine->GetNumberFitPoints()<5) {
          sgood[ir] = 1;
        }
        else {
          if(fFitLine->GetChisquare()/(fFitLine->GetNumberFitPoints()-1.)>2.) {
            sgood[ir] = 2;
          }
        }
      }
      else {sgood[ir] = 3;}

      delete g;
    }    
  }
  MSG("Pulser",Msg::kInfo)<<"FarMeans: ngood = "<<ngood<<endl;
  if( ngood > 50) {
    return 0;
  }
  else {
    return ngood;
  }
}

Int_t PulserGainFit::LineFit(TGraphErrors* g, TF1* fFitLine) {

  Int_t n = g->GetN();
  Double_t* x = g->GetX();
  Double_t* y = g->GetY();

  // Set range of fit
  float xmin = fAdcmin;
  float xmax = fAdcmax;
  float ymin = fAdcmin;
  float ymax = fAdcmax;
  int imin = -1;
  int imax = n-1;
  for (int i=0;i<n;i++) {
    if (x[i] > fAdcmin && y[i] > fAdcmin) {
      imin = i;
      if(i>0) {
        xmin = 0.5*(x[i-1]+x[i]);
        ymin = 0.5*(y[i-1]+y[i]);
      }
      break;
    }
  }
  for (int i=0;i<n;i++) {
    if (x[i] > fAdcmax || y[i] > fAdcmax) {
      imax = i-1;
      if (i>0) {
        xmax = 0.5*(x[i-1]+x[i]);
        ymax = 0.5*(y[i-1]+y[i]);
      }
      break;
    }
  }

  // Initialize slope (to prevent a silly value causing problems for
  // subsequent fits
  fFitLine->SetParameter(0,1.);

  // Perform fit: unless we have no point to fit
  if(imin>=0 && imax>imin) {
    g->Fit(fFitLine,"Q"," ",xmin,xmax);
    return 0;
  }
  return 1;
}

// Fit graph to appropriate function and calculate residuals
Int_t PulserGainFit::Fit(TGraphErrors* g)
{
  // Internal variables
  Int_t ftype = -100;
  Int_t npar = 0;
  Int_t npfit = 0;
  Double_t gain = 0.;
  Double_t gain_frac_err = 0.;
  Double_t k1;
  Double_t k2;
  Double_t chisq = 9999.;
  Double_t par[8];
  for (int j=0; j<8; j++) {par[j] = 0.;}
  Float_t sres = 0.;
  Float_t mres = 0.;

  // Number of points
  Int_t n = g->GetN();
  Double_t* x = g->GetX();
  Double_t* er_x = g->GetEX();
  Double_t* y = g->GetY();
  Double_t* er_y = g->GetEY();

  // First do linear fit up to maximum ADC value
  // If this is good we don't do any more fits
  // Set minimum point using same min ADC as pin fits
  Double_t xmin = fAdcmin;
  Double_t ymin = fAdcmin;
  Int_t imin = -1;
  for (int i=0;i<n;i++) {
    if (x[i] > fAdcmin && y[i] > fAdcmin) {
      imin = i;
      if(i>0) {
        xmin = 0.5*(x[i-1]+x[i]);
        ymin = 0.5*(y[i-1]+y[i]);
      }
      break;
    }
  }

  // If no points are above minimum we can't do any fits.
  // In this case give up
  if(imin<0) {return ftype;}

  // Maximum point
  Int_t imax = n-1;
  Double_t xmax = x[imax]+1.;
  Double_t ymax = y[imax]+1.;

  // Perform fit if at least 1 point in range
  fFitLine->SetParameter(0,5.);
  if(imax>=imin) {
    g->Fit(fFitLine,"Q"," ",xmin,xmax);
    fFit = fFitLine;
    ftype = 1;
    npar = 1;
    npfit = fFitLine->GetNumberFitPoints();
    gain = 1.0/fFitLine->GetParameter(0);
    gain_frac_err = fFitLine->GetParError(0)/fFitLine->GetParameter(0);
    chisq = fFitLine->GetChisquare();
    Float_t slope = fFitLine->GetParameter(0);
    sres = 0.;
    mres = 0.;
    for (int i=imin; i<=imax; i++) {
      Float_t res = y[i] - slope*x[i];
      if(slope*x[i]!=0) {res=100.*res/(slope*x[i]);}
      sres+=fabs(res);
      if(fabs(res)>mres) {mres = fabs(res);}
    }
    sres/=(imax-imin+1);
    fMeanRes = sres;
    fMaxRes = mres;
  }
  else {
    // If no points in full range give up
    return ftype;
  }

  // If we have too few points for Pole fit stop here
  if(npfit < 3) {
    ftype = -1;
    return ftype;
  }

  // If full range gives a good linear fit stop here
  if(chisq/(npfit-npar)<2) {return ftype;}

  // If highest point is still in expected linear region stop 
  if(y[n-1] < 7000) {return ftype;}

  // Now do linear fit in expected linear region
  // Set range of fit: minimum is already set
  // We use same max ADC as pin fits
  xmax = fAdcmax;
  ymax = fAdcmax;
  for (int i=0;i<n;i++) {
    if (x[i] > fAdcmax || y[i] > fAdcmax) {
      imax = i-1;
      if (i>0) {
        xmax = 0.5*(x[i-1]+x[i]);
        ymax = 0.5*(y[i-1]+y[i]);
      }
      break;
    }
  }

  // Perform fit if at least 1 point in range
  fFitLine->SetParameter(0,5.);
  if(imax>=imin) {
    g->Fit(fFitLine,"Q"," ",xmin,xmax);
    fFit = fFitLine;
    ftype = 1;
    npar = 1;
    npfit = fFitLine->GetNumberFitPoints();
    gain = 1.0/fFitLine->GetParameter(0);
    gain_frac_err = fFitLine->GetParError(0)/fFitLine->GetParameter(0);
    chisq = fFitLine->GetChisquare();
    Float_t slope = fFitLine->GetParameter(0);
    sres = 0.;
    mres = 0.;
    for (int i=imin; i<=imax; i++) {
      Float_t res = y[i] - slope*x[i];
      if(slope*x[i]!=0) {res=100.*res/(slope*x[i]);}
      sres+=fabs(res);
      if(fabs(res)>mres) {mres = fabs(res);}
    }
    sres/=(imax-imin+1);
    fMeanRes = sres;
    fMaxRes = mres;
  }
  else {
    // If no points in linear region give up
    ftype = -100;
    return ftype;
  }

  // Scale x values with gain for future fits 
  if( gain != 0) {
    for (int i=0;i<n;i++) {
      x[i] /= gain;
      er_x[i] /= gain;
      g->SetPoint(i,x[i],y[i]);
      g->SetPointError(i,er_x[i],er_y[i]);
    }
    xmin /= gain;
  }
  else {
    MSG("Pulser",Msg::kWarning)<<"Zero gain after linear fit"<<endl;
  }

  // Now fit full range with pole only; 
  // Set range of fit: xmin already set, max is last point
  imax = n-1;
  xmax = x[imax]+1.;
  ymax = y[imax]+1.;

  // If too few points give up, saving anything from linear fit
  if(imax-imin<3) {
    ftype = -1;
    return ftype;
  }

  // Set starting parameters
  fFitPole->SetParameter(1,15000.);
  fFitPole->SetParameter(2,5.);
  if(gain != 0.) {
    fFitPole->FixParameter(7,1.0);
  }
  else {
    // This is not necessary if we don't fit channels with no points in linear region, but you never know
    fFitPole->ReleaseParameter(7);
    fFitPole->SetParameter(7,1.0);
  }

  // Perform fit
  g->Fit(fFitPole,"Q"," ",xmin,xmax);

  // Following is now not necessary, but leave anyway
  // If linear fit used same points and pole does not improve chisq sufficiently, store linear fit and give up
  if(npfit==fFitPole->GetNumberFitPoints()) {
    if(chisq/(npfit-npar) - fFitPole->GetChisquare()/(fFitPole->GetNumberFitPoints()-3) < 0.25) {
      return ftype;
    }
  }

  // Set error on slope from linear fit
  fFitPole->SetParError(7,gain_frac_err);

  // Save Pole results
  ftype = 5;
  npar = 3;
  fFit = fFitPole;
  npfit = fFitPole->GetNumberFitPoints();
  gain = fFitPole->GetParameter(7);
  k1 = fFitPole->GetParameter(1); 
  k2 = fFitPole->GetParameter(2);
  chisq = fFitPole->GetChisquare();
  for (int j=0; j<8; j++) {par[j] = fFitPole->GetParameter(j);}
  sres = 0.;
  mres = 0.;
  for (int i=imin; i<=imax; i++) {
    Double_t fy = CalVaLinearity::FitFunction(&x[i],par);
    Float_t res = y[i] - fy;
    if(fy!=0) {res=100.*res/fy;}
    sres+=fabs(res);
    if(fabs(res)>mres) {mres = fabs(res);}
  }
  sres/=(imax-imin+1);
  fMeanRes = sres;
  fMaxRes = mres;
  // Check and flag problem fits
  par[0] -= 1.;
  Double_t xtest = 200.;
  Double_t ytest = CalVaLinearity::FitFunction(&xtest,par)/par[7];
  if(fabs(ytest-xtest)/xtest > 0.01) ftype = -ftype;
  xtest = 7500.;
  ytest = CalVaLinearity::FitFunction(&xtest,par)/par[7];
  if(fabs(ytest-xtest)/xtest > 0.10) ftype = -ftype;

  // If too few points for pole+kicker stop here
  if(imax-imin<5) {return ftype;}

  // Fit pole + kicker with fixed ped
  // Set starting parameters
  if(chisq < 10*(npfit-npar)) {
    fFitPK1->SetParameter(1,k1);
    fFitPK1->SetParameter(2,k2);
  }
  else {
    fFitPK1->SetParameter(1,15000.);
    fFitPK1->SetParameter(2,5.);
  }
  fFitPK1->SetParameter(3,6000.);
  fFitPK1->SetParameter(4,50000.);
  if(gain != 0.) {
    fFitPK1->FixParameter(7,1.0);
  }
  else {
    fFitPK1->ReleaseParameter(7);
    fFitPK1->SetParameter(7,1.0);
  }

  // Perform fit
  g->Fit(fFitPK1,"Q+"," ",xmin,xmax);

  // Set error on slope from linear fit
  fFitPK1->SetParError(7,gain_frac_err);

  // If this an improvement on pole only save values
  if((chisq/(npfit-npar) > fFitPK1->GetChisquare()/(fFitPK1->GetNumberFitPoints()-5)) || ftype < 0) {
    ftype = 25;
    npar = 5;
    fFit = fFitPK1;
    npfit = fFitPK1->GetNumberFitPoints();
    gain = fFitPK1->GetParameter(7); 
    k1 = fFitPK1->GetParameter(1); 
    k2 = fFitPK1->GetParameter(2); 
    chisq = fFitPK1->GetChisquare();
    for (int j=0; j<8; j++) {par[j] = fFitPK1->GetParameter(j);}
    Float_t sres = 0.;
    Float_t mres = 0.;
    for (int i=imin; i<=imax; i++) {
      Double_t fy = CalVaLinearity::FitFunction(&x[i],par);
      Float_t res = y[i] - fy;
      if(fy!=0) {res=100.*res/fy;}
      sres+=fabs(res);
      if(fabs(res)>mres) {mres = fabs(res);}
    }
    sres/=(imax-imin+1);
    fMeanRes = sres;
    fMaxRes = mres;
    // Check and flag problem fits
    par[0] -= 1.;
    Double_t xtest = 200.;
    Double_t ytest = CalVaLinearity::FitFunction(&xtest,par)/par[7];
    if(fabs(ytest-xtest)/xtest > 0.01) ftype = -ftype;
    xtest = 7500.;
    ytest = CalVaLinearity::FitFunction(&xtest,par)/par[7];
    if(fabs(ytest-xtest)/xtest > 0.10) ftype = -ftype;
  }

  // STOP HERE

  /*    
  // If too few points to free ped, stop here
  if(imax-imin<6) {return ftype;}

  // If pole or pole+kicker fit has good chisq stop;
  // otherwise try freeing pedestal
  if(ftype >=5 && chisq/(npfit-npar) < 2.) {return ftype;}

  // Fit pole + kicker with free ped
  // Set starting parameters
  if(chisq < 10*(npfit-npar)) {
    fFitPK2->SetParameter(1,k1);
    fFitPK2->SetParameter(2,k2);
    fFitPK2->SetParameter(7,1.0);
  }
  else {
    fFitPK2->SetParameter(1,15000.);
    fFitPK2->SetParameter(2,5.);
    fFitPK2->SetParameter(7,1.0);
  }
  fFitPK2->SetParameter(3,5000.);
  fFitPK2->SetParameter(4,40000.);

  // Perform fit
  g->Fit(fFitPK2,"Q+"," ",xmin,xmax);

  // If this an improvement save values
  if((chisq/(npfit-npar) > fFitPK2->GetChisquare()/fFitPK2->GetNumberFreeParameters()) || ftype < 0) {
    ftype = 25;
    fFit = fFitPK2;
    for (int j=0; j<8; j++) {par[j] = fFitPK2->GetParameter(j);}
    Float_t sres = 0.;
    Float_t mres = 0.;
    for (int i=imin; i<=imax; i++) {
      Double_t fy = CalVaLinearity::FitFunction(&x[i],par);
      Float_t res = y[i] - fy;
      if(fy!=0) {res=100.*res/fy;}
      sres+=fabs(res);
      if(fabs(res)>mres) {mres = fabs(res);}
    }
    sres/=(imax-imin+1);
    fMeanRes = sres;
    fMaxRes = mres;
    // Check and flag problem fits
    par[0] -= 1.;
    Double_t xtest = 200.;
    Double_t ytest = (CalVaLinearity::FitFunction(&xtest,par)-par[6])/par[7];
    if(fabs(ytest-xtest)/xtest > 0.01) ftype = -ftype;
    xtest = 7500.;
    ytest = CalVaLinearity::FitFunction(&xtest,par)/par[7];
    if(fabs(ytest-xtest)/xtest > 0.10) ftype = -ftype;
  }
  */

  return ftype;
}

---