StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
L2gammaAlgo.cxx
1 #include <stdio.h>
2 
3 
4 #include <stdio.h>
5 #include "L2gammaAlgo.h"
6 #include "L2gammaResult2006.h"
7 
8 #include "bemcConstants.h"
9 #include "eemcConstants.h"
10 
11 void L2gammaAlgo::setLogFile( const char *fname){ mLogFile=fopen(fname,"w"); }
12 void L2gammaAlgo::setHistFile( const char *fname ){ mHistFile=fopen(fname,"w"); }
13 
14 //#define DEBUG
15 //#define DEBUG_PATCH_THRESHOLDS
16 // ----------------------------------------------------------------------------
17 L2gammaAlgo::L2gammaAlgo(const char* name, L2EmcDb* db, char* outDir, int resOff)
18  : L2VirtualAlgo( name, db, outDir, resOff)
19 {
20 
21  int geom, thresh; // temporary fix,JB
22 #if 1
23  if(strstr(name,"etow_gamma")) {
24  geom= L2gammaAlgo::kEEmcAlgo; thresh = L2gammaAlgo::kThresh1;
25  } else if(strstr(name,"etow_ht2")) {
26  geom= L2gammaAlgo::kEEmcAlgo; thresh = L2gammaAlgo::kThresh2;
27  } else if(strstr(name,"btow_gamma")) {
28  geom= L2gammaAlgo::kBEmcAlgo; thresh = L2gammaAlgo::kThresh1;
29  } else if(strstr(name,"btow_ht2")) {
30  geom= L2gammaAlgo::kBEmcAlgo; thresh = L2gammaAlgo::kThresh2;
31  } else {
32  printf("undefined config of L2gamma-algo ->%s<-, abort L2main\n",name);
33  assert(1==2);
34  }
35 #endif
36 
37 
38  printf("L2gammaEmCall2006 instantiated geom=%d thresh=%d logpath=%s\n",geom,thresh,mOutDir);
39  int I_par[]={
40  1, // 0=bemc, 1=eemc
41  0, // prescale
42  0, // free
43  0, // free
44  0 // free
45  };
46  float F_par[]={
47  2.50, // high tower threshold
48  3.50, // 3x3 patch threshold
49  0., // free
50  0., // free
51  0. // free
52  };
53  mThresholdLevel=thresh;
54  init(geom,I_par,F_par);
55  mL2input=0;
56  mPrescale=0;
57 
58 }
59 
60 
61 // ----------------------------------------------------------------------------
62 //
63 // perform initialization at L2 main program startup not
64 // done in the constructor
65 //
66 void L2gammaAlgo::init(int geom, int I_par[5], float F_par[5] )
67 {
68 
69  printf("+ init geom=%d\n",geom);
70 
71  mLogFile=0;
72  mHistFile=0;
73  mUseOfflineGains = I_par[0];
74  mUseBbc=false;
75  setTowerThreshold( F_par[0] );
76  setPatchThreshold( F_par[1] );
77 
78  for ( int i=0;i<5;i++ )
79  {
80  mDefaultI_par[i]=I_par[i];
81  mDefaultF_par[i]=F_par[i];
82  }
83 
84  if ( I_par[1] ) mPrescale=I_par[1];
85 
86  /*
87  * Select endcap or barrel based on the geom flag.
88  */
89  mEEmc=geom;mBEmc=!geom;
90 
91  mNumEtas = (mEEmc)? kEEmcNumEtas : kBEmcNumEtas; /* number of etabins */
92  mNumPhis = (mEEmc)? kEEmcNumPhis : kBEmcNumPhis; /* number of phi bins */
93  mNumSubs = (mEEmc)? kEEmcNumSubs : kBEmcNumSubs; /* number of subsectors */
94  mNumSecs = (mEEmc)? kEEmcNumSecs : kBEmcNumSecs; /* number of sectors */
95 
96  mNumTower = (mEEmc)? kEEmcNumTower : kBEmcNumTower;
97  mNumClust = mNumTower;
98  mNumRdo = mNumTower;
99  mEtaBins = (mEEmc)? kEEmcEtaBins : kBEmcEtaBins;
100 
101  mMaxADC = (mEEmc)? kEEmcMaxADC : kBEmcMaxADC;
102  mMaxET = (mEEmc)? kEEmcMaxET : kBEmcMaxET;
103  mIdealGainT = (mEEmc)? kEEmcIdealGainT : kBEmcIdealGainT;
104 
105  const char *names[]={"bemc","eemc"};
106 
107  printf("L2gammaAlgo v0.92 (prociutto) \n");
108  printf("registering new threshold for %4s \n",names[geom]);
109 
110 #ifdef DEBUGIT
111  printf(" \n");
112 
113  printf("============================================================================\n");
114  // Jason, I am commenting this for the last time.
115  // there is no benefit of rinting this on the L2 console.
116  // If this printout of any value write it to a file saved on /data/... disk as all other outputs. Then at least you can retrieve it after the run.
117  //Or activate this only in the debug mode - you can have GUI param for debug - if you want
118 
119  printf("algorithm residing at %p\n",this);
120  printf("l2algo used as input at %p\n",mL2input);
121  printf("\n");
122  printf("\n");
123  printf("mNumEtas = %-2d\n",mNumEtas);
124  printf("mNumPhis = %-4d\n",mNumPhis);
125  printf("mNumSubs = %-4d\n",mNumSubs);
126  printf("mNumSecs = %-4d\n",mNumSecs);
127  printf("mNumTower = %-4d\n",mNumTower);
128  printf("mNumClust = %-4d\n",mNumClust);
129  printf("mNumRdo = %-4d\n",mNumRdo);
130  printf("mEtaBins = {");
131  for ( int i=0;i<mNumEtas+1;i++ ) {
132  printf("%4.2f, ",mEtaBins[i]);
133  if ( !((i+1)%4) ) printf("\n ");
134  }
135  printf("X}\n");
136  printf("mMaxADC = %4d\n",mMaxADC);
137  printf("mMaxET = %4.1f\n",mMaxET);
138  printf("mIdealGainT = %5.3f\n",mIdealGainT);
139  printf("\n");
140  printf("============================================================================\n");
141 
142 #endif
143 
144  // book "histograms" for QA (clears from previous run)
145  jbook();
146 
147 }
148 
149 
150 
151 // ----------------------------------------------------------------------------
152 //
153 // runtime initialization
154 //
155 // at this stage, time is not overly critical. so we will take
156 // the time to do most of the hard work. we will do two things:
157 //
158 // 1) establish a high tower threshold for all 720 towers in
159 // terms of adc. This will be the high tower threshold (in
160 // adc) plus pedestal. Thus we will not need to do a ped
161 // subtraction in the event loop.
162 //
163 // 2) We implement a 3x3 tower patch around each tower. For
164 // each rdo channel we determine the number of (physically)
165 // adjacent towers and
166 //
167 // 3) An optional parameter (use offline gains) will correct
168 // the ADC thresholds for gain variations measured offline.
169 // When applied to patch sum thresholds, we will also need
170 // to apply correction to the towers in the patch sums.
171 //
172 int L2gammaAlgo::initRun( int run )
173 {
174  return initRun(run,mDefaultI_par,mDefaultF_par);
175 }
176 int L2gammaAlgo::initRun( char *myname, int run, int I_par[5], float F_par[5] )
177 {
178  return initRun( run, I_par, F_par );
179 }
180 int L2gammaAlgo::initRun( int run, int I_par[5], float F_par[5] )
181 {
182 
183 
184  printf("%s ::initRun() for run=%d\n",mName,run);
185 
186  if ( mDb->initRun(run) ) return -7; // this must be in, if your algos are the only one in the game, _you_ must initialize DB for the new run,JB
187 
188  mRunNumber = run;
189 
190  // recreate histograms (clears)
191  //$$$ jbook();
192  jclear();
193 
194  mUseOfflineGains = I_par[0];
195  if ( I_par[1] ) mPrescale = I_par[1];
196  mUseBbc=false;
197  setTowerThreshold( F_par[0] );
198  setPatchThreshold( F_par[1] );
199 
200 
201  /***********************************************************************
202  * test validity of input parameters
203  */
204  {
205 
206  if ( F_par[1] < 0. ) {
207  printf("L2gammaAlgo::ERROR -- patch threshold < 0\n");
208  return 1;
209  }
210  if ( F_par[1] > 2.0 * F_par[0] ) {
211  printf("L2gammaAlgo::ERROR -- patch threshold > 2x tower threshold\n");
212  return 1;
213  }
214 
215  if ( mL2input )
216  {
217  if ( F_par[0] < mL2input->getPatchThreshold() ) {
218  printf("L2gammaAlgo::WARNING -- l2input to this algorithm has a patch threshold < tower threshold\n");
219  //return 1;
220  }
221  }
222 
223  if ( I_par[0] < 0 || I_par[0] > 1 )
224  {
225  printf("L2gammaAlgo::ERROR -- I_par[0] should be 0 (use ideal gains) or 1 (use online gains)\n");
226  }
227 
228  }/* valid */
229 
230 
231  // close histogram and log files and reopen w/ new paths
232  if ( mLogFile && mLogFile != stdout ) fclose(mLogFile);
233  if ( mHistFile ) fclose(mHistFile);
234 
235  //const char *names[]={"bemc","eemc"};
236 
237  char clog[128];
238  sprintf(clog,"%s/run%d.l2%s.out",mOutDir,run,mName);
239  //printf("%s=\n",clog);
240  //setHistFile(chis);
241  // setLogFile(clog);
242  mLogFile=fopen(clog,"w");
243  // assert(mLogFile);
244  if(mLogFile==0) printf("%s ::initRun() for run=%d UNABLE to open log-file=%s=, it is not fatal, continue initialization\n",mName,run,clog);
245 
246  if(mLogFile) {
247  fprintf(mLogFile,"\n\n====================================================================================================================================\n");
248  fprintf(mLogFile,"L2gammaAlgo start of run %d summary, compiled: %s , %s\n",run,__DATE__,__TIME__);
249  fprintf(mLogFile,"calorimeter: %s\n",mName);
250  fprintf(mLogFile,"run: %d\n",run);
251  fprintf(mLogFile,"use offline gains I[0]: %d\n",mUseOfflineGains);
252  fprintf(mLogFile,"prescaled accept I[1]: %d\n",mPrescale);
253  fprintf(mLogFile,"threshold level: %d\n",mThresholdLevel);
254  fprintf(mLogFile,"tower threshold F[0]: %-5.2f [GeV]\n",mTowerThreshold);
255  fprintf(mLogFile,"patch threshold F[1]: %-5.2f [GeV]\n",mPatchThreshold);
256  fprintf(mLogFile,"patch size: %s\n","3x3"); /* may change in later revisions */
257  fprintf(mLogFile,"logfile: %s\n",clog);
258  // fprintf(mLogFile,"histograms: %s\n",chis);
259  fprintf(mLogFile,"database at: %p\n",mDb);
260  fprintf(mLogFile,"l2 input at: %p\n",mL2input);
261 
262  fprintf(mLogFile,"\nDetector Geometry\n");
263  fprintf(mLogFile,"mEtaBins = {");
264  for ( int i=0;i<mNumEtas+1;i++ ) {
265  fprintf(mLogFile,"%4.2f, ",mEtaBins[i]);
266  if ( !((i+1)%4) ) fprintf(mLogFile,"\n ");
267  }
268  fprintf(mLogFile,"X}\n");
269  fprintf(mLogFile,"mMaxADC = %4d\n",mMaxADC);
270  fprintf(mLogFile,"mMaxET = %4.1f\n",mMaxET);
271  fprintf(mLogFile,"mIdealGainT = %5.3f\n",mIdealGainT);
272  fprintf(mLogFile,"\n");
273 
274  fprintf(mLogFile,"\n");
275  }
276  if ( !mDb ) return 100;
277 
278 
279 
281  mNumberInput = 0;
282  mNumberAcceptHT = 0;
283  mNumberAcceptTP = 0;
284 
285  mNumberLive=0;
286 
288  if(mLogFile) fprintf(mLogFile,"initialize EMC ascii database\n\n");
289  //$$$ assume it's already initialized by a call in main??? gEmcCDb_init();
290  if(mLogFile) fprintf(mLogFile,"\n");
291 
295  //float emc_ideal_gains[mNumEtas];
296  float *emc_ideal_gains = new float[mNumEtas];
297  if(mLogFile) fprintf(mLogFile,"configure ideal gains for %d eta bins\n",mNumEtas);
298  for ( int i=0;i<mNumEtas;i++ )
299  {
300  float eta_mean = mEtaBins[ i ] + mEtaBins[ i+1 ];
301  eta_mean /= 2.0;
302  emc_ideal_gains[i] = mMaxADC / mMaxET / cosh(eta_mean);
303  if(mLogFile) fprintf(mLogFile,"+ etabin=%d eta=%5.2f bemc ideal gain=%5.2f\n", i+1, eta_mean, emc_ideal_gains[i]);
304 
305  }
306 
307  if ( mNumEtas != 12 && mNumEtas != 40 ) {
308  if(mLogFile) fprintf(mLogFile,"L2gammaEmCal::FATAL ERROR -- expect 12 or 40 etabins, got %d\n",mNumEtas);
309  if(mLogFile) fprintf(mLogFile," + likely cause -- invalid runtime parameter\n");
310  delete emc_ideal_gains;
311  return 10; /* EMC not properly configured */
312  }
313 
315  for ( int index=0;index<mNumTower;index++ )
316  {
317  mTowerAdcThreshold[index]=0;
318  mPatchAdcThreshold[index]=mPatchThreshold;
319  mTowerPed[index]=0.;
320  mPatchPed[index]=0.;
321  mTowerFrequency[index]=0;
322  mPatchFrequency[index]=0;
323  mTowerGain[index]=-1.0;
324  }
325 
326  if(mLogFile) fprintf(mLogFile,"\nInitializing tower ADC threshods\n");
327 
329  for ( int index=0; index<EmcDbIndexMax; index++ )
330  {
331 
332  // get the item from the database
333  //$$$struct EmcCDbItem *x = &gEmcCDbByIndex[index];
334  const L2EmcDb::EmcCDbItem *x = mDb->getByIndex(index);
335  if ( x==0 ) continue;
336 
337 #if 1
338  // make sure db item matches the calorimeter we're using
339  if ( mBEmc && !mDb->isBTOW(x) )
340  continue;
341  if ( mEEmc && !mDb->isETOW(x) )
342  continue;
343 #endif
344 
345 
346  /* temporary kludge */
347  //$$$ if ( x->gain < 0 ) x->fail |= 0x0080;
348 
349  // local numbering scheme counts from 0
350  int sec = x->sec - 1;
351  int sub = 8192;
352  if ( mEEmc ) sub = x->sub - 'A';
353  if ( mBEmc ) sub = x->sub - 'a';
354  int eta = x->eta - 1;
355  int phi = phibin(sec,sub); //mNumSubs*sec + sub;
356  int tow = tower(phi,eta); //mNumEtas*phi + eta;
357  int rdo = x->rdo;
358 
359  //fprintf(mLogFile,"sec=%d sub=%d eta=%d tow=%d\n",sec,sub,eta,tow);
360 
361  /* temporary kludge */
362  //if ( eta < 20 ) x->fail |= 0x0080;
363 
364 
365  // equal energy gain
366  //float gain = x->gain;
367  float ped = x->ped;
368  float gain = x->gain;
369  float ideal = emc_ideal_gains[eta];
370 
371 
372 
373  // check fail bit and artificially raise HT threshold
374  unsigned short stat = x->stat;
375  unsigned short fail = x->fail;
376  if ( gain < 0.0 ) {
377  if(mLogFile) fprintf(mLogFile,"L2gammaEmCal::WARN %s gain=%5.2f will be masked\n",x->name,x->gain);
378  fail = 0xffff;
379  }
380 
381  mTowerStat[rdo]=stat; mPatchStat[rdo]=stat;
382  mTowerFail[rdo]=fail; mPatchFail[rdo]=fail;
383  mTowerGain[rdo]=gain;
384  mTowerGainIdeal[rdo]=ideal;
385 
386 
387  if ( !fail )
388  mTowerAdcCorrection[rdo]=(mUseOfflineGains)?ideal/gain:1.0;
389  else
390  mTowerAdcCorrection[rdo]=0.; /* tower is masked out */
391 
392  // calculate high tower threshold for this
393  // rdo channel in ADC
394  if ( !fail ) {
395  mTowerAdcThreshold[rdo] = (unsigned short)(ped + (mIdealGainT * mTowerThreshold)/mTowerAdcCorrection[rdo] + 0.5 );
396  mTowerPed[rdo]=ped;
397  }
398  else {
399  mTowerAdcThreshold[rdo] = 4095;
400  mTowerPed[rdo]=4095;
401  }
402 
403  /* upgrade: verify that tower threshold is above mMinAdc, issue warning and log */
404 
405  mRdo2tower[ rdo ] = tow; // returns tower for given rdo channel
406  mTower2rdo[ tow ] = rdo; // returns rdo channel for given tower
407 
408  if ( !fail ) mNumberLive++;
409 
410  }
411 
412 
413 
414 
415 
416  // fprintf(mLogFile,"L2gammaAlgo::initRun(%d) initialize 3x3 patch thresholds for ET=%5.2f\n",mRunNumber,mPatchThreshold);
417 
418  if(mLogFile) fprintf(mLogFile,"Initializing 3x3 tower cluster thresholds\n");
419 
422  for ( int index=0; index<EmcDbIndexMax; index++ )
423  {
424 
425  // get the item from the database
426  //$$$ struct EmcCDbItem *x = &gEmcCDbByIndex[index];
427  const L2EmcDb::EmcCDbItem *x = mDb->getByIndex(index);
428  if ( x==0 ) continue;
429 
430  // make sure db item matches the calorimeter we're using
431  if ( mBEmc && !mDb->isBTOW(x) )
432  continue;
433  if ( mEEmc && !mDb->isETOW(x) )
434  continue;
435 
436  // local numbering scheme counts from 0
437  int sec = x->sec - 1;
438  int sub = 8192;
439  if ( mEEmc ) sub = x->sub - 'A';
440  if ( mBEmc ) sub = x->sub - 'a';
441  int eta = x->eta - 1;
442  int phi = phibin(sec,sub);
443  //int tow = tower(phi,eta);
444  int rdo = x->rdo;
445 
446  mPatchStat[rdo] |= x->stat;
447  mPatchFail[rdo] |= x->fail;
448 
449  // masked towers will be ignored in the patch thresholds...
450  // our butts are covered in the event that a crate is masked
451  // out, by the requirement of a high tower in coincidence with
452  // a 3x3 tower patch sum.
453  if ( x->fail ) continue;
454  if ( mTowerFail[rdo] ) continue;
455 
456  // equal energy gain
457  //float gain = x->gain;
458  float ped = x->ped;
459 
460 
465  int nn=0;
466  for ( int ieta=eta-1;ieta<=eta+1;ieta++ )
467  for ( int iphi=phi-1;iphi<=phi+1;iphi++ )
468  {
469 
470  if ( ieta < 0 || ieta > mNumEtas-1 ) continue; // out of bounds
471  int jeta=ieta;
472  int jphi=(iphi+mNumPhis)%mNumPhis;
473  int jtow=tower(jphi,jeta);
474  int jrdo=mTower2rdo[jtow];
475 
476  mPatchAdcThreshold[jrdo] += ped/mIdealGainT *
477  mTowerAdcCorrection[rdo];
478 
479  if ( jtow == 12 ) {
480  if(mLogFile) fprintf(mLogFile,"+ name=%s ped=%5.2f gfact=%5.2f thresh=%5.2f\n",x->name,x->ped,mTowerAdcCorrection[rdo],mPatchAdcThreshold[jrdo] );
481  }
482 
483  // add ped for this channel to all patches
484  // which contain it
485  mPatchPed[jrdo]+=ped;
486 
488  mRdoPatch[rdo][nn++]=jrdo;
489 
490 
491  }
492  mNumPatch[rdo]=nn;
493 
494  }
495 
496  //fprintf(mLogFile,"done w/ initialization, writing to log file\n");
497 
498 
499  //
500  // Write to logfile (stdout for now) the thresholds generated
501  //
502 
503  for ( int index=0;index<EmcDbIndexMax;index++ )
504  {
505 
506 
507  // get the item from the database
508  //. struct EmcCDbItem *x = &gEmcCDbByIndex[index];
509  const L2EmcDb::EmcCDbItem *x = mDb->getByIndex(index);
510  if ( x==0 ) continue;
511 
512 #if 1
513  // make sure db item matches the calorimeter we're using
514  if ( mBEmc && !mDb->isBTOW(x) )
515  continue;
516  if ( mEEmc && !mDb->isETOW(x) )
517  continue;
518 #endif
519 
520  float gain=x->gain;
521  int rdo=x->rdo;
522  if ( gain==0. ) gain=-1.;
523 
524  int tow=mRdo2tower[rdo];
525  int eta=tow%mNumEtas;
526 
527 
528  /*
529  Floating point exception occurs in my output????
530  */
531 
532 #if 0
533  fprintf(mLogFile,"tower=%s tow=%d rdo=%d mask=%2x gain=%5.2f ideal=%5.2f ped=%5.2f thr=%d %5.2f GeV patch=%5.2f + %5.2f GeV\n",
534  x->name,
535  mRdo2tower[rdo],
536  rdo,
537  x->fail,
538  x->gain,
539  emc_ideal_gains[eta],
540  x->ped,
541  mTowerAdcThreshold[rdo],
542  (float)((mTowerAdcThreshold[rdo]-x->ped)*mMaxET/mMaxADC),
543  mPatchAdcThreshold[rdo]-mPatchThreshold,
544  mPatchThreshold);
545 #endif
546 
547 #if 1
548  if(mLogFile) fprintf(mLogFile,"tower=%s tow=%d rdo=%d ideal=%5.2f thr=%d patch=%5.2f GeV\n",
549  x->name,
550  mRdo2tower[rdo],
551  rdo,
552  emc_ideal_gains[eta],
553  mTowerAdcThreshold[rdo],
554  mPatchAdcThreshold[rdo]);
555 #endif
556 
557 
558 
559  }
560 
561 #ifdef DEBUG_PATCH_THRESHOLDS
562  for ( int rdo=0;rdo<mNumRdo;rdo++ )
563  printPatchConfig(rdo);
564 #endif
565 
566  if(mLogFile) {
567  fprintf(mLogFile,"\n");
568  fprintf(mLogFile,"total number of towers: %d\n",mNumTower);
569  fprintf(mLogFile,"number of unmasked towers: %d\n",mNumberLive);
570  fprintf(mLogFile,"\n====================================================================================================================================\n\n");
571  }
572 
573  /*
574  * cleanup allocated memory
575  */
576  delete emc_ideal_gains;
577 
578  return 0;
579 
580 }/* ::initRun */
581 
582 
583 
584 // -------------------------------------------------------------------------------------------- eval --
585 //
586 // evaluate trigger conditions
587 //
588 // most of the work has already been done for us. at the start of
589 // the run, we took the user-supplied ET thresholds on the high
590 // tower and associated cluster (patch) and determined two thresholds:
591 //
592 // o mTowerAdcThreshold[rdo] -- the ADC value for which the tower will
593 // exceed the specified ET threshold (i.e. ped*4096/60+E_T).
594 //
595 // o mPatchAdcThreshold[rdo] -- actually a floating point ET threshold
596 // to be calculated ONLY when we have a tower above the specified
597 // ET threshold.
598 //
599 // therefore, the only thing we need to do w/in eval is to loop over
600 // all rdo's and find a high tower, then test the 3x3 tower cluster
601 // against its threshold.
602 //
603 
604 bool L2gammaAlgo::doEvent( int L0trigger, int inputEventID, TrgDataType* trgData,
605  int bemcIn, unsigned short *bemcData,
606  int eemcIn, unsigned short *eemcData )
607 {
608  rdtscl_macro(mEveTimeStart);
609  mAccept=false;
610  if ( mEEmc ) {
611  mAccept=doEvent( inputEventID, trgData, eemcIn, eemcData );
612  }
613  else {
614  mAccept=doEvent( inputEventID, trgData, bemcIn, bemcData );
615  }
616 
617  rdtscl_macro(mEveTimeStop);
618  mEveTimeDiff=mEveTimeStop-mEveTimeStart;
619  int kTick=mEveTimeDiff/1000;
620  // printf("gg=%f t1=%d t2=%d \n",mEveTimeDiff/1000.,mEveTimeStart,mEveTimeStop);
621  mhT->fill(kTick);
622 
623  return mAccept;
624 }
625 
626 
627 bool L2gammaAlgo::doEvent( int inpEveId, TrgDataType* trgData,
628  int emcIn, unsigned short *emcData )
629 {
630  //I think you will need to use inpEveId for second copy, if not - ignore it,JB
631  clear();
632  // we got nothin'
633  if ( !emcIn ) return false;
634 
635  mNumberInput++;
636 
637  mHistos[0].fill(0);
638 
639  //
640  // first section of the code evaluates the trigger decision
641  //
642 
643  // track execution time of the algorithm (in cpu ticks)
644  unsigned long eval_time_start;
645  unsigned long eval_time_stop;
646  rdtscl_macro(eval_time_start);
647 
648  bool eht = false;
649  bool trig = false;
650 
651  float sum=0.; /* tower patch sum, 3x3 cluster around high tower */
652  unsigned short htrdo=8192; /* rdo of high tower */
653  unsigned short tprdo=8192; /* rdo of high cluster */
654 
655  nRdosHT=0;
656  nRdosTP=0;
657 
658  // ------------------------------------------------------------------------------------ preprocess --
659 #ifdef OPT_PREPROCESS
660 
661  // If another "lower threshold" is being run as input to this algorithm,
662  // then we process the results of that algorithm and jump to the QA
663  // section.
664 
665  if ( mL2input ) {
666  for ( unsigned short i=0;i<mL2input->getNumberOfHighTowers();i++ )
667  {
668  unsigned short rdo=mL2input->mListOfRdosHT[i];
669 #ifdef DEBUG
670  mHistos[3].fill( emcData[rdo] );
671 #endif
672  if ( emcData[rdo] < mTowerAdcThreshold[rdo]) continue;
673  mADCofHT[ nRdosHT ] = emcData[rdo];
674  mListOfRdosHT[ nRdosHT++ ] = rdo; /* store rdo of high tower */
675  mHistos[4].fill( emcData[rdo] );
676 
677  eht=true;
678  htrdo=rdo;
679  sum = 0.;
680  for ( int i=0;i<mNumPatch[rdo];i++ ) {
681  int prdo=mRdoPatch[rdo][i];
682  if ( mTowerFail[prdo] ) continue; /* skip masked towers */
683  sum += emcData[prdo] / mIdealGainT * mTowerAdcCorrection[rdo];
684  }
685  unsigned short sumx2=(unsigned short)(2.0*sum);
686 #ifdef DEBUG
687  mHistos[7].fill( sumx2 );
688 #endif
689 
690  if ( sum > mPatchAdcThreshold[rdo] ) {
691  trig=true;
692  tprdo=rdo;
693  mETofTP[ nRdosTP ] = sum-mPatchAdcThreshold[rdo]+mPatchThreshold;
694  mListOfRdosTP[ nRdosTP++ ] = rdo; /* store rdo of cluster */
695  mHistos[5].fill( emcData[rdo] );
696  mHistos[8].fill( sumx2 );
697  }
698  }
699  goto QA;
700  }
701 #endif
702 
703  // ------------------------------------------------------------------------- standalone processing --
704 
705  // Here we loop over all RDO channels and determine whether thresholds have
706  // been met. We store results for higher thresholds/other algorithms to
707  // process at a later point.
708 
709  for ( unsigned short rdo=0; rdo<mNumRdo; rdo++ )
710  {
711 
712  // check if tower is above ADC=ET*gain+ped threshold
713  // these two lines consume about 1/2 of the time of the
714  // algorithm for any given event.
715 #ifdef DEBUG
716  mHistos[3].fill( emcData[rdo] );
717 #endif
718  if ( emcData[rdo] < mTowerAdcThreshold[rdo]) continue;
719 
720  // past here, we have a high tower trigger. when we
721  // have such a trigger, the following code takes about
722  // the same time as the preceding two lines.
723 
724  mADCofHT[ nRdosHT ] = emcData[rdo];
725  mListOfRdosHT[ nRdosHT++ ] = rdo; /* store rdo of high tower for 2nd stage */
726  mHistos[4].fill( emcData[rdo] );
727 
728  eht=true;
729  htrdo=rdo;
730 
731 
732  // Sum the 3x3 patch of towers around the high tower
733  sum = 0.;
734  for ( int i=0;i<mNumPatch[rdo];i++ ) {
735  int prdo=mRdoPatch[rdo][i];
736  if ( mTowerFail[prdo] ) continue; /* skip masked towers */
737  sum += emcData[prdo] / mIdealGainT * mTowerAdcCorrection[rdo];
738  }
739  unsigned short sumx2=(unsigned short)(2.0*sum);
740 #ifdef DEBUG
741  mHistos[7].fill( sumx2 );
742 #endif
743 
744  // we have found a high tower in coincidence
745  // with a 3x3 patch of towers above specified
746  // thresholds
747  if ( sum > mPatchAdcThreshold[rdo] ) {
748  trig=true;
749  tprdo=rdo;
750  mETofTP[ nRdosTP ] = sum-mPatchAdcThreshold[rdo]+mPatchThreshold;
751  mListOfRdosTP[ nRdosTP++ ] = rdo; /* store rdo of cluster */
752  mHistos[5].fill( emcData[rdo] );
753  mHistos[8].fill( sumx2 );
754  }
755 
756  }
757 
758  // ------------------------------------------------------------------------ process QA histograms --
759 
760  // Fill QA histograms, set bits in the L2 result, and return the result
761 
762 #ifdef OPT_PREPROCESS
763  QA:
764 #endif
765 
766  bool prescaleAccept = false;
767  if ( mPrescale>0 )
768  prescaleAccept = !((mNumberInput) % mPrescale);
769 
770  mHistos[1].fill(nRdosHT); // Number of high towers > threshold
771  mHistos[2].fill(nRdosTP); // Number of clusters > threshold
772 
773 
774  if ( eht ) {
775  mNumberAcceptHT++;
776  mTowerFrequency[mRdo2tower[htrdo]]++;
777  mHistos[0].fill(1);
778  }
779 
780  if ( trig ) {
781  mNumberAcceptTP++;
782  mPatchFrequency[mRdo2tower[tprdo]]++;
783  mHistos[0].fill(2);
784  }
785  rdtscl_macro(eval_time_stop);
786 
787  /* algorithm selects fraction of events to do QA on */
788  unsigned long qa_time_start;
789  unsigned long qa_time_stop;
790  rdtscl_macro(qa_time_start);
791 
792  // ------------------------------------------------------------------------- save trigger decision --
793 
794 
795  L2gammaResult mResult;
796  // initialize "result"
797  mResult.result_version = LEVEL2_GAMMA_RESULT_VERSION;
798  mResult.threshold = (mThresholdLevel==kThresh1)? (0x1|0x2) : (0x4|0x8);
799  // reset monitor portion of L2gammaResult
800  mResult.elapsed=0x00;
801  // clear L2 trigger result
802  mResult.trigger = 0x00;
803 
804  float max=0.;
805  unsigned short maxrdo=8192;
806 
807  // set high tower bit
808  if ( eht )
809  mResult.trigger |= 0x1;
810 
811 
812  // set cluster bit and
813  for ( unsigned short i=0;i<nRdosTP;i++ )
814  {
815  if ( mETofTP[i] > max ) {
816  max=mETofTP[i];
817  maxrdo=mListOfRdosTP[i];
818  }
819  }
820  if ( maxrdo < mNumRdo )
821  {
822  mResult.trigger |= 0x2;
823  unsigned short tow=mRdo2tower[maxrdo];
824  mResult.phibin=tow/mNumEtas;
825  mResult.etabin=tow%mNumEtas;
826  max-=mPatchAdcThreshold[maxrdo];
827  max+=mPatchThreshold;
828  if ( max < 127.5 )
829  mResult.ptclusterx2 = ((unsigned short)(2.0*max));
830  else
831  mResult.ptclusterx2 = 255;
832  float pttow=emcData[maxrdo]-mTowerPed[maxrdo];
833  pttow*=mTowerAdcCorrection[maxrdo]/mIdealGainT;
834  if ( pttow < 127.5 )
835  mResult.pttowerx2 = (unsigned short)(2.0*pttow);
836  else
837  mResult.pttowerx2 = 255;
838  if ( mEEmc ) mResult.phibin |= 0x8;
839  }
840 
841 
842  // set prescale bit
843  if ( prescaleAccept )
844  mResult.trigger |= 0x4;
845 
846  // set trigger bit
847  if ( trig || prescaleAccept )
848  mResult.trigger |= 0x8;
849 
850  rdtscl_macro(qa_time_stop);
851  int dt=(int)(qa_time_stop-eval_time_start)/1000;
852 
853  mEvalTime += dt;
854  //mResult.mon.elapsed = 0;
855 
859  unsigned int *pResult
860 // = &trgData->TrgSum.L2Result[L2RESULTS_OFFSET_PIG +2*mEEmc+mThresholdLevel];
861  = &trgData->TrgSum.L2Result[mResultOffset+2*mEEmc+mThresholdLevel];
862 
866  memcpy(pResult,&mResult,sizeof(L2gammaResult));
867  //#define __L2GAMMA_PRINT_RESULT__
868 #ifdef __L2GAMMA_PRINT_RESULT__
869  print_L2gammaResult(mResult);
870 #endif
871 
877  if ( trig || prescaleAccept )
878  {
879  mHistos[0].fill(3);
880  for ( unsigned short i=0;i<nRdosHT;i++ ) {
881  int tower=(int)mRdo2tower[mListOfRdosHT[i]];
882  mHistos[6].fill( (mADCofHT[i] - (unsigned short)mTowerPed[ mListOfRdosHT[i] ])/5,tower );
883  mHistos[13].fill(tower);
884  }
885 
886  for ( unsigned short i=0;i<nRdosTP;i++ ) {
887  int tower=(int)mRdo2tower[mListOfRdosTP[i]];
888  mHistos[9].fill( (int)(mETofTP[i]*2),tower );
889  mHistos[14].fill( tower );
890  }
891 
892  mHistos[18].fill( dt );
893 
894  }
895 
896  mHistos[15].fill(dt);
897  if ( eht ) mHistos[16].fill(dt);
898  if ( trig ) mHistos[17].fill(dt);
899 
900  return ( trig || prescaleAccept );
901 
902 }
903 
904 // ----------------------------------------------------------------------------
905 //
906 
907 void L2gammaAlgo::clear()
908 {
909 
910  nRdosHT=0;
911  nRdosTP=0;
912 
913 }
914 
915 // ----------------------------------------------------------------------------
916 
917 void L2gammaAlgo::finishRun()
918 {
919 
920  int run=mRunNumber;
921 
922 
923  char chis[128];
924  printf("%s/run%d.l2%s.hist.bin\n",mOutDir,run,mName);
925  sprintf(chis,"%s/run%d.l2%s.hist.bin",mOutDir,run,mName);
926  printf("L2gamma:saving '%s'\n",chis);
927  // mLogFile=fopen(clog,"a");
928  // mLogFile=stdout;
929 
930  mHistFile=fopen(chis,"w");
931  if(mHistFile==0) printf("%s ::finishRun() for run=%d UNABLE to open histo-file=%s=, continue initialization\n",mName,run,chis);
932  //assert( mHistFile);
933 
934  //
935  // trap some potential divide by zeros but still do output
936  //
937  if ( !mNumberInput ) mNumberInput=-1; /* short/nonexistent run */
938  if ( !mNumberLive ) mNumberLive=-1; /* db is ft^ */
939 
940 
941  // if ( !mLogFile )
942  // mLogFile = stdout;
943  if(mLogFile) {
944  fprintf(mLogFile,"\n\n===================================================================================================================================\n");
945  fprintf(mLogFile,"L2gammaAlgo end of run %d summary\n",run);
946  fprintf(mLogFile,"run: %d\n",run);
947  fprintf(mLogFile,"tower threshold: %5.2f [GeV]\n",mTowerThreshold);
948  fprintf(mLogFile,"patch threshold: %5.2f [GeV]\n",mPatchThreshold);
949  fprintf(mLogFile,"patch size: 3x3\n"); /* may change in later revisions */
950  fprintf(mLogFile,"eval time: %d [kTicks]\n",mEvalTime);
951  fprintf(mLogFile,"avg time/event: %d [kTicks]\n",mEvalTime/mNumberInput);
952  fprintf(mLogFile,"# input: %d\n",mNumberInput);
953  fprintf(mLogFile,"# ht accept: %d\n",mNumberAcceptHT);
954  fprintf(mLogFile,"# ht+tp accept: %d\n",mNumberAcceptTP);
955  fprintf(mLogFile,"\n");
956 
957  //
958  // first section, loop over the tower and patch frequencies
959  // to look for hot towers
960  //
961  float avgt=0,avgp=0;
962  int maxtf=0,maxpf=0;
963 
964  for ( int tow=0;tow<mNumTower;tow++ )
965  {
966  avgt+=mTowerFrequency[tow];
967  avgp+=mPatchFrequency[tow];
968  if ( mTowerFrequency[tow]>maxtf ) maxtf=mTowerFrequency[tow];
969  if ( mPatchFrequency[tow]>maxpf ) maxpf=mTowerFrequency[tow];
970  }
971  avgt=( ((float)avgt)/(float)mNumberLive );
972  avgp=( ((float)avgp)/(float)mNumberLive );
973 
974  // fprintf(mLogFile,"sumt=%5.2f sump=%5.2f nlive=%d\n",avgt,avgp,mNumberLive);
975  fprintf(mLogFile,"max # ht accept in one tower: %d\n",maxtf);
976  fprintf(mLogFile,"max # ht+tp accept in one tower: %d\n",maxpf);
977  fprintf(mLogFile,"avg # ht accept: %6f\n",avgt);
978  fprintf(mLogFile,"avg # ht+tp accept: %6f\n",avgp);
979  fprintf(mLogFile,"expected # ht accept in one tower: %6f\n",((float)mNumberAcceptHT)/mNumberLive);
980  fprintf(mLogFile,"expected # ht+tp accept in one tower: %6f\n",((float)mNumberAcceptTP)/mNumberLive);
981  fprintf(mLogFile,"\n");
982 
983 
984  //#define DEBUG_PATCH_THRESHOLDS
985 
986 #ifndef DEBUG_PATCH_THRESHOLDS
987  fprintf(mLogFile,"Summary of hot towers (# ht accept > 3*average or # ht+tp accept > 3*average):\n\n");
988 
989  for ( int tow=0;tow<mNumTower;tow++ )
990  {
991  int phi=tow/mNumEtas;
992  int eta=tow%mNumEtas;
993  int sec=phi/mNumSubs;
994  int sub=phi%mNumSubs;
995  int myrdo=mTower2rdo[tow];
996 
997 
998  bool hott = (mTowerFrequency[tow] > (3.0*avgt));
999  bool hotp = 0 && (mPatchFrequency[tow] > (3.0*avgp));
1000 
1001  /* ^^^^^^^^^^^^^ need a statistically robust sample or we will get noise! */
1002 
1003  if ( hott||hotp ) fprintf(mLogFile,"\n");
1004 
1005  if ( mBEmc )
1006  {
1007  if ( hott )
1008  fprintf(mLogFile,"L2gammaAlgo::WARNING -- possible hot tower %02dt%c%02d freq=%d avg=%5.2f\n",sec+1,sub+'a',eta+1,mTowerFrequency[tow],avgt);
1009  if ( hotp )
1010  fprintf(mLogFile,"L2gammaAlgo::WARNING -- possible hot patch %02dt%c%02d freq=%d avg=%5.2f\n",sec+1,sub+'a',eta+1,mPatchFrequency[tow],avgp);
1011  }
1012  else
1013  {
1014  if ( hott )
1015  fprintf(mLogFile,"L2gammaAlgo::WARNING -- possible hot tower %02dT%c%02d freq=%d avg=%5.2f\n",sec+1,sub+'A',eta+1,mTowerFrequency[tow],avgt);
1016  if ( hotp )
1017  fprintf(mLogFile,"L2gammaAlgo::WARNING -- possible hot patch %02dT%c%02d freq=%d avg=%5.2f\n",sec+1,sub+'A',eta+1,mPatchFrequency[tow],avgp);
1018  }
1019 
1020  if ( hott || hotp )
1021  {
1022 
1023  printPatchConfig( myrdo );
1024 
1025  }
1026 
1027  }
1028  }
1029 #endif
1030 
1031  if ( mLogFile ) {
1032  fprintf(mLogFile,"\n====================================================================================================================================\n\n");
1033  finishCommonHistos();
1034 
1035  }
1036 
1037  finish(); /* output histogram */
1038 
1039 }
1040 
1041 void L2gammaAlgo::finish()
1042 {
1043  // Write histograms at end of run
1044  if ( mHistFile ) {
1045  for ( unsigned int i=0;i<sizeof(mHistos)/sizeof(L2Histo);i++)
1046  {
1047  mHistos[i].write(mHistFile);
1048  }
1049  fflush(mHistFile);
1050  }
1051 }
1052 
1053 // ----------------------------------------------------------------------------
1054 
1055 void L2gammaAlgo::jclear()
1056 {
1057  for ( int i=0;i<19;i++ ) mHistos[i].reset();
1058 }
1059 void L2gammaAlgo::jbook()
1060 {
1061  // Book L2 histograms for QA
1062 
1063  // on first call book histograms
1064  if ( 1 ) {
1065 
1066  /* Event counters */
1067 
1068  mHistos[0]=L2Histo(100, "Counter. 0=N_{input} 1=N_{ht} 2=N_{clust} 3=N_{accept}", 10);
1069  mHistos[1]=L2Histo(101, "N high towers > threshold", 10);
1070  mHistos[2]=L2Histo(102, "N clusters > threshold", 10);
1071 
1072  /* ADC spectra */
1073 
1074  mHistos[3]=L2Histo(103, "Raw ADC of ht (if DEBUG)", 1024);
1075  mHistos[4]=L2Histo(104, "Raw ADC of ht, ht > threshold", 1024);
1076  mHistos[5]=L2Histo(105, "Raw ADC of ht, cluster > threshold", 1024);
1077  mHistos[6]=L2Histo(106, "HT event accepted; X: (ADC-ped)/5; Y: tower index", 128,mNumTower);
1078 
1079  /* Cluster spectra */
1080 
1081  mHistos[7]=L2Histo(107, "2*ET sum of cluster, ht > threshold", 32);
1082  mHistos[8]=L2Histo(108, "2*ET sum of cluster, cluster > threshold", 32);
1083  mHistos[9]=L2Histo(109, "2*ET sum of cluster - ped/gain, event accepted", 32,mNumTower);
1084 
1085  /* BBC difference (not implemented) */
1086 
1087  mHistos[10]=L2Histo(110, "Raw BBC TAC difference+256 (not implemented)", 1);
1088  mHistos[11]=L2Histo(111, "Raw BBC TAC difference+256, ht > threshold (not implemented)", 1);
1089  mHistos[12]=L2Histo(112, "Raw BBC TAC difference+256, cluster > threshold (not implemented)" , 1);
1090 
1091  /* Tower and cluster trigger frequencies */
1092 
1093  mHistos[13]=L2Histo(113, "High tower > threshold frequency, X: ieta+mNumEta*iphi", mNumTower);
1094  mHistos[14]=L2Histo(114, "Cluster > threshold frequency, X: ieta+mNumEta*iphi", mNumTower);
1095 
1096  /* Event timing */
1097 
1098  mHistos[15]=L2Histo(115, "kTicks / input event", 128);
1099  mHistos[16]=L2Histo(116, "kTicks / ht event", 128);
1100  mHistos[17]=L2Histo(117, "kTicks / cluster event", 128);
1101  mHistos[18]=L2Histo(118, "kTicks / trigger", 128);
1102 
1103  }
1104 
1105 
1106 }
1107 
1108 // ----------------------------------------------------------------------------
1109 #ifndef DEBUG_PATCH_THRESHOLDS
1110 void L2gammaAlgo::printPatchConfig( int rdo )
1111 {
1112 
1113  // get the tower id
1114  int tow = mRdo2tower[rdo];
1115  int phi = tow/mNumEtas;
1116  int eta = tow%mNumEtas;
1117 
1118  int sec = phi/mNumSubs;
1119  int sub = phi%mNumSubs;
1120 
1121  /*
1122 
1123  The goal here is to output several quantities in an ascii
1124  format which matches the physical layout of the towers in
1125  order to help us debug problems online
1126 
1127  as least as far as the endcap is concerned (I can switch the
1128  order for the barrel if need be. Specifically, we want the
1129  following output:
1130 
1131  tower: 06TC05
1132  gain: 22.64
1133  ped: 20.90
1134  thr: 212 [2.8 GeV]
1135 
1136  pedestals gains, stat, fail, high tower frequency, etc...
1137 
1138  18.7 | -3.0 | 22.5
1139  -----+-----+-----
1140  22.9 | 20.9 | 15.9
1141  -----+------+-----
1142  11.1 | 16.7 | 32.0
1143 
1144  */
1145 
1146  if ( mEEmc )
1147  fprintf(mLogFile,"etow: %02dT%c%02d\n",sec+1,sub+'A',eta+1);
1148  if ( mBEmc )
1149  fprintf(mLogFile,"btow: %02dt%c%02d\n",sec+1,sub+'a',eta+1);
1150 
1151  fprintf(mLogFile,"index: %4d\n",tow);
1152  fprintf(mLogFile,"gain: %6.3f [ideal=%6.3f]\n",mTowerGain[rdo],mTowerGainIdeal[rdo]);
1153  fprintf(mLogFile,"ped: %6.3f\n",mTowerPed[rdo]);
1154  fprintf(mLogFile,"thr: %3d [%4.2f GeV]\n",mTowerAdcThreshold[rdo],mTowerAdcThreshold[rdo]/mIdealGainT);
1155  fprintf(mLogFile,"\n");
1156 
1157  // 14.3 | 14.5 | 8.8 22.6 | 22.9 | 21.6 80 | 80 | 80 0 | 0 | 0 1 | 0 | 0
1158  fprintf(mLogFile,"loaded peds gain factors status bits mask bits frequency > thr\n");
1159 
1160  for ( int i_eta=eta+1;i_eta>=eta-1;i_eta-- ) {
1161 
1162  int phi_l = (phi+mNumPhis-1)%mNumPhis;
1163  int phi_r = (phi+mNumPhis+1)%mNumPhis;
1164  int phi_c = phi;
1165 
1166  //int phis[]={phi_l,phi,phi_r};
1167  float peds[]={0,0,0};
1168  float gains[]={0,0,0};
1169  unsigned short stats[]={0xff,0xff,0xff};
1170  unsigned short fails[]={0xff,0xff,0xff};
1171  int freqs[]={0,0,0};
1172 
1173  // get numbers from db
1174  if ( i_eta < mNumEtas && i_eta >= 0 ) {
1175 
1176  int tow_l = tower(phi_l,i_eta);
1177  int tow_r = tower(phi_r,i_eta);
1178  int tow_c = tower(phi_c,i_eta);
1179 
1180  int rdo_l = mTower2rdo[tow_l];
1181  int rdo_r = mTower2rdo[tow_r];
1182  int rdo_c = mTower2rdo[tow_c];
1183 
1184  peds[0]=mTowerPed[rdo_l];peds[1]=mTowerPed[rdo_c];peds[2]=mTowerPed[rdo_r];
1185  gains[0]=mTowerAdcCorrection[rdo_l];gains[1]=mTowerAdcCorrection[rdo_c];gains[2]=mTowerAdcCorrection[rdo_r];
1186  stats[0]=mTowerStat[rdo_l];stats[1]=mTowerStat[rdo_c];stats[2]=mTowerStat[rdo_r];
1187  fails[0]=mTowerFail[rdo_l];fails[1]=mTowerFail[rdo_c];fails[2]=mTowerFail[rdo_r];
1188  freqs[0]=mTowerFrequency[tow_l];freqs[1]=mTowerFrequency[tow_c];freqs[2]=mTowerFrequency[tow_r];
1189 
1190 
1191  }
1192 
1193  fprintf(mLogFile, "%4.1f | %4.1f | %4.1f\t", peds[0],peds[1],peds[2] );
1194  fprintf(mLogFile, "%4.2f | %4.2f | %4.2f\t", gains[0],gains[1],gains[2] );
1195  fprintf(mLogFile, "%2x | %2x | %2x \t", stats[0],stats[1],stats[2] );
1196  fprintf(mLogFile, "%2x | %2x | %2x \t", fails[0],fails[1],fails[2] );
1197  fprintf(mLogFile, "%4.0f | %4.0f | %4.0f\t", (float)freqs[0],(float)freqs[1],(float)freqs[2] );
1198 
1199  fprintf(mLogFile,"\n");
1200 
1201  if ( i_eta > eta-1 ) {
1202  fprintf(mLogFile, "-----+------+-----\t");
1203  fprintf(mLogFile, "-----+------+-----\t");
1204  fprintf(mLogFile, "-----+------+-----\t");
1205  fprintf(mLogFile, "-----+------+-----\t");
1206  fprintf(mLogFile, "-----+------+-----\t");
1207  fprintf(mLogFile, "\n" );
1208  }
1209 
1210  }
1211 
1212 }
1213 #endif
1214 #ifdef DEBUG_PATCH_THRESHOLDS
1215 void L2gammaAlgo::printPatchConfig( int rdo )
1216 {
1217 
1218  // get the tower id
1219  int tow = mRdo2tower[rdo];
1220  int phi = tow/mNumEtas;
1221  int eta = tow%mNumEtas;
1222 
1223  int sec = phi/mNumSubs;
1224  int sub = phi%mNumSubs;
1225 
1226  /*
1227 
1228  The goal here is to output several quantities in an ascii
1229  format which matches the physical layout of the towers in
1230  order to help us debug problems online
1231 
1232  as least as far as the endcap is concerned (I can switch the
1233  order for the barrel if need be. Specifically, we want the
1234  following output:
1235 
1236  tower: 06TC05
1237  gain: 22.64
1238  ped: 20.90
1239  thr: 212 [2.8 GeV]
1240 
1241  pedestals gains, stat, fail, high tower frequency, etc...
1242 
1243  18.7 | -3.0 | 22.5
1244  -----+------+-----
1245  22.9 | 20.9 | 15.9
1246  -----+------+-----
1247  11.1 | 16.7 | 32.0
1248 
1249  */
1250 
1251  if ( mEEmc )
1252  fprintf(mLogFile,"etow: %02dT%c%02d\n",sec+1,sub+'A',eta+1);
1253  if ( mBEmc )
1254  fprintf(mLogFile,"btow: %02dt%c%02d\n",sec+1,sub+'a',eta+1);
1255 
1256  fprintf(mLogFile,"index: %4d\n",tow);
1257  fprintf(mLogFile,"gain: %6.3f [ideal=%6.3f]\n",mTowerGain[rdo],mTowerGainIdeal[rdo]);
1258  fprintf(mLogFile,"ped: %6.3f\n",mTowerPed[rdo]);
1259  fprintf(mLogFile,"thr: %3d [%4.2f GeV]\n",mTowerAdcThreshold[rdo],mTowerAdcThreshold[rdo]/mIdealGainT);
1260  fprintf(mLogFile,"patch: %6.3f\n",mPatchAdcThreshold[rdo]);
1261  fprintf(mLogFile,"\n");
1262 
1263  // 14.3 | 14.5 | 8.8 22.6 | 22.9 | 21.6 80 | 80 | 80 0 | 0 | 0 1 | 0 | 0
1264  fprintf(mLogFile,"loaded peds gain factors status bits mask bits frequency > thr\n");
1265 
1266  for ( int i_eta=eta+1;i_eta>=eta-1;i_eta-- ) {
1267 
1268  int phi_l = (phi+mNumPhis-1)%mNumPhis;
1269  int phi_r = (phi+mNumPhis+1)%mNumPhis;
1270  int phi_c = phi;
1271 
1272  int phis[]={phi_l,phi,phi_r};
1273  float peds[]={0,0,0};
1274  float gains[]={0,0,0};
1275  unsigned short stats[]={0xff,0xff,0xff};
1276  unsigned short fails[]={0xff,0xff,0xff};
1277  int freqs[]={0,0,0};
1278 
1279  float corr[]={0.,0.,0.};
1280 
1281  // get numbers from db
1282  if ( i_eta < mNumEtas && i_eta >= 0 ) {
1283 
1284  int tow_l = tower(phi_l,i_eta);
1285  int tow_r = tower(phi_r,i_eta);
1286  int tow_c = tower(phi_c,i_eta);
1287 
1288  int rdo_l = mTower2rdo[tow_l];
1289  int rdo_r = mTower2rdo[tow_r];
1290  int rdo_c = mTower2rdo[tow_c];
1291 
1292  peds[0]=mTowerPed[rdo_l];peds[1]=mTowerPed[rdo_c];peds[2]=mTowerPed[rdo_r];
1293  //gains[0]=mTowerAdcCorrection[rdo_l];gains[1]=mTowerAdcCorrection[rdo_c];gains[2]=mTowerAdcCorrection[rdo_r];
1294  gains[0]=(float)mIdealGainT;
1295  gains[1]=(float)mIdealGainT;
1296  gains[2]=(float)mIdealGainT;
1297  //$$$ gains[0]=mTowerGainIdeal[rdo_l]/m
1298  stats[0]=mTowerStat[rdo_l];stats[1]=mTowerStat[rdo_c];stats[2]=mTowerStat[rdo_r];
1299  fails[0]=mTowerFail[rdo_l];fails[1]=mTowerFail[rdo_c];fails[2]=mTowerFail[rdo_r];
1300  freqs[0]=mTowerFrequency[tow_l];freqs[1]=mTowerFrequency[tow_c];freqs[2]=mTowerFrequency[tow_r];
1301 
1302  corr[0]=mTowerAdcCorrection[rdo_l];
1303  corr[1]=mTowerAdcCorrection[rdo_c];
1304  corr[2]=mTowerAdcCorrection[rdo_r];
1305 
1306 
1307  }
1308 
1309  fprintf(mLogFile, "%4.1f | %4.1f | %4.1f\t", peds[0],peds[1],peds[2] );
1310  fprintf(mLogFile, "%4.0f | %4.0f | %4.0f\t", gains[0],gains[1],gains[2] );
1311  fprintf(mLogFile, "%2x | %2x | %2x \t", stats[0],stats[1],stats[2] );
1312  fprintf(mLogFile, "%2x | %2x | %2x \t", fails[0],fails[1],fails[2] );
1313  fprintf(mLogFile, "%4.2f | %4.2f | %4.2f\t",
1314  corr[0]*peds[0]/gains[0],
1315  corr[1]*peds[1]/gains[1],
1316  corr[2]*peds[2]/gains[2]);
1317 
1318  fprintf(mLogFile,"\n");
1319 
1320  if ( i_eta > eta-1 ) {
1321  fprintf(mLogFile, "-----+------+-----\t");
1322  fprintf(mLogFile, "-----+------+-----\t");
1323  fprintf(mLogFile, "-----+------+-----\t");
1324  fprintf(mLogFile, "-----+------+-----\t");
1325  fprintf(mLogFile, "-----+------+-----\t");
1326  fprintf(mLogFile, "\n" );
1327  }
1328  else
1329  fprintf(mLogFile, "\n\n" );
1330 
1331  }
1332 
1333 }
1334 
1335 #endif
L2gammaAlgo::mEEmc
unsigned short mEEmc
Definition: L2gammaAlgo.h:70
L2gammaAlgo::doEvent
bool doEvent(int inpEveId, TrgDataType *trgData, int emcIn, unsigned short *emcData)
Definition: L2gammaAlgo.cxx:627
L2gammaAlgo::L2gammaAlgo
L2gammaAlgo(const char *name, L2EmcDb *db, char *outDir, int resOff)
class constructor
Definition: L2gammaAlgo.cxx:17
L2gammaAlgo::mPatchAdcThreshold
float mPatchAdcThreshold[MAX_TOWERS]
Definition: L2gammaAlgo.h:237
L2gammaAlgo::mNumberInput
int mNumberInput
event counters
Definition: L2gammaAlgo.h:221
L2gammaAlgo::finishRun
void finishRun()
cleanup for runs. output of warnings (eg hot channels) etc...
Definition: L2gammaAlgo.cxx:917
L2gammaAlgo::printPatchConfig
void printPatchConfig(int rdo)
Definition: L2gammaAlgo.cxx:1110
L2gammaAlgo::mPatchThreshold
float mPatchThreshold
Patch threshold.
Definition: L2gammaAlgo.h:211
L2gammaAlgo::mTowerAdcThreshold
unsigned short mTowerAdcThreshold[MAX_TOWERS]
Definition: L2gammaAlgo.h:231
L2gammaAlgo::mRdo2tower
unsigned short mRdo2tower[MAX_TOWERS]
Map to go from Rdo to Tower and Tower to Rdo.
Definition: L2gammaAlgo.h:260
L2gammaAlgo::setLogFile
void setLogFile(const char *fname="./bsqueal.log")
Set logfile for summary output (overrides constructor option)
Definition: L2gammaAlgo.cxx:11
L2gammaAlgo::setPatchThreshold
void setPatchThreshold(float pt)
set the L2 trigger patch threshold. F_par[1]
Definition: L2gammaAlgo.h:306
L2gammaAlgo::mHistos
L2Histo mHistos[19]
L2 histograms.
Definition: L2gammaAlgo.h:287
L2gammaAlgo::mPatchStat
unsigned short mPatchStat[MAX_TOWERS]
Definition: L2gammaAlgo.h:253
L2gammaAlgo::mNumberLive
int mNumberLive
Counters for error monitoring.
Definition: L2gammaAlgo.h:283
L2gammaAlgo::mL2input
L2gammaAlgo * mL2input
Definition: L2gammaAlgo.h:203
L2gammaAlgo::initRun
int initRun(int run)
initialize the data structures for the next run
Definition: L2gammaAlgo.cxx:172
L2gammaAlgo::init
void init(int run, int I_par[5], float F_par[5])
One time initialization.
Definition: L2gammaAlgo.cxx:66
L2gammaAlgo::clear
void clear()
clear data structures for ...
Definition: L2gammaAlgo.cxx:907
L2gammaAlgo::setHistFile
void setHistFile(const char *fname="./bsqueal.dat")
Set filename for histogram output (overrides constructor option)
Definition: L2gammaAlgo.cxx:12
L2gammaAlgo::mNumPatch
unsigned short mNumPatch[MAX_TOWERS]
Number and list of up to 9 neighboring towers.
Definition: L2gammaAlgo.h:264
L2gammaAlgo::finish
void finish()
final output of statistics, histograms, etc...
Definition: L2gammaAlgo.cxx:1041
L2gammaAlgo::phibin
unsigned short phibin(unsigned short sec, unsigned short sub)
Log file.
Definition: L2gammaAlgo.h:295
L2gammaAlgo::mTowerThreshold
float mTowerThreshold
High tower threshold.
Definition: L2gammaAlgo.h:209
L2gammaAlgo::setTowerThreshold
void setTowerThreshold(float pt)
set the L2 high tower threshold. F_par[0]
Definition: L2gammaAlgo.h:302
L2gammaAlgo::mUseBbc
int mUseBbc
Option to correct for vertex effects (not in use yet)
Definition: L2gammaAlgo.h:215
L2gammaAlgo::mRunNumber
int mRunNumber
Run number for the current run.
Definition: L2gammaAlgo.h:206
L2gammaAlgo::mUseOfflineGains
int mUseOfflineGains
Option to use offline gains.
Definition: L2gammaAlgo.h:213
Generated by   doxygen 1.8.5