165 TF1* digiGausDist =
new TF1(
"digiGausDist",
"gaus",-5.,5.);
168 if ( fVerbose > 2 ) {
169 cout << endl <<
"-I- " << fName <<
": executing event" << endl;
170 cout <<
"----------------------------------------------" << endl;
173 Double_t nPoints = 0.;
174 Double_t nOutside = 0.;
175 Double_t nDigisF = 0.;
176 Double_t nDigisB = 0.;
180 cerr <<
"-W- " << fName <<
"::Exec: No input array (STSPoint) "
182 cout <<
"- " << fName << endl;
186 map<CbmStsSensor*, set<Int_t> >::iterator mapIt;
187 for (mapIt=fPointMap.begin(); mapIt!=fPointMap.end(); mapIt++)
188 ((*mapIt).second).clear();
190 for (Int_t iPoint=0; iPoint<fPoints->GetEntriesFast(); iPoint++) {
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();
201 if ( fPointMap.find(sensor) == fPointMap.end() ) {
202 cerr <<
"-E- " << fName <<
"::Exec:: sensor " << curNode->GetName()
203 <<
" not found in digi scheme!" << endl;
206 fPointMap[sensor].insert(iPoint);
209 for (Int_t iStation=fDigiScheme->
GetNStations(); iStation > 0 ; ) {
211 for (Int_t iSector=station->
GetNSectors(); iSector > 0 ; ) {
214 map<Int_t, set<Int_t> >::iterator mapCh;
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();
225 for (Int_t iChannel=nChannels ; iChannel > 0 ; ) {
230 fStripSignalF[--iChannel] = TMath::Abs(gRandom->Gaus(0.,fFNoiseWidth));
231 fStripSignalB[ iChannel] = TMath::Abs(gRandom->Gaus(0.,fBNoiseWidth));
234 for (Int_t iSensor=sector->
GetNSensors(); iSensor > 0 ; ) {
244 for ( Int_t ifstr = 0 ; ifstr < nChannels ; ifstr++ ) {
245 if ( fStripSignalF[ifstr] < fFThreshold )
continue;
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,
254 set<Int_t>::iterator it1;
255 set<Int_t> chPnt = fFChannelPointsMap[ifstr];
258 if ( chPnt.size() == 0 ) {
264 for (it1=chPnt.begin(); it1!=chPnt.end(); it1++) {
266 if ( it1==chPnt.begin() )
269 match->AddPoint(pnt);
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,
285 set<Int_t>::iterator it1;
286 set<Int_t> chPnt = fBChannelPointsMap[ibstr];
289 if ( chPnt.size() == 0 ) {
295 for (it1=chPnt.begin(); it1!=chPnt.end(); it1++) {
297 if ( it1==chPnt.begin() )
300 match->AddPoint(pnt);
311 cout <<
"+ " << flush;
312 cout << setw(15) << left << fName <<
": " << setprecision(4) << setw(8)
313 << fixed << right << fTimer.RealTime()
314 <<
" s, digis " << nDigisF <<
" / " << nDigisB << endl;
317 fNPoints += Double_t(nPoints);
318 fNOutside += Double_t(nOutside);
319 fNDigisFront += Double_t(nDigisF);
320 fNDigisBack += Double_t(nDigisB);
321 fTime += fTimer.RealTime();
329 if ( fPointMap.find(sensor) == fPointMap.end() ) {
330 cerr <<
"-E- " << fName <<
"::Exec:: sensor"
331 <<
" not found in digi scheme!" << endl;
334 pSet = fPointMap[sensor];
339 set<Int_t>::iterator it1;
341 for (it1=pSet.begin(); it1!=pSet.end(); it1++) {
345 Double_t xin = point->
GetXIn();
346 Double_t yin = point->
GetYIn();
347 Double_t zin = point->
GetZIn();
349 Double_t xvec = point->
GetXOut()-xin;
350 Double_t yvec = point->
GetYOut()-yin;
351 Double_t zvec = point->
GetZOut()-zin;
353 Int_t nofSteps = (Int_t)(TMath::Sqrt(xvec*xvec+yvec*yvec+zvec*zvec)/fStep+1);
355 Double_t stepEL = fEnergyLossToSignal*point->GetEnergyLoss()/(nofSteps+1);
357 xvec = xvec/((Double_t)nofSteps);
358 yvec = yvec/((Double_t)nofSteps);
359 zvec = zvec/((Double_t)nofSteps);
361 for ( Int_t istep = 0 ; istep <= nofSteps ; istep++ ) {
366 if ( iIChan != -1 ) {
367 fStripSignalF[iIChan] += stepEL;
368 fFChannelPointsMap[iIChan].insert(iPoint);
375 if ( iIChan != -1 ) {
376 fStripSignalB[iIChan] += stepEL;
377 fBChannelPointsMap[iIChan].insert(iPoint);
395 Double_t xin = pnt->
GetXIn();
396 Double_t yin = pnt->
GetYIn();
397 Double_t zin = pnt->
GetZIn();
399 gGeoManager->FindNode(xin,yin,zin);
400 TGeoNode* curNode = gGeoManager->GetCurrentNode();
404 Double_t xvec = pnt->
GetXOut()-xin;
405 Double_t yvec = pnt->
GetYOut()-yin;
406 Double_t zvec = pnt->
GetZOut()-zin;
408 Int_t nofSteps = (Int_t)(TMath::Sqrt(xvec*xvec+yvec*yvec+zvec*zvec)/fStep+1);
410 Double_t stepEL = fEnergyLossToSignal*pnt->GetEnergyLoss()/(nofSteps+1);
412 xvec = xvec/((Double_t)nofSteps);
413 yvec = yvec/((Double_t)nofSteps);
414 zvec = zvec/((Double_t)nofSteps);
416 for ( Int_t istep = 0 ; istep <= nofSteps ; istep++ ) {
418 Int_t iIChan = (Int_t)iFChan;
424 if ( iIChan == -1 )
continue;
426 for ( Int_t ifstr = 0 ; ifstr < nofStr ; ifstr++ ) {
427 if ( strips[ifstr] == iIChan ) {
428 signals[ifstr] += stepEL;
433 if ( iIChan == -1 )
continue;
435 strips [nofStr] = iIChan;
436 signals[nofStr] = stepEL;
580 cout <<
"============================================================"
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 <<
"============================================================"