Fix Coverity 10974-82
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemtoUser / AliFemtoBPLCMS3DCorrFctnEMCIC.cxx
CommitLineData
fa35b516 1/***************************************************************************
2 *
3 * $Id: AliFemtoBPLCMS3DCorrFctnEMCIC.cxx $
4 *
5 * Author: Nicolas Bock, Ohio State University, bock@mps.ohio-state.edu
6 ***************************************************************************
7 *
8 * Description: Calculates of the 3D Correlation Function, and also
9 * produces histograms to calculate Energy Momentum Conservation
10 * Induced Correlations (EMCICs)
11 *
12 * This Class produces the following histograms as function of Qinv
13 * (for both real and mixed pairs):
14 * 1) E1 + E2
15 * 2) E1 * E2
16 * 3) Pt1*Pt2
17 * 4) Pz1*Pz2
18 *
19 * The class is derived from AliFemtoBPLCMS3DCorrFctn, therefore it produces
20 * also the histograms in that class.
21 *
22 * NOTE: The EMCIC histograms are not averaged in this class, to obtain
23 * the average, the user needs to divide the real pair histograms by
24 * the numerator, and the mixed pair histograms by the denominator
25 *
26 ***************************************************************************
27 *
28 **************************************************************************/
29
30//////////////////////////////////////////////////////////////////////////
31// //
32// AliFemtoBPLCMS3DCorrFctn: a class to calculate 3D correlation //
33// for pairs of identical particles. //
34// It also stored the weighted qinv per bin histogram for the coulomb //
35// correction. //
36// In analysis the function should be first created in a macro, then //
37// added to the analysis, and at the end of the macro the procedure to //
38// write out histograms should be called. //
39// //
40///////////////////////////////////////////////////////////////////////////
41
42#include "AliFemtoBPLCMS3DCorrFctnEMCIC.h"
43#include "AliFemtoKTPairCut.h"
44#include "AliFemtoAnalysisReactionPlane.h"
45//#include "AliFemtoHisto.h"
46#include <cstdio>
47#include <TVector2.h>
48
49#ifdef __ROOT__
50ClassImp(AliFemtoBPLCMS3DCorrFctnEMCIC)
51#endif
52
53//____________________________
54AliFemtoBPLCMS3DCorrFctnEMCIC::AliFemtoBPLCMS3DCorrFctnEMCIC(char* title, const int& nbins, const float& QLo, const float& QHi)
55:
746bf34d 56AliFemtoCorrFctn(),
943d9a85 57//fEnergyTotalReal(0),
58// fEnergyMultReal(0),
59//fPzMultReal(0),
746bf34d 60//fPtMultReal(0),
61 fNumerator(0),
62 fDenominator(0),
fa35b516 63 fEnergyTotalMix(0),
64 fEnergyMultMix(0),
65 fPzMultMix(0),
746bf34d 66 fPtMultMix(0),
67 fUseRPSelection(0)
fa35b516 68{
69
746bf34d 70 // set up numerator
71 char tTitNum[100] = "Num";
72 strcat(tTitNum,title);
73 fNumerator = new TH3D(tTitNum,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
74 // set up denominator
75 char tTitDen[100] = "Den";
76 strcat(tTitDen,title);
77 fDenominator = new TH3D(tTitDen,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
78
fa35b516 79 //Setup EnergyTotalReal
943d9a85 80 /*char tTitNum1[100] = "ESumReal";
fa35b516 81 strcat(tTitNum1,title);
746bf34d 82 fEnergyTotalReal = new TH3D(tTitNum1,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
fa35b516 83
84 //Setup EnergyMultReal
85 char tTitNum2[100] = "EMultReal";
86 strcat(tTitNum2,title);
746bf34d 87 fEnergyMultReal = new TH3D(tTitNum2,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
fa35b516 88
89 //Setup Pz MultReal
90 char tTitNum3[100] = "PzMultReal";
91 strcat(tTitNum3,title);
746bf34d 92 fPzMultReal = new TH3D(tTitNum3,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
fa35b516 93
94 //Setup Pt MultReal
95 char tTitNum4[100] = "PtMultReal";
96 strcat(tTitNum4,title);
746bf34d 97 fPtMultReal = new TH3D(tTitNum4,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi); */
fa35b516 98
99 //Setup EnergyTotalMix
100 char tTitNum5[100] = "ESumMix";
101 strcat(tTitNum5,title);
746bf34d 102 fEnergyTotalMix = new TH3D(tTitNum5,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
fa35b516 103
104 //Setup EnergyMultMix
105 char tTitNum6[100] = "EMultMix";
106 strcat(tTitNum6,title);
746bf34d 107 fEnergyMultMix = new TH3D(tTitNum6,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
fa35b516 108
109 //Setup Pz MultMix
110 char tTitNum7[100] = "PzMultMix";
111 strcat(tTitNum7,title);
746bf34d 112 fPzMultMix = new TH3D(tTitNum7,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
fa35b516 113
114 //Setup Pt MultMix
115 char tTitNum8[100] = "PtMultMix";
116 strcat(tTitNum8,title);
746bf34d 117 fPtMultMix = new TH3D(tTitNum8,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
fa35b516 118 // To enable error bar calculation
119
943d9a85 120 /*fEnergyTotalReal->Sumw2();
fa35b516 121 fEnergyMultReal->Sumw2();
122 fPzMultReal->Sumw2();
943d9a85 123 fPtMultReal->Sumw2(); */
746bf34d 124 fNumerator->Sumw2();
125 fDenominator->Sumw2();
fa35b516 126 fEnergyTotalMix->Sumw2();
127 fEnergyMultMix->Sumw2();
128 fPzMultMix->Sumw2();
129 fPtMultMix->Sumw2();
130
131}
132
746bf34d 133
134// Variable bin size constructor :
135//qBins array of low-edges for each bin. This is an array of size nbins+1
136AliFemtoBPLCMS3DCorrFctnEMCIC::AliFemtoBPLCMS3DCorrFctnEMCIC(char* title, const int& nbins, const float* qBins):
137AliFemtoCorrFctn(),
138 fNumerator(0),
139 fDenominator(0),
140 fEnergyTotalMix(0),
141 fEnergyMultMix(0),
142 fPzMultMix(0),
143 fPtMultMix(0),
144 fUseRPSelection(0)
145{
146
147 // set up numerator
148 char tTitNum[100] = "Num";
149 strcat(tTitNum,title);
150 fNumerator = new TH3D(tTitNum,title,nbins,qBins,nbins,qBins,nbins,qBins);
151 // set up denominator
152 char tTitDen[100] = "Den";
153 strcat(tTitDen,title);
154 fDenominator = new TH3D(tTitDen,title,nbins,qBins,nbins,qBins,nbins,qBins);
155
156 //Setup EnergyTotalReal
157 /*char tTitNum1[100] = "ESumReal";
158 strcat(tTitNum1,title);
159 fEnergyTotalReal = new TH3D(tTitNum1,title,nbins,qBins,nbins,qBins,nbins,qBins);
160
161 //Setup EnergyMultReal
162 char tTitNum2[100] = "EMultReal";
163 strcat(tTitNum2,title);
164 fEnergyMultReal = new TH3D(tTitNum2,title,nbins,qBins,nbins,qBins,nbins,qBins);
165
166 //Setup Pz MultReal
167 char tTitNum3[100] = "PzMultReal";
168 strcat(tTitNum3,title);
169 fPzMultReal = new TH3D(tTitNum3,title,nbins,qBins,nbins,qBins,nbins,qBins);
170
171 //Setup Pt MultReal
172 char tTitNum4[100] = "PtMultReal";
173 strcat(tTitNum4,title);
174 fPtMultReal = new TH3D(tTitNum4,title,nbins,qBins,nbins,qBins,nbins,qBins); */
175
176 //Setup EnergyTotalMix
177 char tTitNum5[100] = "ESumMix";
178 strcat(tTitNum5,title);
179 fEnergyTotalMix = new TH3D(tTitNum5,title,nbins,qBins,nbins,qBins,nbins,qBins);
180
181 //Setup EnergyMultMix
182 char tTitNum6[100] = "EMultMix";
183 strcat(tTitNum6,title);
184 fEnergyMultMix = new TH3D(tTitNum6,title,nbins,qBins,nbins,qBins,nbins,qBins);
185
186 //Setup Pz MultMix
187 char tTitNum7[100] = "PzMultMix";
188 strcat(tTitNum7,title);
189 fPzMultMix = new TH3D(tTitNum7,title,nbins,qBins,nbins,qBins,nbins,qBins);
190
191 //Setup Pt MultMix
192 char tTitNum8[100] = "PtMultMix";
193 strcat(tTitNum8,title);
194 fPtMultMix = new TH3D(tTitNum8,title,nbins,qBins,nbins,qBins,nbins,qBins);
195 // To enable error bar calculation
196
197 /*fEnergyTotalReal->Sumw2();
198 fEnergyMultReal->Sumw2();
199 fPzMultReal->Sumw2();
200 fPtMultReal->Sumw2(); */
201 fNumerator->Sumw2();
202 fDenominator->Sumw2();
203 fEnergyTotalMix->Sumw2();
204 fEnergyMultMix->Sumw2();
205 fPzMultMix->Sumw2();
206 fPtMultMix->Sumw2();
207}
208
209
210
211
212
213
214
fa35b516 215AliFemtoBPLCMS3DCorrFctnEMCIC::AliFemtoBPLCMS3DCorrFctnEMCIC(const AliFemtoBPLCMS3DCorrFctnEMCIC& aCorrFctn) :
746bf34d 216 AliFemtoCorrFctn(aCorrFctn),
943d9a85 217 /*fEnergyTotalReal(0),
fa35b516 218 fEnergyMultReal(0),
219 fPzMultReal(0),
746bf34d 220 fPtMultReal(0),*/
221 fNumerator(0),
222 fDenominator(0),
fa35b516 223 fEnergyTotalMix (0),
224 fEnergyMultMix (0),
225 fPzMultMix(0),
746bf34d 226 fPtMultMix(0),
227 fUseRPSelection(0)
fa35b516 228{
229 // Copy constructor
746bf34d 230 fNumerator = new TH3D(*aCorrFctn.fNumerator);
231 fDenominator = new TH3D(*aCorrFctn.fDenominator);
943d9a85 232 /* fEnergyTotalReal = new TH3D(*aCorrFctn.fEnergyTotalReal);
fa35b516 233 fEnergyMultReal = new TH3D(*aCorrFctn.fEnergyMultReal);
234 fPzMultReal = new TH3D(*aCorrFctn.fPzMultReal);
943d9a85 235 fPtMultReal = new TH3D(*aCorrFctn.fPtMultReal); */
fa35b516 236 fEnergyTotalMix = new TH3D(*aCorrFctn.fEnergyTotalMix);
237 fEnergyMultMix = new TH3D(*aCorrFctn.fEnergyMultMix);
238 fPzMultMix = new TH3D(*aCorrFctn.fPzMultMix);
239 fPtMultMix = new TH3D(*aCorrFctn.fPtMultMix);
240}
241//____________________________
242AliFemtoBPLCMS3DCorrFctnEMCIC::~AliFemtoBPLCMS3DCorrFctnEMCIC(){
243 // Destructor
943d9a85 244 /* delete fEnergyTotalReal;
fa35b516 245 delete fEnergyMultReal;
246 delete fPzMultReal;
943d9a85 247 delete fPtMultReal; */
746bf34d 248 delete fNumerator;
249 delete fDenominator;
fa35b516 250 delete fEnergyTotalMix;
251 delete fEnergyMultMix;
252 delete fPzMultMix;
253 delete fPtMultMix;
254}
255//_________________________
256AliFemtoBPLCMS3DCorrFctnEMCIC& AliFemtoBPLCMS3DCorrFctnEMCIC::operator=(const AliFemtoBPLCMS3DCorrFctnEMCIC& aCorrFctn)
257{
258 // assignment operator
259 if (this == &aCorrFctn)
260 return *this;
746bf34d 261 if (fNumerator) delete fNumerator;
262 fNumerator = new TH3D(*aCorrFctn.fNumerator);
263 if (fDenominator) delete fDenominator;
264 fDenominator = new TH3D(*aCorrFctn.fDenominator);
fa35b516 265 //Emcics
943d9a85 266 /* if (fEnergyTotalReal) delete fEnergyTotalReal;
fa35b516 267 fEnergyTotalReal = new TH3D(*aCorrFctn.fEnergyTotalReal);
268 if (fEnergyMultReal) delete fEnergyMultReal;
269 fEnergyMultReal = new TH3D(*aCorrFctn.fEnergyMultReal);
270 if (fPzMultReal) delete fPzMultReal;
271 fPzMultReal = new TH3D(*aCorrFctn.fPzMultReal);
272 if (fPtMultReal) delete fPtMultReal;
943d9a85 273 fPtMultReal = new TH3D(*aCorrFctn.fPtMultReal); */
fa35b516 274 if (fEnergyTotalMix) delete fEnergyTotalMix;
275 fEnergyTotalMix = new TH3D(*aCorrFctn.fEnergyTotalMix);
276 if (fEnergyMultMix) delete fEnergyMultMix;
277 fEnergyMultMix = new TH3D(*aCorrFctn.fEnergyMultMix);
278 if (fPzMultMix) delete fPzMultMix;
279 fPzMultMix = new TH3D(*aCorrFctn.fPzMultMix);
280 if (fPtMultMix) delete fPtMultMix;
281 fPtMultMix = new TH3D(*aCorrFctn.fPtMultMix);
282
746bf34d 283 fUseRPSelection = aCorrFctn.fUseRPSelection;
284
fa35b516 285 return *this;
286}
287
288//_________________________
289void AliFemtoBPLCMS3DCorrFctnEMCIC::WriteOutHistos(){
746bf34d 290
291 fNumerator->Write();
292 fDenominator->Write();
943d9a85 293 //fEnergyTotalReal->Write();
294 //fEnergyMultReal->Write();
295 //fPzMultReal->Write();
296 //fPtMultReal->Write();
fa35b516 297 fEnergyTotalMix->Write();
298 fEnergyMultMix->Write();
299 fPzMultMix->Write();
300 fPtMultMix->Write();
746bf34d 301 //cout << "write histos emcics" << endl;
fa35b516 302}
303//______________________________
304TList* AliFemtoBPLCMS3DCorrFctnEMCIC::GetOutputList()
305{
306 // Prepare the list of objects to be written to the output
307 TList *tOutputList = new TList();
308 cout << "Getting list from CFemcic" << endl;
746bf34d 309 tOutputList->Add(fNumerator);
310 tOutputList->Add(fDenominator);
943d9a85 311 /*tOutputList->Add(fEnergyTotalReal);
fa35b516 312 tOutputList->Add(fEnergyMultReal);
313 tOutputList->Add(fPzMultReal);
943d9a85 314 tOutputList->Add(fPtMultReal); */
fa35b516 315 tOutputList->Add(fEnergyTotalMix );
316 tOutputList->Add(fEnergyMultMix );
317 tOutputList->Add(fPzMultMix);
318 tOutputList->Add(fPtMultMix);
319 return tOutputList;
320}
321
322
323
324//____________________________
325void AliFemtoBPLCMS3DCorrFctnEMCIC::AddRealPair( AliFemtoPair* pair){
326 // perform operations on real pairs
327
746bf34d 328 if (fPairCut){
fa35b516 329 if (fUseRPSelection) {
330 AliFemtoKTPairCut *ktc = dynamic_cast<AliFemtoKTPairCut *>(fPairCut);
331 if (!ktc) {
332 cout << "RP aware cut requested, but not connected to the CF" << endl;
333 if (!(fPairCut->Pass(pair))) return;
334 }
335 else {
336 AliFemtoAnalysisReactionPlane *arp = dynamic_cast<AliFemtoAnalysisReactionPlane *> (HbtAnalysis());
337 if (!arp) {
338 cout << "RP aware cut requested, but not connected to the CF" << endl;
339 if (!(fPairCut->Pass(pair))) return;
340 }
1f4d1e46 341 else if (!(ktc->Pass(pair, arp->GetCurrentReactionPlane()))) return;
fa35b516 342 }
343 }
344 else
345 if (!(fPairCut->Pass(pair))) return;
746bf34d 346 }
fa35b516 347
746bf34d 348 double qOut = fabs(pair->QOutCMS());
349 double qSide = fabs( pair->QSideCMS());
350 double qLong = fabs(pair->QLongCMS());
fa35b516 351
746bf34d 352 fNumerator->Fill(qOut,qSide,qLong);
353
943d9a85 354 /*AliFemtoLorentzVector tMom1 = pair->Track1()->FourMomentum();
fa35b516 355 AliFemtoLorentzVector tMom2 = pair->Track2()->FourMomentum();
356 double tE1 = tMom1.e();
357 double tE2 = tMom2.e();
358 double tPz1 = tMom1.pz();
359 double tPz2 = tMom2.pz();
943d9a85 360 TVector2 tPt1;
361 TVector2 tPt2;
fa35b516 362 tPt1.Set(tMom1.px(),tMom1.py());
363 tPt2.Set(tMom2.px(),tMom2.py());
943d9a85 364
365
366 double tPt1DotPt2 = tPt1*tPt2;
fa35b516 367
368 fEnergyTotalReal->Fill(qOut,qSide,qLong,tE1+tE2);
369 fEnergyMultReal->Fill(qOut,qSide,qLong,tE1*tE2);
370 fPzMultReal->Fill(qOut,qSide,qLong,tPz1*tPz2);
943d9a85 371 fPtMultReal->Fill(qOut,qSide,qLong,tPt1DotPt2); */
fa35b516 372}
746bf34d 373
374
fa35b516 375//____________________________
376void AliFemtoBPLCMS3DCorrFctnEMCIC::AddMixedPair( AliFemtoPair* pair){
377 // perform operations on mixed pairs
378// if (fPairCut){
379// if (!(fPairCut->Pass(pair))) return;
380// }
381 if (fPairCut){
382 if (fUseRPSelection) {
383 AliFemtoKTPairCut *ktc = dynamic_cast<AliFemtoKTPairCut *>(fPairCut);
384 if (!ktc) {
385 cout << "RP aware cut requested, but not connected to the CF" << endl;
386 if (!(fPairCut->Pass(pair))) return;
387 }
388 else {
389 AliFemtoAnalysisReactionPlane *arp = dynamic_cast<AliFemtoAnalysisReactionPlane *> (HbtAnalysis());
390 if (!arp) {
391 cout << "RP aware cut requested, but not connected to the CF" << endl;
392 if (!(fPairCut->Pass(pair))) return;
393 }
1f4d1e46 394 else if (!(ktc->Pass(pair, arp->GetCurrentReactionPlane()))) return;
fa35b516 395 }
396 }
397 else
398 if (!(fPairCut->Pass(pair))) return;
399 }
400
401
746bf34d 402 double qOut = fabs(pair->QOutCMS());
403 double qSide = fabs(pair->QSideCMS());
404 double qLong = fabs(pair->QLongCMS());
fa35b516 405
746bf34d 406 fDenominator->Fill(qOut,qSide,qLong);
fa35b516 407
408 AliFemtoLorentzVector tMom1 = pair->Track1()->FourMomentum();
409 AliFemtoLorentzVector tMom2 = pair->Track2()->FourMomentum();
410 double tE1 = tMom1.e();
411 double tE2 = tMom2.e();
412 double tPz1 = tMom1.pz();
413 double tPz2 = tMom2.pz();
943d9a85 414 TVector2 tPt1;
415 TVector2 tPt2;
fa35b516 416 tPt1.Set(tMom1.px(),tMom1.py());
417 tPt2.Set(tMom2.px(),tMom2.py());
943d9a85 418 double tPt1DotPt2 = tPt1*tPt2;
fa35b516 419
420 fEnergyTotalMix->Fill(qOut,qSide,qLong,tE1+tE2);
421 fEnergyMultMix->Fill(qOut,qSide,qLong,tE1*tE2);
422 fPzMultMix->Fill(qOut,qSide,qLong,tPz1*tPz2);
423 fPtMultMix->Fill(qOut,qSide,qLong,tPt1DotPt2);
424
425}
426
746bf34d 427
428void AliFemtoBPLCMS3DCorrFctnEMCIC::SetUseRPSelection(unsigned short aRPSel)
429{
430 fUseRPSelection = aRPSel;
431}