StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
L2pedAlgo.cxx
1 #include <stdio.h>
2 #include <assert.h>
3 #include <string.h>
4 #include <stdlib.h>
5 
6 /*********************************************************************
7  * $Id: L2pedAlgo.cxx,v 1.11 2009/11/19 15:48:46 balewski Exp $
8  * \author Jan Balewski, IUCF, 2006
9  *********************************************************************
10  * Descripion:
11  * pedestal algo in L2 , for BTOW & ETOW
12  *********************************************************************
13  */
14 
15 
16 #ifdef IS_REAL_L2 //in l2-ana environment
17  #include "../L2algoUtil/L2EmcDb.h"
18  #include "../L2algoUtil/L2Histo.h"
19 #else
20  #include "StTriggerUtilities/L2Emulator/L2algoUtil/L2EmcDb.h"
21  #include "StTriggerUtilities/L2Emulator/L2algoUtil/L2Histo.h"
22 #endif
23 
24 #include "L2pedAlgo.h"
25 #include "L2pedResults2006.h"
26 
27 //=================================================
28 //=================================================
29 L2pedAlgo::L2pedAlgo(const char* name, L2EmcDb* db, char* outDir, int resOff)
30  : L2VirtualAlgo( name, db, outDir, resOff) {
31  /* called one per days
32  all memory allocation must be done here
33  */
34 
35  par_pedSubtr=false;
36  par_saveBinary=false;
37  par_dbg=0;
38  par_prescAccept=0;
39 
40  int i;
41  for(i=0;i<MaxBtowRdo;i++) {
42  char tit[100];
43  sprintf(tit,"BTOW ADC for rdo=%d; ADC+%d",i,-minAdc);
44  btowAdc[i]=new L2Histo(20000+i,tit,maxAdc-minAdc+1);
45  }
46 
47  for(i=0;i<MaxEtowRdo;i++) {
48  char tit[100];
49  sprintf(tit,"ETOW ADC for rdo=%d; ADC ",i);
50  etowAdc[i]=new L2Histo(10000+i,tit,maxAdc-minAdc+1);
51  }
52 
53  // aux histos
54  memset(hA,0,sizeof(hA));
55  hA[10]=new L2Histo(10, (char*)"total event counter; x=cases",6);
56  hA[11]=new L2Histo(11, (char*)"L2 time used per input event; x: time (CPU 20*kTics); y: events ",500);
57 
58  // BTOW raw spectra
59  hA[20]=new L2Histo(20, (char*)"BTOW pedRes Y=ADC-DBped ; x: chan + 160*crate", 4800);
60  hA[21]=new L2Histo(21, (char*)"BTOW pedRes Z=ADC-DBped, saturated @ |3|; x: etaBin ,[-1,+1]; y: phi bin ~sector",40,120);
61 
62  // ETOW raw spectra
63  hA[30]=new L2Histo(30, (char*)"ETOW pedRes Y=ADC-DBped ; x: chan + 128*crate", 768);
64  hA[31]=new L2Histo(31, (char*)"ETOW pedRes Z=ADC-DBped, saturated @ |3|; x: 12 - Endcap etaBin ,[+1,+2]; y: phi bin ~sector",12,60);
65 
66 
67  printf("L2pedAlgo instantiated, logPath='%s'\n",mOutDir);
68 
69 }
70 
71 /*========================================
72  ======================================== */
73 int
74 L2pedAlgo::initRun(int runNo, int *rc_ints, float *rc_floats) {
75  //myName is not used.
76  // update DB if run # has changed
77  //printf("aaa L2pedAlgo::initRun runNo=%d\n",runNo);
78  if(mDb->initRun(runNo)) return -27;
79  // DB must be initialized prior to lookup tables
80 
81  // unpack input params
82 
83 
84  par_pedSubtr =rc_ints[0]!=0;
85  par_speedFact =rc_ints[1];
86  par_saveBinary=rc_ints[2]!=0;
87  par_dbg =rc_ints[3];
88  par_prescAccept=rc_ints[4];
89 
90  if(par_prescAccept<0) par_prescAccept=0; // prescale can't be negative
91  //note speedFactor can be only powers of 2, range [1-256]
92  if(par_speedFact<1) par_speedFact=1;
93  if(par_speedFact>2) {// ASSURE ONLY POWERS OF 2, round down
94  if(par_speedFact<4) par_speedFact=2;
95  else if(par_speedFact<8) par_speedFact=4;
96  else if(par_speedFact<16) par_speedFact=8;
97  else if(par_speedFact<32) par_speedFact=16;
98  else if(par_speedFact<64) par_speedFact=32;
99  else if(par_speedFact<192) par_speedFact=64;
100  else par_speedFact=192;
101  }
102  // can't be bigger, since there is only 6 etow crates
103 
104  s_lastB=s_lastE=0;
105  s_stepB=MaxBtowRdo/ par_speedFact ;
106  s_stepE=MaxEtowRdo/ par_speedFact ;
107 
108 
109  /* .... clear content */
110  memset(db_btowPed, 0,sizeof(db_btowPed));
111  memset(db_etowPed, 0,sizeof(db_etowPed));
112 
113 
114  int i;
115  int nBtowOk=0, nEtowOk=0;
116  for(i=0; i<EmcDbIndexMax; i++) {
117  const L2EmcDb::EmcCDbItem *x=mDb->getByIndex(i);
118  if(mDb->isEmpty(x)) continue; // dropped not mapped channels
119  if (mDb->isBTOW(x) ) {
120  db_btowPed[x->rdo]=(int) (x->ped);
121  nBtowOk++;
122  } else if (mDb->isETOW(x) ) {
123  db_etowPed[x->rdo]=(int) (x->ped);
124  nEtowOk++;
125  }
126  }
127 
128  printf("L2ped algorithm init()... params:\n dbg=%d, pedSubtr=%d saveBinHist=%d speedFact=%d prescAccept=%d\n mapped channels: nBtow=%d nEtow=%d\n",par_dbg,par_pedSubtr,par_saveBinary,par_speedFact,par_prescAccept,nBtowOk,nEtowOk);
129 
130  for(i=0;i<MaxBtowRdo;i++) btowAdc[i]->reset();
131  for(i=0;i<MaxEtowRdo;i++) etowAdc[i]->reset();
132  for(i=0;i<mxHA;i++) if(hA[i]) hA[i]->reset();
133 
134  nInp=0;
135  run_number=runNo;
136  return 0; // OK
137 }
138 
139 
140 /*========================================
141  ======================================== */
142 bool
143 L2pedAlgo::doEvent(int L0trg, int inpEveId, TrgDataType* trgData,
144  int bemcIn, unsigned short *bemcData,
145  int eemcIn, unsigned short *eemcData){
146  // not used: L0trg, inpEveId
147  mAccept=true; // it never aborts
148  /* STRICT TIME BUDGET START ...., well a bit relaxed for this algo*/
149  rdtscl_macro(mEveTimeStart);
150  nInp++;
151  hA[10]->fill(0);
152  if(par_prescAccept>0) { // value=1 accepts 100% at maximal speed
153  if((rand()>>4) % par_prescAccept ) return false;
154  hA[10]->fill(5);
155  return true;
156  }
157 
158  myTrigData=trgData;
159 
160  /* *****************************************
161  the code below is NOT optimized for speed
162  one would need to make loop indexed by RDO to speed it up
163  not worth for pedestal calculation
164  This version takes 2,500 kTicks/eve --> ~1.5msec
165  ******************************************** */
166 
167  short rdo;
168  if( bemcIn ) {
169  short first=s_lastB%MaxBtowRdo;
170  s_lastB=first+s_stepB;
171  if(first==0) hA[10]->fill(1);
172  // printf("B: f=%d l=%d\n",first,s_lastB);
173 
174  for(rdo=first; rdo<s_lastB; rdo++){
175  int adc=bemcData[rdo];
176  if(par_pedSubtr) adc-=db_btowPed[rdo];
177  btowAdc[rdo]->fill(adc-minAdc);
178  // printf("B: rdo=%d raw=%d ped=%d val=%d\n",rdo,bemcData[rdo],db_etowPed[rdo],adc-minAdc);
179  }
180  }
181 
182 
183  if( eemcIn ) {
184  short first=s_lastE%MaxEtowRdo;
185  s_lastE=first+s_stepE;
186  if(first==0) hA[10]->fill(2);
187  //printf("E: f=%d l=%d\n",first,s_lastE);
188 
189  for(rdo=first; rdo<s_lastE; rdo++){
190  int adc=eemcData[rdo];
191  if(par_pedSubtr) adc-=db_etowPed[rdo];
192  etowAdc[rdo]->fill(adc-minAdc);
193  // printf("E: rdo=%d raw=%d ped=%d val=%d\n",rdo,eemcData[rdo],db_etowPed[rdo],adc-minAdc);
194  }
195  }
196 
197  //...this event will be accepted
198  L2pedResults2006 out; // all output bits lump together
199  memset(&out,0,sizeof(out));
200 
201  out.int0.decision=
202  ( (bemcIn>0) <<3 ) +
203  ( (eemcIn>0) <<4 ) +
204 
205  ( par_pedSubtr <<6 ) ;
206 
207  rdtscl_macro(mEveTimeStop);
208  mEveTimeDiff=mEveTimeStop-mEveTimeStart;
209  int kTick=mEveTimeDiff/1000;
210  hA[11]->fill(kTick/20);
211 
212  // printf("jkTick=%d\n",kTick);
213  const unsigned short maxKT=30000;
214  out.int0.kTick= kTick>maxKT ? maxKT : (int)kTick;
215 
216  unsigned int *outPlace=myTrigData->TrgSum.L2Result+ mResultOffset;
217  memcpy(outPlace,&out,sizeof( L2pedResults2006));
218 
219  if(par_dbg) L2pedResults2006_print(&out);
220 
221  return mAccept;
222 }
223 
224 /*========================================
225  ======================================== */
226 void
227 L2pedAlgo::finishRun() {/* called once at the end of the run */
228  if(run_number<=0) return; // algo not used for current run
229  int nBtowLow=0, nBtowHigh=0 ,nEtowLow=0, nEtowHigh=0 ;
230 
231  char fname[1000];
232  sprintf(fname,"%s/run%d.l2ped.out",mOutDir,run_number); // actual L2 output destination
233  printf("L2ped_finish('%s') , finding pedestals...\n",fname);
234 
235  FILE *fd=fopen(fname,"w");
236  if(fd==0) {printf("failed to open output %s file,skip ped_finish()\n",fname); return;}
237  fprintf(fd,"#L2-ped algorithm finishRun(%d), compiled: %s , %s\n",run_number,__DATE__,__TIME__);
238  fprintf(fd,"#params: pedSubtr=%d speedFact=%d saveBin=%d debug=%d prescAccept=%d\n",par_pedSubtr,par_speedFact,par_saveBinary,par_dbg,par_prescAccept);
239  hA[10]->printCSV(fd); // event accumulated
240  int iMax=-3, iFWHM=-4;
241  hA[11]->findMax( &iMax, &iFWHM);
242  fprintf(fd,"#L2ped CPU/eve MPV %d kTicks, FWHM=%d, seen eve=%d\n",iMax, iFWHM,nInp);
243  printf("L2ped CPU/eve MPV %d kTicks, FWHM=%d, seen eve=%d\n",iMax, iFWHM,nInp);
244  if(par_saveBinary) fprintf(fd,"#L2ped will save full spectra for all towers\n");
245 
246  int par_topAdc=100;
247  int maxPedDeviation=5;
248 //int iadcHigh=maxPedDeviation-minAdc;
249 // int iadcLow =maxPedDeviation+minAdc;
250 
251 // fixed sign on July 23, 2007, JanB
252  int iadcHigh=maxAdc - maxPedDeviation;
253  int iadcLow =minAdc + maxPedDeviation;
254 
255  char xAxis[100];
256  sprintf(xAxis,"raw ADC + %d",-minAdc);
257  if(par_pedSubtr) sprintf(xAxis,"ADC - ped + %d",-minAdc);
258 
259 
260  fprintf(fd,"# L2ped-Adc spectra, run=%d, Z-scale is ln(yield), only first digit shown; maxPedDev=%d table format:\n# name, ped, sigPed, crate, chan, softID-m-s-e, RDO_ID;\n# ADC spectrum: [%d ... <=-10 ... *=0 ... >=+10 ... :=+20 ... %d], Xaxis=%s\n",run_number,maxPedDeviation,minAdc,par_topAdc,xAxis);
261 
262  int i;
263  int nB=0;
264  for(i=0; i<EmcDbIndexMax; i++) {
265  const L2EmcDb::EmcCDbItem *x=mDb->getByIndex(i);
266  if(mDb->isEmpty(x)) continue; /* dropped not mapped channels */
267  if (mDb->isBTOW(x) ||mDb->isETOW(x) ) { //
268  // mDb->printItem(x);
269  nB++;
270  /* if(nB>30) return; */
271  int iMax=-3, iFWHM=-4;
272  char pedQA='?';
273  L2Histo *h=0;
274  if(mDb->isBTOW(x)) h= btowAdc[x->rdo];
275  else if(mDb->isETOW(x)) h= etowAdc[x->rdo];
276  else continue;
277 
278 
279  int pedRes=999;
280  int maxRes=3;
281  if(h->findMax( &iMax, &iFWHM)) {
282  pedQA='0';
283  if(iMax<iadcLow) pedQA='-';
284  else if(iMax>iadcHigh) pedQA='+';
285  pedRes=iMax+minAdc;
286  if(!par_pedSubtr) pedRes=int(pedRes - x->ped); // this looks funny but is right
287  }
288 
289  if(mDb->isBTOW(x)) { //BTOW ...............
290  if(pedQA=='-' ) nBtowLow++;
291  else if(pedQA=='+') nBtowHigh++;
292  //........ residual monito histos ....
293  int ieta= (x->eta-1);
294  int iphi= (x->sec-1)*10 + x->sub-'a' ;
295  int ihard=x->chan+(x->crate-1)*160;
296  if(x->fail) pedRes =0;
297  hA[20]->fillW(ihard,pedRes);
298  if(x->fail) {
299  pedRes =-100;
300  } else {
301  if(pedRes<-maxRes) pedRes=-maxRes;
302  if(pedRes>maxRes) pedRes=maxRes;
303  }
304  hA[21]->fillW(ieta, iphi,pedRes);
305  } else { // ETOW ...........
306  if(pedQA=='-' ) nEtowLow++;
307  else if(pedQA=='+') nEtowHigh++;
308  //........ residual monito histos ....
309  int ieta= 12-x->eta;
310  int iphi= (x->sec-1)*5 + x->sub-'A' ;
311  int ihard=x->chan+(x->crate-1)*128;
312  if(x->fail) pedRes =0;
313  hA[30]->fillW(ihard,pedRes);
314  if(x->fail) {
315  pedRes =-100;
316  } else {
317  if(pedRes<-maxRes) pedRes=-maxRes;
318  if(pedRes>maxRes) pedRes=maxRes;
319  }
320  hA[31]->fillW(ieta, iphi,pedRes);
321  }
322 
323 
324  char okC=' ';
325  if(x->fail) okC='#';
326  fprintf(fd,"%c%s %3d %4.1f 0x%02x 0x%02x %15s %4d ",okC,x->name,iMax+minAdc,iFWHM/2.3,x->crate, x->chan,x->tube, x->rdo);
327  h->printPed(fd,minAdc,par_topAdc,' ');
328  fprintf(fd,"qa=%c\n",pedQA);
329  // exit(1);
330  } /* end of BTOW & ETOW */
331  } /* end of pixels */
332 
333 
334  fprintf(fd,"#L2ped_finishRun() # of towers with |ped-pedDB| >10 chan\n# BTOW: nLow=%d nHigh=%d ; ETOW nLow=%d, nHigh=%d\n",nBtowLow, nBtowHigh,nEtowLow, nEtowHigh);
335  fprintf(fd,"# found peds for nB+E=%d , seen events=%d\n",nB,nInp);
336  printf("l2ped_finish() found peds for nB+E=%d , seen events=%d\n",nB,nInp);
337 
338  fclose(fd);
339 
340 
341  sprintf(fname,"%s/run%d.l2ped.hist.bin",mOutDir,run_number); // full spectra in binary form
342 
343  fd=fopen(fname,"w");
344  if(fd==0) {
345  printf("failed to open output %s file,skip ped_finish()\n",fname);
346  goto end;
347  }
348 
349  for(int j=0;j<mxHA;j++) { // save auxil histos
350  if(hA[j]==0) continue;
351  hA[j]->write(fd);
352  }
353 
354  if(par_saveBinary) {
355  printf("l2ped_finish('%s') , save FULL spectra binary...\n",fname);
356 
357  // .................. SAVE FULL SPECTRA BINARY ...........
358  for(i=0;i<mxHA;i++) if(hA[i]) hA[i]->write(fd);
359  for(i=0; i<EmcDbIndexMax; i++) {
360  const L2EmcDb::EmcCDbItem *x=mDb->getByIndex(i);
361  if(mDb->isEmpty(x)) continue; /* dropped not mapped channels */
362  L2Histo *h=0;
363  char tit[200];
364  if(mDb->isBTOW(x) ){
365  h= btowAdc[x->rdo];
366  //L2EmcDb::printItem(x);
367  sprintf(tit,"BTOW=%s cr/ch=%03d/%03d stat/0x=%04x+%04x soft=%s; %s",x->name, x->crate, x->chan, x->stat, x->fail,x->tube,xAxis);
368  //printf("tit=%s=\n",tit);
369  } else if(mDb->isETOW(x) ) {
370  h= etowAdc[x->rdo];
371  //L2EmcDb::printItem(x);
372  sprintf(tit,"ETOW=%s cr/ch=%03d/%03d stat/0x=%04x+%04x pname=%s; %s",x->name, x->crate, x->chan, x->stat, x->fail,x->tube,xAxis);
373  }
374  if(h==0) continue; //just in case
375  h->setTitle(tit);
376  h->write(fd); // change title, add cr/chan/name, pedSubtrFlag
377  // break;
378  }
379  printf("l2ped_finish() binary full spectra saved\n");
380  }
381  fclose(fd);
382  end:
383  run_number=-2;
384 }
385 
386 
387 
388 /**********************************************************************
389  $Log: L2pedAlgo.cxx,v $
390  Revision 1.11 2009/11/19 15:48:46 balewski
391  add (char*) to many strings to make SL5 happ, few other adjustments
392 
393  Revision 1.10 2007/11/19 22:18:31 balewski
394  most L2algos provide triggerID's
395 
396  Revision 1.9 2007/11/18 21:58:58 balewski
397  L2algos triggerId list fixed
398 
399  Revision 1.8 2007/11/13 23:06:09 balewski
400  toward more unified L2-algos
401 
402  Revision 1.7 2007/11/13 00:12:38 balewski
403  added offline triggerID, take1
404 
405  Revision 1.6 2007/11/08 04:02:33 balewski
406  run on l2ana as well
407 
408  Revision 1.5 2007/11/02 20:44:54 balewski
409  cleanup
410 
411  Revision 1.4 2007/11/02 17:43:11 balewski
412  cleanup & it started to work w/ L2upsilon
413 
414  Revision 1.3 2007/11/02 03:03:50 balewski
415  modified L2VirtualAlgo
416 
417  Revision 1.2 2007/10/25 02:07:06 balewski
418  added L2upsilon & binary event dump
419 
420  Revision 1.1 2007/10/11 00:33:24 balewski
421  L2algo added
422 
423  Revision 1.5 2006/03/28 19:46:51 balewski
424  ver16b, in l2new
425 
426  Revision 1.4 2006/03/11 17:08:35 balewski
427  now CVS comments should work
428 
429 */
430 
Generated by   doxygen 1.8.5