From 8a2d754c4bda48570ba5c832e60c51a0141694e2 Mon Sep 17 00:00:00 2001
From: "theodore.efremov" <theodore.efremov@cea.fr>
Date: Thu, 19 Dec 2024 15:53:36 +0100
Subject: [PATCH] Added automatic DE detection project
---
.../macro/chio/XYCalibration/CutAutoDEE.C | 110 ++++++++++++++++++
.../2024/macro/chio/XYCalibration/DECorr.C | 33 +++++-
2 files changed, 142 insertions(+), 1 deletion(-)
create mode 100644 Projects/AlPhaPha/2024/macro/chio/XYCalibration/CutAutoDEE.C
diff --git a/Projects/AlPhaPha/2024/macro/chio/XYCalibration/CutAutoDEE.C b/Projects/AlPhaPha/2024/macro/chio/XYCalibration/CutAutoDEE.C
new file mode 100644
index 000000000..9c7f11b3a
--- /dev/null
+++ b/Projects/AlPhaPha/2024/macro/chio/XYCalibration/CutAutoDEE.C
@@ -0,0 +1,110 @@
+#include <TCanvas.h>
+#include <TFile.h>
+#include <TGraph.h>
+#include <TH2.h>
+#include <TSpectrum.h>
+#include <vector>
+
+using namespace std;
+
+void PeakFinder(vector<TSpectrum*> *sPeak , vector<int> *vPeak, TH2F *h , vector<vector<Double_t>> *LinePos);
+vector<vector<vector<Double_t>>> PeakLinker(vector<vector<Double_t>> *LinePos);
+
+void CutAutoDEE(){
+
+ //retrieve initial histogram
+ TFile *inFile = new TFile("Output/DE_E.root","open");
+ TH2F *hDE_E = (TH2F*)inFile->Get("DE0234_Ecorr");
+
+ //inverse data to use tspectrum
+ for (int xBin=1; xBin<= hDE_E->GetNbinsX();xBin++){
+ for (int yBin=1; yBin<= hDE_E->GetNbinsY();yBin++){
+ hDE_E->SetBinContent(xBin,yBin , -hDE_E->GetBinContent(xBin,yBin));
+ }
+ }
+ //Store result of peak finding
+ vector<int> *vPeaks = new vector<int>() ;
+ vector<TSpectrum*> *vPeakFinder = new vector<TSpectrum*>();
+ vector<vector<Double_t>> *LinePos = new vector<vector<Double_t>> ;
+
+ PeakFinder(vPeakFinder,vPeaks,hDE_E,LinePos);
+ vector<vector<vector<Double_t>>> CutLine = PeakLinker(LinePos);
+}
+
+void PeakFinder(vector<TSpectrum*> *sPeak , vector<int> *vPeak, TH2F *h, vector<vector<Double_t>> *LinePos){
+
+ int nbinsX = h->GetNbinsX();
+ int step = 10;
+ int iter = 0;
+
+
+ for (int i = 1; i <= nbinsX; i += step) {
+ iter +=1;
+ int x1 = i;
+ int x2 = std::min(i + step - 1, nbinsX);
+
+ TString projName = Form("projY_%d_to_%d", x1, x2);
+ TH1D* projY = h->ProjectionY(projName, x1, x2);
+
+ projY->SetTitle(Form("Projection Y: bins %d to %d", x1, x2));
+
+ TSpectrum *PeakFinder = new TSpectrum(55);
+ int nPeaks = PeakFinder->Search(projY,2,"",0.02);
+
+ vPeak->push_back(nPeaks);
+ sPeak->push_back(PeakFinder);
+
+ if (nPeaks <= 0) {
+ cout << "No peaks found!" << endl;
+ return;
+ }
+
+ //Save the position of each peak
+ Double_t * x = PeakFinder->GetPositionX();
+
+ for(int peaks = 0 ; peaks < nPeaks ; peaks++){
+ vector<Double_t> Position(2);
+ Position.at(0) = (double(i) - double(step)/2.0) * h->GetXaxis()->GetBinWidth(4) ; //Retrieve E position
+ Position.at(1) = *PeakFinder->GetPositionX();
+ LinePos->push_back(Position);
+ }
+
+
+ //Test
+ if (true && iter ==58){
+ // Draw the histogram
+ TCanvas* c = new TCanvas(Form("c%d",i), Form("Histogram with Peaks%d",i), 800, 600);
+ projY->Draw();
+
+
+ }
+ }
+}
+vector<vector<vector<Double_t>>> PeakLinker(vector<vector<Double_t>> *LinePos){
+
+ vector<vector<vector<Double_t>>> Res ;
+ vector<Double_t> X , Y ;
+
+ //Retrieve position of peaks
+ for (int ipeak; ipeak<LinePos->size() ; ipeak ++){
+ X.push_back(LinePos->at(ipeak).at(0));
+ Y.push_back(LinePos->at(ipeak).at(1));
+ }
+
+ // Start from the lowest position in Y
+ // Corresponding to the lowest line in DE-E
+ //
+ auto minElementY = min_element(Y.begin(),Y.end());
+ if (minElementY != Y.end())
+ {
+ int indexMin = distance(Y.begin(), minElementY);
+ //cout << Y[indexMin] << endl;
+
+ //Now we need to find the good threshold between each line in the DE-E
+ //To do so we will get the peaks from the same X and get their distance in Y
+ //This will be our threshold
+
+ }
+
+ return Res;
+}
diff --git a/Projects/AlPhaPha/2024/macro/chio/XYCalibration/DECorr.C b/Projects/AlPhaPha/2024/macro/chio/XYCalibration/DECorr.C
index b5bc34a45..cfdf531e9 100644
--- a/Projects/AlPhaPha/2024/macro/chio/XYCalibration/DECorr.C
+++ b/Projects/AlPhaPha/2024/macro/chio/XYCalibration/DECorr.C
@@ -16,6 +16,7 @@
using namespace std;
vector<double> TxtToVector(const char *Path);
int GetNumberKey(TFile *infile ,const char* ClassName);
+void TH2toCSV(TH2F* h2, const char *PathOut) ;
TSpline3* MakeSpline(TH2F* hInput, Int_t NCallee, double XMin , double
XMax, double YMin, double YMax);
@@ -204,7 +205,9 @@ void DECorr(bool cut = false, bool spline = false) {
auto start = std::chrono::high_resolution_clock::now();
int NSegment = 11;
- for (int e = 0; e < Nentries; e++) {
+ //save data for clustering
+
+ for (int e = 0; e < Nentries; e++) {
if (e % 100000 == 0 && e > 0 ) {
auto now = std::chrono::high_resolution_clock::now();
std::chrono::duration<double> elapsed = now - start;
@@ -541,6 +544,17 @@ void DECorr(bool cut = false, bool spline = false) {
hDE_Ebis->Draw("colz");
hDE_Ebis->GetXaxis()->SetTitle("E");
hDE_Ebis->GetYaxis()->SetTitle("DE");
+
+
+//===========================================================================================================
+// Save for further analysis
+//===========================================================================================================
+
+TFile *oSaved = new TFile("Output/DE_E.root","recreate");
+hDE0234_E_corr->GetXaxis()->SetRangeUser(2000,34000);
+hDE0234_E_corr->Write();
+TH2toCSV(hDE0234_E_corr,"Output/DE_E.csv");
+
} // End spline chio XY
@@ -661,3 +675,20 @@ vector<double> TxtToVector(const char *Path){
return values;
}
+
+// Example: Export 2D data for external clustering
+void TH2toCSV(TH2F* h2, const char *PathOut) {
+ std::ofstream outFile(PathOut);
+ outFile << "E,DE,content\n"; // CSV header
+ for (int i = 1; i <= h2->GetNbinsX(); ++i) {
+ for (int j = 1; j <= h2->GetNbinsY(); ++j) {
+ double content = h2->GetBinContent(i, j);
+ if (content > 0) {
+ double x = h2->GetXaxis()->GetBinCenter(i);
+ double y = h2->GetYaxis()->GetBinCenter(j);
+ outFile << x << "," << y << "," << content << std::endl;
+ }
+ }
+ }
+ outFile.close();
+}
--
GitLab