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