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   fCharge(0),
  72   fNtrack(balance.fNtrack),
  73   fV(0),
  74   fNumberOfBins(balance.fNumberOfBins),
  75   fAnalysisType(balance.fAnalysisType),
  76   fAnalyzedEvents(balance.fAnalyzedEvents),
  77   fP2Start(balance.fP2Start),
  78   fP2Stop(balance.fP2Stop),
  79   fP2Step(balance.fP2Step),
  80   fNn(balance.fNn),
  81   fNp(balance.fNp)
  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 }