diff --git a/macro/psd/build_events.cpp b/macro/psd/build_events.cpp
index e640531619985e65e674ab2b49bb9d10c9e24c5d..f5ecb210fbf21a081abfb389987937c7214940e9 100644
--- a/macro/psd/build_events.cpp
+++ b/macro/psd/build_events.cpp
@@ -1,20 +1,20 @@
-#include <string.h>
-#include <math.h>
+#include <TBrowser.h>
+#include <TCanvas.h>
+#include <TFile.h>
+#include <TGraph.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TMath.h>
+#include <TObjArray.h>
+#include <TObject.h>
+#include <TString.h>
+#include <TTree.h>
-#include <iostream>
#include <fstream>
-#include<TTree.h>
-#include<TString.h>
-#include<TMath.h>
-#include<TObject.h>
-#include<TFile.h>
-#include<TH1F.h>
-#include<TH2F.h>
-#include<TCanvas.h>
-#include<TBrowser.h>
-#include<TObjArray.h>
-#include<TGraph.h>
+#include <iostream>
+#include <math.h>
+#include <string.h>
#include "event_data_struct.h"
@@ -22,52 +22,52 @@ using namespace std;
// root -l -b -q 'GBTreaderEventBuild.cpp("DspTree.root", "EventTree.root")'
-int build_events(TString inFileName, TString outFileName, int dsps_to_read = -1) {
- TString srcDir = gSystem->Getenv("VMCWORKDIR");
- const int total_channels = 32;
- const int gate_length = 60; //in ns
-
- TFile* inFile = TFile::Open(inFileName.Data(), "READONLY");
- TTree* inTree = (TTree*) inFile->Get("cbmsim");
- if (inTree == nullptr) {
- cout << "tree not found" << endl;
- return -1;
- }
-
- std::vector<CbmPsdDsp>* fPsdDspVect = new std::vector<CbmPsdDsp>;
- inTree->SetBranchAddress("PsdDsp", &fPsdDspVect);
- if (fPsdDspVect == nullptr) {
- cout << "no psd dsps vector in tree" << endl;
- return -2;
- }
- Int_t Nhits = 0;
- ULong64_t EventTime = 0;
-
-
- printf("\n\n====================================================\n");
- printf("================= TSA DATA READOUT =================\n");
- printf("====================================================\n\n\n");
-
- TFile *result_file = new TFile(outFileName.Data(), "RECREATE");
- event_data_struct *event_struct = new event_data_struct[total_channels];
- TTree *event_tree = new TTree("EventTree", "EventTree");
- event_tree->Branch("Nhits",&Nhits,"Nhits/I");
- event_tree->Branch("EventTime",&EventTime,"EventTime/l");
- for(Int_t ch = 0; ch < total_channels; ch++)
- {
- (event_struct+ch)->CreateBranch(event_tree, ch);
- (event_struct+ch)->reset();
- }
-
- bool first_all_sync = false;
- double hit_time = 0.;
- double ev_start_time = 0.;
- double prev_hit_time = 0.;
- double timestump = 0.;
- double prev_timestump = 0.;
- std::vector<CbmPsdDsp> event_dsp_vect;
-
-/*
+int build_events(TString inFileName, TString outFileName, int dsps_to_read = -1)
+{
+ TString srcDir = gSystem->Getenv("VMCWORKDIR");
+ const int total_channels = 32;
+ const int gate_length = 60; //in ns
+
+ TFile* inFile = TFile::Open(inFileName.Data(), "READONLY");
+ TTree* inTree = (TTree*) inFile->Get("cbmsim");
+ if (inTree == nullptr) {
+ cout << "tree not found" << endl;
+ return -1;
+ }
+
+ std::vector<CbmPsdDsp>* fPsdDspVect = new std::vector<CbmPsdDsp>;
+ inTree->SetBranchAddress("PsdDsp", &fPsdDspVect);
+ if (fPsdDspVect == nullptr) {
+ cout << "no psd dsps vector in tree" << endl;
+ return -2;
+ }
+ Int_t Nhits = 0;
+ ULong64_t EventTime = 0;
+
+
+ printf("\n\n====================================================\n");
+ printf("================= TSA DATA READOUT =================\n");
+ printf("====================================================\n\n\n");
+
+ TFile* result_file = new TFile(outFileName.Data(), "RECREATE");
+ event_data_struct* event_struct = new event_data_struct[total_channels];
+ TTree* event_tree = new TTree("EventTree", "EventTree");
+ event_tree->Branch("Nhits", &Nhits, "Nhits/I");
+ event_tree->Branch("EventTime", &EventTime, "EventTime/l");
+ for (Int_t ch = 0; ch < total_channels; ch++) {
+ (event_struct + ch)->CreateBranch(event_tree, ch);
+ (event_struct + ch)->reset();
+ }
+
+ bool first_all_sync = false;
+ double hit_time = 0.;
+ double ev_start_time = 0.;
+ double prev_hit_time = 0.;
+ double timestump = 0.;
+ double prev_timestump = 0.;
+ std::vector<CbmPsdDsp> event_dsp_vect;
+
+ /*
{
int dsps_for_time_btw = (int) 1e5; //1e5
TH1F *h1_time_btw_hits = new TH1F("h1_time_btw", "h1_time_btw; time [ns]; Counts", 500, -50, 300);
@@ -99,142 +99,115 @@ int build_events(TString inFileName, TString outFileName, int dsps_to_read = -1)
}
*/
- time_t start_time = time(NULL);
- prev_hit_time = 0.;
- CbmPsdDsp dsp;
- for( int entr_iter = 0; entr_iter < inTree->GetEntries(); entr_iter++ )
- {
- inTree->GetEntry(entr_iter);
- Int_t nPsdDsps = fPsdDspVect->size();
- for (Int_t iDsp = 0; iDsp < nPsdDsps; iDsp++) {
- dsp = (*fPsdDspVect)[iDsp];
- if(dsp.GetSectionID() > 9) continue;
- hit_time = dsp.fdTime;
-
- if( hit_time < ev_start_time) cout<<"negative time "<< hit_time << " " << ev_start_time<<endl;
- if( hit_time < prev_hit_time) cout<<"negative time"<<endl;
-
- if( (hit_time - ev_start_time) < gate_length)
- //if( (hit_time-prev_hit_time) < gate_length )
- {
- event_dsp_vect.push_back(dsp);
- prev_hit_time = dsp.fdTime;
- }
- else
- {
- Nhits = event_dsp_vect.size();
- EventTime = (ULong64_t) ev_start_time;
-
- for (int ch_iter = 0; ch_iter < total_channels; ch_iter++) event_struct[ch_iter].reset();
-
- // Check channel duplicates in dsp vector
- std::map<uint32_t, int> countMap;
- for (auto & elem : event_dsp_vect)
- {
- auto result = countMap.insert(std::pair<uint32_t, int>(elem.GetSectionID(), 1));
- if (result.second == false) result.first->second++;
- }
-
- // Iterate over the map
- for (auto & elem : countMap)
- {
- // If frequency count is greater than 1 then its a duplicate element
- if (elem.second > 1)
- {
- std::cout << "Replacement in channel " << elem.first << " :: total times " << elem.second << " :: Event vector size "<< event_dsp_vect.size() << std::endl;
-
- int counter = 0;
- cout.precision(32);
- for (auto & PsdDsp : event_dsp_vect)
- {
- cout<<"# "<<counter<<" dsp_time-event_begin_time "<<PsdDsp.fdTime-ev_start_time<<endl;
- cout<< " prev_timestump " << prev_timestump << " timestamp "<< timestump <<endl;
- cout<<PsdDsp.GetSectionID()<<" "<<PsdDsp.fdEdep<<" "<<PsdDsp.fdTime<<endl;
- counter++;
- }
-
- cout<<"End of microslice"<<endl;
-
-
- }
- }
-
-
- for (auto & PsdDsp : event_dsp_vect)
- {
- int hit_channel = PsdDsp.GetSectionID();
- event_data_struct &result_event = event_struct[hit_channel];
-
- result_event.EdepFPGA = PsdDsp.fdEdep;
- result_event.EdepWfm = PsdDsp.fdEdepWfm;
- result_event.Ampl = PsdDsp.fdAmpl;
- result_event.Minimum = PsdDsp.fuMinimum;
- result_event.ZL = PsdDsp.fuZL;
- result_event.Time = PsdDsp.fdTime;
- result_event.TimeMax = PsdDsp.fuTimeMax;
- result_event.TimeInEvent = PsdDsp.fdTime-ev_start_time;
-
- result_event.FitEdep = PsdDsp.fdFitEdep;
- result_event.FitAmpl = PsdDsp.fdFitAmpl;
- result_event.FitZL = PsdDsp.fdFitZL;
- result_event.FitR2 = PsdDsp.fdFitR2;
- result_event.FitTimeMax = PsdDsp.fdFitTimeMax;
-
- result_event.EventTime = (ULong64_t) ev_start_time;
-
- }
-
- if(Nhits>0) event_tree->Fill(); //skipping first dsp
- event_dsp_vect.clear();
- //countMap.clear();
-
- event_dsp_vect.push_back(dsp);
- ev_start_time = hit_time;
- prev_hit_time = hit_time;
-
-
-
- }
-
- }
-
- if(dsps_to_read>0 && event_tree->GetEntries() >= dsps_to_read) break;
-
-
- if((entr_iter%100) == 0)
- {
- time_t current_time = time(NULL);
- Int_t proc_sec = difftime(current_time,start_time);
- Float_t percents = (float)entr_iter/(inTree->GetEntries());
- Int_t time_est = (percents == 0)?0:proc_sec/percents*(1. - percents);
- Float_t proc_rate = (float)proc_sec/entr_iter*1000000./60.;
-
- printf("Processed entries: %i (%5.1f%%); [pas %3.0dm %2.0is] [est %3.0dm %2.0is] [%.1f min [%.1fh]/1M ev]\r",
- entr_iter, (percents*100.),
- (proc_sec/60), proc_sec%60,
- (time_est/60), time_est%60,
- proc_rate, (proc_rate/60.));
- cout<<flush;
- }
-
- prev_timestump = timestump;
-
- }
-
-
- event_tree->Write();
- result_file->Close();
-
- return 0;
-
-
-
-
-
-
+ time_t start_time = time(NULL);
+ prev_hit_time = 0.;
+ CbmPsdDsp dsp;
+ for (int entr_iter = 0; entr_iter < inTree->GetEntries(); entr_iter++) {
+ inTree->GetEntry(entr_iter);
+ Int_t nPsdDsps = fPsdDspVect->size();
+ for (Int_t iDsp = 0; iDsp < nPsdDsps; iDsp++) {
+ dsp = (*fPsdDspVect)[iDsp];
+ if (dsp.GetSectionID() > 9) continue;
+ hit_time = dsp.fdTime;
+
+ if (hit_time < ev_start_time) cout << "negative time " << hit_time << " " << ev_start_time << endl;
+ if (hit_time < prev_hit_time) cout << "negative time" << endl;
+
+ if ((hit_time - ev_start_time) < gate_length)
+ //if( (hit_time-prev_hit_time) < gate_length )
+ {
+ event_dsp_vect.push_back(dsp);
+ prev_hit_time = dsp.fdTime;
+ }
+ else {
+ Nhits = event_dsp_vect.size();
+ EventTime = (ULong64_t) ev_start_time;
+
+ for (int ch_iter = 0; ch_iter < total_channels; ch_iter++)
+ event_struct[ch_iter].reset();
+
+ // Check channel duplicates in dsp vector
+ std::map<uint32_t, int> countMap;
+ for (auto& elem : event_dsp_vect) {
+ auto result = countMap.insert(std::pair<uint32_t, int>(elem.GetSectionID(), 1));
+ if (result.second == false) result.first->second++;
+ }
+
+ // Iterate over the map
+ for (auto& elem : countMap) {
+ // If frequency count is greater than 1 then its a duplicate element
+ if (elem.second > 1) {
+ std::cout << "Replacement in channel " << elem.first << " :: total times " << elem.second
+ << " :: Event vector size " << event_dsp_vect.size() << std::endl;
+
+ int counter = 0;
+ cout.precision(32);
+ for (auto& PsdDsp : event_dsp_vect) {
+ cout << "# " << counter << " dsp_time-event_begin_time " << PsdDsp.fdTime - ev_start_time << endl;
+ cout << " prev_timestump " << prev_timestump << " timestamp " << timestump << endl;
+ cout << PsdDsp.GetSectionID() << " " << PsdDsp.fdEdep << " " << PsdDsp.fdTime << endl;
+ counter++;
+ }
+
+ cout << "End of microslice" << endl;
+ }
+ }
+
+
+ for (auto& PsdDsp : event_dsp_vect) {
+ int hit_channel = PsdDsp.GetSectionID();
+ event_data_struct& result_event = event_struct[hit_channel];
+
+ result_event.EdepFPGA = PsdDsp.fdEdep;
+ result_event.EdepWfm = PsdDsp.fdEdepWfm;
+ result_event.Ampl = PsdDsp.fdAmpl;
+ result_event.Minimum = PsdDsp.fuMinimum;
+ result_event.ZL = PsdDsp.fuZL;
+ result_event.Time = PsdDsp.fdTime;
+ result_event.TimeMax = PsdDsp.fuTimeMax;
+ result_event.TimeInEvent = PsdDsp.fdTime - ev_start_time;
+
+ result_event.FitEdep = PsdDsp.fdFitEdep;
+ result_event.FitAmpl = PsdDsp.fdFitAmpl;
+ result_event.FitZL = PsdDsp.fdFitZL;
+ result_event.FitR2 = PsdDsp.fdFitR2;
+ result_event.FitTimeMax = PsdDsp.fdFitTimeMax;
+
+ result_event.EventTime = (ULong64_t) ev_start_time;
+ }
+
+ if (Nhits > 0) event_tree->Fill(); //skipping first dsp
+ event_dsp_vect.clear();
+ //countMap.clear();
+
+ event_dsp_vect.push_back(dsp);
+ ev_start_time = hit_time;
+ prev_hit_time = hit_time;
+ }
+ }
+
+ if (dsps_to_read > 0 && event_tree->GetEntries() >= dsps_to_read) break;
+
+
+ if ((entr_iter % 100) == 0) {
+ time_t current_time = time(NULL);
+ Int_t proc_sec = difftime(current_time, start_time);
+ Float_t percents = (float) entr_iter / (inTree->GetEntries());
+ Int_t time_est = (percents == 0) ? 0 : proc_sec / percents * (1. - percents);
+ Float_t proc_rate = (float) proc_sec / entr_iter * 1000000. / 60.;
+
+ printf("Processed entries: %i (%5.1f%%); [pas %3.0dm %2.0is] [est %3.0dm %2.0is] [%.1f min [%.1fh]/1M ev]\r",
+ entr_iter, (percents * 100.), (proc_sec / 60), proc_sec % 60, (time_est / 60), time_est % 60, proc_rate,
+ (proc_rate / 60.));
+ cout << flush;
+ }
+
+ prev_timestump = timestump;
+ }
+
+
+ event_tree->Write();
+ result_file->Close();
+
+ return 0;
}
-
-
-
-
-
diff --git a/macro/psd/compare_to_sim.cpp b/macro/psd/compare_to_sim.cpp
index ff050c703b40b69f34bf60ffb149091e8f53580a..55245274f5325e7f5420cbccb32dd4e474a79c33 100644
--- a/macro/psd/compare_to_sim.cpp
+++ b/macro/psd/compare_to_sim.cpp
@@ -1,26 +1,26 @@
-#include <string.h>
-#include <math.h>
-
-#include <TROOT.h>
-#include <iostream>
-#include <fstream>
-#include<TTree.h>
-#include<TString.h>
-#include<TMath.h>
-#include<TObject.h>
-#include<TFile.h>
-#include<TH1F.h>
-#include<TH2F.h>
-#include<TCanvas.h>
-#include<TBrowser.h>
-#include<TObjArray.h>
-#include<TGraph.h>
-#include <TLatex.h>
+#include <TBrowser.h>
+#include <TCanvas.h>
#include <TChain.h>
-
-#include <TGraphErrors.h>
#include <TF1.h>
+#include <TFile.h>
+#include <TGraph.h>
+#include <TGraphErrors.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TLatex.h>
+#include <TMath.h>
+#include <TObjArray.h>
+#include <TObject.h>
+#include <TROOT.h>
+#include <TString.h>
#include <TStyle.h>
+#include <TTree.h>
+
+#include <fstream>
+#include <iostream>
+
+#include <math.h>
+#include <string.h>
#include "event_data_struct.h"
@@ -30,231 +30,220 @@ using namespace std;
// root -l -b -q 'GBTtreeShowSimple_fin.cpp("/mnt/c/Users/dmitr/cernbox/SOURCE/CALO/INR_ADC/16_cosmic_TRG09_44V_LC/Readout.root", "/mnt/c/Users/dmitr/cernbox/SOURCE/CALO/INR_ADC/16_cosmic_TRG09_44V_LC/Readout_ss", 1000)'
+Double_t SpectraCenter(TH1* hist, float nrms = 6.)
+{
+ Double_t hist_start_Xmin = hist->GetXaxis()->GetXmin();
+ Double_t hist_start_Xmax = hist->GetXaxis()->GetXmax();
+ hist->SetAxisRange(hist_start_Xmin, hist_start_Xmax);
+ Double_t histMean = hist->GetMean();
+ Double_t histRMS = hist->GetRMS();
+ Double_t histBin = hist->GetBinWidth(1);
+ if (histRMS < histBin) histRMS = histBin;
+ hist->SetAxisRange(histMean - histRMS * nrms, histMean + histRMS * nrms);
+ printf("%s m %f r %f\n", hist->GetName(), histMean, histRMS);
+ return histMean;
+}
+
+
+int compare_to_sim(TString file_name, TString result_file_name, TString calc_mode = "EdepFPGA", int events_to_read = -1,
+ int nhits = 0, int cut_option = 0)
+{
- Double_t SpectraCenter(TH1* hist, float nrms = 6.){
- Double_t hist_start_Xmin = hist->GetXaxis()->GetXmin();
- Double_t hist_start_Xmax = hist->GetXaxis()->GetXmax();
- hist->SetAxisRange(hist_start_Xmin, hist_start_Xmax);
-
- Double_t histMean = hist->GetMean();
- Double_t histRMS = hist->GetRMS();
- Double_t histBin = hist->GetBinWidth(1);
- if (histRMS < histBin) histRMS = histBin;
-
- hist->SetAxisRange(histMean-histRMS*nrms, histMean+histRMS*nrms);
- printf("%s m %f r %f\n", hist->GetName(), histMean, histRMS);
- return histMean;
+ const int total_channels = 10;
+ std::vector<int> hist_pars = {200, 0, 100};
+ result_file_name = result_file_name + Form("_by_%s", calc_mode.Data());
+ switch (cut_option) {
+ case 0: result_file_name = result_file_name + "_noCut"; break;
+ case 1: result_file_name = result_file_name + "_tCut"; break;
+ case 2: result_file_name = result_file_name + "_tR2Cut"; break;
}
+ TFile* file_sim =
+ TFile::Open("/home/nikolay/git_projects/mPSD_test_bench/data/ONi2AGeV_sim/O2AGeV_Ni_Hists_10kk.root");
+ gDirectory->cd("Rint:/");
+
+ TObjArray* canvas_array = new TObjArray;
+ canvas_array->SetName("Canvas");
+
+ TObjArray* result_array = new TObjArray;
+ result_array->SetName("calo_results");
+ result_array->SetOwner();
+
+ TChain* data_tree = new TChain;
+ data_tree->AddFile((file_name + "/EventTree").Data());
+ if (events_to_read == -1) events_to_read = data_tree->GetEntries();
+ event_data_struct* event_struct = new event_data_struct[total_channels];
+ for (Int_t ch = 0; ch < total_channels; ch++)
+ (event_struct + ch)->SetBranch(data_tree, ch);
+ TH1* th1_hist_ptr = nullptr;
+ TH2* th2_hist_ptr = nullptr;
+ th1_hist_ptr = new TH1F("EdepProfile", Form("EdepProfile"), total_channels, 0, total_channels);
+ th1_hist_ptr->SetTitle(Form("Energy profile; Chan; Energy [MeV]"));
+ result_array->Add(th1_hist_ptr);
+ TH1F* temp = new TH1F("temp", Form("temp"), total_channels, 0, total_channels);
+ data_tree->Draw(Form("channel_%i.%s>>temp", 0, calc_mode.Data()), Form("channel_%i.%s > 1", 0, calc_mode.Data()), "",
+ events_to_read);
-int compare_to_sim(TString file_name, TString result_file_name, TString calc_mode = "EdepFPGA", int events_to_read = -1, int nhits = 0, int cut_option = 0) {
-
- const int total_channels = 10;
- std::vector<int> hist_pars = {200, 0, 100};
- result_file_name = result_file_name+Form("_by_%s", calc_mode.Data());
- switch(cut_option){
- case 0:
- result_file_name = result_file_name+"_noCut";
- break;
- case 1:
- result_file_name = result_file_name+"_tCut";
- break;
- case 2:
- result_file_name = result_file_name+"_tR2Cut";
- break;
- }
-
- TFile* file_sim = TFile::Open("/home/nikolay/git_projects/mPSD_test_bench/data/ONi2AGeV_sim/O2AGeV_Ni_Hists_10kk.root");
- gDirectory->cd("Rint:/");
-
- TObjArray *canvas_array = new TObjArray;
- canvas_array->SetName ("Canvas");
-
- TObjArray *result_array = new TObjArray;
- result_array->SetName("calo_results");
- result_array->SetOwner();
-
- TChain *data_tree = new TChain;
- data_tree->AddFile( (file_name + "/EventTree").Data() );
- if(events_to_read==-1) events_to_read = data_tree->GetEntries();
-
- event_data_struct *event_struct = new event_data_struct[total_channels];
- for(Int_t ch = 0; ch < total_channels; ch++)
- (event_struct+ch)->SetBranch(data_tree, ch);
-
- TH1 *th1_hist_ptr = nullptr;
- TH2 *th2_hist_ptr = nullptr;
-
- th1_hist_ptr = new TH1F( "EdepProfile", Form("EdepProfile"), total_channels,0,total_channels);
- th1_hist_ptr->SetTitle( Form("Energy profile; Chan; Energy [MeV]") );
- result_array->Add(th1_hist_ptr);
-
- TH1F* temp = new TH1F( "temp", Form("temp"), total_channels,0,total_channels);
- data_tree->Draw( Form("channel_%i.%s>>temp", 0, calc_mode.Data()), Form("channel_%i.%s > 1", 0, calc_mode.Data()), "", events_to_read );
-
- for(int ch_iter = 0; ch_iter < total_channels; ch_iter++)
- {
- th1_hist_ptr = new TH1F( Form("Edep_ch%i", ch_iter), Form("Edep_ch%i", ch_iter), hist_pars[0],hist_pars[1],hist_pars[2] );
- th1_hist_ptr->SetTitle( Form("Charge channel %i; Edep [MeV]", ch_iter) );
-
- th2_hist_ptr = new TH2F( Form("timeMax_vs_Edep_ch%i", ch_iter), Form("timeMax_vs_Edep_ch%i", ch_iter), 32,0,32, hist_pars[0],hist_pars[1],hist_pars[2] );
- th2_hist_ptr->SetTitle( Form("TimeMax vs Edep channel %i; TimeMax [12,5ns]; Charge [mV]", ch_iter) );
-
- th2_hist_ptr = new TH2F( Form("FitR2_vs_Edep_ch%i", ch_iter), Form("FitR2_vs_Edep_ch%i", ch_iter), hist_pars[0],hist_pars[1],hist_pars[2], 200,0,2 );
- th2_hist_ptr->SetTitle( Form("FitR2 vs EdepFPGA channel %i; Charge [mV]; R2 []", ch_iter) );
- }
-
- TString common_selection = Form("(Nhits > %i)", nhits);
- //draw selection
- for(int ch_iter = 0; ch_iter < total_channels; ch_iter++)
- {
- TString ich_selection = "";
- //ich_selection += Form("&&(channel_%i.Edep>1.5) ", ch_iter);//dont touch (not filled channels removal)
- switch(cut_option){
- case 0:
- break;
- case 1:
- ich_selection += Form("&&(channel_%i.TimeMax>=4 && channel_%i.TimeMax<=7)", ch_iter, ch_iter);
- break;
- case 2:
- ich_selection += Form("&&(channel_%i.TimeMax>=4 && channel_%i.TimeMax<=7)", ch_iter, ch_iter);
- ich_selection += Form("&&(channel_%i.FitR2<0.4)", ch_iter);
- break;
- }
-
- TString hist_name = Form("Edep_ch%i", ch_iter);
- data_tree->Draw( Form("channel_%i.%s>>%s", ch_iter, calc_mode.Data(), hist_name.Data()), (common_selection+ich_selection).Data(), "", events_to_read );
-
- hist_name = Form("timeMax_vs_Edep_ch%i", ch_iter);
- data_tree->Draw( Form("channel_%i.%s:channel_%i.TimeMax>>%s", ch_iter, calc_mode.Data(), ch_iter, hist_name.Data()),
- (common_selection+ich_selection).Data(), "", events_to_read );
-
- hist_name = Form("FitR2_vs_Edep_ch%i", ch_iter);
- data_tree->Draw( Form("channel_%i.FitR2:channel_%i.EdepFPGA>>%s", ch_iter, ch_iter, /*calc_mode.Data(),*/ hist_name.Data()),
- (common_selection+ich_selection).Data(), "", events_to_read );
-
- printf("selection for ch %i: [%s]\n", ch_iter, (common_selection+ich_selection).Data());
- }
-
-
- //drawing histogramms
- for(int ch_iter = 0; ch_iter < total_channels; ch_iter++)
- {
- double mip = 5.0; //MeV
- TH1* th1_ch_ptr = ((TH1*)(gDirectory->FindObjectAny( Form("Edep_ch%i", ch_iter) )));
-
- th1_hist_ptr = ((TH1*)(gDirectory->FindObjectAny( "EdepProfile" )));
- //th1_hist_ptr->SetBinContent(ch_iter+1, th1_ch_ptr->GetMean());//*th1_ch_ptr->GetEntries()/data_tree->GetEntries());
- th1_hist_ptr->SetBinContent(ch_iter+1, th1_ch_ptr->GetMean()*th1_ch_ptr->GetEntries()/temp->GetEntries());
-
- }
-
-
-
- TCanvas *canv_charge_Edep = new TCanvas("canv_charge_Edep", "canv_charge_Edep");
- canvas_array->Add(canv_charge_Edep);
- canv_charge_Edep->DivideSquare(total_channels);
- for (int ch_iter = 0; ch_iter < total_channels; ch_iter++){
- canv_charge_Edep->cd(ch_iter+1)->SetLogy(1);
- th1_hist_ptr = ((TH1*)(gDirectory->FindObjectAny( Form("Edep_ch%i", ch_iter) )));
- th1_hist_ptr->Draw();
- result_array->Add(th1_hist_ptr);
- int start = 8;
- int stop = 20;
- Int_t binx_min = th1_hist_ptr->GetXaxis()->FindBin(start);
- Int_t binx_max = th1_hist_ptr->GetXaxis()->FindBin(stop);
- double data_integral = th1_hist_ptr->Integral(binx_min, binx_max);
-
- th1_hist_ptr = ((TH1F*)(file_sim->FindObjectAny( Form("hf4_%i", ch_iter+1) )));
- if(th1_hist_ptr == nullptr) {cout<<"no hist"<<endl; break;}
- th1_hist_ptr->SetLineColor(kMagenta);
- th1_hist_ptr->SetLineWidth(2);
- binx_min = th1_hist_ptr->GetXaxis()->FindBin(start);
- binx_max = th1_hist_ptr->GetXaxis()->FindBin(stop);
- th1_hist_ptr->Scale(data_integral/th1_hist_ptr->Integral(binx_min, binx_max));
- th1_hist_ptr->Draw("same hist");
-
- }
-
-
- TCanvas *canv_time_Edep = new TCanvas("canv_time_Edep", "canv_time_Edep");
- canvas_array->Add(canv_time_Edep);
- canv_time_Edep->DivideSquare(total_channels);
- for (int ch_iter = 0; ch_iter < total_channels; ch_iter++){
- canv_time_Edep->cd(ch_iter+1)->SetLogz(1);
- th2_hist_ptr = ((TH2*)(gDirectory->FindObjectAny( Form("timeMax_vs_Edep_ch%i", ch_iter) )));
- th2_hist_ptr->Draw("colz");
- result_array->Add(th2_hist_ptr);
- }
-
- TCanvas *canv_R2_Edep = new TCanvas("canv_R2_Edep", "canv_R2_Edep");
- canvas_array->Add(canv_R2_Edep);
- canv_R2_Edep->DivideSquare(total_channels);
- for (int ch_iter = 0; ch_iter < total_channels; ch_iter++){
- canv_R2_Edep->cd(ch_iter+1)->SetLogz(1);
- th2_hist_ptr = ((TH2*)(gDirectory->FindObjectAny( Form("FitR2_vs_Edep_ch%i", ch_iter) )));
- th2_hist_ptr->Draw("colz");
- result_array->Add(th2_hist_ptr);
- }
-
- TCanvas *canv_Profile = new TCanvas("canv_Profile", "canv_Profile");
- canvas_array->Add(canv_Profile);
- th1_hist_ptr = ((TH1*)(gDirectory->FindObjectAny( "EdepProfile" )));
- th1_hist_ptr->GetYaxis()->SetRangeUser(0, 1.1*th1_hist_ptr->GetBinContent(1));
- th1_hist_ptr->Draw("hist");
-
- TFile* file_simm = TFile::Open("/home/nikolay/git_projects/mPSD_test_bench/data/ONi2AGeV_sim/profile_4.00MeV_trshd_10kk.root");
- TCanvas* canv_sim = (TCanvas*) file_simm->Get("c1");
- th1_hist_ptr = (TH1F*) canv_sim->GetPrimitive("hEnergyProfile");
- th1_hist_ptr->SetLineColor(kMagenta);
- th1_hist_ptr->SetLineWidth(2);
- canv_Profile->cd();
- th1_hist_ptr->Draw("same hist");
-
- // saving results #####################################################################
-
- for (Int_t i = 0; i < canvas_array->GetLast (); i++)
- ((TCanvas*) canvas_array->At(i))->SaveAs ((result_file_name + ".pdf(").Data ());
- ((TCanvas*) canvas_array->At (canvas_array->GetLast ()))->SaveAs ((result_file_name + ".pdf)").Data ());
- printf("file %s was written\n", (result_file_name + ".pdf").Data());
-
-
- result_array->Add(canvas_array);
- TFile *result_file = new TFile((result_file_name + ".root").Data(), "RECREATE");
- TString start_path = result_array->GetName();
- result_file->mkdir(start_path.Data());
- result_file->cd(start_path.Data());
- TIter nx_iter((TCollection*)(result_array));
- TObject* obj_ptr;
- while ( (obj_ptr=(TObjArray*)nx_iter()) )
- {
- if(obj_ptr == nullptr) continue;
- printf("Writing %s object\n", obj_ptr->GetName());
- obj_ptr->Write();
- }
- printf("file %s was written\n", (result_file_name + ".root").Data());
- result_file->Close();
- delete result_file;
-
-
- return 0;
+ for (int ch_iter = 0; ch_iter < total_channels; ch_iter++) {
+ th1_hist_ptr =
+ new TH1F(Form("Edep_ch%i", ch_iter), Form("Edep_ch%i", ch_iter), hist_pars[0], hist_pars[1], hist_pars[2]);
+ th1_hist_ptr->SetTitle(Form("Charge channel %i; Edep [MeV]", ch_iter));
+ th2_hist_ptr = new TH2F(Form("timeMax_vs_Edep_ch%i", ch_iter), Form("timeMax_vs_Edep_ch%i", ch_iter), 32, 0, 32,
+ hist_pars[0], hist_pars[1], hist_pars[2]);
+ th2_hist_ptr->SetTitle(Form("TimeMax vs Edep channel %i; TimeMax [12,5ns]; Charge [mV]", ch_iter));
+
+ th2_hist_ptr = new TH2F(Form("FitR2_vs_Edep_ch%i", ch_iter), Form("FitR2_vs_Edep_ch%i", ch_iter), hist_pars[0],
+ hist_pars[1], hist_pars[2], 200, 0, 2);
+ th2_hist_ptr->SetTitle(Form("FitR2 vs EdepFPGA channel %i; Charge [mV]; R2 []", ch_iter));
+ }
+
+ TString common_selection = Form("(Nhits > %i)", nhits);
+ //draw selection
+ for (int ch_iter = 0; ch_iter < total_channels; ch_iter++) {
+ TString ich_selection = "";
+ //ich_selection += Form("&&(channel_%i.Edep>1.5) ", ch_iter);//dont touch (not filled channels removal)
+ switch (cut_option) {
+ case 0: break;
+ case 1: ich_selection += Form("&&(channel_%i.TimeMax>=4 && channel_%i.TimeMax<=7)", ch_iter, ch_iter); break;
+ case 2:
+ ich_selection += Form("&&(channel_%i.TimeMax>=4 && channel_%i.TimeMax<=7)", ch_iter, ch_iter);
+ ich_selection += Form("&&(channel_%i.FitR2<0.4)", ch_iter);
+ break;
+ }
+
+ TString hist_name = Form("Edep_ch%i", ch_iter);
+ data_tree->Draw(Form("channel_%i.%s>>%s", ch_iter, calc_mode.Data(), hist_name.Data()),
+ (common_selection + ich_selection).Data(), "", events_to_read);
+
+ hist_name = Form("timeMax_vs_Edep_ch%i", ch_iter);
+ data_tree->Draw(Form("channel_%i.%s:channel_%i.TimeMax>>%s", ch_iter, calc_mode.Data(), ch_iter, hist_name.Data()),
+ (common_selection + ich_selection).Data(), "", events_to_read);
+
+ hist_name = Form("FitR2_vs_Edep_ch%i", ch_iter);
+ data_tree->Draw(
+ Form("channel_%i.FitR2:channel_%i.EdepFPGA>>%s", ch_iter, ch_iter, /*calc_mode.Data(),*/ hist_name.Data()),
+ (common_selection + ich_selection).Data(), "", events_to_read);
+
+ printf("selection for ch %i: [%s]\n", ch_iter, (common_selection + ich_selection).Data());
+ }
+
+
+ //drawing histogramms
+ for (int ch_iter = 0; ch_iter < total_channels; ch_iter++) {
+ double mip = 5.0; //MeV
+ TH1* th1_ch_ptr = ((TH1*) (gDirectory->FindObjectAny(Form("Edep_ch%i", ch_iter))));
+
+ th1_hist_ptr = ((TH1*) (gDirectory->FindObjectAny("EdepProfile")));
+ //th1_hist_ptr->SetBinContent(ch_iter+1, th1_ch_ptr->GetMean());//*th1_ch_ptr->GetEntries()/data_tree->GetEntries());
+ th1_hist_ptr->SetBinContent(ch_iter + 1, th1_ch_ptr->GetMean() * th1_ch_ptr->GetEntries() / temp->GetEntries());
+ }
+
+
+ TCanvas* canv_charge_Edep = new TCanvas("canv_charge_Edep", "canv_charge_Edep");
+ canvas_array->Add(canv_charge_Edep);
+ canv_charge_Edep->DivideSquare(total_channels);
+ for (int ch_iter = 0; ch_iter < total_channels; ch_iter++) {
+ canv_charge_Edep->cd(ch_iter + 1)->SetLogy(1);
+ th1_hist_ptr = ((TH1*) (gDirectory->FindObjectAny(Form("Edep_ch%i", ch_iter))));
+ th1_hist_ptr->Draw();
+ result_array->Add(th1_hist_ptr);
+ int start = 8;
+ int stop = 20;
+ Int_t binx_min = th1_hist_ptr->GetXaxis()->FindBin(start);
+ Int_t binx_max = th1_hist_ptr->GetXaxis()->FindBin(stop);
+ double data_integral = th1_hist_ptr->Integral(binx_min, binx_max);
+
+ th1_hist_ptr = ((TH1F*) (file_sim->FindObjectAny(Form("hf4_%i", ch_iter + 1))));
+ if (th1_hist_ptr == nullptr) {
+ cout << "no hist" << endl;
+ break;
+ }
+ th1_hist_ptr->SetLineColor(kMagenta);
+ th1_hist_ptr->SetLineWidth(2);
+ binx_min = th1_hist_ptr->GetXaxis()->FindBin(start);
+ binx_max = th1_hist_ptr->GetXaxis()->FindBin(stop);
+ th1_hist_ptr->Scale(data_integral / th1_hist_ptr->Integral(binx_min, binx_max));
+ th1_hist_ptr->Draw("same hist");
+ }
+
+
+ TCanvas* canv_time_Edep = new TCanvas("canv_time_Edep", "canv_time_Edep");
+ canvas_array->Add(canv_time_Edep);
+ canv_time_Edep->DivideSquare(total_channels);
+ for (int ch_iter = 0; ch_iter < total_channels; ch_iter++) {
+ canv_time_Edep->cd(ch_iter + 1)->SetLogz(1);
+ th2_hist_ptr = ((TH2*) (gDirectory->FindObjectAny(Form("timeMax_vs_Edep_ch%i", ch_iter))));
+ th2_hist_ptr->Draw("colz");
+ result_array->Add(th2_hist_ptr);
+ }
+
+ TCanvas* canv_R2_Edep = new TCanvas("canv_R2_Edep", "canv_R2_Edep");
+ canvas_array->Add(canv_R2_Edep);
+ canv_R2_Edep->DivideSquare(total_channels);
+ for (int ch_iter = 0; ch_iter < total_channels; ch_iter++) {
+ canv_R2_Edep->cd(ch_iter + 1)->SetLogz(1);
+ th2_hist_ptr = ((TH2*) (gDirectory->FindObjectAny(Form("FitR2_vs_Edep_ch%i", ch_iter))));
+ th2_hist_ptr->Draw("colz");
+ result_array->Add(th2_hist_ptr);
+ }
+
+ TCanvas* canv_Profile = new TCanvas("canv_Profile", "canv_Profile");
+ canvas_array->Add(canv_Profile);
+ th1_hist_ptr = ((TH1*) (gDirectory->FindObjectAny("EdepProfile")));
+ th1_hist_ptr->GetYaxis()->SetRangeUser(0, 1.1 * th1_hist_ptr->GetBinContent(1));
+ th1_hist_ptr->Draw("hist");
+
+ TFile* file_simm =
+ TFile::Open("/home/nikolay/git_projects/mPSD_test_bench/data/ONi2AGeV_sim/profile_4.00MeV_trshd_10kk.root");
+ TCanvas* canv_sim = (TCanvas*) file_simm->Get("c1");
+ th1_hist_ptr = (TH1F*) canv_sim->GetPrimitive("hEnergyProfile");
+ th1_hist_ptr->SetLineColor(kMagenta);
+ th1_hist_ptr->SetLineWidth(2);
+ canv_Profile->cd();
+ th1_hist_ptr->Draw("same hist");
+
+ // saving results #####################################################################
+
+ for (Int_t i = 0; i < canvas_array->GetLast(); i++)
+ ((TCanvas*) canvas_array->At(i))->SaveAs((result_file_name + ".pdf(").Data());
+ ((TCanvas*) canvas_array->At(canvas_array->GetLast()))->SaveAs((result_file_name + ".pdf)").Data());
+ printf("file %s was written\n", (result_file_name + ".pdf").Data());
+
+
+ result_array->Add(canvas_array);
+ TFile* result_file = new TFile((result_file_name + ".root").Data(), "RECREATE");
+ TString start_path = result_array->GetName();
+ result_file->mkdir(start_path.Data());
+ result_file->cd(start_path.Data());
+ TIter nx_iter((TCollection*) (result_array));
+ TObject* obj_ptr;
+ while ((obj_ptr = (TObjArray*) nx_iter())) {
+ if (obj_ptr == nullptr) continue;
+ printf("Writing %s object\n", obj_ptr->GetName());
+ obj_ptr->Write();
+ }
+ printf("file %s was written\n", (result_file_name + ".root").Data());
+ result_file->Close();
+ delete result_file;
+
+
+ return 0;
}
#ifdef __CINT__
-
+
#pragma link off all globals;
#pragma link off all classes;
#pragma link off all functions;
#pragma link C++ nestedclasses;
-
-#pragma link C++ struct data_header_struct+;
-#pragma link C++ struct digi_struct+;
-#endif
+#pragma link C++ struct data_header_struct + ;
+#pragma link C++ struct digi_struct + ;
+#endif
diff --git a/macro/psd/event_data_struct.h b/macro/psd/event_data_struct.h
index 0cc92e485c348aef0c638c29a4cc9dc50ceb2c2f..cb1e28cb8cb570d3007cc48e0eedd1be9e4e933a 100755
--- a/macro/psd/event_data_struct.h
+++ b/macro/psd/event_data_struct.h
@@ -1,85 +1,80 @@
-#ifndef EVENT_DATA_STRUCT
-#define EVENT_DATA_STRUCT
-
-#include<TTree.h>
-#include<TString.h>
-#include<TMath.h>
-#include<TObject.h>
-#include<fstream>
-
-
-//data processed from waveform
-struct event_data_struct
-{
-
- Float_t EdepFPGA;
- Float_t EdepWfm;
- Float_t Ampl;
- Float_t Minimum;
- Float_t ZL;
- Float_t Time;
- Int_t TimeMax;
- Float_t TimeHalf;
- Float_t TimeInEvent;
-
- Float_t FitEdep;
- Float_t FitAmpl;
- Float_t FitZL;
- Float_t FitR2;
- Float_t FitTimeMax;
-
- ULong64_t EventTime;
-
-
- void reset()
- {
- EdepFPGA = 0.;
- EdepWfm = 0.;
- Ampl = 0.;
- Minimum = 0.;
- ZL = 0.;
- Time = 0.;
- TimeMax = -1;
- TimeHalf = -1.;
- TimeInEvent = -1.;
-
- FitEdep = 0.;
- FitAmpl = 0.;
- FitZL = 0.;
- FitR2 = 999.;
- FitTimeMax = -1.;
-
- EventTime = 0;
- }
-
- void Print()
- {
- printf("[EdepFPGA %.0f; EdepWfm %.0f; Ampl %.0f; Minumum %.0f; ZL %.2f; Time %.2f; TimeMax %i; TimeHalf %.2f; TimeInEvent %.2f; FitEdep %.2f; FitAmpl %.2f; FitZL %.2f; FitR2 %.2f; FitTimeMax %.2f; EventTime %llu]",
- EdepFPGA, EdepWfm, Ampl, Minimum, ZL, Time, TimeMax, TimeHalf, TimeInEvent, FitEdep, FitAmpl, FitZL, FitR2, FitTimeMax, EventTime);
- }
-
- static TString GetChName(Int_t channel_num)
- {
- return TString::Format("channel_%i", channel_num);
- }
-
- TBranch* CreateBranch(TTree *tree, Int_t channel_num)
- {
- return tree->Branch(GetChName(channel_num).Data(), this,
- "EdepFPGA/F:EdepWfm/F:Ampl/F:Minimum/F:ZL/F:Time/F:TimeMax/I:TimeHalf/F:TimeInEvent/F:FitEdep/F:FitAmpl/F:FitZL/F:FitR2/F:FitTimeMax/F:EventTime/l");
- }
-
-
- Int_t SetBranch(TTree *tree, Int_t channel_num)
- {
- return tree->SetBranchAddress(GetChName(channel_num).Data(), this);
- }
-
-};
-
-
-
-
-
-#endif // EVENT_DATA_STRUCT
-
+#ifndef EVENT_DATA_STRUCT
+#define EVENT_DATA_STRUCT
+
+#include <TMath.h>
+#include <TObject.h>
+#include <TString.h>
+#include <TTree.h>
+
+#include <fstream>
+
+
+//data processed from waveform
+struct event_data_struct {
+
+ Float_t EdepFPGA;
+ Float_t EdepWfm;
+ Float_t Ampl;
+ Float_t Minimum;
+ Float_t ZL;
+ Float_t Time;
+ Int_t TimeMax;
+ Float_t TimeHalf;
+ Float_t TimeInEvent;
+
+ Float_t FitEdep;
+ Float_t FitAmpl;
+ Float_t FitZL;
+ Float_t FitR2;
+ Float_t FitTimeMax;
+
+ ULong64_t EventTime;
+
+
+ void reset()
+ {
+ EdepFPGA = 0.;
+ EdepWfm = 0.;
+ Ampl = 0.;
+ Minimum = 0.;
+ ZL = 0.;
+ Time = 0.;
+ TimeMax = -1;
+ TimeHalf = -1.;
+ TimeInEvent = -1.;
+
+ FitEdep = 0.;
+ FitAmpl = 0.;
+ FitZL = 0.;
+ FitR2 = 999.;
+ FitTimeMax = -1.;
+
+ EventTime = 0;
+ }
+
+ void Print()
+ {
+ printf("[EdepFPGA %.0f; EdepWfm %.0f; Ampl %.0f; Minumum %.0f; ZL %.2f; Time %.2f; TimeMax %i; TimeHalf %.2f; "
+ "TimeInEvent %.2f; FitEdep %.2f; FitAmpl %.2f; FitZL %.2f; FitR2 %.2f; FitTimeMax %.2f; EventTime %llu]",
+ EdepFPGA, EdepWfm, Ampl, Minimum, ZL, Time, TimeMax, TimeHalf, TimeInEvent, FitEdep, FitAmpl, FitZL, FitR2,
+ FitTimeMax, EventTime);
+ }
+
+ static TString GetChName(Int_t channel_num) { return TString::Format("channel_%i", channel_num); }
+
+ TBranch* CreateBranch(TTree* tree, Int_t channel_num)
+ {
+ return tree->Branch(GetChName(channel_num).Data(), this,
+ "EdepFPGA/F:EdepWfm/F:Ampl/F:Minimum/F:ZL/F:Time/F:TimeMax/I:TimeHalf/F:TimeInEvent/F:FitEdep/"
+ "F:FitAmpl/F:FitZL/F:FitR2/F:FitTimeMax/F:EventTime/l");
+ }
+
+
+ Int_t SetBranch(TTree* tree, Int_t channel_num)
+ {
+ return tree->SetBranchAddress(GetChName(channel_num).Data(), this);
+ }
+};
+
+
+#endif // EVENT_DATA_STRUCT