StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
L2jetAlgo2009.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: L2jetAlgo2009.cxx,v 1.1 2010/04/17 16:42:14 pibero Exp $
9  * \author Jan Balewski, IUCF, 2006
10  ***********************************************************
11  * Descripion:
12  Reco of mono- & di-jets in L2 using BTOW+ETOW
13  depends on L2-DB class
14  ***********************************************************
15  * version updated 2/17/2006 W.W. Jacobs
16  * 2008 to 2009 change by B.S. Page 11/12/2008
17  */
18 
19 
20 
21 #ifdef IS_REAL_L2 //in l2-ana environment
22  #include "../L2algoUtil/L2EmcDb.h"
23  #include "../L2algoUtil/L2Histo.h"
24 #else
25  #include "StTriggerUtilities/L2Emulator/L2algoUtil/L2EmcDb.h"
26  #include "StTriggerUtilities/L2Emulator/L2algoUtil/L2Histo.h"
27 #endif
28 
29 #include "L2jetAlgo2009.h"
30 #include "L2jetResults2009.h"
31 #include "Map_DeltaPhiJets.h"
32 
33 //=================================================
34 //=================================================
35 L2jetAlgo2009::L2jetAlgo2009(const char* name, L2EmcDb* db, char* outDir, int resOff, bool writeHighResult)
36  : L2VirtualAlgo2009( name, db, outDir, true, true, resOff) {
37  /* called one per days
38  all memory allocation must be done here
39  */
40  namechar='L';
41  if (writeHighResult) namechar='H';
42 
43  createHisto();
44  run_number=-1;
45  printf("L2jetAlgo2009 instantiated, logPath='%s'\n",mOutDir1.c_str());
46  eve_Jet[0]=new L2Jet;
47  eve_Jet[1]=new L2Jet;
48 }
49 
50 /*========================================
51  ======================================== */
52 bool
53 L2jetAlgo2009::paramsChanged( int *rc_ints, float *rc_floats) {
54  int i;
55  for(i=0;i<5;i++)
56  if(rc_ints[i]!=raw_ints[i]) goto foundProblem;
57 
58  for(i=0;i<5;i++)
59  if(fabs(rc_floats[i]-raw_floats[i])>0.00001) goto foundProblem;
60 
61  return false;
62 
63  foundProblem:
64  if (mLogFile) fprintf(mLogFile,"L2jet-ghost initRun - inconsistent params, ABORT initialization\n");
65  return true;
66 }
67 
68 /*========================================
69  ======================================== */
70 int
71 L2jetAlgo2009::initRunUser( int runNo, int *rc_ints, float *rc_floats) {
72 
73  // clear content of all histograms
74  int i;
75  for (i=0; i<mxHA;i++) if(hA[i])hA[i]->reset();
76 
77  /* .... clear content, set threshold @ max as default */
78  memset(db_btowL2PhiBin, 0, sizeof(db_btowL2PhiBin));
79  memset(db_btowL2PatchBin,0, sizeof(db_btowL2PatchBin));
80 
81  memset(db_etowL2PhiBin, 0 ,sizeof(db_etowL2PhiBin));
82  memset(db_etowL2PatchBin,0 ,sizeof(db_etowL2PatchBin));
83 
84  run_startUnix=time(0);
85  run_number =runNo; // serves as a flag this run is initialized
86  raw_ints =rc_ints;
87  raw_floats =rc_floats;
88  run_nEventOneJet=run_nEventDiJet= run_nEventRnd=0;
89 
90  mEventsInRun=0;
91 
92  if( mLogFile==0) printf(" L2jetAlgo2009() UNABLE to open run summary log file, continue anyhow\n");
93 
94  // unpack params from run control GUI
95  par_useBtowEast= (rc_ints[0]&1 )>0;
96  par_useBtowWest= (rc_ints[0]&2)>0;
97  par_useEndcap = rc_ints[1]&1;
98  int par_adcThr = rc_ints[2]; // needed only in initRun()
99  par_dbg = rc_ints[3];
100  par_RndAcceptPrescale = rc_ints[4];
101  par_minPhiBinDiff= 7; //for 1X1 JP was 5 for .6x.6 JP BP
102 
103  par_oneJetThr = rc_floats[0];
104  par_diJetThrHigh= rc_floats[1];
105  par_diJetThrLow = rc_floats[2];
106  par_diJetThr_2 = rc_floats[3]; // wwj 2/11
107  par_diJetThr_3 = rc_floats[4]; // wwj 2/11
108  par_diJetThr_4 = rc_floats[5]; // wwj 2/11
109  par_diJetThr_5 = rc_floats[6]; // wwj 2/11
110  par_diJetEtaHigh = rc_floats[7]; // wwj 2/11
111  par_diJetEtaLow = rc_floats[8]; // wwj 2/11
112 
113  // set other hardcoded or calculated params
114 
115  // to monitor hot towers in E+B Emc
116  float monTwEtThr=2.0; // (GeV) Et threshold,WARN, edit histo title by hand
117  par_hotTwEtThr = monTwEtThr;
118 
119  if(par_dbg) printf("Brian version is running!\n");
120 
121  if (mLogFile) {
122  fprintf(mLogFile,"L2jet algorithm initRun(%d), compiled: %s , %s\n params:\n",run_number,__DATE__,__TIME__);
123  fprintf(mLogFile," - use BTOW: East=%d West=%d, Endcap=%d L2ResOffset=%d\n", par_useBtowEast, par_useBtowWest,par_useEndcap ,mResultOffset);
124  fprintf(mLogFile," - threshold: ADC-ped> %d \n", par_adcThr);
125  fprintf(mLogFile," - min phi opening angle Jet1<->Jet2: %d in L2phiBins\n",par_minPhiBinDiff);
126  fprintf(mLogFile," - diJet Et thrHigh= %.2f (GeV) thrLow= %.2f (GeV)\n", par_diJetThrHigh, par_diJetThrLow);
127  fprintf(mLogFile," - new diJet Et thr_2= %.2f (GeV) ; thr_3= %.2f (Gev) ; thr_4= %.2f (GeV) ; thr_5= %.2f (GeV)\n",par_diJetThr_2, par_diJetThr_3, par_diJetThr_4, par_diJetThr_5); // wwj 2/10
128  fprintf(mLogFile," - oneJet Et thr = %.2f (GeV) ; rndAccPrescale=%d \n",par_oneJetThr,par_RndAcceptPrescale);
129  fprintf(mLogFile," - new diJet Eta_low= %.2f ; diJet Eta_High= %.2f \n", par_diJetEtaLow, par_diJetEtaHigh); // wwj 2/10
130  fprintf(mLogFile," - debug=%d, hot tower threshold: Et> %.1f GeV ( only monitoring)\n",par_dbg, monTwEtThr);
131  }
132 
133  // verify consistency of input params
134  int kBad=0;
135  kBad+=0x0001 * ( !par_useBtowEast & !par_useBtowWest & !par_useEndcap);
136  // par_pedOff no longer exists: kBad+=0x0002 * (par_adcThr<par_pedOff);
137  kBad+=0x0004 * (par_adcThr>16);
138  kBad+=0x0008 * (par_minPhiBinDiff<5);
139  kBad+=0x0010 * (par_minPhiBinDiff>=15);
140  kBad+=0x0020 * (par_oneJetThr<3.);
141  kBad+=0x0040 * (par_oneJetThr>12.);
142  kBad+=0x0080 * (par_diJetThrLow<1.9); // LOOK: wwj 2/11 using lowest threshold
143  kBad+=0x0100 * (par_diJetThrHigh<par_diJetThrLow);
144  kBad+=0x0200 * (par_diJetThrHigh>12.);
145  //0x0400 not used
146  kBad+=0x0800 * ( par_RndAcceptPrescale<0 );
147  kBad+=0x1000 * (par_diJetThr_2<par_diJetThrLow || par_diJetThr_5>par_diJetThrHigh); // wwj 2/16
148  kBad+=0x2000 * (par_diJetThr_2>par_diJetThr_3 || par_diJetThr_3>par_diJetThr_4 || par_diJetThr_4>par_diJetThr_5); // wwj 2/16
149  kBad+=0x4000 * (par_diJetEtaLow>par_diJetEtaHigh); // wwj 2/16
150  if (mLogFile) {
151  fprintf(mLogFile,"L2jet initRun() params checked for consistency, Error flag=0x%04x\n",kBad);
152  if(kBad) fprintf(mLogFile,"L2jet initRun() ABORT\n");
153  }
154 
155  if(kBad) {
156  run_number=-66;
157  if (mLogFile) {
158  fprintf(mLogFile,"L2jet algorithm init: crashB in internal logic\n");
159  fclose(mLogFile);
160  return kBad;
161  }
162  }
163 
164  char tit[100];
165  sprintf(tit,"# BTOW towers>ped+%d (input); x: # of towers/event",par_adcThr);
166  hA[47]->setTitle(tit);
167 
168  sprintf(tit,"# ETOW towers>ped+%d (input); x: # of towers/event",par_adcThr);
169  hA[48]->setTitle(tit);
170 
171 
172  const int mxEtaBinsE=12,mxEtaBinsB=40;
173 
174  // rebuild local lookup tables
175 
176  int etowEtaBin2Patch[mxEtaBinsE]={14,14,13,13,12,12,11,11,11,10,10,10};
177 
178  int nB=0, nE=0; /* counts # of unmasekd towers */
179  int nBg=0, nEg=0; /* counts # of reasonable calibrated towers */
180 
181  for(i=0; i<EmcDbIndexMax; i++) {
182  const L2EmcDb::EmcCDbItem *x=mDb->getByIndex(i);
183  if(mDb->isEmpty(x)) continue; /* dropped not mapped channels */
184  if(x->fail) continue; /* dropped masked channels */
185  if(x->gain<=0) continue; /* dropped uncalibrated towers , tmp */
186  /* if(x->sec!=1) continue; tmp, to test patch mapping */
187 
188  /* WARN, calculate index to match RDO order in the ADC data banks */
189  int ietaP, iphiP;
190  if (mDb->isBTOW(x) ) {
191  /*....... B A R R E L .................*/
192  nB++;
193  if(x->eta<0 || x->eta>mxEtaBinsB) goto crashIt_1;
194  if(!par_useBtowEast && x->eta<=20) continue;
195  if(!par_useBtowWest && x->eta>=21) continue;
196  ietaP= (x->eta-1)/4; /* correct */
197  int iphiTw=(x->sec-1)*10 + x->sub-'a';
198  // allign in phi TP @ L2 w/ L0
199  iphiTw--;
200  if(iphiTw<0) iphiTw=119;
201  // now cr0x1e, mod1, subm2, is beginning of the first BTOW TP
202  iphiP= iphiTw/4 ; /* correct */
203  // if(ietaP==0 && iphiP==5)
204  //if( (iphiTw==21 || iphiTw==22) && x->eta<=4) printf("%s %s ietaP=%d iPhiP=%d cr0x%02x ch=%03d\n",x->name, x->tube,ietaP, iphiP,x->crate,x->chan);
205  db_btowL2PhiBin[x->rdo]=iphiP;
206  db_btowL2PatchBin[x->rdo]=ietaP+ iphiP*cl2jetMaxEtaBins;
207  nBg++;
208  } else if(mDb->isETOW(x) && par_useEndcap) {
209  /*....... E N D C A P ........................*/
210  nE++;
211  int iphiTw= (x->sec-1)*5 + x->sub-'A';
212  // allign in phi TP @ L2 w/ L0
213  iphiTw--;
214  if(iphiTw<0) iphiTw=59;
215  // now subsector 01TB is beginning of the first(==0) ETOW TP in phi
216  iphiP= iphiTw/2 ; /* correct */
217  if(x->eta<0 || x->eta>mxEtaBinsE) goto crashIt_1;
218  ietaP=etowEtaBin2Patch[x->eta-1];
219  //printf("%s %s ietaP=%d iPhiP=%d\n",x->name, x->tube,ietaP, iphiP);
220  db_etowL2PhiBin[x->rdo]=iphiP;
221  db_etowL2PatchBin[x->rdo]=ietaP+ iphiP*cl2jetMaxEtaBins;
222  nEg++;
223  }
224 
225  }
226 
227  if (mLogFile) {
228  fprintf(mLogFile,"L2jet algorithm: found working/calibrated: %d/%d=ETOW & %d/%d=BTOW, based on ASCII DB\n",nE,nEg,nB,nBg);
229  }
230 
231  return 0; //OK
232 
233  crashIt_1: /* fatal initialization error */
234  run_number=-55;
235  if (mLogFile) {
236  fprintf(mLogFile,"L2jet algorithm init: crashC in internal logic\n");
237  fclose(mLogFile);
238  }
239  return -6;
240 
241 }
242 
243 /*========================================
244  ======================================== */
245 void L2jetAlgo2009::computeUser(int token){
246  //put in to satisfy l2virtualalgo2009
247 }
248 
249 /*========================================
250  ======================================== */
251 bool
252 L2jetAlgo2009::decisionUser(int token, int *myL2Result){
253  /* STRICT TIME BUDGET START ....*/
254 
255  //Ross says: Eventually you should take these out - timing is done in the virtual algo class as well.
256  unsigned long mEveTimeStart,mEveTimeStop, mEveTimeDiff;
257 rdtscl_macro(mEveTimeStart);
258 
259  //Ross says: I will check how the following line is done in L2bemcGamma2009
260  //following if statement was eve_ID!=inpEveId BP
261  if(1) {//UUUU
262  //.................... this event has NOT been processed
263 
264  /*
265  Chris doesn't want us to write anything out
266  during event processing ...
267  */
268 
269  bool bemcIn = true;
270  bool eemcIn = true;
271 
272  //eve_ID=inpEveId; // every events is processed only once
273  mEventsInRun++;
274  clearEvent(); /* price=13 kTicks */
275  int runTimeSec=time(0)- run_startUnix;
276  hA[10]->fill(0);
277  hA[12]->fill(runTimeSec);
278 
279  if(par_dbg>1) printf("\n......... in L2Jet_doEvent(ID=%d)... bIn=%d eIn=%d\n",eve_ID,bemcIn,eemcIn);//bemcIn and eemcIn were passed to the old doEvent
280 
281 
282  //Ross asks: will the jet code run on just barrel or just endcap?
283  // if that's the case, will it save time to move those sections of the
284  // code to computeUser() rather than decisionUser()?
285  if (bemcIn || eemcIn){//VVVV this algo has nothing to do w/o any Ecal data
286 
287  //===== step 1: unpack raw data blocks ==============
288 
289  int nBtowTw=0, nEtowTw=0;
290  //Ross says: I've changed projectAdc to use the L2eventStream2009 class
291  /*......... BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB */
292 
293  //we can move this to compute.
294  if(bemcIn==1 && (par_useBtowEast||par_useBtowWest) ) {
295 
296  nBtowTw=projectAdc( mEveStream_btow[token].get_hits(), mEveStream_btow[token].get_hitSize(),
297  db_btowL2PhiBin, db_btowL2PatchBin,
298  hA[20] );
299  }
300 
301  /*........... EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE */
302  //and this as well.
303  if(eemcIn==1 && par_useEndcap ) {
304  nEtowTw=projectAdc( mEveStream_etow[token].get_hits(), mEveStream_etow[token].get_hitSize(),
305  db_etowL2PhiBin, db_etowL2PatchBin,
306  hA[30] );
307  }
308 
309  //============Do 2D scan==========
310  float totEneGeV=true2Dscan();
311  int itotEneGeV=(int)totEneGeV;
312  int iK;
313  for(iK=0; iK<mxJ; iK++)
314  {
315  weightedEtaPhi(iK);
316  L2Jet *K=eve_Jet[iK];
317  K->eneGeV=K->iene;
318  K->phiRad=0.21*(6.0-K->fphiBin);
319  while(K->phiRad<0) K->phiRad+=6.2832;
320  while(K->phiRad>6.2832) K->phiRad-=6.2832;
321  }
322 
323  if(eve_Jet[0]->eneGeV <eve_Jet[1]->eneGeV) {// swap jets fo E1>E2
324  L2Jet *Jx=eve_Jet[0];
325  eve_Jet[0]=eve_Jet[1];
326  eve_Jet[1]=Jx;
327  }
328 
329  if(par_dbg>2) printf("doEvent iphiBin1=%d iene1=%f , iphiBin2=%d iene2=%f rms1PhiBin=%f \n",eve_Jet[0]->iphiBin,eve_Jet[0]->iene,eve_Jet[1]->iphiBin,eve_Jet[1]->iene,eve_Jet[0]->rmsPhiBin);
330 
331  //====== step 4: make trigger decisions====
332 
333  // now put in eta vs. Jet Threshold "contour" cuts ... this vesion cuts on eta_1 + eta_2 .. wwj 2/16
334  // looks like fetaBin goes 0.0->15.0 for eta -1. -> +2. ... so "eta" = fetaBin*0.2 -1.0 ... wwj 2/10
335 
336  float rjetEta_0=eve_Jet[0]->fetaBin*0.2 -1.0; // wwj 2/10
337  float rjetEta_1=eve_Jet[1]->fetaBin*0.2 -1.0; // wwj 2/10
338 
339  // printf(" fetabin_0=%f fetabin_1=%f \n",eve_Jet[0]->fetaBin,eve_Jet[1]->fetaBin); // wwj 2/15
340  // printf(" rjetEta_0=%f rjetEta_1=%f \n",rjetEta_0, rjetEta_1); // wwj 2/15
341 
342  bool acceptDiJet_EE=((rjetEta_0+rjetEta_1) > par_diJetEtaHigh && eve_Jet[0]->eneGeV > par_diJetThr_4 && eve_Jet[1]->eneGeV > par_diJetThrLow); // wwj 2/10
343 
344  bool acceptDiJet_EB=((rjetEta_0+rjetEta_1) > par_diJetEtaLow && eve_Jet[0]->eneGeV > par_diJetThr_5 && eve_Jet[1]->eneGeV > par_diJetThr_2); // wwj 2/10
345 
346 
347  // bool acceptDiJet=( eve_Jet[0]->eneGeV > par_diJetThrHigh) && ( eve_Jet[1]->eneGeV > par_diJetThrLow);
348  bool acceptDiJet_BB=(eve_Jet[0]->eneGeV > par_diJetThrHigh && eve_Jet[1]->eneGeV > par_diJetThr_3); // wwj 2/10
349 
350  bool acceptDiJet=acceptDiJet_EE || acceptDiJet_EB || acceptDiJet_BB; // wwj 2/10
351 
352  bool acceptOneJet=( eve_Jet[0]->eneGeV > par_oneJetThr) ;
353 
354  bool acceptRnd=mRandomAccept>0; // provided by Virtual algo
355  mAccept=acceptDiJet || acceptOneJet || acceptRnd;
356 
357  //====== step 5: update various monitorig histos
358 
359  // histogramming reco Et1-Et2, no cuts
360  int iet1 =(int)eve_Jet[0]->eneGeV;
361  int iet2 =(int)eve_Jet[1]->eneGeV;
362  int ieta1=(int)eve_Jet[0]->fetaBin;
363  int ieta2=(int)eve_Jet[1]->fetaBin;
364  int iphi1=(int)eve_Jet[0]->fphiBin;
365  int iphi2=(int)eve_Jet[1]->fphiBin;
366 
367  hA[40]->fill(iet1,iet2);
368 
369  hA[41]->fill(ieta1,iphi1);
370  hA[42]->fill(ieta2,iphi2);
371  hA[43]->fill(iphi1,iphi2);
372  hA[44]->fill(iet1);
373  hA[45]->fill(iet2);
374  hA[46]->fill(itotEneGeV);
375  hA[47]->fill(nBtowTw);
376  hA[48]->fill(nEtowTw);
377 
378  // sivers delta zeta, the map is still worng
379  int kphi1=int(eve_Jet[0]->phiRad*10.);
380  int kphi2=int(eve_Jet[1]->phiRad*10.);
381  int idelZeta=map_DelPhiJets[kphi1*MxPhiRad10 + kphi2];
382 
383  if( mAccept) hA[10]->fill(8);
384 
385  if(acceptOneJet ){
386  hA[10]->fill(4);
387  run_nEventOneJet++;
388  hA[13]->fill(runTimeSec);
389  hA[50]->fill(iet1);
390  hA[51]->fill(ieta1,iphi1);
391  hA[52]->fill(ieta1);
392  hA[53]->fill(iphi1);
393  }
394 
395  if(acceptDiJet ){
396  hA[10]->fill(5);
397  run_nEventDiJet++;
398  hA[14]->fill(runTimeSec);
399  hA[60]->fill(iet1,iet2);
400  hA[61]->fill(ieta1,iphi1);
401  hA[62]->fill(ieta2,iphi2);
402  hA[63]->fill(iphi1,iphi2);
403  hA[64]->fill(iet1);
404  hA[65]->fill(iet2);
405  hA[66]->fill(ieta1);
406  hA[67]->fill(ieta2);
407  hA[68]->fill(iphi1);
408  hA[69]->fill(iphi2);
409  hA[70]->fill(idelZeta);
410  hA[71]->fill(ieta1,idelZeta);
411  hA[72]->fill(ieta1,ieta2);
412  hA[73]->fill((iphi1+iphi2)/2,idelZeta);
413  hA[74]->fill(itotEneGeV);
414  }
415 
416  if(acceptDiJet_EE ){
417  run_nEventDiJet++; // all these histos #80-#94 added ... wwj 2/10
418  hA[80]->fill(iet1,iet2);
419  hA[81]->fill(ieta1,iphi1);
420  hA[82]->fill(ieta2,iphi2);
421  hA[83]->fill(iphi1,iphi2);
422  hA[84]->fill(iet1);
423  hA[85]->fill(iet2);
424  hA[86]->fill(ieta1);
425  hA[87]->fill(ieta2);
426  hA[88]->fill(iphi1);
427  hA[89]->fill(iphi2);
428  hA[90]->fill(idelZeta);
429  hA[91]->fill(ieta1,idelZeta);
430  hA[92]->fill(ieta1,ieta2);
431  hA[93]->fill((iphi1+iphi2)/2,idelZeta);
432  hA[94]->fill(itotEneGeV);
433  }
434 
435  if(acceptDiJet_EB ){
436  run_nEventDiJet++; // all these histos #100-#114 added ... wwj 2/10
437  hA[100]->fill(iet1,iet2);
438  hA[101]->fill(ieta1,iphi1);
439  hA[102]->fill(ieta2,iphi2);
440  hA[103]->fill(iphi1,iphi2);
441  hA[104]->fill(iet1);
442  hA[105]->fill(iet2);
443  hA[106]->fill(ieta1);
444  hA[107]->fill(ieta2);
445  hA[108]->fill(iphi1);
446  hA[109]->fill(iphi2);
447  hA[110]->fill(idelZeta);
448  hA[111]->fill(ieta1,idelZeta);
449  hA[112]->fill(ieta1,ieta2);
450  hA[113]->fill((iphi1+iphi2)/2,idelZeta);
451  hA[114]->fill(itotEneGeV);
452  }
453 
454  if(acceptDiJet_BB ){
455  run_nEventDiJet++; // all these histos #120-#134 added ... wwj 2/10
456  hA[120]->fill(iet1,iet2);
457  hA[121]->fill(ieta1,iphi1);
458  hA[122]->fill(ieta2,iphi2);
459  hA[123]->fill(iphi1,iphi2);
460  hA[124]->fill(iet1);
461  hA[125]->fill(iet2);
462  hA[126]->fill(ieta1);
463  hA[127]->fill(ieta2);
464  hA[128]->fill(iphi1);
465  hA[129]->fill(iphi2);
466  hA[130]->fill(idelZeta);
467  hA[131]->fill(ieta1,idelZeta);
468  hA[132]->fill(ieta1,ieta2);
469  hA[133]->fill((iphi1+iphi2)/2,idelZeta);
470  hA[134]->fill(itotEneGeV);
471  }
472 
473  if(mRandomAccept){
474  hA[10]->fill(6);
475  run_nEventRnd++;
476  hA[15]->fill(runTimeSec);
477  }
478 
479 
480  //====== step 6: fill L2Result (except time)
481  L2jetResults2009 out; // all output bits lump together
482  memset(&out,0,sizeof(out)); // clear content
483 
484 
485  // add acceptDiJet_EE, acceptDiJet_EB and acceptDiJet_BB to below ... wwj 2/10
486  out.int0.version=L2JET_RESULTS_VERSION;
487  out.int0.decision=
488  ( par_useBtowEast <<0 ) +
489  ( par_useBtowWest <<1 ) +
490  ( par_useEndcap <<2 ) +
491  ( (bemcIn>0) <<3 ) +
492  ( (eemcIn>0) <<4 ) +
493  ( mRandomAccept <<5 ) +
494  ( acceptOneJet <<6 ) +
495  ( acceptDiJet <<7 ) +
496  ( acceptDiJet_EE <<8 ) +
497  ( acceptDiJet_EB <<9 ) +
498  ( acceptDiJet_BB <<10 ) ; // e.g., added last three lines here ... wwj 2/10
499  out.int0.dumm=namechar;
500 
501  out.int1.iTotEne=(unsigned short)(totEneGeV*100.); // now 1=10 MeV
502  out.int2.nBtowTw=nBtowTw;
503  out.int2.nEtowTw=nEtowTw;
504 
505  out.jet1.jPhi=(int)(eve_Jet[0]->phiRad*28.65); //so phi/deg=2*jPhi
506  out.jet1.jEta=(int)(eve_Jet[0]->fetaBin*10.);
507  out.jet1.iEne=(unsigned short)(eve_Jet[0]->eneGeV*100.); // now 1=10 MeV
508 
509  out.jet2.jPhi=(int)(eve_Jet[1]->phiRad*28.65);
510  out.jet2.jEta=(int)(eve_Jet[1]->fetaBin*10.);
511  out.jet2.iEne=(unsigned short)(eve_Jet[1]->eneGeV*100.); // now 1=10 MeV
512 
513  out.jet1.rmsEta=(unsigned short)(eve_Jet[0]->rmsEtaBin*20.); //Fix: 1=0.01 physical eta
514  out.jet1.rmsPhi=(unsigned short)(eve_Jet[0]->rmsPhiBin*120.); //1=0.1 deg
515 
516  out.jet2.rmsEta=(unsigned short)(eve_Jet[1]->rmsEtaBin*20.); //Fix: 1=0.01 physical eta
517  out.jet2.rmsPhi=(unsigned short)(eve_Jet[1]->rmsPhiBin*120.);// 1=0.1 deg
518 
519 
520  rdtscl_macro(mEveTimeStop);
521  mEveTimeDiff=mEveTimeStop-mEveTimeStart;
522  int kTick=mEveTimeDiff/1000;
523  // printf("jj=%f t1=%d t2=%d \n",mEveTimeDiff/1000.,mEveTimeStart,mEveTimeStop);
524  //mhT->fill(kTick); taken out 12/11/08 BP
525 
526  out.int0.kTick= kTick>255 ? 255 : kTick;
527 
528  //calculate and match the check sum
529  out.int1.checkSum=-L2jetResults2009_doCheckSum(&out);
530  // unsigned char cSum=L2jetResults2009_doCheckSum(&out); printf("cSum2=%d\n",cSum);
531 
532 
533  //===== step 7: write L2Result
534  memcpy(myL2Result,&out,sizeof(L2jetResults2009));
535 
536  // dirty tests, clean it up before real use
537 
538  if(par_dbg){//WWWW
539  L2jetResults2009_print(&out);
540  printf(" phiRad1=%f phiRad2=%f \n",eve_Jet[0]->phiRad,eve_Jet[1]->phiRad);
541  printf("idelZeta=%d delZeta/deg=%.1f \n\n",idelZeta,idelZeta/31.416*180);
542 
543  //tmp printouts of errors:
544  if( out.jet1.iEne+out.jet2.iEne > out.int1.iTotEne) {
545  //printf("L2jet-fatal error, eve=%d, iEtot=%d < iEJ1=%d + iEJ2=%d, continue\n",inpEveId, out.int1.iTotEne,out.jet1.iEne,out.jet2.iEne);//no inpEveId BP
546  }
547  if(iphi1==iphi2) {
548  printf("L2jet-fatal error,neveId=%d, phi1,2=%d,%d\n",mEventsInRun,iphi1,iphi2);
549  dumpPatchEneA();
550  }
551 
552  if( L2jetResults2009_doCheckSum(&out)) {
553  printf("L2jet-fatal error, wrong cSum=%d\n", L2jetResults2009_doCheckSum(&out));
554  L2jetResults2009_print(&out);
555  }
556 
557  } // end of WWWW
558 
559  }// end of VVVV (etow or btow has some data)
560  }// end of UUUU event processing
561 
562 
563  return mAccept;
564 }
565 
566 
567 /*========================================
568  ======================================== */
569 void
570 L2jetAlgo2009::finishRunUser() { /* called once at the end of the run */
571  if(run_number<0) return; // already finished
572 
573  if (mLogFile) {
574  fprintf(mLogFile,"L2-jet algorithm finishRun(%d)\n",run_number);
575  fprintf(mLogFile," - %d events seen by L2 di-jet\n",mEventsInRun);
576  fprintf(mLogFile," - accepted: rnd=%d oneJet=%d diJet=%d \n", run_nEventRnd, run_nEventOneJet, run_nEventDiJet);
577 
578  // print few basic histos
579  //printf("Attempting to write to 10:\n");
580  hA[10]->printCSV(mLogFile); // event accumulated
581 
582  }
583  finishRunHisto(); // still needs current DB
584  if( mHistFile==0) {
585  printf(" L2jetAlgo2009: finishRun() UNABLE to open run summary log file, continue anyhow\n");
586  if (mLogFile)
587  fprintf(mLogFile,"L2 di-jet histos NOT saved, I/O error\n");
588  } else { // save histos
589  int j;
590  int nh=0;
591  for(j=0;j<mxHA;j++) {
592  if(hA[j]==0) continue;
593  hA[j]->write(mHistFile);
594  nh++;
595  }
596  }
597 
598  run_number=-2; // clear run #
599 
600 }
601 
602 //=======================================
603 //=======================================
604 void
605 L2jetAlgo2009::createHisto() {
606  memset(hA,0,sizeof(hA));
607 
608  hA[10]=new L2Histo(10,"total event counter; x=cases",9);
609  hA[11]=new L2Histo(11,"L2 time used per input event; x: time (CPU kTics), range=100muSec; y: events ",160);
610 
611  int mxRunDration=2500;
612  hA[12]=new L2Histo(12,"rate of input events; x: time in this run (seconds); y: rate (Hz)", mxRunDration);
613 
614  hA[13]=new L2Histo(13,"rate of accepted one-Jet; x: time in this run (seconds); y: rate (Hz)", mxRunDration);
615  hA[14]=new L2Histo(14,"rate of accepted di-Jet ; x: time in this run (seconds); y: rate (Hz)", mxRunDration);
616  hA[15]=new L2Histo(15,"rate of random accepted ; x: time in this run (seconds); y: rate (Hz)", mxRunDration);
617 
618  // BTOW raw spectra
619  hA[20]=new L2Histo(20,"BTOW tower, Et>2.0 GeV (input); x: BTOW RDO index=chan*30+fiber; y: counts", 4800);
620  hA[21]=new L2Histo(21,"BTOW tower, Et>2.0 GeV (input); x: BTOW softID", 4800);
621  hA[22]=new L2Histo(22,"BTOW tower, Et>2.0 GeV (input); x: eta bin, [-1,+1]; y: phi bin ~sector",40,120);
622 
623  // ETOW raw spectra
624  hA[30]=new L2Histo(30,"ETOW tower, Et>2.0 GeV (input); x: ETOW RDO index=chan*6+fiber; y: counts", 720 );
625  hA[31]=new L2Histo(31,"ETOW tower, Et>2.0 GeV (input); x: i=chan+128*crate", 768);
626  hA[32]=new L2Histo(32,"ETOW tower, Et>2.0 GeV (input); x: 12 - Endcap etaBin ,[+1,+2]; y: phi bin ~sector",12,60);
627 
628  // Di-Jet raw yields
629  hA[40]=new L2Histo(40,"Et Jet1-Jet2 (input); x: Jet1 Et/GeV ; Jet2 Et/GeV",30,30);
630  hA[41]=new L2Histo(41,"diJet1 eta-phi (input); x: iEta [-1,+2] ; y: iPhi ~sector ",15,30);
631  hA[42]=new L2Histo(42,"diJet2 eta-phi (input); x: iEta [-1,+2] ; y: iPhi ~sector",15,30);
632 
633  hA[43]=new L2Histo(43,"diJet phi1-phi2 (input); x: iPhi1 ~sector ; y: iPhi2 ~sector ",30,30);
634 
635  hA[44]=new L2Histo(44,"Jet1 Et (input); x: Et (GeV)", 60);
636  hA[45]=new L2Histo(45,"Jet2 Et (input); x: Et (GeV)", 60);
637  hA[46]=new L2Histo(46,"total Et (input); x: Et (GeV)", 100);
638  hA[47]=new L2Histo(47,"# BTOW towers>thrXX (input); x: # of towers/event", 200);
639  hA[48]=new L2Histo(48,"# ETOW towers>thrXX (input); x: # of towers/event", 100);
640 
641  // ........accepted one-jet events
642  hA[50]=new L2Histo(50,"one-Jet Et (accepted); x: jet Et (GeV)", 60);
643  hA[51]=new L2Histo(51,"one-Jet eta-phi (accepted); x: iEta [-1,+2] ; y: iPhi ~sector ",15,30);
644  hA[52]=new L2Histo(52,"one-Jet eta (accepted); x: iEta [-1,+2]", 15);
645  hA[53]=new L2Histo(53,"one-Jet phi (accepted); x: iPhi ~sector", 30);
646 
647  // Di-Jet accepted
648  hA[60]=new L2Histo(60,"Et of Jet1 vs. Jet2 (accepted); x: Jet1/GeV ; Jet2/GeV",30,30);
649  hA[61]=new L2Histo(61,"diJet1 eta-phi (accepted); x: iEta [-1,+2] ; y: iPhi ~sector ",15,30);
650  hA[62]=new L2Histo(62,"diJet2 eta-phi (accepted); x: iEta [-1,+2] ; y: iPhi ~sector",15,30);
651 
652  hA[63]=new L2Histo(63,"diJet phi1-phi2 (accepted); x: iPhi1 ~sector ; y: iPhi2 ~sector ",30,30);
653 
654  hA[64]=new L2Histo(64,"diJet1 Et (accepted); x: Et (GeV)", 60);
655  hA[65]=new L2Histo(65,"diJet2 Et (accepted); x: Et (GeV)", 60);
656 
657  hA[66]=new L2Histo(66,"diJet1 eta (accepted); x: i Eta [-1,+2]", 15);
658  hA[67]=new L2Histo(67,"diJet2 eta (accepted); x: i Eta [-1,+2]", 15);
659  hA[68]=new L2Histo(68,"diJet1 phi (accepted); x: iPhi ~sector", 30);
660  hA[69]=new L2Histo(69,"diJet2 phi (accepted); x: iPhi ~sector", 30);
661  hA[70]=new L2Histo(70,"diJet delZeta (accepted); x: delta zeta (rad*10)", MxPhiRad10);
662  hA[71]=new L2Histo(71,"diJet delZeta vs. eta1 (accepted); x: iEta1 [-1,+2] ; y: delta zeta (rad*10)",15, MxPhiRad10);
663  hA[72]=new L2Histo(72,"diJet eta2 vs. eta1 (accepted); x: iEta1 [-1,+2] ;x: iEta2 [-1,+2] ",15,15);
664  hA[73]=new L2Histo(73,"diJet delZeta vs. avrPhi (accepted); x: (iphi1+iphi2)/2 (12 deg/bin); y: delta zeta (rad*10)",30, MxPhiRad10);
665  hA[74]=new L2Histo(74,"total Et diJet (accepted); x: Et (GeV)", 100);
666 
667  // Di-Jet_EE accepted ... all these histos #80-#94 added ... wwj 2/10
668  hA[80]=new L2Histo(80,"Et of Jet1 vs. Jet2 (accepted_EE); x: Jet1/GeV ; Jet2/GeV",30,30);
669  hA[81]=new L2Histo(81,"diJet1 eta-phi (accepted_EE); x: iEta [-1,+2] ; y: iPhi ~sector ",15,30);
670  hA[82]=new L2Histo(82,"diJet2 eta-phi (accepted_EE); x: iEta [-1,+2] ; y: iPhi ~sector",15,30);
671 
672  hA[83]=new L2Histo(83,"diJet phi1-phi2 (accepted_EE); x: iPhi1 ~sector ; y: iPhi2 ~sector ",30,30);
673 
674  hA[84]=new L2Histo(84,"diJet1 Et (accepted_EE); x: Et (GeV)", 60);
675  hA[85]=new L2Histo(85,"diJet2 Et (accepted_EE); x: Et (GeV)", 60);
676 
677  hA[86]=new L2Histo(86,"diJet1 eta (accepted_EE); x: i Eta [-1,+2]", 15);
678  hA[87]=new L2Histo(87,"diJet2 eta (accepted_EE); x: i Eta [-1,+2]", 15);
679  hA[88]=new L2Histo(88,"diJet1 phi (accepted_EE); x: iPhi ~sector", 30);
680  hA[89]=new L2Histo(89,"diJet2 phi (accepted_EE); x: iPhi ~sector", 30);
681  hA[90]=new L2Histo(90,"diJet delZeta (accepted_EE); x: delta zeta (rad*10)", MxPhiRad10);
682  hA[91]=new L2Histo(91,"diJet delZeta vs. eta1 (accepted_EE); x: iEta1 [-1,+2] ; y: delta zeta (rad*10)",15, MxPhiRad10);
683  hA[92]=new L2Histo(92,"diJet eta2 vs. eta1 (accepted_EE); x: iEta1 [-1,+2] ;x: iEta2 [-1,+2] ",15,15);
684  hA[93]=new L2Histo(93,"diJet delZeta vs. avrPhi (accepted_EE); x: (iphi1+iphi2)/2 (12 deg/bin); y: delta zeta (rad*10)",30, MxPhiRad10);
685  hA[94]=new L2Histo(94,"total Et diJet (accepted_EE); x: Et (GeV)", 60);
686 
687 // Di-Jet_EB accepted ... all these histos #100-#114 added ... wwj 2/10
688  hA[100]=new L2Histo(100,"Et of Jet1 vs. Jet2 (accepted_EB); x: Jet1/GeV ; Jet2/GeV",30,30);
689  hA[101]=new L2Histo(101,"diJet1 eta-phi (accepted_EB); x: iEta [-1,+2] ; y: iPhi ~sector ",15,30);
690  hA[102]=new L2Histo(102,"diJet2 eta-phi (accepted_EB); x: iEta [-1,+2] ; y: iPhi ~sector",15,30);
691 
692  hA[103]=new L2Histo(103,"diJet phi1-phi2 (accepted_EB); x: iPhi1 ~sector ; y: iPhi2 ~sector ",30,30);
693 
694  hA[104]=new L2Histo(104,"diJet1 Et (accepted_EB); x: Et (GeV)", 60);
695  hA[105]=new L2Histo(105,"diJet2 Et (accepted_EB); x: Et (GeV)", 60);
696 
697  hA[106]=new L2Histo(106,"diJet1 eta (accepted_EB); x: i Eta [-1,+2]", 15);
698  hA[107]=new L2Histo(107,"diJet2 eta (accepted_EB); x: i Eta [-1,+2]", 15);
699  hA[108]=new L2Histo(108,"diJet1 phi (accepted_EB); x: iPhi ~sector", 30);
700  hA[109]=new L2Histo(109,"diJet2 phi (accepted_EB); x: iPhi ~sector", 30);
701  hA[110]=new L2Histo(110,"diJet delZeta (accepted_EB); x: delta zeta (rad*10)", MxPhiRad10);
702  hA[111]=new L2Histo(111,"diJet delZeta vs. eta1 (accepted_EB); x: iEta1 [-1,+2] ; y: delta zeta (rad*10)",15, MxPhiRad10);
703  hA[112]=new L2Histo(112,"diJet eta2 vs. eta1 (accepted_EB); x: iEta1 [-1,+2] ;x: iEta2 [-1,+2] ",15,15);
704  hA[113]=new L2Histo(113,"diJet delZeta vs. avrPhi (accepted_EB); x: (iphi1+iphi2)/2 (12 deg/bin); y: delta zeta (rad*10)",30, MxPhiRad10);
705  hA[114]=new L2Histo(114,"total Et diJet (accepted_EB); x: Et (GeV)", 60);
706 
707 // Di-Jet_BB accepted ... all these histos #120-#134 added ... wwj 2/10
708  hA[120]=new L2Histo(120,"Et of Jet1 vs. Jet2 (accepted_BB); x: Jet1/GeV ; Jet2/GeV",30,30);
709  hA[121]=new L2Histo(121,"diJet1 eta-phi (accepted_BB); x: iEta [-1,+2] ; y: iPhi ~sector ",15,30);
710  hA[122]=new L2Histo(122,"diJet2 eta-phi (accepted_BB); x: iEta [-1,+2] ; y: iPhi ~sector",15,30);
711 
712  hA[123]=new L2Histo(123,"diJet phi1-phi2 (accepted_BB); x: iPhi1 ~sector ; y: iPhi2 ~sector ",30,30);
713 
714  hA[124]=new L2Histo(124,"diJet1 Et (accepted_BB); x: Et (GeV)", 60);
715  hA[125]=new L2Histo(125,"diJet2 Et (accepted_BB); x: Et (GeV)", 60);
716 
717  hA[126]=new L2Histo(126,"diJet1 eta (accepted_BB); x: i Eta [-1,+2]", 15);
718  hA[127]=new L2Histo(127,"diJet2 eta (accepted_BB); x: i Eta [-1,+2]", 15);
719  hA[128]=new L2Histo(128,"diJet1 phi (accepted_BB); x: iPhi ~sector", 30);
720  hA[129]=new L2Histo(129,"diJet2 phi (accepted_BB); x: iPhi ~sector", 30);
721  hA[130]=new L2Histo(130,"diJet delZeta (accepted_BB); x: delta zeta (rad*10)", MxPhiRad10);
722  hA[131]=new L2Histo(131,"diJet delZeta vs. eta1 (accepted_BB); x: iEta1 [-1,+2] ; y: delta zeta (rad*10)",15, MxPhiRad10);
723  hA[132]=new L2Histo(132,"diJet eta2 vs. eta1 (accepted_BB); x: iEta1 [-1,+2] ;x: iEta2 [-1,+2] ",15,15);
724  hA[133]=new L2Histo(133,"diJet delZeta vs. avrPhi (accepted_BB); x: (iphi1+iphi2)/2 (12 deg/bin); y: delta zeta (rad*10)",30, MxPhiRad10);
725  hA[134]=new L2Histo(134,"total Et diJet (accepted_BB); x: Et (GeV)", 100);
726 
727 }
728 
729 //=======================================
730 //=======================================
731 void
732 L2jetAlgo2009::clearEvent(){
733  /* printf("clearEvent_L2jet() executed\n"); */
734 
735  mAccept=false;
736  memset(eve_patchEne,0,sizeof(eve_patchEne));
737  memset(eve_phiEne,0,sizeof(eve_phiEne));
738  eve_Jet[0]->clear();
739  eve_Jet[1]->clear();
740 }
741 
742 
743 
744 //=======================================
745 //=======================================
746 int
747 L2jetAlgo2009::projectAdc(const HitTower1 *hit,const int hitSize,
748  unsigned short *phiBin, unsigned short *patchBin,
749  L2Histo *hHot ){
750  int tmpNused=0; /* counts mapped & used ADC channels */
751  short rdo;
752  float low_noise_et;
753 
754  for(int i=0;i< hitSize;i++,hit++) {
755  rdo=hit->rdo;
756  low_noise_et=hit->low_noise_et;
757  eve_patchEne[patchBin[rdo]]+=low_noise_et;
758  eve_phiEne[phiBin[rdo]]+=low_noise_et;
759  tmpNused++;
760 
761  if(low_noise_et > par_hotTwEtThr) hHot->fill(rdo);//do we want rdo or i?
762 
763  }
764  return tmpNused;
765 }
766 
767 
768 //========================================
769 //========================================
770 float L2jetAlgo2009::true2Dscan(){
771  //implament a true 2-D scan to test the performance of the current algo
772  int iphiBinEdge = 0;
773  int ietaBinEdge = 0;
774  float maxPatchEt = 0.0;
775  float sumTot = 0.0;
776 
777  float *totalCaloEt=eve_phiEne;
778  float eneA[cl2jetMaxEtaBins];
779  float secondPatchArray[cl2jetMaxPhiBins][cl2jetMaxEtaBins-cl2jet_par_mxEtaBin+1];
780  memset(secondPatchArray,0,sizeof(secondPatchArray));
781 
782  int i;//loops over all phi bin. phi edge
783  int j;//collapses phi bins
784  int j1;//phi wrapped phi bin
785  int j2;//loops over eta bins and sums energy
786  int k;//loops over eta bins. eta edge
787 
788  for(i=0;i<cl2jetMaxPhiBins;i++)
789  {
790  sumTot+=totalCaloEt[i];
791  memset(eneA,0,sizeof(eneA));
792  for(j=0;j<cl2jet_par_mxPhiBin;j++)
793  {
794  j1 = (i+j)%cl2jetMaxPhiBins;//must wrap up phi
795  float *patchEneA=eve_patchEne+(j1*cl2jetMaxEtaBins);
796  for(j2=0;j2<cl2jetMaxEtaBins;j2++,patchEneA++)
797  {
798  eneA[j2]+=*patchEneA;
799  }
800  }
801  float *eneAp=eneA;
802  for(k=0;k<cl2jetMaxEtaBins-cl2jet_par_mxEtaBin+1;k++,eneAp++)
803  {
804  // printf(" J1 ieta=%d E=%.2f\n", ix,eneA[0]/l2jet_par_energyScale);
805  // printf(" E=%.2f\n",etaEneA[ix]/l2jet_par_energyScale);
806  float sum=eneAp[0]+eneAp[1]+eneAp[2]+eneAp[3]+eneAp[4];
807  // printf("ix=%d sum=%d sumMax=%d\n",ix,sum,sumMax);
808  secondPatchArray[i][k]=sum;
809  if(maxPatchEt>sum) continue;
810  //float *eneSave=eneA;
811  maxPatchEt=sum;
812  iphiBinEdge=i;
813  ietaBinEdge=k;
814  }
815  }
816 
817  eve_Jet[0]->iphiBin=iphiBinEdge;
818  eve_Jet[0]->ietaBin=ietaBinEdge;
819  eve_Jet[0]->iene=maxPatchEt;
820 
821  //now find second highest patch
822  int iphiBinEdge2=-1;
823  int ietaBinEdge2=-1;
824  float maxPatchEt2=0.0;
825  char doWrap=0;
826  int a1=iphiBinEdge-par_minPhiBinDiff;
827  int a2=iphiBinEdge+par_minPhiBinDiff;
828  if (a1<0) { a1+=cl2jetMaxPhiBins; doWrap+=1; }
829  if (a2>=cl2jetMaxPhiBins) { a2-=cl2jetMaxPhiBins; doWrap+=2; }
830 
831  int b;
832  int b1;
833  if (!doWrap)//masking range does not wrap up
834  {
835  for(b=0;b<cl2jetMaxPhiBins;b++)
836  {
837  if(b>=a1 && b<=a2) continue;
838  for(b1=0;b1<cl2jetMaxEtaBins-cl2jet_par_mxEtaBin+1;b1++)
839  {
840  if(maxPatchEt2>secondPatchArray[b][b1]) continue;
841  maxPatchEt2=secondPatchArray[b][b1];
842  iphiBinEdge2=b;
843  ietaBinEdge2=b1;
844  }
845  }
846  }
847  else
848  {//masking range wraps up
849  for(b=a2;b<a1;b++)
850  {
851  for(b1=0;b1<cl2jetMaxEtaBins-cl2jet_par_mxEtaBin+1;b1++)
852  {
853  if(maxPatchEt2>secondPatchArray[b][b1]) continue;
854  maxPatchEt2=secondPatchArray[b][b1];
855  iphiBinEdge2=b;
856  ietaBinEdge2=b1;
857  }
858  }
859  }
860 
861  eve_Jet[1]->iphiBin=iphiBinEdge2;
862  eve_Jet[1]->ietaBin=ietaBinEdge2;
863  eve_Jet[1]->iene=maxPatchEt2;
864 
865  return sumTot;
866 }
867 
868 //========================================
869 //========================================
870 void L2jetAlgo2009::weightedEtaPhi(int iK){
871  L2Jet *J=eve_Jet[iK];
872 
873  // empty calo protection
874  if(J->iene<=0.1)
875  {
876  J->fphiBin=J->iphiBin+.333;
877  J->fetaBin=J->ietaBin+.333;
878  J->rmsPhiBin=0.0;
879  J->rmsEtaBin=0.0;
880  return;
881  }
882 
883  int iphi0=J->iphiBin;
884  int ieta0=J->ietaBin;
885  float iene0=J->iene;
886 
887  //variables for weighted phi
888  float sumY=0.0;
889  float sum1;//wrk variable
890 
891  //these variables for phi RMS:
892  double sumYY = 0.0;
893  int nEnePhi=0; // counts bins with energy for RMS computation
894 
895  //variable for weighted eta
896  float sumX=0.0;
897 
898  //these variables for eta RMS:
899  double sumXX = 0.0;
900  int nEneEta = 0;
901 
902  float etaEneA[cl2jet_par_mxEtaBin];
903  memset(etaEneA,0,sizeof(etaEneA));
904 
905  int ix,iy;
906  // pick 5x5 eta-phi patch and calculate weighted phi==Y
907  for(iy=iphi0;iy<iphi0+cl2jet_par_mxPhiBin;iy++)
908  {
909  int jy=iy % cl2jetMaxPhiBins;
910  float *patchEneA=eve_patchEne+(jy*cl2jetMaxEtaBins);// phi bins are consecutive
911 
912  sum1=0.0;
913  for(ix=ieta0;ix<ieta0+cl2jet_par_mxEtaBin;ix++)
914  {
915  sum1+=patchEneA[ix];
916  etaEneA[ix-ieta0]+=patchEneA[ix];
917  }
918  if(sum1>0) nEnePhi++;
919  sumY+=sum1*iy;
920  sumYY+=sum1*iy*iy;
921  }
922 
923  // pick 5x5 eta-phi patch and calculate weighted eta==X
924  int jx;
925  for(jx=ieta0;jx<ieta0+cl2jet_par_mxEtaBin;jx++)
926  {
927  if(etaEneA[jx-ieta0]>0) nEneEta++;
928  sumX+=jx*etaEneA[jx-ieta0];
929  sumXX+=(jx*jx*etaEneA[jx-ieta0]);
930  }
931 
932  //calculate phi mean and rms
933  float phiMean= 1.*sumY/iene0;//was phiSum
934  float fphiBinMax=0.5 +phiMean;
935  float rmsPhi=0.0; // for single bin RMS is 0
936  if(nEnePhi>1)
937  { //here I calculate RMS
938  rmsPhi=sqrt(nEnePhi/(nEnePhi-1.) *(sumYY/iene0 - phiMean*phiMean));
939  }
940  if(fphiBinMax>cl2jetMaxPhiBins) fphiBinMax-=cl2jetMaxPhiBins;
941 
942  //calculate eta mean and rms
943  float etaMean = 1.*sumX/iene0;
944  float fetaBin=0.5 + etaMean;
945  float rmsEta=0.0; // for single bin RMS is 0
946  if(nEneEta>1)
947  { //here I calculate RMS
948  rmsEta=sqrt(nEneEta/(nEneEta-1.) *(sumXX/iene0 - etaMean*etaMean));
949  }
950 
951  J->fphiBin=fphiBinMax;
952  J->fetaBin=fetaBin;
953  J->rmsPhiBin=rmsPhi;
954  J->rmsEtaBin=rmsEta;
955 
956 }
957 
958 //========================================
959 //========================================
960 void
961 L2jetAlgo2009:: dumpPatchEneA(){
962  // dump L2 array with energy
963  int ix,iy;
964  for(iy=0;iy<cl2jetMaxPhiBins;iy++) {
965  float *patchEneA=eve_patchEne+(iy*cl2jetMaxEtaBins);// phi bins are consecutive
966 
967  for(ix=0;ix<cl2jetMaxEtaBins;ix++,patchEneA++){
968  printf(" %6f",*patchEneA);
969  }
970  printf(" iPhi=%d\n",iy);
971  }
972 
973 }
974 
975 
976 //=======================================
977 //=======================================
978 void
979 L2jetAlgo2009::computeE(int token){
980 }
981 //=======================================
982 //=======================================
983 void
984 L2jetAlgo2009::finishRunHisto(){
985  // auxialiary operations on histograms at the end of the run
986 
987  const int *data20=hA[20]->getData();
988  const int *data30=hA[30]->getData();
989 
990  int bHotSum=1,bHotId=-1;
991  int eHotSum=1;
992 
993  const L2EmcDb::EmcCDbItem *xE=mDb->getByIndex(402), *xB=mDb->getByIndex(402);
994  int i;
995  for(i=0; i<EmcDbIndexMax; i++) {
996  const L2EmcDb::EmcCDbItem *x=mDb->getByIndex(i);
997  if(mDb->isEmpty(x)) continue;
998  if (mDb->isBTOW(x) ) {
999  int softId=atoi(x->tube+2);
1000  int ieta= (x->eta-1);
1001  int iphi= (x->sec-1)*10 + x->sub-'a' ;
1002  //mDb->printItem(x); printf("softID=%d\n",softId);
1003  hA[21]->fillW(softId,data20[x->rdo]);
1004  hA[22]->fillW(ieta, iphi,data20[x->rdo]);
1005  if(bHotSum<data20[x->rdo]) {
1006  bHotSum=data20[x->rdo];
1007  bHotId=softId;
1008  xB=x;
1009  }
1010  }// end of BTOW
1011  else if (mDb->isETOW(x) ) {
1012  int ihard=x->chan+(x->crate-1)*128;
1013  int ieta= 12-x->eta;
1014  int iphi= (x->sec-1)*5 + x->sub-'A' ;
1015  hA[31]->fillW(ihard,data30[x->rdo]);
1016  hA[32]->fillW(ieta, iphi,data30[x->rdo]);
1017  if(eHotSum<data30[x->rdo]) {
1018  eHotSum=data30[x->rdo];
1019  xE=x;
1020  }
1021 
1022  }// end of BTOW
1023  }
1024  if (mLogFile){
1025  fprintf(mLogFile,"L2jet::finishRun()\n");
1026  fprintf(mLogFile,"#BTOW_hot tower _candidate_ (bHotSum=%d) :, softID %d , crate %d , chan %d , name %s\n",bHotSum,bHotId,xB->crate,xB->chan,xB->name);
1027  fprintf(mLogFile,"#ETOW_hot tower _candidate_ (eHotSum=%d) :, name %s , crate %d , chan %d\n",eHotSum,xE->name,xE->crate,xE->chan);
1028  }
1029 
1030 }
1031 
1032 
1033 /**********************************************************************
1034  $Log: L2jetAlgo2009.cxx,v $
1035  Revision 1.1 2010/04/17 16:42:14 pibero
1036  *** empty log message ***
1037 
1038  Revision 1.1 2007/11/19 22:18:27 balewski
1039  most L2algos provide triggerID's
1040 
1041  Revision 1.8 2007/11/14 03:58:14 balewski
1042  cleanup of common timing measurement
1043 
1044  Revision 1.7 2007/11/13 23:06:07 balewski
1045  toward more unified L2-algos
1046 
1047  Revision 1.6 2007/11/13 00:12:36 balewski
1048  added offline triggerID, take1
1049 
1050  Revision 1.5 2007/11/08 04:02:31 balewski
1051  run on l2ana as well
1052 
1053  Revision 1.4 2007/11/02 17:43:08 balewski
1054  cleanup & it started to work w/ L2upsilon
1055 
1056  Revision 1.3 2007/11/02 03:03:47 balewski
1057  modified L2VirtualAlgo
1058 
1059  Revision 1.2 2007/10/25 02:07:02 balewski
1060  added L2upsilon & binary event dump
1061 
1062  Revision 1.1 2007/10/11 00:33:19 balewski
1063  L2algo added
1064 
1065  Revision 1.5 2006/03/28 19:46:49 balewski
1066  ver16b, in l2new
1067 
1068  Revision 1.4 2006/03/11 17:08:33 balewski
1069  now CVS comments should work
1070 
1071 */
1072 
1073 
Generated by   doxygen 1.8.5