StRoot: StTriggerUtilities/L2Emulator/L2wAlgo/L2wBemc2012.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: L2wBemc2012.cxx,v 1.4 2012/03/21 18:18:04 jml 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/L2EmcDb2012.h"
00026   #include "../L2algoUtil/L2Histo.h"
00027 #else
00028   #include "../L2algoUtil/L2EmcDb2012.h"
00029   #include "../L2algoUtil/L2Histo.h"
00030   #include "../L2algoUtil/L2EmcGeom2012.h"
00031 #endif
00032 
00033 #include "L2wBemc2012.h"
00034 
00035 //=================================================
00036 //=================================================
00037 L2wBemc2012::L2wBemc2012(const char* name, const char *uid, L2EmcDb2012* db, L2EmcGeom2012 *geoX, char* outDir, int resOff)  :  L2VirtualAlgo2012( name,uid,  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(L2wResult2012)!= L2wResult2012::mySizeChar) 
00051     criticalError((char*)"L2wBemc has failed consistency check. sizeof(L2wResult2012)!= L2wResult2012::mySizeChar");
00052   
00053 }
00054 
00055 /* ========================================
00056   ======================================== */
00057 int 
00058 L2wBemc2012::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 L2EmcDb2012::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 L2wBemc2012::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   tow+=BtowGeom::mxEtaBin;
00140   tow%=maxTowers;
00141   
00142   sum+=wrkBtow_et[tow];
00143   
00144   towPlusOne = tow+1;
00145   towPlusOne%= maxTowers;
00146   sum+=wrkBtow_et[towPlusOne];
00147   
00148   return sum;
00149 }
00150   
00151 
00152 /* ========================================
00153   ======================================== */
00154 void 
00155 L2wBemc2012::computeUser(int token){
00156 
00157   clearEvent(token);
00158 
00159   // ------ PROJECT INPUT LIST TO 2D ARRAY AND SCAN FOR SEED TOWERS ----
00160   int i;
00161 
00162   // access inpute list
00163   const HitTower1 *hit=mEveStream_btow[token].get_hits();
00164   const int hitSize=mEveStream_btow[token].get_hitSize();
00165   // preapre output list
00166   L2wBemcEvent2012 *btowEve=mBtow+token;
00167 
00168   for(i=0;i< hitSize;i++,hit++) {
00169     int tower=mRdo2tower[hit->rdo];
00170     wrkBtow_et[tower]=hit->et;
00171     if(hit->et<par_seedEtThres)continue;
00172     wrkBtow_tower_seed[wrkBtow_tower_seed_size++]=tower;
00173   }
00174   hA[2]->fill(hitSize);
00175   int seedTow=-1,seedEta=-1,seedPhi=-1;
00176   float clustET=0;
00177   btowEve->isFresh=L2wBemcEvent2012::kDataFresh;
00178 
00179   // ----------- FIND 2x2 CLUSTER AROUND EVERY SEED -----
00180   for(i=0; i<wrkBtow_tower_seed_size;i++) {
00181     seedTow=wrkBtow_tower_seed[i];
00182     seedEta=seedTow%BtowGeom::mxEtaBin;
00183     seedPhi=seedTow/BtowGeom::mxEtaBin;
00184     
00185     //.... find first 2x2 above cluster thresh
00186     if (seedEta < BtowGeom::mxEtaBin) {
00187       clustET = sumET(seedPhi,seedEta);
00188       if(clustET>par_clustEtThres) goto ACCEPT;
00189       clustET = sumET(seedPhi-1,seedEta);
00190       if(clustET>par_clustEtThres) goto ACCEPT;
00191     }
00192     if (seedEta > 0 ) {
00193       clustET = sumET(seedPhi-1,seedEta-1);
00194       if(clustET>par_clustEtThres) goto ACCEPT;
00195       clustET = sumET(seedPhi,seedEta-1);
00196       if(clustET>par_clustEtThres) goto ACCEPT;
00197     }
00198   }
00199   //.... ABORT
00200   btowEve->resultBlob.seedEt=0;
00201   btowEve->resultBlob.clusterEt=0;
00202   btowEve->resultBlob.seedEtaBin=0;
00203   btowEve->resultBlob.seedPhiBin=0;
00204   btowEve->resultBlob.trigger=0;
00205   btowEve->seedET=0;
00206   btowEve->clusterET=0;
00207   btowEve->tkCompute=0;
00208   return;
00209   
00210  ACCEPT:
00211   btowEve->seedET=wrkBtow_et[seedTow];
00212   btowEve->clusterET=clustET;
00213   btowEve->resultBlob.seedEt   =(unsigned char)(wrkBtow_et[seedTow]*256.0/60.0);
00214   btowEve->resultBlob.clusterEt=(unsigned char)(clustET*256.0/60.0);
00215   btowEve->resultBlob.seedEtaBin=seedEta;
00216   btowEve->resultBlob.seedPhiBin=seedPhi;
00217   btowEve->resultBlob.trigger=2;
00218   rdtscll( btowEve->tkCompute);
00219 
00220   return;
00221 }
00222 
00223 
00224 
00225 
00226 /* ========================================
00227   ======================================== */
00228 bool 
00229 L2wBemc2012::decisionUser(int token, int *myL2Result){
00230   // INPUT: token + comput() results stored internally
00231   // OUTPUT: yes/now + pointer to  L2Result
00232 
00233   // get pointers to internal private event storage
00234   L2wBemcEvent2012 *btowEve=mBtow+token;
00235 
00236   // prescaling decison ws already made in Virtual:
00237   // if (par_RndAcceptPrescale>0 && mRandomAccept) btowEve->resultBlob.trigger+=1;
00238   if (mRandomAccept) btowEve->resultBlob.trigger+=1;
00239   
00240   //...... some histos just for QA
00241   if(btowEve->isFresh>L2wBemcEvent2012::kDataFresh) mhN->fill(6); // stale data, should never happen
00242 
00243   btowEve->isFresh++; // mark the data as  stale
00244 
00245   if(btowEve->resultBlob.trigger&2) {
00246     unsigned long long tkDecision;
00247     rdtscll(tkDecision);
00248     int tkDelta=tkDecision-btowEve->tkCompute;                                  
00249     hA[1]->fill(tkDelta/1000);
00250 
00251     hA[3]->fill((int)btowEve->seedET);
00252     hA[4]->fill((int)btowEve->clusterET);
00253     hA[5]->fill((int)(100.*btowEve->seedET/btowEve->clusterET));
00254     int jET=(int)(btowEve->clusterET/stepETH);
00255     if( jET<=1) jET=1;
00256     if( jET>=9) jET=9;
00257     hA[10+jET]->fill(btowEve->resultBlob.seedEtaBin,btowEve->resultBlob.seedPhiBin);
00258     hA[10]->fill(btowEve->resultBlob.seedEtaBin,btowEve->resultBlob.seedPhiBin);
00259   }
00260   // store final content in TRG data
00261   memcpy(myL2Result,&(btowEve->resultBlob),sizeof(L2wResult2012));
00262   if(par_dbg) L2wResult2012_print((L2wResult2012 *)&(btowEve->resultBlob));
00263   //return btowEve->resultBlob.trigger;
00264   return btowEve->resultBlob.trigger&2;
00265 } 
00266 
00267 
00268 /* ========================================
00269   ======================================== */
00270 void
00271 L2wBemc2012::finishRunUser() {  /* called once at the end of the run */
00272   // do whatever you want, log-file & histo-file are still open
00273   
00274   if (mLogFile){ 
00275     fprintf(mLogFile,"finishRunUser-%s start\n",getName());
00276   }
00277   
00278   // test1:
00279   //  hA[9]->fillW(21,1234);
00280   
00281   //......... scan for hot towers ....& repack 2D--> 1D(softId)
00282 
00283   for(int jh=10;jh<19;jh++) 
00284   {
00285     const L2EmcDb2012::EmcCDbItem *xB=mDb->getByIndex(402);
00286     const int *hiData = hA[jh]->getData();
00287     int hotY=0,totY=0;
00288     for(int i=0; i<EmcDbIndexMax; i++) {
00289       const L2EmcDb2012::EmcCDbItem *x=mDb->getByIndex(i);
00290       if(mDb->isEmpty(x)) continue;
00291       if (!mDb->isBTOW(x) ) continue;
00292       int ieta= (x->eta-1);
00293       int iphi= (x->sec-1)*10 + x->sub-'a' ;
00294       int softId=atoi(x->tube+2);
00295       int ix=ieta+iphi*BtowGeom::mxEtaBin;
00296       int yield=hiData[ix];
00297       totY+=yield;
00298       if(jh==10) hA[9]->fillW(softId,yield);// do it only for ET-integrated
00299       if(hotY>yield)continue;
00300       hotY=yield;
00301       xB=x;
00302     }
00303 
00304     int ieta= (xB->eta-1);
00305     int iphi= (xB->sec-1)*10 + xB->sub-'a' ;
00306     int softId=atoi(xB->tube+2);
00307     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);
00308   }
00309 
00310     
00311 
00312 }
00313 
00314 
00315 //=======================================
00316 //=======================================
00317 void 
00318 L2wBemc2012::createHisto() {
00319   hA[1]=new L2Histo(1,(char*)"W Btow delTime (decision-compute); kTicks",300);
00320   
00321   hA[2]=new L2Histo(2,(char*)"W Btow-compute: # btow towers w/ energy /event; x: # BTOW towers; y: counts", 100); 
00322   hA[3]=new L2Histo(3,(char*)"W Btow-accepted: seeds ; Seed ET (GeV)", 70); 
00323   hA[4]=new L2Histo(4,(char*)"W Btow-accepted: clusters ; Cluster ET ET(GeV)", 70);
00324   hA[5]=new L2Histo(5,(char*)"W Btow-accepted: cluster shape ;100*Seed Et/Cluster Et", 105);
00325   // 6-8 free
00326 
00327   hA[9]=new L2Histo(9,(char*)"W Btow ...9 .....", 5000); // title in initRun
00328   hA[10]=new L2Histo(10,(char*)"W Btow ...10 ...", BtowGeom::mxEtaBin, BtowGeom::mxPhiBin); // title in initRun
00329 
00330 
00331   char tit[100];
00332   for(int j=1;j<=9;j++){
00333     int hid=10+j;
00334     float x1=j*stepETH;
00335     float x2=x1+stepETH;
00336     if(j==1) x1=0;
00337     if(j==9) x2=999;
00338     sprintf(tit,"W Btow-accept: cluster  %.0f <ET< %.0f GeV ; seed eta bin [-1,+1];  seed phi bin ~ TCP sector", x1,x2);
00339 
00340     hA[hid]=new L2Histo(hid,tit, BtowGeom::mxEtaBin, BtowGeom::mxPhiBin);
00341   }
00342 
00343   //20-32 free
00344 
00345 }
00346 
00347 //=======================================
00348 //=======================================
00349 void 
00350 L2wBemc2012::clearEvent(int token){
00351   memset(wrkBtow_et,0,sizeof(wrkBtow_et)); 
00352   wrkBtow_tower_seed_size=0; 
00353   memset(&(mBtow[token].resultBlob),0, sizeof(L2wResult2012));
00354 }
00355   
00356 /* ========================================
00357   ======================================== */
00358 void 
00359 L2wBemc2012::print2(){ // full , local ADC array
00360   int i;
00361   printf("pr2-%s: ---BTOW ADC 2D array, only non-zero\n",getName());
00362 
00363   for(i=0;i<mxBtow;i++) {
00364     if(wrkBtow_et[i]<=0) continue;
00365     int rdo=mTower2rdo[i];
00366     float et=wrkBtow_et[i];
00367     printf("  btow: tower=%4d  rdo=%4d   et=%.3f \n",i,rdo,et);
00368   }
00369 
00370 }
00371 
00372 /* ========================================
00373   ======================================== */
00374 void 
00375 L2wBemc2012::print3(){ // seed list
00376   int i;
00377   printf("pr3-%s: ---seed list, size=%d\n",getName(),wrkBtow_tower_seed_size);
00378 
00379   for(i=0;i<wrkBtow_tower_seed_size;i++) {
00380     int tower=wrkBtow_tower_seed[i];
00381     float et=wrkBtow_et[tower];
00382     printf("  btow: i=%4d  tower=%4d   et=%.3f \n",i,tower,et);
00383   }
00384 
00385 }
00386 
00387 #if 0
00388 
00389 /* ========================================
00390   ======================================== */
00391 void 
00392 L2wBemc2012::print4(int token, int hitSize){ // L2-algo input list
00393   int i;
00394   printf("print4 IS NOT Fully FUNCTIONAL **********************\n");
00395   printf("pr1-%s: ---BTOW Sorted ADC list--- size=%d\n",getName(),hitSize);
00396   //const HitTower *hit=globEve_btow_hit;
00397   for(i=0;i< hitSize;i++) {
00398     int adc=0;//(mEveStream_btow[token].get_hits()[wrkBtow_tower_index[i]]).adc;
00399     int rdo=0;//(mEveStream_btow[token].get_hits()[wrkBtow_tower_index[i]]).rdo;
00400     float et=wrkBtow_et[wrkBtow_tower_index[i]];
00401     float ene=0;//(mEveStream_btow[token].get_hits()[wrkBtow_tower_index[i]]).ene;
00402     printf("  tower=%2d ",wrkBtow_tower_index[i]);
00403     printf("  btow: i=%2d rdo=%4d  adc=%d  et=%.3f  ene=%.3f\n",i,rdo,adc,et,ene);
00404   }
00405 }
00406 
00407 #endif
00408 /**********************************************************************
00409   $Log: L2wBemc2012.cxx,v $
00410   Revision 1.4  2012/03/21 18:18:04  jml
00411   got rid of printfs from 2012 files
00412 
00413   Revision 1.3  2011/10/19 16:12:12  jml
00414   more 2012 stuff
00415 
00416   Revision 1.2  2011/10/19 15:39:45  jml
00417   2012
00418 
00419   Revision 1.1  2011/10/18 15:11:44  jml
00420   adding 2012 algorithms
00421 
00422   Revision 1.2  2009/11/19 15:48:49  balewski
00423   add (char*) to many strings to make SL5 happ, few other adjustments
00424 
00425   Revision 1.1  2009/03/28 19:43:53  balewski
00426   2009 code
00427 
00428 
00429  
00430 */
00431 
00432