StRoot: StTriggerUtilities/L2Emulator/L2algoUtil/L2btowCalAlgo12.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   $Id: L2btowCalAlgo12.cxx,v 1.6 2012/03/21 18:18:03 jml Exp $
00009   \author Jan Balewski, MIT, 2009 
00010  *****************************************************
00011   Descripion: 
00012   calibrates Barrel towers, result is used by other L2-algo
00013  *****************************************************/
00014 
00015 
00016 #ifdef  IS_REAL_L2  //in l2-ana  environment
00017   #include "../L2algoUtil/L2EmcDb2012.h"
00018   #include "../L2algoUtil/L2Histo.h"
00019 #else
00020   #include "L2EmcDb2012.h"
00021   #include "L2Histo.h"
00022   #include "L2EmcGeom2012.h"
00023 #endif
00024 
00025 #include "L2btowCalAlgo12.h"
00026 
00027 /************************
00028 the _only_ copy of L2event for all L2-algos
00029 **************************/
00030 L2eventStream2012  globL2eventStream2012;
00031 /*************************/
00032 
00033 //=================================================
00034 //=================================================
00035 L2btowCalAlgo12::L2btowCalAlgo12(const char* name, const char *uid, L2EmcDb2012* db, L2EmcGeom2012 *geoX, char* outDir, int resOff)  :  L2VirtualAlgo2012( name, uid,  db,  outDir, true, false, resOff) { 
00036   /* called one per days
00037      all memory allocation must be done here
00038   */
00039 
00040   geom=geoX; 
00041   if (!geom)
00042     criticalError("L2btowCalAlgo is broken -- can't find geom.");
00043 
00044   setMaxHist(32);
00045   createHisto();
00046 
00047   // initialize BTOW-Calibrated-data
00048   int k;
00049   for(k=0;k<L2eventStream2012::mxToken;k++){
00050     L2BtowCalibData12 & btowCalibData=globL2eventStream2012.btow[k];  
00051     btowCalibData.nInputBlock=0;
00052     btowCalibData.hitSize=0;
00053   }
00054 
00055   LOG("JEFF", "sizeof L2 event stream: %d", sizeof(L2eventStream2012));
00056 
00057  }
00058 
00059 /* ========================================
00060   ======================================== */
00061 int 
00062 L2btowCalAlgo12::initRunUser( int runNo, int *rc_ints, float *rc_floats) {
00063 
00064   par_adcMask=(unsigned short)(-0x10); 
00065   par_pedOff=0x10/2; //hex must match with above 
00066  
00067   // unpack params from run control GUI
00068   par_dbg       =  rc_ints[0];
00069   par_gainType  =  rc_ints[1];
00070   par_nSigPed   =  rc_ints[2];
00071 
00072   par_twEneThres = rc_floats[0];
00073   par_hotEtThres = rc_floats[1];;
00074 
00075   // verify consistency of input params
00076   int kBad=0;
00077   kBad+=0x00001 * (par_gainType<kGainZero || par_gainType>kGainOffline);
00078   kBad+=0x00002 * (par_nSigPed<2 || par_nSigPed>5);
00079   kBad+=0x00004 * (par_twEneThres<0.1 ||  par_twEneThres>1.5);
00080 
00081   if (mLogFile) { 
00082     fprintf(mLogFile,"L2%s algorithm initRun(R=%d), compiled: %s , %s\n params:\n",getName(),mRunNumber,__DATE__,__TIME__);
00083     fprintf(mLogFile," - use BTOW=%d,  gain Ideal=%d or Offline=%d, debug=%d\n",
00084             par_gainType>=kGainIdeal, par_gainType==kGainIdeal, par_gainType==kGainOffline, par_dbg);
00085     fprintf(mLogFile," - thresholds: ADC-ped> %d*sigPed .AND. energy>%.2f GeV \n", par_nSigPed, par_twEneThres);
00086 
00087     fprintf(mLogFile," - hot tower thresholds:  ET/GeV=%.2f\n",par_hotEtThres);
00088     fprintf(mLogFile,"initRun() params checked for consistency, Error flag=0x%04x\n",kBad);
00089   }
00090 
00091   // initialize BTOW-Calibrated-data                                                                                                         
00092   int k;
00093   for(k=0;k<L2eventStream2012::mxToken;k++){
00094     L2BtowCalibData12 & btowCalibData=globL2eventStream2012.btow[k];
00095     btowCalibData.nInputBlock=0;
00096     btowCalibData.hitSize=0;
00097   }
00098 
00099   
00100   if(kBad) return -kBad;
00101 
00102   // clear content of all histograms
00103   int i;
00104   for (i=0; i<mxHA;i++) if(hA[i])hA[i]->reset();
00105 
00106   // update title of histos
00107   char txt[1000];
00108   sprintf(txt,"BTOW tower, E_T>%.2f GeV (input); x: BTOW RDO index=chan*30+fiber; y: counts",par_hotEtThres);
00109   hA[10]->setTitle(txt);
00110   
00111   sprintf(txt,"BTOW tower, Et>%.2f GeV (input); x: BTOW softID",par_hotEtThres);
00112   hA[11]->setTitle(txt);
00113   sprintf(txt,"BTOW tower, Et>%.2f GeV (input); x: eta bin, [-1,+1];  y: phi bin ~ TPC sector",par_hotEtThres);
00114   hA[12] ->setTitle(txt);
00115   
00116   sprintf(txt,"#BTOW towers / event , Et>%.2f GeV; x: # BTOW towers; y: counts",par_hotEtThres);
00117   hA[14] ->setTitle(txt);
00118   
00119   // re-caluclate geometry properties
00120   geom->btow.clear(); 
00121   int nB=0; /* counts # of unmasekd towers */ 
00122   int nBg=0; /* counts # of reasonable calibrated towers */ 
00123   int nEneThr=0, nPedThr=0; //BTOW count # of towers above & below threshold
00124   if(par_gainType>=kGainIdeal)  // this disables the whole loop below
00125   for(i=0; i<EmcDbIndexMax; i++) {
00126     const L2EmcDb2012::EmcCDbItem *x=mDb->getByIndex(i);
00127     if(mDb->isEmpty(x)) continue;  /* dropped not mapped  channels */
00128     /*....... B A R R E L  .................*/
00129     if (!mDb->isBTOW(x) ) continue; /* drop if not BTOW */
00130     if(x->fail) continue;          /* dropped masked channels */
00131     if(x->gain<=0) continue;       /* dropped uncalibrated towers , tmp */
00132     nB++;   
00133 
00134     float adcThres=x->ped+par_nSigPed* fabs(x->sigPed);
00135     float otherThr=x->ped+par_twEneThres*x->gain;
00136  
00137     if(adcThres<otherThr) { //use energy threshold if higher
00138       adcThres=otherThr;
00139       nEneThr++;
00140     } else {
00141       nPedThr++;
00142     }
00143     
00144     /* use rdo index to match RDO order in the ADC data banks */    
00145     if(x->eta<=0 || x->eta>BtowGeom::mxEtaBin) return -90;
00146     int ietaTw= (x->eta-1); /* correct */
00147 
00148     // use ideal gains for now, hardcoded
00149 
00150     if (par_gainType!=kGainIdeal) return -102;    
00151     geom->btow.gain2Ene_rdo[x->rdo]=geom->btow.idealGain2Ene[ietaTw];
00152     geom->btow.gain2ET_rdo[x->rdo]=geom->getIdealAdc2ET();
00153     
00154     geom->btow.thr_rdo[x->rdo]=(int) (adcThres);
00155     geom->btow.ped_rdo[x->rdo]=(int) (x->ped);
00156     geom->btow.ped_shifted_rdo[x->rdo]=(unsigned short)(par_pedOff - x->ped);
00157     nBg++;      
00158   }
00159   
00160   if (mLogFile) {
00161     fprintf(mLogFile,"  found  towers working=%d calibrated=%d, based on ASCII DB\n",nB,nBg);
00162     fprintf(mLogFile,"  thresh defined by energy=%d  or NsigPed=%d \n",nEneThr, nPedThr);
00163   }
00164 
00165   return 0; //OK
00166 
00167 
00168 }               
00169 
00170 /* ========================================
00171   ======================================== */
00172 void 
00173 L2btowCalAlgo12::calibrateBtow(int token, int bemcIn, ushort *rawAdc){
00174   // Btow calibration is a special case,  must have one exit  at the end
00175   computeStart();
00176 
00177   //LOG("JEFF", "Token=%d",token);
00178   token&=L2eventStream2012::tokenMask; // only protect against bad token, Gerard's trick
00179   //LOG("JEFF", "Token=%d",token);
00180 
00181   //...... token is valid no w ........
00182   L2BtowCalibData12 & btowCalibData=globL2eventStream2012.btow[token];  
00183   btowCalibData.nInputBlock++;
00184   
00185   // clear data for this token from previous event
00186   btowCalibData.hitSize=0;
00187 
00188   int nTower=0; /* counts mapped & used ADC channels */
00189   int nHotTower=0;
00190   if(bemcIn && par_gainType>kGainZero) { // EVEVEVEVEVE
00191     // ............process this event ...............
00192     short rdo;
00193     int low_noise_adc;
00194     int adc; // pedestal subtracted
00195     float et;
00196     float low_noise_et; // bit-chopped
00197     ushort *thr=geom->btow.thr_rdo;
00198     ushort *ped=geom->btow.ped_rdo;
00199     ushort *ped_shifted=geom->btow.ped_shifted_rdo;
00200     float *gain2ET=geom->btow.gain2ET_rdo;
00201     float *gain2Ene=geom->btow.gain2Ene_rdo;
00202     HitTower1 *hit=btowCalibData.hit;
00203     for(rdo=0; rdo<BtowGeom::mxRdo; rdo++){
00204       if(rawAdc[rdo]<thr[rdo])continue;
00205 
00206       //LOG("JEFF", "rawAdc[%d] = %d  thr=%d  ped=%d",rdo,rawAdc[rdo],thr[rdo],ped[rdo],BtowGeom::mxRdo);
00207 
00208       if(nTower>=L2BtowCalibData12::mxListSize) break; // overflow protection
00209       adc=rawAdc[rdo]-ped[rdo];  //do NOT correct for common pedestal noise
00210       et=adc/gain2ET[rdo]; 
00211       low_noise_adc=(rawAdc[rdo]+ped_shifted[rdo]) & par_adcMask;//note that ped_shifted is defined as pedOffset-ped.
00212       low_noise_et=low_noise_adc/gain2ET[rdo];
00213       hit->rdo=rdo;
00214       hit->adc=adc;
00215       hit->et=et;
00216       hit->low_noise_et=low_noise_et;
00217       hit->ene=adc/gain2Ene[rdo]; 
00218       hit++;
00219       nTower++; 
00220       // only monitoring
00221       if(et >par_hotEtThres) {
00222         hA[10]->fill(rdo);
00223         nHotTower++;
00224       }
00225     }
00226     btowCalibData.hitSize=nTower;
00227     
00228     // QA histos
00229     hA[13]->fill(nTower);
00230     hA[14]->fill(nHotTower);
00231     if(nTower>=L2BtowCalibData12::mxListSize) mhN->fill(5); // was overflow
00232   
00233   } // EVEVEVEVEVE
00234 
00235   LOG(DBG, "hits = %d   hot=%d",nTower,nHotTower);
00236 
00237   // debugging should be off for any time critical computation
00238   if(par_dbg>0){
00239     LOG(DBG,"L2-%s-compute: set adcL size=%d, get=%d\n",getName(),nTower,999); //tmp
00240     LOG(DBG,"dbg=%s: found  nTw=%d\n",getName(),nTower);
00241     if(par_dbg>0)   print0();
00242     printCalibratedData(token);
00243   } 
00244   
00245   computeStop(token);
00246 
00247 } 
00248 
00249 /* ======================================== 
00250    ======================================== */
00251 
00252 void
00253 L2btowCalAlgo12::clear(int token){
00254   token&=L2eventStream2012::tokenMask; // only protect against bad token, Gerard's trick
00255   //...... token is valid no w ........
00256   L2BtowCalibData12 & btowCalibData=globL2eventStream2012.btow[token];
00257   btowCalibData.hitSize=0;
00258   return;
00259 }
00260 
00261 /* ========================================
00262   ======================================== */
00263 void 
00264 L2btowCalAlgo12::computeUser(int token ){
00265 
00266   LOG(CRIT,"L2btowCalAlgo12::computeUser has been called and should not have been.  Serious problem in L2");
00267 }
00268 
00269 
00270 /* ========================================
00271   ======================================== */
00272 void
00273 L2btowCalAlgo12::finishRunUser() {  /* called once at the end of the run */
00274   // do whatever you want, log-file & histo-file are still open
00275 
00276   //report mean and RMS of #towers/event:
00277   int mean=0, RMS=0;
00278   hA[13]->findMean(&mean,&RMS);
00279   if (mLogFile){
00280     fprintf(mLogFile,"#BTOW_nonzero_towers_per_event: mean=%d, rms=%d\n",mean,RMS);
00281   }
00282   // search for hot tower, re-project histos vs. other representations
00283  int bHotSum=1,bHotId=-1;
00284  const int *data20=hA[10]->getData();
00285  const L2EmcDb2012::EmcCDbItem *xB=mDb->getByIndex(402); // some wired default?
00286   
00287   int i;
00288   for(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 softId=atoi(x->tube+2);
00293     int ieta= (x->eta-1);
00294     int iphi= (x->sec-1)*BtowGeom::mxSubs + x->sub-'a' ;
00295     hA[11]->fillW(softId,data20[x->rdo]);
00296     hA[12]->fillW(ieta, iphi,data20[x->rdo]);
00297     if(bHotSum<data20[x->rdo]) {
00298       bHotSum=data20[x->rdo];
00299       bHotId=softId;
00300       xB=x;
00301     }
00302   }
00303   
00304   if (mLogFile){
00305     fprintf(mLogFile,"#BTOW_hot tower _candidate_ (bHotSum=%d of %d eve) :, softID %d , crate %d , chan %d , name %s\n",bHotSum,mEventsInRun,bHotId,xB->crate,xB->chan,xB->name);
00306   }
00307   
00308   //...... QA tokens .....
00309   int tkn1=99999, tkn2=0; // min/max token
00310   int nTkn=0;
00311   int tkn3=-1, nTkn3=-1; // most often used token
00312   
00313   int k;
00314   for(k=0;k<L2eventStream2012::mxToken;k++){
00315     L2BtowCalibData12 & btowCalibData=globL2eventStream2012.btow[k];  
00316     if(btowCalibData.nInputBlock==0) continue;
00317     hA[1]->fillW(k,btowCalibData.nInputBlock);
00318     if(nTkn3<btowCalibData.nInputBlock){
00319       nTkn3=btowCalibData.nInputBlock;
00320       tkn3=k;
00321     }
00322 
00323     nTkn++;
00324     if(tkn1>k) tkn1=k;
00325     if(tkn2<k) tkn2=k;
00326   }
00327   if (mLogFile){
00328     fprintf(mLogFile,"#BTOW_token_QA:  _candidate_ hot token=%d used %d for %d events, token range [%d, %d], used %d tokens\n",tkn3,nTkn3,mEventsInRun,tkn1,tkn2,nTkn);
00329   }
00330 
00331 }
00332 
00333 
00334 //=======================================
00335 //=======================================
00336 void 
00337 L2btowCalAlgo12::createHisto() {
00338   memset(hA,0,mxHA*sizeof(L2Histo*));
00339   //token related spectra
00340   hA[1]=new  L2Histo(1,"L2-btow-calib: seen tokens;  x:  token value; y: events ",L2eventStream2012::mxToken);
00341   
00342   // BTOW  raw spectra
00343   hA[10]=new L2Histo(10,"btow hot tower 1", BtowGeom::mxRdo); // title upadted in initRun
00344   hA[11]=new L2Histo(11,"btow hot tower 2", BtowGeom::mxRdo); // title upadted in initRun
00345   hA[12]=new L2Histo(12,"btow hot tower 3", BtowGeom::mxEtaBin,BtowGeom::mxPhiBin); // title upadted in initRun  
00346   hA[13]=new L2Histo(13,"BTOW #tower w/ energy /event; x: # BTOW towers; y: counts", 200); 
00347   hA[14]=new L2Histo(14,"# hot towers/event", 100); 
00348 
00349 }
00350 
00351 
00352 
00353 /* ========================================
00354   ======================================== */
00355 void 
00356 L2btowCalAlgo12::print0(){ // full raw input  ADC array
00357   // empty
00358  }
00359 
00360 
00361 /****************************************************
00362   $Log: L2btowCalAlgo12.cxx,v $
00363   Revision 1.6  2012/03/21 18:18:03  jml
00364   got rid of printfs from 2012 files
00365 
00366   Revision 1.5  2011/10/19 16:12:10  jml
00367   more 2012 stuff
00368 
00369   Revision 1.4  2011/10/19 15:58:06  jml
00370   more compile offline
00371 
00372   Revision 1.3  2011/10/19 15:39:42  jml
00373   2012
00374 
00375   Revision 1.2  2011/10/19 14:34:23  jml
00376   added fakeRtsLog.h to turn log statements into printfs
00377 
00378   Revision 1.1  2011/10/18 15:11:41  jml
00379   adding 2012 algorithms
00380 
00381   Revision 1.3  2010/04/18 02:53:35  pibero
00382   Fixed memset bug.
00383 
00384   Revision 1.2  2010/04/17 16:42:09  pibero
00385   *** empty log message ***
00386 
00387   Revision 1.5  2008/02/01 00:16:40  balewski
00388   add mxListSize to BTOW/ETOW calibration
00389 
00390   Revision 1.4  2008/01/30 00:47:15  balewski
00391   Added L2-Etow-calib
00392 
00393   Revision 1.3  2008/01/18 23:29:13  balewski
00394   now L2result is exported
00395 
00396   Revision 1.2  2008/01/16 23:32:34  balewski
00397   toward token dependent compute()
00398 
00399   Revision 1.1  2007/12/19 02:30:18  balewski
00400   new L2-btow-calib-2008
00401 
00402   Revision 1.1  2007/11/19 22:18:25  balewski
00403   most L2algos provide triggerID's
00404 
00405  
00406 */
00407 
00408