From 631beb7ce903b772070ca0a499ed36b6be63ad2c Mon Sep 17 00:00:00 2001 From: ddobrigk Date: Sat, 8 Nov 2025 17:48:37 +0100 Subject: [PATCH 1/9] Add code for glauber fits --- Common/Tools/Multiplicity/README.md | 17 + Common/Tools/Multiplicity/macros/README.md | 69 ++++ .../Multiplicity/macros/runCalibration.C | 150 ++++++++ .../Tools/Multiplicity/macros/runGlauberFit.C | 352 ++++++++++++++++++ .../Multiplicity/macros/saveCorrelation.C | 56 +++ 5 files changed, 644 insertions(+) create mode 100644 Common/Tools/Multiplicity/README.md create mode 100644 Common/Tools/Multiplicity/macros/README.md create mode 100644 Common/Tools/Multiplicity/macros/runCalibration.C create mode 100755 Common/Tools/Multiplicity/macros/runGlauberFit.C create mode 100644 Common/Tools/Multiplicity/macros/saveCorrelation.C diff --git a/Common/Tools/Multiplicity/README.md b/Common/Tools/Multiplicity/README.md new file mode 100644 index 00000000000..30a2d725fda --- /dev/null +++ b/Common/Tools/Multiplicity/README.md @@ -0,0 +1,17 @@ +# Multiplicity/centrality tools in O2/O2Physics + +This directory serves to aggregate all files necessary for multiplicity +and centrality calibration going from pp to Pb-Pb. It offers functionality +to perform simple slicing in percentiles, such as what is done in +proton-proton collisions, as well as the Glauber + analytical NBD fitter +used to perform anchoring and hadronic event distribution estimates in +nucleus-nucleus collisions. + +## Files +* `README.md` this readme +* `CMakeLists.txt` definition of source files that need compiling +* `multCalibrator.cxx/h` a class to do percentile slicing of a given histogram. Used for all systems. +* `multMCCalibrator.cxx/h` a class to perform data-to-mc matching of average Nch. +* `multGlauberNBDFitter.cxx/h` a class to do glauber fits. +* `multModule.h` a class to perform calculations of multiplicity and centrality tables for analysis. Meant to be used inside the main core service wagon 'multcenttable'. + diff --git a/Common/Tools/Multiplicity/macros/README.md b/Common/Tools/Multiplicity/macros/README.md new file mode 100644 index 00000000000..0054d665d71 --- /dev/null +++ b/Common/Tools/Multiplicity/macros/README.md @@ -0,0 +1,69 @@ +# Example multiplicity calibration macros + +You will find some example macros in this directory that will allow for the calculation of a glauber fit and the corresponding calibration histograms. + +A simplified description of the procedure necessary to generate a +Glauber + NBD fit to a certain histogram is described below. + +## Performing a Glauber + NBD fit + +### First step: calculation of Glauber MC sample + +The Glauber + NBD model assumes that the multiplicity / amplitude +distribution that one intends to fit can be described as coming +from a certain number N_{ancestor} of particle-emitting sources called 'ancestors'. Each ancestor is assumed to produce particles +according to a negative binominal distribution. Further, +the number of ancestors is assumed to be related to the basic +glauber quantities N_{part} and N_{coll} via the formula: + +N_{ancestor} = f * N_{part} + (1-f) N_{coll} + +Usually, the value f is fixed at 0.8, and thus the range of +N_{ancestors} is typically 0-700 in Pb-Pb collisions. + +In order to allow for Glauber + NBD fitting, the correlation of +(N_{part}, N_{coll}) needs to be known. For that purpose, +a tree of Glauber MC needs to be generated using TGlauberMC using the relevant nuclei, which also involves choosing an appropriate nuclear profile. + +Once TGlauberMC has been used to produce a tree with N_{part} and +N_{coll} values, their correlation needs to be saved to a 2D histogram that serves as input to the Glauber + NBD fitter used in +O2/O2Physics. This is done using the macro called `saveCorrelation.C` contained in this directory, which produces a file named `baseHistos.root`. The file `saveCorrelation.C` serves as +an example for the Pb-Pb case and minor adaptation may be needed +in case other nuclei are to be used. + +### Second step: execution of Glauber + NBD fitter + +The fitting procedure is done via the macro ``. Because +the numerical fitting utilizes a convolution of a N_{ancestor} +distribution and NBD distributions, it will not be as fast +as typical one-dimensional fits: typical times of 10-100 seconds +are not unusual. The macro will produce an output file +that contains: + +* the original fitted histogram +* the actual glauber fit distribution, plotted all the way +to zero multiplicity + +This output can then be used to extract percentiles. + +### Third step: calculation of percentile boundaries + +Once both the data and glauber fit distributions are known, +the next step is to create a 'stitched' data/glauber distribution in +which the distribution at low multiplicities follows +the glauber fit and the distribution at higher multiplicities +follows the data. The multiplicity value in which the switch from +glauber to data takes place is called 'anchor point'. + +Because of the fact that this 'stitching' procedure may need to be tuned and the actual glauber fit is slow, the stitching is done +in a separate macro called `runCalibration.C`. It is provided +in this directory as well and it is the third and last step in calculating percentile boundaries. The macro will printout some +key boundaries as well as save an output file with the calibration. + +*Bonus*: at the stage in which the glauber fit has been done +and all information is available, a de-convolution process +can be used to calculate the average N_{part} and N_{coll} +in centrality slices. This functionality is also provided +in O2Physics as part of the `multGlauberNBDFitter` and the +`runCalibration.C` macro can optionally also perform that +deconvolution. *Warning*: this procedure might take a mooment. \ No newline at end of file diff --git a/Common/Tools/Multiplicity/macros/runCalibration.C b/Common/Tools/Multiplicity/macros/runCalibration.C new file mode 100644 index 00000000000..ca029570cb7 --- /dev/null +++ b/Common/Tools/Multiplicity/macros/runCalibration.C @@ -0,0 +1,150 @@ +void runCalibration( TString lInputFileName = "results/AR_544122_glauberNBD_ancestorMode2_hFT0C_BCs.root", double anchorPointPercentage = 90.0, double matchRange = 200.0, bool doNpartNcoll = false){ + TFile *file = new TFile(lInputFileName.Data(), "READ"); + file->ls(); + + TH1F *hData = (TH1F*) file->Get("hV0MUltraFine"); + TH1F *hGlauberParameters = (TH1F*) file->Get("hGlauberParameters"); + TH1F *hGlauberFitRange = (TH1F*) file->Get("hGlauberFitRange"); + hData->SetName("hData"); + TH1F *hStitched = (TH1F*) hData->Clone("hStitched"); + TH1F *hFit = (TH1F*) file->Get("hGlauber"); + + TCanvas *c1 = new TCanvas("c1", "", 800, 600); + c1->SetLeftMargin(0.17); + c1->SetBottomMargin(0.17); + c1->SetRightMargin(0.15); + c1->SetTopMargin(0.05); + c1->SetTicks(1,1); + c1->SetLogz(); + c1->SetFrameFillStyle(0); + c1->SetFillStyle(0); + + + cout<<"Data bin width: "<GetBinWidth(1)<GetBinWidth(1)<GetNbinsX()+1; ii++){ + // renormalize data curve + int bin1 = ii+1; + int bin2 = hData->FindBin( hData->GetBinLowEdge(ii+1) + matchRange + 1e-3 ); + double matchRangeData = hData -> Integral( bin1, bin2); + double matchRangeFit = hFit -> Integral( bin1, bin2); + + // rescale fit to match in the vicinity of the region we're at + hFit->Scale(matchRangeData/matchRangeFit); + + double integralFit = hFit->Integral(1,ii); + double integralData = hData->Integral(ii+1,hData->GetNbinsX()+1); + double integralAll = integralFit+integralData; + + cout<<"at bin #"<GetBinLowEdge(ii+1)<<" fraction above this value is: "<GetBinLowEdge(ii+1); + + if(integralData/integralAllGetNbinsX()+1; ii++){ + // renormalize data curve + if( hData->GetBinCenter(ii) < anchorPoint) hStitched->SetBinContent(ii, hFit->GetBinContent(ii)); + } + + + cout<<"Anchor point determined to be: "<SetLineColor(kRed); + hStitched->SetLineColor(kBlue); + + hData->GetYaxis()->SetTitleSize(0.055); + hData->GetXaxis()->SetTitleSize(0.055); + hData->GetYaxis()->SetLabelSize(0.04); + hData->GetXaxis()->SetLabelSize(0.04); + hData->SetTitle(""); + hData->Draw("hist"); + hFit->Draw("hist same"); + hStitched->Draw("hist same"); + + //All fine, let's try the calibrator + multCalibrator *lCalib = new multCalibrator("lCalib"); + lCalib->SetAnchorPointPercentage(100.0f); + lCalib->SetAnchorPointRaw(-1e-6); + + //Set standard Pb-Pb boundaries + lCalib->SetStandardOnePercentBoundaries(); + + TString calibFileName = lInputFileName.Data(); + calibFileName.ReplaceAll("glauberNBD", "calibration"); + TFile *fileCalib = new TFile(calibFileName.Data(), "RECREATE"); + + TH1F *hCalib = lCalib->GetCalibrationHistogram(hStitched, "hCalib"); + + TCanvas *c2 = new TCanvas("c2", "", 800, 600); + c2->SetLeftMargin(0.17); + c2->SetBottomMargin(0.17); + c2->SetRightMargin(0.15); + c2->SetTopMargin(0.05); + c2->SetTicks(1,1); + // c2->SetLogz(); + c2->SetFrameFillStyle(0); + c2->SetFillStyle(0); + + hCalib->GetYaxis()->SetTitleSize(0.055); + hCalib->GetXaxis()->SetTitleSize(0.055); + hCalib->GetYaxis()->SetLabelSize(0.04); + hCalib->GetXaxis()->SetLabelSize(0.04); + hCalib->SetTitle(""); + hCalib->Draw(); + + + fileCalib->cd(); + + hData->Write(); + hCalib->Write(); + hStitched->Write(); + hFit->Write(); + + if(doNpartNcoll){ + cout<<"Will now attempt to calculate % -> Np, Nc map..."<GetGlauberNBD(); + + //Step 1: open the (Npart, Ncoll) pair information, provide + TFile *fbasefile = new TFile("basehistos.root","READ"); + TH2D *hNpNc = (TH2D*) fbasefile->Get("hNpNc"); + g->SetNpartNcollCorrelation(hNpNc); + g->InitializeNpNc(); + + fitfunc->SetParameter(0, hGlauberParameters->GetBinContent(1)); + fitfunc->SetParameter(1, hGlauberParameters->GetBinContent(2)); + fitfunc->SetParameter(2, hGlauberParameters->GetBinContent(3)); + fitfunc->SetParameter(3, hGlauberParameters->GetBinContent(4)); + fitfunc->SetParameter(4, hGlauberParameters->GetBinContent(5)); + + Double_t lMax = hData->GetBinLowEdge( hData->GetNbinsX() + 1); + + // uncomment if Np Nc needed -> Warning, slow! + g->CalculateAvNpNc( hProfileNpart, hProfileNcoll, h2dNpart, h2dNcoll, hCalib, 0, lMax ); + + hProfileNpart->Write(); + hProfileNcoll->Write(); + h2dNpart->Write(); + h2dNcoll->Write(); + } + + fileCalib->Write(); +} diff --git a/Common/Tools/Multiplicity/macros/runGlauberFit.C b/Common/Tools/Multiplicity/macros/runGlauberFit.C new file mode 100755 index 00000000000..df1041da9fe --- /dev/null +++ b/Common/Tools/Multiplicity/macros/runGlauberFit.C @@ -0,0 +1,352 @@ +#include "multCalibrator.h" +#include "multGlauberNBDFitter.h" +#include "TFile.h" +#include "TH1F.h" +#include "TH2F.h" +#include "TF1.h" +#include "TCanvas.h" +#include "TLegend.h" +#include "TLatex.h" +#include "TLine.h" +#include "TDirectory.h" +#include "TStyle.h" +#include "TSystem.h" +#include "TGraph.h" +#include "TTree.h" +#include "TStopwatch.h" + +//________________________________________________________________ +Double_t FastIntegrate(TF1 *f1, Double_t a, Double_t b, Int_t n = 5){ + //Do fast integration with N sampling points + const Int_t nc = n; + Double_t x[nc], y[nc]; + Double_t lWidth = (b-a)/((double)(n-1)); + for(Int_t ii=0; iiEval( x[ii] ); + } + //Now go via trapezoids, please (this probably has a name) + Double_t lIntegral = 0; + for(Int_t ii=0; iiGetNbinsX(); + Double_t lCountDesired = lPercentile * histo->GetEntries()/100; + Long_t lCount = 0; + for(Long_t ibin=1;ibinGetBinContent(ibin); + if(lCount >= lCountDesired){ + //Found bin I am looking for! + Double_t lWidth = histo->GetBinWidth(ibin); + Double_t lLeftPercentile = 100.*(lCount - histo->GetBinContent(ibin))/histo->GetEntries(); + Double_t lRightPercentile = 100.*lCount / histo->GetEntries(); + + Double_t lProportion = (lPercentile - lLeftPercentile)/(lRightPercentile-lLeftPercentile); + + lReturnValue = histo->GetBinLowEdge(ibin) + lProportion*lWidth; + break; + } + } + return lReturnValue; +} + +// master fit function +// parameters: +// --- input file name (from the grid, typically centrality-studies task) +// --- histogram name: histogram to use within the input file +// --- ancestor mode: 0: truncation, 1: rounding, 2: effective / non-integer (default: 2) +// --- free k: keep k value free (default Pb-Pb: fixed at 1.5) +// --- use dMu/dNanc: allow for a varying production of Nch vs ancestor if Nancestor is large. Models detector saturation in effective manner. +// --- free f: keep f value free (default Pb-Pb: fixed at 0.8) +// --- f value: the value to use for fixed f + +int runGlauberFit(TString lInputFileName = "AnalysisResultsLHC24ar.root", TString histogramName = "hFT0C_BCs", int ancestorMode = 2, Bool_t lFreek = kFALSE, Bool_t use_dMu_dNanc = kFALSE, Bool_t lFreef = kFALSE, Float_t lfvalue = 0.800, TString outputFile = "output.root") { + gStyle->SetLineScalePS(1); + gStyle->SetOptStat(0); + + cout<<"Starting!"< Get(Form("centrality-study/%s", histogramName.Data())); + + // disregard bin zero + cout<<"Received bin zero content: "<< hV0Mfine ->GetBinContent(0)<<", will set to zero..."<SetBinContent(0, 0); + + if(!hV0Mfine) cout<<"Problem with histogram!"<GetEntries()<GetNbinsX()<GetBinLowEdge(hV0Mfine->GetNbinsX()+1)<Clone("hV0M"); + TH1F *hV0MUltraFine = (TH1F*) hV0Mfine->Clone("hV0MUltraFine"); + hV0M->SetName("hV0M"); + hV0M->SetTitle(""); + + //____________________________________________ + // maximum fit range estimate (avoid tails) + // may need adjusting + Double_t lFitRangeMax = GetBoundaryForPercentile(hV0Mfine, 0.008); + cout<<"Fit range max estimated from histogram: "<Rebin(rebinFactor); + + //____________________________________________ + // simple plots for inspection + TCanvas *c1 = new TCanvas("c1", "", 1300,900); + c1->SetFrameFillStyle(0); + c1->SetFillStyle(0); + c1->Divide(1,2); + c1->cd(1)->SetFrameFillStyle(0); + c1->cd(1)->SetFillStyle(0); + c1->cd(2)->SetFrameFillStyle(0); + c1->cd(2)->SetFillStyle(0); + + c1->cd(1); + c1->cd(1)->SetLogy(); + c1->cd(1)->SetTicks(1,1); + c1->cd(1)->SetPad(0,0.5,1,1); + c1->cd(2)->SetPad(0,0.0,1,.5); + + c1->cd(1)->SetBottomMargin(0.001); + c1->cd(1)->SetRightMargin(0.25); + c1->cd(1)->SetTopMargin(0.02); + c1->cd(1)->SetLeftMargin(0.07); + + c1->cd(2)->SetBottomMargin(0.14); + c1->cd(2)->SetRightMargin(0.25); + c1->cd(2)->SetTopMargin(0.001); + c1->cd(2)->SetLeftMargin(0.07); + c1->cd(2)->SetTicks(1,1); + c1->cd(1); + + hV0M->GetXaxis()->SetRangeUser(0,lFitRangeMax); + hV0M->GetYaxis()->SetRangeUser(0.25,hV0M->GetMaximum()*3); + hV0M->SetLineColor(kBlack); + hV0M->SetMarkerStyle(20); + hV0M->SetMarkerColor(kBlack); + hV0M->SetMarkerSize(0.5); + hV0M->GetYaxis()->SetTitleSize(0.07); + hV0M->GetYaxis()->SetLabelSize(0.05); + hV0M->GetYaxis()->SetTitle("Count"); + hV0M->GetYaxis()->SetTitleOffset(0.5); + hV0M->GetXaxis()->SetLabelSize(0.05); + hV0M->GetXaxis()->SetTitleSize(0.06); + hV0M->GetXaxis()->SetTitle("FT0A+C Amplitude"); + hV0M->GetYaxis()->SetTickLength(0.015); + hV0M->SetStats(0); + hV0M->Draw("E"); + + //Stand back! Imma gonna do GLAUBER FITTIN' + multGlauberNBDFitter *g = new multGlauberNBDFitter("lglau"); + g->SetAncestorMode(ancestorMode); + + //Step 1: open the (Npart, Ncoll) pair information, provide + TFile *fbasefile = new TFile("basehistos.root","READ"); + TH2D *hNpNc = (TH2D*) fbasefile->Get("hNpNc"); + // return to proper scope + fOutput->cd(); + g->SetNpartNcollCorrelation(hNpNc); + g->SetInputV0M(hV0M); + g->SetFitRange(lFitRange,lFitRangeMax); + //Step 3: go for it ... + g->SetNorm(1.53527e+08); + TString lString = "REM0"; + g->SetFitOptions(lString.Data()); + g->SetFitNpx(100000); + + TF1 *fitfunc = g->GetGlauberNBD(); + + //____________________________________________ + // + // set initial fit parameters here + // may require manual tuning depending on data! + + Double_t guessedMu = lFitRangeMax/53968.4 * 0.175*3.53971e+02; + cout<<"Guessed GlauberNBD mu value: "<SetParameter(0,guessedMu); //mu value + fitfunc->SetParLimits(0,0.25*guessedMu,guessedMu*2); + + if(!lFreek){ + fitfunc->FixParameter(1,1.5); //k value + }else{ + fitfunc->SetParLimits(1,0.01,35); + fitfunc->SetParameter(1,1.5); + } + if(!lFreef){ + fitfunc->FixParameter(2,lfvalue); // f value + }else{ + fitfunc->SetParLimits(2,0.55,0.97); + fitfunc->SetParameter(2, 0.800); + } + fitfunc->SetParLimits(3,0.1e+5, 800e+10); // normalization + fitfunc->SetParameter(3, 0.80832e+08); + + //dMu/dNanc + fitfunc->SetParameter(4,0); + if( !use_dMu_dNanc ){ + fitfunc->FixParameter(4,0); + } + + //fitfunc->SetParameter(4,-1.15443e-01); + //fitfunc->SetParameter(4,-4.55957e-02); + + //dk/dNanc + fitfunc->FixParameter(5,0); + //fitfunc->SetParameter(5,1.63590e-03); + + //d2Mu/dNanc2 + fitfunc->FixParameter(6,0.0); + //fitfunc->SetParameter(6,4.02271e-05); + + fitfunc->FixParameter(7,0.0); + //fitfunc->SetParameter(7,-1.24349e-06); + + //____________________________________________ + // handle internals for fitting: needs to + // be done before the fit is attempted! + g->InitializeNpNc(); + g->InitAncestor(); + + cout<<"WILL NOW ATTEMPT GLAUBER FIT"<DoFit(); + Int_t lAttempts = 1; + while(lAttempts < 10 && lFitStatus==0){ + // insist on fitting until it works + cout<<"Attempting fit again ("<DoFit(); + } + cout<<"Final fit status: "<SetOptStat(0); + //Do a ratio plot + TH1D *hGlauber = (TH1D*) hV0MUltraFine->Clone("hGlauber"); + TH1D *hRatio = (TH1D*) hV0MUltraFine->Clone("hRatio"); + hGlauber->Reset(); + + c1->cd(1); + fitfunc->SetLineColor(kRed); + fitfunc->SetLineWidth(2); + fitfunc->Draw("same"); + + cout<<"Calculating glauber function histogram with the same binning as data input... please wait..."<GetNbinsX()+1; ii++){ + Double_t lFuncVal = FastIntegrate( fitfunc, hGlauber->GetBinLowEdge(ii), hGlauber->GetBinLowEdge(ii+1), 4); + hGlauber->SetBinContent(ii, lFuncVal); + if(ii%500==0){ + cout<<"At integration #"<GetNbinsX()+1<<"..."<cd(2); + Float_t lLoRangeRatio = 0.35; + Float_t lHiRangeRatio = 1.65; + hRatio->Divide(hGlauber); + hRatio->GetYaxis()->SetTitle("Data/Fit"); + hRatio->GetXaxis()->SetTitle("FT0C amplitude"); + hRatio->GetYaxis()->SetTitleSize(0.055); + hRatio->GetYaxis()->SetTitleOffset(0.7); + hRatio->GetXaxis()->SetTitleSize(0.055); + hRatio->GetYaxis()->SetLabelSize(0.045); + hRatio->GetXaxis()->SetLabelSize(0.045); + hRatio->GetYaxis()->SetRangeUser(lLoRangeRatio,lHiRangeRatio); + hRatio->GetXaxis()->SetRangeUser(0,lFitRangeMax); + hRatio->SetMarkerStyle(20); + hRatio->SetMarkerColor(kGray+2); + hRatio->SetLineColor(kGray+2); + //hRatio->SetMarkerSize(1.0); + hRatio->SetMarkerSize(.7); + hRatio->SetStats(0); + //hRatioWide->SetStats(0); + + hRatio->Draw("hist"); + + TLine *line = new TLine(0,1,lFitRangeMax,1); + line->SetLineStyle(7); + line->SetLineColor(kGray+1); + line->Draw(); + + TLine *lFitRangeLine = new TLine( lFitRange, lLoRangeRatio, lFitRange, 0.9) ; + lFitRangeLine->SetLineColor(kBlue); + lFitRangeLine->SetLineWidth(1); + lFitRangeLine->SetLineStyle(2); + lFitRangeLine->Draw(); + + TH1D *hRatioGrayed = (TH1D*) hRatio->Clone("hRatioGrayed"); + hRatioGrayed->SetMarkerColor(kGray+2); + hRatioGrayed->SetLineColor(kGray+2); + hRatioGrayed->Draw("same"); + + hRatio->SetLineWidth(1); + hRatio->Draw("same hist"); + //hRatioWide->Draw("same"); + + c1->cd(1); + TLatex *lat = new TLatex(); + lat->SetNDC(); + Float_t lPosText = 0.76; + Float_t lYShift = 0.25; + lat->SetTextSize(0.042); + + // save the glauber parameters explicitly + TH1D *hGlauberParameters = new TH1D("hGlauberParameters", "", 10,0,10); + TH1D *hGlauberFitRange = new TH1D("hGlauberFitRange", "", 10,0,10); + + //fitfunc + hGlauberParameters -> SetBinContent( 1, fitfunc -> GetParameter(0)); + hGlauberParameters -> SetBinContent( 2, fitfunc -> GetParameter(1)); + hGlauberParameters -> SetBinContent( 3, fitfunc -> GetParameter(2)); + hGlauberParameters -> SetBinContent( 4, fitfunc -> GetParameter(3)); + hGlauberParameters -> SetBinContent( 5, fitfunc -> GetParameter(4)); + hGlauberParameters -> Write(); + + Double_t lLoRangeGlauber, lHiRangeGlauber; + fitfunc->GetRange(lLoRangeGlauber, lHiRangeGlauber); + hGlauberFitRange->SetBinContent(1, lLoRangeGlauber); + hGlauberFitRange->SetBinContent(2, lHiRangeGlauber); + hGlauberFitRange->Write(); + + hRatio->Write(); + fOutput->Write(); + + return 0; +} diff --git a/Common/Tools/Multiplicity/macros/saveCorrelation.C b/Common/Tools/Multiplicity/macros/saveCorrelation.C new file mode 100644 index 00000000000..9e9e5eb5776 --- /dev/null +++ b/Common/Tools/Multiplicity/macros/saveCorrelation.C @@ -0,0 +1,56 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +// +// This file provides a simple macro to convert TGlauberMC output tuples +// into the 2D correlation histogram necessary for the ALICE machinery +// that performs Glauber + NBD fits. + +void saveCorrelation(TString filename = "gmc-PbPb-snn68.21-md0.40-nd-1.0-rc1-smax99.0.root", TString outputFile = "basehistos.root"){ + TFile *fin = new TFile(filename.Data(),"READ"); + TNtuple *ntup = (TNtuple*) fin->Get("nt_Pb_Pb"); + + // try other Pb nuclear profiles in case "Pb" - "Pb" not found + if(!ntup){ + ntup = (TNtuple*) fin->Get("nt_Pbpn_Pbpn"); + } + if(!ntup){ + ntup = (TNtuple*) fin->Get("nt_PbpnVar1_PbpnVar1"); + } + if(!ntup){ + ntup = (TNtuple*) fin->Get("nt_PbpnVar2_PbpnVar2"); + } + if(!ntup){ + ntup = (TNtuple*) fin->Get("nt_PbpnVar3_PbpnVar3"); + } + if(!ntup){ + ntup = (TNtuple*) fin->Get("nt_PbpnVar4_PbpnVar4"); + } + + if(!ntup){ + cout<<"No tree found!"<GetEntries()<SetTicks(1,1); + c1->SetLogz(); + ntup->Draw("Ncoll:Npart>>hNpNc", "", "colz"); + fout->Write(); +} From 547767367bfeecd45ba8fd7c8b529d725a73e356 Mon Sep 17 00:00:00 2001 From: ALICE Builder Date: Sat, 8 Nov 2025 17:50:16 +0100 Subject: [PATCH 2/9] Please consider the following formatting changes (#495) --- .../Multiplicity/macros/runCalibration.C | 160 ++++---- .../Tools/Multiplicity/macros/runGlauberFit.C | 350 +++++++++--------- .../Multiplicity/macros/saveCorrelation.C | 55 +-- 3 files changed, 288 insertions(+), 277 deletions(-) diff --git a/Common/Tools/Multiplicity/macros/runCalibration.C b/Common/Tools/Multiplicity/macros/runCalibration.C index ca029570cb7..a8611c9a5ec 100644 --- a/Common/Tools/Multiplicity/macros/runCalibration.C +++ b/Common/Tools/Multiplicity/macros/runCalibration.C @@ -1,68 +1,69 @@ -void runCalibration( TString lInputFileName = "results/AR_544122_glauberNBD_ancestorMode2_hFT0C_BCs.root", double anchorPointPercentage = 90.0, double matchRange = 200.0, bool doNpartNcoll = false){ - TFile *file = new TFile(lInputFileName.Data(), "READ"); +void runCalibration(TString lInputFileName = "results/AR_544122_glauberNBD_ancestorMode2_hFT0C_BCs.root", double anchorPointPercentage = 90.0, double matchRange = 200.0, bool doNpartNcoll = false) +{ + TFile* file = new TFile(lInputFileName.Data(), "READ"); file->ls(); - - TH1F *hData = (TH1F*) file->Get("hV0MUltraFine"); - TH1F *hGlauberParameters = (TH1F*) file->Get("hGlauberParameters"); - TH1F *hGlauberFitRange = (TH1F*) file->Get("hGlauberFitRange"); + + TH1F* hData = (TH1F*)file->Get("hV0MUltraFine"); + TH1F* hGlauberParameters = (TH1F*)file->Get("hGlauberParameters"); + TH1F* hGlauberFitRange = (TH1F*)file->Get("hGlauberFitRange"); hData->SetName("hData"); - TH1F *hStitched = (TH1F*) hData->Clone("hStitched"); - TH1F *hFit = (TH1F*) file->Get("hGlauber"); - - TCanvas *c1 = new TCanvas("c1", "", 800, 600); + TH1F* hStitched = (TH1F*)hData->Clone("hStitched"); + TH1F* hFit = (TH1F*)file->Get("hGlauber"); + + TCanvas* c1 = new TCanvas("c1", "", 800, 600); c1->SetLeftMargin(0.17); c1->SetBottomMargin(0.17); c1->SetRightMargin(0.15); c1->SetTopMargin(0.05); - c1->SetTicks(1,1); + c1->SetTicks(1, 1); c1->SetLogz(); c1->SetFrameFillStyle(0); c1->SetFillStyle(0); - - - cout<<"Data bin width: "<GetBinWidth(1)<GetBinWidth(1)<GetBinWidth(1) << endl; + cout << "Fit bin width: " << hFit->GetBinWidth(1) << endl; + cout << "Match range to use: " << matchRange << endl; + //____________________________________________ - double anchorPointFraction = anchorPointPercentage/100.f; + double anchorPointFraction = anchorPointPercentage / 100.f; double anchorPoint = -1; // the anchor point value in raw - + //____________________________________________ // doing partial integration up to certain point for finding anchor point bin - for(int ii=1; iiGetNbinsX()+1; ii++){ + for (int ii = 1; ii < hData->GetNbinsX() + 1; ii++) { // renormalize data curve - int bin1 = ii+1; - int bin2 = hData->FindBin( hData->GetBinLowEdge(ii+1) + matchRange + 1e-3 ); - double matchRangeData = hData -> Integral( bin1, bin2); - double matchRangeFit = hFit -> Integral( bin1, bin2); - + int bin1 = ii + 1; + int bin2 = hData->FindBin(hData->GetBinLowEdge(ii + 1) + matchRange + 1e-3); + double matchRangeData = hData->Integral(bin1, bin2); + double matchRangeFit = hFit->Integral(bin1, bin2); + // rescale fit to match in the vicinity of the region we're at - hFit->Scale(matchRangeData/matchRangeFit); - - double integralFit = hFit->Integral(1,ii); - double integralData = hData->Integral(ii+1,hData->GetNbinsX()+1); - double integralAll = integralFit+integralData; - - cout<<"at bin #"<GetBinLowEdge(ii+1)<<" fraction above this value is: "<GetBinLowEdge(ii+1); - - if(integralData/integralAllScale(matchRangeData / matchRangeFit); + + double integralFit = hFit->Integral(1, ii); + double integralData = hData->Integral(ii + 1, hData->GetNbinsX() + 1); + double integralAll = integralFit + integralData; + + cout << "at bin #" << ii << ", integrated up to " << hData->GetBinLowEdge(ii + 1) << " fraction above this value is: " << integralData / integralAll << endl; + anchorPoint = hData->GetBinLowEdge(ii + 1); + + if (integralData / integralAll < anchorPointFraction) + break; } - + //____________________________________________ - for(int ii=1; iiGetNbinsX()+1; ii++){ + for (int ii = 1; ii < hData->GetNbinsX() + 1; ii++) { // renormalize data curve - if( hData->GetBinCenter(ii) < anchorPoint) hStitched->SetBinContent(ii, hFit->GetBinContent(ii)); + if (hData->GetBinCenter(ii) < anchorPoint) + hStitched->SetBinContent(ii, hFit->GetBinContent(ii)); } - - - cout<<"Anchor point determined to be: "<SetLineColor(kRed); hStitched->SetLineColor(kBlue); - + hData->GetYaxis()->SetTitleSize(0.055); hData->GetXaxis()->SetTitleSize(0.055); hData->GetYaxis()->SetLabelSize(0.04); @@ -71,75 +72,74 @@ void runCalibration( TString lInputFileName = "results/AR_544122_glauberNBD_ance hData->Draw("hist"); hFit->Draw("hist same"); hStitched->Draw("hist same"); - - //All fine, let's try the calibrator - multCalibrator *lCalib = new multCalibrator("lCalib"); + + // All fine, let's try the calibrator + multCalibrator* lCalib = new multCalibrator("lCalib"); lCalib->SetAnchorPointPercentage(100.0f); lCalib->SetAnchorPointRaw(-1e-6); - - //Set standard Pb-Pb boundaries + + // Set standard Pb-Pb boundaries lCalib->SetStandardOnePercentBoundaries(); - + TString calibFileName = lInputFileName.Data(); calibFileName.ReplaceAll("glauberNBD", "calibration"); - TFile *fileCalib = new TFile(calibFileName.Data(), "RECREATE"); - - TH1F *hCalib = lCalib->GetCalibrationHistogram(hStitched, "hCalib"); - - TCanvas *c2 = new TCanvas("c2", "", 800, 600); + TFile* fileCalib = new TFile(calibFileName.Data(), "RECREATE"); + + TH1F* hCalib = lCalib->GetCalibrationHistogram(hStitched, "hCalib"); + + TCanvas* c2 = new TCanvas("c2", "", 800, 600); c2->SetLeftMargin(0.17); c2->SetBottomMargin(0.17); c2->SetRightMargin(0.15); c2->SetTopMargin(0.05); - c2->SetTicks(1,1); + c2->SetTicks(1, 1); // c2->SetLogz(); c2->SetFrameFillStyle(0); c2->SetFillStyle(0); - + hCalib->GetYaxis()->SetTitleSize(0.055); hCalib->GetXaxis()->SetTitleSize(0.055); hCalib->GetYaxis()->SetLabelSize(0.04); hCalib->GetXaxis()->SetLabelSize(0.04); hCalib->SetTitle(""); hCalib->Draw(); - - + fileCalib->cd(); hData->Write(); hCalib->Write(); hStitched->Write(); hFit->Write(); - - if(doNpartNcoll){ - cout<<"Will now attempt to calculate % -> Np, Nc map..."< Np, Nc map..." << endl; + + TProfile* hProfileNpart = new TProfile("hProfileNpart", "", 100, 0, 100); + TProfile* hProfileNcoll = new TProfile("hProfileNcoll", "", 100, 0, 100); + TH2F* h2dNpart = new TH2F("h2dNpart", "", 100, 0, 100, 500, -0.5f, 499.5f); + TH2F* h2dNcoll = new TH2F("h2dNcoll", "", 100, 0, 100, 3000, -0.5f, 2999.5); + // Replay - multGlauberNBDFitter *g = new multGlauberNBDFitter("lglau"); - TF1 *fitfunc = g->GetGlauberNBD(); - - //Step 1: open the (Npart, Ncoll) pair information, provide - TFile *fbasefile = new TFile("basehistos.root","READ"); - TH2D *hNpNc = (TH2D*) fbasefile->Get("hNpNc"); + multGlauberNBDFitter* g = new multGlauberNBDFitter("lglau"); + TF1* fitfunc = g->GetGlauberNBD(); + + // Step 1: open the (Npart, Ncoll) pair information, provide + TFile* fbasefile = new TFile("basehistos.root", "READ"); + TH2D* hNpNc = (TH2D*)fbasefile->Get("hNpNc"); g->SetNpartNcollCorrelation(hNpNc); g->InitializeNpNc(); - + fitfunc->SetParameter(0, hGlauberParameters->GetBinContent(1)); fitfunc->SetParameter(1, hGlauberParameters->GetBinContent(2)); fitfunc->SetParameter(2, hGlauberParameters->GetBinContent(3)); fitfunc->SetParameter(3, hGlauberParameters->GetBinContent(4)); fitfunc->SetParameter(4, hGlauberParameters->GetBinContent(5)); - - Double_t lMax = hData->GetBinLowEdge( hData->GetNbinsX() + 1); - + + Double_t lMax = hData->GetBinLowEdge(hData->GetNbinsX() + 1); + // uncomment if Np Nc needed -> Warning, slow! - g->CalculateAvNpNc( hProfileNpart, hProfileNcoll, h2dNpart, h2dNcoll, hCalib, 0, lMax ); - + g->CalculateAvNpNc(hProfileNpart, hProfileNcoll, h2dNpart, h2dNcoll, hCalib, 0, lMax); + hProfileNpart->Write(); hProfileNcoll->Write(); h2dNpart->Write(); diff --git a/Common/Tools/Multiplicity/macros/runGlauberFit.C b/Common/Tools/Multiplicity/macros/runGlauberFit.C index df1041da9fe..0d20a15e637 100755 --- a/Common/Tools/Multiplicity/macros/runGlauberFit.C +++ b/Common/Tools/Multiplicity/macros/runGlauberFit.C @@ -1,57 +1,60 @@ #include "multCalibrator.h" #include "multGlauberNBDFitter.h" + +#include "TCanvas.h" +#include "TDirectory.h" +#include "TF1.h" #include "TFile.h" +#include "TGraph.h" #include "TH1F.h" #include "TH2F.h" -#include "TF1.h" -#include "TCanvas.h" -#include "TLegend.h" #include "TLatex.h" +#include "TLegend.h" #include "TLine.h" -#include "TDirectory.h" +#include "TStopwatch.h" #include "TStyle.h" #include "TSystem.h" -#include "TGraph.h" #include "TTree.h" -#include "TStopwatch.h" //________________________________________________________________ -Double_t FastIntegrate(TF1 *f1, Double_t a, Double_t b, Int_t n = 5){ - //Do fast integration with N sampling points +Double_t FastIntegrate(TF1* f1, Double_t a, Double_t b, Int_t n = 5) +{ + // Do fast integration with N sampling points const Int_t nc = n; Double_t x[nc], y[nc]; - Double_t lWidth = (b-a)/((double)(n-1)); - for(Int_t ii=0; iiEval( x[ii] ); + Double_t lWidth = (b - a) / ((double)(n - 1)); + for (Int_t ii = 0; ii < n; ii++) { + x[ii] = a + ((double)(ii)) * lWidth; + y[ii] = f1->Eval(x[ii]); } - //Now go via trapezoids, please (this probably has a name) + // Now go via trapezoids, please (this probably has a name) Double_t lIntegral = 0; - for(Int_t ii=0; iiGetNbinsX(); - Double_t lCountDesired = lPercentile * histo->GetEntries()/100; + Double_t lCountDesired = lPercentile * histo->GetEntries() / 100; Long_t lCount = 0; - for(Long_t ibin=1;ibinGetBinContent(ibin); - if(lCount >= lCountDesired){ - //Found bin I am looking for! + if (lCount >= lCountDesired) { + // Found bin I am looking for! Double_t lWidth = histo->GetBinWidth(ibin); - Double_t lLeftPercentile = 100.*(lCount - histo->GetBinContent(ibin))/histo->GetEntries(); - Double_t lRightPercentile = 100.*lCount / histo->GetEntries(); - - Double_t lProportion = (lPercentile - lLeftPercentile)/(lRightPercentile-lLeftPercentile); - - lReturnValue = histo->GetBinLowEdge(ibin) + lProportion*lWidth; + Double_t lLeftPercentile = 100. * (lCount - histo->GetBinContent(ibin)) / histo->GetEntries(); + Double_t lRightPercentile = 100. * lCount / histo->GetEntries(); + + Double_t lProportion = (lPercentile - lLeftPercentile) / (lRightPercentile - lLeftPercentile); + + lReturnValue = histo->GetBinLowEdge(ibin) + lProportion * lWidth; break; } } @@ -68,40 +71,45 @@ Double_t GetBoundaryForPercentile( TH1 *histo, Double_t lPercentileRequested ) { // --- free f: keep f value free (default Pb-Pb: fixed at 0.8) // --- f value: the value to use for fixed f -int runGlauberFit(TString lInputFileName = "AnalysisResultsLHC24ar.root", TString histogramName = "hFT0C_BCs", int ancestorMode = 2, Bool_t lFreek = kFALSE, Bool_t use_dMu_dNanc = kFALSE, Bool_t lFreef = kFALSE, Float_t lfvalue = 0.800, TString outputFile = "output.root") { +int runGlauberFit(TString lInputFileName = "AnalysisResultsLHC24ar.root", TString histogramName = "hFT0C_BCs", int ancestorMode = 2, Bool_t lFreek = kFALSE, Bool_t use_dMu_dNanc = kFALSE, Bool_t lFreef = kFALSE, Float_t lfvalue = 0.800, TString outputFile = "output.root") +{ gStyle->SetLineScalePS(1); gStyle->SetOptStat(0); - - cout<<"Starting!"< Get(Form("centrality-study/%s", histogramName.Data())); + cout << "Starting!" << endl; + TFile* file = new TFile(lInputFileName.Data(), "READ"); + if (!file) + cout << "Problem with file!" << endl; + TH1F* hV0Mfine = 0x0; + + hV0Mfine = (TH1F*)file->Get(Form("centrality-study/%s", histogramName.Data())); // disregard bin zero - cout<<"Received bin zero content: "<< hV0Mfine ->GetBinContent(0)<<", will set to zero..."<SetBinContent(0, 0); + cout << "Received bin zero content: " << hV0Mfine->GetBinContent(0) << ", will set to zero..." << endl; + hV0Mfine->SetBinContent(0, 0); + + if (!hV0Mfine) + cout << "Problem with histogram!" << endl; - if(!hV0Mfine) cout<<"Problem with histogram!"<GetEntries()<GetNbinsX()<GetBinLowEdge(hV0Mfine->GetNbinsX()+1)<GetEntries() << endl; + cout << "NbinsX: " << hV0Mfine->GetNbinsX() << endl; + cout << "MaxX: " << hV0Mfine->GetBinLowEdge(hV0Mfine->GetNbinsX() + 1) << endl; - cout<<"Creating output file..."<Clone("hV0M"); - TH1F *hV0MUltraFine = (TH1F*) hV0Mfine->Clone("hV0MUltraFine"); + TFile* fOutput = new TFile(outputFile.Data(), "RECREATE"); + + TH1F* hV0M = (TH1F*)hV0Mfine->Clone("hV0M"); + TH1F* hV0MUltraFine = (TH1F*)hV0Mfine->Clone("hV0MUltraFine"); hV0M->SetName("hV0M"); hV0M->SetTitle(""); @@ -109,54 +117,56 @@ int runGlauberFit(TString lInputFileName = "AnalysisResultsLHC24ar.root", TStrin // maximum fit range estimate (avoid tails) // may need adjusting Double_t lFitRangeMax = GetBoundaryForPercentile(hV0Mfine, 0.008); - cout<<"Fit range max estimated from histogram: "<Rebin(rebinFactor); //____________________________________________ // simple plots for inspection - TCanvas *c1 = new TCanvas("c1", "", 1300,900); + TCanvas* c1 = new TCanvas("c1", "", 1300, 900); c1->SetFrameFillStyle(0); c1->SetFillStyle(0); - c1->Divide(1,2); + c1->Divide(1, 2); c1->cd(1)->SetFrameFillStyle(0); c1->cd(1)->SetFillStyle(0); c1->cd(2)->SetFrameFillStyle(0); c1->cd(2)->SetFillStyle(0); - + c1->cd(1); c1->cd(1)->SetLogy(); - c1->cd(1)->SetTicks(1,1); - c1->cd(1)->SetPad(0,0.5,1,1); - c1->cd(2)->SetPad(0,0.0,1,.5); - + c1->cd(1)->SetTicks(1, 1); + c1->cd(1)->SetPad(0, 0.5, 1, 1); + c1->cd(2)->SetPad(0, 0.0, 1, .5); + c1->cd(1)->SetBottomMargin(0.001); c1->cd(1)->SetRightMargin(0.25); c1->cd(1)->SetTopMargin(0.02); c1->cd(1)->SetLeftMargin(0.07); - + c1->cd(2)->SetBottomMargin(0.14); c1->cd(2)->SetRightMargin(0.25); c1->cd(2)->SetTopMargin(0.001); c1->cd(2)->SetLeftMargin(0.07); - c1->cd(2)->SetTicks(1,1); + c1->cd(2)->SetTicks(1, 1); c1->cd(1); - - hV0M->GetXaxis()->SetRangeUser(0,lFitRangeMax); - hV0M->GetYaxis()->SetRangeUser(0.25,hV0M->GetMaximum()*3); + + hV0M->GetXaxis()->SetRangeUser(0, lFitRangeMax); + hV0M->GetYaxis()->SetRangeUser(0.25, hV0M->GetMaximum() * 3); hV0M->SetLineColor(kBlack); hV0M->SetMarkerStyle(20); hV0M->SetMarkerColor(kBlack); @@ -171,112 +181,112 @@ int runGlauberFit(TString lInputFileName = "AnalysisResultsLHC24ar.root", TStrin hV0M->GetYaxis()->SetTickLength(0.015); hV0M->SetStats(0); hV0M->Draw("E"); - - //Stand back! Imma gonna do GLAUBER FITTIN' - multGlauberNBDFitter *g = new multGlauberNBDFitter("lglau"); + + // Stand back! Imma gonna do GLAUBER FITTIN' + multGlauberNBDFitter* g = new multGlauberNBDFitter("lglau"); g->SetAncestorMode(ancestorMode); - - //Step 1: open the (Npart, Ncoll) pair information, provide - TFile *fbasefile = new TFile("basehistos.root","READ"); - TH2D *hNpNc = (TH2D*) fbasefile->Get("hNpNc"); + + // Step 1: open the (Npart, Ncoll) pair information, provide + TFile* fbasefile = new TFile("basehistos.root", "READ"); + TH2D* hNpNc = (TH2D*)fbasefile->Get("hNpNc"); // return to proper scope - fOutput->cd(); + fOutput->cd(); g->SetNpartNcollCorrelation(hNpNc); g->SetInputV0M(hV0M); - g->SetFitRange(lFitRange,lFitRangeMax); - //Step 3: go for it ... + g->SetFitRange(lFitRange, lFitRangeMax); + // Step 3: go for it ... g->SetNorm(1.53527e+08); TString lString = "REM0"; g->SetFitOptions(lString.Data()); g->SetFitNpx(100000); - - TF1 *fitfunc = g->GetGlauberNBD(); + + TF1* fitfunc = g->GetGlauberNBD(); //____________________________________________ // // set initial fit parameters here // may require manual tuning depending on data! - - Double_t guessedMu = lFitRangeMax/53968.4 * 0.175*3.53971e+02; - cout<<"Guessed GlauberNBD mu value: "<SetParameter(0,guessedMu); //mu value - fitfunc->SetParLimits(0,0.25*guessedMu,guessedMu*2); - - if(!lFreek){ - fitfunc->FixParameter(1,1.5); //k value - }else{ - fitfunc->SetParLimits(1,0.01,35); - fitfunc->SetParameter(1,1.5); + + Double_t guessedMu = lFitRangeMax / 53968.4 * 0.175 * 3.53971e+02; + cout << "Guessed GlauberNBD mu value: " << guessedMu << endl; + + fitfunc->SetParameter(0, guessedMu); // mu value + fitfunc->SetParLimits(0, 0.25 * guessedMu, guessedMu * 2); + + if (!lFreek) { + fitfunc->FixParameter(1, 1.5); // k value + } else { + fitfunc->SetParLimits(1, 0.01, 35); + fitfunc->SetParameter(1, 1.5); } - if(!lFreef){ - fitfunc->FixParameter(2,lfvalue); // f value - }else{ - fitfunc->SetParLimits(2,0.55,0.97); + if (!lFreef) { + fitfunc->FixParameter(2, lfvalue); // f value + } else { + fitfunc->SetParLimits(2, 0.55, 0.97); fitfunc->SetParameter(2, 0.800); } - fitfunc->SetParLimits(3,0.1e+5, 800e+10); // normalization + fitfunc->SetParLimits(3, 0.1e+5, 800e+10); // normalization fitfunc->SetParameter(3, 0.80832e+08); - - //dMu/dNanc - fitfunc->SetParameter(4,0); - if( !use_dMu_dNanc ){ - fitfunc->FixParameter(4,0); + + // dMu/dNanc + fitfunc->SetParameter(4, 0); + if (!use_dMu_dNanc) { + fitfunc->FixParameter(4, 0); } - - //fitfunc->SetParameter(4,-1.15443e-01); - //fitfunc->SetParameter(4,-4.55957e-02); - - //dk/dNanc - fitfunc->FixParameter(5,0); - //fitfunc->SetParameter(5,1.63590e-03); - - //d2Mu/dNanc2 - fitfunc->FixParameter(6,0.0); - //fitfunc->SetParameter(6,4.02271e-05); - - fitfunc->FixParameter(7,0.0); - //fitfunc->SetParameter(7,-1.24349e-06); - + + // fitfunc->SetParameter(4,-1.15443e-01); + // fitfunc->SetParameter(4,-4.55957e-02); + + // dk/dNanc + fitfunc->FixParameter(5, 0); + // fitfunc->SetParameter(5,1.63590e-03); + + // d2Mu/dNanc2 + fitfunc->FixParameter(6, 0.0); + // fitfunc->SetParameter(6,4.02271e-05); + + fitfunc->FixParameter(7, 0.0); + // fitfunc->SetParameter(7,-1.24349e-06); + //____________________________________________ // handle internals for fitting: needs to // be done before the fit is attempted! g->InitializeNpNc(); g->InitAncestor(); - - cout<<"WILL NOW ATTEMPT GLAUBER FIT"<DoFit(); Int_t lAttempts = 1; - while(lAttempts < 10 && lFitStatus==0){ - // insist on fitting until it works - cout<<"Attempting fit again ("<DoFit(); } - cout<<"Final fit status: "<SetOptStat(0); - //Do a ratio plot - TH1D *hGlauber = (TH1D*) hV0MUltraFine->Clone("hGlauber"); - TH1D *hRatio = (TH1D*) hV0MUltraFine->Clone("hRatio"); + // Do a ratio plot + TH1D* hGlauber = (TH1D*)hV0MUltraFine->Clone("hGlauber"); + TH1D* hRatio = (TH1D*)hV0MUltraFine->Clone("hRatio"); hGlauber->Reset(); - + c1->cd(1); fitfunc->SetLineColor(kRed); fitfunc->SetLineWidth(2); fitfunc->Draw("same"); - - cout<<"Calculating glauber function histogram with the same binning as data input... please wait..."<GetNbinsX()+1; ii++){ - Double_t lFuncVal = FastIntegrate( fitfunc, hGlauber->GetBinLowEdge(ii), hGlauber->GetBinLowEdge(ii+1), 4); + + cout << "Calculating glauber function histogram with the same binning as data input... please wait..." << endl; + for (Int_t ii = 1; ii < hGlauber->GetNbinsX() + 1; ii++) { + Double_t lFuncVal = FastIntegrate(fitfunc, hGlauber->GetBinLowEdge(ii), hGlauber->GetBinLowEdge(ii + 1), 4); hGlauber->SetBinContent(ii, lFuncVal); - if(ii%500==0){ - cout<<"At integration #"<GetNbinsX()+1<<"..."<GetNbinsX() + 1 << "..." << endl; } } - cout<<"Glauber function evaluated. Should go quickly now."<cd(2); Float_t lLoRangeRatio = 0.35; Float_t lHiRangeRatio = 1.65; @@ -288,65 +298,65 @@ int runGlauberFit(TString lInputFileName = "AnalysisResultsLHC24ar.root", TStrin hRatio->GetXaxis()->SetTitleSize(0.055); hRatio->GetYaxis()->SetLabelSize(0.045); hRatio->GetXaxis()->SetLabelSize(0.045); - hRatio->GetYaxis()->SetRangeUser(lLoRangeRatio,lHiRangeRatio); - hRatio->GetXaxis()->SetRangeUser(0,lFitRangeMax); + hRatio->GetYaxis()->SetRangeUser(lLoRangeRatio, lHiRangeRatio); + hRatio->GetXaxis()->SetRangeUser(0, lFitRangeMax); hRatio->SetMarkerStyle(20); - hRatio->SetMarkerColor(kGray+2); - hRatio->SetLineColor(kGray+2); - //hRatio->SetMarkerSize(1.0); + hRatio->SetMarkerColor(kGray + 2); + hRatio->SetLineColor(kGray + 2); + // hRatio->SetMarkerSize(1.0); hRatio->SetMarkerSize(.7); hRatio->SetStats(0); - //hRatioWide->SetStats(0); - + // hRatioWide->SetStats(0); + hRatio->Draw("hist"); - - TLine *line = new TLine(0,1,lFitRangeMax,1); + + TLine* line = new TLine(0, 1, lFitRangeMax, 1); line->SetLineStyle(7); - line->SetLineColor(kGray+1); + line->SetLineColor(kGray + 1); line->Draw(); - - TLine *lFitRangeLine = new TLine( lFitRange, lLoRangeRatio, lFitRange, 0.9) ; + + TLine* lFitRangeLine = new TLine(lFitRange, lLoRangeRatio, lFitRange, 0.9); lFitRangeLine->SetLineColor(kBlue); lFitRangeLine->SetLineWidth(1); lFitRangeLine->SetLineStyle(2); lFitRangeLine->Draw(); - - TH1D *hRatioGrayed = (TH1D*) hRatio->Clone("hRatioGrayed"); - hRatioGrayed->SetMarkerColor(kGray+2); - hRatioGrayed->SetLineColor(kGray+2); + + TH1D* hRatioGrayed = (TH1D*)hRatio->Clone("hRatioGrayed"); + hRatioGrayed->SetMarkerColor(kGray + 2); + hRatioGrayed->SetLineColor(kGray + 2); hRatioGrayed->Draw("same"); - + hRatio->SetLineWidth(1); hRatio->Draw("same hist"); - //hRatioWide->Draw("same"); - + // hRatioWide->Draw("same"); + c1->cd(1); - TLatex *lat = new TLatex(); + TLatex* lat = new TLatex(); lat->SetNDC(); Float_t lPosText = 0.76; Float_t lYShift = 0.25; lat->SetTextSize(0.042); // save the glauber parameters explicitly - TH1D *hGlauberParameters = new TH1D("hGlauberParameters", "", 10,0,10); - TH1D *hGlauberFitRange = new TH1D("hGlauberFitRange", "", 10,0,10); - - //fitfunc - hGlauberParameters -> SetBinContent( 1, fitfunc -> GetParameter(0)); - hGlauberParameters -> SetBinContent( 2, fitfunc -> GetParameter(1)); - hGlauberParameters -> SetBinContent( 3, fitfunc -> GetParameter(2)); - hGlauberParameters -> SetBinContent( 4, fitfunc -> GetParameter(3)); - hGlauberParameters -> SetBinContent( 5, fitfunc -> GetParameter(4)); - hGlauberParameters -> Write(); - + TH1D* hGlauberParameters = new TH1D("hGlauberParameters", "", 10, 0, 10); + TH1D* hGlauberFitRange = new TH1D("hGlauberFitRange", "", 10, 0, 10); + + // fitfunc + hGlauberParameters->SetBinContent(1, fitfunc->GetParameter(0)); + hGlauberParameters->SetBinContent(2, fitfunc->GetParameter(1)); + hGlauberParameters->SetBinContent(3, fitfunc->GetParameter(2)); + hGlauberParameters->SetBinContent(4, fitfunc->GetParameter(3)); + hGlauberParameters->SetBinContent(5, fitfunc->GetParameter(4)); + hGlauberParameters->Write(); + Double_t lLoRangeGlauber, lHiRangeGlauber; fitfunc->GetRange(lLoRangeGlauber, lHiRangeGlauber); - hGlauberFitRange->SetBinContent(1, lLoRangeGlauber); + hGlauberFitRange->SetBinContent(1, lLoRangeGlauber); hGlauberFitRange->SetBinContent(2, lHiRangeGlauber); hGlauberFitRange->Write(); hRatio->Write(); fOutput->Write(); - + return 0; } diff --git a/Common/Tools/Multiplicity/macros/saveCorrelation.C b/Common/Tools/Multiplicity/macros/saveCorrelation.C index 9e9e5eb5776..9f0a3b35c2b 100644 --- a/Common/Tools/Multiplicity/macros/saveCorrelation.C +++ b/Common/Tools/Multiplicity/macros/saveCorrelation.C @@ -11,45 +11,46 @@ // // This file provides a simple macro to convert TGlauberMC output tuples // into the 2D correlation histogram necessary for the ALICE machinery -// that performs Glauber + NBD fits. +// that performs Glauber + NBD fits. + +void saveCorrelation(TString filename = "gmc-PbPb-snn68.21-md0.40-nd-1.0-rc1-smax99.0.root", TString outputFile = "basehistos.root") +{ + TFile* fin = new TFile(filename.Data(), "READ"); + TNtuple* ntup = (TNtuple*)fin->Get("nt_Pb_Pb"); -void saveCorrelation(TString filename = "gmc-PbPb-snn68.21-md0.40-nd-1.0-rc1-smax99.0.root", TString outputFile = "basehistos.root"){ - TFile *fin = new TFile(filename.Data(),"READ"); - TNtuple *ntup = (TNtuple*) fin->Get("nt_Pb_Pb"); - // try other Pb nuclear profiles in case "Pb" - "Pb" not found - if(!ntup){ - ntup = (TNtuple*) fin->Get("nt_Pbpn_Pbpn"); + if (!ntup) { + ntup = (TNtuple*)fin->Get("nt_Pbpn_Pbpn"); } - if(!ntup){ - ntup = (TNtuple*) fin->Get("nt_PbpnVar1_PbpnVar1"); + if (!ntup) { + ntup = (TNtuple*)fin->Get("nt_PbpnVar1_PbpnVar1"); } - if(!ntup){ - ntup = (TNtuple*) fin->Get("nt_PbpnVar2_PbpnVar2"); + if (!ntup) { + ntup = (TNtuple*)fin->Get("nt_PbpnVar2_PbpnVar2"); } - if(!ntup){ - ntup = (TNtuple*) fin->Get("nt_PbpnVar3_PbpnVar3"); + if (!ntup) { + ntup = (TNtuple*)fin->Get("nt_PbpnVar3_PbpnVar3"); } - if(!ntup){ - ntup = (TNtuple*) fin->Get("nt_PbpnVar4_PbpnVar4"); + if (!ntup) { + ntup = (TNtuple*)fin->Get("nt_PbpnVar4_PbpnVar4"); } - - if(!ntup){ - cout<<"No tree found!"<GetEntries()<GetEntries() << endl; + + TFile* fout = new TFile(outputFile.Data(), "RECREATE"); + // 2D correlation plot necessary for Glauber + NBD fitting // The provided range should be enough for Pb-Pb - TH2D *hNpNc = new TH2D("hNpNc","",500,-0.5,499.5,2500,-0.5,2499.5); - + TH2D* hNpNc = new TH2D("hNpNc", "", 500, -0.5, 499.5, 2500, -0.5, 2499.5); + // let's draw this on screen for inspection - TCanvas *c1 = new TCanvas("c1", "", 800, 600); - c1->SetTicks(1,1); + TCanvas* c1 = new TCanvas("c1", "", 800, 600); + c1->SetTicks(1, 1); c1->SetLogz(); ntup->Draw("Ncoll:Npart>>hNpNc", "", "colz"); fout->Write(); From ef2cc90354598144eef930d3185940ed0b2f5f34 Mon Sep 17 00:00:00 2001 From: ddobrigk Date: Sat, 8 Nov 2025 17:51:58 +0100 Subject: [PATCH 3/9] Add copyright notice --- .../Multiplicity/macros/runCalibration.C | 22 +++++++++++-- .../Tools/Multiplicity/macros/runGlauberFit.C | 31 +++++++++++++------ .../Multiplicity/macros/saveCorrelation.C | 3 ++ 3 files changed, 43 insertions(+), 13 deletions(-) diff --git a/Common/Tools/Multiplicity/macros/runCalibration.C b/Common/Tools/Multiplicity/macros/runCalibration.C index a8611c9a5ec..c93584ddf2b 100644 --- a/Common/Tools/Multiplicity/macros/runCalibration.C +++ b/Common/Tools/Multiplicity/macros/runCalibration.C @@ -1,6 +1,22 @@ -void runCalibration(TString lInputFileName = "results/AR_544122_glauberNBD_ancestorMode2_hFT0C_BCs.root", double anchorPointPercentage = 90.0, double matchRange = 200.0, bool doNpartNcoll = false) -{ - TFile* file = new TFile(lInputFileName.Data(), "READ"); +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +// + +/// @brief function to calibrate centrality +/// @param lInputFileName name of input file. +/// @param anchorPointPercentage anchor point percentage to use +/// @param matchRange width of region in which data/glauber matching is to be done in rolling anchoring test +/// @param doNpartNcoll wether or not to attempt calculating Npart, Ncoll in centrality bins +void runCalibration( TString lInputFileName = "results/AR_544122_glauberNBD_ancestorMode2_hFT0C_BCs.root", double anchorPointPercentage = 90.0, double matchRange = 200.0, bool doNpartNcoll = false){ + TFile *file = new TFile(lInputFileName.Data(), "READ"); file->ls(); TH1F* hData = (TH1F*)file->Get("hV0MUltraFine"); diff --git a/Common/Tools/Multiplicity/macros/runGlauberFit.C b/Common/Tools/Multiplicity/macros/runGlauberFit.C index 0d20a15e637..e89b77f2400 100755 --- a/Common/Tools/Multiplicity/macros/runGlauberFit.C +++ b/Common/Tools/Multiplicity/macros/runGlauberFit.C @@ -1,3 +1,15 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +// + #include "multCalibrator.h" #include "multGlauberNBDFitter.h" @@ -61,16 +73,15 @@ Double_t GetBoundaryForPercentile(TH1* histo, Double_t lPercentileRequested) return lReturnValue; } -// master fit function -// parameters: -// --- input file name (from the grid, typically centrality-studies task) -// --- histogram name: histogram to use within the input file -// --- ancestor mode: 0: truncation, 1: rounding, 2: effective / non-integer (default: 2) -// --- free k: keep k value free (default Pb-Pb: fixed at 1.5) -// --- use dMu/dNanc: allow for a varying production of Nch vs ancestor if Nancestor is large. Models detector saturation in effective manner. -// --- free f: keep f value free (default Pb-Pb: fixed at 0.8) -// --- f value: the value to use for fixed f - +/// @brief master glauber fit function +/// @param lInputFileName input file name (from the grid, typically centrality-studies task) +/// @param histogramName histogram name: histogram to use within the input file +/// @param ancestorMode ancestor mode: 0: truncation, 1: rounding, 2: effective / non-integer (default: 2) +/// @param lFreek free k: keep k value free (default Pb-Pb: fixed at 1.5) +/// @param use_dMu_dNanc use dMu/dNanc: allow for a varying production of Nch vs ancestor if Nancestor is large. Models detector saturation in effective manner. +/// @param lFreef free f: keep f value free (default Pb-Pb: fixed at 0.8) +/// @param lfvalue f value: the value to use for fixed f +/// @param outputFile name of output file int runGlauberFit(TString lInputFileName = "AnalysisResultsLHC24ar.root", TString histogramName = "hFT0C_BCs", int ancestorMode = 2, Bool_t lFreek = kFALSE, Bool_t use_dMu_dNanc = kFALSE, Bool_t lFreef = kFALSE, Float_t lfvalue = 0.800, TString outputFile = "output.root") { gStyle->SetLineScalePS(1); diff --git a/Common/Tools/Multiplicity/macros/saveCorrelation.C b/Common/Tools/Multiplicity/macros/saveCorrelation.C index 9f0a3b35c2b..a011c65d32a 100644 --- a/Common/Tools/Multiplicity/macros/saveCorrelation.C +++ b/Common/Tools/Multiplicity/macros/saveCorrelation.C @@ -13,6 +13,9 @@ // into the 2D correlation histogram necessary for the ALICE machinery // that performs Glauber + NBD fits. +/// @brief function to save Npart x Ncoll correlation to file for glauber fits +/// @param filename input TGlauberMC ntuple file +/// @param outputFile output file for Npart x Ncoll correlation TH2D void saveCorrelation(TString filename = "gmc-PbPb-snn68.21-md0.40-nd-1.0-rc1-smax99.0.root", TString outputFile = "basehistos.root") { TFile* fin = new TFile(filename.Data(), "READ"); From 187e22dd153b5830cbf6cba51f1211df4f98d2f4 Mon Sep 17 00:00:00 2001 From: ALICE Builder Date: Sat, 8 Nov 2025 17:58:00 +0100 Subject: [PATCH 4/9] Please consider the following formatting changes (#496) --- Common/Tools/Multiplicity/macros/runCalibration.C | 11 ++++++----- Common/Tools/Multiplicity/macros/saveCorrelation.C | 4 ++-- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/Common/Tools/Multiplicity/macros/runCalibration.C b/Common/Tools/Multiplicity/macros/runCalibration.C index c93584ddf2b..0e3fbfdfec8 100644 --- a/Common/Tools/Multiplicity/macros/runCalibration.C +++ b/Common/Tools/Multiplicity/macros/runCalibration.C @@ -10,13 +10,14 @@ // or submit itself to any jurisdiction. // -/// @brief function to calibrate centrality -/// @param lInputFileName name of input file. -/// @param anchorPointPercentage anchor point percentage to use +/// @brief function to calibrate centrality +/// @param lInputFileName name of input file. +/// @param anchorPointPercentage anchor point percentage to use /// @param matchRange width of region in which data/glauber matching is to be done in rolling anchoring test /// @param doNpartNcoll wether or not to attempt calculating Npart, Ncoll in centrality bins -void runCalibration( TString lInputFileName = "results/AR_544122_glauberNBD_ancestorMode2_hFT0C_BCs.root", double anchorPointPercentage = 90.0, double matchRange = 200.0, bool doNpartNcoll = false){ - TFile *file = new TFile(lInputFileName.Data(), "READ"); +void runCalibration(TString lInputFileName = "results/AR_544122_glauberNBD_ancestorMode2_hFT0C_BCs.root", double anchorPointPercentage = 90.0, double matchRange = 200.0, bool doNpartNcoll = false) +{ + TFile* file = new TFile(lInputFileName.Data(), "READ"); file->ls(); TH1F* hData = (TH1F*)file->Get("hV0MUltraFine"); diff --git a/Common/Tools/Multiplicity/macros/saveCorrelation.C b/Common/Tools/Multiplicity/macros/saveCorrelation.C index a011c65d32a..160d0177c69 100644 --- a/Common/Tools/Multiplicity/macros/saveCorrelation.C +++ b/Common/Tools/Multiplicity/macros/saveCorrelation.C @@ -13,9 +13,9 @@ // into the 2D correlation histogram necessary for the ALICE machinery // that performs Glauber + NBD fits. -/// @brief function to save Npart x Ncoll correlation to file for glauber fits +/// @brief function to save Npart x Ncoll correlation to file for glauber fits /// @param filename input TGlauberMC ntuple file -/// @param outputFile output file for Npart x Ncoll correlation TH2D +/// @param outputFile output file for Npart x Ncoll correlation TH2D void saveCorrelation(TString filename = "gmc-PbPb-snn68.21-md0.40-nd-1.0-rc1-smax99.0.root", TString outputFile = "basehistos.root") { TFile* fin = new TFile(filename.Data(), "READ"); From 8a537e8f3aeb94302fbadf6d59dcea5d587b0c63 Mon Sep 17 00:00:00 2001 From: ddobrigk Date: Sat, 8 Nov 2025 18:01:13 +0100 Subject: [PATCH 5/9] Remove trailing whitespace --- Common/Tools/Multiplicity/README.md | 22 +++--- Common/Tools/Multiplicity/macros/README.md | 84 +++++++++++----------- 2 files changed, 53 insertions(+), 53 deletions(-) diff --git a/Common/Tools/Multiplicity/README.md b/Common/Tools/Multiplicity/README.md index 30a2d725fda..c38f81a4f6e 100644 --- a/Common/Tools/Multiplicity/README.md +++ b/Common/Tools/Multiplicity/README.md @@ -1,17 +1,17 @@ -# Multiplicity/centrality tools in O2/O2Physics +# Multiplicity/centrality tools in O2/O2Physics -This directory serves to aggregate all files necessary for multiplicity -and centrality calibration going from pp to Pb-Pb. It offers functionality -to perform simple slicing in percentiles, such as what is done in -proton-proton collisions, as well as the Glauber + analytical NBD fitter -used to perform anchoring and hadronic event distribution estimates in -nucleus-nucleus collisions. +This directory serves to aggregate all files necessary for multiplicity +and centrality calibration going from pp to Pb-Pb. It offers functionality +to perform simple slicing in percentiles, such as what is done in +proton-proton collisions, as well as the Glauber + analytical NBD fitter +used to perform anchoring and hadronic event distribution estimates in +nucleus-nucleus collisions. ## Files * `README.md` this readme * `CMakeLists.txt` definition of source files that need compiling -* `multCalibrator.cxx/h` a class to do percentile slicing of a given histogram. Used for all systems. -* `multMCCalibrator.cxx/h` a class to perform data-to-mc matching of average Nch. -* `multGlauberNBDFitter.cxx/h` a class to do glauber fits. -* `multModule.h` a class to perform calculations of multiplicity and centrality tables for analysis. Meant to be used inside the main core service wagon 'multcenttable'. +* `multCalibrator.cxx/h` a class to do percentile slicing of a given histogram. Used for all systems. +* `multMCCalibrator.cxx/h` a class to perform data-to-mc matching of average Nch. +* `multGlauberNBDFitter.cxx/h` a class to do glauber fits. +* `multModule.h` a class to perform calculations of multiplicity and centrality tables for analysis. Meant to be used inside the main core service wagon 'multcenttable'. diff --git a/Common/Tools/Multiplicity/macros/README.md b/Common/Tools/Multiplicity/macros/README.md index 0054d665d71..1044342cd8b 100644 --- a/Common/Tools/Multiplicity/macros/README.md +++ b/Common/Tools/Multiplicity/macros/README.md @@ -1,69 +1,69 @@ # Example multiplicity calibration macros -You will find some example macros in this directory that will allow for the calculation of a glauber fit and the corresponding calibration histograms. +You will find some example macros in this directory that will allow for the calculation of a glauber fit and the corresponding calibration histograms. -A simplified description of the procedure necessary to generate a -Glauber + NBD fit to a certain histogram is described below. +A simplified description of the procedure necessary to generate a +Glauber + NBD fit to a certain histogram is described below. ## Performing a Glauber + NBD fit -### First step: calculation of Glauber MC sample +### First step: calculation of Glauber MC sample The Glauber + NBD model assumes that the multiplicity / amplitude -distribution that one intends to fit can be described as coming -from a certain number N_{ancestor} of particle-emitting sources called 'ancestors'. Each ancestor is assumed to produce particles -according to a negative binominal distribution. Further, -the number of ancestors is assumed to be related to the basic -glauber quantities N_{part} and N_{coll} via the formula: +distribution that one intends to fit can be described as coming +from a certain number N_{ancestor} of particle-emitting sources called 'ancestors'. Each ancestor is assumed to produce particles +according to a negative binominal distribution. Further, +the number of ancestors is assumed to be related to the basic +glauber quantities N_{part} and N_{coll} via the formula: -N_{ancestor} = f * N_{part} + (1-f) N_{coll} +N_{ancestor} = f * N_{part} + (1-f) N_{coll} Usually, the value f is fixed at 0.8, and thus the range of -N_{ancestors} is typically 0-700 in Pb-Pb collisions. +N_{ancestors} is typically 0-700 in Pb-Pb collisions. -In order to allow for Glauber + NBD fitting, the correlation of -(N_{part}, N_{coll}) needs to be known. For that purpose, -a tree of Glauber MC needs to be generated using TGlauberMC using the relevant nuclei, which also involves choosing an appropriate nuclear profile. +In order to allow for Glauber + NBD fitting, the correlation of +(N_{part}, N_{coll}) needs to be known. For that purpose, +a tree of Glauber MC needs to be generated using TGlauberMC using the relevant nuclei, which also involves choosing an appropriate nuclear profile. -Once TGlauberMC has been used to produce a tree with N_{part} and -N_{coll} values, their correlation needs to be saved to a 2D histogram that serves as input to the Glauber + NBD fitter used in -O2/O2Physics. This is done using the macro called `saveCorrelation.C` contained in this directory, which produces a file named `baseHistos.root`. The file `saveCorrelation.C` serves as +Once TGlauberMC has been used to produce a tree with N_{part} and +N_{coll} values, their correlation needs to be saved to a 2D histogram that serves as input to the Glauber + NBD fitter used in +O2/O2Physics. This is done using the macro called `saveCorrelation.C` contained in this directory, which produces a file named `baseHistos.root`. The file `saveCorrelation.C` serves as an example for the Pb-Pb case and minor adaptation may be needed -in case other nuclei are to be used. +in case other nuclei are to be used. -### Second step: execution of Glauber + NBD fitter +### Second step: execution of Glauber + NBD fitter -The fitting procedure is done via the macro ``. Because -the numerical fitting utilizes a convolution of a N_{ancestor} -distribution and NBD distributions, it will not be as fast +The fitting procedure is done via the macro ``. Because +the numerical fitting utilizes a convolution of a N_{ancestor} +distribution and NBD distributions, it will not be as fast as typical one-dimensional fits: typical times of 10-100 seconds -are not unusual. The macro will produce an output file -that contains: +are not unusual. The macro will produce an output file +that contains: -* the original fitted histogram -* the actual glauber fit distribution, plotted all the way -to zero multiplicity +* the original fitted histogram +* the actual glauber fit distribution, plotted all the way +to zero multiplicity -This output can then be used to extract percentiles. +This output can then be used to extract percentiles. -### Third step: calculation of percentile boundaries +### Third step: calculation of percentile boundaries -Once both the data and glauber fit distributions are known, +Once both the data and glauber fit distributions are known, the next step is to create a 'stitched' data/glauber distribution in -which the distribution at low multiplicities follows -the glauber fit and the distribution at higher multiplicities -follows the data. The multiplicity value in which the switch from -glauber to data takes place is called 'anchor point'. +which the distribution at low multiplicities follows +the glauber fit and the distribution at higher multiplicities +follows the data. The multiplicity value in which the switch from +glauber to data takes place is called 'anchor point'. -Because of the fact that this 'stitching' procedure may need to be tuned and the actual glauber fit is slow, the stitching is done -in a separate macro called `runCalibration.C`. It is provided +Because of the fact that this 'stitching' procedure may need to be tuned and the actual glauber fit is slow, the stitching is done +in a separate macro called `runCalibration.C`. It is provided in this directory as well and it is the third and last step in calculating percentile boundaries. The macro will printout some -key boundaries as well as save an output file with the calibration. +key boundaries as well as save an output file with the calibration. -*Bonus*: at the stage in which the glauber fit has been done +*Bonus*: at the stage in which the glauber fit has been done and all information is available, a de-convolution process -can be used to calculate the average N_{part} and N_{coll} -in centrality slices. This functionality is also provided +can be used to calculate the average N_{part} and N_{coll} +in centrality slices. This functionality is also provided in O2Physics as part of the `multGlauberNBDFitter` and the -`runCalibration.C` macro can optionally also perform that -deconvolution. *Warning*: this procedure might take a mooment. \ No newline at end of file +`runCalibration.C` macro can optionally also perform that +deconvolution. *Warning*: this procedure might take a mooment. \ No newline at end of file From d7b5bf416bacc9f6ecb235a667bc11c8e6d97fe3 Mon Sep 17 00:00:00 2001 From: ddobrigk Date: Sat, 8 Nov 2025 18:09:47 +0100 Subject: [PATCH 6/9] Megalinter issue fix --- Common/Tools/Multiplicity/macros/runCalibration.C | 2 ++ Common/Tools/Multiplicity/macros/runGlauberFit.C | 2 ++ Common/Tools/Multiplicity/macros/saveCorrelation.C | 2 ++ 3 files changed, 6 insertions(+) diff --git a/Common/Tools/Multiplicity/macros/runCalibration.C b/Common/Tools/Multiplicity/macros/runCalibration.C index 0e3fbfdfec8..f15a5e2199f 100644 --- a/Common/Tools/Multiplicity/macros/runCalibration.C +++ b/Common/Tools/Multiplicity/macros/runCalibration.C @@ -10,6 +10,8 @@ // or submit itself to any jurisdiction. // +#include + /// @brief function to calibrate centrality /// @param lInputFileName name of input file. /// @param anchorPointPercentage anchor point percentage to use diff --git a/Common/Tools/Multiplicity/macros/runGlauberFit.C b/Common/Tools/Multiplicity/macros/runGlauberFit.C index e89b77f2400..817791fb3eb 100755 --- a/Common/Tools/Multiplicity/macros/runGlauberFit.C +++ b/Common/Tools/Multiplicity/macros/runGlauberFit.C @@ -28,6 +28,8 @@ #include "TSystem.h" #include "TTree.h" +#include + //________________________________________________________________ Double_t FastIntegrate(TF1* f1, Double_t a, Double_t b, Int_t n = 5) { diff --git a/Common/Tools/Multiplicity/macros/saveCorrelation.C b/Common/Tools/Multiplicity/macros/saveCorrelation.C index 160d0177c69..69f3e481c41 100644 --- a/Common/Tools/Multiplicity/macros/saveCorrelation.C +++ b/Common/Tools/Multiplicity/macros/saveCorrelation.C @@ -13,6 +13,8 @@ // into the 2D correlation histogram necessary for the ALICE machinery // that performs Glauber + NBD fits. +#include + /// @brief function to save Npart x Ncoll correlation to file for glauber fits /// @param filename input TGlauberMC ntuple file /// @param outputFile output file for Npart x Ncoll correlation TH2D From 2cb2ff932829227b6fb720692b8e1233f4b197cf Mon Sep 17 00:00:00 2001 From: ddobrigk Date: Sat, 8 Nov 2025 18:20:05 +0100 Subject: [PATCH 7/9] Fixes from O2 linter --- .../Tools/Multiplicity/macros/runCalibration.C | 3 +++ .../Tools/Multiplicity/macros/runGlauberFit.C | 17 +++++++++++++---- .../Tools/Multiplicity/macros/saveCorrelation.C | 10 ++++++---- 3 files changed, 22 insertions(+), 8 deletions(-) diff --git a/Common/Tools/Multiplicity/macros/runCalibration.C b/Common/Tools/Multiplicity/macros/runCalibration.C index f15a5e2199f..5ce4b52ecd4 100644 --- a/Common/Tools/Multiplicity/macros/runCalibration.C +++ b/Common/Tools/Multiplicity/macros/runCalibration.C @@ -9,6 +9,9 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. // +/// \file saveCorrelation.C +/// \brief +/// \author ALICE #include diff --git a/Common/Tools/Multiplicity/macros/runGlauberFit.C b/Common/Tools/Multiplicity/macros/runGlauberFit.C index 817791fb3eb..18c116eb086 100755 --- a/Common/Tools/Multiplicity/macros/runGlauberFit.C +++ b/Common/Tools/Multiplicity/macros/runGlauberFit.C @@ -9,6 +9,9 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. // +/// \file saveCorrelation.C +/// \brief +/// \author ALICE #include "multCalibrator.h" #include "multGlauberNBDFitter.h" @@ -134,16 +137,22 @@ int runGlauberFit(TString lInputFileName = "AnalysisResultsLHC24ar.root", TStrin // minimum fit range estimate (guess region that may be unfittable) // may need adjusting - Double_t lFitRange = 0.012 * GetBoundaryForPercentile(hV0Mfine, 0.01); - if (lFitRangeMax < 10000) - lFitRange = 0.02 * GetBoundaryForPercentile(hV0Mfine, 0.01); + Double_t maxPercent = 0.01; + Double_t fractionOfMax = 0.012; + Double_t lFitRange = fractionOfMax * GetBoundaryForPercentile(hV0Mfine, maxPercent); + + // adjust if low mult (Ntracks, typically) + Double_t maxRangeForTracks= 10000; + Double_t fractionOfMaxBroader = 0.02; + if (lFitRangeMax < maxRangeForTracks) + lFitRange = fractionOfMaxBroader * GetBoundaryForPercentile(hV0Mfine, maxPercent); cout << "Fit range min estimated from histogram: " << lFitRange << endl; //____________________________________________ // rebinning matters int rebinFactor = 20; - if (lFitRangeMax < 10000) + if (lFitRangeMax < maxRangeForTracks) rebinFactor = 1; cout << "Creating rebinned histogram with rebin factor: " << rebinFactor << endl; diff --git a/Common/Tools/Multiplicity/macros/saveCorrelation.C b/Common/Tools/Multiplicity/macros/saveCorrelation.C index 69f3e481c41..ca0f3540342 100644 --- a/Common/Tools/Multiplicity/macros/saveCorrelation.C +++ b/Common/Tools/Multiplicity/macros/saveCorrelation.C @@ -8,10 +8,12 @@ // In applying this license CERN does not waive the privileges and immunities // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -// -// This file provides a simple macro to convert TGlauberMC output tuples -// into the 2D correlation histogram necessary for the ALICE machinery -// that performs Glauber + NBD fits. +// +/// \file saveCorrelation.C +/// \brief This file provides a simple macro to convert TGlauberMC output tuples +/// into the 2D correlation histogram necessary for the ALICE machinery +/// that performs Glauber + NBD fits. +/// \author ALICE #include From dc2dd3d9f5d43a74a432633725beb1d35a63c197 Mon Sep 17 00:00:00 2001 From: ALICE Builder Date: Sat, 8 Nov 2025 18:31:35 +0100 Subject: [PATCH 8/9] Please consider the following formatting changes (#497) --- Common/Tools/Multiplicity/macros/runCalibration.C | 2 +- Common/Tools/Multiplicity/macros/runGlauberFit.C | 6 +++--- Common/Tools/Multiplicity/macros/saveCorrelation.C | 2 +- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/Common/Tools/Multiplicity/macros/runCalibration.C b/Common/Tools/Multiplicity/macros/runCalibration.C index 5ce4b52ecd4..7a5e76c12d0 100644 --- a/Common/Tools/Multiplicity/macros/runCalibration.C +++ b/Common/Tools/Multiplicity/macros/runCalibration.C @@ -10,7 +10,7 @@ // or submit itself to any jurisdiction. // /// \file saveCorrelation.C -/// \brief +/// \brief /// \author ALICE #include diff --git a/Common/Tools/Multiplicity/macros/runGlauberFit.C b/Common/Tools/Multiplicity/macros/runGlauberFit.C index 18c116eb086..f92d4169bc8 100755 --- a/Common/Tools/Multiplicity/macros/runGlauberFit.C +++ b/Common/Tools/Multiplicity/macros/runGlauberFit.C @@ -10,7 +10,7 @@ // or submit itself to any jurisdiction. // /// \file saveCorrelation.C -/// \brief +/// \brief /// \author ALICE #include "multCalibrator.h" @@ -141,8 +141,8 @@ int runGlauberFit(TString lInputFileName = "AnalysisResultsLHC24ar.root", TStrin Double_t fractionOfMax = 0.012; Double_t lFitRange = fractionOfMax * GetBoundaryForPercentile(hV0Mfine, maxPercent); - // adjust if low mult (Ntracks, typically) - Double_t maxRangeForTracks= 10000; + // adjust if low mult (Ntracks, typically) + Double_t maxRangeForTracks = 10000; Double_t fractionOfMaxBroader = 0.02; if (lFitRangeMax < maxRangeForTracks) lFitRange = fractionOfMaxBroader * GetBoundaryForPercentile(hV0Mfine, maxPercent); diff --git a/Common/Tools/Multiplicity/macros/saveCorrelation.C b/Common/Tools/Multiplicity/macros/saveCorrelation.C index ca0f3540342..b23415639ba 100644 --- a/Common/Tools/Multiplicity/macros/saveCorrelation.C +++ b/Common/Tools/Multiplicity/macros/saveCorrelation.C @@ -8,7 +8,7 @@ // In applying this license CERN does not waive the privileges and immunities // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -// +// /// \file saveCorrelation.C /// \brief This file provides a simple macro to convert TGlauberMC output tuples /// into the 2D correlation histogram necessary for the ALICE machinery From 02c70fc20e5b9951513a393f28e7fd56f78d9c8c Mon Sep 17 00:00:00 2001 From: ddobrigk Date: Sat, 8 Nov 2025 18:33:27 +0100 Subject: [PATCH 9/9] O2 linter fixes --- Common/Tools/Multiplicity/macros/runGlauberFit.C | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/Common/Tools/Multiplicity/macros/runGlauberFit.C b/Common/Tools/Multiplicity/macros/runGlauberFit.C index f92d4169bc8..fe33c209308 100755 --- a/Common/Tools/Multiplicity/macros/runGlauberFit.C +++ b/Common/Tools/Multiplicity/macros/runGlauberFit.C @@ -281,7 +281,8 @@ int runGlauberFit(TString lInputFileName = "AnalysisResultsLHC24ar.root", TStrin Int_t lFitStatus = 0; lFitStatus = g->DoFit(); Int_t lAttempts = 1; - while (lAttempts < 10 && lFitStatus == 0) { + Int_t lMaxAttempts = 10; + while (lAttempts < lMaxAttempts && lFitStatus == 0) { // insist on fitting until it works cout << "Attempting fit again (" << lAttempts << " attempt)..." << endl; lFitStatus = g->DoFit(); @@ -303,7 +304,8 @@ int runGlauberFit(TString lInputFileName = "AnalysisResultsLHC24ar.root", TStrin for (Int_t ii = 1; ii < hGlauber->GetNbinsX() + 1; ii++) { Double_t lFuncVal = FastIntegrate(fitfunc, hGlauber->GetBinLowEdge(ii), hGlauber->GetBinLowEdge(ii + 1), 4); hGlauber->SetBinContent(ii, lFuncVal); - if (ii % 500 == 0) { + Int_t printEveryThisManyBins = 500; + if (ii % printEveryThisManyBins == 0) { cout << "At integration #" << ii << "/" << hGlauber->GetNbinsX() + 1 << "..." << endl; } }