20 #ifdef IS_REAL_L2 //in l2-ana environment
21 #include "../L2algoUtil/L2EmcDb.h"
22 #include "../L2algoUtil/L2Histo.h"
24 #include "StTriggerUtilities/L2Emulator/L2algoUtil/L2EmcDb.h"
25 #include "StTriggerUtilities/L2Emulator/L2algoUtil/L2Histo.h"
28 #include "L2adc2energyAlgo.h"
29 #include "Map_DeltaPhiJets.h"
33 L2adc2energyAlgo::L2adc2energyAlgo(
const char* name,
L2EmcDb* db,
char* outDir,
int resOff)
40 par_adcMask= (
unsigned short) (-0x10);
44 printf(
"L2adc2energyAlgo instantiated, logPath='%s'\n",mOutDir);
50 L2adc2energyAlgo::paramsChanged(
int *rc_ints,
float *rc_floats) {
53 if(rc_ints[i]!=raw_ints[i])
goto foundProblem;
56 if(fabs(rc_floats[i]-raw_floats[i])>0.00001)
goto foundProblem;
61 if (mLogFile) fprintf(mLogFile,
"L2jet-ghost initRun - inconsistent params, ABORT initialization\n");
68 L2adc2energyAlgo::initRun(
int runNo,
int *rc_ints,
float *rc_floats) {
71 if(mDb->initRun(runNo))
return -7;
74 if(run_number==runNo) {
75 if (mLogFile) fprintf(mLogFile,
"L2jet::initRun-%s(%d)=ghost already initilized, only check params\n",mName, runNo);
76 printf(
"L2jet::initRun-%s(%d)=ghost already initilized, only checking params\n",mName, runNo);
77 int ret= paramsChanged(rc_ints, rc_floats);
82 fprintf(mLogFile,
"L2jet algorithm init: crashA in internal logic\n");
91 for (i=0; i<mxHA;i++)
if(hA[i])hA[i]->reset();
94 memset(db_btowThr, 0xFFFF,
sizeof(db_btowThr));
95 memset(db_btowPedS, 0,
sizeof(db_btowPedS));
96 memset(db_btowGainCorr, 0,
sizeof(db_btowGainCorr));
97 memset(db_btowL2PhiBin, 0,
sizeof(db_btowL2PhiBin));
98 memset(db_btowL2PatchBin,0,
sizeof(db_btowL2PatchBin));
100 memset(db_etowThr, 0xFFFF,
sizeof(db_etowThr));
101 memset(db_etowPedS, 0 ,
sizeof(db_etowPedS));
102 memset(db_etowGainCorr, 0 ,
sizeof(db_etowGainCorr));
103 memset(db_etowL2PhiBin, 0 ,
sizeof(db_etowL2PhiBin));
104 memset(db_etowL2PatchBin,0 ,
sizeof(db_etowL2PatchBin));
112 int par_IgainCorrOne=30;
113 int par_IgainCorrMin=5;
114 int par_IgainCorrMax=60;
116 run_startUnix=time(0);
119 raw_floats =rc_floats;
120 run_nEventOneJet=run_nEventDiJet= run_nEventRnd=0;
123 sprintf(Fname,
"%s/run%d.l2jet.out",mOutDir,run_number);
124 printf(
"L2jet::initRun-%s('%s') ...\n",mName,Fname);
127 mLogFile = fopen(Fname,
"w");
128 if( mLogFile==0) printf(
" L2adc2energyAlgo() UNABLE to open run summary log file, continue anyhow\n");
132 par_cutTag = rc_ints[0];
133 par_useBtowEast= (rc_ints[1]&1 )>0;
134 par_useBtowWest= (rc_ints[1]&2)>0;
135 par_useEndcap = rc_ints[2]&1;
136 int par_adcThr = rc_ints[3];
137 par_minPhiBinDiff=rc_ints[4];
139 par_oneJetThr = rc_floats[0];
140 par_diJetThrHigh= rc_floats[1];
141 par_diJetThrLow = rc_floats[2];
142 par_rndAccProb = rc_floats[3];
143 par_dbg =(int)rc_floats[4];
146 par_energyScale=par_maxADC*par_IgainCorrOne/par_maxEt;
148 float monTwEtThr=2.0;
149 par_hotTwEtThr= (int)(monTwEtThr*par_energyScale);
151 par_rndAccThr= int(par_rndAccProb* RAND_MAX);
152 if(par_rndAccProb<0) {
155 }
else if (par_rndAccProb>0.9999) {
156 par_rndAccThr=RAND_MAX;
160 fprintf(mLogFile,
"L2jet algorithm initRun(%d), compiled: %s , %s\n params:\n",run_number,__DATE__,__TIME__);
161 fprintf(mLogFile,
" - use BTOW: East=%d West=%d, Endcap=%d L2ResOffset=%d\n", par_useBtowEast, par_useBtowWest,par_useEndcap ,mResultOffset);
162 fprintf(mLogFile,
" - threshold: ADC-ped> %d \n", par_adcThr);
163 fprintf(mLogFile,
" - min phi opening angle Jet1<->Jet2: %d in L2phiBins\n",par_minPhiBinDiff);
164 fprintf(mLogFile,
" - diJet Et thrHigh= %.2f (GeV) thrLow= %.2f (GeV)\n", par_diJetThrHigh, par_diJetThrLow);
165 fprintf(mLogFile,
" - oneJet Et thr = %.2f (GeV) ; rndAccProb=%f; cutTag=%d \n",par_oneJetThr,par_rndAccProb,par_cutTag);
166 fprintf(mLogFile,
" - debug=%d, hot tower threshold: Et> %.1f GeV ( only monitoring)\n",par_dbg, monTwEtThr);
171 kBad+=0x0001 * ( !par_useBtowEast & !par_useBtowWest & !par_useEndcap);
172 kBad+=0x0002 * (par_adcThr<par_pedOff);
173 kBad+=0x0004 * (par_adcThr>16);
174 kBad+=0x0008 * (par_minPhiBinDiff<5);
175 kBad+=0x0010 * (par_minPhiBinDiff>=15);
176 kBad+=0x0020 * (par_oneJetThr<3.);
177 kBad+=0x0040 * (par_oneJetThr>12.);
178 kBad+=0x0080 * (par_diJetThrLow<2.9);
179 kBad+=0x0100 * (par_diJetThrHigh<par_diJetThrLow);
180 kBad+=0x0200 * (par_diJetThrHigh>12.);
181 kBad+=0x0400 * (par_cutTag<=0 || par_cutTag>255);
182 kBad+=0x0800 * (par_rndAccProb<0. || par_rndAccProb>1.);
184 fprintf(mLogFile,
"L2jet initRun() params checked for consistency, Error flag=0x%04x\n",kBad);
185 if(kBad) fprintf(mLogFile,
"L2jet initRun() ABORT\n");
191 fprintf(mLogFile,
"L2jet algorithm init: crashB in internal logic\n");
198 sprintf(tit,
"# BTOW towers>ped+%d (input); x: # of towers/event",par_adcThr);
199 hA[47]->setTitle(tit);
201 sprintf(tit,
"# ETOW towers>ped+%d (input); x: # of towers/event",par_adcThr);
202 hA[48]->setTitle(tit);
207 const float edgeEtaBinEtow[] = {
209 1.9008 , 1.8065 , 1.7168 , 1.6317 , 1.5507 , 1.4738 ,
210 1.4007 , 1.3312 , 1.2651 , 1.2023 , 1.1427 , 1.086 ,
214 const int mxEtaBinsE=12,mxEtaBinsB=40;
215 float idealGainEtow[mxEtaBinsE], idealGainBtow[mxEtaBinsB];
216 float coshEtow[mxEtaBinsE],coshBtow[mxEtaBinsB];
218 for(i=0;i<mxEtaBinsE;i++ ){
219 float avrEta=(edgeEtaBinEtow[i]+edgeEtaBinEtow[i+1])/2.;
220 coshEtow[i]=cosh(avrEta);
221 idealGainEtow[i]=par_maxADC/par_maxEt/coshEtow[i];
226 for(i=0;i<mxEtaBinsB;i++ ){
227 float avrEta=-0.975 +i*0.05;
228 coshBtow[i]=cosh(avrEta);
229 idealGainBtow[i]=par_maxADC/par_maxEt/coshBtow[i];
236 int etowEtaBin2Patch[mxEtaBinsE]={14,14,13,13,12,12,11,11,11,10,10,10};
241 for(i=0; i<EmcDbIndexMax; i++) {
243 if(mDb->isEmpty(x))
continue;
244 if(x->fail)
continue;
245 if(x->gain<=0)
continue;
250 if (mDb->isBTOW(x) ) {
253 if(x->eta<0 || x->eta>mxEtaBinsB)
goto crashIt_1;
254 if(!par_useBtowEast && x->eta<=20)
continue;
255 if(!par_useBtowWest && x->eta>=21)
continue;
257 int iphiTw=(x->sec-1)*10 + x->sub-
'a';
260 if(iphiTw<0) iphiTw=119;
265 int IgainCor=int(par_IgainCorrOne*idealGainBtow[x->eta-1]/x->gain);
268 if(IgainCor <par_IgainCorrMin)
continue;
269 if(IgainCor >par_IgainCorrMax)
continue;
270 db_btowGainCorr[x->rdo]=IgainCor;
272 db_btowL2PhiBin[x->rdo]=iphiP;
273 db_btowL2PatchBin[x->rdo]=ietaP+ iphiP*cl2jetMaxEtaBins;
274 db_btowThr[x->rdo]=(int) (x->ped+par_adcThr);
275 db_btowPedS[x->rdo]=(
unsigned short) (par_pedOff-x->ped);
277 }
else if(mDb->isETOW(x) && par_useEndcap) {
280 int iphiTw= (x->sec-1)*5 + x->sub-
'A';
283 if(iphiTw<0) iphiTw=59;
286 if(x->eta<0 || x->eta>mxEtaBinsE)
goto crashIt_1;
287 ietaP=etowEtaBin2Patch[x->eta-1];
289 int IgainCor=int(par_IgainCorrOne*idealGainEtow[x->eta-1]/x->gain);
292 if(IgainCor <par_IgainCorrMin)
continue;
293 if(IgainCor >par_IgainCorrMax)
continue;
294 db_etowGainCorr[x->rdo]=IgainCor;
296 db_etowL2PhiBin[x->rdo]=iphiP;
297 db_etowL2PatchBin[x->rdo]=ietaP+ iphiP*cl2jetMaxEtaBins;
298 db_etowThr[x->rdo]=(int) (x->ped+par_adcThr);
299 db_etowPedS[x->rdo]=(
unsigned short) (par_pedOff-x->ped);
306 fprintf(mLogFile,
"L2jet algorithm: found working/calibrated: %d/%d=ETOW & %d/%d=BTOW, based on ASCII DB\n",nE,nEg,nB,nBg);
314 fprintf(mLogFile,
"L2jet algorithm init: crashC in internal logic\n");
325 L2adc2energyAlgo::doEvent(
int L0trg,
int inpEveId,
TrgDataType* trgData,
326 int bemcIn,
unsigned short *bemcData,
327 int eemcIn,
unsigned short *eemcData){
329 rdtscl_macro(mEveTimeStart);
331 if(L0trg==1) hA[10]->fill(1);
332 else if(L0trg==2) hA[10]->fill(2);
334 if(eve_ID!=inpEveId) {
345 int runTimeSec=time(0)- run_startUnix;
347 hA[12]->fill(runTimeSec);
349 if(par_dbg>1) printf(
"\n......... in L2adc2E_doEvent(ID=%d)... bIn=%d eIn=%d\n",eve_ID,bemcIn,eemcIn);
351 if (bemcIn || eemcIn){
356 int nBtowTw=0, nEtowTw=0;
358 if(bemcIn==1 && (par_useBtowEast||par_useBtowWest) ) {
359 nBtowTw=projectAdc( bemcData, MaxBtowRdo,
360 db_btowThr, db_btowPedS, db_btowGainCorr,
361 db_btowL2PhiBin, db_btowL2PatchBin,
366 if(eemcIn==1 && par_useEndcap ) {
367 nEtowTw=projectAdc( eemcData, 720,
368 db_etowThr, db_etowPedS, db_etowGainCorr,
369 db_etowL2PhiBin, db_etowL2PatchBin,
376 if( mAccept) hA[10]->fill(8);
382 rdtscl_macro(mEveTimeStop);
383 mEveTimeDiff=mEveTimeStop-mEveTimeStart;
384 int kTick=mEveTimeDiff/1000;
404 L2adc2energyAlgo::finishRun() {
405 if(run_number<0)
return;
408 sprintf(Fname,
"%s/run%d.l2adc2E.hist.bin",mOutDir,run_number);
409 printf(
"L2jet::finishRun('%s') , save histo ...\n",Fname);
410 mHistFile = fopen(Fname,
"w");
413 fprintf(mLogFile,
"L2-jet algorithm finishRun(%d)\n",run_number);
414 fprintf(mLogFile,
" - %d events seen by L2 di-jet\n",mEventsInRun);
415 fprintf(mLogFile,
" - accepted: rnd=%d oneJet=%d diJet=%d \n", run_nEventRnd, run_nEventOneJet, run_nEventDiJet);
419 hA[10]->printCSV(mLogFile);
425 printf(
" L2adc2energyAlgo: finishRun() UNABLE to open run summary log file, continue anyhow\n");
427 fprintf(mLogFile,
"L2 di-jet histos NOT saved, I/O error\n");
431 for(j=0;j<mxHA;j++) {
432 if(hA[j]==0)
continue;
433 hA[j]->write(mHistFile);
436 finishCommonHistos();
440 fprintf(mLogFile,
"L2 di-jet: %d histos saved to '%s'\n",nh,Fname);
446 if (mLogFile && mLogFile!=stdout) {
457 L2adc2energyAlgo::createHisto() {
458 memset(hA,0,
sizeof(hA));
460 hA[10]=
new L2Histo(10, (
char*)
"total event counter; x=cases",9);
461 hA[11]=
new L2Histo(11, (
char*)
"L2 time used per input event; x: time (CPU kTics), range=100muSec; y: events ",160);
463 int mxRunDration=2500;
464 hA[12]=
new L2Histo(12, (
char*)
"rate of input events; x: time in this run (seconds); y: rate (Hz)", mxRunDration);
466 hA[13]=
new L2Histo(13, (
char*)
"rate of accepted one-Jet; x: time in this run (seconds); y: rate (Hz)", mxRunDration);
467 hA[14]=
new L2Histo(14, (
char*)
"rate of accepted di-Jet ; x: time in this run (seconds); y: rate (Hz)", mxRunDration);
468 hA[15]=
new L2Histo(15, (
char*)
"rate of random accepted ; x: time in this run (seconds); y: rate (Hz)", mxRunDration);
471 hA[20]=
new L2Histo(20, (
char*)
"BTOW tower, Et>2.0 GeV (input); x: BTOW RDO index=chan*30+fiber; y: counts", 4800);
472 hA[21]=
new L2Histo(21, (
char*)
"BTOW tower, Et>2.0 GeV (input); x: BTOW softID", 4800);
473 hA[22]=
new L2Histo(22, (
char*)
"BTOW tower, Et>2.0 GeV (input); x: eta bin, [-1,+1]; y: phi bin ~sector",40,120);
476 hA[30]=
new L2Histo(30, (
char*)
"ETOW tower, Et>2.0 GeV (input); x: ETOW RDO index=chan*6+fiber; y: counts", 720 );
477 hA[31]=
new L2Histo(31, (
char*)
"ETOW tower, Et>2.0 GeV (input); x: i=chan+128*crate", 768);
478 hA[32]=
new L2Histo(32, (
char*)
"ETOW tower, Et>2.0 GeV (input); x: 12 - Endcap etaBin ,[+1,+2]; y: phi bin ~sector",12,60);
481 hA[40]=
new L2Histo(40, (
char*)
"Et Jet1-Jet2 (input); x: Jet1 Et/GeV ; Jet2 Et/GeV",12,12);
482 hA[41]=
new L2Histo(41, (
char*)
"diJet1 eta-phi (input); x: iEta [-1,+2] ; y: iPhi ~sector ",15,30);
483 hA[42]=
new L2Histo(42, (
char*)
"diJet2 eta-phi (input); x: iEta [-1,+2] ; y: iPhi ~sector",15,30);
485 hA[43]=
new L2Histo(43, (
char*)
"diJet phi1-phi2 (input); x: iPhi1 ~sector ; y: iPhi2 ~sector ",30,30);
487 hA[44]=
new L2Histo(44, (
char*)
"Jet1 Et (input); x: Et (GeV)", 60);
488 hA[45]=
new L2Histo(45, (
char*)
"Jet2 Et (input); x: Et (GeV)", 60);
489 hA[46]=
new L2Histo(46, (
char*)
"total Et (input); x: Et (GeV)", 60);
490 hA[47]=
new L2Histo(47, (
char*)
"# BTOW towers>thrXX (input); x: # of towers/event", 200);
491 hA[48]=
new L2Histo(48, (
char*)
"# ETOW towers>thrXX (input); x: # of towers/event", 100);
494 hA[50]=
new L2Histo(50, (
char*)
"one-Jet Et (accepted); x: jet Et (GeV)", 60);
495 hA[51]=
new L2Histo(51, (
char*)
"one-Jet eta-phi (accepted); x: iEta [-1,+2] ; y: iPhi ~sector ",15,30);
496 hA[52]=
new L2Histo(52, (
char*)
"one-Jet eta (accepted); x: iEta [-1,+2]", 15);
497 hA[53]=
new L2Histo(53, (
char*)
"one-Jet phi (accepted); x: iPhi ~sector", 30);
500 hA[60]=
new L2Histo(60, (
char*)
"Et of Jet1 vs. Jet2 (accepted); x: Jet1/GeV ; Jet2/GeV",12,12);
501 hA[61]=
new L2Histo(61, (
char*)
"diJet1 eta-phi (accepted); x: iEta [-1,+2] ; y: iPhi ~sector ",15,30);
502 hA[62]=
new L2Histo(62, (
char*)
"diJet2 eta-phi (accepted); x: iEta [-1,+2] ; y: iPhi ~sector",15,30);
504 hA[63]=
new L2Histo(63, (
char*)
"diJet phi1-phi2 (accepted); x: iPhi1 ~sector ; y: iPhi2 ~sector ",30,30);
506 hA[64]=
new L2Histo(64, (
char*)
"diJet1 Et (accepted); x: Et (GeV)", 60);
507 hA[65]=
new L2Histo(65, (
char*)
"diJet2 Et (accepted); x: Et (GeV)", 60);
509 hA[66]=
new L2Histo(66, (
char*)
"diJet1 eta (accepted); x: i Eta [-1,+2]", 15);
510 hA[67]=
new L2Histo(67, (
char*)
"diJet2 eta (accepted); x: i Eta [-1,+2]", 15);
511 hA[68]=
new L2Histo(68, (
char*)
"diJet1 phi (accepted); x: iPhi ~sector", 30);
512 hA[69]=
new L2Histo(69, (
char*)
"diJet2 phi (accepted); x: iPhi ~sector", 30);
513 hA[70]=
new L2Histo(70, (
char*)
"diJet delZeta (accepted); x: delta zeta (rad*10)", MxPhiRad10);
514 hA[71]=
new L2Histo(71, (
char*)
"diJet delZeta vs. eta1 (accepted); x: iEta1 [-1,+2] ; y: delta zeta (rad*10)",15, MxPhiRad10);
515 hA[72]=
new L2Histo(72, (
char*)
"diJet eta2 vs. eta1 (accepted); x: iEta1 [-1,+2] ;x: iEta2 [-1,+2] ",15,15);
516 hA[73]=
new L2Histo(73, (
char*)
"diJet delZeta vs. avrPhi (accepted); x: (iphi1+iphi2)/2 (12 deg/bin); y: delta zeta (rad*10)",30, MxPhiRad10);
517 hA[74]=
new L2Histo(74, (
char*)
"total Et diJet (accepted); x: Et (GeV)", 60);
524 L2adc2energyAlgo::clearEvent(){
529 memset(eve_patchEne,0,
sizeof(eve_patchEne));
530 memset(eve_phiEne,0,
sizeof(eve_phiEne));
538 L2adc2energyAlgo::projectAdc(
unsigned short *rawAdc,
int nRdo,
539 unsigned short *thr,
unsigned short *pedS,
unsigned short *gainCorr,
540 unsigned short *phiBin,
unsigned short *patchBin,
547 for(rdo=0; rdo<nRdo; rdo++){
548 if(rawAdc[rdo]<thr[rdo])
continue;
550 adc=(rawAdc[rdo]+pedS[rdo]) & par_adcMask ;
551 adc4=adc*gainCorr[rdo];
552 eve_patchEne[patchBin[rdo]]+=adc4;
553 eve_phiEne[phiBin[rdo]]+=adc4;
557 if(adc4 >par_hotTwEtThr) hHot->fill(rdo);
567 L2adc2energyAlgo:: dumpPatchEneA(){
570 for(iy=0;iy<cl2jetMaxPhiBins;iy++) {
571 int *patchEneA=eve_patchEne+(iy*cl2jetMaxEtaBins);
573 for(ix=0;ix<cl2jetMaxEtaBins;ix++,patchEneA++){
574 printf(
" %6d",*patchEneA);
576 printf(
" iPhi=%d\n",iy);
585 L2adc2energyAlgo::finishRunHisto(){
588 const int *data20=hA[20]->getData();
589 const int *data30=hA[30]->getData();
591 int bHotSum=1,bHotId=-1;
596 for(i=0; i<EmcDbIndexMax; i++) {
598 if(mDb->isEmpty(x))
continue;
599 if (mDb->isBTOW(x) ) {
600 int softId=atoi(x->tube+2);
601 int ieta= (x->eta-1);
602 int iphi= (x->sec-1)*10 + x->sub-
'a' ;
604 hA[21]->fillW(softId,data20[x->rdo]);
605 hA[22]->fillW(ieta, iphi,data20[x->rdo]);
606 if(bHotSum<data20[x->rdo]) {
607 bHotSum=data20[x->rdo];
612 else if (mDb->isETOW(x) ) {
613 int ihard=x->chan+(x->crate-1)*128;
615 int iphi= (x->sec-1)*5 + x->sub-
'A' ;
616 hA[31]->fillW(ihard,data30[x->rdo]);
617 hA[32]->fillW(ieta, iphi,data30[x->rdo]);
618 if(eHotSum<data30[x->rdo]) {
619 eHotSum=data30[x->rdo];
626 fprintf(mLogFile,
"L2jet::finishRun()\n");
627 fprintf(mLogFile,
"#BTOW_hot tower _candidate_ (bHotSum=%d) :, softID %d , crate %d , chan %d , name %s\n",bHotSum,bHotId,xB->crate,xB->chan,xB->name);
628 fprintf(mLogFile,
"#ETOW_hot tower _candidate_ (eHotSum=%d) :, name %s , crate %d , chan %d\n",eHotSum,xE->name,xE->crate,xE->chan);