From 3bb40a3df4af5302f131c0ac1acf9eadcabdf0f3 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Sat, 8 Nov 2025 16:49:54 +0000 Subject: [PATCH] Please consider the following formatting changes --- .../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();