CbmRoot
Loading...
Searching...
No Matches
CbmEventBuilderQa.cxx
Go to the documentation of this file.
1/* Copyright (C) 2017-2020 IKF-UFra, GSI
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Valentina Akishina , Maksym Zyzak, Valentina Akishina [committer] */
4
10// Cbm Headers ----------------------
11#include "CbmEventBuilderQa.h"
12
13#include "CbmEbEventEfficiencies.h"
14#include "CbmEbEventMatch.h"
15#include "CbmEbMCEvent.h"
16#include "CbmEvent.h"
17#include "CbmGlobalTrack.h"
18#include "CbmMCDataArray.h"
19#include "CbmMCDataManager.h"
20#include "CbmMCEventList.h"
21#include "CbmMCTrack.h"
22#include "CbmMuchTrack.h"
23#include "CbmRichRing.h"
24#include "CbmStsDigi.h"
25#include "CbmStsHit.h"
26#include "CbmStsPoint.h"
27#include "CbmStsSensor.h"
28#include "CbmStsSetup.h"
29#include "CbmStsStation.h"
30#include "CbmStsTrack.h"
31#include "CbmTrack.h"
32#include "CbmTrackMatchNew.h"
33
34#include <FairRunAna.h>
35
36#include <TClonesArray.h>
37#include <TDatabasePDG.h>
38#include <TDirectory.h>
39#include <TFile.h>
40#include <TH1F.h>
41#include <TH2F.h>
42#include <TMath.h>
43#include <TObject.h>
44
45#include <cmath>
46#include <iomanip>
47#include <iostream>
48
49using std::cout;
50using std::endl;
51using std::map;
52using std::vector;
53
54struct TH1FParameters {
55 TString name, title;
56 int nbins;
57 float xMin, xMax;
58};
59
60CbmEventBuilderQa::CbmEventBuilderQa(const char* name, Int_t iVerbose, TString outFileName)
61 : FairTask(name, iVerbose)
62 , fPointsInTracks()
63 , fStsTrackBranchName("StsTrack")
64 , fGlobalTrackBranchName("GlobalTrack")
65 , fRichBranchName("RichRing")
66 , fTrdBranchName("TrdTrack")
67 , fTofBranchName("TofHit")
68 , fMuchTrackBranchName("MuchTrack")
69 , fMCTracksBranchName("MCTrack")
70 , fStsTrackMatchBranchName("StsTrackMatch")
71 , fRichRingMatchBranchName("RichRingMatch")
72 , fTrdTrackMatchBranchName("TrdTrackMatch")
73 , fTofHitMatchBranchName("TofHitMatch")
74 , fMuchTrackMatchBranchName("MuchTrackMatch")
75 , fStsDigis(0)
76 , fStsTracks(0)
77 , fMCTracks(0)
78 , fStsHits(0)
79 , fMvdPoints(0)
80 , fStsPoints(0)
81 , fEvents(0)
82 , fStsTrackMatchArray(0)
83 , fStsHitMatch(0)
84 , fEventList(nullptr)
85 , fOutFileName(outFileName)
86 , fOutFile(0)
87 , fHistoDir(0) //, //fNEvents(0),
88 // fPDGtoIndexMap()
89{
90 TFile* curFile = gFile;
91 TDirectory* curDirectory = gDirectory;
92
93 if (!(fOutFileName == ""))
94 fOutFile = new TFile(fOutFileName.Data(), "RECREATE");
95 else
96 fOutFile = gFile;
97
98 fOutFile->cd();
99 fHistoDir = fOutFile->mkdir("4DQA");
100 fHistoDir->cd();
101
102 gDirectory->mkdir("TimeHisto");
103 gDirectory->cd("TimeHisto");
104 {
105 TH1FParameters timeHisto[fNTimeHistos] = {
106 {"NTracksVsTime", "NTracksVsTime", 12100, -100.f, 12000.f},
107 {"TracksTimeResidual", "TracksTimeResidual", 300, -30.f, 30.f},
108 {"TracksTimePull", "TracksTimePull", 100, -10.f, 10.f},
109 {"NHitsVsTime", "NHitsVsTime", 12100, -100.f, 12000.f},
110 {"HitsTimeResidual", "HitsTimeResidual", 300, -30.f, 30.f},
111 {"HitsTimePull", "HitsTimePull", 100, -10.f, 10.f},
112 {"NTracksInEventsVsTime", "NTracksInEventsVsTime", 12100, -100.f, 12000.f},
113 {"NHitsInEventsVsTime", "NHitsInEventsVsTime", 12100, -100.f, 12000.f},
114 {"NTracksVsMCTime", "NTracksVsMCTime", 12100, -100.f, 12000.f},
115 {"NHitsVsMCTime", "NHitsVsMCTime", 12100, -100.f, 12000.f},
116 {"dtDistribution", "dtDistribution", 100, 0.f, 50.f},
117 {"NTracksInEventsVsTime1", "NTracksInEventsVsTime1", 120100, -100.f, 12000.f},
118 {"NTracksInEventsVsTime2", "NTracksInEventsVsTime2", 120100, -100.f, 12000.f},
119 {"NTracksInEventsVsTime3", "NTracksInEventsVsTime3", 120100, -100.f, 12000.f},
120 {"NTracksInEventsVsTime4", "NTracksInEventsVsTime4", 120100, -100.f, 12000.f},
121 {"NTracksInEventsVsTime5", "NTracksInEventsVsTime5", 120100, -100.f, 12000.f},
122
123 {"NHitsInEventsVsTime1", "NHitsInEventsVsTime1", 12100, -100.f, 12000.f},
124 {"NHitsInEventsVsTime2", "NHitsInEventsVsTime2", 12100, -100.f, 12000.f},
125 {"NHitsInEventsVsTime3", "NHitsInEventsVsTime3", 12100, -100.f, 12000.f},
126 {"NHitsInEventsVsTime4", "NHitsInEventsVsTime4", 12100, -100.f, 12000.f},
127 {"NHitsInEventsVsTime5", "NHitsInEventsVsTime5", 12100, -100.f, 12000.f},
128
129 {"Number of merged events", "Number of merged events", 6, -0.5f, 5.5f},
130 {"Event length", "Event length", 100, 0.f, 50.f},
131 {"NTracksPerEvent", "NTracksPerEvent", 250, 0.f, 250.f},
132 {"TrackRecoTimePull", "TrackRecoTimePull", 100, -10.f, 10.f},
133 {"TrackTimeEvent", "TrackTimeEvent", 100, -10.f, 10.f},
134 {"TrackTimeNoEvent", "TrackTimeNoEvent", 100, -10.f, 10.f}};
135 for (int iH = 0; iH < fNTimeHistos; iH++)
136 fTimeHisto[iH] = new TH1F(timeHisto[iH].name.Data(), timeHisto[iH].title.Data(), timeHisto[iH].nbins,
137 timeHisto[iH].xMin, timeHisto[iH].xMax);
138
139 fTimeHisto[0]->GetXaxis()->SetTitle("Time, ns");
140 fTimeHisto[0]->GetYaxis()->SetTitle("Number of tracks");
141 fTimeHisto[16]->GetYaxis()->SetTitle("Number of Events");
142 fTimeHisto[16]->GetXaxis()->SetTitle("Number of MCEvents in Event");
143
144 fTimeHisto[11]->SetLineColor(5);
145 fTimeHisto[11]->SetFillColor(5);
146 fTimeHisto[12]->SetLineColor(2);
147 fTimeHisto[12]->SetFillColor(2);
148 fTimeHisto[13]->SetLineColor(3);
149 fTimeHisto[13]->SetFillColor(3);
150 fTimeHisto[14]->SetLineColor(7);
151 fTimeHisto[14]->SetFillColor(7);
152 fTimeHisto[15]->SetLineColor(6);
153 fTimeHisto[15]->SetFillColor(6);
154
155 fTimeHisto[16]->SetLineColor(5);
156 fTimeHisto[16]->SetFillColor(5);
157 fTimeHisto[17]->SetLineColor(2);
158 fTimeHisto[17]->SetFillColor(2);
159 fTimeHisto[18]->SetLineColor(3);
160 fTimeHisto[18]->SetFillColor(3);
161 fTimeHisto[19]->SetLineColor(7);
162 fTimeHisto[19]->SetFillColor(7);
163 fTimeHisto[20]->SetLineColor(6);
164 fTimeHisto[20]->SetFillColor(6);
165 }
166 gDirectory->cd(".."); //4DQA
167
168
169 gFile = curFile;
170 gDirectory = curDirectory;
171}
172
173CbmEventBuilderQa::~CbmEventBuilderQa() {}
174
175
176struct CbmBuildEventMCTrack {
177 CbmBuildEventMCTrack() : fMCFileId(-1), fMCEventId(-1), fMCTrackId(-1), fRecoTrackId(), fRecoEventId() {}
178
179 int fMCFileId;
180 int fMCEventId;
181 int fMCTrackId;
182
183 vector<int> fRecoTrackId;
184 vector<int> fRecoEventId;
185};
186
187InitStatus CbmEventBuilderQa::Init()
188{
189 //Get ROOT Manager
190 FairRootManager* ioman = FairRootManager::Instance();
191
192 if (ioman == 0) {
193 Warning("CbmEventBuilderQa::Init", "RootManager not instantiated!");
194 return kERROR;
195 }
196
197 CbmMCDataManager* mcManager = (CbmMCDataManager*) ioman->GetObject("MCDataManager");
198 if (mcManager == nullptr) LOG(fatal) << GetName() << ": No CbmMCDataManager!";
199
200 fMCTracks = (CbmMCDataArray*) mcManager->InitBranch("MCTrack");
201 if (fMCTracks == nullptr) LOG(fatal) << GetName() << ": No MCTrack data!";
202
203 fEventList = (CbmMCEventList*) ioman->GetObject("MCEventList.");
204 if (fEventList == 0) {
205 LOG(error) << GetName() << "MC Event List not found!";
206 return kERROR;
207 }
208
209 // open MCTrack array
210 fStsTracks = (TClonesArray*) ioman->GetObject("StsTrack");
211 assert(fStsTracks);
212
213 fStsTrackMatchArray = (TClonesArray*) ioman->GetObject("StsTrackMatch");
214 if (fStsTrackMatchArray == 0) {
215 LOG(error) << GetName() << "track match array not found!";
216 return kERROR;
217 }
218
219 fStsHits = (TClonesArray*) ioman->GetObject("StsHit");
220 assert(fStsHits);
221
222 fStsPoints = mcManager->InitBranch("StsPoint");
223 assert(fStsPoints);
224
225 fMvdPoints = mcManager->InitBranch("StsPoint");
226 assert(fStsPoints);
227
228 // --- Get input array (CbmStsDigi)
229 fStsHitMatch = (TClonesArray*) ioman->GetObject("StsHitMatch");
230 assert(fStsHitMatch);
231
232 // --- Get input array (CbmEvent)
233 fEvents = dynamic_cast<TClonesArray*>(ioman->GetObject("CbmEvent"));
234 if (nullptr == fEvents) {
235 LOG(fatal) << GetName() << "Init: No CbmEvent TClonesArray found!";
236 }
237
238 return kSUCCESS;
239}
240
241void CbmEventBuilderQa::Exec(Option_t* /*opt*/)
242{
243 fPointsInTracks.clear();
244
245 int nMCEvents = fEventList->GetNofEvents();
246
247
248 vector<vector<int>> mcHitMap(nMCEvents);
249 vector<CbmEbMCEvent> fMCEvents(nMCEvents);
250
251 //Fill Histo for Hits
252
253 UInt_t nHits = fStsHits->GetEntriesFast();
254
255
256 for (UInt_t iHit = 0; iHit < nHits; iHit++) {
257 CbmStsHit* hit = (CbmStsHit*) fStsHits->At(iHit);
258
259 float hit_time = hit->GetTime();
260
261 fTimeHisto[3]->Fill(hit_time);
262
263
264 CbmMatch* stsHitMatch = (CbmMatch*) fStsHitMatch->At(iHit);
265 if (stsHitMatch->GetNofLinks() == 0) continue;
266 Float_t bestWeight = 0.f;
267 Float_t totalWeight = 0.f;
268 // Int_t mcHitId = -1;
269 int mcEvent = -1;
270 int iMCPoint = -1;
271 CbmLink link;
272
273 for (int iLink = 0; iLink < stsHitMatch->GetNofLinks(); iLink++) {
274 totalWeight += stsHitMatch->GetLink(iLink).GetWeight();
275 if (stsHitMatch->GetLink(iLink).GetWeight() > bestWeight) {
276 bestWeight = stsHitMatch->GetLink(iLink).GetWeight();
277 iMCPoint = stsHitMatch->GetLink(iLink).GetIndex();
278 link = stsHitMatch->GetLink(iLink);
279 mcEvent = link.GetEntry();
280 }
281 }
282 if (bestWeight / totalWeight < 0.7 || iMCPoint < 0) continue;
283
284
285 int iCol = (mcEvent) % 5;
286
287 if (iCol >= 0) fTimeHisto[16 + iCol]->Fill(hit_time);
288
289 CbmStsPoint* stsMcPoint = (CbmStsPoint*) fStsPoints->Get(link.GetFile(), link.GetEntry(), link.GetIndex());
290 double mcTime = stsMcPoint->GetTime() + fEventList->GetEventTime(link.GetEntry() + 1, link.GetFile());
291
292 fTimeHisto[7]->Fill(mcTime);
293
294 double residual = hit_time - mcTime;
295 double pull = residual / 5.;
296
297 fTimeHisto[4]->Fill(residual);
298 fTimeHisto[5]->Fill(pull);
299 }
300
302 UInt_t nTracks = fStsTracks->GetEntriesFast();
303
304 vector<bool> IsTrackReconstructed(nTracks, 0);
305 vector<bool> IsTrackReconstructable(nTracks, 0);
306
307 const int iMCFile = 0;
308 if (fStsPoints) {
309 fPointsInTracks.resize(nMCEvents);
310 for (int iMCEvent = 0; iMCEvent < nMCEvents; iMCEvent++) {
311 const int nMCTracksInEvent = fMCTracks->Size(iMCFile, iMCEvent);
312 fPointsInTracks[iMCEvent].resize(nMCTracksInEvent);
313
314 unsigned int nStsPoints = fStsPoints->Size(iMCFile, iMCEvent);
315
316 for (unsigned int iStsPoint = 0; iStsPoint < nStsPoints; iStsPoint++) {
317 CbmStsPoint* point = (CbmStsPoint*) fStsPoints->Get(iMCFile, iMCEvent, iStsPoint);
318 const int iMCTrack = point->GetTrackID();
319 fPointsInTracks[iMCEvent][iMCTrack].push_back(iStsPoint);
320 }
321 }
322 }
323
324 for (unsigned int iTrack = 0; iTrack < nTracks; iTrack++) {
325
326 CbmStsTrack* track = (CbmStsTrack*) fStsTracks->At(iTrack);
327
328 fTimeHisto[0]->Fill(track->GetStartTime());
329
330
331 UInt_t NHits = track->GetNofStsHits();
332
333 for (UInt_t iHit = 0; iHit < NHits; iHit++) {
334
335 CbmMatch* stsHitMatch = (CbmMatch*) fStsHitMatch->At(iHit);
336 if (stsHitMatch->GetNofLinks() == 0) continue;
337 Float_t bestWeight = 0.f;
338 Float_t totalWeight = 0.f;
339 // Int_t mcHitId = -1;
340 // int mcEvent = -1;
341 int iMCPoint = -1;
342 CbmLink link;
343 for (int iLink = 0; iLink < stsHitMatch->GetNofLinks(); iLink++) {
344 totalWeight += stsHitMatch->GetLink(iLink).GetWeight();
345 if (stsHitMatch->GetLink(iLink).GetWeight() > bestWeight) {
346 bestWeight = stsHitMatch->GetLink(iLink).GetWeight();
347 iMCPoint = stsHitMatch->GetLink(iLink).GetIndex();
348 link = stsHitMatch->GetLink(iLink);
349 }
350 }
351 if (bestWeight / totalWeight < 0.7 || iMCPoint < 0) continue;
352
353
354 // CbmStsPoint* stsMcPoint = (CbmStsPoint*) fStsPoints->Get(link.GetFile(),link.GetEntry(),link.GetIndex());
355 // double mcTime = stsMcPoint->GetTime() + fEventList->GetEventTime(link.GetEntry()+1, link.GetFile());
356
357 //fTimeHisto[8]->Fill(mcTime);
358 }
359
360 CbmTrackMatchNew* stsTrackMatch = (CbmTrackMatchNew*) fStsTrackMatchArray->At(iTrack);
361 if (stsTrackMatch->GetNofLinks() == 0) continue;
362 Float_t bestWeight = 0.f;
363 Float_t totalWeight = 0.f;
364 Int_t mcTrackId = -1;
365 int mcEvent = -1;
366 CbmLink link;
367 for (int iLink = 0; iLink < stsTrackMatch->GetNofLinks(); iLink++) {
368 totalWeight += stsTrackMatch->GetLink(iLink).GetWeight();
369 if (stsTrackMatch->GetLink(iLink).GetWeight() > bestWeight) {
370 bestWeight = stsTrackMatch->GetLink(iLink).GetWeight();
371 int iMCTrack = stsTrackMatch->GetLink(iLink).GetIndex();
372 link = stsTrackMatch->GetLink(iLink);
373
374
375 mcEvent = link.GetEntry();
376 mcTrackId = iMCTrack;
377 }
378 }
379 if (bestWeight / totalWeight < 0.7 || mcTrackId < 0) continue;
380
381 IsTrackReconstructed[iTrack] = 1;
382
383 CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->Get(0, mcEvent, mcTrackId);
384
385 double mcTime = mcTrack->GetStartT() + fEventList->GetEventTime(mcEvent + 1, link.GetFile());
386 double residual = track->GetStartTime() - mcTime;
387 double pull = residual / track->GetStartTimeError();
388 fTimeHisto[1]->Fill(residual);
389 fTimeHisto[2]->Fill(pull);
390 }
391
392
393 vector<vector<int>> mcTrackMap(nMCEvents);
394 vector<CbmBuildEventMCTrack> mcTracks;
395
396 float tEvent = 0.f;
397 // int eventId = 0;
398
399
400 for (int iMCEvent = 0; iMCEvent < nMCEvents; iMCEvent++) {
401 unsigned int nMCTracks = fMCTracks->Size(0, iMCEvent);
402
403 mcTrackMap[iMCEvent].resize(nMCTracks);
404
405 fTimeHisto[10]->Fill(fEventList->GetEventTime(iMCEvent + 1, 0)
406 - tEvent); // +fEventList->GetEventTime(iMCEvent+1, 0) );
407 tEvent = fEventList->GetEventTime(iMCEvent + 1, 0);
408
409 int nReconstructableMCTracks = 0;
410 for (unsigned int iMC = 0; iMC < nMCTracks; iMC++) {
411 CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->Get(0, iMCEvent, iMC);
412
413 if (CalculateIsReconstructable(0, iMCEvent, iMC)) nReconstructableMCTracks++;
414
415 fTimeHisto[8]->Fill(mcTrack->GetStartT() + fEventList->GetEventTime(iMCEvent + 1, 0));
416
417 CbmBuildEventMCTrack newMCTrack;
418 newMCTrack.fMCFileId = 0;
419 newMCTrack.fMCEventId = iMCEvent;
420 newMCTrack.fMCTrackId = iMC;
421
422 mcTrackMap[iMCEvent][iMC] = mcTracks.size();
423 vector<int>& trackIds = fMCEvents[iMCEvent].GetMCTrackIds();
424 trackIds.push_back(mcTracks.size());
425
426 mcTracks.push_back(newMCTrack);
427
428 // if(mcTrack->GetMotherId()<eventId)
429 // {
430 // fTimeHisto[10]->Fill( mcTrack->GetStartT() - tEvent);// +fEventList->GetEventTime(iMCEvent+1, 0) );
431 // tEvent = mcTrack->GetStartT();
432 // eventId = mcTrack->GetMotherId();
433 // }
434 }
435 if (nReconstructableMCTracks > 2) fMCEvents[iMCEvent].SetReconstructable(1);
436 }
437
438 for (int iTrack = 0; iTrack < fStsTrackMatchArray->GetEntriesFast(); iTrack++) {
439 CbmTrackMatchNew* stsTrackMatch = (CbmTrackMatchNew*) fStsTrackMatchArray->At(iTrack);
440 if (stsTrackMatch->GetNofLinks() == 0) continue;
441 Float_t bestWeight = 0.f;
442 Float_t totalWeight = 0.f;
443 Int_t mcTrackId = -1;
444 int mcEvent = -1;
445 CbmLink link;
446 for (int iLink = 0; iLink < stsTrackMatch->GetNofLinks(); iLink++) {
447 totalWeight += stsTrackMatch->GetLink(iLink).GetWeight();
448 if (stsTrackMatch->GetLink(iLink).GetWeight() > bestWeight) {
449 bestWeight = stsTrackMatch->GetLink(iLink).GetWeight();
450 int iMCTrack = stsTrackMatch->GetLink(iLink).GetIndex();
451 link = stsTrackMatch->GetLink(iLink);
452 mcEvent = link.GetEntry();
453 mcTrackId = iMCTrack;
454 }
455 }
456 if (bestWeight / totalWeight < 0.7 || mcTrackId < 0) continue;
457
458 fMCEvents[mcEvent].GetRecoTrackIds().push_back(iTrack);
459 }
460
461 vector<CbmEbEventMatch> eventMatches;
462
463 for (int iEvent = 0; iEvent < fEvents->GetEntriesFast(); iEvent++) {
464 CbmEvent* event = (CbmEvent*) fEvents->At(iEvent);
465
466
467 vector<int> tracksId;
468 Int_t nEventTracks = event->GetNofData(ECbmDataType::kStsTrack);
469
470 fTimeHisto[23]->Fill(nEventTracks);
471
472 CbmEbEventMatch EventMatch;
473 EventMatch.SetNEventTracks(nEventTracks);
474 vector<int> tracks;
475
476 float tLastTrack = 0;
477 float tFirstTrack = 100000000000000;
478
479 for (int iTr = 0; iTr < nEventTracks; iTr++) {
480 const int stsTrackIndex = event->GetStsTrackIndex(iTr);
481
482 tracks.push_back(stsTrackIndex);
483
484 CbmTrackMatchNew* stsTrackMatch = (CbmTrackMatchNew*) fStsTrackMatchArray->At(stsTrackIndex);
485 if (stsTrackMatch->GetNofLinks() == 0) continue;
486 Float_t bestWeight = 0.f;
487 Float_t totalWeight = 0.f;
488 Int_t mcTrackId = -1;
489 int mcEvent = -1;
490 CbmLink link;
491 for (int iLink = 0; iLink < stsTrackMatch->GetNofLinks(); iLink++) {
492 totalWeight += stsTrackMatch->GetLink(iLink).GetWeight();
493 if (stsTrackMatch->GetLink(iLink).GetWeight() > bestWeight) {
494 bestWeight = stsTrackMatch->GetLink(iLink).GetWeight();
495 int iMCTrack = stsTrackMatch->GetLink(iLink).GetIndex();
496 link = stsTrackMatch->GetLink(iLink);
497
498
499 mcEvent = link.GetEntry();
500 mcTrackId = iMCTrack;
501 }
502 }
503 if (bestWeight / totalWeight < 0.7 || mcTrackId < 0) continue;
504
505 CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->Get(0, mcEvent, mcTrackId);
506
507 CbmStsTrack* track = (CbmStsTrack*) fStsTracks->At(stsTrackIndex);
508
509 if ((track->GetStartTime()) > tLastTrack) tLastTrack = (track->GetStartTime());
510 if ((track->GetStartTime()) < tFirstTrack) tFirstTrack = (track->GetStartTime());
511
512 double mcTime = mcTrack->GetStartT() + fEventList->GetEventTime(mcEvent + 1, link.GetFile());
513
514 fTimeHisto[24]->Fill(fabs(track->GetStartTime() - mcTime) / track->GetStartTimeError());
515
516 if (fMCEvents[mcEvent].IsReconstructable()) EventMatch.AddTrack(mcEvent);
517
518 for (int iTr1 = iTr + 1; iTr1 < nEventTracks; iTr1++) {
519 if (iTr1 == iTr) continue;
520 const int stsTrackIndex1 = event->GetStsTrackIndex(iTr1);
521 CbmStsTrack* track1 = (CbmStsTrack*) fStsTracks->At(stsTrackIndex1);
522
523 fTimeHisto[25]->Fill(fabs(track->GetStartTime() - track1->GetStartTime())
524 / sqrt(track->GetStartTimeError() + track1->GetStartTimeError()));
525 }
526
527 for (int iEvent1 = 0; iEvent1 < fEvents->GetEntriesFast(); iEvent1++) {
528
529 if (iEvent == iEvent1) continue;
530 CbmEvent* event1 = (CbmEvent*) fEvents->At(iEvent1);
531
532 Int_t nTracks1 = event1->GetNofData(ECbmDataType::kStsTrack);
533 for (int iTr1 = 0; iTr1 < nTracks1; iTr1++) {
534 const int stsTrackIndex1 = event1->GetStsTrackIndex(iTr1);
535 CbmStsTrack* track1 = (CbmStsTrack*) fStsTracks->At(stsTrackIndex1);
536
537 if (fabs(track->GetStartTime() - track1->GetStartTime()) > 8.5) continue;
538
539 fTimeHisto[26]->Fill(fabs(track->GetStartTime() - track1->GetStartTime())
540 / sqrt(track->GetStartTimeError() + track1->GetStartTimeError()));
541 }
542 }
543 }
544
545 EventMatch.SetTracks(tracks);
546 fTimeHisto[22]->Fill(tLastTrack - tFirstTrack);
547 eventMatches.push_back(EventMatch);
548
549 for (std::map<int, int>::iterator it = EventMatch.GetMCEvents().begin(); it != EventMatch.GetMCEvents().end();
550 ++it) {
551
552 fMCEvents[it->first].AddRecoEvent(iEvent);
553 }
554 }
555
556
557 vector<SortEvents> Ev;
558
559 for (int iMCEvent = 0; iMCEvent < nMCEvents; iMCEvent++) {
560 int nMcPoints = fStsPoints->Size(0, iMCEvent);
561 for (int i = 0; i < nMcPoints; i++) {
562 CbmStsPoint* stsMcPoint = (CbmStsPoint*) fStsPoints->Get(0, iMCEvent, i);
563 double mcTime = stsMcPoint->GetTime() + fEventList->GetEventTime(iMCEvent + 1, 0);
564 fTimeHisto[9]->Fill(mcTime);
565 }
566 }
567
568
569 int nEvents = fEvents->GetEntriesFast();
570
571 for (int k = 0; k < nEvents; k++) {
572 SortEvents Event1;
573
574 CbmEvent* event = (CbmEvent*) fEvents->At(k);
575
576 Event1.Event = event;
577 const int stsTrackIndex = event->GetStsTrackIndex(0);
578 CbmStsTrack* track = (CbmStsTrack*) fStsTracks->At(stsTrackIndex);
579 Event1.track = *track;
580 Ev.push_back(Event1);
581 }
582
583
584 std::sort(Ev.begin(), Ev.end(), CompareTrackTime);
585
586 for (int k = 0; k < nEvents; k++) {
587 CbmEvent* event = Ev[k].Event; //(CbmEvent*) fEvents->At(k);
588
589 int iCol = k % 5;
590 for (size_t i = 0; i < event->GetNofData(ECbmDataType::kStsTrack); i++) {
591 const int stsTrackIndex = event->GetStsTrackIndex(i);
592 CbmStsTrack* track = (CbmStsTrack*) fStsTracks->At(stsTrackIndex);
593
594 fTimeHisto[6]->Fill(track->GetStartTime());
595
596
597 fTimeHisto[11 + iCol]->Fill(track->GetStartTime());
598
599 UInt_t NHits = track->GetNofStsHits();
600
601 for (UInt_t iHit = 0; iHit < NHits; iHit++) {
602 CbmStsHit* hit = (CbmStsHit*) fStsHits->At(track->GetStsHitIndex(iHit));
603 fTimeHisto[7]->Fill(hit->GetTime());
604 }
605 }
606 }
607
608
609 CbmEbEventEfficiencies fEventEfficiency;
610
611 CbmEbEventEfficiencies eventEfficiency; // efficiencies for current event
612 for (unsigned int iEvent = 0; iEvent < eventMatches.size(); iEvent++)
613 eventEfficiency.AddGhost(eventMatches[iEvent].IsGhost());
614
615 for (unsigned int iMCEvent = 0; iMCEvent < fMCEvents.size(); iMCEvent++) {
616 if (fMCEvents[iMCEvent].IsReconstructable() && fMCEvents[iMCEvent].NMCTracks() > 1) {
617 const vector<int>& recoEvents = fMCEvents[iMCEvent].GetRecoEvents(); // for length calculations
618 bool reco = 0;
619 if (recoEvents.size() == 1) reco = 1;
620 // number of cloned events
621 int nclones = 0;
622 if (fMCEvents[iMCEvent].NClones() > 0) nclones = 1;
623 eventEfficiency.Inc(reco, nclones, "reconstructable");
624 }
625 // if( fMCEvents[iMCEvent].IsReconstructed() )
626 if (fMCEvents[iMCEvent].GetRecoTrackIds().size() > 2) {
627 const vector<int>& recoEvents = fMCEvents[iMCEvent].GetRecoEvents(); // for length calculations
628 bool reco = 0;
629 if (recoEvents.size() == 1) reco = 1;
630 // number of cloned events
631 int nclones = 0;
632 if (fMCEvents[iMCEvent].NClones() > 0) nclones = 1;
633 eventEfficiency.Inc(reco, nclones, "reconstructed");
634 }
635 } // for mcTracks
636 eventEfficiency.IncNEvents();
637 fEventEfficiency += eventEfficiency;
638 eventEfficiency.CalcEff();
639 fEventEfficiency.CalcEff();
640
641 std::cout << "Event reconstruction efficiency" << std::endl;
642 fEventEfficiency.PrintEff();
643 cout << endl;
644
645
646 // int iDMcTrack;
647 for (unsigned int iM = 0; iM < eventMatches.size(); iM++) {
648 // cout<<eventMatches[iM].NMCEvents()<<" Nevents"<<endl;
649
650 fTimeHisto[21]->Fill(eventMatches[iM].NMCEvents());
651 }
652}
653
654bool CbmEventBuilderQa::CalculateIsReconstructable(const int iMCFile, const int iMCEvent, const int iMCTrack)
655{
656
657 CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->Get(iMCFile, iMCEvent, iMCTrack);
658
659 bool f = 1;
660
661 // reject very slow tracks from analysis
662 f &= (mcTrack->GetP() > 0.1);
663 // detected at least in 4 stations
664 // f &= (nMCContStations >= 4);
665
666 int maxNStaMC = 0; // max number of mcPoints on station
667 int maxNSensorMC = 0; // max number of mcPoints with same z
668 int nMCStations = 0;
669
670 maxNStaMC = 0;
671 maxNSensorMC = 0;
672 nMCStations = 0;
673 int lastSta = -1;
674 float lastz = -1;
675 int nMCContStations = 0;
676 int istaold = -1, ncont = 0;
677 int cur_maxNStaMC = 0, cur_maxNSensorMC = 0;
678 for (unsigned int iH = 0; iH < fPointsInTracks[iMCEvent][iMCTrack].size(); iH++) {
679
680 int iStsPoint = fPointsInTracks[iMCEvent][iMCTrack][iH];
681 CbmStsPoint* point = (CbmStsPoint*) fStsPoints->Get(iMCFile, iMCEvent, iStsPoint);
682
683 int currentStation = -1;
684 for (int iStation = 0; iStation < CbmStsSetup::Instance()->GetNofStations(); iStation++) {
685 CbmStsStation* station = CbmStsSetup::Instance()->GetStation(iStation);
686 const float zStation = station->GetZ();
687
688 if (fabs(point->GetZ() - zStation) < 2.5) {
689 currentStation = iStation;
690 break;
691 }
692 }
693 if (currentStation < 0) continue;
694
695 if (currentStation == lastSta)
696 cur_maxNStaMC++;
697 else { // new station
698 if (cur_maxNStaMC > maxNStaMC) maxNStaMC = cur_maxNStaMC;
699 cur_maxNStaMC = 1;
700 lastSta = currentStation;
701 nMCStations++;
702 }
703
704 if (point->GetZ() == lastz) // TODO: works incorrect because of different z
705 cur_maxNSensorMC++;
706 else { // new z
707 if (cur_maxNSensorMC > maxNSensorMC) maxNSensorMC = cur_maxNSensorMC;
708 cur_maxNSensorMC = 1;
709 lastz = point->GetZ();
710 }
711
712 int ista = currentStation;
713 if (ista - istaold == 1)
714 ncont++;
715 else if (ista - istaold > 1) {
716 if (nMCContStations < ncont) nMCContStations = ncont;
717 ncont = 1;
718 }
719 if (ista <= istaold) continue; // backward direction
720 istaold = ista;
721 };
722
723 if (nMCContStations < ncont) nMCContStations = ncont;
724
725 if (cur_maxNStaMC > maxNStaMC) maxNStaMC = cur_maxNStaMC;
726 if (cur_maxNSensorMC > maxNSensorMC) maxNSensorMC = cur_maxNSensorMC;
727 // cout << pdg << " " << p << " " << Points.size() << " > " << maxNStaMC << " " << maxNSensorMC << endl;
728
729
730 // maximul 4 layers for a station.
731 // f &= (maxNStaHits <= 4);
732 f &= (maxNStaMC <= 4);
733
734 bool isReconstructableTrack = f & (nMCContStations >= 4); // L1-MC
735
736 return (isReconstructableTrack);
737};
738
739void CbmEventBuilderQa::Finish()
740{
741 TDirectory* curr = gDirectory;
742 TFile* currentFile = gFile;
743 // Open output file and write histograms
744
745 fOutFile->cd();
746 WriteHistosCurFile(fHistoDir);
747 if (!(fOutFileName == "")) {
748 fOutFile->Close();
749 fOutFile->Delete();
750 }
751 gFile = currentFile;
752 gDirectory = curr;
753}
754
755void CbmEventBuilderQa::WriteHistosCurFile(TObject* obj)
756{
757
758 if (!obj->IsFolder())
759 obj->Write();
760 else {
761 TDirectory* cur = gDirectory;
762 TFile* currentFile = gFile;
763
764 TDirectory* sub = cur->GetDirectory(obj->GetName());
765 sub->cd();
766 TList* listSub = (static_cast<TDirectory*>(obj))->GetList();
767 TIter it(listSub);
768 while (TObject* obj1 = it())
769 WriteHistosCurFile(obj1);
770 cur->cd();
771 gFile = currentFile;
772 gDirectory = cur;
773 }
774}
775
776
777ClassImp(CbmEventBuilderQa);
tracks
TClonesArray * tracks
Definition Analyze_matching.h:21
ClassImp
ClassImp(CbmEventBuilderQa)
nMCTracks
Int_t nMCTracks
Definition CbmHadronAnalysis.cxx:72
CbmStsHit.h
Data class for a reconstructed hit in the STS.
CbmStsTrack.h
Data class for STS tracks.
mcTracks
static vector< vector< QAMCTrack > > mcTracks
Definition CbmTofHitFinderTBQa.cxx:106
size
static constexpr size_t size()
Definition KfSimdPseudo.h:2
CbmEbEventEfficiencies::Inc
void Inc(bool isReco, int _nclones, std::string name)
Definition CbmEbEventEfficiencies.h:70
CbmEbEventMatch::SetTracks
void SetTracks(vector< int > tracks)
Definition CbmEbEventMatch.h:30
CbmEbEventMatch::SetNEventTracks
void SetNEventTracks(int ntracks)
Definition CbmEbEventMatch.h:29
CbmEbEventMatch::GetMCEvents
map< int, int > & GetMCEvents()
Definition CbmEbEventMatch.h:36
CbmEbEventMatch::AddTrack
void AddTrack(int mcEventId)
Definition CbmEbEventMatch.h:28
CbmEventBuilderQa::Exec
virtual void Exec(Option_t *opt)
Definition CbmEventBuilderQa.cxx:241
CbmEventBuilderQa::Init
virtual InitStatus Init()
Definition CbmEventBuilderQa.cxx:187
CbmEventBuilderQa::fMvdPoints
CbmMCDataArray * fMvdPoints
Input array (class CbmStsDigi)
Definition CbmEventBuilderQa.h:98
CbmEventBuilderQa::CalculateIsReconstructable
bool CalculateIsReconstructable(const int iMCFile, const int iMCEvent, const int iMCTrack)
Definition CbmEventBuilderQa.cxx:654
CbmEventBuilderQa::fMCTracks
CbmMCDataArray * fMCTracks
Input array (class CbmStsDigi)
Definition CbmEventBuilderQa.h:96
CbmEventBuilderQa::fNTimeHistos
static const int fNTimeHistos
Definition CbmEventBuilderQa.h:90
CbmEventBuilderQa::fStsTrackMatchArray
TClonesArray * fStsTrackMatchArray
Definition CbmEventBuilderQa.h:101
CbmEventBuilderQa::fStsHitMatch
TClonesArray * fStsHitMatch
Definition CbmEventBuilderQa.h:102
CbmEventBuilderQa::WriteHistosCurFile
void WriteHistosCurFile(TObject *obj)
Definition CbmEventBuilderQa.cxx:755
CbmEventBuilderQa::CbmEventBuilderQa
CbmEventBuilderQa(const char *name="CbmEventBuilderQa", Int_t iVerbose=0, TString outFileName="CbmEventBuilderQa.root")
Definition CbmEventBuilderQa.cxx:60
CbmEventBuilderQa::fEventList
CbmMCEventList * fEventList
Definition CbmEventBuilderQa.h:104
CbmEventBuilderQa::fEvents
TClonesArray * fEvents
Output array (class CbmEvent)
Definition CbmEventBuilderQa.h:100
CbmEventBuilderQa::fStsTracks
TClonesArray * fStsTracks
Input array (class CbmStsDigi)
Definition CbmEventBuilderQa.h:95
CbmEventBuilderQa::fPointsInTracks
std::vector< std::vector< std::vector< int > > > fPointsInTracks
Definition CbmEventBuilderQa.h:56
CbmEventBuilderQa::CompareTrackTime
static bool CompareTrackTime(const SortEvents &a, const SortEvents &b)
Definition CbmEventBuilderQa.h:66
CbmEventBuilderQa::fStsPoints
CbmMCDataArray * fStsPoints
Input array (class CbmStsDigi)
Definition CbmEventBuilderQa.h:99
CbmEventBuilderQa::fTimeHisto
TH1F * fTimeHisto[fNTimeHistos]
Definition CbmEventBuilderQa.h:91
CbmEventBuilderQa::fStsHits
TClonesArray * fStsHits
Input array (class CbmStsDigi)
Definition CbmEventBuilderQa.h:97
CbmEvent
Class characterising one event by a collection of links (indices) to data objects,...
Definition CbmEvent.h:34
CbmEvent::GetNofData
size_t GetNofData() const
Definition CbmEvent.cxx:53
CbmEvent::GetStsTrackIndex
int32_t GetStsTrackIndex(int32_t iTrack)
Definition CbmEvent.h:128
CbmHit::GetTime
double GetTime() const
Definition CbmHit.h:76
CbmMCDataArray
Access to a MC data branch for time-based analysis.
Definition CbmMCDataArray.h:39
CbmMCDataArray::Get
TObject * Get(const CbmLink *lnk)
Definition CbmMCDataArray.h:51
CbmMCDataArray::Size
Int_t Size(Int_t fileNumber, Int_t eventNumber)
Definition CbmMCDataArray.cxx:137
CbmMCDataManager
Task class creating and managing CbmMCDataArray objects.
Definition CbmMCDataManager.h:31
CbmMCDataManager::GetObject
CbmMCDataObject * GetObject(const char *name)
Definition CbmMCDataManager.cxx:147
CbmMCDataManager::InitBranch
CbmMCDataArray * InitBranch(const char *name)
Definition CbmMCDataManager.cxx:115
CbmMCEventList
Container class for MC events with number, file and start time.
Definition CbmMCEventList.h:43
CbmMCEventList::GetNofEvents
std::size_t GetNofEvents() const
Number of events in the list @value Number of events.
Definition CbmMCEventList.h:116
CbmMCEventList::GetEventTime
double GetEventTime(uint32_t event, uint32_t file)
Event start time.
Definition CbmMCEventList.cxx:89
CbmMCTrack::GetStartT
double GetStartT() const
Definition CbmMCTrack.h:76
CbmMCTrack::GetP
double GetP() const
Definition CbmMCTrack.h:98
CbmMatch::GetLink
const CbmLink & GetLink(int32_t i) const
Definition CbmMatch.h:39
CbmMatch::GetNofLinks
int32_t GetNofLinks() const
Definition CbmMatch.h:42
CbmStsHit
data class for a reconstructed 3-d hit in the STS
Definition CbmStsHit.h:35
CbmStsSetup::GetNofStations
Int_t GetNofStations() const
Definition CbmStsSetup.h:94
CbmStsSetup::Instance
static CbmStsSetup * Instance()
Definition CbmStsSetup.cxx:303
CbmStsSetup::GetStation
CbmStsStation * GetStation(Int_t stationId) const
Definition CbmStsSetup.h:101
CbmStsStation
Class representing a station of the StsSystem.
Definition CbmStsStation.h:33
CbmStsStation::GetZ
Double_t GetZ() const
Definition CbmStsStation.h:129
CbmStsTrack::GetStsHitIndex
int32_t GetStsHitIndex(int32_t iHit) const
Definition CbmStsTrack.h:101
CbmStsTrack::GetNofStsHits
int32_t GetNofStsHits() const
Definition CbmStsTrack.h:93
CbmTrack::GetStartTime
double GetStartTime() const
Definition CbmTrack.h:71
CbmTrack::GetStartTimeError
double GetStartTimeError() const
Definition CbmTrack.h:72