StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StHbtCorrFctnDirectYlm.cxx
1 /****************************************************************************
2  * $Id: StHbtCorrFctnDirectYlm.cxx,v 1.2 2015/11/02 20:11:06 perev Exp $
3  *
4  * Author: Yan Yang, hce137@gmail.com
5  * **************************************************************************
6  * Description: Correlation function that is binned in Ylms directly
7  * Provides a way to store the numerator and denominator
8  * in Ylms directly and correctly calculate the correlation
9  * function from them.
10  * The original author is not known, I guess Kisiel would be.
11  * *************************************************************************/
12 
13 #include "StHbtCorrFctnDirectYlm.h"
14 
15 #include <TMath.h>
16 #include "gsl/gsl_blas.h"
17 #include "gsl/gsl_linalg.h"
18 #include "gsl/gsl_complex.h"
19 #include <iostream>
20 #include <cmath>
21 #include "TMatrixD.h"
22 
23 using namespace std;
24 
25 //#define YY_DEBUG
26 #define MAXJM 49
27 
28 #ifdef __ROOT__
29 ClassImp(StHbtCorrFctnDirectYlm)
30 #endif
31 
32 StHbtCorrFctnDirectYlm::StHbtCorrFctnDirectYlm(const char* name, int maxl, int ibin = 30, double vmin = 0.0, double vmax = 0.3)
33  : cftnreal(0), cftnimag(0), covnum(0), covden(0), fR2factor(0),
34  mDoEMCIC(false), emcicP1P2T(NULL), emcicP1P2Z(NULL), emcicE1E2(NULL), emcicE1plusE2(NULL),
35  mName(name), mNbins(ibin), mkmin(vmin), mkmax(vmax),
36  mNormRadius(0.0), mNormBohr(0.0), mNormBinMin(0), mNormBinMax(0)
37 {
38  fMaxL = maxl;
39  maxjm = (maxl + 1) * (maxl + 1);
40 
41  //Fill in factorials table
42  factorials = new double[( 4 * (maxl + 1) )];
43  int fac = 1;
44  factorials[0] = 1;
45  for (int iter = 1; iter < 4 * (maxl + 1); iter++) {
46  fac *= iter;
47  factorials[iter] = fac;
48  }
49 
50  //Fill in els and ems table
51  int el = 0;
52  int em = 0;
53  int il = 0;
54  els = new double[maxjm];
55  ems = new double[maxjm];
56  elsi = new int[maxjm];
57  emsi = new int[maxjm];
58  do {
59  els[il] = el;
60  ems[il] = em;
61  elsi[il] = (int) el;
62  emsi[il] = (int) em;
63 
64  //cout << "il el em " << il << " " << elsi[il] << " " << emsi[il] << endl;
65  em++;
66  il++;
67  if (em > el) {
68  el++;
69  em = -el;
70  }
71  } while (el <= maxl);
72 
73  //Create numerator and denominator histograms
74  numsreal = new TH1D*[maxjm];
75  numsimag = new TH1D*[maxjm];
76  densreal = new TH1D*[maxjm];
77  densimag = new TH1D*[maxjm];
78 
79  char bufname[200];
80  for (int ihist = 0; ihist < maxjm; ihist++) {
81  int em = emsi[ihist] < 0 ? elsi[ihist] - emsi[ihist] : emsi[ihist];
82 
83  sprintf(bufname, "NumReYlm%i%i%s", elsi[ihist], em, name);
84  numsreal[ihist] = new TH1D(bufname, bufname, ibin, vmin, vmax);
85  sprintf(bufname, "NumImYlm%i%i%s", elsi[ihist], em, name);
86  numsimag[ihist] = new TH1D(bufname, bufname, ibin, vmin, vmax);
87  sprintf(bufname, "DenReYlm%i%i%s", elsi[ihist], em, name);
88  densreal[ihist] = new TH1D(bufname, bufname, ibin, vmin, vmax);
89  sprintf(bufname, "DenImYlm%i%i%s", elsi[ihist], em, name);
90  densimag[ihist] = new TH1D(bufname, bufname, ibin, vmin, vmax);
91 
92  numsreal[ihist]->Sumw2();
93  numsimag[ihist]->Sumw2();
94  densreal[ihist]->Sumw2();
95  densimag[ihist]->Sumw2();
96  }
97 
98  sprintf(bufname, "BinCountNum%s", name);
99  binctn = new TH1D(bufname, bufname, ibin, vmin, vmax);
100 
101  sprintf(bufname, "BinCountDen%s", name);
102  binctd = new TH1D(bufname, bufname, ibin, vmin, vmax);
103 
104  ylmbuffer = new complex<double>[maxjm];
105 
106  //Covariance matrices
107  covmnum = new double[maxjm * maxjm * 4 * ibin];
108  covmden = new double[maxjm * maxjm * 4 * ibin];
109 }
110 
111 //_________________________________________________________________________
112 void StHbtCorrFctnDirectYlm::DoEMCIC()
113 {
114  mDoEMCIC = true;
115  emcicP1P2T = new TH1D*[maxjm];
116  emcicP1P2Z = new TH1D*[maxjm];
117  emcicE1E2 = new TH1D*[maxjm];
118  emcicE1plusE2 = new TH1D*[maxjm];
119  char bufname[200];
120  for (int ihist = 0; ihist < maxjm; ihist++) {
121  int em = emsi[ihist] < 0 ? elsi[ihist] - emsi[ihist] : emsi[ihist];
122  sprintf( bufname, "emcicP1P2T%i%i%s", elsi[ihist], em, mName.c_str() );
123  emcicP1P2T[ihist] = new TH1D(bufname, bufname, mNbins, mkmin, mkmax);
124  emcicP1P2T[ihist]->Sumw2();
125  sprintf( bufname, "emcicP1P2Z%i%i%s", elsi[ihist], em, mName.c_str() );
126  emcicP1P2Z[ihist] = new TH1D(bufname, bufname, mNbins, mkmin, mkmax);
127  emcicP1P2Z[ihist]->Sumw2();
128  sprintf( bufname, "emcicE1E2%i%i%s", elsi[ihist], em, mName.c_str() );
129  emcicE1E2[ihist] = new TH1D(bufname, bufname, mNbins, mkmin, mkmax);
130  emcicE1E2[ihist]->Sumw2();
131  sprintf( bufname, "emcicE1plusE2%i%i%s", elsi[ihist], em, mName.c_str() );
132  emcicE1plusE2[ihist] = new TH1D(bufname, bufname, mNbins, mkmin, mkmax);
133  emcicE1plusE2[ihist]->Sumw2();
134  }
135 }
136 
137 //__________________________________________________________________________
138 StHbtCorrFctnDirectYlm::StHbtCorrFctnDirectYlm()
139 {
140  StHbtCorrFctnDirectYlm("StHbtCorrFctnDirectYlm", 2);
141 }
142 
143 //___________________________________________________________________________
144 StHbtCorrFctnDirectYlm::~StHbtCorrFctnDirectYlm()
145 {
146 //cout << "StHbtCorrFctnDirectYlm::~StHbtCorrFctnDirectYlm()\n";
147  for (int ihist = 0; ihist < maxjm; ihist++) {
148  if (numsreal[ihist]) delete numsreal[ihist];
149  if (numsimag[ihist]) delete numsimag[ihist];
150  if (densreal[ihist]) delete densreal[ihist];
151  if (densimag[ihist]) delete densimag[ihist];
152  if (cftnreal && cftnreal[ihist]) delete cftnreal[ihist];
153  if (cftnimag && cftnimag[ihist]) delete cftnimag[ihist];
154  if (mDoEMCIC) {
155  delete emcicP1P2T[ihist];
156  delete emcicP1P2Z[ihist];
157  delete emcicE1E2[ihist];
158  delete emcicE1plusE2[ihist];
159  }
160  }
161 
162  if (mDoEMCIC) {
163  delete []emcicP1P2T;
164  delete []emcicP1P2Z;
165  delete []emcicE1E2;
166  delete []emcicE1plusE2;
167  }
168 
169  delete binctn;
170  delete binctd;
171 
172  delete []numsreal;
173  delete []numsimag;
174  delete []densreal;
175  delete []densimag;
176  if (cftnreal) delete []cftnreal;
177  if (cftnimag) delete []cftnimag;
178 
179  delete []factorials;
180  delete []els;
181  delete []ems;
182  delete []elsi;
183  delete []emsi;
184  delete []ylmbuffer;
185 
186  delete []covmnum;
187  delete []covmden;
188 
189  delete covnum;
190  delete covden;
191 
192 }
193 
194 //____________________________________________________________________________
195 double StHbtCorrFctnDirectYlm::ClebschGordan(double aJot1, double aEm1, double aJot2, double aEm2, double aJot, double aEm)
196 {
197  int mint, maxt;
198  double cgc = 0.0;
199  int titer;
200  double coef;
201 
202  maxt = lrint(aJot1 + aJot2 - aJot);
203  mint = 0;
204  if (lrint(aJot1 - aEm1) < maxt) maxt = lrint(aJot1 - aEm1);
205  if (lrint(aJot2 + aEm2) < maxt) maxt = lrint(aJot2 + aEm2);
206  if (lrint( -(aJot - aJot2 + aEm1) ) > mint) mint = lrint( -(aJot - aJot2 + aEm1) );
207  if (lrint( -(aJot - aJot1 - aEm2) ) > mint) mint = lrint( -(aJot - aJot1 - aEm2) );
208 
209  for (titer = mint; titer <= maxt; titer++) {
210  coef = TMath::Power(-1, titer);
211  coef *= TMath::Sqrt( (2 * aJot + 1) *
212  factorials[lrint(aJot1 + aEm1)] *
213  factorials[lrint(aJot1 - aEm1)] *
214  factorials[lrint(aJot2 + aEm2)] *
215  factorials[lrint(aJot2 - aEm2)] *
216  factorials[lrint(aJot + aEm)] *
217  factorials[lrint(aJot - aEm)] );
218  coef /= (factorials[titer] *
219  factorials[lrint(aJot1 + aJot2 - aJot - titer)] *
220  factorials[lrint(aJot1 - aEm1 - titer)] *
221  factorials[lrint(aJot2 + aEm2 - titer)] *
222  factorials[lrint(aJot - aJot2 + aEm1 + titer)] *
223  factorials[lrint(aJot - aJot1 - aEm2 + titer)]);
224 
225  cgc += coef;
226  }
227 
228  cgc *= DeltaJ(aJot1, aJot2, aJot);
229 
230  return (cgc);
231 }
232 
233 //____________________________________________________________________________
234 double StHbtCorrFctnDirectYlm::DeltaJ(double aJot1, double aJot2, double aJot)
235 {
236  if ( (aJot1 + aJot2 - aJot) < 0 ) {
237  //cout << "J1+J2-J3 < 0 !!!" << " " << aJot1 << " " << aJot2 << " " << aJot << endl;
238  return (0);
239  }
240  if ( (aJot1 - aJot2 + aJot) < 0 ) {
241  //cout << "J1-J2+J3 < 0 !!!" << " " << aJot1 << " " << aJot2 << " " << aJot << endl;
242  return (0);
243  }
244  if ( (-aJot1 + aJot2 + aJot) < 0 ) {
245  //cout << "-J1+J2+J3 < 0 !!!" << " " << aJot1 << " " << aJot2 << " " << aJot << endl;
246  return (0);
247  }
248  if ( (aJot1 + aJot2 + aJot + 1) < 0 ) {
249  //cout << "J1+J2+J3+1 < 0 !!!" << " " << aJot1 << " " << aJot2 << " " << aJot << endl;
250  return (0);
251  }
252  double res = TMath::Sqrt(1.0 *
253  factorials[lrint(aJot1 + aJot2 - aJot)] *
254  factorials[lrint(aJot1 - aJot2 + aJot)] *
255  factorials[lrint(-aJot1 + aJot2 + aJot)] /
256  factorials[lrint(aJot1 + aJot2 + aJot + 1)]);
257 
258  return (res);
259 }
260 
261 //___________________________________________________________________________
262 double StHbtCorrFctnDirectYlm::WignerSymbol(double aJot1, double aEm1, double aJot2, double aEm2, double aJot, double aEm)
263 {
264  if (lrint(aEm1 + aEm2 + aEm) != 0.0)
265  return (0.0);
266  double cge = ClebschGordan(aJot1, aEm1, aJot2, aEm2, aJot, -aEm);
267  if (lrint( abs(aJot1 - aJot2 - aEm) ) % 2)
268  cge *= -1.0;
269  cge /= sqrt(2 * aJot + 1);
270 
271  if (cge == -0.0) cge = 0.0;
272 
273  return (cge);
274 }
275 
276 //____________________________________________________________________________
277 void StHbtCorrFctnDirectYlm::GetMtilde(complex<double>* aMat, double* aMTilde)
278 {
279  //Create the Mtilde for a given q bin
280  double lzero, mzero;
281  double lprim, mprim;
282  double lbis, mbis;
283 
284  int lzeroi, mzeroi;
285  int lprimi, mprimi;
286  int lbisi, mbisi;
287 
288  for (int iz = 0; iz < GetMaxJM() * 2; iz++)
289  for (int ip = 0; ip < GetMaxJM() * 2; ip++)
290  aMTilde[iz * GetMaxJM() * 2 + ip] = 0.0;
291 
292  for (int izero = 0; izero < GetMaxJM(); izero++) {
293  GetElEmForIndex(izero, &lzero, &mzero);
294  GetElEmForIndex(izero, &lzeroi, &mzeroi);
295  for (int ibis = 0; ibis < GetMaxJM(); ibis++) {
296  GetElEmForIndex(ibis, &lbis, &mbis);
297  GetElEmForIndex(ibis, &lbisi, &mbisi);
298  complex<double> val = complex<double>(0.0, 0.0);
299  complex<double> mcomp[MAXJM];
300  for (int iprim = 0; iprim < GetMaxJM(); iprim++) {
301  GetElEmForIndex(iprim, &lprim, &mprim);
302  GetElEmForIndex(iprim, &lprimi, &mprimi);
303  //(-1)^m
304  if (abs(mzeroi) % 2) mcomp[iprim] = complex<double>(-1.0, 0.0);
305  else mcomp[iprim] = complex<double>(1.0, 0.0);
306 
307  //P1
308  mcomp[iprim] *= sqrt( (2 * lzero + 1) * (2 * lprim + 1) * (2 * lbis + 1) );
309  //W1
310  mcomp[iprim] *= WignerSymbol(lzero, 0, lprim, 0, lbis, 0);
311  //W2
312  mcomp[iprim] *= WignerSymbol(lzero, -mzero, lprim, mprim, lbis, mbis);
313  mcomp[iprim] *= aMat[iprim];
314  val += mcomp[iprim];
315  }
316 
317  aMTilde[(izero * 2) * ( 2 * GetMaxJM() ) + (ibis * 2)] = real(val);
318  aMTilde[(izero * 2 + 1) * ( 2 * GetMaxJM() ) + (ibis * 2)] = imag(val);
319  if (imag(val) != 0.0)
320  aMTilde[(izero * 2) * ( 2 * GetMaxJM() ) + (ibis * 2 + 1)] = -imag(val);
321  else
322  aMTilde[(izero * 2) * ( 2 * GetMaxJM() ) + (ibis * 2 + 1)] = 0.0;
323  aMTilde[(izero * 2 + 1) * ( 2 * GetMaxJM() ) + (ibis * 2 + 1)] = real(val);
324  }
325  }
326 }
327 
328 int StHbtCorrFctnDirectYlm::GetMaxJM()
329 { return (maxjm); }
330 
331 void StHbtCorrFctnDirectYlm::GetElEmForIndex(int aIndex, double* aEl, double* aEm)
332 {
333  *aEl = els[aIndex];
334  *aEm = ems[aIndex];
335 }
336 
337 void StHbtCorrFctnDirectYlm::GetElEmForIndex(int aIndex, int* aEl, int* aEm)
338 {
339  *aEl = elsi[aIndex];
340  *aEm = emsi[aIndex];
341 }
342 
343 int StHbtCorrFctnDirectYlm::GetBin(int qbin, int ilmzero, int zeroimag, int ilmprim, int primimag)
344 {
345  return (qbin * GetMaxJM() * GetMaxJM() * 4 +
346  (ilmprim * 2 + primimag) * GetMaxJM() * 2 +
347  ilmzero * 2 + zeroimag);
348 }
349 
350 //_____________________________________________________________________________
351 void StHbtCorrFctnDirectYlm::AddRealPair(const StHbtPair* aPair, double weight)
352 {
353  double qout = aPair->dKOut();
354  double qside = aPair->dKSide();
355  double qlong = aPair->dKLong();
356  double kv = sqrt(qout * qout + qside * qside + qlong * qlong);
357  if (kv >= mkmax || kv < mkmin) return; //yyang
358  int nqbin = binctn->GetXaxis()->FindFixBin(kv) - 1;
359  if (fR2factor) weight *= ( 1.0 / (kv * kv) );
360 
361  StHbtYlm::YlmUpToL(elsi[GetMaxJM() - 1], qout, qside, qlong, ylmbuffer);
362  //check is it not a number
363  for (int ilm = 0; ilm < GetMaxJM(); ilm++) {
364  if ( TMath::IsNaN(real(ylmbuffer[ilm]) * weight) ) return;
365  if ( TMath::IsNaN(imag(ylmbuffer[ilm]) * weight) ) return;
366  if ( ::isinf(real(ylmbuffer[ilm]) * weight) ) return;
367  if ( ::isinf(imag(ylmbuffer[ilm]) * weight) ) return;
368  }
369  if ( nqbin < binctn->GetNbinsX() )
370  for (int ilmzero = 0; ilmzero < GetMaxJM(); ilmzero++)
371  for (int ilmprim = 0; ilmprim < GetMaxJM(); ilmprim++) {
372  if ( TMath::IsNaN(real(ylmbuffer[ilmzero]) * real(ylmbuffer[ilmprim]) * weight * weight) ) return;
373  if ( TMath::IsNaN(real(ylmbuffer[ilmzero]) * -imag(ylmbuffer[ilmprim]) * weight * weight) ) return;
374  if ( TMath::IsNaN(imag(ylmbuffer[ilmzero]) * real(ylmbuffer[ilmprim]) * weight * weight) ) return;
375  if ( TMath::IsNaN(imag(ylmbuffer[ilmzero]) * -imag(ylmbuffer[ilmprim]) * weight * weight) ) return;
376  if ( ::isinf(real(ylmbuffer[ilmzero]) * real(ylmbuffer[ilmprim]) * weight * weight) ) return;
377  if ( ::isinf(real(ylmbuffer[ilmzero]) * -imag(ylmbuffer[ilmprim]) * weight * weight) ) return;
378  if ( ::isinf(imag(ylmbuffer[ilmzero]) * real(ylmbuffer[ilmprim]) * weight * weight) ) return;
379  if ( ::isinf(imag(ylmbuffer[ilmzero]) * -imag(ylmbuffer[ilmprim]) * weight * weight) ) return;
380  }
381 
382  for (int ilm = 0; ilm < GetMaxJM(); ilm++) {
383  numsreal[ilm]->Fill(kv, real(ylmbuffer[ilm]) * weight);
384  numsimag[ilm]->Fill(kv, -imag(ylmbuffer[ilm]) * weight);
385  binctn->Fill(kv, 1.0);
386  }
387  if (mDoEMCIC) {
388  double p1p2t = sqrt( aPair->track1()->FourMomentum().px() * aPair->track1()->FourMomentum().py() + aPair->track2()->FourMomentum().px() * aPair->track2(
389  )->FourMomentum().py() );
390  double p1p2z = aPair->track1()->FourMomentum().pz() * aPair->track2()->FourMomentum().pz();
391  double e1e2 = aPair->track1()->FourMomentum().e() * aPair->track2()->FourMomentum().e();
392  double e1pluse2 = aPair->track1()->FourMomentum().e() + aPair->track2()->FourMomentum().e();
393  if ( !( ::isnan(p1p2t) || ::isnan(p1p2z) || ::isnan(e1e2) || ::isnan(e1pluse2) ||
394  ::isinf(p1p2t) || ::isinf(p1p2z) || ::isinf(e1e2) || ::isinf(e1pluse2) ) )
395  {
396  for (int ilm = 0; ilm < GetMaxJM(); ilm++) {
397  emcicP1P2T[ilm]->Fill(kv, p1p2t);
398  emcicP1P2Z[ilm]->Fill(kv, p1p2z);
399  emcicE1plusE2[ilm]->Fill(kv, e1pluse2);
400  emcicE1E2[ilm]->Fill(kv, e1e2);
401  }
402  }
403  } //if mDoEMCIC
404 
405  //Fill in the error matrix
406  //int tabshift = nqbin*GetMaxJM()*GetMaxJM()*4;
407  //the following code needs to be revised!! // yyang
408  if ( nqbin < binctn->GetNbinsX() )
409  for (int ilmzero = 0; ilmzero < GetMaxJM(); ilmzero++)
410  for (int ilmprim = 0; ilmprim < GetMaxJM(); ilmprim++) {
411  covmnum[GetBin(nqbin, ilmzero, 0, ilmprim, 0)] += real(ylmbuffer[ilmzero]) * real(ylmbuffer[ilmprim]) * weight * weight;
412  covmnum[GetBin(nqbin, ilmzero, 0, ilmprim, 1)] += real(ylmbuffer[ilmzero]) * -imag(ylmbuffer[ilmprim]) * weight * weight;
413  covmnum[GetBin(nqbin, ilmzero, 1, ilmprim, 0)] += -imag(ylmbuffer[ilmzero]) * real(ylmbuffer[ilmprim]) * weight * weight;
414  covmnum[GetBin(nqbin, ilmzero, 1, ilmprim, 1)] += -imag(ylmbuffer[ilmzero]) * -imag(ylmbuffer[ilmprim]) * weight * weight;
415  }
416 }
417 
418 //____________________________________________________________________________
419 void StHbtCorrFctnDirectYlm::AddMixedPair(const StHbtPair* aPair, double weight)
420 {
421  double qout = aPair->dKOut();
422  double qside = aPair->dKSide();
423  double qlong = aPair->dKLong();
424  double kv = sqrt(qout * qout + qside * qside + qlong * qlong);
425  if (kv >= mkmax || kv < mkmin) return; //yyang
426 
427  if (fR2factor) weight *= ( 1.0 / (kv * kv) );
428 
429  StHbtYlm::YlmUpToL(elsi[GetMaxJM() - 1], qout, qside, qlong, ylmbuffer);
430  //check is not a number
431  for (int ilm = 0; ilm < GetMaxJM(); ilm++) {
432  if ( TMath::IsNaN(real(ylmbuffer[ilm]) * weight) ) return;
433  if ( TMath::IsNaN(imag(ylmbuffer[ilm]) * weight) ) return;
434  }
435  for (int ilmzero = 0; ilmzero < GetMaxJM(); ilmzero++)
436  for (int ilmprim = 0; ilmprim < GetMaxJM(); ilmprim++) {
437  if ( TMath::IsNaN( real(ylmbuffer[ilmzero]) * real(ylmbuffer[ilmprim]) ) ) return;
438  if ( TMath::IsNaN( real(ylmbuffer[ilmzero]) * -imag(ylmbuffer[ilmprim]) ) ) return;
439  if ( TMath::IsNaN( imag(ylmbuffer[ilmzero]) * real(ylmbuffer[ilmprim]) ) ) return;
440  if ( TMath::IsNaN( imag(ylmbuffer[ilmzero]) * -imag(ylmbuffer[ilmprim]) ) ) return;
441  }
442  for (int ilm = 0; ilm < GetMaxJM(); ilm++) {
443  densreal[ilm]->Fill(kv, real(ylmbuffer[ilm]) * weight);
444  densimag[ilm]->Fill(kv, -imag(ylmbuffer[ilm]) * weight);
445  binctd->Fill(kv, 1.0);
446  }
447 
448  //Fill in the error matrix
449  int nqbin = binctn->GetXaxis()->FindFixBin(kv) - 1;
450  //int tabshift = nqbin*GetMaxJM()*GetMaxJM()*4;
451  if ( nqbin < binctn->GetNbinsX() )
452  for (int ilmzero = 0; ilmzero < GetMaxJM(); ilmzero++)
453  for (int ilmprim = 0; ilmprim < GetMaxJM(); ilmprim++) {
454  covmden[GetBin(nqbin, ilmzero, 0, ilmprim, 0)] += real(ylmbuffer[ilmzero]) * real(ylmbuffer[ilmprim]);
455  covmden[GetBin(nqbin, ilmzero, 0, ilmprim, 1)] += real(ylmbuffer[ilmzero]) * -imag(ylmbuffer[ilmprim]);
456  covmden[GetBin(nqbin, ilmzero, 1, ilmprim, 0)] += -imag(ylmbuffer[ilmzero]) * real(ylmbuffer[ilmprim]);
457  covmden[GetBin(nqbin, ilmzero, 1, ilmprim, 1)] += -imag(ylmbuffer[ilmzero]) * -imag(ylmbuffer[ilmprim]);
458  }
459 }
460 
461 //________________________________________________________
462 void StHbtCorrFctnDirectYlm::Finish()
463 {
464  PackCovariances();
465 //CalcCorrFctn();
466 
467  for (int ilm = 0; ilm < maxjm; ilm++) {
468  numsreal[ilm]->Write();
469  numsimag[ilm]->Write();
470  densreal[ilm]->Write();
471  densimag[ilm]->Write();
472  if (cftnreal && cftnreal[ilm]) cftnreal[ilm]->Write();
473  if (cftnimag && cftnimag[ilm]) cftnimag[ilm]->Write();
474  if (mDoEMCIC) {
475  emcicP1P2T[ilm]->Write();
476  emcicP1P2Z[ilm]->Write();
477  emcicE1E2[ilm]->Write();
478  emcicE1plusE2[ilm]->Write();
479  }
480  }
481  if (covnum) covnum->Write();
482  if (covden) covden->Write();
483 }
484 
485 //_________________________________________________________________
486 void StHbtCorrFctnDirectYlm::Write(TFile* rfile)
487 {
488  rfile->cd();
489  for (int ilm = 0; ilm < maxjm; ilm++) {
490  numsreal[ilm]->Write();
491  numsimag[ilm]->Write();
492  densreal[ilm]->Write();
493  densimag[ilm]->Write();
494  if (cftnreal && cftnreal[ilm]) cftnreal[ilm]->Write();
495  if (cftnimag && cftnimag[ilm]) cftnimag[ilm]->Write();
496  if (mDoEMCIC) {
497  emcicP1P2T[ilm]->Write();
498  emcicP1P2Z[ilm]->Write();
499  emcicE1E2[ilm]->Write();
500  emcicE1plusE2[ilm]->Write();
501  }
502  }
503  if (covnum) covnum->Write();
504  if (covden) covden->Write();
505  rfile->Write();
506 }
507 
508 void StHbtCorrFctnDirectYlm::ReadFromFile(TFile* ifile)
509 {
510  /* yyang
511  * this function is used for read hadded */
512 
513  //cout << "Reading in numerators and denominators" << endl;
514  char bufname[200];
515  for (int ihist = 0; ihist < maxjm; ihist++) {
516  int em = emsi[ihist] < 0 ? elsi[ihist] - emsi[ihist] : emsi[ihist];
517  sprintf(bufname, "NumReYlm%i%i%s", elsi[ihist], em, mName.c_str());
518  if (numsreal[ihist]) delete numsreal[ihist];
519  numsreal[ihist] = new TH1D( *( (TH1D*) ifile->Get(bufname) ) );
520 
521  sprintf(bufname, "NumImYlm%i%i%s", elsi[ihist], em, mName.c_str());
522  if (numsimag[ihist]) delete numsimag[ihist];
523  numsimag[ihist] = new TH1D( *( (TH1D*) ifile->Get(bufname) ) );
524 
525  sprintf(bufname, "DenReYlm%i%i%s", elsi[ihist], em, mName.c_str());
526  if (densreal[ihist]) delete densreal[ihist];
527  densreal[ihist] = new TH1D( *( (TH1D*) ifile->Get(bufname) ) );
528 
529  sprintf(bufname, "DenImYlm%i%i%s", elsi[ihist], em, mName.c_str());
530  if (densimag[ihist]) delete densimag[ihist];
531  densimag[ihist] = new TH1D( *( (TH1D*) ifile->Get(bufname) ) );
532 
533  sprintf(bufname, "CfnReYlm%i%i%s", elsi[ihist], em, mName.c_str());
534  string keyname(bufname);
535  if (ifile->FindKey( (keyname + ";*").c_str() ) ) ifile->Delete( (keyname + ";*").c_str() );
536 
537  sprintf(bufname, "CfnImYlm%i%i%s", elsi[ihist], em, mName.c_str());
538  keyname = string(bufname);
539  if (ifile->FindKey( (keyname + ";*").c_str() ) ) ifile->Delete( (keyname + ";*").c_str() );
540  }
541 
542  if (covnum) delete covnum;
543  sprintf(bufname, "CovNum%s", mName.c_str());
544  covnum = new TH3D ( *( (TH3D*) ifile->Get(bufname) ) );
545 
546  if (covden) delete covden;
547  sprintf(bufname, "CovDen%s", mName.c_str() );
548  covden = new TH3D ( *( (TH3D*) ifile->Get(bufname) ) );
549 
550  if ( (covnum) && (covden) ) {
551  //cout << "Unpacking covariance matrices from file " << endl;
552  UnpackCovariances();
553  }
554  else {
555  cout << "Creating fake covariance matrices" << endl;
556 
557  for (int ibin = 1; ibin <= numsreal[0]->GetNbinsX(); ibin++) {
558  double nent = numsreal[0]->GetEntries();
559  double nentd = densreal[0]->GetEntries();
560  for (int ilmx = 0; ilmx < GetMaxJM(); ilmx++) {
561  for (int ilmy = 0; ilmy < GetMaxJM(); ilmy++) {
562  double t1t2rr = numsreal[ilmx]->GetBinContent(ibin) * numsreal[ilmy]->GetBinContent(ibin) / nent / nent;
563  double t1t2ri = numsreal[ilmx]->GetBinContent(ibin) * numsimag[ilmy]->GetBinContent(ibin) / nent / nent;
564  double t1t2ir = numsimag[ilmx]->GetBinContent(ibin) * numsreal[ilmy]->GetBinContent(ibin) / nent / nent;
565  double t1t2ii = numsimag[ilmx]->GetBinContent(ibin) * numsimag[ilmy]->GetBinContent(ibin) / nent / nent;
566  if (ilmx == ilmy) {
567  covmnum[GetBin(ibin - 1, ilmx, 0, ilmy, 0)] = nent * (TMath::Power(numsreal[ilmx]->GetBinError(ibin) / nent, 2) * (nent - 1) + t1t2rr);
568  covmnum[GetBin(ibin - 1, ilmx, 0, ilmy, 1)] = nent * t1t2ri;
569  covmnum[GetBin(ibin - 1, ilmx, 1, ilmy, 0)] = nent * t1t2ir;
570  covmnum[GetBin(ibin - 1, ilmx, 1, ilmy, 1)] = nent * (TMath::Power(numsimag[ilmx]->GetBinError(ibin) / nent, 2) * (nent - 1) + t1t2rr);
571  }
572  else {
573  covmnum[GetBin(ibin - 1, ilmx, 0, ilmy, 0)] = nent * t1t2rr;
574  covmnum[GetBin(ibin - 1, ilmx, 0, ilmy, 1)] = nent * t1t2ri;
575  covmnum[GetBin(ibin - 1, ilmx, 1, ilmy, 0)] = nent * t1t2ir;
576  covmnum[GetBin(ibin - 1, ilmx, 1, ilmy, 1)] = nent * t1t2ii;
577  }
578  t1t2rr = densreal[ilmx]->GetBinContent(ibin) * densreal[ilmy]->GetBinContent(ibin) / nentd / nentd;
579  t1t2ri = densreal[ilmx]->GetBinContent(ibin) * densimag[ilmy]->GetBinContent(ibin) / nentd / nentd;
580  t1t2ir = densimag[ilmx]->GetBinContent(ibin) * densreal[ilmy]->GetBinContent(ibin) / nentd / nentd;
581  t1t2ii = densimag[ilmx]->GetBinContent(ibin) * densimag[ilmy]->GetBinContent(ibin) / nentd / nentd;
582 
583  covmden[GetBin(ibin - 1, ilmx, 0, ilmy, 0)] = nentd * t1t2rr;
584  covmden[GetBin(ibin - 1, ilmx, 0, ilmy, 1)] = nentd * t1t2ri;
585  covmden[GetBin(ibin - 1, ilmx, 1, ilmy, 0)] = nentd * t1t2ir;
586  covmden[GetBin(ibin - 1, ilmx, 1, ilmy, 1)] = nentd * t1t2ii;
587  }
588  }
589  }
590  }
591 
592  //Recalculating the correlation functions
593  CalcCorrFctn();
594 
595  //save cf
596  ifile->cd();
597  for (int ilm = 0; ilm < maxjm; ilm++) {
598  cftnreal[ilm]->Write();
599  cftnimag[ilm]->Write();
600  }
601 }
602 
603 int StHbtCorrFctnDirectYlm::PackYlmVector(double* invec, double* outvec)
604 {
605  int ioutcount = 0;
606  int em, el;
607  for (int ilm = 0; ilm < GetMaxJM(); ilm++) {
608  GetElEmForIndex(ilm, &el, &em);
609  outvec[ioutcount++] = invec[ilm * 2];
610  if (em == 0) continue;
611  outvec[ioutcount++] = invec[ilm * 2 + 1];
612  }
613 
614  return (ioutcount);
615 }
616 
617 int StHbtCorrFctnDirectYlm::PackYlmMatrix(double* inmat, double* outmat)
618 {
619  int ioutcountz = 0;
620  int ioutcountp = 0;
621  int emz, elz;
622  int emp, elp;
623  int finalsize = 0;
624 
625  for (int ilm = 0; ilm < GetMaxJM(); ilm++) {
626  GetElEmForIndex(ilm, &elz, &emz);
627  finalsize++;
628  if (emz == 0) continue;
629  finalsize++;
630  }
631 
632  for (int ilmz = 0; ilmz < GetMaxJM(); ilmz++) {
633  GetElEmForIndex(ilmz, &elz, &emz);
634  ioutcountp = 0;
635  for (int ilmp = 0; ilmp < GetMaxJM(); ilmp++) {
636  GetElEmForIndex(ilmp, &elp, &emp);
637  outmat[ioutcountz * finalsize + ioutcountp] = inmat[GetBin(0, ilmz, 0, ilmp, 0)];
638  ioutcountp++;
639  if (emp == 0) continue;
640  outmat[ioutcountz * finalsize + ioutcountp] = inmat[GetBin(0, ilmz, 0, ilmp, 1)];
641  ioutcountp++;
642  }
643  ioutcountz++;
644 
645  if (emz == 0) continue;
646  ioutcountp = 0;
647  for (int ilmp = 0; ilmp < GetMaxJM(); ilmp++) {
648  GetElEmForIndex(ilmp, &elp, &emp);
649  outmat[ioutcountz * finalsize + ioutcountp] = inmat[GetBin(0, ilmz, 1, ilmp, 0)];
650  ioutcountp++;
651  if (emp == 0) continue;
652  outmat[ioutcountz * finalsize + ioutcountp] = inmat[GetBin(0, ilmz, 1, ilmp, 1)];
653  ioutcountp++;
654  }
655  ioutcountz++;
656  }
657 
658  return (ioutcountz);
659 }
660 
661 void StHbtCorrFctnDirectYlm::PackCovariances()
662 {
663  char bufname[200];
664 
665  sprintf(bufname, "CovNum%s", numsreal[0]->GetName() + 10);
666  if (covnum) delete covnum;
667  covnum = new TH3D( bufname, bufname,
668  numsreal[0]->GetNbinsX(), numsreal[0]->GetXaxis()->GetXmin(), numsreal[0]->GetXaxis()->GetXmax(),
669  GetMaxJM() * 2, -0.5, GetMaxJM() * 2 - 0.5,
670  GetMaxJM() * 2, -0.5, GetMaxJM() * 2 - 0.5);
671 
672  for (int ibin = 1; ibin <= covnum->GetNbinsX(); ibin++)
673  for (int ilmz = 0; ilmz < GetMaxJM() * 2; ilmz++)
674  for (int ilmp = 0; ilmp < GetMaxJM() * 2; ilmp++)
675  covnum->SetBinContent(ibin, ilmz + 1, ilmp + 1, covmnum[GetBin(ibin - 1, ilmz / 2, ilmz % 2, ilmp / 2, ilmp % 2)]);
676 
677  sprintf(bufname, "CovDen%s", numsreal[0]->GetName() + 10);
678  if (covden) delete covden;
679  covden = new TH3D( bufname, bufname,
680  densreal[0]->GetNbinsX(), densreal[0]->GetXaxis()->GetXmin(), densreal[0]->GetXaxis()->GetXmax(),
681  GetMaxJM() * 2, -0.5, GetMaxJM() * 2 - 0.5,
682  GetMaxJM() * 2, -0.5, GetMaxJM() * 2 - 0.5);
683 
684  for (int ibin = 1; ibin <= covden->GetNbinsX(); ibin++)
685  for (int ilmz = 0; ilmz < GetMaxJM() * 2; ilmz++)
686  for (int ilmp = 0; ilmp < GetMaxJM() * 2; ilmp++)
687  covden->SetBinContent(ibin, ilmz + 1, ilmp + 1, covmden[GetBin(ibin - 1, ilmz / 2, ilmz % 2, ilmp / 2, ilmp % 2)]);
688 
689 }
690 
691 void StHbtCorrFctnDirectYlm::UnpackCovariances()
692 {
693  if (covnum) {
694  for (int ibin = 1; ibin <= covnum->GetNbinsX(); ibin++)
695  for (int ilmz = 0; ilmz < GetMaxJM() * 2; ilmz++)
696  for (int ilmp = 0; ilmp < GetMaxJM() * 2; ilmp++)
697  covmnum[GetBin(ibin - 1, ilmz / 2, ilmz % 2, ilmp / 2, ilmp % 2)] = covnum->GetBinContent(ibin, ilmz + 1, ilmp + 1);
698  }
699  if (covden) {
700  for (int ibin = 1; ibin <= covden->GetNbinsX(); ibin++)
701  for (int ilmz = 0; ilmz < GetMaxJM() * 2; ilmz++)
702  for (int ilmp = 0; ilmp < GetMaxJM() * 2; ilmp++)
703  covmden[GetBin(ibin - 1, ilmz / 2, ilmz % 2, ilmp / 2, ilmp % 2)] = covden->GetBinContent(ibin, ilmz + 1, ilmp + 1);
704  }
705 }
706 
707 int StHbtCorrFctnDirectYlm::GetIndexForLM(int el, int em)
708 {
709  for (int iter = 0; iter < maxjm; iter++)
710  if ( (el == elsi[iter]) && (em == emsi[iter]) )
711  return (iter);
712  return (-1);
713 }
714 
715 TH1D* StHbtCorrFctnDirectYlm::GetNumRealHist(int el, int em)
716 {
717  if (GetIndexForLM(el, em) >= 0)
718  return (numsreal[GetIndexForLM(el, em)]);
719  else
720  return (0);
721 }
722 
723 TH1D* StHbtCorrFctnDirectYlm::GetNumImagHist(int el, int em)
724 {
725  if (GetIndexForLM(el, em) >= 0)
726  return (numsimag[GetIndexForLM(el, em)]);
727  else
728  return (0);
729 }
730 
731 TH1D* StHbtCorrFctnDirectYlm::GetDenRealHist(int el, int em)
732 {
733  if (GetIndexForLM(el, em) >= 0)
734  return (densreal[GetIndexForLM(el, em)]);
735  else
736  return (0);
737 }
738 
739 TH1D* StHbtCorrFctnDirectYlm::GetDenImagHist(int el, int em)
740 {
741  if (GetIndexForLM(el, em) >= 0)
742  return (densimag[GetIndexForLM(el, em)]);
743  else
744  return (0);
745 }
746 
747 StHbtString StHbtCorrFctnDirectYlm::Report()
748 {
749  return ("StHbtCorrFctnDirectYlm::Finish");
750 }
751 
752 void StHbtCorrFctnDirectYlm::SetR2Factor(bool aFactor)
753 {
754  fR2factor = aFactor;
755 }
756 
757 //__________________________________________________________________________
758 void StHbtCorrFctnDirectYlm::CalcCorrFctn()
759 {
760  cftnreal = new TH1D*[maxjm];
761  cftnimag = new TH1D*[maxjm];
762 
763  char bufname[200];
764  for (int ihist = 0; ihist < maxjm; ihist++) {
765  int em = emsi[ihist] < 0 ? elsi[ihist] - emsi[ihist] : emsi[ihist];
766 
767  sprintf( bufname, "CfnReYlm%i%i%s", elsi[ihist], em, mName.c_str() );
768  cftnreal[ihist] = new TH1D(bufname, bufname, mNbins, mkmin, mkmax);
769  sprintf( bufname, "CfnImYlm%i%i%s", elsi[ihist], em, mName.c_str() );
770  cftnimag[ihist] = new TH1D(bufname, bufname, mNbins, mkmin, mkmax);
771 
772  cftnreal[ihist]->Sumw2();
773  cftnimag[ihist]->Sumw2();
774  }
775 
776  covmcfc = new double[maxjm * maxjm * 4 * mNbins];
777 
778  complex<double> tMq0[maxjm];
779  complex<double> tTq0[maxjm];
780  double tMTilde[maxjm * maxjm * 4];
781  complex<double> tCq0[maxjm];
782 
783  int recalccov = 1;
784  if ( (covnum) && (covnum->GetBinContent(0, 0, 0) > 0.0) ) {
785  cout << "Detected calculated covariance matrix. Do not recalculate !!!\n";
786  recalccov = 0;
787  }
788 
789  //If requested calculate the advanced normalization factor
790  //*** WARNING !!! ***
791  //This is the factor that assumes that the correlation function
792  //calculated is the non-identical particle correlation function
793  //dominated by the asymptotic behaviour of the Coulmb interaction
794  //in the normalization region
795  //For any other function use normal normalization !!!
796  //*** WARNING !!! ***
797 
798  double normfactor = 1.0;
799  if (mNormBinMax > 0) {
800  double sksum = 0.0;
801  double wksum = 0.0;
802 
803  double sk, wk, ks;
804  if (mNormBinMin < 1) mNormBinMin = 1;
805  if ( mNormBinMax > densreal[0]->GetNbinsX() )
806  mNormBinMax = densreal[0]->GetNbinsX();
807 
808  for (int ib = mNormBinMin; ib <= mNormBinMax; ib++) {
809  ks = densreal[0]->GetXaxis()->GetBinCenter(ib);
810  sk = numsreal[0]->GetBinContent(ib) / ( densreal[0]->GetBinContent(ib) * ( 1.0 - mNormPurity / (mNormRadius * mNormBohr * ks * ks) ) );
811  wk = numsreal[0]->GetBinContent(ib);
812  sksum += sk * wk;
813  wksum += wk;
814  }
815  normfactor *= sksum / wksum;
816  normfactor /= numsreal[0]->GetEntries() / densreal[0]->GetEntries();
817  }
818 
819  for (int ibin = 1; ibin <= numsreal[0]->GetNbinsX(); ibin++) {
820  for (int ilm = 0; ilm < maxjm; ilm++) {
821  if (recalccov) {
822  tMq0[ilm] = complex<double>( densreal[ilm]->GetBinContent(ibin) / (densreal[0]->GetEntries() / normfactor),
823  densimag[ilm]->GetBinContent(ibin) / (densreal[0]->GetEntries() / normfactor) );
824  tTq0[ilm] = complex<double>( numsreal[ilm]->GetBinContent(ibin) / numsreal[0]->GetEntries(),
825  numsimag[ilm]->GetBinContent(ibin) / numsreal[0]->GetEntries() );
826  }
827  else {
828  tMq0[ilm] = complex<double>(densreal[ilm]->GetBinContent(ibin) / normfactor,
829  densimag[ilm]->GetBinContent(ibin) / normfactor);
830  tTq0[ilm] = complex<double>( numsreal[ilm]->GetBinContent(ibin),
831  numsimag[ilm]->GetBinContent(ibin) );
832  }
833  }
834 
835  //Calculate the proper error matrix for T
836  //from the temporary covariance matrices
837  if (recalccov) {
838  for (int ilmzero = 0; ilmzero < GetMaxJM(); ilmzero++) {
839  for (int ilmprim = 0; ilmprim < GetMaxJM(); ilmprim++) {
840  if ( ::isnan(covmnum[GetBin(ibin - 1, ilmzero, 0, ilmprim, 0)]) ) {
841  cout << "NaN !!!! RR " << ilmzero << " " << ilmprim << endl;
842  }
843  if ( ::isnan(covmnum[GetBin(ibin - 1, ilmzero, 0, ilmprim, 1)]) ) {
844  cout << "NaN !!!! RI " << ilmzero << " " << ilmprim << endl;
845  }
846  if ( ::isnan(covmnum[GetBin(ibin - 1, ilmzero, 1, ilmprim, 0)]) ) {
847  cout << "NaN !!!! IR " << ilmzero << " " << ilmprim << endl;
848  }
849  if ( ::isnan(covmnum[GetBin(ibin - 1, ilmzero, 1, ilmprim, 1)]) ) {
850  cout << "NaN !!!! II " << ilmzero << " " << ilmprim << endl;
851  }
852 
853  covmnum[GetBin(ibin - 1, ilmzero, 0, ilmprim, 0)] /= numsreal[0]->GetEntries();
854  covmnum[GetBin(ibin - 1, ilmzero, 0, ilmprim, 1)] /= numsreal[0]->GetEntries();
855  covmnum[GetBin(ibin - 1, ilmzero, 1, ilmprim, 0)] /= numsreal[0]->GetEntries();
856  covmnum[GetBin(ibin - 1, ilmzero, 1, ilmprim, 1)] /= numsreal[0]->GetEntries();
857 
858  if ( ::isnan(covmnum[GetBin(ibin - 1, ilmzero, 0, ilmprim, 0)]) ) {
859  cout << "NaN !!!! RR" << endl;
860  }
861  if ( ::isnan(covmnum[GetBin(ibin - 1, ilmzero, 0, ilmprim, 1)]) ) {
862  cout << "NaN !!!! RI" << endl;
863  }
864  if ( ::isnan(covmnum[GetBin(ibin - 1, ilmzero, 1, ilmprim, 0)]) ) {
865  cout << "NaN !!!! IR" << endl;
866  }
867  if ( ::isnan(covmnum[GetBin(ibin - 1, ilmzero, 1, ilmprim, 1)]) ) {
868  cout << "NaN !!!! II" << endl;
869  }
870 
871  covmnum[GetBin(ibin - 1, ilmzero, 0, ilmprim, 0)] -= real(tTq0[ilmzero]) * real(tTq0[ilmprim]);
872  covmnum[GetBin(ibin - 1, ilmzero, 0, ilmprim, 1)] -= real(tTq0[ilmzero]) * imag(tTq0[ilmprim]);
873  covmnum[GetBin(ibin - 1, ilmzero, 1, ilmprim, 0)] -= imag(tTq0[ilmzero]) * real(tTq0[ilmprim]);
874  covmnum[GetBin(ibin - 1, ilmzero, 1, ilmprim, 1)] -= imag(tTq0[ilmzero]) * imag(tTq0[ilmprim]);
875 
876  if ( ::isnan(covmnum[GetBin(ibin - 1, ilmzero, 0, ilmprim, 0)]) ) {
877  cout << "NaN !!!! RR" << endl;
878  }
879  if ( ::isnan(covmnum[GetBin(ibin - 1, ilmzero, 0, ilmprim, 1)]) ) {
880  cout << "NaN !!!! RI" << endl;
881  }
882  if ( ::isnan(covmnum[GetBin(ibin - 1, ilmzero, 1, ilmprim, 0)]) ) {
883  cout << "NaN !!!! IR" << endl;
884  }
885  if ( ::isnan(covmnum[GetBin(ibin - 1, ilmzero, 1, ilmprim, 1)]) ) {
886  cout << "NaN !!!! II" << endl;
887  }
888 
889  covmnum[GetBin(ibin - 1, ilmzero, 0, ilmprim, 0)] /= (numsreal[0]->GetEntries() - 1);
890  covmnum[GetBin(ibin - 1, ilmzero, 0, ilmprim, 1)] /= (numsreal[0]->GetEntries() - 1);
891  covmnum[GetBin(ibin - 1, ilmzero, 1, ilmprim, 0)] /= (numsreal[0]->GetEntries() - 1);
892  covmnum[GetBin(ibin - 1, ilmzero, 1, ilmprim, 1)] /= (numsreal[0]->GetEntries() - 1);
893 
894  if ( ::isnan(covmnum[GetBin(ibin - 1, ilmzero, 0, ilmprim, 0)]) ) {
895  cout << "NaN !!!! RR" << endl;
896  }
897  if ( ::isnan(covmnum[GetBin(ibin - 1, ilmzero, 0, ilmprim, 1)]) ) {
898  cout << "NaN !!!! RI" << endl;
899  }
900  if ( ::isnan(covmnum[GetBin(ibin - 1, ilmzero, 1, ilmprim, 0)]) ) {
901  cout << "NaN !!!! IR" << endl;
902  }
903  if ( ::isnan(covmnum[GetBin(ibin - 1, ilmzero, 1, ilmprim, 1)]) ) {
904  cout << "NaN !!!! II" << endl;
905  }
906  }
907  }
908  }
909  GetMtilde(tMq0, tMTilde);
910 
911  //Perform the solution for the correlation function itself and the errors
912  double mDeltaT[maxjm * maxjm * 4];
913  for (int ilmzero = 0; ilmzero < GetMaxJM() * 2; ilmzero++)
914  for (int ilmprim = 0; ilmprim < GetMaxJM() * 2; ilmprim++)
915  mDeltaT[(ilmzero * maxjm * 2) + ilmprim] = (covmnum[GetBin(ibin - 1, ilmzero / 2, ilmzero % 2, ilmprim / 2, ilmprim % 2)]);
916 
917 #ifdef YY_DEBUG
918  cout << "Delta T matrix " << endl;
919  for (int ilmz = 0; ilmz < GetMaxJM() * 2; ilmz++) {
920  for (int ilmp = 0; ilmp < GetMaxJM() * 2; ilmp++) {
921  cout.precision(3);
922  cout.width(10);
923  cout << mDeltaT[ilmz * GetMaxJM() * 2 + ilmp];
924  }
925  cout << endl;
926  }
927 #endif
928 
929  double mDeltaTPacked[maxjm * maxjm * 4];
930  int msize = PackYlmMatrixIndependentOnly(mDeltaT, mDeltaTPacked);
931 
932 #ifdef YY_DEBUG
933  cout << "Delta T matrix packed " << endl;
934  for (int ilmz = 0; ilmz < msize; ilmz++) {
935  for (int ilmp = 0; ilmp < msize; ilmp++) {
936  cout.precision(3);
937  cout.width(10);
938  cout << mDeltaTPacked[ilmz * msize + ilmp];
939  }
940  cout << endl;
941  }
942 #endif
943 
944  //(1) Solve (DeltaT)^1 Mtilde = Q
945  //Prepare halper matrices
946  double mM[maxjm * maxjm * 4];
947  double mMPacked[maxjm * maxjm * 4];
948  for (int iter = 0; iter < maxjm * maxjm * 4; iter++)
949  mM[iter] = tMTilde[iter];
950  PackYlmMatrixIndependentOnly(mM, mMPacked);
951 
952  gsl_matrix_view matM = gsl_matrix_view_array(mMPacked, msize, msize);
953 
954 #ifdef YY_DEBUG
955  cout << "Mtilde matrix " << endl;
956  for (int ilmz = 0; ilmz < GetMaxJM() * 2; ilmz++) {
957  for (int ilmp = 0; ilmp < GetMaxJM() * 2; ilmp++) {
958  cout.precision(3);
959  cout.width(10);
960  cout << mM[ilmz * GetMaxJM() * 2 + ilmp];
961  }
962  cout << endl;
963  }
964  cout << "Mtilde matrix packed " << endl;
965  for (int ilmz = 0; ilmz < msize; ilmz++) {
966  for (int ilmp = 0; ilmp < msize; ilmp++) {
967  cout.precision(3);
968  cout.width(10);
969  cout << mMPacked[ilmz * msize + ilmp];
970  }
971  cout << endl;
972  }
973 #endif
974 
975  //Inverting matrix DeltaT.
976  double mU[maxjm * maxjm * 4];
977  InvertYlmIndependentMatrix(mDeltaT, mU);
978 
979  double mDTInvertedPacked[maxjm * maxjm * 4];
980  PackYlmMatrixIndependentOnly(mU, mDTInvertedPacked);
981 
982  gsl_matrix_view matDTI = gsl_matrix_view_array(mDTInvertedPacked, msize, msize);
983 #ifdef YY_DEBUG
984  cout << "Delta T matrix inverted packed " << endl;
985  for (int ilmz = 0; ilmz < msize; ilmz++) {
986  for (int ilmp = 0; ilmp < msize; ilmp++) {
987  cout.precision(3);
988  cout.width(10);
989  cout << mDTInvertedPacked[ilmz * msize + ilmp];
990  }
991  cout << endl;
992  }
993 #endif
994 
995  //(2) Multiply DeltaT^1 M = Q
996  double mQ[maxjm * maxjm * 4];
997  for (int iter = 0; iter < msize * msize; iter++)
998  mQ[iter] = 0.0;
999  gsl_matrix_view matQ = gsl_matrix_view_array(mQ, msize, msize);
1000  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, &matDTI.matrix, &matM.matrix, 0.0, &matQ.matrix);
1001 
1002  double mTest[maxjm * maxjm * 4];
1003  gsl_matrix_view matTest = gsl_matrix_view_array(mTest, msize, msize);
1004 
1005  double mF[maxjm * maxjm * 4];
1006  for (int iter = 0; iter < maxjm * maxjm * 4; iter++)
1007  mF[iter] = mDeltaTPacked[iter];
1008  gsl_matrix_view matF = gsl_matrix_view_array(mF, msize, msize);
1009  gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, &matF.matrix, &matQ.matrix, 0.0, &matTest.matrix);
1010 #ifdef YY_DEBUG
1011  cout << "Test matrix packed - compare to Mtilde" << endl;
1012  for (int ilmz = 0; ilmz < msize; ilmz++) {
1013  for (int ilmp = 0; ilmp < msize; ilmp++) {
1014  cout.precision(3);
1015  cout.width(10);
1016  cout << mTest[ilmz * msize + ilmp];
1017  }
1018  cout << endl;
1019  }
1020 #endif
1021  //(2) Multiply Mtilde^T Q = P
1022  double mP[maxjm * maxjm * 4];
1023  for (int iter = 0; iter < maxjm * maxjm * 4; iter++)
1024  mP[iter] = 0;
1025 
1026  gsl_matrix_view matP = gsl_matrix_view_array(mP, msize, msize);
1027  gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, &matM.matrix, &matQ.matrix, 0.0, &matP.matrix);
1028 
1029 #ifdef YY_DEBUG
1030  cout << "P matrix packed" << endl;
1031  for (int ilmz = 0; ilmz < msize; ilmz++) {
1032  for (int ilmp = 0; ilmp < msize; ilmp++) {
1033  cout.precision(3);
1034  cout.width(10);
1035  cout << mP[ilmz * msize + ilmp];
1036  }
1037  cout << endl;
1038  }
1039 #endif
1040 
1041  //(3) Solve P^-1 Mtilde^T = R
1042  double mPUnpacked[maxjm * maxjm * 4];
1043  UnPackYlmMatrixIndependentOnly(mP, mPUnpacked, msize);
1044 
1045 #ifdef YY_DEBUG
1046  cout << "P matrix unpacked" << endl;
1047  for (int ilmz = 0; ilmz < GetMaxJM() * 2; ilmz++) {
1048  for (int ilmp = 0; ilmp < GetMaxJM() * 2; ilmp++) {
1049  cout.precision(3);
1050  cout.width(10);
1051  cout << mPUnpacked[ilmz * GetMaxJM() * 2 + ilmp];
1052  }
1053  cout << endl;
1054  }
1055 #endif
1056 
1057  //Invert the P matrix
1058  double mPInverted[maxjm * maxjm * 4];
1059  InvertYlmIndependentMatrix(mPUnpacked, mPInverted);
1060 
1061  double mPInvertedPacked[maxjm * maxjm * 4];
1062  PackYlmMatrixIndependentOnly(mPInverted, mPInvertedPacked);
1063 
1064  gsl_matrix_view matPI = gsl_matrix_view_array(mPInvertedPacked, msize, msize);
1065 
1066 #ifdef YY_DEBUG
1067  cout << "P matrix inverted packed" << endl;
1068  for (int ilmz = 0; ilmz < msize; ilmz++) {
1069  for (int ilmp = 0; ilmp < msize; ilmp++) {
1070  cout.precision(3);
1071  cout.width(10);
1072  cout << mPInvertedPacked[ilmz * msize + ilmp];
1073  }
1074  cout << endl;
1075  }
1076 #endif
1077 
1078  double mR[maxjm * maxjm * 4];
1079  for (int ir = 0; ir < maxjm * maxjm * 4; ir++)
1080  mR[ir] = 0.0;
1081  gsl_matrix_view matR = gsl_matrix_view_array(mR, msize, msize);
1082 
1083 
1084  //(2) Multiply P^-1 M (Trans) = R
1085 #ifdef YY_DEBUG
1086  cout << "Matrix M Packed " << endl;
1087  for (int ilmz = 0; ilmz < msize; ilmz++) {
1088  for (int ilmp = 0; ilmp < msize; ilmp++) {
1089  cout.precision(3);
1090  cout.width(10);
1091  cout << mMPacked[ilmz * msize + ilmp];
1092  }
1093  cout << endl;
1094  }
1095 #endif
1096 
1097  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, &matPI.matrix, &matM.matrix, 1.0, &matR.matrix);
1098 
1099 #ifdef YY_DEBUG
1100  cout << "R matrix packed " << endl;
1101  for (int ilmz = 0; ilmz < msize; ilmz++) {
1102  for (int ilmp = 0; ilmp < msize; ilmp++) {
1103  cout.precision(3);
1104  cout.width(10);
1105  cout << mR[ilmz * msize + ilmp];
1106  }
1107  cout << endl;
1108  }
1109 #endif
1110 
1111  //(4) Solve DeltaT^-1 T = L
1112  double vL[maxjm * 2];
1113  gsl_vector_view vecL = gsl_vector_view_array(vL, msize);
1114 
1115  //Decomposing the M matrix
1116  //gsl_linalg_SV_decomp(&matF.matrix, &matS.matrix, &vecST.vector, &vecWT.vector);
1117 
1118  double vB[maxjm * 2];
1119  for (int iter = 0; iter < GetMaxJM(); iter++) {
1120  vB[iter * 2] = real(tTq0[iter]);
1121  vB[iter * 2 + 1] = imag(tTq0[iter]);
1122  }
1123 
1124  double vBPacked[maxjm * 2];
1125  PackYlmVectorIndependentOnly(vB, vBPacked);
1126 
1127  gsl_vector_view vecB = gsl_vector_view_array(vBPacked, msize);
1128 
1129 
1130  //Solving the problem
1131  //gsl_linalg_SV_solve(&matF.matrix, &matS.matrix, &vecST.vector, &vecB.vector, &vecL.vector);
1132 
1133 #ifdef YY_DEBUG
1134  cout << "L vector packed " << endl;
1135  for (int ilmp = 0; ilmp < msize; ilmp++) {
1136  cout.precision(3);
1137  cout.width(10);
1138  cout << vL[ilmp];
1139  }
1140  cout << endl;
1141 #endif
1142 
1143  //Multiply DeltaT^-1 T = L
1144  gsl_blas_dgemv(CblasNoTrans, 1.0, &matDTI.matrix, &vecB.vector, 0.0, &vecL.vector);
1145 
1146  //(5) Multiply R L = C
1147  double vY[maxjm * 2];
1148  for (int iter = 0; iter < GetMaxJM() * 2; iter++) {
1149  vY[iter] = 0.0;
1150  }
1151 
1152  //Prepare inputs for solving the problem
1153  gsl_vector_view vecY = gsl_vector_view_array(vY, msize);
1154  gsl_blas_dgemv (CblasNoTrans, 1.0, &matR.matrix, &vecL.vector, 0.0, &vecY.vector);
1155 
1156 #ifdef YY_DEBUG
1157  cout << "C vector packed" << endl;
1158  for (int ilmp = 0; ilmp < msize; ilmp++) {
1159  cout.precision(3);
1160  cout.width(10);
1161  cout << vY[ilmp];
1162  }
1163  cout << endl;
1164 #endif
1165 
1166 
1167  int mpack = 0;
1168  int el, em;
1169  for (int ilm = 0; ilm < maxjm; ilm++) {
1170  //cftnreal[ilm]->SetBinContent(ibin, vC[mpack++]);
1171  GetElEmForIndex(ilm, &el, &em);
1172  if (em < 0) {
1173  cftnreal[ilm]->SetBinContent(ibin, 0.0);
1174  cftnimag[ilm]->SetBinContent(ibin, 0.0);
1175  }
1176  else {
1177  cftnreal[ilm]->SetBinContent(ibin, vY[mpack++]);
1178  if (em == 0)
1179  cftnimag[ilm]->SetBinContent(ibin, 0);
1180  else
1181  //cftnimag[ilm]->SetBinContent(ibin, vC[mpack++]);
1182  cftnimag[ilm]->SetBinContent(ibin, vY[mpack++]);
1183  }
1184  }
1185 
1186  //Invert V
1187  mpack = 0;
1188  for (int ilm = 0; ilm < maxjm; ilm++) {
1189  GetElEmForIndex(ilm, &el, &em);
1190  if (em < 0 ) {
1191  cftnreal[ilm]->SetBinError(ibin, 0);
1192  cftnimag[ilm]->SetBinError(ibin, 0);
1193  }
1194  else {
1195  cftnreal[ilm]->SetBinError( ibin, sqrt( fabs(mPInvertedPacked[mpack * msize + mpack]) ) );
1196  mpack++;
1197  if (em == 0)
1198  cftnimag[ilm]->SetBinError(ibin, 0);
1199  else {
1200  cftnimag[ilm]->SetBinError( ibin, sqrt( fabs(mPInvertedPacked[mpack * msize + mpack]) ) );
1201  mpack++;
1202  }
1203  }
1204  }
1205 
1206  for (int ilmz = 0; ilmz < GetMaxJM() * 2; ilmz++) {
1207  for (int ilmp = 0; ilmp < GetMaxJM() * 2; ilmp++) {
1208  if (ilmp > ilmz)
1209  covmcfc[GetBin(ibin - 1, ilmz / 2, ilmz % 2, ilmp / 2, ilmp % 2)] = mPInverted[ilmz * GetMaxJM() * 2 + ilmp];
1210  else
1211  covmcfc[GetBin(ibin - 1, ilmz / 2, ilmz % 2, ilmp / 2, ilmp % 2)] = mPInverted[ilmp * GetMaxJM() * 2 + ilmz];
1212  }
1213  }
1214  } //for (int ibin=1; ibin<=numsreal[0]->GetNbinsX(); ibin++)
1215 
1216 
1217 //PackCfcCovariance()
1218 }
1219 
1220 //______________________________________________________________________
1221 int StHbtCorrFctnDirectYlm::PackYlmMatrixIndependentOnly(double* inmat, double* outmat)
1222 {
1223  int ioutcountz = 0;
1224  int ioutcountp = 0;
1225  int emz, elz;
1226  int emp, elp;
1227  int finalsize = 0;
1228 
1229  for (int ilm = 0; ilm < GetMaxJM(); ilm++) {
1230  GetElEmForIndex(ilm, &elz, &emz);
1231  if (emz < 0) continue;
1232  finalsize++;
1233  if (emz == 0) continue;
1234  finalsize++;
1235  }
1236  //cout << "Final size " << finalsize << endl;
1237 
1238  for (int ilmz = 0; ilmz < GetMaxJM(); ilmz++) {
1239  GetElEmForIndex(ilmz, &elz, &emz);
1240  ioutcountp = 0;
1241 
1242  if (emz < 0) continue;
1243  for (int ilmp = 0; ilmp < GetMaxJM(); ilmp++) {
1244  GetElEmForIndex(ilmp, &elp, &emp);
1245  if (emp < 0) continue;
1246  outmat[ioutcountz * finalsize + ioutcountp] = inmat[GetBin(0, ilmz, 0, ilmp, 0)];
1247  ioutcountp++;
1248  if (emp == 0) continue;
1249  outmat[ioutcountz * finalsize + ioutcountp] = inmat[GetBin(0, ilmz, 0, ilmp, 1)];
1250  ioutcountp++;
1251  }
1252  ioutcountz++;
1253 
1254  if (emz == 0) continue;
1255  ioutcountp = 0;
1256  for (int ilmp = 0; ilmp < GetMaxJM(); ilmp++) {
1257  GetElEmForIndex(ilmp, &elp, &emp);
1258  if (emp < 0) continue;
1259  outmat[ioutcountz * finalsize + ioutcountp] = inmat[GetBin(0, ilmz, 1, ilmp, 0)];
1260  ioutcountp++;
1261  if (emp == 0) continue;
1262  outmat[ioutcountz * finalsize + ioutcountp] = inmat[GetBin(0, ilmz, 1, ilmp, 1)];
1263  ioutcountp++;
1264  }
1265  ioutcountz++;
1266  }
1267  return (ioutcountz);
1268 }
1269 
1270 //_____________________________________________________________________
1271 void StHbtCorrFctnDirectYlm::InvertYlmIndependentMatrix(double* inmat, double* outmat)
1272 {
1273  //Invert the Ylm matrix by inverting only the matrix
1274  //with independent elements and filling in the rest
1275  //according to sign rules
1276 
1277  double mU[maxjm * maxjm * 4];
1278  int isize = PackYlmMatrixIndependentOnly(inmat, mU);
1279  //cout << "Independent count " << isize << endl;
1280 
1281  gsl_matrix_view matU = gsl_matrix_view_array(mU, isize, isize);
1282 
1283 #ifdef YY_DEBUG
1284  cout << "Input matrix independent only " << endl;
1285  for (int ilmz = 0; ilmz < isize; ilmz++) {
1286  for (int ilmp = 0; ilmp < isize; ilmp++) {
1287  cout.precision(3);
1288  cout.width(10);
1289  cout << mU[ilmz * isize + ilmp];
1290  }
1291  cout << endl;
1292  }
1293 #endif
1294 
1295  //Identity matrix helper for inversion
1296  double mI[maxjm * maxjm * 4];
1297  for (int iterm = 0; iterm < isize; iterm++)
1298  for (int iterp = 0; iterp < isize; iterp++)
1299  if (iterm == iterp)
1300  mI[iterm * isize + iterp] = 1.0;
1301  else
1302  mI[iterm * isize + iterp] = 0.0;
1303 
1304  gsl_matrix_view matI = gsl_matrix_view_array(mI, isize, isize);
1305 
1306  //Invert the matrix
1307  gsl_blas_dtrsm(CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, 1.0, &matU.matrix, &matI.matrix);
1308 
1309  UnPackYlmMatrixIndependentOnly(mI, outmat, isize);
1310 }
1311 
1312 //_____________________________________________________________________
1313 void StHbtCorrFctnDirectYlm::UnPackYlmMatrixIndependentOnly(double* inmat, double* outmat, int insize)
1314 {
1315  int lmax = sqrt(insize) - 1;
1316  //cout << "lmax is " << lmax << endl;
1317  int tmax = (lmax + 1) * (lmax + 1) * 2;
1318  int indexfrom[tmax];
1319  int multfrom[tmax];
1320 
1321  int el, em;
1322  for (int iter = 0; iter < tmax; iter++) {
1323  int im = iter % 2;
1324  GetElEmForIndex(iter / 2, &el, &em);
1325  if (em == 0) {
1326  if (im == 1) {
1327  indexfrom[iter] = 0;
1328  multfrom[iter] = 0;
1329  }
1330  else {
1331  indexfrom[iter] = el * el;
1332  multfrom[iter] = 1;
1333  }
1334  }
1335  else if (em < 0) {
1336  indexfrom[iter] = (el * el) + (-em) * 2 - 1;
1337  if (im) indexfrom[iter]++;
1338  if ( (-em) % 2 )
1339  if (im) multfrom[iter] = 1;
1340  else multfrom[iter] = -1;
1341  else
1342  if (im) multfrom[iter] = -1;
1343  else multfrom[iter] = 1;
1344  }
1345  else if (em > 0) {
1346  indexfrom[iter] = (el * el) + (em) * 2 - 1;
1347  if (im) indexfrom[iter]++;
1348  multfrom[iter] = 1;
1349  }
1350  }
1351 
1352  for (int ilmz = 0; ilmz < GetMaxJM() * 2; ilmz++)
1353  for (int ilmp = 0; ilmp < GetMaxJM() * 2; ilmp++)
1354  outmat[ilmz * GetMaxJM() * 2 + ilmp] = inmat[(indexfrom[ilmz] * insize) + indexfrom[ilmp]] * multfrom[ilmz] * multfrom[ilmp];
1355 }
1356 
1357 //____________________________________________________________________
1358 int StHbtCorrFctnDirectYlm::PackYlmVectorIndependentOnly(double* invec, double* outvec)
1359 {
1360  int ioutcount = 0;
1361  int em, el;
1362  for (int ilm = 0; ilm < GetMaxJM(); ilm++) {
1363  GetElEmForIndex(ilm, &el, &em);
1364  if (em < 0) continue;
1365  outvec[ioutcount++] = invec[ilm * 2];
1366  if (em == 0)
1367  continue;
1368  outvec[ioutcount++] = invec[ilm * 2 + 1];
1369  }
1370  return (ioutcount);
1371 }
1372 
1373 /***************************************************************************
1374  * $Log: StHbtCorrFctnDirectYlm.cxx,v $
1375  * Revision 1.2 2015/11/02 20:11:06 perev
1376  * isnan for new compiler
1377  *
1378  * Revision 1.1 2013/01/18 14:46:02 yyang
1379  * Add ultilities for SHD of CF
1380  *
1381  * ***********************************************************************/
Generated by   doxygen 1.8.5