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