StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
L2wBemc2009.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: L2wBemc2009.cxx,v 1.2 2009/11/19 15:48:49 balewski Exp $
9  * \author Jan Balewski,MIT , 2009
10  *********************************************************************
11  * Descripion: see .h
12  *********************************************************************
13  */
14 // super smart macro from Pibero, needs only final output 'long long' as argument
15 #define rdtscll(val) do { \
16  unsigned int __a,__d; \
17  __asm__ __volatile__("rdtsc" : "=a" (__a), "=d" (__d)); \
18  (val) = ((unsigned long long)__a) | (((unsigned long long)__d)<<32); \
19 } while(0)
20 
21 const float stepETH=5;// for QA histo
22 
23 
24 #ifdef IS_REAL_L2 //in l2-ana environment
25  #include "../L2algoUtil/L2EmcDb.h"
26  #include "../L2algoUtil/L2Histo.h"
27 #else
28  #include "../L2algoUtil/L2EmcDb.h"
29  #include "../L2algoUtil/L2Histo.h"
30  #include "../L2algoUtil/L2EmcGeom.h"
31 #endif
32 
33 #include "L2wBemc2009.h"
34 
35 //=================================================
36 //=================================================
37 L2wBemc2009::L2wBemc2009(const char* name, L2EmcDb* db, L2EmcGeom *geoX, char* outDir, int resOff) : L2VirtualAlgo2009( name, db, outDir, true, false, resOff) {
38  /* called one per day
39  all memory allocation must be done here
40  */
41 
42  mGeom=geoX;
43  if (!mGeom)
44  criticalError((char*)"L2wBemc is broken -- can't find geom.");
45 
46  setMaxHist(32); // set upper range, I uses only 2^N -it is easier to remember
47  createHisto();
48 
49  //------- self-consistency checks, should never fail
50  if (sizeof(L2wResult2009)!= L2wResult2009::mySizeChar)
51  criticalError((char*)"L2wBemc has failed consistency check. sizeof(L2wResult2009)!= L2wResult2009::mySizeChar");
52 
53 }
54 
55 /* ========================================
56  ======================================== */
57 int
58 L2wBemc2009::initRunUser( int runNo, int *rc_ints, float *rc_floats) {
59 
60  // unpack params from run control GUI
61  par_dbg = rc_ints[0];
62  par_RndAcceptPrescale = rc_ints[1];
63  par_seedEtThres = rc_floats[0];
64  par_clustEtThres = rc_floats[1];
65 
66  // verify consistency of input params
67  int kBad=0;
68  kBad+=0x00002 * (par_clustEtThres<par_seedEtThres);
69 
70  // put notes about configuration into the log file
71  if (mLogFile) {
72  fprintf(mLogFile,"L2%s algorithm initRun(R=%d), compiled: %s , %s\n params:\n",
73  getName(),mRunNumber,__DATE__,__TIME__);
74  fprintf(mLogFile," - use seedThres=%.2f (GeV), debug=%d, prescale=%d (0=off,1=100prc, 2=50proc, etc)\n",
75  par_seedEtThres,par_dbg, par_RndAcceptPrescale);
76  fprintf(mLogFile," - accept event cluster Thres=%.2f (GeV)\n",
77  par_clustEtThres);
78  fprintf(mLogFile,"initRun() params checked for consistency, Error flag=0x%04x\n",
79  kBad);
80  }
81 
82  // clear content of all histograms & token-dependent memory
83  int i;
84  for (i=0; i<mxHA;i++) if(hA[i])hA[i]->reset();
85  memset(mBtow,0,sizeof(mBtow));
86 
87  if(kBad>0) return -1*kBad;
88  else if(kBad<0) return kBad;
89 
90 
91  // update titles of histos
92  char txt[1000];
93 
94  sprintf(txt,"W-BTOW-accepted: acc seed tower ET>%.2f GeV; BTOW softID",par_seedEtThres);
95  hA[9]->setTitle(txt);
96 
97  sprintf(txt,"W Btow accepted, seed tower Et>%.2f GeV; eta bin [-1,+1]; y: phi bin ~ TCP sector",par_seedEtThres);
98  hA[10]->setTitle(txt);
99 
100 
101  for ( int index=0; index<EmcDbIndexMax; index++ )
102  {
103  const L2EmcDb::EmcCDbItem *x = mDb->getByIndex(index);
104  if ( x==0 ) continue;
105  if ( !mDb->isBTOW(x) ) continue;
106  int sec = x->sec - 1;
107  int sub = 8192;
108  sub = x->sub - 'a';
109  int eta = x->eta - 1;
110  int phi = BtowGeom::mxSubs *sec + sub;
111  int tow = BtowGeom::mxEtaBin *phi + eta; // phi- changes faster
112  int rdo = x->rdo;
113  if (tow<0 || tow>mxBtow || rdo<0 || rdo>mxBtow) return -101;
114 
115  mTower2rdo[ tow ] = rdo; // returns rdo channel for given tower
116  mRdo2tower[ rdo ] = tow; // returns tower for given rdo channel
117  }
118  return 0; //OK
119 
120 }
121 
122 
123 
124 /* ========================================
125  ======================================== */
126 float
127 L2wBemc2009::sumET(int phi, int eta) {
128  int tow = BtowGeom::mxEtaBin *((phi+BtowGeom::mxPhiBin)%BtowGeom::mxPhiBin) + ((eta+BtowGeom::mxEtaBin)%BtowGeom::mxEtaBin); // phi- changes faster
129 
130  const int maxTowers = BtowGeom::mxEtaBin * BtowGeom::mxPhiBin;
131  int towPlusOne;
132  float sum;
133  sum=0;
134  sum=wrkBtow_et[tow];
135  towPlusOne = tow+1;
136  towPlusOne%= maxTowers;
137  sum+=wrkBtow_et[towPlusOne];
138 
139  //if(tow==0 || towPlusOne==0) {
140  //printf("tow : %d, %d --> %f %f \n",tow, towPlusOne,wrkBtow_et[tow],wrkBtow_et[towPlusOne]);
141  //}
142  tow+=BtowGeom::mxEtaBin;
143  tow%=maxTowers;
144 
145  sum+=wrkBtow_et[tow];
146 
147  towPlusOne = tow+1;
148  towPlusOne%= maxTowers;
149  sum+=wrkBtow_et[towPlusOne];
150 
151  //if(tow==0 || towPlusOne==0) {
152  //printf("tow : %d, %d --> %f %f \n",tow, towPlusOne,wrkBtow_et[tow],wrkBtow_et[towPlusOne]);
153  //}
154  // printf("B sumET=%f\n",sum);
155  return sum;
156 }
157 
158 
159 /* ========================================
160  ======================================== */
161 void
162 L2wBemc2009::computeUser(int token){
163 
164  clearEvent(token);
165 
166  // ------ PROJECT INPUT LIST TO 2D ARRAY AND SCAN FOR SEED TOWERS ----
167  int i;
168  // printf("L2-%s-compute: ---BTOW ADC list--- size=%d\n",getName(),*globEve_btow_hitSize);
169 
170  // access inpute list
171  const HitTower1 *hit=mEveStream_btow[token].get_hits();
172  const int hitSize=mEveStream_btow[token].get_hitSize();
173  // preapre output list
174  L2wBemcEvent2009 *btowEve=mBtow+token;
175 
176  for(i=0;i< hitSize;i++,hit++) {
177  int tower=mRdo2tower[hit->rdo];
178  wrkBtow_et[tower]=hit->et;
179  if(hit->et<par_seedEtThres)continue;
180  wrkBtow_tower_seed[wrkBtow_tower_seed_size++]=tower;
181  //printf("A seed TWID=%d \n",tower);
182  }
183  hA[2]->fill(hitSize);
184  int seedTow=-1,seedEta=-1,seedPhi=-1;
185  float clustET=0;
186  btowEve->isFresh=L2wBemcEvent2009::kDataFresh;
187 
188  // printf("B nseed=%d\n",wrkBtow_tower_seed_size);
189  // ----------- FIND 2x2 CLUSTER AROUND EVERY SEED -----
190  for(i=0; i<wrkBtow_tower_seed_size;i++) {
191  seedTow=wrkBtow_tower_seed[i];
192  seedEta=seedTow%BtowGeom::mxEtaBin;
193  seedPhi=seedTow/BtowGeom::mxEtaBin;
194 
195  //.... find first 2x2 above cluster thresh
196  if (seedEta < BtowGeom::mxEtaBin) {
197  clustET = sumET(seedPhi,seedEta);
198  if(clustET>par_clustEtThres) goto ACCEPT;
199  clustET = sumET(seedPhi-1,seedEta);
200  if(clustET>par_clustEtThres) goto ACCEPT;
201  }
202  if (seedEta > 0 ) {
203  clustET = sumET(seedPhi-1,seedEta-1);
204  if(clustET>par_clustEtThres) goto ACCEPT;
205  clustET = sumET(seedPhi,seedEta-1);
206  if(clustET>par_clustEtThres) goto ACCEPT;
207  }
208  }
209  //.... ABORT
210  btowEve->resultBlob.seedEt=0;
211  btowEve->resultBlob.clusterEt=0;
212  btowEve->resultBlob.seedEtaBin=0;
213  btowEve->resultBlob.seedPhiBin=0;
214  btowEve->resultBlob.trigger=0;
215  btowEve->seedET=0;
216  btowEve->clusterET=0;
217  btowEve->tkCompute=0;
218  return;
219 
220  ACCEPT:
221  // printf("B acc clustET=%f\n",clustET);
222  btowEve->seedET=wrkBtow_et[seedTow];
223  btowEve->clusterET=clustET;
224  btowEve->resultBlob.seedEt =(unsigned char)(wrkBtow_et[seedTow]*256.0/60.0);
225  btowEve->resultBlob.clusterEt=(unsigned char)(clustET*256.0/60.0);
226  btowEve->resultBlob.seedEtaBin=seedEta;
227  btowEve->resultBlob.seedPhiBin=seedPhi;
228  btowEve->resultBlob.trigger=2;
229  rdtscll( btowEve->tkCompute);
230  // printf("TTT1 %dl\n",ticks);
231 
232  return;
233 }
234 
235 
236 
237 
238 /* ========================================
239  ======================================== */
240 bool
241 L2wBemc2009::decisionUser(int token, int *myL2Result){
242  // INPUT: token + comput() results stored internally
243  // OUTPUT: yes/now + pointer to L2Result
244 
245  // get pointers to internal private event storage
246  L2wBemcEvent2009 *btowEve=mBtow+token;
247 
248  // prescaling decison ws already made in Virtual:
249  // if (par_RndAcceptPrescale>0 && mRandomAccept) btowEve->resultBlob.trigger+=1;
250  if (mRandomAccept) btowEve->resultBlob.trigger+=1;
251 
252  //...... some histos just for QA
253  if(btowEve->isFresh>L2wBemcEvent2009::kDataFresh) mhN->fill(6); // stale data, should never happen
254 
255  btowEve->isFresh++; // mark the data as stale
256 
257  if(btowEve->resultBlob.trigger&2) {
258  unsigned long long tkDecision;
259  rdtscll(tkDecision);
260  int tkDelta=tkDecision-btowEve->tkCompute;
261  //printf("TTT2 %d\n",tkDelta);
262  hA[1]->fill(tkDelta/1000);
263 
264  hA[3]->fill((int)btowEve->seedET);
265  hA[4]->fill((int)btowEve->clusterET);
266  hA[5]->fill((int)(100.*btowEve->seedET/btowEve->clusterET));
267  int jET=(int)(btowEve->clusterET/stepETH);
268  if( jET<=1) jET=1;
269  if( jET>=9) jET=9;
270  hA[10+jET]->fill(btowEve->resultBlob.seedEtaBin,btowEve->resultBlob.seedPhiBin);
271  hA[10]->fill(btowEve->resultBlob.seedEtaBin,btowEve->resultBlob.seedPhiBin);
272  }
273  // store final content in TRG data
274  memcpy(myL2Result,&(btowEve->resultBlob),sizeof(L2wResult2009));
275  if(par_dbg) L2wResult2009_print((L2wResult2009 *)&(btowEve->resultBlob));
276  //return btowEve->resultBlob.trigger;
277  return btowEve->resultBlob.trigger&2;
278 }
279 
280 
281 /* ========================================
282  ======================================== */
283 void
284 L2wBemc2009::finishRunUser() { /* called once at the end of the run */
285  // do whatever you want, log-file & histo-file are still open
286 
287  if (mLogFile){
288  fprintf(mLogFile,"finishRunUser-%s start\n",getName());
289  }
290 
291  // test1:
292  // hA[9]->fillW(21,1234);
293 
294  //......... scan for hot towers ....& repack 2D--> 1D(softId)
295 
296  for(int jh=10;jh<19;jh++)
297  {
298  const L2EmcDb::EmcCDbItem *xB=mDb->getByIndex(402);
299  const int *hiData = hA[jh]->getData();
300  int hotY=0,totY=0;
301  for(int i=0; i<EmcDbIndexMax; i++) {
302  const L2EmcDb::EmcCDbItem *x=mDb->getByIndex(i);
303  if(mDb->isEmpty(x)) continue;
304  if (!mDb->isBTOW(x) ) continue;
305  int ieta= (x->eta-1);
306  int iphi= (x->sec-1)*10 + x->sub-'a' ;
307  int softId=atoi(x->tube+2);
308  int ix=ieta+iphi*BtowGeom::mxEtaBin;
309  int yield=hiData[ix];
310  totY+=yield;
311  if(jh==10) hA[9]->fillW(softId,yield);// do it only for ET-integrated
312  if(hotY>yield)continue;
313  hotY=yield;
314  xB=x;
315  }
316 
317  int ieta= (xB->eta-1);
318  int iphi= (xB->sec-1)*10 + xB->sub-'a' ;
319  int softId=atoi(xB->tube+2);
320  fprintf(mLogFile,"#BTOW hot *candidate*, hist: %d, yield: %d, totYield: %d, softID: %d, crate: %d, chan: %d, name: %s, ieta: %d, iphi: %d\n", hA[jh]->getId(),hotY,totY,softId,xB->crate,xB->chan,xB->name,ieta,iphi);
321  }
322 
323 
324 
325 }
326 
327 
328 //=======================================
329 //=======================================
330 void
331 L2wBemc2009::createHisto() {
332  hA[1]=new L2Histo(1,(char*)"W Btow delTime (decision-compute); kTicks",300);
333 
334  hA[2]=new L2Histo(2,(char*)"W Btow-compute: # btow towers w/ energy /event; x: # BTOW towers; y: counts", 100);
335  hA[3]=new L2Histo(3,(char*)"W Btow-accepted: seeds ; Seed ET (GeV)", 70);
336  hA[4]=new L2Histo(4,(char*)"W Btow-accepted: clusters ; Cluster ET ET(GeV)", 70);
337  hA[5]=new L2Histo(5,(char*)"W Btow-accepted: cluster shape ;100*Seed Et/Cluster Et", 105);
338  // 6-8 free
339 
340  hA[9]=new L2Histo(9,(char*)"W Btow ...9 .....", 5000); // title in initRun
341  hA[10]=new L2Histo(10,(char*)"W Btow ...10 ...", BtowGeom::mxEtaBin, BtowGeom::mxPhiBin); // title in initRun
342 
343 
344  char tit[100];
345  for(int j=1;j<=9;j++){
346  int hid=10+j;
347  float x1=j*stepETH;
348  float x2=x1+stepETH;
349  if(j==1) x1=0;
350  if(j==9) x2=999;
351  sprintf(tit,"W Btow-accept: cluster %.0f <ET< %.0f GeV ; seed eta bin [-1,+1]; seed phi bin ~ TCP sector", x1,x2);
352 
353  hA[hid]=new L2Histo(hid,tit, BtowGeom::mxEtaBin, BtowGeom::mxPhiBin);
354  }
355 
356  //20-32 free
357 
358  // printf("L2-%s::createHisto() done\n",getName());
359 }
360 
361 //=======================================
362 //=======================================
363 void
364 L2wBemc2009::clearEvent(int token){
365  memset(wrkBtow_et,0,sizeof(wrkBtow_et));
366  wrkBtow_tower_seed_size=0;
367  memset(&(mBtow[token].resultBlob),0, sizeof(L2wResult2009));
368 }
369 
370 /* ========================================
371  ======================================== */
372 void
373 L2wBemc2009::print2(){ // full , local ADC array
374  int i;
375  printf("pr2-%s: ---BTOW ADC 2D array, only non-zero\n",getName());
376 
377  for(i=0;i<mxBtow;i++) {
378  if(wrkBtow_et[i]<=0) continue;
379  int rdo=mTower2rdo[i];
380  float et=wrkBtow_et[i];
381  printf(" btow: tower=%4d rdo=%4d et=%.3f \n",i,rdo,et);
382  }
383 
384 }
385 
386 /* ========================================
387  ======================================== */
388 void
389 L2wBemc2009::print3(){ // seed list
390  int i;
391  printf("pr3-%s: ---seed list, size=%d\n",getName(),wrkBtow_tower_seed_size);
392 
393  for(i=0;i<wrkBtow_tower_seed_size;i++) {
394  int tower=wrkBtow_tower_seed[i];
395  float et=wrkBtow_et[tower];
396  printf(" btow: i=%4d tower=%4d et=%.3f \n",i,tower,et);
397  }
398 
399 }
400 
401 #if 0
402 
403 /* ========================================
404  ======================================== */
405 void
406 L2wBemc2009::print4(int token, int hitSize){ // L2-algo input list
407  int i;
408  printf("print4 IS NOT Fully FUNCTIONAL **********************\n");
409  printf("pr1-%s: ---BTOW Sorted ADC list--- size=%d\n",getName(),hitSize);
410  //const HitTower *hit=globEve_btow_hit;
411  for(i=0;i< hitSize;i++) {
412  int adc=0;//(mEveStream_btow[token].get_hits()[wrkBtow_tower_index[i]]).adc;
413  int rdo=0;//(mEveStream_btow[token].get_hits()[wrkBtow_tower_index[i]]).rdo;
414  float et=wrkBtow_et[wrkBtow_tower_index[i]];
415  float ene=0;//(mEveStream_btow[token].get_hits()[wrkBtow_tower_index[i]]).ene;
416  printf(" tower=%2d ",wrkBtow_tower_index[i]);
417  printf(" btow: i=%2d rdo=%4d adc=%d et=%.3f ene=%.3f\n",i,rdo,adc,et,ene);
418  }
419 }
420 
421 #endif
422 /**********************************************************************
423  $Log: L2wBemc2009.cxx,v $
424  Revision 1.2 2009/11/19 15:48:49 balewski
425  add (char*) to many strings to make SL5 happ, few other adjustments
426 
427  Revision 1.1 2009/03/28 19:43:53 balewski
428  2009 code
429 
430 
431 
432 */
433 
434 
Generated by   doxygen 1.8.5