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