Added automatic DE detection project

Projects/AlPhaPha/2024/macro/chio/XYCalibration/CutAutoDEE.C

+#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;
+}

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();
+}