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