Making the directory structure of AliFemtoUser flat. All files go into one common...
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemtoUser / AliFemtoModelBPLCMSCorrFctn.cxx
CommitLineData
65423af9 1////////////////////////////////////////////////////////////////////////////////
2/// ///
3/// AliFemtoModelBPLCMSCorrFctn - the class for correlation function which ///
4/// uses the model framework and weight generation and calculated the 3D ///
5/// correlation function in the Bertsh-Pratt LCMS system ///
6/// Authors: Adam Kisiel, kisiel@mps.ohio-state.edu ///
7/// ///
8////////////////////////////////////////////////////////////////////////////////
9#include "AliFemtoModelBPLCMSCorrFctn.h"
10#include <cstdio>
11
12#ifdef __ROOT__
13ClassImp(AliFemtoModelBPLCMSCorrFctn)
14#endif
15
16//____________________________
17AliFemtoModelBPLCMSCorrFctn::AliFemtoModelBPLCMSCorrFctn(char* title, const int& nbins, const float& QLo, const float& QHi)
18 :
19 fNumerator(0),
20 fDenominator(0),
21 fQinvHisto(0)
22{
23
24 // set up numerator
25 char TitNum[100] = "Num";
26 strcat(TitNum,title);
27 fNumerator = new TH3D(TitNum,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
28 // set up denominator
29 char TitDen[100] = "Den";
30 strcat(TitDen,title);
31 fDenominator = new TH3D(TitDen,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
32 // set up ave qInv
33 char TitQinv[100] = "Qinv";
34 strcat(TitQinv,title);
35 fQinvHisto = new TH3D(TitQinv,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
36
37 // to enable error bar calculation...
38 fNumerator->Sumw2();
39 fDenominator->Sumw2();
40}
41
42AliFemtoModelBPLCMSCorrFctn::AliFemtoModelBPLCMSCorrFctn(const AliFemtoModelBPLCMSCorrFctn& aCorrFctn) :
43 fNumerator(0),
44 fDenominator(0),
45 fQinvHisto(0)
46{
47 fNumerator = new TH3D(*aCorrFctn.fNumerator);
48 fDenominator = new TH3D (*aCorrFctn.fDenominator);
49 fQinvHisto = new TH3D(*aCorrFctn.fQinvHisto);
50}
51//____________________________
52AliFemtoModelBPLCMSCorrFctn::~AliFemtoModelBPLCMSCorrFctn(){
53 delete fNumerator;
54 delete fDenominator;
55 delete fQinvHisto;
56}
57//_________________________
58AliFemtoModelBPLCMSCorrFctn& AliFemtoModelBPLCMSCorrFctn::operator=(const AliFemtoModelBPLCMSCorrFctn& aCorrFctn)
59{
60 if (this == &aCorrFctn)
61 return *this;
62 if (fNumerator) delete fNumerator;
63 fNumerator = new TH3D(*aCorrFctn.fNumerator);
64 if (fDenominator) delete fDenominator;
65 fDenominator = new TH3D(*aCorrFctn.fDenominator);
66 if (fQinvHisto) delete fQinvHisto;
67 fQinvHisto = new TH3D(*aCorrFctn.fQinvHisto);
68
69 return *this;
70}
71
72//_________________________
73void AliFemtoModelBPLCMSCorrFctn::WriteOutHistos(){
74
75 fNumerator->Write();
76 fDenominator->Write();
77 fQinvHisto->Write();
78}
79
80//_________________________
81void AliFemtoModelBPLCMSCorrFctn::Finish(){
82 fQinvHisto->Divide(fDenominator);
83}
84
85//____________________________
86AliFemtoString AliFemtoModelBPLCMSCorrFctn::Report(){
87 string stemp = "LCMS Frame Bertsch-Pratt 3D Correlation Function Report:\n";
88 char ctemp[100];
89 sprintf(ctemp,"Number of entries in numerator:\t%E\n",fNumerator->GetEntries());
90 stemp += ctemp;
91 sprintf(ctemp,"Number of entries in denominator:\t%E\n",fDenominator->GetEntries());
92 stemp += ctemp;
93 sprintf(ctemp,"Number of entries in ratio:\t%E\n",fRatio->GetEntries());
94 stemp += ctemp;
95 sprintf(ctemp,"Normalization region in Qinv was:\t%E\t%E\n",fQinvNormLo,fQinvNormHi);
96 stemp += ctemp;
97 sprintf(ctemp,"Number of pairs in Normalization region was:\n");
98 stemp += ctemp;
99 sprintf(ctemp,"In numerator:\t%lu\t In denominator:\t%lu\n",fNumRealsNorm,fNumMixedNorm);
100 stemp += ctemp;
101 /* if (fCorrection)
102 {
103 float radius = fCorrection->GetRadius();
104 sprintf(ctemp,"Coulomb correction used radius of\t%E\n",radius);
105 }
106 else
107 {
108 sprintf(ctemp,"No Coulomb Correction applied to this CorrFctn\n");
109 }
110 stemp += ctemp;
111 */
112
113 if (fPairCut){
114 sprintf(ctemp,"Here is the PairCut specific to this CorrFctn\n");
115 stemp += ctemp;
116 stemp += fPairCut->Report();
117 }
118 else{
119 sprintf(ctemp,"No PairCut specific to this CorrFctn\n");
120 stemp += ctemp;
121 }
122
123 //
124 AliFemtoString returnThis = stemp;
125 return returnThis;
126}
127//____________________________
128void AliFemtoModelBPLCMSCorrFctn::AddRealPair( AliFemtoPair* pair){
129
130 if (fPairCut){
131 if (!(fPairCut->Pass(pair))) return;
132 }
133
134 double Qinv = fabs(pair->qInv()); // note - qInv() will be negative for identical pairs...
135 if ((Qinv < fQinvNormHi) && (Qinv > fQinvNormLo)) fNumRealsNorm++;
136 double qOut = fabs(pair->qOutCMS());
137 double qSide = fabs(pair->qSideCMS());
138 double qLong = fabs(pair->qLongCMS());
139
140 fNumerator->Fill(qOut,qSide,qLong);
141}
142//____________________________
143void AliFemtoModelBPLCMSCorrFctn::AddMixedPair( AliFemtoPair* pair){
144
145 if (fPairCut){
146 if (!(fPairCut->Pass(pair))) return;
147 }
148
149 // double CoulombWeight = (fCorrection ? fCorrection->CoulombCorrect(pair) : 1.0);
150 double CoulombWeight = 1.0;
151
152 double Qinv = fabs(pair->qInv()); // note - qInv() will be negative for identical pairs...
153 if ((Qinv < fQinvNormHi) && (Qinv > fQinvNormLo)) fNumMixedNorm++;
154 double qOut = fabs(pair->qOutCMS());
155 double qSide = fabs(pair->qSideCMS());
156 double qLong = fabs(pair->qLongCMS());
157
158 fDenominator->Fill(qOut,qSide,qLong,CoulombWeight);
159 // fUncorrectedDenominator->Fill(qOut,qSide,qLong,1.0);
160 fQinvHisto->Fill(qOut,qSide,qLong,Qinv);
161
162 /*
163 // now for the momentum resolution stuff...
164 if (fSmearPair){
165 double CorrWeight = 1.0 +
166 fLambda*exp((-qOut*qOut*fRout2 -qSide*qSide*fRside2 -qLong*qLong*fRlong2)/0.038936366329);
167 CorrWeight *= CoulombWeight; // impt.
168
169 fIDNumHisto->Fill(qOut,qSide,qLong,CorrWeight);
170 fIDDenHisto->Fill(qOut,qSide,qLong,CoulombWeight);
171
172 fSmearPair->SetUnsmearedPair(pair);
173 double qOut_prime = fabs(fSmearPair->SmearedPair().qOutCMS());
174 double qSide_prime = fabs(fSmearPair->SmearedPair().qSideCMS());
175 double qLong_prime = fabs(fSmearPair->SmearedPair().qLongCMS());
176
177 fSMNumHisto->Fill(qOut_prime,qSide_prime,qLong_prime,CorrWeight);
178
179 double SmearedCoulombWeight = ( fCorrection ?
180 fCorrection->CoulombCorrect(&(fSmearPair->SmearedPair())) :
181 1.0);
182
183 fSMDenHisto->Fill(qOut_prime,qSide_prime,qLong_prime,SmearedCoulombWeight);
184 }
185 */
186}
187
188