gwm17/SPS_SABRE_EventBuilder
1
0
Fork 0
mirror of https://github.com/sesps/SPS_SABRE_EventBuilder.git synced 2024-11-10 20:38:51 -05:00
Code Issues Projects Releases Wiki Activity

Code review, with a focus on enforcing style and purging extraneous code. Changed name of GWMEventBuilder to EVBApp to better reflect class intent.

Browse Source
This commit is contained in:
Gordon McCann 2021-12-13 12:28:56 -05:00
parent 7c28e252cb
commit afcdd3bb8c
32 changed files with 1536 additions and 1258 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

39
include/ChannelMap.h
View File

@ -9,35 +9,38 @@
#ifndef CHANNELMAP_H
#define CHANNELMAP_H
struct Channel {
struct Channel
{
int detectorType; //What kind of detector we're looking at
int detectorID; //Which specific detector we're looking at
int detectorPart; //Which specific part we're looking at
};
//Detector part/type identifiers for use in the code
enum DetAttribute {
FOCALPLANE,
SCINTLEFT,
SCINTRIGHT,
ANODEFRONT,
ANODEBACK,
DELAYFR,
DELAYFL,
DELAYBR,
DELAYBL,
CATHODE,
MONITOR,
SABRERING = 88, //These are offset to avoid interference at the variable mapping phase
SABREWEDGE = 99 //Just don't add any new attributes with values greater than 88
enum class DetAttribute
{
FocalPlane,
ScintLeft,
ScintRight,
AnodeFront,
AnodeBack,
DelayFR,
DelayFL,
DelayBR,
DelayBL,
Cathode,
Monitor,
SabreRing = 88, //These are offset to avoid interference at the variable mapping phase
SabreWedge = 99 //Just don't add any new attributes with values greater than 88
};
class ChannelMap {
class ChannelMap
{
public:
typedef std::unordered_map<int, Channel> Containter;
typedef std::unordered_map<int, Channel>::iterator Iterator;
typedef std::unordered_map<DetAttribute, Channel> Containter;
typedef std::unordered_map<DetAttribute, Channel>::iterator Iterator;
ChannelMap();
ChannelMap(const std::string& filename);

27
include/CompassFile.h
View File

@ -15,7 +15,9 @@
#include "ShiftMap.h"
#include <memory>
class CompassFile {
class CompassFile
{
public:
CompassFile();
CompassFile(const std::string& filename);
@ -23,17 +25,18 @@ public:
~CompassFile();
void Open(const std::string& filename);
void Close();
inline bool IsOpen() { return m_file->is_open(); };
bool GetNextHit();
inline CompassHit GetCurrentHit() const { return m_currentHit; };
inline std::string GetName() { return m_filename; };
inline bool CheckHitHasBeenUsed() { return hitUsedFlag; }; //query to find out if we've used the current hit
inline void SetHitHasBeenUsed() { hitUsedFlag = true; }; //flip the flag to indicate the current hit has been used
inline bool IsEOF() { return eofFlag; }; //see if we've read all available data
inline bool* GetUsedFlagPtr() { return &hitUsedFlag; };
inline void AttachShiftMap(ShiftMap* map) { m_smap = map; };
inline unsigned int GetSize() { return m_size; };
inline unsigned int GetNumberOfHits() { return m_nHits; };
inline bool IsOpen() const { return m_file->is_open(); };
inline CompassHit GetCurrentHit() const { return m_currentHit; }
inline std::string GetName() const { return m_filename; }
inline bool CheckHitHasBeenUsed() const { return hitUsedFlag; } //query to find out if we've used the current hit
inline void SetHitHasBeenUsed() { hitUsedFlag = true; } //flip the flag to indicate the current hit has been used
inline bool IsEOF() const { return eofFlag; } //see if we've read all available data
inline bool* GetUsedFlagPtr() { return &hitUsedFlag; }
inline void AttachShiftMap(ShiftMap* map) { m_smap = map; }
inline unsigned int GetSize() const { return m_size; }
inline unsigned int GetNumberOfHits() const { return m_nHits; }
private:
@ -50,10 +53,12 @@ private:
char* bufferIter;
char* bufferEnd;
ShiftMap* m_smap; //NOT owned by CompassFile. DO NOT delete
bool hitUsedFlag;
int bufsize = 200000; //size of the buffer in hits
int hitsize = 24; //size of a CompassHit in bytes (without alignment padding)
int m_buffersize;
CompassHit m_currentHit;
FilePointer m_file;
bool eofFlag;

17
include/CompassHit.h
View File

@ -1,14 +1,15 @@
#ifndef COMPASS_HIT_H
#define COMPASS_HIT_H
struct CompassHit {
UShort_t board = 400;
UShort_t channel = 400;
ULong64_t timestamp = 0;
UShort_t lgate = 0;
UShort_t sgate = 0;
UInt_t flags = 0;
UInt_t Ns = 0;
struct CompassHit
{
uint16_t board = 400;
uint16_t channel = 400;
uint64_t timestamp = 0;
uint16_t lgate = 0;
uint16_t sgate = 0;
uint32_t flags = 0;
uint32_t Ns = 0;
};
#endif

14
include/CompassRun.h
View File

@ -18,15 +18,17 @@
#include <TGProgressBar.h>
#include <TSystem.h>
class CompassRun {
class CompassRun
{
public:
CompassRun();
CompassRun(const std::string& dir);
~CompassRun();
inline void SetDirectory(const std::string& dir) { directory = dir; };
inline void SetScalerInput(const std::string& filename) { m_scalerinput = filename; };
inline void SetRunNumber(int n) { runNum = n; };
inline void SetShiftMap(const std::string& filename) { m_smap.SetFile(filename); };
inline void SetDirectory(const std::string& dir) { directory = dir; }
inline void SetScalerInput(const std::string& filename) { m_scalerinput = filename; }
inline void SetRunNumber(int n) { runNum = n; }
inline void SetShiftMap(const std::string& filename) { m_smap.SetFile(filename); }
void Convert2RawRoot(const std::string& name);
void Convert2SortedRoot(const std::string& name, const std::string& mapfile, double window);
void Convert2FastSortedRoot(const std::string& name, const std::string& mapfile, double window, double fsi_window, double fic_window);
@ -35,7 +37,7 @@ public:
void Convert2FastAnalyzedRoot(const std::string& name, const std::string& mapfile, double window, double fsi_window, double fic_window,
int zt, int at, int zp, int ap, int ze, int ae, double bke, double b, double theta);
inline void AttachProgressBar(TGProgressBar* pb) { m_pb = pb; };
inline void AttachProgressBar(TGProgressBar* pb) { m_pb = pb; }
private:
bool GetBinaryFiles();

4
include/CutHandler.h
View File

@ -9,9 +9,9 @@ public:
CutHandler(const std::string& filename);
~CutHandler();
void SetCuts(const std::string& filename);
bool IsValid() { return validFlag; };
bool IsValid() { return validFlag; }
bool IsInside(ProcessedEvent* eaddress);
std::vector<TCutG*> GetCuts() { return cut_array; };
std::vector<TCutG*> GetCuts() { return cut_array; }
private:
void InitVariableMap();

83
include/DataStructs.h
View File

@ -6,61 +6,66 @@
#ifndef DATA_STRUCTS_H
#define DATA_STRUCTS_H
using namespace std;
struct DPPChannel {
Double_t Timestamp;
Int_t Channel, Board, Energy, EnergyShort;
Int_t Flags;
struct DPPChannel
{
double Timestamp;
int Channel, Board, Energy, EnergyShort;
int Flags;
};
struct DetectorHit {
Double_t Long=-1, Short=-1, Time=-1;
Int_t Ch=-1;
struct DetectorHit
{
double Long=-1, Short=-1, Time=-1;
int Ch=-1;
};
struct SabreDetector {
vector<DetectorHit> rings;
vector<DetectorHit> wedges;
struct SabreDetector
{
std::vector<DetectorHit> rings;
std::vector<DetectorHit> wedges;
};
struct FPDetector {
vector<DetectorHit> delayFL, delayFR, delayBL, delayBR;
vector<DetectorHit> anodeF, anodeB, scintL, scintR, cathode;
vector<DetectorHit> monitor;
struct FPDetector
{
std::vector<DetectorHit> delayFL, delayFR, delayBL, delayBR;
std::vector<DetectorHit> anodeF, anodeB, scintL, scintR, cathode;
std::vector<DetectorHit> monitor;
};
struct CoincEvent {
struct CoincEvent
{
FPDetector focalPlane;
SabreDetector sabreArray[5]; //index = ChannelMap Id# -1
};
struct ProcessedEvent {
Double_t fp1_tdiff = -1e6, fp2_tdiff = -1e6, fp1_tsum = -1, fp2_tsum = -1,
struct ProcessedEvent
{
double fp1_tdiff = -1e6, fp2_tdiff = -1e6, fp1_tsum = -1, fp2_tsum = -1,
fp1_tcheck = -1, fp2_tcheck = -1;
Double_t fp1_y=-1, fp2_y=-1;
Double_t anodeFront = -1, anodeBack = -1, scintRight = -1, scintLeft = -1;
Double_t scintRightShort = -1, scintLeftShort = -1;
Double_t cathode = -1;
Double_t xavg = -1e6, x1 = -1e6, x2 = -1e6;
Double_t theta = -1e6;
Double_t sabreRingE[5] = {-1,-1,-1,-1,-1}, sabreWedgeE[5] = {-1,-1,-1,-1,-1};
Double_t sabreRingChannel[5] = {-1,-1,-1,-1,-1}, sabreWedgeChannel[5] = {-1,-1,-1,-1,-1};
Double_t sabreRingTime[5] = {-1,-1,-1,-1,-1}, sabreWedgeTime[5] = {-1,-1,-1,-1,-1};
double fp1_y=-1, fp2_y=-1;
double anodeFront = -1, anodeBack = -1, scintRight = -1, scintLeft = -1;
double scintRightShort = -1, scintLeftShort = -1;
double cathode = -1;
double xavg = -1e6, x1 = -1e6, x2 = -1e6;
double theta = -1e6;
double sabreRingE[5] = {-1,-1,-1,-1,-1}, sabreWedgeE[5] = {-1,-1,-1,-1,-1};
double sabreRingChannel[5] = {-1,-1,-1,-1,-1}, sabreWedgeChannel[5] = {-1,-1,-1,-1,-1};
double sabreRingTime[5] = {-1,-1,-1,-1,-1}, sabreWedgeTime[5] = {-1,-1,-1,-1,-1};
Double_t delayFrontRightE = -1, delayFrontLeftE = -1;
Double_t delayBackRightE = -1, delayBackLeftE = -1;
Double_t delayFrontRightShort = -1, delayFrontLeftShort = -1;
Double_t delayBackRightShort = -1, delayBackLeftShort = -1;
Double_t anodeFrontTime = -1, anodeBackTime = -1;
Double_t scintRightTime = -1, scintLeftTime = -1;
Double_t delayFrontMaxTime = -1, delayBackMaxTime = -1;
Double_t delayFrontLeftTime = -1, delayFrontRightTime = -1;
Double_t delayBackLeftTime = -1, delayBackRightTime = -1;
Double_t cathodeTime = -1;
double delayFrontRightE = -1, delayFrontLeftE = -1;
double delayBackRightE = -1, delayBackLeftE = -1;
double delayFrontRightShort = -1, delayFrontLeftShort = -1;
double delayBackRightShort = -1, delayBackLeftShort = -1;
double anodeFrontTime = -1, anodeBackTime = -1;
double scintRightTime = -1, scintLeftTime = -1;
double delayFrontMaxTime = -1, delayBackMaxTime = -1;
double delayFrontLeftTime = -1, delayFrontRightTime = -1;
double delayBackLeftTime = -1, delayBackRightTime = -1;
double cathodeTime = -1;
Double_t monitorE = -1, monitorShort = -1;
Double_t monitorTime = -1;
double monitorE = -1, monitorShort = -1;
double monitorTime = -1;
SabreDetector sabreArray[5]; //index = ChannelMap Id# -1

78
include/GWMEventBuilder.h → include/EVBApp.h
View File

@ -1,21 +1,21 @@
/*
GWMEventBuilder.h
EVBApp.h
Class which represents the API of the event building environment. Wraps together the core concepts
of the event builder, from conversion to plotting. Even intended to be able to archive data.
Currently under development.
Written by G.W. McCann Oct. 2020
*/
#ifndef GWMEVENTBUILDER_H
#define GWMEVENTBUILDER_H
#ifndef EVBAPP_H
#define EVBAPP_H
#include "RunCollector.h"
#include <TGProgressBar.h>
class GWMEventBuilder {
class EVBApp {
public:
GWMEventBuilder();
~GWMEventBuilder();
EVBApp();
~EVBApp();
bool ReadConfigFile(const std::string& filename);
void WriteConfigFile(const std::string& filename);
@ -28,9 +28,6 @@ public:
void Convert2SlowAnalyzedRoot();
void Convert2FastAnalyzedRoot();
inline void SetAnalysisType(int type) { m_analysisType = type;};
inline int GetAnalysisType() { return m_analysisType; };
inline void SetRunRange(int rmin, int rmax) { m_rmin = rmin; m_rmax = rmax; };
inline void SetWorkDirectory(const std::string& fullpath) { m_workspace = fullpath; };
inline void SetChannelMap(const std::string& name) { m_mapfile = name; };
@ -42,38 +39,39 @@ public:
inline void SetScalerFile(const std::string& fullpath) { m_scalerfile = fullpath; };
bool SetKinematicParameters(int zt, int at, int zp, int ap, int ze, int ae, double b, double theta, double bke);
inline int GetRunMin() {return m_rmin;};
inline int GetRunMax() {return m_rmax;};
inline std::string GetWorkDirectory() {return m_workspace;};
inline int GetTargetZ() {return m_ZT;};
inline int GetTargetA() {return m_AT;};
inline int GetProjectileZ() {return m_ZP;};
inline int GetProjectileA() {return m_AP;};
inline int GetEjectileZ() {return m_ZE;};
inline int GetEjectileA() {return m_AE;};
inline int GetResidualZ() {return m_ZR;};
inline int GetResidualA() {return m_AR;};
inline double GetBField() {return m_B;};
inline double GetBeamKE() {return m_BKE;};
inline double GetTheta() {return m_Theta;};
inline double GetSlowCoincidenceWindow() { return m_SlowWindow; };
inline double GetFastWindowIonChamber() { return m_FastWindowIonCh; };
inline double GetFastWindowSABRE() { return m_FastWindowSABRE; };
inline std::string GetChannelMap() { return m_mapfile; };
inline std::string GetBoardShiftFile() { return m_shiftfile; };
inline std::string GetCutList() { return m_cutList; };
inline std::string GetScalerFile() { return m_scalerfile; };
inline int GetRunMin() const {return m_rmin;}
inline int GetRunMax() const {return m_rmax;}
inline std::string GetWorkDirectory() const {return m_workspace;}
inline int GetTargetZ() const {return m_ZT;}
inline int GetTargetA() const {return m_AT;}
inline int GetProjectileZ() const {return m_ZP;}
inline int GetProjectileA() const {return m_AP;}
inline int GetEjectileZ() const {return m_ZE;}
inline int GetEjectileA() const {return m_AE;}
inline int GetResidualZ() const {return m_ZR;}
inline int GetResidualA() const {return m_AR;}
inline double GetBField() const {return m_B;}
inline double GetBeamKE() const {return m_BKE;}
inline double GetTheta() const {return m_Theta;}
inline double GetSlowCoincidenceWindow() const { return m_SlowWindow; }
inline double GetFastWindowIonChamber() const { return m_FastWindowIonCh; }
inline double GetFastWindowSABRE() const { return m_FastWindowSABRE; }
inline std::string GetChannelMap() const { return m_mapfile; }
inline std::string GetBoardShiftFile() const { return m_shiftfile; }
inline std::string GetCutList() const { return m_cutList; }
inline std::string GetScalerFile() const { return m_scalerfile; }
inline void AttachProgressBar(TGProgressBar* pb) { m_pb = pb; };
inline void AttachProgressBar(TGProgressBar* pb) { m_pb = pb; }
enum BuildType {
CONVERT,
CONVERT_S,
CONVERT_SA,
CONVERT_F,
CONVERT_FA,
MERGE,
PLOT
enum Operation
{
Convert,
ConvertSlow,
ConvertSlowA,
ConvertFast,
ConvertFastA,
Merge,
Plot
};
private:
@ -91,8 +89,6 @@ private:
double m_FastWindowIonCh;
double m_FastWindowSABRE;
int m_analysisType;
RunCollector grabber;
TGProgressBar* m_pb;

31
include/FastSort.h
View File

@ -9,24 +9,25 @@
#include "DataStructs.h"
#include <TH2.h>
class FastSort {
class FastSort
{
public:
FastSort(float si_windowSize, float ion_windowSize);
~FastSort();
std::vector<CoincEvent> GetFastEvents(CoincEvent& event);
public:
FastSort(float si_windowSize, float ion_windowSize);
~FastSort();
std::vector<CoincEvent> GetFastEvents(CoincEvent& event);
private:
void ResetSABRE();
void ResetFocalPlane();
void ProcessSABRE(unsigned int scint_index);
void ProcessFocalPlane(unsigned int scint_index, unsigned int ionch_index);
private:
void ResetSABRE();
void ResetFocalPlane();
void ProcessSABRE(unsigned int scint_index);
void ProcessFocalPlane(unsigned int scint_index, unsigned int ionch_index);
float si_coincWindow, ion_coincWindow;
CoincEvent *event_address, slowEvent;
CoincEvent fastEvent, blank;
SabreDetector sblank;
FPDetector fpblank;
float si_coincWindow, ion_coincWindow;
CoincEvent *event_address, slowEvent;
CoincEvent fastEvent, blank;
SabreDetector sblank;
FPDetector fpblank;
};

40
include/FlagHandler.h
View File

@ -3,7 +3,8 @@
#include <map>
struct FlagCount {
struct FlagCount
{
long total_counts=0;
long dead_time=0;
long time_roll=0;
@ -24,30 +25,31 @@ struct FlagCount {
long adc_shutdown=0;
};
class FlagHandler {
class FlagHandler
{
public:
FlagHandler();
FlagHandler(const std::string& filename);
~FlagHandler();
void CheckFlag(int board, int channel, int flag);
const int DEAD_TIME = 0x00000001;
const int TIME_ROLLOVER = 0x00000002;
const int TIME_RESET = 0x00000004;
const int FAKE_EVENT = 0x00000008;
const int MEM_FULL = 0x00000010;
const int TRIG_LOST = 0x00000020;
const int N_TRIG_LOST = 0x00000040;
const int SATURATING_IN_GATE = 0x00000080;
const int TRIG_1024_COUNTED = 0x00000100;
const int SATURATING_INPUT = 0x00000400;
const int N_TRIG_COUNTED = 0x00000800;
const int EVENT_NOT_MATCHED = 0x00001000;
const int FINE_TIME = 0x00004000;
const int PILE_UP = 0x00008000;
const int PLL_LOCK_LOSS = 0x00080000;
const int OVER_TEMP = 0x00100000;
const int ADC_SHUTDOWN = 0x00200000;
const int DeadTime = 0x00000001;
const int TimeRollover = 0x00000002;
const int TimeReset = 0x00000004;
const int FakeEvent = 0x00000008;
const int MemFull = 0x00000010;
const int TrigLost = 0x00000020;
const int NTrigLost = 0x00000040;
const int SaturatingInGate = 0x00000080;
const int Trig1024Counted = 0x00000100;
const int SaturatingInput = 0x00000400;
const int NTrigCounted = 0x00000800;
const int EventNotMatched = 0x00001000;
const int FineTime = 0x00004000;
const int PileUp = 0x00008000;
const int PLLLockLoss = 0x00080000;
const int OverTemp = 0x00100000;
const int ADCShutdown = 0x00200000;
private:
std::ofstream log;

27
include/MassLookup.h
View File

@ -11,22 +11,23 @@ Written by G.W. McCann Aug. 2020
#ifndef MASS_LOOKUP_H
#define MASS_LOOKUP_H
class MassLookup {
class MassLookup
{
public:
MassLookup();
~MassLookup();
double FindMass(int Z, int A);
std::string FindSymbol(int Z, int A);
public:
MassLookup();
~MassLookup();
double FindMass(int Z, int A);
std::string FindSymbol(int Z, int A);
private:
std::unordered_map<std::string, double> massTable;
std::unordered_map<int, std::string> elementTable;
private:
std::unordered_map<std::string, double> massTable;
std::unordered_map<int, std::string> elementTable;
//constants
static constexpr double u_to_mev = 931.4940954;
static constexpr double electron_mass = 0.000548579909;
//constants
static constexpr double u_to_mev = 931.4940954;
static constexpr double electron_mass = 0.000548579909;
};
//static instance for use throught program

3
include/OrderChecker.h
View File

@ -8,7 +8,8 @@
#ifndef ORDERCHECKER_H
#define ORDERCHECKER_H
class OrderChecker {
class OrderChecker
{
public:
OrderChecker();
~OrderChecker();

62
include/RunCollector.h
View File

@ -10,41 +10,35 @@
#ifndef RUNCOLLECTOR_H
#define RUNCOLLECTOR_H
#include <TSystemDirectory.h>
#include <TSystemFile.h>
#include <TCollection.h>
#include <TList.h>
#include <cstdlib>
#include <cstdio>
class RunCollector
{
public:
RunCollector();
RunCollector(const std::string& dirname, const std::string& prefix, const std::string& suffix);
RunCollector(const std::string& dirname, const std::string& prefix, const std::string& suffix, int min, int max);
~RunCollector();
void SetSearchParams(const std::string& dirname, const std::string& prefix, const std::string& suffix, int min, int max);
bool Merge_hadd(const std::string& outname);
bool Merge_TChain(const std::string& outname);
bool GrabAllFiles();
bool GrabFilesInRange();
std::string GrabFile(int runNum);
inline std::string GetSearchDir() { return directory; }
inline std::string GetSearchPrefix() { return m_prefix; }
inline std::string GetSearchSuffix() { return m_suffix; }
inline int GetRunMin() { return m_minRun; }
inline int GetRunMax() { return m_maxRun; }
inline const std::vector<std::string>& GetFileList() { return m_filelist; }
using namespace std;
class RunCollector {
public:
RunCollector();
RunCollector(const string& dirname, const string& prefix, const string& suffix);
RunCollector(const string& dirname, const string& prefix, const string& suffix, int min, int max);
~RunCollector();
void SetSearchParams(const string& dirname, const string& prefix, const string& suffix, int min, int max);
int Merge_hadd(const string& outname);
int Merge_TChain(const string& outname);
int GrabAllFiles();
int GrabFilesInRange();
std::string GrabFile(int runNum);
inline std::string GetSearchDir() {return dir.Data();};
inline std::string GetSearchPrefix() {return run.Data();};
inline std::string GetSearchSuffix() {return end.Data();};
inline int GetRunMin() {return MinRun;};
inline int GetRunMax() {return MaxRun;};
vector<TString> filelist;
private:
bool initFlag;
TString dir;
TString run;
TString end;
Int_t MaxRun, MinRun; //user run limits
const Int_t LITERALMAX = 1000; //class run limit
private:
bool m_initFlag;
std::string m_directory;
std::string m_prefix;
std::string m_suffix;
int m_minRun, m_maxRun; //user run limits
const int m_maxAllowedRuns = 1000; //class run limit
std::vector<std::string> m_filelist;
};
#endif

43
include/SFPAnalyzer.h
View File

@ -13,34 +13,33 @@
#include "DataStructs.h"
#include "FP_kinematics.h"
using namespace std;
class SFPAnalyzer {
class SFPAnalyzer
{
public:
SFPAnalyzer(int zt, int at, int zp, int ap, int ze, int ae, double ep, double angle,
double b);
~SFPAnalyzer();
ProcessedEvent GetProcessedEvent(CoincEvent& event);
inline void ClearHashTable() { rootObj->Clear(); }
inline THashTable* GetHashTable() { return rootObj; }
public:
SFPAnalyzer(int zt, int at, int zp, int ap, int ze, int ae, double ep, double angle,
double b);
~SFPAnalyzer();
ProcessedEvent GetProcessedEvent(CoincEvent& event);
inline void ClearHashTable() { rootObj->Clear(); };
inline THashTable* GetHashTable() { return rootObj; };
private:
void Reset(); //Sets ouput structure back to "zero"
void GetWeights(); //weights for xavg
void AnalyzeEvent(CoincEvent& event);
private:
void Reset(); //Sets ouput structure back to "zero"
void GetWeights(); //weights for xavg
void AnalyzeEvent(CoincEvent& event);
/*Fill wrappers for use with THashTable*/
void MyFill(const std::string& name, int binsx, double minx, double maxx, double valuex,
int binsy, double miny, double maxy, double valuey);
void MyFill(const std::string& name, int binsx, double minx, double maxx, double valuex);
/*Fill wrappers for use with THashTable*/
void MyFill(const string& name, int binsx, double minx, double maxx, double valuex,
int binsy, double miny, double maxy, double valuey);
void MyFill(const string& name, int binsx, double minx, double maxx, double valuex);
CoincEvent *event_address; //Input branch address
ProcessedEvent pevent, blank; //output branch and reset
CoincEvent *event_address; //Input branch address
ProcessedEvent pevent, blank; //output branch and reset
double w1, w2, zfp;
Double_t w1, w2, zfp;
THashTable *rootObj; //root storage
THashTable *rootObj; //root storage
};
#endif

51
include/SFPPlotter.h
View File

@ -13,39 +13,32 @@
#include "CutHandler.h"
#include <TGProgressBar.h>
using namespace std;
class SFPPlotter
{
public:
SFPPlotter();
~SFPPlotter();
inline void AttachProgressBar(TGProgressBar* pb) { m_pb = pb; }
inline void ApplyCutlist(const std::string& listname) { cutter.SetCuts(listname); }
void Run(const std::vector<std::string>& files, const string& output);
class SFPPlotter {
public:
SFPPlotter();
SFPPlotter(bool tf);
~SFPPlotter();
inline void AttachProgressBar(TGProgressBar* pb) { m_pb = pb; };
void ApplyCutlist(const string& listname);
void Run(vector<TString> files, const string& output);
private:
void SetProgressBar(long total);
void Chain(const std::vector<std::string>& files); //Form TChain
void MakeUncutHistograms(const ProcessedEvent& ev, THashTable* table);
void MakeCutHistograms(const ProcessedEvent& ev, THashTable* table);
private:
void SetProgressBar(long total);
void Chain(vector<TString> files); //Form TChain
void MakeUncutHistograms(ProcessedEvent ev);
void MakeCutHistograms(ProcessedEvent ev);
/*Histogram fill wrapper functions*/
void MyFill(THashTable* table, const std::string& name, int binsx, double minx, double maxx, double valuex,
int binsy, double miny, double maxy, double valuey);
void MyFill(THashTable* table, const std::string& name, int binsx, double minx, double maxx, double valuex);
/*Histogram fill wrapper functions*/
void MyFill(const string& name, int binsx, double minx, double maxx, double valuex,
int binsy, double miny, double maxy, double valuey);
void MyFill(const string& name, int binsx, double minx, double maxx, double valuex);
ProcessedEvent *event_address;
ProcessedEvent *event_address;
/*ROOT Storage*/
THashTable *rootObj;
/*Cuts*/
CutHandler cutter;
TChain *chain;
TGProgressBar* m_pb; //GUI progress
/*Cuts*/
CutHandler cutter;
TGProgressBar* m_pb; //GUI progress
};

11
include/ShiftMap.h
View File

@ -12,23 +12,24 @@
#ifndef SHIFTMAP_H
#define SHIFTMAP_H
class ShiftMap {
class ShiftMap
{
public:
ShiftMap();
ShiftMap(const std::string& filename);
~ShiftMap();
void SetFile(const std::string& filename);
inline bool IsSet() { return is_set; };
inline std::string GetFilename() { return m_filename; };
inline bool IsValid() { return m_validFlag; }
inline std::string GetFilename() { return m_filename; }
Long64_t GetShift(int gchan);
private:
void ParseFile();
std::string m_filename;
bool is_set;
bool m_validFlag;
std::unordered_map<int, Long64_t> m_map;
std::unordered_map<int, uint64_t> m_map;
};

4
include/Stopwatch.h
View File

@ -10,8 +10,8 @@
#include <chrono>
class Stopwatch {
class Stopwatch
{
public:
Stopwatch();
~Stopwatch();

44
src/evb/ChannelMap.cpp
View File

@ -22,9 +22,11 @@ ChannelMap::ChannelMap(const std::string& name) :
ChannelMap::~ChannelMap() {}
bool ChannelMap::FillMap(const std::string& name) {
bool ChannelMap::FillMap(const std::string& name)
{
std::ifstream input(name);
if(!input.is_open()) {
if(!input.is_open())
{
is_valid = false;
return is_valid;
}
@ -34,33 +36,39 @@ bool ChannelMap::FillMap(const std::string& name) {
std::getline(input, junk);
std::getline(input, junk);
Channel this_chan;
while(input>>gchan) {
while(input>>gchan)
{
//Set default values
this_chan.detectorType = -1;
this_chan.detectorID = -1;
this_chan.detectorPart = -1;
input>>id>>type>>partname;
if(type == "SABRERING") {
this_chan.detectorType = SABRERING;
if(type == "SABRERING")
{
this_chan.detectorType = DetAttribute::SabreRing;
this_chan.detectorID = id;
this_chan.detectorPart = std::stoi(partname);
} else if(type == "SABREWEDGE") {
this_chan.detectorType = SABREWEDGE;
}
else if(type == "SABREWEDGE")
{
this_chan.detectorType = DetAttribute::SabreWedge;
this_chan.detectorID = id;
this_chan.detectorPart = std::stoi(partname);
} else if (type == "FOCALPLANE") {
}
else if (type == "FOCALPLANE")
{
this_chan.detectorType = FOCALPLANE;
this_chan.detectorID = id;
if(partname == "SCINTRIGHT") this_chan.detectorPart = SCINTRIGHT;
else if(partname == "SCINTLEFT") this_chan.detectorPart = SCINTLEFT;
else if(partname == "DELAYFR") this_chan.detectorPart = DELAYFR;
else if(partname == "DELAYFL") this_chan.detectorPart = DELAYFL;
else if(partname == "DELAYBR") this_chan.detectorPart = DELAYBR;
else if(partname == "DELAYBL") this_chan.detectorPart = DELAYBL;
else if(partname == "CATHODE") this_chan.detectorPart = CATHODE;
else if(partname == "ANODEFRONT") this_chan.detectorPart = ANODEFRONT;
else if(partname == "ANODEBACK") this_chan.detectorPart = ANODEBACK;
else if(partname == "MONITOR") this_chan.detectorPart = MONITOR;
if(partname == "SCINTRIGHT") this_chan.detectorPart = DetAttribute::ScintRight;
else if(partname == "SCINTLEFT") this_chan.detectorPart = DetAttribute::ScintLeft;
else if(partname == "DELAYFR") this_chan.detectorPart = DetAttribute::DelayFR;
else if(partname == "DELAYFL") this_chan.detectorPart = DetAttribute::DelayFL;
else if(partname == "DELAYBR") this_chan.detectorPart = DetAttribute::DelayBR;
else if(partname == "DELAYBL") this_chan.detectorPart = DetAttribute::DelayBL;
else if(partname == "CATHODE") this_chan.detectorPart = DetAttribute::Cathode;
else if(partname == "ANODEFRONT") this_chan.detectorPart = DetAttribute::AnodeFront;
else if(partname == "ANODEBACK") this_chan.detectorPart = DetAttribute::AnodeBack;
else if(partname == "MONITOR") this_chan.detectorPart = DetAttribute::Monitor;
}
cmap[gchan] = this_chan;

60
src/evb/CompassFile.cpp
View File

@ -35,11 +35,13 @@ CompassFile::CompassFile(const std::string& filename, int bsize) :
Open(filename);
}
CompassFile::~CompassFile() {
CompassFile::~CompassFile()
{
Close();
}
void CompassFile::Open(const std::string& filename) {
void CompassFile::Open(const std::string& filename)
{
eofFlag = false;
hitUsedFlag = true;
m_filename = filename;
@ -48,21 +50,28 @@ void CompassFile::Open(const std::string& filename) {
m_file->seekg(0, std::ios_base::end);
m_size = m_file->tellg();
m_nHits = m_size/24;
if(m_size == 0) {
if(m_size == 0)
{
eofFlag = true;
} else {
}
else
{
m_file->seekg(0, std::ios_base::beg);
}
}
void CompassFile::Close() {
if(IsOpen()) {
void CompassFile::Close()
{
if(IsOpen())
{
m_file->close();
}
}
int CompassFile::GetHitSize() {
if(!IsOpen()) {
int CompassFile::GetHitSize()
{
if(!IsOpen())
{
std::cerr<<"Unable to get hit size due to file not being open!"<<std::endl;
return 0;
}
@ -90,14 +99,17 @@ int CompassFile::GetHitSize() {
If the file cannot be opened, signals as though file is EOF
*/
bool CompassFile::GetNextHit() {
bool CompassFile::GetNextHit()
{
if(!IsOpen()) return true;
if((bufferIter == nullptr || bufferIter == bufferEnd) && !IsEOF()) {
if((bufferIter == nullptr || bufferIter == bufferEnd) && !IsEOF())
{
GetNextBuffer();
}
if(!IsEOF()) {
if(!IsEOF())
{
ParseNextHit();
hitUsedFlag = false;
}
@ -112,9 +124,11 @@ bool CompassFile::GetNextHit() {
bit upon pulling the last buffer, but this class waits until that entire
last buffer is read to singal EOF (the true end of file).
*/
void CompassFile::GetNextBuffer() {
void CompassFile::GetNextBuffer()
{
if(m_file->eof()) {
if(m_file->eof())
{
eofFlag = true;
return;
}
@ -126,24 +140,26 @@ void CompassFile::GetNextBuffer() {
}
void CompassFile::ParseNextHit() {
void CompassFile::ParseNextHit()
{
m_currentHit.board = *((UShort_t*)bufferIter);
m_currentHit.board = *((uint16_t*)bufferIter);
bufferIter += 2;
m_currentHit.channel = *((UShort_t*)bufferIter);
m_currentHit.channel = *((uint16_t*)bufferIter);
bufferIter += 2;
m_currentHit.timestamp = *((ULong64_t*)bufferIter);
m_currentHit.timestamp = *((uint64_t*)bufferIter);
bufferIter += 8;
m_currentHit.lgate = *((UShort_t*)bufferIter);
m_currentHit.lgate = *((uint16_t*)bufferIter);
bufferIter += 2;
m_currentHit.sgate = *((UShort_t*)bufferIter);
m_currentHit.sgate = *((uint16_t*)bufferIter);
bufferIter += 2;
m_currentHit.flags = *((UInt_t*)bufferIter);
m_currentHit.flags = *((uint32_t*)bufferIter);
bufferIter += 4;
m_currentHit.Ns = *((UInt_t*)bufferIter);
m_currentHit.Ns = *((uint32_t*)bufferIter);
bufferIter += 4;
if(m_smap != nullptr) { //memory safety
if(m_smap != nullptr)
{ //memory safety
int gchan = m_currentHit.channel + m_currentHit.board*16;
m_currentHit.timestamp += m_smap->GetShift(gchan);
}

292
src/evb/CompassRun.cpp
View File

@ -33,9 +33,11 @@ CompassRun::~CompassRun() {}
/*Load em into a map*/
void CompassRun::SetScalers() {
void CompassRun::SetScalers()
{
std::ifstream input(m_scalerinput);
if(!input.is_open()) return;
if(!input.is_open())
return;
m_scaler_flag = true;
std::string junk, filename, varname;
@ -43,7 +45,8 @@ void CompassRun::SetScalers() {
std::getline(input, junk);
std::getline(input, junk);
m_scaler_map.clear();
while(input>>filename) {
while(input>>filename)
{
input>>varname;
filename = directory+filename+"_run_"+to_string(runNum)+".bin";
m_scaler_map[filename] = TParameter<Long64_t>(varname.c_str(), init);
@ -51,7 +54,8 @@ void CompassRun::SetScalers() {
input.close();
}
bool CompassRun::GetBinaryFiles() {
bool CompassRun::GetBinaryFiles()
{
std::string prefix = "";
std::string suffix = ".bin"; //binaries
RunCollector grabber(directory, prefix, suffix);
@ -62,26 +66,31 @@ bool CompassRun::GetBinaryFiles() {
bool scalerd;
m_totalHits = 0; //reset total run size
for(auto& entry : grabber.filelist) {
for(auto& entry : grabber.filelist)
{
//Handle scaler files, if they exist
if(m_scaler_flag) {
if(m_scaler_flag)
{
scalerd = false;
for(auto& scaler_pair : m_scaler_map) {
if(std::string(entry.Data()) == scaler_pair.first) {
for(auto& scaler_pair : m_scaler_map)
{
if(std::string(entry.Data()) == scaler_pair.first)
{
ReadScalerData(entry.Data());
scalerd = true;
break;
}
}
if(scalerd) continue;
if(scalerd)
continue;
}
m_datafiles.emplace_back(entry.Data());
m_datafiles[m_datafiles.size()-1].AttachShiftMap(&m_smap);
//Any time we have a file that fails to be found, we terminate the whole process
if(!m_datafiles[m_datafiles.size() - 1].IsOpen()) {
if(!m_datafiles[m_datafiles.size() - 1].IsOpen())
return false;
}
m_totalHits += m_datafiles[m_datafiles.size()-1].GetNumberOfHits();
}
@ -92,16 +101,20 @@ bool CompassRun::GetBinaryFiles() {
Pure counting of scalers. Potential upgrade path to something like
average count rate etc.
*/
void CompassRun::ReadScalerData(const std::string& filename) {
if(!m_scaler_flag) return;
void CompassRun::ReadScalerData(const std::string& filename)
{
if(!m_scaler_flag)
return;
Long64_t count;
count = 0;
CompassFile file(filename);
auto& this_param = m_scaler_map[file.GetName()];
while(true) {
while(true)
{
file.GetNextHit();
if(file.IsEOF()) break;
if(file.IsEOF())
break;
count++;
}
this_param.SetVal(count);
@ -116,25 +129,34 @@ void CompassRun::ReadScalerData(const std::string& filename) {
of a rolling start index. Once a file has gone EOF, we no longer need it. If this is the first file in the list, we can just skip
that index all together. In this way, the loop can go from N times to N-1 times.
*/
bool CompassRun::GetHitsFromFiles() {
bool CompassRun::GetHitsFromFiles()
{
std::pair<CompassHit, bool*> earliestHit = make_pair(CompassHit(), nullptr);
for(unsigned int i=startIndex; i<m_datafiles.size(); i++) {
if(m_datafiles[i].CheckHitHasBeenUsed()) {
for(unsigned int i=startIndex; i<m_datafiles.size(); i++)
{
if(m_datafiles[i].CheckHitHasBeenUsed())
{
m_datafiles[i].GetNextHit();
}
if(m_datafiles[i].IsEOF()) {
if(m_datafiles[i].IsEOF())
{
if(i == startIndex) startIndex++;
continue;
} else if(i == startIndex) {
}
else if(i == startIndex)
{
earliestHit = make_pair(m_datafiles[i].GetCurrentHit(), m_datafiles[i].GetUsedFlagPtr());
} else if(m_datafiles[i].GetCurrentHit().timestamp < earliestHit.first.timestamp) {
}
else if(m_datafiles[i].GetCurrentHit().timestamp < earliestHit.first.timestamp)
{
earliestHit = make_pair(m_datafiles[i].GetCurrentHit(), m_datafiles[i].GetUsedFlagPtr());
}
}
if(earliestHit.second == nullptr) return false; //Make sure that there actually was a hit
if(earliestHit.second == nullptr)
return false; //Make sure that there actually was a hit
hit = earliestHit.first;
*earliestHit.second = true;
return true;
@ -151,96 +173,118 @@ void CompassRun::Convert2RawRoot(const std::string& name) {
outtree->Branch("Timestamp", &hit.timestamp);
outtree->Branch("Flags", &hit.flags);
if(!m_smap.IsSet()) {
if(!m_smap.IsSet())
{
std::cerr<<"Bad shift map at CompassRun::Convert()."<<std::endl;
std::cerr<<"Shifts will be locked to 0"<<std::endl;
}
SetScalers();
if(!GetBinaryFiles()) {
if(!GetBinaryFiles())
{
std::cerr<<"Unable to open a file!"<<std::endl;
return;
}
if(m_pb) SetProgressBar();
if(m_pb)
SetProgressBar();
startIndex = 0; //Reset the startIndex
unsigned int count = 0, flush = m_totalHits*0.01, flush_count = 0;
if(flush == 0) flush = 1;
while(true) {
if(flush == 0)
flush = 1;
while(true)
{
count++;
if(count == flush) { //Progress Log
if(m_pb) {
if(count == flush)
{ //Progress Log
if(m_pb)
{
m_pb->Increment(count);
gSystem->ProcessEvents();
count=0;
} else {
}
else
{
count = 0;
flush_count++;
std::cout<<"\rPercent of run built: "<<flush_count*10<<"%"<<std::flush;
}
}
if(!GetHitsFromFiles()) break;
if(!GetHitsFromFiles())
break;
outtree->Fill();
}
output->cd();
outtree->Write(outtree->GetName(), TObject::kOverwrite);
for(auto& entry : m_scaler_map) {
for(auto& entry : m_scaler_map)
entry.second.Write();
}
output->Close();
}
void CompassRun::Convert2SortedRoot(const std::string& name, const std::string& mapfile, double window) {
void CompassRun::Convert2SortedRoot(const std::string& name, const std::string& mapfile, double window)
{
TFile* output = TFile::Open(name.c_str(), "RECREATE");
TTree* outtree = new TTree("SortTree", "SortTree");
outtree->Branch("event", &event);
if(!m_smap.IsSet()) {
if(!m_smap.IsSet())
{
std::cerr<<"Bad shift map at CompassRun::Convert()."<<std::endl;
std::cerr<<"Shifts will be locked to 0"<<std::endl;
}
SetScalers();
if(!GetBinaryFiles()) {
if(!GetBinaryFiles())
{
std::cerr<<"Unable to open a file!"<<std::endl;
return;
}
if(m_pb) SetProgressBar();
if(m_pb)
SetProgressBar();
startIndex = 0;
SlowSort coincidizer(window, mapfile);
bool killFlag = false;
unsigned int count = 0, flush = m_totalHits*0.01, flush_count = 0;
if(flush == 0) flush = 1;
while(true) {
if(flush == 0)
flush = 1;
while(true)
{
count++;
if(count == flush) {
if(m_pb) {
if(count == flush)
{
if(m_pb)
{
m_pb->Increment(count);
gSystem->ProcessEvents();
count=0;
} else {
}
else
{
count = 0;
flush_count++;
std::cout<<"\rPercent of run built: "<<flush_count*10<<"%"<<std::flush;
}
}
if(!GetHitsFromFiles()) {
if(!GetHitsFromFiles())
{
coincidizer.FlushHitsToEvent();
killFlag = true;
} else {
}
else
coincidizer.AddHitToEvent(hit);
}
if(coincidizer.IsEventReady()) {
if(coincidizer.IsEventReady())
{
event = coincidizer.GetEvent();
outtree->Fill();
if(killFlag) break;
@ -249,32 +293,36 @@ void CompassRun::Convert2SortedRoot(const std::string& name, const std::string&
output->cd();
outtree->Write(outtree->GetName(), TObject::kOverwrite);
for(auto& entry : m_scaler_map) {
for(auto& entry : m_scaler_map)
entry.second.Write();
}
coincidizer.GetEventStats()->Write();
output->Close();
}
void CompassRun::Convert2FastSortedRoot(const std::string& name, const std::string& mapfile, double window, double fsi_window, double fic_window) {
void CompassRun::Convert2FastSortedRoot(const std::string& name, const std::string& mapfile, double window, double fsi_window, double fic_window)
{
TFile* output = TFile::Open(name.c_str(), "RECREATE");
TTree* outtree = new TTree("SortTree", "SortTree");
outtree->Branch("event", &event);
if(!m_smap.IsSet()) {
if(!m_smap.IsSet())
{
std::cerr<<"Bad shift map at CompassRun::Convert()."<<std::endl;
std::cerr<<"Shifts will be locked to 0"<<std::endl;
}
SetScalers();
if(!GetBinaryFiles()) {
if(!GetBinaryFiles())
{
std::cerr<<"Unable to open a file!"<<std::endl;
return;
}
if(m_pb) SetProgressBar();
if(m_pb)
SetProgressBar();
startIndex = 0;
CoincEvent this_event;
@ -286,72 +334,88 @@ void CompassRun::Convert2FastSortedRoot(const std::string& name, const std::stri
bool killFlag = false;
unsigned int count = 0, flush = m_totalHits*0.01, flush_count = 0;
if(flush == 0) flush = 1;
while(true) {
if(flush == 0)
flush = 1;
while(true)
{
count++;
if(count == flush) {
if(m_pb) {
if(count == flush)
{
if(m_pb)
{
m_pb->Increment(count);
gSystem->ProcessEvents();
count=0;
} else {
}
else
{
count = 0;
flush_count++;
std::cout<<"\rPercent of run built: "<<flush_count*10<<"%"<<std::flush;
}
}
if(!GetHitsFromFiles()) {
if(!GetHitsFromFiles())
{
coincidizer.FlushHitsToEvent();
killFlag = true;
} else {
}
else
{
flagger.CheckFlag(hit.board, hit.channel, hit.flags);
coincidizer.AddHitToEvent(hit);
}
if(coincidizer.IsEventReady()) {
if(coincidizer.IsEventReady())
{
this_event = coincidizer.GetEvent();
fast_events = speedyCoincidizer.GetFastEvents(this_event);
for(auto& entry : fast_events) {
for(auto& entry : fast_events)
{
event = entry;
outtree->Fill();
}
if(killFlag) break;
if(killFlag)
break;
}
}
output->cd();
outtree->Write(outtree->GetName(), TObject::kOverwrite);
for(auto& entry : m_scaler_map) {
for(auto& entry : m_scaler_map)
entry.second.Write();
}
coincidizer.GetEventStats()->Write();
output->Close();
}
void CompassRun::Convert2SlowAnalyzedRoot(const std::string& name, const std::string& mapfile, double window,
int zt, int at, int zp, int ap, int ze, int ae, double bke, double b, double theta) {
int zt, int at, int zp, int ap, int ze, int ae, double bke, double b, double theta)
{
TFile* output = TFile::Open(name.c_str(), "RECREATE");
TTree* outtree = new TTree("SPSTree", "SPSTree");
outtree->Branch("event", &pevent);
if(!m_smap.IsSet()) {
if(!m_smap.IsSet())
{
std::cerr<<"Bad shift map at CompassRun::Convert()."<<std::endl;
std::cerr<<"Shifts will be locked to 0"<<std::endl;
}
SetScalers();
if(!GetBinaryFiles()) {
if(!GetBinaryFiles())
{
std::cerr<<"Unable to open a file!"<<std::endl;
return;
}
if(m_pb) SetProgressBar();
if(m_pb)
SetProgressBar();
startIndex = 0;
CoincEvent this_event;
@ -372,44 +436,55 @@ void CompassRun::Convert2SlowAnalyzedRoot(const std::string& name, const std::st
bool killFlag = false;
unsigned int count = 0, flush = m_totalHits*0.01, flush_count = 0;
if(flush == 0) flush = 1;
while(true) {
if(flush == 0)
flush = 1;
while(true)
{
count++;
if(count == flush) {
if(m_pb) {
if(count == flush)
{
if(m_pb)
{
m_pb->Increment(count);
gSystem->ProcessEvents();
count=0;
} else {
}
else
{
count = 0;
flush_count++;
std::cout<<"\rPercent of run built: "<<flush_count*10<<"%"<<std::flush;
}
}
if(!GetHitsFromFiles()) {
if(!GetHitsFromFiles())
{
coincidizer.FlushHitsToEvent();
killFlag = true;
} else {
}
else
{
coincidizer.AddHitToEvent(hit);
}
if(coincidizer.IsEventReady()) {
if(coincidizer.IsEventReady())
{
this_event = coincidizer.GetEvent();
pevent = analyzer.GetProcessedEvent(this_event);
outtree->Fill();
if(killFlag) break;
if(killFlag)
break;
}
}
output->cd();
outtree->Write(outtree->GetName(), TObject::kOverwrite);
for(auto& entry : m_scaler_map) {
for(auto& entry : m_scaler_map)
entry.second.Write();
}
for(auto& entry : parvec) {
for(auto& entry : parvec)
entry.Write();
}
coincidizer.GetEventStats()->Write();
analyzer.GetHashTable()->Write();
analyzer.ClearHashTable();
@ -417,26 +492,30 @@ void CompassRun::Convert2SlowAnalyzedRoot(const std::string& name, const std::st
}
void CompassRun::Convert2FastAnalyzedRoot(const std::string& name, const std::string& mapfile, double window, double fsi_window, double fic_window,
int zt, int at, int zp, int ap, int ze, int ae, double bke, double b, double theta) {
int zt, int at, int zp, int ap, int ze, int ae, double bke, double b, double theta)
{
TFile* output = TFile::Open(name.c_str(), "RECREATE");
TTree* outtree = new TTree("SPSTree", "SPSTree");
outtree->Branch("event", &pevent);
if(!m_smap.IsSet()) {
if(!m_smap.IsSet())
{
std::cerr<<"Bad shift map at CompassRun::Convert()."<<std::endl;
std::cerr<<"Shifts will be locked to 0"<<std::endl;
}
SetScalers();
if(!GetBinaryFiles()) {
if(!GetBinaryFiles())
{
std::cerr<<"Unable to open a file!"<<std::endl;
return;
}
if(m_pb) SetProgressBar();
if(m_pb)
SetProgressBar();
startIndex = 0;
CoincEvent this_event;
@ -461,56 +540,69 @@ void CompassRun::Convert2FastAnalyzedRoot(const std::string& name, const std::st
bool killFlag = false;
unsigned int count = 0, flush = m_totalHits*0.01, flush_count = 0;
if(flush == 0) flush = 1;
while(true) {
if(flush == 0)
flush = 1;
while(true)
{
count++;
if(count == flush) {
if(m_pb) {
if(count == flush)
{
if(m_pb)
{
m_pb->Increment(count);
gSystem->ProcessEvents();
count=0;
} else {
}
else
{
count = 0;
flush_count++;
std::cout<<"\rPercent of run built: "<<flush_count*10<<"%"<<std::flush;
}
}
if(!GetHitsFromFiles()) {
if(!GetHitsFromFiles())
{
coincidizer.FlushHitsToEvent();
killFlag = true;
} else {
}
else
{
flagger.CheckFlag(hit.board, hit.channel, hit.flags);
coincidizer.AddHitToEvent(hit);
}
if(coincidizer.IsEventReady()) {
if(coincidizer.IsEventReady())
{
this_event = coincidizer.GetEvent();
fast_events = speedyCoincidizer.GetFastEvents(this_event);
for(auto& entry : fast_events) {
for(auto& entry : fast_events)
{
pevent = analyzer.GetProcessedEvent(entry);
outtree->Fill();
}
if(killFlag) break;
if(killFlag)
break;
}
}
output->cd();
outtree->Write(outtree->GetName(), TObject::kOverwrite);
for(auto& entry : m_scaler_map) {
for(auto& entry : m_scaler_map)
entry.second.Write();
}
for(auto& entry : parvec) {
for(auto& entry : parvec)
entry.Write();
}
coincidizer.GetEventStats()->Write();
analyzer.GetHashTable()->Write();
analyzer.ClearHashTable();
output->Close();
}
void CompassRun::SetProgressBar() {
void CompassRun::SetProgressBar()
{
m_pb->SetMax(m_totalHits);
m_pb->SetMin(0);
m_pb->SetPosition(0);

46
src/evb/CutHandler.cpp
View File

@ -14,16 +14,19 @@ CutHandler::CutHandler(const std::string& filename) :
InitVariableMap();
}
CutHandler::~CutHandler() {
for(unsigned int i=0; i<file_array.size(); i++) {
if(file_array[i]->IsOpen()) file_array[i]->Close();
}
CutHandler::~CutHandler()
{
for(unsigned int i=0; i<file_array.size(); i++)
if(file_array[i]->IsOpen())
file_array[i]->Close();
}
void CutHandler::SetCuts(const std::string& filename) {
void CutHandler::SetCuts(const std::string& filename)
{
std::ifstream cutlist(filename);
if(!cutlist.is_open()) {
if(!cutlist.is_open())
{
validFlag = false;
}
@ -33,17 +36,21 @@ void CutHandler::SetCuts(const std::string& filename) {
cut_array.clear();
file_array.clear();
while(cutlist>>name) {
while(cutlist>>name)
{
cutlist>>fname>>varx>>vary;
TFile* file = TFile::Open(fname.c_str(), "READ");
TCutG* cut = (TCutG*) file->Get("CUTG");
if(cut) {
if(cut)
{
cut->SetVarX(varx.c_str());
cut->SetVarY(vary.c_str());
cut->SetName(name.c_str());
cut_array.push_back(cut);
file_array.push_back(file);
} else {
}
else
{
validFlag = false;
std::cerr<<"CutHandler has encountered a bad cut at file: "<<file<<"."<<std::endl;
std::cerr<<"The file either does not exist or does not contain a TCutG named CUTG"<<std::endl;
@ -52,17 +59,17 @@ void CutHandler::SetCuts(const std::string& filename) {
}
}
if(cut_array.size() > 0) {
if(cut_array.size() > 0)
validFlag = true;
} else {
else
validFlag = false;
}
}
/*
ADD MORE VARIABLES HERE!
*/
void CutHandler::InitVariableMap() {
void CutHandler::InitVariableMap()
{
varmap["x1"] = &m_event.x1;
varmap["x2"] = &m_event.x2;
varmap["xavg"] = &m_event.xavg;
@ -71,23 +78,26 @@ void CutHandler::InitVariableMap() {
varmap["cathode"] = &m_event.cathode;
}
bool CutHandler::IsInside(ProcessedEvent* eaddress) {
bool CutHandler::IsInside(ProcessedEvent* eaddress)
{
m_event = *eaddress;
std::string x, y;
for(unsigned int i=0; i<cut_array.size(); i++) {
for(unsigned int i=0; i<cut_array.size(); i++)
{
TCutG* cut = cut_array[i];
x = cut->GetVarX();
y = cut->GetVarY();
auto xentry = varmap.find(x);
auto yentry = varmap.find(y);
if(xentry == varmap.end() || yentry == varmap.end()) {
if(xentry == varmap.end() || yentry == varmap.end())
{
std::cerr<<"Unmapped variable called in CutHandler::IsInside()! Var names: "<<xentry->first<<" , "<<yentry->first<<std::endl;
return false;
}
if(!cut->IsInside(*(xentry->second), *(yentry->second))) {
if(!cut->IsInside(*(xentry->second), *(yentry->second)))
return false;
}
}
return true;

103
src/evb/GWMEventBuilder.cpp → src/evb/EVBApp.cpp
View File

@ -1,5 +1,5 @@
/*
GWMEventBuilder.cpp
EVBApp.cpp
Class which represents the API of the event building environment. Wraps together the core concepts
of the event builder, from conversion to plotting. Even intended to be able to archive data.
Currently under development.
@ -8,7 +8,7 @@
*/
#include "EventBuilder.h"
#include <cstdlib>
#include "GWMEventBuilder.h"
#include "EVBApp.h"
#include "RunCollector.h"
#include "CompassRun.h"
#include "SlowSort.h"
@ -16,21 +16,23 @@
#include "SFPAnalyzer.h"
#include "SFPPlotter.h"
GWMEventBuilder::GWMEventBuilder() :
EVBApp::EVBApp() :
m_rmin(0), m_rmax(0), m_ZT(0), m_AT(0), m_ZP(0), m_AP(0), m_ZE(0), m_AE(0), m_ZR(0), m_AR(0),
m_B(0), m_Theta(0), m_BKE(0), m_workspace("none"), m_mapfile("none"), m_shiftfile("none"),
m_cutList("none"), m_SlowWindow(0), m_FastWindowIonCh(0), m_FastWindowSABRE(0), m_pb(nullptr)
{
}
GWMEventBuilder::~GWMEventBuilder()
EVBApp::~EVBApp()
{
}
bool GWMEventBuilder::ReadConfigFile(const std::string& fullpath) {
bool EVBApp::ReadConfigFile(const std::string& fullpath)
{
std::cout<<"Reading in configuration from file: "<<fullpath<<std::endl;
std::ifstream input(fullpath);
if(!input.is_open()) {
if(!input.is_open())
{
std::cout<<"Read failed! Unable to open input file!"<<std::endl;
return false;
}
@ -70,11 +72,13 @@ bool GWMEventBuilder::ReadConfigFile(const std::string& fullpath) {
return true;
}
void GWMEventBuilder::WriteConfigFile(const std::string& fullpath) {
void EVBApp::WriteConfigFile(const std::string& fullpath)
{
std::cout<<"Writing out configuration to file: "<<fullpath<<std::endl;
std::ofstream output(fullpath);
if(!output.is_open()) {
if(!output.is_open())
{
std::cout<<"Write failed! Unable to open output file!"<<std::endl;
return;
}
@ -113,7 +117,8 @@ void GWMEventBuilder::WriteConfigFile(const std::string& fullpath) {
}
void GWMEventBuilder::PlotHistograms() {
void EVBApp::PlotHistograms()
{
std::string analyze_dir = m_workspace+"/analyzed/";
std::string plot_file = m_workspace+"/histograms/run_"+to_string(m_rmin)+"_"+to_string(m_rmax)+".root";
SFPPlotter grammer;
@ -125,20 +130,25 @@ void GWMEventBuilder::PlotHistograms() {
std::cout<<"Cut List File: "<<m_cutList<<std::endl;
std::cout<<"Histogram File: "<<plot_file<<std::endl;
if(m_pb) grammer.AttachProgressBar(m_pb);
if(m_pb)
grammer.AttachProgressBar(m_pb);
grabber.SetSearchParams(analyze_dir, "", ".root", m_rmin, m_rmax);
if(grabber.GrabFilesInRange()) {
if(grabber.GrabFilesInRange())
{
std::cout<<"Working...";
grammer.Run(grabber.filelist, plot_file);
std::cout<<" Complete."<<std::endl;
} else {
}
else
{
std::cout<<"Unable to find files at PlotHistograms"<<std::endl;
}
std::cout<<"-------------------------------------------"<<std::endl;
}
void GWMEventBuilder::Convert2RawRoot() {
void EVBApp::Convert2RawRoot()
{
std::string rawroot_dir = m_workspace+"/raw_root/";
std::string unpack_dir = m_workspace+"/temp_binary/";
std::string binary_dir = m_workspace+"/raw_binary/";
@ -160,13 +170,16 @@ void GWMEventBuilder::Convert2RawRoot() {
CompassRun converter(unpack_dir);
converter.SetShiftMap(m_shiftfile);
converter.SetScalerInput(m_scalerfile);
if(m_pb) converter.AttachProgressBar(m_pb);
if(m_pb)
converter.AttachProgressBar(m_pb);
std::cout<<"Beginning conversion..."<<std::endl;
for(int i=m_rmin; i<=m_rmax; i++) {
for(int i=m_rmin; i<=m_rmax; i++)
{
binfile = grabber.GrabFile(i);
if(binfile == "") continue;
if(binfile == "")
continue;
converter.SetRunNumber(i);
std::cout<<"Converting file: "<<binfile<<std::endl;
@ -184,7 +197,8 @@ void GWMEventBuilder::Convert2RawRoot() {
}
void GWMEventBuilder::MergeROOTFiles() {
void EVBApp::MergeROOTFiles()
{
std::string merge_file = m_workspace+"/merged/run_"+to_string(m_rmin)+"_"+to_string(m_rmax)+".root";
std::string file_dir = m_workspace+"/analyzed/";
std::cout<<"-------------GWM Event Builder-------------"<<std::endl;
@ -196,7 +210,8 @@ void GWMEventBuilder::MergeROOTFiles() {
std::string suffix = ".root";
grabber.SetSearchParams(file_dir, prefix, suffix,m_rmin,m_rmax);
std::cout<<"Beginning the merge...";
if(!grabber.Merge_TChain(merge_file)) {
if(!grabber.Merge_TChain(merge_file))
{
std::cout<<"Unable to find files at MergeROOTFiles"<<std::endl;
return;
}
@ -204,7 +219,8 @@ void GWMEventBuilder::MergeROOTFiles() {
std::cout<<"-------------------------------------------"<<std::endl;
}
void GWMEventBuilder::Convert2SortedRoot() {
void EVBApp::Convert2SortedRoot()
{
std::string sortroot_dir = m_workspace+"/sorted/";
std::string unpack_dir = m_workspace+"/temp_binary/";
std::string binary_dir = m_workspace+"/raw_binary/";
@ -227,13 +243,16 @@ void GWMEventBuilder::Convert2SortedRoot() {
CompassRun converter(unpack_dir);
converter.SetShiftMap(m_shiftfile);
converter.SetScalerInput(m_scalerfile);
if(m_pb) converter.AttachProgressBar(m_pb);
if(m_pb)
converter.AttachProgressBar(m_pb);
std::cout<<"Beginning conversion..."<<std::endl;
for(int i=m_rmin; i<= m_rmax; i++) {
for(int i=m_rmin; i<= m_rmax; i++)
{
binfile = grabber.GrabFile(i);
if(binfile == "") continue;
if(binfile == "")
continue;
converter.SetRunNumber(i);
std::cout<<"Converting file: "<<binfile<<std::endl;
@ -250,7 +269,7 @@ void GWMEventBuilder::Convert2SortedRoot() {
std::cout<<"-------------------------------------------"<<std::endl;
}
void GWMEventBuilder::Convert2FastSortedRoot() {
void EVBApp::Convert2FastSortedRoot() {
std::string sortroot_dir = m_workspace+"/fast/";
std::string unpack_dir = m_workspace+"/temp_binary/";
std::string binary_dir = m_workspace+"/raw_binary/";
@ -273,13 +292,16 @@ void GWMEventBuilder::Convert2FastSortedRoot() {
CompassRun converter(unpack_dir);
converter.SetShiftMap(m_shiftfile);
converter.SetScalerInput(m_scalerfile);
if(m_pb) converter.AttachProgressBar(m_pb);
if(m_pb)
converter.AttachProgressBar(m_pb);
std::cout<<"Beginning conversion..."<<std::endl;
for(int i=m_rmin; i<=m_rmax; i++) {
for(int i=m_rmin; i<=m_rmax; i++)
{
binfile = grabber.GrabFile(i);
if(binfile == "") continue;
if(binfile == "")
continue;
converter.SetRunNumber(i);
std::cout<<"Converting file: "<<binfile<<std::endl;
@ -296,7 +318,7 @@ void GWMEventBuilder::Convert2FastSortedRoot() {
std::cout<<"-------------------------------------------"<<std::endl;
}
void GWMEventBuilder::Convert2SlowAnalyzedRoot() {
void EVBApp::Convert2SlowAnalyzedRoot() {
std::string sortroot_dir = m_workspace+"/analyzed/";
std::string unpack_dir = m_workspace+"/temp_binary/";
std::string binary_dir = m_workspace+"/raw_binary/";
@ -318,13 +340,16 @@ void GWMEventBuilder::Convert2SlowAnalyzedRoot() {
CompassRun converter(unpack_dir);
converter.SetShiftMap(m_shiftfile);
converter.SetScalerInput(m_scalerfile);
if(m_pb) converter.AttachProgressBar(m_pb);
if(m_pb)
converter.AttachProgressBar(m_pb);
std::cout<<"Beginning conversion..."<<std::endl;
for(int i=m_rmin; i<=m_rmax; i++) {
for(int i=m_rmin; i<=m_rmax; i++)
{
binfile = grabber.GrabFile(i);
if(binfile == "") continue;
if(binfile == "")
continue;
converter.SetRunNumber(i);
std::cout<<"Converting file: "<<binfile<<std::endl;
@ -341,7 +366,8 @@ void GWMEventBuilder::Convert2SlowAnalyzedRoot() {
std::cout<<"-------------------------------------------"<<std::endl;
}
void GWMEventBuilder::Convert2FastAnalyzedRoot() {
void EVBApp::Convert2FastAnalyzedRoot()
{
std::string sortroot_dir = m_workspace+"/analyzed/";
std::string unpack_dir = m_workspace+"/temp_binary/";
std::string binary_dir = m_workspace+"/raw_binary/";
@ -366,13 +392,16 @@ void GWMEventBuilder::Convert2FastAnalyzedRoot() {
CompassRun converter(unpack_dir);
converter.SetShiftMap(m_shiftfile);
converter.SetScalerInput(m_scalerfile);
if(m_pb) converter.AttachProgressBar(m_pb);
if(m_pb)
converter.AttachProgressBar(m_pb);
std::cout<<"Beginning conversion..."<<std::endl;
for(int i=m_rmin; i<=m_rmax; i++) {
for(int i=m_rmin; i<=m_rmax; i++)
{
binfile = grabber.GrabFile(i);
if(binfile == "") continue;
if(binfile == "")
continue;
converter.SetRunNumber(i);
std::cout<<"Converting file: "<<binfile<<std::endl;
@ -389,9 +418,11 @@ void GWMEventBuilder::Convert2FastAnalyzedRoot() {
std::cout<<"-------------------------------------------"<<std::endl;
}
bool GWMEventBuilder::SetKinematicParameters(int zt, int at, int zp, int ap, int ze, int ae, double b, double theta, double bke) {
bool EVBApp::SetKinematicParameters(int zt, int at, int zp, int ap, int ze, int ae, double b, double theta, double bke)
{
if((at + ap - ae) < 0 || (zt + zp - ze) < 0) {
if((at + ap - ae) < 0 || (zt + zp - ze) < 0)
{
std::cout<<"Invalid Parameters at SetKinematicParameters"<<std::endl;
return false;
}

88
src/evb/FP_kinematics.cpp
View File

@ -39,63 +39,65 @@
//requires (Z,A) for T, P, and E, as well as energy of P,
// spectrograph angle of interest, and field value
double Delta_Z(int ZT, int AT, int ZP, int AP, int ZE, int AE,
double EP, double angle, double B) {
double EP, double angle, double B)
{
/* CONSTANTS */
const double UTOMEV = 931.4940954; //MeV per u;
const double MEVTOJ = 1.60218E-13; //J per MeV
const double RESTMASS_ELECTRON = 0.000548579909; //amu
const double UNIT_CHARGE = 1.602E-19; //Coulombs
const double C = 2.9979E8; //m/s
/* CONSTANTS */
const double UTOMEV = 931.4940954; //MeV per u;
const double MEVTOJ = 1.60218E-13; //J per MeV
const double RESTMASS_ELECTRON = 0.000548579909; //amu
const double UNIT_CHARGE = 1.602E-19; //Coulombs
const double C = 2.9979E8; //m/s
/* SESPS-SPECIFIC */
const double DISP = 1.96; //dispersion (x/rho)
const double MAG = 0.39; //magnification in x
const double DEGTORAD = M_PI/180.;
/* SESPS-SPECIFIC */
const double DISP = 1.96; //dispersion (x/rho)
const double MAG = 0.39; //magnification in x
const double DEGTORAD = M_PI/180.;
int ZR = ZT + ZP - ZE, AR = AT + AP - AE;
double EE=0; //ejectile energy
int ZR = ZT + ZP - ZE, AR = AT + AP - AE;
double EE=0; //ejectile energy
double MT=0, MP=0, ME=0, MR=0; //masses (MeV)
double MT=0, MP=0, ME=0, MR=0; //masses (MeV)
B /= 10000; //convert to tesla
angle *= DEGTORAD;
B /= 10000; //convert to tesla
angle *= DEGTORAD;
MT = MASS.FindMass(ZT, AT) - ZT*RESTMASS_ELECTRON*UTOMEV;
MP = MASS.FindMass(ZP, AP) - ZP*RESTMASS_ELECTRON*UTOMEV;
ME = MASS.FindMass(ZE, AE) - ZE*RESTMASS_ELECTRON*UTOMEV;
MR = MASS.FindMass(ZR, AR) - ZR*RESTMASS_ELECTRON*UTOMEV;
if (MT*MP*ME*MR == 0) {
std::cerr << "***WARNING: error loading one or more masses; returning 0\n";
return 0;
}
MT = MASS.FindMass(ZT, AT) - ZT*RESTMASS_ELECTRON*UTOMEV;
MP = MASS.FindMass(ZP, AP) - ZP*RESTMASS_ELECTRON*UTOMEV;
ME = MASS.FindMass(ZE, AE) - ZE*RESTMASS_ELECTRON*UTOMEV;
MR = MASS.FindMass(ZR, AR) - ZR*RESTMASS_ELECTRON*UTOMEV;
if (MT*MP*ME*MR == 0)
{
std::cerr << "***WARNING: error loading one or more masses; returning 0\n";
return 0;
}
double Q = MT + MP - ME - MR; //Q-value
//kinematics a la Iliadis p.590
double term1 = sqrt(MP*ME*EP)/(ME + MR)*cos(angle);
double term2 = (EP*(MR - MP) + MR*Q)/(ME + MR);
double Q = MT + MP - ME - MR; //Q-value
//kinematics a la Iliadis p.590
double term1 = sqrt(MP*ME*EP)/(ME + MR)*cos(angle);
double term2 = (EP*(MR - MP) + MR*Q)/(ME + MR);
EE = term1 + sqrt(term1*term1 + term2);
EE *= EE;
EE = term1 + sqrt(term1*term1 + term2);
EE *= EE;
//momentum
double PE = sqrt(EE*(EE+2*ME));
//momentum
double PE = sqrt(EE*(EE+2*ME));
//calculate rho from B a la B*rho = (proj. momentum)/(proj. charge)
double rho = (PE*MEVTOJ)/(ZE*UNIT_CHARGE*C*B)*100; //in cm
//calculate rho from B a la B*rho = (proj. momentum)/(proj. charge)
double rho = (PE*MEVTOJ)/(ZE*UNIT_CHARGE*C*B)*100; //in cm
double K;
double K;
K = sqrt(MP*ME*EP/EE);
K *= sin(angle);
K = sqrt(MP*ME*EP/EE);
K *= sin(angle);
double denom = ME + MR - sqrt(MP*ME*EP/EE)*cos(angle);
double denom = ME + MR - sqrt(MP*ME*EP/EE)*cos(angle);
K /= denom;
std::cout<<"Delta Z= "<<-1*rho*DISP*MAG*K<<std::endl;
return -1*rho*DISP*MAG*K; //delta-Z in cm
K /= denom;
std::cout<<"Delta Z= "<<-1*rho*DISP*MAG*K<<std::endl;
return -1*rho*DISP*MAG*K; //delta-Z in cm
}

183
src/evb/FastSort.cpp
View File

@ -1,24 +1,26 @@
#include "EventBuilder.h"
#include "FastSort.h"
FastSort::FastSort(float si_windowSize, float ion_windowSize) {
si_coincWindow = si_windowSize/1.0e3; //given in pico s, want in nano s
ion_coincWindow = ion_windowSize/1.0e3;
event_address = NULL;
//windows given in picoseconds, converted to nanoseconds
FastSort::FastSort(float si_windowSize, float ion_windowSize)
si_coincWindow(si_windowSize/1.0e3), ion_coincWindow(ion_windowSize/1.0e3), event_address(nullptr)
{
}
FastSort::~FastSort() {
delete event_address;
FastSort::~FastSort()
{
delete event_address;
}
void FastSort::ResetSABRE() {
for(int i=0; i<5; i++) {
fastEvent.sabreArray[i] = sblank;
}
void FastSort::ResetSABRE()
{
for(int i=0; i<5; i++)
fastEvent.sabreArray[i] = sblank;
}
void FastSort::ResetFocalPlane() {
fastEvent.focalPlane = fpblank;
void FastSort::ResetFocalPlane()
{
fastEvent.focalPlane = fpblank;
}
/*Assign a set of ion chamber data to the scintillator*/
@ -28,89 +30,96 @@ void FastSort::ProcessFocalPlane(unsigned int scint_index, unsigned int ionch_in
*In this case, I chose one of the anodes. But in principle you could also choose any other part of the ion
*chamber
*/
if(slowEvent.focalPlane.anodeB.size() > ionch_index) { //Back anode required to move on`
if(slowEvent.focalPlane.anodeB.size() > ionch_index)
{ //Back anode required to move on`
float anodeRelTime = fabs(slowEvent.focalPlane.anodeB[ionch_index].Time - slowEvent.focalPlane.scintL[scint_index].Time);
if(anodeRelTime > ion_coincWindow) return; //Window check
float anodeRelTime = fabs(slowEvent.focalPlane.anodeB[ionch_index].Time - slowEvent.focalPlane.scintL[scint_index].Time);
if(anodeRelTime > ion_coincWindow)
return; //Window check
fastEvent.focalPlane.anodeB.push_back(slowEvent.focalPlane.anodeB[ionch_index]);
fastEvent.focalPlane.scintL.push_back(slowEvent.focalPlane.scintL[scint_index]);
if(slowEvent.focalPlane.delayFL.size() > ionch_index) {
fastEvent.focalPlane.delayFL.push_back(slowEvent.focalPlane.delayFL[ionch_index]);
}
if(slowEvent.focalPlane.delayFR.size() > ionch_index) {
fastEvent.focalPlane.delayFR.push_back(slowEvent.focalPlane.delayFR[ionch_index]);
}
if(slowEvent.focalPlane.delayBR.size() > ionch_index) {
fastEvent.focalPlane.delayBR.push_back(slowEvent.focalPlane.delayBR[ionch_index]);
}
if(slowEvent.focalPlane.delayBL.size() > ionch_index) {
fastEvent.focalPlane.delayBL.push_back(slowEvent.focalPlane.delayBL[ionch_index]);
}
if(slowEvent.focalPlane.scintR.size() > ionch_index) {
fastEvent.focalPlane.scintR.push_back(slowEvent.focalPlane.scintR[ionch_index]);
}
if(slowEvent.focalPlane.anodeF.size() > ionch_index) {
fastEvent.focalPlane.anodeF.push_back(slowEvent.focalPlane.anodeF[ionch_index]);
}
if(slowEvent.focalPlane.cathode.size() > ionch_index) {
fastEvent.focalPlane.cathode.push_back(slowEvent.focalPlane.cathode[ionch_index]);
}
}
fastEvent.focalPlane.anodeB.push_back(slowEvent.focalPlane.anodeB[ionch_index]);
fastEvent.focalPlane.scintL.push_back(slowEvent.focalPlane.scintL[scint_index]);
if(slowEvent.focalPlane.delayFL.size() > ionch_index)
fastEvent.focalPlane.delayFL.push_back(slowEvent.focalPlane.delayFL[ionch_index]);
if(slowEvent.focalPlane.delayFR.size() > ionch_index)
fastEvent.focalPlane.delayFR.push_back(slowEvent.focalPlane.delayFR[ionch_index]);
if(slowEvent.focalPlane.delayBR.size() > ionch_index)
fastEvent.focalPlane.delayBR.push_back(slowEvent.focalPlane.delayBR[ionch_index]);
if(slowEvent.focalPlane.delayBL.size() > ionch_index)
fastEvent.focalPlane.delayBL.push_back(slowEvent.focalPlane.delayBL[ionch_index]);
if(slowEvent.focalPlane.scintR.size() > ionch_index)
fastEvent.focalPlane.scintR.push_back(slowEvent.focalPlane.scintR[ionch_index]);
if(slowEvent.focalPlane.anodeF.size() > ionch_index)
fastEvent.focalPlane.anodeF.push_back(slowEvent.focalPlane.anodeF[ionch_index]);
if(slowEvent.focalPlane.cathode.size() > ionch_index)
fastEvent.focalPlane.cathode.push_back(slowEvent.focalPlane.cathode[ionch_index]);
}
}
/*Assign a set of SABRE data that falls within the coincidence window*/
void FastSort::ProcessSABRE(unsigned int scint_index) {
for(int i=0; i<5; i++) { //loop over SABRE silicons
std::vector<DetectorHit> rings;
std::vector<DetectorHit> wedges;
void FastSort::ProcessSABRE(unsigned int scint_index)
{
for(int i=0; i<5; i++)
{ //loop over SABRE silicons
std::vector<DetectorHit> rings;
std::vector<DetectorHit> wedges;
if(slowEvent.sabreArray[i].rings.size() == 0 || slowEvent.sabreArray[i].wedges.size() == 0) continue; //save some time on empties
if(slowEvent.sabreArray[i].rings.size() == 0 || slowEvent.sabreArray[i].wedges.size() == 0)
continue; //save some time on empties
/*Dump sabre data that doesnt fall within the fast coincidence window with the scint*/
for(unsigned int j=0; j<slowEvent.sabreArray[i].rings.size(); j++) {
float sabreRelTime = fabs(slowEvent.sabreArray[i].rings[j].Time - slowEvent.focalPlane.scintL[scint_index].Time);
if(sabreRelTime < si_coincWindow) {
rings.push_back(slowEvent.sabreArray[i].rings[j]);
}
}
for(unsigned int j=0; j<slowEvent.sabreArray[i].wedges.size(); j++) {
float sabreRelTime = fabs(slowEvent.sabreArray[i].wedges[j].Time - slowEvent.focalPlane.scintL[scint_index].Time);
if(sabreRelTime < si_coincWindow) {
wedges.push_back(slowEvent.sabreArray[i].wedges[j]);
}
}
fastEvent.sabreArray[i].rings = rings;
fastEvent.sabreArray[i].wedges = wedges;
}
/*Dump sabre data that doesnt fall within the fast coincidence window with the scint*/
for(unsigned int j=0; j<slowEvent.sabreArray[i].rings.size(); j++)
{
float sabreRelTime = fabs(slowEvent.sabreArray[i].rings[j].Time - slowEvent.focalPlane.scintL[scint_index].Time);
if(sabreRelTime < si_coincWindow)
rings.push_back(slowEvent.sabreArray[i].rings[j]);
}
for(unsigned int j=0; j<slowEvent.sabreArray[i].wedges.size(); j++)
{
float sabreRelTime = fabs(slowEvent.sabreArray[i].wedges[j].Time - slowEvent.focalPlane.scintL[scint_index].Time);
if(sabreRelTime < si_coincWindow)
wedges.push_back(slowEvent.sabreArray[i].wedges[j]);
}
fastEvent.sabreArray[i].rings = rings;
fastEvent.sabreArray[i].wedges = wedges;
}
}
std::vector<CoincEvent> FastSort::GetFastEvents(CoincEvent& event) {
slowEvent = event;
std::vector<CoincEvent> fast_events;
std::vector<CoincEvent> FastSort::GetFastEvents(CoincEvent& event)
{
slowEvent = event;
std::vector<CoincEvent> fast_events;
unsigned int sizeArray[7];
sizeArray[0] = slowEvent.focalPlane.delayFL.size();
sizeArray[1] = slowEvent.focalPlane.delayFR.size();
sizeArray[2] = slowEvent.focalPlane.delayBL.size();
sizeArray[3] = slowEvent.focalPlane.delayBR.size();
sizeArray[4] = slowEvent.focalPlane.anodeF.size();
sizeArray[5] = slowEvent.focalPlane.anodeB.size();
sizeArray[6] = slowEvent.focalPlane.cathode.size();
unsigned int maxSize = *max_element(sizeArray, sizeArray+7);
//loop over scints
for(unsigned int i=0; i<slowEvent.focalPlane.scintL.size(); i++) {
ResetSABRE();
ProcessSABRE(i);
//loop over ion chamber
//NOTE: as written, this dumps data that does not have an ion chamber hit!
//If you want scint/SABRE singles, move the fill outside of this loop
for(unsigned int j=0; j<maxSize; j++) {
ResetFocalPlane();
ProcessFocalPlane(i, j);
fast_events.push_back(fastEvent);
}
}
return fast_events;
unsigned int sizeArray[7];
sizeArray[0] = slowEvent.focalPlane.delayFL.size();
sizeArray[1] = slowEvent.focalPlane.delayFR.size();
sizeArray[2] = slowEvent.focalPlane.delayBL.size();
sizeArray[3] = slowEvent.focalPlane.delayBR.size();
sizeArray[4] = slowEvent.focalPlane.anodeF.size();
sizeArray[5] = slowEvent.focalPlane.anodeB.size();
sizeArray[6] = slowEvent.focalPlane.cathode.size();
unsigned int maxSize = *max_element(sizeArray, sizeArray+7);
//loop over scints
for(unsigned int i=0; i<slowEvent.focalPlane.scintL.size(); i++)
{
ResetSABRE();
ProcessSABRE(i);
//loop over ion chamber
//NOTE: as written, this dumps data that does not have an ion chamber hit!
//If you want scint/SABRE singles, move the fill outside of this loop
for(unsigned int j=0; j<maxSize; j++)
{
ResetFocalPlane();
ProcessFocalPlane(i, j);
fast_events.push_back(fastEvent);
}
}
return fast_events;
}

64
src/evb/FlagHandler.cpp
View File

@ -2,65 +2,85 @@
#include "FlagHandler.h"
FlagHandler::FlagHandler() :
log("./event_log.txt")
log("./event_log.txt")
{
}
FlagHandler::FlagHandler(const std::string& filename) :
log(filename)
log(filename)
{
}
FlagHandler::~FlagHandler() {
FlagHandler::~FlagHandler()
{
WriteLog();
log.close();
}
void FlagHandler::CheckFlag(int board, int channel, int flag) {
void FlagHandler::CheckFlag(int board, int channel, int flag)
{
int gchan = channel + board*16;
FlagCount& counter = event_count_map[gchan]; //yikes
counter.total_counts++;
if(flag & DEAD_TIME) counter.dead_time++;
if(flag & DeadTime)
counter.dead_time++;
if(flag & TIME_ROLLOVER) counter.time_roll++;
if(flag & TimeRollover)
counter.time_roll++;
if(flag & TIME_RESET) counter.time_reset++;
if(flag & TimeReset)
counter.time_reset++;
if(flag & FAKE_EVENT) counter.fake_event++;
if(flag & FakeEvent)
counter.fake_event++;
if(flag & MEM_FULL) counter.mem_full++;
if(flag & MemFull)
counter.mem_full++;
if(flag & TRIG_LOST) counter.trig_lost++;
if(flag & TrigLost)
counter.trig_lost++;
if(flag & N_TRIG_LOST) counter.n_trig_lost++;
if(flag & NTrigLost)
counter.n_trig_lost++;
if(flag & SATURATING_IN_GATE) counter.sat_in_gate++;
if(flag & SaturatingInGate)
counter.sat_in_gate++;
if(flag & TRIG_1024_COUNTED) counter.trig_1024++;;
if(flag & Trig1024Counted)
counter.trig_1024++;;
if(flag & SATURATING_INPUT) counter.sat_input++;
if(flag & SaturatingInput)
counter.sat_input++;
if(flag & N_TRIG_COUNTED) counter.n_trig_count++;
if(flag & NTrigCounted)
counter.n_trig_count++;
if(flag & EVENT_NOT_MATCHED) counter.event_not_matched++;
if(flag & EventNotMatched)
counter.event_not_matched++;
if(flag & PILE_UP) counter.pile_up++;
if(flag & PileUp)
counter.pile_up++;
if(flag & PLL_LOCK_LOSS) counter.pll_lock_loss++;
if(flag & PLLLockLoss)
counter.pll_lock_loss++;
if(flag & OVER_TEMP) counter.over_temp++;
if(flag & OverTemp)
counter.over_temp++;
if(flag & ADC_SHUTDOWN) counter.adc_shutdown++;
if(flag & ADCShutdown)
counter.adc_shutdown++;
}
void FlagHandler::WriteLog() {
void FlagHandler::WriteLog()
{
log<<"Event Flag Log"<<std::endl;
log<<"-----------------------------"<<std::endl;
for(auto& counter : event_count_map) {
for(auto& counter : event_count_map)
{
log<<"-----------------------------"<<std::endl;
log<<"GLOBAL CHANNEL No.: "<<counter.first<<std::endl;
log<<"Total number of events: "<<counter.second.total_counts<<std::endl;

74
src/evb/MassLookup.cpp
View File

@ -16,46 +16,52 @@ Written by G.W. McCann Aug. 2020
Read in AMDC mass file, preformated to remove excess info. Here assumes that by default
the file is in a local directory etc/
*/
MassLookup::MassLookup() {
std::ifstream massfile("./etc/mass.txt");
if(massfile.is_open()) {
std::string junk, A, element;
int Z;
double atomicMassBig, atomicMassSmall, isotopicMass;
getline(massfile,junk);
getline(massfile,junk);
while(massfile>>junk) {
massfile>>Z>>A>>element>>atomicMassBig>>atomicMassSmall;
isotopicMass = (atomicMassBig + atomicMassSmall*1e-6 - Z*electron_mass)*u_to_mev;
std::string key = "("+std::to_string(Z)+","+A+")";
massTable[key] = isotopicMass;
elementTable[Z] = element;
}
} else {
std::cerr<<"Unable to open mass.txt at MassLookup! Prepare for errors."<<std::endl;
}
MassLookup::MassLookup()
{
std::ifstream massfile("./etc/mass.txt");
if(massfile.is_open())
{
int Z;
double atomicMassBig, atomicMassSmall, isotopicMass;
getline(massfile,junk);
getline(massfile,junk);
while(massfile>>junk)
{
massfile>>Z>>A>>element>>atomicMassBig>>atomicMassSmall;
isotopicMass = (atomicMassBig + atomicMassSmall*1e-6 - Z*electron_mass)*u_to_mev;
std::string key = "("+std::to_string(Z)+","+A+")";
massTable[key] = isotopicMass;
elementTable[Z] = element;
}
}
else
std::cerr<<"Unable to open mass.txt at MassLookup! Prepare for errors."<<std::endl;
}
MassLookup::~MassLookup() {}
//Returns nuclear mass in MeV
double MassLookup::FindMass(int Z, int A) {
std::string key = "("+std::to_string(Z)+","+std::to_string(A)+")";
auto data = massTable.find(key);
if(data == massTable.end()) {
std::cerr<<"Invaild nucleus at MassLookup! Returning mass of 0"<<std::endl;
return 0;
}
return data->second;
double MassLookup::FindMass(int Z, int A)
{
std::string key = "("+std::to_string(Z)+","+std::to_string(A)+")";
auto data = massTable.find(key);
if(data == massTable.end())
{
std::cerr<<"Invaild nucleus at MassLookup! Returning mass of 0"<<std::endl;
return 0;
}
return data->second;
}
//returns element symbol
std::string MassLookup::FindSymbol(int Z, int A) {
auto data = elementTable.find(Z);
if(data == elementTable.end()) {
std::cerr<<"Invaild nucleus at MassLookup! Returning empty symbol"<<std::endl;
return "";
}
std::string fullsymbol = std::to_string(A) + data->second;
return fullsymbol;
std::string MassLookup::FindSymbol(int Z, int A)
{
auto data = elementTable.find(Z);
if(data == elementTable.end())
{
std::cerr<<"Invaild nucleus at MassLookup! Returning empty symbol"<<std::endl;
return "";
}
std::string fullsymbol = std::to_string(A) + data->second;
return fullsymbol;
}

21
src/evb/OrderChecker.cpp
View File

@ -8,25 +8,28 @@
#include "EventBuilder.h"
#include "OrderChecker.h"
OrderChecker::OrderChecker() {
OrderChecker::OrderChecker()
{
}
OrderChecker::~OrderChecker() {
OrderChecker::~OrderChecker()
{
}
bool OrderChecker::IsOrdered(const std::string& filename) {
bool OrderChecker::IsOrdered(const std::string& filename)
{
TFile* file = TFile::Open(filename.c_str(), "READ");
TTree* tree = (TTree*) file->Get("Data");
ULong64_t ts;
uint64_t ts;
tree->SetBranchAddress("Timestamp", &ts);
ULong64_t prevStamp = 0;
uint64_t prevStamp = 0;
for(Long64_t i=0; i<tree->GetEntries(); i++) {
for(Long64_t i=0; i<tree->GetEntries(); i++)
{
tree->GetEntry();
if(prevStamp >= ts) {
if(prevStamp >= ts)
{
std::cerr<<"Bad order at entry "<<i<<" out of "<<tree->GetEntries()<<std::endl;
return false;
}

371
src/evb/RunCollector.cpp
View File

@ -1,196 +1,233 @@
#include "EventBuilder.h"
#include "RunCollector.h"
using namespace std;
#include <TSystemDirectory.h>
#include <TSystemFile.h>
#include <TCollection.h>
#include <TList.h>
#include <cstdlib>
#include <cstdio>
RunCollector::RunCollector():
initFlag(false), dir(""), run(""), end(""), MaxRun(0), MinRun(0)
m_initFlag(false), m_directory(""), m_prefix(""), m_suffix(""), m_maxRun(0), m_minRun(0)
{
}
RunCollector::RunCollector(const string& dirname, const string& prefix, const string& suffix) {
dir = dirname.c_str();
run = prefix.c_str();
end = suffix.c_str();
MinRun = 0; MaxRun = LITERALMAX;
initFlag = true;
RunCollector::RunCollector(const std::string& dirname, const std::string& prefix, const std::string& suffix)
m_initFlag(true), m_directory(dirname), m_prefix(prefix), m_suffix(suffix), m_minRun(0), m_maxRun(0)
{
}
RunCollector::RunCollector(const string& dirname, const string& prefix, const string& suffix, int min, int max) {
dir = dirname.c_str();
run = prefix.c_str();
end = suffix.c_str();
MinRun = min; MaxRun = max;
initFlag = true;
RunCollector::RunCollector(const std::string& dirname, const std::string& prefix, const std::string& suffix, int min, int max)
m_initFlag(true), m_directory(dirname), m_prefix(prefix), m_suffix(suffix), m_minRun(min), m_maxRun(max)
{
}
RunCollector::~RunCollector() {}
void RunCollector::SetSearchParams(const string& dirname, const string& prefix, const string& suffix, int min, int max) {
dir = dirname.c_str();
run = prefix.c_str();
end = suffix.c_str();
MinRun = min; MaxRun = max;
initFlag = true;
void RunCollector::SetSearchParams(const std::string& dirname, const std::string& prefix, const std::string& suffix, int min, int max)
{
m_directory = dirname.c_str();
m_prefix = prefix.c_str();
m_suffix = suffix.c_str();
m_minRun = min; m_maxRun = max;
m_initFlag = true;
}
int RunCollector::GrabAllFiles() {
if(!initFlag) {return 0;}
TSystemDirectory sysdir(dir.Data(), dir.Data());
TList *flist = sysdir.GetListOfFiles();
filelist.clear();
int counter = 0;
if(flist) { //Make sure list is real. If not, means no directory
TSystemFile *file;
TString fname, temp;
TIter next_element(flist); //List iterator
while((file = (TSystemFile*)next_element())) {
temp = file->GetName();
if(!file->IsDirectory() && temp.BeginsWith(run.Data()) && temp.EndsWith(end.Data())) {
counter++;
fname = dir+temp;
filelist.push_back(fname);
}
}
if(counter>0) {
delete flist;
return 1;
} else {
cerr<<"Unable to find files with matching run name in directory; check input.txt"<<endl;
delete flist;
return 0;
}
} else {
cerr<<"Unable to find any files in directory; check name given to the input.txt"<<endl;
delete flist;
return 0;
}
bool RunCollector::GrabAllFiles()
{
if(!m_initFlag)
return false;
TSystemDirectory sysdir(m_directory.c_str(), m_directory.c_str());
TList *flist = sysdir.GetListOfFiles();
m_filelist.clear();
if(!flist) //Make sure list is real. If not, means no directory
{
cerr<<"Unable to find any files in directory; check name given to the input.txt"<<endl;
return false;
}
TSystemFile *file;
std::string fname, temp;
TIter next_element(flist); //List iterator
while((file = (TSystemFile*)next_element()))
{
temp = file->GetName();
if(temp.size() < prefix.size() || temp.size() < suffix.size())
continue;
else if(!file->IsDirectory() && !temp.compare(0,prefix.size(),prefix) &&
!temp.compare(temp.size()-(1+suffix.size()), suffix.size(), suffix))
{
fname = dir+temp;
m_filelist.push_back(fname);
}
}
delete flist;
if(m_filelist.size()>0)
return true;
else
{
cerr<<"Unable to find files with matching run name in directory; check input.txt"<<endl;
return false;
}
}
std::string RunCollector::GrabFile(int runNum) {
if(!initFlag) return "";
TSystemDirectory sysdir(dir.Data(), dir.Data());
TList* flist = sysdir.GetListOfFiles();
if(!m_initFlag)
return "";
TSystemDirectory sysdir(dir.Data(), dir.Data());
TList* flist = sysdir.GetListOfFiles();
if(!flist) return "";
if(!flist)
return "";
TSystemFile *file;
TString fname = "", temp;
string runno = "_"+to_string(runNum)+end.Data();
TIter next_element(flist);
while((file = (TSystemFile*)next_element())) {
temp = file->GetName();
if(!file->IsDirectory() && temp.BeginsWith(run.Data()) && temp.EndsWith(runno.c_str())) {
fname = dir+temp;
break;
}
}
TSystemFile *file;
std::string fname = "", temp;
std::string runno = "_"+to_string(runNum)+m_suffix;
TIter next_element(flist);
while((file = (TSystemFile*)next_element()))
{
temp = file->GetName();
if(temp.size() < m_prefix.size() || temp.size() < runno.size())
continue;
else if(!file->IsDirectory() && !temp.compare(0,m_prefix.size(),m_prefix) &&
!temp.compare(temp.size()-(1+runno.size()),runno.size(), runno))
{
fname = dir+temp;
break;
}
}
return fname.Data();
delete flist;
return fname;
}
/*Grabs all files within a specified run range*/
int RunCollector::GrabFilesInRange() {
if(!initFlag) {return 0;}
TSystemDirectory sysdir(dir.Data(), dir.Data());
TList *flist = sysdir.GetListOfFiles();
filelist.clear();
int counter = 0;
if(flist) {
TSystemFile *file;
TString fname, temp;
string runno;
for(int i=MinRun; i<=MaxRun; i++) {//loop over range
TIter next_element(flist);//list iterator
runno = "_"+to_string(i) + end.Data(); //suffix is now _#.endData
while((file = (TSystemFile*)next_element())) {//look through directory until file found
temp = file->GetName();
if(!file->IsDirectory()&&temp.BeginsWith(run.Data())&&temp.EndsWith(runno.c_str())){
counter++;
fname = dir+temp;
filelist.push_back(fname);
break; //if we find the file, break out of iterator loop
}
}
}
if(counter>0) {
delete flist;
return 1;
} else {
cerr<<"Unable to find files with matching run name in directory; check input.txt"<<endl;
delete flist;
return 0;
}
} else {
cerr<<"Unable to find any files in directory; check name given to input.txt"<<endl;
delete flist;
return 0;
}
int RunCollector::GrabFilesInRange()
{
if(!m_initFlag)
return false;
TSystemDirectory sysdir(dir.Data(), dir.Data());
TList *flist = sysdir.GetListOfFiles();
m_filelist.clear();
int counter = 0;
if(!flist)
{
cerr<<"Unable to find any files in directory; check name given to input.txt"<<endl;
return false;
}
TSystemFile *file;
std::string fname, temp;
std::string runno;
for(int i=m_minRun; i<=m_maxRun; i++) //loop over range
{
TIter next_element(flist);//list iterator
runno = "_"+to_string(i) + m_suffix; //suffix is now _#.suffix
while((file = (TSystemFile*)next_element())) //look through directory until file found
{
temp = file->GetName();
if(temp.size() < m_prefix.size() || temp.size() < runno.size())
continue;
else if(!file->IsDirectory() && !temp.compare(0,m_prefix.size(),m_prefix) &&
!temp.compare(temp.size()-(1+runno.size()),runno.size(), runno))
{
counter++;
fname = dir+temp;
m_filelist.push_back(fname);
break; //if we find the file, break out of iterator loop
}
}
}
delete flist;
if(counter>0)
return true;
else
{
cerr<<"Unable to find files with matching run name in directory; check input.txt"<<endl;
return false;
}
}
int RunCollector::Merge_hadd(const string& outname) {
if(!initFlag) {return 0;}
if(MaxRun == LITERALMAX) {
if(GrabAllFiles()) {
TString clump = "hadd "+outname;
for(unsigned int i=0; i<filelist.size(); i++) {
clump += " "+filelist[i];
}
cout<<"Merging runs into single file..."<<endl;
system(clump.Data());
cout<<"Finished merging"<<endl;
return 1;
} else {
return 0;
}
} else {
if(GrabFilesInRange()) {
TString clump = "hadd "+outname;
for(unsigned int i=0; i<filelist.size(); i++) {
clump += " "+filelist[i];
}
cout<<"Merging runs "<<MinRun<<" to "<<MaxRun<<" into a single file..."<<endl;
system(clump.Data());
cout<<"Finished merging"<<endl;
return 1;
} else {
return 0;
}
}
bool RunCollector::Merge_hadd(const std::string& outname)
{
if(!m_initFlag)
false;
if(m_maxRun == 0)
{
if(GrabAllFiles())
{
std::string clump = "hadd "+outname;
for(unsigned int i=0; i<filelist.size(); i++)
clump += " "+filelist[i];
cout<<"Merging runs into single file..."<<endl;
std::system(clump);
cout<<"Finished merging"<<endl;
return true;
}
else
return false;
}
else
{
if(GrabFilesInRange())
{
std::string clump = "hadd "+outname;
for(unsigned int i=0; i<filelist.size(); i++)
clump += " "+filelist[i];
cout<<"Merging runs "<<m_minRun<<" to "<<m_maxRun<<" into a single file..."<<endl;
system(clump);
cout<<"Finished merging"<<endl;
return true;
}
else
return false;
}
}
int RunCollector::Merge_TChain(const string& outname) {
if(!initFlag) {return 0;}
TFile *output = new TFile(outname.c_str(), "RECREATE");
TChain *chain = new TChain("SPSTree", "SPSTree");
if(MaxRun == LITERALMAX) {
if(GrabAllFiles()) {
for(unsigned int i=0; i<filelist.size(); i++) {
chain->Add(filelist[i].Data());
}
cout<<"Merging runs into single file..."<<endl;
chain->Merge(output,0,"fast");
cout<<"Finished merging"<<endl;
return 1;
} else {
return 0;
}
} else {
if(GrabFilesInRange()) {
for(unsigned int i=0; i<filelist.size(); i++) {
chain->Add(filelist[i]);
}
cout<<"Merging runs "<<MinRun<<" to "<<MaxRun<<" into a single file..."<<endl;
chain->Merge(output,0,"fast");
cout<<"Finished merging"<<endl;
return 1;
} else {
return 0;
}
}
if(output->IsOpen()) output->Close();
return 0;
bool RunCollector::Merge_TChain(const std::string& outname)
{
if(!m_initFlag)
return false;
TFile *output = new TFile(outname.c_str(), "RECREATE");
TChain *chain = new TChain("SPSTree", "SPSTree");
if(m_maxRun == 0)
{
if(GrabAllFiles())
{
for(unsigned int i=0; i<m_filelist.size(); i++)
chain->Add(m_filelist[i].Data());
cout<<"Merging runs into single file..."<<endl;
chain->Merge(output,0,"fast");
cout<<"Finished merging"<<endl;
return true;
}
else
return false;
}
else
{
if(GrabFilesInRange())
{
for(unsigned int i=0; i<m_filelist.size(); i++)
chain->Add(m_filelist[i]);
cout<<"Merging runs "<<m_minRun<<" to "<<m_maxRun<<" into a single file..."<<endl;
chain->Merge(output,0,"fast");
cout<<"Finished merging"<<endl;
return true;
} else
return false;
}
if(output->IsOpen())
output->Close();
return false;
}

306
src/evb/SFPAnalyzer.cpp
View File

@ -12,25 +12,27 @@
#include "EventBuilder.h"
#include "SFPAnalyzer.h"
using namespace std;
/*Constructor takes in kinematic parameters for generating focal plane weights*/
SFPAnalyzer::SFPAnalyzer(int zt, int at, int zp, int ap, int ze, int ae, double ep,
double angle, double b) {
zfp = Delta_Z(zt, at, zp, ap, ze, ae, ep, angle, b);
event_address = new CoincEvent();
rootObj = new THashTable();
GetWeights();
double angle, double b)
{
zfp = Delta_Z(zt, at, zp, ap, ze, ae, ep, angle, b);
event_address = new CoincEvent();
rootObj = new THashTable();
GetWeights();
}
SFPAnalyzer::~SFPAnalyzer() {
rootObj->Clear();
delete rootObj;
delete event_address;
SFPAnalyzer::~SFPAnalyzer()
{
rootObj->Clear();
delete rootObj;
delete event_address;
}
void SFPAnalyzer::Reset() {
pevent = blank; //set output back to blank
void SFPAnalyzer::Reset()
{
pevent = blank; //set output back to blank
}
/*Use functions from FP_kinematics to calculate weights for xavg
@ -38,146 +40,166 @@ void SFPAnalyzer::Reset() {
*from the two focal plane points and finding the intersection with
*the kinematic focal plane
*/
void SFPAnalyzer::GetWeights() {
w1 = (Wire_Dist()/2.0-zfp)/Wire_Dist();
w2 = 1.0-w1;
cout<<"w1: "<<w1<<" w2: "<<w2<<endl;
void SFPAnalyzer::GetWeights()
{
w1 = (Wire_Dist()/2.0-zfp)/Wire_Dist();
w2 = 1.0-w1;
cout<<"w1: "<<w1<<" w2: "<<w2<<endl;
}
/*2D histogram fill wrapper for use with THashTable (faster)*/
void SFPAnalyzer::MyFill(const string& name, int binsx, double minx, double maxx, double valuex,
int binsy, double miny, double maxy, double valuey) {
TH2F *histo = (TH2F*) rootObj->FindObject(name.c_str());
if(histo != NULL) {
histo->Fill(valuex, valuey);
} else {
TH2F *h = new TH2F(name.c_str(), name.c_str(), binsx, minx, maxx, binsy, miny, maxy);
h->Fill(valuex, valuey);
rootObj->Add(h);
}
void SFPAnalyzer::MyFill(const std::string& name, int binsx, double minx, double maxx, double valuex,
int binsy, double miny, double maxy, double valuey)
{
TH2F *histo = (TH2F*) rootObj->FindObject(name.c_str());
if(histo != nullptr)
histo->Fill(valuex, valuey);
else
{
TH2F *h = new TH2F(name.c_str(), name.c_str(), binsx, minx, maxx, binsy, miny, maxy);
h->Fill(valuex, valuey);
rootObj->Add(h);
}
}
/*1D histogram fill wrapper for use with THashTable (faster)*/
void SFPAnalyzer::MyFill(const string& name, int binsx, double minx, double maxx, double valuex) {
TH1F *histo = (TH1F*) rootObj->FindObject(name.c_str());
if(histo != NULL) {
histo->Fill(valuex);
} else {
TH1F *h = new TH1F(name.c_str(), name.c_str(), binsx, minx, maxx);
h->Fill(valuex);
rootObj->Add(h);
}
void SFPAnalyzer::MyFill(const std::string& name, int binsx, double minx, double maxx, double valuex)
{
TH1F *histo = (TH1F*) rootObj->FindObject(name.c_str());
if(histo != nullptr)
histo->Fill(valuex);
else
{
TH1F *h = new TH1F(name.c_str(), name.c_str(), binsx, minx, maxx);
h->Fill(valuex);
rootObj->Add(h);
}
}
void SFPAnalyzer::AnalyzeEvent(CoincEvent& event) {
Reset();
if(!event.focalPlane.anodeF.empty()) {
pevent.anodeFront = event.focalPlane.anodeF[0].Long;
pevent.anodeFrontTime = event.focalPlane.anodeF[0].Time;
}
if(!event.focalPlane.anodeB.empty()) {
pevent.anodeBack = event.focalPlane.anodeB[0].Long;
pevent.anodeBackTime = event.focalPlane.anodeB[0].Time;
}
if(!event.focalPlane.scintL.empty()) {
pevent.scintLeft = event.focalPlane.scintL[0].Long;
pevent.scintLeftShort = event.focalPlane.scintL[0].Short;
pevent.scintLeftTime = event.focalPlane.scintL[0].Time;
}
if(!event.focalPlane.scintR.empty()) {
pevent.scintRight = event.focalPlane.scintR[0].Long;
pevent.scintRightShort = event.focalPlane.scintR[0].Short;
pevent.scintRightTime = event.focalPlane.scintR[0].Time;
}
if(!event.focalPlane.cathode.empty()) {
pevent.cathode = event.focalPlane.cathode[0].Long;
pevent.cathodeTime = event.focalPlane.cathode[0].Time;
}
if(!event.focalPlane.monitor.empty()) {
pevent.monitorE = event.focalPlane.monitor[0].Long;
pevent.monitorShort = event.focalPlane.monitor[0].Short;
pevent.monitorTime = event.focalPlane.monitor[0].Time;
}
void SFPAnalyzer::AnalyzeEvent(CoincEvent& event)
{
Reset();
if(!event.focalPlane.anodeF.empty())
{
pevent.anodeFront = event.focalPlane.anodeF[0].Long;
pevent.anodeFrontTime = event.focalPlane.anodeF[0].Time;
}
if(!event.focalPlane.anodeB.empty())
{
pevent.anodeBack = event.focalPlane.anodeB[0].Long;
pevent.anodeBackTime = event.focalPlane.anodeB[0].Time;
}
if(!event.focalPlane.scintL.empty())
{
pevent.scintLeft = event.focalPlane.scintL[0].Long;
pevent.scintLeftShort = event.focalPlane.scintL[0].Short;
pevent.scintLeftTime = event.focalPlane.scintL[0].Time;
}
if(!event.focalPlane.scintR.empty())
{
pevent.scintRight = event.focalPlane.scintR[0].Long;
pevent.scintRightShort = event.focalPlane.scintR[0].Short;
pevent.scintRightTime = event.focalPlane.scintR[0].Time;
}
if(!event.focalPlane.cathode.empty())
{
pevent.cathode = event.focalPlane.cathode[0].Long;
pevent.cathodeTime = event.focalPlane.cathode[0].Time;
}
if(!event.focalPlane.monitor.empty())
{
pevent.monitorE = event.focalPlane.monitor[0].Long;
pevent.monitorShort = event.focalPlane.monitor[0].Short;
pevent.monitorTime = event.focalPlane.monitor[0].Time;
}
/*Delay lines and all that*/
if(!event.focalPlane.delayFR.empty()) {
pevent.delayFrontRightE = event.focalPlane.delayFR[0].Long;
pevent.delayFrontRightTime = event.focalPlane.delayFR[0].Time;
pevent.delayFrontRightShort = event.focalPlane.delayFR[0].Short;
}
if(!event.focalPlane.delayFL.empty()) {
pevent.delayFrontLeftE = event.focalPlane.delayFL[0].Long;
pevent.delayFrontLeftTime = event.focalPlane.delayFL[0].Time;
pevent.delayFrontLeftShort = event.focalPlane.delayFL[0].Short;
}
if(!event.focalPlane.delayBR.empty()) {
pevent.delayBackRightE = event.focalPlane.delayBR[0].Long;
pevent.delayBackRightTime = event.focalPlane.delayBR[0].Time;
pevent.delayBackRightShort = event.focalPlane.delayBR[0].Short;
}
if(!event.focalPlane.delayBL.empty()) {
pevent.delayBackLeftE = event.focalPlane.delayBL[0].Long;
pevent.delayBackLeftTime = event.focalPlane.delayBL[0].Time;
pevent.delayBackLeftShort = event.focalPlane.delayBL[0].Short;
}
if(!event.focalPlane.delayFL.empty() && !event.focalPlane.delayFR.empty()) {
pevent.fp1_tdiff = (event.focalPlane.delayFL[0].Time-event.focalPlane.delayFR[0].Time)*0.5;
pevent.fp1_tsum = (event.focalPlane.delayFL[0].Time+event.focalPlane.delayFR[0].Time);
pevent.fp1_tcheck = (pevent.fp1_tsum)/2.0-pevent.anodeFrontTime;
pevent.delayFrontMaxTime = max(event.focalPlane.delayFL[0].Time, event.focalPlane.delayFR[0].Time);
pevent.x1 = pevent.fp1_tdiff*1.0/2.10; //position from time, based on total delay
MyFill("x1",1200,-300,300,pevent.x1);
MyFill("x1 vs anodeBack",600,-300,300,pevent.x1,512,0,4096,pevent.anodeBack);
}
if(!event.focalPlane.delayBL.empty() && !event.focalPlane.delayBR.empty()) {
pevent.fp2_tdiff = (event.focalPlane.delayBL[0].Time-event.focalPlane.delayBR[0].Time)*0.5;
pevent.fp2_tsum = (event.focalPlane.delayBL[0].Time+event.focalPlane.delayBR[0].Time);
pevent.fp2_tcheck = (pevent.fp2_tsum)/2.0-pevent.anodeBackTime;
pevent.delayBackMaxTime = max(event.focalPlane.delayBL[0].Time, event.focalPlane.delayBR[0].Time);
pevent.x2 = pevent.fp2_tdiff*1.0/1.98; //position from time, based on total delay
MyFill("x2",1200,-300,300,pevent.x2);
MyFill("x2 vs anodeBack",600,-300,300,pevent.x2,512,0,4096,pevent.anodeBack);
}
/*SABRE data*/
for(int j=0; j<5; j++) {
if(!event.sabreArray[j].rings.empty()) {
pevent.sabreRingE[j] = event.sabreArray[j].rings[0].Long;
pevent.sabreRingChannel[j] = event.sabreArray[j].rings[0].Ch;
pevent.sabreRingTime[j] = event.sabreArray[j].rings[0].Time;
}
if(!event.sabreArray[j].wedges.empty()) {
pevent.sabreWedgeE[j] = event.sabreArray[j].wedges[0].Long;
pevent.sabreWedgeChannel[j] = event.sabreArray[j].wedges[0].Ch;
pevent.sabreWedgeTime[j] = event.sabreArray[j].wedges[0].Time;
}
/*Aaaand passes on all of the rest. 4/24/20 GWM*/
pevent.sabreArray[j] = event.sabreArray[j];
}
/*Delay lines and all that*/
if(!event.focalPlane.delayFR.empty())
{
pevent.delayFrontRightE = event.focalPlane.delayFR[0].Long;
pevent.delayFrontRightTime = event.focalPlane.delayFR[0].Time;
pevent.delayFrontRightShort = event.focalPlane.delayFR[0].Short;
}
if(!event.focalPlane.delayFL.empty())
{
pevent.delayFrontLeftE = event.focalPlane.delayFL[0].Long;
pevent.delayFrontLeftTime = event.focalPlane.delayFL[0].Time;
pevent.delayFrontLeftShort = event.focalPlane.delayFL[0].Short;
}
if(!event.focalPlane.delayBR.empty())
{
pevent.delayBackRightE = event.focalPlane.delayBR[0].Long;
pevent.delayBackRightTime = event.focalPlane.delayBR[0].Time;
pevent.delayBackRightShort = event.focalPlane.delayBR[0].Short;
}
if(!event.focalPlane.delayBL.empty())
{
pevent.delayBackLeftE = event.focalPlane.delayBL[0].Long;
pevent.delayBackLeftTime = event.focalPlane.delayBL[0].Time;
pevent.delayBackLeftShort = event.focalPlane.delayBL[0].Short;
}
if(!event.focalPlane.delayFL.empty() && !event.focalPlane.delayFR.empty())
{
pevent.fp1_tdiff = (event.focalPlane.delayFL[0].Time-event.focalPlane.delayFR[0].Time)*0.5;
pevent.fp1_tsum = (event.focalPlane.delayFL[0].Time+event.focalPlane.delayFR[0].Time);
pevent.fp1_tcheck = (pevent.fp1_tsum)/2.0-pevent.anodeFrontTime;
pevent.delayFrontMaxTime = max(event.focalPlane.delayFL[0].Time, event.focalPlane.delayFR[0].Time);
pevent.x1 = pevent.fp1_tdiff*1.0/2.10; //position from time, based on total delay
MyFill("x1",1200,-300,300,pevent.x1);
MyFill("x1 vs anodeBack",600,-300,300,pevent.x1,512,0,4096,pevent.anodeBack);
}
if(!event.focalPlane.delayBL.empty() && !event.focalPlane.delayBR.empty())
{
pevent.fp2_tdiff = (event.focalPlane.delayBL[0].Time-event.focalPlane.delayBR[0].Time)*0.5;
pevent.fp2_tsum = (event.focalPlane.delayBL[0].Time+event.focalPlane.delayBR[0].Time);
pevent.fp2_tcheck = (pevent.fp2_tsum)/2.0-pevent.anodeBackTime;
pevent.delayBackMaxTime = max(event.focalPlane.delayBL[0].Time, event.focalPlane.delayBR[0].Time);
pevent.x2 = pevent.fp2_tdiff*1.0/1.98; //position from time, based on total delay
MyFill("x2",1200,-300,300,pevent.x2);
MyFill("x2 vs anodeBack",600,-300,300,pevent.x2,512,0,4096,pevent.anodeBack);
}
/*SABRE data*/
for(int j=0; j<5; j++)
{
if(!event.sabreArray[j].rings.empty())
{
pevent.sabreRingE[j] = event.sabreArray[j].rings[0].Long;
pevent.sabreRingChannel[j] = event.sabreArray[j].rings[0].Ch;
pevent.sabreRingTime[j] = event.sabreArray[j].rings[0].Time;
}
if(!event.sabreArray[j].wedges.empty())
{
pevent.sabreWedgeE[j] = event.sabreArray[j].wedges[0].Long;
pevent.sabreWedgeChannel[j] = event.sabreArray[j].wedges[0].Ch;
pevent.sabreWedgeTime[j] = event.sabreArray[j].wedges[0].Time;
}
/*Aaaand passes on all of the rest. 4/24/20 GWM*/
pevent.sabreArray[j] = event.sabreArray[j];
}
/*Make some histograms and xavg*/
MyFill("anodeBack vs scintLeft",512,0,4096,pevent.scintLeft,512,0,4096,pevent.anodeBack);
if(pevent.x1 != -1e6 && pevent.x2 != -1e6) {
pevent.xavg = pevent.x1*w1+pevent.x2*w2;
MyFill("xavg",1200,-300,300,pevent.xavg);
if((pevent.x2-pevent.x1) > 0) {
pevent.theta = atan((pevent.x2-pevent.x1)/36.0);
} else if((pevent.x2-pevent.x1) < 0) {
pevent.theta = TMath::Pi() + atan((pevent.x2-pevent.x1)/36.0);
} else {
pevent.theta = TMath::Pi()/2.0;
}
MyFill("xavg vs theta",600,-300,300,pevent.xavg,314,0,3.14,pevent.theta);
MyFill("x1 vs x2",600,-300,300,pevent.x1,600,-300,300,pevent.x2);
}
if(pevent.anodeFrontTime != -1 && pevent.scintRightTime != -1) {
pevent.fp1_y = pevent.anodeFrontTime-pevent.scintRightTime;
}
if(pevent.anodeBackTime != -1 && pevent.scintRightTime != -1) {
pevent.fp2_y = pevent.anodeBackTime-pevent.scintRightTime;
}
/*Make some histograms and xavg*/
MyFill("anodeBack vs scintLeft",512,0,4096,pevent.scintLeft,512,0,4096,pevent.anodeBack);
if(pevent.x1 != -1e6 && pevent.x2 != -1e6)
{
pevent.xavg = pevent.x1*w1+pevent.x2*w2;
MyFill("xavg",1200,-300,300,pevent.xavg);
if((pevent.x2-pevent.x1) > 0)
pevent.theta = std::atan((pevent.x2-pevent.x1)/36.0);
else if((pevent.x2-pevent.x1) < 0)
pevent.theta = TMath::Pi() + std::atan((pevent.x2-pevent.x1)/36.0);
else
pevent.theta = TMath::Pi()/2.0;
MyFill("xavg vs theta",600,-300,300,pevent.xavg,314,0,3.14,pevent.theta);
MyFill("x1 vs x2",600,-300,300,pevent.x1,600,-300,300,pevent.x2);
}
if(pevent.anodeFrontTime != -1 && pevent.scintRightTime != -1)
pevent.fp1_y = pevent.anodeFrontTime-pevent.scintRightTime;
if(pevent.anodeBackTime != -1 && pevent.scintRightTime != -1)
pevent.fp2_y = pevent.anodeBackTime-pevent.scintRightTime;
}
ProcessedEvent SFPAnalyzer::GetProcessedEvent(CoincEvent& event) {
AnalyzeEvent(event);
return pevent;
ProcessedEvent SFPAnalyzer::GetProcessedEvent(CoincEvent& event)
{
AnalyzeEvent(event);
return pevent;
}

497
src/evb/SFPPlotter.cpp
View File

@ -11,284 +11,289 @@
#include <TSystem.h>
/*Generates storage and initializes pointers*/
SFPPlotter::SFPPlotter() {
rootObj = new THashTable();
rootObj->SetOwner(false);//THashTable doesnt own members; avoid double delete
event_address = new ProcessedEvent();
chain = new TChain("SPSTree");
m_pb = NULL;
SFPPlotter::SFPPlotter()
{
rootObj = new THashTable();
rootObj->SetOwner(false);//THashTable doesnt own members; avoid double delete
event_address = new ProcessedEvent();
chain = new TChain("SPSTree");
m_pb = nullptr;
}
SFPPlotter::~SFPPlotter() {
delete event_address;
SFPPlotter::~SFPPlotter()
{
delete event_address;
}
/*2D histogram fill wrapper*/
void SFPPlotter::MyFill(const string& name, int binsx, double minx, double maxx, double valuex,
int binsy, double miny, double maxy, double valuey) {
TH2F *histo = (TH2F*) rootObj->FindObject(name.c_str());
if(histo != NULL) {
histo->Fill(valuex, valuey);
} else {
TH2F *h = new TH2F(name.c_str(), name.c_str(), binsx, minx, maxx, binsy, miny, maxy);
h->Fill(valuex, valuey);
rootObj->Add(h);
}
void SFPPlotter::MyFill(THashTable* table, const std::string& name, int binsx, double minx, double maxx, double valuex,
int binsy, double miny, double maxy, double valuey)
{
TH2F *histo = (TH2F*) rootObj->FindObject(name.c_str());
if(histo != nullptr)
histo->Fill(valuex, valuey);
else
{
TH2F *h = new TH2F(name.c_str(), name.c_str(), binsx, minx, maxx, binsy, miny, maxy);
h->Fill(valuex, valuey);
rootObj->Add(h);
}
}
/*1D histogram fill wrapper*/
void SFPPlotter::MyFill(const string& name, int binsx, double minx, double maxx, double valuex) {
TH1F *histo = (TH1F*) rootObj->FindObject(name.c_str());
if(histo != NULL) {
histo->Fill(valuex);
} else {
TH1F *h = new TH1F(name.c_str(), name.c_str(), binsx, minx, maxx);
h->Fill(valuex);
rootObj->Add(h);
}
}
void SFPPlotter::ApplyCutlist(const string& listname) {
cutter.SetCuts(listname);
void SFPPlotter::MyFill(THashTable* table, const std::string& name, int binsx, double minx, double maxx, double valuex)
{
TH1F *histo = (TH1F*) rootObj->FindObject(name.c_str());
if(histo != nullptr)
histo->Fill(valuex);
else
{
TH1F *h = new TH1F(name.c_str(), name.c_str(), binsx, minx, maxx);
h->Fill(valuex);
rootObj->Add(h);
}
}
/*Makes histograms where only rejection is unset data*/
void SFPPlotter::MakeUncutHistograms(ProcessedEvent ev) {
MyFill("x1NoCuts_bothplanes",600,-300,300,ev.x2);
MyFill("x2NoCuts_bothplanes",600,-300,300,ev.x2);
MyFill("xavgNoCuts_bothplanes",600,-300,300,ev.xavg);
MyFill("xavgNoCuts_theta_bothplanes",600,-300,300,ev.xavg,100,0,TMath::Pi()/2.,ev.theta);
MyFill("x1_delayBackRightE_NoCuts",600,-300,300,ev.x1,512,0,4096,ev.delayBackRightE);
MyFill("x2_delayBackRightE_NoCuts",600,-300,300,ev.x2,512,0,4096,ev.delayBackRightE);
MyFill("xavg_delayBackRightE_NoCuts",600,-300,300,ev.xavg,512,0,4096,ev.delayBackRightE);
MyFill("x1_x2_NoCuts",600,-300,300,ev.x1,600,-300,300,ev.x2);
void SFPPlotter::MakeUncutHistograms(const ProcessedEvent& ev, THashTable* table)
{
MyFill(table,"x1NoCuts_bothplanes",600,-300,300,ev.x2);
MyFill(table,"x2NoCuts_bothplanes",600,-300,300,ev.x2);
MyFill(table,"xavgNoCuts_bothplanes",600,-300,300,ev.xavg);
MyFill(table,"xavgNoCuts_theta_bothplanes",600,-300,300,ev.xavg,100,0,TMath::Pi()/2.,ev.theta);
MyFill(table,"x1_delayBackRightE_NoCuts",600,-300,300,ev.x1,512,0,4096,ev.delayBackRightE);
MyFill(table,"x2_delayBackRightE_NoCuts",600,-300,300,ev.x2,512,0,4096,ev.delayBackRightE);
MyFill(table,"xavg_delayBackRightE_NoCuts",600,-300,300,ev.xavg,512,0,4096,ev.delayBackRightE);
MyFill(table,"x1_x2_NoCuts",600,-300,300,ev.x1,600,-300,300,ev.x2);
Double_t delayBackAvgE = (ev.delayBackRightE+ev.delayBackLeftE)/2.0;
MyFill("x1_delayBackAvgE_NoCuts",600,-300,300,ev.x1,512,0,4096,delayBackAvgE);
MyFill("x2_delayBackAvgE_NoCuts",600,-300,300,ev.x2,512,0,4096,delayBackAvgE);
MyFill("xavg_delayBackAvgE_NoCuts",600,-300,300,ev.xavg,512,0,4096,delayBackAvgE);
Double_t delayFrontAvgE = (ev.delayFrontRightE+ev.delayFrontLeftE)/2.0;
MyFill("x1_delayFrontAvgE_NoCuts",600,-300,300,ev.x1,512,0,4096,delayFrontAvgE);
MyFill("x2_delayFrontAvgE_NoCuts",600,-300,300,ev.x2,512,0,4096,delayFrontAvgE);
MyFill("xavg_delayFrontAvgE_NoCuts",600,-300,300,ev.xavg,512,0,4096,delayFrontAvgE);
Double_t delayBackAvgE = (ev.delayBackRightE+ev.delayBackLeftE)/2.0;
MyFill(table,"x1_delayBackAvgE_NoCuts",600,-300,300,ev.x1,512,0,4096,delayBackAvgE);
MyFill(table,"x2_delayBackAvgE_NoCuts",600,-300,300,ev.x2,512,0,4096,delayBackAvgE);
MyFill(table,"xavg_delayBackAvgE_NoCuts",600,-300,300,ev.xavg,512,0,4096,delayBackAvgE);
Double_t delayFrontAvgE = (ev.delayFrontRightE+ev.delayFrontLeftE)/2.0;
MyFill(table,"x1_delayFrontAvgE_NoCuts",600,-300,300,ev.x1,512,0,4096,delayFrontAvgE);
MyFill(table,"x2_delayFrontAvgE_NoCuts",600,-300,300,ev.x2,512,0,4096,delayFrontAvgE);
MyFill(table,"xavg_delayFrontAvgE_NoCuts",600,-300,300,ev.xavg,512,0,4096,delayFrontAvgE);
MyFill("scintLeft_anodeBack_NoCuts",512,0,4096,ev.scintLeft,512,0,4096,ev.anodeBack);
MyFill("scintLeft_anodeFront_NoCuts",512,0,4096,ev.scintLeft,512,0,4096,ev.anodeFront);
MyFill("scintLeft_cathode_NoCuts",512,0,4096,ev.scintLeft,512,0,4096,ev.cathode);
MyFill(table,"scintLeft_anodeBack_NoCuts",512,0,4096,ev.scintLeft,512,0,4096,ev.anodeBack);
MyFill(table,"scintLeft_anodeFront_NoCuts",512,0,4096,ev.scintLeft,512,0,4096,ev.anodeFront);
MyFill(table,"scintLeft_cathode_NoCuts",512,0,4096,ev.scintLeft,512,0,4096,ev.cathode);
MyFill("x1_scintLeft_NoCuts",600,-300,300,ev.x1,512,0,4096,ev.scintLeft);
MyFill("x2_scintLeft_NoCuts",600,-300,300,ev.x2,512,0,4096,ev.scintLeft);
MyFill("xavg_scintLeft_NoCuts",600,-300,300,ev.xavg,512,0,4096,ev.scintLeft);
MyFill(table,"x1_scintLeft_NoCuts",600,-300,300,ev.x1,512,0,4096,ev.scintLeft);
MyFill(table,"x2_scintLeft_NoCuts",600,-300,300,ev.x2,512,0,4096,ev.scintLeft);
MyFill(table,"xavg_scintLeft_NoCuts",600,-300,300,ev.xavg,512,0,4096,ev.scintLeft);
MyFill("x1_anodeBack_NoCuts",600,-300,300,ev.x1,512,0,4096,ev.anodeBack);
MyFill("x2_anodeBack_NoCuts",600,-300,300,ev.x2,512,0,4096,ev.anodeBack);
MyFill("xavg_anodeBack_NoCuts",600,-300,300,ev.xavg,512,0,4096,ev.anodeBack);
MyFill(table,"x1_anodeBack_NoCuts",600,-300,300,ev.x1,512,0,4096,ev.anodeBack);
MyFill(table,"x2_anodeBack_NoCuts",600,-300,300,ev.x2,512,0,4096,ev.anodeBack);
MyFill(table,"xavg_anodeBack_NoCuts",600,-300,300,ev.xavg,512,0,4096,ev.anodeBack);
MyFill("x1_anodeFront_NoCuts",600,-300,300,ev.x1,512,0,4096,ev.anodeFront);
MyFill("x2_anodeFront_NoCuts",600,-300,300,ev.x2,512,0,4096,ev.anodeFront);
MyFill("xavg_anodeFront_NoCuts",600,-300,300,ev.xavg,512,0,4096,ev.anodeFront);
MyFill(table,"x1_anodeFront_NoCuts",600,-300,300,ev.x1,512,0,4096,ev.anodeFront);
MyFill(table,"x2_anodeFront_NoCuts",600,-300,300,ev.x2,512,0,4096,ev.anodeFront);
MyFill(table,"xavg_anodeFront_NoCuts",600,-300,300,ev.xavg,512,0,4096,ev.anodeFront);
MyFill("x1_cathode_NoCuts",600,-300,300,ev.x1,512,0,4096,ev.cathode);
MyFill("x2_cathode_NoCuts",600,-300,300,ev.x2,512,0,4096,ev.cathode);
MyFill("xavg_cathode_NoCuts",600,-300,300,ev.xavg,512,0,4096,ev.cathode);
MyFill(table,"x1_cathode_NoCuts",600,-300,300,ev.x1,512,0,4096,ev.cathode);
MyFill(table,"x2_cathode_NoCuts",600,-300,300,ev.x2,512,0,4096,ev.cathode);
MyFill(table,"xavg_cathode_NoCuts",600,-300,300,ev.xavg,512,0,4096,ev.cathode);
/****Timing relative to back anode****/
if(ev.anodeBackTime != -1 && ev.scintLeftTime != -1) {
Double_t anodeRelFT = ev.anodeFrontTime - ev.anodeBackTime;
Double_t delayRelFT = ev.delayFrontMaxTime - ev.anodeBackTime;
Double_t delayRelBT = ev.delayBackMaxTime - ev.anodeBackTime;
Double_t anodeRelBT = ev.anodeBackTime - ev.scintLeftTime;
Double_t delayRelFT_toScint = ev.delayFrontMaxTime - ev.scintLeftTime;
Double_t delayRelBT_toScint = ev.delayBackMaxTime - ev.scintLeftTime;
MyFill("anodeRelFrontTime_NoCuts",1000,-3000,3500, anodeRelFT);
MyFill("delayRelFrontTime_NoCuts",1000,-3000,-3500,delayRelFT);
MyFill("delayRelBackTime_NoCuts",1000,-3000,-3500,delayRelBT);
for(int i=0; i<5; i++) {
if(ev.sabreRingE[i] != -1) {
Double_t sabreRelRT = ev.sabreRingTime[i] - ev.anodeBackTime;
Double_t sabreRelWT = ev.sabreWedgeTime[i] - ev.anodeBackTime;
Double_t sabreRelRT_toScint = ev.sabreRingTime[i] - ev.scintLeftTime;
Double_t sabreRelWT_toScint = ev.sabreWedgeTime[i] - ev.scintLeftTime;
MyFill("xavg_sabrefcoinc_NoCuts",600,-300,300, ev.xavg);
MyFill("sabreRelRingTime_NoCuts",1000,-3000,3500, sabreRelRT);
MyFill("sabreRelWedgeTime_NoCuts",1000,-3000,3500, sabreRelWT);
MyFill("sabreRelRingTime_toScint",1000,-3000,3500,sabreRelRT_toScint);
MyFill("sabreRelWedgeTime_toScint",1000,-3000,3500,sabreRelWT_toScint);
MyFill("sabreRelRTScint_sabreRelRTAnode",500,-3000,3500,sabreRelRT_toScint,500,-3000,3500,sabreRelRT);
MyFill("sabreRelRTScint_sabreRingChannel",500,-3000,3500,sabreRelRT_toScint,144,0,144,ev.sabreRingChannel[i]);
MyFill("sabreRelRTAnode_sabreRingChannel",500,-3000,3500,sabreRelRT,144,0,144,ev.sabreRingChannel[i]);
MyFill("sabreRelWTScint_sabreWedgeChannel",500,-3000,3500,sabreRelWT_toScint,144,0,144,ev.sabreWedgeChannel[i]);
MyFill("sabreRelRT_sabreRelWT",500,-3000,3500,sabreRelRT,500,-3000,3500,sabreRelWT);
MyFill("sabreRelRT_sabreRelWT_scint",500,-3000,3500,sabreRelRT_toScint,500,-3000,3500,sabreRelWT_toScint);
MyFill("sabreRelRTScint_anodeRelT",500,-3000,3500,sabreRelRT_toScint,500,-3000,3500,anodeRelBT);
}
}
MyFill("anodeBackRelTime_toScint",1000,-3000,3500,anodeRelBT);
MyFill("delayRelBackTime_toScint",1000,-3000,3500,delayRelBT_toScint);
MyFill("delayRelFrontTime_toScint",1000,-3000,3500,delayRelFT_toScint);
} else {
MyFill("noscinttime_counter_NoCuts",2,0,1,1);
}
int count = 0;
for(int i=0; i<5; i++) {
if(ev.sabreRingE[i] != -1) { //Again, at this point front&back are required
MyFill("sabreRingE_NoCuts",2000,0,20,ev.sabreRingE[i]);
MyFill("sabreRingChannel_sabreRingE_NoCuts",144,0,144,ev.sabreRingChannel[i],200,0,20,ev.sabreRingE[i]);
MyFill("sabreWedgeE_NoCuts",2000,0,20,ev.sabreWedgeE[i]);
MyFill("sabreWedgeChannel_sabreWedgeE_NoCuts",144,0,144,ev.sabreWedgeChannel[i],200,0,20,ev.sabreWedgeE[i]);
} else {
count++;
}
}
if(count == 80) {
MyFill("xavg_bothplanes_sabreanticoinc_NoCuts",600,-300,300,ev.xavg);
}
if(ev.x1 != -1e6 && ev.x2 == -1e6) {
MyFill("x1NoCuts_only1plane",600,-300,300,ev.x1);
} else if(ev.x2 != -1e6 && ev.x1 == -1e6) {
MyFill("x2NoCuts_only1plane",600,-300,300,ev.x2);
} else if(ev.x1 == -1e6 && ev.x2 == -1e6) {
MyFill("nopos_counter",2,0,1,1);
}
/****Timing relative to back anode****/
if(ev.anodeBackTime != -1 && ev.scintLeftTime != -1)
{
Double_t anodeRelFT = ev.anodeFrontTime - ev.anodeBackTime;
Double_t delayRelFT = ev.delayFrontMaxTime - ev.anodeBackTime;
Double_t delayRelBT = ev.delayBackMaxTime - ev.anodeBackTime;
Double_t anodeRelBT = ev.anodeBackTime - ev.scintLeftTime;
Double_t delayRelFT_toScint = ev.delayFrontMaxTime - ev.scintLeftTime;
Double_t delayRelBT_toScint = ev.delayBackMaxTime - ev.scintLeftTime;
MyFill(table,"anodeRelFrontTime_NoCuts",1000,-3000,3500, anodeRelFT);
MyFill(table,"delayRelFrontTime_NoCuts",1000,-3000,-3500,delayRelFT);
MyFill(table,"delayRelBackTime_NoCuts",1000,-3000,-3500,delayRelBT);
for(int i=0; i<5; i++)
{
if(ev.sabreRingE[i] != -1)
{
Double_t sabreRelRT = ev.sabreRingTime[i] - ev.anodeBackTime;
Double_t sabreRelWT = ev.sabreWedgeTime[i] - ev.anodeBackTime;
Double_t sabreRelRT_toScint = ev.sabreRingTime[i] - ev.scintLeftTime;
Double_t sabreRelWT_toScint = ev.sabreWedgeTime[i] - ev.scintLeftTime;
MyFill(table,"xavg_sabrefcoinc_NoCuts",600,-300,300, ev.xavg);
MyFill(table,"sabreRelRingTime_NoCuts",1000,-3000,3500, sabreRelRT);
MyFill(table,"sabreRelWedgeTime_NoCuts",1000,-3000,3500, sabreRelWT);
MyFill(table,"sabreRelRingTime_toScint",1000,-3000,3500,sabreRelRT_toScint);
MyFill(table,"sabreRelWedgeTime_toScint",1000,-3000,3500,sabreRelWT_toScint);
MyFill(table,"sabreRelRTScint_sabreRelRTAnode",500,-3000,3500,sabreRelRT_toScint,500,-3000,3500,sabreRelRT);
MyFill(table,"sabreRelRTScint_sabreRingChannel",500,-3000,3500,sabreRelRT_toScint,144,0,144,ev.sabreRingChannel[i]);
MyFill(table,"sabreRelRTAnode_sabreRingChannel",500,-3000,3500,sabreRelRT,144,0,144,ev.sabreRingChannel[i]);
MyFill(table,"sabreRelWTScint_sabreWedgeChannel",500,-3000,3500,sabreRelWT_toScint,144,0,144,ev.sabreWedgeChannel[i]);
MyFill(table,"sabreRelRT_sabreRelWT",500,-3000,3500,sabreRelRT,500,-3000,3500,sabreRelWT);
MyFill(table,"sabreRelRT_sabreRelWT_scint",500,-3000,3500,sabreRelRT_toScint,500,-3000,3500,sabreRelWT_toScint);
MyFill(table,"sabreRelRTScint_anodeRelT",500,-3000,3500,sabreRelRT_toScint,500,-3000,3500,anodeRelBT);
}
}
MyFill(table,"anodeBackRelTime_toScint",1000,-3000,3500,anodeRelBT);
MyFill(table,"delayRelBackTime_toScint",1000,-3000,3500,delayRelBT_toScint);
MyFill(table,"delayRelFrontTime_toScint",1000,-3000,3500,delayRelFT_toScint);
}
else
MyFill(table,"noscinttime_counter_NoCuts",2,0,1,1);
for(int i=0; i<5; i++)
{
if(ev.sabreRingE[i] != -1) //Again, at this point front&back are required
{
MyFill(table,"sabreRingE_NoCuts",2000,0,20,ev.sabreRingE[i]);
MyFill(table,"sabreRingChannel_sabreRingE_NoCuts",144,0,144,ev.sabreRingChannel[i],200,0,20,ev.sabreRingE[i]);
MyFill(table,"sabreWedgeE_NoCuts",2000,0,20,ev.sabreWedgeE[i]);
MyFill(table,"sabreWedgeChannel_sabreWedgeE_NoCuts",144,0,144,ev.sabreWedgeChannel[i],200,0,20,ev.sabreWedgeE[i]);
}
}
if(ev.x1 != -1e6 && ev.x2 == -1e6)
MyFill(table,"x1NoCuts_only1plane",600,-300,300,ev.x1);
else if(ev.x2 != -1e6 && ev.x1 == -1e6)
MyFill(table,"x2NoCuts_only1plane",600,-300,300,ev.x2);
else if(ev.x1 == -1e6 && ev.x2 == -1e6)
MyFill(table,"nopos_counter",2,0,1,1);
}
/*Makes histograms with cuts & gates implemented*/
void SFPPlotter::MakeCutHistograms(ProcessedEvent ev) {
if(cutter.IsInside(&ev)) {
MyFill("x1_bothplanes_Cut",600,-300,300,ev.x1);
MyFill("x2_bothplanes_Cut",600,-300,300,ev.x2);
MyFill("xavg_bothplanes_Cut",600,-300,300,ev.xavg);
MyFill("x1_x2_Cut",600,-300,300,ev.x1, 600,-300,300,ev.x2);
MyFill("xavg_theta_Cut_bothplanes",600,-300,300,ev.xavg,100,0,TMath::Pi()/2.,ev.theta);
MyFill("x1_delayBackRightE_Cut",600,-300,300,ev.x1,512,0,4096,ev.delayBackRightE);
MyFill("x2_delayBackRightE_Cut",600,-300,300,ev.x2,512,0,4096,ev.delayBackRightE);
MyFill("xavg_delayBackRightE_Cut",600,-300,300,ev.xavg,512,0,4096,ev.delayBackRightE);
void SFPPlotter::MakeCutHistograms(const ProcessedEvent& ev, THashTable* table)
{
if(!cutter.IsInside(&ev))
return;
Double_t delayBackAvgE = (ev.delayBackRightE+ev.delayBackLeftE)/2.0;
MyFill("x1_delayBackAvgE_Cut",600,-300,300,ev.x1,512,0,4096,delayBackAvgE);
MyFill("x2_delayBackAvgE_Cut",600,-300,300,ev.x2,512,0,4096,delayBackAvgE);
MyFill("xavg_delayBackAvgE_Cut",600,-300,300,ev.xavg,512,0,4096,delayBackAvgE);
Double_t delayFrontAvgE = (ev.delayFrontRightE+ev.delayFrontLeftE)/2.0;
MyFill("x1_delayFrontAvgE_Cut",600,-300,300,ev.x1,512,0,4096,delayFrontAvgE);
MyFill("x2_delayFrontAvgE_Cut",600,-300,300,ev.x2,512,0,4096,delayFrontAvgE);
MyFill("xavg_delayFrontAvgE_Cut",600,-300,300,ev.xavg,512,0,4096,delayFrontAvgE);
MyFill(table,"x1_bothplanes_Cut",600,-300,300,ev.x1);
MyFill(table,"x2_bothplanes_Cut",600,-300,300,ev.x2);
MyFill(table,"xavg_bothplanes_Cut",600,-300,300,ev.xavg);
MyFill(table,"x1_x2_Cut",600,-300,300,ev.x1, 600,-300,300,ev.x2);
MyFill(table,"xavg_theta_Cut_bothplanes",600,-300,300,ev.xavg,100,0,TMath::Pi()/2.,ev.theta);
MyFill(table,"x1_delayBackRightE_Cut",600,-300,300,ev.x1,512,0,4096,ev.delayBackRightE);
MyFill(table,"x2_delayBackRightE_Cut",600,-300,300,ev.x2,512,0,4096,ev.delayBackRightE);
MyFill(table,"xavg_delayBackRightE_Cut",600,-300,300,ev.xavg,512,0,4096,ev.delayBackRightE);
MyFill("scintLeft_anodeBack_Cut",512,0,4096,ev.scintLeft,512,0,4096,ev.anodeBack);
MyFill("scintLeft_anodeFront_Cut",512,0,4096,ev.scintLeft,512,0,4096,ev.anodeFront);
MyFill("scintLeft_cathode_Cut",512,0,4096,ev.scintLeft,512,0,4096,ev.cathode);
Double_t delayBackAvgE = (ev.delayBackRightE+ev.delayBackLeftE)/2.0;
MyFill(table,"x1_delayBackAvgE_Cut",600,-300,300,ev.x1,512,0,4096,delayBackAvgE);
MyFill(table,"x2_delayBackAvgE_Cut",600,-300,300,ev.x2,512,0,4096,delayBackAvgE);
MyFill(table,"xavg_delayBackAvgE_Cut",600,-300,300,ev.xavg,512,0,4096,delayBackAvgE);
Double_t delayFrontAvgE = (ev.delayFrontRightE+ev.delayFrontLeftE)/2.0;
MyFill(table,"x1_delayFrontAvgE_Cut",600,-300,300,ev.x1,512,0,4096,delayFrontAvgE);
MyFill(table,"x2_delayFrontAvgE_Cut",600,-300,300,ev.x2,512,0,4096,delayFrontAvgE);
MyFill(table,"xavg_delayFrontAvgE_Cut",600,-300,300,ev.xavg,512,0,4096,delayFrontAvgE);
MyFill("x1_scintLeft_Cut",600,-300,300,ev.x1,512,0,4096,ev.scintLeft);
MyFill("x2_scintLeft_Cut",600,-300,300,ev.x2,512,0,4096,ev.scintLeft);
MyFill("xavg_scintLeft_Cut",600,-300,300,ev.xavg,512,0,4096,ev.scintLeft);
MyFill(table,"scintLeft_anodeBack_Cut",512,0,4096,ev.scintLeft,512,0,4096,ev.anodeBack);
MyFill(table,"scintLeft_anodeFront_Cut",512,0,4096,ev.scintLeft,512,0,4096,ev.anodeFront);
MyFill(table,"scintLeft_cathode_Cut",512,0,4096,ev.scintLeft,512,0,4096,ev.cathode);
MyFill("x1_anodeBack_Cut",600,-300,300,ev.x1,512,0,4096,ev.anodeBack);
MyFill("x2_anodeBack_Cut",600,-300,300,ev.x2,512,0,4096,ev.anodeBack);
MyFill("xavg_anodeBack_Cut",600,-300,300,ev.xavg,512,0,4096,ev.anodeBack);
MyFill("x1_anodeFront_Cut",600,-300,300,ev.x1,512,0,4096,ev.anodeFront);
MyFill("x2_anodeFront_Cut",600,-300,300,ev.x2,512,0,4096,ev.anodeFront);
MyFill("xavg_anodeFront_Cut",600,-300,300,ev.xavg,512,0,4096,ev.anodeFront);
MyFill("x1_cathode_Cut",600,-300,300,ev.x1,512,0,4096,ev.cathode);
MyFill("x2_cathode_Cut",600,-300,300,ev.x2,512,0,4096,ev.cathode);
MyFill("xavg_cathode_Cut",600,-300,300,ev.xavg,512,0,4096,ev.cathode);
MyFill(table,"x1_scintLeft_Cut",600,-300,300,ev.x1,512,0,4096,ev.scintLeft);
MyFill(table,"x2_scintLeft_Cut",600,-300,300,ev.x2,512,0,4096,ev.scintLeft);
MyFill(table,"xavg_scintLeft_Cut",600,-300,300,ev.xavg,512,0,4096,ev.scintLeft);
/****Timing relative to back anode****/
if(ev.anodeBackTime != -1 && ev.scintLeftTime != -1) {
Double_t anodeRelFT = ev.anodeFrontTime - ev.anodeBackTime;
Double_t anodeRelBT = ev.anodeBackTime - ev.anodeBackTime;
Double_t anodeRelFT_toScint = ev.anodeFrontTime-ev.scintLeftTime;
MyFill("anodeRelBackTime_Cut",1000,-3000,3500, anodeRelBT);
MyFill("anodeRelFrontTime_Cut",1000,-3000,3500, anodeRelFT);
MyFill("anodeRelTime_toScint_Cut",1000,-3000,3500,anodeRelFT_toScint);
for(int i=0; i<5; i++) {
if(ev.sabreRingE[i] != -1) {
Double_t sabreRelRT = ev.sabreRingTime[i] - ev.anodeBackTime;
Double_t sabreRelWT = ev.sabreWedgeTime[i] - ev.anodeBackTime;
MyFill("sabreRelRingTime_Cut",1000,-3000,3500, sabreRelRT);
MyFill("sabreRelWedgeTime_Cut",1000,-3000,3500, sabreRelWT);
}
}
} else {
MyFill("noscinttime_counter_Cut",2,0,1,1);
}
int count = 0;
for(int i=0; i<5; i++) {
if(ev.sabreRingE[i] != -1) {
MyFill("sabreRingE_Cut",2000,0,20,ev.sabreRingE[i]);
MyFill("xavg_Cut_sabrefcoinc",600,-300,300,ev.xavg);
MyFill("xavg_sabreRingE_Cut",600,-300,300,ev.xavg,200,0,20,ev.sabreRingE[i]);
MyFill("sabreWedgeE_Cut",2000,0,20,ev.sabreWedgeE[i]);
MyFill("xavg_sabreWedgeE_Cut",600,-300,300,ev.xavg,200,0,20,ev.sabreWedgeE[i]);
} else {
count++;
}
}
if(count == 80) {
MyFill("xavg_bothplanes_sabreanticoinc_Cut",600,-300,300,ev.xavg);
}
}
MyFill(table,"x1_anodeBack_Cut",600,-300,300,ev.x1,512,0,4096,ev.anodeBack);
MyFill(table,"x2_anodeBack_Cut",600,-300,300,ev.x2,512,0,4096,ev.anodeBack);
MyFill(table,"xavg_anodeBack_Cut",600,-300,300,ev.xavg,512,0,4096,ev.anodeBack);
MyFill(table,"x1_anodeFront_Cut",600,-300,300,ev.x1,512,0,4096,ev.anodeFront);
MyFill(table,"x2_anodeFront_Cut",600,-300,300,ev.x2,512,0,4096,ev.anodeFront);
MyFill(table,"xavg_anodeFront_Cut",600,-300,300,ev.xavg,512,0,4096,ev.anodeFront);
MyFill(table,"x1_cathode_Cut",600,-300,300,ev.x1,512,0,4096,ev.cathode);
MyFill(table,"x2_cathode_Cut",600,-300,300,ev.x2,512,0,4096,ev.cathode);
MyFill(table,"xavg_cathode_Cut",600,-300,300,ev.xavg,512,0,4096,ev.cathode);
/****Timing relative to back anode****/
if(ev.anodeBackTime != -1 && ev.scintLeftTime != -1)
{
Double_t anodeRelFT = ev.anodeFrontTime - ev.anodeBackTime;
Double_t anodeRelBT = ev.anodeBackTime - ev.anodeBackTime;
Double_t anodeRelFT_toScint = ev.anodeFrontTime-ev.scintLeftTime;
MyFill(table,"anodeRelBackTime_Cut",1000,-3000,3500, anodeRelBT);
MyFill(table,"anodeRelFrontTime_Cut",1000,-3000,3500, anodeRelFT);
MyFill(table,"anodeRelTime_toScint_Cut",1000,-3000,3500,anodeRelFT_toScint);
for(int i=0; i<5; i++)
{
if(ev.sabreRingE[i] != -1)
{
Double_t sabreRelRT = ev.sabreRingTime[i] - ev.anodeBackTime;
Double_t sabreRelWT = ev.sabreWedgeTime[i] - ev.anodeBackTime;
MyFill(table,"sabreRelRingTime_Cut",1000,-3000,3500, sabreRelRT);
MyFill(table,"sabreRelWedgeTime_Cut",1000,-3000,3500, sabreRelWT);
}
}
}
else
{
MyFill(table,"noscinttime_counter_Cut",2,0,1,1);
}
for(int i=0; i<5; i++)
{
if(ev.sabreRingE[i] != -1)
{
MyFill(table,"sabreRingE_Cut",2000,0,20,ev.sabreRingE[i]);
MyFill(table,"xavg_Cut_sabrefcoinc",600,-300,300,ev.xavg);
MyFill(table,"xavg_sabreRingE_Cut",600,-300,300,ev.xavg,200,0,20,ev.sabreRingE[i]);
MyFill(table,"sabreWedgeE_Cut",2000,0,20,ev.sabreWedgeE[i]);
MyFill(table,"xavg_sabreWedgeE_Cut",600,-300,300,ev.xavg,200,0,20,ev.sabreWedgeE[i]);
}
}
}
/*Runs a list of files given from a RunMusher/Collector class*/
void SFPPlotter::Run(vector<TString> files, const string& output) {
Chain(files);
chain->SetBranchAddress("event", &event_address);
TFile *outfile = new TFile(output.c_str(), "RECREATE");
/*Runs a list of files given from a RunCollector class*/
void SFPPlotter::Run(const std::vector<std::string>& files, const std::string& output)
{
TFile *outfile = TFile::Open(output.c_str(), "RECREATE");
TChain* chain = new TChain("SPSTree");
for(unsigned int i=0; i<files.size(); i++)
chain->Add(files[i].c_str());
chain->SetBranchAddress("event", &event_address);
THastTable* table = new THashTable();
long blentries = chain->GetEntries();
if(m_pb) SetProgressBar(blentries);
cout<<"Total number of events: "<<blentries<<endl;
long blentries = chain->GetEntries();
if(m_pb)
SetProgressBar(blentries);
cout<<"Total number of events: "<<blentries<<endl;
long count=0, flush=blentries*0.01, nflushes=0;
if(flush == 0) flush = 1;
long count=0, flush_val=blentries*0.01, flush_count=0;
for(long double i=0; i<chain->GetEntries(); i++) {
count++;
if(count == flush) {
if(m_pb) {
m_pb->Increment(count);
gSystem->ProcessEvents();
count = 0;
} else {
nflushes++;
count=0;
std::cout<<"\rPercent of data processed: "<<nflushes*10<<"%"<<std::flush;
}
}
chain->GetEntry(i);
MakeUncutHistograms(*event_address);
if(cutter.IsValid()) MakeCutHistograms(*event_address);
}
cout<<endl;
outfile->cd();
rootObj->Write();
if(cutter.IsValid()) {
auto clist = cutter.GetCuts();
for(unsigned int i=0; i<clist.size(); i++) {
clist[i]->Write();
}
}
delete rootObj;
outfile->Close();
delete outfile;
}
/*Link all files*/
void SFPPlotter::Chain(vector<TString> files) {
for(unsigned int i=0; i<files.size(); i++) {
chain->Add(files[i]);
}
for(long i=0; i<chain->GetEntries(); i++)
{
count++;
if(count == flush_val)
{
if(m_pb) {
m_pb->Increment(count);
gSystem->ProcessEvents();
count = 0;
} else {
flush_count++;
count=0;
std::cout<<"\rPercent of data processed: "<<flush_count*10<<"%"<<std::flush;
}
}
chain->GetEntry(i);
MakeUncutHistograms(*event_address);
if(cutter.IsValid()) MakeCutHistograms(*event_address);
}
cout<<endl;
outfile->cd();
rootObj->Write();
if(cutter.IsValid())
{
auto clist = cutter.GetCuts();
for(unsigned int i=0; i<clist.size(); i++)
clist[i]->Write();
}
delete rootObj;
outfile->Close();
delete outfile;
}
void SFPPlotter::SetProgressBar(long total) {
m_pb->SetMax(total);
m_pb->SetMin(0);
m_pb->SetPosition(0);
gSystem->ProcessEvents();
void SFPPlotter::SetProgressBar(long total)
{
m_pb->SetMax(total);
m_pb->SetMin(0);
m_pb->SetPosition(0);
gSystem->ProcessEvents();
}

44
src/evb/ShiftMap.cpp
View File

@ -13,57 +13,69 @@
#include "ShiftMap.h"
ShiftMap::ShiftMap() :
m_filename(""), is_set(false)
m_filename(""), m_validFlag(false)
{
}
ShiftMap::ShiftMap(const std::string& filename) :
m_filename(filename), is_set(false)
m_filename(filename), m_validFlag(false)
{
ParseFile();
}
ShiftMap::~ShiftMap() {}
void ShiftMap::SetFile(const std::string& filename) {
void ShiftMap::SetFile(const std::string& filename)
{
m_filename = filename;
ParseFile();
}
Long64_t ShiftMap::GetShift(int gchan) {
if(!is_set) return 0.0;
uint64_t ShiftMap::GetShift(int gchan)
{
if(!m_validFlag)
return 0;
auto iter = m_map.find(gchan);
if(iter == m_map.end()) {
return 0.0;
} else return iter->second;
if(iter == m_map.end())
return 0;
else
return iter->second;
}
void ShiftMap::ParseFile() {
void ShiftMap::ParseFile()
{
m_validFlag = false;
std::ifstream input(m_filename);
if(!input.is_open()) return;
if(!input.is_open())
return;
int board, channel, gchan;
Long64_t shift;
uint64_t shift;
std::string junk, temp;
std::getline(input, junk);
std::getline(input, junk);
while(input>>board) {
while(input>>board)
{
input>>temp;
input>>shift;
if(temp == "all") { //keyword to set all channels in this board to same shift
for(int i=0; i<16; i++) {
if(temp == "all") //keyword to set all channels in this board to same shift
{
for(int i=0; i<16; i++)
{
gchan = board*16 + i;
m_map[gchan] = shift;
}
} else {
}
else
{
channel = stoi(temp);
gchan = channel + board*16;
m_map[gchan] = shift;
}
}
is_set = true;
m_validFlag = true;
}

15
src/evb/Stopwatch.cpp
View File

@ -8,25 +8,30 @@
#include "EventBuilder.h"
#include "Stopwatch.h"
Stopwatch::Stopwatch() {
Stopwatch::Stopwatch()
{
start_time = Clock::now();
stop_time = start_time;
}
Stopwatch::~Stopwatch() {}
void Stopwatch::Start() {
void Stopwatch::Start()
{
start_time = Clock::now();
}
void Stopwatch::Stop() {
void Stopwatch::Stop()
{
stop_time = Clock::now();
}
double Stopwatch::GetElapsedSeconds() {
double Stopwatch::GetElapsedSeconds()
{
return std::chrono::duration_cast<std::chrono::duration<double>>(stop_time-start_time).count();
}
double Stopwatch::GetElapsedMilliseconds() {
double Stopwatch::GetElapsedMilliseconds()
{
return std::chrono::duration_cast<std::chrono::duration<double>>(stop_time-start_time).count()*1000.0;
}

52
src/main.cpp
View File

@ -1,9 +1,11 @@
#include "EventBuilder.h"
#include "GWMEventBuilder.h"
#include "EVBApp.h"
#include "Stopwatch.h"
int main(int argc, char** argv) {
if(argc != 3) {
int main(int argc, char** argv)
{
if(argc != 3)
{
std::cerr<<"Incorrect number of command line arguments!"<<std::endl;
std::cerr<<"Need to specify type of operation (buildSlow, buildFast, etc.) and input file."<<std::endl;
return 1;
@ -14,46 +16,40 @@ int main(int argc, char** argv) {
/* DEFAULT Operation Types:
convert (convert binary archive to root data)
convertSlow (convert binary arhcive to event slow data)
convertFast (convert binary archive to event fast data)
convertSlowA (convert binary archive to analyzed slow event data)
convertFastA (convert binary archive to analyzed fast event data)
merge (combine root files)
plot (generate a default histogram file from analyzed data)
Convert (convert binary archive to root data)
ConvertSlow (convert binary arhcive to event slow data)
ConvertFast (convert binary archive to event fast data)
ConvertSlowA (convert binary archive to analyzed slow event data)
ConvertFastA (convert binary archive to analyzed fast event data)
Merge (combine root files)
Plot (generate a default histogram file from analyzed data)
*/
GWMEventBuilder theBuilder;
if(!theBuilder.ReadConfigFile(filename)) {
EVBApp theBuilder;
return 1;
}
Stopwatch timer;
timer.Start();
if(operation == "convert") {
theBuilder.SetAnalysisType(GWMEventBuilder::CONVERT);
if(operation == "Convert")
theBuilder.Convert2RawRoot();
} else if(operation == "merge") {
theBuilder.SetAnalysisType(GWMEventBuilder::MERGE);
else if(operation == "Merge")
theBuilder.MergeROOTFiles();
} else if(operation == "plot") {
theBuilder.SetAnalysisType(GWMEventBuilder::PLOT);
else if(operation == "Plot")
theBuilder.PlotHistograms();
} else if (operation == "convertSlow"){
theBuilder.SetAnalysisType(GWMEventBuilder::CONVERT_S);
else if (operation == "ConvertSlow")
theBuilder.Convert2SortedRoot();
} else if (operation == "convertFast"){
theBuilder.SetAnalysisType(GWMEventBuilder::CONVERT_F);
else if (operation == "ConvertFast")
theBuilder.Convert2FastSortedRoot();
} else if (operation == "convertSlowA"){
theBuilder.SetAnalysisType(GWMEventBuilder::CONVERT_SA);
else if (operation == "ConvertSlowA")
theBuilder.Convert2SlowAnalyzedRoot();
} else if (operation == "convertFastA"){
theBuilder.SetAnalysisType(GWMEventBuilder::CONVERT_FA);
else if (operation == "ConvertFastA")
theBuilder.Convert2FastAnalyzedRoot();
} else {
else
{
std::cerr<<"Unidentified type of operation! Check your first argument."<<std::endl;
return 1;
}
timer.Stop();
std::cout<<"Elapsed time (ms): "<<timer.GetElapsedMilliseconds()<<std::endl;