]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/EBYE/AliBalance.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
First step in moving the GRP reco-params from the reconstruction macros to OCDB.
[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>
5c33329d 28#include <TObjArray.h>
6c178944 29#include <TGraphErrors.h>
9d1f0df5 30#include <TString.h>
6c178944 31
5c33329d 32#include "AliESDtrack.h"
6c178944 33#include "AliBalance.h"
34
35ClassImp(AliBalance)
36
7f0257ea 37//____________________________________________________________________//
5d534bd4 38AliBalance::AliBalance() :
39 TObject(),
9d1f0df5 40 fAnalysisLevel("ESD"), fNumberOfBins(0),
5d534bd4 41 fAnalysisType(0), fAnalyzedEvents(0), fP2Start(0),
42 fP2Stop(0), fP2Step(0), fNn(0), fNp(0) {
6c178944 43 // Default constructor
7f0257ea 44 for(Int_t i = 0; i < MAXIMUM_NUMBER_OF_STEPS; i++) {
45 fNpp[i] = .0;
46 fNnn[i] = .0;
47 fNpn[i] = .0;
48 fB[i] = 0.0;
49 ferror[i] = 0.0;
50 }
6c178944 51}
52
7f0257ea 53//____________________________________________________________________//
5d534bd4 54AliBalance::AliBalance(Double_t p2Start, Double_t p2Stop, Int_t p2Bins) :
9d1f0df5 55 TObject(), fAnalysisLevel("ESD"),
5d534bd4 56 fNumberOfBins(p2Bins), fAnalysisType(0),
57 fAnalyzedEvents(0), fP2Start(p2Start), fP2Stop(p2Stop),
58 fP2Step(TMath::Abs(fP2Start - fP2Stop) / (Double_t)fNumberOfBins),
59 fNn(0), fNp(0) {
6c178944 60 // Constructor
7f0257ea 61 for(Int_t i = 0; i < MAXIMUM_NUMBER_OF_STEPS; i++) {
62 fNpp[i] = .0;
63 fNnn[i] = .0;
64 fNpn[i] = .0;
65 fB[i] = 0.0;
66 ferror[i] = 0.0;
67 }
6c178944 68}
69
7f0257ea 70//____________________________________________________________________//
71AliBalance::AliBalance(const AliBalance& balance):
9d1f0df5 72 TObject(balance), fAnalysisLevel(balance.fAnalysisLevel),
7f0257ea 73 fNumberOfBins(balance.fNumberOfBins),
74 fAnalysisType(balance.fAnalysisType),
75 fAnalyzedEvents(balance.fAnalyzedEvents),
76 fP2Start(balance.fP2Start),
77 fP2Stop(balance.fP2Stop),
78 fP2Step(balance.fP2Step),
79 fNn(balance.fNn),
5c33329d 80 fNp(balance.fNp) {
7f0257ea 81 //copy constructor
7f0257ea 82 for(Int_t i = 0; i < MAXIMUM_NUMBER_OF_STEPS; i++) {
83 fNpp[i] = .0;
84 fNnn[i] = .0;
85 fNpn[i] = .0;
86 fB[i] = 0.0;
87 ferror[i] = 0.0;
88 }
7f0257ea 89}
90
91//____________________________________________________________________//
92AliBalance::~AliBalance() {
6c178944 93 // Destructor
94}
95
7f0257ea 96//____________________________________________________________________//
97void AliBalance::SetNnn(Double_t *nn) {
98 // Setter of the Nnn term
99 for(Int_t i = 0; i < fNumberOfBins; i++) fNnn[i] = nn[i];
100}
101
102//____________________________________________________________________//
103void AliBalance::SetNpp(Double_t *pp) {
104 // Setter of the Npp term
105 for(Int_t i = 0; i < fNumberOfBins; i++) fNpp[i] = pp[i];
106}
107
108//____________________________________________________________________//
109void AliBalance::SetNpn(Double_t *pn) {
110 // Setter of the Npn term
111 for(Int_t i = 0; i < fNumberOfBins; i++) fNpn[i] = pn[i];
112}
113
114//____________________________________________________________________//
115void AliBalance::SetNumberOfBins(Int_t ibins) {
6c178944 116 // Sets the number of bins for the analyzed interval
7f0257ea 117 fNumberOfBins = ibins;
6c178944 118}
119
7f0257ea 120//____________________________________________________________________//
121void AliBalance::SetInterval(Double_t p2Start, Double_t p2Stop) {
6c178944 122 // Sets the analyzed interval.
123 // The analysis variable is set by SetAnalysisType
124 fP2Start = p2Start;
125 fP2Stop = p2Stop;
126 fP2Step = TMath::Abs(p2Start - p2Stop) / (Double_t)fNumberOfBins;
127}
128
7f0257ea 129//____________________________________________________________________//
130void AliBalance::SetAnalysisType(Int_t iType) {
6c178944 131 //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
132 this->fAnalysisType = iType;
7f0257ea 133 if(fAnalysisType==0) {
134 cout<<" ====================== "<<endl;
135 cout<<"||Analysis selected: y||"<<endl;
136 cout<<" ====================== "<<endl;
137 }
138 else if(fAnalysisType==1) {
139 cout<<" ======================== "<<endl;
140 cout<<"||Analysis selected: eta||"<<endl;
141 cout<<" ======================== "<<endl;
142 }
143 else if(fAnalysisType==2) {
144 cout<<" ========================== "<<endl;
145 cout<<"||Analysis selected: Qlong||"<<endl;
146 cout<<" ========================== "<<endl;
147 }
148 else if(fAnalysisType==3) {
149 cout<<" ========================= "<<endl;
150 cout<<"||Analysis selected: Qout||"<<endl;
151 cout<<" ========================= "<<endl;
152 }
153 else if(fAnalysisType==4) {
154 cout<<" ========================== "<<endl;
155 cout<<"||Analysis selected: Qside||"<<endl;
156 cout<<" ========================== "<<endl;
157 }
158 else if(fAnalysisType==5) {
159 cout<<" ========================= "<<endl;
160 cout<<"||Analysis selected: Qinv||"<<endl;
161 cout<<" ========================= "<<endl;
162 }
163 else if(fAnalysisType==6) {
164 cout<<" ======================== "<<endl;
165 cout<<"||Analysis selected: phi||"<<endl;
166 cout<<" ======================== "<<endl;
167 }
168 else {
169 cout<<"Selection of analysis mode failed!!!"<<endl;
170 cout<<"Choices are: 0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi"<<endl;
171 abort();
172 }
6c178944 173}
174
898d5b28 175//____________________________________________________________________//
176const char* AliBalance::GetAnalysisType() {
177 //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
178 TString analysisType;
179 switch(fAnalysisType) {
180 case 0:
181 analysisType = "Rapidity";
182 break;
183 case 1:
184 analysisType = "Pseudo-rapidity";
185 break;
186 case 2:
187 analysisType = "Qlong";
188 break;
189 case 3:
190 analysisType = "Qout";
191 break;
192 case 4:
193 analysisType = "Qside";
194 break;
195 case 5:
196 analysisType = "Qinv";
197 break;
198 case 6:
199 analysisType = "Phi";
200 break;
201 default:
202 break;
203 }
204 analysisType += "\nInterval: ";
205 analysisType += fP2Start; analysisType += " - "; analysisType += fP2Stop;
206 analysisType += "\nSteps: "; analysisType += fP2Step;
207 analysisType += "\nBins: "; analysisType += fNumberOfBins;
208
209 return analysisType.Data();
210}
211
7f0257ea 212//____________________________________________________________________//
5c33329d 213void AliBalance::CalculateBalance(TObjArray *gTrackArray) {
6c178944 214 // Calculates the balance function
215 fAnalyzedEvents++;
216 Int_t i = 0 , j = 0;
217 Int_t ibin = 0;
5c33329d 218
9d1f0df5 219 AliESDtrack* track = 0;
220 AliESDtrack* track1 = 0;
221 AliESDtrack* track2 = 0;
222
5c33329d 223 Int_t gNtrack = gTrackArray->GetEntries();
224 for(i = 0; i < gNtrack; i++) {
9d1f0df5 225 if(fAnalysisLevel == "ESD")
226 track = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
5c33329d 227 Short_t charge = track->Charge();
228 if(charge > 0) fNp += 1.;
229 if(charge < 0) fNn += 1.;
7f0257ea 230 }
5c33329d 231 //Printf("Np: %lf - Nn: %lf",fNp,fNn);
6c178944 232
233 //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
7f0257ea 234 if(fAnalysisType==0) {
5c33329d 235 for(i = 1; i < gNtrack; i++) {
9d1f0df5 236 if(fAnalysisLevel == "ESD")
237 track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
5c33329d 238 Short_t charge1 = track1->Charge();
239 Double_t pZ1 = track1->Pz();
240 Double_t energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
241 TMath::Power(track1->M(),2));
7f0257ea 242 for(j = 0; j < i; j++) {
9d1f0df5 243 if(fAnalysisLevel == "ESD")
244 track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
5c33329d 245 Short_t charge2 = track2->Charge();
246 Double_t pZ2 = track2->Pz();
247 Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
248 TMath::Power(track2->M(),2));
249
250 Double_t rap1 = 0.5*log((energy1 + pZ1)/(energy1 - pZ1));
251 Double_t rap2 = 0.5*log((energy2 + pZ2)/(energy2 - pZ2));
7f0257ea 252 Double_t dy = TMath::Abs(rap1 - rap2);
253 ibin = Int_t(dy/fP2Step);
5c33329d 254 if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
255 if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
256 if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
257 if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
7f0257ea 258 }
259 }
260 }//case 0
261 if(fAnalysisType==1) {
5c33329d 262 for(i = 1; i < gNtrack; i++) {
9d1f0df5 263 if(fAnalysisLevel == "ESD")
264 track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
5c33329d 265 Short_t charge1 = track1->Charge();
266 Double_t pZ1 = track1->Pz();
267 Double_t p1 = track1->P();
7f0257ea 268 for(j = 0; j < i; j++) {
9d1f0df5 269 if(fAnalysisLevel == "ESD")
270 track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
5c33329d 271 Short_t charge2 = track2->Charge();
272 Double_t pZ2 = track2->Pz();
273 Double_t p2 = track2->P();
274 Double_t eta1 = 0.5*log((p1 + pZ1)/(p1 - pZ1));
275 Double_t eta2 = 0.5*log((p2 + pZ2)/(p2 - pZ2));
7f0257ea 276 Double_t deta = TMath::Abs(eta1 - eta2);
277 ibin = Int_t(deta/fP2Step);
5c33329d 278
279 if((charge1 > 0.)&&(charge2 > 0.)) fNpp[ibin] += 1.;
280 if((charge1 < 0.)&&(charge2 < 0.)) fNnn[ibin] += 1.;
281 if((charge1 > 0.)&&(charge2 < 0.)) fNpn[ibin] += 1.;
282 if((charge1 < 0.)&&(charge2 > 0.)) fNpn[ibin] += 1.;
283 //Printf("charge1: %d - eta1: %lf - charge2: %d - eta2: %lf - deta: %lf - ibin: %d - fNpp: %lf - fNnn: %lf - fNpn: %lf",charge1,eta1,charge2,eta2,deta,ibin,fNpp[ibin],fNnn[ibin],fNpn[ibin]);
7f0257ea 284 }
285 }
286 }//case 1
287 if(fAnalysisType==2) {
5c33329d 288 for(i = 1; i < gNtrack; i++) {
9d1f0df5 289 if(fAnalysisLevel == "ESD")
290 track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
5c33329d 291 Short_t charge1 = track1->Charge();
292 Double_t pX1 = track1->Px();
293 Double_t pY1 = track1->Py();
294 Double_t pZ1 = track1->Pz();
295 Double_t energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
296 TMath::Power(track1->M(),2));
7f0257ea 297 for(j = 0; j < i; j++) {
9d1f0df5 298 if(fAnalysisLevel == "ESD")
299 track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
5c33329d 300 Short_t charge2 = track2->Charge();
301 Double_t pX2 = track2->Px();
302 Double_t pY2 = track2->Py();
303 Double_t pZ2 = track2->Pz();
304 Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
305 TMath::Power(track2->M(),2));
306 Double_t eTot = energy1 + energy2;
307 Double_t pxTot = pX1 + pX2;
308 Double_t pyTot = pY1 + pY2;
309 Double_t pzTot = pZ1 + pZ2;
310 Double_t q0Tot = energy1 - energy2;
311 Double_t qzTot = pZ1 - pZ2;
312 Double_t snn = TMath::Power(eTot,2) - TMath::Power(pxTot,2) - TMath::Power(pyTot,2) - TMath::Power(pzTot,2);
313 Double_t ptTot = TMath::Sqrt( TMath::Power(pxTot,2) + TMath::Power(pyTot,2));
314 Double_t qLong = TMath::Abs(eTot*qzTot - pzTot*q0Tot)/TMath::Sqrt(snn + TMath::Power(ptTot,2));
7f0257ea 315 ibin = Int_t(qLong/fP2Step);
316 //cout<<ibin<<endl;
5c33329d 317 if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
318 if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
319 if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
320 if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
7f0257ea 321 }
322 }
323 }//case 2
324 if(fAnalysisType==3) {
5c33329d 325 for(i = 1; i < gNtrack; i++) {
9d1f0df5 326 if(fAnalysisLevel == "ESD")
327 track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
5c33329d 328 Short_t charge1 = track1->Charge();
329 Double_t pX1 = track1->Px();
330 Double_t pY1 = track1->Py();
331 Double_t pZ1 = track1->Pz();
332 Double_t energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
333 TMath::Power(track1->M(),2));
7f0257ea 334 for(j = 0; j < i; j++) {
9d1f0df5 335 if(fAnalysisLevel == "ESD")
336 track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
5c33329d 337 Short_t charge2 = track2->Charge();
338 Double_t pX2 = track2->Px();
339 Double_t pY2 = track2->Py();
340 Double_t pZ2 = track2->Pz();
341 Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
342 TMath::Power(track2->M(),2));
343 Double_t eTot = energy1 + energy2;
344 Double_t pxTot = pX1 + pX2;
345 Double_t pyTot = pY1 + pY2;
346 Double_t pzTot = pZ1 + pZ2;
347 Double_t qxTot = pX1 - pX2;
348 Double_t qyTot = pY1 - pY2;
349 Double_t snn = TMath::Power(eTot,2) - TMath::Power(pxTot,2) - TMath::Power(pyTot,2) - TMath::Power(pzTot,2);
350 Double_t ptTot = TMath::Sqrt( TMath::Power(pxTot,2) + TMath::Power(pyTot,2));
351 Double_t qOut = TMath::Sqrt(snn/(snn + TMath::Power(ptTot,2))) * TMath::Abs(pxTot*qxTot + pyTot*qyTot)/ptTot;
7f0257ea 352 ibin = Int_t(qOut/fP2Step);
353 //cout<<ibin<<endl;
5c33329d 354 if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
355 if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
356 if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
357 if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
7f0257ea 358 }
359 }
360 }//case 3
361 if(fAnalysisType==4) {
5c33329d 362 for(i = 1; i < gNtrack; i++) {
9d1f0df5 363 if(fAnalysisLevel == "ESD")
364 track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
5c33329d 365 Short_t charge1 = track1->Charge();
366 Double_t pX1 = track1->Px();
367 Double_t pY1 = track1->Py();
368 //Double_t energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
369 //TMath::Power(track1->M(),2));
7f0257ea 370 for(j = 0; j < i; j++) {
9d1f0df5 371 if(fAnalysisLevel == "ESD")
372 track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
5c33329d 373 Short_t charge2 = track2->Charge();
374 Double_t pX2 = track2->Px();
375 Double_t pY2 = track2->Py();
376 //Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
377 //TMath::Power(track2->M(),2));
378 //Double_t eTot = energy1 + energy2;
379 Double_t pxTot = pX1 + pX2;
380 Double_t pyTot = pY1 + pY2;
381 Double_t qxTot = pX1 - pX2;
382 Double_t qyTot = pY1 - pY2;
383 Double_t ptTot = TMath::Sqrt( TMath::Power(pxTot,2) + TMath::Power(pyTot,2));
7f0257ea 384 Double_t qSide = TMath::Abs(pxTot*qyTot - pyTot*qxTot)/ptTot;
385 ibin = Int_t(qSide/fP2Step);
386 //cout<<ibin<<endl;
5c33329d 387 if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
388 if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
389 if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
390 if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
7f0257ea 391 }
392 }
393 }//case 4
394 if(fAnalysisType==5) {
5c33329d 395 for(i = 1; i < gNtrack; i++) {
9d1f0df5 396 if(fAnalysisLevel == "ESD")
397 track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
5c33329d 398 Short_t charge1 = track1->Charge();
399 Double_t pX1 = track1->Px();
400 Double_t pY1 = track1->Py();
401 Double_t pZ1 = track1->Pz();
402 Double_t energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
403 TMath::Power(track1->M(),2));
7f0257ea 404 for(j = 0; j < i; j++) {
9d1f0df5 405 if(fAnalysisLevel == "ESD")
406 track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
5c33329d 407 Short_t charge2 = track2->Charge();
408 Double_t pX2 = track2->Px();
409 Double_t pY2 = track2->Py();
410 Double_t pZ2 = track2->Pz();
411 Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
412 TMath::Power(track2->M(),2));
413 Double_t q0Tot = energy1 - energy2;
414 Double_t qxTot = pX1 - pX2;
415 Double_t qyTot = pY1 - pY2;
416 Double_t qzTot = pZ1 - pZ2;
417 Double_t qInv = TMath::Sqrt(TMath::Abs(-TMath::Power(q0Tot,2) +TMath::Power(qxTot,2) +TMath::Power(qyTot,2) +TMath::Power(qzTot,2)));
7f0257ea 418 ibin = Int_t(qInv/fP2Step);
419 //cout<<ibin<<endl;
5c33329d 420 if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
421 if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
422 if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
423 if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
7f0257ea 424 }
425 }
426 }//case 5
427 if(fAnalysisType==6) {
5c33329d 428 for(i = 1; i < gNtrack; i++) {
9d1f0df5 429 if(fAnalysisLevel == "ESD")
430 track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
5c33329d 431 Short_t charge1 = track1->Charge();
432 Double_t phi1 = track1->Phi();
7f0257ea 433 for(j = 0; j < i; j++) {
9d1f0df5 434 if(fAnalysisLevel == "ESD")
435 track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
5c33329d 436 Short_t charge2 = track2->Charge();
437 Double_t phi2 = track2->Phi();
7f0257ea 438 Double_t dphi = TMath::Abs(phi1 - phi2);
439 ibin = Int_t(dphi/fP2Step);
5c33329d 440 if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
441 if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
442 if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
443 if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
7f0257ea 444 }
445 }
446 }//case 6
5c33329d 447
448 /*for(Int_t i = 0; i < fNumberOfBins; i++)
449 Printf("bin: %d - Npp: %lf - Nnn: %lf - Nnp: %lf - Npn: %lf",i,fNpp[i],fNnn[i],fNpn[i],fNpn[i]);*/
6c178944 450}
451
7f0257ea 452//____________________________________________________________________//
453Double_t AliBalance::GetBalance(Int_t p2) {
6c178944 454 // Returns the value of the balance function in bin p2
455 fB[p2] = 0.5*(((fNpn[p2] - 2.0*fNnn[p2])/fNn) + ((fNpn[p2] - 2.0*fNpp[p2])/fNp))/fP2Step;
7f0257ea 456
6c178944 457 return fB[p2];
458}
459
7f0257ea 460//____________________________________________________________________//
461Double_t AliBalance::GetError(Int_t p2) {
6c178944 462 // Returns the error on the BF value for bin p2
5c33329d 463 ferror[p2] = TMath::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])*TMath::Power((0.5/Double_t(fNp) + 0.5/Double_t(fNn)),2))/fP2Step;
6c178944 464
465 return ferror[p2];
466}
467
7f0257ea 468//____________________________________________________________________//
469void AliBalance::PrintResults() {
6c178944 470 // Prints the results
471 Double_t x[MAXIMUM_NUMBER_OF_STEPS];
472 Double_t fSumXi = 0.0, fSumBi = 0.0, fSumBiXi = 0.0;
473 Double_t fSumBiXi2 = 0.0, fSumBi2Xi2 = 0.0;
474 Double_t fSumDeltaBi2 = 0.0, fSumXi2DeltaBi2 = 0.0;
475 Double_t deltaBalP2 = 0.0, integral = 0.0;
476 Double_t deltaErrorNew = 0.0;
477
478 cout<<"=================================================="<<endl;
7f0257ea 479 for(Int_t i = 0; i < fNumberOfBins; i++) {
5c33329d 480 //x[i] = fP2Start + fP2Step*i + fP2Step/2;
481 x[i] = fP2Step*i + fP2Step/2;
7f0257ea 482 cout<<"B: "<<fB[i]<<"\t Error: "<<ferror[i]<<"\t bin: "<<x[i]<<endl;
483 }
6c178944 484 cout<<"=================================================="<<endl;
7f0257ea 485 for(Int_t i = 1; i < fNumberOfBins; i++) {
486 fSumXi += x[i];
487 fSumBi += fB[i];
488 fSumBiXi += fB[i]*x[i];
5c33329d 489 fSumBiXi2 += fB[i]*TMath::Power(x[i],2);
490 fSumBi2Xi2 += TMath::Power(fB[i],2)*TMath::Power(x[i],2);
491 fSumDeltaBi2 += TMath::Power(ferror[i],2);
492 fSumXi2DeltaBi2 += TMath::Power(x[i],2) * TMath::Power(ferror[i],2);
7f0257ea 493
5c33329d 494 deltaBalP2 += fP2Step*TMath::Power(ferror[i],2);
7f0257ea 495 integral += fP2Step*fB[i];
496 }
5c33329d 497 for(Int_t i = 1; i < fNumberOfBins; i++) deltaErrorNew += ferror[i]*(x[i]*fSumBi - fSumBiXi)/TMath::Power(fSumBi,2);
7f0257ea 498
5c33329d 499 Double_t integralError = TMath::Sqrt(deltaBalP2);
6c178944 500
7f0257ea 501 Double_t delta = fSumBiXi / fSumBi;
5c33329d 502 Double_t deltaError = (fSumBiXi / fSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(fSumXi2DeltaBi2)/fSumBiXi),2) + TMath::Power((fSumDeltaBi2/fSumBi),2) );
6c178944 503
504 cout<<"Analyzed events: "<<fAnalyzedEvents<<endl;
505 cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
506 cout<<"New error: "<<deltaErrorNew<<endl;
507 cout<<"Interval: "<<integral<<"\t Error: "<<integralError<<endl;
508 cout<<"=================================================="<<endl;
509}
510
7f0257ea 511//____________________________________________________________________//
512TGraphErrors *AliBalance::DrawBalance() {
6c178944 513 // Draws the BF
514 Double_t x[MAXIMUM_NUMBER_OF_STEPS];
515 Double_t xer[MAXIMUM_NUMBER_OF_STEPS];
516 Double_t b[MAXIMUM_NUMBER_OF_STEPS];
517 Double_t ber[MAXIMUM_NUMBER_OF_STEPS];
518
7f0257ea 519 if((fNp == 0)||(fNn == 0)) {
520 cout<<"Couldn't find any particles in the analyzed interval!!!"<<endl;
521 cout<<"Aborting....."<<endl;
522 abort();
523 }
524
525 for(Int_t i = 0; i < fNumberOfBins; i++) {
526 b[i] = GetBalance(i);
527 ber[i] = GetError(i);
5c33329d 528 //x[i] = fP2Start + fP2Step*i + fP2Step/2;
529 x[i] = fP2Step*i + fP2Step/2;
7f0257ea 530 xer[i] = 0.0;
531 }
532
533 TGraphErrors *gr = new TGraphErrors(fNumberOfBins,x,b,xer,ber);
534 gr->SetMarkerStyle(25);
6c178944 535 gr->GetXaxis()->SetTitleColor(1);
7f0257ea 536 if(fAnalysisType==0) {
537 gr->GetXaxis()->SetTitle("#Delta y");
538 gr->GetYaxis()->SetTitle("B(#Delta y)");
539 }
540 if(fAnalysisType==1) {
541 gr->GetXaxis()->SetTitle("#Delta #eta");
542 gr->GetYaxis()->SetTitle("B(#Delta #eta)");
543 }
544 if(fAnalysisType==2) {
545 gr->GetXaxis()->SetTitle("Q_{long} [GeV]");
546 gr->GetYaxis()->SetTitle("B(Q_{long})");
547 }
548 if(fAnalysisType==3) {
549 gr->GetXaxis()->SetTitle("Q_{out} [GeV]");
550 gr->GetYaxis()->SetTitle("B(Q_{out})");
551 }
552 if(fAnalysisType==4) {
553 gr->GetXaxis()->SetTitle("Q_{side} [GeV]");
554 gr->GetYaxis()->SetTitle("B(Q_{side})");
555 }
556 if(fAnalysisType==5) {
557 gr->GetXaxis()->SetTitle("Q_{inv} [GeV]");
558 gr->GetYaxis()->SetTitle("B(Q_{inv})");
559 }
560 if(fAnalysisType==6) {
561 gr->GetXaxis()->SetTitle("#Delta #phi");
562 gr->GetYaxis()->SetTitle("B(#Delta #phi)");
563 }
6c178944 564
565 return gr;
566}
898d5b28 567
568//____________________________________________________________________//
569void AliBalance::Merge(AliBalance *b) {
570 //Merging function to be used for proof and grid
571 fNp += b->GetNp();
572 fNn += b->GetNn();
573 for(Int_t i = 0; i < fNumberOfBins; i++) {
574 fNnn[i] += b->GetNnn(i);
575 fNpp[i] += b->GetNpp(i);
576 fNpn[i] += b->GetNpn(i);
577 }
578}