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