LIGainCurve.cxx

Go to the documentation of this file.

#define LIGainCurve_cxx

#include <LIGainCurve.h>
#include <TF1.h>
#include <TMath.h>

#include <Validity/VldContext.h>
#include <Plex/PlexStripEndId.h>
#include <Conventions/Detector.h>
#include <Conventions/StripEnd.h>

using namespace std;

static const Short_t MASTERMAX  = 20;
static const Short_t MINDERMAX  = 8;
static const Short_t MENUMAX    = 16;

static const Short_t VARCMAX    = 3;
static const Short_t VMMMAX     = 6;
static const Short_t VFBMAX     = 2;
static const Short_t VAMAX      = 3;
static const Short_t CHANNELMAX = 22;

Bool_t LIGainCurve::fCheckDB = true;


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


//-------------------------------------------------

LIGainCurve::LIGainCurve() {

  fInTree = 0;
  fInFile = 0;

  gErrorIgnoreLevel=1;

  gStyle->SetOptFit(0111);
  gStyle->SetOptStat(0); //iourmen


  //
  // Input File -> A file with a TTree with the GainCurve information.
  // LI-TTree obtained from LIReco.C
  fFileNameIn = "~/DATA/Comp2003/LI_data/LIDataC00071335_0023.root";

  // Output file:
  fFileNameOut = getenv("SRT_PRIVATE_CONTEXT");
  fFileNameOut += "/LISummary/liGainCurveOutput/LIGainCurveOut.ps";

  //
  //Initialised Analysis objects:
  //

  //HG
  fSlopeND_HG_h     = new TH1F("fSlopeND_HG_h","",100,0,0.5);
  fSlopeFD_HG_h     = new TH1F("fSlopeFD_HG_h","",100,0,0.5);
  fInterceptND_HG_h = new TH1F("fInterceptND_HG_h","",100,-75,25);
  fInterceptFD_HG_h = new TH1F("fInterceptFD_HG_h","",100,-75,25);

  fResponseND_HG_p = new TProfile("fResponseND_HG_p","",120,0,12000);
  fResponseFD_HG_p = new TProfile("fResponseFD_HG_p","",120,0,12000);
  fResponseND_HG_h = new TH2F("fResponseND_HG_h","",120,0,12000,100,0,500);
  fResponseFD_HG_h = new TH2F("fResponseFD_HG_h","",120,0,12000,100,0,500);

  //LG
  fSlopeND_LG_h     = new TH1F("fSlopeND_LG_h","",100,0,0.5);
  fSlopeFD_LG_h     = new TH1F("fSlopeFD_LG_h","",100,0,0.5);
  fInterceptND_LG_h = new TH1F("fInterceptND_LG_h","",100,-75,25);
  fInterceptFD_LG_h = new TH1F("fInterceptFD_LG_h","",100,-75,25);

  fResponseND_LG_p = new TProfile("fResponseND_LG_p","",120,0,12000);
  fResponseFD_LG_p = new TProfile("fResponseFD_LG_p","",120,0,12000);
  fResponseND_LG_h = new TH2F("fResponseND_LG_h","",120,0,12000,125,0,250);
  fResponseFD_LG_h = new TH2F("fResponseFD_LG_h","",120,0,12000,125,0,250);

  //Response per LED:
  fResponseND_LG1_p = new TProfile("fResponseND_LG1_p","",40,0,4000);
  fResponseFD_LG1_p = new TProfile("fResponseFD_LG1_p","",40,0,4000);
  fResponseND_LG2_p = new TProfile("fResponseND_LG2_p","",40,0,4000);
  fResponseFD_LG2_p = new TProfile("fResponseFD_LG2_p","",40,0,4000);
  fResponseND_LG3_p = new TProfile("fResponseND_LG3_p","",40,0,4000);
  fResponseFD_LG3_p = new TProfile("fResponseFD_LG3_p","",40,0,4000);
  fResponseND_LG4_p = new TProfile("fResponseND_LG4_p","",40,0,4000);
  fResponseFD_LG4_p = new TProfile("fResponseFD_LG4_p","",40,0,4000);
  fResponseND_LG5_p = new TProfile("fResponseND_LG5_p","",40,0,4000);
  fResponseFD_LG5_p = new TProfile("fResponseFD_LG5_p","",40,0,4000);
  fResponseND_LG6_p = new TProfile("fResponseND_LG6_p","",40,0,4000);
  fResponseFD_LG6_p = new TProfile("fResponseFD_LG6_p","",40,0,4000);

  //Response per LED:
  fResponseND_HG1_p = new TProfile("fResponseND_HG1_p","",50,0,10000);
  fResponseFD_HG1_p = new TProfile("fResponseFD_HG1_p","",50,0,10000);
  fResponseND_HG2_p = new TProfile("fResponseND_HG2_p","",50,0,10000);
  fResponseFD_HG2_p = new TProfile("fResponseFD_HG2_p","",50,0,10000);
  fResponseND_HG3_p = new TProfile("fResponseND_HG3_p","",50,0,10000);
  fResponseFD_HG3_p = new TProfile("fResponseFD_HG3_p","",50,0,10000);
  fResponseND_HG4_p = new TProfile("fResponseND_HG4_p","",50,0,10000);
  fResponseFD_HG4_p = new TProfile("fResponseFD_HG4_p","",50,0,10000);
  fResponseND_HG5_p = new TProfile("fResponseND_HG5_p","",50,0,10000);
  fResponseFD_HG5_p = new TProfile("fResponseFD_HG5_p","",50,0,10000);
  fResponseND_HG6_p = new TProfile("fResponseND_HG6_p","",50,0,10000);
  fResponseFD_HG6_p = new TProfile("fResponseFD_HG6_p","",50,0,10000);



  //Residuals:
  fResidualND_p = new TProfile("fResidualND_p","",50,0,100);
  fResidualFD_p = new TProfile("fResidualFD_p","",50,0,100);
  fResidualFDPMT_p = new TProfile("fResidualFDPMT_p","",50,0,100);

  fResidualFD300_p = new TProfile("fResidualFD300_p","",150,0,300);
  fResidualND300_p = new TProfile("fResidualND300_p","",150,0,300);
  fResidualFD300PMT_p = new TProfile("fResidualFD300PMT_p","",150,0,300);

  fResidualFDADC_p = new TProfile("fResidualFDADC_p","",150,0,15000);
  fResidualNDADC_p = new TProfile("fResidualNDADC_p","",150,0,25000);
  fResidualFDADCPMT_p = new TProfile("fResidualFDADCPMT_p","",150,0,15000);

  fResidualND_h = new TH2F("fResidualND_h","",300,0,300,100,-0.3,0.1);
  fResidualFD_h = new TH2F("fResidualFD_h","",300,0,300,100,-0.3,0.1);
  fResidualFDPMT_h = new TH2F("fResidualFDPMT_h","",300,0,300,100,-0.15,0.1);

  //fResidualN_F_p = new TH1D("fResidualN-F_p","",55,0,1100);
  fResidualN_F_p = new TProfile("fResidualN_F_p","",150,0,300);
  //fResidualN_F_h = new TH2F("fResidualN_F_h","",200,0,200,100,-0.3,0.1);


  fResVsGainND_p = new TProfile("fResVsGainND_p","",25,30,130);
  fResVsGainFD_p = new TProfile("fResVsGainFD_p","",25,30,130);
}

//-------------------------------------------------

LIGainCurve::~LIGainCurve() {

  ;
}

//-------------------------------------------------

void LIGainCurve::ExecuteMe(Int_t nentries) {

  //cout << "Size (empty)" << mapLinResp.size() << endl;
  //mapLinResp[4];
  //LinResp* t = new LinResp();
  //t->addressKey = 40;
  //mapLinResp[4] = t;
  //cout << "Size (one guy!)" << mapLinResp.size() << endl;

  Bool_t success = true;


  MSG("LIGainCurve",Msg::kInfo) << "\tLoading TTree from input file" << endl;

  success = this->LoadTreeLI();

  MSG("LIGainCurve",Msg::kInfo) << "\tAdquiring data from TTree" << endl;

  if(success) success = this->GetData(nentries);

  MSG("LIGainCurve",Msg::kInfo) << "\n\tProcessing the data" << endl;

  if(success) success = this->ProcessData();

  MSG("LIGainCurve",Msg::kInfo) << "\tPlotting results" << endl;

  if(success) this->PlotData();

  if(!success)
    MSG("LIGainCurve",Msg::kFatal) << "There was an error somewhere!" << endl;

  return;
}

//-------------------------------------------------

Bool_t LIGainCurve::LoadTreeLI() {
  
  if(!fInTree) {

    if(!fFileNameIn.empty()) {

      fInFile = TFile::Open( fFileNameIn.c_str() );
      fInTree = static_cast<TTree*>( fInFile->Get("li_tree") );
    }
    else {

      cerr << endl
           << "\tThe Input file name was not given in the ROOT macro 'makeLIGainCurve.C'"
           << endl;
      return false;
    }
    
    if(!fInTree) {

      cerr << "ERROR: TTree is still empty" << endl;
      return false;
    }

    //fInTree->SetBranchAddress("ashtray",&ashtray);
    //fInTree->SetBranchAddress("calibPoint",&calibPoint);
    //fInTree->SetBranchAddress("calibType",&calibType);
    //fInTree->SetBranchAddress("chAdd",&chAdd);
    fInTree->SetBranchAddress("channel",&channel);
    fInTree->SetBranchAddress("chip",&va);
    fInTree->SetBranchAddress("correlatedHit",&correlatedHit);
    fInTree->SetBranchAddress("crate",&crate);
    //fInTree->SetBranchAddress("detectorType",&detectorType);
    fInTree->SetBranchAddress("eastWest",&eastWest);
    fInTree->SetBranchAddress("elecType",&elecType);
    //fInTree->SetBranchAddress("farLed",&farLed);
    //fInTree->SetBranchAddress("farPulserBox",&farPulserBox);
    fInTree->SetBranchAddress("geoAdd",&master);
    //fInTree->SetBranchAddress("inRack",&inRack);
    fInTree->SetBranchAddress("led",&led);
    fInTree->SetBranchAddress("masterCh",&minder);
    fInTree->SetBranchAddress("mean",&mean);
    fInTree->SetBranchAddress("minderCh",&menu);
    fInTree->SetBranchAddress("nearLed",&nearLed);
    fInTree->SetBranchAddress("nearPulserBox",&nearPulserBox);
    fInTree->SetBranchAddress("numEntries",&numEntries);
    //fInTree->SetBranchAddress("numericMuxBox",&numericMuxBox);
    //fInTree->SetBranchAddress("pe",&pe);
    //fInTree->SetBranchAddress("period",&period);
    fInTree->SetBranchAddress("pinGain",&pinGain);
    //fInTree->SetBranchAddress("pinInBox",&pinInBox);
    //fInTree->SetBranchAddress("pixel",&pixel);
    fInTree->SetBranchAddress("plane",&plane);
    fInTree->SetBranchAddress("pulseHeight",&pulseHeight);
    fInTree->SetBranchAddress("pulserBox",&pulserBox);
    //fInTree->SetBranchAddress("pulses",&pulses);
    fInTree->SetBranchAddress("pulseWidth",&pulseWidth);
    //fInTree->SetBranchAddress("rackBay",&rackBay);
    //fInTree->SetBranchAddress("rackLevel",&rackLevel);
    fInTree->SetBranchAddress("readoutType",&readoutType);
    fInTree->SetBranchAddress("rms",&rms);
    //fInTree->SetBranchAddress("runNumber",&runNumber);
    //fInTree->SetBranchAddress("runNumberSub",&runNumberSub);
    //fInTree->SetBranchAddress("runType",&runType);
    fInTree->SetBranchAddress("strip",&strip);
    fInTree->SetBranchAddress("stripEnd",&stripEnd);
    //fInTree->SetBranchAddress("summaryCounter",&summaryCounter);
    fInTree->SetBranchAddress("timestamp",&timestamp);
    fInTree->SetBranchAddress("varc",&varc);
    fInTree->SetBranchAddress("vfb",&vfb);
    fInTree->SetBranchAddress("vmm",&vmm);

  }
  else MSG("LIGainCurve",Msg::kError) << "What's wrong with it? Just continue!" << endl;

  return true;
}

//-------------------------------------------------

Int_t LIGainCurve::EncryptMeKey(const Int_t crate,
                            const Int_t varc,
                            const Int_t vmm,
                            const Int_t vfb,
                            const Int_t va,
                            const Int_t channel) {

  Int_t key = CHANNELMAX*va+channel;
  key      += CHANNELMAX*VAMAX*vfb;
  key      += CHANNELMAX*VAMAX*VFBMAX*vmm;
  key      += CHANNELMAX*VAMAX*VFBMAX*VMMMAX*varc;
  key      += CHANNELMAX*VAMAX*VFBMAX*VMMMAX*VARCMAX*crate;
  
  return key;
}

//-------------------------------------------------

Int_t LIGainCurve:: EncryptMeKey(const Int_t crate,
                             const Int_t master,
                             const Int_t minder,
                             const Int_t menu ) {

  Int_t key = MENUMAX*minder+menu;
  key      += MENUMAX*MINDERMAX*master;
  key      += MENUMAX*MINDERMAX*MASTERMAX*crate;
  key      += 2500; //offset
    
  return key;
}

//-------------------------------------------------

Bool_t LIGainCurve::GetData(Int_t nentries) {

  if (nentries == -1) 
    nentries = Int_t(fInTree->GetEntries());
 
  //LOOP over TTree entries:
  for (Int_t jentry=0; jentry<nentries; jentry++) {

    Int_t ientry = fInTree->GetEntry(jentry);

    if (ientry < 0) {

      cerr  << endl << "ERROR with the TChain entry" << endl;
      break;
    }

    //
    //Collect data into the map of LILinResp's:
    //
    if( correlatedHit == 1           //in agreement with the PLEX
        &&
        pulserBox==nearPulserBox     //Make sure the pulser box is the correcte one
        &&
        led == nearLed ) {           //make sure we use the near LED only


      //
      //PMT data - channels connected to strip-ends:
      //
      if( readoutType == 0x0001 ) {

        PlexStripEndId* seid = 0;
        Int_t key = -1;

        //My local pointer to a LILinResp object 
        //(to be stored in map) 
        LILinResp* myLR = 0;

        //create iter. Initialise to end of map (was =0 but that's
        //not ansi c++ and gcc 4.1 compile fails) This should work ok...
        map<Int_t,LILinResp*>::iterator iter = mapLinResp.end();

        //FD
        if(crate == 0) {
          
          if(!seid) {
            seid = new PlexStripEndId(Detector::kCalDet,plane,strip,StripEnd::kNegative);

            key = this->EncryptMeKey(crate,varc,vmm,vfb,va,channel);

            //Check this channel has already been added or not:
            //if(mapLinResp.size()>0) 
            iter = mapLinResp.find(key);

            //if((*iter).second == 0) {
            if(iter == mapLinResp.end()) {

              myLR = new LILinResp();

              //Give address:
              myLR->typeFEE    = ElecType::kVA;
              myLR->addressKey = key;
              //LED: link to corresponding PIN:
              myLR->led        = led;
            }
            else {

              myLR = (*iter).second;
            }

            if(!myLR) cerr << "This should not happen!" << endl;
          }
          else cerr << "This should not happen!" << endl;
        }
        //ND
        else {
          
          if(!seid) {
            seid = new PlexStripEndId(Detector::kCalDet,plane,strip,StripEnd::kPositive);

            key = this->EncryptMeKey(crate,master,minder,menu);

            //Check this channel has already been added or not:
            //if(mapLinResp.size()>0) 
            iter = mapLinResp.find(key);

            //if((*iter).second == 0) {
            if(iter == mapLinResp.end()) {

              myLR = new LILinResp();

              //Once created initialised:
              //Address:
              myLR->addressKey = key;
              myLR->typeFEE    = ElecType::kQIE;
              //LED: link to corresponding PIN:
              myLR->led        = led;
            }
            else {

              myLR = (*iter).second;
            }

            if(!myLR) cerr << "This should not happen!" << endl;
          }
          else cerr << "This should not happen!" << endl;
        }
        
        //
        // RMS:
        // 
        //  saturation effects on the VA FEE->
        //  RMS seems ok until pulseHeight ~400.
        // 
        //  QIE's RMS looks ok all the time.
        //
        //
                
        //if(key==6218) continue;

        //Get drifted gain of PMT -> PE:
        if(seid) {
          
          Float_t meanCal = 0;
          Float_t mean_e  = rms / sqrt( Float_t(numEntries) );
          
          //cout << jentry << "\t" << key << "\t" << led << "\t" << mapLinResp.size() << endl;

          Float_t gain  = -1;
          Float_t drift = -1;
          
          //If values still not stored into LinResp object -> get them.
          if( (*myLR).gain == -1 && (*myLR).drift == -1 ) {
            
            if(crate == 0 )
              gain = 65;    //(PE/ADC) As measured at CalDet
            if(crate > 0 )
              gain  = 81.;  //(PE/ADC) As measured at CalDet

            //No drift:
            drift = 1.;

            //
            // There was a time that I did all this calculation with the real gains: that was not right!
            //

            //Store into the LinResp objects:
            myLR->gain  = gain;
            myLR->drift = drift;
          }
          else {
            gain  = myLR->gain;
            drift = myLR->drift;
          }

          //Get charge in PEs taken into account the drift:
          if(gain != -1 && drift != -1 )
            meanCal = mean * drift / gain * numEntries/2000.; //Note the total number of entries expected is hardwire here
          else
            MSG("LIGainCurve",Msg::kError) << "No gain/drift numbers available!" << endl;

          //Correct for attenuation of each view:
          //NIY

          

          if(key>-1) {
          
            //Sanity check (>90% of entries):
            if( numEntries > 1800 ) {
              
              myLR->fX.push_back(meanCal);
              myLR->fX_error.push_back(mean_e);
            }
            //Data should not be used -> label with -1 and then removed.
            else {
              
              myLR->fX.push_back(-1);
              myLR->fX_error.push_back(-1);
            }
            mapLinResp[key] = myLR;
          }
        }

        delete seid;
      }
    
      //
      //PIN data:
      //
      if( readoutType == 0x0002 ) {
        
        //NOTE:*******************************************
        //
        //HERE the link between LED <--> PIN is hard-wired 
        //************************************************
        //
        // @ CALDET 2003:
        //
        // LED  |  VMM  |  VFB
        //  1   |   0   |   0
        //  2   |   0   |   1
        //  3   |   1   |   
        //  4   |   2   |   
        //  5   |   4   |   
        //  6   |   3   |   
        //

        MSG("LIGainCurve",Msg::kDebug) << crate << "\t"
                                   << varc << "\t"
                                   << vmm << "\t"
                                   << vfb << "\t"
                                   << va << "\t"
                                   << channel << "\t"
                                   << pinGain << "\t"
                                   << led << "\t"
                                   << nearLed << "\t" 
                                   << numEntries << "\t"
                                   << mean << "\t" 
                                   << endl;
        

        //This swich account for the mapping of the PINs in FEE-space:
        switch(vmm) {
        case 0:
          if(vfb==0) {
            if( numEntries > 1800 ) {
              //cout << "led 1" << endl; 
              if(va==0) LILinResp::fXrefPINLG1.push_back(mean);
              if(va==1) LILinResp::fXrefPINHG1.push_back(mean);
            }
            //Data should not be used:
            else {
              //cout << "led 1" << endl; 
              if(va==0) LILinResp::fXrefPINLG1.push_back(-1);
              if(va==1) LILinResp::fXrefPINHG1.push_back(-1);
            }
          }
          else{
            if( numEntries > 1800 ) {
              //cout << "led 2" << endl; 
              if(va==0) LILinResp::fXrefPINLG2.push_back(mean);
              if(va==1) LILinResp::fXrefPINHG2.push_back(mean);
            }
            //Data should not be used:
            else {
              //cout << "led 2" << endl; 
              if(va==0) LILinResp::fXrefPINLG2.push_back(-1);
              if(va==1) LILinResp::fXrefPINHG2.push_back(-1);
            }
          }
          break;
        case 1:
          if( numEntries > 1800 ) {
              //cout << "led 3" << endl; 
            if(va==0) LILinResp::fXrefPINLG3.push_back(mean);
            if(va==1) LILinResp::fXrefPINHG3.push_back(mean);
          }
          //Data should not be used:
          else {
              //cout << "led 3" << endl; 
            if(va==0) LILinResp::fXrefPINLG3.push_back(-1);
            if(va==1) LILinResp::fXrefPINHG3.push_back(-1);
          }
          break;
        case 2:
          if( numEntries > 1800 ) {
              //cout << "led 4" << endl; 
            if(va==0) LILinResp::fXrefPINLG4.push_back(mean);
            if(va==1) LILinResp::fXrefPINHG4.push_back(mean);
          }
          //Data should not be used:
        else {
              //cout << "led 4" << endl; 
          if(va==0) LILinResp::fXrefPINLG4.push_back(-1);
          if(va==1) LILinResp::fXrefPINHG4.push_back(-1);
        }
          break;
        case 3:
          if( numEntries > 1800 ) {
              //cout << "led 6" << endl; 
            if(va==0) LILinResp::fXrefPINLG6.push_back(mean);
            if(va==1) LILinResp::fXrefPINHG6.push_back(mean);
          }
          //Data should not be used:
          else {
              //cout << "led 6" << endl; 
            if(va==0) LILinResp::fXrefPINLG6.push_back(-1);
            if(va==1) LILinResp::fXrefPINHG6.push_back(-1);
          }
          break;
        case 4:
          if( numEntries > 1800 ) {
              //cout << "led 5" << endl; 
            if(va==0) LILinResp::fXrefPINLG5.push_back(mean);
            if(va==1) LILinResp::fXrefPINHG5.push_back(mean);
          }
          //Data should not be used:
          else {
              //cout << "led 5" << endl; 
            if(va==0) LILinResp::fXrefPINLG5.push_back(-1);
            if(va==1) LILinResp::fXrefPINHG5.push_back(-1);
          }
          break;
        }
      }
      
    } //End of IF => selected response

    else {

      MSG("LILinResp",Msg::kVerbose) << crate << "\t"
                                   << numEntries << "\t"
                                   << correlatedHit << "\t"
                                   << led << "\t"
                                   << nearLed << "\t" 
                                   << channel << "\t"
                                   << endl;
    }

  } //End of LOOP over Input TTREE
  

  MSG("LILinResp",Msg::kInfo) << "\t\t\t\t...Number of LILinResp created: " << mapLinResp.size() << endl;

  //fInFile->Close();

  return true;
}


//-------------------------------------------------


Bool_t LIGainCurve::ProcessData() {

  //Int_t count  = 0;

  Bool_t success = false;
  map<Int_t,LILinResp*>::iterator linRespIter = mapLinResp.begin();

  //
  //LOOP over the map of LinResp
  //
  while( linRespIter != mapLinResp.end() ) {
    
    //Tell LinResp to estimate their linearity
    success = (*linRespIter).second->GetLinearity();

    //(*linRespIter).second->PrintMe();

    if(success) {

      //
      //Fill GLOBAL histograms with slopes & intercepts:
      //

      if((*linRespIter).second->typeFEE == ElecType::kQIE) {

        //HG
        fSlopeND_HG_h->Fill( (*linRespIter).second->slopeHG );
        fInterceptND_HG_h->Fill( (*linRespIter).second->interceptHG );
        //LG
        fSlopeND_LG_h->Fill( (*linRespIter).second->slopeLG );
        fInterceptND_LG_h->Fill( (*linRespIter).second->interceptLG );
      }
      else {

        //HG
        fSlopeFD_HG_h->Fill( (*linRespIter).second->slopeHG );
        fInterceptFD_HG_h->Fill( (*linRespIter).second->interceptHG );
        //LG
        fSlopeFD_LG_h->Fill( (*linRespIter).second->slopeLG );
        fInterceptFD_LG_h->Fill( (*linRespIter).second->interceptLG );
      }

      //
      // Iterate over the vectors in each LinResp to get residuals
      //

      vector<Float_t>* X    = (*linRespIter).second->GetX();
      //vector<Float_t>* X_e  = (*linRespIter).second->GetX_e();
      vector<Float_t>* Res     = (*linRespIter).second->GetResiduals();
      vector<Float_t>* ResPMT  = (*linRespIter).second->GetPMTResiduals();

      vector<Float_t>* XrefPINHG = 0;
      vector<Float_t>* XrefPINLG = 0;

      //Ask each channel for what LED should be injected by.
      //Once the LED is known -> use the LED to find the corresponding Xref.
      switch( (*linRespIter).second->led ) {
      case 1:
        XrefPINHG = ( &LILinResp::fXrefPINHG1 );
        XrefPINLG = ( &LILinResp::fXrefPINLG1 );
        break;
      case 2:
        XrefPINHG = ( &LILinResp::fXrefPINHG2 );
        XrefPINLG = ( &LILinResp::fXrefPINLG2 );
        break;
      case 3:
        XrefPINHG = ( &LILinResp::fXrefPINHG3 );
        XrefPINLG = ( &LILinResp::fXrefPINLG3 );
        break;
      case 4:
        XrefPINHG = ( &LILinResp::fXrefPINHG4 );
        XrefPINLG = ( &LILinResp::fXrefPINLG4 );
        break;
      case 5:
        XrefPINHG = ( &LILinResp::fXrefPINHG5 );
        XrefPINLG = ( &LILinResp::fXrefPINLG5 );
        break;
      case 6:
        XrefPINHG = ( &LILinResp::fXrefPINHG6 );
        XrefPINLG = ( &LILinResp::fXrefPINLG6 );
        break;
      }

      vector<Float_t>::iterator xIter    = X->begin();
      //vector<Float_t>::iterator xeIter   = X_e->begin();
      vector<Float_t>::iterator resIter     = Res->begin();
      vector<Float_t>::iterator resIterPMT  = ResPMT->begin();
      vector<Float_t>::iterator xRefPINHGIter = XrefPINHG->begin();
      vector<Float_t>::iterator xRefPINLGIter = XrefPINLG->begin();


      //
      //LOOP over the responses of each LinResp object:
      //
      while( xIter != X->end() ) {

        if( (*xIter) != -1 && (*xRefPINHGIter) != -1 && (*xRefPINLGIter) != -1 && (*resIter) != -1 ) {

          if( (*linRespIter).second->typeFEE == ElecType::kQIE ) {
            
            if( (*linRespIter).second->led <= 3 ) {

              //Fill GLOBAL histograms with linearity
              //HG:       
              fResponseND_HG_p->Fill( (*xRefPINHGIter),(*xIter) );
              fResponseND_HG_h->Fill( (*xRefPINHGIter),(*xIter) );
              
              //LG:
              fResponseND_LG_p->Fill( (*xRefPINLGIter),(*xIter) );
              fResponseND_LG_h->Fill( (*xRefPINLGIter),(*xIter) );
              
              
              switch( (*linRespIter).second->led ) {
              case 1:
                fResponseND_LG1_p->Fill( (*xRefPINLGIter),(*xIter) );
                fResponseND_HG1_p->Fill( (*xRefPINHGIter),(*xIter) );
                break;
              case 2:
                fResponseND_LG2_p->Fill( (*xRefPINLGIter),(*xIter) );
                fResponseND_HG2_p->Fill( (*xRefPINHGIter),(*xIter) );
                break;
              case 3:
                fResponseND_LG3_p->Fill( (*xRefPINLGIter),(*xIter) );
                fResponseND_HG3_p->Fill( (*xRefPINHGIter),(*xIter) );
                break;
              case 4:
                fResponseND_LG4_p->Fill( (*xRefPINLGIter),(*xIter) );
                fResponseND_HG4_p->Fill( (*xRefPINHGIter),(*xIter) );
                break;
              case 5:
                fResponseND_LG5_p->Fill( (*xRefPINLGIter),(*xIter) );
                fResponseND_HG5_p->Fill( (*xRefPINHGIter),(*xIter) );
                break;
              case 6:
                fResponseND_LG6_p->Fill( (*xRefPINLGIter),(*xIter) );
                fResponseND_HG6_p->Fill( (*xRefPINHGIter),(*xIter) );
                break;
              }
              
              
              //HGPIN was not used
              Float_t xLinearised = (*linRespIter).second->slopeLG * (*xRefPINLGIter) + (*linRespIter).second->interceptLG;

              //Float_t xLinearised = (*linRespIter).second->slopeHG * (*xRefPINHGIter) + (*linRespIter).second->interceptHG;
              

              //Fill GLOBAL histograms with residuals
              if(xLinearised>5 //Low charge structure not necesarily very well describe with GainCurve, there cut, not affect fit for parametrisation
                 && 
                 (*linRespIter).second->slopeLG > 0 
                 &&
                 (*linRespIter).second->slopeLG < 0.3 
                 && 
                 TMath::Abs( (*resIter) ) < 0.2 ) { //For some reason the are some entries screwed up! -> to be investigated
                //fResidualND_p->Fill( (*xIter),(*resIter) );
                //fResidualND_h->Fill( (*xIter),(*resIter) );
                fResidualND_p->Fill( xLinearised,(*resIter) );
                fResidualND_h->Fill( xLinearised,(*resIter) );
                fResidualND300_p->Fill( xLinearised,(*resIter) );
                fResidualNDADC_p->Fill( xLinearised*(*linRespIter).second->gain,(*resIter) );
              }
              
              if(xLinearised > 140 && xLinearised < 160)
                fResVsGainND_p->Fill( (*linRespIter).second->gain, (*resIter) );
            }
          }
          else {
          
            //Fill GLOBAL histograms with linearity
            //HG:
            fResponseFD_HG_p->Fill( (*xRefPINHGIter),(*xIter) );
            fResponseFD_HG_h->Fill( (*xRefPINHGIter),(*xIter) );
          
            //LG:
            fResponseFD_LG_p->Fill( (*xRefPINLGIter),(*xIter) );
            fResponseFD_LG_h->Fill( (*xRefPINLGIter),(*xIter) );
          

            switch( (*linRespIter).second->led ) {
            case 1:
              fResponseFD_LG1_p->Fill( (*xRefPINLGIter),(*xIter) );
              fResponseFD_HG1_p->Fill( (*xRefPINHGIter),(*xIter) );
              break;
            case 2:
              fResponseFD_LG2_p->Fill( (*xRefPINLGIter),(*xIter) );
              fResponseFD_HG2_p->Fill( (*xRefPINHGIter),(*xIter) );
              break;
            case 3:
              fResponseFD_LG3_p->Fill( (*xRefPINLGIter),(*xIter) );
              fResponseFD_HG3_p->Fill( (*xRefPINHGIter),(*xIter) );
              break;
            case 4:
              fResponseFD_LG4_p->Fill( (*xRefPINLGIter),(*xIter) );
              fResponseFD_HG4_p->Fill( (*xRefPINHGIter),(*xIter) );
              break;
            case 5:
              fResponseFD_LG5_p->Fill( (*xRefPINLGIter),(*xIter) );
              fResponseFD_HG5_p->Fill( (*xRefPINHGIter),(*xIter) );
              break;
            case 6:
              fResponseFD_LG6_p->Fill( (*xRefPINLGIter),(*xIter) );
              fResponseFD_HG6_p->Fill( (*xRefPINHGIter),(*xIter) );
              break;
            }


            //HGPIN was not used
            Float_t xLinearised = (*linRespIter).second->slopeLG * (*xRefPINLGIter) + (*linRespIter).second->interceptLG;
            //Float_t xLinearised = (*linRespIter).second->slopeHG * (*xRefPINHGIter) + (*linRespIter).second->interceptHG;


            //Fill GLOBAL histograms with residuals
            if(xLinearised>5 //Low charge structure not necesarily very well describe with GainCurve, there cut, not affect fit for parametrisation
               && 
               (*linRespIter).second->slopeLG > 0 
               &&
               (*linRespIter).second->slopeLG < 0.3 
               && 
               TMath::Abs( (*resIter) ) < 0.4 ) { //For some reason the are some entries screwed up! -> to be investigated

              //fResidualFD_p->Fill( (*xIter),(*resIter) );
              //fResidualFD_h->Fill( (*xIter),(*resIter) );
              fResidualFD_p->Fill( xLinearised,(*resIter) );
              fResidualFD_h->Fill( xLinearised,(*resIter) );
              fResidualFD300_p->Fill( xLinearised,(*resIter) );
              fResidualFDADC_p->Fill( xLinearised*(*linRespIter).second->gain,(*resIter) );

              if(xLinearised <=220) {
                fResidualFDPMT_p->Fill( xLinearised,(*resIterPMT) );
                fResidualFDPMT_h->Fill( xLinearised,(*resIterPMT) );
                fResidualFD300PMT_p->Fill( xLinearised,(*resIterPMT) );
                fResidualFDADCPMT_p->Fill( xLinearised*(*linRespIter).second->gain,(*resIterPMT) );
              }
            }

            if(xLinearised > 140 && xLinearised < 160)
              fResVsGainFD_p->Fill( (*linRespIter).second->gain, (*resIter) );
          }
        }

        xIter++;
        resIter++;
        resIterPMT++;
        xRefPINLGIter++;
        xRefPINHGIter++;
      }
    }

    //count ++;
    linRespIter++;
  }

  //
  //Residuals' difference:
  //
  fResidualN_F_p->Add(fResidualND300_p,       
                      fResidualFD300_p,
                      1,-1);

  //Clear map - not needed any more!
  mapLinResp.clear();

  return true;
}

//-------------------------------------------------

void LIGainCurve::PlotData() {

  //Plot the Summaries:
  Int_t plotCounter = 0;

  fplotFile_c = new TCanvas("plotFile_c","",1200,800);

  string fFileOutStart = fFileNameOut + "(";
  string fFileOutEnd   = fFileNameOut + ")";

  fplotFile_c->cd();
  fplotFile_c->Clear();
  fplotFile_c->Print(fFileOutStart.c_str());

  //
  //Response:
  //

  //ND & HG:
  fplotFile_c->cd();
  fplotFile_c->Divide(2,1);
  fplotFile_c->cd(1);
  gPad->SetLogz(1);
  fResponseND_HG_h->SetXTitle("X_{referece HG}");
  fResponseND_HG_h->SetYTitle("MEAN ND [PE]");
  fResponseND_HG_h->Draw("colz");
  fplotFile_c->cd(2);
  gPad->SetLogz(0);
  fResponseND_HG_p->SetXTitle("X_{referece HG}");
  fResponseND_HG_p->SetYTitle("<MEAN ND [PE]>");
  fResponseND_HG_p->SetLineColor(9);
  fResponseND_HG_p->SetLineWidth(4);
  fResponseND_HG_p->Draw();
  fplotFile_c->cd();
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  plotCounter++;
  fplotFile_c->Clear();


  //ND & LG:
  fplotFile_c->cd();
  fplotFile_c->Divide(2,1);
  fplotFile_c->cd(1);
  gPad->SetLogz(1);
  fResponseND_LG_h->SetXTitle("X_{referece LG}");
  fResponseND_LG_h->SetYTitle("MEAN ND [PE]");
  fResponseND_LG_h->Draw("colz");
  fplotFile_c->cd(2);
  gPad->SetLogz(0);
  fResponseND_LG_p->SetXTitle("X_{referece LG}");
  fResponseND_LG_p->SetYTitle("<MEAN ND [PE]>");
  fResponseND_LG_p->SetLineColor(9);
  fResponseND_LG_p->SetLineWidth(4);
  fResponseND_LG_p->Draw();
  fplotFile_c->cd();
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  plotCounter++;
  fplotFile_c->Clear();



  //FD & HG:
  fplotFile_c->cd();
  fplotFile_c->Divide(2,1);
  fplotFile_c->cd(1);
  gPad->SetLogz(1);
  fResponseFD_HG_h->SetXTitle("X_{referece HG}");
  fResponseFD_HG_h->SetYTitle("MEAN FD [PE]");
  fResponseFD_HG_h->Draw("colz");
  fplotFile_c->cd(2);
  gPad->SetLogz(0);
  fResponseFD_HG_p->SetXTitle("X_{referece HG}");
  fResponseFD_HG_p->SetYTitle("<MEAN FD [PE]>");
  fResponseFD_HG_p->SetLineColor(9);
  fResponseFD_HG_p->SetLineWidth(4);
  fResponseFD_HG_p->Draw();
  fplotFile_c->cd();
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  plotCounter++;
  fplotFile_c->Clear();


  //FD & LG:
  fplotFile_c->cd();
  fplotFile_c->Divide(2,1);
  fplotFile_c->cd(1);
  gPad->SetLogz(1);
  fResponseFD_LG_h->SetXTitle("X_{referece LG}");
  fResponseFD_LG_h->SetYTitle("MEAN FD [PE]");
  fResponseFD_LG_h->Draw("colz");
  fplotFile_c->cd(2);
  gPad->SetLogz(0);
  fResponseFD_LG_p->SetXTitle("X_{referece LG}");
  fResponseFD_LG_p->SetYTitle("<MEAN FD [PE]>");
  fResponseFD_LG_p->SetLineColor(9);
  fResponseFD_LG_p->SetLineWidth(4);
  fResponseFD_LG_p->Draw();
  fplotFile_c->cd();
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  plotCounter++;
  fplotFile_c->Clear();



  //ND & LG:
  fplotFile_c->cd();
  fplotFile_c->Divide(2,3);
  //
  fplotFile_c->cd(1);
  gPad->SetLogz(0);
  fResponseND_LG1_p->SetXTitle("X_{referece LG}(1)");
  fResponseND_LG1_p->SetYTitle("<MEAN ND [PE]>");
  fResponseND_LG1_p->Draw();
  //
  fplotFile_c->cd(2);
  gPad->SetLogz(0);
  fResponseND_LG2_p->SetXTitle("X_{referece LG}(2)");
  fResponseND_LG2_p->SetYTitle("<MEAN ND [PE]>");
  fResponseND_LG2_p->Draw();
  //
  fplotFile_c->cd(3);
  gPad->SetLogz(0);
  fResponseND_LG3_p->SetXTitle("X_{referece LG}(3)");
  fResponseND_LG3_p->SetYTitle("<MEAN ND [PE]>");
  fResponseND_LG3_p->Draw();
  //
  fplotFile_c->cd(4);
  gPad->SetLogz(0);
  fResponseND_LG4_p->SetXTitle("X_{referece LG}(4)");
  fResponseND_LG4_p->SetYTitle("<MEAN ND [PE]>");
  fResponseND_LG4_p->Draw();
  //
  fplotFile_c->cd(5);
  gPad->SetLogz(0);
  fResponseND_LG5_p->SetXTitle("X_{referece LG}(5)");
  fResponseND_LG5_p->SetYTitle("<MEAN ND [PE]>");
  fResponseND_LG5_p->Draw();
  //
  fplotFile_c->cd(6);
  gPad->SetLogz(0);
  fResponseND_LG6_p->SetXTitle("X_{referece LG}(6)");
  fResponseND_LG6_p->SetYTitle("<MEAN ND [PE]>");
  fResponseND_LG6_p->Draw();

  fplotFile_c->cd();
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  plotCounter++;
  fplotFile_c->Clear();



  //FD & LG:
  fplotFile_c->cd();
  fplotFile_c->Divide(2,3);
  //
  fplotFile_c->cd(1);
  gPad->SetLogz(0);
  fResponseFD_LG1_p->SetXTitle("X_{referece LG}(1)");
  fResponseFD_LG1_p->SetYTitle("<MEAN FD [PE]>");
  fResponseFD_LG1_p->Draw();
  //
  fplotFile_c->cd(2);
  gPad->SetLogz(0);
  fResponseFD_LG2_p->SetXTitle("X_{referece LG}(2)");
  fResponseFD_LG2_p->SetYTitle("<MEAN FD [PE]>");
  fResponseFD_LG2_p->Draw();
  //
  fplotFile_c->cd(3);
  gPad->SetLogz(0);
  fResponseFD_LG3_p->SetXTitle("X_{referece LG}(3)");
  fResponseFD_LG3_p->SetYTitle("<MEAN FD [PE]>");
  fResponseFD_LG3_p->Draw();
  //
  fplotFile_c->cd(4);
  gPad->SetLogz(0);
  fResponseFD_LG4_p->SetXTitle("X_{referece LG}(4)");
  fResponseFD_LG4_p->SetYTitle("<MEAN FD [PE]>");
  fResponseFD_LG4_p->Draw();
  //
  fplotFile_c->cd(5);
  gPad->SetLogz(0);
  fResponseFD_LG5_p->SetXTitle("X_{referece LG}(5)");
  fResponseFD_LG5_p->SetYTitle("<MEAN FD [PE]>");
  fResponseFD_LG5_p->Draw();
  //
  fplotFile_c->cd(6);
  gPad->SetLogz(0);
  fResponseFD_LG6_p->SetXTitle("X_{referece LG}(6)");
  fResponseFD_LG6_p->SetYTitle("<MEAN FD [PE]>");
  fResponseFD_LG6_p->Draw();

  fplotFile_c->cd();
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  plotCounter++;
  fplotFile_c->Clear();




  //ND & HG:
  fplotFile_c->cd();
  fplotFile_c->Divide(2,3);
  //
  fplotFile_c->cd(1);
  gPad->SetLogz(0);
  fResponseND_HG1_p->SetXTitle("X_{referece HG}(1)");
  fResponseND_HG1_p->SetYTitle("<MEAN ND [PE]>");
  fResponseND_HG1_p->Draw();
  //
  fplotFile_c->cd(2);
  gPad->SetLogz(0);
  fResponseND_HG2_p->SetXTitle("X_{referece HG}(2)");
  fResponseND_HG2_p->SetYTitle("<MEAN ND [PE]>");
  fResponseND_HG2_p->Draw();
  //
  fplotFile_c->cd(3);
  gPad->SetLogz(0);
  fResponseND_HG3_p->SetXTitle("X_{referece HG}(3)");
  fResponseND_HG3_p->SetYTitle("<MEAN ND [PE]>");
  fResponseND_HG3_p->Draw();
  //
  fplotFile_c->cd(4);
  gPad->SetLogz(0);
  fResponseND_HG4_p->SetXTitle("X_{referece HG}(4)");
  fResponseND_HG4_p->SetYTitle("<MEAN ND [PE]>");
  fResponseND_HG4_p->Draw();
  //
  fplotFile_c->cd(5);
  gPad->SetLogz(0);
  fResponseND_HG5_p->SetXTitle("X_{referece HG}(5)");
  fResponseND_HG5_p->SetYTitle("<MEAN ND [PE]>");
  fResponseND_HG5_p->Draw();
  //
  fplotFile_c->cd(6);
  gPad->SetLogz(0);
  fResponseND_HG6_p->SetXTitle("X_{referece HG}(6)");
  fResponseND_HG6_p->SetYTitle("<MEAN ND [PE]>");
  fResponseND_HG6_p->Draw();

  fplotFile_c->cd();
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  plotCounter++;
  fplotFile_c->Clear();



  //FD & HG:
  fplotFile_c->cd();
  fplotFile_c->Divide(2,3);
  //
  fplotFile_c->cd(1);
  gPad->SetLogz(0);
  fResponseFD_HG1_p->SetXTitle("X_{referece HG}(1)");
  fResponseFD_HG1_p->SetYTitle("<MEAN FD [PE]>");
  fResponseFD_HG1_p->Draw();
  //
  fplotFile_c->cd(2);
  gPad->SetLogz(0);
  fResponseFD_HG2_p->SetXTitle("X_{referece HG}(2)");
  fResponseFD_HG2_p->SetYTitle("<MEAN FD [PE]>");
  fResponseFD_HG2_p->Draw();
  //
  fplotFile_c->cd(3);
  gPad->SetLogz(0);
  fResponseFD_HG3_p->SetXTitle("X_{referece HG}(3)");
  fResponseFD_HG3_p->SetYTitle("<MEAN FD [PE]>");
  fResponseFD_HG3_p->Draw();
  //
  fplotFile_c->cd(4);
  gPad->SetLogz(0);
  fResponseFD_HG4_p->SetXTitle("X_{referece HG}(4)");
  fResponseFD_HG4_p->SetYTitle("<MEAN FD [PE]>");
  fResponseFD_HG4_p->Draw();
  //
  fplotFile_c->cd(5);
  gPad->SetLogz(0);
  fResponseFD_HG5_p->SetXTitle("X_{referece HG}(5)");
  fResponseFD_HG5_p->SetYTitle("<MEAN FD [PE]>");
  fResponseFD_HG5_p->Draw();
  //
  fplotFile_c->cd(6);
  gPad->SetLogz(0);
  fResponseFD_HG6_p->SetXTitle("X_{referece HG}(6)");
  fResponseFD_HG6_p->SetYTitle("<MEAN FD [PE]>");
  fResponseFD_HG6_p->Draw();

  fplotFile_c->cd();
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  plotCounter++;
  fplotFile_c->Clear();


  //
  //Residuals:
  //

  //ND:
  fplotFile_c->cd();
  gPad->SetLogz(1);
  gPad->SetGrid();
  //fResidualND_h->SetMaximum(0.1);
  //fResidualND_h->SetStats(0);
  fResidualND_h->SetXTitle("Linearised Scale (PE)");
  fResidualND_h->SetYTitle("(NonLinear - Linear) / Linear [ND]");
  fResidualND_h->Draw("colz");
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  plotCounter++;
  fplotFile_c->Clear();



  TF1* lin = new TF1("lin","[0]+[1]*x+[2]*x**2+[3]*x**3+[4]*x**4+[5]*x**5");
  //lin->FixParameter(0,0);
  lin->SetParameter(0,0);
  lin->SetParameter(1,0);
  lin->SetParameter(2,0);
  lin->SetParameter(3,0);
  lin->SetParameter(4,0);
  lin->SetParameter(5,0);


  TF1* lin2 = new TF1("lin2","[0]+[1]*x");
  lin2->SetParameter(0,0);
  lin2->SetParameter(1,0);


  //fResidualND300_p->Fit(lin,"Q","",0,300);


  fplotFile_c->cd();
  gPad->SetLogy(0);
  gPad->SetGrid();
  //fResidualND_p->SetStats(0);
  fResidualND_p->SetXTitle("Linearised Scale (PE)");
  fResidualND_p->SetYTitle("(NonLinear - Linear) / Linear [ND]");
  fResidualND_p->SetLineColor(4);
  fResidualND_p->SetMarkerColor(4);
  fResidualND_p->SetLineWidth(4);
  fResidualND_p->Draw();
  lin->Draw("same");
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  fplotFile_c->Print("ND_nonLin-LI.C");
  plotCounter++;
  fplotFile_c->Clear();


  fplotFile_c->cd();
  gPad->SetLogy(0);
  gPad->SetGrid();
  //fResidualND300_p->SetStats(0);
  //fResidualND300_p->Fit(lin,"Q","",0,300);
  fResidualND300_p->SetXTitle("Linearised Scale (PE)");
  fResidualND300_p->SetYTitle("(NonLinear - Linear) / Linear [ND]");
  fResidualND300_p->SetLineColor(4);
  fResidualND300_p->SetMarkerColor(4);
  fResidualND300_p->SetLineWidth(4);
  fResidualND300_p->Draw();
  //lin1->Draw("same");
  //lin2->Draw("same");
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  fplotFile_c->Print("ND300_nonLin-LI.C");
  plotCounter++;
  fplotFile_c->Clear();


  fplotFile_c->cd();
  gPad->SetLogy(0);
  gPad->SetGrid();
  //fResidualNDADC_p->SetStats(0);
  fResidualNDADC_p->Fit(lin2,"Q","",2500,25000);
  fResidualNDADC_p->Fit(lin,"Q","",0,5000);
  fResidualNDADC_p->SetXTitle("Linearised Scale (ADC)");
  fResidualNDADC_p->SetYTitle("(NonLinear - Linear) / Linear [ND]");
  fResidualNDADC_p->SetLineColor(4);
  fResidualNDADC_p->SetMarkerColor(4);
  fResidualNDADC_p->SetLineWidth(4);
  fResidualNDADC_p->Draw();
  lin->Draw("same");
  lin2->Draw("same");
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  fplotFile_c->Print("NDADC_nonLin-LI.C");
  plotCounter++;
  fplotFile_c->Clear();


  delete lin;
  lin=0;

  //FD:
  fplotFile_c->cd();
  gPad->SetLogz(1);
  gPad->SetGrid();
  //fResidualFD_h->SetMaximum(0.1);
  fResidualFD_h->SetStats(0);
  fResidualFD_h->SetXTitle("Linearised Scale (PE)");
  fResidualFD_h->SetYTitle("(NonLinear - Linear) / Linear [FD]");
  fResidualFD_h->Draw("colz");
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  plotCounter++;
  fplotFile_c->Clear();


  //FD M16 parametrisation (Mike's)
  TF1* g = new TF1("f","(5e-2-1.4e-3*x+5.6e-6*x*x-7.5e-9*x*x*x)",0,300);
  g->SetLineColor(6);



  if(!lin) lin = new TF1("lin","[0]+[1]*x+[2]*x**2+[3]*x**3+[4]*x**4+[5]*x**5");
  lin->FixParameter(0,0);
  lin->SetParameter(1,0);
  lin->SetParameter(2,0);
  lin->SetParameter(3,0);
  lin->SetParameter(4,0);
  lin->SetParameter(5,0);


  fplotFile_c->cd();
  gPad->SetLogy(0);
  gPad->SetGrid();
  //fResidualFD_p->SetStats(0);
  //fResidualFD_h->SetMinimum(-0.12);
  //fResidualFD_p->Fit(lin,"Q","",0,300);
  fResidualFD_p->SetXTitle("Linearised Scale (PE)");
  fResidualFD_p->SetYTitle("(NonLinear - Linear) / Linear [FD]");
  fResidualFD_p->SetLineColor(4);
  fResidualFD_p->SetMarkerColor(4);
  fResidualFD_p->SetLineWidth(4);
  fResidualFD_p->Draw();
  //g->Draw("same");
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  fplotFile_c->Print("FD_nonLin-LI.C");
  plotCounter++;
  fplotFile_c->Clear();



  fplotFile_c->cd();
  gPad->SetLogy(0);
  gPad->SetGrid();
  //fResidualFD300_p->SetStats(0);
  //fResidualFD300_h->SetMinimum(-0.12);
  //fResidualFD300_p->Fit(lin,"Q","",0,300);
  fResidualFD300_p->SetXTitle("Linearised Scale (PE)");
  fResidualFD300_p->SetYTitle("(NonLinear - Linear) / Linear [FD]");
  fResidualFD300_p->SetLineColor(4);
  fResidualFD300_p->SetMarkerColor(4);
  fResidualFD300_p->SetLineWidth(4);
  fResidualFD300_p->Draw();
  //g->Draw("same");
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  fplotFile_c->Print("FD300_nonLin-LI.C");
  plotCounter++;
  fplotFile_c->Clear();



  fplotFile_c->cd();
  gPad->SetLogy(0);
  gPad->SetGrid();
  //fResidualFDADC_p->SetStats(0);
  //fResidualFDADC_h->SetMinimum(-0.12);
  fResidualFDADC_p->Fit(lin,"Q","",0,15000);
  fResidualFDADC_p->SetXTitle("Linearised Scale (ADC)");
  fResidualFDADC_p->SetYTitle("(NonLinear - Linear) / Linear [FD]");
  fResidualFDADC_p->SetLineColor(4);
  fResidualFDADC_p->SetMarkerColor(4);
  fResidualFDADC_p->SetLineWidth(4);
  fResidualFDADC_p->Draw();
  //g->Draw("same");
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  fplotFile_c->Print("FDADC_nonLin-LI.C");
  plotCounter++;
  fplotFile_c->Clear();



  delete lin;
  lin = 0;

  //
  //
  //PMT non linearity: M16
  //
  //


  //FD:
  fplotFile_c->cd();
  gPad->SetLogz(1);
  gPad->SetGrid();
  fResidualFDPMT_h->SetStats(0);
  fResidualFDPMT_h->SetXTitle("Linearised Scale (PE)");
  fResidualFDPMT_h->SetYTitle("(NonLinear - Linear) / Linear [FD]");
  //fResidualFDPMT_h->SetMinimum(-0.15);
  fResidualFDPMT_h->Draw("colz");
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  plotCounter++;
  fplotFile_c->Clear();


  //
  //FD M16 parametrisation (Mike's)
  //TF1* g = new TF1("f","(5e-2-1.4e-3*x+5.6e-6*x*x-7.5e-9*x*x*x)",0,300);
  //g->SetLineColor(6);



  if(!lin) lin = new TF1("lin","[0]+[1]*x"); //+[2]*x**2+[3]*x**3+[4]*x**4+[5]*x**5");
  //lin->FixParameter(0,0);
  lin->SetParameter(0,0);
  lin->SetParameter(1,0);
  //lin->SetParameter(2,0);
  //lin->SetParameter(3,0);
  //lin->SetParameter(4,0);
  //lin->SetParameter(5,0);


  fplotFile_c->cd();
  gPad->SetLogy(0);
  gPad->SetGrid();
  //fResidualFDPMT_p->SetStats(0);
  //fResidualFDPMT_p->Fit(lin,"Q","",55,200);
  fResidualFDPMT_p->SetXTitle("Linearised Scale (PE)");
  fResidualFDPMT_p->SetYTitle("(NonLinear - Linear) / Linear [FD]");
  fResidualFDPMT_p->SetLineColor(4);
  fResidualFDPMT_p->SetMarkerColor(4);
  fResidualFDPMT_p->SetLineWidth(4);
  fResidualFDPMT_p->Draw();
  //g->Draw("same");
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  fplotFile_c->Print("M16_nonLin-LI.C");
  plotCounter++;
  fplotFile_c->Clear();



  fplotFile_c->cd();
  gPad->SetLogy(0);
  gPad->SetGrid();
  //fResidualFD300PMT_p->SetStats(0);
  //fResidualFD300PMT_p->Fit(lin,"Q","",55,200);
  fResidualFD300PMT_p->SetXTitle("Linearised Scale (PE)");
  fResidualFD300PMT_p->SetYTitle("(NonLinear - Linear) / Linear [FD]");
  fResidualFD300PMT_p->SetLineColor(4);
  fResidualFD300PMT_p->SetMarkerColor(4);
  fResidualFD300PMT_p->SetLineWidth(4);
  fResidualFD300PMT_p->Draw();
  //g->Draw("same");
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  fplotFile_c->Print("M16-300_nonLin-LI.C");
  plotCounter++;
  fplotFile_c->Clear();


  fplotFile_c->cd();
  gPad->SetLogy(0);
  gPad->SetGrid();
  //fResidualFDADCPMT_p->SetStats(0);
  fResidualFDADCPMT_p->Fit(lin,"Q","",3600,13000);
  fResidualFDADCPMT_p->SetXTitle("Linearised Scale (ADC)");
  fResidualFDADCPMT_p->SetYTitle("(NonLinear - Linear) / Linear [FD]");
  fResidualFDADCPMT_p->SetLineColor(4);
  fResidualFDADCPMT_p->SetMarkerColor(4);
  fResidualFDADCPMT_p->SetLineWidth(4);
  fResidualFDADCPMT_p->Draw();
  //g->Draw("same");
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  fplotFile_c->Print("M16-ADC_nonLin-LI.C");
  plotCounter++;
  fplotFile_c->Clear();



  delete lin;
  lin = 0;


  //
  //Both in the same plot:
  //


  fplotFile_c->cd();
  gPad->SetLogy(0);
  gPad->SetGrid();
  fResidualND300_p->SetStats(0);
  fResidualFD300PMT_p->SetStats(0);
  fResidualND300_p->Draw();
  fResidualFD300PMT_p->Draw("same");
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  fplotFile_c->Print("PMT-300_nonLin-LI.C");
  plotCounter++;
  fplotFile_c->Clear();



  fplotFile_c->cd();
  gPad->SetLogy(0);
  gPad->SetGrid();
  fResidualND300_p->SetStats(0);
  fResidualFD300_p->SetStats(0);
  fResidualND300_p->Draw();
  fResidualFD300_p->Draw("same");
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  plotCounter++;
  fplotFile_c->Clear();




  //
  //Difference between ND and FD:
  //

  if(!lin) lin = new TF1("lin","[0]+[1]*x+[2]*x**2+[3]*x**3+[4]*x**4+[5]*x**5");
  lin->SetParameter(0,0);
  lin->SetParameter(1,0);
  lin->SetParameter(2,0);
  lin->SetParameter(3,0);
  lin->SetParameter(4,0);
  lin->SetParameter(5,0);


  fplotFile_c->cd();
  gPad->SetLogy(0);
  gPad->SetGrid();
  fResidualN_F_p->SetXTitle("Linearised Scale (PE)");
  fResidualN_F_p->SetYTitle("<ND - FD>");
  fResidualN_F_p->SetAxisRange(0,150,"X");
  fResidualN_F_p->Fit(lin,"Q","",0,150);
  fResidualN_F_p->Draw();
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  fplotFile_c->Print("RelNonLin-LI.C");
  plotCounter++;
  fplotFile_c->Clear();
  //TF1* g = new TF1("f","-4.*(-3.63e-4*x*4+1.67e-5*pow(x*4,2)-2.27e-7*pow(x*4,3)+1.30e-9*pow(x*4,4)-3.02e-12*pow(x*4,5))",0,200);



  fplotFile_c->cd();
  gPad->SetLogy(0);
  gPad->SetGrid();
  fResVsGainFD_p->SetXTitle("Gain (PE/ADC)");
  //fResVsGainND_p->SetXTitle("");
  fResVsGainFD_p->SetYTitle("(NonLinear - Linear) / Linear [@ 150PE]");
  fResVsGainFD_p->SetMinimum(-0.1);
  fResVsGainFD_p->SetLineColor(2);
  fResVsGainND_p->SetLineColor(4);
  fResVsGainFD_p->SetMarkerColor(2);
  fResVsGainND_p->SetMarkerColor(4);
  fResVsGainFD_p->Draw();
  fResVsGainND_p->Draw("same");
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  plotCounter++;
  fplotFile_c->Clear();




  //
  //Slopes & Intercepts:
  //


  //ND & HG:
  fplotFile_c->cd();
  gPad->SetLogy(0);
  fSlopeND_HG_h->SetXTitle("Response Slope ND");
  fSlopeND_HG_h->Draw();
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  plotCounter++;
  fplotFile_c->Clear();

  fplotFile_c->cd();
  gPad->SetLogy(0);
  fInterceptND_HG_h->SetXTitle("Response Intercept ND");
  fInterceptND_HG_h->Draw();
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  plotCounter++;
  fplotFile_c->Clear();


  //ND & LG:
  fplotFile_c->cd();
  gPad->SetLogy(0);
  fSlopeND_LG_h->SetXTitle("Response Slope ND");
  fSlopeND_LG_h->Draw();
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  plotCounter++;
  fplotFile_c->Clear();

  fplotFile_c->cd();
  gPad->SetLogy(0);
  fInterceptND_LG_h->SetXTitle("Response Intercept ND");
  fInterceptND_LG_h->Draw();
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  plotCounter++;
  fplotFile_c->Clear();


  //FD & HG:
  fplotFile_c->cd();
  gPad->SetLogy(0);
  fSlopeFD_HG_h->SetXTitle("Response Slope FD");
  fSlopeFD_HG_h->Draw();
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  plotCounter++;
  fplotFile_c->Clear();

  fplotFile_c->cd();
  gPad->SetLogy(0);
  fInterceptFD_HG_h->SetXTitle("Response Intercept FD");
  fInterceptFD_HG_h->Draw();
  fplotFile_c->Update();
  fplotFile_c->Print(fFileNameOut.c_str());
  plotCounter++;
  fplotFile_c->Clear();


  //FD & LG:
  fplotFile_c->cd();
  gPad->SetLogy(0);
  fSlopeFD_LG_h->SetXTitle("Response Slope FD");
  fSlopeFD_LG_h->Draw();
  fplotFile_c->Update();
  //fplotFile_c->Print(fFileNameOut.c_str());
  plotCounter++;
  fplotFile_c->Clear();

  fplotFile_c->cd();
  gPad->SetLogy(0);
  fInterceptFD_LG_h->SetXTitle("Response Intercept FD");
  fInterceptFD_LG_h->Draw();
  fplotFile_c->Update();
  fplotFile_c->Print(fFileOutEnd.c_str());
  //fplotFile_c->Print(fFileNameOut.c_str());
  plotCounter++;
  fplotFile_c->Clear();


  //Close the Booklet:
  fplotFile_c->cd();
  fplotFile_c->Clear();
  fplotFile_c->Update();
  //fplotFile_c->Print(fFileOutEnd.c_str());

  delete fplotFile_c;
 
  return; 
}

//-------------------------------------------------

//EOF

---