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 <TObjArray.h>
  29 #include <TGraphErrors.h>
  30 #include <TString.h>
  31 #include <TH1F.h>
  32
  33 #include "AliVParticle.h"
  34 #include "AliMCParticle.h"
  35 #include "AliESDtrack.h"
  36 #include "AliAODTrack.h"
  37
  38 #include "AliBalance.h"
  39
  40 ClassImp(AliBalance)
  41
  42 //____________________________________________________________________//
  43 AliBalance::AliBalance() :
  44   TObject(),
  45   fAnalysisLevel("ESD"), fNumberOfBins(0),
  46   fAnalysisType(0), fAnalyzedEvents(0), fP2Start(0),
  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)) {
  54   // Default constructor
  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   }
  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   }
 102 }
 103
 104 //____________________________________________________________________//
 105 AliBalance::AliBalance(Double_t p2Start, Double_t p2Stop, Int_t p2Bins) :
 106   TObject(), fAnalysisLevel("ESD"),
 107   fNumberOfBins(p2Bins), fAnalysisType(0),
 108   fAnalyzedEvents(0), fP2Start(p2Start), fP2Stop(p2Stop),
 109   fP2Step(TMath::Abs(fP2Start - fP2Stop) / (Double_t)fNumberOfBins),
 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)) {
 117   // Constructor
 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   }
 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   }
 165 }
 166
 167 //____________________________________________________________________//
 168 AliBalance::AliBalance(const AliBalance& balance):
 169   TObject(balance), fAnalysisLevel(balance.fAnalysisLevel),
 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),
 177   fNp(balance.fNp),
 178   fHistfNnn(balance.fHistfNnn),
 179   fHistfNpp(balance.fHistfNpp),
 180   fHistfNpn(balance.fHistfNpn) {
 181   //copy constructor
 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   }
 189 }
 190
 191 //____________________________________________________________________//
 192 AliBalance::~AliBalance() {
 193   // Destructor
 194   if(fHistfNnn) delete fHistfNnn;
 195   if(fHistfNpp) delete fHistfNpp;
 196   if(fHistfNpn) delete fHistfNpn;
 197 }
 198
 199 //____________________________________________________________________//
 200 void 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 //____________________________________________________________________//
 206 void 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 //____________________________________________________________________//
 212 void 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 //____________________________________________________________________//
 218 void AliBalance::SetNumberOfBins(Int_t ibins) {
 219   // Sets the number of bins for the analyzed interval
 220   fNumberOfBins = ibins;
 221 }
 222
 223 //____________________________________________________________________//
 224 void AliBalance::SetInterval(Double_t p2Start, Double_t p2Stop) {
 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
 232 //____________________________________________________________________//
 233 void AliBalance::SetAnalysisType(Int_t iType) {
 234   //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
 235   this->fAnalysisType = iType;
 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   }
 276 }
 277
 278 //____________________________________________________________________//
 279 void AliBalance::PrintAnalysisSettings() {
 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   }
 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("======================================");
 316 }
 317
 318 //____________________________________________________________________//
 319 void AliBalance::CalculateBalance(TObjArray *gTrackArray) {
 320   // Calculates the balance function
 321   fAnalyzedEvents++;
 322   Int_t i = 0 , j = 0;
 323   Int_t ibin = 0;
 324
 325   AliVParticle* track = 0;
 326   AliVParticle* track1 = 0;
 327   AliVParticle* track2 = 0;
 328
 329   //Printf("(AliBalance) Number of tracks: %d",gTrackArray->GetEntries());
 330   Int_t gNtrack = gTrackArray->GetEntries();
 331   for(i = 0; i < gNtrack; i++) {
 332     if(fAnalysisLevel == "ESD")
 333       track = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
 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));
 338     Short_t charge = track->Charge();
 339     if(charge > 0) fNp += 1.;
 340     if(charge < 0) fNn += 1.;
 341   }
 342   //Printf("Np: %lf - Nn: %lf",fNp,fNn);
 343
 344   //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
 345   if(fAnalysisType==0) {
 346     for(i = 1; i < gNtrack; i++) {
 347       if(fAnalysisLevel == "ESD")
 348         track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
 349       else if(fAnalysisLevel == "AOD")
 350         track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
 351       else if(fAnalysisLevel == "MC")
 352         track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
 353       Short_t charge1 = track1->Charge();
 354       Double_t pZ1 = track1->Pz();
 355       Double_t energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
 356                                      TMath::Power(track1->M(),2));
 357       for(j = 0; j < i; j++) {
 358         if(fAnalysisLevel == "ESD")
 359           track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
 360         else if(fAnalysisLevel == "AOD")
 361           track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
 362         else if(fAnalysisLevel == "MC")
 363           track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
 364         Short_t charge2 = track2->Charge();
 365         Double_t pZ2 = track2->Pz();
 366         Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
 367                                        TMath::Power(track2->M(),2));
 368
 369         Double_t rap1 = 0.5*log((energy1 + pZ1)/(energy1 - pZ1));
 370         Double_t rap2 = 0.5*log((energy2 + pZ2)/(energy2 - pZ2));
 371         Double_t dy = TMath::Abs(rap1 - rap2);
 372         ibin = Int_t(dy/fP2Step);
 373         if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
 374         if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
 375         if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
 376         if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
 377       }
 378     }
 379   }//case 0
 380   if(fAnalysisType==1) {
 381     for(i = 1; i < gNtrack; i++) {
 382       if(fAnalysisLevel == "ESD")
 383         track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
 384       else if(fAnalysisLevel == "AOD")
 385         track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
 386       else if(fAnalysisLevel == "MC")
 387         track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
 388       Short_t charge1 = track1->Charge();
 389       Double_t pZ1 = track1->Pz();
 390       Double_t p1 = track1->P();
 391       for(j = 0; j < i; j++) {
 392         if(fAnalysisLevel == "ESD")
 393           track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
 394         else if(fAnalysisLevel == "AOD")
 395           track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
 396         else if(fAnalysisLevel == "MC")
 397           track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
 398         Short_t charge2 = track2->Charge();
 399         Double_t pZ2 = track2->Pz();
 400         Double_t p2 = track2->P();
 401         Double_t eta1 = 0.5*log((p1 + pZ1)/(p1 - pZ1));
 402         Double_t eta2 = 0.5*log((p2 + pZ2)/(p2 - pZ2));
 403         Double_t deta = TMath::Abs(eta1 - eta2);
 404         ibin = Int_t(deta/fP2Step);
 405
 406         if((charge1 > 0.)&&(charge2 > 0.)) fNpp[ibin] += 1.;
 407         if((charge1 < 0.)&&(charge2 < 0.)) fNnn[ibin] += 1.;
 408         if((charge1 > 0.)&&(charge2 < 0.)) fNpn[ibin] += 1.;
 409         if((charge1 < 0.)&&(charge2 > 0.)) fNpn[ibin] += 1.;
 410         //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]);
 411       }
 412     }
 413   }//case 1
 414   if(fAnalysisType==2) {
 415     for(i = 1; i < gNtrack; i++) {
 416       if(fAnalysisLevel == "ESD")
 417         track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
 418       else if(fAnalysisLevel == "AOD")
 419         track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
 420       else if(fAnalysisLevel == "MC")
 421         track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
 422       Short_t charge1 = track1->Charge();
 423       Double_t pX1 = track1->Px();
 424       Double_t pY1 = track1->Py();
 425       Double_t pZ1 = track1->Pz();
 426       Double_t energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
 427                                      TMath::Power(track1->M(),2));
 428       for(j = 0; j < i; j++) {
 429         if(fAnalysisLevel == "ESD")
 430           track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
 431         else if(fAnalysisLevel == "AOD")
 432           track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
 433         else if(fAnalysisLevel == "MC")
 434           track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
 435         Short_t charge2 = track2->Charge();
 436         Double_t pX2 = track2->Px();
 437         Double_t pY2 = track2->Py();
 438         Double_t pZ2 = track2->Pz();
 439         Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
 440                                        TMath::Power(track2->M(),2));
 441         Double_t eTot = energy1 + energy2;
 442         Double_t pxTot = pX1 + pX2;
 443         Double_t pyTot = pY1 + pY2;
 444         Double_t pzTot = pZ1 + pZ2;
 445         Double_t q0Tot = energy1 - energy2;
 446         Double_t qzTot = pZ1 - pZ2;
 447         Double_t snn = TMath::Power(eTot,2) - TMath::Power(pxTot,2) - TMath::Power(pyTot,2) - TMath::Power(pzTot,2);
 448         Double_t ptTot = TMath::Sqrt( TMath::Power(pxTot,2) + TMath::Power(pyTot,2));
 449         Double_t qLong = TMath::Abs(eTot*qzTot - pzTot*q0Tot)/TMath::Sqrt(snn + TMath::Power(ptTot,2));
 450         ibin = Int_t(qLong/fP2Step);
 451         //cout<<ibin<<endl;
 452         if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
 453         if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
 454         if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
 455         if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
 456       }
 457     }
 458   }//case 2
 459   if(fAnalysisType==3) {
 460     for(i = 1; i < gNtrack; i++) {
 461       if(fAnalysisLevel == "ESD")
 462         track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
 463       else if(fAnalysisLevel == "AOD")
 464         track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
 465       else if(fAnalysisLevel == "MC")
 466         track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
 467       Short_t charge1 = track1->Charge();
 468       Double_t pX1 = track1->Px();
 469       Double_t pY1 = track1->Py();
 470       Double_t pZ1 = track1->Pz();
 471       Double_t energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
 472                                      TMath::Power(track1->M(),2));
 473       for(j = 0; j < i; j++) {
 474         if(fAnalysisLevel == "ESD")
 475           track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
 476         else if(fAnalysisLevel == "AOD")
 477           track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
 478         else if(fAnalysisLevel == "MC")
 479           track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
 480         Short_t charge2 = track2->Charge();
 481         Double_t pX2 = track2->Px();
 482         Double_t pY2 = track2->Py();
 483         Double_t pZ2 = track2->Pz();
 484         Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
 485                                        TMath::Power(track2->M(),2));
 486         Double_t eTot = energy1 + energy2;
 487         Double_t pxTot = pX1 + pX2;
 488         Double_t pyTot = pY1 + pY2;
 489         Double_t pzTot = pZ1 + pZ2;
 490         Double_t qxTot = pX1 - pX2;
 491         Double_t qyTot = pY1 - pY2;
 492         Double_t snn = TMath::Power(eTot,2) - TMath::Power(pxTot,2) - TMath::Power(pyTot,2) - TMath::Power(pzTot,2);
 493         Double_t ptTot = TMath::Sqrt( TMath::Power(pxTot,2) + TMath::Power(pyTot,2));
 494         Double_t qOut = TMath::Sqrt(snn/(snn + TMath::Power(ptTot,2))) * TMath::Abs(pxTot*qxTot + pyTot*qyTot)/ptTot;
 495         ibin = Int_t(qOut/fP2Step);
 496         //cout<<ibin<<endl;
 497         if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
 498         if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
 499         if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
 500         if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
 501       }
 502     }
 503   }//case 3
 504   if(fAnalysisType==4) {
 505     for(i = 1; i < gNtrack; i++) {
 506       if(fAnalysisLevel == "ESD")
 507         track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
 508       else if(fAnalysisLevel == "AOD")
 509         track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
 510       else if(fAnalysisLevel == "MC")
 511         track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
 512       Short_t charge1 = track1->Charge();
 513       Double_t pX1 = track1->Px();
 514       Double_t pY1 = track1->Py();
 515       //Double_t energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
 516       //TMath::Power(track1->M(),2));
 517       for(j = 0; j < i; j++) {
 518         if(fAnalysisLevel == "ESD")
 519           track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
 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));
 524         Short_t charge2 = track2->Charge();
 525         Double_t pX2 = track2->Px();
 526         Double_t pY2 = track2->Py();
 527         //Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
 528         //TMath::Power(track2->M(),2));
 529         //Double_t eTot = energy1 + energy2;
 530         Double_t pxTot = pX1 + pX2;
 531         Double_t pyTot = pY1 + pY2;
 532         Double_t qxTot = pX1 - pX2;
 533         Double_t qyTot = pY1 - pY2;
 534         Double_t ptTot = TMath::Sqrt( TMath::Power(pxTot,2) + TMath::Power(pyTot,2));
 535         Double_t qSide = TMath::Abs(pxTot*qyTot - pyTot*qxTot)/ptTot;
 536         ibin = Int_t(qSide/fP2Step);
 537         //cout<<ibin<<endl;
 538         if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
 539         if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
 540         if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
 541         if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
 542       }
 543     }
 544   }//case 4
 545   if(fAnalysisType==5) {
 546     for(i = 1; i < gNtrack; i++) {
 547       if(fAnalysisLevel == "ESD")
 548         track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
 549       else if(fAnalysisLevel == "AOD")
 550         track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
 551       else if(fAnalysisLevel == "MC")
 552         track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
 553       Short_t charge1 = track1->Charge();
 554       Double_t pX1 = track1->Px();
 555       Double_t pY1 = track1->Py();
 556       Double_t pZ1 = track1->Pz();
 557       Double_t energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
 558                                      TMath::Power(track1->M(),2));
 559       for(j = 0; j < i; j++) {
 560         if(fAnalysisLevel == "ESD")
 561           track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
 562         else if(fAnalysisLevel == "AOD")
 563           track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
 564         else if(fAnalysisLevel == "MC")
 565           track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
 566         Short_t charge2 = track2->Charge();
 567         Double_t pX2 = track2->Px();
 568         Double_t pY2 = track2->Py();
 569         Double_t pZ2 = track2->Pz();
 570         Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
 571                                        TMath::Power(track2->M(),2));
 572         Double_t q0Tot = energy1 - energy2;
 573         Double_t qxTot = pX1 - pX2;
 574         Double_t qyTot = pY1 - pY2;
 575         Double_t qzTot = pZ1 - pZ2;
 576         Double_t qInv = TMath::Sqrt(TMath::Abs(-TMath::Power(q0Tot,2) +TMath::Power(qxTot,2) +TMath::Power(qyTot,2) +TMath::Power(qzTot,2)));
 577         ibin = Int_t(qInv/fP2Step);
 578         //cout<<ibin<<endl;
 579         if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
 580         if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
 581         if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
 582         if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
 583       }
 584     }
 585   }//case 5
 586   if(fAnalysisType==6) {
 587     for(i = 1; i < gNtrack; i++) {
 588       if(fAnalysisLevel == "ESD")
 589         track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
 590       else if(fAnalysisLevel == "AOD")
 591         track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
 592       else if(fAnalysisLevel == "MC")
 593         track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
 594       Short_t charge1 = track1->Charge();
 595       Double_t phi1 = track1->Phi();
 596       for(j = 0; j < i; j++) {
 597         if(fAnalysisLevel == "ESD")
 598           track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
 599         else if(fAnalysisLevel == "AOD")
 600           track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
 601         else if(fAnalysisLevel == "MC")
 602           track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
 603         Short_t charge2 = track2->Charge();
 604         Double_t phi2 = track2->Phi();
 605         Double_t dphi = TMath::Abs(phi1 - phi2);
 606         ibin = Int_t(dphi/fP2Step);
 607         if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
 608         if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
 609         if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
 610         if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
 611       }
 612     }
 613   }//case 6
 614
 615   /*for(Int_t i = 0; i < fNumberOfBins; i++)
 616     Printf("bin: %d - Npp: %lf - Nnn: %lf - Nnp: %lf - Npn: %lf",i,fNpp[i],fNnn[i],fNpn[i],fNpn[i]);*/
 617 }
 618
 619 //____________________________________________________________________//
 620 TH1F *AliBalance::GetHistNnn() {
 621   //Return the control histogram of the -- component
 622   for(Int_t iBin = 1; iBin <= fNumberOfBins; iBin++)
 623     fHistfNnn->SetBinContent(iBin,GetNnn(iBin-1));
 624
 625   return fHistfNnn;
 626 }
 627
 628 //____________________________________________________________________//
 629 TH1F *AliBalance::GetHistNpp() {
 630   //Return the control histogram of the ++ component
 631   for(Int_t iBin = 1; iBin <= fNumberOfBins; iBin++)
 632     fHistfNpp->SetBinContent(iBin,GetNpp(iBin-1));
 633
 634   return fHistfNpp;
 635 }
 636
 637 //____________________________________________________________________//
 638 TH1F *AliBalance::GetHistNpn() {
 639   //Return the control histogram of the +- component
 640   for(Int_t iBin = 1; iBin <= fNumberOfBins; iBin++)
 641     fHistfNpn->SetBinContent(iBin,GetNpn(iBin-1));
 642
 643   return fHistfNpn;
 644 }
 645
 646 //____________________________________________________________________//
 647 Double_t AliBalance::GetBalance(Int_t p2) {
 648   // Returns the value of the balance function in bin p2
 649   fB[p2] = 0.5*(((fNpn[p2] - 2.0*fNnn[p2])/fNn) + ((fNpn[p2] - 2.0*fNpp[p2])/fNp))/fP2Step;
 650
 651   return fB[p2];
 652 }
 653
 654 //____________________________________________________________________//
 655 Double_t AliBalance::GetError(Int_t p2) {
 656   // Returns the error on the BF value for bin p2
 657   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;
 658
 659   return ferror[p2];
 660 }
 661
 662 //____________________________________________________________________//
 663 void AliBalance::PrintResults() {
 664   // Prints the results
 665   Double_t x[MAXIMUM_NUMBER_OF_STEPS];
 666   Double_t fSumXi = 0.0, fSumBi = 0.0, fSumBiXi = 0.0;
 667   Double_t fSumBiXi2 = 0.0, fSumBi2Xi2 = 0.0;
 668   Double_t fSumDeltaBi2 = 0.0, fSumXi2DeltaBi2 = 0.0;
 669   Double_t deltaBalP2 = 0.0, integral = 0.0;
 670   Double_t deltaErrorNew = 0.0;
 671
 672   cout<<"=================================================="<<endl;
 673   for(Int_t i = 0; i < fNumberOfBins; i++) {
 674     //x[i] = fP2Start + fP2Step*i + fP2Step/2;
 675     x[i] = fP2Step*i + fP2Step/2;
 676     cout<<"B: "<<fB[i]<<"\t Error: "<<ferror[i]<<"\t bin: "<<x[i]<<endl;
 677   }
 678   cout<<"=================================================="<<endl;
 679   for(Int_t i = 1; i < fNumberOfBins; i++) {
 680     fSumXi += x[i];
 681     fSumBi += fB[i];
 682     fSumBiXi += fB[i]*x[i];
 683     fSumBiXi2 += fB[i]*TMath::Power(x[i],2);
 684     fSumBi2Xi2 += TMath::Power(fB[i],2)*TMath::Power(x[i],2);
 685     fSumDeltaBi2 +=  TMath::Power(ferror[i],2);
 686     fSumXi2DeltaBi2 += TMath::Power(x[i],2) * TMath::Power(ferror[i],2);
 687
 688     deltaBalP2 += fP2Step*TMath::Power(ferror[i],2);
 689     integral += fP2Step*fB[i];
 690   }
 691   for(Int_t i = 1; i < fNumberOfBins; i++) deltaErrorNew += ferror[i]*(x[i]*fSumBi - fSumBiXi)/TMath::Power(fSumBi,2);
 692
 693   Double_t integralError = TMath::Sqrt(deltaBalP2);
 694
 695   Double_t delta = fSumBiXi / fSumBi;
 696   Double_t deltaError = (fSumBiXi / fSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(fSumXi2DeltaBi2)/fSumBiXi),2) + TMath::Power((fSumDeltaBi2/fSumBi),2) );
 697
 698   cout<<"Analyzed events: "<<fAnalyzedEvents<<endl;
 699   cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
 700   cout<<"New error: "<<deltaErrorNew<<endl;
 701   cout<<"Interval: "<<integral<<"\t Error: "<<integralError<<endl;
 702   cout<<"=================================================="<<endl;
 703 }
 704
 705 //____________________________________________________________________//
 706 TGraphErrors *AliBalance::DrawBalance() {
 707   // Draws the BF
 708   Double_t x[MAXIMUM_NUMBER_OF_STEPS];
 709   Double_t xer[MAXIMUM_NUMBER_OF_STEPS];
 710   Double_t b[MAXIMUM_NUMBER_OF_STEPS];
 711   Double_t ber[MAXIMUM_NUMBER_OF_STEPS];
 712
 713   if((fNp == 0)||(fNn == 0)) {
 714     cout<<"Couldn't find any particles in the analyzed interval!!!"<<endl;
 715     cout<<"Aborting....."<<endl;
 716     abort();
 717   }
 718
 719   for(Int_t i = 0; i < fNumberOfBins; i++) {
 720     b[i] = GetBalance(i);
 721     ber[i] = GetError(i);
 722     //x[i] = fP2Start + fP2Step*i + fP2Step/2;
 723     x[i] = fP2Step*i + fP2Step/2;
 724     xer[i] = 0.0;
 725   }
 726
 727   TGraphErrors *gr = new TGraphErrors(fNumberOfBins,x,b,xer,ber);
 728   gr->SetMarkerStyle(25);
 729   gr->GetXaxis()->SetTitleColor(1);
 730   if(fAnalysisType==0) {
 731     gr->GetXaxis()->SetTitle("#Delta y");
 732     gr->GetYaxis()->SetTitle("B(#Delta y)");
 733   }
 734   if(fAnalysisType==1) {
 735     gr->GetXaxis()->SetTitle("#Delta #eta");
 736     gr->GetYaxis()->SetTitle("B(#Delta #eta)");
 737   }
 738   if(fAnalysisType==2) {
 739     gr->GetXaxis()->SetTitle("q_{long} (GeV/c)");
 740     gr->GetYaxis()->SetTitle("B(q_{long}) [(GeV/c)^{-1}]");
 741   }
 742   if(fAnalysisType==3) {
 743     gr->GetXaxis()->SetTitle("q_{out} (GeV/c)");
 744     gr->GetYaxis()->SetTitle("B(q_{out}) [(GeV/c)^{-1}]");
 745   }
 746   if(fAnalysisType==4) {
 747     gr->GetXaxis()->SetTitle("q_{side} (GeV/c)");
 748     gr->GetYaxis()->SetTitle("B(q_{side}) [(GeV/c)^{-1}]");
 749   }
 750   if(fAnalysisType==5) {
 751     gr->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
 752     gr->GetYaxis()->SetTitle("B(q_{inv}) [(GeV/c)^{-1}]");
 753   }
 754   if(fAnalysisType==6) {
 755     gr->GetXaxis()->SetTitle("#Delta #phi");
 756     gr->GetYaxis()->SetTitle("B(#Delta #phi)");
 757   }
 758
 759   return gr;
 760 }
 761
 762 //____________________________________________________________________//
 763 void AliBalance::Merge(AliBalance *b) {
 764   //Merging function to be used for proof and grid
 765   fNp += b->GetNp();
 766   fNn += b->GetNn();
 767   for(Int_t i = 0; i < fNumberOfBins; i++) {
 768     fNnn[i] += b->GetNnn(i);
 769     fNpp[i] += b->GetNpp(i);
 770     fNpn[i] += b->GetNpn(i);
 771   }
 772 }