MajorityCurvature.cxx

Go to the documentation of this file.

// $Id: MajorityCurvature.cxx,v 1.14 2009/12/08 18:09:26 ahimmel Exp $
//
// Class to calculate the Majority Curvature of a track.
//
// Justin Evans
// j.evans2@physics.ox.ac.uk
//
// $Log: MajorityCurvature.cxx,v $
// Revision 1.14  2009/12/08 18:09:26  ahimmel
//
//
// Removing dependence on NuBarPID now that it's maintainer has moved on.  Let me know ASAP if you see unexpected problems.
//
// Revision 1.13  2009/09/20 15:20:40  ahimmel
//
//
// Now use a NuGeneral object directly rather than  NuLibrary singleton.  This should avoid the chicken-and-egg problem.  I have tested that it runs normally, however I cannot test the file-handling behavior until this change is picked up by the development build.
//
// Revision 1.12  2009/09/20 12:24:54  nickd
// Reverting to version 1.10 to fix chicken-and-egg bug
//
// MajorityCurvature.cxx instanced NuLibrary in it's constructor - NuLibrary also did the same.
// Rather than fixing this myself now, I am reverting it so that the code 'works'
//
// Revision 1.10  2007/10/17 19:17:09  evans
// Removing an unused variable.
//
// Revision 1.9  2007/10/17 19:12:18  evans
// Discovered Google Perftools and made my majority curvature code at
// least a factor of 5 faster. No change in the user interface or
// results.
//
// Revision 1.8  2007/08/14 12:04:17  evans
// Changing the calulation of the smoothMajC variable to be the most up to
// date.
//
// Revision 1.7  2007/08/13 12:47:31  hartnell
//
// Set the gDirectory back to what it was before opening the PID file.
//
// Revision 1.6  2007/07/26 10:32:14  evans
// Adding an implementation of code to make a beam matrix for extrapolation.
//
// NuFluxHelper makes the beam matrix. It's configurable to make the matrix
// for either numus or numubars, and even to make a non-square beam matrix if
// the mood takes you.
//
// NuFluxChain is a wrapper around a fluka file to give the interface required
// by NuFluxHelper.
//
// macros/MakeMMFluxHelpers.C is an example macro of how to make the beam
// matrix.
//
// Revision 1.5  2007/04/13 14:50:35  evans
// Adding a second version of the track jitter using the square of the
// differences between the curvatures of the successive track segments.
//
// Revision 1.4  2007/03/27 13:31:20  evans
// New code to calculate my charge PID (requires pdfs in /data).
// Also includes code to calculate track jitter given an NtpSRTrack object.
//
// Revision 1.3  2007/03/24 20:28:59  evans
// Slight change to the track jitter variable to make sure it's well-defined
// for all tracks.
//
// Revision 1.2  2007/03/09 18:02:11  evans
// Adding extra functionality to MajorityCurvature to allow me to experiment
// with it's capabilities. This requires the addition of the (very
// lightweight) class MajCInfo.
//
// None of the changes affect the function
// MajorityCurvature::Curvature(NtpSRTrack), so anyone using this function
// will be unaffected.
//
// Revision 1.1  2007/02/05 14:45:51  evans
// Adding the code to calculate the majority curvature for a track.
//
// To calculate the majority curvature for an NtpSRTrack object, pass it to the function MajorityCurvature::Curvature(const NtpSRTrack* track) const.
// The return value is the majority curvature of the track.
//
// The class SRMom is a wrapper to the MomNavigator that has been submitted as MajorityCurvature uses some of its methods.
//
//
//
#include <cmath>

#include "TCanvas.h"
#include "TH1D.h"
#include "TF1.h"
#include "TFile.h"
#include "TGraphErrors.h"
#include "TMath.h"

#include <map>
#include "NtupleUtils/NuGeneral.h"
#include "CandNtupleSR/NtpSRTrack.h"
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"

#include "NtupleUtils/SRMom.h"
#include "NtupleUtils/MajorityCurvature.h"

ClassImp(MajorityCurvature)

CVSID("$Id: MajorityCurvature.cxx,v 1.14 2009/12/08 18:09:26 ahimmel Exp $");

TH1F* MajorityCurvature::hvariances = new TH1F("hvariances","",200,-5,5);

//____________________________________________________________________72
MajorityCurvature::MajorityCurvature()
{

//     topDir="MCReweight/data";
//     std::string base="";
//     base=getenv("SRT_PRIVATE_CONTEXT");
//     if(base=="" || base==".") base=getenv("SRT_PUBLIC_CONTEXT");
//     if(base=="") {
//       cout<<"No SRT_PUBLIC_CONTEXT set"<<std::endl;
//       assert(false);
//     }
//     topDir = base+ "/" + topDir;


  TDirectory* tmpd = gDirectory;
  MSG("MajorityCurvature",Msg::kDebug)
    <<"TDirectory before opening file is:"<<endl;
  //tmpd->Print();
    
  const NuGeneral& general=NuGeneral();
  TString filename=general.GetReleaseDirToUse("NtupleUtils");
  filename += "/NtupleUtils/data/pid_evans_cedar_daikon_LE185.root";
  fPDFFile = new
    TFile(filename, "READ");

  fhMajCBa = dynamic_cast<TH1D*> (fPDFFile->Get("hMajCBa"));
  fhMajCBa->SetDirectory(0);
  if (!fhMajCBa){
    cout << "No MajCBa PDF. Expect a seg fault at any moment!" << endl;
  }
  fhMajCNC = dynamic_cast<TH1D*> (fPDFFile->Get("hMajCNC"));
  fhMajCNC->SetDirectory(0);
  if (!fhMajCNC){
    cout << "No MajCNC PDF. Expect a seg fault at any moment!" << endl;
  }
  fhMajCNu = dynamic_cast<TH1D*> (fPDFFile->Get("hMajCNu"));
  fhMajCNu->SetDirectory(0);
  if (!fhMajCNu){
    cout << "No MajCNu PDF. Expect a seg fault at any moment!" << endl;
  }
  fhQOverPBa = dynamic_cast<TH1D*> (fPDFFile->Get("hQOverPBa"));
  fhQOverPBa->SetDirectory(0);
  if (!fhQOverPBa){
    cout << "No QOverPBa PDF. Expect a seg fault at any moment!" << endl;
  }
  fhQOverPNC = dynamic_cast<TH1D*> (fPDFFile->Get("hQOverPNC"));
  fhQOverPNC->SetDirectory(0);
  if (!fhQOverPNC){
    cout << "No QOverPNC PDF. Expect a seg fault at any moment!" << endl;
  }
  fhQOverPNu = dynamic_cast<TH1D*> (fPDFFile->Get("hQOverPNu"));
  fhQOverPNu->SetDirectory(0);
  if (!fhQOverPNu){
    cout << "No QOverPNu PDF. Expect a seg fault at any moment!" << endl;
  }
  fPDFFile->Close();

  MSG("MajorityCurvature",Msg::kDebug)
    <<"TDirectory after using file is:"<<endl;
  //gDirectory->Print();
  MSG("MajorityCurvature",Msg::kDebug)
    <<"After reseting gDirectory its location is:"<<endl;
  gDirectory=tmpd;
  //gDirectory->Print();
  
  for (Int_t row=0; row<500; ++row){
    for (Int_t column=0; column<500; ++column){
      fitLengths[row][column] = 0;
    }
  }
   
  //Iteration 1
  fitLengths[0][0] = 0;
  fitLengths[1][0] = 0;
  fitLengths[2][0] = 0;

  for(Int_t tl=3; tl<8; ++tl){
    Int_t i=0;
    for(Int_t fl=tl; fl>1; --fl, ++i){
      fitLengths[tl][i] = fl;
      if (2==fl){fitLengths[tl][i] = 0;}
    }
  }

  for(Int_t tl=8; tl<11; ++tl){
    Int_t i=0;
    for(Int_t fl = 7; fl>3; --fl,++i){
      fitLengths[tl][i] = fl;
    }
    for(Int_t fl = 8; fl<=tl; ++fl,++i){
      fitLengths[tl][i] = fl;
    }
    fitLengths[tl][i] = 3;
    ++i;
    fitLengths[tl][i] = 0;
  }

  for (Int_t tl=11; tl<500; ++tl){
    fitLengths[tl][0] = 7;
    fitLengths[tl][1] = 8;
    fitLengths[tl][2] = 6;
    fitLengths[tl][3] = 5;
    fitLengths[tl][4] = 4;
    fitLengths[tl][5] = 0;
  }
  
  /*  
  //Iteration 2
  fitLengths[0][0] = 0;
  fitLengths[1][0] = 0;
  fitLengths[2][0] = 0;

  for(Int_t tl=3; tl<8; ++tl){
    Int_t i=0;
    for(Int_t fl=tl; fl>1; --fl, ++i){
      fitLengths[tl][i] = fl;
      if (2==fl){fitLengths[tl][i] = 0;}
    }
  }

  for (Int_t tl = 8; tl < 500; ++tl){
    Int_t i=0;
    for (Int_t fl = 7; fl <= tl; ++fl, ++i){
      fitLengths[tl][i] = fl;
    }
    ++i;
    for (Int_t fl=6; fl > 2; --fl, ++i){
      fitLengths[tl][i] = fl;
    }
    ++i;
    fitLengths[tl][i] = 0;
  }
  */
}

//____________________________________________________________________72
MajorityCurvature::~MajorityCurvature()
{
  if (fhMajCBa) {delete fhMajCBa; fhMajCBa = 0;}
  if (fhMajCNu) {delete fhMajCNu; fhMajCNu = 0;}
  if (fhMajCNC) {delete fhMajCNC; fhMajCNC = 0;}
  if (fhQOverPBa) {delete fhQOverPBa; fhQOverPBa = 0;}
  if (fhQOverPNu) {delete fhQOverPNu; fhQOverPNu = 0;}
  if (fhQOverPNC) {delete fhQOverPNC; fhQOverPNC = 0;}
  fPDFFile->Close();
  if (fPDFFile) {delete fPDFFile; fPDFFile = 0;}
}

//____________________________________________________________________72
const Float_t MajorityCurvature::CalcCombCurv(const MomNavigator* mom,
                                              const Int_t event,
                                              const Int_t track,
                                              const Int_t fitLength) const
{
  //   if (snarl > 11){return JobCResult::kFailed;}
  SRMom srMom(mom);
  TH1F hCurv("hCurv","",100,-5,5);

  const multimap<Float_t, Float_t> mu
    = srMom.ArrayUAlongTrack(event,track);
  const Int_t sizeu = mu.size();
  Float_t* xerr = new Float_t[fitLength];
  Float_t* zerr = new Float_t[fitLength];
  Bool_t smallu = false;

  if (sizeu < fitLength){smallu = true;}
  else{
    multimap<Float_t,Float_t>::const_iterator itu = mu.begin();
    for (Int_t k=0; k<(sizeu-fitLength); ++k){
      Float_t* z = new Float_t[fitLength];
      Float_t* x = new Float_t[fitLength];
      Float_t* absx = new Float_t[fitLength];
      Bool_t changesign = false;
      Bool_t positive = false;
      if (itu->second > 0){positive = true;}
      for (Int_t l=0; l<fitLength; ++l)
        {
          Bool_t thispos = false;
          z[l] = itu->first;
          zerr[l] = 0.0594;
          x[l] = itu->second;
          xerr[l] = 0.041;
          absx[l] = fabs(itu->second);
          if (x[l] > 0){thispos = true;}
          if (thispos != positive){changesign = true;}
          ++itu;
        }
      if (!changesign){
        TGraphErrors g(fitLength,z,absx,zerr,xerr);
        //      TCanvas cd("cd","",0,0,500,500);
        g.Draw("A*");
        g.Fit("pol2","Q");
        //   cd.Print(Form("~/public_html/AASnarl%dEvent%dSegment%d.gif",snarl,event,k));
        Float_t curviness = 2*(g.GetFunction("pol2")->GetParameter(2));
        //       for (Double_t i = -99999999; i < 99999999; i+=0.1){
        //      Float_t j = 1; ++j;
        //       }
        hCurv.Fill(curviness);
      }
      for (Int_t l=0; l<fitLength; ++l)
        {
          --itu;
        }
      ++itu;
      if(z){delete z; z=0;}
      if(x){delete x; x=0;}
      if(absx){delete absx; absx=0;}
    }
  }
  
  const multimap<Float_t, Float_t> mv
    = srMom.ArrayVAlongTrack(event,track);
  const Int_t sizev = mv.size();
  if (sizev < fitLength){
    if (smallu){return 0.0;}
  }
  else{
    multimap<Float_t,Float_t>::const_iterator itv = mv.begin();
    for (Int_t k=0; k<(sizev-fitLength); ++k){
      Float_t* z = new Float_t[fitLength];
      Float_t* x = new Float_t[fitLength];
      Float_t* absx = new Float_t[fitLength];
      Bool_t changesign = false;
      Bool_t positive = false;
      if (itv->second > 0){positive = true;}
      for (Int_t l=0; l<fitLength; ++l)
        {
          Bool_t thispos = false;
          z[l] = itv->first;
          zerr[l] = 0.0594;
          x[l] = itv->second;
          xerr[l] = 0.041;
          absx[l] = fabs(itv->second);
          if (x[l] > 0){thispos = true;}
          if (thispos != positive){changesign = true;}
          ++itv;
        }
      if (!changesign){
        TGraphErrors g(fitLength,z,absx,zerr,xerr);
        //      TCanvas cd("cd","",0,0,500,500);
        g.Draw("A*");
        g.Fit("pol2","Q");
        //   cd.Print(Form("~/public_html/AASnarl%dEvent%dSegment%d.gif",snarl,event,k));
        Float_t curviness = 2*(g.GetFunction("pol2")->GetParameter(2));
        //       for (Double_t i = -99999999; i < 99999999; i+=0.1){
        //      Float_t j = 1; ++j;
        //       }
        hCurv.Fill(curviness);
      }
      for (Int_t l=0; l<fitLength; ++l)
        {
          --itv;
        }
      ++itv;
      if(z){delete z; z=0;}
      if(x){delete x; x=0;}
      if(absx){delete absx; absx=0;}
    }
  }
  TCanvas cCurv("cCurv","",0,0,500,500);
  hCurv.Draw();
  //   width = hCurv.GetRMS();
  //   cCurv.Print(Form("~/public_html/ASnarl%dEvent%d.gif",snarl,event));
  Float_t neg = hCurv.Integral(0,50);
  Float_t pos = hCurv.Integral(51,101);
  if(xerr){delete xerr; xerr = 0;}
  if(zerr){delete zerr; zerr = 0;}
  if (pos > neg) {if (neg) {return pos/neg;} else {return 10.0;}}
  if (neg > pos) {if (pos) {return -1.0*neg/pos;} else {return -10.0;}}
  if (neg == pos) {return 0.0;}
  return 0.0;
}

//____________________________________________________________________72
const Float_t MajorityCurvature::CalcCurv(const multimap<Float_t,
                                          Float_t> md,
                                          const Int_t fitLength) const
{
  //Set up some arrays.
  const Int_t size = md.size();
  Float_t* xerr = new Float_t[fitLength];
  Float_t* zerr = new Float_t[fitLength];
  if (size<fitLength){return 0.0;}

  multimap<Float_t, Float_t> mfits; //<abs(parabolic fit),(parabolic fit)>
  multimap<Float_t,Float_t>::const_iterator it = md.begin();
  //k loops over the fit segments.
  for (Int_t k=0; k<(size-fitLength); ++k){
    Float_t* z = new Float_t[fitLength];
    Float_t* x = new Float_t[fitLength];
    Float_t* absx = new Float_t[fitLength];
    Bool_t changesign = false;
    Bool_t positive = false;
    if (it->second > 0){positive = true;}
    //l loops over the strips in each fit segment
    //The variable 'positive' says whether the first strip in a
    //segment has a positive co-ordinate.
    for (Int_t l=0; l<fitLength; ++l)
      {
        Bool_t thispos = false;
        z[l] = it->first;
        zerr[l] = 0.0594;
        x[l] = it->second;
        xerr[l] = 0.041;
        absx[l] = fabs(it->second);
        if (x[l] > 0){thispos = true;}
        if (thispos != positive){changesign = true;}
        ++it;
      }
    if (!changesign){
      TGraphErrors g(fitLength,z,absx,zerr,xerr);
      //      TCanvas cd("cd","",0,0,500,500);
      g.Draw("A*");
      g.Fit("pol2","Q");
      Float_t curviness = 2*(g.GetFunction("pol2")->GetParameter(2));
      //      hCurv.Fill(curviness);
      mfits.insert(pair<Float_t,Float_t>(fabs(curviness),
                                         curviness));
    }
    for (Int_t l=0; l<fitLength; ++l)
      {
        --it;
      }
    ++it;
    if(z){delete z; z=0;}
    if(x){delete x; x=0;}
    if(absx){delete absx; absx=0;}
  }

  multimap<Float_t,Float_t>::const_iterator itfits = mfits.end();
  Int_t numRemove = 0;
  if (((Int_t) mfits.size()) >= numRemove){
    for (Int_t rem = 0; rem < numRemove; ++rem){
      itfits--;
      mfits.erase(itfits->first);
      itfits = mfits.end();
    }
  }
  TH1F hCurv("hCurv","",100,-5,5);
  for (itfits = mfits.begin();
       itfits != mfits.end();
       ++itfits){
    hCurv.Fill(itfits->second);
  }

  TCanvas cCurv("cCurv","",0,0,500,500);
  hCurv.Draw();
//   width = hCurv.GetRMS();
//   cCurv.Print(Form("~/public_html/ASnarl%dEvent%d.gif",snarl,event));
  Float_t neg = hCurv.Integral(0,50);
  Float_t pos = hCurv.Integral(51,101);
  hvariances->Fill(hCurv.GetRMS());
  //  flastVariance = hCurv.GetRMS();
  if (pos > neg) {if (neg) {return pos/neg;} else {return 10.0;}}
  if (neg > pos) {if (pos) {return -1.0*neg/pos;} else {return -10.0;}}
  if (neg == pos) {return 0.0;}
  else {return 0.0;}
}

//____________________________________________________________________72
const MajCInfo MajorityCurvature::CalcCurvObject(const multimap<Float_t,
                                                 Float_t> md,
                                                 const Int_t fitLength)
  const
{
  //Object to return:
  MajCInfo majCInfo;

  //Set up some arrays.
  const Int_t size = md.size();
//   Float_t* xerr = new Float_t[fitLength];
//   Float_t* zerr = new Float_t[fitLength];
  if (size<fitLength){return majCInfo;}

  multimap<Float_t, Float_t> mfits; //<abs(parabolic fit),(parabolic fit)>
  multimap<Float_t,Float_t>::const_iterator it = md.begin();
  Double_t totJitter = 0.0;
  Double_t totSqJitter = 0.0;
  Double_t prevCurvature = 0.0;
  Int_t numSegments = 0;
  /*  Int_t endChop = 3;*///Don't use the final n segments in jitter
  //k loops over the fit segments.
  TGraphErrors gQuick(fitLength);
  for (Int_t k=0; k<(size-fitLength); ++k){
//     Float_t* z = new Float_t[fitLength];
//     Float_t* x = new Float_t[fitLength];
//     Float_t* absx = new Float_t[fitLength];
    Bool_t changesign = false;
    Bool_t positive = false;
    if (it->second > 0){positive = true;}
    //l loops over the strips in each fit segment
    //The variable 'positive' says whether the first strip in a
    //segment has a positive co-ordinate.
    for (Int_t l=0; l<fitLength; ++l)
      {
        Bool_t thispos = false;
        gQuick.SetPoint(l,it->first,fabs(it->second));
        gQuick.SetPointError(l,0.0594,0.041);
//      z[l] = it->first;
//      zerr[l] = 0.0594;
//      x[l] = it->second;
//      xerr[l] = 0.041;
//      absx[l] = fabs(it->second);
        if (it->second > 0){thispos = true;}
        if (thispos != positive){changesign = true;}
        ++it;
      }
    if (!changesign){
//       gQuick.Draw("A*");
      gQuick.Fit("pol2","Q");
      Float_t curviness = 2*(gQuick.GetFunction("pol2")->GetParameter(2));
      if ((numSegments/*-endChop*/)/* && (k<(size-fitLength-endChop))*/){
        totJitter += fabs(curviness - prevCurvature);
        totSqJitter +=
          (curviness - prevCurvature)*(curviness - prevCurvature);
      }
      ++numSegments;
      prevCurvature = curviness;
      //      hCurv.Fill(curviness);
      mfits.insert(pair<Float_t,Float_t>(fabs(curviness),
                                         curviness));
    }
    for (Int_t l=0; l<fitLength; ++l)
      {
        --it;
      }
    ++it;
//     if(z){delete z; z=0;}
//     if(x){delete x; x=0;}
//     if(absx){delete absx; absx=0;}

    majCInfo.jitter = -1.0;
    majCInfo.sqJitter = -1.0;
//     if (numSegments>1){
//       majCInfo.jitter = totJitter/((Double_t) (numSegments-1));
//      majCInfo.sqJitter = totSqJitter/((Double_t) (numSegments-1));
//     }
    if (numSegments/*-endChop*/){
      majCInfo.jitter = totJitter/((Double_t) (numSegments/*-endChop*/));
      majCInfo.sqJitter = totSqJitter/((Double_t) (numSegments/*-endChop*/));
    }
    gQuick.Clear();
  }

  multimap<Float_t,Float_t>::const_iterator itfits = mfits.end();
  Int_t numRemove = 0;
  if (((Int_t) mfits.size()) >= numRemove){
    for (Int_t rem = 0; rem < numRemove; ++rem){
      itfits--;
      mfits.erase(itfits->first);
      itfits = mfits.end();
    }
  }
  TH1F hCurv("hCurv","",100,-5,5);
  Float_t smallestCurv=100.0;
  Float_t largestCurv = -100.0;
  for (itfits = mfits.begin();
       itfits != mfits.end();
       ++itfits){
    hCurv.Fill(itfits->second);
    if (itfits->second<smallestCurv){smallestCurv=itfits->second;}
    if (itfits->second>largestCurv){largestCurv=itfits->second;}
  }

//   TCanvas cCurv("cCurv","",0,0,500,500);
//   hCurv.Draw();
  majCInfo.rms = hCurv.GetRMS();
  majCInfo.totWidth = largestCurv-smallestCurv;
//   cCurv.Print(Form("~/public_html/ASnarl%dEvent%d.gif",snarl,event));
  Float_t neg = hCurv.Integral(0,50);
  Float_t pos = hCurv.Integral(51,101);
  hvariances->Fill(hCurv.GetRMS());
  //  flastVariance = hCurv.GetRMS();
  majCInfo.simpleMajC = pos-neg;
  if (neg){majCInfo.majCRatio = pos/neg;}
  else{majCInfo.majCRatio = 10.0;}

  if (!pos){
    majCInfo.smoothMajC = -10.0;
  }
  else {
    if (!neg) {
      majCInfo.smoothMajC = 10.0;
    }
    else{
      majCInfo.smoothMajC = TMath::Log(pos/neg);
    }
  }

  if (pos > neg) {
    if (neg) {majCInfo.majC = pos/neg; return majCInfo;}
    else {majCInfo.majC = 10.0; return majCInfo;}
  }
  if (neg > pos) {
    if (pos) {majCInfo.majC = -1.0*neg/pos; return majCInfo;}
    else {majCInfo.majC = -10.0; return majCInfo;}
  }
  if (neg == pos) {majCInfo.majC = 0.0; return majCInfo;}
  else {majCInfo.majC = 0.0; return majCInfo;}
  return majCInfo;
}

//____________________________________________________________________72
const MajCInfo MajorityCurvature::CalcCurvObjectFast(const multimap<Float_t,
                                                     Float_t> md,
                                                     const Int_t fitLength)
  const
{
  //Object to return:
  MajCInfo majCInfo;

  //Set up some arrays.
  const Int_t size = md.size();
//   Float_t* xerr = new Float_t[fitLength];
//   Float_t* zerr = new Float_t[fitLength];
  if (size<fitLength){return majCInfo;}

  multimap<Float_t,Float_t>::const_iterator it = md.begin();
  Int_t numPositiveSegments = 0;
  Int_t numNegativeSegments = 0;
  //k loops over the fit segments.
  for (Int_t k=0; k<(size-fitLength); ++k){
    Bool_t changesign = false;
    Bool_t positive = false;
    if (it->second > 0){positive = true;}
    //l loops over the strips in each fit segment
    //The variable 'positive' says whether the first strip in a
    //segment has a positive co-ordinate.
    ++it;
    Int_t numPositiveIncrements = 0;
    Int_t numNegativeIncrements = 0;
    for (Int_t l=1; l<fitLength-1; ++l)
      {
        Bool_t thispos = false;
        --it;
        Double_t prevZ = it->first;
        Double_t prevX = it->second;
        ++it;
        Double_t thisZ = it->first;
        Double_t thisX = it->second;
        ++it;
        Double_t nextZ = it->first;
        Double_t nextX = it->second;
        if ((nextZ-thisZ) && (thisZ-prevZ)){
          Double_t forwardGrad = (nextX-thisX)/(nextZ-thisZ);
          Double_t backGrad = (thisX-prevX)/(thisZ-prevZ);
          if (positive){
            if (forwardGrad>backGrad){
              ++numPositiveIncrements;
            }
            if (backGrad>forwardGrad){
              ++numNegativeIncrements;
            }
          }
          else{
            if (forwardGrad>backGrad){
              ++numNegativeIncrements;
            }
            if (backGrad>forwardGrad){
              ++numPositiveIncrements;
            }
          }
        }
        --it;
        if (it->second > 0){thispos = true;}
        if (thispos != positive){changesign = true;}
        ++it;
      }
    if (!changesign){
      if (numPositiveIncrements>numNegativeIncrements){
        ++numPositiveSegments;
      }
      if (numNegativeIncrements>numPositiveIncrements){
        ++numNegativeSegments;
      }
    }
    for (Int_t l=1; l<fitLength-1; ++l)
      {
        --it;
      }
    ++it;

    majCInfo.jitter = -1.0;
    majCInfo.sqJitter = -1.0;
  }

  majCInfo.simpleMajC = numPositiveSegments-numNegativeSegments;
  if (numNegativeSegments){majCInfo.majCRatio = numPositiveSegments/numNegativeSegments;}
  else{majCInfo.majCRatio = 10.0;}

  if (!numPositiveSegments){
    majCInfo.smoothMajC = -10.0;
  }
  else {
    if (!numNegativeSegments) {
      majCInfo.smoothMajC = 10.0;
    }
    else{
      majCInfo.smoothMajC = TMath::Log(numPositiveSegments/numNegativeSegments);
    }
  }

  if (numPositiveSegments > numNegativeSegments) {
    if (numNegativeSegments) {majCInfo.majC = numPositiveSegments/numNegativeSegments; return majCInfo;}
    else {majCInfo.majC = 10.0; return majCInfo;}
  }
  if (numNegativeSegments > numPositiveSegments) {
    if (numPositiveSegments) {majCInfo.majC = -1.0*numNegativeSegments/numPositiveSegments; return majCInfo;}
    else {majCInfo.majC = -10.0; return majCInfo;}
  }
  if (numNegativeSegments == numPositiveSegments) {majCInfo.majC = 0.0; return majCInfo;}
  else {majCInfo.majC = 0.0; return majCInfo;}
  return majCInfo;
}

//____________________________________________________________________72
const Float_t MajorityCurvature::CalcXCurv(const MomNavigator* mom,
                                           const Int_t event,
                                           const Int_t track,
                                           const Int_t fitLength) const
{
  SRMom srMom(mom);
  const multimap<Float_t, Float_t> md
    = srMom.ArrayXAlongTrack(event,track);
  return CalcCurv(md, fitLength);
}

//____________________________________________________________________72
const MajCInfo MajorityCurvature::
CalcXCurvObject(const MomNavigator* mom,
                const Int_t event,
                const Int_t track,
                const Int_t fitLength) const
{
  SRMom srMom(mom);
  const multimap<Float_t, Float_t> md
    = srMom.ArrayXAlongTrack(event,track);
  return CalcCurvObject(md, fitLength);
}

//____________________________________________________________________72
const MajCInfo MajorityCurvature::
CalcXCurvObject(const NtpSRTrack* track,
                const Int_t fitLength) const
{
  SRMom srMom;
  const multimap<Float_t, Float_t> md
    = srMom.ArrayXAlongTrack(track);
  return CalcCurvObject(md, fitLength);
}

//____________________________________________________________________72
const Float_t MajorityCurvature::CalcXCurv(const NtpSRTrack* track,
                                           const Int_t fitLength) const
{
  SRMom srMom;
  const multimap<Float_t, Float_t> md
    = srMom.ArrayXAlongTrack(track);
  return CalcCurv(md, fitLength);
}

//____________________________________________________________________72
const Float_t MajorityCurvature::CalcYCurv(const MomNavigator* mom,
                                           const Int_t event,
                                           const Int_t track,
                                           const Int_t fitLength) const
{
  SRMom srMom(mom);
  const multimap<Float_t, Float_t> md
    = srMom.ArrayYAlongTrack(event,track);
  return CalcCurv(md, fitLength);
}

//____________________________________________________________________72
const MajCInfo MajorityCurvature::
CalcYCurvObject(const MomNavigator* mom,
                const Int_t event,
                const Int_t track,
                const Int_t fitLength) const
{
  SRMom srMom(mom);
  const multimap<Float_t, Float_t> md
    = srMom.ArrayYAlongTrack(event,track);
  return CalcCurvObject(md, fitLength);
}

//____________________________________________________________________72
const MajCInfo MajorityCurvature::
CalcYCurvObject(const NtpSRTrack* track,
                const Int_t fitLength) const
{
  SRMom srMom;
  const multimap<Float_t, Float_t> md
    = srMom.ArrayYAlongTrack(track);
  return CalcCurvObject(md, fitLength);
}

//____________________________________________________________________72
const Float_t MajorityCurvature::CalcYCurv(const NtpSRTrack* track,
                                           const Int_t fitLength) const
{
  SRMom srMom;
  const multimap<Float_t, Float_t> md
    = srMom.ArrayYAlongTrack(track);
  return CalcCurv(md, fitLength);
}

//____________________________________________________________________72
const Float_t MajorityCurvature::Curvature(const NtpSRTrack* track) const
{
  Float_t xCurv = XCurvature(track);
  Float_t yCurv = YCurvature(track);
  if (fabs(yCurv) > fabs(xCurv)){return yCurv;}
  else {return xCurv;}
}

//____________________________________________________________________72
const MajCInfo MajorityCurvature::
CurvatureImproved(const NtpSRTrack* track) const
{
  MajCInfo xCurv = XCurvatureObject(track);
  MajCInfo yCurv = YCurvatureObject(track);
  MajCInfo* bestCurv;
  if (fabs(yCurv.majC) > fabs(xCurv.majC)){bestCurv = &yCurv;}
  else {bestCurv = &xCurv;}
  SRMom srMom;
  Double_t qp = srMom.TrackVtxQOvP(track);
  bestCurv->jPID = this->JPID(qp,bestCurv->majC);
  return *bestCurv;
}

//____________________________________________________________________72
const MajCInfo MajorityCurvature::
CurvatureImproved(const MomNavigator* mom,
                  const Int_t event,
                  const Int_t track) const
{
  MajCInfo xCurv = XCurvatureObject(mom,event,track);
  MajCInfo yCurv = YCurvatureObject(mom,event,track);
  MajCInfo* bestCurv;
  if (fabs(yCurv.majC) > fabs(xCurv.majC)){bestCurv = &yCurv;}
  else {bestCurv = &xCurv;}
  SRMom srMom(mom);
  Double_t qp = srMom.TrackVtxQOvP(event,track);
  bestCurv->jPID = this->JPID(qp,bestCurv->majC);
  return *bestCurv;
}

//____________________________________________________________________72
const Float_t MajorityCurvature::Curvature(const MomNavigator* mom,
                                           const Int_t event,
                                           const Int_t track) const
{
  Float_t xCurv = XCurvature(mom,event,track);
  Float_t yCurv = YCurvature(mom,event,track);
  if (fabs(yCurv) > fabs(xCurv)){return yCurv;}
  else {return xCurv;}
}

//____________________________________________________________________72
const Float_t MajorityCurvature::CurvatureComb(const MomNavigator* mom,
                                               const Int_t event,
                                               const Int_t track) const
{
  SRMom srMom(mom);
  const Int_t trkLthPln = srMom.TrackLengthPlanes(event,track);
  Int_t fitLength = FitLength(trkLthPln);
  if (!fitLength) {return 0.0;}
  Float_t xCurv = CalcCombCurv(mom,event,track,fitLength);
  if (xCurv) {return xCurv;}
  while (!xCurv){
    fitLength = NextFitLength(trkLthPln,fitLength);
    if (!fitLength){return 0.0;}
    xCurv = CalcCombCurv(mom,event,track,fitLength);
  }
  return xCurv;
}

//____________________________________________________________________72
const Int_t MajorityCurvature::FitLength(const Int_t trkLthPln) const
{
  return fitLengths[trkLthPln][0];
}

//____________________________________________________________________72
const Double_t MajorityCurvature::JPID(const Double_t qOverP,
                                       const Double_t majC) const
{
  Double_t pBackMajC =
    fhMajCNu->GetBinContent(fhMajCNu->GetXaxis()->FindBin(majC));
  pBackMajC +=
    fhMajCNC->GetBinContent(fhMajCNC->GetXaxis()->FindBin(majC));
  Double_t pSigMajC =
    fhMajCBa->GetBinContent(fhMajCBa->GetXaxis()->FindBin(majC));

  Double_t pBackqOvp =
    fhQOverPNu->GetBinContent(fhQOverPNu->GetXaxis()->FindBin(qOverP));
  pBackqOvp +=
    fhQOverPNC->GetBinContent(fhQOverPNC->GetXaxis()->FindBin(qOverP));
  Double_t pSigqOvp =
    fhQOverPBa->GetBinContent(fhQOverPBa->GetXaxis()->FindBin(qOverP));

  Float_t numerator =
    pBackMajC*pBackqOvp;
  if (numerator<0.0001){numerator = 0.0001;}
  Float_t denominator =
    pSigMajC*pSigqOvp;
  if (denominator<0.0001){denominator = 0.0001;}
  
  Double_t jPID = log(numerator/denominator);
  return jPID;
}

//____________________________________________________________________72
const Int_t MajorityCurvature::NextFitLength(const Int_t trkLthPln,
                                             const Int_t lastFitLength)
  const
{
  Int_t i=0;
  Int_t fitLength = 0;
  for (i=0; fitLength != lastFitLength; ++i){
    fitLength = fitLengths[trkLthPln][i];
    if (!fitLength){return fitLength;}
  }
  return fitLengths[trkLthPln][i];
}

//____________________________________________________________________72
const Float_t MajorityCurvature::XCurvature(const MomNavigator* mom,
                                            const Int_t event,
                                            const Int_t track) const
{
  SRMom srMom(mom);
  const Int_t trkLthPln = srMom.TrackLengthPlanes(event,track);
  Int_t fitLength = FitLength(trkLthPln);
  if (!fitLength) {return 0.0;}
  Float_t xCurv = CalcXCurv(mom,event,track,fitLength);
  if (xCurv) {return xCurv;}
  while (!xCurv){
    fitLength = NextFitLength(trkLthPln,fitLength);
    if (!fitLength){return 0.0;}
    xCurv = CalcXCurv(mom,event,track,fitLength);
  }
  return xCurv;
}

//____________________________________________________________________72
const MajCInfo MajorityCurvature::
XCurvatureObject(const MomNavigator* mom,
                 const Int_t event,
                 const Int_t track) const
{
  SRMom srMom(mom);
  const Int_t trkLthPln = srMom.TrackLengthPlanes(event,track);
  Int_t fitLength = FitLength(trkLthPln);
  if (!fitLength) {MajCInfo empty; return empty;}
  MajCInfo xCurv = CalcXCurvObject(mom,event,track,fitLength);
  if (xCurv.majC) {return xCurv;}
  while (!xCurv.majC){
    fitLength = NextFitLength(trkLthPln,fitLength);
    if (!fitLength){MajCInfo empty; return empty;}
    xCurv = CalcXCurvObject(mom,event,track,fitLength);
  }
  return xCurv;
}

//____________________________________________________________________72
const MajCInfo MajorityCurvature::
XCurvatureObject(const NtpSRTrack* track) const
{
  SRMom srMom;
  const Int_t trkLthPln = srMom.TrackLengthPlanes(track);
  Int_t fitLength = FitLength(trkLthPln);
  if (!fitLength) {MajCInfo empty; return empty;}
  MajCInfo xCurv = CalcXCurvObject(track,fitLength);
  if (xCurv.majC) {return xCurv;}
  while (!xCurv.majC){
    fitLength = NextFitLength(trkLthPln,fitLength);
    if (!fitLength){MajCInfo empty; return empty;}
    xCurv = CalcXCurvObject(track,fitLength);
  }
  return xCurv;
}

//____________________________________________________________________72
const Float_t MajorityCurvature::XCurvature(const NtpSRTrack* track)
  const
{
  SRMom srMom;
  const Int_t trkLthPln = srMom.TrackLengthPlanes(track);
  Int_t fitLength = FitLength(trkLthPln);
  if (!fitLength) {return 0.0;}
  Float_t xCurv = CalcXCurv(track,fitLength);
  if (xCurv) {return xCurv;}
  while (!xCurv){
    fitLength = NextFitLength(trkLthPln,fitLength);
    if (!fitLength){return 0.0;}
    xCurv = CalcXCurv(track,fitLength);
  }
  return xCurv;
}

//____________________________________________________________________72
const Float_t MajorityCurvature::YCurvature(const MomNavigator* mom,
                                            const Int_t event,
                                            const Int_t track) const
{
  SRMom srMom(mom);
  const Int_t trkLthPln = srMom.TrackLengthPlanes(event,track);
  Int_t fitLength = FitLength(trkLthPln);
  if (!fitLength) {return 0.0;}
  Float_t yCurv = CalcYCurv(mom,event,track,fitLength);
  if (yCurv) {return yCurv;}
  while (!yCurv){
    fitLength = NextFitLength(trkLthPln,fitLength);
    if (!fitLength){return 0.0;}
    yCurv = CalcYCurv(mom,event,track,fitLength);
  }
  return yCurv;
}

//____________________________________________________________________72
const MajCInfo MajorityCurvature::
YCurvatureObject(const MomNavigator* mom,
                 const Int_t event,
                 const Int_t track) const
{
  SRMom srMom(mom);
  const Int_t trkLthPln = srMom.TrackLengthPlanes(event,track);
  Int_t fitLength = FitLength(trkLthPln);
  if (!fitLength) {MajCInfo empty; return empty;}
  MajCInfo yCurv = CalcYCurvObject(mom,event,track,fitLength);
  if (yCurv.majC) {return yCurv;}
  while (!yCurv.majC){
    fitLength = NextFitLength(trkLthPln,fitLength);
    if (!fitLength){MajCInfo empty; return empty;}
    yCurv = CalcYCurvObject(mom,event,track,fitLength);
  }
  return yCurv;
}

//____________________________________________________________________72
const MajCInfo MajorityCurvature::
YCurvatureObject(const NtpSRTrack* track) const
{
  SRMom srMom;
  const Int_t trkLthPln = srMom.TrackLengthPlanes(track);
  Int_t fitLength = FitLength(trkLthPln);
  if (!fitLength) {MajCInfo empty; return empty;}
  MajCInfo yCurv = CalcYCurvObject(track,fitLength);
  if (yCurv.majC) {return yCurv;}
  while (!yCurv.majC){
    fitLength = NextFitLength(trkLthPln,fitLength);
    if (!fitLength){MajCInfo empty; return empty;}
    yCurv = CalcYCurvObject(track,fitLength);
  }
  return yCurv;
}

//____________________________________________________________________72
const Float_t MajorityCurvature::YCurvature(const NtpSRTrack* track) const
{
  SRMom srMom;
  const Int_t trkLthPln = srMom.TrackLengthPlanes(track);
  Int_t fitLength = FitLength(trkLthPln);
  if (!fitLength) {return 0.0;}
  Float_t yCurv = CalcYCurv(track,fitLength);
  if (yCurv) {return yCurv;}
  while (!yCurv){
    fitLength = NextFitLength(trkLthPln,fitLength);
    if (!fitLength){return 0.0;}
    yCurv = CalcYCurv(track,fitLength);
  }
  return yCurv;
}


// //____________________________________________________________________72
// const MajCInfo MajorityCurvature::CalcCurvObject(const multimap<Float_t,
//                                               Float_t> md,
//                                               const Int_t fitLength)
//   const
// {
//   //Object to return:
//   MajCInfo majCInfo;

//   //Set up some arrays.
//   const Int_t size = md.size();
//   Float_t* xerr = new Float_t[fitLength];
//   Float_t* zerr = new Float_t[fitLength];
//   if (size<fitLength){return majCInfo;}

//   multimap<Float_t, Float_t> mfits; //<abs(parabolic fit),(parabolic fit)>
//   multimap<Float_t,Float_t>::const_iterator it = md.begin();
//   Double_t totJitter = 0.0;
//   Double_t totSqJitter = 0.0;
//   Double_t prevCurvature = 0.0;
//   Int_t numSegments = 0;
//   /*  Int_t endChop = 3;*///Don't use the final n segments in jitter
//   //k loops over the fit segments.
//   for (Int_t k=0; k<(size-fitLength); ++k){
//     Float_t* z = new Float_t[fitLength];
//     Float_t* x = new Float_t[fitLength];
//     Float_t* absx = new Float_t[fitLength];
//     Bool_t changesign = false;
//     Bool_t positive = false;
//     if (it->second > 0){positive = true;}
//     //l loops over the strips in each fit segment
//     //The variable 'positive' says whether the first strip in a
//     //segment has a positive co-ordinate.
//     for (Int_t l=0; l<fitLength; ++l)
//       {
//      Bool_t thispos = false;
//      z[l] = it->first;
//      zerr[l] = 0.0594;
//      x[l] = it->second;
//      xerr[l] = 0.041;
//      absx[l] = fabs(it->second);
//      if (x[l] > 0){thispos = true;}
//      if (thispos != positive){changesign = true;}
//      ++it;
//       }
//     if (!changesign){
//       TGraphErrors g(fitLength,z,absx,zerr,xerr);
// //       TCanvas cd("cd","",0,0,500,500);
//       g.Draw("A*");
//       g.Fit("pol2","Q");
//       Float_t curviness = 2*(g.GetFunction("pol2")->GetParameter(2));
//       if ((numSegments/*-endChop*/)/* && (k<(size-fitLength-endChop))*/){
//      totJitter += fabs(curviness - prevCurvature);
//      totSqJitter +=
//        (curviness - prevCurvature)*(curviness - prevCurvature);
//       }
//       ++numSegments;
//       prevCurvature = curviness;
//       //      hCurv.Fill(curviness);
//       mfits.insert(pair<Float_t,Float_t>(fabs(curviness),
//                                       curviness));
//     }
//     for (Int_t l=0; l<fitLength; ++l)
//       {
//      --it;
//       }
//     ++it;
//     if(z){delete z; z=0;}
//     if(x){delete x; x=0;}
//     if(absx){delete absx; absx=0;}

//     majCInfo.jitter = -1.0;
//     majCInfo.sqJitter = -1.0;
// //     if (numSegments>1){
// //       majCInfo.jitter = totJitter/((Double_t) (numSegments-1));
// //      majCInfo.sqJitter = totSqJitter/((Double_t) (numSegments-1));
// //     }
//     if (numSegments/*-endChop*/){
//       majCInfo.jitter = totJitter/((Double_t) (numSegments/*-endChop*/));
//       majCInfo.sqJitter = totSqJitter/((Double_t) (numSegments/*-endChop*/));
//     }
//   }

//   multimap<Float_t,Float_t>::const_iterator itfits = mfits.end();
//   Int_t numRemove = 0;
//   if (((Int_t) mfits.size()) >= numRemove){
//     for (Int_t rem = 0; rem < numRemove; ++rem){
//       itfits--;
//       mfits.erase(itfits->first);
//       itfits = mfits.end();
//     }
//   }
//   TH1F hCurv("hCurv","",100,-5,5);
//   Float_t smallestCurv=100.0;
//   Float_t largestCurv = -100.0;
//   for (itfits = mfits.begin();
//        itfits != mfits.end();
//        ++itfits){
//     hCurv.Fill(itfits->second);
//     if (itfits->second<smallestCurv){smallestCurv=itfits->second;}
//     if (itfits->second>largestCurv){largestCurv=itfits->second;}
//   }

//   TCanvas cCurv("cCurv","",0,0,500,500);
//   hCurv.Draw();
//   majCInfo.rms = hCurv.GetRMS();
//   majCInfo.totWidth = largestCurv-smallestCurv;
// //   cCurv.Print(Form("~/public_html/ASnarl%dEvent%d.gif",snarl,event));
//   Float_t neg = hCurv.Integral(0,50);
//   Float_t pos = hCurv.Integral(51,101);
//   hvariances->Fill(hCurv.GetRMS());
//   //  flastVariance = hCurv.GetRMS();
//   majCInfo.simpleMajC = pos-neg;
//   if (neg){majCInfo.majCRatio = pos/neg;}
//   else{majCInfo.majCRatio = 10.0;}

//   if (!pos){
//     majCInfo.smoothMajC = -10.0;
//   }
//   else {
//     if (!neg) {
//       majCInfo.smoothMajC = 10.0;
//     }
//     else{
//       majCInfo.smoothMajC = TMath::Log(pos/neg);
//     }
//   }

//   if (pos > neg) {
//     if (neg) {majCInfo.majC = pos/neg; return majCInfo;}
//     else {majCInfo.majC = 10.0; return majCInfo;}
//   }
//   if (neg > pos) {
//     if (pos) {majCInfo.majC = -1.0*neg/pos; return majCInfo;}
//     else {majCInfo.majC = -10.0; return majCInfo;}
//   }
//   if (neg == pos) {majCInfo.majC = 0.0; return majCInfo;}
//   else {majCInfo.majC = 0.0; return majCInfo;}
//   return majCInfo;
// }

---