LITuning.cxx

Go to the documentation of this file.

// Program name: LITuning.cxx                
//
// Package: LISummary          
//
// Coded by Jeff Hartnell Feb/2003                             
//
// Purpose: To provide algorithms for the tuning of the LI system
//
// Contact: jeffrey.hartnell@physics.ox.ac.uk                         

#include "TString.h"

#include "MessageService/MsgService.h"

#include "LISummary/LITuning.h"
#include "LISummary/LIRun.h"
#include <cmath>

ClassImp(LITuning)
  
  CVSID("$Id: LITuning.cxx,v 1.27 2007/11/11 08:00:53 rhatcher Exp $");

//......................................................................

LITuning::LITuning()
{
  MSG("LITuning", Msg::kDebug) 
    <<"Running LITuning constructor..."<<endl;

  this->InitialiseDataMembers(Detector::kUnknown);

  //some good runs to test the code on:
  //15853 - has lots of really low points
  //17681 - a drift point run
  //17661 - has a very saturated led and a bad pin
  //      - also has some missing info at low ph; 0 adc and 0 ph
  //17616 - led 1 has two bad pins and doesn't saturate!!!!!

  MSG("LITuning", Msg::kDebug) 
    <<"Constructed LITuning object"<<endl;
}

//......................................................................

LITuning::~LITuning()
{
  MSG("LITuning", Msg::kDebug) 
    <<"Running LITuning destructor..."<<endl;

  MSG("LITuning", Msg::kDebug) 
    <<"Finished LITuning destructor"<<endl;
}

//......................................................................

void LITuning::InitialiseDataMembers(Detector::Detector_t det)
{
  MSG("LITuning", Msg::kDebug) 
    <<"Running InitialiseDataMembers method..."<<endl;

  //initialise data members
 
  if (det==Detector::kFar){
    fAdcAtStartOfSaturation=8000;
    fdADCdPHAtSat=1;
    fDetector=det;
    fFirstGcPointAdc=350;
    fGainCurve=0;//set pointer to zero
    fIdealAdc=3500;
    fLastGcPointAdc=13800;
    fPinHighestMaxHG=8000;
    fPinHighestMaxLG=8000;
    fPinLowestMaxHG=400;
    fPinLowestMaxLG=200;
    fPinLowestMinHG=80;
    fPinLowestMinLG=80;
    fMaxPossAdc=15000;
    fMaxPossPh=1023;
    fMaxPossPin=11000;
    fMinHGdPINdPH=3;
    fMinLGdPINdPH=0.5;
    fMinPossAdc=0;
    fMinPossPh=0;
    fMinPossPin=0;
    fNumCalibPoints=15;
    fNumGainPoints=20;
    fS="";
  }
  else if (det==Detector::kNear){
    //these have not been set properly
    //jjh added best guesses on 19/Dec/04
    fAdcAtStartOfSaturation=40000;
    fdADCdPHAtSat=1;
    fDetector=det;
    fFirstGcPointAdc=350;
    fGainCurve=0;//set pointer to zero
    fIdealAdc=5000;
    fLastGcPointAdc=30000;
    fPinHighestMaxHG=60000;
    fPinHighestMaxLG=60000;
    fPinLowestMaxHG=100;
    fPinLowestMaxLG=100;
    fPinLowestMinHG=70;
    fPinLowestMinLG=70;
    fMaxPossAdc=60000;
    fMaxPossPh=1023;
    fMaxPossPin=60000;
    fMinHGdPINdPH=3;
    fMinLGdPINdPH=0.5;
    fMinPossAdc=0;
    fMinPossPh=0;
    fMinPossPin=0;
    fNumCalibPoints=15;
    fNumGainPoints=20;
    fS="";
  }
  else if (det==Detector::kUnknown){//make default same as FD
    fAdcAtStartOfSaturation=8000;
    fdADCdPHAtSat=1;
    fDetector=Detector::kFar;
    fFirstGcPointAdc=350;
    fGainCurve=0;//set pointer to zero
    fIdealAdc=3500;
    fLastGcPointAdc=13800;
    fPinHighestMaxHG=8000;
    fPinHighestMaxLG=8000;
    fPinLowestMaxHG=400;
    fPinLowestMaxLG=200;
    fPinLowestMinHG=80;
    fPinLowestMinLG=80;
    fMaxPossAdc=15000;
    fMaxPossPh=1023;
    fMaxPossPin=11000;
    fMinHGdPINdPH=3;
    fMinLGdPINdPH=0.5;
    fMinPossAdc=0;
    fMinPossPh=0;
    fMinPossPin=0;
    fNumCalibPoints=15;
    fNumGainPoints=20;
    fS="";
  }
  else {
    MSG("LITuning", Msg::kFatal) 
      <<"Detector Type "<<Detector::AsString(fDetector)
      <<" not implemented"<<endl
      <<"Exiting now..."<<endl<<endl;
    exit(0);
  }

  MSG("LITuning", Msg::kDebug) 
    <<"Finished InitialiseDataMembers method"<<endl;
}

//......................................................................

void LITuning::InputDataGc(vector<LIRun>& gc) 
{
  this->InitialiseDataMembers(gc.begin()->GetDetector());

  //set data member according to input
  fGainCurve=&gc;//get the address of the vector
}

//......................................................................

void LITuning::SetNumCalibPoints(Int_t i) 
{
  if (i<=1){
    MSG("LITuning", Msg::kFatal) 
      <<endl
      <<" *** Warning: Number of calib points <= 1 *** "<<endl
      <<"Data input is most likely a drift point run!"<<endl
      <<"Exiting now..."<<endl<<endl;
    exit(0);
  }

  fNumCalibPoints=i;
}

//......................................................................

void LITuning::SetNumGainPoints(Int_t i)
{
  if (i<=1){
    MSG("LITuning", Msg::kFatal) 
      <<endl
      <<" *** Warning: Number of gain curve points <= 1 *** "<<endl
      <<"Exiting now..."<<endl<<endl;
    exit(0);
  }

  fNumGainPoints=i;
}

//......................................................................

void LITuning::PrintConfig()
{
  MSG("LITuning", Msg::kInfo) 
    <<endl
    <<"LITuning object has the following configuration:"<<endl
    <<"fDetector="<<Detector::AsString(fDetector)<<endl
    <<"fNumGainPoints="<<fNumGainPoints<<endl
    <<"fIdealAdc="<<fIdealAdc<<endl
    <<"fFirstGcPointAdc="<<fFirstGcPointAdc
    <<", fLastGcPointAdc="<<fLastGcPointAdc<<endl
    <<"fdADCdPHAtSat="<<fdADCdPHAtSat<<endl
    <<"fAdcAtStartOfSaturation="<<fAdcAtStartOfSaturation<<endl
    <<"fNumCalibPoints="<<fNumCalibPoints<<endl
    <<"fPinHighestMaxLG="<<fPinHighestMaxLG
    <<", fPinHighestMaxHG="<<fPinHighestMaxHG<<endl
    <<"fPinLowestMaxLG="<<fPinLowestMaxLG
    <<", fPinLowestMaxHG="<<fPinLowestMaxHG<<endl
    <<"fPinLowestMinLG="<<fPinLowestMinLG
    <<", fPinLowestMinHG="<<fPinLowestMinHG<<endl
    <<"fMinPossAdc="<<fMinPossAdc
    <<", fMaxPossAdc="<<fMaxPossAdc<<endl
    <<"fMinPossPh="<<fMinPossPh
    <<", fMaxPossPh="<<fMaxPossPh<<endl
    <<"fMinPossPin="<<fMinPossPin
    <<", fMaxPossPin="<<fMaxPossPin<<endl
    <<"fMinLGdPINdPH="<<fMinLGdPINdPH
    <<", fMinHGdPINdPH="<<fMinHGdPINdPH<<endl
    <<endl;
}

//......................................................................

void LITuning::PrintPhPinMsg(vector<Double_t>::iterator adc,
                             vector<Double_t>::iterator ph,
                             vector<Double_t>::iterator pin,
                             Double_t interpPh,
                             Double_t interpPin)
{
  MSG("LITuning",Msg::kInfo)
    <<"    Lower point: adc="<<static_cast<Int_t>(*(adc-1))
    <<", PH="<<static_cast<Int_t>(*(ph-1))
    <<endl<<"    Upper point: adc="<<static_cast<Int_t>(*adc)
    <<", PH="<<static_cast<Int_t>(*ph)
    <<", Interpolated PH="<<static_cast<Int_t>(interpPh)
    <<endl<<"In PIN space:"
    <<endl<<"    Lower point: pinAdc="<<static_cast<Int_t>(*(pin-1))
    <<", PH="<<static_cast<Int_t>(*(ph-1))
    <<endl<<"    Upper point: pinAdc="<<static_cast<Int_t>(*pin)
    <<", PH="<<static_cast<Int_t>(*ph)
    <<", Interpolated Pin="<<static_cast<Int_t>(interpPin)<<endl;
}

//......................................................................

void LITuning::PrintPhMsg(vector<Double_t>::iterator adc,
                          vector<Double_t>::iterator ph,
                          Double_t interpPh)
{
  MSG("LITuning",Msg::kInfo)
    <<"    Lower point: adc="<<static_cast<Int_t>(*(adc-1))
    <<", PH="<<static_cast<Int_t>(*(ph-1))
    <<endl<<"    Upper point: adc="<<static_cast<Int_t>(*adc)
    <<", PH="<<static_cast<Int_t>(*ph)
    <<", Interpolated PH="<<static_cast<Int_t>(interpPh)<<endl;
}

//......................................................................

void LITuning::PrintPinMsg(vector<Double_t>::iterator pin,
                           vector<Double_t>::iterator ph,
                           Double_t interpAdc,
                           Double_t interpPh)
{
  MSG("LITuning",Msg::kInfo)
    <<"    Lower point: pin="<<static_cast<Int_t>(*(pin-1))
    <<", PH="<<static_cast<Int_t>(*(ph-1))
    <<endl<<"    Upper point: pin="<<static_cast<Int_t>(*pin)
    <<", PH="<<static_cast<Int_t>(*ph)
    <<", Interpolated PH="<<static_cast<Int_t>(interpPh)<<endl
    <<"    Corresponding adc="<<interpAdc<<endl; 
}

//......................................................................

void LITuning::PrintFirstLastMsg(Double_t firstPh,Double_t lastPh,
                                 Double_t firstPin,Double_t lastPin)
{
    MSG("LITuning_TuneGc",Msg::kInfo) 
      <<"Input calculated first/last GC data points are:"<<endl
      <<static_cast<Int_t>(firstPh)<<" < PH < "
      <<static_cast<Int_t>(lastPh)<<endl
      <<static_cast<Int_t>(firstPin)<<" < PIN < "
      <<static_cast<Int_t>(lastPin)<<endl
      <<"Number of required GC points = "<<fNumGainPoints<<endl;
}


//......................................................................

Float_t LITuning::FindLowerFraction(map<Float_t,Int_t>& inputMap,
                                      Float_t fraction,
                                      Int_t totalNumHits)
{ 
  MSG("LITuning",Msg::kDebug)  
    <<endl<<" ** Running FindLowerFraction method... ** "<<endl; 

  if (totalNumHits<=0){
    MSG("LITuning",Msg::kWarning)
      <<"Can't find lower fraction: Total number of hits is zero"<<endl;
    MSG("LITuning",Msg::kDebug)
      <<endl<<" ** Finished FindLowerFraction method ** "<<endl; 
    return -1;
  }

  //loop over all the hits for a calib point
  map<Float_t,Int_t>::iterator hit=inputMap.begin();

  Int_t totalHitsIteratedOver=0;
  
  //find the lower mean at X% of hits
  while (hit!=inputMap.end()){
    
    totalHitsIteratedOver+=hit->second;
    
    MSG("LITuning",Msg::kVerbose)
      <<"mean="<<hit->first
      <<", num="<<hit->second
      <<", runTotal="<<totalHitsIteratedOver
      <<endl;
    
    if (totalHitsIteratedOver>totalNumHits*fraction){
      MSG("LITuning",Msg::kDebug)
        <<"Found lower bound, mean="<<hit->first
        <<", num="<<totalHitsIteratedOver<<"/"<<totalNumHits
        <<" (="<<100.*totalHitsIteratedOver/totalNumHits
        <<"%)"<<endl;

      MSG("LITuning",Msg::kDebug)
        <<endl<<" ** Finished FindLowerFraction method ** "<<endl; 
      return hit->first;
    }
    
    hit++;
  }
  
  MSG("LITuning",Msg::kDebug)
    <<endl<<" ** Finished FindLowerFraction method ** "<<endl; 
  return -1;
}

//......................................................................

Float_t LITuning::FindUpperFraction(map<Float_t,Int_t>& inputMap,
                                      Float_t fraction,
                                      Int_t totalNumHits)
{ 
  MSG("LITuning",Msg::kDebug)  
    <<endl<<" ** Running FindUpperFraction method... ** "<<endl; 

  if (totalNumHits<=0){
    MSG("LITuning",Msg::kWarning)
      <<"Can't find upper fraction: Total number of hits is zero"<<endl;
    MSG("LITuning",Msg::kDebug)  
      <<endl<<" ** Finished FindUpperFraction method... ** "<<endl; 
    return -1;
  }

  //loop over all the hits for a calib point
  map<Float_t,Int_t>::iterator hit=inputMap.end();
  hit--;

  Int_t totalHitsIteratedOver=0;
  
  //find the upper mean at X% of hits
  while (hit!=inputMap.begin()){
    
    totalHitsIteratedOver+=hit->second;
    
    MSG("LITuning",Msg::kVerbose)
      <<"mean="<<hit->first
      <<", num="<<hit->second
      <<", runTotal="<<totalHitsIteratedOver
      <<endl;
    
    if (totalHitsIteratedOver>totalNumHits*fraction){
      MSG("LITuning",Msg::kDebug)
        <<"Found upper bound, mean="<<hit->first
        <<", num="<<totalHitsIteratedOver<<"/"<<totalNumHits
        <<" (="<<100.*totalHitsIteratedOver/totalNumHits
        <<"%)"<<endl;

      MSG("LITuning",Msg::kDebug)
        <<endl<<" ** Finished FindUpperFraction method ** "<<endl; 
      return hit->first;
    }
    
    hit--;
  }
  
  MSG("LITuning",Msg::kDebug)
    <<endl<<" ** Finished FindUpperFraction method ** "<<endl; 
  return -1;
}

//......................................................................

Float_t LITuning::FindMidFractionAv(map<Float_t,Int_t>& inputMap,
                                      Float_t lowerFraction,
                                      Float_t upperFraction,
                                      Int_t totalNumHits)
{ 

  MSG("LITuning",Msg::kDebug)  
    <<endl<<" ** Running FindMidFractionAv method... ** "<<endl; 
  
  if (totalNumHits<=0){
    MSG("LITuning",Msg::kWarning)
      <<"Can't find mid fraction average: Total number of hits is zero"
      <<endl;
    MSG("LITuning",Msg::kDebug)  
      <<endl<<" ** Finished FindMidFractionAv method... ** "<<endl; 
    return -1;
  }
  
  map<Float_t,Int_t>::iterator hit=inputMap.begin();
  Int_t totalHitsIteratedOver=0;
  Float_t midFractionAv=0;  
  Int_t numMidFractionHits=0;

  //find the upper mean at X% of hits
  while (hit!=inputMap.end()){
    
    totalHitsIteratedOver+=hit->second;
    
    //ignore the lower X% and the upper Y%
    if (totalHitsIteratedOver<totalNumHits*lowerFraction ||
        totalHitsIteratedOver>totalNumHits*(1-upperFraction)){

      MSG("LITuning",Msg::kVerbose)
        <<"mean="<<hit->first
        <<", num="<<hit->second
        <<", runTotal="<<totalHitsIteratedOver
        <<endl;
      hit++;
      continue;
    }

    //add up the total mean (there might be more than one)
    midFractionAv+=hit->second*hit->first;
    //count the number of hits added
    numMidFractionHits+=hit->second;

    MSG("LITuning",Msg::kVerbose)
        <<"Using hit, mean="<<hit->first
        <<", num="<<totalHitsIteratedOver<<"/"<<totalNumHits
        <<" ("<<100.*totalHitsIteratedOver/totalNumHits
        <<"%)"
        <<", run total="<<midFractionAv
        <<", num added="<<numMidFractionHits
        <<endl;

    hit++;
  }
  
  //calculate the average
  if (numMidFractionHits>0){
    midFractionAv/=numMidFractionHits;
    MSG("LITuning",Msg::kDebug)
      <<"Mid fraction average="<<midFractionAv<<endl;
  }
  else if (totalNumHits==1){
    MSG("LITuning",Msg::kDebug)
      <<"Only 1 hit for this led point"
      <<", returning mid fraction average of zero"<<endl;
  }
  else{
    MSG("LITuning",Msg::kWarning)
      <<endl
      <<"Avoided FPE: Number of hits in mid fraction average is zero"
      <<", returning zero mean"
      <<endl<<"Total number of hits for this led is "<<totalNumHits
      <<endl;
  }

  MSG("LITuning",Msg::kDebug)
    <<endl<<" ** Finished FindMidFractionAv method ** "<<endl; 
  return midFractionAv;
}

//......................................................................

Double_t LITuning::CalcAv_dPINdPH(vector<LIRun>::iterator gcData,
                                  Int_t pinNumber)
{
  MSG("LITuning", Msg::kDebug) 
    <<"Running CalcAv_dPINdPH..."<<endl;

  //get iterators to gc components
  vector<Double_t> ph;
  ph=gcData->GetPh();
  vector<Double_t> pin;
  pin=gcData->GetPin(pinNumber);
  
  Double_t av_dPINdPH=0;
  Int_t counter=0;

  Double_t pinMin=-1;
  Double_t pinMax=-1;

  if (pinNumber==1){
    pinMin=fPinLowestMinHG;
    pinMax=fPinHighestMaxHG;
  }
  else if (pinNumber==2){
    pinMin=fPinLowestMinLG;
    pinMax=fPinHighestMaxLG;
  }

  //loop and add up gradients    
  for (UInt_t i=0;i<ph.size()-1;i++){
    //don't include points above max
    //don't include points with zero pulse height
    if ((ph[i+1]-ph[i])!=0 && pin[i]<pinMax && ph[i]>0){
      av_dPINdPH+=(pin[i+1]-pin[i])/(ph[i+1]-ph[i]);
      counter++;
    }
  }
  
  //calc average
  if (counter!=0) av_dPINdPH/=counter;

  MSG("LITuning",Msg::kDebug) 
    <<"CalcAv_dPINdPH method finished"<<endl;
  return av_dPINdPH;
}

//......................................................................

Double_t LITuning::CalcLastGcPointAdc(vector<LIRun>::iterator gcData)
{
  MSG("LITuning", Msg::kDebug) 
    <<"Running CalcLastGcPointAdc method..."<<endl;

  //get iterators to gc components
  vector<Double_t> ph;
  ph=gcData->GetPh();
  vector<Double_t> adc;
  adc=gcData->GetAdc();
  
  Int_t led=gcData->GetLed();
  Int_t pb=gcData->GetPb();

  Double_t dADCdPH=-1;
  Double_t adcAfterSat=-1;
  UInt_t pointAfterSat=0;
  Double_t dADCdPHAfterSat=-1;
  
  Double_t adcAtSat=1;
  Bool_t saturatedPointFound=false;

  //reset the general purpose string
  fS="";

  for (UInt_t i=0;i<ph.size()-1;i++){
    //avoid fpe and only check for points up near saturation
    if ((ph[i+1]-ph[i])!=0 && adc[i]>fAdcAtStartOfSaturation){

      //calc grad of adc vs ph
      dADCdPH=(adc[i+1]-adc[i])/(ph[i+1]-ph[i]);
      string s=Form("%f",static_cast<Float_t>(dADCdPH));
      fS+=s+",";
      if (dADCdPH<fdADCdPHAtSat){
        if (MsgService::Instance()->IsActive("LITuning_LastGc",
                                             Msg::kDebug)){
          MSG("LITuning",Msg::kInfo)
            <<"dADCdPH is bad ("<<dADCdPH<<") for point "<<i
            <<", adc="<<adc[i]<<", ph="<<ph[i]<<endl;
        }
        if (!saturatedPointFound) {
          adcAfterSat=adc[i];
          pointAfterSat=i;
          dADCdPHAfterSat=dADCdPH;
        }
        saturatedPointFound=true;
      }
    }
  }

  if (MsgService::Instance()->IsActive("LITuning_LastGc",
                                       Msg::kDebug)){
    MSG("LITuning",Msg::kInfo) 
      <<"Gradients are: "<<fS<<endl;
  }

  if (!saturatedPointFound){
    MSG("LITuning",Msg::kWarning) 
      <<endl
      <<"("<<pb<<":"<<led<<")"
      <<" Saturation point never reached in data gain curve"
      <<", last dADCdPH="<<dADCdPH<<endl;
  }
  else{
    if (MsgService::Instance()->IsActive("LITuning_LastGc",
                                         Msg::kDebug)){
      MSG("LITuning",Msg::kInfo)
        <<"("<<pb<<":"<<led<<")"
        <<" adcAfterSat="<<static_cast<Int_t>(adcAfterSat)
        <<", fLastGcPointAdc"<<static_cast<Int_t>(fLastGcPointAdc)
        <<", diff="<<static_cast<Int_t>(fLastGcPointAdc-adcAfterSat)
        <<endl
        <<"       pointAfterSat="<<pointAfterSat
        <<", dADCdPHAfterSat="<<dADCdPHAfterSat<<endl;
    }

    Int_t i=pointAfterSat;
    
    if (pointAfterSat>0){
      if ((ph[i]-ph[i-1])!=0 && (ph[i+1]-ph[i])!=0){
        //interpolate (y-y2)*m+x2;
        //to find adc corresponding to fdADCdPHAtSat
        Double_t dADCdPH1=(adc[i]-adc[i-1])/(ph[i]-ph[i-1]);
        Double_t dADCdPH2=(adc[i+1]-adc[i])/(ph[i+1]-ph[i]);
        //perhaps I should average the two adcs!!! rather than extremes
        Double_t m=(adc[i+1]-adc[i-1])/(dADCdPH2-dADCdPH1);
        adcAtSat=this->Interpolate(fdADCdPHAtSat,dADCdPH2,m,adc[i+1]);

        if (MsgService::Instance()->IsActive("LITuning_LastGc",
                                             Msg::kDebug)){
          MSG("LITuning",Msg::kInfo)
            <<"dADCdPH1="<<dADCdPH1<<", dADCdPH2="<<dADCdPH2
            <<",   m="<<m<<", adcAtSat="<<adcAtSat<<endl
            <<"adc[i-1]="<<static_cast<Int_t>(adc[i-1])
            <<",   adc[i+1]="<<static_cast<Int_t>(adc[i+1])<<endl
            <<" ph[i-1]="<<ph[i-1]<<",      ph[i+1]="<<ph[i+1]
            <<endl;
        }
      }
    }
  }
  
  MSG("LITuning",Msg::kDebug) 
    <<"CalcLastGcPointAdc method finished"<<endl;
  return adcAtSat;
}

//......................................................................

void LITuning::CalcFirstGcPoint(vector<LIRun>::iterator gcData,
                                Int_t goodPin,Double_t& firstPh,
                                Double_t& firstPin)
{
  //need to find where led switches on
  //this method sets firstPin and firstPh

  MSG("LITuning", Msg::kDebug) 
    <<"Running CalcFirstGcPoint..."<<endl;

  Bool_t pointFound=false;
  Bool_t firstPointTooHigh=false;

  //get iterators to gc components
  vector<Double_t> phVector;
  vector<Double_t>::iterator ph;
  phVector=gcData->GetPh();
  ph=phVector.begin();
  vector<Double_t> adcVector;
  vector<Double_t>::iterator adc;
  adcVector=gcData->GetAdc();
  adc=adcVector.begin();
  vector<Double_t>pinVector;
  vector<Double_t>::iterator pin;
  pinVector=gcData->GetPin(goodPin);
  pin=pinVector.begin();
  
  Int_t led=gcData->GetLed();
  Int_t pb=gcData->GetPb();

  MSG("LITuning_FirstGc", Msg::kDebug) 
    <<"("<<pb<<":"<<led<<")"
    <<"Looking for lowest gain curve point..."<<endl;

  Int_t counter=0;//already done one point
  //loop over the points
  while(ph!=phVector.end()){
    //counter
    counter++;

    //check if the first data point has an adc value greater than the
    //first required point
    if (*adc>fFirstGcPointAdc && adc==adcVector.begin()){
        
      //set control variables
      pointFound=false;
      firstPointTooHigh=true;
        
      //jump out of the loop
      break;
    }
    
    //check if the current point has an adc value greater than the
    //first required point
    //and whether a point has been found yet
    if (*adc>fFirstGcPointAdc && !pointFound){

      //set control variable
      pointFound=true;

      //it's best to use two known points than one effectively
      //unknown point
      if (*(adc-1)==0 || *(adc-1)==1){//==1 is a hack, remove eventually
        if ((adc+1)!=phVector.end()){
          MSG("LITuning_FirstGc", Msg::kInfo) 
            <<"("<<pb<<":"<<led<<")"
            <<" Wanted to interpolate using a zero adc point"<<endl
            <<"       but using two higher points instead"<<endl;

          //interpolate in pmt adc space (y-y2)*m+x2;
          Double_t m=(*(ph+1)-*ph)/(*(adc+1)-*adc);
          firstPh=this->Interpolate(fFirstGcPointAdc,*(adc+1),m,
                                    *(ph+1));
          
          //work out the pin value that corresponds to the led
          //just having switched on
          //calculate gradient in pin adc space
          m=(*(pin+1)-*pin)/(*(ph+1)-*ph);
          firstPin=this->Interpolate(firstPh,*(ph+1),m,*(pin+1));

          //correct the values if necessary       
          this->CorrectInterpValues(adc+1,ph+1,pin+1,led,pb,
                                    fFirstGcPointAdc,firstPh,firstPin);
          //jump out of the loop
          break;
        }
        else{
          //else just continue and use the zero point - rare case
          MSG("LITuning_FirstGc", Msg::kInfo) 
            <<" ("<<pb<<":"<<led<<")"
            <<"Wanted to avoid interpolating using a zero adc point"
            <<endl
            <<", but no higher point"<<endl;
        }
      }

      //interpolate in pmt adc space (y-y2)*m+x2;
      Double_t m=(*ph-*(ph-1))/(*adc-*(adc-1));
      firstPh=this->Interpolate(fFirstGcPointAdc,*adc,m,*ph);
      
      //calculate gradient in pin adc space
      m=(*pin-*(pin-1))/(*ph-*(ph-1));
      firstPin=this->Interpolate(firstPh,*ph,m,*pin);

      if (MsgService::Instance()->IsActive("LITuning_FirstGc",
                                           Msg::kDebug)){
        MSG("LITuning_FirstGc",Msg::kInfo)
          <<"("<<pb<<":"<<led<<")"
          <<" Lowest GC point found"
          <<" ("<<counter<<"/"<<adcVector.size()<<")"
          <<", Required Point="<<fFirstGcPointAdc<<endl;
        this->PrintPhPinMsg(adc,ph,pin,firstPh,firstPin);
      }

      //jump out of the loop when point is found
      break;
    }
    
    //go to the next point
    ph++;
    adc++;
    pin++; 
  }

  //check on success of finding point
  if (firstPointTooHigh){
    this->InterpolateAboveOrBelow(kGcUnknown,
                                  adcVector.begin()+1,
                                  phVector.begin()+1,
                                  pinVector.begin()+1,led,pb,
                                  fFirstGcPointAdc,firstPh,firstPin);
    //set control variable
    pointFound=true;
  }
  else if (!pointFound){//i.e. first data point was too low
    this->InterpolateAboveOrBelow(kGcUnknown,
                                  adcVector.end()-1,phVector.end()-1,
                                  pinVector.end()-1,led,pb,
                                  fFirstGcPointAdc,firstPh,firstPin);
    //set control variable
    pointFound=true;
  }

  //this is a bit of a hack
  if (firstPin==0){
    MSG("LITuning_FirstGc",Msg::kWarning)
      <<endl
      <<"("<<pb<<":"<<led<<")"
      <<" Looking for first GC point but first PIN value has"
      <<" been set to zero"<<endl;
    
    ph=phVector.begin()+1;
    adc=adcVector.begin()+1;
    pin=pinVector.begin()+1;
    
    Bool_t pinFound=false;
    Bool_t lookForSecondPinNow=false;

    Int_t counter=0;//already done one point
    //loop over the points
    while(pin!=pinVector.end()){
      //counter
      counter++;

      MSG("LITuning_FirstGc",Msg::kVerbose)
        <<"("<<pb<<":"<<led<<")"<<" Looping...   counter="
        <<counter<<endl;
      
      //check if the current point has an adc value greater than the
      //first required point
      //and whether a point has been found yet
      if (*pin>fMinPossPin && lookForSecondPinNow && !pinFound){

        if (MsgService::Instance()->IsActive("LITuning_FirstGc",
                                             Msg::kDebug)){     
          MSG("LITuning_FirstGc",Msg::kInfo)
            <<"("<<pb<<":"<<led<<")"
            <<" Found 2nd non-zero PIN"<<endl;
        }

        pinFound=true;

        //recalculate gradient in pin adc space
        Double_t m=0;
        if (*ph-*(ph-1)!=0) m=(*pin-*(pin-1))/(*ph-*(ph-1));
        firstPin=this->Interpolate(firstPh,*ph,m,*pin); 

        if (MsgService::Instance()->IsActive("LITuning_FirstGc",
                                             Msg::kDebug)){   
          MSG("LITuning_FirstGc",Msg::kInfo)
            <<"("<<pb<<":"<<led<<")"<<" Interpolated ok, counter="
            <<counter<<endl;
          this->PrintPhPinMsg(adc,ph,pin,firstPh,firstPin);
        }

        //jump out of loop
        break;
      }

      if (*pin>fMinPossPin){
        lookForSecondPinNow=true;
        if (MsgService::Instance()->IsActive("LITuning_FirstGc",
                                             Msg::kDebug)){
          MSG("LITuning_FirstGc",Msg::kInfo)
            <<"("<<pb<<":"<<led<<")"
            <<" Found first non-zero pin, looking for second..."<<endl;
        }
      }
      pin++;
      adc++;
      ph++;
    }

    if (MsgService::Instance()->IsActive("LITuning_FirstGc",
                                         Msg::kDebug)){
      MSG("LITuning_FirstGc",Msg::kInfo)
        <<endl
        <<"("<<pb<<":"<<led<<")"
        <<" After recalculating pin value:"
        <<endl;
      this->PrintPhPinMsg(adc,ph,pin,firstPh,firstPin);
    }

    //correct the values if needed
    this->CorrectInterpValues(adc,ph,pin,led,pb,fFirstGcPointAdc,
                              firstPh,firstPin);
  }

  MSG("LITuning",Msg::kDebug) 
    <<"CalcFirstGcPoint method finished"<<endl;
}

//......................................................................

void LITuning::CalcLastGcPoint(vector<LIRun>::iterator gcData,
                               Int_t goodPin,Double_t& lastPh,
                               Double_t& lastPin)
{
  //need to find where led saturates pmt
  //this method sets lastPin and lastPh

  MSG("LITuning", Msg::kVerbose) 
    <<"Running CalcLastGcPoint..."<<endl;

  Bool_t pointFound=false;
  Bool_t firstPointTooHigh=false;

  //get iterators to gc components
  vector<Double_t> phVector;
  vector<Double_t>::iterator ph;
  phVector=gcData->GetPh();
  ph=phVector.begin();
  vector<Double_t> adcVector;
  vector<Double_t>::iterator adc;
  adcVector=gcData->GetAdc();
  adc=adcVector.begin();
  vector<Double_t>pinVector;
  vector<Double_t>::iterator pin;
  pinVector=gcData->GetPin(goodPin);
  pin=pinVector.begin();
  
  Int_t led=gcData->GetLed();
  Int_t pb=gcData->GetPb();

  MSG("LITuning_LastGc", Msg::kDebug) 
    <<" ("<<pb<<":"<<led<<")"
    <<"Looking for highest gain curve point..."<<endl;

  //declare the variable used to find the last/highest GC point
  //initialise to default
  Double_t lastGcPointAdc=fLastGcPointAdc;

  //only try and calculate last point adc if far detector
  if (fDetector==-1){
    Double_t tempAdc=this->CalcLastGcPointAdc(gcData);
    
    if (fLastGcPointAdc-lastGcPointAdc>4000 ||
        fLastGcPointAdc-lastGcPointAdc<-4000){
      MSG("LITuning_LastGc",Msg::kWarning)
        <<endl
        <<"("<<pb<<":"<<led<<")"
        <<" Calculated lastGcPointAdc="
        <<static_cast<Int_t>(tempAdc)
        <<", fLastGcPointAdc="<<static_cast<Int_t>(fLastGcPointAdc)<<endl
        <<"       Difference="
        <<static_cast<Int_t>(fLastGcPointAdc-tempAdc)
        <<". Note: Using fLastGcPointAdc"<<endl;
      lastGcPointAdc=tempAdc;
    }
  }
  else{
    //just use the one-size-fits-all maximum fLastGcPointAdc
  }

  Int_t counter=0;
  //loop over the points
  while(ph!=phVector.end()){
    //counter
    counter++;

    //check if the first data point has an adc value greater than the
    //last required point
    if (*adc>lastGcPointAdc && adc==adcVector.begin()){
        
      //set control variables
      pointFound=false;
      firstPointTooHigh=true;
        
      //jump out of the loop
      break;
    }
    
    //check if the current point has an adc value greater than the
    //last required point
    //and whether a point has been found yet
    if (*adc>lastGcPointAdc && !pointFound){
        
      //set control variable
      pointFound=true;

      //interpolate in pmt adc space (y-y2)*m+x2;
      Double_t m=(*ph-*(ph-1))/(*adc-*(adc-1));
      lastPh=this->Interpolate(lastGcPointAdc,*adc,m,*ph);
      
      //calculate gradient in pin adc space
      m=(*pin-*(pin-1))/(*ph-*(ph-1));
      lastPin=this->Interpolate(lastPh,*ph,m,*pin);

      if (MsgService::Instance()->IsActive("LITuning_LastGc",
                                           Msg::kDebug)){
        MSG("LITuning_LastGc",Msg::kInfo)
          <<"("<<pb<<":"<<led<<")"
          <<" Highest GC point found"
          <<" ("<<counter<<"/"<<adcVector.size()<<")"
          <<", Required Point="<<lastGcPointAdc<<endl;
        this->PrintPhPinMsg(adc,ph,pin,lastPh,lastPin);
      }

      //jump out of the loop when point is found
      break;
    }
    
    //go to the next point
    ph++;
    adc++;
    pin++; 
  }

  //check on success of finding point
  if (firstPointTooHigh){
    this->InterpolateAboveOrBelow(kGcUnknown,
                                  adcVector.begin()+1,
                                  phVector.begin()+1,
                                  pinVector.begin()+1,led,pb,
                                  lastGcPointAdc,lastPh,lastPin);
    //set control variable
    pointFound=true;
  }
  else if (!pointFound){//i.e. last data point was too low
    this->InterpolateAboveOrBelow(kGcUnknown,
                                  adcVector.end()-1,phVector.end()-1,
                                  pinVector.end()-1,led,pb,
                                  lastGcPointAdc,lastPh,lastPin);
    //set control variable
    pointFound=true;
  }

  //need to check if the pulse height was set to the maximum
  //because if it was then the pin value might be way to high
  //and mess up the linear with pin adc tuning
  //Need to correct the pin value when ph is the maximum
  if (lastPh==fMaxPossPh){
    //get iterators to the end of the vector
    ph=phVector.end()-1;
    adc=adcVector.end()-1;
    pin=pinVector.end()-1;
    

    MSG("LITuning_LastGc",Msg::kWarning)
      <<endl
      <<"("<<pb<<":"<<led<<")"
      <<" Looking for last GC point but PH has been set to max possible"
      <<endl
      <<"lastGcPointAdc="<<static_cast<Int_t>(lastGcPointAdc)
      <<", fLastGcPointAdc="<<static_cast<Int_t>(fLastGcPointAdc)
      <<", Difference="
      <<static_cast<Int_t>(fLastGcPointAdc-lastGcPointAdc)
      <<endl;
    this->PrintPhPinMsg(adc,ph,pin,lastPh,lastPin);

    //interpolate in pmt adc space (y-y2)*m+x2;
    Double_t m=(*adc-*(adc-1))/(*ph-*(ph-1));
    lastGcPointAdc=this->Interpolate(lastPh,*ph,m,*adc);
      
    //calculate gradient in pin adc space
    m=(*pin-*(pin-1))/(*ph-*(ph-1));
    lastPin=this->Interpolate(lastPh,*ph,m,*pin);

    //correct the values if needed
    this->CorrectInterpValues(adc,ph,pin,led,pb,lastGcPointAdc,
                              lastPh,lastPin);

    MSG("LITuning_LastGc",Msg::kWarning)
      <<endl
      <<"("<<pb<<":"<<led<<")"
      <<" After interpolation and potential correction:"
      <<"lastGcPointAdc="<<static_cast<Int_t>(lastGcPointAdc)
      <<endl;
    this->PrintPhPinMsg(adc,ph,pin,lastPh,lastPin);
  }

  MSG("LITuning",Msg::kVerbose)
    <<"CalcLastGcPoint method finished"<<endl;
}

//......................................................................

Int_t LITuning::GoodPin(vector<LIRun>::iterator currentGc)
{

  MSG("LITuning",Msg::kDebug)<<"Running GoodPin method..."<<endl;
  
  //get the two maximum pin values for the input gain curve  
  Double_t maxPin1=currentGc->GetMaxPin(1);
  Double_t maxPin2=currentGc->GetMaxPin(2);

  Int_t led=currentGc->GetLed();
  Int_t pb=currentGc->GetPb();

  Double_t avPin1dPINdPH=this->CalcAv_dPINdPH(currentGc,1);
  Double_t avPin2dPINdPH=this->CalcAv_dPINdPH(currentGc,2);
  
  if (MsgService::Instance()->IsActive("LITuning_Pin",Msg::kDebug)){
    MSG("LITuning_Pin", Msg::kInfo)
      <<"("<<pb<<":"<<led<<")"
      <<" Av gradient Pin1 = "<<avPin1dPINdPH
      <<", Av gradient Pin2 = "<<avPin2dPINdPH
      <<endl;
  }

  Bool_t pin1dPINdPHIsOk=false;
  Bool_t pin2dPINdPHIsOk=false;
  Bool_t pin1RangeIsOk=false;
  Bool_t pin2RangeIsOk=false;

  //work out why the different pins are good or bad
  if (maxPin1>fPinLowestMaxHG && maxPin1<fPinHighestMaxHG){
    pin1RangeIsOk=true;
  }
  if (maxPin2>fPinLowestMaxLG && maxPin2<fPinHighestMaxLG){
    pin2RangeIsOk=true;
  }
  if (avPin1dPINdPH>fMinHGdPINdPH){
    pin1dPINdPHIsOk=true;
  }
  if (avPin2dPINdPH>fMinLGdPINdPH){
    pin2dPINdPHIsOk=true;
  }
  
  if (MsgService::Instance()->IsActive("LITuning_Pin",Msg::kDebug)){
    MSG("LITuning_Pin",Msg::kInfo)
      <<"("<<pb<<":"<<led<<")"
      <<" Good/Bad: Pin1 range="<<pin1RangeIsOk
      <<", gradient="<<pin1dPINdPHIsOk
      <<endl<<"                  Pin2 range="<<pin2RangeIsOk
      <<", gradient="<<pin2dPINdPHIsOk<<endl
      <<"        maxPin1="<<maxPin1
      <<", maxPin2="<<maxPin2
      <<", avPin1dPINdPH="<<avPin1dPINdPH
      <<", avPin2dPINdPH="<<avPin2dPINdPH<<endl;
  }

  //for neardet, always return low gain pin first for the time being
  if (currentGc->GetDetector()==Detector::kNear
      && pin2RangeIsOk && pin2dPINdPHIsOk){
    MSG("LITuning_Pin",Msg::kDebug)
      <<"("<<pb<<":"<<led<<")"
      <<" It's ND, using second pin in calculations"<<endl;
    return 2;
  }


  //if pin1 good
  if (pin1RangeIsOk && pin1dPINdPHIsOk){
    MSG("LITuning_Pin",Msg::kDebug)
      <<"("<<pb<<":"<<led<<")"
      <<" Using first pin in calculations"<<endl;
    return 1;
  }

  //if pin2 good
  if (pin2RangeIsOk && pin2dPINdPHIsOk){
    MSG("LITuning_Pin",Msg::kDebug)
      <<"("<<pb<<":"<<led<<")"
      <<" Pin1 is bad, using second pin in calculations"<<endl;
    return 2;
  }

  //both have range and gradient bad to get this far

  //if both too high
  if (maxPin1>fPinHighestMaxHG && maxPin2>fPinHighestMaxLG){
    //print all details if not already done so
    if (!MsgService::Instance()->IsActive("LITuning_Pin",Msg::kDebug)){
      MSG("LITuning_Pin",Msg::kInfo)
        <<"("<<pb<<":"<<led<<")"
        <<" Good/Bad: Pin1 range="<<pin1RangeIsOk
        <<", gradient="<<pin1dPINdPHIsOk
        <<endl<<"                  Pin2 range="<<pin2RangeIsOk
        <<", gradient="<<pin2dPINdPHIsOk<<endl
        <<"        maxPin1="<<maxPin1
        <<", maxPin2="<<maxPin2
        <<", avPin1dPINdPH="<<avPin1dPINdPH
        <<", avPin2dPINdPH="<<avPin2dPINdPH<<endl;
    }
    MSG("LITuning_Pin",Msg::kWarning)
      <<endl<<"("<<pb<<":"<<led<<")"
      <<" Both Pin1 and pin2 are too high"
      <<endl<<"Selecting the lower of the two with good gradient"<<endl;

    if (maxPin1<maxPin2 && pin1dPINdPHIsOk){
      return 1;
    }
    else if (maxPin2<maxPin1 && pin2dPINdPHIsOk){
      return 2;
    }
  }

  if (!MsgService::Instance()->IsActive("LITuning_Pin",Msg::kDebug)){
    MSG("LITuning_Pin",Msg::kInfo)
      <<"("<<pb<<":"<<led<<")"
      <<" Good/Bad: Pin1 range="<<pin1RangeIsOk
      <<", gradient="<<pin1dPINdPHIsOk
      <<endl<<"                  Pin2 range="<<pin2RangeIsOk
      <<", gradient="<<pin2dPINdPHIsOk<<endl
      <<"        maxPin1="<<maxPin1
      <<", maxPin2="<<maxPin2
      <<", avPin1dPINdPH="<<avPin1dPINdPH
      <<", avPin2dPINdPH="<<avPin2dPINdPH<<endl;
  }

  //if you get this far print out a mesage anyway
  MSG("LITuning_Pin",Msg::kError)
    <<"("<<pb<<":"<<led<<")"<<" ** Both pins appear bad **"<<endl;

  MSG("LITuning",Msg::kDebug)<<"GoodPin method finished"<<endl;
  return 0;
}
/*

//early simple algorithm results for 17661
box06led01#I=62,96,136,195,254,300,343,383,423,462;
box06led02#I=43,405,625,844,1023,1023,1023,1023,1023,1023;
box06led03#I=60,97,157,233,308,380,449,518,587,656;
box06led04#I=64,116,181,246,308,366,422,476,529,580;
box06led05#I=67,103,142,187,233,279,329,375,420,456;
box06led06#I=63,100,147,202,256,313,367,418,464,509;
box06led07#I=63,104,157,212,264,310,351,385,418,445;
box06led08#I=63,99,139,187,236,287,338,387,434,478;
box06led09#I=63,101,149,202,253,301,347,394,438,480;

//before not allowing zero adc points to be used when interp
box06led10#I=29,87,127,180,239,287,332,375,416,458;
//after (57 is a much more realistic value from looking by eye)
box06led10#I=57,95,150,224,284,339,391,442,494,545;
//wohoooooooo

box06led11#I=18,84,129,178,226,268,306,343,377,411;
box06led12#I=62,88,115,148,184,218,252,287,322,358;
box06led13#I=63,85,112,143,175,211,250,292,333,367;
box06led14#I=63,172,240,279,310,341,371,401,431,461;

//before vetoing pin because of bad gradient:
box06led15#I=58,67,74,81,87,93,99,105,112,470;
//after 
box06led15#I=58,90,132,182,234,293,347,390,431,470;
//wohoooooooo

box06led16#I=61,89,122,166,209,249,288,327,363,399;
box06led17#I=66,116,181,249,312,371,427,482,536,589;
box06led18#I=62,101,151,207,263,316,368,419,470,521;
box06led19#I=65,126,191,255,315,374,431,485,538,591;
box06led20#I=63,104,150,201,249,296,341,387,433,478;

//after selecting last point based on pmt adc / ph gradient
//in preference to the fixed value
//and other changes inserted in above
box06led01#I=62,105,169,247,307,361,413,463,514,564;
box06led02#I=43,183,296,385,459,526,593,659,726,794;
box06led03#I=60,101,171,255,336,413,489,563,639,715;
box06led04#I=64,114,177,241,301,358,413,466,517,567;
box06led05#I=67,116,174,236,300,365,425,473,520,567;
box06led06#I=63,118,196,275,356,429,494,560,624,688;
box06led07#I=63,126,205,279,343,390,434,472,511,549;
box06led08#I=63,108,165,227,292,356,419,474,528,581;
box06led09#I=63,123,201,275,345,414,477,532,583,632;
box06led10#I=57,95,150,224,284,339,391,442,494,545;
box06led11#I=52,109,182,249,308,361,412,457,500,540;
box06led12#I=62,96,136,183,229,274,321,368,415,456;
box06led13#I=63,92,129,169,214,265,317,364,406,444;
box06led14#I=63,210,275,320,364,408,451,494,539,584;
box06led15#I=60,100,159,223,298,362,416,466,515,564;
box06led16#I=61,105,169,234,293,351,406,454,500,543;
box06led17#I=66,128,208,286,358,426,491,555,618,680;
box06led18#I=62,108,171,238,302,364,424,484,545,606;
box06led19#I=65,139,217,291,362,431,496,559,623,687;
box06led20#I=63,109,165,223,277,329,382,434,486,534;

*/
//......................................................................

void LITuning::CalculateDriftPoints()
{
  MSG("LITuning", Msg::kInfo) 
    <<"Running CalculateDriftPoints method..."<<endl;

  //get iterators to the gain curves and drift points
  vector<LIRun>::iterator gcData;
  gcData=fGainCurve->begin();
  vector<LIRun>::iterator dpTuned;
  //note can't get begin() before vector contains any elements
  //so get end()-1 for dpTuned below 

  MSG("LITuning_TuneDp",Msg::kDebug)
    <<"Starting the loop over all input data gain curves"<<endl;
    
  //loop over all the gain curves
  while (gcData!=fGainCurve->end()){

    if (!this->DataIsOk(gcData,2)){
      gcData++;
      continue;
    }

    //get iterators to the vectors of gain curve points
    vector<Double_t> phVector;
    vector<Double_t>::iterator ph;
    phVector=gcData->GetPh();
    ph=phVector.begin();
    vector<Double_t> adcVector;
    vector<Double_t>::iterator adc;
    adcVector=gcData->GetAdc();
    adc=adcVector.begin();

    //get the pulserbox and led associated with real data gain curve
    Int_t pulserBox=gcData->GetPb();
    Int_t led=gcData->GetLed();

    MSG("LITuning_TuneDp",Msg::kDebug)
      <<"("<<pulserBox<<":"<<led<<") First point has ph="<<*ph
      <<", adc="<<*adc<<endl;

    //add another DP to vector to hold tuned drift points
    fBestDp.push_back(LIRun(pulserBox,led,LIRun::kDriftPoint));
    //get iterator to the newest (last) element
    dpTuned=fBestDp.end()-1;
    MSG("LITuning_TuneDp",Msg::kDebug)
      <<"("<<pulserBox<<":"<<led<<") Created DP"<<endl;

    //check that the first adc is not too high
    if (*adc>fIdealAdc){  
      //if there is not enough data at low pulse heights
      //it is better to set too high than too low
      //since the led might not have turned on at low PH
      dpTuned->AddPoint(*ph,*adc);
      
      MSG("LITuning_TuneDp",Msg::kInfo)
        <<"("<<pulserBox<<":"<<led<<") DP"
        <<" Not enough points at low PHs"
        <<" for DP. Pulse Height set to lowest point taken"
        <<endl<<"PH="<<*ph<<", ideal ADC="<<fIdealAdc
        <<", actual adc="<<*adc
        <<endl;
      //once it's found then continue to the next GC
      gcData++;
      continue;
    }
    //iterate to take you to the second point
    ph++;
    adc++;

    MSG("LITuning_TuneDp",Msg::kDebug)
      <<"Starting loop over the given gain curve"<<endl;
    
    //loop over the points 2 -> end
    while(ph!=gcData->GetPh().end()){
      
      if (*adc>fIdealAdc){
        //interpolate to find the best PH
        Double_t m=(*ph-*(ph-1))/(*adc-*(adc-1));
        Double_t bestPhInterp=this->Interpolate(fIdealAdc,*adc,m,*ph);
        //(fIdealAdc-*adc)*
        //((*ph-*(ph-1))/(*adc-*(adc-1)))+*ph;
       
        if (bestPhInterp<0){ 
          MSG("LITuning_TuneDp",Msg::kWarning)
            <<endl<<"("<<pulserBox<<":"<<led<<") DP"
            <<" Interpolated point is negative ("
            <<bestPhInterp<<"). Setting to one"
            <<", ideal="<<fIdealAdc
            <<endl<<"       Lower point: adc="<<*(adc-1)
            <<", PH="<<*(ph-1)
            <<endl<<"       Upper point: adc="<<*adc
            <<", PH="<<*ph<<endl<<endl;
          
          //set to 0 if negative
          //this should never happen!!
          bestPhInterp=0;
        }
        else{
          MSG("LITuning_TuneDp",Msg::kInfo)
            <<"("<<pulserBox<<":"<<led<<") DP"
            <<" Interpolated point found, fIdealAdc="<<fIdealAdc
            <<endl<<"       Lower point: adc="<<*(adc-1)
            <<", PH="<<*(ph-1)
            <<endl<<"       Upper point: adc="<<*adc
            <<", PH="<<*ph<<"  Best Ph="<<bestPhInterp
            <<endl;
        }

        //set the interpolated DP
        dpTuned->AddPoint(bestPhInterp,fIdealAdc);
        
        //once point is found then stop loop    
        MSG("LITuning_TuneDp",Msg::kDebug)
          <<"("<<pulserBox<<":"<<led
          <<") Point found. Breaking out of loop"<<endl;
        break;
      }
      
      //iterate the gain curve point
      adc++;
      ph++;
      MSG("LITuning_TuneDp",Msg::kVerbose)
        <<"("<<pulserBox<<":"<<led<<") Iterated to next ph"<<endl;
    }//end of while

    //check if a DP was found and print a message
    if (*adc>10 && ph==gcData->GetPh().end() &&
        dpTuned->GetPh().empty()){
      MSG("LITuning_TuneDp",Msg::kInfo)
        <<"("<<pulserBox<<":"<<led<<") DP"
        <<" Not enough points at high PHs"
        <<" for DP, ideal ADC="<<fIdealAdc
        <<", PH="<<*ph<<", adc="<<*adc
        <<endl;
    }
    
    //iterate to next gain curve
    MSG("LITuning_TuneDp",Msg::kDebug)
      <<"("<<pulserBox<<":"<<led<<") Iterated to next dp"<<endl;
    gcData++;
  }

  MSG("LITuning_TuneDp",Msg::kDebug)
    <<"Finished the loop over all input data gain curves"<<endl;
    
  MSG("LITuning",Msg::kInfo) 
    <<"CalculateDriftPoints method finished"<<endl;
}

//......................................................................

void LITuning::InterpolateAboveOrBelow(EGcTuningType tuningType,
                                       vector<Double_t>::iterator adc,
                                       vector<Double_t>::iterator ph,
                                       vector<Double_t>::iterator pin,
                                       Int_t led,Int_t pb,
                                       Double_t& interpAdc,
                                       Double_t& interpPh,
                                       Double_t& interpPin)
{
  MSG("LITuning",Msg::kVerbose) 
    <<"Running InterpolateAboveOrBelow method..."<<endl;


  if (tuningType!=kGcLinearInPmtAdc && tuningType!=kGcLinearInPinAdc &&
      tuningType!=kGcUnknown){
    MSG("LITuning",Msg::kFatal)
      <<"Tuning type = "<<tuningType<<" not defined."<<endl
      <<"Will exit programme here."<<endl;
    exit(0);
  }
  
  //get iterators to the last points in the vectors
  //vector<Double_t>::iterator ph=phVector.end()-1;
  //vector<Double_t>::iterator adc=adcVector.end()-1;
  //vector<Double_t>::iterator pin=pinVector.end()-1;

  //check the tuning type since it dictates whether to
  //interpolate in pin or pmt adc space
  if (tuningType==kGcUnknown){
    //interpolate in pmt adc space (y-y2)*m+x2;
    Double_t m=(*ph-*(ph-1))/(*adc-*(adc-1));
    interpPh=this->Interpolate(interpAdc,*adc,m,*ph);
    
    //work out the pin value at the interpolated ph
    //calculate gradient in pin adc space
    m=(*pin-*(pin-1))/(*ph-*(ph-1));
    interpPin=this->Interpolate(interpPh,*ph,m,*pin);
  }
  else if (tuningType==kGcLinearInPmtAdc){
    //interpolate in pmt adc space (y-y2)*m+x2;
    Double_t m=(*ph-*(ph-1))/(*adc-*(adc-1));
    interpPh=this->Interpolate(interpAdc,*adc,m,*ph);    

    //work out the pin value at the interpolated ph
    //calculate gradient in pin adc space
    m=(*pin-*(pin-1))/(*ph-*(ph-1));
    interpPin=this->Interpolate(interpPh,*ph,m,*pin);
  }
  else if (tuningType==kGcLinearInPinAdc){
    //interpolate in pin space (y-y2)*m+x2;
    Double_t m=(*ph-*(ph-1))/(*pin-*(pin-1));
    interpPh=this->Interpolate(interpPin,*pin,m,*ph);

    //work out the pmt adc value at the interpolated ph
    //calculate gradient in pmt space 
    m=(*adc-*(adc-1))/(*ph-*(ph-1));
    interpAdc=this->Interpolate(interpPh,*ph,m,*adc);
  }

  Bool_t valuesCorrected=false;

  if (interpPh<fMinPossPh || interpPin<fMinPossPin ||
      interpPin>fMaxPossPin || interpPh>fMaxPossPh){
    this->CorrectInterpValues(adc,ph,pin,led,pb,interpAdc,
                              interpPh,interpPin);
    valuesCorrected=true;
  }

  if (valuesCorrected) {
    fS="Interpolated values were out of range and corrected ";
    fS+=" to be within allowed min/max";
  }
  else {
    fS="Interpolated values stayed within allowed min/max";
  }
  
  if (MsgService::Instance()->IsActive("LITuning",
                                       Msg::kDebug)){
    MSG("LITuning",Msg::kInfo)
      <<endl
      <<"("<<pb<<":"<<led<<")"
      <<" Looking for GC point."
      <<" Interpolated above/below data"
      <<endl<<fS<<endl
      <<"PMT ADC="<<interpAdc<<endl;
    this->PrintPhPinMsg(adc,ph,pin,interpPh,interpPin);
  }

  MSG("LITuning",Msg::kVerbose) 
    <<"InterpolateAboveOrBelow method finished"<<endl;
}

//......................................................................

void LITuning::CorrectInterpValues(vector<Double_t>::iterator adc,
                                   vector<Double_t>::iterator ph,
                                   vector<Double_t>::iterator pin,
                                   Int_t led,Int_t pb,
                                   Double_t interpAdc,
                                   Double_t& interpPh,
                                   Double_t& interpPin)
{
  MSG("LITuning",Msg::kVerbose) 
    <<"Running CorrectInterpValues method..."<<endl;

  //check points are above minimum
  if (interpPh<fMinPossPh || interpPin<fMinPossPin){
    MSG("LITuning_TuneHl",Msg::kWarning)
      <<endl
      <<"("<<pb<<":"<<led<<")"
      <<" Looking for GC point."
      <<" Interpolated but overshot allowed minimum"
      <<endl
      <<"PMT ADC="<<interpAdc<<endl;
    this->PrintPhPinMsg(adc,ph,pin,interpPh,interpPin);

    MSG("LITuning_TuneHl",Msg::kWarning)
      <<endl
      <<"After setting overshoot points to minimum:"<<endl;
    //set points to minimum if required
    if (interpPh<fMinPossPh) interpPh=fMinPossPh;
    if (interpPin<fMinPossPin) interpPin=fMinPossPin;
    //print out their new values
    this->PrintPhPinMsg(adc,ph,pin,interpPh,interpPin);
  }

  //check points are below maximum
  if (interpPin>fMaxPossPin || interpPh>fMaxPossPh){
    MSG("LITuning_TuneHl",Msg::kWarning)
      <<endl
      <<"("<<pb<<":"<<led<<")"
      <<" Looking for GC point."
      <<" Interpolated but overshot allowed maximum"
      <<endl
      <<"PMT ADC="<<interpAdc<<endl;
    this->PrintPhPinMsg(adc,ph,pin,interpPh,interpPin);

    MSG("LITuning_TuneHl",Msg::kWarning)
      <<endl
      <<"After setting overshoot points to maximum:"<<endl;
    //set points to maximum if required
    if (interpPh>fMaxPossPh) interpPh=fMaxPossPh;
    if (interpPin>fMaxPossPin) interpPin=fMaxPossPin;
    //print out their new values
    this->PrintPhPinMsg(adc,ph,pin,interpPh,interpPin);
  }
  
  MSG("LITuning",Msg::kVerbose) 
    <<"CorrectInterpValues method finished"<<endl;
}

//......................................................................

Double_t LITuning::Interpolate(Double_t y,Double_t y2,Double_t m,
                               Double_t x2)
{
  MSG("LITuning",Msg::kVerbose) 
    <<"Running Interpolate method..."<<endl;
    
  //interpolation formula
  Double_t x=(y-y2)*m+x2;

  MSG("LITuning",Msg::kVerbose) 
    <<"Interpolate method finished"<<endl;
  return x;
}

//......................................................................

Bool_t LITuning::DataIsOk(vector<LIRun>::iterator gcData,
                          Int_t minNumPoints,Int_t maxAllowedAdc)
{
  MSG("LITuning",Msg::kVerbose) 
    <<"Running DataIsOk method..."<<endl;
    
  Bool_t sizeOk=false;
  Bool_t adcRangeOk=false;

  Int_t phSize=gcData->GetPh().size();
  Int_t adcSize=gcData->GetAdc().size();

  //DO NOT use a gain curve with less than minNumPoints points
  if (phSize>=minNumPoints && adcSize>=minNumPoints &&
      phSize==adcSize) sizeOk=true;

  //check the range of the adc values
  //this is no good for ND
  //hack alert!!!!
  maxAllowedAdc=0;//just for unused parameter warning
  //if (gcData->GetMaxAdc()<maxAllowedAdc && gcData->GetMaxAdc()>=1 && 
  //  gcData->GetMinAdc()<maxAllowedAdc && 
  //  gcData->GetMinAdc()>=0) adcRangeOk=true; 

  //check the range of the adc values
  if (gcData->GetMaxAdc()>=1 && gcData->GetMinAdc()>=0) adcRangeOk=true;

  MSG("LITuning",Msg::kDebug) 
    <<"DataIsOk: max ADC="<<gcData->GetMaxAdc()
    <<", min ADC="<<gcData->GetMinAdc()
    <<", PH size="<<phSize<<", ADC size="<<adcSize
    <<endl;

  MSG("LITuning",Msg::kVerbose) 
    <<"DataIsOk method finished"<<endl;
  return sizeOk*adcRangeOk;//only ok if both ok
}

//......................................................................

void LITuning::GcLinearInPmtAdc(vector<LIRun>::iterator gcData,
                                Int_t goodPin,
                                vector<LIRun>::iterator& gcTuned,
                                Double_t firstPh,Double_t lastPh,
                                Double_t firstPin,Double_t lastPin)
{
  MSG("LITuning",Msg::kVerbose)
    <<"Running GcLinearInPmtAdc method..."<<endl;

  //get iterators to the vectors of gain curve points
  vector<Double_t> phVector;
  vector<Double_t>::iterator ph;
  phVector=gcData->GetPh();
  ph=phVector.begin();
  vector<Double_t> adcVector;
  vector<Double_t>::iterator adc;
  adcVector=gcData->GetAdc();
  adc=adcVector.begin();
  vector<Double_t>pinVector;
  vector<Double_t>::iterator pin;
  pinVector=gcData->GetPin(goodPin);
  pin=pinVector.begin();

  //get the pulserbox and led associated with real data gain curve
  Int_t pb=gcData->GetPb();
  Int_t led=gcData->GetLed();

  //set the first point as it has already been found
  gcTuned->AddPoint(firstPh,fFirstGcPointAdc,firstPin);

  //loop to find gain curve points (second through to last but one)
  for (Int_t g=1;g<fNumGainPoints-1;g++){
    //calculate the required adc for this point
    Double_t adcSeparation=(fLastGcPointAdc-fFirstGcPointAdc)/
      (fNumGainPoints-1);
    Double_t reqAdc=fFirstGcPointAdc+g*adcSeparation;
      
    if (reqAdc<fMinPossAdc || reqAdc>fMaxPossAdc){
      MSG("LITuning_TuneGc",Msg::kFatal)
        <<endl<<"("<<pb<<":"<<led<<") GC g="<<g+1
        <<" Aborted finding gain curve point, reqAdc="<<reqAdc
        <<", adcSeparation="<<adcSeparation<<endl;
      continue;
    }
      
    if (MsgService::Instance()->IsActive("LITuning_TuneGc",
                                         Msg::kDebug)){
      MSG("LITuning_TuneGc",Msg::kInfo)
        <<endl<<"Gain Curve point "<<g+1<<" ("<<pb<<":"<<led
        <<") reqAdc="<<reqAdc<<", sep="<<adcSeparation<<endl;
    }

    //reset all the variables for this loop
    Bool_t pointFound=false;
    Bool_t firstPointTooHigh=false;
    Int_t counter=0;//already done one point
    ph=phVector.begin();
    adc=adcVector.begin();
    Double_t interpPh=-1;
    Double_t pinAtReqAdc=-1;

    //loop over the points
    while(ph!=phVector.end()){
      //counter
      counter++;

      //check if the first data point has an adc value greater than
      //the last required point
      if (*adc>reqAdc && adc==adcVector.begin()){
            
        //set control variables
        pointFound=false;
        firstPointTooHigh=true;
          
        //jump out of the loop
        break;
      }
    
      //check if the current point has an adc value greater than the
      //required point
      //and whether a point has been found yet
      if (*adc>reqAdc && !pointFound){
        
        //set control variable
        pointFound=true;

        //interpolate in pmt adc space (y-y2)*m+x2;
        Double_t m=(*ph-*(ph-1))/(*adc-*(adc-1));
        interpPh=this->Interpolate(reqAdc,*adc,m,*ph);
      
        if (MsgService::Instance()->IsActive("LITuning_TuneGc",
                                             Msg::kDebug)){
          MSG("LITuning_TuneGc",Msg::kInfo)
            <<"("<<pb<<":"<<led<<")"
            <<" GC point "<<g+1<<" found"
            <<" ("<<counter<<"/"<<adcVector.size()<<")"
            <<", Required Point="<<reqAdc<<endl;
          this->PrintPhMsg(adc,ph,interpPh);
        }

        //jump out of the loop when point is found
        break;
      }
    
      //go to the next point
      ph++;
      adc++;
      pin++;
    }

    //check on success of finding point
    if (firstPointTooHigh){
      this->InterpolateAboveOrBelow(kGcLinearInPmtAdc,
                                    adcVector.begin()+1,
                                    phVector.begin()+1,
                                    pinVector.begin()+1,led,pb,
                                    reqAdc,interpPh,pinAtReqAdc);
      //set control variable
      pointFound=true;
    }
    else if (!pointFound){//i.e. last data point was too low
      this->InterpolateAboveOrBelow(kGcLinearInPmtAdc,
                                    adcVector.end()-1,phVector.end()-1,
                                    pinVector.end()-1,led,pb,
                                    reqAdc,interpPh,pinAtReqAdc);
      //set control variable
      pointFound=true;
    }
      
    //add the point to the tuned GC
    if (pointFound) gcTuned->AddPoint(interpPh,reqAdc,pinAtReqAdc);
    else MSG("LITuning_TuneGc",Msg::kWarning)<<"No GC Point found"
                                             <<endl;
  }//end of for

  //set the last point as it has already been found
  gcTuned->AddPoint(lastPh,fLastGcPointAdc,lastPin);

  MSG("LITuning",Msg::kVerbose) 
    <<"GcLinearInPmtAdc method finished"<<endl;
}

//......................................................................

void LITuning::GcLinearInPinAdc(vector<LIRun>::iterator gcData,
                                Int_t goodPin,
                                vector<LIRun>::iterator& gcTuned,
                                Double_t firstPh,Double_t lastPh,
                                Double_t firstPin,Double_t lastPin)
{
  MSG("LITuning",Msg::kVerbose) 
    <<"Running GcLinearInPinAdc method..."<<endl;

  //get iterators to the vectors of gain curve points
  vector<Double_t> phVector;
  vector<Double_t>::iterator ph;
  phVector=gcData->GetPh();
  ph=phVector.begin();
  vector<Double_t> adcVector;
  vector<Double_t>::iterator adc;
  adcVector=gcData->GetAdc();
  adc=adcVector.begin();
  vector<Double_t>pinVector;
  vector<Double_t>::iterator pin;
  pinVector=gcData->GetPin(goodPin);
  pin=pinVector.begin();

  //get the pulserbox and led associated with real data gain curve
  Int_t pb=gcData->GetPb();
  Int_t led=gcData->GetLed();

  //set the first point as it has already been found
  gcTuned->AddPoint(firstPh,fFirstGcPointAdc,firstPin);

  //loop to find gain curve points (second through to last but one)
  for (Int_t g=1;g<fNumGainPoints-1;g++){
    //calculate the required pin for this point
    Double_t pinSeparation=(lastPin-firstPin)/(fNumGainPoints-1);
    Double_t reqPin=firstPin+g*pinSeparation;
    
    if (reqPin<fMinPossPin || reqPin>fMaxPossPin){
      MSG("LITuning_TuneGc",Msg::kFatal)
        <<endl<<"("<<pb<<":"<<led<<") GC g="<<g+1
        <<" Aborted finding gain curve point, reqPin="<<reqPin
        <<", pinSeparation="<<pinSeparation<<endl;
      continue;
    }
    if (MsgService::Instance()->IsActive("LITuning_TuneGc",
                                         Msg::kDebug)){
      MSG("LITuning_TuneGc",Msg::kInfo)
        <<endl<<"Gain Curve point "<<g+1<<" ("<<pb<<":"<<led
        <<") reqPin="<<reqPin<<", sep="<<pinSeparation<<endl;
      if (MsgService::Instance()->IsActive("LITuning_TuneGc",
                                           Msg::kVerbose)){
        this->PrintFirstLastMsg(firstPh,lastPh,firstPin,lastPin);
      }
    }

    //reset all the variables for this loop
    Bool_t pointFound=false;
    Bool_t firstPointTooHigh=false;
    Int_t counter=0;//already done one point
    adc=adcVector.begin();
    ph=phVector.begin();
    pin=pinVector.begin();
    Double_t interpPh=-1;
    Double_t adcAtReqPin=-1;
    
    //loop over the points
    while(ph!=phVector.end()){
      //counter
      counter++;
      
      //check if the first data point has an pin value greater than
      //the last required point
      if (*pin>reqPin && pin==pinVector.begin()){
        
        //set control variables
        pointFound=false;
        firstPointTooHigh=true;
          
        //jump out of the loop
        break;
      }
    
      //check if the current point has an pin value greater than the
      //required point
      //and whether a point has been found yet
      if (*pin>reqPin && !pointFound){
        
        //set control variable
        pointFound=true;

        //interpolate in pin space (y-y2)*m+x2;
        Double_t m=(*ph-*(ph-1))/(*pin-*(pin-1));
        interpPh=this->Interpolate(reqPin,*pin,m,*ph);
        
        //calculate gradient in pmt space 
        m=(*adc-*(adc-1))/(*ph-*(ph-1));
        //interpolate in pmt space (y-y2)*m+x2;
        adcAtReqPin=this->Interpolate(interpPh,*ph,m,*adc);
     
        if (MsgService::Instance()->IsActive("LITuning_TuneGc",
                                             Msg::kDebug)){
          MSG("LITuning_TuneGc",Msg::kInfo)
            <<"("<<pb<<":"<<led<<")"
            <<" GC point "<<g+1<<" found"
            <<" ("<<counter<<"/"<<pinVector.size()<<")"
            <<", Required PIN ADC="<<reqPin<<endl;
          this->PrintPinMsg(pin,ph,adcAtReqPin,interpPh);
        }
          
        //jump out of the loop when point is found
        break;
      }
    
      //go to the next point
      ph++;
      pin++;
      adc++;
    }

    //check on success of finding point
    if (firstPointTooHigh){
      this->InterpolateAboveOrBelow(kGcLinearInPinAdc,
                                    adcVector.begin()+1,//was pinVector?
                                    phVector.begin()+1,
                                    pinVector.begin()+1,led,pb,
                                    adcAtReqPin,interpPh,reqPin);
      //set control variable
      pointFound=true;
    }
    else if (!pointFound){//i.e. last data point was too low
      this->InterpolateAboveOrBelow(kGcLinearInPinAdc,
                                    adcVector.end()-1,phVector.end()-1,
                                    pinVector.end()-1,led,pb,
                                    adcAtReqPin,interpPh,reqPin);
      //set control variable
      pointFound=true;
    }
      
    //add the point to the tuned GC
    if (pointFound) gcTuned->AddPoint(interpPh,adcAtReqPin,reqPin);
    else MSG("LITuning_TuneGc",Msg::kWarning)<<"No GC point found"
                                             <<endl;

  }//end of for

  //set the last point as it has already been found
  gcTuned->AddPoint(lastPh,fLastGcPointAdc,lastPin);

  MSG("LITuning",Msg::kVerbose) 
    <<"GcLinearInPinAdc method finished"<<endl;
}

//......................................................................

void LITuning::CalculateGainCurve(EGcTuningType tuningType)
{
  MSG("LITuning", Msg::kInfo) 
    <<"Running CalculateGainCurve method..."<<endl;

  if (tuningType!=kGcLinearInPmtAdc && tuningType!=kGcLinearInPinAdc){
    MSG("LITuning_TuneGc",Msg::kFatal)
      <<"Tuning type = "<<tuningType<<" not defined."<<endl
      <<"Will exit programme here."<<endl;
    exit(0);
  }
  
  //get iterators to the gain curve
  vector<LIRun>::iterator gcData;
  gcData=fGainCurve->begin();
  vector<LIRun>::iterator gcTuned;
  //note can't get begin() before vector contains any elements
  //so get end()-1 for gcTuned below 

  //loop over all the gain curves for the different pulser boxes
  //and leds
  while (gcData!=fGainCurve->end()){

    //get the pulserbox and led associated with real data gain curve
    Int_t pulserBox=gcData->GetPb();
    Int_t led=gcData->GetLed();

    if (led!=-1 && pulserBox!=-1){
      if (MsgService::Instance()->IsActive("LITuning_TuneGc",
                                           Msg::kDebug)){
        MSG("LITuning_TuneGc", Msg::kInfo)
          <<endl<<"("<<pulserBox<<":"<<led
          <<") **** Starting GC tuning for this led... **** "
          <<endl<<endl;
      }
    }

    if (!this->DataIsOk(gcData)){
      MSG("LITuning_TuneGc", Msg::kInfo) 
        <<"This gain curve data ("<<pulserBox<<":"<<led<<") is not "
        <<"consistent for calcuations"<<endl
        <<"Either it has less than the minimum number of points or"
        <<" the data is nonsense"<<endl;
      gcData++;
      continue;
    }

    //get iterators to the vectors of gain curve points
    vector<Double_t> phVector;
    vector<Double_t>::iterator ph;
    phVector=gcData->GetPh();
    ph=phVector.begin();
    vector<Double_t> adcVector;
    vector<Double_t>::iterator adc;
    adcVector=gcData->GetAdc();
    adc=adcVector.begin();
    vector<Double_t>pinVector;
    vector<Double_t>::iterator pin;
    Int_t goodPin=this->GoodPin(gcData);
    if (goodPin==1 || goodPin==2){
      pinVector=gcData->GetPin(goodPin);
      pin=pinVector.begin();
    }
    else {
      MSG("LITuning_TuneGc", Msg::kWarning)
        <<endl 
        <<"This gain curve ("<<pulserBox<<":"<<led<<") does not have "
        <<"a good pin diode associated with it (goodPin="<<goodPin<<")"
        <<endl;
      
      //set pin iterator to point to pulse height!!!
      //pin=ph;
      //don't like this, best to skip

      //find out the first and last points before continuing
      
      //variables for first and last points
      Double_t firstPh=-1;
      Double_t firstPin=-1;
      Double_t lastPh=-1;
      Double_t lastPin=-1;
      //get the first and last gain curve points
      this->CalcFirstGcPoint(gcData,1,firstPh,firstPin);
      this->CalcLastGcPoint(gcData,1,lastPh,lastPin);
      
      MSG("LITuning_TuneGc", Msg::kWarning)
        <<"("<<pulserBox<<":"<<led<<") Skipping this LED"<<endl
        <<"Calcuated first and last points anyway:"<<endl
        <<"first Ph="<<firstPh<<", last Ph="<<lastPh<<endl
        <<"first Pin="<<firstPin<<", last Pin="<<lastPin<<endl<<endl;

      gcData++;
      continue;
    }
    
    MSG("LITuning_TuneGc",Msg::kDebug)
      <<"("<<pulserBox<<":"<<led<<"), ph="<<*ph<<", adc="<<*adc
      <<", good pin="<<goodPin<<endl;

    //variables to hold the ph and pin values to find GC points between
    Double_t firstPh=-1;
    Double_t firstPin=-1;
    Double_t lastPh=-1;
    Double_t lastPin=-1;

    //get the first and last gain curve points
    this->CalcFirstGcPoint(gcData,goodPin,firstPh,firstPin);
    this->CalcLastGcPoint(gcData,goodPin,lastPh,lastPin);
    
    if (MsgService::Instance()->IsActive("LITuning_TuneGc",
                                         Msg::kDebug)){
      this->PrintFirstLastMsg(firstPh,lastPh,firstPin,lastPin);
    }

    //add another GC to vector to hold tuned drift points
    fBestGc.push_back(LIRun(pulserBox,led,LIRun::kGainCurve));
    //get iterator to the newest (last) element
    gcTuned=fBestGc.end()-1;

    //Run the tuning algorithm on this particular pulser box and
    //led to generate the points
    if (tuningType==kGcLinearInPmtAdc){
      this->GcLinearInPmtAdc(gcData,goodPin,gcTuned,firstPh,lastPh,
                             firstPin,lastPin);
    }
    else if (tuningType==kGcLinearInPinAdc){
      this->GcLinearInPinAdc(gcData,goodPin,gcTuned,firstPh,lastPh,
                             firstPin,lastPin);
    }

    MSG("LITuning_TuneGc",Msg::kDebug)
      <<"("<<pulserBox<<":"<<led<<") Created tuned GC Point"<<endl;
    
    gcData++;
  }

  MSG("LITuning",Msg::kInfo) 
    <<"CalculateGainCurve method finished"<<endl;
}

//......................................................................

Int_t LITuning::PrintLIConfig(LIRun::ELIRunType liRun)
{
  MSG("LITuning",Msg::kDebug) 
    <<"Running PrintLIConfig method..."<<endl;

  vector<LIRun>::iterator tunedRun;
  vector<LIRun>::iterator tunedRunEnd;

  //select appropriate vectors depending if the run to 
  //be printed is a GC or a DP
  if (liRun==LIRun::kGainCurve){
    if (fBestGc.empty()){
      MSG("LITuning",Msg::kWarning)
        <<"There is nothing to print, no tuned values have been"
        <<" computed for a gain curve"<<endl;
      return 0;
    }
    //select the tuned gain curve
    tunedRun=fBestGc.begin();
    tunedRunEnd=fBestGc.end();
    MSG("LITuning",Msg::kInfo)<<endl<<"<fullCalib>"<<endl;
  }
  else if (liRun==LIRun::kDriftPoint){
    if (fBestDp.empty()){
      MSG("LITuning",Msg::kWarning)
        <<"There is nothing to print, no tuned values have been"
        <<" computed for a drift point"<<endl;
      return 0;
    }
    //select the tuned drift points
    tunedRun=fBestDp.begin();
    tunedRunEnd=fBestDp.end();
    MSG("LITuning",Msg::kInfo)<<endl<<"<driftCalib>"<<endl;
  }
  
  //variable to hold the last pulser box
  Int_t lastPulserBox=-1;

  //loop and print out the list of LI run pulse heights 
  //in the desired format
  while (tunedRun!=tunedRunEnd){

    //get iterators to the vectors of li run points
    vector<Double_t> phVector;
    vector<Double_t>::iterator ph;
    phVector=tunedRun->GetPh();
    ph=phVector.begin();

    //check that there is a point present
    if (phVector.empty()){
      MSG("LITuning",Msg::kWarning)
        <<"("<<tunedRun->GetPb()<<":"<<tunedRun->GetLed()
        <<") DP is empty"<<endl;
      tunedRun++;
      continue;
    }
    else{
      //leave a line between pulser boxes
      if (tunedRun->GetPb()!=lastPulserBox){
        MSG("LITuning",Msg::kInfo)<<endl;
      }
    }

    //keep a record of last pulser box
    lastPulserBox=tunedRun->GetPb();

    //print only on debug
    if (MsgService::Instance()->IsActive("LITuning",Msg::kDebug)){
      tunedRun->PrintAll();
    }

    string sBox="";
    sBox=Form("%d",tunedRun->GetPb());
    if (tunedRun->GetPb()<10) sBox="0"+sBox;
    
    string sLed="";
    sLed=Form("%d",tunedRun->GetLed());
    if (tunedRun->GetLed()<10) sLed="0"+sLed;
    
    //form the base string
    fS="box"+sBox+"led"+sLed+"#I=";

    //loop over all the points in the vector (1 or more)
    while(ph!=phVector.end()){

      MSG("LITuning",Msg::kVerbose)
        <<"Looping over tuned points"<<endl;

      //round up pulse height
      Double_t rem=0;
      Double_t iph=0;
      rem=modf(*ph,&iph);
      if (rem>=0.5) iph++;//round up

      string sPh="";
      if (iph<1e6){
        sPh=Form("%d",static_cast<Int_t>(iph));
      }
      else{
        sPh="Really high PH!";
      }

      MSG("LITuning",Msg::kVerbose)
        <<"iph="<<iph
        <<", sPh="<<sPh<<endl;

      fS+=sPh;
      if (ph<phVector.end()-1) fS+=",";
      //iterate the next point
      ph++;
    }
    fS+=";";
    //print out the string in the required format
    MSG("LITuning",Msg::kInfo)<<fS<<endl; 

    //iterate the li run that you are printing out
    tunedRun++;
  }

  //print out end tag
  if (liRun==LIRun::kGainCurve){
    MSG("LITuning",Msg::kInfo)<<endl<<"</fullCalib>"<<endl;
  }
  else if (liRun==LIRun::kDriftPoint){
    MSG("LITuning",Msg::kInfo)<<endl<<"</driftCalib>"<<endl;
  }
  
  //This is the format under which the parameters are specified
  //in the config file

  //Label under which we specify LI sequence for drift calib:
  //<driftCalib>
  //box00led01#I=400;
  //box00led02#I=440;
  //box00led03#I=510;
  //[etc., etc.]
  //box15led20#I=390;
  //</driftCalib>
  
  //Label under which we specify LI sequence for gain-curve calib:
  //<fullCalib>
  //box00led01#I=100,200,300,400,500,600,700,800,900,1000;
  //box00led02#I=110,220,330,440,550,660,770,880,990,1023;
  //[etc., etc.]
  //box15led20#I=98,196,294,392,490,588,686,784,882,980;
  //</fullCalib>
  
  MSG("LITuning",Msg::kDebug) 
    <<"PrintLIConfig method finished"<<endl;
  return 1;
}

//......................................................................

void LITuning::PrintLedCheckGrid(LIRun::ELIRunType liRun,
                                 Int_t firstPb,Int_t lastPb,
                                 Int_t firstLed,Int_t lastLed)
{
  MSG("LITuning",Msg::kDebug) 
    <<"Running PrintLedCheckGrid method..."<<endl;
 
  vector<LIRun>::iterator tunedRun;
  vector<LIRun>::iterator tunedRunEnd;

  Int_t totalNumPbs=lastPb-firstPb+1;
  Int_t numLedsPerPb=lastLed-firstLed+1;
  
  //select appropriate vectors depending if the run to 
  //be checked is a GC or a DP
  if (liRun==LIRun::kGainCurve){
    if (fBestGc.empty()){
      MSG("LITuning",Msg::kWarning)
        <<"There is nothing to print, no tuned values have been"
        <<" computed for a gain curve"<<endl;
    }
    //select the tuned gain curve
    tunedRun=fBestGc.begin();
    tunedRunEnd=fBestGc.end();
  }
  else if (liRun==LIRun::kDriftPoint){
    if (fBestDp.empty()){
      MSG("LITuning",Msg::kWarning)
        <<"There is nothing to print, no tuned values have been"
        <<" computed for a drift point"<<endl;
    }
    //select the tuned drift points
    tunedRun=fBestDp.begin();
    tunedRunEnd=fBestDp.end();
  }
  
  //vector to hold which leds were tuned
  vector<Int_t> totalLISystem(totalNumPbs*numLedsPerPb,0);

  //loop over tuned leds
  while (tunedRun!=tunedRunEnd){

    //get iterator to the vector of li run points
    vector<Double_t> phVector;
    phVector=tunedRun->GetPh();

    //check that there is a point present
    if (phVector.empty()){
      MSG("LITuning",Msg::kWarning)
        <<"("<<tunedRun->GetPb()<<":"<<tunedRun->GetLed()
        <<") DP is empty"<<endl;
      tunedRun++;
      continue;
    }

    Int_t pb=tunedRun->GetPb();
    Int_t led=tunedRun->GetLed();
    Int_t ledIndex=pb*numLedsPerPb+led-1;

    if (pb<=lastPb && led<=lastLed &&
        pb>=firstPb && led>=firstLed){
      totalLISystem[ledIndex]++;
    }
    else{
      MSG("LITuning",Msg::kInfo)
        <<"Pulser box and or led are out of expected range:"
        <<" PB="<<pb<<", LED="<<led<<endl;
    }
    tunedRun++;
  }
  
  //variable to hold the output string
  string sDataProcessed="";

  //loop over whole li system
  for (Int_t pb=firstPb;pb<=lastPb;pb++){
    for (Int_t led=firstLed;led<=lastLed;led++){
      Int_t ledIndex=pb*numLedsPerPb+led-1;
      MSG("LITuning",Msg::kVerbose)
        <<"PB="<<pb<<", LED="<<led<<endl;
      
      string sBox="";
      sBox=Form("%d",pb);
      if (pb<10) sBox="0"+sBox;
      
      string sLed="";
      sLed=Form("%d",led);
      if (led<10) sLed="0"+sLed;
      
      string sTotalLISystem=Form("%d",totalLISystem[ledIndex]);
      static Int_t lastPulserBox=-1;

      if (pb!=lastPulserBox){
        MSG("LITuning",Msg::kDebug)
          <<"New PB="<<pb<<", LED="<<led<<endl;
        
        //add prefix to each line in string
        if (pb==firstPb){
          sDataProcessed+="PB "+sBox+" LEDs: ";
        }
        else{
          sDataProcessed+="\nPB "+sBox+" LEDs: ";
        }
      }
      //keep a record of last pulser box
      lastPulserBox=pb;

      //add whether point was tuned      
      sDataProcessed+=sTotalLISystem;
      
      //add comma separations
      if (led!=lastLed){
        sDataProcessed+=",";
      }
      
      //add spaces every five leds
      if (led%5==0){
        sDataProcessed+=" ";
      }
    }
  }

  string sLIRunType="non-specified type";
  if (liRun==LIRun::kGainCurve){
    sLIRunType="Gain Curve";
  }
  else if (liRun==LIRun::kDriftPoint){
    sLIRunType="Drift Point";
  }

  //print string
  MSG("LITuning",Msg::kInfo)
    <<endl<<"Of the "<<totalNumPbs*numLedsPerPb
    <<" leds the following were tuned for a "<<sLIRunType<<endl
    <<sDataProcessed<<endl;
  
  MSG("LITuning",Msg::kDebug) 
    <<"Finished PrintLedCheckGrid method"<<endl;
}

//......................................................................

Bool_t LITuning::GetTunedDpGraphs(TGraph*& gAdcVsLed,TGraph*& gPhVsLed,
                                  TGraph*& gPinVsLed,
                                  Int_t numLedsPerPb)
{
  MSG("LITuning",Msg::kInfo) 
    <<"Running GetTunedDpGraphs method..."<<endl;

  //check if a gain curve has been calculated
  if (fBestDp.empty()){
    MSG("LITuning",Msg::kWarning)
      <<"There is nothing to print, no tuned values have been"
      <<" computed for a drift point"<<endl;
    return false;
  }

  vector<LIRun>::iterator tunedRun;
  vector<LIRun>::iterator tunedRunEnd;
 
  //select the tuned gain curve
  tunedRun=fBestDp.begin();
  tunedRunEnd=fBestDp.end();
  
  vector<Double_t> phVector;
  vector<Double_t> pinVector;
  vector<Double_t> adcVector;
  vector<Double_t> ledIndexVector;

  //loop over tuned DPs and fill vectors
  while (tunedRun!=tunedRunEnd){

    //get the pulserbox and led 
    Int_t pulserBox=tunedRun->GetPb();
    Int_t led=tunedRun->GetLed();
    //calculate the ledIndex
    Int_t ledIndex=pulserBox*numLedsPerPb+led-1;

    phVector.push_back(tunedRun->GetPh(0));
    pinVector.push_back(tunedRun->GetPin(2,0));//pin2
    adcVector.push_back(tunedRun->GetAdc(0));
    ledIndexVector.push_back(ledIndex);

   //print only on debug
    if (MsgService::Instance()->IsActive("LITuning",Msg::kDebug)){
      tunedRun->PrintAll();
    }

    //iterate the li run that you are looking at
    tunedRun++;
  }

  MSG("LITuning",Msg::kDebug) 
    <<"Making graphs..."<<endl;
  gAdcVsLed=this->TGraphVect(ledIndexVector,adcVector);
  gPhVsLed=this->TGraphVect(ledIndexVector,phVector);
  gPinVsLed=this->TGraphVect(ledIndexVector,pinVector);
    
  //do pmt adc graph
  fS="Theoretical Tuned PMT ADC of LEDs";
  gAdcVsLed->SetTitle(fS.c_str());
  gAdcVsLed->SetFillColor(0);
  gAdcVsLed->SetMinimum(-2);

  //do ph graph
  fS="Tuned Pulse Height of LEDs";
  gPhVsLed->SetTitle(fS.c_str());
  gPhVsLed->SetFillColor(0);
  gPhVsLed->SetMinimum(-2);

  //do pin adc graph
  fS="PIN ADC at Tuned Drift Point Light Levels";
  gPinVsLed->SetTitle(fS.c_str());
  gPinVsLed->SetFillColor(0);
  gPinVsLed->SetMinimum(-2);

  MSG("LITuning",Msg::kDebug) 
    <<"Finished GetTunedDpGraphs method"<<endl;
  return true;
}

//......................................................................

Bool_t LITuning::GetDataGcGraphs(vector<TGraph*>& gAdcVsPh,
                                 vector<TGraph*>& gAdcFVsPh,
                                 vector<TGraph*>& gPinVsPh,
                                 vector<TGraph*>& gPin2VsPh,
                                 vector<TGraphAsymmErrors*>& 
                                 gAdcErrVsPh,
                                 vector<TGraphAsymmErrors*>& 
                                 gAdcFErrVsPh)
{
  MSG("LITuning",Msg::kInfo) 
    <<"Running GetDataGcGraphs method..."<<endl;

  //check if data has been input
  if (!fGainCurve){//check that it's not a null pointer
    MSG("LITuning",Msg::kWarning)
      <<"There is nothing to print, no data has been input"<<endl;
    return false;
  }

  //get iterators for the histos
  vector<TGraph*>::iterator gAdcVsPhIter;
  vector<TGraph*>::iterator gAdcFVsPhIter;
  vector<TGraph*>::iterator gPinVsPhIter;
  vector<TGraph*>::iterator gPin2VsPhIter;
  vector<TGraphAsymmErrors*>::iterator gAdcErrVsPhIter;
  vector<TGraphAsymmErrors*>::iterator gAdcFErrVsPhIter;

  vector<LIRun>::iterator gcData;
  vector<LIRun>::iterator gcDataEnd;
  
  //select the tuned gain curve
  gcData=fGainCurve->begin();
  gcDataEnd=fGainCurve->end();

  //loop over input data and create graphs
  while (gcData!=gcDataEnd){

    //get iterators to the vectors of li run points
    vector<Double_t> phVector;
    phVector=gcData->GetPh();
    vector<Double_t> pinVector;
    pinVector=gcData->GetPin(1);
    vector<Double_t> pin2Vector;
    pin2Vector=gcData->GetPin(2);
    vector<Double_t> adcVector;
    adcVector=gcData->GetAdc();
    vector<Double_t> adcLowVector;
    adcLowVector=gcData->GetAdcLow();
    vector<Double_t> adcHighVector;
    adcHighVector=gcData->GetAdcHigh();
    vector<Double_t> adcFLowVector;
    adcFLowVector=gcData->GetAdcLowF();
    vector<Double_t> adcFHighVector;
    adcFHighVector=gcData->GetAdcHighF();
    vector<Double_t> adcFVector;
    adcFVector=gcData->GetAdcF();

    //get the pulserbox and led 
    Int_t pulserBox=gcData->GetPb();
    Int_t led=gcData->GetLed();
    Detector::Detector_t detectorType=gcData->GetDetector();

    //check that there is a point present
    if (phVector.empty()){
      MSG("LITuning",Msg::kDebug)
        <<"("<<pulserBox<<":"<<led<<") Input data vectors are empty."
        <<" Skipping this point..."
        <<endl;
      gcData++;
      continue;
    }

    string sPb=Form("%d",pulserBox);
    string sLed=Form("%d",led);
    string sEndBit=", Pulser Box "+sPb+", LED "+sLed;
    MSG("LITuning",Msg::kVerbose) 
      <<"("<<pulserBox<<":"<<led<<") Making graphs..."<<endl;
    gAdcVsPh.push_back(this->TGraphVect(phVector,adcVector));
    gAdcFVsPh.push_back(this->TGraphVect(phVector,adcFVector));
    gPinVsPh.push_back(this->TGraphVect(phVector,pinVector));
    gPin2VsPh.push_back(this->TGraphVect(phVector,pin2Vector));
    gAdcErrVsPh.push_back
      (this->TGraphAsymmErrorsVectEY(phVector,adcVector,adcLowVector,
                                     adcHighVector));
    gAdcFErrVsPh.push_back
      (this->TGraphAsymmErrorsVectEY(phVector,adcFVector,adcFLowVector,
                                     adcFHighVector));
    
    gAdcVsPhIter=gAdcVsPh.end()-1;
    gAdcFVsPhIter=gAdcFVsPh.end()-1;
    gPinVsPhIter=gPinVsPh.end()-1;
    gPin2VsPhIter=gPin2VsPh.end()-1;
    gAdcErrVsPhIter=gAdcErrVsPh.end()-1;
    gAdcFErrVsPhIter=gAdcFErrVsPh.end()-1;

    //do pmt adc graph
    fS="Response of PMTs to LI Pulse Height"+sEndBit; 
    (*gAdcVsPhIter)->SetTitle(fS.c_str());
    (*gAdcVsPhIter)->SetFillColor(0);
    (*gAdcVsPhIter)->SetMinimum(-2);
    if (detectorType==Detector::kFar){
      (*gAdcVsPhIter)->SetMaximum(15000);
    }
    else if (detectorType==Detector::kNear){
      Double_t maxAdc=gcData->GetMaxAdc();
      if (maxAdc>30000){
        (*gAdcVsPhIter)->SetMaximum(maxAdc+0.1*maxAdc);
      }
      else{
        (*gAdcVsPhIter)->SetMaximum(30000);
      }
    }
    else if (detectorType==Detector::kCalDet){
      //don't set max!
    }
    else if (detectorType==Detector::kUnknown){
      //don't set max!
    }
    else {
      MSG("LITuning",Msg::kDebug)
        <<"("<<pulserBox<<":"<<led<<") Detector Type not supported."
        <<"Detector Type=";
      Detector::AsString(detectorType);
      cout<<endl;
    }

    //do far end of strip pmt adc graph
    fS="Response of 'Far End of Strip PMTs' to LI Pulse Height"+
      sEndBit; 
    (*gAdcFVsPhIter)->SetTitle(fS.c_str());
    (*gAdcFVsPhIter)->SetFillColor(0);
    (*gAdcFVsPhIter)->SetMinimum(-2);
    if (detectorType==Detector::kFar){
      (*gAdcFVsPhIter)->SetMaximum(15000);
    }
    else if (detectorType==Detector::kNear){
      Double_t maxAdcF=gcData->GetMaxAdcF();
      if (maxAdcF>30000){
        (*gAdcFVsPhIter)->SetMaximum(maxAdcF+0.1*maxAdcF);
      }
      else{
        (*gAdcFVsPhIter)->SetMaximum(30000);
      }
    }
    else if (detectorType==Detector::kCalDet){
      //don't set max!
    }
    else if (detectorType==Detector::kUnknown){
      //don't set max!
    }
    else {
      MSG("LITuning",Msg::kDebug)
        <<"("<<pulserBox<<":"<<led<<") Detector Type not supported."
        <<"Detector Type=";
      Detector::AsString(detectorType);
      cout<<endl;
    }

    //do pin graph
    fS="Response of PIN to LI Pulse Height"+sEndBit; 
    (*gPinVsPhIter)->SetTitle(fS.c_str());
    (*gPinVsPhIter)->SetFillColor(0);
    (*gPinVsPhIter)->SetMinimum(-2);
    if (detectorType==Detector::kFar){
      (*gPinVsPhIter)->SetMaximum(12500);
    }
    else if (detectorType==Detector::kNear){
      Double_t maxPin=gcData->GetMaxPin();
      if (maxPin>30000){
        (*gPinVsPhIter)->SetMaximum(maxPin+0.1*maxPin);
      }
      else{
        (*gPinVsPhIter)->SetMaximum(30000);
      }
    }
    else if (detectorType==Detector::kCalDet){
      //don't set max!
    }
    else if (detectorType==Detector::kUnknown){
      //don't set max!
    }
    else {
      MSG("LITuning",Msg::kDebug)
        <<"("<<pulserBox<<":"<<led<<") Detector Type not supported."
        <<"Detector Type=";
      Detector::AsString(detectorType);
      cout<<endl;
    }

    //do pin2 graph
    fS="Response of Low Gain PIN to LI Pulse Height"+sEndBit; 
    (*gPin2VsPhIter)->SetTitle(fS.c_str());
    (*gPin2VsPhIter)->SetFillColor(0);
    (*gPin2VsPhIter)->SetMinimum(-2);
    if (detectorType==Detector::kFar){
      (*gPin2VsPhIter)->SetMaximum(12500);
    }
    else if (detectorType==Detector::kNear){
      Double_t maxPin=gcData->GetMaxPin();
      if (maxPin>30000){
        (*gPin2VsPhIter)->SetMaximum(maxPin+0.1*maxPin);
      }
      else{
        (*gPin2VsPhIter)->SetMaximum(30000);
      }
    }
    else if (detectorType==Detector::kCalDet){
      //don't set max!
    }
    else if (detectorType==Detector::kUnknown){
      //don't set max!
    }
    else {
      MSG("LITuning",Msg::kDebug)
        <<"("<<pulserBox<<":"<<led<<") Detector Type not supported."
        <<"Detector Type=";
      Detector::AsString(detectorType);
      cout<<endl;
    }

    //do adcErr graph
    fS="Response of PMTs to LI Pulse Height (Error is 90% range)"+
      sEndBit; 
    (*gAdcErrVsPhIter)->SetTitle(fS.c_str());
    (*gAdcErrVsPhIter)->SetFillColor(0);
    if (detectorType==Detector::kFar){
      (*gAdcErrVsPhIter)->SetMaximum(15000);
    }
    else if (detectorType==Detector::kNear){
      Double_t maxAdc=gcData->GetMaxAdc();
      if (maxAdc>30000){
        (*gAdcErrVsPhIter)->SetMaximum(maxAdc+0.1*maxAdc);
      }
      else{
        (*gAdcErrVsPhIter)->SetMaximum(30000);
      }
    }
    else if (detectorType==Detector::kCalDet){
      //don't set max!
    }
    else if (detectorType==Detector::kUnknown){
      //don't set max!
    }
    else {
      MSG("LITuning",Msg::kDebug)
        <<"("<<pulserBox<<":"<<led<<") Detector Type not supported."
        <<"Detector Type=";
      Detector::AsString(detectorType);
      cout<<endl;
    }

    //do adcErr graph for far end of strip
    fS="Response of PMTs to LI Pulse Height (Error is 90% range)"+
      sEndBit; 
    (*gAdcFErrVsPhIter)->SetTitle(fS.c_str());
    (*gAdcFErrVsPhIter)->SetFillColor(0);
    if (detectorType==Detector::kFar){
      (*gAdcFErrVsPhIter)->SetMaximum(15000);
    }
    else if (detectorType==Detector::kNear){
      Double_t maxAdcF=gcData->GetMaxAdcF();
      if (maxAdcF>30000){
        (*gAdcFErrVsPhIter)->SetMaximum(maxAdcF+0.1*maxAdcF);
      }
      else{
        (*gAdcFErrVsPhIter)->SetMaximum(30000);
      }
    }
    else if (detectorType==Detector::kCalDet){
      //don't set max!
    }
    else if (detectorType==Detector::kUnknown){
      //don't set max!
    }
    else {
      MSG("LITuning",Msg::kDebug)
        <<"("<<pulserBox<<":"<<led<<") Detector Type not supported."
        <<"Detector Type=";
      Detector::AsString(detectorType);
      cout<<endl;
    }

    //print only on debug
    if (MsgService::Instance()->IsActive("LITuning",Msg::kDebug)){
      gcData->PrintAll();
    }

    //iterate the li run that you are looking at
    gcData++;
  }

  MSG("LITuning",Msg::kDebug) 
    <<"Finished GetDataGcGraphs method"<<endl;
  return true;
}

//......................................................................

Bool_t LITuning::GetTunedGcPlots(vector<TH2F*>& hGcAdc,
                                 vector<TH2F*>& hGcPin,
                                 vector<TH2F*>& hGcPh,
                                 Int_t firstLed,Int_t lastLed)
{
  MSG("LITuning",Msg::kInfo) 
    <<"Running GetTunedGcPlots method..."<<endl;

  //check if a gain curve has been calculated
  if (fBestGc.empty()){
    MSG("LITuning",Msg::kWarning)
      <<"There is nothing to print, no tuned values have been"
      <<" computed for a gain curve"<<endl;
    return false;
  }

  //get iterators for the histos
  vector<TH2F*>::iterator hGcAdcIter;
  vector<TH2F*>::iterator hGcPinIter;
  vector<TH2F*>::iterator hGcPhIter;

  vector<LIRun>::iterator tunedRun;
  vector<LIRun>::iterator tunedRunEnd;
 
  //select the tuned gain curve
  tunedRun=fBestGc.begin();
  tunedRunEnd=fBestGc.end();
  
  //variable to hold the last pulser box
  Int_t lastPulserBox=-1;

  //loop over tuned GCs and create the histos
  while (tunedRun!=tunedRunEnd){

    //get iterators to the vectors of li run points
    vector<Double_t> phVector;
    vector<Double_t>::iterator ph;
    phVector=tunedRun->GetPh();
    ph=phVector.begin();
    vector<Double_t> pinVector;
    vector<Double_t>::iterator pin;
    pinVector=tunedRun->GetPin(1);
    pin=pinVector.begin();
    vector<Double_t> adcVector;
    vector<Double_t>::iterator adc;
    adcVector=tunedRun->GetAdc();
    adc=adcVector.begin();

    //get the pulserbox and led 
    Int_t pulserBox=tunedRun->GetPb();
    Int_t led=tunedRun->GetLed();

    //check that there is a point present
    if (phVector.empty()){
      MSG("LITuning",Msg::kInfo)
        <<"("<<pulserBox<<":"<<led<<") Tuned vectors are empty"<<endl;
      tunedRun++;
      continue;
    }

    Int_t xMax=lastLed+1;
    Int_t xMin=firstLed-1;
    Int_t numXbins=xMax-xMin;

    if (pulserBox!=lastPulserBox){
      MSG("LITuning",Msg::kDebug)<<"Next PB..."<<endl;
      string sPb=Form("%d",pulserBox);
      fS="Gain Curve Points in PMT ADCs, Pulser Box "+sPb; 
      hGcAdc.push_back(new TH2F(fS.c_str(),fS.c_str(),
                                numXbins,xMin,xMax,200,0,15000));
      fS="Gain Curve Points in PIN ADCs, Pulser Box "+sPb; 
      hGcPin.push_back(new TH2F(fS.c_str(),fS.c_str(),
                                numXbins,xMin,xMax,200,0,9000));
      fS="Gain Curve PHs, Pulser Box "+sPb; 
      hGcPh.push_back(new TH2F(fS.c_str(),fS.c_str(),
                               numXbins,xMin,xMax,200,0,1050));

      hGcAdcIter=hGcAdc.end()-1;
      hGcPinIter=hGcPin.end()-1;
      hGcPhIter=hGcPh.end()-1;

      (*hGcAdcIter)->GetXaxis()->SetTitle("LED");
      (*hGcAdcIter)->GetXaxis()->CenterTitle();
      (*hGcAdcIter)->GetYaxis()->SetTitle("PMT ADCs");
      (*hGcAdcIter)->GetYaxis()->CenterTitle();
      (*hGcAdcIter)->SetFillColor(0);
      (*hGcAdcIter)->SetBit(TH1::kCanRebin);
      
      (*hGcPinIter)->GetXaxis()->SetTitle("LED");
      (*hGcPinIter)->GetXaxis()->CenterTitle();
      (*hGcPinIter)->GetYaxis()->SetTitle("PIN ADCs");
      (*hGcPinIter)->GetYaxis()->CenterTitle();
      (*hGcPinIter)->SetFillColor(0);
      (*hGcPinIter)->SetBit(TH1::kCanRebin);
      
      (*hGcPhIter)->GetXaxis()->SetTitle("LED");
      (*hGcPhIter)->GetXaxis()->CenterTitle();
      (*hGcPhIter)->GetYaxis()->SetTitle("Pulse Height");
      (*hGcPhIter)->GetYaxis()->CenterTitle();
      (*hGcPhIter)->SetFillColor(0);
      (*hGcPhIter)->SetBit(TH1::kCanRebin);
    }
    //keep a record of last pulser box
    lastPulserBox=pulserBox;

    //print only on debug
    if (MsgService::Instance()->IsActive("LITuning",Msg::kDebug)){
      tunedRun->PrintAll();
    }

    //loop over all the points in the vector (1 or more)
    while(ph!=phVector.end()){

      MSG("LITuning",Msg::kVerbose)
        <<"Looping over tuned points"<<endl;

      (*hGcPinIter)->Fill(led,*pin);
      (*hGcAdcIter)->Fill(led,*adc);
      (*hGcPhIter)->Fill(led,*ph);
      
      adc++;
      pin++;
      ph++;
    }

    //iterate the li run that you are looking at
    tunedRun++;
  }

  MSG("LITuning",Msg::kDebug) 
    <<"GetTunedGcPlots method finished"<<endl;
  return true;
}

//......................................................................

void LITuning::PrintAll()
{
  MSG("LITuning",Msg::kInfo) 
    <<"Running PrintAll method..."<<endl;

  //create iterators
  vector<LIRun>::iterator run;
  vector<LIRun>::iterator runEnd;

  MSG("LITuning",Msg::kDebug) 
    <<"Getting iterators to input gain curve"<<endl;
  //get the begining and end of the vector of gain curves
  run=fGainCurve->begin();
  runEnd=fGainCurve->end();
  MSG("LITuning",Msg::kDebug)
    <<"Got iterators to input gain curve"<<endl;
  
  //loop and print out the data stored in the LIRun objects
  while (run!=runEnd){
    //print out the info stored in the object
    run->PrintAll();
    //iterate the li run that you are printing out
    run++;
  }
  
  MSG("LITuning",Msg::kInfo) 
    <<"PrintAll method finished"<<endl;
}

//......................................................................

---