Skip to content

Commit cb41ddc

Browse files
iouribelikovsawenzel
iouribelikov
authored and
sawenzel
committed
Adding some basic comparison at the outermost layer, and checking the availability of clusters
1 parent 96163fe commit cb41ddc

1 file changed

Lines changed: 131 additions & 54 deletions

File tree

‎Detectors/ITSMFT/ITS/macros/test/CheckTracks.C‎

Lines changed: 131 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -2,60 +2,88 @@
22
/// \brief Simple macro to check ITSU tracks
33

44
#if !defined(__CLING__) || defined(__ROOTCLING__)
5-
#include <array>
6-
7-
#include <TFile.h>
8-
#include <TTree.h>
9-
#include <TClonesArray.h>
10-
#include <TH2F.h>
11-
#include <TNtuple.h>
12-
#include <TCanvas.h>
13-
#include <TMath.h>
14-
#include <TString.h>
15-
16-
#include "SimulationDataFormat/MCTrack.h"
17-
#include "SimulationDataFormat/MCCompLabel.h"
18-
#include "SimulationDataFormat/MCTruthContainer.h"
5+
#include <array>
6+
7+
#include <TFile.h>
8+
#include <TTree.h>
9+
#include <TClonesArray.h>
10+
#include <TH2F.h>
11+
#include <TNtuple.h>
12+
#include <TCanvas.h>
13+
#include <TMath.h>
14+
#include <TString.h>
15+
16+
#include "SimulationDataFormat/TrackReference.h"
17+
#include "SimulationDataFormat/MCTrack.h"
18+
#include "SimulationDataFormat/MCCompLabel.h"
19+
#include "SimulationDataFormat/MCTruthContainer.h"
1920
#include "DataFormatsITSMFT/Cluster.h"
2021
#include "DataFormatsITS/TrackITS.h"
2122
#endif
2223

2324
void CheckTracks(Int_t nEvents = 10, TString mcEngine = "TGeant3") {
25+
using namespace o2::ITSMFT;
2426
using namespace o2::ITS;
2527

26-
TFile *f=TFile::Open("CheckTracks.root","recreate");
27-
TNtuple *nt=new TNtuple("ntt","track ntuple",
28-
"mcPhi:mcLam:mcPt:recPhi:recLam:recPt:ipD:ipZ:label");
28+
TFile* f = TFile::Open("CheckTracks.root", "recreate");
29+
TNtuple* nt = new TNtuple("ntt", "track ntuple",
30+
//"mcYOut:recYOut:"
31+
"mcZOut:recZOut:"
32+
"mcPhiOut:recPhiOut:"
33+
"mcThetaOut:recThetaOut:"
34+
"mcPhi:recPhi:"
35+
"mcLam:recLam:"
36+
"mcPt:recPt:"
37+
"ipD:ipZ:label");
2938

3039
char filename[100];
3140

3241
// MC tracks
3342
sprintf(filename, "AliceO2_%s.mc_%i_event.root", mcEngine.Data(), nEvents);
34-
TFile *file0 = TFile::Open(filename);
35-
TTree *mcTree=(TTree*)gFile->Get("o2sim");
43+
TFile* file0 = TFile::Open(filename);
44+
TTree* mcTree = (TTree*)gFile->Get("o2sim");
3645
std::vector<o2::MCTrack>* mcArr = nullptr;
37-
mcTree->SetBranchAddress("MCTrack",&mcArr);
46+
mcTree->SetBranchAddress("MCTrack", &mcArr);
47+
std::vector<o2::TrackReference>* mcTrackRefs = nullptr;
48+
mcTree->SetBranchAddress("TrackRefs", &mcTrackRefs);
49+
50+
// Clusters
51+
sprintf(filename, "AliceO2_%s.clus_%i_event.root", mcEngine.Data(), nEvents);
52+
TFile::Open(filename);
53+
TTree* clusTree = (TTree*)gFile->Get("o2sim");
54+
std::vector<Cluster>* clusArr = nullptr;
55+
auto* branch = clusTree->GetBranch("ITSCluster");
56+
if (!branch) {
57+
std::cout << "No clusters !" << std::endl;
58+
return;
59+
}
60+
branch->SetAddress(&clusArr);
61+
// Cluster MC labels
62+
o2::dataformats::MCTruthContainer<o2::MCCompLabel>* clusLabArr = nullptr;
63+
clusTree->SetBranchAddress("ITSClusterMCTruth", &clusLabArr);
3864

3965
// Reconstructed tracks
4066
sprintf(filename, "AliceO2_%s.trac_%i_event.root", mcEngine.Data(), nEvents);
41-
TFile *file1 = TFile::Open(filename);
42-
TTree *recTree=(TTree*)gFile->Get("o2sim");
67+
TFile* file1 = TFile::Open(filename);
68+
TTree* recTree = (TTree*)gFile->Get("o2sim");
4369
std::vector<TrackITS>* recArr = nullptr;
44-
recTree->SetBranchAddress("ITSTrack",&recArr);
70+
recTree->SetBranchAddress("ITSTrack", &recArr);
4571
// Track MC labels
46-
o2::dataformats::MCTruthContainer<o2::MCCompLabel> *trkLabArr=nullptr;
47-
recTree->SetBranchAddress("ITSTrackMCTruth",&trkLabArr);
48-
49-
Int_t tf=0, nrec=0;
50-
Int_t lastEventID=-1;
51-
Int_t nev=mcTree->GetEntries();
52-
for (Int_t n=0;n<nev;n++) {
53-
std::cout<<"\nMC event "<<n<<'/'<<nev<<std::endl;
54-
Int_t nGen=0, nGoo=0;
72+
o2::dataformats::MCTruthContainer<o2::MCCompLabel>* trkLabArr = nullptr;
73+
recTree->SetBranchAddress("ITSTrackMCTruth", &trkLabArr);
74+
75+
Int_t tf = 0, nrec = 0;
76+
Int_t lastEventID = -1;
77+
Int_t nev = mcTree->GetEntries();
78+
for (Int_t n = 0; n < nev; n++) {
79+
std::cout << "\nMC event " << n << '/' << nev << std::endl;
80+
Int_t nGen = 0, nGoo = 0;
5581
mcTree->GetEvent(n);
5682
Int_t nmc=mcArr->size();
83+
Int_t nmcrefs=mcTrackRefs->size();
5784

5885
while ((n>lastEventID) && (tf<recTree->GetEntries())) { // Cache a new reconstructed TF
86+
clusTree->GetEvent(tf);
5987
recTree->GetEvent(tf);
6088
nrec=recArr->size();
6189
for (int i=0; i<nrec; i++) { // Find the last MC event within this reconstructed TF
@@ -68,28 +96,48 @@ void CheckTracks(Int_t nEvents = 10, TString mcEngine = "TGeant3") {
6896
tf++;
6997
}
7098

71-
while(nmc--) {
99+
while (nmc--) {
72100
const auto& mcTrack = (*mcArr)[nmc];
73101
Int_t mID = mcTrack.getMotherTrackId();
74-
if (mID >= 0) continue;// Select primary particles
102+
if (mID >= 0)
103+
continue; // Select primary particles
75104
Int_t pdg = mcTrack.GetPdgCode();
76105
if (TMath::Abs(pdg) != 211) continue; // Select pions
77106

107+
int ok=0;
108+
//Check the availability of clusters
109+
for (int i=0; i < clusArr->size(); i++) {
110+
const Cluster &c = (*clusArr)[i];
111+
auto lab = (clusLabArr->getLabels(i))[0];
112+
if (lab.getEventID() != n) continue;
113+
if (lab.getTrackID() != nmc) continue;
114+
auto r= c.getX();
115+
if (TMath::Abs(r-2.2) <0.5) ok |= 0b1;
116+
if (TMath::Abs(r-3.0) <0.5) ok |= 0b10;
117+
if (TMath::Abs(r-3.8) <0.5) ok |= 0b100;
118+
if (TMath::Abs(r-19.5)<0.5) ok |= 0b1000;
119+
if (TMath::Abs(r-24.5)<0.5) ok |= 0b10000;
120+
if (TMath::Abs(r-34.5)<0.5) ok |= 0b100000;
121+
if (TMath::Abs(r-39.5)<0.5) ok |= 0b1000000;
122+
}
123+
if (ok != 0b1111111) continue;
124+
78125
nGen++;// Generated tracks for the efficiency calculation
79126

127+
//Float_t mcYOut=-1., recYOut=-1.;
128+
Float_t mcZOut=-1., recZOut=-1.;
129+
Float_t mcPhiOut=-1., recPhiOut=-1.;
130+
Float_t mcThetaOut=-1., recThetaOut=-1.;
80131
Float_t mcPx = mcTrack.GetStartVertexMomentumX();
81132
Float_t mcPy = mcTrack.GetStartVertexMomentumY();
82133
Float_t mcPz = mcTrack.GetStartVertexMomentumZ();
83-
Float_t mcPt = mcTrack.GetPt();
84-
Float_t mcPhi= TMath::ATan2(mcPy,mcPx);
85-
Float_t mcLam= TMath::ATan2(mcPz,mcPt);
86-
Float_t recPhi=-1.;
87-
Float_t recLam=-1.;
88-
Float_t recPt=-1.;
89-
Float_t ip[2]{0.,0.};
90-
Float_t label=-123456789.;
91-
92-
for (Int_t i=0;i<nrec;i++) {
134+
Float_t mcPhi = TMath::ATan2(mcPy, mcPx), recPhi = -1.;
135+
Float_t mcPt = mcTrack.GetPt(), recPt = -1.;
136+
Float_t mcLam = TMath::ATan2(mcPz, mcPt), recLam = -1.;
137+
Float_t ip[2]{ 0., 0. };
138+
Float_t label = -123456789.;
139+
140+
for (Int_t i = 0; i < nrec; i++) {
93141
const TrackITS& recTrack = (*recArr)[i];
94142
auto mclab = (trkLabArr->getLabels(i))[0];
95143
auto id = mclab.getEventID();
@@ -98,6 +146,25 @@ void CheckTracks(Int_t nEvents = 10, TString mcEngine = "TGeant3") {
98146
Int_t lab = mclab.getTrackID();
99147
if (TMath::Abs(lab) != nmc)
100148
continue;
149+
150+
for (auto& ref : *mcTrackRefs) {
151+
if (ref.getUserId() != 6)
152+
continue;
153+
if (ref.getTrackID() != nmc)
154+
continue;
155+
// mcYOut=ref.LocalY();
156+
mcZOut = ref.Z();
157+
mcPhiOut = ref.Phi();
158+
mcThetaOut = ref.Theta();
159+
break;
160+
}
161+
162+
auto out = recTrack.getParamOut();
163+
// recYOut = out.getY();
164+
recZOut = out.getZ();
165+
recPhiOut = out.getPhi();
166+
recThetaOut = out.getTheta();
167+
101168
std::array<float, 3> p;
102169
recTrack.getPxPyPzGlo(p);
103170
recPt = recTrack.getPt();
@@ -112,18 +179,28 @@ void CheckTracks(Int_t nEvents = 10, TString mcEngine = "TGeant3") {
112179
nGoo++; // Good found tracks for the efficiency calculation
113180
}
114181

115-
nt->Fill(mcPhi,mcLam,mcPt,recPhi,recLam,recPt,ip[0],ip[1],label);
116-
182+
nt->Fill( // mcYOut,recYOut,
183+
mcZOut, recZOut, mcPhiOut, recPhiOut, mcThetaOut, recThetaOut, mcPhi, recPhi, mcLam, recLam, mcPt, recPt, ip[0],
184+
ip[1], label);
117185
}
118-
Float_t eff = (nGen > 0) ? nGoo/Float_t(nGen) : -1.;
119-
std::cout<<"Good found tracks: "<<nGoo<<", efficiency: "<<eff<<std::endl;
186+
Float_t eff = (nGen > 0) ? nGoo / Float_t(nGen) : -1.;
187+
std::cout << "Good found tracks: " << nGoo << ", efficiency: " << eff << std::endl;
120188
}
121-
122-
// "recPt>0" means "found tracks only"
123-
// "label>0" means "found good tracks only"
124-
new TCanvas;nt->Draw("ipD","recPt>0 && label>0");
125-
new TCanvas;nt->Draw("mcLam-recLam","recPt>0 && label>0");
126-
new TCanvas;nt->Draw("mcPt-recPt","recPt>0 && label>0");
189+
190+
// "recPt>0" means "found tracks only"
191+
// "label>0" means "found good tracks only"
192+
new TCanvas;
193+
nt->Draw("ipD", "recPt>0 && label>0");
194+
new TCanvas;
195+
nt->Draw("mcLam-recLam", "recPt>0 && label>0");
196+
new TCanvas;
197+
nt->Draw("mcPt-recPt", "recPt>0 && label>0");
198+
new TCanvas;
199+
nt->Draw("mcZOut-recZOut", "recPt>0 && label>0 && abs(mcZOut-recZOut)<0.025");
200+
new TCanvas;
201+
nt->Draw("mcPhiOut-recPhiOut", "recPt>0 && label>0");
202+
new TCanvas;
203+
nt->Draw("mcThetaOut-recThetaOut", "recPt>0 && label>0");
127204
f->Write();
128205
f->Close();
129206
}

0 commit comments

Comments
 (0)