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