00001 #include <stdio.h>
00002 #include <string.h>
00003 #include <stdlib.h>
00004 #include <time.h>
00005 #include <math.h>
00006 #include <fakeRtsLog.h>
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #ifdef IS_REAL_L2 //in l2-ana environment
00018 #include "../L2algoUtil/L2EmcDb2012.h"
00019 #include "../L2algoUtil/L2Histo.h"
00020 #else //full path needed for cvs'd code
00021 #include "StTriggerUtilities/L2Emulator/L2algoUtil/L2EmcDb2012.h"
00022 #include "StTriggerUtilities/L2Emulator/L2algoUtil/L2Histo.h"
00023 #include "StTriggerUtilities/L2Emulator/L2algoUtil/L2EmcGeom2012.h"
00024 #endif
00025
00026 #include "L2eemcGamma2012.h"
00027
00028
00029
00030
00031 void L2eemcGamma2012::quickSort(int array[], int start, int end)
00032 {
00033 int i = start;
00034 int k = end;
00035
00036 if (end - start >= 1)
00037 {
00038 float pivot = wrkEtow_et[array[start]];
00039
00040 while (k > i)
00041 {
00042 while (wrkEtow_et[array[i]] <= pivot && i <= end && k > i)
00043 i++;
00044 while (wrkEtow_et[array[k]] > pivot && k >= start && k >= i)
00045 k--;
00046 if (k > i)
00047 swap(array, i, k);
00048 }
00049 swap(array, start, k);
00050
00051 quickSort(array, start, k - 1);
00052 quickSort(array, k + 1, end);
00053 }
00054 else
00055 {
00056 return;
00057 }
00058 }
00059
00060 void L2eemcGamma2012::swap(int array[], int index1, int index2)
00061
00062
00063 {
00064 int temp = array[index1];
00065 array[index1] = array[index2];
00066 array[index2] = temp;
00067 }
00068
00069
00070
00071
00072
00073 L2eemcGamma2012::L2eemcGamma2012(const char* name, const char *uid, L2EmcDb2012* db, L2EmcGeom2012 *geoX, char* outDir, int resOff) : L2VirtualAlgo2012( name, uid, db, outDir, false,true, resOff) {
00074
00075
00076
00077 mGeom=geoX;
00078 if (!mGeom)
00079 criticalError("L2eemcGamma is broken -- can't find geom.");
00080
00081 setMaxHist(16);
00082 createHisto();
00083
00084
00085 if (sizeof(L2gammaResult2012)!= L2gammaResult2012::mySizeChar)
00086
00087 criticalError("L2eemcGamma has failed consistency check. sizeof(L2gammaResult2012)!= L2gammaResult2012::mySizeChar");
00088
00089
00090
00091 }
00092
00093
00094
00095 int
00096 L2eemcGamma2012::initRunUser( int runNo, int *rc_ints, float *rc_floats) {
00097
00098
00099 par_dbg = rc_ints[0];
00100 par_RndAcceptPrescale = rc_ints[1];
00101 par_seedEtThres = rc_floats[0];
00102 par_clusterEtThres= rc_floats[1];
00103 par_eventEtThres = rc_floats[2];
00104
00105
00106 int kBad=0;
00107 kBad+=0x00001 * (par_seedEtThres<1.0);
00108 kBad+=0x00002 * (par_clusterEtThres<par_seedEtThres);
00109
00110 if (mLogFile) {
00111 fprintf(mLogFile,"L2%s algorithm initRun(R=%d), compiled: %s , %s\n params:\n",getName(),mRunNumber,__DATE__,__TIME__);
00112 fprintf(mLogFile," - use Thres/GeV seed=%.2f, clusterList=%.2f debug=%d, prescale=%d\n", par_seedEtThres,par_clusterEtThres, par_dbg,par_RndAcceptPrescale);
00113
00114 fprintf(mLogFile," - accept event cluster Thres/GeV=%.2f\n",par_eventEtThres);
00115 fprintf(mLogFile,"initRun() params checked for consistency, Error flag=0x%04x\n",kBad);
00116 }
00117
00118
00119
00120 int i;
00121 for (i=0; i<mxHA;i++) if(hA[i])hA[i]->reset();
00122
00123 if(kBad) return -1*kBad;
00124
00125 memset(mEtow,0,sizeof(mEtow));
00126
00127
00128 char txt[1000];
00129 sprintf(txt,"ETOW-compute: #seed towers ET>%.2f GeV / event; x: # ETOW towers; y: counts",par_seedEtThres);
00130 hA[3]->setTitle(txt);
00131
00132 sprintf(txt,"ETOW-decision: #cluster ET>%.2f GeV / event ; x: # ETOW towers; y: counts",par_clusterEtThres);
00133 hA[4]->setTitle(txt);
00134
00135 sprintf(txt,"ETOW-decision: acc cluster ET>%.2f GeV; x:cluster ET(GeV) ; y: counts",par_eventEtThres);
00136 hA[4]->setTitle(txt);
00137
00138
00139
00140 for ( int index=0; index<EmcDbIndexMax; index++ )
00141 {
00142 const L2EmcDb2012::EmcCDbItem *x = mDb->getByIndex(index);
00143 if ( x==0 ) continue;
00144 if ( !mDb->isETOW(x) ) continue;
00145 int sec = x->sec - 1;
00146 int sub = 8192;
00147 sub = x->sub - 'A';
00148 int eta = x->eta - 1;
00149 int phi = EtowGeom::mxSubs *sec + sub;
00150 int tow = EtowGeom::mxEtaBin *phi + eta;
00151 int rdo = x->rdo;
00152 if (tow<0 || tow>mxEtow || rdo<0 || rdo>mxEtow) return -101;
00153
00154 mTower2rdo[ tow ] = rdo;
00155 mRdo2tower[ rdo ] = tow;
00156 }
00157 return 0;
00158
00159 }
00160
00161
00162
00163 int
00164 L2eemcGamma2012::countNonZeroTow(int phi, int eta) {
00165 int tow = EtowGeom::mxEtaBin *((phi+EtowGeom::mxPhiBin)%EtowGeom::mxPhiBin) + ((eta+EtowGeom::mxEtaBin)%EtowGeom::mxEtaBin);
00166 const int maxTowers = EtowGeom::mxEtaBin * EtowGeom::mxPhiBin;
00167 int towPlusOne;
00168 int count=0;
00169
00170 if (etow_used[tow]==0) {
00171 if(wrkEtow_et[tow]!=0) count++;
00172 }
00173
00174 towPlusOne = tow+1;
00175 towPlusOne%= maxTowers;
00176 if (etow_used[towPlusOne]==0) {
00177 if(wrkEtow_et[towPlusOne]!=0) count++;
00178 }
00179
00180 tow+=EtowGeom::mxEtaBin;
00181 tow%=maxTowers;
00182 if (etow_used[tow]==0) {
00183 if(wrkEtow_et[tow]!=0) count++;
00184 }
00185
00186 towPlusOne = tow+1;
00187 towPlusOne%= maxTowers;
00188 if (etow_used[towPlusOne]==0) {
00189 if(wrkEtow_et[towPlusOne]!=0) count++;
00190 }
00191 return count;
00192 }
00193
00194
00195
00196
00197 void
00198 L2eemcGamma2012::flagUsed(int phi, int eta) {
00199 int tow = EtowGeom::mxEtaBin *((phi+EtowGeom::mxPhiBin)%EtowGeom::mxPhiBin) + ((eta+EtowGeom::mxEtaBin)%EtowGeom::mxEtaBin);
00200 const int maxTowers = EtowGeom::mxEtaBin * EtowGeom::mxPhiBin;
00201 int towPlusOne;
00202 etow_used[tow] = 1;
00203
00204 towPlusOne = tow+1;
00205 towPlusOne%= maxTowers;
00206 etow_used[towPlusOne] = 1;
00207
00208 tow+=EtowGeom::mxEtaBin;
00209 tow%=maxTowers;
00210 etow_used[tow] = 1;
00211
00212 towPlusOne = tow+1;
00213 towPlusOne%= maxTowers;
00214 etow_used[towPlusOne] = 1;
00215 }
00216
00217
00218
00219 void
00220 L2eemcGamma2012::averageEtaPhi(int phi, int eta, float *avePhi, float *aveEta) {
00221 int tow = EtowGeom::mxEtaBin *((phi+EtowGeom::mxPhiBin)%EtowGeom::mxPhiBin) + ((eta+EtowGeom::mxEtaBin)%EtowGeom::mxEtaBin);
00222 const int maxTowers = EtowGeom::mxEtaBin * EtowGeom::mxPhiBin;
00223 int towPlusOne;
00224 float ETsum;
00225 float etaSum;
00226 float phiSum;
00227
00228 ETsum=wrkEtow_et[tow];
00229 etaSum = wrkEtow_et[tow]*(float)eta;
00230 phiSum = wrkEtow_et[tow]*(float)phi;
00231
00232 towPlusOne = tow+1;
00233 towPlusOne%= maxTowers;
00234 ETsum+=wrkEtow_et[towPlusOne];
00235 etaSum += wrkEtow_et[towPlusOne]*(float)(eta+1);
00236 phiSum += wrkEtow_et[towPlusOne]*(float)phi;
00237
00238 tow+=EtowGeom::mxEtaBin;
00239 tow%=maxTowers;
00240 ETsum+=wrkEtow_et[tow];
00241 etaSum += wrkEtow_et[tow]*(float)eta;
00242 phiSum += wrkEtow_et[tow]*(float)(phi+1);
00243
00244 towPlusOne = tow+1;
00245 towPlusOne%= maxTowers;
00246 ETsum+=wrkEtow_et[towPlusOne];
00247 etaSum += wrkEtow_et[towPlusOne]*(float)(eta+1);
00248 phiSum += wrkEtow_et[towPlusOne]*(float)(phi+1);
00249
00250
00251 *aveEta = etaSum/ETsum;
00252 *avePhi = phiSum/ETsum;
00253 }
00254
00255
00256
00257
00258 float
00259 L2eemcGamma2012::sumET(int phi, int eta) {
00260 int tow = EtowGeom::mxEtaBin *((phi+EtowGeom::mxPhiBin)%EtowGeom::mxPhiBin) + ((eta+EtowGeom::mxEtaBin)%EtowGeom::mxEtaBin);
00261
00262 const int maxTowers = EtowGeom::mxEtaBin * EtowGeom::mxPhiBin;
00263 int towPlusOne;
00264 float sum;
00265 sum=0;
00266 if (etow_used[tow]==0) {
00267 sum=wrkEtow_et[tow];
00268 }
00269 towPlusOne = tow+1;
00270 towPlusOne%= maxTowers;
00271 if (etow_used[towPlusOne]==0) {
00272 sum+=wrkEtow_et[towPlusOne];
00273 }
00274
00275 tow+=EtowGeom::mxEtaBin;
00276 tow%=maxTowers;
00277 if (etow_used[tow]==0) {
00278 sum+=wrkEtow_et[tow];
00279 }
00280 towPlusOne = tow+1;
00281 towPlusOne%= maxTowers;
00282 if (etow_used[towPlusOne]==0) {
00283 sum+=wrkEtow_et[towPlusOne];
00284 }
00285 return sum;
00286 }
00287
00288
00289
00290
00291 void
00292 L2eemcGamma2012::computeUser(int token){
00293
00294
00295
00296
00297
00298
00299
00300 clearEvent(token);
00301
00302
00303 int i;
00304
00305 L2eemcGammaEvent2012 *etowEve=mEtow+token;
00306 const HitTower1 *hit=mEveStream_etow[token].get_hits();
00307 const int hitSize=mEveStream_etow[token].get_hitSize();
00308 for(i=0;i< hitSize;i++,hit++) {
00309
00310 int tower=mRdo2tower[hit->rdo];
00311 wrkEtow_tower_index[i]=tower;
00312 wrkEtow_et[tower]=hit->et;
00313
00314
00315 }
00316 hA[2]->fill(hitSize);
00317
00318
00319 quickSort(wrkEtow_tower_index,0,hitSize-1);
00320
00321
00322
00323 for(i=hitSize-1;i >= 0; i--) {
00324 if(wrkEtow_et[wrkEtow_tower_index[i]] < par_seedEtThres) {
00325 i = -1;
00326 continue;
00327 }
00328 int seedTow=wrkEtow_tower_index[i];
00329 int seedEta=seedTow%EtowGeom::mxEtaBin;
00330 int seedPhi=seedTow/EtowGeom::mxEtaBin;
00331
00332
00334
00335
00336
00337
00338
00339 float maxET;
00340 float sum;
00341
00342 if (etow_used[seedTow]!=0) continue;
00343 wrkEtow_tower_seed_size++;
00344
00345 int high_quadrant;
00346 high_quadrant = 0;
00347 maxET = 0;
00348 if(seedEta < EtowGeom::mxEtaBin) {
00349 maxET = sumET(seedPhi,seedEta);
00350 sum=sumET(seedPhi-1,seedEta);
00351 if(maxET<sum) {
00352 maxET=sum;
00353 high_quadrant = 1;
00354 }
00355 }
00356 if (seedEta > 0) {
00357 sum=sumET(seedPhi-1,seedEta-1);
00358 if(maxET<sum) {
00359 maxET=sum;
00360 high_quadrant = 2;
00361 }
00362 sum=sumET(seedPhi,seedEta-1);
00363 if(maxET<sum) {
00364 maxET=sum;
00365 high_quadrant = 3;
00366 }
00367 }
00368
00369 if(maxET<par_clusterEtThres)continue;
00370
00371
00372 hA[6]->fill((int)(maxET));
00373 hA[7]->fill(seedTow);
00374 hA[8]->fill(hitSize-i);
00375
00376
00377 hA[11]->fill((int)wrkEtow_et[seedTow]);
00378 hA[12]->fill((int)(10.*wrkEtow_et[seedTow]/maxET));
00379
00380 if(etowEve->size>=L2eemcGammaEvent2012::mxClust) continue;
00381 etowEve->clusterET[etowEve->size]=maxET;
00382 etowEve->clusterQuad[etowEve->size]=high_quadrant;
00383 etowEve->clusterSeedTow[etowEve->size]=seedTow;
00384
00385
00386
00387
00388
00389 int numNonZeroTow=0;
00390 float clusterEta;
00391 float clusterPhi;
00392 if (high_quadrant==0) {
00393 numNonZeroTow = countNonZeroTow(seedPhi, seedEta);
00394 flagUsed(seedPhi, seedEta);
00395 averageEtaPhi(seedPhi, seedEta,&clusterPhi, &clusterEta);
00396 } else if (high_quadrant==1) {
00397 numNonZeroTow = countNonZeroTow(seedPhi-1, seedEta);
00398 flagUsed(seedPhi-1, seedEta);
00399 averageEtaPhi(seedPhi-1, seedEta,&clusterPhi, &clusterEta);
00400 } else if (high_quadrant==2) {
00401 numNonZeroTow = countNonZeroTow(seedPhi-1, seedEta-1);
00402 flagUsed(seedPhi-1, seedEta-1);
00403 averageEtaPhi(seedPhi-1, seedEta-1,&clusterPhi, &clusterEta);
00404 } else if (high_quadrant==3) {
00405 numNonZeroTow = countNonZeroTow(seedPhi, seedEta-1);
00406 flagUsed(seedPhi, seedEta-1);
00407 averageEtaPhi(seedPhi, seedEta-1,&clusterPhi, &clusterEta);
00408 }
00409 hA[13]->fill(numNonZeroTow);
00410 etowEve->clusterEta[etowEve->size]=clusterEta;
00411 etowEve->clusterPhi[etowEve->size]=clusterPhi;
00412
00413 etowEve->size++;
00414
00415 }
00416
00417 hA[3]->fill(wrkEtow_tower_seed_size);
00418 hA[4]->fill(etowEve->size);
00419
00420
00421 if(par_dbg>0){
00422 LOG(DBG,"dbg=%s, etow-adcL-size=%d\n",getName(),hitSize);
00423
00424 print2();
00425 print3();
00426 }
00427 }
00428
00429
00430
00431
00432 bool
00433 L2eemcGamma2012::decisionUser(int token, int *myL2Result){
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444 L2eemcGammaEvent2012 *etowEve=mEtow+token;
00445
00446 int ic;
00447
00448
00449 if(etowEve->size>= L2eemcGammaEvent2012::mxClust) mhN->fill(5);
00450 if(etowEve->isFresh>L2eemcGammaEvent2012::kDataFresh) mhN->fill(6);
00451 etowEve->isFresh++;
00452
00453 hA[4]->fill(etowEve->size);
00454
00455 for(ic=0;ic<etowEve->size;ic++) {
00456 float clustET=etowEve->clusterET[ic];
00457 hA[5]->fill((int)clustET);
00458 if(clustET<par_eventEtThres) continue;
00459 hA[6]->fill((int)clustET);
00460 }
00461
00462
00463
00464
00465
00466
00467 float maxClusterEt;
00468 maxClusterEt=-1;
00469 int maxClusterID;
00470 maxClusterID =-1;
00471
00472
00473 for(ic=0;ic<etowEve->size;ic++) {
00474 if(etowEve->clusterET[ic]>maxClusterEt) {
00475 maxClusterEt = etowEve->clusterET[ic];
00476 maxClusterID = ic;
00477 }
00478 }
00479 if(maxClusterEt>par_eventEtThres) {
00480 etowEve->resultBlob.clusterEt=(unsigned char)(etowEve->clusterET[maxClusterID]*256.0/60.0);
00481 etowEve->resultBlob.TowerID=(unsigned short)etowEve->clusterSeedTow[maxClusterID]+(etowEve->clusterQuad[maxClusterID]<<13);
00482 etowEve->resultBlob.meanEtaBin=(unsigned char)((int)(etowEve->clusterEta[maxClusterID]*256.0/EtowGeom::mxEtaBin)%256);
00483 etowEve->resultBlob.meanPhiBin=(unsigned char)((int)(etowEve->clusterPhi[maxClusterID]*256.0/EtowGeom::mxPhiBin)%256);
00484 etowEve->resultBlob.isolationSum=(unsigned char)1;
00485 etowEve->resultBlob.Time=(unsigned char)(mComputeTimeDiff[token]/1000);
00486 etowEve->resultBlob.trigger=(unsigned char)5;
00487 if (mRandomAccept) {
00488 etowEve->resultBlob.trigger=(unsigned char)7;
00489 }
00490
00491 hA[9]->fill(etowEve->clusterQuad[maxClusterID]);
00492 hA[10]->fill((int)etowEve->clusterET[maxClusterID]);
00493 if(etowEve->size>= L2eemcGammaEvent2012::mxClust) mhN->fill(15);
00494 mhN->fill(15);
00495 memcpy(myL2Result,&(etowEve->resultBlob),sizeof(L2gammaResult2012));
00496 return true;
00497 }
00498
00499
00500 if (mRandomAccept) {
00501 etowEve->resultBlob.trigger=(unsigned char)6;
00502 etowEve->resultBlob.Time=(unsigned char)(mComputeTimeDiff[token]/1000);
00503 mhN->fill(16);
00504 memcpy(myL2Result,&(etowEve->resultBlob),sizeof(L2gammaResult2012));
00505 return false;
00506 }
00507
00508 memcpy(myL2Result,&(etowEve->resultBlob),sizeof(L2gammaResult2012));
00509 return false;
00510
00511 }
00512
00513
00514
00515
00516 void
00517 L2eemcGamma2012::finishRunUser() {
00518
00519
00520 if (mLogFile){
00521 fprintf(mLogFile,"finishRunUser-%s bhla bhla\n",getName());
00522 }
00523
00524 }
00525
00526
00527
00528
00529 void
00530 L2eemcGamma2012::createHisto() {
00531 setMaxHist(15);
00532
00533
00534 hA[2]=new L2Histo(2,"ETOW-compute: #towers w/ energy /event; x: # ETOW towers; y: counts", 100);
00535 hA[3]=new L2Histo(3,"ETOW-compute: #seed ....... ", 100);
00536 hA[4]=new L2Histo(4,"ETOW-decision: #clust ....... ", 50);
00537
00538 hA[5]=new L2Histo(5,"ETOW-decision: any cluster ; x: ET(GeV)", 30);
00539 hA[6]=new L2Histo(6,"ETOW-decision: accepted clust ... ; x: ET(GeV)", 30);
00540
00541 hA[7]=new L2Histo(7,"ETOW: Seed Tower Number", 5000);
00542 hA[8]=new L2Histo(8,"ETOW: Seed Tower Rank", 20);
00543 hA[9]=new L2Histo(9,"ETOW: Seed Tower Quadrant", 5);
00544 hA[10]=new L2Histo(10,"ETOW: Cluster Et Sum GeV", 30);
00545 hA[11]=new L2Histo(11,"ETOW: Cluster Seed Et GeV", 30);
00546 hA[12]=new L2Histo(12,"ETOW: 10*Seed Et/Cluster Et GeV", 30);
00547 hA[13]=new L2Histo(13,"ETOW: # non-zero towers per cluster", 5);
00548
00549 }
00550
00551
00552
00553 void
00554 L2eemcGamma2012::clearEvent(int token){
00555 memset(wrkEtow_et,0,sizeof(wrkEtow_et));
00556 memset(etow_used,0,sizeof(etow_used));
00557 memset(wrkEtow_tower_seed,0,sizeof(wrkEtow_tower_seed));
00558 wrkEtow_tower_seed_size=0;
00559
00560 mEtow[token].size=0;
00561
00562 mEtow[token].isFresh=L2eemcGammaEvent2012::kDataFresh;
00563
00564 memset(&(mEtow[token].resultBlob),0, sizeof(L2gammaResult2012));
00565 }
00566
00567
00568
00569 void
00570 L2eemcGamma2012::print2(){
00571 int i;
00572 printf("pr2-%s: ---ETOW ADC 2D array, only non-zero\n",getName());
00573
00574 for(i=0;i<mxEtow;i++) {
00575 if(wrkEtow_et[i]<=0) continue;
00576 int rdo=mTower2rdo[i];
00577 float et=wrkEtow_et[i];
00578 printf(" etow: tower=%4d rdo=%4d et=%.3f \n",i,rdo,et);
00579 }
00580
00581 }
00582
00583
00584
00585 void
00586 L2eemcGamma2012::print3(){
00587 int i;
00588 printf("pr3-%s: ---seed list, size=%d\n",getName(),wrkEtow_tower_seed_size);
00589
00590 for(i=0;i<wrkEtow_tower_seed_size;i++) {
00591 int tower=wrkEtow_tower_seed[i];
00592 float et=wrkEtow_et[tower];
00593 printf(" etow: i=%4d tower=%4d et=%.3f \n",i,tower,et);
00594 }
00595
00596 }
00597
00598
00599
00600 void
00601 L2eemcGamma2012::print4(int token, int hitSize){
00602 int i;
00603 printf("print4 IS NOT Fully FUNCTIONAL **********************\n");
00604 printf("pr1-%s: ---ETOW Sorted ADC list--- size=%d\n",getName(),hitSize);
00605
00606 for(i=0;i< hitSize;i++) {
00607 int adc=0;
00608 int rdo=0;
00609 float et=wrkEtow_et[wrkEtow_tower_index[i]];
00610 float ene=0;
00611 printf(" tower=%2d ",wrkEtow_tower_index[i]);
00612 printf(" etow: i=%2d rdo=%4d adc=%d et=%.3f ene=%.3f\n",i,rdo,adc,et,ene);
00613 }
00614 }
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655