Project
Loading...
Searching...
No Matches
ImpactParameter.cxx
Go to the documentation of this file.
1// Copyright 2019-2020 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// Skeleton derived from RS's code in ITSOffStudy
13
14#include <set>
15#include <vector>
16#include <utility>
17#include <string>
18#include "ITSStudies/ImpactParameter.h"
19#include "ITSStudies/TrackCuts.h"
20#include "ITSStudies/Helpers.h"
21#include "Framework/CCDBParamSpec.h"
22#include "Framework/Task.h"
23#include "DetectorsCommonDataFormats/DetID.h"
24#include "ReconstructionDataFormats/VtxTrackRef.h"
25#include "ReconstructionDataFormats/PrimaryVertex.h"
26#include "DataFormatsGlobalTracking/RecoContainer.h"
27#include "DetectorsVertexing/PVertexer.h"
28#include "DetectorsBase/Propagator.h"
29#include "CommonUtils/TreeStreamRedirector.h"
30#include "DetectorsBase/GRPGeomHelper.h"
31#include "DataFormatsParameters/GRPECSObject.h"
32#include "ITSMFTBase/DPLAlpideParam.h"
33#include "DetectorsCommonDataFormats/DetID.h"
34#include "Framework/DeviceSpec.h"
35#include "CommonUtils/ConfigurableParam.h"
36#include "Framework/CompletionPolicy.h"
37#include "Framework/ConfigParamSpec.h"
38#include "Framework/CompletionPolicyHelpers.h"
39#include "Framework/CallbacksPolicy.h"
40#include "Framework/ConfigParamRegistry.h"
41#include <TH1F.h>
42#include <TH2F.h>
43#include <TMath.h>
44#include <TROOT.h>
45#include <TSystem.h>
46#include <TString.h>
47#include "TGeoGlobalMagField.h"
48
49namespace o2
50{
51namespace its
52{
53namespace study
54{
55using namespace o2::framework;
56using namespace o2::globaltracking;
57using namespace o2::its::study;
58
59using DetID = o2::detectors::DetID;
60using PVertex = o2::dataformats::PrimaryVertex;
61using GTrackID = o2::dataformats::GlobalTrackID;
62
63class ImpactParameterStudy : public Task
64{
65 public:
66 ImpactParameterStudy(std::shared_ptr<DataRequest> dr,
67 std::shared_ptr<o2::base::GRPGeomRequest> gr,
68 mask_t src,
69 bool useMC) : mDataRequest(dr),
70 mGGCCDBRequest(gr),
71 mTracksSrc(src),
72 mUseMC(useMC) {}
73 ~ImpactParameterStudy() final = default;
74 void init(InitContext& ic) final;
75 void run(ProcessingContext&) final;
76 void endOfStream(EndOfStreamContext&) final;
77 void finaliseCCDB(ConcreteDataMatcher&, void*) final;
78 void process(o2::globaltracking::RecoContainer&);
79
80 private:
81 void updateTimeDependentParams(ProcessingContext& pc);
82 void saveHistograms();
83 void plotHistograms();
84 std::shared_ptr<o2::base::GRPGeomRequest> mGGCCDBRequest;
85 GTrackID::mask_t mTracksSrc{};
86 bool mUseMC{false};
87 o2::vertexing::PVertexer mVertexer;
88 float mITSROFrameLengthMUS = 0.;
89 float mITSROFBiasMUS = 0.;
90 std::unique_ptr<o2::utils::TreeStreamRedirector> mDBGOut;
91 // output histograms
92 std::unique_ptr<TH1F> mHistoContributorsPV{};
93 std::unique_ptr<TH1F> mHistoTrackType;
94 std::unique_ptr<TH1F> mHistoTrackTypeRej;
95 std::unique_ptr<TH1F> mHistoTrackTypeAcc;
96 std::unique_ptr<TH2F> mHistoXPvVsRefitted{};
97 std::unique_ptr<TH2F> mHistoYPvVsRefitted{};
98 std::unique_ptr<TH2F> mHistoZPvVsRefitted{};
99 std::unique_ptr<TH1F> mHistoXDeltaPVrefit{};
100 std::unique_ptr<TH1F> mHistoYDeltaPVrefit{};
101 std::unique_ptr<TH1F> mHistoZDeltaPVrefit{};
102 std::unique_ptr<TH2F> mHistoImpParXy{};
103 std::unique_ptr<TH2F> mHistoImpParZ{};
104 std::unique_ptr<TH2F> mHistoImpParXyPhi{};
105 std::unique_ptr<TH2F> mHistoImpParZPhi{};
106 std::unique_ptr<TH2F> mHistoImpParXyTop{};
107 std::unique_ptr<TH2F> mHistoImpParZTop{};
108 std::unique_ptr<TH2F> mHistoImpParXyBottom{};
109 std::unique_ptr<TH2F> mHistoImpParZBottom{};
110 std::unique_ptr<TH2F> mHistoImpParXyPositiveCharge{};
111 std::unique_ptr<TH2F> mHistoImpParZPositiveCharge{};
112 std::unique_ptr<TH2F> mHistoImpParXyNegativeCharge{};
113 std::unique_ptr<TH2F> mHistoImpParZNegativeCharge{};
114 std::unique_ptr<TH1F> mHistoImpParXyMeanPhi{};
115 std::unique_ptr<TH1F> mHistoImpParZMeanPhi{};
116 std::unique_ptr<TH1F> mHistoImpParXySigmaPhi{};
117 std::unique_ptr<TH1F> mHistoImpParZSigmaPhi{};
118 std::unique_ptr<TH1F> mHistoImpParXySigma{};
119 std::unique_ptr<TH1F> mHistoImpParZSigma{};
120 std::unique_ptr<TH1F> mHistoImpParXySigmaTop{};
121 std::unique_ptr<TH1F> mHistoImpParZSigmaTop{};
122 std::unique_ptr<TH1F> mHistoImpParXySigmaBottom{};
123 std::unique_ptr<TH1F> mHistoImpParZSigmaBottom{};
124 std::unique_ptr<TH1F> mHistoImpParXyMeanTop{};
125 std::unique_ptr<TH1F> mHistoImpParZMeanTop{};
126 std::unique_ptr<TH1F> mHistoImpParXyMeanBottom{};
127 std::unique_ptr<TH1F> mHistoImpParZMeanBottom{};
128 std::unique_ptr<TH1F> mHistoImpParXySigmaPositiveCharge{};
129 std::unique_ptr<TH1F> mHistoImpParZSigmaPositiveCharge{};
130 std::unique_ptr<TH1F> mHistoImpParXySigmaNegativeCharge{};
131 std::unique_ptr<TH1F> mHistoImpParZSigmaNegativeCharge{};
132
133 // output file
134 TString mOutName{};
135
136 // Data
137 std::shared_ptr<DataRequest> mDataRequest;
138 gsl::span<const PVertex> mPVertices;
139};
140
141void ImpactParameterStudy::init(InitContext& ic)
142{
143 o2::base::GRPGeomHelper::instance().setRequest(mGGCCDBRequest);
144 auto& params = ITSImpactParameterParamConfig::Instance();
145 mOutName = params.outFileName;
146 mDBGOut = std::make_unique<o2::utils::TreeStreamRedirector>(mOutName.Data(), "recreate");
147
148 std::vector<double> logPtBinning = helpers::makeLogBinning(100, 0.1, 10);
149 mHistoTrackType = std::make_unique<TH1F>("trackType", "# Track Type", 32, -0.5, 31.5);
150 mHistoTrackTypeRej = std::make_unique<TH1F>("trackTypeRej", "# Rejected Track Type", 32, -0.5, 31.5);
151 mHistoTrackTypeAcc = std::make_unique<TH1F>("trackTypeAcc", "# Filtered Track Type", 32, -0.5, 31.5);
152 mHistoContributorsPV = std::make_unique<TH1F>("nContribPVrefit", "# Contributors per PV", 100, 0, 100);
153 mHistoXPvVsRefitted = std::make_unique<TH2F>("histo2dXPvVsPVrefit", "#X PV vs PV_{-1}, #mum", 100, -10, 10, 100, -10, 10);
154 mHistoYPvVsRefitted = std::make_unique<TH2F>("histo2dYPvVsPVrefit", "#Y PV vs PV_{-1}, #mum", 100, -10, 10, 100, -10, 10);
155 mHistoZPvVsRefitted = std::make_unique<TH2F>("histo2dZPvVsPVrefit", "#Z PV vs PV_{-1}, #mum", 100, -10, 10, 100, -10, 10);
156 mHistoXDeltaPVrefit = std::make_unique<TH1F>("histoDeltaXPVrefit", "#DeltaX (PV-PV_{-1}), #mum", 300, -15, 15);
157 mHistoYDeltaPVrefit = std::make_unique<TH1F>("histoDeltaYPVrefit", "#DeltaY (PV-PV_{-1}), #mum", 300, -15, 15);
158 mHistoZDeltaPVrefit = std::make_unique<TH1F>("histoDeltaZPVrefit", "#DeltaZ (PV-PV_{-1}), #mum", 300, -15, 15);
159 mHistoImpParXyPhi = std::make_unique<TH2F>("histoImpParXyPhi", "#Phi; #phi; Impact Parameter XY (#mum);", 100, 0., 6.28, 100, -1000, 1000);
160 mHistoImpParZPhi = std::make_unique<TH2F>("histoImpParZPhi", "#Phi; #phi; Impact Parameter Z (#mum);", 100, 0., 6.28, 100, -1000, 1000);
161 mHistoImpParZ = std::make_unique<TH2F>("histoImpParZ", "Impact Parameter Z; #it{p}_{T} (GeV/#it{c}); #mum", logPtBinning.size() - 1, logPtBinning.data(), 100, -1000, 1000);
162 mHistoImpParXy = std::make_unique<TH2F>("histoImpParXy", "Impact Parameter XY; #it{p}_{T} (GeV/#it{c}); #mum", logPtBinning.size() - 1, logPtBinning.data(), 100, -1000, 1000);
163 mHistoImpParXyTop = std::make_unique<TH2F>("histoImpParXyTop", "Impact Parameter XY, #phi(track)<#pi; #it{p}_{T} (GeV/#it{c}); #mum", logPtBinning.size() - 1, logPtBinning.data(), 100, -1000, 1000);
164 mHistoImpParXyBottom = std::make_unique<TH2F>("histoImpParXyBottom", "Impact Parameter XY, #phi(track)>#pi; #it{p}_{T} (GeV/#it{c}); #mum", logPtBinning.size() - 1, logPtBinning.data(), 100, -1000, 1000);
165 mHistoImpParZTop = std::make_unique<TH2F>("histoImpParZTop", "Impact Parameter Z, #phi(track)<#pi; #it{p}_{T} (GeV/#it{c}); #mum", logPtBinning.size() - 1, logPtBinning.data(), 100, -1000, 1000);
166 mHistoImpParZBottom = std::make_unique<TH2F>("histoImpParZBottom", "Impact Parameter Z, #phi(track)>#pi; #it{p}_{T} (GeV/#it{c}); #mum", logPtBinning.size() - 1, logPtBinning.data(), 100, -1000, 1000);
167 mHistoImpParXyNegativeCharge = std::make_unique<TH2F>("histoImpParXyNegativeCharge", "Impact Parameter XY, sign<0; #it{p}_{T} (GeV/#it{c}); #mum", logPtBinning.size() - 1, logPtBinning.data(), 100, -1000, 1000);
168 mHistoImpParXyPositiveCharge = std::make_unique<TH2F>("histoImpParXyPositiveCharge", "Impact Parameter XY, sign>0; #it{p}_{T} (GeV/#it{c}); #mum", logPtBinning.size() - 1, logPtBinning.data(), 100, -1000, 1000);
169 mHistoImpParZNegativeCharge = std::make_unique<TH2F>("histoImpParZNegativeCharge", "Impact Parameter Z, sign<0; #it{p}_{T} (GeV/#it{c}); #mum", logPtBinning.size() - 1, logPtBinning.data(), 100, -1000, 1000);
170 mHistoImpParZPositiveCharge = std::make_unique<TH2F>("histoImpParZPositiveCharge", "Impact Parameter Z, sign>0; #it{p}_{T} (GeV/#it{c}); #mum", logPtBinning.size() - 1, logPtBinning.data(), 100, -1000, 1000);
171
172 mHistoImpParXySigma = std::make_unique<TH1F>("histoImpParXySigma", "Pointing Resolution XY; #it{p}_{T} (GeV/#it{c}); #mum", logPtBinning.size() - 1, logPtBinning.data());
173 mHistoImpParZSigma = std::make_unique<TH1F>("histoImpParZSigma", "Pointing Resolution Z; #it{p}_{T} (GeV/#it{c}); #mum", logPtBinning.size() - 1, logPtBinning.data());
174 mHistoImpParXyMeanPhi = std::make_unique<TH1F>("histoImpParXyMean", "Pointing Resolution XY; #phi; Mean #mum", 100, 0., 6.28);
175 mHistoImpParZMeanPhi = std::make_unique<TH1F>("histoImpParZMean", "Pointing Resolution Z; #phi; Mean #mum", 100, 0., 6.28);
176 mHistoImpParXySigmaPhi = std::make_unique<TH1F>("histoImpParXySigmaPhi", "Pointing Resolution XY; #phi; #sigma #mum", 100, 0., 6.28);
177 mHistoImpParZSigmaPhi = std::make_unique<TH1F>("histoImpParZSigmaPhi", "Pointing Resolution Z; #phi; #sigma #mum", 100, 0., 6.28);
178 mHistoImpParXySigmaTop = std::make_unique<TH1F>("histoImpParXySigmaTop", "Pointing Resolution XY, Top; #it{p}_{T} (GeV/#it{c}); #mum", logPtBinning.size() - 1, logPtBinning.data());
179 mHistoImpParZSigmaTop = std::make_unique<TH1F>("histoImpParZSigmaTop", "Pointing Resolution Z, Top; #it{p}_{T} (GeV/#it{c}); #mum", logPtBinning.size() - 1, logPtBinning.data());
180 mHistoImpParXySigmaBottom = std::make_unique<TH1F>("histoImpParXySigmaBottom", "Pointing Resolution XY, Bottom; #it{p}_{T} (GeV/#it{c}); #mum", logPtBinning.size() - 1, logPtBinning.data());
181 mHistoImpParZSigmaBottom = std::make_unique<TH1F>("histoImpParZSigmaBottom", "Pointing Resolution Z, Bottom; #it{p}_{T} (GeV/#it{c}); #mum", logPtBinning.size() - 1, logPtBinning.data());
182 mHistoImpParXyMeanTop = std::make_unique<TH1F>("histoImpParXyMeanTop", "Pointing Resolution XY, Top; #it{p}_{T} (GeV/#it{c}); #mum", logPtBinning.size() - 1, logPtBinning.data());
183 mHistoImpParZMeanTop = std::make_unique<TH1F>("histoImpParZMeanTop", "Pointing Resolution Z, Top; #it{p}_{T} (GeV/#it{c}); #mum", logPtBinning.size() - 1, logPtBinning.data());
184 mHistoImpParXyMeanBottom = std::make_unique<TH1F>("histoImpParXyMeanBottom", "Mean Pointing Resolution XY, Bottom; #it{p}_{T} (GeV/#it{c}); #mum", logPtBinning.size() - 1, logPtBinning.data());
185 mHistoImpParZMeanBottom = std::make_unique<TH1F>("histoImpParZMeanBottom", "Mean Pointing Resolution Z, Bottom; #it{p}_{T} (GeV/#it{c}); #mum", logPtBinning.size() - 1, logPtBinning.data());
186 mHistoImpParXySigmaPositiveCharge = std::make_unique<TH1F>("histoImpParXySigmaPositiveCharge", "Pointing Resolution XY, sign>0; #it{p}_{T} (GeV/#it{c}); #mum", logPtBinning.size() - 1, logPtBinning.data());
187 mHistoImpParZSigmaPositiveCharge = std::make_unique<TH1F>("histoImpParZSigmaPositiveCharge", "Pointing Resolution Z, sign>0; #it{p}_{T} (GeV/#it{c}); #mum", logPtBinning.size() - 1, logPtBinning.data());
188 mHistoImpParXySigmaNegativeCharge = std::make_unique<TH1F>("histoImpParXySigmaNegativeCharge", "Pointing Resolution XY, sign<0; #it{p}_{T} (GeV/#it{c}); #mum", logPtBinning.size() - 1, logPtBinning.data());
189 mHistoImpParZSigmaNegativeCharge = std::make_unique<TH1F>("histoImpParZSigmaNegativeCharge", "Pointing Resolution Z, sign<0; #it{p}_{T} (GeV/#it{c}); #mum", logPtBinning.size() - 1, logPtBinning.data());
190
191 mHistoTrackType->SetDirectory(nullptr);
192 mHistoTrackTypeRej->SetDirectory(nullptr);
193 mHistoTrackTypeAcc->SetDirectory(nullptr);
194 mHistoContributorsPV->SetDirectory(nullptr);
195 mHistoXPvVsRefitted->SetDirectory(nullptr);
196 mHistoYPvVsRefitted->SetDirectory(nullptr);
197 mHistoZPvVsRefitted->SetDirectory(nullptr);
198 mHistoXDeltaPVrefit->SetDirectory(nullptr);
199 mHistoYDeltaPVrefit->SetDirectory(nullptr);
200 mHistoZDeltaPVrefit->SetDirectory(nullptr);
201 mHistoImpParZ->SetDirectory(nullptr);
202 mHistoImpParXy->SetDirectory(nullptr);
203 mHistoImpParXyPhi->SetDirectory(nullptr);
204 mHistoImpParZPhi->SetDirectory(nullptr);
205 mHistoImpParXyTop->SetDirectory(nullptr);
206 mHistoImpParXyBottom->SetDirectory(nullptr);
207 mHistoImpParZTop->SetDirectory(nullptr);
208 mHistoImpParZBottom->SetDirectory(nullptr);
209 mHistoImpParXyNegativeCharge->SetDirectory(nullptr);
210 mHistoImpParXyPositiveCharge->SetDirectory(nullptr);
211 mHistoImpParZNegativeCharge->SetDirectory(nullptr);
212 mHistoImpParZPositiveCharge->SetDirectory(nullptr);
213 mHistoImpParXyMeanPhi->SetDirectory(nullptr);
214 mHistoImpParZMeanPhi->SetDirectory(nullptr);
215 mHistoImpParXySigmaPhi->SetDirectory(nullptr);
216 mHistoImpParZSigmaPhi->SetDirectory(nullptr);
217 mHistoImpParXySigma->SetDirectory(nullptr);
218 mHistoImpParZSigma->SetDirectory(nullptr);
219 mHistoImpParXySigmaTop->SetDirectory(nullptr);
220 mHistoImpParZSigmaTop->SetDirectory(nullptr);
221 mHistoImpParXySigmaBottom->SetDirectory(nullptr);
222 mHistoImpParZSigmaBottom->SetDirectory(nullptr);
223 mHistoImpParXyMeanTop->SetDirectory(nullptr);
224 mHistoImpParZMeanTop->SetDirectory(nullptr);
225 mHistoImpParXyMeanBottom->SetDirectory(nullptr);
226 mHistoImpParZMeanBottom->SetDirectory(nullptr);
227 mHistoImpParXySigmaPositiveCharge->SetDirectory(nullptr);
228 mHistoImpParZSigmaPositiveCharge->SetDirectory(nullptr);
229 mHistoImpParXySigmaNegativeCharge->SetDirectory(nullptr);
230 mHistoImpParZSigmaNegativeCharge->SetDirectory(nullptr);
231}
232
233void ImpactParameterStudy::run(ProcessingContext& pc)
234{
235 o2::base::GRPGeomHelper::instance().checkUpdates(pc);
236 o2::globaltracking::RecoContainer recoData;
237 recoData.collectData(pc, *mDataRequest.get());
238 updateTimeDependentParams(pc);
239 process(recoData);
240}
241
242void ImpactParameterStudy::process(o2::globaltracking::RecoContainer& recoData)
243{
244 auto& params = ITSImpactParameterParamConfig::Instance();
245
246 o2::base::GRPGeomHelper::instance().getGRPMagField()->print();
247 o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrLUT;
248 std::vector<o2::track::TrackParCov> vecPvContributorTrackParCov;
249 std::vector<int64_t> vec_globID_contr = {};
250 std::vector<o2::track::TrackParCov> trueVecPvContributorTrackParCov;
251 std::vector<int64_t> trueVec_globID_contr = {};
252 float impParRPhi, impParZ;
253 constexpr float toMicrometers = 10000.f; // Conversion from [cm] to [mum]
254 auto trackIndex = recoData.getPrimaryVertexMatchedTracks(); // Global ID's for associated tracks
255 auto vtxRefs = recoData.getPrimaryVertexMatchedTrackRefs(); // references from vertex to these track IDs
256 auto pvertices = recoData.getPrimaryVertices();
257 auto trmTPCTracks = recoData.getTPCTracks();
258
259 TrackCuts cuts; // ITS and TPC(commented) cut implementation
260
261 int nv = vtxRefs.size() - 1; // The last entry is for unassigned tracks, ignore them
262 for (int iv = 0; iv < nv; iv++) { // Loop over PVs
263 LOGP(info, "*** NEW VERTEX {}***", iv);
264 const auto& vtref = vtxRefs[iv];
265 const o2::dataformats::VertexBase& pv = pvertices[iv];
266 int it = vtref.getFirstEntry(), itLim = it + vtref.getEntries();
267 int i = 0;
268
269 for (; it < itLim; it++) {
270 auto tvid = trackIndex[it];
271 // LOGP(info,"ORIGIN: {}", tvid.getSourceName(tvid.getSource()));
272 mHistoTrackType->Fill(tvid.getSource());
273 if (!recoData.isTrackSourceLoaded(tvid.getSource())) {
274 // LOGP(info,"SOURCE Rej: {}", tvid.getSourceName(tvid.getSource()));
275 mHistoTrackTypeRej->Fill(tvid.getSource());
276 continue;
277 }
278 // LOGP(info,"ORIGIN: {} INDEX: {}", tvid.getSourceName(tvid.getSource()), trackIndex[it]);
279 mHistoTrackTypeAcc->Fill(tvid.getSource());
280 const o2::track::TrackParCov& trc = recoData.getTrackParam(tvid); // The actual track
281
282 auto refs = recoData.getSingleDetectorRefs(tvid);
283 if (!refs[GTrackID::ITS].isIndexSet()) { // might be an afterburner track
284 // LOGP(info, " ** AFTERBURN **");
285 continue;
286 }
287 // Apply track selections
288 if (params.applyTrackCuts) {
289 if (!cuts.isSelected(tvid, recoData)) {
290 // LOGP(info, "Fail");
291 continue;
292 }
293 }
294 // Vectors for reconstructing the vertex
295 vec_globID_contr.push_back(trackIndex[it]);
296 vecPvContributorTrackParCov.push_back(trc);
297 // Store vector with index and tracks ITS only -- same order as the GLOBAL vectors
298 int itsTrackID = refs[GTrackID::ITS].getIndex();
299 const o2::track::TrackParCov& trcITS = recoData.getTrackParam(itsTrackID); // The actual ITS track
300 trueVec_globID_contr.push_back(itsTrackID);
301 trueVecPvContributorTrackParCov.push_back(trcITS);
302 // LOGP(info, "SOURCE: {} indexGLOBAL: {} indexITS: {}", tvid.getSourceName(tvid.getSource()), trackIndex[it], trueVec_globID_contr[i]);
303 i++;
304 } // end loop tracks
305 LOGP(info, "************ SIZE INDEX GLOBAL: {} ", vec_globID_contr.size());
306 LOGP(info, "************ SIZE INDEX ITS: {} ", trueVec_globID_contr.size());
307
308 it = vtref.getFirstEntry();
309 // Preparation PVertexer refit
310 o2::conf::ConfigurableParam::updateFromString("pvertexer.useMeanVertexConstraint=false");
311 bool PVrefit_doable = mVertexer.prepareVertexRefit(vecPvContributorTrackParCov, pv);
312 if (!PVrefit_doable) {
313 LOG(info) << "Not enough tracks accepted for the refit --> Skipping vertex";
314 } else {
315 mHistoContributorsPV->Fill(vecPvContributorTrackParCov.size());
316 // Neglect vertices with small number of contributors (configurable)
317 if (vecPvContributorTrackParCov.size() < params.minNumberOfContributors) {
318 continue;
319 }
320 for (it = 0; it < vec_globID_contr.size(); it++) {
321 // vector of booleans to keep track of the track to be skipped
322 std::vector<bool> vec_useTrk_PVrefit(vec_globID_contr.size(), true);
323 auto tvid = vec_globID_contr[it];
324 auto trackIterator = std::find(vec_globID_contr.begin(), vec_globID_contr.end(), vec_globID_contr[it]);
325 if (trackIterator != vec_globID_contr.end()) {
326 // LOGP(info,"************ Trackiterator: {} : {} ", it+1, vec_globID_contr[it]);
327 o2::dataformats::VertexBase PVbase_recalculated;
329 const int entry = std::distance(vec_globID_contr.begin(), trackIterator);
330 if (!params.useAllTracks) {
331 vec_useTrk_PVrefit[entry] = false;
332 }
333 auto Pvtx_refitted = mVertexer.refitVertex(vec_useTrk_PVrefit, pv); // vertex refit
334 // enable the dca recalculation for the current PV contributor, after removing it from the PV refit
335 bool recalc_imppar = true;
336 if (Pvtx_refitted.getChi2() < 0) {
337 // LOG(info) << "---> Refitted vertex has bad chi2 = " << Pvtx_refitted.getChi2();
338 recalc_imppar = false;
339 }
340 vec_useTrk_PVrefit[entry] = true;
341 if (recalc_imppar) {
342 const double DeltaX = pv.getX() - Pvtx_refitted.getX();
343 const double DeltaY = pv.getY() - Pvtx_refitted.getY();
344 const double DeltaZ = pv.getZ() - Pvtx_refitted.getZ();
345 mHistoXPvVsRefitted->Fill(pv.getX(), Pvtx_refitted.getX());
346 mHistoYPvVsRefitted->Fill(pv.getY(), Pvtx_refitted.getY());
347 mHistoZPvVsRefitted->Fill(pv.getZ(), Pvtx_refitted.getZ());
348 mHistoXDeltaPVrefit->Fill(DeltaX);
349 mHistoYDeltaPVrefit->Fill(DeltaY);
350 mHistoZDeltaPVrefit->Fill(DeltaZ);
351
352 // fill the newly calculated PV
353 PVbase_recalculated.setX(Pvtx_refitted.getX());
354 PVbase_recalculated.setY(Pvtx_refitted.getY());
355 PVbase_recalculated.setZ(Pvtx_refitted.getZ());
356 PVbase_recalculated.setCov(Pvtx_refitted.getSigmaX2(), Pvtx_refitted.getSigmaXY(), Pvtx_refitted.getSigmaY2(), Pvtx_refitted.getSigmaXZ(), Pvtx_refitted.getSigmaYZ(), Pvtx_refitted.getSigmaZ2());
357
358 auto trueID = trueVec_globID_contr[it];
359 const o2::track::TrackParCov& trc = recoData.getTrackParam(trueID);
360 auto pt = trc.getPt();
361 o2::gpu::gpustd::array<float, 2> dcaInfo{-999., -999.};
362 // LOGP(info, " ---> Bz={}", o2::base::Propagator::Instance()->getNominalBz());
363 o2::track::TrackPar trcTmp{trc};
364 if (o2::base::Propagator::Instance()->propagateToDCABxByBz({Pvtx_refitted.getX(), Pvtx_refitted.getY(), Pvtx_refitted.getZ()}, trcTmp, 2.f, matCorr, &dcaInfo)) {
365 impParRPhi = dcaInfo[0] * toMicrometers;
366 impParZ = dcaInfo[1] * toMicrometers;
367 mHistoImpParXy->Fill(pt, impParRPhi);
368 mHistoImpParZ->Fill(pt, impParZ);
369 double phi = trcTmp.getPhi();
370 mHistoImpParXyPhi->Fill(phi, impParRPhi);
371 mHistoImpParZPhi->Fill(phi, impParZ);
372 if (phi < TMath::Pi()) {
373 mHistoImpParXyTop->Fill(pt, impParRPhi);
374 mHistoImpParZTop->Fill(pt, impParZ);
375 }
376 if (phi > TMath::Pi()) {
377 mHistoImpParXyBottom->Fill(pt, impParRPhi);
378 mHistoImpParZBottom->Fill(pt, impParZ);
379 }
380 double sign = trcTmp.getSign();
381 if (sign < 0) {
382 mHistoImpParXyNegativeCharge->Fill(pt, impParRPhi);
383 mHistoImpParZNegativeCharge->Fill(pt, impParZ);
384 } else {
385 mHistoImpParXyPositiveCharge->Fill(pt, impParRPhi);
386 mHistoImpParZPositiveCharge->Fill(pt, impParZ);
387 }
388 }
389 } // end recalc impact param
390 }
391 } // end loop tracks in pv refitted
392 } // else pv refit duable */
393 vec_globID_contr.clear();
394 vecPvContributorTrackParCov.clear();
395 trueVec_globID_contr.clear();
396 trueVecPvContributorTrackParCov.clear();
397
398 } // end loop pv
399
400 mHistoImpParXy->FitSlicesY();
401 mHistoImpParZ->FitSlicesY();
402 mHistoImpParXyPhi->FitSlicesY();
403 mHistoImpParZPhi->FitSlicesY();
404 mHistoImpParXyTop->FitSlicesY();
405 mHistoImpParZTop->FitSlicesY();
406 mHistoImpParXyBottom->FitSlicesY();
407 mHistoImpParZBottom->FitSlicesY();
408 mHistoImpParXyNegativeCharge->FitSlicesY();
409 mHistoImpParZNegativeCharge->FitSlicesY();
410 mHistoImpParXyPositiveCharge->FitSlicesY();
411 mHistoImpParZPositiveCharge->FitSlicesY();
412 mHistoImpParXySigma = std::unique_ptr<TH1F>(static_cast<TH1F*>((gROOT->FindObject("histoImpParXy_2"))->Clone()));
413 mHistoImpParZSigma = std::unique_ptr<TH1F>(static_cast<TH1F*>((gROOT->FindObject("histoImpParZ_2"))->Clone()));
414 mHistoImpParXyMeanPhi = std::unique_ptr<TH1F>(static_cast<TH1F*>((gROOT->FindObject("histoImpParXyPhi_1"))->Clone()));
415 mHistoImpParZMeanPhi = std::unique_ptr<TH1F>(static_cast<TH1F*>((gROOT->FindObject("histoImpParZPhi_1"))->Clone()));
416 mHistoImpParXySigmaPhi = std::unique_ptr<TH1F>(static_cast<TH1F*>((gROOT->FindObject("histoImpParXyPhi_2"))->Clone()));
417 mHistoImpParZSigmaPhi = std::unique_ptr<TH1F>(static_cast<TH1F*>((gROOT->FindObject("histoImpParZPhi_2"))->Clone()));
418 mHistoImpParXyMeanTop = std::unique_ptr<TH1F>(static_cast<TH1F*>((gROOT->FindObject("histoImpParXyTop_1"))->Clone()));
419 mHistoImpParZMeanTop = std::unique_ptr<TH1F>(static_cast<TH1F*>((gROOT->FindObject("histoImpParZTop_1"))->Clone()));
420 mHistoImpParXyMeanBottom = std::unique_ptr<TH1F>(static_cast<TH1F*>((gROOT->FindObject("histoImpParXyBottom_1"))->Clone()));
421 mHistoImpParZMeanBottom = std::unique_ptr<TH1F>(static_cast<TH1F*>((gROOT->FindObject("histoImpParZBottom_1"))->Clone()));
422 mHistoImpParXySigmaTop = std::unique_ptr<TH1F>(static_cast<TH1F*>((gROOT->FindObject("histoImpParXyTop_2"))->Clone()));
423 mHistoImpParZSigmaTop = std::unique_ptr<TH1F>(static_cast<TH1F*>((gROOT->FindObject("histoImpParZTop_2"))->Clone()));
424 mHistoImpParXySigmaBottom = std::unique_ptr<TH1F>(static_cast<TH1F*>((gROOT->FindObject("histoImpParXyBottom_2"))->Clone()));
425 mHistoImpParZSigmaBottom = std::unique_ptr<TH1F>(static_cast<TH1F*>((gROOT->FindObject("histoImpParZBottom_2"))->Clone()));
426 mHistoImpParXySigmaNegativeCharge = std::unique_ptr<TH1F>(static_cast<TH1F*>((gROOT->FindObject("histoImpParXyNegativeCharge_2"))->Clone()));
427 mHistoImpParZSigmaNegativeCharge = std::unique_ptr<TH1F>(static_cast<TH1F*>((gROOT->FindObject("histoImpParZNegativeCharge_2"))->Clone()));
428 mHistoImpParXySigmaPositiveCharge = std::unique_ptr<TH1F>(static_cast<TH1F*>((gROOT->FindObject("histoImpParXyPositiveCharge_2"))->Clone()));
429 mHistoImpParZSigmaPositiveCharge = std::unique_ptr<TH1F>(static_cast<TH1F*>((gROOT->FindObject("histoImpParZPositiveCharge_2"))->Clone()));
430}
431// end process
432
433void ImpactParameterStudy::updateTimeDependentParams(ProcessingContext& pc)
434{
435 o2::base::GRPGeomHelper::instance().checkUpdates(pc);
436 static bool initOnceDone = false;
437 if (!initOnceDone) { // this params need to be queried only once
438 initOnceDone = true;
439 auto grp = o2::base::GRPGeomHelper::instance().getGRPECS();
440 mVertexer.init();
441 if (pc.services().get<const o2::framework::DeviceSpec>().inputTimesliceId == 0) {
442 }
443 }
444}
445
446void ImpactParameterStudy::saveHistograms()
447{
448 mDBGOut.reset();
449 TFile fout(mOutName.Data(), "update");
450 fout.WriteTObject(mHistoContributorsPV.get());
451 fout.WriteTObject(mHistoTrackType.get());
452 fout.WriteTObject(mHistoTrackTypeRej.get());
453 fout.WriteTObject(mHistoTrackTypeAcc.get());
454 fout.WriteTObject(mHistoXPvVsRefitted.get());
455 fout.WriteTObject(mHistoYPvVsRefitted.get());
456 fout.WriteTObject(mHistoZPvVsRefitted.get());
457 fout.WriteTObject(mHistoXDeltaPVrefit.get());
458 fout.WriteTObject(mHistoYDeltaPVrefit.get());
459 fout.WriteTObject(mHistoZDeltaPVrefit.get());
460 fout.WriteTObject(mHistoImpParZ.get());
461 fout.WriteTObject(mHistoImpParXy.get());
462 fout.WriteTObject(mHistoImpParXyPhi.get());
463 fout.WriteTObject(mHistoImpParZPhi.get());
464 fout.WriteTObject(mHistoImpParXyTop.get());
465 fout.WriteTObject(mHistoImpParZTop.get());
466 fout.WriteTObject(mHistoImpParXyBottom.get());
467 fout.WriteTObject(mHistoImpParZBottom.get());
468 fout.WriteTObject(mHistoImpParXyPositiveCharge.get());
469 fout.WriteTObject(mHistoImpParZPositiveCharge.get());
470 fout.WriteTObject(mHistoImpParXyNegativeCharge.get());
471 fout.WriteTObject(mHistoImpParZNegativeCharge.get());
472 fout.WriteTObject(mHistoImpParZMeanPhi.get());
473 fout.WriteTObject(mHistoImpParXyMeanPhi.get());
474 fout.WriteTObject(mHistoImpParZSigmaPhi.get());
475 fout.WriteTObject(mHistoImpParXySigmaPhi.get());
476 fout.WriteTObject(mHistoImpParZSigma.get());
477 fout.WriteTObject(mHistoImpParXySigma.get());
478 fout.WriteTObject(mHistoImpParZMeanTop.get());
479 fout.WriteTObject(mHistoImpParXyMeanTop.get());
480 fout.WriteTObject(mHistoImpParZMeanBottom.get());
481 fout.WriteTObject(mHistoImpParXyMeanBottom.get());
482 fout.WriteTObject(mHistoImpParZSigmaTop.get());
483 fout.WriteTObject(mHistoImpParXySigmaTop.get());
484 fout.WriteTObject(mHistoImpParZSigmaBottom.get());
485 fout.WriteTObject(mHistoImpParXySigmaBottom.get());
486 fout.WriteTObject(mHistoImpParZSigmaPositiveCharge.get());
487 fout.WriteTObject(mHistoImpParXySigmaPositiveCharge.get());
488 fout.WriteTObject(mHistoImpParZSigmaNegativeCharge.get());
489 fout.WriteTObject(mHistoImpParXySigmaNegativeCharge.get());
490 LOGP(info, "Impact Parameters histograms stored in {}", mOutName.Data());
491 fout.Close();
492}
493
494void ImpactParameterStudy::plotHistograms()
495{
496 TString directoryName = "./plotsImpactParameter";
497 gSystem->mkdir(directoryName);
498 TCanvas* dcaXyVsdcaZ = helpers::prepareSimpleCanvas2Histograms(*mHistoImpParXySigma, kRed, "#sigma_{XY}", *mHistoImpParZSigma, kBlue, "#sigma_{Z}");
499 TCanvas* dcaXyTopVsBottom = helpers::prepareSimpleCanvas2Histograms(*mHistoImpParXySigmaTop, kRed, "#sigma_{XY}^{TOP}", *mHistoImpParXySigmaBottom, kBlue, "#sigma_{XY}^{BOTTOM}");
500 TCanvas* dcaZTopVsBottom = helpers::prepareSimpleCanvas2Histograms(*mHistoImpParZSigmaTop, kRed, "#sigma_{Z}^{TOP}", *mHistoImpParZSigmaBottom, kBlue, "#sigma_{Z}^{BOTTOM}");
501 TCanvas* dcaXyPosVsNeg = helpers::prepareSimpleCanvas2Histograms(*mHistoImpParXySigmaPositiveCharge, kRed, "#sigma_{XY}^{positive}", *mHistoImpParXySigmaNegativeCharge, kBlue, "#sigma_{XY}^{negative}");
502 TCanvas* dcaZPosVsNeg = helpers::prepareSimpleCanvas2Histograms(*mHistoImpParZSigmaPositiveCharge, kRed, "#sigma_{Z}^{positive}", *mHistoImpParZSigmaNegativeCharge, kBlue, "#sigma_{Z}^{negative}");
503 TCanvas* dcaXYvsPhi = helpers::plot2DwithMeanAndSigma(*mHistoImpParXyPhi, *mHistoImpParXyMeanPhi, *mHistoImpParXySigmaPhi, kRed);
504 TCanvas* dcaZvsPhi = helpers::plot2DwithMeanAndSigma(*mHistoImpParZPhi, *mHistoImpParZMeanPhi, *mHistoImpParZSigmaPhi, kRed);
505 TCanvas* dcaXyPhiMeanTopVsBottom = helpers::prepareSimpleCanvas2DcaMeanValues(*mHistoImpParXyMeanTop, kRed, "mean_{XY}^{TOP}", *mHistoImpParXyMeanBottom, kBlue, "mean_{XY}^{BOTTOM}");
506 TCanvas* dcaZPhiMeanTopVsBottom = helpers::prepareSimpleCanvas2DcaMeanValues(*mHistoImpParZMeanTop, kRed, "mean_{Z}^{TOP}", *mHistoImpParZMeanBottom, kBlue, "mean_{Z}^{BOTTOM}");
507
508 dcaXyVsdcaZ->SaveAs(Form("%s/ComparisonDCAXYvsZ.png", directoryName.Data()));
509 dcaXyTopVsBottom->SaveAs(Form("%s/ComparisonDCAXyTopVsBottom.png", directoryName.Data()));
510 dcaXyPosVsNeg->SaveAs(Form("%s/ComparisonDCAXyPosVsNeg.png", directoryName.Data()));
511 dcaZTopVsBottom->SaveAs(Form("%s/ComparisonDCAZTopVsBottom.png", directoryName.Data()));
512 dcaZPosVsNeg->SaveAs(Form("%s/ComparisonDCAZPosVsNeg.png", directoryName.Data()));
513 dcaXYvsPhi->SaveAs(Form("%s/dcaXYvsPhi.png", directoryName.Data()));
514 dcaZvsPhi->SaveAs(Form("%s/dcaZvsPhi.png", directoryName.Data()));
515 dcaXyPhiMeanTopVsBottom->SaveAs(Form("%s/dcaXyPhiMeanTopVsBottom.png", directoryName.Data()));
516 dcaZPhiMeanTopVsBottom->SaveAs(Form("%s/dcaZPhiMeanTopVsBottom.png", directoryName.Data()));
517}
518
519void ImpactParameterStudy::endOfStream(EndOfStreamContext& ec)
520{
521 auto& params = ITSImpactParameterParamConfig::Instance();
522 if (params.generatePlots) {
523 saveHistograms();
524 plotHistograms();
525 }
526}
527
528void ImpactParameterStudy::finaliseCCDB(ConcreteDataMatcher& matcher, void* obj)
529{
530 if (o2::base::GRPGeomHelper::instance().finaliseCCDB(matcher, obj)) {
531 return;
532 }
533}
534
535DataProcessorSpec getImpactParameterStudy(mask_t srcTracksMask, mask_t srcClustersMask, bool useMC)
536{
537 std::vector<OutputSpec> outputs;
538 auto dataRequest = std::make_shared<DataRequest>();
539 dataRequest->requestTracks(srcTracksMask, useMC);
540 dataRequest->requestPrimaryVertices(useMC);
541
542 auto ggRequest = std::make_shared<o2::base::GRPGeomRequest>(false, // orbitResetTime
543 true, // GRPECS=true
544 true, // GRPLHCIF
545 true, // GRPMagField
546 true, // askMatLUT
547 o2::base::GRPGeomRequest::Aligned, // geometry
548 dataRequest->inputs,
549 true);
550
551 return DataProcessorSpec{
552 "its-study-impactparameter",
553 dataRequest->inputs,
554 outputs,
555 AlgorithmSpec{adaptFromTask<ImpactParameterStudy>(dataRequest, ggRequest, srcTracksMask, useMC)},
556 Options{}};
557}
558
559} // namespace study
560} // namespace its
561} // namespace o2
RecoContainer.h
Wrapper container for different reconstructed object types.
TrackCuts.h
Class to perform track cuts.
i
int32_t i
Definition GPUCommonAlgorithm.h:443
GRPECSObject.h
Header of the AggregatedRunInfo struct.
GRPGeomHelper.h
Helper for geometry and GRP related CCDB requests.
PVertexer.h
Primary vertex finder.
VtxTrackRef.h
Referenc on track indices contributing to the vertex, with possibility chose tracks from specific sou...
o2::base::GRPGeomHelper::checkUpdates
void checkUpdates(o2::framework::ProcessingContext &pc)
Definition GRPGeomHelper.cxx:178
o2::base::GRPGeomHelper::instance
static GRPGeomHelper & instance()
Definition GRPGeomHelper.h:120
o2::base::GRPGeomHelper::setRequest
void setRequest(std::shared_ptr< GRPGeomRequest > req)
Definition GRPGeomHelper.cxx:101
o2::base::PropagatorImpl< float >::Instance
GPUd() value_type estimateLTFast(o2 static GPUd() float estimateLTIncrement(const o2 PropagatorImpl * Instance(bool uninitialized=false)
Definition Propagator.h:143
o2::conf::ConfigurableParam::updateFromString
static void updateFromString(std::string const &)
Definition ConfigurableParam.cxx:433
o2::detectors::DetID
Static class with identifiers, bitmasks and names for ALICE detectors.
Definition DetID.h:58
o2::framework::ProcessingContext::services
ServiceRegistryRef services()
The services registry associated with this processing context.
Definition ProcessingContext.h:39
o2::its::study::ImpactParameterStudy::run
void run(ProcessingContext &) final
Definition ImpactParameter.cxx:233
o2::its::study::ImpactParameterStudy::finaliseCCDB
void finaliseCCDB(ConcreteDataMatcher &, void *) final
Definition ImpactParameter.cxx:528
o2::its::study::ImpactParameterStudy::ImpactParameterStudy
ImpactParameterStudy(std::shared_ptr< DataRequest > dr, std::shared_ptr< o2::base::GRPGeomRequest > gr, mask_t src, bool useMC)
Definition ImpactParameter.cxx:66
o2::its::study::ImpactParameterStudy::process
void process(o2::globaltracking::RecoContainer &)
Definition ImpactParameter.cxx:242
o2::its::study::ImpactParameterStudy::endOfStream
void endOfStream(EndOfStreamContext &) final
This is invoked whenever we have an EndOfStream event.
Definition ImpactParameter.cxx:519
o2::its::study::TrackCuts::isSelected
bool isSelected(GID trackIndex, o2::globaltracking::RecoContainer &data)
Definition TrackCuts.h:107
o2::vertexing::PVertexer::prepareVertexRefit
bool prepareVertexRefit(const TR &tracks, const o2d::VertexBase &vtxSeed)
Definition PVertexer.h:359
o2::vertexing::PVertexer::refitVertex
PVertex refitVertex(const std::vector< bool > useTrack, const o2d::VertexBase &vtxSeed)
Definition PVertexer.cxx:1304
src
GLenum src
Definition glcorearb.h:1767
entry
GLuint entry
Definition glcorearb.h:5735
params
GLenum const GLfloat * params
Definition glcorearb.h:272
o2::framework
Defining PrimaryVertex explicitly as messageable.
Definition TFIDInfo.h:20
o2::framework::Options
std::vector< ConfigParamSpec > Options
Definition DataProcessorSpec.h:31
o2::gpu::gpustd::array
std::array< T, N > array
Definition GPUCommonArray.h:38
o2::its::study::helpers::makeLogBinning
std::vector< double > makeLogBinning(const int nbins, const double min, const double max)
get a vector containing binning info for constant sized bins on a log axis
Definition Helpers.cxx:21
o2::its::study::helpers::prepareSimpleCanvas2Histograms
TCanvas * prepareSimpleCanvas2Histograms(TH1F &h1, int color1, TH1F &h2, int color2)
prepare canvas with two TH1F plots
Definition Helpers.cxx:89
o2::its::study::helpers::prepareSimpleCanvas2DcaMeanValues
TCanvas * prepareSimpleCanvas2DcaMeanValues(TH1F &h1, int color1, TString nameHisto1, TH1F &h2, int color2, TString nameHisto2)
Definition Helpers.cxx:160
o2::its::study::helpers::plot2DwithMeanAndSigma
TCanvas * plot2DwithMeanAndSigma(TH2F &h2D, TH1F &hMean, TH1F &hSigma, int color)
plot canvas with TH2D + TH1D(Mean and Sigma) from slice
Definition Helpers.cxx:178
o2::its::study::getImpactParameterStudy
o2::framework::DataProcessorSpec getImpactParameterStudy(mask_t srcTracksMask, mask_t srcClusMask, bool useMC=false)
Definition ImpactParameter.cxx:535
o2::track::TrackParCov
TrackParCovF TrackParCov
Definition Track.h:33
o2
a couple of static helper functions to create timestamp values for CCDB queries or override obsolete ...
Definition BitstreamReader.h:24
std
Defining DataPointCompositeObject explicitly as copiable.
Definition CcdbObjectInfo.h:121
o2::framework::DeviceSpec::inputTimesliceId
size_t inputTimesliceId
The time pipelining id of this particular device.
Definition DeviceSpec.h:68
o2::globaltracking::RecoContainer::getSingleDetectorRefs
GlobalIDSet getSingleDetectorRefs(GTrackID gidx) const
Definition RecoContainer.cxx:1389
o2::globaltracking::RecoContainer::getTrackParam
const o2::track::TrackParCov & getTrackParam(GTrackID gidx) const
Definition RecoContainer.cxx:1347
o2::globaltracking::RecoContainer::collectData
void collectData(o2::framework::ProcessingContext &pc, const DataRequest &request)
Definition RecoContainer.cxx:609
LOG
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"