StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
YKPLongitudinal.cxx
1 /***************************************************************************
2  *
3  * Author: Dominik Flierl, flierl@bnl.gov
4  ***************************************************************************
5  *
6  * Description:
7  * YKP parametrization in longitudinal direction in kt-Y-Bins
8  *
9  *
10  **************************************************************************/
11 
12 
13 #include "StHbtMaker/CorrFctn/YKPLongitudinal.h"
14 #include <cstdio>
15 
16 #ifdef __ROOT__
17 ClassImp(YKPLongitudinal)
18 #endif
19 
20 //____________________________
21 YKPLongitudinal::YKPLongitudinal( TString ctype, TString frame ,
22  const int& nbins, const float& QLo, const float& QHi,
23  const int& nKtbins, const double& ktMin, const double& ktMax,
24  const int& nYbins, const double& YMin, const double& YMax,
25  const int& nbinsQINV, const float& QLoQINV , const float& QHiQINV )
26 {
27  // set correlation function type
28  if ( !ctype.CompareTo("YKP") || !ctype.CompareTo("BP") )
29  {
30  mCtype = ctype ;
31  cout << "correlation function set to :" << ctype << endl ;
32  }
33  else
34  {
35  cout << "Error: correlationfunction didn't fit BP or YKP -> YKP used !\n" ;
36  mCtype = "YKP" ;
37  }
38 
39  // set frame
40  if ( !frame.CompareTo("LCMS") || !frame.CompareTo("CMS") || !frame.CompareTo("PRF"))
41  {
42  mFrame = frame ;
43  cout << "frame set to :" << frame << endl ;
44  }
45  else
46  {
47  cout << "Error: frame choice didn't fit LCMS,CMS or PRF -> LCMS used !\n" ;
48  mFrame = "LCMS" ;
49  }
50 
51  // get intializers
52  mNumberKtBins = nKtbins ;
53  mNumberYBins = nYbins ;
54 
55  // note : bin = binKt + (binY-1) * mNumberKtBins
56  mNumberBins = mNumberKtBins * mNumberYBins ;
57 
58  // pointer to Coulomb Correction object must be set from macro
59  mCorrection = 0 ;
60 
61  // set up 3d histos
62  mNumerator = new StHbt3DHisto[mNumberBins] ;
63  mDenominator = new StHbt3DHisto[mNumberBins] ;
64  mRatio = new StHbt3DHisto[mNumberBins] ;
65  mQinvNumerator = new StHbt1DHisto[mNumberBins] ;
66  mQinvDenominator = new StHbt1DHisto[mNumberBins] ;
67  mQinvRatio = new StHbt1DHisto[mNumberBins] ;
68 
69  char TitNum[100] ;
70  char TitDen[100] ;
71  char TitRat[100] ;
72  char TitQinvNum[100] ;
73  char TitQinvDen[100] ;
74  char TitQinvRatio[100] ;
75 
76  char xAxisTitle[100] ;
77  char yAxisTitle[100] ;
78  char zAxisTitle[100] ;
79  if (mCtype == "YKP")
80  {
81  sprintf(xAxisTitle,"qLong") ;
82  sprintf(yAxisTitle,"qPerp") ;
83  sprintf(zAxisTitle,"qNull") ;
84  }
85  else if (mCtype == "BP")
86  {
87  sprintf(xAxisTitle,"qSide") ;
88  sprintf(yAxisTitle,"qOut") ;
89  sprintf(zAxisTitle,"qLong") ;
90  }
91 
92  for(int ktindex = 0 ; ktindex < mNumberKtBins ; ktindex++)
93  {
94  for(int yindex = 0 ; yindex < mNumberYBins ; yindex++)
95  {
96  // numerator
97  sprintf(TitNum,"NumKt%dY%d",ktindex+1,yindex+1) ;
98  mNumerator[ktindex+yindex*mNumberKtBins].SetName(TitNum);
99  mNumerator[ktindex+yindex*mNumberKtBins].SetTitle(TitNum);
100  mNumerator[ktindex+yindex*mNumberKtBins].SetBins(nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi) ;
101  mNumerator[ktindex+yindex*mNumberKtBins].SetXTitle(xAxisTitle);
102  mNumerator[ktindex+yindex*mNumberKtBins].SetYTitle(yAxisTitle);
103  mNumerator[ktindex+yindex*mNumberKtBins].SetZTitle(zAxisTitle);
104  // denominator
105  sprintf(TitDen,"DenKt%dY%d",ktindex+1,yindex+1) ;
106  mDenominator[ktindex+yindex*mNumberKtBins].SetName(TitDen);
107  mDenominator[ktindex+yindex*mNumberKtBins].SetTitle(TitDen);
108  mDenominator[ktindex+yindex*mNumberKtBins].SetBins(nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi) ;
109  mDenominator[ktindex+yindex*mNumberKtBins].SetXTitle(xAxisTitle);
110  mDenominator[ktindex+yindex*mNumberKtBins].SetYTitle(yAxisTitle);
111  mDenominator[ktindex+yindex*mNumberKtBins].SetZTitle(zAxisTitle);
112  // ratio
113  sprintf(TitRat,"RatKt%dY%d",ktindex+1,yindex+1) ;
114  mRatio[ktindex+yindex*mNumberKtBins].SetName(TitRat);
115  mRatio[ktindex+yindex*mNumberKtBins].SetTitle(TitRat);
116  mRatio[ktindex+yindex*mNumberKtBins].SetBins(nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi) ;
117  mRatio[ktindex+yindex*mNumberKtBins].SetXTitle(xAxisTitle);
118  mRatio[ktindex+yindex*mNumberKtBins].SetYTitle(yAxisTitle);
119  mRatio[ktindex+yindex*mNumberKtBins].SetZTitle(zAxisTitle);
120  // qinv numerator
121  sprintf(TitQinvNum,"QinvNumeratorKt%dY%d",ktindex+1,yindex+1) ;
122  mQinvNumerator[ktindex+yindex*mNumberKtBins].SetName(TitQinvNum);
123  mQinvNumerator[ktindex+yindex*mNumberKtBins].SetTitle(TitQinvNum);
124  mQinvNumerator[ktindex+yindex*mNumberKtBins].SetBins(nbinsQINV,QLoQINV,QHiQINV) ;
125  mQinvNumerator[ktindex+yindex*mNumberKtBins].SetXTitle("qInv");
126  // qinv denominator
127  sprintf(TitQinvDen,"QinvDenominatorKt%dY%d",ktindex+1,yindex+1) ;
128  mQinvDenominator[ktindex+yindex*mNumberKtBins].SetName(TitQinvDen);
129  mQinvDenominator[ktindex+yindex*mNumberKtBins].SetTitle(TitQinvDen);
130  mQinvDenominator[ktindex+yindex*mNumberKtBins].SetBins(nbinsQINV,QLoQINV,QHiQINV) ;
131  mQinvDenominator[ktindex+yindex*mNumberKtBins].SetXTitle("qInv");
132  // qinv ratio
133  sprintf(TitQinvRatio,"QinvRatioKt%dY%d",ktindex+1,yindex+1) ;
134  mQinvRatio[ktindex+yindex*mNumberKtBins].SetName(TitQinvRatio);
135  mQinvRatio[ktindex+yindex*mNumberKtBins].SetTitle(TitQinvRatio);
136  mQinvRatio[ktindex+yindex*mNumberKtBins].SetBins(nbinsQINV,QLoQINV,QHiQINV) ;
137  mQinvRatio[ktindex+yindex*mNumberKtBins].SetXTitle("qInv") ;
138  // to enable error bar calculation...
139  mNumerator[ktindex+yindex*mNumberKtBins].Sumw2() ;
140  mDenominator[ktindex+yindex*mNumberKtBins].Sumw2() ;
141  mRatio[ktindex+yindex*mNumberKtBins].Sumw2() ;
142  mQinvNumerator[ktindex+yindex*mNumberKtBins].Sumw2() ;
143  mQinvDenominator[ktindex+yindex*mNumberKtBins].Sumw2() ;
144  mQinvRatio[ktindex+yindex*mNumberKtBins].Sumw2() ;
145  }
146  }
147 
148 
150  // set up y pt bins boundaries
152  // kt
153  mktBinsMin = new double[mNumberKtBins] ;
154  mktBinsMax = new double[mNumberKtBins] ;
155  double ktstep = (double)(ktMax-ktMin)/(double) mNumberKtBins ;
156  for(int ktindex = 0 ; ktindex < mNumberKtBins ; ktindex++ )
157  {
158  mktBinsMin[ktindex] = ((double) ktindex) * ktstep + ktMin ;
159  mktBinsMax[ktindex] = ((double) (ktindex+1)) * ktstep + ktMin ;
160  }
161  // rap
162  mYBinsMin = new double[mNumberKtBins] ;
163  mYBinsMax = new double[mNumberYBins] ;
164  double ystep = (double)(YMax-YMin)/(double) mNumberYBins ;
165  for(int yindex = 0 ; yindex < mNumberYBins ; yindex++ )
166  {
167  mYBinsMin[yindex] = ((double) yindex) * ystep + YMin ;
168  mYBinsMax[yindex] = ((double) (yindex+1)) * ystep + YMin ;
169  }
170 
172  // set up normalization
174  // set range normalization
175  mQinvNormLo = 0.15 ;
176  mQinvNormHi = 0.18 ;
177  // set up array with normalization values
178  mNumRealsNorm = new unsigned long int [mNumberBins] ;
179  mNumMixedNorm = new unsigned long int [mNumberBins] ;
180  for( int index = 0 ; index < mNumberBins ; index++ )
181  {
182  mNumRealsNorm[index] = mNumMixedNorm[index] = 0 ;
183  }
184 }
185 //____________________________
186 YKPLongitudinal::~YKPLongitudinal()
187 {
188  delete [] mNumRealsNorm ;
189  delete [] mNumMixedNorm ;
190  delete [] mNumerator ;
191  delete [] mDenominator ;
192  delete [] mRatio ;
193  delete [] mQinvNumerator ;
194  delete [] mQinvDenominator ;
195  delete [] mQinvRatio ;
196 }
197 //_________________________
198 void YKPLongitudinal::Finish()
199 {
200  // here is where we can normalize, fit, etc...
201  double NumFact,DenFact;
202  for(int ktindex = 0 ; ktindex < mNumberKtBins ; ktindex++)
203  {
204  for(int yindex = 0 ; yindex < mNumberYBins ; yindex++)
205  {
206  // get the normalisation
207  if ((mNumRealsNorm[ktindex+yindex*mNumberKtBins] !=0) && (mNumMixedNorm[ktindex+yindex*mNumberKtBins] !=0))
208  {
209  NumFact = double(mNumRealsNorm[ktindex+yindex*mNumberKtBins]) ;
210  DenFact = double(mNumMixedNorm[ktindex+yindex*mNumberKtBins]) ;
211  cout << " Normalizing with Num/Denom = norm : " << NumFact << "\t" << DenFact << "\t" << NumFact/DenFact
212  << " for kt bin : " << ktindex << " y bin : " << yindex << endl ;
213  }
214  // error no normalization filled
215  else
216  {
217  cout << "Warning! - no normalization possible ..." << endl;
218  NumFact = 1.0 ;
219  DenFact = 10.0 ;
220  }
221  // normalize ratio
222  mRatio[ktindex+yindex*mNumberKtBins].Divide(&mNumerator[ktindex+yindex*mNumberKtBins],
223  &mDenominator[ktindex+yindex*mNumberKtBins],DenFact,NumFact) ;
224 
225  // produce qinv ratio
226  mQinvRatio[ktindex+yindex*mNumberKtBins].Divide(&mQinvNumerator[ktindex+yindex*mNumberKtBins],
227  &mQinvDenominator[ktindex+yindex*mNumberKtBins]) ;
228  }
229  }
230 }
231 //____________________________
232 StHbtString YKPLongitudinal::Report()
233 {
234  // report called by manager per standard
235  // get the overall values
236  int mNumeratorEntriesAll = 0 ;
237  int mDenominatorEntriesAll = 0 ;
238  int mRatioEntriesAll = 0 ;
239  int mNumRealsNormAll = 0 ;
240  int mNumMixedNormAll = 0 ;
241 
242  for(int ktindex = 0 ; ktindex < mNumberKtBins ; ktindex++)
243  {
244  for(int yindex = 0 ; yindex < mNumberYBins ; yindex++)
245  {
246  mNumeratorEntriesAll += (int)mNumerator[ktindex+yindex*mNumberKtBins].GetEntries() ;
247  mDenominatorEntriesAll += (int)mDenominator[ktindex+yindex*mNumberKtBins].GetEntries() ;
248  mRatioEntriesAll += (int)mRatio[ktindex+yindex*mNumberKtBins].GetEntries() ;
249  mNumRealsNormAll += mNumRealsNorm[ktindex+yindex*mNumberKtBins] ;
250  mNumMixedNormAll += mNumMixedNorm[ktindex+yindex*mNumberKtBins] ;
251  }
252  }
253  // put the report together
254  string stemp = "YKP Function Report:\n";
255  char ctemp[100];
256  sprintf(ctemp,"Number of entries in numerator in all %d histos:\t%d\n",mNumberBins,mNumeratorEntriesAll);
257  stemp += ctemp;
258  sprintf(ctemp,"Number of entries in denominator in all %d histos:\t%d\n",mNumberBins,mDenominatorEntriesAll);
259  stemp += ctemp;
260  sprintf(ctemp,"Number of entries in ratio in all %d histos:\t%d\n",mNumberBins,mRatioEntriesAll);
261  stemp += ctemp;
262  sprintf(ctemp,"Normalization region in Qinv was:\t%E\t%E\n",mQinvNormLo,mQinvNormHi);
263  stemp += ctemp;
264  sprintf(ctemp,"Number of pairs in Normalization region was:\n");
265  stemp += ctemp;
266  sprintf(ctemp,"In numerator:\t%u\t In denominator:\t%u\n",mNumRealsNormAll,mNumMixedNormAll);
267  stemp += ctemp ;
268  if (mCorrection)
269  {
270  float radius = mCorrection->GetRadius();
271  sprintf(ctemp,"Coulomb correction used radius of\t%E\n",radius) ;
272  }
273  else
274  {
275  sprintf(ctemp,"No Coulomb Correction applied to this CorrFctn\n") ;
276  }
277  stemp += ctemp ;
278 
279  // return it
280  StHbtString returnThis = stemp ;
281  return returnThis ;
282 }
283 //____________________________
284 void YKPLongitudinal::AddRealPair(const StHbtPair* pair)
285 {
286  // get the pair momenta
287  double mKt = pair->kT() ;
288  double mRap = pair->rap() ;
289  double mQinv = fabs(pair->qInv()) ;
290  // sort it into the approriate histo
291  for(int ktindex = 0 ; ktindex < mNumberKtBins ; ktindex++)
292  {
293  if( mKt >= mktBinsMin[ktindex] && mKt < mktBinsMax[ktindex] )
294  {
295  for(int yindex = 0 ; yindex < mNumberYBins ; yindex++)
296  {
297  if( mRap >= mYBinsMin[yindex] && mRap < mYBinsMax[yindex] )
298  {
299  // sum up qinv pairs to do later the normalisation
300  if ((mQinv < mQinvNormHi) && (mQinv > mQinvNormLo)) mNumRealsNorm[ktindex+yindex*mNumberKtBins]++;
301  // fill 1d qinv histo
302  mQinvNumerator[ktindex+yindex*mNumberKtBins].Fill(mQinv) ;
303  // fill 3d : choose correlationfunction type and frame
304  if (mCtype == "YKP")
305  {
306  double qlong, qperp, q0 ;
307  if(mFrame =="LCMS") { pair->qYKPLCMS(qlong,qperp,q0) ;}
308  else if(mFrame =="CMS") { pair->qYKPCMS(qlong,qperp,q0) ;}
309  else if(mFrame =="PRF") { pair->qYKPPF(qlong,qperp,q0) ;}
310  // fill 3d histo
311  mNumerator[ktindex+yindex*mNumberKtBins].Fill(qlong,qperp,q0) ;
312  }
313  else if (mCtype == "BP")
314  {
315  double qside, qout, qlong;
316  if(mFrame =="LCMS") { qside = pair->qSideCMS(); qout = pair->qOutCMS(); qlong = pair->qLongCMS(); }
317  else if(mFrame =="CMS") { qside = pair->qSideBf(); qout = pair->qOutBf(); qlong = pair->qLongBf(); }
318  else if(mFrame =="PRF") { qside = pair->qSidePf(); qout = pair->qOutPf(); qlong = pair->qLongPf(); }
319  // fill 3d histo
320  mNumerator[ktindex+yindex*mNumberKtBins].Fill(fabs(qside),fabs(qout),fabs(qlong)) ;
321  }
322  // end yindex loop
323  break ;
324  }
325  }// for yindex
326  // end kt loop
327  break ;
328  } // if kt
329  } // for kt
330 }
331 //____________________________
332 void YKPLongitudinal::AddMixedPair(const StHbtPair* pair)
333 {
334  // get the pair momenta
335  double mKt = pair->kT() ;
336  double mRap = pair->rap() ;
337  double mQinv = fabs(pair->qInv());
338  // caculate coulomb weight
339  double weight=1.0;
340  if (mCorrection)
341  {
342  weight = mCorrection->CoulombCorrect(pair);
343  }
344 
345  // sort it into the approriate histo
346  for(int ktindex = 0 ; ktindex < mNumberKtBins ; ktindex++)
347  {
348  if( mKt >= mktBinsMin[ktindex] && mKt < mktBinsMax[ktindex] )
349  {
350  for(int yindex = 0 ; yindex < mNumberYBins ; yindex++)
351  {
352  if( mRap >= mYBinsMin[yindex] && mRap < mYBinsMax[yindex] )
353  {
354  // sum up qinv pairs to do later the normalisation
355  if ((mQinv < mQinvNormHi) && (mQinv > mQinvNormLo)) mNumMixedNorm[ktindex+yindex*mNumberKtBins]++;
356  // fill 1d qinv histo
357  mQinvDenominator[ktindex+yindex*mNumberKtBins].Fill(mQinv,weight) ;
358 
359  // fill 3d : choose correlationfunction type and frame
360  if (mCtype == "YKP")
361  {
362  double qlong, qperp, q0 ;
363  if(mFrame =="LCMS") { pair->qYKPLCMS(qlong,qperp,q0) ;}
364  else if(mFrame =="CMS") { pair->qYKPCMS(qlong,qperp,q0) ;}
365  else if(mFrame =="PRF") { pair->qYKPPF(qlong,qperp,q0) ;}
366  // fill 3d histo
367  mDenominator[ktindex+yindex*mNumberKtBins].Fill(qlong,qperp,q0,weight) ;
368  }
369  else if (mCtype == "BP")
370  {
371  double qside, qout, qlong;
372  if(mFrame =="LCMS") { qside = pair->qSideCMS(); qout = pair->qOutCMS(); qlong = pair->qLongCMS(); }
373  else if(mFrame =="CMS") { qside = pair->qSideBf(); qout = pair->qOutBf(); qlong = pair->qLongBf(); }
374  else if(mFrame =="PRF") { qside = pair->qSidePf(); qout = pair->qOutPf(); qlong = pair->qLongPf(); }
375  // fill 3d histo
376  mDenominator[ktindex+yindex*mNumberKtBins].Fill(fabs(qside),fabs(qout),fabs(qlong),weight) ;
377  }
378  // end yindex loop
379  break ;
380  }
381  }// for yindex
382  // end kt loop
383  break ;
384  } // if kt
385  } // for kt
386 }
Generated by   doxygen 1.8.5