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