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