16 #ifdef IS_REAL_L2 //in l2-ana environment
17 #include "../L2algoUtil/L2EmcDb.h"
18 #include "../L2algoUtil/L2Histo.h"
19 #else //full path needed for cvs'd code
20 #include "StTriggerUtilities/L2Emulator/L2algoUtil/L2EmcDb.h"
21 #include "StTriggerUtilities/L2Emulator/L2algoUtil/L2Histo.h"
22 #include "StTriggerUtilities/L2Emulator/L2algoUtil/L2EmcGeom.h"
25 #include "L2eemcGamma2009.h"
30 void L2eemcGamma2009::quickSort(
int array[],
int start,
int end)
37 float pivot = wrkEtow_et[array[start]];
41 while (wrkEtow_et[array[i]] <= pivot && i <= end && k > i)
43 while (wrkEtow_et[array[k]] > pivot && k >= start && k >= i)
48 swap(array, start, k);
50 quickSort(array, start, k - 1);
51 quickSort(array, k + 1, end);
59 void L2eemcGamma2009::swap(
int array[],
int index1,
int index2)
63 int temp = array[index1];
64 array[index1] = array[index2];
72 L2eemcGamma2009::L2eemcGamma2009(
const char* name,
L2EmcDb* db,
L2EmcGeom *geoX,
char* outDir,
int resOff) :
L2VirtualAlgo2009( name, db, outDir, false,true, resOff) {
78 criticalError(
"L2eemcGamma is broken -- can't find geom.");
86 criticalError(
"L2eemcGamma has failed consistency check. sizeof(L2gammaResult2009)!= L2gammaResult2009::mySizeChar");
95 L2eemcGamma2009::initRunUser(
int runNo,
int *rc_ints,
float *rc_floats) {
99 par_RndAcceptPrescale = rc_ints[1];
100 par_seedEtThres = rc_floats[0];
101 par_clusterEtThres= rc_floats[1];
102 par_eventEtThres = rc_floats[2];
106 kBad+=0x00001 * (par_seedEtThres<1.0);
107 kBad+=0x00002 * (par_clusterEtThres<par_seedEtThres);
110 fprintf(mLogFile,
"L2%s algorithm initRun(R=%d), compiled: %s , %s\n params:\n",getName(),mRunNumber,__DATE__,__TIME__);
111 fprintf(mLogFile,
" - use Thres/GeV seed=%.2f, clusterList=%.2f debug=%d, prescale=%d\n", par_seedEtThres,par_clusterEtThres, par_dbg,par_RndAcceptPrescale);
113 fprintf(mLogFile,
" - accept event cluster Thres/GeV=%.2f\n",par_eventEtThres);
114 fprintf(mLogFile,
"initRun() params checked for consistency, Error flag=0x%04x\n",kBad);
120 for (i=0; i<mxHA;i++)
if(hA[i])hA[i]->reset();
122 if(kBad)
return -1*kBad;
124 memset(mEtow,0,
sizeof(mEtow));
128 sprintf(txt,
"ETOW-compute: #seed towers ET>%.2f GeV / event; x: # ETOW towers; y: counts",par_seedEtThres);
129 hA[3]->setTitle(txt);
131 sprintf(txt,
"ETOW-decision: #cluster ET>%.2f GeV / event ; x: # ETOW towers; y: counts",par_clusterEtThres);
132 hA[4]->setTitle(txt);
134 sprintf(txt,
"ETOW-decision: acc cluster ET>%.2f GeV; x:cluster ET(GeV) ; y: counts",par_eventEtThres);
135 hA[4]->setTitle(txt);
139 for (
int index=0; index<EmcDbIndexMax; index++ )
142 if ( x==0 )
continue;
143 if ( !mDb->isETOW(x) )
continue;
144 int sec = x->sec - 1;
147 int eta = x->eta - 1;
148 int phi = EtowGeom::mxSubs *sec + sub;
149 int tow = EtowGeom::mxEtaBin *phi + eta;
151 if (tow<0 || tow>mxEtow || rdo<0 || rdo>mxEtow)
return -101;
153 mTower2rdo[ tow ] = rdo;
154 mRdo2tower[ rdo ] = tow;
163 L2eemcGamma2009::countNonZeroTow(
int phi,
int eta) {
164 int tow = EtowGeom::mxEtaBin *((phi+EtowGeom::mxPhiBin)%EtowGeom::mxPhiBin) + ((eta+EtowGeom::mxEtaBin)%EtowGeom::mxEtaBin);
165 const int maxTowers = EtowGeom::mxEtaBin * EtowGeom::mxPhiBin;
169 if (etow_used[tow]==0) {
170 if(wrkEtow_et[tow]!=0) count++;
174 towPlusOne%= maxTowers;
175 if (etow_used[towPlusOne]==0) {
176 if(wrkEtow_et[towPlusOne]!=0) count++;
179 tow+=EtowGeom::mxEtaBin;
181 if (etow_used[tow]==0) {
182 if(wrkEtow_et[tow]!=0) count++;
186 towPlusOne%= maxTowers;
187 if (etow_used[towPlusOne]==0) {
188 if(wrkEtow_et[towPlusOne]!=0) count++;
197 L2eemcGamma2009::flagUsed(
int phi,
int eta) {
198 int tow = EtowGeom::mxEtaBin *((phi+EtowGeom::mxPhiBin)%EtowGeom::mxPhiBin) + ((eta+EtowGeom::mxEtaBin)%EtowGeom::mxEtaBin);
199 const int maxTowers = EtowGeom::mxEtaBin * EtowGeom::mxPhiBin;
204 towPlusOne%= maxTowers;
205 etow_used[towPlusOne] = 1;
207 tow+=EtowGeom::mxEtaBin;
212 towPlusOne%= maxTowers;
213 etow_used[towPlusOne] = 1;
219 L2eemcGamma2009::averageEtaPhi(
int phi,
int eta,
float *avePhi,
float *aveEta) {
220 int tow = EtowGeom::mxEtaBin *((phi+EtowGeom::mxPhiBin)%EtowGeom::mxPhiBin) + ((eta+EtowGeom::mxEtaBin)%EtowGeom::mxEtaBin);
221 const int maxTowers = EtowGeom::mxEtaBin * EtowGeom::mxPhiBin;
227 ETsum=wrkEtow_et[tow];
228 etaSum = wrkEtow_et[tow]*(float)eta;
229 phiSum = wrkEtow_et[tow]*(float)phi;
232 towPlusOne%= maxTowers;
233 ETsum+=wrkEtow_et[towPlusOne];
234 etaSum += wrkEtow_et[towPlusOne]*(float)(eta+1);
235 phiSum += wrkEtow_et[towPlusOne]*(float)phi;
237 tow+=EtowGeom::mxEtaBin;
239 ETsum+=wrkEtow_et[tow];
240 etaSum += wrkEtow_et[tow]*(float)eta;
241 phiSum += wrkEtow_et[tow]*(float)(phi+1);
244 towPlusOne%= maxTowers;
245 ETsum+=wrkEtow_et[towPlusOne];
246 etaSum += wrkEtow_et[towPlusOne]*(float)(eta+1);
247 phiSum += wrkEtow_et[towPlusOne]*(float)(phi+1);
250 *aveEta = etaSum/ETsum;
251 *avePhi = phiSum/ETsum;
259 L2eemcGamma2009::sumET(
int phi,
int eta) {
260 int tow = EtowGeom::mxEtaBin *((phi+EtowGeom::mxPhiBin)%EtowGeom::mxPhiBin) + ((eta+EtowGeom::mxEtaBin)%EtowGeom::mxEtaBin);
262 const int maxTowers = EtowGeom::mxEtaBin * EtowGeom::mxPhiBin;
266 if (etow_used[tow]==0) {
270 towPlusOne%= maxTowers;
271 if (etow_used[towPlusOne]==0) {
272 sum+=wrkEtow_et[towPlusOne];
277 tow+=EtowGeom::mxEtaBin;
279 if (etow_used[tow]==0) {
280 sum+=wrkEtow_et[tow];
283 towPlusOne%= maxTowers;
284 if (etow_used[towPlusOne]==0) {
285 sum+=wrkEtow_et[towPlusOne];
313 const int hitSize=mEveStream_etow[token].get_hitSize();
314 for(i=0;i< hitSize;i++,hit++) {
316 int tower=mRdo2tower[hit->rdo];
317 wrkEtow_tower_index[i]=tower;
318 wrkEtow_et[tower]=hit->et;
322 hA[2]->fill(hitSize);
325 quickSort(wrkEtow_tower_index,0,hitSize-1);
329 for(i=hitSize-1;i >= 0; i--) {
330 if(wrkEtow_et[wrkEtow_tower_index[i]] < par_seedEtThres) {
334 int seedTow=wrkEtow_tower_index[i];
335 int seedEta=seedTow%EtowGeom::mxEtaBin;
336 int seedPhi=seedTow/EtowGeom::mxEtaBin;
348 if (etow_used[seedTow]!=0)
continue;
349 wrkEtow_tower_seed_size++;
354 if(seedEta < EtowGeom::mxEtaBin) {
355 maxET = sumET(seedPhi,seedEta);
356 sum=sumET(seedPhi-1,seedEta);
364 sum=sumET(seedPhi-1,seedEta-1);
370 sum=sumET(seedPhi,seedEta-1);
379 if(maxET<par_clusterEtThres)
continue;
382 hA[6]->fill((
int)(maxET));
383 hA[7]->fill(seedTow);
384 hA[8]->fill(hitSize-i);
387 hA[11]->fill((
int)wrkEtow_et[seedTow]);
388 hA[12]->fill((
int)(10.*wrkEtow_et[seedTow]/maxET));
390 if(etowEve->size>=L2eemcGammaEvent2009::mxClust)
continue;
391 etowEve->clusterET[etowEve->size]=maxET;
392 etowEve->clusterQuad[etowEve->size]=high_quadrant;
393 etowEve->clusterSeedTow[etowEve->size]=seedTow;
402 if (high_quadrant==0) {
403 numNonZeroTow = countNonZeroTow(seedPhi, seedEta);
404 flagUsed(seedPhi, seedEta);
405 averageEtaPhi(seedPhi, seedEta,&clusterPhi, &clusterEta);
406 }
else if (high_quadrant==1) {
407 numNonZeroTow = countNonZeroTow(seedPhi-1, seedEta);
408 flagUsed(seedPhi-1, seedEta);
409 averageEtaPhi(seedPhi-1, seedEta,&clusterPhi, &clusterEta);
410 }
else if (high_quadrant==2) {
411 numNonZeroTow = countNonZeroTow(seedPhi-1, seedEta-1);
412 flagUsed(seedPhi-1, seedEta-1);
413 averageEtaPhi(seedPhi-1, seedEta-1,&clusterPhi, &clusterEta);
414 }
else if (high_quadrant==3) {
415 numNonZeroTow = countNonZeroTow(seedPhi, seedEta-1);
416 flagUsed(seedPhi, seedEta-1);
417 averageEtaPhi(seedPhi, seedEta-1,&clusterPhi, &clusterEta);
419 hA[13]->fill(numNonZeroTow);
421 etowEve->clusterEta[etowEve->size]=clusterEta;
422 etowEve->clusterPhi[etowEve->size]=clusterPhi;
428 hA[3]->fill(wrkEtow_tower_seed_size);
430 hA[4]->fill(etowEve->size);
439 printf(
"dbg=%s, etow-adcL-size=%d\n",getName(),hitSize);
450 L2eemcGamma2009::decisionUser(
int token,
int *myL2Result){
468 if(etowEve->size>= L2eemcGammaEvent2009::mxClust) mhN->fill(5);
469 if(etowEve->isFresh>L2eemcGammaEvent2009::kDataFresh) mhN->fill(6);
472 hA[4]->fill(etowEve->size);
474 for(ic=0;ic<etowEve->size;ic++) {
475 float clustET=etowEve->clusterET[ic];
476 hA[5]->fill((
int)clustET);
478 if(clustET<par_eventEtThres)
continue;
479 hA[6]->fill((
int)clustET);
494 for(ic=0;ic<etowEve->size;ic++) {
495 if(etowEve->clusterET[ic]>maxClusterEt) {
496 maxClusterEt = etowEve->clusterET[ic];
500 if(maxClusterEt>par_eventEtThres) {
501 etowEve->resultBlob.clusterEt=(
unsigned char)(etowEve->clusterET[maxClusterID]*256.0/60.0);
502 etowEve->resultBlob.TowerID=(
unsigned short)etowEve->clusterSeedTow[maxClusterID]+(etowEve->clusterQuad[maxClusterID]<<13);
503 etowEve->resultBlob.meanEtaBin=(
unsigned char)((
int)(etowEve->clusterEta[maxClusterID]*256.0/EtowGeom::mxEtaBin)%256);
504 etowEve->resultBlob.meanPhiBin=(
unsigned char)((
int)(etowEve->clusterPhi[maxClusterID]*256.0/EtowGeom::mxPhiBin)%256);
505 etowEve->resultBlob.isolationSum=(
unsigned char)1;
506 etowEve->resultBlob.Time=(
unsigned char)(mComputeTimeDiff[token]/1000);
507 etowEve->resultBlob.trigger=(
unsigned char)5;
509 etowEve->resultBlob.trigger=(
unsigned char)7;
520 hA[9]->fill(etowEve->clusterQuad[maxClusterID]);
521 hA[10]->fill((
int)etowEve->clusterET[maxClusterID]);
522 if(etowEve->size>= L2eemcGammaEvent2009::mxClust) mhN->fill(15);
530 etowEve->resultBlob.trigger=(
unsigned char)6;
531 etowEve->resultBlob.Time=(
unsigned char)(mComputeTimeDiff[token]/1000);
546 L2eemcGamma2009::finishRunUser() {
550 fprintf(mLogFile,
"finishRunUser-%s bhla bhla\n",getName());
559 L2eemcGamma2009::createHisto() {
563 hA[2]=
new L2Histo(2,
"ETOW-compute: #towers w/ energy /event; x: # ETOW towers; y: counts", 100);
564 hA[3]=
new L2Histo(3,
"ETOW-compute: #seed ....... ", 100);
565 hA[4]=
new L2Histo(4,
"ETOW-decision: #clust ....... ", 50);
567 hA[5]=
new L2Histo(5,
"ETOW-decision: any cluster ; x: ET(GeV)", 30);
568 hA[6]=
new L2Histo(6,
"ETOW-decision: accepted clust ... ; x: ET(GeV)", 30);
570 hA[7]=
new L2Histo(7,
"ETOW: Seed Tower Number", 5000);
571 hA[8]=
new L2Histo(8,
"ETOW: Seed Tower Rank", 20);
572 hA[9]=
new L2Histo(9,
"ETOW: Seed Tower Quadrant", 5);
573 hA[10]=
new L2Histo(10,
"ETOW: Cluster Et Sum GeV", 30);
574 hA[11]=
new L2Histo(11,
"ETOW: Cluster Seed Et GeV", 30);
575 hA[12]=
new L2Histo(12,
"ETOW: 10*Seed Et/Cluster Et GeV", 30);
576 hA[13]=
new L2Histo(13,
"ETOW: # non-zero towers per cluster", 5);
584 L2eemcGamma2009::clearEvent(
int token){
585 memset(wrkEtow_et,0,
sizeof(wrkEtow_et));
586 memset(etow_used,0,
sizeof(etow_used));
587 memset(wrkEtow_tower_seed,0,
sizeof(wrkEtow_tower_seed));
588 wrkEtow_tower_seed_size=0;
592 mEtow[token].isFresh=L2eemcGammaEvent2009::kDataFresh;
600 L2eemcGamma2009::print2(){
602 printf(
"pr2-%s: ---ETOW ADC 2D array, only non-zero\n",getName());
604 for(i=0;i<mxEtow;i++) {
605 if(wrkEtow_et[i]<=0)
continue;
606 int rdo=mTower2rdo[i];
607 float et=wrkEtow_et[i];
608 printf(
" etow: tower=%4d rdo=%4d et=%.3f \n",i,rdo,et);
616 L2eemcGamma2009::print3(){
618 printf(
"pr3-%s: ---seed list, size=%d\n",getName(),wrkEtow_tower_seed_size);
620 for(i=0;i<wrkEtow_tower_seed_size;i++) {
621 int tower=wrkEtow_tower_seed[i];
622 float et=wrkEtow_et[tower];
623 printf(
" etow: i=%4d tower=%4d et=%.3f \n",i,tower,et);
631 L2eemcGamma2009::print4(
int token,
int hitSize){
633 printf(
"print4 IS NOT Fully FUNCTIONAL **********************\n");
634 printf(
"pr1-%s: ---ETOW Sorted ADC list--- size=%d\n",getName(),hitSize);
636 for(i=0;i< hitSize;i++) {
639 float et=wrkEtow_et[wrkEtow_tower_index[i]];
641 printf(
" tower=%2d ",wrkEtow_tower_index[i]);
642 printf(
" etow: i=%2d rdo=%4d adc=%d et=%.3f ene=%.3f\n",i,rdo,adc,et,ene);
void computeUser(int token)