BmnRoot
Loading...
Searching...
No Matches
CbmStsRealDigitize.cxx
Go to the documentation of this file.
1//* $Id: */
2
3// -------------------------------------------------------------------------
4// ----- CbmStsRealDigitize source file -----
5// ----- Created 08/07/2008 by R. Karabowicz -----
6// -------------------------------------------------------------------------
7
8
9
10
11// Includes from ROOT
12#include "TClonesArray.h"
13#include "TObjArray.h"
14#include "TMath.h"
15#include "TF1.h"
16#include "TRandom3.h"
17
18#include "TGeoManager.h"
19#include "TGeoNode.h"
20
21// Includes from base
22#include "FairRootManager.h"
23#include "FairRunAna.h"
24#include "FairRuntimeDb.h"
25
26// Includes from STS
27#include "CbmGeoStsPar.h"
28#include "CbmStsDigi.h"
29#include "CbmStsDigiMatch.h"
30#include "CbmStsDigiPar.h"
31#include "CbmStsDigiScheme.h"
32#include "CbmStsRealDigitize.h"
33#include "CbmStsPoint.h"
34#include "CbmStsSensor.h"
35#include "CbmStsSector.h"
36#include "CbmStsStation.h"
37
38#include <iostream>
39#include <iomanip>
40#include <map>
41
42using std::cout;
43using std::cerr;
44using std::endl;
45using std::flush;
46using std::pair;
47using std::setw;
48using std::left;
49using std::right;
50using std::fixed;
51using std::setprecision;
52using std::set;
53using std::map;
54using std::ios_base;
55
56
57// ----- Default constructor ------------------------------------------
58CbmStsRealDigitize::CbmStsRealDigitize() : FairTask("STS Digitizer", 1) {
59 fGeoPar = NULL;
60 fDigiPar = NULL;
61 fPoints = NULL;
62 fDigis = NULL;
63 fDigiMatches = NULL;
64 fRealistic = kFALSE;
65 fDigiScheme = new CbmStsDigiScheme();
66 Reset();
67
68 fStep = 0.001;
69
70 fFThreshold = 0.001;
71 fBThreshold = 0.001;
72 fFNoiseWidth = 0.0;
73 fBNoiseWidth = 0.0;
74
75 fFNofBits = 20;
76 fBNofBits = 20;
77 fFMinStep = 0.01;
78 fBMinStep = 0.01;
79
80 fNEvents = 0.;
81}
82// -------------------------------------------------------------------------
83
84// ----- Standard constructor ------------------------------------------
86 : FairTask("STSDigitize", iVerbose) {
87 fGeoPar = NULL;
88 fDigiPar = NULL;
89 fPoints = NULL;
90 fDigis = NULL;
91 fDigiMatches = NULL;
92 fRealistic = kFALSE;
93 fDigiScheme = new CbmStsDigiScheme();
94 Reset();
95
96 fStep = 0.001;
97
98 fFThreshold = 0.001;
99 fBThreshold = 0.001;
100 fFNoiseWidth = 0.0;
101 fBNoiseWidth = 0.0;
102
103 fFNofBits = 20;
104 fBNofBits = 20;
105 fFMinStep = 0.01;
106 fBMinStep = 0.01;
107
108 fNEvents = 0.;
109}
110// -------------------------------------------------------------------------
111
112// ----- Constructor with name -----------------------------------------
113CbmStsRealDigitize::CbmStsRealDigitize(const char* name, Int_t iVerbose)
114 : FairTask(name, iVerbose) {
115 fGeoPar = NULL;
116 fDigiPar = NULL;
117 fPoints = NULL;
118 fDigis = NULL;
119 fDigiMatches = NULL;
120 fRealistic = kFALSE;
121 fDigiScheme = new CbmStsDigiScheme();
122 Reset();
123
124 fStep = 0.001;
125
126 fFThreshold = 0.001;
127 fBThreshold = 0.001;
128 fFNoiseWidth = 0.0;
129 fBNoiseWidth = 0.0;
130
131 fFNofBits = 20;
132 fBNofBits = 20;
133 fFMinStep = 0.01;
134 fBMinStep = 0.01;
135
136 fNEvents = 0.;
137}
138// -------------------------------------------------------------------------
139
140// ----- Destructor ----------------------------------------------------
142 if ( fGeoPar) delete fGeoPar;
143 if ( fDigiPar) delete fDigiPar;
144 if ( fDigis ) {
145 fDigis->Delete();
146 delete fDigis;
147 }
148 if ( fDigiMatches ) {
149 fDigiMatches->Delete();
150 delete fDigiMatches;
151 }
152 if ( fDigiScheme ) delete fDigiScheme;
153 Reset();
154}
155// -------------------------------------------------------------------------
156
157// ----- Public method Exec --------------------------------------------
158void CbmStsRealDigitize::Exec(Option_t* opt) {
159
160// cout << "threshold = " << fFThreshold << endl;
161
162 // Reset all eventwise counters
163 fTimer.Start();
164 Reset();
165 TF1* digiGausDist = new TF1("digiGausDist","gaus",-5.,5.);
166
167 // Verbose screen output
168 if ( fVerbose > 2 ) {
169 cout << endl << "-I- " << fName << ": executing event" << endl;
170 cout << "----------------------------------------------" << endl;
171 }
172
173 Double_t nPoints = 0.;
174 Double_t nOutside = 0.;
175 Double_t nDigisF = 0.;
176 Double_t nDigisB = 0.;
177
178 // Loop over all StsPoints
179 if ( ! fPoints ) {
180 cerr << "-W- " << fName << "::Exec: No input array (STSPoint) "
181 << endl;
182 cout << "- " << fName << endl;
183 return;
184 }
185
186 map<CbmStsSensor*, set<Int_t> >::iterator mapIt;
187 for (mapIt=fPointMap.begin(); mapIt!=fPointMap.end(); mapIt++)
188 ((*mapIt).second).clear();
189
190 for (Int_t iPoint=0; iPoint<fPoints->GetEntriesFast(); iPoint++) {
191 CbmStsPoint* point = (CbmStsPoint*) fPoints->At(iPoint);
192
193 Double_t xin = point->GetXIn();
194 Double_t yin = point->GetYIn();
195 Double_t zin = point->GetZIn();
196 gGeoManager->FindNode(xin,yin,zin);
197 TGeoNode* curNode = gGeoManager->GetCurrentNode();
198
199 CbmStsSensor* sensor = fDigiScheme->GetSensorByName(curNode->GetName());
200
201 if ( fPointMap.find(sensor) == fPointMap.end() ) {
202 cerr << "-E- " << fName << "::Exec:: sensor " << curNode->GetName()
203 << " not found in digi scheme!" << endl;
204 continue;
205 }
206 fPointMap[sensor].insert(iPoint);
207 }
208
209 for (Int_t iStation=fDigiScheme->GetNStations(); iStation > 0 ; ) {
210 CbmStsStation* station = fDigiScheme->GetStation(--iStation);
211 for (Int_t iSector=station->GetNSectors(); iSector > 0 ; ) {
212 CbmStsSector* sector = station->GetSector(--iSector);
213
214 map<Int_t, set<Int_t> >::iterator mapCh;
215
216 for (mapCh=fFChannelPointsMap.begin(); mapCh!=fFChannelPointsMap.end(); mapCh++)
217 ((*mapCh).second).clear();
218 for (mapCh=fBChannelPointsMap.begin(); mapCh!=fBChannelPointsMap.end(); mapCh++)
219 ((*mapCh).second).clear();
220
221 // simulating detector+cables+electronics noise
222 // should be more sophisticated...
223 // the question is: sectorwise or sensorwise???
224 Int_t nChannels = sector->GetNChannelsFront();
225 for (Int_t iChannel=nChannels ; iChannel > 0 ; ) {
226// fStripSignalF[--iChannel] = fGen->Landau(.1,.02);
227// fStripSignalB[ iChannel] = fGen->Landau(.1,.02);
228// fStripSignalF[--iChannel] = 0.;
229// fStripSignalB[ iChannel] = 0.;
230 fStripSignalF[--iChannel] = TMath::Abs(gRandom->Gaus(0.,fFNoiseWidth));
231 fStripSignalB[ iChannel] = TMath::Abs(gRandom->Gaus(0.,fBNoiseWidth));
232 }
233
234 for (Int_t iSensor=sector->GetNSensors(); iSensor > 0 ; ) {
235 CbmStsSensor* sensor = sector->GetSensor(--iSensor);
236
237 ProduceHitResponse(sensor);
238 }
239
240 Int_t stationNr = sector->GetStationNr();
241 Int_t sectorNr = sector->GetSectorNr();
242 Int_t sectorDetId = sector->GetDetectorId();
243
244 for ( Int_t ifstr = 0 ; ifstr < nChannels ; ifstr++ ) {
245 if ( fStripSignalF[ifstr] < fFThreshold ) continue;
246// cout << "digi#" << fNDigis << " -> making fdigi at " << stationNr << "," << sectorNr
247// << " at channel " << ifstr << " with signal " << fStripSignalF[ifstr] << endl;
248 Int_t digiFSignal = 1+(Int_t)((fStripSignalF[ifstr]-fFThreshold)/fFMinStep);
249 if ( digiFSignal >= fFNofSteps ) digiFSignal = fFNofSteps-1;
250 new (( *fDigis)[fNDigis]) CbmStsDigi(0, stationNr, sectorNr,
251 0, ifstr,
252 digiFSignal, 0);
253 // fStripSignalF[ifstr], 0);
254 set<Int_t>::iterator it1;
255 set<Int_t> chPnt = fFChannelPointsMap[ifstr];
256 Int_t pnt;
258 if ( chPnt.size() == 0 ) {
259// cout << "digi#" << fNDigis << " -> making fdigi at " << stationNr << "," << sectorNr
260// << " at channel " << ifstr << " with signal " << fStripSignalF[ifstr] << endl;
261 new ((*fDigiMatches)[fNDigis]) CbmStsDigiMatch(-666);
262 }
263 else {
264 for (it1=chPnt.begin(); it1!=chPnt.end(); it1++) {
265 pnt = (*it1);
266 if ( it1==chPnt.begin() )
267 match = new ((*fDigiMatches)[fNDigis]) CbmStsDigiMatch(pnt);
268 else {
269 match->AddPoint(pnt);
270 fNMulti++;
271 }
272 }
273 }
274 fNDigis++;
275 }
276
277 for ( Int_t ibstr = 0 ; ibstr < nChannels ; ibstr++ ) {
278 if ( fStripSignalB[ibstr] < fBThreshold ) continue;
279 Int_t digiBSignal = 1+(Int_t)((fStripSignalB[ibstr]-fBThreshold)/fBMinStep);
280 if ( digiBSignal >= fBNofSteps ) digiBSignal = fBNofSteps-1;
281 new (( *fDigis)[fNDigis]) CbmStsDigi(0, stationNr, sectorNr,
282 1, ibstr,
283 digiBSignal, 0);
284 // fStripSignalB[ibstr], 0);
285 set<Int_t>::iterator it1;
286 set<Int_t> chPnt = fBChannelPointsMap[ibstr];
287 Int_t pnt;
289 if ( chPnt.size() == 0 ) {
290// cout << "digi#" << fNDigis << " -> making fdigi at " << stationNr << "," << sectorNr
291// << " at channel " << ifstr << " with signal " << fStripSignalF[ifstr] << endl;
292 new ((*fDigiMatches)[fNDigis]) CbmStsDigiMatch(-666);
293 }
294 else {
295 for (it1=chPnt.begin(); it1!=chPnt.end(); it1++) {
296 pnt = (*it1);
297 if ( it1==chPnt.begin() )
298 match = new ((*fDigiMatches)[fNDigis]) CbmStsDigiMatch(pnt);
299 else {
300 match->AddPoint(pnt);
301 fNMulti++;
302 }
303 }
304 }
305 fNDigis++;
306 }
307 }
308 }
309
310 fTimer.Stop();
311 cout << "+ " << flush;
312 cout << setw(15) << left << fName << ": " << setprecision(4) << setw(8)
313 << fixed << right << fTimer.RealTime()
314 << " s, digis " << nDigisF << " / " << nDigisB << endl;
315
316 fNEvents += 1.;
317 fNPoints += Double_t(nPoints);
318 fNOutside += Double_t(nOutside);
319 fNDigisFront += Double_t(nDigisF);
320 fNDigisBack += Double_t(nDigisB);
321 fTime += fTimer.RealTime();
322
323}
324// -------------------------------------------------------------------------
325
326// ----- Private method ProduceHitResponse --------------------------------
328 set <Int_t> pSet;
329 if ( fPointMap.find(sensor) == fPointMap.end() ) {
330 cerr << "-E- " << fName << "::Exec:: sensor"
331 << " not found in digi scheme!" << endl;
332 return;
333 }
334 pSet = fPointMap[sensor];
335
336 Int_t iPoint = -1;
337 CbmStsPoint* point = NULL;
338
339 set<Int_t>::iterator it1;
340
341 for (it1=pSet.begin(); it1!=pSet.end(); it1++) {
342 iPoint = (*it1);
343 point = (CbmStsPoint*) fPoints->At(iPoint);
344
345 Double_t xin = point->GetXIn();
346 Double_t yin = point->GetYIn();
347 Double_t zin = point->GetZIn();
348
349 Double_t xvec = point->GetXOut()-xin;
350 Double_t yvec = point->GetYOut()-yin;
351 Double_t zvec = point->GetZOut()-zin;
352
353 Int_t nofSteps = (Int_t)(TMath::Sqrt(xvec*xvec+yvec*yvec+zvec*zvec)/fStep+1);
354
355 Double_t stepEL = fEnergyLossToSignal*point->GetEnergyLoss()/(nofSteps+1);
356
357 xvec = xvec/((Double_t)nofSteps);
358 yvec = yvec/((Double_t)nofSteps);
359 zvec = zvec/((Double_t)nofSteps);
360
361 for ( Int_t istep = 0 ; istep <= nofSteps ; istep++ ) {
362 // Float_t iFChan = sensor->GetChannelPlus(xin, yin, 0);
363 Int_t iIChan = sensor->GetFrontChannel(xin,yin,zin);
364 //(Int_t)iFChan;
365
366 if ( iIChan != -1 ) {
367 fStripSignalF[iIChan] += stepEL;
368 fFChannelPointsMap[iIChan].insert(iPoint);
369 }
370
371 // iFChan = sensor->GetChannelPlus(xin, yin, 1);
372 // iIChan = (Int_t)iFChan;
373 iIChan = sensor->GetBackChannel (xin,yin,zin);
374
375 if ( iIChan != -1 ) {
376 fStripSignalB[iIChan] += stepEL;
377 fBChannelPointsMap[iIChan].insert(iPoint);
378 }
379
380 xin+=xvec;
381 yin+=yvec;
382 zin+=zvec;
383 }
384
385 }
386
387}
388// -------------------------------------------------------------------------
389
390// ----- Private method FindFiredStrips --------------------------------
391void CbmStsRealDigitize::FindFiredStrips(CbmStsPoint* pnt,Int_t& nofStr,Int_t*& strips,Double_t*& signals, Int_t side) {
392
393 nofStr = 0;
394
395 Double_t xin = pnt->GetXIn();
396 Double_t yin = pnt->GetYIn();
397 Double_t zin = pnt->GetZIn();
398
399 gGeoManager->FindNode(xin,yin,zin);
400 TGeoNode* curNode = gGeoManager->GetCurrentNode();
401
402 CbmStsSensor* sensor = fDigiScheme->GetSensorByName(curNode->GetName());
403
404 Double_t xvec = pnt->GetXOut()-xin;
405 Double_t yvec = pnt->GetYOut()-yin;
406 Double_t zvec = pnt->GetZOut()-zin;
407
408 Int_t nofSteps = (Int_t)(TMath::Sqrt(xvec*xvec+yvec*yvec+zvec*zvec)/fStep+1);
409
410 Double_t stepEL = fEnergyLossToSignal*pnt->GetEnergyLoss()/(nofSteps+1);
411
412 xvec = xvec/((Double_t)nofSteps);
413 yvec = yvec/((Double_t)nofSteps);
414 zvec = zvec/((Double_t)nofSteps);
415
416 for ( Int_t istep = 0 ; istep <= nofSteps ; istep++ ) {
417 Float_t iFChan = sensor->GetChannelPlus(xin, yin, side);
418 Int_t iIChan = (Int_t)iFChan;
419
420 xin+=xvec;
421 yin+=yvec;
422 zin+=zvec;
423
424 if ( iIChan == -1 ) continue;
425
426 for ( Int_t ifstr = 0 ; ifstr < nofStr ; ifstr++ ) {
427 if ( strips[ifstr] == iIChan ) {
428 signals[ifstr] += stepEL;
429 iIChan = -1;
430 break;
431 }
432 }
433 if ( iIChan == -1 ) continue;
434
435 strips [nofStr] = iIChan;
436 signals[nofStr] = stepEL;
437
438 nofStr++;
439
440 }
441}
442// -------------------------------------------------------------------------
443
444// ----- Private method SetParContainers -------------------------------
445void CbmStsRealDigitize::SetParContainers() {
446
447 // Get run and runtime database
448 FairRunAna* run = FairRunAna::Instance();
449 if ( ! run ) Fatal("SetParContainers", "No analysis run");
450
451 FairRuntimeDb* db = run->GetRuntimeDb();
452 if ( ! db ) Fatal("SetParContainers", "No runtime database");
453
454 // Get STS geometry parameter container
455 fGeoPar = (CbmGeoStsPar*) db->getContainer("CbmGeoStsPar");
456
457 // Get STS digitisation parameter container
458 fDigiPar = (CbmStsDigiPar*) db->getContainer("CbmStsDigiPar");
459
460}
461// -------------------------------------------------------------------------
462
463
464
465// ----- Private method Init -------------------------------------------
466InitStatus CbmStsRealDigitize::Init() {
467
468 // Get input array
469 FairRootManager* ioman = FairRootManager::Instance();
470 if ( ! ioman ) Fatal("Init", "No FairRootManager");
471 fPoints = (TClonesArray*) ioman->GetObject("StsPoint");
472
473 // Register output array StsDigi
474 fDigis = new TClonesArray("CbmStsDigi",1000);
475 ioman->Register("StsDigi", "Digital response in STS", fDigis, kTRUE);
476
477 // Register output array StsDigiMatches
478 fDigiMatches = new TClonesArray("CbmStsDigiMatch",1000);
479 ioman->Register("StsDigiMatch", "Digi Match in STS", fDigiMatches, kTRUE);
480
481// fGen = new TRandom3();
482// time_t curtime;
483// time(&curtime);
484// fGen->SetSeed(curtime);
485
486 fStripSignalF = new Double_t[2000];
487 fStripSignalB = new Double_t[2000];
488
489 fEnergyLossToSignal = 200000.;
490
491 fFNofSteps = (Int_t)TMath::Power(2,(Double_t)fFNofBits);
492 fBNofSteps = (Int_t)TMath::Power(2,(Double_t)fBNofBits);
493
494 // Build digitisation scheme
495 if ( fDigiScheme->Init(fGeoPar, fDigiPar) ) {
496 MakeSets();
497
498 if (fVerbose == 1 || fVerbose == 2) fDigiScheme->Print(kFALSE);
499 else if (fVerbose > 2) fDigiScheme->Print(kTRUE);
500 cout << "-I- " << fName << "::Init: "
501 << "STS digitisation scheme succesfully initialised" << endl;
502 cout << " Stations: " << fDigiScheme->GetNStations()
503 << ", Sectors: " << fDigiScheme->GetNSectors() << ", Channels: "
504 << fDigiScheme->GetNChannels() << endl;
505 return kSUCCESS;
506 }
507
508 return kERROR;
509
510}
511// -------------------------------------------------------------------------
512
513
514
515// ----- Private method ReInit -----------------------------------------
516InitStatus CbmStsRealDigitize::ReInit() {
517
518 // Clear digitisation scheme
519 fDigiScheme->Clear();
520
521 // Build new digitisation scheme
522 if ( fDigiScheme->Init(fGeoPar, fDigiPar) ) {
523 MakeSets();
524 return kSUCCESS;
525 }
526
527 return kERROR;
528
529}
530// -------------------------------------------------------------------------
531
532
533// ----- Private method MakeSets ---------------------------------------
534void CbmStsRealDigitize::MakeSets() {
535
536 fPointMap.clear();
537 Int_t nStations = fDigiScheme->GetNStations();
538 for (Int_t iStation=0; iStation<nStations; iStation++) {
539 CbmStsStation* station = fDigiScheme->GetStation(iStation);
540 Int_t nSectors = station->GetNSectors();
541 for (Int_t iSector=0; iSector<nSectors; iSector++) {
542 CbmStsSector* sector = station->GetSector(iSector);
543 Int_t nSensors = sector->GetNSensors();
544 for (Int_t iSensor=0; iSensor<nSensors; iSensor++) {
545 CbmStsSensor* sensor = sector->GetSensor(iSensor);
546 set<Int_t> a;
547 fPointMap[sensor] = a;
548 }
549 }
550 }
551 fFChannelPointsMap.clear();
552 fBChannelPointsMap.clear();
553 for ( Int_t ichan = 2000 ; ichan > 0 ; ) {
554 set<Int_t> a;
555 fFChannelPointsMap[--ichan] = a;
556 set<Int_t> b;
557 fBChannelPointsMap[ ichan] = b;
558 }
559}
560// -------------------------------------------------------------------------
561
562
563// ----- Private method Reset ------------------------------------------
564void CbmStsRealDigitize::Reset() {
565 // fNPoints = fNOutside = fNDigisFront = fNDigisBack = fTime = 0.;
566 fNDigis = fNMulti = 0;
567 fFChannelPointsMap.clear();
568 fBChannelPointsMap.clear();
569 // if ( fDigis ) fDigis->Clear();
570 // if ( fDigiMatches ) fDigiMatches->Clear();
571 if ( fDigis ) fDigis->Delete();
572 if ( fDigiMatches ) fDigiMatches->Delete();
573}
574// -------------------------------------------------------------------------
575
576
577// ----- Virtual method Finish -----------------------------------------
579 cout << endl;
580 cout << "============================================================"
581 << endl;
582 cout << "===== " << fName << ": Run summary " << endl;
583 cout << "===== " << endl;
584 cout << "===== Events processed : " << setw(8) << fNEvents << endl;
585 cout.setf(ios_base::fixed, ios_base::floatfield);
586 cout << "===== Real time per event : "
587 << setw(8) << setprecision(4)
588 << fTime / fNEvents << " s" << endl;
589 cout << "===== StsPoints per event : "
590 << setw(8) << setprecision(2)
591 << fNPoints / fNEvents << endl;
592 cout << "===== Outside hits per event : "
593 << setw(8) << setprecision(2)
594 << fNOutside / fNEvents << " = "
595 << setw(6) << setprecision(2)
596 << fNOutside / fNPoints * 100. << " %" << endl;
597 cout << "===== Front digis per point : "
598 << setw(8) << setprecision(2)
599 << fNDigisFront / (fNPoints-fNOutside) << endl;
600 cout << "===== Back digis per point : "
601 << setw(8) << setprecision(2)
602 << fNDigisBack / (fNPoints-fNOutside) << endl;
603 cout << "============================================================"
604 << endl;
605
606}
607// -------------------------------------------------------------------------
void Print(Bool_t kLong=kFALSE)
CbmStsSensor * GetSensorByName(TString sensorName)
CbmStsStation * GetStation(Int_t iStation)
Double_t GetXIn() const
Definition CbmStsPoint.h:69
Double_t GetZOut() const
Definition CbmStsPoint.h:74
Double_t GetYOut() const
Definition CbmStsPoint.h:73
Double_t GetZIn() const
Definition CbmStsPoint.h:71
Double_t GetYIn() const
Definition CbmStsPoint.h:70
Double_t GetXOut() const
Definition CbmStsPoint.h:72
void ProduceHitResponse(CbmStsSensor *sensor)
virtual void Exec(Option_t *opt)
void FindFiredStrips(CbmStsPoint *pnt, Int_t &nofStr, Int_t *&strips, Double_t *&signals, Int_t side)
Int_t GetSectorNr() const
Int_t GetStationNr() const
Int_t GetDetectorId() const
Int_t GetNChannelsFront() const
CbmStsSensor * GetSensor(Int_t iSensor)
Int_t GetNSensors() const
Float_t GetChannelPlus(Double_t x, Double_t y, Int_t iSide)
Int_t GetFrontChannel(Double_t x, Double_t y, Double_t z, Double_t &dPitch)
Int_t GetBackChannel(Double_t x, Double_t y, Double_t z, Double_t &dPitch)
Int_t GetNSectors() const
CbmStsSector * GetSector(Int_t iSector)
-clang-format