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