CbmRoot
Loading...
Searching...
No Matches
CbmAnaConversionTomography.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
17
18#include "CbmMCTrack.h"
19
20#include "FairRootManager.h"
21
22#include <iostream>
23#include <map>
24
25
26using namespace std;
27
28
30 : fMcTracks(NULL)
31 , fHistoList_tomography()
32 , fhGammaZ(NULL)
33 , fhTomography(NULL)
34 , fhTomography_XZ(NULL)
35 , fhTomography_YZ(NULL)
36 , fhTomography_uptoRICH(NULL)
37 , fhTomography_RICH_complete(NULL)
38 , fhTomography_RICH_beampipe(NULL)
39 , fhTomography_STS_end(NULL)
40 , fhTomography_STS_lastStation(NULL)
41 , fhTomography_RICH_frontplate(NULL)
42 , fhTomography_RICH_backplate(NULL)
43 , fhConversion(NULL)
44 , fhConversion_cut(NULL)
45 , fhConversion_inSTS(NULL)
46 , fhConversion_prob(NULL)
47 , fhConversion_energy(NULL)
48 , fhConversion_p(NULL)
49 , fhConversion_vs_momentum(NULL)
50 , fhTomography_reco(NULL)
51 , fhTomography_reco_XZ(NULL)
52 , fhTomography_reco_YZ(NULL)
53 , fhConversion_reco(NULL)
54 , electronIDs()
55 , electronMotherIDs()
56 , conversionsInDetector()
57 ,
58 //conversionsInDetector_cut(),
59 fhConversionsPerDetector(NULL)
60 , fhConversionsPerDetectorPE(NULL)
61 ,
62 //fhConversionsPerDetectorPE_cut(NULL),
63 timer()
64 , fTime(0.)
65{
66}
67
69
70
72{
73 FairRootManager* ioman = FairRootManager::Instance();
74 if (NULL == ioman) { Fatal("CbmAnaConversionTomography::Init", "RootManager not instantised!"); }
75
76 fMcTracks = (TClonesArray*) ioman->GetObject("MCTrack");
77 if (NULL == fMcTracks) { Fatal("CbmAnaConversion::Init", "No MCTrack array!"); }
78
79
80 InitHistos();
81}
82
83
85{
87
88 fhGammaZ = new TH1D("fhGammaZ", "fhGammaZ;Z in cm;Entries", 300, 0.0, 300.);
90
91 // for all events including gammas and gamma-conversion
92 fhTomography = new TH3D("fhTomography", "Tomography/fhTomography;X in cm;Y in cm;Z in cm", 200, -200., 200., 200,
93 -200., 200., 500, 0., 1000.);
94 fhTomography_XZ = new TH2D("fhTomography_XZ", "fhTomography_XZ;X in cm;Z in cm", 3200, -800., 800., 2200, 0., 1100.);
95 fhTomography_YZ = new TH2D("fhTomography_YZ", "fhTomography_YZ;Y in cm;Z in cm", 2000, -500., 500., 2200, 0., 1100.);
96 //fhTomography_XZ_cut = new TH2D("fhTomography_XZ_cut", "fhTomography_XZ_cut;X [cm];Z [cm]", 2800, -700., 700., 2200, 0., 1100.);
97 //fhTomography_YZ_cut = new TH2D("fhTomography_YZ_cut", "fhTomography_YZ_cut;Y [cm];Z [cm]", 2000, -500., 500., 2200, 0., 1100.);
101 //fHistoList_tomography.push_back(fhTomography_XZ_cut);
102 //fHistoList_tomography.push_back(fhTomography_YZ_cut);
103
105 new TH2D("fhTomography_uptoRICH", "fhTomography_uptoRICH;X in cm;Y in cm", 600, -300., 300., 400, -200., 200.);
106 fhTomography_RICH_complete = new TH2D("fhTomography_RICH_complete", "fhTomography_RICH_complete;X in cm;Y in cm", 600,
107 -300., 300., 500, -250., 250.);
108 fhTomography_RICH_beampipe = new TH2D("fhTomography_RICH_beampipe", "fhTomography_RICH_beampipe;X in cm;Y in cm", 400,
109 -100., 100., 200, -50., 50.);
111 new TH2D("fhTomography_STS_end", "fhTomography_STS_end;X in cm;Y in cm", 400, -200., 200., 400, -200., 200.);
113 "fhTomography_STS_lastStation", "fhTomography_STS_lastStation;X in cm;Y in cm", 200, -100., 100., 200, -100., 100.);
115 "fhTomography_RICH_frontplate", "fhTomography_RICH_frontplate;X in cm;Y in cm", 600, -300., 300., 500, -250., 250.);
116 fhTomography_RICH_backplate = new TH2D("fhTomography_RICH_backplate", "fhTomography_RICH_backplate;X in cm;Y in cm",
117 600, -300., 300., 500, -250., 250.);
125
126
127 fhConversion = new TH1D("fhConversion", "fhConversion;Z in cm;# conversions", 4800, -0.5, 1199.5);
128 fhConversion_cut = new TH1D("fhConversion_cut", "fhConversion_cut;Z in cm;# conversions", 4800, -0.5, 1199.5);
129 fhConversion_inSTS = new TH1D("fhConversion_inSTS", "fhConversion in STS;Z in cm;# conversions", 800, -0.5, 199.5);
130 fhConversion_prob = new TH1D("fhConversion_prob", "fhConversion_prob;Z in cm;# conversions", 1200, -0.5, 1199.5);
135
136 fhConversion_energy = new TH1D("fhConversion_energy", "fhConversion_energy;energy;#", 1000, 0., 100.);
137 fhConversion_p = new TH1D("fhConversion_p", "fhConversion_p;p;#", 1000, 0., 100.);
140
141 fhConversion_vs_momentum = new TH2D("fhConversion_vs_momentum", "fhConversion_vs_momentum;Z in cm;momentum in GeV/c",
142 4800, -0.5, 1199.5, 1000, 0., 10.);
144
145 fhConversionsPerDetector = new TH1I("fhConversionsPerDetector", "fhConversionsPerDetector;;#", 5, 0, 5);
147 fhConversionsPerDetector->GetXaxis()->SetBinLabel(1, "magnet");
148 fhConversionsPerDetector->GetXaxis()->SetBinLabel(2, "STS");
149 fhConversionsPerDetector->GetXaxis()->SetBinLabel(3, "RICH");
150 fhConversionsPerDetector->GetXaxis()->SetBinLabel(4, "TRD");
151 fhConversionsPerDetector->GetXaxis()->SetBinLabel(5, "TOF");
152
154 new TH2I("fhConversionsPerDetectorPE", "fhConversionsPerDetectorPE;;conversions per event", 5, 0, 5, 1500, 0, 1500);
156 fhConversionsPerDetectorPE->GetXaxis()->SetBinLabel(1, "magnet");
157 fhConversionsPerDetectorPE->GetXaxis()->SetBinLabel(2, "STS");
158 fhConversionsPerDetectorPE->GetXaxis()->SetBinLabel(3, "RICH");
159 fhConversionsPerDetectorPE->GetXaxis()->SetBinLabel(4, "TRD");
160 fhConversionsPerDetectorPE->GetXaxis()->SetBinLabel(5, "TOF");
161
162 //fhConversionsPerDetectorPE_cut = new TH2I("fhConversionsPerDetectorPE_cut", "fhConversionsPerDetectorPE_cut;;conversions per event", 5, 0, 5, 1500, 0, 1500);
163 //fHistoList_tomography.push_back(fhConversionsPerDetectorPE_cut);
164 //fhConversionsPerDetectorPE_cut->GetXaxis()->SetBinLabel(1, "magnet");
165 //fhConversionsPerDetectorPE_cut->GetXaxis()->SetBinLabel(2, "STS");
166 //fhConversionsPerDetectorPE_cut->GetXaxis()->SetBinLabel(3, "RICH");
167 //fhConversionsPerDetectorPE_cut->GetXaxis()->SetBinLabel(4, "TRD");
168 //fhConversionsPerDetectorPE_cut->GetXaxis()->SetBinLabel(5, "TOF");
169
170 // tomography from reconstructed tracks
171 fhTomography_reco = new TH3D("fhTomography_reco", "fhTomography_reco;X in cm;Y in cm;Z in cm", 200, -200., 200., 200,
172 -200., 200., 500, 0., 1000.);
174 new TH2D("fhTomography_reco_XZ", "fhTomography_reco_XZ;X in cm;Z in cm", 1600, -400., 400., 2400, 0., 1200.);
176 new TH2D("fhTomography_reco_YZ", "fhTomography_reco_YZ;Y in cm;Z in cm", 1600, -400., 400., 2400, 0., 1200.);
177 fhConversion_reco = new TH1D("fhConversion_reco", "fhConversion_reco;Z in cm;# conversions", 4800, -0.5, 1199.5);
182}
183
184
186{
187 //gDirectory->cd("analysis-conversion");
188 gDirectory->mkdir("Tomography");
189 gDirectory->cd("Tomography");
190 for (UInt_t i = 0; i < fHistoList_tomography.size(); i++) {
191 fHistoList_tomography[i]->Write();
192 }
193 gDirectory->cd("..");
194
195 cout << "CbmAnaConversionTomography: Realtime - " << fTime << endl;
196}
197
198
200{
201 timer.Start();
202
203 electronIDs.clear();
204 electronMotherIDs.clear();
205 conversionsInDetector[0] = 0; // magnet
206 conversionsInDetector[1] = 0; // sts
207 conversionsInDetector[2] = 0; // rich
208 conversionsInDetector[3] = 0; // trd
209 conversionsInDetector[4] = 0; // tof
210
211 Int_t nofMcTracks = fMcTracks->GetEntriesFast();
212 for (int i = 0; i < nofMcTracks; i++) {
213 CbmMCTrack* mctrack = (CbmMCTrack*) fMcTracks->At(i);
214 if (mctrack == NULL) continue;
215
216 if (TMath::Abs(mctrack->GetPdgCode()) == 11) {
217 int motherId = mctrack->GetMotherId();
218 if (motherId == -1) continue;
219 // CbmMCTrack* mother = (CbmMCTrack*) fMcTracks->At(motherId);
220 // int mcMotherPdg = -1;
221 // if (NULL != mother) mcMotherPdg = mother->GetPdgCode();
222
223 electronIDs.push_back(i);
224 electronMotherIDs.push_back(motherId);
225 }
226 }
227
228
229 // int photoncounter = 0;
230 std::multimap<int, int> electronMap;
231 for (unsigned int i = 0; i < electronIDs.size(); i++) {
232 electronMap.insert(std::pair<int, int>(electronMotherIDs[i], electronIDs[i]));
233 }
234
235 // int check = 0;
236 for (std::map<int, int>::iterator it = electronMap.begin(); it != electronMap.end(); ++it) {
237 if (it == electronMap.begin()) { TomographyMC(it->second); }
238 else {
239 std::map<int, int>::iterator zwischen = it;
240 zwischen--;
241 int id = it->first;
242 int id_old = zwischen->first;
243 if (id == id_old) continue;
244 else {
245 TomographyMC(it->second);
246 }
247 }
248 }
249
250 fhConversionsPerDetectorPE->Fill(1.0 * 0, 1.0 * conversionsInDetector[0]);
251 fhConversionsPerDetectorPE->Fill(1.0 * 1, 1.0 * conversionsInDetector[1]);
252 fhConversionsPerDetectorPE->Fill(1.0 * 2, 1.0 * conversionsInDetector[2]);
253 fhConversionsPerDetectorPE->Fill(1.0 * 3, 1.0 * conversionsInDetector[3]);
254 fhConversionsPerDetectorPE->Fill(1.0 * 4, 1.0 * conversionsInDetector[4]);
255
256 //fhConversionsPerDetectorPE_cut->Fill(1.0*0, 1.0*conversionsInDetector_cut[0]);
257 //fhConversionsPerDetectorPE_cut->Fill(1.0*1, 1.0*conversionsInDetector_cut[1]);
258 //fhConversionsPerDetectorPE_cut->Fill(1.0*2, 1.0*conversionsInDetector_cut[2]);
259 //fhConversionsPerDetectorPE_cut->Fill(1.0*3, 1.0*conversionsInDetector_cut[3]);
260 //fhConversionsPerDetectorPE_cut->Fill(1.0*4, 1.0*conversionsInDetector_cut[4]);
261
262
263 timer.Stop();
264 fTime += timer.RealTime();
265}
266
267/*
268void CbmAnaConversionTomography::TomographyMC(CbmMCTrack* mctrack)
269// doing tomography from gamma -> e+ e- decays, MC DATA
270{
271 if (TMath::Abs( mctrack->GetPdgCode()) == 11) {
272 int motherId = mctrack->GetMotherId();
273 if (motherId != -1) {
274 CbmMCTrack* mother = (CbmMCTrack*) fMcTracks->At(motherId);
275 int mcMotherPdg = -1;
276 if (NULL != mother) mcMotherPdg = mother->GetPdgCode();
277
278 if (mcMotherPdg == 22) {
279 TVector3 v;
280 mctrack->GetStartVertex(v);
281 fhGammaZ->Fill(v.Z());
282
283
284 fhTomography->Fill(v.X(), v.Y(), v.Z());
285 fhTomography_XZ->Fill(v.X(), v.Z());
286 fhTomography_YZ->Fill(v.Y(), v.Z());
287 fhConversion_energy->Fill(mctrack->GetEnergy());
288 fhConversion_p->Fill(mctrack->GetP());
289
290 if(v.z() >= 0 && v.Z() <= 170) { // only in region before the RICH
291 fhTomography_uptoRICH->Fill(v.X(), v.Y());
292 }
293 if(v.z() >= 170 && v.Z() <= 400) { // only in region of the RICH detector
294 fhTomography_RICH_complete->Fill(v.X(), v.Y());
295 }
296 if(v.z() >= 220 && v.Z() <= 280) { // only in region of RICH beampipe, without being distorted by mirrors or camera
297 fhTomography_RICH_beampipe->Fill(v.X(), v.Y());
298 }
299 if(v.z() >= 100 && v.Z() <= 170) { // only in a downstream part of STS and magnet
300 fhTomography_STS_end->Fill(v.X(), v.Y());
301 }
302 if(v.z() >= 98 && v.Z() <= 102) { // only in a downstream part of STS and magnet, only last STS station
303 fhTomography_STS_lastStation->Fill(v.X(), v.Y());
304 }
305 if(v.z() >= 179 && v.Z() <= 181) { // only in region of RICH frontplate
306 fhTomography_RICH_frontplate->Fill(v.X(), v.Y());
307 }
308 if(v.z() >= 369 && v.Z() <= 371) { // only in region of RICH backplate
309 fhTomography_RICH_backplate->Fill(v.X(), v.Y());
310 }
311
312 fhConversion->Fill(v.Z());
313
314 if( (TMath::Abs(v.X()) <= 100 && TMath::Abs(v.X()) > 3) && (TMath::Abs(v.Y()) <= 50 && TMath::Abs(v.Y()) > 3) && TMath::Abs(v.Z()) <= 169) {
315 fhConversion_inSTS->Fill(v.Z());
316 }
317 }
318 }
319 }
320}
321*/
322
323
325// doing tomography from gamma -> e+ e- decays, MC DATA
326{
327
328 CbmMCTrack* mctrack = (CbmMCTrack*) fMcTracks->At(electronID);
329
330 TVector3 v;
331 mctrack->GetStartVertex(v);
332 fhGammaZ->Fill(v.Z());
333
334
335 fhTomography->Fill(v.X(), v.Y(), v.Z());
336 fhTomography_XZ->Fill(v.X(), v.Z());
337 fhTomography_YZ->Fill(v.Y(), v.Z());
338 //if(GetNPoints(mctrack)) {
339 // fhTomography_XZ_cut->Fill(v.X(), v.Z());
340 // fhTomography_YZ_cut->Fill(v.Y(), v.Z());
341 //}
342 fhConversion_energy->Fill(mctrack->GetEnergy());
343 fhConversion_p->Fill(mctrack->GetP());
344
345 if (v.z() >= 0 && v.Z() <= 170) { // only in region before the RICH
346 fhTomography_uptoRICH->Fill(v.X(), v.Y());
347 }
348 if (v.z() >= 170 && v.Z() <= 400) { // only in region of the RICH detector
349 fhTomography_RICH_complete->Fill(v.X(), v.Y());
350 }
351 if (v.z() >= 220 && v.Z() <= 280) { // only in region of RICH beampipe, without being distorted by mirrors or camera
352 fhTomography_RICH_beampipe->Fill(v.X(), v.Y());
353 }
354 if (v.z() >= 100 && v.Z() <= 170) { // only in a downstream part of STS and magnet
355 fhTomography_STS_end->Fill(v.X(), v.Y());
356 }
357 if (v.z() >= 98 && v.Z() <= 102) { // only in a downstream part of STS and magnet, only last STS station
358 fhTomography_STS_lastStation->Fill(v.X(), v.Y());
359 }
360 if (v.z() >= 179 && v.Z() <= 181) { // only in region of RICH frontplate
361 fhTomography_RICH_frontplate->Fill(v.X(), v.Y());
362 }
363 if (v.z() >= 369 && v.Z() <= 371) { // only in region of RICH backplate
364 fhTomography_RICH_backplate->Fill(v.X(), v.Y());
365 }
366
367 fhConversion->Fill(v.Z());
368 fhConversion_vs_momentum->Fill(v.Z(), mctrack->GetP());
369
370 if (mctrack->GetP() > 1) { fhConversion_cut->Fill(v.Z()); }
371
372 if ((TMath::Abs(v.X()) <= 100 && TMath::Abs(v.X()) > 3) && (TMath::Abs(v.Y()) <= 50 && TMath::Abs(v.Y()) > 3)
373 && TMath::Abs(v.Z()) <= 169) {
374 fhConversion_inSTS->Fill(v.Z());
375 }
376
377
378 if (v.z() >= 170 && v.Z() <= 380) { // only in region of the RICH detector
381 //if(GetNPoints(mctrack)) { conversionsInDetector_cut[2]++; }
382 }
383 //if(v.z() >= 405 && v.Z() <= 870) { // only in region of the TRD detector, SIS300
384 if (v.z() >= 405 && v.Z() <= 595) { // only in region of the TRD detector
387 //if(GetNPoints(mctrack)) { conversionsInDetector_cut[3]++; }
388 }
389 //if(v.z() > 870 && v.Z() <= 1010) { // only in region of the TOF detector, SIS300
390 if (v.z() > 598 && v.Z() <= 730) { // only in region of the TOF detector
393 //if(GetNPoints(mctrack)) { conversionsInDetector_cut[4]++; }
394 }
395 if ((TMath::Abs(v.X()) <= 100 && TMath::Abs(v.X()) > 3) && (TMath::Abs(v.Y()) <= 65 && TMath::Abs(v.Y()) > 3)
396 && TMath::Abs(v.Z()) <= 105) { // STS
399 //if(GetNPoints(mctrack)) { conversionsInDetector_cut[1]++; }
400 }
401 if ((TMath::Abs(v.X()) > 100 || TMath::Abs(v.Y()) > 65) && TMath::Abs(v.Z()) <= 169) { // magnet
404 //if(GetNPoints(mctrack)) { conversionsInDetector_cut[0]++; }
405 }
406}
407
408
410// doing tomography from gamma -> e+ e- decays, RECONSTRUCTED TRACKS WITH MC DATA
411{
412 timer.Start();
413
414 if (TMath::Abs(mctrack->GetPdgCode()) == 11) {
415 int motherId = mctrack->GetMotherId();
416 if (motherId != -1) {
417 CbmMCTrack* mother = (CbmMCTrack*) fMcTracks->At(motherId);
418 int mcMotherPdg = -1;
419 if (NULL != mother) mcMotherPdg = mother->GetPdgCode();
420
421 if (mcMotherPdg == 22) {
422 TVector3 v;
423 mctrack->GetStartVertex(v);
424
425 fhTomography_reco->Fill(v.X(), v.Y(), v.Z());
426 fhTomography_reco_XZ->Fill(v.X(), v.Z());
427 fhTomography_reco_YZ->Fill(v.Y(), v.Z());
428 fhConversion_reco->Fill(v.Z());
429 }
430 }
431 }
432
433 timer.Stop();
434 fTime += timer.RealTime();
435}
436
437
439{
440 Double_t np_sts = mctrack->GetNPoints(ECbmModuleId::kSts);
441 Double_t np_rich = mctrack->GetNPoints(ECbmModuleId::kRich);
442 Double_t np_trd = mctrack->GetNPoints(ECbmModuleId::kTrd);
443 Double_t np_tof = mctrack->GetNPoints(ECbmModuleId::kTof);
444
445 Bool_t result = false;
446 if (np_sts > 0 || np_rich > 0 || np_trd > 0 || np_tof > 0) { result = true; }
447
448 return result;
449}
@ kTrd
Transition Radiation Detector.
@ kTof
Time-of-flight Detector.
@ kSts
Silicon Tracking System.
@ kRich
Ring-Imaging Cherenkov Detector.
fscal v[fmask::Size]
Definition KfSimdPseudo.h:4
void TomographyReco(CbmMCTrack *mctrack)
Bool_t GetNPoints(CbmMCTrack *mctrack)
double GetP() const
Definition CbmMCTrack.h:98
int32_t GetMotherId() const
Definition CbmMCTrack.h:69
int32_t GetPdgCode() const
Definition CbmMCTrack.h:68
void GetStartVertex(TVector3 &vertex) const
Definition CbmMCTrack.h:176
int32_t GetNPoints(ECbmModuleId detId) const
double GetEnergy() const
Definition CbmMCTrack.h:162
Hash for CbmL1LinkKey.