StRoot: StTriggerUtilities/L2Emulator/L2gammaAlgo/L2bemcGamma2009.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: L2bemcGamma2009.cxx,v 1.1 2011/03/09 16:29:07 pibero Exp $
00009  * \author Jan Balewski,MIT , 2008 
00010  *********************************************************************
00011  * Descripion: see .h
00012   *********************************************************************
00013  */
00014 
00015 
00016 #ifdef  IS_REAL_L2  //in l2-ana  environment
00017   #include "../L2algoUtil/L2EmcDb.h"
00018   #include "../L2algoUtil/L2Histo.h"
00019 #else    //full path needed for cvs'd code
00020   #include "StTriggerUtilities/L2Emulator/L2algoUtil/L2EmcDb.h"
00021   #include "StTriggerUtilities/L2Emulator/L2algoUtil/L2Histo.h"
00022   #include "StTriggerUtilities/L2Emulator/L2algoUtil/L2EmcGeom.h"
00023 #endif
00024 
00025 #include "L2bemcGamma2009.h"
00026 
00027 
00028 //=================================================
00029 //=================================================
00030 L2bemcGamma2009::L2bemcGamma2009(const char* name, L2EmcDb* db, L2EmcGeom *geoX, char* outDir, int resOff)  :  L2VirtualAlgo2009( name,  db, outDir, true, false, resOff) { 
00031   /* called one per days
00032      all memory allocation must be done here
00033   */
00034 
00035   mGeom=geoX; 
00036   if (!mGeom)
00037     criticalError("L2bemcGamma is broken -- can't find geom.");
00038 
00039   setMaxHist(16); // set upper range, I uses only 2^N -it is easier to remember
00040   createHisto();
00041 
00042   //------- self-consistency checks, should never fail
00043   // printf("ZZ %d %d\n", sizeof(L2gammaResult2008), L2gammaResult2008::mySizeChar);
00044   if (sizeof(L2gammaResult2009)!= L2gammaResult2009::mySizeChar) 
00045     criticalError("L2bemcGamma has failed consistency check. sizeof(L2gammaResult2009)!= L2gammaResult2009::mySizeChar");
00046   
00047 }
00048 
00049 /* ========================================
00050   ======================================== */
00051 int 
00052 L2bemcGamma2009::initRunUser( int runNo, int *rc_ints, float *rc_floats) {
00053 
00054   // unpack params from run control GUI
00055   par_dbg           =  rc_ints[0];
00056   par_RndAcceptPrescale      =  rc_ints[1];
00057   par_seedEtThres   =  rc_floats[0];
00058   par_clusterEtThres  =  rc_floats[1];
00059 
00060   // verify consistency of input params
00061   int kBad=0;
00062   kBad+=0x00001 * (par_seedEtThres<1.0);
00063   kBad+=0x00002 * (par_clusterEtThres<par_seedEtThres);
00064 
00065   // put notes about configuration into the log file
00066   if (mLogFile) { 
00067     fprintf(mLogFile,"L2%s algorithm initRun(R=%d), compiled: %s , %s\n params:\n",
00068             getName(),mRunNumber,__DATE__,__TIME__);
00069     fprintf(mLogFile," - use  seedThres=%.2f (GeV),  debug=%d, prescale=%d (0=off,1=100proc, 2=50proc, etc)\n",
00070             par_seedEtThres,par_dbg,par_RndAcceptPrescale);
00071     fprintf(mLogFile," - accept event cluster Thres=%.2f (GeV)\n",
00072             par_clusterEtThres);
00073     fprintf(mLogFile,"initRun() params checked for consistency, Error flag=0x%04x\n",
00074             kBad);
00075   }
00076 
00077   // clear content of all histograms & token-dependent memory
00078   int i;
00079   for (i=0; i<mxHA;i++) if(hA[i])hA[i]->reset();
00080   memset(mBtow,0,sizeof(mBtow));
00081 
00082   if(kBad>0) return -1*kBad;    
00083   else if(kBad<0) return kBad;
00084 
00085  
00086   // update titles of histos
00087   char txt[1000];
00088 
00089   sprintf(txt,"BTOW-decision: acc seed Tw ET>%.2f GeV; x:tower internal ID; y: counts",par_clusterEtThres);
00090   hA[7]->setTitle(txt);
00091   sprintf(txt,"#BTOW decision acc seed Tw, Et>%.2f GeV; eta bin [-1,+1]; y: phi bin ~ TCP sector",par_clusterEtThres);
00092   hA[14]->setTitle(txt);
00093   
00094 
00095   for ( int index=0; index<EmcDbIndexMax; index++ )
00096      {
00097        const L2EmcDb::EmcCDbItem *x = mDb->getByIndex(index);
00098        if ( x==0 ) continue;
00099        if ( !mDb->isBTOW(x) ) continue; 
00100        int sec = x->sec - 1;
00101        int sub = 8192; 
00102        sub = x->sub - 'a';
00103        int eta = x->eta - 1;
00104        int phi = BtowGeom::mxSubs *sec + sub;
00105        int tow = BtowGeom::mxEtaBin *phi + eta; // phi- changes faster
00106        int rdo = x->rdo;
00107        if (tow<0 || tow>mxBtow || rdo<0 || rdo>mxBtow) return -101;
00108        
00109        mTower2rdo[ tow ] = rdo;    // returns rdo channel for given tower
00110        mRdo2tower[ rdo ] = tow;    // returns tower for given rdo channel
00111     }
00112   return 0; //OK
00113 
00114 }               
00115 
00116   
00117 
00118 /* ========================================
00119   ======================================== */
00120 float
00121 L2bemcGamma2009::sumET(int phi, int eta) {
00122   int tow = BtowGeom::mxEtaBin *((phi+BtowGeom::mxPhiBin)%BtowGeom::mxPhiBin) + ((eta+BtowGeom::mxEtaBin)%BtowGeom::mxEtaBin); // phi- changes faster
00123 
00124   const int maxTowers = BtowGeom::mxEtaBin * BtowGeom::mxPhiBin;
00125   int towPlusOne;
00126   float sum;
00127   sum=0;
00128   sum=wrkBtow_et[tow];
00129   towPlusOne = tow+1;
00130   towPlusOne%= maxTowers;
00131   sum+=wrkBtow_et[towPlusOne];
00132   
00133   //if(tow==0 || towPlusOne==0) {
00134     //printf("tow : %d, %d --> %f %f \n",tow, towPlusOne,wrkBtow_et[tow],wrkBtow_et[towPlusOne]);
00135   //}
00136   tow+=BtowGeom::mxEtaBin;
00137   tow%=maxTowers;
00138   
00139   sum+=wrkBtow_et[tow];
00140   
00141   towPlusOne = tow+1;
00142   towPlusOne%= maxTowers;
00143   sum+=wrkBtow_et[towPlusOne];
00144   
00145   //if(tow==0 || towPlusOne==0) {
00146     //printf("tow : %d, %d --> %f %f \n",tow, towPlusOne,wrkBtow_et[tow],wrkBtow_et[towPlusOne]);
00147   //}
00148   // printf("B sumET=%f\n",sum);
00149   return sum;
00150 }
00151   
00152 
00153 /* ========================================
00154   ======================================== */
00155 void 
00156 L2bemcGamma2009::computeUser(int token){
00157   // token range is guaranteed by virtual08-class
00158 
00159   /* 2x2 cluster finder:
00160      - can find clusters through the tower border in Phi
00161      - sorts hit towers by Et and selects Highest Seed Towers first
00162   */
00163   clearEvent(token);
00164 
00165   // ------ PROJECT INPUT LIST TO 2D ARRAY AND SCAN FOR SEED TOWERS ----
00166   int i;
00167   //  printf("L2-%s-compute: ---BTOW ADC list--- size=%d\n",getName(),*globEve_btow_hitSize);
00168 
00169   L2bemcGammaEvent2009 *btowEve=mBtow+token;
00170   const HitTower1 *hit=mEveStream_btow[token].get_hits();
00171   const int hitSize=mEveStream_btow[token].get_hitSize();
00172   for(i=0;i< hitSize;i++,hit++) {
00173     int tower=mRdo2tower[hit->rdo];
00174     wrkBtow_et[tower]=hit->et;
00175     if(hit->et<par_seedEtThres)continue;
00176     wrkBtow_tower_seed[wrkBtow_tower_seed_size++]=tower;
00177     //printf("A  seed TWID=%d \n",tower);
00178   }
00179   hA[2]->fill(hitSize);
00180   int seedTow=-1,seedEta=-1,seedPhi=-1;
00181   float clustET=0;
00182   btowEve->isFresh=L2bemcGammaEvent2009::kDataFresh;
00183 
00184   // printf("B nseed=%d\n",wrkBtow_tower_seed_size);
00185   // ----------- FIND 2x2 CLUSTER AROUND EVERY SEED -----
00186   for(i=0; i<wrkBtow_tower_seed_size;i++) {
00187     seedTow=wrkBtow_tower_seed[i];
00188     seedEta=seedTow%BtowGeom::mxEtaBin;
00189     seedPhi=seedTow/BtowGeom::mxEtaBin;
00190     
00191     //.... find first 2x2 above cluster thresh
00192     if (seedEta < BtowGeom::mxEtaBin) {
00193       clustET = sumET(seedPhi,seedEta);
00194       if(clustET>par_clusterEtThres) goto ACCEPT;
00195       clustET = sumET(seedPhi-1,seedEta);
00196       if(clustET>par_clusterEtThres) goto ACCEPT;
00197     }
00198     if (seedEta > 0 ) {
00199       clustET = sumET(seedPhi-1,seedEta-1);
00200       if(clustET>par_clusterEtThres) goto ACCEPT;
00201       clustET = sumET(seedPhi,seedEta-1);
00202       if(clustET>par_clusterEtThres) goto ACCEPT;
00203     }
00204   }
00205   //.... ABORT
00206   btowEve->resultBlob.clusterEt=0;
00207   btowEve->resultBlob.meanEtaBin=0;
00208   btowEve->resultBlob.meanPhiBin=0;
00209   btowEve->resultBlob.trigger=0;
00210   btowEve->seedTwID=-1;
00211   btowEve->clusterET=0;
00212   btowEve->seedET=0;
00213   return;
00214   
00215  ACCEPT:
00216   // printf("B acc clustET=%f\n",clustET);
00217   btowEve->clusterET=clustET;
00218   btowEve->seedET=wrkBtow_et[seedTow];
00219   btowEve->resultBlob.clusterEt=(unsigned char)(clustET*256.0/60.0);
00220   btowEve->resultBlob.meanEtaBin=seedEta;
00221   btowEve->resultBlob.meanPhiBin=seedPhi;
00222   btowEve->seedTwID=seedTow;
00223   btowEve->resultBlob.trigger=2;
00224   return;
00225 }
00226 
00227 
00228 
00229 
00230 /* ========================================
00231   ======================================== */
00232 bool 
00233 L2bemcGamma2009::decisionUser(int token, int *myL2Result){
00234   // INPUT: token + comput() results stored internally
00235   // OUTPUT: yes/now + pointer to  L2Result
00236 
00237 
00238   // get pointers to internal private event storage
00239   L2bemcGammaEvent2009 *btowEve=mBtow+token;
00240   bool triggerDecision=(btowEve->resultBlob.trigger>0);
00241 
00242   //prescaling is done in Virtual:
00243   if (par_RndAcceptPrescale>0) btowEve->resultBlob.trigger+=mRandomAccept;
00244   
00245   //...... some histos just for fun
00246   if(btowEve->isFresh>L2bemcGammaEvent2009::kDataFresh) mhN->fill(6); // stale data
00247 
00248   btowEve->isFresh++; // mark the data as  stale
00249 
00250   if(btowEve->resultBlob.trigger&2) {
00251     mhN->fill(15);
00252     hA[6]->fill((int)btowEve->clusterET);
00253     hA[7]->fill(btowEve->seedTwID);  
00254     hA[11]->fill((int)btowEve->seedET);
00255     hA[12]->fill((int)(100.*btowEve->seedET/btowEve->clusterET));
00256   }
00257 
00258   if(btowEve->resultBlob.trigger&1)  mhN->fill(16);
00259 
00260   memcpy(myL2Result,&(btowEve->resultBlob),sizeof(L2gammaResult2009));
00261   return triggerDecision;
00262 } 
00263 
00264 
00265 /* ========================================
00266   ======================================== */
00267 void
00268 L2bemcGamma2009::finishRunUser() {  /* called once at the end of the run */
00269   // do whatever you want, log-file & histo-file are still open
00270   
00271   if (mLogFile){ 
00272     fprintf(mLogFile,"finishRunUser-%s bhla bhla\n",getName());
00273   }
00274   
00275   const int *hist15Data = hA[7]->getData();
00276   for (int i = 0; i < BtowGeom::mxEtaBin; i++) {
00277     for (int j = 0; j < BtowGeom::mxPhiBin; j++) {
00278       //printf("i %d, j %d, iXj*n %d, data %d\n", i, j, i+j*BtowGeom::mxEtaBin, hist15Data[i+j*BtowGeom::mxEtaBin]);
00279       hA[14]->fillW(i,j, hist15Data[i+j*BtowGeom::mxEtaBin]);
00280     }
00281   }
00282   
00283 }
00284 
00285 
00286 //=======================================
00287 //=======================================
00288 void 
00289 L2bemcGamma2009::createHisto() {
00290   setMaxHist(16); // PMN added - histogram count does not seem to be initialiazed anywere.
00291   //memset(hA,0,sizeof(hA));
00292 
00293   hA[2]=new L2Histo(2,"BTOW-compute: #towers w/ energy /event; x: # BTOW towers; y: counts", 100); 
00294   hA[6]=new L2Histo(6,"BTOW-decision: accepted clust ...  ; x: ET(GeV)", 30);// title in initRun
00295 
00296   hA[7]=new L2Histo(7,"BTOW: accepted Seed Tower .....", 5000); // title in initRun
00297   hA[11]=new L2Histo(11,"BTOW: decision Cluster Seed Et; ET GeV", 30); // title in initRun
00298   hA[12]=new L2Histo(12,"BTOW: decision ;100*Seed Et/Cluster Et", 100); // title in initRun
00299 
00300   hA[14]=new L2Histo(14,"BTOW: hot tower projection", BtowGeom::mxEtaBin, BtowGeom::mxPhiBin); // title in initRun
00301 
00302   // printf("L2-%s::createHisto() done\n",getName());
00303 }
00304 
00305 //=======================================
00306 //=======================================
00307 void 
00308 L2bemcGamma2009::clearEvent(int token){
00309   memset(wrkBtow_et,0,sizeof(wrkBtow_et)); 
00310   wrkBtow_tower_seed_size=0; 
00311   memset(&(mBtow[token].resultBlob),0, sizeof(L2gammaResult2009));
00312 }
00313   
00314 /* ========================================
00315   ======================================== */
00316 void 
00317 L2bemcGamma2009::print2(){ // full , local ADC array
00318   int i;
00319   printf("pr2-%s: ---BTOW ADC 2D array, only non-zero\n",getName());
00320 
00321   for(i=0;i<mxBtow;i++) {
00322     if(wrkBtow_et[i]<=0) continue;
00323     int rdo=mTower2rdo[i];
00324     float et=wrkBtow_et[i];
00325     printf("  btow: tower=%4d  rdo=%4d   et=%.3f \n",i,rdo,et);
00326   }
00327 
00328 }
00329 
00330 /* ========================================
00331   ======================================== */
00332 void 
00333 L2bemcGamma2009::print3(){ // seed list
00334   int i;
00335   printf("pr3-%s: ---seed list, size=%d\n",getName(),wrkBtow_tower_seed_size);
00336 
00337   for(i=0;i<wrkBtow_tower_seed_size;i++) {
00338     int tower=wrkBtow_tower_seed[i];
00339     float et=wrkBtow_et[tower];
00340     printf("  btow: i=%4d  tower=%4d   et=%.3f \n",i,tower,et);
00341   }
00342 
00343 }
00344 
00345 #if 0
00346 
00347 /* ========================================
00348   ======================================== */
00349 void 
00350 L2bemcGamma2009::print4(int token, int hitSize){ // L2-algo input list
00351   int i;
00352   printf("print4 IS NOT Fully FUNCTIONAL **********************\n");
00353   printf("pr1-%s: ---BTOW Sorted ADC list--- size=%d\n",getName(),hitSize);
00354   //const HitTower *hit=globEve_btow_hit;
00355   for(i=0;i< hitSize;i++) {
00356     int adc=0;//(mEveStream_btow[token].get_hits()[wrkBtow_tower_index[i]]).adc;
00357     int rdo=0;//(mEveStream_btow[token].get_hits()[wrkBtow_tower_index[i]]).rdo;
00358     float et=wrkBtow_et[wrkBtow_tower_index[i]];
00359     float ene=0;//(mEveStream_btow[token].get_hits()[wrkBtow_tower_index[i]]).ene;
00360     printf("  tower=%2d ",wrkBtow_tower_index[i]);
00361     printf("  btow: i=%2d rdo=%4d  adc=%d  et=%.3f  ene=%.3f\n",i,rdo,adc,et,ene);
00362   }
00363 }
00364 
00365 #endif
00366 /**********************************************************************
00367   $Log: L2bemcGamma2009.cxx,v $
00368   Revision 1.1  2011/03/09 16:29:07  pibero
00369   Added L2gamma2009
00370 
00371   Revision 1.6  2008/01/30 21:56:40  balewski
00372   E+B high-enery-filter L2-algo fuly functional
00373 
00374   Revision 1.5  2008/01/30 00:47:17  balewski
00375   Added L2-Etow-calib
00376 
00377   Revision 1.4  2008/01/18 23:29:13  balewski
00378   now L2result is exported
00379 
00380   Revision 1.3  2008/01/17 23:15:51  balewski
00381   bug in token-addressed memory fixed
00382 
00383   Revision 1.2  2008/01/16 23:32:35  balewski
00384   toward token dependent compute()
00385 
00386   Revision 1.1  2007/12/19 02:30:18  balewski
00387   new L2-btow-calib-2008
00388 
00389 
00390  
00391 */
00392 
00393