CbmRoot
Loading...
Searching...
No Matches
CbmAnaConversionRich.cxx
Go to the documentation of this file.
1/* Copyright (C) 2014-2020 Fakultaet fuer Mathematik und Naturwissenschaften, Bergische Universitaet Wuppertal, Wuppertal
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Sascha Reinecke [committer], Florian Uhlig */
4
12#include "CbmAnaConversionRich.h"
13
14#include "CbmGlobalTrack.h"
15#include "CbmMCTrack.h"
16#include "CbmRichPoint.h"
17#include "CbmRichRing.h"
18#include "CbmRichUtil.h"
19#include "CbmStsTrack.h"
20#include "CbmTrackMatchNew.h"
21
22#include "FairRootManager.h"
23#include "Logger.h"
24
25#include <algorithm>
26#include <map>
27
28
29using namespace std;
30
31
32CbmAnaConversionRich::CbmAnaConversionRich()
33 : fRichPoints(NULL)
34 , fRichRings(NULL)
35 , fRichRingMatches(NULL)
36 , fMcTracks(NULL)
37 , fStsTracks(NULL)
38 , fStsTrackMatches(NULL)
39 , fGlobalTracks(NULL)
40 , fPrimVertex(NULL)
41 , fHistoList_richrings()
42 , fTest(NULL)
43 , fRichRings_nofRings(NULL)
44 , fRichRings_motherpdg(NULL)
45 , fRichRings_richpdg(NULL)
46 , fRichRings_electronspE(NULL)
47 , fRichRings_sourcePI0(NULL)
48 , fRichRings_test(NULL)
49 , fRichRings_detectedParticles(NULL)
50 , fRichRings_detParticlesMother(NULL)
51 , fRichRings_Aaxis(NULL)
52 , fRichRings_Aaxis_part1(NULL)
53 , fRichRings_Aaxis_part2(NULL)
54 , fRichRings_Aaxis_part3(NULL)
55 , fRichRings_Aaxis_electrons(NULL)
56 , fRichRings_Baxis(NULL)
57 , fRichRings_Baxis_part1(NULL)
58 , fRichRings_Baxis_part2(NULL)
59 , fRichRings_Baxis_part3(NULL)
60 , fRichRings_Baxis_electrons(NULL)
61 , fRichRings_radius(NULL)
62 , fRichRings_radius_electrons(NULL)
63 , fRichRings_radius_vs_momentum(NULL)
64 , fRichRings_radius_vs_pt(NULL)
65 , fRichRings_distance(NULL)
66 , fRichRings_distance_electrons(NULL)
67 , fhRingtest(NULL)
68 , fhRichRings_AaxisVSmom(NULL)
69 , fhRichRings_BaxisVSmom(NULL)
70 , fhRichRings_test1(NULL)
71 , fhRichRings_test2(NULL)
72 , fhRichRings_test3(NULL)
73 , fhRichRings_test4(NULL)
74 , fhRichRings_test5(NULL)
75 , fhRichRings_test6(NULL)
76 , fhRichRings_pos(NULL)
77 , fhRichRings_protons(NULL)
78 , fhRichRings_protons2(NULL)
79 , fhRichRings_start(NULL)
80 , timer()
81 , fTime(0.)
82{
83}
84
85CbmAnaConversionRich::~CbmAnaConversionRich() {}
86
87
88void CbmAnaConversionRich::Init()
89{
90 FairRootManager* ioman = FairRootManager::Instance();
91 if (NULL == ioman) { Fatal("CbmAnaConversion::Init", "RootManager not instantised!"); }
92
93 fRichPoints = (TClonesArray*) ioman->GetObject("RichPoint");
94 if (NULL == fRichPoints) { Fatal("CbmAnaConversion::Init", "No RichPoint array!"); }
95
96 fMcTracks = (TClonesArray*) ioman->GetObject("MCTrack");
97 if (NULL == fMcTracks) { Fatal("CbmAnaConversion::Init", "No MCTrack array!"); }
98
99 fStsTracks = (TClonesArray*) ioman->GetObject("StsTrack");
100 if (NULL == fStsTracks) { Fatal("CbmAnaConversion::Init", "No StsTrack array!"); }
101
102 fStsTrackMatches = (TClonesArray*) ioman->GetObject("StsTrackMatch");
103 if (NULL == fStsTrackMatches) { Fatal("CbmAnaConversion::Init", "No StsTrackMatch array!"); }
104
105 fGlobalTracks = (TClonesArray*) ioman->GetObject("GlobalTrack");
106 if (NULL == fGlobalTracks) { Fatal("CbmAnaConversion::Init", "No GlobalTrack array!"); }
107
108 // Get pointer to PrimaryVertex object from IOManager if it exists
109 // The old name for the object is "PrimaryVertex" the new one
110 // "PrimaryVertex." Check first for the new name
111 fPrimVertex = dynamic_cast<CbmVertex*>(ioman->GetObject("PrimaryVertex."));
112 if (nullptr == fPrimVertex) { fPrimVertex = dynamic_cast<CbmVertex*>(ioman->GetObject("PrimaryVertex")); }
113 if (nullptr == fPrimVertex) { LOG(fatal) << "No PrimaryVertex array!"; }
114
115 fRichRings = (TClonesArray*) ioman->GetObject("RichRing");
116 if (NULL == fRichRings) { Fatal("CbmAnaConversion::Init", "No RichRing array!"); }
117
118 fRichRingMatches = (TClonesArray*) ioman->GetObject("RichRingMatch");
119 if (NULL == fRichRingMatches) { Fatal("CbmAnaConversion::Init", "No RichRingMatch array!"); }
120
121 InitHistos();
122}
123
124
125void CbmAnaConversionRich::InitHistos()
126{
127 fHistoList_richrings.clear();
128
129 fTest = new TH2D("fTest", "fTest; X; Y", 200, -100, 100, 400, -200, 200);
130 fHistoList_richrings.push_back(fTest);
131
132
133 // #############################################
134 // Histograms related to rich rings
135 fRichRings_nofRings = new TH1D("fRichRings_nofRings", "fRichRings_nofRings;nof rings per Event;#", 101, -0.5, 100.5);
136 fRichRings_motherpdg = new TH1D("fRichRings_motherpdg", "fRichRings_motherpdg;pdg code;#", 5002, -1.5, 5000.5);
137 fRichRings_richpdg = new TH1D("fRichRings_richpdg", "fRichRings_richpdg;pdg code;#", 5001, -0.5, 5000.5);
138 fRichRings_electronspE =
139 new TH1D("fRichRings_electronspE", "fRichRings_electronspE;nof electrons per event;#", 60, -0.5, 59.5);
140 fRichRings_sourcePI0 =
141 new TH1D("fRichRings_sourcePI0", "fRichRings_sourcePI0;nof electrons from pi0;#", 60, -0.5, 59.5);
142 fRichRings_test = new TH1D("fRichRings_test", "fRichRings_test;nof electrons from pi0;#", 20, -0.5, 19.5);
143 fHistoList_richrings.push_back(fRichRings_nofRings);
144 fHistoList_richrings.push_back(fRichRings_motherpdg);
145 fHistoList_richrings.push_back(fRichRings_richpdg);
146 fHistoList_richrings.push_back(fRichRings_electronspE);
147 fHistoList_richrings.push_back(fRichRings_sourcePI0);
148 fHistoList_richrings.push_back(fRichRings_test);
149
150 fRichRings_detectedParticles = new TH1D("fRichRings_detectedParticles", "fRichRings_detectedParticles;;#", 9, 0., 9.);
151 fRichRings_detectedParticles->SetLabelSize(0.06);
152 fHistoList_richrings.push_back(fRichRings_detectedParticles);
153 fRichRings_detectedParticles->GetXaxis()->SetBinLabel(1, "e^{+}/e^{-}"); // pdg code 11
154 fRichRings_detectedParticles->GetXaxis()->SetBinLabel(2, "\\mu^{+}/\\mu^{-}"); // pdg code 13
155 fRichRings_detectedParticles->GetXaxis()->SetBinLabel(3, "\\pi^{+}/\\pi^{-}"); // pdg code 211
156 fRichRings_detectedParticles->GetXaxis()->SetBinLabel(4, "K^{+}/K^{-}"); // pdg code 321
157 fRichRings_detectedParticles->GetXaxis()->SetBinLabel(5,
158 "p"); // pdg code 2212
159 fRichRings_detectedParticles->GetXaxis()->SetBinLabel(6, "\\Sigma^-"); // pdg code 3112
160 fRichRings_detectedParticles->GetXaxis()->SetBinLabel(7, "\\Sigma^+"); // pdg code 3222
161 fRichRings_detectedParticles->GetXaxis()->SetBinLabel(8, "\\Xi^-"); // pdg code 3312
162 fRichRings_detectedParticles->GetXaxis()->SetBinLabel(9, "else"); // everything else
163
164 fRichRings_detParticlesMother =
165 new TH1D("fRichRings_detParticlesMother", "fRichRings_detParticlesMother;;#", 22, 0., 22.);
166 fRichRings_detParticlesMother->SetLabelSize(0.04);
167 fHistoList_richrings.push_back(fRichRings_detParticlesMother);
168 fRichRings_detParticlesMother->GetXaxis()->SetBinLabel(1, "e^{+}/e^{-}"); // pdg code 11
169 fRichRings_detParticlesMother->GetXaxis()->SetBinLabel(2, "\\mu^{+}/\\mu^{-}"); // pdg code 13
170 fRichRings_detParticlesMother->GetXaxis()->SetBinLabel(3, "\\gamma"); // pdg code 22
171 fRichRings_detParticlesMother->GetXaxis()->SetBinLabel(4, "\\pi^{0}"); // pdg code 111
172 fRichRings_detParticlesMother->GetXaxis()->SetBinLabel(5, "K_{L}^{0}"); // pdg code 130
173 fRichRings_detParticlesMother->GetXaxis()->SetBinLabel(6, "\\pi^{+}/\\pi^{-}"); // pdg code 211
174 fRichRings_detParticlesMother->GetXaxis()->SetBinLabel(7, "\\eta"); // pdg code 221
175 fRichRings_detParticlesMother->GetXaxis()->SetBinLabel(8, "K_{S}^{0}"); // pdg code 310
176 fRichRings_detParticlesMother->GetXaxis()->SetBinLabel(9, "K^{+}/K^{-}"); // pdg code 321
177 fRichRings_detParticlesMother->GetXaxis()->SetBinLabel(10,
178 "n"); // pdg code 2112
179 fRichRings_detParticlesMother->GetXaxis()->SetBinLabel(11,
180 "p"); // pdg code 2212
181 fRichRings_detParticlesMother->GetXaxis()->SetBinLabel(12, "\\Sigma^-"); // pdg code 3112
182 fRichRings_detParticlesMother->GetXaxis()->SetBinLabel(13, "\\Lambda^0"); // pdg code 3122
183 fRichRings_detParticlesMother->GetXaxis()->SetBinLabel(14, "\\Sigma^+"); // pdg code 3222
184 fRichRings_detParticlesMother->GetXaxis()->SetBinLabel(15, "\\Xi^-"); // pdg code 3312
185 fRichRings_detParticlesMother->GetXaxis()->SetBinLabel(16, "\\Omega^-"); // pdg code 3334
186 fRichRings_detParticlesMother->GetXaxis()->SetBinLabel(
187 17, "prim."); // no mother particle, i.e. comes from primary vertex
188 fRichRings_detParticlesMother->GetXaxis()->SetBinLabel(
189 18,
190 "prim. E"); // no mother particle, i.e. comes from primary vertex (e+/e-)
191 fRichRings_detParticlesMother->GetXaxis()->SetBinLabel(
192 19,
193 "prim. M"); // no mother particle, i.e. comes from primary vertex (mu+/mu-)
194 fRichRings_detParticlesMother->GetXaxis()->SetBinLabel(
195 20,
196 "prim. P"); // no mother particle, i.e. comes from primary vertex (pi+/pi-)
197 fRichRings_detParticlesMother->GetXaxis()->SetBinLabel(
198 21,
199 "prim. K"); // no mother particle, i.e. comes from primary vertex (K+/K-)
200 fRichRings_detParticlesMother->GetXaxis()->SetBinLabel(
201 22,
202 "prim. else"); // no mother particle, i.e. comes from primary vertex (K+/K-)
203
204
205 fRichRings_Aaxis = new TH1D("fRichRings_Aaxis", "fRichRings_Aaxis;A-axis;#", 300, 0., 30.);
206 fRichRings_Aaxis_part1 = new TH1D("fRichRings_Aaxis_part1", "fRichRings_Aaxis_part1;A-axis;#", 300, 0., 30.);
207 fRichRings_Aaxis_part2 = new TH1D("fRichRings_Aaxis_part2", "fRichRings_Aaxis_part2;A-axis;#", 300, 0., 30.);
208 fRichRings_Aaxis_part3 = new TH1D("fRichRings_Aaxis_part3", "fRichRings_Aaxis_part3;A-axis;#", 300, 0., 30.);
209 fRichRings_Aaxis_electrons =
210 new TH1D("fRichRings_Aaxis_electrons", "fRichRings_Aaxis_electrons;A-axis;#", 300, 0., 30.);
211 fRichRings_Baxis = new TH1D("fRichRings_Baxis", "fRichRings_Baxis;B-axis;#", 300, 0., 30.);
212 fRichRings_Baxis_part1 = new TH1D("fRichRings_Baxis_part1", "fRichRings_Baxis_part1;B-axis;#", 300, 0., 30.);
213 fRichRings_Baxis_part2 = new TH1D("fRichRings_Baxis_part2", "fRichRings_Baxis_part2;B-axis;#", 300, 0., 30.);
214 fRichRings_Baxis_part3 = new TH1D("fRichRings_Baxis_part3", "fRichRings_Baxis_part3;B-axis;#", 300, 0., 30.);
215 fRichRings_Baxis_electrons =
216 new TH1D("fRichRings_Baxis_electrons", "fRichRings_Baxis_electrons;B-axis;#", 300, 0., 30.);
217 fRichRings_radius = new TH1D("fRichRings_radius", "fRichRings_radius;radius;#", 300, 0., 30.);
218 fRichRings_radius_electrons =
219 new TH1D("fRichRings_radius_electrons", "fRichRings_radius_electrons;radius;#", 300, 0., 30.);
220 fRichRings_radius_vs_momentum = new TH2D("fRichRings_radius_vs_momentum",
221 "fRichRings_radius_vs_momentum;momentum;radius", 100, 0., 10., 300, 0., 30.);
222 fRichRings_radius_vs_pt =
223 new TH2D("fRichRings_radius_vs_pt", "fRichRings_radius_vs_pt;pt;radius", 100, 0., 10., 300, 0., 30.);
224 fRichRings_distance = new TH1D("fRichRings_distance", "fRichRings_distance;distance;#", 500, 0., 5.);
225 fRichRings_distance_electrons =
226 new TH1D("fRichRings_distance_electrons", "fRichRings_distance_electrons;distance;#", 500, 0., 5.);
227 fhRingtest = new TH2D("fhRingtest", "fhRingtest;momentum [GeV/c];radius [cm]", 200, 0., 10., 200, 0., 10.);
228 fhRichRings_AaxisVSmom = new TH2D("fhRichRings_AaxisVSmom", "fhRichRings_AaxisVSmom;momentum [GeV/c];a-axis [cm]",
229 200, 0., 10., 200, 0., 10.);
230 fhRichRings_BaxisVSmom = new TH2D("fhRichRings_BaxisVSmom", "fhRichRings_BaxisVSmom;momentum [GeV/c];b-axis [cm]",
231 200, 0., 10., 200, 0., 10.);
232 fHistoList_richrings.push_back(fRichRings_Aaxis);
233 fHistoList_richrings.push_back(fRichRings_Aaxis_part1);
234 fHistoList_richrings.push_back(fRichRings_Aaxis_part2);
235 fHistoList_richrings.push_back(fRichRings_Aaxis_part3);
236 fHistoList_richrings.push_back(fRichRings_Aaxis_electrons);
237 fHistoList_richrings.push_back(fRichRings_Baxis);
238 fHistoList_richrings.push_back(fRichRings_Baxis_part1);
239 fHistoList_richrings.push_back(fRichRings_Baxis_part2);
240 fHistoList_richrings.push_back(fRichRings_Baxis_part3);
241 fHistoList_richrings.push_back(fRichRings_Baxis_electrons);
242 fHistoList_richrings.push_back(fRichRings_radius);
243 fHistoList_richrings.push_back(fRichRings_radius_electrons);
244 fHistoList_richrings.push_back(fRichRings_radius_vs_momentum);
245 fHistoList_richrings.push_back(fRichRings_radius_vs_pt);
246 fHistoList_richrings.push_back(fRichRings_distance);
247 fHistoList_richrings.push_back(fRichRings_distance_electrons);
248 fHistoList_richrings.push_back(fhRingtest);
249 fHistoList_richrings.push_back(fhRichRings_AaxisVSmom);
250 fHistoList_richrings.push_back(fhRichRings_BaxisVSmom);
251
252 fhRichRings_test1 = new TH1D("fRichRings_test1", "fRichRings_test1;A-axis;#", 300, 0., 30.);
253 fhRichRings_test2 = new TH1D("fRichRings_test2", "fRichRings_test2;A-axis;#", 300, 0., 30.);
254 fhRichRings_test3 = new TH1D("fRichRings_test3", "fRichRings_test3;A-axis;#", 300, 0., 5.);
255 fhRichRings_test4 = new TH1D("fRichRings_test4", "fRichRings_test4;A-axis;#", 300, 0., 30.);
256 fhRichRings_test5 = new TH1D("fRichRings_test5", "fRichRings_test5;A-axis;#", 300, 0., 30.);
257 fhRichRings_test6 = new TH1D("fRichRings_test6", "fRichRings_test6;A-axis;#", 300, 0., 30.);
258 fhRichRings_pos = new TH2D("fRichRings_pos", "fRichRings_pos;x;y", 400, -200., 200., 400, -200., 200.);
259 fHistoList_richrings.push_back(fhRichRings_test1);
260 fHistoList_richrings.push_back(fhRichRings_test2);
261 fHistoList_richrings.push_back(fhRichRings_test3);
262 fHistoList_richrings.push_back(fhRichRings_test4);
263 fHistoList_richrings.push_back(fhRichRings_test5);
264 fHistoList_richrings.push_back(fhRichRings_test6);
265 fHistoList_richrings.push_back(fhRichRings_pos);
266
267
268 // for proton analysis
269 fhRichRings_protons = new TH2D("fhRichRings_protons", "fhRichRings_protons;px;py", 1000, -10, 10, 1000, -10, 10);
270 fHistoList_richrings.push_back(fhRichRings_protons);
271 fhRichRings_protons2 = new TH1D("fhRichRings_protons2", "fhRichRings_protons2;nof hits;#", 40, -0.5, 39.5);
272 fHistoList_richrings.push_back(fhRichRings_protons2);
273
274
275 fhRichRings_start = new TH1D("fhRichRings_start", "fhRichRings_start;z;#", 1000, 0, 100);
276 fHistoList_richrings.push_back(fhRichRings_start);
277}
278
279
280void CbmAnaConversionRich::Finish()
281{
282 //gDirectory->cd("analysis-conversion");
283 gDirectory->mkdir("RichRings");
284 gDirectory->cd("RichRings");
285 for (UInt_t i = 0; i < fHistoList_richrings.size(); i++) {
286 fHistoList_richrings[i]->Write();
287 }
288 gDirectory->cd("..");
289
290 cout << "CbmAnaConversionRich: Realtime - " << fTime << endl;
291 //timer.Print();
292}
293
294
295void CbmAnaConversionRich::AnalyseRICHdata()
296{
297 timer.Start();
298
299 Int_t nofMC = fMcTracks->GetEntriesFast();
300 cout << "CbmAnaConversionRich: nof mc tracks " << nofMC << endl;
301
302 // analyse RICH points
303 Int_t nofPoints = fRichPoints->GetEntriesFast();
304 for (int i = 0; i < nofPoints; i++) {
305 CbmRichPoint* point = (CbmRichPoint*) fRichPoints->At(i);
306 if (point == NULL) continue;
307 //cout << point->GetX() << "\t" << point->GetY() << endl;
308 fTest->Fill(point->GetX(), point->GetY());
309 }
310
311 // analyse RICH rings
312 if (fRichRings != NULL) {
313 // Bool_t testbool = true;
314 Int_t nofElectrons = 0;
315 Int_t sourcePI0 = 0;
316 vector<int> pi0ids;
317 pi0ids.clear();
318
319 Int_t primEl = 0;
320 Int_t primMu = 0;
321 Int_t primPi = 0;
322 Int_t primK = 0;
323
324 Int_t nofRings = fRichRings->GetEntriesFast();
325 fRichRings_nofRings->Fill(nofRings);
326 for (int i = 0; i < nofRings; i++) {
327 CbmRichRing* richRing = (CbmRichRing*) fRichRings->At(i);
328 if (richRing == NULL) continue;
329 fRichRings_Aaxis->Fill(richRing->GetAaxis());
330 fRichRings_Baxis->Fill(richRing->GetBaxis());
331 // CbmRichRing::GetDistance() method is no longer supported
332 // If you wan to use cuts update code using CbmRichUtil::GetRingTrackDistance()
333 fRichRings_distance->Fill(1. /*richRing->GetDistance()*/);
334 fRichRings_radius->Fill(richRing->GetRadius());
335 }
336
337 Int_t nofGlobalTracks = fGlobalTracks->GetEntriesFast();
338 cout << "CbmAnaConversionRich: nof global tracks " << nofGlobalTracks << endl;
339 for (int iG = 0; iG < nofGlobalTracks; iG++) {
340 CbmGlobalTrack* gTrack = (CbmGlobalTrack*) fGlobalTracks->At(iG);
341 if (NULL == gTrack) continue;
342 int stsInd = gTrack->GetStsTrackIndex();
343 int richInd = gTrack->GetRichRingIndex();
344 if (richInd < 0) continue; // no RICH segment -> no ring
345
346 //int trdInd = gTrack->GetTrdTrackIndex();
347 //int tofInd = gTrack->GetTofHitIndex();
348
349 if (stsInd < 0) continue;
350 CbmStsTrack* stsTrack = (CbmStsTrack*) fStsTracks->At(stsInd);
351 if (stsTrack == NULL) continue;
352
353 CbmTrackMatchNew* stsMatch = (CbmTrackMatchNew*) fStsTrackMatches->At(stsInd);
354 if (stsMatch == NULL) continue;
355 int stsMcTrackId = stsMatch->GetMatchedLink().GetIndex();
356 if (stsMcTrackId < 0) continue;
357 CbmMCTrack* mcTrack1 = (CbmMCTrack*) fMcTracks->At(stsMcTrackId);
358 if (mcTrack1 == NULL) continue;
359
360 CbmTrackMatchNew* richMatch = (CbmTrackMatchNew*) fRichRingMatches->At(richInd);
361 if (richMatch == NULL) continue;
362 int richMcTrackId = richMatch->GetMatchedLink().GetIndex();
363 if (richMcTrackId < 0) continue;
364 CbmMCTrack* mcTrack2 = (CbmMCTrack*) fMcTracks->At(richMcTrackId);
365 if (mcTrack2 == NULL) continue;
366
367 // stsMcTrackId == richMcTrackId -> track was reconstructed in STS and made a ring in RICH, track matching was correct
368 // in case they are not equal, the ring comes either from a secondary particle or STS track was not reconstructed
369 if (stsMcTrackId != richMcTrackId) continue;
370
371 CbmRichRing* ring = static_cast<CbmRichRing*>(fRichRings->At(richInd));
372 TVector3 momentum;
373 mcTrack1->GetMomentum(momentum);
374 fRichRings_radius_vs_momentum->Fill(momentum.Mag(), ring->GetRadius());
375 fRichRings_radius_vs_pt->Fill(momentum.Perp(), ring->GetRadius());
376
377
378 int pdg = TMath::Abs(mcTrack2->GetPdgCode()); // extract pdg code of particle directly from rich ring
379
380 // get start vertex of track
381 TVector3 start;
382 mcTrack2->GetStartVertex(start);
383 fhRichRings_start->Fill(start.Z());
384
385
386 /* Int_t index = richRing->GetTrackID();
387 if (index < 0) continue;
388 CbmGlobalTrack* gTrack = (CbmGlobalTrack*) fGlobalTracks->At(index);
389 if(NULL == gTrack) continue;
390 int stsInd = gTrack->GetStsTrackIndex();
391 int richInd = gTrack->GetRichRingIndex();
392 int trdInd = gTrack->GetTrdTrackIndex();
393 int tofInd = gTrack->GetTofHitIndex();
394
395 if (stsInd < 0) continue;
396 CbmStsTrack* stsTrack = (CbmStsTrack*) fStsTracks->At(stsInd);
397 if (stsTrack == NULL) continue;
398 CbmTrackMatchNew* stsMatch = (CbmTrackMatchNew*) fStsTrackMatches->At(stsInd);
399 if (stsMatch == NULL) continue;
400 int stsMcTrackId = stsMatch->GetMatchedLink().GetIndex();
401 if (stsMcTrackId < 0) continue;
402 CbmMCTrack* mcTrack1 = (CbmMCTrack*) fMcTracks->At(stsMcTrackId);
403 if (mcTrack1 == NULL) continue;
404 int pdg = TMath::Abs(mcTrack1->GetPdgCode());
405*/
406 int motherId = mcTrack2->GetMotherId();
407 int motherpdg = 0;
408 CbmMCTrack* mothermcTrack1;
409
410 if (motherId != -1) {
411 mothermcTrack1 = (CbmMCTrack*) fMcTracks->At(motherId);
412 motherpdg = TMath::Abs(mothermcTrack1->GetPdgCode());
413 }
414 if (motherId == -1) {
415 motherpdg = -1;
416 if (pdg == 11) primEl++;
417 if (pdg == 13) primMu++;
418 if (pdg == 211) primPi++;
419 if (pdg == 321) primK++;
420 }
421
422 if (pdg == 2212) {
423 Protons(mcTrack1);
424 //fhRichRings_protons2->Fill(richRing->GetNofHits());
425 }
426
427 fRichRings_motherpdg->Fill(motherpdg);
428 fRichRings_richpdg->Fill(pdg);
429
430 if (pdg == 11) { // && stsMcTrackId == richMcTrackId) {
431 nofElectrons++;
432 int grandmotherpdg = 0;
433 int grandmotherId = 0;
434 if (motherpdg == 22) {
435 grandmotherId = mothermcTrack1->GetMotherId();
436 if (grandmotherId != -1) {
437 CbmMCTrack* grandmothermcTrack1 = (CbmMCTrack*) fMcTracks->At(grandmotherId);
438 grandmotherpdg = TMath::Abs(grandmothermcTrack1->GetPdgCode());
439 }
440 if (grandmotherId == -1) { grandmotherpdg = -1; }
441 }
442 if (motherpdg == 111 || grandmotherpdg == 111) {
443 sourcePI0++;
444 if (motherpdg == 111) { pi0ids.push_back(motherId); }
445 else if (grandmotherpdg == 111) {
446 pi0ids.push_back(grandmotherId);
447 }
448 }
449 }
450 FillAdditionalPDGhisto(pdg, motherpdg);
451 }
452
453 fRichRings_electronspE->Fill(nofElectrons);
454 fRichRings_sourcePI0->Fill(sourcePI0);
455
456 TH1I* zwischenhisto = new TH1I("zwischenhisto", "zwischenhisto", 1000000, 0, 1000000);
457 for (unsigned int i = 0; i < pi0ids.size(); i++) {
458 zwischenhisto->Fill(pi0ids[i]);
459 }
460 fRichRings_test->Fill(zwischenhisto->GetMaximum());
461 if (zwischenhisto->GetMaximum() >= 4) {
462 CbmMCTrack* mcTrack2 = (CbmMCTrack*) fMcTracks->At(zwischenhisto->GetMaximumBin() - 1);
463 int pdg = mcTrack2->GetPdgCode();
464 cout << "CbmAnaConversionRich: MAXIMUM BIN: " << zwischenhisto->GetMaximum() << "\t"
465 << "bin id: " << zwischenhisto->GetMaximumBin() - 1 << endl;
466 cout << "CbmAnaConversionRich: pdg code: " << pdg << endl;
467 cout << "CbmAnaConversionRich: \t";
468
469 std::sort(pi0ids.begin(), pi0ids.end());
470 for (unsigned int i = 0; i < pi0ids.size(); i++) {
471 cout << pi0ids[i] << "\t";
472 }
473
474 cout << endl;
475 cout << "CbmAnaConversionRich: number of electrons: " << nofElectrons << "\t e from pi0: " << sourcePI0
476 << "\t pi0ids size: " << pi0ids.size() << endl;
477 }
478 zwischenhisto->Delete();
479
480
481 int photoncounter = 0;
482 std::multimap<int, int> electronMap;
483 for (unsigned int i = 0; i < pi0ids.size(); i++) {
484 electronMap.insert(std::pair<int, int>(pi0ids[i], i));
485 }
486
487 int check = 0;
488 for (std::map<int, int>::iterator it = electronMap.begin(); it != electronMap.end(); ++it) {
489 if (it == electronMap.begin()) check = 1;
490 if (it != electronMap.begin()) {
491 std::map<int, int>::iterator zwischen = it;
492 zwischen--;
493 int id = it->first;
494 int id_old = zwischen->first;
495 if (id == id_old) {
496 check++;
497 if (check > 3) {
498 cout << "CbmAnaConversionRich: richmap - photon found " << id << endl;
499 photoncounter++;
500 }
501 }
502 else
503 check = 1;
504 }
505 }
506 }
507
508
509 // plot momentum vs ring radius/Axis+Baxis
510 //cout << "CbmAnaConversionRich: further ring analyses..." << endl;
511 Int_t nGTracks = fGlobalTracks->GetEntriesFast();
512 for (Int_t iGTrack = 0; iGTrack < nGTracks; iGTrack++) {
513 CbmGlobalTrack* gTrack = (CbmGlobalTrack*) fGlobalTracks->At(iGTrack);
514 if (NULL == gTrack) continue;
515 int stsInd = gTrack->GetStsTrackIndex();
516 int richInd = gTrack->GetRichRingIndex();
517 if (richInd < 0) continue;
518 if (stsInd < 0) continue;
519 CbmRichRing* richRing = (CbmRichRing*) fRichRings->At(richInd);
520
521 CbmTrackMatchNew* richMatch = (CbmTrackMatchNew*) fRichRingMatches->At(richInd);
522 if (richMatch == NULL) continue;
523 Int_t richMcTrackId = richMatch->GetMatchedLink().GetIndex();
524
525 CbmMCTrack* mctrack = (CbmMCTrack*) fMcTracks->At(richMcTrackId);
526 if (mctrack == NULL) continue;
527
528 int pdg = TMath::Abs(mctrack->GetPdgCode());
529
530 Double_t momentum = mctrack->GetP();
531 Double_t radius = richRing->GetRadius();
532 Double_t axisA = richRing->GetAaxis();
533 Double_t axisB = richRing->GetBaxis();
534 Double_t distance = CbmRichUtil::GetRingTrackDistance(iGTrack);
535
536 fhRingtest->Fill(momentum, radius);
537 fhRichRings_AaxisVSmom->Fill(momentum, axisA);
538 fhRichRings_BaxisVSmom->Fill(momentum, axisB);
539
540 if (richRing->GetPhi() <= 5) { fhRichRings_test1->Fill(axisA); }
541
542 if (pdg == 11) {
543 fRichRings_Aaxis_electrons->Fill(axisA);
544 fRichRings_Baxis_electrons->Fill(axisB);
545 fRichRings_radius_electrons->Fill(radius);
546 fRichRings_distance_electrons->Fill(distance);
547 }
548
549 Double_t ringX = richRing->GetCenterX();
550 Double_t ringY = richRing->GetCenterY();
551
552 if (sqrt(ringX * ringX + (abs(ringY) - 105) * (abs(ringY) - 105)) <= 50) { fhRichRings_test2->Fill(axisA); }
553 fhRichRings_test3->Fill(richRing->GetAngle());
554 fhRichRings_pos->Fill(ringX, ringY);
555
556 if (abs(ringX) <= 30 && abs(ringY) <= 130) {
557 fRichRings_Aaxis_part1->Fill(axisA);
558 fRichRings_Baxis_part1->Fill(axisB);
559 }
560 if (abs(ringX) > 30 && abs(ringY) > 130 && abs(ringX) <= 70 && abs(ringY) <= 155) {
561 fRichRings_Aaxis_part2->Fill(axisA);
562 fRichRings_Baxis_part2->Fill(axisB);
563 }
564 if (abs(ringX) > 70 && abs(ringY) > 155) {
565 fRichRings_Aaxis_part3->Fill(axisA);
566 fRichRings_Baxis_part3->Fill(axisB);
567 }
568 }
569
570 CheckMC();
571
572 timer.Stop();
573 fTime += timer.RealTime();
574}
575
576
577void CbmAnaConversionRich::FillAdditionalPDGhisto(Int_t pdg, Int_t motherpdg)
578{
579 if (pdg == 11) fRichRings_detectedParticles->Fill(0);
580 if (pdg == 13) fRichRings_detectedParticles->Fill(1);
581 if (pdg == 211) fRichRings_detectedParticles->Fill(2);
582 if (pdg == 321) fRichRings_detectedParticles->Fill(3);
583 if (pdg == 2212) fRichRings_detectedParticles->Fill(4);
584 if (pdg == 3112) fRichRings_detectedParticles->Fill(5);
585 if (pdg == 3222) fRichRings_detectedParticles->Fill(6);
586 if (pdg == 3312) fRichRings_detectedParticles->Fill(7);
587 if (pdg != 11 && pdg != 13 && pdg != 211 && pdg != 321 && pdg != 2212 && pdg != 3112 && pdg != 3222 && pdg != 3312)
588 fRichRings_detectedParticles->Fill(8);
589
590 if (motherpdg == 11) fRichRings_detParticlesMother->Fill(0);
591 if (motherpdg == 13) fRichRings_detParticlesMother->Fill(1);
592 if (motherpdg == 22) fRichRings_detParticlesMother->Fill(2);
593 if (motherpdg == 111) fRichRings_detParticlesMother->Fill(3);
594 if (motherpdg == 130) fRichRings_detParticlesMother->Fill(4);
595 if (motherpdg == 211) fRichRings_detParticlesMother->Fill(5);
596 if (motherpdg == 221) fRichRings_detParticlesMother->Fill(6);
597 if (motherpdg == 310) fRichRings_detParticlesMother->Fill(7);
598 if (motherpdg == 321) fRichRings_detParticlesMother->Fill(8);
599 if (motherpdg == 2112) fRichRings_detParticlesMother->Fill(9);
600 if (motherpdg == 2212) fRichRings_detParticlesMother->Fill(10);
601 if (motherpdg == 3112) fRichRings_detParticlesMother->Fill(11);
602 if (motherpdg == 3122) fRichRings_detParticlesMother->Fill(12);
603 if (motherpdg == 3222) fRichRings_detParticlesMother->Fill(13);
604 if (motherpdg == 3312) fRichRings_detParticlesMother->Fill(14);
605 if (motherpdg == 3334) fRichRings_detParticlesMother->Fill(15);
606 if (motherpdg == -1) {
607 fRichRings_detParticlesMother->Fill(16);
608 if (pdg == 11) fRichRings_detParticlesMother->Fill(17);
609 if (pdg == 13) fRichRings_detParticlesMother->Fill(18);
610 if (pdg == 211) fRichRings_detParticlesMother->Fill(19);
611 if (pdg == 321) fRichRings_detParticlesMother->Fill(20);
612 //if(pdg != 11 && pdg != 13 && pdg != 211 && pdg != 321 && pdg != 2112 && pdg != 2212 && pdg != 3112 && pdg != 3122 && pdg != 3222 && pdg != 3312 && pdg != 3334) fRichRings_detParticlesMother->Fill(21);
613 if (pdg != 11 && pdg != 13 && pdg != 211 && pdg != 321) fRichRings_detParticlesMother->Fill(21);
614 }
615}
616
617
618void CbmAnaConversionRich::Protons(CbmMCTrack* mcTrack)
619{
620 fhRichRings_protons->Fill(mcTrack->GetPx(), mcTrack->GetPy());
621}
622
623
624void CbmAnaConversionRich::CheckMC()
625{
626 vector<int> ids;
627 vector<int> electronids;
628 ids.clear();
629 electronids.clear();
630
631 Int_t nofGlobalTracks = fGlobalTracks->GetEntriesFast();
632 for (int iG = 0; iG < nofGlobalTracks; iG++) {
633 CbmGlobalTrack* gTrack = (CbmGlobalTrack*) fGlobalTracks->At(iG);
634 if (NULL == gTrack) continue;
635 int stsInd = gTrack->GetStsTrackIndex();
636 // int richInd = gTrack->GetRichRingIndex();
637 // if (richInd < 0) continue; // no RICH segment -> no ring
638 if (stsInd < 0) continue;
639
640 CbmStsTrack* stsTrack = (CbmStsTrack*) fStsTracks->At(stsInd);
641 if (stsTrack == NULL) continue;
642 CbmTrackMatchNew* stsMatch = (CbmTrackMatchNew*) fStsTrackMatches->At(stsInd);
643 if (stsMatch == NULL) continue;
644 if (stsMatch->GetNofLinks() <= 0) continue;
645 int stsMcTrackId = stsMatch->GetMatchedLink().GetIndex();
646 if (stsMcTrackId < 0) continue;
647 CbmMCTrack* mcTrack1 = (CbmMCTrack*) fMcTracks->At(stsMcTrackId);
648 if (mcTrack1 == NULL) continue;
649
650 /* CbmTrackMatchNew* richMatch = (CbmTrackMatchNew*)fRichRingMatches->At(richInd);
651 if (richMatch == NULL) continue;
652 int richMcTrackId = richMatch->GetMatchedLink().GetIndex();
653 if (richMcTrackId < 0) continue;
654 CbmMCTrack* mcTrack2 = (CbmMCTrack*) fMcTracks->At(richMcTrackId);
655 if (mcTrack2 == NULL) continue;
656 // stsMcTrackId == richMcTrackId -> track was reconstructed in STS and made a ring in RICH, track matching was correct
657 // in case they are not equal, the ring comes either from a secondary particle or STS track was not reconstructed
658 // if(stsMcTrackId != richMcTrackId) continue;
659*/
660 int pdg = TMath::Abs(mcTrack1->GetPdgCode()); // extract pdg code of particle directly from rich ring
661 // int pdg2 = TMath::Abs(mcTrack2->GetPdgCode());
662
663 // if(pdg != pdg2) {
664 // cout << "PDGs not matching!" << endl;
665 // }
666
667 int motherId = mcTrack1->GetMotherId();
668 // int motherId2 = mcTrack2->GetMotherId();
669
670 // if(motherId != motherId2) {
671 // cout << "MotherIDs not matching!" << endl;
672 // }
673
674 int motherpdg = 0;
675 int grandmotherpdg = 0;
676 CbmMCTrack* mothermcTrack1;
677 CbmMCTrack* grandmothermcTrack1;
678 int grandmotherId = 0;
679
680 if (motherId != -1) {
681 mothermcTrack1 = (CbmMCTrack*) fMcTracks->At(motherId);
682 motherpdg = TMath::Abs(mothermcTrack1->GetPdgCode());
683 }
684 else if (motherId == -1) {
685 motherpdg = -1;
686 }
687
688 if (pdg == 11) {
689 if (motherpdg == 22) {
690 grandmotherId = mothermcTrack1->GetMotherId();
691 if (grandmotherId != -1) {
692 grandmothermcTrack1 = (CbmMCTrack*) fMcTracks->At(grandmotherId);
693 grandmotherpdg = TMath::Abs(grandmothermcTrack1->GetPdgCode());
694 }
695 if (grandmotherId == -1) { grandmotherpdg = -1; }
696 }
697
698 if (motherpdg == 111) {
699 ids.push_back(motherId);
700 electronids.push_back(stsMcTrackId);
701 }
702 if (motherpdg == 22 && grandmotherpdg == 111) {
703 ids.push_back(grandmotherId);
704 electronids.push_back(stsMcTrackId);
705 }
706 }
707 }
708
709 cout << "CbmAnaConversionRich: vector ids contains " << ids.size() << " entries." << endl;
710 cout << "CbmAnaConversionRich: ids entries: ";
711 for (unsigned int i = 0; i < ids.size(); i++) {
712 cout << ids[i] << " / ";
713 }
714 cout << endl;
715 cout << "CbmAnaConversionRich: electronids entries: ";
716 for (unsigned int i = 0; i < electronids.size(); i++) {
717 cout << electronids[i] << " / ";
718 }
719 cout << endl;
720
721 int photoncounter = 0;
722 std::multimap<int, int> electronMap;
723 for (unsigned int i = 0; i < ids.size(); i++) {
724 electronMap.insert(std::pair<int, int>(ids[i], i));
725 }
726
727 int check = 0;
728 for (std::map<int, int>::iterator it = electronMap.begin(); it != electronMap.end(); ++it) {
729 if (it == electronMap.begin()) check = 1;
730 if (it != electronMap.begin()) {
731 std::map<int, int>::iterator zwischen = it;
732 zwischen--;
733 int id = it->first;
734 int id_old = zwischen->first;
735 if (id == id_old) {
736 check++;
737 if (check > 3) {
738 cout << "CbmAnaConversionRich: richmap - photon found " << id << endl;
739 photoncounter++;
740 }
741 }
742 else
743 check = 1;
744 }
745 }
746}
CbmStsTrack.h
Data class for STS tracks.
sqrt
friend fvec sqrt(const fvec &a)
Definition KfSimdPseudo.h:118
CbmAnaConversionRich::fRichRingMatches
TClonesArray * fRichRingMatches
Definition CbmAnaConversionRich.h:49
CbmAnaConversionRich::FillAdditionalPDGhisto
void FillAdditionalPDGhisto(Int_t pdg, Int_t motherpdg)
Definition CbmAnaConversionRich.cxx:577
CbmAnaConversionRich::Protons
void Protons(CbmMCTrack *mcTrack)
Definition CbmAnaConversionRich.cxx:618
CbmAnaConversionRich::fHistoList_richrings
std::vector< TH1 * > fHistoList_richrings
Definition CbmAnaConversionRich.h:56
CbmAnaConversionRich::fStsTrackMatches
TClonesArray * fStsTrackMatches
Definition CbmAnaConversionRich.h:52
CbmGlobalTrack::GetStsTrackIndex
int32_t GetStsTrackIndex() const
Definition CbmGlobalTrack.h:44
CbmGlobalTrack::GetRichRingIndex
int32_t GetRichRingIndex() const
Definition CbmGlobalTrack.h:47
CbmMCTrack::GetP
double GetP() const
Definition CbmMCTrack.h:98
CbmMCTrack::GetPx
double GetPx() const
Definition CbmMCTrack.h:70
CbmMCTrack::GetMomentum
void GetMomentum(TVector3 &momentum) const
Definition CbmMCTrack.h:170
CbmMCTrack::GetMotherId
int32_t GetMotherId() const
Definition CbmMCTrack.h:69
CbmMCTrack::GetPdgCode
int32_t GetPdgCode() const
Definition CbmMCTrack.h:68
CbmMCTrack::GetPy
double GetPy() const
Definition CbmMCTrack.h:71
CbmMCTrack::GetStartVertex
void GetStartVertex(TVector3 &vertex) const
Definition CbmMCTrack.h:176
CbmMatch::GetNofLinks
int32_t GetNofLinks() const
Definition CbmMatch.h:42
CbmMatch::GetMatchedLink
const CbmLink & GetMatchedLink() const
Definition CbmMatch.h:41
CbmRichRing::GetRadius
float GetRadius() const
Definition CbmRichRing.h:79
CbmRichRing::GetAngle
double GetAngle() const
Definition CbmRichRing.h:95
CbmRichRing::GetAaxis
double GetAaxis() const
Definition CbmRichRing.h:80
CbmRichRing::GetCenterX
float GetCenterX() const
Definition CbmRichRing.h:77
CbmRichRing::GetPhi
double GetPhi() const
Definition CbmRichRing.h:84
CbmRichRing::GetCenterY
float GetCenterY() const
Definition CbmRichRing.h:78
CbmRichRing::GetBaxis
double GetBaxis() const
Definition CbmRichRing.h:81
CbmRichUtil::GetRingTrackDistance
static double GetRingTrackDistance(int globalTrackId)
Definition CbmRichUtil.cxx:24
std
Hash for CbmL1LinkKey.
Definition algo/base/Options.cxx:21