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