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