247 #include "StEvent/StTpcHit.h"
249 #include "StTpcHitMaker.h"
251 #include "TDataSetIter.h"
252 #include "StDAQMaker/StDAQReader.h"
256 #include "StEvent/StTpcHitCollection.h"
257 #include "StEvent/StTpcHit.h"
258 #include "RTS/src/DAQ_TPX/tpxFCF_flags.h"
259 #include "StTpcRawData.h"
260 #include "StThreeVectorF.hh"
262 #include "StDaqLib/TPC/trans_table.hh"
263 #include "StRtsTable.h"
264 #include "StDbUtilities/StTpcCoordinateTransform.hh"
265 #include "StTpcDb/StTpcDb.h"
266 #include "StDbUtilities/StCoordinates.hh"
267 #include "StDetectorDbMaker/St_tss_tssparC.h"
268 #include "StDetectorDbMaker/St_tpcSlewingC.h"
269 #ifdef __CORRECT_S_SHAPE__
270 #include "StDetectorDbMaker/St_TpcPadCorrectionC.h"
272 #include "StDetectorDbMaker/St_tpcPadGainT0BC.h"
273 #include "StDetectorDbMaker/St_tpcAnodeHVavgC.h"
274 #include "StDetectorDbMaker/St_tpcMaxHitsC.h"
275 #include "StDetectorDbMaker/StDetectorDbTpcRDOMasks.h"
276 #include "StDetectorDbMaker/St_tpcPadConfigC.h"
277 #include "StDetectorDbMaker/St_tpcStatusC.h"
281 #include "St_tpc_cl.h"
283 #include "St_daq_cld.h"
285 #include "St_daq_sim_cld.h"
287 #include "St_daq_adc_tb.h"
289 #include "St_daq_sim_adc_tb.h"
292 static TNtuple *pulserP = 0;
293 Float_t StTpcHitMaker::fgDp = .1;
294 Float_t StTpcHitMaker::fgDt = .2;
295 Float_t StTpcHitMaker::fgDperp = .1;
296 Bool_t StTpcHitMaker::fgCosmics = kFALSE;
297 static Int_t _debug = 0;
299 #define __NOT_ZERO_SUPPRESSED_DATA__
302 #ifdef __NOT_ZERO_SUPPRESSED_DATA__
303 #include "StDetectorDbMaker/St_tpcPedestalC.h"
305 static const Char_t *Names[StTpcHitMaker::kAll] = {
"undef",
306 "tpc_hits",
"tpx_hits",
"itpc_hits",
307 "TpcPulser",
"TpxPulser",
"iTPCPulser",
308 "TpcRaw",
"TpxRaw",
"iTPCRaw",
309 "TpcAvLaser",
"TpxAvLaser",
"tpc_hitsO"};
311 StTpcHitMaker::StTpcHitMaker(
const char *name) :
StRTSBaseMaker(
"tpc",name), kMode(kUndefined),
312 kReaderType(kUnknown), mQuery(
""), fTpc(0), fAvLaser(0), fSectCounts(0), fThr(0), fSeq(0) {
313 SetAttr(
"minSector",1);
314 SetAttr(
"maxSector",24);
316 SetAttr(
"UseTonkoClusterAnnotation",1);
319 Int_t StTpcHitMaker::Init() {
320 LOG_INFO <<
"StTpcHitMaker::Init as\t" <<
GetName() << endm;
322 for (Int_t k = 1; k < kAll; k++) {
323 if (MkName.CompareTo(Names[k],TString::kIgnoreCase) == 0) {kMode = (EMode) k;
break;}
326 memset(maxHits,0,
sizeof(maxHits));
332 void StTpcHitMaker::InitializeHistograms(Int_t token) {
333 static Int_t oldToken = -1;
334 Int_t newToken = token/10;
335 if (newToken == oldToken)
return;
336 TFile *f = GetTFile();
337 if (! f) {gMessMgr->Error() <<
"with Tpx/Tpc AvLaser you must provide TFile as the 5-th parameter in bfc.C macro" << endm; assert(0);}
340 if (fSectCounts) {fSectCounts->Write();
delete fSectCounts;}
342 for (Int_t s = 1; s <= 24; s++) {
343 if (fAvLaser[s-1]) {fAvLaser[s-1]->Write();
delete fAvLaser[s-1];}
345 delete [] fAvLaser; fAvLaser = 0;
349 const Char_t *NameV[NoDim] = {
"row",
"pad",
"time"};
350 const Double_t xMin[NoDim] = {0.5 , 0.5, -0.5};
351 const Double_t xMax[NoDim] = {0.5+St_tpcPadConfigC::instance()->numberOfRows(20), 182.5, 399.5};
352 Int_t nBins[NoDim] = { St_tpcPadConfigC::instance()->numberOfRows(20), 182, 400};
353 fSectCounts =
new TH1F(Form(
"SectorCounts_%03i",newToken),
"Count no. of sectors",25,-0.5,24.5);
354 #ifdef __USE__THnSparse__
355 fAvLaser =
new THnSparseF *[24];
356 for (Int_t s = 1; s <= 24; s++) {
357 fAvLaser[s-1] =
new THnSparseF(Form(
"AvLaser_%02i_%03i",s,newToken), Form(
"Averaged laser event for sector %02i for token %03i",s,newToken),
358 NoDim, nBins, xMin, xMax);
360 for (Int_t i = 0; i < NoDim; i++) {
361 fAvLaser[s-1]->GetAxis(i)->SetName(NameV[i]);
363 f->Add(fAvLaser[s-1]);
366 fAvLaser =
new TH3F *[24];
367 TH1::SetDefaultSumw2(kTRUE);
368 for (Int_t s = 1; s <= 24; s++) {
369 TString name(Form(
"AvLaser_%02i",s));
370 if (newToken) name += Form(
"_%03i",newToken);
371 fAvLaser[s-1] =
new TH3F(name, Form(
"Averaged laser event for sector %02i for token %03i",s,newToken),
372 nBins[0],xMin[0],xMax[0],
373 nBins[1],xMin[1],xMax[1],
374 nBins[2],xMin[2],xMax[2]);
375 fAvLaser[s-1]->GetXaxis()->SetTitle(NameV[0]);
376 fAvLaser[s-1]->GetYaxis()->SetTitle(NameV[1]);
377 fAvLaser[s-1]->GetZaxis()->SetTitle(NameV[2]);
383 Int_t StTpcHitMaker::InitRun(Int_t runnumber) {
384 SetAttr(
"maxRow",St_tpcPadConfigC::instance()->numberOfRows(20));
385 if (IAttr(
"Cosmics")) SetCosmics();
389 Int_t maxHitsPerSector = St_tpcMaxHitsC::instance()->maxSectorHits();
390 Int_t maxBinZeroHits = St_tpcMaxHitsC::instance()->maxBinZeroHits();
393 Float_t liveFrac = 1;
394 for(Int_t sector=1;sector<=24;sector++) {
395 Int_t liveSecPads = 0;
396 Int_t totalSecPads = 0;
397 if (maxHitsPerSector > 0 || maxBinZeroHits > 0) {
398 for(Int_t row=1;row<=St_tpcPadConfigC::instance()->numberOfRows(sector);row++) {
399 Int_t numPadsAtRow = St_tpcPadConfigC::instance()->padsPerRow(sector,row);
400 totalSecPads += numPadsAtRow;
401 if (StDetectorDbTpcRDOMasks::instance()->isRowOn(sector,row,1) &&
402 St_tpcAnodeHVavgC::instance()->livePadrow(sector,row))
403 liveSecPads += numPadsAtRow;
405 livePads += liveSecPads;
406 totalPads += totalSecPads;
408 if (maxHitsPerSector > 0) {
409 liveFrac = TMath::Max(0.1f,
410 ((Float_t) liveSecPads) / (1e-15f + (Float_t) totalSecPads));
411 maxHits[sector-1] = (Int_t) (liveFrac * maxHitsPerSector);
412 if (Debug()) {LOG_INFO <<
"maxHits in sector " << sector
413 <<
" = " << maxHits[sector-1] << endm;}
415 maxHits[sector-1] = 0;
416 if (Debug()) {LOG_INFO <<
"No maxHits in sector " << sector << endm;}
419 if (maxBinZeroHits > 0) {
420 liveFrac = TMath::Max(0.1f,
421 ((Float_t) livePads) / ((Float_t) totalPads));
422 maxBin0Hits = (Int_t) (liveFrac * maxBinZeroHits);
423 if (Debug()) {LOG_INFO <<
"maxBinZeroHits " << maxBin0Hits << endm;}
426 if (Debug()) {LOG_INFO <<
"No maxBinZeroHits" << endm;}
429 if (kMode == kTpxAvLaser || kMode == kTpcAvLaser) {
432 TFile *f = GetTFile();
433 if (! f) {gMessMgr->Error() <<
"with Tpx/Tpc AvLaser you must provide TFile as the 5-th parameter in bfc.C macro" << endm; assert(0);}
442 if (St_tpcStatusC::instance()->isDead()) {
443 LOG_WARN <<
"TPC status indicates it is unusable for this event. Ignoring hits." << endm;
447 static Int_t minSector = IAttr(
"minSector");
448 static Int_t maxSector = IAttr(
"maxSector");
449 static Int_t minRow = IAttr(
"minRow");
450 static Int_t maxRow = IAttr(
"maxRow");
451 if (kMode == kTpxAvLaser || kMode == kTpcAvLaser) {
453 InitializeHistograms(
Token());
455 InitializeHistograms(0);
458 StMaker* maskMk = GetMakerInheritsFrom(
"StMtdTrackingMaskMaker");
459 unsigned int mask = (maskMk ? maskMk->UAttr(
"TpcSectorsByMtd") : ~0U);
461 if (fSectCounts) fSectCounts->Fill(0);
462 static const Char_t *tpcDataNames[5] = {0,
"tpc/legacy",
"tpx/legacy",
"tpx",
"itpc"};
463 TString cldadc(
"cld");
464 if ( kMode == kTpxRaw || kMode == kTpcRaw || kMode == kiTPCRaw ||
465 kMode == kTpcAvLaser || kMode == kTpxAvLaser) cldadc =
"adc";
466 for (Int_t sector = minSector; sector <= maxSector; sector++) {
467 if (!((1U<<(sector-1)) & mask))
continue;
471 Int_t kMin = kUnknown;
472 for (Int_t k = kStandardiTPC; k > kMin; k--) {
474 mQuery = Form(
"%s/%s[%i]",tpcDataNames[k],cldadc.Data(),sector);
476 mQuery = Form(
"%s[%i]",tpcDataNames[k],sector);
478 if (! daqTpcTable)
continue;
481 kReaderType = (EReaderType) k;
482 if (kReaderType > kLegacyTpx) kMin = kLegacyTpx;
483 while (daqTpcTable) {
484 if (Sector() == sector) {
485 if (Debug()/100 > 0) {
486 daqTpcTable->
Print(0,10);
490 if (kReaderType == kLegacyTpx || kReaderType == kLegacyTpc) fTpc = (
tpc_t*)*
DaqDta()->begin();
491 Int_t row = RowNumber();
492 if (row >= minRow && row <= maxRow) {
494 LOG_INFO <<
"StTpcHitMaker::Make(" << Names[kMode] <<
") => " << tpcDataNames[k] <<
" for sector = " << sector <<
" row = " << row << endm;
500 case kTpx: hitsAdded += UpdateHitCollection(sector);
break;
502 case kTpxPulser:
if (fTpc) DoPulser(sector);
break;
505 if ( fTpc) TpcAvLaser(sector);
506 else TpxAvLaser(sector);
507 fSectCounts->Fill(sector);
512 if ( fTpc) RawTpcData(sector);
513 else RawTpxData(sector);
524 if (maxHits[sector-1] && hitsAdded > maxHits[sector-1]) {
525 LOG_ERROR <<
"Too many hits (" << hitsAdded <<
") in one sector ("
526 << sector <<
"). Skipping event." << endm;
530 if (maxBin0Hits && bin0Hits > maxBin0Hits) {
531 LOG_ERROR <<
"Too many hits (" << bin0Hits
532 <<
") starting at time bin 0. Skipping event." << endm;
535 if (kMode == kTpc || kMode == kTpx || kMode == kTpxO) {
536 StEvent *pEvent =
dynamic_cast<StEvent *
> (GetInputDS(
"StEvent"));
537 if (Debug()) {LOG_INFO <<
"StTpcHitMaker::Make : StEvent has been retrieved " <<pEvent<< endm;}
538 if (! pEvent) {LOG_INFO <<
"StTpcHitMaker::Make : StEvent has not been found " << endm;
return kStWarn;}
540 if (hitCollection && ! IAttr(
"NoTpxAfterBurner")) AfterBurner(hitCollection);
542 if (IAttr(
"CheckThrSeq") && (kMode == kTpcRaw || kMode == kTpxRaw || kMode == kiTPCRaw)) {
548 void StTpcHitMaker::CheckThrSeq() {
549 if (! fThr || !fSeq) {
550 fThr =
new TH2C(
"Thr",
"ADC Thresold value versus sector and row",24,0.5,24.5,72,0.5,72.5); fThr->SetDirectory(0);
551 fSeq =
new TH2C(
"Seq",
"ADC sequnce value versus sector and row",24,0.5,24.5,72,0.5,72.5); fSeq->SetDirectory(0);
552 for (Int_t s = 1; s <= 24; s++)
553 for (Int_t r = 1; r <= 72; r++) {
554 fThr->SetBinContent(s,r,127);
555 fSeq->SetBinContent(s,r,127);
558 for (Int_t s = 1; s <= 24; s++){
560 if (! digitalSector)
continue;
561 Int_t Nrows = digitalSector->numberOfRows();
562 for (Int_t r = 1; r <= Nrows; r++) {
563 Int_t Npads = digitalSector->numberOfPadsInRow(r);
564 for (Int_t p = 1; p <= Npads; p++) {
565 Int_t ntb = digitalSector->numberOfTimeBins(r,p);
567 digitalSector->getTimeAdc(r,p,ADCs,IDTs);
570 Int_t tbF = -1, tbL = -1;
571 for (Int_t tb = 0; tb < __MaxNumberOfTimeBins__; tb++) {
573 if (tbF > -1 && tbL >= tbF) {
574 if (seq > tbL - tbF + 1) seq = tbL - tbF + 1;
579 if (ADCs[tb] < adcMin) {
580 adcMin = TMath::Max(ADCs[tb-2],TMath::Max(ADCs[tb-1],TMath::Max(ADCs[tb],TMath::Max(ADCs[tb+1],ADCs[tb+1]))));
582 if (tbF < 0) tbF = tb;
585 if (adcMin < 127 && seq < 127) {
586 if (seq == 1 && adcMin == 1) {
587 static Int_t ibreak = 0;
590 Char_t th = fThr->GetBinContent(s,r);
591 if (th > adcMin) fThr->SetBinContent(s,r, adcMin);
592 Char_t sq = fSeq->GetBinContent(s,r);
593 if (sq > seq) fSeq->SetBinContent(s,r, seq);
601 #ifdef __USE__THnSparse__
602 if (GetTFile() && fAvLaser) {
603 for (Int_t sector = 1; sector <= 24; sector++) {
604 if (fAvLaser[sector-1]) {
605 THnSparseF *hnew = CompressTHn(fAvLaser[sector-1]);
606 GetTFile()->Remove(fAvLaser[sector-1]);
607 delete fAvLaser[sector-1];
608 fAvLaser[sector-1] = hnew;
609 GetTFile()->Add(fAvLaser[sector-1]);
615 Int_t adcMinFL = 127;
617 Int_t s1 = -1, s2 = -1, r1 = -1, r2 = -1;
618 for (Int_t s = 1; s <= 24; s++){
619 for (Int_t r = 1; r <= 72; r++) {
620 Int_t adcMin = fThr->GetBinContent(s,r);
621 Int_t seq = fSeq->GetBinContent(s,r);
622 if (adcMin == 127 || seq == 127)
continue;
623 if (adcMinFL == 127 && seqFL == 127) {
629 if (adcMinFL == adcMin && seqFL == seq) {
633 LOG_INFO <<
"StTpcHitMaker::Finish CheckThrSeq in sectors [" << s1 <<
"," << s2 <<
"] and rows[" << r1 <<
"," << r2 <<
"] Threshold = " << adcMinFL <<
" and sequence = " << seqFL << endm;
642 LOG_INFO <<
"StTpcHitMaker::Finish CheckThrSeq in sectors [" << s1 <<
"," << s2 <<
"] and rows[" << r1 <<
"," << r2 <<
"] Threshold = " << adcMinFL <<
" and sequence = " << seqFL << endm;
649 Int_t StTpcHitMaker::UpdateHitCollection(Int_t sector) {
651 StEvent *pEvent =
dynamic_cast<StEvent *
> (GetInputDS(
"StEvent"));
652 if (Debug()) {LOG_INFO <<
"StTpcHitMaker::Make : StEvent has been retrieved " <<pEvent<< endm;}
653 if (! pEvent) {LOG_INFO <<
"StTpcHitMaker::Make : StEvent has not been found " << endm;
return 0;}
655 if ( !hitCollection ) {
658 pEvent->setTpcHitCollection(hitCollection);
661 if (NRows <= 0)
return 0;
662 Int_t nhitsBefore = hitCollection->numberOfHits();
663 Int_t sec =
DaqDta()->Sector();
664 Int_t row = RowNumber();
665 if (kReaderType == kLegacyTpc || kReaderType == kLegacyTpx) {
667 for (Int_t l = 0; l < NRows; tpc++) {
668 if ( !tpc->has_clusters )
return 0;
669 for(Int_t padrow=0;padrow<St_tpcPadConfigC::instance()->numberOfRows(sector);padrow++) {
670 tpc_cl *c = &tpc->cl[padrow][0];
671 Int_t ncounts = tpc->cl_counts[padrow];
672 for(Int_t j=0;j<ncounts;j++,c++) {
673 if (! c || ! c->charge)
continue;
675 (c->flags & ~(FCF_ONEPAD | FCF_MERGED | FCF_BIG_CHARGE)))
continue;
676 if (kMode == kTpxO) c->flags |= 256;
677 Int_t row = padrow + 1;
679 Int_t iRdo = StDetectorDbTpcRDOMasks::instance()->rdoForPadrow(sector,row,pad);
680 if ( ! StDetectorDbTpcRDOMasks::instance()->isOn(sector,iRdo))
continue;
681 StTpcHit *tpcHit = CreateTpcHit(*c,sector,row);
682 Int_t iok = hitCollection->addHit(tpcHit);
691 LOG_INFO << Form(
"CLD sec %2d: row %2d: clusters: %3d",sec, row, NRows) << endm;
693 for (Int_t l = 0; l < NRows; l++, cld++) {
695 LOG_INFO << Form(
" pad %f[%d:%d], tb %f[%d:%d], cha %d, fla 0x%X",
706 if (! cld->pad || ! cld->charge)
continue;
707 if (! cld->tb)
continue;
708 if (cld->tb < 0 || cld->tb >= __MaxNumberOfTimeBins__)
continue;
709 if (cld->t1 < 0 || cld->t1 >= __MaxNumberOfTimeBins__)
continue;
710 if (cld->t2 < 0 || cld->t2 >= __MaxNumberOfTimeBins__)
continue;
712 (cld->flags & ~(FCF_ONEPAD | FCF_MERGED | FCF_BIG_CHARGE)))
continue;
713 if (kMode == kTpxO) cld->flags |= 256;
714 Float_t pad = cld->pad;
715 Int_t iRdo = StDetectorDbTpcRDOMasks::instance()->rdoForPadrow(sector,row,pad);
716 if ( ! StDetectorDbTpcRDOMasks::instance()->isOn(sector,iRdo))
continue;
717 StTpcHit *tpcHit = CreateTpcHit(*cld,sector,row);
718 Int_t iok = hitCollection->addHit(tpcHit);
722 Int_t nhits = hitCollection->numberOfHits() - nhitsBefore;
724 LOG_INFO <<
" Total hits in Sector : row " << sector <<
" : " << row <<
" = " << nhits << endm;
729 StTpcHit *StTpcHitMaker::CreateTpcHit(
const tpc_cl &cluster, Int_t sector, Int_t row) {
732 Float_t pad = cluster.p;
733 Float_t time = cluster.t;
734 if (kReaderType == kLegacyTpx) time += 22;
739 transform(padcoord,local,kFALSE);
740 transform(local,global);
747 Double_t gain = (row<=St_tpcPadConfigC::instance()->innerPadRows(sector)) ? St_tss_tssparC::instance()->gain_in() : St_tss_tssparC::instance()->gain_out();
748 Double_t wire_coupling = (row<=St_tpcPadConfigC::instance()->innerPadRows(sector)) ? St_tss_tssparC::instance()->wire_coupling_in() : St_tss_tssparC::instance()->wire_coupling_out();
751 StTpcHit *
hit = StTpcHitFlag(global.position(),hard_coded_errors,hw,q
764 if (hit->minTmbk() == 0) bin0Hits++;
766 LOG_INFO <<
"StTpcHitMaker::CreateTpcHit fromt tpc_cl\t" <<*hit << endm;
771 StTpcHit *StTpcHitMaker::CreateTpcHit(
const daq_cld &cluster, Int_t sector, Int_t row) {
774 Float_t pad = cluster.pad;
775 Float_t time = cluster.tb;
777 Double_t gain = (row<=St_tpcPadConfigC::instance()->innerPadRows(sector)) ? St_tss_tssparC::instance()->gain_in() : St_tss_tssparC::instance()->gain_out();
778 Double_t wire_coupling = (row<=St_tpcPadConfigC::instance()->innerPadRows(sector)) ? St_tss_tssparC::instance()->wire_coupling_in() : St_tss_tssparC::instance()->wire_coupling_out();
779 Double_t q = cluster.charge * ((Double_t)St_tss_tssparC::instance()->ave_ion_pot() * (Double_t)St_tss_tssparC::instance()->scale())/(gain*wire_coupling) ;
785 if (tpcSlewing && tpcSlewing->type() != 1001) tpcSlewing = 0;
788 Double_t freq = gStTpcDb->Electronics()->samplingFrequency();
789 time = freq * tpcSlewing->correctedT(sector,row,q,time/freq);
795 transform(padcoord,local,kFALSE);
796 transform(local,global);
802 UShort_t flag = cluster.flags;
803 if (kMode == kTpxO) flag |= 256;
805 StTpcHit *hit = StTpcHitFlag(global.position(),hard_coded_errors,hw,q
818 if (hit->minTmbk() == 0) bin0Hits++;
820 LOG_INFO <<
"StTpcHitMaker::CreateTpcHit fromt daq_cld\t" <<*hit << endm;
825 void StTpcHitMaker::DoPulser(Int_t sector) {
826 struct Pulser_t {Float_t sector, row, pad, gain, t0, nnoise, noise, npeak;};
827 static const Char_t *names =
"sector:row:pad:gain:t0:nnoise:noise:npeak";
828 static Pulser_t Pulser;
830 TFile *f = GetTFile();
833 pulserP =
new TNtuple(
"pulserP",
"Pulser analysis",names);
835 Int_t r, p, tb, tbmax;
838 if (! fTpc->channels_sector)
return;
839 for(Int_t row = 1; row <= St_tpcPadConfigC::instance()->numberOfRows(sector); row++) {
841 if (! fTpc->cl_counts[r])
continue;
842 for (Int_t pad = 1; pad <= 182; pad++) {
844 Int_t ncounts = fTpc->counts[r][p];
845 if (! ncounts)
continue;
846 static UShort_t adc[512];
847 memset (adc, 0,
sizeof(adc));
850 for (Int_t i = 0; i < ncounts; i++) {
851 tb = fTpc->timebin[r][p][i];
852 adc[tb] = log8to10_table[fTpc->adc[r][p][i]];
853 if (adc[tb] > adcmax) {
858 if (tbmax < 2 || tbmax > 504)
continue;
860 Int_t i1s = TMath::Max( 0, tbmax - 2);
861 Int_t i2s = TMath::Min(511, tbmax + 7);
862 Int_t i1 = TMath::Max(0 ,i1s - 20);
863 Int_t i2 = TMath::Min(511,i2s + 20);
867 for (Int_t i = i1; i <= i2; i++) {
868 if (i >= i1s && i <= i2s)
continue;
872 if (nnoise) noise /= nnoise;
873 for (Int_t i = i1s; i <= i2s; i++) {
875 peak += adc[i] - noise;
876 t0 += i*(adc[i] - noise);
878 if (peak <= 0)
continue;
880 Pulser.sector = sector;
885 Pulser.nnoise = nnoise;
886 Pulser.noise = noise;
887 Pulser.npeak = npeak;
888 pulserP->Fill(&Pulser.sector);
893 void StTpcHitMaker::TpcAvLaser(Int_t sector) {
894 if (! fTpc || !fAvLaser)
return;
897 Double_t sector, row, pad, time;
899 #ifdef __USE__THnSparse__
900 if (fAvLaser[sector-1]->GetNbins() > 500000) {
901 THnSparseF *hnew = CompressTHn(fAvLaser[sector-1]);
902 GetTFile()->Remove(fAvLaser[sector-1]);
903 delete fAvLaser[sector-1];
904 fAvLaser[sector-1] = hnew;
905 GetTFile()->Add(fAvLaser[sector-1]);
909 pixel.sector = sector;
910 for(Int_t r = 0; r < St_tpcPadConfigC::instance()->numberOfRows(sector); r++) {
912 for (Int_t pad = 1; pad <= 182; pad++) {
915 Double_t gain = St_tpcPadGainT0BC::instance()->Gain(pixel.sector,pixel.row,pixel.pad);
916 if (gain <= 0)
continue;
917 Int_t ncounts = fTpc->counts[r][p];
918 if (! ncounts)
continue;
919 for (Int_t i = 0; i < ncounts; i++) {
920 pixel.time = fTpc->timebin[r][p][i];
921 Double_t adc = log8to10_table[fTpc->adc[r][p][i]];
922 #ifdef __USE__THnSparse__
923 fAvLaser[sector-1]->Fill(&pixel.row,adc);
925 fAvLaser[sector-1]->Fill(pixel.row,pixel.pad,pixel.time,adc);
931 LOG_INFO <<
" Total pixels in Sector : " << sector <<
" = " << npixels << endm;
934 void StTpcHitMaker::TpxAvLaser(Int_t sector) {
935 assert(fAvLaser[sector-1]);
936 #ifdef __USE__THnSparse__
937 if (fAvLaser[sector-1]->GetNbins() > 1000000) {
938 THnSparseF *hnew = CompressTHn(fAvLaser[sector-1]);
939 GetTFile()->Remove(fAvLaser[sector-1]);
940 delete fAvLaser[sector-1];
941 fAvLaser[sector-1] = hnew;
942 GetTFile()->Add(fAvLaser[sector-1]);
945 Int_t r=RowNumber() ;
948 Int_t p =
Pad() - 1 ;
949 if (p < 0 || p >= St_tpcPadConfigC::instance()->padsPerRow(sector,r+1))
return;
951 Double_t sector, row, pad, time;
954 pixel.sector = sector;
958 for (;iword !=
DaqDta()->end();++iword) {
960 Int_t tb = daqadc.tb;
962 Double_t adc = daqadc.adc;
963 #ifdef __NOT_ZERO_SUPPRESSED_DATA__
965 if (tb >= 368 && tb <= 383)
967 adc -= St_tpcPedestalC::instance()->Pedestal(sector,r+1,p+1);
970 adc -= St_tpcPedestalC::instance()->Pedestal(sector,r+1,p+1);
975 #ifdef __USE__THnSparse__
976 fAvLaser[sector-1]->Fill(&pixel.row,adc);
978 fAvLaser[sector-1]->Fill(pixel.row,pixel.pad,pixel.time,adc);
983 void StTpcHitMaker::DumpPixels2Ntuple(Int_t sector, Int_t row, Int_t pad) {
985 Float_t event, sector, row, pad, tb, adc, idt;
987 static const Char_t *BName =
"event:sector:row:pad:tb:adc:idt";
988 static TNtuple *adcP = 0;
989 if (! adcP && GetTFile() ) {
991 adcP =
new TNtuple(
"adcP",
"Pulser ADC",BName);
995 P.event = GetEventNumber();
999 for (Int_t i = 0; i < __MaxNumberOfTimeBins__; i++) {
1000 if (! ADCs[i])
continue;
1004 adcP->Fill(&P.event);
1008 void StTpcHitMaker::PrintSpecial(Int_t sector) {
1016 UInt_t tot_pix = 0 ;
1017 UInt_t cl_count = 0 ;
1020 for(r=0;r<St_tpcPadConfigC::instance()->numberOfRows(sector);r++) {
1021 for(p=0;p<182;p++) {
1022 for(t=0;t<fTpc->counts[r][p];t++) {
1023 val = fTpc->adc[r][p][t] ;
1024 Int_t vali = log8to10_table[val];
1028 Int_t timebin = fTpc->timebin[r][p][t] ;
1029 printf(
"%d %d %d %d %d\n",sector,r+1,p+1,timebin,vali) ;
1034 if(fTpc->has_clusters) {
1035 cl_count += fTpc->cl_counts[r] ;
1038 if(fTpc->has_clusters) {
1039 for(i=0;i<fTpc->cl_counts[r];i++) {
1040 tpc_cl *c = &fTpc->cl[r][i] ;
1042 printf(
"%d %d %f %f %d %d %d %d %d %d\n",
1043 sector,r+1,c->p,c->t,c->charge,c->flags,c->p1,c->p2,c->t1,c->t2) ;
1048 LOG_INFO << Form(
"TPC: Sector %d: occupancy %3d %%, charge %d, pixels %u, clusters %d",sector,
1049 (
int)(100.0 *((
double)fTpc->channels_sector/(
double)fTpc->max_channels_sector)),
1050 adc,tot_pix,cl_count) << endm;
1055 TDataSet *
event = GetData(
"Event");
1064 if (! digitalSector) {
1066 data->setSector(sector,digitalSector);
1068 return digitalSector;
1071 Int_t StTpcHitMaker::RawTpxData(Int_t sector) {
1072 static Int_t TonkoAnn = IAttr(
"UseTonkoClusterAnnotation");
1074 #ifndef __TFG__VERSION__
1075 UShort_t IDTs2[512];
1079 memset(ADCs, 0,
sizeof(ADCs));
1080 memset(IDTs, 0,
sizeof(IDTs));
1084 Int_t Total_data = 0;
1085 Int_t r=RowNumber() ;
1087 Int_t p =
Pad() - 1 ;
1089 if (p < 0 || p >= St_tpcPadConfigC::instance()->padsPerRow(sector,r+1))
return 0;
1091 if (iword ==
DaqDta()->end())
return 0;
1092 Int_t some_data = 0;
1095 Total_data += some_data;
1097 if (! digitalSector) digitalSector = GetDigitalSector(sector);
1098 Int_t ntbold = digitalSector->numberOfTimeBins(r_old+1,p_old+1);
1101 LOG_INFO <<
"digitalSector " << sector <<
" row " << r+1 <<
" pad " << p + 1
1102 <<
" already has " << ntbold <<
" time bins at row/pad " << r_old+1 <<
"/" << p_old+1 << endm;
1104 digitalSector->getTimeAdc(r_old+1,p_old+1,ADCs2,IDTs2);
1107 for (Int_t i = 0; i < __MaxNumberOfTimeBins__; i++) {
1108 if (! ADCs2[i])
continue;
1109 if ((IDTs[i] || IDTs2[i]) && ADCs[i] < ADCs2[i]) IDTs[i] = IDTs2[i];
1110 ADCs[i] += ADCs2[i];
1114 digitalSector->putTimeAdc(r_old+1,p_old+1,ADCs,IDTs);
1115 if (IAttr(
"TpxDumpPxls2Nt")) {
1116 DumpPixels2Ntuple(sector,r_old+1,p_old+1);
1118 memset(ADCs, 0,
sizeof(ADCs));
1119 memset(IDTs, 0,
sizeof(IDTs));
1121 if (r_old != r || p_old != p) {
1125 for (;iword !=
DaqDta()->end();++iword) {
1127 Int_t tb = daqadc.tb;
1128 Int_t adc = daqadc.adc;
1129 #ifdef __NOT_ZERO_SUPPRESSED_DATA__
1131 adc -= St_tpcPedestalC::instance()->Pedestal(sector,r+1,p+1);
1132 if (adc <= 0)
continue;
1141 }
while (some_data);
1142 if (Debug() && Total_data > 0) {
1143 LOG_INFO <<
"StTpcHitMaker::RawTpxData \tsector = " << sector <<
"\trow = " << r+1 <<
"\tpad = " << p+1 << endm;
1144 digitalSector->PrintTimeAdc(r+1,p+1);
1149 Int_t StTpcHitMaker::RawTpcData(Int_t sector) {
1150 static Int_t TonkoAnn = IAttr(
"UseTonkoClusterAnnotation");
1151 if (! fTpc)
return 0;
1152 memset(ADCs, 0,
sizeof(ADCs));
1153 memset(IDTs, 0,
sizeof(IDTs));
1155 assert( digitalSector );
1156 Int_t Total_data = 0;
1157 for (Int_t row = 1; row <= digitalSector->numberOfRows(); row++) {
1159 for (Int_t pad = 1; pad <= digitalSector->numberOfPadsAtRow(row); pad++) {
1161 memset(ADCs, 0,
sizeof(ADCs));
1162 memset(IDTs, 0,
sizeof(IDTs));
1163 Int_t ncounts = fTpc->counts[r][p];
1164 if (! ncounts)
continue;
1165 for (Int_t i = 0; i < ncounts; i++) {
1166 Int_t tb = fTpc->timebin[r][p][i];
1167 ADCs[tb] = log8to10_table[fTpc->adc[r][p][i]];
1174 Int_t ntbold = digitalSector->numberOfTimeBins(row,pad);
1177 LOG_INFO <<
"digitalSector " << sector
1178 <<
" already has " << ntbold <<
" at row/pad " << row <<
"/" << pad << endm;
1181 digitalSector->putTimeAdc(row,pad,ADCs,IDTs);
1182 if (IAttr(
"TpxDumpPxls2Nt")) {
1183 DumpPixels2Ntuple(sector,row,pad);
1188 LOG_INFO <<
"Read " << Total_data <<
" pixels from Sector " << sector << endm;
1193 #ifdef __MAKE_NTUPLE__
1194 static TNtuple *tup = 0;
1195 struct TpcHitPair_t {
1197 qK, padK, tbkK, padKmn, padKmx, tbkKmn, tbkKmx, IdTK, QAK,
1198 qL, padL, tbkL, padLmn, padLmx, tbkLmn, tbkLmx, IdTL, QAL,
1201 static const Char_t *vTpcHitMRPair =
"sec:row:"
1202 "qK:padK:tbkK:padKmn:padKmx:tbkKmn:tbkKmx:IdTK:QAK:"
1203 "qL:padL:tbkL:padLmn:padLmx:tbkLmn:tbkLmx:IdTL:QAL:"
1208 return (200*lhs->timeBucket() + lhs->pad()) < (200*rhs->timeBucket() + rhs->pad());
1212 static Float_t padDiff = 2.5, timeBucketDiff = 5.0;
1216 if (! TpcHitCollection)
return;
1218 #ifdef __MAKE_NTUPLE__
1220 if (StChain::GetChain()->GetTFile()) {
1221 StChain::GetChain()->GetTFile()->cd();
1222 tup =
new TNtuple(
"HitT",
"Cluster Pair Info",vTpcHitMRPair);
1227 Char_t mergedHits[65536];
1228 UInt_t numberOfSectors = TpcHitCollection->numberOfSectors();
1229 UInt_t TotNoHits = 0;
1230 UInt_t RejNoHits = 0;
1231 for (UInt_t sec = 1; sec <= numberOfSectors; sec++) {
1233 if (sectorCollection) {
1234 UInt_t numberOfPadrows = sectorCollection->numberOfPadrows();
1235 for (UInt_t row = 1; row <= numberOfPadrows; row++) {
1237 if (rowCollection) {
1238 UInt_t NoHits = rowCollection->hits().size();
1239 TotNoHits += NoHits;
1240 if (NoHits < 2)
continue;
1241 sort(rowCollection->hits().begin(),
1242 rowCollection->hits().end(),
1246 Char_t mergePass = 0;
1247 memset(mergedHits,0,65536*
sizeof(Char_t));
1248 Int_t newlyMerged = 1;
1249 while (newlyMerged > 0) {
1250 if (mergePass == 255) {
1251 LOG_WARN <<
"Too many merging passes! Afterburner aborted for sector="
1252 << sec <<
" row=" << row << endm;
1256 for (UInt_t k = 0; k < NoHits; k++) {
1257 if (mergedHits[k] != mergePass)
continue;
1258 StTpcHit* kHit = rowCollection->hits().at(k);
1259 if (_debug) {cout <<
"k " << k; kHit->Print();}
1260 if (kHit->flag())
continue;
1261 if (kHit->adc() <= 0.0)
continue;
1262 #ifdef __MAKE_NTUPLE__
1265 pairC.qK = kHit->charge();
1266 pairC.padK = kHit->pad();
1267 pairC.tbkK = kHit->timeBucket();
1268 pairC.padKmn = kHit->minPad();
1269 pairC.padKmx = kHit->maxPad();
1270 pairC.tbkKmn = kHit->minTmbk();
1271 pairC.tbkKmx = kHit->maxTmbk();
1272 pairC.IdTK = kHit->idTruth();
1273 pairC.QAK = kHit->qaTruth();
1275 for (UInt_t l = k+1; l < NoHits; l++) {
1276 StTpcHit* lHit = rowCollection->hits().at(l);
1277 if (_debug) {cout <<
"l " << l; lHit->Print();}
1278 if (lHit->flag())
continue;
1279 if (lHit->adc() <= 0.0)
continue;
1282 bool notNear = (TMath::Abs(kHit->pad() - lHit->pad()) > padDiff ||
1283 TMath::Abs(kHit->timeBucket() - lHit->timeBucket()) > timeBucketDiff);
1284 #ifndef __MAKE_NTUPLE__
1285 if (notNear)
continue;
1288 Int_t padOverlap = TMath::Min(kHit->maxPad(),lHit->maxPad())
1289 - TMath::Max(kHit->minPad(),lHit->minPad());
1290 if (padOverlap < 0)
continue;
1291 Int_t tmbkOverlap = TMath::Min(kHit->maxTmbk(),lHit->maxTmbk())
1292 - TMath::Max(kHit->minTmbk(),lHit->minTmbk());
1293 if (tmbkOverlap < 0)
continue;
1294 #ifdef __MAKE_NTUPLE__
1296 pairC.qL = lHit->charge();
1297 pairC.padL = lHit->pad();
1298 pairC.tbkL = lHit->timeBucket();
1299 pairC.padLmn = lHit->minPad();
1300 pairC.padLmx = lHit->maxPad();
1301 pairC.tbkLmn = lHit->minTmbk();
1302 pairC.tbkLmx = lHit->maxTmbk();
1303 pairC.IdTL = lHit->idTruth();
1304 pairC.QAL = lHit->qaTruth();
1305 pairC.padOv = padOverlap;
1306 pairC.tbkOv = tmbkOverlap;
1307 tup->Fill(&pairC.sec);
1309 if (notNear)
continue;
1312 UShort_t flag = lHit->flag() | FCF_CHOPPED;
1314 UShort_t flag = lHit->flag() | 0x080;
1316 lHit->setFlag(flag);
1319 cout <<
"mk" << k; kHit->Print();
1320 cout <<
"ml" << l; lHit->Print();
1322 Float_t q = lHit->charge() + kHit->charge();
1323 Float_t adc = lHit->adc() + kHit->adc();
1324 Float_t pad = (lHit->adc()*lHit->pad() + kHit->adc()*kHit->pad() )/adc;
1325 Float_t timeBucket = (lHit->adc()*lHit->timeBucket() + kHit->adc()*kHit->timeBucket())/adc;
1326 Short_t minPad = TMath::Min(lHit->minPad(), kHit->minPad());
1327 Short_t maxPad = TMath::Max(lHit->maxPad(), kHit->maxPad());
1328 Short_t minTmbk = TMath::Min(lHit->minTmbk(), kHit->minTmbk());
1329 Short_t maxTmbk = TMath::Max(lHit->maxTmbk(), kHit->maxTmbk());
1330 Int_t IdT = lHit->idTruth();
1331 Double_t QA = lHit->adc()*lHit->qaTruth();
1332 if (IdT == kHit->idTruth()) {
1333 QA += kHit->adc()*kHit->qaTruth();
1335 if (lHit->adc()*lHit->qaTruth() < kHit->adc()*kHit->qaTruth()) {
1336 QA = kHit->adc()*kHit->qaTruth();
1337 IdT = kHit->idTruth();
1341 kHit->setIdTruth(IdT, TMath::Nint(QA));
1343 kHit->setExtends(pad,timeBucket,minPad,maxPad,minTmbk,maxTmbk);
1346 transform(padcoord,local,kFALSE);
1347 transform(local,global);
1348 kHit->setPosition(global.position());
1350 cout <<
"m " << k; kHit->Print();
1354 mergedHits[k] = mergePass + 1;
1360 #ifdef __CORRECT_S_SHAPE__
1362 for (UInt_t k = 0; k < NoHits; k++) {
1363 StTpcHit* kHit = TpcHitCollection->sector(sec-1)->padrow(row-1)->hits().at(k);
1364 if (kHit->flag())
continue;
1365 Double_t pad = kHit->pad();
1366 Double_t timeBucket = kHit->timeBucket();
1368 if (row > St_tpcPadConfigC::instance()->innerPadRows(sector)) io = 2;
1369 Int_t np = kHit->padsInHit();
1370 pad += St_TpcPadCorrectionC::instance()->GetCorrection(pad,io,np,0);
1372 transform(padcoord,local,kFALSE);
1373 transform(local,global);
1374 kHit->setPosition(global.position());
1381 LOG_INFO <<
"StTpcHitMaker::AfterBurner from " << TotNoHits <<
" Total no. of hits " << RejNoHits <<
" were rejected" << endm;
1387 UInt_t hw,
float q, UChar_t c,
1388 Int_t idTruth, UShort_t quality,
1390 UShort_t mnpad, UShort_t mxpad, UShort_t mntmbk,
1391 UShort_t mxtmbk, Float_t cl_x, Float_t cl_t, UShort_t adc,
1394 StTpcHit* hit =
new StTpcHit(p,e,hw,q,c,idTruth,quality,
id,mnpad,mxpad,mntmbk,mxtmbk,cl_x,cl_t,adc);
1396 if ( mntmbk<0 || mxtmbk<0 || mntmbk>500 || mxtmbk>500
1397 || mnpad <0 || mxpad <0 || mnpad >500 || mxpad >500
1398 || mxpad-mnpad > 100
1399 || (Float_t) mntmbk>cl_t || (Float_t) mxtmbk<cl_t
1400 || (Float_t) mnpad >cl_x || (Float_t) mxpad <cl_x
1401 ) flag |= FCF_SANITY;
1403 if (fgCosmics && TMath::Abs(p.z()) > 150) flag |= FCF_SANITY;
1406 assert(! ( TMath::IsNaN(hit->position().x()) || TMath::IsNaN(hit->position().y()) || TMath::IsNaN(hit->position().x())));
1409 #ifdef __USE__THnSparse__
1411 THnSparseF *StTpcHitMaker::CompressTHn(THnSparseF *hist, Double_t compress) {
1412 if (! hist)
return 0;
1413 Int_t nd = hist->GetNdimensions();
1414 Int_t *nbins =
new Int_t[nd];
1415 for (Int_t i = 0; i < nd; i++) nbins[i] = hist->GetAxis(i)->GetNbins();
1416 THnSparseF *hnew =
new THnSparseF(hist->GetName(),hist->GetTitle(),nd, nbins, 0, 0);
1417 hnew->CalculateErrors(kTRUE);
1419 for (Int_t i = 0; i < nd; i++) {
1420 TAxis *ax = hist->GetAxis(i);
1421 if (ax->IsVariableBinSize()) hnew->GetAxis(i)->Set(ax->GetNbins(), ax->GetXbins()->GetArray());
1422 else hnew->GetAxis(i)->Set(ax->GetNbins(), ax->GetXmin(), ax->GetXmax());
1424 Int_t *bins =
new Int_t[nd];
1425 Double_t *x =
new Double_t[nd];
1426 Long64_t N = hist->GetNbins(); cout << hist->GetName() <<
" has " << N <<
" bins before compression." << endl;
1427 #ifndef __NOT_ZERO_SUPPRESSED_DATA__
1429 for (Long64_t i = 0; i < N; ++i) {
1430 Double_t cont = hist->GetBinContent(i, bins);
1431 if (cont > max) max = cont;
1433 for (Long64_t i = 0; i < N; ++i) {
1434 Double_t cont = hist->GetBinContent(i, bins);
1435 if (cont < max/compress)
continue;
1437 for (Int_t d = 0; d < nd; ++d) {x[d] = hist->GetAxis(d)->GetBinCenter(bins[d]);}
1441 for (Long64_t i = 0; i < N; ++i) {
1442 Double_t cont = hist->GetBinContent(i, bins);
1443 if (cont < 1.0)
continue;
1445 for (Int_t d = 0; d < nd; ++d) {x[d] = hist->GetAxis(d)->GetBinCenter(bins[d]);}
1451 cout << hnew->GetName() <<
" has " << hnew->GetNbins() <<
" bins after compression." << endl;
1456 Int_t StTpcHitMaker::RowNumber(){
1457 Int_t sector =
DaqDta()->Sector();
1458 Int_t row =
DaqDta()->Row();
1459 if (row == 0 && (kReaderType == kLegacyTpc || kReaderType == kLegacyTpx)) row = 1;
1460 if (kReaderType != kStandardiTPC) {
1461 if (St_tpcPadConfigC::instance()->iTpc(sector) && row > 13) {
StRtsTable * GetNextDaqElement(const char *elementPath)
Query the STAR production chain for the DAQ data.
Class StRTSBaseMaker - is an abstract StMaker to define the interface to access the DAQ data from the...
virtual void AddData(TDataSet *data, const char *dir=".data")
User methods.
virtual Int_t Write(const char *name=0, Int_t option=0, Int_t bufsize=0)
virtual Long_t GetNRows() const
Returns the number of the used rows for the wrapped table.
StRtsTable * DaqDta()
Return the current DAQ data block. This member function is provided for convenience.
virtual const char * GetName() const
special overload
static UInt_t Token()
current token
virtual Char_t * Print(Char_t *buf, Int_t n) const
Create IDL table defintion (to be used for XDF I/O)