More mods to suit ROOT head
[u/mrichter/AliRoot.git] / PWG2 / EBYE / AliBalance.cxx
CommitLineData
6c178944 1/**************************************************************************
2 * Author: Panos Christakoglou. *
3 * Contributors are mentioned in the code where appropriate. *
4 * *
5 * Permission to use, copy, modify and distribute this software and its *
6 * documentation strictly for non-commercial purposes is hereby granted *
7 * without fee, provided that the above copyright notice appears in all *
8 * copies and that both the copyright notice and this permission notice *
9 * appear in the supporting documentation. The authors make no claims *
10 * about the suitability of this software for any purpose. It is *
11 * provided "as is" without express or implied warranty. *
12 **************************************************************************/
13
14/* $Id$ */
15
16//-----------------------------------------------------------------
17// Balance Function class
18// This is the class to deal with the Balance Function analysis
19// Origin: Panos Christakoglou, UOA-CERN, Panos.Christakoglou@cern.ch
20//-----------------------------------------------------------------
21
22
23//ROOT
24#include <Riostream.h>
25#include <TMath.h>
d560b581 26#include <TAxis.h>
6c178944 27#include <TLorentzVector.h>
28#include <TGraphErrors.h>
29
30#include "AliBalance.h"
31
32ClassImp(AliBalance)
33
7f0257ea 34//____________________________________________________________________//
5d534bd4 35AliBalance::AliBalance() :
36 TObject(),
37 fCharge(0), fNtrack(0), fV(0), fNumberOfBins(0),
38 fAnalysisType(0), fAnalyzedEvents(0), fP2Start(0),
39 fP2Stop(0), fP2Step(0), fNn(0), fNp(0) {
6c178944 40 // Default constructor
7f0257ea 41 for(Int_t i = 0; i < MAXIMUM_NUMBER_OF_STEPS; i++) {
42 fNpp[i] = .0;
43 fNnn[i] = .0;
44 fNpn[i] = .0;
45 fB[i] = 0.0;
46 ferror[i] = 0.0;
47 }
6c178944 48}
49
7f0257ea 50//____________________________________________________________________//
5d534bd4 51AliBalance::AliBalance(Double_t p2Start, Double_t p2Stop, Int_t p2Bins) :
52 TObject(),
53 fCharge(0), fNtrack(0), fV(0),
54 fNumberOfBins(p2Bins), fAnalysisType(0),
55 fAnalyzedEvents(0), fP2Start(p2Start), fP2Stop(p2Stop),
56 fP2Step(TMath::Abs(fP2Start - fP2Stop) / (Double_t)fNumberOfBins),
57 fNn(0), fNp(0) {
6c178944 58 // Constructor
7f0257ea 59 for(Int_t i = 0; i < MAXIMUM_NUMBER_OF_STEPS; i++) {
60 fNpp[i] = .0;
61 fNnn[i] = .0;
62 fNpn[i] = .0;
63 fB[i] = 0.0;
64 ferror[i] = 0.0;
65 }
6c178944 66}
67
7f0257ea 68//____________________________________________________________________//
69AliBalance::AliBalance(const AliBalance& balance):
e559055a 70 TObject(balance),
7f0257ea 71 fNtrack(balance.fNtrack),
72 fNumberOfBins(balance.fNumberOfBins),
73 fAnalysisType(balance.fAnalysisType),
74 fAnalyzedEvents(balance.fAnalyzedEvents),
75 fP2Start(balance.fP2Start),
76 fP2Stop(balance.fP2Stop),
77 fP2Step(balance.fP2Step),
78 fNn(balance.fNn),
79 fNp(balance.fNp),
80 fCharge(0),
81 fV(0)
6c178944 82{
7f0257ea 83 //copy constructor
84
85 for(Int_t i = 0; i < MAXIMUM_NUMBER_OF_STEPS; i++) {
86 fNpp[i] = .0;
87 fNnn[i] = .0;
88 fNpn[i] = .0;
89 fB[i] = 0.0;
90 ferror[i] = 0.0;
91 }
92
93 if (balance.fV) fV = new TLorentzVector(*(balance.fV));
94 if (balance.fCharge) fCharge = new Double_t(*(balance.fCharge));
95}
96
97//____________________________________________________________________//
98AliBalance::~AliBalance() {
6c178944 99 // Destructor
7f0257ea 100 delete fV;
101 delete fCharge;
6c178944 102}
103
7f0257ea 104//____________________________________________________________________//
105void AliBalance::SetNnn(Double_t *nn) {
106 // Setter of the Nnn term
107 for(Int_t i = 0; i < fNumberOfBins; i++) fNnn[i] = nn[i];
108}
109
110//____________________________________________________________________//
111void AliBalance::SetNpp(Double_t *pp) {
112 // Setter of the Npp term
113 for(Int_t i = 0; i < fNumberOfBins; i++) fNpp[i] = pp[i];
114}
115
116//____________________________________________________________________//
117void AliBalance::SetNpn(Double_t *pn) {
118 // Setter of the Npn term
119 for(Int_t i = 0; i < fNumberOfBins; i++) fNpn[i] = pn[i];
120}
121
122//____________________________________________________________________//
123void AliBalance::SetNumberOfBins(Int_t ibins) {
6c178944 124 // Sets the number of bins for the analyzed interval
7f0257ea 125 fNumberOfBins = ibins;
6c178944 126}
127
7f0257ea 128//____________________________________________________________________//
129void AliBalance::SetInterval(Double_t p2Start, Double_t p2Stop) {
6c178944 130 // Sets the analyzed interval.
131 // The analysis variable is set by SetAnalysisType
132 fP2Start = p2Start;
133 fP2Stop = p2Stop;
134 fP2Step = TMath::Abs(p2Start - p2Stop) / (Double_t)fNumberOfBins;
135}
136
7f0257ea 137//____________________________________________________________________//
138void AliBalance::SetAnalysisType(Int_t iType) {
6c178944 139 //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
140 this->fAnalysisType = iType;
7f0257ea 141 if(fAnalysisType==0) {
142 cout<<" ====================== "<<endl;
143 cout<<"||Analysis selected: y||"<<endl;
144 cout<<" ====================== "<<endl;
145 }
146 else if(fAnalysisType==1) {
147 cout<<" ======================== "<<endl;
148 cout<<"||Analysis selected: eta||"<<endl;
149 cout<<" ======================== "<<endl;
150 }
151 else if(fAnalysisType==2) {
152 cout<<" ========================== "<<endl;
153 cout<<"||Analysis selected: Qlong||"<<endl;
154 cout<<" ========================== "<<endl;
155 }
156 else if(fAnalysisType==3) {
157 cout<<" ========================= "<<endl;
158 cout<<"||Analysis selected: Qout||"<<endl;
159 cout<<" ========================= "<<endl;
160 }
161 else if(fAnalysisType==4) {
162 cout<<" ========================== "<<endl;
163 cout<<"||Analysis selected: Qside||"<<endl;
164 cout<<" ========================== "<<endl;
165 }
166 else if(fAnalysisType==5) {
167 cout<<" ========================= "<<endl;
168 cout<<"||Analysis selected: Qinv||"<<endl;
169 cout<<" ========================= "<<endl;
170 }
171 else if(fAnalysisType==6) {
172 cout<<" ======================== "<<endl;
173 cout<<"||Analysis selected: phi||"<<endl;
174 cout<<" ======================== "<<endl;
175 }
176 else {
177 cout<<"Selection of analysis mode failed!!!"<<endl;
178 cout<<"Choices are: 0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi"<<endl;
179 abort();
180 }
6c178944 181}
182
7f0257ea 183//____________________________________________________________________//
184void AliBalance::SetParticles(TLorentzVector *P, Double_t *charge, Int_t dim) {
6c178944 185 // Sets a new particle with given 4-momentum and charge.
186 // dim is the size of the array of charges and corresponds
187 // to the number of selected tracks.
188 this->fV = P;
189 this->fCharge = charge;
190 fNtrack = dim;
191}
192
193
7f0257ea 194//____________________________________________________________________//
195void AliBalance::CalculateBalance() {
6c178944 196 // Calculates the balance function
197 fAnalyzedEvents++;
198 Int_t i = 0 , j = 0;
199 Int_t ibin = 0;
200
7f0257ea 201 for(i = 0; i < fNtrack; i++) {
202 if(fCharge[i] > 0) fNp += 1.;
203 if(fCharge[i] < 0) fNn += 1.;
204 }
6c178944 205
206 //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
7f0257ea 207 if(fAnalysisType==0) {
208 for(i = 1; i < fNtrack; i++) {
209 for(j = 0; j < i; j++) {
210 Double_t rap1 = 0.5*log((fV[i].E() + fV[i].Pz())/(fV[i].E() - fV[i].Pz()));
211 Double_t rap2 = 0.5*log((fV[j].E() + fV[j].Pz())/(fV[j].E() - fV[j].Pz()));
212 Double_t dy = TMath::Abs(rap1 - rap2);
213 ibin = Int_t(dy/fP2Step);
214 if((fCharge[i] > 0)&&(fCharge[j] > 0)) fNpp[ibin] += 1.;
215 if((fCharge[i] < 0)&&(fCharge[j] < 0)) fNnn[ibin] += 1.;
216 if((fCharge[i] > 0)&&(fCharge[j] < 0)) fNpn[ibin] += 1.;
217 if((fCharge[i] < 0)&&(fCharge[j] > 0)) fNpn[ibin] += 1.;
218 }
219 }
220 }//case 0
221 if(fAnalysisType==1) {
222 for(i = 1; i < fNtrack; i++) {
223 for(j = 0; j < i; j++) {
224 Double_t p1 = sqrt(pow(fV[i].Px(),2) + pow(fV[i].Py(),2) + pow(fV[i].Pz(),2));
225 Double_t p2 = sqrt(pow(fV[j].Px(),2) + pow(fV[j].Py(),2) + pow(fV[j].Pz(),2));
226 Double_t eta1 = 0.5*log((p1 + fV[i].Pz())/(p1 - fV[i].Pz()));
227 Double_t eta2 = 0.5*log((p2 + fV[j].Pz())/(p2 - fV[j].Pz()));
228 Double_t deta = TMath::Abs(eta1 - eta2);
229 ibin = Int_t(deta/fP2Step);
230 if((fCharge[i] > 0)&&(fCharge[j] > 0)) fNpp[ibin] += 1.;
231 if((fCharge[i] < 0)&&(fCharge[j] < 0)) fNnn[ibin] += 1.;
232 if((fCharge[i] > 0)&&(fCharge[j] < 0)) fNpn[ibin] += 1.;
233 if((fCharge[i] < 0)&&(fCharge[j] > 0)) fNpn[ibin] += 1.;
234 }
235 }
236 }//case 1
237 if(fAnalysisType==2) {
238 for(i = 1; i < fNtrack; i++) {
239 for(j = 0; j < i; j++) {
240 Double_t eTot = fV[i].E() + fV[j].E();
241 Double_t pxTot = fV[i].Px() + fV[j].Px();
242 Double_t pyTot = fV[i].Py() + fV[j].Py();
243 Double_t pzTot = fV[i].Pz() + fV[j].Pz();
244 Double_t q0Tot = fV[i].E() - fV[j].E();
245 Double_t qzTot = fV[i].Pz() - fV[j].Pz();
246 Double_t snn = pow(eTot,2) - pow(pxTot,2) - pow(pyTot,2) - pow(pzTot,2);
247 Double_t ptTot = sqrt( pow(pxTot,2) + pow(pyTot,2));
248 Double_t qLong = TMath::Abs(eTot*qzTot - pzTot*q0Tot)/sqrt(snn + pow(ptTot,2));
249 ibin = Int_t(qLong/fP2Step);
250 //cout<<ibin<<endl;
251 if((fCharge[i] > 0)&&(fCharge[j] > 0)) fNpp[ibin] += 1.;
252 if((fCharge[i] < 0)&&(fCharge[j] < 0)) fNnn[ibin] += 1.;
253 if((fCharge[i] > 0)&&(fCharge[j] < 0)) fNpn[ibin] += 1.;
254 if((fCharge[i] < 0)&&(fCharge[j] > 0)) fNpn[ibin] += 1.;
255 }
256 }
257 }//case 2
258 if(fAnalysisType==3) {
259 for(i = 1; i < fNtrack; i++) {
260 for(j = 0; j < i; j++) {
261 Double_t eTot = fV[i].E() + fV[j].E();
262 Double_t pxTot = fV[i].Px() + fV[j].Px();
263 Double_t pyTot = fV[i].Py() + fV[j].Py();
264 Double_t pzTot = fV[i].Pz() + fV[j].Pz();
265 Double_t qxTot = fV[i].Px() - fV[j].Px();
266 Double_t qyTot = fV[i].Py() - fV[j].Py();
267 Double_t snn = pow(eTot,2) - pow(pxTot,2) - pow(pyTot,2) - pow(pzTot,2);
268 Double_t ptTot = sqrt( pow(pxTot,2) + pow(pyTot,2));
269 Double_t qOut = sqrt(snn/(snn + pow(ptTot,2))) * TMath::Abs(pxTot*qxTot + pyTot*qyTot)/ptTot;
270 ibin = Int_t(qOut/fP2Step);
271 //cout<<ibin<<endl;
272 if((fCharge[i] > 0)&&(fCharge[j] > 0)) fNpp[ibin] += 1.;
273 if((fCharge[i] < 0)&&(fCharge[j] < 0)) fNnn[ibin] += 1.;
274 if((fCharge[i] > 0)&&(fCharge[j] < 0)) fNpn[ibin] += 1.;
275 if((fCharge[i] < 0)&&(fCharge[j] > 0)) fNpn[ibin] += 1.;
276 }
277 }
278 }//case 3
279 if(fAnalysisType==4) {
280 for(i = 1; i < fNtrack; i++) {
281 for(j = 0; j < i; j++) {
282 Double_t pxTot = fV[i].Px() + fV[j].Px();
283 Double_t pyTot = fV[i].Py() + fV[j].Py();
284 Double_t qxTot = fV[i].Px() - fV[j].Px();
285 Double_t qyTot = fV[i].Py() - fV[j].Py();
286 Double_t ptTot = sqrt( pow(pxTot,2) + pow(pyTot,2));
287 Double_t qSide = TMath::Abs(pxTot*qyTot - pyTot*qxTot)/ptTot;
288 ibin = Int_t(qSide/fP2Step);
289 //cout<<ibin<<endl;
290 if((fCharge[i] > 0)&&(fCharge[j] > 0)) fNpp[ibin] += 1.;
291 if((fCharge[i] < 0)&&(fCharge[j] < 0)) fNnn[ibin] += 1.;
292 if((fCharge[i] > 0)&&(fCharge[j] < 0)) fNpn[ibin] += 1.;
293 if((fCharge[i] < 0)&&(fCharge[j] > 0)) fNpn[ibin] += 1.;
294 }
295 }
296 }//case 4
297 if(fAnalysisType==5) {
298 for(i = 1; i < fNtrack; i++) {
299 for(j = 0; j < i; j++) {
300 Double_t q0Tot = fV[i].E() - fV[j].E();
301 Double_t qxTot = fV[i].Px() - fV[j].Px();
302 Double_t qyTot = fV[i].Py() - fV[j].Py();
303 Double_t qzTot = fV[i].Pz() - fV[j].Pz();
304 Double_t qInv = sqrt(TMath::Abs(-pow(q0Tot,2) +pow(qxTot,2) +pow(qyTot,2) +pow(qzTot,2)));
305 ibin = Int_t(qInv/fP2Step);
306 //cout<<ibin<<endl;
307 if((fCharge[i] > 0)&&(fCharge[j] > 0)) fNpp[ibin] += 1.;
308 if((fCharge[i] < 0)&&(fCharge[j] < 0)) fNnn[ibin] += 1.;
309 if((fCharge[i] > 0)&&(fCharge[j] < 0)) fNpn[ibin] += 1.;
310 if((fCharge[i] < 0)&&(fCharge[j] > 0)) fNpn[ibin] += 1.;
311 }
312 }
313 }//case 5
314 if(fAnalysisType==6) {
315 for(i = 1; i < fNtrack; i++) {
316 for(j = 0; j < i; j++) {
317 Double_t phi1 = TMath::ATan(fV[i].Py()/fV[i].Px())*180.0/TMath::Pi();
318 Double_t phi2 = TMath::ATan(fV[j].Py()/fV[j].Px())*180.0/TMath::Pi();
319 Double_t dphi = TMath::Abs(phi1 - phi2);
320 ibin = Int_t(dphi/fP2Step);
321 if((fCharge[i] > 0)&&(fCharge[j] > 0)) fNpp[ibin] += 1.;
322 if((fCharge[i] < 0)&&(fCharge[j] < 0)) fNnn[ibin] += 1.;
323 if((fCharge[i] > 0)&&(fCharge[j] < 0)) fNpn[ibin] += 1.;
324 if((fCharge[i] < 0)&&(fCharge[j] > 0)) fNpn[ibin] += 1.;
325 }
326 }
327 }//case 6
6c178944 328}
329
7f0257ea 330//____________________________________________________________________//
331Double_t AliBalance::GetBalance(Int_t p2) {
6c178944 332 // Returns the value of the balance function in bin p2
333 fB[p2] = 0.5*(((fNpn[p2] - 2.0*fNnn[p2])/fNn) + ((fNpn[p2] - 2.0*fNpp[p2])/fNp))/fP2Step;
7f0257ea 334
6c178944 335 return fB[p2];
336}
337
7f0257ea 338//____________________________________________________________________//
339Double_t AliBalance::GetError(Int_t p2) {
6c178944 340 // Returns the error on the BF value for bin p2
341 ferror[p2] = sqrt( Double_t(fNpp[p2])/(Double_t(fNp)*Double_t(fNp)) + Double_t(fNnn[p2])/(Double_t(fNn)*Double_t(fNn)) + Double_t(fNpn[p2])*pow((0.5/Double_t(fNp) + 0.5/Double_t(fNn)),2))/fP2Step;
342
343 return ferror[p2];
344}
345
7f0257ea 346//____________________________________________________________________//
347void AliBalance::PrintResults() {
6c178944 348 // Prints the results
349 Double_t x[MAXIMUM_NUMBER_OF_STEPS];
350 Double_t fSumXi = 0.0, fSumBi = 0.0, fSumBiXi = 0.0;
351 Double_t fSumBiXi2 = 0.0, fSumBi2Xi2 = 0.0;
352 Double_t fSumDeltaBi2 = 0.0, fSumXi2DeltaBi2 = 0.0;
353 Double_t deltaBalP2 = 0.0, integral = 0.0;
354 Double_t deltaErrorNew = 0.0;
355
356 cout<<"=================================================="<<endl;
7f0257ea 357 for(Int_t i = 0; i < fNumberOfBins; i++) {
358 x[i] = fP2Start + fP2Step*i + fP2Step/2;
359 cout<<"B: "<<fB[i]<<"\t Error: "<<ferror[i]<<"\t bin: "<<x[i]<<endl;
360 }
6c178944 361 cout<<"=================================================="<<endl;
7f0257ea 362 for(Int_t i = 1; i < fNumberOfBins; i++) {
363 fSumXi += x[i];
364 fSumBi += fB[i];
365 fSumBiXi += fB[i]*x[i];
366 fSumBiXi2 += fB[i]*pow(x[i],2);
367 fSumBi2Xi2 += pow(fB[i],2)*pow(x[i],2);
368 fSumDeltaBi2 += pow(ferror[i],2);
369 fSumXi2DeltaBi2 += pow(x[i],2) * pow(ferror[i],2);
370
371 deltaBalP2 += fP2Step*pow(ferror[i],2);
372 integral += fP2Step*fB[i];
373 }
374 for(Int_t i = 1; i < fNumberOfBins; i++) deltaErrorNew += ferror[i]*(x[i]*fSumBi - fSumBiXi)/pow(fSumBi,2);
375
376 Double_t integralError = sqrt(deltaBalP2);
6c178944 377
7f0257ea 378 Double_t delta = fSumBiXi / fSumBi;
379 Double_t deltaError = (fSumBiXi / fSumBi) * sqrt(pow((sqrt(fSumXi2DeltaBi2)/fSumBiXi),2) + pow((fSumDeltaBi2/fSumBi),2) );
6c178944 380
381 cout<<"Analyzed events: "<<fAnalyzedEvents<<endl;
382 cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
383 cout<<"New error: "<<deltaErrorNew<<endl;
384 cout<<"Interval: "<<integral<<"\t Error: "<<integralError<<endl;
385 cout<<"=================================================="<<endl;
386}
387
7f0257ea 388//____________________________________________________________________//
389TGraphErrors *AliBalance::DrawBalance() {
6c178944 390 // Draws the BF
391 Double_t x[MAXIMUM_NUMBER_OF_STEPS];
392 Double_t xer[MAXIMUM_NUMBER_OF_STEPS];
393 Double_t b[MAXIMUM_NUMBER_OF_STEPS];
394 Double_t ber[MAXIMUM_NUMBER_OF_STEPS];
395
7f0257ea 396 if((fNp == 0)||(fNn == 0)) {
397 cout<<"Couldn't find any particles in the analyzed interval!!!"<<endl;
398 cout<<"Aborting....."<<endl;
399 abort();
400 }
401
402 for(Int_t i = 0; i < fNumberOfBins; i++) {
403 b[i] = GetBalance(i);
404 ber[i] = GetError(i);
405 x[i] = fP2Start + fP2Step*i + fP2Step/2;
406 xer[i] = 0.0;
407 }
408
409 TGraphErrors *gr = new TGraphErrors(fNumberOfBins,x,b,xer,ber);
410 gr->SetMarkerStyle(25);
6c178944 411 gr->GetXaxis()->SetTitleColor(1);
7f0257ea 412 if(fAnalysisType==0) {
413 gr->GetXaxis()->SetTitle("#Delta y");
414 gr->GetYaxis()->SetTitle("B(#Delta y)");
415 }
416 if(fAnalysisType==1) {
417 gr->GetXaxis()->SetTitle("#Delta #eta");
418 gr->GetYaxis()->SetTitle("B(#Delta #eta)");
419 }
420 if(fAnalysisType==2) {
421 gr->GetXaxis()->SetTitle("Q_{long} [GeV]");
422 gr->GetYaxis()->SetTitle("B(Q_{long})");
423 }
424 if(fAnalysisType==3) {
425 gr->GetXaxis()->SetTitle("Q_{out} [GeV]");
426 gr->GetYaxis()->SetTitle("B(Q_{out})");
427 }
428 if(fAnalysisType==4) {
429 gr->GetXaxis()->SetTitle("Q_{side} [GeV]");
430 gr->GetYaxis()->SetTitle("B(Q_{side})");
431 }
432 if(fAnalysisType==5) {
433 gr->GetXaxis()->SetTitle("Q_{inv} [GeV]");
434 gr->GetYaxis()->SetTitle("B(Q_{inv})");
435 }
436 if(fAnalysisType==6) {
437 gr->GetXaxis()->SetTitle("#Delta #phi");
438 gr->GetYaxis()->SetTitle("B(#Delta #phi)");
439 }
6c178944 440
441 return gr;
442}