StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StFcsTriggerSimMaker.cxx
1 // \class StFmsTriggerSimMaker
2 // \author Akio Ogawa
3 //
4 // $Id: StFcsTriggerSimMaker.cxx,v 1.2 2021/05/30 21:40:56 akio Exp $
5 // $Log: StFcsTriggerSimMaker.cxx,v $
6 // Revision 1.2 2021/05/30 21:40:56 akio
7 // Many updates from trigger commissioning on Run21 OO data
8 //
9 // Revision 1.1 2021/03/30 13:33:53 akio
10 // Moved from $CVSROOT/offline/upgrade/akio/ to $CVSROOT/StRoot/StSpinPool/
11 //
12 // Revision 1.7 2021/02/25 21:56:10 akio
13 // Int_t -> int
14 //
15 // Revision 1.6 2020/07/24 17:22:39 akio
16 // adding option to reading in EPD masks
17 //
18 // Revision 1.5 2020/06/01 20:33:42 akio
19 // adapt for DAQ_FCS change
20 //
21 // Revision 1.4 2020/05/29 18:55:47 akio
22 // Modiying to make it run with Tonko's wrapper
23 //
24 // Revision 1.3 2019/06/26 18:01:06 akio
25 // assuming first timebin from MC has ADC
26 //
27 // Revision 1.2 2019/05/17 15:58:56 akio
28 // updates
29 //
30 // Revision 1.1 2018/11/12 13:15:58 akio
31 // Initial version
32 //
33 
34 #include "StFcsTriggerSimMaker.h"
35 
36 #include "TTree.h"
37 #include "TFile.h"
38 
39 #include "StMessMgr.h"
40 #include "Stypes.h"
41 #include "StarGenerator/BASE/StarPrimaryMaker.h"
42 #include "StarGenerator/EVENT/StarGenEvent.h"
43 
44 #include "StThreeVectorF.hh"
45 #include "StEvent/StEventTypes.h"
46 #include "StEvent/StFcsHit.h"
47 #include "StFcsDbMaker/StFcsDb.h"
48 
49 #include "RTS/include/rtsLog.h"
50 
51 #include "StRoot/RTS/src/TRG_FCS/fcs_trg_base.h"
52 
53 #include "StMaker.h"
54 #include "StChain.h"
55 #include "StMuDSTMaker/COMMON/StMuDst.h"
56 #include "StMuDSTMaker/COMMON/StMuEvent.h"
57 #include "StMuDSTMaker/COMMON/StMuFcsCollection.h"
58 #include "StMuDSTMaker/COMMON/StMuFcsHit.h"
59 #include "StMuDSTMaker/COMMON/StMuTypes.hh"
60 #include "StMuDSTMaker/COMMON/StMuEvent.h"
61 #include "StMuDSTMaker/COMMON/StMuFcsCollection.h"
62 
63 namespace {
64  enum {kMaxNS=2, kMaxDet=3, kMaxDep=24, kMaxCh=32, kMaxEcalDep=24, kMaxHcalDep=8, kMaxPresDep=4, kMaxLink2=2};
65  uint32_t fcs_trg_sim_adc[kMaxNS][kMaxDet][kMaxDep][kMaxCh] ;
66  float fcs_trg_pt_correction[kMaxNS][kMaxDet][kMaxDep][kMaxCh];
67  float fcs_trg_gain_correction[kMaxNS][kMaxDet][kMaxDep][kMaxCh];
68  uint16_t fcs_trg_pedestal[kMaxNS][kMaxDet][kMaxDep][kMaxCh] ;
69 
70  static const int mNTRG=21;
71  static const char* ctrg[mNTRG]={"JP2", "JPA1", "JPA0", "JPBC1", "JPBC0", "JPDE1", "JPDE0",
72  "DiJP", "DiJPAsy",
73  "DY","JPsi","DYNoEpd","DYAsy",
74  "Had2","Had1","Had0",
75  "EM2","EM1","EM0",
76  "ELE2","EM3"};
77  int NTRG[mNTRG+1];
78 }
79 
80 ClassImp(StFcsTriggerSimMaker);
81 
82 StFcsTriggerSimMaker::StFcsTriggerSimMaker(const char* name): StMaker(name) {}
83 
84 StFcsTriggerSimMaker::~StFcsTriggerSimMaker(){}
85 
86 int StFcsTriggerSimMaker::Init(){
87  LOG_INFO << "StFcsTriggerSimMaker::Init" << endm;
88 
89  mFcsDb=static_cast<StFcsDb*>(GetDataSet("fcsDb"));
90  if(!mFcsDb){
91  LOG_ERROR << "StFcsTriggerSimMaker::Init Failed to get StFcsDb" << endm;
92  return kStFatal;
93  }
94 
95  rtsLogOutput(RTS_LOG_STDERR) ;
96 
97  mTrgSim = new fcs_trg_base();
98  mTrgSim->sim_mode=1;
99  mTrgSim->init(".");
100  mTrgSim->run_start(0);
101  mTrgSim->fcs_trgDebug=mDebug;
102 
103  //trigegr versions
104  if(mTrgSelect==201900){
105  mTrgSim->stage_version[0]=0;
106  mTrgSim->stage_version[1]=0;
107  mTrgSim->stage_version[2]=0;
108  mTrgSim->stage_version[3]=0;
109  }else if(mTrgSelect==202201){
110  mTrgSim->stage_version[0]=0;
111  mTrgSim->stage_version[1]=1;
112  mTrgSim->stage_version[2]=1;
113  mTrgSim->stage_version[3]=1;
114  }else if(mTrgSelect==202203){
115  mTrgSim->stage_version[0]=2;
116  mTrgSim->stage_version[1]=1;
117  mTrgSim->stage_version[2]=3;
118  mTrgSim->stage_version[3]=3;
119  }else if(mTrgSelect==202204){
120  mTrgSim->stage_version[0]=2;
121  mTrgSim->stage_version[1]=1;
122  mTrgSim->stage_version[2]=4;
123  mTrgSim->stage_version[3]=3;
124  }else if(mTrgSelect==202205){
125  mTrgSim->stage_version[0]=2;
126  mTrgSim->stage_version[1]=1;
127  mTrgSim->stage_version[2]=5;
128  mTrgSim->stage_version[3]=3;
129  }else if(mTrgSelect==202206){
130  mTrgSim->stage_version[0]=2;
131  mTrgSim->stage_version[1]=1;
132  mTrgSim->stage_version[2]=6;
133  mTrgSim->stage_version[3]=3;
134  }else if(mTrgSelect==202207){
135  mTrgSim->stage_version[0]=2;
136  mTrgSim->stage_version[1]=1;
137  mTrgSim->stage_version[2]=7;
138  mTrgSim->stage_version[3]=7;
139  }else if(mTrgSelect==202209){
140  mTrgSim->stage_version[0]=3;
141  mTrgSim->stage_version[1]=1;
142  mTrgSim->stage_version[2]=7;
143  mTrgSim->stage_version[3]=7;
144  }
145 
146  //Thresholds
147  if(mThresholdFile){
148  readThresholdFile();
149  }else if(mThresholdDb){
150  readThresholdDb();
151  }
152  //mTrgSim->EM_HERATIO_THR = 32;
153  //mTrgSim->HAD_HERATIO_THR = 32;
154 
155  //Trigger Id names
156  readTrgId();
157 
158  //EPD mask
159  if(mPresMask){
160  printf("Reading PresMask from %s\n",mPresMask);
161  FILE* F=fopen(mPresMask,"r");
162  if(F==NULL){
163  printf("Cannot open %s\n",mPresMask);
164  }else{
165  char line[512];
166  int r,c,m[6];
167  while(fgets(line,sizeof(line),F)){
168  if(line[0]=='#') {
169  printf("%s",line);
170  continue;
171  }
172  sscanf(line,"%d %d %x %x %x %x %x %x",&r,&c,&m[0],&m[1],&m[2],&m[3],&m[4],&m[5]);
173  printf("%2d %1d %08x %08x %08x %08x %08x %08x\n",r,c,m[0],m[1],m[2],m[3],m[4],m[5]);
174  for(int i=0; i<6; i++) mTrgSim->PRES_MASK[r-1][c-1][i]=m[i];
175  }
176  mTrgSim->fcs_readPresMaskFromText=1;
177  }
178  }
179 
180  memset(NTRG,0,sizeof(NTRG));
181  return kStOK;
182 }
183 
184 int StFcsTriggerSimMaker::InitRun(int runNumber){
185  LOG_INFO << "StFcsTriggerSimMaker::InitRun" << endm;
186  //print out 4x4 and JP info
187  if(mDebug>0){
188  print4B4();
189  printJP();
190  }
191 
192  //QA root file
193  if(mQaTreeFilename){
194  mQaTreeFile=new TFile(mQaTreeFilename,"RECREATE");
195  mTree = new TTree("trgsim","trigger sim QA");
196  mTree->Branch("flt",&mFlt,"flt/I");
197  mTree->Branch("trg",&mTrg,"trg/I");
198  }
199  if(mQaHistFilename){
200  mQaHistFile=new TFile(mQaHistFilename,"RECREATE");
201  mTrgRate = new TH1F("FcsTrgRate","FcsTrgRate",mNTRG+1,0,mNTRG+1);
202  }
203 
204  //Write Text event file & gainfile
205  FILE* gainfile=0;
206  FILE* gainfile2=0;
207  if(mFilename){
208  mFile = fopen(mFilename,"w");
209  gainfile=fopen("fcs_et_gain.txt","w");
210  gainfile2=fopen("fcs_et_gain2.txt","w");
211  }
212 
213  //Fill ETgain, GainCorr and Pedestal
214  for(int det=0; det<=kFcsNDet; det++){
215  int nid=mFcsDb->maxId(det);
216  for(int id=0; id<nid; id++){
217  int ehp,ns,crt,sub,dep,ch;
218  mFcsDb->getDepfromId(det,id,ehp,ns,crt,sub,dep,ch);
219  if(det<4){
220  fcs_trg_pt_correction[ns][ehp][dep][ch] = mFcsDb->getEtGain(det,id,mEtFactor);
221  fcs_trg_gain_correction[ns][ehp][dep][ch] = mFcsDb->getGainOnline(det,id);
222  }else{
223  fcs_trg_pt_correction[ns][ehp][dep][ch] = 1.0;
224  fcs_trg_gain_correction[ns][ehp][dep][ch] = 1.0;
225  }
226  fcs_trg_pedestal[ns][ehp][dep][ch] = 0;
227 
228  mTrgSim->p_g[ns][ehp][dep][ch].ped = fcs_trg_pedestal[ns][ehp][dep][ch];
229 
230  if(mOverwriteGain==1){
231  float ggg = fcs_trg_pt_correction[ns][ehp][dep][ch];
232  //float ggg = (fcs_trg_pt_correction[ns][ehp][dep][ch]-1.0)/2.0 + 1.0;
233  float gg = ggg * fcs_trg_gain_correction[ns][ehp][dep][ch];
234  int g = (uint32_t)(gg*256.0+0.5) ;
235  mTrgSim->p_g[ns][ehp][dep][ch].gain = g;
236  }
237 
238  /*
239  printf("AAAGAIN %1d %1d %2d %2d pT=%6.3f corr=%6.3f ped=%4d\n",ns,ehp,dep,ch,
240  fcs_trg_pt_correction[ns][ehp][dep][ch],
241  fcs_trg_gain_correction[ns][ehp][dep][ch],
242  fcs_trg_pedestal[ns][ehp][dep][ch]);
243  */
244 
245  if(gainfile)
246  fprintf(gainfile,"%2d %2d %2d %2d %8.3f\n",ns,ehp,dep,ch,
247  fcs_trg_pt_correction[ns][ehp][dep][ch]);
248  if(gainfile2)
249  fprintf(gainfile2,"%2d %2d %2d %2d %8.3f\n",ns,ehp,dep,ch,
250  (fcs_trg_pt_correction[ns][ehp][dep][ch]-1.0)/2.0 + 1.0);
251  }
252  }
253  if(gainfile) fclose(gainfile);
254  if(gainfile2) fclose(gainfile2);
255  return kStOK;
256 }
257 
258 int StFcsTriggerSimMaker::Finish(){
259  mTrgSim->run_stop();
260  if(mFile) {
261  printf("Closing %s\n",mFilename);
262  fclose(mFile);
263  }
264  if(mQaTreeFile){
265  printf("Closing %s\n",mQaTreeFilename);
266  mTree->Write();
267  mQaTreeFile->Close();
268  }
269  if(mQaHistFile){
270  printf("Closing %s\n",mQaHistFilename);
271  mTrgRate->Write();
272  mQaHistFile->Close();
273  }
274 
275  int tot = NTRG[mNTRG];
276  LOG_INFO << "Triggers counts/"<<tot<<" (rate at 5MHz BBC)"<<endm;
277  for(int i=0; i<mNTRG; i++)
278  LOG_INFO << Form("%8s %9d (%12.2f)",ctrg[i],NTRG[i],double(NTRG[i])/tot*5.0e6)<<endm;
279 
280  return kStOK;
281 }
282 
283 int StFcsTriggerSimMaker::Make(){
284  StEvent* event = nullptr;
285  event = (StEvent*)GetInputDS("StEvent");
286  if(!event){
287  LOG_INFO << "StFcsTriggerSimMaker::Make did not find StEvent"<<endm;}
288  else {
289  mFcsColl = event->fcsCollection();
290  if(!mFcsColl){LOG_INFO << "StFcsTriggerSimMaker::Make did not find StEvent->StFcsCollection"<<endm;}
291  }
292 
293  StMuEvent* muevent = nullptr;
294  if((!event)||(!mFcsColl)){
295  LOG_INFO<<"No StEvent info available for StFcsTriggerSimMaker. "<< endm;
296  muevent = StMuDst::event();
297  mMuFcsColl = StMuDst::muFcsCollection();
298  if(muevent && mMuFcsColl){
299  LOG_INFO <<"Proceeding with StMuDst info to be used with StfcsTriggerSimMaker."<< endm;
300  }else{
301  LOG_ERROR << "StFcsTriggerSimMaker::Make did not find StEvent and MuEvent." << endm;
302  return kStErr;
303  }
304  }
305 
306  mTrgSim->start_event();
307 
308  //Feed ADC
309  static uint16_t data[8];
310  memset(data,0,sizeof(data)) ;
311  memset(fcs_trg_sim_adc,0,sizeof(fcs_trg_sim_adc));
312  int n=0;
313  for(int det=0; det<=kFcsNDet; det++){
314 
315  int ns = mFcsDb->northSouth(det);
316  int ehp = mFcsDb->ecalHcalPres(det);
317 
318  if((event)&&(mFcsColl)){
319  StSPtrVecFcsHit& hits = mFcsColl->hits(det);
320 
321  int nh = mFcsColl->numberOfHits(det);
322 
323  for(int i=0; i<nh; i++){
324  StFcsHit* hit=hits[i];
325  unsigned short ch = hit->channel();
326 
327  if(ehp<0 || ch>=32) continue;
328  feedADC(hit, ns, ehp, data);
329  n++;
330  }
331 
332  }else if((muevent)&&(mMuFcsColl)){ //Use StMuDst info instead of StEvent for Trigger Reconstruction
333 
334  int nh = mMuFcsColl->numberOfHits(det);
335  int det_hit_index = mMuFcsColl->indexOfFirstHit(det);
336 
337  for(int i=0; i<nh; i++){
338 
339  int hit_index = i + det_hit_index;
340  StMuFcsHit* hit = mMuFcsColl->getHit(hit_index);
341  unsigned short ch = hit->channel();
342 
343  if(ehp<0 || ch>=32) continue;
344  feedADC(hit, ns, ehp, data);
345  n++;
346  }
347  }
348  }
349  if(mFile) fprintf(mFile,"%2d %2d %2d %2d %5d\n",-1,0,0,0,0);
350  LOG_INFO << Form("StFcsTriggerSimMaker feeded %d hits",n) << endm;;
351 
352  //Run Trigger Simulation
353  // uint16_t dsm_out = fcs_trg_run(mTrgSelect, mDebug);
354  uint32_t dsm_out = mTrgSim->end_event();
355  LOG_INFO << Form("AAA dsmout=%08x",dsm_out)<<endm;
356 
357  //QA Tree
358  mFlt=0;
359  StarPrimaryMaker* pmkr= static_cast<StarPrimaryMaker*>(GetMaker("PrimaryMaker"));
360  if(pmkr){
361  StarGenEvent *ge = pmkr->event();
362  if(ge){
363  mFlt=ge->GetFilterResult();
364  }
365  }
366  mTcu=dsm_out;
367  if(mQaTreeFile) mTree->Fill();
368 
369  //Results
370  LOG_INFO << Form("Output to TCU = 0x%08x Filter=0x%08x",mTcu,mFlt)<<endm;
371 
372  //TCU Trigger Emulation
373  int trg[mNTRG];
374  mTrg=0;
375  memset(trg,0,sizeof(trg));
376  if((dsm_out >> 6) & 0x1) {trg[ 0]=1; mTrg+=1<<0;} //JP2
377  if((dsm_out >> 7) & 0x1) {trg[ 1]=1; mTrg+=1<<1;} //JPA1
378  if((dsm_out >>10) & 0x1) {trg[ 2]=1; mTrg+=1<<2;} //JPA0
379  if((dsm_out >> 8) & 0x1) {trg[ 3]=1; mTrg+=1<<3;} //JPBC1
380  if((dsm_out >>11) & 0x1) {trg[ 4]=1; mTrg+=1<<4;} //JPBC0
381  if((dsm_out >> 9) & 0x1) {trg[ 5]=1; mTrg+=1<<5;} //JPDE1
382  if((dsm_out >>12) & 0x1) {trg[ 6]=1; mTrg+=1<<6;} //JPDE0
383  if((dsm_out >>13) & 0x1) {trg[ 7]=1; mTrg+=1<<7;} //DiJP
384  if((dsm_out >>14) & 0x1) {trg[ 8]=1; mTrg+=1<<8;} //DiJpAsy
385  if( ((dsm_out >>18) & 0x1) &&
386  ((dsm_out >>26) & 0x1) ) {trg[ 9]=1; mTrg+=1<<9;} //DY
387  if( ((dsm_out >>17) & 0x1) &&
388  ((dsm_out >>25) & 0x1) ) {trg[10]=1; mTrg+=1<<10;} //Jpsi
389  if( ((dsm_out >>19) & 0x1) &&
390  ((dsm_out >>27) & 0x1) ) {trg[11]=1; mTrg+=1<<11;} //DYNoEpd
391  if((dsm_out >>15) & 0x1) {trg[12]=1; mTrg+=1<<12;} //DYAsy
392  if((dsm_out >> 2) & 0x1) {trg[13]=1; mTrg+=1<<13;} //Had2
393  if((dsm_out >> 1) & 0x1) {trg[14]=1; mTrg+=1<<14;} //Had1
394  if((dsm_out >> 0) & 0x1) {trg[15]=1; mTrg+=1<<15;} //Had0
395  if((dsm_out >> 5) & 0x1) {trg[16]=1; mTrg+=1<<16;} //EM2
396  if((dsm_out >> 4) & 0x1) {trg[17]=1; mTrg+=1<<17;} //EM1
397  if((dsm_out >> 3) & 0x1) {trg[18]=1; mTrg+=1<<18;} //EM0
398  if( ((dsm_out >>18) & 0x1) ||
399  ((dsm_out >>26) & 0x1) ) {trg[19]=1; mTrg+=1<<19;} //ELE2
400  if( ((dsm_out >>19) & 0x1) ||
401  ((dsm_out >>27) & 0x1) ) {trg[20]=1; mTrg+=1<<20;} //EM3
402 
403  if(mTrgRate) mTrgRate->Fill(mNTRG);
404  NTRG[mNTRG]++;
405  for(int i=0; i<mNTRG; i++){
406  if(trg[i]) {
407  if(mTrgRate) mTrgRate->Fill(i);
408  NTRG[i]++;
409  }
410  }
411 
412  LOG_INFO << "Triggers = ";
413  for(int i=0; i<mNTRG; i++){ if(trg[i]) LOG_INFO << ctrg[i] << " ";}
414  LOG_INFO << endm;
415 
416  return kStOK;
417 }
418 
419 void StFcsTriggerSimMaker::runStage2(link_t ecal[], link_t hcal[], link_t pres[], geom_t& geo, link_t output[], unsigned short& dsm,
420  int dta[], int dsmout, int sim[], int simdsmout, int iev){
421  unsigned short emu[8];
422  mTrgSim->stage_2(ecal,hcal,pres,geo,output,&dsm);
423  for(int i=0; i<8; i++) emu[i] = output[0].d[i] + (output[1].d[i] << 8);
424  printf("Event#=%3d emuout = ",iev); for(int i=0; i<8; i++) {printf("%04x ",emu[i]);}
425  printf("TCU=%04x\n",dsm);
426  printf("Event#=%3d dtaout = ",iev); for(int i=0; i<8; i++) {printf("%04x ",dta[i]);}
427  printf("TCU=%04x\n",dsmout);
428  printf("Event#=%3d simout = ",iev); for(int i=0; i<8; i++) {printf("%04x ",sim[i]);}
429  printf("TCU=%04x\n",simdsmout);
430  const char* s2bit[3][16]={{"EM0 ","EM1 ","EM2 ","EM3 ","ELE0","ELE1","ELE2","PRS ",
431  "HAD0","HAD1","HAD2","xxxx","EHT ","HHT ","ETOT","HTOT"},
432  {"JPA2","JPB2","JPC2","JPD2","JPE2","xxxx","xxxx","xxxx",
433  "JPA1","JPB1","JPC1","JPD1","JPE1","xxxx","xxxx","xxxx"},
434  {"JPA0","JPB0","JPC0","JPD0","JPE0","xxxx","xxxx","xxxx",
435  "JPAd","JPBd","JPCd","JPDd","JPEd","xxxx","xxxx","xxxx"}};
436  for(int i=0; i<3; i++){
437  for(int j=0; j<16; j++){
438  if( ((dta[i]>>j)&1) != ((sim[i]>>j)&1) ||
439  ((sim[i]>>j)&1) != ((emu[i]>>j)&1) ||
440  ((emu[i]>>j)&1) != ((dta[i]>>j)&1) ){
441  printf("Event#=%3d STG2to3 ns=%1d i=%d j=%2d %s dat=%x sim=%x emu=%x",
442  iev,geo.ns,i,j,s2bit[i][j],(dta[i]>>j)&1,(sim[i]>>j)&1,(emu[i]>>j)&1);
443  if(i==0 && j<7){
444  int maxr=0, maxc=0, max=0;
445  for(int r=0; r<15; r++){
446  for(int c=0; c<9; c++){
447  if(max < mTrgSim->esum[geo.ns][r][c]) {maxr=r; maxc=c; max=mTrgSim->esum[geo.ns][r][c];}
448  }
449  }
450  printf(" EsumMax=%4d ratio=%4.3f",
451  mTrgSim->esum[geo.ns][maxr][maxc],
452  mTrgSim->ratiomax[geo.ns][maxr][maxc]);
453  if(j>=4) printf(" Epd=%1d",mTrgSim->epdcoin[geo.ns][maxr][maxc]);
454  if(j==0) printf(" EMTHR0=%4d",mTrgSim->EMTHR0);
455  if(j==1) printf(" EMTHR1=%4d",mTrgSim->EMTHR1);
456  if(j==2) printf(" EMTHR2=%4d",mTrgSim->EMTHR2);
457  if(j==3) printf(" ELETHR2=%4d",mTrgSim->ELETHR2);
458  if(j==4) printf(" ELETHR0=%4d",mTrgSim->ELETHR0);
459  if(j==5) printf(" ELETHR1=%4d",mTrgSim->ELETHR1);
460  if(j==6) printf(" ELETHR2=%4d",mTrgSim->ELETHR2);
461  }
462  if(i==1 && j== 0){ printf(" JPA=%4d thr2=%4d",mTrgSim->jet[geo.ns][0],mTrgSim->JPATHR2);}
463  if(i==1 && j== 1){ printf(" JPB=%4d thr2=%4d",mTrgSim->jet[geo.ns][1],mTrgSim->JPBCTHR2);}
464  if(i==1 && j== 2){ printf(" JPC=%4d thr2=%4d",mTrgSim->jet[geo.ns][2],mTrgSim->JPBCTHR2);}
465  if(i==1 && j== 3){ printf(" JPD=%4d thr2=%4d",mTrgSim->jet[geo.ns][3],mTrgSim->JPDETHR2);}
466  if(i==1 && j== 4){ printf(" JPE=%4d thr2=%4d",mTrgSim->jet[geo.ns][4],mTrgSim->JPDETHR2);}
467  if(i==1 && j== 8){ printf(" JPA=%4d thr1=%4d",mTrgSim->jet[geo.ns][0],mTrgSim->JPATHR1);}
468  if(i==1 && j== 9){ printf(" JPB=%4d thr1=%4d",mTrgSim->jet[geo.ns][1],mTrgSim->JPBCTHR1);}
469  if(i==1 && j==10){ printf(" JPC=%4d thr1=%4d",mTrgSim->jet[geo.ns][2],mTrgSim->JPBCTHR1);}
470  if(i==1 && j==11){ printf(" JPD=%4d thr1=%4d",mTrgSim->jet[geo.ns][3],mTrgSim->JPDETHR1);}
471  if(i==1 && j==12){ printf(" JPE=%4d thr1=%4d",mTrgSim->jet[geo.ns][4],mTrgSim->JPDETHR1);}
472  if(i==2 && j== 0){ printf(" JPA=%4d thr0=%4d",mTrgSim->jet[geo.ns][0],mTrgSim->JPATHR0);}
473  if(i==2 && j== 1){ printf(" JPB=%4d thr0=%4d",mTrgSim->jet[geo.ns][1],mTrgSim->JPBCTHR0);}
474  if(i==2 && j== 2){ printf(" JPC=%4d thr0=%4d",mTrgSim->jet[geo.ns][2],mTrgSim->JPBCTHR0);}
475  if(i==2 && j== 3){ printf(" JPD=%4d thr0=%4d",mTrgSim->jet[geo.ns][3],mTrgSim->JPDETHR0);}
476  if(i==2 && j== 4){ printf(" JPE=%4d thr0=%4d",mTrgSim->jet[geo.ns][4],mTrgSim->JPDETHR0);}
477  if(i==2 && j== 8){ printf(" JPA=%4d thrD=%4d",mTrgSim->jet[geo.ns][0],-1);}
478  if(i==2 && j== 9){ printf(" JPB=%4d thrD=%4d",mTrgSim->jet[geo.ns][1],mTrgSim->JPBCTHRD);}
479  if(i==2 && j==10){ printf(" JPC=%4d thrD=%4d",mTrgSim->jet[geo.ns][2],mTrgSim->JPBCTHRD);}
480  if(i==2 && j==11){ printf(" JPD=%4d thrD=%4d",mTrgSim->jet[geo.ns][3],mTrgSim->JPDETHRD);}
481  if(i==2 && j==12){ printf(" JPE=%4d thrD=%4d",mTrgSim->jet[geo.ns][4],mTrgSim->JPDETHRD);}
482  printf("\n");
483  }
484  }
485  }
486 }
487 
488 void StFcsTriggerSimMaker::print4B4(){
489  //printout ecal 4x4
490  FILE* f1=fopen("EH4by4.txt","w");
491  FILE* f2=fopen("EH4by4dist.txt","w");
492  FILE* f3=fopen("EH4by4map.txt","w");
493  FILE* f4=fopen("EH2by2dist.txt","w");
494  FILE* f5=fopen("EH2by2map.txt","w");
495  FILE* f6=fopen("EH2by2map2.txt","w");
496 
497  // v1 with hcal top2/bottom2 rows not in trigger
498  //enum {EX2B2=10,EY2B2=16,HX2B2=6,HY2B2=8};
499  //enum {EXOFF=1,EYOFF=1,HXOFF=0,HYOFF=2};
500  //v2 with hcal top2/bottom2 rows in trigger
501  enum {EX2B2=10,EY2B2=16,HX2B2=6,HY2B2=10};
502  enum {EXOFF=1,EYOFF=1,HXOFF=0,HYOFF=0};
503  // v3 with hcal top2/bottom2 rows in trigger, put offset of 1 in HX
504  //enum {EX2B2=10,EY2B2=16,HX2B2=6,HY2B2=10};
505  //enum {EXOFF=1,EYOFF=1,HXOFF=1,HYOFF=0};
506  StThreeVectorF exyz[2][EX2B2-1][EY2B2-1];
507  StThreeVectorF hxyz[2][HX2B2-1][HY2B2-1]; //4x4
508  StThreeVectorF Hxyz[2][HX2B2][HY2B2]; //2x2
509  StThreeVectorF sxyz[2][EX2B2-1][EY2B2-1]; //ecal 4x4 extraprated at hcal
510  StThreeVectorF eoff = mFcsDb->getDetectorOffset(1);
511  StThreeVectorF hoff = mFcsDb->getDetectorOffset(3);
512  float esmx=mFcsDb->getShowerMaxZ(1);
513  float hsmx=mFcsDb->getShowerMaxZ(3);
514  float r1 = sqrt(eoff.x()*eoff.x()+eoff.z()*eoff.z()); //distance from IP to ecal surface
515  float r2 = r1 + esmx; //distance from IP to ecal SMax
516  float r3 = sqrt(hoff.x()*hoff.x()+hoff.z()*hoff.z()); //distance from IP to hcal surface
517  float r4 = r3 + hsmx; //distance from IP to hcal SMax
518  float sf = r4/r2; // scale factor for extraporating ecal point to hcal plane
519  fprintf(f1,"Distannce from IP to EcalSmax=%8.3f HcalSMax=%8.3f Ratio=%6.3f\n",r2,r4,sf);
520 
521  fprintf(f1,"\nHcal 4x4\n");
522  fprintf(f1," C R XY[cell] XYZ[cm]\n");
523  for(int ns=1; ns<2; ns++){
524  for(int j=0; j<HY2B2-1; j++){
525  float y=j*2 + HYOFF + 2;
526  for(int i=0; i<HX2B2-1; i++){
527  float x=i*2 + HXOFF + 2;
528  hxyz[ns][i][j] = mFcsDb->getStarXYZfromColumnRow(ns+2,x,y);
529  fprintf(f1,"H %2d %2d %4.1f %4.1f %8.2f %8.2f %8.2f\n",
530  i,j,x,y,hxyz[ns][i][j].x(),hxyz[ns][i][j].y(),hxyz[ns][i][j].z());
531  }
532  }
533  }
534 
535  fprintf(f1,"\nHcal 2x2\n");
536  fprintf(f1," C R XY[cell] XYZ[cm]\n");
537  for(int ns=1; ns<2; ns++){
538  for(int j=0; j<HY2B2; j++){
539  float y=j*2 + HYOFF + 1;
540  for(int i=0; i<HX2B2; i++){
541  float x=i*2 + HXOFF + 1;
542  Hxyz[ns][i][j] = mFcsDb->getStarXYZfromColumnRow(ns+2,x,y);
543  fprintf(f1,"H %2d %2d %4.1f %4.1f %8.2f %8.2f %8.2f\n",
544  i,j,x,y,Hxyz[ns][i][j].x(),Hxyz[ns][i][j].y(),Hxyz[ns][i][j].z());
545  }
546  }
547  }
548 
549  fprintf(f1,"\nEcal 4x4\n");
550  fprintf(f1," C R XY[cell] XYZ[cm] | XYZ at HCAL | C R distance XYZ of closest Hcal 4x4\n");
551 
552  fprintf(f3,"static const int EtoHmap[%d][%d][2] = {\n",EY2B2-1,EX2B2-1);
553  fprintf(f5,"static const int EtoH2map[%d][%d][2] = {\n",EY2B2-1,EX2B2-1);
554  fprintf(f6,"static const int EtoH3map[%d][%d][4] = {\n",EY2B2-1,EX2B2-1);
555  for(int ns=1; ns<2; ns++){
556  for(int j=0; j<EY2B2-1; j++){
557  float y=j*2 + EYOFF + 2;
558  fprintf(f3," { ");
559  for(int i=0; i<EX2B2-1; i++){
560  float x=i*2 + EXOFF + 2;
561  exyz[ns][i][j] = mFcsDb->getStarXYZfromColumnRow(ns,x,y);
562  sxyz[ns][i][j] = sf * exyz[ns][i][j];
563  //search for closest hcal 4x4
564  float dmin=999.0;
565  int k,l;
566  for(int ii=0; ii<HX2B2-1; ii++){
567  for(int jj=0; jj<HY2B2-1; jj++){
568  StThreeVectorF diff = sxyz[ns][i][j]-hxyz[ns][ii][jj];
569  float d = diff.mag();
570  //printf("%2d %2d %2d %2d d=%6.2f diff=%8.2f %8.2f %8.2f\n",
571  // i,j,ii,jj,d,diff.x(),diff.y(),diff.z());
572  if(d<dmin){
573  dmin=d;
574  k=ii;
575  l=jj;
576  }
577  }
578  }
579  fprintf(f1,"E %2d %2d %4.1f %4.1f %8.2f %8.2f %8.2f | %8.2f %8.2f %8.2f | %2d %2d %6.2f %8.2f %8.2f %8.2f\n",
580  i,j,x,y,
581  exyz[ns][i][j].x(),exyz[ns][i][j].y(),exyz[ns][i][j].z(),
582  sxyz[ns][i][j].x(),sxyz[ns][i][j].y(),sxyz[ns][i][j].z(),
583  k,l,dmin,
584  hxyz[ns][k][l].x(),hxyz[ns][k][l].y(),hxyz[ns][k][l].z());
585  fprintf(f2,"%6.2f ",dmin);
586  if(i==EX2B2-2) fprintf(f2,"\n");
587  fprintf(f3,"{%2d,%2d}",l,k);
588  if(i<EX2B2-2) fprintf(f3,",");
589  if(i==EX2B2-2) fprintf(f3,"}");
590 
591  //hcal 2x2
592  float dmin2=999.0;
593  int kk,ll;
594  for(int ii=0; ii<HX2B2; ii++){
595  for(int jj=0; jj<HY2B2; jj++){
596  StThreeVectorF diff = sxyz[ns][i][j]-Hxyz[ns][ii][jj];
597  float d = diff.mag();
598  //printf("%2d %2d %2d %2d d=%6.2f diff=%8.2f %8.2f %8.2f\n",
599  // i,j,ii,jj,d,diff.x(),diff.y(),diff.z());
600  if(d<dmin2){
601  dmin2=d;
602  kk=ii;
603  ll=jj;
604  }
605  }
606  }
607  fprintf(f1,"E %2d %2d %4.1f %4.1f %8.2f %8.2f %8.2f | %8.2f %8.2f %8.2f | %2d %2d %6.2f %8.2f %8.2f %8.2f\n",
608  i,j,x,y,
609  exyz[ns][i][j].x(),exyz[ns][i][j].y(),exyz[ns][i][j].z(),
610  sxyz[ns][i][j].x(),sxyz[ns][i][j].y(),sxyz[ns][i][j].z(),
611  kk,ll,dmin2,
612  Hxyz[ns][kk][ll].x(),Hxyz[ns][kk][ll].y(),Hxyz[ns][kk][ll].z());
613  fprintf(f4,"%6.2f ",dmin2);
614  if(i==EX2B2-2) fprintf(f4,"\n");
615 
616  fprintf(f5,"{%2d,%2d}",ll,kk);
617  if(i<EX2B2-2) fprintf(f5,",");
618  if(i==EX2B2-2) fprintf(f5,"}");
619 
620  int i1=kk-1, i2=kk, j1=ll-1, j2=ll;
621  int r1=-1, r2=-1, r3=-1, r4=-1;
622  if(i1>=0 && i1<HX2B2-1 && j1>=0 && j1<HY2B2-1) r1=i1+j1*(HX2B2-1);
623  if(i2>=0 && i2<HX2B2-1 && j1>=0 && j1<HY2B2-1) r2=i2+j1*(HX2B2-1);
624  if(i1>=0 && i1<HX2B2-1 && j2>=0 && j2<HY2B2-1) r3=i1+j2*(HX2B2-1);
625  if(i2>=0 && i2<HX2B2-1 && j2>=0 && j1<HY2B2-1) r4=i2+j2*(HX2B2-1);
626  fprintf(f6,"{%2d,%2d,%2d,%2d}",r1,r2,r3,r4);
627  if(i<EX2B2-2) fprintf(f6,",");
628  if(i==EX2B2-2) fprintf(f6,"}");
629  }
630  if(j<EY2B2-2) fprintf(f3,",\n");
631  if(j==EY2B2-2) fprintf(f3,"\n");
632  if(j<EY2B2-2) fprintf(f5,",\n");
633  if(j==EY2B2-2) fprintf(f5,"\n");
634  if(j<EY2B2-2) fprintf(f6,",\n");
635  if(j==EY2B2-2) fprintf(f6,"\n");
636  }
637  fprintf(f3,"}\n");
638  fprintf(f5,"}\n");
639  fprintf(f6,"}\n");
640  }
641  fclose(f1);
642  fclose(f2);
643  fclose(f3);
644  fclose(f4);
645  fclose(f5);
646  fclose(f6);
647  return;
648 }
649 
650 
651 void StFcsTriggerSimMaker::printJP(){
652  //printout ecal 4x4
653  FILE* f1=fopen("EHJP.txt","w");
654  FILE* f2=fopen("EHJPdist.txt","w");
655 
656  enum {EXOFF=1,EYOFF=1,HXOFF=0,HYOFF=2};
657  StThreeVectorF exyz[2][3][5];
658  StThreeVectorF hxyz[2][3][5];
659  StThreeVectorF sxyz[2][3][5]; //ecal 4x4 extraprated at hcal
660  StThreeVectorF eoff = mFcsDb->getDetectorOffset(1);
661  StThreeVectorF hoff = mFcsDb->getDetectorOffset(3);
662  float esmx=mFcsDb->getShowerMaxZ(1);
663  float hsmx=mFcsDb->getShowerMaxZ(3);
664  float r1 = sqrt(eoff.x()*eoff.x()+eoff.z()*eoff.z()); //distance from IP to ecal surface
665  float r2 = r1 + esmx; //distance from IP to ecal SMax
666  float r3 = sqrt(hoff.x()*hoff.x()+hoff.z()*hoff.z()); //distance from IP to hcal surface
667  float r4 = r3 + hsmx; //distance from IP to hcal SMax
668  float sf = r4/r2; // scale factor for extraporating ecal point to hcal plane
669  fprintf(f1,"Distannce from IP to EcalSmax=%8.3f HcalSMax=%8.3f Ratio=%6.3f\n",r2,r4,sf);
670  fprintf(f1,"\nHcal 8x8\n");
671  fprintf(f1," C R XY[cell] XYZ[cm]\n");
672  for(int ns=1; ns<2; ns++){
673  for(int j=0; j<5; j++){
674  float y=j*2 + HYOFF + 4;
675  for(int i=0; i<3; i++){
676  float x=i*2 + HXOFF + 4;
677  hxyz[ns][i][j] = mFcsDb->getStarXYZfromColumnRow(ns+2,x,y);
678  fprintf(f1,"H %2d %2d %4.1f %4.1f %8.2f %8.2f %8.2f\n",
679  i,j,x,y,hxyz[ns][i][j].x(),hxyz[ns][i][j].y(),hxyz[ns][i][j].z());
680  }
681  }
682  }
683  fprintf(f1,"\nEcal 12x16\n");
684  fprintf(f1," C R XY[cell] XYZ[cm] | XYZ at HCAL | distance XYZ of Hcal 8x8\n");
685  for(int ns=1; ns<2; ns++){
686  for(int j=0; j<5; j++){
687  float y=j*4 + EYOFF + 8;
688  for(int i=0; i<3; i++){
689  float x=i*4 + EXOFF + 6;
690  exyz[ns][i][j] = mFcsDb->getStarXYZfromColumnRow(ns,x,y);
691  sxyz[ns][i][j] = sf * exyz[ns][i][j];
692  StThreeVectorF diff = sxyz[ns][i][j]-hxyz[ns][i][j];
693  float d = diff.mag();
694  fprintf(f1,"E %2d %2d %4.1f %4.1f %8.2f %8.2f %8.2f | %8.2f %8.2f %8.2f | %6.2f %8.2f %8.2f %8.2f\n",
695  i,j,x,y,
696  exyz[ns][i][j].x(),exyz[ns][i][j].y(),exyz[ns][i][j].z(),
697  sxyz[ns][i][j].x(),sxyz[ns][i][j].y(),sxyz[ns][i][j].z(),
698  d,
699  hxyz[ns][i][j].x(),hxyz[ns][i][j].y(),hxyz[ns][i][j].z());
700  fprintf(f2,"%6.2f ",d);
701  if(i==2) fprintf(f2,"\n");
702  }
703  }
704  }
705  fclose(f1);
706  fclose(f2);
707  return;
708 }
709 
710 void StFcsTriggerSimMaker::readThresholdFile(){
711  LOG_INFO << Form("Reading Thresholds from %s",mThresholdFile)<<endm;
712  FILE* F=fopen(mThresholdFile,"r");
713  if(F == NULL){
714  LOG_ERROR << Form("Could not open %s",mThresholdFile)<<endm;
715  return;
716  }
717  char f[10],name[100];
718  int i,id,v;
719  while(fscanf(F,"%s %d %d %s %d",f, &i, &id, name, &v) != EOF){
720  if(f[0] == '#') continue;
721  TString trg(name);
722  printf("Reading Threshold %s %d\n",trg.Data(),v);
723  if (trg=="FCS_PHTTHR") mTrgSim->PHTTHR=v;
724  else if(trg=="FCS_EM-HERATIO-THR") mTrgSim->EM_HERATIO_THR=v;
725  else if(trg=="FCS_HAD-HERATIO-THR") mTrgSim->HAD_HERATIO_THR=v;
726  else if(trg=="FCS_EMTHR2") mTrgSim->EMTHR2=v;
727  else if(trg=="FCS_EMTHR1") mTrgSim->EMTHR1=v;
728  else if(trg=="FCS_EMTHR0") mTrgSim->EMTHR0=v;
729  else if(trg=="FCS_ELETHR2") mTrgSim->ELETHR2=v;
730  else if(trg=="FCS_ELETHR1") mTrgSim->ELETHR1=v;
731  else if(trg=="FCS_ELETHR0") mTrgSim->ELETHR0=v;
732  else if(trg=="FCS_HADTHR2") mTrgSim->HADTHR2=v;
733  else if(trg=="FCS_HADTHR1") mTrgSim->HADTHR1=v;
734  else if(trg=="FCS_HADTHR0") mTrgSim->HADTHR0=v;
735  else if(trg=="FCS_JPATHR2") mTrgSim->JPATHR2=v;
736  else if(trg=="FCS_JPATHR1") mTrgSim->JPATHR2=v;
737  else if(trg=="FCS_JPATHR0") mTrgSim->JPATHR0=v;
738  else if(trg=="FCS_JPBCTHR2") mTrgSim->JPBCTHR2=v;
739  else if(trg=="FCS_JPBCTHR1") mTrgSim->JPBCTHR1=v;
740  else if(trg=="FCS_JPBCTHR0") mTrgSim->JPBCTHR0=v;
741  else if(trg=="FCS_JPBCTHRD") mTrgSim->JPBCTHRD=v;
742  else if(trg=="FCS_JPDETHR2") mTrgSim->JPDETHR2=v;
743  else if(trg=="FCS_JPDETHR1") mTrgSim->JPDETHR1=v;
744  else if(trg=="FCS_JPDETHR0") mTrgSim->JPDETHR0=v;
745  else if(trg=="FCS_JPDETHRD") mTrgSim->JPDETHRD=v;
746  else if(trg=="FCS_EHTTHR") mTrgSim->EHTTHR=v;
747  else if(trg=="FCS_HHTTHR") mTrgSim->HHTTHR=v;
748  else if(trg=="FCS_ETOTTHR") mTrgSim->ETOTTHR=v;
749  else if(trg=="FCS_HTOTTHR") mTrgSim->HTOTTHR=v;
750  else{
751  printf("No Threshold called %s %d\n",trg.Data(),v);
752  }
753  }
754  fclose(F);
755 }
756 
757 //read from offline copy of online runlog DB
758 void StFcsTriggerSimMaker::readThresholdDb(){
759  //to be implemented before run22 online DB moves to offline
760 }
761 
762 template<typename T> void StFcsTriggerSimMaker::feedADC(T* hit, int ns, int ehp, uint16_t data_array[]){
763 
764  unsigned short dep = hit->dep();
765  unsigned short ch = hit->channel();
766 
767 // printf("ns=%1d ehp=%1d dep=%2d ch=%2d adc=%4d sim=%d\n",ns,ehp,dep,ch,hit->adc(0),mSimMode);
768  fcs_trg_sim_adc[ns][ehp][dep][ch] = hit->adc(0);
769  if(mSimMode==0){
770  int ntb=hit->nTimeBin();
771  for(int t=0; t<ntb; t++){
772  int tb = hit->timebin(t);
773  if(tb>=mTrgTimebin-3 && tb<=mTrgTimebin+4){
774  data_array[tb-mTrgTimebin+3] = hit->adc(t);
775 // printf("tb=%3d i=%2d adc=%4d\n",tb,tb-mTrgTimebin+3,hit->adc(t));
776  }
777  }
778  mTrgSim->fill_event(ehp,ns,dep,ch,data_array,8) ;
779  }else{
780  data_array[1] = hit->adc(0)-1; //removing 1 to add at tb6
781  data_array[6] = 1; //add this so tb6>tb7
782  mTrgSim->fill_event(ehp,ns,dep,ch,data_array,8) ;
783  }
784  if(mFile) fprintf(mFile,"%2d %2d %2d %2d %5d\n",ns,ehp,dep,ch,hit->adc(0));
785 
786 }
787 
788 void StFcsTriggerSimMaker::readTrgId(){
789  int i;
790  char trgn[200];
791  if(mTrgIdFile){
792  LOG_INFO<<"Reading "<<mTrgIdFile<<endm;
793  FILE* fp=fopen(mTrgIdFile,"r");
794  if(!fp) {LOG_WARN << "Cannot open "<<mTrgIdFile<<endm; return;}
795  while(!feof(fp)) {
796  fscanf(fp,"%d %s",&i,trgn);
797  mTrgIdName[i]=trgn;
798  LOG_INFO << "TRGID "<<i<<"="<<mTrgIdName[i].Data()<<endm;
799  }
800  fclose(fp);
801  }
802 }
StMuDst::muFcsCollection
static StMuFcsCollection * muFcsCollection()
returns pointer to current StMuFcsCollection
Definition: StMuDst.h:395
StFcsDb::getGainOnline
float getGainOnline(int det, int id) const
get the pres valley position for cut
Definition: StFcsDb.cxx:984
StFcsDb::northSouth
int northSouth(int det) const
Ecal=0, Hcal=1, Pres=2.
Definition: StFcsDb.cxx:430
StFcsDb::getDetectorOffset
StThreeVectorD getDetectorOffset(int det, double zdepth=-1) const
Utility functions related to DetectorPosition.
Definition: StFcsDb.cxx:578
StFcsDb::getStarXYZfromColumnRow
StThreeVectorD getStarXYZfromColumnRow(int det, float col, float row, float FcsZ=-1.0) const
get the STAR frame cooridnates from other way
Definition: StFcsDb.cxx:863
link_t
Definition: fcs_trg_base.h:10
StFcsDb::ecalHcalPres
int ecalHcalPres(int det) const
Ecal North=0, Ecal South=1, Hcal North=2, Hcal South=3, Pres=4/5.
Definition: StFcsDb.cxx:423
StMuDst::event
static StMuEvent * event()
returns pointer to current StMuEvent (class holding the event wise information, e.g. event number, run number)
Definition: StMuDst.h:320
StarGenEvent
Base class for event records.
Definition: StarGenEvent.h:81
StarPrimaryMaker::event
StarGenEvent * event()
Return a pointer to the event.
Definition: StarPrimaryMaker.h:169
kStOK
Definition: Stypes.h:40
StarPrimaryMaker
Main steering class for event generation.
Definition: StarPrimaryMaker.h:103
StFcsDb::maxId
int maxId(int det) const
number of column
Definition: StFcsDb.cxx:468
StarGenEvent::GetFilterResult
UInt_t GetFilterResult()
Returns the filter result.
Definition: StarGenEvent.h:206
geom_t
Definition: fcs_trg_base.h:19
kStErr
Definition: Stypes.h:44
fcs_trg_base
Definition: fcs_trg_base.h:33
StFcsDb::getDepfromId
void getDepfromId(int detectorId, int id, int &ehp, int &ns, int &crt, int &slt, int &dep, int &ch) const
print ET gain
Definition: StFcsDb.cxx:1018
Generated by   doxygen 1.8.5