Skip to content

Commit

Permalink
Add xX0 scan macro
Browse files Browse the repository at this point in the history
  • Loading branch information
@mconcas
mconcas committed Jan 21, 2024
1 parent e0e1365 commit d272137
Show file tree
Hide file tree
Showing 3 changed files with 315 additions and 1 deletion.
3 changes: 2 additions & 1 deletion Detectors/Upgrades/ALICE3/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,4 +15,5 @@ add_subdirectory(FT3)
add_subdirectory(FCT)
add_subdirectory(AOD)
add_subdirectory(IOTOF)
add_subdirectory(RICH)
add_subdirectory(RICH)
add_subdirectory(macros)
13 changes: 13 additions & 0 deletions Detectors/Upgrades/ALICE3/macros/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,13 @@
# 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.

o2_add_test_root_macro(scanXX0.C
LABELS alice3)
300 changes: 300 additions & 0 deletions Detectors/Upgrades/ALICE3/macros/scanXX0.C
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,300 @@
// 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.
//
// Original authors: M. Sitta, F. Grosa
// Author: M. Concas

#if !defined(__CLING__) || defined(__ROOTCLING__)

#include "DetectorsBase/GeometryManager.h"
#include "ITSBase/GeometryTGeo.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TGeoManager.h"
#include "TH2F.h"
#include "TMath.h"
#include "TRandom.h"
#include "TSystem.h"
#include <TLegend.h>
#include <TStyle.h>

#include "TGeoMaterial.h"
#include "TGeoMedium.h"
#include <TGeoTube.h>
#include <TGeoVolume.h>
#include <TPaveText.h>
#include <TROOT.h>
#include <TText.h>

// #include "FairGeoParSet.h" // for FairGeoParSet
#include <fairlogger/Logger.h> // for LOG, LOG_IF

#endif

//______________________________________________________________________
void ComputeMaterialBudget(double rmin, double rmax, double etaPos,
double phiMin, double phiMax, TH1F* xOverX0VsPhi,
TH1F* xOverX0VsEta, TH1F* xOverX0VsZvtx, TH2F* xOverX02d = nullptr)
{
// Ancillary function to compute material budget between rmin and rmax

int n = 1e6; // testtracks

float len = 900; // cm
TH1F* nParticlesVsPhi = new TH1F("", "", 90, phiMin, phiMax);
TH1F* nParticlesVsEta = new TH1F("", "", 200, -2, 2);
TH1F* nParticlesVsZvtx = new TH1F("", "", 300, -len / 2, len / 2);
TH2F* nParticlesVsPhiEta = new TH2F("", "", 200, -2, 2, 45, phiMax, phiMax);

double x1, y1, z1, x2, y2, z2;

for (int it = 0; it < n; it++) { // we simulate flat in phi and eta, from Zvtx=0

// PHI VS ETA
double phi = gRandom->Uniform(phiMin, phiMax);
double eta = gRandom->Uniform(-4, 4);
double theta = TMath::ATan(TMath::Exp(-eta)) * 2;
x1 = rmin * TMath::Cos(phi);
y1 = rmin * TMath::Sin(phi);
z1 = rmin / TMath::Tan(theta);
x2 = rmax * TMath::Cos(phi);
y2 = rmax * TMath::Sin(phi);
z2 = rmax / TMath::Tan(theta);

auto mparam = o2::base::GeometryManager::meanMaterialBudget(x1, y1, z1, x2, y2, z2);

if (std::abs(eta) < etaPos) {
xOverX0VsPhi->Fill(phi, mparam.meanX2X0 * 100);
nParticlesVsPhi->Fill(phi);
}
xOverX0VsEta->Fill(eta, mparam.meanX2X0 * 100);
nParticlesVsEta->Fill(eta);

if (xOverX02d != 0) {
xOverX02d->Fill(eta, phi, mparam.meanX2X0 * 100);
nParticlesVsPhiEta->Fill(eta, phi);
}
}

for (int it = 0; it < n; it++) { // then we simulate flat in Zvtx with eta = 0 and flat in phi

double phi = TMath::Pi() / 2;
double eta = 0.;
double theta = TMath::ATan(TMath::Exp(-eta)) * 2;
double z = gRandom->Uniform(-len / 2, len / 2);
x1 = rmin * TMath::Cos(phi);
y1 = rmin * TMath::Sin(phi);
x2 = rmax * TMath::Cos(phi);
y2 = rmax * TMath::Sin(phi);

auto mparam = o2::base::GeometryManager::meanMaterialBudget(x1, y1, z, x2, y2, z);

xOverX0VsZvtx->Fill(z, mparam.meanX2X0 * 100);
nParticlesVsZvtx->Fill(z);
}

// normalization to number of particles in case of phi vs eta
double theta = TMath::ATan(TMath::Exp(-etaPos / 2)) * 2;
LOGP(info, "<η>={} -> Sin(θ) {}", etaPos / 2, TMath::Sin(theta));

for (int ix = 1; ix <= nParticlesVsPhi->GetNbinsX(); ix++) {
if (nParticlesVsPhi->GetBinContent(ix) > 0)
xOverX0VsPhi->SetBinContent(ix, xOverX0VsPhi->GetBinContent(ix) / nParticlesVsPhi->GetBinContent(ix) * TMath::Sin(theta));
}

for (int ix = 1; ix <= nParticlesVsEta->GetNbinsX(); ix++) {
if (nParticlesVsEta->GetBinContent(ix) > 0)
xOverX0VsEta->SetBinContent(ix, xOverX0VsEta->GetBinContent(ix) / nParticlesVsEta->GetBinContent(ix));
}

for (int ix = 1; ix <= nParticlesVsZvtx->GetNbinsX(); ix++) {
if (nParticlesVsZvtx->GetBinContent(ix) > 0)
xOverX0VsZvtx->SetBinContent(ix, xOverX0VsZvtx->GetBinContent(ix) / nParticlesVsZvtx->GetBinContent(ix));
}

if (xOverX02d) {
for (int ix = 1; ix <= nParticlesVsPhiEta->GetNbinsX(); ix++) {
for (int iy = 1; iy <= nParticlesVsPhiEta->GetNbinsY(); iy++) {
if (nParticlesVsPhiEta->GetBinContent(ix, iy) > 0)
xOverX02d->SetBinContent(ix, iy, xOverX02d->GetBinContent(ix, iy) / nParticlesVsPhiEta->GetBinContent(ix, iy));
}
}
}
}

void printMaterialDefinitions(TGeoManager* gman)
{
TGeoMedium* med;
TGeoMaterial* mat;
char mediaName[50], matName[50], shortName[50];

Int_t nMedia = gman->GetListOfMedia()->GetEntries();

LOGP(info, " =================== ALICE 3 Material Properties ================= ");
LOGP(info, " A Z d (g/cm3) RadLen (cm) IntLen (cm)\t Name\n");

// Loop on media, select ITS3 materials, print their characteristics
for (Int_t i = 0; i < nMedia; i++) {
med = (TGeoMedium*)(gman->GetListOfMedia()->At(i));
mat = med->GetMaterial();
LOGP(info, "{:5.1f} {:6.1f} {:8.3f} {:13.3f} {:11.3f}\t {}", mat->GetA(), mat->GetZ(), mat->GetDensity(), mat->GetRadLen(), mat->GetIntLen(), mat->GetName());
}
}

void scanXX0(const string fileName = "o2sim_geometry.root", const string path = "./")
{
gStyle->SetPadTopMargin(0.035);
gStyle->SetPadRightMargin(0.035);
gStyle->SetPadBottomMargin(0.14);
gStyle->SetPadLeftMargin(0.14);
gStyle->SetTitleOffset(1.4, "y");
gStyle->SetPadTickX(1);
gStyle->SetPadTickY(1);
gStyle->SetOptStat(0);
gStyle->SetOptTitle(0);

double etaPos = 1.;

TCanvas* canv = new TCanvas("canv", "canv", 2400, 800);
canv->Divide(3, 1);

TLegend* legVsPhi = new TLegend(0.25, 0.6, 0.85, 0.9);
legVsPhi->SetFillColor(kWhite);
legVsPhi->SetTextSize(0.045);
legVsPhi->SetHeader(Form("ALICE 3, |#it{#eta}| < %0.f, #it{Z}_{vtx} = 0", etaPos));

TLegend* legVsEta = new TLegend(0.25, 0.6, 0.85, 0.9);
legVsEta->SetFillColor(kWhite);
legVsEta->SetTextSize(0.045);
legVsEta->SetHeader("ALICE 3, 0 < #it{#varphi} < #pi, #it{Z}_{vtx} = 0");

TLegend* legVsZvtx = new TLegend(0.25, 0.6, 0.85, 0.9);
legVsZvtx->SetFillColor(kWhite);
legVsZvtx->SetTextSize(0.045);
legVsZvtx->SetHeader("ALICE 3, #it{#varphi} = #pi/2, #it{#eta} = 0");

std::vector<TH1F*> xOverX0VsPhi;
std::vector<TH1F*> xOverX0VsEta;
std::vector<TH1F*> xOverX0VsZvtx;

TGeoManager::Import((path + fileName).c_str());
printMaterialDefinitions(gGeoManager);

const double rmin = 0.2;
const double rmax = 100.;
const double phiMin = 0;
const double phiMax = 2 * TMath::Pi();
const double len = 900.;
std::array<int, 3> colors = {kGray + 2, kAzure + 4, kRed + 1};

for (size_t iMedium{0}; gGeoManager->GetListOfMedia()->GetEntries(); ++iMedium) {
xOverX0VsPhi.emplace_back(new TH1F(Form("xOverX0VsPhi_step%zu", iMedium), "", 90, phiMin, phiMax));
xOverX0VsEta.emplace_back(new TH1F(Form("xOverX0VsEta_step%zu", iMedium), "", 200, -2, 2));
xOverX0VsZvtx.emplace_back(new TH1F(Form("xOverX0VsZvtx_step%zu", iMedium), "", 300, -len / 2, len / 2));

ComputeMaterialBudget(rmin, rmax, etaPos, phiMin, phiMax, xOverX0VsPhi.back(), xOverX0VsEta.back(), xOverX0VsZvtx.back());

double meanX0vsPhi = 0, meanX0vsEta = 0, meanX0vsZvtx = 0;
for (int ix = 1; ix <= xOverX0VsPhi.back()->GetNbinsX(); ix++) {
meanX0vsPhi += xOverX0VsPhi.back()->GetBinContent(ix);
}
meanX0vsPhi /= xOverX0VsPhi.back()->GetNbinsX();

for (int ix = 1; ix <= xOverX0VsEta.back()->GetNbinsX(); ix++) {
meanX0vsEta += xOverX0VsEta.back()->GetBinContent(ix);
}
meanX0vsEta /= xOverX0VsEta.back()->GetNbinsX();

for (int ix = 1; ix <= xOverX0VsZvtx.back()->GetNbinsX(); ix++) {
meanX0vsZvtx += xOverX0VsZvtx.back()->GetBinContent(ix);
}
meanX0vsZvtx /= xOverX0VsZvtx.back()->GetNbinsX();

LOGP(info, "Mean X/X0 vs. phi: {}", meanX0vsPhi);
LOGP(info, "Mean X/X0 vs. eta: {}", meanX0vsEta);
LOGP(info, "Mean X/X0 vs. Zvtx: {}", meanX0vsZvtx);

xOverX0VsPhi.back()->GetXaxis()->SetTitle("#it{#varphi} (rad)");
xOverX0VsPhi.back()->GetXaxis()->SetTitleSize(0.05);
xOverX0VsPhi.back()->GetXaxis()->SetLabelSize(0.045);
xOverX0VsPhi.back()->GetYaxis()->SetTitle("#it{X}/#it{X}_{0} (%)");
xOverX0VsPhi.back()->GetYaxis()->SetTitleSize(0.05);
xOverX0VsPhi.back()->GetYaxis()->SetLabelSize(0.045);
xOverX0VsPhi.back()->GetYaxis()->SetDecimals();
xOverX0VsPhi.back()->SetFillColorAlpha(colors[0], 0.5);
xOverX0VsPhi.back()->SetLineColor(colors[0]);
xOverX0VsPhi.back()->SetLineWidth(2);

xOverX0VsEta.back()->GetXaxis()->SetTitle("#it{#eta}");
xOverX0VsEta.back()->GetXaxis()->SetTitleSize(0.05);
xOverX0VsEta.back()->GetXaxis()->SetLabelSize(0.045);
xOverX0VsEta.back()->GetYaxis()->SetTitle("#it{X}/#it{X}_{0} (%)");
xOverX0VsEta.back()->GetYaxis()->SetTitleSize(0.05);
xOverX0VsEta.back()->GetYaxis()->SetLabelSize(0.045);
xOverX0VsEta.back()->GetYaxis()->SetDecimals();
xOverX0VsEta.back()->SetFillColorAlpha(colors[0], 0.5);
xOverX0VsEta.back()->SetLineColor(colors[0]);
xOverX0VsEta.back()->SetLineWidth(2);

xOverX0VsZvtx.back()->GetXaxis()->SetTitle("#it{Z}_{vtx} (cm)");
xOverX0VsZvtx.back()->GetXaxis()->SetTitleSize(0.05);
xOverX0VsZvtx.back()->GetXaxis()->SetLabelSize(0.045);
xOverX0VsZvtx.back()->GetYaxis()->SetTitle("#it{X}/#it{X}_{0} (%)");
xOverX0VsZvtx.back()->GetYaxis()->SetTitleSize(0.05);
xOverX0VsZvtx.back()->GetYaxis()->SetLabelSize(0.045);
xOverX0VsZvtx.back()->GetYaxis()->SetDecimals();
xOverX0VsZvtx.back()->SetFillColorAlpha(colors[0], 0.5);
xOverX0VsZvtx.back()->SetLineColor(colors[0]);
xOverX0VsZvtx.back()->SetLineWidth(2);

if (xOverX0VsPhi.size() == 1) {
legVsPhi->AddEntry("", Form("#LT #it{X}/#it{X}_{0} #GT = %0.3f %%", meanX0vsPhi), "");
legVsEta->AddEntry("", Form("#LT #it{X}/#it{X}_{0} #GT = %0.3f %%", meanX0vsEta), "");
legVsZvtx->AddEntry("", Form("#LT #it{X}/#it{X}_{0} #GT = %0.3f %%", meanX0vsZvtx), "");
}
legVsPhi->AddEntry(xOverX0VsPhi.back(), "Total", "f");
legVsEta->AddEntry(xOverX0VsPhi.back(), "Total", "f");
legVsZvtx->AddEntry(xOverX0VsZvtx.back(), "Total", "f");

canv->cd(1)->SetGrid();
if (xOverX0VsPhi.size() == 1) {
xOverX0VsPhi.back()->SetMinimum(1.e-4);
xOverX0VsPhi.back()->SetMaximum(20.f);
xOverX0VsPhi.back()->DrawCopy("HISTO");
legVsPhi->Draw();
} else {
xOverX0VsPhi.back()->DrawCopy("HISTO SAME");
}

canv->cd(2)->SetGrid();
if (xOverX0VsEta.size() == 1) {
xOverX0VsEta.back()->SetMinimum(1.e-4);
xOverX0VsEta.back()->SetMaximum(60.f);
xOverX0VsEta.back()->DrawCopy("HISTO");
legVsEta->Draw();
} else {
xOverX0VsEta.back()->DrawCopy("HISTO SAME");
}

canv->cd(3)->SetGrid();
if (xOverX0VsZvtx.size() == 1) {
xOverX0VsZvtx.back()->SetMinimum(1.e-4);
xOverX0VsZvtx.back()->SetMaximum(100.f);
xOverX0VsZvtx.back()->DrawCopy("HISTO");
// firstPlot = 0;
legVsZvtx->Draw();
} else {
xOverX0VsZvtx.back()->DrawCopy("HISTO SAME");
}
break;
}
canv->SaveAs("alice3_material_vsphietaz.pdf");
}

0 comments on commit d272137

Please sign in to comment.