BmnRoot
Loading...
Searching...
No Matches
CbmStsMatchHits.cxx
Go to the documentation of this file.
1//* $Id: */
2
3// -------------------------------------------------------------------------
4// ----- CbmStsMatchHits source file -----
5// ----- Created 01/07/2008 by R. Karabowicz -----
6// -------------------------------------------------------------------------
7
8
9// --- Includes from ROOT
10#include "TClonesArray.h"
11
12// --- Includes from base
13#include "FairRootManager.h"
14#include "FairRunAna.h"
15#include "FairRuntimeDb.h"
16
17// --- Includes from STS
18#include "CbmGeoStsPar.h"
19#include "FairGeoPassivePar.h"
20#include "CbmStsDigiPar.h"
21#include "CbmStsDigiScheme.h"
22#include "CbmStsDigiMatch.h"
23#include "CbmStsCluster.h"
24#include "CbmStsDigi.h"
25#include "CbmStsHit.h"
26#include "CbmStsMatchHits.h"
27#include "CbmStsPoint.h"
28#include "CbmStsSensor.h" //AZ
29#include "CbmStsSector.h"
30#include "CbmStsStation.h"
31#include "BmnSiliconPoint.h" //AZ-280322
32#include "FairGeoVector.h"
33#include "FairGeoNode.h"
34
35#include "TGeoManager.h" //AZ
36#include "TMath.h"
37
38#include <iostream>
39#include <iomanip>
40#include <map>
41
42using std::cout;
43using std::flush;
44using std::endl;
45using std::setw;
46using std::left;
47using std::right;
48using std::fixed;
49using std::setprecision;
50using std::map;
51
52// ----- Default constructor -------------------------------------------
53CbmStsMatchHits::CbmStsMatchHits()
54: FairTask("STSMatchHits", 1),
55 fGeoPar(NULL),
56 fDigiPar(NULL),
57 fDigiScheme(NULL),
58 fPoints(NULL),
59 fPointsSi(nullptr), //AZ-280322
60 fDigis(NULL),
61 fDigiMatches(NULL),
62 fHits(NULL),
63 fTimer(),
64 fNEvents(0),
65 fNEventsFailed(0),
66 fTime(0.),
67 fNHits(0.),
68 fNMatched(0.),
69 fNDistant(0.),
70 fNBackgrd(0.),
71 fPassGeo(NULL),
72 fTargetPos(0.,0.,0.),
73 fNStations(0),
74 fRealistic(kTRUE),
75 fCandMap(),
76 fIter()
77{
78 //AZ fDigiScheme = new CbmStsDigiScheme();
79 fDigiScheme = CbmStsDigiScheme::Instance();
80}
81// -------------------------------------------------------------------------
82
83
84// ----- Standard constructor ------------------------------------------
85CbmStsMatchHits::CbmStsMatchHits(Int_t iVerbose)
86: FairTask("STSMatchHits", iVerbose),
87 fGeoPar(NULL),
88 fDigiPar(NULL),
89 fDigiScheme(NULL),
90 fPoints(NULL),
91 fPointsSi(nullptr), //AZ-280322
92 fDigis(NULL),
93 fDigiMatches(NULL),
94 fHits(NULL),
95 fTimer(),
96 fNEvents(0),
97 fNEventsFailed(0),
98 fTime(0.),
99 fNHits(0.),
100 fNMatched(0.),
101 fNDistant(0.),
102 fNBackgrd(0.),
103 fPassGeo(NULL),
104 fTargetPos(0.,0.,0.),
105 fNStations(0),
106 fRealistic(kTRUE),
107 fCandMap(),
108 fIter()
109{
110 //AZ fDigiScheme = new CbmStsDigiScheme();
111 fDigiScheme = CbmStsDigiScheme::Instance();
112}
113// -------------------------------------------------------------------------
114
115
116// ----- Constructor with name -----------------------------------------
117CbmStsMatchHits::CbmStsMatchHits(const char* name, Int_t iVerbose)
118: FairTask(name, iVerbose),
119 fGeoPar(NULL),
120 fDigiPar(NULL),
121 fDigiScheme(NULL),
122 fPoints(NULL),
123 fPointsSi(nullptr), //AZ-280322
124 fDigis(NULL),
125 fDigiMatches(NULL),
126 fHits(NULL),
127 fTimer(),
128 fNEvents(0),
129 fNEventsFailed(0),
130 fTime(0.),
131 fNHits(0.),
132 fNMatched(0.),
133 fNDistant(0.),
134 fNBackgrd(0.),
135 fPassGeo(NULL),
136 fTargetPos(0.,0.,0.),
137 fNStations(0),
138 fRealistic(kTRUE),
139 fCandMap(),
140 fIter()
141{
142 //AZ fDigiScheme = new CbmStsDigiScheme();
143 fDigiScheme = CbmStsDigiScheme::Instance();
144}
145// -------------------------------------------------------------------------
146
147
148
149// ----- Destructor ----------------------------------------------------
150CbmStsMatchHits::~CbmStsMatchHits() {
151 if ( fPassGeo ) delete fPassGeo;
152 if ( fGeoPar ) delete fGeoPar;
153 if ( fDigiPar ) delete fDigiPar;
154 //AZ if ( fDigiScheme ) delete fDigiScheme;
155}
156// -------------------------------------------------------------------------
157
158
159
160// ----- Public method Exec --------------------------------------------
161void CbmStsMatchHits::Exec(Option_t* opt) {
162
163 if ( fRealistic ) {
164 ExecReal(opt);
165 return;
166 }
167
168 // Timer
169 fTimer.Start();
170 Bool_t warn = kFALSE;
171
172 // Counters
173 Int_t nHits = fHits->GetEntriesFast();
174 Int_t nNoDigi = 0;
175 Int_t nBackgrd = 0;
176 Int_t nDistant = 0;
177 Int_t nMatched = 0;
178
179 // Loop over all StsHits
180 for (Int_t iHit=0; iHit<nHits; iHit++) {
181 CbmStsHit* hit = (CbmStsHit*) fHits->At(iHit);
182 if ( ! hit ) {
183 cout << "-W- " << GetName() << "::Exec: Empty hit at index "
184 << iHit << endl;
185 warn = kTRUE;
186 continue;
187 }
188
189 // Determine sector type and channel numbers
190 Int_t iStation = hit->GetStationNr();
191 Int_t iSector = hit->GetSectorNr();
192 CbmStsStation* station = fDigiScheme->GetStationByNr(iStation);
193 CbmStsSector* sector = fDigiScheme->GetSector(iStation, iSector);
194 Int_t iType = sector->GetType();
195 CbmStsDigiMatch* dMatchF = NULL;
196 CbmStsDigiMatch* dMatchB = NULL;
197
198
199 // Hit coordinates and errors
200 Double_t xH = hit->GetX();
201 Double_t yH = hit->GetY();
202 Double_t dX = hit->GetDx();
203 Double_t dY = hit->GetDy();
204
205 // Get front side DigiMatch corresponding to hit
206 Int_t iMatchF = (Int_t)hit->GetDigi(0);
207 if ( iMatchF >= 0 )
208 dMatchF = (CbmStsDigiMatch*) fDigiMatches->At(iMatchF);
209 if ( ! dMatchF ) {
210 cout << "-E- " << GetName() << "::Exec: "
211 << "No DigiMatchF for hit " << iHit << endl;
212 hit->SetRefIndex(-1);
213 nNoDigi++;
214 warn = kTRUE;
215 continue;
216 }
217
218 // Get back side DigiMatch of hit (for strip sensors only)
219 if ( iType != 1 ) {
220 Int_t iMatchB = (Int_t)hit->GetDigi(1);
221 if ( iMatchB >= 0 )
222 dMatchB = (CbmStsDigiMatch*) fDigiMatches->At(iMatchB);
223 if ( ! dMatchB ) {
224 cout << "-E- " << GetName() << "::Exec: "
225 << "No DigiMatchB for hit " << iHit << endl;
226 hit->SetRefIndex(-1);
227 nNoDigi++;
228 warn = kTRUE;
229 continue;
230 }
231 }
232
233 // Map candidate points to distance to hit
234 fCandMap.clear();
235 Int_t nPointsF = 0;
236 Int_t nPointsB = 0;
237
238 if ( iType == 1 ) {
239 // Case pixel: Candidates are all points of the digi.
240 for (Int_t jMatchF=0; jMatchF<3; jMatchF++) {
241 Int_t iPointF = dMatchF->GetRefIndex(jMatchF);
242 if ( iPointF < 0 ) continue;
243 nPointsF++;
244 // Calculate distance to hit
245 CbmStsPoint* point = (CbmStsPoint*) fPoints->At(iPointF);
246 if ( ! point ) {
247 cout << "-E- " << GetName() << "::Exec: "
248 << "No StsPoint (" << iPointF << ") for pixel hit "
249 << iHit << endl;
250 warn = kTRUE;
251 continue;
252 }
253 // point coordinates in station midplane
254 Double_t xP = point->GetX(station->GetZ());
255 Double_t yP = point->GetY(station->GetZ());
256 Double_t dist = TMath::Sqrt( (xP-xH)*(xP-xH) + (yP-yH)*(yP-yH) );
257 // cout << "candidate dist = " << dist << endl;
258 fCandMap[dist] = iPointF;
259 } // front digi loop
260 } // pixel sensor
261
262 else if ( iType == 2 || iType == 3 ) {
263 // Case strip: Candidates are points corresponding to front
264 // and back side digis
265 for ( Int_t jMatchF=0; jMatchF<3; jMatchF++) {
266 Int_t iPointF = dMatchF->GetRefIndex(jMatchF);
267 if ( iPointF < 0 ) continue;
268 // cout << " got pointf index = " << iPointF << endl;
269 nPointsF++;
270 for ( Int_t iMatchB=0; iMatchB<3; iMatchB++) {
271 Int_t iPointB = dMatchB->GetRefIndex(iMatchB);
272 if ( iPointB < 0 ) continue;
273 // cout << " got pointb index = " << iPointB << endl;
274 if ( iMatchF == 0 ) nPointsB++;
275 if ( iPointB != iPointF ) continue; // chance combination
276 // Calculate distance to hit
277 CbmStsPoint* point = (CbmStsPoint*) fPoints->At(iPointF);
278 if ( ! point ) {
279 cout << "-E- " << GetName() << "::Exec: "
280 << "No StsPoint (" << iPointF << ") for strip hit "
281 << iHit << endl;
282 warn = kTRUE;
283 continue;
284 }
285 // point coordinates in station midplane
286 Double_t xP = point->GetX(station->GetZ());
287 Double_t yP = point->GetY(station->GetZ());
288 Double_t dist = TMath::Sqrt( (xP-xH)*(xP-xH) +
289 (yP-yH)*(yP-yH) );
290 fCandMap[dist] = iPointF;
291 } // back digi loop
292 } // front digi loop
293 } // strip sensor
294
295 else {
296 // Unknown sensor type
297 cout << "-E- " << GetName() << "::Exec: Unknwon sensor type "
298 << iType << endl;
299 Fatal("Exec", "Unknwon sensor type");
300 }
301
302 if ( fVerbose > 1 ) cout << "-I- " << GetName() << ": Hit "
303 << iHit << ", type " << iType
304 << ", points " << nPointsF << " / "
305 << nPointsB << ", candidates "
306 << fCandMap.size();
307
308
309 // Check for at least one candidate point. Else background hit.
310 // This can happen for noise digis or for fake combinations
311 // of strip digis.
312 if ( fCandMap.empty() ) {
313 hit->SetRefIndex(-1);
314 if (fVerbose>1) cout << ", background " << endl;
315 nBackgrd++;
316 continue;
317 }
318
319 // Select closest point from candidates.
320 Double_t distMin = 999999.;
321 Int_t iPoint = -1;
322 for (fIter=fCandMap.begin(); fIter!=fCandMap.end(); fIter++) {
323 if ( (*fIter).first < distMin ) {
324 distMin = (*fIter).first;
325 iPoint = (*fIter).second;
326 }
327 }
328 if ( iPoint < 0 ) {
329 cout << "-E- " << GetName() << "::Exec: "
330 << "No closest point found in candidate map!" << endl;
331 Fatal("Exec", "No closest point");
332 }
333 CbmStsPoint* point = (CbmStsPoint*) fPoints->At(iPoint);
334 if (fVerbose>1) cout << ", matched to " << iPoint << ", distance "
335 << distMin << " cm";
336
337 // Check whether closest point is more than 5 sigma away from hit.
338 // This should not happen in case of pixel or strip OSU sensors,
339 // but can happen for strip GSI sensors due to multiple
340 // crossings of front and back side strips.
341 Double_t xP = point->GetX(station->GetZ());
342 Double_t yP = point->GetY(station->GetZ());
343 if ( TMath::Abs(xP-xH) > 5. * dX ||
344 TMath::Abs(yP-yH) > 5. * dY ) {
345 hit->SetRefIndex(-2);
346 nDistant++;
347 if (fVerbose>1) cout << ", distant" << endl;
348 if ( iType == 1 || iType == 2) {
349 cout << "-W- " << fName << "::Exec: "
350 << "Distant hit in pixel / strip MSU station" << endl;
351 cout << "Hit " << iHit << " coordinates " << xH
352 << ", " << yH << endl;
353 cout << "Matched point " << iPoint << ", distance "
354 << distMin << endl;
355 cout << "Distance y " << TMath::Abs(xP-xH) << " Error " << dX << endl;
356 cout << "Distance y " << TMath::Abs(yP-yH) << " Error " << dY << endl;
357 }
358
359
360 continue;
361 }
362
363 // Match closest StsPoint to hit
364 hit->SetRefIndex(iPoint);
365 nMatched++;
366 if (fVerbose>1) cout << ", good match" << endl;
367
368 // Check whether hit and point are in the same station.
369 // Else something has really gone wrong.
370 if ( TMath::Abs(hit->GetZ() - point->GetZ()) > 1. ) {
371 cout << "-E- " << GetName() << "::Exec: Hit " << iHit
372 << " is at z = " << hit->GetZ() << " cm, but matched point "
373 << iPoint << " is at z = " << point->GetZ() << "cm " << endl;
374 Fatal("Exec", "Different stations for hit and point");
375 }
376
377 } // Loop over StsHits
378
379 // Event statistics
380 fTimer.Stop();
381 if ( warn ) cout << "- ";
382 else cout << "+ ";
383 cout << setw(15) << left << fName << ": " << setprecision(4) << setw(8)
384 << fixed << right << fTimer.RealTime()
385 << " s, hits " << nHits << ", matched " << nMatched << ", distant "
386 << nDistant << ", background " << setw(6) << nBackgrd << endl;
387
388 // Run statistics
389 if ( warn) fNEventsFailed++;
390 else {
391 fNEvents++;
392 fTime += fTimer.RealTime();
393 fNHits += Double_t(nHits);
394 fNMatched += Double_t(nMatched);
395 fNDistant += Double_t(nDistant);
396 fNBackgrd += Double_t(nBackgrd);
397 }
398
399}
400
401// -------------------------------------------------------------------------
402
403// ----- Public method ExecReal --------------------------------------------
404//AZ
405void CbmStsMatchHits::ExecReal(Option_t* opt) {
406
407 // Timer
408 fTimer.Start();
409 Bool_t warn = kFALSE;
410
411 // Counters
412 Int_t nHits = fHits->GetEntriesFast();
413 //Int_t nNoDigi = 0;
414 Int_t nBackgrd = 0;
415 Int_t nDistant = 0;
416 Int_t nMatched = 0;
417
418 Int_t nofStsHits = fHits->GetEntriesFast();
419 Int_t nofStsPoints = fPoints->GetEntriesFast();
420 int nofSiPoints = 0; //AZ-280322
421 if (fPointsSi) nofSiPoints = fPointsSi->GetEntriesFast(); //AZ-280322
422 //cout << "there are " << nofStsPoints << " points and " << nofStsHits << " hits." << endl;
423 Int_t hitStationLimits[2][100];
424
425 // Workaround. TODO.
426 //fNStations = 8;
427 fNStations = 20; //AZ
428
429 for ( Int_t ist = 0 ; ist < fNStations ; ist++ ) {
430 hitStationLimits[0][ist] = -1;
431 hitStationLimits[1][ist] = -1;
432 }
433
434 // check for limits of hit indices on different stations...
435 for ( Int_t ihit = 0 ; ihit < nofStsHits ; ihit++ ) {
436 CbmStsHit *stsHit = (CbmStsHit*)fHits->At(ihit);
437 stsHit->SetRefIndex(-1);
438 if ( hitStationLimits[0][stsHit->GetStationNr()-1] == -1 )
439 hitStationLimits[0][stsHit->GetStationNr()-1] = ihit;
440 CbmStsHit *stsHitBack = (CbmStsHit*)fHits->At(nofStsHits-ihit-1);
441 if ( hitStationLimits[1][stsHitBack->GetStationNr()-1] == -1 ) {
442 hitStationLimits[1][stsHitBack->GetStationNr()-1] = nofStsHits-ihit;
443 }
444 }
445 //for ( Int_t istat = 0 ; istat < fNStations ; istat++ )
446 // cout << "station " << istat << " hits from " << hitStationLimits[0][istat] << " to " << hitStationLimits[1][istat] << endl;
447
448 int nTotPoints = nofSiPoints + nofStsPoints; //AZ-280322
449 CbmStsPoint *stsPoint = nullptr; //AZ-280322
450
451 //AZ-280322 for ( Int_t ipnt = 0 ; ipnt < nofStsPoints ; ipnt++ ) {
452 for ( Int_t iii = 0 ; iii < nTotPoints ; iii++ ) { //AZ-280322
453 //AZ-280322 CbmStsPoint *stsPoint = (CbmStsPoint*)fPoints->At(ipnt);
454 int ipnt = (iii < nofSiPoints) ? iii : iii - nofSiPoints;
455 if (iii < nofSiPoints) stsPoint = (CbmStsPoint*)fPointsSi->At(ipnt); //AZ-280322
456 else stsPoint = (CbmStsPoint*)fPoints->At(ipnt); //AZ-280322
457
458 // This is a workaround. Has to be worked out when new digi scheme is implemented. TODO.
459 //AZ Double_t zPoint = stsPoint->GetZIn();
460 //AZ Int_t stationNr = TMath::Nint((zPoint-30.)/10.);
461 // AZ - Get station number
462 Double_t xin = stsPoint->GetXIn();
463 Double_t yin = stsPoint->GetYIn();
464 Double_t zin = stsPoint->GetZIn();
465 gGeoManager->FindNode(xin,yin,zin);
466 TGeoNode* curNode = gGeoManager->GetCurrentNode();
467 CbmStsSensor* sensor = NULL;
468 //if ( fDigiScheme->IsNewGeometry() ) {
469 TString curPath = fDigiScheme->GetCurrentPath();
470 sensor = fDigiScheme->GetSensorByName(curPath);
471 if ( !sensor ) {
472 //AZ - 160920: Check one level above (for double-sensor modules)
473 gGeoManager->CdUp();
474 curPath = fDigiScheme->GetCurrentPath();
475 sensor = fDigiScheme->GetSensorByName(curPath);
476 //
477 if (!sensor) {
478 std::cerr << "-E- " << fName << "::ExecReal:: sensor " << curNode->GetName()
479 << " not found in digi scheme!" << endl;
480 continue;
481 }
482 }
483 Int_t stationNr = sensor->GetStationNr() - 1;
484 //}
485 //AZ
486
487 // This does not work with new geometries since the StsGeoPar are not filled.
488 // Int_t stationNr = fStationNrFromMcId[stsPoint->GetDetectorID()]];
489
490 Int_t startHit = hitStationLimits[0][stationNr];
491 Int_t finalHit = hitStationLimits[1][stationNr];
492
493 if ( startHit == -1 && finalHit == -1 ) continue;
494 if ( startHit < 0 || finalHit < 0 || startHit > fHits->GetEntriesFast() || finalHit > fHits->GetEntriesFast()) continue;
495
496 Int_t NofMatched = 0;
497 //AZ Double_t GoodHit[100];
498 //AZ Int_t GoodHitIndex[nofStsHits];
499 Double_t GoodHit[50000];
500 Int_t GoodHitIndex[50000];
501 for ( Int_t ihit = startHit ; ihit < finalHit ; ihit++ ) {
502 CbmStsHit *stsHit= (CbmStsHit*)fHits->At(ihit);
503 if (NULL == stsHit) {
504 std::cout << "CbmStsMatchHits::ExecReal: CbmStsHit NULL pointer: ihit=" << ihit
505 << " " << "startHit=" << startHit << " finalHit=" << finalHit << std::endl;
506 continue;
507 }
508
509 // Check one side
510 Int_t iL = 0, nMatch = 0;
511 FairLink link = stsHit->GetLink(iL);
512 CbmStsCluster *stsCluster = dynamic_cast<CbmStsCluster*> (fClusters->UncheckedAt(link.GetIndex()));
513 Int_t nLinks2 = stsCluster->GetNLinks();
514 for (Int_t iL2 = 0; iL2 < nLinks2; ++iL2) {
515 FairLink link2 = stsCluster->GetLink(iL2);
516 //CbmStsDigi *stsDigi = dynamic_cast<CbmStsDigi*> (fDigis->UncheckedAt(link2.GetIndex()));
517 CbmStsDigiMatch *stsDigiMatch = dynamic_cast<CbmStsDigiMatch*> (fDigiMatches->UncheckedAt(link2.GetIndex()));
518 //const Int_t nLinks3 = stsDigi->GetNLinks();
519 const Int_t nLinks3 = 3;
520 for (Int_t iL3 = 0; iL3 < nLinks3; ++iL3) {
521 //FairLink link3 = stsDigi->GetLink(iL3);
522 //Int_t stsPointId = link3.GetIndex();
523 Int_t stsPointId = stsDigiMatch->GetRefIndex(iL3);
524 if (stsPointId == ipnt) { ++nMatch; break; }
525 else if (stsPointId < 0) break;
526 } // for mcPoint
527 if (nMatch) break;
528 } // for digi
529
530 // Check the other side
531 iL = 1;
532 link = stsHit->GetLink(iL);
533 stsCluster = dynamic_cast<CbmStsCluster*> (fClusters->UncheckedAt(link.GetIndex()));
534 nLinks2 = stsCluster->GetNLinks();
535 for (Int_t iL2 = 0; iL2 < nLinks2; ++iL2) {
536 FairLink link2 = stsCluster->GetLink(iL2);
537 //CbmStsDigi *stsDigi = dynamic_cast<CbmStsDigi*> (fDigis->UncheckedAt(link2.GetIndex()));
538 CbmStsDigiMatch *stsDigiMatch = dynamic_cast<CbmStsDigiMatch*> (fDigiMatches->UncheckedAt(link2.GetIndex()));
539 //const Int_t nLinks3 = stsDigi->GetNLinks();
540 const Int_t nLinks3 = 3;
541 for (Int_t iL3 = 0; iL3 < nLinks3; ++iL3) {
542 //FairLink link3 = stsDigi->GetLink(iL3);
543 //Int_t stsPointId = link3.GetIndex();
544 Int_t stsPointId = stsDigiMatch->GetRefIndex(iL3);
545 if (stsPointId == ipnt) { ++nMatch; break; }
546 else if (stsPointId < 0) break;
547 } // for mcPoint
548 if (nMatch == 2) break;
549 } // for digi
550
551 if (nMatch == 2) {
552 Double_t xH = stsHit->GetX();
553 Double_t yH = stsHit->GetY();
554 //AZ-280322 Double_t xP = stsPoint->GetX(stsHit->GetZ());
555 //AZ-280322 Double_t yP = stsPoint->GetY(stsHit->GetZ());
556 Double_t xP = stsPoint->GetXIn(); //AZ-280322
557 Double_t yP = stsPoint->GetYIn(); //AZ-280322
558 Double_t dist = TMath::Sqrt( (xP-xH)*(xP-xH) + (yP-yH)*(yP-yH) );
559 nMatched++;
560 NofMatched++;
561 GoodHitIndex[NofMatched]=ihit;
562 GoodHit[NofMatched]=dist;
563 }
564
565 } // for ( Int_t ihit = startHit ; ihit < finalHit ;
566
567 Int_t k=NofMatched;
568 if (NofMatched>1){
569 for (Int_t i=NofMatched-1; i>0; i--) {
570 if (GoodHit[i]<GoodHit[k]) k=i;
571 }
572 }
573 if (NofMatched>0) {
574 CbmStsHit *stsHit= (CbmStsHit*)fHits->At(GoodHitIndex[k]);
575 stsHit->SetRefIndex(ipnt);
576 }
577 } // for ( Int_t iii = 0 ; iii < nTotPoints ;
578
579 // Event statistics
580 fTimer.Stop();
581 if ( warn ) cout << "- ";
582 else cout << "+ ";
583 cout << setw(15) << left << fName << ": " << setprecision(4) << setw(8)
584 << fixed << right << fTimer.RealTime()
585 << " s, hits " << nHits << ", matched " << nMatched << ", distant "
586 << nDistant << ", background " << setw(6) << nBackgrd << endl;
587
588 // Run statistics
589 if ( warn) fNEventsFailed++;
590 else {
591 fNEvents++;
592 fTime += fTimer.RealTime();
593 fNHits += Double_t(nHits);
594 fNMatched += Double_t(nMatched);
595 fNDistant += Double_t(nDistant);
596 fNBackgrd += Double_t(nBackgrd);
597 }
598
599}
600// -------------------------------------------------------------------------
601
602// ----- Public method ExecReal --------------------------------------------
603/*
604void CbmStsMatchHits::ExecReal(Option_t* opt) {
605
606 // Timer
607 fTimer.Start();
608 Bool_t warn = kFALSE;
609
610 // Counters
611 Int_t nHits = fHits->GetEntriesFast();
612 Int_t nNoDigi = 0;
613 Int_t nBackgrd = 0;
614 Int_t nDistant = 0;
615 Int_t nMatched = 0;
616
617 Int_t nofStsHits = fHits->GetEntriesFast();
618 Int_t nofStsPoints = fPoints->GetEntriesFast();
619 //cout << "there are " << nofStsPoints << " points and " << nofStsHits << " hits." << endl;
620 Int_t hitStationLimits[2][100];
621
622 // Workaround. TODO.
623 fNStations = 8;
624
625 for ( Int_t ist = 0 ; ist < fNStations ; ist++ ) {
626 hitStationLimits[0][ist] = -1;
627 hitStationLimits[1][ist] = -1;
628 }
629
630 // check for limits of hit indices on different stations...
631 for ( Int_t ihit = 0 ; ihit < nofStsHits ; ihit++ ) {
632 CbmStsHit *stsHit = (CbmStsHit*)fHits->At(ihit);
633 stsHit->SetRefIndex(-1);
634 if ( hitStationLimits[0][stsHit->GetStationNr()-1] == -1 )
635 hitStationLimits[0][stsHit->GetStationNr()-1] = ihit;
636 CbmStsHit *stsHitBack = (CbmStsHit*)fHits->At(nofStsHits-ihit-1);
637 if ( hitStationLimits[1][stsHitBack->GetStationNr()-1] == -1 ) {
638 hitStationLimits[1][stsHitBack->GetStationNr()-1] = nofStsHits-ihit;
639 }
640 }
641 //for ( Int_t istat = 0 ; istat < fNStations ; istat++ )
642 // cout << "station " << istat << " hits from " << hitStationLimits[0][istat] << " to " << hitStationLimits[1][istat] << endl;
643
644 for ( Int_t ipnt = 0 ; ipnt < nofStsPoints ; ipnt++ ) {
645 CbmStsPoint *stsPoint = (CbmStsPoint*)fPoints->At(ipnt);
646
647 // This is a workaround. Has to be worked out when new digi scheme is implemented. TODO.
648 Double_t zPoint = stsPoint->GetZIn();
649 Int_t stationNr = TMath::Nint((zPoint-30.)/10.);
650
651 // This does not work with new geometries since the StsGeoPar are not filled.
652 // Int_t stationNr = fStationNrFromMcId[stsPoint->GetDetectorID()]];
653
654 Int_t startHit = hitStationLimits[0][stationNr];
655 Int_t finalHit = hitStationLimits[1][stationNr];
656
657 if ( startHit == -1 && finalHit == -1 ) continue;
658 if ( startHit < 0 || finalHit < 0 || startHit > fHits->GetEntriesFast() || finalHit > fHits->GetEntriesFast()) continue;
659
660
661 Int_t NofMatched = 0;
662 Double_t GoodHit[100];
663 Int_t GoodHitIndex[nofStsHits];
664 for ( Int_t ihit = startHit ; ihit < finalHit ; ihit++ ) {
665 CbmStsHit *stsHit= (CbmStsHit*)fHits->At(ihit);
666 if (NULL == stsHit) {
667 std::cout << "CbmStsMatchHits::ExecReal: CbmStsHit NULL pointer: ihit=" << ihit
668 << " " << "startHit=" << startHit << " finalHit=" << finalHit << std::endl;
669 continue;
670 }
671 Double_t xH = stsHit->GetX();
672 Double_t yH = stsHit->GetY();
673 Double_t xP = stsPoint->GetX(stsHit->GetZ());
674 Double_t yP = stsPoint->GetY(stsHit->GetZ());
675 if ( ( TMath::Abs(stsHit->GetX()-stsPoint->GetX(stsHit->GetZ())) < .01 ) &&
676 ( TMath::Abs(stsHit->GetY()-stsPoint->GetY(stsHit->GetZ())) < .04 ) ) {
677 Double_t dist = TMath::Sqrt( (xP-xH)*(xP-xH) + (yP-yH)*(yP-yH) );
678 // cout << "matching "
679 // << "X: " << stsHit->GetX() << " - " << stsPoint->GetX(stsHit->GetZ())
680 // << "Y: " << stsHit->GetY() << " - " << stsPoint->GetY(stsHit->GetZ()) << endl;
681// cout << "candidate dist = " << dist << endl;
682
683 // cout << "setting ref index = " << stsHit->GetRefIndex() << " (max sts pnt = " << nofStsPoints << ")" << endl;
684 nMatched++;
685 NofMatched++;
686 GoodHitIndex[NofMatched]=ihit;
687 GoodHit[NofMatched]=dist;
688// cout <<" NofMatched = "<<NofMatched<<" dist = "<<dist<<flush;
689 }
690
691 }
692
693 Int_t k=NofMatched;
694 if (NofMatched>1){
695 for (Int_t i=NofMatched-1; i>0; i--) {
696
697 if (GoodHit[i]<GoodHit[k]) {
698 k=i;
699 }
700 }
701 }
702 if (NofMatched>0) {
703 CbmStsHit *stsHit= (CbmStsHit*)fHits->At(GoodHitIndex[k]);
704 stsHit->SetRefIndex(ipnt);
705 }
706 }
707
708 // Event statistics
709 fTimer.Stop();
710 if ( warn ) cout << "- ";
711 else cout << "+ ";
712 cout << setw(15) << left << fName << ": " << setprecision(4) << setw(8)
713 << fixed << right << fTimer.RealTime()
714 << " s, hits " << nHits << ", matched " << nMatched << ", distant "
715 << nDistant << ", background " << setw(6) << nBackgrd << endl;
716
717 // Run statistics
718 if ( warn) fNEventsFailed++;
719 else {
720 fNEvents++;
721 fTime += fTimer.RealTime();
722 fNHits += Double_t(nHits);
723 fNMatched += Double_t(nMatched);
724 fNDistant += Double_t(nDistant);
725 fNBackgrd += Double_t(nBackgrd);
726 }
727
728}
729*/
730// -------------------------------------------------------------------------
731
732// ----- Private method SetParContainers -------------------------------
733void CbmStsMatchHits::SetParContainers() {
734
735 // Get run and runtime database
736 FairRunAna* run = FairRunAna::Instance();
737 if ( ! run ) Fatal("SetParContainers", "No analysis run");
738 FairRuntimeDb* db = run->GetRuntimeDb();
739 if ( ! db ) Fatal("SetParContainers", "No runtime database");
740
741 // Get STS geometry and digitisation parameter container
742 fPassGeo = (FairGeoPassivePar*) db->getContainer("FairGeoPassivePar");
743 fGeoPar = (CbmGeoStsPar*) db->getContainer("CbmGeoStsPar");
744 fDigiPar = (CbmStsDigiPar*) db->getContainer("CbmStsDigiPar");
745}
746// -------------------------------------------------------------------------
747
748
749
750// ----- Private method Init -------------------------------------------
751InitStatus CbmStsMatchHits::Init() {
752
753 // Reset counters
754 fNEvents = fNEventsFailed = 0;
755 fTime = fNHits = fNMatched = fNDistant = fNBackgrd = 0.;
756
757 // Get input arrays
758 FairRootManager* ioman = FairRootManager::Instance();
759 if ( ! ioman ) Fatal("Init", "No FairRootManager");
760 fPoints = (TClonesArray*) ioman->GetObject("StsPoint");
761 if ( ! fPoints ) {
762 //cout << "-E- " << GetName() << "::Init: No STSPoint array!" << endl;
763 return kERROR;
764 }
765 fPointsSi = (TClonesArray*) ioman->GetObject("SiliconPoint"); //AZ-289322
766 fDigis = (TClonesArray*) ioman->GetObject("StsDigi");
767 if ( ! fDigis ) {
768 cout << "-E- " << GetName() << "::Init: No StsDigi array!" << endl;
769 return kERROR;
770 }
771 fDigiMatches = (TClonesArray*) ioman->GetObject("StsDigiMatch");
772 if ( ! fDigiMatches ) {
773 cout << "-E- " << GetName() << "::Init: No StsDigiMatch array!"
774 << endl;
775 return kERROR;
776 }
777 fHits = (TClonesArray*) ioman->GetObject("StsHit");
778 if ( ! fHits ) {
779 cout << "-E- " << GetName() << "::Init: No StsHit array!" << endl;
780 return kERROR;
781 }
782 //AZ
783 fClusters = (TClonesArray*) ioman->GetObject("StsCluster");
784 if ( ! fClusters ) {
785 cout << "-E- " << GetName() << "::Init: No StsCluster array!" << endl;
786 return kERROR;
787 }
788 //AZ
789
790 InitStatus geoStatus = GetGeometry();
791 if ( geoStatus != kSUCCESS ) {
792 cout << "-E- " << GetName() << "::Init: Error in reading geometry!"
793 << endl;
794 return geoStatus;
795 }
796
797 // Build digitisation scheme
798
799 if ( fDigiScheme->Init(fGeoPar, fDigiPar) ) {
800
801 if (fVerbose == 1 || fVerbose == 2) fDigiScheme->Print(kFALSE);
802 else if (fVerbose > 2) fDigiScheme->Print(kTRUE);
803 cout << "-I- " << fName << "::Init: "
804 << "STS digitisation scheme succesfully initialised" << endl;
805 cout << " Stations: " << fDigiScheme->GetNStations()
806 << ", Sectors: " << fDigiScheme->GetNSectors() << ", Channels: "
807 << fDigiScheme->GetNChannels() << endl;
808 return kSUCCESS;
809 }
810
811 return kERROR;
812
813}
814// -------------------------------------------------------------------------
815
816
817
818// ----- Private method ReInit -----------------------------------------
819InitStatus CbmStsMatchHits::ReInit() {
820
821 // Clear digitisation scheme
822 fDigiScheme->Clear();
823
824 // Build new digitisation scheme
825 if ( fDigiScheme->Init(fGeoPar, fDigiPar) ) return kSUCCESS;
826
827 InitStatus geoStatus = GetGeometry();
828 if ( geoStatus != kSUCCESS ) {
829 cout << "-E- " << GetName() << "::Init: Error in reading geometry!"
830 << endl;
831 return geoStatus;
832 }
833
834 return kERROR;
835
836}
837// -------------------------------------------------------------------------
838
839// ----- Private method GetGeometry ------------------------------------
840InitStatus CbmStsMatchHits::GetGeometry() {
841cout << "-W- " << endl;
842 // Get target geometry
843 if ( ! fPassGeo ) {
844 cout << "-W- " << GetName() << "::GetGeometry: No passive geometry!"
845 <<endl;
846 fTargetPos.SetXYZ(0., 0., 0.);
847 return kERROR;
848 }
849 TObjArray* passNodes = fPassGeo->GetGeoPassiveNodes();
850 if ( ! passNodes ) {
851 cout << "-W- " << GetName() << "::GetGeometry: No passive node array"
852 << endl;
853 fTargetPos.SetXYZ(0., 0., 0.);
854 return kERROR;
855 }
856 /* FairGeoNode* target = (FairGeoNode*) passNodes->FindObject("targ");
857 if ( ! target ) {
858 cout << "-E- " << GetName() << "::GetGeometry: No target node"
859 << endl;
860 fTargetPos.SetXYZ(0., 0., 0.);
861 return kERROR;
862 }
863 FairGeoVector targetPos = target->getLabTransform()->getTranslation();
864 FairGeoVector centerPos = target->getCenterPosition().getTranslation();
865 Double_t targetX = targetPos.X() + centerPos.X();
866 Double_t targetY = targetPos.Y() + centerPos.Y();
867 Double_t targetZ = targetPos.Z() + centerPos.Z();
868 fTargetPos.SetXYZ(targetX, targetY, targetZ);*/
869
870 // Get STS geometry
871 if ( ! fGeoPar ) {
872 cout << "-W- " << GetName() << "::GetGeometry: No passive geometry!"
873 <<endl;
874 fNStations = 0;
875 return kERROR;
876 }
877 TObjArray* stsNodes = fGeoPar->GetGeoSensitiveNodes();
878 if ( ! stsNodes ) {
879 cout << "-E- " << GetName() << "::GetGeometry: No STS node array"
880 << endl;
881 fNStations = 0;
882 return kERROR;
883 }
884 Int_t tempNofStations = stsNodes->GetEntries();
885 cout << "Nodes in STS: " << tempNofStations << endl;
886
887 cout << "There are " << tempNofStations << " nodes" << (tempNofStations > 10 ? "!!!" : "" ) << endl;
888
889 TString geoNodeName;
890 fNStations = 0;
891 TString stationNames[1000];
892 for ( Int_t ist = 0 ; ist < tempNofStations ; ist++ ) {
893 FairGeoNode* stsNode = (FairGeoNode*)stsNodes->At(ist);
894 if ( ! stsNode ) {
895 cout << "-W- CbmStsDigiScheme::Init: station#" << ist
896 << " not found among sensitive nodes " << endl;
897 continue;
898 }
899 geoNodeName = stsNode->getName();
900 // TArrayD* params = stsNode->getParameters();
901
902 Bool_t stationKnown = kFALSE;
903 // check if the node belongs to some station, save the MCId and outer radius
904 for ( Int_t ikst = 0 ; ikst < fNStations ; ikst++ )
905 if ( geoNodeName.Contains(stationNames[ikst]) ) {
906 fStationNrFromMcId[stsNode->getMCid()] = ikst;
907 stationKnown = kTRUE;
908 }
909
910 if ( stationKnown ) continue;
911
912 // if not known, register it and save MCId
913 fStationNrFromMcId[stsNode->getMCid()] = fNStations;
914
915 // it will work only if the node name is organized as:
916 // for station name is "stsstationXX", where XX is the station number (f.e. XX=07 for station number 7)
917 // for sector name is "stsstationXXanythingHereToDistinguishDifferentSectors"
918 geoNodeName.Remove(12,geoNodeName.Length()-12);
919 stationNames[fNStations] = geoNodeName.Data();
920 fNStations++;
921
922 cout << "station #" << fNStations << " has MCID = " << stsNode->getMCid() << " and name " << stsNode->GetName() << endl;
923
924 // fStationsMCId[fNStations] = stsNode->getMCid(); // not used
925 }
926 cout << "There are " << fNStations << " stations" << endl;
927
928 return kSUCCESS;
929
930}
931// -------------------------------------------------------------------------
932
933
934
935// ----- Private method Finish -----------------------------------------
936void CbmStsMatchHits::Finish() {
937
938 cout << endl;
939 cout << "============================================================"
940 << endl;
941 cout << "===== " << GetName() << ": Run summary " << endl;
942 cout << "===== " << endl;
943 cout << "===== Good events : " << setw(6) << fNEvents << endl;
944 cout << "===== Failed events : " << setw(6) << fNEventsFailed << endl;
945 cout << "===== Average time : " << setprecision(4) << setw(8) << right
946 << fTime / Double_t(fNEvents) << " s" << endl;
947 cout << "===== " << endl;
948 cout << "===== Hits per event : " << fixed << setprecision(0)
949 << fNHits / Double_t(fNEvents) << endl;
950 cout << setprecision(2);
951 cout << "===== Matched hits : " << fixed << setw(6) << right
952 << fNMatched / fNHits * 100. << " %" << endl;
953 cout << "===== Distant hits : " << fixed << setw(6) << right
954 << fNDistant / fNHits * 100. << " %" << endl;
955 cout << "===== Background hits : " << fixed << setw(6) << right
956 << fNBackgrd / fNHits * 100. << " % " << endl;
957 cout << "============================================================"
958 << endl;
959
960}
961// -------------------------------------------------------------------------
dist
vector< Double_t > dist(vector< Double_t > qp, Double_t mu)
Definition BmnMath.cxx:869
fDigiPar
BmnSsdDigitizeParameters * fDigiPar
Digitisation parameters.
Definition BmnSsdDigitize.cxx:65
i
int i
Definition P4_F32vec4.h:22
CbmGeoStsPar::GetGeoSensitiveNodes
TObjArray * GetGeoSensitiveNodes()
Definition CbmGeoStsPar.h:20
CbmStsDigiMatch::GetRefIndex
Int_t GetRefIndex(Int_t i=0) const
Definition CbmStsDigiMatch.cxx:65
CbmStsDigiScheme::Print
void Print(Bool_t kLong=kFALSE)
Definition CbmStsDigiScheme.cxx:248
CbmStsDigiScheme::Instance
static CbmStsDigiScheme * Instance(int version=1)
Definition CbmStsDigiScheme.cxx:65
CbmStsDigiScheme::GetSensorByName
CbmStsSensor * GetSensorByName(TString sensorName)
Definition CbmStsDigiScheme.cxx:233
CbmStsDigiScheme::GetStationByNr
CbmStsStation * GetStationByNr(Int_t stationNr)
Definition CbmStsDigiScheme.cxx:193
CbmStsDigiScheme::GetSector
CbmStsSector * GetSector(Int_t stationNr, Int_t sectorNr)
Definition CbmStsDigiScheme.cxx:204
CbmStsHit::GetStationNr
virtual Int_t GetStationNr() const
Definition CbmStsHit.h:66
CbmStsHit::GetSectorNr
Int_t GetSectorNr() const
Definition CbmStsHit.h:68
CbmStsHit::GetDigi
Int_t GetDigi(Int_t iSide) const
Definition CbmStsHit.cxx:116
CbmStsMatchHits::Exec
virtual void Exec(Option_t *opt)
Definition CbmStsMatchHits.cxx:161
CbmStsMatchHits::ExecReal
virtual void ExecReal(Option_t *opt)
Definition CbmStsMatchHits.cxx:405
CbmStsPoint::GetY
Double_t GetY(Double_t z) const
Definition CbmStsPoint.cxx:112
CbmStsPoint::GetX
Double_t GetX(Double_t z) const
Definition CbmStsPoint.cxx:101
CbmStsPoint::GetXIn
Double_t GetXIn() const
Definition CbmStsPoint.h:69
CbmStsPoint::GetZIn
Double_t GetZIn() const
Definition CbmStsPoint.h:71
CbmStsPoint::GetYIn
Double_t GetYIn() const
Definition CbmStsPoint.h:70
CbmStsSector::GetType
Int_t GetType() const
Definition CbmStsSector.h:93
CbmStsSensor::GetStationNr
Int_t GetStationNr() const
Definition CbmStsSensor.h:85
CbmStsStation::GetZ
Double_t GetZ(Int_t it=0)
Definition CbmStsStation.cxx:152
FairGeoPassivePar::GetGeoPassiveNodes
TObjArray * GetGeoPassiveNodes()
Definition FairGeoPassivePar.h:32
run
-clang-format