0c68d6d30778a718f39867addc3a335f39055d4f
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemtoUser / AliFemtoModelBPLCMSCorrFctn.cxx
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__ 
13 ClassImp(AliFemtoModelBPLCMSCorrFctn)
14 #endif
15
16 //____________________________
17 AliFemtoModelBPLCMSCorrFctn::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
42 AliFemtoModelBPLCMSCorrFctn::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 //____________________________
52 AliFemtoModelBPLCMSCorrFctn::~AliFemtoModelBPLCMSCorrFctn(){
53   delete fNumerator;
54   delete fDenominator;
55   delete fQinvHisto;
56 }
57 //_________________________
58 AliFemtoModelBPLCMSCorrFctn& 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 //_________________________
73 void AliFemtoModelBPLCMSCorrFctn::WriteOutHistos(){
74
75   fNumerator->Write();
76   fDenominator->Write();
77   fQinvHisto->Write();
78 }
79
80 //_________________________
81 void AliFemtoModelBPLCMSCorrFctn::Finish(){
82   fQinvHisto->Divide(fDenominator);
83 }
84
85 //____________________________
86 AliFemtoString 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 //____________________________
128 void 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 //____________________________
143 void 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