StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
NonIdPurityCorrFctn.cxx
1 /***************************************************************************
2  *
3  * $Id: NonIdPurityCorrFctn.cxx,v 1.2 2003/01/31 19:21:09 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  * a simple Q-invariant correlation function
10  *
11  ***************************************************************************
12  *
13  * $Log: NonIdPurityCorrFctn.cxx,v $
14  * Revision 1.2 2003/01/31 19:21:09 magestro
15  * Cleared up simple compiler warnings on i386_linux24
16  *
17  * Revision 1.1 2002/12/12 17:02:49 kisiel
18  * Use KStar instead of 2*KStar for non-identical particles
19  *
20  * Revision 1.4 2000/01/25 17:34:45 laue
21  * I. In order to run the stand alone version of the StHbtMaker the following
22  * changes have been done:
23  * a) all ClassDefs and ClassImps have been put into #ifdef __ROOT__ statements
24  * b) unnecessary includes of StMaker.h have been removed
25  * c) the subdirectory StHbtMaker/doc/Make has been created including everything
26  * needed for the stand alone version
27  *
28  * II. To reduce the amount of compiler warning
29  * a) some variables have been type casted
30  * b) some destructors have been declared as virtual
31  *
32  * Revision 1.3 1999/07/29 02:47:09 lisa
33  * 1) add OpeningAngle correlation function 2) add StHbtMcEventReader 3) make histos in CorrFctns do errors correctly
34  *
35  * Revision 1.2 1999/07/06 22:33:20 lisa
36  * Adjusted all to work in pro and new - dev itself is broken
37  *
38  * Revision 1.1.1.1 1999/06/29 16:02:57 lisa
39  * Installation of StHbtMaker
40  *
41  **************************************************************************/
42 
43 #include "StHbtMaker/CorrFctn/NonIdPurityCorrFctn.h"
44 #include <cstdio>
45 
46 #ifdef __ROOT__
47 ClassImp(NonIdPurityCorrFctn)
48 #endif
49 
50 //____________________________
51 NonIdPurityCorrFctn::NonIdPurityCorrFctn(char* title, const int& nbins, const float& QinvLo, const float& QinvHi, int p1Type, int p2Type){
52  // set up numerator
53  char TitNumP[100] = "NumP";
54  strcat(TitNumP,title);
55  mNumP = new StHbt1DHisto(TitNumP,title,nbins,QinvLo,QinvHi);
56  // set up denominator
57  char TitDenP[100] = "DenP";
58  strcat(TitDenP,title);
59  mDenP = new StHbt1DHisto(TitDenP,title,nbins,QinvLo,QinvHi);
60  // set up ratio
61  char TitRatP[100] = "RatP";
62  strcat(TitRatP,title);
63  mRatP = new StHbt1DHisto(TitRatP,title,nbins,QinvLo,QinvHi);
64  // set up numerator
65  char TitNumN[100] = "NumN";
66  strcat(TitNumN,title);
67  mNumN = new StHbt1DHisto(TitNumN,title,nbins,QinvLo,QinvHi);
68  // set up denominator
69  char TitDenN[100] = "DenN";
70  strcat(TitDenN,title);
71  mDenN = new StHbt1DHisto(TitDenN,title,nbins,QinvLo,QinvHi);
72  // set up ratio
73  char TitRatN[100] = "RatN";
74  strcat(TitRatN,title);
75  mRatN = new StHbt1DHisto(TitRatN,title,nbins,QinvLo,QinvHi);
76  // set up ratio
77  char TitRat[100] = "Rat";
78  strcat(TitRat,title);
79  mRat = new StHbt1DHisto(TitRat,title,nbins,QinvLo,QinvHi);
80  // set up
81  char TitRatNOverP[100] = "RatNOverP";
82  strcat(TitRatNOverP,title);
83  mRatNOverP = new StHbt1DHisto(TitRatNOverP,title,nbins,QinvLo,QinvHi);
84 
85 // // set up numerator
86 // char TitNumPinvP[100] = "NumPinvP";
87 // strcat(TitNumPinvP,title);
88 // mNumPinvP = new StHbt1DHisto(TitNumPinvP,title,nbins,QinvLo,QinvHi);
89 // // set up denominator
90 // char TitDenPinvP[100] = "DenPinvP";
91 // strcat(TitDenPinvP,title);
92 // mDenPinvP = new StHbt1DHisto(TitDenPinvP,title,nbins,QinvLo,QinvHi);
93 // // set up numerator
94 // char TitNumPinvN[100] = "NumPinvN";
95 // strcat(TitNumPinvN,title);
96 // mNumPinvN = new StHbt1DHisto(TitNumPinvN,title,nbins,QinvLo,QinvHi);
97 // // set up denominator
98 // char TitDenPinvN[100] = "DenPinvN";
99 // strcat(TitDenPinvN,title);
100 // mDenPinvN = new StHbt1DHisto(TitDenPinvN,title,nbins,QinvLo,QinvHi);
101 // // set up ratio
102 // char TitRatPinv[100] = "RatPinv";
103 // strcat(TitRatPinv,title);
104 // mRatPinv = new StHbt1DHisto(TitRatPinv,title,nbins,QinvLo,QinvHi);
105 
106 // // set up ratio
107 // char TitRatPinvNormP[100] = "RatPinvNormP";
108 // strcat(TitRatPinvNormP,title);
109 // mRatPinvNormP = new StHbt1DHisto(TitRatPinvNormP,title,nbins,QinvLo,QinvHi);
110 // // set up ratio
111 // char TitRatPinvNormN[100] = "RatPinvNormN";
112 // strcat(TitRatPinvNormN,title);
113 // mRatPinvNormN = new StHbt1DHisto(TitRatPinvNormN,title,nbins,QinvLo,QinvHi);
114 // // set up ratio
115 // char TitRatPinvNorm[100] = "RatPinvNorm";
116 // strcat(TitRatPinvNorm,title);
117 // mRatPinvNorm = new StHbt1DHisto(TitRatPinvNorm,title,nbins,QinvLo,QinvHi);
118 // // set up ratio
119 // char TitRatPinvNormNOverP[100] = "RatPinvNormNOverP";
120 // strcat(TitRatPinvNormNOverP,title);
121 // mRatPinvNormNOverP = new StHbt1DHisto(TitRatPinvNormNOverP,
122 // title,nbins,QinvLo,QinvHi);
123  mp1Type = p1Type;
124  mp2Type = p2Type;
125  char TitPairPurity[100];
126  strcpy(TitPairPurity,"PairPurityOut");
127  strcat(TitPairPurity,title);
128  mPairPurityOut = new TProfile(TitPairPurity,title,2*nbins,-QinvHi,QinvHi);
129  strcpy(TitPairPurity,"PairPuritySide");
130  strcat(TitPairPurity,title);
131  mPairPuritySide = new TProfile(TitPairPurity,title,2*nbins,-QinvHi,QinvHi);
132  strcpy(TitPairPurity,"PairPurityLong");
133  strcat(TitPairPurity,title);
134  mPairPurityLong = new TProfile(TitPairPurity,title,2*nbins,-QinvHi,QinvHi);
135 
136  // to enable error bar calculation...
137  mNumP->Sumw2();
138  mDenP->Sumw2();
139  mRatP->Sumw2();
140  mNumN->Sumw2();
141  mDenN->Sumw2();
142  mRatN->Sumw2();
143  mRat->Sumw2();
144  mRatNOverP->Sumw2();
145 
146 // mNumPinvP->Sumw2();
147 // mDenPinvP->Sumw2();
148 // mNumPinvN->Sumw2();
149 // mDenPinvN->Sumw2();
150 // mRatPinv->Sumw2();
151 
152 // mRatPinvNormP->Sumw2();
153 // mRatPinvNormN->Sumw2();
154 // mRatPinvNorm->Sumw2();
155 // mRatPinvNormNOverP->Sumw2();
156 
157 // mQinvPt1 = new StHbt2DHisto("QinvPt1","QinvPt1",100,0.,1.,
158 // nbins,QinvLo,QinvHi);
159 // mQinvY1 = new StHbt2DHisto("QinvY1","QinvY1",100,-1.,1.,
160 // nbins,QinvLo,QinvHi);
161 
162 // mHKCVKSame = new StHbt2DHisto("HKCVKSame","HKCVKSame",
163 // 100,-1.,1.,nbins,QinvLo,QinvHi);
164 // mHKCVKDiff = new StHbt2DHisto("HKCVKDiff","HKCVKDiff",
165 // 100,-1.,1.,nbins,QinvLo,QinvHi);
166 }
167 
168 //____________________________
169 NonIdPurityCorrFctn::~NonIdPurityCorrFctn(){
170  delete mNumP;
171  delete mDenP;
172  delete mRatP;
173  delete mNumN;
174  delete mDenN;
175  delete mRatN;
176  delete mRat;
177  delete mRatNOverP;
178 
179 // delete mNumPinvP;
180 // delete mDenPinvP;
181 // delete mNumPinvN;
182 // delete mDenPinvN;
183 // delete mRatPinv;
184 
185 // delete mRatPinvNormP;
186 // delete mRatPinvNormN;
187 // delete mRatPinvNorm;
188 // delete mRatPinvNormNOverP;
189 
190 // delete mQinvPt1;
191 // delete mQinvY1;
192 
193 // delete mHKCVKSame;
194 // delete mHKCVKDiff;
195 
196  delete mPairPurityOut;
197  delete mPairPuritySide;
198  delete mPairPurityLong;
199 
200 }
201 //_________________________
202 void NonIdPurityCorrFctn::Finish(){
203  double tScale;
204  int tLastNormBin = mNumP->GetNbinsX();
205  int tFirstNormBin = tLastNormBin/2+1;
206  // Make cvk dependant correlation function
207  mRatP->Divide(mNumP,mDenP,1.0,1.0);
208  tScale = mRatP->Integral(tFirstNormBin,tLastNormBin);
209  tScale/= (tLastNormBin-tFirstNormBin+1);
210  mRatP->Scale(1./tScale);
211  mRatN->Divide(mNumN,mDenN,1.0,1.0);
212  tScale = mRatN->Integral(tFirstNormBin,tLastNormBin);
213  tScale/= (tLastNormBin-tFirstNormBin+1);
214  mRatN->Scale(1./tScale);
215  mRatNOverP->Divide(mRatN,mRatP,1.0,1.0);
216 
217 // mRatPinvNormP->Divide(mNumP,mNumPinvP,1.0,1.0);
218 // mRatPinvNormN->Divide(mNumN,mNumPinvN,1.0,1.0);
219 // mRatPinvNormNOverP->Divide(mRatPinvNormN,mRatPinvNormP);
220 
221  // Regular correlation function
222  TH1D tHNum(*mNumP);
223  tHNum.SetName("tHNum");
224  tHNum.Add(mNumN);
225  TH1D tHDen(*mDenP);
226  tHDen.SetName("tHDen");
227  tHDen.Add(mDenN);
228  mRat->Divide(&tHNum,&tHDen);
229  tScale = mRat->Integral(tFirstNormBin,tLastNormBin);
230  tScale/= (tLastNormBin-tFirstNormBin+1);
231  mRat->Scale(1./tScale);
232 
233  // Pinv correlation function
234 // TH1D tHNumPinv(*mNumPinvP);
235 // tHNumPinv.SetName("tHNumPinv");
236 // tHNumPinv.Add(mNumPinvN);
237 // TH1D tHDenPinv(*mDenPinvP);
238 // tHDenPinv.SetName("tHDenPinv");
239 // tHDenPinv.Add(mDenPinvN);
240 // mRatPinv->Divide(&tHNumPinv,&tHDenPinv);
241 // tScale = mRatPinv->Integral(tFirstNormBin,tLastNormBin);
242 // tScale/= (tLastNormBin-tFirstNormBin+1);
243 // mRatPinv->Scale(1./tScale);
244 
245  // Pinv normalised correlation function
246 // mRatPinvNorm->Divide(&tHNum,&tHNumPinv);
247 }
248 
249 void NonIdPurityCorrFctn::Write(){
250  mNumP->Write();
251  mDenP->Write();
252  mRatP->Write();
253  mNumN->Write();
254  mDenN->Write();
255  mRatN->Write();
256  mRat->Write();
257  mRatNOverP->Write();
258 
259 // mNumPinvP->Write();
260 // mDenPinvP->Write();
261 // mNumPinvN->Write();
262 // mDenPinvN->Write();
263 // mRatPinv->Write();
264 
265 // mRatPinvNormP->Write();
266 // mRatPinvNormN->Write();
267 // mRatPinvNorm->Write();
268 // mRatPinvNormNOverP->Write();
269 
270 // mHKCVKSame->Write();
271 // mHKCVKDiff->Write();
272  //mQinvPt1->Write();
273  //mQinvY1->Write();
274 
275  mPairPurityOut->Write();
276  mPairPuritySide->Write();
277  mPairPurityLong->Write();
278 }
279 
280 //____________________________
281 StHbtString NonIdPurityCorrFctn::Report(){
282  string stemp = "Qinv Correlation Function Report:\n";
283  char ctemp[100];
284  sprintf(ctemp,"Number of entries in numerator:\t%E\n",mNumP->GetEntries());
285  stemp += ctemp;
286  sprintf(ctemp,"Number of entries in denominator:\t%E\n",mDenP->GetEntries());
287  stemp += ctemp;
288  sprintf(ctemp,"Number of entries in ratio:\t%E\n",mRatP->GetEntries());
289  stemp += ctemp;
290  StHbtString returnThis = stemp;
291  return returnThis;
292 }
293 //____________________________
294 void NonIdPurityCorrFctn::AddRealPair(const StHbtPair* pair){
295  double tKStar = fabs(pair->KStar());
296  //double tPStar = 2*fabs(pair->KStarFlipped());
297  double tCVK = pair->CVK();
298  double pPurity = 0.0;
299  double tKOut = pair->dKOut();
300  double tKSide = pair->dKSide();
301  double tKLong = pair->dKLong();
302 
303  //mQinvPt1->Fill(pair->track1()->FourMomentum().perp(),tKStar,1);
304  //mQinvY1->Fill(pair->track1()->FourMomentum().rapidity(),tKStar,1);
305 
306  //if(pair->track1()->FourMomentum().pz()<0.){
307  // cout << pair->track1()->FourMomentum() << endl;
308  //}
309 // mHKCVKSame->Fill(tCVK,tKStar,1.);
310 
311  if(tCVK>0.){
312  mNumP->Fill(tKStar);
313  }
314  else{
315  mNumN->Fill(tKStar);
316  }
317 // if(pair->CVKFlipped()>0.){
318 // mNumPinvP->Fill(tPStar);
319 // }
320 // else{
321 // mNumPinvN->Fill(tPStar);
322 // }
323 
324  switch (mp1Type)
325  {
326  case 1:
327  pPurity = pair->track1()->GetPionPurity();
328  break;
329  case 2:
330  pPurity = pair->track1()->GetKaonPurity();
331  break;
332  case 3:
333  pPurity = pair->track1()->GetProtonPurity();
334  break;
335  case 4:
336  pPurity = pair->track1()->Track()->PidProbElectron();
337  break;
338  }
339 
340  // cout << "First particle purity: " << pPurity << endl;
341 
342  switch (mp2Type)
343  {
344  case 1:
345  pPurity *= pair->track2()->GetPionPurity();
346  break;
347  case 2:
348  pPurity *= pair->track2()->GetKaonPurity();
349  break;
350  case 3:
351  pPurity *= pair->track2()->GetProtonPurity();
352  break;
353  case 4:
354  pPurity *= pair->track2()->Track()->PidProbElectron();
355  break;
356  }
357  // cout << "Pair purity: " << pPurity << endl;
358 
359  mPairPurityOut->Fill(((tKOut>0)*2-1) *tKStar, pPurity);
360  mPairPuritySide->Fill(((tKSide>0)*2-1)*tKStar, pPurity);
361  mPairPurityLong->Fill(((tKLong>0)*2-1)*tKStar, pPurity);
362 
363 }
364 //____________________________
365 void NonIdPurityCorrFctn::AddMixedPair(const StHbtPair* pair){
366  double tKStar = 2*fabs(pair->KStar());
367  //double tPStar = 2*fabs(pair->KStarFlipped());
368  double tCVK = pair->CVK();
369 // mHKCVKDiff->Fill(tCVK,tKStar,1.);
370  if(tCVK>0.){
371  mDenP->Fill(tKStar);
372  }
373  else{
374  mDenN->Fill(tKStar);
375  }
376 // if(pair->CVKFlipped()>0.){
377 // mDenPinvP->Fill(tPStar);
378 // }
379 // else{
380 // mDenPinvN->Fill(tPStar);
381 // }
382 }
383 
384 
Generated by   doxygen 1.8.5