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