[AlPhaPha] Added optimisation of AoQ macro

Projects/AlPhaPha/2024/macro/AoQ/DrawAoQ.C

@@ -5,7 +5,7 @@ void DrawAoQ(int sec=8, float time_offset=0){
   chain = new TChain("PhysicsTree");
   //chain->Add("../../root/analysis/run_48_calib.root");
   //chain->Add("../../root/analysis/run_48.root");
-  chain->Add("../../root/analysis/Run247AoQ.root");
+  chain->Add("../../root/analysis/Run247NoAlex.root");
 
   // T13 //
   TString to_draw = Form("FF_Brho/3.10761* (FF_T13+%f)/FF_D13*29.9792*sqrt(1-pow(FF_D13,2)/(pow((FF_T13+%f)*29.9792,2)))>>h(1000,2.5,4.5)",time_offset,time_offset) ;

Projects/AlPhaPha/2024/macro/AoQ/MergeToff.C

+#include <iostream>
+#include <fstream>
+#include <vector>
+#include <string>
+#include <sstream>
+#include <algorithm>
+
+struct Entry {
+    int secNumber;
+    double value;
+};
+
+int MergeToff() {
+    std::vector<Entry> entries;
+    
+    for (int sec = 0; sec <= 21; ++sec) {
+        std::ostringstream filename;
+        filename << "Output/ToffAoQ_Sec_" << sec << ".cal";
+        std::ifstream infile(filename.str());
+        
+        if (!infile) {
+            std::cerr << "Warning: Could not open " << filename.str() << "\n";
+            continue;
+        }
+        
+        std::string line;
+        while (std::getline(infile, line)) {
+            std::istringstream iss(line);
+            std::string header;
+            double value;
+            
+            if (iss >> header >> value) {
+                std::ostringstream expectedHeader;
+                expectedHeader << "ToffAoQ_Sec_" << sec;
+                if (header == expectedHeader.str()) {
+                    entries.push_back({sec, value});
+                    break;
+                }
+            }
+        }
+    }
+    
+    std::sort(entries.begin(), entries.end(), [](const Entry &a, const Entry &b) {
+        return a.secNumber < b.secNumber;
+    });
+    
+    std::ofstream outfile("Output/ToffAoQ.cal");
+    for (const auto &entry : entries) {
+        outfile << "ToffAoQ_Sec_" << entry.secNumber << " " << entry.value << "\n";
+    }
+    
+    std::cout << "Merged file created: ToffAoQ.cal\n";
+    return 0;
+}
+

Projects/AlPhaPha/2024/macro/AoQ/Minimization.C

+TChain *chain;
+double* parameter;
+const int nb_parameter = 11;
+int iteration = 0;
+int m_section = 9;
+
+int NumericalMinimization(const char* minName = "Minuit", const char* algoName = "");
+double ConstantFactor(const double* parameter);
+
+TH2F* h2 = new TH2F("h2","h2",1000,2.5,4.5,500,30,40);
+
+///////////////////////////////////////////////////
+void Minimization(int section=8){
+  m_section = section;
+  chain = new TChain("tree"); 
+  chain->Add("SelectTree.root");
+
+  //double buffer[nb_parameter] = {1,1,1,1,1,1,1,1,1,1,1};
+  double buffer[nb_parameter] = {0.02411,0.8686,0.7199,0.6233,0.4697,0.9787,0.9892,2.1038,1.9429,1.754,2.5};
+  parameter = new double[nb_parameter];
+  parameter = buffer;
+
+  ConstantFactor(parameter);
+  NumericalMinimization("Minuit","Migrad");
+
+  TCanvas* c1 = new TCanvas("c1","c1",800,800);
+  c1->cd();
+  h2->Draw("colz");
+}
+
+
+
+////////////////////////////////////////////////////
+int NumericalMinimization(const char* minName, const char* algoName){
+
+  ROOT::Math::Minimizer* min = ROOT::Math::Factory::CreateMinimizer(minName, algoName);
+
+  min->SetMaxFunctionCalls(1000000);
+  min->SetMaxIterations(10000);
+  min->SetTolerance(0.01);
+  min->SetPrecision(0.01);
+  min->SetPrintLevel(1);
+
+  ROOT::Math::Functor f(&ConstantFactor,nb_parameter);
+
+  min->SetFunction(f);
+
+  min->SetLimitedVariable(0,"scalor",parameter[0],0.01,parameter[0]*0.8,parameter[0]*1.2);
+  for(int i=1; i<nb_parameter; i++){
+    string seg_name = "seg"+to_string(i);
+    min->SetLimitedVariable(i,seg_name,parameter[i],0.01,parameter[i]*0.8,parameter[i]*1.2);
+  }
+  //min->SetLimitedVariable(nb_parameter-1,"scalor",parameter[nb_parameter-1],0.01,0,2);
+
+  min->Minimize();
+
+  const double* xs = min->X();
+
+  ofstream ofile;
+  string filename = "Chio_E_sec" + to_string(m_section) + ".cal";
+  ofile.open(filename.c_str());
+  double init_par[nb_parameter] = {0.02411,0.8686,0.7199,0.6233,0.4697,0.9787,0.9892,2.1038,1.9429,1.754,2.5};
+  cout << "**********************" << endl;
+  cout << "Minimum : " << endl;
+  string parname = "IC_ETOT_SCALING_SEC" + to_string(m_section);
+  ofile << parname << " "  << xs[0] << endl;
+  for(int i=1; i<nb_parameter; i++){ 
+    cout << "par" << i << " = " << xs[i] << " -> " << xs[i] << endl;
+    parname = "IC_SEC"+to_string(m_section)+"_SEG"+to_string(i)+"_ALIGN";
+    ofile << parname << " " << xs[i] << endl;
+  }
+  ofile.close();
+  cout << "MinValue = " << min->MinValue() << endl;
+  ConstantFactor(xs);
+
+  return 0;
+
+}
+
+////////////////////////////////////////////////////
+double ConstantFactor(const double* parameter){
+
+  double distance=0;
+  iteration++;
+
+  double FF_Gamma;
+  double FF_AoQ;
+  double FF_Q;
+  int FPMW_Section;
+
+
+  TICPhysics* IC = new TICPhysics();
+
+  chain->SetBranchStatus("FF_Gamma","true");
+  chain->SetBranchAddress("FF_Gamma",&FF_Gamma);
+
+  chain->SetBranchStatus("FF_AoQ","true");
+  chain->SetBranchAddress("FF_AoQ",&FF_AoQ);
+
+  chain->SetBranchStatus("FF_Q","true");
+  chain->SetBranchAddress("FF_Q",&FF_Q);
+
+  chain->SetBranchStatus("FPMW_Section","true");
+  chain->SetBranchAddress("FPMW_Section",&FPMW_Section);
+
+  chain->SetBranchStatus("fIC","true");
+  chain->SetBranchAddress("IC",&IC);
+
+
+  int nentries = chain->GetEntries();
+
+  h2->Reset();
+  double Etot=0;
+  double M1;
+  double Q;
+  for(int i=0; i<nentries; i++){
+    chain->GetEntry(i);
+    Etot = 0;
+    M1 = 0;
+    Q = 0;
+    if(FPMW_Section==m_section && FF_AoQ>2){
+      for(int j=1; j<nb_parameter; j++){
+        Etot += parameter[j]*IC->fIC[j-1];
+      }
+      Etot = parameter[0]*Etot;
+      M1 = Etot/931.5016/(FF_Gamma-1);
+      Q = M1/FF_AoQ;
+
+      distance += pow((Q - 36),2);
+
+      h2->Fill(FF_AoQ,Q);
+      /*if(cut1->IsInside(FF_AoQ,FF_Q)){
+        if(Q>0) distance += pow((Q - 1410),2);
+      }
+      else if(cut2->IsInside(FF_AoQ,FF_Q)){
+        if(Q>0) distance += pow((Q - 1605),2);
+      }
+      else if(cut3->IsInside(FF_AoQ,FF_Q)){
+        if(Q>0) distance += pow((Q - 1800),2);
+      }*/
+
+      //cout << distance << " " << Etot << " " << M1 << " " << Q << endl;
+    }
+  }
+  if(iteration%1==0){
+    cout << "Number of iterations: " << iteration  << endl;
+    cout << "distance= " << distance << endl;
+  }
+
+  return distance;
+}

Projects/AlPhaPha/2024/macro/AoQ/Output/ReadMe.md

+# Output
+

Projects/AlPhaPha/2024/macro/AoQ/ToffAllSec.sh

+#!/bin/bash
+
+
+#launch job in parralel
+MAX_JOBS=$(( $(nproc) -4 ))   # Maximum parallel jobs
+SEMAPHORE="job_semaphore"
+
+# Create a named pipe (FIFO) for job control
+mkfifo $SEMAPHORE
+exec 3<> $SEMAPHORE
+rm $SEMAPHORE
+
+# Pre-fill the semaphore with tokens equal to MAX_JOBS
+for ((i=0; i<MAX_JOBS; i++)); do
+    echo "" >&3
+done
+
+
+for section in {0..20}; do
+    read -u 3  # Wait for an available slot
+    {
+        root -l -q "./ToffHighRes.C($section,0.0)" 
+        echo "" >&3  # Release the slot
+    } & 
+done
+
+wait
+
+exec 3>&- #close semaphore 
+
+wait 
+
+root -l -q "./MergeToff.C()"

Projects/AlPhaPha/2024/macro/AoQ/ToffHighRes.C

+#include <TCanvas.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TGraph.h>
+#include <TH1.h>
+#include <TSpectrum.h>
+#include <TFitter.h>
+#include <TStopwatch.h>
+#include <iostream>
+#include <TSystem.h>
+#include <TPolyMarker.h>
+#include <fstream>
+
+using namespace std;
+
+
+//===========================================================================================================
+//                              Init mini
+//===========================================================================================================
+
+int iteration_count = 0;
+const int nb_parameter = 1;
+double* parameter;
+int m_section = 9;
+const int Bin = 1500;
+//===========================================================================================================
+
+
+double DistCalc(const double *parameter);
+long double NumericalMinimization(const char* minName, const char* algoName);
+
+void ToffHighRes(int section =8, double initToff = 0.0){
+
+    m_section = section;
+    //===========================================================================================================
+
+    double buffer[nb_parameter] = {initToff};
+    parameter = new double[nb_parameter];
+    parameter = buffer;
+    //===========================================================================================================
+
+    DistCalc(parameter);
+    long double ToffSection = NumericalMinimization("Minuit","Migrad");
+
+    //===========================================================================================================
+    TChain *chain = new TChain("PhysicsTree");
+    chain->Add("../../root/analysis/Run247NoAlex.root");
+ 
+    TH1F * hOriginal= new TH1F ("hOriginal",Form("Section_%d Original",m_section),Bin,2.5,4.5);
+    TH1F *hCorrected = new TH1F("hCorrected" , Form("Section_%d Corrected",m_section),Bin,2.5,4.5) ;
+    TString to_draw = Form("FF_Brho/3.10761* (FF_T13+%f)/FF_D13*29.9792*sqrt(1-pow(FF_D13,2)"
+            "/(pow((FF_T13+%f)*29.9792,2)))>>hOriginal"
+            ,parameter[0],parameter[0]) ;
+
+    TString to_draw2 = Form("FF_Brho/3.10761* (FF_T13+%Lf)/FF_D13*29.9792*sqrt(1-pow(FF_D13,2)"
+            "/(pow((FF_T13+%Lf)*29.9792,2)))>>hCorrected"
+            ,ToffSection,ToffSection) ;
+
+    TString cond = Form("FPMW_Section==%d",m_section);
+    //===========================================================================================================
+
+    TCanvas* c1 = new TCanvas("c2","c2",800,800);
+    c1->Divide(2);
+    c1->cd(1);
+    chain->Draw(to_draw,cond,"");
+    gPad->SetGrid();
+    c1->cd(2);
+    chain->Draw(to_draw2,cond,"");
+    gPad->SetGrid();
+
+    //===========================================================================================================
+    //                                      Saving
+    //===========================================================================================================
+
+    string filename = "Output/ToffAoQ_Sec_" +to_string(m_section)+".cal";
+    ofstream outfile(Form("Output/ToffAoQ_Sec_%d.cal",m_section)) ;
+    string target_prefix = "ToffAoQ_Sec_" + to_string(m_section);
+       
+    if(outfile.is_open()){
+        outfile << target_prefix << " " << ToffSection;
+        outfile.close();
+    }
+}
+
+long double NumericalMinimization(const char* minName, const char* algoName){
+
+    ROOT::Math::Minimizer* min = ROOT::Math::Factory::CreateMinimizer(minName, algoName);
+
+    min->SetMaxFunctionCalls(1000000);
+    min->SetMaxIterations(10000);
+    min->SetTolerance(0.01);
+    min->SetPrecision(0.01);
+    min->SetPrintLevel(1);
+
+    ROOT::Math::Functor f(&DistCalc,nb_parameter);
+
+    min->SetFunction(f);
+
+    double step = 0.1 , lowerBound = -1.2 , higherBound = 1.2;
+    min->SetLimitedVariable(0,Form("ToffSection_%d",m_section),parameter[0],step,lowerBound,higherBound);
+    min->Minimize();
+
+    const double* xs = min->X();
+
+    DistCalc(xs);
+    return xs[0];
+
+}
+
+
+double DistCalc(const  double* parameter){
+
+    long double distance = 0;
+    iteration_count++;  // Increment iteration count
+
+    //===========================================================================================================
+    //                                      LoadTree
+    //===========================================================================================================
+    //===========================================================================================================
+
+    TChain *chain = new TChain("PhysicsTree");
+    chain->Add("../../root/analysis/Run247NoAlex.root");
+    //===========================================================================================================
+    //                              Draw hist
+    //===========================================================================================================
+
+    TH1F *hist = (new TH1F("hist" , "hist" , Bin , 2.5 , 4.5));
+
+    TString to_draw = Form("FF_Brho/3.10761* (FF_T13+%f)/FF_D13*29.9792*sqrt(1-pow(FF_D13,2)/(pow((FF_T13+%f)*29.9792,2)))>>hist"
+            ,parameter[0],parameter[0]) ;
+    TString cond = Form("FPMW_Section==%d",m_section);
+    chain->Draw(to_draw,cond,"");
+
+    //===========================================================================================================
+    TSpectrum PeakFinder(30);
+    Double_t source[Bin], dest[Bin];
+    // Load hist in a array
+    for (int i=0 ; i< Bin ; i++){
+        source[i] = hist->GetBinContent(i+1);
+    }
+    //===========================================================================================================
+    // find peaks
+    int nPeaks = PeakFinder.SearchHighRes(source,dest,Bin,3,30,kTRUE,20,kTRUE,10);
+    Double_t *xPeaks = PeakFinder.GetPositionX(); //peak position in bin number
+    vector<Double_t> fPositionX(nPeaks), fPositionY(nPeaks); // Vector to hold peak position in X and Y with size number of peak
+    Double_t fPositionXarr[nPeaks], fPositionYarr[nPeaks]; // Vector to hold peak position in X and Y with size number of peak
+    vector<double> peakClose3, peakClose4;
+
+
+    for (int i = 0; i < nPeaks; i++) {
+        Double_t a = xPeaks[i];
+        int bin = 1 + Int_t(a + 0.5);
+        fPositionX[i] = hist->GetBinCenter(bin);
+        fPositionXarr[i] = hist->GetBinCenter(bin);
+        fPositionY[i] = hist->GetBinContent(bin);
+        fPositionYarr[i] = hist->GetBinContent(bin);
+
+        // find if closest to 3 or 4
+        (round(fPositionX[i]) == 3 ? peakClose3 : peakClose4).push_back(fPositionX[i]) ;
+    }
+
+
+    //===========================================================================================================
+    // find peak nearest to 3 if there is at least 2 peak near 3
+    long double bestPeak3 = 1e6;
+    long double bestPeak4 = 1e6;
+
+    if (peakClose3.size()>2){
+        for (int j=0; j<peakClose3.size(); j++){
+            //cout << "Curr dist " << fabs(xPeaks[j] - 3.0) << endl;
+            //cout << "Curr position " << peakClose3[j]  << endl;
+            //cout << "Best dist" <<fabs(bestPeak - 3.0) << endl;
+            if (fabs(peakClose3[j] - 3.0 ) < fabs(bestPeak3 -3.0)){
+                bestPeak3 = peakClose3[j];
+            } // end cond nearer
+        }// end loop peaks
+    }
+
+    if (peakClose4.size()>2){
+        for (int j=0; j<peakClose4.size(); j++){
+            //cout << "Curr dist " << fabs(xPeaks[j] - 3.0) << endl;
+            //cout << "Curr position " << xPeaks[j]  << endl;
+            //cout << "Best dist" <<fabs(bestPeak - 3.0) << endl;
+            if (fabs(peakClose4[j] - 4.0 ) < fabs(bestPeak4 -4.0)){
+                bestPeak4 = peakClose4[j];
+            } // end cond nearer
+        }// end loop peaks
+    }
+
+
+    //===========================================================================================================
+
+    cout << " ************************ " << iteration_count  << " *************************** "  << endl;
+    cout << "Toff " << parameter[0] << endl;
+    if (peakClose3.size()>2){
+        cout << "Best peak3 " << bestPeak3 << endl;
+    }
+    if (peakClose4.size()>2){
+        cout << "Best peak4 " << bestPeak4 << endl;
+    }
+    //===========================================================================================================
+    distance += fabs( (peakClose3.size() > 2 ?  bestPeak3 : 3.0) - 3.0) 
+        +fabs( (peakClose4.size() > 2 ?  bestPeak4 : 4.0) - 4.0); 
+    /*
+    //if result check is needed
+    TPolyMarker *pm = (TPolyMarker *)hist->GetListOfFunctions()->FindObject("TPolyMarker");
+    if (pm) {
+    hist->GetListOfFunctions()->Remove(pm);
+    delete pm;
+    }
+    pm = new TPolyMarker(nPeaks, fPositionXarr, fPositionYarr);
+    hist->GetListOfFunctions()->Add(pm);
+    pm->SetMarkerStyle(23);
+    pm->SetMarkerColor(kRed);
+    pm->SetMarkerSize(1.3);
+
+    TH1F *d = new TH1F("d", "", Bin, 2.5, 4.5);
+    for (int i = 0; i < Bin; i++){
+    d->SetBinContent(i + 1, dest[i]);
+    }
+    d->SetLineColor(kRed);
+
+    hist->Draw();
+    d->Draw("SAME");
+    //d->Delete();
+    */
+
+    hist->Delete();
+
+    // Garbage collector
+    chain->Reset();
+    delete chain;
+
+    return distance;
+
+}
+
+