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