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