From 938286ffb01fe4588a9eea8b61a4d7bc19e1470f Mon Sep 17 00:00:00 2001
From: Nicolas de Sereville <deserevi@ipno.in2p3.fr>
Date: Fri, 4 Dec 2015 11:23:28 +0100
Subject: [PATCH] + Add TSplitPoleNMR class holding the magnetic field
information
+ Add support for using star/stop run times
---
NPAnalysis/SPcoincW1/RunToTreat.txt | 4 +-
NPLib/SplitPole/CMakeLists.txt | 7 +-
NPLib/SplitPole/TSplitPoleNMR.cxx | 289 ++++++++++++++++++++++++++
NPLib/SplitPole/TSplitPoleNMR.h | 83 ++++++++
NPLib/SplitPole/TSplitPolePhysics.cxx | 61 +++++-
NPLib/SplitPole/TSplitPolePhysics.h | 20 ++
6 files changed, 459 insertions(+), 5 deletions(-)
create mode 100644 NPLib/SplitPole/TSplitPoleNMR.cxx
create mode 100644 NPLib/SplitPole/TSplitPoleNMR.h
diff --git a/NPAnalysis/SPcoincW1/RunToTreat.txt b/NPAnalysis/SPcoincW1/RunToTreat.txt
index 62739fccd..3fd8d3e5e 100644
--- a/NPAnalysis/SPcoincW1/RunToTreat.txt
+++ b/NPAnalysis/SPcoincW1/RunToTreat.txt
@@ -6,5 +6,7 @@ RootFileName
% /scratch/gypaos/data/al27pp/oct2015/midas/root/R50_0_SPonly.root
% /scratch/gypaos/data/al27pp/oct2015/midas/root/R50_0.root
% /scratch/gypaos/data/al27pp/oct2015/midas/root/R131_0_test.root
- /scratch/gypaos/data/al27pp/oct2015/midas/root/R108_0_test.root
+% /scratch/gypaos/data/al27pp/oct2015/midas/root/R108_0_test.root
+ /scratch/gypaos/data/al27pp/oct2015/midas/root/R74_0.root
+ /scratch/gypaos/data/al27pp/oct2015/midas/root/R72_0.root
diff --git a/NPLib/SplitPole/CMakeLists.txt b/NPLib/SplitPole/CMakeLists.txt
index 796cbb1cd..079ef004c 100644
--- a/NPLib/SplitPole/CMakeLists.txt
+++ b/NPLib/SplitPole/CMakeLists.txt
@@ -1,6 +1,7 @@
add_custom_command(OUTPUT TSplitPolePhysicsDict.cxx COMMAND ../scripts/build_dict.sh TSplitPolePhysics.h TSplitPolePhysicsDict.cxx TSplitPolePhysics.rootmap libNPSplitPole.dylib DEPENDS TSplitPolePhysics.h)
add_custom_command(OUTPUT TSplitPoleDataDict.cxx COMMAND ../scripts/build_dict.sh TSplitPoleData.h TSplitPoleDataDict.cxx TSplitPoleData.rootmap libNPSplitPole.dylib DEPENDS TSplitPoleData.h)
-add_library(NPSplitPole SHARED TSplitPoleSpectra.cxx TSplitPoleData.cxx TSplitPolePhysics.cxx TSplitPoleDataDict.cxx TSplitPolePhysicsDict.cxx )
-target_link_libraries(NPSplitPole ${ROOT_LIBRARIES} NPCore)
-install(FILES TSplitPoleData.h TSplitPolePhysics.h TSplitPoleSpectra.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
+add_custom_command(OUTPUT TSplitPoleNMRDict.cxx COMMAND ../scripts/build_dict.sh TSplitPoleNMR.h TSplitPoleNMRDict.cxx TSplitPoleNMR.rootmap libNPSplitPole.dylib DEPENDS TSplitPoleNMR.h)
+add_library(NPSplitPole SHARED TSplitPoleSpectra.cxx TSplitPoleData.cxx TSplitPolePhysics.cxx TSplitPoleNMR.cxx TSplitPoleDataDict.cxx TSplitPolePhysicsDict.cxx TSplitPoleNMRDict.cxx)
+target_link_libraries(NPSplitPole ${ROOT_LIBRARIES} NPCore)
+install(FILES TSplitPoleData.h TSplitPolePhysics.h TSplitPoleSpectra.h TSplitPoleNMR.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
diff --git a/NPLib/SplitPole/TSplitPoleNMR.cxx b/NPLib/SplitPole/TSplitPoleNMR.cxx
new file mode 100644
index 000000000..ecf588ae6
--- /dev/null
+++ b/NPLib/SplitPole/TSplitPoleNMR.cxx
@@ -0,0 +1,289 @@
+/*****************************************************************************
+ * Original Author: N. de Sereville contact address: deserevi@ipno.in2p3.fr *
+ * *
+ * Creation Date : april 2012 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: This class reads the rmn file and provides some simple *
+ * operations: apply delay, calculate mean value for subset of *
+ * the data, .... *
+ * Operations support both absolute and relative time *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * + Based on cSplitRun class from Iulian Stefan *
+ * *
+ *****************************************************************************/
+
+// class header
+#include "TSplitPoleNMR.h"
+
+// ROOT header
+#include "TMath.h"
+
+// C++ header
+#include <sstream>
+using namespace std;
+
+//#define DEBUG
+
+
+ClassImp(TSplitPoleNMR)
+
+
+TSplitPoleNMR::TSplitPoleNMR()
+ : fOpenFileTime(0),
+ fDelay(6797),
+ fFileName("Data_run.dat"),
+ fRmnRelativeTime(new TGraph()),
+ fRmn(new TGraph()),
+ fIsRelativeTime(0),
+ fIsLargeField(0),
+ fMean(-1),
+ fMin(10),
+ fMax(-1)
+{
+}
+
+
+
+TSplitPoleNMR::TSplitPoleNMR(const char* fileName)
+ : fOpenFileTime(0),
+ fDelay(6797),
+ fFileName(fileName),
+ fRmnRelativeTime(new TGraph()),
+ fRmn(new TGraph()),
+ fIsRelativeTime(0),
+ fIsLargeField(0),
+ fMean(-1),
+ fMin(10),
+ fMax(-1)
+{
+ ReadRmnFile();
+}
+
+
+
+TSplitPoleNMR::TSplitPoleNMR(const char* fileName, Double_t delay)
+ : fOpenFileTime(0),
+ fDelay(delay),
+ fFileName(fileName),
+ fRmnRelativeTime(new TGraph()),
+ fRmn(new TGraph()),
+ fIsRelativeTime(0),
+ fIsLargeField(0),
+ fMean(-1),
+ fMin(10),
+ fMax(-1)
+{
+ ReadRmnFile();
+}
+
+
+
+TSplitPoleNMR::~TSplitPoleNMR()
+{
+ if (fRmnRelativeTime) delete fRmnRelativeTime;
+ if (fRmn) delete fRmn;
+}
+
+
+
+
+Int_t TSplitPoleNMR::ReadRmnFile()
+{
+ ifstream in_rmn;
+ Double_t old_x=-10.;
+ std::cout<<"\nReading Rmn data file: "<<fFileName<<"..."<<std::flush;
+// DeleteRmnGraph();
+// fRmn= new TGraph();
+ in_rmn.open(fFileName);
+
+ if(!in_rmn.is_open()){
+ std::cout<<"\n Cant open Rmn data file "<<fFileName<<" ...\n"<<std::flush;
+ std::cout<<"\n Nor fRmnFile "<<fFileName<<" ...\n"<<std::flush;
+
+ return -1;
+ }
+
+ char tmp[1000];
+ in_rmn.getline(tmp,sizeof(tmp),'\n');
+ // TTimeStamp a;
+ in_rmn.getline(tmp,sizeof(tmp),':');
+ Int_t year,month,day,hour,minutes,seconds;
+
+ char ch_tmp;
+ in_rmn>>year;in_rmn>>ch_tmp;in_rmn>>month;in_rmn>>ch_tmp;in_rmn>>day;
+ in_rmn.getline(tmp,sizeof(tmp),'-');
+ in_rmn>>hour;in_rmn>>ch_tmp;in_rmn>>minutes;in_rmn>>ch_tmp;in_rmn>>seconds;
+
+ fOpenFileTime.Set(year,month,day,hour,minutes,seconds, 0, 1, 0);
+ fOpenFileTime.Print();
+ Double_t time=fOpenFileTime.AsDouble(),x=1.,y;
+ Int_t i=0;
+
+ Double_t ymin = 10;
+ Double_t ymax = -1;
+ Double_t ymean = 0;
+
+#ifdef DEBUG
+// std::cout << "time 0 = " << time << std::endl;
+// std::cout<<"\nx="<<x<<"\ty="<<y<<std::flush;
+#endif
+
+ in_rmn.getline(tmp,sizeof(tmp),'\n');
+ in_rmn.getline(tmp,sizeof(tmp),'\n');
+ for(;;){
+ in_rmn>>x;in_rmn>>ch_tmp;in_rmn>>y;
+ y*=1.e-7;
+ // apply large field
+ if (fIsLargeField) y += 1;
+#ifdef DEBUG
+// std::cout<<"\nx="<<x<<"\ty="<<y<<std::flush;
+#endif
+ if(x > old_x) {
+ // TGraph
+ fRmnRelativeTime -> SetPoint(i, x, y);
+ fRmn -> SetPoint(i++, x+time+fDelay, y);
+#ifdef DEBUG
+// std::cout<<"\nx="<<x+time+fDelay<<"\ty="<<y<<std::flush;
+#endif
+ // min, max
+ if (y < ymin) ymin = y;
+ if (y > ymax) ymax = y;
+ // mean
+ ymean += y;
+ }
+ else break;
+ old_x=x;
+ x=0.;
+
+ }
+ in_rmn.close();
+ std::cout<<" Ok!"<<std::flush;
+ std::cout << std::endl;
+
+ // set mean, min and max field values
+ fMin = ymin;
+ fMax = ymax;
+ fMean = ymean / fRmnRelativeTime->GetN();
+
+ return 0;
+}
+
+
+
+void TSplitPoleNMR::SetIsRelativeTime(bool isRelativeTime)
+{
+ fIsRelativeTime = isRelativeTime;
+}
+
+
+
+void TSplitPoleNMR::SetLargeField(bool isLargeField)
+{
+ fIsLargeField = isLargeField;
+ ReadRmnFile();
+}
+
+
+
+TGraph* TSplitPoleNMR::GetRmnGraph() const
+{
+ if (!fIsRelativeTime) return fRmn;
+ else return fRmnRelativeTime;
+}
+
+
+
+void TSplitPoleNMR::Draw(const char* option) const
+{
+ if (fRmn->GetN() > 0) {
+ if (!fIsRelativeTime) fRmn->Draw(option);
+ else fRmnRelativeTime->Draw(option);
+ }
+}
+
+
+
+void TSplitPoleNMR::ApplyDelay()
+{
+ if (!fRmnRelativeTime) {
+ cout << "TGraph containing rmn data is missing" << endl;
+ }
+ else {
+ Double_t timeRef = fOpenFileTime.AsDouble();
+ Int_t dim = fRmnRelativeTime->GetN();
+ Double_t *x = fRmnRelativeTime->GetX();
+ Double_t *y = fRmnRelativeTime->GetY();
+ for (Int_t i = 0; i < dim; ++i) { // loop on TGraph size
+ fRmn->SetPoint(i, x[i] + fDelay + timeRef, y[i]);
+ } // end loop on TGraph size
+ }
+}
+
+
+
+void TSplitPoleNMR::Dump()
+{
+ std::cout<<"\nPrinting TSplitPoleNMR object "<<this<<"\n"<<std::flush;
+
+
+
+ if(fRmn){
+ std::cout<<"\nfRmn defined "<<fRmn->GetName();
+ std::cout<<"\nfrom file "<<fFileName;
+ }
+ else
+ std::cout<<"\nNo fRmn object";
+
+}
+
+
+
+Double_t TSplitPoleNMR::EvalB(Double_t Time) const
+{
+ if(!fRmn) return 0.;
+ Double_t *x;
+ x=fRmn->GetX();
+ Int_t n= fRmn->GetN();
+
+ if( Time>=x[0] && Time<=x[n-1])
+ return fRmn->Eval(Time);
+ else
+ return fRmn->GetMean(2);
+
+}
+
+
+
+Double_t TSplitPoleNMR::Mean(Double_t t1, Double_t t2) const
+{
+ // get dimension and x array
+ Int_t dim;
+ Double_t *x, *y;
+ if (fIsRelativeTime) {
+ dim = fRmnRelativeTime->GetN();
+ x = fRmnRelativeTime->GetX();
+ y = fRmnRelativeTime->GetY();
+ }
+ else {
+ dim = fRmn->GetN();
+ x = fRmn->GetX();
+ y = fRmn->GetY();
+ }
+
+ // find index for t1 and t2
+ Double_t index1 = TMath::BinarySearch(dim, x, t1);
+ Double_t index2 = TMath::BinarySearch(dim, x, t2);
+ // if t1=0, BinarySearh returns index = -1, we force it to 0
+ if (index1 < 0) index1 = 0;
+
+ // calculate mean
+ Double_t mean = 0;
+ for (Int_t i = index1; i < index2+1; ++i) { // loop on field TGraph
+ mean += y[i];
+ } // end loop on field TGraph
+
+
+ return mean / (index2 - index1 + 1);
+}
diff --git a/NPLib/SplitPole/TSplitPoleNMR.h b/NPLib/SplitPole/TSplitPoleNMR.h
new file mode 100644
index 000000000..0e78cbbe5
--- /dev/null
+++ b/NPLib/SplitPole/TSplitPoleNMR.h
@@ -0,0 +1,83 @@
+#ifndef __SPLITPOLENMR__
+#define __SPLITPOLENMR__
+/*****************************************************************************
+ * Copyright (C) 2009-2014 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: N. de Sereville contact address: deserevi@ipno.in2p3.fr *
+ * *
+ * Creation Date : april 2012 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: This class reads the rmn file and provides some simple *
+ * operations: apply delay, calculate mean value for subset of *
+ * the data, .... *
+ * Operations support both absolute and relative time *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * + Based on cSplitRun class from Iulian Stefan *
+ * *
+ *****************************************************************************/
+
+// ROOT headers
+#include "TObject.h"
+#include "TTimeStamp.h"
+#include "TGraph.h"
+#include "TString.h"
+
+// C++ headers
+#include <iostream>
+#include <fstream>
+using namespace std;
+
+
+class TSplitPoleNMR : public TObject
+{
+ private:
+ TTimeStamp fOpenFileTime; //
+ Double_t fDelay; // Trun-Trmn in sec
+ TString fFileName;
+ TGraph* fRmnRelativeTime; // in sec wrt file beginning
+ TGraph* fRmn; // absolute time taking into account fDelay
+ bool fIsRelativeTime; // specify if returned graph is in relative of absolute units
+ bool fIsLargeField; // specifiy if field is greater than 1 T
+ Double_t fMean; // mean value for field
+ Double_t fMin; // min value for field
+ Double_t fMax; // max value for field
+
+
+ public:
+ TSplitPoleNMR();
+ TSplitPoleNMR(const char* fileName);
+ TSplitPoleNMR(const char* fileName, Double_t delay);
+ virtual ~TSplitPoleNMR();
+
+ void SetDelay(Double_t delay) {fDelay = delay; ApplyDelay();}
+ void SetIsRelativeTime(bool isRelativeTime);
+ void SetLargeField(bool isLargeField);
+ void SetFileName(TString fileName) {fFileName = fileName;}
+
+ Double_t GetDelay() {return fDelay;}
+ Double_t GetMean() {return fMean;}
+ Double_t GetMin() {return fMin;}
+ Double_t GetMax() {return fMax;}
+ TString GetFileName() {return fFileName;}
+
+ Int_t ReadRmnFile();
+ TGraph* GetRmnGraph() const;
+ void Draw(const char* option) const;
+ void ApplyDelay();
+
+ void Dump();
+
+ Double_t EvalB(Double_t Time) const;
+ Double_t Mean(Double_t t1, Double_t t2) const;
+
+ ClassDef(TSplitPoleNMR,1);
+};
+
+#endif
diff --git a/NPLib/SplitPole/TSplitPolePhysics.cxx b/NPLib/SplitPole/TSplitPolePhysics.cxx
index fe18fe437..d81569a9f 100644
--- a/NPLib/SplitPole/TSplitPolePhysics.cxx
+++ b/NPLib/SplitPole/TSplitPolePhysics.cxx
@@ -12,7 +12,7 @@
* Last update : *
*---------------------------------------------------------------------------*
* Decription: *
- * This class hold SplitPole Physics *
+ * This class hold SplitPole Physics *
* *
*---------------------------------------------------------------------------*
* Comment: *
@@ -48,12 +48,21 @@ TSplitPolePhysics::TSplitPolePhysics()
: m_EventData(new TSplitPoleData),
m_PreTreatedData(new TSplitPoleData),
m_EventPhysics(this),
+ m_RunStart(2015, 10, 6, 0, 0, 0),
+ m_RunStop(2015, 10, 7, 0, 0, 0),
+ m_RunLength(0),
+ m_FrequenceClock(2.03),
+ m_TickMin(0),
+ m_TickMax(0),
+ m_RunNumber(0),
+ m_CurrentRunNumber(0),
m_MagneticFieldCorrection(0),
m_TimeDelay(6500),
m_CalibP0(0.65),
m_CalibP1(6e-6)
{
Clear();
+ ReadTimeTable();
}
@@ -133,6 +142,44 @@ void TSplitPolePhysics::ReadConfiguration(string Path)
}
+
+///////////////////////////////////////////////////////////////////////////
+void TSplitPolePhysics::ReadTimeTable()
+{
+ ifstream file("TimeTable.txt");
+
+ Int_t run, date1, date2, time1, time2;
+ TTimeStamp start, stop;
+ pair<TTimeStamp, TTimeStamp> time;
+ while (!file.eof()) {
+ file >> run >> date1 >> time1 >> date2 >> time2;
+ cout << run << "\t" << date1 << "\t" << time1 << "\t" << date2 << "\t" << time2 << endl;
+ start.Set(date1, time1, 0, 1, 0);
+ stop.Set(date2, time2, 0, 1, 0);
+ cout << "dt = " << stop.GetSec() - start.GetSec() << endl;
+ time.first = start;
+ time.second = stop;
+ m_TimeTable[run] = time;
+ }
+
+ cout << m_TimeTable[74].first << endl;
+ cout << m_TimeTable.size() << endl;
+
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////
+Bool_t TSplitPolePhysics::IsSameRun()
+{
+ Bool_t isSameRun = true;
+ if (m_CurrentRunNumber != m_RunNumber) isSameRun = false;
+
+ return isSameRun;
+}
+
+
+
///////////////////////////////////////////////////////////////////////////
void TSplitPolePhysics::AddParameterToCalibrationManager()
{
@@ -151,6 +198,7 @@ void TSplitPolePhysics::InitializeRootInputRaw()
inputChain->SetBranchStatus("fPlasticP", true);
inputChain->SetBranchStatus("fPlasticG", true);
inputChain->SetBranchAddress("SplitPole", &m_EventData);
+ inputChain->SetBranchAddress("RunNumber", &m_RunNumber);
}
@@ -209,6 +257,17 @@ void TSplitPolePhysics::BuildSimplePhysicalEvent()
fTime2.push_back(m_EventData->GetTime2(i));
} // end loop on multiplicity
+
+ // Magnetic field correction
+ // store localy run number, run start and stop times and rmn data
+// cout << "rrrrrrrrrrun number = " << m_RunNumber << endl;
+ if (!IsSameRun()) {
+ m_CurrentRunNumber = m_RunNumber;
+ m_RunStart = m_TimeTable[m_CurrentRunNumber].first;
+ m_RunStop = m_TimeTable[m_CurrentRunNumber].second;
+ m_RunLength = m_RunStop.AsDouble() - m_RunStart.AsDouble();
+// cout << m_CurrentRunNumber << "\t" << m_RunStart << "\t" << m_RunStop << "\t" << m_RunLength << endl;
+ }
// Correct for magnetic field variation
// fBrho = (m_CalibP0 + m_CalibP1*m_EventData->GetPosition()) * 0.5;
fBrho = (m_CalibP0 + m_CalibP1*m_EventData->GetPlasticG()) * 0.5;
diff --git a/NPLib/SplitPole/TSplitPolePhysics.h b/NPLib/SplitPole/TSplitPolePhysics.h
index e4f02655f..541938284 100644
--- a/NPLib/SplitPole/TSplitPolePhysics.h
+++ b/NPLib/SplitPole/TSplitPolePhysics.h
@@ -23,12 +23,14 @@
*****************************************************************************/
// STL
#include <vector>
+#include <utility>
#include <map>
using namespace std;
// ROOT
#include "TObject.h"
#include "TH1.h"
+#include "TTimeStamp.h"
// NPL
#include "TSplitPoleData.h"
@@ -162,12 +164,30 @@ class TSplitPolePhysics : public TObject, public NPL::VDetector
// map< int, vector<bool> > m_BackChannelStatus; //!
+ private: // parameters needed for magnetic field correction
+ map<Int_t, pair<TTimeStamp, TTimeStamp> > m_TimeTable; //!
+ TTimeStamp m_RunStart; //!
+ TTimeStamp m_RunStop; //!
+ Double_t m_RunLength; //! // in sec
+ Double_t m_FrequenceClock; //! // in Hz
+ Double_t m_TickMin; //!
+ Double_t m_TickMax; //!
+ Int_t m_RunNumber; //! // read event by event from TTree
+ Int_t m_CurrentRunNumber; //!
+
private: // Parameters used in the analysis
Bool_t m_MagneticFieldCorrection; //!
Double_t m_TimeDelay; //!
Double_t m_CalibP0; //!
Double_t m_CalibP1; //!
+ // methods for magnetic field correction
+ public: // called once
+ void ReadTimeTable();
+
+ public: // called event by event
+ Bool_t IsSameRun();
+
private: // Spectra Class
TSplitPoleSpectra* m_Spectra; // !
--
GitLab