Fix AliFemtoModelFreezeOutGenerator undefined references
[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$
0215f606 14 * Revision 1.1.1.1 2007/04/25 15:38:41 panos
15 * Importing the HBT code dir
16 *
67427ff7 17 * Revision 1.1.1.1 2007/03/07 10:14:49 mchojnacki
18 * First version on CVS
19 *
20 * Revision 1.7 2003/01/31 19:20:54 magestro
21 * Cleared up simple compiler warnings on i386_linux24
22 *
23 * Revision 1.6 2002/06/07 22:51:38 lisa
24 * Widely used AliFemtoBPLCMS3DCorrFctn class now accumulates UNcorrected denominator and has a WriteOutHistos method
25 *
26 * Revision 1.5 2001/05/23 00:19:04 lisa
27 * Add in Smearing classes and methods needed for momentum resolution studies and correction
28 *
29 * Revision 1.4 2000/10/26 19:48:49 rcwells
30 * Added functionality for Coulomb correction of <qInv> in 3D correltions
31 *
32 * Revision 1.3 2000/09/14 18:36:53 lisa
33 * Added Qinv and ExitSep pair cuts and AliFemtoBPLCMS3DCorrFctn_SIM CorrFctn
34 *
35 * Revision 1.2 2000/08/23 19:43:43 lisa
36 * added alternate normalization algorithm to 3d CorrFctns in case normal one fails
37 *
38 * Revision 1.1 2000/08/17 20:48:39 lisa
39 * Adding correlationfunction in LCMS frame
40 *
41 *
42 **************************************************************************/
43
44#include "CorrFctn/AliFemtoBPLCMS3DCorrFctn.h"
45//#include "Infrastructure/AliFemtoHisto.h"
46#include <cstdio>
47
48#ifdef __ROOT__
49ClassImp(AliFemtoBPLCMS3DCorrFctn)
50#endif
51
52//____________________________
0215f606 53AliFemtoBPLCMS3DCorrFctn::AliFemtoBPLCMS3DCorrFctn(char* title, const int& nbins, const float& QLo, const float& QHi)
54 :
55 fIDNumHisto(0),
56 fIDDenHisto(0),
57 fIDRatHisto(0),
58 fSMNumHisto(0),
59 fSMDenHisto(0),
60 fSMRatHisto(0),
61 fCorrectionHisto(0),
62 fCorrCFHisto(0),
63 fNumerator(0),
64 fDenominator(0),
65 fRatio(0),
66 fQinvHisto(0),
67 fLambda(0),
68 fRout2(0),
69 fRside2(0),
70 fRlong2(0),
71 fPairCut(0),
72 fQinvNormLo(0),
73 fQinvNormHi(0),
74 fNumRealsNorm(0),
75 fNumMixedNorm(0)
76{
67427ff7 77
78 // set some stuff...
79 fQinvNormLo = 0.15;
80 fQinvNormHi = 0.18;
81 fNumRealsNorm = 0;
82 fNumMixedNorm = 0;
83 // fCorrection = 0; // pointer to Coulomb Correction object
84
85 fPairCut = 0; // added Sept2000 - CorrFctn-specific PairCut
86
87 // fSmearPair = 0; // no resolution correction unless user sets SmearPair
88
89 // set up numerator
90 char TitNum[100] = "Num";
91 strcat(TitNum,title);
92 fNumerator = new TH3D(TitNum,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
93 // set up denominator
94 char TitDen[100] = "Den";
95 strcat(TitDen,title);
96 fDenominator = new TH3D(TitDen,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
97 // set up uncorrected denominator
98 char TitDenUncoul[100] = "DenNoCoul";
99 strcat(TitDenUncoul,title);
100 // fUncorrectedDenominator = new TH3D(TitDenUncoul,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
101 // set up ratio
102 char TitRat[100] = "Rat";
103 strcat(TitRat,title);
104 fRatio = new TH3D(TitRat,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
105 // set up ave qInv
106 char TitQinv[100] = "Qinv";
107 strcat(TitQinv,title);
108 fQinvHisto = new TH3D(TitQinv,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
109
110 // to enable error bar calculation...
111 fNumerator->Sumw2();
112 fDenominator->Sumw2();
113 // fUncorrectedDenominator->Sumw2();
114 fRatio->Sumw2();
115
116 // Following histos are for the momentum resolution correction
117 // they are filled only if a AliFemtoSmear object is plugged in
118 // here comes the "idea" numerator and denominator and ratio...
119 char TitNumID[100] = "IDNum";
120 strcat(TitNumID,title);
121 fIDNumHisto = new TH3D(TitNumID,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
122 char TitDenID[100] = "IDDen";
123 strcat(TitDenID,title);
124 fIDDenHisto = new TH3D(TitDenID,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
125 char TitRatID[100] = "IDRat";
126 strcat(TitRatID,title);
127 fIDRatHisto = new TH3D(TitRatID,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
128
129 fIDNumHisto->Sumw2();
130 fIDDenHisto->Sumw2();
131 fIDRatHisto->Sumw2();
132
133 //
134 // here comes the "smeared" numerator and denominator...
135 char TitNumSM[100] = "SMNum";
136 strcat(TitNumSM,title);
137 fSMNumHisto = new TH3D(TitNumSM,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
138 char TitDenSM[100] = "SMDen";
139 strcat(TitDenSM,title);
140 fSMDenHisto = new TH3D(TitDenSM,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
141 char TitRatSM[100] = "SMRat";
142 strcat(TitRatSM,title);
143 fSMRatHisto = new TH3D(TitRatSM,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
144 //
145 fSMNumHisto->Sumw2();
146 fSMDenHisto->Sumw2();
147 fSMRatHisto->Sumw2();
148 //
149 // here comes the correction factor (which is just ratio of ideal ratio to smeared ratio)
150 char TitCorrection[100] = "CorrectionFactor";
151 strcat(TitCorrection,title);
152 fCorrectionHisto = new TH3D(TitCorrection,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
153 fCorrectionHisto->Sumw2();
154 // here comes the fully corrected correlation function
155 char TitCorrCF[100] = "CorrectedCF";
156 strcat(TitCorrCF,title);
157 fCorrCFHisto = new TH3D(TitCorrCF,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
158 fCorrCFHisto->Sumw2();
159
160 // user can (and should) override these defaults...
161 fLambda = 0.6;
162 fRout2 = 6.0*6.0;
163 fRside2 = 6.0*6.0;
164 fRlong2 = 7.0*7.0;
165
166}
167
0215f606 168AliFemtoBPLCMS3DCorrFctn::AliFemtoBPLCMS3DCorrFctn(const AliFemtoBPLCMS3DCorrFctn& aCorrFctn) :
169 fIDNumHisto(0),
170 fIDDenHisto(0),
171 fIDRatHisto(0),
172 fSMNumHisto(0),
173 fSMDenHisto(0),
174 fSMRatHisto(0),
175 fCorrectionHisto(0),
176 fCorrCFHisto(0),
177 fNumerator(0),
178 fDenominator(0),
179 fRatio(0),
180 fQinvHisto(0),
181 fLambda(0),
182 fRout2(0),
183 fRside2(0),
184 fRlong2(0),
185 fPairCut(0),
186 fQinvNormLo(0),
187 fQinvNormHi(0),
188 fNumRealsNorm(0),
189 fNumMixedNorm(0)
190{
191 fIDNumHisto = aCorrFctn.fIDNumHisto;
192 fIDDenHisto = aCorrFctn.fIDDenHisto;
193 fIDRatHisto = aCorrFctn.fIDRatHisto;
194 fSMNumHisto = aCorrFctn.fSMNumHisto;
195 fSMDenHisto = aCorrFctn.fSMDenHisto;
196 fSMRatHisto = aCorrFctn.fSMRatHisto;
197 fCorrectionHisto = aCorrFctn.fCorrectionHisto;
198 fCorrCFHisto = aCorrFctn.fCorrCFHisto;
199 fNumerator = aCorrFctn.fNumerator;
200 fDenominator = aCorrFctn.fDenominator;
201 fRatio = aCorrFctn.fRatio;
202 fQinvHisto = aCorrFctn.fQinvHisto;
203 fLambda = aCorrFctn.fLambda;
204 fRout2 = aCorrFctn.fRout2;
205 fRside2 = aCorrFctn.fRside2;
206 fRlong2 = aCorrFctn.fRlong2;
207 fPairCut = aCorrFctn.fPairCut;
208 fQinvNormLo = aCorrFctn.fQinvNormLo;
209 fQinvNormHi = aCorrFctn.fQinvNormHi;
210 fNumRealsNorm = aCorrFctn.fNumRealsNorm;
211 fNumMixedNorm = aCorrFctn.fNumMixedNorm;
212}
67427ff7 213//____________________________
214AliFemtoBPLCMS3DCorrFctn::~AliFemtoBPLCMS3DCorrFctn(){
215 delete fNumerator;
216 delete fDenominator;
217 delete fRatio;
218 delete fQinvHisto;
219 delete fIDNumHisto;
220 delete fIDDenHisto;
221 delete fIDRatHisto;
222 delete fSMNumHisto;
223 delete fSMDenHisto;
224 delete fSMRatHisto;
225 delete fCorrectionHisto;
226 delete fCorrCFHisto;
227}
0215f606 228//_________________________
229AliFemtoBPLCMS3DCorrFctn& AliFemtoBPLCMS3DCorrFctn::operator=(const AliFemtoBPLCMS3DCorrFctn& aCorrFctn)
230{
231 if (this == &aCorrFctn)
232 return *this;
233 if (fIDNumHisto) delete fIDNumHisto;
234 fIDNumHisto = new TH3D(*aCorrFctn.fIDNumHisto);
235 if (fIDDenHisto) delete fIDDenHisto;
236 fIDDenHisto = new TH3D(*aCorrFctn.fIDDenHisto);
237 if (fIDRatHisto) delete fIDRatHisto;
238 fIDRatHisto = new TH3D(*aCorrFctn.fIDRatHisto);
239 if (fSMNumHisto) delete fSMNumHisto;
240 fSMNumHisto = new TH3D(*aCorrFctn.fSMNumHisto);
241 if (fSMDenHisto) delete fSMDenHisto;
242 fSMDenHisto = new TH3D(*aCorrFctn.fSMDenHisto);
243 if (fSMRatHisto) delete fSMRatHisto;
244 fSMRatHisto = new TH3D(*aCorrFctn.fSMRatHisto);
245
246 if (fCorrectionHisto) delete fCorrectionHisto;
247 fCorrectionHisto = new TH3D(*aCorrFctn.fCorrectionHisto);
248 if (fCorrCFHisto) delete fCorrCFHisto;
249 fCorrCFHisto = new TH3D(*aCorrFctn.fCorrCFHisto);
250 if (fNumerator) delete fNumerator;
251 fNumerator = new TH3D(*aCorrFctn.fNumerator);
252 if (fDenominator) delete fDenominator;
253 fDenominator = new TH3D(*aCorrFctn.fDenominator);
254 if (fRatio) delete fRatio;
255 fRatio = new TH3D(*aCorrFctn.fRatio);
256 if (fQinvHisto) delete fQinvHisto;
257 fQinvHisto = new TH3D(*aCorrFctn.fQinvHisto);
258
259 fLambda = aCorrFctn.fLambda;
260 fRout2 = aCorrFctn.fRout2;
261 fRside2 = aCorrFctn.fRside2;
262 fRlong2 = aCorrFctn.fRlong2;
263 fPairCut = aCorrFctn.fPairCut;
264 fQinvNormLo = aCorrFctn.fQinvNormLo;
265 fQinvNormHi = aCorrFctn.fQinvNormHi;
266 fNumRealsNorm = aCorrFctn.fNumRealsNorm;
267 fNumMixedNorm = aCorrFctn.fNumMixedNorm;
268
269 return *this;
270}
67427ff7 271
272//_________________________
273void AliFemtoBPLCMS3DCorrFctn::WriteOutHistos(){
274
275 fNumerator->Write();
276 fDenominator->Write();
277 // fUncorrectedDenominator->Write();
278 fRatio->Write();
279 fQinvHisto->Write();
280
281 /*
282 if (fSmearPair){
283 fIDNumHisto->Write();
284 fIDDenHisto->Write();
285 fIDRatHisto->Write();
286 //
287 fSMNumHisto->Write();
288 fSMDenHisto->Write();
289 fSMRatHisto->Write();
290 //
291 fCorrectionHisto->Write();
292 fCorrCFHisto->Write();
293 }
294 */
295}
296
297//_________________________
298void AliFemtoBPLCMS3DCorrFctn::Finish(){
299 // here is where we should normalize, fit, etc...
300 double NumFact,DenFact;
301 if ((fNumRealsNorm !=0) && (fNumMixedNorm !=0)){
302 NumFact = double(fNumRealsNorm);
303 DenFact = double(fNumMixedNorm);
304 }
305 // can happen that the fNumRealsNorm and fNumMixedNorm = 0 if you do non-standard
306 // things like making a new CorrFctn and just setting the Numerator and Denominator
307 // from OTHER CorrFctns which you read in (like when doing parallel processing)
308 else{
309 cout << "Warning! - no normalization constants defined - I do the best I can..." << endl;
310 int nbins = fNumerator->GetNbinsX();
311 int half_way = nbins/2;
312 NumFact = fNumerator->Integral(half_way,nbins,half_way,nbins,half_way,nbins);
313 DenFact = fDenominator->Integral(half_way,nbins,half_way,nbins,half_way,nbins);
314 }
315
316 fRatio->Divide(fNumerator,fDenominator,DenFact,NumFact);
317 // fQinvHisto->Divide(fUncorrectedDenominator);
318 fQinvHisto->Divide(fDenominator);
319
320 /*
321 // now do all the resolution correction stuff..
322 if (fSmearPair){ // but only do it if we have been working with a SmearPair
323 fIDRatHisto->Divide(fIDNumHisto,fIDDenHisto);
324 fSMRatHisto->Divide(fSMNumHisto,fSMDenHisto);
325 fCorrectionHisto->Divide(fIDRatHisto,fSMRatHisto);
326 fCorrCFHisto->Multiply(fRatio,fCorrectionHisto);
327 }
328 */
329
330}
331
332//____________________________
333AliFemtoString AliFemtoBPLCMS3DCorrFctn::Report(){
334 string stemp = "LCMS Frame Bertsch-Pratt 3D Correlation Function Report:\n";
335 char ctemp[100];
336 sprintf(ctemp,"Number of entries in numerator:\t%E\n",fNumerator->GetEntries());
337 stemp += ctemp;
338 sprintf(ctemp,"Number of entries in denominator:\t%E\n",fDenominator->GetEntries());
339 stemp += ctemp;
340 sprintf(ctemp,"Number of entries in ratio:\t%E\n",fRatio->GetEntries());
341 stemp += ctemp;
342 sprintf(ctemp,"Normalization region in Qinv was:\t%E\t%E\n",fQinvNormLo,fQinvNormHi);
343 stemp += ctemp;
344 sprintf(ctemp,"Number of pairs in Normalization region was:\n");
345 stemp += ctemp;
346 sprintf(ctemp,"In numerator:\t%lu\t In denominator:\t%lu\n",fNumRealsNorm,fNumMixedNorm);
347 stemp += ctemp;
348 /* if (fCorrection)
349 {
350 float radius = fCorrection->GetRadius();
351 sprintf(ctemp,"Coulomb correction used radius of\t%E\n",radius);
352 }
353 else
354 {
355 sprintf(ctemp,"No Coulomb Correction applied to this CorrFctn\n");
356 }
357 stemp += ctemp;
358 */
359
360 if (fPairCut){
361 sprintf(ctemp,"Here is the PairCut specific to this CorrFctn\n");
362 stemp += ctemp;
363 stemp += fPairCut->Report();
364 }
365 else{
366 sprintf(ctemp,"No PairCut specific to this CorrFctn\n");
367 stemp += ctemp;
368 }
369
370 //
371 AliFemtoString returnThis = stemp;
372 return returnThis;
373}
374//____________________________
375void AliFemtoBPLCMS3DCorrFctn::AddRealPair(const AliFemtoPair* pair){
376
377 if (fPairCut){
378 if (!(fPairCut->Pass(pair))) return;
379 }
380
381 double Qinv = fabs(pair->qInv()); // note - qInv() will be negative for identical pairs...
382 if ((Qinv < fQinvNormHi) && (Qinv > fQinvNormLo)) fNumRealsNorm++;
383 double qOut = fabs(pair->qOutCMS());
384 double qSide = fabs(pair->qSideCMS());
385 double qLong = fabs(pair->qLongCMS());
386
387 fNumerator->Fill(qOut,qSide,qLong);
388}
389//____________________________
390void AliFemtoBPLCMS3DCorrFctn::AddMixedPair(const AliFemtoPair* pair){
391
392 if (fPairCut){
393 if (!(fPairCut->Pass(pair))) return;
394 }
395
396 // double CoulombWeight = (fCorrection ? fCorrection->CoulombCorrect(pair) : 1.0);
397 double CoulombWeight = 1.0;
398
399 double Qinv = fabs(pair->qInv()); // note - qInv() will be negative for identical pairs...
400 if ((Qinv < fQinvNormHi) && (Qinv > fQinvNormLo)) fNumMixedNorm++;
401 double qOut = fabs(pair->qOutCMS());
402 double qSide = fabs(pair->qSideCMS());
403 double qLong = fabs(pair->qLongCMS());
404
405 fDenominator->Fill(qOut,qSide,qLong,CoulombWeight);
406 // fUncorrectedDenominator->Fill(qOut,qSide,qLong,1.0);
407 fQinvHisto->Fill(qOut,qSide,qLong,Qinv);
408
409 /*
410 // now for the momentum resolution stuff...
411 if (fSmearPair){
412 double CorrWeight = 1.0 +
413 fLambda*exp((-qOut*qOut*fRout2 -qSide*qSide*fRside2 -qLong*qLong*fRlong2)/0.038936366329);
414 CorrWeight *= CoulombWeight; // impt.
415
416 fIDNumHisto->Fill(qOut,qSide,qLong,CorrWeight);
417 fIDDenHisto->Fill(qOut,qSide,qLong,CoulombWeight);
418
419 fSmearPair->SetUnsmearedPair(pair);
420 double qOut_prime = fabs(fSmearPair->SmearedPair().qOutCMS());
421 double qSide_prime = fabs(fSmearPair->SmearedPair().qSideCMS());
422 double qLong_prime = fabs(fSmearPair->SmearedPair().qLongCMS());
423
424 fSMNumHisto->Fill(qOut_prime,qSide_prime,qLong_prime,CorrWeight);
425
426 double SmearedCoulombWeight = ( fCorrection ?
427 fCorrection->CoulombCorrect(&(fSmearPair->SmearedPair())) :
428 1.0);
429
430 fSMDenHisto->Fill(qOut_prime,qSide_prime,qLong_prime,SmearedCoulombWeight);
431 }
432 */
433}
434
435