]> 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
Fix for segfault if methods called in wrong order.
[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));
c511520f 338
339 if(track) {
340 Short_t charge = track->Charge();
341 if(charge > 0) fNp += 1.;
342 if(charge < 0) fNn += 1.;
343 }
344 else continue;
7f0257ea 345 }
5c33329d 346 //Printf("Np: %lf - Nn: %lf",fNp,fNn);
6c178944 347
348 //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
7f0257ea 349 if(fAnalysisType==0) {
5c33329d 350 for(i = 1; i < gNtrack; i++) {
9d1f0df5 351 if(fAnalysisLevel == "ESD")
352 track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
cd54a838 353 else if(fAnalysisLevel == "AOD")
354 track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
355 else if(fAnalysisLevel == "MC")
356 track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
c511520f 357
358 Short_t charge1 = 0;
359 Double_t pZ1 = 0., energy1 = 0.;
360 if(track1) {
361 charge1 = track1->Charge();
362 pZ1 = track1->Pz();
363 energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
364 TMath::Power(track1->M(),2));
365 }
366 else continue;
367
7f0257ea 368 for(j = 0; j < i; j++) {
9d1f0df5 369 if(fAnalysisLevel == "ESD")
370 track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
cd54a838 371 else if(fAnalysisLevel == "AOD")
372 track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
373 else if(fAnalysisLevel == "MC")
374 track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
c511520f 375
376 if(track2) {
377 Short_t charge2 = track2->Charge();
378 Double_t pZ2 = track2->Pz();
379 Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
380 TMath::Power(track2->M(),2));
381
382 Double_t rap1 = 0.5*log((energy1 + pZ1)/(energy1 - pZ1));
383 Double_t rap2 = 0.5*log((energy2 + pZ2)/(energy2 - pZ2));
384 Double_t dy = TMath::Abs(rap1 - rap2);
385 ibin = Int_t(dy/fP2Step);
386 if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
387 if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
388 if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
389 if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
390 }
391 else continue;
7f0257ea 392 }
393 }
394 }//case 0
395 if(fAnalysisType==1) {
5c33329d 396 for(i = 1; i < gNtrack; i++) {
9d1f0df5 397 if(fAnalysisLevel == "ESD")
398 track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
cd54a838 399 else if(fAnalysisLevel == "AOD")
400 track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
401 else if(fAnalysisLevel == "MC")
402 track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
c511520f 403
404 Short_t charge1 = 0;
405 Double_t pZ1 = 0., p1 = 0.;
406 if(track1) {
407 charge1 = track1->Charge();
408 pZ1 = track1->Pz();
409 p1 = track1->P();
410 }
411 else continue;
412
7f0257ea 413 for(j = 0; j < i; j++) {
9d1f0df5 414 if(fAnalysisLevel == "ESD")
415 track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
cd54a838 416 else if(fAnalysisLevel == "AOD")
417 track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
418 else if(fAnalysisLevel == "MC")
419 track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
c511520f 420
421 if(track2) {
422 Short_t charge2 = track2->Charge();
423 Double_t pZ2 = track2->Pz();
424 Double_t p2 = track2->P();
425 Double_t eta1 = 0.5*log((p1 + pZ1)/(p1 - pZ1));
426 Double_t eta2 = 0.5*log((p2 + pZ2)/(p2 - pZ2));
427 Double_t deta = TMath::Abs(eta1 - eta2);
428 ibin = Int_t(deta/fP2Step);
429
430 if((charge1 > 0.)&&(charge2 > 0.)) fNpp[ibin] += 1.;
431 if((charge1 < 0.)&&(charge2 < 0.)) fNnn[ibin] += 1.;
432 if((charge1 > 0.)&&(charge2 < 0.)) fNpn[ibin] += 1.;
433 if((charge1 < 0.)&&(charge2 > 0.)) fNpn[ibin] += 1.;
434 //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]);
435 }
436 else continue;
7f0257ea 437 }
438 }
439 }//case 1
440 if(fAnalysisType==2) {
5c33329d 441 for(i = 1; i < gNtrack; i++) {
9d1f0df5 442 if(fAnalysisLevel == "ESD")
443 track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
cd54a838 444 else if(fAnalysisLevel == "AOD")
445 track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
446 else if(fAnalysisLevel == "MC")
447 track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
c511520f 448
449 Short_t charge1 = 0;
450 Double_t pX1 = 0., pY1 = 0., pZ1 = 0., energy1 = 0.;
451 if(track1) {
452 charge1 = track1->Charge();
453 pX1 = track1->Px();
454 pY1 = track1->Py();
455 pZ1 = track1->Pz();
456 energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
457 TMath::Power(track1->M(),2));
458 }
459 else continue;
460
7f0257ea 461 for(j = 0; j < i; j++) {
9d1f0df5 462 if(fAnalysisLevel == "ESD")
463 track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
cd54a838 464 else if(fAnalysisLevel == "AOD")
465 track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
466 else if(fAnalysisLevel == "MC")
467 track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
c511520f 468
469 if(track2) {
470 Short_t charge2 = track2->Charge();
471 Double_t pX2 = track2->Px();
472 Double_t pY2 = track2->Py();
473 Double_t pZ2 = track2->Pz();
474 Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
475 TMath::Power(track2->M(),2));
476 Double_t eTot = energy1 + energy2;
477 Double_t pxTot = pX1 + pX2;
478 Double_t pyTot = pY1 + pY2;
479 Double_t pzTot = pZ1 + pZ2;
480 Double_t q0Tot = energy1 - energy2;
481 Double_t qzTot = pZ1 - pZ2;
482 Double_t snn = TMath::Power(eTot,2) - TMath::Power(pxTot,2) - TMath::Power(pyTot,2) - TMath::Power(pzTot,2);
483 Double_t ptTot = TMath::Sqrt( TMath::Power(pxTot,2) + TMath::Power(pyTot,2));
484 Double_t qLong = TMath::Abs(eTot*qzTot - pzTot*q0Tot)/TMath::Sqrt(snn + TMath::Power(ptTot,2));
485 ibin = Int_t(qLong/fP2Step);
486 //cout<<ibin<<endl;
487 if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
488 if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
489 if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
490 if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
491 }
492 else continue;
7f0257ea 493 }
494 }
495 }//case 2
496 if(fAnalysisType==3) {
5c33329d 497 for(i = 1; i < gNtrack; i++) {
9d1f0df5 498 if(fAnalysisLevel == "ESD")
499 track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
cd54a838 500 else if(fAnalysisLevel == "AOD")
501 track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
502 else if(fAnalysisLevel == "MC")
503 track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
c511520f 504
505 Short_t charge1 = 0;
506 Double_t pX1 = 0., pY1 = 0., pZ1 = 0., energy1 = 0.;
507 if(track1) {
508 charge1 = track1->Charge();
509 pX1 = track1->Px();
510 pY1 = track1->Py();
511 pZ1 = track1->Pz();
512 energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
513 TMath::Power(track1->M(),2));
514 }
515 else continue;
516
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));
c511520f 524
525 if(track2) {
526 Short_t charge2 = track2->Charge();
527 Double_t pX2 = track2->Px();
528 Double_t pY2 = track2->Py();
529 Double_t pZ2 = track2->Pz();
530 Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
531 TMath::Power(track2->M(),2));
532 Double_t eTot = energy1 + energy2;
533 Double_t pxTot = pX1 + pX2;
534 Double_t pyTot = pY1 + pY2;
535 Double_t pzTot = pZ1 + pZ2;
536 Double_t qxTot = pX1 - pX2;
537 Double_t qyTot = pY1 - pY2;
538 Double_t snn = TMath::Power(eTot,2) - TMath::Power(pxTot,2) - TMath::Power(pyTot,2) - TMath::Power(pzTot,2);
539 Double_t ptTot = TMath::Sqrt( TMath::Power(pxTot,2) + TMath::Power(pyTot,2));
540 Double_t qOut = TMath::Sqrt(snn/(snn + TMath::Power(ptTot,2))) * TMath::Abs(pxTot*qxTot + pyTot*qyTot)/ptTot;
541 ibin = Int_t(qOut/fP2Step);
542 //cout<<ibin<<endl;
543 if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
544 if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
545 if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
546 if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
547 }
548 else continue;
7f0257ea 549 }
550 }
551 }//case 3
552 if(fAnalysisType==4) {
5c33329d 553 for(i = 1; i < gNtrack; i++) {
9d1f0df5 554 if(fAnalysisLevel == "ESD")
555 track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
cd54a838 556 else if(fAnalysisLevel == "AOD")
557 track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
558 else if(fAnalysisLevel == "MC")
559 track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
c511520f 560
561 Short_t charge1 = 0;
562 Double_t pX1 = 0., pY1 = 0.;
563 if(track1) {
564 charge1 = track1->Charge();
d25c18a7 565 pX1 = track1->Px();
566 pY1 = track1->Py();
c511520f 567 //Double_t energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
568 //TMath::Power(track1->M(),2));
569 }
570 else continue;
571
7f0257ea 572 for(j = 0; j < i; j++) {
9d1f0df5 573 if(fAnalysisLevel == "ESD")
574 track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
cd54a838 575 else if(fAnalysisLevel == "AOD")
576 track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
577 else if(fAnalysisLevel == "MC")
578 track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
c511520f 579
580 if(track2) {
581 Short_t charge2 = track2->Charge();
582 Double_t pX2 = track2->Px();
583 Double_t pY2 = track2->Py();
584 //Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
585 //TMath::Power(track2->M(),2));
586 //Double_t eTot = energy1 + energy2;
587 Double_t pxTot = pX1 + pX2;
588 Double_t pyTot = pY1 + pY2;
589 Double_t qxTot = pX1 - pX2;
590 Double_t qyTot = pY1 - pY2;
591 Double_t ptTot = TMath::Sqrt( TMath::Power(pxTot,2) + TMath::Power(pyTot,2));
592 Double_t qSide = TMath::Abs(pxTot*qyTot - pyTot*qxTot)/ptTot;
593 ibin = Int_t(qSide/fP2Step);
594 //cout<<ibin<<endl;
595 if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
596 if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
597 if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
598 if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
599 }
600 else continue;
7f0257ea 601 }
602 }
603 }//case 4
604 if(fAnalysisType==5) {
5c33329d 605 for(i = 1; i < gNtrack; i++) {
9d1f0df5 606 if(fAnalysisLevel == "ESD")
607 track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
cd54a838 608 else if(fAnalysisLevel == "AOD")
609 track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
610 else if(fAnalysisLevel == "MC")
611 track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
c511520f 612
613 Short_t charge1 = 0;
614 Double_t pX1 = 0., pY1 = 0., pZ1 = 0., energy1 = 0.;
615 if(track1) {
616 charge1 = track1->Charge();
617 pX1 = track1->Px();
618 pY1 = track1->Py();
619 pZ1 = track1->Pz();
620 energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
621 TMath::Power(track1->M(),2));
622 }
623 else continue;
624
7f0257ea 625 for(j = 0; j < i; j++) {
9d1f0df5 626 if(fAnalysisLevel == "ESD")
627 track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
cd54a838 628 else if(fAnalysisLevel == "AOD")
629 track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
630 else if(fAnalysisLevel == "MC")
631 track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
c511520f 632
633 if(track2) {
634 Short_t charge2 = track2->Charge();
635 Double_t pX2 = track2->Px();
636 Double_t pY2 = track2->Py();
637 Double_t pZ2 = track2->Pz();
638 Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
639 TMath::Power(track2->M(),2));
640 Double_t q0Tot = energy1 - energy2;
641 Double_t qxTot = pX1 - pX2;
642 Double_t qyTot = pY1 - pY2;
643 Double_t qzTot = pZ1 - pZ2;
644 Double_t qInv = TMath::Sqrt(TMath::Abs(-TMath::Power(q0Tot,2) +TMath::Power(qxTot,2) +TMath::Power(qyTot,2) +TMath::Power(qzTot,2)));
645 ibin = Int_t(qInv/fP2Step);
646 //cout<<ibin<<endl;
647 if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
648 if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
649 if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
650 if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
651 }
652 else continue;
7f0257ea 653 }
654 }
655 }//case 5
656 if(fAnalysisType==6) {
5c33329d 657 for(i = 1; i < gNtrack; i++) {
9d1f0df5 658 if(fAnalysisLevel == "ESD")
659 track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
cd54a838 660 else if(fAnalysisLevel == "AOD")
661 track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
662 else if(fAnalysisLevel == "MC")
663 track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
c511520f 664
665 Short_t charge1 = 0;
666 Double_t phi1 = 0.;
667 if(track1) {
668 charge1 = track1->Charge();
669 phi1 = track1->Phi();
670 }
671 else continue;
672
7f0257ea 673 for(j = 0; j < i; j++) {
9d1f0df5 674 if(fAnalysisLevel == "ESD")
675 track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
cd54a838 676 else if(fAnalysisLevel == "AOD")
677 track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
678 else if(fAnalysisLevel == "MC")
679 track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
c511520f 680
681 if(track2) {
682 Short_t charge2 = track2->Charge();
683 Double_t phi2 = track2->Phi();
684 Double_t dphi = TMath::Abs(phi1 - phi2);
685 ibin = Int_t(dphi/fP2Step);
686 if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
687 if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
688 if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
689 if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
690 }
691 else continue;
7f0257ea 692 }
693 }
694 }//case 6
5c33329d 695
696 /*for(Int_t i = 0; i < fNumberOfBins; i++)
697 Printf("bin: %d - Npp: %lf - Nnn: %lf - Nnp: %lf - Npn: %lf",i,fNpp[i],fNnn[i],fNpn[i],fNpn[i]);*/
6c178944 698}
699
1d8792e7 700//____________________________________________________________________//
701TH1F *AliBalance::GetHistNnn() {
702 //Return the control histogram of the -- component
703 for(Int_t iBin = 1; iBin <= fNumberOfBins; iBin++)
704 fHistfNnn->SetBinContent(iBin,GetNnn(iBin-1));
705
706 return fHistfNnn;
707}
708
709//____________________________________________________________________//
710TH1F *AliBalance::GetHistNpp() {
711 //Return the control histogram of the ++ component
712 for(Int_t iBin = 1; iBin <= fNumberOfBins; iBin++)
713 fHistfNpp->SetBinContent(iBin,GetNpp(iBin-1));
714
715 return fHistfNpp;
716}
717
718//____________________________________________________________________//
719TH1F *AliBalance::GetHistNpn() {
720 //Return the control histogram of the +- component
721 for(Int_t iBin = 1; iBin <= fNumberOfBins; iBin++)
722 fHistfNpn->SetBinContent(iBin,GetNpn(iBin-1));
723
724 return fHistfNpn;
725}
726
7f0257ea 727//____________________________________________________________________//
728Double_t AliBalance::GetBalance(Int_t p2) {
6c178944 729 // Returns the value of the balance function in bin p2
730 fB[p2] = 0.5*(((fNpn[p2] - 2.0*fNnn[p2])/fNn) + ((fNpn[p2] - 2.0*fNpp[p2])/fNp))/fP2Step;
7f0257ea 731
6c178944 732 return fB[p2];
733}
734
7f0257ea 735//____________________________________________________________________//
736Double_t AliBalance::GetError(Int_t p2) {
6c178944 737 // Returns the error on the BF value for bin p2
5c33329d 738 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 739
740 return ferror[p2];
741}
742
7f0257ea 743//____________________________________________________________________//
744void AliBalance::PrintResults() {
6c178944 745 // Prints the results
746 Double_t x[MAXIMUM_NUMBER_OF_STEPS];
747 Double_t fSumXi = 0.0, fSumBi = 0.0, fSumBiXi = 0.0;
748 Double_t fSumBiXi2 = 0.0, fSumBi2Xi2 = 0.0;
749 Double_t fSumDeltaBi2 = 0.0, fSumXi2DeltaBi2 = 0.0;
750 Double_t deltaBalP2 = 0.0, integral = 0.0;
751 Double_t deltaErrorNew = 0.0;
752
753 cout<<"=================================================="<<endl;
7f0257ea 754 for(Int_t i = 0; i < fNumberOfBins; i++) {
5c33329d 755 //x[i] = fP2Start + fP2Step*i + fP2Step/2;
756 x[i] = fP2Step*i + fP2Step/2;
7f0257ea 757 cout<<"B: "<<fB[i]<<"\t Error: "<<ferror[i]<<"\t bin: "<<x[i]<<endl;
758 }
6c178944 759 cout<<"=================================================="<<endl;
7f0257ea 760 for(Int_t i = 1; i < fNumberOfBins; i++) {
761 fSumXi += x[i];
762 fSumBi += fB[i];
763 fSumBiXi += fB[i]*x[i];
5c33329d 764 fSumBiXi2 += fB[i]*TMath::Power(x[i],2);
765 fSumBi2Xi2 += TMath::Power(fB[i],2)*TMath::Power(x[i],2);
766 fSumDeltaBi2 += TMath::Power(ferror[i],2);
767 fSumXi2DeltaBi2 += TMath::Power(x[i],2) * TMath::Power(ferror[i],2);
7f0257ea 768
5c33329d 769 deltaBalP2 += fP2Step*TMath::Power(ferror[i],2);
7f0257ea 770 integral += fP2Step*fB[i];
771 }
5c33329d 772 for(Int_t i = 1; i < fNumberOfBins; i++) deltaErrorNew += ferror[i]*(x[i]*fSumBi - fSumBiXi)/TMath::Power(fSumBi,2);
7f0257ea 773
5c33329d 774 Double_t integralError = TMath::Sqrt(deltaBalP2);
6c178944 775
7f0257ea 776 Double_t delta = fSumBiXi / fSumBi;
5c33329d 777 Double_t deltaError = (fSumBiXi / fSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(fSumXi2DeltaBi2)/fSumBiXi),2) + TMath::Power((fSumDeltaBi2/fSumBi),2) );
6c178944 778
779 cout<<"Analyzed events: "<<fAnalyzedEvents<<endl;
780 cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
781 cout<<"New error: "<<deltaErrorNew<<endl;
782 cout<<"Interval: "<<integral<<"\t Error: "<<integralError<<endl;
783 cout<<"=================================================="<<endl;
784}
785
7f0257ea 786//____________________________________________________________________//
787TGraphErrors *AliBalance::DrawBalance() {
6c178944 788 // Draws the BF
789 Double_t x[MAXIMUM_NUMBER_OF_STEPS];
790 Double_t xer[MAXIMUM_NUMBER_OF_STEPS];
791 Double_t b[MAXIMUM_NUMBER_OF_STEPS];
792 Double_t ber[MAXIMUM_NUMBER_OF_STEPS];
793
7f0257ea 794 if((fNp == 0)||(fNn == 0)) {
795 cout<<"Couldn't find any particles in the analyzed interval!!!"<<endl;
796 cout<<"Aborting....."<<endl;
797 abort();
798 }
799
800 for(Int_t i = 0; i < fNumberOfBins; i++) {
801 b[i] = GetBalance(i);
802 ber[i] = GetError(i);
5c33329d 803 //x[i] = fP2Start + fP2Step*i + fP2Step/2;
804 x[i] = fP2Step*i + fP2Step/2;
7f0257ea 805 xer[i] = 0.0;
806 }
807
808 TGraphErrors *gr = new TGraphErrors(fNumberOfBins,x,b,xer,ber);
809 gr->SetMarkerStyle(25);
6c178944 810 gr->GetXaxis()->SetTitleColor(1);
7f0257ea 811 if(fAnalysisType==0) {
812 gr->GetXaxis()->SetTitle("#Delta y");
813 gr->GetYaxis()->SetTitle("B(#Delta y)");
814 }
815 if(fAnalysisType==1) {
816 gr->GetXaxis()->SetTitle("#Delta #eta");
817 gr->GetYaxis()->SetTitle("B(#Delta #eta)");
818 }
819 if(fAnalysisType==2) {
09bb7bf4 820 gr->GetXaxis()->SetTitle("q_{long} (GeV/c)");
821 gr->GetYaxis()->SetTitle("B(q_{long}) [(GeV/c)^{-1}]");
7f0257ea 822 }
823 if(fAnalysisType==3) {
09bb7bf4 824 gr->GetXaxis()->SetTitle("q_{out} (GeV/c)");
825 gr->GetYaxis()->SetTitle("B(q_{out}) [(GeV/c)^{-1}]");
7f0257ea 826 }
827 if(fAnalysisType==4) {
09bb7bf4 828 gr->GetXaxis()->SetTitle("q_{side} (GeV/c)");
829 gr->GetYaxis()->SetTitle("B(q_{side}) [(GeV/c)^{-1}]");
7f0257ea 830 }
831 if(fAnalysisType==5) {
09bb7bf4 832 gr->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
833 gr->GetYaxis()->SetTitle("B(q_{inv}) [(GeV/c)^{-1}]");
7f0257ea 834 }
835 if(fAnalysisType==6) {
836 gr->GetXaxis()->SetTitle("#Delta #phi");
837 gr->GetYaxis()->SetTitle("B(#Delta #phi)");
838 }
6c178944 839
840 return gr;
841}
898d5b28 842
843//____________________________________________________________________//
844void AliBalance::Merge(AliBalance *b) {
845 //Merging function to be used for proof and grid
846 fNp += b->GetNp();
847 fNn += b->GetNn();
848 for(Int_t i = 0; i < fNumberOfBins; i++) {
849 fNnn[i] += b->GetNnn(i);
850 fNpp[i] += b->GetNpp(i);
851 fNpn[i] += b->GetNpn(i);
852 }
853}