StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
L2exampleAlgo08.cxx
1 #include <stdio.h>
2 #include <assert.h>
3 #include <string.h>
4 #include <stdlib.h>
5 #include <time.h>
6 #include <math.h>
7 
8 /*********************************************************************
9  * $Id: L2exampleAlgo08.cxx,v 1.6 2008/01/30 21:56:40 balewski Exp $
10  * \author Jan Balewski,MIT , 2008
11  *********************************************************************
12  * Descripion: see .h
13  *********************************************************************
14  */
15 
16 
17 #ifdef IS_REAL_L2 //in l2-ana environment
18  #include "../L2algoUtil/L2EmcDb.h"
19  #include "../L2algoUtil/L2Histo.h"
20 #else
21  #include "../L2algoUtil/L2EmcDb.h"
22  #include "../L2algoUtil/L2Histo.h"
23  #include "../L2algoUtil/L2EmcGeom.h"
24 #endif
25 
26 #include "L2exampleAlgo08.h"
27 
28 
29 //=================================================
30 //=================================================
31 L2exampleAlgo08::L2exampleAlgo08(const char* name, L2EmcDb* db, L2EmcGeom *geoX, char* outDir) : L2VirtualAlgo2008( name, db, outDir) {
32  /* called one per days
33  all memory allocation must be done here
34  */
35 
36  mGeom=geoX; assert(mGeom); // avaliable but not used in this example
37  setMaxHist(16); // set upper range, I uses only 2^N -it is easier to remember
38  createHisto();
39 
40  //------- self-consistency checks, should never fail
41  // printf("ZZ %d %d\n", sizeof(L2exampleResult08), L2exampleResult08::mySizeChar);
42  assert(sizeof(L2exampleResult08)== L2exampleResult08::mySizeChar);
43 }
44 
45 /* ========================================
46  ======================================== */
47 int
48 L2exampleAlgo08::initRunUser( int runNo, int *rc_ints, float *rc_floats) {
49 
50  // unpack params from run control GUI
51  par_dbg = rc_ints[0];
52  par_seedEtThres = rc_floats[0];
53  par_clusterEtThres= rc_floats[1];
54  par_eventEtThres = rc_floats[2];
55 
56  // verify consistency of input params
57  int kBad=0;
58  kBad+=0x00001 * (par_seedEtThres<1.0);
59  kBad+=0x00002 * (par_clusterEtThres<par_seedEtThres);
60 
61  if (mLogFile) {
62  fprintf(mLogFile,"L2%s algorithm initRun(R=%d), compiled: %s , %s\n params:\n",getName(),mRunNumber,__DATE__,__TIME__);
63  fprintf(mLogFile," - use Thres/GeV seed=%.2f, clusterList=%.2f debug=%d\n", par_seedEtThres,par_clusterEtThres, par_dbg);
64  fprintf(mLogFile," - accept event cluster Thres/GeV=%.2f\n",par_eventEtThres);
65  fprintf(mLogFile,"initRun() params checked for consistency, Error flag=0x%04x\n",kBad);
66  }
67 
68  if(kBad) return kBad;
69 
70 
71  // clear content of all histograms & token-dependet memory
72  int i;
73  for (i=0; i<mxHA;i++) if(hA[i])hA[i]->reset();
74  memset(mBtow,0,sizeof(mBtow));
75 
76  // update tiltles of histos
77  char txt[1000];
78  sprintf(txt,"BTOW-compute: #seed towers ET>%.2f GeV / event; x: # BTOW towers; y: counts",par_seedEtThres);
79  hA[3]->setTitle(txt);
80 
81  sprintf(txt,"BTOW-decision: #cluster ET>%.2f GeV / event ; x: # BTOW towers; y: counts",par_clusterEtThres);
82  hA[4]->setTitle(txt);
83 
84  sprintf(txt,"BTOW-decision: acc cluster ET>%.2f GeV; x:cluster ET(GeV) ; y: counts",par_eventEtThres);
85  hA[4]->setTitle(txt);
86 
87 
88 
89  for ( int index=0; index<EmcDbIndexMax; index++ )
90  {
91  const L2EmcDb::EmcCDbItem *x = mDb->getByIndex(index);
92  if ( x==0 ) continue;
93  if ( !mDb->isBTOW(x) ) continue;
94  int sec = x->sec - 1;
95  int sub = 8192;
96  sub = x->sub - 'a';
97  int eta = x->eta - 1;
98  int phi = BtowGeom::mxSubs *sec + sub;
99  int tow = BtowGeom::mxEtaBin *phi + eta; // phi- changes faster
100  int rdo = x->rdo;
101  assert(tow>=0); assert(tow<=mxBtow);
102  assert(rdo>=0); assert(rdo<=mxBtow);
103 
104  mTower2rdo[ tow ] = rdo; // returns rdo channel for given tower
105  mRdo2tower[ rdo ] = tow; // returns tower for given rdo channel
106  }
107  return 0; //OK
108 
109 }
110 
111 /* ========================================
112  ======================================== */
113 float
114 L2exampleAlgo08::sumET(int phi, int eta) {
115  int tow = BtowGeom::mxEtaBin *phi + eta; // phi- changes faster
116  float sum=wrkBtow_et[tow]+wrkBtow_et[tow+1];
117  // printf("tow : %d, %d --> %f %f \n",tow, tow+1,wrkBtow_et[tow],wrkBtow_et[tow+1]);
118  tow+=BtowGeom::mxEtaBin;
119  // printf("tow : %d, %d --> %f %f \n",tow, tow+1,wrkBtow_et[tow],wrkBtow_et[tow+1]);
120  sum+=wrkBtow_et[tow]+wrkBtow_et[tow+1];
121  return sum;
122 }
123 
124 
125 /* ========================================
126  ======================================== */
127 void
128 L2exampleAlgo08::computeUser(int token){
129  // token range is guaranteed by virtual08-class
130 
131  /* primitive 2x2 cluster finder:
132  - ignores seed towers at eta or phi boundary
133  - doublcounts if 2 seeds are neighbours
134  */
135 
136  clearEvent(token);
137 
138  // ------ PROJECT INPUT LIST TO 2D ARRAY AND SCAN FOR SEED TOWERS ----
139  int i;
140  // printf("L2-%s-compute: ---BTOW ADC list--- size=%d\n",getName(),*globEve_btow_hitSize);
141 
142  const HitTower1 *hit=mEveStream_btow[token].get_hits();
143  const int hitSize=mEveStream_btow[token].get_hitSize();
144  for(i=0;i< hitSize;i++,hit++) {
145  int tower=mRdo2tower[hit->rdo];
146  wrkBtow_et[tower]=hit->et;
147  if(hit->et<par_seedEtThres)continue;
148  wrkBtow_tower_seed[wrkBtow_tower_seed_size++]=tower;
149  }
150  hA[2]->fill(hitSize);
151  hA[3]->fill(wrkBtow_tower_seed_size);
152 
153  // ----------- FIND 2x2 CLUSTER AROUND EVERY SEED -----
154 
155  // get pointers to internal private event storage
156  L2exampleEvent08 *btowEve=mBtow+token;
157 
158  // compute & store values
159  for(i=0;i<wrkBtow_tower_seed_size;i++) {
160  int seedTow=wrkBtow_tower_seed[i];
161  int seedEta=seedTow%BtowGeom::mxEtaBin;
162  int seedPhi=seedTow/BtowGeom::mxEtaBin;
163  // float seedET= wrkBtow_et[seedTow];
164  // printf("sumE seed ET=%.3f tow=%d phiBin=%d etaBin=%d\n",seedET,seedTow,seedPhi,seedEta);
165 
166  // ........drop seeds close to boundaries, fix it
167  if(seedEta==0 || seedEta==BtowGeom::mxEtaBin-1) continue;
168  if(seedPhi==0 || seedPhi==BtowGeom::mxPhiBin-1) continue;
169  // now every seed has 4 2x2 clusters
170 
171  /* Pibero's idea to avoid double counting:
172  "drop seed if a neighbour tower has more energy"
173  has not been implemented - sounds great!
174  */
175 
176  //.... find max of 4 possible clusters
177  float maxET=sumET(seedPhi,seedEta);
178  float sum=sumET(seedPhi-1,seedEta);
179  if(maxET<sum) maxET=sum;
180  sum=sumET(seedPhi-1,seedEta-1);
181  if(maxET<sum) maxET=sum;
182  sum=sumET(seedPhi,seedEta-1);
183  if(maxET<sum) maxET=sum;
184  // printf("CCC clustET=%f \n",maxET);
185  if(maxET<par_clusterEtThres)continue;
186  if(btowEve->size>=L2exampleEvent08::mxClust) continue; // overflow protection
187  //........record largest cluster for this seed....
188  btowEve->clusterET[btowEve->size++]=maxET;
189  }// end of cluster search
190 
191  btowEve->isFresh=L2exampleEvent08::kDataFresh;
192  // printf("compuzzzzzzzzzzzzzzzzz s=%d tkn=%d\n",*btow_clusterET_size,token);
193 
194  // debugging should be off for any time critical computation
195  if(par_dbg>0){
196  printf("dbg=%s, btow-adcL-size=%d\n",getName(),hitSize);
197  // if(par_dbg>0) print1();
198  print2();
199  print3();// tmp, must have token
200  }
201 }
202 
203 
204 /* ========================================
205  ======================================== */
206 bool
207 L2exampleAlgo08::decisionUser(int token, void **myL2Result){
208  // INPUT: token + comput() results stored internally
209  // OUTPUT: yes/now + pointer to L2Result
210 
211  // get pointers to internal private event storage
212  L2exampleEvent08 *btowEve=mBtow+token;
213  (*myL2Result)=&(btowEve->resultBlob);
214  // printf("eee %p\n",*myL2Result);
215  int ic;
216 
217  //...... some histos just for fun
218  if(btowEve->size>= L2exampleEvent08::mxClust) mhN->fill(5); // was overflow
219  if(btowEve->isFresh>L2exampleEvent08::kDataFresh) mhN->fill(6); // stale data
220  btowEve->isFresh++; // mark the data as stale
221 
222  hA[4]->fill(btowEve->size);
223 
224  for(ic=0;ic<btowEve->size;ic++) {
225  float clustET=btowEve->clusterET[ic];
226  hA[5]->fill((int)clustET);
227  // printf("DDD clust ET=%.3f\n", clustET);
228  if(clustET<par_eventEtThres) continue;
229  hA[6]->fill((int)clustET);
230  }
231  // printf("clustzzzzzzzzzzzzzzzzz s=%d tkn=%d\n",btow_clusterET_size,token);
232 
233  //........ fill (some) L2Result
234  btowEve->resultBlob.kTicksCompute=mComputeTimeDiff[token]/1000;
235  btowEve->resultBlob.decision=0;
236  btowEve->resultBlob.numberOfL2Clust=btowEve->size;
237 
238  //........ compute the final decision
239  for(ic=0;ic<btowEve->size;ic++) {
240  if(btowEve->clusterET[ic]<par_eventEtThres) continue;
241  btowEve->resultBlob.decision=1;
242  btowEve->resultBlob.clusterET=btowEve->clusterET[ic];
243  goto accepted;
244  }
245  return false;
246 
247  /****************/
248  accepted:
249  /****************/
250  if(btowEve->size>= L2exampleEvent08::mxClust) mhN->fill(15);
251  return true;
252 
253 }
254 
255 
256 /* ========================================
257  ======================================== */
258 void
259 L2exampleAlgo08::finishRunUser() { /* called once at the end of the run */
260  // do whatever you want, log-file & histo-file are still open
261 
262  if (mLogFile){
263  fprintf(mLogFile,"finishRunUser-%s bhla bhla\n",getName());
264  }
265 
266 }
267 
268 
269 //=======================================
270 //=======================================
271 void
272 L2exampleAlgo08::createHisto() {
273  memset(hA,0,sizeof(hA));
274 
275  hA[2]=new L2Histo(2,"BTOW-compute: #towers w/ energy /event; x: # BTOW towers; y: counts", 100);
276  hA[3]=new L2Histo(3,"BTOW-compute: #seed ....... ", 100); // title in initRun
277  hA[4]=new L2Histo(4,"BTOW-decision: #clust ....... ", 50); // title in initRun
278 
279  hA[5]=new L2Histo(5,"BTOW-decision: any cluster ; x: ET(GeV)", 30);
280  hA[6]=new L2Histo(6,"BTOW-decision: accepted clust ... ; x: ET(GeV)", 30);// title in initRun
281 
282  // printf("L2-%s::createHisto() done\n",getName());
283 }
284 
285 //=======================================
286 //=======================================
287 void
288 L2exampleAlgo08::clearEvent(int token){
289  memset(wrkBtow_et,0,sizeof(wrkBtow_et));
290  memset(wrkBtow_tower_seed,0,sizeof(wrkBtow_tower_seed));
291  wrkBtow_tower_seed_size=0;
292  mBtow[token].size=0; // do not set 'fresh-flag here, only after compute() finishes
293  memset(&(mBtow[token].resultBlob),0, sizeof(L2exampleResult08));
294 }
295 
296 /* ========================================
297  ======================================== */
298 void
299 L2exampleAlgo08::print2(){ // full , local ADC array
300  int i;
301  printf("pr2-%s: ---BTOW ADC 2D array, only non-zero\n",getName());
302 
303  for(i=0;i<mxBtow;i++) {
304  if(wrkBtow_et[i]<=0) continue;
305  int rdo=mTower2rdo[i];
306  float et=wrkBtow_et[i];
307  printf(" btow: tower=%4d rdo=%4d et=%.3f \n",i,rdo,et);
308  }
309 
310 }
311 
312 /* ========================================
313  ======================================== */
314 void
315 L2exampleAlgo08::print3(){ // seed list
316  int i;
317  printf("pr3-%s: ---seed list, size=%d\n",getName(),wrkBtow_tower_seed_size);
318 
319  for(i=0;i<wrkBtow_tower_seed_size;i++) {
320  int tower=wrkBtow_tower_seed[i];
321  float et=wrkBtow_et[tower];
322  printf(" btow: i=%4d tower=%4d et=%.3f \n",i,tower,et);
323  }
324 
325 }
326 
327 /**********************************************************************
328  $Log: L2exampleAlgo08.cxx,v $
329  Revision 1.6 2008/01/30 21:56:40 balewski
330  E+B high-enery-filter L2-algo fuly functional
331 
332  Revision 1.5 2008/01/30 00:47:17 balewski
333  Added L2-Etow-calib
334 
335  Revision 1.4 2008/01/18 23:29:13 balewski
336  now L2result is exported
337 
338  Revision 1.3 2008/01/17 23:15:51 balewski
339  bug in token-addressed memory fixed
340 
341  Revision 1.2 2008/01/16 23:32:35 balewski
342  toward token dependent compute()
343 
344  Revision 1.1 2007/12/19 02:30:18 balewski
345  new L2-btow-calib-2008
346 
347 
348 
349 */
350 
351 
Generated by   doxygen 1.8.5