From 4cf0e572d074e88cf85b15cf0effe410654a2882 Mon Sep 17 00:00:00 2001
From: Nikolay Karpushkin <karpushkin@inr.ru>
Date: Sat, 10 Jul 2021 18:54:24 +0200
Subject: [PATCH] psd macros

---
 macro/psd/analyse_events.sh   |  27 ++++
 macro/psd/build_events.cpp    | 240 +++++++++++++++++++++++++++++++
 macro/psd/compare_to_sim.cpp  | 260 ++++++++++++++++++++++++++++++++++
 macro/psd/event_data_struct.h |  85 +++++++++++
 macro/psd/make_results.sh     |  19 +++
 5 files changed, 631 insertions(+)
 create mode 100755 macro/psd/analyse_events.sh
 create mode 100644 macro/psd/build_events.cpp
 create mode 100644 macro/psd/compare_to_sim.cpp
 create mode 100755 macro/psd/event_data_struct.h
 create mode 100755 macro/psd/make_results.sh

diff --git a/macro/psd/analyse_events.sh b/macro/psd/analyse_events.sh
new file mode 100755
index 0000000000..c54a612cad
--- /dev/null
+++ b/macro/psd/analyse_events.sh
@@ -0,0 +1,27 @@
+#!/bin/bash
+
+export FILELIST=$1
+export FILEINDX=$2
+export MAXTSA=$3
+export MAXEVENT=$4
+
+while read F  ; do
+	echo "making events from ${F}/unp_mcbm_00$FILEINDX.root"
+    root -l -b -q "PsdBuildEvents.cpp(\"${F}/unp_mcbm_00$FILEINDX.root\", \"${F}/EventTree$FILEINDX.root\", $MAXTSA)"
+
+	echo "processing ${F}/EventTree$FILEINDX.root"
+	mkdir -p ${F}/results
+	root -l -b -q "compare_to_sim.cpp(\"${F}/EventTree$FILEINDX.root\", \"${F}/results/Results\", \"EdepFPGA\", $MAXEVENT, 0, 0)"
+	root -l -b -q "compare_to_sim.cpp(\"${F}/EventTree$FILEINDX.root\", \"${F}/results/Results\", \"EdepFPGA\", $MAXEVENT, 0, 1)"
+	root -l -b -q "compare_to_sim.cpp(\"${F}/EventTree$FILEINDX.root\", \"${F}/results/Results\", \"EdepFPGA\", $MAXEVENT, 0, 2)"
+
+	root -l -b -q "compare_to_sim.cpp(\"${F}/EventTree$FILEINDX.root\", \"${F}/results/Results\", \"EdepWfm\", $MAXEVENT, 0, 0)"
+	root -l -b -q "compare_to_sim.cpp(\"${F}/EventTree$FILEINDX.root\", \"${F}/results/Results\", \"EdepWfm\", $MAXEVENT, 0, 1)"
+	root -l -b -q "compare_to_sim.cpp(\"${F}/EventTree$FILEINDX.root\", \"${F}/results/Results\", \"EdepWfm\", $MAXEVENT, 0, 2)"
+
+	root -l -b -q "compare_to_sim.cpp(\"${F}/EventTree$FILEINDX.root\", \"${F}/results/Results\", \"FitEdep\", $MAXEVENT, 0, 0)"
+	root -l -b -q "compare_to_sim.cpp(\"${F}/EventTree$FILEINDX.root\", \"${F}/results/Results\", \"FitEdep\", $MAXEVENT, 0, 1)"
+	root -l -b -q "compare_to_sim.cpp(\"${F}/EventTree$FILEINDX.root\", \"${F}/results/Results\", \"FitEdep\", $MAXEVENT, 0, 2)"
+done <$FILELIST
+
+
diff --git a/macro/psd/build_events.cpp b/macro/psd/build_events.cpp
new file mode 100644
index 0000000000..e640531619
--- /dev/null
+++ b/macro/psd/build_events.cpp
@@ -0,0 +1,240 @@
+#include <string.h>
+#include <math.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 "event_data_struct.h"
+
+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 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);
+		hit_time = 0.;
+
+ 		CbmPsdDsp dsp;
+		for( int entr_iter = 0; entr_iter < dsps_for_time_btw; entr_iter++ )
+		{
+			inTree->GetEntry(entr_iter);
+			Int_t nPsdDsps = fPsdDspVect->size();
+			for (Int_t iDsp = 0; iDsp < nPsdDsps; iDsp++) {
+				dsp = (*fPsdDspVect)[iDsp];
+
+				hit_time = dsp.fdTime;
+				int hit_channel = dsp.GetSectionID();
+				if(hit_channel>9) continue;
+				h1_time_btw_hits->Fill(hit_time-prev_hit_time);
+				prev_hit_time = hit_time;
+			}
+
+		}
+
+		TCanvas *cc = new TCanvas();
+		gPad->SetLogy();
+		h1_time_btw_hits->Draw();
+
+		hit_time = 0.;
+		prev_hit_time = 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
new file mode 100644
index 0000000000..ff050c703b
--- /dev/null
+++ b/macro/psd/compare_to_sim.cpp
@@ -0,0 +1,260 @@
+#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 <TChain.h>
+
+#include <TGraphErrors.h>
+#include <TF1.h>
+#include <TStyle.h>
+
+#include "event_data_struct.h"
+
+
+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) {
+
+	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;
+
+}
+
+#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
+
diff --git a/macro/psd/event_data_struct.h b/macro/psd/event_data_struct.h
new file mode 100755
index 0000000000..0cc92e485c
--- /dev/null
+++ b/macro/psd/event_data_struct.h
@@ -0,0 +1,85 @@
+#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
+
diff --git a/macro/psd/make_results.sh b/macro/psd/make_results.sh
new file mode 100755
index 0000000000..a3447f1a52
--- /dev/null
+++ b/macro/psd/make_results.sh
@@ -0,0 +1,19 @@
+#!/bin/bash
+
+while read F  ; do
+	echo "processing ${F}/EventTree.root"
+	mkdir -p ${F}/results
+	root -l -b -q "compare_to_sim.cpp(\"${F}/EventTree.root\", \"${F}/results/Results\", \"EdepFPGA\", -1, 0, 0)"
+	root -l -b -q "compare_to_sim.cpp(\"${F}/EventTree.root\", \"${F}/results/Results\", \"EdepFPGA\", -1, 0, 1)"
+	root -l -b -q "compare_to_sim.cpp(\"${F}/EventTree.root\", \"${F}/results/Results\", \"EdepFPGA\", -1, 0, 2)"
+
+	root -l -b -q "compare_to_sim.cpp(\"${F}/EventTree.root\", \"${F}/results/Results\", \"EdepWfm\", -1, 0, 0)"
+	root -l -b -q "compare_to_sim.cpp(\"${F}/EventTree.root\", \"${F}/results/Results\", \"EdepWfm\", -1, 0, 1)"
+	root -l -b -q "compare_to_sim.cpp(\"${F}/EventTree.root\", \"${F}/results/Results\", \"EdepWfm\", -1, 0, 2)"
+
+	root -l -b -q "compare_to_sim.cpp(\"${F}/EventTree.root\", \"${F}/results/Results\", \"FitEdep\", -1, 0, 0)"
+	root -l -b -q "compare_to_sim.cpp(\"${F}/EventTree.root\", \"${F}/results/Results\", \"FitEdep\", -1, 0, 1)"
+	root -l -b -q "compare_to_sim.cpp(\"${F}/EventTree.root\", \"${F}/results/Results\", \"FitEdep\", -1, 0, 2)"
+done <filelist.txt
+
+
-- 
GitLab