StRoot: StTriggerUtilities/L2Emulator/L2wAlgo/L2wBemc2009.cxx Source File

00001 #include <stdio.h>
00002 #include <string.h>
00003 #include <stdlib.h>
00004 #include <time.h>
00005 #include <math.h>
00006 
00007 /*********************************************************************
00008  * $Id: L2wBemc2009.cxx,v 1.2 2009/11/19 15:48:49 balewski Exp $
00009  * \author Jan Balewski,MIT , 2009 
00010  *********************************************************************
00011  * Descripion: see .h
00012   *********************************************************************
00013  */
00014 // super smart macro from Pibero, needs only final output 'long long' as argument
00015 #define rdtscll(val) do { \
00016     unsigned int __a,__d; \
00017     __asm__ __volatile__("rdtsc" : "=a" (__a), "=d" (__d)); \
00018     (val) = ((unsigned long long)__a) | (((unsigned long long)__d)<<32); \
00019 } while(0)
00020 
00021 const float stepETH=5;// for QA histo
00022 
00023 
00024 #ifdef  IS_REAL_L2  //in l2-ana  environment
00025   #include "../L2algoUtil/L2EmcDb.h"
00026   #include "../L2algoUtil/L2Histo.h"
00027 #else
00028   #include "../L2algoUtil/L2EmcDb.h"
00029   #include "../L2algoUtil/L2Histo.h"
00030   #include "../L2algoUtil/L2EmcGeom.h"
00031 #endif
00032 
00033 #include "L2wBemc2009.h"
00034 
00035 //=================================================
00036 //=================================================
00037 L2wBemc2009::L2wBemc2009(const char* name, L2EmcDb* db, L2EmcGeom *geoX, char* outDir, int resOff)  :  L2VirtualAlgo2009( name,  db, outDir, true, false, resOff) { 
00038   /* called one per day
00039      all memory allocation must be done here
00040   */
00041 
00042   mGeom=geoX; 
00043   if (!mGeom)
00044     criticalError((char*)"L2wBemc is broken -- can't find geom.");
00045 
00046   setMaxHist(32); // set upper range, I uses only 2^N -it is easier to remember
00047   createHisto();
00048 
00049   //------- self-consistency checks, should never fail
00050    if (sizeof(L2wResult2009)!= L2wResult2009::mySizeChar) 
00051     criticalError((char*)"L2wBemc has failed consistency check. sizeof(L2wResult2009)!= L2wResult2009::mySizeChar");
00052   
00053 }
00054 
00055 /* ========================================
00056   ======================================== */
00057 int 
00058 L2wBemc2009::initRunUser( int runNo, int *rc_ints, float *rc_floats) {
00059 
00060   // unpack params from run control GUI
00061   par_dbg           =  rc_ints[0];
00062   par_RndAcceptPrescale      =  rc_ints[1];
00063   par_seedEtThres   =  rc_floats[0];
00064   par_clustEtThres  =  rc_floats[1];
00065 
00066   // verify consistency of input params
00067   int kBad=0;
00068   kBad+=0x00002 * (par_clustEtThres<par_seedEtThres);
00069 
00070   // put notes about configuration into the log file
00071   if (mLogFile) { 
00072     fprintf(mLogFile,"L2%s algorithm initRun(R=%d), compiled: %s , %s\n params:\n",
00073             getName(),mRunNumber,__DATE__,__TIME__);
00074     fprintf(mLogFile," - use  seedThres=%.2f (GeV), debug=%d, prescale=%d (0=off,1=100prc, 2=50proc, etc)\n",
00075             par_seedEtThres,par_dbg, par_RndAcceptPrescale);
00076     fprintf(mLogFile," - accept event cluster Thres=%.2f (GeV)\n",
00077             par_clustEtThres);
00078     fprintf(mLogFile,"initRun() params checked for consistency, Error flag=0x%04x\n",
00079             kBad);
00080   }
00081 
00082   // clear content of all histograms & token-dependent memory
00083   int i;
00084   for (i=0; i<mxHA;i++) if(hA[i])hA[i]->reset();
00085   memset(mBtow,0,sizeof(mBtow));
00086 
00087   if(kBad>0) return -1*kBad;    
00088   else if(kBad<0) return kBad;
00089 
00090  
00091   // update titles of histos
00092   char txt[1000];
00093 
00094   sprintf(txt,"W-BTOW-accepted: acc seed tower ET>%.2f GeV; BTOW softID",par_seedEtThres);
00095   hA[9]->setTitle(txt);
00096   
00097   sprintf(txt,"W Btow accepted, seed tower Et>%.2f GeV; eta bin [-1,+1]; y: phi bin ~ TCP sector",par_seedEtThres);
00098   hA[10]->setTitle(txt);
00099   
00100 
00101   for ( int index=0; index<EmcDbIndexMax; index++ )
00102      {
00103        const L2EmcDb::EmcCDbItem *x = mDb->getByIndex(index);
00104        if ( x==0 ) continue;
00105        if ( !mDb->isBTOW(x) ) continue; 
00106        int sec = x->sec - 1;
00107        int sub = 8192; 
00108        sub = x->sub - 'a';
00109        int eta = x->eta - 1;
00110        int phi = BtowGeom::mxSubs *sec + sub;
00111        int tow = BtowGeom::mxEtaBin *phi + eta; // phi- changes faster
00112        int rdo = x->rdo;
00113        if (tow<0 || tow>mxBtow || rdo<0 || rdo>mxBtow) return -101;
00114        
00115        mTower2rdo[ tow ] = rdo;    // returns rdo channel for given tower
00116        mRdo2tower[ rdo ] = tow;    // returns tower for given rdo channel
00117     }
00118   return 0; //OK
00119 
00120 }               
00121 
00122   
00123 
00124 /* ========================================
00125   ======================================== */
00126 float
00127 L2wBemc2009::sumET(int phi, int eta) {
00128   int tow = BtowGeom::mxEtaBin *((phi+BtowGeom::mxPhiBin)%BtowGeom::mxPhiBin) + ((eta+BtowGeom::mxEtaBin)%BtowGeom::mxEtaBin); // phi- changes faster
00129 
00130   const int maxTowers = BtowGeom::mxEtaBin * BtowGeom::mxPhiBin;
00131   int towPlusOne;
00132   float sum;
00133   sum=0;
00134   sum=wrkBtow_et[tow];
00135   towPlusOne = tow+1;
00136   towPlusOne%= maxTowers;
00137   sum+=wrkBtow_et[towPlusOne];
00138   
00139   //if(tow==0 || towPlusOne==0) {
00140     //printf("tow : %d, %d --> %f %f \n",tow, towPlusOne,wrkBtow_et[tow],wrkBtow_et[towPlusOne]);
00141   //}
00142   tow+=BtowGeom::mxEtaBin;
00143   tow%=maxTowers;
00144   
00145   sum+=wrkBtow_et[tow];
00146   
00147   towPlusOne = tow+1;
00148   towPlusOne%= maxTowers;
00149   sum+=wrkBtow_et[towPlusOne];
00150   
00151   //if(tow==0 || towPlusOne==0) {
00152     //printf("tow : %d, %d --> %f %f \n",tow, towPlusOne,wrkBtow_et[tow],wrkBtow_et[towPlusOne]);
00153   //}
00154   // printf("B sumET=%f\n",sum);
00155   return sum;
00156 }
00157   
00158 
00159 /* ========================================
00160   ======================================== */
00161 void 
00162 L2wBemc2009::computeUser(int token){
00163 
00164   clearEvent(token);
00165 
00166   // ------ PROJECT INPUT LIST TO 2D ARRAY AND SCAN FOR SEED TOWERS ----
00167   int i;
00168   //  printf("L2-%s-compute: ---BTOW ADC list--- size=%d\n",getName(),*globEve_btow_hitSize);
00169 
00170   // access inpute list
00171   const HitTower1 *hit=mEveStream_btow[token].get_hits();
00172   const int hitSize=mEveStream_btow[token].get_hitSize();
00173   // preapre output list
00174   L2wBemcEvent2009 *btowEve=mBtow+token;
00175 
00176   for(i=0;i< hitSize;i++,hit++) {
00177     int tower=mRdo2tower[hit->rdo];
00178     wrkBtow_et[tower]=hit->et;
00179     if(hit->et<par_seedEtThres)continue;
00180     wrkBtow_tower_seed[wrkBtow_tower_seed_size++]=tower;
00181     //printf("A  seed TWID=%d \n",tower);
00182   }
00183   hA[2]->fill(hitSize);
00184   int seedTow=-1,seedEta=-1,seedPhi=-1;
00185   float clustET=0;
00186   btowEve->isFresh=L2wBemcEvent2009::kDataFresh;
00187 
00188   // printf("B nseed=%d\n",wrkBtow_tower_seed_size);
00189   // ----------- FIND 2x2 CLUSTER AROUND EVERY SEED -----
00190   for(i=0; i<wrkBtow_tower_seed_size;i++) {
00191     seedTow=wrkBtow_tower_seed[i];
00192     seedEta=seedTow%BtowGeom::mxEtaBin;
00193     seedPhi=seedTow/BtowGeom::mxEtaBin;
00194     
00195     //.... find first 2x2 above cluster thresh
00196     if (seedEta < BtowGeom::mxEtaBin) {
00197       clustET = sumET(seedPhi,seedEta);
00198       if(clustET>par_clustEtThres) goto ACCEPT;
00199       clustET = sumET(seedPhi-1,seedEta);
00200       if(clustET>par_clustEtThres) goto ACCEPT;
00201     }
00202     if (seedEta > 0 ) {
00203       clustET = sumET(seedPhi-1,seedEta-1);
00204       if(clustET>par_clustEtThres) goto ACCEPT;
00205       clustET = sumET(seedPhi,seedEta-1);
00206       if(clustET>par_clustEtThres) goto ACCEPT;
00207     }
00208   }
00209   //.... ABORT
00210   btowEve->resultBlob.seedEt=0;
00211   btowEve->resultBlob.clusterEt=0;
00212   btowEve->resultBlob.seedEtaBin=0;
00213   btowEve->resultBlob.seedPhiBin=0;
00214   btowEve->resultBlob.trigger=0;
00215   btowEve->seedET=0;
00216   btowEve->clusterET=0;
00217   btowEve->tkCompute=0;
00218   return;
00219   
00220  ACCEPT:
00221   // printf("B acc clustET=%f\n",clustET);
00222   btowEve->seedET=wrkBtow_et[seedTow];
00223   btowEve->clusterET=clustET;
00224   btowEve->resultBlob.seedEt   =(unsigned char)(wrkBtow_et[seedTow]*256.0/60.0);
00225   btowEve->resultBlob.clusterEt=(unsigned char)(clustET*256.0/60.0);
00226   btowEve->resultBlob.seedEtaBin=seedEta;
00227   btowEve->resultBlob.seedPhiBin=seedPhi;
00228   btowEve->resultBlob.trigger=2;
00229   rdtscll( btowEve->tkCompute);
00230   //  printf("TTT1 %dl\n",ticks);
00231 
00232   return;
00233 }
00234 
00235 
00236 
00237 
00238 /* ========================================
00239   ======================================== */
00240 bool 
00241 L2wBemc2009::decisionUser(int token, int *myL2Result){
00242   // INPUT: token + comput() results stored internally
00243   // OUTPUT: yes/now + pointer to  L2Result
00244 
00245   // get pointers to internal private event storage
00246   L2wBemcEvent2009 *btowEve=mBtow+token;
00247 
00248   // prescaling decison ws already made in Virtual:
00249   // if (par_RndAcceptPrescale>0 && mRandomAccept) btowEve->resultBlob.trigger+=1;
00250   if (mRandomAccept) btowEve->resultBlob.trigger+=1;
00251   
00252   //...... some histos just for QA
00253   if(btowEve->isFresh>L2wBemcEvent2009::kDataFresh) mhN->fill(6); // stale data, should never happen
00254 
00255   btowEve->isFresh++; // mark the data as  stale
00256 
00257   if(btowEve->resultBlob.trigger&2) {
00258     unsigned long long tkDecision;
00259     rdtscll(tkDecision);
00260     int tkDelta=tkDecision-btowEve->tkCompute;                                  
00261     //printf("TTT2 %d\n",tkDelta);
00262     hA[1]->fill(tkDelta/1000);
00263 
00264     hA[3]->fill((int)btowEve->seedET);
00265     hA[4]->fill((int)btowEve->clusterET);
00266     hA[5]->fill((int)(100.*btowEve->seedET/btowEve->clusterET));
00267     int jET=(int)(btowEve->clusterET/stepETH);
00268     if( jET<=1) jET=1;
00269     if( jET>=9) jET=9;
00270     hA[10+jET]->fill(btowEve->resultBlob.seedEtaBin,btowEve->resultBlob.seedPhiBin);
00271     hA[10]->fill(btowEve->resultBlob.seedEtaBin,btowEve->resultBlob.seedPhiBin);
00272   }
00273   // store final content in TRG data
00274   memcpy(myL2Result,&(btowEve->resultBlob),sizeof(L2wResult2009));
00275   if(par_dbg) L2wResult2009_print((L2wResult2009 *)&(btowEve->resultBlob));
00276   //return btowEve->resultBlob.trigger;
00277   return btowEve->resultBlob.trigger&2;
00278 } 
00279 
00280 
00281 /* ========================================
00282   ======================================== */
00283 void
00284 L2wBemc2009::finishRunUser() {  /* called once at the end of the run */
00285   // do whatever you want, log-file & histo-file are still open
00286   
00287   if (mLogFile){ 
00288     fprintf(mLogFile,"finishRunUser-%s start\n",getName());
00289   }
00290   
00291   // test1:
00292   //  hA[9]->fillW(21,1234);
00293   
00294   //......... scan for hot towers ....& repack 2D--> 1D(softId)
00295 
00296   for(int jh=10;jh<19;jh++) 
00297   {
00298     const L2EmcDb::EmcCDbItem *xB=mDb->getByIndex(402);
00299     const int *hiData = hA[jh]->getData();
00300     int hotY=0,totY=0;
00301     for(int i=0; i<EmcDbIndexMax; i++) {
00302       const L2EmcDb::EmcCDbItem *x=mDb->getByIndex(i);
00303       if(mDb->isEmpty(x)) continue;
00304       if (!mDb->isBTOW(x) ) continue;
00305       int ieta= (x->eta-1);
00306       int iphi= (x->sec-1)*10 + x->sub-'a' ;
00307       int softId=atoi(x->tube+2);
00308       int ix=ieta+iphi*BtowGeom::mxEtaBin;
00309       int yield=hiData[ix];
00310       totY+=yield;
00311       if(jh==10) hA[9]->fillW(softId,yield);// do it only for ET-integrated
00312       if(hotY>yield)continue;
00313       hotY=yield;
00314       xB=x;
00315     }
00316 
00317     int ieta= (xB->eta-1);
00318     int iphi= (xB->sec-1)*10 + xB->sub-'a' ;
00319     int softId=atoi(xB->tube+2);
00320     fprintf(mLogFile,"#BTOW hot *candidate*, hist: %d,  yield: %d, totYield: %d, softID: %d, crate: %d, chan: %d, name: %s, ieta: %d, iphi: %d\n", hA[jh]->getId(),hotY,totY,softId,xB->crate,xB->chan,xB->name,ieta,iphi);
00321   }
00322 
00323     
00324 
00325 }
00326 
00327 
00328 //=======================================
00329 //=======================================
00330 void 
00331 L2wBemc2009::createHisto() {
00332   hA[1]=new L2Histo(1,(char*)"W Btow delTime (decision-compute); kTicks",300);
00333   
00334   hA[2]=new L2Histo(2,(char*)"W Btow-compute: # btow towers w/ energy /event; x: # BTOW towers; y: counts", 100); 
00335   hA[3]=new L2Histo(3,(char*)"W Btow-accepted: seeds ; Seed ET (GeV)", 70); 
00336   hA[4]=new L2Histo(4,(char*)"W Btow-accepted: clusters ; Cluster ET ET(GeV)", 70);
00337   hA[5]=new L2Histo(5,(char*)"W Btow-accepted: cluster shape ;100*Seed Et/Cluster Et", 105);
00338   // 6-8 free
00339 
00340   hA[9]=new L2Histo(9,(char*)"W Btow ...9 .....", 5000); // title in initRun
00341   hA[10]=new L2Histo(10,(char*)"W Btow ...10 ...", BtowGeom::mxEtaBin, BtowGeom::mxPhiBin); // title in initRun
00342 
00343 
00344   char tit[100];
00345   for(int j=1;j<=9;j++){
00346     int hid=10+j;
00347     float x1=j*stepETH;
00348     float x2=x1+stepETH;
00349     if(j==1) x1=0;
00350     if(j==9) x2=999;
00351     sprintf(tit,"W Btow-accept: cluster  %.0f <ET< %.0f GeV ; seed eta bin [-1,+1];  seed phi bin ~ TCP sector", x1,x2);
00352 
00353     hA[hid]=new L2Histo(hid,tit, BtowGeom::mxEtaBin, BtowGeom::mxPhiBin);
00354   }
00355 
00356   //20-32 free
00357 
00358   // printf("L2-%s::createHisto() done\n",getName());
00359 }
00360 
00361 //=======================================
00362 //=======================================
00363 void 
00364 L2wBemc2009::clearEvent(int token){
00365   memset(wrkBtow_et,0,sizeof(wrkBtow_et)); 
00366   wrkBtow_tower_seed_size=0; 
00367   memset(&(mBtow[token].resultBlob),0, sizeof(L2wResult2009));
00368 }
00369   
00370 /* ========================================
00371   ======================================== */
00372 void 
00373 L2wBemc2009::print2(){ // full , local ADC array
00374   int i;
00375   printf("pr2-%s: ---BTOW ADC 2D array, only non-zero\n",getName());
00376 
00377   for(i=0;i<mxBtow;i++) {
00378     if(wrkBtow_et[i]<=0) continue;
00379     int rdo=mTower2rdo[i];
00380     float et=wrkBtow_et[i];
00381     printf("  btow: tower=%4d  rdo=%4d   et=%.3f \n",i,rdo,et);
00382   }
00383 
00384 }
00385 
00386 /* ========================================
00387   ======================================== */
00388 void 
00389 L2wBemc2009::print3(){ // seed list
00390   int i;
00391   printf("pr3-%s: ---seed list, size=%d\n",getName(),wrkBtow_tower_seed_size);
00392 
00393   for(i=0;i<wrkBtow_tower_seed_size;i++) {
00394     int tower=wrkBtow_tower_seed[i];
00395     float et=wrkBtow_et[tower];
00396     printf("  btow: i=%4d  tower=%4d   et=%.3f \n",i,tower,et);
00397   }
00398 
00399 }
00400 
00401 #if 0
00402 
00403 /* ========================================
00404   ======================================== */
00405 void 
00406 L2wBemc2009::print4(int token, int hitSize){ // L2-algo input list
00407   int i;
00408   printf("print4 IS NOT Fully FUNCTIONAL **********************\n");
00409   printf("pr1-%s: ---BTOW Sorted ADC list--- size=%d\n",getName(),hitSize);
00410   //const HitTower *hit=globEve_btow_hit;
00411   for(i=0;i< hitSize;i++) {
00412     int adc=0;//(mEveStream_btow[token].get_hits()[wrkBtow_tower_index[i]]).adc;
00413     int rdo=0;//(mEveStream_btow[token].get_hits()[wrkBtow_tower_index[i]]).rdo;
00414     float et=wrkBtow_et[wrkBtow_tower_index[i]];
00415     float ene=0;//(mEveStream_btow[token].get_hits()[wrkBtow_tower_index[i]]).ene;
00416     printf("  tower=%2d ",wrkBtow_tower_index[i]);
00417     printf("  btow: i=%2d rdo=%4d  adc=%d  et=%.3f  ene=%.3f\n",i,rdo,adc,et,ene);
00418   }
00419 }
00420 
00421 #endif
00422 /**********************************************************************
00423   $Log: L2wBemc2009.cxx,v $
00424   Revision 1.2  2009/11/19 15:48:49  balewski
00425   add (char*) to many strings to make SL5 happ, few other adjustments
00426 
00427   Revision 1.1  2009/03/28 19:43:53  balewski
00428   2009 code
00429 
00430 
00431  
00432 */
00433 
00434