CbmRoot
Loading...
Searching...
No Matches
CbmTofExtendTracks.cxx
Go to the documentation of this file.
1/* Copyright (C) 2020-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Norbert Herrmann [committer] */
4
5// -------------------------------------------------------------------------
6// ----- CbmTofExtendTracks source file -----
7// ----- Created 22/12/20 by N. Herrmann -----
8// -------------------------------------------------------------------------
9
10#include "CbmTofExtendTracks.h"
11
12#include "CbmDefs.h"
13#include "CbmEvent.h"
14#include "CbmMatch.h"
15#include "CbmTofAddress.h" // in cbmdata/tof
16#include "CbmTofCalibrator.h"
17#include "CbmTofCell.h" // in tof/TofData
18#include "CbmTofCreateDigiPar.h" // in tof/TofTools
19#include "CbmTofDetectorId_v12b.h" // in cbmdata/tof
20#include "CbmTofDetectorId_v14a.h" // in cbmdata/tof
21#include "CbmTofDetectorId_v21a.h" // in cbmdata/tof
22#include "CbmTofDigiBdfPar.h" // in tof/TofParam
23#include "CbmTofDigiPar.h" // in tof/TofParam
24#include "CbmTofGeoHandler.h" // in tof/TofTools
25#include "CbmTofHit.h"
26#include "CbmTofTrackFinder.h"
27#include "CbmTofTrackFinderNN.h"
28#include "CbmTofTracklet.h"
29#include "CbmTofTrackletParam.h"
30#include "CbmTofTrackletTools.h"
31#include "CbmVertex.h"
32#include "FairRootFileSink.h"
33#include "FairRootManager.h"
34#include "FairRunAna.h"
35#include "FairRuntimeDb.h"
36#include "LKFMinuit.h"
37#include "Logger.h"
38#include "TClonesArray.h"
39#include "TDirectory.h"
40#include "TF1.h"
41#include "TFile.h"
42#include "TFitResult.h"
43#include "TGeoManager.h"
44#include "TH1.h"
45#include "TH2.h"
46#include "TH3.h"
47#include "TMath.h"
48#include "TProfile.h"
49#include "TROOT.h"
50#include "TRandom.h"
51#include "TSpectrum.h"
52#include "TString.h"
53
54#include <iostream>
55#include <vector>
56
57using std::cout;
58using std::endl;
59using std::vector;
60
61const Int_t NDefSetup = 100;
62const double dStDZ = 1.;
64static Int_t fiTS = 0;
65static Double_t dSUT_z = 0.;
66
68
70
71// ----- Default constructor -------------------------------------------
73{
74 if (!fInstance) fInstance = this;
75}
76// -------------------------------------------------------------------------
77
78
79// ----- Standard constructor ------------------------------------------
80CbmTofExtendTracks::CbmTofExtendTracks(const char* name, const char* /*title*/, CbmTofTrackFinder* finder)
81 : FairTask(name)
82 , fFinder(finder)
83 , fTrackletTools(NULL)
84 , fTofCalibrator(NULL)
85 , fEventsColl(NULL)
86 , fTofHitArrayIn(NULL)
87 , fStsHitArrayIn(NULL)
88 , fMuchHitArrayIn(NULL)
89 , fRichHitArrayIn(NULL)
90 , fTofMatchArrayIn(NULL)
91 , fTofHitArray(NULL)
92 , fTofTrackArrayIn(NULL)
93 , fTrackArrayOut(nullptr)
94 , fvTofHitIndex()
95 , fvTofTrackIndex()
96 , fvStsHitIndex()
97 , fvMuchHitIndex()
98 , fvRichHitIndex()
99 , fvTofStationZ()
100 , fvStsStationZ()
101 , fvMuchStationZ()
102 , fvRichStationZ()
103 , fvAllHitPointer()
104 , fvTrkCalHits()
105 , fvTrkPar()
106 , fvToff()
107 , fvXoff()
108 , fvYoff()
109 , fvZoff()
110 , fvTsig()
111 , fvXsig()
112 , fvYsig()
113 , fvZsig()
114 , fhMulCorTrkTof(NULL)
115 , fhMulCorTrkSts(NULL)
116 , fhMulCorTrkMuch(NULL)
117 , fhMulCorTrkRich(NULL)
118 , fhPosCorTrkTof(NULL)
119 , fhPosCorTrkSts(NULL)
120 , fhPosCorTrkMuch(NULL)
121 , fhPosCorTrkRich(NULL)
122 , fhTrkStationDX()
123 , fhTrkStationDY()
124 , fhTrkStationDZ()
125 , fhTrkStationDT()
126 , fhTrkStationNHits()
127 , fhTrkStationDXDY()
128 , fhTrkPullDX()
129 , fhTrkPullDY()
130 , fhTrkPullDT()
131 , fhExt_Toff(NULL)
132 , fhExt_Xoff(NULL)
133 , fhExt_Yoff(NULL)
134 , fhExt_Zoff(NULL)
135 , fhExt_Tsig(NULL)
136 , fhExt_Xsig(NULL)
137 , fhExt_Ysig(NULL)
138 , fhExt_Zsig(NULL)
139 , fhExt_TrkSizVel()
140 , fhExt_TrkSizChiSq()
141 , fhExtSutXY_Found(NULL)
142 , fhExtSutXY_Missed(NULL)
143 , fhExtSutXY_DX(NULL)
144 , fhExtSutXY_DY(NULL)
145 , fhExtSutXY_DT(NULL)
146 , fVTXNorm(0.)
147 , fVTX_T(0.)
148 , fVTX_X(0.)
149 , fVTX_Y(0.)
150 , fVTX_Z(0.)
151 , fT0MAX(0.5)
152 , fiTrkHitsMin(0)
153 , fdTrkCutDX(100.)
154 , fdTrkCutDY(100.)
155 , fdTrkCutDT(100.)
156 , fdChi2Max(5.)
157 , fiCorSrc(0)
158 , fiCorMode(0)
159 , fiAddStations(0)
160 , fiReqStations(-1)
161 , fiStationUT(-1)
162 , fiCutStationMaxHitMul(1000)
163 , fiNTrkTofMax(1000)
164 , fiEvent(0)
165{
166 if (!fInstance) fInstance = this;
167}
168// -------------------------------------------------------------------------
169
170
171// ----- Destructor ----------------------------------------------------
176// -------------------------------------------------------------------------
177
178
179// ----- Public method Init (abstract in base class) --------------------
181{
182
183 fTrackletTools = new CbmTofTrackletTools(); // initialize tools
184
185 // Get and check FairRootManager
186 FairRootManager* ioman = FairRootManager::Instance();
187 if (!ioman) {
188 cout << "-E- CbmTofExtendTracks::Init: "
189 << "RootManager not instantiated!" << endl;
190 return kFATAL;
191 }
192
193 ioman->InitSink();
194
195 fEventsColl = dynamic_cast<TClonesArray*>(ioman->GetObject("Event"));
196 if (!fEventsColl) fEventsColl = dynamic_cast<TClonesArray*>(ioman->GetObject("CbmEvent"));
197 if (!fEventsColl) {
198 LOG(fatal) << "CbmEvent not found in input file";
199 }
200
201 // Get TOF hit Array
202 fTofHitArrayIn = (TClonesArray*) ioman->GetObject("TofCalHit");
203 if (!fTofHitArrayIn) {
204 fTofHitArrayIn = (TClonesArray*) ioman->GetObject("TofHit");
205 if (!fTofHitArrayIn) {
206 LOG(fatal) << "-W- CbmTofExtendTracks::Init: No TofHit array!";
207 return kERROR;
208 }
209 }
210 else {
211 LOG(info) << "-I- CbmTofExtendTracks::Init: TofCalHit array!";
212 }
213
214 // Get TOF Tracklet Array
215 fTofTrackArrayIn = (TClonesArray*) ioman->GetObject("TofTracklets");
216 if (!fTofTrackArrayIn) {
217 LOG(fatal) << "-W- CbmTofExtendTracks::Init: No TofTrack array!";
218 return kERROR;
219 }
220
221 // Get Sts hit Array
222 fStsHitArrayIn = (TClonesArray*) ioman->GetObject("StsHit");
223 if (!fStsHitArrayIn) {
224 LOG(warn) << "-W- CbmTofExtendTracks::Init: No StsHit array!";
225 //return kERROR;
226 }
227
228 // Get Much hit Array
229 fMuchHitArrayIn = (TClonesArray*) ioman->GetObject("MuchPixelHit");
230 if (!fMuchHitArrayIn) {
231 LOG(warn) << "-W- CbmTofExtendTracks::Init: No MuchHit array!";
232 //return kERROR;
233 }
234
235 // Get Rich hit Array
236 fRichHitArrayIn = (TClonesArray*) ioman->GetObject("RichHit");
237 if (!fRichHitArrayIn) {
238 LOG(warn) << "-W- CbmTofExtendTracks::Init: No RichHit array!";
239 //return kERROR;
240 }
241
242 if (kFALSE == InitParameters()) return kFATAL;
243
245
247
248 LOG(info) << "CbmTofExtendTracks initialized";
249 return kSUCCESS;
250}
251// -------------------------------------------------------------------------
252/************************************************************************************/
254{
255 UInt_t NSt = fMapStationZ.size();
256 fvToff.resize(NSt);
257 for (uint i = 0; i < NSt; i++)
258 fvToff[i] = 0.;
259 fvXoff.resize(NSt);
260 for (uint i = 0; i < NSt; i++)
261 fvXoff[i] = 0.;
262 fvYoff.resize(NSt);
263 for (uint i = 0; i < NSt; i++)
264 fvYoff[i] = 0.;
265 fvZoff.resize(NSt);
266 for (uint i = 0; i < NSt; i++)
267 fvZoff[i] = 0.;
268
269 if (fCalParFileName.IsNull()) return kTRUE;
270
272 TFile* oldFile = gFile;
273 TDirectory* oldDir = gDirectory;
274
275 fCalParFile = new TFile(fCalParFileName, "");
276 if (NULL == fCalParFile) {
277 LOG(error) << "CbmTofExtendTracks::LoadCalParameter: "
278 << "file " << fCalParFileName << " does not exist!";
279 return kTRUE;
280 }
281
282 LOG(info) << "CbmTofExtendTracks::LoadCalParameter: "
283 << " read from file " << fCalParFileName;
284
285 TH1D* fht = (TH1D*) gDirectory->FindObjectAny(Form("hExt_Toff"));
286 TH1D* fhx = (TH1D*) gDirectory->FindObjectAny(Form("hExt_Xoff"));
287 TH1D* fhy = (TH1D*) gDirectory->FindObjectAny(Form("hExt_Yoff"));
288 TH1D* fhz = (TH1D*) gDirectory->FindObjectAny(Form("hExt_Zoff"));
289
290 TH1D* fhts = (TH1D*) gDirectory->FindObjectAny(Form("hExt_Tsig"));
291 TH1D* fhxs = (TH1D*) gDirectory->FindObjectAny(Form("hExt_Xsig"));
292 TH1D* fhys = (TH1D*) gDirectory->FindObjectAny(Form("hExt_Ysig"));
293 TH1D* fhzs = (TH1D*) gDirectory->FindObjectAny(Form("hExt_Zsig"));
294
296 gFile = oldFile;
297 gDirectory = oldDir;
298
299 if (NULL != fht) {
300 fhExt_Toff = (TH1D*) fht->Clone();
301 for (UInt_t iSt = 0; iSt < NSt; iSt++) {
302 fvToff[iSt] = fhExt_Toff->GetBinContent(iSt + 1);
303 }
304 }
305 if (NULL != fhx) {
306 fhExt_Xoff = (TH1D*) fhx->Clone();
307 for (UInt_t iSt = 0; iSt < NSt; iSt++) {
308 fvXoff[iSt] = fhExt_Xoff->GetBinContent(iSt + 1);
309 }
310 }
311 if (NULL != fhy) {
312 fhExt_Yoff = (TH1D*) fhy->Clone();
313 for (UInt_t iSt = 0; iSt < NSt; iSt++) {
314 fvYoff[iSt] = fhExt_Yoff->GetBinContent(iSt + 1);
315 }
316 }
317 if (NULL != fhz) {
318 fhExt_Zoff = (TH1D*) fhz->Clone();
319 for (UInt_t iSt = 0; iSt < NSt; iSt++) {
320 fvZoff[iSt] = fhExt_Zoff->GetBinContent(iSt + 1);
321 }
322 }
323
324 if (NULL != fhts) {
325 fhExt_Tsig = (TH1D*) fhts->Clone();
326 for (UInt_t iSt = 0; iSt < NSt; iSt++) {
327 fvTsig[iSt] = fhExt_Tsig->GetBinContent(iSt + 1);
328 }
329 }
330 if (NULL != fhxs) {
331 fhExt_Xsig = (TH1D*) fhxs->Clone();
332 for (UInt_t iSt = 0; iSt < NSt; iSt++) {
333 fvXsig[iSt] = fhExt_Xsig->GetBinContent(iSt + 1);
334 }
335 }
336 if (NULL != fhys) {
337 fhExt_Ysig = (TH1D*) fhys->Clone();
338 for (UInt_t iSt = 0; iSt < NSt; iSt++) {
339 fvYsig[iSt] = fhExt_Ysig->GetBinContent(iSt + 1);
340 }
341 }
342 if (NULL != fhzs) {
343 fhExt_Zsig = (TH1D*) fhzs->Clone();
344 for (UInt_t iSt = 0; iSt < NSt; iSt++) {
345 fvZsig[iSt] = fhExt_Zsig->GetBinContent(iSt + 1);
346 }
347 }
348
349
350 fCalParFile->Close();
351
352 return kTRUE;
353}
354//-------------------------------------------------------------------------------------------------
355Bool_t CbmTofExtendTracks::InitParameters() { return kTRUE; }
356// ----- SetParContainers -------------------------------------------------
358{
359 /*
360 FairRunAna* ana = FairRunAna::Instance();
361 FairRuntimeDb* rtdb = ana->GetRuntimeDb();
362 */
363}
364// -------------------------------------------------------------------------
365
367{
368 if (fiEvent <= NDefSetup) return kTRUE;
369
370
371 while (fiCorMode > 0) {
372 Int_t iCorMode = fiCorMode % 10;
373 fiCorMode /= 10;
374 Int_t iLev = fiCorSrc / 10;
375 LOG(info) << "UpdateCalHist on level " << iLev << " from src " << fiCorSrc << " in mode " << iCorMode;
376 switch (iCorMode) {
377 case 0: { // T
378 TH2* hCorDT = NULL;
379 if (fiCorSrc % 10 == 0)
380 hCorDT = fhTrkStationDT[iLev];
381 else
382 hCorDT = fhTrkPullDT[iLev];
383 if (NULL != hCorDT) {
384 Double_t nx = hCorDT->GetNbinsX();
385 if (NULL == fhExt_Toff) {
386 fhExt_Toff = new TH1D("hExt_Toff", ";station;Toff (ns)", nx, 0, nx);
387 LOG(warn) << "Created " << fhExt_Toff->GetName();
388 }
389 if (NULL == fhExt_Tsig) {
390 fhExt_Tsig = new TH1D("hExt_Tsig", ";station;Tsig (ns)", nx, 0, nx);
391 LOG(warn) << "Created " << fhExt_Tsig->GetName();
392 }
393 for (Int_t ix = 0; ix < nx; ix++) {
394 Double_t dVal = fhExt_Toff->GetBinContent(ix + 1);
395 TH1D* hpy = hCorDT->ProjectionY(Form("%s_py%d", hCorDT->GetName(), ix), ix + 1, ix + 1, "");
396 Double_t dFMean = 0.;
397 LOG(warn) << "TOff Entries for station " << ix << ": " << hpy->GetEntries();
398 if (hpy->GetEntries() > 100) {
399 Int_t iBmax = hpy->GetMaximumBin();
400 TAxis* xaxis = hpy->GetXaxis();
401 Double_t dMean = xaxis->GetBinCenter(iBmax); //X-value of bin with maximal content
402 Double_t dLim = 1. * hpy->GetRMS();
403 //Double_t dLim = 5. * hpy->GetBinWidth(1);
404 if (dLim > 0.) {
405 //TF1 * f = new TF1 ("f","gaus");
406 //f->SetParLimits(2,dMean - dLim, dMean + dLim);
407 //f->SetParLimits(3,0.,dLim);
408 TFitResultPtr fRes =
409 // hpy->Fit(f,"B");
410 hpy->Fit("gaus", "SQM", "", dMean - dLim, dMean + dLim);
411 Int_t iFitStatus = fRes;
412 //LOG(warn)<<Form("TRes 0x%08x ",(UInt_t)iFitStatus);
413
414 if (iFitStatus != -1) {
415 dFMean = fRes->Parameter(1);
416 dVal -= dFMean;
417 fhExt_Tsig->SetBinContent(ix + 1, fRes->Parameter(2));
418 }
419 else
420 dVal -= dMean;
421 LOG(warn) << "Update hExt_Toff Ind " << ix << ": Old " << fhExt_Toff->GetBinContent(ix + 1)
422 << ", FitMean " << dFMean << " => " << dVal << ", width "
423 << fhExt_Tsig->GetBinContent(ix + 1);
424 }
425 }
426 else {
427 LOG(warn) << "Update hExt_Toff " << ix << ": insufficient counts: " << hpy->GetEntries();
428 }
429 fhExt_Toff->SetBinContent(ix + 1, dVal);
430 }
431 }
432 } break;
433 case 1: { // X
434 TH2* hCorDX = NULL;
435 if (fiCorSrc % 10 == 0)
436 hCorDX = fhTrkStationDX[iLev];
437 else
438 hCorDX = fhTrkPullDX[iLev];
439 if (NULL != hCorDX) {
440 Double_t nx = hCorDX->GetNbinsX();
441 if (NULL == fhExt_Xoff) {
442 fhExt_Xoff = new TH1D("hExt_Xoff", ";station;Xoff (cm)", nx, 0, nx);
443 LOG(warn) << "Created " << fhExt_Xoff->GetName();
444 }
445 if (NULL == fhExt_Xsig) {
446 fhExt_Xsig = new TH1D("hExt_Xsig", ";station;Xsig (cm)", nx, 0, nx);
447 LOG(warn) << "Created " << fhExt_Xsig->GetName();
448 }
449 for (Int_t ix = 0; ix < nx; ix++) {
450 Double_t dVal = fhExt_Xoff->GetBinContent(ix + 1);
451 TH1D* hpy = hCorDX->ProjectionY(Form("%s_py%d", hCorDX->GetName(), ix), ix + 1, ix + 1, "");
452 Double_t dFMean = 0.;
453 LOG(warn) << "XOff Entries for station " << ix << ": " << hpy->GetEntries();
454 if (hpy->GetEntries() > 100) {
455 Int_t iBmax = hpy->GetMaximumBin();
456 TAxis* xaxis = hpy->GetXaxis();
457 Double_t dMean = xaxis->GetBinCenter(iBmax); //X-value of bin with maximal content
458 //Double_t dLim = 1. * hpy->GetRMS();
459 Double_t dLim = 5. * hpy->GetBinWidth(1);
460 if (dLim > 0.) {
461 //TF1 * f = new TF1 ("f","gaus");
462 //f->SetParLimits(2,dMean - dLim, dMean + dLim);
463 //f->SetParLimits(3,0.,dLim);
464 TFitResultPtr fRes =
465 // hpy->Fit(f,"B");
466 hpy->Fit("gaus", "SQM", "", dMean - dLim, dMean + dLim);
467 Int_t iFitStatus = fRes;
468 //LOG(warn)<<Form("XRes 0x%08x ",(UInt_t)iFitStatus);
469 if (iFitStatus != -1) { // check validity of fit
470 dFMean = fRes->Parameter(1);
471 dVal -= dFMean;
472 fhExt_Xsig->SetBinContent(ix + 1, fRes->Parameter(2));
473 }
474 else
475 dVal -= dMean;
476 LOG(warn) << "Update hExt_Xoff Ind " << ix << ": Old " << fhExt_Xoff->GetBinContent(ix + 1)
477 << ", FitMean " << dFMean << " => " << dVal << ", width "
478 << fhExt_Xsig->GetBinContent(ix + 1);
479 }
480 }
481 else {
482 LOG(warn) << "Update hExt_Xoff " << ix << ": insufficient counts: " << hpy->GetEntries();
483 }
484 fhExt_Xoff->SetBinContent(ix + 1, dVal);
485 }
486 }
487 } break;
488 case 2: { // Y
489 TH2* hCorDY = NULL;
490 if (fiCorSrc % 10 == 0)
491 hCorDY = fhTrkStationDY[iLev];
492 else
493 hCorDY = fhTrkPullDY[iLev];
494 if (NULL != hCorDY) {
495 Double_t nx = hCorDY->GetNbinsX();
496 if (NULL == fhExt_Yoff) {
497 fhExt_Yoff = new TH1D("hExt_Yoff", ";station;Yoff (cm)", nx, 0, nx);
498 LOG(warn) << "Created " << fhExt_Yoff->GetName();
499 }
500 if (NULL == fhExt_Ysig) {
501 fhExt_Ysig = new TH1D("hExt_Ysig", ";station;Ysig (cm)", nx, 0, nx);
502 LOG(warn) << "Created " << fhExt_Ysig->GetName();
503 }
504 for (Int_t ix = 0; ix < nx; ix++) {
505 Double_t dVal = fhExt_Yoff->GetBinContent(ix + 1);
506 TH1D* hpy = hCorDY->ProjectionY(Form("%s_py%d", hCorDY->GetName(), ix), ix + 1, ix + 1, "");
507 Double_t dFMean = 0.;
508 LOG(warn) << "YOff Entries for station " << ix << ": " << hpy->GetEntries();
509 if (hpy->GetEntries() > 100) {
510 Int_t iBmax = hpy->GetMaximumBin();
511 TAxis* xaxis = hpy->GetXaxis();
512 Double_t dMean = xaxis->GetBinCenter(iBmax); //X-value of bin with maximal content
513 //Double_t dLim = 1. * hpy->GetRMS();
514 Double_t dLim = 5. * hpy->GetBinWidth(1);
515 if (dLim > 0.) {
516 //TF1 * f = new TF1 ("f","gaus");
517 //f->SetParLimits(2,dMean - dLim, dMean + dLim);
518 //f->SetParLimits(3,0.,dLim);
519 // hpy->Draw("");
520 LOG(warn) << "Fit gaus with " << dMean << ", " << dLim;
521 TFitResultPtr fRes =
522 // hpy->Fit(f,"B");
523 hpy->Fit("gaus", "SQM", "", dMean - dLim, dMean + dLim);
524 Int_t iFitStatus = fRes;
525 //LOG(warn)<<Form("YRes 0x%08x ",iFitStatus);
526 if (iFitStatus != -1) { // check validity of fit
527 dFMean = fRes->Parameter(1);
528 dVal -= dFMean;
529 fhExt_Ysig->SetBinContent(ix + 1, fRes->Parameter(2));
530 }
531 else
532 dVal -= dMean;
533
534 LOG(warn) << "Update hExt_Yoff Ind " << ix << ": Old " << fhExt_Yoff->GetBinContent(ix + 1)
535 << ", FitMean " << dFMean << " => " << dVal << ", width "
536 << fhExt_Ysig->GetBinContent(ix + 1);
537 }
538 }
539 else {
540 LOG(warn) << "Update hExt_Yoff " << ix << ": insufficient counts: " << hpy->GetEntries();
541 }
542 fhExt_Yoff->SetBinContent(ix + 1, dVal);
543 }
544 }
545 } break;
546 default: LOG(info) << "Correction Mode not implemented";
547 } // switch end
548 }
549 return kTRUE;
550}
552{
553 if (fiCorMode < 0) return kTRUE;
554 // Write calibration histogramms to the file
555 TFile* oldFile = gFile;
556 TDirectory* oldDir = gDirectory;
557 TFile* fHist = new TFile(fCalOutFileName, "RECREATE");
558 fHist->cd();
559
560 if (NULL != fhExt_Toff) fhExt_Toff->Write();
561 if (NULL != fhExt_Xoff) fhExt_Xoff->Write();
562 if (NULL != fhExt_Yoff) fhExt_Yoff->Write();
563 if (NULL != fhExt_Zoff) fhExt_Zoff->Write();
564
565 gFile = oldFile;
566 gDirectory = oldDir;
567
568 fHist->Close();
569 return kTRUE;
570}
571
572
573// ----- Public method Exec --------------------------------------------
574void CbmTofExtendTracks::Exec(Option_t* opt)
575{
576
577 if (!fEventsColl) {
578 // fTofHitArray = (TClonesArray*)fTofHitArrayIn->Clone();
579 LOG(fatal) << "Analysis needs EventsColl ";
580 // ExecExtend(opt);
581 }
582 else {
583 LOG(debug) << "ExtExec TS " << fiTS << " with " << fEventsColl->GetEntriesFast() << " evts";
584 fiTS++;
585 for (Int_t iEvent = 0; iEvent < fEventsColl->GetEntriesFast(); iEvent++) {
586 CbmEvent* tEvent = dynamic_cast<CbmEvent*>(fEventsColl->At(iEvent));
587 LOG(debug) << "Process TS event " << iEvent << " with " << tEvent->GetNofData(ECbmDataType::kBmonHit) << " T0, "
588 << tEvent->GetNofData(ECbmDataType::kTofHit) << " TOF, " << tEvent->GetNofData(ECbmDataType::kStsHit)
589 << " STS, " << tEvent->GetNofData(ECbmDataType::kMuchPixelHit) << " MUCH hits";
590
591 ExecExtend(opt, tEvent);
592 }
593 }
594}
595
596void CbmTofExtendTracks::ExecExtend(Option_t* /*opt*/, CbmEvent* tEvent)
597{
598 fiEvent++;
600
601 UInt_t iNbTofStations = 0;
602 UInt_t iNbStsStations = 0;
603 UInt_t iNbMuchStations = 0;
604 UInt_t iNbRichStations = 0;
605 const Double_t STATION_TOF_ZWIDTH = 10.;
606 const Double_t STATION_STS_ZWIDTH = 1.5;
607 const Double_t STATION_MUCH_ZWIDTH = 1.;
608 const Double_t STATION_RICH_ZWIDTH = 3.;
609 fvTofHitIndex.clear();
610 fvStsHitIndex.clear();
611 fvMuchHitIndex.clear();
612 fvRichHitIndex.clear();
613
614 LOG(debug) << "Ev " << fiEvent << " has hits "
615 << " Tof: " << tEvent->GetNofData(ECbmDataType::kTofHit)
616 << " Sts: " << tEvent->GetNofData(ECbmDataType::kStsHit)
617 << " Much: " << tEvent->GetNofData(ECbmDataType::kMuchPixelHit)
618 << " Rich: " << tEvent->GetNofData(ECbmDataType::kRichHit);
619
620 // cleanup
621 /*
622 if (NULL != fvAllHitPointer )
623 for(UInt_t iSt=0; iSt<fvAllHitPointer.size(); iSt++) {
624 for (UInt_t iHit=0; iHit<fvAllHitPointer[iSt].size(); iHit++)
625 delete(fvAllHitPointer[iSt][iHit]);
626 fvAllHitPointer[iSt].clear();
627 }
628 */
629 fvAllHitPointer.clear();
630 UInt_t iNbAllStations = fMapStationZ.size();
631 fvAllHitPointer.resize(iNbAllStations);
632
633 if (tEvent->GetNofData(ECbmDataType::kTofHit) > 0) {
634 // init to 1 tracking station
635 iNbTofStations = 1;
636 fvTofStationZ.resize(iNbTofStations);
637 fvTofHitIndex.resize(iNbTofStations);
638 for (size_t iHit = 0; iHit < tEvent->GetNofData(ECbmDataType::kTofHit); iHit++) {
639 Int_t iHitIndex = static_cast<Int_t>(tEvent->GetIndex(ECbmDataType::kTofHit, iHit));
640 CbmTofHit* tHit = dynamic_cast<CbmTofHit*>(fTofHitArrayIn->At(iHitIndex));
641 LOG(debug) << Form("Inspect Ev %d, TofHit %lu, Ind %d at %6.1f in %lu (%u) stations", fiEvent, iHit, iHitIndex,
642 tHit->GetZ(), fvTofStationZ.size(), iNbAllStations);
643
644 Int_t iStZ = (Int_t)(tHit->GetZ() / dStDZ);
645 itMapStationZ = fMapStationZ.find(iStZ);
646
647 UInt_t iSt = 0;
648 for (; iSt < fvTofStationZ.size(); iSt++) {
649 if (iSt == 0 && fvTofHitIndex[iSt].size() == 1) {
650 fvTofStationZ[iSt] = tHit->GetZ();
651 fvTofHitIndex[iSt].resize(1);
652 fvTofHitIndex[iSt][0] = iHitIndex;
653 LOG(debug) << Form("Ev %d, init TofSt %d, Mul %lu at z %8.1f from Ind %d", fiEvent, iSt,
654 fvTofHitIndex[iSt].size(), fvTofStationZ[iSt], iHitIndex);
655 break;
656 }
657 else {
658 if (TMath::Abs(fvTofStationZ[iSt] - tHit->GetZ()) < STATION_TOF_ZWIDTH) {
659 // update average z position of station
660 fvTofStationZ[iSt] =
661 (fvTofStationZ[iSt] * fvTofHitIndex[iSt].size() + tHit->GetZ()) / (fvTofHitIndex[iSt].size() + 1);
662 fvTofHitIndex[iSt].resize(fvTofHitIndex[iSt].size() + 1);
663 fvTofHitIndex[iSt][fvTofHitIndex[iSt].size() - 1] = iHitIndex;
664 LOG(debug) << Form("Ev %d, upd TofSt %d, Mul %lu at z %8.1f from Ind %d", fiEvent, iSt,
665 fvTofHitIndex[iSt].size(), fvTofStationZ[iSt], iHitIndex);
666 break;
667 }
668 }
669 } // Station loop end
670
671 if (iSt == iNbTofStations) { // add new station
672 iNbTofStations++;
673 fvTofStationZ.resize(iNbTofStations);
674 fvTofHitIndex.resize(iNbTofStations);
675 fvTofStationZ[iSt] = tHit->GetZ();
676 fvTofHitIndex[iSt].resize(1);
677 fvTofHitIndex[iSt][0] = iHitIndex;
678 LOG(debug) << Form("Ev %d, add TofSt %d (%d), Mul %lu at z %8.1f from Ind %d", fiEvent, iSt, iNbTofStations,
679 fvTofHitIndex[iSt].size(), fvTofStationZ[iSt], iHitIndex);
680 }
681 if (fiEvent < NDefSetup) {
682 if (itMapStationZ == fMapStationZ.end()) {
683 LOG(info) << "Insert new tracking station " << Int_t(ECbmModuleId::kTof) * 100 + iSt << " at z=" << iStZ;
684 fMapStationZ[iStZ] = Int_t(ECbmModuleId::kTof) * 100 + iSt;
685 itMapStationZ = fMapStationZ.begin();
686 Int_t iStId = Int_t(ECbmModuleId::kTof) * 100;
687 for (; itMapStationZ != fMapStationZ.end(); ++itMapStationZ) {
688 Int_t iSysId = itMapStationZ->second / 100;
689 if (iSysId == Int_t(ECbmModuleId::kTof)) {
690 itMapStationZ->second = iStId++;
691 }
692 LOG(info) << "MapStationZ: " << itMapStationZ->first << " " << itMapStationZ->second;
693 }
694 }
695 }
696 else { // Define Setup end
697 // sort hits into stations
698 if (itMapStationZ != fMapStationZ.end()) {
699 Int_t iAllSt = (itMapStationZ->second) % 100;
700 /* TOF recalibration disabled
701 pHit->SetX(pHit->GetX() + fvXoff[iAllSt]);
702 pHit->SetY(pHit->GetY() + fvYoff[iAllSt]);
703 pHit->SetZ(pHit->GetZ() + fvZoff[iAllSt]);
704 pHit->SetTime(pHit->GetTime() + fvToff[iAllSt]);
705 */
706 //CbmTofHit* pHit = new CbmTofHit(*tHit); // copy construction
707 // apply recalibration if necessary
708 fvAllHitPointer[iAllSt].push_back(dynamic_cast<CbmPixelHit*>(tHit));
709 //LOG(info)<<"TpHit: "<<pHit->ToString();
710 }
711 }
712 } // TofHit loop end
713 } // TofHit condition end
714
715
716 if (tEvent->GetNofData(ECbmDataType::kStsHit) > 0) {
717 // init to 1 tracking station
718 iNbStsStations = 1;
719 fvStsStationZ.resize(iNbStsStations);
720 fvStsHitIndex.resize(iNbStsStations);
721
722 for (size_t iHit = 0; iHit < tEvent->GetNofData(ECbmDataType::kStsHit); iHit++) {
723 Int_t iHitIndex = static_cast<Int_t>(tEvent->GetIndex(ECbmDataType::kStsHit, iHit));
724 CbmPixelHit* tHit = dynamic_cast<CbmPixelHit*>(fStsHitArrayIn->At(iHitIndex));
725
726 Int_t iStZ = (Int_t)(tHit->GetZ() / dStDZ);
727 itMapStationZ = fMapStationZ.find(iStZ);
728
729 UInt_t iSt = 0;
730 for (; iSt < fvStsStationZ.size(); iSt++) {
731 if (iSt == 0 && fvStsHitIndex[iSt].size() == 1) {
732 fvStsStationZ[iSt] = tHit->GetZ();
733 fvStsHitIndex[iSt].resize(1);
734 fvStsHitIndex[iSt][0] = iHitIndex;
735 LOG(debug) << Form("Ev %d, init StsSt %d at z %8.1f from Ind %d", fiEvent, iSt, fvStsStationZ[iSt],
736 iHitIndex);
737 break;
738 }
739 else {
740 if (TMath::Abs(fvStsStationZ[iSt] - tHit->GetZ()) < STATION_STS_ZWIDTH) {
741 // update average z position of station
742 fvStsStationZ[iSt] =
743 (fvStsStationZ[iSt] * fvStsHitIndex[iSt].size() + tHit->GetZ()) / (fvStsHitIndex[iSt].size() + 1);
744 fvStsHitIndex[iSt].resize(fvStsHitIndex[iSt].size() + 1);
745 fvStsHitIndex[iSt][fvStsHitIndex[iSt].size() - 1] = iHitIndex;
746 LOG(debug) << Form("Ev %d, upd StsSt %d at z %8.1f from Ind %d", fiEvent, iSt, fvStsStationZ[iSt],
747 iHitIndex);
748 break;
749 }
750 }
751 } // Station loop end
752 if (iSt == iNbStsStations) { // add new station
753 iNbStsStations++;
754 fvStsStationZ.resize(iNbStsStations);
755 fvStsHitIndex.resize(iNbStsStations);
756 fvStsStationZ[iSt] = tHit->GetZ();
757 fvStsHitIndex[iSt].resize(fvStsHitIndex[iSt].size() + 1);
758 fvStsHitIndex[iSt][fvStsHitIndex[iSt].size() - 1] = iHitIndex;
759 LOG(debug) << Form("Ev %d, add StsSt %d (%d) at z %8.1f from Ind %d", fiEvent, iSt, iNbStsStations,
760 fvStsStationZ[iSt], iHitIndex);
761 }
762 if (fiEvent < NDefSetup) {
763 if (itMapStationZ == fMapStationZ.end()) {
764 LOG(info) << "Insert new tracking station " << Int_t(ECbmModuleId::kSts) * 100 + iSt << " at z=" << iStZ;
765 fMapStationZ[iStZ] = Int_t(ECbmModuleId::kSts) * 100 + iSt;
766 itMapStationZ = fMapStationZ.begin();
767 Int_t iStId = Int_t(ECbmModuleId::kSts) * 100;
768 for (; itMapStationZ != fMapStationZ.end(); ++itMapStationZ) {
769 Int_t iSysId = itMapStationZ->second / 100;
770 if (iSysId == Int_t(ECbmModuleId::kSts)) {
771 itMapStationZ->second = iStId++;
772 }
773 LOG(info) << "MapStationZ: " << itMapStationZ->first << " " << itMapStationZ->second;
774 }
775 }
776 }
777 else { // Define Setup end
778 // sort hits into stations
779 if (itMapStationZ != fMapStationZ.end()) {
780 Int_t iAllSt = (itMapStationZ->second) % 100;
781 UInt_t iH = fvAllHitPointer[iAllSt].size();
782 fvAllHitPointer[iAllSt].push_back(dynamic_cast<CbmPixelHit*>(tHit));
783 CbmPixelHit* pHit = fvAllHitPointer[iAllSt][iH];
784 pHit->SetX(pHit->GetX() + fvXoff[iAllSt]);
785 pHit->SetY(pHit->GetY() + fvYoff[iAllSt]);
786 pHit->SetZ(pHit->GetZ() + fvZoff[iAllSt]);
787 pHit->SetTime(pHit->GetTime() + fvToff[iAllSt]);
788 }
789 }
790 } // Sts Hit loop end
791 } //Sts condition end
792
793 if (tEvent->GetNofData(ECbmDataType::kMuchPixelHit) > 0) {
794 // init to 1 tracking station
795 iNbMuchStations = 1;
796 fvMuchStationZ.resize(iNbMuchStations);
797 fvMuchHitIndex.resize(iNbMuchStations);
798 for (size_t iHit = 0; iHit < tEvent->GetNofData(ECbmDataType::kMuchPixelHit); iHit++) {
799 Int_t iHitIndex = static_cast<Int_t>(tEvent->GetIndex(ECbmDataType::kMuchPixelHit, iHit));
800 CbmPixelHit* tHit = dynamic_cast<CbmPixelHit*>(fMuchHitArrayIn->At(iHitIndex));
801
802 Int_t iStZ = (Int_t)(tHit->GetZ() / dStDZ);
803 itMapStationZ = fMapStationZ.find(iStZ);
804
805 UInt_t iSt = 0;
806 for (; iSt < fvMuchStationZ.size(); iSt++) {
807 if (iSt == 0 && fvMuchHitIndex[iSt].size() == 1) {
808 fvMuchStationZ[iSt] = tHit->GetZ();
809 fvMuchHitIndex[iSt].resize(1);
810 fvMuchHitIndex[iSt][0] = iHitIndex;
811 LOG(info) << Form("Ev %d, init MuchSt %d at z %8.1f from Ind %d", fiEvent, iSt, fvMuchStationZ[iSt],
812 iHitIndex);
813 break;
814 }
815 else {
816 if (TMath::Abs(fvMuchStationZ[iSt] - tHit->GetZ()) < STATION_MUCH_ZWIDTH) {
817 // update average z position of station
818 fvMuchStationZ[iSt] =
819 (fvMuchStationZ[iSt] * fvMuchHitIndex[iSt].size() + tHit->GetZ()) / (fvMuchHitIndex[iSt].size() + 1);
820 fvMuchHitIndex[iSt].resize(fvMuchHitIndex[iSt].size() + 1);
821 fvMuchHitIndex[iSt][fvMuchHitIndex[iSt].size() - 1] = iHitIndex;
822 LOG(debug) << Form("Ev %d, upd MuchSt %d at z %8.1f from Ind %d", fiEvent, iSt, fvMuchStationZ[iSt],
823 iHitIndex);
824 break;
825 }
826 }
827 } // Station loop end
828 if (iSt == iNbMuchStations) { // add new station
829 iNbMuchStations++;
830 fvMuchStationZ.resize(iNbMuchStations);
831 fvMuchHitIndex.resize(iNbMuchStations);
832 fvMuchStationZ[iSt] = tHit->GetZ();
833 fvMuchHitIndex[iSt].resize(fvMuchHitIndex[iSt].size() + 1);
834 fvMuchHitIndex[iSt][fvMuchHitIndex[iSt].size() - 1] = iHitIndex;
835 LOG(debug) << Form("Ev %d, add MuchSt %d (%d) at z %8.1f from Ind %d", fiEvent, iSt, iNbMuchStations,
836 fvMuchStationZ[iSt], iHitIndex);
837 }
838
839 if (fiEvent < NDefSetup) {
840 if (itMapStationZ == fMapStationZ.end()) {
841 LOG(info) << "Insert new tracking station " << Int_t(ECbmModuleId::kMuch) * 100 + iSt << " at z=" << iStZ;
842 fMapStationZ[iStZ] = Int_t(ECbmModuleId::kMuch) * 100 + iSt;
843 itMapStationZ = fMapStationZ.begin();
844 Int_t iStId = Int_t(ECbmModuleId::kMuch) * 100;
845 for (; itMapStationZ != fMapStationZ.end(); ++itMapStationZ) {
846 Int_t iSysId = itMapStationZ->second / 100;
847 if (iSysId == Int_t(ECbmModuleId::kMuch)) {
848 itMapStationZ->second = iStId++;
849 }
850 LOG(info) << "MapStationZ: " << itMapStationZ->first << " " << itMapStationZ->second;
851 }
852 }
853 }
854 else { // Define Setup end
855 // sort hits into stations
856 if (itMapStationZ != fMapStationZ.end()) {
857 Int_t iAllSt = (itMapStationZ->second) % 100;
858 UInt_t iH = fvAllHitPointer[iAllSt].size();
859 fvAllHitPointer[iAllSt].push_back(dynamic_cast<CbmPixelHit*>(tHit));
860 CbmPixelHit* pHit = fvAllHitPointer[iAllSt][iH];
861 Double_t Xin = pHit->GetX();
862 pHit->SetX(pHit->GetX() + fvXoff[iAllSt]);
863 Double_t Xout = pHit->GetX();
864 // LOG(info) << "CalMuchX in "<<Xin<<", out "<<Xout<<", diff: "<<Xout-Xin;
865 if (fvXoff[iAllSt] != 0.) assert(Xin != Xout);
866
867 pHit->SetY(pHit->GetY() + fvYoff[iAllSt]);
868 pHit->SetZ(pHit->GetZ() + fvZoff[iAllSt]);
869 pHit->SetTime(pHit->GetTime() + fvToff[iAllSt]);
870 //assign errors
871 const TVector3 MuchHitError = {0.5, 0.5, 0.5};
872 pHit->SetPositionError(MuchHitError);
873
874 assert(iAllSt == fMapStationZ[(Int_t)(pHit->GetZ())] % 100);
875
876 LOG(info) << Form("Proc ev %d, MuchIn St %d, H %d, MpHit: ", fiEvent, iAllSt, iH) << tHit->ToString();
877
878 //delete(pHit);
879 }
880 else {
881 LOG(warn) << "Undefined station for Much ";
882 }
883 }
884 } // MuchHit loop end
885 } // Much Hit condition end
886
887 if (tEvent->GetNofData(ECbmDataType::kRichHit) > 0) {
888 // init to 1 tracking station
889 iNbRichStations = 1;
890 fvRichStationZ.resize(iNbRichStations);
891 fvRichHitIndex.resize(iNbRichStations);
892 for (size_t iHit = 0; iHit < tEvent->GetNofData(ECbmDataType::kRichHit); iHit++) {
893 Int_t iHitIndex = static_cast<Int_t>(tEvent->GetIndex(ECbmDataType::kRichHit, iHit));
894 CbmPixelHit* tHit = dynamic_cast<CbmPixelHit*>(fRichHitArrayIn->At(iHitIndex));
895
896 Int_t iStZ = (Int_t)(tHit->GetZ() / dStDZ);
897 itMapStationZ = fMapStationZ.find(iStZ);
898
899 UInt_t iSt = 0;
900 for (; iSt < fvRichStationZ.size(); iSt++) {
901 if (iSt == 0 && fvRichHitIndex[iSt].size() == 1) {
902 fvRichStationZ[iSt] = tHit->GetZ();
903 fvRichHitIndex[iSt].resize(1);
904 fvRichHitIndex[iSt][0] = iHitIndex;
905 LOG(debug) << Form("Ev %d, init RichSt %d at z %8.1f from Ind %d", fiEvent, iSt, fvRichStationZ[iSt],
906 iHitIndex);
907 break;
908 }
909 else {
910 if (TMath::Abs(fvRichStationZ[iSt] - tHit->GetZ()) < STATION_RICH_ZWIDTH) {
911 // update average z position of station
912 fvRichStationZ[iSt] =
913 (fvRichStationZ[iSt] * fvRichHitIndex[iSt].size() + tHit->GetZ()) / (fvRichHitIndex[iSt].size() + 1);
914 fvRichHitIndex[iSt].resize(fvRichHitIndex[iSt].size() + 1);
915 fvRichHitIndex[iSt][fvRichHitIndex[iSt].size() - 1] = iHitIndex;
916 LOG(debug) << Form("Ev %d, upd RichSt %d at z %8.1f from Ind %d", fiEvent, iSt, fvRichStationZ[iSt],
917 iHitIndex);
918 break;
919 }
920 }
921 } // Station loop end
922 if (iSt == iNbRichStations) { // add new station
923 iNbRichStations++;
924 fvRichStationZ.resize(iNbRichStations);
925 fvRichHitIndex.resize(iNbRichStations);
926 fvRichStationZ[iSt] = tHit->GetZ();
927 fvRichHitIndex[iSt].resize(fvRichHitIndex[iSt].size() + 1);
928 fvRichHitIndex[iSt][fvRichHitIndex[iSt].size() - 1] = iHitIndex;
929 LOG(debug) << Form("Ev %d, add RichSt %d (%d) at z %8.1f from Ind %d", fiEvent, iSt, iNbRichStations,
930 fvRichStationZ[iSt], iHitIndex);
931 }
932
933 if (fiEvent < NDefSetup) {
934 if (itMapStationZ == fMapStationZ.end()) {
935 LOG(info) << "Insert new tracking station " << Int_t(ECbmModuleId::kRich) * 100 + iSt << " at z=" << iStZ;
936 fMapStationZ[iStZ] = Int_t(ECbmModuleId::kRich) * 100 + iSt;
937 itMapStationZ = fMapStationZ.begin();
938 Int_t iStId = Int_t(ECbmModuleId::kRich) * 100;
939 for (; itMapStationZ != fMapStationZ.end(); ++itMapStationZ) {
940 Int_t iSysId = itMapStationZ->second / 100;
941 if (iSysId == Int_t(ECbmModuleId::kRich)) {
942 itMapStationZ->second = iStId++;
943 }
944 LOG(info) << "MapStationZ: " << itMapStationZ->first << " " << itMapStationZ->second;
945 }
946 }
947 }
948 else { // Define Setup end
949 // sort hits into stations
950 if (itMapStationZ != fMapStationZ.end()) {
951 Int_t iAllSt = (itMapStationZ->second) % 100;
952 //CbmPixelHit* pHit = new CbmPixelHit(*tHit); // copy construction
953 // apply recalibration if necessary
954 UInt_t iH = fvAllHitPointer[iAllSt].size();
955 fvAllHitPointer[iAllSt].push_back(dynamic_cast<CbmPixelHit*>(tHit));
956 CbmPixelHit* pHit = fvAllHitPointer[iAllSt][iH];
957 pHit->SetX(pHit->GetX() + fvXoff[iAllSt]);
958 pHit->SetY(pHit->GetY() + fvYoff[iAllSt]);
959 pHit->SetZ(pHit->GetZ() + fvZoff[iAllSt]);
960 pHit->SetTime(pHit->GetTime() + fvToff[iAllSt]);
961
962 //LOG(info)<<"RpHit: "<<pHit->ToString();
963 //delete(pHit);
964 }
965 }
966 } //RichHit loop end
967 } // Rich condition end
968 /*
969 LOG(info) << Form(
970 "Ev %d, found %d Tof, %d Sts, %d Much, %d Rich - stations, %d tracks",
971 fiEvent,
972 iNbTofStations,
973 iNbStsStations,
974 iNbMuchStations,
975 iNbRichStations,
976 tEvent->GetNofData(ECbmDataType::kTofTracklet));
977
978 LOG(info) << "Station Multipliplicities:";
979 for(UInt_t iSt=0; iSt<fvAllHitPointer.size();iSt++)
980 LOG(info)<<" "<<iSt<<": "<<fvAllHitPointer[iSt].size();
981 */
982
983 // loop on tof tracklet
984 size_t iNTrks = tEvent->GetNofData(ECbmDataType::kTofTracklet);
985
986 fvTrkCalHits.clear();
987 fvTrkCalHits.resize(iNTrks);
988 fvTrkPar.clear();
989 for (size_t iTr = 0; iTr < iNTrks; iTr++) {
990 Int_t iTrkIndex = static_cast<Int_t>(tEvent->GetIndex(ECbmDataType::kTofTracklet, iTr));
991 CbmTofTracklet* pTrk = dynamic_cast<CbmTofTracklet*>(fTofTrackArrayIn->At(iTrkIndex));
992 //CbmTofTrackletParam* pTrkPar=new CbmTofTrackletParam(*pTrk->GetTrackParameter());
993 //fvTrkPar.push_back(pTrkPar);
994 fvTrkPar.push_back((CbmTofTrackletParam*) (pTrk->GetTrackParameter()));
995
996 Int_t iNHits = pTrk->GetNofHits();
997 fvTrkCalHits[iTr].resize(0);
998 for (Int_t iHit = 0; iHit < iNHits; iHit++) {
999 Int_t iHitInd = pTrk->GetTofHitIndex(iHit);
1000 CbmTofHit* pHitIn = dynamic_cast<CbmTofHit*>(fTofHitArrayIn->At(iHitInd));
1001 if (NULL == pHitIn) {
1002 LOG(warn) << Form("Hit %d not found at index %d ", iHit, iHitInd);
1003 continue;
1004 }
1005 //CbmTofHit* pHit = new CbmTofHit(*pHitIn); // copy construction
1006 // apply recalibration if necessary
1007 fvTrkCalHits[iTr].push_back(dynamic_cast<CbmPixelHit*>(pHitIn));
1008
1009 LOG(debug) << Form("Added PixHit %d, ind %d: x %6.3f, y %6.3f, z %6.3f, t %9.2f "
1010 ", dx %6.3f, dy %6.3f, dz %6.3f ",
1011 iHit, pTrk->GetTofHitIndex(iHit), (fvTrkCalHits[iTr][iHit])->GetX(),
1012 (fvTrkCalHits[iTr][iHit])->GetY(), (fvTrkCalHits[iTr][iHit])->GetZ(),
1013 (fvTrkCalHits[iTr][iHit])->GetTime(), (fvTrkCalHits[iTr][iHit])->GetDx(),
1014 (fvTrkCalHits[iTr][iHit])->GetDy(), (fvTrkCalHits[iTr][iHit])->GetDz());
1015 }
1016 Line3Dfit(fvTrkCalHits[iTr], fvTrkPar[iTr]);
1017 fvTrkPar[iTr]->SetTt(pTrk->GetTt());
1018 fvTrkPar[iTr]->SetT(pTrk->GetT0());
1019
1020 // compare to input
1021 LOG(debug) << "CompareTrk " << iTr
1022 << Form(": DTx %10.8f, DTy %10.8f", pTrk->GetTrackTx() - fvTrkPar[iTr]->GetTx(),
1023 pTrk->GetTrackTy() - fvTrkPar[iTr]->GetTy());
1024 } // Tracklet loop end
1025 if (fiEvent > NDefSetup) {
1026 Int_t iAddStations = fiAddStations;
1027 while (iAddStations > 0) {
1028 Int_t iAddNextStation = iAddStations % 100;
1029 iAddStations /= 100;
1030 TrkAddStation(iAddNextStation);
1031 }
1032
1033 //FindVertex();
1034 FillHistograms(tEvent);
1035 //assert(fiEvent<1);
1036 }
1037}
1038// -------------------------------------------------------------------------
1039
1040
1041// ----- Public method Finish ------------------------------------------
1043{
1044
1046
1047 // save everything in output file
1048
1049 FairRootManager* ioman = FairRootManager::Instance();
1050 TList* tList = gROOT->GetList();
1051
1052 ((FairRootFileSink*) ioman->GetSink())->GetRootFile()->cd();
1053
1054 TIter next(tList);
1055 TObject* obj;
1056 while ((obj = (TObject*) next())) {
1057 if (obj->InheritsFrom(TH1::Class())) obj->Write();
1058 }
1059
1060 WriteHistos();
1061
1062 LOG(info) << Form(" CbmTofExtendTracks::Finished ");
1063}
1064// -------------------------------------------------------------------------
1065
1067{
1068 if (fiEvent < NDefSetup) return;
1069 TDirectory* oldir = gDirectory; // <= To prevent histos from being sucked in by the param file of the TRootManager!
1070 gROOT->cd(); // <= To prevent histos from being sucked in by the param file of the TRootManager !
1071
1072 // define histos here
1073 // Correlation with Tof Tracklets
1074 //
1076 new TH2F("hMulCorTrkTof", Form("Multiplicity correlation ; N_{Track}; N_{TofHits}"), 50, 0, 50, 150, 0, 150);
1077
1079 new TH2F("hMulCorTrkSts", Form("Multiplicity correlation ; N_{Track}; N_{StsHits}"), 50, 0, 50, 150, 0, 150);
1081 new TH2F("hMulCorTrkMuch", Form("Multiplicity correlation ; N_{Track}; N_{MuchHits}"), 50, 0, 50, 150, 0, 150);
1083 new TH2F("hMulCorTrkRich", Form("Multiplicity correlation ; N_{Track}; N_{RichHits}"), 50, 0, 50, 150, 0, 150);
1084 fhPosCorTrkTof = new TH2F("hPosCorTrkTof", Form("Tof position correlation ; #DeltaX (cm); #DeltaY (cm)"), 100, -10,
1085 10, 100, -10, 10);
1086 fhPosCorTrkSts = new TH2F("hPosCorTrkSts", Form("Sts position correlation ; #DeltaX (cm); #DeltaY (cm)"), 100, -10,
1087 10, 100, -10, 10);
1088 fhPosCorTrkMuch = new TH2F("hPosCorTrkMuch", Form("Much position correlation ; #DeltaX (cm); #DeltaY (cm)"), 100,
1089 -10, 10, 100, -10, 10);
1090 fhPosCorTrkRich = new TH2F("hPosCorTrkRich", Form("Rich position correlation ; #DeltaX (cm); #DeltaY (cm)"), 100,
1091 -10, 10, 100, -10, 10);
1092
1093
1094 // Correlation with extended tracks
1095 // Book Histos for all tracking stations found in fMap
1096 LOG(info) << "Book histos for the following stations: ";
1097 Int_t iStNum = 0;
1098 itMapStationZ = fMapStationZ.begin();
1099 for (; itMapStationZ != fMapStationZ.end(); ++itMapStationZ) {
1100 Int_t iSysId = itMapStationZ->second / 100;
1101 itMapStationZ->second = iSysId * 100 + iStNum++;
1102 LOG(info) << Form(" station %2d, z %3d, Id %4d ", iStNum, itMapStationZ->first, itMapStationZ->second);
1103 }
1104 UInt_t NSt = fMapStationZ.size();
1105
1106 Int_t NLev = 2;
1107 //if (fiAddStations > 0 ) NLev=2;
1108
1109 fhTrkStationDX.resize(NLev);
1110 fhTrkStationDY.resize(NLev);
1111 fhTrkStationDZ.resize(NLev);
1112 fhTrkStationDT.resize(NLev);
1113 fhTrkStationNHits.resize(NLev);
1114 fhTrkStationDXDY.resize(NLev);
1115 fhTrkPullDX.resize(NLev);
1116 fhTrkPullDY.resize(NLev);
1117 fhTrkPullDT.resize(NLev);
1118 fhExt_TrkSizVel.resize(NLev);
1119 fhExt_TrkSizChiSq.resize(NLev);
1120
1121 for (Int_t iLev = 0; iLev < NLev; iLev++) {
1122 fhExt_TrkSizVel[iLev] = new TH2F(Form("hExt_TrkSizVel%d", iLev), ";TrkSize;v (cm/ns)", 15, 0, 15, 100, 0., 50.);
1123 fhExt_TrkSizChiSq[iLev] = new TH2F(Form("hExt_TrkSizChiSq%d", iLev), ";TrkSize;#chi^2", 15, 0, 15, 50, 0., 5.);
1124
1125 fhTrkStationDX[iLev] = new TH2F(Form("hTrkStationDX%d", iLev), Form("TrkStationDX; StationNr ; #DeltaX (cm)"), NSt,
1126 0, NSt, 100, -10., 10.);
1127 fhTrkStationDY[iLev] = new TH2F(Form("hTrkStationDY%d", iLev), Form("TrkStationDY; StationNr ; #DeltaY (cm)"), NSt,
1128 0, NSt, 100, -10., 10.);
1129 fhTrkStationDZ[iLev] = new TH2F(Form("hTrkStationDZ%d", iLev), Form("TrkStationDZ; StationNr ; #DeltaZ (cm)"), NSt,
1130 0, NSt, 100, -20., 20.);
1131 fhTrkStationDT[iLev] = new TH2F(Form("hTrkStationDT%d", iLev), Form("TrkStationDT; StationNr ; #DeltaT (ns)"), NSt,
1132 0, NSt, 100, -50., 50.);
1133 fhTrkStationNHits[iLev] =
1134 new TH2F(Form("hTrkStationNHits%d", iLev), Form("TrkStationNHits; StationNr ; Number of Hits"), NSt, 0, NSt, 100,
1135 0., 100.);
1136
1137 fhTrkPullDX[iLev] =
1138 new TH2F(Form("hTrkPullDX%d", iLev), Form("TrkPullDX; StationNr ; #DeltaX (cm)"), NSt, 0, NSt, 100, -10., 10.);
1139 fhTrkPullDY[iLev] =
1140 new TH2F(Form("hTrkPullDY%d", iLev), Form("TrkPullDY; StationNr ; #DeltaY (cm)"), NSt, 0, NSt, 100, -10., 10.);
1141 fhTrkPullDT[iLev] =
1142 new TH2F(Form("hTrkPullDT%d", iLev), Form("TrkPullDT; StationNr ; #DeltaT (ns)"), NSt, 0, NSt, 100, -50., 50.);
1143
1144 fhTrkStationDXDY[iLev].resize(NSt);
1145 for (UInt_t iSt = 0; iSt < NSt; iSt++) {
1146 fhTrkStationDXDY[iLev][iSt] = new TH2F(Form("hTrkPosCor%d_St%d", iLev, iSt),
1147 Form("Lev%d-St%d position correlation ; "
1148 "#DeltaX (cm); #DeltaY (cm)",
1149 iLev, iSt),
1150 100, -10, 10, 100, -10, 10);
1151 }
1152 }
1153 if (fiStationUT > -1) {
1154 // histos for Station Under Test
1155 itMapStationZ = fMapStationZ.begin();
1156 for (Int_t iSt = 0; iSt < fiStationUT; iSt++) {
1157 itMapStationZ++;
1158 dSUT_z = itMapStationZ->first;
1159 LOG(info) << "Create SUT histos for station " << fiStationUT << " at distance " << dSUT_z;
1160 const Double_t dSUT_RefDx = 0.15;
1161 const Double_t dSUT_RefDy = 0.3;
1162 const Double_t dNbinX = 100;
1163 const Double_t dNbinY = 100;
1164 const Double_t dNbinZ = 50;
1165
1166 Double_t dSUT_dx = dSUT_RefDx * dSUT_z;
1167 Double_t dSUT_dy = dSUT_RefDy * dSUT_z;
1169 new TH2F("hExtSutXY_Found", Form("StationUnderTest %d found hits ; X (cm); Y (cm)", fiStationUT), dNbinX,
1170 -dSUT_dx, dSUT_dx, dNbinY, -dSUT_dy, dSUT_dy);
1172 new TH2F("hExtSutXY_Missed", Form("StationUnderTest %d missed hits ; X (cm); Y (cm)", fiStationUT), dNbinX,
1173 -dSUT_dx, dSUT_dx, dNbinY, -dSUT_dy, dSUT_dy);
1175 new TH3F("hExtSutXY_DX", Form("StationUnderTest %d #DeltaX ; X (cm); Y (cm); #DeltaX (cm)", fiStationUT),
1176 dNbinX, -dSUT_dx, dSUT_dx, dNbinY, -dSUT_dy, dSUT_dy, dNbinZ, -10., 10.);
1178 new TH3F("hExtSutXY_DY", Form("StationUnderTest %d #DeltaY ; X (cm); Y (cm); #DeltaY (cm)", fiStationUT),
1179 dNbinX, -dSUT_dx, dSUT_dx, dNbinY, -dSUT_dy, dSUT_dy, dNbinZ, -10., 10.);
1181 new TH3F("hExtSutXY_DT", Form("StationUnderTest %d #DeltaT ; X (cm); Y (cm); #DeltaT (ns)", fiStationUT),
1182 dNbinX, -dSUT_dx, dSUT_dx, dNbinY, -dSUT_dy, dSUT_dy, dNbinZ, -50., 50.);
1183 }
1184 } // StationUT end
1185
1186 gDirectory->cd(oldir->GetPath()); // <= To prevent histos from being sucked in by the param file of the TRootManager!
1187
1189}
1190
1192{
1193 fVTX_T = 0.; //reset
1194 fVTX_X = 0.;
1195 fVTX_Y = 0.;
1196 fVTX_Z = 0.;
1197 fVTXNorm = 0.;
1198 Int_t fMinNofHits = 3;
1199
1200 for (Int_t iTrk = 0; iTrk < fTofTrackArrayIn->GetEntriesFast(); iTrk++) {
1201 CbmTofTracklet* pTrk = (CbmTofTracklet*) fTofTrackArrayIn->At(iTrk);
1202 if (NULL == pTrk) continue;
1203 Double_t w = pTrk->GetNofHits();
1204
1205 if (w > (Double_t) fMinNofHits) { // for further analysis request minimum number of hits
1206 fVTXNorm += w;
1207 fVTX_T += w * pTrk->GetFitT(0.);
1208 fVTX_X += w * pTrk->GetFitX(0.);
1209 fVTX_Y += w * pTrk->GetFitY(0.);
1210 }
1211 }
1212 if (fVTXNorm > 0.) {
1213 fVTX_T /= fVTXNorm;
1214 fVTX_X /= fVTXNorm;
1215 fVTX_Y /= fVTXNorm;
1216 fVTX_Z /= fVTXNorm;
1217 }
1218 LOG(debug) << Form("CbmTofExtendTracks::FindVertex: N %3.0f, T %6.2f, "
1219 "X=%6.2f, Y=%6.2f Z=%6.2f ",
1221}
1222
1224{
1225 // event selection: limit number of hits per added station
1226 Int_t iAddStations = fiAddStations;
1227 while (iAddStations > 0) {
1228 Int_t iAddNextStation = iAddStations % 100;
1229 iAddStations /= 100;
1230 if ((Int_t) fvAllHitPointer[iAddNextStation].size() > fiCutStationMaxHitMul) return;
1231 }
1232
1233 // limit track multiplicity
1234 Int_t NTrkTof = static_cast<Int_t>(tEvent->GetNofData(ECbmDataType::kTofTracklet));
1235 if (NTrkTof > fiNTrkTofMax) return;
1236
1237 Int_t iLev = 0;
1238 //
1240
1242
1245
1247
1248
1249 // correlation with TOF Tracklets
1250 for (Int_t iTr = 0; iTr < NTrkTof; iTr++) {
1251 Int_t iTrkIndex = static_cast<Int_t>(tEvent->GetIndex(ECbmDataType::kTofTracklet, iTr));
1252 CbmTofTracklet* tTrk = dynamic_cast<CbmTofTracklet*>(fTofTrackArrayIn->At(iTrkIndex));
1253
1254 // tTrk->PrintInfo();
1255
1256 if (tTrk->GetNofHits() < fiTrkHitsMin) continue;
1257
1258 // correlation with all TofHits
1259 for (UInt_t iSt = 0; iSt < fvTofHitIndex.size(); iSt++) {
1260 for (UInt_t iHit = 0; iHit < fvTofHitIndex[iSt].size(); iHit++) {
1261 CbmTofHit* tHit = dynamic_cast<CbmTofHit*>(fTofHitArrayIn->At(fvTofHitIndex[iSt][iHit]));
1262 if (NULL == tHit) {
1263 LOG(warn) << Form("Invalid TofHit %d(%lu) ind %d in station %d(%lu)", iHit, fvTofHitIndex[iSt].size(),
1264 fvTofHitIndex[iSt][iHit], iSt, fvTofHitIndex.size());
1265 assert(tHit);
1266 }
1267 Int_t iTrSt = fMapStationZ[(Int_t)(tHit->GetZ() / dStDZ)] % 100;
1268 if (iTrSt < 0 || iTrSt >= (Int_t) fMapStationZ.size()) {
1269 LOG(error) << "Invalid station index " << iTrSt;
1270 continue;
1271 }
1272
1273 Double_t dDX = tHit->GetX() + fvXoff[iTrSt] - tTrk->GetFitX(tHit->GetZ());
1274 if (TMath::Abs(dDX) > fdTrkCutDX) continue;
1275 Double_t dDY = tHit->GetY() + fvYoff[iTrSt] - tTrk->GetFitY(tHit->GetZ());
1276 if (TMath::Abs(dDY) > fdTrkCutDY) continue;
1277 Double_t dDT = tHit->GetTime() + fvToff[iTrSt] - tTrk->GetFitT(tHit->GetZ());
1278 if (TMath::Abs(dDT) > fdTrkCutDT) continue;
1279
1280 fhPosCorTrkTof->Fill(dDX, dDY);
1281 // fill tracking station histos
1282 fhTrkStationDX[iLev]->Fill(iTrSt, dDX);
1283 fhTrkStationDY[iLev]->Fill(iTrSt, dDY);
1284 fhTrkStationDT[iLev]->Fill(iTrSt, dDT);
1285 fhTrkStationNHits[iLev]->Fill(iTrSt, (Double_t) fvTofHitIndex[iSt].size());
1286 }
1287 } // Tof Station loop end
1288
1289 // correlation with all StsHits
1290 for (UInt_t iSt = 0; iSt < fvStsHitIndex.size(); iSt++) {
1291 for (UInt_t iHit = 0; iHit < fvStsHitIndex[iSt].size(); iHit++) {
1292 CbmPixelHit* tHit = dynamic_cast<CbmPixelHit*>(fStsHitArrayIn->At(fvStsHitIndex[iSt][iHit]));
1293 if (NULL == tHit) {
1294 LOG(warn) << Form("Invalid StsHit %d(%lu) ind %d in station %d(%lu)", iHit, fvStsHitIndex[iSt].size(),
1295 fvStsHitIndex[iSt][iHit], iSt, fvStsHitIndex.size());
1296 assert(tHit);
1297 }
1298 Int_t iTrSt = fMapStationZ[(Int_t)(tHit->GetZ() / dStDZ)] % 100;
1299
1300 Double_t dDX = tHit->GetX() + fvXoff[iTrSt] - tTrk->GetFitX(tHit->GetZ());
1301 if (TMath::Abs(dDX) > fdTrkCutDX) continue;
1302 Double_t dDY = tHit->GetY() + fvYoff[iTrSt] - tTrk->GetFitY(tHit->GetZ());
1303 if (TMath::Abs(dDY) > fdTrkCutDY) continue;
1304 Double_t dDT = tHit->GetTime() + fvToff[iTrSt] - tTrk->GetFitT(tHit->GetZ());
1305 if (TMath::Abs(dDT) > fdTrkCutDT) continue;
1306
1307 fhPosCorTrkSts->Fill(dDX, dDY);
1308 // fill tracking station histos
1309 fhTrkStationDX[iLev]->Fill(iTrSt, dDX);
1310 fhTrkStationDY[iLev]->Fill(iTrSt, dDY);
1311 fhTrkStationDT[iLev]->Fill(iTrSt, dDT);
1312 fhTrkStationNHits[iLev]->Fill(iTrSt, (Double_t) fvStsHitIndex[iSt].size());
1313 }
1314 } // Sts Station loop end
1315
1316 // correlation with all MuchHits
1317 for (UInt_t iSt = 0; iSt < fvMuchHitIndex.size(); iSt++) {
1318 for (UInt_t iHit = 0; iHit < fvMuchHitIndex[iSt].size(); iHit++) {
1319 CbmPixelHit* tHit = dynamic_cast<CbmPixelHit*>(fMuchHitArrayIn->At(fvMuchHitIndex[iSt][iHit]));
1320 if (NULL == tHit) {
1321 LOG(warn) << Form("Invalid MuchHit %d(%lu) ind %d in station %d(%lu)", iHit, fvMuchHitIndex[iSt].size(),
1322 fvMuchHitIndex[iSt][iHit], iSt, fvMuchHitIndex.size());
1323 assert(tHit);
1324 }
1325 Int_t iTrSt = fMapStationZ[(Int_t)(tHit->GetZ() / dStDZ)] % 100;
1326
1327 LOG(info) << Form("Proc ev %d, Trk %d, St %d, iLe%d MtHit: ", fiEvent, iTr, iTrSt, iLev) << tHit->ToString();
1328
1329 Double_t dDX = tHit->GetX() + fvXoff[iTrSt] - tTrk->GetFitX(tHit->GetZ());
1330 if (TMath::Abs(dDX) > fdTrkCutDX) continue;
1331 Double_t dDY = tHit->GetY() + fvYoff[iTrSt] - tTrk->GetFitY(tHit->GetZ());
1332 if (TMath::Abs(dDY) > fdTrkCutDY) continue;
1333 Double_t dDT = tHit->GetTime() + fvToff[iTrSt] - tTrk->GetFitT(tHit->GetZ());
1334 if (TMath::Abs(dDT) > fdTrkCutDT) continue;
1335 /*
1336 LOG(info)<<"Got MuchCorI Tr "<<iTr
1337 <<", St "<<iTrSt
1338 <<", Hit "<< iHit <<", ind "<<fvMuchHitIndex[iSt][iHit]
1339 <<", poi "<< tHit
1340 <<", dx " << dDX
1341 <<", dy " << dDY;
1342 */
1343 LOG(info) << Form("Proc ev %d, Trk %d, St %d, iLe%d, DX %6.3f, dY %6.3f, MtHit: ", fiEvent, iTr, iTrSt, iLev,
1344 dDX, dDY)
1345 << tHit->ToString();
1346
1347 fhPosCorTrkMuch->Fill(dDX, dDY);
1348 // fill tracking station histos
1349 fhTrkStationDX[iLev]->Fill(iTrSt, dDX);
1350 fhTrkStationDY[iLev]->Fill(iTrSt, dDY);
1351 fhTrkStationDT[iLev]->Fill(iTrSt, dDT);
1352 fhTrkStationNHits[iLev]->Fill(iTrSt, (Double_t) fvMuchHitIndex[iSt].size());
1353 }
1354 } // Much Station loop end
1355
1356 // correlation with all RichHits
1357 for (UInt_t iSt = 0; iSt < fvRichHitIndex.size(); iSt++) {
1358 for (UInt_t iHit = 0; iHit < fvRichHitIndex[iSt].size(); iHit++) {
1359 CbmPixelHit* tHit = dynamic_cast<CbmPixelHit*>(fRichHitArrayIn->At(fvRichHitIndex[iSt][iHit]));
1360 if (NULL == tHit) {
1361 LOG(warn) << Form("Invalid RichHit %d(%lu) ind %d in station %d(%lu)", iHit, fvRichHitIndex[iSt].size(),
1362 fvRichHitIndex[iSt][iHit], iSt, fvRichHitIndex.size());
1363 assert(tHit);
1364 }
1365 Int_t iTrSt = fMapStationZ[(Int_t)(tHit->GetZ() / dStDZ)] % 100;
1366
1367 Double_t dDX = tHit->GetX() + fvXoff[iTrSt] - tTrk->GetFitX(tHit->GetZ());
1368 if (TMath::Abs(dDX) > fdTrkCutDX) continue;
1369 Double_t dDY = tHit->GetY() + fvYoff[iTrSt] - tTrk->GetFitY(tHit->GetZ());
1370 if (TMath::Abs(dDY) > fdTrkCutDY) continue;
1371 Double_t dDT = tHit->GetTime() + fvToff[iTrSt] - tTrk->GetFitT(tHit->GetZ());
1372 if (TMath::Abs(dDT) > fdTrkCutDT) continue;
1373
1374 fhPosCorTrkRich->Fill(dDX, dDY);
1375 // fill tracking station histos
1376 fhTrkStationDX[iLev]->Fill(iTrSt, dDX);
1377 fhTrkStationDY[iLev]->Fill(iTrSt, dDY);
1378 fhTrkStationDT[iLev]->Fill(iTrSt, dDT);
1379 fhTrkStationNHits[iLev]->Fill(iTrSt, (Double_t) fvRichHitIndex[iSt].size());
1380 }
1381 } // Rich Station loop end
1382
1383 } // TOF Tracklet loop end
1384
1385 iLev = 1;
1386 // Fully assembled tracks
1387 UInt_t NSt = fMapStationZ.size();
1388 UInt_t NTrk = fvTrkPar.size();
1389 if (NTrkTof < 0) NTrkTof = 0;
1390 assert(NTrk == (UInt_t) NTrkTof);
1391
1392 for (UInt_t iTr = 0; iTr < NTrk; iTr++) {
1393 //select tracks
1394 // minimal number of TOF hits
1395 Int_t NTofHits = 0;
1396 for (UInt_t iH = 0; iH < fvTrkCalHits[iTr].size(); iH++) {
1397 if (fMapStationZ[(Int_t)(fvTrkCalHits[iTr][iH]->GetZ())] / 100 == Int_t(ECbmModuleId::kTof)) NTofHits++;
1398 }
1399 if (NTofHits < fiTrkHitsMin) continue;
1400 // contain requested stations
1401 if (fiReqStations > -1) {
1402 Int_t iReqStations = fiReqStations;
1403 Int_t iHit = (Int_t) fvTrkCalHits[iTr].size() - 1;
1404 Int_t iReqStation = 0;
1405 while (iReqStations > 0) {
1406 iReqStation = iReqStations % 100;
1407 iHit = (Int_t) fvTrkCalHits[iTr].size() - 1;
1408 // LOG(info)<<"Check tr "<<iTr<<" for station "<<iReqStation
1409 // <<" starting from hit "<<iHit<<" in station "
1410 // <<fMapStationZ[(Int_t)(fvTrkCalHits[iTr][iHit]->GetZ())]%100;
1411 iReqStations /= 100;
1412 for (; iHit > 0; iHit--)
1413 if (iReqStation == fMapStationZ[(Int_t)(fvTrkCalHits[iTr][iHit]->GetZ())] % 100) {
1414 // if (iReqStations ==0 )
1415 // LOG(info)<<"Check tr "<<iTr<<" accepted for analysis";
1416 break;
1417 }
1418 // LOG(info)<<"Check last Hit "<<iHit<<", ReqSt "<<iReqStations
1419 // <<"("<<fiReqStations<<") ";
1420 if (iReqStation != 0 && iHit == 0) break; // station not found
1421 // LOG(info)<<"CheckI with "<<iReqStations<<","<<iHit;
1422 }
1423 // LOG(info)<<"CheckII with "<<iReqStations<<","<<iHit;
1424 if (iReqStation != 0 && iHit == 0) continue; //skip this track
1425 }
1426 /*
1427 LOG(info) << "CheckIII accepted tr " << iTr << " with " << fiReqStations;
1428 for (UInt_t iHit = 0; iHit < fvTrkCalHits[iTr].size(); iHit++) {
1429 LOG(info) << " Hit " << iHit << ", station " << fMapStationZ[(Int_t)(fvTrkCalHits[iTr][iHit]->GetZ())] % 100;
1430 }
1431 */
1432 fhExt_TrkSizChiSq[iLev]->Fill((Double_t) fvTrkCalHits[iTr].size(), fvTrkPar[iTr]->GetChiSq());
1433
1434 fhExt_TrkSizVel[iLev]->Fill((Double_t) fvTrkCalHits[iTr].size(), 1. / (fvTrkPar[iTr]->GetTt()));
1435
1436 // Pulls
1437 CbmTofTrackletParam* tTrkPar = fvTrkPar[iTr];
1438 Bool_t bSUT_Found = kFALSE;
1439 for (UInt_t iHit = 0; iHit < fvTrkCalHits[iTr].size(); iHit++) {
1440 CbmPixelHit* tHit = fvTrkCalHits[iTr][iHit];
1441 Int_t iTrSt = fMapStationZ[(Int_t)(tHit->GetZ() / dStDZ)] % 100;
1442
1443 Double_t dDX = tHit->GetX() - GetFitX(tHit->GetZ(), tTrkPar);
1444 Double_t dDY = tHit->GetY() - GetFitY(tHit->GetZ(), tTrkPar);
1445 Double_t dDT = tHit->GetTime() - GetFitT(tHit->GetZ(), tTrkPar);
1446
1447 if (TMath::Abs(dDX) > fdTrkCutDX) continue;
1448 if (TMath::Abs(dDY) > fdTrkCutDY) continue;
1449 if (TMath::Abs(dDT) > fdTrkCutDT) continue;
1450
1451 fhTrkPullDX[iLev]->Fill(iTrSt, dDX);
1452 fhTrkPullDY[iLev]->Fill(iTrSt, dDY);
1453 fhTrkPullDT[iLev]->Fill(iTrSt, dDT);
1454
1455 if (iTrSt == fiStationUT) { // hit in SUT
1456 bSUT_Found = kTRUE;
1457 fhExtSutXY_Found->Fill(tHit->GetX(), tHit->GetY());
1458 fhExtSutXY_DX->Fill(tHit->GetX(), tHit->GetY(), dDX);
1459 fhExtSutXY_DY->Fill(tHit->GetX(), tHit->GetY(), dDY);
1460 fhExtSutXY_DT->Fill(tHit->GetX(), tHit->GetY(), dDT);
1461 }
1462
1463 // pulls against original tracklets
1464 Int_t iTrkInd = static_cast<Int_t>(tEvent->GetIndex(ECbmDataType::kTofTracklet, iTr));
1465 CbmTofTracklet* pTrk = dynamic_cast<CbmTofTracklet*>(fTofTrackArrayIn->At(iTrkInd));
1466 dDX = tHit->GetX() - pTrk->GetFitX(tHit->GetZ());
1467 dDY = tHit->GetY() - pTrk->GetFitY(tHit->GetZ());
1468 dDT = tHit->GetTime() - pTrk->GetFitT(tHit->GetZ());
1469 fhTrkPullDX[0]->Fill(iTrSt, dDX);
1470 fhTrkPullDY[0]->Fill(iTrSt, dDY);
1471 fhTrkPullDT[0]->Fill(iTrSt, dDT);
1472 fhTrkStationDXDY[0][iTrSt]->Fill(dDX, dDY);
1473 LOG(debug) << Form("Proc ev %d, Trk %d, St %d, Lev%d, DX %6.3f, dY %6.3f, MtHit: ", fiEvent, iTr, iTrSt, 0, dDX,
1474 dDY)
1475 << tHit->ToString();
1476 } // fvTrkCalHits loop end
1477
1478 if (bSUT_Found == kFALSE) {
1479 fhExtSutXY_Missed->Fill(GetFitX(dSUT_z, tTrkPar), GetFitY(dSUT_z, tTrkPar));
1480 }
1481
1482 //Deviations to all hits in all stations
1483 for (UInt_t iSt = 0; iSt < NSt; iSt++) {
1484 for (UInt_t iHit = 0; iHit < fvAllHitPointer[iSt].size(); iHit++) {
1485 CbmPixelHit* tHit = fvAllHitPointer[iSt][iHit];
1486 /*
1487 LOG(info)<<"Got AllHit for Tr "<<iTr
1488 <<", St "<<iSt
1489 <<", StZ "<<fMapStationZ[(Int_t)(tHit->GetZ()/dStDZ)]%100
1490 <<", Hit "<< iHit <<", poi "<< tHit;
1491 LOG(info)<<"tHit: "<<tHit->ToString();
1492 */
1493 assert((Int_t) iSt == fMapStationZ[(Int_t)(tHit->GetZ() / dStDZ)] % 100);
1494
1495 Double_t dDX = tHit->GetX() - GetFitX(tHit->GetZ(), tTrkPar);
1496 if (TMath::Abs(dDX) > fdTrkCutDX) continue;
1497 Double_t dDY = tHit->GetY() - GetFitY(tHit->GetZ(), tTrkPar);
1498 if (TMath::Abs(dDY) > fdTrkCutDY) continue;
1499 Double_t dDT = tHit->GetTime() - GetFitT(tHit->GetZ(), tTrkPar);
1500 if (TMath::Abs(dDT) > fdTrkCutDT) continue;
1501
1502 LOG(debug) << Form("Proc ev %d, Trk %d, St %d, Lev%d, DX %6.3f, dY %6.3f, MtHit: ", fiEvent, iTr, iSt, iLev,
1503 dDX, dDY)
1504 << tHit->ToString();
1505
1506 // fill tracking station histos
1507 fhTrkStationDX[iLev]->Fill(iSt, dDX);
1508 fhTrkStationDY[iLev]->Fill(iSt, dDY);
1509 fhTrkStationDT[iLev]->Fill(iSt, dDT);
1510 fhTrkStationNHits[iLev]->Fill(iSt, (Double_t) fvAllHitPointer[iSt].size());
1511 fhTrkStationDXDY[iLev][iSt]->Fill(dDX, dDY);
1512 }
1513 }
1514 } // loop on extended tracks end
1515
1516} // FillHist end
1517
1518void CbmTofExtendTracks::Line3Dfit(std::vector<CbmPixelHit*> pHits, CbmTofTrackletParam* pTrkPar)
1519{
1520 TGraph2DErrors* gr = new TGraph2DErrors();
1521 // Fill the 2D graph
1522 // generate graph with the 3d points
1523 for (UInt_t N = 0; N < pHits.size(); N++) {
1524 double x, y, z = 0;
1525 x = pHits[N]->GetX();
1526 y = pHits[N]->GetY();
1527 z = pHits[N]->GetZ();
1528 gr->SetPoint(N, x, y, z);
1529 double dx, dy, dz = 0.;
1530 dx = pHits[N]->GetDx();
1531 dy = pHits[N]->GetDy();
1532 dz = pHits[N]->GetDz(); //FIXME
1533 gr->SetPointError(N, dx, dy, dz);
1534 LOG(debug) << "Line3Dfit add N = " << N << ",\t" << x << ",\t" << y << ",\t" << z << ",\t" << dx << ",\t" << dy
1535 << ",\t" << dz;
1536 }
1537 // fit the graph now
1538
1539 Double_t pStart[4] = {0., 0., 0., 0.};
1540 pStart[0] = pTrkPar->GetX();
1541 pStart[1] = pTrkPar->GetTx();
1542 pStart[2] = pTrkPar->GetY();
1543 pStart[3] = pTrkPar->GetTy();
1544 LOG(debug) << "Line3Dfit init: X0 " << pStart[0] << ", TX " << pStart[1] << ", Y0 " << pStart[2] << ", TY "
1545 << pStart[3];
1546
1547 fMinuit.DoFit(gr, pStart);
1548 //gr->Draw("err p0");
1549 gr->Delete();
1550 Double_t* dRes;
1551 dRes = fMinuit.GetParFit();
1552 LOG(debug) << "Line3Dfit result(" << pHits.size() << ") " << gMinuit->fCstatu << " : X0 "
1553 << Form(" %7.4f, TX %7.4f, Y0 %7.4f, TY %7.4f, Chi2DoF %7.4f ", dRes[0], dRes[1], dRes[2], dRes[3],
1555
1556 pTrkPar->SetX(dRes[0]);
1557 pTrkPar->SetY(dRes[2]);
1558 pTrkPar->SetZ(0.);
1559 pTrkPar->SetTx(dRes[1]);
1560 pTrkPar->SetTy(dRes[3]);
1561 pTrkPar->SetQp(1.);
1562 pTrkPar->SetChiSq(fMinuit.GetChi2DoF()); // / (Double_t)pHits.size());
1563}
1564
1566{
1567 return fTrkPar->GetX() + fTrkPar->GetTx() * (dZ - fTrkPar->GetZ());
1568}
1569
1571{
1572 return fTrkPar->GetY() + fTrkPar->GetTy() * (dZ - fTrkPar->GetZ());
1573}
1574
1576{
1577 return fTrkPar->GetT()
1578 + fTrkPar->GetTt() * (dZ - fTrkPar->GetZ())
1579 * TMath::Sqrt(1. + fTrkPar->GetTx() * fTrkPar->GetTx() + fTrkPar->GetTy() * fTrkPar->GetTy());
1580}
1581
1583{
1584
1585 LOG(debug) << "Add " << fvAllHitPointer[iStation].size() << " hits from station " << iStation << " to "
1586 << fvTrkPar.size() << " tracks ";
1587
1588 const Int_t NCand = 100;
1589 // tbd: get matching width from histos !!
1590 Double_t tSIGX = 1.;
1591 Double_t tSIGY = 1.;
1592 Double_t tSIGT = 1000.;
1593 UInt_t NTr = fvTrkPar.size();
1594 Double_t MatchChi2Min[NCand] = {NCand * fdChi2Max};
1595 Int_t MatchHit[NCand] = {NCand * -1};
1596 Int_t MatchTrk[NCand] = {NCand * -1};
1597 // find list of best matches for each track, max NCand
1598 Int_t iMatch = 0;
1599 for (UInt_t iTr = 0; iTr < NTr; iTr++) {
1600 for (UInt_t iHit = 0; iHit < fvAllHitPointer[iStation].size(); iHit++) {
1601 CbmPixelHit* tHit = fvAllHitPointer[iStation][iHit];
1602 LOG(debug) << " Hit " << iHit << ": X " << tHit->GetX() << ", Y " << tHit->GetY() << ", Z " << tHit->GetZ()
1603 << ", T " << tHit->GetTime();
1604 Double_t MatchChi2 =
1605 TMath::Power((tHit->GetX() + fvXoff[iStation] - GetFitX(tHit->GetZ(), fvTrkPar[iTr])) / tSIGX, 2)
1606 + TMath::Power((tHit->GetY() + fvYoff[iStation] - GetFitY(tHit->GetZ(), fvTrkPar[iTr])) / tSIGY, 2);
1607 TMath::Power((tHit->GetTime() + fvToff[iStation] - GetFitT(tHit->GetZ(), fvTrkPar[iTr])) / tSIGT, 2);
1608 MatchChi2 /= 3.;
1609 LOG(debug) << "Match Tr " << iTr << ", Hit " << iHit << ": " << MatchChi2;
1610 if (MatchChi2 < fdChi2Max) {
1611 Int_t iCand = 0;
1612 for (; iCand < iMatch + 1; iCand++) {
1613 if (MatchChi2 < MatchChi2Min[iCand]) {
1614 LOG(debug) << "Better Match found for iTr " << iTr << " Chi2 " << MatchChi2;
1615 for (Int_t i = iMatch - 1; i >= iCand; i--) {
1616 MatchChi2Min[i + 1] = MatchChi2Min[i];
1617 MatchHit[i + 1] = MatchHit[i];
1618 MatchTrk[i + 1] = MatchTrk[i];
1619 }
1620 MatchChi2Min[iCand] = MatchChi2;
1621 MatchHit[iCand] = (Int_t) iHit;
1622 MatchTrk[iCand] = (Int_t) iTr;
1623 break;
1624 }
1625 } // end candidate loop
1626 iMatch++;
1627 if (iMatch == NCand) iMatch--; // cutoff
1628 LOG(debug) << "New Match " << iMatch << " stored as candidate " << iCand << " for Tr " << iTr << " with hit "
1629 << iHit;
1630 }
1631 } // station hit loop end
1632 } // track loop end
1633
1634 // append best matches to track
1635 for (Int_t iCand = 0; iCand < iMatch; iCand++) {
1636 if (MatchHit[iCand] > -1) {
1637 CbmPixelHit* tHit = fvAllHitPointer[iStation][MatchHit[iCand]];
1638 /*
1639 LOG(info)<<"Add hit "<<MatchHit[iCand]
1640 <<", station "<<fMapStationZ[(Int_t)(tHit->GetZ()/dStDZ)]%100
1641 <<" to track "<<MatchTrk[iCand]
1642 <<" with chi2 "<<MatchChi2Min[iCand];
1643 */
1644 //LOG(info)<<"MatHit: "<<tHit->ToString();
1645
1646 fvTrkCalHits[MatchTrk[iCand]].push_back(dynamic_cast<CbmPixelHit*>(tHit));
1647 // prevent match with other tracks
1648 for (Int_t i = iCand + 1; i < iMatch; i++) {
1649 if (MatchHit[i] == MatchHit[iCand] || MatchTrk[i] == MatchTrk[iCand]) MatchHit[i] = -1;
1650 }
1651 // refit track, update TrkPar
1652 fvTrkPar[MatchTrk[iCand]]->SetX(0.);
1653 fvTrkPar[MatchTrk[iCand]]->SetY(0.);
1654
1655 Line3Dfit(fvTrkCalHits[MatchTrk[iCand]], fvTrkPar[MatchTrk[iCand]]);
1656 // Inspect pulls (debugging)
1657 /*
1658 Double_t dDX=tHit->GetX() - GetFitX(tHit->GetZ(),fvTrkPar[MatchTrk[iCand]]);
1659 Double_t dDY=tHit->GetY() - GetFitY(tHit->GetZ(),fvTrkPar[MatchTrk[iCand]]);
1660 LOG(info)<<"MatchPulls X "<<dDX<<", Y "<<dDY;
1661 */
1662 }
1663 } // end of loop on match candidates
1664}
@ kTof
Time-of-flight Detector.
@ kSts
Silicon Tracking System.
@ kMuch
Muon detection system.
@ kRich
Ring-Imaging Cherenkov Detector.
Char_t * gr
static TFile * fHist
TClonesArray * fEventsColl
static Int_t fiTS
const Int_t NDefSetup
ClassImp(CbmTofExtendTracks)
static LKFMinuit fMinuit
const double dStDZ
static Double_t dSUT_z
Data class for track parameters.
static constexpr size_t size()
Definition KfSimdPseudo.h:2
Class characterising one event by a collection of links (indices) to data objects,...
Definition CbmEvent.h:34
size_t GetNofData() const
Definition CbmEvent.cxx:53
uint32_t GetIndex(ECbmDataType type, uint32_t iData)
Definition CbmEvent.cxx:42
double GetTime() const
Definition CbmHit.h:76
void SetZ(double z)
Definition CbmHit.h:80
double GetZ() const
Definition CbmHit.h:71
void SetTime(double time)
Definition CbmHit.h:85
void SetX(double x)
Definition CbmPixelHit.h:92
double GetY() const
Definition CbmPixelHit.h:74
virtual std::string ToString() const
Inherited from CbmBaseHit.
void SetY(double y)
Definition CbmPixelHit.h:93
void SetPositionError(const TVector3 &dpos)
double GetX() const
Definition CbmPixelHit.h:73
std::vector< Double_t > fvZsig
std::vector< Double_t > fvStsStationZ
std::vector< Double_t > fvXsig
std::vector< Double_t > fvToff
virtual void Line3Dfit(std::vector< CbmPixelHit * >, CbmTofTrackletParam *)
virtual InitStatus Init()
std::vector< Double_t > fvMuchStationZ
std::vector< CbmTofTrackletParam * > fvTrkPar
std::vector< std::vector< CbmPixelHit * > > fvTrkCalHits
virtual void CreateHistograms()
virtual void TrkAddStation(Int_t iStation)
std::vector< TH2 * > fhTrkStationDY
virtual void ExecExtend(Option_t *opt, CbmEvent *tEvent=NULL)
std::map< Int_t, Int_t >::iterator itMapStationZ
std::vector< TH2 * > fhTrkPullDY
std::vector< TH2 * > fhTrkStationDT
std::vector< TH2 * > fhTrkPullDX
std::vector< std::vector< TH2 * > > fhTrkStationDXDY
virtual void FillHistograms(CbmEvent *tEvent=NULL)
std::vector< std::vector< CbmPixelHit * > > fvAllHitPointer
std::vector< std::vector< Int_t > > fvTofHitIndex
std::vector< Double_t > fvXoff
TClonesArray * fMuchHitArrayIn
std::vector< Double_t > fvTofStationZ
TClonesArray * fTofHitArrayIn
std::vector< TH2 * > fhTrkStationNHits
Double_t GetFitY(Double_t, CbmTofTrackletParam *)
std::vector< Double_t > fvTsig
virtual void SetParContainers()
std::vector< Double_t > fvYsig
std::vector< TH2 * > fhTrkStationDX
std::map< Int_t, Int_t > fMapStationZ
TClonesArray * fTofTrackArrayIn
TClonesArray * fRichHitArrayIn
std::vector< TH2 * > fhTrkPullDT
Double_t GetFitT(Double_t, CbmTofTrackletParam *)
std::vector< TH2 * > fhExt_TrkSizVel
Double_t GetFitX(Double_t, CbmTofTrackletParam *)
std::vector< std::vector< Int_t > > fvMuchHitIndex
std::vector< Double_t > fvYoff
static CbmTofExtendTracks * fInstance
TClonesArray * fStsHitArrayIn
TClonesArray * fEventsColl
std::vector< Double_t > fvRichStationZ
std::vector< Double_t > fvZoff
CbmTofTrackletTools * fTrackletTools
std::vector< std::vector< Int_t > > fvRichHitIndex
std::vector< TH2 * > fhExt_TrkSizChiSq
std::vector< std::vector< Int_t > > fvStsHitIndex
std::vector< TH2 * > fhTrkStationDZ
virtual void Exec(Option_t *opt)
contains fits and resolution functions
Provides information on attaching a TofHit to a TofTrack.
double GetFitY(double Z)
int32_t GetTofHitIndex(int32_t ind) const
double GetTt() const
double GetTrackTx() const
int32_t GetNofHits() const
double GetFitT(double R)
double GetTrackTy() const
double GetT0() const
CbmTofTrackletParam * GetTrackParameter()
double GetFitX(double Z)
int DoFit(TGraph2DErrors *gr, double pStart[])
Definition LKFMinuit.cxx:36
int Initialize()
Definition LKFMinuit.cxx:23
double * GetParFit()
Definition LKFMinuit.h:37
double GetChi2DoF()
Definition LKFMinuit.h:39