StRoot: StTriggerUtilities/L2Emulator/L2upsilon/L2Upsilon2012.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 #ifdef  IS_REAL_L2  //in l2-ana  environment
00009   #include "../L2algoUtil/L2EmcDb2012.h"
00010   #include "../L2algoUtil/L2Histo.h"
00011 #else    //full path needed for cvs'd code
00012   #include "StTriggerUtilities/L2Emulator/L2algoUtil/L2EmcDb2012.h"
00013   #include "StTriggerUtilities/L2Emulator/L2algoUtil/L2Histo.h"
00014   #include "StTriggerUtilities/L2Emulator/L2algoUtil/L2EmcGeom2012.h"
00015 #endif
00016 
00017 #include "L2Upsilon2012.h"
00018 #include "read_in_geog.h"
00019 //extern inline void *memset(void *s, int c, size_t n);
00020 
00021 
00022 //=================================================
00023 //=================================================
00024 L2Upsilon2012::L2Upsilon2012(const char* name, const char *uid, L2EmcDb2012* db, L2EmcGeom2012 *geoX, char* outDir, int resOff)  :  L2VirtualAlgo2012( name, uid, db, outDir, true, false, resOff) { 
00025   /* called one per days
00026      all memory allocation must be done here
00027   */
00028 
00029   mGeom=geoX;
00030   if (!mGeom)
00031     criticalError("L2Upsilon2012 is broken -- can't find geom.");
00032 
00033   setMaxHist(32); // set upper range, I uses only 2^N -it is easier to remember
00034   createHisto();
00035 
00036   //------- self-consistency checks, should never fail
00037   if (!(sizeof(L2UpsilonResult2012)== L2UpsilonResult2012::mySizeChar))
00038     {
00039       criticalError("L2Upsilon2012 has failed consistency check 'sizeof(L2UpsilonResult2012)== L2UpsilonResult2012::mySizeChar'");
00040     }
00041 
00042   prescale = 1;
00043 }
00044 
00045 /* ========================================
00046   ======================================== */
00047 int 
00048 L2Upsilon2012::initRunUser( int runNo, int *rc_ints, float *rc_floats) {
00049         
00050   // unpack params from run control GUI
00051 
00052   par_prescale      =  rc_ints[0];
00053   fMaxDynamicMaskTowers   =rc_ints[1];
00054   fHowManyEventPerUpdateDynamicMask=rc_ints[2];
00055   par_RndAcceptPrescale=rc_ints[3]; //added by Ross.
00056 
00057   fL0SeedThreshold = rc_floats[0];
00058   fL2SeedThreshold = rc_floats[1];
00059   fMinL0ClusterEnergy = rc_floats[2];
00060   fMinL2ClusterEnergy = rc_floats[3];
00061   fMinInvMass         = rc_floats[4];
00062   fMaxInvMass         = rc_floats[5];
00063   fMaxCosTheta        = rc_floats[6];
00064   fHotTowerThreshold=rc_floats[7];
00065   fThresholdRatioOfHotTower=rc_floats[8];
00066 
00067   fHotTowerSeenTimesThreshold=int(fHowManyEventPerUpdateDynamicMask*fThresholdRatioOfHotTower);
00068 
00069 /*   suggested pars
00070   par_prescale      =  1;
00071   fMaxDynamicMaskTowers   =25;
00072   fHowManyEventPerUpdateDynamicMask=50;
00073 
00074   fL0SeedThreshold = 4.;
00075   fL2SeedThreshold = 2.5;
00076   fMinL0ClusterEnergy = 4.5;
00077   fMinL2ClusterEnergy = 3.;
00078   fMinInvMass         = 5;
00079   fMaxInvMass         = 20;
00080   fMaxCosTheta        = 0;
00081   fHotTowerThreshold=2.5;
00082   fThresholdRatioOfHotTower=0.5;
00083   
00084   fHotTowerSeenTimesThreshold=25;
00085 */
00086 
00087 /*
00088 ---  temp par values for test -----
00089   par_prescale =2;
00090   fL0SeedThreshold = 1.;
00091   fL2SeedThreshold = 1;
00092   fMinL0ClusterEnergy = 1;
00093   fMinL2ClusterEnergy = 1.;
00094   fMinInvMass         = 1;
00095   fMaxInvMass         = 20;
00096   fMaxCosTheta        = 1;
00097 //--------------------------------   
00098 */
00099 
00100         EventSeen=0;
00101         event_accept=0;
00102         memset(wrkDynamicMask_tower_stat,0,sizeof(wrkDynamicMask_tower_stat)); 
00103         wrkNumberOfMasked=0;
00104 
00105     // verify consistency of input params
00106   int kBad=0;
00107   kBad+=0x00001 * (fMaxCosTheta>1.0||fMaxCosTheta<-1.0);
00108   kBad+=0x00002 * (fL0SeedThreshold<fL2SeedThreshold);
00109   kBad+=0x00004 * (fMinL0ClusterEnergy<fMinL2ClusterEnergy);
00110   kBad+=0x00008 * (fMinInvMass>fMaxInvMass);
00111   kBad+=0x00010 * (fHowManyEventPerUpdateDynamicMask<fHotTowerSeenTimesThreshold);
00112   kBad+=0x00020 * (fMaxDynamicMaskTowers>100);
00113   kBad+=0x00040 * (fThresholdRatioOfHotTower>1);
00114 
00115   if (mLogFile) { 
00116     fprintf(mLogFile," L2%s upsilon algorithm initRun(R=%d), compiled: %s , %s\n params:\n",getName(),mRunNumber,__DATE__,__TIME__);
00117     fprintf(mLogFile," initRun() params checked for consistency, Error flag=%d\n",kBad);
00118     fprintf(mLogFile," prescale=%d \n",  par_prescale);
00119     fprintf(mLogFile," rndAccept prescale=%d \n",  par_RndAcceptPrescale);
00120     fprintf(mLogFile," L0SeedThreshold= %.2f  (GeV)\n",fL0SeedThreshold);
00121     fprintf(mLogFile," L2SeedThreshold= %.2f  (GeV)\n",fL2SeedThreshold);
00122     fprintf(mLogFile," MinL0ClusterEnergy= %.2f  (GeV)\n",fMinL0ClusterEnergy);
00123     fprintf(mLogFile," MinL2ClusterEnergy= %.2f  (GeV)\n",fMinL2ClusterEnergy);
00124     fprintf(mLogFile," MinInvMass= %.2f  (GeV)\n",fMinInvMass);
00125     fprintf(mLogFile," MaxInvMass= %.2f  (GeV)\n",fMaxInvMass);
00126     fprintf(mLogFile," MaxCosTheta= %.2f  \n",fMaxCosTheta);
00127     fprintf(mLogFile," MaxDynamicMaskedTowers= %d (should not greater than 100) \n",fMaxDynamicMaskTowers);
00128     fprintf(mLogFile," HotTowerThreshold= %f  (GeV) \n",fHotTowerThreshold);
00129     fprintf(mLogFile," HowManyEventsPerUpdateDynamicMask= %d \n",fHowManyEventPerUpdateDynamicMask);
00130     fprintf(mLogFile," fHotTowerSeenTimesThreshold= %d   \n",fHotTowerSeenTimesThreshold);
00131     fprintf(mLogFile," fThresholdRatioOfHotTower= %f  (count_hot / count_check) \n",fThresholdRatioOfHotTower);
00132   }
00133 
00134 
00135   // clear content of all histograms & token-dependet memory
00136   int i;
00137   for (i=0; i<mxHA;i++) if(hA[i])hA[i]->reset();
00138   memset(mBtow,0,sizeof(mBtow));
00139 
00140   if(kBad>0) return -1*kBad;    
00141   else if(kBad<0) return kBad;
00142  
00143   for ( int index=0; index<EmcDbIndexMax; index++ )
00144      {
00145        const L2EmcDb2012::EmcCDbItem *x = mDb->getByIndex(index);
00146        if ( x==0 ) continue;
00147        if ( !mDb->isBTOW(x) ) continue; 
00148        int sec = x->sec - 1;
00149        int sub = 8192; 
00150        sub = x->sub - 'a';
00151        int eta = x->eta - 1;
00152        int phi = BtowGeom::mxSubs *sec + sub;
00153        int tow = BtowGeom::mxEtaBin *phi + eta; // phi- changes faster
00154        int rdo = x->rdo;
00155        if (tow<0 || tow>mxBtow || rdo<0 || rdo>mxBtow) return -101;
00156        
00157        mTower2rdo[ tow ] = rdo;    // returns rdo channel for given tower
00158        mRdo2tower[ rdo ] = tow;    // returns tower for given rdo channel
00159 
00160                 int a,b,c,d;
00161                 sscanf(x->tube,"id%d-%d-%d-%d",&a,&b,&c,&d);
00162                 rdo2softID[x->rdo]=a;
00163 //              hA[20]->fill(a,x->rdo);
00164         }
00165   return 0; //OK
00166 
00167 }               
00168 
00169 
00170 /* ========================================
00171   ======================================== */
00172 void 
00173 L2Upsilon2012::computeUser(int token){
00174   clearEvent(token);
00175 
00176   EventSeen++;
00177   if(EventSeen%fHowManyEventPerUpdateDynamicMask==0) update_DynamicMask();
00178 
00179   L2UpsilonEvent2012 *btowEve=mBtow+token;
00180   const HitTower1 *hit=mEveStream_btow[token].get_hits();
00181   const int hitSize=mEveStream_btow[token].get_hitSize();
00182 
00183   for(int i=0;i< hitSize;i++,hit++) {
00184 //    int tower=mRdo2tower[hit->rdo];   
00185     int tower=rdo2softID[hit->rdo];
00186         wrkBtow_tower[i]=tower;
00187     wrkBtow_ene[tower]=hit->ene;     //    fill energy array    (all towers)
00188         if(wrkBtow_ene[tower] > fHotTowerThreshold)   wrkDynamicMask_tower_stat[tower]++;  //  wrkDynamicMask_tower_stat is needed to judge hot towers
00189         hA[2]->fill(int(wrkBtow_ene[tower]));
00190   }
00191   hA[1]->fill(hitSize);
00192 
00193 //  Dynamic Masking 
00194         for(int j=0;j<wrkNumberOfMasked;j++) 
00195                 wrkBtow_ene[wrkDynamicMasked_tower[j]]=0;
00196 
00197 // ---- store all L0 L2 candidates
00198   wrkNumberOfL0=0;
00199   wrkNumberOfL2=0;
00200   for(int i=0;i< hitSize;i++) {
00201         int tower=wrkBtow_tower[i];
00202         if(wrkBtow_ene[tower] > fL2SeedThreshold){
00203                 float clusterE=wrkBtow_ene[tower];
00204                 for(int k=0;k<geog_Nneighbors[tower];k++){
00205                         float Ek=wrkBtow_ene[geog_neighbor[tower][k]];
00206                         int Ecount=0;
00207                         for(int m=0;m<geog_Nneighbors[tower];m++){
00208                                 if(m==k) continue;
00209                                 if(Ek<wrkBtow_ene[geog_neighbor[tower][m]])Ecount++;
00210                         }
00211                         if(Ecount<2)clusterE+=Ek;  // top 2 cluster
00212                 }
00213 
00214                 if(clusterE>fMinL2ClusterEnergy) {   //   tower seed good && cluster good  for L2 threshold
00215                         wrkL2_seed_tower[wrkNumberOfL2]=tower;
00216                         wrkL2_seed_ClusterE[wrkNumberOfL2]=clusterE;
00217                         hA[6]->fill(int(wrkBtow_ene[tower]));
00218                         hA[8]->fill(int(clusterE));
00219                         wrkNumberOfL2++;
00220                         if(wrkBtow_ene[tower] > fL0SeedThreshold&&clusterE>fMinL0ClusterEnergy) {   //   tower seed good && cluster good  for L0 threshold
00221                                 wrkL0_seed_tower[wrkNumberOfL0]=tower;
00222                                 wrkL0_seed_ClusterE[wrkNumberOfL0]=clusterE;
00223                                 hA[5]->fill(int(wrkBtow_ene[tower]));
00224                                 hA[7]->fill(int(clusterE));
00225                                 wrkNumberOfL0++;
00226                         }
00227                 }
00228         }
00229   }
00230   hA[3]->fill(wrkNumberOfL2,wrkNumberOfL0);
00231 
00232 //  ------  pair L0 L2, calculate invMass
00233         for(int i=0;i<wrkNumberOfL0;i++){
00234                 int idL2,idL0;
00235                 float cosTheta,invMass;
00236                 idL0= wrkL0_seed_tower[i];
00237                 for(int j=0;j<wrkNumberOfL2;j++){
00238                         idL2= wrkL2_seed_tower[j];
00239                         if(idL0==idL2)continue;
00240                         cosTheta=geog_x[idL0]*geog_x[idL2]+geog_y[idL0]*geog_y[idL2]+geog_z[idL0]*geog_z[idL2];
00241                         if(cosTheta>fMaxCosTheta)  continue;
00242                         invMass = sqrt(2 *wrkL0_seed_ClusterE[i] * wrkL2_seed_ClusterE[j] * (1 - cosTheta));
00243                         if(fMinInvMass>invMass || invMass>fMaxInvMass) continue;
00244                         hA[10]->fill(int(invMass));
00245                         // accept !!!
00246                         btowEve->resultBlob.L0SeedTowerID=idL0;
00247                         btowEve->resultBlob.L2SeedTowerID=idL2;
00248                         if(wrkL0_seed_ClusterE[i]>25.5)wrkL0_seed_ClusterE[i]=25.5;
00249                         if(wrkL0_seed_ClusterE[j]>25.5)wrkL0_seed_ClusterE[j]=25.5;
00250                         if(invMass>25.5)invMass=25.5;
00251                         btowEve->resultBlob.energyOfL0Cluster=(unsigned char)(wrkL0_seed_ClusterE[i]*10);  //   256/25.6=10
00252                         btowEve->resultBlob.energyOfL2Cluster=(unsigned char)(wrkL2_seed_ClusterE[j]*10);  //   256/25.6=10
00253                         btowEve->resultBlob.invMass=(unsigned char)(invMass*10);  //   256/25.6=10
00254                         btowEve->resultBlob.trigger=true;
00255 //                      btowEve->resultBlob.numberOfL0SeedTowers=wrkNumberOfL0;
00256 //                      btowEve->resultBlob.numberOfL2SeedTowers=wrkNumberOfL2;
00257                         return;
00258                 }
00259         }
00261         btowEve->resultBlob.L0SeedTowerID=5555;
00262         btowEve->resultBlob.L2SeedTowerID=5555;
00263         btowEve->resultBlob.energyOfL0Cluster=0;
00264         btowEve->resultBlob.energyOfL2Cluster=0;
00265         btowEve->resultBlob.invMass=0;
00266         btowEve->resultBlob.trigger=false;
00267 //      btowEve->resultBlob.numberOfL0SeedTowers=wrkNumberOfL0;
00268 //      btowEve->resultBlob.numberOfL2SeedTowers=wrkNumberOfL2;
00269   return;
00270 }
00271 
00272 
00273 /* ========================================
00274   ======================================== */
00275 bool 
00276 L2Upsilon2012::decisionUser(int token, int *myL2Result){
00277   // INPUT: token + comput() results stored internally
00278   // OUTPUT: yes/no + pointer to  L2Result
00279 
00280   // get pointers to internal private event storage
00281   L2UpsilonEvent2012 *btowEve=mBtow+token;
00282 //  (*myL2Result)=&(btowEve->resultBlob);
00283   memcpy(myL2Result,&(btowEve->resultBlob),sizeof(L2UpsilonResult2012));
00284   if (par_prescale > 1) {
00285     prescale++;
00286         prescale%=par_prescale;
00287     if(prescale!=0)    btowEve->resultBlob.trigger=false;
00288   }
00289   btowEve->isFresh++; // mark the data as  stale
00290         hA[13]->fill(1);
00291 
00292         if(btowEve->resultBlob.trigger){
00293                 event_accept++;
00294                 hA[13]->fill(2);
00295                 hA[14]->fill(btowEve->resultBlob.L0SeedTowerID);
00296                 hA[15]->fill(btowEve->resultBlob.L2SeedTowerID);
00297                 hA[16]->fill(int((btowEve->resultBlob.energyOfL0Cluster)*0.1));
00298                 hA[17]->fill(int((btowEve->resultBlob.energyOfL2Cluster)*0.1));
00299                 hA[18]->fill(int((btowEve->resultBlob.invMass)*0.1));
00300         }
00301         else hA[13]->fill(3);
00302 
00303   return btowEve->resultBlob.trigger;
00304 } 
00305 
00306 
00307 /* ========================================
00308   ======================================== */
00309 void
00310 L2Upsilon2012::finishRunUser() {  /* called once at the end of the run */
00311   // do whatever you want, log-file & histo-file are still open
00312   
00313   if (mLogFile){ 
00314     fprintf(mLogFile,"finishRunUser - %s \n",getName());
00315     fprintf(mLogFile,"EventSeen = %d \n",EventSeen);
00316     fprintf(mLogFile,"event accept = %d \n",event_accept);
00317   }
00318   
00319 }
00320 
00321 
00322 //=======================================
00323 //=======================================
00324 void 
00325 L2Upsilon2012::createHisto() {
00326   setMaxHist(32); // PMN added - histogram count does not seem to be initialiazed anywere.
00327   //memset(hA,0,sizeof(hA));
00328 
00329         hA[1]=new L2Histo(1,"Upsilon:# of hits (no L0 required)", 200);
00330         hA[2]=new L2Histo(2,"Upsilon:energy for all hits (GeV)  (no L0 required)", 20);
00331         hA[3]=new L2Histo(3,"Upsilon: # of L0 candidates vs. # of L2 candidates (no L0 required)", 50,50);
00332 
00333         hA[5]=new L2Histo(5,"Upsilon: L0 seeds Energy (GeV) (no L0 required)", 20);
00334         hA[6]=new L2Histo(6,"Upsilon: L2 seeds Energy (GeV) (no L0 required)", 20);
00335 
00336         hA[7]=new L2Histo(7,"Upsilon: L0 cluster Energy (GeV) (no L0 required)", 25);
00337         hA[8]=new L2Histo(8,"Upsilon: L2 cluster Energy (GeV) (no L0 required)", 25);
00338 
00339         hA[10]=  new L2Histo(10,"Upsilon: InvMass of accepted pair (GeV) (no L0 required)", 25);
00340         hA[11]=new L2Histo(11,"Upsilon: # of masked towers  (no L0 required)", 100);
00341         hA[12]=new L2Histo(12,"Upsilon: masked tower softid  (no L0 required) ", 4801);
00342 
00343         hA[13] = new L2Histo(13,"Upsilon: L0 counter (1=L0 fired, 2=L0 fired+L2 accepted, 3=L0 fired+L2 NOT accepted)",5);
00344         hA[14] = new L2Histo(14,"Upsilon: L0 seed tower id (L0 fired+L2 accepted)",4801);
00345         hA[15] = new L2Histo(15,"Upsilon: L2 seed tower id (L0 fired+L2 accepted)",4801);
00346         hA[16]=  new L2Histo(16,"Upsilon: L0 cluster Energy (GeV) (L0 fired+L2 accepted)", 25);
00347         hA[17]=  new L2Histo(17,"Upsilon: L2 cluster Energy (GeV) (L0 fired+L2 accepted)", 25);
00348         hA[18]=  new L2Histo(18,"Upsilon: InvMass of accepted pair (GeV) (L0 fired+L2 accepted)", 25);
00349 
00350 //      hA[20]=new L2Histo(20,"id",4801,4801);
00351 
00352 }
00353 
00354 //=======================================
00355 //=======================================
00356 void 
00357 L2Upsilon2012::clearEvent(int token){
00358   memset(wrkBtow_ene,0,sizeof(wrkBtow_ene)); 
00359   memset(wrkBtow_tower,0,sizeof(wrkBtow_tower)); 
00360   memset(wrkL2_seed_tower,0,sizeof(wrkL2_seed_tower)); 
00361   memset(wrkL2_seed_ClusterE,0,sizeof(wrkL2_seed_ClusterE)); 
00362   memset(wrkL0_seed_tower,0,sizeof(wrkL0_seed_tower)); 
00363   memset(wrkL0_seed_ClusterE,0,sizeof(wrkL0_seed_ClusterE)); 
00364   memset(&(mBtow[token].resultBlob),0, sizeof(L2UpsilonResult2012));
00365   wrkNumberOfL2=wrkNumberOfL0=0;
00366 }
00367 
00368 void
00369 L2Upsilon2012::update_DynamicMask(){
00370         int count=0;
00371         for(int i=0;i<mxBtow;i++){
00372                 if(wrkDynamicMask_tower_stat[i]>fHotTowerSeenTimesThreshold && count<fMaxDynamicMaskTowers){   // tower i is hot
00373                         wrkDynamicMasked_tower[count]=i;
00374                         hA[12]->fill(i);
00375                         count++;
00376                 }
00377         }
00378         memset(wrkDynamicMask_tower_stat,0,sizeof(wrkDynamicMask_tower_stat)); 
00379         wrkNumberOfMasked=count;
00380         hA[11]->fill(count);
00381 }