StRoot: StTriggerUtilities/L2Emulator/L2jetAlgo/L2jetAlgo2012.cxx Source File

00001 #include <stdio.h>
00002 #include <string.h>
00003 #include <stdlib.h>
00004 #include <time.h>
00005 #include <math.h>
00006 #include <fakeRtsLog.h>
00007 
00008 /***********************************************************
00009  * $Id: L2jetAlgo2012.cxx,v 1.4 2012/03/21 18:18:03 jml Exp $
00010  * \author Jan Balewski, IUCF, 2006 
00011  ***********************************************************
00012  * Descripion:
00013   Reco of mono- & di-jets in L2 using BTOW+ETOW
00014   depends on L2-DB class 
00015  ***********************************************************
00016  * version updated 2/17/2006 W.W. Jacobs
00017  * 2008 to 2009 change by B.S. Page 11/12/2008
00018  */
00019 
00020 
00021 
00022 #ifdef  IS_REAL_L2  //in l2-ana  environment
00023   #include "../L2algoUtil/L2EmcDb2012.h"
00024   #include "../L2algoUtil/L2Histo.h"
00025 #else
00026   #include "StTriggerUtilities/L2Emulator/L2algoUtil/L2EmcDb2012.h"
00027   #include "StTriggerUtilities/L2Emulator/L2algoUtil/L2Histo.h"
00028 #endif
00029 
00030 #include "L2jetAlgo2012.h"
00031 #include "L2jetResults2012.h" 
00032 #include "Map_DeltaPhiJets.h"
00033 
00034 //=================================================
00035 //=================================================
00036 L2jetAlgo2012::L2jetAlgo2012(const char* name, const char *uid, L2EmcDb2012* db, char* outDir, int resOff, bool writeHighResult) 
00037   :  L2VirtualAlgo2012( name, uid,  db,  outDir, true, true, resOff) { 
00038   /* called one per days
00039      all memory allocation must be done here
00040   */
00041   namechar='L';
00042   if (writeHighResult) namechar='H';
00043 
00044   createHisto();
00045   run_number=-1;
00046   LOG(DBG,"L2jetAlgo2012 instantiated, logPath='%s'\n",mOutDir1.c_str());
00047   eve_Jet[0]=new L2Jet;
00048   eve_Jet[1]=new L2Jet;
00049 }
00050 
00051 /*========================================
00052   ======================================== */
00053 bool
00054 L2jetAlgo2012::paramsChanged( int *rc_ints, float *rc_floats) {
00055   int i;
00056   for(i=0;i<5;i++)
00057     if(rc_ints[i]!=raw_ints[i]) goto  foundProblem;
00058 
00059   for(i=0;i<5;i++)
00060     if(fabs(rc_floats[i]-raw_floats[i])>0.00001)  goto  foundProblem;
00061 
00062   return false;
00063 
00064  foundProblem:
00065       if (mLogFile) fprintf(mLogFile,"L2jet-ghost initRun - inconsistent params, ABORT initialization\n");
00066   return true;
00067 }
00068 
00069 /*========================================
00070   ======================================== */
00071 int 
00072 L2jetAlgo2012::initRunUser( int runNo, int *rc_ints, float *rc_floats) {
00073 
00074   // clear content of all histograms
00075   int i;
00076   for (i=0; i<mxHA;i++) if(hA[i])hA[i]->reset();
00077 
00078   /* .... clear content, set threshold @ max as default */
00079   memset(db_btowL2PhiBin,  0,  sizeof(db_btowL2PhiBin));
00080   memset(db_btowL2PatchBin,0,  sizeof(db_btowL2PatchBin));
00081 
00082   memset(db_etowL2PhiBin,  0  ,sizeof(db_etowL2PhiBin));
00083   memset(db_etowL2PatchBin,0  ,sizeof(db_etowL2PatchBin));
00084 
00085   run_startUnix=time(0);
00086   run_number  =runNo;  // serves as a flag this run is initialized
00087   raw_ints    =rc_ints;
00088   raw_floats  =rc_floats;
00089   run_nEventOneJet=run_nEventDiJet= run_nEventRnd=0;
00090 
00091   mEventsInRun=0;
00092 
00093   if( mLogFile==0) LOG(ERR," L2jetAlgo2012() UNABLE to open run summary log file, continue anyhow\n");
00094 
00095   // unpack params from run control GUI
00096   par_useBtowEast= (rc_ints[0]&1 )>0;
00097   par_useBtowWest= (rc_ints[0]&2)>0;
00098   par_useEndcap  =  rc_ints[1]&1;
00099   int par_adcThr =  rc_ints[2]; // needed only in initRun()
00100   par_dbg = rc_ints[3];
00101   par_RndAcceptPrescale = rc_ints[4];
00102   par_minPhiBinDiff= 7; //for 1X1 JP was 5 for .6x.6 JP BP
00103 
00104   par_oneJetThr   =  rc_floats[0];
00105   par_diJetThrHigh=  rc_floats[1];
00106   par_diJetThrLow =  rc_floats[2];
00107   par_diJetThr_2 =   rc_floats[3];  // wwj 2/11
00108   par_diJetThr_3 =   rc_floats[4];  // wwj 2/11
00109   par_diJetThr_4 =   rc_floats[5];  // wwj 2/11
00110   par_diJetThr_5 =   rc_floats[6];  // wwj 2/11
00111   par_diJetEtaHigh =  rc_floats[7];  // wwj 2/11
00112   par_diJetEtaLow = rc_floats[8];  // wwj 2/11
00113 
00114   // set other hardcoded or calculated  params
00115 
00116   // to monitor hot towers  in E+B Emc
00117   float monTwEtThr=2.0; // (GeV) Et threshold,WARN,  edit histo title by hand
00118   par_hotTwEtThr = monTwEtThr;
00119 
00120   if(par_dbg) LOG(DBG,"Brian version is running!\n");  
00121 
00122   if (mLogFile) { 
00123     fprintf(mLogFile,"L2jet algorithm initRun(%d), compiled: %s , %s\n params:\n",run_number,__DATE__,__TIME__);
00124     fprintf(mLogFile," - use BTOW: East=%d West=%d, Endcap=%d  L2ResOffset=%d\n", par_useBtowEast, par_useBtowWest,par_useEndcap ,mResultOffset);
00125     fprintf(mLogFile," - threshold: ADC-ped> %d \n", par_adcThr);
00126     fprintf(mLogFile," - min phi opening angle Jet1<->Jet2: %d in L2phiBins\n",par_minPhiBinDiff);   
00127     fprintf(mLogFile," - diJet  Et thrHigh= %.2f (GeV)   thrLow= %.2f  (GeV)\n", par_diJetThrHigh, par_diJetThrLow); 
00128     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 
00129     fprintf(mLogFile," - oneJet Et thr = %.2f (GeV) ; rndAccPrescale=%d \n",par_oneJetThr,par_RndAcceptPrescale);
00130     fprintf(mLogFile," - new diJet Eta_low= %.2f ; diJet Eta_High= %.2f \n", par_diJetEtaLow, par_diJetEtaHigh); // wwj 2/10
00131     fprintf(mLogFile," - debug=%d, hot tower threshold: Et> %.1f GeV ( only monitoring)\n",par_dbg, monTwEtThr);
00132   }
00133 
00134   // verify consistency of input params
00135   int kBad=0;
00136   kBad+=0x0001 * ( !par_useBtowEast & !par_useBtowWest & !par_useEndcap);
00137   // par_pedOff no longer exists:  kBad+=0x0002 * (par_adcThr<par_pedOff); 
00138   kBad+=0x0004 * (par_adcThr>16); 
00139   kBad+=0x0008 * (par_minPhiBinDiff<5);
00140   kBad+=0x0010 * (par_minPhiBinDiff>=15);
00141   kBad+=0x0020 * (par_oneJetThr<3.);
00142   kBad+=0x0040 * (par_oneJetThr>12.);
00143   kBad+=0x0080 * (par_diJetThrLow<1.9);  //     LOOK: wwj 2/11  using lowest threshold 
00144   kBad+=0x0100 * (par_diJetThrHigh<par_diJetThrLow); 
00145   kBad+=0x0200 * (par_diJetThrHigh>12.); 
00146   //0x0400 not used
00147   kBad+=0x0800 * ( par_RndAcceptPrescale<0 );
00148   kBad+=0x1000 * (par_diJetThr_2<par_diJetThrLow || par_diJetThr_5>par_diJetThrHigh); // wwj 2/16 
00149   kBad+=0x2000 * (par_diJetThr_2>par_diJetThr_3 || par_diJetThr_3>par_diJetThr_4 || par_diJetThr_4>par_diJetThr_5); // wwj 2/16
00150   kBad+=0x4000 * (par_diJetEtaLow>par_diJetEtaHigh); // wwj 2/16
00151   if (mLogFile) {
00152     fprintf(mLogFile,"L2jet initRun() params checked for consistency, Error flag=0x%04x\n",kBad);
00153     if(kBad)   fprintf(mLogFile,"L2jet initRun()  ABORT\n");
00154   }
00155 
00156   if(kBad) {
00157     run_number=-66;
00158     if (mLogFile) { 
00159       fprintf(mLogFile,"L2jet algorithm init: crashB in internal logic\n");
00160       fclose(mLogFile);
00161       return kBad;
00162     }
00163   }
00164 
00165   char tit[100];
00166   sprintf(tit,"# BTOW towers>ped+%d (input); x: # of towers/event",par_adcThr);
00167   hA[47]->setTitle(tit);
00168   
00169   sprintf(tit,"# ETOW towers>ped+%d (input); x: # of towers/event",par_adcThr);
00170   hA[48]->setTitle(tit);
00171 
00172  
00173   const int mxEtaBinsE=12,mxEtaBinsB=40;
00174 
00175   // rebuild local lookup tables
00176   
00177   int etowEtaBin2Patch[mxEtaBinsE]={14,14,13,13,12,12,11,11,11,10,10,10};
00178 
00179   int nB=0, nE=0; /* counts # of unmasekd towers */ 
00180   int nBg=0, nEg=0; /* counts # of reasonable calibrated towers */ 
00181 
00182   for(i=0; i<EmcDbIndexMax; i++) {
00183     const L2EmcDb2012::EmcCDbItem *x=mDb->getByIndex(i);
00184     if(mDb->isEmpty(x)) continue;  /* dropped not mapped  channels */
00185     if(x->fail) continue; /* dropped masked channels */
00186     if(x->gain<=0) continue; /* dropped uncalibrated towers , tmp */
00187     /* if(x->sec!=1) continue;   tmp, to test patch mapping */
00188 
00189     /* WARN, calculate index to match RDO order in the ADC data banks */
00190     int ietaP, iphiP;
00191     if (mDb->isBTOW(x) ) {
00192       /*....... B A R R E L  .................*/
00193       nB++;   
00194       if(x->eta<0 || x->eta>mxEtaBinsB) goto crashIt_1;
00195       if(!par_useBtowEast && x->eta<=20) continue;
00196       if(!par_useBtowWest && x->eta>=21) continue;
00197       ietaP= (x->eta-1)/4; /* correct */
00198       int iphiTw=(x->sec-1)*10 + x->sub-'a';
00199       // allign in phi  TP @ L2 w/ L0 
00200       iphiTw--;
00201       if(iphiTw<0) iphiTw=119;
00202       // now cr0x1e, mod1, subm2, is beginning of the first BTOW TP
00203       iphiP= iphiTw/4 ; /* correct */
00204       //      if(ietaP==0 && iphiP==5) 
00205       db_btowL2PhiBin[x->rdo]=iphiP;
00206       db_btowL2PatchBin[x->rdo]=ietaP+ iphiP*cl2jetMaxEtaBins;
00207       nBg++;
00208     } else if(mDb->isETOW(x) &&  par_useEndcap) {
00209       /*....... E N D C A P ........................*/
00210       nE++;   
00211       int iphiTw= (x->sec-1)*5 + x->sub-'A';
00212       // allign in phi  TP @ L2 w/ L0 
00213       iphiTw--;
00214       if(iphiTw<0) iphiTw=59;
00215       // now subsector 01TB is beginning of the first(==0) ETOW TP in phi
00216       iphiP= iphiTw/2 ; /* correct */
00217       if(x->eta<0 || x->eta>mxEtaBinsE) goto crashIt_1;
00218       ietaP=etowEtaBin2Patch[x->eta-1];
00219       db_etowL2PhiBin[x->rdo]=iphiP;
00220       db_etowL2PatchBin[x->rdo]=ietaP+ iphiP*cl2jetMaxEtaBins;
00221       nEg++;
00222     }
00223     
00224   }
00225 
00226   if (mLogFile) {
00227     fprintf(mLogFile,"L2jet algorithm: found  working/calibrated: %d/%d=ETOW  & %d/%d=BTOW, based on ASCII DB\n",nE,nEg,nB,nBg);
00228   }
00229 
00230   return 0; //OK
00231   
00232  crashIt_1: /* fatal initialization error */
00233   run_number=-55;
00234   if (mLogFile) { 
00235     fprintf(mLogFile,"L2jet algorithm init: crashC in internal logic\n");
00236     fclose(mLogFile);
00237   }
00238  return -6;
00239 
00240 }               
00241 
00242 /*========================================
00243   ======================================== */
00244 void L2jetAlgo2012::computeUser(int token){
00245   //put in to satisfy l2virtualalgo2012
00246 }
00247 
00248 /*========================================
00249   ======================================== */
00250 bool
00251 L2jetAlgo2012::decisionUser(int token, int *myL2Result){
00252   /* STRICT TIME BUDGET  START ....*/
00253   
00254   //Ross says:  Eventually you should take these out - timing is done in the virtual algo class as well.
00255   unsigned long mEveTimeStart,mEveTimeStop, mEveTimeDiff;
00256 rdtscl_macro(mEveTimeStart);
00257 
00258   //Ross says:  I will check how the following line is done in L2bemcGamma2012
00259   //following if statement was eve_ID!=inpEveId BP
00260   if(1) {//UUUU
00261     //.................... this event has NOT been processed
00262 
00263   /*
00264     Chris doesn't want us to write anything out 
00265     during event processing ...
00266   */
00267 
00268     bool bemcIn = true;
00269     bool eemcIn = true;
00270   
00271     //eve_ID=inpEveId; // every events is processed only once
00272     mEventsInRun++;
00273     clearEvent(); /* price=13 kTicks */  
00274     int runTimeSec=time(0)- run_startUnix;
00275     hA[10]->fill(0);
00276     hA[12]->fill(runTimeSec);
00277     
00278     if(par_dbg>1) LOG(DBG,"\n......... in  L2Jet_doEvent(ID=%d)... bIn=%d eIn=%d\n",eve_ID,bemcIn,eemcIn);//bemcIn and eemcIn were passed to the old doEvent
00279   
00280 
00281     //Ross asks: will the jet code run on just barrel or just endcap?
00282     //  if that's the case, will it save time to move those sections of the
00283     //  code to computeUser() rather than decisionUser()? 
00284     if (bemcIn || eemcIn){//VVVV this algo has nothing to do w/o any Ecal data 
00285       
00286   //===== step 1: unpack raw data blocks ==============
00287    
00288   int nBtowTw=0, nEtowTw=0;
00289   //Ross says:  I've changed projectAdc to use the L2eventStream2012 class
00290   /*......... BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB */
00291 
00292   //we can move this to compute.
00293   if(bemcIn==1 && (par_useBtowEast||par_useBtowWest) ) {
00294 
00295     nBtowTw=projectAdc( mEveStream_btow[token].get_hits(), mEveStream_btow[token].get_hitSize(),
00296                         db_btowL2PhiBin, db_btowL2PatchBin,
00297                         hA[20] );
00298   }
00299 
00300   /*........... EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE */
00301   //and this as well.
00302   if(eemcIn==1 && par_useEndcap ) {
00303     nEtowTw=projectAdc(  mEveStream_etow[token].get_hits(), mEveStream_etow[token].get_hitSize(),
00304                         db_etowL2PhiBin, db_etowL2PatchBin,
00305                         hA[30] );
00306   }
00307 
00308   //============Do 2D scan==========
00309   float totEneGeV=true2Dscan();
00310   int itotEneGeV=(int)totEneGeV;
00311   int iK;
00312   for(iK=0; iK<mxJ; iK++)
00313     {
00314       weightedEtaPhi(iK);
00315       L2Jet *K=eve_Jet[iK];
00316       K->eneGeV=K->iene;
00317       K->phiRad=0.21*(6.0-K->fphiBin);
00318       while(K->phiRad<0) K->phiRad+=6.2832;
00319       while(K->phiRad>6.2832) K->phiRad-=6.2832;
00320     }
00321 
00322   if(eve_Jet[0]->eneGeV <eve_Jet[1]->eneGeV) {// swap jets fo E1>E2
00323     L2Jet *Jx=eve_Jet[0];
00324     eve_Jet[0]=eve_Jet[1];
00325     eve_Jet[1]=Jx;
00326   }
00327 
00328   if(par_dbg>2) LOG(DBG,"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);
00329   
00330   //====== step 4: make trigger decisions====
00331   
00332   // now put in eta vs. Jet Threshold "contour" cuts ... this vesion cuts on eta_1 + eta_2 .. wwj 2/16 
00333   // looks like fetaBin goes 0.0->15.0 for eta -1. -> +2. ... so "eta" = fetaBin*0.2 -1.0 ... wwj 2/10
00334 
00335   float rjetEta_0=eve_Jet[0]->fetaBin*0.2 -1.0; //  wwj 2/10 
00336   float rjetEta_1=eve_Jet[1]->fetaBin*0.2 -1.0; //  wwj 2/10 
00337 
00338   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
00339 
00340   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
00341   
00342   
00343   // bool acceptDiJet=( eve_Jet[0]->eneGeV > par_diJetThrHigh) && ( eve_Jet[1]->eneGeV > par_diJetThrLow);
00344   bool acceptDiJet_BB=(eve_Jet[0]->eneGeV > par_diJetThrHigh && eve_Jet[1]->eneGeV > par_diJetThr_3); // wwj 2/10 
00345 
00346   bool acceptDiJet=acceptDiJet_EE || acceptDiJet_EB || acceptDiJet_BB; // wwj 2/10
00347 
00348   bool acceptOneJet=( eve_Jet[0]->eneGeV > par_oneJetThr) ;
00349 
00350   bool acceptRnd=mRandomAccept>0; // provided by Virtual algo
00351   mAccept=acceptDiJet || acceptOneJet || acceptRnd;
00352  
00353   //====== step 5: update various monitorig histos
00354 
00355   // histogramming reco Et1-Et2, no cuts
00356   int iet1 =(int)eve_Jet[0]->eneGeV;
00357   int iet2 =(int)eve_Jet[1]->eneGeV;
00358   int ieta1=(int)eve_Jet[0]->fetaBin;
00359   int ieta2=(int)eve_Jet[1]->fetaBin;
00360   int iphi1=(int)eve_Jet[0]->fphiBin;
00361   int iphi2=(int)eve_Jet[1]->fphiBin;
00362 
00363   hA[40]->fill(iet1,iet2);
00364 
00365   hA[41]->fill(ieta1,iphi1);
00366   hA[42]->fill(ieta2,iphi2);
00367   hA[43]->fill(iphi1,iphi2);
00368   hA[44]->fill(iet1);
00369   hA[45]->fill(iet2);
00370   hA[46]->fill(itotEneGeV);
00371   hA[47]->fill(nBtowTw);
00372   hA[48]->fill(nEtowTw);
00373 
00374   // sivers delta zeta, the map is still worng
00375   int kphi1=int(eve_Jet[0]->phiRad*10.);
00376   int kphi2=int(eve_Jet[1]->phiRad*10.);
00377   int idelZeta=map_DelPhiJets[kphi1*MxPhiRad10 + kphi2];
00378 
00379   if( mAccept)  hA[10]->fill(8);
00380 
00381   if(acceptOneJet ){ 
00382     hA[10]->fill(4);
00383     run_nEventOneJet++;
00384     hA[13]->fill(runTimeSec);
00385     hA[50]->fill(iet1);
00386     hA[51]->fill(ieta1,iphi1);
00387     hA[52]->fill(ieta1);
00388     hA[53]->fill(iphi1);
00389   }
00390   
00391   if(acceptDiJet  ){
00392     hA[10]->fill(5);
00393     run_nEventDiJet++;
00394     hA[14]->fill(runTimeSec);
00395     hA[60]->fill(iet1,iet2);
00396     hA[61]->fill(ieta1,iphi1);
00397     hA[62]->fill(ieta2,iphi2);
00398     hA[63]->fill(iphi1,iphi2);
00399     hA[64]->fill(iet1);
00400     hA[65]->fill(iet2);
00401     hA[66]->fill(ieta1);
00402     hA[67]->fill(ieta2);
00403     hA[68]->fill(iphi1);
00404     hA[69]->fill(iphi2);
00405     hA[70]->fill(idelZeta);
00406     hA[71]->fill(ieta1,idelZeta);
00407     hA[72]->fill(ieta1,ieta2);
00408     hA[73]->fill((iphi1+iphi2)/2,idelZeta);
00409     hA[74]->fill(itotEneGeV);
00410    }
00411 
00412   if(acceptDiJet_EE  ){
00413     run_nEventDiJet++;        // all these histos #80-#94 added ... wwj 2/10
00414     hA[80]->fill(iet1,iet2);
00415     hA[81]->fill(ieta1,iphi1);
00416     hA[82]->fill(ieta2,iphi2);
00417     hA[83]->fill(iphi1,iphi2);
00418     hA[84]->fill(iet1);
00419     hA[85]->fill(iet2);
00420     hA[86]->fill(ieta1);
00421     hA[87]->fill(ieta2);
00422     hA[88]->fill(iphi1);
00423     hA[89]->fill(iphi2);
00424     hA[90]->fill(idelZeta);
00425     hA[91]->fill(ieta1,idelZeta);
00426     hA[92]->fill(ieta1,ieta2);
00427     hA[93]->fill((iphi1+iphi2)/2,idelZeta);
00428     hA[94]->fill(itotEneGeV);
00429    }
00430 
00431   if(acceptDiJet_EB  ){
00432     run_nEventDiJet++;        // all these histos #100-#114 added ... wwj 2/10
00433     hA[100]->fill(iet1,iet2);
00434     hA[101]->fill(ieta1,iphi1);
00435     hA[102]->fill(ieta2,iphi2);
00436     hA[103]->fill(iphi1,iphi2);
00437     hA[104]->fill(iet1);
00438     hA[105]->fill(iet2);
00439     hA[106]->fill(ieta1);
00440     hA[107]->fill(ieta2);
00441     hA[108]->fill(iphi1);
00442     hA[109]->fill(iphi2);
00443     hA[110]->fill(idelZeta);
00444     hA[111]->fill(ieta1,idelZeta);
00445     hA[112]->fill(ieta1,ieta2);
00446     hA[113]->fill((iphi1+iphi2)/2,idelZeta);
00447     hA[114]->fill(itotEneGeV);
00448    }
00449 
00450   if(acceptDiJet_BB  ){
00451     run_nEventDiJet++;        // all these histos #120-#134 added ... wwj 2/10
00452     hA[120]->fill(iet1,iet2);
00453     hA[121]->fill(ieta1,iphi1);
00454     hA[122]->fill(ieta2,iphi2);
00455     hA[123]->fill(iphi1,iphi2);
00456     hA[124]->fill(iet1);
00457     hA[125]->fill(iet2);
00458     hA[126]->fill(ieta1);
00459     hA[127]->fill(ieta2);
00460     hA[128]->fill(iphi1);
00461     hA[129]->fill(iphi2);
00462     hA[130]->fill(idelZeta);
00463     hA[131]->fill(ieta1,idelZeta);
00464     hA[132]->fill(ieta1,ieta2);
00465     hA[133]->fill((iphi1+iphi2)/2,idelZeta);
00466     hA[134]->fill(itotEneGeV);
00467    }
00468  
00469   if(mRandomAccept){  
00470     hA[10]->fill(6);
00471     run_nEventRnd++;
00472     hA[15]->fill(runTimeSec);
00473   }
00474   
00475   
00476   //====== step 6:  fill L2Result  (except time)
00477   L2jetResults2012 out; // all output bits lump together
00478   memset(&out,0,sizeof(out)); // clear content
00479   
00480   
00481   // add acceptDiJet_EE, acceptDiJet_EB and acceptDiJet_BB to below ... wwj 2/10
00482   out.int0.version=L2JET_RESULTS_VERSION;
00483   out.int0.decision=
00484     ( par_useBtowEast <<0 ) +
00485     ( par_useBtowWest <<1 ) +
00486     ( par_useEndcap   <<2 ) +
00487     ( (bemcIn>0)      <<3 ) + 
00488     ( (eemcIn>0)      <<4 ) + 
00489     (  mRandomAccept  <<5 ) +
00490     (  acceptOneJet   <<6 ) +
00491     (  acceptDiJet    <<7 ) +
00492     ( acceptDiJet_EE  <<8 ) +
00493     ( acceptDiJet_EB  <<9 ) +
00494     ( acceptDiJet_BB <<10 ) ;  // e.g., added last three lines here  ... wwj 2/10
00495   out.int0.dumm=namechar;
00496 
00497   out.int1.iTotEne=(unsigned short)(totEneGeV*100.); // now 1=10 MeV
00498   out.int2.nBtowTw=nBtowTw;
00499   out.int2.nEtowTw=nEtowTw;
00500   
00501   out.jet1.jPhi=(int)(eve_Jet[0]->phiRad*28.65);  //so phi/deg=2*jPhi
00502   out.jet1.jEta=(int)(eve_Jet[0]->fetaBin*10.);
00503   out.jet1.iEne=(unsigned short)(eve_Jet[0]->eneGeV*100.); // now 1=10 MeV
00504   
00505   out.jet2.jPhi=(int)(eve_Jet[1]->phiRad*28.65);
00506   out.jet2.jEta=(int)(eve_Jet[1]->fetaBin*10.);
00507   out.jet2.iEne=(unsigned short)(eve_Jet[1]->eneGeV*100.); // now 1=10 MeV
00508 
00509   out.jet1.rmsEta=(unsigned short)(eve_Jet[0]->rmsEtaBin*20.); //Fix: 1=0.01 physical eta
00510   out.jet1.rmsPhi=(unsigned short)(eve_Jet[0]->rmsPhiBin*120.); //1=0.1 deg
00511 
00512   out.jet2.rmsEta=(unsigned short)(eve_Jet[1]->rmsEtaBin*20.); //Fix: 1=0.01 physical eta
00513   out.jet2.rmsPhi=(unsigned short)(eve_Jet[1]->rmsPhiBin*120.);// 1=0.1 deg
00514 
00515   
00516   rdtscl_macro(mEveTimeStop);
00517   mEveTimeDiff=mEveTimeStop-mEveTimeStart;
00518   int  kTick=mEveTimeDiff/1000;
00519 
00520   //mhT->fill(kTick); taken out 12/11/08 BP
00521 
00522   out.int0.kTick=  kTick>255 ? 255 : kTick;
00523  
00524   //calculate and match the check sum
00525   out.int1.checkSum=-L2jetResults2012_doCheckSum(&out);
00526 
00527   
00528 
00529   //===== step 7: write L2Result
00530   memcpy(myL2Result,&out,sizeof(L2jetResults2012));
00531 
00532   // dirty tests, clean it up before real use
00533   
00534   if(par_dbg){//WWWW
00535     L2jetResults2012_print(&out);
00536     LOG(DBG," phiRad1=%f  phiRad2=%f \n",eve_Jet[0]->phiRad,eve_Jet[1]->phiRad);
00537     LOG(DBG,"idelZeta=%d  delZeta/deg=%.1f \n\n",idelZeta,idelZeta/31.416*180);
00538     
00539     //tmp printouts of errors:
00540     if( out.jet1.iEne+out.jet2.iEne > out.int1.iTotEne) {
00541     }
00542     if(iphi1==iphi2) {
00543       LOG(ERR,"L2jet-fatal error,neveId=%d, phi1,2=%d,%d\n",mEventsInRun,iphi1,iphi2);
00544       dumpPatchEneA();    
00545     }
00546     
00547     if( L2jetResults2012_doCheckSum(&out)) {
00548       LOG(ERR,"L2jet-fatal error, wrong cSum=%d\n", L2jetResults2012_doCheckSum(&out));
00549       L2jetResults2012_print(&out);
00550     }
00551     
00552   } // end of WWWW
00553   
00554     }// end of VVVV (etow or btow has some data)
00555   }// end of UUUU event processing
00556   
00557   
00558   return   mAccept;
00559 }
00560 
00561 
00562 /*========================================
00563   ======================================== */
00564 void 
00565 L2jetAlgo2012::finishRunUser() {  /* called once at the end of the run */
00566   if(run_number<0) return; // already finished
00567 
00568   if (mLogFile) {
00569     fprintf(mLogFile,"L2-jet algorithm finishRun(%d)\n",run_number);
00570     fprintf(mLogFile," - %d events seen by L2 di-jet\n",mEventsInRun);
00571     fprintf(mLogFile," - accepted: rnd=%d  oneJet=%d diJet=%d \n", run_nEventRnd,  run_nEventOneJet, run_nEventDiJet);
00572 
00573     // print few basic histos
00574     hA[10]->printCSV(mLogFile); // event accumulated
00575 
00576   }
00577   finishRunHisto(); // still needs current DB
00578   if( mHistFile==0) {
00579     LOG(ERR," L2jetAlgo2012: finishRun() UNABLE to open run summary log file, continue anyhow\n");
00580     if (mLogFile)
00581       fprintf(mLogFile,"L2 di-jet histos NOT saved, I/O error\n");
00582   } else { // save histos  
00583     int j;
00584     int nh=0;
00585     for(j=0;j<mxHA;j++) {
00586       if(hA[j]==0) continue;
00587       hA[j]->write(mHistFile);
00588       nh++;
00589     }    
00590   }
00591   
00592   run_number=-2; // clear run #
00593  
00594 }
00595 
00596 //=======================================
00597 //=======================================
00598 void 
00599 L2jetAlgo2012::createHisto() {
00600   memset(hA,0,sizeof(hA));
00601 
00602   hA[10]=new   L2Histo(10,"total event counter; x=cases",9);
00603   hA[11]=new   L2Histo(11,"L2 time used per input event;  x: time (CPU kTics), range=100muSec; y: events ",160);
00604 
00605   int mxRunDration=2500;
00606   hA[12]=new   L2Histo(12,"rate of input events; x: time in this run (seconds); y: rate (Hz)", mxRunDration);
00607   
00608   hA[13]=new   L2Histo(13,"rate of  accepted one-Jet; x: time in this run (seconds); y: rate (Hz)", mxRunDration);
00609   hA[14]=new   L2Histo(14,"rate of  accepted di-Jet ; x: time in this run (seconds); y: rate (Hz)", mxRunDration);
00610   hA[15]=new   L2Histo(15,"rate of  random accepted  ; x: time in this run (seconds); y: rate (Hz)", mxRunDration);
00611 
00612   // BTOW  raw spectra
00613   hA[20]=new   L2Histo(20,"BTOW tower, Et>2.0 GeV (input); x: BTOW RDO index=chan*30+fiber; y: counts", 4800);
00614   hA[21]=new   L2Histo(21,"BTOW tower, Et>2.0 GeV (input); x: BTOW softID", 4800);
00615   hA[22]=new   L2Histo(22,"BTOW tower, Et>2.0 GeV (input); x: eta bin, [-1,+1];  y: phi bin ~sector",40,120);
00616   
00617   // ETOW  raw spectra
00618   hA[30]=new   L2Histo(30,"ETOW tower, Et>2.0 GeV (input); x: ETOW RDO index=chan*6+fiber; y: counts", 720 );
00619   hA[31]=new   L2Histo(31,"ETOW tower, Et>2.0 GeV (input); x: i=chan+128*crate", 768);
00620   hA[32]=new   L2Histo(32,"ETOW tower, Et>2.0 GeV (input); x: 12 - Endcap etaBin ,[+1,+2];  y: phi bin ~sector",12,60);
00621   
00622  // Di-Jet raw yields
00623  hA[40]=new   L2Histo(40,"Et Jet1-Jet2 (input); x: Jet1 Et/GeV ; Jet2 Et/GeV",30,30);
00624  hA[41]=new   L2Histo(41,"diJet1 eta-phi (input); x: iEta [-1,+2] ; y: iPhi ~sector ",15,30);
00625  hA[42]=new   L2Histo(42,"diJet2 eta-phi (input); x: iEta [-1,+2]  ; y: iPhi ~sector",15,30);
00626 
00627  hA[43]=new  L2Histo(43,"diJet phi1-phi2 (input); x: iPhi1 ~sector ; y: iPhi2 ~sector ",30,30);
00628 
00629  hA[44]=new  L2Histo(44,"Jet1 Et (input); x: Et (GeV)", 60);
00630  hA[45]=new  L2Histo(45,"Jet2 Et (input); x: Et (GeV)", 60);
00631  hA[46]=new  L2Histo(46,"total Et (input); x: Et (GeV)", 100);
00632  hA[47]=new  L2Histo(47,"# BTOW towers>thrXX (input); x: # of towers/event", 200);
00633  hA[48]=new  L2Histo(48,"# ETOW towers>thrXX (input); x: # of towers/event", 100);
00634 
00635  // ........accepted one-jet events
00636  hA[50]=new  L2Histo(50,"one-Jet Et (accepted); x: jet Et (GeV)", 60);
00637  hA[51]=new  L2Histo(51,"one-Jet eta-phi (accepted); x: iEta [-1,+2] ; y: iPhi ~sector ",15,30);
00638  hA[52]=new  L2Histo(52,"one-Jet eta (accepted); x: iEta [-1,+2]", 15);
00639  hA[53]=new  L2Histo(53,"one-Jet phi (accepted); x: iPhi ~sector", 30);
00640 
00641  // Di-Jet accepted
00642  hA[60]=new   L2Histo(60,"Et of Jet1 vs. Jet2  (accepted); x: Jet1/GeV ; Jet2/GeV",30,30);
00643  hA[61]=new   L2Histo(61,"diJet1 eta-phi   (accepted); x: iEta [-1,+2] ; y: iPhi ~sector ",15,30);
00644  hA[62]=new   L2Histo(62,"diJet2 eta-phi   (accepted); x: iEta [-1,+2]  ; y: iPhi ~sector",15,30);
00645 
00646  hA[63]=new  L2Histo(63,"diJet phi1-phi2   (accepted); x: iPhi1 ~sector ; y: iPhi2 ~sector ",30,30);
00647 
00648  hA[64]=new  L2Histo(64,"diJet1 Et  (accepted); x: Et (GeV)", 60);
00649  hA[65]=new  L2Histo(65,"diJet2 Et   (accepted); x: Et (GeV)", 60);
00650 
00651  hA[66]=new  L2Histo(66,"diJet1  eta (accepted); x: i Eta [-1,+2]", 15);
00652  hA[67]=new  L2Histo(67,"diJet2  eta (accepted); x: i Eta [-1,+2]", 15);
00653  hA[68]=new  L2Histo(68,"diJet1 phi (accepted); x: iPhi ~sector", 30);
00654  hA[69]=new  L2Histo(69,"diJet2 phi (accepted); x: iPhi ~sector", 30);
00655  hA[70]=new  L2Histo(70,"diJet delZeta  (accepted); x: delta zeta  (rad*10)", MxPhiRad10);
00656  hA[71]=new  L2Histo(71,"diJet delZeta vs. eta1 (accepted); x: iEta1 [-1,+2] ; y: delta zeta  (rad*10)",15, MxPhiRad10);
00657  hA[72]=new  L2Histo(72,"diJet eta2 vs. eta1  (accepted); x: iEta1 [-1,+2] ;x: iEta2 [-1,+2] ",15,15);
00658  hA[73]=new  L2Histo(73,"diJet   delZeta vs. avrPhi (accepted); x: (iphi1+iphi2)/2  (12 deg/bin); y: delta zeta  (rad*10)",30, MxPhiRad10);
00659  hA[74]=new  L2Histo(74,"total Et diJet (accepted); x: Et (GeV)", 100);
00660 
00661  // Di-Jet_EE accepted  ... all these histos #80-#94 added  ... wwj 2/10
00662  hA[80]=new   L2Histo(80,"Et of Jet1 vs. Jet2  (accepted_EE); x: Jet1/GeV ; Jet2/GeV",30,30);
00663  hA[81]=new   L2Histo(81,"diJet1 eta-phi   (accepted_EE); x: iEta [-1,+2] ; y: iPhi ~sector ",15,30);
00664  hA[82]=new   L2Histo(82,"diJet2 eta-phi   (accepted_EE); x: iEta [-1,+2]  ; y: iPhi ~sector",15,30);
00665 
00666  hA[83]=new  L2Histo(83,"diJet phi1-phi2   (accepted_EE); x: iPhi1 ~sector ; y: iPhi2 ~sector ",30,30);
00667 
00668  hA[84]=new  L2Histo(84,"diJet1 Et  (accepted_EE); x: Et (GeV)", 60);
00669  hA[85]=new  L2Histo(85,"diJet2 Et   (accepted_EE); x: Et (GeV)", 60);
00670 
00671  hA[86]=new  L2Histo(86,"diJet1  eta (accepted_EE); x: i Eta [-1,+2]", 15);
00672  hA[87]=new  L2Histo(87,"diJet2  eta (accepted_EE); x: i Eta [-1,+2]", 15);
00673  hA[88]=new  L2Histo(88,"diJet1 phi (accepted_EE); x: iPhi ~sector", 30);
00674  hA[89]=new  L2Histo(89,"diJet2 phi (accepted_EE); x: iPhi ~sector", 30);
00675  hA[90]=new  L2Histo(90,"diJet delZeta  (accepted_EE); x: delta zeta  (rad*10)", MxPhiRad10);
00676  hA[91]=new  L2Histo(91,"diJet delZeta vs. eta1 (accepted_EE); x: iEta1 [-1,+2] ; y: delta zeta  (rad*10)",15, MxPhiRad10);
00677  hA[92]=new  L2Histo(92,"diJet eta2 vs. eta1  (accepted_EE); x: iEta1 [-1,+2] ;x: iEta2 [-1,+2] ",15,15);
00678  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);
00679  hA[94]=new  L2Histo(94,"total Et diJet (accepted_EE); x: Et (GeV)", 60);
00680 
00681 // Di-Jet_EB accepted ... all these histos #100-#114 added  ... wwj 2/10
00682  hA[100]=new   L2Histo(100,"Et of Jet1 vs. Jet2  (accepted_EB); x: Jet1/GeV ; Jet2/GeV",30,30);
00683  hA[101]=new   L2Histo(101,"diJet1 eta-phi   (accepted_EB); x: iEta [-1,+2] ; y: iPhi ~sector ",15,30);
00684  hA[102]=new   L2Histo(102,"diJet2 eta-phi   (accepted_EB); x: iEta [-1,+2]  ; y: iPhi ~sector",15,30);
00685 
00686  hA[103]=new  L2Histo(103,"diJet phi1-phi2   (accepted_EB); x: iPhi1 ~sector ; y: iPhi2 ~sector ",30,30);
00687 
00688  hA[104]=new  L2Histo(104,"diJet1 Et  (accepted_EB); x: Et (GeV)", 60);
00689  hA[105]=new  L2Histo(105,"diJet2 Et   (accepted_EB); x: Et (GeV)", 60);
00690 
00691  hA[106]=new  L2Histo(106,"diJet1  eta (accepted_EB); x: i Eta [-1,+2]", 15);
00692  hA[107]=new  L2Histo(107,"diJet2  eta (accepted_EB); x: i Eta [-1,+2]", 15);
00693  hA[108]=new  L2Histo(108,"diJet1 phi (accepted_EB); x: iPhi ~sector", 30);
00694  hA[109]=new  L2Histo(109,"diJet2 phi (accepted_EB); x: iPhi ~sector", 30);
00695  hA[110]=new  L2Histo(110,"diJet delZeta  (accepted_EB); x: delta zeta  (rad*10)", MxPhiRad10);
00696  hA[111]=new  L2Histo(111,"diJet delZeta vs. eta1 (accepted_EB); x: iEta1 [-1,+2] ; y: delta zeta  (rad*10)",15, MxPhiRad10);
00697  hA[112]=new  L2Histo(112,"diJet eta2 vs. eta1  (accepted_EB); x: iEta1 [-1,+2] ;x: iEta2 [-1,+2] ",15,15);
00698  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);
00699  hA[114]=new  L2Histo(114,"total Et diJet (accepted_EB); x: Et (GeV)", 60);
00700 
00701 // Di-Jet_BB accepted ... all these histos #120-#134 added ... wwj 2/10 
00702  hA[120]=new   L2Histo(120,"Et of Jet1 vs. Jet2  (accepted_BB); x: Jet1/GeV ; Jet2/GeV",30,30);
00703  hA[121]=new   L2Histo(121,"diJet1 eta-phi   (accepted_BB); x: iEta [-1,+2] ; y: iPhi ~sector ",15,30);
00704  hA[122]=new   L2Histo(122,"diJet2 eta-phi   (accepted_BB); x: iEta [-1,+2]  ; y: iPhi ~sector",15,30);
00705 
00706  hA[123]=new  L2Histo(123,"diJet phi1-phi2   (accepted_BB); x: iPhi1 ~sector ; y: iPhi2 ~sector ",30,30);
00707 
00708  hA[124]=new  L2Histo(124,"diJet1 Et  (accepted_BB); x: Et (GeV)", 60);
00709  hA[125]=new  L2Histo(125,"diJet2 Et   (accepted_BB); x: Et (GeV)", 60);
00710 
00711  hA[126]=new  L2Histo(126,"diJet1  eta (accepted_BB); x: i Eta [-1,+2]", 15);
00712  hA[127]=new  L2Histo(127,"diJet2  eta (accepted_BB); x: i Eta [-1,+2]", 15);
00713  hA[128]=new  L2Histo(128,"diJet1 phi (accepted_BB); x: iPhi ~sector", 30);
00714  hA[129]=new  L2Histo(129,"diJet2 phi (accepted_BB); x: iPhi ~sector", 30);
00715  hA[130]=new  L2Histo(130,"diJet delZeta  (accepted_BB); x: delta zeta  (rad*10)", MxPhiRad10);
00716  hA[131]=new  L2Histo(131,"diJet delZeta vs. eta1 (accepted_BB); x: iEta1 [-1,+2] ; y: delta zeta  (rad*10)",15, MxPhiRad10);
00717  hA[132]=new  L2Histo(132,"diJet eta2 vs. eta1  (accepted_BB); x: iEta1 [-1,+2] ;x: iEta2 [-1,+2] ",15,15);
00718  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);
00719  hA[134]=new  L2Histo(134,"total Et diJet (accepted_BB); x: Et (GeV)", 100);
00720 
00721 }
00722 
00723 //=======================================
00724 //=======================================
00725 void 
00726 L2jetAlgo2012::clearEvent(){
00727   mAccept=false;
00728   memset(eve_patchEne,0,sizeof(eve_patchEne));
00729   memset(eve_phiEne,0,sizeof(eve_phiEne));
00730   eve_Jet[0]->clear();
00731   eve_Jet[1]->clear();
00732 }
00733 
00734 
00735 
00736 //=======================================
00737 //=======================================
00738 int 
00739 L2jetAlgo2012::projectAdc(const HitTower1 *hit,const int hitSize,
00740              ushort *phiBin, ushort *patchBin,
00741              L2Histo *hHot       ){
00742   int tmpNused=0; /* counts mapped & used ADC channels */
00743   short rdo;
00744   float low_noise_et;
00745 
00746   for(int i=0;i< hitSize;i++,hit++) {
00747     rdo=hit->rdo;
00748     low_noise_et=hit->low_noise_et;
00749     eve_patchEne[patchBin[rdo]]+=low_noise_et;
00750     eve_phiEne[phiBin[rdo]]+=low_noise_et;
00751     tmpNused++;
00752 
00753     if(low_noise_et > par_hotTwEtThr) hHot->fill(rdo);//do we want rdo or i?
00754 
00755   }
00756   return tmpNused;
00757 }
00758 
00759 
00760 //========================================
00761 //========================================
00762 float L2jetAlgo2012::true2Dscan(){
00763   //implament a true 2-D scan to test the performance of the current algo
00764   int iphiBinEdge = 0;
00765   int ietaBinEdge = 0;
00766   float maxPatchEt = 0.0;
00767   float sumTot = 0.0;
00768 
00769   float *totalCaloEt=eve_phiEne;
00770   float  eneA[cl2jetMaxEtaBins];
00771   float  secondPatchArray[cl2jetMaxPhiBins][cl2jetMaxEtaBins-cl2jet_par_mxEtaBin+1];
00772   memset(secondPatchArray,0,sizeof(secondPatchArray));
00773 
00774   int i;//loops over all phi bin. phi edge
00775   int j;//collapses phi bins
00776   int j1;//phi wrapped phi bin
00777   int j2;//loops over eta bins and sums energy
00778   int k;//loops over eta bins. eta edge
00779 
00780   for(i=0;i<cl2jetMaxPhiBins;i++)
00781     {
00782       sumTot+=totalCaloEt[i];
00783       memset(eneA,0,sizeof(eneA));
00784       for(j=0;j<cl2jet_par_mxPhiBin;j++)
00785         {
00786           j1 = (i+j)%cl2jetMaxPhiBins;//must wrap up phi
00787           float *patchEneA=eve_patchEne+(j1*cl2jetMaxEtaBins);
00788           for(j2=0;j2<cl2jetMaxEtaBins;j2++,patchEneA++)
00789             {
00790               eneA[j2]+=*patchEneA;
00791             }
00792         }
00793       float *eneAp=eneA;
00794       for(k=0;k<cl2jetMaxEtaBins-cl2jet_par_mxEtaBin+1;k++,eneAp++) 
00795         {
00796           float sum=eneAp[0]+eneAp[1]+eneAp[2]+eneAp[3]+eneAp[4];
00797           secondPatchArray[i][k]=sum;
00798           if(maxPatchEt>sum) continue;
00799           //float *eneSave=eneA;
00800           maxPatchEt=sum;
00801           iphiBinEdge=i;
00802           ietaBinEdge=k;
00803         }
00804     }
00805 
00806   eve_Jet[0]->iphiBin=iphiBinEdge;
00807   eve_Jet[0]->ietaBin=ietaBinEdge;
00808   eve_Jet[0]->iene=maxPatchEt;
00809 
00810   //now find second highest patch
00811   int iphiBinEdge2=-1;
00812   int ietaBinEdge2=-1;
00813   float maxPatchEt2=0.0;
00814   char doWrap=0;
00815   int a1=iphiBinEdge-par_minPhiBinDiff;
00816   int a2=iphiBinEdge+par_minPhiBinDiff;
00817   if (a1<0) { a1+=cl2jetMaxPhiBins; doWrap+=1; }
00818   if (a2>=cl2jetMaxPhiBins) { a2-=cl2jetMaxPhiBins; doWrap+=2; }
00819     
00820   int b;
00821   int b1;
00822   if (!doWrap)//masking range does not wrap up
00823     {
00824       for(b=0;b<cl2jetMaxPhiBins;b++)
00825         {
00826           if(b>=a1 && b<=a2) continue;
00827           for(b1=0;b1<cl2jetMaxEtaBins-cl2jet_par_mxEtaBin+1;b1++)
00828             {
00829               if(maxPatchEt2>secondPatchArray[b][b1]) continue;
00830               maxPatchEt2=secondPatchArray[b][b1];
00831               iphiBinEdge2=b;
00832               ietaBinEdge2=b1;
00833             }
00834         }
00835     } 
00836   else
00837     {//masking range wraps up
00838       for(b=a2;b<a1;b++)
00839         {
00840           for(b1=0;b1<cl2jetMaxEtaBins-cl2jet_par_mxEtaBin+1;b1++)
00841             {
00842               if(maxPatchEt2>secondPatchArray[b][b1]) continue;
00843               maxPatchEt2=secondPatchArray[b][b1];
00844               iphiBinEdge2=b;
00845               ietaBinEdge2=b1;
00846             }
00847         }
00848     }
00849 
00850   eve_Jet[1]->iphiBin=iphiBinEdge2;
00851   eve_Jet[1]->ietaBin=ietaBinEdge2;
00852   eve_Jet[1]->iene=maxPatchEt2;
00853 
00854   return sumTot;
00855 }
00856 
00857 //========================================
00858 //========================================
00859 void L2jetAlgo2012::weightedEtaPhi(int iK){
00860   L2Jet *J=eve_Jet[iK];
00861 
00862   // empty calo protection
00863   if(J->iene<=0.1) 
00864     {
00865       J->fphiBin=J->iphiBin+.333;
00866       J->fetaBin=J->ietaBin+.333;
00867       J->rmsPhiBin=0.0;
00868       J->rmsEtaBin=0.0;
00869       return;
00870     }
00871 
00872   int iphi0=J->iphiBin;
00873   int ieta0=J->ietaBin;
00874   float iene0=J->iene;
00875 
00876   //variables for weighted phi
00877   float sumY=0.0;
00878   float sum1;//wrk variable
00879 
00880   //these variables for phi RMS:
00881   double sumYY = 0.0; 
00882   int nEnePhi=0; //  counts bins with energy for RMS computation
00883 
00884   //variable for weighted eta
00885   float sumX=0.0;
00886 
00887   //these variables for eta RMS:
00888   double sumXX = 0.0;
00889   int nEneEta = 0;
00890 
00891   float  etaEneA[cl2jet_par_mxEtaBin];
00892   memset(etaEneA,0,sizeof(etaEneA));
00893 
00894   int ix,iy;
00895   // pick 5x5 eta-phi patch and calculate weighted phi==Y
00896   for(iy=iphi0;iy<iphi0+cl2jet_par_mxPhiBin;iy++)
00897     {
00898       int jy=iy % cl2jetMaxPhiBins;
00899       float  *patchEneA=eve_patchEne+(jy*cl2jetMaxEtaBins);// phi bins are consecutive
00900     
00901       sum1=0.0;
00902       for(ix=ieta0;ix<ieta0+cl2jet_par_mxEtaBin;ix++)
00903         {
00904           sum1+=patchEneA[ix];
00905           etaEneA[ix-ieta0]+=patchEneA[ix];
00906         }
00907       if(sum1>0) nEnePhi++;
00908       sumY+=sum1*iy;
00909       sumYY+=sum1*iy*iy;
00910     }
00911 
00912   // pick 5x5 eta-phi patch and calculate weighted eta==X
00913   int jx;
00914   for(jx=ieta0;jx<ieta0+cl2jet_par_mxEtaBin;jx++)
00915     {
00916       if(etaEneA[jx-ieta0]>0) nEneEta++;
00917       sumX+=jx*etaEneA[jx-ieta0];
00918       sumXX+=(jx*jx*etaEneA[jx-ieta0]);
00919     }
00920 
00921   //calculate phi mean and rms
00922   float phiMean= 1.*sumY/iene0;//was phiSum
00923   float fphiBinMax=0.5 +phiMean; 
00924   float rmsPhi=0.0; // for single bin RMS is 0
00925   if(nEnePhi>1)
00926     {   //here I calculate RMS
00927       rmsPhi=sqrt(nEnePhi/(nEnePhi-1.) *(sumYY/iene0 - phiMean*phiMean));
00928     }
00929   if(fphiBinMax>cl2jetMaxPhiBins) fphiBinMax-=cl2jetMaxPhiBins;
00930 
00931   //calculate eta mean and rms
00932   float etaMean = 1.*sumX/iene0;
00933   float fetaBin=0.5 + etaMean; 
00934   float rmsEta=0.0; // for single bin RMS is 0
00935   if(nEneEta>1)
00936     {   //here I calculate RMS
00937       rmsEta=sqrt(nEneEta/(nEneEta-1.) *(sumXX/iene0 - etaMean*etaMean));
00938     }
00939 
00940   J->fphiBin=fphiBinMax;
00941   J->fetaBin=fetaBin;
00942   J->rmsPhiBin=rmsPhi;
00943   J->rmsEtaBin=rmsEta;
00944 
00945 }
00946 
00947 //========================================
00948 //========================================
00949 void 
00950 L2jetAlgo2012:: dumpPatchEneA(){
00951   // dump L2 array with energy 
00952   int ix,iy;
00953   for(iy=0;iy<cl2jetMaxPhiBins;iy++) {
00954     float  *patchEneA=eve_patchEne+(iy*cl2jetMaxEtaBins);// phi bins are consecutive
00955     
00956     for(ix=0;ix<cl2jetMaxEtaBins;ix++,patchEneA++){
00957       printf(" %6f",*patchEneA);
00958     }
00959     printf(" iPhi=%d\n",iy);
00960   }
00961   
00962 }
00963 
00964 
00965 //=======================================
00966 //=======================================
00967 void 
00968 L2jetAlgo2012::computeE(int token){
00969 }
00970 //=======================================
00971 //=======================================
00972 void 
00973 L2jetAlgo2012::finishRunHisto(){
00974   // auxialiary operations on histograms at the end of the run
00975 
00976   const int *data20=hA[20]->getData();
00977   const int *data30=hA[30]->getData();
00978 
00979   int bHotSum=1,bHotId=-1;
00980   int eHotSum=1;
00981 
00982   const L2EmcDb2012::EmcCDbItem *xE=mDb->getByIndex(402), *xB=mDb->getByIndex(402);
00983   int i;
00984   for(i=0; i<EmcDbIndexMax; i++) {
00985     const L2EmcDb2012::EmcCDbItem *x=mDb->getByIndex(i);
00986     if(mDb->isEmpty(x)) continue; 
00987     if (mDb->isBTOW(x) ) {
00988       int softId=atoi(x->tube+2);
00989       int ieta= (x->eta-1);
00990       int iphi= (x->sec-1)*10 + x->sub-'a' ;
00991       hA[21]->fillW(softId,data20[x->rdo]);
00992       hA[22]->fillW(ieta, iphi,data20[x->rdo]);
00993       if(bHotSum<data20[x->rdo]) {
00994         bHotSum=data20[x->rdo];
00995         bHotId=softId;
00996         xB=x;
00997       }
00998    }// end of BTOW
00999     else  if (mDb->isETOW(x) ) {
01000       int ihard=x->chan+(x->crate-1)*128;
01001       int ieta= 12-x->eta;
01002       int iphi= (x->sec-1)*5 + x->sub-'A' ;
01003       hA[31]->fillW(ihard,data30[x->rdo]);
01004       hA[32]->fillW(ieta, iphi,data30[x->rdo]);
01005       if(eHotSum<data30[x->rdo]) {
01006         eHotSum=data30[x->rdo];
01007         xE=x;
01008       }
01009 
01010     }// end of BTOW
01011   }
01012   if (mLogFile){
01013     fprintf(mLogFile,"L2jet::finishRun()\n");
01014     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); 
01015    fprintf(mLogFile,"#ETOW_hot tower _candidate_ (eHotSum=%d) :, name %s , crate %d , chan %d\n",eHotSum,xE->name,xE->crate,xE->chan); 
01016   }
01017 
01018 }
01019 
01020 
01021 /**********************************************************************
01022   $Log: L2jetAlgo2012.cxx,v $
01023   Revision 1.4  2012/03/21 18:18:03  jml
01024   got rid of printfs from 2012 files
01025 
01026   Revision 1.3  2011/10/19 16:12:11  jml
01027   more 2012 stuff
01028 
01029   Revision 1.2  2011/10/19 15:39:44  jml
01030   2012
01031 
01032   Revision 1.1  2011/10/18 15:11:43  jml
01033   adding 2012 algorithms
01034 
01035   Revision 1.1  2010/04/17 16:42:14  pibero
01036   *** empty log message ***
01037 
01038   Revision 1.1  2007/11/19 22:18:27  balewski
01039   most L2algos provide triggerID's
01040 
01041   Revision 1.8  2007/11/14 03:58:14  balewski
01042   cleanup of common timing measurement
01043 
01044   Revision 1.7  2007/11/13 23:06:07  balewski
01045   toward more unified L2-algos
01046 
01047   Revision 1.6  2007/11/13 00:12:36  balewski
01048   added offline triggerID, take1
01049 
01050   Revision 1.5  2007/11/08 04:02:31  balewski
01051   run on l2ana as well
01052 
01053   Revision 1.4  2007/11/02 17:43:08  balewski
01054   cleanup & it started to work w/ L2upsilon
01055 
01056   Revision 1.3  2007/11/02 03:03:47  balewski
01057   modified L2VirtualAlgo
01058 
01059   Revision 1.2  2007/10/25 02:07:02  balewski
01060   added L2upsilon & binary event dump
01061 
01062   Revision 1.1  2007/10/11 00:33:19  balewski
01063   L2algo added
01064 
01065   Revision 1.5  2006/03/28 19:46:49  balewski
01066   ver16b, in l2new
01067 
01068   Revision 1.4  2006/03/11 17:08:33  balewski
01069   now CVS comments should work
01070 
01071 */
01072 
01073