Commit 002aafbf authored by Riccardo Longo's avatar Riccardo Longo
Browse files
parents b734dd89 a4ed697f
......@@ -32,7 +32,10 @@ class Channel {
* [4,4] [4,3] [4,2] [4,1]
* 2021 mapping will come with the dedicated classes
*/
/** Type of detector - ZDC or RPD **/
/*
* Hardware channel information
*/
/** Type of detector - ZDC or RPD **/
std::string detector;
/** Channel name - CXX with XX in [1,20] */
std::string name;
......@@ -52,6 +55,11 @@ class Channel {
int Vop;
/** PMT reference code **/
std::string PMTcode;
/** Polarity of the channel */
bool pos_polarity;
/*
* Waveform storage
*/
/** Raw waveform for a particular event */
std::vector < float > WF;
/** First derivative of the waveform for a particular event */
......@@ -68,50 +76,54 @@ class Channel {
TH1D* FirstDerivative;
/** Fitting function for the PMT pulse **/
TF1* FitFunc;
/** RMS value of the first derrivative of the waveform **/
double FirstDerivativeRMS;
/** Bin number of the peak center*/
int Peak_center;
/** Bin number of the derivative peak */
int Diff_Peak_center;
/** Calibrated time of the peak center in ns*/
double Peak_time;
/** Calibrated time of max slope in ns*/
double Diff_Peak_time;
/** Height of the peak */
double Peak_max;
/** Max value of the derivative */
double Diff_max;
/** Bin number of 1/3 peak value on the rising edge*/
int DiscriminatorTime;
/** Start and end of the hit window*/
std::pair< int, int > hit_window;
/** Flag if the channel was hit */
bool was_hit;
/** Flag if the waveform saturated */
bool saturated;
/** Polarity of the channel */
bool pos_polarity;
/** Flag if the waveform has a DRS4 spike present */
bool spike;
/*
* Waveform processing variables
*/
/** Flag if the low pass filtering is to be performed on this channel */
bool doLPF;
/** Flag if DRS4 non-linearity compensation is to be performed on this channel */
bool DRS4NLcomp;
/** Pedestal of raw waveform */
double PedMean;
/** Pedestal RMS of raw waveform */
double PedRMS;
/** Integral of PWF_histo in the hit window in pC */
double Charge;
/** Flag if median filtering is to be performed on this channel */
bool doMedianFilter;
/** Hit threshold multiplier for this channel */
double Threshold;
/** Smoothing factor (number of bins to be summed) in calculation of waveform first derivative */
int diffSmoothing;
/** Low pass filter frequency */
int LPFfreq;
/** Crossing zero points - dummy - to be checked by Sheng **/
std::vector < int > CrossZeroPoints;
/** Median filtering half window size */
int medFiltHalfWindow;
/*
* Variables produced by WFAnalysis
*/
/** Flag if the channel was hit */
bool was_hit;
/** Number of hits detected in a single event */
int nHits;
/** Pair of vectors containing the start and end of each hit window*/
std::pair< std::vector< int >, std::vector< int > > hit_window;
/** RMS value of the first derrivative of the waveform **/
double FirstDerivativeRMS;
/** Vector of bin number of the peak location*/
std::vector< int > Peak_center;
/** Vector of bin number of the derivative peaks */
std::vector< int > Diff_Peak_center;
/** Vector of calibrated time of the peak in ns*/
std::vector< double > Peak_time;
/** Vector of calibrated time of max slope for each hit in ns*/
std::vector< double > Diff_Peak_time;
/** Vector of height of the peaks */
std::vector< double > Peak_max;
/** Vector of max value of the derivative in each hit window */
std::vector< double > Diff_max;
/** Vector of integral of PWF_histo in each hit window in pC */
std::vector< double > Charge;
/** Pedestal of raw waveform */
double PedMean;
/** Pedestal RMS of raw waveform */
double PedRMS;
/** Flag if the waveform saturated */
bool saturated;
};
......
......@@ -50,8 +50,8 @@ class DataReader2021{
void ReadListOfFiles( std::string listname );
void LoadAlignmentFile (std::string _inFile = std::getenv("JZCaPA") + std::string("/Utils/Alignment_2021.xml"));
void LoadConfigurationFile (std::string _inFile = std::getenv("JZCaPA") + std::string("/Utils/ConfigFile2021.xml"));
void LoadAlignmentFile (std::string _inFile = "");
void LoadConfigurationFile (std::string _inFile = "");
void LoadTimingFile (std::string _inFile = "" );
void SetDebugMode ( ) { m_debug = true; }
void SetVerbosity ( int _level ){ m_verbose = _level; }
......
......@@ -15,6 +15,7 @@
#include <string>
#include "Containers.h"
#include "Visualizer.h"
#include "TTree.h"
class Detector{
......@@ -33,20 +34,24 @@ class Detector{
virtual Alignment* GetAlignment( ) { return m_Alignment; }
virtual Alignment2021* GetAlignment2021( ) { return m_Alignment2021; }
virtual void SetNSamples ( int _nSamples ) { m_nSamp = _nSamples; }
virtual void SetElement ( Channel* _entry) { m_Element.push_back(_entry); }
virtual void SetPosition (double x, double y, double z) { m_Position[0] = x; m_Position[1] = y; m_Position[2] = z; }
virtual void SetAngle (double _cosx = 0, double _cosy = 0, double _cosz = 0) { m_Angle[0] = _cosx; m_Angle[1] = _cosy; m_Angle[2] = _cosz; }
virtual void SetBranches ( TTree* _dataTree );
virtual void SetAlignment( Alignment* _alignment ){ m_Alignment = _alignment; }
virtual void SetAlignment2021( Alignment2021* _alignment ){ m_Alignment2021 = _alignment; }
virtual void SetWFAthresholds( double _threshold, int _diffSmoothing ){ m_Threshold = _threshold; m_diffSmoothing = _diffSmoothing; }
virtual void SetFFTlowPass( int _freq ){ m_LPFfreq = _freq; m_doLPF = true; }
virtual void DoDRS4NLcomp( ){ m_doDRS4NLcomp = true; }
virtual void DeclareHistograms ( );
virtual void FillHistograms ( );
virtual void SetNSamples ( int _nSamples ) { m_nSamp = _nSamples; }
virtual void SetElement ( Channel* _entry) { m_Element.push_back(_entry); }
virtual void SetPosition ( double x, double y, double z) { m_Position[0] = x; m_Position[1] = y; m_Position[2] = z; }
virtual void SetAngle ( double _cosx = 0, double _cosy = 0, double _cosz = 0) { m_Angle[0] = _cosx; m_Angle[1] = _cosy; m_Angle[2] = _cosz; }
virtual void SetBranches ( TTree* _dataTree );
virtual void SetAlignment ( Alignment* _alignment ){ m_Alignment = _alignment; }
virtual void SetAlignment2021 ( Alignment2021* _alignment ){ m_Alignment2021 = _alignment; }
virtual void SetWFAthresholds ( double _threshold, int _diffSmoothing ){ m_Threshold = _threshold; m_diffSmoothing = _diffSmoothing; }
virtual void SetFFTlowPass ( int _freq ){ m_LPFfreq = _freq; m_doLPF = true; }
virtual void SetMedianFilter ( int _window ){ m_medHW = _window; m_doMedianFilter = true; }
virtual void DoDRS4NLcomp ( ){ m_doDRS4NLcomp = true; }
virtual void DeclareHistograms( );
virtual void FillHistograms ( );
virtual void PrintMap ( ) = 0;
virtual void PrintMap ( ) = 0;
virtual void PlotWF ( int eventNum );
virtual void PlotWFandDerivative( int eventNum, double _yMin1=0, double _yMax1=0, double _yMin2=0, double _yMax2=0 );
virtual void PlotWFandPWF ( int eventNum, double _yMin1=0, double _yMax1=0, double _yMin2=0, double _yMax2=0 );
void SetName ( std::string _name ) { m_name = _name; }
std::string GetName ( ) { return m_name; }
......@@ -56,6 +61,16 @@ class Detector{
std::string m_name;
/** Vector of channels associated to the dector **/
std::vector< Channel* > m_Element;
/** 1d vector of waveform histograms used with Visualizer for low level analysis */
std::vector < TH1* > hWFvec;
/** 1d vector of processed waveform histograms used with Visualizer for low level analysis */
std::vector < TH1* > hPWFvec;
/** 1d vector of waveform derivative histograms used with Visualizer for low level analysis */
std::vector < TH1* > hDiffvec;
/** 2d vector of TMarkers or TLines to be drawn on the pads. 1st index is channel number, 2nd index is mark to be drawn */
std::vector< std::vector < TObject* >* >* hMarksFront;
/** 2d vector of TMarkers or TLines to be drawn on the pads. 1st index is channel number, 2nd index is mark to be drawn */
std::vector< std::vector < TObject* >* >* hMarksBack;
/** Three element array with x, y, and z of some pre-defined point on the detector **/
double m_Position[3];
/** Three element array of angle about the x, y, and z axis **/
......@@ -66,12 +81,18 @@ class Detector{
int m_diffSmoothing = 25;
/** Frequency cutoff for low pass filtering of waveforms if selected */
int m_LPFfreq = 50;
/** Half window size for median filtering */
int m_medHW = 0;
/** Flag if low pass filtering is to be performed on this detector's waveforms */
bool m_doLPF = false;
/** Perform DRS4 non-linearity compensation on channels of this detector */
bool m_doDRS4NLcomp = false;
/** Perform median filtering on channels of this detector */
bool m_doMedianFilter = false;
/** Number of samples per channel **/
int m_nSamp = 1024;
/** Visualizer for plots **/
Visualizer* m_viz = 0;
/** Alignment of the 2018 Testbeam */
Alignment* m_Alignment = 0;
/** Alignment of the 2021 Testbeam */
......
......@@ -41,15 +41,15 @@ class WFAnalysis : public Analysis{
virtual void AnalyzeEvent ( const std::vector< TH1* >& );
virtual void AnalyzeEvent ( const std::vector< std::vector< float > >& );
virtual void AnalyzeEvent ( const std::vector< Channel* > vCh );
virtual int GetMaximumBin ( TH1* h, int _start, int _end );
virtual double GetCharge ( Channel* ch, int start, int end );
virtual void GetDifferential ( Channel* Ch );
virtual double GetRMS ( Channel* Ch );
virtual void GetPedestal ( Channel* ch );
virtual void SubtractPedestal( Channel* ch );
virtual void GetCharge ( Channel* ch );
virtual void FindHitWindow ( Channel* ch );
virtual void FindHitWindowNew( Channel* ch );
virtual void ZeroSuppress ( Channel* ch );
virtual void SpikeSupress ( Channel* ch );
virtual void MedianFilter ( Channel* ch );
virtual void LowPassFilter ( Channel* ch, TH1D* hIn = 0 );
virtual void DRS4Cal ( Channel* ch );
......
......@@ -170,7 +170,7 @@ void DataReader2021::ReadListOfFiles( std::string listname ){
* @param _inFile
*/
void DataReader2021::LoadAlignmentFile(std::string _inFile ){
if( _inFile == "" ) _inFile = std::getenv("JZCaPA") + std::string("/Utils/Alignment_2021.xml");
m_XMLparser = new XMLSettingsReader();
if (!m_XMLparser->parseFile(_inFile)) {
......@@ -215,7 +215,7 @@ void DataReader2021::LoadAlignmentFile(std::string _inFile ){
* @param _inFile
*/
void DataReader2021::LoadConfigurationFile(std::string _inFile ){
if( _inFile == "" ) _inFile = std::getenv("JZCaPA") + std::string("/Utils/ConfigFile2021.xml");
m_XMLparser = new XMLSettingsReader();
if (!m_XMLparser->parseFile(_inFile)) {
......
......@@ -12,6 +12,11 @@
#include "Detector.h"
#include "Containers.h"
#include "TLine.h"
#include "TMarker.h"
#include "TSystem.h"
#include <vector>
......@@ -82,6 +87,8 @@ void Detector::SetBranches( TTree *_dataTree ){
* @brief Declare histograms to be filled with the raw waveform
*/
void Detector::DeclareHistograms(){
hMarksFront = new std::vector< std::vector< TObject* >* >;
hMarksBack = new std::vector< std::vector< TObject* >* >;
for( uint ch = 0; ch < m_Element.size(); ch++ ){
m_Element.at(ch)->WF_histo = new TH1D( m_Element.at(ch)->name.c_str(), (m_Element.at(ch)->name + ", " + m_Element.at(ch)->detector).c_str(), m_nSamp, 0, m_nSamp);
......@@ -97,6 +104,15 @@ void Detector::DeclareHistograms(){
m_Element.at(ch)->doLPF = m_doLPF;
m_Element.at(ch)->LPFfreq = m_LPFfreq;
m_Element.at(ch)->DRS4NLcomp = m_doDRS4NLcomp;
m_Element.at(ch)->doMedianFilter = m_doMedianFilter;
m_Element.at(ch)->medFiltHalfWindow = m_medHW;
hWFvec.push_back(m_Element.at(ch)->WF_histo);
hPWFvec.push_back(m_Element.at(ch)->PWF_histo);
hDiffvec.push_back(m_Element.at(ch)->FirstDerivative);
hMarksFront->push_back( new std::vector< TObject* > );
hMarksBack->push_back( new std::vector< TObject* > );
}
}
......@@ -111,7 +127,6 @@ void Detector::FillHistograms(){
m_Element.at(ch)->FirstDerivative->Reset();
m_Element.at(ch)->vDiff.clear();
m_Element.at(ch)->FirstDerivativeRMS = 0;
m_Element.at(ch)->CrossZeroPoints.clear();
m_Element.at(ch)->FitFunc->SetParameters(25000,0.5,2,500);
// Loop over samples in each channel
for( uint samp = 0; samp < m_nSamp; samp++ ){
......@@ -119,3 +134,178 @@ void Detector::FillHistograms(){
} // End loop over samples in each channel
}
}
/** @brief Generic implementation of PlotWF. Draw the waveforms for all Channels for a given event
* WARNING: This is cpu intensive and intended for low level inspection of a few events at a time
*
* @param - Event number for output file naming
*/
void Detector::PlotWF( int eventNum ){
if(eventNum == 0){//If this is the first event create an output directory
gSystem->Exec( Form("mkdir -p %s/results/plots",std::getenv("JZCaPA") ) );
}
for(int ch = 0; ch < m_Element.size(); ch++){
//Delete and clear old markers
for(int i = 0; i < hMarksFront->at(ch)->size(); i++){
if( !hMarksFront->at(ch)->at(i)->InheritsFrom("TF1") ){
delete hMarksFront->at(ch)->at(i);
}
}
hMarksFront->at(ch)->clear();
//Draw the pedestal no matter what
hMarksFront->at(ch)->push_back( new TLine( 0, m_Element[ch]->PedMean, m_nSamp, m_Element[ch]->PedMean) );
((TLine*)hMarksFront->at(ch)->back())->SetLineColor(kGreen);
//Draw hit window lines, peak marker, and set the range of the fit function to the full hit window
if(m_Element[ch]->is_on){
m_Element[ch]->WF_histo->ResetStats();
if(m_Element[ch]->was_hit){
for(int hit = 0; hit < m_Element[ch]->nHits; hit++){
//Draw vertical red lines around the hit window
hMarksFront->at(ch)->push_back( new TLine( m_Element[ch]->hit_window.first[hit], m_Element[ch]->WF_histo->GetMinimum(), m_Element[ch]->hit_window.first[hit], m_Element[ch]->WF_histo->GetMaximum()) );
((TLine*)hMarksFront->at(ch)->back())->SetLineColor(kRed);
hMarksFront->at(ch)->push_back( new TLine( m_Element[ch]->hit_window.second[hit], m_Element[ch]->WF_histo->GetMinimum(), m_Element[ch]->hit_window.second[hit], m_Element[ch]->WF_histo->GetMaximum()) );
((TLine*)hMarksFront->at(ch)->back())->SetLineColor(kRed);
hMarksFront->at(ch)->push_back( new TMarker( m_Element[ch]->Peak_center[hit], m_Element[ch]->Peak_max[hit] + m_Element[ch]->PedMean, 5) );
((TMarker*)hMarksFront->at(ch)->back())->SetMarkerColor(kMagenta+2);
((TMarker*)hMarksFront->at(ch)->back())->SetMarkerSize(3);
}//end hit loop
}//end if channel was hit
}//end if channel is on
}//end channels loop
if(m_viz == NULL){
m_viz = new Visualizer( "ATLAS" );
m_viz->SetTestBeamLabel(m_Alignment2021->runNumber, m_Alignment2021 );
}
int pads = ceil(sqrt(m_Element.size()));
m_viz->SimplePadsPlot( hWFvec, pads, pads, "Time [bins]", "Amplitude [ADC]", //Background histos, foreground histos, nCol, nRow, Detector label
Form("%swf%d.png", m_name.c_str(), eventNum), m_name.c_str(), "", // Output file name, chule name, plot type
Form("%s/results/plots/", std::getenv("JZCaPA") ), // Output directry
true, hMarksFront ); // Autoscale, Plot markers
}
void Detector::PlotWFandDerivative( int eventNum, double _yMin1, double _yMax1, double _yMin2, double _yMax2 ){
if(eventNum == 0){//If this is the first event create an output directory
gSystem->Exec( Form("mkdir -p %s/results/plots",std::getenv("JZCaPA") ) );
}
for(int ch = 0; ch < m_Element.size(); ch++){
//Delete and clear old markers
for(int i = 0; i < hMarksBack->at(ch)->size(); i++){
delete hMarksBack->at(ch)->at(i);
}
hMarksBack->at(ch)->clear();
for(int i = 0; i < hMarksFront->at(ch)->size(); i++){
if( !hMarksFront->at(ch)->at(i)->InheritsFrom("TF1") ){ // Don't delete the fit functions
delete hMarksFront->at(ch)->at(i);
}
}
hMarksFront->at(ch)->clear();
//Draw the pedestal no matter what
hMarksBack->at(ch)->push_back( new TLine( 0, m_Element[ch]->PedMean, m_nSamp, m_Element[ch]->PedMean) );
((TLine*)hMarksBack->at(ch)->back())->SetLineColor(kGreen);
//Plot the threshold lines on the derivative no matter what
hMarksFront->at(ch)->push_back( new TLine( 0, m_Element[ch]->FirstDerivativeRMS*m_Element[ch]->Threshold, m_nSamp, m_Element[ch]->FirstDerivativeRMS*m_Element[ch]->Threshold) );
((TLine*)hMarksFront->at(ch)->back())->SetLineColor(kMagenta);
hMarksFront->at(ch)->push_back( new TLine( 0, -1*m_Element[ch]->FirstDerivativeRMS*m_Element[ch]->Threshold, m_nSamp, -1*m_Element[ch]->FirstDerivativeRMS*m_Element[ch]->Threshold) );
((TLine*)hMarksFront->at(ch)->back())->SetLineColor(kMagenta);
//Draw hit window lines, peak marker, and set the range of the fit function to the full hit window
if(m_Element[ch]->is_on){
m_Element[ch]->WF_histo->GetYaxis()->SetRange(0,0);
m_Element[ch]->WF_histo->ResetStats();
if(m_Element[ch]->was_hit){
for(int hit = 0; hit < m_Element[ch]->nHits; hit++){
//Draw vertical red lines around the hit window
hMarksBack->at(ch)->push_back( new TLine( m_Element[ch]->hit_window.first[hit], -2950, m_Element[ch]->hit_window.first[hit], -2600) );
((TLine*)hMarksBack->at(ch)->back())->SetLineColor(kRed);
hMarksBack->at(ch)->push_back( new TLine( m_Element[ch]->hit_window.second[hit], -2950, m_Element[ch]->hit_window.second[hit], -2600) );
((TLine*)hMarksBack->at(ch)->back())->SetLineColor(kRed);
//Draw an X on the peak
hMarksBack->at(ch)->push_back( new TMarker( m_Element[ch]->Peak_center[hit], m_Element[ch]->Peak_max[hit] + m_Element[ch]->PedMean, 5) );
((TMarker*)hMarksBack->at(ch)->back())->SetMarkerColor(kMagenta+2);
((TMarker*)hMarksBack->at(ch)->back())->SetMarkerSize(3);
}//end hit loop
}//end if channel was hit
}//end if channel is on
}//end channel loop
if(m_viz == NULL){
m_viz = new Visualizer( "ATLAS" );
m_viz->SetTestBeamLabel(m_Alignment2021->runNumber, m_Alignment2021 );
}
int pads = ceil(sqrt(m_Element.size()));
m_viz->ManyPadsPlot( hWFvec, hDiffvec, pads, pads, m_name.c_str(), //Background histos, foreground histos, nCol, nRow, Detector label
Form("%s/results/plots/%swfDiff_ev%d", std::getenv("JZCaPA"), m_name.c_str(), eventNum ), // Output file name
"Time [bins]", "Amplitude [ADC]", "overlay", //X-axis label, Y-axis label, plot type
hMarksBack, hMarksFront, //Background plot markers, Foreground plot markers
_yMin1, _yMax1, _yMin2, _yMax2); //Y axis ranges
}
void Detector::PlotWFandPWF( int eventNum, double _yMin1, double _yMax1, double _yMin2, double _yMax2 ){
if(eventNum == 0){//If this is the first event create an output directory
gSystem->Exec( Form("mkdir -p %s/results/plots",std::getenv("JZCaPA") ) );
}
for(int ch = 0; ch < m_Element.size(); ch++){
//Delete and clear old markers
for(int i = 0; i < hMarksBack->at(ch)->size(); i++){
delete hMarksBack->at(ch)->at(i);
}
hMarksBack->at(ch)->clear();
for(int i = 0; i < hMarksFront->at(ch)->size(); i++){
if( !hMarksFront->at(ch)->at(i)->InheritsFrom("TF1") ){
delete hMarksFront->at(ch)->at(i);
}
}
hMarksFront->at(ch)->clear();
//Draw the pedestal no matter what
hMarksBack->at(ch)->push_back( new TLine( 0, m_Element[ch]->PedMean, m_nSamp, m_Element[ch]->PedMean) );
((TLine*)hMarksBack->at(ch)->back())->SetLineColor(kGreen);
//Draw hit window lines, peak marker, and set the range of the fit function to the full hit window
if(m_Element[ch]->is_on){
m_Element[ch]->WF_histo->ResetStats();
if(m_Element[ch]->was_hit){
for(int hit = 0; hit < m_Element[ch]->nHits; hit++){
//Draw vertical red lines around the hit window
hMarksBack->at(ch)->push_back( new TLine( m_Element[ch]->hit_window.first[hit], m_Element[ch]->WF_histo->GetMinimum(), m_Element[ch]->hit_window.first[hit], m_Element[ch]->WF_histo->GetMaximum()) );
((TLine*)hMarksBack->at(ch)->back())->SetLineColor(kRed);
hMarksBack->at(ch)->push_back( new TLine( m_Element[ch]->hit_window.second[hit],m_Element[ch]->WF_histo->GetMinimum(), m_Element[ch]->hit_window.second[hit], m_Element[ch]->WF_histo->GetMaximum()) );
((TLine*)hMarksBack->at(ch)->back())->SetLineColor(kRed);
hMarksBack->at(ch)->push_back( new TMarker( m_Element[ch]->Peak_center[hit], m_Element[ch]->Peak_max[hit] + m_Element[ch]->PedMean, 5) );
((TMarker*)hMarksBack->at(ch)->back())->SetMarkerColor(kMagenta+2);
((TMarker*)hMarksBack->at(ch)->back())->SetMarkerSize(3);
// hMarksFront->at(ch)->push_back( m_Element[ch]->FitFunc );
// ((TF1*)hMarksFront->at(ch)->back())->SetRange(m_Element[ch]->hit_window.first, m_Element[ch]->hit_window.second);
// ((TF1*)hMarksFront->at(ch)->back())->SetLineColor(kCyan);
}//end hit loop
}//end if channel was hit
}//end if channel is on
}//end channel loop
if(m_viz == NULL){
m_viz = new Visualizer( "ATLAS" );
m_viz->SetTestBeamLabel(m_Alignment2021->runNumber, m_Alignment2021 );
}
int pads = ceil(sqrt(m_Element.size()));
m_viz->ManyPadsPlot( hWFvec, hPWFvec, pads, pads, m_name.c_str(), //Background histos, foreground histos, nCol, nRow, Detector label
Form("%s/results/plots/%swfPwfPlot_ev%d",std::getenv("JZCaPA"), m_name.c_str(), eventNum ), // Output file name
"Time [bins]", "Amplitude [ADC]", "overlay", //X-axis label, Y-axis label, plot type
hMarksBack, hMarksFront, //Background plot markers, Foreground plot markers
_yMin1, _yMax1, _yMin2, _yMax2); //Y axis ranges
}
......@@ -113,17 +113,17 @@ void EMAnalysis::SetBranches( TTree* _tree ){
for(int row = 0; row < m_numRows; row++){
for(int col = 0; col < m_numCols; col++){
m_AnalysisTree->Branch( Form("em%d_%d_Charge", row, col), &em[row][col]->Charge, Form("em%d_%d_Charge/D", row, col) );
m_AnalysisTree->Branch( Form("em%d_%d_Peak_max", row, col), &em[row][col]->Peak_max, Form("em%d_%d_Peak_max/D", row, col) );
m_AnalysisTree->Branch( Form("em%d_%d_Diff_max", row, col), &em[row][col]->Diff_max, Form("em%d_%d_Diff_max/D", row, col) );
m_AnalysisTree->Branch( Form("em%d_%d_Peak_center", row, col), &em[row][col]->Peak_center, Form("em%d_%d_Peak_center/I", row, col) );
m_AnalysisTree->Branch( Form("em%d_%d_Diff_Peak_center", row, col), &em[row][col]->Diff_Peak_center, Form("em%d_%d_Diff_Peak_center/I", row, col) );
m_AnalysisTree->Branch( Form("em%d_%d_Charge", row, col), "std::vector<double>", &em[row][col]->Charge );
m_AnalysisTree->Branch( Form("em%d_%d_Peak_max", row, col), "std::vector<double>", &em[row][col]->Peak_max );
m_AnalysisTree->Branch( Form("em%d_%d_Diff_max", row, col), "std::vector<double>", &em[row][col]->Diff_max );
m_AnalysisTree->Branch( Form("em%d_%d_Peak_center", row, col), "std::vector<int>", &em[row][col]->Peak_center );
m_AnalysisTree->Branch( Form("em%d_%d_Diff_Peak_center", row, col), "std::vector<int>", &em[row][col]->Diff_Peak_center );
}
}
m_AnalysisTree->Branch("em_Charge_sum", &ChargeSum, "ChargeSum/D" );
m_AnalysisTree->Branch("em_Peak_sum", &PeakSum, "PeakSum/D" );
m_AnalysisTree->Branch("em_Diff_Peak_sum", &DiffPeakSum, "DiffPeakSum/D" );
m_AnalysisTree->Branch("em_Charge_sum", &ChargeSum, "ChargeSum/D" );
m_AnalysisTree->Branch("em_Peak_sum", &PeakSum, "PeakSum/D" );
m_AnalysisTree->Branch("em_Diff_Peak_sum",&DiffPeakSum,"DiffPeakSum/D" );
}
......@@ -140,16 +140,18 @@ void EMAnalysis::AnalyzeEvent( ){
for(int row = 0; row < m_numRows; row++){
for(int col = 0; col < m_numCols; col++){
if( em[row][col]->is_on ){
ChargeSum += em[row][col]->Charge;
PeakSum += em[row][col]->Peak_max;
DiffPeakSum += em[row][col]->Diff_max;
hChargeArr[row][col]->Fill(em[row][col]->Charge);
hPeakArr[row][col]->Fill(em[row][col]->Peak_max);
hDPeakArr[row][col]->Fill(em[row][col]->Diff_max);
}
}
}
for(int hit = 0; hit < em[row][col]->nHits; hit++){
ChargeSum += em[row][col]->Charge[hit];
PeakSum += em[row][col]->Peak_max[hit];
DiffPeakSum += em[row][col]->Diff_max[hit];
hChargeArr[row][col]->Fill(em[row][col]->Charge[hit]);
hPeakArr[row][col]->Fill(em[row][col]->Peak_max[hit]);
hDPeakArr[row][col]->Fill(em[row][col]->Diff_max[hit]);
}//end hits loop
}//end if channel is on
}//end column loop
}//end row loop
//This loop finds the center of mass per event (should probably be it's own function)
......@@ -159,9 +161,11 @@ void EMAnalysis::AnalyzeEvent( ){
for(int i = 0; i < m_numRows; i++){
for(int j = 0; j < m_numCols; j++){
if( em[i][j]->is_on && em[i][j]->was_hit ){
totalCharge += em[i][j]->Charge;
weightedRow += em[i][j]->Charge * zPos[i];
weightedCol += em[i][j]->Charge * xPos[j];
for(int hit = 0; hit < em[i][j]->nHits; hit++){
totalCharge += em[i][j]->Charge[hit];
weightedRow += em[i][j]->Charge[hit] * zPos[i];
weightedCol += em[i][j]->Charge[hit] * xPos[j];
}
}
}
}
......@@ -170,6 +174,11 @@ void EMAnalysis::AnalyzeEvent( ){
hPeakSum->Fill(PeakSum);
hDiffPeakSum->Fill(DiffPeakSum);
// m_EM->PlotWF(m_counter);
// m_EM->PlotWFandDerivative(m_counter, -3300, -2550, -400, 900);
// m_EM->PlotWFandPWF(m_counter, -3300, -2550, -10, 700);
m_counter++;
}
/** @brief Finalize method for EMAnalysis
......
......@@ -117,11 +117,11 @@ void RPDAnalysis::SetBranches( TTree* _tree ){
for(int row = 0; row < 4; row++){
for(int col = 0; col < 4; col++){
m_AnalysisTree->Branch( Form("rpd%d_%d_Charge", row, col), &rpd[row][col]->Charge, Form("rpd%d_%d_Charge/D", row, col) );
m_AnalysisTree->Branch( Form("rpd%d_%d_Peak_max", row, col), &rpd[row][col]->Peak_max, Form("rpd%d_%d_Peak_max/D", row, col) );
m_AnalysisTree->Branch( Form("rpd%d_%d_Diff_max", row, col), &rpd[row][col]->Diff_max, Form("rpd%d_%d_Diff_max/D", row, col) );
m_AnalysisTree->Branch( Form("rpd%d_%d_Peak_center", row, col), &rpd[row][col]->Peak_center, Form("rpd%d_%d_Peak_center/I", row, col) );
m_AnalysisTree->Branch( Form("rpd%d_%d_Diff_Peak_center", row, col), &rpd[row][col]->Diff_Peak_center, Form("rpd%d_%d_Diff_Peak_center/I", row, col) );
m_AnalysisTree->Branch( Form("rpd%d_%d_Charge", row, col), "std::vector<double>", &rpd[row][col]->Charge );
m_AnalysisTree->Branch( Form("rpd%d_%d_Peak_max", row, col), "std::vector<double>", &rpd[row][col]->Peak_max );
m_AnalysisTree->Branch( Form("rpd%d_%d_Diff_max", row, col), "std::vector<double>", &rpd[row][col]->Diff_max );
m_AnalysisTree->Branch( Form("rpd%d_%d_Peak_center", row, col), "std::vector<int>", &rpd[row][col]->Peak_center );
m_AnalysisTree->Branch( Form("rpd%d_%d_Diff_Peak_center", row, col), "std::vector<int>", &rpd[row][col]->Diff_Peak_center );
}
}
......@@ -146,16 +146,18 @@ void RPDAnalysis::AnalyzeEvent( ){
for(int row = 0; row < 4; row++){
for(int col = 0; col < 4; col++){
if( rpd[row][col]->is_on ){
ChargeSum += rpd[row][col]->Charge;
PeakSum += rpd[row][col]->Peak_max;