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  //EPD mask
156  if(mPresMask){
157  printf("Reading PresMask from %s\n",mPresMask);
158  FILE* F=fopen(mPresMask,"r");
159  if(F==NULL){
160  printf("Cannot open %s\n",mPresMask);
161  }else{
162  char line[512];
163  int r,c,m[6];
164  while(fgets(line,sizeof(line),F)){
165  if(line[0]=='#') {
166  printf("%s",line);
167  continue;
168  }
169  sscanf(line,"%d %d %x %x %x %x %x %x",&r,&c,&m[0],&m[1],&m[2],&m[3],&m[4],&m[5]);
170  printf("%2d %1d %08x %08x %08x %08x %08x %08x\n",r,c,m[0],m[1],m[2],m[3],m[4],m[5]);
171  for(int i=0; i<6; i++) mTrgSim->PRES_MASK[r-1][c-1][i]=m[i];
172  }
173  mTrgSim->fcs_readPresMaskFromText=1;
174  }
175  }
176 
177  memset(NTRG,0,sizeof(NTRG));
178  return kStOK;
179 }
180 
181 int StFcsTriggerSimMaker::InitRun(int runNumber){
182  LOG_INFO << "StFcsTriggerSimMaker::InitRun" << endm;
183  //print out 4x4 and JP info
184  if(mDebug>0){
185  print4B4();
186  printJP();
187  }
188 
189  //QA root file
190  if(mQaTreeFilename){
191  mQaTreeFile=new TFile(mQaTreeFilename,"RECREATE");
192  mTree = new TTree("trgsim","trigger sim QA");
193  mTree->Branch("flt",&mFlt,"flt/I");
194  mTree->Branch("trg",&mTrg,"trg/I");
195  }
196  if(mQaHistFilename){
197  mQaHistFile=new TFile(mQaHistFilename,"RECREATE");
198  mTrgRate = new TH1F("FcsTrgRate","FcsTrgRate",mNTRG+1,0,mNTRG+1);
199  }
200 
201  //Write Text event file & gainfile
202  FILE* gainfile=0;
203  FILE* gainfile2=0;
204  if(mFilename){
205  mFile = fopen(mFilename,"w");
206  gainfile=fopen("fcs_et_gain.txt","w");
207  gainfile2=fopen("fcs_et_gain2.txt","w");
208  }
209 
210  //Fill ETgain, GainCorr and Pedestal
211  for(int det=0; det<=kFcsNDet; det++){
212  int nid=mFcsDb->maxId(det);
213  for(int id=0; id<nid; id++){
214  int ehp,ns,crt,sub,dep,ch;
215  mFcsDb->getDepfromId(det,id,ehp,ns,crt,sub,dep,ch);
216  if(det<4){
217  fcs_trg_pt_correction[ns][ehp][dep][ch] = mFcsDb->getEtGain(det,id,mEtFactor);
218  fcs_trg_gain_correction[ns][ehp][dep][ch] = mFcsDb->getGainCorrection(det,id);
219  }else{
220  fcs_trg_pt_correction[ns][ehp][dep][ch] = 1.0;
221  fcs_trg_gain_correction[ns][ehp][dep][ch] = 1.0;
222  }
223  fcs_trg_pedestal[ns][ehp][dep][ch] = 0;
224 
225  mTrgSim->p_g[ns][ehp][dep][ch].ped = fcs_trg_pedestal[ns][ehp][dep][ch];
226 
227  float ggg = fcs_trg_pt_correction[ns][ehp][dep][ch];
228  //float ggg = (fcs_trg_pt_correction[ns][ehp][dep][ch]-1.0)/2.0 + 1.0;
229  float gg = ggg * fcs_trg_gain_correction[ns][ehp][dep][ch];
230  int g = (uint32_t)(gg*256.0+0.5) ;
231  mTrgSim->p_g[ns][ehp][dep][ch].gain = g;
232 
233  /*
234  printf("AAAGAIN %1d %1d %2d %2d pT=%6.3f corr=%6.3f ped=%4d\n",ns,ehp,dep,ch,
235  fcs_trg_pt_correction[ns][ehp][dep][ch],
236  fcs_trg_gain_correction[ns][ehp][dep][ch],
237  fcs_trg_pedestal[ns][ehp][dep][ch]);
238  */
239 
240  if(gainfile)
241  fprintf(gainfile,"%2d %2d %2d %2d %8.3f\n",ns,ehp,dep,ch,
242  fcs_trg_pt_correction[ns][ehp][dep][ch]);
243  if(gainfile2)
244  fprintf(gainfile2,"%2d %2d %2d %2d %8.3f\n",ns,ehp,dep,ch,
245  (fcs_trg_pt_correction[ns][ehp][dep][ch]-1.0)/2.0 + 1.0);
246  }
247  }
248  if(gainfile) fclose(gainfile);
249  if(gainfile2) fclose(gainfile2);
250  return kStOK;
251 }
252 
253 int StFcsTriggerSimMaker::Finish(){
254  mTrgSim->run_stop();
255  if(mFile) {
256  printf("Closing %s\n",mFilename);
257  fclose(mFile);
258  }
259  if(mQaTreeFile){
260  printf("Closing %s\n",mQaTreeFilename);
261  mTree->Write();
262  mQaTreeFile->Close();
263  }
264  if(mQaHistFile){
265  printf("Closing %s\n",mQaHistFilename);
266  mTrgRate->Write();
267  mQaHistFile->Close();
268  }
269 
270  int tot = NTRG[mNTRG];
271  LOG_INFO << "Triggers counts/"<<tot<<" (rate at 5MHz BBC)"<<endm;
272  for(int i=0; i<mNTRG; i++)
273  LOG_INFO << Form("%8s %9d (%12.2f)",ctrg[i],NTRG[i],double(NTRG[i])/tot*5.0e6)<<endm;
274 
275  return kStOK;
276 }
277 
278 int StFcsTriggerSimMaker::Make(){
279  StEvent* event = nullptr;
280  event = (StEvent*)GetInputDS("StEvent");
281  if(!event) {LOG_INFO << "StFcsTriggerSimMaker::Make did not find StEvent"<<endm;}
282  mFcsColl = event->fcsCollection();
283  if(!mFcsColl) {LOG_INFO << "StFcsTriggerSimMaker::Make did not find StEvent->StFcsCollection"<<endm;}
284 
285  StMuEvent* muevent = nullptr;
286  if((!event)||(!mFcsColl)){
287 
288  LOG_INFO<<"No StEvent info available for StFcsTriggerSimMaker. "<< endm;
289 
290  muevent = StMuDst::event();
291  mMuFcsColl = StMuDst::muFcsCollection();
292  if(muevent && mMuFcsColl){
293  LOG_INFO <<"Proceeding with StMuDst info to be used with StfcsTriggerSimMaker."<< endm;
294  }else{
295  LOG_ERROR << "StFcsTriggerSimMaker::Make did not find StEvent and MuEvent." << endm;
296  return kStErr;
297  }
298  }
299 
300  mTrgSim->start_event();
301 
302  //Feed ADC
303  static uint16_t data[8];
304  memset(data,0,sizeof(data)) ;
305  memset(fcs_trg_sim_adc,0,sizeof(fcs_trg_sim_adc));
306  int n=0;
307  for(int det=0; det<=kFcsNDet; det++){
308 
309  int ns = mFcsDb->northSouth(det);
310  int ehp = mFcsDb->ecalHcalPres(det);
311 
312  if((event)&&(mFcsColl)){
313  StSPtrVecFcsHit& hits = mFcsColl->hits(det);
314 
315  int nh = mFcsColl->numberOfHits(det);
316 
317  for(int i=0; i<nh; i++){
318  StFcsHit* hit=hits[i];
319  unsigned short ch = hit->channel();
320 
321  if(ehp<0 || ch>=32) continue;
322  feedADC(hit, ns, ehp, data);
323  n++;
324  }
325 
326  }else if((muevent)&&(mMuFcsColl)){ //Use StMuDst info instead of StEvent for Trigger Reconstruction
327 
328  int nh = mMuFcsColl->numberOfHits(det);
329  int det_hit_index = mMuFcsColl->indexOfFirstHit(det);
330 
331  for(int i=0; i<nh; i++){
332 
333  int hit_index = i + det_hit_index;
334  StMuFcsHit* hit = mMuFcsColl->getHit(hit_index);
335  unsigned short ch = hit->channel();
336 
337  if(ehp<0 || ch>=32) continue;
338  feedADC(hit, ns, ehp, data);
339  n++;
340  }
341  }
342  }
343  if(mFile) fprintf(mFile,"%2d %2d %2d %2d %5d\n",-1,0,0,0,0);
344  LOG_INFO << Form("StFcsTriggerSimMaker feeded %d hits",n) << endm;;
345 
346  //Run Trigger Simulation
347  // uint16_t dsm_out = fcs_trg_run(mTrgSelect, mDebug);
348  uint32_t dsm_out = mTrgSim->end_event();
349 
350  //QA Tree
351  mFlt=0;
352  StarPrimaryMaker* pmkr= static_cast<StarPrimaryMaker*>(GetMaker("PrimaryMaker"));
353  if(pmkr){
354  StarGenEvent *ge = pmkr->event();
355  if(ge){
356  mFlt=ge->GetFilterResult();
357  }
358  }
359  mTcu=dsm_out;
360  if(mQaTreeFile) mTree->Fill();
361 
362  //Results
363  LOG_INFO << Form("Output to TCU = 0x%08x Filter=0x%08x",mTcu,mFlt)<<endm;
364 
365  //TCU Trigger Emulation
366  int trg[mNTRG];
367  mTrg=0;
368  memset(trg,0,sizeof(trg));
369  if((dsm_out >> 6) & 0x1) {trg[ 0]=1; mTrg+=1<<0;} //JP2
370  if((dsm_out >> 7) & 0x1) {trg[ 1]=1; mTrg+=1<<1;} //JPA1
371  if((dsm_out >>10) & 0x1) {trg[ 2]=1; mTrg+=1<<2;} //JPA0
372  if((dsm_out >> 8) & 0x1) {trg[ 3]=1; mTrg+=1<<3;} //JPBC1
373  if((dsm_out >>11) & 0x1) {trg[ 4]=1; mTrg+=1<<4;} //JPBC0
374  if((dsm_out >> 9) & 0x1) {trg[ 5]=1; mTrg+=1<<5;} //JPDE1
375  if((dsm_out >>12) & 0x1) {trg[ 6]=1; mTrg+=1<<6;} //JPDE0
376  if((dsm_out >>13) & 0x1) {trg[ 7]=1; mTrg+=1<<7;} //DiJP
377  if((dsm_out >>14) & 0x1) {trg[ 8]=1; mTrg+=1<<8;} //DiJpAsy
378  if( ((dsm_out >>18) & 0x1) &&
379  ((dsm_out >>26) & 0x1) ) {trg[ 9]=1; mTrg+=1<<9;} //DY
380  if( ((dsm_out >>17) & 0x1) &&
381  ((dsm_out >>25) & 0x1) ) {trg[10]=1; mTrg+=1<<10;} //Jpsi
382  if( ((dsm_out >>19) & 0x1) &&
383  ((dsm_out >>27) & 0x1) ) {trg[11]=1; mTrg+=1<<11;} //DYNoEpd
384  if((dsm_out >>15) & 0x1) {trg[12]=1; mTrg+=1<<12;} //DYAsy
385  if((dsm_out >> 2) & 0x1) {trg[13]=1; mTrg+=1<<13;} //Had2
386  if((dsm_out >> 1) & 0x1) {trg[14]=1; mTrg+=1<<14;} //Had1
387  if((dsm_out >> 0) & 0x1) {trg[15]=1; mTrg+=1<<15;} //Had0
388  if((dsm_out >> 5) & 0x1) {trg[16]=1; mTrg+=1<<16;} //EM2
389  if((dsm_out >> 4) & 0x1) {trg[17]=1; mTrg+=1<<17;} //EM1
390  if((dsm_out >> 3) & 0x1) {trg[18]=1; mTrg+=1<<18;} //EM0
391  if( ((dsm_out >>18) & 0x1) ||
392  ((dsm_out >>26) & 0x1) ) {trg[19]=1; mTrg+=1<<19;} //ELE2
393  if( ((dsm_out >>19) & 0x1) ||
394  ((dsm_out >>27) & 0x1) ) {trg[20]=1; mTrg+=1<<20;} //EM3
395 
396  if(mTrgRate) mTrgRate->Fill(mNTRG);
397  NTRG[mNTRG]++;
398  for(int i=0; i<mNTRG; i++){
399  if(trg[i]) {
400  if(mTrgRate) mTrgRate->Fill(i);
401  NTRG[i]++;
402  }
403  }
404 
405  LOG_INFO << "Triggers = ";
406  for(int i=0; i<mNTRG; i++){ if(trg[i]) LOG_INFO << ctrg[i] << " ";}
407  LOG_INFO << endm;
408 
409  return kStOK;
410 }
411 
412 void StFcsTriggerSimMaker::runStage2(link_t ecal[], link_t hcal[], link_t pres[], geom_t& geo, link_t output[]){
413  uint16_t s2_to_dsm;
414  mTrgSim->stage_2(ecal,hcal,pres,geo,output,&s2_to_dsm);
415 }
416 
417 void StFcsTriggerSimMaker::print4B4(){
418  //printout ecal 4x4
419  FILE* f1=fopen("EH4by4.txt","w");
420  FILE* f2=fopen("EH4by4dist.txt","w");
421  FILE* f3=fopen("EH4by4map.txt","w");
422  FILE* f4=fopen("EH2by2dist.txt","w");
423  FILE* f5=fopen("EH2by2map.txt","w");
424  FILE* f6=fopen("EH2by2map2.txt","w");
425 
426  // v1 with hcal top2/bottom2 rows not in trigger
427  //enum {EX2B2=10,EY2B2=16,HX2B2=6,HY2B2=8};
428  //enum {EXOFF=1,EYOFF=1,HXOFF=0,HYOFF=2};
429  //v2 with hcal top2/bottom2 rows in trigger
430  enum {EX2B2=10,EY2B2=16,HX2B2=6,HY2B2=10};
431  enum {EXOFF=1,EYOFF=1,HXOFF=0,HYOFF=0};
432  // v3 with hcal top2/bottom2 rows in trigger, put offset of 1 in HX
433  //enum {EX2B2=10,EY2B2=16,HX2B2=6,HY2B2=10};
434  //enum {EXOFF=1,EYOFF=1,HXOFF=1,HYOFF=0};
435  StThreeVectorF exyz[2][EX2B2-1][EY2B2-1];
436  StThreeVectorF hxyz[2][HX2B2-1][HY2B2-1]; //4x4
437  StThreeVectorF Hxyz[2][HX2B2][HY2B2]; //2x2
438  StThreeVectorF sxyz[2][EX2B2-1][EY2B2-1]; //ecal 4x4 extraprated at hcal
439  StThreeVectorF eoff = mFcsDb->getDetectorOffset(1);
440  StThreeVectorF hoff = mFcsDb->getDetectorOffset(3);
441  float esmx=mFcsDb->getShowerMaxZ(1);
442  float hsmx=mFcsDb->getShowerMaxZ(3);
443  float r1 = sqrt(eoff.x()*eoff.x()+eoff.z()*eoff.z()); //distance from IP to ecal surface
444  float r2 = r1 + esmx; //distance from IP to ecal SMax
445  float r3 = sqrt(hoff.x()*hoff.x()+hoff.z()*hoff.z()); //distance from IP to hcal surface
446  float r4 = r3 + hsmx; //distance from IP to hcal SMax
447  float sf = r4/r2; // scale factor for extraporating ecal point to hcal plane
448  fprintf(f1,"Distannce from IP to EcalSmax=%8.3f HcalSMax=%8.3f Ratio=%6.3f\n",r2,r4,sf);
449 
450  fprintf(f1,"\nHcal 4x4\n");
451  fprintf(f1," C R XY[cell] XYZ[cm]\n");
452  for(int ns=1; ns<2; ns++){
453  for(int j=0; j<HY2B2-1; j++){
454  float y=j*2 + HYOFF + 2;
455  for(int i=0; i<HX2B2-1; i++){
456  float x=i*2 + HXOFF + 2;
457  hxyz[ns][i][j] = mFcsDb->getStarXYZfromColumnRow(ns+2,x,y);
458  fprintf(f1,"H %2d %2d %4.1f %4.1f %8.2f %8.2f %8.2f\n",
459  i,j,x,y,hxyz[ns][i][j].x(),hxyz[ns][i][j].y(),hxyz[ns][i][j].z());
460  }
461  }
462  }
463 
464  fprintf(f1,"\nHcal 2x2\n");
465  fprintf(f1," C R XY[cell] XYZ[cm]\n");
466  for(int ns=1; ns<2; ns++){
467  for(int j=0; j<HY2B2; j++){
468  float y=j*2 + HYOFF + 1;
469  for(int i=0; i<HX2B2; i++){
470  float x=i*2 + HXOFF + 1;
471  Hxyz[ns][i][j] = mFcsDb->getStarXYZfromColumnRow(ns+2,x,y);
472  fprintf(f1,"H %2d %2d %4.1f %4.1f %8.2f %8.2f %8.2f\n",
473  i,j,x,y,Hxyz[ns][i][j].x(),Hxyz[ns][i][j].y(),Hxyz[ns][i][j].z());
474  }
475  }
476  }
477 
478  fprintf(f1,"\nEcal 4x4\n");
479  fprintf(f1," C R XY[cell] XYZ[cm] | XYZ at HCAL | C R distance XYZ of closest Hcal 4x4\n");
480 
481  fprintf(f3,"static const int EtoHmap[%d][%d][2] = {\n",EY2B2-1,EX2B2-1);
482  fprintf(f5,"static const int EtoH2map[%d][%d][2] = {\n",EY2B2-1,EX2B2-1);
483  fprintf(f6,"static const int EtoH3map[%d][%d][4] = {\n",EY2B2-1,EX2B2-1);
484  for(int ns=1; ns<2; ns++){
485  for(int j=0; j<EY2B2-1; j++){
486  float y=j*2 + EYOFF + 2;
487  fprintf(f3," { ");
488  for(int i=0; i<EX2B2-1; i++){
489  float x=i*2 + EXOFF + 2;
490  exyz[ns][i][j] = mFcsDb->getStarXYZfromColumnRow(ns,x,y);
491  sxyz[ns][i][j] = sf * exyz[ns][i][j];
492  //search for closest hcal 4x4
493  float dmin=999.0;
494  int k,l;
495  for(int ii=0; ii<HX2B2-1; ii++){
496  for(int jj=0; jj<HY2B2-1; jj++){
497  StThreeVectorF diff = sxyz[ns][i][j]-hxyz[ns][ii][jj];
498  float d = diff.mag();
499  //printf("%2d %2d %2d %2d d=%6.2f diff=%8.2f %8.2f %8.2f\n",
500  // i,j,ii,jj,d,diff.x(),diff.y(),diff.z());
501  if(d<dmin){
502  dmin=d;
503  k=ii;
504  l=jj;
505  }
506  }
507  }
508  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",
509  i,j,x,y,
510  exyz[ns][i][j].x(),exyz[ns][i][j].y(),exyz[ns][i][j].z(),
511  sxyz[ns][i][j].x(),sxyz[ns][i][j].y(),sxyz[ns][i][j].z(),
512  k,l,dmin,
513  hxyz[ns][k][l].x(),hxyz[ns][k][l].y(),hxyz[ns][k][l].z());
514  fprintf(f2,"%6.2f ",dmin);
515  if(i==EX2B2-2) fprintf(f2,"\n");
516  fprintf(f3,"{%2d,%2d}",l,k);
517  if(i<EX2B2-2) fprintf(f3,",");
518  if(i==EX2B2-2) fprintf(f3,"}");
519 
520  //hcal 2x2
521  float dmin2=999.0;
522  int kk,ll;
523  for(int ii=0; ii<HX2B2; ii++){
524  for(int jj=0; jj<HY2B2; jj++){
525  StThreeVectorF diff = sxyz[ns][i][j]-Hxyz[ns][ii][jj];
526  float d = diff.mag();
527  //printf("%2d %2d %2d %2d d=%6.2f diff=%8.2f %8.2f %8.2f\n",
528  // i,j,ii,jj,d,diff.x(),diff.y(),diff.z());
529  if(d<dmin2){
530  dmin2=d;
531  kk=ii;
532  ll=jj;
533  }
534  }
535  }
536  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",
537  i,j,x,y,
538  exyz[ns][i][j].x(),exyz[ns][i][j].y(),exyz[ns][i][j].z(),
539  sxyz[ns][i][j].x(),sxyz[ns][i][j].y(),sxyz[ns][i][j].z(),
540  kk,ll,dmin2,
541  Hxyz[ns][kk][ll].x(),Hxyz[ns][kk][ll].y(),Hxyz[ns][kk][ll].z());
542  fprintf(f4,"%6.2f ",dmin2);
543  if(i==EX2B2-2) fprintf(f4,"\n");
544 
545  fprintf(f5,"{%2d,%2d}",ll,kk);
546  if(i<EX2B2-2) fprintf(f5,",");
547  if(i==EX2B2-2) fprintf(f5,"}");
548 
549  int i1=kk-1, i2=kk, j1=ll-1, j2=ll;
550  int r1=-1, r2=-1, r3=-1, r4=-1;
551  if(i1>=0 && i1<HX2B2-1 && j1>=0 && j1<HY2B2-1) r1=i1+j1*(HX2B2-1);
552  if(i2>=0 && i2<HX2B2-1 && j1>=0 && j1<HY2B2-1) r2=i2+j1*(HX2B2-1);
553  if(i1>=0 && i1<HX2B2-1 && j2>=0 && j2<HY2B2-1) r3=i1+j2*(HX2B2-1);
554  if(i2>=0 && i2<HX2B2-1 && j2>=0 && j1<HY2B2-1) r4=i2+j2*(HX2B2-1);
555  fprintf(f6,"{%2d,%2d,%2d,%2d}",r1,r2,r3,r4);
556  if(i<EX2B2-2) fprintf(f6,",");
557  if(i==EX2B2-2) fprintf(f6,"}");
558  }
559  if(j<EY2B2-2) fprintf(f3,",\n");
560  if(j==EY2B2-2) fprintf(f3,"\n");
561  if(j<EY2B2-2) fprintf(f5,",\n");
562  if(j==EY2B2-2) fprintf(f5,"\n");
563  if(j<EY2B2-2) fprintf(f6,",\n");
564  if(j==EY2B2-2) fprintf(f6,"\n");
565  }
566  fprintf(f3,"}\n");
567  fprintf(f5,"}\n");
568  fprintf(f6,"}\n");
569  }
570  fclose(f1);
571  fclose(f2);
572  fclose(f3);
573  fclose(f4);
574  fclose(f5);
575  fclose(f6);
576  return;
577 }
578 
579 
580 void StFcsTriggerSimMaker::printJP(){
581  //printout ecal 4x4
582  FILE* f1=fopen("EHJP.txt","w");
583  FILE* f2=fopen("EHJPdist.txt","w");
584 
585  enum {EXOFF=1,EYOFF=1,HXOFF=0,HYOFF=2};
586  StThreeVectorF exyz[2][3][5];
587  StThreeVectorF hxyz[2][3][5];
588  StThreeVectorF sxyz[2][3][5]; //ecal 4x4 extraprated at hcal
589  StThreeVectorF eoff = mFcsDb->getDetectorOffset(1);
590  StThreeVectorF hoff = mFcsDb->getDetectorOffset(3);
591  float esmx=mFcsDb->getShowerMaxZ(1);
592  float hsmx=mFcsDb->getShowerMaxZ(3);
593  float r1 = sqrt(eoff.x()*eoff.x()+eoff.z()*eoff.z()); //distance from IP to ecal surface
594  float r2 = r1 + esmx; //distance from IP to ecal SMax
595  float r3 = sqrt(hoff.x()*hoff.x()+hoff.z()*hoff.z()); //distance from IP to hcal surface
596  float r4 = r3 + hsmx; //distance from IP to hcal SMax
597  float sf = r4/r2; // scale factor for extraporating ecal point to hcal plane
598  fprintf(f1,"Distannce from IP to EcalSmax=%8.3f HcalSMax=%8.3f Ratio=%6.3f\n",r2,r4,sf);
599  fprintf(f1,"\nHcal 8x8\n");
600  fprintf(f1," C R XY[cell] XYZ[cm]\n");
601  for(int ns=1; ns<2; ns++){
602  for(int j=0; j<5; j++){
603  float y=j*2 + HYOFF + 4;
604  for(int i=0; i<3; i++){
605  float x=i*2 + HXOFF + 4;
606  hxyz[ns][i][j] = mFcsDb->getStarXYZfromColumnRow(ns+2,x,y);
607  fprintf(f1,"H %2d %2d %4.1f %4.1f %8.2f %8.2f %8.2f\n",
608  i,j,x,y,hxyz[ns][i][j].x(),hxyz[ns][i][j].y(),hxyz[ns][i][j].z());
609  }
610  }
611  }
612  fprintf(f1,"\nEcal 12x16\n");
613  fprintf(f1," C R XY[cell] XYZ[cm] | XYZ at HCAL | distance XYZ of Hcal 8x8\n");
614  for(int ns=1; ns<2; ns++){
615  for(int j=0; j<5; j++){
616  float y=j*4 + EYOFF + 8;
617  for(int i=0; i<3; i++){
618  float x=i*4 + EXOFF + 6;
619  exyz[ns][i][j] = mFcsDb->getStarXYZfromColumnRow(ns,x,y);
620  sxyz[ns][i][j] = sf * exyz[ns][i][j];
621  StThreeVectorF diff = sxyz[ns][i][j]-hxyz[ns][i][j];
622  float d = diff.mag();
623  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",
624  i,j,x,y,
625  exyz[ns][i][j].x(),exyz[ns][i][j].y(),exyz[ns][i][j].z(),
626  sxyz[ns][i][j].x(),sxyz[ns][i][j].y(),sxyz[ns][i][j].z(),
627  d,
628  hxyz[ns][i][j].x(),hxyz[ns][i][j].y(),hxyz[ns][i][j].z());
629  fprintf(f2,"%6.2f ",d);
630  if(i==2) fprintf(f2,"\n");
631  }
632  }
633  }
634  fclose(f1);
635  fclose(f2);
636  return;
637 }
638 
639 void StFcsTriggerSimMaker::readThresholdFile(){
640  LOG_INFO << Form("Reading Thresholds from %s",mThresholdFile)<<endm;
641  FILE* F=fopen(mThresholdFile,"r");
642  if(F == NULL){
643  LOG_ERROR << Form("Could not open %s",mThresholdFile)<<endm;
644  return;
645  }
646  char f[10],name[100];
647  int i,id,v;
648  while(fscanf(F,"%s %d %d %s %d",f, &i, &id, name, &v) != EOF){
649  if(f[0] == '#') continue;
650  TString trg(name);
651  printf("Reading Threshold %s %d\n",trg.Data(),v);
652  if (trg=="FCS_PHTTHR") mTrgSim->PHTTHR=v;
653  else if(trg=="FCS_EM-HERATIO-THR") mTrgSim->EM_HERATIO_THR=v;
654  else if(trg=="FCS_HAD-HERATIO-THR") mTrgSim->HAD_HERATIO_THR=v;
655  else if(trg=="FCS_EMTHR2") mTrgSim->EMTHR2=v;
656  else if(trg=="FCS_EMTHR1") mTrgSim->EMTHR1=v;
657  else if(trg=="FCS_EMTHR0") mTrgSim->EMTHR0=v;
658  else if(trg=="FCS_ELETHR2") mTrgSim->ELETHR2=v;
659  else if(trg=="FCS_ELETHR1") mTrgSim->ELETHR1=v;
660  else if(trg=="FCS_ELETHR0") mTrgSim->ELETHR0=v;
661  else if(trg=="FCS_HADTHR2") mTrgSim->HADTHR2=v;
662  else if(trg=="FCS_HADTHR1") mTrgSim->HADTHR1=v;
663  else if(trg=="FCS_HADTHR0") mTrgSim->HADTHR0=v;
664  else if(trg=="FCS_JPATHR2") mTrgSim->JPATHR2=v;
665  else if(trg=="FCS_JPATHR1") mTrgSim->JPATHR2=v;
666  else if(trg=="FCS_JPATHR0") mTrgSim->JPATHR0=v;
667  else if(trg=="FCS_JPBCTHR2") mTrgSim->JPBCTHR2=v;
668  else if(trg=="FCS_JPBCTHR1") mTrgSim->JPBCTHR1=v;
669  else if(trg=="FCS_JPBCTHR0") mTrgSim->JPBCTHR0=v;
670  else if(trg=="FCS_JPBCTHRD") mTrgSim->JPBCTHRD=v;
671  else if(trg=="FCS_JPDETHR2") mTrgSim->JPDETHR2=v;
672  else if(trg=="FCS_JPDETHR1") mTrgSim->JPDETHR1=v;
673  else if(trg=="FCS_JPDETHR0") mTrgSim->JPDETHR0=v;
674  else if(trg=="FCS_JPDETHRD") mTrgSim->JPDETHRD=v;
675  else if(trg=="FCS_EHTTHR") mTrgSim->EHTTHR=v;
676  else if(trg=="FCS_HHTTHR") mTrgSim->HHTTHR=v;
677  else if(trg=="FCS_ETOTTHR") mTrgSim->ETOTTHR=v;
678  else if(trg=="FCS_HTOTTHR") mTrgSim->HTOTTHR=v;
679  else{
680  printf("No Threshold called %s %d\n",trg.Data(),v);
681  }
682  }
683  fclose(F);
684 }
685 
686 //read from offline copy of online runlog DB
687 void StFcsTriggerSimMaker::readThresholdDb(){
688  //to be implemented before run22 online DB moves to offline
689 }
690 
691 template<typename T> void StFcsTriggerSimMaker::feedADC(T* hit, int ns, int ehp, uint16_t data_array[]){
692 
693  unsigned short dep = hit->dep();
694  unsigned short ch = hit->channel();
695 
696 // printf("ns=%1d ehp=%1d dep=%2d ch=%2d adc=%4d sim=%d\n",ns,ehp,dep,ch,hit->adc(0),mSimMode);
697  fcs_trg_sim_adc[ns][ehp][dep][ch] = hit->adc(0);
698  if(mSimMode==0){
699  int ntb=hit->nTimeBin();
700  for(int t=0; t<ntb; t++){
701  int tb = hit->timebin(t);
702  if(tb>=mTrgTimebin-3 && tb<=mTrgTimebin+4){
703  data_array[tb-mTrgTimebin+3] = hit->adc(t);
704 // printf("tb=%3d i=%2d adc=%4d\n",tb,tb-mTrgTimebin+3,hit->adc(t));
705  }
706  }
707  mTrgSim->fill_event(ehp,ns,dep,ch,data_array,8) ;
708  }else{
709  data_array[1] = hit->adc(0)-1; //removing 1 to add at tb6
710  data_array[6] = 1; //add this so tb6>tb7
711  mTrgSim->fill_event(ehp,ns,dep,ch,data_array,8) ;
712  }
713  if(mFile) fprintf(mFile,"%2d %2d %2d %2d %5d\n",ns,ehp,dep,ch,hit->adc(0));
714 
715 }
StMuDst::muFcsCollection
static StMuFcsCollection * muFcsCollection()
returns pointer to current StMuFcsCollection
Definition: StMuDst.h:395
StFcsDb::getGainCorrection
float getGainCorrection(int det, int id) const
get the gain for the channel for 8 timebin sum
Definition: StFcsDb.cxx:910
StFcsDb::northSouth
int northSouth(int det) const
Ecal=0, Hcal=1, Pres=2.
Definition: StFcsDb.cxx:388
StFcsDb::getDetectorOffset
StThreeVectorD getDetectorOffset(int det, double zdepth=-1) const
Utility functions related to DetectorPosition.
Definition: StFcsDb.cxx:536
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:821
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:381
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:426
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:946
Generated by   doxygen 1.8.5