StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StFstCalibrationMaker.cxx
1 /***************************************************************************
2  * $Id: StFstCalibrationMaker.cxx$
3  *
4  * Author: Shenghui Zhang, Oct. 2021
5  ****************************************************************************
6  * Description:
7  * See header file.
8  ****************************************************************************
9  ****************************************************************************
10  * StFstCalibrationMaker.cxx,v 1.0
11  * Revision 1.0 2021/10/15 Shenghui Zhang
12  * Initial version
13  ****************************************************************************/
14 
15 #include <string>
16 #include "StEvent.h"
17 #include "StEvent/StEnumerations.h"
18 #include "StEvent/StFstConsts.h"
19 #include "StIOMaker/StIOMaker.h"
20 #include "StFstCalibrationMaker.h"
21 #include "StRoot/StFstUtil/StFstCollection.h"
22 #include "StRoot/StFstUtil/StFstRawHitCollection.h"
23 #include "StEvent/StFstRawHit.h"
24 #include "StRoot/StFstDbMaker/StFstDb.h"
25 
26 #include "tables/St_fstMapping_Table.h"
27 #include "tables/St_fstControl_Table.h"
28 
29 #include "TMath.h"
30 #include "TH2S.h"
31 
32 const string StFstCalibrationMaker::sectionLabel[72]={ "1I","1O","2I","2O","3I","3O","4I","4O","5I","5O","6I","6O","7I","7O","8I","8O","9I","9O","10I","10O","11I","11O","12I","12O","13I","13O","14I","14O","15I","15O","16I","16O","17I","17O","18I","18O","19I","19O","20I","20O","21I","21O","22I","22O","23I","23O","24I","24O","25I","25O","26I","26O","27I","27O","28I","28O","29I","29O","30I","30O","31I","31O","32I","32O","33I","33O","34I","34O","35I","35O","36I","36O"};
33 
34 // constructor
35 StFstCalibrationMaker::StFstCalibrationMaker( const char* name ) : StMaker( name ), mTimeBinMask(0xFF), mRunHist(true), mDoPedCut(true), evtIdx(0), mHasFinished(0), mFstDb(0), mDoOutput(false) {
36  for(unsigned char iTb=0; iTb<kFstNumTimeBins; iTb++) {
37  hist_meanPed[iTb] = NULL;
38  hist_rmsPed[iTb] = NULL;
39  hist_cmNoise[iTb] = NULL;
40  hist_ranNoise[iTb] = NULL;
41  hist_sumPed[iTb] = NULL;
42  hist_sumRms[iTb] = NULL;
43  hist_sumCmn[iTb] = NULL;
44  hist_sumRan[iTb] = NULL;
45  hist_adcSpectrum[iTb] = NULL;
46  }
47 
48  mPedVec1stLoop.resize(kFstNumTimeBins*kFstNumElecIds);
49  mRmsVec1stLoop.resize(kFstNumTimeBins*kFstNumElecIds);
50  mRanVec1stLoop.resize(kFstNumTimeBins*kFstNumElecIds);
51 
52  mMathPedVec.resize(kFstNumTimeBins*kFstNumElecIds);
53  mMathRmsVec.resize(kFstNumTimeBins*kFstNumElecIds);
54  mMathCouVec.resize(kFstNumTimeBins*kFstNumElecIds);
55 
56  mMathRanVec.resize(kFstNumTimeBins*kFstNumElecIds);
57  mMathPedRanVec.resize(kFstNumTimeBins*kFstNumElecIds);
58  mMathRmsRanVec.resize(kFstNumTimeBins*kFstNumElecIds);
59  mMathCouRanVec.resize(kFstNumTimeBins*kFstNumElecIds);
60 
61  if(mRunHist)
62  {
63  mHistPedVec.assign( kFstNumTimeBins*kFstNumElecIds, (TH1S*)0 );
64  mHistRanVec.assign( kFstNumTimeBins*kFstNumElecIds, (TH1F*)0 );
65  }
66 
67  mHistCmnVec.assign( kFstNumTimeBins*kFstNumApvs*kFstNumRStripsPerSensor, (TH1F*)0 );
68  mPedVec.resize( kFstNumTimeBins*kFstNumElecIds ); // set to all zeros
69  mCmnVec.resize( kFstNumTimeBins*kFstNumApvs*kFstNumRStripsPerSensor ); // set to all zeros
70 
71  mMappingVec.resize( kFstNumElecIds );
72 };
73 
74 // deconstructor
75 StFstCalibrationMaker::~StFstCalibrationMaker()
76 {
77  mPedVec1stLoop.clear();
78  mRmsVec1stLoop.clear();
79  mRanVec1stLoop.clear();
80 
81  mMathPedVec.clear();
82  mMathRmsVec.clear();
83  mMathCouVec.clear();
84 
85  mMathRanVec.clear();
86  mMathPedRanVec.clear();
87  mMathRmsRanVec.clear();
88  mMathCouRanVec.clear();
89 
90  if(mRunHist)
91  {
92  while (!mHistPedVec.empty()) delete mHistPedVec.back(), mHistPedVec.pop_back();
93  while (!mHistRanVec.empty()) delete mHistRanVec.back(), mHistRanVec.pop_back();
94  }
95 
96  while (!mHistCmnVec.empty()) delete mHistCmnVec.back(), mHistCmnVec.pop_back();
97 
98  mPedVec.clear();
99  mCmnVec.clear();
100  mMappingVec.clear();
101 };
102 
103 // initialize
104 Int_t StFstCalibrationMaker::Init()
105 {
106  Int_t ierr = kStOk;
107 
108  char buffer[100];
109  for(unsigned char iTB=0; iTB<kFstNumTimeBins; iTB++) {
110  sprintf(buffer,"hist_meanPedestal_TimeBin%d",iTB);
111  hist_meanPed[iTB] = new TH1F(buffer, buffer, kFstNumElecIds, 0, kFstNumElecIds);
112  hist_meanPed[iTB]->SetStats(false);
113  hist_meanPed[iTB]->GetXaxis()->SetTitle("Channel Geometry Index");
114  hist_meanPed[iTB]->GetXaxis()->SetNdivisions(-72,false); //sections
115  hist_meanPed[iTB]->GetYaxis()->SetTitle("ADC count");
116 
117  sprintf(buffer,"hist_rmsPedestal_TimeBin%d",iTB);
118  hist_rmsPed[iTB] = new TH1F(buffer, buffer, kFstNumElecIds, 0, kFstNumElecIds);
119  hist_rmsPed[iTB]->SetStats(false);
120  hist_rmsPed[iTB]->GetXaxis()->SetTitle("Channel Geometry Index");
121  hist_rmsPed[iTB]->GetXaxis()->SetNdivisions(-72,false);//sections
122  hist_rmsPed[iTB]->GetYaxis()->SetTitle("ADC count");
123 
124  sprintf(buffer,"hist_cmNoise_TimeBin%d",iTB);
125  hist_cmNoise[iTB] = new TH1F(buffer, buffer, kFstNumApvs*kFstNumRStripsPerSensor, 0, kFstNumApvs*kFstNumRStripsPerSensor);
126  hist_cmNoise[iTB]->SetStats(false);
127  hist_cmNoise[iTB]->GetXaxis()->SetTitle("APV Geometry Index");
128  hist_cmNoise[iTB]->GetXaxis()->SetNdivisions(-72,false);//sections
129  hist_cmNoise[iTB]->GetYaxis()->SetTitle("ADC count");
130 
131  sprintf(buffer,"hist_ranNoise_TimeBin%d",iTB);
132  hist_ranNoise[iTB] = new TH1F(buffer, buffer, kFstNumElecIds, 0, kFstNumElecIds);
133  hist_ranNoise[iTB]->SetStats(false);
134  hist_ranNoise[iTB]->GetXaxis()->SetTitle("Channel Geometry Index");
135  hist_ranNoise[iTB]->GetXaxis()->SetNdivisions(-72,false);//sections
136  hist_ranNoise[iTB]->GetYaxis()->SetTitle("ADC count");
137 
138  for(int i=0; i<72; i++) {
139  TString binBuffer = "";
140  binBuffer = sectionLabel[i];
141  hist_meanPed[iTB]->GetXaxis()->SetBinLabel(i*512+256, binBuffer);
142  hist_rmsPed[iTB]->GetXaxis()->SetBinLabel(i*512+256, binBuffer);
143  hist_cmNoise[iTB]->GetXaxis()->SetBinLabel(i*4+4, binBuffer);
144  hist_ranNoise[iTB]->GetXaxis()->SetBinLabel(i*512+256, binBuffer);
145  }
146 
147  sprintf(buffer,"hist_sumPedestal_TimeBin%d",iTB);
148  hist_sumPed[iTB] = new TH1F(buffer, buffer, 128, 0, 4096);
149  hist_sumPed[iTB]->SetStats(kTRUE);
150  hist_sumPed[iTB]->GetXaxis()->SetTitle("Pedestal [ADC counts]");
151  hist_sumPed[iTB]->GetYaxis()->SetTitle("Counts");
152 
153  sprintf(buffer,"hist_sumRmsNoise_TimeBin%d",iTB);
154  hist_sumRms[iTB] = new TH1F(buffer, buffer, 128, 0, 256);
155  hist_sumRms[iTB]->SetStats(kTRUE);
156  hist_sumRms[iTB]->GetXaxis()->SetTitle("RMS Noise [ADC counts]");
157  hist_sumRms[iTB]->GetYaxis()->SetTitle("Counts");
158 
159  sprintf(buffer,"hist_sumCmNoise_TimeBin%d",iTB);
160  hist_sumCmn[iTB] = new TH1F(buffer, buffer, 128, 0, 256);
161  hist_sumCmn[iTB]->SetStats(kTRUE);
162  hist_sumCmn[iTB]->GetXaxis()->SetTitle("CM Noise [ADC counts]");
163  hist_sumCmn[iTB]->GetYaxis()->SetTitle("Counts");
164 
165  sprintf(buffer,"hist_sumRanNoise_TimeBin%d",iTB);
166  hist_sumRan[iTB] = new TH1F(buffer, buffer, 128, 0, 256);
167  hist_sumRan[iTB]->SetStats(kTRUE);
168  hist_sumRan[iTB]->GetXaxis()->SetTitle("RAN Noise [ADC counts]");
169  hist_sumRan[iTB]->GetYaxis()->SetTitle("Counts");
170 
171  sprintf(buffer,"hist_adcSpectrum_TimeBin%d",iTB);
172  hist_adcSpectrum[iTB] = new TH2S(buffer, buffer, kFstNumElecIds, 0, kFstNumElecIds, 512, 0, 4096);
173  hist_adcSpectrum[iTB]->SetStats(false);
174  hist_adcSpectrum[iTB]->GetXaxis()->SetTitle("Channel Electronics Index");
175  hist_adcSpectrum[iTB]->GetYaxis()->SetTitle("ADC counts");
176  }
177 
178  //read the pedestal/rms value from 1st loop for signal-like channel excluding in 2nd loop
179  if(mDoPedCut) {
180  TString infile = "fstPedNoiseTable.dat";
181  std::ifstream in(infile);
182  if (!in.is_open()) {
183  mDoPedCut=false;
184  LOG_WARN << "Could not find fstPedNoiseTable.dat! Set mDoPedCut to false" << endm;
185  } else {
186  LOG_INFO << "Read Pedestal and RMS from fstPedNoiseTable.dat!"<< endm;
187  int chElecId, rdo, arm, apv, chan, chTimeBin, chCode;
188  float chPed, chRms, chRan;
189  while(!in.eof()) {
190  in >> chElecId >> rdo >> arm >> apv >> chan >> chTimeBin >> chPed >> chRms >> chRan;
191 
192  chCode = kFstNumTimeBins * chElecId + chTimeBin;
193  mPedVec1stLoop[chCode] = chPed;
194  mRmsVec1stLoop[chCode] = chRms;
195  mRanVec1stLoop[chCode] = chRan;
196  }
197  }
198  in.close();
199  }
200 
201  return ierr;
202 };
203 
204 Int_t StFstCalibrationMaker::InitRun(Int_t runnumber)
205 {
206  Int_t ierr = kStOk;
207 
208  TObjectSet *fstDbDataSet = (TObjectSet *)GetDataSet("fst_db");
209  if (fstDbDataSet) {
210  mFstDb = (StFstDb *)fstDbDataSet->GetObject();
211  assert(mFstDb);
212  }
213  else {
214  LOG_ERROR << "InitRun : no fstDb" << endm;
215  return kStErr;
216  }
217 
218  //control parameter
219  const fstControl_st *fstControlTable = mFstDb->getControl();
220  if (!fstControlTable) {
221  LOG_ERROR << "Pointer to FST control table is null" << endm;
222  ierr = kStErr;
223  }
224  else
225  mPedCut = fstControlTable[0].kFstPedCutDefault;
226 
227  //channel mapping
228  const fstMapping_st *gM = mFstDb->getMapping();
229  if( !gM ) {
230  LOG_ERROR << "Pointer to FST mapping table is null" << endm;
231  ierr = kStErr;
232  }
233  else {
234  for(int i=0; i<kFstNumElecIds; i++) {
235  LOG_DEBUG<<Form(" Print entry %d : geoId=%d ",i,gM[0].mapping[i])<<endm;
236  mMappingVec[i] = gM[0].mapping[i];
237  }
238  }
239 
240  return ierr;
241 };
242 
243 Int_t StFstCalibrationMaker::Make()
244 {
245  Int_t ierr = kStOk;
246  //obtain raw hit information from temporary dataset
247  TObjectSet* fstDataSet = (TObjectSet*)GetDataSet("fstRawHitAndCluster");
248  if (! fstDataSet) {
249  LOG_WARN << "Make() - there is no fstDataSet (raw hits container) " << endm;
250  ierr = kStWarn;
251  }
252 
253  StFstCollection* fstCollectionPtr = (StFstCollection*)fstDataSet->GetObject();
254  if(!fstCollectionPtr) {
255  LOG_WARN << "Make() - no fstCollection."<<endm;
256  ierr = kStWarn;
257  }
258 
259  if( !ierr ){
260  if( evtIdx%10 == 0 )
261  cout << "event index: " << evtIdx << endl;
262 
263  std::stringstream ss, sc;
264  int sumAdcPerEvent[kFstApvsPerWedge][kFstNumRStripsPerSensor][kFstNumTimeBins];
265  int channelCountsPerEvent[kFstApvsPerWedge][kFstNumRStripsPerSensor][kFstNumTimeBins];
266  double cmNoisePerEvent[kFstApvsPerWedge][kFstNumRStripsPerSensor][kFstNumTimeBins];
267 
268  for( unsigned char wedgeIdx=0; wedgeIdx<kFstNumWedges; ++wedgeIdx ){
269  StFstRawHitCollection *rawHitCollectionPtr = fstCollectionPtr->getRawHitCollection( wedgeIdx );
270 
271  for(int i=0; i<kFstApvsPerWedge; i++) {
272  for(int j=0; j<kFstNumRStripsPerSensor; j++) {
273  for(int k=0; k<kFstNumTimeBins; k++) {
274  sumAdcPerEvent[i][j][k] = 0;
275  channelCountsPerEvent[i][j][k] = 0;
276  cmNoisePerEvent[i][j][k] = 0.;
277  }
278  }
279  }
280 
281  if( rawHitCollectionPtr ){
282  std::vector<StFstRawHit*>& rawHitVec = rawHitCollectionPtr->getRawHitVec();
283  std::vector< StFstRawHit* >::iterator rawHitIter;
284 
285  for( rawHitIter = rawHitVec.begin(); rawHitIter != rawHitVec.end(); ++rawHitIter ){
286  int elecId = (*rawHitIter)->getChannelId();
287  int geoId = (*rawHitIter)->getGeoId(); //channel geometry Id, counting from 0 to 36863
288  int apvId = elecId/kFstNumApvChannels; //APV chip geometry Id, counting from 0 to 287
289  int ChPerRdo = kFstNumArmsPerRdo*kFstNumChanPerArm; // 6144: 128 channels * 8 APVs * 6 modules
290  int rdoIdx = elecId/ChPerRdo + 1; // 1-6
291  int armIdx = (elecId - (rdoIdx-1)*ChPerRdo)/kFstNumChanPerArm; // 0-2
292  int refElecChanId = elecId - (rdoIdx-1)*ChPerRdo - armIdx*kFstNumChanPerArm; // 0-2047
293  int refApvIdx = refElecChanId/kFstNumApvChannels; // 0-15
294  int portIdx = refApvIdx/8; // 0-1
295  int lclApvIdx = refApvIdx-portIdx*8; // 0-7
296 
297  if(elecId >= kFstNumElecIds || geoId >= kFstNumElecIds || apvId >= kFstNumApvs)
298  continue;
299 
300  for( unsigned char timeBin = 0; timeBin < kFstNumTimeBins; ++timeBin ){
301  Int_t adc = (*rawHitIter)->getCharge( timeBin );
302  if(adc) {
303  int t = (int)timeBin;
304  if(mTimeBinMask==0) t=0;
305  int code = kFstNumTimeBins * elecId + t;
306 
307  //exclude signal-like channels time bin by time bin in current event
308  bool pass = kTRUE;
309  if(mDoPedCut) {
310  // pass = (adc > (mPedVec1stLoop[code]-mPedCut*mRmsVec1stLoop[code])) && (adc < (mPedVec1stLoop[code]+mPedCut*mRmsVec1stLoop[code]));
311  pass = (adc < (mPedVec1stLoop[code]+mPedCut*mRmsVec1stLoop[code]));
312  hist_adcSpectrum[t]->Fill(elecId, adc);
313  }
314  if(!pass) continue;
315 
316  mMathPedVec[code] += adc;
317  mMathRmsVec[code] += adc*adc;
318  mMathCouVec[code] ++;
319 
320  if(mRunHist) {
321  TH1S* histPed = mHistPedVec[ code ];
322  if( !histPed ){
323  ss.str("");
324  ss.clear();
325  ss << "hist_Pedestal_Ch" << code / kFstNumTimeBins << "_TB" << code % kFstNumTimeBins;
326  histPed = new TH1S( ss.str().data(), "", 512, 0, kFstMaxAdc );
327  mHistPedVec[ code ] = histPed;
328  }
329  histPed->Fill( adc );
330  }
331 
332  if(mDoPedCut) { // need the third run to be accurate
333  int rstrip = (geoId % (kFstNumInnerSensorsPerWedge * kFstNumStripsPerInnerSensor + kFstNumOuterSensorsPerWedge * kFstNumStripsPerOuterSensor))/kFstNumPhiSegPerWedge;
334  if(rstrip > 3) rstrip = rstrip-4;
335  sumAdcPerEvent[lclApvIdx][rstrip][t] += adc - mPedVec1stLoop[code];
336  channelCountsPerEvent[lclApvIdx][rstrip][t] ++;
337  }
338  }//adc cut
339  }//time bin loop
340  }//raw hits loop
341 
342  if(mDoPedCut) {
343  //common mode calculation per event
344  for(int i=0; i<kFstApvsPerWedge; i++) {
345  for(int j=0; j<kFstNumRStripsPerSensor; j++) {
346  for(int k=0; k<kFstNumTimeBins; k++) {
347  if( channelCountsPerEvent[i][j][k] ) {
348  cmNoisePerEvent[i][j][k] = (double)sumAdcPerEvent[i][j][k] / channelCountsPerEvent[i][j][k];
349 
350  int apvId = wedgeIdx*kFstApvsPerWedge + i;
351  int groupId = apvId*kFstNumRStripsPerSensor + j;
352  int code = groupId*kFstNumTimeBins + k;
353 
354  TH1F* histCmn = mHistCmnVec[ code ];
355 
356  if( !histCmn ){
357  sc.str("");
358  sc.clear();
359  sc << "hist_CMNoise_APV" << apvId << "_R" << j << "_TB" << k;
360  histCmn = new TH1F( sc.str().data(), "", -1024, -kFstMaxAdc, kFstMaxAdc);
361  mHistCmnVec[ code ] = histCmn;
362  }
363  histCmn->Fill( cmNoisePerEvent[i][j][k] );
364  }
365  }
366  }
367  }
368 
369  for( rawHitIter = rawHitVec.begin(); rawHitIter != rawHitVec.end(); ++rawHitIter ){
370  int elecId = (*rawHitIter)->getChannelId();
371  int geoId = (*rawHitIter)->getGeoId(); //channel geometry Id, counting from 0 to 36863
372  int apvId = elecId/kFstNumApvChannels; //APV chip geometry Id, counting from 0 to 287
373  int ChPerRdo = kFstNumArmsPerRdo*kFstNumChanPerArm; // 6144: 128 channels * 8 APVs * 6 modules
374  int rdoIdx = elecId/ChPerRdo + 1; // 1-6
375  int armIdx = (elecId - (rdoIdx-1)*ChPerRdo)/kFstNumChanPerArm; // 0-2
376  int refElecChanId = elecId - (rdoIdx-1)*ChPerRdo - armIdx*kFstNumChanPerArm; // 0-2047
377  int refApvIdx = refElecChanId/kFstNumApvChannels; // 0-15
378  int portIdx = refApvIdx/8; // 0-1
379  int lclApvIdx = refApvIdx-portIdx*8; // 0-7
380 
381  if(elecId >= kFstNumElecIds || geoId >= kFstNumElecIds || apvId >= kFstNumApvs)
382  continue;
383 
384  for( unsigned char timeBin = 0; timeBin < kFstNumTimeBins; ++timeBin ){
385  Int_t adc = (*rawHitIter)->getCharge( timeBin );
386  if(adc) {
387  int t = (int)timeBin;
388  if(mTimeBinMask==0) t=0;
389  int code = kFstNumTimeBins * elecId + t;
390 
391  //exclude signal-like channels time bin by time bin in current event
392  bool pass = (adc < (mPedVec1stLoop[code]+mPedCut*mRmsVec1stLoop[code]));
393  if(!pass) continue;
394 
395  int rstrip = (geoId % (kFstNumInnerSensorsPerWedge * kFstNumStripsPerInnerSensor + kFstNumOuterSensorsPerWedge * kFstNumStripsPerOuterSensor))/kFstNumPhiSegPerWedge;
396  if(rstrip > 3) rstrip = rstrip-4;
397 
398  mMathPedRanVec[code] += adc-mPedVec1stLoop[code]-cmNoisePerEvent[lclApvIdx][rstrip][t];
399  mMathRmsRanVec[code] += (adc-mPedVec1stLoop[code]-cmNoisePerEvent[lclApvIdx][rstrip][t])*(adc-mPedVec1stLoop[code]-cmNoisePerEvent[lclApvIdx][rstrip][t]);
400  mMathCouRanVec[code] ++;
401 
402  if(mRunHist) {
403  TH1F* histRan = mHistRanVec[ code ];
404  if( !histRan ){
405  ss.str("");
406  ss.clear();
407  ss << "hist_RanNoise_Ch" << code / kFstNumTimeBins << "_TB" << code % kFstNumTimeBins;
408  histRan = new TH1F( ss.str().data(), "", 1024, -kFstMaxAdc, kFstMaxAdc );
409  mHistRanVec[ code ] = histRan;
410  }
411  histRan->Fill( adc-mPedVec1stLoop[code]-cmNoisePerEvent[lclApvIdx][rstrip][t] );
412  }
413  }//adc cut
414  }//time bin loop
415  }//raw hits loop
416  }// calculate random noise
417 
418  //clear raw hit objects in vector
419  while (!rawHitVec.empty()) delete rawHitVec.back(), rawHitVec.pop_back();
420  }//raw hit collection cut
421  }//wedgeIdx loop
422  }//ierr cut
423  evtIdx++;
424 
425  return ierr;
426 };
427 
428 // save as needed
429 Int_t StFstCalibrationMaker::Finish()
430 {
431  Int_t ierr = kStOk;
432 
433  if( !mHasFinished ){
434  mHasFinished = 1;
435  cout << "StFstCalibrationMaker::Finish()" << endl;
436 
437  //calculate pedestal/RMS with mathematical method
438  for(int i=0; i<kFstNumTimeBins*kFstNumElecIds; i++) {
439  double mathPed, mathRms;
440  if(mMathCouVec[i] < 1) {
441  mathPed = 0.;
442  mathRms = 100.;
443  }
444  else {
445  mathPed = mMathPedVec[i]/mMathCouVec[i];
446  mathRms = mMathRmsVec[i]/mMathCouVec[i];
447  }
448 
449  double variance = mathRms - mathPed*mathPed;
450  if(variance > 0)
451  mathRms = sqrt( variance );
452  else
453  mathRms = 100.;
454 
455  mMathPedVec[i] = mathPed;
456  mMathRmsVec[i] = mathRms;
457  }
458 
459  // calculate random noise with mathematical method
460  for(int i=0; i<kFstNumTimeBins*kFstNumElecIds; i++) {
461  double mathPedRan, mathRmsRan;
462  if(mMathCouRanVec[i] < 1) {
463  mathPedRan = 0.;
464  mathRmsRan = 10000.;
465  }
466  else {
467  mathPedRan = mMathPedRanVec[i]/mMathCouRanVec[i];
468  mathRmsRan = mMathRmsRanVec[i]/mMathCouRanVec[i];
469  }
470 
471  double variance = mathRmsRan - mathPedRan*mathPedRan;
472  if(variance > 0)
473  mathRmsRan = sqrt( variance );
474  else
475  mathRmsRan = 100.;
476 
477  mMathRanVec[i] = mathRmsRan;
478  }
479 
480  //calculate pedestal/RMS with histogram method
481  if(mRunHist) {
482  std::vector< TH1S* >::iterator mHistPedVecIter;
483  int elecIdx = 0;
484 
485  for( mHistPedVecIter = mHistPedVec.begin(); mHistPedVecIter != mHistPedVec.end(); ++mHistPedVecIter ){
486  TH1S *histPed = *mHistPedVecIter;
487  elecIdx = std::distance( mHistPedVec.begin(), mHistPedVecIter );
488  short timebin = elecIdx % kFstNumTimeBins;
489  int chanIdx = elecIdx / kFstNumTimeBins;
490  if(histPed) {
491  int entries = histPed->GetEntries();
492  float meanPed = histPed->GetMean();
493  float rmsPed = histPed->GetRMS();
494 
495  if( entries == 0 ) {
496  meanPed = 0.;
497  rmsPed = 100.; //marked as dead channel
498  }
499 
500  pedNoiseData_t &data = mPedVec[ elecIdx ];
501  data.n = entries;
502  data.ped = meanPed;
503  data.rms = rmsPed;
504 
505  if(mDoPedCut) {//only fill histograms in the 2nd loop
506  int geomIdx = mMappingVec[chanIdx];
507  hist_meanPed[timebin]->SetBinContent(geomIdx+1, meanPed);
508  hist_meanPed[timebin]->SetBinError(geomIdx+1, 0.);
509  hist_rmsPed[timebin]->SetBinContent(geomIdx+1, rmsPed);
510  hist_rmsPed[timebin]->SetBinError(geomIdx+1, 0.);
511  hist_sumPed[timebin]->Fill( meanPed );
512  hist_sumRms[timebin]->Fill( rmsPed );
513  }
514  }
515  }
516  }
517 
518  //calculate random noise with histogram method
519  if(mRunHist) {
520  std::vector< TH1F* >::iterator mHistRanVecIter;
521  int elecIdx = 0;
522 
523  for( mHistRanVecIter = mHistRanVec.begin(); mHistRanVecIter != mHistRanVec.end(); ++mHistRanVecIter ){
524  TH1F *histRan = *mHistRanVecIter;
525  elecIdx = std::distance( mHistRanVec.begin(), mHistRanVecIter );
526  short timebin = elecIdx % kFstNumTimeBins;
527  int chanIdx = elecIdx / kFstNumTimeBins;
528  if(histRan) {
529  int entries = histRan->GetEntries();
530  float rmsRan = histRan->GetRMS();
531 
532  if( entries == 0 ) {
533  rmsRan = 100.; //marked as dead channel
534  }
535 
536  pedNoiseData_t &data = mPedVec[ elecIdx ];
537  data.ran = rmsRan;
538 
539  if(mDoPedCut) {//only fill histograms in the 2nd loop
540  int geomIdx = mMappingVec[chanIdx];
541  hist_ranNoise[timebin]->SetBinContent(geomIdx+1, rmsRan);
542  hist_ranNoise[timebin]->SetBinError(geomIdx+1, 0.);
543  hist_sumRan[timebin]->Fill( rmsRan );
544  }
545  }
546  }
547  }
548 
549  std::vector< TH1F* >::iterator mHistCmnVecIter;
550  int groupIdx = 0;
551  for( mHistCmnVecIter = mHistCmnVec.begin(); mHistCmnVecIter != mHistCmnVec.end(); ++mHistCmnVecIter ){
552  TH1F *histCmn = *mHistCmnVecIter;
553  groupIdx = std::distance( mHistCmnVec.begin(), mHistCmnVecIter );
554 
555  short timebin = groupIdx % kFstNumTimeBins;
556  int gouIdx = groupIdx / kFstNumTimeBins;
557 
558  if(histCmn) {
559  int entries = histCmn->GetEntries();
560  float cmNoise = histCmn->GetRMS();
561 
562  if(entries == 0) {
563  cmNoise = 100.; //marked as dead chip
564  }
565 
566  cmNoiseData_t &data = mCmnVec[ groupIdx ];
567  data.n = entries;
568  data.cmn = cmNoise;
569 
570  if(mDoPedCut) {
571  hist_cmNoise[timebin]->SetBinContent(gouIdx+1, cmNoise);
572  hist_cmNoise[timebin]->SetBinError(gouIdx+1, 0.);
573  hist_sumCmn[timebin]->Fill(cmNoise);
574  }
575  }
576  }
577 
578  ierr = saveToFile();
579  }
580 
581  return ierr;
582 };
583 
584 // functions that actually do the saving
585 Int_t StFstCalibrationMaker::saveToFile()
586 {
587  Int_t ierr = kStOk;
588 
589  if(mDoOutput)
590  {
591  //create output file
592  StIOMaker *ioMaker = (StIOMaker * )GetMaker("inputStream");
593  if (!ioMaker) {
594  LOG_WARN << "StFstCalibrationMaker::Init(): No StIOMaker" << endm;
595  }
596 
597  TString filename = TString(ioMaker->GetFile());
598  int found = filename.Last('/');
599  if(found >= 0){
600  filename.Replace(0, found + 1, "");
601  }
602  found = filename.First(".");
603  if(found == 0) found = filename.Length();
604 
605  TString mRootFilename = filename.Data();
606  TString mPedNoiseFilename_hist = filename.Data();
607  TString mCmNoiseFilename = filename.Data();
608  TString mPedNoiseFilename_math = filename.Data();
609 
610  mRootFilename.Replace(found, mRootFilename.Length() - found, ".fstCaliQa.root");
611  LOG_INFO << "FST Calibration QA File Name: " << mRootFilename << endm;
612 
613  //create QA file
614  TFile *myRootFile = new TFile(mRootFilename.Data(),"RECREATE");
615  if( !myRootFile ) {
616  LOG_WARN << "Error recreating file '" << mRootFilename << "'" << endl;
617  ierr = kStWarn;
618  }
619  //save calibration QA histograms
620  for(unsigned char iTB=0; iTB<kFstNumTimeBins; iTB++) {
621  myRootFile->WriteTObject(hist_meanPed[iTB]);
622  myRootFile->WriteTObject(hist_rmsPed[iTB]);
623  myRootFile->WriteTObject(hist_cmNoise[iTB]);
624  myRootFile->WriteTObject(hist_ranNoise[iTB]);
625  myRootFile->WriteTObject(hist_sumPed[iTB]);
626  myRootFile->WriteTObject(hist_sumRms[iTB]);
627  myRootFile->WriteTObject(hist_sumCmn[iTB]);
628  myRootFile->WriteTObject(hist_sumRan[iTB]);
629  myRootFile->WriteTObject(hist_adcSpectrum[iTB]);
630  }
631  myRootFile->Close();
632 
633  //save pedestal and rms noise with mathematical method
634  mPedNoiseFilename_math.Replace(found, mPedNoiseFilename_math.Length() - found, ".fstPedNoise_math.dat");
635  LOG_INFO << "FST Pedestal and RMS File Name using math method: " << mPedNoiseFilename_math << endm;
636 
637  std::ofstream fout_ped_math( mPedNoiseFilename_math.Data(), std::ios_base::out & std::ios_base::trunc );
638  if( !fout_ped_math ){
639  LOG_ERROR << "Error opening file '" << mPedNoiseFilename_math << "'" << endm;
640  ierr = kStFatal;
641  }
642  fout_ped_math.setf(std::ios::fixed, std::ios::floatfield);
643  fout_ped_math.precision(5);
644 
645  for(int i=0; i<kFstNumTimeBins*kFstNumElecIds; i++) {
646  //obtain rdo/arm/apv/chan
647  int rdo = 0, arm = -1, apv = -1, chan = -1;
648  short timebin = i % kFstNumTimeBins;
649  int elecId = i / kFstNumTimeBins;
650  rdo = 1 + elecId/(kFstNumArmsPerRdo*kFstNumChanPerArm);
651  arm = (elecId%(kFstNumArmsPerRdo*kFstNumChanPerArm))/kFstNumChanPerArm;
652  apv = ((elecId%(kFstNumArmsPerRdo*kFstNumChanPerArm))%kFstNumChanPerArm)/kFstNumApvChannels;
653  chan = ((elecId%(kFstNumArmsPerRdo*kFstNumChanPerArm))%kFstNumChanPerArm)%kFstNumApvChannels;
654 
655  fout_ped_math << elecId << ' ' << rdo << ' ' << arm << ' ' << apv << ' ' << chan << ' ' << timebin << ' ' << mMathPedVec[i] << ' ' << mMathRmsVec[i] << ' ' << mMathRanVec[i] << endl;
656  }
657  fout_ped_math.close();
658 
659  //create a link as fstPedNoiseTable.dat in 1st loop
660  if(!mDoPedCut) {
661  char cmd[128];
662  sprintf(cmd, "/bin/ln -f -s %s fstPedNoiseTable.dat", mPedNoiseFilename_math.Data());
663  system(cmd);
664  }
665 
666  //save pedestal and rms noise with histogram method
667  if(mRunHist) {
668  mPedNoiseFilename_hist.Replace(found, mPedNoiseFilename_hist.Length() - found, ".fstPedNoise_hist.dat");
669  LOG_INFO << "FST Pedestal and RMS File Name using histogram method: " << mPedNoiseFilename_hist << endm;
670 
671  std::ofstream fout_ped_hist( mPedNoiseFilename_hist.Data(), std::ios_base::out & std::ios_base::trunc );
672  if( !fout_ped_hist ){
673  LOG_ERROR << "Error opening file '" << mPedNoiseFilename_hist << "'" << endm;
674  ierr = kStFatal;
675  }
676  fout_ped_hist.setf(std::ios::fixed, std::ios::floatfield);
677  fout_ped_hist.precision(5);
678 
679  pedNoiseDataVec_t::iterator pedDataVecIter;
680  int idx = 0;
681  for( pedDataVecIter = mPedVec.begin(); pedDataVecIter != mPedVec.end() && !ierr; ++pedDataVecIter, ++idx ){
682  //obtain rdo/arm/apv/chan
683  int rdo = 0, arm = -1, apv = -1, chan = -1;
684  short timebin = idx % kFstNumTimeBins;
685  int elecId = idx / kFstNumTimeBins;
686  rdo = 1 + elecId/(kFstNumArmsPerRdo*kFstNumChanPerArm);
687  arm = (elecId%(kFstNumArmsPerRdo*kFstNumChanPerArm))/kFstNumChanPerArm;
688  apv = ((elecId%(kFstNumArmsPerRdo*kFstNumChanPerArm))%kFstNumChanPerArm)/kFstNumApvChannels;
689  chan = ((elecId%(kFstNumArmsPerRdo*kFstNumChanPerArm))%kFstNumChanPerArm)%kFstNumApvChannels;
690 
691  fout_ped_hist << elecId << ' ' << rdo << ' ' << arm << ' ' << apv << ' ' << chan << ' ' << timebin << ' ' << pedDataVecIter->ped << ' ' << pedDataVecIter->rms << ' ' << pedDataVecIter->ran << endl;
692  }
693  fout_ped_hist.close();
694  }
695 
696  //save common mode noise
697  mCmNoiseFilename.Replace(found, mCmNoiseFilename.Length() - found, ".fstCmNoise.dat");
698  LOG_INFO << "FST Common Mode Noise File Name: " << mCmNoiseFilename << endm;
699 
700  std::ofstream fout_cmn( mCmNoiseFilename.Data(), std::ios_base::out & std::ios_base::trunc );
701  if( !fout_cmn ){
702  LOG_ERROR << "Error opening file '" << mCmNoiseFilename << "'" << endm;
703  ierr = kStFatal;
704  }
705  fout_cmn.setf(std::ios::fixed, std::ios::floatfield);
706  fout_cmn.precision(5);
707 
708  cmNoiseDataVec_t::iterator cmnDataVecIter;
709  int idx = 0;
710  for( cmnDataVecIter = mCmnVec.begin(); cmnDataVecIter != mCmnVec.end() && !ierr; ++cmnDataVecIter, ++idx ){
711  short timebin = idx % kFstNumTimeBins;
712  int groupId = idx / kFstNumTimeBins;
713  int rId = groupId % kFstNumRStripsPerSensor;
714  int apvId = groupId / kFstNumRStripsPerSensor;
715  int wedgeGeomId = 1 + apvId/kFstApvsPerWedge;
716  int rdo = (wedgeGeomId-1)/kFstNumWedsPerRdo + 1; //1-6
717  int arm = ((wedgeGeomId-1)%kFstNumWedsPerRdo)/kFstNumWedsPerArm; //0-2
718  // int apv = apvId%(kFstNumArmsPerRdo*kFstNumRdos); //0-15
719  int apv = apvId%kFstNumApvsPerArm; //0-15
720 
721  fout_cmn << apvId << ' ' << rdo << ' ' << arm << ' '<< apv << ' ' << rId << ' ' << timebin << ' ' << cmnDataVecIter->cmn << endl;
722  }
723  fout_cmn.close();
724  }
725 
726  return ierr;
727 };
728 
729 ClassImp( StFstCalibrationMaker );
StFstCalibrationMaker::cmNoiseData_t
Definition: StFstCalibrationMaker.h:56
StFstCalibrationMaker::pedNoiseData_t
Definition: StFstCalibrationMaker.h:45
kStWarn
Definition: Stypes.h:42
TObjectSet::GetObject
virtual TObject * GetObject() const
The depricated method (left here for the sake of the backward compatibility)
Definition: TObjectSet.h:56
kStErr
Definition: Stypes.h:44
kStOk
Definition: Stypes.h:41
Generated by   doxygen 1.8.5