8d384c6bff8931ee0ba62f26d7cf86fb8b64a5ea
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemto / AliFemtoBPLCMS3DCorrFctn.cxx
1 ///////////////////////////////////////////////////////////////////////////
2 //                                                                       //
3 // AliFemtoBPLCMS3DCorrFctn: a class to calculate 3D correlation         //
4 // for pairs of identical particles.                                     //
5 // It also stored the weighted qinv per bin histogram for the coulomb    //
6 // correction.                                                           //
7 // In analysis the function should be first created in a macro, then     //
8 // added to the analysis, and at the end of the macro the procedure to   //
9 // write out histograms should be called.                                //
10 //                                                                       //
11 ///////////////////////////////////////////////////////////////////////////
12
13 #include "AliFemtoBPLCMS3DCorrFctn.h"
14 //#include "AliFemtoHisto.h"
15 #include <cstdio>
16
17 #ifdef __ROOT__ 
18 ClassImp(AliFemtoBPLCMS3DCorrFctn)
19 #endif
20
21 //____________________________
22 AliFemtoBPLCMS3DCorrFctn::AliFemtoBPLCMS3DCorrFctn(char* title, const int& nbins, const float& QLo, const float& QHi)
23   :
24   fIDNumHisto(0),
25   fIDDenHisto(0),
26   fIDRatHisto(0),
27   fSMNumHisto(0),
28   fSMDenHisto(0),
29   fSMRatHisto(0),
30   fCorrectionHisto(0),
31   fCorrCFHisto(0),
32   fNumerator(0),
33   fDenominator(0),
34   fRatio(0),
35   fQinvHisto(0),
36   fLambda(0),
37   fRout2(0),
38   fRside2(0),
39   fRlong2(0),
40   fPairCut(0), 
41   fQinvNormLo(0),
42   fQinvNormHi(0),
43   fNumRealsNorm(0),
44   fNumMixedNorm(0)
45 {
46   // Basic constructor
47   // set some stuff...
48   fQinvNormLo = 0.15;
49   fQinvNormHi = 0.18;
50   fNumRealsNorm = 0;
51   fNumMixedNorm = 0;
52   //  fCorrection = 0;  // pointer to Coulomb Correction object
53
54   fPairCut = 0; // added Sept2000 - CorrFctn-specific PairCut
55
56   //  fSmearPair = 0; // no resolution correction unless user sets SmearPair
57
58   // set up numerator
59   char tTitNum[100] = "Num";
60   strcat(tTitNum,title);
61   fNumerator = new TH3D(tTitNum,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
62   // set up denominator
63   char tTitDen[100] = "Den";
64   strcat(tTitDen,title);
65   fDenominator = new TH3D(tTitDen,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
66   // set up uncorrected denominator
67   char tTitDenUncoul[100] = "DenNoCoul";
68   strcat(tTitDenUncoul,title);
69   //  fUncorrectedDenominator = new TH3D(tTitDenUncoul,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
70   // set up ratio
71   char tTitRat[100] = "Rat";
72   strcat(tTitRat,title);
73   fRatio = new TH3D(tTitRat,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
74   // set up ave qInv
75   char tTitQinv[100] = "Qinv";
76   strcat(tTitQinv,title);
77   fQinvHisto = new TH3D(tTitQinv,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
78
79   // to enable error bar calculation...
80   fNumerator->Sumw2();
81   fDenominator->Sumw2();
82   //  fUncorrectedDenominator->Sumw2();
83   fRatio->Sumw2();
84
85   // Following histos are for the momentum resolution correction
86   // they are filled only if a AliFemtoSmear object is plugged in
87   // here comes the "idea" numerator and denominator and ratio...
88   char tTitNumID[100] = "IDNum";
89   strcat(tTitNumID,title);
90   fIDNumHisto = new TH3D(tTitNumID,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
91   char tTitDenID[100] = "IDDen";
92   strcat(tTitDenID,title);
93   fIDDenHisto = new TH3D(tTitDenID,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
94   char tTitRatID[100] = "IDRat";
95   strcat(tTitRatID,title);
96   fIDRatHisto = new TH3D(tTitRatID,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
97
98   fIDNumHisto->Sumw2();
99   fIDDenHisto->Sumw2();
100   fIDRatHisto->Sumw2();
101
102   //
103   // here comes the "smeared" numerator and denominator...
104   char tTitNumSM[100] = "SMNum";
105   strcat(tTitNumSM,title);
106   fSMNumHisto = new TH3D(tTitNumSM,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
107   char tTitDenSM[100] = "SMDen";
108   strcat(tTitDenSM,title);
109   fSMDenHisto = new TH3D(tTitDenSM,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
110   char tTitRatSM[100] = "SMRat";
111   strcat(tTitRatSM,title);
112   fSMRatHisto = new TH3D(tTitRatSM,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
113   //
114   fSMNumHisto->Sumw2();
115   fSMDenHisto->Sumw2();
116   fSMRatHisto->Sumw2();
117   //
118   // here comes the correction factor (which is just ratio of ideal ratio to smeared ratio)
119   char tTitCorrection[100] = "CorrectionFactor";
120   strcat(tTitCorrection,title);
121   fCorrectionHisto = new TH3D(tTitCorrection,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);  
122   fCorrectionHisto->Sumw2();
123   // here comes the fully corrected correlation function
124   char tTitCorrCF[100] = "CorrectedCF";
125   strcat(tTitCorrCF,title);
126   fCorrCFHisto = new TH3D(tTitCorrCF,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
127   fCorrCFHisto->Sumw2();
128
129   // user can (and should) override these defaults...
130   fLambda = 0.6;
131   fRout2 = 6.0*6.0;
132   fRside2 = 6.0*6.0;
133   fRlong2 = 7.0*7.0;
134
135 }
136
137 AliFemtoBPLCMS3DCorrFctn::AliFemtoBPLCMS3DCorrFctn(const AliFemtoBPLCMS3DCorrFctn& aCorrFctn) :
138   AliFemtoCorrFctn(),
139   fIDNumHisto(0),
140   fIDDenHisto(0),
141   fIDRatHisto(0),
142   fSMNumHisto(0),
143   fSMDenHisto(0),
144   fSMRatHisto(0),
145   fCorrectionHisto(0),
146   fCorrCFHisto(0),
147   fNumerator(0),
148   fDenominator(0),
149   fRatio(0),
150   fQinvHisto(0),
151   fLambda(0),
152   fRout2(0),
153   fRside2(0),
154   fRlong2(0),
155   fPairCut(0), 
156   fQinvNormLo(0),
157   fQinvNormHi(0),
158   fNumRealsNorm(0),
159   fNumMixedNorm(0)
160 {
161   // Copy constructor
162   fIDNumHisto = new TH3D(*aCorrFctn.fIDNumHisto);
163   fIDDenHisto = new TH3D(*aCorrFctn.fIDDenHisto);
164   fIDRatHisto = new TH3D(*aCorrFctn.fIDRatHisto);
165   fSMNumHisto = new TH3D(*aCorrFctn.fSMNumHisto);
166   fSMDenHisto = new TH3D(*aCorrFctn.fSMDenHisto);
167   fSMRatHisto = new TH3D(*aCorrFctn.fSMRatHisto);
168   fCorrectionHisto = new TH3D(*aCorrFctn.fCorrectionHisto);
169   fCorrCFHisto = new TH3D(*aCorrFctn.fCorrCFHisto);
170   fNumerator = new TH3D(*aCorrFctn.fNumerator);
171   fDenominator = new TH3D(*aCorrFctn.fDenominator);
172   fRatio = new TH3D(*aCorrFctn.fRatio);
173   fQinvHisto = new TH3D(*aCorrFctn.fQinvHisto);
174   fLambda = aCorrFctn.fLambda;
175   fRout2 = aCorrFctn.fRout2;
176   fRside2 = aCorrFctn.fRside2;
177   fRlong2 = aCorrFctn.fRlong2;
178   fPairCut = aCorrFctn.fPairCut; 
179   fQinvNormLo = aCorrFctn.fQinvNormLo;
180   fQinvNormHi = aCorrFctn.fQinvNormHi;
181   fNumRealsNorm = aCorrFctn.fNumRealsNorm;
182   fNumMixedNorm = aCorrFctn.fNumMixedNorm;
183 }
184 //____________________________
185 AliFemtoBPLCMS3DCorrFctn::~AliFemtoBPLCMS3DCorrFctn(){
186   // Destructor
187   delete fNumerator;
188   delete fDenominator;
189   delete fRatio;
190   delete fQinvHisto;
191   delete fIDNumHisto;
192   delete fIDDenHisto;
193   delete fIDRatHisto;
194   delete fSMNumHisto;
195   delete fSMDenHisto;
196   delete fSMRatHisto;
197   delete fCorrectionHisto;
198   delete fCorrCFHisto;
199 }
200 //_________________________
201 AliFemtoBPLCMS3DCorrFctn& AliFemtoBPLCMS3DCorrFctn::operator=(const AliFemtoBPLCMS3DCorrFctn& aCorrFctn)
202 {
203   // assignment operator
204   if (this == &aCorrFctn)
205     return *this;
206   if (fIDNumHisto) delete fIDNumHisto;
207   fIDNumHisto = new TH3D(*aCorrFctn.fIDNumHisto);
208   if (fIDDenHisto) delete fIDDenHisto;
209   fIDDenHisto = new TH3D(*aCorrFctn.fIDDenHisto);
210   if (fIDRatHisto) delete fIDRatHisto;
211   fIDRatHisto = new TH3D(*aCorrFctn.fIDRatHisto);
212   if (fSMNumHisto) delete fSMNumHisto;
213   fSMNumHisto = new TH3D(*aCorrFctn.fSMNumHisto);
214   if (fSMDenHisto) delete fSMDenHisto;
215   fSMDenHisto = new TH3D(*aCorrFctn.fSMDenHisto);
216   if (fSMRatHisto) delete fSMRatHisto;
217   fSMRatHisto = new TH3D(*aCorrFctn.fSMRatHisto);
218
219   if (fCorrectionHisto) delete fCorrectionHisto;
220   fCorrectionHisto = new TH3D(*aCorrFctn.fCorrectionHisto);
221   if (fCorrCFHisto) delete fCorrCFHisto;
222   fCorrCFHisto = new TH3D(*aCorrFctn.fCorrCFHisto);
223   if (fNumerator) delete fNumerator;
224   fNumerator = new TH3D(*aCorrFctn.fNumerator);
225   if (fDenominator) delete fDenominator;
226   fDenominator = new TH3D(*aCorrFctn.fDenominator);
227   if (fRatio) delete fRatio;
228   fRatio = new TH3D(*aCorrFctn.fRatio);
229   if (fQinvHisto) delete fQinvHisto;
230   fQinvHisto = new TH3D(*aCorrFctn.fQinvHisto);
231
232   fLambda = aCorrFctn.fLambda;
233   fRout2 = aCorrFctn.fRout2;
234   fRside2 = aCorrFctn.fRside2;
235   fRlong2 = aCorrFctn.fRlong2;
236   fPairCut = aCorrFctn.fPairCut; 
237   fQinvNormLo = aCorrFctn.fQinvNormLo;
238   fQinvNormHi = aCorrFctn.fQinvNormHi;
239   fNumRealsNorm = aCorrFctn.fNumRealsNorm;
240   fNumMixedNorm = aCorrFctn.fNumMixedNorm;
241   
242   return *this;
243 }
244
245 //_________________________
246 void AliFemtoBPLCMS3DCorrFctn::WriteOutHistos(){
247   // Write out all histograms to file
248   fNumerator->Write();
249   fDenominator->Write();
250   //  fUncorrectedDenominator->Write();
251   fRatio->Write();
252   fQinvHisto->Write();
253
254   /*
255     if (fSmearPair){
256     fIDNumHisto->Write();
257     fIDDenHisto->Write();
258     fIDRatHisto->Write();
259     //
260     fSMNumHisto->Write();
261     fSMDenHisto->Write();
262     fSMRatHisto->Write();
263     //
264     fCorrectionHisto->Write();
265     fCorrCFHisto->Write();
266     }
267   */
268 }
269 //______________________________
270 TList* AliFemtoBPLCMS3DCorrFctn::GetOutputList()
271 {
272   // Prepare the list of objects to be written to the output
273   TList *tOutputList = new TList();
274
275   tOutputList->Add(fNumerator); 
276   tOutputList->Add(fDenominator);  
277   tOutputList->Add(fQinvHisto);  
278
279   return tOutputList;
280 }
281
282 //_________________________
283 void AliFemtoBPLCMS3DCorrFctn::Finish(){
284   // here is where we should normalize, fit, etc...
285   double tNumFact,tDenFact;
286   if ((fNumRealsNorm !=0) && (fNumMixedNorm !=0)){
287     tNumFact = double(fNumRealsNorm);
288     tDenFact = double(fNumMixedNorm);
289   }
290   // can happen that the fNumRealsNorm and fNumMixedNorm = 0 if you do non-standard
291   //   things like making a new CorrFctn and just setting the Numerator and Denominator
292   //   from OTHER CorrFctns which you read in (like when doing parallel processing) 
293   else{
294     cout << "Warning! - no normalization constants defined - I do the best I can..." << endl;
295     int nbins = fNumerator->GetNbinsX();
296     int half_way = nbins/2;
297     tNumFact = fNumerator->Integral(half_way,nbins,half_way,nbins,half_way,nbins);
298     tDenFact = fDenominator->Integral(half_way,nbins,half_way,nbins,half_way,nbins);
299   }
300
301   fRatio->Divide(fNumerator,fDenominator,tDenFact,tNumFact);
302   //  fQinvHisto->Divide(fUncorrectedDenominator);
303   fQinvHisto->Divide(fDenominator);
304
305   /*
306   // now do all the resolution correction stuff..
307   if (fSmearPair){  // but only do it if we have been working with a SmearPair
308   fIDRatHisto->Divide(fIDNumHisto,fIDDenHisto);
309   fSMRatHisto->Divide(fSMNumHisto,fSMDenHisto);
310   fCorrectionHisto->Divide(fIDRatHisto,fSMRatHisto);
311   fCorrCFHisto->Multiply(fRatio,fCorrectionHisto);
312   }
313   */
314
315 }
316
317 //____________________________
318 AliFemtoString AliFemtoBPLCMS3DCorrFctn::Report(){
319   // Construct the report
320   string stemp = "LCMS Frame Bertsch-Pratt 3D Correlation Function Report:\n";
321   char ctemp[100];
322   sprintf(ctemp,"Number of entries in numerator:\t%E\n",fNumerator->GetEntries());
323   stemp += ctemp;
324   sprintf(ctemp,"Number of entries in denominator:\t%E\n",fDenominator->GetEntries());
325   stemp += ctemp;
326   sprintf(ctemp,"Number of entries in ratio:\t%E\n",fRatio->GetEntries());
327   stemp += ctemp;
328   sprintf(ctemp,"Normalization region in Qinv was:\t%E\t%E\n",fQinvNormLo,fQinvNormHi);
329   stemp += ctemp;
330   sprintf(ctemp,"Number of pairs in Normalization region was:\n");
331   stemp += ctemp;
332   sprintf(ctemp,"In numerator:\t%lu\t In denominator:\t%lu\n",fNumRealsNorm,fNumMixedNorm);
333   stemp += ctemp;
334   /*  if (fCorrection)
335       {
336       float radius = fCorrection->GetRadius();
337       sprintf(ctemp,"Coulomb correction used radius of\t%E\n",radius);
338       }
339       else
340       {
341       sprintf(ctemp,"No Coulomb Correction applied to this CorrFctn\n");
342       }
343       stemp += ctemp;
344   */
345
346   if (fPairCut){
347     sprintf(ctemp,"Here is the PairCut specific to this CorrFctn\n");
348     stemp += ctemp;
349     stemp += fPairCut->Report();
350   }
351   else{
352     sprintf(ctemp,"No PairCut specific to this CorrFctn\n");
353     stemp += ctemp;
354   }
355
356   //  
357   AliFemtoString returnThis = stemp;
358   return returnThis;
359 }
360 //____________________________
361 void AliFemtoBPLCMS3DCorrFctn::AddRealPair( AliFemtoPair* pair){
362   // perform operations on real pairs
363   if (fPairCut){
364     if (!(fPairCut->Pass(pair))) return;
365   }
366
367   double tQinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
368   if ((tQinv < fQinvNormHi) && (tQinv > fQinvNormLo)) fNumRealsNorm++;
369   double qOut = fabs(pair->QOutCMS());
370   double qSide = fabs(pair->QSideCMS());
371   double qLong = fabs(pair->QLongCMS());
372
373   fNumerator->Fill(qOut,qSide,qLong);
374 }
375 //____________________________
376 void AliFemtoBPLCMS3DCorrFctn::AddMixedPair( AliFemtoPair* pair){
377   // perform operations on mixed pairs
378   if (fPairCut){
379     if (!(fPairCut->Pass(pair))) return;
380   }
381
382   //  double CoulombWeight = (fCorrection ? fCorrection->CoulombCorrect(pair) : 1.0);
383   double tCoulombWeight = 1.0;
384
385   double tQinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
386   if ((tQinv < fQinvNormHi) && (tQinv > fQinvNormLo)) fNumMixedNorm++;
387   double qOut = fabs(pair->QOutCMS());
388   double qSide = fabs(pair->QSideCMS());
389   double qLong = fabs(pair->QLongCMS());
390
391   fDenominator->Fill(qOut,qSide,qLong,tCoulombWeight);
392   //  fUncorrectedDenominator->Fill(qOut,qSide,qLong,1.0);
393   fQinvHisto->Fill(qOut,qSide,qLong,tQinv);
394
395   /*
396   // now for the momentum resolution stuff...
397   if (fSmearPair){
398       double CorrWeight =  1.0 + 
399       fLambda*exp((-qOut*qOut*fRout2 -qSide*qSide*fRside2 -qLong*qLong*fRlong2)/0.038936366329);
400     CorrWeight *= CoulombWeight;  // impt.
401
402     fIDNumHisto->Fill(qOut,qSide,qLong,CorrWeight);
403     fIDDenHisto->Fill(qOut,qSide,qLong,CoulombWeight);
404
405     fSmearPair->SetUnsmearedPair(pair);
406     double qOut_prime = fabs(fSmearPair->SmearedPair().qOutCMS());
407     double qSide_prime = fabs(fSmearPair->SmearedPair().qSideCMS());
408     double qLong_prime = fabs(fSmearPair->SmearedPair().qLongCMS());
409
410     fSMNumHisto->Fill(qOut_prime,qSide_prime,qLong_prime,CorrWeight);
411
412     double SmearedCoulombWeight = ( fCorrection ? 
413                                     fCorrection->CoulombCorrect(&(fSmearPair->SmearedPair())) : 
414                                     1.0);
415
416     fSMDenHisto->Fill(qOut_prime,qSide_prime,qLong_prime,SmearedCoulombWeight);
417   }
418   */
419 }
420
421