CutAutoDEE_FineTuning.C

#include "CutAutoDEE_FineTuning.h"
#include <numeric>
#include <set>
#include <iostream>
#include <unordered_map>

// Important variable
#define STEP  2;        // Step to parse through the histogram in X
#define BINSIZE 10;    // Size of the bin in X on which to projectY the plot
#define THRESHOLDY 75.0;
#define THRESHOLDX 200.0;

//Main
void CutAutoDEE_FineTuning(){

    // Retrieve initial histogram
    TFile *inFile = new TFile("Output/DE_E_merged.root","read");
    TH2F *hDE_E = (TH2F*)inFile->Get("DE_E");
    // 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
    //This is not a good way to do this : We should normalize the structure of the code by
    //directly using a PeakData class to store each of the peak.
    vector<int> *vPeaks = new vector<int>() ;
    vector<TSpectrum*> *vPeakFinder = new vector<TSpectrum*>();
    vector<vector<Double_t>> *LinePos = new vector<vector<Double_t>> ;

    // Tspectrum used to find the peak : If you have problem with peak finding
    // adjust the setting in this function 
    // Step and Binsize are used in this function
    Double_t sigma = 3;
    Double_t threshold = 0.0001;
    PeakFinder(vPeakFinder,vPeaks,hDE_E,LinePos,sigma,threshold);


    // Link the peak, this is where the most critical variable are.
    // Specifically the maximum distance between each peak and the minimum number of peak in a line 
    vector<TCutG*> CutLine =  PeakLinker(LinePos);

    // Inverse data for plotting
    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));
        }
    }
    
    //Plot
    gStyle->SetPalette(kRainBow);
    TCanvas* c2 = new TCanvas("c2","c2");
    hDE_E->Draw("colz");

    TCanvas* c3 = new TCanvas("c3","c3");
    hDE_E->Draw("colz");
    for (auto cut : CutLine){
        cut->Draw("SAME");
    }
    
    //Save the cut !
    TFile *ofile = new TFile("Output/CutAutoZtest.root","RECREATE");
    for (auto cut : CutLine){
        cut->Write();
    }
    ofile->Close();
}

/** 
  *  Find the peak of an histogram by subdivinding the histogram on it's x axis.
  *  Store the number of peak and their position in different vector
  *
  *  input and Output:
  *     vector<TSpectrum*> *sPeak : Vector of TSpectrum to output
  *     vector<int> *vPeak        : Vector of number of peak for each bin range
  *     vector<vector<Double_t>> *LinePos : Vector of position of each peak for each bin range 
  */
void PeakFinder(vector<TSpectrum*> *sPeak , vector<int> *vPeak, TH2F *h, vector<vector<Double_t>> *LinePos, Double_t sigma, Double_t threshold){

    // Initialising
    int nbinsX = h->GetNbinsX();
    int iter = 0;
    int stepParsing = STEP;
    int SizeBin = BINSIZE;

    // Looping on the histogram x axis.
    for (int i = 1; i <= nbinsX; i += stepParsing) {

        // Limit in binnumber
        iter +=1;
        int x1 = i;
        int x2 = std::min(i + SizeBin - 1, nbinsX);

        // Projection on Y between the binnumber
        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));

        // PeakFinder
		TSpectrum *PeakFinder = new TSpectrum(600); 
		int nPeaks = PeakFinder->Search(projY,sigma,"",threshold);
      
        //Store npeaks and tspectrum
        vPeak->push_back(nPeaks); // number of peak for this bin range
        sPeak->push_back(PeakFinder);

        //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(SizeBin)/2.0) * h->GetXaxis()->GetBinWidth(4) ; //Retrieve E position as the center
                                                                                              //of the bin
            Position.at(1) = PeakFinder->GetPositionX()[peaks]; // Position in DE
            LinePos->push_back(Position);
        }


        //Test
        if (true && iter ==50){
            // Draw the histogram
            TCanvas* c = new TCanvas(Form("c%d",i), Form("Histogram with E = %f",(h->GetXaxis()->GetBinCenter(i))), 800, 600);
            projY->Draw();
        }

    }
}

/** 
  *  Link the different point to create TCUT from the peak found previously
  *  To do so a max range is given and then all point separated by a distance inferior to it are linked.
  *
  *  Input :
  *     vec<vec<Double>> *LinePos : pointer to a vector of length number of peak holding position in x and y for each peak
  *  Output:
  *    TCUTG : a cut for each charge in vamos
  */
vector<TCutG*> PeakLinker(vector<vector<Double_t>> *LinePos){

    //Retrieve position of peaks
    vector<Double_t> X , Y ;
    for (int ipeak = 0; ipeak<LinePos->size() ; ipeak ++){
        X.push_back(LinePos->at(ipeak).at(0));  
        Y.push_back(LinePos->at(ipeak).at(1));
    }

    //===========================================================================================================
    //                             1. Find the threshold between two point
    //===========================================================================================================
  
    // Find the distance between two point in X.
    //Sort the vec
    vector<Double_t> XSub(X); 
    sort(XSub.begin(),XSub.end());    

	//Take only the unique element 
    set<Double_t> SetX(XSub.begin() , XSub.end());
    vector<Double_t> XThresh ;
    for (Double_t elem : SetX) {
        XThresh.push_back(elem);
    } 

	//Find the average dist in X
    vector<Double_t> DiffX = diff(XThresh);
    Double_t ThresholdXestim ;
    Double_t ThresholdX ;
	for (int i =0 ; i<DiffX.size() ; i++) { ThresholdXestim += DiffX.at(i);}
	ThresholdXestim = ThresholdXestim / DiffX.size();	


    // In Y the threshold is not easily findable, However it's critical that the threshold 
    // is superior to the average distance in X. If not line can't be found by construction
    // The threshold should also be inferior to the distance between each line in Y
    // To do so adjust the bin STEP at the beginning of the code.

    //Very important set
    ThresholdX = THRESHOLDX;
    Double_t ThresholdY = THRESHOLDY;
    cout << "Threshold X : "  << ThresholdX  << endl;
    cout << "Threshold X approx : "  << ThresholdXestim  << endl;
    cout <<   "N of X "       << SetX.size() << endl;
    cout << "ThresholdY " << ThresholdY << endl;

    //===========================================================================================================
    //                                           2.LINK POINT
    //===========================================================================================================
   
    // Now that the threshold is found we will link the point by comparing each point 
    // Notice that if two point have the same X they can't be linked 

    // Create a vector holding the id of each peak
    int ID = 0 ; // Initial Id of the first line
    vector<int> vID(LinePos->size(), -1000);

    // We want to loop on the UNIDENTIFIED peak.
    for (int ipeak = 0 ; ipeak < LinePos->size() ; ipeak ++){

        //Here do ID manipulation
        //Create a set of ID without the unidentified peak
        set<int> SetId(vID.begin() , vID.end());
        if (SetId.size() > 1 ) {SetId.erase(-1000) ;  ID = *prev(SetId.end()) + 1 ;}

        //If the current peak is identified set newID to it's ID
        if (vID.at(ipeak) != -1000) ID = vID.at(ipeak); // If the point is identified propagate it!

        Double_t XPosPeak = LinePos->at(ipeak).at(0);
        Double_t YPosPeak = LinePos->at(ipeak).at(1);


        // Parse through all the next peak
        // This could be done faster by sorting the list by X and Y before and stoping after a threshold
        // For now it'll do just fine
        for (int ipeakCompare = ipeak + 1  ; ipeakCompare < LinePos->size() ; ipeakCompare ++){
           
            if ((vID.at(ipeakCompare) != -1000) || (XPosPeak == LinePos->at(ipeakCompare).at(0)))  continue; // If the point is identified skip it !
            
            //Calculate the distance between the two peak
            Double_t XPosPeakComp = LinePos->at(ipeakCompare).at(0);
            Double_t YPosPeakComp = LinePos->at(ipeakCompare).at(1);
            Double_t dist = sqrt( pow(abs(XPosPeak-XPosPeakComp),2) +  pow(abs(YPosPeak-YPosPeakComp),2) );
            Double_t distX = sqrt( pow(abs(XPosPeak-XPosPeakComp),2) );
            Double_t distY = sqrt(  pow(abs(YPosPeak-YPosPeakComp),2) );

            //If the distance is inferior to the threshold then set both ID to the new one !
            if ((distX < ThresholdX )&& (distY <ThresholdY )) {
                vID.at(ipeak) = ID ;
                vID.at(ipeakCompare) = ID ;                    
            } // end threshold if 
        } // end compare
    }   // end loop peak

    //Check if needed
    /*
       for (int i=0 ; i<vID.size() ; i++){
       if (vID.at(i)!=-1000){
       cout << vID.at(i) << " , ";
       }
       }
       */




    // Here we flush all the id that have less element than a threshold value
    //===========================================================================================================
    //                                 3. Flush line with not enough value
    //===========================================================================================================
   
    //This is a pretty sensitive variable 
    int MinElem = SetX.size()/4;
    cout << "Min Elem " << MinElem << endl ;


    set<int> SetIdTemp(vID.begin() , vID.end());
    for (int Id : SetIdTemp){
        int NId = 0 ;
        for (int ipeak = 0 ; ipeak < vID.size() ; ipeak ++) {
            if (vID.at(ipeak) == Id) NId += 1 ;
        }
        if ( NId < MinElem ) {

            for (int ipeak = 0 ; ipeak < vID.size() ; ipeak ++) {
                if (vID.at(ipeak) == Id) vID.at(ipeak) = -1000 ;
            }

        }
    }

    //===========================================================================================================
    //                                    4. TCutG Creator
    //===========================================================================================================

    // We now have identified each peak to a id that correspond to a line in the DE E
    // Now we need to link the ADJACENT line and then transform those point into tcutg

    set<int> SetId(vID.begin() , vID.end()); // This is all the unique element of vID 
    vector<int> IndexPeak(LinePos->size());
    for (int ipeak = 0 ; ipeak < LinePos->size() ; ipeak ++){
        IndexPeak.at(ipeak) = ipeak ;
    }

    //===========================================================================================================
    //                      4.1 Sort by id
    //===========================================================================================================
    // The next part of the code SHOULD sort by meanY but does not.
    // It therefore only serve as a convoluted way to sort by ID.
    // We need to sort by id either way so no need to rewrite i
    
    vector<Double_t> MeanYId(SetId.size()) ;
    vector<Double_t> NelemYId(SetId.size()) ;

    //Loop on the element of the set
    for (int Id = 0 ; Id < SetId.size() ; Id ++){
        Double_t Sum = 0 , Npeak = 1;
        auto itSet = SetId.begin();
        advance(itSet,Id);

        //Loop on all the peak
        for (int ipeak = 0 ; ipeak < LinePos->size() ; ipeak ++){

            // If the ID of the peak correspond to one of the set then increment his Y and NPeak value
            if (vID.at(ipeak) == *itSet )
            {
                Npeak += 1 ; 
                Sum += LinePos->at(ipeak).at(1);
            } // End if Id 
        } // Loop on Peak

        NelemYId.at(Id) = Npeak ;
        MeanYId.at(Id) = Sum /Npeak ;
        // cout << Id << "  " << MeanYId.at(Id) << endl;
    } // Loop on Id

    //Now we will sort the MeanY vector and return a vector holding the ID sorted 
    vector<size_t> IndexMeanY(MeanYId.size());
    iota(IndexMeanY.begin(),IndexMeanY.end(),0);
    sort(IndexMeanY.begin(),IndexMeanY.end() , [&](size_t i, size_t j) {
            return MeanYId.at(i) < MeanYId.at(j);
            });

    // Now we have to sort the SetId according to this index ordering
    // Since set can't be sorted we need to copy it to a vector
    vector<int> vSetId(SetId.begin(), SetId.end());
    vector<int> SortedSetId(SetId.begin(), SetId.end());
    for (int Index = 0  ; Index < IndexMeanY.size() ; Index ++){ SortedSetId.at(IndexMeanY.at(Index)) = vSetId.at(Index) ;}

    /*
       cout << endl << " SORTED" <<endl ; 
       for (int i = 0 ; i<IndexMeanY.size() ; i++){
       if (IndexMeanY.at(i) == 49 || IndexMeanY.at(i) == 32){
       cout << NelemYId.at(IndexMeanY.at(i)) << endl; 
       }
       cout << IndexMeanY.at(i) << "  " << SortedSetId.at(i) << "   "  << MeanYId.at(IndexMeanY.at(i)) << endl;
       }

       cout << endl << " *********************" <<endl ; 
       */
    
    //===========================================================================================================
    //                         4.2 Recompacting data into a struct
    //===========================================================================================================
    

    // We now have the different lines sorted by Id we just need to link them together by group of two.
    // to do so we will sort each id by x and then link them

    // 1st step is to sort the original IndexPeak according to the MeanY sort
    // To do so we need to sort the vID vector according to meanY and apply this transformation to IndexPeak too

    // vID as the ID of each peak -> correspondance of peak to a line
    // IndexPeak                  -> index of each of these peak of size number of peak
    // SortedSetId                -> Set of sorted ID of peak according to mean y
    
    //Create a map for the sort of ID
    unordered_map<int,int> order_map ;
    for (size_t i = 0 ; i < SortedSetId.size() ; i++  ) {
        order_map[SortedSetId[i]] = i ; 
    }

    // Combine Index and ID into a single vector of pair
    vector<pair<int,int>> PeakId ;       
    for(size_t i = 0 ; i < IndexPeak.size() ; i++){
        PeakId.emplace_back(IndexPeak.at(i), vID.at(i));
    }

    // Sort the PeakId on the custom order of IDs 
    sort(PeakId.begin(), PeakId.end(), [&order_map](const auto& a, const auto& b){
            return order_map[a.second] < order_map[b.second];
            } );

    //2nd step is to take the original Position vector and sort it thanks to PeakId

    vector<vector<Double_t>> LinePosSorted = *LinePos ;

    //Swap LinePosSorted elemn according to Index of PeakID
    for (int i = 0 ; i < PeakId.size() ; i++ ){
        LinePosSorted.at(i) = LinePos->at(PeakId.at(i).first);
    }    

    // 3rd step is to compile x y id and index in one structure

    if (PeakId.size() != LinePos->size() ){
        cerr << "Vector not same size" << endl; 
    }

    vector<PeakData> PeakStruct ;
    for (size_t i=0 ; i<LinePosSorted.size() ; i++){
        PeakData TempPeak = {LinePosSorted.at(i).at(0), LinePosSorted.at(i).at(1) , PeakId.at(i).second , PeakId.at(i).first};
        PeakStruct.push_back(TempPeak);
    }

    //===========================================================================================================
    //                          4.3 Cleanupi and sortY
    //===========================================================================================================
    // We drop all non clustered element and sort the whole struct by the mean value of Y for their ID.

    // 4th step is to drop all element = -1000
    //

    auto groupBegin = find_if(PeakStruct.begin(),PeakStruct.end(),[&](const PeakData& elem){
            return elem.LineId == -1000 ;
            });

    auto groupEnd = find_if(groupBegin,PeakStruct.end(),[&](const PeakData& elem){
            return elem.LineId != -1000 ;
            });

    PeakStruct.erase(groupBegin,groupEnd);

    //4.1 step we will resort by mean y because it didnt work before , probably due to data structure 

    auto ItMeanY = PeakStruct.begin();
    vector<int> SortedId ;
    vector<Indexed2D> MeanY ;
    int index = 0;

    while (ItMeanY != PeakStruct.end()){

        //Get position of end of the current id
        auto IdEnd = find_if(ItMeanY, PeakStruct.end() , [&](const PeakData& elem){
                return elem.LineId != ItMeanY->LineId; 
                });

        // Fill MeanY of the current Id 

        Double_t Sum  = 0;
        Double_t Npts = 0;

        for (auto i = ItMeanY ; i!=IdEnd ; i++){
            Sum += i->Y;        
            Npts += 1;        
        }

        MeanY.push_back({Sum/Npts ,ItMeanY->LineId  , index });

        index +=1;
        ItMeanY = IdEnd;
    }

    // We now have a all the mean Y of each ID and their index
    // Let's now sort the id from this meanY

    sort(MeanY.begin(),MeanY.end(),[&](Indexed2D a , Indexed2D b){
            return a.x < b.x;
            });

    unordered_map<int,size_t> ordered_map ;
    for (size_t i=0 ; i<MeanY.size() ; i++){
        ordered_map[MeanY.at(i).id]=i;
    }


    sort(PeakStruct.begin(),PeakStruct.end(), [&ordered_map](const PeakData& a , const PeakData &b) {
            return ordered_map[a.LineId]<ordered_map[b.LineId];
            });

    //===========================================================================================================
    //              4.4  Sort by X alternating order
    //===========================================================================================================
    // We sort by X in alternating order because a TCutG need the point to be aligned to link them together

    // Iterate and sort groups with the same LineId
    auto it = PeakStruct.begin();// create an iterator that point to the elements of the PeakData
    bool Oddness = true;             // to know in which way sort
    while (it != PeakStruct.end()){

        // Find range where id is the same
        // This return the first iterator where the condition is fulfilled 
        auto groupEnd = find_if(it,PeakStruct.end(),[lineId = it->LineId](const PeakData& elem){
                return elem.LineId != lineId ;
                });

        // Now we sort this whole group by x

        sort(it,groupEnd,[Oddness](const PeakData& a , const PeakData& b){
                if (Oddness) return a.X < b.X ;         
                if (!Oddness) return a.X > b.X ;        
                else return a.X < b.X;
                });

        // Move to next group
        it = groupEnd;
        Oddness = not(Oddness) ;
    }

    // works ! 
    /*
    for (int i : SortedSetId) {cout << i << endl;}

    int tempId= -1000;
    for (PeakData i : PeakStruct) {
        if (tempId != i.LineId  ){
            tempId = i.LineId;
            cout << " Id " << i.LineId << " Y  " << i.Y << endl;
        }
    } 
    */
    //===========================================================================================================
    //              4.5 Create TCut 
    //===========================================================================================================
    // We need to make sublist for each 2 lines

    // We will iterate through our data and create substructure for each adjacent pair of LineID
    // and store them in a vector of TCutG

    vector<TCutG*> Res ;

    it = PeakStruct.begin();
    int iteration = 0;


    while(it != PeakStruct.end()){

        // 1st line
        int LineIdBottom = it->LineId;
        auto BottomLine = find_if(it,PeakStruct.end(), [&] (const PeakData& elem){
                return elem.LineId != LineIdBottom ;
                }) ;

        //2nd line
        int LineIdTop = BottomLine->LineId;
        auto TopLine = find_if(BottomLine,PeakStruct.end(), [&] (const PeakData& elem){
                return elem.LineId != LineIdTop ;
                }) ;

        //cout << LineIdBottom << "  " << LineIdTop << endl;

        // first we need to find the shortest line
        int SizeBottom = distance(it,BottomLine);
        int SizeTop = distance(BottomLine,TopLine);

        //Subvector with the longest line first then the sortest
        vector<PeakData> CutVec;
        vector<PeakData> LongVec;
        vector<PeakData> ShortVec;

        if(SizeBottom >= SizeTop ){
            CutVec.insert(CutVec.end(),it,TopLine);
            LongVec.insert(LongVec.end(),it,BottomLine);
            ShortVec.insert(ShortVec.end(),BottomLine,TopLine);
        }

        else {
            CutVec.insert(CutVec.end(),BottomLine,TopLine);
            CutVec.insert(CutVec.end(),it,BottomLine);
            LongVec.insert(LongVec.end(),BottomLine,TopLine);
            ShortVec.insert(ShortVec.end(),it,BottomLine);
        }
        


        //nbr of point
        int Npoints = CutVec.size(); 
        //Close the loop by finding the closest point to the last one
        int closestIndexLast = 0;
        int closestIndexBegin = LongVec.size()-1;
        int SizeArray = 0;
        if (ShortVec.size() == 0) {
            closestIndexLast = CutVec.size(); 
            SizeArray = CutVec.size();
        }
        else{
            PeakData lastPoint = ShortVec.back();
            double minDistance = distancePeak(lastPoint,LongVec.at(0));
            for (size_t i = 1 ; i < LongVec.size(); i++){
                double dist = distancePeak(lastPoint,LongVec.at(i));
                if(dist<minDistance){
                    minDistance = dist;
                    closestIndexLast = i;
                }
            }
            PeakData firstPoint = ShortVec.at(0);
            double minDistanceFirst = distancePeak(firstPoint,LongVec.at(0));
            for (size_t i = 1 ; i < LongVec.size(); i++){
                double dist = distancePeak(firstPoint,LongVec.at(i));
                if(dist<minDistanceFirst){
                    minDistanceFirst = dist;
                    closestIndexBegin = i;
                }
            }
            SizeArray = ShortVec.size() + (closestIndexBegin - closestIndexLast);
        }
        //Create array to hold position
        double *x = new double[SizeArray + 1 ];
        double *y = new double[SizeArray + 1 ];

        if (SizeTop == 0){

            for(int i = 0 ; i<CutVec.size() ; i++) {
                x[i] = CutVec.at(i).X ;
                y[i] = CutVec.at(i).Y ;
            }
            x[SizeArray ] = CutVec.at(0).X;
            y[SizeArray ] = CutVec.at(0).Y;

        }

        else {
            for(int i = 0 ; i<ShortVec.size() ; i++) {
                x[i] = ShortVec.at(i).X ;
                y[i] = ShortVec.at(i).Y ;
            }

            for(int i = closestIndexLast ; i<=closestIndexBegin ; i++) {
                x[i + ShortVec.size() - closestIndexLast] = LongVec.at(i).X ;
                y[i + ShortVec.size() - closestIndexLast ] = LongVec.at(i).Y ;
            }



            x[SizeArray] = ShortVec.at(0).X;
            y[SizeArray] = ShortVec.at(0).Y;

        }
        //Saving and moving on

        if (Npoints >10 && (TopLine != PeakStruct.end())){ //Dont link last line to itself 
            TCutG* Cut = new TCutG(Form("Cut %d",iteration) , SizeArray+1, x , y);
            Cut->SetLineWidth(2);
            Cut->SetLineColor(kRed);
            Res.push_back(Cut);
        }
        iteration += 1 ;
        it = BottomLine ;
    } // end while

    return Res;
}


vector<Double_t> diff(const vector<Double_t>& input){
    vector<Double_t> result;
    if (input.size() < 2) return result; // Return empty vector if size < 2

    for (size_t i = 1; i < input.size(); ++i) {
        result.push_back(input[i] - input[i - 1]);
    }
    return result;
}

void SortAndRenameTCutG() {
    TIter next(gROOT->GetListOfSpecials());
    std::vector<std::pair<int, TCutG*>> cuts;

    TObject *obj;
    while ((obj = next())) {
        TCutG *cut = dynamic_cast<TCutG*>(obj);
        if (cut) {
            std::string name = cut->GetName();
            if (name.rfind("Cut", 0) == 0) { // Check if name starts with "Cut"
                try {
                    int num = std::stoi(name.substr(3));
                    cuts.emplace_back(num, cut);
                } catch (...) {
                    continue; // Ignore malformed names
                }
            }
        }
    }

    // Sort by extracted number
    std::sort(cuts.begin(), cuts.end(), [](const auto &a, const auto &b) {
            return a.first < b.first;
            });

    // Rename sequentially
    for (size_t i = 0; i < cuts.size(); ++i) {
        cuts[i].second->SetName(Form("Cut %zu", i));
    }

    std::cout << "TCutG objects renamed from original numbering to 0:N.\n";
}
void SaveTCutGToFile() {
    TFile file("Output/CutAutoZ.root", "RECREATE");
    if (!file.IsOpen()) {
        std::cerr << "Error: Could not open Output/CutAutoZ.root for writing.\n";
        return;
    }

    TIter next(gROOT->GetListOfSpecials());
    TObject *obj;
    while ((obj = next())) {
        TCutG *cut = dynamic_cast<TCutG*>(obj);
        if (cut) {
            cut->Write();
        }
    }

    file.Close();
    std::cout << "All TCutG objects saved to Output/CutAutoZ.root.\n";
}

double distancePeak(const PeakData& p1, const PeakData& p2) {
    return abs(p1.X - p2.X);
}