NtpMaker.cxx

Go to the documentation of this file.

#include "NtpMaker.h"
#include "AtmosCalculator.h"

#include "AtNuEvent/AtmosEvent.h"
#include "AtNuEvent/AtmosScintHit.h"
#include "AtNuEvent/AtmosStdHEP.h"
#include "AtNuEvent/AtmosDeadChip.h"
#include "AtNuEvent/AtmosShieldPlank.h"
#include "AtNuEvent/AtmosStrip.h"
#include "AtNuEvent/AtmosShower.h"
#include "AtNuEvent/AtmosTrack.h"

#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "JobControl/JobCModuleRegistry.h"

#include "RerootExodus/RerootExodus.h"
#include "REROOT_Classes/REROOT_NeuKin.h"
#include "Record/SimSnarlRecord.h"
#include "Record/SimSnarlHeader.h"
#include "Digitization/DigiSignal.h"
#include "Digitization/DigiScintHit.h"

#include "RawData/RawRecord.h"
#include "RawData/RawSnarlHeaderBlock.h"
#include "RawData/RawDigit.h"
#include "RawData/RawDigitDataBlock.h"
#include "RawData/RawDaqSnarlHeader.h"

#include "CandData/CandRecord.h"
#include "CandData/CandHeader.h"
#include "CandMorgue/CandDataQualityHandle.h"
#include "CandMorgue/CandDeadChipHandle.h"
#include "CandDigit/CandDeMuxDigitHandle.h"
#include "CandDigit/CandDeMuxDigitListHandle.h"
#include "RecoBase/CandStripHandle.h"
#include "RecoBase/CandStripListHandle.h"
#include "RecoBase/CandSliceHandle.h"
#include "RecoBase/CandSliceListHandle.h"
#include "RecoBase/CandShowerHandle.h"
#include "RecoBase/CandShowerListHandle.h"
#include "RecoBase/CandTrackListHandle.h"
#include "RecoBase/CandTrackHandle.h"
#include "RecoBase/CandFitTrackListHandle.h"
#include "RecoBase/CandFitTrackHandle.h"

#include "Calibrator/Calibrator.h"
#include "OnlineUtil/mdSpillData.h"
#include "SpillTiming/SpillTimeFinder.h"
#include "BeamDataUtil/BeamMonSpill.h"
#include "BeamDataUtil/BDSpillAccessor.h"
#include "BeamDataUtil/BMSpillAna.h"
#include "UgliGeometry/UgliGeomHandle.h"
#include "UgliGeometry/UgliStripHandle.h"
#include "Plex/PlexHandle.h"
#include "Plex/PlexStripEndId.h"
#include "Plex/PlexPlaneId.h"

#include "FarDetShieldPlankListHandle.h"
#include "FarDetShieldPlankHandle.h"
#include "FarDetStripListHandle.h"
#include "FarDetStripHandle.h"
#include "FarDetSliceListHandle.h"
#include "FarDetSliceHandle.h"
#include "FarDetEventListHandle.h"
#include "FarDetEventHandle.h"

#include "TClonesArray.h"
#include <vector>
#include "TParticle.h"

#include <cmath>

using std::vector;

JOBMODULE(NtpMaker,"NtpMaker","Make Ntps");

CVSID("$Id: NtpMaker.cxx,v 1.30 2009/02/13 15:03:01 blake Exp $");

NtpMaker::NtpMaker()
{
  MSG("NtpMaker",Msg::kDebug) << " *** NtpMaker::NtpMaker() *** " << endl;
  
  fFile = 0;
  fTree = 0;
  fNtpName = "ntp.root";
  fMaxFileSizeBytes = 0.0;
  fFileNumber = 1;
  fWriteStrips = true;
  fWriteScintHits = true;
  fWriteStdHEPs = true;
  fEvent = new AtmosEvent();
}

//======================================================================

NtpMaker::~NtpMaker()
{
  MSG("NtpMaker",Msg::kDebug) << " *** NtpMaker::~NtpMaker() *** " << endl;
}

//======================================================================

void NtpMaker::BeginJob()
{
  MSG("NtpMaker",Msg::kDebug) << " *** NtpMaker::BeginJob() *** " << endl;

}

//======================================================================

JobCResult NtpMaker::Ana(const MomNavigator* mom)
{
  MSG("NtpMaker",Msg::kDebug) << " *** NtpMaker::Ana(...) *** " << endl;

  JobCResult result(JobCResult::kPassed);  
  Int_t simflag=-1;
  Int_t run=-1, subrun=-1, snarl=-1;  
  Int_t runtype=-1, trigsrc=-1, timeframe=-1;
  Int_t date=-1, time=-1, microsec=-1, nanosec=-1, unixtime=-1;
  Int_t ndigits=-1;
  fEvent->Reset();

  // Process the SimSnarlRecord
  SimSnarlRecord* simrec = dynamic_cast<SimSnarlRecord *>(mom->GetFragment("SimSnarlRecord"));
  if(simrec)
    {
      MSG("NtpMaker",Msg::kDebug) << " *** PROCESS SIMSNARLRECORD *** " << endl;
      
      if(simflag<0) simflag = simrec->GetVldContext()->GetSimFlag();

      if(run<0) run=RerootExodus::GetRunNo();
      if(snarl<0) snarl=RerootExodus::GetEventNo();
      
      const SimSnarlHeader* hdr = dynamic_cast<const SimSnarlHeader*>(simrec->GetSimSnarlHeader());
      if(hdr)
        {
          if(run<0) run = hdr->GetRun();
          if(subrun<0) subrun = hdr->GetSubRun();
          if(snarl<0) snarl = hdr->GetSnarl();
          if(runtype<0) runtype = hdr->GetRunType();
        }
      
      this->FillMCInfo(simrec);  
    }

  // Process the RawRecord
  RawRecord* rawrec = dynamic_cast<RawRecord*>(mom->GetFragment("RawRecord"));
  if(rawrec)
    {
      MSG("NtpMaker",Msg::kDebug) << " *** PROCESS RAWRECORD *** " << endl;

      if(simflag<0) simflag = rawrec->GetVldContext()->GetSimFlag();

      const RawDaqSnarlHeader* hdr = dynamic_cast<const RawDaqSnarlHeader*>(rawrec->GetRawHeader());
      if(hdr)
        {
      
          if(run<0) run = hdr->GetRun(); 
          if(subrun<0) subrun = hdr->GetSubRun();
          if(snarl<0) snarl = hdr->GetSnarl();
          if(runtype<0) runtype = hdr->GetRunType();
          if(trigsrc<0) trigsrc = hdr->GetTrigSrc();
          if(timeframe<0) timeframe = hdr->GetTimeFrameNum();
          if(ndigits<0) ndigits = hdr->GetNumRawDigits();
          if(unixtime<0) unixtime = ((VldTimeStamp)(hdr->GetVldContext().GetTimeStamp())).GetSec();
          if(date<0) date = (((VldTimeStamp)(hdr->GetVldContext().GetTimeStamp())).GetSec()-1059696000)/(3600*24);
          if(time<0) time = (((VldTimeStamp)(hdr->GetVldContext().GetTimeStamp())).GetSec()-1059696000)%(3600*24);
          if(nanosec<0) nanosec = (Int_t)(1.000*(((VldTimeStamp)(hdr->GetVldContext().GetTimeStamp())).GetNanoSec()));
          if(microsec<0) microsec = (Int_t)(0.001*(((VldTimeStamp)(hdr->GetVldContext().GetTimeStamp())).GetNanoSec()));
          if(fEvent->SpillInfo.spilltype<0) fEvent->SpillInfo.spilltype =  hdr->GetRemoteSpillType();
        }      
      this->FillRawInfo(rawrec);
    }

  // Process the CandRecord
  CandRecord* candrec = dynamic_cast<CandRecord*>(mom->GetFragment("CandRecord"));
  if(candrec)
    {
      MSG("NtpMaker",Msg::kDebug) << " *** PROCESS CANDRECORD *** " << endl;
      
      if(simflag<0) simflag = candrec->GetVldContext()->GetSimFlag();

      const CandHeader* hdr  = dynamic_cast<const CandHeader*>(candrec->GetHeader());
      if(hdr)
        {
          Calibrator& cal = Calibrator::Instance();
          cal.ReInitialise(hdr->GetVldContext());
          
          if(run<0) run = hdr->GetRun();
          if(snarl<0) snarl = hdr->GetSnarl();
          if(unixtime<0) unixtime = ((VldTimeStamp)(hdr->GetVldContext().GetTimeStamp())).GetSec();
          if(date<0) date = (((VldTimeStamp)(hdr->GetVldContext().GetTimeStamp())).GetSec()-1059696000)/(3600*24);
          if(time<0) time = (((VldTimeStamp)(hdr->GetVldContext().GetTimeStamp())).GetSec()-1059696000)%(3600*24);
          if(nanosec<0) nanosec = (Int_t)(1.000*(((VldTimeStamp)(hdr->GetVldContext().GetTimeStamp())).GetNanoSec()));
          if(microsec<0) microsec = (Int_t)(0.001*(((VldTimeStamp)(hdr->GetVldContext().GetTimeStamp())).GetNanoSec()));
        }
      
      this->FillCandInfo(candrec);
    }

  fEvent->Run = run;
  fEvent->SubRun = subrun;  
  fEvent->Snarl = snarl;            
  fEvent->RunType = runtype;          
  fEvent->TrigSrc = trigsrc;           
  fEvent->TimeFrame = timeframe;        
  fEvent->Date = date;   
  fEvent->Time = time;       
  fEvent->MicroSec = microsec;     
  fEvent->NanoSec = nanosec;      
  fEvent->UnixTime = unixtime;     
  fEvent->Ndigits = ndigits;
  fEvent->SimFlag = simflag;

  MSG("NtpMaker",Msg::kDebug) << " *** EVENT SUMMARY *** " << endl
                              << "   Run = " << fEvent->Run << endl
                              << "   SubRun = " << fEvent->SubRun << endl
                              << "   Snarl = " << fEvent->Snarl << endl
                              << "   RunType = " << fEvent->RunType << endl 
                              << "   TrigSrc = " << fEvent->TrigSrc << endl
                              << "   TimeFrame = " << fEvent->TimeFrame << endl 
                              << "   Date = " << fEvent->Date << endl
                              << "   Time = " << fEvent->Time << endl
                              << "   MicroSec = " << fEvent->MicroSec << endl
                              << "   NanoSec = " << fEvent->NanoSec << endl
                              << "   UnixTime = " << fEvent->UnixTime << endl
                              << "   Ndigits = " << fEvent->Ndigits << endl
                              << "   SimFlag = " << fEvent->SimFlag << endl;

  
  MSG("NtpMaker",Msg::kInfo) << " *** NtpMaker: " << fEvent->Run << "/" << fEvent->SubRun << " " << fEvent->Snarl << " DONE *** " << endl;

  //Check if user wants strips written to file.
  if (!fWriteStrips)
    {
      fEvent->NStrips = 0;
      fEvent->StripList->Clear();
    }

  //Check if user wants strips written to file.
  if (!fWriteScintHits)
    {
      fEvent->NScintHits = 0;
      fEvent->ScintHitList->Clear();
    }

  //Check if user wants strips written to file.
  if (!fWriteStdHEPs)
    {
      fEvent->NStdHEPs = 0;
      fEvent->StdHEPList->Clear();
    }

  // create output file
  if(!fFile)
    {
      MSG("NtpMaker",Msg::kInfo) << " CREATING OUTPUT FILE: " << fNtpName.Data() << endl;

      // create the output file
      TDirectory *tmpd = gDirectory;
      fFile = TFile::Open(fNtpName.Data(), "recreate");
      fTree = new TTree("ntp", "FarDet Atmos Events");
      fTree->SetDirectory(fFile);
      fTree->Branch("evt", "AtmosEvent", &fEvent);
      fTree->SetAutoSave(100);
      tmpd->cd();
    }

  // re-create output file (at maximum file size)
  if( fFile
   && fMaxFileSizeBytes>0 
   && fFile->GetBytesWritten()>fMaxFileSizeBytes )
    {
      MSG("NtpMaker",Msg::kWarning) << " FILE SIZE EXCEEDS LIMIT (" << 1.0e-6*fFile->GetBytesWritten() << ">" << 1.0e-6*fMaxFileSizeBytes << ")" << endl;
      MSG("NtpMaker",Msg::kWarning) << "  ... SWITCHING TO NEW FILE" << endl;

      // write out the current file
      TDirectory *tmpd = gDirectory;
      fFile->cd();
      fFile->Write();
      fFile->Close();
      gDirectory = tmpd;

      // adjust the file name
      fFileNumber++;
      TDirectory *tmpd2 = gDirectory;
      TString filename = fNtpName;
      TString tempfilename = ".";
      tempfilename += fFileNumber;

      if( filename.EndsWith(".root") ){
        filename.Insert(filename.Length()-5,tempfilename);
      }
      else{
        filename.Append(tempfilename);
      }

      MSG("NtpMaker",Msg::kWarning) << "  ... CREATING NEW OUTPUT FILE: " << filename.Data() << endl;

      // create the new output file
      fFile = TFile::Open(filename.Data(), "recreate");
      fTree = new TTree("ntp", "FarDet Atmos Events");
      fTree->SetDirectory(fFile);
      fTree->Branch("evt", "AtmosEvent", &fEvent);
      fTree->SetAutoSave(100);
      tmpd2->cd();   
    }

  // check that both output file and tree exist
  if(!fFile)
    {
      MSG("NtpMaker",Msg::kError) << " Can't find output file: " << fNtpName.Data() << endl;
      MSG("NtpMaker",Msg::kFatal) << " Aaaggghhh... " << endl;
      assert(0);
    }
  if(!fTree)
    {
      MSG("NtpMaker",Msg::kError) << " Can't find output tree: " << fNtpName.Data() << endl;
      MSG("NtpMaker",Msg::kFatal) << " Aaaggghhh... " << endl;
      assert(0);
    }

  // write current event to file
  if(fFile)
    {
      MSG("NtpMaker",Msg::kDebug) << " writing event to file... " << endl;
      Int_t filesize(0),filesizeNew(0); 

      TDirectory* tmpd = gDirectory;
      filesize = fFile->GetBytesWritten();
      fFile->cd();
      fTree->Fill();
      filesizeNew = fFile->GetBytesWritten();
      gDirectory = tmpd;
      
      MSG("NtpMaker",Msg::kDebug) << "  ... bytes written: " << filesizeNew-filesize << endl;
      MSG("NtpMaker",Msg::kDebug) << "  ... current file size: " << filesizeNew << endl;
    }

  return result;
}

//======================================================================

const Registry& NtpMaker::DefaultConfig() const
{
  MSG("NtpMaker",Msg::kDebug) << " *** NtpMaker::DefaultConfig() *** " << endl;

  // Supply default configuration for the module
  static Registry r; 
  std::string name = this->GetName();
  name += ".config.default";
  r.SetName(name.c_str());
  r.UnLockKeys();
  r.UnLockValues();
  r.Set("NtpFileName","ntp.root");
  r.Set("MaxFileSize",0);
  r.Set("WriteStrips",1);
  r.Set("WriteScintHits",1);
  r.Set("WriteStdHEPs",1);
  r.LockValues();
  r.LockKeys();
  return r;
}

//======================================================================

void NtpMaker::Config(const Registry& r)
{
  MSG("NtpMaker",Msg::kDebug) << " *** NtpMaker::Config(...) *** " << endl;

  // Configure the module given the Registry r
  Int_t tmpint; const char* tmpchar = 0;
  if(r.Get("NtpFileName",tmpchar)) { fNtpName = tmpchar; }
  if(r.Get("MaxFileSize",tmpint)) { fMaxFileSizeBytes = 1.0e6*tmpint; }
  if(r.Get("WriteStrips",tmpint)) { fWriteStrips = (Bool_t)tmpint; }
  if(r.Get("WriteScintHits",tmpint)) { fWriteScintHits = (Bool_t)tmpint; }
  if(r.Get("WriteStdHEPs",tmpint)) { fWriteStdHEPs = (Bool_t)tmpint; }
}

//======================================================================

void NtpMaker::EndJob()
{
  MSG("NtpMaker",Msg::kDebug) << " *** NtpMaker::EndJob(...) *** " << endl;
  if( !fFile )
    {
      //Create the output file if it doesn't already exist.
      MSG("NtpMaker",Msg::kWarning) << " *** NTPMAKER: CREATING EMPTY OUTPUT FILE *** " << endl;
      TDirectory *tmpd = gDirectory;
      fFile = TFile::Open(fNtpName.Data(), "recreate");
      fTree = new TTree("ntp", "FarDet Atmos Events");
      fTree->SetDirectory(fFile);
      fTree->Branch("evt", "AtmosEvent", &fEvent);
      fTree->SetAutoSave(100);
      tmpd->cd();
    }
  // write out file
  if( fFile )
    {
      TDirectory* tmpd = gDirectory;
      fFile->cd();
      MSG("NtpMaker",Msg::kDebug) << " WRITING FILE: " << fFile->GetBytesWritten() << " bytes (before)" << endl;
      fFile->Write();
      MSG("NtpMaker",Msg::kDebug) << " WRITING FILE: " << fFile->GetBytesWritten() << " bytes (after)" << endl;
      fFile->Close();
      gDirectory = tmpd;
    }
}

//======================================================================

void NtpMaker::FillCandInfo(CandRecord* candrec)
{
  MSG("NtpMaker",Msg::kDebug) << " *** NtpMaker::FillCandInfo(...) *** " << endl;

  // Fill information from the candidates
  TIter candbit(candrec->GetCandHandleIter());
  while(TObject* tob = (TObject*)(candbit()))
    {
      MSG("NtpMaker",Msg::kDebug) << tob->GetName() << endl;
  
      if(tob->InheritsFrom("CandDataQualityHandle"))
        {
          FillDataQualityInfo(tob); continue;
        }

      if(tob->InheritsFrom("CandDeMuxDigitListHandle"))
        {
          FillDeMuxInfo(tob); continue;
        }
      if(tob->InheritsFrom("FarDetShieldPlankListHandle"))
        {
          FillShieldInfo(tob); continue;
        }
  
      if(tob->InheritsFrom("FarDetStripListHandle"))
        {
          FillStripInfo(tob); continue;
        }
 
      if(tob->InheritsFrom("FarDetSliceListHandle"))
        {
          FillSliceInfo(tob); continue;
        }

      if(tob->InheritsFrom("CandFitTrackListHandle"))
        {
          FillFitTrackInfo(tob); continue;
        }      

      if(tob->InheritsFrom("CandTrackListHandle"))
        {
          FillTrackInfo(tob); continue;
        }
  
      if(tob->InheritsFrom("CandShowerListHandle"))
        {
          FillShowerInfo(tob); continue;
        }
  
      if(tob->InheritsFrom("FarDetEventListHandle"))
        {
          FillEventInfo(tob); continue;
        }
      MSG("NtpMaker",Msg::kDebug) << "   NtpMaker has no use for " << tob->GetName() <<endl;
    }

  // now we've got all the information from the candidates, we can calculate 
  // the atmospheric neutrino analysis variables using the AtmosCalculator tool

  MSG("NtpMaker",Msg::kDebug) << " Running AtmosCalculator..." << endl;
  AtmosCalculator AtmosCalc;
  TClonesArray& MyFitList = *(fEvent->TrackList);
  TClonesArray* MyStrpList = (fEvent->StripList);
  TClonesArray& MyShowerList = *(fEvent->ShowerList);

  int trklimit = MyFitList.GetLast()+1;
  for(int i=0; i<trklimit; ++i)
    {
      AtmosTrack* MyFit = (AtmosTrack*)(MyFitList[i]);
      AtmosCalc.TrackProperties(MyFit,MyStrpList);
    }

  int shwlimit = MyShowerList.GetLast()+1;
  for(int i=0; i<shwlimit; ++i)
    {
      AtmosShower* MyShower = (AtmosShower*)(MyShowerList[i]);
      AtmosCalc.ShowerProperties(MyShower,MyStrpList);
    }

  AtmosCalc.EventProperties(fEvent,MyStrpList);

  MSG("NtpMaker",Msg::kDebug) << " ... AtmosCalculator Finished" << endl;

  return;
}

//======================================================================

void NtpMaker::FillDataQualityInfo(TObject* tob)
{
  // Fill data quality information
  // (uses CandDataQualityHandle from CandMorgue package)
  CandDataQualityHandle* Data = (CandDataQualityHandle*)(tob);
  MSG("NtpMaker",Msg::kDebug) << " ... found CandDataQualityHandle " << endl;
  TClonesArray& MyDeadChipList = *(fEvent->DeadChipList);

  fEvent->DataInfo.TriggerSource =     Data->GetTriggerSource();
  fEvent->DataInfo.TriggerTime =       Data->GetTriggerTime();
  fEvent->DataInfo.ErrorCode =         Data->GetErrorCode();
  fEvent->DataInfo.CrateMask =         Data->GetCrateMask();
  fEvent->DataInfo.PreTriggerDigits =  Data->GetPreTriggerDigits();
  fEvent->DataInfo.PostTriggerDigits = Data->GetPostTriggerDigits();
  fEvent->DataInfo.SnarlMultiplicity = Data->GetSnarlMultiplicity();
  fEvent->DataInfo.SpillStatus =       Data->GetSpillStatus();
  fEvent->DataInfo.SpillType =         Data->GetSpillType();
  fEvent->DataInfo.SpillTimeError =    Data->GetSpillTimeError();
  fEvent->DataInfo.LiTrigger =         Data->GetLiTrigger();
  fEvent->DataInfo.LiTime =            Data->GetLiTime();
  fEvent->DataInfo.LiSubtractedTime =  Data->GetLiSubtractedTime();
  fEvent->DataInfo.LiRelativeTime =    Data->GetLiRelativeTime();
  fEvent->DataInfo.LiCalibPoint =      Data->GetLiCalibPoint();
  fEvent->DataInfo.LiCalibType =       Data->GetLiCalibType();
  fEvent->DataInfo.LiPulserBox =       Data->GetLiPulserBox();
  fEvent->DataInfo.LiPulserLed =       Data->GetLiPulserLed();
  fEvent->DataInfo.LiPulseHeight =     Data->GetLiPulseHeight();
  fEvent->DataInfo.LiPulseWidth =      Data->GetLiPulseWidth();
  fEvent->DataInfo.ColdChips =         Data->GetColdChips();
  fEvent->DataInfo.HotChips =          Data->GetHotChips();  
  fEvent->DataInfo.BusyChips =         Data->GetBusyChips();
  fEvent->DataInfo.ReadoutErrors =     Data->GetReadoutErrors();
  fEvent->DataInfo.DataQualityBitMap = (Int_t)Data->GetDataQuality();

  //Get Details of Dead Chips
  PlexHandle plex(*(Data->GetVldContext()));
  TIter DataItr(Data->GetDaughterIterator());
  while(CandDeadChipHandle* DeadChip = dynamic_cast<CandDeadChipHandle*>(DataItr()))
    {
      AtmosDeadChip* MyDeadChip = new( MyDeadChipList[(fEvent->NDeadChips)++] ) AtmosDeadChip();
      RawChannelId rcid = DeadChip->GetChannelId();
      
      MyDeadChip->Crate = rcid.GetCrate();
      MyDeadChip->Varc = rcid.GetVarcId();
      MyDeadChip->Vmm = rcid.GetVmm();
      MyDeadChip->VaAdc = rcid.GetVaAdcSel();
      MyDeadChip->VaChip = rcid.GetVaChip();
      RawChannelId rcidCheck = RawChannelId(Detector::kFar,ElecType::kVA,
                                            MyDeadChip->Crate,  MyDeadChip->Varc,
                                            MyDeadChip->Vmm, MyDeadChip->VaAdc, 
                                            MyDeadChip->VaChip, 10);

      if( plex.GetSEIdAltL(rcidCheck).IsValid() && plex.GetSEIdAltL(rcidCheck).GetPlane()>=0 && plex.GetSEIdAltL(rcidCheck).GetPlane()<1000 )
        {
          MyDeadChip->InReadout = 1;

          // Veto Shield
          if(plex.GetSEIdAltL(rcidCheck).IsVetoShield())
            {
              MyDeadChip->Shield = 1;
            }

          // Main Detector
          if(!plex.GetSEIdAltL(rcidCheck).IsVetoShield())
            {//Only set the plane numbers if we're in the main detector 
          
              RawChannelId rcidLow = RawChannelId(Detector::kFar,ElecType::kVA,
                                                  MyDeadChip->Crate,  MyDeadChip->Varc,
                                                  MyDeadChip->Vmm, MyDeadChip->VaAdc, 
                                                  0,2);

              RawChannelId rcidHigh = RawChannelId(Detector::kFar,ElecType::kVA,
                                                   MyDeadChip->Crate,  MyDeadChip->Varc,
                                                   MyDeadChip->Vmm, MyDeadChip->VaAdc, 
                                                   1,17);

              if(MyDeadChip->VaChip!=1) MyDeadChip->Plane[0] = plex.GetSEIdAltL(rcidLow).GetPlane();
              if(MyDeadChip->VaChip!=0) MyDeadChip->Plane[1] = plex.GetSEIdAltL(rcidHigh).GetPlane();
            }
        }
      
      MyDeadChip->ErrorCode = DeadChip->GetErrorCode();
      MyDeadChip->Status = (Int_t)DeadChip->GetChipStatus();
    }

  return;
}

//======================================================================

void NtpMaker::FillDeMuxInfo(TObject* tob)
{
  // Get output from demuxing 
  CandDeMuxDigitListHandle* DigitList = (CandDeMuxDigitListHandle*)(tob);
  MSG("NtpMaker",Msg::kDebug) << " ... found CandDeMuxDigitListHandle " << endl;
  fEvent->RecoInfo.Ndigits = DigitList->GetNDaughters();
  fEvent->RecoInfo.DeMuxFlagWord = DigitList->GetDeMuxDigitListFlagWord();
  fEvent->RecoInfo.MultipleMuonFlag = (Int_t)((DigitList->GetDeMuxDigitListFlagWord())&(CandDeMuxDigitList::kMultipleMuonEvent));
  fEvent->RecoInfo.AbsTime = DigitList->GetAbsTime();

  fEvent->RecoInfo.Nplndigits = 0;
  TIter digitr(DigitList->GetDaughterIterator());
  while(CandDigitHandle* cdh = dynamic_cast<CandDigitHandle*>(digitr()))
    {
      if( cdh )
        {
          if( !cdh->GetPlexSEIdAltL().IsVetoShield() 
           && cdh->GetPlexSEIdAltL().GetPlane()>=0
           && cdh->GetPlexSEIdAltL().GetPlane()<500 
           && cdh->GetSubtractedTime()>=0 )
            {
              ++fEvent->RecoInfo.Nplndigits;
            }
        }
    }


  return;
}

//======================================================================

void NtpMaker::FillEventInfo(TObject* tob)
{
  FarDetEventListHandle* EventList = (FarDetEventListHandle*)(tob);
  MSG("NtpMaker",Msg::kDebug) << " ... found FarDetEventListHandle " << endl;

  if(EventList->GetNDaughters()<1) 
    {
      MSG("NtpMaker",Msg::kWarning) << " ... FarDetEventListHandle has no daughters. " << endl;
      MSG("NtpMaker",Msg::kWarning) << " ... failed to find FarDetEventHandle " << endl;
      return;
    }
  // for atmospheric neutrinos in the Far Detector, assume only 1 event per snarl
  const FarDetEventHandle* Event = dynamic_cast<const FarDetEventHandle*>(EventList->GetDaughter(0));
  
  if(Event==0) 
    {
      MSG("NtpMaker",Msg::kWarning) << " ... failed to find FarDetEventHandle " << endl;
      return;
    }
  MSG("NtpMaker",Msg::kDebug) << " ... found FarDetEventHandle " << endl;

  fEvent->RecoInfo.EMCharge = Event->GetEMCharge();

  fEvent->RecoInfo.MuReco = Event->GetMuReco();
  fEvent->RecoInfo.MuMomentumRange = Event->GetMuMomentumRange();
  fEvent->RecoInfo.MuMomentumCurve = Event->GetMuMomentumCurve();
  fEvent->RecoInfo.MuDirCosU = Event->GetMuDirCosU();
  fEvent->RecoInfo.MuDirCosV = Event->GetMuDirCosV();
  fEvent->RecoInfo.MuDirCosZ = Event->GetMuDirCosZ();
  fEvent->RecoInfo.MuDirCosX = -999;
  fEvent->RecoInfo.MuDirCosY = -999;
  if( fEvent->RecoInfo.MuReco )
    {
      fEvent->RecoInfo.MuDirCosX = 0.7071*(Event->GetMuDirCosU()-Event->GetMuDirCosV());
      fEvent->RecoInfo.MuDirCosY = 0.7071*(Event->GetMuDirCosU()+Event->GetMuDirCosV());
    }
  
  fEvent->RecoInfo.ShwReco = Event->GetShwReco();
  fEvent->RecoInfo.ShwMipsLinear = Event->GetShwMipsLinear();
  fEvent->RecoInfo.ShwMipsDeweighted = Event->GetShwMipsDeweighted();
  fEvent->RecoInfo.ShwEnergyGeVLinear = Event->GetShwEnergyGeVLinear(); 
  fEvent->RecoInfo.ShwEnergyGeVDeweighted = Event->GetShwEnergyGeVDeweighted(); 
  fEvent->RecoInfo.ShwDirCosU = Event->GetShwDirCosU();
  fEvent->RecoInfo.ShwDirCosV = Event->GetShwDirCosV();
  fEvent->RecoInfo.ShwDirCosZ = Event->GetShwDirCosZ();
  fEvent->RecoInfo.ShwDirCosX = -999;
  fEvent->RecoInfo.ShwDirCosY = -999;
  if( fEvent->RecoInfo.ShwReco )
    {
      fEvent->RecoInfo.ShwDirCosX = 0.7071*(Event->GetShwDirCosU()-Event->GetShwDirCosV());
      fEvent->RecoInfo.ShwDirCosY = 0.7071*(Event->GetShwDirCosU()+Event->GetShwDirCosV());
    }

  fEvent->RecoInfo.VtxTime = Event->GetVtxTime();
  fEvent->RecoInfo.VtxU = Event->GetVtxU();
  fEvent->RecoInfo.VtxV = Event->GetVtxV();
  fEvent->RecoInfo.VtxX = 0.7071*(Event->GetVtxU()-Event->GetVtxV());
  fEvent->RecoInfo.VtxY = 0.7071*(Event->GetVtxU()+Event->GetVtxV());
  fEvent->RecoInfo.VtxZ = Event->GetVtxZ();
  fEvent->RecoInfo.VtxPlane = Event->GetVtxPlane();

  fEvent->RecoInfo.MinPlane = Event->GetMinPlane();
  fEvent->RecoInfo.MaxPlane = Event->GetMaxPlane();
  fEvent->RecoInfo.MinZ = Event->GetMinZ();
  fEvent->RecoInfo.MaxZ = Event->GetMaxZ();

  return;

}

//======================================================================

void NtpMaker::FillFitTrackInfo(TObject* tob)
{
  CandFitTrackListHandle* FitList = (CandFitTrackListHandle*)(tob);
  MSG("NtpMaker",Msg::kDebug) << " ... found CandFitTrackListHandle containing "<<FitList->GetNDaughters()<<" tracks."<< endl;
  float nstripcalib(0.0),stripsintrack(0.0);
  Calibrator& cal = Calibrator::Instance(); 
      
  //Reset the track flags
  TClonesArray* MyStrpList = (TClonesArray*)(fEvent->StripList);
  for(Int_t myint=0;myint<fEvent->NStrips;myint++)
    {
      AtmosStrip& MyStrp = *((AtmosStrip*)(MyStrpList->At(myint)));
      MyStrp.Trk = 0;
      MyStrp.dS = -999.9;
    }

  // Set information from CandFitTrackHandles
  // assume same number of CandTracks and CandFitTracks
  // (the CandFitTracks entirely overwrite the CandTracks)
  TClonesArray* MyFitList = (TClonesArray*)(fEvent->TrackList);  
  TIter FitItr(FitList->GetDaughterIterator());
  fEvent->NFits = 0; 
  while(CandFitTrackHandle* Fit = (CandFitTrackHandle*)(FitItr()))
    {
      if( fEvent->NFits < fEvent->NTracks )
        {
          AtmosTrack* MyFit = (AtmosTrack*)(MyFitList->At(fEvent->NFits++));
          MyFit->VtxPlane = Fit->GetVtxPlane();
          MyFit->VtxTime = Fit->GetVtxT();
          MyFit->VtxU = Fit->GetVtxU();
          MyFit->VtxV = Fit->GetVtxV();
          MyFit->VtxR = sqrt((MyFit->VtxU*MyFit->VtxU)+(MyFit->VtxV*MyFit->VtxV));
          MyFit->VtxX = 0.7071*(Fit->GetVtxU()-Fit->GetVtxV());
          MyFit->VtxY = 0.7071*(Fit->GetVtxU()+Fit->GetVtxV());
          MyFit->VtxZ = Fit->GetVtxZ();
          MyFit->VtxTrace = Fit->GetVtxTrace();
          MyFit->VtxTraceZ = Fit->GetVtxTraceZ();
          MyFit->VtxDirCosU = Fit->GetDirCosU();
          MyFit->VtxDirCosV = Fit->GetDirCosV();
          MyFit->VtxDirCosX = 0.7071*(Fit->GetDirCosU()-Fit->GetDirCosV());
          MyFit->VtxDirCosY = 0.7071*(Fit->GetDirCosU()+Fit->GetDirCosV());
          MyFit->VtxDirCosZ = Fit->GetDirCosZ();
          MyFit->EndPlane = Fit->GetEndPlane();
          MyFit->EndTime = Fit->GetEndT();
          MyFit->EndU = Fit->GetEndU();
          MyFit->EndV = Fit->GetEndV();
          MyFit->EndR = sqrt((MyFit->EndU*MyFit->EndU)+(MyFit->EndV*MyFit->EndV)); 
          MyFit->EndX = 0.7071*(Fit->GetEndU()-Fit->GetEndV());
          MyFit->EndY = 0.7071*(Fit->GetEndU()+Fit->GetEndV());
          MyFit->EndZ = Fit->GetEndZ();
          MyFit->EndTrace = Fit->GetEndTrace();
          MyFit->EndTraceZ = Fit->GetEndTraceZ();
          MyFit->EndDirCosU = Fit->GetEndDirCosU();
          MyFit->EndDirCosV = Fit->GetEndDirCosV();
          MyFit->EndDirCosX = 0.7071*(Fit->GetEndDirCosU()-Fit->GetEndDirCosV());
          MyFit->EndDirCosY = 0.7071*(Fit->GetEndDirCosU()+Fit->GetEndDirCosV());
          MyFit->EndDirCosZ = Fit->GetEndDirCosZ();
          MyFit->RangeGCM2 = Fit->GetRange(Fit->GetVtxPlane());//NB this will lose the extrapolation at vtx and end
          MyFit->RangeMetres = Fit->GetdS(Fit->GetVtxPlane());//NB this will lose the extrapolation at vtx and end
          MyFit->Nstrips = Fit->GetNStrip(); // AtNu Trk->GetNStrips();
          MyFit->Nplanes = abs(MyFit->EndPlane-MyFit->VtxPlane)+1;//Trk->GetNPlane() // AtNu Trk->GetNPlanes();

          MyFit->TimingFitChi2 = Fit->GetTimeFitChi2();
          MyFit->TimingFitNdf = Fit->GetNTimeFitDigit()-2;
          MyFit->TimeSlope = Fit->GetTimeSlope();
          MyFit->TimeOffset = Fit->GetTimeOffset(); 

          if(MyFit->EndPlane-MyFit->VtxPlane>0)
            {//Forward and Backward variables defined relative to increasing Z
              MyFit->VtxDirTimeFitRMS = Fit->GetTimeForwardFitRMS();  
              MyFit->EndDirTimeFitRMS = Fit->GetTimeBackwardFitRMS(); 
              MyFit->VtxDirTimeFitNdf = Fit->GetTimeForwardFitNDOF(); 
              MyFit->EndDirTimeFitNdf = Fit->GetTimeBackwardFitNDOF();
            }
          else
            {
              MyFit->VtxDirTimeFitRMS = Fit->GetTimeBackwardFitRMS(); 
              MyFit->EndDirTimeFitRMS = Fit->GetTimeForwardFitRMS();  
              MyFit->VtxDirTimeFitNdf = Fit->GetTimeBackwardFitNDOF();
              MyFit->EndDirTimeFitNdf = Fit->GetTimeForwardFitNDOF(); 
            }

          MyFit->Momentum = Fit->GetMomentumCurve();
          MyFit->MomentumRange = Fit->GetMomentumRange();
          MyFit->MomentumCurve = Fit->GetMomentumCurve();

          MyFit->FitPass = Fit->GetPass();    
          MyFit->EMcharge = Fit->GetEMCharge();
          MyFit->QPvtxChi2 = Fit->GetChi2();
          MyFit->QPvtxNdf = Fit->GetNDOF();
          MyFit->QPvtx = Fit->GetPlaneQP(MyFit->VtxPlane);
          MyFit->QPvtxTweaked = (MyFit->Momentum!=0.0)?MyFit->EMcharge/MyFit->Momentum:0.0;
          MyFit->QPvtxErr = Fit->GetVtxQPError();

          MyFit->Index = (int)pow(2.,fEvent->NFits-1);

          // (don't) reset track strips
          //for(Int_t myint=0;myint<fEvent->NStrips;myint++)
          //  {
          //    AtmosStrip& MyStrp = *((AtmosStrip*)(MyStrpList->At(myint)));
          //    MyStrp.Trk = 0;
          //  }

          //set information for strips that are part of this track
          nstripcalib=(0.0); stripsintrack = Fit->GetNDaughters();
          TIter StrpItr(Fit->GetDaughterIterator());
          PlexStripEndId* seidP = 0;
          PlexStripEndId* seidN = 0;
          while(CandStripHandle* TrkStrp = (CandStripHandle*)(StrpItr()))
            {
              for(Int_t myint=0;myint<fEvent->NStrips;myint++)
                {
                  AtmosStrip& MyStrp = *((AtmosStrip*)(MyStrpList->At(myint)));

                  if( TrkStrp->GetPlane()==MyStrp.Plane && TrkStrp->GetStrip()==MyStrp.Strip )
                    {
                      if( MyStrp.Trk==0 && Fit->IsTPosValid(MyStrp.Plane) )
                        { // (using MyStrp.Trk==0 as primary track info should dominate)
                          if( MyStrp.View==0 )
                            {
                              MyStrp.L = Fit->GetV(MyStrp.Plane);
                              MyStrp.dS = Fit->GetdS(MyStrp.Plane);
                              MyStrp.GreenFibre[0] = MyStrp.HalfLength - MyStrp.L;
                              MyStrp.GreenFibre[1] = MyStrp.HalfLength + MyStrp.L;
                              seidN = new PlexStripEndId(Detector::kFar, MyStrp.Plane, MyStrp.Strip, StripEnd::kNegative, StripEnd::kWhole, PlaneView::kU, PlaneCoverage::kTotal);
                              seidP = new PlexStripEndId(Detector::kFar, MyStrp.Plane, MyStrp.Strip, StripEnd::kPositive, StripEnd::kWhole, PlaneView::kU, PlaneCoverage::kTotal);
                            }
                          if( MyStrp.View==1 )
                            {
                              MyStrp.L = Fit->GetU(MyStrp.Plane);
                              MyStrp.dS = Fit->GetdS(MyStrp.Plane);
                              MyStrp.GreenFibre[0] = MyStrp.HalfLength + MyStrp.L;
                              MyStrp.GreenFibre[1] = MyStrp.HalfLength - MyStrp.L;
                              seidN = new PlexStripEndId(Detector::kFar, MyStrp.Plane, MyStrp.Strip, StripEnd::kNegative, StripEnd::kWhole, PlaneView::kV, PlaneCoverage::kTotal);
                              seidP = new PlexStripEndId(Detector::kFar, MyStrp.Plane, MyStrp.Strip, StripEnd::kPositive, StripEnd::kWhole, PlaneView::kV, PlaneCoverage::kTotal);
                            }
                          double lsign = pow(-1.0, (1.0+double(MyStrp.View)) );
                          if(seidN && MyStrp.Sigcorr[0]>0.0)
                            {
                              MyStrp.Sigmap[0] = cal.GetAttenCorrected(MyStrp.Sigcorr[0], lsign*MyStrp.L, *(seidN) );
                              MyStrp.MIP[0] = cal.GetMIP(MyStrp.Sigmap[0], *(seidN) );
                            }
                          if(seidP && MyStrp.Sigcorr[1]>0.0)
                            {
                              MyStrp.Sigmap[1] = cal.GetAttenCorrected(MyStrp.Sigcorr[1], lsign*MyStrp.L, *(seidP) );
                              MyStrp.MIP[1] = cal.GetMIP(MyStrp.Sigmap[1], *(seidP) );
                            }
                          if(MyStrp.Sigcorr[0]>0.0 || MyStrp.Sigcorr[1]>0.0) nstripcalib+=1.0;
                        }
                      MyStrp.Trk += MyFit->Index;

                      if(seidP) delete seidP; seidP = 0;
                      if(seidN) delete seidN; seidN = 0;

                      break; // break out of loop
                             // faster - this should be ok as long as the
                             // same strip doesn't exist twice in a snarl
                    }
                } 
            }

          MSG("NtpMaker",Msg::kVerbose) <<"   "<<(int)stripsintrack<<" strips in the fitted track of which "<<(int)nstripcalib<<" were calibrated."<<endl<<flush;
        }
    }
  return;
}

//======================================================================

void NtpMaker::FillMCInfo(SimSnarlRecord* simrec)
{
  MSG("NtpMaker",Msg::kDebug) << " *** NtpMaker::FillMCInfo(...) *** " << endl;

  UgliGeomHandle ugh = (*(simrec->GetVldContext())); 
  const TObjArray arr(simrec->GetComponents());
  TIter iter(arr.MakeIterator());
  while(TObject* tob = (TObject*)(iter()))
    {
      MSG("NtpMaker",Msg::kDebug) << tob->GetName() << endl;

      if(tob->GetName()==TString("REROOT_NeuKin"))
        {
          REROOT_NeuKin* NeuKin = dynamic_cast<REROOT_NeuKin*>(tob);
          MSG("NtpMaker",Msg::kDebug) << " ... found REROOT_NeuKin " << endl;
          fEvent->MCInfo.IDnu = NeuKin->INu();
          fEvent->MCInfo.IDnunoosc = NeuKin->INuNoOsc();
          fEvent->MCInfo.IDact = NeuKin->IAction();
          fEvent->MCInfo.IDres = NeuKin->IResonance();
          fEvent->MCInfo.IDboson = NeuKin->IBoson();
          fEvent->MCInfo.IDtarget = NeuKin->ITg();
          fEvent->MCInfo.A = NeuKin->A();
          fEvent->MCInfo.Z = NeuKin->Z();
          fEvent->MCInfo.Xsection = NeuKin->Sigma();
          fEvent->MCInfo.EMfrac = NeuKin->EMFrac();
          fEvent->MCInfo.W2 = NeuKin->W2();
          fEvent->MCInfo.Q2 = NeuKin->Q2();
          fEvent->MCInfo.x = NeuKin->X();
          fEvent->MCInfo.y = NeuKin->Y();
          fEvent->MCInfo.PnuX = NeuKin->P4Neu()[0];
          fEvent->MCInfo.PnuY = NeuKin->P4Neu()[1];
          fEvent->MCInfo.PnuZ = NeuKin->P4Neu()[2];
          fEvent->MCInfo.Enu = NeuKin->P4Neu()[3];
          fEvent->MCInfo.PmuX = NeuKin->P4Mu1()[0];
          fEvent->MCInfo.PmuY = NeuKin->P4Mu1()[1];
          fEvent->MCInfo.PmuZ = NeuKin->P4Mu1()[2];
          fEvent->MCInfo.Emu = NeuKin->P4Mu1()[3];
          fEvent->MCInfo.PelX = NeuKin->P4El1()[0];
          fEvent->MCInfo.PelY = NeuKin->P4El1()[1];
          fEvent->MCInfo.PelZ = NeuKin->P4El1()[2];
          fEvent->MCInfo.Eel = NeuKin->P4El1()[3];
          fEvent->MCInfo.PhadX = NeuKin->P4Shw()[0];
          fEvent->MCInfo.PhadY = NeuKin->P4Shw()[1];
          fEvent->MCInfo.PhadZ = NeuKin->P4Shw()[2];
          fEvent->MCInfo.Ehad = NeuKin->P4Shw()[3];
          fEvent->MCInfo.PtargX = NeuKin->P4Tgt()[0];
          fEvent->MCInfo.PtargY = NeuKin->P4Tgt()[1];
          fEvent->MCInfo.PtargZ = NeuKin->P4Tgt()[2];
          fEvent->MCInfo.Etarg = NeuKin->P4Tgt()[3];
        }

      if(tob->GetName()==TString("StdHep"))
        {
          TClonesArray* ParticleList = (TClonesArray*)(tob);
          MSG("NtpMaker",Msg::kDebug) << " ... found StdHep " << endl;        
          TClonesArray& MyStdHEPList = *(fEvent->StdHEPList);

          TParticle* Particle = dynamic_cast<TParticle*>(ParticleList->At(0));
          fEvent->MCInfo.VtxU = 0.7071*(Particle->Vy()+Particle->Vx()); 
          fEvent->MCInfo.VtxV = 0.7071*(Particle->Vy()-Particle->Vx()); 
          fEvent->MCInfo.VtxX = Particle->Vx(); 
          fEvent->MCInfo.VtxY = Particle->Vy(); 
          fEvent->MCInfo.VtxZ = Particle->Vz();
          TIter ParticleItr(ParticleList->MakeIterator());
          while( const TParticle* StdHEP = dynamic_cast<const TParticle*>(ParticleItr()))
            {
              AtmosStdHEP* MyStdHEP = new( MyStdHEPList[(fEvent->NStdHEPs)++] ) AtmosStdHEP();
              MyStdHEP->Index     = fEvent->NStdHEPs - 1;

              MyStdHEP->Parent[0] = StdHEP->GetMother(0);
              MyStdHEP->Parent[1] = StdHEP->GetMother(1);
              MyStdHEP->Child[0]  = StdHEP->GetDaughter(0);
              MyStdHEP->Child[1]  = StdHEP->GetDaughter(1);
              MyStdHEP->Id        = StdHEP->GetPdgCode();
              MyStdHEP->Status    = StdHEP->GetStatusCode();
              MyStdHEP->Mass      = StdHEP->GetCalcMass();
              MyStdHEP->PX        = StdHEP->Px();
              MyStdHEP->PY        = StdHEP->Py();
              MyStdHEP->PZ        = StdHEP->Pz();
              MyStdHEP->E         = StdHEP->Energy();
              MyStdHEP->VtxX      = StdHEP->Vx();
              MyStdHEP->VtxY      = StdHEP->Vy();
              MyStdHEP->VtxX      = StdHEP->Vz();
              MyStdHEP->VtxT      = StdHEP->T();
            }
        }

      if(tob->GetName()==TString("DigiScintHits"))
        {
          TClonesArray* ScintHitList = (TClonesArray*)(tob); 
          MSG("NtpMaker",Msg::kDebug) << " ... found DigitScintHits " << endl;
          TClonesArray& MyScintHitList = *(fEvent->ScintHitList);

          Int_t plane,strip;
          Int_t idnu,ipdg,trkid;
          Double_t u,v,x,y,z,r;
          Double_t xm,xp,ym,yp,um,up,vm,vp;
          Double_t muvtxU,muvtxV;
          Double_t muvtxX,muvtxY,muvtxZ;
          Double_t muvtxR,muendR;
          Double_t muendU,muendV;
          Double_t muendX,muendY,muendZ;
          Int_t muvtxpln,muendpln;
          Int_t muvtxstrp,muendstrp;
          Double_t pmu,maxpmu,minpmu;

          muvtxU=0.0; muvtxV=0.0; 
          muvtxX=0.0; muvtxY=0.0; muvtxZ=0.0; 
          muendU=0.0; muendV=0.0; 
          muendX=0.0; muendY=0.0; muendZ=0.0; 
          muvtxR=999.9; muendR=999.9;
          muvtxpln=-999; muendpln=-999;
          muvtxstrp=-999; muendstrp=-999; 
          maxpmu=-999.9; minpmu=-999.9;

          PlexStripEndId seid;
          UgliStripHandle ush;
          TVector3 myglobal;
          TIter ScintHitItr(ScintHitList->MakeIterator());
          while( const DigiScintHit* ScintHit = dynamic_cast<const DigiScintHit*>(ScintHitItr()))
            {
              if( ScintHit->DE()<1.0E-4 && abs(ScintHit->ParticleId())!=13 ) continue; 
                 //this is to prevent huge numbers of scinthits being stored in the ntuple

              AtmosScintHit* MyScintHit = new( MyScintHitList[(fEvent->NScintHits)++] ) AtmosScintHit();
              MyScintHit->Id = ScintHit->ParticleId();
              MyScintHit->TrkId = ScintHit->TrackId();        
              MyScintHit->Plane = ScintHit->Plane();
              MyScintHit->Strip = ScintHit->Strip();
                
              seid = ScintHit->StripEndId();
              MyScintHit->View = -1;
              if(seid.GetPlaneView()==PlaneView::kU) MyScintHit->View = 0;
              if(seid.GetPlaneView()==PlaneView::kV) MyScintHit->View = 1;

              ush = ugh.GetStripHandle(seid);
      
              TVector3 mylocal( ScintHit->X1(), ScintHit->Y1(), ScintHit->Z1() );
              myglobal = ush.LocalToGlobal(mylocal);
            //Note: doubles are all being converted to floats here
            // SHOULD PROBABLY DO CHECKS HERE TO PREVENT DOUBLES LARGER 
            // THAN FLT_MAX FROM SCREWING THINGS UP
              MyScintHit->U[0]  = float(0.7071*(myglobal.Y()+myglobal.X())); 
              MyScintHit->V[0]  = float(0.7071*(myglobal.Y()-myglobal.X()));
              MyScintHit->X[0]  = float(myglobal.X()); 
              MyScintHit->Y[0]  = float(myglobal.Y());
              MyScintHit->Z[0]  = float(myglobal.Z());
              MyScintHit->T[0]  = float(ScintHit->T1());

              mylocal.SetXYZ( ScintHit->X2(), ScintHit->Y2(), ScintHit->Z2() );
              myglobal = ush.LocalToGlobal(mylocal);
       
              MyScintHit->U[1]  = float(0.7071*(myglobal.Y()+myglobal.X())); 
              MyScintHit->V[1]  = float(0.7071*(myglobal.Y()-myglobal.X()));
              MyScintHit->X[1]  = float(myglobal.X()); 
              MyScintHit->Y[1]  = float(myglobal.Y());
              MyScintHit->Z[1]  = float(myglobal.Z());
              MyScintHit->T[1]  = float(ScintHit->T2());

              MyScintHit->DS = float(ScintHit->DS());
              MyScintHit->DE = float(ScintHit->DE());
              MyScintHit->ParticleE = float(ScintHit->ParticleEnergy());

              idnu = fEvent->MCInfo.IDnu;
              trkid = ScintHit->TrackId();
              ipdg = ScintHit->ParticleId();

              if( trkid>=0 
               && ( (idnu==0 && (ipdg==13 || ipdg==-13) ) 
                 || (idnu==14 && ipdg==13) || (idnu==-14 && ipdg==-13) ) )
                {
                  plane = ScintHit->Plane();
                  strip = ScintHit->Strip();
                  x = myglobal.X(); y = myglobal.Y(); z = myglobal.Z();
                  u = 0.7071*(y + x); v = 0.7071*(y - x);  
                  pmu = ScintHit->ParticleEnergy();

                  r=999.9;
                  up=4.0-u; if(up<r) r=up;
                  um=4.0+u; if(um<r) r=um;
                  vp=4.0-v; if(vp<r) r=vp;
                  vm=4.0+v; if(vm<r) r=vm;
                  xp=4.0-x; if(xp<r) r=xp;
                  xm=4.0+x; if(xm<r) r=xm;
                  yp=4.0-y; if(yp<r) r=yp;
                  ym=4.0+y; if(ym<r) r=ym;

                  if( minpmu<0 || pmu<minpmu )
                    {
                      muendU=u; muendV=v;
                      muendX=x; muendY=y; muendZ=z; muendR=r;
                      muendstrp=strip; muendpln=plane; minpmu=pmu;
                    }
                  if( maxpmu<0 || pmu>maxpmu )
                    {
                      muvtxU=u; muvtxV=v;
                      muvtxX=x; muvtxY=y; muvtxZ=z; muvtxR=r;
                      muvtxstrp=strip; muvtxpln=plane; maxpmu=pmu;
                    }
                }
            }

          fEvent->MCInfo.MuVtxU = muvtxU;              
          fEvent->MCInfo.MuVtxV = muvtxV;    
          fEvent->MCInfo.MuVtxX = muvtxX;              
          fEvent->MCInfo.MuVtxY = muvtxY;               
          fEvent->MCInfo.MuVtxZ = muvtxZ;               
          fEvent->MCInfo.MuVtxDistToEdge = muvtxR;               
          fEvent->MCInfo.MuVtxPlane = muvtxpln; 
          fEvent->MCInfo.MuVtxStrip = muvtxstrp; 
          fEvent->MCInfo.MuEndU = muendU;            
          fEvent->MCInfo.MuEndV = muendV;      
          fEvent->MCInfo.MuEndX = muendX;            
          fEvent->MCInfo.MuEndY = muendY;            
          fEvent->MCInfo.MuEndZ = muendZ;            
          fEvent->MCInfo.MuEndDistToEdge = muendR;            
          fEvent->MCInfo.MuEndPlane = muendpln;          
          fEvent->MCInfo.MuEndStrip = muendstrp;      
          fEvent->MCInfo.MuPbeg = maxpmu;               
          fEvent->MCInfo.MuPend = minpmu;               

        }
    }

  return;
}

//======================================================================

void NtpMaker::FillRawInfo(RawRecord* rawrec)
{
  MSG("NtpMaker",Msg::kDebug) << " *** NtpMaker::FillRawInfo(...) *** " << endl;
  
  const VldContext& vldc = *(rawrec->GetVldContext());
  this->FillSpillInfo(vldc);

  return;
}

//=======================================================================

void NtpMaker::FillShieldInfo(TObject* tob)
{
  FarDetShieldPlankListHandle* ShldList = (FarDetShieldPlankListHandle*)(tob);
  MSG("NtpMaker",Msg::kDebug) << " ... found FarDetShieldPlankListHandle " << endl;
  TClonesArray& MyShldList = *(fEvent->ShieldPlankList);

  TIter ShldItr(ShldList->GetDaughterIterator());
  while(FarDetShieldPlankHandle* Shld = (FarDetShieldPlankHandle*)(ShldItr()))
    {
      AtmosShieldPlank* MyShld = new( MyShldList[(fEvent->NShieldPlanks)++] ) AtmosShieldPlank();

      MyShld->Section = Shld->GetSection();       
      MyShld->SubSection = Shld->GetSubSection();
      MyShld->Plane = Shld->GetPlane();
      MyShld->Plank = Shld->GetPlank();
      MyShld->Nstrips = Shld->GetNStrips();
      MyShld->Nerrors = Shld->GetGeomErrors();
      MyShld->Ndigits = Shld->GetNDaughters();    
      MyShld->X = Shld->GetX();
      MyShld->Y = Shld->GetY();  
      MyShld->Z[0] = Shld->GetZ(StripEnd::kNegative);
      MyShld->Z[1] = Shld->GetZ(StripEnd::kPositive);
      MyShld->Traw[0] = Shld->GetTime(StripEnd::kNegative,CalTimeType::kNone);
      MyShld->Traw[1] = Shld->GetTime(StripEnd::kPositive,CalTimeType::kNone);
      MyShld->Tcal[0] = Shld->GetTime(StripEnd::kNegative,CalTimeType::kT0);
      MyShld->Tcal[1] = Shld->GetTime(StripEnd::kPositive,CalTimeType::kT0);
      MyShld->Qadc[0] = Shld->GetCharge(StripEnd::kNegative,CalDigitType::kNone);
      MyShld->Qadc[1] = Shld->GetCharge(StripEnd::kPositive,CalDigitType::kNone);
      MyShld->QPE[0] = Shld->GetCharge(StripEnd::kNegative,CalDigitType::kPE);
      MyShld->QPE[1] = Shld->GetCharge(StripEnd::kPositive,CalDigitType::kPE);
      //Do we add info for sigmap and mips?
      MyShld->GreenFibre[0] = Shld->GetGreenFibre(StripEnd::kNegative);
      MyShld->GreenFibre[1] = Shld->GetGreenFibre(StripEnd::kPositive);
      MyShld->WlsPigtail[0] = Shld->GetWlsPigtail(StripEnd::kNegative);
      MyShld->WlsPigtail[1] = Shld->GetWlsPigtail(StripEnd::kPositive);
      MyShld->ClearFibre[0] = Shld->GetClearFibre(StripEnd::kNegative);
      MyShld->ClearFibre[1] = Shld->GetClearFibre(StripEnd::kPositive);

    }
  return;
}

//======================================================================

void NtpMaker::FillShowerInfo(TObject* tob)
{
  CandShowerListHandle* ShwList = (CandShowerListHandle*)(tob);
  MSG("NtpMaker",Msg::kDebug) << " ... found CandShowerListHandle containing "<<ShwList->GetNDaughters()<<" showers."<< endl;
  float nstripcalib(0.0), nshwstrips(0.0);
  Calibrator& cal = Calibrator::Instance();
  TClonesArray* MyStrpList = (TClonesArray*)(fEvent->StripList);

  // Set information from CandShowerHandle
  TClonesArray& MyShwList = *(fEvent->ShowerList);
  TIter ShwItr(ShwList->GetDaughterIterator());
  while(CandShowerHandle* Shw = (CandShowerHandle*)(ShwItr()))
    {
      AtmosShower* MyShw = new( MyShwList[(fEvent->NShowers)++] ) AtmosShower();
        
      MyShw->VtxPlane = Shw->GetVtxPlane();
      MyShw->VtxTime = Shw->GetVtxT();
      MyShw->VtxU = Shw->GetVtxU();
      MyShw->VtxV = Shw->GetVtxV();
      MyShw->VtxR = sqrt((MyShw->VtxU*MyShw->VtxU)+(MyShw->VtxV*MyShw->VtxV));
      MyShw->VtxX = 0.7071*(Shw->GetVtxU()-Shw->GetVtxV());
      MyShw->VtxY = 0.7071*(Shw->GetVtxU()+Shw->GetVtxV());
      MyShw->VtxZ = Shw->GetVtxZ();

      MyShw->VtxDirCosU = Shw->GetDirCosU();
      MyShw->VtxDirCosV = Shw->GetDirCosV();
      MyShw->VtxDirCosX = 0.7071*(Shw->GetDirCosU()-Shw->GetDirCosV());
      MyShw->VtxDirCosY = 0.7071*(Shw->GetDirCosU()+Shw->GetDirCosV());
      MyShw->VtxDirCosZ = Shw->GetDirCosZ();   
      MyShw->TimeSlope = Shw->GetTimeSlope();
      MyShw->TimeOffset = Shw->GetTimeOffset(); 
      MyShw->Energy = Shw->GetEnergy();
      
      MyShw->Nstrips = Shw->GetNStrip(); //AtNu Shw->GetNStrips();
      MyShw->Nplanes = Shw->GetNPlane(); //AtNu Shw->GetNPlanes();
      MyShw->Index = (int)pow(2.,fEvent->NShowers-1);

      // Set information for strips that are part of this shower      
      nshwstrips=Shw->GetNDaughters(); nstripcalib=0.0;
      TIter StrpItr(Shw->GetDaughterIterator());
      PlexStripEndId* seidP = 0;
      PlexStripEndId* seidN = 0;
      while(CandStripHandle* ShwStrp = (CandStripHandle*)(StrpItr()))
        {
          for(Int_t myint=0;myint<fEvent->NStrips;myint++)
            {
              AtmosStrip& MyStrp = *((AtmosStrip*)(MyStrpList->At(myint)));
              if( ShwStrp->GetPlane()==MyStrp.Plane && ShwStrp->GetStrip()==MyStrp.Strip )
                {
                  if( MyStrp.Trk==0 && MyStrp.Shw==0 && Shw->GetU(MyStrp.Plane)>-999.9 && Shw->GetV(MyStrp.Plane)>-999.9 )
                    {
                      
                      if( MyStrp.View==0 )
                        {
                          MyStrp.L = Shw->GetV(MyStrp.Plane);
                          MyStrp.GreenFibre[0] = MyStrp.HalfLength - MyStrp.L;
                          MyStrp.GreenFibre[1] = MyStrp.HalfLength + MyStrp.L;
                          seidN = new PlexStripEndId(Detector::kFar, MyStrp.Plane, MyStrp.Strip, StripEnd::kNegative, StripEnd::kWhole, PlaneView::kU, PlaneCoverage::kTotal);
                          seidP = new PlexStripEndId(Detector::kFar, MyStrp.Plane, MyStrp.Strip, StripEnd::kPositive, StripEnd::kWhole, PlaneView::kU, PlaneCoverage::kTotal);
                        }
                      if( MyStrp.View==1 )
                        {
                          MyStrp.L = Shw->GetU(MyStrp.Plane);
                          MyStrp.GreenFibre[0] = MyStrp.HalfLength + MyStrp.L;
                          MyStrp.GreenFibre[1] = MyStrp.HalfLength - MyStrp.L;
                          seidN = new PlexStripEndId(Detector::kFar, MyStrp.Plane, MyStrp.Strip, StripEnd::kNegative, StripEnd::kWhole, PlaneView::kV, PlaneCoverage::kTotal);
                          seidP = new PlexStripEndId(Detector::kFar, MyStrp.Plane, MyStrp.Strip, StripEnd::kPositive, StripEnd::kWhole, PlaneView::kV, PlaneCoverage::kTotal);
                        }
                      double lsign = pow(-1.0, (1.0+double(MyStrp.View)) );
                      if(seidN && MyStrp.Sigcorr[0]>0.0)
                        {
                          MyStrp.Sigmap[0] = cal.GetAttenCorrected(MyStrp.Sigcorr[0], lsign*MyStrp.L, *(seidN) );
                          MyStrp.MIP[0] = cal.GetMIP(MyStrp.Sigmap[0], *(seidN) );
                        }
                      if(seidP && MyStrp.Sigcorr[1]>0.0)
                        {
                          MyStrp.Sigmap[1] = cal.GetAttenCorrected(MyStrp.Sigcorr[1], lsign*MyStrp.L, *(seidP) );
                          MyStrp.MIP[1] = cal.GetMIP(MyStrp.Sigmap[1], *(seidP) );
                        }
                    }
                  if(MyStrp.Sigcorr[0]>0.0 || MyStrp.Sigcorr[1]>0.0) nstripcalib+=1.0;
                  MyStrp.Shw += MyShw->Index;
                }

              if(seidP) delete seidP; seidP = 0;
              if(seidN) delete seidN; seidN = 0;
            }
        }
      MSG("NtpMaker",Msg::kVerbose) <<"   "<<(int)nshwstrips<<" strips in the shower of which "<<(int)nstripcalib<<" were calibrated."<<endl<<flush;
    }
  return;
}

//======================================================================

void NtpMaker::FillSliceInfo(TObject* tob)
{
  FarDetSliceListHandle* SliceList = (FarDetSliceListHandle*)(tob);
  MSG("NtpMaker",Msg::kDebug) << " ... found FarDetSliceListHandle " << endl;

  TIter SliceItr = SliceList->GetDaughterIterator();
  FarDetSliceHandle* Slice = dynamic_cast<FarDetSliceHandle*>(SliceItr()); 
  if( Slice )
    {
      fEvent->FilterInfo.EventId = Slice->GetEventId();
      fEvent->FilterInfo.EventIdSM1 = Slice->GetEventIdSM1();
      fEvent->FilterInfo.EventIdSM2 = Slice->GetEventIdSM2();
      fEvent->FilterInfo.EventIdX = Slice->GetEventIdX();
      fEvent->FilterInfo.Edges=Slice->GetEdges();
      fEvent->FilterInfo.EdgesSM1=Slice->GetEdgesSM1();
      fEvent->FilterInfo.EdgesSM2=Slice->GetEdgesSM2();
      fEvent->FilterInfo.GoodUPlanes = Slice->GetGoodUPlanes();
      fEvent->FilterInfo.GoodVPlanes = Slice->GetGoodVPlanes();
      fEvent->FilterInfo.GoodPlanes = Slice->GetGoodPlanes();
      fEvent->FilterInfo.GoodStrips = Slice->GetGoodStrips();
      fEvent->FilterInfo.FidCharge = Slice->GetFidCharge();
      fEvent->FilterInfo.TotalCharge = Slice->GetTotalCharge();
      fEvent->FilterInfo.MaxPlaneCharge = Slice->GetMaxPlaneCharge();

      for(Int_t myint=0; myint<11; myint++)
        {
          fEvent->FilterInfo.FidChargeSM1[myint] = Slice->GetFidChargeSM1(myint);
          fEvent->FilterInfo.FidChargeSM2[myint] = Slice->GetFidChargeSM2(myint);
          fEvent->FilterInfo.MeanXPosSM1[myint] = Slice->GetMeanXPosSM1(myint);
          fEvent->FilterInfo.MeanXPosSM2[myint] = Slice->GetMeanXPosSM2(myint);
          fEvent->FilterInfo.MeanYPosSM1[myint] = Slice->GetMeanYPosSM1(myint);
          fEvent->FilterInfo.MeanYPosSM2[myint] = Slice->GetMeanYPosSM2(myint);
          fEvent->FilterInfo.MeanZPosSM1[myint] = Slice->GetMeanZPosSM1(myint);
          fEvent->FilterInfo.MeanZPosSM2[myint] = Slice->GetMeanZPosSM2(myint);
        }
    } 
  return;
}

//======================================================================

void NtpMaker::FillSpillInfo(const VldContext& vldc)
{
  MSG("NtpMaker",Msg::kDebug) << " *** NtpMaker::FillSpillInfo(...) *** " << endl;
  double dt1(0.0),dt2(0.0),dt(0.0); 

  // only apply to data
  if(vldc.GetSimFlag()==SimFlag::kData)
    {
      
      SpillTimeFinder &stf = SpillTimeFinder::Instance();
      fEvent->SpillInfo.stfDataAvailable = stf.DataIsAvailable(vldc);
      fEvent->SpillInfo.TimeToNearestSpill_stf = stf.GetTimeToNearestSpill(vldc);

      const SpillTimeND& ndspill = stf.GetNearestSpill(vldc);
     
      const BeamMonSpill* bmspill = BDSpillAccessor::Get().LoadSpill(vldc.GetTimeStamp());

      // fill members with info from the BeamMonSpill
      if( bmspill )
        {    
          VldTimeStamp ndspilltime = ndspill.GetTimeStamp();
          VldTimeStamp bmspilltime = bmspill->SpillTime();
          dt1 = bmspilltime.GetSec() - ndspilltime.GetSec();
          dt2 = ( bmspilltime.GetNanoSec() - ndspilltime.GetNanoSec() )*1e-9;
          dt = dt1 + dt2;

          fEvent->SpillInfo.bmsDataAvailable = false;
          fEvent->SpillInfo.TimeToNearestSpill_bms = fEvent->SpillInfo.TimeToNearestSpill_stf + dt;

          if( abs(fEvent->SpillInfo.TimeToNearestSpill_stf)<86400.0 )
            {
              fEvent->SpillInfo.bmsDataAvailable = true;

              // Beam Status Info
              fEvent->SpillInfo.beam_type = bmspill->BeamType();
              fEvent->SpillInfo.bmsStatus = bmspill->GetStatusInt();

              BeamMonSpill::StatusBits status_bits = bmspill->GetStatusBits();
              fEvent->SpillInfo.horn_on = status_bits.horn_on;
              fEvent->SpillInfo.target_in = status_bits.target_in;
              fEvent->SpillInfo.n_batches = status_bits.n_batches;

              //Beam Monitoring Info
              fEvent->SpillInfo.tor101 = bmspill->fTor101;           
              fEvent->SpillInfo.tr101d = bmspill->fTr101d;           
              fEvent->SpillInfo.tortgt = bmspill->fTortgt;           
              fEvent->SpillInfo.trtgtd = bmspill->fTrtgtd;
              fEvent->SpillInfo.horncur = bmspill->fHornCur;     
              for (Int_t i=0;i<6;++i)
                {
                  fEvent->SpillInfo.bposx[i]=bmspill->fTargBpmX[i];
                  fEvent->SpillInfo.bposy[i]=bmspill->fTargBpmY[i];
                  fEvent->SpillInfo.bpmint[i]=bmspill->fBpmInt[i];
                }
              fEvent->SpillInfo.bwidx = bmspill->fProfWidX;  
              fEvent->SpillInfo.bwidy = bmspill->fProfWidY;  
              fEvent->SpillInfo.hadint = bmspill->fHadInt; 
              fEvent->SpillInfo.muint1 = bmspill->fMuInt1;  
              fEvent->SpillInfo.muint2 = bmspill->fMuInt2;  
              fEvent->SpillInfo.muint3 = bmspill->fMuInt3; 

              VldTimeStamp BMsubtime = bmspilltime-ndspill.GetTimeStamp();
              BMSpillAna BMana;
              BMana.UseDatabaseCuts();
              BMana.SetSpill(*bmspill);
              BMana.SetTimeDiff(BMsubtime.GetSeconds());
              fEvent->SpillInfo.BeamSelectSpill = BMana.SelectSpill(); 
            }
        }
      else
        { 
          MSG("NtpMaker",Msg::kDebug)<< "No BeamMonSpill found!" <<endl;
        }
  
    }

  return;
}

//======================================================================

void NtpMaker::FillStripInfo(TObject* tob)
{
  FarDetStripListHandle* StrpList = (FarDetStripListHandle*)(tob);
  MSG("NtpMaker",Msg::kDebug) << " ... found FarDetStripListHandle " << endl;
  TClonesArray& MyStrpList = *(fEvent->StripList);
  Calibrator& cal = Calibrator::Instance();
  float nstripcalib(0.0);
  Int_t pln;
  Double_t totq,totqo,opos;

  TIter StrpItr(StrpList->GetDaughterIterator());        
  while(FarDetStripHandle* Strp = dynamic_cast<FarDetStripHandle*>(StrpItr()))
    {
      AtmosStrip* MyStrp = new( MyStrpList[ (fEvent->NStrips)++] ) AtmosStrip();

      MyStrp->Plane = Strp->GetPlane();
      MyStrp->Strip = Strp->GetStrip();
      if(Strp->GetPlaneView()==PlaneView::kU) MyStrp->View = 0;
      if(Strp->GetPlaneView()==PlaneView::kV) MyStrp->View = 1;
      MyStrp->L = 0.0;
      MyStrp->T = Strp->GetTPos();
      MyStrp->Z = Strp->GetZPos();
      MyStrp->Ndigits = Strp->GetNDaughters();
      MyStrp->Traw[0] = Strp->GetTime(StripEnd::kNegative,CalTimeType::kNone);
      MyStrp->Traw[1] = Strp->GetTime(StripEnd::kPositive,CalTimeType::kNone);
      MyStrp->Tcal[0] = Strp->GetTime(StripEnd::kNegative,CalTimeType::kT0);
      MyStrp->Tcal[1] = Strp->GetTime(StripEnd::kPositive,CalTimeType::kT0);
      MyStrp->Qadc[0] = Strp->GetCharge(StripEnd::kNegative,CalDigitType::kNone);
      MyStrp->Qadc[1] = Strp->GetCharge(StripEnd::kPositive,CalDigitType::kNone);
      MyStrp->QPE[0] = Strp->GetCharge(StripEnd::kNegative,CalDigitType::kPE);
      MyStrp->QPE[1] = Strp->GetCharge(StripEnd::kPositive,CalDigitType::kPE);
      MyStrp->Sigcorr[0] = Strp->GetCharge(StripEnd::kNegative,CalDigitType::kSigCorr);
      MyStrp->Sigcorr[1] = Strp->GetCharge(StripEnd::kPositive,CalDigitType::kSigCorr);
      MyStrp->QPEcorr[0] = MyStrp->Sigcorr[0]/65.0; // PEcorr=Sigcorr/<Gain>
      MyStrp->QPEcorr[1] = MyStrp->Sigcorr[1]/65.0; // hack the conversion
      MyStrp->GreenFibre[0] = Strp->GetGreenFibre(StripEnd::kNegative);
      MyStrp->GreenFibre[1] = Strp->GetGreenFibre(StripEnd::kPositive);     
      MyStrp->WlsPigtail[0] = Strp->GetWlsPigtail(StripEnd::kNegative);
      MyStrp->WlsPigtail[1] = Strp->GetWlsPigtail(StripEnd::kPositive);
      MyStrp->ClearFibre[0] = Strp->GetClearFibre(StripEnd::kNegative);
      MyStrp->ClearFibre[1] = Strp->GetClearFibre(StripEnd::kPositive);
      MyStrp->HalfLength = Strp->GetHalfLength();
      MyStrp->Xtalk = Strp->IsXtalk(StripEnd::kWhole);
      
      if( MyStrp->Plane>0 && MyStrp->Plane<500 && MyStrp->Xtalk==0 && MyStrp->QPE[0]+MyStrp->QPE[1]>2.0 )
        {
          fStrpList[MyStrp->Plane].push_back(*MyStrp);
        }
    }
  vector<AtmosStrip>::iterator PlaneStrpIter;
  vector<AtmosStrip>::iterator PlaneListEnd;
  PlexStripEndId* seidP = 0;
  PlexStripEndId* seidN = 0;
  for(Int_t myint=0;myint<fEvent->NStrips;myint++)
    {
      AtmosStrip& MyStrp = *((AtmosStrip*)(MyStrpList.At(myint)));
      pln = MyStrp.Plane;

      totqo=0.0; totq=0.0; opos=0.0;
      if( (pln>1 && pln<249) || (pln>250 && pln<498) )
        {
          PlaneStrpIter = fStrpList[pln-1].begin();
          PlaneListEnd = fStrpList[pln-1].end();
          while(PlaneStrpIter!=PlaneListEnd)
            {
              AtmosStrip& mystrip = *(PlaneStrpIter);  
              totqo+=(mystrip.QPE[0]+mystrip.QPE[1])*(mystrip.T);
              totq+=(mystrip.QPE[0]+mystrip.QPE[1]);
              ++PlaneStrpIter;
            } 
        }        
      if( (pln>0 && pln<248) || (pln>249 && pln<497) )
        {
          PlaneStrpIter = fStrpList[pln+1].begin();
          PlaneListEnd = fStrpList[pln+1].end();
          while(PlaneStrpIter!=PlaneListEnd)
            {
              AtmosStrip& mystrip = *(PlaneStrpIter);  
              totqo+=(mystrip.QPE[0]+mystrip.QPE[1])*(mystrip.T);
              totq+=(mystrip.QPE[0]+mystrip.QPE[1]);
              ++PlaneStrpIter;
            }
        }
      if(totq>0.0){ opos=totqo/totq; } else{ opos=0.0; }

      MyStrp.L = opos;
      if( MyStrp.View==0 )
        {           
          MyStrp.GreenFibre[0] = MyStrp.HalfLength - MyStrp.L;
          MyStrp.GreenFibre[1] = MyStrp.HalfLength + MyStrp.L;
          seidN = new PlexStripEndId(Detector::kFar, MyStrp.Plane, MyStrp.Strip, StripEnd::kNegative, StripEnd::kWhole, PlaneView::kU, PlaneCoverage::kTotal);
          seidP = new PlexStripEndId(Detector::kFar, MyStrp.Plane, MyStrp.Strip, StripEnd::kPositive, StripEnd::kWhole, PlaneView::kU, PlaneCoverage::kTotal);
        }
      if( MyStrp.View==1 )
        {
          seidN = new PlexStripEndId(Detector::kFar, MyStrp.Plane, MyStrp.Strip, StripEnd::kNegative, StripEnd::kWhole, PlaneView::kV, PlaneCoverage::kTotal);
          seidP = new PlexStripEndId(Detector::kFar, MyStrp.Plane, MyStrp.Strip, StripEnd::kPositive, StripEnd::kWhole, PlaneView::kV, PlaneCoverage::kTotal);
          MyStrp.GreenFibre[0] = MyStrp.HalfLength + MyStrp.L;
          MyStrp.GreenFibre[1] = MyStrp.HalfLength - MyStrp.L;
        }
      
      double lsign = pow(-1.0, (1.0+double(MyStrp.View)) );
      if(seidN && MyStrp.Sigcorr[0]>0.0)
        {
          MyStrp.Sigmap[0] = cal.GetAttenCorrected(MyStrp.Sigcorr[0], lsign*opos, *(seidN) );
          MyStrp.MIP[0] = cal.GetMIP(MyStrp.Sigmap[0], *(seidN) );
        }
      if(seidP && MyStrp.Sigcorr[1]>0.0)
        {
          MyStrp.Sigmap[1] = cal.GetAttenCorrected(MyStrp.Sigcorr[1], lsign*opos, *(seidP) );
          MyStrp.MIP[1] = cal.GetMIP(MyStrp.Sigmap[1], *(seidP) );
        }
      if(MyStrp.Sigcorr[0]>0.0 || MyStrp.Sigcorr[1]>0.0) nstripcalib+=1.0;

      if(seidP) delete seidP; seidP = 0;
      if(seidN) delete seidN; seidN = 0;
    }
  MSG("NtpMaker",Msg::kVerbose) <<"   "<<(int)(fEvent->NStrips)<<" strips in the event of which "<<(int)nstripcalib<< " were calibrated."<<endl<<flush;

  for(Int_t myi=0;myi<500;myi++)
    {
      fStrpList[myi].clear();        
    }
  return;
}

//======================================================================

void NtpMaker::FillTrackInfo(TObject* tob)
{
  CandTrackListHandle* TrkList = (CandTrackListHandle*)(tob);
  MSG("NtpMaker",Msg::kDebug) << " ... found CandTrackListHandle containing "<<TrkList->GetNDaughters()<<" tracks."<< endl;
  float nstripcalib(0.0),stripsintrack(0.0);
  Calibrator& cal = Calibrator::Instance();
  TClonesArray* MyStrpList = (TClonesArray*)(fEvent->StripList);

  // Set information from CandTrackHandle
  TClonesArray& MyTrkList = *(fEvent->TrackList);
  TIter TrkItr = TrkList->GetDaughterIterator();
  while(CandTrackHandle* Trk = dynamic_cast<CandTrackHandle*>(TrkItr()))
    {
      AtmosTrack* MyTrk = new( MyTrkList[(fEvent->NTracks)++] ) AtmosTrack();
      MyTrk->VtxPlane = Trk->GetVtxPlane();
      MyTrk->VtxTime = Trk->GetVtxT();
      MyTrk->VtxU = Trk->GetVtxU();
      MyTrk->VtxV = Trk->GetVtxV();
      MyTrk->VtxR = sqrt( (MyTrk->VtxU*MyTrk->VtxU) + (MyTrk->VtxV*MyTrk->VtxV) );
      MyTrk->VtxX = 0.7071*(Trk->GetVtxU()-Trk->GetVtxV());
      MyTrk->VtxY = 0.7071*(Trk->GetVtxU()+Trk->GetVtxV());
      MyTrk->VtxZ = Trk->GetVtxZ();
      MyTrk->VtxTrace = Trk->GetVtxTrace();
      MyTrk->VtxTraceZ = Trk->GetVtxTraceZ();
      MyTrk->VtxDirCosU = Trk->GetDirCosU();
      MyTrk->VtxDirCosV = Trk->GetDirCosV();
      MyTrk->VtxDirCosX = 0.7071*(Trk->GetDirCosU()-Trk->GetDirCosV());
      MyTrk->VtxDirCosY = 0.7071*(Trk->GetDirCosU()+Trk->GetDirCosV());
      MyTrk->VtxDirCosZ = Trk->GetDirCosZ();
      MyTrk->EndPlane = Trk->GetEndPlane();
      MyTrk->EndTime = Trk->GetEndT();
      MyTrk->EndU = Trk->GetEndU();
      MyTrk->EndV = Trk->GetEndV();
      MyTrk->EndR = sqrt((MyTrk->EndU*MyTrk->EndU)+(MyTrk->EndV*MyTrk->EndV));
      MyTrk->EndX = 0.7071*(Trk->GetEndU()-Trk->GetEndV());
      MyTrk->EndY = 0.7071*(Trk->GetEndU()+Trk->GetEndV());
      MyTrk->EndZ = Trk->GetEndZ();
      MyTrk->EndTrace = Trk->GetEndTrace();
      MyTrk->EndTraceZ = Trk->GetEndTraceZ();
      MyTrk->EndDirCosU = Trk->GetEndDirCosU();
      MyTrk->EndDirCosV = Trk->GetEndDirCosV();
      MyTrk->EndDirCosX = 0.7071*(Trk->GetEndDirCosU()-Trk->GetEndDirCosV());
      MyTrk->EndDirCosY = 0.7071*(Trk->GetEndDirCosU()+Trk->GetEndDirCosV());
      MyTrk->EndDirCosZ = Trk->GetEndDirCosZ();
      MyTrk->RangeGCM2 = Trk->GetRange(Trk->GetEndPlane());
       // NB this will lose the extrapolation at vtx and end
      MyTrk->RangeMetres = Trk->GetdS(Trk->GetEndPlane());
       // NB this will lose the extrapolation at vtx and end
      MyTrk->Nstrips = Trk->GetNStrip(); 
      MyTrk->Nplanes = abs(MyTrk->EndPlane-MyTrk->VtxPlane)+1;
      MyTrk->AtNuNplanes = MyTrk->Nplanes;

      MyTrk->TimingFitChi2 = Trk->GetTimeFitChi2();
      MyTrk->TimingFitNdf = Trk->GetNTimeFitDigit()-2;
      MyTrk->TimeSlope = Trk->GetTimeSlope();
      MyTrk->TimeOffset = Trk->GetTimeOffset(); 

      if(MyTrk->EndPlane-MyTrk->VtxPlane>0)
        {//Forward and Backward variables defined relative to increasing Z
          MyTrk->VtxDirTimeFitRMS = Trk->GetTimeForwardFitRMS();  
          MyTrk->EndDirTimeFitRMS = Trk->GetTimeBackwardFitRMS(); 
          MyTrk->VtxDirTimeFitNdf = Trk->GetTimeForwardFitNDOF(); 
          MyTrk->EndDirTimeFitNdf = Trk->GetTimeBackwardFitNDOF();
        }
      else
        {
          MyTrk->VtxDirTimeFitRMS = Trk->GetTimeBackwardFitRMS(); 
          MyTrk->EndDirTimeFitRMS = Trk->GetTimeForwardFitRMS();  
          MyTrk->VtxDirTimeFitNdf = Trk->GetTimeBackwardFitNDOF();
          MyTrk->EndDirTimeFitNdf = Trk->GetTimeForwardFitNDOF(); 
        }

      MyTrk->Momentum = Trk->GetMomentum();
      MyTrk->MomentumRange = Trk->GetMomentum();

      MyTrk->Index = (int)pow(2.,fEvent->NTracks-1);

      //set information for strips that are part of this track
      TIter StrpItr(Trk->GetDaughterIterator());
      nstripcalib=(0.0); stripsintrack = Trk->GetNDaughters();
      PlexStripEndId* seidP = 0;
      PlexStripEndId* seidN = 0;
      while(CandStripHandle* TrkStrp = (CandStripHandle*)(StrpItr()))
        {
          for(Int_t myint=0;myint<fEvent->NStrips;myint++)
            {
              AtmosStrip& MyStrp = *((AtmosStrip*)(MyStrpList->At(myint)));
              
              if( TrkStrp->GetPlane()==MyStrp.Plane && TrkStrp->GetStrip()==MyStrp.Strip )
                {
                  if( MyStrp.Trk==0 && Trk->IsTPosValid(MyStrp.Plane) )
                    {
                      if( MyStrp.View==0 )
                        {
                          MyStrp.L = Trk->GetV(MyStrp.Plane);
                          MyStrp.dS = Trk->GetdS(MyStrp.Plane);
                          MyStrp.GreenFibre[0] = MyStrp.HalfLength - MyStrp.L;
                          MyStrp.GreenFibre[1] = MyStrp.HalfLength + MyStrp.L;
                          seidN = new PlexStripEndId(Detector::kFar, MyStrp.Plane, MyStrp.Strip, StripEnd::kNegative, StripEnd::kWhole, PlaneView::kU, PlaneCoverage::kTotal);
                          seidP = new PlexStripEndId(Detector::kFar, MyStrp.Plane, MyStrp.Strip, StripEnd::kPositive, StripEnd::kWhole, PlaneView::kU, PlaneCoverage::kTotal);
                        }
                      if( MyStrp.View==1 )
                        {
                          MyStrp.L = Trk->GetU(MyStrp.Plane);
                          MyStrp.dS = Trk->GetdS(MyStrp.Plane);
                          MyStrp.GreenFibre[0] = MyStrp.HalfLength + MyStrp.L;
                          MyStrp.GreenFibre[1] = MyStrp.HalfLength - MyStrp.L;
                          seidN = new PlexStripEndId(Detector::kFar, MyStrp.Plane, MyStrp.Strip, StripEnd::kNegative, StripEnd::kWhole, PlaneView::kV, PlaneCoverage::kTotal);
                          seidP = new PlexStripEndId(Detector::kFar, MyStrp.Plane, MyStrp.Strip, StripEnd::kPositive, StripEnd::kWhole, PlaneView::kV, PlaneCoverage::kTotal);
                        }
                      double lsign = pow(-1.0, (1.0+double(MyStrp.View)) );
                      if(seidN && MyStrp.Sigcorr[0]>0.0)
                        {
                          MyStrp.Sigmap[0] = cal.GetAttenCorrected(MyStrp.Sigcorr[0], lsign*MyStrp.L, *(seidN) );
                          MyStrp.MIP[0] = cal.GetMIP(MyStrp.Sigmap[0], *(seidN) );
                        }
                      if(seidP && MyStrp.Sigcorr[1]>0.0)
                        {
                          MyStrp.Sigmap[1] = cal.GetAttenCorrected(MyStrp.Sigcorr[1], lsign*MyStrp.L, *(seidP) );
                          MyStrp.MIP[1] = cal.GetMIP(MyStrp.Sigmap[1], *(seidP) );
                        }
                      if(MyStrp.Sigcorr[0]>0.0 || MyStrp.Sigcorr[1]>0.0) nstripcalib+=1.0;
                    }
                  MyStrp.Trk += MyTrk->Index;
                  
                  if(seidP) delete seidP; seidP = 0;
                  if(seidN) delete seidN; seidN = 0;

                  break; // break out of loop
                         // faster - this should be ok as long as the
                         // same strip doesn't exist twice in a snarl
                }
            }
        }
      MSG("NtpMaker",Msg::kVerbose) <<"   "<<(int)stripsintrack<<" strips in the track of which "<<(int)nstripcalib<< " were calibrated."<<endl<<flush;
    }
  return;
}

//======================================================================

---