StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StTpcHitMaker.cxx
1 /***************************************************************************
2  *
3  * $Id: StTpcHitMaker.cxx,v 1.82 2021/05/10 21:13:19 fisyak Exp $
4  *
5  * Author: Valeri Fine, BNL Feb 2007
6  ***************************************************************************
7  *
8  * Description: Fill the StEvent from the DAQ clusters
9  *
10  * Input: DAQReader
11  * Output: StTpcHit collection added to StEvent
12  *
13  ***************************************************************************
14  *
15  * $Log: StTpcHitMaker.cxx,v $
16  * Revision 1.82 2021/05/10 21:13:19 fisyak
17  * Clean up
18  *
19  * Revision 1.81 2020/03/20 02:25:35 genevb
20  * Only look for StEvent when needed, e.g. not for raw modes such as in embedding
21  *
22  * Revision 1.80 2020/02/24 23:15:06 genevb
23  * Faster Afterburner(), Remove compiler ambiguity on double vs. float in TMath::Min,Max()
24  *
25  * Revision 1.79 2020/02/20 18:36:54 genevb
26  * Reduce time-expenseive calls to GetInputDS(), some coverity cleanup
27  *
28  * Revision 1.78 2019/05/11 02:20:34 genevb
29  * Add TPC-dead status handling
30  *
31  * Revision 1.77 2019/03/22 18:08:46 fisyak
32  * recover skip of legacy if exists Tpx or iTpc. Checked by Irakli
33  *
34  * Revision 1.76 2019/03/01 15:50:20 fisyak
35  * Fix bug. 3385
36  *
37  * Revision 1.75 2019/02/14 17:40:39 fisyak
38  * Fix selection for row number (Thanks Iraklii for checking)
39  *
40  * Revision 1.74 2018/10/17 20:45:27 fisyak
41  * Restore update for Run XVIII dE/dx calibration removed by Gene on 08/07/2018
42  *
43  * Revision 1.72 2018/06/22 18:35:19 perev
44  * Merging with TPC group code
45  *
46  * Revision 1.55 2018/02/18 23:35:33 perev
47  * Remove iTPC update
48  *
49  * Revision 1.53 2015/04/09 19:54:03 genevb
50  * Introduce sector masking (only MTD-based so far)
51  *
52  * Revision 1.52 2015/03/02 21:10:15 genevb
53  * pad and timebucket sanity units were off by x64, mistakenly using units of StTpcHit data members (RT 3507)
54  *
55  * Revision 1.51 2014/06/26 21:31:41 fisyak
56  * New Tpc Alignment, v632
57  *
58  * Revision 1.49 2013/04/07 21:58:36 fisyak
59  * Move selection of sector range from InitRun to Init, add cluster averaging based on TH3
60  *
61  * Revision 1.48 2013/01/29 23:28:16 fisyak
62  * comment out occupancy print outs
63  *
64  * Revision 1.47 2013/01/28 20:26:50 fisyak
65  * Simplify loop over clusters
66  *
67  * Revision 1.46 2012/12/06 14:33:47 fisyak
68  * Keep only clusters with flag == 0 or FCF_ONEPAD | FCF_MERGED | FCF_BIG_CHARGE. Tonko's instruction
69  *
70  * Revision 1.45 2012/11/20 22:55:56 fisyak
71  * Set the same cuts for online cluster maker as for offline one
72  *
73  * Revision 1.44 2012/10/24 13:36:06 fisyak
74  * Increase no. of pad rows
75  *
76  * Revision 1.43 2012/09/13 21:00:04 fisyak
77  * Corrections for iTpx, clean up
78  *
79  * Revision 1.42 2012/05/07 15:51:01 fisyak
80  * Remove hard coded TPC numbers
81  *
82  * Revision 1.41 2011/06/09 20:52:08 genevb
83  * Set sanity flag
84  *
85  * Revision 1.40 2011/04/07 23:33:12 genevb
86  * Restore to version before previous commit
87  *
88  * Revision 1.38 2011/03/31 19:31:12 fisyak
89  * Add adc to Tpc hit
90  *
91  * Revision 1.37 2011/03/08 18:20:44 genevb
92  * Limit on number of hits starting at time bin 0
93  *
94  * Revision 1.36 2011/01/21 18:35:25 fisyak
95  * change flag type from UChar_t to UShort_t
96  *
97  * Revision 1.35 2011/01/20 18:26:30 genevb
98  * Add FCF_flags include, and exclude any flagged hit from AfterBurner()
99  *
100  * Revision 1.34 2010/11/05 16:25:19 genevb
101  * No longer include hits found on dead padrows
102  *
103  * Revision 1.33 2010/11/04 19:39:12 genevb
104  * Maintain backward reproducibility
105  *
106  * Revision 1.32 2010/11/04 18:30:47 genevb
107  * Typo correction
108  *
109  * Revision 1.31 2010/11/04 18:29:58 genevb
110  * Max hits scaling does not need to use PadGainT0 table
111  *
112  * Revision 1.30 2010/10/04 19:06:56 fisyak
113  * Use FCF flag definition
114  *
115  * Revision 1.29 2010/09/08 15:44:41 genevb
116  * Slightly better arrangement for limiting excessive TPC events
117  *
118  * Revision 1.28 2010/09/01 21:14:33 fisyak
119  * Add codes for S-shape correction (disactivated)
120  *
121  * Revision 1.27 2010/08/31 15:19:36 genevb
122  * Lower bound on reduced hit maxima
123  *
124  * Revision 1.26 2010/08/31 14:16:30 genevb
125  * Correct mistake from prev commit of location of TPC cluster check
126  *
127  * Revision 1.25 2010/08/30 18:02:01 genevb
128  * Introduce hit maxima for tracking
129  *
130  * Revision 1.24 2010/08/02 23:06:15 fisyak
131  * Fix format
132  *
133  * Revision 1.23 2010/03/25 15:05:54 fisyak
134  * Add AfterBurner
135  *
136  * Revision 1.22 2010/02/19 23:36:08 fisyak
137  * Add hit Id
138  *
139  * Revision 1.21 2010/01/12 22:54:36 fisyak
140  * Propagate flags from online clustering into StEvent
141  *
142  * Revision 1.20 2009/11/18 14:29:02 fisyak
143  * Restore slewing correction
144  *
145  * Revision 1.19 2009/11/10 21:05:08 fisyak
146  * Add attributes for sector and pad row selections
147  *
148  * Revision 1.18 2009/09/11 22:11:58 genevb
149  * Introduce TPC slewing corrections
150  *
151  * Revision 1.17 2009/03/18 14:21:06 fisyak
152  * Move sector check under condition that there is some TPC data
153  *
154  * Revision 1.16 2009/03/17 19:19:21 fisyak
155  * Account new Valery's interface for adc values
156  *
157  * Revision 1.15 2009/03/16 13:41:45 fisyak
158  * Switch to new scheme (avoid legacy) for TPX cluster reading
159  *
160  * Revision 1.14 2009/03/11 18:38:20 fisyak
161  * Add 22 time bins to account subtracted by Tonko, clean up
162  *
163  * Revision 1.13 2009/02/20 22:06:15 fisyak
164  * Restore access to TPX
165  *
166  * Revision 1.12 2008/12/29 23:58:06 fine
167  * Optimize the DAQ data access
168  *
169  * Revision 1.11 2008/12/29 21:14:41 fine
170  * Sort out the tps/tpc data handling
171  *
172  * Revision 1.10 2008/12/29 18:23:48 fine
173  * avoid using the dummy data
174  *
175  * Revision 1.9 2008/12/18 20:20:25 fine
176  * access two different detectors tpx/tpc
177  *
178  * Revision 1.8 2008/12/17 23:27:04 fine
179  * Clean up
180  *
181  * Revision 1.7 2008/12/17 23:26:00 fine
182  * Adjust the sector number
183  *
184  * Revision 1.6 2008/12/17 02:04:28 fine
185  * fix the sector number to make the new interface happy
186  *
187  * Revision 1.5 2008/12/16 20:43:25 fine
188  * add the DAQ_READER compliant access to the tpx sector
189  *
190  * Revision 1.4 2008/12/15 21:04:01 fine
191  * For for the NEW_DAQ_READER
192  *
193  * Revision 1.3 2008/07/31 20:45:26 fisyak
194  * Add TpcMixer
195  *
196  * Revision 1.2 2008/06/23 20:13:53 fisyak
197  * Add real data pixel annotation
198  *
199  * Revision 1.1.1.1 2008/05/27 14:22:41 fisyak
200  * Maker to access TPC DAQ information via EVP_READER
201  *
202  * Revision 1.3 2008/05/27 14:18:18 fisyak
203  * Freeze before moving to STAR repository
204  *
205  * Revision 1.2 2008/04/28 14:37:15 fisyak
206  * Rearrage TpcHitMaker to make it run for parallel taks, add the first version of online clustering
207  *
208  * Revision 1.1.1.1 2008/04/03 20:16:41 fisyak
209  * Initial version
210  *
211  * Revision 1.9 2008/01/29 02:44:38 fine
212  * INFO
213  *
214  * Revision 1.8 2008/01/29 02:42:31 fine
215  * remove 16th sector constarin. EVP_READER can read all of them alone now.
216  *
217  * Revision 1.7 2008/01/28 23:48:39 fine
218  * use the new base class
219  *
220  * Revision 1.6 2008/01/10 01:12:49 fine
221  * makr to use the full TPC + TPX
222  *
223  * Revision 1.5 2008/01/09 15:16:48 fine
224  * Correct the sector number
225  *
226  * Revision 1.4 2008/01/09 00:43:29 fine
227  * Working version. It can be used as the protopty for anither maker that calles RTS_READER to fill the 16-th TPX sector
228  *
229  * Revision 1.3 2008/01/07 19:04:07 fine
230  * Add the interface to access the DAQ clusters
231  *
232  * Revision 1.2 2008/01/07 17:37:39 fine
233  * check for tpcHitCollection and new StTpcHit object
234  *
235  * Revision 1.1 2008/01/04 17:52:32 fine
236  * New maker to populate the StEvent from the tpc structure filled by the new EVP_READER package
237  *
238  *
239  * StTpcHitMaker - class to fille the StEvewnt from DAQ reader
240  *
241  **************************************************************************/
242 //#define __MAKE_NTUPLE__
243 //#define __CORRECT_S_SHAPE__
244 //#define __TOKENIZED__
245 //#define __USE__THnSparse__
246 #include <assert.h>
247 #include "StEvent/StTpcHit.h"
248 #include <algorithm>
249 #include "StTpcHitMaker.h"
250 
251 #include "TDataSetIter.h"
252 #include "StDAQMaker/StDAQReader.h"
253 #include "TError.h"
254 #include "string.h"
255 #include "StEvent.h"
256 #include "StEvent/StTpcHitCollection.h"
257 #include "StEvent/StTpcHit.h"
258 #include "RTS/src/DAQ_TPX/tpxFCF_flags.h" // for FCF flag definition
259 #include "StTpcRawData.h"
260 #include "StThreeVectorF.hh"
261 
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"
271 #endif /* __CORRECT_S_SHAPE__ */
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"
278 #include "TFile.h"
279 #include "TNtuple.h"
280 #include "TH2.h"
281 #include "St_tpc_cl.h"
282 TableClassImpl(St_tpc_cl,tcl_cl);
283 #include "St_daq_cld.h"
284 TableClassImpl(St_daq_cld,tcl_cl);
285 #include "St_daq_sim_cld.h"
286 TableClassImpl(St_daq_sim_cld,tcl_cl);
287 #include "St_daq_adc_tb.h"
288 TableClassImpl(St_daq_adc_tb,daq_adc_tb);
289 #include "St_daq_sim_adc_tb.h"
290 TableClassImpl(St_daq_sim_adc_tb,daq_sim_adc_tb);
291 ClassImp(StTpcHitMaker);
292 static TNtuple *pulserP = 0;
293 Float_t StTpcHitMaker::fgDp = .1; // hardcoded errors
294 Float_t StTpcHitMaker::fgDt = .2;
295 Float_t StTpcHitMaker::fgDperp = .1;
296 Bool_t StTpcHitMaker::fgCosmics = kFALSE;
297 static Int_t _debug = 0;
298 #ifdef __TOKENIZED__
299 #define __NOT_ZERO_SUPPRESSED_DATA__
300 #endif
301 //#define __NOT_ZERO_SUPPRESSED_DATA__
302 #ifdef __NOT_ZERO_SUPPRESSED_DATA__
303 #include "StDetectorDbMaker/St_tpcPedestalC.h"
304 #endif
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"};
310 //_____________________________________________________________
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);
315  SetAttr("minRow",1);
316  SetAttr("UseTonkoClusterAnnotation",1);
317 }
318 //_____________________________________________________________
319 Int_t StTpcHitMaker::Init() {
320  LOG_INFO << "StTpcHitMaker::Init as\t" << GetName() << endm;
321  TString MkName(GetName());
322  for (Int_t k = 1; k < kAll; k++) {
323  if (MkName.CompareTo(Names[k],TString::kIgnoreCase) == 0) {kMode = (EMode) k; break;}
324  }
325  assert(kMode);
326  memset(maxHits,0,sizeof(maxHits));
327  maxBin0Hits = 0;
328  bin0Hits = 0;
329  return kStOK ;
330 }
331 //________________________________________________________________________________
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);}
338  f->cd();
339  if (oldToken >= 0) {
340  if (fSectCounts) {fSectCounts->Write(); delete fSectCounts;}
341  if (fAvLaser) {
342  for (Int_t s = 1; s <= 24; s++) {
343  if (fAvLaser[s-1]) {fAvLaser[s-1]->Write(); delete fAvLaser[s-1];}
344  }
345  delete [] fAvLaser; fAvLaser = 0;
346  }
347  }
348  enum {NoDim = 3};
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);
359  // fAvLaser[s-1]->CalculateErrors(kTRUE);
360  for (Int_t i = 0; i < NoDim; i++) {
361  fAvLaser[s-1]->GetAxis(i)->SetName(NameV[i]);
362  }
363  f->Add(fAvLaser[s-1]);
364  }
365 #else /* ! __USE__THnSparse__ */
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]);
378  }
379 #endif /* __USE__THnSparse__ */
380  oldToken = newToken;
381 }
382 //________________________________________________________________________________
383 Int_t StTpcHitMaker::InitRun(Int_t runnumber) {
384  SetAttr("maxRow",St_tpcPadConfigC::instance()->numberOfRows(20));
385  if (IAttr("Cosmics")) SetCosmics();
386 // Prepare scaled hit maxima
387 
388  // No hit maxima if these DB params are 0
389  Int_t maxHitsPerSector = St_tpcMaxHitsC::instance()->maxSectorHits();
390  Int_t maxBinZeroHits = St_tpcMaxHitsC::instance()->maxBinZeroHits();
391  Int_t livePads = 0;
392  Int_t totalPads = 0;
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;
404  }
405  livePads += liveSecPads;
406  totalPads += totalSecPads;
407  }
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;}
414  } else {
415  maxHits[sector-1] = 0;
416  if (Debug()) {LOG_INFO << "No maxHits in sector " << sector << endm;}
417  }
418  }
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;}
424  } else {
425  maxBin0Hits = 0;
426  if (Debug()) {LOG_INFO << "No maxBinZeroHits" << endm;}
427  }
428  // write event header for AvLaser
429  if (kMode == kTpxAvLaser || kMode == kTpcAvLaser) {
430  StEvtHddr *header = GetEvtHddr();
431  if (header) {
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);}
434  f->cd();
435  header->Write();
436  }
437  }
438  return kStOK;
439 }
440 //_____________________________________________________________
442  if (St_tpcStatusC::instance()->isDead()) {
443  LOG_WARN << "TPC status indicates it is unusable for this event. Ignoring hits." << endm;
444  return kStOK;
445  }
446 
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) {
452 #ifdef __TOKENIZED__
453  InitializeHistograms(Token());
454 #else
455  InitializeHistograms(0);
456 #endif
457  }
458  StMaker* maskMk = GetMakerInheritsFrom("StMtdTrackingMaskMaker");
459  unsigned int mask = (maskMk ? maskMk->UAttr("TpcSectorsByMtd") : ~0U); // 24 bit masking for sectors 1..24
460  bin0Hits = 0;
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; // sector masking
468  fId = 0;
469  // invoke tpcReader to fill the TPC DAQ sector structure
470  Int_t hitsAdded = 0;
471  Int_t kMin = kUnknown;
472  for (Int_t k = kStandardiTPC; k > kMin; k--) {
473  if (k > kLegacyTpx)
474  mQuery = Form("%s/%s[%i]",tpcDataNames[k],cldadc.Data(),sector);
475  else
476  mQuery = Form("%s[%i]",tpcDataNames[k],sector);
477  StRtsTable *daqTpcTable = GetNextDaqElement(mQuery);
478  if (! daqTpcTable) continue;
479  // Int_t Nrows = daqTpcTable->GetNRows();
480  // if (! Nrows) 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);
487  }
488  if (daqTpcTable->GetNRows()) {
489  fTpc = 0;
490  if (kReaderType == kLegacyTpx || kReaderType == kLegacyTpc) fTpc = (tpc_t*)*DaqDta()->begin();
491  Int_t row = RowNumber();
492  if (row >= minRow && row <= maxRow) {
493  if (Debug()) {
494  LOG_INFO << "StTpcHitMaker::Make(" << Names[kMode] << ") => " << tpcDataNames[k] << " for sector = " << sector << " row = " << row << endm;
495  }
496  switch (kMode) {
497  case kTpc:
498  case kiTPC:
499  case kTpxO:
500  case kTpx: hitsAdded += UpdateHitCollection(sector); break;
501  case kTpcPulser:
502  case kTpxPulser: if (fTpc) DoPulser(sector); break;
503  case kTpcAvLaser:
504  case kTpxAvLaser:
505  if ( fTpc) TpcAvLaser(sector);
506  else TpxAvLaser(sector);
507  fSectCounts->Fill(sector);
508  break;
509  case kTpcRaw:
510  case kTpxRaw:
511  case kiTPCRaw:
512  if ( fTpc) RawTpcData(sector);
513  else RawTpxData(sector);
514  break;
515  default:
516  break;
517  }
518  }
519  }
520  }
521  daqTpcTable = GetNextDaqElement(mQuery);
522  }
523  } // Loop over ReaderType
524  if (maxHits[sector-1] && hitsAdded > maxHits[sector-1]) {
525  LOG_ERROR << "Too many hits (" << hitsAdded << ") in one sector ("
526  << sector << "). Skipping event." << endm;
527  return kStSkip;
528  }
529  }
530  if (maxBin0Hits && bin0Hits > maxBin0Hits) {
531  LOG_ERROR << "Too many hits (" << bin0Hits
532  << ") starting at time bin 0. Skipping event." << endm;
533  return kStSkip;
534  }
535  if (kMode == kTpc || kMode == kTpx || kMode == kTpxO) { // || kMode == kiTPC) { --> no after burner for iTpc
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;}
539  StTpcHitCollection *hitCollection = pEvent->tpcHitCollection();
540  if (hitCollection && ! IAttr("NoTpxAfterBurner")) AfterBurner(hitCollection);
541  }
542  if (IAttr("CheckThrSeq") && (kMode == kTpcRaw || kMode == kTpxRaw || kMode == kiTPCRaw)) {
543  CheckThrSeq();
544  }
545  return kStOK;
546 }
547 //_____________________________________________________________
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);
556  }
557  }
558  for (Int_t s = 1; s <= 24; s++){
559  StTpcDigitalSector *digitalSector = GetDigitalSector(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);
566  if (! ntb) continue;
567  digitalSector->getTimeAdc(r,p,ADCs,IDTs);
568  Int_t adcMin = 127;
569  Int_t seq = 127;
570  Int_t tbF = -1, tbL = -1;
571  for (Int_t tb = 0; tb < __MaxNumberOfTimeBins__; tb++) {
572  if (! ADCs[tb]) {
573  if (tbF > -1 && tbL >= tbF) {
574  if (seq > tbL - tbF + 1) seq = tbL - tbF + 1;
575  }
576  tbF = tbL = -1;
577  continue;
578  }
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]))));
581  }
582  if (tbF < 0) tbF = tb;
583  tbL = tb;
584  }
585  if (adcMin < 127 && seq < 127) {
586  if (seq == 1 && adcMin == 1) {
587  static Int_t ibreak = 0;
588  ibreak++;
589  }
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);
594  }
595  }
596  }
597  }
598 }
599 //_____________________________________________________________
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]);
610  }
611  }
612  }
613 #endif /* __USE__THnSparse__ */
614  if (fThr && fSeq) {
615  Int_t adcMinFL = 127;
616  Int_t seqFL = 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) {
624  adcMinFL = adcMin;
625  seqFL = seq;
626  s1 = s2 = s;
627  r1 = r2 = r;
628  } else {
629  if (adcMinFL == adcMin && seqFL == seq) {
630  r2 = r;
631  s2 = s;
632  } else {
633  LOG_INFO << "StTpcHitMaker::Finish CheckThrSeq in sectors [" << s1 << "," << s2 << "] and rows[" << r1 << "," << r2 << "] Threshold = " << adcMinFL << " and sequence = " << seqFL << endm;
634  adcMinFL = adcMin;
635  seqFL = seq;
636  s1 = s2 = s;
637  r1 = r2 = r;
638  }
639  }
640  }
641  }
642  LOG_INFO << "StTpcHitMaker::Finish CheckThrSeq in sectors [" << s1 << "," << s2 << "] and rows[" << r1 << "," << r2 << "] Threshold = " << adcMinFL << " and sequence = " << seqFL << endm;
643  }
644  SafeDelete(fThr);
645  SafeDelete(fSeq);
646  return StMaker::Finish();
647 }
648 //_____________________________________________________________
649 Int_t StTpcHitMaker::UpdateHitCollection(Int_t sector) {
650  // Populate StEvent with StTpcHit collection
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;}
654  StTpcHitCollection *hitCollection = pEvent->tpcHitCollection();
655  if ( !hitCollection ) {
656  // Save the hit collection to StEvent...if needed
657  hitCollection = new StTpcHitCollection();
658  pEvent->setTpcHitCollection(hitCollection);
659  }
660  Int_t NRows = DaqDta()->GetNRows();
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) {
666  tpc_t *tpc = (tpc_t *) DaqDta()->GetTable();
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;
674  if (c->flags &&
675  (c->flags & ~(FCF_ONEPAD | FCF_MERGED | FCF_BIG_CHARGE))) continue;
676  if (kMode == kTpxO) c->flags |= 256; // mark cluster if it is coming from extra online maker
677  Int_t row = padrow + 1;
678  Float_t pad = c->p;
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);
683  assert(iok);
684  }
685  }
686  }
687  } else {
688  // kReaderType == kStandardTpx
689  daq_cld *cld = (daq_cld *) DaqDta()->GetTable();
690  if (Debug() > 1) {
691  LOG_INFO << Form("CLD sec %2d: row %2d: clusters: %3d",sec, row, NRows) << endm;
692  }
693  for (Int_t l = 0; l < NRows; l++, cld++) {
694  if (Debug() > 1) {
695  LOG_INFO << Form(" pad %f[%d:%d], tb %f[%d:%d], cha %d, fla 0x%X",//, Id %d, Q %d ",
696  cld->pad,
697  cld->p1,
698  cld->p2,
699  cld->tb,
700  cld->t1,
701  cld->t2,
702  cld->charge,
703  cld->flags
704  ) << endm;
705  }
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;
711  if (cld->flags &&
712  (cld->flags & ~(FCF_ONEPAD | FCF_MERGED | FCF_BIG_CHARGE))) continue;
713  if (kMode == kTpxO) cld->flags |= 256; // mark cluster if it is coming from extra online maker
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);
719  assert(iok);
720  }
721  }
722  Int_t nhits = hitCollection->numberOfHits() - nhitsBefore;
723  if (Debug()) {
724  LOG_INFO << " Total hits in Sector : row " << sector << " : " << row << " = " << nhits << endm;
725  }
726  return nhits;
727 }
728 //_____________________________________________________________
729 StTpcHit *StTpcHitMaker::CreateTpcHit(const tpc_cl &cluster, Int_t sector, Int_t row) {
730  // Create an instance of the StTpcHit from the tpcReader data
731 
732  Float_t pad = cluster.p;
733  Float_t time = cluster.t;
734  if (kReaderType == kLegacyTpx) time += 22; // remove Tonko's offset
735  static StTpcCoordinateTransform transform(gStTpcDb);
736  static StTpcLocalSectorCoordinate local;
737  static StTpcLocalCoordinate global;
738  StTpcPadCoordinate padcoord(sector, row, pad, time);
739  transform(padcoord,local,kFALSE);
740  transform(local,global);
741 
742  UInt_t hw = 1; // detid_tpc
743  hw += sector << 4; // (row/100 << 4); // sector
744  hw += row << 9; // (row%100 << 9); // row
745  static StThreeVector<double> hard_coded_errors(fgDp,fgDt,fgDperp);
746 #if 0
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();
749 #endif
750  Double_t q = 0; //cluster.charge * ((Double_t)St_tss_tssparC::instance()->ave_ion_pot() * (Double_t)St_tss_tssparC::instance()->scale())/(gain*wire_coupling) ;
751  StTpcHit *hit = StTpcHitFlag(global.position(),hard_coded_errors,hw,q
752  , (UChar_t ) 0 // c
753  , (Int_t) 0 // idTruth=0
754  , (UShort_t) 0 // quality=0,
755  , ++fId // id
756  , cluster.p1 // mnpad
757  , cluster.p2 // mxpad
758  , cluster.t1 // mntmbk
759  , cluster.t2 // mxtmbk
760  , pad
761  , time
762  , cluster.charge
763  ,cluster.flags);
764  if (hit->minTmbk() == 0) bin0Hits++;
765  if (Debug()) {
766  LOG_INFO << "StTpcHitMaker::CreateTpcHit fromt tpc_cl\t" <<*hit << endm;
767  }
768  return hit;
769 }
770 //_____________________________________________________________
771 StTpcHit *StTpcHitMaker::CreateTpcHit(const daq_cld &cluster, Int_t sector, Int_t row) {
772  // Create an instance of the StTpcHit from the tpcReader data
773 
774  Float_t pad = cluster.pad;
775  Float_t time = cluster.tb;
776 #if 1
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) ;
780 #else /* used in TFG till 07/31/20 */
781  Double_t q = 0;
782 #endif
783  // Check that slewing is active
784  static St_tpcSlewingC *tpcSlewing = St_tpcSlewingC::instance();
785  if (tpcSlewing && tpcSlewing->type() != 1001) tpcSlewing = 0;
786  if (tpcSlewing) {
787  // Correct for slewing (needs corrected q, and time in microsec)
788  Double_t freq = gStTpcDb->Electronics()->samplingFrequency();
789  time = freq * tpcSlewing->correctedT(sector,row,q,time/freq);
790  }
791  static StTpcCoordinateTransform transform(gStTpcDb);
792  static StTpcLocalSectorCoordinate local;
793  static StTpcLocalCoordinate global;
794  StTpcPadCoordinate padcoord(sector, row, pad, time);
795  transform(padcoord,local,kFALSE);
796  transform(local,global);
797 
798  UInt_t hw = 1; // detid_tpc
799  hw += sector << 4; // (row/100 << 4); // sector
800  hw += row << 9; // (row%100 << 9); // row
801  static StThreeVector<double> hard_coded_errors(fgDp,fgDt,fgDperp);
802  UShort_t flag = cluster.flags;
803  if (kMode == kTpxO) flag |= 256; // mark cluster if it is coming from extra online maker
804 
805  StTpcHit *hit = StTpcHitFlag(global.position(),hard_coded_errors,hw,q
806  , (UChar_t ) 0 // c
807  , (Int_t) 0 // idTruth=0
808  , (UShort_t) 0 // quality=0,
809  , ++fId // id =0
810  , cluster.p1 // mnpad
811  , cluster.p2 // mxpad
812  , cluster.t1 // mntmbk
813  , cluster.t2 // mxtmbk
814  , pad
815  , time
816  , cluster.charge
817  , flag);
818  if (hit->minTmbk() == 0) bin0Hits++;
819  if (Debug()) {
820  LOG_INFO << "StTpcHitMaker::CreateTpcHit fromt daq_cld\t" <<*hit << endm;
821  }
822  return hit;
823 }
824 //________________________________________________________________________________
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;
829  if (! pulserP) {
830  TFile *f = GetTFile();
831  assert(f);
832  f->cd();
833  pulserP = new TNtuple("pulserP","Pulser analysis",names);
834  }
835  Int_t r, p, tb, tbmax;
836  Int_t npeak, nnoise;
837  if (! fTpc) return;
838  if (! fTpc->channels_sector) return;
839  for(Int_t row = 1; row <= St_tpcPadConfigC::instance()->numberOfRows(sector); row++) {
840  r = row - 1;
841  if (! fTpc->cl_counts[r]) continue;
842  for (Int_t pad = 1; pad <= 182; pad++) {
843  p = pad - 1;
844  Int_t ncounts = fTpc->counts[r][p];
845  if (! ncounts) continue;
846  static UShort_t adc[512];
847  memset (adc, 0, sizeof(adc));
848  tbmax = 513;
849  UShort_t adcmax = 0;
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) {
854  tbmax = tb;
855  adcmax = adc[tb];
856  }
857  }
858  if (tbmax < 2 || tbmax > 504) continue;
859  npeak = nnoise = 0;
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);
864  Double_t peak = 0;
865  Double_t noise = 0;
866  Double_t t0 = 0;
867  for (Int_t i = i1; i <= i2; i++) {
868  if (i >= i1s && i <= i2s) continue;
869  nnoise++;
870  noise += adc[i];
871  }
872  if (nnoise) noise /= nnoise;
873  for (Int_t i = i1s; i <= i2s; i++) {
874  npeak++;
875  peak += adc[i] - noise;
876  t0 += i*(adc[i] - noise);
877  }
878  if (peak <= 0) continue;
879  t0 /= peak;
880  Pulser.sector = sector;
881  Pulser.row = row;
882  Pulser.pad = pad;
883  Pulser.gain = peak;
884  Pulser.t0 = t0;
885  Pulser.nnoise = nnoise;
886  Pulser.noise = noise;
887  Pulser.npeak = npeak;
888  pulserP->Fill(&Pulser.sector);
889  }
890  }
891 }
892 //________________________________________________________________________________
893 void StTpcHitMaker::TpcAvLaser(Int_t sector) {
894  if (! fTpc || !fAvLaser) return;
895  Int_t npixels = 0;
896  struct pixl_t {
897  Double_t sector, row, pad, time;
898  };
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]);
906  }
907 #endif /* __USE__THnSparse__ */
908  pixl_t pixel;
909  pixel.sector = sector;
910  for(Int_t r = 0; r < St_tpcPadConfigC::instance()->numberOfRows(sector); r++) {
911  pixel.row = r+1;
912  for (Int_t pad = 1; pad <= 182; pad++) {
913  pixel.pad = pad;
914  Int_t p = pad - 1;
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);
924 #else /* ! __USE__THnSparse__ */
925  fAvLaser[sector-1]->Fill(pixel.row,pixel.pad,pixel.time,adc);
926 #endif /* __USE__THnSparse__ */
927  npixels++;
928  }
929  }
930  }
931  LOG_INFO << " Total pixels in Sector : " << sector << " = " << npixels << endm;
932 }
933 //________________________________________________________________________________
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]);
943  }
944 #endif /* __USE__THnSparse__ */
945  Int_t r=RowNumber() ; // I count from 1
946  if(r==0) return; // TPC does not support unphy. rows so we skip em
947  r-- ; // TPC wants from 0
948  Int_t p = Pad() - 1 ; // ibid.
949  if (p < 0 || p >= St_tpcPadConfigC::instance()->padsPerRow(sector,r+1)) return;
950  struct pixl_t {
951  Double_t sector, row, pad, time;
952  };
953  pixl_t pixel;
954  pixel.sector = sector;
955  pixel.row = r+1;
956  pixel.pad = p+1;
957  TGenericTable::iterator iword = DaqDta()->begin();
958  for (;iword != DaqDta()->end();++iword) {
959  daq_adc_tb &daqadc = (*(daq_adc_tb *)*iword);
960  Int_t tb = daqadc.tb;
961  pixel.time = tb;
962  Double_t adc = daqadc.adc;
963 #ifdef __NOT_ZERO_SUPPRESSED_DATA__
964 #ifdef __TOKENIZED__
965  if (tb >= 368 && tb <= 383)
966  // adc -= St_tpcPedestalC::instance()->Pedestal(sector,r+1,p+1,tb);
967  adc -= St_tpcPedestalC::instance()->Pedestal(sector,r+1,p+1);
968 #else /* ! __TOKENIZED__ */
969  // adc -= St_tpcPedestalC::instance()->Pedestal(sector,r+1,p+1,tb);
970  adc -= St_tpcPedestalC::instance()->Pedestal(sector,r+1,p+1);
971 #endif /* __TOKENIZED__ */
972 #else
973  // if (adc < 6) continue;
974 #endif /* __NOT_ZERO_SUPPRESSED_DATA__ */
975 #ifdef __USE__THnSparse__
976  fAvLaser[sector-1]->Fill(&pixel.row,adc);
977 #else /* ! __USE__THnSparse__ */
978  fAvLaser[sector-1]->Fill(pixel.row,pixel.pad,pixel.time,adc);
979 #endif /* __USE__THnSparse__ */
980  }
981 }
982 //________________________________________________________________________________
983 void StTpcHitMaker::DumpPixels2Ntuple(Int_t sector, Int_t row, Int_t pad) {
984  struct BPoint_t {
985  Float_t event, sector, row, pad, tb, adc, idt;
986  };
987  static const Char_t *BName = "event:sector:row:pad:tb:adc:idt";
988  static TNtuple *adcP = 0;
989  if (! adcP && GetTFile() ) {
990  GetTFile()->cd();
991  adcP = new TNtuple("adcP","Pulser ADC",BName);
992  }
993  if (! adcP) return;
994  static BPoint_t P;
995  P.event = GetEventNumber();
996  P.sector = sector;
997  P.row = row;
998  P.pad = pad;
999  for (Int_t i = 0; i < __MaxNumberOfTimeBins__; i++) {
1000  if (! ADCs[i]) continue;
1001  P.tb = i;
1002  P.adc = ADCs[i];
1003  P.idt = IDTs[i];
1004  adcP->Fill(&P.event);
1005  }
1006 }
1007 //________________________________________________________________________________
1008 void StTpcHitMaker::PrintSpecial(Int_t sector) {
1009  // example usage: calculate total charge and
1010  // print occupancy
1011  Int_t r,p,t ;
1012  UInt_t adc = 0;
1013  UChar_t val ;
1014  if (! fTpc) return;
1015  if(fTpc->mode==0) { // normal event
1016  UInt_t tot_pix = 0 ;
1017  UInt_t cl_count = 0 ;
1018  Int_t i ;
1019 
1020  for(r=0;r<St_tpcPadConfigC::instance()->numberOfRows(sector);r++) { // padrow
1021  for(p=0;p<182;p++) { // pad
1022  for(t=0;t<fTpc->counts[r][p];t++) {
1023  val = fTpc->adc[r][p][t] ;
1024  Int_t vali = log8to10_table[val];
1025  adc += val ;
1026  if(val) tot_pix++ ;
1027  if (Debug() > 1) {
1028  Int_t timebin = fTpc->timebin[r][p][t] ;
1029  printf("%d %d %d %d %d\n",sector,r+1,p+1,timebin,vali) ;
1030  }
1031  }
1032  }
1033 
1034  if(fTpc->has_clusters) {
1035  cl_count += fTpc->cl_counts[r] ;
1036  }
1037  if (Debug() > 1) {
1038  if(fTpc->has_clusters) {
1039  for(i=0;i<fTpc->cl_counts[r];i++) {
1040  tpc_cl *c = &fTpc->cl[r][i] ;
1041 
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) ;
1044  }
1045  }
1046  }
1047  }
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;
1051  }
1052 }
1053 //________________________________________________________________________________
1054 StTpcDigitalSector *StTpcHitMaker::GetDigitalSector(Int_t sector) {
1055  TDataSet *event = GetData("Event");
1056  StTpcRawData *data = 0;
1057  if (! event) {
1058  data = new StTpcRawData(24);
1059  event = new TObjectSet("Event", data);
1060  AddData(event);
1061  } else data = (StTpcRawData *) event->GetObject();
1062  assert(data);
1063  StTpcDigitalSector *digitalSector = data->GetSector(sector);
1064  if (! digitalSector) {
1065  digitalSector = new StTpcDigitalSector(sector);
1066  data->setSector(sector,digitalSector);
1067  }
1068  return digitalSector;
1069 }
1070 //________________________________________________________________________________
1071 Int_t StTpcHitMaker::RawTpxData(Int_t sector) {
1072  static Int_t TonkoAnn = IAttr("UseTonkoClusterAnnotation");
1073  Short_t ADCs2[512];
1074 #ifndef __TFG__VERSION__
1075  UShort_t IDTs2[512];
1076 #else
1077  Int_t IDTs2[512];
1078 #endif
1079  memset(ADCs, 0, sizeof(ADCs));
1080  memset(IDTs, 0, sizeof(IDTs));
1081  StTpcDigitalSector *digitalSector = 0;
1082  Int_t r_old = -1;
1083  Int_t p_old = -1;
1084  Int_t Total_data = 0;
1085  Int_t r=RowNumber() ; // I count from 1
1086  if(r==0) return 0 ; // TPC does not support unphysical rows so we skip them
1087  Int_t p = Pad() - 1 ; // ibid.
1088  r-- ; // TPC wants from 0
1089  if (p < 0 || p >= St_tpcPadConfigC::instance()->padsPerRow(sector,r+1)) return 0;
1090  TGenericTable::iterator iword = DaqDta()->begin();
1091  if (iword == DaqDta()->end()) return 0;
1092  Int_t some_data = 0;
1093  do {
1094  if (some_data) {
1095  Total_data += some_data;
1096  some_data = 0;
1097  if (! digitalSector) digitalSector = GetDigitalSector(sector);
1098  Int_t ntbold = digitalSector->numberOfTimeBins(r_old+1,p_old+1);
1099  if (ntbold) {
1100  if (Debug()) {
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;
1103  }
1104  digitalSector->getTimeAdc(r_old+1,p_old+1,ADCs2,IDTs2);
1105  // if (IAttr("UseTonkoClusterAnnotation")) {
1106  if (TonkoAnn) {
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];
1111  }
1112  }
1113  }
1114  digitalSector->putTimeAdc(r_old+1,p_old+1,ADCs,IDTs);
1115  if (IAttr("TpxDumpPxls2Nt")) {
1116  DumpPixels2Ntuple(sector,r_old+1,p_old+1);
1117  }
1118  memset(ADCs, 0, sizeof(ADCs));
1119  memset(IDTs, 0, sizeof(IDTs));
1120  }
1121  if (r_old != r || p_old != p) {
1122  r_old = r;
1123  p_old = p;
1124  }
1125  for (;iword != DaqDta()->end();++iword) {
1126  daq_adc_tb &daqadc = (*(daq_adc_tb *)*iword);
1127  Int_t tb = daqadc.tb;
1128  Int_t adc = daqadc.adc;
1129 #ifdef __NOT_ZERO_SUPPRESSED_DATA__
1130  // adc -= St_tpcPedestalC::instance()->Pedestal(sector,r+1,p+1,tb);
1131  adc -= St_tpcPedestalC::instance()->Pedestal(sector,r+1,p+1);
1132  if (adc <= 0) continue;
1133 #endif
1134  ADCs[tb] = adc;
1135  // if (IAttr("UseTonkoClusterAnnotation")) {
1136  if (TonkoAnn) {
1137  IDTs[tb] = 65535;
1138  }
1139  some_data++ ; // I don't know the bytecount but I'll return something...
1140  }
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);
1145  }
1146  return Total_data;
1147 }
1148 //________________________________________________________________________________
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));
1154  StTpcDigitalSector *digitalSector = GetDigitalSector(sector);
1155  assert( digitalSector );
1156  Int_t Total_data = 0;
1157  for (Int_t row = 1; row <= digitalSector->numberOfRows(); row++) {
1158  Int_t r = row - 1;
1159  for (Int_t pad = 1; pad <= digitalSector->numberOfPadsAtRow(row); pad++) {
1160  Int_t p = pad - 1;
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]];
1168  // if (IAttr("UseTonkoClusterAnnotation")) {
1169  if (TonkoAnn) {
1170  IDTs[tb] = 65535;
1171  }
1172  Total_data++;
1173  }
1174  Int_t ntbold = digitalSector->numberOfTimeBins(row,pad);
1175  if (ntbold) {
1176 #if 0
1177  LOG_INFO << "digitalSector " << sector
1178  << " already has " << ntbold << " at row/pad " << row << "/" << pad << endm;
1179 #endif
1180  }
1181  digitalSector->putTimeAdc(row,pad,ADCs,IDTs);
1182  if (IAttr("TpxDumpPxls2Nt")) {
1183  DumpPixels2Ntuple(sector,row,pad);
1184  }
1185  }
1186  }
1187  if (Total_data) {
1188  LOG_INFO << "Read " << Total_data << " pixels from Sector " << sector << endm;
1189  }
1190  return Total_data;
1191 }
1192 //________________________________________________________________________________
1193 #ifdef __MAKE_NTUPLE__
1194 static TNtuple *tup = 0;
1195 struct TpcHitPair_t {
1196  Float_t sec, row,
1197  qK, padK, tbkK, padKmn, padKmx, tbkKmn, tbkKmx, IdTK, QAK,
1198  qL, padL, tbkL, padLmn, padLmx, tbkLmn, tbkLmx, IdTL, QAL,
1199  padOv, tbkOv;
1200 };
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:"
1204  "padOv:tbkOv";
1205 #endif
1206 //________________________________________________________________________________
1207 Bool_t TpcHitLess(const StTpcHit *lhs, const StTpcHit *rhs) {
1208  return (200*lhs->timeBucket() + lhs->pad()) < (200*rhs->timeBucket() + rhs->pad());
1209 };
1210 //________________________________________________________________________________
1211 void StTpcHitMaker::AfterBurner(StTpcHitCollection *TpcHitCollection) {
1212  static Float_t padDiff = 2.5, timeBucketDiff = 5.0;
1213  static StTpcCoordinateTransform transform(gStTpcDb);
1214  static StTpcLocalSectorCoordinate local;
1215  static StTpcLocalCoordinate global;
1216  if (! TpcHitCollection) return;
1217 
1218 #ifdef __MAKE_NTUPLE__
1219  if (! tup) {
1220  if (StChain::GetChain()->GetTFile()) {
1221  StChain::GetChain()->GetTFile()->cd();
1222  tup = new TNtuple("HitT","Cluster Pair Info",vTpcHitMRPair);
1223  }
1224  }
1225  TpcHitPair_t pairC;
1226 #endif
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++) {
1232  StTpcSectorHitCollection* sectorCollection = TpcHitCollection->sector(sec-1);
1233  if (sectorCollection) {
1234  UInt_t numberOfPadrows = sectorCollection->numberOfPadrows();
1235  for (UInt_t row = 1; row <= numberOfPadrows; row++) {
1236  StTpcPadrowHitCollection *rowCollection = sectorCollection->padrow(row-1);
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(),
1243  TpcHitLess);
1244  // Merge splitted clusters
1245  Int_t merged = 0;
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;
1253  break;
1254  }
1255  newlyMerged = 0;
1256  for (UInt_t k = 0; k < NoHits; k++) {
1257  if (mergedHits[k] != mergePass) continue; // skip any hits not merged in previous pass
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__
1263  pairC.sec = sec;
1264  pairC.row = row;
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();
1274 #endif
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;
1280  // Most stringent tests first to reduce calls
1281  // check hits near by
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;
1286 #endif
1287  // Are extents overlapped ?
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__
1295  if (tup) {
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);
1308  }
1309  if (notNear) continue;
1310 #endif
1311 #ifdef FCF_CHOPPED
1312  UShort_t flag = lHit->flag() | FCF_CHOPPED;
1313 #else
1314  UShort_t flag = lHit->flag() | 0x080;
1315 #endif
1316  lHit->setFlag(flag);
1317  RejNoHits++;
1318  if (_debug) {
1319  cout << "mk" << k; kHit->Print();
1320  cout << "ml" << l; lHit->Print();
1321  }
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();
1334  } else {
1335  if (lHit->adc()*lHit->qaTruth() < kHit->adc()*kHit->qaTruth()) {
1336  QA = kHit->adc()*kHit->qaTruth();
1337  IdT = kHit->idTruth();
1338  }
1339  }
1340  QA = QA/adc;
1341  kHit->setIdTruth(IdT, TMath::Nint(QA));
1342  kHit->setCharge(q);
1343  kHit->setExtends(pad,timeBucket,minPad,maxPad,minTmbk,maxTmbk);
1344  kHit->setAdc(adc);
1345  StTpcPadCoordinate padcoord(sec, row, pad, timeBucket);
1346  transform(padcoord,local,kFALSE);
1347  transform(local,global);
1348  kHit->setPosition(global.position());
1349  if (_debug) {
1350  cout << "m " << k; kHit->Print();
1351  }
1352  merged++;
1353  newlyMerged++;
1354  mergedHits[k] = mergePass + 1;
1355  break; // done with hit k for now, will look again later
1356  } // l hit loop
1357  } // k hit loop
1358  mergePass++;
1359  } // merging pass
1360 #ifdef __CORRECT_S_SHAPE__
1361  // Correct S - shape in pad direction
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();
1367  Int_t io = 1;
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);
1371  StTpcPadCoordinate padcoord(sec, row, pad, timeBucket);
1372  transform(padcoord,local,kFALSE);
1373  transform(local,global);
1374  kHit->setPosition(global.position());
1375  } // k hit loop
1376 #endif
1377  } // row collection exists
1378  } // rows
1379  } // sector collection exists
1380  } // sectors
1381  LOG_INFO << "StTpcHitMaker::AfterBurner from " << TotNoHits << " Total no. of hits " << RejNoHits << " were rejected" << endm;
1382  return;
1383 };
1384 //________________________________________________________________________________
1385 StTpcHit* StTpcHitMaker::StTpcHitFlag(const StThreeVectorF& p,
1386  const StThreeVectorF& e,
1387  UInt_t hw, float q, UChar_t c,
1388  Int_t idTruth, UShort_t quality,
1389  UShort_t id,
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,
1392  UShort_t flag) {
1393  // New hit
1394  StTpcHit* hit = new StTpcHit(p,e,hw,q,c,idTruth,quality,id,mnpad,mxpad,mntmbk,mxtmbk,cl_x,cl_t,adc);
1395  // Check for sanity
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;
1402 #if 0
1403  if (fgCosmics && TMath::Abs(p.z()) > 150) flag |= FCF_SANITY;
1404 #endif
1405  hit->setFlag(flag);
1406  assert(! ( TMath::IsNaN(hit->position().x()) || TMath::IsNaN(hit->position().y()) || TMath::IsNaN(hit->position().x())));
1407  return hit;
1408 }
1409 #ifdef __USE__THnSparse__
1410 //________________________________________________________________________________
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);//, hist->GetChunkSize());
1417  hnew->CalculateErrors(kTRUE);
1418  delete [] nbins;
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());
1423  }
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__
1428  Double_t max = -1;
1429  for (Long64_t i = 0; i < N; ++i) {
1430  Double_t cont = hist->GetBinContent(i, bins);
1431  if (cont > max) max = cont;
1432  }
1433  for (Long64_t i = 0; i < N; ++i) {
1434  Double_t cont = hist->GetBinContent(i, bins);
1435  if (cont < max/compress) continue;
1436  // Long64_t bin = hnew->GetBin(bins);
1437  for (Int_t d = 0; d < nd; ++d) {x[d] = hist->GetAxis(d)->GetBinCenter(bins[d]);}
1438  hnew->Fill(x,cont);
1439  }
1440 #else
1441  for (Long64_t i = 0; i < N; ++i) {
1442  Double_t cont = hist->GetBinContent(i, bins);
1443  if (cont < 1.0) continue;
1444  // Long64_t bin = hnew->GetBin(bins);
1445  for (Int_t d = 0; d < nd; ++d) {x[d] = hist->GetAxis(d)->GetBinCenter(bins[d]);}
1446  hnew->Fill(x,cont);
1447  }
1448 #endif
1449  delete [] bins;
1450  delete [] x;
1451  cout << hnew->GetName() << " has " << hnew->GetNbins() << " bins after compression." << endl;
1452  return hnew;
1453 }
1454 #endif /* __USE__THnSparse__ */
1455 //________________________________________________________________________________
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) {
1462  row += 41-14;
1463  }
1464  }
1465  return row;
1466 }
virtual Int_t Finish()
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.
Definition: StMaker.cxx:332
virtual Int_t Write(const char *name=0, Int_t option=0, Int_t bufsize=0)
Definition: TDataSet.cxx:893
Definition: Stypes.h:49
virtual Long_t GetNRows() const
Returns the number of the used rows for the wrapped table.
Definition: TTable.cxx:1388
StRtsTable * DaqDta()
Return the current DAQ data block. This member function is provided for convenience.
Definition: daq_tpc.h:22
Definition: daq_tpc.h:11
virtual const char * GetName() const
special overload
Definition: StMaker.cxx:237
Definition: Stypes.h:42
Definition: Stypes.h:40
static UInt_t Token()
current token
virtual Int_t Finish()
Definition: StMaker.cxx:776
virtual Char_t * Print(Char_t *buf, Int_t n) const
Create IDL table defintion (to be used for XDF I/O)
Definition: TTable.cxx:1548