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