BmnRoot
Loading...
Searching...
No Matches
BmnSiliconHitMaker.cxx
Go to the documentation of this file.
1#include "BmnSiliconHitMaker.h"
2
3#include "BmnSiliconModule.h"
4#include "FairRootManager.h"
5#include "TClonesArray.h"
6#include "TSystem.h"
7#include "UniDetectorParameter.h"
8
9#include <TStopwatch.h>
10
11static Double_t workTime = 0.0;
12
13BmnSiliconHitMaker::BmnSiliconHitMaker()
14 : fHitMatching(kTRUE)
15{
16
17 fInputPointsBranchName = "SiliconPoint";
18 fInputDigitsBranchName = "BmnSiliconDigit";
19 fInputDigitMatchesBranchName = "BmnSiliconDigitMatch";
20
21 fOutputHitsBranchName = "BmnSiliconHit";
22
23 fCurrentConfig = BmnSiliconConfiguration::None;
24 StationSet = nullptr;
25 TransfSet = nullptr;
26
27 fBmnSiliconHitsArray = nullptr;
28 fBmnSiliconUpperClustersArray = nullptr;
29 fBmnSiliconLowerClustersArray = nullptr;
30}
31
32BmnSiliconHitMaker::BmnSiliconHitMaker(Int_t run_period, Int_t run_number, Bool_t isExp, Bool_t isSrc)
33 : fHitMatching(kTRUE)
34{
35 fIsExp = isExp;
36 fIsSrc = isSrc;
37
38 fRunPeriod = run_period;
39 fRunNumber = run_number;
40
41 fInputPointsBranchName = "SiliconPoint";
42 fInputDigitsBranchName = (!isExp) ? "BmnSiliconDigit" : "SILICON";
43 fInputDigitMatchesBranchName = "BmnSiliconDigitMatch";
44
45 fOutputHitsBranchName = "BmnSiliconHit";
46
47 fCurrentConfig = BmnSiliconConfiguration::None;
48 StationSet = nullptr;
49 TransfSet = nullptr;
50
51 fBmnSiliconHitsArray = nullptr;
52 fBmnSiliconUpperClustersArray = nullptr;
53 fBmnSiliconLowerClustersArray = nullptr;
54
55 fSignalLow = 0.;
56 fSignalUp = DBL_MAX;
57
58 switch (fRunPeriod) {
59 case 6: // BM@N RUN-6
60 fCurrentConfig = BmnSiliconConfiguration::RunSpring2017;
61 break;
62 case 7: // BM@N RUN-7 (and SRC)
63 if (fIsSrc) {
64 fCurrentConfig = BmnSiliconConfiguration::RunSRCSpring2018;
65 } else {
66 fCurrentConfig = BmnSiliconConfiguration::RunSpring2018;
67 }
68 break;
69 case 8: // BM@N RUN-8
70 fCurrentConfig = BmnSiliconConfiguration::Run8;
71 break;
72 case 9: // BM@N RUN-9
73 fCurrentConfig = BmnSiliconConfiguration::Run9;
74 break;
75
76 default:
77 fCurrentConfig = BmnSiliconConfiguration::None;
78 }
79}
80
81void BmnSiliconHitMaker::createSiliconDetector()
82{
83 TString gPathSiliconConfig = gSystem->Getenv("VMCWORKDIR");
84 gPathSiliconConfig += "/parameters/silicon/XMLConfigs/";
85
86 // Create SILICON detector --------------------------------------------------
87 switch (fCurrentConfig) {
88 case BmnSiliconConfiguration::RunSpring2017:
89 XMLConfigFile = gPathSiliconConfig + "SiliconRunSpring2017.xml";
90 if (fVerbose)
91 cout << " Current SILICON Configuration : RunSpring2017"
92 << "\n";
93 break;
94
95 case BmnSiliconConfiguration::RunSpring2018:
96 XMLConfigFile = gPathSiliconConfig + "SiliconRunSpring2018.xml";
97 if (fVerbose)
98 cout << " Current SILICON Configuration : RunSpring2018"
99 << "\n";
100 break;
101
102 case BmnSiliconConfiguration::RunSRCSpring2018:
103 XMLConfigFile = gPathSiliconConfig + "SiliconRunSRCSpring2018.xml";
104 if (fVerbose)
105 cout << " Current SILICON Configuration : RunSRCSpring2018"
106 << "\n";
107 break;
108
109 case BmnSiliconConfiguration::Run8_3stations:
110 XMLConfigFile = gPathSiliconConfig + "SiliconRun8_3stations.xml";
111 if (fVerbose)
112 cout << " Current SILICON Configuration : SiliconRun8_3stations"
113 << "\n";
114 break;
115
116 case BmnSiliconConfiguration::Run8_4stations:
117 XMLConfigFile = gPathSiliconConfig + "SiliconRun8_4stations.xml";
118 if (fVerbose)
119 cout << " Current SILICON Configuration : SiliconRun8_4stations"
120 << "\n";
121 break;
122
123 case BmnSiliconConfiguration::Run8_5stations:
124 XMLConfigFile = gPathSiliconConfig + "SiliconRun8_5stations.xml";
125 if (fVerbose)
126 cout << " Current SILICON Configuration : SiliconRun8_5stations"
127 << "\n";
128 break;
129
130 case BmnSiliconConfiguration::Run8_mods_6_10_14_18:
131 XMLConfigFile = gPathSiliconConfig + "SiliconRun8_mods_6_10_14_18.xml";
132 if (fVerbose)
133 cout << " Current SILICON Configuration : SiliconRun8_mods_6_10_14_18"
134 << "\n";
135 break;
136
137 case BmnSiliconConfiguration::Run8_mods_10_14rot_18:
138 XMLConfigFile = gPathSiliconConfig + "SiliconRun8_mods_10_14rot_18.xml";
139 if (fVerbose)
140 cout << " Current SILICON Configuration : SiliconRun8_mods_10_14rot_18"
141 << "\n";
142 break;
143
144 case BmnSiliconConfiguration::Run8:
145 XMLConfigFile = gPathSiliconConfig + "SiliconRun8.xml";
146 if (fVerbose)
147 cout << " Current SILICON Configuration : Run8"
148 << "\n";
149 break;
150
151 case BmnSiliconConfiguration::Run9:
152 XMLConfigFile = gPathSiliconConfig + "SiliconRun9.xml";
153 if (fVerbose)
154 cout << " Current SILICON Configuration : Run9"
155 << "\n";
156 break;
157
158 case BmnSiliconConfiguration::Run_5stations_2026:
159 XMLConfigFile = gPathSiliconConfig + "Silicon_5stations_2026.xml";
160 if (fVerbose)
161 cout << " Current SILICON Configuration : 5stations_2026"
162 << "\n";
163 break;
164
165 default:;
166 }
167
168 if (!gSystem->AccessPathName(XMLConfigFile)) { // returns FALSE if the file exists
169 if (fVerbose)
170 cout << " XMLConfigFile : " << XMLConfigFile << "\n";
171 StationSet = new BmnSiliconStationSet(XMLConfigFile);
172 TransfSet = new BmnSiliconTransform();
173 TransfSet->LoadFromXMLFile(XMLConfigFile);
174 } else {
175 Fatal("BmnSiliconHitMaker::createSiliconDetector()", " !!! Current configuration is not set !!! ");
176 }
177
178 if (fIsExp) {
179 const Int_t nStat = StationSet->GetNStations();
180 UniDetectorParameter* coeffAlignCorrs =
181 UniDetectorParameter::GetDetectorParameter("Silicon", "alignment_shift", fRunPeriod, fRunNumber);
182 vector<UniValue*> algnShifts;
183 if (coeffAlignCorrs)
184 coeffAlignCorrs->GetValue(algnShifts);
185 fAlignCor = new Double_t**[nStat];
186 for (Int_t iStat = 0; iStat < nStat; iStat++) {
187 const Int_t nParams = 3;
188 Int_t nModul = StationSet->GetSiliconStation(iStat)->GetNModules();
189 fAlignCor[iStat] = new Double_t*[nModul];
190 for (Int_t iMod = 0; iMod < nModul; iMod++) {
191 fAlignCor[iStat][iMod] = new Double_t[nParams];
192 for (Int_t iPar = 0; iPar < nParams; iPar++) {
193 fAlignCor[iStat][iMod][iPar] = 0.0;
194 }
195 }
196 }
197
198 for (size_t i = 0; i < algnShifts.size(); ++i) {
199 // cout << ((AlignmentValue*)algnShifts[i])->value[0] << endl;
200 // cout << ((AlignmentValue*)algnShifts[i])->value[1] << endl;
201 // cout << ((AlignmentValue*)algnShifts[i])->value[2] << endl;
202 Int_t st = ((AlignmentValue*)algnShifts[i])->station;
203 Int_t mod = ((AlignmentValue*)algnShifts[i])->module;
204 fAlignCor[st][mod][0] = ((AlignmentValue*)algnShifts[i])->value[0];
205 fAlignCor[st][mod][1] = ((AlignmentValue*)algnShifts[i])->value[1];
206 fAlignCor[st][mod][2] = ((AlignmentValue*)algnShifts[i])->value[2];
207 }
208 }
209}
210
211BmnSiliconHitMaker::~BmnSiliconHitMaker()
212{
213
214 if (fIsExp) {
215 for (Int_t iStat = 0; iStat < StationSet->GetNStations(); iStat++) {
216 for (Int_t iMod = 0; iMod < StationSet->GetSiliconStation(iStat)->GetNModules(); iMod++)
217 delete[] fAlignCor[iStat][iMod];
218 delete[] fAlignCor[iStat];
219 }
220 delete[] fAlignCor;
221 }
222
223 if (StationSet) {
224 delete StationSet;
225 }
226
227 if (TransfSet) {
228 delete TransfSet;
229 }
230
231 if (fBmnSiliconHitsArray) {
232 fBmnSiliconHitsArray->Delete();
233 delete fBmnSiliconHitsArray;
234 };
235 if (fBmnSiliconUpperClustersArray) {
236 fBmnSiliconUpperClustersArray->Delete();
237 delete fBmnSiliconUpperClustersArray;
238 };
239 if (fBmnSiliconLowerClustersArray) {
240 fBmnSiliconLowerClustersArray->Delete();
241 delete fBmnSiliconLowerClustersArray;
242 }
243}
244
245InitStatus BmnSiliconHitMaker::Init()
246{
247
248 if (fVerbose > 1)
249 cout << "=================== BmnSiliconHitMaker::Init() started ================" << endl;
250
251 createSiliconDetector();
252
253 FairRootManager* ioman = FairRootManager::Instance();
254
255 fBmnSiliconDigitsArray = (TClonesArray*)ioman->GetObject(fInputDigitsBranchName);
256 if (!fBmnSiliconDigitsArray) {
257 cout << "BmnSiliconHitMaker::Init(): branch " << fInputDigitsBranchName
258 << " not found! Task will be deactivated" << endl;
259 SetActive(kFALSE);
260 return kERROR;
261 }
262 fBmnSiliconDigitMatchesArray = (TClonesArray*)ioman->GetObject(fInputDigitMatchesBranchName);
263
264 if (fVerbose > 1) {
265 if (fBmnSiliconDigitMatchesArray)
266 cout << " Strip matching information exists!\n";
267 else
268 cout << " Strip matching information doesn`t exist!\n";
269 }
270
271 fBmnSiliconHitsArray = new TClonesArray(fOutputHitsBranchName);
272 ioman->Register(fOutputHitsBranchName, "SILICON", fBmnSiliconHitsArray, kTRUE);
273 fBmnSiliconUpperClustersArray = new TClonesArray("StripCluster");
274 ioman->Register("BmnSiliconUpperCluster", "SILICON", fBmnSiliconUpperClustersArray, kTRUE);
275 fBmnSiliconLowerClustersArray = new TClonesArray("StripCluster");
276 ioman->Register("BmnSiliconLowerCluster", "SILICON", fBmnSiliconLowerClustersArray, kTRUE);
277
278 if (fVerbose > 1)
279 cout << "=================== BmnSiliconHitMaker::Init() finished ===============" << endl;
280
281 return kSUCCESS;
282}
283
284void BmnSiliconHitMaker::Exec(Option_t* opt)
285{
286
287 TStopwatch sw;
288 sw.Start();
289
290 if (!IsActive())
291 return;
292
293 fBmnSiliconHitsArray->Delete();
294 fBmnSiliconUpperClustersArray->Delete();
295 fBmnSiliconLowerClustersArray->Delete();
296
297 BmnSiliconLayer::SetLowerUniqueID(0);
298 BmnSiliconLayer::SetUpperUniqueID(0);
299
300 if (fVerbose > 1)
301 cout << "=================== BmnSiliconHitMaker::Exec() started ================" << endl;
302 if (fVerbose > 1)
303 cout << " BmnSiliconHitMaker::Exec(), Number of BmnSiliconDigits = " << fBmnSiliconDigitsArray->GetEntriesFast()
304 << "\n";
305
306 ProcessDigits();
307
308 if (fVerbose > 1)
309 cout << "=================== BmnSiliconHitMaker::Exec() finished ===============" << endl;
310
311 sw.Stop();
312 workTime += sw.RealTime();
313}
314
315void BmnSiliconHitMaker::ProcessDigits()
316{
317
318 // FairMCPoint* MCPoint;
319 BmnSiliconDigit* digit;
320 BmnMatch* strip_match; // MC-information for a strip
321
322 BmnSiliconStation* station;
323 BmnSiliconModule* module;
324
325 // Loading digits ---------------------------------------------------------------
326 Int_t AddedDigits = 0;
327 Int_t AddedStripDigitMatches = 0;
328
329 for (Int_t idigit = 0; idigit < fBmnSiliconDigitsArray->GetEntriesFast(); idigit++) {
330 digit = (BmnSiliconDigit*)fBmnSiliconDigitsArray->At(idigit);
331 if (!digit->IsGoodDigit())
332 continue;
333
334 if (digit->GetStripSignal() < fSignalLow || digit->GetStripSignal() > fSignalUp)
335 continue;
336
337 station = StationSet->GetSiliconStation(digit->GetStation());
338 module = station->GetModule(digit->GetModule());
339 if (!module)
340 return;
341
342 if (module->SetStripSignalInLayerByZoneId(digit->GetStripLayer(), digit->GetStripNumber(),
343 digit->GetStripSignal()))
344 AddedDigits++;
345
346 // Add a MC-match to the current strip if this MC-match array exists
347 if (fBmnSiliconDigitMatchesArray) {
348 strip_match = (BmnMatch*)fBmnSiliconDigitMatchesArray->At(idigit);
349 if (module->SetStripMatchInLayerByZoneId(digit->GetStripLayer(), digit->GetStripNumber(), *strip_match))
350 AddedStripDigitMatches++;
351 }
352
353 // Add a digit number match to the current strip
354 BmnMatch stripDigitNumberMatch; // digit number information for the current strip
355 stripDigitNumberMatch.AddLink(1.0, idigit);
356 module->SetStripDigitNumberMatchInLayerByZoneId(digit->GetStripLayer(), digit->GetStripNumber(),
357 stripDigitNumberMatch);
358 }
359
360 if (fVerbose > 1)
361 cout << " Processed strip digits : " << AddedDigits << "\n";
362 if (fVerbose > 1 && fBmnSiliconDigitMatchesArray)
363 cout << " Added strip digit MC-matches : " << AddedStripDigitMatches << "\n";
364 //------------------------------------------------------------------------------
365
366 // Processing digits
367 StationSet->ProcessPointsInDetector();
368
369 Int_t NCalculatedPoints = StationSet->CountNProcessedPointsInDetector();
370 if (fVerbose > 1)
371 cout << " Calculated points : " << NCalculatedPoints << "\n";
372 if (fVerbose == 1)
373 cout << "BmnSiliconHitMaker: " << NCalculatedPoints << " hits\n";
374
375 Int_t clear_matched_points_cnt = 0; // points with the only one match-indexes
376
377 map<Int_t, StripCluster> UniqueUpperClusters;
378 map<Int_t, StripCluster> UniqueLowerClusters;
379
380 for (Int_t iStation = 0; iStation < StationSet->GetNStations(); ++iStation) {
381 station = StationSet->GetSiliconStation(iStation);
382
383 for (Int_t iModule = 0; iModule < station->GetNModules(); ++iModule) {
384 module = station->GetModule(iModule);
385
386 for (Int_t iLayer = 0; iLayer < module->GetNStripLayers(); ++iLayer) {
387 BmnSiliconLayer layer = module->GetStripLayer(iLayer);
388
389 for (size_t iCluster = 0; iCluster < layer.GetStripClusters().size(); ++iCluster) {
390 StripCluster cluster = layer.GetStripClusters()[iCluster];
391 cluster.SetModule(iModule);
392 cluster.SetStation(iStation);
393
394 // Saving all strip clusters into branches (each cluster have its own unique id in the current
395 // event)
396 if (cluster.GetType() == UpperStripLayer) {
397 new ((*fBmnSiliconUpperClustersArray)[fBmnSiliconUpperClustersArray->GetEntriesFast()])
398 StripCluster(cluster);
399 } else {
400 new ((*fBmnSiliconLowerClustersArray)[fBmnSiliconLowerClustersArray->GetEntriesFast()])
401 StripCluster(cluster);
402 }
403 }
404 }
405
406 Int_t NIntersectionPointsInModule = module->GetNIntersectionPoints();
407
408 // Int_t NPseudoIntersectionsInModule = module->GetNPseudoIntersections();
409 Int_t NPseudoIntersectionsInModule = 0; // !!! IGNORE PSEUDO INTERSECTIONS !!!
410
411 for (Int_t iPoint = 0; iPoint < NIntersectionPointsInModule + NPseudoIntersectionsInModule; ++iPoint) {
412
413 // Double_t threshold = 2000000;// ADC
414
415 // Double_t sigL;
416 // Double_t sigU;
417 Double_t x, y, z;
418 Double_t x_err, y_err, z_err;
419
420 if (iPoint < NIntersectionPointsInModule) {
421 // sigL = module->GetIntersectionPoint_LowerLayerSripTotalSignal(iPoint);
422 // sigU = module->GetIntersectionPoint_UpperLayerSripTotalSignal(iPoint);
423 x = module->GetIntersectionPointX(iPoint);
424 y = module->GetIntersectionPointY(iPoint);
425 z = module->GetZPositionRegistered();
426 x_err = module->GetIntersectionPointXError(iPoint);
427 y_err = module->GetIntersectionPointYError(iPoint);
428 z_err = 0.0;
429 } else {
430 // sigL =
431 // module->GetPseudoIntersection_LowerLayerSripTotalSignal(iPoint-NIntersectionPointsInModule); sigU
432 // = module->GetPseudoIntersection_UpperLayerSripTotalSignal(iPoint-NIntersectionPointsInModule);
433 x = module->GetPseudoIntersectionX(iPoint - NIntersectionPointsInModule);
434 y = module->GetPseudoIntersectionY(iPoint - NIntersectionPointsInModule);
435 z = module->GetZPositionRegistered();
436 x_err = module->GetPseudoIntersectionXError(iPoint - NIntersectionPointsInModule);
437 y_err = module->GetPseudoIntersectionYError(iPoint - NIntersectionPointsInModule);
438 z_err = 0.0;
439 }
440
441 z += fIsExp ? fAlignCor[iStation][iModule][2] : 0.; // alignment shift
442
443 // Transform hit coordinates from local coordinate system to global
444 if (TransfSet) {
445 Plane3D::Point glob_point = TransfSet->ApplyTransforms(Plane3D::Point(-x, y, z), iStation, iModule);
446 x = -glob_point.X();
447 y = glob_point.Y();
448 z = glob_point.Z();
449 }
450
451 Int_t RefMCIndex = -1;
452
453 // MC-matching for the current hit (define RefMCIndex)) ---------
454 BmnMatch mc_match_hit;
455 if (iPoint < NIntersectionPointsInModule) {
456 mc_match_hit = module->GetIntersectionPointMatch(iPoint);
457 } else {
458 mc_match_hit = module->GetPseudoIntersectionMatch(iPoint - NIntersectionPointsInModule);
459 }
460
461 Int_t most_probably_index = -1;
462 Double_t max_weight = 0;
463
464 Int_t n_links = mc_match_hit.GetNofLinks();
465 if (n_links == 1)
466 clear_matched_points_cnt++;
467 for (Int_t ilink = 0; ilink < n_links; ilink++) {
468 Int_t index = mc_match_hit.GetLink(ilink).GetIndex();
469 Double_t weight = mc_match_hit.GetLink(ilink).GetWeight();
470 if (weight > max_weight) {
471 max_weight = weight;
472 most_probably_index = index;
473 }
474 }
475
476 RefMCIndex = most_probably_index;
477 //--------------------------------------------------------------
478
479 // Add hit ------------------------------------------------------
480 x *= -1; // invert to global X
481 Double_t deltaX = fIsExp ? fAlignCor[iStation][iModule][0] : 0.;
482 Double_t deltaY = fIsExp ? fAlignCor[iStation][iModule][1] : 0.;
483
484 x += deltaX;
485 y += deltaY;
486
487 new ((*fBmnSiliconHitsArray)[fBmnSiliconHitsArray->GetEntriesFast()])
488 BmnSiliconHit(0, TVector3(x, y, z), TVector3(x_err, y_err, z_err), RefMCIndex);
489
490 BmnSiliconHit* hit =
491 (BmnSiliconHit*)fBmnSiliconHitsArray->At(fBmnSiliconHitsArray->GetEntriesFast() - 1);
492 hit->SetStation(iStation);
493 hit->SetModule(iModule);
494 hit->SetIndex(fBmnSiliconHitsArray->GetEntriesFast() - 1);
495
496 if (iPoint < NIntersectionPointsInModule) {
497 hit->SetDigitNumberMatch(
498 module->GetIntersectionPointDigitNumberMatch(iPoint)); // digit number match for the hit
499 hit->SetUpperClusterIndex(module->GetUpperCluster(iPoint).GetUniqueID());
500 hit->SetLowerClusterIndex(module->GetLowerCluster(iPoint).GetUniqueID());
501 hit->SetPseudo(false);
502 } else {
503 hit->SetDigitNumberMatch(module->GetPseudoIntersectionDigitNumberMatch(
504 iPoint - NIntersectionPointsInModule)); // digit number match for the hit
505 hit->SetUpperClusterIndex(
506 module->GetUpperCluster_PseudoIntersections(iPoint - NIntersectionPointsInModule)
507 .GetUniqueID());
508 hit->SetLowerClusterIndex(
509 module->GetLowerCluster_PseudoIntersections(iPoint - NIntersectionPointsInModule)
510 .GetUniqueID());
511 hit->SetPseudo(true);
512 }
513 //--------------------------------------------------------------
514 /*
515 StripCluster ucls = module->GetUpperCluster(iPoint);
516 StripCluster lcls = module->GetLowerCluster(iPoint);
517 ucls.SetModule(iModule);
518 lcls.SetModule(iModule);
519 ucls.SetStation(iStation);
520 lcls.SetStation(iStation);
521
522 Int_t upperClusterIndex = ucls.GetUniqueID();
523 Int_t lowerClusterIndex = lcls.GetUniqueID();
524
525 //Filtration of duplicate clusters in the current module for the same type of strips
526 Bool_t duplicateUpperCluster = false;
527 for(auto it : UniqueUpperClusters) {
528 StripCluster cluster = it.second;
529
530 if(cluster.GetClusterSize() != ucls.GetClusterSize()) continue;
531
532 Int_t identic_cnt = 0;
533 for(Int_t istrip = 0; istrip < cluster.GetClusterSize(); ++istrip) {
534 Int_t strip_num = cluster.Strips.at(istrip);
535 Double_t strip_signal = cluster.Signals.at(istrip);
536
537 if( (strip_num == ucls.Strips.at(istrip)) && (strip_signal ==
538 ucls.Signals.at(istrip)) ) identic_cnt++; else break;
539 }
540 if(identic_cnt == ucls.GetClusterSize()) {
541 upperClusterIndex = cluster.GetUniqueID();
542 duplicateUpperCluster = true;
543 break;
544 }
545 }
546
547 Bool_t duplicateLowerCluster = false;
548 for(auto it : UniqueLowerClusters) {
549 StripCluster cluster = it.second;
550
551 if(cluster.GetClusterSize() != lcls.GetClusterSize()) continue;
552
553 Int_t identic_cnt = 0;
554 for(Int_t istrip = 0; istrip < cluster.GetClusterSize(); ++istrip) {
555 Int_t strip_num = cluster.Strips.at(istrip);
556 Double_t strip_signal = cluster.Signals.at(istrip);
557
558 if( (strip_num == lcls.Strips.at(istrip)) && (strip_signal ==
559 lcls.Signals.at(istrip)) ) identic_cnt++; else break;
560 }
561 if(identic_cnt == lcls.GetClusterSize()) {
562 lowerClusterIndex = cluster.GetUniqueID();
563 duplicateLowerCluster = true;
564 break;
565 }
566 }
567 //--------------------------------------------------------------
568
569 if(!duplicateUpperCluster) UniqueUpperClusters[upperClusterIndex] = ucls;
570 if(!duplicateLowerCluster) UniqueLowerClusters[lowerClusterIndex] = lcls;
571 hit->SetUpperClusterIndex(upperClusterIndex);
572 hit->SetLowerClusterIndex(lowerClusterIndex);
573 */
574 if (fHitMatching) {
575 BmnMatch digiMatch;
576 Int_t idx0, idx1;
577 if (iPoint < NIntersectionPointsInModule) {
578 digiMatch = module->GetIntersectionPointDigitNumberMatch(iPoint);
579 idx0 = digiMatch.GetLink(0).GetIndex();
580 idx1 = digiMatch.GetLink(1).GetIndex();
581 } else {
582 digiMatch = module->GetPseudoIntersectionDigitNumberMatch(iPoint - NIntersectionPointsInModule);
583 idx0 = digiMatch.GetLink(0).GetIndex();
584 idx1 = digiMatch.GetLink(0).GetIndex();
585 }
586 BmnMatch* digiMatch0 = (BmnMatch*)fBmnSiliconDigitMatchesArray->At(idx0);
587 BmnMatch* digiMatch1 = (BmnMatch*)fBmnSiliconDigitMatchesArray->At(idx1);
588
589 Bool_t hitOk = kFALSE;
590 for (Int_t ilink = 0; ilink < digiMatch0->GetNofLinks(); ilink++) {
591 Int_t iindex = digiMatch0->GetLink(ilink).GetIndex();
592 for (Int_t jlink = 0; jlink < digiMatch1->GetNofLinks(); jlink++) {
593 Int_t jindex = digiMatch1->GetLink(jlink).GetIndex();
594 if (iindex == jindex) {
595 hitOk = kTRUE;
596 break;
597 }
598 }
599 if (hitOk)
600 break;
601 }
602
603 hit->SetType(hitOk);
604 if (!hitOk)
605 hit->SetRefIndex(-1);
606
607 // hit MC-matching ----------------------------------------------
608 FairRootManager::Instance()->SetUseFairLinks(kTRUE);
609 BmnMatch hitMatch;
610 if (iPoint < NIntersectionPointsInModule) {
611 hitMatch = module->GetIntersectionPointMatch(iPoint);
612 } else {
613 hitMatch = module->GetPseudoIntersectionMatch(iPoint - NIntersectionPointsInModule);
614 }
615
616 for (BmnLink lnk : hitMatch.GetLinks())
617 hit->AddLink(FairLink(-1, lnk.GetIndex(), lnk.GetWeight()));
618 FairRootManager::Instance()->SetUseFairLinks(kFALSE);
619 }
620 //--------------------------------------------------------------
621 }
622 }
623 }
624
625 /*
626 for (auto it : UniqueUpperClusters) {
627 for (Int_t i = 0; i < fBmnSiliconHitsArray->GetEntriesFast(); i++) {
628 BmnSiliconHit* hit = (BmnSiliconHit*)fBmnSiliconHitsArray->At(i);
629 if (hit->GetUpperClusterIndex() != it.first) continue;
630 hit->SetUpperClusterIndex(fBmnSiliconUpperClustersArray->GetEntriesFast());
631 }
632 it.second.SetUniqueID(fBmnSiliconUpperClustersArray->GetEntriesFast());
633 new ((*fBmnSiliconUpperClustersArray)[fBmnSiliconUpperClustersArray->GetEntriesFast()])
634 StripCluster(it.second);
635 }
636 for (auto it : UniqueLowerClusters) {
637 for (Int_t i = 0; i < fBmnSiliconHitsArray->GetEntriesFast(); i++) {
638 BmnSiliconHit* hit = (BmnSiliconHit*)fBmnSiliconHitsArray->At(i);
639 if (hit->GetLowerClusterIndex() != it.first) continue;
640 hit->SetLowerClusterIndex(fBmnSiliconLowerClustersArray->GetEntriesFast());
641 }
642 it.second.SetUniqueID(fBmnSiliconLowerClustersArray->GetEntriesFast());
643 new ((*fBmnSiliconLowerClustersArray)[fBmnSiliconLowerClustersArray->GetEntriesFast()])
644 StripCluster(it.second);
645 }
646 */
647
648 if (fVerbose > 1)
649 cout << " N clear matches with MC-points = " << clear_matched_points_cnt << "\n";
650 //------------------------------------------------------------------------------
651 StationSet->Reset();
652}
653
654void BmnSiliconHitMaker::OnlineWrite(const std::unique_ptr<TTree>& resultTree)
655{
656 if (!IsActive())
657 return;
658
659 resultTree->Branch(fOutputHitsBranchName, &fBmnSiliconHitsArray);
660 resultTree->Branch("BmnSiliconUpperCluster", &fBmnSiliconUpperClustersArray);
661 resultTree->Branch("BmnSiliconLowerCluster", &fBmnSiliconLowerClustersArray);
662 resultTree->Fill();
663}
664
665void BmnSiliconHitMaker::Finish()
666{
667 printf("Work time of BmnSiliconHitMaker: %4.2f sec.\n", workTime);
668}
i
int i
Definition P4_F32vec4.h:22
UpperStripLayer
@ UpperStripLayer
Definition BmnStripData.h:13
BmnHit::SetModule
void SetModule(Int_t mod)
Definition BmnHit.h:73
BmnHit::SetType
void SetType(Int_t type)
Definition BmnHit.h:81
BmnHit::SetUpperClusterIndex
void SetUpperClusterIndex(Int_t idx)
Definition BmnHit.h:133
BmnHit::SetDigitNumberMatch
void SetDigitNumberMatch(BmnMatch match)
Definition BmnHit.h:121
BmnHit::SetIndex
void SetIndex(Int_t id)
Definition BmnHit.h:57
BmnHit::SetPseudo
void SetPseudo(Bool_t pseudo)
Definition BmnHit.h:141
BmnHit::SetStation
void SetStation(Short_t st)
Definition BmnHit.h:69
BmnHit::SetLowerClusterIndex
void SetLowerClusterIndex(Int_t idx)
Definition BmnHit.h:137
BmnMatch::GetLink
const BmnLink & GetLink(Int_t i) const
Definition BmnMatch.h:37
BmnMatch::AddLink
void AddLink(const BmnMatch &match)
Definition BmnMatch.cxx:43
BmnMatch::GetLinks
const vector< BmnLink > & GetLinks() const
Definition BmnMatch.h:38
BmnMatch::GetNofLinks
Int_t GetNofLinks() const
Definition BmnMatch.h:40
BmnSiliconHitMaker::Exec
virtual void Exec(Option_t *opt)
Definition BmnSiliconHitMaker.cxx:284
BmnSiliconHitMaker::OnlineWrite
virtual void OnlineWrite(const std::unique_ptr< TTree > &resultTree)
Write task resul to tree.
Definition BmnSiliconHitMaker.cxx:654
BmnSiliconHitMaker::Init
virtual InitStatus Init()
Definition BmnSiliconHitMaker.cxx:245
BmnSiliconLayer::SetLowerUniqueID
static void SetLowerUniqueID(Int_t id)
Definition BmnSiliconLayer.h:111
BmnSiliconLayer::GetStripClusters
vector< StripCluster > GetStripClusters()
Definition BmnSiliconLayer.h:110
BmnSiliconLayer::SetUpperUniqueID
static void SetUpperUniqueID(Int_t id)
Definition BmnSiliconLayer.h:112
BmnSiliconStationSet::GetSiliconStation
BmnSiliconStation * GetSiliconStation(Int_t station_num)
Definition BmnSiliconStationSet.cxx:89
BmnSiliconTransform::LoadFromXMLFile
Bool_t LoadFromXMLFile(TString xml_config_file)
Definition BmnSiliconTransform.cxx:95
BmnSiliconTransform::ApplyTransforms
Plane3D::Point ApplyTransforms(Plane3D::Point point, Int_t station, Int_t module)
Definition BmnSiliconTransform.cxx:225
BmnStripDigit::IsGoodDigit
Bool_t IsGoodDigit()
Definition BmnStripDigit.h:32
BmnStripDigit::GetStripNumber
Int_t GetStripNumber()
Definition BmnStripDigit.h:29
BmnStripDigit::GetStripLayer
Int_t GetStripLayer()
Definition BmnStripDigit.h:28
BmnStripDigit::GetStation
Int_t GetStation()
Definition BmnStripDigit.h:26
BmnStripDigit::GetStripSignal
Double_t GetStripSignal()
Definition BmnStripDigit.h:30
StripCluster::SetStation
void SetStation(Int_t st)
Definition BmnStripData.h:215
StripCluster::GetType
Int_t GetType()
Definition BmnStripData.h:211
StripCluster::SetModule
void SetModule(Int_t m)
Definition BmnStripData.h:217
UniDetectorParameter::GetValue
int GetValue(vector< UniValue * > &parameter_value)
get value of detector parameter presented by an array
Definition UniDetectorParameter.cxx:1107
UniDetectorParameter::GetDetectorParameter
static UniDetectorParameter * GetDetectorParameter(int value_id)
get detector parameter from the database
Definition UniDetectorParameter.cxx:163