CDAnalysis.cxx

Go to the documentation of this file.

// Program name: CDAnalysis.cxx
//                                                                     
// Package: CalDetTracker   
//                                                                    
// Coded by Jeff Hartnell Jul/2003 onwards
//                                                                    
// Purpose: To make plots and analyse the Tracker*.root files
//                                                                    
// Contact: jeffrey.hartnell@physics.ox.ac.uk                         

#include <fstream>
#include <cmath>
#include <cassert>

#include "TCanvas.h"
#include "TChain.h"
#include "TError.h"
#include "TFile.h"
#include "TF1.h"
#include "TGraphErrors.h"
#include "TH2F.h"
#include "TLegend.h"
#include "TProfile.h"
#include "TProfile2D.h"
#include "TStyle.h"
#include "TTree.h"
#include "TVector3.h"
#include "TClonesArray.h"

#include "CalDetPID/CalDetParticleType.h"
#include "Calibrator/Calibrator.h"
#include "Digitization/DigiSignal.h"
#include "MessageService/MsgService.h"
#include "MuELoss/MuEnergyLoss.h"
#include "Conventions/Munits.h"
#include "Conventions/StripEnd.h"
#include "Validity/VldContext.h"

#include "CalDetTracker/CDAnalysis.h"
#include "CalDetTracker/CDSimpleMC.h"
#include "CalDetTracker/CDTruthHitInfo.h"
#include "CalDetTracker/CDPIDInfo.h"
#include "CalDetTracker/CDTrackedHitInfo.h"
#include "CalDetTracker/CDTrackInfo.h"
#include "CalDetTracker/CDTrackerOptions.h"
#include "CalDetTracker/CDXTalkHitInfo.h"

using std::endl;
using std::cout;
using std::map;
using std::vector;

string CDAnalysis::fInputFileName="";

ClassImp(CDAnalysis)

CVSID("$Id: CDAnalysis.cxx,v 1.74 2007/07/03 16:01:27 hartnell Exp $");

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

CDAnalysis::CDAnalysis()
{
  MSG("CDAnalysis",Msg::kInfo)
      <<"Running CDAnalysis Constructor..."<<endl;

  //initialise variables
  fOutFile=0;//initialise output file pointer
  fTrackerTree=0;
  fOptTree=0;
  fTrkOpt=0;//tracking options
  fPIDInfo=0;
  fTrackInfo=0;
  fTrkHitInfo=0;
  fUnTrkHitInfo=0;
  fXTalkHits=0;
  fTruthHitInfo=0;
  
  fSnarl=-1;
  fSec=-1;
  fNumDeadChips=-1;
  fTemperature=-1;
  fBeamMomentum=0;
  fSimFlag=SimFlag::kUnknown;
  fRunNumber=0;
  fSubRun=0;
  fStartSnarl=0;
  fEndSnarl=0;
  fStartTime=0;
  fEndTime=0;
  fMainParticleEnergy=0;
  fEvents=-1;

  this->InitialiseLoopVariables();
  this->InitialiseHitInfoVariables();
  this->InitialisePidVariables();

  //set a default value of 1 so it stays as sigcor
  fSigCorsPerMip=1;
  fSigCorPerMEU=1;
  fScEnDepPerMEU=1;
  fRatioScEnToBeamEn=1;
  
  fMeuVsEvLenCutM=0;
  fSigCorVsDistM=0;
  fPeakLastPlane=0;
  fPeakPlaneSigCors=0;
  fMyPiMuCounter=0;
  fPidPiMuCounter=0;
  fMyElecCounter=0;
  fPidElecCounter=0;
  fBothPiMuCounter=0;
  fBothElecCounter=0;
  fAvTemperature=0;

  fS="";//general purpose string
  fDoReCalibration=false;//don't do recalibration by default

  //all data members are initialised by this point...

  this->MakeChain();
  this->SetRootFileObjects();

  //get the main particle energy
  fMainParticleEnergy=this->GetGreatestMainParticleEnergy();
  if (fMainParticleEnergy>0){
    MSG("CDAnalysis",Msg::kInfo)
      <<"Found main energy="<<fMainParticleEnergy<<endl;
    MSG("CDAnalysis",Msg::kInfo)
      <<"Beam Momentum from file="<<fBeamMomentum<<endl;
  }

  //setup some root stuff
  //include the under and overflow counts
  gStyle->SetOptStat(1111111);
  gStyle->SetOptFit(1111);
  gStyle->SetPalette(1,(Int_t*)0);

  //do these last in the constructor
  //set the values appropriately if possible
  this->SetSigCorPerMEU();
  this->SetScEnDepPerMEU();
  this->SetRatioScEnToBeamEn();

  MSG("CDAnalysis",Msg::kInfo)
      <<"Finished CDAnalysis Constructor"<<endl;
}

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

CDAnalysis::~CDAnalysis()
{
  MSG("CDAnalysis",Msg::kDebug)
      <<"Running CDAnalysis Destructor..."<<endl;
  
  if (fOutFile){
    fOutFile->Close();
  }

  MSG("CDAnalysis",Msg::kDebug)
      <<"Finished CDAnalysis Destructor"<<endl;
}

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

void CDAnalysis::InitialiseLoopVariables()
{
  MSG("CDAnalysis",Msg::kDebug)<<"Initialising loop variables..."<<endl;

  fFirstPlane=69;
  fLastPlane=-1;
  fLastPlaneOdd=-1;
  fLastPlaneEven=-1;

  fNumPlanesHitAll.clear();
  fNumPlanesHit25.clear();
  fNumPlanesHitTrk.clear();
  fNumPlanesHitPeCut.clear();
  fNumPlanesHitPeCut10.clear();

  fNumHitsPerPlane=-1;
  fNumHitsPerPlane25=-1;
  fNumHitsPerPlanePeCut10=-1;
  fStripsFromCentrePeCut=-1;

  //variables to calc av strip for ps muon cut
  fAvStrip=-1;
  fStCount=0;
  fAvStrip1=-1;
  fStCount1=0;
  fAvStrip2=-1;
  fStCount2=0;

  fLowScint_Tof=1e9;
  fUpScint_Tof=-1e9;

  MSG("CDAnalysis",Msg::kDebug)<<"Initialisation complete"<<endl;
}

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

void CDAnalysis::InitialisePidVariables()
{
  //these are the SI variables
  fTriggerPMT=-1;
  fFafErr=-1;
  fSparseErr=-1;
  fTrigSource=-1;

  fKovADC1=-1;
  fKovTimeStamp1=-1;
  fKovADC2=-1;
  fKovTimeStamp2=-1;
  fKovADC3=-1;
  fKovTimeStamp3=-1;

  fSnarlTimeFrame=-1;
  fSnarlMinTimeStamp=0;
  fSnarlMaxTimeStamp=0;

  fTofTDC0=-1;
  fTofTDC1=-1;
  fTofTDC2=-1;
  fTofADC0=-1;
  fTofADC1=-1;
  fTofADC2=-1;
  fTofADCTimeStamp0=-1;
  fTofADCTimeStamp1=-1;
  fTofADCTimeStamp2=-1;
  fTofTimeStamp=-1;//don't use this, use fTofADCTimeStamp1==TTAG

  fTickSinceLast=-1;

  //these are the PID variables
  fNoOverlap=kFALSE;
  fInCERTime=kFALSE;
  fPIDType=0;
  fNoOverlapBits=0; 
  fInCERTimeBits=0; 
  fOLChi2=0;
}

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

void CDAnalysis::InitialiseHitInfoVariables()
{
  fTime=0;

  fPlane=-1;
  fResultOdd=-1;
  fResultEven=-1;
  fStrip=-1;
  fStripend=-1;
  fTransPos=-1;

  fO1=false;
  fO2=false;
  fE1=false;
  fE2=false;

  fChargeAdc=-1;
  fChargeMip=-1;
  fChargeSigCor=-1;
  fChargeSigLin=-1;
  fChargePe=-1;
  fGain=-1;
}

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

void CDAnalysis::WriteOutHistos() const
{
  //write out the histos to the file, if it's open
  if (fOutFile){
    if (fOutFile->IsWritable()){
      fOutFile->cd();

      //create a tree and set branches      
      TTree runInfo("runInfo","runInfo");
      //have to create local variables here
      //I think it's something to do with the fact that they are
      //read in on other trees in this class as well???
      //No, actually not because fSigCorsPerMip is not in a tree
      //don't understand this crazyness!
      Int_t lSimFlag=fSimFlag;
      runInfo.Branch("SimFlag",&lSimFlag,"SimFlag/I");
      Float_t lBeamMomentum=fBeamMomentum;
      runInfo.Branch("BeamMomentum",&lBeamMomentum,"BeamMomentum/F");
      Int_t lRunNumber=fRunNumber;
      runInfo.Branch("RunNumber",&lRunNumber,"RunNumber/I");
      Int_t lSubRun=fSubRun;
      runInfo.Branch("SubRun",&lSubRun,"SubRun/I");
      Int_t lRunPeriod=this->RunNumber2RunPeriod(fRunNumber);
      runInfo.Branch("RunPeriod",&lRunPeriod,"RunPeriod/I");
      Int_t lStartTime=fStartTime;
      runInfo.Branch("StartTime",&lStartTime,"StartTime/I");
      Int_t lEndTime=fEndTime;
      runInfo.Branch("EndTime",&lEndTime,"EndTime/I");
      Int_t lStartSnarl=fStartSnarl;
      runInfo.Branch("StartSnarl",&lStartSnarl,"StartSnarl/I");
      Int_t lEndSnarl=fEndSnarl;
      runInfo.Branch("EndSnarl",&lEndSnarl,"EndSnarl/I");
      Float_t lAvTemperature=fAvTemperature;
      runInfo.Branch("AvTemperature",&lAvTemperature,"AvTemperature/F");

      Float_t lSigCorsPerMip=fSigCorsPerMip;
      runInfo.Branch("SigCorsPerMip",&lSigCorsPerMip,"SigCorsPerMip/F");
      Float_t lScEnDepPerMEU=fScEnDepPerMEU;
      runInfo.Branch("ScEnDepPerMEU",&lScEnDepPerMEU,"ScEnDepPerMEU/F");
      Float_t lRatioScEnToBeamEn=fRatioScEnToBeamEn;
      runInfo.Branch("RatioScEnToBeamEn",&lRatioScEnToBeamEn,
                     "RatioScEnToBeamEn/F");
      Float_t lPeakLastPlane=fPeakLastPlane;
      runInfo.Branch("PeakLastPlane",&lPeakLastPlane,
                     "PeakLastPlane/F");
      Float_t lPeakPlaneSigCors=fPeakPlaneSigCors;
      runInfo.Branch("PeakPlaneSigCors",&lPeakPlaneSigCors,
                     "PeakPlaneSigCors/F");
      Float_t lSigCorVsDistM=fSigCorVsDistM;
      runInfo.Branch("SigCorVsDistM",&lSigCorVsDistM,
                     "SigCorVsDistM/F");
      Float_t lMeuVsEvLenCutM=fMeuVsEvLenCutM;
      runInfo.Branch("MeuVsEvLenCutM",&lMeuVsEvLenCutM,
                     "MeuVsEvLenCutM/F");
      Int_t lEvLenCutMin=-1;
      Int_t lEvLenCutMax=-1;
      this->GetEvLenMinMax(lEvLenCutMin,lEvLenCutMax);
      runInfo.Branch("EvLenCutMin",&lEvLenCutMin,"EvLenCutMin/I");
      runInfo.Branch("EvLenCutMax",&lEvLenCutMax,"EvLenCutMax/I");

      //pid counters
      Int_t lMyPiMuCounter=fMyPiMuCounter;
      runInfo.Branch("MyPiMuCounter",&lMyPiMuCounter,"MyPiMuCounter/I");
      Int_t lPidPiMuCounter=fPidPiMuCounter;
      runInfo.Branch("PidPiMuCounter",&lPidPiMuCounter,
                     "PidPiMuCounter/I");
      Int_t lMyElecCounter=fMyElecCounter;
      runInfo.Branch("MyElecCounter",&lMyElecCounter,"MyElecCounter/I");
      Int_t lPidElecCounter=fPidElecCounter;
      runInfo.Branch("PidElecCounter",&lPidElecCounter,
                     "PidElecCounter/I");
      Int_t lBothPiMuCounter=fBothPiMuCounter;
      runInfo.Branch("BothPiMuCounter",&lBothPiMuCounter,
                     "BothPiMuCounter/I");
      Int_t lBothElecCounter=fBothElecCounter;
      runInfo.Branch("BothElecCounter",&lBothElecCounter,
                     "BothElecCounter/I");

      //fill and write      
      runInfo.Fill();
      runInfo.Write();

      MSG("CDAnalysis",Msg::kInfo)
        <<"Writing histos to: "<<fOutFile->GetName()<<" ..."<<endl;
      fOutFile->Write();
      //fOutFile->Close();//this makes histos disappear from
      //the canvases for some strange reason!
    }
    else {
      MSG("CDAnalysis",Msg::kWarning)
        <<"File not writable!"<<endl;
    }
  }
}

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

void CDAnalysis::InputFileName(string f)
{
  MSG("CDAnalysis",Msg::kInfo)
    <<"Running with input file name="<<f<<endl;
  fInputFileName=f;
}

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

vector<string> CDAnalysis::MakeFileList()
{

  vector<string> fileList;
  
  if (fInputFileName!=""){
    ifstream inputFile(fInputFileName.c_str());

    //check if file exists
    if (inputFile){
      //variables to hold input from file
      string file="";
   
      //read in from the text file and fill objects
      while(inputFile>>file) {
        MSG("CDAnalysis",Msg::kDebug)
          <<"Found input file name="<<file<<endl;
        fileList.push_back(file);
      }
      MSG("CDAnalysis",Msg::kDebug)
        <<"Files names found in txt file="<<fileList.size()<<endl;
    }
    else{
      MSG("CDAnalysis",Msg::kFatal)
        <<endl<<endl
        <<"***********************************************************"
        <<endl<<"Input txt file of file names does not exist!"<<endl
        <<"InputFileName="<<fInputFileName<<endl
        <<"***********************************************************"
        <<endl<<endl<<"Program will exit here"<<endl;
      exit(0);
    }
  }
  //return the fileList if files were found
  if (fileList.size()>0) return fileList;
  
  //Check the env variable to find files
  char* envVariable=getenv("CDDATA");
  if (envVariable==NULL){
    MSG("CDAnalysis",Msg::kFatal)
      <<endl<<endl
      <<"*************************************************************"
      <<endl<<"Environmental variable CDDATA not set!"<<endl
      <<"Please set CDDATA to the directory containing the"
      <<" Tracker*.root files"<<endl
      <<"Note: If more than one file is found they will be"
      <<" concatenated and treated as one"<<endl
      <<"*************************************************************"
      <<endl<<endl<<"Program will exit here"<<endl;
    exit(0);
  }
  string sEnv=envVariable;
  MSG("CDAnalysis",Msg::kInfo)
    <<"Looking for Tracker*.root files using the env variable"<<endl
    <<"CDDATA="<<sEnv<<endl;
  sEnv+="/Tracker*.root";
  fileList.push_back(sEnv);

  return fileList;
}

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

void CDAnalysis::MakeChain()
{
  //get the files to open
  vector<string> fileList=this->MakeFileList();

  //create a chain with TrackerTree
  fTrackerTree=new TChain("TrackerTree");
  fOptTree=new TChain("TrackerOptions");

  Int_t nf=0;
  //add the files to the chain
  for (vector<string>::iterator file=fileList.begin();
       file!=fileList.end();++file){
    
    //test if file already exists
    ifstream openOk((*file).c_str()); 

    //check if a wildcard was used because ifstream can't open wildcards
    Int_t stars=(*file).find("*");
    Int_t quest=(*file).find("?");

    //check if file existed
    if (!openOk && !(quest>=0 || stars>=0)){
      MSG("CDAnalysis",Msg::kInfo)
        <<endl<<endl
        <<"***********************************************************"
        <<endl<<"Can't find file="<<*file<<endl
        <<"Note: you can't use '~/'. It has to be the full path"<<endl
        <<"***********************************************************"
        <<endl<<endl
        <<"Exiting here!"<<endl;
      exit(0);
    }
    
    MSG("CDAnalysis",Msg::kInfo)<<"Adding file="<<*file<<endl;
    nf+=fTrackerTree->Add((*file).c_str());
    fOptTree->Add((*file).c_str());
  }

  if(nf==0){
    MSG("CDAnalysis",Msg::kFatal)
      <<endl<<endl
      <<"*************************************************************"
      <<endl<<"No Tracker*.root files found"<<endl
      <<"Please set CDDATA to the directory containing the"
      <<" Tracker*.root files"<<endl
      <<"Or give the txt file containing the files to be input"<<endl
      <<"Note: If more than one file is found they will be"
      <<" concatenated in a TChain and treated as one"<<endl
      <<"*************************************************************"
      <<endl<<endl<<"Program will exit here"<<endl;
    exit(0);
  }
    
  MSG("CDAnalysis",Msg::kInfo) 
    <<"Tracker Tree information:"<<endl;
  fTrackerTree->Show(0);
  MSG("CDAnalysis",Msg::kInfo) 
    <<"Options Tree information:"<<endl;
  fOptTree->Show(0);
  //fTrackerTree->Print()
  
  MSG("CDAnalysis",Msg::kInfo)
    <<endl<<"Analysing "<<nf<<" file(s). Reading from disk..."<<endl;
}

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

void CDAnalysis::SetRootFileObjects()
{
  MSG("CDAnalysis",Msg::kDebug)
    <<"Running the SetRootFileObjects method..."<<endl;

  //tracker options tree
  fOptTree->SetBranchAddress("SimFlag",&fSimFlag);
  fOptTree->SetBranchAddress("Run",&fRunNumber);
  fOptTree->SetBranchAddress("SubRun",&fSubRun);
  fOptTree->SetBranchAddress("TrackerOptions",fTrkOpt);
  fOptTree->SetBranchAddress("StartTime",&fStartTime);
  fOptTree->SetBranchAddress("EndTime",&fEndTime);
  fOptTree->SetBranchAddress("StartSnarl",&fStartSnarl);
  fOptTree->SetBranchAddress("EndSnarl",&fEndSnarl);
  fOptTree->SetBranchAddress("BeamMomentum",&fBeamMomentum);

  //set the number of events in the tree
  Int_t events=static_cast<Int_t>(fOptTree->GetEntries());
  MSG("CDAnalysis",Msg::kDebug)
    <<"OptTree has "<<events<<" events"<<endl;
  if (events>0){
    //get the snarl/event
    fOptTree->GetEvent(0); 
    MSG("CDAnalysis",Msg::kInfo)
      <<"File has run="<<fRunNumber<<", subrun="<<fSubRun
      <<", SimFlag="<<SimFlag::AsString(fSimFlag)
      <<", beam momentum="<<fBeamMomentum<<endl
      <<"Time="<<fStartTime<<"->"<<fEndTime
      <<", snarl="<<fStartSnarl<<"->"<<fEndSnarl<<endl; 
  }

  //set the branch to point at the track info object and pid
  fTrackerTree->SetBranchAddress("TrackInfo",&fTrackInfo);
  fTrackerTree->SetBranchAddress("Snarl",&fSnarl);
  fTrackerTree->SetBranchAddress("NumDeadChips",&fNumDeadChips);
  fTrackerTree->SetBranchAddress("Sec",&fSec);
  fTrackerTree->SetBranchAddress("Temperature",&fTemperature);

  //check and set pid info branch
  if (fTrackerTree->GetBranchStatus("PIDInfo")){
    fTrackerTree->SetBranchAddress("PIDInfo",&fPIDInfo);
  }
  else{
    MSG("CDAnalysis",Msg::kInfo)
      <<"PIDInfo branch does not exist, not setting tree address"<<endl;
  }

  //create the TClonesArrays
  fTrkHitInfo=new TClonesArray("CDTrackedHitInfo",1000);
  fUnTrkHitInfo=new TClonesArray("CDTrackedHitInfo",1000);
  fXTalkHits=new TClonesArray("CDXTalkHitInfo",1000);

  //set the branches to point to the tclonesarrays
  fTrackerTree->SetBranchAddress("StraightTrackedHitInfo",
                                 &fTrkHitInfo);
  fTrackerTree->SetBranchAddress("UnTrackedHitInfo",
                                 &fUnTrkHitInfo);
  fTrackerTree->SetBranchAddress("XTalkHitInfo",&fXTalkHits);

  //check and set truth hit info branch
  if (fTrackerTree->GetBranchStatus("TruthHitInfo")){
    fTruthHitInfo=new TClonesArray("CDTruthHitInfo",1000);
    fTrackerTree->SetBranchAddress("TruthHitInfo",&fTruthHitInfo);
  }
  else{
    MSG("CDAnalysis",Msg::kInfo)
      <<"TruthHitInfo branch does not exist"
      <<", not setting tree address"<<endl;
  }

  //set the number of events in the tree
  fEvents=static_cast<Int_t>(fTrackerTree->GetEntries());

  //get the first snarl/event
  if (fTrackerTree) fTrackerTree->GetEvent(0); 
  MSG("CDAnalysis",Msg::kInfo)
    <<"First snarl="<<fSnarl<<" (sec="<<fSec<<")"<<endl;

  //this is a really ugly, ugly, ugly place to do this...
  Calibrator& cal=Calibrator::Instance();  
  if (fRunNumber>=70000 && fRunNumber<80000){
    cout<<"Setting DriftCalibrator=PulserDriftCalScheme"<<endl;
    cal.Set("DriftCalibrator=PulserDriftCalScheme");
    cal.GetDriftCalibrator().Set("VATemperatureCalibration=on");
    cal.PrintConfig(cout);
  }
  else{//no tables for 2002 for using PulserDriftCalScheme
    cout<<"Setting DriftCalibrator=PulserSigLinCalScheme"<<endl;
    cal.Set("DriftCalibrator=PulserSigLinCalScheme");
    cal.GetDriftCalibrator().Set("VATemperatureCalibration=on");
    cal.PrintConfig(cout);
  }

  MSG("CDAnalysis",Msg::kDebug)
      <<"Finished the SetRootFileObjects method"<<endl;
}
//......................................................................

Int_t CDAnalysis::RunNumber2RunPeriod(Int_t run) const
{
  //Key to Run Periods (as defined by CBS)
  //1   T11 2002 Near-Far        (30000s)
  //2   T11 2002 Far-Far         (40000s)
  //3   T11 2002 Far-Far   15deg (late 40000s)
  //4   T7  2002 Far-Far         (50000s)
  //5   T7  2003 Near-Far        (70000s)
  //6   T7  2003 Near Only 3m    (80000s)
  //7   T7  2003 Near Only 1m    (90000s)
  //8   T11 2003 Near Only norm  (100000s)
  //9   T11 2003 Near Only 45deg (101000s)
  //10  T11 2003 Near Only 30deg (late 101000s)
  //11  T11 2003 Near Only 15deg (102000s)

  Int_t runPeriod=0;
  if (run>=30000 && run<40000){
    runPeriod=1;
  }
  if (run>=40000 && run<50000){
    runPeriod=2;
  }
  else if (run>=50000 && run<60000){
    runPeriod=4;
  }
  else if (run>=60000 && run<70000){
    runPeriod=-1;//this is not defined
  }
  else if (run>=70000 && run<80000){
    runPeriod=5;
  }
  else if (run>=80000 && run<90000){
    runPeriod=6;
  }
  else if (run>=90000 && run<100000){
    runPeriod=7;
  }
  else if (run>=100000 && run<101000){
    runPeriod=8;
  }
  else if (run>=101000 && run<102000){
    runPeriod=9;//or 10, need run number
  }
  else if (run>=102000 && run<120000){
    runPeriod=11;
  }
  
  return runPeriod;
}

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

void CDAnalysis::SetSigCorPerMEU()
{
  if (fSimFlag==SimFlag::kData){
    if (fRunNumber>70000 && fRunNumber<80000){
      fSigCorPerMEU=591.6;
    }
    else if (fRunNumber==100161){
      fSigCorPerMEU=377.2;
    }
    else{
      MSG("CDAnalysis",Msg::kWarning) 
        <<"no run number for sigcor per meu"<<endl;
    }
  }
  else if(fSimFlag==SimFlag::kMC){
    fSigCorPerMEU=667.104;
  }
  else{
    MSG("CDAnalysis",Msg::kWarning) 
      <<"SetSigCorPerMEU: SimFlag not set"<<endl;
  }

  MSG("CDAnalysis",Msg::kInfo) 
    <<"Using sigcor per MEU="<<fSigCorPerMEU<<endl;
}

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

void CDAnalysis::SetScEnDepPerMEU()
{
  //using same selection of tracks as in data
  fScEnDepPerMEU=0.00195253*Munits::gigaelectronvolt;//new MK dEdx
  //fScEnDepPerMEU=0.00194761*Munits::gigaelectronvolt;//new MK dEdx
  //fScEnDepPerMEU=0.0296718*Munits::gigaelectronvolt;//orig number

  MSG("CDAnalysis",Msg::kInfo) 
    <<"Using ScEnDep per MEU="<<fScEnDepPerMEU<<endl;
}

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

void CDAnalysis::SetRatioScEnToBeamEn()
{
  //using same selection of tracks as in data
  fRatioScEnToBeamEn=0.056064;//with new MK dEdx
  //fRatioScEnToBeamEn=0.0559147;//with new MK dEdx
  //fRatioScEnToBeamEn=0.05638;//orig number
  
  MSG("CDAnalysis",Msg::kInfo) 
    <<"Using Ratio ScEn/BeamEn="<<fRatioScEnToBeamEn<<endl;
}

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

TFile* CDAnalysis::OpenFile(Int_t runNumber,std::string prefix)
{
  //create the tfile pointer to be returned
  TFile* outputFile=0;
  
  //get the environmental variable
  char *anaDir=getenv("CDANA_DIR");
  
  //use a string to hold env instead 
  string sAnaDir="";
  
  if (anaDir!=NULL) {
    sAnaDir=anaDir;
  }
  else {
    MSG("CDAnalysis",Msg::kInfo) 
      <<"Environmental variable $CDANA_DIR not set." 
      <<" Writing file(s) to current directory"<<endl;
    sAnaDir=".";
  }
  
  //convert variables to string
  string sRunNumber=Form("%d",runNumber);
  string sDetector="C";
  string sPrefix="h";//default
  if (prefix!="") sPrefix+=prefix;
  //start all filenames with an h
  string sBase=sAnaDir+"/"+sPrefix+sDetector+sRunNumber;
  string sFileName=sBase+".root";
  
  //test if file already exists
  ifstream Test(sFileName.c_str());
  
  //open the appropriate file
  if(!Test){
    outputFile=new TFile(sFileName.c_str(),"RECREATE");
  }
  else {
    //Need new filename
    Int_t fred=1;
    while(Test) {
      Test.close();
      string sAppendage=Form("%d",fred);
      sFileName=sBase+"_"+sAppendage+".root";
      Test.open(sFileName.c_str());
      fred++;
    }
    outputFile=new TFile(sFileName.c_str(),"NEW");
    outputFile->SetCompressionLevel(9);
  }
  
  string sTmp="No File!";
  if (outputFile) sTmp=outputFile->GetName();

  MSG("CDAnalysis",Msg::kInfo) 
    <<"Output file opened: "<<sTmp<<endl;
  return outputFile;
}

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

ofstream* CDAnalysis::OpenTxtFile(Int_t runNumber,std::string prefix)
{
  //create the tfile pointer to be returned
  ofstream* outputFile=0;
  
  //convert variables to string
  string sRunNumber=Form("%d",runNumber);
  string sDetector="C";
  string sPrefix="";//default
  string sAnaDir=".";
  if (prefix!="") sPrefix+=prefix;
  string sBase=sAnaDir+"/"+sPrefix+sDetector+sRunNumber;
  string sFileName=sBase+".txt";

  outputFile=new ofstream(sFileName.c_str());
  
  if (outputFile){
    MSG("CDAnalysis",Msg::kInfo) 
      <<"Output txt file opened: "<<sFileName<<endl;
  }
  else{
    MSG("CDAnalysis",Msg::kInfo)
      <<"Txt file NOT opened: "<<sFileName<<endl;
  }

  return outputFile;
}

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

TGraph* CDAnalysis::TGraphVect(const std::vector<Double_t>& vX,
                               const std::vector<Double_t>& vY) const
{
  MSG("CDAnalysis",Msg::kDebug) 
    <<" ** Running TGraphVect method... **"<<endl;

  TGraph* g=new TGraph(vX.size());

  for (UInt_t i=0;i<vX.size();i++){
    g->SetPoint(i,vX[i],vY[i]);
  }

  MSG("CDAnalysis",Msg::kDebug) 
    <<" ** Finished TGraphVect method **"<<endl;
  return g;
}
//......................................................................

void CDAnalysis::TGraphMinMax(vector<TGraph*>& v) const
{
  Float_t max=-1e200;
  Float_t min=+1e200;

  for (vector<TGraph*>::iterator vIt=v.begin();vIt!=v.end();++vIt){
    for (Int_t i=0;i<(*vIt)->GetN();i++){
      Double_t x=0;
      Double_t y=0;
      (*vIt)->GetPoint(i,x,y);

      //find the min and max
      if (y>max) max=y;
      if (y<min) min=y;
    }
  }

  MSG("CDAnalysis",Msg::kDebug)
    <<"Found min="<<min<<", max="<<max<<endl;

  //now set the min and max for 10% above and below size of range
  for (vector<TGraph*>::iterator vIt=v.begin();vIt!=v.end();++vIt){
    (*vIt)->SetMinimum(min-0.1*(max-min));
    (*vIt)->SetMaximum(max+0.1*(max-min));
  }
}

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

TGraph* CDAnalysis::TGraphMap(const std::map<Int_t,Float_t>& m) const
{
  MSG("CDAnalysis",Msg::kDebug) 
    <<" ** Running TGraphMap method... **"<<endl;

  TGraph* g=new TGraph(m.size());

  map<Int_t,Float_t>::const_iterator mIter=m.begin();

  Int_t i=0;

  while (mIter!=m.end()){
    g->SetPoint(i,mIter->first,mIter->second);

    i++;
    mIter++;
  }

  MSG("CDAnalysis",Msg::kDebug) 
    <<" ** Finished TGraphMap method **"<<endl;
  return g;
}

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

Int_t CDAnalysis::SumMapValues(const std::map<Int_t,Int_t>& mapIn) const
{
  Int_t sum=0;
  
  //loop and sum the values stored in the map
  for (map<Int_t,Int_t>::const_iterator it=mapIn.begin();
       it!=mapIn.end();++it){
    sum+=it->second;
  }
  
  return sum;
}

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

void CDAnalysis::NormaliseVector(std::vector<Double_t>& v,
                                 Int_t mode) const
{
  MSG("CDAnalysis",Msg::kDebug) 
    <<" ** Running NormaliseVector method... **"<<endl;
  
  if (v.size()>0){

    if (mode==1){//simple average
      //find average
      Double_t average=0;
      for (UInt_t i=0;i<v.size();i++) average+=v[i];
      average/=v.size();
      
      MSG("CDAnalysis",Msg::kVerbose)<<"Average="<<average<<endl; 
      
      //normalise
      for (UInt_t i=0;i<v.size();i++) v[i]/=average;
    }
    else if (mode==2){//fit the function to remove effect of outliers
      
      Double_t vMax=-9e50;
      Double_t vMin=9e50;

      //work out mins and maxs
      for (UInt_t i=0;i<v.size();i++){
        if (v[i]>vMax) vMax=v[i];
        if (v[i]<vMin) vMin=v[i];
      }
      
      //limit the number of bins
      Int_t numBins=300;
      if (4*v.size()<10000) numBins=4*v.size();

      MSG("LIAnalysis",Msg::kInfo)
        <<"vMax="<<vMax <<", vMin="<<vMin<<", numBins="<<numBins<<endl;

      //create the histo
      TH1F* h=new TH1F("h","h",numBins,vMin-abs(0.1*vMin),
                       vMax+abs(0.1*vMin));
      h->SetBit(TH1::kCanRebin);

      Double_t average=0;
      for (UInt_t i=0;i<v.size();i++){
        //find average
        average+=v[i];

        //fill the histo
        h->Fill(v[i]);
      }
      average/=v.size();

      //fit but use 0 (zero!) otherwise it stamps on the current canvas!
      h->Fit("gaus","q0");
      TF1* func=h->GetFunction("gaus");
      Double_t p1=func->GetParameter(1);
      Double_t p2=func->GetParameter(2);

      //don't use the fit if it is screwed up
      Double_t avToNormWith=p1;
      if (p1<vMin || p1>vMax) avToNormWith=average;

      MSG("LIAnalysis",Msg::kInfo)
        <<"Fitted mean="<<p1
        <<" (average="<<average
        <<") fitted rms="<<p2
        <<", using avToNormWith="<<avToNormWith<<endl;

      //normalise
      for (UInt_t i=0;i<v.size();i++) v[i]/=avToNormWith;

      delete h;
    }
  }

  MSG("CDAnalysis",Msg::kDebug) 
    <<" ** Finished NormaliseVector method **"<<endl;
}

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

void CDAnalysis::SetLoopVariables(Int_t event)
{
  Float_t fract=ceil(fEvents/20.);
  
  if (ceil(((Float_t)event)/fract)==((Float_t)event)/fract){
    MSG("CDAnalysis",Msg::kInfo) 
      <<"Fraction of loop complete: "<<event 
      <<"/"<<fEvents<<"  ("
      <<(Int_t)(100.*event/fEvents)<<"%)"<<endl;
  }

  //get the snarl/event
  fTrackerTree->GetEvent(event); 

  //initialise snarl level variables
  this->InitialiseLoopVariables();

  this->InitialisePidVariables();

  //get the pid info
  if (fPIDInfo){
    fTriggerPMT=fPIDInfo->GetTriggerPMT();
    fFafErr=fPIDInfo->GetFafErr() ;
    fSparseErr=fPIDInfo->GetSparseErr();
    fTrigSource=fPIDInfo->GetTrigSource();
    fKovADC1=fPIDInfo->GetKovADC1() ;
    fKovTimeStamp1=fPIDInfo->GetKovTimeStamp1();
    fKovADC2=fPIDInfo->GetKovADC2();
    fKovTimeStamp2=fPIDInfo->GetKovTimeStamp2();
    fKovADC3=fPIDInfo->GetKovADC3();
    fKovTimeStamp3=fPIDInfo->GetKovTimeStamp3(); 
    fSnarlTimeFrame=fPIDInfo->GetSnarlTimeFrame();
    fSnarlMinTimeStamp=fPIDInfo->GetSnarlMinTimeStamp();
    fSnarlMaxTimeStamp=fPIDInfo->GetSnarlMaxTimeStamp();
    fTofTDC0=fPIDInfo->GetTofTDC0();//actual time of flight
    fTofTDC1=fPIDInfo->GetTofTDC1();
    fTofTDC2=fPIDInfo->GetTofTDC2();
    fTofADC0=fPIDInfo->GetTofADC0();
    fTofADC1=fPIDInfo->GetTofADC1();
    fTofADC2=fPIDInfo->GetTofADC2();
    /*
      kTofADCChannel0=ph.GetSpecialChannelContains("TADC0");
      kTofADCChannel1=ph.GetSpecialChannelContains("TTAG");
      if(kTofADCChannel1==0){
         kTofADCChannel1=ph.GetSpecialChannelContains("TADC1");
      }
      kTofADCChannel2=ph.GetSpecialChannelContains("TimingFid");
      if(kTofADCChannel2==0){
         kTofADCChannel2=ph.GetSpecialChannelContains("TADC2");
      }
      kTofTimeStampChannel=ph.GetSpecialChannelContains("TTIME");
    */
    fTofADCTimeStamp0=fPIDInfo->GetTofADCTimeStamp0();//TADC0
    fTofADCTimeStamp1=fPIDInfo->GetTofADCTimeStamp1();//TTAG || TADC1
    fTofADCTimeStamp2=fPIDInfo->GetTofADCTimeStamp2();//TimingFid |TADC2
    //TofTimeStamp=fPIDInfo->GetTofTimeStamp();//don't use
    fTickSinceLast=fPIDInfo->GetTickSinceLast();
    
    //these are the PID variables
    fNoOverlap=fPIDInfo->GetNoOverlap();
    fInCERTime=fPIDInfo->GetInCERTime();
    fPIDType=fPIDInfo->GetPIDType();
    fNoOverlapBits=fPIDInfo->GetNoOverlapBits();
    fInCERTimeBits=fPIDInfo->GetInCERTimeBits();
    fOLChi2=fPIDInfo->GetOLChi2();
  }
}  

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

void CDAnalysis::PrintBasicInfo(std::string s) const
{
  MSG("CDAnalysis",Msg::kInfo)
    <<s<<"(Pl,St,End)=("<<fPlane<<","<<fStrip<<","<<fStripend<<")"
    <<", adc="<<fChargeAdc
    <<", sigL="<<fChargeSigLin
    <<", sigC="<<fChargeSigCor
    <<", mip="<<fChargeMip
    <<", pe="<<fChargePe
    <<endl;
}

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

void CDAnalysis::PrintEvent(Int_t event)
{
  if (event>=fEvents) {
    MSG("CDAnalysis",Msg::kFatal)
      <<"No such event in tree="<<event<<", fEvents="<<fEvents<<endl;
    return;
  }

  MSG("CDAnalysis",Msg::kInfo)
    <<"Printing for event="<<event<<", fEvents="<<fEvents<<endl;

  //get the snarl/event
  fTrackerTree->GetEvent(event); 

  //get tclones arrays for this snarl
  TClonesArray &cTrk=*fTrkHitInfo;
  Int_t numTrkHits=fTrkHitInfo->GetEntries();
  TClonesArray &cUnTrk = *fUnTrkHitInfo;
  Int_t numUnTrkHits = fUnTrkHitInfo->GetEntries();
  TClonesArray &cXTalk = *fXTalkHits;
  Int_t numXTalkHits=fXTalkHits->GetEntries();
  
  //loop over the untracked hits in the snarl
  MSG("CDAnalysis",Msg::kInfo)<<"Untracked hits..."<<endl;
  for (Int_t hit=0;hit<numUnTrkHits;hit++){
    //cast the tclonesarray up to a trackedhitinfo object
    CDTrackedHitInfo *hitInfo=
      dynamic_cast<CDTrackedHitInfo*>(cUnTrk[hit]);
    
    this->ReadInHitInfo(hitInfo);
    
    this->PrintBasicInfo();
  }
  
  //loop over xtalk hits in the snarl
  MSG("CDAnalysis",Msg::kInfo)<<"XTalk hits..."<<endl;
  for (Int_t hit=0;hit<numXTalkHits;hit++){
    CDXTalkHitInfo *hitInfo=
      dynamic_cast<CDXTalkHitInfo*>(cXTalk[hit]);
    
    this->ReadInHitInfo(hitInfo);

    this->PrintBasicInfo();
  }

  //loop over the tracked hits
  MSG("CDAnalysis",Msg::kInfo)<<"Tracked hits..."<<endl;
  for (Int_t hit=0;hit<numTrkHits;hit++){
    CDTrackedHitInfo* hitInfo=
      dynamic_cast<CDTrackedHitInfo*>(cTrk[hit]);

    this->ReadInHitInfo(hitInfo);
    
    this->PrintBasicInfo();
  }
}

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

void CDAnalysis::CalcNumPlanesHit(Int_t plane,Float_t pe)
{
  //calc the num planes hit
  fNumPlanesHitAll[plane]++;

  //calc number of planes hit in first 25
  if (plane<25) fNumPlanesHit25[plane]++;

  //calc if pe cut satisfied
  if (pe>1.5) {
    fNumPlanesHitPeCut[plane]++;

    //only calc planes hit in first 10 planes
    if (plane<10) fNumPlanesHitPeCut10[plane]++;
  }
}

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

void CDAnalysis::CalcNumPlanesHitTrk(Int_t plane)
{
  //calc the num planes hit
  fNumPlanesHitTrk[plane]++;
}

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

void CDAnalysis::CalcAvStripForPSMuonCut(Int_t plane,Int_t strip,
                                         Int_t end)
{
  if (plane<4){
    //calc weighted average strip
    fAvStrip=fAvStrip*(fStCount/(fStCount+1))+strip*(1/(fStCount+1));
    fStCount++;
    MAXMSG("CDAnalysis",Msg::kVerbose,1000)
      <<"  strip="<<strip<<", plane="<<plane
      <<", fEnd="<<fStripend<<", fAvStrip="<<fAvStrip
      <<", fStCount="<<fStCount<<endl;

    //now do each end separately
    if (end==1 && plane>0) {
      fAvStrip1=fAvStrip1*(fStCount1/(fStCount1+1))+
        strip*(1/(fStCount1+1));
      fStCount1++;
    }
    else if (end==2 && plane>0) {
      fAvStrip2=fAvStrip2*(fStCount2/(fStCount2+1))+
        strip*(1/(fStCount2+1));
      fStCount2++;
    }
  }
}

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

void CDAnalysis::CalcLowUpScint_TofDiff(Double_t time)
{
  Double_t tickInNs=25./16;//ns, va clock tick
  Float_t timeDiff=(time*1e9)-(fTofADCTimeStamp1*tickInNs);//ns

  if (fTofADCTimeStamp1<1 || fTofADCTimeStamp1>1e9/tickInNs){
    MAXMSG("CDAnalysis",Msg::kWarning,10)
      <<"Bad Tof (TTAG) time="<<fTofADCTimeStamp1
      <<" ticks, ("<<fTofADCTimeStamp1*tickInNs<<" ns)"<<endl;
  }

  if (timeDiff>1000 || timeDiff<-1000){
    MAXMSG("CDAnalysis",Msg::kWarning,10)
      <<"Bad Scint-Tof (TTAG) time difference="<<timeDiff<<" ns"
      <<", time="<<time<<", tof time="<<fTofADCTimeStamp1
      <<" ticks, ("<<fTofADCTimeStamp1*tickInNs<<" ns)"<<endl;
  }

  //calc highest and lowest scint tof time diff.
  if (timeDiff>fUpScint_Tof) fUpScint_Tof=timeDiff;
  if (timeDiff<fLowScint_Tof) fLowScint_Tof=timeDiff;
}

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

void CDAnalysis::CalcFirstLastPlane(Int_t plane)
{
  if (plane<fFirstPlane) fFirstPlane=plane;
  if (plane>fLastPlane) fLastPlane=plane;
  if (plane>fLastPlaneOdd && plane%2!=0) fLastPlaneOdd=plane;
  if (plane>fLastPlaneEven && plane%2==0) fLastPlaneEven=plane;
}

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

Int_t CDAnalysis::CalcLastPlaneTruth() const
{  
  if (!fTruthHitInfo) return 0;

  TClonesArray &cTruth = *fTruthHitInfo;
  Int_t numTruthHits=fTruthHitInfo->GetEntries();

  Int_t lastPlane=0;

  //loop over the truth hits
  for (Int_t hit=0;hit<numTruthHits;hit++){
    //cast the tclonesarray up to a truthhitinfo object
    CDTruthHitInfo *hitInfo=
      dynamic_cast<CDTruthHitInfo*>(cTruth[hit]);
    
    Int_t plane=hitInfo->GetPlane();
    Int_t pmtTruth1=hitInfo->GetPmtTruth1();
    Int_t pmtTruth2=hitInfo->GetPmtTruth2();

    //check either side
    Int_t sumTruth=pmtTruth1|pmtTruth2;

    //check if genuine
    if ((sumTruth & DigiSignal::kGenuine)==DigiSignal::kGenuine){
      //see if bigger
      if (plane>lastPlane) lastPlane=plane;
    }
    MSG("CDAnalysis",Msg::kVerbose)
      <<"sumTruth="<<sumTruth<<endl;
  }

  MSG("CDAnalysis",Msg::kVerbose)
    <<"CalcLastPlaneTruth: LastPlane="<<lastPlane
    <<", fLastPlane="<<fLastPlane<<endl; 

  return lastPlane;
}

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

void CDAnalysis::CalcLastPlaneOnTrk()
{

  TClonesArray &cTrk=*fTrkHitInfo;
  Int_t numTrkHits=fTrkHitInfo->GetEntries();

  //loop over the tracked hits
  for (Int_t hit=0;hit<numTrkHits;hit++){
    CDTrackedHitInfo *trackedHitInfo=
      dynamic_cast<CDTrackedHitInfo*>(cTrk[hit]);
    
    this->ReadInHitInfo(trackedHitInfo);
    
    //cut out the bad channels
    if (this->CutOnBadPedestals(fPlane,fStrip,fStripend)) continue;
    
    //check the channel's sanity
    this->StandardSanityChecks();
    
    //calculate the first and last planes
    this->CalcFirstLastPlane(fPlane);
    
  }
  //end of loop over the tracked hits
  
  MSG("CDAnalysis",Msg::kVerbose)
    <<"CalcLastPlaneOnTrk: fLastPlane="<<fLastPlane<<endl; 
}

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

Int_t CDAnalysis::CalcLastPlaneOnTrkNoXTalk() const
{
  //get tclones arrays for this snarl
  TClonesArray &cTrk=*fTrkHitInfo;
  Int_t numTrkHits=fTrkHitInfo->GetEntries();
  TClonesArray &cUnTrk = *fUnTrkHitInfo;
  Int_t numUnTrkHits = fUnTrkHitInfo->GetEntries();
  TClonesArray &cXTalk = *fXTalkHits;
  Int_t numXTalkHits=fXTalkHits->GetEntries();
  
  map<Int_t,Float_t> planeGreatestPe;
  Int_t lastPlane=0;  

  //the untracked and xtalk hits are just used to find the greatest
  //pe in each plane; they are not included in the last plane 
  //directly

  //loop over the untracked hits in the snarl
  for (Int_t hit=0;hit<numUnTrkHits;hit++){
    //cast the tclonesarray up to a trackedhitinfo object
    CDTrackedHitInfo *hitInfo=
      dynamic_cast<CDTrackedHitInfo*>(cUnTrk[hit]);
    
    //get the variables
    Int_t plane=hitInfo->GetPlane();
    Int_t strip=hitInfo->GetStrip();
    Int_t stripend=hitInfo->GetEnd();
    Float_t chargePe=hitInfo->GetCharge(CDTrackedHitInfo::kPe);

    //cut out the bad channels
    if (this->CutOnBadPedestals(plane,strip,stripend)) continue;
    
    if (chargePe>planeGreatestPe[fPlane]) planeGreatestPe[plane]=
                                            chargePe;
  }
  
  //loop over xtalk hits in the snarl
  for (Int_t hit=0;hit<numXTalkHits;hit++){
    CDXTalkHitInfo *hitInfo=
      dynamic_cast<CDXTalkHitInfo*>(cXTalk[hit]);
    
    //get the variables
    Int_t plane=hitInfo->GetPlane();
    Int_t strip=hitInfo->GetStrip();
    Int_t stripend=hitInfo->GetEnd();
    Float_t chargePe=hitInfo->GetCharge(CDTrackedHitInfo::kPe);

    //cut out the bad channels
    if (this->CutOnBadPedestals(plane,strip,stripend)) continue;
    
    if (chargePe>planeGreatestPe[fPlane]) planeGreatestPe[plane]=
                                            chargePe;
  }

  //loop over the tracked hits
  for (Int_t hit=0;hit<numTrkHits;hit++){
    CDTrackedHitInfo* hitInfo=
      dynamic_cast<CDTrackedHitInfo*>(cTrk[hit]);
    
    //get the variables
    Int_t plane=hitInfo->GetPlane();
    Int_t strip=hitInfo->GetStrip();
    Int_t stripend=hitInfo->GetEnd();
    Float_t chargePe=hitInfo->GetCharge(CDTrackedHitInfo::kPe);

    if (plane>lastPlane) lastPlane=plane;

    //cut out the bad channels
    if (this->CutOnBadPedestals(plane,strip,stripend)) continue;
    
    if (chargePe>planeGreatestPe[fPlane]) planeGreatestPe[plane]=
                                            chargePe;
  }
  //end of loop over the tracked hits

  Int_t newLastPlane=lastPlane;
  Double_t xtalkLevel=1.7;

  if (fLastPlane>=3){//make sure track is at least 4 planes
    if (planeGreatestPe[lastPlane]>xtalkLevel){
      //do nothing
    }
    else if (planeGreatestPe[lastPlane-1]>xtalkLevel){
      newLastPlane=lastPlane-1;
    }
    else if (planeGreatestPe[lastPlane-2]>xtalkLevel){
      newLastPlane=lastPlane-2;
    }
    else{
      if (planeGreatestPe[lastPlane-3]<xtalkLevel){
        MSG("CDAnalysis",Msg::kDebug)
          <<"4th plane with low pe hit"
          <<", pe0="<<planeGreatestPe[lastPlane-0]
          <<", pe1="<<planeGreatestPe[lastPlane-1]
          <<", pe2="<<planeGreatestPe[lastPlane-2]
          <<", pe3="<<planeGreatestPe[lastPlane-3]
          <<endl;
      }
      newLastPlane=lastPlane-3;
    }
  }
  
  MSG("CDAnalysis",Msg::kVerbose)
    <<"CalcLastPlaneOnTrkNoXTalk: lastPlane="<<lastPlane
    <<", newLastPlane="<<newLastPlane<<endl; 
  
  return newLastPlane;
}

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

void CDAnalysis::ReadInHitInfo(CDTrackedHitInfo* hitInfo)
{
  //initialise variables to be on the safe side
  this->InitialiseHitInfoVariables();

  fTime=hitInfo->GetTime();

  //get the variables
  fPlane=hitInfo->GetPlane();
  fStrip=hitInfo->GetStrip();
  fStripend=hitInfo->GetEnd();
  fTransPos=hitInfo->GetTransPos();
  
  //decide which stripend is which; odd, even, 1 or 2.
  fO1=(fStripend==1 && fPlane%2==1);
  fO2=(fStripend==2 && fPlane%2==1);
  fE1=(fStripend==1 && fPlane%2==0);
  fE2=(fStripend==2 && fPlane%2==0);

  //look at the calorimetry info
  fChargeAdc=hitInfo->GetCharge(CDTrackedHitInfo::kAdc);
  fChargeMip=hitInfo->GetCharge(CDTrackedHitInfo::kMip);
  fChargeSigCor=hitInfo->GetCharge(CDTrackedHitInfo::kSigCorr);
  fChargeSigLin=hitInfo->GetCharge(CDTrackedHitInfo::kSigLin);
  fChargePe=hitInfo->GetCharge(CDTrackedHitInfo::kPe);
  if (fChargePe!=0) fGain=fChargeSigLin/fChargePe;
  else fGain=-1;

  static Bool_t firstTime=true;
  if (fChargeMip>0 && fChargeSigCor>0 && firstTime){
    firstTime=false;
    fSigCorsPerMip=fChargeSigCor/fChargeMip;
    MAXMSG("CDAnalysis",Msg::kVerbose,100)
      <<"fSigCorsPerMip="<<fSigCorsPerMip<<endl;
  }

  //get the tracker evaluation of the event
  fResultEven=fTrackInfo->GetResult(0);
  fResultOdd=fTrackInfo->GetResult(1);
}

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

void CDAnalysis::ReadInHitInfo(CDXTalkHitInfo* hitInfo)
{
  //initialise variables to be on the safe side
  this->InitialiseHitInfoVariables();

  fTime=hitInfo->GetTime();

  //get the variables
  fPlane=hitInfo->GetPlane();
  fStrip=hitInfo->GetStrip();
  fStripend=hitInfo->GetEnd();
  //fTransPos=hitInfo->GetTransPos();
  fTransPos=-1;//npt possible to calculate for xtalk

  //decide which stripend is which; odd, even, 1 or 2.
  fO1=(fStripend==1 && fPlane%2==1);
  fO2=(fStripend==2 && fPlane%2==1);
  fE1=(fStripend==1 && fPlane%2==0);
  fE2=(fStripend==2 && fPlane%2==0);

  //look at the calorimetry info
  fChargeAdc=hitInfo->GetCharge(CDTrackedHitInfo::kAdc);
  fChargeMip=hitInfo->GetCharge(CDTrackedHitInfo::kMip);
  fChargeSigCor=hitInfo->GetCharge(CDTrackedHitInfo::kSigCorr);
  fChargeSigLin=hitInfo->GetCharge(CDTrackedHitInfo::kSigLin);
  fChargePe=hitInfo->GetCharge(CDTrackedHitInfo::kPe);
  if (fChargePe!=0) fGain=fChargeSigLin/fChargePe;
  else fGain=-1;
}

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

Float_t CDAnalysis::CalcPathLengthCorection(const TVector3& end1,
                                            const TVector3& end2) const
{
  //calc pathlength between two points
  Float_t x=end1.X()-end2.X();
  Float_t y=end1.Y()-end2.Y();
  Float_t z=end1.Z()-end2.Z();
  Float_t pathLength=sqrt(pow(x,2)+pow(y,2)+pow(z,2));
  //calc path length correction
  Float_t pathLengthCor=-1;
  if (z!=0) pathLengthCor=pathLength/fabs(z);
  else cout<<"pathlength cor problem: z="<<z
           <<", x="<<x<<", y="<<y<<endl;
  return pathLengthCor;
}

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

void CDAnalysis::CalcXorYandZ(Int_t strip,Int_t plane,
                              Float_t& x,Float_t& y,Float_t& z) const
{
  Float_t planePitch=5.94*Munits::cm;
  Float_t stripWidth=4.1*Munits::cm;
  Float_t stripsFromCentre=strip-11.5;//middle 2 strips 11&12 are 0.5 st
  Float_t distFromCentre=stripsFromCentre*stripWidth;
  
  //calc z
  z=plane*planePitch;

  if (plane%2==0){//even, v-view (measure v and y), horizontal strips
    x=-1;
    y=distFromCentre;
  }
  else if (plane%2==1){//odd, u-view (measure u and x), vertical strips
    x=distFromCentre;;
    y=-1;
  }
}

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

void CDAnalysis::TrueHighLowEn(Float_t& highEn,Float_t& lowEn) const
{
  if (!fTruthHitInfo) return;

  //initialise
  lowEn=999999;
  highEn=-1;

  TClonesArray &cTruth = *fTruthHitInfo;
  Int_t numTruthHits=fTruthHitInfo->GetEntries();
  
  //loop over the truth hits
  for (Int_t hit=0;hit<numTruthHits;hit++){
    //cast the tclonesarray up to a truthhitinfo object
    CDTruthHitInfo *hitInfo=
      dynamic_cast<CDTruthHitInfo*>(cTruth[hit]);
    
    Double_t energy=hitInfo->GetMainPartEn();
    Int_t particle=hitInfo->GetMainParticle();
    Int_t muon=13;
    
    if (abs(particle)==muon) {
      if (energy>highEn) highEn=energy;
      if (energy<lowEn) lowEn=energy; 
    }
  }

  //reset to -1 if not found to be consistent with highEn
  if (lowEn==999999) lowEn=-1;
}

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

void CDAnalysis::TruePLCor(std::map<Int_t,Float_t>& mPLCor,
                           Int_t event) const
{
  Int_t lastPlane=this->GetEventLength();
  if (lastPlane<0) return;
  Float_t lastPLCor=-1;
  Float_t lastPL=-1;

  for (Int_t pl=0;pl<=lastPlane;pl++){
    Float_t pathLen=this->TrueMainPL(pl,event);
    //get a zero if no dsh in a plane
    //also know that PL must be at least 0.95 cm
    if (pathLen<0.95*Munits::cm) {
      if (lastPL>=0.95*Munits::cm) pathLen=lastPL;
      else pathLen=0.95*Munits::cm;
    }
    //record the last PL
    lastPL=pathLen;

    Float_t pathLengthCor=pathLen/(0.95*Munits::cm);
    if (pathLengthCor<0.999) {
      MAXMSG("CDAnalysis",Msg::kWarning,500)
        <<"bad PLCor="<<pathLengthCor<<endl;
      if (lastPLCor>=1) pathLengthCor=lastPLCor;
      else pathLengthCor=1;
    }
    else if (pathLengthCor>0.999 && pathLengthCor<1){
      MAXMSG("CDAnalysis",Msg::kVerbose,500)
        <<"fractionally wrong PLCor="<<pathLengthCor<<endl;
      pathLengthCor=1;
    }

    //insert the value into the map
    mPLCor[pl]=pathLengthCor;

    //record the last PLCor
    lastPLCor=pathLengthCor;
  }
}

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

Bool_t CDAnalysis::CalcPLCor(const std::map<Int_t,TVector3>& mPlPos,
                             std::map<Int_t,Float_t>& mPLCor) const
{
  map<Int_t,TVector3>::const_iterator plPosEndIt=mPlPos.end();

  //Int_t lastPlaneNoXTalk=this->CalcLastPlaneOnTrkNoXTalk();
  Int_t lastPlane=this->GetEventLength();
  Float_t initValue=-999;
 
  for (Int_t pl=0;pl<=lastPlane;pl++){
    Int_t plVtxSide=pl-2;
    Int_t plStopSide=pl+2;
    
    //protect against going outside allowed ranges
    if (plVtxSide<0) plVtxSide=0;
    if (plStopSide>lastPlane) plStopSide=lastPlane;

    MAXMSG("CDAnalysis",Msg::kDebug,500)
      <<"pl="<<pl<<", plStopSide="<<plStopSide
      <<", plVtxSide="<<plVtxSide<<endl;

    //fit track position
    TGraph gx_z(3);//min of 3 points
    TGraph gy_z(3);
    Int_t tmpPl=plStopSide;
    Int_t tmpPlVtxSide=plVtxSide+1;
    Bool_t posDir=true;//always forward going in caldet
    if (posDir) tmpPlVtxSide=plVtxSide-1;
    Int_t i=0;
    //loop over planes to fit
    while (tmpPl!=tmpPlVtxSide){
      map<Int_t,TVector3>::const_iterator plPosIt=mPlPos.find(tmpPl);
      if (plPosIt!=plPosEndIt){
        Float_t x=(plPosIt->second).x();
        Float_t y=(plPosIt->second).y();
        Float_t z=(plPosIt->second).z();
        gx_z.SetPoint(i,z,x);
        gy_z.SetPoint(i,z,y);
        MAXMSG("CDAnalysis",Msg::kVerbose,500)
          <<"Adding point: x="<<x<<", y="<<y<<", z="<<z<<endl;
        //count the number of points added
        i++;
      }
      //increment the plane
      if (posDir) tmpPl--;
      else tmpPl++;
    }
  
    Float_t newPLCor=initValue;
    if (i>=3){
      gx_z.Fit("pol1","q");
      gy_z.Fit("pol1","q");
      Float_t mx_z=gx_z.GetFunction("pol1")->GetParameter(1);
      Float_t my_z=gy_z.GetFunction("pol1")->GetParameter(1);
      //Float_t mx_zErr=gx_z.GetFunction("pol1")->GetParError(1);
      //Float_t my_zErr=gy_z.GetFunction("pol1")->GetParError(1);
      newPLCor=sqrt(pow(mx_z,2)+pow(my_z,2)+1);
      MAXMSG("CDAnalysis",Msg::kDebug,500)
        <<"pl="<<pl<<", found plcor="<<newPLCor<<endl;
    }
    else {
      MSG("CDAnalysis",kWarning)
        <<"Less than 3 points, can't do local grad. fit"<<endl;
      return false;
    }

    if (newPLCor!=initValue) mPLCor[pl]=newPLCor;
    else MSG("CDAnalysis",kInfo)<<"Bad PLCor"<<endl;

  } 

  return true;
}

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

Bool_t CDAnalysis::CalcXYZ(std::map<Int_t,TVector3>& mPlPos) const
{
  //get tclones arrays for this snarl
  TClonesArray &cTrk=*fTrkHitInfo;
  Int_t numTrkHits=fTrkHitInfo->GetEntries();
  
  if (numTrkHits<2){
    MAXMSG("CDAnalysis",Msg::kWarning,100)
      <<"CalcXYZ::numTrkHits is bad="<<numTrkHits<<endl;
    return false;
  }

  //Int_t lastPlaneNoXTalk=this->CalcLastPlaneOnTrkNoXTalk();
  Int_t lastPlane=this->GetEventLength();
  if (lastPlane<2 || lastPlane>59){
    MAXMSG("CDAnalysis",Msg::kWarning,100)
      <<"CalcXYZ::last plane is bad="<<lastPlane<<endl;
    return false;
  }

  Bool_t foundFirstFew=false;
  Bool_t unFilledGap=false;
  Float_t initValue=-999;
  vector<TVector3> plPosition(60);
  for (Int_t pl=0;pl<60;pl++){
    plPosition[pl]=TVector3(initValue,initValue,initValue);
  }

  //loop over the tracked hits
  for (Int_t hit=0;hit<numTrkHits;hit++){
    CDTrackedHitInfo* hitInfo=
      dynamic_cast<CDTrackedHitInfo*>(cTrk[hit]);
    
    //get the variables
    Int_t plane=hitInfo->GetPlane();
    Int_t strip=hitInfo->GetStrip();
    Int_t stripend=hitInfo->GetEnd();
    //Float_t chargePe=hitInfo->GetCharge(CDTrackedHitInfo::kPe);

    //cut out the bad channels
    if (this->CutOnBadPedestals(plane,strip,stripend)) continue;

    Float_t x=initValue;
    Float_t y=initValue;
    Float_t z=initValue;
    this->CalcXorYandZ(strip,plane,x,y,z);

    plPosition[plane]=TVector3(x,y,z);

  }
  //end of loop over the tracked hits

  //calc the x and y of the first two planes
  //loop and look at more planes if can't find any hits in first few
  for (Int_t np=2;np<=5;np++){
    //look at first few even planes
    Float_t yAvFirstFew=0;
    Float_t yAvFirstFewCounter=0;
    for (Int_t pl=0;pl<np*2;pl+=2){//even measures v or y
      if (plPosition[pl].y()!=initValue){
        yAvFirstFew+=plPosition[pl].y();
        yAvFirstFewCounter++;
      }
    }
    if (yAvFirstFewCounter>0) yAvFirstFew/=yAvFirstFewCounter;
    else yAvFirstFew=initValue;
    
    //look at first few odd planes
    Float_t xAvFirstFew=0;
    Float_t xAvFirstFewCounter=0;
    for (Int_t pl=1;pl<np*2;pl+=2){//odd measures u or x
      if (plPosition[pl].x()!=initValue){
        xAvFirstFew+=plPosition[pl].x();
        xAvFirstFewCounter++;
      }
    }
    if (xAvFirstFewCounter>0) xAvFirstFew/=xAvFirstFewCounter;
    else xAvFirstFew=initValue;
    
    MAXMSG("CDAnalysis",Msg::kDebug,1000)
      <<"num planes="<<np
      <<", avX="<<xAvFirstFew<<" ("<<xAvFirstFewCounter<<")"
      <<", avY="<<yAvFirstFew<<" ("<<yAvFirstFewCounter<<")"<<endl;

    //now set the TVectors if an average was obtained
    if (xAvFirstFew==initValue || yAvFirstFew==initValue){
      MAXMSG("CDAnalysis",Msg::kVerbose,100)
        <<endl<<"Not found tracked hits with num planes="<<np
        <<", looping again..."<<endl<<endl;
    }
    else{
      //set the first 2 planes (what is not measured by that plane)
      plPosition[0].SetX(xAvFirstFew);
      plPosition[1].SetY(yAvFirstFew);
      
      //sometimes have to set what should have been measured if 
      //it is missing
      if (plPosition[0].Y()==initValue || plPosition[1].X()==initValue){
        MAXMSG("CDAnalysis",Msg::kDebug,100)
          <<endl<<"Having to manually set x or y in first few planes"
          <<endl<<endl;
        //plPosition[0].Print();
        //plPosition[1].Print();
      }
      if (plPosition[0].Y()==initValue) plPosition[0].SetY(yAvFirstFew);
      if (plPosition[1].X()==initValue) plPosition[1].SetX(xAvFirstFew);
      foundFirstFew=true;
      break;
    }
  }

  if (!foundFirstFew){
    MAXMSG("CDAnalysis",Msg::kDebug,100)
      <<endl<<"Not found hits in first few planes ignoring this event"
      <<endl<<endl;
    return false;
  }

  //probably best to do the +/- 2 planes first to fill the gaps
  //then do the interpolation into the planes that don't measure
  //what is needed.

  //loop and fill the gaps
  for (Int_t pl=2;pl<=lastPlane;pl++){

    Bool_t foundGap=false;
    
    //check for gaps
    if (pl%2==0){//even planes, measures v and y
      if (plPosition[pl].Y()==initValue) foundGap=true;
    }
    else if (pl%2==1){//odd planes, measures u and x
      if (plPosition[pl].X()==initValue) foundGap=true;
    }

    if (foundGap) {
      Int_t plPrev=pl-2;
      Int_t plNext=pl+2;
      Int_t plPrev4=pl-4;
      Int_t plNext4=pl+4;
      Float_t posPrev=initValue;
      Float_t posNext=initValue;
      Float_t pos=initValue;
      Float_t posPrev4=initValue;
      Float_t posNext4=initValue;
      Float_t pos4=initValue;

      //protect against going outside allowed ranges
      if (plPrev<0) {
        if (abs(plPrev%2)==0) plPrev=0;//even
        else if (abs(plPrev%2)==1) plPrev=1;//odd
      }
      if (plNext>59){
        if (plNext%2==0) plNext=58;//even
        else if (plNext%2==1) plNext=59;//odd
      }
      if (plPrev4<0) {
        if (abs(plPrev4%2)==0) plPrev4=0;//even
        else if (abs(plPrev4%2)==1) plPrev4=1;//odd
      }
      if (plNext4>59){
        if (plNext4%2==0) plNext4=58;//even
        else if (plNext4%2==1) plNext4=59;//odd
      }
      
      if (pl%2==0){//even planes, measures v and y
        //get position
        posPrev=plPosition[plPrev].Y();
        posNext=plPosition[plNext].Y();
        posPrev4=plPosition[plPrev4].Y();
        posNext4=plPosition[plNext4].Y();
      }
      else if (pl%2==1){//odd planes, measures u and x
        //get position
        posPrev=plPosition[plPrev].X();
        posNext=plPosition[plNext].X();
        posPrev4=plPosition[plPrev4].X();
        posNext4=plPosition[plNext4].X();
      }
      
      if (posPrev==initValue || posNext==initValue){
        MAXMSG("CDAnalysis",Msg::kDebug,500)
          <<"pl="<<pl<<", not found position"
          <<", prev="<<posPrev<<", next="<<posNext
          <<endl;
      }
      else{
        //calc the new position
        pos=(posPrev+posNext)*0.5;
        //MAXMSG("CDAnalysis",Msg::kDebug,1000)
        //<<"pl="<<pl
        //<<", found positions, prev="<<posPrev<<", next="<<posNext
        //<<", pos="<<pos<<endl;
      }

      if (posPrev4==initValue || posNext4==initValue){
        MAXMSG("CDAnalysis",Msg::kDebug,500)
          <<"pl="<<pl<<", not found position"
          <<", prev4="<<posPrev4<<", next4="<<posNext
          <<endl;
      }
      else{
        //calc the new position
        pos4=(posPrev4+posNext4)*0.5;
        //MAXMSG("CDAnalysis",Msg::kDebug,1000)
        //<<"pl="<<pl
        //<<", found positions, prev4="<<posPrev4<<", next4="<<posNext4
        //<<", pos4="<<pos4<<endl;
      }
    
      Float_t posToUse=pos;
      if (posToUse==initValue) posToUse=pos4;

      if (posToUse!=initValue){
        if (pl%2==0){//even planes, measures v and y
          //set position
          //MAXMSG("CDAnalysis",Msg::kDebug,1000)
          //<<"pl="<<pl<<", setting pos="<<posToUse
          //<<", y was="<<plPosition[pl].Y()<<endl;
          plPosition[pl].SetY(posToUse);
        }
        else if (pl%2==1){//odd planes, measures u and x
          //set position
          //MAXMSG("CDAnalysis",Msg::kDebug,1000)
          //<<"pl="<<pl<<", setting pos="<<posToUse
          //<<", x was="<<plPosition[pl].X()<<endl;
          plPosition[pl].SetX(posToUse);
        }
      }
      else{
        MAXMSG("CDAnalysis",Msg::kDebug,3000)
          <<"pl="<<pl<<", fLastPlane="<<fLastPlane
          <<", lastPl="<<lastPlane
          <<", not managed to fill gap with first method"<<endl;
        unFilledGap=true;
      }
    }
    else{
      MAXMSG("CDAnalysis",Msg::kDebug,3000)
        <<"pl="<<pl<<", no gap"<<endl;
    }
  }

  //fill any gaps that failed to get filled above
  //simple algorithm to just use the previous x/y value in the track
  Float_t lastX=initValue;
  Float_t lastY=initValue;
  if (unFilledGap){
    for (Int_t pl=0;pl<=lastPlane;pl++){
      if (pl%2==0){//even planes, measures v and y
        if (plPosition[pl].Y()==initValue){
          plPosition[pl].SetY(lastY);
          //MAXMSG("CDAnalysis",Msg::kVerbose,1000)
          //<<"pl="<<pl<<", fLastPlane="<<fLastPlane
          //<<", lastPl="<<lastPlane
          //<<", filled gap (y="<<lastY<<")"<<endl;
        }
        //get last position in track
        lastY=plPosition[pl].Y();
      }
      else if (pl%2==1){//odd planes, measures u and x
        if (plPosition[pl].X()==initValue){
          plPosition[pl].SetX(lastX);
          //MAXMSG("CDAnalysis",Msg::kVerbose,1000)
          //<<"pl="<<pl<<", fLastPlane="<<fLastPlane
          //<<", lastPl="<<lastPlane
          //<<", filled gap (x="<<lastX<<")"<<endl;
        }
        //get last position in track
        lastX=plPosition[pl].X();
      }
    }
  }
  
  //loop and just do the simple calculations 
  //+/- 1 plane to calc what is not measured
  for (Int_t pl=2;pl<=lastPlane;pl++){
    Int_t plPrev=pl-1;
    Int_t plNext=pl+1;
    Float_t posPrev=initValue;
    Float_t posNext=initValue;
    Float_t pos=initValue;
    
    //protect against going outside allowed ranges
    if (plPrev<0) {
      if (abs(plPrev%2)==0) plPrev=0;//even
      else if (abs(plPrev%2)==1) plPrev=1;//odd
    }
    if (plNext>59){
      if (plNext%2==0) plNext=58;//even
      else if (plNext%2==1) plNext=59;//odd
    }
    
    if (pl%2==0){//even planes, measures v and y
      //get position
      posPrev=plPosition[plPrev].X();
      posNext=plPosition[plNext].X();
    }
    else if (pl%2==1){//odd planes, measures u and x
      //get position
      posPrev=plPosition[plPrev].Y();
      posNext=plPosition[plNext].Y();
    }
    
    if (posPrev==initValue || posNext==initValue){
      MAXMSG("CDAnalysis",Msg::kDebug,500)
        <<"pl="<<pl<<", fLastPlane="<<fLastPlane
        <<", lastPl="<<lastPlane
        <<", not found position"
        <<", prev="<<posPrev<<", next="<<posNext
        <<endl;
    }
    else{
      //calc the new position
      pos=(posPrev+posNext)*0.5;
      MAXMSG("CDAnalysis",Msg::kDebug,1000)
        <<"pl="<<pl
        <<", found positions, prev="<<posPrev<<", next="<<posNext
        <<", pos="<<pos<<endl;
    }
    
    if (pl%2==0){//even planes, measures v and y
      //set position
      plPosition[pl].SetX(pos);
    }
    else if (pl%2==1){//odd planes, measures u and x
      //set position
      plPosition[pl].SetY(pos);
    }
  }
  
  //need to do a 3D vector type correction to work out the last plane
  //but hack for now - use same x/y as previous plane
  Int_t pl=lastPlane;
  if (pl%2==0){//even planes, measures v and y
    //set position
    Float_t lastPlX=plPosition[pl-1].X();
    plPosition[pl].SetX(lastPlX);
  }
  else if (pl%2==1){//odd planes, measures u and x
    //set position
    Float_t lastPlY=plPosition[pl-1].Y();
    plPosition[pl].SetY(lastPlY);
  }
  
  //MAXMSG("CDAnalysis",Msg::kInfo,200)
  //<<"fLastPlane="<<fLastPlane
  //<<", lastPl="<<lastPlane
  //<<endl;
  //loop and print positions
  //for (Int_t pl=0;pl<60;pl++){
  //MAXMSG("CDAnalysis",Msg::kInfo,5000)
  //  <<"pl="<<pl<<" (x,y,z)=("<<plPosition[pl].X()
  //  <<","<<plPosition[pl].Y()
  //  <<","<<plPosition[pl].Z()<<")"<<endl;
  //}
  
  //make map to return (not nice)
  for (Int_t pl=0;pl<60;pl++){
    mPlPos[pl]=plPosition[pl];
  }
  
  return true;
}

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

Float_t CDAnalysis::CalcDistToPlaneCentre(Int_t strip,
                                          Float_t transPos) const
{
  Float_t stripSize=0.041*Munits::m;
  Float_t distToPlaneCentre=-1;
  //note that strips are numbered from 0  

  //cutting on transPos==0 is ugly since it's an allowed value
  //but this is what it is initialised to and it 
  //should only cut real values rarely
  if (strip<0 || strip>23 || transPos<0 || transPos==0 || 
      transPos>23) return distToPlaneCentre;

  Float_t dist1=-1;
  Float_t dist2=-1;
  Float_t stripsFromCentre1=-1;
  Float_t transPosFromCentre1=-1;

  //get distances measured from plane centre in each view
  if (strip<=11) stripsFromCentre1=11-strip;
  else if (strip>11) stripsFromCentre1=strip-12;
  if (transPos<=11) transPosFromCentre1=11-transPos;
  else if (transPos>11) transPosFromCentre1=transPos-11;

  //calc real distances
  dist1=stripsFromCentre1*stripSize;
  dist2=transPosFromCentre1*stripSize;
  
  //move it to the middle of a strip so 
  //you know which strip it's actually in
  dist1+=stripSize*0.5;
  //transPos is not quantised so don't worry about it here

  //simple bit of trig
  if (dist1>0 && dist2>0){
    distToPlaneCentre=sqrt(pow(dist1,2)+pow(dist2,2));
  }
  
  return distToPlaneCentre;
}

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

Float_t CDAnalysis::TrueDistToPlaneCentre(Int_t strip,
                                          Float_t x) const
{
  Float_t stripSize=0.041*Munits::m;
  Float_t distToPlaneCentre=-1;
  //note that strips are numbered from 0  

  if (strip<0 || strip>23 || x>0.5 || x<-0.5) return distToPlaneCentre;

  Float_t dist1=-1;
  Float_t dist2=-1;
  Float_t stripsFromCentre1=-1;

  //get distances measured from plane centre in each view
  if (strip<=11) stripsFromCentre1=11-strip;
  else if (strip>11) stripsFromCentre1=strip-12;
  if (x<0) x=x-2*x;//convert to positive for later check

  //calc real distances
  dist1=stripsFromCentre1*stripSize;
  dist2=x;//don't need to do it for x  

  //move it to the middle of a strip so 
  //you know which strip it's actually in
  dist1+=stripSize*0.5;

  //simple bit of trig
  if (dist1>0 && dist2>0){
    distToPlaneCentre=sqrt(pow(dist1,2)+pow(dist2,2));
  }
  
  return distToPlaneCentre;
}

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

void CDAnalysis::StandardSanityChecks() const
{
  //new check to not allow any capid errors
  if (fChargeAdc<1){
    MSG("CDAnalysis",Msg::kFatal)
      <<"Run has bad adcs!"<<endl
      <<"adc="<<fChargeAdc<<" ("<<fPlane<<":"<<fStrip<<":"<<fStripend
      <<")"<<endl;
    assert(false);
  }
}

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

void CDAnalysis::ScaleVector(std::vector<Double_t>& v,
                             Double_t scaleFactor) const
{
  MSG("CDAnalysis",Msg::kDebug) 
    <<" ** Running ScaleVector method... **"<<endl;

  for (UInt_t i=0;i<v.size();i++) v[i]*=scaleFactor;

  MSG("CDAnalysis",Msg::kDebug) 
    <<" ** Finished ScaleVector method **"<<endl;
}

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

void CDAnalysis::FlipVector(std::vector<Double_t>& v) const
{
  MSG("CDAnalysis",Msg::kDebug) 
    <<" ** Running FlipVector method... **"<<endl;

  vector<Double_t> tempV=v;

  for (UInt_t i=0;i<v.size();i++) v[i]=tempV[v.size()-1-i];

  MSG("CDAnalysis",Msg::kDebug) 
    <<" ** Finished FlipVector method **"<<endl;
}

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

Int_t CDAnalysis::Rnd(Double_t x) const
{
  Double_t xi=0;

  Double_t floatPart=modf(x,&xi); 

  //round up if necessary
  if (floatPart>0.5) xi++;

  return static_cast<Int_t>(xi);
}

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

void CDAnalysis::ScaleMap(std::map<Int_t,Float_t>& m,
                          const std::map<Int_t,Float_t>& scaleFactor) const
{
  MSG("CDAnalysis",Msg::kDebug) 
    <<" ** Running ScaleMap method... **"<<endl;

  if (m.size()==scaleFactor.size()){
    
    map<Int_t,Float_t>::iterator mIter=m.begin();
    map<Int_t,Float_t>::const_iterator scale=scaleFactor.begin();

    //loop over map and scale it
    while (mIter!=m.end()){

      //divide by the scale factor
      if (scale->second!=0){
        MSG("CDAnalysis",Msg::kVerbose)
          <<"mIter="<<mIter->second<<", scale="<<scale->second<<endl;
        mIter->second/=scale->second;
      }
      else MSG("CDAnalysis",Msg::kWarning)<<"Scale factor zero!"<<endl;

      mIter++;
      scale++;
    }
  }
  else{
    MSG("CDAnalysis",Msg::kWarning)
      <<"Not scaling map, vectors are different sizes!"<<endl;
  }

  MSG("CDAnalysis",Msg::kDebug) 
    <<" ** Finished ScaleMap method **"<<endl;
}

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

Bool_t CDAnalysis::CutOnFiducialVolume() const
{
  Bool_t cutThisOne=false;

  Int_t stripsFromCentreCut=10;

  if (fStripsFromCentrePeCut>=stripsFromCentreCut) cutThisOne=true;

  if (cutThisOne){
    MAXMSG("CDAnalysis",Msg::kInfo,1)
      <<"Applying cut on fiducial volume:"<<endl
      <<" stripsFromCentre<"<<stripsFromCentreCut
      <<" (fStripsFromCentrePeCut="<<fStripsFromCentrePeCut<<")"
      <<endl;
  }
  
  return cutThisOne;
}

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

Bool_t CDAnalysis::CutOnNumHitsPerPlane() const
{
  Bool_t cutThisOne=false;

  Float_t hitsPerPlLowCut=1.7;
  Float_t hitsPerPlCut=2.9;

  if (fNumHitsPerPlane<hitsPerPlLowCut ||
      fNumHitsPerPlane>hitsPerPlCut) cutThisOne=true;

  if (cutThisOne){
    MAXMSG("CDAnalysis",Msg::kInfo,1)
      <<"Applying cut on hits per plane:"
      <<" numHitsPerPlane<"<<hitsPerPlCut<<" and >"<<hitsPerPlLowCut
      <<" (fNumHitsPerPlane="<<fNumHitsPerPlane<<")"<<endl
      <<endl;
  }
  
  return cutThisOne;
}

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

Bool_t CDAnalysis::CutOnNumHitsPerPlane10() const
{
  Bool_t cutThisOne=false;

  Float_t hitsPerPlPeCut10LowCut=1.1;
  Float_t hitsPerPlPeCut10Cut=2.7;

  if (fNumHitsPerPlanePeCut10<hitsPerPlPeCut10LowCut ||
      fNumHitsPerPlanePeCut10>hitsPerPlPeCut10Cut) cutThisOne=true;

  if (cutThisOne){
    MAXMSG("CDAnalysis",Msg::kInfo,1)
      <<"Applying cut on hits per plane:"
      <<" numHitsPerPlanePeCut10<"<<hitsPerPlPeCut10Cut
      <<" and >"<<hitsPerPlPeCut10LowCut
      <<" (fNumHitsPerPlanePeCut10="<<fNumHitsPerPlanePeCut10<<")"
      <<endl;
  }
  
  return cutThisOne;
}

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

Bool_t CDAnalysis::CutOnDeadChips() const
{
  Bool_t cutThisOne=false;

  //check if dead chips
  if (fNumDeadChips>0) cutThisOne=true;

  if (cutThisOne){
    MAXMSG("CDAnalysis",Msg::kInfo,1)
      <<"Applying cut on dead chips:"
      <<" numDeadChips="<<fNumDeadChips<<endl;
  }
  
  return cutThisOne;
}

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

Bool_t CDAnalysis::CutOnBadPedestals(Int_t plane,Int_t strip,
                                     Int_t end) const
{
  //only have bad channels in real data so return if not data
  if (fSimFlag!=SimFlag::kData) return false;

  Bool_t cutThisOne=false;

  //cut out the strips with bad pedestals
  //not sure exactly for which run number the problem was fixed
  //had gone away by 71266
  if (fRunNumber>=70000 && fRunNumber<=70888){
    if ((plane==21 || plane==41) && strip==19) cutThisOne=true;
  }
  else if (fRunNumber>40000 && fRunNumber<60000){
    //cut out the strips with bad pedestals
    if ((plane==21 || plane==41) && strip==19 && end==1){
      cutThisOne=true;
    }
    if ((plane==20 || plane==40) && strip==4 && end==2){
      cutThisOne=true;
    }
  }    
  
  //not sure about this, MK reported it bad
  //in the 30000s and 40000s (i.e. T11) 38,10,2 was bad
  
  if (cutThisOne){
    MAXMSG("CDAnalysis",Msg::kInfo,6)
      <<"Applying cut to bad pedestal channel:"
      <<" (pl="<<plane<<", str="<<strip<<", end="<<end<<")"<<endl;
  }
  
  return cutThisOne;
}

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

Bool_t CDAnalysis::CutOnBadCalorimetry(Int_t plane,Int_t strip,
                                       Int_t end) const
{
  //only have bad channels in real data so return if not data
  if (fSimFlag!=SimFlag::kData) return false;

  Bool_t cutThisOne=false;

  //cut out the strips with bad calorimetry
  //do this for all run periods to be consistent
  if (fPlane==0 || fPlane==35) cutThisOne=true;
  
  if (cutThisOne){
    MAXMSG("CDAnalysis",Msg::kInfo,6)
      <<"Applying cut to bad calorimetry:"
      <<" (pl="<<plane<<", str="<<strip<<", end="<<end<<")"<<endl;
  }
  
  return cutThisOne;
}

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

Bool_t CDAnalysis::CutOnTrackQuality(Bool_t requireGoodTrack) const
{
  //set default to not cut
  Bool_t cutThisOne=false;

  //look if a good first plane was found
  //this is a good indication of tracking success
  Bool_t goodFirstPlane=(fFirstPlane>=0 && fFirstPlane<59);
  if (requireGoodTrack){
    if (!goodFirstPlane) cutThisOne=true;
  }

  //don't do this as not necessary for MEU calibration
  //cut on the trackers opinion of what's a good track
  //Bool_t goodInBothViews=(fResultEven==1 && fResultOdd==1);
  //if (requireGoodTrack){
  //  if (!goodInBothViews) cutThisOne=true;
  //}

  if (cutThisOne){
    MAXMSG("CDAnalysis",Msg::kInfo,1)
      <<"Applying cut to track quality"<<endl;
  }

  return cutThisOne;
}

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

Bool_t CDAnalysis::CutOnOverlap() const
{
  //set default to not cut
  Bool_t cutThisOne=false;

  //only have pid for real data so return if not data
  if (fSimFlag!=SimFlag::kData) return false;

  //if (!fNoOverlap){ //not no overlap=overlap
  //cutThisOne=true;
  //}
  Float_t olChi2Threshold=1.1;
  if (fOLChi2>olChi2Threshold) cutThisOne=true;

  if (cutThisOne){
    MAXMSG("CDAnalysis",Msg::kInfo,1)
      <<"Applying cut on overlap, OLChi2<"<<olChi2Threshold
      <<endl;
  }

  return cutThisOne;
}

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

Bool_t CDAnalysis::CutOnOverlapForElec() const
{
  //set default to not cut
  Bool_t cutThisOne=false;

  //only have pid for real data so return if not data
  if (fSimFlag!=SimFlag::kData) return false;

  Float_t olChi2Threshold=1.1;
  if (fOLChi2>olChi2Threshold) cutThisOne=true;

  if (cutThisOne){
    MAXMSG("CDAnalysis",Msg::kInfo,1)
      <<"Applying cut on overlap for elec, OLChi2<"<<olChi2Threshold
      <<endl;
  }

  return cutThisOne;
}

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

Bool_t CDAnalysis::CutOnPidForElec() const
{
  //set default to not cut
  Bool_t cutThisOne=false;

  //only have pid for real data so return if not data
  if (fSimFlag!=SimFlag::kData) return false;

  if (fPIDType!=CalDetParticleType::kElectron){ 
    cutThisOne=true;
  }
  if (!fInCERTime){ 
    cutThisOne=true;
  }

  if (cutThisOne){
    MAXMSG("CDAnalysis",Msg::kInfo,1)
      <<"Applying cut on pid to pass only electrons"<<endl;
  }

  return cutThisOne;
}

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

Bool_t CDAnalysis::CutByHandOnPidForElec() const
{
  //set default to not cut
  Bool_t cutThisOne=false;

  //only have pid for real data so return if not data
  if (fSimFlag!=SimFlag::kData) return false;

  string particleSelection="for electrons";
  MAXMSG("CDAnalysis",Msg::kInfo,1)
    <<"Applying cut: "<<particleSelection<<endl;

  Bool_t elecKovADC1=false;
  Bool_t elecKovADC3=false;
  Bool_t goodTof=false;
  Bool_t elecInTime=false;

  if (fRunNumber==70579){//a 1 GeV run

    if (fTofTDC2-fTofTDC0>340 && fTofTDC2-fTofTDC0<500) goodTof=true;

    //don't analyse if not a elec
    if (!elecKovADC1) cutThisOne=true;
    if (!elecKovADC3) cutThisOne=true;
    
    //don't analyse if bad tof
    if (!goodTof) cutThisOne=true;
  }
  else if (fRunNumber==70709){
    MAXMSG("CDAnalysis",Msg::kInfo,1)
      <<"Applying cut: "<<particleSelection<<endl;

    if (fKovADC3>2000 && fKovADC3<3500) elecKovADC3=true;//elec ADC
    if (fKovADC1>1000 && fKovADC1<2500) elecKovADC1=true;
    if (fTofTDC2-fTofTDC0>385 && fTofTDC2-fTofTDC0<450) goodTof=true;
    //Tof in time with elec triggering kov
    if ((fTofADCTimeStamp1-fKovTimeStamp3)<-80 &&
        (fTofADCTimeStamp1-fKovTimeStamp3)>-85 &&
        (fTofADCTimeStamp1-fKovTimeStamp1)<-30 &&
        (fTofADCTimeStamp1-fKovTimeStamp1)>-40) elecInTime=true;
    
    //don't analyse if not a elec
    if (!elecKovADC1) cutThisOne=true;
    if (!elecKovADC3) cutThisOne=true;
    
    //don't analyse if bad tof
    if (!goodTof) cutThisOne=true;
    
    //don't analyse if tof not in time with kov set to elecs
    if (!elecInTime) cutThisOne=true;

  }
  else if (fRunNumber==71266){

  }
  else if (fRunNumber==100131){

  }
  else if (fRunNumber==100161){

    MAXMSG("CDAnalysis",Msg::kInfo,1)
      <<"Applying cut: "<<particleSelection<<endl;

    if (fKovADC3>1600 && fKovADC3<2100) elecKovADC3=true;//elec ADC
    if (fKovADC1>1100 && fKovADC1<3100) elecKovADC1=true;
    if (fTofTDC2-fTofTDC0>385 && fTofTDC2-fTofTDC0<470) goodTof=true;
    //Tof in time with elec triggering kov
    if ((fTofADCTimeStamp1-fKovTimeStamp3)<-20 &&
        (fTofADCTimeStamp1-fKovTimeStamp3)>-30 &&
        (fTofADCTimeStamp1-fKovTimeStamp1)<-50 &&
        (fTofADCTimeStamp1-fKovTimeStamp1)>-62) elecInTime=true;
    
    //don't analyse if not a elec
    if (!elecKovADC1) cutThisOne=true;
    if (!elecKovADC3) cutThisOne=true;
    
    //don't analyse if bad tof
    if (!goodTof) cutThisOne=true;
    
    //don't analyse if tof not in time with kov set to elecs
    if (!elecInTime) cutThisOne=true;


  }
  else if (fRunNumber==100243){

    
  }
  else if (fRunNumber==100270){

  }
  
  if (cutThisOne){
    MAXMSG("CDAnalysis",Msg::kInfo,1)
      <<"Applying 'by hand' cut on pid to pass only electrons"<<endl;
  }

  return cutThisOne;
}

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

Bool_t CDAnalysis::CutOnPid() const
{
  //set default to not cut
  Bool_t cutThisOne=false;

  //only have pid for real data so return if not data
  if (fSimFlag!=SimFlag::kData) return false;

  Bool_t muonOrPion=fPIDType & (CalDetParticleType::kMuon |
                                CalDetParticleType::kPion);

  if (!muonOrPion){ 
    cutThisOne=true;
  }
  if (!fInCERTime){ 
    cutThisOne=true;
  }

  if (cutThisOne){
    MAXMSG("CDAnalysis",Msg::kInfo,1)
      <<"Applying cut on pid to pass only pions and muons"<<endl;
  }

  return cutThisOne;
}

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

Bool_t CDAnalysis::CutByHandOnPid() const
{
  //set default to not cut
  Bool_t cutThisOne=false;

  //only have pid for real data so return if not data
  if (fSimFlag!=SimFlag::kData) return false;

  Bool_t muonKovADC1=false;
  Bool_t muonKovADC3=false;
  Bool_t goodTof=false;
  Bool_t muonInTime=false;

  if (fRunNumber==70579){//a 1 GeV run
    if (fKovADC3<500) muonKovADC3=true;
    if (fKovADC1<300) muonKovADC1=true;//not an electron ADC
    if (fTofTDC2-fTofTDC0>340 && fTofTDC2-fTofTDC0<500) goodTof=true;

    //don't analyse if not a muon
    if (!muonKovADC1) cutThisOne=true;
    if (!muonKovADC3) cutThisOne=true;
    
    //don't analyse if bad tof
    if (!goodTof) cutThisOne=true;
  }
  else if (fRunNumber==70709){
    //string particleSelection="muons";
    string particleSelection="notElectrons";
    //string particleSelection="notElectronsNoCkv";

    if (particleSelection=="muons"){
      MAXMSG("CDAnalysis",Msg::kWarning,1)
        <<"Applying cut: "<<particleSelection<<endl;
      
      if (fKovADC3<1400 && fKovADC3>0) muonKovADC3=true;//muon ADC
      if (fKovADC1<500) muonKovADC1=true;//not an electron ADC
      if (fTofTDC2-fTofTDC0>385 && fTofTDC2-fTofTDC0<450) goodTof=true;
      //Tof in time with muon triggering kov
      if ((fTofADCTimeStamp1-fKovTimeStamp3)<-78 &&
          (fTofADCTimeStamp1-fKovTimeStamp3)>-96) muonInTime=true;
      
      //don't analyse if not a muon
      if (!muonKovADC1) cutThisOne=true;
      if (!muonKovADC3) cutThisOne=true;
      
      //don't analyse if bad tof
      if (!goodTof) cutThisOne=true;
      
      //don't analyse if tof not in time with kov set to muons
      if (!muonInTime) cutThisOne=true;
    }
    else if (particleSelection=="notElectrons"){
      MAXMSG("CDAnalysis",Msg::kWarning,1)
        <<"Applying cut: "<<particleSelection<<endl;
      
      if (fKovADC3<1400) muonKovADC3=true;//not an electron ADC
      if (fKovADC1<500) muonKovADC1=true;//not an electron ADC
      if (fTofTDC2-fTofTDC0>385 && fTofTDC2-fTofTDC0<450) goodTof=true;
      
      //don't analyse if not a muon
      if (!muonKovADC1) cutThisOne=true;
      if (!muonKovADC3) cutThisOne=true;
      
      //don't analyse if bad tof
      if (!goodTof) cutThisOne=true;

      //don't want the muon in time cut
    }
  }
  else if (fRunNumber==71266){
    if (fKovADC3<1400) muonKovADC3=true;//not an electron ADC
    if (fKovADC1<1000) muonKovADC1=true;//not an electron ADC
    if (fTofTDC2-fTofTDC0>385 && fTofTDC2-fTofTDC0<450) goodTof=true;

    //don't analyse if not a muon
    if (!muonKovADC1) cutThisOne=true;
    if (!muonKovADC3) cutThisOne=true;
    
    //don't analyse if bad tof
    if (!goodTof) cutThisOne=true;
  }
  else if (fRunNumber==100131){
    if (fKovADC3<1000) muonKovADC3=true;
    if (fKovADC1<500) muonKovADC1=true;//not an electron ADC
    if (fTofTDC2-fTofTDC0>370 && fTofTDC2-fTofTDC0<525) goodTof=true;

    //don't analyse if not a muon
    if (!muonKovADC1) cutThisOne=true;
    if (!muonKovADC3) cutThisOne=true;
    
    //don't analyse if bad tof
    if (!goodTof) cutThisOne=true;
  }
  else if (fRunNumber==100161){

    //string particleSelection="muons";
    string particleSelection="notElectrons";
    //string particleSelection="notElectronsNoCkv";

    if (particleSelection=="muons"){
      MAXMSG("CDAnalysis",Msg::kWarning,1)
        <<"Applying cut: "<<particleSelection<<endl;

      if (fKovADC3<850 && fKovADC3>0) muonKovADC3=true;
      if (fKovADC1<400) muonKovADC1=true;//not an electron ADC
      if (fTofTDC2-fTofTDC0>385 && fTofTDC2-fTofTDC0<490) goodTof=true;
      //Tof in time with muon triggering kov
      if ((fTofADCTimeStamp1-fKovTimeStamp3)<-24 &&
          (fTofADCTimeStamp1-fKovTimeStamp3)>-45) muonInTime=true;
      
      //don't analyse if not a muon
      if (!muonKovADC1) cutThisOne=true;
      if (!muonKovADC3) cutThisOne=true;
      
      //don't analyse if bad tof
      if (!goodTof) cutThisOne=true;
      
      //don't analyse if tof not in time with kov set to muons
      if (!muonInTime) cutThisOne=true;
    }
    else if (particleSelection=="notElectrons"){
      MAXMSG("CDAnalysis",Msg::kWarning,1)
        <<"Applying cut: "<<particleSelection<<endl;

      if (fKovADC3<850) muonKovADC3=true;
      if (fKovADC1<400) muonKovADC1=true;//not an electron ADC
      if (fTofTDC2-fTofTDC0>385 && fTofTDC2-fTofTDC0<490) goodTof=true;

      //don't analyse if not a muon
      if (!muonKovADC1) cutThisOne=true;
      if (!muonKovADC3) cutThisOne=true;
      
      //don't analyse if bad tof
      if (!goodTof) cutThisOne=true;

      //don't want the muon in time cut
    }
    else if (particleSelection=="notElectronsNoCkv"){
      MAXMSG("CDAnalysis",Msg::kWarning,1)
        <<"Applying cut: "<<particleSelection<<endl;

      if (fKovADC3<100) muonKovADC3=true;
      if (fKovADC1<100) muonKovADC1=true;//not an electron ADC
      if (fTofTDC2-fTofTDC0>385 && fTofTDC2-fTofTDC0<490) goodTof=true;

      //don't analyse if not a muon
      if (!muonKovADC1) cutThisOne=true;
      if (!muonKovADC3) cutThisOne=true;
      
      //don't analyse if bad tof
      if (!goodTof) cutThisOne=true;

      //don't want the muon in time cut
    }
  }
  else if (fRunNumber==100243){
    if (fKovADC3<900) muonKovADC3=true;
    if (fKovADC1<1000) muonKovADC1=true;//not an electron ADC
    if (fTofTDC2-fTofTDC0>350 && fTofTDC2-fTofTDC0<490) goodTof=true;

    //don't analyse if not a muon
    if (!muonKovADC1) cutThisOne=true;
    if (!muonKovADC3) cutThisOne=true;
    
    //don't analyse if bad tof
    if (!goodTof) cutThisOne=true;
    
  }
  else if (fRunNumber==100270){
    if (fKovADC3<1100) muonKovADC3=true;
    if (fKovADC1<500) muonKovADC1=true;//not an electron ADC
    if (fTofTDC2-fTofTDC0>360 && fTofTDC2-fTofTDC0<525) goodTof=true;

    //don't analyse if not a muon
    if (!muonKovADC1) cutThisOne=true;
    if (!muonKovADC3) cutThisOne=true;
    
    //don't analyse if bad tof
    if (!goodTof) cutThisOne=true;    
  }

  if (cutThisOne){
    MAXMSG("CDAnalysis",Msg::kInfo,1)
      <<"Applying 'by hand' cut on pid"<<endl;
  }

  return cutThisOne;
}

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

Bool_t CDAnalysis::CutOnScint_TofDiff() const
{
  //set default to not cut
  Bool_t cutThisOne=false;

  Float_t lowTime=-1000000;
  Float_t upTime=1000000;
  
  this->GetLowUpTimeCuts(lowTime,upTime);
  
  if (fLowScint_Tof<lowTime || fUpScint_Tof>upTime) cutThisOne=true;

  if (cutThisOne){
    MAXMSG("CDAnalysis",Msg::kInfo,10)
      <<"Using TOF-scint diff. cut: "
      <<lowTime<<" -> "<<upTime<<" ns"
      <<" (t="<<fLowScint_Tof<<" -> "<<fUpScint_Tof<<" ns)"
      <<endl;
  }

  return cutThisOne;
}

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

void CDAnalysis::GetLowUpTimeCuts(Float_t& lowTime,
                                  Float_t& upTime) const
{
  //this works in ns rather naughtily!
  lowTime=-1000000;
  upTime=1000000;

  Int_t run=fRunNumber;

  //set the time gate here so that you get the same for each run 
  //period
  Int_t lowTimeLimit=35;
  Int_t upTimeLimit=295;
  Int_t peak=-1;
  
  //only have timing issues for real data
  if (fSimFlag==SimFlag::kData) {
    if (run>=40000 && run<50000){
      //just for interest: the peak is 180 ns into the 600 ns varc gate
      peak=-125;//ns
      lowTime=peak-lowTimeLimit;//was -175;
      upTime=peak+upTimeLimit;//was 135;
    }
    else if (run>=50000 && run<60000){
      //just for interest: the peak is 280 ns into the 600 ns varc gate
      peak=30;//ns
      lowTime=peak-lowTimeLimit;
      upTime=peak+upTimeLimit;
    }
    else if (run>=70000 && run<80000){
      //just for interest: the peak is 220 ns into the 600 ns varc gate
      peak=5;//ns
      lowTime=peak-lowTimeLimit;//-30;//-50;
      upTime=peak+upTimeLimit;//300;//260;
    }
    else{
      MAXMSG("CDAnalysis",Msg::kInfo,10)
        <<"No time cuts cut for this run="<<fRunNumber<<endl;
    }
  }
  else{
    MAXMSG("CDAnalysis",Msg::kInfo,1)
      <<"MC data so no sensible timing cuts to set"<<endl;
  }

  MAXMSG("CDAnalysis",Msg::kInfo,1)
    <<"Set time cuts:  low="<<lowTime<<" ns, upper="<<upTime<<" ns"
    <<endl;
}

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

Bool_t CDAnalysis::CutOnEventLengthForElec() const
{
  //set default to not cut
  Bool_t cutThisOne=false;

  //Int_t minPlane=0;
  Int_t maxPlane=20;

  //old cut
  //if (fLastPlane<minPlane || fLastPlane>maxPlane) cutThisOne=true;

  //new cut - number of planes hit allows for a noise hit in the
  //back of the detector but will cut out muons sucessfully
  Int_t numPlanes=fNumPlanesHitAll.size();
  if (numPlanes>maxPlane) cutThisOne=true;

  if (cutThisOne){
    MAXMSG("CDAnalysis",Msg::kInfo,1)
      <<"Using event length cuts for electrons: "
      <<"Num hit planes<="<<maxPlane<<endl;
  }

  return cutThisOne;
}

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

void CDAnalysis::GetEvLenMinMax(Int_t& minPlane,Int_t& maxPlane) const
{
  //I don't understand abs. Sometimes it just returns an int!!
  //e.g. in a root macro
  Float_t aBeamP=fabs(fBeamMomentum);
  
  //new cuts are now inclusive so it's 31<=Pl<=45

  if (fabs(fBeamMomentum)>1.3 && fabs(fBeamMomentum)<1.5){
    minPlane=31;//minPlane=35;
    maxPlane=45;//maxPlane=44;
  }
  else if (fabs(fBeamMomentum)>1.5 && fabs(fBeamMomentum)<1.7){
    minPlane=35;//minPlane=40;
    maxPlane=51;//maxPlane=49;
  }
  else if (fabs(fBeamMomentum)>1.7 && fabs(fBeamMomentum)<1.9){
    minPlane=39;//minPlane=42;
    maxPlane=56;//maxPlane=57;
  }
  else if (fabs(fBeamMomentum)>1.9){
    minPlane=44;//minPlane=44;
    maxPlane=57;//maxPlane=58;
  }
  else{
    MAXMSG("CDAnalysis",Msg::kInfo,10)
      <<"No event length cut for this particle energy="
      <<fBeamMomentum<<endl;
    assert(false);
  }
  MAXMSG("CDAnalysis",Msg::kInfo,1)
    <<"Set Ev. Len min="<<minPlane<<", max="<<maxPlane
    <<", aBeamP="<<aBeamP<<endl;
}

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

Bool_t CDAnalysis::CutOnNumPlanesHit(Int_t planesHit) const
{
  //set default to not cut
  Bool_t cutThisOne=false;

  Int_t minPlane=0;
  Int_t maxPlane=60;

  this->GetEvLenMinMax(minPlane,maxPlane);

  //new cuts are now inclusive so it's 31<=Pl<=45
  if (planesHit<minPlane || planesHit>maxPlane) cutThisOne=true;

  if (cutThisOne){
    MAXMSG("CDAnalysis",Msg::kInfo,1)
      <<"Using num planes hit cut to pass muons: "
      <<minPlane<<"<=Plane<="<<maxPlane
      <<", fBeamMomentum="<<fBeamMomentum<<endl;
  }

  MSG("CDAnalysis",Msg::kDebug)
    <<"fLastPlane="<<fLastPlane<<", cut this one="<<cutThisOne<<endl;

  return cutThisOne;
}

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

Bool_t CDAnalysis::CutOnPSMuons() const
{
  //set default to not cut
  Bool_t cutThisOne=false;

  if (fAvStrip<9|| fAvStrip>13){
    cutThisOne=true;
  }

  return cutThisOne;
}

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

Int_t CDAnalysis::GetEventLength() const
{
  Int_t lastPlane=fLastPlane;
  
  //if there was no track then a last plane is not found
  //in this case use the number of planes hit
  //this will be fairly rare for muons
  if (lastPlane==-1) {
    MAXMSG("CDAnalysis",Msg::kInfo,1)
      <<"GetEventLength::Using num planes hit in event length:"
      <<" fLastPlane="<<fLastPlane
      <<", numPlanesHit="<<fNumPlanesHitAll.size()<<endl;
    lastPlane=fNumPlanesHitAll.size();
  }
  
  return lastPlane;
}

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

Bool_t CDAnalysis::CutOnEventLength() const
{
  //set default to not cut
  Bool_t cutThisOne=false;

  Int_t lastPlane=this->GetEventLength();
  
  Int_t minPlane=0;
  Int_t maxPlane=60;
  this->GetEvLenMinMax(minPlane,maxPlane);

  //new cuts are now inclusive so it's 31<=Pl<=45
  if (lastPlane<minPlane || lastPlane>maxPlane) cutThisOne=true;

  if (cutThisOne){
    MAXMSG("CDAnalysis",Msg::kInfo,1)
      <<"Using event length cuts to pass muons: "
      <<minPlane<<"<=Plane<="<<maxPlane
      <<" (lastPlane="<<lastPlane<<")"
      <<", fBeamMomentum="<<fBeamMomentum<<endl;
  }

  MSG("CDAnalysis",Msg::kDebug)
    <<"fLastPlane="<<fLastPlane<<", cut this one="<<cutThisOne<<endl;

  return cutThisOne;
}

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

Bool_t CDAnalysis::CutByHandOnEventLength() const
{
  //set default to not cut
  Bool_t cutThisOne=false;

  Int_t minPlane=0;
  Int_t maxPlane=60;

  if (fSimFlag==SimFlag::kData){
    
    //might want to tighten these up to boost the signal to noise ratio
    
    if (fRunNumber==70709 ||
        fRunNumber==71266 ||
        fRunNumber==100161 ||
        fRunNumber==100270){//1.8
      //cut out through going muons, no 58,59
      minPlane=42;
      maxPlane=57;
      if (fLastPlane<minPlane || fLastPlane>maxPlane) cutThisOne=true;
    }
    else if (fRunNumber==100243){//1.6
      minPlane=40;
      maxPlane=49;
      if (fLastPlane<minPlane || fLastPlane>maxPlane) cutThisOne=true;
    }
    else if (fRunNumber==100131){//1.4
      minPlane=35;
      maxPlane=44;
      if (fLastPlane<minPlane || fLastPlane>maxPlane) cutThisOne=true;
    }
    else if (fRunNumber==70579){//1.0 GeV run
      //hack - these are 1.8 cuts
      if (fLastPlane<42 || fLastPlane>57) cutThisOne=true;
    }
  }
  else if (fSimFlag==SimFlag::kMC){
    
    if (fMainParticleEnergy>1.3 && fMainParticleEnergy<1.5){
      minPlane=35;
      maxPlane=44;
      if (fLastPlane<minPlane || fLastPlane>maxPlane) cutThisOne=true;
    }
    else if (fMainParticleEnergy>1.5 && fMainParticleEnergy<1.7){
      minPlane=40;
      maxPlane=49;
      if (fLastPlane<minPlane || fLastPlane>maxPlane) cutThisOne=true;
    }
    else if (fMainParticleEnergy>1.7 && fMainParticleEnergy<1.9){
      minPlane=42;
      maxPlane=57;
      if (fLastPlane<minPlane || fLastPlane>maxPlane) cutThisOne=true;
    }
    else if (fMainParticleEnergy>1.9){
      minPlane=44;
      maxPlane=58;//last but one, so didn't exit
      if (fLastPlane<minPlane || fLastPlane>maxPlane) cutThisOne=true;
    }
    else{
      MAXMSG("CDAnalysis",Msg::kInfo,10)
        <<"No event length cut for this particle energy="
        <<fMainParticleEnergy<<endl;
    }
  }

  if (cutThisOne){
    MAXMSG("CDAnalysis",Msg::kInfo,1)
      <<"Using 'by hand' event length cuts to pass muons: "
      <<minPlane<<"<Plane<"<<maxPlane
      <<", fBeamMomentum="<<fBeamMomentum<<endl;
  }

  MSG("CDAnalysis",Msg::kDebug)
    <<"fLastPlane="<<fLastPlane<<", cut this one="<<cutThisOne<<endl;

  return cutThisOne;
}

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

void CDAnalysis::FillProfHisto(TProfile* prof,
                               std::map<Int_t,Float_t>& counter,
                               const std::map<Int_t,Float_t>& addMe,
                               Int_t offset) const
{
  if (prof){
    for (map<Int_t,Float_t>::const_iterator it=addMe.begin();
         it!=addMe.end();++it){
      Int_t key=it->first;
      if (offset!=0) key=offset-(it->first);

      prof->Fill(key,it->second);
      counter[key]++;
    }
  }
  else{
    MSG("CDAnalysis",Msg::kWarning)
      <<"TProfile not initialised"<<endl;
  }
}

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

Bool_t CDAnalysis::IsPlaneXTalkOnly(Int_t inplane) const
{
  if (!fTruthHitInfo) return false;

  TClonesArray &cTruth = *fTruthHitInfo;
  Int_t numTruthHits=fTruthHitInfo->GetEntries();
  
  Int_t sumTruth=DigiSignal::kUnknown;

  //loop over the truth hits
  for (Int_t hit=0;hit<numTruthHits;hit++){
    //cast the tclonesarray up to a truthhitinfo object
    CDTruthHitInfo *hitInfo=
      dynamic_cast<CDTruthHitInfo*>(cTruth[hit]);
    
    Int_t plane=hitInfo->GetPlane();

    if (plane!=inplane) continue;

    Int_t pmtTruth1=hitInfo->GetPmtTruth1();
    Int_t pmtTruth2=hitInfo->GetPmtTruth2();

    //add up the bits
    sumTruth|=pmtTruth1;
    sumTruth|=pmtTruth2;

    MSG("CDAnalysis",Msg::kVerbose)
      <<"sumTruth="<<sumTruth<<endl;
  }

  //check if genuine or not
  Bool_t genuineBit=((sumTruth & DigiSignal::kGenuine)==
                     DigiSignal::kGenuine);
  Bool_t optXTalkBit=((sumTruth & DigiSignal::kCrosstalkOptical)==
                      DigiSignal::kCrosstalkOptical);
  Bool_t xTalkBit=((sumTruth & DigiSignal::kCrosstalk)==
                   DigiSignal::kCrosstalk);
  
  if (genuineBit || sumTruth==0) return false;
  else if (optXTalkBit) return true;
  else if (xTalkBit) return true;
  else return false;
}

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

Bool_t CDAnalysis::IsPlaneGenuine(Int_t inplane) const
{
  if (!fTruthHitInfo) return true;

  MAXMSG("CDAnalysis",Msg::kInfo,1)
    <<"Running IsPlaneGenuine(), fTruthHitInfo="<<fTruthHitInfo<<endl;

  TClonesArray &cTruth = *fTruthHitInfo;
  Int_t numTruthHits=fTruthHitInfo->GetEntries();
  
  Int_t sumTruth=DigiSignal::kUnknown;

  //loop over the truth hits
  for (Int_t hit=0;hit<numTruthHits;hit++){
    //cast the tclonesarray up to a truthhitinfo object
    CDTruthHitInfo *hitInfo=
      dynamic_cast<CDTruthHitInfo*>(cTruth[hit]);
    
    Int_t plane=hitInfo->GetPlane();

    if (plane!=inplane) continue;

    Int_t pmtTruth1=hitInfo->GetPmtTruth1();
    Int_t pmtTruth2=hitInfo->GetPmtTruth2();

    //add up the bits
    sumTruth|=pmtTruth1;
    sumTruth|=pmtTruth2;

    MSG("CDAnalysis",Msg::kVerbose)
      <<"sumTruth="<<sumTruth<<endl;
  }

  //check if genuine or not
  Bool_t genuineBit=((sumTruth & DigiSignal::kGenuine)==
                     DigiSignal::kGenuine);
  
  if (genuineBit) return true;
  else return false;
}

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

Bool_t CDAnalysis::IsSharedPmtHit(Int_t inplane) const
{
  if (!fTruthHitInfo) return true;

  TClonesArray &cTruth = *fTruthHitInfo;
  Int_t numTruthHits=fTruthHitInfo->GetEntries();
  
  Bool_t sharedPmtHit=false;

  //loop over the truth hits
  for (Int_t hit=0;hit<numTruthHits;hit++){
    //cast the tclonesarray up to a truthhitinfo object
    CDTruthHitInfo *hitInfo=
      dynamic_cast<CDTruthHitInfo*>(cTruth[hit]);
    
    Int_t plane=hitInfo->GetPlane();
    
    if (plane!=inplane) continue;
    
    Int_t vaChip1=hitInfo->GetVaChip1();
    Int_t vaChip2=hitInfo->GetVaChip2();

    if (vaChip1==1 || vaChip2==1) sharedPmtHit=true;

    MSG("CDAnalysis",Msg::kVerbose)
      <<"sharedPmtHit="<<sharedPmtHit<<endl;
  }

  return sharedPmtHit;
}

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

Bool_t CDAnalysis::IsPlaneSideHit(Int_t inplane,Int_t stripEnd) const
{
  if (!fTruthHitInfo) return true;

  TClonesArray &cTruth = *fTruthHitInfo;
  Int_t numTruthHits=fTruthHitInfo->GetEntries();
  
  Bool_t planeSideHit=false;

  //loop over the truth hits
  for (Int_t hit=0;hit<numTruthHits;hit++){
    //cast the tclonesarray up to a truthhitinfo object
    CDTruthHitInfo *hitInfo=
      dynamic_cast<CDTruthHitInfo*>(cTruth[hit]);
    
    Int_t plane=hitInfo->GetPlane();
    
    if (plane!=inplane) continue;

    Int_t pmtTruth1=hitInfo->GetPmtTruth1();
    Int_t pmtTruth2=hitInfo->GetPmtTruth2();
    

    if (pmtTruth1 && stripEnd==StripEnd::kEast){
      planeSideHit=true;
      break;//only need the one hit
    }
    if (pmtTruth2 && stripEnd==StripEnd::kWest){
      planeSideHit=true;
      break;//only need the one hit
    }
  }

  MSG("CDAnalysis",Msg::kVerbose)
    <<"planeSideHit="<<planeSideHit<<endl;

  return planeSideHit;
}

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

Bool_t CDAnalysis::IsDoubleEnded(Int_t inplane) const
{
  if (!fTruthHitInfo) return true;

  TClonesArray &cTruth = *fTruthHitInfo;
  Int_t numTruthHits=fTruthHitInfo->GetEntries();
  
  Bool_t doubleEnded=false;

  //loop over the truth hits
  for (Int_t hit=0;hit<numTruthHits;hit++){
    //cast the tclonesarray up to a truthhitinfo object
    CDTruthHitInfo *hitInfo=
      dynamic_cast<CDTruthHitInfo*>(cTruth[hit]);
    
    Int_t plane=hitInfo->GetPlane();
    
    if (plane!=inplane) continue;
    
    Int_t pmtTruth1=hitInfo->GetPmtTruth1();
    Int_t pmtTruth2=hitInfo->GetPmtTruth2();

    if (pmtTruth1 && pmtTruth2) doubleEnded=true;

    MSG("CDAnalysis",Msg::kVerbose)
      <<"doubleEnded="<<doubleEnded<<endl;
  }

  return doubleEnded;
}

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

Bool_t CDAnalysis::IsGenuineOrXTalk(Int_t inplane) const
{
  if (!fTruthHitInfo) return true;

  TClonesArray &cTruth = *fTruthHitInfo;
  Int_t numTruthHits=fTruthHitInfo->GetEntries();
  
  Int_t sumTruth=DigiSignal::kUnknown;

  //loop over the truth hits
  for (Int_t hit=0;hit<numTruthHits;hit++){
    //cast the tclonesarray up to a truthhitinfo object
    CDTruthHitInfo *hitInfo=
      dynamic_cast<CDTruthHitInfo*>(cTruth[hit]);
    
    Int_t plane=hitInfo->GetPlane();
    
    if (plane!=inplane) continue;
    
    Int_t pmtTruth1=hitInfo->GetPmtTruth1();
    Int_t pmtTruth2=hitInfo->GetPmtTruth2();

    //add up the bits
    sumTruth|=pmtTruth1;
    sumTruth|=pmtTruth2;
  }
  
  MSG("CDAnalysis",Msg::kVerbose)
    <<"sumTruth(GenOrXTalk)="<<sumTruth<<endl;

  Bool_t genuineBit=((sumTruth & DigiSignal::kGenuine)==
                     DigiSignal::kGenuine);
  Bool_t optXTalkBit=((sumTruth & DigiSignal::kCrosstalkOptical)==
                      DigiSignal::kCrosstalkOptical);
  Bool_t xTalkBit=((sumTruth & DigiSignal::kCrosstalk)==
                   DigiSignal::kCrosstalk);

  if (genuineBit || optXTalkBit) return true;
  else if (xTalkBit){
    MSG("CDAnalysis",Msg::kDebug)
      <<" **** Only non-optXTalk, sumTruth="<<sumTruth<<endl;
    return false;
  }
  else{
    MSG("CDAnalysis",Msg::kDebug)
      <<"Not genuine or xtalk, sumTruth="<<sumTruth<<endl;
    return false;
  }
}

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

Bool_t CDAnalysis::StraightTrack_Radius(Double_t radiusCut) const
{
  //check if positive
  if (radiusCut<0) radiusCut=radiusCut-2*radiusCut;
  
  //set default return value to true
  Bool_t straight=true;

  //get tclones arrays for this snarl
  TClonesArray &cTrk=*fTrkHitInfo;
  Int_t numTrkHits=fTrkHitInfo->GetEntries();
  
  for (Int_t hit=0;hit<numTrkHits;hit++){
    CDTrackedHitInfo *trackedHitInfo=
      dynamic_cast<CDTrackedHitInfo*>(cTrk[hit]);
    
    Int_t strip=trackedHitInfo->GetStrip();
    Double_t transPos=trackedHitInfo->GetTransPos();

    //calc dist to plane centre
    Float_t radius=this->CalcDistToPlaneCentre(strip,transPos);

    //test for straight tracks
    if (radius>radiusCut){
      straight=false;
      break;
    } 
  }

  //return result
  return straight;
}

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

Bool_t CDAnalysis::IsStraightTrack_Radius(Double_t radiusCut) const
{
  if (!fTruthHitInfo) return true;

  //this method will get called once per event so use event to satisfy
  //the TrueXInStripFrame method used below
  //bit of a hack but it should work fine
  static Int_t event=0;
  event++;

  //check if positive
  if (radiusCut<0) radiusCut=radiusCut-2*radiusCut;

  TClonesArray &cTruth = *fTruthHitInfo;
  Int_t numTruthHits=fTruthHitInfo->GetEntries();
  
  //set default return value to true
  Bool_t straight=true;

  //loop over the truth hits
  for (Int_t hit=0;hit<numTruthHits;hit++){
    //cast the tclonesarray up to a truthhitinfo object
    CDTruthHitInfo *hitInfo=
      dynamic_cast<CDTruthHitInfo*>(cTruth[hit]);
    
    Int_t mainParticle=hitInfo->GetMainParticle();
    Int_t plane=hitInfo->GetPlane();
    Int_t strip=hitInfo->GetStrip();
    Int_t muon=13;

    if (abs(mainParticle)!=muon) continue;

    Float_t trueX=this->TrueXInStripFrame(plane,strip,event);   
    Float_t radiusTrue=this->TrueDistToPlaneCentre(strip,trueX);
    MAXMSG("CDAnalysis",Msg::kVerbose,300)
      <<"trueX="<<trueX<<", trueRadius="<<radiusTrue
      <<", striaght="<<straight<<", radiusCut="<<radiusCut<<endl;

    //test for straight tracks
    if (radiusTrue>radiusCut){
      straight=false;
      break;
    }
  }

  MAXMSG("CDAnalysis",Msg::kVerbose,100)
    <<"Returning straight="<<straight<<endl;

  //return result
  return straight;
}

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

Bool_t CDAnalysis::IsStraightTrack(Double_t mainX1Cut) const
{
  if (!fTruthHitInfo) return true;

  TClonesArray &cTruth = *fTruthHitInfo;
  Int_t numTruthHits=fTruthHitInfo->GetEntries();
  
  //set default return value to true
  Bool_t straight=true;

  //loop over the truth hits
  for (Int_t hit=0;hit<numTruthHits;hit++){
    //cast the tclonesarray up to a truthhitinfo object
    CDTruthHitInfo *hitInfo=
      dynamic_cast<CDTruthHitInfo*>(cTruth[hit]);
    
    Int_t mainParticle=hitInfo->GetMainParticle();
    Double_t mainX1=hitInfo->GetMainX1();
    Int_t muon=13;

    //test for straight tracks
    if (abs(mainParticle)==muon && 
        (mainX1>mainX1Cut || mainX1<-1*mainX1Cut)){
      straight=false;
      break;
    }
  }

  //return result
  return straight;
}

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

Int_t CDAnalysis::TrueMuonHunter(Int_t event)
{
  if (!fTruthHitInfo) return -1;

  TClonesArray &cTruth = *fTruthHitInfo;
  Int_t numTruthHits=fTruthHitInfo->GetEntries();
  
  Bool_t foundFirstMuon=false;
  Int_t firstMuonPlane=-1;

  //loop over the truth hits
  for (Int_t hit=0;hit<numTruthHits;hit++){
    //cast the tclonesarray up to a truthhitinfo object
    CDTruthHitInfo *hitInfo=
      dynamic_cast<CDTruthHitInfo*>(cTruth[hit]);
    
    Int_t plane=hitInfo->GetPlane();
    Double_t enDep=hitInfo->GetTotalEnDep();
    //Double_t earliestT1=hitInfo->GetEarliestT1();
    Int_t mainParticle=hitInfo->GetMainParticle();
    Double_t mainPartEn=hitInfo->GetMainPartEn();
    Int_t muon=13;

    if (abs(mainParticle)==muon){
      
      if (!foundFirstMuon){
        MSG("CDAnalysis",Msg::kInfo)
          <<"Event="<<event<<": First plane with muon="<<plane<<endl;
        firstMuonPlane=plane;
        foundFirstMuon=true;
      }

      if (plane>32){
        MSG("CDAnalysis",Msg::kInfo)
          <<"Plane="<<plane<<", enDep="<<enDep
          <<", muon energy="<<mainPartEn<<endl;
      }
    }

  }

  return firstMuonPlane;
}

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

Float_t CDAnalysis::TrueXInStripFrame(Int_t plane,Int_t strip,
                                      Int_t event) const
{
  if (!fTruthHitInfo) return -1;

  static Int_t lastEvent=-1;
  static map<Int_t,Double_t> trueX;
  static map<Int_t,Int_t> counter;
  Int_t inStripKey=500*plane+strip;

  //if a new event then have to rebuild the map of x
  //as a function of globalstrip key
  if (event!=lastEvent){
    //set the lastevent
    lastEvent=event;

    //clear the map
    trueX.clear();
    counter.clear();

    //now build the new map for this event

    TClonesArray &cTruth = *fTruthHitInfo;
    Int_t numTruthHits=fTruthHitInfo->GetEntries();
    
    //loop over the truth hits
    for (Int_t hit=0;hit<numTruthHits;hit++){
      //cast the tclonesarray up to a truthhitinfo object
      CDTruthHitInfo *hitInfo=
        dynamic_cast<CDTruthHitInfo*>(cTruth[hit]);
    
      Int_t pmtTruth1=hitInfo->GetPmtTruth1();
      Int_t pmtTruth2=hitInfo->GetPmtTruth2();
      Int_t sumTruth=0;
      
      //add up the bits
      sumTruth|=pmtTruth1;
      sumTruth|=pmtTruth2;

      //check if genuine bit is set
      Bool_t genuineBit=((sumTruth & DigiSignal::kGenuine)==
                         DigiSignal::kGenuine);

      //if not genuine then skip this strip
      if (!genuineBit) {
        continue;
      }

      Int_t plane=hitInfo->GetPlane();
      Int_t strip=hitInfo->GetStrip();
      Int_t stripKey=500*plane+strip;

      Double_t mainX1=hitInfo->GetMainX1();
      Double_t mainX2=hitInfo->GetMainX2();
      Double_t avX=(mainX1+mainX2)*0.5;

      if (avX==0){
        MSG("CDAnalysis",Msg::kDebug)
          <<"CDAnalysis::TrueXInStripFrame::average x="<<avX
          <<", enDep="<<hitInfo->GetTotalEnDep()<<endl;
      }

      trueX[stripKey]=avX;
      counter[stripKey]++;
    }
  }

  if (counter[inStripKey]>0) return trueX[inStripKey];
  else return -1;
}

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

Double_t CDAnalysis::TrueEnDep(Int_t inplane,Int_t event) const
{
  if (!fTruthHitInfo) return 0;

  static Int_t lastEvent=-1;
  static map<Int_t,Double_t> plEnDep;

  //if a new event then have to rebuild the map of energy dep
  //as a function of plane
  if (event!=lastEvent){
    //set the lastevent
    lastEvent=event;

    //clear the map
    plEnDep.clear();

    //now build the new map for this event

    TClonesArray &cTruth = *fTruthHitInfo;
    Int_t numTruthHits=fTruthHitInfo->GetEntries();
    
    //loop over the truth hits
    for (Int_t hit=0;hit<numTruthHits;hit++){
      //cast the tclonesarray up to a truthhitinfo object
      CDTruthHitInfo *hitInfo=
        dynamic_cast<CDTruthHitInfo*>(cTruth[hit]);
    
      Int_t plane=hitInfo->GetPlane();
      Double_t enDep=hitInfo->GetTotalEnDep();
      //Double_t earliestT1=hitInfo->GetEarliestT1();

      //cut out the decay electron
      //if (earliestT1>15e-9) continue;

      plEnDep[plane]+=enDep;
    }
  }

  return plEnDep[inplane];
}

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

Double_t CDAnalysis::TrueMainPL(Int_t inplane,Int_t event) const
{
  if (!fTruthHitInfo) return 0;

  static Int_t lastEvent=-1;
  static map<Int_t,Double_t> plPL;
  //static map<Int_t,Int_t> plStrip;

  //if a new event then have to rebuild the map of energy dep
  //as a function of plane
  if (event!=lastEvent){
    MAXMSG("CDAnalysis",Msg::kDebug,50)
      <<"Building map..."<<endl;
    //set the lastevent
    lastEvent=event;

    //clear the map
    plPL.clear();
    //plStrip.clear();

    //now build the new map for this event

    TClonesArray &cTruth = *fTruthHitInfo;
    Int_t numTruthHits=fTruthHitInfo->GetEntries();
    
    //loop over the truth hits
    for (Int_t hit=0;hit<numTruthHits;hit++){
      //cast the tclonesarray up to a truthhitinfo object
      CDTruthHitInfo *hitInfo=
        dynamic_cast<CDTruthHitInfo*>(cTruth[hit]);
    
      Int_t plane=hitInfo->GetPlane();
      Int_t strip=hitInfo->GetStrip();
      Double_t pathLen=hitInfo->GetMainPathLength();
      Int_t particle=hitInfo->GetMainParticle();
      Int_t muon=13;

      if (abs(particle)!=muon) {
        MAXMSG("CDAnalysis",Msg::kDebug,100)
          <<"Not a muon so ignoring..."<<endl;
        continue;
      }

      MAXMSG("CDAnalysis",Msg::kDebug,500)
        <<"plane="<<plane//<<", lastSt="<<plStrip[plane]
        <<" (currSt="<<strip<<")"<<", lastPL="<<plPL[plane]*100
        <<", currentPL="<<pathLen*100<<endl;
      plPL[plane]+=pathLen;
      //plStrip[plane]=strip;
    }
    MAXMSG("CDAnalysis",Msg::kDebug,50)
      <<"Finished building map"<<endl;
  }

  return plPL[inplane];
}

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

Double_t CDAnalysis::GetMainParticleEnergy() const
{
  if (!fTruthHitInfo) return 0;

  Double_t mainEnergy=0;

  TClonesArray &cTruth = *fTruthHitInfo;
  Int_t numTruthHits=fTruthHitInfo->GetEntries();
  
  //loop over the truth hits
  for (Int_t hit=0;hit<numTruthHits;hit++){
    //cast the tclonesarray up to a truthhitinfo object
    CDTruthHitInfo *hitInfo=
      dynamic_cast<CDTruthHitInfo*>(cTruth[hit]);
    
    Double_t en=hitInfo->GetMainPartEn();

    MSG("CDAnalysis",Msg::kDebug)
      <<"en="<<en<<endl;

    //find the greatest energy
    if (en>mainEnergy) {
      mainEnergy=en;
      MSG("CDAnalysis",Msg::kDebug)
        <<"mainEnergy="<<mainEnergy<<endl;
    }
  }
 
  return mainEnergy;
}

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

Double_t CDAnalysis::GetGreatestMainParticleEnergy(Int_t nEvents) const
{
  Double_t mainEnergy=0;

  if (nEvents>fEvents) nEvents=fEvents;

  for(Int_t event=0;event<nEvents;event++){
    
    //get the snarl/event
    fTrackerTree->GetEvent(event); 
    
    Double_t mEn=this->GetMainParticleEnergy();
    
    //select the greatest one
    if (mEn>mainEnergy) mainEnergy=mEn;
  }

  return mainEnergy;
}

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

Double_t CDAnalysis::CalcMeuPLCorrected(std::map<Int_t,Float_t> plEnDep,
                                        std::map<Int_t,Float_t> plPLCor,
                                        Int_t event,
                                        Float_t* GeVPerMeu
                                        ) const
{
  if (plEnDep.size()==0 || plPLCor.size()==0){
    MAXMSG("CDAnalysis",Msg::kWarning,10000)
      <<"Can't CalcMeuPLCorrected, map has size zero"<<endl;
    return -1;
  }

  Float_t local14_16=0;
  Float_t localSigCor=0;//for calc meu without fractions of planes
  Float_t localEnDep=0;//for calc meu without fractions of planes
  Int_t plCounter=0;
  Float_t initValue=-999;
  Float_t materialFromTrkEnd=0;
  Float_t materialInWindow=0;
  const Float_t material01Pl=(2.5+1)*Munits::cm;
  const Float_t material16Pl=16*material01Pl;
  const Float_t material14Pl=14*material01Pl;
  MAXMSG("CDAnalysis",Msg::kInfo,1)
    <<"Using 14 pls material="<<material14Pl
    <<", 16 plns="<<material16Pl
    <<", 1 plns="<<material01Pl<<endl;

  Int_t endPlane=this->GetEventLength();

  //calculate energy deposition in the window
  if (endPlane>=0){
    MAXMSG("CDAnalysis",Msg::kDebug,100)
      <<"Analysing track window, vtxPl=0"
      <<", endPl="<<endPlane<<endl;
    
    //loop over the track
    Int_t pl=endPlane;
    while (pl!=0){
      
      Float_t pathLengthCor=initValue;
      pathLengthCor=plPLCor[pl];
      if (pathLengthCor<1){
        MAXMSG("CDAnalysis",Msg::kWarning,1000)
          <<"pl="<<pl<<", endPlane="<<endPlane
          <<", path len cor wrong="<<pathLengthCor<<endl;
        //set default to 1
        pathLengthCor=1;
      }

      //sum up energy in window
      if ((materialFromTrkEnd+
           pathLengthCor*material01Pl)>material16Pl && 
          //check if window starts within or before this plane
          //&& that don't already have whole window of material
          materialInWindow<material14Pl){

        //alternatively start window after the plane in which 
        //the window should officially start part way through
        //if (materialFromTrkEnd>material16Pl && 
        //materialInWindow<material14Pl){
        
        //work out meu without using a fraction of a plane
        Float_t plCorEnDep=0;
        if (pathLengthCor) plCorEnDep=plEnDep[pl]/pathLengthCor;
        localSigCor+=plCorEnDep;
        plCounter++;
        //do it for the true en dep as well
        Float_t trueEnDep=this->TrueEnDep(pl,event);
        MAXMSG("CDAnalysis",Msg::kInfo,100)
          <<"trueEnDep="<<trueEnDep
          <<", PLCor="<<pathLengthCor
          <<", sum="<<localEnDep<<endl;
        if (pathLengthCor) trueEnDep/=pathLengthCor;
        localEnDep+=trueEnDep;

        //check if adding this much energy takes you over
        //the required amount of material
        if (materialInWindow+pathLengthCor*material01Pl>material14Pl){

          Float_t materialNeeded=material14Pl-materialInWindow;
          //calc the fraction of material/energy needed from this plane
          Float_t fractNeeded=materialNeeded/
            (pathLengthCor*material01Pl);
          
          MAXMSG("CDAnalysis",Msg::kDebug,100)
            <<"Last Chunk: matNeed="<<materialNeeded
            <<", matPl="<<pathLengthCor*material01Pl
            <<", 14Pl="<<material14Pl
            <<", matWin="<<materialInWindow
            <<", matTrk="<<materialFromTrkEnd<<endl
            <<"fractNeeded="<<fractNeeded<<", en in plane="<<plEnDep[pl]
            <<", en used="<<fractNeeded*plEnDep[pl]<<endl;

          //increment the material in window
          materialInWindow+=fractNeeded*(pathLengthCor*material01Pl);

          //add the required fraction of energy in this plane to
          //the total
          local14_16+=fractNeeded*plEnDep[pl];
        }
        else{
          //add the energy in this plane to the total
          materialInWindow+=pathLengthCor*material01Pl;

          MAXMSG("CDAnalysis",Msg::kDebug,500)
            <<"Window:p="<<pl
            <<", matPl="<<pathLengthCor*material01Pl
            <<", matWin="<<materialInWindow
            <<", matTrk="<<materialFromTrkEnd
            <<", plEn="<<plEnDep[pl]
            <<endl;

          //increment the energy deposited
          local14_16+=plEnDep[pl];
        }
      }

      MAXMSG("CDAnalysis",Msg::kDebug,1000)
        <<"local: matPl="<<pathLengthCor*material01Pl
        <<", matTrk="<<materialFromTrkEnd
        <<", matWin="<<materialInWindow<<", sumEnDep="<<local14_16
        <<"PLCor="<<pathLengthCor
        <<endl;

      //sum up material traversed
      if (pathLengthCor>0) materialFromTrkEnd+=pathLengthCor*
                             material01Pl;
      else {
        MAXMSG("CDAnalysis",Msg::kInfo,1000)
          <<"Ahhh path len cor wrong="<<pathLengthCor<<endl;
      }

      //increment the plane
      pl--;
    }
  }

  Float_t winSigCor=0;
  if (plCounter){
    winSigCor=localSigCor/plCounter;
    localEnDep/=plCounter;
    //set the output variable
    if (GeVPerMeu) *GeVPerMeu=localEnDep;
  }

  //meu from fraction of last plane in window
  local14_16/=14;

  //decide which to use
  Float_t meu=local14_16;
  Bool_t useNoFraction=true;
  if (useNoFraction) meu=winSigCor;

  //calc fraction of material in window actually obtained
  Float_t fractWin=materialInWindow/material14Pl;
  
  MAXMSG("CDAnalysis",Msg::kInfo,500)
    <<"14_16="<<local14_16<<", winSigCor="<<winSigCor
    <<", plCounter="<<plCounter<<", fractWin="<<fractWin<<endl;

  //don't return meu when only a fraction of window was obtained
  if (fractWin<1) meu=-1;

  return meu;
}

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

void CDAnalysis::FillEnVsDist(std::map<Int_t,Float_t> plEnDep,
                              std::map<Int_t,Float_t> plPLCor,
                              Int_t event) const
{
  if (plEnDep.size()==0 || plPLCor.size()==0){
    MAXMSG("CDAnalysis",Msg::kWarning,10000)
      <<"Can't FillEnVsDist, map has size zero"<<endl;
    return;
  }

  Float_t materialFromTrkEnd=0;
  Float_t trueMaterialFromTrkEnd=0;
  Int_t planesFromTrkEnd=0;
  const Float_t material01Pl=(5.94)*Munits::cm;//including air gap
  Float_t initValue=-999;
  
  
  //static Int_t event=-999;
  map<Int_t,Float_t> truePLCor;
  this->TruePLCor(truePLCor,event);
  //event++;

  static TProfile* pEnVsDist=0;
  static TProfile* pEnVsDistAll=0;
  static TProfile* pEnVsDistAll2=0;
  static TProfile* pEnVsDist0=0;
  static TProfile* pTrueEnVsDist0=0;
  static TProfile* pTrueEnPLCVsDist0=0;
  static TProfile* pTrueEnPLCVsDist0270=0;
  static TProfile* pTrueEnPLCVsDist0290=0;
  static TProfile* pTrueEnPLCVsDist0310=0;
  static TProfile* pEnVsPlane0=0;
  static TProfile* pTrueEnVsPlane0=0;
  static TProfile* pEnVsDistNoEnd=0;
  static TProfile* pEnVsDistNoEnd4951=0;
  static TProfile* pEnVsDistNoEnd40=0;
  static TProfile* pEnVsDistNoEnd250=0;
  static TProfile* pEnVsDistNoEnd270=0;
  static TProfile* pEnVsDistNoEnd290=0;
  static TProfile* pEnVsDistNoEnd290a=0;
  static TProfile* pEnVsDistNoEnd290b=0;
  static TProfile* pEnVsDistNoEnd310=0;
  static TProfile* pEnVsDistNoEndTrue=0;
  static TProfile* pEnVsDistNoEndTrue2=0;
  static TProfile* pPLCorVsPlane290=0;
  static TProfile* pPLCorVsPlaneTrue=0;
  static TH1F* hMatFromTrkEnd=0;
  static TH1F* hMatFromTrkEndTrue=0;
  static TH1F* hSigCorPLC=0;

  static vector<TH1F*> meuSigCorNoPLCorHistos(60);
  static vector<TH1F*> meuSigCorPLCorHistos(60);
  static vector<TH1F*> meuSigCorNoPLCorPerpHistos(60);
  static vector<TH1F*> meuSigCorPLCorPerpHistos(60);
  static vector<TH1F*> meuSigCorNoPLCorZerosHistos(60);
  static vector<TH1F*> meuSigCorPLCorZerosHistos(60);

  if (pEnVsDist==0) {
    pEnVsDist=new TProfile("pEnVsDist","pEnVsDist",60,0,
                           (60*material01Pl)/Munits::cm);
    pEnVsDistAll=new TProfile("pEnVsDistAll","pEnVsDistAll",60,0,
                              (60*material01Pl)/Munits::cm);
    pEnVsDistAll2=new TProfile("pEnVsDistAll2","pEnVsDistAll2",60,0,
                               (60*material01Pl)/Munits::cm);
    pEnVsDistNoEnd=new TProfile("pEnVsDistNoEnd","pEnVsDistNoEnd",60,0,
                                (60*material01Pl)/Munits::cm);
    pEnVsDist0=new TProfile("pEnVsDist0","pEnVsDist0",60,0,
                            (60*material01Pl)/Munits::cm);
    pEnVsPlane0=new TProfile("pEnVsPlane0","pEnVsPlane0",60,0,60);
    pTrueEnVsDist0=new TProfile("pTrueEnVsDist0","pTrueEnVsDist0",60,0,
                                (60*material01Pl)/Munits::cm);
    pTrueEnPLCVsDist0=new TProfile("pTrueEnPLCVsDist0",
                                   "pTrueEnPLCVsDist0",
                                   60,0,(60*material01Pl)/Munits::cm);
    pTrueEnPLCVsDist0270=new TProfile("pTrueEnPLCVsDist0270",
                                   "pTrueEnPLCVsDist0270",
                                   60,0,(60*material01Pl)/Munits::cm);
    pTrueEnPLCVsDist0290=new TProfile("pTrueEnPLCVsDist0290",
                                   "pTrueEnPLCVsDist0290",
                                   60,0,(60*material01Pl)/Munits::cm);
    pTrueEnPLCVsDist0310=new TProfile("pTrueEnPLCVsDist0310",
                                   "pTrueEnPLCVsDist0310",
                                   60,0,(60*material01Pl)/Munits::cm);
    pTrueEnVsPlane0=new TProfile("pTrueEnVsPlane0","pTrueEnVsPlane0",
                                 60,0,60);
    pEnVsDistNoEnd4951=new TProfile("pEnVsDistNoEnd4951",
                                    "pEnVsDistNoEnd4951",60,0,
                                    (60*material01Pl)/Munits::cm);
    pEnVsDistNoEnd40=new TProfile("pEnVsDistNoEnd40",
                                  "pEnVsDistNoEnd40",60,0,
                                  (60*material01Pl)/Munits::cm);
    pEnVsDistNoEnd250=new TProfile("pEnVsDistNoEnd250",
                                   "pEnVsDistNoEnd250",60,0,
                                   (60*material01Pl)/Munits::cm);
    pEnVsDistNoEnd270=new TProfile("pEnVsDistNoEnd270",
                                   "pEnVsDistNoEnd270",60,0,
                                   (60*material01Pl)/Munits::cm);
    pEnVsDistNoEnd290=new TProfile("pEnVsDistNoEnd290",
                                   "pEnVsDistNoEnd290",60,0,
                                   (60*material01Pl)/Munits::cm);
    pEnVsDistNoEnd290a=new TProfile("pEnVsDistNoEnd290a",
                                   "pEnVsDistNoEnd290a",60,0,
                                  (60*material01Pl)/Munits::cm);
    pEnVsDistNoEnd290b=new TProfile("pEnVsDistNoEnd290b",
                                   "pEnVsDistNoEnd290b",60,0,
                                  (60*material01Pl)/Munits::cm);
    pEnVsDistNoEnd310=new TProfile("pEnVsDistNoEnd310",
                                   "pEnVsDistNoEnd310",60,0,
                                   (60*material01Pl)/Munits::cm);
    pEnVsDistNoEndTrue=new TProfile("pEnVsDistNoEndTrue",
                                    "pEnVsDistNoEndTrue",60,0,
                                    (60*material01Pl)/Munits::cm);
    pEnVsDistNoEndTrue2=new TProfile("pEnVsDistNoEndTrue2",
                                     "pEnVsDistNoEndTrue2",60,0,
                                     (60*material01Pl)/Munits::cm);
    hMatFromTrkEnd=new TH1F("hMatFromTrkEnd","hMatFromTrkEnd",200,0,4);
    hMatFromTrkEndTrue=new TH1F("hMatFromTrkEndTrue",
                                "hMatFromTrkEndTrue",200,0,4);
    hSigCorPLC=new TH1F("hSigCorPLC","hSigCorPLC",5000,-100,20000);
    pPLCorVsPlane290=new TProfile("pPLCorVsPlane290",
                                  "pPLCorVsPlane290",60,0,60);
    pPLCorVsPlaneTrue=new TProfile("pPLCorVsPlaneTrue",
                                   "pPLCorVsPlaneTrue",60,0,60);

    for (Int_t i=0;i<60;i++){
      string sNameSigCorNoPLCor="hMeuSigCorNoPLCorPlane";
      string sNameSigCorPLCor="hMeuSigCorPLCorPlane";
      string sNameSigCorNoPLCorPerp="hMeuSigCorNoPLCorPerpPlane";
      string sNameSigCorPLCorPerp="hMeuSigCorPLCorPerpPlane";
      string sNameSigCorNoPLCorZeros="hMeuSigCorNoPLCorZerosPlane";
      string sNameSigCorPLCorZeros="hMeuSigCorPLCorZerosPlane";
      //string sNameAdc="hAdcNoPLCorPlane";
      string sNum=Form("%d",i);
      sNameSigCorNoPLCor+=sNum;
      sNameSigCorPLCor+=sNum;
      sNameSigCorNoPLCorPerp+=sNum;
      sNameSigCorPLCorPerp+=sNum;
      sNameSigCorNoPLCorZeros+=sNum;
      sNameSigCorPLCorZeros+=sNum;
      
      meuSigCorNoPLCorHistos[i]=new TH1F
        (sNameSigCorNoPLCor.c_str(),sNameSigCorNoPLCor.c_str(),
         300,-2,3000);
      meuSigCorPLCorHistos[i]=new TH1F
        (sNameSigCorPLCor.c_str(),sNameSigCorPLCor.c_str(),
         300,-2,3000);

      meuSigCorNoPLCorPerpHistos[i]=new TH1F
        (sNameSigCorNoPLCorPerp.c_str(),sNameSigCorNoPLCorPerp.c_str(),
         300,-2,3000);
      meuSigCorPLCorPerpHistos[i]=new TH1F
        (sNameSigCorPLCorPerp.c_str(),sNameSigCorPLCorPerp.c_str(),
         300,-2,3000);
      
      meuSigCorNoPLCorZerosHistos[i]=new TH1F
        (sNameSigCorNoPLCorZeros.c_str(),
         sNameSigCorNoPLCorZeros.c_str(),
         300,-2,3000);
      meuSigCorPLCorZerosHistos[i]=new TH1F
        (sNameSigCorPLCorZeros.c_str(),
         sNameSigCorPLCorZeros.c_str(),
         300,-2,3000);
    }

    MSG("CDAnalysis",Msg::kInfo)
      <<"Creating TProfile, pEnVsDist, max bin="
      <<(60*material01Pl)/Munits::cm<<endl;
  }

  //Int_t endPlane=this->GetEventLength();
  Int_t endPlane=this->CalcLastPlaneOnTrkNoXTalk();
  Int_t planeToStopBefore=0;//have to change length cuts too!!!

  //calculate energy deposition in the window
  if (endPlane>=0){
    MAXMSG("CDAnalysis",Msg::kInfo,100)
      <<"Filling histo, vtxPl=0"
      <<", endPl="<<endPlane<<endl;
    
    //loop over the track
    Int_t pl=endPlane;
    while (pl!=planeToStopBefore){//cut out the first x planes (en dep is low, !d-ray)
      
      Float_t pathLengthCor=initValue;
      pathLengthCor=plPLCor[pl];
      Float_t truePathLengthCor=truePLCor[pl];
      //might be beyond plane end so use reco PLCor
      if (truePathLengthCor<1) truePathLengthCor=pathLengthCor;
      if (pl==endPlane){
        //cap the very last plane
        if (truePathLengthCor>3) {
          MAXMSG("CDAnalysis",Msg::kWarning,1000)
            <<"Capping the last plane PLCor"<<endl;
          truePathLengthCor=3;
        }
      }

      if (pathLengthCor<1){
        MAXMSG("CDAnalysis",Msg::kWarning,1000)
          <<"pl="<<pl<<", endPlane="<<endPlane
          <<", path len cor wrong="<<pathLengthCor<<endl;
        //set default to 1
        pathLengthCor=1;
      }

      Float_t meuSigCor=0;
      if (pathLengthCor) meuSigCor=plEnDep[pl]/pathLengthCor;
      Float_t meuSigCorTrue=0;
      if (truePathLengthCor) meuSigCorTrue=plEnDep[pl]/
                               truePathLengthCor;
      Float_t plMaterial=pathLengthCor*material01Pl;
      Float_t truePlMaterial=truePathLengthCor*material01Pl;

      //plot the value half way through the range of dE/dx it 
      //effectively integrates over
      Float_t toPlotMaterial=plMaterial*0.5+materialFromTrkEnd;
      Float_t trueToPlotMaterial=truePlMaterial*0.5+
        trueMaterialFromTrkEnd;

      Float_t trueEnDep=this->TrueEnDep(pl,event);
      Float_t trueEnDepPLC=-1;
      if (truePathLengthCor) {
        trueEnDepPLC=trueEnDep/truePathLengthCor;
      }
      
      //fill the histograms for each individual plane
      if (pl<60 && pl>=0) {
        //fill both with and without the pathlength correction
        //fill as a function of planes from track end point
        meuSigCorNoPLCorHistos[planesFromTrkEnd]->Fill(plEnDep[pl]);
        meuSigCorPLCorHistos[planesFromTrkEnd]->Fill(meuSigCor);

        //fill as a function of distance from track end point
        //distance in units of "perpendicular planes"
        Int_t planesFromTrkEndPerp=
          static_cast<Int_t>(materialFromTrkEnd/material01Pl);
        MAXMSG("CDAnalysis",Msg::kInfo,1000)
          <<"pl="<<pl<<", endPlane="<<endPlane
          <<", plFrTrkEnd="<<planesFromTrkEnd
          <<", plFrTrkEndPerp="<<planesFromTrkEndPerp
          <<", matFrTrkEnd="<<materialFromTrkEnd
          <<", material01Pl="<<material01Pl<<endl;
        if (planesFromTrkEndPerp<60 && planesFromTrkEndPerp>=0) {
          meuSigCorNoPLCorPerpHistos[planesFromTrkEndPerp]->
            Fill(plEnDep[pl]);
          meuSigCorPLCorPerpHistos[planesFromTrkEndPerp]->
            Fill(meuSigCor);
        }
        else cout<<"ahhh, planesFromTrkEndPerp="
                 <<planesFromTrkEndPerp<<endl;

        //fill as a function of plane number
        meuSigCorNoPLCorZerosHistos[pl]->Fill(plEnDep[pl]);
        meuSigCorPLCorZerosHistos[pl]->Fill(meuSigCor);
      }
      else cout<<"ahhh, pl="<<pl<<endl;

      //fill plot of PLC energy deposition
      hSigCorPLC->Fill(meuSigCor);
      pEnVsDistAll->Fill((materialFromTrkEnd/Munits::cm),meuSigCor);
      if (meuSigCor<7000) pEnVsDistAll2->Fill
                            ((materialFromTrkEnd/Munits::cm),meuSigCor);

      //put a check on the crazy brems
      //can get high landau-tail events in individual planes
      if (meuSigCor<5000){
        pEnVsDist->Fill((toPlotMaterial/Munits::cm),meuSigCor);
        
        //fill using 0 for the first pl, 
        //don't try and get an average position
        pEnVsDist0->Fill((materialFromTrkEnd/Munits::cm),meuSigCor);
        pEnVsPlane0->Fill(pl,meuSigCor);
        pTrueEnVsDist0->Fill((materialFromTrkEnd/Munits::cm),trueEnDep);
        pTrueEnPLCVsDist0->Fill((materialFromTrkEnd/Munits::cm),
                                trueEnDepPLC);
        pTrueEnVsPlane0->Fill(pl,trueEnDep);
        
        if (pl!=endPlane) {
          pEnVsDistNoEnd->Fill((toPlotMaterial/Munits::cm),meuSigCor);
          if (truePLCor.size()>0){
            pEnVsDistNoEndTrue->Fill((trueToPlotMaterial/Munits::cm),
                                     meuSigCorTrue);
            //make one with same positions as reco but true enDep
            pEnVsDistNoEndTrue2->Fill((toPlotMaterial/Munits::cm),
                                      meuSigCorTrue);
          }
        }
      }

      if (meuSigCor<5000 && endPlane>=49 && endPlane<=51){
        if (pl!=endPlane) {
          pEnVsDistNoEnd4951->Fill((toPlotMaterial/Munits::cm),
                                   meuSigCor);
        }
      }

      if (meuSigCor<5000 && endPlane>40){
        if (pl!=endPlane) {
          pEnVsDistNoEnd40->Fill((toPlotMaterial/Munits::cm),
                                 meuSigCor);
        }
      }

      //sum up material traversed
      materialFromTrkEnd+=plMaterial;
      trueMaterialFromTrkEnd+=truePlMaterial;
      planesFromTrkEnd++;

      MAXMSG("CDAnalysis",Msg::kInfo,500)
        <<"Window:p="<<pl
        <<", matPl="<<pathLengthCor*material01Pl
        <<", matTrk="<<materialFromTrkEnd
        <<", plEn="<<plEnDep[pl]
        <<endl;

      //increment the plane
      pl--;
    }
    hMatFromTrkEnd->Fill(materialFromTrkEnd);
    hMatFromTrkEndTrue->Fill(trueMaterialFromTrkEnd);
    Float_t totalMaterial=materialFromTrkEnd;
    materialFromTrkEnd=0;
    trueMaterialFromTrkEnd=0;
    
    //2nd loop over the track
    pl=endPlane;
    while (pl!=planeToStopBefore){//cut out the first x planes (en dep is low, !d-ray)
      
      Float_t pathLengthCor=initValue;
      pathLengthCor=plPLCor[pl];
      Float_t truePathLengthCor=truePLCor[pl];
      //might be beyond plane end so use reco PLCor
      if (truePathLengthCor<1) truePathLengthCor=pathLengthCor;
      if (pl==endPlane){
        //cap the very last plane
        if (truePathLengthCor>3) {
          MAXMSG("CDAnalysis",Msg::kWarning,1000)
            <<"Capping the last plane PLCor"<<endl;
          truePathLengthCor=3;
        }
      }

      if (pathLengthCor<1){
        MAXMSG("CDAnalysis",Msg::kWarning,1000)
          <<"pl="<<pl<<", endPlane="<<endPlane
          <<", path len cor wrong="<<pathLengthCor<<endl;
        //set default to 1
        pathLengthCor=1;
      }

      Float_t meuSigCor=0;
      if (pathLengthCor) meuSigCor=plEnDep[pl]/pathLengthCor;
      Float_t meuSigCorTrue=0;
      if (truePathLengthCor) meuSigCorTrue=plEnDep[pl]/
                               truePathLengthCor;
      Float_t plMaterial=pathLengthCor*material01Pl;
      Float_t truePlMaterial=truePathLengthCor*material01Pl;

      //plot the value half way through the range of dE/dx it 
      //effectively integrates over
      Float_t toPlotMaterial=plMaterial*0.5+materialFromTrkEnd;
      //Float_t trueToPlotMaterial=truePlMaterial*0.5+
      //trueMaterialFromTrkEnd;

      Float_t trueEnDep=this->TrueEnDep(pl,event);
      Float_t trueEnDepPLC=-1;
      if (truePathLengthCor) {
        trueEnDepPLC=trueEnDep/truePathLengthCor;
      }

      //put a check on the crazy brems and cut out the really steep ones
      //can also get high landau-tail events in individual planes
      if (meuSigCor<5000 && totalMaterial>2.49 && totalMaterial<2.61){
        if (pl!=endPlane) {
          pEnVsDistNoEnd250->Fill((toPlotMaterial/Munits::cm),
                                  meuSigCor);
        }
      }
      if (meuSigCor<5000 && totalMaterial>2.7 && totalMaterial<2.9){
        if (pl!=endPlane) {
          pEnVsDistNoEnd270->Fill((toPlotMaterial/Munits::cm),
                                  meuSigCor);
          pTrueEnPLCVsDist0270->Fill((materialFromTrkEnd/Munits::cm),
                                     trueEnDepPLC);
        }
      }
      //was 2.8->3.0 when pl!=1
      if (meuSigCor<5000 && totalMaterial>2.9 && totalMaterial<3.1){
        if (pl!=endPlane) {
          pPLCorVsPlane290->Fill(pl,pathLengthCor);
          pTrueEnPLCVsDist0290->Fill((materialFromTrkEnd/Munits::cm),
                                     trueEnDepPLC);
          pPLCorVsPlaneTrue->Fill(pl,truePathLengthCor);
          pEnVsDistNoEnd290->Fill((toPlotMaterial/Munits::cm),
                                  meuSigCor);
        }
      }
      if (meuSigCor<5000 && totalMaterial>3.1 && totalMaterial<3.3){
        if (pl!=endPlane) {
          pEnVsDistNoEnd310->Fill((toPlotMaterial/Munits::cm),
                                  meuSigCor);
          pTrueEnPLCVsDist0310->Fill((materialFromTrkEnd/Munits::cm),
                                     trueEnDepPLC);
        }
      }

      if (meuSigCor<5000 && totalMaterial>2.95 && totalMaterial<3.05){
        if (pl!=endPlane) {
          pEnVsDistNoEnd290a->Fill((toPlotMaterial/Munits::cm),
                                   meuSigCor);
        }
      }


      if (meuSigCor<5000 && totalMaterial>2.96 && totalMaterial<3.03){
        if (pl!=endPlane) {
          pEnVsDistNoEnd290b->Fill((toPlotMaterial/Munits::cm),
                                  meuSigCor);
        }
      }

      //sum up material traversed
      materialFromTrkEnd+=plMaterial;
      trueMaterialFromTrkEnd+=truePlMaterial;

      //increment the plane
      pl--;
    }
  }
}

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

Double_t CDAnalysis::TrueEnDepNotSc(Int_t inplane,Int_t event) const
{
  if (!fTruthHitInfo) return 0;

  static Int_t lastEvent=-1;
  static map<Int_t,Double_t> plEnDep;
  static map<Int_t,Double_t> plMainPartEn;
  static Double_t totalEnDep=0;

  //if a new event then have to rebuild the map of energy dep
  //as a function of plane
  if (event!=lastEvent){

    MSG("CDAnalysis",Msg::kVerbose) 
      <<"Clearing the maps..."<<endl;

    //set the lastevent
    lastEvent=event;

    //clear the map
    plEnDep.clear();
    plMainPartEn.clear();
    totalEnDep=0;

    //now build the new map for this event

    TClonesArray &cTruth = *fTruthHitInfo;
    Int_t numTruthHits=fTruthHitInfo->GetEntries();
    
    //loop over the truth hits
    for (Int_t hit=0;hit<numTruthHits;hit++){
      //cast the tclonesarray up to a truthhitinfo object
      CDTruthHitInfo *hitInfo=
        dynamic_cast<CDTruthHitInfo*>(cTruth[hit]);
    
      Int_t plane=hitInfo->GetPlane();
      Int_t mainParticle=hitInfo->GetMainParticle();
      Double_t enDep=hitInfo->GetTotalEnDep();
      Double_t mainPartEn=hitInfo->GetMainPartEn();
      //Double_t earliestT1=hitInfo->GetEarliestT1();
      Int_t muon=13;

      //cut out the decay electron
      //if (earliestT1>15e-9) continue;

      plEnDep[plane]+=enDep;
      //select the biggest mainPartEn
      if (mainPartEn>plMainPartEn[plane] && abs(mainParticle)==muon){
        plMainPartEn[plane]=mainPartEn;      
      }
    }
  }

  Double_t plEnDepNotSc=0;
    
  if (!(plMainPartEn[inplane]==0 || plMainPartEn[inplane+1]==0)) {
    plEnDepNotSc=plMainPartEn[inplane]-plEnDep[inplane]-
      plMainPartEn[inplane+1];
    totalEnDep+=plEnDep[inplane]+plEnDepNotSc;
  }

  MSG("CDAnalysis",Msg::kVerbose) 
    <<"pl="<<inplane<<", pEn="<<plMainPartEn[inplane]
    <<", pEn+1="<<plMainPartEn[inplane+1]
    <<", eDep="<<1000*plEnDep[inplane]
    <<", eDepNotSc="<<1000*plEnDepNotSc
    <<", tE="<<totalEnDep
    <<endl;

  return plEnDepNotSc;
}

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

Int_t CDAnalysis::TrueNumDigiScintHits(Int_t inplane,
                                          Int_t event) const
{
  if (!fTruthHitInfo) return 0;

  static Int_t lastEvent=-1;
  static map<Int_t,Int_t> numHits;

  //if a new event then have to rebuild the map of hits
  //as a function of plane
  if (event!=lastEvent){
    //set the lastevent
    lastEvent=event;

    //clear the map
    numHits.clear();

    //now build the new map for this event

    TClonesArray &cTruth = *fTruthHitInfo;
    Int_t numTruthHits=fTruthHitInfo->GetEntries();
    
    //loop over the truth hits
    for (Int_t hit=0;hit<numTruthHits;hit++){
      //cast the tclonesarray up to a truthhitinfo object
      CDTruthHitInfo *hitInfo=
        dynamic_cast<CDTruthHitInfo*>(cTruth[hit]);
    
      Int_t plane=hitInfo->GetPlane();
      Int_t numDsh=hitInfo->GetNumDigiScintHits();
      //Double_t earliestT1=hitInfo->GetEarliestT1();

      //cut out the decay electron
      //if (earliestT1>15e-9) continue;

      numHits[plane]+=numDsh;
    }
  }

  return numHits[inplane];
}

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

void CDAnalysis::TruthEnDepFe()
{
  map<Int_t,Double_t> plEnDep;
  map<Int_t,Double_t> plMainPartEn;
  Double_t totalEnDep=0;
  vector<TH1F*> hEnDepFe;

  for (Int_t pl=0;pl<60;pl++){
    fS="Energy Deposition in Iron, plane=";
    string sPl=Form("%d",pl);
    fS+=sPl;
    hEnDepFe.push_back(new TH1F(fS.c_str(),fS.c_str(),600,-2,200));
  }

  
  this->InitialiseLoopVariables();  
  
  for(Int_t event=0;event<fEvents;event++){
    
    this->SetLoopVariables(event);
    
    //clear the map
    plEnDep.clear();
    plMainPartEn.clear();
    totalEnDep=0;
    
    //now build the new map for this event

    if (!fTruthHitInfo) break;

    TClonesArray &cTruth = *fTruthHitInfo;
    Int_t numTruthHits=fTruthHitInfo->GetEntries();
    
    //loop over the truth hits
    for (Int_t hit=0;hit<numTruthHits;hit++){
      //cast the tclonesarray up to a truthhitinfo object
      CDTruthHitInfo *hitInfo=
        dynamic_cast<CDTruthHitInfo*>(cTruth[hit]);
    
      Int_t plane=hitInfo->GetPlane();
      Int_t mainParticle=hitInfo->GetMainParticle();
      Double_t enDep=hitInfo->GetTotalEnDep();
      Double_t mainPartEn=hitInfo->GetMainPartEn();
      Double_t earliestT1=hitInfo->GetEarliestT1();
      Int_t muon=13;

      //cut out the decay electron
      if (earliestT1>15e-9) continue;
      cout<<"Cutting out the decay electron"<<endl;

      plEnDep[plane]+=enDep;
      //select the biggest mainPartEn
      if (mainPartEn>plMainPartEn[plane] && abs(mainParticle)==muon){
        plMainPartEn[plane]=mainPartEn;      
      }
    }

    
    //vector<Double_t> plEnDepFe(60);

    //loop 60-1 planes
    for (Int_t p=0;p<59;p++){

      Double_t plEnDepFe=plMainPartEn[p]-plEnDep[p]-
        plMainPartEn[p+1];
      totalEnDep+=plEnDep[p]+plEnDepFe;
      
      //fill the histo
      hEnDepFe[p]->Fill(1000*plEnDepFe);
      
      MSG("CDAnalysis",Msg::kVerbose) 
        <<"pl="<<p<<", pEn="<<plMainPartEn[p]
        <<", pEn+1="<<plMainPartEn[p+1]
        <<", eDep="<<1000*plEnDep[p]
        <<", eDepFe="<<1000*plEnDepFe
        <<", tE="<<totalEnDep
        <<endl;
      
    }
  }//end of for

  TCanvas *cEnDep=new TCanvas("cEnDep","EnDep",0,0,1200,800);
  cEnDep->SetFillColor(0);
  cEnDep->Divide(2,3);
  cEnDep->cd(1);
  hEnDepFe[0]->Draw();
  cEnDep->cd(2);
  hEnDepFe[1]->Draw();
  cEnDep->cd(3);
  hEnDepFe[2]->Draw();
  cEnDep->cd(4);
  hEnDepFe[3]->Draw();
  cEnDep->cd(5);
  hEnDepFe[4]->Draw();
  cEnDep->cd(6);
  hEnDepFe[5]->Draw();

}

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

void CDAnalysis::TruthEventLength()
{
  //only have pid for real data so return if not data
  //if (fSimFlag!=SimFlag::kData) abort(true);

  //open the output file for the histograms
  fOutFile=this->OpenFile(fRunNumber,"EvLen");

  TH1F *hEventLength=new TH1F("hEventLength","EventLength hit",
                              77,-2,75);
  hEventLength->GetXaxis()->SetTitle("Event Length (planes)");
  hEventLength->GetXaxis()->CenterTitle();
  hEventLength->GetYaxis()->SetTitle("");
  hEventLength->GetYaxis()->CenterTitle();
  hEventLength->SetFillColor(0);
  hEventLength->SetBit(TH1::kCanRebin);

  TH1F *hEventLength2=new TH1F("hEventLength2","EventLength",77,-2,75);
  hEventLength2->GetXaxis()->SetTitle("Event Length (planes)");
  hEventLength2->GetXaxis()->CenterTitle();
  hEventLength2->GetYaxis()->SetTitle("");
  hEventLength2->GetYaxis()->CenterTitle();
  hEventLength2->SetFillColor(0);
  hEventLength2->SetBit(TH1::kCanRebin);

  TH1F *hEventLength3=new TH1F("hEventLength3","EventLength",77,-2,75);
  hEventLength3->GetXaxis()->SetTitle("Event Length (planes)");
  hEventLength3->GetXaxis()->CenterTitle();
  hEventLength3->GetYaxis()->SetTitle("");
  hEventLength3->GetYaxis()->CenterTitle();
  hEventLength3->SetFillColor(0);
  hEventLength3->SetBit(TH1::kCanRebin);

  TH1F *hEventLength4=new TH1F("hEventLength4","EventLength",77,-2,75);
  hEventLength4->GetXaxis()->SetTitle("Event Length (planes)");
  hEventLength4->GetXaxis()->CenterTitle();
  hEventLength4->GetYaxis()->SetTitle("");
  hEventLength4->GetYaxis()->CenterTitle();
  hEventLength4->SetFillColor(0);
  hEventLength4->SetBit(TH1::kCanRebin);

  
  this->InitialiseLoopVariables();  
  
  for(Int_t event=0;event<fEvents;event++){
    
    this->SetLoopVariables(event);
    
    Int_t overallLastPlane=0;
    Int_t notNoiseLastPlane=0;
    Int_t particleLastPlane=0;
    //Int_t inTimeLastPlane=0;
    Int_t muonLastPlane=0;
    
    if (!fTruthHitInfo) break;
    TClonesArray &cTruth = *fTruthHitInfo;
    Int_t numTruthHits=fTruthHitInfo->GetEntries();

    //loop over the truth hits
    for (Int_t hit=0;hit<numTruthHits;hit++){
      //cast the tclonesarray up to a truthhitinfo object
      CDTruthHitInfo *hitInfo=
        dynamic_cast<CDTruthHitInfo*>(cTruth[hit]);
    
      Int_t plane=hitInfo->GetPlane();
      Int_t mainParticle=hitInfo->GetMainParticle();
      Int_t pmtTruth1=hitInfo->GetPmtTruth1();
      Int_t pmtTruth2=hitInfo->GetPmtTruth2();
      Double_t earliestT1=hitInfo->GetEarliestT1();
      Int_t muon=13;

      if (plane>overallLastPlane) overallLastPlane=plane;

      Int_t sumTruth=pmtTruth1 | pmtTruth2;
      Bool_t genuine=(sumTruth & DigiSignal::kGenuine)==
        DigiSignal::kGenuine;
      Bool_t fibreLight=(sumTruth & DigiSignal::kFibreLight)==
        DigiSignal::kFibreLight;
      Bool_t darkNoise=(sumTruth & DigiSignal::kDarkNoise)==
        DigiSignal::kDarkNoise;

      if (!genuine && (fibreLight || darkNoise)) continue;
      
      if (plane>notNoiseLastPlane) notNoiseLastPlane=plane;

      if (!genuine) continue;

      //cut out the late decay electron, ND is alive for ~400 ns
      if (earliestT1>500e-9) continue;
      
      if (plane>particleLastPlane) particleLastPlane=plane;

      //in time but not necessarily a muon
      //if (plane>inTimeLastPlane) inTimeLastPlane=plane;

      //only look at muons
      if (abs(mainParticle)!=muon) continue;      

      if (plane>muonLastPlane) muonLastPlane=plane;
    }

    hEventLength->Fill(overallLastPlane+1);
    hEventLength2->Fill(notNoiseLastPlane+1);
    hEventLength3->Fill(particleLastPlane+1);
    //hEventLength3->Fill(inTimeLastPlane+1);
    hEventLength4->Fill(muonLastPlane+1);

  }//end of for

  TCanvas *cGeom=new TCanvas("cGeom","Geom",0,0,1200,800);
  cGeom->SetFillColor(0);
  cGeom->Divide(2,3);
  cGeom->cd(1);
  hEventLength->Draw();
  cGeom->cd(2);
  hEventLength2->Draw();
  hEventLength2->SetLineColor(2);
  cGeom->cd(3);
  hEventLength3->Draw();
  hEventLength3->SetLineColor(3);
  cGeom->cd(4);
  hEventLength4->Draw();
  hEventLength->SetLineColor(4);
  cGeom->cd(5);
  hEventLength4->Draw();
  hEventLength2->Draw("sames");
  cGeom->cd(6);
  hEventLength4->Draw();
  hEventLength3->Draw("sames");
  hEventLength2->Draw("sames");
  hEventLength->Draw("sames");
}

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

void CDAnalysis::MuonResponse()
{
  MSG("CDAnalysis",Msg::kInfo) 
    <<" ** Running MuonResponse method... **"<<endl;

  //open the output file for the histograms
  fOutFile=this->OpenFile(fRunNumber,"MuonRes");

  TH2F *hPeVsRange=new TH2F("hPeVsRange","PeVsRange",
                            20,40,60,50,0,10);
  hPeVsRange->GetXaxis()->SetTitle("Last Plane");
  hPeVsRange->GetXaxis()->CenterTitle();
  hPeVsRange->GetYaxis()->SetTitle("Greatest PE Value");
  hPeVsRange->GetYaxis()->CenterTitle();
  hPeVsRange->SetFillColor(0);
  //hPeVsRange->SetBit(TH1::kCanRebin);

  TH2F *hPeVsRangeGen=new TH2F("hPeVsRangeGen","PeVsRangeGen",
                            20,40,60,50,0,10);
  hPeVsRangeGen->GetXaxis()->SetTitle("Last Plane");
  hPeVsRangeGen->GetXaxis()->CenterTitle();
  hPeVsRangeGen->GetYaxis()->SetTitle("Greatest PE Value");
  hPeVsRangeGen->GetYaxis()->CenterTitle();
  hPeVsRangeGen->SetFillColor(0);
  //hPeVsRangeGen->SetBit(TH1::kCanRebin);

  TH2F *hPeVsRangeGen2=new TH2F("hPeVsRangeGen2","PeVsRangeGen2",
                            20,40,60,50,0,10);
  hPeVsRangeGen2->GetXaxis()->SetTitle("Last Plane");
  hPeVsRangeGen2->GetXaxis()->CenterTitle();
  hPeVsRangeGen2->GetYaxis()->SetTitle("Greatest PE Value");
  hPeVsRangeGen2->GetYaxis()->CenterTitle();
  hPeVsRangeGen2->SetFillColor(0);
  //hPeVsRangeGen2->SetBit(TH1::kCanRebin);

  TH2F *hPeVsRangeGen3=new TH2F("hPeVsRangeGen3","PeVsRangeGen3",
                            20,40,60,50,0,10);
  hPeVsRangeGen3->GetXaxis()->SetTitle("Last Plane");
  hPeVsRangeGen3->GetXaxis()->CenterTitle();
  hPeVsRangeGen3->GetYaxis()->SetTitle("Greatest PE Value");
  hPeVsRangeGen3->GetYaxis()->CenterTitle();
  hPeVsRangeGen3->SetFillColor(0);
  //hPeVsRangeGen3->SetBit(TH1::kCanRebin);

  TH2F *hPeVsRangeXTalk=new TH2F("hPeVsRangeXTalk","PeVsRangeXTalk",
                            20,40,60,50,0,10);
  hPeVsRangeXTalk->GetXaxis()->SetTitle("Last Plane");
  hPeVsRangeXTalk->GetXaxis()->CenterTitle();
  hPeVsRangeXTalk->GetYaxis()->SetTitle("Greatest PE Value");
  hPeVsRangeXTalk->GetYaxis()->CenterTitle();
  hPeVsRangeXTalk->SetFillColor(0);
  //hPeVsRangeXTalk->SetBit(TH1::kCanRebin);

  TH2F *hPeVsRangeXTalk2=new TH2F("hPeVsRangeXTalk2","PeVsRangeXTalk2",
                            20,40,60,50,0,10);
  hPeVsRangeXTalk2->GetXaxis()->SetTitle("Last Plane");
  hPeVsRangeXTalk2->GetXaxis()->CenterTitle();
  hPeVsRangeXTalk2->GetYaxis()->SetTitle("Greatest PE Value");
  hPeVsRangeXTalk2->GetYaxis()->CenterTitle();
  hPeVsRangeXTalk2->SetFillColor(0);
  //hPeVsRangeXTalk2->SetBit(TH1::kCanRebin);

  TH2F *hPeVsRangeXTalk3=new TH2F("hPeVsRangeXTalk3","PeVsRangeXTalk3",
                            20,40,60,50,0,10);
  hPeVsRangeXTalk3->GetXaxis()->SetTitle("Last Plane");
  hPeVsRangeXTalk3->GetXaxis()->CenterTitle();
  hPeVsRangeXTalk3->GetYaxis()->SetTitle("Greatest PE Value");
  hPeVsRangeXTalk3->GetYaxis()->CenterTitle();
  hPeVsRangeXTalk3->SetFillColor(0);
  //hPeVsRangeXTalk3->SetBit(TH1::kCanRebin);

  TH1F *hNumPlanes=new TH1F("hNumPlanes","NumPlanes hit",77,-2,75);
  hNumPlanes->GetXaxis()->SetTitle("Number of Planes Hit");
  hNumPlanes->GetXaxis()->CenterTitle();
  hNumPlanes->GetYaxis()->SetTitle("");
  hNumPlanes->GetYaxis()->CenterTitle();
  hNumPlanes->SetFillColor(0);
  hNumPlanes->SetBit(TH1::kCanRebin);

  TH1F *h15_18All=new TH1F("h15_18All","h15_18All",1000,-2,75);
  h15_18All->GetXaxis()->SetTitle("SigCors");
  h15_18All->GetXaxis()->CenterTitle();
  h15_18All->GetYaxis()->SetTitle("");
  h15_18All->GetYaxis()->CenterTitle();
  h15_18All->SetFillColor(0);
  h15_18All->SetBit(TH1::kCanRebin);

  TH1F *h15_18=new TH1F("h15_18","h15_18",1000,-2,75);
  h15_18->GetXaxis()->SetTitle("SigCors");
  h15_18->GetXaxis()->CenterTitle();
  h15_18->GetYaxis()->SetTitle("");
  h15_18->GetYaxis()->CenterTitle();
  h15_18->SetFillColor(0);
  h15_18->SetBit(TH1::kCanRebin);

  TH1F *hSigCorTotal=new TH1F("hSigCorTotal","hSigCorTotal",1000,-2,75);
  hSigCorTotal->GetXaxis()->SetTitle("SigCors");
  hSigCorTotal->GetXaxis()->CenterTitle();
  hSigCorTotal->GetYaxis()->SetTitle("");
  hSigCorTotal->GetYaxis()->CenterTitle();
  hSigCorTotal->SetFillColor(0);
  hSigCorTotal->SetBit(TH1::kCanRebin);

  TH1F *hSigCorTotalTight=new TH1F("hSigCorTotalTight",
                                   "hSigCorTotalTight",1000,-2,75);
  hSigCorTotalTight->GetXaxis()->SetTitle("SigCors");
  hSigCorTotalTight->GetXaxis()->CenterTitle();
  hSigCorTotalTight->GetYaxis()->SetTitle("");
  hSigCorTotalTight->GetYaxis()->CenterTitle();
  hSigCorTotalTight->SetFillColor(0);
  hSigCorTotalTight->SetBit(TH1::kCanRebin);

  TH1F *hSigCorTotalPeak=new TH1F("hSigCorTotalPeak",
                                  "hSigCorTotalPeak",1000,-2,75);
  hSigCorTotalPeak->GetXaxis()->SetTitle("SigCors");
  hSigCorTotalPeak->GetXaxis()->CenterTitle();
  hSigCorTotalPeak->GetYaxis()->SetTitle("");
  hSigCorTotalPeak->GetYaxis()->CenterTitle();
  hSigCorTotalPeak->SetFillColor(0);
  hSigCorTotalPeak->SetBit(TH1::kCanRebin);

  TH1F *hNumPlanesAll=new TH1F("hNumPlanesAll","NumPlanesAll hit",
                               77,-2,75);
  hNumPlanesAll->GetXaxis()->SetTitle("Number of Planes Hit");
  hNumPlanesAll->GetXaxis()->CenterTitle();
  hNumPlanesAll->GetYaxis()->SetTitle("");
  hNumPlanesAll->GetYaxis()->CenterTitle();
  hNumPlanesAll->SetFillColor(0);
  hNumPlanesAll->SetBit(TH1::kCanRebin);

  TH1F *hNumPlanesMu=new TH1F("hNumPlanesMu","hNumPlanesMu hit",
                              77,-2,75);
  hNumPlanesMu->GetXaxis()->SetTitle("Number of Planes Hit");
  hNumPlanesMu->GetXaxis()->CenterTitle();
  hNumPlanesMu->GetYaxis()->SetTitle("");
  hNumPlanesMu->GetYaxis()->CenterTitle();
  hNumPlanesMu->SetFillColor(0);
  hNumPlanesMu->SetLineColor(2);
  hNumPlanesMu->SetBit(TH1::kCanRebin);
  
  TH1F *hFirstPlane=new TH1F("hFirstPlane","FirstPlane hit",77,-2,75);
  hFirstPlane->GetXaxis()->SetTitle("First Plane Hit");
  hFirstPlane->GetXaxis()->CenterTitle();
  hFirstPlane->GetYaxis()->SetTitle("");
  hFirstPlane->GetYaxis()->CenterTitle();
  hFirstPlane->SetFillColor(0);
  hFirstPlane->SetBit(TH1::kCanRebin);

  TH1F *hEventLength=new TH1F("hEventLength","EventLength hit",
                              77,-2,75);
  hEventLength->GetXaxis()->SetTitle("Event Length (planes)");
  hEventLength->GetXaxis()->CenterTitle();
  hEventLength->GetYaxis()->SetTitle("");
  hEventLength->GetYaxis()->CenterTitle();
  hEventLength->SetFillColor(0);
  hEventLength->SetBit(TH1::kCanRebin);

  TH1F *hEventLength2=new TH1F("hEventLength2","EventLength",77,-2,75);
  hEventLength2->GetXaxis()->SetTitle("Event Length (planes)");
  hEventLength2->GetXaxis()->CenterTitle();
  hEventLength2->GetYaxis()->SetTitle("");
  hEventLength2->GetYaxis()->CenterTitle();
  hEventLength2->SetFillColor(0);
  hEventLength2->SetBit(TH1::kCanRebin);

  TH1F *hEventLength3=new TH1F("hEventLength3","EventLength",77,-2,75);
  hEventLength3->GetXaxis()->SetTitle("Event Length (planes)");
  hEventLength3->GetXaxis()->CenterTitle();
  hEventLength3->GetYaxis()->SetTitle("");
  hEventLength3->GetYaxis()->CenterTitle();
  hEventLength3->SetFillColor(0);
  hEventLength3->SetBit(TH1::kCanRebin);

  TH2F *hCkvAdcVsRange=new TH2F("hCkvAdcVsRange","CkvAdcVsRange",
                                20,40,60,100,-1,1500);
  hCkvAdcVsRange->GetXaxis()->SetTitle("Last Plane");
  hCkvAdcVsRange->GetXaxis()->CenterTitle();
  hCkvAdcVsRange->GetYaxis()->SetTitle("Ckv ADC");
  hCkvAdcVsRange->GetYaxis()->CenterTitle();
  hCkvAdcVsRange->SetFillColor(0);
  hCkvAdcVsRange->SetBit(TH1::kCanRebin);

  //sum of all strip ends, normal and offset planes
  //tracked
  TProfile* pSigCorVsPlane=new TProfile("pSigCorVsPlane",
                                        "pSigCorVsPlane",
                                        62,0,62);
  TProfile* pSigCorVsOsPlane=new TProfile("pSigCorVsOsPlane",
                                          "pSigCorVsOsPlane",
                                          62,0,62);

  //xtalk + untracked hits
  TProfile* pSigCorVsPlaneX=new TProfile("pSigCorVsPlaneX",
                                         "pSigCorVsPlaneX",
                                         62,0,62);
  TProfile* pSigCorVsOsPlaneX=new TProfile("pSigCorVsOsPlaneX",
                                           "pSigCorVsOsPlaneX",
                                           62,0,62);

  //total
  TProfile* pSigCorVsPlaneT=new TProfile("pSigCorVsPlaneT",
                                         "pSigCorVsPlaneT",
                                         62,0,62);
  TProfile* pSigCorVsOsPlaneT=new TProfile("pSigCorVsOsPlaneT",
                                           "pSigCorVsOsPlaneT",
                                           62,0,62);

  //individual strip ends
  //tracked
  TProfile* pSigCorVsPlaneO1=new TProfile("pSigCorVsPlaneO1",
                                          "pSigCorVsPlaneO1",
                                          62,0,62);
  TProfile* pSigCorVsPlaneO2=new TProfile("pSigCorVsPlaneO2",
                                          "pSigCorVsPlaneO2",
                                          62,0,62);
  TProfile* pSigCorVsPlaneE1=new TProfile("pSigCorVsPlaneE1",
                                          "pSigCorVsPlaneE1",
                                          62,0,62);
  TProfile* pSigCorVsPlaneE2=new TProfile("pSigCorVsPlaneE2",
                                          "pSigCorVsPlaneE2",
                                          62,0,62);
  //xtalk + untracked hits
  TProfile* pSigCorVsPlaneO1X=new TProfile("pSigCorVsPlaneO1X",
                                          "pSigCorVsPlaneO1X",
                                          62,0,62);
  TProfile* pSigCorVsPlaneO2X=new TProfile("pSigCorVsPlaneO2X",
                                          "pSigCorVsPlaneO2X",
                                          62,0,62);
  TProfile* pSigCorVsPlaneE1X=new TProfile("pSigCorVsPlaneE1X",
                                          "pSigCorVsPlaneE1X",
                                          62,0,62);
  TProfile* pSigCorVsPlaneE2X=new TProfile("pSigCorVsPlaneE2X",
                                          "pSigCorVsPlaneE2X",
                                          62,0,62);
  //total
  TProfile* pSigCorVsPlaneO1T=new TProfile("pSigCorVsPlaneO1T",
                                          "pSigCorVsPlaneO1T",
                                          62,0,62);
  TProfile* pSigCorVsPlaneO2T=new TProfile("pSigCorVsPlaneO2T",
                                          "pSigCorVsPlaneO2T",
                                          62,0,62);
  TProfile* pSigCorVsPlaneE1T=new TProfile("pSigCorVsPlaneE1T",
                                          "pSigCorVsPlaneE1T",
                                          62,0,62);
  TProfile* pSigCorVsPlaneE2T=new TProfile("pSigCorVsPlaneE2T",
                                          "pSigCorVsPlaneE2T",
                                          62,0,62);

  //individual strip ends - offset planes
  //tracked
  TProfile* pSigCorVsOsPlaneO1=new TProfile("pSigCorVsOsPlaneO1",
                                            "pSigCorVsOsPlaneO1",
                                            62,0,62);
  TProfile* pSigCorVsOsPlaneO2=new TProfile("pSigCorVsOsPlaneO2",
                                            "pSigCorVsOsPlaneO2",
                                            62,0,62);
  TProfile* pSigCorVsOsPlaneE1=new TProfile("pSigCorVsOsPlaneE1",
                                            "pSigCorVsOsPlaneE1",
                                            62,0,62);
  TProfile* pSigCorVsOsPlaneE2=new TProfile("pSigCorVsOsPlaneE2",
                                            "pSigCorVsOsPlaneE2",
                                            62,0,62);

  //xtalk + untracked hits
  TProfile* pSigCorVsOsPlaneO1X=new TProfile("pSigCorVsOsPlaneO1X",
                                             "pSigCorVsOsPlaneO1X",
                                             62,0,62);
  TProfile* pSigCorVsOsPlaneO2X=new TProfile("pSigCorVsOsPlaneO2X",
                                            "pSigCorVsOsPlaneO2X",
                                            62,0,62);
  TProfile* pSigCorVsOsPlaneE1X=new TProfile("pSigCorVsOsPlaneE1X",
                                            "pSigCorVsOsPlaneE1X",
                                            62,0,62);
  TProfile* pSigCorVsOsPlaneE2X=new TProfile("pSigCorVsOsPlaneE2X",
                                            "pSigCorVsOsPlaneE2X",
                                            62,0,62);
  //total
  TProfile* pSigCorVsOsPlaneO1T=new TProfile("pSigCorVsOsPlaneO1T",
                                            "pSigCorVsOsPlaneO1T",
                                            62,0,62);
  TProfile* pSigCorVsOsPlaneO2T=new TProfile("pSigCorVsOsPlaneO2T",
                                            "pSigCorVsOsPlaneO2T",
                                            62,0,62);
  TProfile* pSigCorVsOsPlaneE1T=new TProfile("pSigCorVsOsPlaneE1T",
                                            "pSigCorVsOsPlaneE1T",
                                            62,0,62);
  TProfile* pSigCorVsOsPlaneE2T=new TProfile("pSigCorVsOsPlaneE2T",
                                            "pSigCorVsOsPlaneE2T",
                                            62,0,62);

  TProfile* pEnVsRange=new TProfile("pEnVsRange","pEnVsRange",
                                    62,0,62);

  vector<TProfile*> pSigCorVsDist;
  vector<TProfile*> pSigCorVsDistN;
  vector<Int_t> vWindowSize;
  vector<Float_t> bigAv;
  vector<Int_t> bigAvCounter;
  Float_t sum15_18=0;
  Float_t sum15_18Counter=0;
  Float_t sum15_18Mip=0;
  Float_t sum15_18MipCounter=0;
  Float_t sumSigCorTotal=0;
  Float_t sumSigCorTotalCounter=0;
  Float_t sumMipTotal=0;
  Float_t sumMipTotalCounter=0;

  Int_t lastPlaneOddCounter=0;
  Int_t lastPlaneEvenCounter=0;
  Int_t lastPlaneO1Counter=0;
  Int_t lastPlaneO2Counter=0;
  Int_t lastPlaneE1Counter=0;
  Int_t lastPlaneE2Counter=0;
  Int_t lastPlaneOBothCounter=0;
  Int_t lastPlaneEBothCounter=0;

  Int_t totalCounter=0;
  Int_t firstPlaneInPairCounter=0;
  Int_t sharedPmtCounter=0;
  Int_t sharedPmtXTalkCounter=0;
  Int_t side1Counter=0;
  Int_t side1XTalkCounter=0;
  Int_t side2Counter=0;
  Int_t side2XTalkCounter=0;
  Int_t doubleCounter=0;
  Int_t doubleXTalkCounter=0;
  Int_t notGenuineOrXTalkCounter=0;

  //num maps for all stripends
  map<Int_t,Float_t> num;
  map<Int_t,Float_t> numX;
  map<Int_t,Float_t> numT;
  map<Int_t,Float_t> numOs;
  map<Int_t,Float_t> numOsX;
  map<Int_t,Float_t> numOsT;

  //non-offset num maps for individual stripends
  map<Int_t,Float_t> numO1;
  map<Int_t,Float_t> numO1X;
  map<Int_t,Float_t> numO1T;
  map<Int_t,Float_t> numO2;
  map<Int_t,Float_t> numO2X;
  map<Int_t,Float_t> numO2T;
  map<Int_t,Float_t> numE1;
  map<Int_t,Float_t> numE1X;
  map<Int_t,Float_t> numE1T;
  map<Int_t,Float_t> numE2;
  map<Int_t,Float_t> numE2X;
  map<Int_t,Float_t> numE2T;

  //offset num maps for individual stripends
  map<Int_t,Float_t> numOsO1;
  map<Int_t,Float_t> numOsO1X;
  map<Int_t,Float_t> numOsO1T;
  map<Int_t,Float_t> numOsO2;
  map<Int_t,Float_t> numOsO2X;
  map<Int_t,Float_t> numOsO2T;
  map<Int_t,Float_t> numOsE1;
  map<Int_t,Float_t> numOsE1X;
  map<Int_t,Float_t> numOsE1T;
  map<Int_t,Float_t> numOsE2;
  map<Int_t,Float_t> numOsE2X;
  map<Int_t,Float_t> numOsE2T;

  Int_t noTrackCounter=0;
  Int_t passPidCounter=0;
 
  
  this->InitialiseLoopVariables();  
  
  for(Int_t event=0;event<fEvents;event++){
    
    this->SetLoopVariables(event);

    if (this->CutOnDeadChips()) continue;
    if (this->CutOnPid()) continue;
    //cut on not straight tracks
    //Double_t mainX1Cut=0.5;
    //if (!this->IsStraightTrack(mainX1Cut)){
    //continue;
    //}

    //count the number of events passing the pid
    passPidCounter++;

    map<Int_t,Int_t> numPlanesHit;
    map<Int_t,Int_t> numPlanesHit1;
    map<Int_t,Int_t> numPlanesHit2;
    map<Int_t,Int_t> numPlanesHitAll;
    map<Int_t,Float_t> planeSigCor;
    map<Int_t,Float_t> planeSigCorX;
    map<Int_t,Float_t> planeSigCorT;
    map<Int_t,Float_t> planeSigCorO1;    
    map<Int_t,Float_t> planeSigCorO1X;
    map<Int_t,Float_t> planeSigCorO1T;
    map<Int_t,Float_t> planeSigCorO2;
    map<Int_t,Float_t> planeSigCorO2X;
    map<Int_t,Float_t> planeSigCorO2T;
    map<Int_t,Float_t> planeSigCorE1;
    map<Int_t,Float_t> planeSigCorE1X;
    map<Int_t,Float_t> planeSigCorE1T;
    map<Int_t,Float_t> planeSigCorE2;
    map<Int_t,Float_t> planeSigCorE2X;
    map<Int_t,Float_t> planeSigCorE2T;

    map<Int_t,Float_t> planeMip;
    map<Int_t,Float_t> planeMipX;
    map<Int_t,Float_t> planeMipT;

    map<Int_t,Float_t> planePe;
    map<Int_t,Float_t> planeGreatestPe;

    Double_t sigCorTotal=0;

    //get tclones arrays for this snarl
    TClonesArray &cTrk=*fTrkHitInfo;
    Int_t numTrkHits=fTrkHitInfo->GetEntries();
    TClonesArray &cUnTrk = *fUnTrkHitInfo;
    Int_t numUnTrkHits = fUnTrkHitInfo->GetEntries();
    //CDTrackedHitInfo *trackedHitInfo=(CDTrackedHitInfo*) c[hit];
    TClonesArray &cXTalk = *fXTalkHits;
    Int_t numXTalkHits=fXTalkHits->GetEntries();
    //TClonesArray &cTruth = *fTruthHitInfo;
    Int_t numTruthHits=0;
    if (fTruthHitInfo) numTruthHits=fTruthHitInfo->GetEntries();

    //look for untracked big events
    if (numTrkHits<2 && numXTalkHits+numUnTrkHits>1){
      MSG("CDAnalysis",Msg::kDebug)
        <<"Event="<<event
        <<", Num hits: trk="<<numTrkHits<<", unTrk="<<numUnTrkHits
        <<", xtalk="<<numXTalkHits<<", truth="<<numTruthHits
        <<endl;
      noTrackCounter++;
      //continue;
    }

    //loop over the untracked hits in the snarl
    for (Int_t hit=0;hit<numUnTrkHits;hit++){
      //cast the tclonesarray up to a trackedhitinfo object
      CDTrackedHitInfo *hitInfo=
        dynamic_cast<CDTrackedHitInfo*>(cUnTrk[hit]);

      this->ReadInHitInfo(hitInfo);

      if (fPlane==0 && fChargeSigCor>1000){
        MSG("CDAnalysis",Msg::kDebug)
          <<"XTalk Event="<<event
          <<", fPlane="<<fPlane
          <<", fStrip="<<fStrip
          <<", fStripend="<<fStripend
          <<", fChargeAdc="<<fChargeAdc
          <<", fSigCor="<<fChargeSigCor<<endl;
      }

      //cut out the bad channels
      if (this->CutOnBadPedestals(fPlane,fStrip,fStripend)) continue;

      //count the number of planes hit
      numPlanesHitAll[fPlane]++;

      //check the channel's sanity
      this->StandardSanityChecks();

      //add up the charge of the hits in each plane
      planeSigCorX[fPlane]+=fChargeSigCor;
      planeSigCorT[fPlane]+=fChargeSigCor;
      planeMipX[fPlane]+=fChargeMip;
      planeMipT[fPlane]+=fChargeMip;
      if (fO1) planeSigCorO1X[fPlane]+=fChargeSigCor;
      if (fO2) planeSigCorO2X[fPlane]+=fChargeSigCor;
      if (fE1) planeSigCorE1X[fPlane]+=fChargeSigCor;
      if (fE2) planeSigCorE2X[fPlane]+=fChargeSigCor;
      if (fO1) planeSigCorO1T[fPlane]+=fChargeSigCor;
      if (fO2) planeSigCorO2T[fPlane]+=fChargeSigCor;
      if (fE1) planeSigCorE1T[fPlane]+=fChargeSigCor;
      if (fE2) planeSigCorE2T[fPlane]+=fChargeSigCor;

      if (fChargePe>planeGreatestPe[fPlane]) planeGreatestPe[fPlane]=
                                               fChargePe;

      if (fPlane==0) continue;
      //sum the sigcor in the event
      //make up for no plane 0 with 2 times plane 1
      else if (fPlane==1) sigCorTotal+=2*fChargeSigCor;
      else sigCorTotal+=fChargeSigCor;
    }

    //loop over xtalk hits in the snarl
    for (Int_t hit=0;hit<numXTalkHits;hit++){
      CDXTalkHitInfo *hitInfo=
        dynamic_cast<CDXTalkHitInfo*>(cXTalk[hit]);

      this->ReadInHitInfo(hitInfo);
      
      //cut out the bad channels
      if (this->CutOnBadPedestals(fPlane,fStrip,fStripend)) continue;

      //count the number of planes hit
      numPlanesHitAll[fPlane]++;

      //check the channel's sanity
      this->StandardSanityChecks();

      //add up the charge of the hits in each plane
      planeSigCorX[fPlane]+=fChargeSigCor;
      planeSigCorT[fPlane]+=fChargeSigCor;
      planeMipX[fPlane]+=fChargeMip;
      planeMipT[fPlane]+=fChargeMip;
      if (fO1) planeSigCorO1X[fPlane]+=fChargeSigCor;
      if (fO2) planeSigCorO2X[fPlane]+=fChargeSigCor;
      if (fE1) planeSigCorE1X[fPlane]+=fChargeSigCor;
      if (fE2) planeSigCorE2X[fPlane]+=fChargeSigCor;
      if (fO1) planeSigCorO1T[fPlane]+=fChargeSigCor;
      if (fO2) planeSigCorO2T[fPlane]+=fChargeSigCor;
      if (fE1) planeSigCorE1T[fPlane]+=fChargeSigCor;
      if (fE2) planeSigCorE2T[fPlane]+=fChargeSigCor;

      if (fChargePe>planeGreatestPe[fPlane]) planeGreatestPe[fPlane]=
                                               fChargePe;

      if (fPlane==0) continue;
      //sum the sigcor in the event
      //make up for no plane 0 with 2 times plane 1
      else if (fPlane==1) sigCorTotal+=2*fChargeSigCor;
      else sigCorTotal+=fChargeSigCor;
    }

    Bool_t lastPlaneO1Hit=false;
    Bool_t lastPlaneO2Hit=false;
    Bool_t lastPlaneE1Hit=false;
    Bool_t lastPlaneE2Hit=false;
    Int_t previousfLastPlane=-1;

    //loop over the tracked hits
    for (Int_t hit=0;hit<numTrkHits;hit++){
      CDTrackedHitInfo *trackedHitInfo=
        dynamic_cast<CDTrackedHitInfo*>(cTrk[hit]);
          
      this->ReadInHitInfo(trackedHitInfo);

      //cut out the bad channels
      if (this->CutOnBadPedestals(fPlane,fStrip,fStripend)) continue;

      //check the channel's sanity
      this->StandardSanityChecks();

      //calculate the first and last planes
      this->CalcFirstLastPlane(fPlane);

      //count the number of planes hit
      numPlanesHit[fPlane]++;
      numPlanesHitAll[fPlane]++;

      if (fLastPlane==fPlane && fPlane%2==0){

        //check if a new last plane has been found
        if (fLastPlane!=previousfLastPlane){
          lastPlaneE1Hit=false;
          lastPlaneE2Hit=false;
        }
        previousfLastPlane=fLastPlane;

        //check which stripend to set to true
        if (fStripend==1) lastPlaneE1Hit=true;
        if (fStripend==2) lastPlaneE2Hit=true;

        //set the hits in the other plane to be false
        lastPlaneO1Hit=false;
        lastPlaneO2Hit=false;
        MSG("CDAnalysis",Msg::kVerbose)
          <<"ev="<<event
          <<", last="<<fLastPlane
          <<", pl="<<fPlane
          <<", se="<<fStripend
          <<", E1="<<lastPlaneE1Hit
          <<", E2="<<lastPlaneE2Hit
          <<", O1="<<lastPlaneO1Hit
          <<", O2="<<lastPlaneO2Hit<<endl;
      }
      else if (fLastPlane==fPlane && fPlane%2==1){

        //check if a new last plane has been found
        if (fLastPlane!=previousfLastPlane){
          lastPlaneO1Hit=false;
          lastPlaneO2Hit=false;
        }
        previousfLastPlane=fLastPlane;

        //check which stripend to set to true
        if (fStripend==1) lastPlaneO1Hit=true;
        if (fStripend==2) lastPlaneO2Hit=true;

        //set the hits in the other plane to be false
        lastPlaneE1Hit=false;
        lastPlaneE2Hit=false;

        MSG("CDAnalysis",Msg::kVerbose)
          <<"ev="<<event
          <<", last="<<fLastPlane
          <<", pl="<<fPlane
          <<", se="<<fStripend
          <<", E1="<<lastPlaneE1Hit
          <<", E2="<<lastPlaneE2Hit
          <<", O1="<<lastPlaneO1Hit
          <<", O2="<<lastPlaneO2Hit<<endl;
      }

      //add up the charge of the hits in each plane
      planeSigCor[fPlane]+=fChargeSigCor;
      planeSigCorT[fPlane]+=fChargeSigCor;
      planeMip[fPlane]+=fChargeMip;
      planeMipT[fPlane]+=fChargeMip;
      if (fChargePe>planeGreatestPe[fPlane]) planeGreatestPe[fPlane]=
                                               fChargePe;
      planePe[fPlane]+=fChargePe;
      if (fO1) planeSigCorO1[fPlane]+=fChargeSigCor;
      if (fO2) planeSigCorO2[fPlane]+=fChargeSigCor;
      if (fE1) planeSigCorE1[fPlane]+=fChargeSigCor;
      if (fE2) planeSigCorE2[fPlane]+=fChargeSigCor;
      if (fO1) planeSigCorO1T[fPlane]+=fChargeSigCor;
      if (fO2) planeSigCorO2T[fPlane]+=fChargeSigCor;
      if (fE1) planeSigCorE1T[fPlane]+=fChargeSigCor;
      if (fE2) planeSigCorE2T[fPlane]+=fChargeSigCor;

      if (fPlane==0) continue;
      //sum the sigcor in the event
      //make up for no plane 0 with 2 times plane 1
      else if (fPlane==1) sigCorTotal+=2*fChargeSigCor;
      else sigCorTotal+=fChargeSigCor;
    }
    //end of loop over the tracked hits

    //fill histos for all data
    hNumPlanes->Fill(numPlanesHit.size());
    hNumPlanesAll->Fill(numPlanesHitAll.size());
    hEventLength->Fill(fLastPlane+1);

    //make sure that the first plane makes sense and was actually set
    if (fFirstPlane<0 || fFirstPlane>59) continue;

    //calculate the event energy
    Double_t eventEn=(fScEnDepPerMEU/fSigCorPerMEU)*sigCorTotal*
      (1/fRatioScEnToBeamEn);
    //calc the last plane on the track with no xtalk
    Int_t lastPlaneNoXTalk=this->CalcLastPlaneOnTrkNoXTalk();

    //fill the prof
    pEnVsRange->Fill(lastPlaneNoXTalk+1,eventEn);

    if (fLastPlane%2==1) lastPlaneOddCounter++;
    if (fLastPlane%2==0) lastPlaneEvenCounter++;

    //work out which stripends in the last plane were hit
    if (lastPlaneE1Hit || lastPlaneE2Hit){
      if (lastPlaneE1Hit && !lastPlaneE2Hit) lastPlaneE1Counter++;
      else if (!lastPlaneE1Hit && lastPlaneE2Hit) lastPlaneE2Counter++;
      else lastPlaneEBothCounter++;
    }
    else if (lastPlaneO1Hit || lastPlaneO2Hit){
      if (lastPlaneO1Hit && !lastPlaneO2Hit) lastPlaneO1Counter++;
      else if (!lastPlaneO1Hit && lastPlaneO2Hit) lastPlaneO2Counter++;
      else lastPlaneOBothCounter++;
    }
    else{
      MSG("CDAnalysis",Msg::kInfo)
        <<"Nothing hit!"<<endl;

      MSG("CDAnalysis",Msg::kInfo)
        <<"ev="<<event
        <<", Num hits: trk="<<numTrkHits<<", unTrk="<<numUnTrkHits
        <<", xtalk="<<numXTalkHits
        <<", E1="<<lastPlaneE1Hit
        <<", E2="<<lastPlaneE2Hit
        <<", O1="<<lastPlaneO1Hit
        <<", O2="<<lastPlaneO2Hit<<endl;
    }

    MSG("CDAnalysis",Msg::kVerbose)<<"New event="<<event<<endl;

    totalCounter++;
    //look for the second of a pair of shared pmts
    if ((fLastPlane==26 || fLastPlane==27 ||
         fLastPlane==30 || fLastPlane==31 ||
         fLastPlane==34 || fLastPlane==35 ||
         fLastPlane==38 || fLastPlane==39 ||
         fLastPlane==42 || fLastPlane==43 ||
         fLastPlane==46 || fLastPlane==47 ||
         fLastPlane==50 || fLastPlane==51 ||
         fLastPlane==54 || fLastPlane==55 ||
         fLastPlane==58 || fLastPlane==59)){
      //fLastPlaneOdd
      //fLastPlaneEven
      
      if (this->IsDoubleEnded(fLastPlane)){
        doubleCounter++;
        if (this->IsPlaneXTalkOnly(fLastPlane)){
          doubleXTalkCounter++;
        }
      }
      else if (this->IsPlaneSideHit(fLastPlane,StripEnd::kEast)){
        side1Counter++;

        if (this->IsPlaneXTalkOnly(fLastPlane)){
          side1XTalkCounter++;
        }

        if (this->IsSharedPmtHit(fLastPlane)){
          sharedPmtCounter++;
          if (this->IsPlaneXTalkOnly(fLastPlane)){
            sharedPmtXTalkCounter++;
          }
        }
      }
      else if (this->IsPlaneSideHit(fLastPlane,StripEnd::kWest)){
        side2Counter++;
         
        if (this->IsPlaneXTalkOnly(fLastPlane)){
          side2XTalkCounter++;
        }
 
        if (this->IsSharedPmtHit(fLastPlane)){
          sharedPmtCounter++;
          if (this->IsPlaneXTalkOnly(fLastPlane)){
            sharedPmtXTalkCounter++;
          }
          MSG("CDAnalysis",Msg::kWarning)
            <<"This should not happen at CalDet!"<<endl;
        }
      }
      
      MSG("CDAnalysis",Msg::kVerbose)
        <<"Shared="<<sharedPmtCounter
        <<" ("<<sharedPmtXTalkCounter<<")"
        <<", 2end="<<doubleCounter
        <<" ("<<doubleXTalkCounter<<")"
        <<", end1="<<side1Counter
        <<" ("<<side1XTalkCounter<<")"
        <<", end2="<<side2Counter
        <<" ("<<side2XTalkCounter<<")"
        <<endl;
    }
    else{
      firstPlaneInPairCounter++;
    }

    if (!this->IsGenuineOrXTalk(fLastPlane)){
      notGenuineOrXTalkCounter++;
    }

    //make event length cuts
    if (this->CutOnEventLength()) continue;
    //make track quality cuts
    if (this->CutOnTrackQuality()) continue;

    MAXMSG("CDAnalysis",Msg::kInfo,10)
      <<"Event passing plane cuts="<<event<<endl;

    //fill the total sigcor histo
    hSigCorTotal->Fill(sigCorTotal);

    //fill event length plots for those passing cuts
    hFirstPlane->Fill(fFirstPlane);
    hEventLength2->Fill(fLastPlane+1);
    hEventLength3->Fill(lastPlaneNoXTalk+1);

    //look at a certain part of the spectrum
    if (lastPlaneNoXTalk+1>50 && lastPlaneNoXTalk+1<53){
      //fill the total sigcor histo
      hSigCorTotalTight->Fill(sigCorTotal);
    }
    
    if (lastPlaneNoXTalk+1==51){
      //fill the total sigcor histo
      hSigCorTotalPeak->Fill(sigCorTotal);
    }

    //look at pe versus range
    hPeVsRange->Fill(fLastPlane+1,planeGreatestPe[fLastPlane]);
    if (this->IsPlaneGenuine(fLastPlane)){
      hPeVsRangeGen->Fill(fLastPlane+1,planeGreatestPe[fLastPlane]);
    }
    else{
      hPeVsRangeXTalk->Fill(fLastPlane+1,planeGreatestPe[fLastPlane]);

      //now look at second from last plane
      if (this->IsPlaneGenuine(fLastPlane-1)){
        hPeVsRangeGen2->Fill(fLastPlane+1-1,
                               planeGreatestPe[fLastPlane-1]);
      }
      else{
        hPeVsRangeXTalk2->Fill(fLastPlane+1-1,
                               planeGreatestPe[fLastPlane-1]);

        //now look at third from last plane
        if (this->IsPlaneGenuine(fLastPlane-2)){
          hPeVsRangeGen3->Fill(fLastPlane+1-2,
                               planeGreatestPe[fLastPlane-2]);
        }
        else{
          hPeVsRangeXTalk3->Fill(fLastPlane+1-2,
                                 planeGreatestPe[fLastPlane-2]);
        }
      }
    }

    hCkvAdcVsRange->Fill(fLastPlane+1,fKovADC3);

    //section for response

    //const Int_t reNormPlane=70;
    //Int_t planeOffset=reNormPlane-fLastPlane;
    const Int_t os=fLastPlane;

    //do nasty hack for plane 40 - very nasty!!!!
    //if (fRunNumber==70709 || fRunNumber==70705){
    //planeSigCor[40]=planeSigCor[38];
      //don't do xtalk since if pl 38 has zero xtalk then you will make
      //the map bigger and will average a zero!!!!
      //planeSigCorX[40]=planeSigCorX[38];
      //planeSigCorT[40]=planeSigCorT[38];
    //}
    //do nasty hack for plane 35 - very nasty!!!!
    if (fRunNumber>40000 && fRunNumber<60000){
      planeSigCor[35]=planeSigCor[37];
      planeSigCorT[35]=planeSigCorT[37];
    }
    //do nasty hack for plane 35 - very nasty!!!!
    else if (fRunNumber>100000 && fRunNumber<110000){
      //map<Int_t,Float_t> tmp=planeSigCor;
      //map<Int_t,Float_t> tmpT=planeSigCor;
      //if (tmp[15]!=0){
        //planeSigCor[13]=planeSigCor[15];
      //}

      //if (tmpT[15]!=0){
        //planeSigCorT[13]=planeSigCorT[15];
      //}
    }

    //fill profs for all stripends
    this->FillProfHisto(pSigCorVsPlane,num,planeSigCor);
    this->FillProfHisto(pSigCorVsPlaneX,numX,planeSigCorX);
    this->FillProfHisto(pSigCorVsPlaneT,numT,planeSigCorT);
    //offset
    this->FillProfHisto(pSigCorVsOsPlane,numOs,planeSigCor,os);
    this->FillProfHisto(pSigCorVsOsPlaneX,numOsX,planeSigCorX,os);
    this->FillProfHisto(pSigCorVsOsPlaneT,numOsT,planeSigCorT,os);

    //fill non-offset profs for individual stripends
    this->FillProfHisto(pSigCorVsPlaneO1,numO1,planeSigCorO1);
    this->FillProfHisto(pSigCorVsPlaneO1X,numO1X,planeSigCorO1X);
    this->FillProfHisto(pSigCorVsPlaneO1T,numO1T,planeSigCorO1T);
    this->FillProfHisto(pSigCorVsPlaneO2,numO2,planeSigCorO2);
    this->FillProfHisto(pSigCorVsPlaneO2X,numO2X,planeSigCorO2X);
    this->FillProfHisto(pSigCorVsPlaneO2T,numO2T,planeSigCorO2T);
    this->FillProfHisto(pSigCorVsPlaneE1,numE1,planeSigCorE1);
    this->FillProfHisto(pSigCorVsPlaneE1X,numE1X,planeSigCorE1X);
    this->FillProfHisto(pSigCorVsPlaneE1T,numE1T,planeSigCorE1T);
    this->FillProfHisto(pSigCorVsPlaneE2,numE2,planeSigCorE2);
    this->FillProfHisto(pSigCorVsPlaneE2X,numE2X,planeSigCorE2X);
    this->FillProfHisto(pSigCorVsPlaneE2T,numE2T,planeSigCorE2T);

    //fill offset profs for individual stripends
    this->FillProfHisto(pSigCorVsOsPlaneO1,numOsO1,planeSigCorO1,os);
    this->FillProfHisto(pSigCorVsOsPlaneO1X,numOsO1X,planeSigCorO1X,os);
    this->FillProfHisto(pSigCorVsOsPlaneO1T,numOsO1T,planeSigCorO1T,os);
    this->FillProfHisto(pSigCorVsOsPlaneO2,numOsO2,planeSigCorO2,os);
    this->FillProfHisto(pSigCorVsOsPlaneO2X,numOsO2X,planeSigCorO2X,os);
    this->FillProfHisto(pSigCorVsOsPlaneO2T,numOsO2T,planeSigCorO2T,os);
    this->FillProfHisto(pSigCorVsOsPlaneE1,numOsE1,planeSigCorE1,os);
    this->FillProfHisto(pSigCorVsOsPlaneE1X,numOsE1X,planeSigCorE1X,os);
    this->FillProfHisto(pSigCorVsOsPlaneE1T,numOsE1T,planeSigCorE1T,os);
    this->FillProfHisto(pSigCorVsOsPlaneE2,numOsE2,planeSigCorE2,os);
    this->FillProfHisto(pSigCorVsOsPlaneE2X,numOsE2X,planeSigCorE2X,os);
    this->FillProfHisto(pSigCorVsOsPlaneE2T,numOsE2T,planeSigCorE2T,os);

    //section for total event sigcor and mips
    Float_t localSumSigCorTotal=0;
    Float_t localSumSigCorTotalCounter=0;
    for (map<Int_t,Float_t>::iterator sig=planeSigCorT.begin();
         sig!=planeSigCorT.end();++sig){
      localSumSigCorTotal+=sig->second;
      localSumSigCorTotalCounter++;
    }

    //add total of this event to the sum of all events 
    if (localSumSigCorTotalCounter!=0){
      sumSigCorTotal+=localSumSigCorTotal;
      sumSigCorTotalCounter++;
    }

    Float_t localSumMipTotal=0;
    Float_t localSumMipTotalCounter=0;
    for (map<Int_t,Float_t>::iterator sig=planeMipT.begin();
         sig!=planeMipT.end();++sig){
      localSumMipTotal+=sig->second;
      localSumMipTotalCounter++;
    }

    //add total of this event to the sum of all events 
    if (localSumMipTotalCounter!=0){
      sumMipTotal+=localSumMipTotal;
      sumMipTotalCounter++;
    }

    //section for sliding window
    Int_t minEventLength=35;
    if (fLastPlane<minEventLength) continue;

    //initialise the profile histos
    static Bool_t firstTime=true;
    if (firstTime){
      for (Int_t p=15;p>1;p--){
        vWindowSize.push_back(p);
        bigAv.push_back(0);
        bigAvCounter.push_back(0);
      }

      for (vector<Int_t>::iterator windowSize=vWindowSize.begin();
           windowSize!=vWindowSize.end();++windowSize){
        
        string sName=Form("%d",*windowSize);
        pSigCorVsDist.push_back
          (new TProfile(("pSigCorVsDist"+sName).c_str(),
                        ("pSigCorVsDist"+sName).c_str(),100,-2,40));
      }
    }
    firstTime=false;

    vector<TProfile*>::iterator prof=pSigCorVsDist.begin();
    vector<Float_t>::iterator av=bigAv.begin();
    vector<Int_t>::iterator count=bigAvCounter.begin();
    Float_t localSum15_18=0;
    Float_t localSum15_18Counter=0;
    Float_t localSum15_18Mip=0;
    Float_t localSum15_18MipCounter=0;

    //loop over the different window sizes
    for (vector<Int_t>::iterator windowSize=vWindowSize.begin();
         windowSize!=vWindowSize.end();++windowSize){
      
      //slide the window along from the end to 40-windowSize
      for (Int_t p=0;p<minEventLength-*windowSize;p++){
        Int_t windowStart=fLastPlane-p;

        MSG("CDAnalysis",Msg::kVerbose)
          <<"Window size="<<*windowSize
          <<", window start="<<windowStart
          <<endl;

        //loop through the hits in the window
        for (Int_t i=0;i<*windowSize;i++){
          Int_t plane=windowStart-i;
          
          //cut out planes that are greater than the last plane
          //in the view - ie the zeros at the end
          Bool_t doZeroCorrection=false;
          if (doZeroCorrection && 
              ((plane%2==0 && plane>fLastPlaneEven) || 
               (plane%2!=0 && plane>fLastPlaneOdd))){
            MSG("CDAnalysis",Msg::kVerbose)
              <<"Not using: Event="<<event
              <<", pl="<<plane
              <<", sigCor="<<planeSigCorT[plane]
              <<", fLastPlane="<<fLastPlane
              <<", fLastPlSig="<<planeSigCorT[fLastPlane]
              <<", fLastPlOdd="<<fLastPlaneOdd
              <<", fLastPlEven="<<fLastPlaneEven
              <<endl;
          }
          else{

            //calculate the best window value
            //if fLP=48 then the first plane included is 30
            //this is actually the 19th from the end of the track
            if (*windowSize==15 && windowStart==fLastPlane-18){
              localSum15_18+=planeSigCorT[plane];//total including xtalk
              localSum15_18Counter++;
              localSum15_18Mip+=planeMipT[plane];//total including xtalk
              localSum15_18MipCounter++;

              if (plane%2==0){//even
                
              }
              else if (plane%2==1){//odd

              }

            }

            (*prof)->Fill(p,planeSigCorT[plane]);
            (*av)+=planeSigCorT[plane];
            (*count)++;
          }
        }
      }
      av++;
      count++;
      prof++;
    }

    if (localSum15_18Counter!=0){
      h15_18All->Fill(localSum15_18);
      h15_18->Fill(localSum15_18/15);
      sum15_18+=localSum15_18;//add this event's energy on to total
      sum15_18Counter++;
      sum15_18Mip+=localSum15_18Mip;
      sum15_18MipCounter++;
    }

    if (localSum15_18Counter==0 || localSum15_18<200){
      MAXMSG("CDAnalysis",Msg::kInfo,40)
        <<"Event="<<event<<", strange 15_18="<<localSum15_18
        <<", counter="<<localSum15_18Counter<<endl;
    }

  }//end of for                                       
  

  MSG("CDAnalysis",Msg::kInfo)<<"Finished main loop"<<endl;

  //turn off the stats box printing
  gStyle->SetOptStat(0);

  //all stripends
  string sTitle="Energy Deposition (SigCor, all Stripends)";
  string sXTitle="Distance from end of track (planes)";
  this->DrawResponsePlot(sTitle,pSigCorVsPlane,pSigCorVsPlaneX,
                         pSigCorVsPlaneT,num,numX,numT);
  sTitle="Energy Deposition Vs Distance (SigCor, all Stripends)";
  this->DrawResponsePlot(sTitle,pSigCorVsOsPlane,pSigCorVsOsPlaneX,
                         pSigCorVsOsPlaneT,numOs,numOsX,numOsT,sXTitle);

  //non-offset planes
  sTitle="Energy Deposition (SigCor in Odd Planes FD stripend)";
  this->DrawResponsePlot(sTitle,pSigCorVsPlaneO1,pSigCorVsPlaneO1X,
                         pSigCorVsPlaneO1T,numO1,numO1X,numO1T);

  sTitle="Energy Deposition (SigCor in Odd Planes ND stripend)";
  this->DrawResponsePlot(sTitle,pSigCorVsPlaneO2,pSigCorVsPlaneO2X,
                         pSigCorVsPlaneO2T,numO2,numO2X,numO2T);

  sTitle="Energy Deposition (SigCor in Even Planes FD stripend)";
  this->DrawResponsePlot(sTitle,pSigCorVsPlaneE1,pSigCorVsPlaneE1X,
                         pSigCorVsPlaneE1T,numE1,numE1X,numE1T);

  sTitle="Energy Deposition (SigCor in Even Planes ND stripend)";
  this->DrawResponsePlot(sTitle,pSigCorVsPlaneE2,pSigCorVsPlaneE2X,
                         pSigCorVsPlaneE2T,numE2,numE2X,numE2T);

  //offset planes
  sTitle="Energy Deposition Vs Distance (SigCor in Odd Planes FD stripend)";
  this->DrawResponsePlot(sTitle,pSigCorVsOsPlaneO1,pSigCorVsOsPlaneO1X,
                         pSigCorVsOsPlaneO1T,numOsO1,numOsO1X,numOsO1T,
                         sXTitle);

  sTitle="Energy Deposition Vs Distance (SigCor in Odd Planes ND stripend)";
  this->DrawResponsePlot(sTitle,pSigCorVsOsPlaneO2,pSigCorVsOsPlaneO2X,
                         pSigCorVsOsPlaneO2T,numOsO2,numOsO2X,numOsO2T,
                         sXTitle);

  sTitle="Energy Deposition Vs Distance (SigCor in Even Planes FD stripend)";
  this->DrawResponsePlot(sTitle,pSigCorVsOsPlaneE1,pSigCorVsOsPlaneE1X,
                         pSigCorVsOsPlaneE1T,numOsE1,numOsE1X,numOsE1T,
                         sXTitle);

  sTitle="Energy Deposition Vs Distance (SigCor in Even Planes ND stripend)";
  this->DrawResponsePlot(sTitle,pSigCorVsOsPlaneE2,pSigCorVsOsPlaneE2X,
                         pSigCorVsOsPlaneE2T,numOsE2,numOsE2X,numOsE2T,
                         sXTitle);


  //turn on the stats box printing
  gStyle->SetOptStat(11);

  TCanvas *cProf=new TCanvas("cProf","Prof",0,0,1200,800);
  cProf->SetFillColor(0);
  string sName="SlidingWindow.ps";
  cProf->Print((sName+"[").c_str());
  gErrorIgnoreLevel=1;
  //plot the unnormalised prof first
  for (vector<TProfile*>::iterator prof=pSigCorVsDist.begin();
       prof!=pSigCorVsDist.end();++prof){      
    cProf->Clear();
    (*prof)->Draw();
    (*prof)->SetTitle("Average SigCor in Sliding Window");
    (*prof)->GetXaxis()->SetTitle("Start of sliding window (planes from end of track)");
    (*prof)->GetYaxis()->SetTitle("Average SigCor");
    (*prof)->GetXaxis()->CenterTitle();
    (*prof)->GetYaxis()->CenterTitle();
    cProf->Print(sName.c_str());
  }

  //now normalise the prof
  vector<Float_t>::iterator av=bigAv.begin();
  vector<Int_t>::iterator count=bigAvCounter.begin();
  for (vector<TProfile*>::iterator prof=pSigCorVsDist.begin();
       prof!=pSigCorVsDist.end();++prof){
      
    //normalise the prof and scale the axes
    MSG("CDAnalysis",Msg::kDebug)
      <<"bigCounter="<<*count<<", bigAv="<<*av<<endl;

    //pSigCorVsDistN.push_back
    //(new TProfile(("pSigCorVsDistN"+sName).c_str(),
    //      ("pSigCorVsDistN"+sName).c_str(),100,-2,40));
    if (*av>0) (*prof)->Scale((*count)/(*av));
    (*prof)->SetMaximum(1.2);
    (*prof)->SetMinimum(0.8);
      
    cProf->Clear();
    (*prof)->Draw();
    (*prof)->SetTitle("Average SigCor in Sliding Window");
    (*prof)->GetXaxis()->SetTitle("Start of sliding window (planes from end of track)");
    (*prof)->GetYaxis()->SetTitle("Average SigCor");
    (*prof)->GetXaxis()->CenterTitle();
    (*prof)->GetYaxis()->CenterTitle();
    cProf->Print(sName.c_str());
    av++;
    count++;
  }
  //close the file
  gErrorIgnoreLevel=0;
  cProf->Print((sName+"]").c_str());

  //now do a few profs on the same plot
  TCanvas *cProf2=new TCanvas("cProf2","Prof",0,0,1200,800);
  cProf2->SetFillColor(0);
  sName="SlidingWindowAll.ps";
  //cProf2->Print((sName+"[").c_str());
  Int_t counter=0;
  for (vector<TProfile*>::iterator prof=pSigCorVsDist.begin();
       prof!=pSigCorVsDist.end();++prof){

    (*prof)->SetMaximum(1.09);
    (*prof)->SetMinimum(0.96);

    //count which prof got to
    counter++;

    //decide which profs to draw
    static Bool_t firstTime=true;
    if (counter==1 || counter==6 || counter==11){ 
      if (firstTime) {
        (*prof)->Draw();
        firstTime=false;
      }
      else  (*prof)->Draw("sames");

      (*prof)->SetTitle("Average SigCor in Sliding Window");
      (*prof)->GetXaxis()->SetTitle("Start of sliding window (planes from end of track)");
      (*prof)->GetYaxis()->SetTitle("Average SigCor");
      (*prof)->GetXaxis()->CenterTitle();
      (*prof)->GetYaxis()->CenterTitle();
      static Int_t color=2;
      (*prof)->SetLineColor(color);
      color+=2;
    }
  }

  //turn off the stats box printing
  //gStyle->SetOptStat(0);
  gStyle->SetOptStat(1111111);

  TCanvas *cGeom=new TCanvas("cGeom","Geom",0,0,1200,800);
  cGeom->SetFillColor(0);
  cGeom->Divide(2,3);
  cGeom->cd(1);
  hEventLength2->Draw();
  cGeom->cd(2);
  hEventLength3->Draw();
  cGeom->cd(3);
  hFirstPlane->Draw();
  cGeom->cd(4);
  hEventLength->Draw();
  cGeom->cd(5);
  hNumPlanes->Draw();
  cGeom->cd(6);
  hNumPlanesAll->Draw();
  hNumPlanesMu->Draw("same");
 
  /*
  TCanvas *cCkv=new TCanvas("cCkv","Ckv",0,0,1200,800);
  cCkv->SetFillColor(0);
  cCkv->Divide(2,2);
  cCkv->cd(1);
  hCkvAdcVsRange->Draw("colz");
  hCkvAdcVsRange->Draw("profsame");
  cCkv->cd(2);
  TProfile* p=hCkvAdcVsRange->ProfileX();
  p->Draw();
  cCkv->cd(3);
  hCkvAdcVsRange->Draw("colz");
  p->Draw("sames");
  */

  TCanvas *c15_18=new TCanvas("c15_18","15_18",0,0,1200,800);
  c15_18->SetFillColor(0);
  c15_18->Divide(2,3);
  c15_18->cd(1);
  h15_18->Draw();
  c15_18->cd(2);
  h15_18All->Draw();
  c15_18->cd(3);
  hSigCorTotal->Draw();
  c15_18->cd(4);
  hSigCorTotalTight->Draw();
  c15_18->cd(5);
  hSigCorTotalPeak->Draw();
  c15_18->cd(6);
  pEnVsRange->Draw();

  MSG("CDAnalysis",Msg::kInfo) 
    <<endl
    <<"Mean sigCorTotal="<<hSigCorTotal->GetMean()<<endl
    <<"Ratio ScEn to BeamEn="<<fRatioScEnToBeamEn
    <<" ("<<1/fRatioScEnToBeamEn<<")"<<endl
    <<"GeV per SigCor="<<fScEnDepPerMEU/fSigCorPerMEU<<endl
    <<"Beam energy from Calorimetry (lose)="
    <<((fScEnDepPerMEU/fSigCorPerMEU)*
    hSigCorTotal->GetMean()*(1/fRatioScEnToBeamEn))/
    Munits::gigaelectronvolt<<" "<<Munits::base_energy_name<<endl
    <<"Beam energy from Calorimetry (tight)="
    <<((fScEnDepPerMEU/fSigCorPerMEU)*
    hSigCorTotalTight->GetMean()*(1/fRatioScEnToBeamEn))/
    Munits::gigaelectronvolt<<" "<<Munits::base_energy_name<<endl
    <<"Beam energy from Calorimetry (peak)="
    <<((fScEnDepPerMEU/fSigCorPerMEU)*
    hSigCorTotalPeak->GetMean()*(1/fRatioScEnToBeamEn))/
    Munits::gigaelectronvolt<<" "<<Munits::base_energy_name<<endl<<endl;

  /*
  TCanvas *cXtalk=new TCanvas("cXtalk","Xtalk",0,0,1200,800);
  cXtalk->SetFillColor(0);
  cXtalk->Divide(3,3);
  cXtalk->cd(1);
  hPeVsRange->Draw("colz");
  cXtalk->cd(2);
  hPeVsRangeGen->Draw("colz");
  cXtalk->cd(3);
  hPeVsRangeXTalk->Draw("colz");

  cXtalk->cd(5);
  hPeVsRangeGen2->Draw("colz");
  cXtalk->cd(6);
  hPeVsRangeXTalk2->Draw("colz");
  cXtalk->cd(8);
  hPeVsRangeGen3->Draw("colz");
  cXtalk->cd(9);
  hPeVsRangeXTalk3->Draw("colz");
  */

  Int_t oddEvenSum=lastPlaneEvenCounter+lastPlaneOddCounter;
  Float_t oddPercent=-1;
  Float_t evenPercent=-1;
  if (oddEvenSum) oddPercent=100.*lastPlaneOddCounter/oddEvenSum;
  if (oddEvenSum) evenPercent=100.*lastPlaneEvenCounter/oddEvenSum;

  MSG("CDAnalysis",Msg::kInfo) 
    <<"Last plane, odd="<<lastPlaneOddCounter<<" ("<<oddPercent<<"%)"
    <<", even="<<lastPlaneEvenCounter<<" ("<<evenPercent<<"%)"<<endl;

  Int_t totalSum=lastPlaneE1Counter+lastPlaneE2Counter+
    lastPlaneO1Counter+lastPlaneO2Counter+
    lastPlaneEBothCounter+lastPlaneOBothCounter;
  Float_t e1Percent=-1;
  Float_t e2Percent=-1;
  Float_t o1Percent=-1;
  Float_t o2Percent=-1;
  Float_t eBothPercent=-1;
  Float_t oBothPercent=-1;
  if (totalSum) o1Percent=100.*lastPlaneO1Counter/totalSum;
  if (totalSum) o2Percent=100.*lastPlaneO2Counter/totalSum;
  if (totalSum) e1Percent=100.*lastPlaneE1Counter/totalSum;
  if (totalSum) e2Percent=100.*lastPlaneE2Counter/totalSum;
  if (totalSum) eBothPercent=100.*lastPlaneEBothCounter/totalSum;
  if (totalSum) oBothPercent=100.*lastPlaneOBothCounter/totalSum;

  Float_t totalPercent=o1Percent+o2Percent+e1Percent+e2Percent+
    oBothPercent+eBothPercent;

  MSG("CDAnalysis",Msg::kInfo) 
    <<"Last Plane:"<<endl
    <<"  Total number hits="<<totalSum<<endl
    <<"  Hits on both sides of plane:"<<endl
    <<"    Odd="<<lastPlaneOBothCounter<<" ("<<oBothPercent<<"%)"<<endl
    <<"    Even="<<lastPlaneEBothCounter<<" ("<<eBothPercent<<"%)"<<endl
    <<"  Hits on only one side of plane:"<<endl
    <<"    Odd-ND="<<lastPlaneO2Counter<<" ("<<o2Percent<<"%)"<<endl
    <<"    Odd-FD="<<lastPlaneO1Counter<<" ("<<o1Percent<<"%)"<<endl
    <<"    Even-ND="<<lastPlaneE2Counter<<" ("<<e2Percent<<"%)"<<endl
    <<"    Even-FD="<<lastPlaneE1Counter<<" ("<<e1Percent<<"%)"<<endl
    <<"  Sum of 6 percentages="<<totalPercent<<endl;

  Float_t sharedXTalkPercent=-1;
  Float_t side1XTalkPercent=-1;
  Float_t side2XTalkPercent=-1;
  Float_t doubleXTalkPercent=-1;
  Float_t notGenuineOrXTalkPercent=-1;
  if (sharedPmtCounter) sharedXTalkPercent=100.*sharedPmtXTalkCounter/
                          sharedPmtCounter;
  if (side1Counter) side1XTalkPercent=100.*side1XTalkCounter/
                      side1Counter;
  if (side2Counter) side2XTalkPercent=100.*side2XTalkCounter/
                      side2Counter;
  if (doubleCounter) doubleXTalkPercent=100.*doubleXTalkCounter/
                      doubleCounter;
  if (totalCounter) notGenuineOrXTalkPercent=100.*
                      notGenuineOrXTalkCounter/totalCounter;

  MSG("CDAnalysis",Msg::kInfo) 
    <<endl
    <<"totalCounter="<<totalCounter<<endl
    <<"firstPlaneInPairCounter="<<firstPlaneInPairCounter<<endl
    <<"double="<<doubleCounter
    <<", xtalk="<<doubleXTalkCounter
    <<" ("<<doubleXTalkPercent<<"%)"<<endl
    <<"side1="<<side1Counter
    <<", xtalk="<<side1XTalkCounter
    <<" ("<<side1XTalkPercent<<"%)"<<endl
    <<"side2="<<side2Counter
    <<", xtalk="<<side2XTalkCounter
    <<" ("<<side2XTalkPercent<<"%)"<<endl
    <<"Number of times shared pmt on second of"
    <<" two planes was last plane="
    <<sharedPmtCounter<<endl
    <<"Number of times above was XTalk only="<<sharedPmtXTalkCounter
    <<" ("<<sharedXTalkPercent<<"%)"<<endl
    <<"notGenuineOrXTalkCounter="<<notGenuineOrXTalkCounter
    <<" ("<<notGenuineOrXTalkPercent<<"%)"<<endl;

  MSG("CDAnalysis",Msg::kInfo) 
    <<endl<<"passPidCounter="<<passPidCounter<<endl
    <<"noTrackCounter="<<noTrackCounter
    <<" ("<<Pct(noTrackCounter,passPidCounter)<<"%)"<<endl;
   
  if (sum15_18Counter!=0){
    MSG("CDAnalysis",Msg::kInfo) 
      <<endl
      <<" ** Calibration Summary **"<<endl
      <<"SigCor in sliding window 15_18:"<<endl
      <<"    SigCor dep="<<sum15_18<<" (N="<<sum15_18Counter<<")"<<endl
      <<"    Av SigCor dep="<<sum15_18/sum15_18Counter
      <<" +/- "<<100*(h15_18->GetRMS()/sqrt(h15_18->GetEntries()))/
      h15_18->GetMean()<<"%"<<endl
      <<"    Av SigCor dep per plane="<<(sum15_18/sum15_18Counter)/15
      <<" +/- "<<100*(h15_18->GetRMS()/sqrt(h15_18->GetEntries()))/
      h15_18->GetMean()<<"%"<<endl
      <<"MIPs in sliding window 15_18:"<<endl
      <<"    MIPs dep="<<sum15_18Mip
      <<" (N="<<sum15_18MipCounter<<")"<<endl
      <<"    Av MIPs dep="<<sum15_18Mip/sum15_18MipCounter
      <<endl;
    if (sumSigCorTotalCounter!=0){
      MSG("CDAnalysis",Msg::kInfo) 
        <<"SigCor of whole event:"<<endl
        <<"    Sum of SigCor dep="<<sumSigCorTotal
        <<" (N="<<sumSigCorTotalCounter<<")"<<endl
        <<"    Av sigCor dep="
        <<sumSigCorTotal/sumSigCorTotalCounter<<endl
        <<"    1 sigCor="
        <<1800/(sumSigCorTotal/sumSigCorTotalCounter)<<" MeV"
        <<endl;
    }
    if (sumMipTotalCounter!=0){
      MSG("CDAnalysis",Msg::kInfo) 
        <<"MIPs of whole event:"<<endl
        <<"    Sum of MIPs dep="<<sumMipTotal
        <<" (N="<<sumMipTotalCounter<<")"<<endl
        <<"    Av MIPs dep="
        <<sumMipTotal/sumMipTotalCounter<<endl
        <<"    1 MIP="
        <<1800/(sumMipTotal/sumMipTotalCounter)<<" MeV"<<endl
        <<"    1 GeV="
        <<(sumMipTotal/sumMipTotalCounter)/1.8<<" MIPs"
        <<endl<<endl;
    }
  }

  MSG("CDAnalysis",Msg::kInfo) 
    <<" ** Finished MuonResponse method **"<<endl;
}

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

void CDAnalysis::MuonNearFar()
{
  MSG("CDAnalysis",Msg::kInfo) 
    <<" ** Running MuonNearFar method... **"<<endl;

  //open the output file for the histograms
  fOutFile=this->OpenFile(fRunNumber,"MuonNF");

  TH1F *h15_18All=new TH1F("h15_18All","h15_18All",1000,-2,75);
  h15_18All->GetXaxis()->SetTitle("SigCors");
  h15_18All->GetXaxis()->CenterTitle();
  h15_18All->GetYaxis()->SetTitle("");
  h15_18All->GetYaxis()->CenterTitle();
  h15_18All->SetFillColor(0);
  h15_18All->SetBit(TH1::kCanRebin);

  TH1F *h15_18=new TH1F("h15_18","h15_18",1000,-2,75);
  h15_18->GetXaxis()->SetTitle("SigCors");
  h15_18->GetXaxis()->CenterTitle();
  h15_18->GetYaxis()->SetTitle("");
  h15_18->GetYaxis()->CenterTitle();
  h15_18->SetFillColor(0);
  h15_18->SetBit(TH1::kCanRebin);

  TH1F *h15_18_O1=new TH1F("h15_18_O1","h15_18_O1",1000,-2,75);
  h15_18_O1->GetXaxis()->SetTitle("SigCors");
  h15_18_O1->GetXaxis()->CenterTitle();
  h15_18_O1->GetYaxis()->SetTitle("");
  h15_18_O1->GetYaxis()->CenterTitle();
  h15_18_O1->SetFillColor(0);
  h15_18_O1->SetBit(TH1::kCanRebin);

  TH1F *h15_18_E1=new TH1F("h15_18_E1","h15_18_E1",1000,-2,75);
  h15_18_E1->GetXaxis()->SetTitle("SigCors");
  h15_18_E1->GetXaxis()->CenterTitle();
  h15_18_E1->GetYaxis()->SetTitle("");
  h15_18_E1->GetYaxis()->CenterTitle();
  h15_18_E1->SetFillColor(0);
  h15_18_E1->SetBit(TH1::kCanRebin);

  TH1F *h15_18_O2=new TH1F("h15_18_O2","h15_18_O2",1000,-2,75);
  h15_18_O2->GetXaxis()->SetTitle("SigCors");
  h15_18_O2->GetXaxis()->CenterTitle();
  h15_18_O2->GetYaxis()->SetTitle("");
  h15_18_O2->GetYaxis()->CenterTitle();
  h15_18_O2->SetFillColor(0);
  h15_18_O2->SetBit(TH1::kCanRebin);

  TH1F *h15_18_E2=new TH1F("h15_18_E2","h15_18_E2",1000,-2,75);
  h15_18_E2->GetXaxis()->SetTitle("SigCors");
  h15_18_E2->GetXaxis()->CenterTitle();
  h15_18_E2->GetYaxis()->SetTitle("");
  h15_18_E2->GetYaxis()->CenterTitle();
  h15_18_E2->SetFillColor(0);
  h15_18_E2->SetBit(TH1::kCanRebin);

  //sum of all strip ends, normal and offset planes
  //tracked
  TProfile* pSigCorVsPlane=new TProfile("pSigCorVsPlane",
                                        "pSigCorVsPlane",
                                        62,0,62);
  TProfile* pSigCorVsOsPlane=new TProfile("pSigCorVsOsPlane",
                                          "pSigCorVsOsPlane",
                                          62,0,62);

  //xtalk + untracked hits
  TProfile* pSigCorVsPlaneX=new TProfile("pSigCorVsPlaneX",
                                         "pSigCorVsPlaneX",
                                         62,0,62);
  TProfile* pSigCorVsOsPlaneX=new TProfile("pSigCorVsOsPlaneX",
                                           "pSigCorVsOsPlaneX",
                                           62,0,62);

  //total
  TProfile* pSigCorVsPlaneT=new TProfile("pSigCorVsPlaneT",
                                         "pSigCorVsPlaneT",
                                         62,0,62);
  TProfile* pSigCorVsOsPlaneT=new TProfile("pSigCorVsOsPlaneT",
                                           "pSigCorVsOsPlaneT",
                                           62,0,62);

  //individual strip ends
  //tracked
  TProfile* pSigCorVsPlaneO1=new TProfile("pSigCorVsPlaneO1",
                                          "pSigCorVsPlaneO1",
                                          62,0,62);
  TProfile* pSigCorVsPlaneO2=new TProfile("pSigCorVsPlaneO2",
                                          "pSigCorVsPlaneO2",
                                          62,0,62);
  TProfile* pSigCorVsPlaneE1=new TProfile("pSigCorVsPlaneE1",
                                          "pSigCorVsPlaneE1",
                                          62,0,62);
  TProfile* pSigCorVsPlaneE2=new TProfile("pSigCorVsPlaneE2",
                                          "pSigCorVsPlaneE2",
                                          62,0,62);
  //xtalk + untracked hits
  TProfile* pSigCorVsPlaneO1X=new TProfile("pSigCorVsPlaneO1X",
                                          "pSigCorVsPlaneO1X",
                                          62,0,62);
  TProfile* pSigCorVsPlaneO2X=new TProfile("pSigCorVsPlaneO2X",
                                          "pSigCorVsPlaneO2X",
                                          62,0,62);
  TProfile* pSigCorVsPlaneE1X=new TProfile("pSigCorVsPlaneE1X",
                                          "pSigCorVsPlaneE1X",
                                          62,0,62);
  TProfile* pSigCorVsPlaneE2X=new TProfile("pSigCorVsPlaneE2X",
                                          "pSigCorVsPlaneE2X",
                                          62,0,62);
  //total
  TProfile* pSigCorVsPlaneO1T=new TProfile("pSigCorVsPlaneO1T",
                                          "pSigCorVsPlaneO1T",
                                          62,0,62);
  TProfile* pSigCorVsPlaneO2T=new TProfile("pSigCorVsPlaneO2T",
                                          "pSigCorVsPlaneO2T",
                                          62,0,62);
  TProfile* pSigCorVsPlaneE1T=new TProfile("pSigCorVsPlaneE1T",
                                          "pSigCorVsPlaneE1T",
                                          62,0,62);
  TProfile* pSigCorVsPlaneE2T=new TProfile("pSigCorVsPlaneE2T",
                                          "pSigCorVsPlaneE2T",
                                          62,0,62);

  //individual strip ends - offset planes
  //tracked
  TProfile* pSigCorVsOsPlaneO1=new TProfile("pSigCorVsOsPlaneO1",
                                            "pSigCorVsOsPlaneO1",
                                            62,0,62);
  TProfile* pSigCorVsOsPlaneO2=new TProfile("pSigCorVsOsPlaneO2",
                                            "pSigCorVsOsPlaneO2",
                                            62,0,62);
  TProfile* pSigCorVsOsPlaneE1=new TProfile("pSigCorVsOsPlaneE1",
                                            "pSigCorVsOsPlaneE1",
                                            62,0,62);
  TProfile* pSigCorVsOsPlaneE2=new TProfile("pSigCorVsOsPlaneE2",
                                            "pSigCorVsOsPlaneE2",
                                            62,0,62);

  //xtalk + untracked hits
  TProfile* pSigCorVsOsPlaneO1X=new TProfile("pSigCorVsOsPlaneO1X",
                                             "pSigCorVsOsPlaneO1X",
                                             62,0,62);
  TProfile* pSigCorVsOsPlaneO2X=new TProfile("pSigCorVsOsPlaneO2X",
                                            "pSigCorVsOsPlaneO2X",
                                            62,0,62);
  TProfile* pSigCorVsOsPlaneE1X=new TProfile("pSigCorVsOsPlaneE1X",
                                            "pSigCorVsOsPlaneE1X",
                                            62,0,62);
  TProfile* pSigCorVsOsPlaneE2X=new TProfile("pSigCorVsOsPlaneE2X",
                                            "pSigCorVsOsPlaneE2X",
                                            62,0,62);
  //total
  TProfile* pSigCorVsOsPlaneO1T=new TProfile("pSigCorVsOsPlaneO1T",
                                            "pSigCorVsOsPlaneO1T",
                                            62,0,62);
  TProfile* pSigCorVsOsPlaneO2T=new TProfile("pSigCorVsOsPlaneO2T",
                                            "pSigCorVsOsPlaneO2T",
                                            62,0,62);
  TProfile* pSigCorVsOsPlaneE1T=new TProfile("pSigCorVsOsPlaneE1T",
                                            "pSigCorVsOsPlaneE1T",
                                            62,0,62);
  TProfile* pSigCorVsOsPlaneE2T=new TProfile("pSigCorVsOsPlaneE2T",
                                            "pSigCorVsOsPlaneE2T",
                                            62,0,62);


  //near/far differences
  TProfile* pNearFarOdd=new TProfile("pNearFarOdd",
                                     "pNearFarOdd",62,0,62);
  TProfile* pNearFarEven=new TProfile("pNearFarEven",
                                      "pNearFarEven",62,0,62);


  vector<TProfile*> pSigCorVsDist;
  vector<Int_t> vWindowSize;
  vector<Float_t> bigAv;
  vector<Int_t> bigAvCounter;
  Float_t sum15_18=0;
  Float_t sum15_18Counter=0;
  Float_t sum15_18Mip=0;
  Float_t sum15_18MipCounter=0;
  Float_t sumSigCorTotal=0;
  Float_t sumSigCorTotalCounter=0;
  Float_t sumMipTotal=0;
  Float_t sumMipTotalCounter=0;

  //variables to calculate the energy dep for different stripends
  Float_t sum15_18_O1=0;
  Float_t sum15_18CounterO1=0;
  Float_t sum15_18_O2=0;
  Float_t sum15_18CounterO2=0;
  Float_t sum15_18_E1=0;
  Float_t sum15_18CounterE1=0;
  Float_t sum15_18_E2=0;
  Float_t sum15_18CounterE2=0;

  //num maps for all stripends
  map<Int_t,Float_t> num;
  map<Int_t,Float_t> numX;
  map<Int_t,Float_t> numT;
  map<Int_t,Float_t> numOs;
  map<Int_t,Float_t> numOsX;
  map<Int_t,Float_t> numOsT;

  //non-offset num maps for individual stripends
  map<Int_t,Float_t> numO1;
  map<Int_t,Float_t> numO1X;
  map<Int_t,Float_t> numO1T;
  map<Int_t,Float_t> numO2;
  map<Int_t,Float_t> numO2X;
  map<Int_t,Float_t> numO2T;
  map<Int_t,Float_t> numE1;
  map<Int_t,Float_t> numE1X;
  map<Int_t,Float_t> numE1T;
  map<Int_t,Float_t> numE2;
  map<Int_t,Float_t> numE2X;
  map<Int_t,Float_t> numE2T;

  //offset num maps for individual stripends
  map<Int_t,Float_t> numOsO1;
  map<Int_t,Float_t> numOsO1X;
  map<Int_t,Float_t> numOsO1T;
  map<Int_t,Float_t> numOsO2;
  map<Int_t,Float_t> numOsO2X;
  map<Int_t,Float_t> numOsO2T;
  map<Int_t,Float_t> numOsE1;
  map<Int_t,Float_t> numOsE1X;
  map<Int_t,Float_t> numOsE1T;
  map<Int_t,Float_t> numOsE2;
  map<Int_t,Float_t> numOsE2X;
  map<Int_t,Float_t> numOsE2T;

  Int_t noTrackCounter=0;
  Int_t passPidCounter=0;
 
  
  this->InitialiseLoopVariables();  
  
  for(Int_t event=0;event<fEvents;event++){
    
    this->SetLoopVariables(event);

    if (this->CutOnDeadChips()) continue;
    if (this->CutOnPid()) continue;

    //cut on not straight tracks
    //Double_t mainX1Cut=0.5;
    //if (!this->IsStraightTrack(mainX1Cut)){
    //continue;
    //}

    //count the number of events passing the pid
    passPidCounter++;

    map<Int_t,Int_t> numPlanesHit;
    map<Int_t,Int_t> numPlanesHit1;
    map<Int_t,Int_t> numPlanesHit2;
    map<Int_t,Int_t> numPlanesHitAll;
    map<Int_t,Float_t> planeSigCor;
    map<Int_t,Float_t> planeSigCorX;
    map<Int_t,Float_t> planeSigCorT;
    map<Int_t,Float_t> planeSigCorO1;    
    map<Int_t,Float_t> planeSigCorO1X;
    map<Int_t,Float_t> planeSigCorO1T;
    map<Int_t,Float_t> planeSigCorO2;
    map<Int_t,Float_t> planeSigCorO2X;
    map<Int_t,Float_t> planeSigCorO2T;
    map<Int_t,Float_t> planeSigCorE1;
    map<Int_t,Float_t> planeSigCorE1X;
    map<Int_t,Float_t> planeSigCorE1T;
    map<Int_t,Float_t> planeSigCorE2;
    map<Int_t,Float_t> planeSigCorE2X;
    map<Int_t,Float_t> planeSigCorE2T;

    map<Int_t,Float_t> planeMip;
    map<Int_t,Float_t> planeMipX;
    map<Int_t,Float_t> planeMipT;

    map<Int_t,Float_t> planePe;
    map<Int_t,Float_t> planeGreatestPe;

    //get tclones arrays for this snarl
    TClonesArray &cTrk=*fTrkHitInfo;
    Int_t numTrkHits=fTrkHitInfo->GetEntries();
    TClonesArray &cUnTrk = *fUnTrkHitInfo;
    Int_t numUnTrkHits = fUnTrkHitInfo->GetEntries();
    //CDTrackedHitInfo *trackedHitInfo=(CDTrackedHitInfo*) c[hit];
    TClonesArray &cXTalk = *fXTalkHits;
    Int_t numXTalkHits=fXTalkHits->GetEntries();
    //TClonesArray &cTruth = *fTruthHitInfo;
    Int_t numTruthHits=0;
    if (fTruthHitInfo) numTruthHits=fTruthHitInfo->GetEntries();

    //look for untracked big events
    if (numTrkHits<2 && numXTalkHits+numUnTrkHits>1){
      MSG("CDAnalysis",Msg::kDebug)
        <<"Event="<<event
        <<", Num hits: trk="<<numTrkHits<<", unTrk="<<numUnTrkHits
        <<", xtalk="<<numXTalkHits<<", truth="<<numTruthHits
        <<endl;
      noTrackCounter++;
      //continue;
    }

    //loop over the untracked hits in the snarl
    for (Int_t hit=0;hit<numUnTrkHits;hit++){
      //cast the tclonesarray up to a trackedhitinfo object
      CDTrackedHitInfo *hitInfo=
        dynamic_cast<CDTrackedHitInfo*>(cUnTrk[hit]);

      this->ReadInHitInfo(hitInfo);

      if (fPlane==0 && fChargeSigCor>1000){
        MSG("CDAnalysis",Msg::kDebug)
          <<"XTalk Event="<<event
          <<", fPlane="<<fPlane
          <<", fStrip="<<fStrip
          <<", fStripend="<<fStripend
          <<", fChargeAdc="<<fChargeAdc
          <<", fSigCor="<<fChargeSigCor<<endl;
      }

      //cut out the bad channels
      if (this->CutOnBadPedestals(fPlane,fStrip,fStripend)) continue;

      //count the number of planes hit
      numPlanesHitAll[fPlane]++;

      //check the channel's sanity
      this->StandardSanityChecks();

      //add up the charge of the hits in each plane
      planeSigCorX[fPlane]+=fChargeSigCor;
      planeSigCorT[fPlane]+=fChargeSigCor;
      planeMipX[fPlane]+=fChargeMip;
      planeMipT[fPlane]+=fChargeMip;
      if (fO1) planeSigCorO1X[fPlane]+=fChargeSigCor;
      if (fO2) planeSigCorO2X[fPlane]+=fChargeSigCor;
      if (fE1) planeSigCorE1X[fPlane]+=fChargeSigCor;
      if (fE2) planeSigCorE2X[fPlane]+=fChargeSigCor;
      if (fO1) planeSigCorO1T[fPlane]+=fChargeSigCor;
      if (fO2) planeSigCorO2T[fPlane]+=fChargeSigCor;
      if (fE1) planeSigCorE1T[fPlane]+=fChargeSigCor;
      if (fE2) planeSigCorE2T[fPlane]+=fChargeSigCor;

      if (fChargePe>planeGreatestPe[fPlane]) planeGreatestPe[fPlane]=
                                               fChargePe;
    }

    //loop over xtalk hits in the snarl
    for (Int_t hit=0;hit<numXTalkHits;hit++){
      CDXTalkHitInfo *hitInfo=
        dynamic_cast<CDXTalkHitInfo*>(cXTalk[hit]);

      this->ReadInHitInfo(hitInfo);
      
      //cut out the bad channels
      if (this->CutOnBadPedestals(fPlane,fStrip,fStripend)) continue;

      //count the number of planes hit
      numPlanesHitAll[fPlane]++;

      //check the channel's sanity
      this->StandardSanityChecks();

      //add up the charge of the hits in each plane
      planeSigCorX[fPlane]+=fChargeSigCor;
      planeSigCorT[fPlane]+=fChargeSigCor;
      planeMipX[fPlane]+=fChargeMip;
      planeMipT[fPlane]+=fChargeMip;
      if (fO1) planeSigCorO1X[fPlane]+=fChargeSigCor;
      if (fO2) planeSigCorO2X[fPlane]+=fChargeSigCor;
      if (fE1) planeSigCorE1X[fPlane]+=fChargeSigCor;
      if (fE2) planeSigCorE2X[fPlane]+=fChargeSigCor;
      if (fO1) planeSigCorO1T[fPlane]+=fChargeSigCor;
      if (fO2) planeSigCorO2T[fPlane]+=fChargeSigCor;
      if (fE1) planeSigCorE1T[fPlane]+=fChargeSigCor;
      if (fE2) planeSigCorE2T[fPlane]+=fChargeSigCor;

      if (fChargePe>planeGreatestPe[fPlane]) planeGreatestPe[fPlane]=
                                               fChargePe;
    }

    Bool_t lastPlaneO1Hit=false;
    Bool_t lastPlaneO2Hit=false;
    Bool_t lastPlaneE1Hit=false;
    Bool_t lastPlaneE2Hit=false;
    Int_t previousfLastPlane=-1;

    //loop over the tracked hits
    for (Int_t hit=0;hit<numTrkHits;hit++){
      CDTrackedHitInfo *trackedHitInfo=
        dynamic_cast<CDTrackedHitInfo*>(cTrk[hit]);
          
      this->ReadInHitInfo(trackedHitInfo);

      //cut out the bad channels
      if (this->CutOnBadPedestals(fPlane,fStrip,fStripend)) continue;

      //check the channel's sanity
      this->StandardSanityChecks();

      //calculate the first and last planes
      this->CalcFirstLastPlane(fPlane);

      //count the number of planes hit
      numPlanesHit[fPlane]++;
      numPlanesHitAll[fPlane]++;

     if (fLastPlane==fPlane && fPlane%2==0){

        //check if a new last plane has been found
        if (fLastPlane!=previousfLastPlane){
          lastPlaneE1Hit=false;
          lastPlaneE2Hit=false;
        }
        previousfLastPlane=fLastPlane;

        //check which stripend to set to true
        if (fStripend==1) lastPlaneE1Hit=true;
        if (fStripend==2) lastPlaneE2Hit=true;

        //set the hits in the other plane to be false
        lastPlaneO1Hit=false;
        lastPlaneO2Hit=false;
        MSG("CDAnalysis",Msg::kVerbose)
          <<"ev="<<event
          <<", last="<<fLastPlane
          <<", pl="<<fPlane
          <<", se="<<fStripend
          <<", E1="<<lastPlaneE1Hit
          <<", E2="<<lastPlaneE2Hit
          <<", O1="<<lastPlaneO1Hit
          <<", O2="<<lastPlaneO2Hit<<endl;
      }
      else if (fLastPlane==fPlane && fPlane%2==1){

        //check if a new last plane has been found
        if (fLastPlane!=previousfLastPlane){
          lastPlaneO1Hit=false;
          lastPlaneO2Hit=false;
        }
        previousfLastPlane=fLastPlane;

        //check which stripend to set to true
        if (fStripend==1) lastPlaneO1Hit=true;
        if (fStripend==2) lastPlaneO2Hit=true;

        //set the hits in the other plane to be false
        lastPlaneE1Hit=false;
        lastPlaneE2Hit=false;

        MSG("CDAnalysis",Msg::kVerbose)
          <<"ev="<<event
          <<", last="<<fLastPlane
          <<", pl="<<fPlane
          <<", se="<<fStripend
          <<", E1="<<lastPlaneE1Hit
          <<", E2="<<lastPlaneE2Hit
          <<", O1="<<lastPlaneO1Hit
          <<", O2="<<lastPlaneO2Hit<<endl;
      }

      //add up the charge of the hits in each plane
      planeSigCor[fPlane]+=fChargeSigCor;
      planeSigCorT[fPlane]+=fChargeSigCor;
      planeMip[fPlane]+=fChargeMip;
      planeMipT[fPlane]+=fChargeMip;
      if (fChargePe>planeGreatestPe[fPlane]) planeGreatestPe[fPlane]=
                                               fChargePe;
      planePe[fPlane]+=fChargePe;
      if (fO1) planeSigCorO1[fPlane]+=fChargeSigCor;
      if (fO2) planeSigCorO2[fPlane]+=fChargeSigCor;
      if (fE1) planeSigCorE1[fPlane]+=fChargeSigCor;
      if (fE2) planeSigCorE2[fPlane]+=fChargeSigCor;
      if (fO1) planeSigCorO1T[fPlane]+=fChargeSigCor;
      if (fO2) planeSigCorO2T[fPlane]+=fChargeSigCor;
      if (fE1) planeSigCorE1T[fPlane]+=fChargeSigCor;
      if (fE2) planeSigCorE2T[fPlane]+=fChargeSigCor;
    }
    //end of loop over the tracked hits

    //make sure that the first plane makes sense and was actually set
    if (fFirstPlane<0 || fFirstPlane>59) continue;

    //calc the last plane on the track with no xtalk
    //Int_t lastPlaneNoXTalk=this->CalcLastPlaneOnTrkNoXTalk();

    //make event length cuts
    if (this->CutOnEventLength()) continue;
    //make track quality cuts
    if (this->CutOnTrackQuality()) continue;

    //section for response
    const Int_t os=fLastPlane;

    //must fill the profs before looping over the maps because
    //zeros will get added to the maps at the window stage!
    //fill profs for all stripends
    this->FillProfHisto(pSigCorVsPlane,num,planeSigCor);
    this->FillProfHisto(pSigCorVsPlaneX,numX,planeSigCorX);
    this->FillProfHisto(pSigCorVsPlaneT,numT,planeSigCorT);
    //offset
    this->FillProfHisto(pSigCorVsOsPlane,numOs,planeSigCor,os);
    this->FillProfHisto(pSigCorVsOsPlaneX,numOsX,planeSigCorX,os);
    this->FillProfHisto(pSigCorVsOsPlaneT,numOsT,planeSigCorT,os);

    //fill non-offset profs for individual stripends
    this->FillProfHisto(pSigCorVsPlaneO1,numO1,planeSigCorO1);
    this->FillProfHisto(pSigCorVsPlaneO1X,numO1X,planeSigCorO1X);
    this->FillProfHisto(pSigCorVsPlaneO1T,numO1T,planeSigCorO1T);
    this->FillProfHisto(pSigCorVsPlaneO2,numO2,planeSigCorO2);
    this->FillProfHisto(pSigCorVsPlaneO2X,numO2X,planeSigCorO2X);
    this->FillProfHisto(pSigCorVsPlaneO2T,numO2T,planeSigCorO2T);
    this->FillProfHisto(pSigCorVsPlaneE1,numE1,planeSigCorE1);
    this->FillProfHisto(pSigCorVsPlaneE1X,numE1X,planeSigCorE1X);
    this->FillProfHisto(pSigCorVsPlaneE1T,numE1T,planeSigCorE1T);
    this->FillProfHisto(pSigCorVsPlaneE2,numE2,planeSigCorE2);
    this->FillProfHisto(pSigCorVsPlaneE2X,numE2X,planeSigCorE2X);
    this->FillProfHisto(pSigCorVsPlaneE2T,numE2T,planeSigCorE2T);

    //fill offset profs for individual stripends
    this->FillProfHisto(pSigCorVsOsPlaneO1,numOsO1,planeSigCorO1,os);
    this->FillProfHisto(pSigCorVsOsPlaneO1X,numOsO1X,planeSigCorO1X,os);
    this->FillProfHisto(pSigCorVsOsPlaneO1T,numOsO1T,planeSigCorO1T,os);
    this->FillProfHisto(pSigCorVsOsPlaneO2,numOsO2,planeSigCorO2,os);
    this->FillProfHisto(pSigCorVsOsPlaneO2X,numOsO2X,planeSigCorO2X,os);
    this->FillProfHisto(pSigCorVsOsPlaneO2T,numOsO2T,planeSigCorO2T,os);
    this->FillProfHisto(pSigCorVsOsPlaneE1,numOsE1,planeSigCorE1,os);
    this->FillProfHisto(pSigCorVsOsPlaneE1X,numOsE1X,planeSigCorE1X,os);
    this->FillProfHisto(pSigCorVsOsPlaneE1T,numOsE1T,planeSigCorE1T,os);
    this->FillProfHisto(pSigCorVsOsPlaneE2,numOsE2,planeSigCorE2,os);
    this->FillProfHisto(pSigCorVsOsPlaneE2X,numOsE2X,planeSigCorE2X,os);
    this->FillProfHisto(pSigCorVsOsPlaneE2T,numOsE2T,planeSigCorE2T,os);

    //section for total event sigcor and mips
    Float_t localSumSigCorTotal=0;
    Float_t localSumSigCorTotalCounter=0;
    for (map<Int_t,Float_t>::iterator sig=planeSigCorT.begin();
         sig!=planeSigCorT.end();++sig){
      localSumSigCorTotal+=sig->second;
      localSumSigCorTotalCounter++;
    }

    //add total of this event to the sum of all events 
    if (localSumSigCorTotalCounter!=0){
      sumSigCorTotal+=localSumSigCorTotal;
      sumSigCorTotalCounter++;
    }

    Float_t localSumMipTotal=0;
    Float_t localSumMipTotalCounter=0;
    for (map<Int_t,Float_t>::iterator sig=planeMipT.begin();
         sig!=planeMipT.end();++sig){
      localSumMipTotal+=sig->second;
      localSumMipTotalCounter++;
    }

    //add total of this event to the sum of all events 
    if (localSumMipTotalCounter!=0){
      sumMipTotal+=localSumMipTotal;
      sumMipTotalCounter++;
    }

    //section for sliding window
    Int_t minEventLength=35;
    if (fLastPlane<minEventLength) continue;

    //initialise the profile histos
    static Bool_t firstTime=true;
    if (firstTime){
      for (Int_t p=15;p>1;p--){
        vWindowSize.push_back(p);
        bigAv.push_back(0);
        bigAvCounter.push_back(0);
      }

      for (vector<Int_t>::iterator windowSize=vWindowSize.begin();
           windowSize!=vWindowSize.end();++windowSize){
        
        string sName=Form("%d",*windowSize);
        pSigCorVsDist.push_back
          (new TProfile(("pSigCorVsDist"+sName).c_str(),
                        ("pSigCorVsDist"+sName).c_str(),100,-2,40));
      }
    }
    firstTime=false;

    vector<TProfile*>::iterator prof=pSigCorVsDist.begin();
    vector<Float_t>::iterator av=bigAv.begin();
    vector<Int_t>::iterator count=bigAvCounter.begin();
    Float_t localSum15_18=0;
    Float_t localSum15_18Counter=0;
    Float_t localSum15_18Mip=0;
    Float_t localSum15_18MipCounter=0;

    Float_t localSum15_18_O1=0;
    Float_t localSum15_18CounterO1=0;
    Float_t localSum15_18_O2=0;
    Float_t localSum15_18CounterO2=0;
    Float_t localSum15_18_E1=0;
    Float_t localSum15_18CounterE1=0;
    Float_t localSum15_18_E2=0;
    Float_t localSum15_18CounterE2=0;
    Bool_t oneMoreEvenPlane=false;

    //loop over the different window sizes
    for (vector<Int_t>::iterator windowSize=vWindowSize.begin();
         windowSize!=vWindowSize.end();++windowSize){
      
      Float_t localSumO1=0;
      Float_t localSumCounterO1=0;
      Float_t localSumO2=0;
      Float_t localSumCounterO2=0;
      Float_t localSumE1=0;
      Float_t localSumCounterE1=0;
      Float_t localSumE2=0;
      Float_t localSumCounterE2=0;

      //slide the window along from the end to 40-windowSize
      for (Int_t p=0;p<minEventLength-*windowSize;p++){
        Int_t windowStart=fLastPlane-p;

        MSG("CDAnalysis",Msg::kVerbose)
          <<"Window size="<<*windowSize
          <<", window start="<<windowStart
          <<endl;

        //loop through the hits in the window
        for (Int_t i=0;i<*windowSize;i++){
          Int_t plane=windowStart-i;
          
          //cut out planes that are greater than the last plane
          //in the view - ie the zeros at the end
          Bool_t doZeroCorrection=false;
          if (doZeroCorrection && 
              ((plane%2==0 && plane>fLastPlaneEven) || 
               (plane%2!=0 && plane>fLastPlaneOdd))){
            MSG("CDAnalysis",Msg::kVerbose)
              <<"Not using: Event="<<event
              <<", pl="<<plane
              <<", sigCor="<<planeSigCorT[plane]
              <<", fLastPlane="<<fLastPlane
              <<", fLastPlSig="<<planeSigCorT[fLastPlane]
              <<", fLastPlOdd="<<fLastPlaneOdd
              <<", fLastPlEven="<<fLastPlaneEven
              <<endl;
          }
          else{

            //calculate the best window value
            //if fLP=48 then the first plane included is 30
            //this is actually the 19th from the end of the track
            if (*windowSize==15 && windowStart==fLastPlane-18){
              localSum15_18+=planeSigCorT[plane];//total including xtalk
              localSum15_18Counter++;
              localSum15_18Mip+=planeMipT[plane];//total including xtalk
              localSum15_18MipCounter++;

              //a 15 plane window means 7 even and 8 odd or vice versa!
              if (windowStart%2==0) oneMoreEvenPlane=true;

              //look at the individual stripends
              if (plane%2==0){//even
                localSum15_18_E1+=planeSigCorE1T[plane];
                localSum15_18CounterE1++;
                localSum15_18_E2+=planeSigCorE2T[plane];
                localSum15_18CounterE2++;
              }
              else if (plane%2==1){//odd
                localSum15_18_O1+=planeSigCorO1T[plane];
                localSum15_18CounterO1++;
                localSum15_18_O2+=planeSigCorO2T[plane];
                localSum15_18CounterO2++;
              }
            }


            if (*windowSize==15){
              //look at the individual stripends
              if (plane%2==0){//even
                localSumE1+=planeSigCorE1T[plane];
                localSumCounterE1++;
                localSumE2+=planeSigCorE2T[plane];
                localSumCounterE2++;
              }
              else if (plane%2==1){//odd
                localSumO1+=planeSigCorO1T[plane];
                localSumCounterO1++;
                localSumO2+=planeSigCorO2T[plane];
                localSumCounterO2++;
              }
            }

            (*prof)->Fill(p,planeSigCorT[plane]);
            (*av)+=planeSigCorT[plane];
            (*count)++;
          }
        }

        //calculate the near/far differences
        if (*windowSize==15){

          //a 15 plane window means 7 even and 8 odd or vice versa!
          Bool_t moreEven=false;
          if (windowStart%2==0) moreEven=true;
          Float_t divEven=1;
          Float_t divOdd=1;

          if (moreEven) {divEven=8; divOdd=7;}
          else {divEven=7; divOdd=8;}

          MAXMSG("CDAnalysis",Msg::kVerbose,70)
            <<"WindowStart="<<windowStart<<", moreEven="<<moreEven
            <<", divEven="<<divEven<<", divOdd="<<divOdd
            <<endl;
    
          Float_t largestOdd=localSumO1/divOdd;
          Float_t largestEven=localSumE1/divEven;
          
          if (localSumO2/divOdd>localSumO1/divOdd){
            largestOdd=localSumO2/divOdd;
          }
          if (localSumE2/divEven>localSumE1/divEven){
            largestEven=localSumE2/divEven;
          }

          //protect fpe
          if (largestOdd>0){
            Float_t diffNearFarOdd=100*(localSumO1/divOdd-
                                        localSumO2/divOdd)/largestOdd;
            pNearFarOdd->Fill(p,diffNearFarOdd);
          }
          if (largestEven>0){
            Float_t diffNearFarEven=100*(localSumE1/divEven-
                                         localSumE2/divEven
                                         )/largestEven;
            pNearFarEven->Fill(p,diffNearFarEven);
          }      
        }
      }
      av++;
      count++;
      prof++;
    }

    if (localSum15_18Counter!=0){
      h15_18All->Fill(localSum15_18);
      h15_18->Fill(localSum15_18/15);
      sum15_18+=localSum15_18;//add this event's energy on to total
      sum15_18Counter++;
      sum15_18Mip+=localSum15_18Mip;
      sum15_18MipCounter++;
    }
    else{
      MAXMSG("CDAnalysis",Msg::kInfo,40)
        <<"Event="<<event<<", strange 15_18="<<localSum15_18
        <<", counter="<<localSum15_18Counter<<endl;
    }

    Float_t divEven=1;
    Float_t divOdd=1;
    
    if (oneMoreEvenPlane) {divEven=8; divOdd=7;}
    else {divEven=7; divOdd=8;}
    
    h15_18_O1->Fill(localSum15_18_O1/divOdd);
    sum15_18_O1+=localSum15_18_O1/divOdd;
    sum15_18CounterO1++;
    h15_18_O2->Fill(localSum15_18_O2/divOdd);
    sum15_18_O2+=localSum15_18_O2/divOdd;
    sum15_18CounterO2++;
    
    h15_18_E1->Fill(localSum15_18_E1/divEven);
    sum15_18_E1+=localSum15_18_E1/divEven;
    sum15_18CounterE1++;
    h15_18_E2->Fill(localSum15_18_E2/divEven);
    sum15_18_E2+=localSum15_18_E2/divEven;
    sum15_18CounterE2++;

  }//end of for                                       
  

  MSG("CDAnalysis",Msg::kInfo)<<"Finished main loop"<<endl;

  //turn off the stats box printing
  gStyle->SetOptStat(0);

  //all stripends
  string sTitle="Energy Deposition (SigCor, all Stripends)";
  string sXTitle="Distance from end of track (planes)";
  this->DrawResponsePlot(sTitle,pSigCorVsPlane,pSigCorVsPlaneX,
                         pSigCorVsPlaneT,num,numX,numT);
  sTitle="Energy Deposition Vs Distance (SigCor, all Stripends)";
  this->DrawResponsePlot(sTitle,pSigCorVsOsPlane,pSigCorVsOsPlaneX,
                         pSigCorVsOsPlaneT,numOs,numOsX,numOsT,sXTitle);

  //non-offset planes
  sTitle="Energy Deposition (SigCor in Odd Planes FD stripend)";
  this->DrawResponsePlot(sTitle,pSigCorVsPlaneO1,pSigCorVsPlaneO1X,
                         pSigCorVsPlaneO1T,numO1,numO1X,numO1T);

  sTitle="Energy Deposition (SigCor in Odd Planes ND stripend)";
  this->DrawResponsePlot(sTitle,pSigCorVsPlaneO2,pSigCorVsPlaneO2X,
                         pSigCorVsPlaneO2T,numO2,numO2X,numO2T);

  sTitle="Energy Deposition (SigCor in Even Planes FD stripend)";
  this->DrawResponsePlot(sTitle,pSigCorVsPlaneE1,pSigCorVsPlaneE1X,
                         pSigCorVsPlaneE1T,numE1,numE1X,numE1T);

  sTitle="Energy Deposition (SigCor in Even Planes ND stripend)";
  this->DrawResponsePlot(sTitle,pSigCorVsPlaneE2,pSigCorVsPlaneE2X,
                         pSigCorVsPlaneE2T,numE2,numE2X,numE2T);

  //offset planes
  sTitle="Energy Deposition Vs Distance (SigCor in Odd Planes FD stripend)";
  this->DrawResponsePlot(sTitle,pSigCorVsOsPlaneO1,pSigCorVsOsPlaneO1X,
                         pSigCorVsOsPlaneO1T,numOsO1,numOsO1X,numOsO1T,
                         sXTitle);

  sTitle="Energy Deposition Vs Distance (SigCor in Odd Planes ND stripend)";
  this->DrawResponsePlot(sTitle,pSigCorVsOsPlaneO2,pSigCorVsOsPlaneO2X,
                         pSigCorVsOsPlaneO2T,numOsO2,numOsO2X,numOsO2T,
                         sXTitle);

  sTitle="Energy Deposition Vs Distance (SigCor in Even Planes FD stripend)";
  this->DrawResponsePlot(sTitle,pSigCorVsOsPlaneE1,pSigCorVsOsPlaneE1X,
                         pSigCorVsOsPlaneE1T,numOsE1,numOsE1X,numOsE1T,
                         sXTitle);

  sTitle="Energy Deposition Vs Distance (SigCor in Even Planes ND stripend)";
  this->DrawResponsePlot(sTitle,pSigCorVsOsPlaneE2,pSigCorVsOsPlaneE2X,
                         pSigCorVsOsPlaneE2T,numOsE2,numOsE2X,numOsE2T,
                         sXTitle);

  //turn on the stats box printing
  gStyle->SetOptStat(111111);

  TCanvas *c15_18=new TCanvas("c15_18","15_18",0,0,1200,800);
  c15_18->SetFillColor(0);
  c15_18->Divide(2,3);
  c15_18->cd(1);
  h15_18->Draw();
  c15_18->cd(2);
  h15_18All->Draw();
  c15_18->cd(3);
  h15_18_O1->Draw();
  c15_18->cd(4);
  h15_18_O2->Draw();
  c15_18->cd(5);
  h15_18_E1->Draw();
  c15_18->cd(6);
  h15_18_E2->Draw();

  TCanvas *cNearFar=new TCanvas("cNearFar","NearFar",0,0,1200,800);
  cNearFar->SetFillColor(0);
  cNearFar->Divide(2,3);
  cNearFar->cd(1);
  pNearFarOdd->Draw();
  cNearFar->cd(2);
  pNearFarEven->Draw();

  if (sum15_18Counter!=0){
    MSG("CDAnalysis",Msg::kInfo) 
      <<endl
      <<" ** Calibration Summary **"<<endl
      <<"SigCor in sliding window 15_18:"<<endl
      <<"    SigCor dep="<<sum15_18<<" (N="<<sum15_18Counter<<")"<<endl
      <<"    Av SigCor dep="<<sum15_18/sum15_18Counter
      <<" +/- "<<100*(h15_18->GetRMS()/sqrt(h15_18->GetEntries()))/
      h15_18->GetMean()<<"%"<<endl
      <<"    Av SigCor dep per plane="<<(sum15_18/sum15_18Counter)/15
      <<" +/- "<<100*(h15_18->GetRMS()/sqrt(h15_18->GetEntries()))/
      h15_18->GetMean()<<"%"<<endl
      <<"MIPs in sliding window 15_18:"<<endl
      <<"    MIPs dep="<<sum15_18Mip
      <<" (N="<<sum15_18MipCounter<<")"<<endl
      <<"    Av MIPs dep="<<sum15_18Mip/sum15_18MipCounter
      <<endl;
    if (sumSigCorTotalCounter!=0){
      MSG("CDAnalysis",Msg::kInfo) 
        <<"SigCor of whole event:"<<endl
        <<"    Sum of SigCor dep="<<sumSigCorTotal
        <<" (N="<<sumSigCorTotalCounter<<")"<<endl
        <<"    Av sigCor dep="
        <<sumSigCorTotal/sumSigCorTotalCounter<<endl
        <<"    1 sigCor="
        <<1800/(sumSigCorTotal/sumSigCorTotalCounter)<<" MeV"
        <<endl;
    }
    if (sumMipTotalCounter!=0){
      MSG("CDAnalysis",Msg::kInfo) 
        <<"MIPs of whole event:"<<endl
        <<"    Sum of MIPs dep="<<sumMipTotal
        <<" (N="<<sumMipTotalCounter<<")"<<endl
        <<"    Av MIPs dep="
        <<sumMipTotal/sumMipTotalCounter<<endl
        <<"    1 MIP="
        <<1800/(sumMipTotal/sumMipTotalCounter)<<" MeV"<<endl
        <<"    1 GeV="
        <<(sumMipTotal/sumMipTotalCounter)/1.8<<" MIPs"
        <<endl<<endl;
    }
  }

  if (sum15_18CounterO1!=0 && sum15_18CounterE1!=0){
    MSG("CDAnalysis",Msg::kInfo) 
      <<endl
      <<"SigCor in sliding window 15_18 in stripend 1:"<<endl
      <<" Odd planes:"<<endl
      <<"    SigCor dep="<<sum15_18_O1
      <<" (N="<<sum15_18CounterO1<<")"<<endl
      <<"    Av SigCor dep="<<sum15_18_O1/sum15_18CounterO1
      <<" +/- "<<h15_18_O1->GetRMS()/sqrt(h15_18_O1->GetEntries())
      <<" ("<<100*((h15_18_O1->GetRMS()/
                   sqrt(h15_18_O1->GetEntries()))/
                   h15_18_O1->GetMean())<<"%)"
      <<endl
      <<" Even planes:"<<endl
      <<"    SigCor dep="<<sum15_18_E1
      <<" (N="<<sum15_18CounterE1<<")"<<endl
      <<"    Av SigCor dep="<<sum15_18_E1/sum15_18CounterE1
      <<" +/- "<<h15_18_E1->GetRMS()/sqrt(h15_18_E1->GetEntries())
      <<" ("<<100*((h15_18_E1->GetRMS()/
                   sqrt(h15_18_E1->GetEntries()))/
                   h15_18_E1->GetMean())<<"%)"<<endl
      <<endl;
  }
  if (sum15_18CounterO2!=0 && sum15_18CounterE2!=0){
    MSG("CDAnalysis",Msg::kInfo) 
      <<endl
      <<"SigCor in sliding window 15_18 in stripend 2:"<<endl
      <<" Odd planes:"<<endl
      <<"    SigCor dep="<<sum15_18_O2
      <<" (N="<<sum15_18CounterO2<<")"<<endl
      <<"    Av SigCor dep="<<sum15_18_O2/sum15_18CounterO2
      <<" +/- "<<h15_18_O2->GetRMS()/sqrt(h15_18_O2->GetEntries())
      <<" ("<<100*((h15_18_O2->GetRMS()/
                   sqrt(h15_18_O2->GetEntries()))/
                   h15_18_O2->GetMean())<<"%)"
      <<endl
      <<" Even planes:"<<endl
      <<"    SigCor dep="<<sum15_18_E2
      <<" (N="<<sum15_18CounterE2<<")"<<endl
      <<"    Av SigCor dep="<<sum15_18_E2/sum15_18CounterE2
      <<" +/- "<<h15_18_E2->GetRMS()/sqrt(h15_18_E2->GetEntries())
      <<" ("<<100*((h15_18_E2->GetRMS()/
                   sqrt(h15_18_E2->GetEntries()))/
                   h15_18_E2->GetMean())<<"%)"<<endl
      <<endl;
  }

  //calculate the differences between near/far
  if (sum15_18CounterO2!=0 && sum15_18CounterE2!=0 &&
      sum15_18CounterO1!=0 && sum15_18CounterE1!=0){

    Float_t largestOdd=sum15_18_O1/sum15_18CounterO1;
    Float_t largestEven=sum15_18_E1/sum15_18CounterE1;
    Float_t diffOdd=sum15_18_O1/sum15_18CounterO1-
      sum15_18_O2/sum15_18CounterO2;
    Float_t diffEven=sum15_18_E1/sum15_18CounterE1-
      sum15_18_E2/sum15_18CounterE2;

    if (sum15_18_O2/sum15_18CounterO2>sum15_18_O1/sum15_18CounterO1) {
      largestOdd=sum15_18_O2/sum15_18CounterO2;
    }
    if (sum15_18_E2/sum15_18CounterE2>sum15_18_E1/sum15_18CounterE1) {
      largestEven=sum15_18_E2/sum15_18CounterE2;
    }
    
    MSG("CDAnalysis",Msg::kInfo) 
      <<endl
      <<"Difference between stripend 1 and 2 (1-2/largest):"<<endl
      <<" Odd planes:"<<endl
      <<"    Difference="<<diffOdd<<" ("<<100*diffOdd/largestOdd<<"%)"
      <<endl
      <<" Even planes:"<<endl
      <<"    Difference="<<diffEven<<" ("
      <<100*diffEven/largestEven<<"%)"
      <<endl;
  }

  MSG("CDAnalysis",Msg::kInfo) 
    <<" ** Finished MuonNearFar method **"<<endl;

}

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

Float_t CDAnalysis::ReCalibrate(Float_t chargeAdc,Int_t plane,
                                Int_t strip,Int_t stripend) const
{
  if (fDoReCalibration && fSimFlag==SimFlag::kData){
    MAXMSG("CDAnalysis",Msg::kInfo,10)
      <<"Running ReCalibrate: fSec="<<fSec
      <<", simFlag="<<fSimFlag<<endl;
    
    //get a calibrator
    Calibrator& cal=Calibrator::Instance();
    //VldTimeStamp (const time_t &t, const Int_t nsec)
    VldTimeStamp vldts(fSec,0);
    VldContext vc(Detector::kCalDet,fSimFlag,vldts);
    cal.ReInitialise(vc);
    MAXMSG("CDAnalysis",Msg::kVerbose,100)
      <<"Calibrator reinitialised sucessfully"<<endl;
    
    //get the seid
    PlexStripEndId seid(Detector::kCalDet,plane,strip,
                        static_cast<StripEnd::EStripEnd>(stripend));
    
    MAXMSG("CDAnalysis",Msg::kVerbose,100)
      <<"Now calibrating drift..."<<endl;
    //now calibrate: first drift, then sts
    Float_t tmpCharge=chargeAdc;
    tmpCharge=cal.GetDriftCorrected(tmpCharge,seid);
    MAXMSG("CDAnalysis",Msg::kVerbose,100)
    <<"Now calibrating sts..."<<endl;
    tmpCharge=cal.GetStripToStripCorrected(tmpCharge,seid);
    
    MAXMSG("CDAnalysis",Msg::kInfo,50)
      <<"ReCalibrate: fChargeSigCor="<<fChargeSigCor
      <<", newSigCor="<<tmpCharge
      <<", temperature="<<cal.GetTemperature()<<endl;
    
    //return the newly calibrated number
    return tmpCharge;
  }
  else {
    //just return the original value of sigcor
    return fChargeSigCor;
  }
}

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

void CDAnalysis::StoppingMuonCalibration()
{
  //open the output file for the histograms
  fOutFile=this->OpenFile(fRunNumber,"StopMuCal");

  //set the recalibration flag
  fDoReCalibration=false;

  VldTimeStamp ts;
  cout<<"Time in secs right now is:"<<ts.GetSec()<<endl;
  Int_t minRunTime=4;//is:1112976714
  if (ts.GetSec()>minRunTime){
    //analyse the electrons and fill histos
    this->ElectronResponse();
  }

  MSG("CDAnalysis",Msg::kInfo) 
    <<" ** Running StoppingMuonCalibration method... **"<<endl;

  vector<TH1F*> meuHistos(60);
  vector<TH1F*> sigCorNoPLCorHistos(60);
  //vector<TH1F*> adcNoPLCorHistos(60);//no adc info stored: June 2007
  for (Int_t i=0;i<60;i++){
    string sName="hMeuPlane";
    string sNameSigCor="hSigCorNoPLCorPlane";
    //string sNameAdc="hAdcNoPLCorPlane";
    string sNum=Form("%d",i);
    sName+=sNum;
    sNameSigCor+=sNum;
    //sNameAdc+=sNum;

    meuHistos[i]=new TH1F(sName.c_str(),sName.c_str(),300,-2,3000);
    meuHistos[i]->SetBit(TH1::kCanRebin);//leave this as is June/07
    //meuHistos[i]->Print();

    //added sigcor in June 2007
    sigCorNoPLCorHistos[i]=new TH1F
      (sNameSigCor.c_str(),sNameSigCor.c_str(),300,-2,3000);
    //adcNoPLCorHistos[i]=new TH1F
    //(sName.c_str(),sName.c_str(),300,-2,3000);
  }
  
  TH1F *hMeu=new TH1F("hMeu","hMeu",1000,-2,1000);
  hMeu->GetXaxis()->SetTitle("SigCors");
  hMeu->GetXaxis()->CenterTitle();
  hMeu->GetYaxis()->SetTitle("");
  hMeu->GetYaxis()->CenterTitle();
  hMeu->SetFillColor(0);
  hMeu->SetBit(TH1::kCanRebin);

  TH1F* hGeVPerMeu=new TH1F("hGeVPerMeu","hGeVPerMeu",3000,-0.001,0.3);
  hGeVPerMeu->SetFillColor(0);
  hGeVPerMeu->SetTitle("GeV/MEU");
  hGeVPerMeu->GetXaxis()->SetTitle("GeV/MEU");
  hGeVPerMeu->GetXaxis()->CenterTitle();

  TH1F *hMeuPLCor=new TH1F("hMeuPLCor","hMeuPLCor",1000,-2,1000);
  hMeuPLCor->GetXaxis()->SetTitle("SigCors");
  hMeuPLCor->GetXaxis()->CenterTitle();
  hMeuPLCor->GetYaxis()->SetTitle("");
  hMeuPLCor->GetYaxis()->CenterTitle();
  hMeuPLCor->SetFillColor(0);
  hMeuPLCor->SetBit(TH1::kCanRebin);

  TH1F *hMeuTruePLCor=new TH1F("hMeuTruePLCor","hMeuTruePLCor",
                               1000,-2,1000);
  hMeuTruePLCor->GetXaxis()->SetTitle("SigCors");
  hMeuTruePLCor->GetXaxis()->CenterTitle();
  hMeuTruePLCor->GetYaxis()->SetTitle("");
  hMeuTruePLCor->GetYaxis()->CenterTitle();
  hMeuTruePLCor->SetFillColor(0);
  hMeuTruePLCor->SetBit(TH1::kCanRebin);

  TH1F *hMeuOnlyEventLengthCut=new TH1F("hMeuOnlyEventLengthCut",
                                        "hMeuOnlyEventLengthCut",
                                        1000,-2,1000);
  hMeuOnlyEventLengthCut->GetXaxis()->SetTitle("SigCors");
  hMeuOnlyEventLengthCut->GetXaxis()->CenterTitle();
  hMeuOnlyEventLengthCut->GetYaxis()->SetTitle("");
  hMeuOnlyEventLengthCut->GetYaxis()->CenterTitle();
  hMeuOnlyEventLengthCut->SetFillColor(0);
  hMeuOnlyEventLengthCut->SetBit(TH1::kCanRebin);

  TH1F *hMeuOnlyEL_FidCut=new TH1F("hMeuOnlyEL_FidCut",
                                   "hMeuOnlyEL_FidCut",
                                   1000,-2,1000);
  hMeuOnlyEL_FidCut->GetXaxis()->SetTitle("SigCors");
  hMeuOnlyEL_FidCut->GetXaxis()->CenterTitle();
  hMeuOnlyEL_FidCut->GetYaxis()->SetTitle("");
  hMeuOnlyEL_FidCut->GetYaxis()->CenterTitle();
  hMeuOnlyEL_FidCut->SetFillColor(0);
  hMeuOnlyEL_FidCut->SetBit(TH1::kCanRebin);

  TH1F *hMeuOnlyEL_Fid_TimeCut=new TH1F("hMeuOnlyEL_Fid_TimeCut",
                                        "hMeuOnlyEL_Fid_TimeCut",
                                        1000,-2,1000);
  hMeuOnlyEL_Fid_TimeCut->GetXaxis()->SetTitle("SigCors");
  hMeuOnlyEL_Fid_TimeCut->GetXaxis()->CenterTitle();
  hMeuOnlyEL_Fid_TimeCut->GetYaxis()->SetTitle("");
  hMeuOnlyEL_Fid_TimeCut->GetYaxis()->CenterTitle();
  hMeuOnlyEL_Fid_TimeCut->SetFillColor(0);
  hMeuOnlyEL_Fid_TimeCut->SetBit(TH1::kCanRebin);

  TH1F *hMeuOnlyEL_TimeCut=new TH1F("hMeuOnlyEL_TimeCut",
                                   "hMeuOnlyEL_TimeCut",
                                   1000,-2,1000);
  hMeuOnlyEL_TimeCut->GetXaxis()->SetTitle("SigCors");
  hMeuOnlyEL_TimeCut->GetXaxis()->CenterTitle();
  hMeuOnlyEL_TimeCut->GetYaxis()->SetTitle("");
  hMeuOnlyEL_TimeCut->GetYaxis()->CenterTitle();
  hMeuOnlyEL_TimeCut->SetFillColor(0);
  hMeuOnlyEL_TimeCut->SetBit(TH1::kCanRebin);

  TH1F *hMeuOnlyEL_HppCut=new TH1F("hMeuOnlyEL_HppCut",
                                   "hMeuOnlyEL_HppCut",
                                   1000,-2,1000);
  hMeuOnlyEL_HppCut->GetXaxis()->SetTitle("SigCors");
  hMeuOnlyEL_HppCut->GetXaxis()->CenterTitle();
  hMeuOnlyEL_HppCut->GetYaxis()->SetTitle("");
  hMeuOnlyEL_HppCut->GetYaxis()->CenterTitle();
  hMeuOnlyEL_HppCut->SetFillColor(0);
  hMeuOnlyEL_HppCut->SetBit(TH1::kCanRebin);

  TH1F *hMeuOnlyPID_EventLengthCut=new TH1F("hMeuOnlyPID_EventLengthCut",
                                        "hMeuOnlyPID_EventLengthCut",
                                        1000,-2,1000);
  hMeuOnlyPID_EventLengthCut->GetXaxis()->SetTitle("SigCors");
  hMeuOnlyPID_EventLengthCut->GetXaxis()->CenterTitle();
  hMeuOnlyPID_EventLengthCut->GetYaxis()->SetTitle("");
  hMeuOnlyPID_EventLengthCut->GetYaxis()->CenterTitle();
  hMeuOnlyPID_EventLengthCut->SetFillColor(0);
  hMeuOnlyPID_EventLengthCut->SetBit(TH1::kCanRebin);

  TH1F *hMeuOnlyPID_EL_FidCut=new TH1F("hMeuOnlyPID_EL_FidCut",
                                   "hMeuOnlyPID_EL_FidCut",
                                   1000,-2,1000);
  hMeuOnlyPID_EL_FidCut->GetXaxis()->SetTitle("SigCors");
  hMeuOnlyPID_EL_FidCut->GetXaxis()->CenterTitle();
  hMeuOnlyPID_EL_FidCut->GetYaxis()->SetTitle("");
  hMeuOnlyPID_EL_FidCut->GetYaxis()->CenterTitle();
  hMeuOnlyPID_EL_FidCut->SetFillColor(0);
  hMeuOnlyPID_EL_FidCut->SetBit(TH1::kCanRebin);

  TH1F *hMeuOnlyPID_EL_Fid_TimeCut=new TH1F("hMeuOnlyPID_EL_Fid_TimeCut",
                                        "hMeuOnlyPID_EL_Fid_TimeCut",
                                        1000,-2,1000);
  hMeuOnlyPID_EL_Fid_TimeCut->GetXaxis()->SetTitle("SigCors");
  hMeuOnlyPID_EL_Fid_TimeCut->GetXaxis()->CenterTitle();
  hMeuOnlyPID_EL_Fid_TimeCut->GetYaxis()->SetTitle("");
  hMeuOnlyPID_EL_Fid_TimeCut->GetYaxis()->CenterTitle();
  hMeuOnlyPID_EL_Fid_TimeCut->SetFillColor(0);
  hMeuOnlyPID_EL_Fid_TimeCut->SetBit(TH1::kCanRebin);

  TH1F *hMeuOnlyPID_EL_TimeCut=new TH1F("hMeuOnlyPID_EL_TimeCut",
                                   "hMeuOnlyPID_EL_TimeCut",
                                   1000,-2,1000);
  hMeuOnlyPID_EL_TimeCut->GetXaxis()->SetTitle("SigCors");
  hMeuOnlyPID_EL_TimeCut->GetXaxis()->CenterTitle();
  hMeuOnlyPID_EL_TimeCut->GetYaxis()->SetTitle("");
  hMeuOnlyPID_EL_TimeCut->GetYaxis()->CenterTitle();
  hMeuOnlyPID_EL_TimeCut->SetFillColor(0);
  hMeuOnlyPID_EL_TimeCut->SetBit(TH1::kCanRebin);

  TH1F *hMeuOnlyPID_EL_HppCut=new TH1F("hMeuOnlyPID_EL_HppCut",
                                   "hMeuOnlyPID_EL_HppCut",
                                   1000,-2,1000);
  hMeuOnlyPID_EL_HppCut->GetXaxis()->SetTitle("SigCors");
  hMeuOnlyPID_EL_HppCut->GetXaxis()->CenterTitle();
  hMeuOnlyPID_EL_HppCut->GetYaxis()->SetTitle("");
  hMeuOnlyPID_EL_HppCut->GetYaxis()->CenterTitle();
  hMeuOnlyPID_EL_HppCut->SetFillColor(0);
  hMeuOnlyPID_EL_HppCut->SetBit(TH1::kCanRebin);

  TH1F *hMeuNoPidCut=new TH1F("hMeuNoPidCut","hMeuNoPidCut",
                              1000,-2,1000);
  hMeuNoPidCut->GetXaxis()->SetTitle("SigCors");
  hMeuNoPidCut->GetXaxis()->CenterTitle();
  hMeuNoPidCut->GetYaxis()->SetTitle("");
  hMeuNoPidCut->GetYaxis()->CenterTitle();
  hMeuNoPidCut->SetFillColor(0);
  hMeuNoPidCut->SetBit(TH1::kCanRebin);

  TH1F *hMeuNoEventLengthCut=new TH1F("hMeuNoEventLengthCut",
                                      "hMeuNoEventLengthCut",
                                      1000,-2,1000);
  hMeuNoEventLengthCut->GetXaxis()->SetTitle("SigCors");
  hMeuNoEventLengthCut->GetXaxis()->CenterTitle();
  hMeuNoEventLengthCut->GetYaxis()->SetTitle("");
  hMeuNoEventLengthCut->GetYaxis()->CenterTitle();
  hMeuNoEventLengthCut->SetFillColor(0);
  hMeuNoEventLengthCut->SetBit(TH1::kCanRebin);

  TH1F *hMeuNoFidVolCut=new TH1F("hMeuNoFidVolCut","hMeuNoFidVolCut",
                                 1000,-2,1000);
  hMeuNoFidVolCut->GetXaxis()->SetTitle(