1+ // Copyright 2019-2026 CERN and copyright holders of ALICE O2.
2+ // See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3+ // All rights not expressly granted are reserved.
4+ //
5+ // This software is distributed under the terms of the GNU General Public
6+ // License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7+ //
8+ // In applying this license CERN does not waive the privileges and immunities
9+ // granted to it by virtue of its status as an Intergovernmental Organization
10+ // or submit itself to any jurisdiction.
11+
12+ /// \file postClusterSizeVsEta.C
13+ /// \brief A post-processing macro to draw average cluster size vs eta
14+
15+ #if !defined(__CLING__ ) || defined(__ROOTCLING__ )
16+ #include <iostream>
17+ #include <TCanvas.h>
18+ #include <TFile.h>
19+ #include <TH1F.h>
20+ #include <TH2F.h>
21+ #include <TNtuple.h>
22+ #include <TString.h>
23+ #include <TTree.h>
24+ #include <TROOT.h>
25+ #include <TStyle.h>
26+ #include <TLegend.h>
27+ #include <TProfile.h>
28+ #endif
29+
30+ using namespace std ;
31+
32+ // ### required input file: CheckClusters.root, which is the output of CheckClusters.C macro
33+ void postClusterSizeVsEta (const std ::string & strFileInput = "CheckClusters.root" )
34+ {
35+ gStyle -> SetOptStat (0 );
36+
37+ TFile * fileInput = new TFile (strFileInput .c_str ());
38+ TTree * tree = (TTree * )fileInput -> Get ("ntc" );
39+ std ::cout << "Opened tree: " << tree -> GetName () << ", entries = " << tree -> GetEntries () << std ::endl ;
40+
41+ // set branch addresses
42+ Float_t event ;
43+ Float_t mcTrackID ;
44+ Float_t hitLocX , hitLocZ ;
45+ Float_t hitGlobX , hitGlobY , hitGlobZ ;
46+ Float_t clusGlobX , clusGlobY , clusGlobZ ;
47+ Float_t clusLocX , clusLocZ ;
48+ Float_t rofFrame ;
49+ Float_t clusSize ;
50+ Float_t chipID ;
51+ Float_t layer ;
52+ Float_t disk ;
53+ Float_t subdet ;
54+ Float_t row , col ;
55+ Float_t pt ;
56+
57+ // set branch addresses
58+ tree -> SetBranchAddress ("event" , & event );
59+ tree -> SetBranchAddress ("mcTrackID" , & mcTrackID );
60+ tree -> SetBranchAddress ("hitLocX" , & hitLocX );
61+ tree -> SetBranchAddress ("hitLocZ" , & hitLocZ );
62+ tree -> SetBranchAddress ("hitGlobX" , & hitGlobX );
63+ tree -> SetBranchAddress ("hitGlobY" , & hitGlobY );
64+ tree -> SetBranchAddress ("hitGlobZ" , & hitGlobZ );
65+ tree -> SetBranchAddress ("clusGlobX" , & clusGlobX );
66+ tree -> SetBranchAddress ("clusGlobY" , & clusGlobY );
67+ tree -> SetBranchAddress ("clusGlobZ" , & clusGlobZ );
68+ tree -> SetBranchAddress ("clusLocX" , & clusLocX );
69+ tree -> SetBranchAddress ("clusLocZ" , & clusLocZ );
70+ tree -> SetBranchAddress ("rofFrame" , & rofFrame );
71+ tree -> SetBranchAddress ("clusSize" , & clusSize );
72+ tree -> SetBranchAddress ("chipID" , & chipID );
73+ tree -> SetBranchAddress ("layer" , & layer );
74+ tree -> SetBranchAddress ("disk" , & disk );
75+ tree -> SetBranchAddress ("subdet" , & subdet );
76+ tree -> SetBranchAddress ("row" , & row );
77+ tree -> SetBranchAddress ("col" , & col );
78+ tree -> SetBranchAddress ("pt" , & pt );
79+
80+ // Some QA histograms
81+ TH1F * hPt = new TH1F ("hPt" , "p_{T};p_{T};Entries" , 100 , 0. , 10. );
82+ TH1F * hClusSize = new TH1F ("hClusSize" , "Cluster size;clusSize;Entries" , 20 , 0. , 20. );
83+ TH1F * hLayer = new TH1F ("hLayer" , "Layer;layer;Entries" , 20 , -0.5 , 19.5 );
84+ TH1F * hDxGlob = new TH1F ("hDxGlob" , "clusGlobX - hitGlobX;#DeltaX [global];Entries" , 200 , -1. , 1. );
85+ TH1F * hDzGlob = new TH1F ("hDzGlob" , "clusGlobZ - hitGlobZ;#DeltaZ [global];Entries" , 200 , -1. , 1. );
86+ TH2F * hHitXY = new TH2F ("hHitXY" , "Hit global XY;hitGlobX;hitGlobY" , 200 , -20. , 20. , 200 , -20. , 20. );
87+ TH2F * hClusVsHitX = new TH2F ("hClusVsHitX" , "clusGlobX vs hitGlobX;hitGlobX;clusGlobX" , 200 , -20. , 20. , 200 , -20. , 20. );
88+
89+ // histograms for cluster size vs eta for each barrel layer:
90+ const int nLayers = 11 ;
91+ TH2F * hClustSizePerLayerVsEta [nLayers ];
92+ for (int i = 0 ; i < nLayers ; i ++ ) {
93+ hClustSizePerLayerVsEta [i ] = new TH2F (Form ("hClustSizePerLayerVsEta_Lay%d" , i ), Form ("Cluster size vs eta for layer %d;#eta;Cluster size" , i ), 200 , -5 , 5 , 101 , -0.5 , 100.5 );
94+ }
95+
96+ // Loop over entries
97+ const Long64_t nEntries = tree -> GetEntries ();
98+ for (Long64_t i = 0 ; i < nEntries ; ++ i ) {
99+ tree -> GetEntry (i );
100+
101+ // Fill QA histograms
102+ float dXGlob = clusGlobX - hitGlobX ;
103+ float dZGlob = clusGlobZ - hitGlobZ ;
104+ hPt -> Fill (pt );
105+ hClusSize -> Fill (clusSize );
106+ hLayer -> Fill (layer );
107+ hDxGlob -> Fill (dXGlob );
108+ hDzGlob -> Fill (dZGlob );
109+ hHitXY -> Fill (hitGlobX , hitGlobY );
110+ hClusVsHitX -> Fill (hitGlobX , clusGlobX );
111+
112+ // cls size vs eta:
113+ float clustR = sqrt (clusGlobX * clusGlobX + clusGlobY * clusGlobY );
114+ float clustPhi = atan2 (clusGlobY , clusGlobX );
115+ float clustTheta = atan2 (clustR , clusGlobZ );
116+ float clustEta = - log (tan (clustTheta / 2 ));
117+
118+ // !!! important: to avoid VD layers (numeration for ML starts from 0, while VD layers are also numbered as 0,1,2)
119+ if (clustR > 5 ) // cm
120+ hClustSizePerLayerVsEta [(int )layer + 3 ]-> Fill (clustEta , clusSize );
121+ else if (layer < 3 ) // VD layers
122+ hClustSizePerLayerVsEta [(int )layer ]-> Fill (clustEta , clusSize );
123+
124+ // progress print
125+ if ((i + 1 ) % 200000 == 0 ) {
126+ std ::cout << "Processed " << (i + 1 ) << " / " << nEntries << " entries" << std ::endl ;
127+ }
128+ }
129+
130+ // Save histograms to file
131+ TFile * fout = TFile ::Open ("clusterSizes_vs_eta.root" , "RECREATE" );
132+ hPt -> Write ();
133+ hClusSize -> Write ();
134+ hLayer -> Write ();
135+ hDxGlob -> Write ();
136+ hDzGlob -> Write ();
137+ hHitXY -> Write ();
138+ hClusVsHitX -> Write ();
139+
140+ // draw some QA histograms
141+ TCanvas * c1 = new TCanvas ("canv_clusters_QA" , "Clusters QA" , 1200 , 800 );
142+ c1 -> Divide (2 , 2 );
143+ c1 -> cd (1 );
144+ hPt -> Draw ();
145+ c1 -> cd (2 );
146+ hClusSize -> Draw ();
147+ c1 -> cd (3 );
148+ hDxGlob -> Draw ();
149+ c1 -> cd (4 );
150+ hHitXY -> Draw ("COLZ" );
151+
152+ int colors [] = {kRed , kBlue + 1 , kMagenta + 1 ,
153+ kRed , kBlue + 1 , kMagenta + 1 ,
154+ kCyan + 1 , kGray + 2 , kRed , kBlue , kMagenta + 1 , kCyan , kAzure + 1 , kOrange - 9 , kRed + 2 , kBlue + 2 , kMagenta + 2 };
155+
156+ TCanvas * canv_clsSize_vs_eta [nLayers ];
157+ TProfile * profPerLayerVsEta [nLayers ];
158+ for (int i = 0 ; i < nLayers ; i ++ ) {
159+ canv_clsSize_vs_eta [i ] = new TCanvas (Form ("canv_clsSize_vs_eta_Lay%d" , i ), Form ("Cluster size vs eta for layer %d" , i ), 800 , 600 );
160+ hClustSizePerLayerVsEta [i ]-> Draw ("COLZ" );
161+ gPad -> SetLogz ();
162+ profPerLayerVsEta [i ] = hClustSizePerLayerVsEta [i ]-> ProfileX ();
163+ profPerLayerVsEta [i ]-> SetLineColor (colors [i ]);
164+ profPerLayerVsEta [i ]-> SetMarkerColor (colors [i ]);
165+ profPerLayerVsEta [i ]-> SetMarkerStyle (i < 8 ? 20 : 24 );
166+ profPerLayerVsEta [i ]-> SetTitle (";#eta;#LTcluster size#GT" );
167+ profPerLayerVsEta [i ]-> DrawCopy ("same" );
168+
169+ hClustSizePerLayerVsEta [i ]-> Write ();
170+ profPerLayerVsEta [i ]-> Write ();
171+ }
172+
173+ // ### canvas with profiles for 3 VD layers
174+ TCanvas * canv_av_clsSize_vs_eta_VD_layers = new TCanvas ("canv_clsSize_vs_eta_VD_layers" , "Cluster size vs eta for VD layers" , 800 , 600 );
175+ TLegend * legLayersVD = new TLegend (0.3 , 0.72 , 0.65 , 0.89 );
176+ for (int i = 0 ; i < 3 ; i ++ ) {
177+ profPerLayerVsEta [i ]-> GetYaxis ()-> SetRangeUser (0. , 60. );
178+ profPerLayerVsEta [i ]-> DrawCopy (i == 0 ? "P" : "P same" );
179+ legLayersVD -> AddEntry (profPerLayerVsEta [i ], Form ("VD layer %d" , i ), "P" );
180+ }
181+ legLayersVD -> Draw ();
182+ gPad -> SetGrid ();
183+ canv_av_clsSize_vs_eta_VD_layers -> SaveAs ("clsSize_vs_eta_VD_layers.png" );
184+ canv_av_clsSize_vs_eta_VD_layers -> Write ();
185+
186+ // ### canvas with profiles for MLOT layers
187+ TCanvas * canv_av_clsSize_vs_eta_MLOT_layers = new TCanvas ("canv_clsSize_vs_eta_MLOT_layers" , "Cluster size vs eta for MLOT layers" , 800 , 600 );
188+ TLegend * legLayersMLOT = new TLegend (0.3 , 0.52 , 0.65 , 0.89 );
189+ for (int i = 3 ; i < nLayers ; i ++ ) {
190+ profPerLayerVsEta [i ]-> GetYaxis ()-> SetRangeUser (0. , 12.5 );
191+ profPerLayerVsEta [i ]-> GetXaxis ()-> SetRangeUser (-3.5 , 3.5 );
192+ profPerLayerVsEta [i ]-> DrawCopy (i == 3 ? "P" : "P same" );
193+ legLayersMLOT -> AddEntry (profPerLayerVsEta [i ], Form ("MLOT layer %d" , i ), "P" );
194+ }
195+ legLayersMLOT -> Draw ();
196+ gPad -> SetGrid ();
197+ canv_av_clsSize_vs_eta_MLOT_layers -> SaveAs ("clsSize_vs_eta_MLOT_layers.png" );
198+ canv_av_clsSize_vs_eta_MLOT_layers -> Write ();
199+ }
0 commit comments