Fix AliFemtoModelFreezeOutGenerator undefined references
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemto / CorrFctn / AliFemtoBPLCMS3DCorrFctn.cxx
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/04/25 15:38:41  panos
15  * Importing the HBT code dir
16  *
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__ 
49 ClassImp(AliFemtoBPLCMS3DCorrFctn)
50 #endif
51
52 //____________________________
53 AliFemtoBPLCMS3DCorrFctn::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 {
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
168 AliFemtoBPLCMS3DCorrFctn::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 }
213 //____________________________
214 AliFemtoBPLCMS3DCorrFctn::~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 }
228 //_________________________
229 AliFemtoBPLCMS3DCorrFctn& 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 }
271
272 //_________________________
273 void 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 //_________________________
298 void 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 //____________________________
333 AliFemtoString 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 //____________________________
375 void 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 //____________________________
390 void 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