Lines getting the matched track moved to a method in AliCalorimeterUtils. Lines copie...
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemto / AliFemtoBPLCMS3DCorrFctn.cxx
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__ 
20 ClassImp(AliFemtoBPLCMS3DCorrFctn)
21 #endif
22
23 //____________________________
24 AliFemtoBPLCMS3DCorrFctn::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
138 AliFemtoBPLCMS3DCorrFctn::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 //____________________________
186 AliFemtoBPLCMS3DCorrFctn::~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 //_________________________
202 AliFemtoBPLCMS3DCorrFctn& 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 //_________________________
247 void 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 //______________________________
271 TList* 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 //_________________________
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);
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 //____________________________
319 AliFemtoString 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 //____________________________
362 void 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 //____________________________
393 void 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
458 void AliFemtoBPLCMS3DCorrFctn::SetUseRPSelection(unsigned short aRPSel)
459 {
460   fUseRPSelection = aRPSel;
461 }