CbmRoot
Loading...
Searching...
No Matches
CbmTofHitMaker.cxx
Go to the documentation of this file.
1/* Copyright (C) 2018-2021 PI-UHd/GSI, Heidelberg/Darmstadt
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Norbert Herrmann [committer] */
4
16#include "CbmTofHitMaker.h"
17
18// TOF Classes and includes
19#include "CbmBmonDigi.h" // in cbmdata/bmon
20#include "CbmDigiManager.h"
21#include "CbmEvent.h"
22#include "CbmMatch.h"
23#include "CbmTofAddress.h" // in cbmdata/tof
24#include "CbmTofCell.h" // in tof/TofData
25#include "CbmTofCreateDigiPar.h" // in tof/TofTools
26#include "CbmTofDetectorId_v12b.h" // in cbmdata/tof
27#include "CbmTofDetectorId_v14a.h" // in cbmdata/tof
28#include "CbmTofDetectorId_v21a.h" // in cbmdata/tof
29#include "CbmTofDigi.h" // in cbmdata/tof
30#include "CbmTofDigiBdfPar.h" // in tof/TofParam
31#include "CbmTofDigiPar.h" // in tof/TofParam
32#include "CbmTofGeoHandler.h" // in tof/TofTools
33#include "CbmTofHit.h" // in cbmdata/tof
34#include "CbmTofPoint.h" // in cbmdata/tof
35#include "CbmVertex.h"
36#include "TTrbHeader.h"
37
38// CBMroot classes and includes
39#include "CbmMCTrack.h"
40
41// FAIR classes and includes
42#include "FairEventHeader.h"
43#include "FairRootFileSink.h"
44#include "FairRootManager.h"
45#include "FairRunAna.h"
46#include "FairRuntimeDb.h"
47#include "Logger.h"
48
49// ROOT Classes and includes
50#include "TClonesArray.h"
51#include "TDirectory.h"
52#include "TF1.h"
53#include "TF2.h"
54#include "TFitResult.h"
55#include "TGeoManager.h"
56#include "TH1.h"
57#include "TH2.h"
58#include "TH3.h"
59#include "TLine.h"
60#include "TMath.h"
61#include "TMinuit.h"
62#include "TProfile.h"
63#include "TROOT.h"
64#include "TRandom.h"
65#include "TRandom3.h"
66#include "TVector3.h"
67
68// Constants definitions
70
71// C++ Classes and includes
72// Globals
73#include <vector>
74
75static Int_t iNbTs = 0;
76
77static Bool_t bAddBeamCounterSideDigi = kTRUE;
78static TRandom3* fRndm = new TRandom3();
79
80// std::vector< CbmTofPoint* > vPtsRef;
81
83
84/************************************************************************************/
86{
87 // if ( !fInstance ) fInstance = this;
88}
89
90CbmTofHitMaker::CbmTofHitMaker(const char* name, Int_t verbose, Bool_t writeDataInOut)
91 : FairTask(TString(name), verbose)
92 , fGeoHandler(new CbmTofGeoHandler())
93 , fTofId(NULL)
94 , fDigiPar(NULL)
95 , fChannelInfo(NULL)
96 , fDigiBdfPar(NULL)
97 , fTrbHeader(NULL)
98 , fTofPointsColl(NULL)
99 , fMcTracksColl(NULL)
100 , fDigiMan(nullptr)
101 , fEventsColl(nullptr)
102 , fbWriteHitsInOut(writeDataInOut)
103 , fbWriteDigisInOut(writeDataInOut)
104 , fTofHitsColl(NULL)
105 , fTofDigiMatchColl(NULL)
106 , fTofHitsCollOut(NULL)
107 , fTofDigiMatchCollOut(NULL)
108 , fiNbHits(0)
109 , fVerbose(verbose)
110 , fStorDigi()
111 , fStorDigiInd()
112 , vDigiIndRef()
113 , fviClusterMul()
114 , fviClusterSize()
115 , fviTrkMul()
116 , fvdX()
117 , fvdY()
118 , fvdDifX()
119 , fvdDifY()
120 , fvdDifCh()
121 , fhClustBuildTime(NULL)
122 , fhHitsPerTracks(NULL)
123 , fhPtsPerHit(NULL)
124 , fhTimeResSingHits(NULL)
125 , fhTimeResSingHitsB(NULL)
126 , fhTimePtVsHits(NULL)
127 , fhClusterSize(NULL)
128 , fhClusterSizeType(NULL)
129 , fhTrackMul(NULL)
130 , fhClusterSizeMulti(NULL)
131 , fhTrk1MulPos(NULL)
132 , fhHiTrkMulPos(NULL)
133 , fhAllTrkMulPos(NULL)
134 , fhMultiTrkProbPos(NULL)
135 , fhDigSpacDifClust(NULL)
136 , fhDigTimeDifClust(NULL)
137 , fhDigDistClust(NULL)
138 , fhClustSizeDifX(NULL)
139 , fhClustSizeDifY(NULL)
140 , fhChDifDifX(NULL)
141 , fhChDifDifY(NULL)
142 , fhCluMulCorDutSel(NULL)
143 , fhEvCluMul(NULL)
144 , fhRpcDigiCor()
145 , fhRpcDigiMul()
146 , fhRpcDigiStatus()
147 , fhRpcDigiDTLD()
148 , fhRpcDigiDTFD()
149 , fhRpcDigiDTMul()
150 , fhRpcCluMul()
151 , fhRpcCluRate()
152 , fhRpcCluRate10s()
153 , fhRpcCluPosition()
154 , fhRpcCluPositionEvol()
155 , fhRpcCluTimeEvol()
156 , fhRpcCluDelPos()
157 , fhRpcCluDelMatPos()
158 , fhRpcCluTOff()
159 , fhRpcCluDelTOff()
160 , fhRpcCluDelMatTOff()
161 , fhRpcCluTrms()
162 , fhRpcCluTot()
163 , fhRpcCluSize()
164 , fhRpcCluAvWalk()
165 , fhRpcCluAvLnWalk()
166 , fhRpcCluWalk()
167 , fhSmCluPosition()
168 , fhSmCluTOff()
169 , fhSmCluSvel()
170 , fhSmCluFpar()
171 , fhRpcDTLastHits()
172 , fhRpcDTLastHits_Tot()
173 , fhRpcDTLastHits_CluSize()
174 , fhTRpcCluMul()
175 , fhTRpcCluPosition()
176 , fhTRpcCluTOff()
177 , fhTRpcCluTofOff()
178 , fhTRpcCluTot()
179 , fhTRpcCluSize()
180 , fhTRpcCluAvWalk()
181 , fhTRpcCluDelTof()
182 , fhTRpcCludXdY()
183 , fhTRpcCluWalk()
184 , fhTRpcCluWalk2()
185 , fhTSmCluPosition()
186 , fhTSmCluTOff()
187 , fhTSmCluTRun()
188 , fhTRpcCluTOffDTLastHits()
189 , fhTRpcCluTotDTLastHits()
190 , fhTRpcCluSizeDTLastHits()
191 , fhTRpcCluMemMulDTLastHits()
192 , fhSeldT()
193 , fvCPDelTof()
194 , fvCPTOff()
195 , fvCPTotGain()
196 , fvCPTotOff()
197 , fvCPWalk()
198 , fvLastHits()
199 , fvDeadStrips()
200 , fvTimeLastDigi()
201 , fiNbSameSide(0)
202 , fhNbSameSide(NULL)
203 , fhNbDigiPerChan(NULL)
204 , fStart()
205 , fStop()
206 , dTRef(0.)
207 , fdTRefMax(0.)
208 , fCalMode(0)
209 , fCalSel(0)
210 , fCalSmAddr(0)
211 , fdCaldXdYMax(0.)
212 , fiCluMulMax(0)
213 , fTRefMode(0)
214 , fTRefHits(0)
215 , fIdMode(0)
216 , fDutId(0)
217 , fDutSm(0)
218 , fDutRpc(0)
219 , fDutAddr(0)
220 , fSelId(0)
221 , fSelSm(0)
222 , fSelRpc(0)
223 , fSelAddr(0)
224 , fiBeamRefType(0)
225 , fiBeamRefSm(0)
226 , fiBeamRefDet(0)
227 , fiBeamRefAddr(0)
228 , fiBeamRefMulMax(1)
229 , fiBeamAddRefMul(0)
230 , fSel2Id(-1)
231 , fSel2Sm(0)
232 , fSel2Rpc(0)
233 , fSel2Addr(0)
234 , fSel2MulMax(1)
235 , fDetIdIndexMap()
236 , fviDetId()
237 , fPosYMaxScal(0.)
238 , fTRefDifMax(0.)
239 , fTotMax(0.)
240 , fTotMin(0.)
241 , fTotOff(0.)
242 , fTotMean(0.)
243 , fdDelTofMax(60.)
244 , fTotPreRange(0.)
245 , fMaxTimeDist(0.)
246 , fdChannelDeadtime(0.)
247 , fdMemoryTime(0.)
248 , fdYFitMin(1.E6)
249 , fdToDAv(0.033) // in ns/cm
250 , fEnableMatchPosScaling(kTRUE)
251 , fEnableAvWalk(kFALSE)
252 , fbPs2Ns(kFALSE)
253 , fCalParFileName("")
254 , fOutHstFileName("")
255 , fCalParFile(NULL)
256 , fiNevtBuild(0)
257 , fiMsgCnt(100)
258 , fdTOTMax(50.)
259 , fdTOTMin(0.)
260 , fdTTotMean(2.)
261 , fdMaxTimeDist(0.)
262 , fdMaxSpaceDist(0.)
263 , fdEvent(0)
264 , fbSwapChannelSides(kFALSE)
265 , fiOutputTreeEntry(0)
266 , fiFileIndex(0)
267 , fbAlternativeBranchNames(kFALSE)
268{
269 if (!fInstance) fInstance = this;
270}
271
273{
274 if (fGeoHandler) delete fGeoHandler;
275 if (fInstance == this) fInstance = 0;
276 // DeleteHistos(); // <-- if needed ?
277}
278
279/************************************************************************************/
280// FairTasks inherited functions
282{
283 LOG(info) << "CbmTofHitMaker initializing... expect Digis in ns units! ";
284
285 if (kFALSE == RegisterInputs()) return kFATAL;
286
287 if (kFALSE == RegisterOutputs()) return kFATAL;
288
289 if (kFALSE == InitParameters()) return kFATAL;
290
291 if (kFALSE == LoadGeometry()) return kFATAL;
292
293 if (kFALSE == InitCalibParameter()) return kFATAL;
294
295 if (kFALSE == CreateHistos()) return kFATAL;
296 return kSUCCESS;
297}
298
300{
301 LOG(info) << "=> Get the digi parameters for tof";
302 //LOG(warning)<<"Return without action";
303 //return;
304 // Get Base Container
305 FairRunAna* ana = FairRunAna::Instance();
306 FairRuntimeDb* rtdb = ana->GetRuntimeDb();
307
308 fDigiPar = (CbmTofDigiPar*) (rtdb->getContainer("CbmTofDigiPar"));
309
310 LOG(info) << "found " << fDigiPar->GetNrOfModules() << " cells ";
311 fDigiBdfPar = (CbmTofDigiBdfPar*) (rtdb->getContainer("CbmTofDigiBdfPar"));
312}
313
314void CbmTofHitMaker::Exec(Option_t* option)
315{
316
318 if (fEventsColl) {
319 LOG(info) << "CbmTofHitMaker::Exec => New timeslice " << iNbTs << " with " << fEventsColl->GetEntriesFast()
320 << " events, " << fDigiMan->GetNofDigis(ECbmModuleId::kTof) << " TOF digis + "
321 << fDigiMan->GetNofDigis(ECbmModuleId::kBmon) << " Bmon digis ";
322 iNbTs++;
323
324 Int_t iNbHits = 0;
325 Int_t iNbCalDigis = 0;
326 fTofDigiMatchCollOut->Delete(); // costly, FIXME
327 fTofHitsCollOut->Delete(); // costly, FIXME
328 //fTofDigiMatchCollOut->Clear("C"); // not sufficient, memory leak
329 for (Int_t iEvent = 0; iEvent < fEventsColl->GetEntriesFast(); iEvent++) {
330 CbmEvent* tEvent = dynamic_cast<CbmEvent*>(fEventsColl->At(iEvent));
331 fTofDigiVec.clear();
332 //if (fTofDigisColl) fTofDigisColl->Clear("C");
333 //Int_t iNbDigis=0; (VF) not used
334 LOG(debug) << "TS event " << iEvent << " with " << tEvent->GetNofData(ECbmDataType::kBmonDigi) << " Bmon and "
335 << tEvent->GetNofData(ECbmDataType::kTofDigi) << " Tof digis ";
336
337 for (size_t iDigi = 0; iDigi < tEvent->GetNofData(ECbmDataType::kBmonDigi); iDigi++) {
338 Int_t iDigiIndex = static_cast<Int_t>(tEvent->GetIndex(ECbmDataType::kBmonDigi, iDigi));
339 CbmTofDigi tDigi(fDigiMan->Get<CbmBmonDigi>(iDigiIndex));
340 if (tDigi.GetType() != 5)
341 LOG(fatal) << "Wrong T0 type " << tDigi.GetType() << ", Addr 0x" << std::hex << tDigi.GetAddress();
342 fTofDigiVec.push_back(tDigi);
343 }
344 for (size_t iDigi = 0; iDigi < tEvent->GetNofData(ECbmDataType::kTofDigi); iDigi++) {
345 Int_t iDigiIndex = static_cast<Int_t>(tEvent->GetIndex(ECbmDataType::kTofDigi, iDigi));
346 const CbmTofDigi* tDigi = fDigiMan->Get<CbmTofDigi>(iDigiIndex);
347 fTofDigiVec.push_back(CbmTofDigi(*tDigi));
348 //new((*fTofDigisColl)[iNbDigis++]) CbmTofDigi(*tDigi);
349 }
350
351 ExecEvent(option);
352
353 // --- In event-by-event mode: copy caldigis, hits and matches to output array and register them to event
354
355 // Int_t iDigi0=iNbCalDigis; //starting index of current event (VF) not used
356 for (UInt_t index = 0; index < fTofCalDigiVec->size(); index++) {
357 // for (Int_t index = 0; index < fTofCalDigisColl->GetEntriesFast(); index++){
358 CbmTofDigi* tDigi = &(fTofCalDigiVec->at(index));
359 //CbmTofDigi* tDigi = dynamic_cast<CbmTofDigi*>(fTofCalDigisColl->At(index));
360 tEvent->AddData(ECbmDataType::kTofCalDigi, iNbCalDigis);
361 fTofCalDigiVecOut->push_back(CbmTofDigi(*tDigi));
362 iNbCalDigis++;
363 //new((*fTofCalDigisCollOut)[iNbCalDigis++]) CbmTofDigi(*tDigi);
364 }
365
366 for (Int_t index = 0; index < fTofHitsColl->GetEntriesFast(); index++) {
367 CbmTofHit* pHit = (CbmTofHit*) fTofHitsColl->At(index);
368 new ((*fTofHitsCollOut)[iNbHits]) CbmTofHit(*pHit);
369 tEvent->AddData(ECbmDataType::kTofHit, iNbHits);
370
371 CbmMatch* pDigiMatch = (CbmMatch*) fTofDigiMatchColl->At(index);
372 // update content of match object, not necessary if event definition is kept !
373 /*
374 for (Int_t iLink=0; iLink<pDigiMatch->GetNofLinks(); iLink++) { // loop over digis
375 CbmLink Link = pDigiMatch->GetLink(iLink);
376 Link.SetIndex(Link.GetIndex()+iDigi0);
377 }
378 */
379 new ((*fTofDigiMatchCollOut)[iNbHits]) CbmMatch(*pDigiMatch);
380
381 iNbHits++;
382 }
383 //fTofDigisColl->Delete();
384 fTofDigiVec.clear();
385 //fTofCalDigi->Delete();//Clear("C"); //otherwise memoryleak! FIXME
386 fTofCalDigiVec->clear();
387 fTofHitsColl->Clear("C");
388 fTofDigiMatchColl->Delete(); //Clear("C");
389 }
390 }
391 else {
392 // fTofDigisColl=fTofRawDigisColl;
393 // (VF) This does not work here. The digi manager does not foresee to add
394 // new data to the input array. So, I here copy the input digis into
395 // the array fTofDigisColl. Not very efficient, but temporary only, until
396 // also the internal data representations are changed to std::vectors.
397
398 fTofDigiVec.clear();
399 //if (fTofDigisColl) fTofDigisColl->Clear("C");
400 // Int_t iNbDigis=0; (VF) not used
401 for (Int_t iDigi = 0; iDigi < fDigiMan->GetNofDigis(ECbmModuleId::kTof); iDigi++) {
402 const CbmTofDigi* tDigi = fDigiMan->Get<CbmTofDigi>(iDigi);
403 fTofDigiVec.push_back(CbmTofDigi(*tDigi));
404 //new((*fTofDigisColl)[iNbDigis++]) CbmTofDigi(*tDigi);
405 }
406 ExecEvent(option);
407 }
408}
409
410void CbmTofHitMaker::ExecEvent(Option_t* /*option*/)
411{
412 // Clear output arrays
413 //fTofCalDigisColl->Delete(); //otherwise memoryleak if 'CbmDigi::fMatch' points to valid MC match objects (simulation)! FIXME
414 fTofCalDigiVec->clear();
415 fTofHitsColl->Clear("C");
416 //fTofHitsColl->Delete(); // Computationally costly!, but hopefully safe
417 //for (Int_t i=0; i<fTofDigiMatchColl->GetEntriesFast(); i++) ((CbmMatch *)(fTofDigiMatchColl->At(i)))->ClearLinks(); // FIXME, try to tamper memory leak (did not help)
418 //fTofDigiMatchColl->Clear("C+L"); // leads to memory leak
419 fTofDigiMatchColl->Delete();
420 FairRootFileSink* bla = (FairRootFileSink*) FairRootManager::Instance()->GetSink();
421 if (bla) fiOutputTreeEntry = ((FairRootFileSink*) FairRootManager::Instance()->GetSink())->GetOutTree()->GetEntries();
422
423 fiNbHits = 0;
424
425 fStart.Set();
426
428
430
431 fStop.Set();
432
433 fdEvent++;
434 FillHistos();
435
436 // fTofDigisColl->RemoveAll();
437}
438
439/************************************************************************************/
441{
442 if (fdEvent < 100) return; // don't save histos with insufficient statistics
443 WriteHistos();
444 // Prevent them from being sucked in by the CbmHadronAnalysis WriteHistograms method
445 // DeleteHistos();
446 if (fdMemoryTime > 0.) CleanLHMemory();
447}
448
449void CbmTofHitMaker::Finish(Double_t calMode)
450{
451 if (fdEvent < 100) return; // don't save histos with insufficient statistics
452 SetCalMode(calMode);
453 WriteHistos();
454}
455
456/************************************************************************************/
457// Functions common for all clusters approximations
459{
460 FairRootManager* fManager = FairRootManager::Instance();
461
462 if (NULL == fManager) {
463 LOG(error) << "CbmTofHitMaker::RegisterInputs => Could not find "
464 "FairRootManager!!!";
465 return kFALSE;
466 } // if( NULL == fTofDigisColl)
467
468 fEventsColl = dynamic_cast<TClonesArray*>(fManager->GetObject("Event"));
469 if (NULL == fEventsColl) fEventsColl = dynamic_cast<TClonesArray*>(fManager->GetObject("CbmEvent"));
470
471 if (NULL == fEventsColl) LOG(info) << "CbmEvent not found in input file, assume eventwise input";
472
474 fDigiMan->Init();
476 LOG(error) << GetName() << ": No Tof digi input!";
477 return kFALSE;
478 }
480 LOG(info) << GetName() << ": found separate Bmon digi input!";
481 }
482 else {
483 LOG(info) << "No separate Bmon digi input found.";
484 } // if( ! fT0DigiVec )
485
486 fTrbHeader = (TTrbHeader*) fManager->GetObject("TofTrbHeader.");
487 if (NULL == fTrbHeader) {
488 LOG(info) << "CbmTofHitMaker::RegisterInputs => Could not get "
489 "TofTrbHeader Object";
490 }
491
492 if (NULL == fEventsColl) {
493 //fTofDigisColl = new TClonesArray("CbmTofDigi");
494 }
495 else {
496 // time based input
497 LOG(info) << "CbmEvent found in input file, assume time based input";
498 //fTofDigisColl = new TClonesArray("CbmTofDigi");
499 }
500
501 return kTRUE;
502}
504{
505 FairRootManager* rootMgr = FairRootManager::Instance();
506 // FairRunAna* ana = FairRunAna::Instance(); (VF) not used
507
508 rootMgr->InitSink();
509
510 //fTofCalDigisColl = new TClonesArray("CbmTofDigi");
511 fTofCalDigiVec = new std::vector<CbmTofDigi>();
512
513 fTofHitsColl = new TClonesArray("CbmTofHit");
514
515 fTofDigiMatchColl = new TClonesArray("CbmMatch", 100);
516
517 TString tHitBranchName;
518 TString tHitDigiMatchBranchName;
519
521 tHitBranchName = "ATofHit";
522 tHitDigiMatchBranchName = "ATofDigiMatch";
523 }
524 else {
525 tHitBranchName = "TofHit";
526 tHitDigiMatchBranchName = "TofHitCalDigiMatch";
527 }
528
529 if (NULL == fEventsColl) {
530 // Flag check to control whether digis are written in output root file
531 //rootMgr->Register( "TofCalDigi","Tof", fTofCalDigisColl, fbWriteDigisInOut);
532 rootMgr->RegisterAny("TofCalDigi", fTofCalDigiVec, fbWriteDigisInOut);
533
534 // Flag check to control whether digis are written in output root file
535 rootMgr->Register(tHitBranchName, "Tof", fTofHitsColl, fbWriteHitsInOut);
536
537 rootMgr->Register(tHitDigiMatchBranchName, "Tof", fTofDigiMatchColl, fbWriteHitsInOut);
538 }
539 else { // CbmEvent - mode
540 //fTofCalDigisCollOut = new TClonesArray("CbmTofDigi");
541 fTofCalDigiVecOut = new std::vector<CbmTofDigi>();
542 fTofHitsCollOut = new TClonesArray("CbmTofHit");
543 fTofDigiMatchCollOut = new TClonesArray("CbmMatch", 100);
544 //rootMgr->Register( "TofCalDigi","Tof", fTofCalDigisCollOut, fbWriteDigisInOut);
545 rootMgr->RegisterAny("TofCalDigi", fTofCalDigiVecOut, fbWriteDigisInOut);
546 rootMgr->Register(tHitBranchName, "Tof", fTofHitsCollOut, fbWriteHitsInOut);
547 rootMgr->Register(tHitDigiMatchBranchName, "Tof", fTofDigiMatchCollOut, fbWriteHitsInOut);
548 }
549 LOG(info) << "out branches: " << tHitBranchName << ", " << tHitDigiMatchBranchName;
550 return kTRUE;
551}
553{
554
555 // Initialize the TOF GeoHandler
556 Bool_t isSimulation = kFALSE;
557 LOG(info) << "CbmTofHitMaker::InitParameters - Geometry, Mapping, ... ??";
558
559 // Get Base Container
560 FairRun* ana = FairRun::Instance();
561 FairRuntimeDb* rtdb = ana->GetRuntimeDb();
562
563 Int_t iGeoVersion = fGeoHandler->Init(isSimulation);
564 if (k14a > iGeoVersion) {
565 LOG(error) << "CbmTofHitMaker::InitParameters => Only compatible "
566 "with geometries after v14a !!!";
567 return kFALSE;
568 }
569 if (iGeoVersion == k14a)
571 else
573
574 // create digitization parameters from geometry file
575 CbmTofCreateDigiPar* tofDigiPar = new CbmTofCreateDigiPar("TOF Digi Producer", "TOF task");
576 LOG(info) << "Create DigiPar ";
577 tofDigiPar->Init();
578
579 fDigiPar = (CbmTofDigiPar*) (rtdb->getContainer("CbmTofDigiPar"));
580 if (0 == fDigiPar) {
581 LOG(error) << "CbmTofHitMaker::InitParameters => Could not obtain "
582 "the CbmTofDigiPar ";
583 return kFALSE;
584 }
585
586 fDigiBdfPar = (CbmTofDigiBdfPar*) (rtdb->getContainer("CbmTofDigiBdfPar"));
587 if (0 == fDigiBdfPar) {
588 LOG(error) << "CbmTofHitMaker::InitParameters => Could not obtain "
589 "the CbmTofDigiBdfPar ";
590 return kFALSE;
591 }
592
593 rtdb->initContainers(ana->GetRunId());
594
595 LOG(info) << "CbmTofHitMaker::InitParameter: currently " << fDigiPar->GetNrOfModules() << " digi cells ";
596
599
601 fdMaxTimeDist = fMaxTimeDist; // modify default
603 * 0.5; // cut consistently on positions (with default signal velocity)
604 }
605
606 LOG(info) << " BuildCluster with MaxTimeDist " << fdMaxTimeDist << ", MaxSpaceDist " << fdMaxSpaceDist;
607
608 if (fiCluMulMax == 0) fiCluMulMax = 100;
609 if (fOutHstFileName == "") {
610 fOutHstFileName = "./tofEventClust.hst.root";
611 }
612
613 LOG(info) << " Hst Output filename = " << fOutHstFileName;
614
615 LOG(info) << "<I> BeamRefType = " << fiBeamRefType << ", Sm " << fiBeamRefSm << ", Det " << fiBeamRefDet
616 << ", MulMax " << fiBeamRefMulMax;
617
618 return kTRUE;
619}
620/************************************************************************************/
622{
623 // dimension and initialize calib parameter
624 // Int_t iNbDet = fDigiBdfPar->GetNbDet(); (VF) not used
625 Int_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes();
626
627 if (fTotMean != 0.) fdTTotMean = fTotMean; // adjust target mean for TOT
628
629 fvCPTOff.resize(iNbSmTypes);
630 fvCPTotGain.resize(iNbSmTypes);
631 fvCPTotOff.resize(iNbSmTypes);
632 fvCPWalk.resize(iNbSmTypes);
633 fvCPDelTof.resize(iNbSmTypes);
634 for (Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
635 Int_t iNbSm = fDigiBdfPar->GetNbSm(iSmType);
636 Int_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
637 fvCPTOff[iSmType].resize(iNbSm * iNbRpc);
638 fvCPTotGain[iSmType].resize(iNbSm * iNbRpc);
639 fvCPTotOff[iSmType].resize(iNbSm * iNbRpc);
640 fvCPWalk[iSmType].resize(iNbSm * iNbRpc);
641 fvCPDelTof[iSmType].resize(iNbSm * iNbRpc);
642 for (Int_t iSm = 0; iSm < iNbSm; iSm++) {
643 for (Int_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
644 // LOG(info)<<Form(" fvCPDelTof resize for SmT %d, R %d, B %d ",iSmType,iNbSm*iNbRpc,nbClDelTofBinX)
645 // ;
646 fvCPDelTof[iSmType][iSm * iNbRpc + iRpc].resize(nbClDelTofBinX);
647 for (Int_t iBx = 0; iBx < nbClDelTofBinX; iBx++) {
648 // LOG(info)<<Form(" fvCPDelTof for SmT %d, R %d, B %d",iSmType,iSm*iNbRpc+iRpc,iBx);
649 fvCPDelTof[iSmType][iSm * iNbRpc + iRpc][iBx].resize(iNSel);
650 for (Int_t iSel = 0; iSel < iNSel; iSel++)
651 fvCPDelTof[iSmType][iSm * iNbRpc + iRpc][iBx][iSel] = 0.; // initialize
652 }
653
654 Int_t iNbChan = fDigiBdfPar->GetNbChan(iSmType, iRpc);
655 fvCPTOff[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
656 fvCPTotGain[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
657 fvCPTotOff[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
658 fvCPWalk[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
659 Int_t nbSide = 2 - fDigiBdfPar->GetChanType(iSmType, iRpc);
660 for (Int_t iCh = 0; iCh < iNbChan; iCh++) {
661 fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh].resize(nbSide);
662 fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh].resize(nbSide);
663 fvCPTotOff[iSmType][iSm * iNbRpc + iRpc][iCh].resize(nbSide);
664 fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh].resize(nbSide);
665 for (Int_t iSide = 0; iSide < nbSide; iSide++) {
666 fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][iSide] = 0.; //initialize
667 fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][iSide] = 1.; //initialize
668 fvCPTotOff[iSmType][iSm * iNbRpc + iRpc][iCh][iSide] = 0.; //initialize
669 fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][iSide].resize(nbClWalkBinX);
670 for (Int_t iWx = 0; iWx < nbClWalkBinX; iWx++) {
671 fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][iSide][iWx] = 0.;
672 }
673 }
674 }
675 }
676 }
677 }
678 LOG(info) << "CbmTofHitMaker::InitCalibParameter: defaults set";
679
681 // <= To prevent histos from being sucked in by the param file of the TRootManager!
682 TFile* oldFile = gFile;
683 TDirectory* oldDir = gDirectory;
684
685 if (0 < fCalMode) {
686 LOG(info) << "CbmTofHitMaker::InitCalibParameter: read histos from "
687 << "file " << fCalParFileName;
688
689 // read parameter from histos
690 if (fCalParFileName.IsNull()) return kTRUE;
691
692 fCalParFile = new TFile(fCalParFileName, "");
693 if (NULL == fCalParFile) {
694 LOG(fatal) << "CbmTofHitMaker::InitCalibParameter: "
695 << "file " << fCalParFileName << " does not exist!";
696 return kTRUE;
697 }
698 /*
699 gDirectory->Print();
700 fCalParFile->cd();
701 fCalParFile->ls();
702 */
703 for (Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
704 Int_t iNbSm = fDigiBdfPar->GetNbSm(iSmType);
705 Int_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
706 TProfile* hSvel = (TProfile*) gDirectory->FindObjectAny(Form("cl_SmT%01d_Svel", iSmType));
707
708 // copy Histo to memory
709 TDirectory* curdir = gDirectory;
710 if (NULL != hSvel) {
711 gDirectory->cd(oldDir->GetPath());
712 // TProfile *hSvelmem = (TProfile *)hSvel->Clone(); (VF) not used
713 gDirectory->cd(curdir->GetPath());
714 }
715 else {
716 LOG(info) << "Svel histogram not found for module type " << iSmType;
717 }
718
719 for (Int_t iPar = 0; iPar < 4; iPar++) {
720 TProfile* hFparcur = (TProfile*) gDirectory->FindObjectAny(Form("cl_SmT%01d_Fpar%1d", iSmType, iPar));
721 if (NULL != hFparcur) {
722 gDirectory->cd(oldDir->GetPath());
723 // TProfile *hFparmem = (TProfile *)hFparcur->Clone(); (VF) not used
724 gDirectory->cd(curdir->GetPath());
725 }
726 }
727
728 for (Int_t iSm = 0; iSm < iNbSm; iSm++)
729 for (Int_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
730
731 // update default parameter
732 if (NULL != hSvel) {
733 Double_t Vscal = 1.; //hSvel->GetBinContent(iSm*iNbRpc+iRpc+1);
734 if (Vscal == 0.) Vscal = 1.;
735 fDigiBdfPar->SetSigVel(iSmType, iSm, iRpc, fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc) * Vscal);
736 LOG(info) << "Modify " << iSmType << iSm << iRpc << " Svel by " << Vscal << " to "
737 << fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc);
738 }
739 TH2F* htempPos_pfx =
740 (TH2F*) gDirectory->FindObjectAny(Form("cl_CorSmT%01d_sm%03d_rpc%03d_Pos_pfx", iSmType, iSm, iRpc));
741 TH2F* htempTOff_pfx =
742 (TH2F*) gDirectory->FindObjectAny(Form("cl_CorSmT%01d_sm%03d_rpc%03d_TOff_pfx", iSmType, iSm, iRpc));
743 TH1D* htempTot_Mean =
744 (TH1D*) gDirectory->FindObjectAny(Form("cl_CorSmT%01d_sm%03d_rpc%03d_Tot_Mean", iSmType, iSm, iRpc));
745 TH1D* htempTot_Off =
746 (TH1D*) gDirectory->FindObjectAny(Form("cl_CorSmT%01d_sm%03d_rpc%03d_Tot_Off", iSmType, iSm, iRpc));
747 if (NULL != htempPos_pfx && NULL != htempTOff_pfx && NULL != htempTot_Mean && NULL != htempTot_Off) {
748 Int_t iNbCh = fDigiBdfPar->GetNbChan(iSmType, iRpc);
749 Int_t iNbinTot = htempTot_Mean->GetNbinsX();
750 for (Int_t iCh = 0; iCh < iNbCh; iCh++) {
751
752 for (Int_t iSide = 0; iSide < 2; iSide++) {
753 Double_t TotMean = htempTot_Mean->GetBinContent(iCh * 2 + 1 + iSide); //nh +1 empirical(?)
754 if (0.001 < TotMean) {
755 fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][iSide] *= fdTTotMean / TotMean;
756 }
757 fvCPTotOff[iSmType][iSm * iNbRpc + iRpc][iCh][iSide] = htempTot_Off->GetBinContent(iCh * 2 + 1 + iSide);
758 }
759
760 Double_t YMean = ((TProfile*) htempPos_pfx)->GetBinContent(iCh + 1);
761 Double_t TMean = ((TProfile*) htempTOff_pfx)->GetBinContent(iCh + 1);
762 //Double_t dTYOff=YMean/fDigiBdfPar->GetSignalSpeed() ;
763 Double_t dTYOff = YMean / fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc);
764 fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][0] += -dTYOff + TMean;
765 fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][1] += +dTYOff + TMean;
766
767 if (5 == iSmType || 8 == iSmType) { // for PAD counters
768 fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][1] = fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][0];
769 fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][1] = fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][0];
770 fvCPTotOff[iSmType][iSm * iNbRpc + iRpc][iCh][1] = fvCPTotOff[iSmType][iSm * iNbRpc + iRpc][iCh][0];
771 }
772
773 LOG(debug) << "CbmTofHitMaker::InitCalibParameter:"
774 << " SmT " << iSmType << " Sm " << iSm << " Rpc " << iRpc << " Ch " << iCh
775 << Form(": YMean %f, TMean %f", YMean, TMean) << " -> "
776 << Form(" %f, %f, %f, %f ", fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][0],
777 fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][1],
778 fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][0],
779 fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][1])
780 << ", NbinTot " << iNbinTot;
781
782 TH1D* htempWalk0 = (TH1D*) gDirectory->FindObjectAny(
783 Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S0_Walk_px", iSmType, iSm, iRpc, iCh));
784 TH1D* htempWalk1 = (TH1D*) gDirectory->FindObjectAny(
785 Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S1_Walk_px", iSmType, iSm, iRpc, iCh));
786 if (NULL != htempWalk0 && NULL != htempWalk1) { // reinitialize Walk array
787 LOG(debug) << "Initialize Walk correction for "
788 << Form(" SmT%01d_sm%03d_rpc%03d_Ch%03d", iSmType, iSm, iRpc, iCh);
789 if (htempWalk0->GetNbinsX() != nbClWalkBinX)
790 LOG(error) << "CbmTofHitMaker::InitCalibParameter: "
791 "Inconsistent Walk histograms";
792 for (Int_t iBin = 0; iBin < nbClWalkBinX; iBin++) {
793 fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][0][iBin] = htempWalk0->GetBinContent(iBin + 1);
794 fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][1][iBin] = htempWalk1->GetBinContent(iBin + 1);
795 if (iCh == 5 && iBin == 10) // debugging
796 LOG(info) << Form("Read New SmT%01d_sm%03d_rpc%03d_Ch%03d bin %d walk %f ", iSmType, iSm, iRpc, iCh,
797 iBin, fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][0][iBin]);
798 if (5 == iSmType || 8 == iSmType) { // Pad structure
799 fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][1][iBin] =
800 fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][0][iBin];
801 }
802 }
803 }
804 else {
805 LOG(info) << "No Walk histograms for TSRC " << iSmType << iSm << iRpc << iCh;
806 }
807 }
808 }
809 else {
810 LOG(warning) << " Calibration histos " << Form("cl_SmT%01d_sm%03d_rpc%03d_XXX", iSmType, iSm, iRpc)
811 << " not found. ";
812 }
813 for (Int_t iSel = 0; iSel < iNSel; iSel++) {
814 TH1D* htmpDelTof = (TH1D*) gDirectory->FindObjectAny(
815 Form("cl_CorSmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof", iSmType, iSm, iRpc, iSel));
816 if (NULL == htmpDelTof) {
817 LOG(debug) << " Histos " << Form("cl_CorSmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof", iSmType, iSm, iRpc, iSel)
818 << " not found. ";
819 continue;
820 }
821 LOG(debug) << " Load DelTof from histos "
822 << Form("cl_CorSmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof", iSmType, iSm, iRpc, iSel) << ".";
823 for (Int_t iBx = 0; iBx < nbClDelTofBinX; iBx++) {
824 fvCPDelTof[iSmType][iSm * iNbRpc + iRpc][iBx][iSel] += htmpDelTof->GetBinContent(iBx + 1);
825 }
826
827 // copy Histo to memory
828 // TDirectory * curdir = gDirectory;
829 gDirectory->cd(oldDir->GetPath());
830 TH1D* h1DelTof =
831 (TH1D*) htmpDelTof->Clone(Form("cl_CorSmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof", iSmType, iSm, iRpc, iSel));
832
833 LOG(debug) << " copy histo " << h1DelTof->GetName() << " to directory " << oldDir->GetName();
834
835 gDirectory->cd(curdir->GetPath());
836 }
837 }
838 }
839 }
840 // fCalParFile->Delete();
842 // <= To prevent histos from being sucked in by the param file of the TRootManager!
843 gFile = oldFile;
844 gDirectory = oldDir;
845 LOG(info) << "CbmTofHitMaker::InitCalibParameter: initialization done";
846 return kTRUE;
847}
848/************************************************************************************/
850{
851 LOG(info) << "CbmTofHitMaker::LoadGeometry starting for " << fDigiBdfPar->GetNbDet() << " described detectors, "
852 << fDigiPar->GetNrOfModules() << " geometrically known cells ";
853
854 Int_t iNrOfCells = fDigiPar->GetNrOfModules();
855 LOG(info) << "Digi Parameter container contains " << iNrOfCells << " cells.";
856 for (Int_t icell = 0; icell < iNrOfCells; ++icell) {
857
858 Int_t cellId = fDigiPar->GetCellId(icell); // cellId is assigned in CbmTofCreateDigiPar
859 fChannelInfo = fDigiPar->GetCell(cellId);
860
861 Int_t smtype = fGeoHandler->GetSMType(cellId);
862 Int_t smodule = fGeoHandler->GetSModule(cellId);
863 Int_t module = fGeoHandler->GetCounter(cellId);
864 Int_t cell = fGeoHandler->GetCell(cellId);
865
866 Double_t x = fChannelInfo->GetX();
867 Double_t y = fChannelInfo->GetY();
868 Double_t z = fChannelInfo->GetZ();
869 Double_t dx = fChannelInfo->GetSizex();
870 Double_t dy = fChannelInfo->GetSizey();
871 LOG(debug) << "-I- InitPar " << icell << " Id: " << Form("0x%08x", cellId) << " " << cell << " tmcs: " << smtype
872 << " " << smodule << " " << module << " " << cell << " x=" << Form("%6.2f", x)
873 << " y=" << Form("%6.2f", y) << " z=" << Form("%6.2f", z) << " dx=" << dx << " dy=" << dy;
874
875 TGeoNode* fNode = // prepare local->global trafo
876 gGeoManager->FindNode(fChannelInfo->GetX(), fChannelInfo->GetY(), fChannelInfo->GetZ());
877 LOG(debug2) << Form(" Node at (%6.1f,%6.1f,%6.1f) : 0x%p", fChannelInfo->GetX(), fChannelInfo->GetY(),
878 fChannelInfo->GetZ(), fNode);
879 if (icell == 0) {
880 TGeoHMatrix* cMatrix = gGeoManager->GetCurrentMatrix();
881 fNode->Print();
882 fDigiPar->GetNode(cellId)->Print();
883 cMatrix->Print();
884 }
885 }
886
887 Int_t iNbDet = fDigiBdfPar->GetNbDet();
888 fvDeadStrips.resize(iNbDet);
889 fvTimeLastDigi.resize(iNbDet);
890 fvTimeFirstDigi.resize(iNbDet);
891 fvMulDigi.resize(iNbDet);
892
893 for (Int_t iDetIndx = 0; iDetIndx < iNbDet; iDetIndx++) {
894 Int_t iUniqueId = fDigiBdfPar->GetDetUId(iDetIndx);
895 Int_t iSmType = CbmTofAddress::GetSmType(iUniqueId);
896 Int_t iSmId = CbmTofAddress::GetSmId(iUniqueId);
897 Int_t iRpcId = CbmTofAddress::GetRpcId(iUniqueId);
898 LOG(info) << " DetIndx " << iDetIndx << "(" << iNbDet << "), SmType " << iSmType << ", SmId " << iSmId << ", RpcId "
899 << iRpcId << " => UniqueId " << Form("0x%08x ", iUniqueId)
900 << Form(" Svel %6.6f, DeadStrips 0x%08x ", fDigiBdfPar->GetSigVel(iSmType, iSmId, iRpcId),
901 fvDeadStrips[iDetIndx]);
902 Int_t iNbChan = fDigiBdfPar->GetNbChan(iSmType, iRpcId);
903 fvTimeLastDigi[iDetIndx].resize(iNbChan * 2);
904 for (Int_t iCh = 0; iCh < iNbChan * 2; iCh++)
905 fvTimeLastDigi[iDetIndx][iCh] = 0.;
906 fvTimeFirstDigi[iDetIndx].resize(iNbChan * 2);
907 fvMulDigi[iDetIndx].resize(iNbChan * 2);
908 for (Int_t iCh = 0; iCh < iNbChan * 2; iCh++) {
909 fvTimeFirstDigi[iDetIndx][iCh] = 0.;
910 fvMulDigi[iDetIndx][iCh] = 0.;
911 }
912
913 Int_t iCell = -1;
914 while (kTRUE) {
915 Int_t iUCellId = CbmTofAddress::GetUniqueAddress(iSmId, iRpcId, ++iCell, 0, iSmType);
916 fChannelInfo = fDigiPar->GetCell(iUCellId);
917 if (NULL == fChannelInfo) break;
918 LOG(debug3) << " Cell " << iCell << Form(" 0x%08x ", iUCellId) << Form(", fCh 0x%p ", fChannelInfo)
919 << ", x: " << fChannelInfo->GetX() << ", y: " << fChannelInfo->GetY()
920 << ", z: " << fChannelInfo->GetZ();
921 }
922 }
923
924 // return kTRUE;
925
926 Int_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes();
927
928 if (kTRUE == fDigiBdfPar->UseExpandedDigi()) {
929 fStorDigi.resize(iNbSmTypes);
930 fStorDigiInd.resize(iNbSmTypes);
931 fviClusterSize.resize(iNbSmTypes);
932 fviTrkMul.resize(iNbSmTypes);
933 fvdX.resize(iNbSmTypes);
934 fvdY.resize(iNbSmTypes);
935 fvdDifX.resize(iNbSmTypes);
936 fvdDifY.resize(iNbSmTypes);
937 fvdDifCh.resize(iNbSmTypes);
938 fviClusterMul.resize(iNbSmTypes);
939 fvLastHits.resize(iNbSmTypes);
940
941 for (Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
942 Int_t iNbSm = fDigiBdfPar->GetNbSm(iSmType);
943 Int_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
944 fStorDigi[iSmType].resize(iNbSm * iNbRpc);
945 fStorDigiInd[iSmType].resize(iNbSm * iNbRpc);
946 fviClusterSize[iSmType].resize(iNbRpc);
947 fviTrkMul[iSmType].resize(iNbRpc);
948 fvdX[iSmType].resize(iNbRpc);
949 fvdY[iSmType].resize(iNbRpc);
950 fvdDifX[iSmType].resize(iNbRpc);
951 fvdDifY[iSmType].resize(iNbRpc);
952 fvdDifCh[iSmType].resize(iNbRpc);
953 for (Int_t iSm = 0; iSm < iNbSm; iSm++) {
954 fviClusterMul[iSmType].resize(iNbSm);
955 fvLastHits[iSmType].resize(iNbSm);
956 for (Int_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
957 fviClusterMul[iSmType][iSm].resize(iNbRpc);
958 fvLastHits[iSmType][iSm].resize(iNbRpc);
959 Int_t iNbChan = fDigiBdfPar->GetNbChan(iSmType, iRpc);
960 if (iNbChan == 0) {
961 LOG(warning) << "CbmTofHitMaker::LoadGeometry: StoreDigi "
962 "without channels "
963 << Form("SmTy %3d, Sm %3d, NbRpc %3d, Rpc, %3d ", iSmType, iSm, iNbRpc, iRpc);
964 }
965 LOG(debug1) << "CbmTofHitMaker::LoadGeometry: StoreDigi with "
966 << Form(" %3d %3d %3d %3d %5d ", iSmType, iSm, iNbRpc, iRpc, iNbChan);
967 fStorDigi[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
968 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
969 fvLastHits[iSmType][iSm][iRpc].resize(iNbChan);
970 } // for( Int_t iRpc = 0; iRpc < iNbRpc; iRpc++ )
971 } // for( Int_t iSm = 0; iSm < iNbSm; iSm++ )
972 } // for( Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++ )
973 } // if( kTRUE == fDigiBdfPar->UseExpandedDigi() )
974
975 return kTRUE;
976}
978{
979 LOG(info) << "CbmTofHitMaker::DeleteGeometry starting";
980 return kTRUE;
981 Int_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes();
982 if (kTRUE == fDigiBdfPar->UseExpandedDigi()) {
983 for (Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
984 Int_t iNbSm = fDigiBdfPar->GetNbSm(iSmType);
985 Int_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
986 for (Int_t iSm = 0; iSm < iNbSm; iSm++) {
987 for (Int_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
988 fStorDigi[iSmType][iSm * iNbRpc + iRpc].clear();
989 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc].clear();
990 }
991 } // for( Int_t iSm = 0; iSm < iNbSm; iSm++ )
992 fStorDigi[iSmType].clear();
993 fStorDigiInd[iSmType].clear();
994 } // for( Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++ )
995 fStorDigi.clear();
996 fStorDigiInd.clear();
997 } // if( kTRUE == fDigiBdfPar->UseExpandedDigi() )
998 return kTRUE;
999}
1000/************************************************************************************/
1001// Histogramming functions
1003{
1004 // <= To prevent histos from being sucked in by the param file of the TRootManager!
1005 TDirectory* oldir = gDirectory;
1006 gROOT->cd();
1008 new TH1I("TofClustBuildTime", "Time needed to build clusters in each event; Time [s]", 100, 0.0, 4.0);
1009 gDirectory->cd(oldir->GetPath());
1010 // <= To prevent histos from being sucked in by the param file of the TRootManager!
1011
1012 return kTRUE;
1013}
1014
1016{
1017 fhClustBuildTime->Fill(fStop.GetSec() - fStart.GetSec() + (fStop.GetNanoSec() - fStart.GetNanoSec()) / 1e9);
1018 return kTRUE;
1019}
1020
1022{
1024 TFile* oldFile = gFile;
1025 TDirectory* oldDir = gDirectory;
1026
1027 TFile* fHist;
1028 fHist = new TFile(fOutHstFileName, "RECREATE");
1029 fHist->cd();
1030 fhClustBuildTime->Write();
1031
1033 gFile = oldFile;
1034 gDirectory = oldDir;
1035
1036 fHist->Close();
1037
1038 return kTRUE;
1039}
1041{
1042 delete fhClustBuildTime;
1043 return kTRUE;
1044}
1045/************************************************************************************/
1047{
1048 //gGeoManager->SetTopVolume( gGeoManager->FindVolumeFast("tof_v14a") );
1049 gGeoManager->CdTop();
1050
1051 // if(NULL == fTofDigisColl) {
1052 if (fTofDigiVec.empty()) {
1053 LOG(info) << " No RawDigis defined ! Check! ";
1054 return kFALSE;
1055 }
1056 fiNevtBuild++;
1057 LOG(debug) << "Build clusters from "
1058 // <<fTofDigisColl->GetEntriesFast()<<" digis in event "<<fiNevtBuild;
1059 << fTofDigiVec.size() << " digis in event " << fiNevtBuild;
1060
1061 fTRefHits = 0.;
1062
1063 Int_t iNbTofDigi = fTofDigiVec.size();
1064 //Int_t iNbTofDigi = fTofDigisColl->GetEntriesFast();
1065 if (iNbTofDigi > 100000) {
1066 LOG(warning) << "Too many digis in event " << fiNevtBuild;
1067 return kFALSE;
1068 }
1070 // Duplicate type "5" - digis
1071 // Int_t iNbDigi=iNbTofDigi;
1072 for (Int_t iDigInd = 0; iDigInd < iNbTofDigi; iDigInd++) {
1073 CbmTofDigi* pDigi = &(fTofDigiVec.at(iDigInd));
1074 //CbmTofDigi *pDigi = (CbmTofDigi*) fTofDigisColl->At( iDigInd );
1075 if (pDigi->GetType() == 5) {
1076 if (pDigi->GetSide() == 1) {
1077 bAddBeamCounterSideDigi = kFALSE; // disable for current data set
1078 LOG(info) << "Start counter digi duplication disabled";
1079 break;
1080 }
1081 fTofDigiVec.push_back(CbmTofDigi(*pDigi));
1082 CbmTofDigi* pDigiN = &(fTofDigiVec.back());
1083 // CbmTofDigi *pDigiN = new((*fTofDigisColl)[iNbDigi++]) CbmTofDigi( *pDigi );
1084 pDigiN->SetAddress(pDigi->GetSm(), pDigi->GetRpc(), pDigi->GetChannel(), (0 == pDigi->GetSide()) ? 1 : 0,
1085 pDigi->GetType());
1086 LOG(debug) << "Duplicated digi " << fTofDigiVec.size() << " with address 0x" << std::hex
1087 << pDigiN->GetAddress();
1088 }
1089 }
1090 iNbTofDigi = fTofDigiVec.size();
1091 //iNbTofDigi = fTofDigisColl->GetEntriesFast(); // Update
1092 }
1093
1094 if (kTRUE) {
1095 for (UInt_t iDetIndx = 0; iDetIndx < fvTimeFirstDigi.size(); iDetIndx++)
1096 for (UInt_t iCh = 0; iCh < fvTimeFirstDigi[iDetIndx].size(); iCh++) {
1097 fvTimeFirstDigi[iDetIndx][iCh] = 0.;
1098 fvMulDigi[iDetIndx][iCh] = 0.;
1099 }
1100
1101 for (Int_t iDigInd = 0; iDigInd < iNbTofDigi; iDigInd++) {
1102 //CbmTofDigi *pDigi = (CbmTofDigi*) fTofDigisColl->At( iDigInd );
1103 CbmTofDigi* pDigi = &(fTofDigiVec.at(iDigInd));
1104 Int_t iDetIndx = fDigiBdfPar->GetDetInd(pDigi->GetAddress());
1105
1106 LOG(debug) << iDigInd << " " << pDigi << Form(" Address : 0x%08x ", pDigi->GetAddress()) << " SmT "
1107 << pDigi->GetType() << " Sm " << pDigi->GetSm() << " Rpc " << pDigi->GetRpc() << " Ch "
1108 << pDigi->GetChannel() << " S " << pDigi->GetSide() << ", DetIndx " << iDetIndx << " : "
1109 << pDigi->ToString()
1110 // <<" Time "<<pDigi->GetTime()
1111 // <<" Tot " <<pDigi->GetTot()
1112 ;
1113
1114 if (fDigiBdfPar->GetNbDet() - 1 < iDetIndx || iDetIndx < 0) {
1115 LOG(debug) << Form(" Wrong DetIndx %d >< %d ", iDetIndx, fDigiBdfPar->GetNbDet());
1116 break;
1117 }
1118
1119 size_t iDigiCh = pDigi->GetChannel() * 2 + pDigi->GetSide();
1120 if (iDigiCh < fvTimeLastDigi[iDetIndx].size()) {
1121 fvTimeLastDigi[iDetIndx][iDigiCh] = pDigi->GetTime();
1122
1123 if (fvTimeFirstDigi[iDetIndx][iDigiCh] != 0.) {
1124 fvMulDigi[iDetIndx][iDigiCh]++;
1125 }
1126 else {
1127 fvTimeFirstDigi[iDetIndx][iDigiCh] = pDigi->GetTime();
1128 fvMulDigi[iDetIndx][iDigiCh]++;
1129 }
1130 }
1131 }
1132 for (UInt_t iDetIndx = 0; iDetIndx < fvTimeFirstDigi.size(); iDetIndx++)
1133 for (UInt_t iCh = 0; iCh < fvTimeFirstDigi[iDetIndx].size(); iCh++) {
1134 if (fvTimeFirstDigi[iDetIndx][iCh] != 0.) fhRpcDigiDTMul[iDetIndx]->Fill(iCh, fvMulDigi[iDetIndx][iCh]);
1135 }
1136 } // kTRUE end
1137
1138 // Calibrate RawDigis
1139 if (kTRUE == fDigiBdfPar->UseExpandedDigi()) {
1140 CbmTofDigi* pDigi;
1141 // CbmTofDigi *pCalDigi=NULL; (VF) not used
1142 CalibRawDigis();
1143
1144 // Then loop over the digis array and store the Digis in separate vectors for
1145 // each RPC modules
1146
1147 // iNbTofDigi = fTofCalDigisColl->GetEntriesFast();
1148 iNbTofDigi = fTofCalDigiVec->size();
1149 for (Int_t iDigInd = 0; iDigInd < iNbTofDigi; iDigInd++) {
1150 // pDigi = (CbmTofDigi*) fTofCalDigisColl->At( iDigInd );
1151 pDigi = &(fTofCalDigiVec->at(iDigInd));
1152 LOG(debug1) << "AC " // After Calibration
1153 << Form("0x%08x", pDigi->GetAddress()) << " TSRC " << pDigi->GetType() << pDigi->GetSm()
1154 << pDigi->GetRpc() << Form("%2d", (Int_t) pDigi->GetChannel()) << " " << pDigi->GetSide() << " "
1155 << Form("%f", pDigi->GetTime()) << " " << pDigi->GetTot();
1156
1157 if (fDigiBdfPar->GetNbSmTypes() > pDigi->GetType() // prevent crash due to misconfiguration
1158 && fDigiBdfPar->GetNbSm(pDigi->GetType()) > pDigi->GetSm()
1159 && fDigiBdfPar->GetNbRpc(pDigi->GetType()) > pDigi->GetRpc()
1160 && fDigiBdfPar->GetNbChan(pDigi->GetType(), pDigi->GetRpc()) > pDigi->GetChannel()) {
1161 fStorDigi[pDigi->GetType()][pDigi->GetSm() * fDigiBdfPar->GetNbRpc(pDigi->GetType()) + pDigi->GetRpc()]
1162 [pDigi->GetChannel()]
1163 .push_back(pDigi);
1164 fStorDigiInd[pDigi->GetType()][pDigi->GetSm() * fDigiBdfPar->GetNbRpc(pDigi->GetType()) + pDigi->GetRpc()]
1165 [pDigi->GetChannel()]
1166 .push_back(iDigInd);
1167 }
1168 else {
1169 LOG(info) << "Skip2 Digi "
1170 << " Type " << pDigi->GetType() << " " << fDigiBdfPar->GetNbSmTypes() << " Sm " << pDigi->GetSm()
1171 << " " << fDigiBdfPar->GetNbSm(pDigi->GetType()) << " Rpc " << pDigi->GetRpc() << " "
1172 << fDigiBdfPar->GetNbRpc(pDigi->GetType()) << " Ch " << pDigi->GetChannel() << " "
1173 << fDigiBdfPar->GetNbChan(pDigi->GetType(), 0);
1174 }
1175 } // for( Int_t iDigInd = 0; iDigInd < nTofDigi; iDigInd++ )
1176
1177 // inspect digi array
1178 Int_t iNbDet = fDigiBdfPar->GetNbDet();
1179 for (Int_t iDetIndx = 0; iDetIndx < iNbDet; iDetIndx++) {
1180 Int_t iDetId = fviDetId[iDetIndx];
1181 Int_t iSmType = CbmTofAddress::GetSmType(iDetId);
1182 Int_t iSm = CbmTofAddress::GetSmId(iDetId);
1183 Int_t iRpc = CbmTofAddress::GetRpcId(iDetId);
1184 Int_t iNbStrips = fDigiBdfPar->GetNbChan(iSmType, iRpc);
1185 for (Int_t iStrip = 0; iStrip < iNbStrips; iStrip++) {
1186 Int_t iDigiMul = fStorDigi[iSmType][iSm * fDigiBdfPar->GetNbRpc(iSmType) + iRpc][iStrip].size();
1187 //LOG(info)<<"Inspect TSRC "<<iSmType<<iSm<<iRpc<<iStrip<<" with "<<iNbStrips<<" strips: Mul "<<iDigiMul;
1188 if (iDigiMul > 0) {
1189 fhRpcDigiMul[iDetIndx]->Fill(iStrip, iDigiMul);
1190 if (iDigiMul == 1) {
1191 fhRpcDigiStatus[iDetIndx]->Fill(iStrip, 0);
1192 if (iStrip > 0)
1193 if (fStorDigi[iSmType][iSm * fDigiBdfPar->GetNbRpc(iSmType) + iRpc][iStrip - 1].size() > 1) {
1194 fhRpcDigiStatus[iDetIndx]->Fill(iStrip, 1);
1195 if (TMath::Abs(
1196 fStorDigi[iSmType][iSm * fDigiBdfPar->GetNbRpc(iSmType) + iRpc][iStrip][0]->GetTime()
1197 - fStorDigi[iSmType][iSm * fDigiBdfPar->GetNbRpc(iSmType) + iRpc][iStrip - 1][0]->GetTime())
1198 < fMaxTimeDist)
1199 fhRpcDigiStatus[iDetIndx]->Fill(iStrip, 3);
1200 }
1201 if (iStrip < iNbStrips - 2) {
1202 if (fStorDigi[iSmType][iSm * fDigiBdfPar->GetNbRpc(iSmType) + iRpc][iStrip + 1].size() > 1) {
1203 fhRpcDigiStatus[iDetIndx]->Fill(iStrip, 2);
1204 if (TMath::Abs(
1205 fStorDigi[iSmType][iSm * fDigiBdfPar->GetNbRpc(iSmType) + iRpc][iStrip][0]->GetTime()
1206 - fStorDigi[iSmType][iSm * fDigiBdfPar->GetNbRpc(iSmType) + iRpc][iStrip + 1][0]->GetTime())
1207 < fMaxTimeDist)
1208 fhRpcDigiStatus[iDetIndx]->Fill(iStrip, 4);
1209 }
1210 }
1211 }
1212 }
1213 }
1214 }
1215
1216 BuildHits();
1217
1218 } // if( kTRUE == fDigiBdfPar->UseExpandedDigi() )
1219 else {
1220 LOG(error) << " Compressed Digis not implemented ... ";
1221 return kFALSE; // not implemented properly yet
1222 }
1223 return kTRUE;
1224}
1225
1227{
1228 // Merge clusters from neigbouring Rpc within a (Super)Module
1229 if (NULL == fTofHitsColl) {
1230 LOG(info) << " No Hits defined ! Check! ";
1231 return kFALSE;
1232 }
1233 // inspect hits
1234 for (Int_t iHitInd = 0; iHitInd < fTofHitsColl->GetEntriesFast(); iHitInd++) {
1235 CbmTofHit* pHit = (CbmTofHit*) fTofHitsColl->At(iHitInd);
1236 if (NULL == pHit) continue;
1237
1238 Int_t iDetId = (pHit->GetAddress() & DetMask);
1239 Int_t iSmType = CbmTofAddress::GetSmType(iDetId);
1240 Int_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
1241 if (iSmType != 5 && iSmType != 8) continue; // only merge diamonds and Pad
1242 LOG(debug) << "MergeClusters: in SmT " << iSmType << " for " << iNbRpc << " Rpcs";
1243
1244 if (iNbRpc > 1) { // check for possible mergers
1245 Int_t iSm = CbmTofAddress::GetSmId(iDetId);
1246 Int_t iRpc = CbmTofAddress::GetRpcId(iDetId);
1247 Int_t iChId = pHit->GetAddress();
1248 fChannelInfo = fDigiPar->GetCell(iChId);
1249 Int_t iCh = CbmTofAddress::GetChannelId(iChId);
1250 LOG(debug) << "MergeClusters: Check for mergers in "
1251 << Form(" SmT %d, Sm %d, Rpc %d, Ch %d - hit %d", iSmType, iSm, iRpc, iCh, iHitInd);
1252 for (Int_t iHitInd2 = iHitInd + 1; iHitInd2 < fTofHitsColl->GetEntriesFast(); iHitInd2++) {
1253 CbmTofHit* pHit2 = (CbmTofHit*) fTofHitsColl->At(iHitInd2);
1254 if (NULL == pHit2) continue;
1255 Int_t iDetId2 = (pHit2->GetAddress() & DetMask);
1256 Int_t iSmType2 = CbmTofAddress::GetSmType(iDetId2);
1257 if (iSmType2 == iSmType) {
1258 Int_t iSm2 = CbmTofAddress::GetSmId(iDetId2);
1259 if (iSm2 == iSm || iSmType == 5) {
1260 Int_t iRpc2 = CbmTofAddress::GetRpcId(iDetId2);
1261 if (TMath::Abs(iRpc - iRpc2) == 1 || iSm2 != iSm) { // Found neighbour
1262 Int_t iChId2 = pHit2->GetAddress();
1263 // CbmTofCell *fChannelInfo2 = fDigiPar->GetCell( iChId2 ); (VF) not used
1264 Int_t iCh2 = CbmTofAddress::GetChannelId(iChId2);
1265 Double_t xPos = pHit->GetX();
1266 Double_t yPos = pHit->GetY();
1267 Double_t tof = pHit->GetTime();
1268 Double_t xPos2 = pHit2->GetX();
1269 Double_t yPos2 = pHit2->GetY();
1270 Double_t tof2 = pHit2->GetTime();
1271 LOG(debug) << "MergeClusters: Found hit in neighbour "
1272 << Form(" SmT %d, Sm %d, Rpc %d, Ch %d - hit %d", iSmType2, iSm2, iRpc2, iCh2, iHitInd2)
1273 << Form(" DX %6.1f, DY %6.1f, DT %6.1f", xPos - xPos2, yPos - yPos2, tof - tof2);
1274
1275 if (TMath::Abs(xPos - xPos2) < fdCaldXdYMax * 2. && TMath::Abs(yPos - yPos2) < fdCaldXdYMax * 2.
1276 && TMath::Abs(tof - tof2) < fMaxTimeDist) {
1277
1278 CbmMatch* digiMatch = (CbmMatch*) fTofDigiMatchColl->At(iHitInd);
1279 Double_t dTot = 0;
1280 for (Int_t iLink = 0; iLink < digiMatch->GetNofLinks(); iLink += 2) { // loop over digis
1281 CbmLink L0 = digiMatch->GetLink(iLink);
1282 UInt_t iDigInd0 = L0.GetIndex();
1283 UInt_t iDigInd1 = (digiMatch->GetLink(iLink + 1)).GetIndex();
1284 // if (iDigInd0 < fTofCalDigisColl->GetEntriesFast() && iDigInd1 < fTofCalDigisColl->GetEntriesFast()){
1285 if (iDigInd0 < fTofCalDigiVec->size() && iDigInd1 < fTofCalDigiVec->size()) {
1286 // CbmTofDigi *pDig0 = (CbmTofDigi*) (fTofCalDigisColl->At(iDigInd0));
1287 // CbmTofDigi *pDig1 = (CbmTofDigi*) (fTofCalDigisColl->At(iDigInd1));
1288 CbmTofDigi* pDig0 = &(fTofCalDigiVec->at(iDigInd0));
1289 CbmTofDigi* pDig1 = &(fTofCalDigiVec->at(iDigInd1));
1290 dTot += pDig0->GetTot();
1291 dTot += pDig1->GetTot();
1292 }
1293 }
1294
1295 CbmMatch* digiMatch2 = (CbmMatch*) fTofDigiMatchColl->At(iHitInd2);
1296 Double_t dTot2 = 0;
1297 for (Int_t iLink = 0; iLink < digiMatch2->GetNofLinks(); iLink += 2) { // loop over digis
1298 CbmLink L0 = digiMatch2->GetLink(iLink);
1299 UInt_t iDigInd0 = L0.GetIndex();
1300 UInt_t iDigInd1 = (digiMatch2->GetLink(iLink + 1)).GetIndex();
1301 // if (iDigInd0 < fTofCalDigisColl->GetEntriesFast() && iDigInd1 < fTofCalDigisColl->GetEntriesFast()){
1302 if (iDigInd0 < fTofCalDigiVec->size() && iDigInd1 < fTofCalDigiVec->size()) {
1303 // CbmTofDigi *pDig0 = (CbmTofDigi*) (fTofCalDigisColl->At(iDigInd0));
1304 // CbmTofDigi *pDig1 = (CbmTofDigi*) (fTofCalDigisColl->At(iDigInd1));
1305 CbmTofDigi* pDig0 = &(fTofCalDigiVec->at(iDigInd0));
1306 CbmTofDigi* pDig1 = &(fTofCalDigiVec->at(iDigInd1));
1307 dTot2 += pDig0->GetTot();
1308 dTot2 += pDig1->GetTot();
1309 digiMatch->AddLink(CbmLink(pDig0->GetTot(), iDigInd0, fiOutputTreeEntry, fiFileIndex));
1310 digiMatch->AddLink(CbmLink(pDig1->GetTot(), iDigInd1, fiOutputTreeEntry, fiFileIndex));
1311 }
1312 }
1313 LOG(debug) << "MergeClusters: Found merger in neighbour "
1314 << Form(" SmT %d, Sm %d, Rpc %d, Ch %d - hit %d(%d)", iSmType2, iSm2, iRpc2, iCh2, iHitInd2,
1315 fTofHitsColl->GetEntriesFast())
1316 << Form(" DX %6.1f, DY %6.1f, DT %6.1f", xPos - xPos2, yPos - yPos2, tof - tof2)
1317 << Form(" Tots %6.1f - %6.1f", dTot, dTot2);
1318 Double_t dTotSum = dTot + dTot2;
1319 Double_t dxPosM = (xPos * dTot + xPos2 * dTot2) / dTotSum;
1320 Double_t dyPosM = (yPos * dTot + yPos2 * dTot2) / dTotSum;
1321 Double_t dtofM = (tof * dTot + tof2 * dTot2) / dTotSum;
1322 pHit->SetX(dxPosM);
1323 pHit->SetY(dyPosM);
1324 pHit->SetTime(dtofM);
1325
1326 // remove merged hit at iHitInd2 and update digiMatch
1327
1328 fTofHitsColl->RemoveAt(iHitInd2);
1329 fTofDigiMatchColl->RemoveAt(iHitInd2);
1330 fTofDigiMatchColl->Compress();
1331 fTofHitsColl->Compress();
1332 LOG(debug) << "MergeClusters: Compress TClonesArrays to " << fTofHitsColl->GetEntriesFast() << ", "
1333 << fTofDigiMatchColl->GetEntriesFast();
1334 /*
1335 for(Int_t i=iHitInd2; i<fTofHitsColl->GetEntriesFast(); i++){ // update RefLinks
1336 CbmTofHit *pHiti = (CbmTofHit*) fTofHitsColl->At( i );
1337 pHiti->SetRefId(i);
1338 }
1339 */
1340 //check merged hit (cluster)
1341 //pHit->Print();
1342 }
1343 }
1344 }
1345 }
1346 }
1347 }
1348 }
1349 return kTRUE;
1350}
1351
1352static Double_t f1_xboxe(double* x, double* par)
1353{
1354 double xx = x[0];
1355 double wx = 1. - par[4] * TMath::Power(xx + par[5], 2);
1356 double xboxe = par[0] * 0.25 * (1. + TMath::Erf((xx + par[1] - par[3]) / par[2]))
1357 * (1. + TMath::Erf((-xx + par[1] + par[3]) / par[2]));
1358 return xboxe * wx;
1359}
1360
1361void CbmTofHitMaker::fit_ybox(const char* hname)
1362{
1363 TH1* h1;
1364 h1 = (TH1*) gROOT->FindObjectAny(hname);
1365 if (NULL != h1) {
1366 fit_ybox(h1, 0.);
1367 }
1368}
1369
1370void CbmTofHitMaker::fit_ybox(TH1* h1, Double_t ysize)
1371{
1372 Double_t* fpar = NULL;
1373 fit_ybox(h1, ysize, fpar);
1374}
1375
1376void CbmTofHitMaker::fit_ybox(TH1* h1, Double_t ysize, Double_t* fpar = NULL)
1377{
1378 TAxis* xaxis = h1->GetXaxis();
1379 Double_t Ymin = xaxis->GetXmin();
1380 Double_t Ymax = xaxis->GetXmax();
1381 TF1* f1 = new TF1("YBox", f1_xboxe, Ymin, Ymax, 6);
1382 Double_t yini = (h1->GetMaximum() + h1->GetMinimum()) * 0.5;
1383 if (ysize == 0.) ysize = Ymax * 0.8;
1384 f1->SetParameters(yini, ysize * 0.5, 1., 0., 0., 0.);
1385 // f1->SetParLimits(1,ysize*0.8,ysize*1.2);
1386 f1->SetParLimits(2, 0.2, 3.);
1387 f1->SetParLimits(3, -4., 4.);
1388 if (fpar != NULL) {
1389 Double_t fp[4];
1390 for (Int_t i = 0; i < 4; i++)
1391 fp[i] = *fpar++;
1392 for (Int_t i = 0; i < 4; i++)
1393 f1->SetParameter(2 + i, fp[i]);
1394 LOG(debug) << "Ini Fpar for " << h1->GetName() << " with "
1395 << Form(" %6.3f %6.3f %6.3f %6.3f ", fp[0], fp[1], fp[2], fp[3]);
1396 }
1397
1398 h1->Fit("YBox", "QM0");
1399
1400 double res[10];
1401 double err[10];
1402 res[9] = f1->GetChisquare();
1403
1404 for (int i = 0; i < 6; i++) {
1405 res[i] = f1->GetParameter(i);
1406 err[i] = f1->GetParError(i);
1407 //cout << " FPar "<< i << ": " << res[i] << ", " << err[i] << endl;
1408 }
1409 LOG(debug) << "YBox Fit of " << h1->GetName() << " ended with chi2 = " << res[9]
1410 << Form(", strip length %7.2f +/- %5.2f, position resolution "
1411 "%7.2f +/- %5.2f at y_cen = %7.2f +/- %5.2f",
1412 2. * res[1], 2. * err[1], res[2], err[2], res[3], err[3]);
1413}
1414
1416{
1417 if (fvLastHits.size() != static_cast<size_t>(fDigiBdfPar->GetNbSmTypes()))
1418 LOG(fatal) << Form("Inconsistent LH Smtype size %lu, %d ", fvLastHits.size(), fDigiBdfPar->GetNbSmTypes());
1419
1420 for (Int_t iSmType = 0; iSmType < fDigiBdfPar->GetNbSmTypes(); iSmType++) {
1421 if (fvLastHits[iSmType].size() != static_cast<size_t>(fDigiBdfPar->GetNbSm(iSmType)))
1422 LOG(fatal) << Form("Inconsistent LH Sm size %lu, %d T %d", fvLastHits[iSmType].size(),
1423 fDigiBdfPar->GetNbSm(iSmType), iSmType);
1424 for (Int_t iSm = 0; iSm < fDigiBdfPar->GetNbSm(iSmType); iSm++) {
1425 if (fvLastHits[iSmType][iSm].size() != static_cast<size_t>(fDigiBdfPar->GetNbRpc(iSmType)))
1426 LOG(fatal) << Form("Inconsistent LH Rpc size %lu, %d TS %d%d ", fvLastHits[iSmType][iSm].size(),
1427 fDigiBdfPar->GetNbRpc(iSmType), iSmType, iSm);
1428 for (Int_t iRpc = 0; iRpc < fDigiBdfPar->GetNbRpc(iSmType); iRpc++) {
1429 if (fvLastHits[iSmType][iSm][iRpc].size() != static_cast<size_t>(fDigiBdfPar->GetNbChan(iSmType, iRpc)))
1430 LOG(fatal) << Form("Inconsistent LH RpcChannel size %lu, %d TSR %d%d%d ",
1431 fvLastHits[iSmType][iSm][iRpc].size(), fDigiBdfPar->GetNbChan(iSmType, iRpc), iSmType, iSm,
1432 iRpc);
1433 for (Int_t iCh = 0; iCh < fDigiBdfPar->GetNbChan(iSmType, iRpc); iCh++)
1434 if (fvLastHits[iSmType][iSm][iRpc][iCh].size() > 0) {
1435 CbmTofDetectorInfo xDetInfo(ECbmModuleId::kTof, iSmType, iSm, iRpc, 0, iCh);
1436 Int_t iAddr = fTofId->SetDetectorInfo(xDetInfo);
1437 if (fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetAddress() != iAddr)
1438 LOG(fatal) << Form("Inconsistent address for Ev %8.0f in list of size %lu for "
1439 "TSRC %d%d%d%d: 0x%08x, time %f",
1440 fdEvent, fvLastHits[iSmType][iSm][iRpc][iCh].size(), iSmType, iSm, iRpc, iCh,
1441 fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetAddress(),
1442 fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetTime());
1443 }
1444 }
1445 }
1446 }
1447 LOG(debug) << Form("LH check passed for event %8.0f", fdEvent);
1448}
1449
1451{
1452 if (fvLastHits.size() != static_cast<size_t>(fDigiBdfPar->GetNbSmTypes()))
1453 LOG(fatal) << Form("Inconsistent LH Smtype size %lu, %d ", fvLastHits.size(), fDigiBdfPar->GetNbSmTypes());
1454 for (Int_t iSmType = 0; iSmType < fDigiBdfPar->GetNbSmTypes(); iSmType++) {
1455 if (fvLastHits[iSmType].size() != static_cast<size_t>(fDigiBdfPar->GetNbSm(iSmType)))
1456 LOG(fatal) << Form("Inconsistent LH Sm size %lu, %d T %d", fvLastHits[iSmType].size(),
1457 fDigiBdfPar->GetNbSm(iSmType), iSmType);
1458 for (Int_t iSm = 0; iSm < fDigiBdfPar->GetNbSm(iSmType); iSm++) {
1459 if (fvLastHits[iSmType][iSm].size() != static_cast<size_t>(fDigiBdfPar->GetNbRpc(iSmType)))
1460 LOG(fatal) << Form("Inconsistent LH Rpc size %lu, %d TS %d%d ", fvLastHits[iSmType][iSm].size(),
1461 fDigiBdfPar->GetNbRpc(iSmType), iSmType, iSm);
1462 for (Int_t iRpc = 0; iRpc < fDigiBdfPar->GetNbRpc(iSmType); iRpc++) {
1463 if (fvLastHits[iSmType][iSm][iRpc].size() != static_cast<size_t>(fDigiBdfPar->GetNbChan(iSmType, iRpc)))
1464 LOG(fatal) << Form("Inconsistent LH RpcChannel size %lu, %d TSR %d%d%d ",
1465 fvLastHits[iSmType][iSm][iRpc].size(), fDigiBdfPar->GetNbChan(iSmType, iRpc), iSmType, iSm,
1466 iRpc);
1467 for (Int_t iCh = 0; iCh < fDigiBdfPar->GetNbChan(iSmType, iRpc); iCh++)
1468 while (fvLastHits[iSmType][iSm][iRpc][iCh].size() > 0) {
1469 CbmTofDetectorInfo xDetInfo(ECbmModuleId::kTof, iSmType, iSm, iRpc, 0, iCh);
1470 Int_t iAddr = fTofId->SetDetectorInfo(xDetInfo);
1471 if (fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetAddress() != iAddr)
1472 LOG(fatal) << Form("Inconsistent address for Ev %8.0f in list of size %lu for "
1473 "TSRC %d%d%d%d: 0x%08x, time %f",
1474 fdEvent, fvLastHits[iSmType][iSm][iRpc][iCh].size(), iSmType, iSm, iRpc, iCh,
1475 fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetAddress(),
1476 fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetTime());
1477 fvLastHits[iSmType][iSm][iRpc][iCh].front()->Delete();
1478 fvLastHits[iSmType][iSm][iRpc][iCh].pop_front();
1479 }
1480 }
1481 }
1482 }
1483 LOG(info) << Form("LH cleaning done after %8.0f events", fdEvent);
1484}
1485
1486Bool_t CbmTofHitMaker::AddNextChan(Int_t iSmType, Int_t iSm, Int_t iRpc, Int_t iLastChan, Double_t dLastPosX,
1487 Double_t dLastPosY, Double_t dLastTime, Double_t dLastTotS)
1488{
1489 // Int_t iNbSm = fDigiBdfPar->GetNbSm( iSmType); (VF) not used
1490 Int_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
1491 Int_t iNbCh = fDigiBdfPar->GetNbChan(iSmType, iRpc);
1492 // Int_t iChType = fDigiBdfPar->GetChanType( iSmType, iRpc ); (VF) not used
1493 Int_t iDetId = CbmTofAddress::GetUniqueAddress(iSm, iRpc, 0, 0, iSmType);
1494 Int_t iDetIndx = fDetIdIndexMap[iDetId]; // Detector Index
1495
1496 Int_t iCh = iLastChan + 1;
1497 Int_t iChId = CbmTofAddress::GetUniqueAddress(iSm, iRpc, iCh, 0, iSmType);
1498
1499 while (fvDeadStrips[iDetIndx] & (1 << iCh)) {
1500 LOG(debug) << "Skip channel " << iCh << " of detector " << Form("0x%08x", iDetId);
1501 iCh++;
1502 iLastChan++;
1503 if (iCh >= iNbCh) return kFALSE;
1504 }
1505 LOG(debug1) << Form("Inspect channel TSRC %d%d%d%d at time %f, pos %f, size ", iSmType, iSm, iRpc, iCh, dLastTime,
1506 dLastPosY)
1507 << fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size();
1508 if (iCh == iNbCh) return kFALSE;
1509 if (0 == fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size()) return kFALSE;
1510 if (0 < fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size())
1511 fhNbDigiPerChan->Fill(fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size());
1512 if (1 < fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size()) {
1513 Bool_t AddedHit = kFALSE;
1514 for (size_t i1 = 0; i1 < fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size() - 1; i1++) {
1515 if (AddedHit) break;
1516 size_t i2 = i1 + 1;
1517 while (!AddedHit && i2 < fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size()) {
1518 LOG(debug1) << "check digi pair " << i1 << "," << i2 << " with size "
1519 << fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size();
1520
1521 if ((fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][i1])->GetSide()
1522 == (fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][i2])->GetSide()) {
1523 i2++;
1524 continue;
1525 } // endif same side
1526 // 2 Digis, both sides present
1527 CbmTofDigi* xDigiA = fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][i1];
1528 CbmTofDigi* xDigiB = fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][i2];
1529 Double_t dTime = 0.5 * (xDigiA->GetTime() + xDigiB->GetTime());
1530 if (TMath::Abs(dTime - dLastTime) < fdMaxTimeDist) {
1531 CbmTofDetectorInfo xDetInfo(ECbmModuleId::kTof, iSmType, iSm, iRpc, 0, iCh);
1532 iChId = fTofId->SetDetectorInfo(xDetInfo);
1533 fChannelInfo = fDigiPar->GetCell(iChId);
1534 gGeoManager->FindNode(fChannelInfo->GetX(), fChannelInfo->GetY(), fChannelInfo->GetZ());
1535
1536 Double_t dTimeDif = xDigiA->GetTime() - xDigiB->GetTime();
1537 Double_t dPosY = 0.;
1538 if (1 == xDigiA->GetSide())
1539 dPosY = fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc) * dTimeDif * 0.5;
1540 else
1541 dPosY = -fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc) * dTimeDif * 0.5;
1542
1543 if (TMath::Abs(dPosY - dLastPosY) < fdMaxSpaceDist) { // append digi pair to current cluster
1544
1545 Double_t dNClHits = (Double_t)(vDigiIndRef.size() / 2);
1546 Double_t dPosX = ((Double_t)(-iNbCh / 2 + iCh) + 0.5) * fChannelInfo->GetSizex();
1547 Double_t dTotS = xDigiA->GetTot() + xDigiB->GetTot();
1548 Double_t dNewTotS = (dLastTotS + dTotS);
1549 dLastPosX = (dLastPosX * dLastTotS + dPosX * dTotS) / dNewTotS;
1550 dLastPosY = (dLastPosY * dLastTotS + dPosY * dTotS) / dNewTotS;
1551 dLastTime = (dLastTime * dLastTotS + dTime * dTotS) / dNewTotS;
1552 dLastTotS = dNewTotS;
1553 // attach selected digis from pool
1554 Int_t Ind1 = fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][i1];
1555 Int_t Ind2 = fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][i2];
1556 vDigiIndRef.push_back(Ind1);
1557 vDigiIndRef.push_back(Ind2);
1558 // remove selected digis from pool
1559 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin()
1560 + i1);
1561 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1562 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin() + i1);
1563
1564 std::vector<int>::iterator it;
1565 it = find(fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin(),
1566 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].end(), Ind2);
1567 if (it != fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].end()) {
1568 auto ipos = it - fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin();
1569 LOG(debug1) << "Found i2 " << i2 << " with Ind2 " << Ind2 << " at position " << ipos;
1570 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin()
1571 + ipos);
1572 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1573 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin() + ipos);
1574 }
1575 else {
1576 LOG(fatal) << " Did not find i2 " << i2 << " with Ind2 " << Ind2;
1577 }
1578
1579 //if(iCh == iNbCh-1) break; //Last strip reached
1580 if (iCh != (iNbCh - 1)
1581 && AddNextChan(iSmType, iSm, iRpc, iCh, dLastPosX, dLastPosY, dLastTime, dLastTotS)) {
1582 LOG(debug1) << "Added Strip " << iCh << " to cluster of size " << dNClHits;
1583 return kTRUE; // signal hit was already added
1584 }
1585 AddedHit = kTRUE;
1586 } //TMath::Abs(dPosY - dLastPosY) < fdMaxSpaceDist
1587 } //TMath::Abs(dTime-dLastTime)<fdMaxTimeDist)
1588 i2++;
1589 } // while(i2 < fStorDigi[iSmType][iSm*iNbRpc+iRpc][iCh].size()-1 )
1590 } // end for i1
1591 } // end if size
1592 Double_t hitpos_local[3] = {3 * 0.};
1593 hitpos_local[0] = dLastPosX;
1594 hitpos_local[1] = dLastPosY;
1595 hitpos_local[2] = 0.;
1596 /*
1597 if( 5 == iSmType || 8 == iSmType) { // for PAD counters
1598 hitpos_local[0] = (gRandom->Rndm()-0.5)*fChannelInfo->GetSizex()*0.5;
1599 hitpos_local[1] = (gRandom->Rndm()-0.5)*fChannelInfo->GetSizey()*0.5;
1600 }
1601 */
1602 Double_t hitpos[3] = {3 * 0.};
1603 if (5 != iSmType) { // Diamond beam counter always at (0,0,0)
1604 /*TGeoNode* cNode = */ gGeoManager->GetCurrentNode();
1605 /*TGeoHMatrix* cMatrix = */ gGeoManager->GetCurrentMatrix();
1606 gGeoManager->LocalToMaster(hitpos_local, hitpos);
1607 }
1608 TVector3 hitPos(hitpos[0], hitpos[1], hitpos[2]);
1609 TVector3 hitPosErr(0.5, 0.5, 0.5); // FIXME including positioning uncertainty
1610 Int_t iChm = floor(dLastPosX / fChannelInfo->GetSizex()) + iNbCh / 2;
1611 if (iChm < 0) iChm = 0;
1612 if (iChm > iNbCh - 1) iChm = iNbCh - 1;
1613 iDetId = CbmTofAddress::GetUniqueAddress(iSm, iRpc, iChm, 0, iSmType);
1614
1615 Int_t iNbChanInHit = vDigiIndRef.size() / 2;
1616
1617 TString cstr = "Save A-Hit ";
1618 cstr += Form(" %3d %3d 0x%08x %3d 0x%08x %8.2f %6.2f", // %3d %3d
1619 fiNbHits, iNbChanInHit, iDetId, iLastChan,
1620 0, //vPtsRef.size(),vPtsRef[0])
1621 dLastTime, dLastPosY);
1622 cstr += Form(", DigiSize: %lu ", vDigiIndRef.size());
1623 cstr += ", DigiInds: ";
1624
1625 fviClusterMul[iSmType][iSm][iRpc]++;
1626
1627 for (UInt_t i = 0; i < vDigiIndRef.size(); i++) {
1628 cstr += Form(" %d (M,%d)", vDigiIndRef.at(i), fviClusterMul[iSmType][iSm][iRpc]);
1629 }
1630 LOG(debug) << cstr;
1631
1632 CbmTofHit* pHit = new CbmTofHit(iDetId, hitPos,
1633 hitPosErr, //local detector coordinates
1634 fiNbHits, // this number is used as reference!!
1635 dLastTime,
1636 vDigiIndRef.size(), // number of linked digis = 2*CluSize
1637 //vPtsRef.size(), // flag = number of TofPoints generating the cluster
1638 Int_t(dLastTotS * 10.)); //channel -> Tot
1639 // output hit
1640 new ((*fTofHitsColl)[fiNbHits]) CbmTofHit(*pHit);
1641 if (fdMemoryTime > 0.) { // memorize hit
1642 LH_store(iSmType, iSm, iRpc, iChm, pHit);
1643 }
1644 else {
1645 pHit->Delete();
1646 }
1647 CbmMatch* digiMatch = new ((*fTofDigiMatchColl)[fiNbHits]) CbmMatch();
1648 for (size_t i = 0; i < vDigiIndRef.size(); i++) {
1649 Double_t dTot = (fTofCalDigiVec->at(vDigiIndRef.at(i))).GetTot();
1650 digiMatch->AddLink(CbmLink(dTot, vDigiIndRef.at(i), fiOutputTreeEntry, fiFileIndex));
1651 }
1652 fiNbHits++;
1653 vDigiIndRef.clear();
1654
1655 return kTRUE;
1656}
1657
1658void CbmTofHitMaker::LH_store(Int_t iSmType, Int_t iSm, Int_t iRpc, Int_t iChm, CbmTofHit* pHit)
1659{
1660
1661 if (fvLastHits[iSmType][iSm][iRpc][iChm].size() == 0)
1662 fvLastHits[iSmType][iSm][iRpc][iChm].push_back(pHit);
1663 else {
1664 Double_t dLastTime = pHit->GetTime();
1665 if (dLastTime >= fvLastHits[iSmType][iSm][iRpc][iChm].back()->GetTime()) {
1666 fvLastHits[iSmType][iSm][iRpc][iChm].push_back(pHit);
1667 LOG(debug) << Form(" Store LH from Ev %8.0f for TSRC %d%d%d%d, size %lu, addr 0x%08x, "
1668 "time %f, dt %f",
1669 fdEvent, iSmType, iSm, iRpc, iChm, fvLastHits[iSmType][iSm][iRpc][iChm].size(),
1670 pHit->GetAddress(), dLastTime,
1671 dLastTime - fvLastHits[iSmType][iSm][iRpc][iChm].front()->GetTime());
1672 }
1673 else {
1674 if (dLastTime
1675 >= fvLastHits[iSmType][iSm][iRpc][iChm].front()->GetTime()) { // hit has to be inserted in the proper place
1676 std::list<CbmTofHit*>::iterator it;
1677 for (it = fvLastHits[iSmType][iSm][iRpc][iChm].begin(); it != fvLastHits[iSmType][iSm][iRpc][iChm].end(); ++it)
1678 if ((*it)->GetTime() > dLastTime) break;
1679 fvLastHits[iSmType][iSm][iRpc][iChm].insert(--it, pHit);
1680 Double_t deltaTime = dLastTime - (*it)->GetTime();
1681 LOG(debug) << Form("Hit inserted into LH from Ev %8.0f for TSRC "
1682 "%d%d%d%d, size %lu, addr 0x%08x, delta time %f ",
1683 fdEvent, iSmType, iSm, iRpc, iChm, fvLastHits[iSmType][iSm][iRpc][iChm].size(),
1684 pHit->GetAddress(), deltaTime);
1685 }
1686 else { // this hit is first
1687 Double_t deltaTime = dLastTime - fvLastHits[iSmType][iSm][iRpc][iChm].front()->GetTime();
1688 LOG(debug) << Form("first LH from Ev %8.0f for TSRC %d%d%d%d, size "
1689 "%lu, addr 0x%08x, delta time %f ",
1690 fdEvent, iSmType, iSm, iRpc, iChm, fvLastHits[iSmType][iSm][iRpc][iChm].size(),
1691 pHit->GetAddress(), deltaTime);
1692 if (deltaTime == 0.) {
1693 // remove hit, otherwise double entry?
1694 pHit->Delete();
1695 }
1696 else {
1697 fvLastHits[iSmType][iSm][iRpc][iChm].push_front(pHit);
1698 }
1699 }
1700 }
1701 }
1702}
1703
1705{
1706 // Then build clusters inside each RPC module
1707 // Assume only 0 or 1 Digi per channel/side in each event
1708 // Use simplest method possible, scan direction independent:
1709 // a) Loop over channels in the RPC starting from 0
1710 // * If strips
1711 // i) Loop over Digis to check if both ends of the channel have a Digi
1712 // ii) Reconstruct a mean channel time and a mean position
1713 // + If a Hit is currently filled & the mean position (space, time) is less than XXX from last channel position
1714 // iii) Add the mean channel time and the mean position to the ones of the hit
1715 // + else
1716 // iii) Use nb of strips in cluster to cal. the hit mean time and pos (charge/tot weighting)
1717 // iv) Save the hit
1718 // v) Start a new hit with current channel
1719 // * else (pads)
1720 // i) Loop over Digis to find if this channel fired
1721 // ii) FIXME: either scan all other channels to check for matching Digis or have more than 1 hit open
1722 Int_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes();
1723 // Hit variables
1724 Double_t dWeightedTime = 0.0;
1725 Double_t dWeightedPosX = 0.0;
1726 Double_t dWeightedPosY = 0.0;
1727 Double_t dWeightedPosZ = 0.0;
1728 Double_t dWeightsSum = 0.0;
1729 //vPtsRef.clear();
1730 vDigiIndRef.clear();
1731 // CbmTofCell *fTrafoCell=NULL; (VF) not used
1732 // Int_t iTrafoCell=-1; (VF) not used
1733 Int_t iNbChanInHit = 0;
1734 // Last Channel Temp variables
1735 Int_t iLastChan = -1;
1736 Double_t dLastPosX = 0.0; // -> Comment to remove warning because set but never used
1737 Double_t dLastPosY = 0.0;
1738 Double_t dLastTime = 0.0;
1739 // Channel Temp variables
1740 Double_t dPosX = 0.0;
1741 Double_t dPosY = 0.0;
1742 Double_t dPosZ = 0.0;
1743 Double_t dTime = 0.0;
1744 Double_t dTimeDif = 0.0;
1745 Double_t dTotS = 0.0;
1746 fiNbSameSide = 0;
1747 if (kTRUE == fDigiBdfPar->UseExpandedDigi()) {
1748 for (Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
1749 Int_t iNbSm = fDigiBdfPar->GetNbSm(iSmType);
1750 Int_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
1751 for (Int_t iSm = 0; iSm < iNbSm; iSm++)
1752 for (Int_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
1753 Int_t iNbCh = fDigiBdfPar->GetNbChan(iSmType, iRpc);
1754 Int_t iChType = fDigiBdfPar->GetChanType(iSmType, iRpc);
1755 LOG(debug2) << "RPC - Loop " << Form(" %3d %3d %3d %3d ", iSmType, iSm, iRpc, iChType);
1756 fviClusterMul[iSmType][iSm][iRpc] = 0;
1757 Int_t iChId = 0;
1758 Int_t iDetId = CbmTofAddress::GetUniqueAddress(iSm, iRpc, 0, 0, iSmType);
1759 ;
1760 Int_t iDetIndx = fDetIdIndexMap[iDetId]; // Detector Index
1761 if (0 == iChType) {
1762 // Don't spread clusters over RPCs!!!
1763 dWeightedTime = 0.0;
1764 dWeightedPosX = 0.0;
1765 dWeightedPosY = 0.0;
1766 dWeightedPosZ = 0.0;
1767 dWeightsSum = 0.0;
1768 iNbChanInHit = 0;
1769 //vPtsRef.clear();
1770 // For safety reinitialize everything
1771 iLastChan = -1;
1772 // dLastPosX = 0.0; // -> Comment to remove warning because set but never used
1773 dLastPosY = 0.0;
1774 dLastTime = 0.0;
1775 LOG(debug2) << "ChanOrient "
1776 << Form(" %3d %3d %3d %3d %3d ", iSmType, iSm, iRpc, fDigiBdfPar->GetChanOrient(iSmType, iRpc),
1777 iNbCh);
1778
1779 if (1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc)) {
1780 // Horizontal strips => X comes from left right time difference
1781 } // if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
1782 else { // Vertical strips => Y comes from bottom top time difference
1783 for (Int_t iCh = 0; iCh < iNbCh; iCh++) {
1784 LOG(debug3) << "VDigisize "
1785 << Form(" T %3d Sm %3d R %3d Ch %3d Size %3lu ", iSmType, iSm, iRpc, iCh,
1786 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size());
1787 if (0 == fStorDigi[iSmType][iSm * iNbRpc + iRpc].size()) continue;
1788 if (fvDeadStrips[iDetIndx] & (1 << iCh)) continue; // skip over dead channels
1789 if (0 < fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size())
1790 fhNbDigiPerChan->Fill(fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size());
1791
1792 while (1 < fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size()) {
1793
1794 while ((fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][0])->GetSide()
1795 == (fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][1])->GetSide()) {
1796 // Not one Digi of each end!
1797 fiNbSameSide++;
1798 if (fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size() > 2) {
1799 LOG(debug) << "SameSide Digis! on TSRC " << iSmType << iSm << iRpc << iCh << ", Times: "
1800 << Form("%f", (fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][0])->GetTime()) << ", "
1801 << Form("%f", (fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][1])->GetTime())
1802 << ", DeltaT "
1803 << (fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][1])->GetTime()
1804 - (fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][0])->GetTime()
1805 << ", array size: " << fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size();
1806 if (fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][2]->GetSide()
1807 == fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][0]->GetSide()) {
1808 LOG(debug) << "3 consecutive SameSide Digis! on TSRC " << iSmType << iSm << iRpc << iCh
1809 << ", Times: " << (fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][0])->GetTime()
1810 << ", " << (fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][1])->GetTime()
1811 << ", DeltaT "
1812 << (fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][1])->GetTime()
1813 - (fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][0])->GetTime()
1814 << ", array size: " << fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size();
1815 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1816 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
1817 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1818 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
1819 }
1820 else {
1821 if (fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][2]->GetTime()
1822 - fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][0]->GetTime()
1823 > fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][2]->GetTime()
1824 - fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][1]->GetTime()) {
1825 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1826 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
1827 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1828 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
1829 }
1830 else {
1831 LOG(debug) << Form("Ev %8.0f, digis not properly time ordered, TSRCS "
1832 "%d%d%d%d%d ",
1833 fdEvent, iSmType, iSm, iRpc, iCh,
1834 (Int_t) fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][0]->GetSide());
1835 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1836 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin() + 1);
1837 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1838 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin() + 1);
1839 }
1840 }
1841 }
1842 else {
1843 LOG(debug2) << "SameSide Erase fStor entries(d) " << iSmType << ", SR " << iSm * iNbRpc + iRpc
1844 << ", Ch" << iCh;
1845 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1846 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
1847 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1848 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
1849 }
1850 if (2 > fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size()) break;
1851 continue;
1852 } // same condition side end
1853
1854 LOG(debug2) << "digis processing for "
1855 << Form(" SmT %3d Sm %3d Rpc %3d Ch %3d # %3lu ", iSmType, iSm, iRpc, iCh,
1856 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size());
1857 if (2 > fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size()) {
1858 LOG(debug) << Form("Leaving digi processing for TSRC %d%d%d%d, size %3lu", iSmType, iSm, iRpc, iCh,
1859 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size());
1860 break;
1861 }
1862 /* Int_t iLastChId = iChId; // Save Last hit channel*/
1863
1864 // 2 Digis = both sides present
1865 CbmTofDetectorInfo xDetInfo(ECbmModuleId::kTof, iSmType, iSm, iRpc, 0, iCh);
1866 iChId = fTofId->SetDetectorInfo(xDetInfo);
1867 Int_t iUCellId = CbmTofAddress::GetUniqueAddress(iSm, iRpc, iCh, 0, iSmType);
1868 LOG(debug1) << Form(" TSRC %d%d%d%d size %3lu ", iSmType, iSm, iRpc, iCh,
1869 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size())
1870 << Form(" ChId: 0x%08x 0x%08x ", iChId, iUCellId);
1871 fChannelInfo = fDigiPar->GetCell(iChId);
1872
1873 if (NULL == fChannelInfo) {
1874 LOG(error) << "CbmTofHitMaker::BuildClusters: no "
1875 "geometry info! "
1876 << Form(" %3d %3d %3d %3d 0x%08x 0x%08x ", iSmType, iSm, iRpc, iCh, iChId, iUCellId);
1877 break;
1878 }
1879
1880 TGeoNode* fNode = // prepare local->global trafo
1881 gGeoManager->FindNode(fChannelInfo->GetX(), fChannelInfo->GetY(), fChannelInfo->GetZ());
1882 LOG(debug2) << Form(" Node at (%6.1f,%6.1f,%6.1f) : %p", fChannelInfo->GetX(), fChannelInfo->GetY(),
1883 fChannelInfo->GetZ(), fNode);
1884 // fNode->Print();
1885
1886 CbmTofDigi* xDigiA = fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][0];
1887 CbmTofDigi* xDigiB = fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][1];
1888
1889 LOG(debug2) << " " << xDigiA->ToString();
1890 LOG(debug2) << " " << xDigiB->ToString();
1891
1892 dTimeDif = (xDigiA->GetTime() - xDigiB->GetTime());
1893 if (5 == iSmType && dTimeDif != 0.) {
1894 // FIXME -> Overflow treatment in calib/tdc/TMbsCalibTdcTof.cxx
1895 LOG(debug) << "CbmTofHitMaker::BuildClusters: "
1896 "Diamond hit in "
1897 << iSm << " with inconsistent digits " << xDigiA->GetTime() << ", " << xDigiB->GetTime()
1898 << " -> " << dTimeDif;
1899 LOG(debug) << " " << xDigiA->ToString();
1900 LOG(debug) << " " << xDigiB->ToString();
1901 }
1902 if (1 == xDigiA->GetSide())
1903 // 0 is the top side, 1 is the bottom side
1904 dPosY = fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc) * dTimeDif * 0.5;
1905 else
1906 // 0 is the bottom side, 1 is the top side
1907 dPosY = -fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc) * dTimeDif * 0.5;
1908
1909 while (TMath::Abs(dPosY) > fChannelInfo->GetSizey() * fPosYMaxScal
1910 && fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size() > 2) {
1911 LOG(debug) << "Hit candidate outside correlation window, check for "
1912 "better possible digis, "
1913 << " mul " << fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size();
1914
1915 CbmTofDigi* xDigiC = fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][2];
1916 Double_t dPosYN = 0.;
1917 Double_t dTimeDifN = 0;
1918 if (xDigiC->GetSide() == xDigiA->GetSide())
1919 dTimeDifN = xDigiC->GetTime() - xDigiB->GetTime();
1920 else
1921 dTimeDifN = xDigiA->GetTime() - xDigiC->GetTime();
1922
1923 if (1 == xDigiA->GetSide())
1924 dPosYN = fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc) * dTimeDifN * 0.5;
1925 else
1926 dPosYN = -fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc) * dTimeDifN * 0.5;
1927
1928 if (TMath::Abs(dPosYN) < TMath::Abs(dPosY)) {
1929 LOG(debug) << "Replace digi on side " << xDigiC->GetSide() << ", yPosNext " << dPosYN
1930 << " old: " << dPosY;
1931 dTimeDif = dTimeDifN;
1932 dPosY = dPosYN;
1933 if (xDigiC->GetSide() == xDigiA->GetSide()) {
1934 xDigiA = xDigiC;
1935 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1936 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
1937 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1938 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
1939 }
1940 else {
1941 xDigiB = xDigiC;
1942 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1943 ++(fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin() + 1));
1944 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1945 ++(fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin() + 1));
1946 }
1947 }
1948 else
1949 break;
1950 } //while loop end
1951
1952 if (xDigiA->GetSide() == xDigiB->GetSide()) {
1953 LOG(fatal) << "Wrong combinations of digis " << fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][0]
1954 << "," << fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][1];
1955 }
1956
1957 if (TMath::Abs(dPosY) > fChannelInfo->GetSizey() * fPosYMaxScal) { // remove both digis
1958 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1959 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
1960 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
1961 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
1962 continue;
1963 }
1964 // The "Strip" time is the mean time between each end
1965 dTime = 0.5 * (xDigiA->GetTime() + xDigiB->GetTime());
1966
1967 // Weight is the total charge => sum of both ends ToT
1968 dTotS = xDigiA->GetTot() + xDigiB->GetTot();
1969
1970 // use local coordinates, (0,0,0) is in the center of counter ?
1971 dPosX = ((Double_t)(-iNbCh / 2 + iCh) + 0.5) * fChannelInfo->GetSizex();
1972 dPosZ = 0.;
1973
1974 LOG(debug1) << "NbChanInHit "
1975 << Form(" %3d %3d %3d %3d %3d 0x%p %1.0f Time %f PosX %f "
1976 "PosY %f Svel %f ",
1977 iNbChanInHit, iSmType, iRpc, iCh, iLastChan, xDigiA, xDigiA->GetSide(), dTime,
1978 dPosX, dPosY, fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc))
1979 // << Form( ", Offs %f, %f ",fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][0],
1980 // fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][1])
1981 ;
1982
1983 // Now check if a hit/cluster is already started
1984 if (0 < iNbChanInHit) {
1985 if (iLastChan == iCh - 1) {
1986 fhDigTimeDifClust->Fill(dTime - dLastTime);
1987 fhDigSpacDifClust->Fill(dPosY - dLastPosY);
1988 fhDigDistClust->Fill(dPosY - dLastPosY, dTime - dLastTime);
1989 }
1990 // if( iLastChan == iCh - 1 )
1991 // a cluster is already started => check distance in space/time
1992 // For simplicity, just check along strip direction for now
1993 // and break cluster when a not fired strip is found
1994 if (TMath::Abs(dTime - dLastTime) < fdMaxTimeDist && iLastChan == iCh - 1
1995 && TMath::Abs(dPosY - dLastPosY) < fdMaxSpaceDist) {
1996 // Add to cluster/hit
1997 dWeightedTime += dTime * dTotS;
1998 dWeightedPosX += dPosX * dTotS;
1999 dWeightedPosY += dPosY * dTotS;
2000 dWeightedPosZ += dPosZ * dTotS;
2001 dWeightsSum += dTotS;
2002 iNbChanInHit += 1;
2003
2004 vDigiIndRef.push_back((Int_t)(fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][0]));
2005 vDigiIndRef.push_back((Int_t)(fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][1]));
2006
2007 LOG(debug1) << " Add Digi and erase fStor entries(a): NbChanInHit " << iNbChanInHit << ", "
2008 << iSmType << ", SR " << iSm * iNbRpc + iRpc << ", Ch" << iCh;
2009
2010 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
2011 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
2012 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
2013 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
2014 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
2015 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
2016 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
2017 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
2018
2019 } // if current Digis compatible with last fired chan
2020 else {
2021 // Save Hit
2022 dWeightedTime /= dWeightsSum;
2023 dWeightedPosX /= dWeightsSum;
2024 dWeightedPosY /= dWeightsSum;
2025 dWeightedPosZ /= dWeightsSum;
2026 // TVector3 hitPosLocal(dWeightedPosX, dWeightedPosY, dWeightedPosZ);
2027 //TVector3 hitPos;
2028 Double_t hitpos_local[3];
2029 hitpos_local[0] = dWeightedPosX;
2030 hitpos_local[1] = dWeightedPosY;
2031 hitpos_local[2] = dWeightedPosZ;
2032 /*
2033 if( 5 == iSmType || 8 == iSmType) { // for PAD counters
2034 hitpos_local[0] = (gRandom->Rndm()-0.5)*fChannelInfo->GetSizex();
2035 hitpos_local[1] = (gRandom->Rndm()-0.5)*fChannelInfo->GetSizey();
2036 }
2037 */
2038 Double_t hitpos[3] = {3 * 0.};
2039 if (5 != iSmType) {
2040 /*TGeoNode* cNode =*/gGeoManager->GetCurrentNode();
2041 /*TGeoHMatrix* cMatrix =*/gGeoManager->GetCurrentMatrix();
2042 //cNode->Print();
2043 //cMatrix->Print();
2044
2045 gGeoManager->LocalToMaster(hitpos_local, hitpos);
2046 }
2047 LOG(debug1) << Form(" LocalToMaster: (%6.1f,%6.1f,%6.1f) "
2048 "->(%6.1f,%6.1f,%6.1f)",
2049 hitpos_local[0], hitpos_local[1], hitpos_local[2], hitpos[0], hitpos[1],
2050 hitpos[2]);
2051
2052 TVector3 hitPos(hitpos[0], hitpos[1], hitpos[2]);
2053
2054 // Simple errors, not properly done at all for now
2055 // Right way of doing it should take into account the weight distribution
2056 // and real system time resolution
2057 TVector3 hitPosErr(0.5, 0.5, 0.5); // including positioning uncertainty
2058 /*
2059 TVector3 hitPosErr( fChannelInfo->GetSizex()/TMath::Sqrt(12.0), // Single strips approximation
2060 0.5, // Use generic value
2061 1.);
2062
2063 */
2064 //fDigiBdfPar->GetFeeTimeRes() * fDigiBdfPar->GetSigVel(iSmType,iRpc), // Use the electronics resolution
2065 //fDigiBdfPar->GetNbGaps( iSmType, iRpc)*
2066 //fDigiBdfPar->GetGapSize( iSmType, iRpc)/ //10.0 / // Change gap size in cm
2067 //TMath::Sqrt(12.0) ); // Use full RPC thickness as "Channel" Z size
2068 // Int_t iDetId = vPtsRef[0]->GetDetectorID();// detID = pt->GetDetectorID() <= from TofPoint
2069 // calc mean ch from dPosX=((Double_t)(-iNbCh/2 + iCh)+0.5)*fChannelInfo->GetSizex();
2070 Int_t iChm = floor(dWeightedPosX / fChannelInfo->GetSizex()) + iNbCh / 2;
2071 if (iChm < 0) iChm = 0;
2072 if (iChm > iNbCh - 1) iChm = iNbCh - 1;
2073 iDetId = CbmTofAddress::GetUniqueAddress(iSm, iRpc, iChm, 0, iSmType);
2074 Int_t iRefId = 0; // Index of the correspondng TofPoint
2075 // if(NULL != fTofPointsColl) {
2076 //iRefId = fTofPointsColl->IndexOf( vPtsRef[0] );
2077 //}
2078 TString sRef = "";
2079 for (UInt_t i = 0; i < vDigiIndRef.size(); i++) {
2080 sRef += Form(" %d, (M%d)", vDigiIndRef.at(i), fviClusterMul[iSmType][iSm][iRpc]);
2081 }
2082 LOG(debug) << "Save Hit "
2083 << Form(" %3d %3d 0x%08x %3d %3d %3d %f %f", fiNbHits, iNbChanInHit, iDetId, iChm,
2084 iLastChan, iRefId, dWeightedTime, dWeightedPosY)
2085 << ", DigiSize: " << vDigiIndRef.size() << ", DigiInds: " << sRef;
2086
2087 fviClusterMul[iSmType][iSm][iRpc]++;
2088
2089 if (vDigiIndRef.size() < 2) {
2090 LOG(warning) << "Digi refs for Hit " << fiNbHits << ": vDigiIndRef.size()";
2091 }
2092 if (fiNbHits > 0) {
2093 CbmTofHit* pHitL = (CbmTofHit*) fTofHitsColl->At(fiNbHits - 1);
2094 if (iDetId == pHitL->GetAddress() && dWeightedTime == pHitL->GetTime()) {
2095 LOG(debug) << "Store Hit twice? "
2096 << " fiNbHits " << fiNbHits << ", "
2097 << Form("0x%08x, MatchCollSize %d, IndRefSize %lu ", iDetId,
2098 fTofDigiMatchColl->GetEntriesFast(), vDigiIndRef.size());
2099
2100 for (UInt_t i = 0; i < vDigiIndRef.size(); i++) {
2101 if (vDigiIndRef.at(i) < (Int_t) fTofCalDigiVec->size()) {
2102 CbmTofDigi* pDigiC = &(fTofCalDigiVec->at(vDigiIndRef.at(i)));
2103 LOG(debug) << " Digi " << i << " " << pDigiC->ToString();
2104 }
2105 else {
2106 LOG(fatal) << "Insufficient CalDigiVec size for i " << i << ", Ind " << vDigiIndRef.at(i);
2107 }
2108 }
2109
2110 if (NULL == fTofDigiMatchColl) assert("No DigiMatchColl");
2111 CbmMatch* digiMatchL = NULL;
2112 if (fTofDigiMatchColl->GetEntriesFast() >= fiNbHits - 1) {
2113 digiMatchL = (CbmMatch*) fTofDigiMatchColl->At(fiNbHits - 1);
2114 }
2115 else {
2116 LOG(fatal) << "DigiMatchColl has insufficient size " << fTofDigiMatchColl->GetEntriesFast();
2117 }
2118
2119 if (NULL != digiMatchL)
2120 for (Int_t i = 0; i < digiMatchL->GetNofLinks(); i++) {
2121 CbmLink L0 = digiMatchL->GetLink(i);
2122 LOG(debug) << "report link " << i << "(" << digiMatchL->GetNofLinks() << "), ind "
2123 << L0.GetIndex();
2124 Int_t iDigIndL = L0.GetIndex();
2125 if (iDigIndL >= (Int_t) vDigiIndRef.size()) {
2126 if (iDetId != fiBeamRefAddr) {
2127 LOG(warn) << Form("Invalid DigiRefInd for det 0x%08x", iDetId);
2128 continue;
2129 }
2130 }
2131 if (vDigiIndRef.at(iDigIndL) >= (Int_t) fTofCalDigiVec->size()) {
2132 LOG(warn) << "Invalid CalDigiInd";
2133 continue;
2134 }
2135 CbmTofDigi* pDigiC = &(fTofCalDigiVec->at(vDigiIndRef.at(iDigIndL)));
2136 LOG(debug) << " DigiL " << pDigiC->ToString();
2137 }
2138 else {
2139 LOG(warn) << "Invalid digMatch Link at Index " << fiNbHits - 1;
2140 }
2141 }
2142 LOG(debug) << "Current HitsColl length " << fTofHitsColl->GetEntriesFast();
2143 }
2144 CbmTofHit* pHit =
2145 new CbmTofHit(iDetId, hitPos,
2146 hitPosErr, //local detector coordinates
2147 fiNbHits, // this number is used as reference!!
2148 dWeightedTime,
2149 vDigiIndRef.size(), // number of linked digis = 2*CluSize
2150 //vPtsRef.size(), // flag = number of TofPoints generating the cluster
2151 Int_t(dWeightsSum * 10.)); //channel -> Tot
2152 //0) ; //channel
2153 // output hit
2154 new ((*fTofHitsColl)[fiNbHits]) CbmTofHit(*pHit);
2155 // memorize hit
2156 if (fdMemoryTime > 0.) {
2157 LH_store(iSmType, iSm, iRpc, iChm, pHit);
2158 }
2159 else {
2160 pHit->Delete();
2161 }
2162 /*
2163 new((*fTofDigiMatchColl)[fiNbHits]) CbmMatch();
2164 CbmMatch* digiMatch = (CbmMatch *)fTofDigiMatchColl->At(fiNbHits);
2165 */
2166 CbmMatch* digiMatch = new ((*fTofDigiMatchColl)[fiNbHits]) CbmMatch();
2167 for (size_t i = 0; i < vDigiIndRef.size(); i++) {
2168 Double_t dTot = (fTofCalDigiVec->at(vDigiIndRef.at(i))).GetTot();
2169 digiMatch->AddLink(CbmLink(dTot, vDigiIndRef.at(i), fiOutputTreeEntry, fiFileIndex));
2170 }
2171
2172 fiNbHits++;
2173 // For Histogramming
2174 fviClusterSize[iSmType][iRpc].push_back(iNbChanInHit);
2175 //fviTrkMul[iSmType][iRpc].push_back( vPtsRef.size() );
2176 fvdX[iSmType][iRpc].push_back(dWeightedPosX);
2177 fvdY[iSmType][iRpc].push_back(dWeightedPosY);
2178 /* no TofPoint available for data!
2179 fvdDifX[iSmType][iRpc].push_back( vPtsRef[0]->GetX() - dWeightedPosX);
2180 fvdDifY[iSmType][iRpc].push_back( vPtsRef[0]->GetY() - dWeightedPosY);
2181 fvdDifCh[iSmType][iRpc].push_back( fGeoHandler->GetCell( vPtsRef[0]->GetDetectorID() ) -1 -iLastChan );
2182 */
2183 //vPtsRef.clear();
2184 vDigiIndRef.clear();
2185
2186 // Start a new hit
2187 dWeightedTime = dTime * dTotS;
2188 dWeightedPosX = dPosX * dTotS;
2189 dWeightedPosY = dPosY * dTotS;
2190 dWeightedPosZ = dPosZ * dTotS;
2191 dWeightsSum = dTotS;
2192 iNbChanInHit = 1;
2193 // Save pointer on CbmTofPoint
2194 // vPtsRef.push_back( (CbmTofPoint*)(xDigiA->GetLinks()) );
2195 // Save next digi address
2196 LOG(debug2) << " Next fStor Digi " << iSmType << ", SR " << iSm * iNbRpc + iRpc << ", Ch" << iCh
2197 << ", Dig0 " << (fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][0]) << ", Dig1 "
2198 << (fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][1]);
2199
2200 vDigiIndRef.push_back((Int_t)(fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][0]));
2201 vDigiIndRef.push_back((Int_t)(fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][1]));
2202 LOG(debug2) << " Erase fStor entries(b) " << iSmType << ", SR " << iSm * iNbRpc + iRpc << ", Ch"
2203 << iCh;
2204 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
2205 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
2206 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
2207 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
2208 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
2209 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
2210 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
2211 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
2212
2213 if (kTRUE == fDigiBdfPar->ClustUseTrackId()) {
2214 // if( ((CbmTofPoint*)(xDigiA->GetLinks()))->GetTrackID() !=
2215 // ((CbmTofPoint*)(xDigiB->GetLinks()))->GetTrackID() )
2216 // if other side come from a different Track,
2217 // also save the pointer on CbmTofPoint
2218 // vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) );
2219 } // if( kTRUE == fDigiBdfPar->ClustUseTrackId() )
2220 //else if( xDigiA->GetLinks() != xDigiB->GetLinks() )
2221 // if other side come from a different TofPoint,
2222 // also save the pointer on CbmTofPoint
2223 // vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) );
2224 } // else of if current Digis compatible with last fired chan
2225 } // if( 0 < iNbChanInHit)
2226 else {
2227 LOG(debug) << Form("1.Hit on TSRC %d%d%d%d, time: %f, PosY %f, Tdif %f ", iSmType, iSm, iRpc, iCh,
2228 dTime, dPosY, dTimeDif);
2229
2230 // first fired strip in this RPC
2231 dWeightedTime = dTime * dTotS;
2232 dWeightedPosX = dPosX * dTotS;
2233 dWeightedPosY = dPosY * dTotS;
2234 dWeightedPosZ = dPosZ * dTotS;
2235 dWeightsSum = dTotS;
2236 iNbChanInHit = 1;
2237 // Save pointer on CbmTofPoint
2238 //if(NULL != fTofPointsColl)
2239 // vPtsRef.push_back( (CbmTofPoint*)(xDigiA->GetLinks()) );
2240 vDigiIndRef.push_back((Int_t)(fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][0]));
2241 vDigiIndRef.push_back((Int_t)(fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][1]));
2242
2243 LOG(debug2) << " Erase fStor entries(c) " << iSmType << ", SR " << iSm * iNbRpc + iRpc << ", Ch"
2244 << iCh;
2245 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
2246 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
2247 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
2248 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
2249 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
2250 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
2251 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].erase(
2252 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
2253
2254 if (kTRUE == fDigiBdfPar->ClustUseTrackId()) {
2255 // if( ((CbmTofPoint*)(xDigiA->GetLinks()))->GetTrackID() !=
2256 // ((CbmTofPoint*)(xDigiB->GetLinks()))->GetTrackID() )
2257 // if other side come from a different Track,
2258 // also save the pointer on CbmTofPoint
2259 // vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) );
2260 } // if( kTRUE == fDigiBdfPar->ClustUseTrackId() )
2261 // else if( xDigiA->GetLinks() != xDigiB->GetLinks() )
2262 // if other side come from a different TofPoint,
2263 // also save the pointer on CbmTofPoint
2264 // vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) );
2265 } // else of if( 0 < iNbChanInHit)
2266 iLastChan = iCh;
2267 dLastPosX = dPosX;
2268 dLastPosY = dPosY;
2269 dLastTime = dTime;
2270 if (AddNextChan(iSmType, iSm, iRpc, iLastChan, dLastPosX, dLastPosY, dLastTime, dWeightsSum)) {
2271 iNbChanInHit = 0; // cluster already stored
2272 }
2273 } // while( 1 < fStorDigi[iSmType][iSm*iNbRpc+iRpc][iCh].size() )
2274 fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].clear();
2275 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].clear();
2276 } // for( Int_t iCh = 0; iCh < iNbCh; iCh++ )
2277 LOG(debug2) << "finished V-RPC"
2278 << Form(" %3d %3d %3d %d %f %fx", iSmType, iSm, iRpc, fTofHitsColl->GetEntriesFast(),
2279 dLastPosX, dLastPosY);
2280 } // else of if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
2281 } // if( 0 == iChType)
2282 else {
2283 LOG(error) << "=> Cluster building "
2284 << "from digis to hits not implemented for pads, Sm type " << iSmType << " Rpc " << iRpc;
2285 return kFALSE;
2286 } // else of if( 0 == iChType)
2287
2288 // Now check if another hit/cluster is started
2289 // and save it if it's the case
2290 if (0 < iNbChanInHit) {
2291 LOG(debug1) << "Process cluster " << iNbChanInHit;
2292
2293 // Check orientation to properly assign errors
2294 if (1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc)) {
2295 LOG(debug1) << "H-Hit ";
2296 } // if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
2297 else {
2298 LOG(debug2) << "V-Hit ";
2299 // Save Hit
2300 dWeightedTime /= dWeightsSum;
2301 dWeightedPosX /= dWeightsSum;
2302 dWeightedPosY /= dWeightsSum;
2303 dWeightedPosZ /= dWeightsSum;
2304 //TVector3 hitPos(dWeightedPosX, dWeightedPosY, dWeightedPosZ);
2305
2306 Double_t hitpos_local[3] = {3 * 0.};
2307 hitpos_local[0] = dWeightedPosX;
2308 hitpos_local[1] = dWeightedPosY;
2309 hitpos_local[2] = dWeightedPosZ;
2310 /*
2311 if( 5 == iSmType || 8 == iSmType) { // for PAD counters
2312 hitpos_local[0] = (gRandom->Rndm()-0.5)*fChannelInfo->GetSizex();
2313 hitpos_local[1] = (gRandom->Rndm()-0.5)*fChannelInfo->GetSizey();
2314 }
2315 */
2316 Double_t hitpos[3] = {3 * 0.};
2317 if (5 != iSmType) {
2318 /*TGeoNode* cNode=*/gGeoManager->GetCurrentNode();
2319 /*TGeoHMatrix* cMatrix =*/gGeoManager->GetCurrentMatrix();
2320 //cNode->Print();
2321 //cMatrix->Print();
2322 gGeoManager->LocalToMaster(hitpos_local, hitpos);
2323 }
2324 LOG(debug1) << Form(" LocalToMaster for V-node: "
2325 "(%6.1f,%6.1f,%6.1f) ->(%6.1f,%6.1f,%6.1f)",
2326 hitpos_local[0], hitpos_local[1], hitpos_local[2], hitpos[0], hitpos[1], hitpos[2]);
2327
2328 TVector3 hitPos(hitpos[0], hitpos[1], hitpos[2]);
2329 // Event errors, not properly done at all for now
2330 // Right way of doing it should take into account the weight distribution
2331 // and real system time resolution
2332 TVector3 hitPosErr(0.5, 0.5, 0.5); // including positioning uncertainty
2333 /*
2334 TVector3 hitPosErr( fChannelInfo->GetSizex()/TMath::Sqrt(12.0), // Single strips approximation
2335 0.5, // Use generic value
2336 1.);
2337 */
2338 // fDigiBdfPar->GetFeeTimeRes() * fDigiBdfPar->GetSigVel(iSmType,iRpc), // Use the electronics resolution
2339 // fDigiBdfPar->GetNbGaps( iSmType, iRpc)*
2340 // fDigiBdfPar->GetGapSize( iSmType, iRpc)/10.0 / // Change gap size in cm
2341 // TMath::Sqrt(12.0) ); // Use full RPC thickness as "Channel" Z size
2342 // cout<<"a "<<vPtsRef.size()<<endl;
2343 // cout<<"b "<<vPtsRef[0]<<endl;
2344 // cout<<"c "<<vPtsRef[0]->GetDetectorID()<<endl;
2345 // Int_t iDetId = vPtsRef[0]->GetDetectorID();// detID = pt->GetDetectorID() <= from TofPoint
2346 // Int_t iDetId = iChId;
2347 Int_t iChm = floor(dWeightedPosX / fChannelInfo->GetSizex()) + iNbCh / 2;
2348 if (iChm < 0) iChm = 0;
2349 if (iChm > iNbCh - 1) iChm = iNbCh - 1;
2350 iDetId = CbmTofAddress::GetUniqueAddress(iSm, iRpc, iChm, 0, iSmType);
2351 Int_t iRefId = 0; // Index of the correspondng TofPoint
2352 //if(NULL != fTofPointsColl) iRefId = fTofPointsColl->IndexOf( vPtsRef[0] );
2353 TString cstr = "Save V-Hit ";
2354 cstr += Form(" %3d %3d 0x%08x %3d 0x%08x %8.2f %6.2f", // %3d %3d
2355 fiNbHits, iNbChanInHit, iDetId, iLastChan,
2356 iRefId, //vPtsRef.size(),vPtsRef[0])
2357 dWeightedTime, dWeightedPosY);
2358
2359 cstr += Form(", DigiSize: %lu ", vDigiIndRef.size());
2360 cstr += ", DigiInds: ";
2361
2362 fviClusterMul[iSmType][iSm][iRpc]++;
2363
2364 for (UInt_t i = 0; i < vDigiIndRef.size(); i++) {
2365 cstr += Form(" %d (M,%d)", vDigiIndRef.at(i), fviClusterMul[iSmType][iSm][iRpc]);
2366 }
2367 LOG(debug) << cstr;
2368
2369 if (vDigiIndRef.size() < 2) {
2370 LOG(warning) << "Digi refs for Hit " << fiNbHits << ": vDigiIndRef.size()";
2371 }
2372 if (fiNbHits > 0) {
2373 CbmTofHit* pHitL = (CbmTofHit*) fTofHitsColl->At(fiNbHits - 1);
2374 if (iDetId == pHitL->GetAddress() && dWeightedTime == pHitL->GetTime())
2375 LOG(debug) << "Store Hit twice? "
2376 << " fiNbHits " << fiNbHits << ", " << Form("0x%08x", iDetId);
2377 }
2378
2379 CbmTofHit* pHit = new CbmTofHit(iDetId, hitPos,
2380 hitPosErr, //local detector coordinates
2381 fiNbHits, // this number is used as reference!!
2382 dWeightedTime,
2383 vDigiIndRef.size(), // number of linked digis = 2*CluSize
2384 //vPtsRef.size(), // flag = number of TofPoints generating the cluster
2385 Int_t(dWeightsSum * 10.)); //channel -> Tot
2386 // 0) ; //channel
2387 // vDigiIndRef);
2388 // output hit
2389 new ((*fTofHitsColl)[fiNbHits]) CbmTofHit(*pHit);
2390 // memorize hit
2391 if (fdMemoryTime > 0.) {
2392 LH_store(iSmType, iSm, iRpc, iChm, pHit);
2393 }
2394 else {
2395 pHit->Delete();
2396 }
2397 /*
2398 new((*fTofDigiMatchColl)[fiNbHits]) CbmMatch();
2399 CbmMatch* digiMatch = (CbmMatch *)fTofDigiMatchColl->At(fiNbHits);
2400 */
2401 CbmMatch* digiMatch = new ((*fTofDigiMatchColl)[fiNbHits]) CbmMatch();
2402
2403 for (size_t i = 0; i < vDigiIndRef.size(); i++) {
2404 Double_t dTot = fTofCalDigiVec->at(vDigiIndRef.at(i)).GetTot();
2405 digiMatch->AddLink(CbmLink(dTot, vDigiIndRef.at(i), fiOutputTreeEntry, fiFileIndex));
2406 }
2407
2408 fiNbHits++;
2409 // For Histogramming
2410 fviClusterSize[iSmType][iRpc].push_back(iNbChanInHit);
2411 //fviTrkMul[iSmType][iRpc].push_back( vPtsRef.size() );
2412 fvdX[iSmType][iRpc].push_back(dWeightedPosX);
2413 fvdY[iSmType][iRpc].push_back(dWeightedPosY);
2414 /*
2415 fvdDifX[iSmType][iRpc].push_back( vPtsRef[0]->GetX() - dWeightedPosX);
2416 fvdDifY[iSmType][iRpc].push_back( vPtsRef[0]->GetY() - dWeightedPosY);
2417 fvdDifCh[iSmType][iRpc].push_back( fGeoHandler->GetCell( vPtsRef[0]->GetDetectorID() ) -1 -iLastChan );
2418 */
2419 //vPtsRef.clear();
2420 vDigiIndRef.clear();
2421 } // else of if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
2422 } // if( 0 < iNbChanInHit)
2423 LOG(debug2) << " Fini-A " << Form(" %3d %3d %3d M%3d", iSmType, iSm, iRpc, fviClusterMul[iSmType][iSm][iRpc]);
2424 } // for each sm/rpc pair
2425 LOG(debug2) << " Fini-B " << Form(" %3d ", iSmType);
2426 } // for( Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++ )
2427 }
2428 return kTRUE;
2429}
2430
2432{
2433 CbmTofDigi* pDigi;
2434 CbmTofDigi* pCalDigi = NULL;
2435 Int_t iDigIndCal = -1;
2436 // channel deadtime map
2437 std::map<Int_t, Double_t> mChannelDeadTime;
2438
2439 Int_t iNbTofDigi = fTofDigiVec.size();
2440 //Int_t iNbTofDigi = fTofDigisColl->GetEntriesFast();
2441 for (Int_t iDigInd = 0; iDigInd < iNbTofDigi; iDigInd++) {
2442 pDigi = &(fTofDigiVec.at(iDigInd));
2443 //pDigi = (CbmTofDigi*) fTofDigisColl->At( iDigInd );
2444
2445 if (fbSwapChannelSides && 5 != pDigi->GetType() && 8 != pDigi->GetType()) {
2446 pDigi->SetAddress(pDigi->GetSm(), pDigi->GetRpc(), pDigi->GetChannel(), (0 == pDigi->GetSide()) ? 1 : 0,
2447 pDigi->GetType());
2448 }
2449
2450 Int_t iAddr = pDigi->GetAddress();
2451
2452 LOG(debug1) << "BC " // Before Calibration
2453 << Form("0x%08x", pDigi->GetAddress()) << " TSRC " << pDigi->GetType() << pDigi->GetSm()
2454 << pDigi->GetRpc() << Form("%2d", (Int_t) pDigi->GetChannel()) << " " << pDigi->GetSide() << " "
2455 << Form("%f", pDigi->GetTime()) << " " << pDigi->GetTot();
2456 /*
2457 if (pDigi->GetType() == 5
2458 || pDigi->GetType()
2459 == 8) // for Pad counters generate fake digi to mockup a strip
2460 if (pDigi->GetSide() == 1)
2461 continue; // skip one side to avoid double entries
2462 */
2463 Bool_t bValid = kTRUE;
2464 std::map<Int_t, Double_t>::iterator it;
2465 it = mChannelDeadTime.find(iAddr);
2466 if (it != mChannelDeadTime.end()) {
2467 LOG(debug1) << "CCT found valid ChannelDeadtime entry " << mChannelDeadTime[iAddr] << ", DeltaT "
2468 << pDigi->GetTime() - mChannelDeadTime[iAddr];
2469 if ((bValid = (pDigi->GetTime() > mChannelDeadTime[iAddr] + fdChannelDeadtime))) {
2470
2471 // pCalDigi = new((*fTofCalDigisColl)[++iDigIndCal]) CbmTofDigi( *pDigi );
2472 fTofCalDigiVec->push_back(CbmTofDigi(*pDigi));
2473 pCalDigi = &(fTofCalDigiVec->back());
2474 iDigIndCal++;
2475 }
2476 }
2477 else {
2478 fTofCalDigiVec->push_back(CbmTofDigi(*pDigi));
2479 pCalDigi = &(fTofCalDigiVec->back());
2480 iDigIndCal++;
2481 // pCalDigi = new((*fTofCalDigisColl)[++iDigIndCal]) CbmTofDigi( *pDigi );
2482 }
2483 mChannelDeadTime[iAddr] = pDigi->GetTime();
2484 if (!bValid) continue;
2485
2486 LOG(debug1) << "DC " // After deadtime check. before Calibration
2487 << Form("0x%08x", pDigi->GetAddress()) << " TSRC " << pDigi->GetType() << pDigi->GetSm()
2488 << pDigi->GetRpc() << Form("%2d", (Int_t) pDigi->GetChannel()) << " " << pDigi->GetSide() << " "
2489 << Form("%f", pDigi->GetTime()) << " " << pDigi->GetTot();
2490
2491 if (fbPs2Ns) {
2492 pCalDigi->SetTime(pCalDigi->GetTime() / 1000.); // for backward compatibility
2493 pCalDigi->SetTot(pCalDigi->GetTot() / 1000.); // for backward compatibility
2494 }
2495
2496 if (fDigiBdfPar->GetNbSmTypes() > pDigi->GetType() // prevent crash due to misconfiguration
2497 && fDigiBdfPar->GetNbSm(pDigi->GetType()) > pDigi->GetSm()
2498 && fDigiBdfPar->GetNbRpc(pDigi->GetType()) > pDigi->GetRpc()
2499 && fDigiBdfPar->GetNbChan(pDigi->GetType(), pDigi->GetRpc()) > pDigi->GetChannel()) {
2500
2501 LOG(debug2) << " CluCal-Init: " << pDigi->ToString();
2502 // apply calibration vectors
2503 pCalDigi->SetTime(pCalDigi->GetTime() - // calibrate Digi Time
2504 fvCPTOff[pDigi->GetType()][pDigi->GetSm() * fDigiBdfPar->GetNbRpc(pDigi->GetType())
2505 + pDigi->GetRpc()][pDigi->GetChannel()][pDigi->GetSide()]);
2506 LOG(debug2) << " CluCal-TOff: " << pCalDigi->ToString();
2507
2508 Double_t dTot = pCalDigi->GetTot() + fRndm->Uniform(0, 1) - // subtract Offset
2509 fvCPTotOff[pDigi->GetType()][pDigi->GetSm() * fDigiBdfPar->GetNbRpc(pDigi->GetType())
2510 + pDigi->GetRpc()][pDigi->GetChannel()][pDigi->GetSide()];
2511 if (dTot < 0.001) dTot = 0.001;
2512 pCalDigi->SetTot(dTot * // calibrate Digi ToT
2513 fvCPTotGain[pDigi->GetType()][pDigi->GetSm() * fDigiBdfPar->GetNbRpc(pDigi->GetType())
2514 + pDigi->GetRpc()][pDigi->GetChannel()][pDigi->GetSide()]);
2515
2516 // walk correction
2517 Double_t dTotBinSize = (fdTOTMax - fdTOTMin) / nbClWalkBinX;
2518 Int_t iWx = (Int_t)((pCalDigi->GetTot() - fdTOTMin) / dTotBinSize);
2519 if (0 > iWx) iWx = 0;
2520 if (iWx >= nbClWalkBinX) iWx = nbClWalkBinX - 1;
2521 Double_t dDTot = (pCalDigi->GetTot() - fdTOTMin) / dTotBinSize - (Double_t) iWx - 0.5;
2522 Double_t dWT =
2523 fvCPWalk[pCalDigi->GetType()][pCalDigi->GetSm() * fDigiBdfPar->GetNbRpc(pCalDigi->GetType())
2524 + pCalDigi->GetRpc()][pCalDigi->GetChannel()][pCalDigi->GetSide()][iWx];
2525 if (dDTot > 0) { // linear interpolation to next bin
2526 if (iWx < nbClWalkBinX - 1) { // linear interpolation to next bin
2527
2528 dWT +=
2529 dDTot
2530 * (fvCPWalk[pCalDigi->GetType()][pCalDigi->GetSm() * fDigiBdfPar->GetNbRpc(pCalDigi->GetType())
2531 + pCalDigi->GetRpc()][pCalDigi->GetChannel()][pCalDigi->GetSide()][iWx + 1]
2532 - fvCPWalk[pCalDigi->GetType()]
2533 [pCalDigi->GetSm() * fDigiBdfPar->GetNbRpc(pCalDigi->GetType()) + pCalDigi->GetRpc()]
2534 [pCalDigi->GetChannel()][pCalDigi->GetSide()][iWx]); //memory leak???
2535 }
2536 }
2537 else // dDTot < 0, linear interpolation to next bin
2538 {
2539 if (0 < iWx) { // linear interpolation to next bin
2540 dWT -=
2541 dDTot
2542 * (fvCPWalk[pCalDigi->GetType()][pCalDigi->GetSm() * fDigiBdfPar->GetNbRpc(pCalDigi->GetType())
2543 + pCalDigi->GetRpc()][pCalDigi->GetChannel()][pCalDigi->GetSide()][iWx - 1]
2544 - fvCPWalk[pCalDigi->GetType()]
2545 [pCalDigi->GetSm() * fDigiBdfPar->GetNbRpc(pCalDigi->GetType()) + pCalDigi->GetRpc()]
2546 [pCalDigi->GetChannel()][pCalDigi->GetSide()][iWx]); //memory leak???
2547 }
2548 }
2549
2550 pCalDigi->SetTime(pCalDigi->GetTime() - dWT); // calibrate Digi Time
2551 LOG(debug2) << " CluCal-Walk: " << pCalDigi->ToString();
2552
2553 if (1) { //pDigi->GetType()==7 && pDigi->GetSm()==0){
2554 LOG(debug)
2555 << "BuildClusters: CalDigi " << Form("%02d TSRCS ", iDigIndCal) << pCalDigi->GetType() << pCalDigi->GetSm()
2556 << pCalDigi->GetRpc() << Form("%02d", Int_t(pCalDigi->GetChannel())) << pCalDigi->GetSide()
2557 << Form(", T %15.3f", pCalDigi->GetTime()) << ", Tot " << pCalDigi->GetTot() << ", TotGain "
2558 << fvCPTotGain[pCalDigi->GetType()][pCalDigi->GetSm() * fDigiBdfPar->GetNbRpc(pCalDigi->GetType())
2559 + pCalDigi->GetRpc()][pCalDigi->GetChannel()][pCalDigi->GetSide()]
2560 << ", TotOff "
2561 << fvCPTotOff[pCalDigi->GetType()][pCalDigi->GetSm() * fDigiBdfPar->GetNbRpc(pCalDigi->GetType())
2562 + pCalDigi->GetRpc()][pCalDigi->GetChannel()][pCalDigi->GetSide()]
2563 << ", Walk " << iWx << ": "
2564 << fvCPWalk[pCalDigi->GetType()][pCalDigi->GetSm() * fDigiBdfPar->GetNbRpc(pCalDigi->GetType())
2565 + pCalDigi->GetRpc()][pCalDigi->GetChannel()][pCalDigi->GetSide()][iWx];
2566
2567 LOG(debug1)
2568 << " dDTot " << dDTot << " BinSize: " << dTotBinSize << ", CPWalk "
2569 << fvCPWalk[pCalDigi->GetType()][pCalDigi->GetSm() * fDigiBdfPar->GetNbRpc(pCalDigi->GetType())
2570 + pCalDigi->GetRpc()][pCalDigi->GetChannel()][pCalDigi->GetSide()][iWx - 1]
2571 << ", "
2572 << fvCPWalk[pCalDigi->GetType()][pCalDigi->GetSm() * fDigiBdfPar->GetNbRpc(pCalDigi->GetType())
2573 + pCalDigi->GetRpc()][pCalDigi->GetChannel()][pCalDigi->GetSide()][iWx]
2574 << ", "
2575 << fvCPWalk[pCalDigi->GetType()][pCalDigi->GetSm() * fDigiBdfPar->GetNbRpc(pCalDigi->GetType())
2576 + pCalDigi->GetRpc()][pCalDigi->GetChannel()][pCalDigi->GetSide()][iWx + 1]
2577 << " -> dWT = " << dWT;
2578 }
2579 }
2580 else {
2581 LOG(info) << "Skip1 Digi "
2582 << " Type " << pDigi->GetType() << " " << fDigiBdfPar->GetNbSmTypes() << " Sm " << pDigi->GetSm() << " "
2583 << fDigiBdfPar->GetNbSm(pDigi->GetType()) << " Rpc " << pDigi->GetRpc() << " "
2584 << fDigiBdfPar->GetNbRpc(pDigi->GetType()) << " Ch " << pDigi->GetChannel() << " "
2585 << fDigiBdfPar->GetNbChan(pDigi->GetType(), 0);
2586 }
2587
2588 if (0) // (bAddBeamCounterSideDigi)
2589 if (pCalDigi->GetType() == 5
2590 || pCalDigi->GetType() == 8) { // for Pad counters generate fake digi to mockup a strip
2591 LOG(debug) << "add Pad counter 2. Side digi for 0x" << std::hex << pCalDigi->GetAddress();
2592 fTofCalDigiVec->push_back(CbmTofDigi(*pCalDigi));
2593 CbmTofDigi* pCalDigi2 = &(fTofCalDigiVec->back());
2594 iDigIndCal++;
2595 // CbmTofDigi *pCalDigi2 = new((*fTofCalDigisColl)[++iDigIndCal]) CbmTofDigi( *pCalDigi );
2596 if (pCalDigi->GetSide() == 0)
2597 pCalDigi2->SetAddress(pCalDigi->GetSm(), pCalDigi->GetRpc(), pCalDigi->GetChannel(), 1, pCalDigi->GetType());
2598 else
2599 pCalDigi2->SetAddress(pCalDigi->GetSm(), pCalDigi->GetRpc(), pCalDigi->GetChannel(), 0, pCalDigi->GetType());
2600 }
2601
2602 } // for( Int_t iDigInd = 0; iDigInd < nTofDigi; iDigInd++ )
2603
2604 // iNbTofDigi = fTofCalDigisColl->GetEntriesFast(); // update because of added duplicted digis
2605 iNbTofDigi = fTofCalDigiVec->size(); // update because of added duplicted digis
2606 //if(fTofCalDigisColl->IsSortable())
2607 // LOG(debug)<<"CbmTofHitMaker::BuildClusters: Sort "<<fTofCalDigisColl->GetEntriesFast()<<" calibrated digis ";
2608 LOG(debug) << "CbmTofHitMaker::BuildClusters: Sort " << fTofCalDigiVec->size() << " calibrated digis ";
2609 if (iNbTofDigi > 1) {
2610 // fTofCalDigisColl->Sort(iNbTofDigi); // Time order again, in case modified by the calibration
2612 std::sort(fTofCalDigiVec->begin(), fTofCalDigiVec->end(),
2613 [](const CbmTofDigi& a, const CbmTofDigi& b) -> bool { return a.GetTime() < b.GetTime(); });
2614 // std::sort(fTofCalDigiVec->begin(), fTofCalDigiVec->end());
2615 // if(!fTofCalDigisColl->IsSorted()){
2616 // if ( ! std::is_sorted(fTofCalDigiVec->begin(), fTofCalDigiVec->end()))
2617 if (!std::is_sorted(fTofCalDigiVec->begin(), fTofCalDigiVec->end(),
2618 [](const CbmTofDigi& a, const CbmTofDigi& b) -> bool { return a.GetTime() < b.GetTime(); }))
2619 LOG(warning) << "CbmTofHitMaker::BuildClusters: Sorting not successful ";
2620 }
2621 // }
2622
2623 return kTRUE;
2624}
2625
2626void CbmTofHitMaker::SetDeadStrips(Int_t iDet, Int_t ival)
2627{
2628 if (fvDeadStrips.size() < static_cast<size_t>(iDet + 1)) fvDeadStrips.resize(iDet + 1);
2629 fvDeadStrips[iDet] = ival;
2630}
@ kTof
Time-of-flight Detector.
static FairRootManager * rootMgr
static Int_t iNbTs
static TFile * fHist
TClonesArray * fTofHitsColl
TClonesArray * fEventsColl
const Int_t DetMask
static Double_t fdMemoryTime
CbmDigiManager * fDigiMan
std::vector< TH1 * > hSvel
const Int_t iNSel
const Int_t nbClWalkBinX
const Int_t nbClDelTofBinX
static TRandom3 * fRndm
static Bool_t bAddBeamCounterSideDigi
@ k14a
static Int_t iNbTs
static Double_t f1_xboxe(double *x, double *par)
static TRandom3 * fRndm
static Bool_t bAddBeamCounterSideDigi
std::vector< Int_t > vDigiIndRef
static constexpr size_t size()
Definition KfSimdPseudo.h:2
Data class for a signal in the t-zero detector.
Definition CbmBmonDigi.h:30
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
size_t GetNofData() const
Definition CbmEvent.cxx:53
uint32_t GetIndex(ECbmDataType type, uint32_t iData)
Definition CbmEvent.cxx:42
void AddData(ECbmDataType type, uint32_t index)
Definition CbmEvent.cxx:33
double GetTime() const
Definition CbmHit.h:76
int32_t GetAddress() const
Definition CbmHit.h:74
void SetTime(double time)
Definition CbmHit.h:85
const CbmLink & GetLink(int32_t i) const
Definition CbmMatch.h:39
int32_t GetNofLinks() const
Definition CbmMatch.h:42
void AddLink(const CbmLink &newLink)
Definition CbmMatch.cxx:47
void SetX(double x)
Definition CbmPixelHit.h:92
double GetY() const
Definition CbmPixelHit.h:74
void SetY(double y)
Definition CbmPixelHit.h:93
double GetX() const
Definition CbmPixelHit.h:73
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)
static int32_t GetSmId(uint32_t address)
static int32_t GetRpcId(uint32_t address)
static int32_t GetSmType(uint32_t address)
static int32_t GetChannelId(uint32_t address)
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.
Int_t GetNbSmTypes() const
void SetSigVel(Int_t iSmType, Int_t iSm, Int_t iRpc, Double_t dvel)
Int_t GetNbSm(Int_t iSmType) const
Double_t GetMaxDistAlongCh() const
Int_t GetNbChan(Int_t iSmType, Int_t iRpc) const
Int_t GetDetInd(Int_t iAddr)
Bool_t ClustUseTrackId() 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
Int_t GetDetUId(Int_t iDet)
Double_t GetSignalSpeed() const
Double_t GetSigVel(Int_t iSmType, Int_t iSm, Int_t iRpc) const
CbmTofCell * GetCell(Int_t i)
TGeoNode * GetNode(Int_t iCell)
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
std::string ToString() const
void SetTot(double tot)
Definition CbmTofDigi.h:166
double GetChannel() const
Channel .
Definition CbmTofDigi.h:156
int32_t GetAddress() const
Inherited from CbmDigi.
Definition CbmTofDigi.h:112
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 SetAddress(int32_t address)
Definition CbmTofDigi.h:163
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)
std::vector< std::vector< std::vector< std::vector< Double_t > > > > fvCPTOff
TClonesArray * fEventsColl
Double_t fdChannelDeadtime
CbmTofDigiBdfPar * fDigiBdfPar
virtual InitStatus Init()
Inherited from FairTask.
virtual void Exec(Option_t *option)
Inherited from FairTask.
void SetDeadStrips(Int_t iDet, Int_t ival)
TClonesArray * fTofDigiMatchCollOut
Double_t fdCaldXdYMax
Bool_t RegisterInputs()
Recover pointer on input TClonesArray: TofPoints, TofDigis...
std::map< UInt_t, UInt_t > fDetIdIndexMap
TClonesArray * fTofDigiMatchColl
TClonesArray * fTofHitsCollOut
// Calibrated TOF Digis
std::vector< std::vector< std::vector< std::vector< Int_t > > > > fStorDigiInd
virtual void fit_ybox(const char *hname)
virtual void CheckLHMemory()
Bool_t fbSwapChannelSides
void SetCalMode(Int_t iMode)
Bool_t LoadGeometry()
Load the geometry: for now just resizing the Digis temporary vectors.
Double_t fPosYMaxScal
TTrbHeader * fTrbHeader
TString fOutHstFileName
std::vector< std::vector< std::vector< Int_t > > > fviTrkMul
CbmTofDetectorId * fTofId
std::vector< std::vector< std::vector< Double_t > > > fvdDifY
std::vector< std::vector< std::vector< std::vector< Double_t > > > > fvCPDelTof
std::vector< Int_t > fvDeadStrips
std::vector< Int_t > fviDetId
CbmTofCell * fChannelInfo
Bool_t InitParameters()
Initialize other parameters not included in parameter classes.
std::vector< std::vector< std::vector< Double_t > > > fvdY
virtual ~CbmTofHitMaker()
Destructor.
std::vector< CbmTofDigi > * fTofCalDigiVecOut
std::vector< std::vector< std::vector< std::vector< std::list< CbmTofHit * > > > > > fvLastHits
std::vector< Int_t > vDigiIndRef
std::vector< TH2 * > fhRpcDigiMul
std::vector< std::vector< std::vector< std::vector< std::vector< Double_t > > > > > fvCPWalk
TString fCalParFileName
std::vector< std::vector< std::vector< Double_t > > > fvdDifCh
std::vector< std::vector< Double_t > > fvTimeFirstDigi
std::vector< CbmTofDigi > * fTofCalDigiVec
CbmTofHitMaker()
Constructor.
CbmDigiManager * fDigiMan
TOF Digis.
std::vector< TH2 * > fhRpcDigiDTMul
virtual void SetParContainers()
Inherited from FairTask.
virtual void CleanLHMemory()
std::vector< std::vector< Double_t > > fvTimeLastDigi
std::vector< std::vector< std::vector< Int_t > > > fviClusterSize
virtual Bool_t AddNextChan(Int_t iSmType, Int_t iSm, Int_t iRpc, Int_t iLastChan, Double_t dLastPosX, Double_t dLastPosY, Double_t dLastTime, Double_t dLastTot)
CbmTofDigiPar * fDigiPar
static CbmTofHitMaker * fInstance
Bool_t fbAlternativeBranchNames
virtual void Finish()
Inherited from FairTask.
TTimeStamp fStart
Double_t fMaxTimeDist
Double_t fdMaxSpaceDist
std::vector< CbmTofDigi > fTofDigiVec
std::vector< std::vector< std::vector< Double_t > > > fvdDifX
std::vector< std::vector< Double_t > > fvMulDigi
std::vector< std::vector< std::vector< Double_t > > > fvdX
virtual void LH_store(Int_t iSmType, Int_t iSm, Int_t iRpc, Int_t iChm, CbmTofHit *pHit)
std::vector< std::vector< std::vector< std::vector< CbmTofDigi * > > > > fStorDigi
CbmTofGeoHandler * fGeoHandler
std::vector< TH2 * > fhRpcDigiStatus
std::vector< std::vector< std::vector< Int_t > > > fviClusterMul
TClonesArray * fTofHitsColl
// Calibrated TOF Digis
virtual void ExecEvent(Option_t *option)
std::vector< std::vector< std::vector< std::vector< Double_t > > > > fvCPTotOff
std::vector< std::vector< std::vector< std::vector< Double_t > > > > fvCPTotGain
Bool_t BuildClusters()
Build clusters out of ToF Digis and store the resulting info in a TofHit.
Bool_t InitCalibParameter()
Initialize other parameters not included in parameter classes.
Bool_t RegisterOutputs()
Create and register output TClonesArray of Tof Hits.
Double_t fdMaxTimeDist
Bool_t DeleteGeometry()
Delete the geometry related arrays: for now just clearing the Digis temporary vectors.
Double_t fdMemoryTime
Double_t f1_xboxe(double *x, double *par)
Definition fit_ybox.h:6