NuGeneral.cxx

Go to the documentation of this file.

// Coded by Jeff Hartnell Jul/2005       
//
// Contact: j.j.hartnell@rl.ac.uk

#include <cmath>
#include <cassert>
#include <map>
#include <sys/stat.h>

#include "TDirectory.h"
#include "TGraph.h"
#include "TF1.h"
#include "TH1F.h"
#include "TMath.h"
#include "TProfile2D.h"
#include "TROOT.h"

#include "MessageService/MsgService.h"

#include "NtupleUtils/NuGeneral.h"

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

CVSID("$Id: NuGeneral.cxx,v 1.8 2010/02/01 11:39:48 evansj Exp $");

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

NuGeneral::NuGeneral()
{
  MSG("NuGeneral",Msg::kDebug)
    <<"Running NuGeneral Constructor..."<<endl;


  MSG("NuGeneral",Msg::kDebug)
    <<"Finished NuGeneral Constructor"<<endl;
}

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

NuGeneral::~NuGeneral()
{
  MSG("NuGeneral",Msg::kDebug)
    <<"Running NuGeneral Destructor..."<<endl;
  

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

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

std::string NuGeneral::GetReleaseDirToUse(std::string sPackageDir) const
{

  const string sPub=this->GetBaseReleaseDir();
  const string sPriv=this->GetTestReleaseDir();

  //set default return value to sPriv
  string sAnaDir=sPriv;
  
  //try private context first
  string sTmp=sPriv+"/"+sPackageDir;
  struct stat ss;
  MAXMSG("NuGeneral",Msg::kDebug,10)<<"sTmp="<<sTmp<<endl;
  if(stat(sTmp.c_str(),&ss)==-1) {
    MAXMSG("NuGeneral",Msg::kInfo,10)
      <<"Package="<<sPackageDir
      <<" not present in SRT_PRIVATE_CONTEXT="<<sPriv
      <<", trying SRT_PUBLIC_CONTEXT="<<sPub<<"..."<<endl;
    
    //now test public context
    sTmp=sPub+"/"+sPackageDir;
    MAXMSG("NuGeneral",Msg::kDebug,10)<<"sTmp2="<<sTmp<<endl;
    if(stat(sTmp.c_str(),&ss)==-1) {
      MAXMSG("NuGeneral",Msg::kWarning,10)
        <<"Package="<<sPackageDir
        <<" not present in SRT_PUBLIC_CONTEXT="<<sPub
        <<" or SRT_PRIVATE_CONTEXT="<<sPriv
        <<", will abort here..."<<endl;
      assert(false);
    }
    //public context exists so use that
    else {
      MAXMSG("NuGeneral",Msg::kInfo,10)
        <<"Using SRT_PUBLIC_CONTEXT="<<sPub
        <<" for package="<<sPackageDir<<endl;
      sAnaDir=sPub;
    }
  }
  else {
    MAXMSG("NuGeneral",Msg::kInfo,10)
      <<"Using SRT_PRIVATE_CONTEXT="<<sPriv
      <<" for package="<<sPackageDir<<endl;
  }
  
  return sAnaDir;
}

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

std::string NuGeneral::GetTestReleaseDir() const
{
  //get the environmental variable
  char* anaDir=getenv("SRT_PRIVATE_CONTEXT");
  
  //use a string to hold env instead 
  string sAnaDir="";
  
  if (anaDir!=NULL) {
    sAnaDir=anaDir;
  }
  else {
    MSG("NuGeneral",Msg::kError) 
      <<"Environmental variable $SRT_PRIVATE_CONTEXT not set"<<endl
      <<"Will abort here..."<<endl;
    assert(false);
  }

  return sAnaDir;
}

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

std::string NuGeneral::GetBaseReleaseDir() const
{
  //get the environmental variable
  char* anaDir=getenv("SRT_PUBLIC_CONTEXT");
  
  //use a string to hold env instead 
  string sAnaDir="";
  
  if (anaDir!=NULL) {
    sAnaDir=anaDir;
  }
  else {
    MSG("NuGeneral",Msg::kError) 
      <<"Environmental variable $SRT_PUBLIC_CONTEXT not set"<<endl
      <<"Will abort here..."<<endl;
    assert(false);
  }
  
  return sAnaDir;
}

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

void NuGeneral::CopyAcrossObjects(TDirectory* dirInput,
                                  TDirectory* dirOutput,
                                  const std::vector<std::string>& 
                                  vObjectNames) const
{
  cout<<"NuGeneral: copying across objects"<<endl;

  TDirectory* directoryOrig=gDirectory;
  //cout<<"Original directory is:"<<endl;
  //directoryOrig->Print();
  
  //change to the directory where the files are kept
  gDirectory=dirInput;
  //cout<<"Changed to input directory:"<<endl;
  //gDirectory->Print();
  
  //vector to store the pointers to the histograms to copy across
  vector<TObject*> vObjects;
  
  //get the pointers
  for (vector<string>::const_iterator iter=vObjectNames.begin();
       iter!=vObjectNames.end();++iter){
    TObject* h=gROOT->FindObject(iter->c_str());
    vObjects.push_back(h);
    if (!h) {
      MAXMSG("NuGeneral",Msg::kInfo,10)
        <<"CopyAcrossObjects: Can't find object="<<*iter<<endl;
    }
  }
  
  //change to the directory where the objects are to be written
  gDirectory=dirOutput;
  //cout<<"Changed to output directory:"<<endl;
  //gDirectory->Print();

  //loop over the objects and write them out
  for (vector<TObject*>::const_iterator iter=vObjects.begin();
       iter!=vObjects.end();++iter){
    if (*iter) {
      (*iter)->Write();
    }
  }

  //reset the global directory to the original directory
  gDirectory=directoryOrig;
}

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

void NuGeneral::SetGraphAxisEtc(TGraph* g,Int_t startTimeSecs,
                                Int_t endTimeSecs) const 
{
  MSG("NuGeneral",Msg::kVerbose)
    <<endl<<" ** Running the SetGraphAxisEtc method... ** "<<endl;

  TDatime startDatime;
  startDatime.Set(startTimeSecs+788918400);
  TDatime endDatime;
  endDatime.Set(endTimeSecs+788918400);

  //play around with the formatting a little
  string sStartHour=Form("%d",startDatime.GetHour());
  string sStartMinute=Form("%d",startDatime.GetMinute());
  string sStartSecond=Form("%d",startDatime.GetSecond());
  string sStartYear=Form("%d",startDatime.GetYear());
  string sStartMonth=Form("%d",startDatime.GetMonth());
  string sStartDay=Form("%d",startDatime.GetDay());
  if (startDatime.GetHour()<10) sStartHour="0"+sStartHour;
  if (startDatime.GetMinute()<10) sStartMinute="0"+sStartMinute;
  if (startDatime.GetSecond()<10) sStartSecond="0"+sStartSecond;
  if (startDatime.GetMonth()<10) sStartMonth="0"+sStartMonth;
  if (startDatime.GetDay()<10) sStartDay="0"+sStartDay;

  //play around with the formatting a little
  string sEndHour=Form("%d",endDatime.GetHour());
  string sEndMinute=Form("%d",endDatime.GetMinute());
  string sEndSecond=Form("%d",endDatime.GetSecond());
  string sEndYear=Form("%d",endDatime.GetYear());
  string sEndMonth=Form("%d",endDatime.GetMonth());
  string sEndDay=Form("%d",endDatime.GetDay());
  if (endDatime.GetHour()<10) sEndHour="0"+sEndHour;
  if (endDatime.GetMinute()<10) sEndMinute="0"+sEndMinute;
  if (endDatime.GetSecond()<10) sEndSecond="0"+sEndSecond;
  if (endDatime.GetMonth()<10) sEndMonth="0"+sEndMonth;
  if (endDatime.GetDay()<10) sEndDay="0"+sEndDay;

  static Bool_t firstTime=true;
  if (firstTime){
    MSG("NuGeneral",Msg::kInfo) 
      <<"Setting graph axis info:"<<endl;
    MSG("NuGeneral",Msg::kInfo) 
      <<"  Start time "<<sStartHour
      <<":"<<sStartMinute
      <<":"<<sStartSecond
      <<" on "<<startDatime.GetYear()
      <<"/"<<sStartMonth
      <<"/"<<sStartDay
      <<endl;
    
    MSG("NuGeneral",Msg::kInfo) 
      <<"  End time "<<sEndHour
      <<":"<<sEndMinute
      <<":"<<sEndSecond
      <<" on "<<endDatime.GetYear()
      <<"/"<<sEndMonth
      <<"/"<<sEndDay
      <<endl;
  }
  firstTime=false;

  Int_t timeRange=endTimeSecs-startTimeSecs;

  //check graph pointer is not null
  if (g){
    if (sStartDay==sEndDay && timeRange<1*24*60*60){
      string title="Time ("+sStartYear+"/"+sStartMonth+"/"+sStartDay+")";
      g->GetXaxis()->SetTitle(title.c_str());
      g->GetXaxis()->SetTimeFormat("%H:%M");
    }
    else if (timeRange>1*24*60*60 && timeRange<2*24*60*60){//1-2 days
      string title="Time ("+sStartYear+"/"+sStartMonth+"/"+sStartDay
        +" - "+sEndYear+"/"+sEndMonth+"/"+sEndDay+")";
      g->GetXaxis()->SetTitle(title.c_str());
      g->GetXaxis()->SetTimeFormat("%H:%M");
    }
    else if (timeRange>2*24*60*60 && timeRange<5*24*60*60){//2-5 days
      string title="Time ("+sStartYear+"/"+sStartMonth+"/"+sStartDay
        +" - "+sEndYear+"/"+sEndMonth+"/"+sEndDay+")";
      g->GetXaxis()->SetTitle(title.c_str());
      g->GetXaxis()->SetTimeFormat("%H:%M-%d/%m");
    }
    else{
      string title="Time ("+sStartYear+"/"+sStartMonth+"/"+sStartDay
        +" - "+sEndYear+"/"+sEndMonth+"/"+sEndDay+")";
      g->GetXaxis()->SetTitle(title.c_str());
      g->GetXaxis()->SetTimeFormat("%H:%M-%d/%m");
      //stop it cramming loads of days on the axis!
      g->GetXaxis()->SetNdivisions(506);
      //n = N1 + 100*N2 + 10000*N3
      //N1=number of primary divisions.
      //N2=number of secondary divisions.
      //N3=number of 3rd divisions.
    }
    
    //other x-axis defaults
    g->GetXaxis()->SetTimeDisplay(1);
    g->GetXaxis()->CenterTitle();
    
    //other defaults
    g->SetMarkerStyle(3);
    g->SetMarkerColor(2);
    g->SetMarkerSize(0.35);
    g->SetLineColor(46);
  }
  else{
    MSG("NuGeneral",Msg::kError)
      <<"Input graph pointer is null!"
      <<" Will do nothing."<<endl;
  }

  MSG("NuGeneral",Msg::kVerbose)
    <<endl<<" ** Finished the SetGraphAxisEtc method... **"<<endl;
}

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

void NuGeneral::SetGraphAxis(TAxis* a,Int_t startTimeSecs,
                             Int_t endTimeSecs) const 
{
  if (!a){
    MSG("NuGeneral",Msg::kWarning)
      <<"Input axis pointer is null! Will do nothing"<<endl;
    return;
  }

  //get datimes
  TDatime startDatime;
  startDatime.Set(startTimeSecs+788918400);
  TDatime endDatime;
  endDatime.Set(endTimeSecs+788918400);

  //play around with the formatting a little
  string sStartHour=Form("%d",startDatime.GetHour());
  string sStartMinute=Form("%d",startDatime.GetMinute());
  string sStartSecond=Form("%d",startDatime.GetSecond());
  string sStartYear=Form("%d",startDatime.GetYear());
  string sStartMonth=Form("%d",startDatime.GetMonth());
  string sStartDay=Form("%d",startDatime.GetDay());
  if (startDatime.GetHour()<10) sStartHour="0"+sStartHour;
  if (startDatime.GetMinute()<10) sStartMinute="0"+sStartMinute;
  if (startDatime.GetSecond()<10) sStartSecond="0"+sStartSecond;
  if (startDatime.GetMonth()<10) sStartMonth="0"+sStartMonth;
  if (startDatime.GetDay()<10) sStartDay="0"+sStartDay;

  //play around with the formatting a little
  string sEndHour=Form("%d",endDatime.GetHour());
  string sEndMinute=Form("%d",endDatime.GetMinute());
  string sEndSecond=Form("%d",endDatime.GetSecond());
  string sEndYear=Form("%d",endDatime.GetYear());
  string sEndMonth=Form("%d",endDatime.GetMonth());
  string sEndDay=Form("%d",endDatime.GetDay());
  if (endDatime.GetHour()<10) sEndHour="0"+sEndHour;
  if (endDatime.GetMinute()<10) sEndMinute="0"+sEndMinute;
  if (endDatime.GetSecond()<10) sEndSecond="0"+sEndSecond;
  if (endDatime.GetMonth()<10) sEndMonth="0"+sEndMonth;
  if (endDatime.GetDay()<10) sEndDay="0"+sEndDay;

  static Bool_t firstTime=true;
  if (firstTime){
    MSG("NuGeneral",Msg::kInfo) 
      <<"Setting axis info:"<<endl;
    MSG("NuGeneral",Msg::kInfo) 
      <<"  Start time "<<sStartHour
      <<":"<<sStartMinute
      <<":"<<sStartSecond
      <<" on "<<startDatime.GetYear()
      <<"/"<<sStartMonth
      <<"/"<<sStartDay
      <<endl;

    MSG("NuGeneral",Msg::kInfo) 
      <<"  End time   "<<sEndHour
      <<":"<<sEndMinute
      <<":"<<sEndSecond
      <<" on "<<endDatime.GetYear()
      <<"/"<<sEndMonth
      <<"/"<<sEndDay
      <<endl;
  }

  Int_t timeRange=endTimeSecs-startTimeSecs;
   
  if (sStartDay==sEndDay && timeRange<1*24*60*60){
    string title="Time ("+sStartYear+"/"+sStartMonth+"/"+sStartDay+")";
    a->SetTitle(title.c_str());
    a->SetTimeFormat("%H:%M");
  }
  else if (timeRange>=1*24*60*60 && timeRange<2*24*60*60){//1-2 days
    string title="Time ("+sStartYear+"/"+sStartMonth+"/"+sStartDay
      +" - "+sEndYear+"/"+sEndMonth+"/"+sEndDay+")";
    a->SetTitle(title.c_str());
    a->SetTimeFormat("%H:%M");
  }
  else if (timeRange>=2*24*60*60 && timeRange<5*24*60*60){//2-5 days
    string title="Time ("+sStartYear+"/"+sStartMonth+"/"+sStartDay
      +" - "+sEndYear+"/"+sEndMonth+"/"+sEndDay+")";
    a->SetTitle(title.c_str());
    a->SetTimeFormat("%H:%M-%d/%m");
  }
  else if (timeRange>=5*24*60*60 && timeRange<20*24*60*60){//5-20 days
    string title="Time ("+sStartYear+"/"+sStartMonth+"/"+sStartDay
      +" - "+sEndYear+"/"+sEndMonth+"/"+sEndDay+")";
    a->SetTitle(title.c_str());
    a->SetTimeFormat("%d/%m");
  }
  else if (timeRange>=20*24*60*60 && timeRange<90*24*60*60){//20-90 days
    string title="Time ("+sStartYear+"/"+sStartMonth+"/"+sStartDay
      +" - "+sEndYear+"/"+sEndMonth+"/"+sEndDay+")";
    a->SetTitle(title.c_str());
    a->SetTimeFormat("%Y/%m/%d");
  }
  else if (timeRange>=90*24*60*60 && 
           timeRange<2*365*24*60*60){//50-2*365 days
    string title="Time ("+sStartYear+"/"+sStartMonth+"/"+sStartDay
      +" - "+sEndYear+"/"+sEndMonth+"/"+sEndDay+")";
    a->SetTitle(title.c_str());
    a->SetTimeFormat("%Y/%m");//Y for year with century
  }
  else{
    string title="Time ("+sStartYear+"/"+sStartMonth+"/"+sStartDay
      +" - "+sEndYear+"/"+sEndMonth+"/"+sEndDay+")";
    a->SetTitle(title.c_str());
    a->SetTimeFormat("%H:%M-%d/%m");
    //stop it cramming loads of days on the axis!
    a->SetNdivisions(506);
    //n = N1 + 100*N2 + 10000*N3
    //N1=number of primary divisions.
    //N2=number of secondary divisions.
    //N3=number of 3rd divisions.
  }
   
  //TDatime da(2003,02,28,12,00,00);
  //          gStyle->SetTimeOffset(da.Convert());

  if (firstTime) MSG("NuGeneral",Msg::kInfo) 
    <<"Set axis title to be '"<<a->GetTitle()<<"'"<<endl;

  //other x-axis defaults
  a->SetTimeDisplay(1);
  firstTime=false;
}

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

void NuGeneral::TH1FFill(TH1F* h,const std::vector<Double_t>& vX) const
{
  if (vX.empty()) {
    MSG("NuGeneral",Msg::kWarning) 
      <<"One or more vectors is empty"<<endl;
    return;
  }

  if (!h){
    MSG("NuGeneral",Msg::kWarning) 
      <<"pointer==0"<<endl;
    return;
  }

  UInt_t vSize=vX.size();
  for (UInt_t i=0;i<vSize;i++){
    h->Fill(vX[i]);
  }
}

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

void NuGeneral::TH1FFill(TH1F* h,const std::vector<Double_t>& vX,
                       const std::vector<Double_t>& vY) const
{

  if (vX.empty() || vY.empty()) {
    MSG("NuGeneral",Msg::kWarning) 
      <<"One or more vectors is empty"<<endl;
    return;
  }

  if (!h){
    MSG("NuGeneral",Msg::kWarning) 
      <<"pointer==0"<<endl;
    return;
  }

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

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

void NuGeneral::TProfileFill(TProfile* p,
                             std::vector<Double_t>& vX,
                             std::vector<Double_t>& vY) const
{
  MSG("NuGeneral",Msg::kDebug) 
    <<" ** Running FillTProfile method... **"<<endl;

  if (vX.empty()) {
    MSG("NuGeneral",Msg::kWarning) 
      <<"FillTProfile: One or more vectors is empty"<<endl;
    return;
  }

  if (!p){
    MSG("NuGeneral",Msg::kWarning) 
      <<"FillTProfile: TProfile pointer==0"<<endl;
    return;
  }

  if (vX.size()!=vY.size()) {
    MSG("NuGeneral",Msg::kWarning) 
      <<"FillTProfile: vectors different sizes; returning..."<<endl;
    return;
  }

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

  MSG("NuGeneral",Msg::kDebug) 
    <<" ** Finished FillTProfile method **"<<endl;
}

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

void NuGeneral::TProfile2DFill(TProfile2D* p,
                               std::vector<Double_t>& vX,
                               std::vector<Double_t>& vY,
                               std::vector<Double_t>& vZ) const
{
  if (vX.empty()) {
    MSG("NuGeneral",Msg::kWarning) 
      <<"FillTProfile2D: One or more vectors is empty"<<endl;
    return;
  }

  if (!p){
    MSG("NuGeneral",Msg::kWarning) 
      <<"FillTProfile2D: TProfile2D pointer==0"<<endl;
    return;
  }

  if (vX.size()!=vY.size() || vX.size()!=vZ.size()) {
    MSG("NuGeneral",Msg::kWarning) 
      <<"FillTProfile2D: vectors different sizes; returning..."<<endl;
    return;
  }

  UInt_t vSize=vX.size();
  for (UInt_t i=0;i<vSize;i++){
    p->Fill(vX[i],vY[i],vZ[i]);
  }
}

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

void NuGeneral::EpochTo1995(std::vector<Double_t>& t) const
{
  if (t.empty()) {
    MSG("NuGeneral",Msg::kWarning) 
      <<"EpochTo1995: One or more vectors is empty"<<endl;
    return;
  }

  UInt_t vSize=t.size();
  for (UInt_t i=0;i<vSize;i++){
    t[i]-=788918400;// seconds between 1970 and 1995
  }
}

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

void NuGeneral::EpochTo1995(Int_t& t) const
{
  t-=788918400;// seconds between 1970 and 1995
}

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

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

  if (vX.empty()) {
    MSG("NuGeneral",Msg::kWarning) 
      <<"TGraphVect: One or more vectors is empty"<<endl;
    return 0;
  }

  if (vX.size()!=vY.size()) {
    MSG("NuGeneral",Msg::kWarning) 
      <<"TGraphVect: vectors different sizes; returning..."<<endl;
    return 0;
  }

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

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

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

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

void NuGeneral::TGraphOffset(std::vector<TGraph*>& v,
                             Double_t factor) const
{
  Double_t max=-1e200;
  Double_t min=+1e200;

  //find the min and max
  for (vector<TGraph*>::iterator vIt=v.begin();vIt!=v.end();++vIt){
    for (Int_t i=0;i<(*vIt)->GetN();i++){
      Double_t x=0;
      Double_t y=0;
      (*vIt)->GetPoint(i,x,y);
      MSG("NuGeneral",Msg::kVerbose)
        <<"i="<<i<<", x="<<x<<", y="<<y<<endl;

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

  //change the scale
  for (vector<TGraph*>::iterator vIt=v.begin();vIt!=v.end();++vIt){
    for (Int_t i=0;i<(*vIt)->GetN();i++){
      Double_t x=0;
      Double_t y=0;
      (*vIt)->GetPoint(i,x,y);
      MSG("NuGeneral",Msg::kVerbose)
        <<"i="<<i<<", x="<<x<<", y="<<y<<endl;

      //scale the y value
      y-=min;
      y*=factor;//scale by factor, e.g. convert to nanosec
      (*vIt)->SetPoint(i,x,y);
    }
  }

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

  //convert min and max
  max=(max-min)*factor;  
  min=(min-min)*factor;
  MSG("NuGeneral",Msg::kDebug)
    <<"Converted range to min="<<min<<", max="<<max<<endl;
  
  //now set the min and max for 10% above and below size of range
  for (vector<TGraph*>::iterator vIt=v.begin();vIt!=v.end();++vIt){
    (*vIt)->SetMinimum(min-0.1*(max-min));
    (*vIt)->SetMaximum(max+0.1*(max-min));
  }
}

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

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

  for (vector<TGraph*>::iterator vIt=v.begin();vIt!=v.end();++vIt){
    for (Int_t i=0;i<(*vIt)->GetN();i++){
      Double_t x=0;
      Double_t y=0;
      (*vIt)->GetPoint(i,x,y);
      MSG("NuGeneral",Msg::kVerbose)
        <<"i="<<i<<", x="<<x<<", y="<<y<<endl;
  
      //find the min and max
      if (y>max) max=y;
      if (y<min) min=y;
    }
  }

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

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

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

TGraph* NuGeneral::TGraphMap(const std::map<Int_t,Float_t>& m) const
{
  TGraph* g=new TGraph(m.size());

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

  Int_t i=0;

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

    i++;
    mIter++;
  }
  return g;
}

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

void NuGeneral::NormaliseVector(std::vector<Double_t>& v,Int_t mode) const
{
  if (v.size()>0){

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

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

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

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

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

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

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

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

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

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

      delete h;
    }
  }
}

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

void NuGeneral::ScaleVector(std::vector<Double_t>& v,
                            Double_t scaleFactor) const
{
  for (UInt_t i=0;i<v.size();i++) v[i]*=scaleFactor;
}

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

void NuGeneral::FlipVector(std::vector<Double_t>& v) const
{
  vector<Double_t> tempV=v;
  for (UInt_t i=0;i<v.size();i++) v[i]=tempV[v.size()-1-i];
}

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

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

  Double_t floatPart=modf(x,&xi); 

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

  return static_cast<Int_t>(xi);
}

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

void NuGeneral::ScaleMap(std::map<Int_t,Float_t>& m,
                         const std::map<Int_t,Float_t>& scaleFactor) const
{
  if (m.size()==scaleFactor.size()){
    
    map<Int_t,Float_t>::iterator mIter=m.begin();
    map<Int_t,Float_t>::const_iterator scale=scaleFactor.begin();

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

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

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

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

// TODO CJB - this helper function is also in NuUtilities
Double_t NuGeneral::OscWeight(Double_t DeltaM2, Double_t Sin22Theta,
                              Double_t Enu) const
{
  //calculate oscillation probability
  const Double_t Distance = 735;
  Double_t weight=1;
  Double_t sinm=0;
  
  sinm=sin((1.27*DeltaM2*Distance)/Enu);
  weight=1.- Sin22Theta * sinm * sinm;

  return weight;
}

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

std::string NuGeneral::IdHEPAsString(Int_t IdHEP) const
{
  string sIdHEP="";
  
  Bool_t shortName=true;
  if (shortName) {
    if (IdHEP==22) sIdHEP="photon";
    
    else if (IdHEP==11) sIdHEP="e-";
    else if (IdHEP==-11) sIdHEP="e+";
    else if (IdHEP==13) sIdHEP="mu-";
    else if (IdHEP==-13) sIdHEP="mu+";
    else if (IdHEP==15) sIdHEP="tau-";
    else if (IdHEP==-15) sIdHEP="tau+";

    else if (IdHEP==12) sIdHEP="nue";
    else if (IdHEP==-12) sIdHEP="nuebar";
    else if (IdHEP==14) sIdHEP="numu";
    else if (IdHEP==-14) sIdHEP="numubar";
    else if (IdHEP==16) sIdHEP="nutaubar";
    else if (IdHEP==-16) sIdHEP="nutaubar";

    else if (IdHEP==211) sIdHEP="pi+";
    else if (IdHEP==-211) sIdHEP="pi-";
    else if (IdHEP==221) sIdHEP="eta+";
    else if (IdHEP==-221) sIdHEP="eta-";
    else if (IdHEP==111) sIdHEP="pi0";
    
    else if (IdHEP==213) sIdHEP="rho+";
    else if (IdHEP==-213) sIdHEP="rho-";
    else if (IdHEP==223) sIdHEP="omega+";
    else if (IdHEP==-223) sIdHEP="omega-";
    else if (IdHEP==113) sIdHEP="rho0";
    
    else if (IdHEP==321) sIdHEP="k+";
    else if (IdHEP==-321) sIdHEP="k-";
    else if (IdHEP==310) sIdHEP="k-short";
    else if (IdHEP==130) sIdHEP="k-long";
    else if (IdHEP==311) sIdHEP="k0";

    else if (IdHEP==411) sIdHEP="D+";
    else if (IdHEP==421) sIdHEP="D0";
    else if (IdHEP==431) sIdHEP="D??";
    else if (IdHEP==4122) sIdHEP="Lambda_c+";

    else if (IdHEP==2212) sIdHEP="p+";
    else if (IdHEP==-2212) sIdHEP="p-";
    else if (IdHEP==2112) sIdHEP="n";

    else if (IdHEP==1114) sIdHEP="delta-";
    else if (IdHEP==-2114) sIdHEP="delta0";
    else if (IdHEP==2214) sIdHEP="delta+";
    else if (IdHEP==2224) sIdHEP="delta++";
    
    else if (IdHEP==28) sIdHEP="geantino";
  }
  else {

  }
  
  if (sIdHEP==""){
    string sForm=Form("%d",IdHEP);
    sIdHEP="IdHEP="+sForm+"=?";
  }

  return sIdHEP;
}

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

---