BmnRoot 26.10.0-18-g4e3414f9
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
CbmStsFindHits.cxx
Go to the documentation of this file.
1//* $Id: */
2
3// -------------------------------------------------------------------------
4// ----- CbmStsFindHits source file -----
5// ----- Created 26/06/2008 by R. Karabowicz -----
6// -------------------------------------------------------------------------
7
8#include "CbmStsFindHits.h"
9
10#include "CbmGeoStsPar.h"
11#include "CbmStsCluster.h"
12#include "CbmStsDigiPar.h"
13#include "CbmStsDigiScheme.h"
14#include "CbmStsHit.h"
15#include "CbmStsSector.h"
16#include "CbmStsStation.h"
17#include "FairEventHeader.h"
18#include "FairField.h" //AZ
19#include "FairRootManager.h"
20#include "FairRunAna.h"
21#include "FairRuntimeDb.h"
22#include "TClonesArray.h"
23#include "TMath.h"
24
25#include <TRandom.h> //AZ
26#include <iomanip>
27
28using std::cerr;
29using std::cout;
30using std::endl;
31using std::fixed;
32using std::flush;
33using std::left;
34using std::map;
35using std::right;
36using std::set;
37using std::setprecision;
38using std::setw;
39
40// ----- Default constructor ------------------------------------------
41CbmStsFindHits::CbmStsFindHits()
42 : FairTask("STS Hit Finder", 1)
43 , fGeoPar(NULL)
44 , fDigiPar(NULL)
45 , fDigiScheme(CbmStsDigiScheme::Instance())
46 , fClusters(NULL)
47 , fHits(NULL)
48 , fClusterMapF()
49 , fClusterMapB()
50 , fTimer()
51 , fNHits(0)
52{}
53
54// ----- Standard constructor ------------------------------------------
55CbmStsFindHits::CbmStsFindHits(Int_t iVerbose)
56 : FairTask("STSRealFindHits", iVerbose)
57 , fGeoPar(NULL)
58 , fDigiPar(NULL)
59 , fDigiScheme(CbmStsDigiScheme::Instance())
60 , fClusters(NULL)
61 , fHits(NULL)
62 , fClusterMapF()
63 , fClusterMapB()
64 , fTimer()
65 , fNHits(0)
66{}
67
68// ----- Constructor with name -----------------------------------------
69CbmStsFindHits::CbmStsFindHits(const char* name, Int_t iVerbose)
70 : FairTask(name, iVerbose)
71 , fGeoPar(NULL)
72 , fDigiPar(NULL)
73 , fDigiScheme(CbmStsDigiScheme::Instance())
74 , fClusters(NULL)
75 , fHits(NULL)
76 , fClusterMapF()
77 , fClusterMapB()
78 , fTimer()
79 , fNHits(0)
80{}
81
82// ----- Destructor ----------------------------------------------------
83CbmStsFindHits::~CbmStsFindHits()
84{
85 if (fHits) {
86 fHits->Delete();
87 delete fHits;
88 }
89}
90// -------------------------------------------------------------------------
91
92// ----- Public method Exec --------------------------------------------
93void CbmStsFindHits::Exec(Option_t* opt)
94{
95
96 fTimer.Start();
97 Bool_t warn = kFALSE;
98
99 // Check for digi scheme
100 if (!fDigiScheme) {
101 cerr << "-E- " << fName << "::Exec: No digi scheme!" << endl;
102 return;
103 }
104
105 // Clear output array
106 // fHits->Clear();
107 fHits->Delete();
108
109 // Sort STS digis with respect to sectors
110 SortClusters();
111
112 // Find hits in sectors
113 Int_t nClustersF = 0;
114 Int_t nClustersB = 0;
115 Int_t nHits = 0;
116 Int_t nStations = fDigiScheme->GetNStations();
117 CbmStsStation* station = NULL;
118 for (Int_t iStation = 0; iStation < nStations; iStation++) {
119 station = fDigiScheme->GetStation(iStation);
120 // if (station->GetStationNr() > 4) continue; //AZ-220322 - only 4 stations
121
122 Int_t nClustersFInStation = 0;
123 Int_t nClustersBInStation = 0;
124 Int_t nHitsInStation = 0;
125 Int_t nSectors = station->GetNSectors();
126 for (Int_t iSector = 0; iSector < nSectors; iSector++) {
127 CbmStsSector* sector = station->GetSector(iSector);
128 set<Int_t> fSet, bSet;
129 if (fClusterMapF.find(sector) == fClusterMapF.end()) {
130 cout << "-E- " << fName << "::Exec: sector " << sector->GetSectorNr() << " of station "
131 << station->GetStationNr() << "not found in front map!" << endl;
132 warn = kTRUE;
133 continue;
134 }
135 fSet = fClusterMapF[sector];
136 if (sector->GetType() == 2 || sector->GetType() == 3) {
137 if (fClusterMapB.find(sector) == fClusterMapB.end()) {
138 cout << "-E- " << fName << "::Exec: sector " << sector->GetSectorNr() << " of station "
139 << station->GetStationNr() << "not found in back map!" << endl;
140 warn = kTRUE;
141 continue;
142 }
143 }
144 bSet = fClusterMapB[sector];
145 Int_t nClustersFInSector = fSet.size();
146 Int_t nClustersBInSector = bSet.size();
147 Int_t nHitsInSector = FindHits(station, sector, fSet, bSet);
148 if (fVerbose > 2)
149 cout << "Sector " << sector->GetSectorNr() << ", Clusters front " << nClustersFInSector
150 << ", Clusters Back " << nClustersBInSector << ", Hits " << nHitsInSector << endl;
151 nHitsInStation += nHitsInSector;
152 nClustersFInStation += nClustersFInSector;
153 nClustersBInStation += nClustersBInSector;
154 } // Sector loop
155
156 if (fVerbose > 1)
157 cout << "Total for station " << station->GetStationNr() << ": Clusters front " << nClustersFInStation
158 << ", clusters back " << nClustersBInStation << ", hits " << nHitsInStation << endl;
159 nClustersB += nClustersBInStation;
160 nClustersF += nClustersFInStation;
161 nHits += nHitsInStation;
162
163 } // Station loop
164
165 fTimer.Stop();
166 if (fVerbose) {
167 cout << endl;
168 // AZ cout << "-I- " << fName << ":Event summary" << endl;
169 cout << "-I- " << fName << ":Event " << FairRunAna::Instance()->GetEventHeader()->GetMCEntryNumber()
170 << " summary" << endl;
171 cout << " Clusters front side : " << nClustersF << endl;
172 cout << " Clusters back side : " << nClustersB << endl;
173 cout << " Hits created : " << nHits << endl;
174 cout << " Real time : " << fTimer.RealTime() << endl;
175 } else {
176 if (warn)
177 cout << "- ";
178 else
179 cout << "+ ";
180 cout << setw(15) << left << fName << ": " << setprecision(4) << setw(8) << fixed << right << fTimer.RealTime()
181 << " s, clusters " << nClustersF << " / " << nClustersB << ", hits: " << nHits << endl;
182 }
183
184 fNHits += nHits;
185}
186// -------------------------------------------------------------------------
187
188// ----- Private method SetParContainers -------------------------------
189void CbmStsFindHits::SetParContainers()
190{
191
192 // Get run and runtime database
193 FairRunAna* run = FairRunAna::Instance();
194 if (!run)
195 Fatal("SetParContainers", "No analysis run");
196
197 FairRuntimeDb* db = run->GetRuntimeDb();
198 if (!db)
199 Fatal("SetParContainers", "No runtime database");
200
201 // Get STS geometry parameter container
202 fGeoPar = (CbmGeoStsPar*)db->getContainer("CbmGeoStsPar");
203
204 // Get STS digitisation parameter container
205 fDigiPar = (CbmStsDigiPar*)db->getContainer("CbmStsDigiPar");
206}
207// -------------------------------------------------------------------------
208
209// ----- Private method Init -------------------------------------------
210InitStatus CbmStsFindHits::Init()
211{
212
213 // Get input array
214 FairRootManager* ioman = FairRootManager::Instance();
215 if (!ioman)
216 Fatal("Init", "No FairRootManager");
217 fClusters = (TClonesArray*)ioman->GetObject("StsCluster");
218
219 // Register output array
220 fHits = new TClonesArray("CbmStsHit", 1000);
221 ioman->Register("StsHit", "Hit in STS", fHits, kTRUE);
222
223 // Build digitisation scheme
224 Bool_t success = fDigiScheme->Init(fGeoPar, fDigiPar);
225 if (!success)
226 return kERROR;
227
228 // Create sectorwise cluster sets
229 MakeSets();
230
231 // Control output
232
233 if (fVerbose == 1 || fVerbose == 2)
234 fDigiScheme->Print(kFALSE);
235 else if (fVerbose > 2)
236 fDigiScheme->Print(kTRUE);
237 cout << "-I- " << fName << "::Init: "
238 << "STS digitisation scheme succesfully initialised" << endl;
239 cout << " Stations: " << fDigiScheme->GetNStations() << ", Sectors: " << fDigiScheme->GetNSectors()
240 << ", Channels: " << fDigiScheme->GetNChannels() << endl;
241
242 return kSUCCESS;
243}
244// -------------------------------------------------------------------------
245
246// ----- Private method ReInit -----------------------------------------
247InitStatus CbmStsFindHits::ReInit()
248{
249
250 // Clear digitisation scheme
251 fDigiScheme->Clear();
252
253 // Build new digitisation scheme
254 Bool_t success = fDigiScheme->Init(fGeoPar, fDigiPar);
255 if (!success)
256 return kERROR;
257
258 // Create sectorwise digi sets
259 MakeSets();
260
261 // Control output
262 if (fVerbose == 1 || fVerbose == 2)
263 fDigiScheme->Print(kFALSE);
264 else if (fVerbose > 2)
265 fDigiScheme->Print(kTRUE);
266 cout << "-I- " << fName << "::Init: "
267 << "STS digitisation scheme succesfully reinitialised" << endl;
268 cout << " Stations: " << fDigiScheme->GetNStations() << ", Sectors: " << fDigiScheme->GetNSectors()
269 << ", Channels: " << fDigiScheme->GetNChannels() << endl;
270
271 return kSUCCESS;
272}
273// -------------------------------------------------------------------------
274
275// ----- Private method MakeSets ---------------------------------------
276void CbmStsFindHits::MakeSets()
277{
278
279 fClusterMapF.clear();
280 fClusterMapB.clear();
281 Int_t nStations = fDigiScheme->GetNStations();
282 for (Int_t iStation = 0; iStation < nStations; iStation++) {
283 CbmStsStation* station = fDigiScheme->GetStation(iStation);
284 Int_t nSectors = station->GetNSectors();
285 for (Int_t iSector = 0; iSector < nSectors; iSector++) {
286 CbmStsSector* sector = station->GetSector(iSector);
287 set<Int_t> a;
288 fClusterMapF[sector] = a;
289 if (sector->GetType() == 2 || sector->GetType() == 3) {
290 set<Int_t> b;
291 fClusterMapB[sector] = b;
292 }
293 }
294 }
295}
296// -------------------------------------------------------------------------
297
298// ----- Private method SortClusters --------------------------------------
299void CbmStsFindHits::SortClusters()
300{
301
302 // Check input array
303 if (!fClusters) {
304 cout << "-E- " << fName << "::SortClusters: No input array!" << endl;
305 return;
306 }
307
308 // Clear sector cluster sets
309 map<CbmStsSector*, set<Int_t>>::iterator mapIt;
310 for (mapIt = fClusterMapF.begin(); mapIt != fClusterMapF.end(); mapIt++)
311 ((*mapIt).second).clear();
312 for (mapIt = fClusterMapB.begin(); mapIt != fClusterMapB.end(); mapIt++)
313 ((*mapIt).second).clear();
314
315 // Fill clusters into sets
316 CbmStsCluster* cluster = NULL;
317 CbmStsSector* sector = NULL;
318 Int_t stationNr = -1;
319 Int_t sectorNr = -1;
320 Int_t iSide = -1;
321 Int_t nClusters = fClusters->GetEntriesFast();
322 for (Int_t iClus = 0; iClus < nClusters; iClus++) {
323 cluster = (CbmStsCluster*)fClusters->At(iClus);
324 stationNr = cluster->GetStationNr();
325 sectorNr = cluster->GetSectorNr();
326 iSide = cluster->GetSide();
327 if (fDigiScheme->GetStationByNr(stationNr) == NULL)
328 continue; // AZ - reduced setup
329 sector = fDigiScheme->GetSector(stationNr, sectorNr);
330 if (iSide == 0) {
331 if (fClusterMapF.find(sector) == fClusterMapF.end()) {
332 cerr << "-E- " << fName << "::SortClusters:: sector " << sectorNr << " of station " << stationNr
333 << " not found in digi scheme (F)!" << endl;
334 continue;
335 }
336 fClusterMapF[sector].insert(iClus);
337 } else if (iSide == 1) {
338 if (fClusterMapB.find(sector) == fClusterMapB.end()) {
339 cerr << "-E- " << fName << "::SortClusters:: sector " << sectorNr << " of station " << stationNr
340 << " not found in digi scheme (B)!" << endl;
341 continue;
342 }
343 fClusterMapB[sector].insert(iClus);
344 }
345 }
346}
347// -------------------------------------------------------------------------
348
349// ----- Private method FindHits ---------------------------------------
350Int_t CbmStsFindHits::FindHits(CbmStsStation* station, CbmStsSector* sector, set<Int_t>& fSet, set<Int_t>& bSet)
351{
352
353 // AZ - get magnetic field (for Lorentz shifts in GEMs)
354 static Int_t ifield = -1;
355 if (ifield < 0) {
356 FairRunAna* run = FairRunAna::Instance();
357 FairField* magField = run->GetField();
358 ifield = TMath::Nint(TMath::Abs(magField->GetBy(0, 0, 135.0)));
359 // ifield = 8; //AZ-281021 Dirty trick to use fast digitizer
360 }
361
362 // Counter
363 Int_t nNew = 0;
364
365 // Get sector parameters
366 // Int_t detId = sector->GetDetectorId();
367 Int_t iType = sector->GetType();
368
369 Double_t rot = sector->GetRotation();
370 Double_t dx = sector->GetDx();
371 // Double_t dy = sector->GetDy();
372 Double_t stereoB = sector->GetStereoB();
373 Double_t stereoF = sector->GetStereoF();
374
375 // Double_t z = station->GetZ();
376 Int_t stationNr = station->GetStationNr();
377 Int_t sectorNr = sector->GetSectorNr();
378
379 // Some auxiliary values
380 Double_t sinrot = TMath::Sin(rot);
381 Double_t cosrot = TMath::Cos(rot);
382 // Double_t tanstr = TMath::Tan(stereoB);
383
384 // Calculate error matrix in sector system
385 // Double_t vX, vY, vXY;
386 // if ( iType == 1 ) {
387 // vX = dx / TMath::Sqrt(12.);
388 // vY = dy / TMath::Sqrt(12.);
389 // vXY = 0.;
390 // }
391 // else if ( iType == 2 || iType == 3 ) {
392 //
393 // if (stereoF==0.&&stereoB*180/TMath::Pi()<80) {
394 // vX = dx / TMath::Sqrt(12.);
395 // vY = dx / TMath::Sqrt(6.) / TMath::Abs(TMath::Tan(stereoB));
396 // vXY = -1. * dx * dx / 12. / TMath::Tan(stereoB);
397 // }
398 // else if (stereoF==0.&&stereoB*180/TMath::Pi()>80) {
399 // vX = dx / TMath::Sqrt(12.);
400 // vY = dx / TMath::Sqrt(12.);
401 // vXY = 0.;
402 // }
403 // else {
404 // vX = dx / TMath::Sqrt(24.);
405 // vY = dx / TMath::Sqrt(24.) / TMath::Abs(TMath::Tan(stereoB));
406 // vXY = 0.;
407 // }
408 // }
409 // else {
410 // cerr << "-E- " << fName << "::FindHits: Illegal sector type "
411 // << iType << endl;
412 // return 0;
413 // }
414
415 // Transform variances into global c.s.
416 // Double_t wX = vX * vX * cosrot * cosrot
417 // - 2. * vXY * cosrot * sinrot
418 // + vY * vY * sinrot * sinrot;
419 // Double_t wY = vX * vX * sinrot * sinrot
420 // + 2. * vXY * cosrot * sinrot
421 // + vY * vY * cosrot * cosrot;
422 // Double_t wXY = (vX*vX - vY*vY) * cosrot * sinrot
423 // + vXY * ( cosrot*cosrot - sinrot*sinrot );
424 // Double_t sigmaX = TMath::Sqrt(wX);
425 // Double_t sigmaY = TMath::Sqrt(wY);
426
427 // Now perform the loop over active channels
428 set<Int_t>::iterator it1;
429 set<Int_t>::iterator it2;
430
431 // ----- Type 1 : Pixel sector ---------------------------------------
432 if (iType == 1) {
433 Fatal("FindHits", "Sorry, not implemented yet");
434 } // Pixel sensor
435 // ---------------------------------------------------------------------
436
437 // ----- Type 2: Strip sector OSU -----------------------------------
438 /*AZ
439 else if ( iType == 2 ) {
440 Fatal("FindHits","Sorry, not implemented yet");
441 } // Strip OSU
442 */
443 // ---------------------------------------------------------------------
444
445 // ----- Type 3: Strip sector GSI -----------------------------------
446 // AZ else if (iType == 3 ) {
447 else if (iType == 2 || iType == 3)
448 {
449 Int_t iClusF = -1;
450 Int_t iClusB = -1;
451 Double_t chanF = -1;
452 Double_t chanB = -1;
453 Int_t nHits = fHits->GetEntriesFast();
454 Double_t xHit;
455 Double_t yHit;
456 Double_t zHit;
457 TVector3 pos, dpos;
458 CbmStsCluster* clusterF = NULL;
459 CbmStsCluster* clusterB = NULL;
460
461 Double_t vX, vY, vXY;
462 if (stereoF == 0. && stereoB * 180 / TMath::Pi() < 80) {
463 vX = dx / TMath::Sqrt(12.);
464 vY = dx / TMath::Sqrt(6.) / TMath::Abs(TMath::Tan(stereoB));
465 vXY = -1. * dx * dx / 12. / TMath::Tan(stereoB);
466 } else if (stereoF == 0. && stereoB * 180 / TMath::Pi() > 80) {
467 vX = dx / TMath::Sqrt(12.);
468 vY = dx / TMath::Sqrt(12.);
469 vXY = 0.;
470 } else {
471 vX = dx / TMath::Sqrt(12.);
472 vY = dx / TMath::Sqrt(12.) / TMath::Abs(TMath::Tan(stereoB));
473 vXY = 0.;
474 }
475
476 for (it1 = fSet.begin(); it1 != fSet.end(); it1++) {
477 iClusF = (*it1);
478 clusterF = (CbmStsCluster*)fClusters->At(iClusF);
479 if (!clusterF) {
480 cout << "-W- " << GetName() << "::FindHits: Invalid cluster index " << iClusF
481 << " in front set of sector " << sector->GetSectorNr() << ", station " << station->GetStationNr()
482 << endl;
483 continue;
484 }
485 chanF = clusterF->GetMean();
486 for (it2 = bSet.begin(); it2 != bSet.end(); it2++) {
487 iClusB = (*it2);
488 clusterB = (CbmStsCluster*)fClusters->At(iClusB);
489 if (!clusterB) {
490 cout << "-W- " << GetName() << "::FindHits: Invalid cluster index " << iClusB
491 << " in back set of sector " << sector->GetSectorNr() << ", station "
492 << station->GetStationNr() << endl;
493 continue;
494 }
495 chanB = clusterB->GetMean();
496
497 // AZ Add extra smearing
498 /* Wrong - no correlation ???
499 Double_t sigma = 0.3; // 800*0.3=240 um
500 chanF += gRandom->Gaus(0.0,sigma);
501 chanB += gRandom->Gaus(0.0,sigma);
502 */
503 // AZ
504
505 Int_t sensorDetId = sector->IntersectClusters(chanF, chanB, xHit, yHit, zHit);
506
507 if (sensorDetId == -1)
508 continue;
509
510 // AZ Add extra smearing
511 /*
512 Double_t sigma = 0.06; // 400 um
513 Double_t dx = gRandom->Gaus(0.0,sigma);
514 xHit += dx;
515 yHit += gRandom->Gaus(0.0,sigma*3.5);
516 */
517 // AZ
518 // AZ - Lorentz shift (Ar:CO2 70:30)
519 ///*
520 // Double_t lorShift = 0.1, zShift = -0.3; // (Ar:CO2 70:30)
521 Double_t lorShift = 0.1575, zShift = -0.27; // (Ar:C4H10 80:20 - 0.8T)
522 if (ifield == 4) {
523 lorShift = 0.0765;
524 zShift = -0.27;
525 } // (Ar:C4H10 80:20 - 0.4T - Au 1.5 GeV)
526 else if (ifield == 6)
527 {
528 lorShift = 0.1145;
529 zShift = -0.27;
530 } // (Ar:C4H10 80:20 - 0.6T - Au 2.5 GeV)
531 // else if (ifield == 6) { lorShift = 0.1190; zShift = -0.27; } // (Ar:C4H10 80:20 - 0.6T - Au 2.5 GeV)
532 else if (ifield == 7)
533 {
534 lorShift = 0.1360;
535 zShift = -0.27;
536 } // (Ar:C4H10 80:20 - 0.7T - Au 3.5 GeV)
537 else if (ifield == 9)
538 {
539 lorShift = 0.1800;
540 zShift = -0.27;
541 } // (Ar:C4H10 80:20 - 0.9T - Au 4.5 GeV)
542
543 // if (station->GetStationNr() == 1) lorShift = zShift = 0.0; //Si - run6
544 // if (station->GetStationNr() <= 3) lorShift = zShift = 0.0; //Si - run7
545 if (sector->GetDx() < 0.02)
546 lorShift = zShift = 0.0; // Si
547 else if (station->GetStationNr() % 2 == 0) {
548 lorShift *= -1;
549 zShift *= -1;
550 }
551 // zShift = (station->GetStationNr() % 2 == 0) ? -0.27 : 0.27; //AZ-220322
552 xHit += lorShift;
553 // xHit += gRandom->Gaus(0.0,0.15); //AZ-220322
554 // AZ zHit += zShift; //???
555 //*/
556 // AZ
557
558 pos.SetXYZ(xHit, yHit, zHit);
559
560 Double_t vXTemp, vYTemp, vXYTemp;
561
562 Double_t wX = vX * vX * cosrot * cosrot - 2. * vXY * cosrot * sinrot + vY * vY * sinrot * sinrot;
563 Double_t wY = vX * vX * sinrot * sinrot + 2. * vXY * cosrot * sinrot + vY * vY * cosrot * cosrot;
564 Double_t wXY = (vX * vX - vY * vY) * cosrot * sinrot + vXY * (cosrot * cosrot - sinrot * sinrot);
565 // Double_t sigmaX = TMath::Sqrt(wX);
566 // Double_t sigmaY = TMath::Sqrt(wY);
567
568 if (stereoF == 0. && stereoB * 180 / TMath::Pi() < 80) { // 0&15 case
569
570 vXTemp = vX * clusterF->GetMeanError();
571 vYTemp = (vX / (TMath::Tan(stereoB)))
572 * TMath::Sqrt(clusterF->GetMeanError() * clusterF->GetMeanError()
573 + clusterB->GetMeanError() * clusterB->GetMeanError());
574 vXYTemp = -vXTemp * vXTemp / (TMath::Tan(stereoB));
575
576 wX = vXTemp * vXTemp * cosrot * cosrot - 2. * vXYTemp * cosrot * sinrot
577 + vYTemp * vYTemp * sinrot * sinrot;
578 wY = vXTemp * vXTemp * sinrot * sinrot + 2. * vXYTemp * cosrot * sinrot
579 + vYTemp * vYTemp * cosrot * cosrot;
580 wXY = (vXTemp * vXTemp - vYTemp * vYTemp) * cosrot * sinrot
581 + vXYTemp * (cosrot * cosrot - sinrot * sinrot);
582 // dpos.SetXYZ(TMath::Sqrt(wX),TMath::Sqrt(wY), 0.);
583 // dpos.SetXYZ(vX, vX, 0.); //AZ - errors on each side
584 // if (zHit > 10.0) dpos.SetXYZ(vX*6, vX*6, 0.); //AZ - errors on each side
585 // if (zHit > 10.0) dpos.SetXYZ(vX*2, vX*2, 0.); //AZ - errors on each side
586 // else dpos.SetXYZ(vX*2, vX*10, 0.); //AZ - errors on each side
587 /*run6
588 if (zHit > 10.0) dpos.SetXYZ(0.015, 0.02, 0.); //AZ - errors on each side
589 else dpos.SetXYZ(0.002, 0.04, 0.); //AZ - errors on each side
590 */
591 // if (zHit > 35.0) dpos.SetXYZ(0.015, 0.058, 0.); //AZ - run7 - /sin(15)
592 // else dpos.SetXYZ(0.01, 0.24, 0.); //AZ - errors on each side
593 // if (zHit > 85.0) dpos.SetXYZ(0.015, 0.058, 0.); //AZ - CBM+BM@N
594 // else dpos.SetXYZ(0.01, 0.077, 0.); //AZ - errors on each side
595 if (sector->GetDx() < 0.02)
596 dpos.SetXYZ(sector->GetDx() / TMath::Sqrt(12.),
597 sector->GetDx() / TMath::Sqrt(12.) / TMath::Sin(7.5 * TMath::DegToRad()), 0.);
598 else
599 dpos.SetXYZ(0.015, 0.058, 0.); // GEMS
600 // else dpos.SetXYZ(vX*10, vX*100, 0.); //AZ - errors on each side
601 vXTemp = clusterF->GetMeanError(); // AZ
602 vYTemp = clusterB->GetMeanError(); // AZ
603 // dpos.SetXYZ(vXTemp, vYTemp, 0.); //AZ - errors on each side
604 } else if (stereoF == 0. && stereoB * 180 / TMath::Pi() > 80) { // 0&90 case
605
606 // vXTemp = vX * clusterF->GetMeanError();
607 // vYTemp = vX * clusterB->GetMeanError();
608 // vXTemp = vX * 1.; // * clusterF->GetMeanError(); //AZ
609 // vYTemp = vX * 1.; // * clusterB->GetMeanError(); //AZ
610 vXTemp = clusterF->GetMeanError(); // AZ
611 vYTemp = clusterB->GetMeanError(); // AZ
612 vXYTemp = 0.;
613
614 wX = vXTemp * vXTemp * cosrot * cosrot - 2. * vXYTemp * cosrot * sinrot
615 + vYTemp * vYTemp * sinrot * sinrot;
616 wY = vXTemp * vXTemp * sinrot * sinrot + 2. * vXYTemp * cosrot * sinrot
617 + vYTemp * vYTemp * cosrot * cosrot;
618 wXY = (vXTemp * vXTemp - vYTemp * vYTemp) * cosrot * sinrot
619 + vXYTemp * (cosrot * cosrot - sinrot * sinrot);
620 // dpos.SetXYZ(TMath::Sqrt(wX),TMath::Sqrt(wY), 0.);
621 dpos.SetXYZ(vX * 6, vX * 6, 0.); // AZ - errors on each side
622 // dpos.SetXYZ(vXTemp, vYTemp, 0.); //AZ - errors on each side
623 }
624
625 else
626 { // 7.5&-7.5 case
627
628 vXTemp = (vX / 2.)
629 * TMath::Sqrt(clusterF->GetMeanError() * clusterF->GetMeanError()
630 + clusterB->GetMeanError() * clusterB->GetMeanError());
631 vYTemp = (vX / (2. * (TMath::Tan(stereoB))))
632 * TMath::Sqrt(clusterF->GetMeanError() * clusterF->GetMeanError()
633 + clusterB->GetMeanError() * clusterB->GetMeanError());
634 vXYTemp = (vX * vX / (4. * (TMath::Tan(stereoB))))
635 * (clusterB->GetMeanError() * clusterB->GetMeanError()
636 - clusterF->GetMeanError() * clusterF->GetMeanError());
637
638 wX = vXTemp * vXTemp * cosrot * cosrot - 2. * vXYTemp * cosrot * sinrot
639 + vYTemp * vYTemp * sinrot * sinrot;
640 wY = vXTemp * vXTemp * sinrot * sinrot + 2. * vXYTemp * cosrot * sinrot
641 + vYTemp * vYTemp * cosrot * cosrot;
642 wXY = (vXTemp * vXTemp - vYTemp * vYTemp) * cosrot * sinrot
643 + vXYTemp * (cosrot * cosrot - sinrot * sinrot);
644 dpos.SetXYZ(TMath::Sqrt(wX), TMath::Sqrt(wY), 0.);
645 }
646
647 Int_t statLayer = -1;
648 for (Int_t istatL = station->GetNofZ(); istatL > 0; istatL--)
649 if (TMath::Abs(zHit - station->GetZ(istatL - 1)) < 0.00001) {
650 statLayer = istatL - 1;
651 break;
652 }
653
654 if (statLayer == -1)
655 cout << "unknown layer for hit at z = " << zHit << endl;
656
657 if (sector->GetDx() > 0.03)
658 pos[2] += zShift; // AZ - effective position (middle of GEM drift region)
659 // AZ new ((*fHits)[nHits++]) CbmStsHit(sensorDetId, pos, dpos, wXY,
660 CbmStsHit* hhNew = new ((*fHits)[nHits++])
661 CbmStsHit(sensorDetId, pos, dpos, wXY, iClusF, iClusB, (Int_t)chanF, (Int_t)chanB, statLayer);
662 // max(iClusF-10,0), iClusB+1, max((Int_t)chanF-1,0), (Int_t)chanB+1, statLayer); //AZ-230322 - test
663 hhNew->SetSignalDiv(TMath::Min(clusterF->GetQtot(), clusterB->GetQtot())
664 / TMath::Max(clusterF->GetQtot(), clusterB->GetQtot())); // AZ
665 nNew++;
666 if (fVerbose > 3)
667 cout << "New StsHit at (" << xHit << ", " << yHit << ", " << zHit << "), station " << stationNr
668 << ", sector " << sectorNr << ", channel " << chanF << " / " << chanB << endl;
669
670 } // back side strip loop
671 } // front side strip loop
672
673 } // strip GSI
674 // ---------------------------------------------------------------------
675
676 return nNew;
677}
678// -------------------------------------------------------------------------
679
680// ----- Virtual method Finish -----------------------------------------
681void CbmStsFindHits::Finish()
682{
683 cout << endl;
684 cout << "============================================================" << endl;
685 cout << "===== " << fName << ": Run summary " << endl;
686 cout << "===== " << endl;
687 cout << "===== Number of hits : " << setw(8) << setprecision(2) << fNHits << endl;
688 cout << "============================================================" << endl;
689}
690// -------------------------------------------------------------------------
Bool_t
bool Bool_t
Int_t
int Int_t
kFALSE
constexpr Bool_t kFALSE
Double_t
double Double_t
kTRUE
constexpr Bool_t kTRUE
Option_t
const char Option_t
Fatal
void Fatal(const char *location, const char *msgfmt,...)
CbmStsCluster
CbmStsCluster.h.
Definition CbmStsCluster.h:32
CbmStsCluster::GetStationNr
Int_t GetStationNr() const
Definition CbmStsCluster.h:59
CbmStsCluster::GetSectorNr
Int_t GetSectorNr() const
Definition CbmStsCluster.h:60
CbmStsCluster::GetMean
Double_t GetMean() const
Definition CbmStsCluster.h:68
CbmStsCluster::GetMeanError
Double_t GetMeanError() const
Definition CbmStsCluster.h:69
CbmStsCluster::GetQtot
Double_t GetQtot() const
Definition CbmStsCluster.h:70
CbmStsCluster::GetSide
Int_t GetSide() const
Definition CbmStsCluster.h:64
CbmStsFindHits::MakeSets
void MakeSets()
Make sectorwise sets for sigis.
Definition CbmStsFindHits.cxx:276
CbmStsFindHits::SortClusters
void SortClusters()
Sort digis sectorwise.
Definition CbmStsFindHits.cxx:299
CbmStsFindHits::FindHits
Int_t FindHits(CbmStsStation *station, CbmStsSector *sector, std::set< Int_t > &fSet, std::set< Int_t > &bSet)
Find hits in one sector.
CbmStsFindHits::~CbmStsFindHits
virtual ~CbmStsFindHits()
Destructor.
Definition CbmStsFindHits.cxx:83
CbmStsFindHits::fClusterMapB
std::map< CbmStsSector *, std::set< Int_t > > fClusterMapB
sector clusters (front)
Definition CbmStsFindHits.h:75
CbmStsFindHits::CbmStsFindHits
CbmStsFindHits()
Default constructor.
Definition CbmStsFindHits.cxx:41
CbmStsFindHits::ReInit
virtual InitStatus ReInit()
Reinitialisation.
Definition CbmStsFindHits.cxx:247
CbmStsFindHits::fClusterMapF
std::map< CbmStsSector *, std::set< Int_t > > fClusterMapF
Output array of CbmStsHit.
Definition CbmStsFindHits.h:74
CbmStsFindHits::fGeoPar
CbmGeoStsPar * fGeoPar
Definition CbmStsFindHits.h:66
CbmStsFindHits::fHits
TClonesArray * fHits
SM - Input array of BMN Silicon hits.
Definition CbmStsFindHits.h:73
CbmStsFindHits::Finish
virtual void Finish()
Virtual method Finish.
Definition CbmStsFindHits.cxx:681
CbmStsFindHits::fClusters
TClonesArray * fClusters
Digitisation scheme.
Definition CbmStsFindHits.h:69
CbmStsFindHits::fDigiPar
CbmStsDigiPar * fDigiPar
Geometry parameters.
Definition CbmStsFindHits.h:67
CbmStsFindHits::fTimer
TStopwatch fTimer
sector clusters (back)
Definition CbmStsFindHits.h:76
CbmStsFindHits::Exec
virtual void Exec(Option_t *opt)
Execution.
Definition CbmStsFindHits.cxx:93
CbmStsFindHits::Init
virtual InitStatus Init()
Intialisation.
Definition CbmStsFindHits.cxx:210
CbmStsFindHits::fDigiScheme
CbmStsDigiScheme * fDigiScheme
Digitisation parameters.
Definition CbmStsFindHits.h:68
CbmStsFindHits::SetParContainers
virtual void SetParContainers()
Get parameter containers.
Definition CbmStsFindHits.cxx:189
CbmStsHit::SetSignalDiv
void SetSignalDiv(Double_t sigDiv)
Definition CbmStsHit.h:91
CbmStsSector::GetType
Int_t GetType() const
Definition CbmStsSector.h:81
CbmStsSector::GetSectorNr
Int_t GetSectorNr() const
Definition CbmStsSector.h:63
CbmStsSector::GetDx
Double_t GetDx() const
Definition CbmStsSector.h:83
CbmStsSector::IntersectClusters
Int_t IntersectClusters(Double_t fChan, Double_t bChan, Double_t &xCross, Double_t &yCross, Double_t &zCross)
Definition CbmStsSector.cxx:314
CbmStsSector::GetStereoF
Double_t GetStereoF() const
Definition CbmStsSector.h:85
CbmStsSector::GetRotation
Double_t GetRotation() const
Definition CbmStsSector.h:82
CbmStsSector::GetStereoB
Double_t GetStereoB() const
Definition CbmStsSector.h:86
CbmStsStation::GetNSectors
Int_t GetNSectors() const
Definition CbmStsStation.h:62
CbmStsStation::GetZ
Double_t GetZ(Int_t it=0)
Definition CbmStsStation.cxx:151
CbmStsStation::GetStationNr
Int_t GetStationNr() const
Definition CbmStsStation.h:56
CbmStsStation::GetSector
CbmStsSector * GetSector(Int_t iSector)
Definition CbmStsStation.h:65
TVector3::SetXYZ
void SetXYZ(Double_t x, Double_t y, Double_t z)
TMath::Min
Double_t Min(Double_t a, Double_t b)
TMath::Nint
Int_t Nint(T x)
TMath::DegToRad
constexpr Double_t DegToRad()
TMath::Sqrt
Double_t Sqrt(Double_t x)
TMath::Cos
Double_t Cos(Double_t)
TMath::Pi
constexpr Double_t Pi()
TMath::Abs
Double_t Abs(Double_t d)
TMath::Sin
Double_t Sin(Double_t)
TMath::Max
Double_t Max(Double_t a, Double_t b)
TMath::Tan
Double_t Tan(Double_t)