MeuAnalysis.cxx

Go to the documentation of this file.

// Coded by Jeff Hartnell Jul/2005 onwards
//                                                                    
// Purpose: To provide a framework for analysis of stopping muon data 
//                                                                    
// Contact: j.j.hartnell@rl.ac.uk                         

#include <set>

#include "TChain.h"
#include "TClonesArray.h"
#include "TF1.h"
#include "TFile.h"
#include "TGraph.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TObjString.h"
#include "TProfile.h"
#include "TRegexp.h"
#include "TStyle.h"
#include "TSystem.h"

#include "AtNuEvent/AtmosEvent.h"
#include "AtNuEvent/AtmosStrip.h"
#include "AtNuEvent/AtmosTrack.h"
#include "BeamDataNtuple/BeamDataLiteHeader.h"
#include "BeamDataUtil/BeamMonSpill.h"
#include "BeamDataUtil/BMSpillAna.h"
#include "Calibrator/Calibrator.h"
#include "MessageService/MsgService.h"
#include "BeamDataNtuple/NtpBDLiteRecord.h"
#include "CandNtupleSR/NtpSREvent.h"
#include "CandNtupleSR/NtpSRTrack.h"
#include "CandNtupleSR/NtpSRShower.h"
#include "CandNtupleSR/NtpSRSlice.h"
#include "CandNtupleSR/NtpSRStrip.h"
#include "StandardNtuple/NtpStRecord.h"
#include "Plex/PlexHandle.h"
#include "Record/RecCandHeader.h"
#include "UgliGeometry/UgliGeomHandle.h"

#include "MeuCal/MeuAnalysis.h"
#include "MeuCal/MeuCounter.h"
#include "MeuCal/MeuCuts.h"
#include "MeuCal/MeuHitInfo.h"
#include "MeuCal/MeuHistos.h"
#include "MeuCal/MeuReco.h"
#include "MeuCal/MeuSummary.h"
#include "MeuCal/MeuSummaryWriter.h"

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

//initialise the static variables
string MeuAnalysis::fInputFileName="";//use hard coded place by default
Bool_t MeuAnalysis::fLoadTrees=false;//don't load trees by default
Bool_t MeuAnalysis::fRecalibrateSigLin=false;//no recalibrate by default
Bool_t MeuAnalysis::fRecalibrateSigCor=false;//no recalibrate by default
Bool_t MeuAnalysis::fRecalibratePE=false;//no recalibrate by default

CVSID("$Id: MeuAnalysis.cxx,v 1.33 2007/06/06 05:07:12 rhatcher Exp $");

ClassImp(MeuAnalysis)
  
//......................................................................

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

  fChain=0;
  fChainBD=0;
  fev=0;//atnu one
  fOutFile=0;
  fRec=0;//sr
  fRecBD=0;//sr beam data
  fUseAtNu=false;
  
  fEntries=0;
  fEntriesBD=0;
  fRun=100;//default
  fSubrun=0;

  if (fLoadTrees) {
    MSG("MeuAnalysis",Msg::kInfo)<<"Loading trees"<<endl;
    this->MakeChain();
    this->SetRootFileObjects();
    //this->InitialiseHitInfoVariables();
  
    //get basic run number stuff
    if (fRec && fChain){//sr
      //get the first event
      fChain->GetEvent(0); 
      NtpStRecord& ntp=(*fRec);
      const RecCandHeader& rec=ntp.GetHeader();
      MSG("MeuAnalysis",Msg::kInfo)
        <<"Found: run="<<rec.GetRun()
        <<", subrun="<<rec.GetSubRun()<<endl;
      fRun=rec.GetRun();
      fSubrun=rec.GetSubRun();
    }
    else if (fev && fChain){//atnu
      //get the first event
      fChain->GetEvent(0); 
      MSG("MeuAnalysis",Msg::kInfo)
        <<"Found: run="<<fev->Run
        <<", subrun="<<fev->SubRun<<endl;
      fRun=fev->Run;
      fSubrun=fev->SubRun;
    }
    else cout<<"Ahhh, no root files"<<endl;
  
  }
  else  {
    MSG("MeuAnalysis",Msg::kInfo)
      <<"Not loading trees in constructor"<<endl;
  }
  
  //setup some root stuff
  //include the under and overflow counts
  gStyle->SetOptStat(1111111);
  gStyle->SetOptFit(1111);
  gStyle->SetPalette(1);
  
  MSG("MeuAnalysis",Msg::kInfo)
    <<"Finished MeuAnalysis Constructor"<<endl;
}

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

MeuAnalysis::~MeuAnalysis()
{
  MSG("MeuAnalysis",Msg::kInfo)
    <<"Running MeuAnalysis Destructor..."<<endl;
  
  if (fOutFile){
    //this makes histos disappear from the canvases 
    //so do it at the very end
    fOutFile->Close();
  }
  
  MSG("MeuAnalysis",Msg::kDebug)
    <<"Finished MeuAnalysis Destructor"<<endl;
}

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

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

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

void MeuAnalysis::LoadTrees(Bool_t loadTrees)
{
  MSG("MeuAnalysis",Msg::kInfo)
    <<"Setting fLoadTrees="<<loadTrees<<endl;
  fLoadTrees=loadTrees;
}

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

void MeuAnalysis::RecalibrateSigLin(Bool_t recalibrate)
{
  MSG("MeuAnalysis",Msg::kInfo)
    <<"Setting fRecalibrateSigLin="<<recalibrate<<endl;
  fRecalibrateSigLin=recalibrate;
  if (recalibrate) {
    MSG("MeuAnalysis",Msg::kInfo)
      <<"Note that the calibrator has to be correctly configured"
      <<" for the recalibration to work!"<<endl;
  }
}

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

void MeuAnalysis::RecalibrateSigCor(Bool_t recalibrate)
{
  MSG("MeuAnalysis",Msg::kInfo)
    <<"Setting fRecalibrateSigCor="<<recalibrate<<endl;
  fRecalibrateSigCor=recalibrate;
  if (recalibrate) {
    MSG("MeuAnalysis",Msg::kInfo)
      <<"Note that the calibrator has to be correctly configured"
      <<" for the recalibration to work!"<<endl;
  }
}

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

void MeuAnalysis::RecalibratePE(Bool_t recalibrate)
{
  MSG("MeuAnalysis",Msg::kInfo)
    <<"Setting fRecalibratePE="<<recalibrate<<endl;
  fRecalibratePE=recalibrate;
  if (recalibrate) {
    MSG("MeuAnalysis",Msg::kInfo)
      <<"Note that the calibrator has to be correctly configured"
      <<" for the recalibration to work!"<<endl;
  }
}

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

void MeuAnalysis::UseAtNu()
{
  //fUseAtNu=true;
}

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

void MeuAnalysis::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=0;
      //runInfo.Branch("SimFlag",&lSimFlag,"SimFlag/I");


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

      MSG("MeuAnalysis",Msg::kInfo)
        <<"Writing histos to: "<<fOutFile->GetName()<<" ..."<<endl;
      fOutFile->Write();
      //fOutFile->Close();//this makes histos disappear from canvases
      //so do it in the destructor (need to make MeuAnalysis on heap)
    }
    else {
      MSG("MeuAnalysis",Msg::kWarning)
        <<"File not writable!"<<endl;
    }
  }
}

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

vector<string> MeuAnalysis::MakeFileList()
{

  vector<string> fileList;
  
  if (fInputFileName!=""){

    Int_t findGives=fInputFileName.find(".root");
    MSG("MeuAnalysis",Msg::kInfo)
      <<"Checking input string for \"*.root\", position in string="
      <<findGives<<endl;

    if (findGives>0){
      //add the file direct to the list
      fileList.push_back(fInputFileName);
    }
    else{//treat the file as a list of root files 
      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("MeuAnalysis",Msg::kDebug)
            <<"Found input file name="<<file<<endl;
          fileList.push_back(file);
        }
        MSG("MeuAnalysis",Msg::kDebug)
          <<"Files names found in txt file="<<fileList.size()<<endl;
      }
      else{
        MSG("MeuAnalysis",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("MEUDATA");
  if (envVariable==NULL){
    MSG("MeuAnalysis",Msg::kFatal)
      <<endl<<endl
      <<"*************************************************************"
      <<endl<<"Environmental variable MEUDATA not set!"<<endl
      <<"Please set MEUDATA to the directory containing the"
      <<" ntuple 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="";
  sEnv=envVariable;
  sEnv+="/*.root";
  MSG("MeuAnalysis",Msg::kInfo)
    <<"Looking for *.root files using the env variable"<<endl
    <<"MEUDATA="<<sEnv<<endl;    
  fileList.push_back(sEnv);

  return fileList;
}


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

Bool_t MeuAnalysis::ObjectExistsInFile(TFile* file,
                                        std::string objectName) const
{
  if (!file) {
    MSG("MeuAnalysis",Msg::kWarning)
      <<"File does not exist so can't test if object exists!"<<endl;
    return false;
  }

  TList* listOfKeys=file->GetListOfKeys();

  static TFile* lastFile=0;
  if (lastFile!=file){
    MSG("MeuAnalysis",Msg::kInfo)
      <<"File contains these keys:"<<endl;
    listOfKeys->Print();
  }
  lastFile=file;

  TObject* ob=listOfKeys->FindObject(objectName.c_str());
  if (ob) {
    MSG("MeuAnalysis",Msg::kInfo)
      <<"Requested object with name="<<objectName
      <<" exists in file:"<<endl;
    ob->Print();
    return true;
  }
  else {
    MSG("MeuAnalysis",Msg::kInfo)
      <<"Requested object with name="<<objectName
      <<" does NOT exist in file"<<endl;
    return false;
  }
}

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

std::string MeuAnalysis::GetFirstFileName(std::string wildcardString) const
{
  //check if there is actually a wildcard
  if (!TString(wildcardString.c_str()).MaybeWildcard()) {
    return wildcardString;
  }
  
  //this was ripped off from ROOT's TChain::AddFile

  vector<string> fileList;

  //wildcarding used in name
  TString basename(wildcardString.c_str());

  Int_t dotslashpos = basename.Index(".root/");
  TString behind_dot_root;
  if (dotslashpos>=0) {
    // Copy the tree name specification
    behind_dot_root = basename(dotslashpos+6,basename.Length()-dotslashpos+6);
    // and remove it from basename
    basename.Remove(dotslashpos+5);
  }

  Int_t slashpos = basename.Last('/');
  TString directory;
  if (slashpos>=0) {
    directory = basename(0,slashpos); // Copy the directory name
    basename.Remove(0,slashpos+1);      // and remove it from basename
  } else {
    directory = gSystem->WorkingDirectory();
  }

  const char *file;
  void *dir = gSystem->OpenDirectory(gSystem->ExpandPathName(directory.Data()));

  if (dir) {
    //create a TList to store the file names (not yet sorted)
    TList l;
    TRegexp re(basename,kTRUE);
    while ((file = gSystem->GetDirEntry(dir))) {
      if (!strcmp(file,".") || !strcmp(file,"..")) continue;
      TString s = file;
      if ( (basename!=file) && s.Index(re) == kNPOS) continue;
      l.Add(new TObjString(file));
    }
    gSystem->FreeDirectory(dir);
    //sort the files in alphanumeric order
    l.Sort();
    TIter next(&l);
    TObjString *obj;
    while ((obj = (TObjString*)next())) {
      file = obj->GetName();
      if (behind_dot_root.Length()!=0){
        string fileName=Form("%s/%s/%s",directory.Data(),
                             file,behind_dot_root.Data());
        fileList.push_back(fileName);
        //nf += AddFile(Form("%s/%s/%s",directory.Data(),file,behind_dot_root.Data()),kBigNumber);
      }
      else {
        string fileName=Form("%s/%s",directory.Data(),file);
        fileList.push_back(fileName);
        //nf += AddFile(Form("%s/%s",directory.Data(),file),kBigNumber);
      }
    }
    l.Delete();
  }

  //check if any files were found
  if (fileList.begin()!=fileList.end()){
    cout<<"Used wildcard expansion to find first file name="
        <<*fileList.begin()<<endl;
    return *fileList.begin();
  }
  else {
    return "";//null string
  }
}

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

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

  TFile* firstFile=0;
  // case with one single file
  if (!TString(*fileList.begin()).MaybeWildcard()) {
    firstFile=TFile::Open((*fileList.begin()).c_str(),"READ");
  }
  else {
    MSG("MeuAnalysis",Msg::kInfo)
      <<"Found a wildcard present in list of file names: "
      <<*fileList.begin()
      <<endl<<"Looking for first filename using wildcard..."<<endl;
    string fileName=this->GetFirstFileName(*fileList.begin());
    if (fileName!="") {
      firstFile=TFile::Open(fileName.c_str(),"READ");
    }
    else {
      MSG("MeuAnalysis",Msg::kInfo)
        <<"No files found matching:"<<endl
        <<*fileList.begin()<<", will exit here..."<<endl;
      exit(1);
    }
  }
  
  //create a chain
  if (false && fUseAtNu) {
    if (this->ObjectExistsInFile(firstFile,"ntp")) {
      fChain=new TChain("ntp");
    }
  }
  else {
    //check if the AtmosEvent ntuple exists in the file
    if (this->ObjectExistsInFile(firstFile,"ntp")) {
      fChain=new TChain("ntp");
      fUseAtNu=true;//set this to true
    }
    //else check if the SR ntuple exists
    else if (this->ObjectExistsInFile(firstFile,"NtpSt")) {
      fChain=new TChain("NtpSt");
    }
    else {
      MSG("MeuAnalysis",Msg::kInfo)
        <<"Not creating NtpSt TChain because it does not exist in file"
        <<endl;
    }

    //check in BD ntuple exists
    if (this->ObjectExistsInFile(firstFile,"NtpBDLite")) {
      fChainBD=new TChain("NtpBDLite");//have to load libraries too
    }
    else {
      MSG("MeuAnalysis",Msg::kInfo)
        <<"Not creating NtpBDLite TChain because it does not exist"
        <<" in file"<<endl;
    }
  }
  
  firstFile->Close();
  if (firstFile) delete firstFile;
  firstFile=0;

  Int_t nf=0;
  Int_t nfBD=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("MeuAnalysis",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("MeuAnalysis",Msg::kInfo)<<"Adding file(s)="<<*file<<endl;
    nf+=fChain->Add((*file).c_str());
    if (fChainBD) nfBD+=fChainBD->Add((*file).c_str());
    //fChain->Print();//very verbose
  }

  if(nf==0){
    MSG("MeuAnalysis",Msg::kFatal)
      <<endl<<endl
      <<"*************************************************************"
      <<endl<<"No *.root files found"<<endl
      <<"Please set MEUDATA to the directory containing the"
      <<" *.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("MeuAnalysis",Msg::kInfo)
    <<endl<<"Added "<<nf<<" file(s) to TChain"<<endl;
  MSG("MeuAnalysis",Msg::kInfo)
    <<endl<<"Added "<<nfBD<<" beam data (BD) file(s) to TChain"
    <<endl;
  
  if (fChain) {
    MSG("MeuAnalysis",Msg::kInfo) 
      <<"Ntuple information:"<<endl;
    fChain->Show(0);
    
    if (fChainBD) {
      MSG("MeuAnalysis",Msg::kInfo) 
        <<"NtpBDLite information:"<<endl;
      fChainBD->Show(0);
    }
  }
  
  MSG("MeuAnalysis",Msg::kInfo)
    <<endl<<"Added "<<nfBD<<" beam data (BD) file(s) to TChain."
    <<endl;
  MSG("MeuAnalysis",Msg::kInfo)
    <<endl<<"Analysing "<<nf<<" file(s). Reading from disk..."<<endl;
}

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

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

  if (fUseAtNu) {
    fChain->SetBranchAddress("evt", &fev);
    
    //set the number of events in the tree
    fEntries=static_cast<Int_t>(fChain->GetEntries());
    MSG("MeuAnalysis",Msg::kInfo)
      <<"AtNuTree has "<<fEntries<<" entries"<<endl;
    if (fEntries>0){
      //get the snarl/event
      fChain->GetEvent(0); 
      //fRunNumber=(*fev).Run;
      MSG("MeuAnalysis",Msg::kInfo)
        <<"SetRootFileObjects: first UnixTime="<<(*fev).UnixTime<<endl;
    }
  }
  else {
    //set up tree
    fChain->SetBranchAddress("NtpStRecord",&fRec);
    if (fChainBD) fChainBD->SetBranchAddress("NtpBDLiteRecord",&fRecBD);
    //fBDtortgt
    
    //set the number of events in the tree
    fEntries=static_cast<Int_t>(fChain->GetEntries());
    if (fChainBD) fEntriesBD=static_cast<Int_t>(fChainBD->GetEntries());
    MSG("MeuAnalysis",Msg::kInfo)
      <<"NtpSt tree has "<<fEntries<<" entries"<<endl;
    MSG("MeuAnalysis",Msg::kInfo)
      <<"NtpBDLite tree has "<<fEntriesBD<<" entries"<<endl;
    
    if (fEntries>0){
      //get the snarl/event
      fChain->GetEvent(0); 
      //MSG("MeuAnalysis",Msg::kInfo)
      //<<"First run number="<<fRec->fHeader.fRun<<endl;
      MSG("MeuAnalysis",Msg::kInfo)
        <<"SetRootFileObjects: first snarl nslice="
        <<fRec->evthdr.nslice<<endl;
    }

    if (fEntriesBD>0){
      //get the snarl/event
      fChainBD->GetEvent(0); 
      MSG("MeuAnalysis",Msg::kInfo)
        <<"First tortgt="<<fRecBD->tortgt<<" E12 POT"<<endl;
    }
  }
  MSG("MeuAnalysis",Msg::kDebug)
    <<"Finished the SetRootFileObjects method"<<endl;
}

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

TFile* MeuAnalysis::OpenFile(Int_t runNumber,string prefix) const
{
  //create the tfile pointer to be returned
  TFile* outputFile=0;
  
  //get the environmental variable
  char *anaDir=getenv("MEUANA_DIR");
  
  //use a string to hold env instead 
  string sAnaDir="";
  
  if (anaDir!=NULL) {
    sAnaDir=anaDir;
  }
  else {
    MSG("MeuAnalysis",Msg::kInfo) 
      <<"Environmental variable $MEUANA_DIR not set." 
      <<" Writing file(s) to current directory"<<endl;
    sAnaDir=".";
  }
  
  //convert variables to string
  string sRunNumber=Form("%d",runNumber);
  //string sDetector="C";
  string sDetector="";
  //string sPrefix="h";//default
  string sPrefix="";//default
  if (prefix!="") sPrefix+=prefix;
  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("MeuAnalysis",Msg::kInfo) 
    <<"Output file opened: "<<sTmp<<endl;
  return outputFile;
}

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

ofstream* MeuAnalysis::OpenTxtFile(Int_t runNumber,
                                    std::string prefix) const
{
  //create the tfile pointer to be returned
  ofstream* outputFile=0;
  
  //convert variables to string
  string sRunNumber=Form("%d",runNumber);
  string sDetector="F";
  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("MeuAnalysis",Msg::kInfo) 
      <<"Output txt file opened: "<<sFileName<<endl;
  }
  else{
    MSG("MeuAnalysis",Msg::kInfo)
      <<"Txt file NOT opened: "<<sFileName<<endl;
  }

  return outputFile;
}

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

void MeuAnalysis::SetLoopVariables(Int_t entry,Int_t printMode)
{
  if (printMode==1){
    Float_t fract=ceil(fEntries/20.);  
    if (ceil(((Float_t)entry)/fract)==((Float_t)entry)/fract){
      MSG("MeuAnalysis",Msg::kInfo) 
        <<"Fraction of loop complete: "<<entry 
        <<"/"<<fEntries<<"  ("
        <<(Int_t)(100.*entry/fEntries)<<"%)"<<endl;
    }
  }

  //get the snarl/entry
  fChain->GetEvent(entry);
  if (fChainBD) fChainBD->GetEvent(entry); 
}

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

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

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

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

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

  
  //this->InitialiseLoopVariables();  
  
  //  for(Int_t event=0;event<fEntries;event++){
    
  //this->SetLoopVariables(event);

  for (int i=0;i<5;i++) {
    fChain->GetEntry(i);
    NtpStRecord& ntpst=(*fRec);   // Make a reference instead of a pointer
    
    TClonesArray& tcaTk=(*ntpst.trk);
    Int_t numTrks=tcaTk.GetEntries();
   
    MSG("MeuAnalysis",Msg::kInfo)
      <<"numTrks="<<numTrks
      <<", ndigit="<<fRec->evthdr.ndigit
      <<", nTrack="<<fRec->evthdr.ntrack<<endl;

    //Loop over tracks
    for (Int_t itrk=0;itrk<numTrks;itrk++){
      const NtpSRTrack* ptrk=
        dynamic_cast<NtpSRTrack*>(tcaTk[itrk]);
      const NtpSRTrack& trk=(*ptrk);
      cout<<"range="<<trk.range<<" nstrip="<<trk.nstrip<<" ds="<<trk.ds
          <<" fidend.dr,dz="<<trk.fidend.dr<<" "<<trk.fidend.dz<<" end-vtx.z:"
          <<trk.end.z-trk.vtx.z<<" "<<trk.stpz[1]-trk.stpz[0]<<endl;


      TClonesArray& tcaStp=(*ntpst.stp);
      Int_t numStps=tcaStp.GetEntries();
      MSG("MeuAnalysis",Msg::kInfo)
        <<"numStps="<<numStps<<endl;
      
      for (Int_t i=0;i<trk.nstrip;i++){

        //check for bug where strip index is -1
        if (trk.stp[i]<0) {
          MAXMSG("MeuAnalysis",Msg::kInfo,500)
            <<"Skipping strip with trk.stp[i]="<<trk.stp[i]<<endl;
          continue;
        }

        //const NtpSRStrip* pstp=dynamic_cast<NtpSRStrip*>(tcaStp[i]);
        const NtpSRStrip* pstp=dynamic_cast<NtpSRStrip*>(tcaStp[trk.stp[i]]);

        const NtpSRStrip& stp=(*pstp);
        MSG("MeuAnalysis",Msg::kInfo)
          <<"Strip="<<stp.strip<<", pl="<<stp.plane
          <<", sigCor="<<stp.ph1.sigcor<<", tpos="<<stp.tpos<<endl;
      }

    }


    TClonesArray& tcaStp=(*ntpst.stp);
    Int_t numStps=tcaStp.GetEntries();
    MSG("MeuAnalysis",Msg::kInfo)
      <<endl<<"2nd: numStps="<<numStps<<endl;
    
    for (Int_t i=0;i<numStps;i++){      
      const NtpSRStrip* pstp=dynamic_cast<NtpSRStrip*>(tcaStp[i]);
      const NtpSRStrip& stp=(*pstp);
      MSG("MeuAnalysis",Msg::kInfo)
        <<"Strip="<<stp.strip<<", pl="<<stp.plane<<endl;
    }
    
    //    for (Int_t nsloop=0;nsloop<trk.nstrip;nsloop++) {
    //const NtpSRStrip* pstp=
    //dynamic_cast<NtpSRStrip*>(tcaStp[trk.stp[ns]]);
    //const NtpSRStrip& stp=(*pstp);

    //      if (stp.plane) ;   // Avoid compiler message ifndef PRINTLOTS
    //float du = trk.stpu[ns]-uprev;
    //float dv = trk.stpv[ns]-vprev;
    //float dz = trk.stpz[ns]-zprev;
    //cout<<" "<<ns<<" "<<du<<" "<<dv<<" "<<dz
    //    <<" uv:"<<trk.stpu[ns] <<" "<<trk.stpv[ns]
    //    <<" xyz:"<<trk.stpx[ns] <<" "<<trk.stpy[ns]<<" "<<trk.stpz[ns]
    //    <<" ds:"<<trk.stpds[ns]<<" stp: "<<trk.stp[ns]
    //    <<" "<<stp.plane<<endl;
    //  printf("%2d%8.4f%8.4f%8.4f%8.4f%2d uv:%8.4f%8.4f"
    //       " xyz:%8.4f%8.4f%8.4f%8.4f"
    //       " stp: %4d%4d\n"
    //       ,ns,du,dv,dz,stp.tpos,stp.planeview
    //       ,trk.stpu[ns],trk.stpv[ns]
    //       ,trk.stpx[ns],trk.stpy[ns],trk.stpz[ns]
    //       ,trk.stpds[ns],trk.stp[ns],stp.plane);
      
    //uprev = trk.stpu[ns];
    //vprev = trk.stpv[ns];
    //zprev = trk.stpz[ns];
    //}// stp loop

  }//end of for                                       
  

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

  MSG("MeuAnalysis",Msg::kInfo) 
    <<" ** Finished Test method **"<<endl;
}

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

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

  MeuCuts cuts;
  
  
  this->InitialiseLoopVariables();  
  
  for(Int_t entry=0;entry<fEntries;entry++) {
  //for(Int_t entry=0;entry<5000;entry++) {
    
    this->SetLoopVariables(entry);

    NtpStRecord& ntp=(*fRec);

    TClonesArray& evtTca=(*ntp.evt);
    const Int_t numEvts=evtTca.GetEntriesFast();
    
    TClonesArray& trkTca=(*ntp.trk);
    const Int_t numTrks=trkTca.GetEntriesFast();

    TClonesArray& shwTca=(*ntp.shw);
    const Int_t numShws=shwTca.GetEntriesFast();

    TClonesArray& stpTca=(*ntp.stp);
    const Int_t numStps=stpTca.GetEntries();

    TClonesArray& slcTca=(*ntp.slc);
    const Int_t numSlcs=slcTca.GetEntries();
    
    Float_t totalSnarlAdc=0;
    Float_t totalEvtAdc=0;
    Float_t totalSlcAdc=0;
    Float_t thisSlcAdc=0;
    
    for (Int_t i=0;i<numStps;i++) {
      const NtpSRStrip* pstp=
        dynamic_cast<NtpSRStrip*>(stpTca[i]);
      const NtpSRStrip& stp=(*pstp);
      totalSnarlAdc+=stp.ph0.raw+stp.ph1.raw;
    }
    
    if (numEvts==0 && numTrks>0) {
      cout<<"Ahhh, num events in tca="<<numEvts
          <<", trks="<<numTrks<<endl;
    }
    
    if (numEvts>1) {
      MAXMSG("MeuAnalysis",Msg::kInfo,10)
        <<"High num evts="<<numEvts<<", trks="<<numTrks<<endl;
    }

    if (numSlcs!=1) {
      MAXMSG("MeuAnalysis",Msg::kInfo,10)
        <<"High num slices in tca="<<numSlcs
        <<", trks="<<numTrks<<endl;
      //continue;
    }
    
    MAXMSG("MeuAnalysis",Msg::kInfo,100)
      <<"Num events in tca="<<numEvts
      <<", slcs="<<numSlcs<<", trks="<<numTrks<<endl;
    
    Int_t slcStps=0;
    
    //loop over the slices
    for (Int_t islc=0;islc<numSlcs;++islc) {
      const NtpSRSlice* pslc=
        dynamic_cast<NtpSRSlice*>(slcTca[islc]);
      const NtpSRSlice& slc=(*pslc);

      MAXMSG("MeuAnalysis",Msg::kDebug,200)
        <<"  slice "<<islc+1<<" of "<<numSlcs<<endl;
      
      slcStps+=slc.nstrip;

      //loop over strips in slc
      for (Int_t i=0;i<slc.nstrip;++i) {
        const NtpSRStrip* pstp=
          dynamic_cast<NtpSRStrip*>(stpTca[slc.stp[i]]);
        const NtpSRStrip& stp=(*pstp);
        totalSlcAdc+=stp.ph0.raw+stp.ph1.raw;
      }

      //loop over tracks in slc  ---- can't do this!!!
      //for(int i=0;i<slc.ntrack;i++){    
      //const NtpSRTrack* ptrk=
      //  dynamic_cast<NtpSRTrack*>(trkTca[slc.trk[i]]);
      //const NtpSRTrack& trk=(*ptrk);

      //trkStps=trk.nstrip;
      //MAXMSG("MeuAnalysis",Msg::kInfo,300)
      //  <<"    track "<<i+1<<" of "<<slc.ntrack
      //  <<", trkStps="<<trkStps<<endl;
      //}
    }

    map<Int_t,Int_t> evtPerSlc;
    
    //loop over the events
    for (Int_t ntpevt=0;ntpevt<numEvts;++ntpevt) {
      const NtpSREvent* pevt=
        dynamic_cast<NtpSREvent*>(evtTca[ntpevt]);
      const NtpSREvent& evt=(*pevt);
      
      MAXMSG("MeuAnalysis",Msg::kInfo,100)
        <<"Entry "<<entry<<" has tracks="<<evt.ntrack
        <<", shws="<<numShws<<", slcs="<<numSlcs<<endl;

      //get the slice associated with this event
      const NtpSRSlice* pslc=
        dynamic_cast<NtpSRSlice*>(slcTca[evt.slc]);
      const NtpSRSlice& slc=(*pslc);

      evtPerSlc[evt.slc]++;

      //check that the filter works
      Int_t evts=-1;
      Bool_t goodEvt=!cuts.FilterBadEvtPerSlc(ntp,evt.slc,entry,&evts);
      if (goodEvt){
        MAXMSG("MeuAnalysis",Msg::kInfo,100)
          <<"TEST::entry="<<entry
          <<", evt="<<ntpevt<<" is good, evtsPerSlc="<<evts
          <<" (slc="<<evt.slc<<")"<<endl;
      }
      else{
        MAXMSG("MeuAnalysis",Msg::kInfo,100)
          <<"TEST::entry="<<entry
          <<", evt="<<ntpevt<<" is bad, evtsPerSlc="<<evts
          <<" (slc="<<evt.slc<<")"<<endl;
      }
      
      if (evt.ntrack!=1 || evt.nshower>1) continue;

      Int_t trkStps=0;
      Int_t shwStps=0;
      
      //loop over strips in slc
      for (Int_t i=0;i<slc.nstrip;++i) {
        const NtpSRStrip* pstp=
          dynamic_cast<NtpSRStrip*>(stpTca[slc.stp[i]]);
        const NtpSRStrip& stp=(*pstp);
        thisSlcAdc+=stp.ph0.raw+stp.ph1.raw;
      }
      
      //loop over tracks in evt
      for(int i=0;i<evt.ntrack;i++) {    
        const NtpSRTrack* ptrk=
          dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[i]]);
        const NtpSRTrack& trk=(*ptrk);

        trkStps=trk.nstrip;
        
        for (Int_t i=0;i<trk.nstrip;i++) {

          //check for bug where strip index is -1
          if (trk.stp[i]<0) {
            MAXMSG("MeuAnalysis",Msg::kInfo,500)
              <<"Skipping strip with trk.stp[i]="<<trk.stp[i]<<endl;
            continue;
          }

          const NtpSRStrip* pstp=
            dynamic_cast<NtpSRStrip*>(stpTca[trk.stp[i]]);
          const NtpSRStrip& stp=(*pstp);
            MAXMSG("MeuAnalysis",Msg::kDebug,100)
              <<"strip time="<<stp.time0<<endl;
        }
      }

      //loop over showers in evt
      for(int i=0;i<evt.nshower;i++) {    
        const NtpSRShower* pshw=
          dynamic_cast<NtpSRShower*>(shwTca[evt.shw[i]]);
        const NtpSRShower& shw=(*pshw);
        
        shwStps=shw.nstrip;
      }

      //loop over strips in evt
      for (Int_t i=0;i<evt.nstrip;i++) {
        const NtpSRStrip* pstp=
          dynamic_cast<NtpSRStrip*>(stpTca[evt.stp[i]]);
        const NtpSRStrip& stp=(*pstp);
        totalEvtAdc+=stp.ph0.raw+stp.ph1.raw;
      }
      MAXMSG("MeuAnalysis",Msg::kInfo,100)
        <<"e="<<entry<<", evt="<<ntpevt
        <<", ADCs: snl="<<totalSnarlAdc
        <<", evt="<<totalEvtAdc
        <<", slc="<<totalSlcAdc
        <<", tslc="<<thisSlcAdc
        <<", #st="<<numStps<<","<<evt.nstrip<<","<<slcStps
        //<<", trk="<<trkStps<<", shw="<<shwStps
        <<endl;
    }

    //make this into a function
    //Bool_t FilterMultiEvtPerSlc(ntp,slc,entry)
    typedef map<Int_t,Int_t>::const_iterator slcIt;
    MAXMSG("MeuAnalysis",Msg::kInfo,30)
      <<"MAIN::Events per slice for entry="<<entry<<endl;
    for (slcIt it=evtPerSlc.begin();it!=evtPerSlc.end();++it) {
      MAXMSG("MeuAnalysis",Msg::kInfo,200)
        <<"  slice "<<it->first<<" has "<<it->second<<" event(s)"<<endl;
    }
    
  }//end of for                                       
    
  
  MSG("MeuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl;

  MSG("MeuAnalysis",Msg::kInfo) 
    <<" ** Finished TestEventLoop method **"<<endl;
}

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

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

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

  TH1F* hTrkPl=new TH1F("hTrkPl","hTrkPl",201,-1,200);
  hTrkPl->SetTitle("hTrkPl");
  hTrkPl->GetXaxis()->SetTitle("hTrkPl");
  hTrkPl->GetXaxis()->CenterTitle();
  hTrkPl->SetFillColor(0);
  hTrkPl->SetLineWidth(3);
  //hTrkPl->SetBit(TH1::kCanRebin);

  TH1F* hTrkPlN_1=new TH1F("hTrkPlN_1","hTrkPlN_1",201,-1,200);
  hTrkPlN_1->SetTitle("hTrkPlN_1");
  hTrkPlN_1->GetXaxis()->SetTitle("hTrkPlN_1");
  hTrkPlN_1->GetXaxis()->CenterTitle();
  hTrkPlN_1->SetLineColor(2);//red for N-1 cuts
  hTrkPlN_1->SetFillColor(0);
  hTrkPlN_1->SetLineWidth(3);
  //hTrkPlN_1->SetBit(TH1::kCanRebin);

  TH1F* hTrkPlN=new TH1F("hTrkPlN","hTrkPlN",201,-1,200);
  hTrkPlN->SetTitle("hTrkPlN");
  hTrkPlN->GetXaxis()->SetTitle("hTrkPlN");
  hTrkPlN->GetXaxis()->CenterTitle();
  hTrkPlN->SetLineColor(4);//blue for N cuts
  hTrkPlN->SetFillColor(0);
  hTrkPlN->SetLineWidth(3);
  //hTrkPlN->SetBit(TH1::kCanRebin);

  TH1F* hDistEdge=new TH1F("hDistEdge","hDistEdge",201,-3,3);
  hDistEdge->SetTitle("hDistEdge");
  hDistEdge->GetXaxis()->SetTitle("hDistEdge");
  hDistEdge->GetXaxis()->CenterTitle();
  hDistEdge->SetFillColor(0);
  hDistEdge->SetLineWidth(3);
  //hDistEdge->SetBit(TH1::kCanRebin);

  TH1F* hDistEdgeN_1=new TH1F("hDistEdgeN_1","hDistEdgeN_1",201,-3,3);
  hDistEdgeN_1->SetTitle("hDistEdgeN_1");
  hDistEdgeN_1->GetXaxis()->SetTitle("hDistEdgeN_1");
  hDistEdgeN_1->GetXaxis()->CenterTitle();
  hDistEdgeN_1->SetLineColor(2);//red for N-1 cuts
  hDistEdgeN_1->SetFillColor(0);
  hDistEdgeN_1->SetLineWidth(3);
  //hDistEdgeN_1->SetBit(TH1::kCanRebin);

  TH1F* hDistEdgeN=new TH1F("hDistEdgeN","hDistEdgeN",201,-3,3);
  hDistEdgeN->SetTitle("hDistEdgeN");
  hDistEdgeN->GetXaxis()->SetTitle("hDistEdgeN");
  hDistEdgeN->GetXaxis()->CenterTitle();
  hDistEdgeN->SetLineColor(4);//blue for N cuts
  hDistEdgeN->SetFillColor(0);
  hDistEdgeN->SetLineWidth(3);
  //hDistEdgeN->SetBit(TH1::kCanRebin);

  TH1F* hViewsVtx=new TH1F("hViewsVtx","hViewsVtx",101,-1,100);
  hViewsVtx->SetTitle("hViewsVtx");
  hViewsVtx->GetXaxis()->SetTitle("hViewsVtx");
  hViewsVtx->GetXaxis()->CenterTitle();
  hViewsVtx->SetFillColor(0);
  hViewsVtx->SetLineWidth(3);
  //hViewsVtx->SetBit(TH1::kCanRebin);

  TH1F* hViewsVtxN_1=new TH1F("hViewsVtxN_1","hViewsVtxN_1",101,-1,100);
  hViewsVtxN_1->SetTitle("hViewsVtxN_1");
  hViewsVtxN_1->GetXaxis()->SetTitle("hViewsVtxN_1");
  hViewsVtxN_1->GetXaxis()->CenterTitle();
  hViewsVtxN_1->SetLineColor(2);//red for N-1 cuts
  hViewsVtxN_1->SetFillColor(0);
  hViewsVtxN_1->SetLineWidth(3);
  //hViewsVtxN_1->SetBit(TH1::kCanRebin);

  TH1F* hViewsVtxN=new TH1F("hViewsVtxN","hViewsVtxN",101,-1,100);
  hViewsVtxN->SetTitle("hViewsVtxN");
  hViewsVtxN->GetXaxis()->SetTitle("hViewsVtxN");
  hViewsVtxN->GetXaxis()->CenterTitle();
  hViewsVtxN->SetLineColor(4);//blue for N cuts
  hViewsVtxN->SetFillColor(0);
  hViewsVtxN->SetLineWidth(3);
  //hViewsVtxN->SetBit(TH1::kCanRebin);

  TH1F* hViewsEnd=new TH1F("hViewsEnd","hViewsEnd",101,-1,100);
  hViewsEnd->SetTitle("hViewsEnd");
  hViewsEnd->GetXaxis()->SetTitle("hViewsEnd");
  hViewsEnd->GetXaxis()->CenterTitle();
  hViewsEnd->SetFillColor(0);
  hViewsEnd->SetLineWidth(3);
  //hViewsEnd->SetBit(TH1::kCanRebin);

  TH1F* hViewsEndN_1=new TH1F("hViewsEndN_1","hViewsEndN_1",101,-1,100);
  hViewsEndN_1->SetTitle("hViewsEndN_1");
  hViewsEndN_1->GetXaxis()->SetTitle("hViewsEndN_1");
  hViewsEndN_1->GetXaxis()->CenterTitle();
  hViewsEndN_1->SetLineColor(2);//red for N-1 cuts
  hViewsEndN_1->SetFillColor(0);
  hViewsEndN_1->SetLineWidth(3);
  //hViewsEndN_1->SetBit(TH1::kCanRebin);

  TH1F* hViewsEndN=new TH1F("hViewsEndN","hViewsEndN",101,-1,100);
  hViewsEndN->SetTitle("hViewsEndN");
  hViewsEndN->GetXaxis()->SetTitle("hViewsEndN");
  hViewsEndN->GetXaxis()->CenterTitle();
  hViewsEndN->SetLineColor(4);//blue for N cuts
  hViewsEndN->SetFillColor(0);
  hViewsEndN->SetLineWidth(3);
  //hViewsEndN->SetBit(TH1::kCanRebin);

  TH1F* hDistEndSt=new TH1F("hDistEndSt","hDistEndSt",200,-0.1,2);
  hDistEndSt->SetTitle("hDistEndSt");
  hDistEndSt->GetXaxis()->SetTitle("hDistEndSt");
  hDistEndSt->GetXaxis()->CenterTitle();
  hDistEndSt->SetFillColor(0);
  hDistEndSt->SetLineWidth(3);
  //hDistEndSt->SetBit(TH1::kCanRebin);

  TH1F* hDistEndStN_1=new TH1F("hDistEndStN_1","hDistEndStN_1",
                               200,-0.1,2);
  hDistEndStN_1->SetTitle("hDistEndStN_1");
  hDistEndStN_1->GetXaxis()->SetTitle("hDistEndStN_1");
  hDistEndStN_1->GetXaxis()->CenterTitle();
  hDistEndStN_1->SetLineColor(2);//red for N-1 cuts
  hDistEndStN_1->SetFillColor(0);
  hDistEndStN_1->SetLineWidth(3);
  //hDistEndStN_1->SetBit(TH1::kCanRebin);

  TH1F* hDistEndStN=new TH1F("hDistEndStN","hDistEndStN",200,-0.1,2);
  hDistEndStN->SetTitle("hDistEndStN");
  hDistEndStN->GetXaxis()->SetTitle("hDistEndStN");
  hDistEndStN->GetXaxis()->CenterTitle();
  hDistEndStN->SetLineColor(4);//blue for N cuts
  hDistEndStN->SetFillColor(0);
  hDistEndStN->SetLineWidth(3);
  //hDistEndStN->SetBit(TH1::kCanRebin);

  TH1F* hSigBeyondEnd=new TH1F("hSigBeyondEnd","hSigBeyondEnd",1000,-1,50000);
  hSigBeyondEnd->SetTitle("hSigBeyondEnd");
  hSigBeyondEnd->GetXaxis()->SetTitle("hSigBeyondEnd");
  hSigBeyondEnd->GetXaxis()->CenterTitle();
  hSigBeyondEnd->SetFillColor(0);
  hSigBeyondEnd->SetLineWidth(3);
  //hSigBeyondEnd->SetBit(TH1::kCanRebin);
  
  TH1F* hSigBeyondEndN_1=new TH1F("hSigBeyondEndN_1","hSigBeyondEndN_1",
                                  1000,-1,50000);
  hSigBeyondEndN_1->SetTitle("hSigBeyondEndN_1");
  hSigBeyondEndN_1->GetXaxis()->SetTitle("hSigBeyondEndN_1");
  hSigBeyondEndN_1->GetXaxis()->CenterTitle();
  hSigBeyondEndN_1->SetLineColor(2);//red for N-1 cuts
  hSigBeyondEndN_1->SetFillColor(0);
  hSigBeyondEndN_1->SetLineWidth(3);
  //hSigBeyondEndN_1->SetBit(TH1::kCanRebin);

  TH1F* hSigBeyondEndN=new TH1F("hSigBeyondEndN","hSigBeyondEndN",
                                1000,-1,50000);
  hSigBeyondEndN->SetTitle("hSigBeyondEndN");
  hSigBeyondEndN->GetXaxis()->SetTitle("hSigBeyondEndN");
  hSigBeyondEndN->GetXaxis()->CenterTitle();
  hSigBeyondEndN->SetLineColor(4);//blue for N cuts
  hSigBeyondEndN->SetFillColor(0);
  hSigBeyondEndN->SetLineWidth(3);
  //hSigBeyondEndN->SetBit(TH1::kCanRebin);


  TH1F* hSigTrkEndGaps=new TH1F("hSigTrkEndGaps","hSigTrkEndGaps",1000,-1,50000);
  hSigTrkEndGaps->SetTitle("hSigTrkEndGaps");
  hSigTrkEndGaps->GetXaxis()->SetTitle("hSigTrkEndGaps");
  hSigTrkEndGaps->GetXaxis()->CenterTitle();
  hSigTrkEndGaps->SetFillColor(0);
  hSigTrkEndGaps->SetLineWidth(3);
  //hSigTrkEndGaps->SetBit(TH1::kCanRebin);
  
  TH1F* hSigTrkEndGapsN_1=new TH1F("hSigTrkEndGapsN_1","hSigTrkEndGapsN_1",
                                  1000,-1,50000);
  hSigTrkEndGapsN_1->SetTitle("hSigTrkEndGapsN_1");
  hSigTrkEndGapsN_1->GetXaxis()->SetTitle("hSigTrkEndGapsN_1");
  hSigTrkEndGapsN_1->GetXaxis()->CenterTitle();
  hSigTrkEndGapsN_1->SetLineColor(2);//red for N-1 cuts
  hSigTrkEndGapsN_1->SetFillColor(0);
  hSigTrkEndGapsN_1->SetLineWidth(3);
  //hSigTrkEndGapsN_1->SetBit(TH1::kCanRebin);

  TH1F* hSigTrkEndGapsN=new TH1F("hSigTrkEndGapsN","hSigTrkEndGapsN",
                                1000,-1,50000);
  hSigTrkEndGapsN->SetTitle("hSigTrkEndGapsN");
  hSigTrkEndGapsN->GetXaxis()->SetTitle("hSigTrkEndGapsN");
  hSigTrkEndGapsN->GetXaxis()->CenterTitle();
  hSigTrkEndGapsN->SetLineColor(4);//blue for N cuts
  hSigTrkEndGapsN->SetFillColor(0);
  hSigTrkEndGapsN->SetLineWidth(3);
  //hSigTrkEndGapsN->SetBit(TH1::kCanRebin);

  TH1F* hMatTrav=new TH1F("hMatTrav","hMatTrav",400,0,10);
  hMatTrav->SetTitle("hMatTrav");
  hMatTrav->GetXaxis()->SetTitle("hMatTrav");
  hMatTrav->GetXaxis()->CenterTitle();
  hMatTrav->SetFillColor(0);
  hMatTrav->SetLineWidth(3);
  //hMatTrav->SetBit(TH1::kCanRebin);

  TH1F* hMatTravN_1=new TH1F("hMatTravN_1","hMatTravN_1",400,0,10);
  hMatTravN_1->SetTitle("hMatTravN_1");
  hMatTravN_1->GetXaxis()->SetTitle("hMatTravN_1");
  hMatTravN_1->GetXaxis()->CenterTitle();
  hMatTravN_1->SetLineColor(2);//red for N-1 cuts
  hMatTravN_1->SetFillColor(0);
  hMatTravN_1->SetLineWidth(3);
  //hMatTravN_1->SetBit(TH1::kCanRebin);

  TH1F* hMatTravN=new TH1F("hMatTravN","hMatTravN",400,0,10);
  hMatTravN->SetTitle("hMatTravN");
  hMatTravN->GetXaxis()->SetTitle("hMatTravN");
  hMatTravN->GetXaxis()->CenterTitle();
  hMatTravN->SetLineColor(4);//blue for N cuts
  hMatTravN->SetFillColor(0);
  hMatTravN->SetLineWidth(3);
  //hMatTravN->SetBit(TH1::kCanRebin);
        
  //reconstruction object
  MeuReco reco;
  MeuCuts cuts;
  MeuSummary s;

  
  this->InitialiseLoopVariables();  
  
  //for(Int_t entry=0;entry<10000;entry++){
  for(Int_t entry=0;entry<fEntries;entry++){
      
    this->SetLoopVariables(entry);
      
    NtpStRecord& ntp=(*fRec);//Make a reference instead of a pointer

    //ensure there is only 1 track in the entry
    TClonesArray& trkTca=(*ntp.trk);
    Int_t numTrks=trkTca.GetEntries();
    if (numTrks!=1) continue;
    const NtpSRTrack* ptrk=dynamic_cast<NtpSRTrack*>(trkTca[0]);
    const NtpSRTrack& trk=(*ptrk);
   
    //cut on number of planes
    Int_t numPlanes=abs(trk.vtx.plane-trk.end.plane);
    if (numPlanes<8) continue;

    s.Reset();
    s.RFid=1.0;
    s.DistToEdgeFid=0.4;
    s.Event=entry;
    cuts.FillSTSumDetails(ntp,s,-1,-1);
    map<Int_t,MeuHitInfo> plInfo;

    //check if the track is in the spectrometer
    if (trk.vtx.plane>120 && trk.end.plane>120) continue;

    TClonesArray& evtTca=(*ntp.evt);
    const NtpSREvent* pevt=dynamic_cast<NtpSREvent*>(evtTca[0]);
    cuts.ExtractPlInfo(ntp,s,plInfo,*pevt);
    
    reco.CalcVtx(s,plInfo);

    //check if track stops in spectrometer
    if (s.EndPlane>120) continue;

    //calculate the containment (before doing the vtxSpecial)
    reco.CalcFCPC(s,plInfo);

    //recalculate the vtx of the tracks coming in from the spectrometer
    reco.CalcVtxSpecialND(s,plInfo);

    //cut on number of planes with potential new vertex
    
    numPlanes=abs(s.VtxPlane-s.EndPlane);
    if (numPlanes<8) continue;

    reco.CalcStripDists(s,plInfo);
    
    //calc distance to detector edge
    Float_t dist=reco.CalcSmallestDeepDistToEdge(plInfo[s.EndPlane]);

    Int_t diffAtEnd=-1;
    Int_t diffAtVtx=-1;
    cuts.CheckTrkViews(s,&diffAtVtx,&diffAtEnd);

    Float_t minDist=-1;
    Float_t sigBeyondEnd=-1;
    Float_t adcInGap=-1;
    
    Bool_t goodContainment=false;
    Bool_t goodReco=false;
    if (s.PC) {
      goodReco=reco.Reconstruct(s,plInfo);
      if (goodReco) {

        TClonesArray& evtTca=(*ntp.evt);
        const NtpSREvent* pevt=dynamic_cast<NtpSREvent*>(evtTca[0]);
        this->Recalibrate(ntp,s,plInfo,*pevt);
        reco.CalcWindow(s,plInfo);

        //now check that window is deeply contained
        if (s.WinSigMap>0) {
          goodContainment=reco.CheckWinContainment(s,plInfo);
        }

        //now do some of the calculations that require good PC
        if (goodContainment){
          cuts.FilterBadDistEndStrip(s,plInfo,&minDist);
          cuts.FilterBadTrackEnd(s,plInfo,&sigBeyondEnd);
          cuts.FilterBadEndGaps(s,plInfo,&adcInGap);
        }
      }
    }


    //if (s.WinSigMap<=0) { //matTrav
    
    
    //start of "no cuts"
    hTrkPl->Fill(numPlanes);
    hDistEdge->Fill(dist);
    hViewsVtx->Fill(diffAtVtx);
    hViewsEnd->Fill(diffAtEnd);
    if (goodContainment){ //only have a value in this case
      hDistEndSt->Fill(minDist);
      hSigBeyondEnd->Fill(sigBeyondEnd);
      hSigTrkEndGaps->Fill(adcInGap);
    }
    //only have a value in this case (requires calcWindow)
    if (goodReco) hMatTrav->Fill(s.TotalMatTraversed);
    
    //cuts
    Bool_t goodDist=dist>0.4;
    Bool_t goodTrkPl=numPlanes>10;
    Bool_t goodViewsVtx=diffAtVtx<=9;//for ND full inst. planes
    Bool_t goodViewsEnd=diffAtEnd<=5;
    Bool_t goodDistEndSt=minDist>0.2;
    Bool_t goodSigBeyondEnd=sigBeyondEnd<500 && sigBeyondEnd>=0;
    Bool_t goodSigTrkEndGaps=adcInGap<500 && adcInGap>=0;
    Bool_t goodMatTrav=!cuts.FilterLowMatTraversed(s);
  
    Bool_t goodN=
      goodContainment && 
      goodDist && 
      goodTrkPl && 
      goodViewsVtx &&
      goodViewsEnd &&
      goodDistEndSt &&
      goodSigBeyondEnd &&
      goodSigTrkEndGaps &&
      goodMatTrav &&
      true;//for coding convenience
    
    //start of N-1 cuts plots
    if (goodContainment && 
        goodDist &&
        goodViewsVtx &&
        goodViewsEnd &&
        goodDistEndSt &&
        goodSigBeyondEnd &&
        goodSigTrkEndGaps &&
        goodMatTrav &&
        true) {
      hTrkPlN_1->Fill(numPlanes);

      if (numPlanes<10){
        MAXMSG("MeuAnalysis",Msg::kInfo,10)
          <<"vtxPl="<<s.VtxPlane
          <<", endPl="<<s.EndPlane
          <<", stopPl="<<s.WinVtxSidePl
          <<", startPl="<<s.WinStopSidePl
          <<", meu="<<s.WinSigMap
          <<", PLC="<<s.WinAvPLCor
          <<endl;
      }
    }
    if (goodContainment && 
        goodTrkPl &&
        goodViewsVtx &&
        goodViewsEnd &&
        goodDistEndSt &&
        goodSigBeyondEnd &&
        goodSigTrkEndGaps &&
        goodMatTrav &&
        true) {
      hDistEdgeN_1->Fill(dist);
    }
    if (goodContainment && 
        goodDist &&
        goodTrkPl &&
        goodViewsEnd &&
        goodDistEndSt &&
        goodSigBeyondEnd &&
        goodSigTrkEndGaps &&
        goodMatTrav &&
        true) {
      hViewsVtxN_1->Fill(diffAtVtx);
    }
    if (goodContainment && 
        goodDist &&
        goodTrkPl &&
        goodViewsVtx &&
        goodDistEndSt &&
        goodSigBeyondEnd &&
        goodSigTrkEndGaps &&
        goodMatTrav &&
        true) {
      hViewsEndN_1->Fill(diffAtEnd);
    }
    if (goodContainment && 
        goodDist &&
        goodTrkPl &&
        goodViewsVtx &&
        goodViewsEnd &&
        goodSigBeyondEnd &&
        goodSigTrkEndGaps &&
        goodMatTrav &&
        true) {
      hDistEndStN_1->Fill(minDist);
    }
    if (goodContainment && 
        goodDist &&
        goodTrkPl &&
        goodViewsVtx &&
        goodViewsEnd &&
        goodDistEndSt &&
        goodSigTrkEndGaps &&
        goodMatTrav &&
        true) {
      hSigBeyondEndN_1->Fill(sigBeyondEnd);
    }
    if (goodContainment && 
        goodDist &&
        goodTrkPl &&
        goodViewsVtx &&
        goodViewsEnd &&
        goodDistEndSt &&
        goodSigBeyondEnd &&
        goodMatTrav &&
        true) {
      hSigTrkEndGapsN_1->Fill(adcInGap);
    }
    if (goodContainment && 
        goodDist &&
        goodTrkPl &&
        goodViewsVtx &&
        goodViewsEnd &&
        goodDistEndSt &&
        goodSigBeyondEnd &&
        goodSigTrkEndGaps &&
        true) {
      hMatTravN_1->Fill(s.TotalMatTraversed);
    }    

    //start of N cuts plots
    if (goodN) {
      hTrkPlN->Fill(numPlanes);
      hDistEdgeN->Fill(dist);
      hViewsVtxN->Fill(diffAtVtx);
      hViewsEndN->Fill(diffAtEnd);
      hDistEndStN->Fill(minDist);
      hSigBeyondEndN->Fill(sigBeyondEnd);
      hSigTrkEndGapsN->Fill(adcInGap);
      hMatTravN->Fill(s.TotalMatTraversed);
    }
  }//end of for                                       
  

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

  MSG("MeuAnalysis",Msg::kInfo) 
    <<" ** Finished N_1Plots method **"<<endl;
}

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

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

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


  TH1F* hNTrack=new TH1F("hNTrack","hNTrack",10,-1,9);
  hNTrack->SetTitle("NTrack");
  hNTrack->GetXaxis()->SetTitle("NTrack");
  hNTrack->GetXaxis()->CenterTitle();
  hNTrack->GetYaxis()->SetTitle("");
  hNTrack->GetYaxis()->CenterTitle();
  hNTrack->SetFillColor(0);
  //hNTrack->SetBit(TH1::kCanRebin);

  TH1F* hStPerTrk=new TH1F("hStPerTrk","hStPerTrk",201,-1,200);
  hStPerTrk->SetTitle("StPerTrk");
  hStPerTrk->GetXaxis()->SetTitle("StPerTrk");
  hStPerTrk->GetXaxis()->CenterTitle();
  hStPerTrk->GetYaxis()->SetTitle("");
  hStPerTrk->GetYaxis()->CenterTitle();
  hStPerTrk->SetFillColor(0);
  //hStPerTrk->SetBit(TH1::kCanRebin);

  TH1F* hTPos=new TH1F("hTPos","hTPos",200,-4,4);
  hTPos->SetTitle("TPos");
  hTPos->GetXaxis()->SetTitle("TPos");
  hTPos->GetXaxis()->CenterTitle();
  hTPos->GetYaxis()->SetTitle("");
  hTPos->GetYaxis()->CenterTitle();
  hTPos->SetFillColor(0);
  //hTPos->SetBit(TH1::kCanRebin);

  TH1F* hSigCor=new TH1F("hSigCor","hSigCor",3000,-1,15000);
  hSigCor->SetTitle("SigCor");
  hSigCor->GetXaxis()->SetTitle("SigCor");
  hSigCor->GetXaxis()->CenterTitle();
  hSigCor->GetYaxis()->SetTitle("");
  hSigCor->GetYaxis()->CenterTitle();
  hSigCor->SetFillColor(0);
  //hSigCor->SetBit(TH1::kCanRebin);

  TH1F* hView=new TH1F("hView","hView",6,-1,5);
  hView->SetTitle("View");
  hView->GetXaxis()->SetTitle("View");
  hView->GetXaxis()->CenterTitle();
  hView->GetYaxis()->SetTitle("");
  hView->GetYaxis()->CenterTitle();
  hView->SetFillColor(0);
  //hView->SetBit(TH1::kCanRebin);

  TH1F* hPlane=new TH1F("hPlane","hPlane",301,-1,300);
  hPlane->SetTitle("Plane");
  hPlane->GetXaxis()->SetTitle("Plane");
  hPlane->GetXaxis()->CenterTitle();
  hPlane->GetYaxis()->SetTitle("");
  hPlane->GetYaxis()->CenterTitle();
  hPlane->SetFillColor(0);
  //hPlane->SetBit(TH1::kCanRebin);

  TH1F* hStrip=new TH1F("hStrip","hStrip",101,-1,100);
  hStrip->SetTitle("Strip");
  hStrip->GetXaxis()->SetTitle("Strip");
  hStrip->GetXaxis()->CenterTitle();
  hStrip->GetYaxis()->SetTitle("");
  hStrip->GetYaxis()->CenterTitle();
  hStrip->SetFillColor(0);
  //hStrip->SetBit(TH1::kCanRebin);

  TH1F* hZ=new TH1F("hZ","hZ",100,-1,30);
  hZ->SetTitle("Z");
  hZ->GetXaxis()->SetTitle("Z");
  hZ->GetXaxis()->CenterTitle();
  hZ->GetYaxis()->SetTitle("");
  hZ->GetYaxis()->CenterTitle();
  hZ->SetFillColor(0);
  //hZ->SetBit(TH1::kCanRebin);

  TH1F* hR=new TH1F("hR","hR",100,-1,10);
  hR->SetTitle("R");
  hR->GetXaxis()->SetTitle("R");
  hR->GetXaxis()->CenterTitle();
  hR->GetYaxis()->SetTitle("");
  hR->GetYaxis()->CenterTitle();
  hR->SetFillColor(0);
  //hR->SetBit(TH1::kCanRebin);

  TH1F* hMaxR=new TH1F("hMaxR","hMaxR",100,-1,10);
  hMaxR->SetTitle("MaxR");
  hMaxR->GetXaxis()->SetTitle("MaxR");
  hMaxR->GetXaxis()->CenterTitle();
  hMaxR->GetYaxis()->SetTitle("");
  hMaxR->GetYaxis()->CenterTitle();
  hMaxR->SetFillColor(0);
  //hMaxR->SetBit(TH1::kCanRebin);

  TH1F* hMinR=new TH1F("hMinR","hMinR",100,-1,10);
  hMinR->SetTitle("MinR");
  hMinR->GetXaxis()->SetTitle("MinR");
  hMinR->GetXaxis()->CenterTitle();
  hMinR->GetYaxis()->SetTitle("");
  hMinR->GetYaxis()->CenterTitle();
  hMinR->SetFillColor(0);
  //hMinR->SetBit(TH1::kCanRebin);

  TH2F* hTPosVsPlane=new TH2F("hTPosVsPlane","hTPosVsPlane",
                              301,-1,300,100,-3,3);
  hTPosVsPlane->SetTitle("TPos vs Plane");
  hTPosVsPlane->GetXaxis()->SetTitle("Plane");
  hTPosVsPlane->GetXaxis()->CenterTitle();
  hTPosVsPlane->GetYaxis()->SetTitle("TPos");
  hTPosVsPlane->GetYaxis()->CenterTitle();
  hTPosVsPlane->SetFillColor(0);
  //hTPosVsPlane->SetBit(TH1::kCanRebin);

  TH2F* hYvsX=new TH2F("hYvsX","hYvsX",
                       100,-5,5,100,-5,5);
  hYvsX->SetTitle("Y vs X");
  hYvsX->GetXaxis()->SetTitle("X");
  hYvsX->GetXaxis()->CenterTitle();
  hYvsX->GetYaxis()->SetTitle("Y");
  hYvsX->GetYaxis()->CenterTitle();
  hYvsX->SetFillColor(0);
  //hYvsX->SetBit(TH1::kCanRebin);

  TH2F* hYvsXVtx=new TH2F("hYvsXVtx","hYvsXVtx",
                          100,-5,5,100,-5,5);
  hYvsXVtx->SetTitle("Y vs X");
  hYvsXVtx->GetXaxis()->SetTitle("X");
  hYvsXVtx->GetXaxis()->CenterTitle();
  hYvsXVtx->GetYaxis()->SetTitle("Y");
  hYvsXVtx->GetYaxis()->CenterTitle();
  hYvsXVtx->SetFillColor(0);
  //hYvsXVtx->SetBit(TH1::kCanRebin);

  TH2F* hYvsXEnd=new TH2F("hYvsXEnd","hYvsXEnd",
                          100,-5,5,100,-5,5);
  hYvsXEnd->SetTitle("Y vs X");
  hYvsXEnd->GetXaxis()->SetTitle("X");
  hYvsXEnd->GetXaxis()->CenterTitle();
  hYvsXEnd->GetYaxis()->SetTitle("Y");
  hYvsXEnd->GetYaxis()->CenterTitle();
  hYvsXEnd->SetFillColor(0);
  //hYvsXEnd->SetBit(TH1::kCanRebin);

  TH2F* hYvsXVtxFid=new TH2F("hYvsXVtxFid","hYvsXVtxFid",
                             100,-5,5,100,-5,5);
  hYvsXVtxFid->SetTitle("Y vs X");
  hYvsXVtxFid->GetXaxis()->SetTitle("X");
  hYvsXVtxFid->GetXaxis()->CenterTitle();
  hYvsXVtxFid->GetYaxis()->SetTitle("Y");
  hYvsXVtxFid->GetYaxis()->CenterTitle();
  hYvsXVtxFid->SetFillColor(0);
  //hYvsXVtxFid->SetBit(TH1::kCanRebin);

  TH2F* hYvsXEndFid=new TH2F("hYvsXEndFid","hYvsXEndFid",
                             100,-5,5,100,-5,5);
  hYvsXEndFid->SetTitle("Y vs X");
  hYvsXEndFid->GetXaxis()->SetTitle("X");
  hYvsXEndFid->GetXaxis()->CenterTitle();
  hYvsXEndFid->GetYaxis()->SetTitle("Y");
  hYvsXEndFid->GetYaxis()->CenterTitle();
  hYvsXEndFid->SetFillColor(0);
  //hYvsXEndFid->SetBit(TH1::kCanRebin);

  TH2F* hYvsZ=new TH2F("hYvsZ","hYvsZ",400,-1,30,100,-5,5);
  hYvsZ->SetTitle("Y vs Z");
  hYvsZ->GetXaxis()->SetTitle("Z");
  hYvsZ->GetXaxis()->CenterTitle();
  hYvsZ->GetYaxis()->SetTitle("Y");
  hYvsZ->GetYaxis()->CenterTitle();
  hYvsZ->SetFillColor(0);
  //hYvsZ->SetBit(TH1::kCanRebin);

  TH2F* hXvsZ=new TH2F("hXvsZ","hXvsZ",400,-1,30,100,-5,5);
  hXvsZ->SetTitle("X vs Z");
  hXvsZ->GetXaxis()->SetTitle("Z");
  hXvsZ->GetXaxis()->CenterTitle();
  hXvsZ->GetYaxis()->SetTitle("Y");
  hXvsZ->GetYaxis()->CenterTitle();
  hXvsZ->SetFillColor(0);
  //hXvsZ->SetBit(TH1::kCanRebin);

  TH2F* hYvsPlane=new TH2F("hYvsPlane","hYvsPlane",301,-1,300,100,-5,5);
  hYvsPlane->SetTitle("Y vs Z");
  hYvsPlane->GetXaxis()->SetTitle("Z");
  hYvsPlane->GetXaxis()->CenterTitle();
  hYvsPlane->GetYaxis()->SetTitle("Y");
  hYvsPlane->GetYaxis()->CenterTitle();
  hYvsPlane->SetFillColor(0);
  //hYvsPlane->SetBit(TH1::kCanRebin);

  TH2F* hXvsPlane=new TH2F("hXvsPlane","hXvsPlane",301,-1,300,100,-5,5);
  hXvsPlane->SetTitle("X vs Z");
  hXvsPlane->GetXaxis()->SetTitle("Z");
  hXvsPlane->GetXaxis()->CenterTitle();
  hXvsPlane->GetYaxis()->SetTitle("Y");
  hXvsPlane->GetYaxis()->CenterTitle();
  hXvsPlane->SetFillColor(0);
  //hXvsPlane->SetBit(TH1::kCanRebin);

  vector<TH1F*> vhStrips(282);
  vector<TH1F*> vhTPos(282);
  for (Int_t i=0;i<282;i++) {
    string sName="hStripsPlane";
    string sPl=Form("%d",i);
    sName+=sPl;
    //cout<<"name="<<sName<<endl;
    TH1F* hStrip=new TH1F(sName.c_str(),sName.c_str(),
                          101,-1,100);
    hStrip->SetTitle(sName.c_str());
    hStrip->GetXaxis()->SetTitle("Strip");
    hStrip->GetXaxis()->CenterTitle();
    hStrip->GetYaxis()->SetTitle("");
    hStrip->GetYaxis()->CenterTitle();
    hStrip->SetFillColor(0);
    //hStrip->SetBit(TH1::kCanRebin);

    vhStrips[i]=hStrip;
    //strips are 0-95(96), 4-67(64), 0-63(64)
      
    sName="hTPosPlane";
    sName+=sPl;
    TH1F* hTPos=new TH1F(sName.c_str(),sName.c_str(),200,-4,4);
    hTPos->SetTitle(sName.c_str());
    hTPos->GetXaxis()->SetTitle("TPos");
    hTPos->GetXaxis()->CenterTitle();
    hTPos->GetYaxis()->SetTitle("");
    hTPos->GetYaxis()->CenterTitle();
    hTPos->SetFillColor(0);
    //hTPos->SetBit(TH1::kCanRebin);

    vhTPos[i]=hTPos;
  }

  vector<TH1F*> vhStripsTrk(282);
  vector<TH1F*> vhTPosTrk(282);
  for (Int_t i=0;i<282;i++) {
    string sName="hStripsTrkPlane";
    string sPl=Form("%d",i);
    sName+=sPl;
    //cout<<"name="<<sName<<endl;
    TH1F* hStrip=new TH1F(sName.c_str(),sName.c_str(),
                          101,-1,100);
    hStrip->SetTitle(sName.c_str());
    hStrip->GetXaxis()->SetTitle("Strip");
    hStrip->GetXaxis()->CenterTitle();
    hStrip->GetYaxis()->SetTitle("");
    hStrip->GetYaxis()->CenterTitle();
    hStrip->SetFillColor(0);
    //hStrip->SetBit(TH1::kCanRebin);

    vhStripsTrk[i]=hStrip;
    //strips are 0-95(96), pl 8,10: 4-67(64), pl 7,9: 0-63(64)
    //plane 6 has low strips in PI part (I think 0-67)
    //plane 11 has high strips in PI part (I think -95)

    //strips TPos:
    //plane 6 goes -2.64 -> 1.32 m
    //plane 11 goes -1.32 -> 2.64 m

    //plane 4,8,10 goes -2.40 -> 0.24 m
    //plane 5,7,9 goes -0.24 -> 2.40 m

    //2.64 - 2.40 = 24 cm = 5.85 strips
    //the FI planes "stick out" by ~6 strips

    sName="hTPosTrkPlane";
    sName+=sPl;
    TH1F* hTPos=new TH1F(sName.c_str(),sName.c_str(),200,-4,4);
    hTPos->SetTitle(sName.c_str());
    hTPos->GetXaxis()->SetTitle("TPos");
    hTPos->GetXaxis()->CenterTitle();
    hTPos->GetYaxis()->SetTitle("");
    hTPos->GetYaxis()->CenterTitle();
    hTPos->SetFillColor(0);
    //hTPos->SetBit(TH1::kCanRebin);

    vhTPosTrk[i]=hTPos;
  }

  //reconstruction object
  MeuReco reco;
  MeuSummary s;

  
  this->InitialiseLoopVariables();  
  
  //for(Int_t entry=0;entry<10000;entry++){
  for(Int_t entry=0;entry<fEntries;entry++){
      
    this->SetLoopVariables(entry);
      
    NtpStRecord& ntpst=(*fRec);//Make a reference instead of a pointer
    
    //
    this->FilterLI(ntpst,s);
    continue;

    TClonesArray& tcaTk=(*ntpst.trk);
    Int_t numTrks=tcaTk.GetEntries();
    
    MAXMSG("MeuAnalysis",Msg::kInfo,10)
      <<"numTrks="<<numTrks
      <<", ndigit="<<fRec->evthdr.ndigit
      <<", nTrack="<<fRec->evthdr.ntrack<<endl;

    MeuSummary s;
    s.RFid=1.0;
    map<Int_t,MeuHitInfo> cp;
    
    hNTrack->Fill(numTrks);
    if (numTrks>1) continue;

    //Loop over tracks
    for (Int_t itrk=0;itrk<numTrks;itrk++){
      const NtpSRTrack* ptrk=
        dynamic_cast<NtpSRTrack*>(tcaTk[itrk]);
      const NtpSRTrack& trk=(*ptrk);

      MAXMSG("MeuAnalysis",Msg::kInfo,10)
        <<"range="<<trk.range<<" nstrip="<<trk.nstrip<<" ds="<<trk.ds
        <<" fidend.dr,dz="<<trk.fidend.dr<<" "<<trk.fidend.dz<<" end-vtx.z:"
        <<trk.end.z-trk.vtx.z<<" "<<trk.stpz[1]-trk.stpz[0]<<endl;
      
      hStPerTrk->Fill(trk.nstrip);
      hYvsXVtx->Fill(trk.vtx.x,trk.vtx.y);
      hYvsXEnd->Fill(trk.end.x,trk.end.y);
      if (trk.end.plane>5 && trk.end.plane<115){
        hYvsXVtxFid->Fill(trk.vtx.x,trk.vtx.y);
        hYvsXEndFid->Fill(trk.end.x,trk.end.y);
      }

      TClonesArray& tcaStp=(*ntpst.stp);
      Int_t numStps=tcaStp.GetEntries();
      MAXMSG("MeuAnalysis",Msg::kDebug,10)
        <<"numStps="<<numStps<<endl;
      
      Float_t minR=999;
      Float_t maxR=0;

      for (Int_t i=0;i<trk.nstrip;i++){

        //check for bug where strip index is -1
        if (trk.stp[i]<0) {
          MAXMSG("MeuAnalysis",Msg::kInfo,500)
            <<"Skipping strip with trk.stp[i]="<<trk.stp[i]<<endl;
          continue;
        }

        //const NtpSRStrip* pstp=dynamic_cast<NtpSRStrip*>(tcaStp[i]);
        const NtpSRStrip* pstp=
          dynamic_cast<NtpSRStrip*>(tcaStp[trk.stp[i]]);

        const NtpSRStrip& stp=(*pstp);
        MAXMSG("MeuAnalysis",Msg::kInfo,20)
          <<"Strip="<<stp.strip<<", pl="<<stp.plane
          <<", sigCor="<<stp.ph1.sigcor<<", tpos="<<stp.tpos<<endl;
        
        vhStripsTrk[stp.plane]->Fill(stp.strip);
        vhTPosTrk[stp.plane]->Fill(stp.tpos);
 
        if (stp.plane>120) continue;

        hTPos->Fill(stp.tpos);
        hSigCor->Fill(stp.ph1.sigcor);
        hView->Fill(stp.planeview);
        hPlane->Fill(stp.plane);
        hStrip->Fill(stp.strip);
        hZ->Fill(stp.z);

        //fill the calibpos
        MeuHitInfo& c=cp[stp.plane];
        c.Plane=stp.plane;
        c.View=stp.planeview;
        c.Strip=stp.strip;
        c.TPos=stp.tpos;
        c.SigCor2=stp.ph1.sigcor;
        c.X=trk.stpx[i];
        c.Y=trk.stpy[i];
        c.Z=trk.stpz[i];
        
        hTPosVsPlane->Fill(stp.plane,stp.tpos);
        MAXMSG("MeuAnalysis",Msg::kInfo,10)
          <<"x="<<trk.stpx[i]
          <<", y="<<trk.stpy[i]
          <<", z="<<trk.stpz[i]<<endl;
        hYvsX->Fill(trk.stpx[i],trk.stpy[i]);
        hYvsZ->Fill(trk.stpz[i],trk.stpy[i]);
        hXvsZ->Fill(trk.stpz[i],trk.stpx[i]);
        hYvsPlane->Fill(stp.plane,trk.stpy[i]);
        hXvsPlane->Fill(stp.plane,trk.stpx[i]);
        
        Float_t r=sqrt(pow((trk.stpx[i]-1.5),2)+pow(trk.stpy[i],2));
        if (r<minR) minR=r;
        if (r>maxR) maxR=r;
        hR->Fill(r);

      }

      //fill r histos
      if (maxR!=0 && minR!=999) {
        hMaxR->Fill(maxR);
        hMinR->Fill(minR);
      }
    }

    TClonesArray& tcaStp=(*ntpst.stp);
    Int_t numStps=tcaStp.GetEntries();
    for (Int_t i=0;i<numStps;i++){      
      const NtpSRStrip* pstp=dynamic_cast<NtpSRStrip*>(tcaStp[i]);
      const NtpSRStrip& stp=(*pstp);
      //MAXMSG("MeuAnalysis",Msg::kInfo,200)
      //  <<"Strip="<<stp.strip<<", pl="<<stp.plane<<endl;
        
      vhStrips[stp.plane]->Fill(stp.strip);
      vhTPos[stp.plane]->Fill(stp.tpos);
    }

  }//end of for                                       
  

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

  MSG("MeuAnalysis",Msg::kInfo) 
    <<" ** Finished BasicPlots method **"<<endl;
}

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

Bool_t MeuAnalysis::PlaneIsWithinTrack(const MeuSummary& s,
                                        Int_t plane,Int_t vtxPl,
                                        Int_t endPl) const
{
  Bool_t withinTrk=false;
  
  if (s.Detector==Detector::kNear){
    if (plane<0 || plane>120) {
      MSG("MeuAnalysis",Msg::kVerbose)
        <<"MeuAnalysis::PlaneIsWithinTrack."
        <<" Plane is out of range="<<plane
        <<", det="<<s.Detector<<endl;
      return withinTrk;
    }
  }
  else if (s.Detector==Detector::kFar){
    if (plane<0 || plane>485) {
      MSG("MeuAnalysis",Msg::kVerbose)
        <<"MeuAnalysis::PlaneIsWithinTrack."
        <<" Plane is out of range="<<plane
        <<", det="<<s.Detector<<endl;
      return withinTrk;
    }
  }
  else {
    MSG("MeuAnalysis",Msg::kWarning)
      <<"MeuAnalysis::PlaneIsWithinTrack."
      <<" Detector not supported"<<plane
      <<endl; 
  }


  //count the end planes as being "within" the track
  if ((plane>=vtxPl && plane<=endPl) ||
      (plane<=vtxPl && plane>=endPl)) withinTrk=true;
  
  return withinTrk;
}

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

void MeuAnalysis::CheckTrackGaps2(const NtpStRecord& ntp) const
{
  TClonesArray& tcaTk=(*ntp.trk);
  Int_t numTrks=tcaTk.GetEntries();
  
  //Loop over tracks
  for (Int_t itrk=0;itrk<numTrks;itrk++){
    const NtpSRTrack* ptrk=
      dynamic_cast<NtpSRTrack*>(tcaTk[itrk]);
    const NtpSRTrack& trk=(*ptrk);
      
    TClonesArray& tcaStp=(*ntp.stp);
      
    //this loop is just over the tracked strips
    for (Int_t i=0;i<trk.nstrip;i++){
        
      //check for bug where strip index is -1
      if (trk.stp[i]<0) {
        MAXMSG("MeuAnalysis",Msg::kInfo,500)
          <<"Skipping strip with trk.stp[i]="<<trk.stp[i]<<endl;
        continue;
      }

      const NtpSRStrip* pstp=
        dynamic_cast<NtpSRStrip*>(tcaStp[trk.stp[i]]);
      const NtpSRStrip& stp=(*pstp);
      //MAXMSG("MeuAnalysis",Msg::kInfo,200)
      //  <<"Strip="<<stp.strip<<", pl="<<stp.plane
      //  <<", sigCor="<<stp.ph1.sigcor<<", tpos="<<stp.tpos<<endl;

      Int_t pl=stp.plane;

      if (pl>120) continue;
    }
  }
}

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

void MeuAnalysis::Recalibrate(const NtpStRecord& ntp,MeuSummary& s,
                              std::map<Int_t,MeuHitInfo>& plInfo,
                              const NtpSREvent& evt) const
{
  
  //variable to turn on all the useful messages if required
  Msg::LogLevel_t logLevel=Msg::kDebug;
  //Msg::LogLevel_t logLevelSts=Msg::kDebug;
  //Msg::LogLevel_t logLevelLin=Msg::kInfo;

  //make plots of the calibration constants for the recalibration
  //this is also done elsewhere for ntuple values of the cal constants
  //if the recalibration is not done then they are the same as sntp
  static TH1F* hAdcPerSigLinReCal1=0;
  static TH1F* hSigLinPerSigCorReCal1=0;
  static TH1F* hSigCorPerSigMapReCal1=0;
  static TH1F* hSigMapPerMipReCal1=0;
  static TH1F* hMipPerGeVReCal1=0;
  static TH1F* hAdcPerPeReCal1=0;
  static TH1F* hAdcPerSigDrfReCal1=0;
  static TH1F* hAdcPerSigLinOnlyReCal1=0;

  static TH1F* hAdcPerSigLinReCal2=0;
  static TH1F* hSigLinPerSigCorReCal2=0;
  static TH1F* hSigCorPerSigMapReCal2=0;
  static TH1F* hSigMapPerMipReCal2=0;
  static TH1F* hMipPerGeVReCal2=0;
  static TH1F* hAdcPerPeReCal2=0;
  static TH1F* hAdcPerSigDrfReCal2=0;
  static TH1F* hAdcPerSigLinOnlyReCal2=0;

  if (hAdcPerSigLinReCal1==0) {
    //stripend 1
    hAdcPerSigLinReCal1=new TH1F("hAdcPerSigLinReCal1",
                                 "hAdcPerSigLinReCal1",
                                 1100,-1,10);
    hAdcPerSigLinReCal1->SetFillColor(0);
    hAdcPerSigLinReCal1->SetTitle("AdcPerSigLin");
    hAdcPerSigLinReCal1->GetXaxis()->SetTitle("AdcPerSigLin");
    hAdcPerSigLinReCal1->GetXaxis()->CenterTitle();
    //hAdcPerSigLinReCal1->SetBit(TH1::kCanRebin);
    
    hSigLinPerSigCorReCal1=new TH1F("hSigLinPerSigCorReCal1",
                                 "hSigLinPerSigCorReCal1",
                                 1100,-1,10);
    hSigLinPerSigCorReCal1->SetFillColor(0);
    hSigLinPerSigCorReCal1->SetTitle("SigLinPerSigCor");
    hSigLinPerSigCorReCal1->GetXaxis()->SetTitle("SigLinPerSigCor");
    hSigLinPerSigCorReCal1->GetXaxis()->CenterTitle();
    //hSigLinPerSigCorReCal1->SetBit(TH1::kCanRebin);
    
    hSigCorPerSigMapReCal1=new TH1F("hSigCorPerSigMapReCal1",
                                 "hSigCorPerSigMapReCal1",
                                 10000,-1,100);
    hSigCorPerSigMapReCal1->SetFillColor(0);
    hSigCorPerSigMapReCal1->SetTitle("SigCorPerSigMap");
    hSigCorPerSigMapReCal1->GetXaxis()->SetTitle("SigCorPerSigMap");
    hSigCorPerSigMapReCal1->GetXaxis()->CenterTitle();
    //hSigCorPerSigMapReCal1->SetBit(TH1::kCanRebin);
    
    hSigMapPerMipReCal1=new TH1F("hSigMapPerMipReCal1",
                              "hSigMapPerMipReCal1",
                              10100,-10,1000);
    hSigMapPerMipReCal1->SetFillColor(0);
    hSigMapPerMipReCal1->SetTitle("SigMapPerMip");
    hSigMapPerMipReCal1->GetXaxis()->SetTitle("SigMapPerMip");
    hSigMapPerMipReCal1->GetXaxis()->CenterTitle();
    //hSigMapPerMipReCal1->SetBit(TH1::kCanRebin);
    
    hMipPerGeVReCal1=new TH1F("hMipPerGeVReCal1","hMipPerGeVReCal1",
                           10100,-10,1000);
    hMipPerGeVReCal1->SetFillColor(0);
    hMipPerGeVReCal1->SetTitle("MipPerGeV");
    hMipPerGeVReCal1->GetXaxis()->SetTitle("MipPerGeV");
    hMipPerGeVReCal1->GetXaxis()->CenterTitle();
    //hMipPerGeVReCal1->SetBit(TH1::kCanRebin);

    hAdcPerPeReCal1=new TH1F("hAdcPerPeReCal1","hAdcPerPeReCal1",
                           8080,-20,2000);
    hAdcPerPeReCal1->SetFillColor(0);
    hAdcPerPeReCal1->SetTitle("AdcPerPe");
    hAdcPerPeReCal1->GetXaxis()->SetTitle("AdcPerPe");
    hAdcPerPeReCal1->GetXaxis()->CenterTitle();
    //hAdcPerPeReCal1->SetBit(TH1::kCanRebin);

    hAdcPerSigDrfReCal1=new TH1F("hAdcPerSigDrfReCal1",
                                 "hAdcPerSigDrfReCal1",
                                 1100,-1,10);
    hAdcPerSigDrfReCal1->SetFillColor(0);
    hAdcPerSigDrfReCal1->SetTitle("AdcPerSigDrf");
    hAdcPerSigDrfReCal1->GetXaxis()->SetTitle("AdcPerSigDrf");
    hAdcPerSigDrfReCal1->GetXaxis()->CenterTitle();
    //hAdcPerSigDrfReCal1->SetBit(TH1::kCanRebin);

    hAdcPerSigLinOnlyReCal1=new TH1F("hAdcPerSigLinOnlyReCal1",
                                     "hAdcPerSigLinOnlyReCal1",
                                     2100,-1,20);
    hAdcPerSigLinOnlyReCal1->SetFillColor(0);
    hAdcPerSigLinOnlyReCal1->SetTitle("AdcPerSigLinOnly");
    hAdcPerSigLinOnlyReCal1->GetXaxis()->SetTitle("AdcPerSigLinOnly");
    hAdcPerSigLinOnlyReCal1->GetXaxis()->CenterTitle();
    //hAdcPerSigLinOnlyReCal1->SetBit(TH1::kCanRebin);



    //stripend 2
    hAdcPerSigLinReCal2=new TH1F("hAdcPerSigLinReCal2","hAdcPerSigLinReCal2",
                              1100,-1,10);
    hAdcPerSigLinReCal2->SetFillColor(0);
    hAdcPerSigLinReCal2->SetTitle("AdcPerSigLin");
    hAdcPerSigLinReCal2->GetXaxis()->SetTitle("AdcPerSigLin");
    hAdcPerSigLinReCal2->GetXaxis()->CenterTitle();
    //hAdcPerSigLinReCal2->SetBit(TH1::kCanRebin);
    
    hSigLinPerSigCorReCal2=new TH1F("hSigLinPerSigCorReCal2",
                                 "hSigLinPerSigCorReCal2",
                                 1100,-1,10);
    hSigLinPerSigCorReCal2->SetFillColor(0);
    hSigLinPerSigCorReCal2->SetTitle("SigLinPerSigCor");
    hSigLinPerSigCorReCal2->GetXaxis()->SetTitle("SigLinPerSigCor");
    hSigLinPerSigCorReCal2->GetXaxis()->CenterTitle();
    //hSigLinPerSigCorReCal2->SetBit(TH1::kCanRebin);
    
    hSigCorPerSigMapReCal2=new TH1F("hSigCorPerSigMapReCal2",
                                 "hSigCorPerSigMapReCal2",
                                 10000,-1,100);
    hSigCorPerSigMapReCal2->SetFillColor(0);
    hSigCorPerSigMapReCal2->SetTitle("SigCorPerSigMap");
    hSigCorPerSigMapReCal2->GetXaxis()->SetTitle("SigCorPerSigMap");
    hSigCorPerSigMapReCal2->GetXaxis()->CenterTitle();
    //hSigCorPerSigMapReCal2->SetBit(TH1::kCanRebin);
    
    hSigMapPerMipReCal2=new TH1F("hSigMapPerMipReCal2",
                              "hSigMapPerMipReCal2",
                              10100,-10,1000);
    hSigMapPerMipReCal2->SetFillColor(0);
    hSigMapPerMipReCal2->SetTitle("SigMapPerMip");
    hSigMapPerMipReCal2->GetXaxis()->SetTitle("SigMapPerMip");
    hSigMapPerMipReCal2->GetXaxis()->CenterTitle();
    //hSigMapPerMipReCal2->SetBit(TH1::kCanRebin);
    
    hMipPerGeVReCal2=new TH1F("hMipPerGeVReCal2","hMipPerGeVReCal2",
                           10100,-10,1000);
    hMipPerGeVReCal2->SetFillColor(0);
    hMipPerGeVReCal2->SetTitle("MipPerGeV");
    hMipPerGeVReCal2->GetXaxis()->SetTitle("MipPerGeV");
    hMipPerGeVReCal2->GetXaxis()->CenterTitle();
    //hMipPerGeVReCal2->SetBit(TH1::kCanRebin);

    hAdcPerPeReCal2=new TH1F("hAdcPerPeReCal2","hAdcPerPeReCal2",
                           8080,-20,2000);
    hAdcPerPeReCal2->SetFillColor(0);
    hAdcPerPeReCal2->SetTitle("AdcPerPe");
    hAdcPerPeReCal2->GetXaxis()->SetTitle("AdcPerPe");
    hAdcPerPeReCal2->GetXaxis()->CenterTitle();
    //hAdcPerPeReCal2->SetBit(TH1::kCanRebin);

    hAdcPerSigDrfReCal2=new TH1F("hAdcPerSigDrfReCal2",
                                 "hAdcPerSigDrfReCal2",
                                 1100,-1,10);
    hAdcPerSigDrfReCal2->SetFillColor(0);
    hAdcPerSigDrfReCal2->SetTitle("AdcPerSigLinDrf");
    hAdcPerSigDrfReCal2->GetXaxis()->SetTitle("AdcPerSigLinDrf");
    hAdcPerSigDrfReCal2->GetXaxis()->CenterTitle();
    //hAdcPerSigDrfReCal2->SetBit(TH1::kCanRebin);
    
    hAdcPerSigLinOnlyReCal2=new TH1F("hAdcPerSigLinOnlyReCal2",
                                     "hAdcPerSigLinOnlyReCal2",
                                     2100,-1,20);
    hAdcPerSigLinOnlyReCal2->SetFillColor(0);
    hAdcPerSigLinOnlyReCal2->SetTitle("AdcPerSigLinOnly");
    hAdcPerSigLinOnlyReCal2->GetXaxis()->SetTitle("AdcPerSigLinOnly");
    hAdcPerSigLinOnlyReCal2->GetXaxis()->CenterTitle();
    //hAdcPerSigLinOnlyReCal2->SetBit(TH1::kCanRebin);
  }

  Calibrator& cal=Calibrator::Instance();
  //VldTimeStamp (const time_t &t, const Int_t nsec)
  //VldTimeStamp vldts;
  VldTimeStamp vldts(s.TimeSec,0);
  VldContext vc(static_cast<Detector::Detector_t>(s.Detector),
                static_cast<SimFlag::ESimFlag>(s.SimFlag),vldts);
  cal.ReInitialise(vc);//same as reset
  //vc.Print();
  //cal.PrintConfig(cout);

  MAXMSG("MeuAnalysis",logLevel,1000)
    <<"Recalibrate: track "<<s.VtxPlane<<"->"<<s.EndPlane<<endl; 
  
  //zero all the current values
  for (map<Int_t,MeuHitInfo>::iterator it=plInfo.begin();
       it!=plInfo.end();++it){
    //cache the current cp
    MeuHitInfo& cp=it->second;
    
    MAXMSG("MeuAnalysis",logLevel,1000)
      <<"pl="<<cp.Plane      
      <<", st="<<cp.Strip
      <<", sigCor="<<cp.SigCor
      <<", sigMap="<<cp.SigMap
      <<", lpos="<<cp.LPos<<endl;

    cp.SigMap=0;
    cp.SigMap1=0;
    cp.SigMap2=0;
    cp.MCSigMap=0;

    if (fRecalibrateSigCor) {
      MAXMSG("MeuAnalysis",Msg::kInfo,1)
        <<"Reapplying strip-to-strip calibration"
        <<", s.SimFlag="<<s.SimFlag<<endl;
      cp.SigCor=0;
      cp.SigCor1=0;
      cp.SigCor2=0;
      cp.SigCorTrk1=0;
      cp.SigCorTrk2=0;
    }
    else {
      MAXMSG("MeuAnalysis",Msg::kInfo,1)
        <<"NOT Reapplying strip-to-strip calibration"<<endl;
    }

    if (fRecalibrateSigLin) {
      MAXMSG("MeuAnalysis",Msg::kInfo,1)
        <<"Reapplying SigLin calibration"<<endl;
      cp.SigLin=0;
      cp.SigLin1=0;
      cp.SigLin2=0;
    }
    else {
      MAXMSG("MeuAnalysis",Msg::kInfo,1)
        <<"NOT Reapplying SigLin calibration (hence no drift)"<<endl;
    }

    if (fRecalibratePE) {
      MAXMSG("MeuAnalysis",Msg::kInfo,1)
        <<"Reapplying PE calibration"<<endl;
      cp.Pe=0;
      cp.Pe1=0;
      cp.Pe2=0;
    }
    else {
      MAXMSG("MeuAnalysis",Msg::kInfo,1)
        <<"NOT Reapplying PE calibration"<<endl;
    }

    //these should already be zero... but to be on the safe side:
    cp.SigLinOnly=0;
    cp.SigLinOnly1=0;
    cp.SigLinOnly2=0;
    cp.SigDrf=0;
    cp.SigDrf1=0;
    cp.SigDrf2=0;
  }
  
  TClonesArray& slcTca=(*ntp.slc);
  TClonesArray& stpTca=(*ntp.stp);

  //get the slice associated with this event
  const NtpSRSlice* pslc=
    dynamic_cast<NtpSRSlice*>(slcTca[evt.slc]);
  const NtpSRSlice& slc=(*pslc);
  
  //loop over strips in slc
  for (Int_t i=0;i<slc.nstrip;++i) {
    const NtpSRStrip* pstp=
      dynamic_cast<NtpSRStrip*>(stpTca[slc.stp[i]]);
    const NtpSRStrip& stp=(*pstp);
    
    Int_t pl=stp.plane;
    //check if the plane is in the track
    //this cut is not in NuCuts::ExtractPlInfo but it should
    //never matter since the track window is always within the track
    if (!this->PlaneIsWithinTrack(s,pl,
                                  s.VtxPlane,s.EndPlane)) continue;

    //filter out the times that are way off to make the 
    //Near and Far detectors more similar
    //get the time0, if between 0 and 1 switch to time1
    Double_t time=stp.time0;
    if (time<0 || time>1) time=stp.time1;
    if (time<0 || time>1) cout<<"Ahhhhh, time"<<endl;;
    
    if (time>s.MedianTime+(400*Munits::ns) ||
        time<s.MedianTime-(400*Munits::ns)) {
      MAXMSG("MeuCuts",Msg::kInfo,10)
        <<"Cut strip: time="<<time
        <<", pl="<<stp.plane<<", st="<<stp.strip
        <<", sigCor="<<stp.ph0.sigcor+stp.ph1.sigcor
        <<", dT="<<(time-s.MedianTime)/Munits::ns<<endl;
      continue;//messes up beg/end planes (if not filtered bad trks)!
    }

    //cache the current cp
    MeuHitInfo& cp=plInfo[pl];

    //get the raw adcs
    Float_t adc1=stp.ph0.raw;//east
    Float_t adc2=stp.ph1.raw;//west

    PlexStripEndId seidE(static_cast<Detector::Detector_t>
                         (s.Detector),pl,stp.strip,
                         StripEnd::kEast);
    PlexStripEndId seidW(static_cast<Detector::Detector_t>
                         (s.Detector),pl,stp.strip,
                         StripEnd::kWest);
    
    //sigDrf and sigLinOnly values are NOT used to get the MEU number
    //get the sigLinOnly value (applied to strip)
    //useful for systematic error evaluation
    //NOTE: in macro have to set:
    //cal.LinearityCalScheme().Set("BucketMode=0"); 
    //for the results to have any meaning
    Float_t sigLinOnly1=0;
    Float_t sigLinOnly2=0;
    if (adc1>0) sigLinOnly1=cal.GetLinearized(adc1,seidE);
    else sigLinOnly1=0;
    if (adc2>0) sigLinOnly2=cal.GetLinearized(adc2,seidW);
    else sigLinOnly2=0;

    //now sum up everything in the plane
    cp.SigLinOnly+=sigLinOnly1+sigLinOnly2;
    cp.SigLinOnly1+=sigLinOnly1;
    cp.SigLinOnly2+=sigLinOnly2;
    
    //get the drift correction applied to adcs
    Float_t sigDrf1=0;
    Float_t sigDrf2=0;
    if (adc1>0) sigDrf1=cal.GetDriftCorrected(adc1,seidE);
    else sigDrf1=0;
    if (adc2>0) sigDrf2=cal.GetDriftCorrected(adc2,seidW);
    else sigDrf2=0;

    //now sum up everything in the plane
    cp.SigDrf+=sigDrf1+sigDrf2;
    cp.SigDrf1+=sigDrf1;
    cp.SigDrf2+=sigDrf2;


    //section to recalibrate PEs if required
    Float_t pe1=stp.ph0.pe;//east
    Float_t pe2=stp.ph1.pe;//west
    if (fRecalibratePE){
      if (adc1>0) pe1=cal.GetPhotoElectrons(adc1,seidE);
      else pe1=0;
      if (adc2>0) pe2=cal.GetPhotoElectrons(adc2,seidW);
      else pe2=0;

      //now sum up everything in the plane
      cp.Pe+=pe1+pe2;
      cp.Pe1+=pe1;
      cp.Pe2+=pe2;
    }
    
    //section to recalibrate siglin if required
    Float_t sigLin1=stp.ph0.siglin;//east
    Float_t sigLin2=stp.ph1.siglin;//west
    if (fRecalibrateSigLin){
      //linearity calibration is done on digits (not strips)
      //this matters for the ND!
      //It is not totally correct to recalibrate with the standard
      //ntuples because the digit information is not available for ND
      //however the calibrator can implement an accurate fudge:
      //cal.LinearityCalScheme().Set("BucketMode=0"); 

      MAXMSG("MeuAnalysis",Msg::kInfo,1)
        <<"Recalibrating drift and linearity for detector="
        <<s.Detector<<endl;
      
      //first apply the linearity 
      if (adc1>0) sigLin1=cal.GetLinearized(adc1,seidE);
      else sigLin1=0;
      if (adc2>0) sigLin2=cal.GetLinearized(adc2,seidW);
      else sigLin2=0;
      
      //then do the drift correction (to sigLin - driftCor this time)
      if (sigLin1>0) sigLin1=cal.GetDriftCorrected(sigLin1,seidE);
      else sigLin1=0;
      if (sigLin2>0) sigLin2=cal.GetDriftCorrected(sigLin2,seidW);
      else sigLin2=0;   
      
      //now sum up everything in the plane
      cp.SigLin+=sigLin1+sigLin2;
      cp.SigLin1+=sigLin1;
      cp.SigLin2+=sigLin2;
    }
    
    //section to recalibrate sigcor if required
    Float_t sigCor1=stp.ph0.sigcor;//east
    Float_t sigCor2=stp.ph1.sigcor;//west
    if (fRecalibrateSigCor){
      MAXMSG("MeuAnalysis",Msg::kInfo,1)
        <<"Recalibrating strip-to-strip..."<<endl;
      if (s.Detector!=Detector::kNear) {
        if (sigLin1>0) sigCor1=cal.GetStripToStripCorrected(sigLin1,
                                                            seidE);
      }
      //else sigCor1=0;
      if (sigLin2>0) sigCor2=cal.GetStripToStripCorrected(sigLin2,
                                                          seidW);
      else sigCor2=0;
      
      //sum both stripends
      Float_t sigCorBoth=sigCor1+sigCor2;
      
      //sum the sigcor in the object
      cp.SigCor+=sigCorBoth;
      cp.SigCor1+=sigCor1;
      cp.SigCor2+=sigCor2;
      
      //store the sigcor of both ends of the tracked strip
      //this means that fSigCor!=fSigCor1+fSigCor2 always
      if (stp.strip==cp.Strip){
        cp.SigCorTrk1=sigCor1;
        cp.SigCorTrk2=sigCor2;
      }
    }
    else if (fRecalibrateSigLin) {
      MAXMSG("MeuAnalysis",Msg::kInfo,1)
        <<"Using ntp strip-to-strip calibration on recalibrated"
        <<" SigLin values..."<<endl;
      //if siglin was recalibrated then have to feed the value through
      //the value of the sigcor correction in the ntp is just
      //stp.ph0.sigcor/stp.ph0.siglin so multiply by that:
      sigCor1=sigLin1*stp.ph0.sigcor/stp.ph0.siglin;//east
      sigCor2=sigLin2*stp.ph1.sigcor/stp.ph1.siglin;//west
    }

    //find the lpos for this plane then calibrate
    Float_t lpos=cp.LPos;
    Float_t sigMap1=0;
    Float_t sigMap2=0;

    //only calibrate the east side for !ND
    //also only do the calibration if there is a hit on stripend
    if (s.Detector!=Detector::kNear) {
      if (sigCor1>0) sigMap1=cal.GetAttenCorrectedTpos(sigCor1,
                                                       lpos,seidE);
    }
    if (sigCor2>0) sigMap2=cal.GetAttenCorrectedTpos(sigCor2,
                                                     lpos,seidW);
    
    Float_t sigMapBoth=sigMap1+sigMap2;
    
    //sum the sigmap
    cp.SigMap+=sigMapBoth;
    cp.SigMap1+=sigMap1;
    cp.SigMap2+=sigMap2;
    
    MAXMSG("MeuAnalysis",logLevel,1000)
      <<"strip="<<stp.strip<<", pl="<<pl<<", lpos="<<lpos<<endl
      <<"  sigCor1="<<sigCor1<<", sigMap1="<<sigMap1<<endl
      <<"  sigCor2="<<sigCor2<<", sigMap2="<<sigMap2<<endl
      <<endl;

    //calibrate "true MC" values as well
    if (s.SimFlag==SimFlag::kMC){
      //get the MC lpos
      Float_t lpos=cp.MCLPos;

      Float_t sigMap1=0;
      Float_t sigMap2=0;

      if (sigCor1>0) {
        sigMap1=cal.GetAttenCorrectedTpos(sigCor1,lpos,seidE);
      }
      if (sigCor2>0) {
        sigMap2=cal.GetAttenCorrectedTpos(sigCor2,lpos,seidW);
      }

      //calc sigmap and add to cp
      sigMapBoth=sigMap1+sigMap2;
      cp.MCSigMap+=sigMapBoth;
    }

    //fill calibration value histos
    //can't do mip or gev ones
    if (sigLin1) hAdcPerSigLinReCal1->Fill(adc1/sigLin1);
    if (sigCor1) hSigLinPerSigCorReCal1->Fill(sigLin1/sigCor1);
    if (sigCor1) hSigCorPerSigMapReCal1->Fill(sigCor1/sigMap1);
    if (pe1) hAdcPerPeReCal1->Fill(adc1/pe1);
    if (sigLinOnly1) hAdcPerSigLinOnlyReCal1->Fill(adc1/sigLinOnly1);
    if (sigDrf1) hAdcPerSigDrfReCal1->Fill(adc1/sigDrf1);
    //stripend 2
    if (sigLin2) hAdcPerSigLinReCal2->Fill(adc2/sigLin2);
    if (sigCor2) hSigLinPerSigCorReCal2->Fill(sigLin2/sigCor2);
    if (sigCor2) hSigCorPerSigMapReCal2->Fill(sigCor2/sigMap2);
    if (pe2) hAdcPerPeReCal2->Fill(adc2/pe2);
    if (sigLinOnly2) hAdcPerSigLinOnlyReCal2->Fill(adc2/sigLinOnly2);
    if (sigDrf2) hAdcPerSigDrfReCal2->Fill(adc2/sigDrf2);
    
    //print out
    static Bool_t print=false;
    if (print) {
      if (sigDrf2) {
        if (adc2/sigDrf2>2) {
          MAXMSG("MeuAnalysis",Msg::kInfo,1000)
            <<"adc2/sigDrf2="<<adc2/sigDrf2
            <<", adc2="<<adc2<<", sigDrf2="<<sigDrf2<<endl;
        }
      }
      if (sigLinOnly2) {
        if (adc2/sigLinOnly2>2) {
          MAXMSG("MeuAnalysis",Msg::kInfo,1000)
            <<"adc2/sigLinOnly2="<<adc2/sigLinOnly2
            <<", adc2="<<adc2<<", sigLinOnly2="<<sigLinOnly2
            <<", pl,st="<<cp.Plane<<","<<cp.Strip<<endl;
        }
      }
    }
  }//loop over strips in slc

  //print
  Bool_t printMap=false;
  if (printMap){
    for (map<Int_t,MeuHitInfo>::iterator it=plInfo.begin();
         it!=plInfo.end();++it){
      //cache the current cp
      MeuHitInfo& cp=it->second;
      
      MAXMSG("MeuAnalysis",Msg::kInfo,1000)
        <<"Printing: pl="<<cp.Plane      
        <<", st="<<cp.Strip
        <<", sigCor="<<cp.SigCor
        <<", sigMap="<<cp.SigMap
        <<", lpos="<<cp.LPos<<endl;
    }
  }
}

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

void MeuAnalysis::RecalibrateAtNu
(AtmosEvent& ev,MeuSummary& s,
 std::map<Int_t,MeuHitInfo>& plInfo) const
{
  //variable to turn on all the useful messages if required
  Msg::LogLevel_t logLevel=Msg::kDebug;
  //Msg::LogLevel_t logLevelSts=Msg::kDebug;
  //Msg::LogLevel_t logLevelLin=Msg::kInfo;

  Bool_t fRecalibrateSigLin=true;//have to configure the calibrator too! 
  Bool_t fRecalibrateSigCor=true;//have to configure the calibrator too! 
  Bool_t fRecalibratePE=true;//have to configure the calibrator too! 

  Calibrator& cal=Calibrator::Instance();
  //VldTimeStamp (const time_t &t, const Int_t nsec)
  //VldTimeStamp vldts;
  VldTimeStamp vldts(ev.UnixTime,0);
  VldContext vc(static_cast<Detector::Detector_t>(s.Detector),
                static_cast<SimFlag::ESimFlag>(s.SimFlag),vldts);
  cal.ReInitialise(vc);

  static Bool_t firstTime=false;//switch this off with false
  if (firstTime){
    if (s.Detector==Detector::kFar){
      cal.Set("LinCalibrator=PulserLinearityCalScheme");
    }
    else if (s.Detector==Detector::kNear){
      cal.Set("LinCalibrator=QuadLinearityCalScheme"); 
    }
    else cout<<"Ahh, detector="<<s.Detector<<endl;
    
    cout<<"Actually using calibrator configured as:"<<endl;
    cal.PrintConfig(cout);
    firstTime=false;
  }
  
  MAXMSG("MeuAnalysis",logLevel,1000)
    <<"Recalibrate: track "<<s.VtxPlane<<"->"<<s.EndPlane<<endl; 

  //zero all the current values
  for (map<Int_t,MeuHitInfo>::iterator it=plInfo.begin();
       it!=plInfo.end();++it){
    //cache the current cp
    MeuHitInfo& cp=it->second;
    
    MAXMSG("MeuAnalysis",logLevel,1000)
      <<"pl="<<cp.Plane      
      <<", st="<<cp.Strip
      <<", sigCor="<<cp.SigCor
      <<", sigMap="<<cp.SigMap
      <<", lpos="<<cp.LPos<<endl;

    cp.SigMap=0;
    cp.SigMap1=0;
    cp.SigMap2=0;
    cp.MCSigMap=0;
    if (fRecalibrateSigCor) {
      MAXMSG("MeuAnalysis",Msg::kInfo,1)
        <<"Reapplying strip-to-strip calibration, cp.SimFlag="
        <<s.SimFlag<<endl;
      cp.SigCor=0;
      cp.SigCor1=0;
      cp.SigCor2=0;
      cp.SigCorTrk1=0;
      cp.SigCorTrk2=0;
    }
    else {
      MAXMSG("MeuAnalysis",Msg::kInfo,1)
        <<"NOT Reapplying strip-to-strip calibration"<<endl;
    }

    if (fRecalibrateSigLin) {
      MAXMSG("MeuAnalysis",Msg::kInfo,1)
        <<"Reapplying SigLin calibration"<<endl;
      cp.SigLin=0;
      cp.SigLin1=0;
      cp.SigLin2=0;
    }
    else {
      MAXMSG("MeuAnalysis",Msg::kInfo,1)
        <<"NOT Reapplying SigLin calibration"<<endl;
    }

    if (fRecalibratePE) {
      MAXMSG("MeuAnalysis",Msg::kInfo,1)
        <<"Reapplying PE calibration"<<endl;
      cp.Pe=0;
      cp.Pe1=0;
      cp.Pe2=0;
    }
    else {
      MAXMSG("MeuAnalysis",Msg::kInfo,1)
        <<"NOT Reapplying PE calibration"<<endl;
    }

  }

  TClonesArray& strips=(*ev.StripList);
  Int_t numStrips=strips.GetEntries();
  //loop over the strips in the snarl
  for (Int_t hit=0;hit<numStrips;hit++){
    const AtmosStrip* strip=dynamic_cast<AtmosStrip*>(strips[hit]);
    const AtmosStrip& stp=(*strip);
    
    Int_t pl=stp.Plane;
    //check if the plane is in the track
    if (!this->PlaneIsWithinTrack(s,pl,
                                  s.VtxPlane,s.EndPlane)) continue;

    //cache the current cp
    MeuHitInfo& cp=plInfo[pl];

    Float_t adc1=stp.Qadc[0];//east
    Float_t adc2=stp.Qadc[1];//west
    Float_t pe1=stp.QPE[0];//east
    Float_t pe2=stp.QPE[1];//west
    Float_t sigLin1=stp.Qadc[0];//east//siglin is no longer in ntp
    Float_t sigLin2=stp.Qadc[1];//west//so use adcs instead
    Float_t sigCor1=stp.Qadc[0];//east//siglin is no longer in ntp
    Float_t sigCor2=stp.Qadc[1];//west//so use adcs instead
    //Float_t sigCor1=stp.QPEcorr[0];//east//this is pe size so don't use it!
    //Float_t sigCor2=stp.QPEcorr[1];//west
    
    PlexStripEndId seidE(static_cast<Detector::Detector_t>
                         (s.Detector),pl,stp.Strip,
                         StripEnd::kEast);
    PlexStripEndId seidW(static_cast<Detector::Detector_t>
                         (s.Detector),pl,stp.Strip,
                         StripEnd::kWest);

    //section to recalibrate PEs if required
    if (fRecalibratePE){
      if (adc1>0) pe1=cal.GetPhotoElectrons(adc1,seidE);
      else pe1=0;
      if (adc2>0) pe2=cal.GetPhotoElectrons(adc2,seidW);
      else pe2=0;

      //now sum up everything in the plane
      cp.Pe+=pe1+pe2;
      cp.Pe1+=pe1;
      cp.Pe2+=pe2;
    }

    //section to recalibrate siglin if required
    if (fRecalibrateSigLin){
      //first do the drift correction
      if (adc1>0) sigLin1=cal.GetDriftCorrected(adc1,seidE);
      else sigLin1=0;
      if (adc2>0) sigLin2=cal.GetDriftCorrected(adc2,seidW);
      else sigLin2=0;

      //then apply the linearity (to sigLin - driftCor this time)
      if (adc1>0) sigLin1=cal.GetLinearized(sigLin1,seidE);
      else sigLin1=0;
      if (adc2>0) sigLin2=cal.GetLinearized(sigLin2,seidW);
      else sigLin2=0;

      //now sum up everything in the plane
      cp.SigLin+=sigLin1+sigLin2;
      cp.SigLin1+=sigLin1;
      cp.SigLin2+=sigLin2;
    }

    //section to recalibrate sigcor if required
    //  have to do this for AtNu for MC AND DATA 
    //  since it has pe-sized units of "cor"
    if (fRecalibrateSigCor){// && s.SimFlag==SimFlag::kData){
      if (sigLin1>0) sigCor1=cal.GetStripToStripCorrected(sigLin1,
                                                          seidE);
      else sigCor1=0;
      if (sigLin2>0) sigCor2=cal.GetStripToStripCorrected(sigLin2,
                                                          seidW);
      else sigCor2=0;
      
      //sum the sigcor in the object
      cp.SigCor+=sigCor1+sigCor2;
      cp.SigCor1+=sigCor1;
      cp.SigCor2+=sigCor2;
      
      //store the sigcor of both ends of the tracked strip
      //this means that fSigCor!=fSigCor1+fSigCor2 always
      if (stp.Strip==cp.Strip){
        cp.SigCorTrk1=sigCor1;
        cp.SigCorTrk2=sigCor2;
      }
    }

    //find the lpos for this plane then calibrate
    Float_t lpos=cp.LPos;
    Float_t sigMap1=0;
    Float_t sigMap2=0;

    //only do the look up if non zero
    if (sigCor1>0) sigMap1=cal.GetAttenCorrectedTpos(sigCor1,
                                                     lpos,seidE);
    if (sigCor2>0) sigMap2=cal.GetAttenCorrectedTpos(sigCor2,
                                                     lpos,seidW);
    
    Float_t sigMapBoth=sigMap1+sigMap2;
    
    //sum the sigmap
    cp.SigMap+=sigMapBoth;
    //store the sigmap of both ends of the tracked strip
    //this means that fSigMap!=fSigMap1+fSigMap2 always
    //if (st.Strip==cp.Strip){
    cp.SigMap1+=sigMap1;
    cp.SigMap2+=sigMap2;
    //}

    MAXMSG("MeuAnalysis",logLevel,1000)
      <<"strip="<<stp.Strip<<", pl="<<pl<<", lpos="<<lpos<<endl
      <<"  sigCor1="<<sigCor1<<", sigMap1="<<sigMap1<<endl
      <<"  sigCor2="<<sigCor2<<", sigMap2="<<sigMap2<<endl
      <<endl;

    //do MC as well
    if (s.SimFlag==SimFlag::kMC){
      //get the MC lpos
      Float_t lpos=cp.MCLPos;

      Float_t sigMap1=0;
      Float_t sigMap2=0;

      if (sigCor1>0) {
        sigMap1=cal.GetAttenCorrectedTpos(sigCor1,lpos,seidE);
      }
      if (sigCor2>0) {
        sigMap2=cal.GetAttenCorrectedTpos(sigCor2,lpos,seidW);
      }

      //calc sigmap and add to cp
      sigMapBoth=sigMap1+sigMap2;
      cp.MCSigMap+=sigMapBoth;
    }
  }

  //print
  Bool_t printMap=false;
  if (printMap){
    for (map<Int_t,MeuHitInfo>::iterator it=plInfo.begin();
         it!=plInfo.end();++it){
      //cache the current cp
      MeuHitInfo& cp=it->second;
      
      MAXMSG("MeuAnalysis",Msg::kInfo,1000)
        <<"Printing: pl="<<cp.Plane      
        <<", st="<<cp.Strip
        <<", sigCor="<<cp.SigCor
        <<", sigMap="<<cp.SigMap
        <<", lpos="<<cp.LPos<<endl;
    }
  }
}

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

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

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

  MeuSummaryWriter summaryOut;

  TH1F* hWinSigMap=new TH1F("hWinSigMap","hWinSigMap",3000,-1,4000);
  hWinSigMap->SetTitle("hWinSigMap");
  hWinSigMap->GetXaxis()->SetTitle("hWinSigMap");
  hWinSigMap->GetXaxis()->CenterTitle();
  hWinSigMap->SetFillColor(0);
  hWinSigMap->SetLineWidth(3);
  //hWinSigMap->SetBit(TH1::kCanRebin);

  Int_t FCCounter=0;
  Int_t PCCounter=0;
  Int_t TGCounter=0;
  Int_t badWindowCounter=0;
  Int_t badRecoCounter=0;
  Int_t goodWindowCounter=0;
  Int_t badFidCounter=0;

  ofstream& eventsTxtFile=*(this->OpenTxtFile(100,"myevents"));

  //reconstruction object
  MeuReco reco;
  MeuCuts cuts;
  
  //variable to turn on all the useful messages if required
  Msg::LogLevel_t logLevel=Msg::kDebug;

  
  this->InitialiseLoopVariables();  
  
  //for(Int_t entry=0;entry<10000;entry++){
  for(Int_t entry=0;entry<fEntries;entry++){
      
    this->SetLoopVariables(entry);
    
    MSG("MeuAnalysis",logLevel)
      <<"Entry="<<entry<<endl;

    //make a reference instead of a pointer
    NtpStRecord& ntp=(*fRec);
    const RecCandHeader& rec=ntp.GetHeader();
    MAXMSG("MeuAnalysis",Msg::kInfo,1)
      <<"First entry: run="<<rec.GetRun()
      <<", snarl="<<rec.GetSnarl()<<endl;

    //MeuSummary s;
    MeuSummary& s=summaryOut.GetMeuSummaryToFill();
    s.RFid=1.0;
    s.DistToEdgeFid=0.4;
    s.Event=entry;
    cuts.FillSTSumDetails(ntp,s,-1,-1);
    //a map to hold the info for each plane in the entry
    map<Int_t,MeuHitInfo> plInfo;

    //ensure there is only 1 track in the entry
    TClonesArray& trkTca=(*ntp.trk);
    Int_t numTrks=trkTca.GetEntries();
    if (numTrks!=1) continue;
    const NtpSRTrack* ptrk=dynamic_cast<NtpSRTrack*>(trkTca[0]);
    const NtpSRTrack& trk=(*ptrk);
    
    //cut on number of planes
    if (sqrt(pow(1.*trk.plane.beg-trk.plane.end,2))<15) continue;

    MSG("MeuAnalysis",logLevel)
      <<"Entry="<<entry<<", run="<<rec.GetRun()
      <<", snarl="<<rec.GetSnarl()
      <<", trk.nstrips="<<trk.nstrip<<endl;

    TClonesArray& evtTca=(*ntp.evt);
    const NtpSREvent* pevt=dynamic_cast<NtpSREvent*>(evtTca[0]);
    cuts.ExtractPlInfo(ntp,s,plInfo,*pevt);
    
    reco.CalcVtx(s,plInfo);

    //check if the track extends into the calorimeter
    if (s.VtxPlane>120 && s.EndPlane>120) continue;

    //check if track stops in spectrometer
    if (s.EndPlane>120) continue;

    MSG("MeuAnalysis",logLevel)<<"CalcFCPC..."<<endl;
    reco.CalcFCPC(s,plInfo);

    //recalculate the vtx of the tracks coming in from the spectrometer
    MSG("MeuAnalysis",logLevel)<<"CalcVtxSpecialND..."<<endl;
    reco.CalcVtxSpecialND(s,plInfo);

    //cut on number of planes with potential new vertex
    if (sqrt(pow(1.*s.VtxPlane-s.EndPlane,2))<15) continue;

    Bool_t TG=!s.FC && !s.PC;
    
    if (s.PC) PCCounter++;
    else if (s.FC) FCCounter++;
    else if (TG) TGCounter++;
    if (s.FC && s.PC) MSG("MeuAnalysis",Msg::kError)
      <<"Both FC and PC!"<<endl;
    
    if (s.PC) {
      MSG("MeuAnalysis",logLevel)<<"Found PC event..."<<endl;
      Bool_t goodReco=reco.Reconstruct(s,plInfo);
      if (goodReco){
        MSG("MeuAnalysis",logLevel)<<"Recalibrating..."<<endl;
        TClonesArray& evtTca=(*ntp.evt);
        const NtpSREvent* pevt=dynamic_cast<NtpSREvent*>(evtTca[0]);
        this->Recalibrate(ntp,s,plInfo,*pevt);
          
        MSG("MeuAnalysis",logLevel)<<"CalcWindow..."<<endl;
        reco.CalcWindow(s,plInfo);
          
        MSG("MeuAnalysis",logLevel)<<"CheckContainment..."<<endl;
        //fill the tree once event has been analysed
        if (s.WinSigMap<=0) {
          badWindowCounter++;
        }
        else if (!reco.CheckWinContainment(s,plInfo)){
          badFidCounter++;
          //MSG("MeuAnalysis",Msg::kInfo)<<"Outside fid:"<<endl;
        }
        else {
          goodWindowCounter++;
              
          //Bool_t gapsOk=this->CheckTrackGaps(plInfo); 
          //if (!gapsOk) {
          //  MSG("MeuAnalysis",Msg::kInfo)<<"Large gap:"<<endl;
          //  reco.PrintMeuHitInfo(plInfo);
          //}
              
          hWinSigMap->Fill(s.WinSigMap);
              
          eventsTxtFile<<entry<<endl;
          MAXMSG("MeuAnalysis",Msg::kInfo,10)
            <<"Good window: run="
            <<rec.GetRun()<<", snarl="<<rec.GetSnarl()
            <<endl;
              
          cuts.FillSTSumRecoDetails(s,plInfo);
          summaryOut.SummaryTreeFill(plInfo);
        }
      }
      else {
        badRecoCounter++;
        MSG("MeuAnalysis",logLevel)<<"Bad reco..."<<endl;
      }
    }
    
  }//end of for                                       
  

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

  Double_t quantile=0.5;//quantile to compute
  Double_t meu=-1;
  hWinSigMap->GetQuantiles(1,&meu,&quantile);
  Double_t err=hWinSigMap->GetRMS()/sqrt(hWinSigMap->GetEntries());
  Double_t relErr=err/meu;

  MSG("MeuAnalysis",Msg::kInfo) 
    <<endl<<"*** Run Summary ***"<<endl
    <<"Raw MEU="<<meu<<" +/- "<<err<<" ("<<100.*relErr<<"%)"<<endl
    //<<endl

    <<"Total Entries="<<fEntries<<endl
    /*<<"Not LI="<<notLICounter<<endl
      <<"Not MM="<<notMMCounter<<endl
      <<"Good track="<<goodTrackCounter<<endl
      <<"Good length="<<goodLengthCounter<<endl
      <<"Good end & vtx="<<goodEndAndVtxCounter<<endl
      <<"Good XY reco="<<goodXYCounter<<endl
      <<endl
      <<"Good containmentPC="<<containmentCounter<<endl
      <<"  both SM hit (when PC)="<<bothSMCounter
      <<"   ("<<bothSM<<"% of PC)"<<endl
      <<"  coil hit (when PC)   ="<<coilHitCounter
      <<"   ("<<coil_PC<<"% of PC)"<<endl
      <<"  close LI (when PC)   ="<<closeLICounter
      <<"   ("<<closeLI<<"% of PC)"<<endl
      <<"  dead Chips (when PC) ="<<deadChipsCounter
      <<"   ("<<deadChips<<"% of PC)"<<endl
      <<"  boundary TF (when PC)="<<boundaryTFCounter
      <<"   ("<<boundaryTF<<"% of PC)"<<endl
      <<"  bad views (when PC)  ="<<badViewsCounter
      <<"   ("<<badViews<<"% of PC)"<<endl
      <<"  bad trk end (when PC)="<<badTrkEndCounter
      <<"   ("<<badTrkEnd<<"% of PC)"<<endl
      <<"  bad trk gap (when PC)="<<badTrkGapsCounter
      <<"   ("<<badTrkGaps<<"% of PC)"<<endl
      <<"  bad traceZ (when PC) ="<<badTraceZCounter
      <<"   ("<<badTraceZ<<"% of PC)"<<endl
      <<"  bad window (when PC) ="<<badWindowCounter
      <<"   ("<<badWin<<"% of PC)"<<endl
      <<"Good containmentFC="<<containmentCounterFC<<endl
      <<"  both SM hit (when FC)="<<bothSMCounterFC
      <<"   ("<<bothSMFC<<"% of FC)"<<endl
      <<"  coil hit (when FC)   ="<<coilHitCounterFC
      <<"   ("<<coil_FC<<"% of FC)"<<endl
      <<"Good containmentTG="<<containmentCounterTG<<endl
      <<"  both SM hit (when TG)="<<bothSMCounterTG
      <<"   ("<<bothSMTG<<"% of TG)"<<endl
      <<"  coil hit (when TG)   ="<<coilHitCounterTG
      <<"   ("<<coil_TG<<"% of TG)"<<endl
      <<endl
      <<"Bad charge pos moment="<<momentCounter<<endl
      <<"Too high StripsPerPlane="<<stPerPlCounter<<endl
      <<"  Either of above two="<<stPlOrMomentCounter<<endl
      <<endl
      <<"MC Info:"<<endl
      <<"  Muon lowest energy >0.6 GeV="<<lowEnCounter<<endl
      <<"  Muon lowest energy >5.0 GeV="<<lowEn5Counter<<endl
      <<endl
    */
    <<"Containment summary:"<<endl
    <<"  PC="<<PCCounter<<endl
//    <<"     (PC - all cuts="
    //closeLICounter-deadChipsCounter-boundaryTFCounter-
    //badViewsCounter-badTrkEndCounter-badTrkGapsCounter-badTraceZCounter-
    //badWindowCounter<<")"<<endl
    <<"  FC="<<FCCounter<<endl//<<"     (FC-coil-bothSM="
    //<<FCCounter-coilHitCounterFC-bothSMCounterFC<<")"<<endl
    <<"  TG="<<TGCounter<<endl//<<"     (TG-coil-bothSM="
//    <<TGCounter-coilHitCounterTG-bothSMCounterTG<<")"<<endl
    <<"  FC+PC+TG="<<FCCounter+PCCounter+TGCounter<<endl
    <<endl
    <<"Total PC with good track window = "<<goodWindowCounter<<endl
    <<"              bad track window  = "<<badWindowCounter<<endl
    <<"              bad fiducial vol  = "<<badFidCounter<<endl
    <<"              bad reco          = "<<badRecoCounter<<endl
    <<endl;

  //finish the summary tree and output it
  summaryOut.SummaryTreeFinish();

  if (fRec) delete fRec;

  MSG("MeuAnalysis",Msg::kInfo) 
    <<" ** Finished BasicReco method **"<<endl;
}


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

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

  ofstream& snarlList=*(this->OpenTxtFile(100,"snarlList"));
  
  //reconstruction object
  MeuReco reco;
  MeuCuts cuts;
  MeuSummary s;
  
  //variable to turn on all the useful messages if required
  Msg::LogLevel_t logLevel=Msg::kDebug;

  
  this->InitialiseLoopVariables();  
  
  //for(Int_t entry=55000;entry<fEntries;entry++){
  for(Int_t entry=0;entry<fEntries;entry++){
      
    this->SetLoopVariables(entry);
    
    MSG("MeuAnalysis",logLevel)
      <<"Entry="<<entry<<endl;

    //make a reference instead of a pointer
    NtpStRecord& ntp=(*fRec);
    const RecCandHeader& rec=ntp.GetHeader();
    MAXMSG("MeuAnalysis",Msg::kInfo,1)
      <<"First entry: run="<<rec.GetRun()
      <<", snarl="<<rec.GetSnarl()<<endl;

    //MeuSummary s;
    s.Reset();
    s.RFid=1.0;
    s.DistToEdgeFid=0.4;
    s.Event=entry;
    cuts.FillSTSumDetails(ntp,s,-1,-1);
    //a map to hold the info for each plane in the entry
    map<Int_t,MeuHitInfo> plInfo;

    //ensure there is only 1 track in the entry
    TClonesArray& trkTca=(*ntp.trk);
    Int_t numTrks=trkTca.GetEntries();
    if (numTrks!=1) continue;
    const NtpSRTrack* ptrk=dynamic_cast<NtpSRTrack*>(trkTca[0]);
    const NtpSRTrack& trk=(*ptrk);
    
    //cut on number of planes
    Int_t numPlanes=abs(trk.vtx.plane-trk.end.plane);
    if (numPlanes<6) continue;

    TClonesArray& evtTca=(*ntp.evt);
    const NtpSREvent* pevt=dynamic_cast<NtpSREvent*>(evtTca[0]);
    cuts.ExtractPlInfo(ntp,s,plInfo,*pevt);
    
    reco.CalcVtx(s,plInfo);

    //check if track stops in spectrometer
    if (s.EndPlane>120) continue;

    reco.CalcFCPC(s,plInfo);

    //recalculate the vtx of the tracks coming in from the spectrometer
    reco.CalcVtxSpecialND(s,plInfo);

    if (s.PC) {
      MAXMSG("MeuAnalysis",Msg::kInfo,100)
        <<"Entry="<<entry<<", run="<<rec.GetRun()
        <<", snarl="<<rec.GetSnarl()<<endl;
      snarlList<<s.Run<<"\t"<<s.Snarl<<endl;
    }
  }//end of for                                       
  

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

  MSG("MeuAnalysis",Msg::kInfo) 
    <<" ** Finished SnarlList method **"<<endl;
}

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

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

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

  TH1F* hWinSigMap=new TH1F("hWinSigMap","hWinSigMap",3000,-1,4000);
  hWinSigMap->SetTitle("hWinSigMap");
  hWinSigMap->GetXaxis()->SetTitle("hWinSigMap");
  hWinSigMap->GetXaxis()->CenterTitle();
  hWinSigMap->SetFillColor(0);
  hWinSigMap->SetLineWidth(3);
  //hWinSigMap->SetBit(TH1::kCanRebin);

  TH1F* hWinSigCor=new TH1F("hWinSigCor","hWinSigCor",3000,-1,4000);
  hWinSigCor->SetTitle("hWinSigCor");
  hWinSigCor->GetXaxis()->SetTitle("hWinSigCor");
  hWinSigCor->GetXaxis()->CenterTitle();
  hWinSigCor->SetFillColor(0);
  hWinSigCor->SetLineWidth(3);
  //hWinSigCor->SetBit(TH1::kCanRebin);
  
  TH1F* hWinSigLin=new TH1F("hWinSigLin","hWinSigLin",3000,-1,4000);
  hWinSigLin->SetTitle("hWinSigLin");
  hWinSigLin->GetXaxis()->SetTitle("hWinSigLin");
  hWinSigLin->GetXaxis()->CenterTitle();
  hWinSigLin->SetFillColor(0);
  hWinSigLin->SetLineWidth(3);
  //hWinSigLin->SetBit(TH1::kCanRebin);

  TH1F* hWinAdc=new TH1F("hWinAdc","hWinAdc",3000,-1,4000);
  hWinAdc->SetTitle("hWinAdc");
  hWinAdc->GetXaxis()->SetTitle("hWinAdc");
  hWinAdc->GetXaxis()->CenterTitle();
  hWinAdc->SetFillColor(0);
  hWinAdc->SetLineWidth(3);
  //hWinAdc->SetBit(TH1::kCanRebin);

  TH1F* hWinPe=new TH1F("hWinPe","hWinPe",3000,-1,400);
  hWinPe->SetTitle("hWinPe");
  hWinPe->GetXaxis()->SetTitle("hWinPe");
  hWinPe->GetXaxis()->CenterTitle();
  hWinPe->SetFillColor(0);
  hWinPe->SetLineWidth(3);
  //hWinPe->SetBit(TH1::kCanRebin);

  TH1F* hTrkVtxPl=new TH1F("hTrkVtxPl","hTrkVtxPl",501,-1,500);
  hTrkVtxPl->SetTitle("hTrkVtxPl");
  hTrkVtxPl->GetXaxis()->SetTitle("hTrkVtxPl");
  hTrkVtxPl->GetXaxis()->CenterTitle();
  hTrkVtxPl->SetFillColor(0);
  hTrkVtxPl->SetLineWidth(3);
  //hTrkVtxPl->SetBit(TH1::kCanRebin);

  TH1F* hTrkEndPl=new TH1F("hTrkEndPl","hTrkEndPl",501,-1,500);
  hTrkEndPl->SetTitle("hTrkEndPl");
  hTrkEndPl->GetXaxis()->SetTitle("hTrkEndPl");
  hTrkEndPl->GetXaxis()->CenterTitle();
  hTrkEndPl->SetFillColor(0);
  hTrkEndPl->SetLineWidth(3);
  //hTrkEndPl->SetBit(TH1::kCanRebin);

  TH1F* hTrkLengthPl=new TH1F("hTrkLengthPl","hTrkLengthPl",501,-1,500);
  hTrkLengthPl->SetTitle("hTrkLengthPl");
  hTrkLengthPl->GetXaxis()->SetTitle("hTrkLengthPl");
  hTrkLengthPl->GetXaxis()->CenterTitle();
  hTrkLengthPl->SetFillColor(0);
  hTrkLengthPl->SetLineWidth(3);
  //hTrkLengthPl->SetBit(TH1::kCanRebin);

  TH1F* hWinPlFromShwEnd=new TH1F
    ("hWinPlFromShwEnd","hWinPlFromShwEnd",1000,-500,500);
  hWinPlFromShwEnd->SetTitle("hWinPlFromShwEnd");
  hWinPlFromShwEnd->GetXaxis()->SetTitle("hWinPlFromShwEnd");
  hWinPlFromShwEnd->GetXaxis()->CenterTitle();
  hWinPlFromShwEnd->SetFillColor(0);
  hWinPlFromShwEnd->SetLineWidth(3);
  //hWinPlFromShwEnd->SetBit(TH1::kCanRebin);
  
  TH1F* hShwVtxPl=new TH1F("hShwVtxPl","hShwVtxPl",501,-1,500);
  hShwVtxPl->SetTitle("hShwVtxPl");
  hShwVtxPl->GetXaxis()->SetTitle("hShwVtxPl");
  hShwVtxPl->GetXaxis()->CenterTitle();
  hShwVtxPl->SetFillColor(0);
  hShwVtxPl->SetLineWidth(3);
  //hShwVtxPl->SetBit(TH1::kCanRebin);

  TH1F* hShwEndPl=new TH1F("hShwEndPl","hShwEndPl",501,-1,500);
  hShwEndPl->SetTitle("hShwEndPl");
  hShwEndPl->GetXaxis()->SetTitle("hShwEndPl");
  hShwEndPl->GetXaxis()->CenterTitle();
  hShwEndPl->SetFillColor(0);
  hShwEndPl->SetLineWidth(3);
  //hShwEndPl->SetBit(TH1::kCanRebin);

  TH1F* hShwLengthPl=new TH1F("hShwLengthPl","hShwLengthPl",501,-1,500);
  hShwLengthPl->SetTitle("hShwLengthPl");
  hShwLengthPl->GetXaxis()->SetTitle("hShwLengthPl");
  hShwLengthPl->GetXaxis()->CenterTitle();
  hShwLengthPl->SetFillColor(0);
  hShwLengthPl->SetLineWidth(3);
  //hShwLengthPl->SetBit(TH1::kCanRebin);

  TH1F* hTrkShwVtxDiff=new TH1F("hTrkShwVtxDiff","hTrkShwVtxDiff",
                                1000,-500,500);
  hTrkShwVtxDiff->SetTitle("hTrkShwVtxDiff");
  hTrkShwVtxDiff->GetXaxis()->SetTitle("hTrkShwVtxDiff");
  hTrkShwVtxDiff->GetXaxis()->CenterTitle();
  hTrkShwVtxDiff->SetFillColor(0);
  hTrkShwVtxDiff->SetLineWidth(3);
  //hTrkShwVtxDiff->SetBit(TH1::kCanRebin);
  
  TH1F* hTrkShwEndDiff=new TH1F("hTrkShwEndDiff","hTrkShwEndDiff",
                                1000,-500,500);
  hTrkShwEndDiff->SetTitle("hTrkShwEndDiff");
  hTrkShwEndDiff->GetXaxis()->SetTitle("hTrkShwEndDiff");
  hTrkShwEndDiff->GetXaxis()->CenterTitle();
  hTrkShwEndDiff->SetFillColor(0);
  hTrkShwEndDiff->SetLineWidth(3);
  //hTrkShwEndDiff->SetBit(TH1::kCanRebin);
  
  TH1F* hShwVtxPl2=new TH1F("hShwVtxPl2","hShwVtxPl2",501,-1,500);
  hShwVtxPl2->SetTitle("hShwVtxPl2");
  hShwVtxPl2->GetXaxis()->SetTitle("hShwVtxPl2");
  hShwVtxPl2->GetXaxis()->CenterTitle();
  hShwVtxPl2->SetFillColor(0);
  hShwVtxPl2->SetLineWidth(3);
  //hShwVtxPl2->SetBit(TH1::kCanRebin);

  TH1F* hShwEndPl2=new TH1F("hShwEndPl2","hShwEndPl2",501,-1,500);
  hShwEndPl2->SetTitle("hShwEndPl2");
  hShwEndPl2->GetXaxis()->SetTitle("hShwEndPl2");
  hShwEndPl2->GetXaxis()->CenterTitle();
  hShwEndPl2->SetFillColor(0);
  hShwEndPl2->SetLineWidth(3);
  //hShwEndPl2->SetBit(TH1::kCanRebin);

  TH1F* hShwLengthPl2=new TH1F("hShwLengthPl2","hShwLengthPl2",501,-1,500);
  hShwLengthPl2->SetTitle("hShwLengthPl2");
  hShwLengthPl2->GetXaxis()->SetTitle("hShwLengthPl2");
  hShwLengthPl2->GetXaxis()->CenterTitle();
  hShwLengthPl2->SetFillColor(0);
  hShwLengthPl2->SetLineWidth(3);
  //hShwLengthPl2->SetBit(TH1::kCanRebin);

  TH1F* hTrkShwVtxDiff2=new TH1F("hTrkShwVtxDiff2","hTrkShwVtxDiff2",
                                1000,-500,500);
  hTrkShwVtxDiff2->SetTitle("hTrkShwVtxDiff2");
  hTrkShwVtxDiff2->GetXaxis()->SetTitle("hTrkShwVtxDiff2");
  hTrkShwVtxDiff2->GetXaxis()->CenterTitle();
  hTrkShwVtxDiff2->SetFillColor(0);
  hTrkShwVtxDiff2->SetLineWidth(3);
  //hTrkShwVtxDiff2->SetBit(TH1::kCanRebin);
  
  TH1F* hTrkShwEndDiff2=new TH1F("hTrkShwEndDiff2","hTrkShwEndDiff2",
                                1000,-500,500);
  hTrkShwEndDiff2->SetTitle("hTrkShwEndDiff2");
  hTrkShwEndDiff2->GetXaxis()->SetTitle("hTrkShwEndDiff2");
  hTrkShwEndDiff2->GetXaxis()->CenterTitle();
  hTrkShwEndDiff2->SetFillColor(0);
  hTrkShwEndDiff2->SetLineWidth(3);
  //hTrkShwEndDiff2->SetBit(TH1::kCanRebin);

  TH1F* hDiffYEnds=new TH1F("hDiffYEnds","hDiffYEnds",201,-4.5,4.5);
  hDiffYEnds->SetTitle("hDiffYEnds");
  hDiffYEnds->GetXaxis()->SetTitle("hDiffYEnds");
  hDiffYEnds->GetXaxis()->CenterTitle();
  hDiffYEnds->SetFillColor(0);
  hDiffYEnds->SetLineWidth(3);
  //hDiffYEnds->SetBit(TH1::kCanRebin);

  TH1F* hDiffYEndsSR=new TH1F("hDiffYEndsSR","hDiffYEndsSR",201,-4.5,4.5);
  hDiffYEndsSR->SetTitle("hDiffYEndsSR");
  hDiffYEndsSR->GetXaxis()->SetTitle("hDiffYEndsSR");
  hDiffYEndsSR->GetXaxis()->CenterTitle();
  hDiffYEndsSR->SetFillColor(0);
  hDiffYEndsSR->SetLineWidth(3);
  //hDiffYEndsSR->SetBit(TH1::kCanRebin);

  TH1F* hEvtShws=new TH1F("hEvtShws","hEvtShws",100,-50,50);
  hEvtShws->SetTitle("hEvtShws");
  hEvtShws->GetXaxis()->SetTitle("hEvtShws");
  hEvtShws->GetXaxis()->CenterTitle();
  hEvtShws->SetFillColor(0);
  hEvtShws->SetLineWidth(3);
  //hEvtShws->SetBit(TH1::kCanRebin);

  TProfile* pEnVsPl=new TProfile("pEnVsPl","pEnVsPl",1000,-500,500);
  pEnVsPl->SetTitle("Signal vs Distance from Trk Vtx");
  pEnVsPl->GetXaxis()->SetTitle("Distance (Planes)");
  pEnVsPl->GetXaxis()->CenterTitle();
  pEnVsPl->GetYaxis()->SetTitle("Signal");
  pEnVsPl->GetYaxis()->CenterTitle();
  pEnVsPl->SetLineColor(1);
  pEnVsPl->SetFillColor(0);
  //pEnVsPl->SetBit(TH1::kCanRebin);

  TProfile* pEnVsPlShw=new TProfile
    ("pEnVsPlShw","pEnVsPlShw",1000,-500,500);
  pEnVsPlShw->SetTitle("Signal vs Distance from Trk Vtx (w/ shw)");
  pEnVsPlShw->GetXaxis()->SetTitle("Distance (Planes)");
  pEnVsPlShw->GetXaxis()->CenterTitle();
  pEnVsPlShw->GetYaxis()->SetTitle("Signal");
  pEnVsPlShw->GetYaxis()->CenterTitle();
  pEnVsPlShw->SetLineColor(1);
  pEnVsPlShw->SetFillColor(0);
  //pEnVsPlShw->SetBit(TH1::kCanRebin);

  TProfile* pEnVsPlNoShw=new TProfile
    ("pEnVsPlNoShw","pEnVsPlNoShw",1000,-500,500);
  pEnVsPlNoShw->SetTitle("Signal vs Distance from Trk Vtx (no shw)");
  pEnVsPlNoShw->GetXaxis()->SetTitle("Distance (Planes)");
  pEnVsPlNoShw->GetXaxis()->CenterTitle();
  pEnVsPlNoShw->GetYaxis()->SetTitle("Signal");
  pEnVsPlNoShw->GetYaxis()->CenterTitle();
  pEnVsPlNoShw->SetLineColor(1);
  pEnVsPlNoShw->SetFillColor(0);
  //pEnVsPlNoShw->SetBit(TH1::kCanRebin);

  TProfile* pEnVsPlFromShw=new TProfile
    ("pEnVsPlFromShw","pEnVsPlFromShw",1000,-500,500);
  pEnVsPlFromShw->SetTitle("Signal vs Distance from Shw End");
  pEnVsPlFromShw->GetXaxis()->SetTitle("Distance (Planes)");
  pEnVsPlFromShw->GetXaxis()->CenterTitle();
  pEnVsPlFromShw->GetYaxis()->SetTitle("Signal");
  pEnVsPlFromShw->GetYaxis()->CenterTitle();
  pEnVsPlFromShw->SetLineColor(1);
  pEnVsPlFromShw->SetFillColor(0);
  //pEnVsPlFromShw->SetBit(TH1::kCanRebin);

  Int_t evtCounter=0;
  Int_t oneTrackCounter=0;
  Int_t oneEvtPerSlcCounter=0;
  Int_t longTrackCounter=0;
  Int_t closeTimesCounter=0;
  Int_t endPlaneSpectCounter=0;
  Int_t longTrack2Counter=0;

  Int_t badDirCounter=0;
  Int_t badDirSRCounter=0;
  
  Int_t FCCounter=0;
  Int_t PCCounter=0;
  Int_t TGCounter=0;
  Int_t badWindowCounter=0;
  Int_t badRecoCounter=0;
  Int_t goodWindowCounter=0;
  Int_t badFidCounter=0;
  Int_t badViewsCounter=0;
  Int_t badEndGapsCounter=0;
  Int_t badTrackEndCounter=0;
  Int_t bad2ndContCounter=0;
  Int_t badDistEndCounter=0;
  Int_t badShwDistCounter=0;
  
  string temps="myeventsLowEn";
  if (fRun>10000) temps+="MC";
  ofstream& eventsTxtFile=*(this->OpenTxtFile(100,"myevents"));
  ofstream& txtFileLowEn=*(this->OpenTxtFile(100,temps.c_str()));
  ofstream& txtFileNonPC=*(this->OpenTxtFile(100,"myeventsNonPC"));
  
  //reconstruction object
  MeuReco reco;
  MeuCuts cuts;
  MeuSummary s;
  
  //variable to turn on all the useful messages if required
  Msg::LogLevel_t logLevel=Msg::kDebug;

  
  this->InitialiseLoopVariables();  
  
  //for(Int_t entry=55000;entry<fEntries;entry++){
  for(Int_t entry=0;entry<fEntries;entry++){
      
    this->SetLoopVariables(entry);
    
    MSG("MeuAnalysis",logLevel)
      <<"Entry="<<entry<<endl;

    //make a reference instead of a pointer
    NtpStRecord& ntp=(*fRec);
    const RecCandHeader& rec=ntp.GetHeader();
    MAXMSG("MeuAnalysis",Msg::kInfo,1)
      <<"First entry: run="<<rec.GetRun()
      <<", snarl="<<rec.GetSnarl()<<endl;
    
    TClonesArray& evtTca=(*ntp.evt);
    const Int_t numEvts=evtTca.GetEntriesFast();

    //loop over the events
    for (Int_t ievt=0;ievt<numEvts;++ievt) {
      const NtpSREvent* pevt=
        dynamic_cast<NtpSREvent*>(evtTca[ievt]);
      const NtpSREvent& evt=(*pevt);

      evtCounter++;

      //ensure there is only 1 track in the event
      if (evt.ntrack!=1) continue;
      oneTrackCounter++;
        
      //filter out events occuring in a slice with more than 1 event
      if (cuts.FilterBadEvtPerSlc(ntp,evt.slc,entry)) continue;
      oneEvtPerSlcCounter++;
      
      TClonesArray& trkTca=(*ntp.trk);
      const NtpSRTrack* ptrk=
        dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[0]]);
      const NtpSRTrack& trk=(*ptrk);

      //MeuSummary s;
      s.Clear();
      s.RFid=1.0;
      s.DistToEdgeFid=0.4;
      s.Event=entry;
      cuts.FillSTSumDetails(ntp,s,ievt,evt.slc);
      //a map to hold the info for each plane in the entry
      map<Int_t,MeuHitInfo> plInfo;
      
      //cut on number of planes
      if (sqrt(pow(1.*trk.plane.beg-trk.plane.end,2))<8) continue;
      longTrackCounter++;
      
      MSG("MeuAnalysis",logLevel)
        <<"Entry="<<entry<<", run="<<rec.GetRun()
        <<", snarl="<<rec.GetSnarl()
        <<", trk.nstrips="<<trk.nstrip<<endl;
      
      if (s.Snarl==-1){
        MSG("MeuAnalysis",Msg::kInfo)
          <<"Entry="<<entry<<", run="<<rec.GetRun()
          <<", snarl="<<rec.GetSnarl()
          <<", trk.nstrips="<<trk.nstrip<<endl;
        cuts.PrintNtpSt(ntp);
      }

      if (cuts.FilterBadTrkTimes(ntp,s,evt)) continue;
      closeTimesCounter++;
      
      cuts.ExtractPlInfo(ntp,s,plInfo,evt);
      
      reco.CalcVtx(s,plInfo);
      
      hDiffYEndsSR->Fill(plInfo[s.VtxPlane].Y-
                         plInfo[s.EndPlane].Y);

      cuts.ExtractTruthInfo(ntp,s,plInfo);
      Bool_t mcPosDir=s.MCVtxPlane<s.MCEndPlane;
      Bool_t srPosDir=trk.plane.beg<trk.plane.end;
      Bool_t posDir=s.VtxPlane<s.EndPlane;
      
      if (posDir!=mcPosDir && s.SimFlag==SimFlag::kMC){
        badDirCounter++;
        MAXMSG("MeuAnalysis",Msg::kVerbose,200)
          <<"Direction wrong! Trk: truth "
          <<s.MCVtxPlane<<" -> "<<s.MCEndPlane
          <<", reco "<<s.VtxPlane<<" -> "<<s.EndPlane
          //<<", y="<<s.VtxY<<" -> "<<s.EndY
          <<endl;
      }
      if (srPosDir!=mcPosDir && s.SimFlag==SimFlag::kMC){
        badDirSRCounter++;
        MAXMSG("MeuAnalysis",Msg::kVerbose,100)
          <<"Direction wrong in SR! Trk: truth "
          <<s.MCVtxPlane<<" -> "<<s.MCEndPlane
          <<", SR reco "<<trk.plane.beg<<" -> "<<trk.plane.end
          <<endl;
      }

      //check if track stops in spectrometer
      if (s.EndPlane>120) continue;
      endPlaneSpectCounter++;
      
      MSG("MeuAnalysis",logLevel)<<"CalcFCPC..."<<endl;
      reco.CalcFCPC(s,plInfo);
      
      //recalculate the vtx of the tracks
      //coming in from the spectrometer
      MSG("MeuAnalysis",logLevel)<<"CalcVtxSpecialND..."<<endl;
      reco.CalcVtxSpecialND(s,plInfo);

      //cut on number of planes with potential new vertex
      if (sqrt(pow(1.*s.VtxPlane-s.EndPlane,2))<8) continue;
      longTrack2Counter++;
        
      Bool_t TG=!s.FC && !s.PC;
      
      if (s.PC) PCCounter++;
      else if (s.FC) FCCounter++;
      else if (TG) TGCounter++;
      if (s.FC && s.PC) MSG("MeuAnalysis",Msg::kError)
        <<"Both FC and PC!"<<endl;
      
      if (s.PC) {
        
        Int_t diffAtEnd=-1;
        Int_t diffAtVtx=-1;
        if (!cuts.CheckTrkViews(s,&diffAtVtx,&diffAtEnd)) {
          badViewsCounter++;
          continue;
        }
        
        MSG("MeuAnalysis",logLevel)<<"Found PC event..."<<endl;
        Bool_t goodReco=reco.Reconstruct(s,plInfo);
        if (goodReco){
          MSG("MeuAnalysis",logLevel)<<"Recalibrating..."<<endl;
          this->Recalibrate(ntp,s,plInfo,evt);
          
          MSG("MeuAnalysis",logLevel)<<"CalcWindow..."<<endl;
          reco.CalcWindow(s,plInfo);
          
          //recalculate the containment!
          s.DistToEdgeFid=0.35;//allow a slight shift in position
          reco.CalcFCPC(s,plInfo);
          
          reco.CalcStripDists(s,plInfo);
          
          MSG("MeuAnalysis",logLevel)<<"CheckContainment..."<<endl;
          //fill the tree once event has been analysed
          if (s.WinSigMap<=0) {
            badWindowCounter++;
          }
          else if (!reco.CheckWinContainment(s,plInfo)){
            badFidCounter++;
            //MSG("MeuAnalysis",Msg::kInfo)<<"Outside fid:"<<endl;
          }
          else if (cuts.FilterBadTrackEnd(s,plInfo)){
            badTrackEndCounter++;
          }
          else if (cuts.FilterBadEndGaps(s,plInfo)){
            badEndGapsCounter++;
          }
          else if (!s.PC){
            bad2ndContCounter++;
            txtFileNonPC<<entry<<endl;

            Bool_t print=false;
            if (print){
              MAXMSG("MeuAnalysis",Msg::kDebug,50)
                <<endl<<endl<<"Recalculated containment and not PC!"
                <<" PC="<<s.PC<<", FC="<<s.FC
                <<" entry="<<entry
                <<", snarl="<<s.Snarl<<", run="<<s.Run<<endl
                <<"Low En="<<s.MCLowEn<<" high="<<s.MCHighEn<<endl
                <<endl;
              MAXMSG("MeuAnalysis",Msg::kDebug,50)
                <<"vtx Pl="<<s.VtxPlane<<", end Pl="<<s.EndPlane
                <<", dVtx="<<s.VtxDistToEdge
                <<", dEnd="<<s.EndDistToEdge<<endl;
              
              for (map<Int_t,MeuHitInfo>::const_iterator it=
                     plInfo.begin();
                   it!=plInfo.end();++it){
                const MeuHitInfo& c=it->second;
                MAXMSG("MeuAnalysis",Msg::kInfo,500)
                  <<"pl="<<c.Plane<<", z="<<c.Z
                  <<", x="<<c.X<<", y="<<c.Y
                  <<"  t="<<c.TPos<<", l="<<c.LPos<<endl;
              }          
            }
          }
          else if (cuts.FilterBadDistEndStrip(s,plInfo)){
            badDistEndCounter++;
          }
          else {
            hTrkVtxPl->Fill(s.VtxPlane);
            hTrkEndPl->Fill(s.EndPlane);
            hTrkLengthPl->Fill(s.EndPlane-s.VtxPlane);
            hDiffYEnds->Fill(s.VtxY-s.EndY);
            hEvtShws->Fill(evt.nshower);

            Int_t shwVtxPl=999;
            Int_t shwEndPl=-999;
            
            //loop over showers in evt
            TClonesArray& shwTca=(*ntp.shw);
            for(int i=0;i<evt.nshower;i++) {    
              const NtpSRShower* pshw=
                dynamic_cast<NtpSRShower*>(shwTca[evt.shw[i]]);
              const NtpSRShower& shw=(*pshw);
              
              hShwVtxPl->Fill(shw.plane.beg);
              hShwEndPl->Fill(shw.plane.end);
              hShwLengthPl->Fill(shw.plane.end-shw.plane.beg);

              hTrkShwVtxDiff->Fill(s.VtxPlane-shw.plane.beg);
              hTrkShwEndDiff->Fill(s.EndPlane-shw.plane.end);

              //define a vtx shower to be one that starts within 5 pl
              Bool_t vtxShw=abs(s.VtxPlane-shw.plane.beg)<=5;
              
              if (vtxShw){
                if (shw.plane.beg<shwVtxPl) shwVtxPl=shw.plane.beg;
                if (shw.plane.end>shwEndPl) shwEndPl=shw.plane.end;
              }
            }

            Bool_t vtxShw=false;
            if (shwEndPl!=-999 && shwVtxPl!=999) {
              vtxShw=true;
              hShwVtxPl2->Fill(shwVtxPl);
              hShwEndPl2->Fill(shwEndPl);
              hShwLengthPl2->Fill(shwEndPl-shwVtxPl);
              hTrkShwVtxDiff2->Fill(s.VtxPlane-shwVtxPl);
              hTrkShwEndDiff2->Fill(s.EndPlane-shwEndPl);
            }

            //fill profiles
            typedef map<Int_t,MeuHitInfo>::const_iterator plIt;
            for (plIt it=plInfo.begin();it!=plInfo.end();++it){
              const MeuHitInfo& cp=it->second;

              //check that plane is within track (assume forward going)
              if (cp.Plane<s.VtxPlane ||
                  cp.Plane>s.EndPlane) continue;

              Int_t plFromVtx=cp.Plane-s.VtxPlane;
              pEnVsPl->Fill(plFromVtx,cp.SigMap);
              
              if (vtxShw) {
                pEnVsPlShw->Fill(plFromVtx,cp.SigMap);
                Int_t plFromShwEnd=cp.Plane-shwEndPl;
                if (plFromShwEnd>=0){
                  pEnVsPlFromShw->Fill(plFromShwEnd,cp.SigMap);
                }
              }
              else pEnVsPlNoShw->Fill(plFromVtx,cp.SigMap);
            }

            if (vtxShw){
              Int_t winPlFromShwEnd=s.WinVtxSidePl-shwEndPl;
              hWinPlFromShwEnd->Fill(winPlFromShwEnd);
              
              //require that shw stops 10 planes before window starts
              if (winPlFromShwEnd<10 && s.TrigSrc>100) {
                badShwDistCounter++;
                continue;
              }
            }
            

            goodWindowCounter++;
            
            cuts.ExtractTruthInfo(ntp,s,plInfo);
            
            cuts.FillSTSumRecoDetails(s,plInfo);
            
            hWinSigMap->Fill(s.WinSigMap);
            hWinSigCor->Fill(s.WinSigCor);
            hWinSigLin->Fill(s.WinSigLin);
            hWinAdc->Fill(s.WinAdc);
            hWinPe->Fill(s.WinPe);
            
            eventsTxtFile<<entry<<endl;
            MAXMSG("MeuAnalysis",Msg::kInfo,10)
              <<"Good window: run="
              <<rec.GetRun()<<", snarl="<<rec.GetSnarl()
              <<endl;
            
            if (s.MCLowEn>0.5*Munits::GeV &&
                s.SimFlag==SimFlag::kMC) {
              txtFileLowEn<<entry<<endl;
              MSG("MeuAnalysis",Msg::kInfo)
                <<"Low En="<<s.MCLowEn<<" high="<<s.MCHighEn
                <<", vtxPl="<<s.VtxPlane<<", endPl="<<s.EndPlane
                <<endl<<"entry="<<entry
                <<", snarl="<<s.Snarl<<", run="<<s.Run<<endl;
            }

            cuts.FillEnVsDist(s,plInfo);
            cuts.FillEnVsDist2(s,plInfo);
            cuts.FillTimeHistos(ntp,evt,s);
          }
        }
        else {
          badRecoCounter++;
          MSG("MeuAnalysis",logLevel)<<"Bad reco..."<<endl;
        }
      }
    }//end of loop over events
  }//end of for                                       
    

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

  Double_t quantile=0.5;//quantile to compute

  Double_t meuWinSigMap=-1;
  hWinSigMap->GetQuantiles(1,&meuWinSigMap,&quantile);
  Double_t errWinSigMap=hWinSigMap->GetRMS()/
    sqrt(hWinSigMap->GetEntries());
  Double_t relErrWinSigMap=errWinSigMap/meuWinSigMap;
  
  Double_t meuWinSigCor=-1;
  hWinSigCor->GetQuantiles(1,&meuWinSigCor,&quantile);
  Double_t errWinSigCor=hWinSigCor->GetRMS()/
    sqrt(hWinSigCor->GetEntries());
  Double_t relErrWinSigCor=errWinSigCor/meuWinSigCor;

  Double_t meuWinSigLin=-1;
  hWinSigLin->GetQuantiles(1,&meuWinSigLin,&quantile);
  Double_t errWinSigLin=hWinSigLin->GetRMS()/
    sqrt(hWinSigLin->GetEntries());
  Double_t relErrWinSigLin=errWinSigLin/meuWinSigLin;

  Double_t meuWinAdc=-1;
  hWinAdc->GetQuantiles(1,&meuWinAdc,&quantile);
  Double_t errWinAdc=hWinAdc->GetRMS()/
    sqrt(hWinAdc->GetEntries());
  Double_t relErrWinAdc=errWinAdc/meuWinAdc;

  Double_t meuWinPe=-1;
  hWinPe->GetQuantiles(1,&meuWinPe,&quantile);
  Double_t errWinPe=hWinPe->GetRMS()/
    sqrt(hWinPe->GetEntries());
  Double_t relErrWinPe=errWinPe/meuWinPe;

  hWinPe->Fit("gaus");
  TF1* fPe=hWinPe->GetFunction("gaus");
  Double_t errWinPe2=fPe->GetParameter(2)/sqrt(hWinPe->GetEntries());
  
//  cout<<"RESOLUTIONS (from Gaussian fit):"<<endl
  //    <<"ADCs:   "<<fAdc->GetParameter(2)/fAdc->GetParameter(1)<<endl
  //  <<"PEs:    "<<fPe->GetParameter(2)/fPe->GetParameter(1)<<endl
  
  MSG("MeuAnalysis",Msg::kInfo) 
    <<endl<<"*** Run Summary ***"<<endl
    //<<"Raw MEU="<<meu<<" +/- "<<err<<" ("<<100.*relErr<<"%)"<<endl
    <<"Raw MEU (WinSigMap)="<<meuWinSigMap<<" +/- "<<errWinSigMap
    <<" ("<<100.*relErrWinSigMap<<"%)"<<endl
    <<"Raw MEU (WinSigCor)="<<meuWinSigCor<<" +/- "<<errWinSigCor
    <<" ("<<100.*relErrWinSigCor<<"%)"<<endl
    <<"Raw MEU (WinSigLin)="<<meuWinSigLin<<" +/- "<<errWinSigLin
    <<" ("<<100.*relErrWinSigLin<<"%)"<<endl
    <<"Raw MEU (WinAdc)   ="<<meuWinAdc<<" +/- "<<errWinAdc
    <<" ("<<100.*relErrWinAdc<<"%)"<<endl
    <<"Raw MEU (WinPe)    ="<<meuWinPe<<" +/- "<<errWinPe
    <<" ("<<100.*relErrWinPe<<"%), fitErr="<<errWinPe2<<endl;
    //<<endl
    /*<<"Not LI="<<notLICounter<<endl
      <<"Not MM="<<notMMCounter<<endl
      <<"Good track="<<goodTrackCounter<<endl
      <<"Good length="<<goodLengthCounter<<endl
      <<"Good end & vtx="<<goodEndAndVtxCounter<<endl
      <<"Good XY reco="<<goodXYCounter<<endl
      <<endl
      <<"Good containmentPC="<<containmentCounter<<endl
      <<"  both SM hit (when PC)="<<bothSMCounter
      <<"   ("<<bothSM<<"% of PC)"<<endl
      <<"  coil hit (when PC)   ="<<coilHitCounter
      <<"   ("<<coil_PC<<"% of PC)"<<endl
      <<"  close LI (when PC)   ="<<closeLICounter
      <<"   ("<<closeLI<<"% of PC)"<<endl
      <<"  dead Chips (when PC) ="<<deadChipsCounter
      <<"   ("<<deadChips<<"% of PC)"<<endl
      <<"  boundary TF (when PC)="<<boundaryTFCounter
      <<"   ("<<boundaryTF<<"% of PC)"<<endl
      <<"  bad views (when PC)  ="<<badViewsCounter
      <<"   ("<<badViews<<"% of PC)"<<endl
      <<"  bad trk end (when PC)="<<badTrkEndCounter
      <<"   ("<<badTrkEnd<<"% of PC)"<<endl
      <<"  bad trk gap (when PC)="<<badTrkGapsCounter
      <<"   ("<<badTrkGaps<<"% of PC)"<<endl
      <<"  bad traceZ (when PC) ="<<badTraceZCounter
      <<"   ("<<badTraceZ<<"% of PC)"<<endl
      <<"  bad window (when PC) ="<<badWindowCounter
      <<"   ("<<badWin<<"% of PC)"<<endl
      <<"Good containmentFC="<<containmentCounterFC<<endl
      <<"  both SM hit (when FC)="<<bothSMCounterFC
      <<"   ("<<bothSMFC<<"% of FC)"<<endl
      <<"  coil hit (when FC)   ="<<coilHitCounterFC
      <<"   ("<<coil_FC<<"% of FC)"<<endl
      <<"Good containmentTG="<<containmentCounterTG<<endl
      <<"  both SM hit (when TG)="<<bothSMCounterTG
      <<"   ("<<bothSMTG<<"% of TG)"<<endl
      <<"  coil hit (when TG)   ="<<coilHitCounterTG
      <<"   ("<<coil_TG<<"% of TG)"<<endl
      <<endl
      <<"Bad charge pos moment="<<momentCounter<<endl
      <<"Too high StripsPerPlane="<<stPerPlCounter<<endl
      <<"  Either of above two="<<stPlOrMomentCounter<<endl
      <<endl
      <<"MC Info:"<<endl
      <<"  Muon lowest energy >0.6 GeV="<<lowEnCounter<<endl
      <<"  Muon lowest energy >5.0 GeV="<<lowEn5Counter<<endl
      <<endl
    */

  Int_t totalBadPCEvents=badWindowCounter+badFidCounter+
    badRecoCounter+badViewsCounter+
    badTrackEndCounter+badEndGapsCounter+bad2ndContCounter+
    badDistEndCounter+badShwDistCounter;
  
  MSG("MeuAnalysis",Msg::kInfo)
    <<"Total Entries (snarls)="<<fEntries<<endl
    <<"Total Evts="<<evtCounter<<endl
    <<"  One track      ="<<oneTrackCounter<<endl
    <<"  One evt in slc ="<<oneEvtPerSlcCounter<<endl
    <<"  Long track     ="<<longTrackCounter<<endl
    <<"  Close trk times="<<closeTimesCounter<<endl
    <<"  End plane cal  ="<<endPlaneSpectCounter<<endl
    <<"  Long track 2   ="<<longTrack2Counter<<endl
    <<"  Bad dir (MC)   ="<<badDirCounter<<endl
    <<"  Bad dir SR (MC)="<<badDirSRCounter<<endl
    <<"Containment summary:"<<endl
    <<"  PC="<<PCCounter<<endl
//    <<"     (PC - all cuts="
    //closeLICounter-deadChipsCounter-boundaryTFCounter-
    //badViewsCounter-badTrkEndCounter-badTrkGapsCounter-badTraceZCounter-
    //badWindowCounter<<")"<<endl
    <<"  FC="<<FCCounter<<endl//<<"     (FC-coil-bothSM="
    //<<FCCounter-coilHitCounterFC-bothSMCounterFC<<")"<<endl
    <<"  TG="<<TGCounter<<endl//<<"     (TG-coil-bothSM="
//    <<TGCounter-coilHitCounterTG-bothSMCounterTG<<")"<<endl
    <<"  FC+PC+TG="<<FCCounter+PCCounter+TGCounter<<endl
    <<endl
    <<"Total PC with good track window = "<<goodWindowCounter<<endl
    <<"              bad track window  = "<<badWindowCounter<<endl
    <<"              bad fiducial vol  = "<<badFidCounter<<endl
    <<"              bad reco          = "<<badRecoCounter<<endl
    <<"              bad views         = "<<badViewsCounter<<endl
    <<"              bad track end     = "<<badTrackEndCounter<<endl
    <<"              bad track end gaps= "<<badEndGapsCounter<<endl
    <<"              bad 2nd contain.  = "<<bad2ndContCounter<<endl
    <<"              bad dist strip end= "<<badDistEndCounter<<endl
    <<"              bad dist to shw   = "<<badShwDistCounter<<endl
    <<" Total bad                      = "<<totalBadPCEvents<<endl
    <<" Total bad + good               = "
    <<totalBadPCEvents+goodWindowCounter<<endl
    <<endl;

  if (fRec) delete fRec;

  MSG("MeuAnalysis",Msg::kInfo) 
    <<" ** Finished SpillPlots method **"<<endl;
}

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

Bool_t MeuAnalysis::FilterLI(const NtpStRecord& ntp,
                             const MeuSummary& s) const
{
  //only analyse if FD
  if (s.Detector!=Detector::kFar) return false;

  const Bool_t cleanTheEvent=false;

  Float_t sumSigCorEast=0;
  Float_t sumSigCorWest=0;

  map<Int_t,Int_t> hpp;  // <plane, # of hits>
  multiset<Double_t> times;
  
  Int_t numStrips=ntp.stp->GetEntries();
  if (numStrips<=0) return false;//pass the event, not LI
  //loop over the strips in the snarl
  for (Int_t hit=0;hit<numStrips;hit++){
    NtpSRStrip* strip = dynamic_cast<NtpSRStrip*>((*ntp.stp)[hit]);

    Double_t time=strip->time0;
    if (time<0 || time>1) time=strip->time1;

    if (time>0 && time<=1) {
      times.insert(time);
    }
    //else just don't put the time in the map - clearly rubbish
  }//end of loop over strips

  //get the median time from the map
  multiset<Double_t>::const_iterator it=times.begin();
  advance(it,times.size()/2);
  Double_t medianTime=*it;
  MAXMSG("MeuAnalysis",Msg::kDebug,100)
    <<"FilterLI: Median time="<<medianTime<<endl;

  //second loop
  for (Int_t hit=0;hit<numStrips;hit++){
    NtpSRStrip* strip = dynamic_cast<NtpSRStrip*>((*ntp.stp)[hit]);

    if (cleanTheEvent) {
      if (strip->ph0.sigcor > 0){
        if (strip->time0<medianTime-200e-9 || 
            strip->time0>medianTime+200e-9) continue;
      }
      if (strip->ph1.sigcor > 0){
        if (strip->time1<medianTime-200e-9 || 
            strip->time1>medianTime+200e-9) continue;
      }
   } // if (cleanTheEvent)

    sumSigCorEast+=strip->ph0.sigcor;
    sumSigCorWest+=strip->ph1.sigcor;

    hpp[strip->plane] += 1;
  }//end of loop over strips

  Float_t avHpp=0;
  vector<Float_t> pb(8);//to store fractions hit in each pb region

  //loop over all the planes flashed
  for (map<Int_t,Int_t>::iterator hppIt=hpp.begin();
       hppIt!=hpp.end();++hppIt){
    avHpp+=hppIt->second;

    Int_t pl=hppIt->first;

    if (pl>=1 && pl<=64) pb[0]+=1./64;//0=bookend
    else if (pl>=65 && pl<=128) pb[1]+=1./64;
    else if (pl>=129 && pl<=192) pb[2]+=1./64;
    else if (pl>=193 && pl<=248) pb[3]+=1./56;//249=bookend
    else if (pl>=250 && pl<=313) pb[4]+=1./64;
    else if (pl>=314 && pl<=377) pb[5]+=1./64;
    else if (pl>=378 && pl<=441) pb[6]+=1./64;
    else if (pl>=442 && pl<=485) pb[7]+=1./44;

    MAXMSG("MeuAnalysis",Msg::kDebug,2000)
      <<"plane="<<pl<<", hpp="<<hppIt->second<<endl;
  }
  if (hpp.size()>0) avHpp/=hpp.size();
  Float_t sumSig=sumSigCorEast+sumSigCorWest;
  Float_t minSigCorEW = TMath::Min(sumSigCorEast, sumSigCorWest);
  Float_t asym=0;
  
  if (sumSig) asym = (sumSig - 2 * minSigCorEW) / sumSig;

  Float_t fractFlashed1=0;
  Float_t fractFlashed2=0;

  for (vector<Float_t>::iterator it=pb.begin();
       it!=pb.end();++it){
    Float_t fract=*it;
    if (fract>fractFlashed1){
      //copy the 1 into 2
      fractFlashed2=fractFlashed1;
      //then set new highest fractFlashed
      fractFlashed1=fract;
    }
    //make sure to get fractFlashed2 if fractFlashed2<fract<frFla1
    else if (fract>fractFlashed2) {
      fractFlashed2=fract;
    }

    MAXMSG("MeuAnalysis",Msg::kDebug,200)
      <<"fract="<<fract<<", f1="<<fractFlashed1
      <<", f2="<<fractFlashed2<<endl;
  }

  Float_t ratioFlashed=0;
  if (fractFlashed1>0) ratioFlashed=fractFlashed2/fractFlashed1;

  Bool_t li=false;

  // check if LI
  // don't check the asymmetry if the VA is saturating
  if ( (asym>0.5 || minSigCorEW>1.7e6) && 
       avHpp>3 && ratioFlashed<0.05 &&
       fractFlashed1>0.85 && fractFlashed2<0.1 ){
    li=true;
    
    MAXMSG("MeuAnalysis",Msg::kInfo,10)
      <<"LI Event:"<<endl
      <<"  Asymmetry="<<asym<<endl
      <<"  Av. HPP  ="<<avHpp<<endl
      <<"  Fract Fl1="<<fractFlashed1
      <<", Fract Fl2="<<fractFlashed2
      <<"  Ratio Fla="<<ratioFlashed<<endl;
  }
  else{
    MAXMSG("MeuAnalysis",Msg::kDebug,100)
      <<"Non-LI Event:"<<endl
      <<"  Asymmetry="<<asym<<endl
      <<"  Av. HPP  ="<<avHpp<<endl
      <<"  Fract Fl1="<<fractFlashed1
      <<", Fract Fl2="<<fractFlashed2
      <<"  Ratio Fla="<<ratioFlashed<<endl;
  }
  
  return li;
}

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

Bool_t MeuAnalysis::FilterLI(const AtmosEvent& ev) const
{
  Float_t sumSigCorEast=0;
  Float_t sumSigCorWest=0;
  
  map<Int_t,Int_t> hpp;
  
  TClonesArray& strips=(*ev.StripList);
  Int_t numStrips=strips.GetEntries();
  if (numStrips<=0) return false;//pass the event, not LI
  //loop over the strips in the snarl
  for (Int_t hit=0;hit<numStrips;hit++){
    const AtmosStrip* strip=dynamic_cast<AtmosStrip*>(strips[hit]);
    const AtmosStrip& st=(*strip);
    
    sumSigCorEast+=st.QPEcorr[0];
    sumSigCorWest+=st.QPEcorr[1];
    
    hpp[st.Plane]++;
  }//end of loop over strips
  
  Float_t avHpp=0;
  vector<Float_t> pb(8);//to store fractions hit in each pb region
  
  //loop over all the planes flashed
  for (map<Int_t,Int_t>::iterator hppIt=hpp.begin();
       hppIt!=hpp.end();++hppIt){
    avHpp+=hppIt->second;
    
    Int_t pl=hppIt->first;
    
    //work out the fraction of the planes flashed by each pb
    //are actually hit
    if (pl>=1 && pl<=64) pb[0]+=1./64;//0=bookend
    else if (pl>=65 && pl<=128) pb[1]+=1./64;
    else if (pl>=129 && pl<=192) pb[2]+=1./64;
    else if (pl>=193 && pl<=248) pb[3]+=1./56;//249=bookend
    else if (pl>=250 && pl<=313) pb[4]+=1./64;
    else if (pl>=314 && pl<=377) pb[5]+=1./64;
    else if (pl>=378 && pl<=441) pb[6]+=1./64;
    else if (pl>=442 && pl<=485) pb[7]+=1./44;
  }
  if (hpp.size()>0) avHpp/=hpp.size();
  Float_t sumSig=sumSigCorEast+sumSigCorWest;
  Float_t minSigCorEW = TMath::Min(sumSigCorEast, sumSigCorWest);
  Float_t asym=0;
  if (sumSig) asym = (sumSig - 2 * minSigCorEW) / sumSig;

  Float_t fractFlashed1=0;
  Float_t fractFlashed2=0;
  
  for (vector<Float_t>::iterator it=pb.begin();
       it!=pb.end();++it){
    Float_t fract=*it;
    if (fract>fractFlashed1){
      //copy the 1 into 2
      fractFlashed2=fractFlashed1;
      //then set new highest fractFlashed
      fractFlashed1=fract;
    }
    MAXMSG("MeuAnalysis",Msg::kVerbose,200) 
      <<"fract="<<fract<<", f1="<<fractFlashed1
      <<", f2="<<fractFlashed2<<endl;
  }
  
  Float_t ratioFlashed=0;
  if (fractFlashed1>0) ratioFlashed=fractFlashed2/fractFlashed1;
  
  Bool_t li=false;
  
  //check if LI
  if ( (asym>0.5 || minSigCorEW>1.7e6) && 
       avHpp>3 && ratioFlashed<0.05 &&
       fractFlashed1>0.85 && fractFlashed2<0.1 ){
    li=true;
    
    MAXMSG("MeuAnalysis",Msg::kDebug,200) 
      <<"LI Event:"<<endl
      <<"  Asymmetry="<<asym<<endl
      <<"  Av. HPP  ="<<avHpp<<endl
      <<"  Fract Fl1="<<fractFlashed1<<", Fract Fl2="<<fractFlashed2
      <<"  Ratio Fla="<<ratioFlashed<<endl;
  }
  
  return li;
}

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

Float_t MeuAnalysis::GetTemperature(const MeuSummary& s) const
{
  static Bool_t firstCall=true;
  static Bool_t foundFile=false;
  static vector <Float_t> vTemperatures;
  static Int_t firstTime=-1;
  static Int_t lastTime=-1;

  //only get the temperature for the FD at present
  if (s.Detector!=Detector::kFar) {
    MAXMSG("MeuAnalysis",Msg::kInfo,1)
      <<"Not retrieving the temperature for dt="<<s.Detector<<endl;
    return 0;
  }
  
  if (firstCall){
    char* envVariable=getenv("MEUDATA");
    if (envVariable==NULL){
      MAXMSG("MeuAnalysis",Msg::kWarning,1)
        <<"Can't find ENV variable \"MEUDATA\" to"
        <<" open temperatures file"<<endl;
      firstCall=false;
      return 0;
    }
    string sEnv=envVariable;

    string sTemperatureFileName=sEnv+"/temps.dat";
    MSG("MeuAnalysis",Msg::kInfo)
      <<"GetTemperature: Looking for file: "
      <<sTemperatureFileName<<endl;
    ifstream temperatureFile(sTemperatureFileName.c_str());
    if(!temperatureFile){
      MAXMSG("MeuAnalysis",Msg::kWarning,1)
        <<"Can't find file of temperatures"<<endl;
      firstCall=false;
      foundFile=false;
      return 0;
    }

    foundFile=true;

    map<Int_t,Float_t> temperatures;

    MSG("MeuAnalysis",Msg::kInfo)
      <<"Loading temperatures from txt file..."<<endl;
    if (temperatureFile){
      Float_t Tf=-1;
      Float_t Tc=-1;
      Int_t time=0;
      while(temperatureFile>>time>>Tf) {
        Tc=(Tf-32)/1.8;
        temperatures[time]=Tc;
      }
    }
    
    firstTime=temperatures.begin()->first;
    lastTime=(--temperatures.end())->first;
    MSG("MeuAnalysis",Msg::kInfo)
      <<"Found first time="<<firstTime<<", lastTime="<<lastTime<<endl;
    VldTimeStamp first(firstTime,0);
    first.Print();
    VldTimeStamp last(lastTime,0);
    last.Print();
    
    Int_t lastIndex=(lastTime-firstTime)/300;
    MSG("MeuAnalysis",Msg::kInfo)<<"Size of table="<<lastIndex<<endl;
    vTemperatures.reserve(lastIndex+1);
    for (Int_t i=0;i<lastIndex+1;i++) vTemperatures.push_back(0);
    
    //loop and fill vector
    for (map<Int_t,Float_t>::iterator it=temperatures.begin();
         it!=temperatures.end();++it){
      Int_t time=it->first;
      Int_t index=(time-firstTime)/300;//every 5 mins
      MAXMSG("MeuAnalysis",Msg::kDebug,50)
        <<"time="<<it->first
        <<", mins since first="<<(time-firstTime)/300
        <<", temperature="<<it->second<<endl;
      vTemperatures[index]=it->second;
    } 
  }
  firstCall=false;
  
  if (!foundFile) return 0;

  Int_t ind=(s.TimeSec-firstTime)/300;
  Float_t Tc=0;
  Int_t lastIndex=(lastTime-firstTime)/300;
  Int_t timeSkipped=0;

  if (s.TimeSec>=firstTime && s.TimeSec<=lastTime){
    while (ind<=lastIndex && ind>0){
      MAXMSG("MeuAnalysis",Msg::kVerbose,500)
        <<"500 mins T="<<vTemperatures[ind]<<", ind="<<ind<<endl;
      if (vTemperatures[ind]>0){
        MAXMSG("MeuAnalysis",Msg::kDebug,50)
          <<"Found good T="<<vTemperatures[ind]<<endl;
        Tc=vTemperatures[ind];
        break;
      }
      timeSkipped++;
      ind++;
    }
  }
  else MSG("MeuAnalysis",Msg::kWarning)<<"Time is out of range"<<endl;

  //check that not too big a time (1 hour) was skipped
  if (timeSkipped<=12) return Tc;
  else {
    MAXMSG("MeuAnalysis",Msg::kInfo,50)
      <<"Too big a gap (>1 hour) in temperature reading"
      <<" to be reliable ("<<timeSkipped<<"*5 mins)."
      <<" Returning 0 degC"<<endl;
    return 0;
  }
}

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

void MeuAnalysis::MakeSummaryTree()
{
  //this method just adds a level of indirection so that you only 
  //need a single macro to make the summary tree

  if (fUseAtNu){
    this->MakeSummaryTreeWithAtNu();
  }
  else{
    this->MakeSummaryTreeWithNtpSt();
  }
}

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

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

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

  MeuSummaryWriter summaryOut;
  summaryOut.SummaryTreeSetup(fRun,fSubrun);
 
  //string temps="myeventsLowEn";
  //if (fRun>10000) temps+="MC";
  //ofstream& eventsTxtFile=*(this->OpenTxtFile(100,"myevents"));
  //ofstream& txtFileLowEn=*(this->OpenTxtFile(100,temps.c_str()));
  //ofstream& txtFileNonPC=*(this->OpenTxtFile(100,"myeventsNonPC"));
 
  //reconstruction object
  const MeuReco reco;
  const MeuCuts cuts;
  const MeuHistos histos;
  MeuCounter cnt;

  //variable to turn on all the useful messages if required
  Msg::LogLevel_t logLevel=Msg::kDebug;

  
  this->InitialiseLoopVariables();  
  
  //for(Int_t entry=24196;entry<fEntries;entry++){
  //for(Int_t entry=20307;entry<fEntries;entry++){//crashed on 32998
  //for(Int_t entry=30000;entry<fEntries;entry++){//didn't crash 32998
  //for(Int_t entry=0;entry<2000;entry++){
  for(Int_t entry=0;entry<fEntries;entry++){
    
    this->SetLoopVariables(entry);
    
    MSG("MeuAnalysis",logLevel)
      <<"Entry="<<entry<<endl;
    
    //make a reference instead of a pointer
    AtmosEvent& ev=(*fev);
    MAXMSG("MeuAnalysis",Msg::kInfo,1)
      <<"First entry: run="<<ev.Run
      <<", snarl="<<ev.Snarl<<endl;

    //MSG("MeuAnalysis",Msg::kInfo)
    //<<"run="<<ev.Run
    //<<", snarl="<<ev.Snarl
    //<<", entry="<<entry
    //<<", unixTime="<<ev.UnixTime
    //<<", #scintHits="<<ev.NScintHits
    //<<endl;

    //cut out the LI events
    if (this->FilterLI(ev)){
      //eventsTxtFileLI<<event<<endl;
      continue;
    }
    cnt.notLICounter++;

    //cut out events with multiple tracks
    if (ev.NTracks!=1) continue;
    cnt.oneTrackCounter++;

    const AtmosTrack* track=dynamic_cast<const AtmosTrack*>
      (ev.TrackList->At(0));
    const AtmosTrack& trk=(*track);
    
    //cut out events with less than 20 planes
    Bool_t goodLength=trk.AtNuNplanes>20;
    if (!goodLength){
      continue;
    }
    cnt.longTrackCounter++;

    //MeuSummary s;
    MeuSummary& s=summaryOut.GetMeuSummaryToFill();
    s.RFid=3.0*Munits::m;
    s.RFidCoil=0.5*Munits::m;
    s.DistToEdgeFid=1.0*Munits::m;
    s.Event=entry;
    cuts.FillSTSumDetails(ev,s);
    //a map to hold the info for each plane in the entry
    map<Int_t,MeuHitInfo> plInfo;

    cuts.ExtractPlInfo(ev,s,plInfo);
    
    reco.CalcVtx(s,plInfo);
      
    MSG("MeuAnalysis",logLevel)<<"CalcFCPC..."<<endl;
    reco.CalcFCPC(s,plInfo);

    //check that the XY position is sensible
    if (cuts.FilterBadXY(ev,s)){
      continue;
    }
    cnt.goodXYCounter++;
    
    Bool_t TG=!s.FC && !s.PC;    
    if (s.PC) cnt.PCCounter++;
    else if (s.FC) cnt.FCCounter++;
    else if (TG) cnt.TGCounter++;
    if (s.FC && s.PC) MSG("MeuAnalysis",Msg::kError)
      <<"Both FC and PC!"<<endl;

    Bool_t awayFromCoil=cuts.AnalyseCoilProximity(ev,s);
    if (s.PC && !awayFromCoil) cnt.coilHitCounter++;
    if (s.PC && s.SMBoth && awayFromCoil) cnt.bothSMHitCounter++;

    Bool_t goodContainment=s.PC && awayFromCoil && !s.SMBoth;
    if (!goodContainment) continue;
    MSG("MeuAnalysis",logLevel)<<"Found PC event..."<<endl;
    cnt.goodContainmentCounter++;
    
    Int_t diffAtEnd=-1;
    Int_t diffAtVtx=-1;
    if (!cuts.CheckTrkViews(s,&diffAtVtx,&diffAtEnd)) {
      cnt.badViewsCounter++;
      continue;
    }

    if (cuts.FilterBadTrackEnd(s,plInfo)){
      cnt.badTrackEndCounter++;
      continue;
    }
      
    Bool_t goodReco=reco.Reconstruct(s,plInfo);
    if (!goodReco){
      MSG("MeuAnalysis",logLevel)<<"Bad reco..."<<endl;
      cnt.badRecoCounter++;
      continue;
    }

    if (cuts.FilterBadEndGaps(s,plInfo)){//after reco only
      cnt.badEndGapsCounter++;
      continue;
    }
    
    MSG("MeuAnalysis",logLevel)<<"Recalibrating..."<<endl;
    this->RecalibrateAtNu(ev,s,plInfo);
    
    MSG("MeuAnalysis",logLevel)<<"CalcWindow..."<<endl;
    reco.CalcWindow(s,plInfo);
        
    //MSG("MeuAnalysis",logLevel)<<"CheckContainment..."<<endl;
    //fill the tree once event has been analysed
    if (s.WinSigMap<=0) {
      cnt.badWindowCounter++;
      continue;
    }
    
    if (!reco.CheckWinContainment(s,plInfo)){//nothing for FD
      cnt.badFidCounter++;
      continue;
    }
    
    //recalculate the containment!
    //s.DistToEdgeFid=0.35;//allow a slight shift in position
    //reco.CalcFCPC(s,plInfo);
    if (!s.PC){
      cnt.bad2ndContCounter++;
      static Bool_t print=false;
      if (print){
        MAXMSG("MeuAnalysis",Msg::kInfo,100)
          <<endl<<endl<<"Recalculated containment and not PC!"
          <<" PC="<<s.PC<<", FC="<<s.FC
          <<" entry="<<entry
          <<", snarl="<<s.Snarl<<", run="<<s.Run<<endl
          <<"Low En="<<s.MCLowEn<<" high="<<s.MCHighEn<<endl
          <<endl;
        MAXMSG("MeuAnalysis",Msg::kInfo,100)
          <<"vtx Pl="<<s.VtxPlane<<", end Pl="<<s.EndPlane
          <<", dVtx="<<s.VtxDistToEdge
          <<", dEnd="<<s.EndDistToEdge<<endl;
        
        for (map<Int_t,MeuHitInfo>::const_iterator it=
               plInfo.begin();
             it!=plInfo.end();++it){
          const MeuHitInfo& c=it->second;
          MAXMSG("MeuAnalysis",Msg::kInfo,500)
            <<"pl="<<c.Plane<<", z="<<c.Z
            <<", x="<<c.X<<", y="<<c.Y
            <<"  t="<<c.TPos<<", l="<<c.LPos<<endl;
        }
      }
      continue;
    }
    //txtFileNonPC<<entry<<endl;
      
    //filter out events with a low traceZ
    Float_t trace=999;
    Float_t traceZ=999;
    reco.CalcTrace(s,plInfo,&trace,&traceZ);
    if (traceZ<0.5){
      cnt.badTraceZCounter++;
      //eventsTxtFileTrace<<event<<endl;
      MAXMSG("MeuAnalysis",Msg::kDebug,1000)
        <<"trace="<<trace<<", traceZ="<<traceZ
        <<", trk.EndTraceZ="<<trk.EndTraceZ<<endl;
      continue;
    }

    reco.CalcStripDists(s,plInfo);
    if (cuts.FilterBadDistEndStrip(s,plInfo)){//nothing for FD
      cnt.badDistEndCounter++;
      continue;
    }
    
    //MUST BE GOOD!
    cnt.goodWindowCounter++;
      
    cuts.ExtractTruthInfo(ev,s,plInfo);
    cuts.FillSTSumRecoDetails(s,plInfo);
    histos.FillMeuHistos(s);

    //eventsTxtFile<<entry<<endl;
    MAXMSG("MeuAnalysis",Msg::kInfo,10)
      <<"Good window: run="
      <<ev.Run<<", snarl="<<ev.Snarl<<endl;
          
    if (s.MCLowEn>0.5*Munits::GeV) {
      //txtFileLowEn<<entry<<endl;
      MSG("MeuAnalysis",Msg::kInfo)
        <<"Low En="<<s.MCLowEn<<" high="<<s.MCHighEn
        <<", vtxPl="<<s.VtxPlane<<", endPl="<<s.EndPlane
        <<endl<<"entry="<<entry
        <<", snarl="<<s.Snarl<<", run="<<s.Run<<endl;
    }
            
    cuts.FillEnVsDist(s,plInfo);
    cuts.FillEnVsDist2(s,plInfo);
    //cuts.FillTimeHistos(ntp,evt,s);
    s.Temperature=this->GetTemperature(s);
    summaryOut.SummaryTreeFill(plInfo);
  }//end of for                                       
    

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

  MSG("MeuAnalysis",Msg::kInfo) 
    <<endl<<"*** Run Summary ***"<<endl;
  histos.PrintMeuValues();

  MSG("MeuAnalysis",Msg::kInfo)
    <<"Total Entries (snarls)="<<fEntries<<endl;
  cnt.PrintNtpSt();

  //finish the summary tree and output it
  summaryOut.SummaryTreeFinish();

  MSG("MeuAnalysis",Msg::kInfo) 
    <<" ** Finished MakeSummaryTreeWithAtNu method **"<<endl;
}

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

void MeuAnalysis::MakeSummaryTreeWithNtpStOneSnarl
(const NtpStRecord* pntp,const NtpBDLiteRecord* pntpBD) const
{

  //get a reference (coding preference)
  //don't get a reference to NtpBDLiteRecord because it doesn't exist
  //for MC
  const NtpStRecord& ntp=*pntp;

  //get candheader
  const RecCandHeader& rec=ntp.GetHeader();
  Int_t run=rec.GetRun();
  Int_t subrun=rec.GetSubRun();
  Int_t snarl=rec.GetSnarl();
  static Int_t lastRun=-1;
  static Int_t lastSubrun=-1;
  static Int_t runningTotalOfRuns=0;
  static Int_t runningTotalOfSubruns=0;
  if (run!=lastRun) {
    runningTotalOfRuns++;
    MSG("MeuAnalysis",Msg::kInfo)
      <<"Found new run="<<run
      <<", running total of runs used="<<runningTotalOfRuns<<endl;
  }
  lastRun=run;
  if (subrun!=lastSubrun) {
    runningTotalOfSubruns++;
    MSG("MeuAnalysis",Msg::kInfo)
      <<"Found new subrun="<<subrun
      <<", running total of subruns used="<<runningTotalOfSubruns<<endl;
  }
  lastSubrun=subrun;
  MAXMSG("MeuAnalysis",Msg::kInfo,1)
    <<"First entry: run="<<run<<", snarl="<<snarl<<endl;
  
  static const Int_t planeCut=8;
  //variable to turn on all the useful messages if required
  static const Msg::LogLevel_t logLevel=Msg::kDebug;

  //create static variables to prevent multiple constructions
  static Bool_t firstTime=true;

  //reconstruction object
  static const MeuReco reco;
  static const MeuCuts cuts;
  static const MeuHistos histos;
  static MeuCounter cnt;
  //count all the calls of this loop
  cnt.entriesAll++;
  
  static MeuSummaryWriter* sOut=0;

  static ofstream* peventsTxtFile=0;
  static ofstream* ptxtFileLowEn=0;
  static ofstream* ptxtFileNonPC=0;

  if (firstTime) {
    //set the control variable so code not run again
    firstTime=false;
    
    //create the writer
    sOut=new MeuSummaryWriter();
    sOut->SummaryTreeSetup(run,subrun);
  
    //this stores the objects' pointers for access at job end
    this->StoreOrFinishSummaryTree(sOut,&cnt,false);

    //printout the first element from the tree
    ntp.Print(cout);

    string temps="myeventsLowEn";
    if (run>10000) temps+="MC";
    peventsTxtFile=this->OpenTxtFile(100,"myevents");
    ptxtFileLowEn=this->OpenTxtFile(100,temps.c_str());
    ptxtFileNonPC=this->OpenTxtFile(100,"myeventsNonPC");
  }

  //get references (coding preference)
  MeuSummaryWriter& summaryOut=*sOut;
  ofstream& eventsTxtFile=*peventsTxtFile;
  ofstream& txtFileLowEn=*ptxtFileLowEn;
  ofstream& txtFileNonPC=*ptxtFileNonPC;
  

  MSG("MeuAnalysis",logLevel)
    <<"Entry="<<cnt.entry<<endl;
  
  TClonesArray& evtTca=(*ntp.evt);
  const Int_t numEvts=evtTca.GetEntriesFast();

  //count the number of snarls with zero events
  if (numEvts<=0) cnt.zeroEvtsCounter++;

  //loop over the events
  for (Int_t ievt=0;ievt<numEvts;++ievt) {
    const NtpSREvent* pevt=
      dynamic_cast<NtpSREvent*>(evtTca[ievt]);
    const NtpSREvent& evt=(*pevt);

    cnt.evtCounter++;
      
    //ensure there is only 1 track in the event
    if (evt.ntrack!=1) continue;
    cnt.oneTrackCounter++;
      
    //filter out events occuring in a slice with more than 1 event
    if (cuts.FilterBadEvtPerSlc(ntp,evt.slc,cnt.entry)) continue;
    cnt.oneEvtPerSlcCounter++;

    TClonesArray& trkTca=(*ntp.trk);
    const NtpSRTrack* ptrk=
      dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[0]]);
    const NtpSRTrack& trk=(*ptrk);

    //cut on number of planes
    if (sqrt(pow(1.*trk.plane.beg-
                 trk.plane.end,2))<planeCut) continue;
    cnt.longTrackCounter++;
      
    //MeuSummary s;
    MeuSummary& s=summaryOut.GetMeuSummaryToFill();
    s.Event=cnt.entry;
    cuts.FillSTSumDetails(ntp,s,ievt,evt.slc);
    cuts.SetSpecificCuts(s);
    //map to hold the info for each plane in the entry
    map<Int_t,MeuHitInfo> plInfo;      
      
    if (!cuts.IsCorrectTrigSrc(s)){
      continue;
    }
    cnt.correctTrigSrcCounter++;

    //cut out the LI events
    if (this->FilterLI(ntp,s)){
      //eventsTxtFileLI<<event<<endl;
      continue;
    }
    cnt.notLICounter++;

    if (pntpBD) cuts.GetBDSelectSpillInfo(*pntpBD,s,cnt.entry,
                                          cnt.goodSpillCounter);
    
    //print out info for specific snarl(s)
    if (s.Snarl==-1){
      MSG("MeuAnalysis",Msg::kInfo)
        <<"Entry="<<cnt.entry<<", run="<<run<<", snarl="<<snarl
        <<", trk.nstrips="<<trk.nstrip<<endl;
      cuts.PrintNtpSt(ntp);
    }

    if (cuts.FilterBadTrkTimes(ntp,s,evt)) continue;
    cnt.closeTimesCounter++;
      
    cuts.ExtractPlInfo(ntp,s,plInfo,evt);
    if (plInfo.size()==0) {
      MAXMSG("MeuAnalysis",Msg::kWarning,20)
        <<"No data extracted for entry="<<cnt.entry<<endl;
      continue;
    }
    cnt.notEmptyEventCounter++;
      
    reco.CalcVtx(s,plInfo);
      
    //check if track stops in spectrometer
    if (s.Detector==Detector::kNear && s.EndPlane>120) continue;
    cnt.endPlaneSpectCounter++;

    MSG("MeuAnalysis",logLevel)<<"CalcFCPC..."<<endl;
    reco.CalcFCPC(s,plInfo);
      
    //recalculate the vtx of the tracks
    //coming in from the spectrometer
    MSG("MeuAnalysis",logLevel)<<"CalcVtxSpecialND..."<<endl;
    reco.CalcVtxSpecialND(s,plInfo);

    //cut on number of planes with potential new vertex
    if (sqrt(pow(1.*s.VtxPlane-s.EndPlane,2))<planeCut) {
      MAXMSG("MeuAnalysis",Msg::kDebug,100)
        <<"Cutting on track length, run="<<s.Run
        <<", subrun="<<s.SubRun<<", snarl="<<s.Snarl<<endl
        <<"s.VtxPlane="<<s.VtxPlane
        <<", s.EndPlane="<<s.EndPlane
        <<", trk.beg="<<trk.plane.beg
        <<", trk.end="<<trk.plane.end<<endl;
      //cuts.PrintNtpSt(ntp);
      continue;
    }
    cnt.longTrack2Counter++;
      
    Bool_t TG=!s.FC && !s.PC;
      
    if (s.PC) cnt.PCCounter++;
    else if (s.FC) cnt.FCCounter++;
    else if (TG) cnt.TGCounter++;
    if (s.FC && s.PC) MSG("MeuAnalysis",Msg::kError)
      <<"Both FC and PC!"<<endl;
      
    Bool_t awayFromCoil=cuts.IsAwayFromCoil(s,plInfo);
    if (s.PC && !awayFromCoil) cnt.coilHitCounter++;
    if (s.PC && s.SMBoth && awayFromCoil) cnt.bothSMHitCounter++;
      
    Bool_t goodContainment=s.PC && !s.SMBoth && awayFromCoil;
    if (!goodContainment) continue;
    cnt.goodContainmentCounter++;
    MSG("MeuAnalysis",logLevel)<<"Found Good PC event..."<<endl;
    
    Int_t diffAtEnd=-1;
    Int_t diffAtVtx=-1;
    if (!cuts.CheckTrkViews(s,&diffAtVtx,&diffAtEnd)) {
      cnt.badViewsCounter++;
      continue;
    }

    //check that the hits beyond the end of the track are ok
    if (cuts.FilterBadTrackEnd(s,plInfo)){
      cnt.badTrackEndCounter++;
      continue;
    }

    Bool_t goodReco=reco.Reconstruct(s,plInfo);
    if (!goodReco){
      cnt.badRecoCounter++;
      MSG("MeuAnalysis",logLevel)<<"Bad reco..."<<endl;
      continue;
    }

    //check that there are no big gaps at the end of the track
    if (cuts.FilterBadEndGaps(s,plInfo)){//after reco only
      cnt.badEndGapsCounter++;
      continue;
    }

    //recalculate the containment (using MeuReco reconstruction)
    s.DistToEdgeFid=0.35;//allow a slight shift in position
    reco.CalcFCPC(s,plInfo);
          
    //re-check containment (using MeuReco reconstruction)
    if (!s.PC){
      cnt.bad2ndContCounter++;
      static Bool_t print=false;
      if (print){
        MAXMSG("MeuAnalysis",Msg::kInfo,100)
          <<endl<<endl<<"Recalculated containment and not PC!"
          <<" PC="<<s.PC<<", FC="<<s.FC
          <<" entry="<<cnt.entry
          <<", snarl="<<s.Snarl<<", run="<<s.Run<<endl
          <<"Low En="<<s.MCLowEn<<" high="<<s.MCHighEn<<endl
          <<endl;
        MAXMSG("MeuAnalysis",Msg::kInfo,100)
          <<"vtx Pl="<<s.VtxPlane<<", end Pl="<<s.EndPlane
          <<", dVtx="<<s.VtxDistToEdge
          <<", dEnd="<<s.EndDistToEdge<<endl;
          
        for (map<Int_t,MeuHitInfo>::const_iterator it=
               plInfo.begin();
             it!=plInfo.end();++it){
          const MeuHitInfo& c=it->second;
          MAXMSG("MeuAnalysis",Msg::kInfo,500)
            <<"pl="<<c.Plane<<", z="<<c.Z
            <<", x="<<c.X<<", y="<<c.Y
            <<"  t="<<c.TPos<<", l="<<c.LPos<<endl;
        }          
      }
      txtFileNonPC<<cnt.entry<<endl;
      continue;
    }

    MSG("MeuAnalysis",logLevel)<<"Recalibrating..."<<endl;
    this->Recalibrate(ntp,s,plInfo,evt);
          
    MSG("MeuAnalysis",logLevel)<<"CalcWindow..."<<endl;
    reco.CalcWindow(s,plInfo);
      
    MSG("MeuAnalysis",logLevel)<<"CheckContainment..."<<endl;
    //check that a track window was found
    if (s.WinSigMap<=0) {
      cnt.badWindowCounter++;
      continue;
    }

    //check that the track window is contained within the fid vol
    if (!reco.CheckWinContainment(s,plInfo)){//nothing for FD
      cnt.badFidCounter++;
      continue;
    }

    //filter out events with a low traceZ
    Float_t trace=999;
    Float_t traceZ=999;
    reco.CalcTrace(s,plInfo,&trace,&traceZ);//only FD now
    if (traceZ<0.5){
      cnt.badTraceZCounter++;
      //eventsTxtFileTrace<<event<<endl;
      MAXMSG("MeuAnalysis",Msg::kDebug,1000)
        <<"trace="<<trace<<", traceZ="<<traceZ<<endl;
      continue;
    }
      
    //check that all the hits in the track window are >20 cm from
    //the end of the strips
    reco.CalcStripDists(s,plInfo);
    if (cuts.FilterBadDistEndStrip(s,plInfo)){//requires window
      cnt.badDistEndCounter++;
      continue;
    }
      
    //if beam event check that the track window is away from 
    //the possible hadronic shower
    if (cuts.FilterBadShwDist(ntp,s,evt)){//only for ND spill data
      cnt.badShwDistCounter++;
      continue;
    }
      
    //MUST BE GOOD!
    cnt.goodWindowCounter++;
            
    cuts.ExtractTruthInfo(ntp,s,plInfo);
    cuts.FillSTSumRecoDetails(s,plInfo);
    histos.FillMeuHistos(s);
    if (pntpBD) cuts.FillBeamMonDetails(*pntpBD,s);
    
    eventsTxtFile<<cnt.entry<<endl;
    MAXMSG("MeuAnalysis",Msg::kInfo,10)
      <<"Good track window: run="<<run<<", snarl="<<snarl
      <<", evt="<<ievt+1<<"/"<<numEvts<<endl;
      
    if (s.MCLowEn>0.5*Munits::GeV) {
      txtFileLowEn<<cnt.entry<<endl;
      MSG("MeuAnalysis",Msg::kInfo)
        <<"Low En="<<s.MCLowEn<<" high="<<s.MCHighEn
        <<", vtxPl="<<s.VtxPlane<<", endPl="<<s.EndPlane
        <<endl<<"entry="<<cnt.entry
        <<", snarl="<<s.Snarl<<", run="<<s.Run<<endl;
    }
      
    cuts.FillEnVsDist(s,plInfo);
    cuts.FillEnVsDist2(s,plInfo);
    cuts.FillTimeHistos(ntp,evt,s);//a loop over lots of strips...
    s.Temperature=this->GetTemperature(s);
    summaryOut.SummaryTreeFill(plInfo);
  }//end of loop over events

  //increment ntuple entry (snarl) counter
  cnt.entry++;
}

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

void MeuAnalysis::StoreOrFinishSummaryTree(MeuSummaryWriter* psOut,
                                           MeuCounter* pcnt,
                                           Bool_t finish) const
{
  //this stores the pointers the first time the function is called
  static MeuCounter& cnt=*pcnt;
  static MeuSummaryWriter& summaryOut=*psOut;
  
  //this section runs methods on objects at the end of job
  if (finish){
    MSG("MeuAnalysis",Msg::kInfo) 
      <<endl<<"*** Run Summary ***"<<endl;
    MeuHistos histos;
    histos.PrintMeuValues();//this looks up TObjects by name 
    
    //print out the counter information
    cnt.PrintNtpSt();
    
    //finish the summary tree and output it
    summaryOut.SummaryTreeFinish();    
  }
}

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

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

  
  this->InitialiseLoopVariables();  
  
  //for(Int_t entry=55000;entry<fEntries;entry++){
  for(Int_t entry=0;entry<fEntries;entry++){
      
    this->SetLoopVariables(entry);

    //run the algorithm to calc meu values for events in snarl
    this->MakeSummaryTreeWithNtpStOneSnarl(fRec,fRecBD);

  }//end of for                                       
    

  MSG("MeuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl;
  
  //finish the job
  this->StoreOrFinishSummaryTree(NULL,NULL,1);
  
  MSG("MeuAnalysis",Msg::kInfo) 
    <<" ** Finished MakeSummaryTreeWithNtpSt method **"<<endl;
}

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

---