StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
BPLCMSFrame3DCorrFctn.cxx
1 /***************************************************************************
2  *
3  * $Id: BPLCMSFrame3DCorrFctn.cxx,v 1.7 2003/01/31 19:20:54 magestro Exp $
4  *
5  * Author: Mike Lisa, Ohio State, lisa@mps.ohio-state.edu
6  ***************************************************************************
7  *
8  * Description: part of STAR HBT Framework: StHbtMaker package
9  * 3D Bertsch-Pratt decomposition in the LCMS frame
10  *
11  ***************************************************************************
12  *
13  * $Log: BPLCMSFrame3DCorrFctn.cxx,v $
14  * Revision 1.7 2003/01/31 19:20:54 magestro
15  * Cleared up simple compiler warnings on i386_linux24
16  *
17  * Revision 1.6 2002/06/07 22:51:38 lisa
18  * Widely used BPLCMSFrame3DCorrFctn class now accumulates UNcorrected denominator and has a WriteOutHistos method
19  *
20  * Revision 1.5 2001/05/23 00:19:04 lisa
21  * Add in Smearing classes and methods needed for momentum resolution studies and correction
22  *
23  * Revision 1.4 2000/10/26 19:48:49 rcwells
24  * Added functionality for Coulomb correction of <qInv> in 3D correltions
25  *
26  * Revision 1.3 2000/09/14 18:36:53 lisa
27  * Added Qinv and ExitSep pair cuts and BPLCMSFrame3DCorrFctn_SIM CorrFctn
28  *
29  * Revision 1.2 2000/08/23 19:43:43 lisa
30  * added alternate normalization algorithm to 3d CorrFctns in case normal one fails
31  *
32  * Revision 1.1 2000/08/17 20:48:39 lisa
33  * Adding correlationfunction in LCMS frame
34  *
35  *
36  **************************************************************************/
37 
38 #include "StHbtMaker/CorrFctn/BPLCMSFrame3DCorrFctn.h"
39 //#include "StHbtMaker/Infrastructure/StHbtHisto.hh"
40 #include <cstdio>
41 
42 #ifdef __ROOT__
43 ClassImp(BPLCMSFrame3DCorrFctn)
44 #endif
45 
46 //____________________________
47 BPLCMSFrame3DCorrFctn::BPLCMSFrame3DCorrFctn(char* title, const int& nbins, const float& QLo, const float& QHi){
48 
49  // set some stuff...
50  mQinvNormLo = 0.15;
51  mQinvNormHi = 0.18;
52  mNumRealsNorm = 0;
53  mNumMixedNorm = 0;
54  mCorrection = 0; // pointer to Coulomb Correction object
55 
56  mPairCut = 0; // added Sept2000 - CorrFctn-specific PairCut
57 
58  mSmearPair = 0; // no resolution correction unless user sets SmearPair
59 
60  // set up numerator
61  char TitNum[100] = "Num";
62  strcat(TitNum,title);
63  mNumerator = new StHbt3DHisto(TitNum,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
64  // set up denominator
65  char TitDen[100] = "Den";
66  strcat(TitDen,title);
67  mDenominator = new StHbt3DHisto(TitDen,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
68  // set up uncorrected denominator
69  char TitDenUncoul[100] = "DenNoCoul";
70  strcat(TitDenUncoul,title);
71  mUncorrectedDenominator = new StHbt3DHisto(TitDenUncoul,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
72  // set up ratio
73  char TitRat[100] = "Rat";
74  strcat(TitRat,title);
75  mRatio = new StHbt3DHisto(TitRat,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
76  // set up ave qInv
77  char TitQinv[100] = "Qinv";
78  strcat(TitQinv,title);
79  mQinvHisto = new StHbt3DHisto(TitQinv,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
80 
81  // to enable error bar calculation...
82  mNumerator->Sumw2();
83  mDenominator->Sumw2();
84  mUncorrectedDenominator->Sumw2();
85  mRatio->Sumw2();
86 
87  // Following histos are for the momentum resolution correction
88  // they are filled only if a StHbtSmear object is plugged in
89  // here comes the "idea" numerator and denominator and ratio...
90  char TitNumID[100] = "IDNum";
91  strcat(TitNumID,title);
92  mIDNumHisto = new StHbt3DHisto(TitNumID,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
93  char TitDenID[100] = "IDDen";
94  strcat(TitDenID,title);
95  mIDDenHisto = new StHbt3DHisto(TitDenID,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
96  char TitRatID[100] = "IDRat";
97  strcat(TitRatID,title);
98  mIDRatHisto = new StHbt3DHisto(TitRatID,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
99 
100  mIDNumHisto->Sumw2();
101  mIDDenHisto->Sumw2();
102  mIDRatHisto->Sumw2();
103 
104  //
105  // here comes the "smeared" numerator and denominator...
106  char TitNumSM[100] = "SMNum";
107  strcat(TitNumSM,title);
108  mSMNumHisto = new StHbt3DHisto(TitNumSM,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
109  char TitDenSM[100] = "SMDen";
110  strcat(TitDenSM,title);
111  mSMDenHisto = new StHbt3DHisto(TitDenSM,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
112  char TitRatSM[100] = "SMRat";
113  strcat(TitRatSM,title);
114  mSMRatHisto = new StHbt3DHisto(TitRatSM,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
115  //
116  mSMNumHisto->Sumw2();
117  mSMDenHisto->Sumw2();
118  mSMRatHisto->Sumw2();
119  //
120  // here comes the correction factor (which is just ratio of ideal ratio to smeared ratio)
121  char TitCorrection[100] = "CorrectionFactor";
122  strcat(TitCorrection,title);
123  mCorrectionHisto = new StHbt3DHisto(TitCorrection,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
124  mCorrectionHisto->Sumw2();
125  // here comes the fully corrected correlation function
126  char TitCorrCF[100] = "CorrectedCF";
127  strcat(TitCorrCF,title);
128  mCorrCFHisto = new StHbt3DHisto(TitCorrCF,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
129  mCorrCFHisto->Sumw2();
130 
131  // user can (and should) override these defaults...
132  mLambda = 0.6;
133  mRout2 = 6.0*6.0;
134  mRside2 = 6.0*6.0;
135  mRlong2 = 7.0*7.0;
136 
137 }
138 
139 //____________________________
140 BPLCMSFrame3DCorrFctn::~BPLCMSFrame3DCorrFctn(){
141  delete mNumerator;
142  delete mDenominator;
143  delete mRatio;
144  delete mQinvHisto;
145  delete mIDNumHisto;
146  delete mIDDenHisto;
147  delete mIDRatHisto;
148  delete mSMNumHisto;
149  delete mSMDenHisto;
150  delete mSMRatHisto;
151  delete mCorrectionHisto;
152  delete mCorrCFHisto;
153 }
154 
155 //_________________________
156 void BPLCMSFrame3DCorrFctn::WriteOutHistos(){
157 
158  mNumerator->Write();
159  mDenominator->Write();
160  mUncorrectedDenominator->Write();
161  mRatio->Write();
162  mQinvHisto->Write();
163 
164  if (mSmearPair){
165  mIDNumHisto->Write();
166  mIDDenHisto->Write();
167  mIDRatHisto->Write();
168  //
169  mSMNumHisto->Write();
170  mSMDenHisto->Write();
171  mSMRatHisto->Write();
172  //
173  mCorrectionHisto->Write();
174  mCorrCFHisto->Write();
175  }
176 }
177 
178 //_________________________
179 void BPLCMSFrame3DCorrFctn::Finish(){
180  // here is where we should normalize, fit, etc...
181  double NumFact,DenFact;
182  if ((mNumRealsNorm !=0) && (mNumMixedNorm !=0)){
183  NumFact = double(mNumRealsNorm);
184  DenFact = double(mNumMixedNorm);
185  }
186  // can happen that the mNumRealsNorm and mNumMixedNorm = 0 if you do non-standard
187  // things like making a new CorrFctn and just setting the Numerator and Denominator
188  // from OTHER CorrFctns which you read in (like when doing parallel processing)
189  else{
190  cout << "Warning! - no normalization constants defined - I do the best I can..." << endl;
191  int nbins = mNumerator->GetNbinsX();
192  int half_way = nbins/2;
193  NumFact = mNumerator->Integral(half_way,nbins,half_way,nbins,half_way,nbins);
194  DenFact = mDenominator->Integral(half_way,nbins,half_way,nbins,half_way,nbins);
195  }
196 
197  mRatio->Divide(mNumerator,mDenominator,DenFact,NumFact);
198  mQinvHisto->Divide(mUncorrectedDenominator);
199 
200  // now do all the resolution correction stuff..
201  if (mSmearPair){ // but only do it if we have been working with a SmearPair
202  mIDRatHisto->Divide(mIDNumHisto,mIDDenHisto);
203  mSMRatHisto->Divide(mSMNumHisto,mSMDenHisto);
204  mCorrectionHisto->Divide(mIDRatHisto,mSMRatHisto);
205  mCorrCFHisto->Multiply(mRatio,mCorrectionHisto);
206  }
207 
208 }
209 
210 //____________________________
211 StHbtString BPLCMSFrame3DCorrFctn::Report(){
212  string stemp = "LCMS Frame Bertsch-Pratt 3D Correlation Function Report:\n";
213  char ctemp[100];
214  sprintf(ctemp,"Number of entries in numerator:\t%E\n",mNumerator->GetEntries());
215  stemp += ctemp;
216  sprintf(ctemp,"Number of entries in denominator:\t%E\n",mDenominator->GetEntries());
217  stemp += ctemp;
218  sprintf(ctemp,"Number of entries in ratio:\t%E\n",mRatio->GetEntries());
219  stemp += ctemp;
220  sprintf(ctemp,"Normalization region in Qinv was:\t%E\t%E\n",mQinvNormLo,mQinvNormHi);
221  stemp += ctemp;
222  sprintf(ctemp,"Number of pairs in Normalization region was:\n");
223  stemp += ctemp;
224  sprintf(ctemp,"In numerator:\t%lu\t In denominator:\t%lu\n",mNumRealsNorm,mNumMixedNorm);
225  stemp += ctemp;
226  if (mCorrection)
227  {
228  float radius = mCorrection->GetRadius();
229  sprintf(ctemp,"Coulomb correction used radius of\t%E\n",radius);
230  }
231  else
232  {
233  sprintf(ctemp,"No Coulomb Correction applied to this CorrFctn\n");
234  }
235  stemp += ctemp;
236 
237  if (mPairCut){
238  sprintf(ctemp,"Here is the PairCut specific to this CorrFctn\n");
239  stemp += ctemp;
240  stemp += mPairCut->Report();
241  }
242  else{
243  sprintf(ctemp,"No PairCut specific to this CorrFctn\n");
244  stemp += ctemp;
245  }
246 
247  //
248  StHbtString returnThis = stemp;
249  return returnThis;
250 }
251 //____________________________
252 void BPLCMSFrame3DCorrFctn::AddRealPair(const StHbtPair* pair){
253 
254  if (mPairCut){
255  if (!(mPairCut->Pass(pair))) return;
256  }
257 
258  double Qinv = fabs(pair->qInv()); // note - qInv() will be negative for identical pairs...
259  if ((Qinv < mQinvNormHi) && (Qinv > mQinvNormLo)) mNumRealsNorm++;
260  double qOut = fabs(pair->qOutCMS());
261  double qSide = fabs(pair->qSideCMS());
262  double qLong = fabs(pair->qLongCMS());
263 
264  mNumerator->Fill(qOut,qSide,qLong);
265 }
266 //____________________________
267 void BPLCMSFrame3DCorrFctn::AddMixedPair(const StHbtPair* pair){
268 
269  if (mPairCut){
270  if (!(mPairCut->Pass(pair))) return;
271  }
272 
273  double CoulombWeight = (mCorrection ? mCorrection->CoulombCorrect(pair) : 1.0);
274 
275  double Qinv = fabs(pair->qInv()); // note - qInv() will be negative for identical pairs...
276  if ((Qinv < mQinvNormHi) && (Qinv > mQinvNormLo)) mNumMixedNorm++;
277  double qOut = fabs(pair->qOutCMS());
278  double qSide = fabs(pair->qSideCMS());
279  double qLong = fabs(pair->qLongCMS());
280 
281  mDenominator->Fill(qOut,qSide,qLong,CoulombWeight);
282  mUncorrectedDenominator->Fill(qOut,qSide,qLong,1.0);
283  mQinvHisto->Fill(qOut,qSide,qLong,Qinv);
284 
285  // now for the momentum resolution stuff...
286  if (mSmearPair){
287  double CorrWeight = 1.0 +
288  mLambda*exp((-qOut*qOut*mRout2 -qSide*qSide*mRside2 -qLong*qLong*mRlong2)/0.038936366329);
289  CorrWeight *= CoulombWeight; // impt.
290 
291  mIDNumHisto->Fill(qOut,qSide,qLong,CorrWeight);
292  mIDDenHisto->Fill(qOut,qSide,qLong,CoulombWeight);
293 
294  mSmearPair->SetUnsmearedPair(pair);
295  double qOut_prime = fabs(mSmearPair->SmearedPair().qOutCMS());
296  double qSide_prime = fabs(mSmearPair->SmearedPair().qSideCMS());
297  double qLong_prime = fabs(mSmearPair->SmearedPair().qLongCMS());
298 
299  mSMNumHisto->Fill(qOut_prime,qSide_prime,qLong_prime,CorrWeight);
300 
301  double SmearedCoulombWeight = ( mCorrection ?
302  mCorrection->CoulombCorrect(&(mSmearPair->SmearedPair())) :
303  1.0);
304 
305  mSMDenHisto->Fill(qOut_prime,qSide_prime,qLong_prime,SmearedCoulombWeight);
306  }
307 }
308 
Generated by   doxygen 1.8.5