StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
L2wBemc2012.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: L2wBemc2012.cxx,v 1.4 2012/03/21 18:18:04 jml 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/L2EmcDb2012.h"
26  #include "../L2algoUtil/L2Histo.h"
27 #else
28  #include "../L2algoUtil/L2EmcDb2012.h"
29  #include "../L2algoUtil/L2Histo.h"
30  #include "../L2algoUtil/L2EmcGeom2012.h"
31 #endif
32 
33 #include "L2wBemc2012.h"
34 
35 //=================================================
36 //=================================================
37 L2wBemc2012::L2wBemc2012(const char* name, const char *uid, L2EmcDb2012* db, L2EmcGeom2012 *geoX, char* outDir, int resOff) : L2VirtualAlgo2012( name,uid, 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(L2wResult2012)!= L2wResult2012::mySizeChar)
51  criticalError((char*)"L2wBemc has failed consistency check. sizeof(L2wResult2012)!= L2wResult2012::mySizeChar");
52 
53 }
54 
55 /* ========================================
56  ======================================== */
57 int
58 L2wBemc2012::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 L2EmcDb2012::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 L2wBemc2012::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  tow+=BtowGeom::mxEtaBin;
140  tow%=maxTowers;
141 
142  sum+=wrkBtow_et[tow];
143 
144  towPlusOne = tow+1;
145  towPlusOne%= maxTowers;
146  sum+=wrkBtow_et[towPlusOne];
147 
148  return sum;
149 }
150 
151 
152 /* ========================================
153  ======================================== */
154 void
155 L2wBemc2012::computeUser(int token){
156 
157  clearEvent(token);
158 
159  // ------ PROJECT INPUT LIST TO 2D ARRAY AND SCAN FOR SEED TOWERS ----
160  int i;
161 
162  // access inpute list
163  const HitTower1 *hit=mEveStream_btow[token].get_hits();
164  const int hitSize=mEveStream_btow[token].get_hitSize();
165  // preapre output list
166  L2wBemcEvent2012 *btowEve=mBtow+token;
167 
168  for(i=0;i< hitSize;i++,hit++) {
169  int tower=mRdo2tower[hit->rdo];
170  wrkBtow_et[tower]=hit->et;
171  if(hit->et<par_seedEtThres)continue;
172  wrkBtow_tower_seed[wrkBtow_tower_seed_size++]=tower;
173  }
174  hA[2]->fill(hitSize);
175  int seedTow=-1,seedEta=-1,seedPhi=-1;
176  float clustET=0;
177  btowEve->isFresh=L2wBemcEvent2012::kDataFresh;
178 
179  // ----------- FIND 2x2 CLUSTER AROUND EVERY SEED -----
180  for(i=0; i<wrkBtow_tower_seed_size;i++) {
181  seedTow=wrkBtow_tower_seed[i];
182  seedEta=seedTow%BtowGeom::mxEtaBin;
183  seedPhi=seedTow/BtowGeom::mxEtaBin;
184 
185  //.... find first 2x2 above cluster thresh
186  if (seedEta < BtowGeom::mxEtaBin) {
187  clustET = sumET(seedPhi,seedEta);
188  if(clustET>par_clustEtThres) goto ACCEPT;
189  clustET = sumET(seedPhi-1,seedEta);
190  if(clustET>par_clustEtThres) goto ACCEPT;
191  }
192  if (seedEta > 0 ) {
193  clustET = sumET(seedPhi-1,seedEta-1);
194  if(clustET>par_clustEtThres) goto ACCEPT;
195  clustET = sumET(seedPhi,seedEta-1);
196  if(clustET>par_clustEtThres) goto ACCEPT;
197  }
198  }
199  //.... ABORT
200  btowEve->resultBlob.seedEt=0;
201  btowEve->resultBlob.clusterEt=0;
202  btowEve->resultBlob.seedEtaBin=0;
203  btowEve->resultBlob.seedPhiBin=0;
204  btowEve->resultBlob.trigger=0;
205  btowEve->seedET=0;
206  btowEve->clusterET=0;
207  btowEve->tkCompute=0;
208  return;
209 
210  ACCEPT:
211  btowEve->seedET=wrkBtow_et[seedTow];
212  btowEve->clusterET=clustET;
213  btowEve->resultBlob.seedEt =(unsigned char)(wrkBtow_et[seedTow]*256.0/60.0);
214  btowEve->resultBlob.clusterEt=(unsigned char)(clustET*256.0/60.0);
215  btowEve->resultBlob.seedEtaBin=seedEta;
216  btowEve->resultBlob.seedPhiBin=seedPhi;
217  btowEve->resultBlob.trigger=2;
218  rdtscll( btowEve->tkCompute);
219 
220  return;
221 }
222 
223 
224 
225 
226 /* ========================================
227  ======================================== */
228 bool
229 L2wBemc2012::decisionUser(int token, int *myL2Result){
230  // INPUT: token + comput() results stored internally
231  // OUTPUT: yes/now + pointer to L2Result
232 
233  // get pointers to internal private event storage
234  L2wBemcEvent2012 *btowEve=mBtow+token;
235 
236  // prescaling decison ws already made in Virtual:
237  // if (par_RndAcceptPrescale>0 && mRandomAccept) btowEve->resultBlob.trigger+=1;
238  if (mRandomAccept) btowEve->resultBlob.trigger+=1;
239 
240  //...... some histos just for QA
241  if(btowEve->isFresh>L2wBemcEvent2012::kDataFresh) mhN->fill(6); // stale data, should never happen
242 
243  btowEve->isFresh++; // mark the data as stale
244 
245  if(btowEve->resultBlob.trigger&2) {
246  unsigned long long tkDecision;
247  rdtscll(tkDecision);
248  int tkDelta=tkDecision-btowEve->tkCompute;
249  hA[1]->fill(tkDelta/1000);
250 
251  hA[3]->fill((int)btowEve->seedET);
252  hA[4]->fill((int)btowEve->clusterET);
253  hA[5]->fill((int)(100.*btowEve->seedET/btowEve->clusterET));
254  int jET=(int)(btowEve->clusterET/stepETH);
255  if( jET<=1) jET=1;
256  if( jET>=9) jET=9;
257  hA[10+jET]->fill(btowEve->resultBlob.seedEtaBin,btowEve->resultBlob.seedPhiBin);
258  hA[10]->fill(btowEve->resultBlob.seedEtaBin,btowEve->resultBlob.seedPhiBin);
259  }
260  // store final content in TRG data
261  memcpy(myL2Result,&(btowEve->resultBlob),sizeof(L2wResult2012));
262  if(par_dbg) L2wResult2012_print((L2wResult2012 *)&(btowEve->resultBlob));
263  //return btowEve->resultBlob.trigger;
264  return btowEve->resultBlob.trigger&2;
265 }
266 
267 
268 /* ========================================
269  ======================================== */
270 void
271 L2wBemc2012::finishRunUser() { /* called once at the end of the run */
272  // do whatever you want, log-file & histo-file are still open
273 
274  if (mLogFile){
275  fprintf(mLogFile,"finishRunUser-%s start\n",getName());
276  }
277 
278  // test1:
279  // hA[9]->fillW(21,1234);
280 
281  //......... scan for hot towers ....& repack 2D--> 1D(softId)
282 
283  for(int jh=10;jh<19;jh++)
284  {
285  const L2EmcDb2012::EmcCDbItem *xB=mDb->getByIndex(402);
286  const int *hiData = hA[jh]->getData();
287  int hotY=0,totY=0;
288  for(int 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 ieta= (x->eta-1);
293  int iphi= (x->sec-1)*10 + x->sub-'a' ;
294  int softId=atoi(x->tube+2);
295  int ix=ieta+iphi*BtowGeom::mxEtaBin;
296  int yield=hiData[ix];
297  totY+=yield;
298  if(jh==10) hA[9]->fillW(softId,yield);// do it only for ET-integrated
299  if(hotY>yield)continue;
300  hotY=yield;
301  xB=x;
302  }
303 
304  int ieta= (xB->eta-1);
305  int iphi= (xB->sec-1)*10 + xB->sub-'a' ;
306  int softId=atoi(xB->tube+2);
307  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);
308  }
309 
310 
311 
312 }
313 
314 
315 //=======================================
316 //=======================================
317 void
318 L2wBemc2012::createHisto() {
319  hA[1]=new L2Histo(1,(char*)"W Btow delTime (decision-compute); kTicks",300);
320 
321  hA[2]=new L2Histo(2,(char*)"W Btow-compute: # btow towers w/ energy /event; x: # BTOW towers; y: counts", 100);
322  hA[3]=new L2Histo(3,(char*)"W Btow-accepted: seeds ; Seed ET (GeV)", 70);
323  hA[4]=new L2Histo(4,(char*)"W Btow-accepted: clusters ; Cluster ET ET(GeV)", 70);
324  hA[5]=new L2Histo(5,(char*)"W Btow-accepted: cluster shape ;100*Seed Et/Cluster Et", 105);
325  // 6-8 free
326 
327  hA[9]=new L2Histo(9,(char*)"W Btow ...9 .....", 5000); // title in initRun
328  hA[10]=new L2Histo(10,(char*)"W Btow ...10 ...", BtowGeom::mxEtaBin, BtowGeom::mxPhiBin); // title in initRun
329 
330 
331  char tit[100];
332  for(int j=1;j<=9;j++){
333  int hid=10+j;
334  float x1=j*stepETH;
335  float x2=x1+stepETH;
336  if(j==1) x1=0;
337  if(j==9) x2=999;
338  sprintf(tit,"W Btow-accept: cluster %.0f <ET< %.0f GeV ; seed eta bin [-1,+1]; seed phi bin ~ TCP sector", x1,x2);
339 
340  hA[hid]=new L2Histo(hid,tit, BtowGeom::mxEtaBin, BtowGeom::mxPhiBin);
341  }
342 
343  //20-32 free
344 
345 }
346 
347 //=======================================
348 //=======================================
349 void
350 L2wBemc2012::clearEvent(int token){
351  memset(wrkBtow_et,0,sizeof(wrkBtow_et));
352  wrkBtow_tower_seed_size=0;
353  memset(&(mBtow[token].resultBlob),0, sizeof(L2wResult2012));
354 }
355 
356 /* ========================================
357  ======================================== */
358 void
359 L2wBemc2012::print2(){ // full , local ADC array
360  int i;
361  printf("pr2-%s: ---BTOW ADC 2D array, only non-zero\n",getName());
362 
363  for(i=0;i<mxBtow;i++) {
364  if(wrkBtow_et[i]<=0) continue;
365  int rdo=mTower2rdo[i];
366  float et=wrkBtow_et[i];
367  printf(" btow: tower=%4d rdo=%4d et=%.3f \n",i,rdo,et);
368  }
369 
370 }
371 
372 /* ========================================
373  ======================================== */
374 void
375 L2wBemc2012::print3(){ // seed list
376  int i;
377  printf("pr3-%s: ---seed list, size=%d\n",getName(),wrkBtow_tower_seed_size);
378 
379  for(i=0;i<wrkBtow_tower_seed_size;i++) {
380  int tower=wrkBtow_tower_seed[i];
381  float et=wrkBtow_et[tower];
382  printf(" btow: i=%4d tower=%4d et=%.3f \n",i,tower,et);
383  }
384 
385 }
386 
387 #if 0
388 
389 /* ========================================
390  ======================================== */
391 void
392 L2wBemc2012::print4(int token, int hitSize){ // L2-algo input list
393  int i;
394  printf("print4 IS NOT Fully FUNCTIONAL **********************\n");
395  printf("pr1-%s: ---BTOW Sorted ADC list--- size=%d\n",getName(),hitSize);
396  //const HitTower *hit=globEve_btow_hit;
397  for(i=0;i< hitSize;i++) {
398  int adc=0;//(mEveStream_btow[token].get_hits()[wrkBtow_tower_index[i]]).adc;
399  int rdo=0;//(mEveStream_btow[token].get_hits()[wrkBtow_tower_index[i]]).rdo;
400  float et=wrkBtow_et[wrkBtow_tower_index[i]];
401  float ene=0;//(mEveStream_btow[token].get_hits()[wrkBtow_tower_index[i]]).ene;
402  printf(" tower=%2d ",wrkBtow_tower_index[i]);
403  printf(" btow: i=%2d rdo=%4d adc=%d et=%.3f ene=%.3f\n",i,rdo,adc,et,ene);
404  }
405 }
406 
407 #endif
408 /**********************************************************************
409  $Log: L2wBemc2012.cxx,v $
410  Revision 1.4 2012/03/21 18:18:04 jml
411  got rid of printfs from 2012 files
412 
413  Revision 1.3 2011/10/19 16:12:12 jml
414  more 2012 stuff
415 
416  Revision 1.2 2011/10/19 15:39:45 jml
417  2012
418 
419  Revision 1.1 2011/10/18 15:11:44 jml
420  adding 2012 algorithms
421 
422  Revision 1.2 2009/11/19 15:48:49 balewski
423  add (char*) to many strings to make SL5 happ, few other adjustments
424 
425  Revision 1.1 2009/03/28 19:43:53 balewski
426  2009 code
427 
428 
429 
430 */
431 
432 
Generated by   doxygen 1.8.5