BmnRoot
Loading...
Searching...
No Matches
CbmStsDigitize.cxx
Go to the documentation of this file.
1// -------------------------------------------------------------------------
2// ----- CbmStsDigitize source file -----
3// ----- Created 08/07/2008 by R. Karabowicz -----
4// -------------------------------------------------------------------------
5
6// Includes from ROOT
7#include "TClonesArray.h"
8#include "TF1.h"
9#include "TGeoBBox.h"
10#include "TGeoManager.h"
11#include "TGeoNode.h"
12#include "TMath.h"
13#include "TObjArray.h"
14#include "TRandom3.h"
15
16#include <TFile.h> //AZ-290322
17
18// Includes from base
19#include "FairField.h" //AZ
20#include "FairRootManager.h"
21#include "FairRunAna.h"
22#include "FairRuntimeDb.h"
23
24// Includes from STS
25#include "BmnGemFastDigitize.h" //AZ
26#include "BmnSiliconPoint.h" //AZ-280322
27#include "CbmGeoStsPar.h"
28#include "CbmMCTrack.h" //AZ
29#include "CbmStsDigi.h"
30#include "CbmStsDigiMatch.h"
31#include "CbmStsDigiPar.h"
32#include "CbmStsDigiScheme.h"
33#include "CbmStsDigitize.h"
34#include "CbmStsPoint.h"
35#include "CbmStsSector.h"
36#include "CbmStsSensor.h"
37#include "CbmStsStation.h"
38
39#include <TDatabasePDG.h>
40#include <iomanip>
41#include <iostream>
42#include <map>
43#include <vector>
44using std::cerr;
45using std::cout;
46using std::endl;
47using std::fixed;
48using std::flush;
49using std::ios_base;
50using std::left;
51using std::map;
52using std::pair;
53using std::right;
54using std::set;
55using std::setprecision;
56using std::setw;
57using std::vector;
58
59// ----- Default constructor ------------------------------------------
60CbmStsDigitize::CbmStsDigitize()
61 : FairTask("STS Digitizer", 1)
62 , fGeoPar(NULL)
63 , fDigiPar(NULL)
64 , fDigiScheme(NULL)
65 , fPointsSi(nullptr)
66 , // AZ-280322
67 fPoints(NULL)
68 , fDigis(NULL)
69 , fDigiMatches(NULL)
70 , fNDigis(0)
71 , fNMulti(0)
72 , fNEvents(0)
73 , fNPoints(0)
74 , fNDigisFront(0)
75 , fNDigisBack(0)
76 , fTime(0.)
77 , fStep(0.001)
78 , fTimer()
79 , fRealistic(kFALSE)
80 , fEnergyLossToSignal(0.)
81 , fFThreshold(4.0)
82 , fBThreshold(4.0)
83 , fFNoiseWidth(0.1)
84 , fBNoiseWidth(0.1)
85 , fStripDeadTime(10)
86 , fFNofBits(11)
87 , fBNofBits(11)
88 , fFNofElPerAdc(10.)
89 , fBNofElPerAdc(10.)
90 , fFNofSteps(0)
91 , fBNofSteps(0)
92 , fStripSignalF(NULL)
93 , fStripSignalB(NULL)
94 , fFChannelPointsMap()
95 , fBChannelPointsMap()
96 , fPointMap()
97 , fastDig(NULL) // AZ
98{
99 fDigiScheme = CbmStsDigiScheme::Instance();
100 Reset();
101}
102// -------------------------------------------------------------------------
103
104// ----- Standard constructor ------------------------------------------
105CbmStsDigitize::CbmStsDigitize(Int_t iVerbose)
106 : FairTask("STSDigitize", iVerbose)
107 , fGeoPar(NULL)
108 , fDigiPar(NULL)
109 , fDigiScheme(NULL)
110 , fPointsSi(nullptr)
111 , // AZ-280322
112 fPoints(NULL)
113 , fDigis(NULL)
114 , fDigiMatches(NULL)
115 , fNDigis(0)
116 , fNMulti(0)
117 , fNEvents(0)
118 , fNPoints(0)
119 , fNDigisFront(0)
120 , fNDigisBack(0)
121 , fTime(0.)
122 , fStep(0.001)
123 , fTimer()
124 , fRealistic(kFALSE)
125 , fEnergyLossToSignal(0.)
126 , fFThreshold(4.0)
127 , fBThreshold(4.0)
128 , fFNoiseWidth(0.1)
129 , fBNoiseWidth(0.1)
130 , fStripDeadTime(10)
131 , fFNofBits(11)
132 , fBNofBits(11)
133 , fFNofElPerAdc(10.)
134 , fBNofElPerAdc(10.)
135 , fFNofSteps(0)
136 , fBNofSteps(0)
137 , fStripSignalF(NULL)
138 , fStripSignalB(NULL)
139 , fFChannelPointsMap()
140 , fBChannelPointsMap()
141 , fPointMap()
142 , fastDig(NULL) // AZ
143{
144 fDigiScheme = CbmStsDigiScheme::Instance();
145 Reset();
146}
147// -------------------------------------------------------------------------
148
149// ----- Constructor with name -----------------------------------------
150CbmStsDigitize::CbmStsDigitize(const char* name, Int_t iVerbose)
151 : FairTask(name, iVerbose)
152 , fGeoPar(NULL)
153 , fDigiPar(NULL)
154 , fDigiScheme(NULL)
155 , fPointsSi(nullptr)
156 , // AZ-280322
157 fPoints(NULL)
158 , fDigis(NULL)
159 , fDigiMatches(NULL)
160 , fNDigis(0)
161 , fNMulti(0)
162 , fNEvents(0)
163 , fNPoints(0)
164 , fNDigisFront(0)
165 , fNDigisBack(0)
166 , fTime(0.)
167 , fStep(0.001)
168 , fTimer()
169 , fRealistic(kFALSE)
170 , fEnergyLossToSignal(0.)
171 , fFThreshold(4.0)
172 , fBThreshold(4.0)
173 , fFNoiseWidth(0.1)
174 , fBNoiseWidth(0.1)
175 , fStripDeadTime(10)
176 , fFNofBits(11)
177 , fBNofBits(11)
178 , fFNofElPerAdc(10.)
179 , fBNofElPerAdc(10.)
180 , fFNofSteps(0)
181 , fBNofSteps(0)
182 , fStripSignalF(NULL)
183 , fStripSignalB(NULL)
184 , fFChannelPointsMap()
185 , fBChannelPointsMap()
186 , fPointMap()
187 , fastDig(NULL) // AZ
188{
189 fGeoPar = NULL;
190 fDigiPar = NULL;
191 fPoints = NULL;
192 fDigis = NULL;
193 fDigiMatches = NULL;
194 fRealistic = kFALSE;
195 fDigiScheme = CbmStsDigiScheme::Instance();
196 Reset();
197
198 fStep = 0.001;
199
200 fFThreshold = 4.0;
201 fBThreshold = 4.0;
202 fFNoiseWidth = 0.1;
203 fBNoiseWidth = 0.1;
204
205 fFNofBits = 11;
206 fBNofBits = 11;
207 fFNofElPerAdc = 10.;
208 fBNofElPerAdc = 10.;
209 fStripDeadTime = 10;
210 fNEvents = 0.;
211}
212// -------------------------------------------------------------------------
213
214// ----- Destructor ----------------------------------------------------
215CbmStsDigitize::~CbmStsDigitize()
216{
217 if (fGeoPar)
218 delete fGeoPar;
219 if (fDigiPar)
220 delete fDigiPar;
221 if (fDigis) {
222 fDigis->Delete();
223 delete fDigis;
224 }
225 if (fDigiMatches) {
226 fDigiMatches->Delete();
227 delete fDigiMatches;
228 }
229 // if ( fDigiScheme ) delete fDigiScheme; // M&SG HEPL: Deleting singleton is safe at the end of a whole program
230 // only!
231 Reset();
232}
233// -------------------------------------------------------------------------
234
235// ----- Public method Exec --------------------------------------------
236void CbmStsDigitize::Exec(Option_t* opt)
237{
238
239 // Reset all eventwise counters
240 fTimer.Start();
241 Reset();
242 // Verbose screen output
243 if (fVerbose > 2) {
244 cout << endl << "-I- " << fName << ": executing event" << endl;
245 cout << "----------------------------------------------" << endl;
246 }
247
248 Int_t nPoints = 0;
249 Int_t nDigisF = 0;
250 Int_t nDigisB = 0;
251
252 // Loop over all StsPoints
253 if (!fPoints) {
254 cerr << "-W- " << fName << "::Exec: No input array (STSPoint) " << endl;
255 cout << "- " << fName << endl;
256 return;
257 }
258 if (fFNofBits > CbmStsDigi::GetNofAdcBits() || fBNofBits > CbmStsDigi::GetNofAdcBits()) {
259 cerr << "-W- " << fName << "::Exec: Number of AdcBits(" << fFNofBits
260 << ") during digitization exceeds ADC range(" << CbmStsDigi::GetNofAdcBits() << ") defined in data class "
261 << endl;
262 cout << "- " << fName << endl;
263 return;
264 }
265 // AZ if (fFNofBits<=CbmStsDigi::GetNofAdcBits()&& fFNofElPerAdc!=10*TMath::Power(2,(12-fFNofBits))) {
266 if (0 && fFNofBits <= CbmStsDigi::GetNofAdcBits() && fFNofElPerAdc != 10 * TMath::Power(2, (12 - fFNofBits))) {
267 cerr << "-W- " << fName << "::Exec: Number of electrons per AdcChannel does not match Adc range " << endl;
268 cout << "- " << fName << endl;
269 return;
270 }
271 // AZ if (fBNofBits<=CbmStsDigi::GetNofAdcBits()&& fBNofElPerAdc!=10*TMath::Power(2,(12-fFNofBits))) {
272 if (0 && fBNofBits <= CbmStsDigi::GetNofAdcBits() && fBNofElPerAdc != 10 * TMath::Power(2, (12 - fFNofBits))) {
273 cerr << "-W- " << fName << "::Exec: Number of electrons per AdcChannel does not match Adc range " << endl;
274 cout << "- " << fName << endl;
275 return;
276 }
277 map<CbmStsSensor*, set<Int_t>>::iterator mapIt;
278 for (mapIt = fPointMap.begin(); mapIt != fPointMap.end(); mapIt++)
279 ((*mapIt).second).clear();
280
281 for (int idet = 0; idet < 2; ++idet) { // AZ-280322
282 TClonesArray* points = (idet == 0) ? fPointsSi : fPoints;
283 if (points == nullptr)
284 continue;
285
286 // AZ-280322 for (Int_t iPoint=0; iPoint<fPoints->GetEntriesFast(); iPoint++) {
287 for (Int_t iPoint = 0; iPoint < points->GetEntriesFast(); iPoint++) {
288 // AZ-280322 CbmStsPoint* point = (CbmStsPoint*) fPoints->At(iPoint);
289 CbmStsPoint* point = (CbmStsPoint*)points->At(iPoint);
290 /*AZ!!! Exclude some points !!!
291 CbmMCTrack *mcTr = (CbmMCTrack*) fMCTracks->UncheckedAt(point->GetTrackID());
292 if (TMath::Abs(mcTr->GetPdgCode()) != 211) continue;
293 TVector3 mom;
294 mcTr->GetMomentum(mom);
295 if (mom.Eta() > 2) continue;
296 //AZ
297 */
298 // if (point->GetZIn() > 100.0) continue; //!!!AZ - exclude GEMs (010920)
299
300 Double_t xin = point->GetXIn();
301 Double_t yin = point->GetYIn();
302 Double_t zin = point->GetZIn();
303 gGeoManager->FindNode(xin, yin, zin);
304 TGeoNode* curNode = gGeoManager->GetCurrentNode();
305 if (TString(curNode->GetName()).Contains("part")) {
306 // AZ - Test beam silicon
307 // cout << curNode->GetName() << " " << gGeoManager->GetPath() << endl;
308 gGeoManager->CdUp();
309 curNode = gGeoManager->GetCurrentNode();
310 // cout << curNode->GetName() << " " << gGeoManager->GetPath() << endl;
311 // exit(0);
312 }
313
314 CbmStsSensor* sensor = NULL;
315 if (fDigiScheme->IsNewGeometry()) {
316 TString curPath = fDigiScheme->GetCurrentPath();
317 // if (!curPath.Contains("Si") && CrossSpacer(curNode, point)) continue; //AZ - emulate spacers
318 sensor = fDigiScheme->GetSensorByName(curPath);
319 } else
320 sensor = fDigiScheme->GetSensorByName(curNode->GetName());
321
322 if (fPointMap.find(sensor) == fPointMap.end()) {
323 cerr << "-E- " << fName << "::Exec:: sensor " << curNode->GetName() << " not found in digi scheme!"
324 << endl;
325 continue;
326 }
327 fPointMap[sensor].insert(iPoint);
328 nPoints++;
329 }
330 } // for (int idet = 0; idet < 2;
331
332 for (Int_t iStation = fDigiScheme->GetNStations(); iStation > 0;) {
333 CbmStsStation* station = fDigiScheme->GetStation(--iStation);
334 for (Int_t iSector = station->GetNSectors(); iSector > 0;) {
335 CbmStsSector* sector = station->GetSector(--iSector);
336
337 map<Int_t, set<Int_t>>::iterator mapCh;
338
339 for (mapCh = fFChannelPointsMap.begin(); mapCh != fFChannelPointsMap.end(); mapCh++)
340 ((*mapCh).second).clear();
341 for (mapCh = fBChannelPointsMap.begin(); mapCh != fBChannelPointsMap.end(); mapCh++)
342 ((*mapCh).second).clear();
343
344 // simulating detector+cables+electronics noise
345 // should be more sophisticated...
346 // the question is: sectorwise or sensorwise???
347 Int_t nChannels = sector->GetNChannelsFront();
348
349 //-----aplying noise on every channel-----
350 for (Int_t iChannel = nChannels; iChannel > 0;) {
351 // fStripSignalF[--iChannel] = fGen->Landau(.1,.02);
352 // fStripSignalB[ iChannel] = fGen->Landau(.1,.02);
353 // fStripSignalF[--iChannel] = 0.;
354 // fStripSignalB[ iChannel] = 0.;
355 // AZ fStripSignalF[--iChannel] = TMath::Abs(gRandom->Gaus(0.,fFNoiseWidth));
356 fStripSignalF[--iChannel] = 0;
357 // AZ fStripSignalB[ iChannel] = TMath::Abs(gRandom->Gaus(0.,fBNoiseWidth));
358 }
359 // AZ
360 Int_t nChannelsB = sector->GetNChannelsBack();
361 //-----aplying noise on every channel-----
362 for (Int_t iChannel = nChannelsB; iChannel > 0;) {
363 // AZ fStripSignalB[--iChannel] = TMath::Abs(gRandom->Gaus(0.,fBNoiseWidth));
364 fStripSignalB[--iChannel] = 0;
365 }
366 // AZ
367 // if (sector->GetStationNr() > 3) cout << " yyyy " << sector->GetStationNr() << " " <<
368 // sector->GetSectorNr() << endl; //AZ-290322
369
370 for (Int_t iSensor = sector->GetNSensors(); iSensor > 0;) {
371 CbmStsSensor* sensor = sector->GetSensor(--iSensor);
372
373 if (sensor->GetDx() < 0.02)
374 ProduceHitResponseSi(sensor);
375 else if (fastDig) {
376 if (fPointMap.find(sensor) != fPointMap.end())
377 ((BmnGemFastDigitize*)fastDig)
378 ->ProduceHitResponseFast(sensor, fPointMap[sensor], fStripSignalF, fStripSignalB,
379 fFChannelPointsMap, fBChannelPointsMap);
380 } else
381 ProduceHitResponse(sensor);
382 // ProduceHitResponseAZ(sensor); //AZ
383 }
384
385 Int_t stationNr = sector->GetStationNr();
386 Int_t sectorNr = sector->GetSectorNr();
387 // Int_t sectorDetId = sector->GetDetectorId();
388 Double_t dx = sector->GetDx();
389
390 for (Int_t ifstr = 0; ifstr < nChannels; ifstr++) {
391 if (fStripSignalF[ifstr] < 1)
392 continue; // AZ
393 // AZ if ( fStripSignalF[ifstr] < (fFThreshold*1000.) ) continue;//threshold cut
394 // if (fastDig && dx > 0.02) {
395 // if (fStripSignalF[ifstr] < 13.0*1) continue;
396 // } else if ( fStripSignalF[ifstr] < (fFThreshold*1000.) ) continue;//threshold cut
397
398 //-----random strip inefficiency-----
399 // Double_t generator;
400 // generator = gRandom->Rndm()*100.;
401 // AZ if (generator< (fStripDeadTime/100.)*occupancy [iStation][iSector][ifstr/125]) continue;
402 //-----------------------------------
403
404 Int_t digiFSignal = 0;
405 if (fastDig && dx > 0.02)
406 digiFSignal = fStripSignalF[ifstr] / 1.;
407 else
408 digiFSignal = 1 + (Int_t)((fStripSignalF[ifstr]) / fFNofElPerAdc);
409 if (digiFSignal >= fFNofSteps)
410 digiFSignal = fFNofSteps - 1;
411 if (digiFSignal < 13)
412 continue; // AZ
413 // if (digiFSignal < 1) continue; //AZ - for cluster library
414 // if (digiFSignal < 5) continue; //AZ - for cluster library
415 new ((*fDigis)[fNDigis]) CbmStsDigi(stationNr, sectorNr, 0, ifstr, digiFSignal, 0);
416 set<Int_t>::iterator it1;
417 set<Int_t> chPnt = fFChannelPointsMap[ifstr];
418 Int_t pnt;
419 CbmStsDigiMatch* match;
420 if (chPnt.size() == 0) {
421 new ((*fDigiMatches)[fNDigis]) CbmStsDigiMatch(-666);
422 } else {
423 for (it1 = chPnt.begin(); it1 != chPnt.end(); it1++) {
424 pnt = (*it1);
425 if (it1 == chPnt.begin())
426 match = new ((*fDigiMatches)[fNDigis]) CbmStsDigiMatch(pnt);
427 else {
428 match->AddPoint(pnt);
429 fNMulti++;
430 }
431 }
432 }
433 fNDigis++;
434 nDigisF++;
435 }
436
437 nChannels = sector->GetNChannelsBack(); // AZ
438
439 for (Int_t ibstr = 0; ibstr < nChannels; ibstr++) {
440 if (fStripSignalB[ibstr] < 1)
441 continue; // AZ
442 // AZ if ( fStripSignalB[ibstr] < (fBThreshold*1000.) ) continue; //threshold cut
443 // if (fastDig && dx > 0.02) {
444 // if (fStripSignalB[ibstr] < 13.0*1) continue;
445 // } else if ( fStripSignalB[ibstr] < (fBThreshold*1000.) ) continue; //threshold cut
446
447 //-----random strip inefficiency-----
448 // Double_t generator;
449 // generator = gRandom->Rndm()*100.;
450 // AZ if (generator< (fStripDeadTime/100.)*occupancy [iStation][iSector][ibstr/125]) continue;
451 //-----------------------------------
452
453 Int_t digiBSignal = 0;
454 if (fastDig && dx > 0.02)
455 digiBSignal = fStripSignalB[ibstr] / 1.;
456 else
457 digiBSignal = 1 + (Int_t)((fStripSignalB[ibstr]) / fBNofElPerAdc);
458 if (digiBSignal >= fBNofSteps)
459 digiBSignal = fBNofSteps - 1;
460 if (digiBSignal < 13)
461 continue; // AZ
462 // if (digiBSignal < 1) continue; //AZ - for cluster library
463 // if (digiBSignal < 5) continue; //AZ - for cluster library
464 new ((*fDigis)[fNDigis]) CbmStsDigi(stationNr, sectorNr, 1, ibstr, digiBSignal, 0);
465 set<Int_t>::iterator it1;
466 set<Int_t> chPnt = fBChannelPointsMap[ibstr];
467 Int_t pnt;
468 CbmStsDigiMatch* match;
469 if (chPnt.size() == 0) {
470 new ((*fDigiMatches)[fNDigis]) CbmStsDigiMatch(-666);
471 } else {
472 for (it1 = chPnt.begin(); it1 != chPnt.end(); it1++) {
473 pnt = (*it1);
474 if (it1 == chPnt.begin())
475 match = new ((*fDigiMatches)[fNDigis]) CbmStsDigiMatch(pnt);
476 else {
477 match->AddPoint(pnt);
478 fNMulti++;
479 }
480 }
481 }
482 fNDigis++;
483 nDigisB++;
484 }
485 }
486 }
487
488 fTimer.Stop();
489 cout << "+ " << flush;
490 cout << setw(15) << left << fName << ": " << setprecision(4) << setw(8) << fixed << right << fTimer.RealTime()
491 << " s, digis " << nDigisF << " / " << nDigisB << endl;
492
493 fNEvents += 1.;
494 fNPoints += Double_t(nPoints);
495 fNDigisFront += Double_t(nDigisF);
496 fNDigisBack += Double_t(nDigisB);
497 fTime += fTimer.RealTime();
498}
499// -------------------------------------------------------------------------
500
501// ----- Private method CrossSpacer ------------------------------------
502Bool_t CbmStsDigitize::CrossSpacer(const TGeoNode* node, const CbmStsPoint* point)
503{
504 // AZ - Check if particle goes thru the spacer (dead space)
505
506 TString name = node->GetName();
507 // if (name.Contains("SensorSV")) return kFALSE; // short strip area
508
509 Double_t width = 0.1, pitch = 10.0; // spacer width/2 and pitch
510 Double_t xyzloc[3], xyz[3] = {(point->GetXIn() + point->GetXOut()) / 2, (point->GetYIn() + point->GetYOut()) / 2,
511 (point->GetZIn() + point->GetZOut()) / 2};
512 gGeoManager->MasterToLocal(xyz, xyzloc);
513
514 TGeoVolume* vol = node->GetVolume();
515 TGeoBBox* box = (TGeoBBox*)vol->GetShape();
516 Double_t dx = box->GetDX(), dy = box->GetDY();
517 gGeoManager->CdUp();
518 TGeoNode* nodeM = gGeoManager->GetCurrentNode();
519 const Double_t* transl = nodeM->GetMatrix()->GetTranslation();
520
521 // X-spacers: check distance from the beam pipe corner
522 Double_t dist = 0;
523 if (transl[0] > -1)
524 dist = xyzloc[0] + dx;
525 else
526 dist = dx - xyzloc[0];
527 Double_t dscaled = dist / pitch;
528 Int_t nspacer = TMath::Nint(dscaled);
529 if (nspacer == 0)
530 return kFALSE; // at the border
531 if (TMath::Abs(dist - nspacer * pitch) < width)
532 return kTRUE;
533
534 // Y-spacers: check distance from the beam pipe corner
535 if (transl[1] > 0)
536 dist = xyzloc[1] + dy;
537 else
538 dist = dy - xyzloc[1];
539 dscaled = dist / pitch;
540 nspacer = TMath::Nint(dscaled);
541 if (nspacer == 0)
542 return kFALSE; // at the border
543 if (TMath::Abs(dist - nspacer * pitch) < width)
544 return kTRUE; // Y-spacer
545
546 return kFALSE;
547}
548// -------------------------------------------------------------------------
549
550// ----- Private method ProduceHitResponse --------------------------------
551void CbmStsDigitize::ProduceHitResponseSi(CbmStsSensor* sensor)
552{
553 // Produce response in Silicon
554
555 // AZ set <Int_t> pSet;
556 if (fPointMap.find(sensor) == fPointMap.end()) {
557 cerr << "-E- " << fName << "::ProduceHitResponse:: sensor"
558 << " not found in digi scheme!" << endl;
559 return;
560 }
561
562 // AZ pSet = fPointMap[sensor];
563 set<Int_t>& pSet = fPointMap[sensor];
564
565 Int_t iPoint = -1;
566 CbmStsPoint* point = NULL;
567 BmnSiliconPoint* pointsi = nullptr; // AZ-280322
568
569 set<Int_t>::iterator it1;
570
571 Double_t dPitch, step = fStep; // AZ
572 if (TString(sensor->GetName()).Contains("Si") || sensor->GetDx() < 0.0200)
573 step = 0.001; // AZ - silicon
574 TClonesArray* points = (fPointsSi == nullptr) ? fPoints : fPointsSi; // AZ-280322
575 TVector3 posIn, posOut; // AZ-280322
576
577 for (it1 = pSet.begin(); it1 != pSet.end(); it1++) {
578 iPoint = (*it1);
579 // AZ-280322 point = (CbmStsPoint*) fPoints->At(iPoint);
580 if (fPointsSi) {
581 pointsi = (BmnSiliconPoint*)points->UncheckedAt(iPoint);
582 pointsi->PositionIn(posIn);
583 pointsi->PositionOut(posOut);
584 } else {
585 point = (CbmStsPoint*)points->At(iPoint);
586 point->PositionIn(posIn);
587 point->PositionOut(posOut);
588 }
589
590 /*AZ-280322
591 Double_t xin = point->GetXIn();
592 Double_t yin = point->GetYIn();
593 Double_t zin = point->GetZIn();
594
595 Double_t xvec = point->GetXOut()-xin;
596 Double_t yvec = point->GetYOut()-yin;
597 Double_t zvec = point->GetZOut()-zin;
598 cout << xin << " " << yin << " " << zin << " " << point->GetXOut() << " " << yvec << " " << point->GetZOut() <<
599 endl; //AZ-280322
600 */
601 posOut -= posIn; // AZ-280322
602 Double_t xin = posIn.X(), yin = posIn.Y(), zin = posIn.Z(); // AZ-280322
603 Double_t xvec = posOut.X(), yvec = posOut.Y(), zvec = posOut.Z(); // AZ-280322
604
605 // Int_t nofSteps = (Int_t)(TMath::Sqrt(xvec*xvec+yvec*yvec+zvec*zvec)/fStep+1);
606 Int_t nofSteps = (Int_t)(TMath::Sqrt(xvec * xvec + yvec * yvec + zvec * zvec) / step + 1);
607
608 // AZ-280322 Double_t stepEL = fEnergyLossToSignal*point->GetEnergyLoss()/(nofSteps+1);
609 Double_t eloss = (fPointsSi == nullptr) ? point->GetEnergyLoss() : pointsi->GetEnergyLoss(); // AZ-280322
610 Double_t stepEL = fEnergyLossToSignal * eloss / (nofSteps + 1); // AZ-280322
611 // Double_t stepEL = 0.0; //AZ
612 // if (sensor->GetDx() < 0.0200) stepEL = 280000000 * 1.5 * point->GetEnergyLoss() / (nofSteps+1); //AZ - Si
613 // else stepEL = fEnergyLossToSignal*point->GetEnergyLoss()/(nofSteps+1); //AZ
614
615 xvec = xvec / ((Double_t)nofSteps);
616 yvec = yvec / ((Double_t)nofSteps);
617 zvec = zvec / ((Double_t)nofSteps);
618
619 for (Int_t istep = 0; istep <= nofSteps; istep++) {
620
621 // AZ - introduce smearing
622 Double_t xinS = xin, yinS = yin, sigma = 0.05; // 500um
623 if (step < 0.0009) {
624 // Smearing only for GEMs
625 xinS += gRandom->Gaus(0.0, sigma);
626 yinS += gRandom->Gaus(0.0, sigma);
627 }
628 // AZ
629
630 // AZ Int_t iIChan = sensor->GetFrontChannel(xin,yin,zin);
631 Int_t iIChan = sensor->GetFrontChannel(xinS, yinS, zin, dPitch);
632
633 // AZ!!! Exclude 50% of channels for large stations
634 // if (sensor->GetNChannelsFront() > 900 && sensor->GetSectorNr() == 1 && iIChan > 500) iIChan = -1;
635 // else if (sensor->GetNChannelsFront() > 900 && sensor->GetSectorNr() == 3 && iIChan < 500) iIChan = -1;
636
637 if (iIChan != -1) {
638 fStripSignalF[iIChan] += stepEL;
639 fFChannelPointsMap[iIChan].insert(iPoint);
640 }
641
642 // AZ iIChan = sensor->GetBackChannel (xin,yin,zin);
643 iIChan = sensor->GetBackChannel(xinS, yinS, zin, dPitch);
644
645 // AZ!!! Exclude 50% of channels for large stations
646 // if (sensor->GetNChannelsFront() > 900 && sensor->GetSectorNr() == 1 && iIChan > 600) iIChan = -1;
647 // else if (sensor->GetNChannelsFront() > 900 && sensor->GetSectorNr() == 3 && iIChan < 600) iIChan = -1;
648
649 if (iIChan != -1) {
650 fStripSignalB[iIChan] += stepEL;
651 fBChannelPointsMap[iIChan].insert(iPoint);
652 }
653
654 xin += xvec;
655 yin += yvec;
656 zin += zvec;
657 }
658 }
659}
660// -------------------------------------------------------------------------
661
662Double_t CbmStsDigitize::GetNumberOfClusters(Double_t beta, Double_t gamma, Double_t charge, Double_t p0, Double_t p1)
663{ // ES
664
665 Double_t beta2 = beta * beta;
666 Double_t gamma2 = gamma * gamma;
667 Double_t val = p0 / beta2 * (p1 + TMath::Log(beta2 * gamma2) - beta2);
668 // AZ return val;
669 return TMath::Min(val * charge * charge, 20000.);
670}
671
672// ----- Private method ProduceHitResponse --------------------------------
673void CbmStsDigitize::ProduceHitResponse(CbmStsSensor* sensor)
674{
675 //
676
677 // AZ - get magnetic field (for digitizer parameters in GEMs)
678 static Int_t ifield = -1;
679 if (ifield < 0) {
680 FairRunAna* run = FairRunAna::Instance();
681 FairField* magField = run->GetField();
682 ifield = TMath::Nint(TMath::Abs(magField->GetBy(0, 0, 135.0)));
683 }
684
685 // AZ set <Int_t> pSet;
686
687 if (fPointMap.find(sensor) == fPointMap.end()) {
688 cerr << "-E- " << fName << "::ProduceHitResponse:: sensor"
689 << " not found in digi scheme!" << endl;
690 return;
691 }
692
693 // AZ pSet = fPointMap[sensor];
694 set<Int_t>& pSet = fPointMap[sensor]; // AZ
695
696 Int_t iPoint = -1;
697 CbmStsPoint* point = NULL;
698
699 set<Int_t>::iterator it1;
700
701 Double_t dPitch = fStep /*, dx = sensor->GetDx(), dy = sensor->GetDy()*/; // AZ
702
703 for (it1 = pSet.begin(); it1 != pSet.end(); it1++) {
704
705 // AZ Apply overall efficiency
706 // const Double_t eff0 = 0.65;
707 // if (gRandom->Rndm() > eff0) continue;
708
709 iPoint = (*it1);
710 // AZ point = (CbmStsPoint*)fPoints->At(iPoint);
711 point = (CbmStsPoint*)fPoints->UncheckedAt(iPoint); // AZ
712 TVector3 mom;
713 point->Momentum(mom);
714
715 CbmMCTrack* mcTr = (CbmMCTrack*)fMCTracks->UncheckedAt(point->GetTrackID()); // ES
716 TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(mcTr->GetPdgCode());
717 Double_t mnoc = 1.0, mass = 3.7283999; // He4
718 if (part)
719 mass = part->Mass();
720 if (mass > 0.0001) {
721 TLorentzVector lorv;
722 lorv.SetVectM(mom, mass);
723 Double_t beta = lorv.Beta();
724 Double_t gamma = lorv.Gamma();
725 Double_t charge = (part) ? part->Charge() / 3 : 1;
726 mnoc = GetNumberOfClusters(beta, gamma, charge, 1.787, 12.33); // ES
727 if (mnoc < 0.5)
728 continue; // AZ
729 // cout << " ***clus*** " << mass << " " << beta << " " << gamma << " " << mnoc << endl;
730 }
731
732 Double_t moduleThickness = 0.3 + 0.25 + 0.2 + 0.15; // DriftGapThickness + FirstTransferGapThickness +
733 // SecondTransferGapThickness + InductionGapThickness;
734 Int_t stationNr = Int_t(sensor->GetStationNr());
735
736 Double_t xin = point->GetXIn();
737 Double_t yin = point->GetYIn();
738 Double_t zin = point->GetZIn();
739
740 Double_t xvec = point->GetXOut() - xin;
741 Double_t yvec = point->GetYOut() - yin;
742 // AZ Double_t zvec;
743
744 // AZ if (point->GetPz() > 0.0) zvec = point->GetZOut()-zin;
745 // AZ else zvec = zin - point->GetZOut();
746 Double_t zvec = point->GetZOut() - zin; // AZ
747
748 // length of track
749 Double_t track_length = TMath::Sqrt(xvec * xvec + yvec * yvec + zvec * zvec);
750
751 // Energy loss step for cluster, for electron
752 // Double_t stepEL, stepEl_current;
753
754 // Exponential dependency on a distance between clusters
755 Double_t current_length = 0.0; // traversed length (counter)
756 Double_t current_length_ratio = 0.0; // ratio of the traversed length to common track length (counter)
757
758 Int_t collisions_cnt = 0; // total collision counter
759 Double_t current_step = 0.0; // current distance between two collision points
760
761 // Collection of collision points
762 std::vector<std::vector<Double_t>> collision_points;
764 Double_t mcd = 1 / mnoc; // ES
765
766 mcd = 0.0245; // AZ-190322 - test - pi+, p = 1 GeV/c
767
768 while (current_length < track_length) {
769
770 current_step = gRandom->Exp(mcd);
771 current_length += current_step;
772
773 if (current_length > track_length)
774 break;
775
776 current_length_ratio = current_length / track_length;
777
778 // In terms of distance
779 Double_t current_x = xin + current_length_ratio * xvec;
780 Double_t current_y = yin + current_length_ratio * yvec;
781 Double_t current_z = zin + current_length_ratio * zvec;
782
783 std::vector<Double_t> collPoint;
784 collPoint.push_back(current_x);
785 collPoint.push_back(current_y);
786 collPoint.push_back(current_z);
787 collision_points.push_back(collPoint);
788 collisions_cnt++;
789 }
790
791 // Each level - distance to the readout plane
792 const Double_t level1 = 0.15; // InductionGapThickness;
793 const Double_t level2 = 0.2 + 0.15; // InductionGapThickness+SecondTransferGapThickness;
794 const Double_t level3 =
795 0.2 + 0.15 + 0.25; // InductionGapThickness+SecondTransferGapThickness+FirstTransferGapThickness;
796
797 // Mean electron shift along x-axis (under the influence of the Lorentz force)
798 Double_t xmean, sigma; // the dependence fitted by polynomial: f(x) = (p0 + p1*x + p2*x^2 + p3*x^3)
799 // Double_t p0_xmean, p1_xmean, p2_xmean, p3_xmean;
800 // Double_t p0_sigma, p1_sigma, p2_sigma, p3_sigma, p4_sigma, p5_sigma;
801 Double_t p_xmean4[4] = {+4.45E-04, +0.0518626, +0.0676802, +0.00195376}; // ArC4H10 (80:20, 0.4T - T=1.5 GeV)
802 Double_t p_xmean5[4] = {+5.52E-04, +0.0659937, +0.0809207,
803 +0.00718051}; // ArC4H10 (80:20, 0.5T - Kr T = 2.36 Gev, B=0.5036 T)
804 Double_t p_xmean6[4] = {+0.000602183, +0.0809795, +0.0918447,
805 +0.0157542}; // ArC4H10 (80:20, 0.6T - T=2.5 GeV)
806 Double_t p_xmean7[4] = {+0.000909404, +0.0896112, +0.125486,
807 +0.00723235}; // ArC4H10 (80:20, 0.7T - T=3.5 GeV)
808 Double_t p_xmean8[4] = {+0.000962017, +0.104726, +0.137705, +0.016483}; // ArC4H10 (80:20, 0.8T - T=4.0 GeV)
809 Double_t p_xmean9[4] = {+0.00103555, +0.117826, +0.160241, +0.0178768}; // ArC4H10 (80:20, 0.9T - T=4.5 GeV)
810 // Sigma electron smearing
811 // Double_t sigma; //depends on the distance from current z-position to the readout plane
812 Double_t p_sigma4[6] = {+0.000181581, +0.120391, -0.425589, +0.887732, -0.885041, +0.332821};
813 Double_t p_sigma5[6] = {+1.87E-04, +0.12034, -0.412872, +0.830354, -0.800251, +0.292069};
814 Double_t p_sigma6[6] = {+0.000137712, +0.125844, -0.468967, +1.01775, -1.04621, +0.403336};
815 Double_t p_sigma7[6] = {+0.000179813, +0.122855, -0.442513, +0.950076, -0.978707, +0.380569};
816 Double_t p_sigma8[6] = {+1.86E-04, +0.122116, -0.435668, +0.934093, -0.963362, +0.376114};
817 Double_t p_sigma9[6] = {+0.000190609, +0.121689, -0.424431, +0.881712, -0.871142, +0.323333};
818
819 Double_t *p_xmean = p_xmean8, *p_sigma = p_sigma8;
820
821 if (ifield == 4) {
822 p_xmean = p_xmean4;
823 p_sigma = p_sigma4;
824 } else if (ifield == 5) {
825 p_xmean = p_xmean5;
826 p_sigma = p_sigma5;
827 } else if (ifield == 6) {
828 p_xmean = p_xmean6;
829 p_sigma = p_sigma6;
830 } else if (ifield == 7) {
831 p_xmean = p_xmean7;
832 p_sigma = p_sigma7;
833 } else if (ifield == 9) {
834 p_xmean = p_xmean9;
835 p_sigma = p_sigma9;
836 }
837
838 // stepEL = fEnergyLossToSignal*point->GetEnergyLoss()/(collisions_cnt + 1);
839
840 for (size_t icoll = 0; icoll < collision_points.size(); ++icoll) {
841 Double_t xcoll = collision_points[icoll][0]; // x
842 Double_t ycoll = collision_points[icoll][1]; // y
843 Double_t zcoll = collision_points[icoll][2]; // z
844
845 Double_t zdist; // current z-distance to the readout
846
847 // Find z-distance to the readout depending on the electron drift direction
848 // Consider that electron drift direction equals to forward Z axis drift
849
850 // AZ if(stationNr % 2 != 0) zdist = moduleThickness - (zcoll - zin);
851 // AZ else zdist = zcoll - zin;
852 zdist = TMath::Abs(zcoll - zin); // AZ
853 if (stationNr % 2 != 0 && point->GetPz() > 0)
854 zdist = moduleThickness - zdist; // AZ
855 else if (stationNr % 2 == 0 && point->GetPz() < 0)
856 zdist = moduleThickness - zdist; // AZ
857 // std::cout << "zdist , zvec" << zdist << vec << std::endl;
858 Double_t zdist2 = zdist * zdist, zdist3 = zdist2 * zdist;
859
860 xmean = p_xmean[0] + p_xmean[1] * zdist + p_xmean[2] * zdist2 + p_xmean[3] * zdist3;
861 sigma = p_sigma[0] + p_sigma[1] * zdist + p_sigma[2] * zdist2 + p_sigma[3] * zdist3
862 + p_sigma[4] * zdist3 * zdist + p_sigma[5] * zdist3 * zdist2;
863
864 // Number of electrons in the current collision
865 Int_t n_electrons_cluster = gRandom->Landau(1.027, 0.11);
866 if (n_electrons_cluster < 1)
867 n_electrons_cluster = 1; // min
868 if (n_electrons_cluster > 6)
869 n_electrons_cluster = 6; // max
870
871 for (Int_t ielectron = 0; ielectron < n_electrons_cluster; ++ielectron) {
872
873 // Electron gain in each GEM cascade
874 // Polya distribution is better, but Exponential is good too in our case
875 Double_t gain_gem1 =
876 gRandom->Exp(15); //...->Exp(V), where V is the mean value of the exponential distribution
877 Double_t gain_gem2 = gRandom->Exp(15);
878 Double_t gain_gem3 = gRandom->Exp(15);
879
880 // if(gain_gem1 < 1.0) gain_gem1 = 1.0;
881 // if(gain_gem2 < 1.0) gain_gem2 = 1.0;
882 // if(gain_gem3 < 1.0) gain_gem3 = 1.0;
883
884 int total_gain = 0;
885
886 if (zdist < level1) {
887 total_gain = 1.0;
888 } else if (zdist >= level1 && zdist < level2) {
889 if (gain_gem3 < 1.0)
890 gain_gem3 = 1.0;
891 total_gain = gain_gem3;
892 } else if (zdist >= level2 && zdist < level3) {
893 if (gain_gem3 < 1.0)
894 gain_gem3 = 1.0;
895 if (gain_gem2 < 1.0)
896 gain_gem2 = 1.0;
897 total_gain = gain_gem3 * gain_gem2;
898 } else if (zdist >= level3) {
899 if (gain_gem3 < 1.0)
900 gain_gem3 = 1.0;
901 if (gain_gem2 < 1.0)
902 gain_gem2 = 1.0;
903 if (gain_gem1 < 1.0)
904 gain_gem1 = 1.0;
905 total_gain = gain_gem3 * gain_gem2 * gain_gem1;
906 }
907
908 // Projection of the current electron on the readout (x,y-coordinates)
909 double x_readout, y_readout;
910 // stepEl_current = stepEL/(total_gain + 1);
911 for (int igain = 0; igain < total_gain; ++igain) {
912
913 // x-shift of the electon depends on the electron drift direction
914 // Consider that electron drift direction equals to forward Z axis drift
915
916 if (stationNr % 2 != 0)
917 x_readout = gRandom->Gaus(xcoll - xmean, sigma);
918 else
919 x_readout = gRandom->Gaus(xcoll + xmean, sigma);
920
921 y_readout = gRandom->Gaus(ycoll, sigma);
922
923 Int_t iIChan = sensor->GetFrontChannel(x_readout, y_readout, zcoll, dPitch);
924 // AZ!!! Exclude 50% of channels for large stations
925 // if (sensor->GetNChannelsFront() > 900 && sensor->GetSectorNr() == 1 && iIChan > 500) iIChan = -1;
926 // else if (sensor->GetNChannelsFront() > 900 && sensor->GetSectorNr() == 3 && iIChan < 500) iIChan
927 // = -1;
928 /*
929 if ( iIChan != -1 ) {
930 fStripSignalF[iIChan] += 1.0;//stepEl_current;
931 fFChannelPointsMap[iIChan].insert(iPoint);
932 }
933 */
934 // Double_t dPrel = dPitch / dx;
935 // if ( iIChan != -1 && dPitch/dx < 0.2) { // strip width 0.16 mm
936 // AZ-200322 if ( iIChan != -1 && TMath::Abs(dPitch/dx) < 0.1) { // strip width 0.16 mm
937 if (iIChan != -1) { // AZ-200322
938 // AZ-200322 fStripSignalF[iIChan] += 1.0;//stepEl_current;
939 fStripSignalF[iIChan] += 0.25; // stepEl_current;
940 fFChannelPointsMap[iIChan].insert(iPoint);
941 // AZ-200322 continue; // electron does not reach backside strips
942 }
943
944 // AZ iIChan = sensor->GetBackChannel (xin,yin,zin);
945 iIChan = sensor->GetBackChannel(x_readout, y_readout, zcoll, dPitch);
946
947 // AZ!!! Exclude 50% of channels for large stations
948 // if (sensor->GetNChannelsFront() > 900 && sensor->GetSectorNr() == 1 && iIChan > 600) iIChan = -1;
949 // else if (sensor->GetNChannelsFront() > 900 && sensor->GetSectorNr() == 3 && iIChan < 600) iIChan
950 // = -1;
951
952 // if ( iIChan != -1 && dPitch/dy < 0.85) { // strip width 0.68 mm
953 // AZ-200322 if ( iIChan != -1 && TMath::Abs(dPitch/dy) < 0.425) { // strip width 0.68 mm
954 if (iIChan != -1) { // AZ-200322
955 // AZ-200322 if (gRandom->Rndm() < 0.68) continue; // kill electron to equalize responses from 2
956 // sides AZ-200322 fStripSignalB[iIChan] += 1.0;//stepEl_current;
957 fStripSignalB[iIChan] += 0.25; // stepEl_current;
958 fBChannelPointsMap[iIChan].insert(iPoint);
959 }
960 }
961 }
962 }
963 }
964}
965
966// AZ ----- Private method ProduceHitResponseAZ ----------------------------
967void CbmStsDigitize::ProduceHitResponseAZ(CbmStsSensor* sensor)
968{
969
970 set<Int_t> pSet;
971 if (fPointMap.find(sensor) == fPointMap.end()) {
972 cerr << "-E- " << fName << "::ProduceHitResponseAZ:: sensor"
973 << " not found in digi scheme!" << endl;
974 return;
975 }
976 pSet = fPointMap[sensor];
977
978 Int_t iPoint = -1;
979 CbmStsPoint* point = NULL;
980
981 set<Int_t>::iterator it1;
982
983 const Int_t nclus = 25; // number of primary clusters / cm
984 const Double_t mfpath = 1.0 / nclus; // mean free path
985 const Double_t sigma = 0.05; // charge spread
986 const Int_t maxClus = 102;
987 // Double_t clusSize[maxClus] = {};
988 const Double_t clusSize[maxClus] = {
989 0.63403, 0.19005, 0.07334, 0.03335, 0.017906, 0.010391, 0.0066257, 0.0046042,
990 0.0034062, 0.0025641, 0.0020853, 0.001613, 0.0012888, 0.0010765, 0.00094932, 0.00087466,
991 0.00079435, 0.00072317, 0.00062637, 0.00053044, 0.00055171, 0.00053304, 0.00056429, 0.00058426,
992 0.00061465, 0.00057688, 0.00058166, 0.00050613, 0.00049571, 0.00046055, 0.00039327, 0.00035898,
993 0.00034162, 0.00030255, 0.00028475, 0.00027998, 0.00023961, 0.0002179, 0.0002153, 0.00019012,
994 0.00018188, 0.00015844, 0.000158, 0.00014368, 0.00013413, 0.00012805, 0.00011373, 0.00011112,
995 0.00011633, 0.000097666, 0.000099837, 0.000082474, 0.000095062, 0.000082474, 0.000091155, 0.000077699,
996 0.000074661, 0.000077699, 0.000074661, 0.000077699, 0.000062072, 0.000059468, 0.000054693, 0.000056429,
997 0.000049484, 0.000043407, 0.000054693, 0.000048616, 0.000044275, 0.000041237, 0.000041671, 0.000043407,
998 0.000041671, 0.000034292, 0.000033858, 0.000039501, 0.000027781, 0.000031687, 0.000029083, 0.000033424,
999 0.000030819, 0.000031253, 0.000032555, 0.000027781, 0.000024308, 0.000022138, 0.000026044, 0.000029083,
1000 0.000019099, 0.000023006, 0.000022572, 0.000020401, 0.000027781, 0.000018231, 0.000024308, 0.000026478,
1001 0.00002344, 0.000018665, 0.000018231, 0.000019099, 0.000016061, 0.0014806};
1002 static TH1D* hClSize = NULL;
1003 if (hClSize == NULL) {
1004 hClSize = new TH1D("hClSize", "", maxClus, 0, maxClus);
1005 for (Int_t i = 0; i < maxClus; ++i)
1006 hClSize->Fill(i + 0.5, clusSize[i]);
1007 }
1008
1009 for (it1 = pSet.begin(); it1 != pSet.end(); it1++) {
1010 iPoint = (*it1);
1011 point = (CbmStsPoint*)fPoints->At(iPoint);
1012
1013 TVector3 posIn, posOut, posCur, vec;
1014 point->PositionIn(posIn);
1015 point->PositionOut(posOut);
1016 vec = posOut;
1017 vec -= posIn;
1018 Double_t leng0 = vec.Mag(), leng = 0.0, dPitch;
1019 vec.SetMag(1.0);
1020 vector<TVector3> clusPos;
1021 vector<Double_t> vClusSize;
1022
1023 // Go thru active material
1024 while (1) {
1025 Double_t step = gRandom->Exp(mfpath);
1026 posCur += step * vec;
1027 leng = posCur.Mag();
1028 if (leng <= leng0) {
1029 clusPos.push_back(posCur);
1030 vClusSize.push_back(Int_t(hClSize->GetRandom() + 1));
1031 if (leng == leng0)
1032 break;
1033 } else
1034 break;
1035 }
1036
1037 // Rescale cluster size sample to account for GEANT energy loss
1038 Int_t nofSteps = clusPos.size();
1039 Double_t nEl = 0, nElGeant = fEnergyLossToSignal * point->GetEnergyLoss();
1040 for (Int_t i = 0; i < nofSteps; ++i)
1041 nEl += vClusSize[i];
1042 Double_t scale = nElGeant / nEl;
1043 for (Int_t i = 0; i < nofSteps; ++i)
1044 vClusSize[i] *= scale;
1045 // cout << " steps: " << nofSteps << " " << nElGeant << " " << leng0 << endl;
1046 // for (Int_t i = 0; i < nofSteps; ++i)
1047 // cout << i << " " << clusPos[i].Mag() << " " << vClusSize[i] << endl;
1048
1049 // Double_t stepEL = fEnergyLossToSignal * point->GetEnergyLoss() / nofSteps;
1050
1051 for (Int_t istep = 0; istep < nofSteps; istep++) {
1052 clusPos[istep] += posIn;
1053 Double_t stepEL = vClusSize[istep];
1054
1055 Int_t iIChan = sensor->GetFrontChannel(clusPos[istep].X(), clusPos[istep].Y(), clusPos[istep].Z(), dPitch);
1056
1057 // AZ!!! Exclude 50% of channels for large stations
1058 // if (sensor->GetNChannelsFront() > 900 && sensor->GetSectorNr() == 1 && iIChan > 500) iIChan = -1;
1059 // else if (sensor->GetNChannelsFront() > 900 && sensor->GetSectorNr() == 3 && iIChan < 500) iIChan = -1;
1060
1061 if (iIChan != -1) {
1062 // fStripSignalF[iIChan] += stepEL;
1063 // fFChannelPointsMap[iIChan].insert(iPoint);
1064 for (Int_t jj = -4; jj < 5; ++jj) {
1065 Int_t ich = iIChan + jj;
1066 if (ich < 0 || ich > sensor->GetNChannelsFront())
1067 continue;
1068 Double_t cx1 = TMath::Erf((-dPitch + jj * sensor->GetDx()) / sigma);
1069 Double_t cx2 = TMath::Erf((-dPitch + (jj + 1) * sensor->GetDx()) / sigma);
1070 Double_t frac = 0.5 * TMath::Abs(cx2 - cx1);
1071 if (frac < 1.e-7)
1072 continue;
1073 Double_t dQ = stepEL * frac;
1074 fStripSignalF[ich] += dQ;
1075 fFChannelPointsMap[ich].insert(iPoint);
1076 }
1077 }
1078
1079 iIChan = sensor->GetBackChannel(clusPos[istep].X(), clusPos[istep].Y(), clusPos[istep].Z(), dPitch);
1080
1081 // AZ!!! Exclude 50% of channels for large stations
1082 // if (sensor->GetNChannelsFront() > 900 && sensor->GetSectorNr() == 1 && iIChan > 600) iIChan = -1;
1083 // else if (sensor->GetNChannelsFront() > 900 && sensor->GetSectorNr() == 3 && iIChan < 600) iIChan = -1;
1084
1085 if (iIChan != -1) {
1086 // fStripSignalB[iIChan] += stepEL;
1087 // fBChannelPointsMap[iIChan].insert(iPoint);
1088 for (Int_t jj = -4; jj < 5; ++jj) {
1089 Int_t ich = iIChan + jj;
1090 if (ich < 0 || ich > sensor->GetNChannelsBack())
1091 continue;
1092 Double_t cx1 = TMath::Erf((-dPitch + jj * sensor->GetDx()) / sigma);
1093 Double_t cx2 = TMath::Erf((-dPitch + (jj + 1) * sensor->GetDx()) / sigma);
1094 Double_t frac = 0.5 * TMath::Abs(cx2 - cx1);
1095 if (frac < 1.e-7)
1096 continue;
1097 Double_t dQ = stepEL * frac;
1098 fStripSignalB[ich] += dQ;
1099 fBChannelPointsMap[ich].insert(iPoint);
1100 }
1101 }
1102 }
1103 }
1104}
1105
1106// ----- Private method SetParContainers -------------------------------
1107void CbmStsDigitize::SetParContainers()
1108{
1109
1110 // Get run and runtime database
1111 FairRunAna* run = FairRunAna::Instance();
1112 if (!run)
1113 Fatal("SetParContainers", "No analysis run");
1114
1115 FairRuntimeDb* db = run->GetRuntimeDb();
1116 if (!db)
1117 Fatal("SetParContainers", "No runtime database");
1118
1119 // Get STS geometry parameter container
1120 fGeoPar = (CbmGeoStsPar*)db->getContainer("CbmGeoStsPar");
1121
1122 // Get STS digitisation parameter container
1123 fDigiPar = (CbmStsDigiPar*)db->getContainer("CbmStsDigiPar");
1124}
1125// -------------------------------------------------------------------------
1126
1127// ----- Private method Init -------------------------------------------
1128InitStatus CbmStsDigitize::Init()
1129{
1130
1131 // Get input array
1132 FairRootManager* ioman = FairRootManager::Instance();
1133 if (!ioman)
1134 Fatal("Init", "No FairRootManager");
1135 fPointsSi = (TClonesArray*)ioman->GetObject("SiliconPoint"); // AZ-280322
1136 fPoints = (TClonesArray*)ioman->GetObject("StsPoint");
1137 fMCTracks = (TClonesArray*)ioman->GetObject("MCTrack"); // AZ
1138
1139 // Register output array StsDigi
1140 fDigis = new TClonesArray("CbmStsDigi", 1000);
1141 ioman->Register("StsDigi", "Digital response in STS", fDigis, kTRUE);
1142
1143 // Register output array StsDigiMatches
1144 fDigiMatches = new TClonesArray("CbmStsDigiMatch", 1000);
1145 ioman->Register("StsDigiMatch", "Digi Match in STS", fDigiMatches, kTRUE);
1146
1147 // fGen = new TRandom3();
1148 // time_t curtime;
1149 // time(&curtime);
1150 // fGen->SetSeed(curtime);
1151
1152 fStripSignalF = new Double_t[2000];
1153 fStripSignalB = new Double_t[2000];
1154
1155 // AZ fEnergyLossToSignal = 280000000.;
1156
1157 fFNofSteps = (Int_t)TMath::Power(2, (Double_t)fFNofBits);
1158 fBNofSteps = (Int_t)TMath::Power(2, (Double_t)fBNofBits);
1159
1160 // Build digitisation scheme
1161 if (fDigiScheme->Init(fGeoPar, fDigiPar)) {
1162 MakeSets();
1163
1164 if (fVerbose == 1 || fVerbose == 2)
1165 fDigiScheme->Print(kFALSE);
1166 else if (fVerbose > 2)
1167 fDigiScheme->Print(kTRUE);
1168 cout << "-I- " << fName << "::Init: "
1169 << "STS digitisation scheme succesfully initialised" << endl;
1170 if (fDigiScheme->IsNewGeometry())
1171 cout << "-I- Using new geometry" << endl;
1172 cout << " Stations: " << fDigiScheme->GetNStations() << ", Sectors: " << fDigiScheme->GetNSectors()
1173 << ", Channels: " << fDigiScheme->GetNChannels() << endl;
1174
1175 // AZ
1176 // fastDig = new BmnGemFastDigitize();
1177 // Add(fastDig);
1178 fastDig = (BmnGemFastDigitize*)FairRun::Instance()->GetTask("GemFastDigitize");
1179
1180 return kSUCCESS;
1181 }
1182
1183 return kERROR;
1184}
1185// -------------------------------------------------------------------------
1186
1187// ----- Private method ReInit -----------------------------------------
1188InitStatus CbmStsDigitize::ReInit()
1189{
1190
1191 // Clear digitisation scheme
1192 fDigiScheme->Clear();
1193
1194 // Build new digitisation scheme
1195 if (fDigiScheme->Init(fGeoPar, fDigiPar)) {
1196 MakeSets();
1197 return kSUCCESS;
1198 }
1199
1200 return kERROR;
1201}
1202// -------------------------------------------------------------------------
1203
1204// ----- Private method MakeSets ---------------------------------------
1205void CbmStsDigitize::MakeSets()
1206{
1207
1208 fPointMap.clear();
1209 Int_t nStations = fDigiScheme->GetNStations();
1210
1211 Double_t fSectorWidth = 0.;
1212
1213 for (Int_t iStation = 0; iStation < nStations; iStation++) {
1214 CbmStsStation* station = fDigiScheme->GetStation(iStation);
1215 Int_t nSectors = station->GetNSectors();
1216 for (Int_t iSector = 0; iSector < nSectors; iSector++) {
1217 CbmStsSector* sector = station->GetSector(iSector);
1218 Int_t nSensors = sector->GetNSensors();
1219 for (Int_t iSensor = 0; iSensor < nSensors; iSensor++) {
1220 CbmStsSensor* sensor = sector->GetSensor(iSensor);
1221 set<Int_t> a;
1222 fPointMap[sensor] = a;
1223 fSectorWidth = 10. * sensor->GetLx();
1224
1225 Int_t nofChips =
1226 (Int_t)(TMath::Ceil(fSectorWidth / 7.5)); // fwidth in mm, 7.5mm = 125(channels)*60mum(pitch)
1227 Int_t lastChip = (Int_t)(TMath::Ceil(10. * fSectorWidth));
1228 lastChip = lastChip % 75;
1229 lastChip = (Int_t)(lastChip / .6);
1230 // cout << nofChips << " chips on " << iStation+1 << " " << iSector+1 << endl;
1231 TString addInfo = "";
1232 if (nofChips != 8) {
1233 addInfo = Form(", only %d strips", lastChip);
1234 // cout << fSectorWidth << " -> " << addInfo.Data() << endl;
1235 }
1236
1237 for (Int_t iChip = 0; iChip < nofChips; iChip++) {
1238 occupancy[iStation][iSector][iChip] = 3.;
1239 // cout << "OCCUPANCY [" << iStation+1 << "][" << iSector+1 << "][" << iChip << "] "<<
1240 // occupancy [iStation][iSector][iChip] << "%" << endl;
1241 }
1242 }
1243 }
1244 }
1245 fFChannelPointsMap.clear();
1246 fBChannelPointsMap.clear();
1247 for (Int_t ichan = 2000; ichan > 0;) {
1248 set<Int_t> a;
1249 fFChannelPointsMap[--ichan] = a;
1250 set<Int_t> b;
1251 fBChannelPointsMap[ichan] = b;
1252 }
1253}
1254// -------------------------------------------------------------------------
1255void CbmStsDigitize::MakeSets1()
1256{ // with occupancy file - default not used
1257
1258 fPointMap.clear();
1259 Int_t nStations = fDigiScheme->GetNStations();
1260
1261 TH1F* fhFNofDigisPChip[10][1000][20];
1262 TH1F* fhBNofDigisPChip[10][1000][20];
1263 TString qaFileName;
1264 qaFileName = "occup.sts.reco.root";
1265 cout << "Occupancy read from file: \"" << qaFileName.Data() << "\"" << endl;
1266 TFile* occuF = TFile::Open(qaFileName.Data());
1267
1268 TString directoryName = "STSFindHitsQA";
1269
1270 Double_t fSectorWidth = 0.;
1271
1272 for (Int_t iStation = 0; iStation < nStations; iStation++) {
1273 CbmStsStation* station = fDigiScheme->GetStation(iStation);
1274 Int_t nSectors = station->GetNSectors();
1275
1276 for (Int_t iSector = 0; iSector < nSectors; iSector++) {
1277 CbmStsSector* sector = station->GetSector(iSector);
1278 Int_t nSensors = sector->GetNSensors();
1279 // Int_t nChannels = sector->GetNChannelsFront();
1280
1281 for (Int_t iSensor = 0; iSensor < nSensors; iSensor++) {
1282 CbmStsSensor* sensor = sector->GetSensor(iSensor);
1283 set<Int_t> a;
1284 fPointMap[sensor] = a;
1285 fSectorWidth = 10. * sensor->GetLx();
1286
1287 Int_t nofChips =
1288 (Int_t)(TMath::Ceil(fSectorWidth / 7.5)); // fwidth in mm, 7.5mm = 125(channels)*60mum(pitch)
1289 Int_t lastChip = (Int_t)(TMath::Ceil(10. * fSectorWidth));
1290 lastChip = lastChip % 75;
1291 lastChip = (Int_t)(lastChip / .6);
1292 // cout << nofChips << " chips on " << iStation+1 << " " << iSector+1 << endl;
1293 TString addInfo = "";
1294 if (nofChips != 8) {
1295 addInfo = Form(", only %d strips", lastChip);
1296 // cout << fSectorWidth << " -> " << addInfo.Data() << endl;
1297 }
1298
1299 for (Int_t iChip = 0; iChip < nofChips; iChip++) {
1300 fhFNofDigisPChip[iStation][iSector][iChip] =
1301 (TH1F*)occuF->Get(Form("%s/Station%d/hNofFiredDigisFSt%dSect%dChip%d", directoryName.Data(),
1302 iStation + 1, iStation + 1, iSector + 1, iChip + 1));
1303 fhBNofDigisPChip[iStation][iSector][iChip] =
1304 (TH1F*)occuF->Get(Form("%s/Station%d/hNofFiredDigisBSt%dSect%dChip%d", directoryName.Data(),
1305 iStation + 1, iStation + 1, iSector + 1, iChip + 1));
1306 occupancy[iStation][iSector][iChip] =
1307 100. * fhFNofDigisPChip[iStation][iSector][iChip]->GetMean() / 125.;
1308 occupancy[iStation][iSector][iChip] =
1309 100. * fhBNofDigisPChip[iStation][iSector][iChip]->GetMean() / 125.;
1310 // if ( !occuF ) {
1311 // occupancy [iStation][iSector][iChip] = 3.;
1312 // }
1313 // cout << "OCCUPANCY [" << iStation+1 << "][" << iSector+1 << "][" << iChip << "] "<<
1314 // occupancy [iStation][iSector][iChip] << "%" << endl;
1315 }
1316 }
1317 }
1318 }
1319 fFChannelPointsMap.clear();
1320 fBChannelPointsMap.clear();
1321 for (Int_t ichan = 2000; ichan > 0;) {
1322 set<Int_t> a;
1323 fFChannelPointsMap[--ichan] = a;
1324 set<Int_t> b;
1325 fBChannelPointsMap[ichan] = b;
1326 }
1327}
1328// -------------------------------------------------------------------------
1329
1330// ----- Private method Reset ------------------------------------------
1331void CbmStsDigitize::Reset()
1332{
1333 // fNPoints = fNDigisFront = fNDigisBack = fTime = 0.;
1334 fNDigis = fNMulti = 0;
1335 fFChannelPointsMap.clear();
1336 fBChannelPointsMap.clear();
1337 // if ( fDigis ) fDigis->Clear();
1338 // if ( fDigiMatches ) fDigiMatches->Clear();
1339 if (fDigis)
1340 fDigis->Delete();
1341 if (fDigiMatches)
1342 fDigiMatches->Delete();
1343}
1344// -------------------------------------------------------------------------
1345
1346// ----- Virtual method Finish -----------------------------------------
1347void CbmStsDigitize::Finish()
1348{
1349 cout << endl;
1350 cout << "============================================================" << endl;
1351 cout << "===== " << fName << ": Run summary " << endl;
1352 cout << "===== " << endl;
1353 cout << "===== Events processed : " << setw(8) << fNEvents << endl;
1354 cout.setf(ios_base::fixed, ios_base::floatfield);
1355 cout << "===== Real time per event : " << setw(8) << setprecision(4) << fTime / fNEvents << " s" << endl;
1356 cout << "===== StsPoints per event : " << setw(8) << setprecision(2) << fNPoints / fNEvents << endl;
1357 cout << "===== StsDigis per event : " << setw(8) << setprecision(2)
1358 << (fNDigisFront + fNDigisBack) / fNEvents << endl;
1359 cout << "===== Front digis per point : " << setw(8) << setprecision(2) << fNDigisFront / (fNPoints) << endl;
1360 cout << "===== Back digis per point : " << setw(8) << setprecision(2) << fNDigisBack / (fNPoints) << endl;
1361 cout << "============================================================" << endl;
1362}
1363// -------------------------------------------------------------------------
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
BmnSiliconPoint::PositionOut
void PositionOut(TVector3 &pos)
Definition BmnSiliconPoint.h:44
BmnSiliconPoint::PositionIn
void PositionIn(TVector3 &pos)
Definition BmnSiliconPoint.h:43
CbmMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition CbmMCTrack.h:56
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::IsNewGeometry
Bool_t IsNewGeometry() const
Definition CbmStsDigiScheme.h:65
CbmStsDigiScheme::GetSensorByName
CbmStsSensor * GetSensorByName(TString sensorName)
Definition CbmStsDigiScheme.cxx:233
CbmStsDigiScheme::GetStation
CbmStsStation * GetStation(Int_t iStation)
Definition CbmStsDigiScheme.cxx:184
CbmStsDigi::GetNofAdcBits
static Int_t GetNofAdcBits()
Definition CbmStsDigi.h:68
CbmStsDigitize::Exec
virtual void Exec(Option_t *opt)
Definition CbmStsDigitize.cxx:236
CbmStsDigitize::~CbmStsDigitize
virtual ~CbmStsDigitize()
Definition CbmStsDigitize.cxx:215
CbmStsDigitize::ProduceHitResponseSi
void ProduceHitResponseSi(CbmStsSensor *sensor)
Definition CbmStsDigitize.cxx:551
CbmStsDigitize::ProduceHitResponseAZ
void ProduceHitResponseAZ(CbmStsSensor *sensor)
Definition CbmStsDigitize.cxx:967
CbmStsDigitize::ProduceHitResponse
void ProduceHitResponse(CbmStsSensor *sensor)
Definition CbmStsDigitize.cxx:673
CbmStsDigitize::Finish
virtual void Finish()
Definition CbmStsDigitize.cxx:1347
CbmStsPoint::PositionOut
void PositionOut(TVector3 &pos)
Definition CbmStsPoint.h:79
CbmStsPoint::GetXIn
Double_t GetXIn() const
Definition CbmStsPoint.h:69
CbmStsPoint::PositionIn
void PositionIn(TVector3 &pos)
Definition CbmStsPoint.h:78
CbmStsPoint::GetZOut
Double_t GetZOut() const
Definition CbmStsPoint.h:74
CbmStsPoint::GetYOut
Double_t GetYOut() const
Definition CbmStsPoint.h:73
CbmStsPoint::GetZIn
Double_t GetZIn() const
Definition CbmStsPoint.h:71
CbmStsPoint::GetYIn
Double_t GetYIn() const
Definition CbmStsPoint.h:70
CbmStsPoint::GetXOut
Double_t GetXOut() const
Definition CbmStsPoint.h:72
CbmStsSector::GetSectorNr
Int_t GetSectorNr() const
Definition CbmStsSector.h:76
CbmStsSector::GetStationNr
Int_t GetStationNr() const
Definition CbmStsSector.h:74
CbmStsSector::GetDx
Double_t GetDx() const
Definition CbmStsSector.h:95
CbmStsSector::GetNChannelsFront
Int_t GetNChannelsFront() const
Definition CbmStsSector.h:79
CbmStsSector::GetNChannelsBack
Int_t GetNChannelsBack() const
Definition CbmStsSector.h:80
CbmStsSector::GetSensor
CbmStsSensor * GetSensor(Int_t iSensor)
Definition CbmStsSector.h:85
CbmStsSector::GetNSensors
Int_t GetNSensors() const
Definition CbmStsSector.h:82
CbmStsSensor::GetNChannelsFront
Int_t GetNChannelsFront() const
Definition CbmStsSensor.h:109
CbmStsSensor::GetNChannelsBack
Int_t GetNChannelsBack() const
Definition CbmStsSensor.h:110
CbmStsSensor::GetFrontChannel
Int_t GetFrontChannel(Double_t x, Double_t y, Double_t z, Double_t &dPitch)
Definition CbmStsSensor.cxx:614
CbmStsSensor::GetDx
Double_t GetDx() const
Definition CbmStsSensor.h:101
CbmStsSensor::GetBackChannel
Int_t GetBackChannel(Double_t x, Double_t y, Double_t z, Double_t &dPitch)
Definition CbmStsSensor.cxx:650
CbmStsSensor::GetStationNr
Int_t GetStationNr() const
Definition CbmStsSensor.h:85
CbmStsSensor::GetLx
Double_t GetLx() const
Definition CbmStsSensor.h:98
CbmStsStation::GetNSectors
Int_t GetNSectors() const
Definition CbmStsStation.h:62
CbmStsStation::GetSector
CbmStsSector * GetSector(Int_t iSector)
Definition CbmStsStation.h:65
run
-clang-format