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