Skip to content

Commit 43a0f6e

Browse files
urasantonioAntonio Urasnjacazio
authored
Introducing the MID (Muon Identifier) detector (#63)
Introducing PID information based on the MID (Muon Identifier) detector * Adding MID (Muon Identifier) information in the O2 tables * MID information is written in the AOD tables only if the setup of MIDdetector is a valid one * Creating O2 geometry files automatically * Fix in the MID information for the O2 tables * Fix in MID information for the O2 tables * Further fix to MID information in O2 tables * Further fix for MID information in O2 tables * Clang fix * Applying Nicolo's comments Co-authored-by: Antonio Uras <auras@alicecerno2.cern.ch> Co-authored-by: Nicolò Jacazio <nicolo.jacazio@cern.ch> Co-authored-by: Nicolò Jacazio <njacazio@users.noreply.github.com>
1 parent dc42182 commit 43a0f6e

8 files changed

Lines changed: 337 additions & 7 deletions

File tree

‎examples/aod/createO2tables.C‎

Lines changed: 23 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -2,8 +2,8 @@ R__LOAD_LIBRARY(libDelphes)
22
R__LOAD_LIBRARY(libDelphesO2)
33

44
#include <algorithm> // std::shuffle
5-
#include <random> // std::default_random_engine
6-
#include <chrono> // std::chrono::system_clock
5+
#include <random> // std::default_random_engine
6+
#include <chrono> // std::chrono::system_clock
77

88
// ROOT includes
99
#include "TMath.h"
@@ -26,6 +26,7 @@ R__LOAD_LIBRARY(libDelphesO2)
2626
#include "TrackSmearer.hh"
2727
#include "TOFLayer.hh"
2828
#include "RICHdetector.hh"
29+
#include "MIDdetector.hh"
2930
#include "TrackUtils.hh"
3031

3132
#include "createO2tables.h"
@@ -44,6 +45,8 @@ const double rich_index = 1.03;
4445
const double rich_radiator_length = 2.;
4546
const double rich_efficiency = 0.4;
4647
const double rich_sigma = 7.e-3;
48+
// MID
49+
const char* inputFileAccMuonPID = "muonAccEffPID.root";
4750

4851
// Simulation parameters
4952
const bool do_vertexing = true;
@@ -154,6 +157,11 @@ int createO2tables(const char* inputFile = "delphes.root",
154157
richdetector.setRadiatorLength(rich_radiator_length);
155158
richdetector.setEfficiency(rich_efficiency);
156159
richdetector.setSigma(rich_sigma);
160+
// MID detector
161+
o2::delphes::MIDdetector midDetector;
162+
printf("creating MID detector...\n");
163+
bool isMID = midDetector.setup(inputFileAccMuonPID);
164+
printf("isMID = %d\n", isMID);
157165

158166
// create output
159167
auto fout = TFile::Open(outputFile, "RECREATE");
@@ -163,6 +171,7 @@ int createO2tables(const char* inputFile = "delphes.root",
163171
MakeTreeO2trackCov();
164172
MakeTreeO2trackExtra();
165173
MakeTreeO2rich();
174+
MakeTreeO2mid();
166175
MakeTreeO2collision();
167176
MakeTreeO2collisionExtra();
168177
MakeTreeO2mccollision();
@@ -185,8 +194,7 @@ int createO2tables(const char* inputFile = "delphes.root",
185194

186195
// Random generator for reshuffling tracks when reading them
187196
std::default_random_engine e(std::chrono::system_clock::now().time_since_epoch().count()); // time-based seed:
188-
189-
for (Int_t ientry = 0; ientry < numberOfEntries; ++ientry) { // Loop over events
197+
for (Int_t ientry = 0; ientry < numberOfEntries; ++ientry) { // Loop over events
190198
// Adjust start indices for this event in all trees by adding the number of entries of the previous event
191199
for (auto i = 0; i < kTrees; ++i) {
192200
eventextra.fStart[i] += eventextra.fNentries[i];
@@ -240,6 +248,8 @@ int createO2tables(const char* inputFile = "delphes.root",
240248
std::vector<TrackAlice3> tracks_for_vertexing;
241249
std::vector<o2::InteractionRecord> bcData;
242250
std::vector<Track*> tof_tracks;
251+
const int multiplicity = tracks->GetEntries();
252+
243253
// Build index array of tracks to randomize track writing order
244254
std::vector<int> tracks_indices(tracks->GetEntries()); // vector with tracks->GetEntries()
245255
std::iota(std::begin(tracks_indices), std::end(tracks_indices), 0); // Fill with 0, 1, ...
@@ -348,6 +358,15 @@ int createO2tables(const char* inputFile = "delphes.root",
348358
rich.fRICHNsigmaPr = nsigma[4];
349359
FillTree(kRICH);
350360
}
361+
// check if it is within the acceptance of the MID
362+
if (isMID) {
363+
if (midDetector.hasMID(*track)) {
364+
mid.fIndexCollisions = ientry + eventOffset;
365+
mid.fIndexTracks = fTrackCounter; // Index in the Track table
366+
mid.fMIDIsMuon = midDetector.isMuon(*track, multiplicity);
367+
FillTree(kMID);
368+
}
369+
}
351370
if (do_vertexing) {
352371
o2::InteractionRecord ir(ientry + eventOffset, 0);
353372
const float t = (ir.bc2ns() + gRandom->Gaus(0., 100.)) * 1e-3;

‎examples/aod/createO2tables.h‎

Lines changed: 19 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -27,14 +27,15 @@ enum TreeIndex { // Index of the output trees
2727
kRun2BCInfo,
2828
kHMPID,
2929
kRICH,
30+
kMID,
3031
kTrees
3132
};
3233

3334
const int fBasketSizeEvents = 1000000; // Maximum basket size of the trees for events
3435
const int fBasketSizeTracks = 10000000; // Maximum basket size of the trees for tracks
3536

36-
const TString TreeName[kTrees] = {"O2collision", "DbgEventExtra", "O2track", "O2trackcov", "O2trackextra", "O2calo", "O2calotrigger", "O2muon", "O2muoncluster", "O2zdc", "O2fv0a", "O2fv0c", "O2ft0", "O2fdd", "O2v0", "O2cascade", "O2tof", "O2mcparticle", "O2mccollision", "O2mctracklabel", "O2mccalolabel", "O2mccollisionlabel", "O2bc", "O2run2bcinfo", "O2hmpid", "O2rich"};
37-
const TString TreeTitle[kTrees] = {"Collision tree", "Collision extra", "Barrel tracks Parameters", "Barrel tracks Covariance", "Barrel tracks Extra", "Calorimeter cells", "Calorimeter triggers", "MUON tracks", "MUON clusters", "ZDC", "FV0A", "FV0C", "FT0", "FDD", "V0s", "Cascades", "TOF hits", "Kinematics", "MC collisions", "MC track labels", "MC calo labels", "MC collision labels", "BC info", "Run 2 BC Info", "HMPID info", "RICH info"};
37+
const TString TreeName[kTrees] = {"O2collision", "DbgEventExtra", "O2track", "O2trackcov", "O2trackextra", "O2calo", "O2calotrigger", "O2muon", "O2muoncluster", "O2zdc", "O2fv0a", "O2fv0c", "O2ft0", "O2fdd", "O2v0", "O2cascade", "O2tof", "O2mcparticle", "O2mccollision", "O2mctracklabel", "O2mccalolabel", "O2mccollisionlabel", "O2bc", "O2run2bcinfo", "O2hmpid", "O2rich","O2mid"};
38+
const TString TreeTitle[kTrees] = {"Collision tree", "Collision extra", "Barrel tracks Parameters", "Barrel tracks Covariance", "Barrel tracks Extra", "Calorimeter cells", "Calorimeter triggers", "MUON tracks", "MUON clusters", "ZDC", "FV0A", "FV0C", "FT0", "FDD", "V0s", "Cascades", "TOF hits", "Kinematics", "MC collisions", "MC track labels", "MC calo labels", "MC collision labels", "BC info", "Run 2 BC Info", "HMPID info", "RICH info", "MID info"};
3839

3940
TTree* Trees[kTrees] = {nullptr}; // Array of created TTrees
4041
TTree* CreateTree(TreeIndex t)
@@ -329,6 +330,22 @@ void MakeTreeO2rich()
329330
t->SetBasketSize("*", fBasketSizeTracks);
330331
}
331332

333+
struct {
334+
// MID data
335+
Int_t fIndexCollisions = -1; /// Collision ID
336+
Int_t fIndexTracks = -1; /// Track ID
337+
Bool_t fMIDIsMuon = kFALSE; /// MID response for the muon hypothesis
338+
} mid; //! structure to keep MID info
339+
340+
void MakeTreeO2mid()
341+
{
342+
TTree* t = CreateTree(kMID);
343+
t->Branch("fIndexCollisions", &mid.fIndexCollisions, "fIndexCollisions/I");
344+
t->Branch("fIndexTracks", &mid.fIndexTracks, "fIndexTracks/I");
345+
t->Branch("fMIDIsMuon", &mid.fMIDIsMuon, "fMIDIsMuon/b");
346+
t->SetBasketSize("*", fBasketSizeTracks);
347+
}
348+
332349
struct {
333350
// MC particle
334351

‎examples/scripts/muon.sh‎

Lines changed: 89 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,89 @@
1+
#! /usr/bin/env bash
2+
3+
NJOBS=1 # number of max parallel runs
4+
NRUNS=1 # number of runs
5+
NEVENTS=100 # number of events in a run
6+
7+
LUTTAG="werner" # LUT tag name
8+
PY8CFG="pythia8_ccbar" # pythia8 configuration
9+
BFIELD=5. # magnetic field [kG]
10+
SIGMAT=0.020 # time resolution [ns]
11+
SIGMA0=0.200 # vertex time spread [ns]
12+
TAILLX=1.0 # tail on left [q]
13+
TAILRX=1.0 # tail on right [q]
14+
TOFRAD=100. # TOF radius [cm]
15+
TOFLEN=200. # TOF half length [cm]
16+
TOFETA=1.443 # TOF max pseudorapidity
17+
18+
### calculate max eta from geometry
19+
TOFETA=`awk -v a=$TOFRAD -v b=$TOFLEN 'BEGIN {th=atan2(a,b)*0.5; sth=sin(th); cth=cos(th); print -log(sth/cth)}'`
20+
echo "maxEta = $TOFETA"
21+
22+
### copy relevant files in the working directory
23+
cp $DELPHESO2_ROOT/examples/cards/propagate.2kG.tails.tcl propagate.tcl
24+
cp $DELPHESO2_ROOT/examples/pythia8/$PY8CFG.cfg pythia8.cfg
25+
#cp $DELPHESO2_ROOT/examples/smearing/muon.C .
26+
cp ../smearing/muon.C .
27+
28+
### adjust pythia8 configuration
29+
echo "" >> pythia8.cfg
30+
echo "### run time configuration" >> pythia8.cfg
31+
echo "Main:numberOfEvents $NEVENTS" >> pythia8.cfg
32+
echo "Beams:allowVertexSpread on " >> pythia8.cfg
33+
echo "Beams:sigmaTime 60." >> pythia8.cfg
34+
35+
### set magnetic field
36+
sed -i -e "s/set barrel_Bz .*$/set barrel_Bz ${BFIELD}e\-1/" propagate.tcl
37+
### set TOF radius
38+
sed -i -e "s/set barrel_Radius .*$/set barrel_Radius ${TOFRAD}e\-2/" propagate.tcl
39+
sed -i -e "s/double tof_radius = .*$/double tof_radius = ${TOFRAD}\;/" muon.C
40+
### set TOF length
41+
sed -i -e "s/set barrel_HalfLength .*$/set barrel_HalfLength ${TOFLEN}e\-2/" propagate.tcl
42+
sed -i -e "s/double tof_length = .*$/double tof_length = ${TOFLEN}\;/" muon.C
43+
### set TOF acceptance
44+
sed -i -e "s/set barrel_Acceptance .*$/set barrel_Acceptance \{ 0.0 + 1.0 * fabs(eta) < ${TOFETA} \}/" propagate.tcl
45+
### set TOF time resolution and tails
46+
sed -i -e "s/set barrel_TimeResolution .*$/set barrel_TimeResolution ${SIGMAT}e\-9/" propagate.tcl
47+
sed -i -e "s/set barrel_TailRight .*$/set barrel_TailRight ${TAILRX}/" propagate.tcl
48+
sed -i -e "s/set barrel_TailLeft .*$/set barrel_TailLeft ${TAILLX}/" propagate.tcl
49+
sed -i -e "s/double tof_sigmat = .*$/double tof_sigmat = ${SIGMAT}\;/" muon.C
50+
sed -i -e "s/double tof_sigma0 = .*$/double tof_sigma0 = ${SIGMA0}\;/" muon.C
51+
### set TOF mismatch information
52+
sed -i -e "s/double tof_mismatch.*$/double tof_mismatch = 0.01;/" muon.C
53+
#sed -i -e "s/std::string tof_mismatch_fname.*$/std::string tof_mismatch_fname = \"tof_mismatch_template.root\";/" muon.C
54+
55+
### create LUTs
56+
BFIELDT=`awk -v a=$BFIELD 'BEGIN {print a*0.1}'`
57+
$DELPHESO2_ROOT/examples/scripts/create_luts.sh $LUTTAG $BFIELDT $TOFRAD
58+
59+
### loop over runs
60+
rm -f .running.* delphes.*.root
61+
for I in $(seq 1 $NRUNS); do
62+
63+
### wait for a free slot
64+
while [ $(ls .running.* 2> /dev/null | wc -l) -ge $NJOBS ]; do
65+
echo " --- waiting for a free slot"
66+
sleep 1
67+
done
68+
69+
### book the slot
70+
echo " --- starting run $I"
71+
touch .running.$I
72+
73+
### copy pythia8 configuration and set random seed
74+
cp pythia8.cfg pythia8.$I.cfg
75+
echo "Random:setSeed on" >> pythia8.$I.cfg
76+
echo "Random:seed = $I" >> pythia8.$I.cfg
77+
78+
### run Delphes and analysis
79+
DelphesPythia8 propagate.tcl pythia8.$I.cfg delphes.$I.root &> delphes.$I.log &&
80+
root -b -q -l "muon.C(\"delphes.$I.root\",\"muonAccEffPID.root\",\"muon.$I.root\")" &> muon.$I.log &&
81+
rm -rf delphes.$I.root &&
82+
rm -rf .running.$I &&
83+
echo " --- complete run $I" &
84+
85+
done
86+
87+
### merge runs when all done
88+
wait
89+
hadd -f muon.root muon.*.root && rm -rf muon.*.root

‎examples/smearing/muon.C‎

Lines changed: 78 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,78 @@
1+
R__LOAD_LIBRARY(libDelphes)
2+
R__LOAD_LIBRARY(libDelphesO2)
3+
4+
void muon(const char *inputFile = "delphes.root",
5+
const char *inputFileAccMuonPID = "muonAccEffPID.root",
6+
const char *outputFile = "muon.root") {
7+
8+
// Create chain of root trees
9+
TChain chain("Delphes");
10+
chain.Add(inputFile);
11+
12+
// Create object of class ExRootTreeReader
13+
auto treeReader = new ExRootTreeReader(&chain);
14+
auto numberOfEntries = treeReader->GetEntries();
15+
16+
// Get pointers to branches used in this analysis
17+
auto events = treeReader->UseBranch("Event");
18+
auto tracks = treeReader->UseBranch("Track");
19+
auto particles = treeReader->UseBranch("Particle");
20+
21+
// MID detector
22+
o2::delphes::MIDdetector mid;
23+
mid.setup(inputFileAccMuonPID);
24+
25+
// smearer
26+
// o2::delphes::TrackSmearer smearer;
27+
// smearer.useEfficiency(true);
28+
// smearer.loadTable(11, "lutCovm.el.dat");
29+
// smearer.loadTable(13, "lutCovm.mu.dat");
30+
// smearer.loadTable(211, "lutCovm.pi.dat");
31+
// smearer.loadTable(321, "lutCovm.ka.dat");
32+
// smearer.loadTable(2212, "lutCovm.pr.dat");
33+
34+
// logx binning
35+
const Int_t nbins = 80;
36+
double xmin = 1.e-2;
37+
double xmax = 1.e2;
38+
double logxmin = std::log10(xmin);
39+
double logxmax = std::log10(xmax);
40+
double binwidth = (logxmax - logxmin) / nbins;
41+
double xbins[nbins + 1];
42+
for (Int_t i=0; i<=nbins; ++i) xbins[i] = std::pow(10., logxmin + i * binwidth);
43+
44+
auto hPtMuons = new TH1F("hPtMuons", "hPtMuons;#it{p_{T}} (GeV/#it{c});", nbins, xbins);
45+
auto hPtAll = new TH1F("hPtAll", "hPtAll;#it{p_{T}} (GeV/#it{c});", nbins, xbins);
46+
47+
for (Int_t ientry = 0; ientry < numberOfEntries; ++ientry) {
48+
49+
// Load selected branches with data from specified event
50+
treeReader->ReadEntry(ientry);
51+
52+
// loop over tracks
53+
Int_t multiplicity = tracks->GetEntries();
54+
for (Int_t itrack = 0; itrack < tracks->GetEntries(); ++itrack) {
55+
56+
auto track = (Track*) tracks->At(itrack);
57+
58+
// smear track
59+
// if (!smearer.smearTrack(*track)) continue;
60+
61+
// check if has MID
62+
if (mid.hasMID(*track)) {
63+
hPtAll->Fill(track->PT);
64+
if (mid.isMuon(*track, multiplicity)) hPtMuons->Fill(track->PT);
65+
}
66+
67+
}
68+
69+
}
70+
71+
auto fout = TFile::Open(outputFile, "RECREATE");
72+
73+
hPtMuons ->Write();
74+
hPtAll ->Write();
75+
76+
fout->Close();
77+
78+
}

‎src/CMakeLists.txt‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,14 +6,16 @@ set(SOURCES
66
TrackUtils.cc
77
TOFLayer.cc
88
RICHdetector.cc
9+
MIDdetector.cc
910
)
1011

1112
set(HEADERS
1213
VertexFitter.hh
13-
TrackSmearer.hh
14+
TrackSmearer.hh
1415
TrackUtils.hh
1516
TOFLayer.hh
1617
RICHdetector.hh
18+
MIDdetector.hh
1719
)
1820

1921
get_target_property(DELPHES_INCLUDE_DIRECTORIES

‎src/DelphesO2LinkDef.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,7 @@
1212
#pragma link C++ class o2::delphes::TrackUtils+;
1313
#pragma link C++ class o2::delphes::TOFLayer+;
1414
#pragma link C++ class o2::delphes::RICHdetector+;
15+
#pragma link C++ class o2::delphes::MIDdetector+;
1516
#pragma link C++ struct o2::delphes::Vertex+;
1617

1718
#endif

‎src/MIDdetector.cc‎

Lines changed: 78 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,78 @@
1+
/// @author: Roberto Preghenella
2+
/// @email: preghenella@bo.infn.it
3+
4+
/// @author: Antonio Uras
5+
/// @email: antonio.uras@cern.ch
6+
7+
#include "MIDdetector.hh"
8+
#include "TDatabasePDG.h"
9+
#include "THnSparse.h"
10+
#include "TRandom.h"
11+
#include "TDatime.h"
12+
#include "TFile.h"
13+
#include "TVector3.h"
14+
15+
namespace o2 {
16+
namespace delphes {
17+
18+
//==========================================================================================================
19+
20+
bool MIDdetector::setup(const Char_t *nameInputFile = "muonAccEffPID.root") {
21+
22+
TDatime t;
23+
gRandom->SetSeed(t.GetDate()+t.GetYear()*t.GetHour()*t.GetMinute()*t.GetSecond());
24+
25+
mFileAccEffMuonPID = new TFile(nameInputFile);
26+
if (!mFileAccEffMuonPID) {
27+
printf("File %s not found\n",nameInputFile);
28+
return kFALSE;
29+
}
30+
if (!(mFileAccEffMuonPID->IsOpen())) {
31+
printf("File %s not open\n",nameInputFile);
32+
return kFALSE;
33+
}
34+
35+
for (Int_t iPart=kMuon; iPart<kNPart; iPart++) {
36+
mAccEffMuonPID[iPart] = (THnSparse*) mFileAccEffMuonPID->Get(Form("mAccEffMuonPID_%s",partLabel[iPart]));
37+
if (!mAccEffMuonPID[iPart]) {
38+
printf("Object %s not found, quitting\n",Form("mAccEffMuonPID_%s",partLabel[iPart]));
39+
return kFALSE;
40+
}
41+
}
42+
43+
printf("Setup of MIDdetector successfully completed\n");
44+
return kTRUE;
45+
46+
}
47+
48+
//==========================================================================================================
49+
50+
bool MIDdetector::hasMID(const Track &track) {
51+
52+
TVector3 v(track.XOuter, track.YOuter, track.ZOuter);
53+
return (TMath::Abs(v.Eta()) < mEtaMax);
54+
55+
}
56+
57+
//==========================================================================================================
58+
59+
bool MIDdetector::isMuon(const Track &track, int multiplicity=1) {
60+
61+
auto pdg = std::abs(track.PID);
62+
auto part = pidmap[pdg];
63+
if (part == kElectron) return kFALSE;
64+
65+
auto particle = (GenParticle*) track.Particle.GetObject();
66+
67+
Double_t var[4] = {track.P, track.Eta, particle->Z, double(multiplicity)};
68+
Double_t probMuonPID = mAccEffMuonPID[part]->GetBinContent(mAccEffMuonPID[part]->GetBin(var));
69+
return (gRandom->Uniform() < probMuonPID);
70+
71+
}
72+
73+
//==========================================================================================================
74+
75+
76+
} /** namespace delphes **/
77+
78+
} /** namespace o2 **/

0 commit comments

Comments
 (0)