CbmRoot
Loading...
Searching...
No Matches
CbmTofSimpClusterizer.cxx
Go to the documentation of this file.
1/* Copyright (C) 2013-2021 Facility for Antiproton and Ion Research in Europe, Darmstadt
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Pierre-Alain Loizeau [committer], Norbert Herrmann */
4
11
12// TOF Classes and includes
13#include "CbmTofAddress.h" // in cbmdata/tof
14#include "CbmTofCell.h" // in tof/TofData
15#include "CbmTofCreateDigiPar.h"
16#include "CbmTofDetectorId_v12b.h" // in cbmdata/tof
17#include "CbmTofDetectorId_v14a.h" // in cbmdata/tof
18#include "CbmTofDetectorId_v21a.h" // in cbmdata/tof
19#include "CbmTofDigi.h" // in cbmdata/tof
20#include "CbmTofDigiBdfPar.h" // in tof/TofParam
21#include "CbmTofDigiPar.h" // in tof/TofParam
22#include "CbmTofGeoHandler.h" // in tof/TofTools
23#include "CbmTofHit.h" // in cbmdata/tof
24#include "CbmTofPoint.h" // in cbmdata/tof
25
26// CBMroot classes and includes
27#include "CbmDigiManager.h"
28#include "CbmEvent.h"
29#include "CbmMCTrack.h"
30#include "CbmMatch.h"
31
32// FAIR classes and includes
33#include "FairEventHeader.h" // from CbmStsDigitize, for GetEventInfo
34#include "FairMCEventHeader.h" // from CbmStsDigitize, for GetEventInfo
35#include "FairRootManager.h"
36#include "FairRunAna.h"
37#include "FairRunSim.h" // from CbmStsDigitize, for GetEventInfo
38#include "FairRuntimeDb.h"
39
40#include <Logger.h>
41
42// ROOT Classes and includes
43#include "TClonesArray.h"
44#include "TDirectory.h"
45#include "TF2.h"
46#include "TGeoManager.h"
47#include "TH1.h"
48#include "TH2.h"
49#include "TLine.h"
50#include "TMath.h"
51#include "TProfile.h"
52#include "TROOT.h"
53#include "TRandom3.h"
54#include "TVector3.h"
55
56#include <TFile.h>
57
58// C++ Classes and includes
59#include <iomanip>
60#include <iostream>
61
62using std::fixed;
63using std::left;
64using std::pair;
65using std::right;
66using std::setprecision;
67using std::setw;
68using std::stringstream;
69
70// const Int_t DetMask = 4194303; (VF) not used
71const Int_t nbClWalkBinX = 20;
72//const Int_t nbClWalkBinY=41; (VF) not used // choose odd number to have central bin symmetric around 0
73//const Double_t WalkNHmin=100; (VF) not used // minimal number of hits in bin for walk correction
74Double_t TOTMax = 5.E4;
75Double_t TOTMin = 2.E4;
76const Double_t TTotMean = 2.E4;
77
78const Int_t nbClDelTofBinX = 50;
79//const Int_t nbClDelTofBinY=49; (VF) not used
80//const Double_t DelTofMax=5000.; (VF) not used
81
82//const Int_t nbCldXdYBinX=49; (VF) not used
83//const Int_t nbCldXdYBinY=49; (VF) not used
84//const Double_t dXdYMax=10.; (VF) not used
85
86const Int_t iNTrg = 1;
87
88//const Double_t Zref = 200.; (VF) not used // distance of projection plane to target
89// C++ Classes and includes
90
91/************************************************************************************/
93 : FairTask("CbmTofSimpClusterizer")
94 , fGeoHandler(new CbmTofGeoHandler())
95 , fTofId(NULL)
96 , fDigiPar(NULL)
97 , fChannelInfo(NULL)
98 , fDigiBdfPar(NULL)
99 , fdParFeeTimeRes(0.0)
100 , fdParSystTimeRes(0.0)
101 , fTofPointsColl(NULL)
102 , fMcTracksColl(NULL)
103 , fDigiMan(nullptr)
104 , fTofHitsColl(NULL)
105 , fTofDigiMatchColl(NULL)
106 , fiNbHits(0)
107 , fVerbose(1)
108 , fStorDigiExp()
109 , fStorDigiInd()
110 , fviClusterMul()
111 , fviClusterSize()
112 , fviTrkMul()
113 , fvdX()
114 , fvdY()
115 , fvdDifX()
116 , fvdDifY()
117 , fvdDifCh()
118 , fsHistoOutFilename("")
119 , fhClustBuildTime(NULL)
120 , fhHitsPerTracks(NULL)
121 , fhPtsPerHit(NULL)
122 , fhTimeResSingHits(NULL)
123 , fhTimeResSingHitsB(NULL)
124 , fhTimePtVsHits(NULL)
125 , fhClusterSize(NULL)
126 , fhClusterSizeType(NULL)
127 , fhTrackMul(NULL)
128 , fhClusterSizeMulti(NULL)
129 , fhTrk1MulPos(NULL)
130 , fhHiTrkMulPos(NULL)
131 , fhAllTrkMulPos(NULL)
132 , fhMultiTrkProbPos(NULL)
133 , fhDigSpacDifClust(NULL)
134 , fhDigTimeDifClust(NULL)
135 , fhDigDistClust(NULL)
136 , fhClustSizeDifX(NULL)
137 , fhClustSizeDifY(NULL)
138 , fhChDifDifX(NULL)
139 , fhChDifDifY(NULL)
140 , fhRpcDigiCor()
141 , fhRpcCluMul()
142 , fhRpcSigPropSpeed()
143 , fhRpcCluPosition()
144 , fhRpcCluTOff()
145 , fhRpcCluTrms()
146 , fhRpcCluTot()
147 , fhRpcCluSize()
148 , fhRpcCluAvWalk()
149 , fhRpcCluWalk()
150 , fhTRpcCluMul()
151 , fhTRpcCluPosition()
152 , fhTRpcCluTOff()
153 , fhTRpcCluTot()
154 , fhTRpcCluSize()
155 , fhTRpcCluAvWalk()
156 , fhTRpcCluDelTof()
157 , fhTRpcCludXdY()
158 , fhTRpcCluWalk()
159 , fhTrgdT()
160 , fvCPSigPropSpeed()
161 , fvCPDelTof()
162 , fvCPTOff()
163 , fvCPTotGain()
164 , fvCPWalk()
165 , fiNbSameSide(0)
166 , fhNbSameSide(NULL)
167 , fhNbDigiPerChan(NULL)
168 , fStart()
169 , fStop()
170 , fTimer()
171 , fiNofEvents(0.)
172 , fdNofHitsTot(0.)
173 , fdTimeTot(0.)
174 , dTRef(0.)
175 , fdTRefMax(0.)
176 , fCalMode(0)
177 , fCalTrg(0)
178 , fCalSmType(0)
179 , fdCaldXdYMax(0.)
180 , fTRefMode(0)
181 , fTRefHits(0)
182 , fPosYMaxScal(0.)
183 , fTRefDifMax(0.)
184 , fTotMax(0.)
185 , fTotMin(0.)
186 , fOutTimeFactor(1.)
187 , fCalParFileName("")
188 , fCalParFile(NULL)
189 , fbMcTrkMonitor(kFALSE)
190{
191}
192
193CbmTofSimpClusterizer::CbmTofSimpClusterizer(const char* name, Int_t verbose)
194 : FairTask(TString(name), verbose)
195 , fGeoHandler(new CbmTofGeoHandler())
196 , fTofId(NULL)
197 , fDigiPar(NULL)
198 , fChannelInfo(NULL)
199 , fDigiBdfPar(NULL)
200 , fdParFeeTimeRes(0.0)
201 , fdParSystTimeRes(0.0)
202 , fTofPointsColl(NULL)
203 , fMcTracksColl(NULL)
204 , fDigiMan(nullptr)
205 , fTofHitsColl(NULL)
206 , fTofDigiMatchColl(NULL)
207 , fiNbHits(0)
208 , fVerbose(verbose)
209 , fStorDigiExp()
210 , fStorDigiInd()
211 , fviClusterMul()
212 , fviClusterSize()
213 , fviTrkMul()
214 , fvdX()
215 , fvdY()
216 , fvdDifX()
217 , fvdDifY()
218 , fvdDifCh()
219 , fsHistoOutFilename("")
220 , fhClustBuildTime(NULL)
221 , fhHitsPerTracks(NULL)
222 , fhPtsPerHit(NULL)
223 , fhTimeResSingHits(NULL)
224 , fhTimeResSingHitsB(NULL)
225 , fhTimePtVsHits(NULL)
226 , fhClusterSize(NULL)
227 , fhClusterSizeType(NULL)
228 , fhTrackMul(NULL)
229 , fhClusterSizeMulti(NULL)
230 , fhTrk1MulPos(NULL)
231 , fhHiTrkMulPos(NULL)
232 , fhAllTrkMulPos(NULL)
233 , fhMultiTrkProbPos(NULL)
234 , fhDigSpacDifClust(NULL)
235 , fhDigTimeDifClust(NULL)
236 , fhDigDistClust(NULL)
237 , fhClustSizeDifX(NULL)
238 , fhClustSizeDifY(NULL)
239 , fhChDifDifX(NULL)
240 , fhChDifDifY(NULL)
241 , fhRpcDigiCor()
242 , fhRpcCluMul()
243 , fhRpcSigPropSpeed()
244 , fhRpcCluPosition()
245 , fhRpcCluTOff()
246 , fhRpcCluTrms()
247 , fhRpcCluTot()
248 , fhRpcCluSize()
249 , fhRpcCluAvWalk()
250 , fhRpcCluWalk()
251 , fhTRpcCluMul()
252 , fhTRpcCluPosition()
253 , fhTRpcCluTOff()
254 , fhTRpcCluTot()
255 , fhTRpcCluSize()
256 , fhTRpcCluAvWalk()
257 , fhTRpcCluDelTof()
258 , fhTRpcCludXdY()
259 , fhTRpcCluWalk()
260 , fhTrgdT()
261 , fvCPSigPropSpeed()
262 , fvCPDelTof()
263 , fvCPTOff()
264 , fvCPTotGain()
265 , fvCPWalk()
266 , fiNbSameSide(0)
267 , fhNbSameSide(NULL)
268 , fhNbDigiPerChan(NULL)
269 , fStart()
270 , fStop()
271 , fTimer()
272 , fiNofEvents(0.)
273 , fdNofHitsTot(0.)
274 , fdTimeTot(0.)
275 , dTRef(0.)
276 , fdTRefMax(0.)
277 , fCalMode(0)
278 , fCalTrg(0)
279 , fCalSmType(0)
280 , fdCaldXdYMax(0.)
281 , fTRefMode(0)
282 , fTRefHits(0)
283 , fPosYMaxScal(0.)
284 , fTRefDifMax(0.)
285 , fTotMax(0.)
286 , fTotMin(0.)
287 , fOutTimeFactor(1.)
288 , fCalParFileName("")
289 , fCalParFile(NULL)
290 , fbMcTrkMonitor(kFALSE)
291
292{
293}
294
296{
297 if (fGeoHandler) delete fGeoHandler;
298 // DeleteHistos(); // <-- if needed ?
299}
300
301/************************************************************************************/
302// FairTasks inherited functions
304{
305
307
308 if (kFALSE == RegisterInputs()) return kFATAL;
309
310 if (kFALSE == RegisterOutputs()) return kFATAL;
311
312 if (kFALSE == InitParameters()) return kFATAL;
313
314 if (kFALSE == LoadGeometry()) return kFATAL;
315
316 if (kFALSE == InitCalibParameter()) return kFATAL;
317
318 if (kFALSE == CreateHistos()) return kFATAL;
319
320 return kSUCCESS;
321}
322
324{
325 LOG(info) << " CbmTofSimpClusterizer => Get the digi parameters for tof";
326
327 // Get Base Container
328 FairRunAna* ana = FairRunAna::Instance();
329 FairRuntimeDb* rtdb = ana->GetRuntimeDb();
330
331 fDigiPar = (CbmTofDigiPar*) (rtdb->getContainer("CbmTofDigiPar"));
332
333 LOG(info) << " CbmTofSimpClusterizer::SetParContainers found " << fDigiPar->GetNrOfModules() << " cells ";
334
335 fDigiBdfPar = (CbmTofDigiBdfPar*) (rtdb->getContainer("CbmTofDigiBdfPar"));
336}
337
338void CbmTofSimpClusterizer::Exec(Option_t* /*option*/)
339{
340 // Start timer counter
341 fTimer.Start();
342
343 fTofHitsColl->Clear("C");
344 // fTofDigiMatchColl->Clear("C"); // Not enough => CbmMatch has no Clear functions!!
345 fTofDigiMatchColl->Delete();
346
347 fiNbHits = 0;
348
349 fStart.Set();
350
351 // --- Local variables
353 Int_t nEvents = 0;
354 Int_t nDigis = 0;
355 Int_t nHits = 0;
356 CbmEvent* event = nullptr;
357 pair<Int_t, Int_t> result;
358
359 // --- Time-slice mode: process entire digi array
360 if (!fEvents) {
361 result = BuildClusters(nullptr);
362 nDigis = result.first;
363 nHits = result.second;
364 }
365
366 // --- Event-based mode: read and process event after event
367 else {
368 nEvents = fEvents->GetEntriesFast();
369 for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
370 event = dynamic_cast<CbmEvent*>(fEvents->At(iEvent));
371 assert(event);
372 result = BuildClusters(event);
373 nDigis += result.first;
374 nHits += result.second;
375 } //# events
376 } //? event mode
377
378
379 fStop.Set();
380
381 FillHistos();
382 fTimer.Stop();
383
384 // --- Timeslice log
385 stringstream logOut;
386 logOut << setw(20) << left << GetName() << " [";
387 logOut << fixed << setw(8) << setprecision(1) << right << fTimer.RealTime() * 1000. << " ms] ";
388 logOut << "TS " << fiNofTs;
389 if (fEvents) logOut << ", events " << nEvents;
390 logOut << ", digis " << nDigis << " / " << nDigisAll;
391 logOut << ", hits " << nHits;
392 LOG(info) << logOut.str();
393
394
395 // --- Update Counters
396 fiNofTs++;
397 fiNofEvents += nEvents;
398 fNofDigisAll += nDigisAll;
399 fNofDigisUsed += nDigis;
400 fdNofHitsTot += nHits;
401 fdTimeTot += fTimer.RealTime();
402}
403
405{
406 WriteHistos();
407 // Prevent them from being sucked in by the CbmHadronAnalysis WriteHistograms method
408 DeleteHistos();
409
410 // Screen log
411 std::cout << std::endl;
412 LOG(info) << "=====================================";
413 LOG(info) << GetName() << ": Run summary";
414 LOG(info) << "Time slices : " << fiNofTs;
415 LOG(info) << "Digis / TS : " << fixed << setprecision(2) << fNofDigisAll / Double_t(fiNofTs);
416 LOG(info) << "Hits / TS : " << fixed << setprecision(2) << fdNofHitsTot / Double_t(fiNofTs);
417 LOG(info) << "Time / TS : " << fixed << setprecision(2) << 1000. * fdTimeTot / Double_t(fiNofTs) << " ms";
418 if (fEvents) {
419 Double_t unusedFrac = 100. * (1. - fNofDigisUsed / fNofDigisAll);
420 LOG(info) << "Digis outside events : " << fNofDigisAll - fNofDigisUsed << " = " << unusedFrac << " %";
421 LOG(info) << "Events : " << fiNofEvents;
422 LOG(info) << "Events / TS : " << fixed << setprecision(2) << Double_t(fiNofEvents) / Double_t(fiNofTs);
423 LOG(info) << "Digis / event : " << fixed << setprecision(2) << fNofDigisUsed / Double_t(fiNofEvents);
424 LOG(info) << "Hits / event : " << fixed << setprecision(2) << fdNofHitsTot / Double_t(fiNofEvents);
425 }
426 LOG(info) << "=====================================\n";
427}
428
429/************************************************************************************/
430// Functions common for all clusters approximations
432{
433 FairRootManager* fManager = FairRootManager::Instance();
434
443 fTofPointsColl = nullptr;
444
445
446 // --- Check input branch (TofDigiExp). If not present, set task inactive.
448 LOG(error) << GetName() << ": No TofDigi input array present; "
449 << "task will be inactive.";
450 return kERROR;
451 }
452
453
454 fMcTracksColl = (TClonesArray*) fManager->GetObject("MCTrack");
455 if (NULL == fMcTracksColl) {
456 LOG(info) << "CbmTofSimpClusterizer: No MCTrack array.";
457 } // if( NULL == fMcTracksColl)
458
459
460 // --- Look for event branch
461 fEvents = dynamic_cast<TClonesArray*>(fManager->GetObject("Event"));
462 if (fEvents)
463 LOG(info) << GetName() << ": Found Event branch; run event-based";
464 else {
465 fEvents = dynamic_cast<TClonesArray*>(fManager->GetObject("CbmEvent"));
466 if (fEvents)
467 LOG(info) << GetName() << ": Found CbmEvent branch; run event-based";
468 else
469 LOG(info) << GetName() << ": No event branch found; run time-based";
470 }
471
472
473 return kTRUE;
474}
476{
477 FairRootManager* rootMgr = FairRootManager::Instance();
478
479 fTofHitsColl = new TClonesArray("CbmTofHit");
480
481 // Flag check to control whether digis are written in output root file
482 rootMgr->Register("TofHit", "Tof", fTofHitsColl, IsOutputBranchPersistent("TofHit"));
483
484 fTofDigiMatchColl = new TClonesArray("CbmMatch", 100);
485 rootMgr->Register("TofHitDigiMatch", "Tof", fTofDigiMatchColl, IsOutputBranchPersistent("TofHitDigiMatch"));
486
487 return kTRUE;
488}
490{
491
492 // Initialize the TOF GeoHandler
493 Bool_t isSimulation = kFALSE;
494 Int_t iGeoVersion = fGeoHandler->Init(isSimulation);
495 LOG(info) << "CbmTofSimpClusterizer::InitParameters with GeoVersion " << iGeoVersion;
496
497 if (k12b > iGeoVersion) {
498 LOG(error) << "CbmTofSimpClusterizer::InitParameters => Only compatible "
499 "with geometries after v12b !!!";
500 return kFALSE;
501 }
502
503 if (NULL != fTofId)
504 LOG(info) << "CbmTofSimpClusterizer::InitParameters with GeoVersion " << fGeoHandler->GetGeoVersion();
505 else {
506 switch (iGeoVersion) {
507 case k12b: fTofId = new CbmTofDetectorId_v12b(); break;
508 case k14a: fTofId = new CbmTofDetectorId_v14a(); break;
509 case k21a: fTofId = new CbmTofDetectorId_v21a(); break;
510 default:
511 LOG(error) << "CbmTofSimpClusterizer::InitParameters => Invalid geometry!!!" << iGeoVersion;
512 return kFALSE;
513 }
514 }
515
516 LOG(info) << "=> Get the digi parameters for tof";
517 FairRunAna* ana = FairRunAna::Instance();
518 FairRuntimeDb* rtdb = ana->GetRuntimeDb();
519
520 // create digitization parameters from geometry file
521 CbmTofCreateDigiPar* tofDigiPar = new CbmTofCreateDigiPar("TOF Digi Producer", "TOF task");
522 LOG(info) << "Create DigiPar ";
523 tofDigiPar->Init();
524
525 fDigiPar = (CbmTofDigiPar*) (rtdb->getContainer("CbmTofDigiPar"));
526 if (0 == fDigiPar) {
527 LOG(error) << "CbmTofSimpleClusterizer::InitParameters => Could not obtain "
528 "the CbmTofDigiPar ";
529 return kFALSE;
530 }
531
533 fdParSystTimeRes = 0.080;
534
535 LOG(info) << " CbmTofSimpClusterizer::InitParameters found Tres FEE = " << fdParFeeTimeRes << " ns ";
536 LOG(info) << " CbmTofSimpClusterizer::InitParameters found Tres Syst = " << fdParSystTimeRes << " ns ";
537
538 return kTRUE;
539}
540/************************************************************************************/
542{
543 // dimension and initialize calib parameter
544 //
545 Int_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes();
546
547 fvCPSigPropSpeed.resize(iNbSmTypes);
548 for (Int_t iT = 0; iT < iNbSmTypes; iT++) {
549 Int_t iNbRpc = fDigiBdfPar->GetNbRpc(iT);
550 fvCPSigPropSpeed[iT].resize(iNbRpc);
551 for (Int_t iRpc = 0; iRpc < iNbRpc; iRpc++)
552 if (0.0 < fDigiBdfPar->GetSigVel(iT, 0, iRpc))
553 fvCPSigPropSpeed[iT][iRpc] = fDigiBdfPar->GetSigVel(iT, 0, iRpc);
554 else
556 } // for (Int_t iT=0; iT<iNbSmTypes; iT++)
557
558 fvCPTOff.resize(iNbSmTypes);
559 fvCPTotGain.resize(iNbSmTypes);
560 fvCPWalk.resize(iNbSmTypes);
561 fvCPDelTof.resize(iNbSmTypes);
562 for (Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
563 Int_t iNbSm = fDigiBdfPar->GetNbSm(iSmType);
564 Int_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
565 fvCPTOff[iSmType].resize(iNbSm * iNbRpc);
566 fvCPTotGain[iSmType].resize(iNbSm * iNbRpc);
567 fvCPWalk[iSmType].resize(iNbSm * iNbRpc);
568 fvCPDelTof[iSmType].resize(iNbSm * iNbRpc);
569 for (Int_t iSm = 0; iSm < iNbSm; iSm++) {
570 for (Int_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
571 // LOG(info)<<Form(" fvCPDelTof resize for SmT %d, R %d, B %d ",iSmType,iNbSm*iNbRpc,nbClDelTofBinX)
572 // ;
573 fvCPDelTof[iSmType][iSm * iNbRpc + iRpc].resize(nbClDelTofBinX);
574 for (Int_t iBx = 0; iBx < nbClDelTofBinX; iBx++) {
575 // LOG(info)<<Form(" fvCPDelTof for SmT %d, R %d, B %d",iSmType,iSm*iNbRpc+iRpc,iBx);
576 fvCPDelTof[iSmType][iSm * iNbRpc + iRpc][iBx].resize(iNTrg);
577 for (Int_t iTrg = 0; iTrg < iNTrg; iTrg++)
578 fvCPDelTof[iSmType][iSm * iNbRpc + iRpc][iBx][iTrg] = 0.; // initialize
579 }
580
581 Int_t iNbChan = fDigiBdfPar->GetNbChan(iSmType, iRpc);
582 fvCPTOff[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
583 fvCPTotGain[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
584 fvCPWalk[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
585 Int_t nbSide = 2 - fDigiBdfPar->GetChanType(iSmType, iRpc);
586 for (Int_t iCh = 0; iCh < iNbChan; iCh++) {
587 fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh].resize(nbSide);
588 fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh].resize(nbSide);
589 fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh].resize(nbSide);
590 for (Int_t iSide = 0; iSide < nbSide; iSide++) {
591 fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][iSide] = 0.; //initialize
592 fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][iSide] = 1.; //initialize
593 fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][iSide].resize(nbClWalkBinX);
594 for (Int_t iWx = 0; iWx < nbClWalkBinX; iWx++) {
595 fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][iSide][iWx] = 0.;
596 }
597 }
598 }
599 }
600 }
601 }
602
603 LOG(info) << "CbmTofSimpClusterizer::InitCalibParameter: defaults set";
604
606 // <= To prevent histos from being sucked in by the param file of the TRootManager!
607 TFile* oldFile = gFile;
608 TDirectory* oldDir = gDirectory;
609
610 if (0 < fCalMode) {
611 LOG(info) << "CbmTofSimpClusterizer::InitCalibParameter: read histos from "
612 << "file " << fCalParFileName;
613
614 // read parameter from histos
615 if (fCalParFileName.IsNull()) return kTRUE;
616
617 fCalParFile = new TFile(fCalParFileName, "");
618 if (NULL == fCalParFile) {
619 LOG(error) << "CbmTofSimpClusterizer::InitCalibParameter: "
620 << "file " << fCalParFileName << " does not exist!";
621 return kTRUE;
622 }
623 /*
624 gDirectory->Print();
625 fCalParFile->cd();
626 fCalParFile->ls();
627 */
628
629 for (Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
630 Int_t iNbSm = fDigiBdfPar->GetNbSm(iSmType);
631 Int_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
632 for (Int_t iSm = 0; iSm < iNbSm; iSm++)
633 for (Int_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
634 TH2F* htempPos_pfx =
635 (TH2F*) gDirectory->FindObjectAny(Form("cl_CorSmT%01d_sm%03d_rpc%03d_Pos_pfx", iSmType, iSm, iRpc));
636 TH2F* htempTOff_pfx =
637 (TH2F*) gDirectory->FindObjectAny(Form("cl_CorSmT%01d_sm%03d_rpc%03d_TOff_pfx", iSmType, iSm, iRpc));
638 TH2F* htempTot_pfx =
639 (TH2F*) gDirectory->FindObjectAny(Form("cl_CorSmT%01d_sm%03d_rpc%03d_Tot_pfx", iSmType, iSm, iRpc));
640 if (NULL != htempPos_pfx && NULL != htempTOff_pfx && NULL != htempTot_pfx) {
641 Int_t iNbCh = fDigiBdfPar->GetNbChan(iSmType, iRpc);
642 Int_t iNbinTot = htempTot_pfx->GetNbinsX();
643 for (Int_t iCh = 0; iCh < iNbCh; iCh++) {
644 Double_t YMean = ((TProfile*) htempPos_pfx)->GetBinContent(iCh + 1); //nh +1 empirical(?)
645 Double_t TMean = ((TProfile*) htempTOff_pfx)->GetBinContent(iCh + 1);
646 Double_t dTYOff = YMean / fvCPSigPropSpeed[iSmType][iRpc];
647 fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][0] += -dTYOff + TMean;
648 fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][1] += +dTYOff + TMean;
649
650 Double_t TotMean = ((TProfile*) htempTot_pfx)->GetBinContent(iCh + 1); //nh +1 empirical(?)
651 if (1 < TotMean) {
652 fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][0] *= TTotMean / TotMean;
653 fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][1] *= TTotMean / TotMean;
654 }
655 LOG(debug1) << "CbmTofSimpClusterizer::InitCalibParameter:"
656 << " SmT " << iSmType << " Sm " << iSm << " Rpc " << iRpc << " Ch " << iCh << ": YMean "
657 << YMean << ", TMean " << TMean << " -> " << fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][0]
658 << ", " << fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][1] << ", "
659 << fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][0] << ", NbinTot " << iNbinTot;
660
661 TH1D* htempWalk0 = (TH1D*) gDirectory->FindObjectAny(
662 Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S0_Walk_px", iSmType, iSm, iRpc, iCh));
663 TH1D* htempWalk1 = (TH1D*) gDirectory->FindObjectAny(
664 Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S1_Walk_px", iSmType, iSm, iRpc, iCh));
665 if (NULL != htempWalk0 && NULL != htempWalk1) { // reinitialize Walk array
666 LOG(info) << "Initialize Walk correction for "
667 << Form(" SmT%01d_sm%03d_rpc%03d_Ch%03d", iSmType, iSm, iRpc, iCh);
668 if (htempWalk0->GetNbinsX() != nbClWalkBinX)
669 LOG(error) << "CbmTofSimpClusterizer::InitCalibParameter: "
670 "Inconsistent Walk histograms";
671 for (Int_t iBin = 0; iBin < nbClWalkBinX; iBin++) {
672 fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][0][iBin] = htempWalk0->GetBinContent(iBin + 1);
673 fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][1][iBin] = htempWalk1->GetBinContent(iBin + 1);
674 LOG(debug1) << Form(" SmT%01d_sm%03d_rpc%03d_Ch%03d bin %d walk %f ", iSmType, iSm, iRpc, iCh, iBin,
675 fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][0][iBin]);
676 }
677 }
678 }
679 }
680 else {
681 LOG(error) << " Histos " << Form("cl_SmT%01d_sm%03d_rpc%03d_XXX", iSmType, iSm, iRpc) << " not found. ";
682 }
683 for (Int_t iTrg = 0; iTrg < iNTrg; iTrg++) {
684 TH1D* htmpDelTof = (TH1D*) gDirectory->FindObjectAny(
685 Form("cl_CorSmT%01d_sm%03d_rpc%03d_Trg%02d_DelTof", iSmType, iSm, iRpc, iTrg));
686 if (NULL == htmpDelTof) {
687 LOG(info) << " Histos " << Form("cl_CorSmT%01d_sm%03d_rpc%03d_Trg%02d_DelTof", iSmType, iSm, iRpc, iTrg)
688 << " not found. ";
689 continue;
690 }
691 LOG(info) << " Load DelTof from histos "
692 << Form("cl_CorSmT%01d_sm%03d_rpc%03d_Trg%02d_DelTof", iSmType, iSm, iRpc, iTrg) << ".";
693 for (Int_t iBx = 0; iBx < nbClDelTofBinX; iBx++) {
694 fvCPDelTof[iSmType][iSm * iNbRpc + iRpc][iBx][iTrg] += htmpDelTof->GetBinContent(iBx + 1);
695 }
696
697 // copy Histo to memory
698 TDirectory* curdir = gDirectory;
699 gDirectory->cd(oldDir->GetPath());
700 TH1D* h1DelTof =
701 (TH1D*) htmpDelTof->Clone(Form("cl_CorSmT%01d_sm%03d_rpc%03d_Trg%02d_DelTof", iSmType, iSm, iRpc, iTrg));
702
703 LOG(info) << " copy histo " << h1DelTof->GetName() << " to directory " << oldDir->GetName();
704
705 gDirectory->cd(curdir->GetPath());
706 }
707 }
708 }
709 }
710 // fCalParFile->Delete();
712 // <= To prevent histos from being sucked in by the param file of the TRootManager!
713 gFile = oldFile;
714 gDirectory = oldDir;
715 LOG(info) << "CbmTofSimpClusterizer::InitCalibParameter: initialization done";
716 return kTRUE;
717}
718/************************************************************************************/
720{
721 LOG(debug) << "CbmTofSimpClusterizer::LoadGeometry starting for " << fDigiPar->GetNrOfModules()
722 << " geometrically known modules ";
723
724 Int_t iNrOfCells = fDigiPar->GetNrOfModules();
725 LOG(debug) << "Digi Parameter container contains " << iNrOfCells << " cells"
726 << ", interpret with GeoVersion " << fGeoHandler->GetGeoVersion();
727
729
730 for (Int_t icell = 0; icell < iNrOfCells; ++icell) {
731
732 Int_t cellId = fDigiPar->GetCellId(icell); // cellId is assigned in CbmTofCreateDigiPar
733 fChannelInfo = fDigiPar->GetCell(cellId);
734
735 Int_t smtype = fGeoHandler->GetSMType(cellId);
736 Int_t smodule = fGeoHandler->GetSModule(cellId);
737 Int_t module = fGeoHandler->GetCounter(cellId);
738 Int_t cell = fGeoHandler->GetCell(cellId);
739
740 Double_t x = fChannelInfo->GetX();
741 Double_t y = fChannelInfo->GetY();
742 Double_t z = fChannelInfo->GetZ();
743 Double_t dx = fChannelInfo->GetSizex();
744 Double_t dy = fChannelInfo->GetSizey();
745 LOG(debug1) << "-I- InitPar " << icell << " Id: " << Form("0x%08x", cellId) << " " << cell << " tmcs: " << smtype
746 << " " << smodule << " " << module << " " << cell << " x=" << Form("%6.2f", x)
747 << " y=" << Form("%6.2f", y) << " z=" << Form("%6.2f", z) << " dx=" << dx << " dy=" << dy;
748 }
749
750 Int_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes();
751
752 if (kTRUE == fDigiBdfPar->UseExpandedDigi()) {
753 fStorDigiExp.resize(iNbSmTypes);
754 fStorDigiInd.resize(iNbSmTypes);
755 fviClusterSize.resize(iNbSmTypes);
756 fviTrkMul.resize(iNbSmTypes);
757 fvdX.resize(iNbSmTypes);
758 fvdY.resize(iNbSmTypes);
759 fvdDifX.resize(iNbSmTypes);
760 fvdDifY.resize(iNbSmTypes);
761 fvdDifCh.resize(iNbSmTypes);
762 fviClusterMul.resize(iNbSmTypes);
763 for (Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
764 Int_t iNbSm = fDigiBdfPar->GetNbSm(iSmType);
765 Int_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
766 fStorDigiExp[iSmType].resize(iNbSm * iNbRpc);
767 fStorDigiInd[iSmType].resize(iNbSm * iNbRpc);
768 fviClusterSize[iSmType].resize(iNbRpc);
769 fviTrkMul[iSmType].resize(iNbRpc);
770 fvdX[iSmType].resize(iNbRpc);
771 fvdY[iSmType].resize(iNbRpc);
772 fvdDifX[iSmType].resize(iNbRpc);
773 fvdDifY[iSmType].resize(iNbRpc);
774 fvdDifCh[iSmType].resize(iNbRpc);
775 for (Int_t iSm = 0; iSm < iNbSm; iSm++) {
776 fviClusterMul[iSmType].resize(iNbSm);
777 for (Int_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
778 fviClusterMul[iSmType][iSm].resize(iNbRpc);
779 Int_t iNbChan = fDigiBdfPar->GetNbChan(iSmType, iRpc);
780 LOG(debug) << "CbmTofSimpClusterizer::LoadGeometry: StoreDigi with "
781 << Form(" %3d %3d %3d %3d %5d ", iSmType, iSm, iNbRpc, iRpc, iNbChan);
782 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
783 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
784 } // for( Int_t iRpc = 0; iRpc < iNbRpc; iRpc++ )
785 } // for( Int_t iSm = 0; iSm < iNbSm; iSm++ )
786 } // for( Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++ )
787 } // if( kTRUE == fDigiBdfPar->UseExpandedDigi() )
788 LOG(debug) << "CbmTofSimpClusterizer::LoadGeometry: Done!";
789 return kTRUE;
790}
792{
793 Int_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes();
794 if (kTRUE == fDigiBdfPar->UseExpandedDigi()) {
795 for (Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
796 Int_t iNbSm = fDigiBdfPar->GetNbSm(iSmType);
797 Int_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
798 for (Int_t iSm = 0; iSm < iNbSm; iSm++) {
799 for (Int_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
800 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc].clear();
801 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc].clear();
802 }
803 } // for( Int_t iSm = 0; iSm < iNbSm; iSm++ )
804 fStorDigiExp[iSmType].clear();
805 fStorDigiInd[iSmType].clear();
806 } // for( Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++ )
807 fStorDigiExp.clear();
808 fStorDigiInd.clear();
809 } // if( kTRUE == fDigiBdfPar->UseExpandedDigi() )
810 return kTRUE;
811}
812/************************************************************************************/
813// Histogramming functions
815{
816 TDirectory* oldir = gDirectory; // <= To prevent histos from being sucked in by the param file of the TRootManager!
817 gROOT->cd(); // <= To prevent histos from being sucked in by the param file of the TRootManager !
819 new TH1I("TofSimpClus_ClustBuildTime", "Time needed to build clusters in each event; Time [s]", 4000, 0.0, 4.0);
820 fhHitsPerTracks = new TH1I("TofSimpClus_TofHitPerTrk",
821 "Mean Number of TofHit per Mc Track; Nb TofHits/Nb MC Tracks []", 2000, 0.0, 20.0);
822 if (kFALSE == fDigiBdfPar->ClustUseTrackId())
823 fhPtsPerHit = new TH1I("TofSimpClus_TofPtsPerHit",
824 "Distribution of the Number of MCPoints associated "
825 "to each TofHit; Nb MCPoint []",
826 20, 0.0, 20.0);
827 if (kTRUE == fDigiBdfPar->ClustUseTrackId()) {
828 fhTimeResSingHits = new TH1I("TofSimpClus_TofTimeResClust",
829 "Time resolution for TofHits containing Digis from a single MC "
830 "Track; t(1st Mc Point) -tTofHit [ns]",
831 10000, -25.0, 25.0);
832 fhTimeResSingHitsB = new TH2I("TofSimpClus_TofTimeResClustB",
833 "Time resolution for TofHits containing Digis from a single MC "
834 "Track; (1st Mc Point) -tTofHit [ns]",
835 5000, -25.0, 25.0, 6, 0, 6);
836 fhTimePtVsHits = new TH2I("TofSimpClus_TofTimePtVsHit",
837 "Time resolution for TofHits containing Digis from a single MC "
838 "Track; t(1st Mc Point) -tTofHit [ps]",
839 2000, 0.0, 50000.0, 2000, 0.0, 50000.0);
840 }
841 else {
842 fhTimeResSingHits = new TH1I("TofSimpClus_TofTimeResClust",
843 "Time resolution for TofHits containing Digis from a single "
844 "TofPoint; tMcPoint -tTofHit [ns]",
845 10000, -25.0, 25.0);
846 fhTimeResSingHitsB = new TH2I("TofSimpClus_TofTimeResClustB",
847 "Time resolution for TofHits containing Digis from a single "
848 "TofPoint; tMcPoint -tTofHit [ns]",
849 5000, -25.0, 25.0, 6, 0, 6);
850 fhTimePtVsHits = new TH2I("TofSimpClus_TofTimePtVsHit",
851 "Time resolution for TofHits containing Digis "
852 "from a single TofPoint; tMcPoint -tTofHit [ps]",
853 2000, 0.0, 50000.0, 2000, 0.0, 50000.0);
854 } // else of if( kTRUE == fDigiBdfPar->ClustUseTrackId() )
856 new TH1I("TofSimpClus_ClusterSize", "Cluster Size distribution; Cluster Size [Strips]", 100, 0.5, 100.5);
858 new TH2I("TofSimpClus_ClusterSizeType",
859 "Cluster Size distribution in each (SM type, Rpc) pair; Cluster "
860 "Size [Strips]; 10*SM Type + Rpc Index []",
861 100, 0.5, 100.5, 40 * fDigiBdfPar->GetNbSmTypes(), 0.0, 40 * fDigiBdfPar->GetNbSmTypes());
862 if (kTRUE == fDigiBdfPar->ClustUseTrackId()) {
863 fhTrackMul = new TH1I("TofSimpClus_TrackMul",
864 "Number of MC tracks generating the cluster; MC Tracks multiplicity []", 100, 0.5, 100.5);
865 fhClusterSizeMulti = new TH2I("TofSimpClus_ClusterSizeMulti",
866 "Cluster Size distribution as function of Number of MC tracks generating "
867 "the cluster; Cluster Size [Strips]; MC tracks mul. []",
868 100, 0.5, 100.5, 100, 0.5, 100.5);
869 fhTrk1MulPos = new TH2D("TofSimpClus_Trk1MulPos",
870 "Position of Clusters with only 1 MC tracks "
871 "generating the cluster; X [cm]; Y [cm]",
872 1500, -750, 750, 1000, -500, 500);
873 fhHiTrkMulPos = new TH2D("TofSimpClus_HiTrkMulPos",
874 "Position of Clusters with >1 MC tracks "
875 "generating the cluster; X [cm]; Y [cm]",
876 1500, -750, 750, 1000, -500, 500);
878 new TH2D("TofSimpClus_AllTrkMulPos", "Position of all clusters generating the cluster; X [cm]; Y [cm]", 1500,
879 -750, 750, 1000, -500, 500);
880 fhMultiTrkProbPos = new TH2D("TofSimpClus_MultiTrkProbPos",
881 "Probability of having a cluster with multiple tracks as "
882 "function of position; X [cm]; Y [cm]; Prob. [%]",
883 1500, -750, 750, 1000, -500, 500);
884 } // if( kTRUE == fDigiBdfPar->ClustUseTrackId() )
885
886 fhDigSpacDifClust = new TH1I("TofSimpClus_DigSpacDifClust",
887 "Space difference along channel direction between Digi pairs on "
888 "adjacent channels; PosCh i - Pos Ch i+1 [cm]",
889 5000, -10.0, 10.0);
890 fhDigTimeDifClust = new TH1I("TofSimpClus_DigTimeDifClust",
891 "Time difference between Digi pairs on adjacent channels; "
892 "0.5*(tDigiA + tDigiA)chi - 0.5*(tDigiA + tDigiA)chi+1 [ns]",
893 5000, -5.0, 5.0);
894 fhDigDistClust = new TH2I("TofSimpClus_DigDistClust",
895 "Distance between Digi pairs on adjacent channels; PosCh i - Pos Ch i+1 "
896 "[cm along ch]; 0.5*(tDigiA + tDigiA)chi - 0.5*(tDigiA + tDigiA)chi+1 [ns]",
897 5000, -10.0, 10.0, 2000, -5.0, 5.0);
898
899 fhClustSizeDifX = new TH2I("TofSimpClus_ClustSizeDifX",
900 "Position difference between Point and Hit as function of Cluster "
901 "Size; Cluster Size [Strips]; dX [cm]",
902 100, 0.5, 100.5, 500, -50, 50);
903 fhClustSizeDifY = new TH2I("TofSimpClus_ClustSizeDifY",
904 "Position difference between Point and Hit as function of Cluster "
905 "Size; Cluster Size [Strips]; dY [cm]",
906 100, 0.5, 100.5, 500, -50, 50);
907 fhChDifDifX = new TH2I("TofSimpClus_ChDifDifX",
908 "Position difference between Point and Hit as "
909 "function of Channel dif; Ch Dif [Strips]; dX [cm]",
910 101, -50.5, 50.5, 500, -50, 50);
911 fhChDifDifY = new TH2I("TofSimpClus_ChDifDifY",
912 "Position difference between Point and Hit as "
913 "function of Channel Dif; Ch Dif [Strips]; dY [cm]",
914 101, -50.5, 50.5, 500, -50, 50);
915
916 fhNbSameSide = new TH1I("TofSimpClus_NbSameSide",
917 "How many time per event the 2 digis on a channel "
918 "were of the same side ; Counts/Event []",
919 500, 0.0, 500.0);
920 fhNbDigiPerChan = new TH1I("TofSimpClus_NbDigiPerChan", "Nb of Digis per channel; Nb Digis []", 100, 0.0, 100.0);
921
922 gDirectory->cd(oldir->GetPath()); // <= To prevent histos from being sucked in by the param file of the TRootManager!
923
924 return kTRUE;
925}
927{
928 fhClustBuildTime->Fill(fStop.GetSec() - fStart.GetSec() + (fStop.GetNanoSec() - fStart.GetNanoSec()) / 1e9);
929 Int_t iNbTofHits = fTofHitsColl->GetEntriesFast();
930
932 Int_t iNbTracks = fMcTracksColl->GetEntriesFast();
933
934 // Trakcs Info
935 Int_t iNbTofTracks = 0;
936 Int_t iNbTofTracksPrim = 0;
937 CbmMCTrack* pMcTrk;
938 for (Int_t iTrkInd = 0; iTrkInd < iNbTracks; iTrkInd++) {
939 pMcTrk = (CbmMCTrack*) fMcTracksColl->At(iTrkInd);
940
941 if (0 < pMcTrk->GetNPoints(ECbmModuleId::kTof)) {
942 iNbTofTracks++;
943 }
944 if (0 < pMcTrk->GetNPoints(ECbmModuleId::kTof) && -1 == pMcTrk->GetMotherId()) iNbTofTracksPrim++;
945
946 } // for(Int_t iTrkInd = 0; iTrkInd < nMcTracks; iTrkInd++)
947
948 if (0 < iNbTofTracks) fhHitsPerTracks->Fill((Double_t) iNbTofHits / (Double_t) iNbTofTracks);
949 } // if( fbMcTrkMonitor )
950
951 CbmTofHit* pHit;
952 for (Int_t iHitInd = 0; iHitInd < iNbTofHits; iHitInd++) {
953 pHit = (CbmTofHit*) fTofHitsColl->At(iHitInd);
954 if (kFALSE == fDigiBdfPar->ClustUseTrackId()) fhPtsPerHit->Fill(pHit->GetFlag());
955 if (1 == pHit->GetFlag()) {
956 // CbmTofPoint* pPt = (CbmTofPoint*)pHit->GetRefId();
957 // Using the SetLinks/GetLinks of the TofHit class seems to conflict
958 // with something in littrack QA
959 // CbmTofPoint* pPt = (CbmTofPoint*)(pHit->GetLinks());
960 /*
961 // Use instead the index => WRONG and not consistent with hit creation
962 // Need to loop on digi match object, check the digi to pnt match object for each, etc....
963 // ====> Just comment this code, if anybody needs it, have to implement proper solution
964 CbmTofPoint* pPt = (CbmTofPoint*)fTofPointsColl->At( pHit->GetRefId() );
965 fhTimePtVsHits->Fill( pPt->GetTime(), pHit->GetTime() );
966 fhTimeResSingHits->Fill( pHit->GetTime() - pPt->GetTime() );
967 fhTimeResSingHitsB->Fill( pHit->GetTime() - pPt->GetTime(),
968 fGeoHandler->GetSMType(pPt->GetDetectorID()) );
969*/
970 } // if( 1 == pHit->GetFlag() )
971 } // for( Int_t iHitInd = 0; iHitInd < iNbTofHits; iHitInd++)
972
973 for (Int_t iSmType = 0; iSmType < fDigiBdfPar->GetNbSmTypes(); iSmType++)
974 for (Int_t iRpc = 0; iRpc < fDigiBdfPar->GetNbRpc(iSmType); iRpc++) {
975 for (UInt_t uCluster = 0; uCluster < fviClusterSize[iSmType][iRpc].size(); uCluster++) {
976 fhClusterSize->Fill(fviClusterSize[iSmType][iRpc][uCluster]);
977 fhClusterSizeType->Fill(fviClusterSize[iSmType][iRpc][uCluster], 40 * iSmType + iRpc);
978 if (kTRUE == fDigiBdfPar->ClustUseTrackId()) {
979 fhTrackMul->Fill(fviTrkMul[iSmType][iRpc][uCluster]);
980 fhClusterSizeMulti->Fill(fviClusterSize[iSmType][iRpc][uCluster], fviTrkMul[iSmType][iRpc][uCluster]);
981 if (1 == fviTrkMul[iSmType][iRpc][uCluster])
982 fhTrk1MulPos->Fill(fvdX[iSmType][iRpc][uCluster], fvdY[iSmType][iRpc][uCluster]);
983 if (1 < fviTrkMul[iSmType][iRpc][uCluster])
984 fhHiTrkMulPos->Fill(fvdX[iSmType][iRpc][uCluster], fvdY[iSmType][iRpc][uCluster]);
985 fhAllTrkMulPos->Fill(fvdX[iSmType][iRpc][uCluster], fvdY[iSmType][iRpc][uCluster]);
986 } // if( kTRUE == fDigiBdfPar->ClustUseTrackId() )
987 if (1 == fviTrkMul[iSmType][iRpc][uCluster]) {
988 fhClustSizeDifX->Fill(fviClusterSize[iSmType][iRpc][uCluster], fvdDifX[iSmType][iRpc][uCluster]);
989 fhClustSizeDifY->Fill(fviClusterSize[iSmType][iRpc][uCluster], fvdDifY[iSmType][iRpc][uCluster]);
990 if (1 == fviClusterSize[iSmType][iRpc][uCluster]) {
991 fhChDifDifX->Fill(fvdDifCh[iSmType][iRpc][uCluster], fvdDifX[iSmType][iRpc][uCluster]);
992 fhChDifDifY->Fill(fvdDifCh[iSmType][iRpc][uCluster], fvdDifY[iSmType][iRpc][uCluster]);
993 }
994 }
995 } // for( UInt_t uCluster = 0; uCluster < fviClusterSize[iSmType][iRpc].size(); uCluster++ )
996 fviClusterSize[iSmType][iRpc].clear();
997 fviTrkMul[iSmType][iRpc].clear();
998 fvdX[iSmType][iRpc].clear();
999 fvdY[iSmType][iRpc].clear();
1000 fvdDifX[iSmType][iRpc].clear();
1001 fvdDifY[iSmType][iRpc].clear();
1002 fvdDifCh[iSmType][iRpc].clear();
1003 } // for( Int_t iRpc = 0; iRpc < fDigiBdfPar->GetNbRpc( iSmType); iRpc++ )
1004
1006
1007 return kTRUE;
1008}
1010{
1011 if ("" == fsHistoOutFilename) {
1012 LOG(info) << "CbmTofSimpClusterizer::WriteHistos => Control histograms "
1013 "will not be written to disk!";
1014 LOG(info) << "CbmTofSimpClusterizer::WriteHistos => To change this "
1015 "behavior please provide a full path "
1016 << "with the SetHistoFileName method";
1017 return kTRUE;
1018 } // if( "" == fsHistoOutFilename )
1019
1021 TFile* oldFile = gFile;
1022 TDirectory* oldDir = gDirectory;
1023
1024 TFile* fHist = new TFile(fsHistoOutFilename, "RECREATE");
1025
1026 fHist->cd();
1027 fhClustBuildTime->Write();
1028 fhHitsPerTracks->Write();
1029 if (kFALSE == fDigiBdfPar->ClustUseTrackId()) fhPtsPerHit->Write();
1030 fhTimeResSingHits->Write();
1031 fhTimeResSingHitsB->Write();
1032 fhTimePtVsHits->Write();
1033 fhClusterSize->Write();
1034 fhClusterSizeType->Write();
1035 if (kTRUE == fDigiBdfPar->ClustUseTrackId()) {
1036 fhTrackMul->Write();
1037 fhClusterSizeMulti->Write();
1038 fhTrk1MulPos->Write();
1039 fhHiTrkMulPos->Write();
1040 fhAllTrkMulPos->Write();
1042 fhMultiTrkProbPos->Scale(100.0);
1043 fhMultiTrkProbPos->Write();
1044 } // if( kTRUE == fDigiBdfPar->ClustUseTrackId() )
1045 fhDigSpacDifClust->Write();
1046 fhDigTimeDifClust->Write();
1047 fhDigDistClust->Write();
1048
1049 fhClustSizeDifX->Write();
1050 fhClustSizeDifY->Write();
1051
1052 fhChDifDifX->Write();
1053 fhChDifDifY->Write();
1054
1055 fhNbSameSide->Write();
1056 fhNbDigiPerChan->Write();
1057
1059 gFile = oldFile;
1060 gDirectory = oldDir;
1061
1062 fHist->Close();
1063
1064 return kTRUE;
1065}
1067{
1068 fsHistoOutFilename = sFilenameIn;
1069 return kTRUE;
1070}
1072{
1073 delete fhClustBuildTime;
1074 delete fhHitsPerTracks;
1075 delete fhPtsPerHit;
1076 delete fhTimeResSingHits;
1077 delete fhTimeResSingHitsB;
1078 delete fhTimePtVsHits;
1079 delete fhClusterSize;
1080 delete fhClusterSizeType;
1081
1082 if (kTRUE == fDigiBdfPar->ClustUseTrackId()) {
1083 delete fhTrackMul;
1084 delete fhClusterSizeMulti;
1085 delete fhTrk1MulPos;
1086 delete fhHiTrkMulPos;
1087 delete fhAllTrkMulPos;
1088 delete fhMultiTrkProbPos;
1089 }
1090 delete fhDigSpacDifClust;
1091 delete fhDigTimeDifClust;
1092 delete fhDigDistClust;
1093
1094 delete fhClustSizeDifX;
1095 delete fhClustSizeDifY;
1096
1097 delete fhChDifDifX;
1098 delete fhChDifDifY;
1099
1100 delete fhNbSameSide;
1101 delete fhNbDigiPerChan;
1102
1103 return kTRUE;
1104}
1105/************************************************************************************/
1107{
1108
1109 // --- MC Event info (input file, entry number, start time)
1110 Int_t iInputNr = 0;
1111 Int_t iEventNr = 0;
1112 Double_t dEventTime = 0.;
1113 GetEventInfo(iInputNr, iEventNr, dEventTime);
1114
1115 // Local variables
1116 Int_t nDigis = 0;
1117 Int_t nHits = 0;
1118 Int_t hitIndex = 0;
1119
1120 /*
1121 * FIXME: maybe use the 2D distance (X/Y) as criteria instead of the distance long channel
1122 * direction
1123 */
1124 Double_t dMaxTimeDist = fDigiBdfPar->GetMaxTimeDist();
1125 Double_t dMaxSpaceDist = fDigiBdfPar->GetMaxDistAlongCh();
1126
1127 LOG(debug) << " BuildCluster with MaxTimeDist " << dMaxTimeDist << ", MaxSpaceDist " << dMaxSpaceDist;
1128
1129 // First Time order the Digis array
1130 // fTofDigisColl->Sort();
1131
1132 // Then loop over the digis array and store the Digis in separate vectors for
1133 // each RPC modules
1134
1135 nDigis = (event ? event->GetNofData(ECbmDataType::kTofDigi) : fDigiMan->GetNofDigis(ECbmModuleId::kTof));
1136 LOG(debug) << "Number of TOF digis: " << nDigis;
1137 if (kTRUE == fDigiBdfPar->UseExpandedDigi()) {
1138 CbmTofDigi* pDigi = nullptr;
1139 UInt_t digiIndex = -1;
1140 for (Int_t iDigi = 0; iDigi < nDigis; iDigi++) {
1141 digiIndex = (event ? event->GetIndex(ECbmDataType::kTofDigi, iDigi) : iDigi);
1142 assert(fDigiMan->Get<CbmTofDigi>(digiIndex));
1143 pDigi = new CbmTofDigi(*(fDigiMan->Get<CbmTofDigi>(digiIndex)));
1144 assert(pDigi);
1145 LOG(debug1) << "CbmTofSimpClusterizer::BuildClusters: " << digiIndex << " " << pDigi << " " << pDigi->GetType()
1146 << " " << pDigi->GetSm() << " " << pDigi->GetRpc() << " " << pDigi->GetChannel() << " "
1147 << Form("T %6.2f, Tot %6.1f ", pDigi->GetTime(), pDigi->GetTot());
1148 if (fDigiBdfPar->GetNbSmTypes() > pDigi->GetType() // prevent crash due to misconfiguration
1149 && fDigiBdfPar->GetNbSm(pDigi->GetType()) > pDigi->GetSm()
1150 && fDigiBdfPar->GetNbRpc(pDigi->GetType()) > pDigi->GetRpc()
1151 && fDigiBdfPar->GetNbChan(pDigi->GetType(), 0) > pDigi->GetChannel()) {
1152
1153 fStorDigiExp[pDigi->GetType()][pDigi->GetSm() * fDigiBdfPar->GetNbRpc(pDigi->GetType()) + pDigi->GetRpc()]
1154 [pDigi->GetChannel()]
1155 .push_back(pDigi);
1156 fStorDigiInd[pDigi->GetType()][pDigi->GetSm() * fDigiBdfPar->GetNbRpc(pDigi->GetType()) + pDigi->GetRpc()]
1157 [pDigi->GetChannel()]
1158 .push_back(digiIndex);
1159
1160 // apply calibration vectors
1161 pDigi->SetTime(pDigi->GetTime() - // calibrate Digi Time
1162 fvCPTOff[pDigi->GetType()][pDigi->GetSm() * fDigiBdfPar->GetNbRpc(pDigi->GetType())
1163 + pDigi->GetRpc()][pDigi->GetChannel()][pDigi->GetSide()]);
1164
1165 pDigi->SetTot(pDigi->GetTot() * // calibrate Digi ToT
1166 fvCPTotGain[pDigi->GetType()][pDigi->GetSm() * fDigiBdfPar->GetNbRpc(pDigi->GetType())
1167 + pDigi->GetRpc()][pDigi->GetChannel()][pDigi->GetSide()]);
1168
1169 // walk correction
1170 Double_t dTotBinSize = (TOTMax - TOTMin) / 2. / nbClWalkBinX;
1171 Int_t iWx = (Int_t)((pDigi->GetTot() - TOTMin / 2.) / dTotBinSize);
1172 if (0 > iWx) iWx = 0;
1173 if (iWx > (nbClWalkBinX - 1)) iWx = nbClWalkBinX - 1;
1174 Double_t dDTot = (pDigi->GetTot() - TOTMin / 2.) / dTotBinSize - (Double_t) iWx - 0.5;
1175 Double_t dWT = fvCPWalk[pDigi->GetType()][pDigi->GetSm() * fDigiBdfPar->GetNbRpc(pDigi->GetType())
1176 + pDigi->GetRpc()][pDigi->GetChannel()][pDigi->GetSide()][iWx];
1177 if (dDTot > 0) {
1178 if (iWx < nbClWalkBinX - 1) { // linear interpolation to next bin
1179 dWT += dDTot
1180 * (fvCPWalk[pDigi->GetType()][pDigi->GetSm() * fDigiBdfPar->GetNbRpc(pDigi->GetType())
1181 + pDigi->GetRpc()][pDigi->GetChannel()][pDigi->GetSide()][iWx + 1]
1182 - fvCPWalk[pDigi->GetType()][pDigi->GetSm() * fDigiBdfPar->GetNbRpc(pDigi->GetType())
1183 + pDigi->GetRpc()][pDigi->GetChannel()][pDigi->GetSide()][iWx]);
1184 } // if(iWx < nbClWalkBinX -1)
1185 }
1186 else // dDTot < 0,
1187 {
1188 if (0 < iWx) { // linear interpolation to next bin
1189 dWT -= dDTot
1190 * (fvCPWalk[pDigi->GetType()][pDigi->GetSm() * fDigiBdfPar->GetNbRpc(pDigi->GetType())
1191 + pDigi->GetRpc()][pDigi->GetChannel()][pDigi->GetSide()][iWx - 1]
1192 - fvCPWalk[pDigi->GetType()][pDigi->GetSm() * fDigiBdfPar->GetNbRpc(pDigi->GetType())
1193 + pDigi->GetRpc()][pDigi->GetChannel()][pDigi->GetSide()][iWx]);
1194 } // if(0 < iWx)
1195 }
1196
1197 pDigi->SetTime(pDigi->GetTime() - dWT); // calibrate Digi Time
1198
1199 if (0) { //pDigi->GetType()==2 && pDigi->GetSm()==0){
1200 LOG(info) << "CbmTofSimpClusterizer::BuildClusters: CalDigi " << digiIndex << ", T " << pDigi->GetType()
1201 << ", Sm " << pDigi->GetSm() << ", R " << pDigi->GetRpc() << ", Ch " << pDigi->GetChannel()
1202 << ", S " << pDigi->GetSide() << ", T " << pDigi->GetTime() << ", Tot " << pDigi->GetTot()
1203 << ", TotGain "
1204 << fvCPTotGain[pDigi->GetType()][pDigi->GetSm() * fDigiBdfPar->GetNbRpc(pDigi->GetType())
1205 + pDigi->GetRpc()][pDigi->GetChannel()][pDigi->GetSide()]
1206 << ", Walk " << iWx << ": "
1207 << fvCPWalk[pDigi->GetType()][pDigi->GetSm() * fDigiBdfPar->GetNbRpc(pDigi->GetType())
1208 + pDigi->GetRpc()][pDigi->GetChannel()][pDigi->GetSide()][iWx];
1209
1210 LOG(info) << " dDTot " << dDTot << " BinSize: " << dTotBinSize << ", CPWalk "
1211 << fvCPWalk[pDigi->GetType()][pDigi->GetSm() * fDigiBdfPar->GetNbRpc(pDigi->GetType())
1212 + pDigi->GetRpc()][pDigi->GetChannel()][pDigi->GetSide()][iWx - 1]
1213 << ", "
1214 << fvCPWalk[pDigi->GetType()][pDigi->GetSm() * fDigiBdfPar->GetNbRpc(pDigi->GetType())
1215 + pDigi->GetRpc()][pDigi->GetChannel()][pDigi->GetSide()][iWx]
1216 << ", "
1217 << fvCPWalk[pDigi->GetType()][pDigi->GetSm() * fDigiBdfPar->GetNbRpc(pDigi->GetType())
1218 + pDigi->GetRpc()][pDigi->GetChannel()][pDigi->GetSide()][iWx + 1]
1219 << " -> dWT = " << dWT;
1220 }
1221 }
1222 else {
1223 LOG(debug) << "CbmTofSimpBeamClusterizer::BuildClusters: Skip Digi "
1224 << " Type " << pDigi->GetType() << " " << fDigiBdfPar->GetNbSmTypes() << " Sm " << pDigi->GetSm()
1225 << " " << fDigiBdfPar->GetNbSm(pDigi->GetType()) << " Rpc " << pDigi->GetRpc() << " "
1226 << fDigiBdfPar->GetNbRpc(pDigi->GetType()) << " Ch " << pDigi->GetChannel() << " "
1227 << fDigiBdfPar->GetNbChan(pDigi->GetType(), 0);
1228 }
1229 } // for( Int_t iDigInd = 0; iDigInd < nTofDigi; iDigInd++ )
1230 } // if( kTRUE == fDigiBdfPar->UseExpandedDigi() )
1231
1232 // Then build clusters inside each RPC module
1233 // Assume only 0 or 1 Digi per channel/side in each event
1234 // Use simplest method possible, scan direction independent:
1235 // a) Loop over channels in the RPC starting from 0
1236 // * If strips
1237 // i) Loop over Digis to check if both ends of the channel have a Digi
1238 // ii) Reconstruct a mean channel time and a mean position
1239 // + If a Hit is currently filled & the mean position (space, time) is less than XXX from last channel position
1240 // iii) Add the mean channel time and the mean position to the ones of the hit
1241 // + else
1242 // iii) Use nb of strips in cluster to cal. the hit mean time and pos (charge/tot weighting)
1243 // iv) Save the hit
1244 // v) Start a new hit with current channel
1245 // * else (pads)
1246 // i) Loop over Digis to find if this channel fired
1247 // ii) FIXME: either scan all other channels to check for matching Digis or have more than 1 hit open
1248 Int_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes();
1249 // Hit variables
1250 Double_t dWeightedTime = 0.0;
1251 Double_t dWeightedPosX = 0.0;
1252 Double_t dWeightedPosY = 0.0;
1253 Double_t dWeightedPosZ = 0.0;
1254 Double_t dWeightedTimeErr = 0.0;
1255 Double_t dWeightsSum = 0.0;
1256 std::vector<CbmTofPoint*> vPtsRef;
1257 std::vector<Int_t> vDigiIndRef;
1258 CbmTofCell* fTrafoCell = NULL;
1259 Int_t iTrafoCell = -1;
1260 Int_t iNbChanInHit = 0;
1261 // Last Channel Temp variables
1262 Int_t iLastChan = -1;
1263 // Double_t dLastPosX = 0.0; // -> Comment to remove warning because set but never used
1264 Double_t dLastPosY = 0.0;
1265 Double_t dLastTime = 0.0;
1266 // Channel Temp variables
1267 Double_t dPosX = 0.0;
1268 Double_t dPosY = 0.0;
1269 Double_t dPosZ = 0.0;
1270 Double_t dTime = 0.0;
1271 Double_t dTotS = 0.0;
1272 fiNbSameSide = 0;
1275 gGeoManager->CdTop();
1276 if (kTRUE == fDigiBdfPar->UseExpandedDigi()) {
1277 for (Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
1278 Int_t iNbSm = fDigiBdfPar->GetNbSm(iSmType);
1279 Int_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
1280 for (Int_t iSm = 0; iSm < iNbSm; iSm++)
1281 for (Int_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
1282 Int_t iNbCh = fDigiBdfPar->GetNbChan(iSmType, iRpc);
1283 Int_t iChType = fDigiBdfPar->GetChanType(iSmType, iRpc);
1284 LOG(debug2) << "CbmTofSimpClusterizer::BuildClusters: RPC - Loop "
1285 << Form(" %3d %3d %3d %3d ", iSmType, iSm, iRpc, iChType);
1286 fviClusterMul[iSmType][iSm][iRpc] = 0;
1287 Int_t iChId = 0;
1288 if (0 == iChType) {
1289 // Don't spread clusters over RPCs!!!
1290 dWeightedTime = 0.0;
1291 dWeightedPosX = 0.0;
1292 dWeightedPosY = 0.0;
1293 dWeightedPosZ = 0.0;
1294 dWeightedTimeErr = 0.0;
1295 dWeightsSum = 0.0;
1296 iNbChanInHit = 0;
1297 vPtsRef.clear();
1298 // For safety reinitialize everything
1299 iLastChan = -1;
1300 // dLastPosX = 0.0; // -> Comment to remove warning because set but never used
1301 dLastPosY = 0.0;
1302 dLastTime = 0.0;
1303 fTrafoCell = NULL;
1304 iTrafoCell = -1;
1305 LOG(debug2) << "CbmTofSimpClusterizer::BuildClusters: ChanOrient "
1306 << Form(" %3d %3d %3d %3d %3d ", iSmType, iSm, iRpc, fDigiBdfPar->GetChanOrient(iSmType, iRpc),
1307 iNbCh);
1308
1309 if (1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc)) {
1310 // Horizontal strips => X comes from left right time difference
1311 } // if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
1312 else {
1313 // Vertical strips => Y comes from bottom top time difference
1314 for (Int_t iCh = 0; iCh < iNbCh; iCh++) {
1315 LOG(debug3) << "CbmTofSimpClusterizer::BuildClusters: VDigisize"
1316 << Form(" T %3d Sm %3d R %3d Ch %3d Size %3zu ", iSmType, iSm, iRpc, iCh,
1317 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].size());
1318
1319 if (0 < fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].size())
1320 fhNbDigiPerChan->Fill(fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].size());
1321 while (1 < fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].size()) {
1322 LOG(debug2) << "CbmTofSimpClusterizer::BuildClusters: digis processing "
1323 "for "
1324 << Form(" SmT %3d Sm %3d Rpc %3d Ch %3d # %3zu ", iSmType, iSm, iRpc, iCh,
1325 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].size());
1326 /*
1327 if (3 < fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size()) // needs more work!
1328 LOG(debug) << "CbmTofSimpClusterizer::BuildClusters: digis ignored for "
1329 << Form(" SmT %3d Sm %3d Rpc %3d Ch %3d # %3d ",iSmType,iSm,iRpc,iCh,fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size())
1330 ;
1331 */
1332
1333 while ((fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh][0])->GetSide()
1334 == (fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh][1])->GetSide()) {
1335 // Not one Digi of each end!
1336 fiNbSameSide++;
1337 LOG(debug) << "CbmTofSimpClusterizer::BuildClusters: SameSide Hits! "
1338 "Times: "
1339 << (fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh][0])->GetTime() << ", "
1340 << (fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh][1])->GetTime() << ", DeltaT "
1341 << (fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh][1])->GetTime()
1342 - (fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh][0])->GetTime();
1343 LOG(debug2) << " SameSide Erase fStor entries(d) " << iSmType << ", SR " << iSm * iNbRpc + iRpc
1344 << ", Ch" << iCh;
1345 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1346 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
1347 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1348 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
1349 if (2 > fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].size()) break;
1350 continue;
1351 }
1352
1353 LOG(debug2) << "CbmTofSimpClusterizer::BuildClusters: digis processing "
1354 "for "
1355 << Form(" SmT %3d Sm %3d Rpc %3d Ch %3d # %3zu ", iSmType, iSm, iRpc, iCh,
1356 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].size());
1357 if (2 > fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].size()) break;
1358 /* Int_t iLastChId = iChId; // Save Last hit channel */
1359
1360 // 2 Digis = both sides present
1361 Int_t iCh1 = iCh;
1362 if (fGeoHandler->GetGeoVersion() < k14a) iCh1 = iCh1 + 1; //FIXME
1363 CbmTofDetectorInfo xDetInfo(ECbmModuleId::kTof, iSmType, iSm, iRpc, 0, iCh1);
1364
1365 iChId = fTofId->SetDetectorInfo(xDetInfo);
1366 Int_t iUCellId = CbmTofAddress::GetUniqueAddress(iSm, iRpc, iCh, 0, iSmType);
1367
1368 fChannelInfo = fDigiPar->GetCell(iChId);
1369 if (NULL == fChannelInfo) {
1370 LOG(error) << "CbmTofSimpClusterizer::BuildClusters: no "
1371 "geometry info! "
1372 << Form(" %3d %3d %3d %3d 0x%08x 0x%08x ", iSmType, iSm, iRpc, iCh, iChId, iUCellId);
1373 break;
1374 }
1375 LOG(debug1) << "CbmTofSimpClusterizer::BuildClusters:"
1376 << Form(" T %3d Sm %3d R %3d Ch %3d size %3zu ", iSmType, iSm, iRpc, iCh,
1377 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].size())
1378 << Form(" ChId: 0x%08x 0x%08x %p", iChId, iUCellId, fChannelInfo);
1379
1380 CbmTofDigi* xDigiA = fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh][0];
1381 CbmTofDigi* xDigiB = fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh][1];
1382
1383
1384 // The "Strip" time is the mean time between each end
1385 dTime = 0.5 * (xDigiA->GetTime() + xDigiB->GetTime());
1386
1387 // Weight is the total charge => sum of both ends ToT
1388 dTotS = xDigiA->GetTot() + xDigiB->GetTot();
1389
1390 // X position is just the center of the channel
1391 //dPosX = fChannelInfo->GetX();
1392
1393 // Y position from strip ends time difference
1394 //dPosY = fChannelInfo->GetY();
1395
1396 // For Z always just take the one of the channel itself( in fact its gap one)
1397 //dPosZ = fChannelInfo->GetZ();
1398
1399 TGeoNode* fNode = // prepare local->global trafo
1400 gGeoManager->FindNode(fChannelInfo->GetX(), fChannelInfo->GetY(), fChannelInfo->GetZ());
1401 if (NULL == fTrafoCell) {
1402 fTrafoCell = fChannelInfo;
1403 iTrafoCell = iCh;
1404 }
1405
1406 //fNode->Print();
1407 LOG(debug1) << Form(" Node at (%6.1f,%6.1f,%6.1f) : %p, info %p, %p", fChannelInfo->GetX(),
1408 fChannelInfo->GetY(), fChannelInfo->GetZ(), fNode, fChannelInfo, fTrafoCell);
1409 // fNode->Print();
1410
1411
1412 // switch to local coordinates, (0,0,0) is in the center of counter ?
1413 dPosX = ((Double_t)(-iTrafoCell + iCh)) * fChannelInfo->GetSizex();
1414 dPosY = 0.;
1415 dPosZ = 0.;
1416
1417
1418 if (1 == xDigiA->GetSide())
1419
1420 // 0 is the top side, 1 is the bottom side
1421 dPosY += fvCPSigPropSpeed[iSmType][iRpc] * (xDigiA->GetTime() - xDigiB->GetTime()) / 2.0;
1422
1423 else
1424
1425 // 0 is the bottom side, 1 is the top side
1426 dPosY += fvCPSigPropSpeed[iSmType][iRpc] * (xDigiB->GetTime() - xDigiA->GetTime()) / 2.0;
1427
1428 if (fChannelInfo->GetSizey() / 2.0 < dPosY || -1 * fChannelInfo->GetSizey() / 2.0 > dPosY) {
1429 // if outside of strip limits, the pair is bad => try to remove one end and check the next pair
1430 // (if another possibility exist)
1431 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1432 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
1433 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1434 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
1435 continue;
1436 } // Pair leads to hit oustide of strip limits
1437
1438 LOG(debug1) << "CbmTofSimpClusterizer::BuildClusters: NbChanInHit"
1439 << Form(" %3d %3d %3d %p %2.0f Time %6.1f PosY %5.1f Svel "
1440 "%5.1f",
1441 iNbChanInHit, iCh, iLastChan, xDigiA, xDigiA->GetSide(), dTime, dPosY,
1442 fvCPSigPropSpeed[iSmType][iRpc])
1443 << Form(", Offs %5.1f, %5.1f", fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][0],
1444 fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][1]);
1445
1446 // Now check if a hit/cluster is already started
1447 if (0 < iNbChanInHit) {
1448
1449
1450 if (iLastChan == iCh - 1) {
1451 fhDigTimeDifClust->Fill(dTime - dLastTime);
1452 fhDigSpacDifClust->Fill(dPosY - dLastPosY);
1453 fhDigDistClust->Fill(dPosY - dLastPosY, dTime - dLastTime);
1454 } // if( iLastChan == iCh - 1 )
1455 // a cluster is already started => check distance in space/time
1456 // For simplicity, just check along strip direction for now
1457 // and break cluster when a not fired strip is found
1458 if (TMath::Abs(dTime - dLastTime) < dMaxTimeDist && iLastChan == iCh - 1
1459 && TMath::Abs(dPosY - dLastPosY) < dMaxSpaceDist) {
1460 // Add to cluster/hit
1461 dWeightedTime += dTime * dTotS;
1462 dWeightedPosX += dPosX * dTotS;
1463 dWeightedPosY += dPosY * dTotS;
1464 dWeightedPosZ += dPosZ * dTotS;
1465 dWeightedTimeErr += dTotS * dTotS;
1466 dWeightsSum += dTotS;
1467 iNbChanInHit += 1;
1468 // if one of the side digi comes from a CbmTofPoint not already found
1469 // in this cluster, save its pointer
1470 // Bool_t bFoundA = kFALSE; // -> Comment to remove warning because set but never used
1471 // Bool_t bFoundB = kFALSE; // -> Comment to remove warning because set but never used
1472 if (NULL != fTofPointsColl) { // MC
1473
1475 /*
1476 if( kTRUE == fDigiBdfPar->ClustUseTrackId() )
1477 for( UInt_t uPtRef = 0; uPtRef < vPtsRef.size(); uPtRef++)
1478 {
1479 //if( ((CbmTofPoint*)(vPtsRef[uPtRef]))->GetTrackID() == ((CbmTofPoint*)(xDigiA->GetLinks()))->GetTrackID() )
1480// bFoundA = kTRUE; // -> Comment to remove warning because set but never used
1481 //if( ((CbmTofPoint*)(vPtsRef[uPtRef]))->GetTrackID() == ((CbmTofPoint*)(xDigiB->GetLinks()))->GetTrackID() )
1482// bFoundB = kTRUE; // -> Comment to remove warning because set but never used
1483 } // for( Int
1484 else for( UInt_t uPtRef = 0; uPtRef < vPtsRef.size(); uPtRef++)
1485 {
1486 // if( vPtsRef[uPtRef] == (CbmTofPoint*)xDigiA->GetLinks() )
1487// bFoundA = kTRUE; // -> Comment to remove warning because set but never used
1488 // if( vPtsRef[uPtRef] == (CbmTofPoint*)xDigiB->GetLinks() )
1489// bFoundB = kTRUE; // -> Comment to remove warning because set but never used
1490 } // for( Int_t uPtRef = 0; uPtRef < vPtsRef.size(); uPtRef++)
1491*/
1492
1493 // CbmTofPoint pointer for 1st digi not found => save it
1494 //if( kFALSE == bFoundA )
1495 // vPtsRef.push_back( (CbmTofPoint*)(xDigiA->GetLinks()) );
1496 // CbmTofPoint pointer for 2nd digi not found => save it
1497 // if( kFALSE == bFoundB )
1498 // vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) );
1499 } // MC end
1500 vDigiIndRef.push_back((Int_t)(fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][0]));
1501 vDigiIndRef.push_back((Int_t)(fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][1]));
1502
1503 LOG(debug1) << " Add Digi and erase fStor entries(a): NbChanInHit " << iNbChanInHit << ", "
1504 << iSmType << ", SR " << iSm * iNbRpc + iRpc << ", Ch" << iCh;
1505
1506 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1507 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
1508 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1509 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
1510 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1511 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
1512 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1513 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
1514
1515 } // if current Digis compatible with last fired chan
1516 else {
1517 // Save Hit
1518 dWeightedTime /= dWeightsSum;
1519 dWeightedPosX /= dWeightsSum;
1520 dWeightedPosY /= dWeightsSum;
1521 dWeightedPosZ /= dWeightsSum;
1522 // simple error scaling with TOT
1523 // dWeightedTimeErr = TMath::Sqrt( dWeightedTimeErr ) * fdParSystTimeRes / dWeightsSum;
1524 // OR harcoded value
1525 dWeightedTimeErr = fdParSystTimeRes;
1526
1527 Double_t hitpos_local[3];
1528 Double_t hitpos[3];
1529
1530 TGeoNode* tNode = // prepare local->global trafo
1531 gGeoManager->FindNode(fTrafoCell->GetX(), fTrafoCell->GetY(), fTrafoCell->GetZ());
1532 TGeoNode* cNode = gGeoManager->GetCurrentNode();
1533
1534 TGeoRotation rotMatrix(*gGeoManager->GetCurrentMatrix());
1535
1536 hitpos[0] = fTrafoCell->GetX();
1537 hitpos[1] = fTrafoCell->GetY();
1538 hitpos[2] = fTrafoCell->GetZ();
1539 gGeoManager->MasterToLocal(hitpos, hitpos_local);
1540 LOG(debug1) << Form(" Node0 at (%6.1f,%6.1f,%6.1f) : "
1541 "(%6.1f,%6.1f,%6.1f) pointer %p",
1543 hitpos_local[0], hitpos_local[1], hitpos_local[2], tNode);
1544
1545 hitpos_local[0] += dWeightedPosX;
1546 hitpos_local[1] += dWeightedPosY;
1547 hitpos_local[2] += dWeightedPosZ;
1548
1549 gGeoManager->LocalToMaster(hitpos_local, hitpos);
1550 LOG(debug1) << Form(" LTM for node %p, info %p: "
1551 "(%6.1f,%6.1f,%6.1f) ->(%6.1f,%6.1f,%6.1f)",
1552 cNode, fChannelInfo, hitpos_local[0], hitpos_local[1], hitpos_local[2],
1553 hitpos[0], hitpos[1], hitpos[2]);
1554
1555 TVector3 hitPos(hitpos[0], hitpos[1], hitpos[2]);
1556
1557 // Simple errors, not properly done at all for now
1558 // Right way of doing it should take into account the weight distribution
1559 // and real system time resolution
1560 Double_t hiterr_local[3];
1561 Double_t hiterr[3];
1562
1563 hiterr_local[0] = fChannelInfo->GetSizex() / TMath::Sqrt(12.0); // Single strips approximation
1564 hiterr_local[1] =
1565 fdParFeeTimeRes * fvCPSigPropSpeed[iSmType][iRpc]; // Use the electronics resolution
1566 hiterr_local[2] = fDigiBdfPar->GetNbGaps(iSmType, iRpc) * fDigiBdfPar->GetGapSize(iSmType, iRpc)
1567 / 10.0 // Change gap size in cm
1568 / TMath::Sqrt(12.0); // Use full RPC thickness as "Channel" Z size
1569
1570 rotMatrix.LocalToMaster(hiterr_local, hiterr);
1571 // TVector3 hitPosErr( hiterr_local[0], hiterr_local[1], hiterr_local[2] );
1572 TVector3 hitPosErr(fabs(hiterr[0]), fabs(hiterr[1]), fabs(hiterr[2]));
1573
1574 // Int_t iDetId = vPtsRef[0]->GetDetectorID();// detID = pt->GetDetectorID() <= from TofPoint
1575 // calc mean ch from dPosX=((Double_t)(-iNbCh/2 + iCh)+0.5)*fChannelInfo->GetSizex();
1576
1577 Int_t iChm = iTrafoCell + floor(dWeightedPosX / fChannelInfo->GetSizex());
1578 if (iChm < 0 || iChm > iNbCh) {
1579 LOG(debug) << "CbmTofSimpClusterizer::BuildClusters: "
1580 "Invalid mean channel";
1581 }
1582 Int_t iDetId = CbmTofAddress::GetUniqueAddress(iSm, iRpc, iChm, 0, iSmType);
1583 Int_t iRefId = 0; // Index of the correspondng TofPoint
1584 // if(NULL != fTofPointsColl) { iRefId = fTofPointsColl->IndexOf( vPtsRef[0] ); }
1585 LOG(debug) << "CbmTofSimpClusterizer::BuildClusters: Save Hit "
1586 << Form(" %3d %3d 0x%08x %3d %3d %3d %f %f", fiNbHits, iNbChanInHit, iDetId, iCh,
1587 iLastChan, iRefId, dWeightedTime, dWeightedPosY)
1588 << ", DigiSize: " << vDigiIndRef.size() << ", DigiInds: ";
1589 fviClusterMul[iSmType][iSm][iRpc]++;
1590 for (UInt_t i = 0; i < vDigiIndRef.size(); i++) {
1591 LOG(debug) << " " << vDigiIndRef.at(i);
1592 }
1593 LOG(debug);
1594
1595 hitIndex = fTofHitsColl->GetEntriesFast();
1596 new ((*fTofHitsColl)[hitIndex])
1597 CbmTofHit(iDetId, hitPos,
1598 hitPosErr, //local detector coordinates
1599 fiNbHits, dWeightedTime * fOutTimeFactor, dWeightedTimeErr,
1600 vPtsRef.size(), // flag = number of TofPoints generating the cluster
1601 0); //channel
1602 // vDigiIndRef);
1603 nHits++;
1604 if (event) event->AddData(ECbmDataType::kTofHit, hitIndex);
1605
1606 CbmMatch* digiMatch = new CbmMatch();
1607 for (UInt_t i = 0; i < vDigiIndRef.size(); i++) {
1608 digiMatch->AddLink(CbmLink(0., vDigiIndRef.at(i), iEventNr, iInputNr));
1609 }
1610 new ((*fTofDigiMatchColl)[hitIndex]) CbmMatch(*digiMatch);
1611 delete digiMatch;
1612
1613 // Using the SetLinks/GetLinks of the TofHit class seems to conflict
1614 // with something in littrack QA
1615 // ((CbmTofHit*)fTofHitsColl->At(fiNbHits))->SetLinks(vPtsRef[0]);
1616 //fiNbHits++; is in local variable now (VF)
1617 // For Histogramming
1618 fviClusterSize[iSmType][iRpc].push_back(iNbChanInHit);
1619 fviTrkMul[iSmType][iRpc].push_back(vPtsRef.size());
1620 fvdX[iSmType][iRpc].push_back(dWeightedPosX);
1621 fvdY[iSmType][iRpc].push_back(dWeightedPosY);
1622 /* no TofPoint available for data!
1623 fvdDifX[iSmType][iRpc].push_back( vPtsRef[0]->GetX() - dWeightedPosX);
1624 fvdDifY[iSmType][iRpc].push_back( vPtsRef[0]->GetY() - dWeightedPosY);
1625 fvdDifCh[iSmType][iRpc].push_back( fGeoHandler->GetCell( vPtsRef[0]->GetDetectorID() ) -1 -iLastChan );
1626 */
1627 vPtsRef.clear();
1628 vDigiIndRef.clear();
1629
1630 // Start a new hit
1631 dWeightedTime = dTime * dTotS;
1632 dWeightedPosX = dPosX * dTotS;
1633 dWeightedPosY = dPosY * dTotS;
1634 dWeightedPosZ = dPosZ * dTotS;
1635 dWeightedTimeErr = dTotS * dTotS;
1636 dWeightsSum = dTotS;
1637 iNbChanInHit = 1;
1638 // Save pointer on CbmTofPoint
1639 // vPtsRef.push_back( (CbmTofPoint*)(xDigiA->GetLinks()) );
1640 // Save next digi address
1641 LOG(debug4) << " Next fStor Digi " << iSmType << ", SR " << iSm * iNbRpc + iRpc << ", Ch" << iCh
1642 << ", Dig0 " << (fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][0]) << ", Dig1 "
1643 << (fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][1]);
1644
1645 vDigiIndRef.push_back((Int_t)(fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][0]));
1646 vDigiIndRef.push_back((Int_t)(fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][1]));
1647 LOG(debug3) << " Erase fStor entries(b) " << iSmType << ", SR " << iSm * iNbRpc + iRpc << ", Ch"
1648 << iCh;
1649 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1650 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
1651 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1652 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
1653 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1654 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
1655 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1656 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
1657
1658 if (kTRUE == fDigiBdfPar->ClustUseTrackId()) {
1659 // if( ((CbmTofPoint*)(xDigiA->GetLinks()))->GetTrackID() !=
1660 // ((CbmTofPoint*)(xDigiB->GetLinks()))->GetTrackID() )
1661 // if other side come from a different Track,
1662 // also save the pointer on CbmTofPoint
1663 // vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) );
1664 } // if( kTRUE == fDigiBdfPar->ClustUseTrackId() )
1665 //else if( xDigiA->GetLinks() != xDigiB->GetLinks() )
1666 // if other side come from a different TofPoint,
1667 // also save the pointer on CbmTofPoint
1668 // vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) );
1669 } // else of if current Digis compatible with last fired chan
1670 } // if( 0 < iNbChanInHit)
1671 else {
1672 LOG(debug1) << "CbmTofSimpClusterizer::BuildClusters: 1.Hit on "
1673 "channel "
1674 << iCh << ", time: " << dTime << ", x: " << dPosX << ", y: " << dPosY
1675 << ", Tot: " << dTotS;
1676
1677 // first fired strip in this RPC
1678 dWeightedTime = dTime * dTotS;
1679 dWeightedPosX = dPosX * dTotS;
1680 dWeightedPosY = dPosY * dTotS;
1681 dWeightedPosZ = dPosZ * dTotS;
1682 dWeightedTimeErr = dTotS * dTotS;
1683 dWeightsSum = dTotS;
1684 iNbChanInHit = 1;
1685 // Save pointer on CbmTofPoint
1686 //if(NULL != fTofPointsColl)
1687 // vPtsRef.push_back( (CbmTofPoint*)(xDigiA->GetLinks()) );
1688 vDigiIndRef.push_back((Int_t)(fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][0]));
1689 vDigiIndRef.push_back((Int_t)(fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][1]));
1690
1691 LOG(debug2) << " Erase fStor entries(c) " << iSmType << ", SR " << iSm * iNbRpc + iRpc << ", Ch"
1692 << iCh;
1693 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1694 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
1695 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1696 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
1697 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1698 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
1699 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1700 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
1701
1702 if (kTRUE == fDigiBdfPar->ClustUseTrackId()) {
1703 // if( ((CbmTofPoint*)(xDigiA->GetLinks()))->GetTrackID() !=
1704 // ((CbmTofPoint*)(xDigiB->GetLinks()))->GetTrackID() )
1705 // if other side come from a different Track,
1706 // also save the pointer on CbmTofPoint
1707 // vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) );
1708 } // if( kTRUE == fDigiBdfPar->ClustUseTrackId() )
1709 // else if( xDigiA->GetLinks() != xDigiB->GetLinks() )
1710 // if other side come from a different TofPoint,
1711 // also save the pointer on CbmTofPoint
1712 // vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) );
1713 } // else of if( 0 < iNbChanInHit)
1714 iLastChan = iCh;
1715 // dLastPosX = dPosX; // -> Comment to remove warning because set but never used
1716 dLastPosY = dPosY;
1717 dLastTime = dTime;
1718 } // if( 2 == fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size() )
1719 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].clear();
1720 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].clear();
1721 } // for( Int_t iCh = 0; iCh < iNbCh; iCh++ )
1722 LOG(debug2) << "CbmTofSimpClusterizer::BuildClusters: finished V-RPC"
1723 << Form(" %3d %3d %3d %d", iSmType, iSm, iRpc, fTofHitsColl->GetEntriesFast());
1724 } // else of if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
1725 } // if( 0 == iChType)
1726 else {
1727 LOG(fatal) << "CbmTofSimpClusterizer::BuildClusters => Cluster building "
1728 << "from digis to hits not implemented for pads, Sm type " << iSmType << " Rpc " << iRpc;
1729 return std::make_pair(0, 0);
1730 } // else of if( 0 == iChType)
1731
1732 // Now check if another hit/cluster is started
1733 // and save it if it's the case
1734 if (0 < iNbChanInHit) {
1735 LOG(debug1) << "CbmTofSimpClusterizer::BuildClusters: Process cluster " << iNbChanInHit;
1736
1737 // Check orientation to properly assign errors
1738 if (1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc)) {
1739 LOG(debug1) << "CbmTofSimpClusterizer::BuildClusters: H-Hit ";
1740 } // if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
1741 else {
1742 LOG(debug2) << "CbmTofSimpClusterizer::BuildClusters: V-Hit ";
1743 // Save Hit
1744 dWeightedTime /= dWeightsSum;
1745 dWeightedPosX /= dWeightsSum;
1746 dWeightedPosY /= dWeightsSum;
1747 dWeightedPosZ /= dWeightsSum;
1748 // simple error scaling with TOT
1749 // dWeightedTimeErr = TMath::Sqrt( dWeightedTimeErr ) * fdParSystTimeRes / dWeightsSum;
1750 // OR harcoded value
1751 dWeightedTimeErr = fdParSystTimeRes;
1752
1753 Double_t hitpos_local[3];
1754 Double_t hitpos[3];
1755
1756 TGeoNode* tNode = // prepare local->global trafo
1757 gGeoManager->FindNode(fTrafoCell->GetX(), fTrafoCell->GetY(), fTrafoCell->GetZ());
1758 //TGeoNode *cNode =
1759 gGeoManager->GetCurrentNode();
1760
1761 TGeoRotation rotMatrix(*gGeoManager->GetCurrentMatrix());
1762
1763 hitpos[0] = fTrafoCell->GetX();
1764 hitpos[1] = fTrafoCell->GetY();
1765 hitpos[2] = fTrafoCell->GetZ();
1766 gGeoManager->MasterToLocal(hitpos, hitpos_local);
1767 LOG(debug1) << Form(" Node0 at (%6.1f,%6.1f,%6.1f) : (%6.1f,%6.1f,%6.1f)", fChannelInfo->GetX(),
1768 fChannelInfo->GetY(), fChannelInfo->GetZ(), hitpos_local[0], hitpos_local[1],
1769 hitpos_local[2]);
1770
1771 hitpos_local[0] += dWeightedPosX;
1772 hitpos_local[1] += dWeightedPosY;
1773 hitpos_local[2] += dWeightedPosZ;
1774
1775 gGeoManager->LocalToMaster(hitpos_local, hitpos);
1776 LOG(debug1) << Form(" LTM for V-node %p, info %p, tra %p: (%6.1f,%6.1f,%6.1f) "
1777 "->(%6.1f,%6.1f,%6.1f) [(%6.1f,%6.1f,%6.1f)]",
1778 tNode, fChannelInfo, fTrafoCell, hitpos_local[0], hitpos_local[1], hitpos_local[2],
1779 hitpos[0], hitpos[1], hitpos[2], fTrafoCell->GetX(), fTrafoCell->GetY(),
1780 fTrafoCell->GetZ());
1781
1782 TVector3 hitPos(hitpos[0], hitpos[1], hitpos[2]);
1783 // TestBeam errors, not properly done at all for now
1784 // Right way of doing it should take into account the weight distribution
1785 // and real system time resolution
1786 Double_t hiterr_local[3];
1787 Double_t hiterr[3];
1788
1789 hiterr_local[0] = fTrafoCell->GetSizex() / TMath::Sqrt(12.0); // Single strips approximation
1790 hiterr_local[1] = fdParFeeTimeRes * fvCPSigPropSpeed[iSmType][iRpc]; // Use the electronics resolution
1791 hiterr_local[2] = fDigiBdfPar->GetNbGaps(iSmType, iRpc) * fDigiBdfPar->GetGapSize(iSmType, iRpc)
1792 / 10.0 // Change gap size in cm
1793 / TMath::Sqrt(12.0); // Use full RPC thickness as "Channel" Z size
1794
1795 rotMatrix.LocalToMaster(hiterr_local, hiterr);
1796 // TVector3 hitPosErr( hiterr_local[0], hiterr_local[1], hiterr_local[2] );
1797 TVector3 hitPosErr(fabs(hiterr[0]), fabs(hiterr[1]), fabs(hiterr[2]));
1798 /*
1799 LOG(info)<< " Size X " << fTrafoCell->GetSizex()
1800 ;
1801 LOG(info)<< " Nb gaps " << fDigiBdfPar->GetNbGaps( iSmType, iRpc)
1802 << " gap size " << fDigiBdfPar->GetGapSize( iSmType, iRpc)
1803 ;
1804 LOG(info)<<"CbmTofSimpClusterizer::BuildClusters: Errors in local "
1805 << hiterr_local[0] << " "<< hiterr_local[1] << " " << hiterr_local[2]
1806 ;
1807 LOG(info)<<"CbmTofSimpClusterizer::BuildClusters: Errors in master "
1808 << hiterr[0] << " "<< hiterr[1] << " " << hiterr[2]
1809 ;
1810*/
1811 // cout<<"a "<<vPtsRef.size()<<endl;
1812 // cout<<"b "<<vPtsRef[0]<<endl;
1813 // cout<<"c "<<vPtsRef[0]->GetDetectorID()<<endl;
1814 // Int_t iDetId = vPtsRef[0]->GetDetectorID();// detID = pt->GetDetectorID() <= from TofPoint
1815 // Int_t iDetId = iChId;
1816
1817 Int_t iChm = iTrafoCell + floor(dWeightedPosX / fChannelInfo->GetSizex());
1818 if (iChm < 0 || iChm > iNbCh) {
1819 LOG(debug) << "CbmTofSimpClusterizer::BuildClusters: Invalid "
1820 "mean channel "
1821 << iChm << "(" << iNbCh << ")";
1822 }
1823 Int_t iDetId = CbmTofAddress::GetUniqueAddress(iSm, iRpc, iChm, 0, iSmType);
1824 Int_t iRefId = 0; // Index of the correspondng TofPoint
1825 // if(NULL != fTofPointsColl) iRefId = fTofPointsColl->IndexOf( vPtsRef[0] );
1826 LOG(debug) << "CbmTofSimpClusterizer::BuildClusters: Save V-Hit "
1827 // << Form(" %3d %3d 0x%08x %3d %3d %3d 0x%08x",
1828 << Form(" %3d %3d %3d %3d %3d ", fiNbHits, iNbChanInHit, iDetId, iLastChan,
1829 iRefId) //vPtsRef.size(),vPtsRef[0])
1830 // dWeightedTime,dWeightedPosY)
1831 << ", DigiSize: " << vDigiIndRef.size();
1832 LOG(debug) << ", DigiInds: ";
1833 for (UInt_t i = 0; i < vDigiIndRef.size(); i++) {
1834 LOG(debug) << " " << vDigiIndRef.at(i);
1835 }
1836 LOG(debug);
1837
1838 fviClusterMul[iSmType][iSm][iRpc]++;
1839
1840 hitIndex = fTofHitsColl->GetEntriesFast();
1841 new ((*fTofHitsColl)[hitIndex])
1842 CbmTofHit(iDetId, hitPos, hitPosErr,
1843 fiNbHits, //iRefId
1844 dWeightedTime * fOutTimeFactor, dWeightedTimeErr, vPtsRef.size(), 0);
1845 // vDigiIndRef);
1846 nHits++;
1847 if (event) event->AddData(ECbmDataType::kTofHit, hitIndex);
1848 CbmMatch* digiMatch = new CbmMatch();
1849 for (UInt_t i = 0; i < vDigiIndRef.size(); i++) {
1850 digiMatch->AddLink(CbmLink(0., vDigiIndRef.at(i), iEventNr, iInputNr));
1851 }
1852 new ((*fTofDigiMatchColl)[hitIndex]) CbmMatch(*digiMatch);
1853 delete digiMatch;
1854
1855 // Using the SetLinks/GetLinks of the TofHit class seems to conflict
1856 // with something in littrack QA
1857 // ((CbmTofHit*)fTofHitsColl->At(fiNbHits))->SetLinks(vPtsRef[0]);
1858 // fiNbHits++; is in local variable (VF)
1859 // For Histogramming
1860 fviClusterSize[iSmType][iRpc].push_back(iNbChanInHit);
1861 fviTrkMul[iSmType][iRpc].push_back(vPtsRef.size());
1862 fvdX[iSmType][iRpc].push_back(dWeightedPosX);
1863 fvdY[iSmType][iRpc].push_back(dWeightedPosY);
1864 /*
1865 fvdDifX[iSmType][iRpc].push_back( vPtsRef[0]->GetX() - dWeightedPosX);
1866 fvdDifY[iSmType][iRpc].push_back( vPtsRef[0]->GetY() - dWeightedPosY);
1867 fvdDifCh[iSmType][iRpc].push_back( fGeoHandler->GetCell( vPtsRef[0]->GetDetectorID() ) -1 -iLastChan );
1868 */
1869 vPtsRef.clear();
1870 vDigiIndRef.clear();
1871 } // else of if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
1872 } // if( 0 < iNbChanInHit)
1873 LOG(debug4) << " Fini-A " << Form(" %3d %3d %3d ", iSmType, iSm, iRpc);
1874 } // for each sm/rpc pair
1875 LOG(debug3) << " Fini-B " << Form(" %3d ", iSmType);
1876 } // for( Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++ )
1877 } // if( kTRUE == fDigiBdfPar->UseExpandedDigi() )
1878 else {
1879 LOG(error) << " Compressed Digis not implemented ... ";
1880 }
1881 return std::make_pair(nDigis, nHits);
1882}
1883
1884void CbmTofSimpClusterizer::GetEventInfo(Int_t& inputNr, Int_t& eventNr, Double_t& eventTime)
1885{
1886
1887 // --- In a FairRunAna, take the information from FairEventHeader
1888 if (FairRunAna::Instance()) {
1889 FairEventHeader* event = FairRunAna::Instance()->GetEventHeader();
1890 inputNr = event->GetInputFileId();
1891 eventNr = event->GetMCEntryNumber();
1892 eventTime = event->GetEventTime();
1893 }
1894
1895 // --- In a FairRunSim, the input number and event time are always zero;
1896 // --- only the event number is retrieved.
1897 else {
1898 if (!FairRunSim::Instance()) LOG(fatal) << GetName() << ": neither SIM nor ANA run.";
1899 FairMCEventHeader* event = FairRunSim::Instance()->GetMCEventHeader();
1900 inputNr = 0;
1901 eventNr = event->GetEventID();
1902 eventTime = 0.;
1903 }
1904}
@ kTof
Time-of-flight Detector.
static FairRootManager * rootMgr
static TFile * fHist
TClonesArray * fTofHitsColl
CbmDigiManager * fDigiMan
const Int_t iNTrg
const Int_t nbClWalkBinX
const Int_t nbClDelTofBinX
@ k12b
@ k21a
@ k14a
const Int_t iNTrg
const Int_t nbClWalkBinX
Double_t TOTMax
const Int_t nbClDelTofBinX
Double_t TOTMin
const Double_t TTotMean
std::vector< Int_t > vDigiIndRef
std::vector< CbmTofPoint * > vPtsRef
static constexpr size_t size()
Definition KfSimdPseudo.h:2
static Int_t GetNofDigis(ECbmModuleId systemId)
static Bool_t IsPresent(ECbmModuleId systemId)
Presence of a digi branch.
InitStatus Init()
Initialisation.
const Digi * Get(Int_t index) const
Get a digi object.
static CbmDigiManager * Instance()
Static instance.
Class characterising one event by a collection of links (indices) to data objects,...
Definition CbmEvent.h:34
int32_t GetMotherId() const
Definition CbmMCTrack.h:69
int32_t GetNPoints(ECbmModuleId detId) const
void AddLink(const CbmLink &newLink)
Definition CbmMatch.cxx:47
static uint32_t GetUniqueAddress(uint32_t Sm, uint32_t Rpc, uint32_t Channel, uint32_t Side=0, uint32_t SmType=0, uint32_t RpcType=0)
Double_t GetSizey() const
Definition CbmTofCell.h:40
Double_t GetY() const
Definition CbmTofCell.h:36
Double_t GetSizex() const
Definition CbmTofCell.h:39
Double_t GetX() const
Definition CbmTofCell.h:35
Double_t GetZ() const
Definition CbmTofCell.h:37
virtual InitStatus Init()
virtual int32_t SetDetectorInfo(const CbmTofDetectorInfo detectorInfo)=0
Parameters class for the CBM ToF digitizer using beam data distributions.
Double_t GetFeeTimeRes() const
Int_t GetNbSmTypes() const
Double_t GetGapSize(Int_t iSmType, Int_t iRpc) const
Int_t GetNbSm(Int_t iSmType) const
Double_t GetMaxDistAlongCh() const
Int_t GetNbChan(Int_t iSmType, Int_t iRpc) const
Bool_t ClustUseTrackId() const
Int_t GetNbGaps(Int_t iSmType, Int_t iRpc) const
Int_t GetChanType(Int_t iSmType, Int_t iRpc) const
Int_t GetNbRpc(Int_t iSmType) const
Bool_t UseExpandedDigi() const
Double_t GetMaxTimeDist() const
Int_t GetChanOrient(Int_t iSmType, Int_t iRpc) const
Double_t GetSignalSpeed() const
Double_t GetSigVel(Int_t iSmType, Int_t iSm, Int_t iRpc) const
CbmTofCell * GetCell(Int_t i)
Int_t GetNrOfModules()
Int_t GetCellId(Int_t i)
Data class for expanded digital TOF information.
Definition CbmTofDigi.h:47
double GetSide() const
Channel Side.
Definition CbmTofDigi.h:160
void SetTot(double tot)
Definition CbmTofDigi.h:166
double GetChannel() const
Channel .
Definition CbmTofDigi.h:156
double GetSm() const
Sm.
Definition CbmTofDigi.h:144
double GetTime() const
Inherited from CbmDigi.
Definition CbmTofDigi.h:131
double GetType() const
Sm Type .
Definition CbmTofDigi.h:148
double GetRpc() const
Detector aka Module aka RPC .
Definition CbmTofDigi.h:152
double GetTot() const
Alias for GetCharge.
Definition CbmTofDigi.h:140
void SetTime(double time)
Definition CbmTofDigi.h:165
Int_t GetCell(Int_t uniqueId)
Int_t Init(Bool_t isSimulation=kFALSE)
Int_t GetSModule(Int_t uniqueId)
Int_t GetCounter(Int_t uniqueId)
Int_t GetSMType(Int_t uniqueId)
int32_t GetFlag() const
Definition CbmTofHit.h:75
void GetEventInfo(Int_t &inputNr, Int_t &eventNr, Double_t &eventTime)
Retrieve event info from run manager to properly fill the CbmLink objects.
Bool_t DeleteGeometry()
Delete the geometry related arrays: for now just clearing the Digis temporary vectors.
Bool_t InitCalibParameter()
Initialize other parameters not included in parameter classes.
std::vector< std::vector< std::vector< Int_t > > > fviClusterSize
Bool_t RegisterOutputs()
Create and register output TClonesArray of Tof Hits.
std::vector< std::vector< std::vector< std::vector< Int_t > > > > fStorDigiInd
std::vector< std::vector< std::vector< std::vector< std::vector< Double_t > > > > > fvCPWalk
CbmTofDigiBdfPar * fDigiBdfPar
Double_t fdNofHitsTot
Total number of hits produced.
std::vector< std::vector< Double_t > > fvCPSigPropSpeed
std::vector< std::vector< std::vector< Double_t > > > fvdX
std::vector< std::vector< std::vector< Double_t > > > fvdDifCh
Double_t fNofDigisUsed
Total number of Tof Digis processed.
Int_t fiNofEvents
Total number of events processed.
std::vector< std::vector< std::vector< std::vector< Double_t > > > > fvCPDelTof
std::vector< std::vector< std::vector< Double_t > > > fvdY
Bool_t SetHistoFileName(TString sFilenameIn="./tofSimpClust.hst.root")
std::pair< Int_t, Int_t > BuildClusters(CbmEvent *event)
Build clusters out of ToF Digis and store the resulting info in a TofHit.
std::vector< std::vector< std::vector< std::vector< Double_t > > > > fvCPTotGain
virtual void SetParContainers()
Inherited from FairTask.
Bool_t RegisterInputs()
Recover pointer on input TClonesArray: TofPoints, TofDigis...
virtual ~CbmTofSimpClusterizer()
Destructor.
std::vector< std::vector< std::vector< std::vector< CbmTofDigi * > > > > fStorDigiExp
Bool_t LoadGeometry()
Load the geometry: for now just resizing the Digis temporary vectors.
Double_t fdTimeTot
Total execution time.
Bool_t InitParameters()
Initialize other parameters not included in parameter classes.
std::vector< std::vector< std::vector< Double_t > > > fvdDifX
CbmTofGeoHandler * fGeoHandler
std::vector< std::vector< std::vector< std::vector< Double_t > > > > fvCPTOff
Int_t fiNofTs
Number of processed timeslices.
virtual void Exec(Option_t *option)
Inherited from FairTask.
virtual InitStatus Init()
Inherited from FairTask.
TStopwatch fTimer
ROOT timer.
std::vector< std::vector< std::vector< Int_t > > > fviTrkMul
virtual void Finish()
Inherited from FairTask.
Double_t fNofDigisAll
Total number of TOF digis in input.
std::vector< std::vector< std::vector< Double_t > > > fvdDifY
std::vector< std::vector< std::vector< Int_t > > > fviClusterMul