Lines getting the matched track moved to a method in AliCalorimeterUtils. Lines copie...
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemtoUser / AliFemtoModelCorrFctn3DSpherical.cxx
1 ///////////////////////////////////////////////////////////////////////////
2 //                                                                       //
3 // AliFemtoModelCorrFctn3DSpherical: a class to calculate 3D correlation //
4 // for pairs of identical particles, binned in spherical coordinates.    //
5 // In analysis the function should be first created in a macro, then     //
6 // added to the analysis, and at the end of the macro the procedure to   //
7 // write out histograms should be called.                                //
8 //                                                                       //
9 ///////////////////////////////////////////////////////////////////////////
10
11 #include "AliFemtoModelCorrFctn3DSpherical.h"
12 #include "AliFemtoModelManager.h"
13 #include <TMath.h>
14 #include <cstdio>
15 //#include <Math/SpecFunc.h>
16
17 #ifdef __ROOT__ 
18 ClassImp(AliFemtoModelCorrFctn3DSpherical)
19 #endif
20
21 //____________________________
22 AliFemtoModelCorrFctn3DSpherical::AliFemtoModelCorrFctn3DSpherical(char* title, const int& nqbins, const float& QLo, const float& QHi, const int& nphibins, const int& ncthetabins):
23   AliFemtoModelCorrFctn(),
24   fTrueNumeratorSph(0),
25   fFakeNumeratorSph(0),
26   fDenominatorSph(0),
27   fPairCut(0x0)
28 {
29   // set up numerator
30   char tTitNum[101] = "NumTrue";
31   strncat(tTitNum,title, 100);
32   fTrueNumeratorSph = new TH3D(tTitNum,title,nqbins,QLo,QHi,nphibins,-TMath::Pi(),TMath::Pi(),ncthetabins,-1.0,1.0);
33   // set up numerator
34   char tTitNumF[101] = "NumFake";
35   strncat(tTitNumF,title, 100);
36   fFakeNumeratorSph = new TH3D(tTitNumF,title,nqbins,QLo,QHi,nphibins,-TMath::Pi(),TMath::Pi(),ncthetabins,-1.0,1.0);
37   // set up denominator
38   char tTitDen[101] = "Den";
39   strncat(tTitDen,title, 100);
40   fDenominatorSph = new TH3D(tTitDen,title,nqbins,QLo,QHi,nphibins,-TMath::Pi(),TMath::Pi(),ncthetabins,-1.0,1.0);
41
42   // to enable error bar calculation...
43   fTrueNumeratorSph->Sumw2();
44   fFakeNumeratorSph->Sumw2();
45   fDenominatorSph->Sumw2();
46 }
47
48 AliFemtoModelCorrFctn3DSpherical::AliFemtoModelCorrFctn3DSpherical(const AliFemtoModelCorrFctn3DSpherical& aCorrFctn) :
49   AliFemtoModelCorrFctn(aCorrFctn),
50   fTrueNumeratorSph(0),
51   fFakeNumeratorSph(0),
52   fDenominatorSph(0),
53   fPairCut(0x0)
54 {
55   // Copy constructor
56   fTrueNumeratorSph = new TH3D(*aCorrFctn.fTrueNumeratorSph);
57   fFakeNumeratorSph = new TH3D(*aCorrFctn.fFakeNumeratorSph);
58   fDenominatorSph = new TH3D(*aCorrFctn.fDenominatorSph);
59   fPairCut = aCorrFctn.fPairCut;
60 }
61 //____________________________
62 AliFemtoModelCorrFctn3DSpherical::~AliFemtoModelCorrFctn3DSpherical(){
63   // Destructor
64   delete fTrueNumeratorSph;
65   delete fFakeNumeratorSph;
66   delete fDenominatorSph;
67 }
68 //_________________________
69 AliFemtoModelCorrFctn3DSpherical& AliFemtoModelCorrFctn3DSpherical::operator=(const AliFemtoModelCorrFctn3DSpherical& aCorrFctn)
70 {
71   // assignment operator
72   if (this == &aCorrFctn)
73     return *this;
74
75   if (fTrueNumeratorSph) delete fTrueNumeratorSph;
76   fTrueNumeratorSph = new TH3D(*aCorrFctn.fTrueNumeratorSph);
77   if (fFakeNumeratorSph) delete fFakeNumeratorSph;
78   fFakeNumeratorSph = new TH3D(*aCorrFctn.fFakeNumeratorSph);
79   if (fDenominatorSph) delete fDenominatorSph;
80   fDenominatorSph = new TH3D(*aCorrFctn.fDenominatorSph);
81   
82   fPairCut = aCorrFctn.fPairCut;
83   
84   return *this;
85 }
86
87 //_________________________
88 void AliFemtoModelCorrFctn3DSpherical::WriteOutHistos(){
89   // Write out all histograms to file
90   fTrueNumeratorSph->Write();
91   fFakeNumeratorSph->Write();
92   fDenominatorSph->Write();
93 }
94 //______________________________
95 TList* AliFemtoModelCorrFctn3DSpherical::GetOutputList()
96 {
97   // Prepare the list of objects to be written to the output
98   TList *tOutputList = new TList();
99
100   tOutputList->Add(fTrueNumeratorSph); 
101   tOutputList->Add(fFakeNumeratorSph); 
102   tOutputList->Add(fDenominatorSph);  
103
104   return tOutputList;
105 }
106
107 //_________________________
108 void AliFemtoModelCorrFctn3DSpherical::Finish(){
109   // here is where we should normalize, fit, etc...
110 }
111
112 //____________________________
113 AliFemtoString AliFemtoModelCorrFctn3DSpherical::Report(){
114   // Construct the report
115   string stemp = "PRF Frame Spherical 3D Model Correlation Function Report:\n";
116   char ctemp[100];
117   snprintf(ctemp , 100, "Number of entries in numerator:\t%E\n",fTrueNumeratorSph->GetEntries());
118   stemp += ctemp;
119   snprintf(ctemp , 100, "Number of entries in denominator:\t%E\n",fDenominatorSph->GetEntries());
120   stemp += ctemp;
121
122   if (fPairCut){
123     snprintf(ctemp , 100, "Here is the PairCut specific to this CorrFctn\n");
124     stemp += ctemp;
125     stemp += fPairCut->Report();
126   }
127   else{
128     snprintf(ctemp , 100, "No PairCut specific to this CorrFctn\n");
129     stemp += ctemp;
130   }
131
132   //  
133   AliFemtoString returnThis = stemp;
134   return returnThis;
135 }
136 //____________________________
137 void AliFemtoModelCorrFctn3DSpherical::AddRealPair( AliFemtoPair* pair){
138   // perform operations on real pairs
139   if (fPairCut){
140     if (!(fPairCut->Pass(pair))) return;
141   }
142
143   Double_t weight = fManager->GetWeight(pair);
144
145   double tKO = pair->KOut();
146   double tKS = pair->KSide();
147   double tKL = pair->KLong();
148
149   double tKR = sqrt(tKO*tKO + tKS*tKS + tKL*tKL);
150   double tKC;
151   if ( fabs(tKR) < 1e-10 ) tKC = 0.0;
152   else tKC=tKL/tKR;
153   double tKP=atan2(tKS,tKO);
154
155   fTrueNumeratorSph->Fill(tKR,tKP,tKC,weight);
156 }
157 //____________________________
158 void AliFemtoModelCorrFctn3DSpherical::AddMixedPair( AliFemtoPair* pair){
159   // perform operations on mixed pairs
160   if (fPairCut){
161     if (!(fPairCut->Pass(pair))) return;
162   }
163
164   Double_t weight = fManager->GetWeight(pair);
165
166   double tKO = pair->KOut();
167   double tKS = pair->KSide();
168   double tKL = pair->KLong();
169
170   double tKR = sqrt(tKO*tKO + tKS*tKS + tKL*tKL);
171   double tKC;
172   if ( fabs(tKR) < 1e-10 ) tKC = 0.0;
173   else tKC=tKL/tKR;
174   double tKP=atan2(tKS,tKO);
175
176   fFakeNumeratorSph->Fill(tKR,tKP,tKC,weight);
177   fDenominatorSph->Fill(tKR,tKP,tKC);
178 }
179
180