StRoot: StTriggerUtilities/L2Emulator/L2pedAlgo/L2pedAlgo.cxx Source File

00001 #include <stdio.h>
00002 #include <assert.h>
00003 #include <string.h>
00004 #include <stdlib.h>
00005 
00006 /*********************************************************************
00007  * $Id: L2pedAlgo.cxx,v 1.11 2009/11/19 15:48:46 balewski Exp $
00008  * \author Jan Balewski, IUCF, 2006 
00009  *********************************************************************
00010  * Descripion:
00011  * pedestal algo in L2 , for BTOW & ETOW
00012  *********************************************************************
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 "StTriggerUtilities/L2Emulator/L2algoUtil/L2EmcDb.h"
00021   #include "StTriggerUtilities/L2Emulator/L2algoUtil/L2Histo.h"
00022 #endif
00023 
00024 #include "L2pedAlgo.h"
00025 #include "L2pedResults2006.h"
00026 
00027 //=================================================
00028 //=================================================
00029 L2pedAlgo::L2pedAlgo(const char* name, L2EmcDb* db, char* outDir, int resOff) 
00030   :  L2VirtualAlgo( name,  db,  outDir, resOff) { 
00031   /* called one per days
00032      all memory allocation must be done here
00033   */
00034 
00035   par_pedSubtr=false;
00036   par_saveBinary=false;
00037   par_dbg=0;
00038   par_prescAccept=0;
00039 
00040   int i;
00041   for(i=0;i<MaxBtowRdo;i++) {
00042     char tit[100];
00043     sprintf(tit,"BTOW ADC for rdo=%d; ADC+%d",i,-minAdc);
00044     btowAdc[i]=new L2Histo(20000+i,tit,maxAdc-minAdc+1); 
00045   }
00046   
00047   for(i=0;i<MaxEtowRdo;i++) {
00048     char tit[100];
00049     sprintf(tit,"ETOW ADC for rdo=%d; ADC ",i);
00050     etowAdc[i]=new L2Histo(10000+i,tit,maxAdc-minAdc+1); 
00051  }
00052 
00053   // aux histos
00054   memset(hA,0,sizeof(hA));
00055   hA[10]=new   L2Histo(10, (char*)"total event counter; x=cases",6);
00056   hA[11]=new   L2Histo(11, (char*)"L2 time used per input event;  x: time (CPU 20*kTics); y: events ",500);
00057 
00058   // BTOW  raw spectra
00059   hA[20]=new   L2Histo(20, (char*)"BTOW pedRes  Y=ADC-DBped ; x: chan + 160*crate", 4800);
00060   hA[21]=new   L2Histo(21, (char*)"BTOW pedRes  Z=ADC-DBped, saturated @ |3|; x:  etaBin ,[-1,+1];  y: phi bin ~sector",40,120);
00061 
00062   // ETOW  raw spectra
00063   hA[30]=new   L2Histo(30, (char*)"ETOW pedRes Y=ADC-DBped ; x: chan + 128*crate", 768);
00064   hA[31]=new   L2Histo(31, (char*)"ETOW pedRes Z=ADC-DBped, saturated @ |3|; x: 12 - Endcap etaBin ,[+1,+2];  y: phi bin ~sector",12,60);
00065  
00066 
00067   printf("L2pedAlgo instantiated, logPath='%s'\n",mOutDir);
00068 
00069 }
00070 
00071 /*========================================
00072   ======================================== */
00073 int
00074 L2pedAlgo::initRun(int runNo, int *rc_ints, float *rc_floats) {
00075   //myName is not used.
00076   // update DB if run # has changed
00077   //printf("aaa L2pedAlgo::initRun runNo=%d\n",runNo);
00078   if(mDb->initRun(runNo)) return -27; 
00079   // DB must be initialized prior to lookup tables
00080 
00081   // unpack input params
00082 
00083 
00084   par_pedSubtr  =rc_ints[0]!=0;
00085   par_speedFact =rc_ints[1];
00086   par_saveBinary=rc_ints[2]!=0;
00087   par_dbg       =rc_ints[3];
00088   par_prescAccept=rc_ints[4];
00089 
00090   if(par_prescAccept<0) par_prescAccept=0; // prescale can't be negative
00091   //note speedFactor can be only powers of 2, range [1-256]
00092   if(par_speedFact<1) par_speedFact=1;
00093   if(par_speedFact>2) {// ASSURE ONLY POWERS OF 2, round down
00094     if(par_speedFact<4) par_speedFact=2;
00095     else if(par_speedFact<8) par_speedFact=4;
00096     else if(par_speedFact<16) par_speedFact=8;
00097     else if(par_speedFact<32) par_speedFact=16;
00098     else if(par_speedFact<64) par_speedFact=32;
00099     else if(par_speedFact<192) par_speedFact=64;
00100     else par_speedFact=192;
00101   }
00102   // can't be bigger, since there is only 6 etow crates
00103 
00104   s_lastB=s_lastE=0;
00105   s_stepB=MaxBtowRdo/ par_speedFact ;
00106   s_stepE=MaxEtowRdo/ par_speedFact ;
00107 
00108 
00109   /* .... clear content */
00110   memset(db_btowPed,       0,sizeof(db_btowPed));
00111   memset(db_etowPed,       0,sizeof(db_etowPed));
00112   
00113 
00114   int i;  
00115   int nBtowOk=0, nEtowOk=0;
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     if (mDb->isBTOW(x) ) {
00120           db_btowPed[x->rdo]=(int) (x->ped);
00121           nBtowOk++;
00122     } else if (mDb->isETOW(x) ) {
00123       db_etowPed[x->rdo]=(int) (x->ped);
00124       nEtowOk++;
00125     }
00126   }
00127   
00128   printf("L2ped algorithm init()... params:\n  dbg=%d, pedSubtr=%d saveBinHist=%d  speedFact=%d prescAccept=%d\n  mapped channels: nBtow=%d nEtow=%d\n",par_dbg,par_pedSubtr,par_saveBinary,par_speedFact,par_prescAccept,nBtowOk,nEtowOk);
00129 
00130   for(i=0;i<MaxBtowRdo;i++)  btowAdc[i]->reset();
00131   for(i=0;i<MaxEtowRdo;i++)  etowAdc[i]->reset();
00132   for(i=0;i<mxHA;i++) if(hA[i]) hA[i]->reset();
00133 
00134   nInp=0;
00135   run_number=runNo;
00136   return 0; // OK
00137 }               
00138 
00139 
00140 /*========================================
00141   ======================================== */
00142 bool
00143 L2pedAlgo::doEvent(int L0trg, int inpEveId, TrgDataType* trgData, 
00144                    int bemcIn, ushort *bemcData,
00145                    int eemcIn, ushort *eemcData){
00146   // not used: L0trg, inpEveId
00147   mAccept=true; // it never aborts
00148   /* STRICT TIME BUDGET  START ...., well a bit relaxed for this algo*/
00149   rdtscl_macro(mEveTimeStart);
00150   nInp++;
00151   hA[10]->fill(0);
00152   if(par_prescAccept>0) { // value=1 accepts 100% at maximal speed
00153     if((rand()>>4) % par_prescAccept ) return false;
00154     hA[10]->fill(5);
00155     return true;
00156   }
00157 
00158   myTrigData=trgData;
00159   
00160   /* *****************************************
00161      the code below is NOT optimized for speed
00162      one would need to make loop indexed by RDO to speed it up
00163      not worth for pedestal calculation
00164      This version takes 2,500 kTicks/eve --> ~1.5msec
00165      ******************************************** */
00166 
00167   short rdo;  
00168   if( bemcIn ) {
00169     short first=s_lastB%MaxBtowRdo;
00170     s_lastB=first+s_stepB;
00171     if(first==0) hA[10]->fill(1);
00172     // printf("B: f=%d l=%d\n",first,s_lastB);
00173     
00174     for(rdo=first; rdo<s_lastB; rdo++){
00175       int adc=bemcData[rdo];
00176       if(par_pedSubtr) adc-=db_btowPed[rdo];
00177       btowAdc[rdo]->fill(adc-minAdc);
00178       // printf("B: rdo=%d raw=%d ped=%d val=%d\n",rdo,bemcData[rdo],db_etowPed[rdo],adc-minAdc);
00179     }
00180   }
00181   
00182   
00183   if( eemcIn ) {
00184     short first=s_lastE%MaxEtowRdo;
00185     s_lastE=first+s_stepE;
00186     if(first==0) hA[10]->fill(2);
00187     //printf("E: f=%d l=%d\n",first,s_lastE);
00188     
00189     for(rdo=first; rdo<s_lastE; rdo++){
00190       int adc=eemcData[rdo];
00191       if(par_pedSubtr) adc-=db_etowPed[rdo];
00192       etowAdc[rdo]->fill(adc-minAdc);
00193       //    printf("E: rdo=%d raw=%d ped=%d val=%d\n",rdo,eemcData[rdo],db_etowPed[rdo],adc-minAdc);
00194     }
00195   }
00196   
00197   //...this  event will be accepted
00198   L2pedResults2006 out; // all output bits lump together
00199   memset(&out,0,sizeof(out));
00200   
00201   out.int0.decision=
00202     ( (bemcIn>0)    <<3 ) +
00203     ( (eemcIn>0)    <<4 ) +
00204   
00205     (  par_pedSubtr <<6 ) ;
00206  
00207   rdtscl_macro(mEveTimeStop);
00208   mEveTimeDiff=mEveTimeStop-mEveTimeStart;
00209   int  kTick=mEveTimeDiff/1000;
00210   hA[11]->fill(kTick/20);
00211   
00212   //  printf("jkTick=%d\n",kTick);
00213   const ushort maxKT=30000;
00214   out.int0.kTick=  kTick>maxKT ? maxKT : (int)kTick;
00215   
00216   unsigned int *outPlace=myTrigData->TrgSum.L2Result+ mResultOffset;
00217   memcpy(outPlace,&out,sizeof( L2pedResults2006));
00218   
00219   if(par_dbg) L2pedResults2006_print(&out);
00220   
00221   return  mAccept;
00222 }
00223 
00224 /*========================================
00225   ======================================== */
00226 void 
00227 L2pedAlgo::finishRun() {/* called once at the end of the run */
00228   if(run_number<=0) return; // algo not used for current run
00229   int nBtowLow=0, nBtowHigh=0 ,nEtowLow=0, nEtowHigh=0 ;
00230 
00231   char fname[1000];
00232   sprintf(fname,"%s/run%d.l2ped.out",mOutDir,run_number); // actual L2 output destination 
00233   printf("L2ped_finish('%s') , finding pedestals...\n",fname);
00234   
00235   FILE *fd=fopen(fname,"w");
00236   if(fd==0) {printf("failed to open output %s file,skip ped_finish()\n",fname); return;}  
00237   fprintf(fd,"#L2-ped algorithm finishRun(%d), compiled: %s , %s\n",run_number,__DATE__,__TIME__);
00238   fprintf(fd,"#params: pedSubtr=%d speedFact=%d saveBin=%d debug=%d prescAccept=%d\n",par_pedSubtr,par_speedFact,par_saveBinary,par_dbg,par_prescAccept);
00239   hA[10]->printCSV(fd); // event accumulated
00240   int iMax=-3, iFWHM=-4;
00241   hA[11]->findMax( &iMax, &iFWHM);
00242   fprintf(fd,"#L2ped  CPU/eve MPV %d kTicks,  FWHM=%d, seen eve=%d\n",iMax, iFWHM,nInp);
00243   printf("L2ped  CPU/eve MPV %d kTicks,  FWHM=%d, seen eve=%d\n",iMax, iFWHM,nInp);
00244   if(par_saveBinary)  fprintf(fd,"#L2ped  will save full spectra for all towers\n");
00245 
00246   int par_topAdc=100;
00247   int maxPedDeviation=5;
00248 //int iadcHigh=maxPedDeviation-minAdc;
00249 // int iadcLow =maxPedDeviation+minAdc;
00250 
00251 // fixed sign on July 23, 2007, JanB
00252   int iadcHigh=maxAdc - maxPedDeviation;
00253   int iadcLow =minAdc + maxPedDeviation;
00254 
00255   char xAxis[100];
00256   sprintf(xAxis,"raw ADC + %d",-minAdc);
00257   if(par_pedSubtr) sprintf(xAxis,"ADC - ped + %d",-minAdc);
00258   
00259   
00260   fprintf(fd,"# L2ped-Adc spectra, run=%d, Z-scale is ln(yield), only first digit shown;  maxPedDev=%d  table format:\n# name, ped, sigPed, crate, chan, softID-m-s-e, RDO_ID;\n#                                   ADC spectrum: [%d ...  <=-10 ... *=0  ...  >=+10 ... :=+20 ... %d],  Xaxis=%s\n",run_number,maxPedDeviation,minAdc,par_topAdc,xAxis);
00261   
00262   int i;
00263   int nB=0;
00264   for(i=0; i<EmcDbIndexMax; i++) {
00265     const L2EmcDb::EmcCDbItem *x=mDb->getByIndex(i);
00266     if(mDb->isEmpty(x)) continue;  /* dropped not mapped  channels */
00267     if (mDb->isBTOW(x) ||mDb->isETOW(x) ) { //
00268       // mDb->printItem(x);
00269       nB++;
00270       /* if(nB>30) return; */
00271       int iMax=-3, iFWHM=-4; 
00272       char pedQA='?';
00273       L2Histo *h=0;
00274       if(mDb->isBTOW(x)) h= btowAdc[x->rdo];
00275       else if(mDb->isETOW(x)) h= etowAdc[x->rdo];
00276       else continue;
00277 
00278 
00279       int pedRes=999;
00280       int maxRes=3;
00281       if(h->findMax( &iMax, &iFWHM)) {
00282         pedQA='0';
00283         if(iMax<iadcLow) pedQA='-';
00284         else if(iMax>iadcHigh) pedQA='+';
00285         pedRes=iMax+minAdc;
00286         if(!par_pedSubtr) pedRes=int(pedRes - x->ped); // this looks funny but is right
00287       }
00288 
00289       if(mDb->isBTOW(x)) { //BTOW ...............
00290         if(pedQA=='-' )  nBtowLow++;
00291         else if(pedQA=='+')  nBtowHigh++;
00292         //........ residual monito histos ....
00293         int ieta= (x->eta-1);
00294         int iphi= (x->sec-1)*10 + x->sub-'a' ;
00295         int ihard=x->chan+(x->crate-1)*160;
00296         if(x->fail) pedRes =0;
00297         hA[20]->fillW(ihard,pedRes);
00298         if(x->fail) {
00299           pedRes =-100;
00300         } else {
00301           if(pedRes<-maxRes) pedRes=-maxRes; 
00302           if(pedRes>maxRes) pedRes=maxRes;
00303         }       
00304         hA[21]->fillW(ieta, iphi,pedRes);
00305      } else {             // ETOW  ...........
00306         if(pedQA=='-' )  nEtowLow++;
00307         else if(pedQA=='+')  nEtowHigh++;  
00308         //........ residual monito histos ....
00309         int ieta= 12-x->eta;
00310         int iphi= (x->sec-1)*5 + x->sub-'A' ;
00311         int ihard=x->chan+(x->crate-1)*128;
00312         if(x->fail) pedRes =0;
00313         hA[30]->fillW(ihard,pedRes);
00314         if(x->fail) {
00315           pedRes =-100;
00316         } else {
00317           if(pedRes<-maxRes) pedRes=-maxRes; 
00318           if(pedRes>maxRes) pedRes=maxRes;
00319         }       
00320         hA[31]->fillW(ieta, iphi,pedRes);
00321       }
00322    
00323 
00324       char okC=' ';
00325       if(x->fail) okC='#';
00326       fprintf(fd,"%c%s %3d %4.1f 0x%02x 0x%02x %15s %4d ",okC,x->name,iMax+minAdc,iFWHM/2.3,x->crate, x->chan,x->tube, x->rdo);
00327       h->printPed(fd,minAdc,par_topAdc,' ');
00328       fprintf(fd,"qa=%c\n",pedQA);
00329       // exit(1); 
00330     } /* end of BTOW & ETOW */
00331   } /* end of pixels */
00332 
00333    
00334   fprintf(fd,"#L2ped_finishRun() # of towers with |ped-pedDB| >10 chan\n#    BTOW: nLow=%d nHigh=%d ;   ETOW  nLow=%d, nHigh=%d\n",nBtowLow, nBtowHigh,nEtowLow, nEtowHigh);
00335   fprintf(fd,"#    found peds for nB+E=%d , seen events=%d\n",nB,nInp);
00336   printf("l2ped_finish() found peds for nB+E=%d , seen events=%d\n",nB,nInp);
00337 
00338   fclose(fd);
00339 
00340 
00341   sprintf(fname,"%s/run%d.l2ped.hist.bin",mOutDir,run_number); // full spectra in binary form
00342   
00343   fd=fopen(fname,"w");
00344   if(fd==0) {
00345     printf("failed to open output %s file,skip ped_finish()\n",fname); 
00346     goto end;
00347   }  
00348   
00349   for(int j=0;j<mxHA;j++) { // save auxil histos
00350     if(hA[j]==0) continue;
00351     hA[j]->write(fd);
00352   }
00353   
00354   if(par_saveBinary) {
00355     printf("l2ped_finish('%s') , save FULL spectra binary...\n",fname);
00356 
00357     // .................. SAVE FULL SPECTRA BINARY ...........
00358     for(i=0;i<mxHA;i++) if(hA[i]) hA[i]->write(fd);
00359     for(i=0; i<EmcDbIndexMax; i++) {
00360       const L2EmcDb::EmcCDbItem *x=mDb->getByIndex(i);
00361       if(mDb->isEmpty(x)) continue;  /* dropped not mapped  channels */
00362       L2Histo *h=0;
00363       char tit[200];
00364       if(mDb->isBTOW(x) ){
00365         h= btowAdc[x->rdo];
00366         //L2EmcDb::printItem(x);
00367         sprintf(tit,"BTOW=%s cr/ch=%03d/%03d stat/0x=%04x+%04x soft=%s; %s",x->name, x->crate, x->chan, x->stat, x->fail,x->tube,xAxis);
00368         //printf("tit=%s=\n",tit);      
00369       } else if(mDb->isETOW(x) ) {
00370         h= etowAdc[x->rdo];
00371         //L2EmcDb::printItem(x);
00372         sprintf(tit,"ETOW=%s cr/ch=%03d/%03d stat/0x=%04x+%04x pname=%s; %s",x->name, x->crate, x->chan, x->stat, x->fail,x->tube,xAxis);
00373       }
00374       if(h==0) continue; //just in case
00375       h->setTitle(tit);
00376       h->write(fd); // change title, add cr/chan/name, pedSubtrFlag
00377       // break;
00378     }
00379     printf("l2ped_finish() binary full spectra saved\n");
00380   }
00381   fclose(fd);
00382  end:
00383   run_number=-2;
00384 }
00385 
00386          
00387 
00388 /**********************************************************************
00389   $Log: L2pedAlgo.cxx,v $
00390   Revision 1.11  2009/11/19 15:48:46  balewski
00391   add (char*) to many strings to make SL5 happ, few other adjustments
00392 
00393   Revision 1.10  2007/11/19 22:18:31  balewski
00394   most L2algos provide triggerID's
00395 
00396   Revision 1.9  2007/11/18 21:58:58  balewski
00397   L2algos triggerId list fixed
00398 
00399   Revision 1.8  2007/11/13 23:06:09  balewski
00400   toward more unified L2-algos
00401 
00402   Revision 1.7  2007/11/13 00:12:38  balewski
00403   added offline triggerID, take1
00404 
00405   Revision 1.6  2007/11/08 04:02:33  balewski
00406   run on l2ana as well
00407 
00408   Revision 1.5  2007/11/02 20:44:54  balewski
00409   cleanup
00410 
00411   Revision 1.4  2007/11/02 17:43:11  balewski
00412   cleanup & it started to work w/ L2upsilon
00413 
00414   Revision 1.3  2007/11/02 03:03:50  balewski
00415   modified L2VirtualAlgo
00416 
00417   Revision 1.2  2007/10/25 02:07:06  balewski
00418   added L2upsilon & binary event dump
00419 
00420   Revision 1.1  2007/10/11 00:33:24  balewski
00421   L2algo added
00422 
00423   Revision 1.5  2006/03/28 19:46:51  balewski
00424   ver16b, in l2new
00425 
00426   Revision 1.4  2006/03/11 17:08:35  balewski
00427   now CVS comments should work
00428 
00429 */
00430