StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StRHICfRecoPos.cxx
1 #include "StRHICfRecoPos.h"
2 
3 StRHICfRecoPos::StRHICfRecoPos()
4 {
5  memset(mGSOBarSTable, 0, sizeof(mGSOBarSTable));
6  memset(mGSOBarTTable, 0, sizeof(mGSOBarTTable));
7 
8  for(int it=0; it<kRHICfNtower; it++){
9  for(int il=0; il<kRHICfNlayer; il++){
10  for(int ixy=0; ixy<kRHICfNxy; ixy++){
11  int gsobarSize = checkGSOBarSize(it);
12 
13  mGSOBarHist[it][il][ixy] = new TH1D("", "", gsobarSize, 0, gsobarSize);
14  mGSOBarHistExpend[it][il][ixy] = new TH1D("", "", gsobarSize+4, -2, gsobarSize+2);
15  mGSOBarGraph[it][il][ixy] = new TGraphErrors(gsobarSize);
16 
17  mSingleFit[it][il][ixy] = new TF1("SingleFit", this, &StRHICfRecoPos::getLorenzianSingle, 0., double(gsobarSize), 6);
18  mMultiFit[it][il][ixy] = new TF1("MultiFit", this, &StRHICfRecoPos::getLorenzianMulti, 0., double(gsobarSize), 11);
19  mEachFit[it][il][ixy][0] = new TF1("EachFit1", this, &StRHICfRecoPos::getLorenzianSingle, 0., double(gsobarSize), 6);
20  mEachFit[it][il][ixy][1] = new TF1("EachFit2", this, &StRHICfRecoPos::getLorenzianSingle, 0., double(gsobarSize), 6);
21 
22  mSpectrum[it][il][ixy] = new TSpectrum();
23  }
24  }
25  }
26 }
27 
28 StRHICfRecoPos::~StRHICfRecoPos()
29 {
30  for(int it=0; it<kRHICfNtower; it++){
31  for(int il=0; il<kRHICfNlayer; il++){
32  for(int ixy=0; ixy<kRHICfNxy; ixy++){
33  delete mGSOBarHist[it][il][ixy];
34  delete mGSOBarHistExpend[it][il][ixy];
35  delete mGSOBarGraph[it][il][ixy];
36  delete mSingleFit[it][il][ixy];
37  delete mMultiFit[it][il][ixy];
38  delete mEachFit[it][il][ixy][0];
39  delete mEachFit[it][il][ixy][1];
40 
41  mGSOBarHist[it][il][ixy] = nullptr;
42  mGSOBarHistExpend[it][il][ixy] = nullptr;
43  mGSOBarGraph[it][il][ixy] = nullptr;
44  mSingleFit[it][il][ixy] = nullptr;
45  mMultiFit[it][il][ixy] = nullptr;
46  mEachFit[it][il][ixy][0] = nullptr;
47  mEachFit[it][il][ixy][1] = nullptr;
48  }
49  }
50  }
51 }
52 
53 void StRHICfRecoPos::init()
54 {
55  memset(mWorthy, false, sizeof(mWorthy));
56  memset(mOverlap, false, sizeof(mOverlap));
57  memset(mGSOMaxLayer, 0, sizeof(mGSOMaxLayer));
58  memset(mGSOBarSE, 0, sizeof(mGSOBarSE));
59  memset(mGSOBarTE, 0, sizeof(mGSOBarTE));
60 
61  for(int it=0; it<kRHICfNtower; it++){
62  for(int il=0; il<kRHICfNlayer; il++){
63  for(int ixy=0; ixy<kRHICfNxy; ixy++){
64  mGSOBarHist[it][il][ixy] -> Reset();
65  mGSOBarHistExpend[it][il][ixy] -> Reset();
66  mGSOBarGraph[it][il][ixy] -> Set(0);
67 
68  for(int p=0; p<11; p++){
69  mMultiFit[it][il][ixy] -> SetParameter(p, 0);
70  if(p < 6){
71  mSingleFit[it][il][ixy] -> SetParameter(p, 0);
72  mEachFit[it][il][ixy][0] -> SetParameter(p, 0);
73  mEachFit[it][il][ixy][1] -> SetParameter(p, 0);
74  }
75  }
76  }
77  }
78  }
79 }
80 
81 void StRHICfRecoPos::setGSOBarEnergy(int tower, int layer, int xy, int bar, Double_t val)
82 {
83  if(tower==0){mGSOBarSE[layer][xy][bar] = val;}
84  if(tower==1){mGSOBarTE[layer][xy][bar] = val;}
85 }
86 
87 void StRHICfRecoPos::setGSOBarTable(int tower, int layer, int xy, int bar, Double_t val)
88 {
89  if(tower==0){mGSOBarSTable[layer][xy][bar] = val;}
90  if(tower==1){mGSOBarTTable[layer][xy][bar] = val;}
91 }
92 
93 bool StRHICfRecoPos::fillData()
94 {
95  // For worthy
96  int tmpWorthy[kRHICfNxy][kRHICfNtower];
97  memset(tmpWorthy, 0, sizeof(tmpWorthy));
98 
99  Double_t pedeRMSGeV = mPedADCRMS/(m1000Vto600VFactor*mSpecialFactor*mAvgConvFactor);
100 
101  for(int it=0; it<kRHICfNtower; it++){
102  for(int il=0; il<kRHICfNlayer; il++){
103  for(int ixy=0; ixy<kRHICfNxy; ixy++){
104  mGSOBarHistExpend[it][il][ixy] -> SetBinContent(1, 0.);
105  mGSOBarHistExpend[it][il][ixy] -> SetBinContent(2, 0.);
106 
107  int gsobarSize = checkGSOBarSize(it);
108  for(int ich=0; ich<gsobarSize; ich++){
109  Double_t gsobarE = 0;
110  Double_t gsobarTable = 0;
111  if(it==0){
112  gsobarE = mGSOBarSE[il][ixy][ich];
113  gsobarTable = mGSOBarSTable[il][ixy][ich];
114  }
115  if(it==1){
116  gsobarE = mGSOBarTE[il][ixy][ich];
117  gsobarTable = mGSOBarTTable[il][ixy][ich];
118  }
119 
120  mGSOBarHist[it][il][ixy] -> SetBinContent(ich+1, gsobarE);
121  mGSOBarHistExpend[it][il][ixy] -> SetBinContent(ich+3, gsobarE);
122  mGSOBarGraph[it][il][ixy] -> SetPoint(ich, gsobarTable, gsobarE);
123  mGSOBarGraph[it][il][ixy] -> SetPointError(ich, 0., mGSOBarMapError * gsobarE + sqrt(pedeRMSGeV));
124 
125  // For worthy
126  double edep = mGSOBarHist[it][il][ixy] -> GetBinContent(ich+1);
127  if(edep > mNoiseThreshold){tmpWorthy[ixy][it] += 1;}
128  }
129  mGSOBarHistExpend[it][il][ixy] -> SetBinContent(gsobarSize+3, 0.);
130  mGSOBarHistExpend[it][il][ixy] -> SetBinContent(gsobarSize+4, 0.);
131  }
132  }
133  }
134  // For worthy
135  for(int it=0; it<kRHICfNtower; it++){if(tmpWorthy[0][it]>0 && tmpWorthy[1][it]>0) mWorthy[it] = true;}
136 
137  // Remove the dead channal
138  mGSOBarGraph[0][3][0]->SetPoint(12, mGSOBarSTable[3][0][12], 0.);
139  mGSOBarGraph[0][3][0]->RemovePoint(12);
140 
141  findMaxLayer();
142 
143  return kRHICfOk;
144 }
145 
146 bool StRHICfRecoPos::calculate()
147 {
148  initSetParamaters();
149  fitting();
150 
151  std::cout << "StRHICfRecoPos::calculate() -- Done " << std::endl;
152  return kRHICfOk;
153 }
154 
155 bool StRHICfRecoPos::separateMultiFit(int tower)
156 {
157  for(int il=0; il<kRHICfNlayer; il++){
158  for(int ixy=0; ixy<kRHICfNxy; ixy++){
159  Double_t parMulti[11];
160  mMultiFit[tower][il][ixy] -> GetParameters(parMulti);
161  mEachFit[tower][il][ixy][0] -> SetParameters(parMulti[0], parMulti[1], parMulti[2], parMulti[3], parMulti[4], parMulti[5]/2.);
162  mEachFit[tower][il][ixy][1] -> SetParameters(parMulti[6], parMulti[7], parMulti[8], parMulti[9], parMulti[10], parMulti[5]/2.);
163  }
164  }
165  return kRHICfOk;
166 }
167 
168 int StRHICfRecoPos::getGSOMaxLayer(int tower, int order){return mGSOMaxLayer[tower][order];}
169 
170 int StRHICfRecoPos::getMaximumBin(int tower, int layer, int xy)
171 {
172  int bin = 0;
173  double tmpVal = 0;
174  double tmpE = 0;
175 
176  for(int ich=0; ich<checkGSOBarSize(tower); ich++){
177  if(tower==0){tmpE = mGSOBarSE[layer][xy][ich];}
178  if(tower==1){tmpE = mGSOBarTE[layer][xy][ich];}
179  if(tmpE > tmpVal){
180  tmpVal = tmpE;
181  bin = ich;
182  }
183  }
184  return bin+1;
185 }
186 
187 int StRHICfRecoPos::getEvalHitNum(int tower)
188 {
189  int nHit[kRHICfNlayer][kRHICfNxy];
190  for(int il=0; il<kRHICfNlayer; il++){
191  for(int ixy=0; ixy<kRHICfNxy; ixy++){
192  nHit[il][ixy] = getEvalHitNum(tower, il, ixy);
193  }
194  }
195 
196  bool isNullData = true;
197  for(int il=0; il<4; ++il){
198  if(mGSOMaxLayer[tower][0]==il && nHit[il][0] > 0 && nHit[il][1] > 0){isNullData = false;}
199  }
200 
201  mOverlap[tower][0] = false;
202  mOverlap[tower][1] = false;
203 
204  if(isNullData == true){return 0;}
205 
206  if(nHit[mGSOMaxLayer[tower][0]][0]==2 && nHit[mGSOMaxLayer[tower][0]][1]==2){return 2;}
207 
208  if(mMultiFit[tower][mGSOMaxLayer[tower][0]][0]->GetParameter(1) == mMultiFit[tower][mGSOMaxLayer[tower][0]][0]->GetParameter(7) && nHit[mGSOMaxLayer[tower][0]][1] == 2){
209  mOverlap[tower][0] = true;
210  return 2;
211  }
212  if(mMultiFit[tower][mGSOMaxLayer[tower][0]][1]->GetParameter(1) == mMultiFit[tower][mGSOMaxLayer[tower][0]][1]->GetParameter(7) && nHit[mGSOMaxLayer[tower][0]][0] == 2){
213  mOverlap[tower][1] = true;
214  return 2;
215  }
216  return 1;
217 }
218 
219 int StRHICfRecoPos::getEvalHitNum(int tower, int layer, int xy)
220 {
221  int nHitSingleFit = 0;
222  int nHitMultiFit = 0;
223 
224  // Check nHit to fitting result for single
225  if(getSinglePeakHeight(tower, layer, xy) > mNoiseThreshold){nHitSingleFit = 1;}
226 
227  // Check nHit to fitting result for multi
228  if(getMultiPeakHeight(tower, layer, xy, 0) > mNoiseThreshold){nHitMultiFit++;}
229  if(getMultiPeakHeight(tower, layer, xy, 1) > mNoiseThreshold){nHitMultiFit++;}
230 
231  // Check nHit to TSpectrum
232  if(mSpectrum[tower][layer][xy]->GetNPeaks() < 2){
233  nHitMultiFit = 1;
234  mMultiFit[tower][layer][xy] -> SetParameter(7, mMultiFit[tower][layer][xy]->GetParameter(1));
235  }
236 
237  // Check nHit to distance
238  Double_t parMulti[11];
239  mMultiFit[tower][layer][xy] -> GetParameters(parMulti);
240  if(abs(parMulti[1] - parMulti[7]) < mPeakDistThreshold){
241  nHitMultiFit = 1;
242  mMultiFit[tower][layer][xy] -> SetParameter(7, mMultiFit[tower][layer][xy]->GetParameter(1));
243  }
244 
245  // Check nHit to peak height ratio.
246  double heightRatio = getMultiPeakHeight(tower, layer, xy, 1)/getMultiPeakHeight(tower, layer, xy, 0);
247  if(heightRatio < mTSpecRatioThreshold){nHitMultiFit = 1;}
248 
249  double chi2Ratio = getSingleChi2(tower, layer, xy)/getMultiChi2(tower, layer, xy);
250 
251  if(nHitMultiFit > 1 && chi2Ratio > 1.0){return 2;}
252  if(nHitSingleFit > 0){return 1;}
253  return 0;
254 }
255 
256 double StRHICfRecoPos::getDepositEnergy(int tower, int layer)
257 {
258  double edepX = mGSOBarHist[tower][layer][0] -> Integral();
259  double edepY = mGSOBarHist[tower][layer][1] -> Integral();
260  return edepX+edepY;
261 }
262 
263 double StRHICfRecoPos::getMultiEnergySum(int tower, int layer, int xy, int order){return mEachFit[tower][layer][xy][order] -> Integral(0, checkGSOBarSize(tower));}
264 
265 Double_t StRHICfRecoPos::getSinglePeakPos(int tower, int layer, int xy){return mSingleFit[tower][layer][xy] -> GetParameter(1);}
266 
267 Double_t StRHICfRecoPos::getMultiPeakPos(int tower, int layer, int xy, int order)
268 {
269  if(order==0){return mMultiFit[tower][layer][xy] -> GetParameter(1);}
270  else{return mMultiFit[tower][layer][xy] -> GetParameter(7);}
271 }
272 
273 double StRHICfRecoPos::getSinglePeakHeight(int tower, int layer, int xy){return mSingleFit[tower][layer][xy] ->Eval(getSinglePeakPos(tower, layer, xy));}
274 
275 double StRHICfRecoPos::getMultiPeakHeight(int tower, int layer, int xy, int order)
276 {
277  double peak[2];
278  Double_t multiPar[11];
279  memset(peak, 0, sizeof(peak));
280  memset(multiPar, 0, sizeof(multiPar));
281 
282  mMultiFit[tower][layer][xy] -> GetParameters(multiPar);
283 
284  peak[0] = multiPar[2]*(multiPar[4]/pow(multiPar[0], 0.5) + (1 - multiPar[4])/pow(multiPar[3], 0.5))/2. + multiPar[5];
285  peak[1] = multiPar[8]*(multiPar[10]/pow(multiPar[6], 0.5) + (1 - multiPar[10])/pow(multiPar[9], 0.5))/2. + multiPar[5];
286 
287  if(order==0){return std::max(peak[0], peak[1]);}
288  else{return std::min(peak[0], peak[1]);}
289 }
290 
291 double StRHICfRecoPos::getMultiPeakRaw(int tower, int layer, int xy, int order)
292 {
293  Int_t posBin = mGSOBarHist[tower][layer][xy] -> FindBin(getMultiPeakPos(tower, layer, xy, order));
294  double candidateBins[3] = {mGSOBarHist[tower][layer][xy]->GetBinContent(posBin-1), mGSOBarHist[tower][layer][xy]->GetBinContent(posBin), mGSOBarHist[tower][layer][xy]->GetBinContent(posBin+1)};
295  int maxIdx = TMath::LocMax(3, candidateBins);
296  return candidateBins[maxIdx];
297 }
298 
299 double StRHICfRecoPos::getEvalPeakHeight(int tower, int layer, int xy, int order){return mEachFit[tower][layer][xy][order] -> Eval(getMultiPeakPos(tower, layer, xy, order));}
300 
301 double StRHICfRecoPos::getSingleChi2(int tower, int layer, int xy)
302 {
303  double chi2 = mSingleFit[tower][layer][xy]->GetChisquare();
304  double dof = (double) mSingleFit[tower][layer][xy]->GetNDF();
305  return chi2/dof;
306 }
307 
308 double StRHICfRecoPos::getMultiChi2(int tower, int layer, int xy)
309 {
310  double chi2 = mMultiFit[tower][layer][xy]->GetChisquare();
311  double dof = (double) mMultiFit[tower][layer][xy]->GetNDF();
312  return chi2/dof;
313 }
314 
315 bool StRHICfRecoPos::getWorthy(int tower){return mWorthy[tower];}
316 
317 bool StRHICfRecoPos::getOverlap(int tower, int xy){return mOverlap[tower][xy];}
318 
319 void StRHICfRecoPos::findMaxLayer()
320 {
321  double maxSumBar[kRHICfNtower][kRHICfNlayer];
322  for(int it=0; it<kRHICfNtower; it++){
323  for(int il=0; il<kRHICfNlayer; il++){maxSumBar[it][il] = getDepositEnergy(it, il);}
324  mGSOMaxLayer[it][0] = TMath::LocMax(4, maxSumBar[it]);
325  maxSumBar[it][mGSOMaxLayer[it][0]] = 0.;
326  mGSOMaxLayer[it][1] = TMath::LocMax(4, maxSumBar[it]);
327  }
328 }
329 
330 void StRHICfRecoPos::initSetParamaters()
331 {
332  Double_t parInMaxLayer[2][kRHICfNtower][kRHICfNxy][2];
333  Double_t parSingle[kRHICfNtower][kRHICfNlayer][kRHICfNxy][2];
334 
335  memset(parInMaxLayer, 0, sizeof(parInMaxLayer));
336  memset(parSingle, 0, sizeof(parSingle));
337 
338  for(int it=0; it<kRHICfNtower; it++){
339  for(int il=0; il<kRHICfNlayer; il++){
340  for(int ixy=0; ixy<kRHICfNxy; ixy++){
341  Double_t posPrime = 0.;
342  Double_t valuePrime = 0.;
343  Double_t posSub = 0.;
344  Double_t valueSub = 0.;
345 
346  // Search the peak
347  mSpectrum[it][il][ixy]->Search(mGSOBarHistExpend[it][il][ixy], mTSpecSigma, "nodraw goff", mTSpecRatioThreshold);
348  int numHit = mSpectrum[it][il][ixy]->GetNPeaks();
349  auto tSpecPeak = mSpectrum[it][il][ixy]->GetPositionX();
350 
351  // Setting prime, sub position and hight
352  if(numHit==0){
353  Int_t maxBin = mGSOBarHist[it][il][ixy]->GetMaximumBin();
354  Int_t subBin = checkGSOBarSize(it)/2;
355 
356  posPrime = mGSOBarHist[it][il][ixy]->GetBinCenter(maxBin);
357  valuePrime = mGSOBarHist[it][il][ixy]->GetBinContent(maxBin);
358  posSub = Double_t(subBin);
359  valueSub = mGSOBarHist[it][il][ixy]->GetBinContent(subBin);
360  }
361  else if(numHit==1){
362  Int_t maxBin = mGSOBarHist[it][il][ixy]->GetXaxis()->FindBin(tSpecPeak[0]);
363  Int_t subBin = checkGSOBarSize(it)/2;
364 
365  posPrime = tSpecPeak[0];
366  valuePrime = mGSOBarHist[it][il][ixy]->GetBinContent(maxBin);
367  posSub = Double_t(subBin);
368  valueSub = mGSOBarHist[it][il][ixy]->GetBinContent(subBin);
369  }
370  else if(numHit>1){
371  double bufferPeakValue = 0;
372  int primeIdx = -1;
373 
374  // searching for prime peak
375  for(int in=0; in<numHit; in++){
376  Int_t bin = mGSOBarHist[it][il][ixy]->GetXaxis()->FindBin(tSpecPeak[in]);
377  Double_t value = mGSOBarHist[it][il][ixy]->GetBinContent(bin);
378 
379  if(value > bufferPeakValue){
380  bufferPeakValue = value;
381  posPrime = tSpecPeak[in];
382  valuePrime = value;
383  primeIdx = in;
384  }
385  }
386 
387  bufferPeakValue = 0;
388  //searching for sub peak
389  for(int in=0; in<numHit; in++){
390  if(in == primeIdx){continue;}
391 
392  Int_t bin = mGSOBarHist[it][il][ixy]->GetXaxis()->FindBin(tSpecPeak[in]);
393  Double_t value = mGSOBarHist[it][il][ixy]->GetBinContent(bin);
394 
395  if(value > bufferPeakValue){
396  bufferPeakValue = value;
397  posSub = tSpecPeak[in];
398  valueSub = value;
399  }
400  }
401  }
402 
403  if(valuePrime<0){valuePrime = -valuePrime;}
404  if(valueSub<0){valueSub = -valueSub;}
405 
406  if(getGSOMaxLayer(it, 0) == il){
407  parInMaxLayer[0][it][ixy][0] = posPrime;
408  parInMaxLayer[0][it][ixy][1] = valuePrime;
409  parInMaxLayer[1][it][ixy][0] = posSub;
410  parInMaxLayer[1][it][ixy][1] = valueSub;
411  }
412  parSingle[it][il][ixy][0] = posPrime;
413  parSingle[it][il][ixy][1] = valuePrime;
414  }
415  }
416  }
417 
418  // Parameters setup
419  for(int it=0; it<kRHICfNtower; it++){
420  for(int il=0; il<kRHICfNlayer; il++){
421  for(int ixy=0; ixy<kRHICfNxy; ixy++){
422  Double_t initParSingle[6] = {mParWidth1, parSingle[it][il][ixy][0], parSingle[it][il][ixy][1], mParWidth2, mParRatio, mParBaseLine};
423  mSingleFit[it][il][ixy] -> SetParameters(initParSingle);
424 
425  mSingleFit[it][il][ixy] -> SetParLimits(0, mParWidth1Min, mParWidth1Max);
426  mSingleFit[it][il][ixy] -> SetParLimits(1, -1., 1. + double(checkGSOBarSize(it)));
427  mSingleFit[it][il][ixy] -> SetParLimits(2, mParHeightMin, mParHeightMax);
428  mSingleFit[it][il][ixy] -> SetParLimits(3, mParWidth2Min, mParWidth2Max);
429  mSingleFit[it][il][ixy] -> SetParLimits(4, mParRatioMin, mParRatioMax);
430 
431  Double_t initParMulti[11] = {mParWidth1, parInMaxLayer[0][it][ixy][0], parInMaxLayer[0][it][ixy][1], mParWidth2, mParRatio, mParBaseLine, mParWidth1, parInMaxLayer[1][it][ixy][0], parInMaxLayer[1][it][ixy][1], mParWidth2, mParRatio};
432  mMultiFit[it][il][ixy] -> SetParameters(initParMulti);
433 
434  mMultiFit[it][il][ixy] -> SetParLimits(0, mParWidth1Min, mParWidth1Max);
435  mMultiFit[it][il][ixy] -> SetParLimits(1, parInMaxLayer[0][it][ixy][0]-2, parInMaxLayer[0][it][ixy][0]+2);
436  mMultiFit[it][il][ixy] -> SetParLimits(2, mParHeightMin, mParHeightMax);
437  mMultiFit[it][il][ixy] -> SetParLimits(3, mParWidth2Min, mParWidth2Max);
438  mMultiFit[it][il][ixy] -> SetParLimits(4, mParRatioMin, mParRatioMax);
439  mMultiFit[it][il][ixy] -> SetParLimits(6, mParWidth1Min, mParWidth1Max);
440  mMultiFit[it][il][ixy] -> SetParLimits(7, parInMaxLayer[1][it][ixy][0]-2, parInMaxLayer[1][it][ixy][0]+2);
441  mMultiFit[it][il][ixy] -> SetParLimits(8, mParHeightMin, mParHeightMax);
442  mMultiFit[it][il][ixy] -> SetParLimits(9, mParWidth2Min, mParWidth2Max);
443  mMultiFit[it][il][ixy] -> SetParLimits(10, mParRatioMin, mParRatioMax);
444  }
445  }
446  }
447 }
448 
449 void StRHICfRecoPos::fitting()
450 {
451  for(int it=0; it<kRHICfNtower; it++){
452  for(int il=0; il<kRHICfNlayer; il++){
453  for(int ixy=0; ixy<kRHICfNxy; ixy++){
454  mGSOBarGraph[it][il][ixy] -> Fit(mSingleFit[it][il][ixy], "QR");
455  mGSOBarGraph[it][il][ixy] -> Fit(mMultiFit[it][il][ixy], "QR");
456  }
457  }
458  }
459 }
460 
461 Double_t StRHICfRecoPos::getLorenzianSingle(Double_t* x, Double_t* par)
462 {
463  Double_t val = par[4]*par[0]/pow(((x[0]-par[1])*(x[0]-par[1])+par[0]), 1.5)/2.;
464  Double_t val2 = (1-par[4])*(par[3]/2.)/pow(((x[0]-par[1])*(x[0]-par[1])+par[3]), 1.5);
465 
466  return par[2]*(val+val2) + par[5];
467 }
468 
469 Double_t StRHICfRecoPos::getLorenzianMulti(Double_t* x, Double_t* par)
470 {
471  Double_t p1val = par[4]*par[0]/pow(((x[0]-par[1])*(x[0]-par[1])+par[0]), 1.5)/2.;
472  Double_t p1val2 = (1.-par[4])*par[3]/pow(((x[0]-par[1])*(x[0]-par[1])+par[3]), 1.5)/2.;
473  Double_t peak1 = par[2]*(p1val+p1val2);
474 
475  Double_t p2val = par[10]*par[6]/pow(((x[0]-par[7])*(x[0]-par[7])+par[6]), 1.5)/2.;
476  Double_t p2val2 = (1.-par[10])*par[9]/pow(((x[0]-par[7])*(x[0]-par[7])+par[9]), 1.5)/2.;
477  Double_t peak2 = par[8]*(p2val+p2val2);
478 
479  return peak1 + peak2 + par[5];
480 }
Generated by   doxygen 1.8.5