StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
L2upsilon2006.cc
1 //
2 // Pibero Djawwotho <pibero@iucf.indiana.edu>
3 // Indiana University
4 // October 23, 2007
5 //
6 
7 #include <cmath>
8 #include <cstdlib>
9 #include <algorithm>
10 #include <functional>
11 #include <cstring>
12 #include <cassert>
13 
14 #include "L2upsilonResult2006.h"
15 #include "L2upsilon2006.hh"
16 #ifdef IS_REAL_L2 //in l2-ana environment
17  #include "trgStructures.h"
18  #include "../L2algoUtil/L2EmcDb.h"
19  #include "../L2algoUtil/L2Histo.h"
20 #else
21  #include "StDaqLib/TRG/trgStructures.h"
22  #include "StTriggerUtilities/L2Emulator/L2algoUtil/L2EmcDb.h"
23  #include "StTriggerUtilities/L2Emulator/L2algoUtil/L2Histo.h"
24 #endif
25 
26 #define L0_HT_TRIGGER_BIT_SHIFT 4
27 #define CPU_MHZ 1603.680
28 
29 L2upsilon2006::L2upsilon2006(const char* name, L2EmcDb* db, char* outDir, int resOff) :
30  L2VirtualAlgo(name, db, outDir, resOff), mLogFile(stdout), mRunNumber(0), mUnfinished(false)
31 {
32  createHistograms();
33 }
34 
35 L2upsilon2006::~L2upsilon2006()
36 {
37  deleteHistograms();
38 }
39 
40 int L2upsilon2006::initRun(int runNumber, int* userInt, float* userFloat)
41 {
42  //
43  // Open log file and write out algorithm parameters
44  //
45  if (mUnfinished) {
46  fprintf(mLogFile, "L2upsilon2006: WARNING - finishRun() was not called for last run\n");
47  finishRun();
48  }
49 
50  mRunNumber = runNumber;
51  mUnfinished = true;
52 
53  char filename[FILENAME_MAX];
54  sprintf(filename, "%s/run%d.l2ups.log", mOutDir, runNumber);
55 
56  if (FILE* fp = fopen(filename, "w"))
57  mLogFile = fp;
58  else
59  printf("L2upsilon2006: WARNING - Can't open log file %s\n", filename);
60 
61  fprintf(mLogFile, "L2upsilon2006::initRun(char* name=\"%s\", int runNumber=%d, int* userInt=%p, float* userFloat=%p)\n", mName, runNumber, userInt, userFloat);
62 
63  //
64  // Set float and int run control parameters
65  //
66  mMinL0ClusterEnergy = userFloat[0];
67  mMinL2ClusterEnergy = userFloat[1];
68  mMinInvMass = userFloat[2];
69  mMaxInvMass = userFloat[3];
70  mMaxCosTheta = userFloat[4];
71 
72  mL0SeedThreshold = userInt[0];
73  mL2SeedThreshold = userInt[1];
74  mUseCtb = userInt[2];
75  mUseVertexZ = userInt[3];
76  mNumberOfTowersPerCluster = userInt[4];
77 
78  fprintf(mLogFile, "Run Number: %d\n", runNumber);
79  fprintf(mLogFile, "Start Time: %s", timeString());
80  fprintf(mLogFile, "L0 seed threshold: %d\n", mL0SeedThreshold);
81  fprintf(mLogFile, "L2 seed threshold: %d\n", mL2SeedThreshold);
82  fprintf(mLogFile, "L0 cluster energy threshold [GeV]: %f\n", mMinL0ClusterEnergy);
83  fprintf(mLogFile, "L2 cluster energy threshold [GeV]: %f\n", mMinL2ClusterEnergy);
84  fprintf(mLogFile, "Number of towers per cluster: %d\n", mNumberOfTowersPerCluster);
85  fprintf(mLogFile, "Use CTB: %d\n", mUseCtb);
86  fprintf(mLogFile, "Use vertex Z: %d\n", mUseVertexZ);
87  fprintf(mLogFile, "Minimum invariant mass [GeV]: %f\n", mMinInvMass);
88  fprintf(mLogFile, "Maximum invariant mass [GeV]: %f\n", mMaxInvMass);
89  fprintf(mLogFile, "Maximum cos(theta): %f\n", mMaxCosTheta);
90 
91  //
92  // Reset event counters
93  //
94  mEventsSeen = 0;
95  mEventsAccepted = 0;
96 
97  mL0SeedThreshold <<= L0_HT_TRIGGER_BIT_SHIFT;
98 
99  // Get pedestals, gains, status, ...
100  mDb->initRun(runNumber);
101  memset(bemcTower, 0, sizeof(bemcTower));
102  memset(mSoftIdToRdo, 0, sizeof(mSoftIdToRdo));
103  memset(mPhiEtaToRdo, 0, sizeof(mPhiEtaToRdo));
104  for (int i = 0; i < EmcDbIndexMax; ++i) {
105  const L2EmcDb::EmcCDbItem* x = mDb->getByIndex(i);
106  if (mDb->isBTOW(x)) {
107  int softId = atoi(x->tube + 2); // 1-4800
108  int eta = x->eta - 1; // 0-39
109  int phi = (x->sec - 1) * 10 + x->sub - 'a'; // 0-119
110  bemcTower[x->rdo].softId = softId;
111  bemcTower[x->rdo].gain = x->gain;
112  bemcTower[x->rdo].eta = eta;
113  bemcTower[x->rdo].phi = phi;
114  bemcTower[x->rdo].pedestal = x->ped;
115  bemcTower[x->rdo].stat = x->stat;
116  bemcTower[x->rdo].fail = x->fail;
117  mSoftIdToRdo[softId] = x->rdo;
118  mPhiEtaToRdo[phi][eta] = x->rdo;
119  }
120  }
121 
122  // Get neighbors
123  for (int i = 0; i < EmcDbIndexMax; ++i) {
124  const L2EmcDb::EmcCDbItem* x = mDb->getByIndex(i);
125  if (mDb->isBTOW(x)) {
126  int eta = x->eta - 1; // 0-39
127  int phi = (x->sec - 1) * 10 + x->sub - 'a'; // 0-119
128  int n = 0;
129  for (int deta = -1; deta <= 1; ++deta) {
130  int eta2 = eta + deta;
131  if (eta2 < 0 || eta2 >= 40) continue;
132  for (int dphi = -1; dphi <= 1; ++dphi) {
133  if (deta == 0 && dphi == 0) continue;
134  int phi2 = phi + dphi;
135  phi2 += 120;
136  phi2 %= 120;
137  int rdo2 = mPhiEtaToRdo[phi2][eta2];
138  bemcTower[x->rdo].neighbor[n++] = rdo2;
139  }
140  }
141  bemcTower[x->rdo].numberOfNeighbors = n;
142  }
143  }
144 
145  resetHistograms();
146  timer.start();
147 
148  return 0;
149 }
150 
151 void L2upsilon2006::readGeomXYZ(const char *fname){
152  printf(" Get BTOW towers x, y, z from=%s=\n",fname);
153 
154  FILE* fp = fopen(fname, "r");
155  assert(fp);
156  int id;
157  float x, y, z;
158  int ret;
159  while ((ret = fscanf(fp, "%d %f %f %f", &id, &x, &y, &z)) != EOF) {
160  assert(ret==4);
161  int rdo = mSoftIdToRdo[id];
162  bemcTower[rdo].x = x;
163  bemcTower[rdo].y = y;
164  bemcTower[rdo].z = z;
165  }
166  fclose(fp);
167 }
168 
169 
170 bool L2upsilon2006::doEvent(int L0trg, int eventNumber, TrgDataType* trgData,
171  int bemcIn, unsigned short* bemcData,
172  int eemcIn, unsigned short* eemcData)
173 {
174 
175 #if 0
176  printf("L2upsilon2006::doEvent(int L0trg=%d, int eventNumber=%d, TrgDataType* trgData=%p, int bemcIn=%d, unsigned short* bemcData=%p, int eemcIn=%d, unsigned short* eemcData=%p)\n", L0trg, eventNumber, trgData, bemcIn, bemcData, eemcIn, eemcData);
177 #endif
178  rdtscl_macro(mEveTimeStart);
179  mAccept=false;
180  if (trgData && bemcData) {// UUUU , I inverted logic to positive
181 
182  //
183  // We have a valid event
184  //
185  unsigned int timeStart;
186  unsigned int timeStop;
187 
188  rdtscl_macro(timeStart);
189  hL0rate->fill(timer.time());
190  ++mEventsSeen;
191  this->trgData = trgData;
192  this->bemcData = bemcData;
193 
194  //
195  // Save the trigger data and BEMC data pointers for use in member functions
196  //
197  vector<int> L0Seeds;
198  vector<int> L2Seeds;
199 
200  L0Seeds.reserve(100);
201  L2Seeds.reserve(100);
202  findSeedTowers(L0Seeds, L2Seeds);
203 
204  L2upsilonResult2006 result;
205  for (size_t i = 0; i < L0Seeds.size(); ++i) {
206  const int id1 = L0Seeds[i];
207  calcCluster(id1);
208  if (bemcCluster[id1].E < mMinL0ClusterEnergy || bemcCluster[id1].E == 0) continue;
209  for (size_t j = 0; j < L2Seeds.size(); ++j) {
210  const int id2 = L2Seeds[j];
211  if (id1 == id2) continue;
212  calcCluster(id2);
213  if (bemcCluster[id2].E < mMinL2ClusterEnergy || bemcCluster[id2].E == 0) continue;
214  float cosTheta = (bemcCluster[id1].x * bemcCluster[id2].x +
215  bemcCluster[id1].y * bemcCluster[id2].y +
216  bemcCluster[id1].z * bemcCluster[id2].z);
217  if (cosTheta > mMaxCosTheta) continue;
218  float invMass = sqrt(2 * bemcCluster[id1].E * bemcCluster[id2].E * (1 - cosTheta));
219  if (mMinInvMass < invMass && invMass < mMaxInvMass) {
220  //
221  // We have a candidate
222  //
223  hL2rate->fill(timer.time());
224 
225  result.numberOfL0SeedTowers = L0Seeds.size();
226  result.numberOfL2SeedTowers = L2Seeds.size();
227  result.invariantMass = invMass;
228  result.eventsSeen = mEventsSeen;
229  result.eventsAccepted = ++mEventsAccepted;
230  result.energyOfL0Cluster = bemcCluster[id1].E;
231  result.energyOfL2Cluster = bemcCluster[id2].E;
232  rdtscl_macro(timeStop);
233  result.processingTime = (timeStop > timeStart) ? timeStop - timeStart : 0;
234 
235  //
236  // Fill histograms
237  //
238  hL0SeedTowers->fill(bemcTower[id1].softId-1);
239  hL2SeedTowers->fill(bemcTower[id2].softId-1);
240  hNumberOfL0Seeds->fill(L0Seeds.size());
241  hNumberOfL2Seeds->fill(L2Seeds.size());
242  hInvMass->fill(int(invMass*10));
243  hTime->fill(int(0.5*result.processingTime/CPU_MHZ));
244  hEnergyOfL0Cluster->fill(int(bemcCluster[id1].E*10));
245  hEnergyOfL2Cluster->fill(int(bemcCluster[id2].E*10));
246  hCosTheta->fill(int((cosTheta+1)*50));
247 
248  memcpy(&trgData->TrgSum.L2Result[mResultOffset], &result, sizeof(result));
249 
250 #ifndef IS_REAL_L2
251  print();
252 #endif
253  mAccept=true;
254  goto doEventEnd;
255  }
256  }
257  }
258 
259  //
260  // No candidate found
261  //
262  result.numberOfL0SeedTowers = L0Seeds.size();
263  result.numberOfL2SeedTowers = L2Seeds.size();
264  result.invariantMass = 0;
265  result.eventsSeen = mEventsSeen;
266  result.eventsAccepted = mEventsAccepted;
267  result.energyOfL0Cluster = 0;
268  result.energyOfL2Cluster = 0;
269  rdtscl_macro(timeStop);
270  result.processingTime = (timeStop > timeStart) ? timeStop - timeStart : 0;
271  memcpy(&trgData->TrgSum.L2Result[mResultOffset], &result, sizeof(result));
272 
273 #ifndef IS_REAL_L2
274  print();
275 #endif
276  } // end of UUUU
277  doEventEnd:
278  rdtscl_macro(mEveTimeStop);
279  mEveTimeDiff=mEveTimeStop-mEveTimeStart;
280  int kTick=mEveTimeDiff/1000;
281  // printf("uu=%f t1=%d t2=%d \n",mEveTimeDiff/1000.,mEveTimeStart,mEveTimeStop);
282  mhT->fill(kTick);
283 
284  return mAccept;
285 }
286 
287 void L2upsilon2006::finishRun()
288 {
289  mUnfinished = false;
290  fprintf(mLogFile, (char*)"L2upsilon2006::finishRun()\n");
291  finishCommonHistos() ;
292  //
293  // Write out events summary and close log file
294  //
295  fprintf(mLogFile, "Stop Time: %s", timeString());
296  fprintf(mLogFile, "Events seen: %d\n", mEventsSeen);
297  fprintf(mLogFile, "Events accepted: %d\n", mEventsAccepted);
298 
299  writeHistograms();
300 
301  if (mLogFile != stdout) fclose(mLogFile);
302 }
303 
304 void L2upsilon2006::findSeedTowers(vector<int>& L0Seeds, vector<int>& L2Seeds)
305 {
306  for (int rdo = 0; rdo < 4800; ++rdo) {
307  if (bemcTower[rdo].fail) continue;
308  if (bemcData[rdo] - bemcTower[rdo].pedestal >= mL2SeedThreshold) {
309  hHighTowers->fill(bemcTower[rdo].softId-1);
310  L2Seeds.push_back(rdo);
311  if (bemcData[rdo] - bemcTower[rdo].pedestal >= mL0SeedThreshold)
312  L0Seeds.push_back(rdo);
313  }
314  }
315 }
316 
317 void L2upsilon2006::calcCluster(int rdo)
318 {
319  bemcCluster[rdo].x = 0;
320  bemcCluster[rdo].y = 0;
321  bemcCluster[rdo].z = 0;
322  bemcCluster[rdo].E = 0;
323 
324  //
325  // Calculate the energies of all neighbors
326  //
327  float energies[8];
328  for (int i = 0; i < bemcTower[rdo].numberOfNeighbors; ++i) {
329  int j = bemcTower[rdo].neighbor[i];
330  if (bemcTower[j].fail || bemcData[j] < bemcTower[j].pedestal) {
331  energies[i] = 0;
332  }
333  else {
334  energies[i] = (bemcData[j] - bemcTower[j].pedestal) / bemcTower[j].gain;
335  }
336  }
337 
338  //
339  // Sort the high tower neighbors in order of descending energy.
340  // Use an insertion sort algorithm which is fairly efficient
341  // for small arrays.
342  //
343  for (int i = 1; i < bemcTower[rdo].numberOfNeighbors; ++i) {
344  int j = i;
345  float indexE = energies[j];
346  int indexId = bemcTower[rdo].neighbor[j];
347  while (j > 0 && indexE > energies[j-1]) {
348  energies[j] = energies[j-1];
349  bemcTower[rdo].neighbor[j] = bemcTower[rdo].neighbor[j-1];
350  --j;
351  }
352  energies[j] = indexE;
353  bemcTower[rdo].neighbor[j] = indexId;
354  }
355 
356  //
357  // Calculate energy-weighted centroid of cluster
358  //
359  bemcCluster[rdo].E = (bemcData[rdo] - bemcTower[rdo].pedestal) / bemcTower[rdo].gain;
360  bemcCluster[rdo].x = bemcTower[rdo].x * bemcCluster[rdo].E;
361  bemcCluster[rdo].y = bemcTower[rdo].y * bemcCluster[rdo].E;
362  bemcCluster[rdo].z = bemcTower[rdo].z * bemcCluster[rdo].E;
363 
364  for (int i = 0; i < mNumberOfTowersPerCluster - 1; ++i) {
365  bemcCluster[rdo].E += energies[i];
366  bemcCluster[rdo].x += bemcTower[bemcTower[rdo].neighbor[i]].x * energies[i];
367  bemcCluster[rdo].y += bemcTower[bemcTower[rdo].neighbor[i]].y * energies[i];
368  bemcCluster[rdo].z += bemcTower[bemcTower[rdo].neighbor[i]].z * energies[i];
369  }
370 
371  bemcCluster[rdo].x /= bemcCluster[rdo].E;
372  bemcCluster[rdo].y /= bemcCluster[rdo].E;
373  bemcCluster[rdo].z /= bemcCluster[rdo].E;
374 
375  //
376  // Normalize (x, y, z) vector to unity
377  //
378  float norm = sqrt(bemcCluster[rdo].x * bemcCluster[rdo].x +
379  bemcCluster[rdo].y * bemcCluster[rdo].y +
380  bemcCluster[rdo].z * bemcCluster[rdo].z);
381 
382  if (norm) {
383  bemcCluster[rdo].x /= norm;
384  bemcCluster[rdo].y /= norm;
385  bemcCluster[rdo].z /= norm;
386  }
387 }
388 
389 void L2upsilon2006::createHistograms()
390 {
391  hL0SeedTowers = new L2Histo(100, (char*)"L0 seed towers;softId", 4800);
392  hL2SeedTowers = new L2Histo(101, (char*)"L2 seed towers;softId", 4800);
393  hNumberOfL0Seeds = new L2Histo(102, (char*)";Number of L0 seeds", 10);
394  hNumberOfL2Seeds = new L2Histo(103, (char*) ";Number of L2 seeds", 10);
395  hInvMass = new L2Histo(104, (char*)";m_{ee} [GeV]", 200);
396  hTime = new L2Histo(105, (char*)";processing time [#mus]", 200);
397  hEnergyOfL0Cluster = new L2Histo(106, (char*)";Energy of L0 cluster [GeV]", 200);
398  hEnergyOfL2Cluster = new L2Histo(107, (char*)";Energy of L2 cluster [GeV]", 200);
399  hCosTheta = new L2Histo(108, (char*)";cos(#Theta)", 100);
400  hVertexZ = new L2Histo(109, (char*)";z_{vertex} (cm)", 100);
401  hCtbIndex = new L2Histo(110, (char*)";CTB index", 256);
402  hHighTowers = new L2Histo(111, (char*)";softId", 4800);
403  hL0rate = new L2Histo(112, (char*)"L0;time (sec);rate (Hz)", 1800);
404  hL2rate = new L2Histo(113, (char*)"L2;time (sec);rate (Hz)", 1800);
405 
406  mHistograms.push_back(hL0SeedTowers);
407  mHistograms.push_back(hL2SeedTowers);
408  mHistograms.push_back(hNumberOfL0Seeds);
409  mHistograms.push_back(hNumberOfL2Seeds);
410  mHistograms.push_back(hInvMass);
411  mHistograms.push_back(hTime);
412  mHistograms.push_back(hEnergyOfL0Cluster);
413  mHistograms.push_back(hEnergyOfL2Cluster);
414  mHistograms.push_back(hCosTheta);
415  mHistograms.push_back(hVertexZ);
416  mHistograms.push_back(hCtbIndex);
417  mHistograms.push_back(hHighTowers);
418  mHistograms.push_back(hL0rate);
419  mHistograms.push_back(hL2rate);
420 }
421 
422 void L2upsilon2006::writeHistograms()
423 {
424  char filename[FILENAME_MAX];
425  sprintf(filename, "%s/run%d.l2ups.histo.bin", mOutDir, mRunNumber);
426  FILE* fp = fopen(filename, "w");
427  if (!fp) {
428  fprintf(mLogFile, "L2upsilon2006: WARNING - Can't open histogram file %s\n", filename);
429  return;
430  }
431  for (list<L2Histo*>::iterator i = mHistograms.begin(); i != mHistograms.end(); ++i)
432  (*i)->write(fp);
433  fclose(fp);
434 }
435 
436 void L2upsilon2006::resetHistograms()
437 {
438  for_each(mHistograms.begin(), mHistograms.end(), mem_fun(&L2Histo::reset));
439 }
440 
441 void L2upsilon2006::deleteHistograms()
442 {
443  for (list<L2Histo*>::iterator i = mHistograms.begin(); i != mHistograms.end(); ++i) {
444  delete *i;
445  *i = 0;
446  }
447 }
448 
449 void L2upsilon2006::print()
450 {
451  L2upsilonResult2006* result = (L2upsilonResult2006*)&trgData->TrgSum.L2Result[mResultOffset];
452  fprintf(mLogFile, "-----------------------------------------------\n");
453  fprintf(mLogFile, "Number of L0 seed towers: %d\n", result->numberOfL0SeedTowers);
454  fprintf(mLogFile, "Number of L2 seed towers: %d\n", result->numberOfL2SeedTowers);
455  fprintf(mLogFile, "Invariant mass [GeV]: %f\n", result->invariantMass);
456  fprintf(mLogFile, "Events seen: %d\n", result->eventsSeen);
457  fprintf(mLogFile, "Events accepted: %d\n", result->eventsAccepted);
458  fprintf(mLogFile, "Processing time [us]: %f\n", result->processingTime / CPU_MHZ);
459  fprintf(mLogFile, "Energy of L0 cluster [GeV]: %f\n", result->energyOfL0Cluster);
460  fprintf(mLogFile, "Energy of L2 cluster [GeV]: %f\n", result->energyOfL2Cluster);
461 }
Generated by   doxygen 1.8.5