BmnRoot
Loading...
Searching...
No Matches
CbmStsClusterFinder.cxx
Go to the documentation of this file.
1//* $Id: */
2
3// -------------------------------------------------------------------------
4// ----- CbmStsClusterFinder source fil -----
5// ----- Created 26/06/2008 by R. Karabowicz -----
6// -------------------------------------------------------------------------
7// -------------------------------------------------------------------------
8#include "TClonesArray.h"
9#include "TH1S.h"
10
11#include "FairRootManager.h"
12#include "FairRunAna.h"
13#include "FairRuntimeDb.h"
14
15#include "CbmGeoStsPar.h"
16#include "CbmStsCluster.h"
17#include "CbmStsDigi.h"
18#include "CbmStsDigiPar.h"
19#include "CbmStsDigiScheme.h"
20#include "CbmStsClusterFinder.h"
21#include "CbmStsHit.h"
22#include "CbmStsSector.h"
23#include "CbmStsStation.h"
24#include "TMath.h"
25
26#include <iomanip>
27
28using std::cout;
29using std::cerr;
30using std::endl;
31using std::flush;
32using std::fixed;
33using std::right;
34using std::left;
35using std::setw;
36using std::setprecision;
37using std::set;
38using std::map;
39
40// ----- Default constructor ------------------------------------------
41CbmStsClusterFinder::CbmStsClusterFinder()
42 : FairTask("STS Cluster Finder", 1),
43 fGeoPar(NULL),
44 fDigiPar(NULL),
45 fDigiScheme(CbmStsDigiScheme::Instance()),
46 fDigis(NULL),
47 fClustersCand(NULL),
48 fClusters(NULL),
49 fNofDigis(0),
50 fDigiMapF(),
51 fDigiMapB(),
52 fTimer(),
53 fNofClustersCand(0),
54 fNofClusters(0),
55 fNofClustersGood(0),
56 fNofClustersWP(0),
57 fNofClustersWM(0),
58 fLongestCluster(0),
59 fLongestGoodCluster(0)
60{}
61
62// ----- Standard constructor ------------------------------------------
63CbmStsClusterFinder::CbmStsClusterFinder(Int_t iVerbose)
64 : FairTask("STSClusterFinder", iVerbose),
65 fGeoPar(NULL),
66 fDigiPar(NULL),
67 fDigiScheme(CbmStsDigiScheme::Instance()),
68 fDigis(NULL),
69 fClustersCand(NULL),
70 fClusters(NULL),
71 fNofDigis(0),
72 fDigiMapF(),
73 fDigiMapB(),
74 fTimer(),
75 fNofClustersCand(0),
76 fNofClusters(0),
77 fNofClustersGood(0),
78 fNofClustersWP(0),
79 fNofClustersWM(0),
80 fLongestCluster(0),
81 fLongestGoodCluster(0)
82{}
83
84// ----- Constructor with name -----------------------------------------
85CbmStsClusterFinder::CbmStsClusterFinder(const char* name, Int_t iVerbose)
86 : FairTask(name, iVerbose),
87 fGeoPar(NULL),
88 fDigiPar(NULL),
89 fDigiScheme(CbmStsDigiScheme::Instance()),
90 fDigis(NULL),
91 fClustersCand(NULL),
92 fClusters(NULL),
93 fNofDigis(0),
94 fDigiMapF(),
95 fDigiMapB(),
96 fTimer(),
97 fNofClustersCand(0),
98 fNofClusters(0),
99 fNofClustersGood(0),
100 fNofClustersWP(0),
101 fNofClustersWM(0),
102 fLongestCluster(0),
103 fLongestGoodCluster(0)
104{}
105
106// ----- Destructor ----------------------------------------------------
107CbmStsClusterFinder::~CbmStsClusterFinder() {
108 if ( fClustersCand ) {
109 fClustersCand->Delete();
110 delete fClustersCand;
111 }
112 if ( fClusters ) {
113 fClusters->Delete();
114 delete fClusters;
115 }
116}
117// -------------------------------------------------------------------------
118
119
120
121// ----- Public method Exec --------------------------------------------
122void CbmStsClusterFinder::Exec(Option_t* opt) {
123
124 fTimer.Start();
125 Bool_t warn = kFALSE;
126
127 if ( fVerbose ) {
128 cout << endl << "-I- " << fName << ": executing event with " << fDigis->GetEntriesFast() << " digis." << endl;
129 cout << "----------------------------------------------" << endl;
130 }
131
132 // Check for digi scheme
133 if ( ! fDigiScheme ) {
134 cerr << "-E- " << fName << "::Exec: No digi scheme!" << endl;
135 return;
136 }
137
138 // Clear output array
139 // cout << "before clear: " << fClusters->GetEntriesFast() << endl;
140 fNofClustersCand = 0;
141 fNofClusters = 0;
142 // fClustersCand->Clear();
143 // fClusters->Clear();
144 fClusters->Delete();
145 fClustersCand->Delete();
146 // cout << " after clear: " << fClusters->GetEntriesFast() << endl;
147
148 // Sort STS digis with respect to sectors
149 SortDigis();
150
151 // Find hits in sectors
152 Int_t nDigisF = 0;
153 Int_t nDigisB = 0;
154 Int_t nStations = fDigiScheme->GetNStations();
155 CbmStsStation* station = NULL;
156 for (Int_t iStation=0; iStation<nStations; iStation++) {
157 station = fDigiScheme->GetStation(iStation);
158 Int_t nDigisFInStation = 0;
159 Int_t nDigisBInStation = 0;
160 Int_t nSectors = station->GetNSectors();
161
162 for (Int_t iSector=0; iSector<nSectors; iSector++) {
163 CbmStsSector* sector = station->GetSector(iSector);
164 // cout << "=============================================================================" << endl;
165 // cout << "station " << iStation+1 << " sector " << iSector+1 << " " << endl;
166 set <Int_t> fSet, bSet;
167 if ( fDigiMapF.find(sector) == fDigiMapF.end() ) {
168 cout << "-E- " << fName << "::Exec: sector "
169 << sector->GetSectorNr() << " of station "
170 << station->GetStationNr() << "not found in front map!"
171 << endl;
172 warn = kTRUE;
173 continue;
174 }
175 fSet = fDigiMapF[sector];
176 FindClusters(iStation+1,iSector+1,0, fSet);
177 if ( sector->GetType() == 2 || sector->GetType() == 3 ) {
178 if ( fDigiMapB.find(sector) == fDigiMapB.end() ) {
179 cout << "-E- " << fName << "::Exec: sector "
180 << sector->GetSectorNr() << " of station "
181 << station->GetStationNr() << "not found in back map!"
182 << endl;
183 warn = kTRUE;
184 continue;
185 }
186 }
187 bSet = fDigiMapB[sector];
188 FindClusters(iStation+1,iSector+1,1, bSet);
189 Int_t nDigisFInSector = fSet.size();
190 Int_t nDigisBInSector = bSet.size();
191 nDigisFInStation += nDigisFInSector;
192 nDigisBInStation += nDigisBInSector;
193 } // Sector loop
194
195 nDigisF += nDigisFInStation;
196 nDigisB += nDigisBInStation;
197
198 } // Station loop
199
200 // analyze clusters...
201 AnalyzeClusters();
202
203 fTimer.Stop();
204 if ( fVerbose ) {
205 cout << endl;
206 cout << "-I- " << fName << ":Event summary" << endl;
207 cout << " Active channels front side: " << nDigisF << endl;
208 cout << " Active channels back side : " << nDigisB << endl;
209 cout << " Real time : " << fTimer.RealTime()
210 << endl;
211 }
212 else {
213 if ( warn ) cout << "- ";
214 else cout << "+ ";
215 cout << setw(15) << left << fName << ": " << setprecision(4) << setw(8)
216 << fixed << right << fTimer.RealTime()
217 << " s, " << fNofClustersCand
218 << "(" << fNofClustersGood
219 << "+" << fNofClustersWP
220 << ") clusters, longest till now " << fLongestCluster << endl;
221 // cout << fDigis->GetEntriesFast() << " vs " << fClusters->GetEntriesFast() << endl;
222 }
223
224
225 //AZ - add links to clusters (digits)
226 Int_t nClust = fClusters->GetEntriesFast();
227 for (Int_t i = 0; i < nClust; ++i) {
228 CbmStsCluster *clus = (CbmStsCluster*) fClusters->UncheckedAt(i);
229
230 Int_t nDigis = clus->GetNDigis();
231 for (Int_t j = 0; j < nDigis; ++j) {
232 Int_t idigi = clus->GetDigi(j);
233 clus->AddLink(FairLink(kStsDigi, idigi));
234 }
235 }
236
237}
238// -------------------------------------------------------------------------
239
240
241
242
243// ----- Private method SetParContainers -------------------------------
244void CbmStsClusterFinder::SetParContainers() {
245
246 // Get run and runtime database
247 FairRunAna* run = FairRunAna::Instance();
248 if ( ! run ) Fatal("SetParContainers", "No analysis run");
249
250 FairRuntimeDb* db = run->GetRuntimeDb();
251 if ( ! db ) Fatal("SetParContainers", "No runtime database");
252
253 // Get STS geometry parameter container
254 fGeoPar = (CbmGeoStsPar*) db->getContainer("CbmGeoStsPar");
255
256 // Get STS digitisation parameter container
257 fDigiPar = (CbmStsDigiPar*) db->getContainer("CbmStsDigiPar");
258
259}
260// -------------------------------------------------------------------------
261
262
263
264
265// ----- Private method Init -------------------------------------------
266InitStatus CbmStsClusterFinder::Init() {
267
268 // Get input array
269 FairRootManager* ioman = FairRootManager::Instance();
270 if ( ! ioman ) Fatal("Init", "No FairRootManager");
271 fDigis = (TClonesArray*) ioman->GetObject("StsDigi");
272
273 fClustersCand = new TClonesArray("CbmStsCluster", 30000);
274 ioman->Register("StsClusterCand", "Cluster in STS", fClustersCand, kTRUE);
275
276 fClusters = new TClonesArray("CbmStsCluster", 30000);
277 ioman->Register("StsCluster", "Cluster in STS", fClusters, kTRUE);
278
279 // Build digitisation scheme
280 Bool_t success = fDigiScheme->Init(fGeoPar, fDigiPar);
281 if ( ! success ) return kERROR;
282
283 // Create sectorwise digi sets
284 MakeSets();
285
286 // Control output
287 if (fVerbose == 1 || fVerbose == 2) fDigiScheme->Print(kFALSE);
288 else if (fVerbose > 2) fDigiScheme->Print(kTRUE);
289 cout << "-I- " << fName << "::Init: "
290 << "STS digitisation scheme succesfully initialised" << endl;
291 cout << " Stations: " << fDigiScheme->GetNStations()
292 << ", Sectors: " << fDigiScheme->GetNSectors() << ", Channels: "
293 << fDigiScheme->GetNChannels() << endl;
294
295 fNofClustersCand = 0;
296 fNofClusters = 0;
297 fNofClustersGood = 0;
298 fNofClustersWP = 0;
299 fNofClustersWM = 0;
300
301 fLongestCluster = 0;
302 fLongestGoodCluster = 0;
303
304 return kSUCCESS;
305}
306// -------------------------------------------------------------------------
307
308
309
310
311// ----- Private method ReInit -----------------------------------------
312InitStatus CbmStsClusterFinder::ReInit() {
313
314 // Clear digitisation scheme
315 fDigiScheme->Clear();
316
317 // Build new digitisation scheme
318 Bool_t success = fDigiScheme->Init(fGeoPar, fDigiPar);
319 if ( ! success ) return kERROR;
320
321 // Create sectorwise digi sets
322 MakeSets();
323
324 // Control output
325 if (fVerbose == 1 || fVerbose == 2) fDigiScheme->Print(kFALSE);
326 else if (fVerbose > 2) fDigiScheme->Print(kTRUE);
327 cout << "-I- " << fName << "::Init: "
328 << "STS digitisation scheme succesfully reinitialised" << endl;
329 cout << " Stations: " << fDigiScheme->GetNStations()
330 << ", Sectors: " << fDigiScheme->GetNSectors() << ", Channels: "
331 << fDigiScheme->GetNChannels() << endl;
332
333 fNofClustersCand = 0;
334 fNofClusters = 0;
335 fNofClustersGood = 0;
336 fNofClustersWP = 0;
337 fNofClustersWM = 0;
338
339 fLongestCluster = 0;
340 fLongestGoodCluster = 0;
341
342 return kSUCCESS;
343}
344// -------------------------------------------------------------------------
345
346
347
348
349// ----- Private method MakeSets ---------------------------------------
350void CbmStsClusterFinder::MakeSets() {
351
352 fDigiMapF.clear();
353 fDigiMapB.clear();
354 Int_t nStations = fDigiScheme->GetNStations();
355 for (Int_t iStation=0; iStation<nStations; iStation++) {
356 CbmStsStation* station = fDigiScheme->GetStation(iStation);
357 Int_t nSectors = station->GetNSectors();
358 for (Int_t iSector=0; iSector<nSectors; iSector++) {
359 CbmStsSector* sector = station->GetSector(iSector);
360 set<Int_t> a;
361 fDigiMapF[sector] = a;
362 if ( sector->GetType() == 2 || sector->GetType() ==3 ) {
363 set<Int_t> b;
364 fDigiMapB[sector] = b;
365 }
366 }
367 }
368
369}
370// -------------------------------------------------------------------------
371
372
373
374
375// ----- Private method SortDigis --------------------------------------
376void CbmStsClusterFinder::SortDigis() {
377
378 // Check input array
379 if ( ! fDigis ) {
380 cout << "-E- " << fName << "::SortDigis: No input array!" << endl;
381 return;
382 }
383
384 CbmStsDigi* digi = NULL;
385 const Int_t nDigis = fDigis->GetEntriesFast();
386
387 // Clear sector digi sets
388 map<CbmStsSector*, set<Int_t> >::iterator mapIt;
389 for (mapIt=fDigiMapF.begin(); mapIt!=fDigiMapF.end(); mapIt++)
390 ((*mapIt).second).clear();
391 for (mapIt=fDigiMapB.begin(); mapIt!=fDigiMapB.end(); mapIt++)
392 ((*mapIt).second).clear();
393
394 // Fill digis into sets
395 CbmStsSector* sector = NULL;
396 Int_t stationNr = -1;
397 Int_t sectorNr = -1;
398 Int_t iSide = -1;
399 //Int_t channelNr = -1;
400
401 for (Int_t iDigi=0; iDigi<nDigis; iDigi++) {
402 digi = (CbmStsDigi*) fDigis->At(iDigi);
403
404 stationNr = digi->GetStationNr();
405 sectorNr = digi->GetSectorNr();
406 iSide = digi->GetSide();
407 sector = fDigiScheme->GetSector(stationNr, sectorNr);
408
409 if (iSide == 0 ) {
410 if ( fDigiMapF.find(sector) == fDigiMapF.end() ) {
411 cerr << "-E- " << fName << "::SortDigits:: sector " << sectorNr
412 << " of station " << stationNr
413 << " not found in digi scheme (F)!" << endl;
414 continue;
415 }
416 fDigiMapF[sector].insert(iDigi);
417 }
418 else if (iSide == 1 ) {
419 if ( fDigiMapB.find(sector) == fDigiMapB.end() ) {
420 cerr << "-E- " << fName << "::SortDigits:: sector " << sectorNr
421 << " of station " << stationNr
422 << " not found in digi scheme (B)!" << endl;
423 continue;
424 }
425 fDigiMapB[sector].insert(iDigi);
426 }
427 }
428}
429// -------------------------------------------------------------------------
430
431// ----- Private method RealisticResponse ------------------------------
432Int_t CbmStsClusterFinder::FindClusters(Int_t stationNr, Int_t sectorNr, Int_t iSide, set<Int_t>& digiSet) {
433
434 set<Int_t>::iterator it1;
435
436 Int_t iDigi = -1;
437 CbmStsDigi* digi = NULL;
438 CbmStsCluster* clusterCand = NULL;
439 //CbmStsCluster* cluster = NULL;
440
441 Int_t digiPos = -1;
442 Double_t digiSig = -1.;
443 Int_t lastDigiNr = -1;
444 Int_t lastDigiPos = -1;
445 Double_t lastDigiSig = 100000.;
446 Int_t clusterMaxNr = -1;
447 //Int_t clusterMaxPos = -1;
448 Double_t clusterMaxSig = -1.;
449 Bool_t clusterHasMinimum = kFALSE;
450 Bool_t clusterHasPlateau = kFALSE;
451 Int_t clusterBeginPos = 0;
452 Int_t clusterWidth = 0;
453 //cout << "====================================================================" << digiSet.size() << endl;
454 //cout << " cluster " << flush;
455 // if ( fNofClusters != fNofClustersGood+fNofClustersWP ) cout << fNofClusters-(fNofClustersGood+fNofClustersWP) << " -> " << fNofClusters <<"!="<< fNofClustersGood<<"+"<<fNofClustersWP<<endl;
456 for (it1=digiSet.begin(); it1!=digiSet.end(); it1++) {
457 iDigi = (*it1);
458 digi = (CbmStsDigi*) fDigis->At(iDigi);
459
460 digiPos = digi->GetChannelNr();
461 digiSig = digi->GetAdc();
462
463 if ( lastDigiPos == -1 ) {
464 clusterCand = new ((*fClustersCand)[fNofClustersCand++]) CbmStsCluster(digiSig, stationNr,sectorNr,iSide);
465 // cout << "first cluster ADC " << digiSig << endl;
466 clusterBeginPos = digiPos;
467 }
468
469 if ( digiPos == lastDigiPos+1 ) {
470
471 if ( digiSig == lastDigiSig ) {
472 clusterCand->SetMeanError(digiSig);
473 clusterHasPlateau = kTRUE;
474 }
475
476 // if the signal is larger that last one and the previous one is smaller than maximum
477 if ( digiSig > lastDigiSig && lastDigiSig < clusterMaxSig && digiSig != clusterMaxSig) {
478 clusterCand->SetMean(clusterMaxNr);
479 // cluster->SetMeanError(clusterMaxSig);
480 clusterCand = new ((*fClustersCand)[fNofClustersCand++]) CbmStsCluster(digiSig, stationNr,sectorNr,iSide);
481 clusterCand->AddDigi(lastDigiNr);
482 // cout << " +end cluster " << lastDigiPos << endl;
483
484 // cout << "+new cluster " << digiPos << endl;
485 clusterWidth = lastDigiPos - clusterBeginPos + 1;
486 if ( clusterWidth > fLongestCluster )
487 fLongestCluster = clusterWidth;
488 if ( clusterHasPlateau ) {
489 fNofClustersWP++;
490 }
491 if ( !clusterHasMinimum && !clusterHasPlateau ) {
492 fNofClustersGood++;
493 if ( clusterWidth > fLongestGoodCluster )
494 fLongestGoodCluster = clusterWidth;
495 }
496
497 clusterHasPlateau = kFALSE;
498 //clusterMaxPos = -1;
499 clusterMaxSig = -1.;
500 clusterBeginPos = digiPos;
501 }
502 }
503 else if ( lastDigiPos>=0 ) {
504 clusterCand->SetMean(clusterMaxNr);
505 // cluster->SetMeanError(clusterMaxSig);
506 clusterCand = new ((*fClustersCand)[fNofClustersCand++]) CbmStsCluster(digiSig, stationNr,sectorNr,iSide);
507 // cout << " end cluster " << lastDigiPos << endl;
508 // cout << "new cluster " << digiPos << endl;
509 clusterWidth = lastDigiPos - clusterBeginPos + 1;
510 if ( clusterWidth > fLongestCluster )
511 fLongestCluster = clusterWidth;
512 if ( clusterHasPlateau ) {
513 fNofClustersWP++;
514 }
515 if ( !clusterHasMinimum && !clusterHasPlateau ) {
516 fNofClustersGood++;
517 if ( clusterWidth > fLongestGoodCluster )
518 fLongestGoodCluster = clusterWidth;
519 }
520
521 clusterHasPlateau = kFALSE;
522 //clusterMaxPos = -1;
523 clusterMaxSig = -1.;
524 clusterBeginPos = digiPos;
525 }
526 if ( digiSig > clusterMaxSig ) {
527 clusterMaxNr = iDigi;
528 //clusterMaxPos = digiPos;
529 clusterMaxSig = digiSig;
530 }
531
532 clusterWidth = lastDigiPos - clusterBeginPos + 1;
533 if ( clusterWidth < 90 ) {
534 clusterCand->AddDigi(iDigi);}
535 // clusterCand->AddIndex(iDigi,digiSig);}
536 else {
537 // return 1;
538 clusterCand->SetMean(clusterMaxNr);
539 // cluster->SetMeanError(clusterMaxSig);
540
541 clusterCand = new ((*fClustersCand)[fNofClustersCand++]) CbmStsCluster(digiSig, stationNr,sectorNr,iSide);
542 clusterCand->AddDigi(iDigi);
543 clusterWidth = lastDigiPos - clusterBeginPos + 1;
544 if ( clusterWidth > fLongestCluster )
545 fLongestCluster = clusterWidth;
546 if ( clusterHasPlateau ) {
547 fNofClustersWP++;
548 }
549 if ( !clusterHasMinimum && !clusterHasPlateau ) {
550 fNofClustersGood++;
551 if ( clusterWidth > fLongestGoodCluster )
552 fLongestGoodCluster = clusterWidth;
553 }
554
555 clusterHasPlateau = kFALSE;
556 //clusterMaxPos = -1;
557 clusterMaxSig = -1.;
558 clusterBeginPos = digiPos;
559 }
560
561 lastDigiNr = iDigi;
562 lastDigiPos = digiPos;
563 lastDigiSig = digiSig;
564 }
565 if ( digiSig > 0. ) {
566 clusterCand->SetMean(clusterMaxNr);
567
568 // cout << " last cluster " << lastDigiPos << endl;
569 clusterWidth = lastDigiPos - clusterBeginPos + 1;
570 if ( clusterWidth > fLongestCluster )
571 fLongestCluster = clusterWidth;
572 if ( clusterHasPlateau ) {
573 fNofClustersWP++;
574 }
575 if ( lastDigiPos && !clusterHasPlateau ) {
576 fNofClustersGood++;
577 if ( clusterWidth > fLongestGoodCluster )
578 fLongestGoodCluster = clusterWidth;
579 }
580 }
581 // delete clusterCand;
582
583 return 1;
584
585 map<Int_t , Int_t> channelSorted;
586 for (it1=digiSet.begin(); it1!=digiSet.end(); it1++) {
587 iDigi = (*it1);
588 digi = (CbmStsDigi*) fDigis->At(iDigi);
589 cout << digi->GetChannelNr() << " " << flush;
590 channelSorted[digi->GetChannelNr()] = iDigi+1;
591 }
592 cout << endl;
593 // print channels/signals
594 size_t iFS = 0;
595 for (Int_t iter=0; iter<1100; iter++) {
596 iDigi = channelSorted[iter];
597 if ( iDigi == 0 ) continue;
598 iFS++;
599 digi = (CbmStsDigi*) fDigis->At(iDigi-1);
600 if (iFS >= channelSorted.size()) break;
601 }
602 // cout << endl;
603
604 // cout << "size = " << channelSorted.size() << endl;
605 /* Int_t iFS = 0;
606 Float_t lastSignal = 0.;
607 for (Int_t iter=0; iter<1100; iter++) {
608 iDigi = channelSorted[iter];
609 if ( iDigi == 0 ) {lastSignal=0.;continue;}
610 iFS++;
611 digi = (CbmStsDigi*) fDigis->At(iDigi-1);
612 Float_t signal = digi->GetADC();
613 if ( signal <= lastSignal ) { lastSignal = signal; continue; }
614
615 if ( iFS >= channelSorted.size() ) break;
616
617
618 }*/
619 // cout << endl << "=============================================" << endl;
620
621}
622// -------------------------------------------------------------------------
623
624// ----- Method AnalyzeClusters ----------------------------------------
625void CbmStsClusterFinder::AnalyzeClusters() {
626 for ( Int_t iclus = 0 ; iclus < fNofClustersCand ; iclus++ ) {
627 AnalyzeCluster(iclus);
628 }
629}
630// -------------------------------------------------------------------------
631
632// ----- Method AnalyzeClusters ----------------------------------------
633void CbmStsClusterFinder::AnalyzeCluster(Int_t iCluster) {
634
635 CbmStsCluster* clusterCand = (CbmStsCluster*) fClustersCand->At(iCluster);
636 CbmStsCluster* cluster = (CbmStsCluster*) fClusters->At(iCluster);
637 // CbmStsSector* sector = NULL;
638 Int_t maxDigiNr = (Int_t)clusterCand->GetMean();
639 //Double_t plateau = clusterCand->GetMeanError();
640 Double_t maxDigiSig = 0.;
641 CbmStsDigi* digi = NULL;
642 digi = (CbmStsDigi*)fDigis->At(maxDigiNr);
643 //Int_t maxDigiPos = digi->GetChannelNr();
644 maxDigiSig = digi->GetAdc();
645
646 if (clusterCand->GetNDigis() > 5) {
647 // AZ - split large clusters
648 //SplitCluster(iCluster);
649 //return;
650 }
651 // cout << "Cluster " << iCluster+1 << " has " << cluster->GetNDigis() << " digis, max at " << maxDigiNr << endl;
652 Double_t chanNr = 0;
653 Double_t chanADC = 0.;
654 Double_t sumW = 0;
655 Double_t sumWX = 0;
656 Double_t sumCh = 0;
657 Double_t error = 0;
658 Double_t sumWX2 = 0.0;
659 for ( Int_t itemp = 0 ; itemp < clusterCand->GetNDigis() ; itemp++ ) {
660 digi = (CbmStsDigi*)fDigis->At(clusterCand->GetDigi(itemp));
661 chanNr = (Double_t)digi->GetChannelNr();
662 chanADC = digi->GetAdc();
663 sumW += chanADC;
664 sumWX += chanNr*chanADC;
665 sumCh += chanNr;
666 error += ( chanADC*chanADC );
667 sumWX2 += chanADC * chanNr * chanNr; //AZ
668 // cout << chanADC << " + " << flush;
669 // cout << "channel " << digi->GetChannelNr() << " with signal = " << digi->GetADC() << " (" << sumWX << "/" << sumW << ") - plateau = " << plateau << endl;
670 }
671
672 if (sumW>50.) {
673 //Int_t imean = TMath::Nint (sumWX/sumW); //AZ
674 for ( Int_t itemp = 0 ; itemp < clusterCand->GetNDigis() ; itemp++ ) {
675
676 digi = (CbmStsDigi*)fDigis->At(clusterCand->GetDigi(itemp));
677 if (itemp==0) {
678 cluster = new ((*fClusters)[fNofClusters++]) CbmStsCluster(digi->GetAdc(), digi->GetStationNr(), digi->GetSectorNr(), digi->GetSide());
679 //AZ cluster->AddDigi(itemp);
680 cluster->AddDigi(clusterCand->GetDigi(itemp));
681 //if (digi->GetChannelNr() == imean) cluster->AddDigi(clusterCand->GetDigi(itemp));
682 }
683 else if (itemp>0) {
684 //AZ cluster->AddDigi(itemp);
685 cluster->AddDigi(clusterCand->GetDigi(itemp));
686 //if (digi->GetChannelNr() == imean) cluster->AddDigi(clusterCand->GetDigi(itemp));
687 }
688 }
689 // cout << " mean at = " << sumWX/sumW << endl;
690 cluster->SetMean(sumWX/sumW);
691 cluster->SetQtot(sumW); //AZ
692 cluster->SetMeanError( (1./(sumW)) * TMath::Sqrt(error) );
693 if ( sumW < maxDigiSig ) {
694 cout << " MAX DIGI = " << maxDigiSig << ", while SUMW = " << sumW << endl;
695 for ( Int_t itemp = 0 ; itemp < clusterCand->GetNDigis() ; itemp++ ) {
696 digi = (CbmStsDigi*)fDigis->At(clusterCand->GetDigi(itemp));
697 cout << "digi ADC = " << digi->GetAdc() << " at channel# " << digi->GetChannelNr() << endl;
698 }
699 }
700 // AZ - correct coordinate for 2 adjacent clusters, i.e. the ones sharing 1 channel
701 //if (0) {
702 if (fNofClusters > 1) {
703 CbmStsCluster *clusPrev = (CbmStsCluster*) fClusters->UncheckedAt(fNofClusters-2);
704 if (clusPrev->GetDetectorId() == cluster->GetDetectorId()) {
705
706 Int_t nDigis = clusPrev->GetNDigis();
707 CbmStsDigi *digPrev = (CbmStsDigi*) fDigis->UncheckedAt(clusPrev->GetDigi(nDigis-1));
708 Int_t chan2 = digPrev->GetChannelNr();
709 if (chan2 == ((CbmStsDigi*)fDigis->UncheckedAt(cluster->GetDigi(0)))->GetChannelNr()) {
710
711 Double_t adcPrev = digPrev->GetAdc();
712 Double_t qPrev = clusPrev->GetQtot();
713 Double_t xPrev = clusPrev->GetMean() * qPrev;
714 xPrev -= chan2 * adcPrev * 0.5;
715 xPrev /= (qPrev - adcPrev * 0.5);
716 clusPrev->SetMean(xPrev);
717 clusPrev->SetQtot(qPrev - adcPrev * 0.5);
718 Double_t xCor = sumWX;
719 xCor -= chan2 * adcPrev * 0.5;
720 xCor /= (sumW - adcPrev * 0.5);
721 cluster->SetMean(xCor);
722 cluster->SetQtot(sumW - adcPrev * 0.5);
723 }
724 }
725 }
726 //* AZ - correct coordinate and choose error estimate for previous cluster
727 Int_t chan2 = -1;
728 if (fNofClusters > 1) {
729 CbmStsCluster *clusPrev = (CbmStsCluster*) fClusters->UncheckedAt(fNofClusters-2);
730 if (clusPrev->GetDetectorId() == cluster->GetDetectorId()) {
731 Int_t nDigis = clusPrev->GetNDigis();
732 CbmStsDigi *digPrev = (CbmStsDigi*) fDigis->UncheckedAt(clusPrev->GetDigi(nDigis-1));
733 chan2 = digPrev->GetChannelNr();
734 if (chan2 != ((CbmStsDigi*)fDigis->UncheckedAt(cluster->GetDigi(0)))->GetChannelNr()) chan2 = -1;
735 }
736 EvalErrors(clusPrev, chan2);
737 }
738 //Double_t rms = sumWX2 / sumW - cluster->GetMean() * cluster->GetMean();
739 Double_t rms = (sumWX2 - sumWX * sumWX / sumW) / sumW;
740 cluster->SetMeanError(TMath::Sqrt(rms));
741 //cluster->SetMeanError(0.04/TMath::Sqrt(12.)); //rms);
742 cluster->SetLeft(chan2); // if > 0 - there is a left neighbour
743 if (iCluster == fNofClustersCand-1) EvalErrors(cluster, -1); // last cluster
744 //AZ
745 //*/
746 } // if (sumW>50.)
747 // cout << sumWE/sumW << " mean at " << endl;
748 // cout << "error = " << sumW/maxDigiSig << endl;
749}
750// -------------------------------------------------------------------------
751
752//AZ -----------------------------------------------------------------------
753void CbmStsClusterFinder::SplitCluster(Int_t iCluster)
754{
755 // Split large cluster into subclusters
756
757 const Int_t nDigMax = 13;
758 static Int_t first = 1, nSubs[nDigMax] = {0};
759
760 if (first) {
761 first = 0;
762 nSubs[6] = nSubs[7] = nSubs[8] = nSubs[9] = 2; // number of subclusters
763 nSubs[10] = nSubs[11] = nSubs[12] = 3;
764 }
765
766 CbmStsCluster* clusterCand = (CbmStsCluster*) fClustersCand->At(iCluster);
767 CbmStsCluster* cluster = NULL;
768 CbmStsDigi *digi = NULL;
769
770 Int_t mult = clusterCand->GetNDigis(), nsub = 1;
771 if (mult < nDigMax) nsub = nSubs[mult];
772 else nsub = TMath::Nint (mult * 1.0 / 4 - 0.1);
773 Int_t leng0 = TMath::Nint (mult * 1.0 / nsub);
774 Int_t i0 = 0, itemp = 0;
775 Double_t chanNr = 0.0, chanADC = 0.0, sumW = 0.0, sumWX = 0.0, sumWX2 = 0.0;
776
777 //for ( Int_t itemp = 0; itemp < mult; itemp++ ) {
778 for (Int_t icl = 0; icl < nsub; ++icl) {
779 Int_t leng = leng0;
780 if (icl == 0 && mult - leng0 * nsub == 1) ++leng; // add extra digi to first subcluster
781 else if (icl == nsub - 1 && mult - leng0 * nsub == 2) ++leng; // add extra digi to last subclus.
782 while (1) {
783 digi = (CbmStsDigi*) fDigis->UncheckedAt(clusterCand->GetDigi(itemp));
784 if (i0 == 0) {
785 cluster = new ((*fClusters)[fNofClusters++]) CbmStsCluster(digi->GetAdc(), digi->GetStationNr(), digi->GetSectorNr(), digi->GetSide());
786 cluster->AddDigi(clusterCand->GetDigi(itemp));
787 } else {
788 cluster->AddDigi(clusterCand->GetDigi(itemp));
789 }
790 ++i0;
791 ++itemp;
792 chanNr = digi->GetChannelNr();
793 chanADC = digi->GetAdc();
794 sumW += chanADC;
795 sumWX += chanNr * chanADC;
796 sumWX2 += chanADC * chanNr * chanNr;
797 if (i0 == leng) {
798 // Subcluster
799 cluster->SetMean(sumWX/sumW);
800 cluster->SetQtot(sumW); //AZ
801 cluster->SetDefaultType(mult);
802
803 // AZ - correct coordinate for 2 adjacent clusters, i.e. the ones sharing 1 channel or back to back
804 //if (0) {
805 if (fNofClusters > 1) {
806 CbmStsCluster *clusPrev = (CbmStsCluster*) fClusters->UncheckedAt(fNofClusters-2);
807 if (clusPrev->GetDetectorId() == cluster->GetDetectorId()) {
808
809 Int_t nDigis = clusPrev->GetNDigis();
810 CbmStsDigi *digPrev = (CbmStsDigi*) fDigis->UncheckedAt(clusPrev->GetDigi(nDigis-1));
811 Int_t chan2 = digPrev->GetChannelNr();
812 CbmStsDigi *digThis = (CbmStsDigi*) fDigis->UncheckedAt(cluster->GetDigi(0));
813 Int_t chan1 = digThis->GetChannelNr();
814 if (TMath::Abs(chan2-chan1) < 1) {
815 // Adjacent clusters
816 Double_t adcPrev = digPrev->GetAdc();
817 Double_t qPrev = clusPrev->GetQtot();
818 Double_t xPrev = clusPrev->GetMean() * qPrev;
819 xPrev -= chan2 * adcPrev * 0.5;
820 xPrev /= (qPrev - adcPrev * 0.5);
821 clusPrev->SetMean(xPrev);
822 clusPrev->SetQtot(qPrev - adcPrev * 0.5);
823 adcPrev = digThis->GetAdc();
824 Double_t xCor = sumWX;
825 xCor -= chan1 * adcPrev * 0.5;
826 xCor /= (sumW - adcPrev * 0.5);
827 cluster->SetMean(xCor);
828 cluster->SetQtot(sumW - adcPrev * 0.5);
829 }
830 }
831 }
832 //* AZ - correct coordinate and choose error estimate for previous cluster
833 Int_t chan2 = -1;
834 if (fNofClusters > 1) {
835 CbmStsCluster *clusPrev = (CbmStsCluster*) fClusters->UncheckedAt(fNofClusters-2);
836 if (clusPrev->GetDetectorId() == cluster->GetDetectorId()) {
837 Int_t nDigis = clusPrev->GetNDigis();
838 CbmStsDigi *digPrev = (CbmStsDigi*) fDigis->UncheckedAt(clusPrev->GetDigi(nDigis-1));
839 chan2 = digPrev->GetChannelNr();
840 if (TMath::Abs(chan2-((CbmStsDigi*)fDigis->UncheckedAt(cluster->GetDigi(0)))->GetChannelNr()) > 0) chan2 = -1;
841 }
842 EvalErrors(clusPrev, chan2);
843 }
844 //Double_t rms = sumWX2 / sumW - cluster->GetMean() * cluster->GetMean();
845 Double_t rms = (sumWX2 - sumWX * sumWX / sumW) / sumW;
846 cluster->SetMeanError(TMath::Sqrt(rms));
847 //cluster->SetMeanError(0.04/TMath::Sqrt(12.)); //rms);
848 cluster->SetLeft(chan2); // if > 0 - there is a left neighbour
849 if (iCluster == fNofClustersCand-1) EvalErrors(cluster, -1); // last cluster
850
851 //cluster->SetMeanError( (1./(sumW)) * TMath::Sqrt(error) );
852 i0 = 0;
853 if (itemp + leng * (nsub - icl - 1) > mult) --itemp; // move one digi backward
854 sumW = sumWX = sumWX2 = 0.0;
855 break;
856 }
857 }
858 }
859}
860// -------------------------------------------------------------------------
861
862//AZ -----------------------------------------------------------------------
863void CbmStsClusterFinder::EvalErrors(CbmStsCluster *clus, Int_t chan2) {
864 // Evaluate errors and correct coordinate if necessary
865
866 Int_t mult = clus->GetNDigis();
867 Double_t sigma = 0.0, xcor = clus->GetMean(), rms = clus->GetMeanError();
868
869 if (chan2 >= 0 && clus->GetLeft() >= 0) {
870 // Cluster has 2 neighbours
871 if (mult <= 3) sigma = 0.0298;
872 else if (mult == 4) {
873 if (clus->GetDefaultType()) sigma = 0.0398; // split cluster - check that !
874 else sigma = 0.0135; // fit
875 }
876 else if (mult == 5) {
877 if (clus->GetDefaultType()) sigma = 0.0463; // split cluster - check that !
878 else sigma = 0.0109; // fit
879 }
880 }
881
882 else if (chan2 >= 0 || clus->GetLeft() >= 0) {
883 // 1 neighbour
884 //Int_t ishft = 1;
885 //if (chan2 >= 0) ishft = -1;
886 if (mult <= 2) {
887 sigma = 0.0220;
888 //xcor = xcor - ishft * 0.0122 / 0.0400;
889 }
890 else if (mult == 5) {
894 //sigma = 0.0055; // fit
895 //xcor = xcor + ishft * (0.0070 / 0.0400 + (rms - 1.234) * 0.44 / 4.5 / 0.0400);
896 if (clus->GetDefaultType()) sigma = 0.0487; // split cluster - check that !
897 else sigma = 0.0118; // fit
898 }
899 else if (mult == 3) {
900 //sigma = 0.0080;
901 //xcor = xcor + ishft * (-0.0025 / 0.0400 + (rms - 0.776) * 0.54 / 2.5 / 0.0400);
902 if (clus->GetDefaultType()) sigma = 0.0315; // split cluster - check that !
903 else sigma = 0.0098; // fit
904 }
905 else if (mult == 4) {
906 //sigma = 0.0065;
907 //xcor = xcor + ishft * (0.0045 / 0.0400 + (rms - 1.028) * 0.62 / 3.5 / 0.0400);
908 if (clus->GetDefaultType()) sigma = 0.0366; // split cluster - check that !
909 else sigma = 0.0100; // fit
910 }
911 xcor = TMath::Max (0.0, xcor);
912 }
913
914 else {
915 // No neighbours
916 if (mult == 1) sigma = 0.0122;
917 else if (mult == 2) {
918 if (rms > 0.495) sigma = 0.0071; // RMS
919 else sigma = 0.0224; // edge effect - fix !
920 }
921 else if (mult == 3) {
922 if (clus->GetDefaultType()) sigma = 0.0255; // split cluster - check that !
923 else sigma = 0.0064; // fit
924 }
925 else if (mult == 4) {
926 if (clus->GetDefaultType()) sigma = 0.0337; // split cluster - check that !
927 sigma = 0.0048; // fit
928 }
929 else if (mult == 5) {
930 //if (rms < 1.19) sigma = 0.0108;
931 if (rms < 1.19) sigma = 0.0038; // fit
932 else sigma = 0.0415; // should be split also !
933 }
934 }
935
936 clus->SetMean(xcor);
937 clus->SetMeanError(sigma);
938}
939//AZ -----------------------------------------------------------------------
940
941// ----- Virtual method Finish -----------------------------------------
942void CbmStsClusterFinder::Finish() {
943 fNofClustersCand = fNofClustersGood+fNofClustersWP;
944 cout << endl;
945 cout << "============================================================"
946 << endl;
947 cout << "===== " << fName << ": Run summary " << endl;
948 cout << "===== " << endl;
949 cout << "===== Number of clusters : "
950 << setw(8) << setprecision(2)
951 << fNofClustersCand << endl;
952 cout << "===== Number of good clusters : "
953 << setw(8) << setprecision(2)
954 << fNofClustersGood << " ("
955 << setw(8) << setprecision(2)
956 << 100.*fNofClustersGood/fNofClustersCand << "%)" << endl;
957 cout << "===== Number of plateau clusters : "
958 << setw(8) << setprecision(2)
959 << fNofClustersWP << " ("
960 << setw(8) << setprecision(2)
961 << 100.*fNofClustersWP/fNofClustersCand << "%)" << endl;
962 cout << "===== Number of minimum clusters : "
963 << setw(8) << setprecision(2)
964 << fNofClustersWM << " ("
965 << setw(8) << setprecision(2)
966 << 100.*fNofClustersWM/fNofClustersCand << "%)" << endl;
967 cout << "===== Longest cluster : "
968 << setw(8) << setprecision(2)
969 << fLongestCluster << endl;
970 cout << "===== Longest good cluster : "
971 << setw(8) << setprecision(2)
972 << fLongestGoodCluster << endl;
973 cout << "============================================================"
974 << endl;
975
976}
977// -------------------------------------------------------------------------
fDigiPar
BmnSsdDigitizeParameters * fDigiPar
Digitisation parameters.
Definition BmnSsdDigitize.cxx:65
i
int i
Definition P4_F32vec4.h:22
kStsDigi
@ kStsDigi
Definition BmnDetectorList.h:82
CbmStsClusterFinder::Exec
virtual void Exec(Option_t *opt)
Definition CbmStsClusterFinder.cxx:122
CbmStsCluster::SetMeanError
void SetMeanError(Double_t err)
Definition CbmStsCluster.h:57
CbmStsCluster::GetLeft
Int_t GetLeft() const
Definition CbmStsCluster.h:74
CbmStsCluster::SetMean
void SetMean(Double_t chan)
Definition CbmStsCluster.h:56
CbmStsCluster::GetMean
Double_t GetMean() const
Definition CbmStsCluster.h:71
CbmStsCluster::SetQtot
void SetQtot(Double_t qtot)
Definition CbmStsCluster.h:75
CbmStsCluster::GetNDigis
Int_t GetNDigis() const
Definition CbmStsCluster.h:69
CbmStsCluster::GetDetectorId
Int_t GetDetectorId() const
Definition CbmStsCluster.h:60
CbmStsCluster::GetMeanError
Double_t GetMeanError() const
Definition CbmStsCluster.h:72
CbmStsCluster::GetQtot
Double_t GetQtot() const
Definition CbmStsCluster.h:73
CbmStsCluster::GetDigi
Int_t GetDigi(Int_t inum)
Definition CbmStsCluster.h:70
CbmStsCluster::SetLeft
void SetLeft(Int_t left)
Definition CbmStsCluster.h:76
CbmStsCluster::AddDigi
void AddDigi(Int_t idigi)
Definition CbmStsCluster.h:52
CbmStsDigiScheme::Print
void Print(Bool_t kLong=kFALSE)
Definition CbmStsDigiScheme.cxx:248
CbmStsDigiScheme::GetSector
CbmStsSector * GetSector(Int_t stationNr, Int_t sectorNr)
Definition CbmStsDigiScheme.cxx:204
CbmStsDigiScheme::GetStation
CbmStsStation * GetStation(Int_t iStation)
Definition CbmStsDigiScheme.cxx:184
CbmStsDigi::GetSide
Int_t GetSide() const
Definition CbmStsDigi.h:79
CbmStsDigi::GetAdc
Int_t GetAdc() const
Definition CbmStsDigi.h:87
CbmStsDigi::GetSectorNr
Int_t GetSectorNr() const
Definition CbmStsDigi.h:75
CbmStsDigi::GetChannelNr
Int_t GetChannelNr() const
Definition CbmStsDigi.h:83
CbmStsDigi::GetStationNr
Int_t GetStationNr() const
Definition CbmStsDigi.h:72
CbmStsSector::GetType
Int_t GetType() const
Definition CbmStsSector.h:93
CbmStsSector::GetSectorNr
Int_t GetSectorNr() const
Definition CbmStsSector.h:76
CbmStsStation::GetNSectors
Int_t GetNSectors() const
Definition CbmStsStation.h:62
CbmStsStation::GetStationNr
Int_t GetStationNr() const
Definition CbmStsStation.h:56
CbmStsStation::GetSector
CbmStsSector * GetSector(Int_t iSector)
Definition CbmStsStation.h:65
nlohmann::detail::cbor_tag_handler_t::error
@ error
throw a parse_error exception in case of a tag
run
-clang-format