CondensedNtpModuleAtm.cxx

Go to the documentation of this file.

//$Id: CondensedNtpModuleAtm.cxx,v 1.3 2007/11/11 08:06:36 rhatcher Exp $
//
//CondensedNtpModuleAtm.cxx
//
//Module for making analysis trees from SR ntuples
//
//B. Rebel 10/2004
#include <cassert>

#include "AnalysisNtuples/Module/CondensedNtpModuleAtm.h"
#include "AnalysisNtuples/Module/ANtpRecoNtpManipulator.h"
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "JobControl/JobCModuleRegistry.h" // For JOBMODULE macro
#include "JobControl/JobCommand.h"
#include "Validity/VldTimeStamp.h" 
#include "CandNtupleSR/NtpSRRecord.h"
#include "CandNtupleSR/NtpSREvent.h"
#include "CandNtupleSR/NtpSREventSummary.h"
#include "CandNtupleSR/NtpSRSlice.h"
#include "CandNtupleSR/NtpSRTrack.h"
#include "CandNtupleSR/NtpSRShower.h"
#include "CandNtupleSR/NtpSRShieldStrip.h"
#include "CandNtupleSR/NtpSRShieldSummary.h"
#include "CandNtupleSR/NtpSRCosmicRay.h"
#include "CandNtupleSR/NtpSRDmxStatus.h"
#include "MCNtuple/NtpMCRecord.h"
#include "MCNtuple/NtpMCTruth.h"
#include "TruthHelperNtuple/NtpTHRecord.h"
#include "TruthHelperNtuple/NtpTHEvent.h"
#include "Record/RecCandHeader.h"
#include "CandData/CandRecord.h"
#include "CandData/CandHeader.h"
#include "RecoBase/ArrivalTime.h"
#include "RecoBase/LinearFit.h"
#include "Conventions/Detector.h"
#include "Conventions/SimFlag.h"
#include "Conventions/PlaneView.h"
#include "BeamDataNtuple/NtpBDLiteRecord.h"

#include "TFolder.h"
#include "TDirectory.h"
#include "TObjArray.h"
#include "TMatrix.h"
#include "TH1.h"
#include "TStopwatch.h"
#include "TRandom.h"
#include "Riostream.h"

#include <cassert>
#include <algorithm>

using namespace std;

ClassImp(CondensedNtpModuleAtm)

CVSID("$Id: CondensedNtpModuleAtm.cxx,v 1.3 2007/11/11 08:06:36 rhatcher Exp $");

// Declare this module to JobControll. Arguments are:
//  (1) The class name 
//  (2) The human-readable name 
//  (3) A short, human-readable description of what the module does
JOBMODULE(CondensedNtpModuleAtm, 
          "CondensedNtpModuleAtm",
          "A module used for making analysis tress from SR ntuples");

//......................................................................
CondensedNtpModuleAtm::CondensedNtpModuleAtm() : 
    fFileName("analysisNtuple.root"),
    fTreeName("analysisNtuple"),
    fNtpFile(0),
    fNtuple(0),
    fBeamTree(0),
    fDetector(1),
    fEndPlane(485),
    fEventInfo(0),
    fHeaderInfo(0),
    fTrackInfo(0),
    fTruthInfo(0),
    fCtr(0)
{
    MSG("JobC", Msg::kDebug) << "CondensedNtpModuleAtm::Constructor" << endl;

    fHeaderInfo = new ANtpHeaderInfo();
    fBeamInfo = new ANtpBeamInfo();
    fEventInfo = new ANtpEventInfo();
    fTrackInfo = new ANtpTrackInfoAtm();
    fTruthInfo = new ANtpTruthInfoAtm();
    
    fInfoFiller = new ANtpInfoObjectFillerBeam();

}

//----------------------------------------------------------------------
CondensedNtpModuleAtm::~CondensedNtpModuleAtm()
{

    MSG("JobC", Msg::kDebug) << "CondensedNtpModuleAtm::Destructor" << endl;

    if(fHeaderInfo) delete fHeaderInfo; 
    if(fBeamInfo) delete fBeamInfo;     
    if(fEventInfo) delete fEventInfo;   
    if(fTrackInfo) delete fTrackInfo;   
    if(fTruthInfo) delete fTruthInfo;   
    
}

//......................................................................
void CondensedNtpModuleAtm::BeginJob()
{
  // save the current working directory
  TDirectory* savedir  = gDirectory;  
  
    
  //create the new TFile for holding the electronics tree
  // this changes the value of gDirectory
  fNtpFile = new TFile(fFileName.c_str(),"RECREATE");  
    
  MSG("CondensedNtpModuleAtm", Msg::kDebug) << fFileName << " " << fTreeName 
                                            << " " << fNtpFile->GetName() << endl;
  
  
  // create TTree, these will attach themselves to the current 
  //working directory   
  fNtuple = new TTree(fTreeName.c_str(), "Analysis Tree");
  
  fNtuple->Branch("header.", "ANtpHeaderInfo", &fHeaderInfo, 64000, 2);
  fNtuple->Branch("beam.", "ANtpBeamInfo", &fBeamInfo, 64000, 2);
  fNtuple->Branch("event.", "ANtpEventInfo", &fEventInfo, 64000, 2);
  fNtuple->Branch("track.", "ANtpTrackInfoAtm", &fTrackInfo, 64000, 2);
  if(fDataType==1) fNtuple->Branch("truth.", "ANtpTruthInfoAtm", &fTruthInfo, 64000, 2);
  
  MSG("CondensedNtpModuleAtm", Msg::kDebug) << "got branches" << endl;

  fFailTree = new TTree("failures", "Failed Events");
  fFailTree->Branch("header.", "ANtpHeaderInfo", &fHeaderInfo, 64000, 2);
  fFailTree->Branch("multiple", &fIsMultiple, "multiple/I");
  fFailTree->Branch("validplane", &fValidPlaneFail, "validplane/I");
  fFailTree->Branch("failDeMux", &fFailDeMux, "failDeMux/I");
  
  // change back to current working directory before leaving constructor
  savedir->cd();   
  
  //load the vertical intensity information into some arrays
  //get the slant depths and the number of events for each slant depth
  //read in rock data
  // http://www.hep.umn.ed/~schubert/far/s2rock/slantdepth.data
  // 1-cos(zen); azi(degrees); slantdepth (gmcm**2) * cos(zen); 
  //    density (gm/cm**2)
  //  N.B., units for density are incorrectly labelled in the file
  ifstream densityFile(fRockMapFileName.c_str());
  if( !densityFile.is_open() ){
    cout << "could not open rock map - exit" << endl;
    exit(0);
  }
  
  int ir = 0;
  while( ir<1404 ){
    densityFile >> fPathCosZenDeg[ir] >> fPathAzimuthDeg[ir] >> fColDenCosZen[ir] >> fDensity[ir]; 
    ++ir;
  }
  
  //get the beam info tree
  //check to see if you need to calculate the total Protons on Target
  MSG("BeamTestModule", Msg::kInfo) << "beam info path = " << fBeamPath.c_str()
      << " tree name = " << fBeamTreeName.c_str() << endl;
  
//   if(fDataType < 1){
    
//     TFile *beamFile = new TFile(fBeamPath.c_str());
//     fBeamTree = dynamic_cast<TTree *>(beamFile->Get(fBeamTreeName.c_str()));
    
//     //build the index for the tree to select the utc timestamps
//     fBeamTree->BuildIndex("utc");
//     fBeamTree->SetBranchAddress("intensity", &fCurrentPOT);
//     fBeamTree->SetBranchAddress("hornCurrent", &fHornCurrent);
//     fBeamTree->SetBranchAddress("hpos", &fBeamHPos);
//     fBeamTree->SetBranchAddress("vpos", &fBeamVPos);
//     fBeamTree->SetBranchAddress("tgtpos", &fTargetPos);
//   }
    
  return;
}

//......................................................................
//here is where you would loop over the events in the snarl to fill the 
//analysis tree
JobCResult CondensedNtpModuleAtm::Ana(const MomNavigator *mom)
{
  MSG("CondensedNtpModuleAtm", Msg::kDebug) << "start ana" << endl;
  
  JobCResult result(JobCResult::kPassed);
  
  //get the records from MOM
  assert(mom);
    
  NtpSRRecord *record = dynamic_cast<NtpSRRecord *>(mom->GetFragment("NtpSRRecord"));
  NtpStRecord *stRecord = dynamic_cast<NtpStRecord *>(mom->GetFragment("NtpStRecord"));
  
  if(!record && !stRecord){
    MSG("CondensedNtpModuleAtm", Msg::kWarning) << "Could not get NtpSR or NtpSt Record from MOM" << endl;
    result.SetWarning().SetFailed();
    return result;
  }//end if no NtpSRRecord
  
  NtpMCRecord *mcRecord = 0;
  NtpTHRecord *thRecord = 0;
  
  mcRecord = dynamic_cast<NtpMCRecord *>(mom->GetFragment("NtpMCRecord"));
  thRecord = dynamic_cast<NtpTHRecord *>(mom->GetFragment("NtpTHRecord"));              
  
  //get the record for the beam tree
  NtpBDLiteRecord *bdRecord = dynamic_cast<NtpBDLiteRecord *>(mom->GetFragment("NtpBDLiteRecord"));
  if(!bdRecord) MSG("CondensedNtpModuleAtm", Msg::kDebug) << "NO BEAM RECORD" << endl;

  //fill the header information
  Detector::Detector_t det = Detector::kFar;
  if(fDetector<1) det = Detector::kNear;
  //get a sim flag
  SimFlag::SimFlag_t dataType = SimFlag::kData;
  if(fDataType) dataType = SimFlag::kMC;

  //instantiate a NtpHelper object to help you get the info you want
  ANtpRecoNtpManipulator *ntpManipulator = 0;
  if(record) 
    ntpManipulator = new ANtpRecoNtpManipulator(record, mcRecord, thRecord);
  else if(stRecord) 
    ntpManipulator = new ANtpRecoNtpManipulator(stRecord);
    
  VldTimeStamp timeStamp = stRecord ? stRecord->GetVldContext()->GetTimeStamp() : record->GetVldContext()->GetTimeStamp();

  MSG("CondensedNtpModuleAtm", Msg::kDebug) << timeStamp << endl;

//   TStopwatch stopwatch;
//   stopwatch.Start();

  fInfoFiller->SetStripArray(ntpManipulator->GetStripArray());
  fInfoFiller->FillHeaderInformation(ntpManipulator, det, dataType, fHeaderInfo);

//   stopwatch.Stop();
//   MSG("CondensedNtpModuleAtm", Msg::kDebug) << "took " << stopwatch.RealTime()
//                                         << " to fill header" << endl;
//   stopwatch.Reset();

//   stopwatch.Start();

  if(fHeaderInfo->year > 2004){
    if(bdRecord) fInfoFiller->FillBeamInformation(bdRecord, fBeamInfo);
    else fInfoFiller->FillBeamInformation(timeStamp, det, dataType, fBeamInfo);
  }

//   stopwatch.Stop();
//   MSG("CondensedNtpModuleAtm", Msg::kDebug) << "took " << stopwatch.RealTime()
//                                         << " to fill beam" << endl;
//   stopwatch.Reset();

  //set up which flags you want to use to determine the primary shower or track
  //a value of 0 for a flag means it will not be used
  ntpManipulator->SetPrimaryTrackCriteria(0,1,0); // nplanes, length, total pulse height
  ntpManipulator->SetPrimaryShowerCriteria(0,1); // nplanes, total pulse height
  
  //check that the current event passed demuxing
  bool passedDeMux = true; 
  fGoodTrack = true;
  fIsMultiple = 0;
  fValidPlaneFail = 0;
  fFailDeMux = 0;

  MSG("CondensedNtpModuleAtm", Msg::kDebug) << "run,snarl = " << fHeaderInfo->run 
                                            << "," << fHeaderInfo->snarl << " " 
                                            << ntpManipulator->GetSnarlEventSummary().ntrack 
                                            << " " << passedDeMux << endl;

  if(fHeaderInfo->detector == Detector::kFar){

    if(ntpManipulator->GetDmxStatusInfo().ismultimuon){
      passedDeMux = false;
      fIsMultiple = 1;
    }
    if(ntpManipulator->GetDmxStatusInfo().validplanesfail){
      passedDeMux = false;
      fValidPlaneFail = 1;
    }
    if(fMajorRelease < 18){

      if(ntpManipulator->GetDmxStatusInfo().vertexplanefail
         || ntpManipulator->GetDmxStatusInfo().nonphysicalfail) passedDeMux = false;
      
      UShort_t uValidPlanes = ntpManipulator->GetDmxStatusInfo().uvalidplanes;
      UShort_t uStrayPlanes = ntpManipulator->GetDmxStatusInfo().ustrayplanes;
      UShort_t vValidPlanes = ntpManipulator->GetDmxStatusInfo().vvalidplanes;
      UShort_t vStrayPlanes = ntpManipulator->GetDmxStatusInfo().vstrayplanes;
      Float_t ufrac = 0., vfrac = 0.;
    
      if(uValidPlanes>0) 
        ufrac = (1.*uStrayPlanes)/(1.*uValidPlanes);
      else ufrac = 1.;
      if(vValidPlanes>0) 
        vfrac = (1.*vStrayPlanes)/(1.*vValidPlanes);
      else vfrac = 1.;
      if(ufrac>=0.45) passedDeMux = false;
      else if( vfrac>=0.45 ) passedDeMux = false;
      else if(uStrayPlanes>=3&&uValidPlanes<20) passedDeMux = false;
      else if(vStrayPlanes>=3&&vValidPlanes<20) passedDeMux = false;
      else if(uStrayPlanes>=4&&uValidPlanes>=20 
              && uValidPlanes<=40) passedDeMux = false;
      else if(vStrayPlanes>=4&&vValidPlanes>=20 
              && vValidPlanes<=40) passedDeMux = false;
      else if(ufrac>=0.1 && uValidPlanes>=40) passedDeMux = false;
      else if(vfrac>=0.1 && vValidPlanes>=40) passedDeMux = false;
    
    }//end if before switch over to cambridge demuxing
  }//end if far detector

  //if the event isnt demuxed well, then forget it also forget it
  //if there arent any tracks in the event
  if(!passedDeMux || ntpManipulator->GetSnarlEventSummary().ntrack==0 || fEndPlane<248){
    fFailTree->Fill();
    return result;
  }

  //set up some pointers to typical NtpSR objects
  NtpSRTrack *ntpTrack = 0;
  NtpSREvent *ntpEvent = 0;
  NtpMCTruth *ntpMCTruth = 0;
  NtpMCStdHep *ntpMCStdHep = 0;
  NtpTHEvent *ntpTHEvent = 0;
  
  //loop over the events in the snarl
  for(Int_t i = 0; i < ntpManipulator->GetSnarlEventSummary().nevent; ++i){
    
    MSG("CondensedNtpModuleAtm", Msg::kDebug) << "on event " << i << endl;
    
    ResetTreeVariables();
    
    fEventInfo->index = i;

//     stopwatch.Start();
    //get event i.  the call to ANtpRecoNtpManipulator::GetEventInSnarl sets the 
    //NtpSREvent data member in that object to the current event so that
    //you can later call for the primary track, shower, etc.
    ntpManipulator->GetEventManipulator()->SetEventInSnarl(i);          
    ntpEvent = ntpManipulator->GetEventManipulator()->GetEvent();               
    if(ntpEvent) FillEventInformation(ntpManipulator,ntpEvent);
//     stopwatch.Stop();
//     MSG("CondensedNtpModuleAtm", Msg::kInfo) << "took " << stopwatch.RealTime()
//                                           << " to fill event" << endl;
//     stopwatch.Reset();
    
//     stopwatch.Start();
    //get the primary track for the event - if no track is present it 
    //returns 0
    fGoodTrack = false;
    ntpTrack = ntpManipulator->GetEventManipulator()->GetPrimaryTrack();                
    if(ntpTrack){
      fGoodTrack = true;
      FillTrackInformation(ntpTrack, ntpManipulator->GetStripArray(), 
                           ntpManipulator->GetCosmicRayInfo().azimuth,
                           ntpManipulator->GetCosmicRayInfo().ra,
                           ntpManipulator->GetCosmicRayInfo().dec);
      FillTrackMagneticBendingVariables(ntpTrack, ntpManipulator->GetStripArray());
      if(fDetector>0&&fGoodTrack) FillTrackTimingVariables(ntpTrack, ntpManipulator->GetStripArray());
    } 
//     stopwatch.Stop();
//     MSG("CondensedNtpModuleAtm", Msg::kInfo) << "took " << stopwatch.RealTime()
//                                           << " to fill track" << endl;
//     stopwatch.Reset();

//     stopwatch.Start();
    // Get best neu match from Truthhelper
    ntpTHEvent = ntpManipulator->GetMCManipulator()->GetNtpTHEvent(i);
    MSG("CondensedNtpModuleAtm", Msg::kDebug) << "start filling mc info " << ntpTHEvent->neumc << endl;
    if(ntpTHEvent) ntpMCTruth = ntpManipulator->GetMCManipulator()->GetNtpMCTruth(ntpTHEvent->neumc);
    if(ntpMCTruth){
      //get the NtpMCStdHep object - use the fact that i put the values i want in the parent 
      //place, which means i want the 0th entry
      if(fNewMC) ntpMCStdHep = ntpManipulator->GetMCManipulator()->GetNtpMCStdHep(0);
      FillMCInformation(ntpMCTruth, ntpMCStdHep);       
    }
//     stopwatch.Stop();
//     MSG("CondensedNtpModuleAtm", Msg::kInfo) << "took " << stopwatch.RealTime()
//                                           << " to fill truth" << endl;
//     stopwatch.Reset();

    if(fEventInfo->totalStrips<fTrackInfo->totalStrips) fGoodTrack = false;
    if(fTrackInfo->planes<1) fGoodTrack = false;
    if(fTrackInfo->totalStrips<1) fGoodTrack = false;
    
    MSG("CondensedNtpModuleAtm", Msg::kDebug) << "strips = " << fEventInfo->totalStrips 
                                              << "/" << fTrackInfo->totalStrips << " planes = " 
                                              << fTrackInfo->planes << " good track = " 
                                              << (int)fGoodTrack << endl;
    
    //got all the necessary sr info, now fill the ntuple
    if(!fGoodTrack) continue;
    
    fNtuple->Fill();
  }//end loop over events for this snarl
  
  ++fCtr;

  delete ntpManipulator;

  return result;
}

//......................................................................
void CondensedNtpModuleAtm::Help() 
{
  MSG("JobC", Msg::kInfo) 
    << "NearElectronicsCheck::Help\n"
    <<"CondensedNtpModuleAtm is a module which demultiplexes events "
    <<"in the far detector."
    << endl;
}

//----------------------------------------------------------------------
void CondensedNtpModuleAtm::EndJob() 
{

  MSG("CondensedNtpModuleAtm", Msg::kDebug) << "start end job method" << endl;  
  
  //save current working directory
  TDirectory *savedir = gDirectory;
  
  //cd to the "directory" of the file
  fNtpFile->cd();
  
  MSG("CondensedNtpModuleAtm", Msg::kInfo) << fNtuple->GetEntries() << endl;
  
  fNtuple->Write();
  fFailTree->Write();

  //go back to original working directory
  savedir->cd();  
  
  fNtpFile->Write();
  fNtpFile->Close();
  
  return;
}

//......................................................................
const Registry& CondensedNtpModuleAtm::DefaultConfig() const
{
  
  int itrue = 1;  // work around for lack of bool in registry
  int ifalse = 0; // work around for lack of bool in registry
  
  static Registry r;
  
  r.UnLockValues();
  
  r.Set("FileName",           "analysisTree.root");
  r.Set("RockMapFileName",    "density.txt");
  r.Set("TreeName",           "analysisTree");
  r.Set("BeamPath",           "/afs/fnal.gov/files/data/minos/d04/brebel/beamInfo.root");
  r.Set("BeamTreeName",       "beam");
  r.Set("EndPlane",           485);
  r.Set("Detector",           1);
  r.Set("DataType",           0);
  r.Set("NewMC",              itrue);
  r.Set("MajorRelease",       14);
  r.LockValues();

  itrue = ifalse; //quiet the compiler warnings

  return r;
}

//......................................................................
void CondensedNtpModuleAtm::Config(const Registry& r)
{
  int         tmpb;  // a temp bool. See comment under DefaultConfig...
  int         tmpi;  // a temp int.
  double      tmpd;  // a temp double.
  const char* tmps;  // a temp string
  
  if (r.Get("FileName",           tmps)) fFileName        =  tmps;
  if (r.Get("RockMapFileName",    tmps)) fRockMapFileName =  tmps;
  if (r.Get("TreeName",           tmps)) fTreeName        =  tmps;
  if (r.Get("BeamPath",           tmps)) fBeamPath        = tmps;
  if (r.Get("BeamTreeName",       tmps)) fBeamTreeName    = tmps;
  if (r.Get("EndPlane",           tmpi)) fEndPlane        =  tmpi;
  if (r.Get("Detector",           tmpi)) fDetector        =  tmpi;
  if (r.Get("DataType",           tmpi)) fDataType        =  tmpi;
  if (r.Get("NewMC",              tmpb)) fNewMC           =  tmpb;
  if (r.Get("MajorRelease",       tmpi)) fMajorRelease    =  tmpi;
  
  //quiet the compiler warnings
  tmpb = 0;
  tmpd = 1.;

  return;
}

//----------------------------------------------------------------------
//this is the method where you fill all variables related to a track
void CondensedNtpModuleAtm::FillTrackInformation(NtpSRTrack *ntpTrack, 
                                                 TClonesArray *stripArray, 
                                                 double azimuth, double ra, double dec)
{
  MSG("CondensedNtpModuleAtm", Msg::kDebug) << "filling track info" << endl;
  
  //fill the basics first
  fInfoFiller->FillTrackInformation(ntpTrack, fTrackInfo);

  NtpSRStrip *ntpStrip = 0;
  
  Int_t plane = 0, strip = 0, numDoubleEndedStrips = 0, prevPlane = -1;
  float xPos = 0., yPos = 0.;
  Float_t radius = 0.;
  
  Int_t index = 0;

  //stuff to find the slant depth
  double vec1[39],vec2[36];
  int index1, index2, index3;
  
  for(int i=0;i<39;++i){
    vec1[i] = i*0.02;
    //     cout << vec1[i] << endl;
    if(i<36){
      vec2[i] = i*10. + 5.;
      //       cout << vec2[i] << endl;
    }
  }

  // Kashahara picks events according to where they come from
  //   she uses altitude, not zenith angle
  index1 = TMath::BinarySearch(39,vec1,1.-(-fTrackInfo->dcosYVtx))+1;
  index2 = TMath::BinarySearch(36,vec2,azimuth)+1;
  index3 = (index1)*36 +(index2);
  fTrackInfo->slantDepth = fDensity[index3]*0.01; //density is in 10^5g/cm^2, ie km.w.e, 
                                                  //0.01 converts to m.w.e 10^3m/km/(10^5g/cm^2/km.w.e)
  
  fTrackInfo->azimuth = azimuth;

  fTrackInfo->ra = ra;
  fTrackInfo->dec = dec;
  fTrackInfo->trackLikePlanes = ntpTrack->plane.ntrklike;
  Int_t nVPlanes = ntpTrack->plane.nv;
  Int_t nUPlanes = ntpTrack->plane.nu;
  fTrackInfo->invBeta = ntpTrack->time.cdtds;
  fTrackInfo->planesIn10Meter = 0;
  fTrackInfo->planesIn15Meter = 0;
  fTrackInfo->planesIn20Meter = 0;
  fTrackInfo->planesIn25Meter = 0;
  fTrackInfo->planesIn30Meter = 0;
  fTrackInfo->planesIn35Meter = 0;
  fTrackInfo->planesIn40Meter = 0;
  
  //check if the zenith angle indicates the right direction for the near detector. 
  //if not, swap the ends and adjust the azimuth
  if(fTrackInfo->dcosYVtx > 0. && fDetector<1){
    
    fTrackInfo->vtxX = ntpTrack->end.x; 
    fTrackInfo->vtxY = ntpTrack->end.y; 
    fTrackInfo->vtxZ = ntpTrack->end.z;                 
    fTrackInfo->endX = ntpTrack->vtx.x; 
    fTrackInfo->endY = ntpTrack->vtx.y; 
    fTrackInfo->endZ = ntpTrack->vtx.z; 
    fTrackInfo->dcosXVtx = -ntpTrack->end.dcosx;
    fTrackInfo->dcosYVtx = -ntpTrack->end.dcosy;
    fTrackInfo->dcosZVtx = -ntpTrack->end.dcosz;
    fTrackInfo->dcosXEnd = -ntpTrack->vtx.dcosx;
    fTrackInfo->dcosYEnd = -ntpTrack->vtx.dcosy;
    fTrackInfo->dcosZEnd = -ntpTrack->vtx.dcosz;
    
    const double rotMatrix[] = { 0.894507011,0,0.447053919,
                                 0,          1,          0,
                                 -0.447053919,0,0.894507011};
    
    double dcosx_ideal, dcosy_ideal, dcosz_ideal;
    
    //get the azimuth from AstroUtill
    dcosx_ideal = (rotMatrix[0]*fTrackInfo->dcosXVtx
                   + rotMatrix[1]*fTrackInfo->dcosYVtx
                   + rotMatrix[2]*fTrackInfo->dcosZVtx);
    dcosy_ideal = (rotMatrix[3]*fTrackInfo->dcosXVtx
                   + rotMatrix[4]*fTrackInfo->dcosYVtx
                   + rotMatrix[5]*fTrackInfo->dcosZVtx);
    dcosz_ideal = (rotMatrix[6]*fTrackInfo->dcosXVtx
                   + rotMatrix[7]*fTrackInfo->dcosYVtx
                   + rotMatrix[8]*fTrackInfo->dcosZVtx);
    
    fTrackInfo->azimuth = 180.*TMath::ATan2(-dcosx_ideal,dcosz_ideal)/TMath::Pi();
    if(fTrackInfo->azimuth < 0.) fTrackInfo->azimuth += 360.;
  }
  
  if(ntpTrack->fit.ndof == 0 || ntpTrack->momentum.qp== 0.){
    fGoodTrack = false;
    return;
  }
  
  //the value of trk.fidvtx.dz is not valid for events with the full detector
  //so hack up my own value of fiducialDz
  Float_t fullDetEndZ = 0.0594*485. + 1.;
  
  fTrackInfo->fiducialVtxDz = ntpTrack->fidvtx.dz;
  if(fEndPlane>248) 
    fTrackInfo->fiducialVtxDz = TMath::Min(fullDetEndZ - fTrackInfo->vtxZ, fTrackInfo->vtxZ); 
  fTrackInfo->fiducialVtxDr = ntpTrack->fidvtx.dr;
  
  fTrackInfo->fiducialEndDz = ntpTrack->fidend.dz;
  if(fEndPlane>248) 
    fTrackInfo->fiducialEndDz = TMath::Min(fullDetEndZ - fTrackInfo->endZ, fTrackInfo->endZ); 
  fTrackInfo->fiducialEndDr = ntpTrack->fidend.dr;
  
  fTrackInfo->uvAsymmetry = 1.;
  if(nVPlanes+nUPlanes>0)
    fTrackInfo->uvAsymmetry = TMath::Abs((1.*nUPlanes-1.*nVPlanes)/(1.*nVPlanes+1.*nUPlanes));
  
  fTrackInfo->signalUseFraction = 0.;
  if(fEventInfo->pulseHeight>0.) 
    fTrackInfo->signalUseFraction = fTrackInfo->pulseHeight/(fEventInfo->pulseHeight);
  
  fTrackInfo->planeUseFraction = 0.;
  if(fEventInfo->planes>0) 
    fTrackInfo->planeUseFraction = 1.*fTrackInfo->trackLikePlanes/(1.*fEventInfo->planes);
  
  Int_t numStrips = 0;
  Int_t numDigits = 0;
  
  fTrackInfo->alternateChi2 = 0.;
  double deltaTPos = 0.;
  double stripPosError = 0.041/TMath::Sqrt(12.);

  fTrackInfo->impactParameter = 20.0;
  fTrackInfo->totalStrips = 0;
  //does the pathlength of the track make sense?
  if(fTrackInfo->length<50.){
    
    
    //loop over all strips in the track
    numStrips = (int)ntpTrack->nstrip;
    numDoubleEndedStrips = 0;
    
    MSG("CondensedNtpModuleAtm", Msg::kDebug) << "start looping over track strips" << endl;
    
    for(Int_t ns = 0; ns < numStrips; ++ns){
      
      //get the index for this strip
      index = ntpTrack->stp[ns];
      
      //get the NtpSRStrip object
      ntpStrip = dynamic_cast<NtpSRStrip *>(stripArray->At(index));
      plane = (int)ntpStrip->plane;
      strip = (int)ntpStrip->strip;
      numDigits = (int)ntpStrip->ndigit;
      
      xPos = ntpTrack->stpx[ns];
      yPos = ntpTrack->stpy[ns];
     
      if((int)ntpStrip->planeview == PlaneView::kU)
        deltaTPos = ntpTrack->stpu[ns]-ntpStrip->tpos;
      else if((int)ntpStrip->planeview == PlaneView::kV)
        deltaTPos = ntpTrack->stpv[ns]-ntpStrip->tpos;

      MSG("CondensedNtpModuleAtm", Msg::kDebug) << ntpTrack->stpu[ns] << " " << ntpStrip->tpos
                                                  << " " << ntpTrack->stpv[ns] << " " << deltaTPos << endl;

      fTrackInfo->alternateChi2 += deltaTPos*deltaTPos/(stripPosError*stripPosError);

      if(numDigits >= 2){
        
        ++numDoubleEndedStrips;
        
        radius = TMath::Sqrt(yPos*yPos + xPos*xPos);
        if(radius<fTrackInfo->impactParameter)
          fTrackInfo->impactParameter = TMath::Sqrt(yPos*yPos + xPos*xPos);
        
        //get values to determine the number of plances crossed by the track 
        //in concentric circles in the detector from 0.3 m out to 
        //4.0 m.  the points come in order, so i dont have to test each one 
        //to make sure it is before or after the previous points
        if(radius>0.3 && radius<1.0 && plane != prevPlane) ++fTrackInfo->planesIn10Meter;
        if(radius>0.3 && radius<1.5 && plane != prevPlane) ++fTrackInfo->planesIn15Meter;
        if(radius>0.3 && radius<2.0 && plane != prevPlane) ++fTrackInfo->planesIn20Meter;
        if(radius>0.3 && radius<2.5 && plane != prevPlane) ++fTrackInfo->planesIn25Meter;
        if(radius>0.3 && radius<3.0 && plane != prevPlane) ++fTrackInfo->planesIn30Meter;
        if(radius>0.3 && radius<3.5 && plane != prevPlane) ++fTrackInfo->planesIn35Meter;
        if(radius>0.3 && radius<4.0 && plane != prevPlane) ++fTrackInfo->planesIn40Meter;
        
        ++fTrackInfo->totalStrips;
        
      } // end if the strip had double sided readout
      
      prevPlane = plane;
      
    } // end loop over strips in the track
    
    fTrackInfo->twoEndStripFraction = 0.;
    if(numStrips>0){
      fTrackInfo->twoEndStripFraction = (1.*numDoubleEndedStrips)/(1.*numStrips);
      fTrackInfo->alternateChi2 /= (1.*numStrips);
    }
  }//end if the track length is reasonable
  else fGoodTrack = false;
  
  MSG("CondensedNtpModuleAtm", Msg::kDebug) << "end fill track info - good track = " 
                                            << (int)fGoodTrack << endl;
  
  return;
}

//----------------------------------------------------------------------
void CondensedNtpModuleAtm::FillTrackTimingVariables(NtpSRTrack *ntpTrack, 
                                                     TClonesArray *stripArray)
{    
  MSG("CondensedNtpModuleAtm", Msg::kDebug) << "in FillTrackTimingVariables" << endl;
  
  NtpSRStrip *ntpStrip = 0;
  
  Double_t time[2] = {0.};
  Double_t xPos = 0., yPos = 0., zPos = 0.;
  Double_t yPosition[5000] = {0.};
  Double_t yPosition1[5000] = {0.};
  Double_t weight[5000] = {0.};
  Double_t weight2[5000] = {0.};
  Double_t pathLength[5000] = {0.};
  Double_t weight1[5000] = {0.};
  Double_t timeAv[5000] = {0.};
  Double_t timeAv1[5000]  = {0.};
  Int_t arrayCtr = 0, betaCtr = 0, actr = 0, plane = 0, strip = 0, index = 0;
  Double_t minTime = 1.e20;     
  
  //loop over all strips in the track
  Int_t numStrips = (int)ntpTrack->nstrip;
  Int_t numDigits = 0;
  Double_t smear = 0.;

  MSG("CondensedNtpModuleAtm", Msg::kDebug) << "fill track timing" << endl;
  
  for(Int_t ns = 0; ns < numStrips; ++ns){

    if(fHeaderInfo->dataType == (int)SimFlag::kMC) 
      smear = gRandom->Gaus(0., 0.75);
    else smear = 0.;

    //get the index for this strip
    index = ntpTrack->stp[ns];
    
    //get the NtpSRStrip object
    ntpStrip = dynamic_cast<NtpSRStrip *>(stripArray->At(index));
    plane = (int)ntpStrip->plane;
    strip = (int)ntpStrip->strip;
    numDigits = (int)ntpStrip->ndigit;
    
    xPos = ntpTrack->stpx[ns];
    yPos = ntpTrack->stpy[ns];
    zPos = ntpTrack->stpz[ns];
    
    if(numDigits == 2){
      
      time[0] = ntpTrack->stpt0[ns];
      time[1] = ntpTrack->stpt1[ns];
      
      if(time[0] < minTime) minTime = time[0];
      if(time[1] < minTime) minTime = time[1];
      
      //have a point for each end of the strip
      timeAv[arrayCtr] = 1.e9*time[0] + smear;
      timeAv[arrayCtr+1] = 1.e9*time[1] + smear;
      pathLength[arrayCtr] = fTrackInfo->length - ntpTrack->stpds[ns];
      pathLength[arrayCtr+1] = fTrackInfo->length - ntpTrack->stpds[ns];
      
      //stuff to get timing information
      yPosition[arrayCtr] = yPos;
      weight[arrayCtr] = GetTimeWeight(ntpStrip->ph0.sigcor);
      if (weight[arrayCtr]>0.4) weight[arrayCtr]=0.4;
      yPosition[arrayCtr+1] = yPos;
      weight[arrayCtr+1] = GetTimeWeight(ntpStrip->ph1.sigcor);
      if (weight[arrayCtr+1]>0.4) weight[arrayCtr+1]=0.4;
      MSG("CondensedNtpModuleAtm", Msg::kDebug) << weight[arrayCtr] << " " 
                                                << weight[arrayCtr+1] << endl;
      arrayCtr += 2;
      
      
    } // end if the strip had double sided readout
  } // end loop over strips in the track

  //subtract the minimum time and copy the weights into another array to do the fit to 
  //a flat line later
  
  for(Int_t ti = 0; ti< arrayCtr; ++ti){
    timeAv[ti] -= 1.e9*minTime;
    weight2[ti] = weight[ti];
  }
  
  //cout << "got strip info" << endl;
  
  //array for fit parameters
  Double_t param[] = {-10000., -10000.};
  Double_t eparam[] = {-10000., -10000.};
  
  //find the inverse beta value
  betaCtr = arrayCtr;
  FitMinimizingResidual(betaCtr, pathLength, timeAv, 
                        weight, param, eparam);
  fTrackInfo->invBeta = param[1]*0.3;
  if(betaCtr > 2) fTrackInfo->invBetaChi2 = eparam[0] /(1.*(betaCtr-2));
  
  //find the chi^2 for a fit to a straight line so you 
  //can see if the sloped fit is really better
  FitMinimizingResidual(betaCtr, pathLength, timeAv, 
                        weight2, param, eparam, false);
  if(betaCtr > 1) fTrackInfo->flatLineChi2 = eparam[0] / (1.*betaCtr-1);

  //now loop over the points in the array and fill new arrays 
  //containing only those points that have a nonzero weight
  actr = 0;
  MSG("CondensedNtpModuleAtm", Msg::kDebug) << fHeaderInfo->run << " " 
                                            << fHeaderInfo->subRun
                                            << " " << fHeaderInfo->snarl 
                                            << " " << arrayCtr << endl;
  for(Int_t ii = 0; ii < arrayCtr; ++ii){
    MSG("CondensedNtpModuleAtm", Msg::kDebug) << " " << ii << "/" << arrayCtr << " " 
                                              << yPosition[ii] << " " << weight[ii] 
                                              << " " << timeAv[ii] << " " 
                                              << actr << endl;
    if(weight[ii]>0.){
      yPosition1[actr] = yPosition[ii];
      weight1[actr] = weight[ii];
      timeAv1[actr] = timeAv[ii];
      ++actr;
    }
  }//end loop to select good points along the track
  
  //get the fit for the time in the y direction
  FitMinimizingResidual(actr, yPosition1, timeAv1, weight1, param, eparam);
  fTrackInfo->timeSlopeY = param[1];
  fTrackInfo->timeSlopeYChi2 = eparam[0];
  
  if(actr > 2) fTrackInfo->timeSlopeYChi2 /= 1.*(actr-2);  
  
  MSG("CondensedNtpModuleAtm", Msg::kDebug) << "actr = " << actr 
                                            << "/" << arrayCtr 
                                            << " goodtrack = " 
                                            << (int)fGoodTrack << endl;
  
  if(1.*actr < 0.5*arrayCtr) fGoodTrack = false;
  
  return;
}

//----------------------------------------------------------------------
//this method is to fill information related to the magnetic measurements of
//the track.  
void CondensedNtpModuleAtm::FillTrackMagneticBendingVariables(NtpSRTrack *ntpTrack,  
                                                              TClonesArray *stripArray)
{
    
  MSG("CondensedNtpModuleAtm", Msg::kDebug) << "in FillTrackMagneticBendingVariables" 
                                            << endl;
  
  NtpSRStrip *ntpStrip = 0;
  
  
  //find the direction cosines for the straight line between the vertex and end points
  //of the track
  Float_t deltaX = ntpTrack->end.x - ntpTrack->vtx.x;
  Float_t deltaY = ntpTrack->end.y - ntpTrack->vtx.y;
  Float_t deltaZ = ntpTrack->end.z - ntpTrack->vtx.z;
  Float_t deltaS = TMath::Sqrt(deltaX*deltaX + deltaY*deltaY + deltaZ*deltaZ);
  Float_t dcosX = deltaX/deltaS;
  Float_t dcosY = deltaY/deltaS;
  Float_t dcosZ = deltaZ/deltaS;
  Float_t xPos = 0.;
  Float_t yPos = 0.;
  Float_t zPos = 0.;
  Float_t linXPos = 0.;
  Float_t linYPos = 0.;
  
  //loop over the strips in the track to find the sagitta
  fTrackInfo->sagitta = 0.;
  fTrackInfo->netDistFromLinearFit = 0.;
  fTrackInfo->meanDistFromLinearFit = 0.;
  fTrackInfo->rmsDistFromLinearFit = 0.;
  fTrackInfo->zeroCurvatureChi2 = 0.;

  if(deltaZ == 0. || deltaX == 0. || deltaY == 0.){
    fGoodTrack = false;
    return;
  }
  
  TH1F dist("dist", "", 500, 0., 5.);
  
  Float_t curDistFromLinearFit = 0.;
  
  //loop over all strips in the track
  Int_t numStrips = (int)ntpTrack->nstrip;
  Int_t index = 0;
  double stripPosError = 0.041/TMath::Sqrt(12.);

  for(Int_t ns = 0; ns < numStrips; ++ns){
    
    //get the index for this strip
    index = ntpTrack->stp[ns];
    
    //get the NtpSRStrip object
    ntpStrip = dynamic_cast<NtpSRStrip *>(stripArray->At(index));
    
    xPos = ntpTrack->stpx[ns];
    yPos = ntpTrack->stpy[ns];
    zPos = ntpTrack->stpz[ns];
    
    linXPos = ntpTrack->vtx.x + (dcosX/dcosZ) * (zPos - ntpTrack->vtx.z);
    linYPos = ntpTrack->vtx.y + (dcosY/dcosZ) * (zPos - ntpTrack->vtx.z);
    deltaX = xPos - linXPos;
    deltaY = yPos - linYPos;

    fTrackInfo->zeroCurvatureChi2 += (deltaX*deltaX + deltaY*deltaY)/(stripPosError*stripPosError);
    
    curDistFromLinearFit = TMath::Sqrt(deltaX*deltaX + deltaY*deltaY);
    dist.Fill(curDistFromLinearFit);
    
    if(curDistFromLinearFit > fTrackInfo->sagitta) fTrackInfo->sagitta = curDistFromLinearFit;
    
  }
  
  if(numStrips>1) fTrackInfo->zeroCurvatureChi2 /= 1.*(numStrips-1);
  fTrackInfo->meanDistFromLinearFit = dist.GetMean();
  fTrackInfo->rmsDistFromLinearFit = dist.GetRMS();
  
  //now to figure out the net pt kick to the muon
  
  return;
}

//----------------------------------------------------------------------
//this is the method where you fill all variables related to an event
void CondensedNtpModuleAtm::FillEventInformation(ANtpRecoNtpManipulator *ntpManipulator,
                                                 NtpSREvent *ntpEvent)
{
  MSG("CondensedNtpModuleAtm", Msg::kDebug) << "in FillEventInformation" << endl;
  
  fInfoFiller->FillEventInformation(ntpManipulator,
                                    ntpEvent, 
                                    fEventInfo);
  
  return;
}

//----------------------------------------------------------------------
//this is the method where you fill all variables related to an event
void CondensedNtpModuleAtm::FillMCInformation(NtpMCTruth *ntpMCTruth, NtpMCStdHep *ntpMCStdHep)
{
  MSG("CondensedNtpModuleAtm", Msg::kDebug) << "in FillMCInformation" << endl;
  
  const double rotMatrix[] = { 0.894507011,0,0.447053919,
                               0,          1,          0,
                              -0.447053919,0,0.894507011};
  
  double dcosx_ideal, dcosy_ideal, dcosz_ideal;
  
  fInfoFiller->FillMCTruthInformation(ntpMCTruth, fTruthInfo);
  
  if(fTruthInfo->nuEnergy>0.) fTruthInfo->baseLine = 0.001*CheapBaseline(-ntpMCTruth->p4neu[1]
                                                                         /TMath::Abs(fTruthInfo->nuEnergy));
  
  //get the azimuth from AstroUtill
  dcosx_ideal = (rotMatrix[0]*fTruthInfo->leptonDCosX
                 + rotMatrix[1]*fTruthInfo->leptonDCosY
                 + rotMatrix[2]*fTruthInfo->leptonDCosZ);
  dcosy_ideal = (rotMatrix[3]*fTruthInfo->leptonDCosX
                 + rotMatrix[4]*fTruthInfo->leptonDCosY
                 + rotMatrix[5]*fTruthInfo->leptonDCosZ);
  dcosz_ideal = (rotMatrix[6]*fTruthInfo->leptonDCosX
                 + rotMatrix[7]*fTruthInfo->leptonDCosY
                 + rotMatrix[8]*fTruthInfo->leptonDCosZ);
  fTruthInfo->azimuth = 180.*TMath::ATan2(-dcosx_ideal,dcosz_ideal)/TMath::Pi();
  if(fTruthInfo->azimuth < 0.) fTruthInfo->azimuth += 360.;

  //use the NtpMCStdHep object to figure out the weight for this event based on the 
  //fact that the new MC picks events at the surface.  we plugged the value into
  //the vtx Y location in the stdhep common bloc, gminos multiplies it by 10^-3
  //because it thought it was converting mm to m
  if(ntpMCStdHep){
    fTruthInfo->weight = 1.e3*ntpMCStdHep->vtx[1];
  }

  return;
}

//-------------------------------------------------------------------
void CondensedNtpModuleAtm::ResetTreeVariables()
{
  fEventInfo->Reset();
  fTrackInfo->Reset();
  fTruthInfo->Reset();
  fGoodTrack = false;
  
  return;
}

//-------------------------------------------------------------------
Double_t CondensedNtpModuleAtm::GetTimeWeight(Float_t ph)
{
  //   cout << "in GetTimeWeight" << endl;
  
  Double_t npe = ph/60.;
  if (npe<1.) npe=1.;
  
  // cout << ArrivalTime_Weight(npe) << endl;
  
  return ArrivalTime_Weight(npe);
}

//-------------------------------------------------------------------
Double_t CondensedNtpModuleAtm::ArrivalTime_Weight(Double_t npe) const
{
  //   cout << "in ArrivalTime_Weight" << endl;
  
  Double_t weight = 1./ArrivalTime_Uncertainty(npe);
  
  //    //cout << weight << endl;
  
  return weight*weight;
}

//--------------------------------------------------------------------
Double_t CondensedNtpModuleAtm::ArrivalTime_Uncertainty(Double_t npe) const
{
  //   cout << "in ArrivalTime_Uncertainty" << endl;
  
  Double_t sigma = 0.; // sigma is uncertainty in arrival of first p.e.
  
  // a simple monte carlo was done to calculate the uncertainty in the
  // arrival time of the first p.e. as a function of # of p.e.
  //
  // the only assumption made was about the shape of the pulse, a 2 ns risetime
  // and an 8 ns decay time of the fluor
  //
  // rlee feb 2002
  //
  
  if(npe<10.) sigma = 0.90+8.0/npe;
  else sigma = 5.0*exp(-0.5*log(npe));
  
  //  cout << sigma << endl;
  
  return sigma;
}


//--------------------------------------------------------------------
void CondensedNtpModuleAtm::LinearFit_Weighted(const Int_t n, 
                                               const Double_t *x, 
                                               const Double_t *y, 
                                               const Double_t *w,
                                               Double_t *parm, Double_t *eparm)
{
  //   cout << "in LinearFit_Weighted with errors" << endl;
  
  Double_t sumx=0.;
  Double_t sumx2=0.;
  Double_t sumy=0.;
  Double_t sumy2=0.;
  Double_t sumxy=0.;
  Double_t sumw=0.;
  
  for(Int_t i=0; i<n; ++i){
    sumx += x[i]*w[i];
    sumx2 += x[i]*x[i]*w[i];
    sumy += y[i]*w[i]; 
    sumy2 += y[i]*y[i]*w[i];
    sumxy += x[i]*y[i]*w[i];
    sumw += w[i];
    
    //     cout << x[i] << " " << y[i] << " " << w[i] << endl;
  }
  
  //   cout << sumx2 << " " << sumw << " " << sumx << " " << sumxy << " " << sumy
  //        << endl;
  
  if(sumx2*sumw-sumx*sumx>0.&&sumx2<1.e5){
    parm[0] = (sumy*sumx2-sumx*sumxy)/(sumx2*sumw-sumx*sumx);
    parm[1] = (sumxy*sumw-sumx*sumy)/(sumx2*sumw-sumx*sumx);
    
    eparm[0] = TMath::Sqrt(sumx2*(sumx2*sumw-sumx*sumx))/(sumx2*sumw-sumx*sumx);
    eparm[1] = TMath::Sqrt(sumx2*sumw-sumx*sumx)/(sumx2*sumw-sumx*sumx);
  }
  
  return;
}

//--------------------------------------------------------------------
void CondensedNtpModuleAtm::WeightedAverage(const Int_t n, 
                                            const Double_t *x, 
                                            const Double_t *w,
                                            Double_t *parm, Double_t *eparm)
{
  //   cout << "in LinearFit_Weighted with errors" << endl;

  double sumWeights = 0.;
  double sumXWeights = 0.;
  //find the weighted average of the points
  for(int i = 0; i < n; ++i){
    sumXWeights += x[i]*w[i];
    sumWeights += w[i];
  }
  
  if(sumWeights > 0.) parm[0] = sumXWeights/sumWeights;
  else parm[0] = -9999.;
  parm[1] = 0.;
  
  eparm[0] = 0.;
  eparm[1] = 0.;

  
  return;
}

//-------------------------------------------------------------------
void CondensedNtpModuleAtm::FitMinimizingResidual(Int_t &n, Double_t *x, 
                                                  Double_t *y, 
                                                  Double_t *w,
                                                  Double_t *param, 
                                                  Double_t *eparam, bool findSlope)
{
  //   cout << "in LinearFit_MinimizeResid" << endl;
  
  Int_t nremove = 0;
  Double_t maxresid = 11.;
  Double_t resid = 0.;
  Double_t parm[2] = {0.};
  Double_t eparm[2] = {0.};
  Double_t timefitchi2 = 0.;
  while(maxresid>10.
        &&n-nremove>4
        &&(1.*nremove)/(1.*n)<0.2){
    
    timefitchi2=0.;
    if(findSlope) LinearFit_Weighted(n,x,y,w,parm,eparm);
    else WeightedAverage(n,y,w,parm,eparm);
    maxresid = 0.;
    Int_t imaxresid=0;
    for (Int_t i=0; i<n; i++) {
      if (w[i]>0.) {
        if(findSlope) resid = (parm[1]*x[i]+parm[0]-y[i]);
        else resid = parm[0] - y[i];
        timefitchi2 += w[i]*resid*resid;
        if (TMath::Abs(resid)>maxresid) {
          maxresid = TMath::Abs(resid);
          imaxresid = i;
        }
      }
    }
    if (maxresid>10.) {
      w[imaxresid] = 0.;
      nremove++;
    }
  }//end loop to minimize residuals
  
  //get the chi^2
  timefitchi2 = 0.;
  for (Int_t i=0; i<n; i++) {
    if (w[i]>0.) {
      if(findSlope) resid = (parm[1]*x[i]+parm[0]-y[i]);
      else resid = parm[0] - y[i];
      timefitchi2 += w[i]*resid*resid;
    }
  }
  
  
  n -= nremove;
  param[0] = parm[0];
  param[1] = parm[1];
  eparam[0] = timefitchi2;
  eparam[1] = eparm[1];
  
  return;
}

//----------------------------------------------------------------
void CondensedNtpModuleAtm::FindLinearFit(Double_t *x, Double_t *y, 
                                          Double_t *weight, 
                                          Int_t nPoints, Double_t &a1, 
                                          Double_t &a2, Double_t &chiSq)
{
  
  //use method of determinants to do the least squares fit
  //this is from bevington chapter 6
  //
  //f_1(x) = 1
  //f_2(x) = x
  
  //need to find the sums
  //SIGMA y_i*f_1(x_i)/weight_i
  //SIGMA y_i*f_2(x_i)/weight_i
  //SIGMA f_1(x_i)*f_1(x_i)/weight_i
  //SIGMA f_1(x_i)*f_2(x_i)/weight_i
  //SIGMA f_2(x_i)*f_1(x_i)/weight_i
  //SIGMA f_2(x_i)*f_2(x_i)/weight_i
  //
  //some of these are the same, ie f_2*f_1 = f_1*f_2 so only need
  //3 instead of 4 variables
  
  Double_t yf1divWeight = 0.;
  Double_t yf2divWeight = 0.;
  Double_t f1f1divWeight = 0.;
  Double_t f1f2divWeight = 0.;
  Double_t f2f2divWeight = 0.;
  
  for(Int_t i=0; i<nPoints; i++){
    //     cout << y[i] << " " << x[i] << " " << weight[i] << endl;
    
    yf1divWeight += y[i]/(weight[i]*weight[i]);
    yf2divWeight += (y[i]*x[i])/(weight[i]*weight[i]);
    f1f1divWeight += 1./(weight[i]*weight[i]);
    f1f2divWeight += x[i]/(weight[i]*weight[i]);
    f2f2divWeight += (x[i]*x[i])/(weight[i]*weight[i]);
  }
  
  //make a bunch of matricies to hold these values
  TMatrix deltaMatrix(2,2);
  TMatrix a1Matrix(2,2);
  TMatrix a2Matrix(2,2);
  
  //fill the matricies
  deltaMatrix(0,0) = f1f1divWeight;
  deltaMatrix(0,1) = f1f2divWeight;
  deltaMatrix(1,0) = f1f2divWeight;
  deltaMatrix(1,1) = f2f2divWeight;
  
  a1Matrix(0,0) = yf1divWeight;
  a1Matrix(0,1) = f1f2divWeight;
  a1Matrix(1,0) = yf2divWeight;
  a1Matrix(1,1) = f2f2divWeight;
  
  a2Matrix(0,0) = f1f1divWeight;
  a2Matrix(0,1) = yf1divWeight;
  a2Matrix(1,0) = f1f2divWeight;
  a2Matrix(1,1) = yf2divWeight;
  
  a1 = -10000.;
  a2 = -10000.;
  Double_t delta = deltaMatrix.Determinant();
  
  //   cout << "delta = " << delta << endl;
  
  if(delta != 0.){
    a1 = a1Matrix.Determinant()/delta;
    a2 = a2Matrix.Determinant()/delta;
  }
  
  //find the chi^2 value for the fit
  chiSq = 0.;
  Double_t chi = 0.;
  for(Int_t j = 0; j<nPoints; j++){
    
    chi = (y[j] - (a1 + a2*x[j]))/weight[j];
    
    //     cout << chi*weight[j] << endl;
    
    chiSq += chi*chi;
    
  }  
  
  //   cout << "chi^2 = " << chiSq << " " << a1 << " " << a2 << endl;
  
  //get the chiSq per degree of freedom
  
  chiSq /= (1.*nPoints - 2.);
  
  return;
}

//----------------------------------------------------------------
//method to do a simple straight line fit to track in 2D and return
//a chi^2 for that fit.  the fit is found by using the end points
//of the track only
//y = a1 + a2*x
void CondensedNtpModuleAtm::FindChi2(Double_t *x, Double_t *y, 
                                     Double_t *weight, 
                                     Int_t nPoints, Double_t &a1, 
                                     Double_t &a2, Double_t &chiSq)
{
  
  double vtxX = x[0];
  double vtxY = y[0];
  double endX = x[nPoints-1];
  double endY = y[nPoints-1];
  
  a2 = endY-vtxY;
  if(TMath::Abs(vtxX-endX)>0.) a2 /= (endX-vtxX);
  else{
    cout << "cannot find slope for n = " << nPoints << " points " << endl;
    for(Int_t i = 0; i<nPoints; ++i) 
      cout << x[i] << " " << y[i] << " " << weight[i] << endl;
  }
  
  a1 = endY-a2*endX;
  
  //find the chi^2 value for the fit
  chiSq = 0.;
  Double_t chi = 0.;
  for(Int_t j = 0; j<nPoints; j++){
    
    chi = (y[j] - (a1 + a2*x[j]))/weight[j];
    
    //     cout << y[j] << " " << a1+a2*x[j] << " " << chiSq << endl;
    
    chiSq += chi*chi;
    
  }  
  
  //   cout << "chi^2 = " << chiSq << " " << a1 << " " << a2 << endl;
  
  //get the chiSq per degree of freedom
  
  if(nPoints>2) chiSq /= 1.*(nPoints-2);
  
  return;
}

//----------------------------------------------------------------
//method to do a simple straight line fit to track in 2D and return
//a the sum of the deviations for that fit.  the fit is found by using 
//the end points of the track only.
//y = a1 + a2*x
void CondensedNtpModuleAtm::FindStraightLineDeviation(Double_t *x, 
                                                      Double_t *y, 
                                                      Int_t nPoints, 
                                                      Double_t &a1, 
                                                      Double_t &a2, 
                                                      Double_t &deviation)
{
  
  double vtxX = x[0];
  double vtxY = y[0];
  double endX = x[nPoints-1];
  double endY = y[nPoints-1];
  
  a2 = endY-vtxY;
  if(TMath::Abs(vtxX-endX)>0.) a2 /= (endX-vtxX);
  else{
    cout << "cannot find slope for n = " << nPoints << " points " << endl;
    for(Int_t i = 0; i<nPoints; ++i) cout << x[i] << " " << y[i] << endl;
  }
  
  a1 = endY-a2*endX;
  
  //find the chi^2 value for the fit
  deviation = 0.;
  Double_t dev = 0.;
  for(Int_t j = 0; j<nPoints; j++){
    
    dev = (y[j] - (a1 + a2*x[j]));
    
    //     cout << y[j] << " " << a1+a2*x[j] << " " << chiSq << endl;
    
    deviation += dev*dev;
    
  }  
  
  return;
}

//----------------------------------------------------------------
//find the signal weighted rms of the hits in the track for a given
//plane
void CondensedNtpModuleAtm::FindMeanAndRMS(double *x, double *weight, 
                                           int nPoints,
                                           double &mean, double &rms)
{
  
  //   cout << "rms" << endl;
  
  double sumXXW = 0.;
  double sumXW = 0.;
  double sumW = 0.;
  
  for(Int_t i = 0; i<nPoints; ++i){
    sumXXW += x[i]*x[i]*weight[i];
    sumXW += x[i]*weight[i];
    sumW += weight[i];
    //     cout << sumXW << " " << sumXXW << " " << sumW << " " 
    //   << x[i] << " " << weight[i] << endl;
  }
  
  if(sumW > 0. && nPoints>1){
    mean = sumXW/sumW;
    if(sumXXW >=0.) rms =  TMath::Sqrt((sumXXW/sumW - mean*mean)
                                       *(nPoints*1./(1.*nPoints-1.)));
  }
  else{
    rms = 0.;
    mean = x[0];
  }
  return;
}

//--------------------------------------------------------------
//return the baseline in meters
Float_t CondensedNtpModuleAtm::CheapBaseline(float cosz)
{
  static float Rearth = 6378140.0; /* Equatorial radius of earth */
  static float nuHeight = 1500.0; /* Average Neutrino production height */
  static float altitudeSK = -210.0; /* Altitude of the detector - MINOS is 689ft below sea level*/
  
  static float x1s = (Rearth+nuHeight)*(Rearth+nuHeight);
  static float x2 = Rearth+altitudeSK;
  static float x2s = (Rearth+altitudeSK)*(Rearth+altitudeSK);
  
  return sqrt(x1s+x2s*(cosz*cosz-1.0))-x2*cosz;
}

---