StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
L2Upsilon2009.cxx
1 #include <stdio.h>
2 #include <string.h>
3 #include <stdlib.h>
4 #include <time.h>
5 #include <math.h>
6 
7 
8 #ifdef IS_REAL_L2 //in l2-ana environment
9  #include "../L2algoUtil/L2EmcDb.h"
10  #include "../L2algoUtil/L2Histo.h"
11 #else //full path needed for cvs'd code
12  #include "StTriggerUtilities/L2Emulator/L2algoUtil/L2EmcDb.h"
13  #include "StTriggerUtilities/L2Emulator/L2algoUtil/L2Histo.h"
14  #include "StTriggerUtilities/L2Emulator/L2algoUtil/L2EmcGeom.h"
15 #endif
16 
17 #include "L2Upsilon2009.h"
18 #include "read_in_geog.h"
19 //extern inline void *memset(void *s, int c, size_t n);
20 
21 
22 //=================================================
23 //=================================================
24 L2Upsilon2009::L2Upsilon2009(const char* name, L2EmcDb* db, L2EmcGeom *geoX, char* outDir, int resOff) : L2VirtualAlgo2009( name, db, outDir, true, false, resOff) {
25  /* called one per days
26  all memory allocation must be done here
27  */
28 
29  mGeom=geoX;
30  if (!mGeom)
31  criticalError("L2Upsilon2009 is broken -- can't find geom.");
32 
33  setMaxHist(32); // set upper range, I uses only 2^N -it is easier to remember
34  createHisto();
35 
36  //------- self-consistency checks, should never fail
37  if (!(sizeof(L2UpsilonResult2009)== L2UpsilonResult2009::mySizeChar))
38  {
39  criticalError("L2Upsilon2009 has failed consistency check 'sizeof(L2UpsilonResult2009)== L2UpsilonResult2009::mySizeChar'");
40  }
41 
42  prescale = 1;
43 }
44 
45 /* ========================================
46  ======================================== */
47 int
48 L2Upsilon2009::initRunUser( int runNo, int *rc_ints, float *rc_floats) {
49 
50  // unpack params from run control GUI
51 
52  par_prescale = rc_ints[0];
53  fMaxDynamicMaskTowers =rc_ints[1];
54  fHowManyEventPerUpdateDynamicMask=rc_ints[2];
55  par_RndAcceptPrescale=rc_ints[3]; //added by Ross.
56 
57  fL0SeedThreshold = rc_floats[0];
58  fL2SeedThreshold = rc_floats[1];
59  fMinL0ClusterEnergy = rc_floats[2];
60  fMinL2ClusterEnergy = rc_floats[3];
61  fMinInvMass = rc_floats[4];
62  fMaxInvMass = rc_floats[5];
63  fMaxCosTheta = rc_floats[6];
64  fHotTowerThreshold=rc_floats[7];
65  fThresholdRatioOfHotTower=rc_floats[8];
66 
67  fHotTowerSeenTimesThreshold=int(fHowManyEventPerUpdateDynamicMask*fThresholdRatioOfHotTower);
68 
69 /* suggested pars
70  par_prescale = 1;
71  fMaxDynamicMaskTowers =25;
72  fHowManyEventPerUpdateDynamicMask=50;
73 
74  fL0SeedThreshold = 4.;
75  fL2SeedThreshold = 2.5;
76  fMinL0ClusterEnergy = 4.5;
77  fMinL2ClusterEnergy = 3.;
78  fMinInvMass = 5;
79  fMaxInvMass = 20;
80  fMaxCosTheta = 0;
81  fHotTowerThreshold=2.5;
82  fThresholdRatioOfHotTower=0.5;
83 
84  fHotTowerSeenTimesThreshold=25;
85 */
86 
87 /*
88 --- temp par values for test -----
89  par_prescale =2;
90  fL0SeedThreshold = 1.;
91  fL2SeedThreshold = 1;
92  fMinL0ClusterEnergy = 1;
93  fMinL2ClusterEnergy = 1.;
94  fMinInvMass = 1;
95  fMaxInvMass = 20;
96  fMaxCosTheta = 1;
97 //--------------------------------
98 */
99 
100  EventSeen=0;
101  event_accept=0;
102  memset(wrkDynamicMask_tower_stat,0,sizeof(wrkDynamicMask_tower_stat));
103  wrkNumberOfMasked=0;
104 
105  // verify consistency of input params
106  int kBad=0;
107  kBad+=0x00001 * (fMaxCosTheta>1.0||fMaxCosTheta<-1.0);
108  kBad+=0x00002 * (fL0SeedThreshold<fL2SeedThreshold);
109  kBad+=0x00004 * (fMinL0ClusterEnergy<fMinL2ClusterEnergy);
110  kBad+=0x00008 * (fMinInvMass>fMaxInvMass);
111  kBad+=0x00010 * (fHowManyEventPerUpdateDynamicMask<fHotTowerSeenTimesThreshold);
112  kBad+=0x00020 * (fMaxDynamicMaskTowers>100);
113  kBad+=0x00040 * (fThresholdRatioOfHotTower>1);
114 
115  if (mLogFile) {
116  fprintf(mLogFile," L2%s upsilon algorithm initRun(R=%d), compiled: %s , %s\n params:\n",getName(),mRunNumber,__DATE__,__TIME__);
117  fprintf(mLogFile," initRun() params checked for consistency, Error flag=%d\n",kBad);
118  fprintf(mLogFile," prescale=%d \n", par_prescale);
119  fprintf(mLogFile," rndAccept prescale=%d \n", par_RndAcceptPrescale);
120  fprintf(mLogFile," L0SeedThreshold= %.2f (GeV)\n",fL0SeedThreshold);
121  fprintf(mLogFile," L2SeedThreshold= %.2f (GeV)\n",fL2SeedThreshold);
122  fprintf(mLogFile," MinL0ClusterEnergy= %.2f (GeV)\n",fMinL0ClusterEnergy);
123  fprintf(mLogFile," MinL2ClusterEnergy= %.2f (GeV)\n",fMinL2ClusterEnergy);
124  fprintf(mLogFile," MinInvMass= %.2f (GeV)\n",fMinInvMass);
125  fprintf(mLogFile," MaxInvMass= %.2f (GeV)\n",fMaxInvMass);
126  fprintf(mLogFile," MaxCosTheta= %.2f \n",fMaxCosTheta);
127  fprintf(mLogFile," MaxDynamicMaskedTowers= %d (should not greater than 100) \n",fMaxDynamicMaskTowers);
128  fprintf(mLogFile," HotTowerThreshold= %f (GeV) \n",fHotTowerThreshold);
129  fprintf(mLogFile," HowManyEventsPerUpdateDynamicMask= %d \n",fHowManyEventPerUpdateDynamicMask);
130  fprintf(mLogFile," fHotTowerSeenTimesThreshold= %d \n",fHotTowerSeenTimesThreshold);
131  fprintf(mLogFile," fThresholdRatioOfHotTower= %f (count_hot / count_check) \n",fThresholdRatioOfHotTower);
132  }
133 
134 
135  // clear content of all histograms & token-dependet memory
136  int i;
137  for (i=0; i<mxHA;i++) if(hA[i])hA[i]->reset();
138  memset(mBtow,0,sizeof(mBtow));
139 
140  if(kBad>0) return -1*kBad;
141  else if(kBad<0) return kBad;
142 
143  for ( int index=0; index<EmcDbIndexMax; index++ )
144  {
145  const L2EmcDb::EmcCDbItem *x = mDb->getByIndex(index);
146  if ( x==0 ) continue;
147  if ( !mDb->isBTOW(x) ) continue;
148  int sec = x->sec - 1;
149  int sub = 8192;
150  sub = x->sub - 'a';
151  int eta = x->eta - 1;
152  int phi = BtowGeom::mxSubs *sec + sub;
153  int tow = BtowGeom::mxEtaBin *phi + eta; // phi- changes faster
154  int rdo = x->rdo;
155  if (tow<0 || tow>mxBtow || rdo<0 || rdo>mxBtow) return -101;
156 
157  mTower2rdo[ tow ] = rdo; // returns rdo channel for given tower
158  mRdo2tower[ rdo ] = tow; // returns tower for given rdo channel
159 
160  int a,b,c,d;
161  sscanf(x->tube,"id%d-%d-%d-%d",&a,&b,&c,&d);
162  rdo2softID[x->rdo]=a;
163 // hA[20]->fill(a,x->rdo);
164  }
165  return 0; //OK
166 
167 }
168 
169 
170 /* ========================================
171  ======================================== */
172 void
173 L2Upsilon2009::computeUser(int token){
174  clearEvent(token);
175 
176  EventSeen++;
177  if(EventSeen%fHowManyEventPerUpdateDynamicMask==0) update_DynamicMask();
178 
179  L2UpsilonEvent2009 *btowEve=mBtow+token;
180  const HitTower1 *hit=mEveStream_btow[token].get_hits();
181  const int hitSize=mEveStream_btow[token].get_hitSize();
182 
183  for(int i=0;i< hitSize;i++,hit++) {
184 // int tower=mRdo2tower[hit->rdo];
185  int tower=rdo2softID[hit->rdo];
186  wrkBtow_tower[i]=tower;
187  wrkBtow_ene[tower]=hit->ene; // fill energy array (all towers)
188  if(wrkBtow_ene[tower] > fHotTowerThreshold) wrkDynamicMask_tower_stat[tower]++; // wrkDynamicMask_tower_stat is needed to judge hot towers
189  hA[2]->fill(int(wrkBtow_ene[tower]));
190  }
191  hA[1]->fill(hitSize);
192 
193 // Dynamic Masking
194  for(int j=0;j<wrkNumberOfMasked;j++)
195  wrkBtow_ene[wrkDynamicMasked_tower[j]]=0;
196 
197 // ---- store all L0 L2 candidates
198  wrkNumberOfL0=0;
199  wrkNumberOfL2=0;
200  for(int i=0;i< hitSize;i++) {
201  int tower=wrkBtow_tower[i];
202  if(wrkBtow_ene[tower] > fL2SeedThreshold){
203  float clusterE=wrkBtow_ene[tower];
204  for(int k=0;k<geog_Nneighbors[tower];k++){
205  float Ek=wrkBtow_ene[geog_neighbor[tower][k]];
206  int Ecount=0;
207  for(int m=0;m<geog_Nneighbors[tower];m++){
208  if(m==k) continue;
209  if(Ek<wrkBtow_ene[geog_neighbor[tower][m]])Ecount++;
210  }
211  if(Ecount<2)clusterE+=Ek; // top 2 cluster
212  }
213 
214  if(clusterE>fMinL2ClusterEnergy) { // tower seed good && cluster good for L2 threshold
215  wrkL2_seed_tower[wrkNumberOfL2]=tower;
216  wrkL2_seed_ClusterE[wrkNumberOfL2]=clusterE;
217  hA[6]->fill(int(wrkBtow_ene[tower]));
218  hA[8]->fill(int(clusterE));
219  wrkNumberOfL2++;
220  if(wrkBtow_ene[tower] > fL0SeedThreshold&&clusterE>fMinL0ClusterEnergy) { // tower seed good && cluster good for L0 threshold
221  wrkL0_seed_tower[wrkNumberOfL0]=tower;
222  wrkL0_seed_ClusterE[wrkNumberOfL0]=clusterE;
223  hA[5]->fill(int(wrkBtow_ene[tower]));
224  hA[7]->fill(int(clusterE));
225  wrkNumberOfL0++;
226  }
227  }
228  }
229  }
230  hA[3]->fill(wrkNumberOfL2,wrkNumberOfL0);
231 
232 // ------ pair L0 L2, calculate invMass
233  for(int i=0;i<wrkNumberOfL0;i++){
234  int idL2,idL0;
235  float cosTheta,invMass;
236  idL0= wrkL0_seed_tower[i];
237  for(int j=0;j<wrkNumberOfL2;j++){
238  idL2= wrkL2_seed_tower[j];
239  if(idL0==idL2)continue;
240  cosTheta=geog_x[idL0]*geog_x[idL2]+geog_y[idL0]*geog_y[idL2]+geog_z[idL0]*geog_z[idL2];
241  if(cosTheta>fMaxCosTheta) continue;
242  invMass = sqrt(2 *wrkL0_seed_ClusterE[i] * wrkL2_seed_ClusterE[j] * (1 - cosTheta));
243  if(fMinInvMass>invMass || invMass>fMaxInvMass) continue;
244  hA[10]->fill(int(invMass));
245  // accept !!!
246  btowEve->resultBlob.L0SeedTowerID=idL0;
247  btowEve->resultBlob.L2SeedTowerID=idL2;
248  if(wrkL0_seed_ClusterE[i]>25.5)wrkL0_seed_ClusterE[i]=25.5;
249  if(wrkL0_seed_ClusterE[j]>25.5)wrkL0_seed_ClusterE[j]=25.5;
250  if(invMass>25.5)invMass=25.5;
251  btowEve->resultBlob.energyOfL0Cluster=(unsigned char)(wrkL0_seed_ClusterE[i]*10); // 256/25.6=10
252  btowEve->resultBlob.energyOfL2Cluster=(unsigned char)(wrkL2_seed_ClusterE[j]*10); // 256/25.6=10
253  btowEve->resultBlob.invMass=(unsigned char)(invMass*10); // 256/25.6=10
254  btowEve->resultBlob.trigger=true;
255 // btowEve->resultBlob.numberOfL0SeedTowers=wrkNumberOfL0;
256 // btowEve->resultBlob.numberOfL2SeedTowers=wrkNumberOfL2;
257  return;
258  }
259  }
261  btowEve->resultBlob.L0SeedTowerID=5555;
262  btowEve->resultBlob.L2SeedTowerID=5555;
263  btowEve->resultBlob.energyOfL0Cluster=0;
264  btowEve->resultBlob.energyOfL2Cluster=0;
265  btowEve->resultBlob.invMass=0;
266  btowEve->resultBlob.trigger=false;
267 // btowEve->resultBlob.numberOfL0SeedTowers=wrkNumberOfL0;
268 // btowEve->resultBlob.numberOfL2SeedTowers=wrkNumberOfL2;
269  return;
270 }
271 
272 
273 /* ========================================
274  ======================================== */
275 bool
276 L2Upsilon2009::decisionUser(int token, int *myL2Result){
277  // INPUT: token + comput() results stored internally
278  // OUTPUT: yes/no + pointer to L2Result
279 
280  // get pointers to internal private event storage
281  L2UpsilonEvent2009 *btowEve=mBtow+token;
282 // (*myL2Result)=&(btowEve->resultBlob);
283  memcpy(myL2Result,&(btowEve->resultBlob),sizeof(L2UpsilonResult2009));
284  if (par_prescale > 1) {
285  prescale++;
286  prescale%=par_prescale;
287  if(prescale!=0) btowEve->resultBlob.trigger=false;
288  }
289  btowEve->isFresh++; // mark the data as stale
290  hA[13]->fill(1);
291 
292  if(btowEve->resultBlob.trigger){
293  event_accept++;
294  hA[13]->fill(2);
295  hA[14]->fill(btowEve->resultBlob.L0SeedTowerID);
296  hA[15]->fill(btowEve->resultBlob.L2SeedTowerID);
297  hA[16]->fill(int((btowEve->resultBlob.energyOfL0Cluster)*0.1));
298  hA[17]->fill(int((btowEve->resultBlob.energyOfL2Cluster)*0.1));
299  hA[18]->fill(int((btowEve->resultBlob.invMass)*0.1));
300  }
301  else hA[13]->fill(3);
302 
303  return btowEve->resultBlob.trigger;
304 }
305 
306 
307 /* ========================================
308  ======================================== */
309 void
310 L2Upsilon2009::finishRunUser() { /* called once at the end of the run */
311  // do whatever you want, log-file & histo-file are still open
312 
313  if (mLogFile){
314  fprintf(mLogFile,"finishRunUser - %s \n",getName());
315  fprintf(mLogFile,"EventSeen = %d \n",EventSeen);
316  fprintf(mLogFile,"event accept = %d \n",event_accept);
317  }
318 
319 }
320 
321 
322 //=======================================
323 //=======================================
324 void
325 L2Upsilon2009::createHisto() {
326  setMaxHist(32); // PMN added - histogram count does not seem to be initialiazed anywere.
327  //memset(hA,0,sizeof(hA));
328 
329  hA[1]=new L2Histo(1,"Upsilon:# of hits (no L0 required)", 200);
330  hA[2]=new L2Histo(2,"Upsilon:energy for all hits (GeV) (no L0 required)", 20);
331  hA[3]=new L2Histo(3,"Upsilon: # of L0 candidates vs. # of L2 candidates (no L0 required)", 50,50);
332 
333  hA[5]=new L2Histo(5,"Upsilon: L0 seeds Energy (GeV) (no L0 required)", 20);
334  hA[6]=new L2Histo(6,"Upsilon: L2 seeds Energy (GeV) (no L0 required)", 20);
335 
336  hA[7]=new L2Histo(7,"Upsilon: L0 cluster Energy (GeV) (no L0 required)", 25);
337  hA[8]=new L2Histo(8,"Upsilon: L2 cluster Energy (GeV) (no L0 required)", 25);
338 
339  hA[10]= new L2Histo(10,"Upsilon: InvMass of accepted pair (GeV) (no L0 required)", 25);
340  hA[11]=new L2Histo(11,"Upsilon: # of masked towers (no L0 required)", 100);
341  hA[12]=new L2Histo(12,"Upsilon: masked tower softid (no L0 required) ", 4801);
342 
343  hA[13] = new L2Histo(13,"Upsilon: L0 counter (1=L0 fired, 2=L0 fired+L2 accepted, 3=L0 fired+L2 NOT accepted)",5);
344  hA[14] = new L2Histo(14,"Upsilon: L0 seed tower id (L0 fired+L2 accepted)",4801);
345  hA[15] = new L2Histo(15,"Upsilon: L2 seed tower id (L0 fired+L2 accepted)",4801);
346  hA[16]= new L2Histo(16,"Upsilon: L0 cluster Energy (GeV) (L0 fired+L2 accepted)", 25);
347  hA[17]= new L2Histo(17,"Upsilon: L2 cluster Energy (GeV) (L0 fired+L2 accepted)", 25);
348  hA[18]= new L2Histo(18,"Upsilon: InvMass of accepted pair (GeV) (L0 fired+L2 accepted)", 25);
349 
350 // hA[20]=new L2Histo(20,"id",4801,4801);
351 
352 }
353 
354 //=======================================
355 //=======================================
356 void
357 L2Upsilon2009::clearEvent(int token){
358  memset(wrkBtow_ene,0,sizeof(wrkBtow_ene));
359  memset(wrkBtow_tower,0,sizeof(wrkBtow_tower));
360  memset(wrkL2_seed_tower,0,sizeof(wrkL2_seed_tower));
361  memset(wrkL2_seed_ClusterE,0,sizeof(wrkL2_seed_ClusterE));
362  memset(wrkL0_seed_tower,0,sizeof(wrkL0_seed_tower));
363  memset(wrkL0_seed_ClusterE,0,sizeof(wrkL0_seed_ClusterE));
364  memset(&(mBtow[token].resultBlob),0, sizeof(L2UpsilonResult2009));
365  wrkNumberOfL2=wrkNumberOfL0=0;
366 }
367 
368 void
369 L2Upsilon2009::update_DynamicMask(){
370  int count=0;
371  for(int i=0;i<mxBtow;i++){
372  if(wrkDynamicMask_tower_stat[i]>fHotTowerSeenTimesThreshold && count<fMaxDynamicMaskTowers){ // tower i is hot
373  wrkDynamicMasked_tower[count]=i;
374  hA[12]->fill(i);
375  count++;
376  }
377  }
378  memset(wrkDynamicMask_tower_stat,0,sizeof(wrkDynamicMask_tower_stat));
379  wrkNumberOfMasked=count;
380  hA[11]->fill(count);
381 }
Generated by   doxygen 1.8.5