More mods to suit ROOT head
[u/mrichter/AliRoot.git] / PWG2 / EBYE / AliBalance.cxx
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>
26 #include <TAxis.h>
27 #include <TLorentzVector.h>
28 #include <TGraphErrors.h>
29
30 #include "AliBalance.h"
31
32 ClassImp(AliBalance)
33
34 //____________________________________________________________________//
35 AliBalance::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) {
40   // Default constructor
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   } 
48 }
49
50 //____________________________________________________________________//
51 AliBalance::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) {
58   // Constructor
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   } 
66 }
67
68 //____________________________________________________________________//
69 AliBalance::AliBalance(const AliBalance& balance):
70   TObject(balance),
71   fNtrack(balance.fNtrack),
72   fNumberOfBins(balance.fNumberOfBins),
73   fAnalysisType(balance.fAnalysisType),
74   fAnalyzedEvents(balance.fAnalyzedEvents),
75   fP2Start(balance.fP2Start),
76   fP2Stop(balance.fP2Stop),
77   fP2Step(balance.fP2Step),
78   fNn(balance.fNn),
79   fNp(balance.fNp),
80   fCharge(0),
81   fV(0)
82 {
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 //____________________________________________________________________//
98 AliBalance::~AliBalance() {
99   // Destructor
100   delete fV;
101   delete fCharge;
102 }
103
104 //____________________________________________________________________//
105 void 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 //____________________________________________________________________//
111 void 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 //____________________________________________________________________//
117 void 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 //____________________________________________________________________//
123 void AliBalance::SetNumberOfBins(Int_t ibins) {
124   // Sets the number of bins for the analyzed interval
125   fNumberOfBins = ibins;
126 }
127
128 //____________________________________________________________________//
129 void AliBalance::SetInterval(Double_t p2Start, Double_t p2Stop) {
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
137 //____________________________________________________________________//
138 void AliBalance::SetAnalysisType(Int_t iType) {
139   //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
140   this->fAnalysisType = iType; 
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   }
181 }
182
183 //____________________________________________________________________//
184 void AliBalance::SetParticles(TLorentzVector *P, Double_t *charge, Int_t dim) {
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
194 //____________________________________________________________________//
195 void AliBalance::CalculateBalance() {
196   // Calculates the balance function
197   fAnalyzedEvents++;
198   Int_t i = 0 , j = 0;
199   Int_t ibin = 0;
200
201   for(i = 0; i < fNtrack; i++) {
202     if(fCharge[i] > 0) fNp += 1.;
203     if(fCharge[i] < 0) fNn += 1.;
204   }
205
206   //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
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
328 }
329
330 //____________________________________________________________________//
331 Double_t AliBalance::GetBalance(Int_t p2) {
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;
334   
335   return fB[p2];
336 }
337     
338 //____________________________________________________________________//
339 Double_t AliBalance::GetError(Int_t p2) {
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
346 //____________________________________________________________________//
347 void AliBalance::PrintResults() {
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;
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   } 
361   cout<<"=================================================="<<endl;
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);
377   
378   Double_t delta = fSumBiXi / fSumBi;
379   Double_t deltaError = (fSumBiXi / fSumBi) * sqrt(pow((sqrt(fSumXi2DeltaBi2)/fSumBiXi),2) + pow((fSumDeltaBi2/fSumBi),2) );
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   
388 //____________________________________________________________________//
389 TGraphErrors *AliBalance::DrawBalance() {
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
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);
411   gr->GetXaxis()->SetTitleColor(1);
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   }
440
441   return gr;
442 }