StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
L2btowCalAlgo12.cxx
1 #include <stdio.h>
2 #include <string.h>
3 #include <stdlib.h>
4 #include <time.h>
5 #include <math.h>
6 #include "fakeRtsLog.h"
7 /*********************************************************
8  $Id: L2btowCalAlgo12.cxx,v 1.6 2012/03/21 18:18:03 jml Exp $
9  \author Jan Balewski, MIT, 2009
10  *****************************************************
11  Descripion:
12  calibrates Barrel towers, result is used by other L2-algo
13  *****************************************************/
14 
15 
16 #ifdef IS_REAL_L2 //in l2-ana environment
17  #include "../L2algoUtil/L2EmcDb2012.h"
18  #include "../L2algoUtil/L2Histo.h"
19 #else
20  #include "L2EmcDb2012.h"
21  #include "L2Histo.h"
22  #include "L2EmcGeom2012.h"
23 #endif
24 
25 #include "L2btowCalAlgo12.h"
26 
27 /************************
28 the _only_ copy of L2event for all L2-algos
29 **************************/
30 L2eventStream2012 globL2eventStream2012;
31 /*************************/
32 
33 //=================================================
34 //=================================================
35 L2btowCalAlgo12::L2btowCalAlgo12(const char* name, const char *uid, L2EmcDb2012* db, L2EmcGeom2012 *geoX, char* outDir, int resOff) : L2VirtualAlgo2012( name, uid, db, outDir, true, false, resOff) {
36  /* called one per days
37  all memory allocation must be done here
38  */
39 
40  geom=geoX;
41  if (!geom)
42  criticalError("L2btowCalAlgo is broken -- can't find geom.");
43 
44  setMaxHist(32);
45  createHisto();
46 
47  // initialize BTOW-Calibrated-data
48  int k;
49  for(k=0;k<L2eventStream2012::mxToken;k++){
50  L2BtowCalibData12 & btowCalibData=globL2eventStream2012.btow[k];
51  btowCalibData.nInputBlock=0;
52  btowCalibData.hitSize=0;
53  }
54 
55  LOG("JEFF", "sizeof L2 event stream: %d", sizeof(L2eventStream2012));
56 
57  }
58 
59 /* ========================================
60  ======================================== */
61 int
62 L2btowCalAlgo12::initRunUser( int runNo, int *rc_ints, float *rc_floats) {
63 
64  par_adcMask=(unsigned short)(-0x10);
65  par_pedOff=0x10/2; //hex must match with above
66 
67  // unpack params from run control GUI
68  par_dbg = rc_ints[0];
69  par_gainType = rc_ints[1];
70  par_nSigPed = rc_ints[2];
71 
72  par_twEneThres = rc_floats[0];
73  par_hotEtThres = rc_floats[1];;
74 
75  // verify consistency of input params
76  int kBad=0;
77  kBad+=0x00001 * (par_gainType<kGainZero || par_gainType>kGainOffline);
78  kBad+=0x00002 * (par_nSigPed<2 || par_nSigPed>5);
79  kBad+=0x00004 * (par_twEneThres<0.1 || par_twEneThres>1.5);
80 
81  if (mLogFile) {
82  fprintf(mLogFile,"L2%s algorithm initRun(R=%d), compiled: %s , %s\n params:\n",getName(),mRunNumber,__DATE__,__TIME__);
83  fprintf(mLogFile," - use BTOW=%d, gain Ideal=%d or Offline=%d, debug=%d\n",
84  par_gainType>=kGainIdeal, par_gainType==kGainIdeal, par_gainType==kGainOffline, par_dbg);
85  fprintf(mLogFile," - thresholds: ADC-ped> %d*sigPed .AND. energy>%.2f GeV \n", par_nSigPed, par_twEneThres);
86 
87  fprintf(mLogFile," - hot tower thresholds: ET/GeV=%.2f\n",par_hotEtThres);
88  fprintf(mLogFile,"initRun() params checked for consistency, Error flag=0x%04x\n",kBad);
89  }
90 
91  // initialize BTOW-Calibrated-data
92  int k;
93  for(k=0;k<L2eventStream2012::mxToken;k++){
94  L2BtowCalibData12 & btowCalibData=globL2eventStream2012.btow[k];
95  btowCalibData.nInputBlock=0;
96  btowCalibData.hitSize=0;
97  }
98 
99 
100  if(kBad) return -kBad;
101 
102  // clear content of all histograms
103  int i;
104  for (i=0; i<mxHA;i++) if(hA[i])hA[i]->reset();
105 
106  // update title of histos
107  char txt[1000];
108  sprintf(txt,"BTOW tower, E_T>%.2f GeV (input); x: BTOW RDO index=chan*30+fiber; y: counts",par_hotEtThres);
109  hA[10]->setTitle(txt);
110 
111  sprintf(txt,"BTOW tower, Et>%.2f GeV (input); x: BTOW softID",par_hotEtThres);
112  hA[11]->setTitle(txt);
113  sprintf(txt,"BTOW tower, Et>%.2f GeV (input); x: eta bin, [-1,+1]; y: phi bin ~ TPC sector",par_hotEtThres);
114  hA[12] ->setTitle(txt);
115 
116  sprintf(txt,"#BTOW towers / event , Et>%.2f GeV; x: # BTOW towers; y: counts",par_hotEtThres);
117  hA[14] ->setTitle(txt);
118 
119  // re-caluclate geometry properties
120  geom->btow.clear();
121  int nB=0; /* counts # of unmasekd towers */
122  int nBg=0; /* counts # of reasonable calibrated towers */
123  int nEneThr=0, nPedThr=0; //BTOW count # of towers above & below threshold
124  if(par_gainType>=kGainIdeal) // this disables the whole loop below
125  for(i=0; i<EmcDbIndexMax; i++) {
126  const L2EmcDb2012::EmcCDbItem *x=mDb->getByIndex(i);
127  if(mDb->isEmpty(x)) continue; /* dropped not mapped channels */
128  /*....... B A R R E L .................*/
129  if (!mDb->isBTOW(x) ) continue; /* drop if not BTOW */
130  if(x->fail) continue; /* dropped masked channels */
131  if(x->gain<=0) continue; /* dropped uncalibrated towers , tmp */
132  nB++;
133 
134  float adcThres=x->ped+par_nSigPed* fabs(x->sigPed);
135  float otherThr=x->ped+par_twEneThres*x->gain;
136 
137  if(adcThres<otherThr) { //use energy threshold if higher
138  adcThres=otherThr;
139  nEneThr++;
140  } else {
141  nPedThr++;
142  }
143 
144  /* use rdo index to match RDO order in the ADC data banks */
145  if(x->eta<=0 || x->eta>BtowGeom::mxEtaBin) return -90;
146  int ietaTw= (x->eta-1); /* correct */
147 
148  // use ideal gains for now, hardcoded
149 
150  if (par_gainType!=kGainIdeal) return -102;
151  geom->btow.gain2Ene_rdo[x->rdo]=geom->btow.idealGain2Ene[ietaTw];
152  geom->btow.gain2ET_rdo[x->rdo]=geom->getIdealAdc2ET();
153 
154  geom->btow.thr_rdo[x->rdo]=(int) (adcThres);
155  geom->btow.ped_rdo[x->rdo]=(int) (x->ped);
156  geom->btow.ped_shifted_rdo[x->rdo]=(unsigned short)(par_pedOff - x->ped);
157  nBg++;
158  }
159 
160  if (mLogFile) {
161  fprintf(mLogFile," found towers working=%d calibrated=%d, based on ASCII DB\n",nB,nBg);
162  fprintf(mLogFile," thresh defined by energy=%d or NsigPed=%d \n",nEneThr, nPedThr);
163  }
164 
165  return 0; //OK
166 
167 
168 }
169 
170 /* ========================================
171  ======================================== */
172 void
173 L2btowCalAlgo12::calibrateBtow(int token, int bemcIn, unsigned short *rawAdc){
174  // Btow calibration is a special case, must have one exit at the end
175  computeStart();
176 
177  //LOG("JEFF", "Token=%d",token);
178  token&=L2eventStream2012::tokenMask; // only protect against bad token, Gerard's trick
179  //LOG("JEFF", "Token=%d",token);
180 
181  //...... token is valid no w ........
182  L2BtowCalibData12 & btowCalibData=globL2eventStream2012.btow[token];
183  btowCalibData.nInputBlock++;
184 
185  // clear data for this token from previous event
186  btowCalibData.hitSize=0;
187 
188  int nTower=0; /* counts mapped & used ADC channels */
189  int nHotTower=0;
190  if(bemcIn && par_gainType>kGainZero) { // EVEVEVEVEVE
191  // ............process this event ...............
192  short rdo;
193  int low_noise_adc;
194  int adc; // pedestal subtracted
195  float et;
196  float low_noise_et; // bit-chopped
197  unsigned short *thr=geom->btow.thr_rdo;
198  unsigned short *ped=geom->btow.ped_rdo;
199  unsigned short *ped_shifted=geom->btow.ped_shifted_rdo;
200  float *gain2ET=geom->btow.gain2ET_rdo;
201  float *gain2Ene=geom->btow.gain2Ene_rdo;
202  HitTower1 *hit=btowCalibData.hit;
203  for(rdo=0; rdo<BtowGeom::mxRdo; rdo++){
204  if(rawAdc[rdo]<thr[rdo])continue;
205 
206  //LOG("JEFF", "rawAdc[%d] = %d thr=%d ped=%d",rdo,rawAdc[rdo],thr[rdo],ped[rdo],BtowGeom::mxRdo);
207 
208  if(nTower>=L2BtowCalibData12::mxListSize) break; // overflow protection
209  adc=rawAdc[rdo]-ped[rdo]; //do NOT correct for common pedestal noise
210  et=adc/gain2ET[rdo];
211  low_noise_adc=(rawAdc[rdo]+ped_shifted[rdo]) & par_adcMask;//note that ped_shifted is defined as pedOffset-ped.
212  low_noise_et=low_noise_adc/gain2ET[rdo];
213  hit->rdo=rdo;
214  hit->adc=adc;
215  hit->et=et;
216  hit->low_noise_et=low_noise_et;
217  hit->ene=adc/gain2Ene[rdo];
218  hit++;
219  nTower++;
220  // only monitoring
221  if(et >par_hotEtThres) {
222  hA[10]->fill(rdo);
223  nHotTower++;
224  }
225  }
226  btowCalibData.hitSize=nTower;
227 
228  // QA histos
229  hA[13]->fill(nTower);
230  hA[14]->fill(nHotTower);
231  if(nTower>=L2BtowCalibData12::mxListSize) mhN->fill(5); // was overflow
232 
233  } // EVEVEVEVEVE
234 
235  LOG(DBG, "hits = %d hot=%d",nTower,nHotTower);
236 
237  // debugging should be off for any time critical computation
238  if(par_dbg>0){
239  LOG(DBG,"L2-%s-compute: set adcL size=%d, get=%d\n",getName(),nTower,999); //tmp
240  LOG(DBG,"dbg=%s: found nTw=%d\n",getName(),nTower);
241  if(par_dbg>0) print0();
242  printCalibratedData(token);
243  }
244 
245  computeStop(token);
246 
247 }
248 
249 /* ========================================
250  ======================================== */
251 
252 void
253 L2btowCalAlgo12::clear(int token){
254  token&=L2eventStream2012::tokenMask; // only protect against bad token, Gerard's trick
255  //...... token is valid no w ........
256  L2BtowCalibData12 & btowCalibData=globL2eventStream2012.btow[token];
257  btowCalibData.hitSize=0;
258  return;
259 }
260 
261 /* ========================================
262  ======================================== */
263 void
264 L2btowCalAlgo12::computeUser(int token ){
265 
266  LOG(CRIT,"L2btowCalAlgo12::computeUser has been called and should not have been. Serious problem in L2");
267 }
268 
269 
270 /* ========================================
271  ======================================== */
272 void
273 L2btowCalAlgo12::finishRunUser() { /* called once at the end of the run */
274  // do whatever you want, log-file & histo-file are still open
275 
276  //report mean and RMS of #towers/event:
277  int mean=0, RMS=0;
278  hA[13]->findMean(&mean,&RMS);
279  if (mLogFile){
280  fprintf(mLogFile,"#BTOW_nonzero_towers_per_event: mean=%d, rms=%d\n",mean,RMS);
281  }
282  // search for hot tower, re-project histos vs. other representations
283  int bHotSum=1,bHotId=-1;
284  const int *data20=hA[10]->getData();
285  const L2EmcDb2012::EmcCDbItem *xB=mDb->getByIndex(402); // some wired default?
286 
287  int i;
288  for(i=0; i<EmcDbIndexMax; i++) {
289  const L2EmcDb2012::EmcCDbItem *x=mDb->getByIndex(i);
290  if(mDb->isEmpty(x)) continue;
291  if (!mDb->isBTOW(x) ) continue;
292  int softId=atoi(x->tube+2);
293  int ieta= (x->eta-1);
294  int iphi= (x->sec-1)*BtowGeom::mxSubs + x->sub-'a' ;
295  hA[11]->fillW(softId,data20[x->rdo]);
296  hA[12]->fillW(ieta, iphi,data20[x->rdo]);
297  if(bHotSum<data20[x->rdo]) {
298  bHotSum=data20[x->rdo];
299  bHotId=softId;
300  xB=x;
301  }
302  }
303 
304  if (mLogFile){
305  fprintf(mLogFile,"#BTOW_hot tower _candidate_ (bHotSum=%d of %d eve) :, softID %d , crate %d , chan %d , name %s\n",bHotSum,mEventsInRun,bHotId,xB->crate,xB->chan,xB->name);
306  }
307 
308  //...... QA tokens .....
309  int tkn1=99999, tkn2=0; // min/max token
310  int nTkn=0;
311  int tkn3=-1, nTkn3=-1; // most often used token
312 
313  int k;
314  for(k=0;k<L2eventStream2012::mxToken;k++){
315  L2BtowCalibData12 & btowCalibData=globL2eventStream2012.btow[k];
316  if(btowCalibData.nInputBlock==0) continue;
317  hA[1]->fillW(k,btowCalibData.nInputBlock);
318  if(nTkn3<btowCalibData.nInputBlock){
319  nTkn3=btowCalibData.nInputBlock;
320  tkn3=k;
321  }
322 
323  nTkn++;
324  if(tkn1>k) tkn1=k;
325  if(tkn2<k) tkn2=k;
326  }
327  if (mLogFile){
328  fprintf(mLogFile,"#BTOW_token_QA: _candidate_ hot token=%d used %d for %d events, token range [%d, %d], used %d tokens\n",tkn3,nTkn3,mEventsInRun,tkn1,tkn2,nTkn);
329  }
330 
331 }
332 
333 
334 //=======================================
335 //=======================================
336 void
337 L2btowCalAlgo12::createHisto() {
338  memset(hA,0,mxHA*sizeof(L2Histo*));
339  //token related spectra
340  hA[1]=new L2Histo(1,"L2-btow-calib: seen tokens; x: token value; y: events ",L2eventStream2012::mxToken);
341 
342  // BTOW raw spectra
343  hA[10]=new L2Histo(10,"btow hot tower 1", BtowGeom::mxRdo); // title upadted in initRun
344  hA[11]=new L2Histo(11,"btow hot tower 2", BtowGeom::mxRdo); // title upadted in initRun
345  hA[12]=new L2Histo(12,"btow hot tower 3", BtowGeom::mxEtaBin,BtowGeom::mxPhiBin); // title upadted in initRun
346  hA[13]=new L2Histo(13,"BTOW #tower w/ energy /event; x: # BTOW towers; y: counts", 200);
347  hA[14]=new L2Histo(14,"# hot towers/event", 100);
348 
349 }
350 
351 
352 
353 /* ========================================
354  ======================================== */
355 void
356 L2btowCalAlgo12::print0(){ // full raw input ADC array
357  // empty
358  }
359 
360 
361 /****************************************************
362  $Log: L2btowCalAlgo12.cxx,v $
363  Revision 1.6 2012/03/21 18:18:03 jml
364  got rid of printfs from 2012 files
365 
366  Revision 1.5 2011/10/19 16:12:10 jml
367  more 2012 stuff
368 
369  Revision 1.4 2011/10/19 15:58:06 jml
370  more compile offline
371 
372  Revision 1.3 2011/10/19 15:39:42 jml
373  2012
374 
375  Revision 1.2 2011/10/19 14:34:23 jml
376  added fakeRtsLog.h to turn log statements into printfs
377 
378  Revision 1.1 2011/10/18 15:11:41 jml
379  adding 2012 algorithms
380 
381  Revision 1.3 2010/04/18 02:53:35 pibero
382  Fixed memset bug.
383 
384  Revision 1.2 2010/04/17 16:42:09 pibero
385  *** empty log message ***
386 
387  Revision 1.5 2008/02/01 00:16:40 balewski
388  add mxListSize to BTOW/ETOW calibration
389 
390  Revision 1.4 2008/01/30 00:47:15 balewski
391  Added L2-Etow-calib
392 
393  Revision 1.3 2008/01/18 23:29:13 balewski
394  now L2result is exported
395 
396  Revision 1.2 2008/01/16 23:32:34 balewski
397  toward token dependent compute()
398 
399  Revision 1.1 2007/12/19 02:30:18 balewski
400  new L2-btow-calib-2008
401 
402  Revision 1.1 2007/11/19 22:18:25 balewski
403  most L2algos provide triggerID's
404 
405 
406 */
407 
408 
Generated by   doxygen 1.8.5