--- /dev/null
+#include <stdio.h>
+#include <math.h>
+#include <string>
+#include <sstream>
+#include <fstream>
+#include <vector>
+
+#include "base/EventProc.h"
+#include "base/Event.h"
+#include "hadaq/TdcSubEvent.h"
+#include "hadaq/HldProcessor.h"
+#include <typeinfo>
+
+#include <iostream>
+
+#include "TFile.h"
+#include "TH2.h"
+#include "TH1.h"
+
+using namespace std;
+
+class SecondProc : public base::EventProc
+{
+ protected:
+ /*Analysis flags*/
+ bool makeTdiffCorr = false;
+ bool makeTdiffCorrFunction = false;
+ bool makeTWCorrParameters = false;
+ bool makeTWCorrFunction = false;
+ bool makeRescale = false;
+ bool makeClusterRemovalCorrCh = false;
+ bool makeClusterRemovalBothCh = false;
+ bool makeCutsCorrChannel = false;
+ bool makeCutsBothChannels = false;
+ bool makeCrosstalkAnalysis = false;
+ bool makeClusterAnalysis = false;
+
+
+ double CorrCutMax=10000;
+ double CorrCutMin=0;
+ double RefCutMax=10000;
+ double RefCutMin=0;
+
+ double TDiffCut = 10000000.;
+
+ static const int numOfTDCs = 5; // number of used TDCs +1
+ static const int numOfChannels = 17; // number of channels +1
+ static const int numOfDetectors = 3; // number of detectors + 1
+
+ int eventnumber=0;
+
+ // Correlation scheme for setup: [TDC][channel][{channel_number, correlated_channel_number_in_other_TDC, detector_number, position_in_detector}]
+ int correlationScheme[numOfTDCs][numOfChannels][4] = {
+
+ {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}},
+
+ {{0,0,0,0}, {1,1,1,2}, {2,2,1,4}, {3,3,1,6}, {4,4,1,8}, {5,5,1,10}, {6,6,1,12}, {7,5,1,14}, {8,8,1,16}, {9,9,1,15}, {10,10,1,13}, {11,12,1,11}, {12,12,1,9}, {13,13,1,7}, {14,14,1,5}, {15,15,1,3}, {16,16,1,1}},
+
+ {{0,0,0,0}, {1,2,1,2}, {2,1,1,4}, {3,3,1,6}, {4,4,1,8}, {5,5,1,10}, {6,6,1,12}, {7,5,1,14}, {8,8,1,16}, {9,9,1,15}, {10,10,1,13}, {11,12,1,11}, {12,12,1,9}, {13,13,1,7}, {14,14,1,5}, {15,15,1,3}, {16,16,1,1}},
+
+ {{0,0,0,0}, {1,1,1,2}, {2,2,1,4}, {3,3,1,6}, {4,4,1,8}, {5,5,1,10}, {6,6,1,12}, {7,5,1,14}, {8,8,1,16}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}},
+
+ {{0,0,0,0}, {1,1,2,2}, {2,2,2,4}, {3,3,2,6}, {4,4,2,8}, {5,5,2,10}, {6,6,2,12}, {7,5,2,14}, {8,8,2,16}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}}
+
+ };
+
+
+
+// int correlationScheme[numOfTDCs][numOfChannels][4] = { // Correlation scheme which correlates only with channel 1 of the other detector
+//
+// {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}},
+//
+// {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {9,9,1,15}, {10,10,1,13}, {11,12,1,11}, {12,12,1,9}, {13,13,1,7}, {14,14,1,5}, {15,15,1,3}, {16,16,1,1}},
+//
+// {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {9,9,2,15}, {10,10,2,13}, {11,12,2,11}, {12,12,2,9}, {13,13,2,7}, {14,14,2,5}, {15,15,2,3}, {16,16,2,1}},
+//
+// {{0,0,0,0}, {1,1,1,2}, {2,1,1,4}, {3,1,1,6}, {4,1,1,8}, {5,1,1,10}, {6,1,1,12}, {7,1,1,14}, {8,1,1,16}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}},
+//
+// {{0,0,0,0}, {1,1,2,2}, {2,1,2,4}, {3,1,2,6}, {4,1,2,8}, {5,1,2,10}, {6,1,2,12}, {7,1,2,14}, {8,1,2,16}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}}
+//
+// };
+
+
+ double TWalkCorrFunction [numOfTDCs][numOfChannels][4]={
+
+// {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}},
+//
+// {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}},
+//
+// {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}},
+//
+// {{0,0,0,0}, {-183.685, 2.17905, 0.00191316, 0.160909}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}},
+//
+// {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}},
+
+
+
+ { {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}},
+
+ { {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}},
+
+ { {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}},
+
+ { {0,0,0,0}, {-202.368,2.2181,0.000303494,0.0765854}, {-176.324,2.10958,-0.00205473,0.178147}, {-237.962,2.34708,0.00103295,0.0445063}, {-291.213,2.42026,0.000945137,0.0717607}, {-196.381,2.23578,8.7334e-05,0.0848016}, {-171.153,2.14604,-8.42002e-05,0.0879119}, {-375.62,2.25024,-0.42893,6.41755}, {-163.742,1.82885,-0.00667092,0.472399}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}},
+
+ { {0,0,0,0}, {-214.793,2.35896,0.00166868,0.121898}, {-224.044,2.33503,0.00109374,0.158345}, {-233.527,2.40481,0.00249939,0.030645}, {-263.802,2.42822,0.00152629,0.0869203}, {-268.726,2.44789,0.0015527,0.116454}, {-112.684,2.0367,-0.00324999,0.490074}, {-263.643,2.51534,0.00296977,0.0282881}, {-305.755,2.25468,0.000236872,0.202798}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0}}
+
+
+
+
+ };
+
+
+
+
+// double cutValuesLower [numOfTDCs] [numOfChannels]= {
+//
+// {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,},
+//
+// {0, 0, 0, 0, 0, 0, 0, 0, 0, 4.25668, 5.28533, 4.64967, 5.47919, 4.89474, 4.1154, 5.33308, 6.47866, 0, 0, 0, 0, 0, 0, 0, 0, 5.38717, 5.64863, 0, 5.17477, 5.95315, 5.34243, 5.13881, 6.10115},
+//
+// {0, 0, 0, 0, 0, 0, 0, 0, 0, 6.59853, 6.93174, 6.53145, 6.88821, 6.51382, 6.3149, 0, 6.68963, 0, 0, 0, 0, 0, 0, 0, 0, 16.5065, 6.78918, 0, 15.8425, 0, 16.0851, 0, 7.1868}
+// } ;
+
+ double cutValuesLower [numOfTDCs] [numOfChannels]= { // test to check if ToT cuts after cluster removal are necessary, or really improve resolution
+
+ {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+
+ {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+
+ {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+
+ {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+
+ {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
+ } ;
+
+// double cutValuesUpper [numOfTDCs] [numOfChannels]= {
+//
+// {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,},
+//
+// {0, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 9.70702, 10.2622, 10.1603, 10.3137, 9.33094, 9.10002, 9.26649, 10.2667, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 10.2682, 11.3712, 0, 8.70153, 10.0298, 9.49266, 9.26477, 10.3372,},
+//
+// {0,1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 11.8185, 11.7319, 12.0786, 11.3878, 10.693, 10.2768, 0, 10.856, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 24.4851, 10.942, 0, 25.1402, 0, 25.0127, 0, 10.259 }
+// } ;
+
+ double cutValuesUpper [numOfTDCs] [numOfChannels]= { // cut values set by hand ca. to bins, where TWalk correction is not done due to missing statistics to check if ToT cuts after cluster removal are necessary, or really improve resolution
+
+ {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+
+ {0, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000},
+
+ {0, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000},
+
+ {0, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000},
+
+ {0, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000}
+ } ;
+
+
+
+
+ int detectorScheme[numOfDetectors][17][2];
+
+ /* Variables */
+ std::string fTdcId1; //!< tdc id where channels will be selected
+ std::string fTdcId2; //!< tdc id where channels will be selected
+ std::string fTdcId3; //!< tdc id where channels will be selected
+ std::string fTdcId4; //!< tdc id where channels will be selected
+
+ double fHits[numOfTDCs][numOfChannels][3]; // 4 TDCs, 17 channels, 2 edges + ToT
+ double fHitsCleared[numOfTDCs][numOfChannels][3]; // 4 TDCs, 17 channels, 2 edges + ToT
+ double fHitsDetector[numOfDetectors][18][3]; //3 detectors, 17 channels, 2 edges + ToT rescaled
+
+ std::vector<std::vector<int>> detectorClusters[numOfDetectors];
+
+ int fSubEvents [numOfTDCs][numOfChannels][2];
+ std::vector<double> fSubTdiff[numOfTDCs][numOfChannels];
+
+ std::vector<double> Timecorr;
+ std::vector<double> rescale;
+
+
+ /* Histograms */
+ base::H1handle hNumHits;
+ base::H1handle hNumEntries;
+
+ base::H1handle clusterSizeDetector[numOfDetectors];
+ base::H1handle multiplicityDetector[numOfDetectors];
+
+ base::H2handle ToTAllChannels;
+
+ base:: H2handle hChMap[numOfTDCs];
+
+ //Number of events in subevent per channel for each TDC
+ base::H2handle hChEvents[numOfTDCs];
+
+ //Subevents Tdiff per event per chanel for each TDC
+ base::H2handle hChTdiff1[numOfTDCs];
+ base::H2handle hChTdiff2[numOfTDCs];
+
+ // Create ToTvsTDiff Histograms
+ base::H2handle ToTvsTDiff[numOfTDCs][numOfChannels];
+ base::H2handle TDiffAllChannels[numOfTDCs][numOfChannels];
+
+
+ base::H2handle NeighbourCrosstalk[numOfTDCs][numOfChannels];
+ base::H2handle LocalMaxNeighbourCrosstalk[numOfTDCs][numOfChannels];
+
+ base::H2handle Right_TDiff_Crosstalk[numOfTDCs][numOfChannels];
+ base::H2handle Left_TDiff_Crosstalk[numOfTDCs][numOfChannels];
+
+ base::H2handle LocalMax_Right_TDiff_Crosstalk[numOfTDCs][numOfChannels];
+ base::H2handle LocalMax_Left_TDiff_Crosstalk[numOfTDCs][numOfChannels];
+
+
+ base::H2handle Right_ToT_Crosstalk[numOfTDCs][numOfChannels];
+ base::H2handle Left_ToT_Crosstalk[numOfTDCs][numOfChannels];
+
+ base::H2handle LocalMax_Right_ToT_Crosstalk[numOfTDCs][numOfChannels];
+ base::H2handle LocalMax_Left_ToT_Crosstalk[numOfTDCs][numOfChannels];
+
+ base::H1handle sparkMultiplicity[numOfTDCs];
+
+
+ //Timediff multiplicity
+ base::H2handle TDTF[numOfTDCs];
+
+ base::H2handle hEdges[numOfTDCs][numOfChannels];
+
+ //TW histograms
+ TH2F* TW_histo[5];
+
+ public:
+ SecondProc(const char* procname, const char* _tdcid1, const char* _tdcid2, const char* _tdcid3, const char* _tdcid4) :
+ base::EventProc(procname),
+ fTdcId1(_tdcid1),
+ fTdcId2(_tdcid2),
+ fTdcId3(_tdcid3),
+ fTdcId4(_tdcid4),
+ ToTAllChannels(0)
+ {
+ //Time Walk correction source file reading
+ if(makeTWCorrParameters) {
+ TFile *inFile = new TFile("2M_TW.root");
+ if(!inFile) cerr<<"TW file not found \n"<<endl;
+
+// Select Histogram - pointer to the histogram
+ for(int i = 1; i < numOfTDCs; i++)
+ TW_histo[i] = (TH2F*)inFile->Get(Form("histo%i", i));
+
+ }
+
+
+
+
+ if(makeTdiffCorr || makeTdiffCorrFunction)
+ Timecorr = {
+
+ // TDiff calibration before everything else
+// 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+// 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+//
+// 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.379719,
+// -0.807165, 0, 4.29306, 1.92666, 0.98097, 1.08864, 1.77656, 0, 0, 0,
+//
+// 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+// 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+//
+// 0, 1.07526, 0.138509, 0.779046, 1.56691, 1.46053, 1.59004, 1.00212, 0.464563, 0,
+// 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+//
+// 0
+ // after TWalk correction
+// 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.259997,
+// 0.234026, 3.03417, 5.12116, 2.59994, 1.75335, 1.74567, 2.51493, 0, 0, 0,
+//
+// 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.389521,
+// -0.810575, 0, 4.29027, 1.92367, 0.984565, 1.08809, 1.7731, 0, 0, 0,
+//
+// 0, 0.973676, -1.5956e-05, 0.4992, 1.17251, 1.31455, 1.55055, -1.33833, 0.288798, 0,
+// 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+//
+// 0, 1.0702, 0.140716, 0.7769, 1.56193, 1.45633, 1.59049, 0.995151, 0.458112, 0,
+// 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+//
+// 0
+ // after TWalk correction and double bin (5 ps bin Width)
+// 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.259981,
+// 0.234088, 3.03421, 5.12126, 2.6, 1.75323, 1.74567, 2.51498, 0, 0, 0,
+//
+// 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.389559,
+// -0.810571, 0, 4.29025, 1.92365, 0.984593, 1.08807, 1.7953, 0, 0, 0,
+//
+// 0, 0.973764, -7.69487e-06, 0.499233, 1.17253, 1.3146, 1.55058, -1.3353, 0.288827, 0,
+// 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+//
+// 0, 1.07014, 0.140753, 0.776919, 1.56193, 1.45629, 1.59047, 0.995192, 0.458113, 0,
+// 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+//
+// 0
+
+
+ // TDiff correction parameter after cuts, TWalk correction and cluster removal, 10 ps bin width
+// 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+// 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+//
+// 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.399899,
+// -0.812385, 0, 4.28664, 1.92076, 0.986521, 1.09002, 1.77038, 0, 0, 0,
+//
+// 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+// 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+//
+// 0, 1.06717, 0.140334, 0.776773, 1.55831, 1.45336, 1.58836, 0.992174, 0.456232, 0,
+// 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+// 0
+
+
+ // TDiff correction parameter after cuts, TWalk correction and cluster removal, 5 ps bin width
+
+
+ /*
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.256519,
+ 0.229751, 3.03087, 5.122, 2.6181, 1.75288, 1.75239, 2.51801, 0, 0, 0,
+
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.399977,
+ -0.812364, 0, 4.28658, 1.92077, 0.986573, 1.09, 1.7705, 0, 0, 0,
+
+ 0, 0.972426, -1.15853e-05, 0.497899, 1.17262, 1.31565, 1.54837, 0, 0.289097, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+
+ 0, 1.06703, 0.140375, 0.776698, 1.55825, 1.45338, 1.58829, 0.992229, 0.456329, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 0
+ */
+
+
+
+ // Mean TDiff parameters. Acquired by first making normal TDiff calibration, then finding the mean correction parameters for the other TDC and then running analysis to get mean correction
+ //parameters of the first TDC
+
+// 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.256519,
+// 0.229751, 3.03087, 5.122, 2.6181, 1.75288, 1.75239, 2.51801, 0, 0, 0,
+//
+// 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.399977,
+// -0.812364, 0, 4.28658, 1.92077, 0.986573, 1.09, 1.7705, 0, 0, 0,
+//
+//
+// 0, 0.978305 , 0, 0.50338, 1.17631, 1.31946, 1.55397, 0, 0.293874, 0,
+// 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+//
+// 0, 1.07614, 0.14949, 0.789326, 1.56768, 1.46074, 1.59607, 1.00274, 0.467685, 0,
+// 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+// 0
+
+
+ // New TDiff parameters generated after TWalk with function in one step
+/*
+
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.241436,
+ 0.216216, 3.00324, 5.11242, 2.53487, 1.71291, 1.69343, 2.44475, 0, 0, 0,
+
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.432846,
+ -0.865883, 0, 4.22372, 1.89486, 0.95267, 1.1048, 1.78417, 0, 0, 0,
+
+ 0, 0.956509, -0.00782424 , 0.483118, 1.15808, 1.30818, 1.58146, 0, 0.286464, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+
+ 0, 1.06409, 0.134968, 0.762091, 1.55373, 1.43657, 1.55457, 0.96967, 0.467006, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 0, */
+
+
+ // TDiff Parameters after TWalk correction with a function and Cluster removal and ToT cuts generated in Two steps/two files
+
+// 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.249248,
+// 0.224022, 3.01106, 5.1202, 2.54274, 1.7208, 1.70127, 2.45253, 0, 0, 0,
+//
+// 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.425687,
+// -0.861085, 0, 4.22932, 1.89926, 0.959421, 1.11553, 1.79165, 0, 0, 0,
+//
+// 0, 0.964389, 3.0101e-05, 0.490929, 1.16599, 1.31605, 1.58926, 0, 0.294278, 0,
+// 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+//
+// 0, 1.07204, 0.142791, 0.771455, 1.55964, 1.44587, 1.56574, 0.978453, 0.469983, 0,
+// 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+// 0
+
+ // TDiff Parameters after TWalk correction with a function and Cluster removal and ToT cuts generated in one step/one file
+
+
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.106466,
+ 0.0812532, 2.86827, 4.97746, 2.39991, 1.57795, 1.55847, 2.30977, 0, 0, 0,
+
+
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.425687,
+ -0.861085, 0, 4.22932, 1.89926, 0.959421, 1.11553, 1.79165, 0, 0, 0,
+
+ 0, 0.821543, -0.142791, 0.348155, 1.02312, 1.1732, 1.4465, 0, 0.151474, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+
+ 0, 1.07204, 0.142791, 0.771455, 1.55964, 1.44587, 1.56574, 0.978453, 0.469983, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 0
+
+
+
+
+ };else Timecorr = std::vector<double> (81, 0);
+
+ if(/*makeClusterRemovalCorrCh || makeClusterRemovalBothCh ||*/ makeRescale)
+ rescale = {
+ 1,
+
+ 1, 1, 1, 1, 1, 1, 1, 1, 4.02951, 3.94215, 4.44876, 4.36756, 4.26959, 4.07083, 4.39062, 3.93766, 1, 1, 1, 1,
+ 1, 1, 1, 1, 1, 1, 1, 1, 4.32465, 4.53967, 1, 4.78249, 4.63644, 4.60738, 4.75951, 4.45675, 1, 1, 1, 1,
+
+ // normalize nino to 35
+ 1.09407, 1.12791, 1.09457, 1.20249, 1.10512, 1.09579, 1 , 0.998456, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
+ 1.55471, 1.58341, 1.2849, 1.49527, 1.59509, 2.31983, 1.4235, 1.39305, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
+
+
+
+ };
+ else
+ rescale = std::vector<double> (81, 1);
+
+ hNumHits = MakeH1("NumHits","Num hits", 100, 0, 100, "number");
+ hNumEntries = MakeH1("NumEntries", "Num Entries", 60, 0, 60, "channel; entries");
+
+ // Time over Threshold for all channels histogram
+ ToTAllChannels = MakeH2("ToTAllChannels","ToTAllChannels" , 100, 0, 100, 1000, -10, 200, "TDC_ch;ToT[ns]");
+
+ char histoName[100];
+ for (int i = 1; i < numOfTDCs; i++) {
+ sprintf(histoName,"Spark_Multiplicity_TDC%i", i);
+ sparkMultiplicity[i] = MakeH1(histoName, histoName, 35, 0, 35, "Multiplicity; Entries");
+
+ for (int j=1; j<numOfChannels; j++){
+ if(correlationScheme[i][j][0]){
+
+ if (i<3){
+ sprintf(histoName,"ToT%i_vs_TDiff_ch%i_%i", i,correlationScheme[i][j][0] ,correlationScheme[i][j][1]);
+ ToTvsTDiff[i][j] = MakeH2(histoName, histoName,600,0,60,500,-10,10,"ToT [ns]; TDiff [ns]");
+ }
+
+ else if (i>2){
+ sprintf(histoName,"ToT%i_vs_TDiff_ch%i_%i", i,correlationScheme[i][j][0] ,correlationScheme[i][j][1]);
+ ToTvsTDiff[i][j] = MakeH2(histoName, histoName,180,0,60,500,-10,10,"ToT [ns]; TDiff [ns]");
+ }
+
+ if(j < 9 || (j > 8 && correlationScheme[i][j][0] != correlationScheme[i][j - 8][0] )){
+ sprintf(histoName,"RefTDC%i_ch%i", i,correlationScheme[i][j][0]);
+ TDiffAllChannels[i][j] = MakeH2(histoName, histoName, 80,0,80,4000,-10.0025,9.9975,"Second_TDC_channels; TDiff [ns]");
+ }
+ sprintf(histoName,"NumberOfEdgesTDC_%i_ch_%i", i,correlationScheme[i][j][0]);
+ hEdges[i][j] = MakeH2(histoName, histoName, 20,0,20,20,0,20,"Leading edges; Trailing edges");
+
+ if (makeCrosstalkAnalysis){
+ if (i>2){
+ sprintf(histoName,"NeighbouringCrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]);
+ NeighbourCrosstalk[i][j] = MakeH2(histoName, histoName, 180,0,60,5 ,0,5,"ToT [ns]; ClusterSize");
+
+ sprintf(histoName,"LocalMaxNeighbouringCrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]);
+ LocalMaxNeighbourCrosstalk[i][j] = MakeH2(histoName, histoName, 180,0,60,5 ,0,5,"ToT [ns]; ClusterSize");
+
+ Left_TDiff_Crosstalk[i][j] = MakeH2(Form("Left_TDiff_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), Form("Left_TDiff_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), 180,0,60,250 ,-5,5,"ToT [ns]; TDiff_Left");
+ Right_TDiff_Crosstalk[i][j] = MakeH2(Form("Right_TDiff_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), Form("Right_TDiff_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), 180,0,60,250 ,-5,5,"ToT [ns]; TDiff_Right");
+
+ LocalMax_Left_TDiff_Crosstalk[i][j] = MakeH2(Form("LocalMaxLeft_TDiff_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), Form("LocalMaxLeft_TDiff_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), 180,0,60,250 ,-5,5,"ToT [ns]; TDiff_Left");
+ LocalMax_Right_TDiff_Crosstalk[i][j] = MakeH2(Form("LocalMaxRight_TDiff_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), Form("LocalMaxRight_TDiff_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), 180,0,60,250 ,-5,5,"ToT [ns]; TDiff_Right");
+
+
+ Left_ToT_Crosstalk[i][j] = MakeH2(Form("Left_ToT_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), Form("Left_ToT_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), 180,0,60,180 ,0,18,"ToT [ns]; ToT_Left");
+ Right_ToT_Crosstalk[i][j] = MakeH2(Form("Right_ToT_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), Form("Right_ToT_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), 180,0,60,180 ,0,15,"ToT [ns]; ToT_Right");
+
+ LocalMax_Left_ToT_Crosstalk[i][j] = MakeH2(Form("LocalMaxLeft_ToT_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), Form("LocalMaxLeft_ToT_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), 90,0,60,90 ,0,60,"ToT [ns]; ToT_Left");
+ LocalMax_Right_ToT_Crosstalk[i][j] = MakeH2(Form("LocalMaxRight_ToT_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), Form("LocalMaxRight_ToT_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), 90,0,60,90 ,0,60,"ToT [ns]; ToT_Right");
+ }
+
+
+ else if (i<3){
+ sprintf(histoName,"NeighbouringCrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]);
+ NeighbourCrosstalk[i][j] = MakeH2(histoName, histoName, 180,0,18,5 ,0,5,"ToT [ns]; ClusterSize");
+
+ sprintf(histoName,"LocalMaxNeighbouringCrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]);
+ LocalMaxNeighbourCrosstalk[i][j] = MakeH2(histoName, histoName, 180,0,18,5 ,0,5,"ToT [ns]; ClusterSize");
+
+ Left_TDiff_Crosstalk[i][j] = MakeH2(Form("Left_TDiff_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), Form("Left_TDiff_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), 180,0,18,250 ,-5,5,"ToT [ns]; TDiff_Left");
+ Right_TDiff_Crosstalk[i][j] = MakeH2(Form("Right_TDiff_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), Form("Right_TDiff_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), 180,0,18,250 ,-5,5,"ToT [ns]; TDiff_Right");
+
+ LocalMax_Left_TDiff_Crosstalk[i][j] = MakeH2(Form("LocalMaxLeft_TDiff_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), Form("LocalMaxLeft_TDiff_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), 180,0,18,250 ,-5,5,"ToT [ns]; TDiff_Left");
+ LocalMax_Right_TDiff_Crosstalk[i][j] = MakeH2(Form("LocalMaxRight_TDiff_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), Form("LocalMaxRight_TDiff_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), 180,0,18,250 ,-5,5,"ToT [ns]; TDiff_Right");
+
+
+ Left_ToT_Crosstalk[i][j] = MakeH2(Form("Left_ToT_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), Form("Left_ToT_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), 180,0,18,180 ,0,60,"ToT [ns]; ToT_Left");
+ Right_ToT_Crosstalk[i][j] = MakeH2(Form("Right_ToT_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), Form("Right_ToT_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), 180,0,18,180 ,0,60,"ToT [ns]; ToT_Right");
+
+ LocalMax_Left_ToT_Crosstalk[i][j] = MakeH2(Form("LocalMaxLeft_ToT_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), Form("LocalMaxLeft_ToT_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), 90,0,60,90 ,0,60,"ToT [ns]; ToT_Left");
+ LocalMax_Right_ToT_Crosstalk[i][j] = MakeH2(Form("LocalMaxRight_ToT_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), Form("LocalMaxRight_ToT_CrosstalkTDC_%i_ch%i", i,correlationScheme[i][j][0]), 90,0,60,90 ,0,60,"ToT [ns]; ToT_Right");
+
+
+
+ }
+ }
+
+ }
+ }
+ if (i<numOfDetectors){
+ sprintf(histoName,"MulitplicityDetector_%i", i);
+ multiplicityDetector[i] = MakeH1(histoName, histoName, 17, 0,17,"Multiplicity");
+ sprintf(histoName,"clusterSizeDetector%i", i);
+ clusterSizeDetector[i] = MakeH1(histoName, histoName, 17, 0,17,"ClusterSize");
+ }
+
+ //Number of events in subevent per channel for each TDC
+ sprintf(histoName,"SignalsPerEventTDC%i", i);
+ hChEvents[i] = MakeH2(histoName, histoName, 16, 1,17, 30, 0, 30,"channel;# of edges");
+
+ // Multiplicity for time difference histograms
+ sprintf(histoName, "TDTF%i", i);
+ TDTF[i] = MakeH2(histoName, histoName,4000,0,1000,17,0,17,"Tdiff[ns];TDC1_ch");
+
+ // Channel correlation mapping histograms
+ sprintf(histoName, "ChMapTDC%i_TDCs", i);
+ char axisLabel[50];
+ sprintf(axisLabel, "ch_TDCs;ch_refTDC%i", i);
+ hChMap[i] = MakeH2(histoName, histoName, 80, 0, 80, 80, 0, 80, axisLabel);
+ }
+
+ CreateDetectorScheme();
+ }
+
+ double CalculateTWCorrection(double ToT) {
+ if(ToT >= 9.45)
+ return 0.0106702 * ToT - 0.112834;
+ else
+ return -0.0313777 * ToT * ToT + 0.527936 * ToT - 2.2427;
+
+ }
+
+ long Unpack(hadaq::TdcSubEventFloat* sub, int shift, double hits[][3], int eventsMap[][2], std::vector<double> channelTdiff[numOfChannels])
+ {
+ long msgcnt = 0;
+
+ double t0=500000; //variable to store the time of the first particle
+ int ch_ref=0; //the chid for the first hit
+
+ // Set all passed variables to 0
+ for (unsigned n=0; n<numOfChannels; n++) {
+ hits[n][0] = hits[n][1] = hits[n][2] = 0.;
+ eventsMap[n][0] = 0; // number of leading edges for each channel
+ eventsMap[n][1] = 0; // number of trailing edges
+ channelTdiff[n].clear();
+ }
+
+ if (!sub) return 0;
+
+ for (unsigned cnt=0;cnt<sub->Size();cnt++)
+ {
+ const hadaq::MessageFloat& msg = sub->msg(cnt);
+ unsigned chid = msg.getCh();
+ unsigned edge = msg.getEdge(); // 0 -leading, 1 - trailing
+
+
+ FillH1(hNumEntries, chid+shift);
+
+ double tm = msg.stamp;
+ if (makeTdiffCorr)
+ tm = tm + Timecorr[chid+shift];
+
+// if(shift == 60)
+// tm += Timecorr[42];
+
+// if ((chid > numOfChannels - 1) || (edge > 1)) printf("ERROR !!!! %u %u\n", chid, edge);
+ if ((chid > 26) || (edge > 1)) printf("ERROR !!!! %u %u\n", chid, edge);
+ if ((chid < numOfChannels) && (edge < 2))
+ {
+ hits[chid][edge] = tm; // this is for the last hit
+ msgcnt++;
+
+ if(edge == 0){
+ channelTdiff[chid].push_back(tm);
+ eventsMap[chid][0] += 1;
+ } else if (edge == 1)
+ eventsMap[chid][1] += 1;
+
+ }
+ if ( (chid < numOfChannels) && (edge==0) && (t0>hits[chid][edge]) )
+ {
+ ch_ref=chid;
+ t0=hits[chid][edge]; //save the new reference
+ }
+ }
+
+ // fill ToT
+ int Tot1=0;
+// int Tot2=0;
+
+ double TW_corr = 0;
+
+ //Time Walk correction mulitpliers
+ double TDC12 = 3.0;
+ double TDC34 = 3.0;
+ // ToT values have to be mulitplied by those values to aquire proper bin number in TW correction histogram
+ // TDC 1 and TDC2 have range 0-15 and 150 bins -> mulitplier 10.0
+ // TDC3 and TDC4 have range 0-50 and 150 bins -> mulitplier 3.0
+
+ //TW correction and ToT
+ for (int ch=1; ch < numOfChannels; ch++)
+ {
+ //Calculate ToT
+ hits[ch][2]=hits[ch][1]-hits[ch][0];
+// if(hits[ch][2] < 0)
+// return -999;
+//
+ // Accessing bins with TW correction
+ Tot1 = int(floor(TDC12*(hits[ch][2])));
+// Tot2 = int(floor(TDC34*(hits[ch][2])));
+
+ int index = shift / 20 + 1; // CHANGE IF SHIFT CHANGED
+ if (makeTWCorrParameters && hits[ch][0]!=0 && hits[ch][1]!=0 && correlationScheme[index][ch][1] != 0) // TWalk correction with the generated Parameterfiles if both t1 and t0 are not 0 (recording went fine) and any channel is correlated with it (Not a broken channel)
+ {
+ if(int(floor(hits[ch][1] - hits[ch][0])) < 0) {
+ hits[ch][1] = 0;
+ hits[ch][0] = 0;
+ } else {
+ TW_corr = ((TH2F*)TW_histo[index])->GetBinContent(ch,Tot1+2);
+ hits[ch][1] = hits[ch][1] - TW_corr;
+ hits[ch][0] = hits[ch][0] - TW_corr;
+ }
+ }
+
+
+
+
+
+ if (makeTWCorrFunction && hits[ch][0]!=0 && hits[ch][1]!=0 && correlationScheme[index][ch][1] != 0) //TWalk correction with the generated parameters via fit functions if both t1 and t0 are not 0 (recording went fine) and any channel is correlated with it (Not a broken channel)
+ {
+ if(int(floor(hits[ch][1] - hits[ch][0])) < 0) {
+ hits[ch][1] = 0;
+ hits[ch][0] = 0;
+ } else {
+ double param0=TWalkCorrFunction[index][ch][0];
+ double param1=TWalkCorrFunction[index][ch][1];
+ double param2=TWalkCorrFunction[index][ch][2];
+ double param3=TWalkCorrFunction[index][ch][3];
+ if (makeTdiffCorrFunction){
+ param3=param3 + Timecorr[ ch+shift];
+ if(shift == 60)
+ param3 += Timecorr[42];
+
+ }
+
+
+ if (param0!=0 || param1!=0 || param2!=0 || param3!=0)
+ TW_corr = param0/(pow(hits[ch][2],param1))+param2*hits[ch][2]+param3 ;
+ else TW_corr =0;
+ hits[ch][1] = hits[ch][1] + TW_corr;
+ hits[ch][0] = hits[ch][0] + TW_corr;
+ }
+ }
+
+
+
+
+
+
+
+ if ( (msgcnt>=4) && (ch!=ch_ref) ) //at least two particles
+ {
+ FillH2(TDTF[index], hits[ch][0]-t0, ch_ref);
+ }
+
+ if(hits[ch][0]!=0)
+ FillH2(ToTAllChannels, shift+ch, hits[ch][2]);
+
+ }
+
+ return msgcnt;
+ }
+
+ void FillTDiff(int referenceTDC, int correlatedTDC, int shift)
+ {
+ int refChannel, corrChannel;
+ double *refData, *corrData;
+ double refChCutLower = 0, refChCutUpper = 0,corrChCutLower = 0, corrChCutUpper = 0;
+
+ for (int n1=1; n1<numOfChannels; n1++) {
+ refChannel = correlationScheme[referenceTDC][n1][0];
+
+ if(makeClusterRemovalBothCh || makeRescale)
+ refData = fHitsCleared[referenceTDC][refChannel];
+ else
+ refData = fHits[referenceTDC][refChannel];
+
+ refChCutLower = makeCutsBothChannels ? cutValuesLower[referenceTDC][refChannel] : RefCutMin;
+ refChCutUpper = makeCutsBothChannels ? cutValuesUpper[referenceTDC][refChannel] : RefCutMax;
+
+ if (refChannel != 0 && refData[2] > refChCutLower && refData[2] < refChCutUpper ) {
+ for (int n2=1; n2<numOfChannels; n2++) {
+ corrChannel= correlationScheme[correlatedTDC][n2][0];
+
+ if(!makeClusterRemovalCorrCh && !makeClusterRemovalBothCh && !makeRescale)
+ corrData = fHits[correlatedTDC][corrChannel];
+ else
+ corrData = fHitsCleared[correlatedTDC][corrChannel];
+
+ corrChCutLower = (makeCutsCorrChannel || makeCutsBothChannels) ? cutValuesLower[correlatedTDC][corrChannel] : CorrCutMin,
+ corrChCutUpper = (makeCutsCorrChannel || makeCutsBothChannels) ? cutValuesUpper[correlatedTDC][corrChannel] : CorrCutMax;
+
+ if (corrChannel != 0 && corrData[2] > corrChCutLower && corrData[2] < corrChCutUpper && abs(corrData[0] - refData[0]) <TDiffCut) {
+ if(corrChannel!=refChannel)
+ FillH2(hChMap[referenceTDC], corrChannel+shift, refChannel);
+ if (referenceTDC!= correlatedTDC)
+ FillH2(TDiffAllChannels[referenceTDC][n1], corrChannel+shift, corrData[0] - refData[0]); // n1 because of correlationScheme of TDC3
+ }
+ }
+ }
+ }
+ }
+
+ void FillToT(int referenceTDC, int correlatedTDC)
+ {
+ int refChannel, corrChannel;
+ double *refData, *corrData;
+ double refChCutLower = 0, refChCutUpper = 0,corrChCutLower = 0, corrChCutUpper = 0;
+
+ for (int i=1; i<numOfChannels; i++){
+ refChannel= correlationScheme[referenceTDC][i][0];
+
+ if(makeClusterRemovalBothCh || makeRescale)
+ refData = fHitsCleared[referenceTDC][refChannel];
+ else
+ refData = fHits[referenceTDC][refChannel];
+
+ refChCutLower = makeCutsBothChannels ? cutValuesLower[referenceTDC][refChannel] : RefCutMin;
+ refChCutUpper = makeCutsBothChannels ? cutValuesUpper[referenceTDC][refChannel] : RefCutMax;
+
+ if(refChannel != 0 && refData[2]>refChCutLower && refData[2]<refChCutUpper){
+ corrChannel= correlationScheme[referenceTDC][i][1];
+
+ if(!makeClusterRemovalCorrCh && !makeClusterRemovalBothCh && !makeRescale)
+ corrData = fHits[correlatedTDC][corrChannel];
+ else
+ corrData = fHitsCleared[correlatedTDC][corrChannel];
+
+ corrChCutLower = (makeCutsCorrChannel || makeCutsBothChannels) ? cutValuesLower[correlatedTDC][corrChannel] : CorrCutMin,
+ corrChCutUpper = (makeCutsCorrChannel || makeCutsBothChannels) ? cutValuesUpper[correlatedTDC][corrChannel] : CorrCutMax;
+
+ if(corrData[2]>corrChCutLower && corrData[2]<corrChCutUpper && abs(corrData[0] - refData[0]) < TDiffCut)
+ FillH2(ToTvsTDiff[referenceTDC][i], refData[2], corrData[0] - refData[0]);
+
+// if (refData[2]<28)
+// cout << "TDC"<< referenceTDC<< " channel " << refChannel<< " T0: "<< refData[0] << " T1: " << refData[1]<< " ToT: "<< refData[2]<<endl;
+// FillH2(ToTvsTDiff[referenceTDC][i], fHits[referenceTDC][refChannel][2], corrData[0] - refData[0]);
+ }
+ }
+ }
+
+// void FillSubevents()
+ bool FillSubevents()
+ {
+ int channel, firedChannelsCounter = 0;
+ bool spark = false;
+ bool edgeMissed = false;
+ for(int i = 1; i < numOfTDCs; i++){ // TDC number
+ spark = false;
+
+ for(int j = 1; j < numOfChannels; j++) // channel number
+ if(correlationScheme[i][j][0] != 0){
+ channel = correlationScheme[i][j][0];
+
+ // Number of subevents per channel
+ // FillH2(hChEvents[i], channel, fSubEvents[i][channel]);
+ if (fSubEvents[i][channel][0]!=0 || fSubEvents[i][channel][1]!=0)
+ FillH2(hEdges[i][channel], fSubEvents[i][channel][0], fSubEvents[i][channel][1]);
+
+ //Multiplicity
+ if(fSubEvents[i][channel][0] > 3 || fSubEvents[i][channel][1] > 3 || abs(fSubEvents[i][channel][0] - fSubEvents[i][channel][1]) > 0)
+ spark = true;
+ if (abs(fSubEvents[i][channel][0] - fSubEvents[i][channel][1]) > 0)
+ edgeMissed=true;
+//
+ }
+
+ if (spark){
+ firedChannelsCounter = 0;
+
+ for(int j = 1; j < numOfChannels; j++){
+ channel = correlationScheme[i][j][0];
+
+ if(fSubEvents[i][channel][0] > 0 || fSubEvents[i][channel][1] > 0)
+ firedChannelsCounter += 1;
+ }
+
+ FillH1(sparkMultiplicity[i], firedChannelsCounter);
+// if (i==2 && firedChannelsCounter){
+// for (int j=1; j<33; j++)
+// if (correlationScheme[i][j][0] !=0)
+// cout << fHits[i][j][2] << " ";
+// cout << eventnumber <<endl;
+// }
+ }
+ }
+ return edgeMissed;
+ }
+
+ void FillDetectorHits(int detector)
+ {
+ //Set all to 0
+ for(int j = 0; j < 18; j++)
+ for(int k = 0; k < 3; k++)
+ fHitsDetector[detector][j][k] = 0;
+
+ for(int i = 1; i < numOfTDCs; i++){
+ for(int j = 0; j < numOfChannels; j++){
+ if(correlationScheme[i][j][2] == detector){
+ for(int k = 0; k < 3; k++)
+ if(k<2)
+ fHitsDetector[detector][correlationScheme[i][j][3]][k] = fHits[i][correlationScheme[i][j][0]][k];
+ else
+ fHitsDetector[detector][correlationScheme[i][j][3]][k] = fHits[i][correlationScheme[i][j][0]][k] * rescale[(i-1)*20 + correlationScheme[i][j][0]];
+ }
+ }
+ }
+ }
+
+ std::vector<std::vector<int>> findClusters(int detector)
+ {
+ std::vector<std::vector<int>> result;
+ result.clear();
+
+ int firstPosition = 0;
+ int lastPosition = 0;
+ int maxPosition = 0;
+ int currentPosition = 1;
+
+ std::vector<int> clusterData;
+ while(currentPosition<numOfChannels) {
+ clusterData.clear();
+
+ for (int k=currentPosition; k<17; k++) {
+ if(fHitsDetector[detector][k][2] != 0 ) {
+ if (fHitsDetector[detector][k-1][2] == 0) {
+ firstPosition = k ;
+ }
+
+ if (fHitsDetector[detector][k][2] > fHitsDetector[detector][maxPosition][2]) {
+ maxPosition = k;
+ }
+
+ if ((k < 16) && (fHitsDetector[detector][k+1][2] == 0) ) {
+ lastPosition = k + 1;
+ currentPosition = k;
+ }
+ else if( k == 16) {
+ lastPosition = k+1; // cluster size = last position - first position, index 17 is OK!
+ currentPosition = k;
+ }
+ }
+
+ if( lastPosition!=0) break;
+ }
+
+ if(firstPosition != 0 && lastPosition !=0){
+ clusterData.push_back(firstPosition);
+ clusterData.push_back(maxPosition);
+ clusterData.push_back(lastPosition);
+ result.push_back(clusterData);
+ }
+
+ currentPosition++;
+ maxPosition = 0;
+ lastPosition = 0;
+ firstPosition = 0;
+ }
+
+ return result;
+ }
+
+ void RemoveClusters(int detector)
+ {
+ auto detectorClustersLocal = detectorClusters[detector];
+ for(auto subCluster : detectorClustersLocal )
+ for(int i = subCluster[0]; i < subCluster[2]; i++){
+ if ( i != subCluster[1] ) {
+ fHitsDetector[detector][i][0] = 0;
+ fHitsDetector[detector][i][1] = 0;
+ fHitsDetector[detector][i][2] = 0;
+ }
+ }
+ }
+
+ void createClearedfHits()
+ {
+ int* coordinates;
+ for(int i = 0; i < numOfTDCs; i++)
+ for(int j = 0; j < numOfChannels; j++)
+ for(int k = 0; k < 3; k++)
+ fHitsCleared[i][j][k] = 0.;
+
+ for(int i = 1; i < numOfDetectors; i++) {
+ for(int j = 1; j < 17; j++){
+ if(fHitsDetector[i][j][2] != 0){
+ coordinates = DetectorHitsToHits(i, j);
+
+ for(int k = 0; k < 3; k++){
+ fHitsCleared[coordinates[0]][coordinates[1]][k] = fHits[coordinates[0]][coordinates[1]][k];
+
+ if(makeRescale)
+ fHitsCleared[coordinates[0]][coordinates[1]][2] = fHitsDetector[i][j][2];
+ }
+ }
+ }
+ }
+ }
+
+ int* DetectorHitsToHits(int detector, int position)
+ {
+ return detectorScheme[detector][position];
+ }
+
+ void FillClusterDetector (int detector) {
+ int multiplicity=0;
+ int clusterSize = 0;
+
+ for (unsigned int i=0; i< detectorClusters[detector].size(); i++ ){
+ clusterSize = detectorClusters[detector][i][2]-detectorClusters[detector][i][0];
+ multiplicity += clusterSize ;
+ if(clusterSize>0)
+ FillH1 (clusterSizeDetector[detector], clusterSize);
+ }
+
+ if(multiplicity>0)
+ FillH1 (multiplicityDetector[detector], multiplicity);
+ }
+
+ void FillNeighbouringCrosstalk (int detector){
+ int clusterSize=0;
+ int otherDetector =0;
+ if (detector == 1)
+ otherDetector =2;
+ else otherDetector=1;
+
+ for (int i=1; i< 17; i++){ // 17 is num of channels connected to detector
+
+ auto position = DetectorHitsToHits(detector, i); // [TDC, channel]
+ auto otherPosition = DetectorHitsToHits(otherDetector, i); // output->[(other)TDC, channel]
+
+ int channel= position[1];
+ int TDC=position[0];
+ int otherTDC= otherPosition[0];
+ int otherChannel= correlationScheme[TDC][channel][1];
+
+ if (fHitsDetector[detector][i][2]>0 && fHits[otherTDC][otherChannel][2] >0 && TMath::Abs(fHitsDetector[detector][i][0]-fHits[otherTDC][otherChannel][0]) <15 ){// this checks if signals used for filling are really coming from a particle traversing both detectors, i.e. inducing a signal in both detectors, which a TDiff of smaller than 15 ns
+ clusterSize +=1;
+ if (fHitsDetector[detector][i+1][2]>0 && TMath::Abs(fHitsDetector[detector][i][0] - fHitsDetector[detector][i+1][0]) <10)
+ clusterSize +=1;
+ if (fHitsDetector[detector][i-1][2]>0 && TMath::Abs(fHitsDetector[detector][i][0] - fHitsDetector[detector][i-1][0]) <10)
+ clusterSize +=1;
+
+ FillH2 (NeighbourCrosstalk[position[0]][position[1]], fHitsDetector[detector][i][2] ,clusterSize); // TODO correlation between channels for hit-> TDiff around 10 ns, not more
+
+ if( fHitsDetector[detector][i][2]>fHitsDetector[detector][i+1][2] && fHitsDetector[detector][i][2]>fHitsDetector[detector][i-1][2])
+ FillH2 (LocalMaxNeighbourCrosstalk[position[0]][position[1]], fHitsDetector[detector][i][2] ,clusterSize);
+
+ clusterSize = 0;
+ }
+ }
+ }
+
+ void FillTDiff_Left_Right (int detector){
+ int otherDetector =0;
+ if (detector == 1)
+ otherDetector =2;
+ else otherDetector=1;
+
+ for (int i=1; i< 17; i++){
+
+ auto position = DetectorHitsToHits(detector, i); // [TDC, channel]
+ auto otherPosition = DetectorHitsToHits(otherDetector, i); // output->[(other)TDC, channel]
+
+ int channel= position[1];
+ int TDC=position[0];
+ int otherTDC= otherPosition[0];
+ int otherChannel= correlationScheme[TDC][channel][1];
+
+ if (fHitsDetector[detector][i][2]>0 && fHitsDetector[detector][i+1][2]>0 && fHits[otherTDC][otherChannel][2] >0 && TMath::Abs(fHitsDetector[detector][i][0]-fHits[otherTDC][otherChannel][0]) <15 ) {
+ FillH2(Right_TDiff_Crosstalk[position[0]][position[1]], fHitsDetector[detector][i][2], fHitsDetector[detector][i][0]-fHitsDetector[detector][i+1][0]);
+
+ if( fHitsDetector[detector][i][2]>fHitsDetector[detector][i+1][2] && fHitsDetector[detector][i][2]>fHitsDetector[detector][i-1][2] )
+ FillH2(LocalMax_Right_TDiff_Crosstalk[position[0]][position[1]], fHitsDetector[detector][i][2], fHitsDetector[detector][i][0]-fHitsDetector[detector][i+1][0]);
+ }
+ if (fHitsDetector[detector][i][2]>0 && fHitsDetector[detector][i-1][2]>0 && fHits[otherTDC][otherChannel][2] >0 && TMath::Abs(fHitsDetector[detector][i][0]-fHits[otherTDC][otherChannel][0]) <15 ){
+ FillH2(Left_TDiff_Crosstalk[position[0]][position[1]], fHitsDetector[detector][i][2], fHitsDetector[detector][i][0]-fHitsDetector[detector][i-1][0]);
+ // FillH2(Left_TDiff_Crosstalk[position[0]][position[1]], fHitsDetector[detector][i][2], fHitsDetector[detector][i-1][2]);
+ if( fHitsDetector[detector][i][2]>fHitsDetector[detector][i-1][2] && fHitsDetector[detector][i][2]>fHitsDetector[detector][i+1][2])
+ FillH2(LocalMax_Left_TDiff_Crosstalk[position[0]][position[1]], fHitsDetector[detector][i][2], fHitsDetector[detector][i][0]-fHitsDetector[detector][i-1][0]);
+ }
+
+ channel =0;
+ TDC =0;
+ otherTDC=0;
+ otherChannel=0;
+ }
+ }
+
+ void FillToT_Left_Right (int detector){
+ int otherDetector =0;
+ if (detector == 1)
+ otherDetector =2;
+ else otherDetector=1;
+
+ for (int i=1; i< 17; i++){
+ auto position = DetectorHitsToHits(detector, i); // [TDC, channel]
+ auto otherPosition = DetectorHitsToHits(otherDetector, i); // output->[(other)TDC, channel]
+
+ int channel= position[1];
+ int TDC=position[0];
+ int otherTDC= otherPosition[0];
+ int otherChannel= correlationScheme[TDC][channel][1];
+
+ if (fHitsDetector[detector][i][2]>0 && fHitsDetector[detector][i+1][2]>0 && fHits[otherTDC][otherChannel][2] >0 && TMath::Abs(fHitsDetector[detector][i][0]-fHits[otherTDC][otherChannel][0]) <15 ) {
+ FillH2(Right_ToT_Crosstalk[position[0]][position[1]], fHitsDetector[detector][i][2],fHitsDetector[detector][i+1][2]);
+
+ if( fHitsDetector[detector][i][2]>fHitsDetector[detector][i+1][2] && fHitsDetector[detector][i][2]>fHitsDetector[detector][i-1][2] )
+ FillH2(LocalMax_Right_ToT_Crosstalk[position[0]][position[1]], fHitsDetector[detector][i][2], fHitsDetector[detector][i+1][2]);
+ }
+ if (fHitsDetector[detector][i][2]>0 && fHitsDetector[detector][i-1][2]>0 && fHits[otherTDC][otherChannel][2] >0 && TMath::Abs(fHitsDetector[detector][i][0]-fHits[otherTDC][otherChannel][0]) <15 ){
+ FillH2(Left_ToT_Crosstalk[position[0]][position[1]], fHitsDetector[detector][i][2], fHitsDetector[detector][i-1][2]);
+
+ if( fHitsDetector[detector][i][2]>fHitsDetector[detector][i-1][2] && fHitsDetector[detector][i][2]>fHitsDetector[detector][i+1][2])
+ FillH2(LocalMax_Left_ToT_Crosstalk[position[0]][position[1]], fHitsDetector[detector][i][2], fHitsDetector[detector][i-1][2]);
+ }
+
+ channel =0;
+ TDC =0;
+ otherTDC=0;
+ otherChannel=0;
+ }
+ }
+
+ void CreateDetectorScheme(){
+ // Iterate over detectorScheme and fill it
+ for(int i = 1; i < numOfDetectors; i++)
+ for(int j = 1; j < 17; j++) {
+
+ // Iterate over correlationScheme
+ for(int k = 1; k < numOfTDCs; k++) // TDC number
+ for(int z = 1; z < numOfChannels; z++) // channel
+ if(correlationScheme[k][z][2] == i && correlationScheme[k][z][3] == j){
+ detectorScheme[i][j][0] = k;
+ detectorScheme[i][j][1] = z;
+ }
+ }
+ }
+
+ //Main function of this analysis
+ virtual bool Process(base::Event* ev)
+ {
+ eventnumber++;
+// if(eventnumber > 2000000){
+// cout << "2 000 000 events passed, analysis stopped" << endl;
+// return false;
+// }
+
+ hadaq::TdcSubEventFloat* sub1 =
+ dynamic_cast<hadaq::TdcSubEventFloat*> (ev->GetSubEvent(fTdcId1));
+
+ hadaq::TdcSubEventFloat* sub2 =
+ dynamic_cast<hadaq::TdcSubEventFloat*> (ev->GetSubEvent(fTdcId2));
+
+ hadaq::TdcSubEventFloat* sub3 =
+ dynamic_cast<hadaq::TdcSubEventFloat*> (ev->GetSubEvent(fTdcId3));
+
+ hadaq::TdcSubEventFloat* sub4 =
+ dynamic_cast<hadaq::TdcSubEventFloat*> (ev->GetSubEvent(fTdcId4));
+
+ long num1 = 0;
+ long num2 = 0;
+ long num3 = 0;
+ long num4 = 0;
+ long num = 0;
+
+ num1 = Unpack(sub1, 0, fHits[1], fSubEvents[1], fSubTdiff[1]);
+ num2 = Unpack(sub2, 20, fHits[2], fSubEvents[2], fSubTdiff[2]);
+ num3 = Unpack(sub3, 40, fHits[3], fSubEvents[3], fSubTdiff[3]);
+ num4 = Unpack(sub4, 60, fHits[4], fSubEvents[4], fSubTdiff[4]);
+
+ num=num1+num2+num3+num4;
+
+ // when return false, event processing is cancelled
+ if (num == 0 ) return false;
+
+ FillH1(hNumHits, num);
+
+ for(int i = 1; i < numOfDetectors; i++) {
+ FillDetectorHits(i);
+
+//
+// if(makeCrosstalkAnalysis){
+// FillNeighbouringCrosstalk(i);
+// FillTDiff_Left_Right(i);
+// FillToT_Left_Right(i);
+// }
+//
+ if (makeClusterAnalysis || makeClusterRemovalCorrCh || makeClusterRemovalBothCh){
+ detectorClusters[i] = findClusters(i);
+//
+ if (makeClusterAnalysis )
+ FillClusterDetector(i);
+ if (makeClusterRemovalCorrCh || makeClusterRemovalBothCh)
+ RemoveClusters(i);
+ }
+ }
+//
+ createClearedfHits();
+
+ //Fill ToT
+ FillToT(1,1);
+// FillToT(1,2);
+// FillToT(2,1);
+// FillToT(3,4);
+// FillToT(4,3);
+
+ //Fill TDiff histograms
+ for(int i=1; i < numOfTDCs; i++) // Reference TDC
+ for(int j = 1; j <numOfTDCs; j++) // Correlated TDC
+ // if(i != j)
+ FillTDiff(i, j, (j-1) * 20);
+
+
+
+ return true;
+ };
+};
+
+void second()
+{
+ new SecondProc("TEST", "TDC_2600", "TDC_2601", "TDC_2602", "TDC_2603");
+}
jspc29 / daqtools.git / commitdiff