StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StRHICfRecoEnergy.cxx
1 #include "StRHICfRecoEnergy.h"
2 
3 StRHICfRecoEnergy::StRHICfRecoEnergy()
4 {
5  init();
6 }
7 
8 StRHICfRecoEnergy::~StRHICfRecoEnergy()
9 {
10  for(int it=0; it<kRHICfNtower; it++){
11  delete mLeakOutTableNeutron[it];
12  mLeakOutTableNeutron[it] = nullptr;
13  for(int ip=0; ip<kRHICfNplate; ip++){
14  delete mLeakInTablePhoton[it][ip];
15  delete mLeakOutTablePhoton[it][ip];
16  mLeakInTablePhoton[it][ip] = nullptr;
17  mLeakOutTablePhoton[it][ip] = nullptr;
18  }
19  }
20 }
21 
22 void StRHICfRecoEnergy::init()
23 {
24  mRunType = 999;
25  std::fill(mRecoHitNum, mRecoHitNum+kRHICfNtower, -999);
26 
27  memset(mPlateE, 0, sizeof(mPlateE));
28  memset(mPlateECorr, 0, sizeof(mPlateECorr));
29  memset(mOverlap, 0, sizeof(mOverlap));
30  memset(mRecoHitPos, 0, sizeof(mRecoHitPos));
31  memset(mSumEnergy, 0, sizeof(mSumEnergy));
32  memset(mResultEnergy, 0, sizeof(mResultEnergy));
33  memset(mEnergyRatio, 0, sizeof(mEnergyRatio));
34  memset(mMultiHitPos, 0, sizeof(mMultiHitPos));
35  memset(mMultiPeakHeight, 0, sizeof(mMultiPeakHeight));
36 }
37 
38 bool StRHICfRecoEnergy::calculate()
39 {
40  for(int it=0; it<kRHICfNtower; it++){
41  mSumEnergy[it][0] = getSumEnergy(it, 0, 15, mPlateE[it]);
42  mSumEnergy[it][1] = getSumEnergy(it, 1, 11, mPlateE[it]);
43  double mParticleEnergy = 0.;
44  recoPhotonEnergySimple(it, mParticleEnergy);
45 
46  if(mParticleEnergy < mEnergyThreshold){setResultHitNum(it, 0);}
47  }
48 
49  // Type-I pi0.
50  if(getResultHitNum(0)>0 && getResultHitNum(1)>0){
51  recoPhotonEnergyDouble();
52 
53  for(int it=0; it<2; it++){
54  if(mResultEnergy[it][kRHICfPhoton] < mEnergyThreshold){setResultHitNum(it, 0);}
55  }
56  }
57 
58  // Single photon and Type-II pi0.
59  if(getResultHitNum(0)>0 && getResultHitNum(1)==0){ // for TS tower
60  correctLightYield();
61  recoPhotonEnergySingle(0);
62  }
63  if(getResultHitNum(0)==0 && getResultHitNum(1)>0){ // for TL tower
64  correctLightYield();
65  recoPhotonEnergySingle(1);
66  }
67 
68  // Hadron
69  for(int it=0; it<kRHICfNtower; it++){
70  if(getResultHitNum(it) == 0){continue;}
71  recoHadronEnergy(it);
72  }
73 
74  std::cout << "StRHICfRecoEnergy::calculate() -- Done" << std::endl;
75  return kRHICfOk;
76 }
77 
78 void StRHICfRecoEnergy::setRunType(int runType){mRunType = runType;}
79 
80 void StRHICfRecoEnergy::setPlateEnergy(int tower, int layer, double val){mPlateE[tower][layer] = val;}
81 
82 void StRHICfRecoEnergy::setResultHitPos(int tower, int xy, double val){mRecoHitPos[tower][xy] = val;}
83 
84 void StRHICfRecoEnergy::setResultHitNum(int tower, int val){mRecoHitNum[tower] = val;}
85 
86 void StRHICfRecoEnergy::setMultiHitPos(int tower, int layer, int xy, int order, double val){mMultiHitPos[tower][layer][xy][order] = val;}
87 
88 void StRHICfRecoEnergy::setMultiPeakHeight(int tower, int layer, int xy, int order, double val){mMultiPeakHeight[tower][layer][xy][order] = val;}
89 
90 void StRHICfRecoEnergy::setOverlap(int tower, int xy, bool val){mOverlap[tower][xy] = val;}
91 
92 void StRHICfRecoEnergy::setLeakageInTable(int tower, int layer, TH2D* table){mLeakInTablePhoton[tower][layer] = table;}
93 
94 void StRHICfRecoEnergy::setLeakageOutTablePhoton(int tower, int layer, TH2D* table){mLeakOutTablePhoton[tower][layer] = table;}
95 
96 void StRHICfRecoEnergy::setLeakageOutTableNeutron(int tower, TH2D* table){mLeakOutTableNeutron[tower] = table;}
97 
98 int StRHICfRecoEnergy::getResultHitNum(int tower){return mRecoHitNum[tower];}
99 
100 double StRHICfRecoEnergy::getPlateSumEnergy(int tower, bool all){
101  if(all){return mSumEnergy[tower][0];}
102  else{return mSumEnergy[tower][1];}
103  return 0.;
104 }
105 
106 double StRHICfRecoEnergy::getResultEnergy(int tower, int particle){return mResultEnergy[tower][particle];}
107 
108 double StRHICfRecoEnergy::getEnergyRatio(int tower, int order){return mEnergyRatio[tower][order];}
109 
110 void StRHICfRecoEnergy::recoPhotonEnergySimple(int tower, double& energy)
111 {
112  double posX = mRecoHitPos[tower][0];
113  double posY = mRecoHitPos[tower][1];
114  double caliSumEnergy = 0.;
115 
116  for(int ip=1; ip<=11; ip++){
117  double efficiency = getLeakageOutPhoton(tower, ip, posX, posY);
118  if(efficiency<0.1 || efficiency>2.){
119 
120  double gsobarSize = checkGSOBarSize(tower);
121  double edgeSize = (ip<12 ? 1.0 : 2.0);
122 
123  if(posX < edgeSize){posX = edgeSize;}
124  if(posX > gsobarSize-edgeSize){posX = gsobarSize-edgeSize;}
125  if(posY < edgeSize){posY = edgeSize;}
126  if(posY > gsobarSize-edgeSize){posY = gsobarSize-edgeSize;}
127 
128  efficiency = getLeakageOutPhoton(tower, ip, posX, posY);
129  }
130  caliSumEnergy += mPlateE[tower][ip]/efficiency;
131  }
132  energy = getPhotonEnergyConvert(tower, caliSumEnergy);
133 }
134 
135 void StRHICfRecoEnergy::recoPhotonEnergySingle(int tower)
136 {
137  for(int it=0; it<kRHICfNtower; it++){
138  if(it==tower){
139  mSumEnergy[it][0] = getSumEnergy(it, 0, 15, mPlateECorr[it]);
140  mSumEnergy[it][1] = getSumEnergy(it, 1, 11, mPlateECorr[it]);
141  mResultEnergy[it][kRHICfPhoton] = getPhotonEnergyConvert(it, mSumEnergy[it][1]);
142  }
143  else{
144  mSumEnergy[it][0] = getSumEnergy(it, 0, 15, mPlateE[it]);
145  mSumEnergy[it][1] = 0.;
146  mResultEnergy[it][kRHICfPhoton] = getPhotonEnergyConvert(it, mSumEnergy[it][1]);
147  }
148  }
149 }
150 
151 void StRHICfRecoEnergy::recoPhotonEnergyDouble()
152 {
153  double mLeakInEffi[kRHICfNtower][kRHICfNplate];
154  double mLeakOutEffi[kRHICfNtower][kRHICfNplate];
155 
156  memset(mLeakInEffi, 0, sizeof(mLeakInEffi));
157  memset(mLeakOutEffi, 0, sizeof(mLeakOutEffi));
158 
159  for(int it=0; it<kRHICfNtower; it++){
160  double posX = mRecoHitPos[it][0];
161  double posY = mRecoHitPos[it][1];
162  double gsobarSize = checkGSOBarSize(it);
163 
164  if(posX < 0.){posX = 0.1;}
165  if(posY < 0.){posY = 0.1;}
166  if(posX > gsobarSize){posX = gsobarSize-0.1;}
167  if(posY > gsobarSize){posY = gsobarSize-0.1;}
168 
169  for(int ip=0; ip<kRHICfNplate; ip++){
170  mLeakOutEffi[it][ip] = getLeakageOutPhoton(it, ip, posX, posY);
171  mLeakInEffi[it][ip] = getLeakageInPhoton(it, ip, posX, posY);
172 
173  if(ip >= 11){mLeakOutEffi[it][ip] = mLeakOutEffi[it][10];}
174  if(ip >= 9){mLeakInEffi[it][ip] = mLeakInEffi[it][8];}
175  }
176  }
177 
178  for(int ip=0; ip<kRHICfNplate; ip++){
179  double leakInFactor = mLeakInEffi[1][ip] * mLeakInEffi[0][ip] - 1.;
180  mPlateECorr[0][ip] = ((mLeakInEffi[1][ip] * mPlateE[1][ip] - mPlateE[0][ip]) / leakInFactor) / mLeakOutEffi[0][ip];
181  mPlateECorr[1][ip] = ((mLeakInEffi[0][ip] * mPlateE[0][ip] - mPlateE[1][ip]) / leakInFactor) / mLeakOutEffi[1][ip];
182  }
183 
184  for(int it=0; it<kRHICfNtower; it++){
185  mSumEnergy[it][0] = getSumEnergy(it, 0, 15, mPlateECorr[it]);
186  mSumEnergy[it][1] = getSumEnergy(it, 1, 11, mPlateECorr[it]);
187  mResultEnergy[it][kRHICfPhoton] = getPhotonEnergyConvert(it, mSumEnergy[it][1]);
188  }
189 }
190 
191 void StRHICfRecoEnergy::recoHadronEnergy(int tower)
192 {
193  double posX = mRecoHitPos[tower][0];
194  double posY = mRecoHitPos[tower][1];
195  double efficiency = getLeakageOutNeutron(tower, posX, posY);
196  double sumEnergy = getSumEnergy(tower, 0, 15, mPlateE[tower]);
197  double sumECorr = sumEnergy/efficiency;
198  mResultEnergy[tower][kRHICfHadron] = getHadronEnergyConvert(tower, sumECorr);
199 }
200 
201 void StRHICfRecoEnergy::correctLightYield()
202 {
203  for(int it=0; it<kRHICfNtower; it++){
204  // single photon
205  if(getResultHitNum(it)==1){
206  double posX = mRecoHitPos[it][0];
207  double posY = mRecoHitPos[it][1];
208  getCalibrationEnergy(it, posX, posY);
209  }
210 
211  // Type-II pi0.
212  if(getResultHitNum(it) > 1){
213  double posX1 = mMultiHitPos[it][0][0][0];
214  double posY1 = mMultiHitPos[it][0][1][0];
215  double posX2 = mMultiHitPos[it][0][0][1];
216  double posY2 = mMultiHitPos[it][0][1][1];
217  double firstRatio = 0.;
218  double secondRatio = 0.;
219 
220  if(!mOverlap[it][0]){firstRatio += mMultiPeakHeight[it][0][0][0];}
221  if(!mOverlap[it][1]){firstRatio += mMultiPeakHeight[it][0][1][0];}
222  if(!mOverlap[it][0]){firstRatio += mMultiPeakHeight[it][1][0][0];}
223  if(!mOverlap[it][1]){firstRatio += mMultiPeakHeight[it][1][1][0];}
224  if(!mOverlap[it][0]){secondRatio += mMultiPeakHeight[it][0][0][1];}
225  if(!mOverlap[it][1]){secondRatio += mMultiPeakHeight[it][0][1][1];}
226  if(!mOverlap[it][0]){secondRatio += mMultiPeakHeight[it][1][0][1];}
227  if(!mOverlap[it][1]){secondRatio += mMultiPeakHeight[it][1][1][1];}
228 
229  mEnergyRatio[it][0] = firstRatio/(firstRatio + secondRatio);
230  mEnergyRatio[it][1] = 1. - mEnergyRatio[it][0];
231 
232  getCalibrationEnergy(it, posX1, posY1, posX2, posY2, mEnergyRatio[it][0], mEnergyRatio[it][1]);
233  }
234  }
235 }
236 
237 void StRHICfRecoEnergy::getCalibrationEnergy(int tower, double posX1, double posY1, double posX2, double posY2, double firstRatio, double secondRatio)
238 {
239  for(int ip=0; ip<kRHICfNplate; ip++){
240  double efficiency1 = getLeakageOutPhoton(tower, ip, posX1, posY1);
241  double efficiency2 = getLeakageOutPhoton(tower, ip, posX2, posY2);
242 
243  if(efficiency1 < 0.1 || efficiency1 > 2.){
244  double gsobarSize = checkGSOBarSize(tower);
245  double edgeSize = (ip<12 ? 1.0 : 2.0);
246 
247  if(posX1 < edgeSize){posX1 = edgeSize;}
248  if(posX1 > gsobarSize-edgeSize){posX1 = gsobarSize-edgeSize;}
249  if(posY1 < edgeSize){posY1 = edgeSize;}
250  if(posY1 > gsobarSize-edgeSize){posY1 = gsobarSize-edgeSize;}
251 
252  efficiency1 = getLeakageOutPhoton(tower, ip, posX1, posY1);
253  }
254  if(efficiency2 < 0.1 || efficiency2 > 2.){
255  double gsobarSize = checkGSOBarSize(tower);
256  double edgeSize = (ip<12 ? 1.0 : 2.0);
257 
258  if(posX2 < edgeSize){posX2 = edgeSize;}
259  if(posX2 > gsobarSize-edgeSize){posX2 = gsobarSize-edgeSize;}
260  if(posY2 < edgeSize){posY2 = edgeSize;}
261  if(posY2 > gsobarSize-edgeSize){posY2 = gsobarSize-edgeSize;}
262 
263  efficiency2 = getLeakageOutPhoton(tower, ip, posX2, posY2);
264  }
265 
266  if(firstRatio <= -100. || secondRatio <= -100.){
267  mPlateECorr[tower][ip] = mPlateE[tower][ip]/efficiency1;
268  }
269  if(firstRatio >= 0. || secondRatio >= 0.){
270  double calibFactor = efficiency1*firstRatio + efficiency2*secondRatio;
271  mPlateECorr[tower][ip] = mPlateE[tower][ip]/calibFactor;
272  }
273  }
274 }
275 
276 double StRHICfRecoEnergy::getSumEnergy(int tower, int startIdx, int endIdx, double* energy)
277 {
278  double mSumEnergy = 0.;
279  for(int ip=startIdx; ip<=endIdx; ip++){
280  if(ip<11){mSumEnergy += energy[ip];}
281  else{mSumEnergy += energy[ip]*2.;}
282  }
283  return mSumEnergy;
284 }
285 
286 double StRHICfRecoEnergy::getPhotonEnergyConvert(int tower, double sumEnergy)
287 {
288  if(mRunType==kRHICfTL){
289  if(tower==0) return 31.1323*sumEnergy+1.42462;
290  if(tower==1) return 27.3429*sumEnergy+1.11522;
291  }
292  if(mRunType==kRHICfTS){
293  if(tower==0) return 31.1666*sumEnergy+1.39715;
294  if(tower==1) return 27.3389*sumEnergy+1.12892;
295  }
296  if(mRunType==kRHICfTOP){
297  if(tower==0) return 31.1749*sumEnergy+1.41172;
298  if(tower==1) return 27.3434*sumEnergy+1.12205;
299  }
300  return 0;
301 }
302 
303 double StRHICfRecoEnergy::getHadronEnergyConvert(int tower, double sumEnergy)
304 {
305  if(tower==0) return 87.383*sumEnergy-11.1752;
306  if(tower==1) return 68.472*sumEnergy-15.0983;
307  return 0;
308 }
309 
310 double StRHICfRecoEnergy::getLeakageInPhoton(int tower, int layer, double x, double y)
311 {
312  return mLeakInTablePhoton[tower][layer]->Interpolate(x, y);
313 }
314 
315 double StRHICfRecoEnergy::getLeakageOutPhoton(int tower, int layer, double x, double y)
316 {
317  return mLeakOutTablePhoton[tower][layer]->Interpolate(x, y);
318 }
319 
320 double StRHICfRecoEnergy::getLeakageOutNeutron(int tower, double x, double y)
321 {
322  return mLeakOutTableNeutron[tower]->Interpolate(x, y);
323 }
Generated by   doxygen 1.8.5