StRoot: StTriggerUtilities/L2Emulator/L2algoUtil/L2btowCalAlgo08.cxx Source File

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