StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StFcsFastSimulatorMaker.cxx
1 // $Id: StFcsFastSimulatorMaker.cxx,v 1.2 2021/03/30 13:40:09 akio Exp $
2 //
3 // $Log: StFcsFastSimulatorMaker.cxx,v $
4 // Revision 1.2 2021/03/30 13:40:09 akio
5 // FCS code after peer review and moved from $CVSROOT/offline/upgrades/akio
6 //
7 // Revision 1.10 2021/02/25 21:54:41 akio
8 // Int_t -> int
9 //
10 // Revision 1.9 2021/02/25 19:25:48 akio
11 // Code modified for STAR code review
12 //
13 // Revision 1.8 2021/02/23 16:25:50 akio
14 // Modification to attend comments from STAR code review (Jason)
15 //
16 // Revision 1.7 2020/08/27 22:08:09 akio
17 // fix a continue bug found by Ting Lin
18 //
19 // Revision 1.6 2020/05/29 18:51:02 akio
20 // adding EPD g2t reading as PRES
21 //
22 // Revision 1.5 2019/10/23 17:15:42 akio
23 // *** empty log message ***
24 //
25 // Revision 1.4 2019/07/22 18:56:41 akio
26 // Added LeakyHcal option 2 and 3 for 2d light collection efficiency parametrization
27 //
28 // Revision 1.3 2019/06/26 18:02:28 akio
29 // MOdify StFcsHit creating to match StEvent update
30 //
31 // Revision 1.2 2019/05/16 16:11:55 akio
32 // Adding leaky hcal option
33 //
34 // Revision 1.1 2018/11/14 16:50:16 akio
35 // FCS codes in offline/upgrade/akio
36 //
37 //
38 // Implementation of StFcsFastSimulatorMaker, the FCS fast simulator
39 //
40 
41 #include "StFcsFastSimulatorMaker/StFcsFastSimulatorMaker.h"
42 
43 #include "St_base/StMessMgr.h"
44 #include "StEvent/StEvent.h"
45 #include "StEvent/StFcsCollection.h"
46 #include "StEvent/StFcsHit.h"
47 #include "tables/St_g2t_emc_hit_Table.h"
48 #include "tables/St_g2t_hca_hit_Table.h"
49 #include "tables/St_g2t_epd_hit_Table.h"
50 #include "StFcsDbMaker/StFcsDb.h"
51 
52 namespace{
53  static const char name[kFcsEHP][4]={"wca","hca","pre"};
54 }
55 
56 StFcsFastSimulatorMaker::StFcsFastSimulatorMaker(const char* name) : StMaker(name) {
57  setLeakyHcal(0);
58  setHcalZDepEff(0);
59 }
60 
61 int StFcsFastSimulatorMaker::Init() {
62  mFcsDb = static_cast<StFcsDb*>(GetDataSet("fcsDb"));
63  if (!mFcsDb) {
64  LOG_ERROR << "StFcsFastSimulatorMaker initializing failed due to no StFcsDb" << endm;
65  return kStErr;
66  }
67 
68  memset(mEcalMap,0,sizeof(mEcalMap));
69  memset(mHcalMap,0,sizeof(mHcalMap));
70  memset(mPresMap,0,sizeof(mPresMap));
71  return kStOK;
72 }
73 
74 void StFcsFastSimulatorMaker::Clear(Option_t *option){
75  memset(mEcalMap,0,sizeof(mEcalMap));
76  memset(mHcalMap,0,sizeof(mHcalMap));
77  memset(mPresMap,0,sizeof(mPresMap));
78  StMaker::Clear(option);
79  return;
80 }
81 
83  LOG_DEBUG << "StFcsFastSimulatorMaker::Make" << endm;
84 
85  // Get the existing StEvent, or add one if it doesn't exist.
86  StEvent* event = static_cast<StEvent*>(GetDataSet("StEvent"));
87  if (!event) {
88  event = new StEvent;
89  AddData(event);
90  LOG_DEBUG << "Creating StEvent" << endm;
91  }
92 
93  // Add an FCS collection to the event if one does not already exist.
94  if (!event->fcsCollection()) {
95  event->setFcsCollection(new StFcsCollection);
96  LOG_DEBUG << "Creating StFcsCollection" << endm;
97  }
98 
99  fillStEvent(event);
100  return kStOk;
101 }
102 
103 /* Fill an event with StFcsHits. */
104 void StFcsFastSimulatorMaker::fillStEvent(StEvent* event) {
105  StFcsCollection * fcscollection = event->fcsCollection();
106  int ehp;
107  int ng2thit[kFcsEHP]={};
108  StPtrVecFcsHit hits; //temp storage for hits
109  int leakyHcal = IAttr("FcsLeakyHcal");
110  int hcalZdepEff = IAttr("FcsHcalZdepEff");
111 
112  // Read the g2t table for FCS Wcal
113  ehp=0;
114  St_g2t_emc_hit* hitTable = static_cast<St_g2t_emc_hit*>(GetDataSet(Form("g2t_%3s_hit",name[ehp])));
115  if (!hitTable) {
116  LOG_INFO << Form("g2t_%3s_hit table is empty",name[ehp]) << endm;
117  }else{
118  const int nHits = hitTable->GetNRows();
119  ng2thit[ehp]=nHits;
120  LOG_INFO << Form("g2t_%s_hit table has %d hit",name[ehp],nHits) << endm;
121  if(nHits>0){
122  const g2t_emc_hit_st* hit = hitTable->GetTable();
123  if(!hit){
124  LOG_INFO << Form("g2t_%3s_hit GetTable failed",name[ehp]) << endm;
125  }else{
126  for (int i=0; i < nHits; ++i) {
127  if (!hit) {hit++; continue;}
128  const int ns = hit->volume_id / 1000 - 1;
129  const int id = hit->volume_id % 1000 - 1;
130  const int det = mFcsDb->detectorId(ehp,ns);
131  if(det<0 || det>=kFcsNDet || id<0 || id>=kFcsEcalMaxId){
132  LOG_WARN << Form("ECAL det=%1d id=%3d volid=%5d e=%f out of range (%d)",
133  det,id,hit->volume_id,hit->de,kFcsMaxId) << endm;
134  hit++;
135  continue;
136  }else if(GetDebug()){
137  LOG_INFO << Form("ECAL det=%1d id=%3d volid=%4d e=%f",
138  det,id,hit->volume_id,hit->de) << endm;
139  }
140  float de = hit->de;
141  StFcsHit* fcshit=0;
142  if(mEcalMap[ns][id]==0){ // New hit
143  int ehp=0, rns=0, crt=0, sub=0, dep=0, ch=0;
144  mFcsDb->getDepfromId(det, id, ehp, rns, crt, sub, dep, ch);
145  fcshit = new StFcsHit(1, det, id, rns, ehp, dep, ch, de);
146  hits.push_back(fcshit);
147  mEcalMap[ns][id]=fcshit;
148  }else{ // Adding energy to old hit
149  fcshit = mEcalMap[ns][id];
150  fcshit->setEnergy(fcshit->energy() + de);
151  }
152  fcshit->addGeantTrack(hit->track_p,de);
153  hit++;
154  }
155  }
156  }
157  }
158 
159  // Read the g2t table for FCS Hcal
160  ehp=1;
161  St_g2t_hca_hit* hitTable_h = static_cast<St_g2t_hca_hit*>(GetDataSet(Form("g2t_%3s_hit",name[ehp])));
162  if (!hitTable_h) {
163  LOG_INFO << Form("g2t_%3s_hit table is empty",name[ehp]) << endm;
164  }else{
165  const int nHits = hitTable_h->GetNRows();
166  ng2thit[ehp]=nHits;
167  LOG_INFO << Form("g2t_%s_hit table has %d hit",name[ehp],nHits) << endm;
168  if(nHits>0){
169  const g2t_hca_hit_st* hit = hitTable_h->GetTable();
170  if(!hit){
171  LOG_INFO << Form("g2t_%3s_hit GetTable failed",name[ehp]) << endm;
172  }else{
173  for (int i=0; i < nHits; i++) {
174  if (!hit) {hit++; continue;}
175  const int ns = hit->volume_id / 1000 - 1;
176  const int id = hit->volume_id % 1000 - 1;
177  const int det = mFcsDb->detectorId(ehp,ns);
178  if(det<0 || det>=kFcsNDet || id<0 || id>=kFcsHcalMaxId){
179  LOG_WARN << Form("HCAL det=%d id=%d volid=%5d e=%f out of range (%d)",
180  det,id,hit->volume_id,hit->de,kFcsMaxId) << endm;
181  hit++;
182  continue;
183  }else if(GetDebug()){
184  LOG_INFO << Form("HCAL det=%d id=%d volid=%5d e=%f",
185  det,id,hit->volume_id,hit->de) << endm;
186  }
187  StFcsHit* fcshit=0;
188  int ehp=0, rns=0, crt=0, sub=0, dep=0, ch=0;
189  mFcsDb->getDepfromId(det, id, ehp, rns, crt, sub, dep, ch);
190  if(leakyHcal==0 || leakyHcal==2){
191  float de;
192  de = hit->de;
193  //if(leakyHcal==2) de = hit->de2;
194  if(hcalZdepEff==1) de = hit->deA;
195  if(hcalZdepEff==2) de = hit->deB;
196  if(mHcalMap[ns][id]==0){ // New hit
197  fcshit = new StFcsHit(1, det, id, rns, ehp, dep, ch, de);
198  hits.push_back(fcshit);
199  mHcalMap[ns][id]=fcshit;
200  }else{ // Adding energy to old hit
201  fcshit = mHcalMap[ns][id];
202  fcshit->setEnergy(fcshit->energy() + de);
203  }
204  fcshit->addGeantTrack(hit->track_p,de);
205  hit++;
206  }else{ //leaky hcal with up to 4 WLSP getting lights from a tower
207  float de[4];
208  if(leakyHcal==1){
209  de[0] = hit->deA;
210  de[1] = hit->deB;
211  de[2] = hit->deC;
212  de[3] = hit->deD;
213  }else{ //tested and turned to be minor. Not implemented in official XML
214  //de[0] = hit->de2A;
215  //de[1] = hit->de2B;
216  //de[2] = hit->de2C;
217  //de[3] = hit->de2D;
218  }
219  int col = mFcsDb->getColumnNumber(det,id); // col goes 1 ~ ncol
220  int ncol= mFcsDb->nColumn(det);
221  for(int j=0; j<4; j++){
222  int id2;
223  int jj = j-2;; //jj goes from -2 ~ +1
224  if(col==1 && jj<0) {id2=id;} // if col=1, add leaked light back to col=1
225  else if(col==2 && jj==-2) {id2=id-1;} // if col=2, add leaked light back to col=1
226  else if(col==ncol && jj==1) {id2=id;} // if col=ncol, add leaked light back to col=ncol
227  else {id2=id+jj;} // add leaked lights to its neighbors
228  mFcsDb-> getDepfromId(det, id2, ehp, rns, crt, sub, dep, ch);
229  if(mHcalMap[ns][id2]==0){ // New hit
230  int ehp=0, rns=0, crt=0, sub=0, dep=0, ch=0;
231  mFcsDb-> getDepfromId(det, id2, ehp, rns, crt, sub, dep, ch);
232  fcshit = new StFcsHit(1, det, id2, rns, ehp, dep, ch, de[j]);
233  hits.push_back(fcshit);
234  mHcalMap[ns][id2]=fcshit;
235  }else{ // Adding energy to old hit
236  fcshit = mHcalMap[ns][id2];
237  fcshit->setEnergy(fcshit->energy() + de[j]);
238  }
239  fcshit->addGeantTrack(hit->track_p,de[j]);
240  }
241  hit++;
242  }
243  }
244  }
245  }
246  }
247 
248  // Read the g2t table for FCS Pres == EPD now
249  ehp=2;
250  St_g2t_epd_hit* hitTable_p = static_cast<St_g2t_epd_hit*>(GetDataSet("g2t_epd_hit"));
251  if (!hitTable_p) {
252  LOG_INFO << Form("g2t_epd_hit table is empty") << endm;
253  }else{
254  const int nHits = hitTable_p->GetNRows();
255  ng2thit[ehp]=nHits;
256  LOG_INFO << Form("g2t_epd_hit table has %d hit",nHits) << endm;
257  if(nHits>0){
258  const g2t_epd_hit_st* hit = hitTable_p->GetTable();
259  if(!hit){
260  LOG_INFO << Form("g2t_epd_hit GetTable failed") << endm;
261  }else{
262  for (int i=0; i < nHits; ++i) {
263  if (!hit) {hit++; continue;}
264  const int volume_id = hit->volume_id;
265  const int ew = volume_id/100000;
266  const int pp = (volume_id%100000)/1000;
267  int tt = (volume_id%1000)/10;
268  const float de = hit->de * 1000.0;
269  if(ew==1) {hit++; continue;} //west side only
270  //Hack! reverse TT Even/Odd. To be removed when EPD XML file fixed
271  if(tt>1){
272  if(tt%2==0) {tt+=1;}
273  else {tt-=1;}
274  }
275  int det,id,ehp,ns,crt,slt,dep,ch;
276  mFcsDb->getIdfromEPD(pp,tt,det,id);
277  mFcsDb->getDepfromId(det,id,ehp,ns,crt,slt,dep,ch);
278  if(det<0 || det>=kFcsNDet || id<0 || id>=kFcsPresMaxId){
279  LOG_WARN << Form("Pres det=%1d id=%3d volid=%4d e=%f10.6 id out of range (%d)",
280  det,id,hit->volume_id,1000*hit->de,kFcsPresId) << endm;
281  hit++;
282  continue;
283  }else if(GetDebug()){
284  LOG_INFO << Form("Pres det=%1d id=%3d volid=%4d e=%10.6f",
285  det,id,hit->volume_id,1000*hit->de) << endm;
286  }
287  StFcsHit* fcshit=0;
288  if(mPresMap[ns][id]==0){ // New hit
289  fcshit = new StFcsHit(1, det, id, ns, ehp, dep, ch, de);
290  hits.push_back(fcshit);
291  mPresMap[ns][id]=fcshit;
292  }else{ // Adding energy to old hit
293  fcshit = mPresMap[ns][id];
294  fcshit->setEnergy(fcshit->energy() + de);
295  }
296  hit++;
297  fcshit->addGeantTrack(hit->track_p,de);
298  }
299  }
300  }
301  }
302 
303  int nhittot=hits.size();
304  int zs=0;
305  int nhit[kFcsEHP]={};
306  float etot[kFcsEHP]={};
307  // Loop over hits and digitize & push to StEvent if it survives ZS
308  for(int i=0; i<nhittot; i++){
309  const int det = hits[i]->detectorId();
310  const int id = hits[i]->id();
311  float de = hits[i]->energy();
312  float sf = mFcsDb->getSamplingFraction(det);
313  float gain= mFcsDb->getGain(det, id);
314  float corr= mFcsDb->getGainCorrection(det, id);
315  int adc = static_cast<int>(de / (sf * gain * corr));
316  adc = std::min(adc, 4095*8); // Cap maximum ADC =12bit * 8 tim
317  zs = mFcsDb->getZeroSuppression(det);
318  if(GetDebug()) LOG_INFO << Form("Det=%1d id=%3d dE=%8.3f SF=%6.3f gain=%6.3f corr=%6.3f ADC=%4d ZS=%2d digiE=%8.3f",
319  det,id,de,sf,gain,corr,adc,zs,adc*gain*corr) << endm;
320  if(adc>zs){ // zero suppress low ADC
321  float digi_energy = adc * gain * corr;
322  int ehp = mFcsDb->ecalHcalPres(det);
323  hits[i]->setAdc(0,adc);
324  hits[i]->setEnergy(digi_energy);
325  fcscollection->addHit(det,hits[i]);
326  etot[ehp] += digi_energy;
327  nhit[ehp]++;
328  }else{ // just delete hits bellow ZS threshold
329  delete hits[i];
330  }
331  }
332 
333  if(GetDebug()>0){
334  for(int ehp=0; ehp<kFcsEHP; ehp++){
335  LOG_INFO << Form("%s Found %d g2t hits in %d cells, created %d hits with ADC>ZS(%d) and Etot=%8.3f",
336  name[ehp],ng2thit[ehp],nhit[ehp],
337  fcscollection->numberOfHits(ehp*2) + fcscollection->numberOfHits(ehp*2+1),
338  zs, etot[ehp])
339  << endm;
340  }
341  if(GetDebug()>1) fcscollection->print();
342  }
343 }
int getColumnNumber(int det, int id) const
get the row number for the channel
Definition: StFcsDb.cxx:490
void getIdfromEPD(int pp, int tt, int &det, int &id)
Map header.
Definition: StFcsDb.cxx:1467
int getZeroSuppression(int det) const
fcsGain/GainCorrection related
Definition: StFcsDb.cxx:920
virtual void AddData(TDataSet *data, const char *dir=".data")
User methods.
Definition: StMaker.cxx:332
float getGainCorrection(int det, int id) const
get the gain for the channel for 8 timebin sum
Definition: StFcsDb.cxx:952
virtual void Clear(Option_t *option="")
User defined functions.
Definition: StMaker.cxx:634
int detectorId(int eh, int ns) const
6
Definition: StFcsDb.cxx:418
float getSamplingFraction(int det) const
get zero suppression threshold
Definition: StFcsDb.cxx:913
int ecalHcalPres(int det) const
Ecal North=0, Ecal South=1, Hcal North=2, Hcal South=3, Pres=4/5.
Definition: StFcsDb.cxx:423
float getGain(int det, int id) const
get sampling fraction
Definition: StFcsDb.cxx:928
void Clear(Option_t *option="")
User defined functions.
Definition: Stypes.h:40
int nColumn(int det) const
number of rows
Definition: StFcsDb.cxx:451
Definition: Stypes.h:44
Definition: Stypes.h:41
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