1 ///////////////////////////////////////////////////////////////////////////
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 //
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. //
11 ///////////////////////////////////////////////////////////////////////////
13 #include "AliFemtoBPLCMS3DCorrFctn.h"
14 #include "AliFemtoKTPairCut.h"
15 #include "AliFemtoAnalysisReactionPlane.h"
16 //#include "AliFemtoHisto.h"
20 ClassImp(AliFemtoBPLCMS3DCorrFctn)
23 //____________________________
24 AliFemtoBPLCMS3DCorrFctn::AliFemtoBPLCMS3DCorrFctn(char* title, const int& nbins, const float& QLo, const float& QHi)
33 // fCorrectionHisto(0),
51 fQinvNormLo = (QHi-QLo)*0.8;
52 fQinvNormHi = (QHi-QLo)*0.8;
55 // fCorrection = 0; // pointer to Coulomb Correction object
57 // fSmearPair = 0; // no resolution correction unless user sets SmearPair
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);
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);
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);
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);
80 // to enable error bar calculation...
82 fDenominator->Sumw2();
83 // fUncorrectedDenominator->Sumw2();
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);
99 // fIDNumHisto->Sumw2();
100 // fIDDenHisto->Sumw2();
101 // fIDRatHisto->Sumw2();
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);
115 // fSMNumHisto->Sumw2();
116 // fSMDenHisto->Sumw2();
117 // fSMRatHisto->Sumw2();
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();
130 // user can (and should) override these defaults...
138 AliFemtoBPLCMS3DCorrFctn::AliFemtoBPLCMS3DCorrFctn(const AliFemtoBPLCMS3DCorrFctn& aCorrFctn) :
139 AliFemtoCorrFctn(aCorrFctn),
146 // fCorrectionHisto(0),
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;
185 //____________________________
186 AliFemtoBPLCMS3DCorrFctn::~AliFemtoBPLCMS3DCorrFctn(){
192 // delete fIDNumHisto;
193 // delete fIDDenHisto;
194 // delete fIDRatHisto;
195 // delete fSMNumHisto;
196 // delete fSMDenHisto;
197 // delete fSMRatHisto;
198 // delete fCorrectionHisto;
199 // delete fCorrCFHisto;
201 //_________________________
202 AliFemtoBPLCMS3DCorrFctn& AliFemtoBPLCMS3DCorrFctn::operator=(const AliFemtoBPLCMS3DCorrFctn& aCorrFctn)
204 // assignment operator
205 if (this == &aCorrFctn)
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);
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);
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;
246 //_________________________
247 void AliFemtoBPLCMS3DCorrFctn::WriteOutHistos(){
248 // Write out all histograms to file
250 fDenominator->Write();
251 // fUncorrectedDenominator->Write();
257 fIDNumHisto->Write();
258 fIDDenHisto->Write();
259 fIDRatHisto->Write();
261 fSMNumHisto->Write();
262 fSMDenHisto->Write();
263 fSMRatHisto->Write();
265 fCorrectionHisto->Write();
266 fCorrCFHisto->Write();
270 //______________________________
271 TList* AliFemtoBPLCMS3DCorrFctn::GetOutputList()
273 // Prepare the list of objects to be written to the output
274 TList *tOutputList = new TList();
276 tOutputList->Add(fNumerator);
277 tOutputList->Add(fDenominator);
278 tOutputList->Add(fQinvHisto);
283 //_________________________
284 void 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);
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)
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);
302 fRatio->Divide(fNumerator,fDenominator,tDenFact,tNumFact);
303 // fQinvHisto->Divide(fUncorrectedDenominator);
304 fQinvHisto->Divide(fDenominator);
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);
318 //____________________________
319 AliFemtoString AliFemtoBPLCMS3DCorrFctn::Report(){
320 // Construct the report
321 string stemp = "LCMS Frame Bertsch-Pratt 3D Correlation Function Report:\n";
323 snprintf(ctemp , 100, "Number of entries in numerator:\t%E\n",fNumerator->GetEntries());
325 snprintf(ctemp , 100, "Number of entries in denominator:\t%E\n",fDenominator->GetEntries());
327 snprintf(ctemp , 100, "Number of entries in ratio:\t%E\n",fRatio->GetEntries());
329 snprintf(ctemp , 100, "Normalization region in Qinv was:\t%E\t%E\n",fQinvNormLo,fQinvNormHi);
331 snprintf(ctemp , 100, "Number of pairs in Normalization region was:\n");
333 snprintf(ctemp , 100, "In numerator:\t%lu\t In denominator:\t%lu\n",fNumRealsNorm,fNumMixedNorm);
337 float radius = fCorrection->GetRadius();
338 snprintf(ctemp , 100, "Coulomb correction used radius of\t%E\n",radius);
342 snprintf(ctemp , 100, "No Coulomb Correction applied to this CorrFctn\n");
348 snprintf(ctemp , 100, "Here is the PairCut specific to this CorrFctn\n");
350 stemp += fPairCut->Report();
353 snprintf(ctemp , 100, "No PairCut specific to this CorrFctn\n");
358 AliFemtoString returnThis = stemp;
361 //____________________________
362 void AliFemtoBPLCMS3DCorrFctn::AddRealPair( AliFemtoPair* pair){
363 // perform operations on real pairs
365 if (fUseRPSelection) {
366 AliFemtoKTPairCut *ktc = dynamic_cast<AliFemtoKTPairCut *>(fPairCut);
368 cout << "RP aware cut requested, but not connected to the CF" << endl;
369 if (!(fPairCut->Pass(pair))) return;
372 AliFemtoAnalysisReactionPlane *arp = dynamic_cast<AliFemtoAnalysisReactionPlane *> (HbtAnalysis());
374 cout << "RP aware cut requested, but not connected to the CF" << endl;
375 if (!(fPairCut->Pass(pair))) return;
377 else if (!(ktc->Pass(pair, arp->GetCurrentReactionPlane()))) return;
381 if (!(fPairCut->Pass(pair))) return;
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());
390 fNumerator->Fill(qOut,qSide,qLong);
392 //____________________________
393 void AliFemtoBPLCMS3DCorrFctn::AddMixedPair( AliFemtoPair* pair){
394 // perform operations on mixed pairs
396 // if (!(fPairCut->Pass(pair))) return;
399 if (fUseRPSelection) {
400 AliFemtoKTPairCut *ktc = dynamic_cast<AliFemtoKTPairCut *>(fPairCut);
402 cout << "RP aware cut requested, but not connected to the CF" << endl;
403 if (!(fPairCut->Pass(pair))) return;
406 AliFemtoAnalysisReactionPlane *arp = dynamic_cast<AliFemtoAnalysisReactionPlane *> (HbtAnalysis());
408 cout << "RP aware cut requested, but not connected to the CF" << endl;
409 if (!(fPairCut->Pass(pair))) return;
411 else if (!(ktc->Pass(pair, arp->GetCurrentReactionPlane()))) return;
415 if (!(fPairCut->Pass(pair))) return;
418 // double CoulombWeight = (fCorrection ? fCorrection->CoulombCorrect(pair) : 1.0);
419 double tCoulombWeight = 1.0;
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());
427 fDenominator->Fill(qOut,qSide,qLong,tCoulombWeight);
428 // fUncorrectedDenominator->Fill(qOut,qSide,qLong,1.0);
429 fQinvHisto->Fill(qOut,qSide,qLong,tQinv);
432 // now for the momentum resolution stuff...
434 double CorrWeight = 1.0 +
435 fLambda*exp((-qOut*qOut*fRout2 -qSide*qSide*fRside2 -qLong*qLong*fRlong2)/0.038936366329);
436 CorrWeight *= CoulombWeight; // impt.
438 fIDNumHisto->Fill(qOut,qSide,qLong,CorrWeight);
439 fIDDenHisto->Fill(qOut,qSide,qLong,CoulombWeight);
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());
446 fSMNumHisto->Fill(qOut_prime,qSide_prime,qLong_prime,CorrWeight);
448 double SmearedCoulombWeight = ( fCorrection ?
449 fCorrection->CoulombCorrect(&(fSmearPair->SmearedPair())) :
452 fSMDenHisto->Fill(qOut_prime,qSide_prime,qLong_prime,SmearedCoulombWeight);
458 void AliFemtoBPLCMS3DCorrFctn::SetUseRPSelection(unsigned short aRPSel)
460 fUseRPSelection = aRPSel;