BmnRoot
Loading...
Searching...
No Matches
BmnNdetClusterFinder.cxx
Go to the documentation of this file.
1
7#include "BmnNdetClusterFinder.h"
8
9#include "TMath.h"
10
11#include <iostream>
12
13using namespace std;
14
15BmnNdetClusterFinder::BmnNdetClusterFinder(const Int_t kNColumns, const Int_t kNRows, const Int_t kNLayers)
16 : TObject()
17 , fClusters()
18 , fIndexXYZ(kNColumns, vector<vector<Int_t>>(kNRows, vector<Int_t>(kNLayers, -1)))
19{
20 fNCells = 0;
21 fNClusters = 0;
22}
23
24BmnNdetClusterFinder::~BmnNdetClusterFinder()
25{
26 for (Int_t c = 0; c < GetNColumns(); ++c) {
27 for (Int_t r = 0; r < GetNRows(); ++r) {
28 fIndexXYZ[c][r].clear();
29 }
30 fIndexXYZ[c].clear();
31 }
32 fIndexXYZ.clear();
33}
34
35void BmnNdetClusterFinder::Fill(Int_t cell,
36 Int_t column,
37 Int_t row,
38 Int_t layer,
39 Double_t x,
40 Double_t y,
41 Double_t z,
42 Double_t t,
43 Double_t e_dep)
44{
45 Double_t e2 = e_dep;
46 Double_t t2 = t;
47 if (t2 <= 0.)
48 return;
49 if (e2 <= 0.)
50 return;
51 Int_t index = fIndexXYZ[column][row][layer];
52 if (index >= 0) {
53 Double_t t1 = fTime[index];
54 Double_t e1 = fEdep[index];
55 t2 = (t2 > t1 ? t1 : t2);
56 e2 += e1;
57 }
58 if (column < 0 || column >= GetNColumns())
59 return;
60 if (row < 0 || row >= GetNRows())
61 return;
62 if (layer < 0 || layer >= GetNLayers())
63 return;
64 fCell.push_back(cell);
65 fColumn.push_back(column);
66 fRow.push_back(row);
67 fLayer.push_back(layer);
68 fXcoord.push_back(x);
69 fYcoord.push_back(y);
70 fZcoord.push_back(z);
71 fTime.push_back(t2);
72 fEdep.push_back(e2);
73
74 Double_t distance = TMath::Sqrt(TMath::Sq(x) + TMath::Sq(y) + TMath::Sq(z));
75 Double_t beta = (distance * 1.E-2) / (t2 * 1.E-9) / TMath::C();
76 fBeta.push_back(beta);
77 fIndexTime.push_back(fNCells);
78 fIndexEdep.push_back(fNCells);
79 fIndexBeta.push_back(fNCells);
80 fCellStatus.push_back(kIsNotUsed);
81
82 fIndexXYZ[column][row][layer] = fNCells;
83
84 Int_t temp;
85 // Sort ascending Time
86 for (Int_t i = fNCells; i > 0; --i) {
87 if (fTime[fIndexTime[i - 1]] < fTime[fIndexTime[i]])
88 break;
89 temp = fIndexTime[i - 1];
90 fIndexTime[i - 1] = fIndexTime[i];
91 fIndexTime[i] = temp;
92 }
93 for (Int_t i = fNCells; i > 0; --i) {
94 if (fEdep[fIndexEdep[i - 1]] > fEdep[fIndexEdep[i]])
95 break;
96 temp = fIndexEdep[i - 1];
97 fIndexEdep[i - 1] = fIndexEdep[i];
98 fIndexEdep[i] = temp;
99 }
100 for (Int_t i = fNCells; i > 0; --i) {
101 if (fBeta[fIndexBeta[i - 1]] > fBeta[fIndexBeta[i]])
102 break;
103 temp = fIndexBeta[i - 1];
104 fIndexBeta[i - 1] = fIndexBeta[i];
105 fIndexBeta[i] = temp;
106 }
107 ++fNCells;
108}
109
110Int_t BmnNdetClusterFinder::GetNClusters()
111{
112 return fNClusters;
113}
114
115Int_t BmnNdetClusterFinder::GetNCells()
116{
117 return fNCells;
118}
119
120Int_t BmnNdetClusterFinder::GetNCells(Double_t threshold)
121{
122 Int_t nCells = 0.;
123 for (Int_t iCell = 0; iCell < fNCells; ++iCell) {
124 if (fEdep[iCell] > threshold)
125 ++nCells;
126 }
127 return nCells;
128}
129
130Int_t BmnNdetClusterFinder::GetNColumns()
131{
132 return fIndexXYZ.size();
133}
134
135Int_t BmnNdetClusterFinder::GetNRows()
136{
137 return fIndexXYZ[0].size();
138}
139
140Int_t BmnNdetClusterFinder::GetNLayers()
141{
142 return fIndexXYZ[0][0].size();
143}
144
145vector<Int_t> BmnNdetClusterFinder::GetCells()
146{
147 return fCell;
148}
149vector<Int_t> BmnNdetClusterFinder::GetColumns()
150{
151 return fColumn;
152}
153vector<Int_t> BmnNdetClusterFinder::GetRows()
154{
155 return fRow;
156}
157vector<Int_t> BmnNdetClusterFinder::GetLayers()
158{
159 return fLayer;
160}
161vector<Double_t> BmnNdetClusterFinder::GetX()
162{
163 return fXcoord;
164}
165vector<Double_t> BmnNdetClusterFinder::GetY()
166{
167 return fYcoord;
168}
169vector<Double_t> BmnNdetClusterFinder::GetZ()
170{
171 return fZcoord;
172}
173vector<Double_t> BmnNdetClusterFinder::GetTime()
174{
175 return fTime;
176}
177vector<Double_t> BmnNdetClusterFinder::GetEdep()
178{
179 return fEdep;
180}
181vector<Double_t> BmnNdetClusterFinder::GetBeta()
182{
183 return fBeta;
184}
185
186BmnNdetCluster* BmnNdetClusterFinder::GetCluster(Int_t i_cluster)
187{
188 if (i_cluster < 0 || i_cluster >= fNClusters)
189 return nullptr;
190 return fClusters[i_cluster];
191}
192
193Int_t BmnNdetClusterFinder::GetClusterStatus(Int_t i_cluster)
194{
195 if (i_cluster >= fNClusters || i_cluster < 0)
196 return kNoCluster;
197 return fClusterStatus[i_cluster];
198}
199
200void BmnNdetClusterFinder::SetClusterStatus(Int_t i_cluster, eClusterStatus status)
201{
202 if (i_cluster >= fNClusters || i_cluster < 0)
203 return;
204
205 fClusterStatus[i_cluster] = status;
206}
207
208BmnNdetCluster* BmnNdetClusterFinder::AddCluster()
209{
210 BmnNdetCluster* cl = new BmnNdetCluster();
211 fClusters.push_back(cl);
212 fClusterStatus.push_back(kUnchecked);
213 ++fNClusters;
214 return cl;
215}
216
217void BmnNdetClusterFinder::AddCellToCluster(Int_t i_cell, Int_t i_cluster = -1)
218{
219 if (i_cluster == -1)
220 i_cluster = fNClusters - 1;
221 if (i_cell < 0 || i_cell >= fNCells)
222 return;
223 if (i_cluster < 0 || i_cluster >= fNClusters)
224 return;
225 BmnNdetCluster* cl = GetCluster(i_cluster);
226 Int_t i = fCell[i_cell];
227 Int_t c = fColumn[i_cell];
228 Int_t r = fRow[i_cell];
229 Int_t l = fLayer[i_cell];
230 Double_t x = fXcoord[i_cell];
231 Double_t y = fYcoord[i_cell];
232 Double_t z = fZcoord[i_cell];
233 Double_t t = fTime[i_cell];
234 Double_t e = fEdep[i_cell];
235 Double_t b = fBeta[i_cell];
236 cl->AddCell(i, c, r, l, x, y, z, t, e, b);
237 fCellStatus[i_cell] = kIsUsed;
238 --fNUncheckedCellsInEvent;
239}
240
241Int_t BmnNdetClusterFinder::FindClusters1()
242{
243 fNUncheckedCellsInEvent = fNCells;
244 while (fNUncheckedCellsInEvent > 0) {
245 // Search for the first unchecked cell
246 Int_t iCell1 = -1;
247 for (Int_t iTime = 0; iTime < fNCells; ++iTime) {
248 Int_t iCell = fIndexTime[iTime];
249 if (fCellStatus[iCell] == kIsNotUsed) {
250 iCell1 = iCell;
251 break;
252 }
253 }
254 if (iCell1 < 0)
255 return fNClusters;
256 else {
257 // Start a new cluster
258 BmnNdetCluster* cl = AddCluster();
259 AddCellToCluster(iCell1);
260 // Chain of checkouts of all neighbours and their neighbours
261 Int_t iCellInCluster = 0;
262 while (iCellInCluster < cl->GetNCells()) {
263 // Loop over all neighbours
264 for (Int_t dIx : {-1, 0, 1}) {
265 Int_t iX = cl->GetColumns()[iCellInCluster] + dIx;
266 if (iX < 0 || iX >= GetNColumns())
267 continue;
268 for (Int_t dIy : {-1, 0, 1}) {
269 Int_t iY = cl->GetRows()[iCellInCluster] + dIy;
270 if (iY < 0 || iY >= GetNRows())
271 continue;
272 for (Int_t dIz : {-1, 0, 1}) {
273 Int_t iZ = cl->GetLayers()[iCellInCluster] + dIz;
274 if (iZ < 0 || iZ >= GetNLayers())
275 continue;
276
277 // Check neighbour
278 Int_t iNeighbour = fIndexXYZ[iX][iY][iZ];
279 if (iNeighbour < 0)
280 continue;
281 if (fCellStatus[iNeighbour] == kIsNotUsed)
282 AddCellToCluster(iNeighbour);
283 }
284 }
285 }
286 ++iCellInCluster;
287 }
288 }
289 }
290 return fNClusters;
291}
292
293Int_t BmnNdetClusterFinder::FindClusters1_2()
294{
295 fNUncheckedCellsInEvent = fNCells;
296 while (fNUncheckedCellsInEvent > 0) {
297 // Search for the first unchecked cell
298 Int_t iCell1 = -1;
299 for (Int_t iTime = 0; iTime < fNCells; ++iTime) {
300 Int_t iCell = fIndexTime[iTime];
301 if (fCellStatus[iCell] == kIsNotUsed) {
302 iCell1 = iCell;
303 break;
304 }
305 }
306 if (iCell1 < 0)
307 return fNClusters;
308 else {
309 // Start a new cluster
310 BmnNdetCluster* cl = AddCluster();
311 AddCellToCluster(iCell1);
312 Bool_t checkAllLayers = kFALSE;
313 if (cl->GetLayers()[0] == 0)
314 checkAllLayers = kTRUE;
315 // Chain of checkouts of all neighbours and their neighbours
316 Int_t iCellInCluster = 0;
317 while (iCellInCluster < cl->GetNCells()) {
318 // Loop over all neighbours
319 for (Int_t dIx : {-1, 0, 1}) {
320 Int_t iX = cl->GetColumns()[iCellInCluster] + dIx;
321 if (iX < 0 || iX >= GetNColumns())
322 continue;
323 for (Int_t dIy : {-1, 0, 1}) {
324 Int_t iY = cl->GetRows()[iCellInCluster] + dIy;
325 if (iY < 0 || iY >= GetNRows())
326 continue;
327 if (checkAllLayers) {
328 for (Int_t iZ = 0; iZ < GetNLayers(); ++iZ) {
329 // Check neighbour
330 Int_t iNeighbour = fIndexXYZ[iX][iY][iZ];
331 if (iNeighbour < 0)
332 continue;
333 if (fCellStatus[iNeighbour] == kIsNotUsed)
334 AddCellToCluster(iNeighbour);
335 }
336 } else {
337 for (Int_t dIz : {-1, 0, 1}) {
338 Int_t iZ = cl->GetLayers()[iCellInCluster] + dIz;
339 if (iZ < 0 || iZ >= GetNLayers())
340 continue;
341
342 // Check neighbour
343 Int_t iNeighbour = fIndexXYZ[iX][iY][iZ];
344 if (iNeighbour < 0)
345 continue;
346 if (fCellStatus[iNeighbour] == kIsNotUsed)
347 AddCellToCluster(iNeighbour);
348 }
349 }
350 }
351 }
352 ++iCellInCluster;
353 }
354 }
355 }
356 return fNClusters;
357}
358
359Int_t BmnNdetClusterFinder::FindClusters2()
360{
361 fNUncheckedCellsInEvent = fNCells;
362 while (fNUncheckedCellsInEvent > 0) {
363 // Search for the first unchecked cell
364 Int_t iCell1 = -1;
365 for (Int_t iBeta = 0; iBeta < fNCells; ++iBeta) {
366 Int_t iCell = fIndexBeta[iBeta];
367 if (fCellStatus[iCell] == kIsNotUsed) {
368 iCell1 = iCell;
369 break;
370 }
371 }
372 if (iCell1 < 0)
373 return fNClusters;
374 else {
375 // Start a new cluster
376 BmnNdetCluster* cl = AddCluster();
377 AddCellToCluster(iCell1);
378 // Chain of checkouts of all neighbours and their neighbours
379 Int_t iCellInCluster = 0;
380 while (iCellInCluster < cl->GetNCells()) {
381 // Loop over all neighbours
382 for (Int_t dIx : {-1, 0, 1}) {
383 Int_t iX = cl->GetColumns()[iCellInCluster] + dIx;
384 if (iX < 0 || iX >= GetNColumns())
385 continue;
386 for (Int_t dIy : {-1, 0, 1}) {
387 Int_t iY = cl->GetRows()[iCellInCluster] + dIy;
388 if (iY < 0 || iY >= GetNRows())
389 continue;
390 for (Int_t dIz : {-1, 0, 1}) {
391 Int_t iZ = cl->GetLayers()[iCellInCluster] + dIz;
392 if (iZ < 0 || iZ >= GetNLayers())
393 continue;
394
395 // Check neighbour
396 Int_t iNeighbour = fIndexXYZ[iX][iY][iZ];
397 if (iNeighbour < 0)
398 continue;
399 if (fCellStatus[iNeighbour] == kIsNotUsed)
400 AddCellToCluster(iNeighbour);
401 }
402 }
403 }
404 ++iCellInCluster;
405 }
406 }
407 }
408 return fNClusters;
409}
410
411Int_t BmnNdetClusterFinder::FindClusters2_2()
412{
413 fNUncheckedCellsInEvent = fNCells;
414 while (fNUncheckedCellsInEvent > 0) {
415 // Search for the first unchecked cell
416 Int_t iCell1 = -1;
417 for (Int_t iBeta = 0; iBeta < fNCells; ++iBeta) {
418 Int_t iCell = fIndexBeta[iBeta];
419 if (fCellStatus[iCell] == kIsNotUsed) {
420 iCell1 = iCell;
421 break;
422 }
423 }
424 if (iCell1 < 0)
425 return fNClusters;
426 else {
427 // Start a new cluster
428 BmnNdetCluster* cl = AddCluster();
429 AddCellToCluster(iCell1);
430 // Chain of checkouts of all neighbours and their neighbours
431 Int_t iCellInCluster = 0;
432 while (iCellInCluster < cl->GetNCells()) {
433 // Loop over all neighbours
434 for (Int_t dIx : {-1, 0, 1}) {
435 Int_t iX = cl->GetColumns()[iCellInCluster] + dIx;
436 if (iX < 0 || iX >= GetNColumns())
437 continue;
438 for (Int_t dIy : {-1, 0, 1}) {
439 Int_t iY = cl->GetRows()[iCellInCluster] + dIy;
440 if (iY < 0 || iY >= GetNRows())
441 continue;
442 for (Int_t dIz : {-1, 0, 1, 2}) {
443 Int_t iZ = cl->GetLayers()[iCellInCluster] + dIz;
444 if (iZ < 0 || iZ >= GetNLayers())
445 continue;
446
447 // Check neighbour
448 Int_t iNeighbour = fIndexXYZ[iX][iY][iZ];
449 if (iNeighbour < 0)
450 continue;
451 if (fCellStatus[iNeighbour] == kIsNotUsed)
452 AddCellToCluster(iNeighbour);
453 }
454 }
455 }
456 ++iCellInCluster;
457 }
458 }
459 }
460 return fNClusters;
461}
462
463Int_t BmnNdetClusterFinder::FindClusters3(Double_t cut_beta)
464{
465 fNUncheckedCellsInEvent = fNCells;
466 while (fNUncheckedCellsInEvent > 0) {
467 // Search for the first unchecked cell
468 Int_t iCell1 = -1;
469 for (Int_t iBeta = 0; iBeta < fNCells; ++iBeta) {
470 Int_t iCell = fIndexBeta[iBeta];
471 if (fCellStatus[iCell] == kIsNotUsed) {
472 iCell1 = iCell;
473 break;
474 }
475 }
476 if (iCell1 < 0)
477 return fNClusters;
478 else {
479 // Start a new cluster
480 BmnNdetCluster* cl = AddCluster();
481 AddCellToCluster(iCell1);
482 Double_t beta0 = fBeta[iCell1];
483 // Chain of checkouts of all neighbours and their neighbours
484 Int_t iCellInCluster = 0;
485 while (iCellInCluster < cl->GetNCells()) {
486 // Loop over all neighbours
487 for (Int_t dIx : {-1, 0, 1}) {
488 Int_t iX = cl->GetColumns()[iCellInCluster] + dIx;
489 if (iX < 0 || iX >= GetNColumns())
490 continue;
491 for (Int_t dIy : {-1, 0, 1}) {
492 Int_t iY = cl->GetRows()[iCellInCluster] + dIy;
493 if (iY < 0 || iY >= GetNRows())
494 continue;
495 for (Int_t dIz : {-1, 0, 1}) {
496 Int_t iZ = cl->GetLayers()[iCellInCluster] + dIz;
497 if (iZ < 0 || iZ >= GetNLayers())
498 continue;
499
500 // Check neighbour
501 Int_t iNeighbour = fIndexXYZ[iX][iY][iZ];
502 if (iNeighbour < 0)
503 continue;
504 if (fCellStatus[iNeighbour] == kIsUsed)
505 continue;
506 if (TMath::Abs(fBeta[iNeighbour] - beta0) > cut_beta)
507 continue;
508 AddCellToCluster(iNeighbour);
509 }
510 }
511 }
512 ++iCellInCluster;
513 }
514 }
515 }
516 return fNClusters;
517}
518
519Int_t BmnNdetClusterFinder::FindClusters4(Double_t cut_beta)
520{
521 fNUncheckedCellsInEvent = fNCells;
522 while (fNUncheckedCellsInEvent > 0) {
523 // Search for the first unchecked cell
524 Int_t iCell1 = -1;
525 for (Int_t iTime = 0; iTime < fNCells; ++iTime) {
526 Int_t iCell = fIndexTime[iTime];
527 if (fCellStatus[iCell] == kIsNotUsed) {
528 iCell1 = iCell;
529 break;
530 }
531 }
532 if (iCell1 < 0)
533 return fNClusters;
534 else {
535 // Start a new cluster
536 BmnNdetCluster* cl = AddCluster();
537 AddCellToCluster(iCell1);
538
539 Double_t x0 = cl->GetX()[0];
540 Double_t y0 = cl->GetY()[0];
541 Double_t z0 = cl->GetZ()[0];
542 Double_t t0 = cl->GetTime()[0];
543 Double_t beta0 = cl->BetaFirst();
544 // Chain of checkouts of all neighbours and their neighbours
545 Int_t iCellInCluster = 0;
546 while (iCellInCluster < cl->GetNCells()) {
547 Double_t x1 = cl->GetX()[iCellInCluster];
548 Double_t y1 = cl->GetY()[iCellInCluster];
549 Double_t z1 = cl->GetZ()[iCellInCluster];
550 Double_t t1 = cl->GetTime()[iCellInCluster];
551 // Loop over all neighbours
552 for (Int_t dIx : {-1, 0, 1}) {
553 Int_t iX = cl->GetColumns()[iCellInCluster] + dIx;
554 if (iX < 0 || iX >= GetNColumns())
555 continue;
556 for (Int_t dIy : {-1, 0, 1}) {
557 Int_t iY = cl->GetRows()[iCellInCluster] + dIy;
558 if (iY < 0 || iY >= GetNRows())
559 continue;
560 for (Int_t dIz : {-1, 0, 1}) {
561 Int_t iZ = cl->GetLayers()[iCellInCluster] + dIz;
562 if (iZ < 0 || iZ >= GetNLayers())
563 continue;
564
565 // Check neighbour
566 Int_t iNeighbour = fIndexXYZ[iX][iY][iZ];
567 if (iNeighbour < 0)
568 continue;
569 Double_t x2 = cl->GetX()[iNeighbour];
570 Double_t y2 = cl->GetY()[iNeighbour];
571 Double_t z2 = cl->GetZ()[iNeighbour];
572 Double_t t2 = cl->GetTime()[iNeighbour];
573 Double_t distance21 =
574 TMath::Sqrt(TMath::Sq(x2 - x1) + TMath::Sq(y2 - y1) + TMath::Sq(z2 - z1));
575 Double_t distance20 =
576 TMath::Sqrt(TMath::Sq(x2 - x0) + TMath::Sq(y2 - y0) + TMath::Sq(z2 - z0));
577 Double_t d_time20 = t2 - t0;
578 Double_t d_time21 = t2 - t1;
579 Double_t beta20 =
580 d_time20 <= 0 ? -1. : distance20 * 1.E-2 / (d_time20 * 1.E-9) / TMath::C();
581 Double_t beta21 =
582 d_time21 <= 0 ? -1. : distance21 * 1.E-2 / (d_time21 * 1.E-9) / TMath::C();
583 if (beta20 > 1. || beta20 < 0.)
584 continue;
585 if (beta21 > 1. || beta21 < 0.)
586 continue;
587 if (TMath::Abs(beta21 - beta0) > cut_beta)
588 continue;
589 if (fCellStatus[iNeighbour] == kIsUsed)
590 continue;
591 else
592 AddCellToCluster(iNeighbour);
593 }
594 }
595 }
596 ++iCellInCluster;
597 }
598 }
599 }
600 return fNClusters;
601}
602
603Int_t BmnNdetClusterFinder::MergeClusters(Int_t i_cluster_1, Int_t i_cluster_2)
604{
605 if (i_cluster_1 == i_cluster_2)
606 return fNClusters;
607 Int_t i1 = min(i_cluster_1, i_cluster_2);
608 Int_t i2 = max(i_cluster_1, i_cluster_2);
609
610 BmnNdetCluster* cl1 = GetCluster(i1);
611 BmnNdetCluster* cl2 = GetCluster(i2);
612
613 for (Int_t iCell2 = 0; iCell2 < cl2->GetNCells(); ++iCell2) {
614
615 cl1->AddCell(cl2->GetColumns()[iCell2], cl2->GetRows()[iCell2], cl2->GetLayers()[iCell2], cl2->GetX()[iCell2],
616 cl2->GetY()[iCell2], cl2->GetZ()[iCell2], cl2->GetTime()[iCell2], cl2->GetEdep()[iCell2],
617 cl2->GetBeta()[iCell2]);
618 }
619 fClusters.erase(fClusters.begin() + i2);
620 fClusterStatus.erase(fClusterStatus.begin() + i2);
621 --fNClusters;
622 return fNClusters;
623}
624
625Int_t BmnNdetClusterFinder::SelectNeutrons(Double_t cut_angle, Int_t cut_ncells, Bool_t check_border)
626{
627 Int_t nNeutrons = 0;
628 BmnNdetCluster* cl = nullptr;
629 for (Int_t iCluster = 0; iCluster < fNClusters; ++iCluster) {
630 // if(fClusterStatus[iCluster]!=kUnchecked)continue;
631 cl = GetCluster(iCluster);
632 if (cl->AngleToTarget() < cut_angle && cl->GetNCells() >= cut_ncells
633 && (!cl->StartsOnBorder(GetNColumns(), GetNRows(), GetNLayers()) && check_border))
634 {
635 fClusterStatus[iCluster] = kIsNeutron;
636 ++nNeutrons;
637 } else {
638 fClusterStatus[iCluster] = kNotNeutron;
639 }
640 }
641 cl = nullptr;
642 return nNeutrons;
643}
644
645Int_t BmnNdetClusterFinder::CalculateNumberOfNeutrons(Double_t weight_edep,
646 Double_t weight_ncells,
647 Double_t weight_nclusters,
648 Double_t threshold12)
649{
650 Double_t edep = 0.;
651 Int_t n_cells = 0;
652 Int_t n_clusters = 0;
653 BmnNdetCluster* cl = nullptr;
654 for (Int_t iCluster = 0; iCluster < fNClusters; ++iCluster) {
655 if (fClusterStatus[iCluster] != kIsNeutron)
656 continue;
657 cl = GetCluster(iCluster);
658 edep += cl->EnergyDeposited();
659 n_cells += Double_t(cl->GetNCells());
660 n_clusters += 1.;
661 }
662 cl = nullptr;
663 if (n_clusters < 1)
664 return 0;
665 Double_t separation_parameter = weight_edep * edep + weight_ncells * n_cells + weight_nclusters * n_clusters;
666 if (separation_parameter > threshold12)
667 return 2;
668 return 1;
669}
670
671void BmnNdetClusterFinder::Clear()
672{
673 for (Int_t i = 0; i < fNCells; ++i) {
674 fCell.pop_back();
675 fColumn.pop_back();
676 fRow.pop_back();
677 fLayer.pop_back();
678 fXcoord.pop_back();
679 fYcoord.pop_back();
680 fZcoord.pop_back();
681 fTime.pop_back();
682 fEdep.pop_back();
683 fBeta.pop_back();
684 fIndexTime.pop_back();
685 fIndexEdep.pop_back();
686 fIndexBeta.pop_back();
687 fCellStatus.pop_back();
688 }
689 fNCells = 0;
690
691 for (Int_t c = 0; c < GetNColumns(); ++c) {
692 for (Int_t r = 0; r < GetNRows(); ++r) {
693 for (Int_t l = 0; l < GetNLayers(); ++l) {
694 fIndexXYZ[c][r][l] = -1;
695 }
696 }
697 }
698
699 for (Int_t i = 0; i < fNClusters; ++i) {
700 fClusters.pop_back();
701 fClusterStatus.pop_back();
702 }
703 fNClusters = 0;
704}
705
706void BmnNdetClusterFinder::Print()
707{
708 cout << "----------------------------- BmnNdetClusterFinder -----------------------------" << endl;
709 for (Int_t i = 0; i < fNCells; ++i) {
710 cout << fCell[i] << " cell (" << fColumn[i] << ", " << fRow[i] << ", " << fLayer[i] << "): ";
711 cout << " i: " << fIndexXYZ[fColumn[i]][fRow[i]][fLayer[i]];
712 cout << "\tX: " << fXcoord[i] << "[cm], Y: " << fYcoord[i] << "[cm], Z: " << fZcoord[i] << "[cm]";
713 cout << "\tt: " << fTime[i] << "[ns], Edep: " << 1.E3 * fEdep[i] << "[MeV]" << endl;
714 }
715 cout << endl;
716
717 cout << "Sorted time" << endl;
718 for (Int_t i = 0; i < fNCells; ++i) {
719 cout << fIndexTime[i] << "(" << fTime[fIndexTime[i]] << "), ";
720 }
721 cout << endl;
722
723 cout << "Sorted edep" << endl;
724 for (Int_t i = 0; i < fNCells; ++i) {
725 cout << fIndexEdep[i] << "(" << fEdep[fIndexEdep[i]] << "), ";
726 }
727 cout << endl;
728
729 cout << "Sorted beta" << endl;
730 for (Int_t i = 0; i < fNCells; ++i) {
731 cout << fIndexBeta[i] << "(" << fBeta[fIndexBeta[i]] << "), ";
732 }
733 cout << endl;
734 cout << "------------------------------------------------------------------------------------" << endl;
735}
iZ
@ iZ
Definition BmnAlignDefines.h:87
iX
@ iX
Definition BmnAlignDefines.h:85
iY
@ iY
Definition BmnAlignDefines.h:86
z0
const Float_t z0
Definition BmnMwpcHit.cxx:6
i
int i
Definition P4_F32vec4.h:22
min
friend F32vec4 min(const F32vec4 &a, const F32vec4 &b)
Definition P4_F32vec4.h:30
max
friend F32vec4 max(const F32vec4 &a, const F32vec4 &b)
Definition P4_F32vec4.h:31
BmnNdetClusterFinder::GetZ
vector< Double_t > GetZ()
Definition BmnNdetClusterFinder.cxx:169
BmnNdetClusterFinder::Clear
void Clear()
Clear all vectors, set fNCells=0.
Definition BmnNdetClusterFinder.cxx:671
BmnNdetClusterFinder::GetTime
vector< Double_t > GetTime()
Definition BmnNdetClusterFinder.cxx:173
BmnNdetClusterFinder::GetEdep
vector< Double_t > GetEdep()
Definition BmnNdetClusterFinder.cxx:177
BmnNdetClusterFinder::BmnNdetClusterFinder
BmnNdetClusterFinder(const Int_t kNColumns=11, const Int_t kNRows=11, const Int_t kNLayers=16)
Default constructor.
Definition BmnNdetClusterFinder.cxx:15
BmnNdetClusterFinder::SelectNeutrons
Int_t SelectNeutrons(Double_t cut_angle=45., Int_t cut_ncells=3, Bool_t check_border=kFALSE)
Definition BmnNdetClusterFinder.cxx:625
BmnNdetClusterFinder::CalculateNumberOfNeutrons
Int_t CalculateNumberOfNeutrons(Double_t weight_edep=1., Double_t weight_ncells=0.01, Double_t weight_nclusters=0.05, Double_t threshold12=0.)
Definition BmnNdetClusterFinder.cxx:645
BmnNdetClusterFinder::FindClusters4
Int_t FindClusters4(Double_t cut_beta=0.01)
Definition BmnNdetClusterFinder.cxx:519
BmnNdetClusterFinder::FindClusters3
Int_t FindClusters3(Double_t cut_beta=0.03)
Definition BmnNdetClusterFinder.cxx:463
BmnNdetClusterFinder::GetY
vector< Double_t > GetY()
Definition BmnNdetClusterFinder.cxx:165
BmnNdetClusterFinder::AddCluster
BmnNdetCluster * AddCluster()
Definition BmnNdetClusterFinder.cxx:208
BmnNdetClusterFinder::GetClusterStatus
Int_t GetClusterStatus(Int_t i_cluster)
Definition BmnNdetClusterFinder.cxx:193
BmnNdetClusterFinder::GetBeta
vector< Double_t > GetBeta()
Definition BmnNdetClusterFinder.cxx:181
BmnNdetClusterFinder::GetCluster
BmnNdetCluster * GetCluster(Int_t i_cluster)
Definition BmnNdetClusterFinder.cxx:186
BmnNdetClusterFinder::eClusterStatus
eClusterStatus
Status of the cluster.
Definition BmnNdetClusterFinder.h:37
BmnNdetClusterFinder::SetClusterStatus
void SetClusterStatus(Int_t i_cluster, eClusterStatus status)
Definition BmnNdetClusterFinder.cxx:200
BmnNdetClusterFinder::MergeClusters
Int_t MergeClusters(Int_t i_cluster_1, Int_t i_cluster_2)
Definition BmnNdetClusterFinder.cxx:603
BmnNdetClusterFinder::GetX
vector< Double_t > GetX()
Definition BmnNdetClusterFinder.cxx:161
BmnNdetClusterFinder::Fill
void Fill(Int_t cell, Int_t column, Int_t row, Int_t layer, Double_t x, Double_t y, Double_t z, Double_t t, Double_t e_dep)
Definition BmnNdetClusterFinder.cxx:35
BmnNdetClusterFinder::AddCellToCluster
void AddCellToCluster(Int_t i_cell, Int_t i_cluster)
Definition BmnNdetClusterFinder.cxx:217
BmnNdetClusterFinder::~BmnNdetClusterFinder
~BmnNdetClusterFinder()
Default destructor.
Definition BmnNdetClusterFinder.cxx:24
BmnNdetCluster
Class for description of recognized cluster in Bmn Ndet detector.
Definition BmnNdetCluster.h:26
BmnNdetCluster::GetZ
vector< Double_t > GetZ()
Definition BmnNdetCluster.cxx:106
BmnNdetCluster::EnergyDeposited
Double_t EnergyDeposited(Double_t threshold=0.003)
Definition BmnNdetCluster.cxx:166
BmnNdetCluster::GetBeta
vector< Double_t > GetBeta()
Definition BmnNdetCluster.cxx:118
BmnNdetCluster::GetLayers
vector< Int_t > GetLayers()
Definition BmnNdetCluster.cxx:94
BmnNdetCluster::AddCell
void AddCell(Int_t cell, Int_t column, Int_t row, Int_t layer, Double_t x, Double_t y, Double_t z, Double_t t, Double_t e_dep, Double_t beta=-1.)
Definition BmnNdetCluster.cxx:36
BmnNdetCluster::GetTime
vector< Double_t > GetTime()
Definition BmnNdetCluster.cxx:110
BmnNdetCluster::GetY
vector< Double_t > GetY()
Definition BmnNdetCluster.cxx:102
BmnNdetCluster::GetEdep
vector< Double_t > GetEdep()
Definition BmnNdetCluster.cxx:114
BmnNdetCluster::GetColumns
vector< Int_t > GetColumns()
Definition BmnNdetCluster.cxx:86
BmnNdetCluster::AngleToTarget
Double_t AngleToTarget()
Definition BmnNdetCluster.cxx:255
BmnNdetCluster::GetX
vector< Double_t > GetX()
Definition BmnNdetCluster.cxx:98
BmnNdetCluster::StartsOnBorder
Bool_t StartsOnBorder(Int_t size_x=11, Int_t size_y=11, Int_t size_z=8)
Definition BmnNdetCluster.cxx:132
BmnNdetCluster::GetRows
vector< Int_t > GetRows()
Definition BmnNdetCluster.cxx:90
std
STL namespace.