ANALYSIS/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 <TLorentzVector.h>
  27 #include <TGraphErrors.h>
  28
  29 #include "AliBalance.h"
  30
  31 ClassImp(AliBalance)
  32
  33 //----------------------------------------//
  34 AliBalance::AliBalance()
  35 {
  36   // Default constructor
  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     }
  45   fNp = 0.0;
  46   fNn = 0.0;
  47   fP2Start = 0.0;
  48   fP2Stop = 0.0;
  49   fP2Step = 0.0;
  50   fAnalysisType = 0;
  51   fNumberOfBins = 0;
  52   fNtrack = 0;
  53   fAnalyzedEvents = 0;
  54 }
  55
  56 //----------------------------------------//
  57 AliBalance::AliBalance(Double_t p2Start, Double_t p2Stop, Int_t p2Bins)
  58 {
  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;
  78 }
  79
  80 //----------------------------------------//
  81 AliBalance::~AliBalance()
  82 {
  83   // Destructor
  84 }
  85
  86 //----------------------------------------//
  87 void AliBalance::SetNumberOfBins(Int_t ibins)
  88 {
  89   // Sets the number of bins for the analyzed interval
  90   fNumberOfBins = ibins ;
  91 }
  92
  93 //----------------------------------------//
  94 void AliBalance::SetInterval(Double_t p2Start, Double_t p2Stop)
  95 {
  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;
 101 }
 102
 103 //----------------------------------------//
 104 void 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 //----------------------------------------//
 159 void AliBalance::SetParticles(TLorentzVector *P, Double_t *charge, Int_t dim)
 160 {
 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.
 164   this->fV = P;
 165   this->fCharge = charge;
 166   fNtrack = dim;
 167 }
 168
 169
 170 //----------------------------------------//
 171 void AliBalance::CalculateBalance()
 172 {
 173   // Calculates the balance function
 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             {
 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()));
 195               Double_t dy = TMath::Abs(rap1 - rap2);
 196               ibin = Int_t(dy/fP2Step);
 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             {
 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));
 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()));
 218               Double_t deta = TMath::Abs(eta1 - eta2);
 219               ibin = Int_t(deta/fP2Step);
 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             {
 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);
 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             {
 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);
 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             {
 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);
 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             {
 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);
 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);
 345               ibin = Int_t(dphi/fP2Step);
 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 //----------------------------------------//
 360 Double_t AliBalance::GetBalance(Int_t p2)
 361 {
 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;
 364
 365   return fB[p2];
 366 }
 367
 368 //----------------------------------------//
 369 Double_t AliBalance::GetError(Int_t p2)
 370 {
 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;
 373
 374   return ferror[p2];
 375 }
 376
 377 //----------------------------------------//
 378 void AliBalance::PrintResults()
 379 {
 380   // Prints the results
 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;
 385   Double_t deltaBalP2 = 0.0, integral = 0.0;
 386   Double_t deltaErrorNew = 0.0;
 387
 388   cout<<"=================================================="<<endl;
 389   for(Int_t i = 0; i < fNumberOfBins; i++)
 390     {
 391       x[i] = fP2Start + fP2Step*i + fP2Step/2 ;
 392       cout<<"B: "<<fB[i]<<"\t Error: "<<ferror[i]<<"\t bin: "<<x[i]<<endl;
 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
 405       deltaBalP2 += fP2Step*pow(ferror[i],2) ;
 406       integral += fP2Step*fB[i] ;
 407     }
 408   for(Int_t i = 1; i < fNumberOfBins; i++)
 409     {
 410       deltaErrorNew += ferror[i]*(x[i]*fSumBi - fSumBiXi)/pow(fSumBi,2);
 411     }
 412   Double_t integralError = sqrt(deltaBalP2) ;
 413
 414   Double_t delta = fSumBiXi / fSumBi ;
 415   Double_t deltaError = (fSumBiXi / fSumBi) * sqrt(pow((sqrt(fSumXi2DeltaBi2)/fSumBiXi),2) + pow((fSumDeltaBi2/fSumBi),2) ) ;
 416
 417   cout<<"Analyzed events: "<<fAnalyzedEvents<<endl;
 418   cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
 419   cout<<"New error: "<<deltaErrorNew<<endl;
 420   cout<<"Interval: "<<integral<<"\t Error: "<<integralError<<endl;
 421   cout<<"=================================================="<<endl;
 422 }
 423
 424 //----------------------------------------//
 425 TGraphErrors *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 }