]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ANALYSIS/AliBalance.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Updated list of classes
[u/mrichter/AliRoot.git] / ANALYSIS / AliBalance.cxx
CommitLineData
cd286c84 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
cd286c84 23//ROOT
24#include <Riostream.h>
cd286c84 25#include <TMath.h>
26#include <TLorentzVector.h>
36b803af 27#include <TGraphErrors.h>
cd286c84 28
29#include "AliBalance.h"
30
31ClassImp(AliBalance)
32
33//----------------------------------------//
34AliBalance::AliBalance()
35{
93809d0d 36 // Default constructor
cd286c84 37 for(Int_t i = 0; i < MAXIMUM_NUMBER_OF_STEPS; i++)
38 {
39 fNpp[i] = .0;
40 fNnn[i] = .0;
41 fNpn[i] = .0;
42 fB[i] = 0.0;
43 ferror[i] = 0.0;
44 }
94d49dd3 45 fNp = 0.0;
46 fNn = 0.0;
93809d0d 47 fP2Start = 0.0;
48 fP2Stop = 0.0;
49 fP2Step = 0.0;
94d49dd3 50 fAnalysisType = 0;
51 fNumberOfBins = 0;
52 fNtrack = 0;
cd286c84 53 fAnalyzedEvents = 0;
54}
55
56//----------------------------------------//
93809d0d 57AliBalance::AliBalance(Double_t p2Start, Double_t p2Stop, Int_t p2Bins)
cd286c84 58{
93809d0d 59 // Constructor
60 for(Int_t i = 0; i < MAXIMUM_NUMBER_OF_STEPS; i++)
61 {
62 fNpp[i] = .0;
63 fNnn[i] = .0;
64 fNpn[i] = .0;
65 fB[i] = 0.0;
66 ferror[i] = 0.0;
67 }
68 fNp = 0.0;
69 fNn = 0.0;
70 fAnalysisType = 0;
71 fNtrack = 0;
72 fAnalyzedEvents = 0;
73
74 fP2Start = p2Start;
75 fP2Stop = p2Stop;
76 fNumberOfBins = p2Bins;
77 fP2Step = TMath::Abs(fP2Start - fP2Stop) / (Double_t)fNumberOfBins;
cd286c84 78}
79
80//----------------------------------------//
81AliBalance::~AliBalance()
82{
93809d0d 83 // Destructor
cd286c84 84}
85
86//----------------------------------------//
87void AliBalance::SetNumberOfBins(Int_t ibins)
88{
93809d0d 89 // Sets the number of bins for the analyzed interval
90 fNumberOfBins = ibins ;
cd286c84 91}
92
93//----------------------------------------//
93809d0d 94void AliBalance::SetInterval(Double_t p2Start, Double_t p2Stop)
cd286c84 95{
93809d0d 96 // Sets the analyzed interval.
97 // The analysis variable is set by SetAnalysisType
98 fP2Start = p2Start;
99 fP2Stop = p2Stop;
100 fP2Step = TMath::Abs(p2Start - p2Stop) / (Double_t)fNumberOfBins;
cd286c84 101}
102
103//----------------------------------------//
104void AliBalance::SetAnalysisType(Int_t iType)
105{
106 //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
107 this->fAnalysisType = iType;
108 if(fAnalysisType==0)
109 {
110 cout<<" ====================== "<<endl;
111 cout<<"||Analysis selected: y||"<<endl;
112 cout<<" ====================== "<<endl;
113 }
114 else if(fAnalysisType==1)
115 {
116 cout<<" ======================== "<<endl;
117 cout<<"||Analysis selected: eta||"<<endl;
118 cout<<" ======================== "<<endl;
119 }
120 else if(fAnalysisType==2)
121 {
122 cout<<" ========================== "<<endl;
123 cout<<"||Analysis selected: Qlong||"<<endl;
124 cout<<" ========================== "<<endl;
125 }
126 else if(fAnalysisType==3)
127 {
128 cout<<" ========================= "<<endl;
129 cout<<"||Analysis selected: Qout||"<<endl;
130 cout<<" ========================= "<<endl;
131 }
132 else if(fAnalysisType==4)
133 {
134 cout<<" ========================== "<<endl;
135 cout<<"||Analysis selected: Qside||"<<endl;
136 cout<<" ========================== "<<endl;
137 }
138 else if(fAnalysisType==5)
139 {
140 cout<<" ========================= "<<endl;
141 cout<<"||Analysis selected: Qinv||"<<endl;
142 cout<<" ========================= "<<endl;
143 }
144 else if(fAnalysisType==6)
145 {
146 cout<<" ======================== "<<endl;
147 cout<<"||Analysis selected: phi||"<<endl;
148 cout<<" ======================== "<<endl;
149 }
150 else
151 {
152 cout<<"Selection of analysis mode failed!!!"<<endl;
153 cout<<"Choices are: 0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi"<<endl;
154 abort();
155 }
156}
157
158//----------------------------------------//
159void AliBalance::SetParticles(TLorentzVector *P, Double_t *charge, Int_t dim)
160{
93809d0d 161 // Sets a new particle with given 4-momentum and charge.
162 // dim is the size of the array of charges and corresponds
163 // to the number of selected tracks.
cd286c84 164 this->fV = P;
165 this->fCharge = charge;
166 fNtrack = dim;
167}
168
169
170//----------------------------------------//
171void AliBalance::CalculateBalance()
172{
93809d0d 173 // Calculates the balance function
cd286c84 174 fAnalyzedEvents++;
175 Int_t i = 0 , j = 0;
176 Int_t ibin = 0;
177
178 for(i = 0; i < fNtrack; i++)
179 {
180 if(fCharge[i] > 0)
181 fNp += 1.;
182 if(fCharge[i] < 0)
183 fNn += 1.;
184 }
185
186 //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
187 if(fAnalysisType==0)
188 {
189 for(i = 1; i < fNtrack; i++)
190 {
191 for(j = 0; j < i; j++)
192 {
5a2d5990 193 Double_t rap1 = 0.5*log((fV[i].E() + fV[i].Pz())/(fV[i].E() - fV[i].Pz()));
194 Double_t rap2 = 0.5*log((fV[j].E() + fV[j].Pz())/(fV[j].E() - fV[j].Pz()));
cd286c84 195 Double_t dy = TMath::Abs(rap1 - rap2);
93809d0d 196 ibin = Int_t(dy/fP2Step);
cd286c84 197 if((fCharge[i] > 0)&&(fCharge[j] > 0))
198 fNpp[ibin] += 1.;
199 if((fCharge[i] < 0)&&(fCharge[j] < 0))
200 fNnn[ibin] += 1.;
201 if((fCharge[i] > 0)&&(fCharge[j] < 0))
202 fNpn[ibin] += 1.;
203 if((fCharge[i] < 0)&&(fCharge[j] > 0))
204 fNpn[ibin] += 1.;
205 }
206 }
207 }//case 0
208 if(fAnalysisType==1)
209 {
210 for(i = 1; i < fNtrack; i++)
211 {
212 for(j = 0; j < i; j++)
213 {
93809d0d 214 Double_t p1 = sqrt(pow(fV[i].Px(),2) + pow(fV[i].Py(),2) + pow(fV[i].Pz(),2));
215 Double_t p2 = sqrt(pow(fV[j].Px(),2) + pow(fV[j].Py(),2) + pow(fV[j].Pz(),2));
5a2d5990 216 Double_t eta1 = 0.5*log((p1 + fV[i].Pz())/(p1 - fV[i].Pz()));
217 Double_t eta2 = 0.5*log((p2 + fV[j].Pz())/(p2 - fV[j].Pz()));
cd286c84 218 Double_t deta = TMath::Abs(eta1 - eta2);
93809d0d 219 ibin = Int_t(deta/fP2Step);
cd286c84 220 if((fCharge[i] > 0)&&(fCharge[j] > 0))
221 fNpp[ibin] += 1.;
222 if((fCharge[i] < 0)&&(fCharge[j] < 0))
223 fNnn[ibin] += 1.;
224 if((fCharge[i] > 0)&&(fCharge[j] < 0))
225 fNpn[ibin] += 1.;
226 if((fCharge[i] < 0)&&(fCharge[j] > 0))
227 fNpn[ibin] += 1.;
228 }
229 }
230 }//case 1
231 if(fAnalysisType==2)
232 {
233 for(i = 1; i < fNtrack; i++)
234 {
235 for(j = 0; j < i; j++)
236 {
93809d0d 237 Double_t eTot = fV[i].E() + fV[j].E();
238 Double_t pxTot = fV[i].Px() + fV[j].Px();
239 Double_t pyTot = fV[i].Py() + fV[j].Py();
240 Double_t pzTot = fV[i].Pz() + fV[j].Pz();
241 Double_t q0Tot = fV[i].E() - fV[j].E();
242 Double_t qzTot = fV[i].Pz() - fV[j].Pz();
243 Double_t snn = pow(eTot,2) - pow(pxTot,2) - pow(pyTot,2) - pow(pzTot,2);
244 Double_t ptTot = sqrt( pow(pxTot,2) + pow(pyTot,2));
245 Double_t qLong = TMath::Abs(eTot*qzTot - pzTot*q0Tot)/sqrt(snn + pow(ptTot,2));
246 ibin = Int_t(qLong/fP2Step);
cd286c84 247 //cout<<ibin<<endl;
248 if((fCharge[i] > 0)&&(fCharge[j] > 0))
249 fNpp[ibin] += 1.;
250 if((fCharge[i] < 0)&&(fCharge[j] < 0))
251 fNnn[ibin] += 1.;
252 if((fCharge[i] > 0)&&(fCharge[j] < 0))
253 fNpn[ibin] += 1.;
254 if((fCharge[i] < 0)&&(fCharge[j] > 0))
255 fNpn[ibin] += 1.;
256 }
257 }
258 }//case 2
259 if(fAnalysisType==3)
260 {
261 for(i = 1; i < fNtrack; i++)
262 {
263 for(j = 0; j < i; j++)
264 {
93809d0d 265 Double_t eTot = fV[i].E() + fV[j].E();
266 Double_t pxTot = fV[i].Px() + fV[j].Px();
267 Double_t pyTot = fV[i].Py() + fV[j].Py();
268 Double_t pzTot = fV[i].Pz() + fV[j].Pz();
269 Double_t qxTot = fV[i].Px() - fV[j].Px();
270 Double_t qyTot = fV[i].Py() - fV[j].Py();
271 Double_t snn = pow(eTot,2) - pow(pxTot,2) - pow(pyTot,2) - pow(pzTot,2);
272 Double_t ptTot = sqrt( pow(pxTot,2) + pow(pyTot,2));
273 Double_t qOut = sqrt(snn/(snn + pow(ptTot,2))) * TMath::Abs(pxTot*qxTot + pyTot*qyTot)/ptTot;
274 ibin = Int_t(qOut/fP2Step);
cd286c84 275 //cout<<ibin<<endl;
276 if((fCharge[i] > 0)&&(fCharge[j] > 0))
277 fNpp[ibin] += 1.;
278 if((fCharge[i] < 0)&&(fCharge[j] < 0))
279 fNnn[ibin] += 1.;
280 if((fCharge[i] > 0)&&(fCharge[j] < 0))
281 fNpn[ibin] += 1.;
282 if((fCharge[i] < 0)&&(fCharge[j] > 0))
283 fNpn[ibin] += 1.;
284 }
285 }
286 }//case 3
287 if(fAnalysisType==4)
288 {
289 for(i = 1; i < fNtrack; i++)
290 {
291 for(j = 0; j < i; j++)
292 {
93809d0d 293 Double_t pxTot = fV[i].Px() + fV[j].Px();
294 Double_t pyTot = fV[i].Py() + fV[j].Py();
295 Double_t qxTot = fV[i].Px() - fV[j].Px();
296 Double_t qyTot = fV[i].Py() - fV[j].Py();
297 Double_t ptTot = sqrt( pow(pxTot,2) + pow(pyTot,2));
298 Double_t qSide = TMath::Abs(pxTot*qyTot - pyTot*qxTot)/ptTot;
299 ibin = Int_t(qSide/fP2Step);
cd286c84 300 //cout<<ibin<<endl;
301 if((fCharge[i] > 0)&&(fCharge[j] > 0))
302 fNpp[ibin] += 1.;
303 if((fCharge[i] < 0)&&(fCharge[j] < 0))
304 fNnn[ibin] += 1.;
305 if((fCharge[i] > 0)&&(fCharge[j] < 0))
306 fNpn[ibin] += 1.;
307 if((fCharge[i] < 0)&&(fCharge[j] > 0))
308 fNpn[ibin] += 1.;
309 }
310 }
311 }//case 4
312 if(fAnalysisType==5)
313 {
314 for(i = 1; i < fNtrack; i++)
315 {
316 for(j = 0; j < i; j++)
317 {
93809d0d 318 Double_t q0Tot = fV[i].E() - fV[j].E();
319 Double_t qxTot = fV[i].Px() - fV[j].Px();
320 Double_t qyTot = fV[i].Py() - fV[j].Py();
321 Double_t qzTot = fV[i].Pz() - fV[j].Pz();
322 Double_t qInv = sqrt(TMath::Abs(-pow(q0Tot,2) +pow(qxTot,2) +pow(qyTot,2) +pow(qzTot,2)));
323 ibin = Int_t(qInv/fP2Step);
cd286c84 324 //cout<<ibin<<endl;
325 if((fCharge[i] > 0)&&(fCharge[j] > 0))
326 fNpp[ibin] += 1.;
327 if((fCharge[i] < 0)&&(fCharge[j] < 0))
328 fNnn[ibin] += 1.;
329 if((fCharge[i] > 0)&&(fCharge[j] < 0))
330 fNpn[ibin] += 1.;
331 if((fCharge[i] < 0)&&(fCharge[j] > 0))
332 fNpn[ibin] += 1.;
333 }
334 }
335 }//case 5
336 if(fAnalysisType==6)
337 {
338 for(i = 1; i < fNtrack; i++)
339 {
340 for(j = 0; j < i; j++)
341 {
342 Double_t phi1 = TMath::ATan(fV[i].Py()/fV[i].Px())*180.0/TMath::Pi();
343 Double_t phi2 = TMath::ATan(fV[j].Py()/fV[j].Px())*180.0/TMath::Pi();
344 Double_t dphi = TMath::Abs(phi1 - phi2);
93809d0d 345 ibin = Int_t(dphi/fP2Step);
cd286c84 346 if((fCharge[i] > 0)&&(fCharge[j] > 0))
347 fNpp[ibin] += 1.;
348 if((fCharge[i] < 0)&&(fCharge[j] < 0))
349 fNnn[ibin] += 1.;
350 if((fCharge[i] > 0)&&(fCharge[j] < 0))
351 fNpn[ibin] += 1.;
352 if((fCharge[i] < 0)&&(fCharge[j] > 0))
353 fNpn[ibin] += 1.;
354 }
355 }
356 }//case 6
357}
358
359//----------------------------------------//
360Double_t AliBalance::GetBalance(Int_t p2)
361{
93809d0d 362 // Returns the value of the balance function in bin p2
363 fB[p2] = 0.5*(((fNpn[p2] - 2.0*fNnn[p2])/fNn) + ((fNpn[p2] - 2.0*fNpp[p2])/fNp))/fP2Step;
cd286c84 364
365 return fB[p2];
366}
367
368//----------------------------------------//
369Double_t AliBalance::GetError(Int_t p2)
370{
93809d0d 371 // Returns the error on the BF value for bin p2
372 ferror[p2] = sqrt( Double_t(fNpp[p2])/(Double_t(fNp)*Double_t(fNp)) + Double_t(fNnn[p2])/(Double_t(fNn)*Double_t(fNn)) + Double_t(fNpn[p2])*pow((0.5/Double_t(fNp) + 0.5/Double_t(fNn)),2))/fP2Step;
cd286c84 373
374 return ferror[p2];
375}
376
377//----------------------------------------//
378void AliBalance::PrintResults()
379{
93809d0d 380 // Prints the results
cd286c84 381 Double_t x[MAXIMUM_NUMBER_OF_STEPS];
382 Double_t fSumXi = 0.0, fSumBi = 0.0, fSumBiXi = 0.0;
383 Double_t fSumBiXi2 = 0.0, fSumBi2Xi2 = 0.0;
384 Double_t fSumDeltaBi2 = 0.0, fSumXi2DeltaBi2 = 0.0;
93809d0d 385 Double_t deltaBalP2 = 0.0, integral = 0.0;
386 Double_t deltaErrorNew = 0.0;
cd286c84 387
388 cout<<"=================================================="<<endl;
389 for(Int_t i = 0; i < fNumberOfBins; i++)
390 {
93809d0d 391 x[i] = fP2Start + fP2Step*i + fP2Step/2 ;
94d49dd3 392 cout<<"B: "<<fB[i]<<"\t Error: "<<ferror[i]<<"\t bin: "<<x[i]<<endl;
cd286c84 393 }
394 cout<<"=================================================="<<endl;
395 for(Int_t i = 1; i < fNumberOfBins; i++)
396 {
397 fSumXi += x[i];
398 fSumBi += fB[i];
399 fSumBiXi += fB[i]*x[i];
400 fSumBiXi2 += fB[i]*pow(x[i],2);
401 fSumBi2Xi2 += pow(fB[i],2)*pow(x[i],2);
402 fSumDeltaBi2 += pow(ferror[i],2) ;
403 fSumXi2DeltaBi2 += pow(x[i],2) * pow(ferror[i],2) ;
404
93809d0d 405 deltaBalP2 += fP2Step*pow(ferror[i],2) ;
406 integral += fP2Step*fB[i] ;
cd286c84 407 }
408 for(Int_t i = 1; i < fNumberOfBins; i++)
409 {
93809d0d 410 deltaErrorNew += ferror[i]*(x[i]*fSumBi - fSumBiXi)/pow(fSumBi,2);
cd286c84 411 }
93809d0d 412 Double_t integralError = sqrt(deltaBalP2) ;
cd286c84 413
93809d0d 414 Double_t delta = fSumBiXi / fSumBi ;
415 Double_t deltaError = (fSumBiXi / fSumBi) * sqrt(pow((sqrt(fSumXi2DeltaBi2)/fSumBiXi),2) + pow((fSumDeltaBi2/fSumBi),2) ) ;
cd286c84 416
417 cout<<"Analyzed events: "<<fAnalyzedEvents<<endl;
93809d0d 418 cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
419 cout<<"New error: "<<deltaErrorNew<<endl;
420 cout<<"Interval: "<<integral<<"\t Error: "<<integralError<<endl;
cd286c84 421 cout<<"=================================================="<<endl;
422}
423
36b803af 424//----------------------------------------//
425TGraphErrors *AliBalance::DrawBalance()
426{
427 // Draws the BF
428 Double_t x[MAXIMUM_NUMBER_OF_STEPS];
429 Double_t xer[MAXIMUM_NUMBER_OF_STEPS];
430 Double_t b[MAXIMUM_NUMBER_OF_STEPS];
431 Double_t ber[MAXIMUM_NUMBER_OF_STEPS];
432
433 if((fNp == 0)||(fNn == 0))
434 {
435 cout<<"Couldn't find any particles in the analyzed interval!!!"<<endl;
436 cout<<"Aborting....."<<endl;
437 abort();
438 }
439
440 for(Int_t i = 0; i < fNumberOfBins; i++)
441 {
442 b[i] = GetBalance(i);
443 ber[i] = GetError(i);
444 x[i] = fP2Start + fP2Step*i + fP2Step/2 ;
445 xer[i] = 0.0;
446 }
447
448 TGraphErrors *gr = new TGraphErrors(fNumberOfBins,x,b,xer,ber) ;
449 gr->SetMarkerStyle(25) ;
450 gr->GetXaxis()->SetTitleColor(1);
451 if(fAnalysisType==0)
452 {
453 gr->GetXaxis()->SetTitle("#Delta y");
454 gr->GetYaxis()->SetTitle("B(#Delta y)");
455 }
456 if(fAnalysisType==1)
457 {
458 gr->GetXaxis()->SetTitle("#Delta #eta");
459 gr->GetYaxis()->SetTitle("B(#Delta #eta)");
460 }
461 if(fAnalysisType==2)
462 {
463 gr->GetXaxis()->SetTitle("Q_{long} [GeV]");
464 gr->GetYaxis()->SetTitle("B(Q_{long})");
465 }
466 if(fAnalysisType==3)
467 {
468 gr->GetXaxis()->SetTitle("Q_{out} [GeV]");
469 gr->GetYaxis()->SetTitle("B(Q_{out})");
470 }
471 if(fAnalysisType==4)
472 {
473 gr->GetXaxis()->SetTitle("Q_{side} [GeV]");
474 gr->GetYaxis()->SetTitle("B(Q_{side})");
475 }
476 if(fAnalysisType==5)
477 {
478 gr->GetXaxis()->SetTitle("Q_{inv} [GeV]");
479 gr->GetYaxis()->SetTitle("B(Q_{inv})");
480 }
481 if(fAnalysisType==6)
482 {
483 gr->GetXaxis()->SetTitle("#Delta #phi");
484 gr->GetYaxis()->SetTitle("B(#Delta #phi)");
485 }
486 gr->Draw("AP") ;
487
488 return gr;
489}