This commit was generated by cvs2svn to compensate for changes in r18145,
[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/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__ 
46 ClassImp(AliFemtoBPLCMS3DCorrFctn)
47 #endif
48
49 //____________________________
50 AliFemtoBPLCMS3DCorrFctn::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 //____________________________
143 AliFemtoBPLCMS3DCorrFctn::~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 //_________________________
159 void 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 //_________________________
184 void 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 //____________________________
219 AliFemtoString 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 //____________________________
261 void 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 //____________________________
276 void 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