Balance function analysis (P.Christakoglou)
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 17 Aug 2005 13:31:53 +0000 (13:31 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 17 Aug 2005 13:31:53 +0000 (13:31 +0000)
ANALYSIS/ANALYSISLinkDef.h
ANALYSIS/AliAnalysisEventCuts.cxx [new file with mode: 0644]
ANALYSIS/AliAnalysisEventCuts.h [new file with mode: 0644]
ANALYSIS/AliAnalysisTrackCuts.cxx [new file with mode: 0644]
ANALYSIS/AliAnalysisTrackCuts.h [new file with mode: 0644]
ANALYSIS/AliBalance.cxx [new file with mode: 0644]
ANALYSIS/AliBalance.h [new file with mode: 0644]
ANALYSIS/libANALYSIS.pkg

index f632c2a..2359005 100644 (file)
@@ -93,4 +93,8 @@
     
 #pragma link C++ class AliAODv0+;
 
+#pragma link C++ class AliAnalysisEventCuts+;
+#pragma link C++ class AliAnalysisTrackCuts+;
+#pragma link C++ class AliBalance+;
+
 #endif
diff --git a/ANALYSIS/AliAnalysisEventCuts.cxx b/ANALYSIS/AliAnalysisEventCuts.cxx
new file mode 100644 (file)
index 0000000..32ea4e5
--- /dev/null
@@ -0,0 +1,295 @@
+/**************************************************************************
+ * Author: Panos Christakoglou.                                           *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+//-----------------------------------------------------------------
+//           AliAnalysisEventCuts class
+//   This is the class to deal with the event and track level cuts
+//   Origin: Panos Christakoglou, UOA-CERN, Panos.Christakoglou@cern.ch
+//-----------------------------------------------------------------
+
+
+
+//ROOT
+#include <TROOT.h>
+#include <TObject.h>
+#include <TSystem.h>
+#include <TObject.h>
+
+#include <TPaveText.h>
+#include <TText.h>
+#include <TLine.h>
+#include <TCanvas.h>
+
+#include <TControlBar.h>
+#include <TRootControlBar.h>
+
+
+
+#include "AliLog.h"
+#include "AliESD.h"
+
+#include "AliAnalysisEventCuts.h"
+
+ClassImp(AliAnalysisEventCuts)
+
+//----------------------------------------//
+AliAnalysisEventCuts::AliAnalysisEventCuts()
+{
+  Reset();
+}
+
+//----------------------------------------//
+AliAnalysisEventCuts::~AliAnalysisEventCuts()
+{
+}
+
+//----------------------------------------//
+void AliAnalysisEventCuts::Reset()
+{
+  fVxMin = -1000.0;
+  fVxMax = 1000.0; 
+  fVyMin = -1000.0;
+  fVyMax = 1000.0;  
+  fVzMin = -1000.0;
+  fVzMax = 1000.0;
+  fMultMin = 0;
+  fMultMax = 100000;
+  
+  fMult = 0;  
+  fVx = 0;  
+  fVy = 0; 
+  fVz = 0; 
+  fTotalEvents = 0; 
+  fAcceptedEvents = 0; 
+
+  fFlagMult = 0;  
+  fFlagVx = 0;  
+  fFlagVy = 0; 
+  fFlagVz = 0; 
+}
+
+//----------------------------------------//
+void AliAnalysisEventCuts::SetPrimaryVertexXRange(Float_t r1, Float_t r2)
+{
+  fVxMin = r1;
+  fVxMax = r2; 
+  fFlagVx = 1;  
+}
+
+//----------------------------------------//
+void AliAnalysisEventCuts::SetPrimaryVertexYRange(Float_t r1, Float_t r2)
+{
+  fVyMin = r1;
+  fVyMax = r2; 
+  fFlagVy = 1;
+}
+
+//----------------------------------------//
+void AliAnalysisEventCuts::SetPrimaryVertexZRange(Float_t r1, Float_t r2)
+{
+  fVzMin = r1;
+  fVzMax = r2; 
+  fFlagVz = 1;
+}
+
+//----------------------------------------//
+void AliAnalysisEventCuts::SetMultiplicityRange(Int_t n1, Int_t n2)
+{
+  fMultMin = n1;
+  fMultMax = n2; 
+  fFlagMult = 1;
+}
+
+
+//----------------------------------------//
+Bool_t AliAnalysisEventCuts::IsAccepted(AliESD *esd)
+{
+  fTotalEvents++;
+  if((esd->GetNumberOfTracks() < fMultMin) || (esd->GetNumberOfTracks() > fMultMax))
+    {
+      fMult++;
+      AliInfo(Form("Event rejected due to multiplicity cut"));
+      return kFALSE;
+    }
+  const AliESDVertex *esdvertex = esd->GetVertex();
+  if((esdvertex->GetXv() < fVxMin) || (esdvertex->GetXv() > fVxMax))
+    {
+      fVx++;
+      AliInfo(Form("Event rejected due to Vx cut"));
+      return kFALSE;
+    }
+  if((esdvertex->GetYv() < fVyMin) || (esdvertex->GetYv() > fVyMax))
+    {
+      fVy++;
+      AliInfo(Form("Event rejected due to Vy cut"));
+      return kFALSE;
+    }
+ if((esdvertex->GetZv() < fVzMin) || (esdvertex->GetZv() > fVzMax))
+    {
+      fVz++;
+      AliInfo(Form("Event rejected due to Vz cut"));
+      return kFALSE;
+    }
+ fAcceptedEvents++;
+
+ return kTRUE;
+}
+
+
+//----------------------------------------//
+TPaveText *AliAnalysisEventCuts::GetEventCuts()
+{
+  TCanvas *ccuts = new TCanvas("ccuts","Event cuts",10,10,400,400);
+  ccuts->SetFillColor(10);
+  ccuts->SetHighLightColor(10);
+
+  TPaveText *pave = new TPaveText(0.01,0.01,0.98,0.98);
+  pave->SetFillColor(5);
+  Char_t CutName[256];
+  TLine *l1 = pave->AddLine(0,0.78,1,0.78);
+  l1->SetLineWidth(2);
+  TLine *l2 = pave->AddLine(0,0.58,1,0.58);
+  l2->SetLineWidth(2);
+  TLine *l3 = pave->AddLine(0,0.38,1,0.38);
+  l3->SetLineWidth(2);
+  TLine *l4 = pave->AddLine(0,0.18,1,0.18);
+  l4->SetLineWidth(2);
+  sprintf(CutName,"Total number of events: %d",fTotalEvents);
+  TText *t1 = pave->AddText(CutName);
+  t1->SetTextColor(4);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+  sprintf(CutName,"Total number of accepted events: %d",fAcceptedEvents);
+  t1 = pave->AddText(CutName);
+  t1->SetTextColor(4);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+  sprintf(CutName,"Multiplicity range: [%d,%d]",fMultMin,fMultMax);
+  t1 = pave->AddText(CutName);
+  t1->SetTextColor(4);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+  sprintf(CutName,"Events rejected: %d",fMult);
+  t1 = pave->AddText(CutName);
+  t1->SetTextColor(4);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+  sprintf(CutName,"Vx range: [%f,%f]",fVxMin,fVxMax);
+  t1 = pave->AddText(CutName);
+  t1->SetTextColor(4);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+  sprintf(CutName,"Events rejected: %d",fVx);
+  t1 = pave->AddText(CutName);
+  t1->SetTextColor(4);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+  sprintf(CutName,"Vy range: [%f,%f]",fVyMin,fVyMax);
+  t1 = pave->AddText(CutName);
+  t1->SetTextColor(4);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+  sprintf(CutName,"Events rejected: %d",fVy);
+  t1 = pave->AddText(CutName);
+  t1->SetTextColor(4);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+  sprintf(CutName,"Vz range: [%f,%f]",fVzMin,fVzMax);
+  t1 = pave->AddText(CutName);
+  t1->SetTextColor(4);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+  sprintf(CutName,"Events rejected: %d",fVz);
+  t1 = pave->AddText(CutName);
+  t1->SetTextColor(4);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+  return pave;
+}
+
+//----------------------------------------//
+void AliAnalysisEventCuts::GetEventStats()
+{
+  AliInfo(Form("Total number of events: %d",fTotalEvents));
+  AliInfo(Form("Total number of accepted events: %d",fAcceptedEvents)); 
+}
+
+//----------------------------------------//
+void AliAnalysisEventCuts::GetMultStats()
+{
+  AliInfo(Form("Multiplicity range: [%d,%d]",fMultMin,fMultMax));
+  AliInfo(Form("Events rejected: %d",fMult));
+}
+
+//----------------------------------------//
+void AliAnalysisEventCuts::GetVxStats()
+{
+  AliInfo(Form("Vx range: [%f,%f]",fVxMin,fVxMax));
+  AliInfo(Form("Events rejected: %d",fVx));
+}
+
+//----------------------------------------//
+void AliAnalysisEventCuts::GetVyStats()
+{
+  AliInfo(Form("Vy range: [%f,%f]",fVyMin,fVyMax));
+  AliInfo(Form("Events rejected: %d",fVy));
+}
+
+//----------------------------------------//
+void AliAnalysisEventCuts::GetVzStats()
+{
+  AliInfo(Form("Vz range: [%f,%f]",fVzMin,fVzMax));
+  AliInfo(Form("Events rejected: %d",fVz));
+}
+
+//----------------------------------------//
+void AliAnalysisEventCuts::PrintEventCuts()
+{
+  /*gROOT->Reset();
+  TControlBar *menu1 = new TControlBar("vertical","Event Cuts",10,10);
+  menu1->AddButton("Event statistics","GetEventStats()","Displays the event statistics");
+  menu1->AddButton("Multipicity cut","gAlice->Run()","Events rejected due to multiplicity cut");
+  menu1->AddButton("Vx cut","gAlice->RunLego()","Events rejected due to Vx cut");
+  menu1->AddButton("Vy cut","DrawTopView()","Events rejected due to Vy cut");
+  menu1->AddButton("Vz cut","DrawSideView()","Events rejected due to Vz cut");
+  menu1->Show();
+  gROOT->SaveContext(); */
+
+  //GetEventCuts()->Draw();  
+
+  AliInfo(Form("**************EVENT CUTS**************"));
+  GetEventStats();
+  if(fFlagMult)
+    GetMultStats();
+  if(fFlagVx)
+    GetVxStats();
+  if(fFlagVy)
+    GetVyStats();
+  if(fFlagVz)
+    GetVzStats();
+  AliInfo(Form("**************************************"));
+}
+
+
+     
diff --git a/ANALYSIS/AliAnalysisEventCuts.h b/ANALYSIS/AliAnalysisEventCuts.h
new file mode 100644 (file)
index 0000000..b9ab4c2
--- /dev/null
@@ -0,0 +1,68 @@
+#ifndef ALIANALYSISEVENTCUTS_H
+#define ALIANALYSISEVENTCUTS_H
+/*  See cxx source for full Copyright notice */
+
+
+/* $Id$ */
+
+//-------------------------------------------------------------------------
+//                       Class AliAnalysisEventCuts
+//   This is the class for the cuts in event & track level
+//
+//    Origin: Panos Christakoglou, UOA-CERN, Panos.Christakoglou@cern.ch
+//-------------------------------------------------------------------------
+
+#include <TPaveText.h>
+#include <TROOT.h>
+#include <TObject.h>
+
+#include "AliESD.h"
+
+class AliAnalysisEventCuts : public TObject
+{
+ private:
+  Float_t fVxMin, fVxMax;  //Definition of the range of the Vx
+  Float_t fVyMin, fVyMax;  //Definition of the range of the Vy
+  Float_t fVzMin, fVzMax;  //Definition of the range of the Vz
+  Int_t fMultMin, fMultMax;  //Definition of the range of the multiplicity
+
+  Int_t fMult;  //Number of events rejected due to multiplicity cut
+  Int_t fVx;  //Number of events rejected due to Vx cut
+  Int_t fVy;  //Number of events rejected due to Vy cut
+  Int_t fVz;  //Number of events rejected due to Vz cut
+  Int_t fTotalEvents;  //Total number of events
+  Int_t fAcceptedEvents;  //Total number of events accepted
+
+  Int_t fFlagMult; //Flag that shows if the multiplicity cut was imposed
+  Int_t fFlagVx; //Flag that shows if the Vx cut was imposed
+  Int_t fFlagVy; //Flag that shows if the Vy cut was imposed
+  Int_t fFlagVz; //Flag that shows ifthe Vz cut was imposed
+ public:
+  AliAnalysisEventCuts();
+  
+  ~AliAnalysisEventCuts();
+
+  void Reset();
+  
+  void SetPrimaryVertexXRange(Float_t r1, Float_t r2);
+  void SetPrimaryVertexYRange(Float_t r1, Float_t r2);
+  void SetPrimaryVertexZRange(Float_t r1, Float_t r2);
+  void SetMultiplicityRange(Int_t n1, Int_t n2);
+  Bool_t IsAccepted(AliESD *esd);
+
+  TPaveText *GetEventCuts();
+  void PrintEventCuts(); 
+  void GetEventStats();
+  void GetMultStats();
+  void GetVxStats();
+  void GetVyStats();
+  void GetVzStats();
+  ClassDef(AliAnalysisEventCuts, 1)
+} ;
+
+#endif
diff --git a/ANALYSIS/AliAnalysisTrackCuts.cxx b/ANALYSIS/AliAnalysisTrackCuts.cxx
new file mode 100644 (file)
index 0000000..10edd69
--- /dev/null
@@ -0,0 +1,499 @@
+/**************************************************************************
+ * Author: Panos Christakoglou.                                           *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+//-----------------------------------------------------------------
+//           AliAnalysisTrackCuts class
+//   This is the class to deal with the event and track level cuts
+//   Origin: Panos Christakoglou, UOA-CERN, Panos.Christakoglou@cern.ch
+//-----------------------------------------------------------------
+
+
+
+//ROOT
+#include <TROOT.h>
+#include <TObject.h>
+#include <TSystem.h>
+#include <TObject.h>
+#include <TPaveText.h>
+#include <TText.h>
+#include <TLine.h>
+#include <TCanvas.h>
+
+#include "AliLog.h"
+#include "AliESDtrack.h"
+#include "AliESD.h"
+
+#include "AliAnalysisTrackCuts.h"
+
+ClassImp(AliAnalysisTrackCuts)
+
+//----------------------------------------//
+AliAnalysisTrackCuts::AliAnalysisTrackCuts()
+{
+  Reset();
+}
+
+//----------------------------------------//
+AliAnalysisTrackCuts::~AliAnalysisTrackCuts()
+{
+}
+
+//----------------------------------------//
+void AliAnalysisTrackCuts::Reset()
+{
+  fPxMin = -1000.0;
+  fPxMax = 1000.0; 
+  fPyMin = -1000.0;
+  fPyMax = 1000.0;  
+  fPzMin = -1000.0;
+  fPzMax = 1000.0;
+  fPtMin = 0.0;
+  fPtMax = 1000.0;
+  fPMin = 0.0;
+  fPMax = 1000.0;
+  fBrMin = 0.0;
+  fBrMax = 1000.0;
+  fBzMin = 0.0;
+  fBzMax = 1000.0;
+  fEtaMin = -100.0;
+  fEtaMax = 100.0;
+  fRapMin = -100.0;
+  fRapMax = 100.0;
+  
+  fP = 0;  
+  fPt = 0; 
+  fPx = 0;  
+  fPy = 0; 
+  fPz = 0; 
+  fbr = 0; 
+  fbz = 0; 
+  fEta = 0;  
+  fRap = 0; 
+  fTotalTracks = 0; 
+  fAcceptedTracks = 0; 
+
+  fFlagP = 0;  
+  fFlagPt = 0;  
+  fFlagPx = 0;  
+  fFlagPy = 0;  
+  fFlagPz = 0;  
+  fFlagEta = 0;  
+  fFlagRap = 0;  
+  fFlagbr = 0;  
+  fFlagbz = 0;  
+}
+
+//----------------------------------------//
+void AliAnalysisTrackCuts::SetPxRange(Float_t r1, Float_t r2)
+{
+  fPxMin = r1;
+  fPxMax = r2;
+  fFlagPx = 1;
+}
+
+//----------------------------------------//
+void AliAnalysisTrackCuts::SetPyRange(Float_t r1, Float_t r2)
+{
+  fPyMin = r1;
+  fPyMax = r2; 
+  fFlagPy = 1;
+}
+
+//----------------------------------------//
+void AliAnalysisTrackCuts::SetPzRange(Float_t r1, Float_t r2)
+{
+  fPzMin = r1;
+  fPzMax = r2; 
+  fFlagPy = 1;
+}
+
+//----------------------------------------//
+void AliAnalysisTrackCuts::SetPRange(Float_t r1, Float_t r2)
+{
+  fPMin = r1;
+  fPMax = r2; 
+  fFlagPz = 1;
+}
+
+//----------------------------------------//
+void AliAnalysisTrackCuts::SetPtRange(Float_t r1, Float_t r2)
+{
+  fPtMin = r1;
+  fPtMax = r2; 
+  fFlagPt = 1;
+}
+
+//----------------------------------------//
+void AliAnalysisTrackCuts::SetBrRange(Float_t r1, Float_t r2)
+{
+  fBrMin = r1;
+  fBrMax = r2; 
+  fFlagbr = 1;
+}
+
+//----------------------------------------//
+void AliAnalysisTrackCuts::SetBzRange(Float_t r1, Float_t r2)
+{
+  fBzMin = r1;
+  fBzMax = r2; 
+  fFlagbz = 1;
+}
+
+//----------------------------------------//
+void AliAnalysisTrackCuts::SetEtaRange(Float_t r1, Float_t r2)
+{
+  fEtaMin = r1;
+  fEtaMax = r2; 
+  fFlagEta = 1;
+}
+
+//----------------------------------------//
+void AliAnalysisTrackCuts::SetRapRange(Float_t r1, Float_t r2)
+{
+  fRapMin = r1;
+  fRapMax = r2; 
+  fFlagRap = 1;
+}
+
+//----------------------------------------//
+void AliAnalysisTrackCuts::GetTrackStats()
+{
+  AliInfo(Form("Total number of tracks: %d",fTotalTracks));
+  AliInfo(Form("Total number of accepted tracks: %d",fAcceptedTracks)); 
+}
+
+//----------------------------------------//
+void AliAnalysisTrackCuts::GetPStats()
+{
+  AliInfo(Form("P range: [%f,%f]",fPMin,fPMax));
+  if(fTotalTracks != 0)
+    AliInfo(Form("Tracks rejected: %f",100.0*fP/fTotalTracks)); 
+}
+
+//----------------------------------------//
+void AliAnalysisTrackCuts::GetPtStats()
+{
+  AliInfo(Form("Pt range: [%f,%f]",fPtMin,fPtMax));
+  if(fTotalTracks != 0) 
+    AliInfo(Form("Tracks rejected: %f",100.0*fPt/fTotalTracks)); 
+}
+
+//----------------------------------------//
+void AliAnalysisTrackCuts::GetPxStats()
+{
+  AliInfo(Form("Px range: [%f,%f]",fPxMin,fPxMax));
+  if(fTotalTracks != 0) 
+    AliInfo(Form("Tracks rejected: %f",100.0*fPx/fTotalTracks)); 
+}
+
+//----------------------------------------//
+void AliAnalysisTrackCuts::GetPyStats()
+{
+  AliInfo(Form("Py range: [%f,%f]",fPyMin,fPyMax));
+  if(fTotalTracks != 0) 
+    AliInfo(Form("Tracks rejected: %f",100.0*fPy/fTotalTracks)); 
+}
+
+//----------------------------------------//
+void AliAnalysisTrackCuts::GetPzStats()
+{
+  AliInfo(Form("Pz range: [%f,%f]",fPzMin,fPzMax));
+  if(fTotalTracks != 0) 
+    AliInfo(Form("Tracks rejected: %f",100.0*fPz/fTotalTracks)); 
+}
+
+//----------------------------------------//
+void AliAnalysisTrackCuts::GetEtaStats()
+{
+  AliInfo(Form("eta range: [%f,%f]",fEtaMin,fEtaMax));
+  if(fTotalTracks != 0)
+    AliInfo(Form("Tracks rejected: %f",100.0*fEta/fTotalTracks)); 
+}
+
+//----------------------------------------//
+void AliAnalysisTrackCuts::GetRapStats()
+{
+  AliInfo(Form("y range: [%f,%f]",fRapMin,fRapMax));
+  if(fTotalTracks != 0)
+    AliInfo(Form("Tracks rejected: %f",100.0*fRap/fTotalTracks)); 
+}
+
+//----------------------------------------//
+void AliAnalysisTrackCuts::GetBrStats()
+{
+  AliInfo(Form("br range: [%f,%f]",fBrMin,fBrMax));
+  if(fTotalTracks != 0) 
+    AliInfo(Form("Tracks rejected: %f",100.0*fbr/fTotalTracks)); 
+}
+
+//----------------------------------------//
+void AliAnalysisTrackCuts::GetBzStats()
+{
+  AliInfo(Form("bz range: [%f,%f]",fBzMin,fBzMax));
+  if(fTotalTracks != 0)
+    AliInfo(Form("Tracks rejected: %f",100.0*fbz/fTotalTracks)); 
+}
+
+
+//----------------------------------------//
+Bool_t AliAnalysisTrackCuts::IsAccepted(AliESD *esd ,AliESDtrack *esdtrack)
+{
+  fTotalTracks++;
+  
+  //momentum related calculations
+  Double_t p[3];
+  esdtrack->GetPxPyPz(p);
+  Float_t P = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
+  Float_t Pt = sqrt(pow(p[0],2) + pow(p[1],2));
+  Float_t E = sqrt(pow(esdtrack->GetMass(),2) + pow(P,2));
+
+  //y-eta related calculations
+  Float_t eta = 0.5*log((P - p[2])/(P + p[2]));
+  Float_t y = 0.5*log((E - p[2])/(E + p[2]));
+  //impact parameter related calculations
+  Double_t TrackPosition[3];
+  esdtrack->GetXYZ(TrackPosition);
+  const AliESDVertex * VertexIn = esd->GetVertex();
+  Double_t VertexPosition[3];
+  VertexIn->GetXYZ(VertexPosition);                
+  for (Int_t ii=0; ii<3; ii++) TrackPosition[ii] -= VertexPosition[ii];
+                   
+  Float_t br = Float_t(TMath::Sqrt(pow(TrackPosition[0],2) + pow(TrackPosition[1],2)));
+  Float_t bz = Float_t(TMath::Abs(TrackPosition[2]));
+  if((P < fPMin) || (P > fPMax))
+    {
+      fP++;
+      return kFALSE;
+    }
+  if((Pt < fPtMin) || (Pt > fPtMax))
+    {
+      fPt++;
+      return kFALSE;
+    }
+  if((p[0] < fPxMin) || (p[0] > fPxMax))
+    {
+      fPx++;
+      return kFALSE;
+    }
+  if((p[1] < fPyMin) || (p[1] > fPyMax))
+    {
+      fPy++;
+      return kFALSE;
+    }
+  if((p[2] < fPzMin) || (p[2] > fPzMax))
+    {
+      fPz++;
+      return kFALSE;
+    } 
+  if((br < fBrMin) || (br > fBrMax))
+    {
+      fbr++;
+      return kFALSE;
+    }
+  if((bz < fBzMin) || (bz > fBzMax))
+    {
+      fbz++;
+      return kFALSE;
+    }
+  if((eta < fEtaMin) || (eta > fEtaMax))
+    {
+      fEta++;
+      return kFALSE;
+    }
+  if((y < fRapMin) || (y > fRapMax))
+    {
+      fRap++;
+      return kFALSE;
+    }
+  
+  fAcceptedTracks++;
+  
+  return kTRUE;
+}
+
+
+//----------------------------------------//
+TPaveText *AliAnalysisTrackCuts::GetTrackCuts()
+{
+  TCanvas *ccuts2 = new TCanvas("ccuts2","Track cuts",410,10,400,400);
+  ccuts2->SetFillColor(10);
+  ccuts2->SetHighLightColor(10);
+
+  TPaveText *pave = new TPaveText(0.01,0.01,0.98,0.98);
+  pave->SetFillColor(3);
+  Char_t CutName[256];
+  TLine *l1 = pave->AddLine(0,0.89,1,0.89);
+  l1->SetLineWidth(2);
+  TLine *l2 = pave->AddLine(0,0.79,1,0.79);
+  l2->SetLineWidth(2);
+  TLine *l3 = pave->AddLine(0,0.69,1,0.69);
+  l3->SetLineWidth(2);
+  TLine *l4 = pave->AddLine(0,0.59,1,0.59);
+  l4->SetLineWidth(2);
+  TLine *l5 = pave->AddLine(0,0.49,1,0.49);
+  l5->SetLineWidth(2);
+  TLine *l6 = pave->AddLine(0,0.39,1,0.39);
+  l6->SetLineWidth(2);
+  TLine *l7 = pave->AddLine(0,0.29,1,0.29);
+  l7->SetLineWidth(2);
+  TLine *l8 = pave->AddLine(0,0.19,1,0.19);
+  l8->SetLineWidth(2);
+  TLine *l9 = pave->AddLine(0,0.09,1,0.09);
+  l9->SetLineWidth(2);
+
+  sprintf(CutName,"Total number of tracks: %d",fTotalTracks);
+  TText *t1 = pave->AddText(CutName);
+  t1->SetTextColor(1);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+  sprintf(CutName,"Total number of accepted tracks: %d",fAcceptedTracks);
+  t1 = pave->AddText(CutName);
+  t1->SetTextColor(1);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+  sprintf(CutName,"P range: [%f,%f]",fPMin,fPMax);
+  t1 = pave->AddText(CutName);
+  t1->SetTextColor(1);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+  sprintf(CutName,"Tracks rejected: %f",100.0*fP/fTotalTracks);
+  t1 = pave->AddText(CutName);
+  t1->SetTextColor(1);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+  sprintf(CutName,"Pt range: [%f,%f]",fPtMin,fPtMax);
+  t1 = pave->AddText(CutName);
+  t1->SetTextColor(1);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+  sprintf(CutName,"Tracks rejected: %f",100.0*fPt/fTotalTracks);
+  t1 = pave->AddText(CutName);
+  t1->SetTextColor(1);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+
+  sprintf(CutName,"Px range: [%f,%f]",fPxMin,fPxMax);
+  t1 = pave->AddText(CutName);
+  t1->SetTextColor(1);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+  sprintf(CutName,"Tracks rejected: %f",100.0*fPx/fTotalTracks);
+  t1 = pave->AddText(CutName);
+  t1->SetTextColor(1);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+  sprintf(CutName,"Py range: [%f,%f]",fPyMin,fPyMax);
+  t1 = pave->AddText(CutName);
+  t1->SetTextColor(1);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+  sprintf(CutName,"Tracks rejected: %f",100.0*fPy/fTotalTracks);
+  t1 = pave->AddText(CutName);
+  t1->SetTextColor(1);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+  sprintf(CutName,"Pz range: [%f,%f]",fPzMin,fPzMax);
+  t1 = pave->AddText(CutName);
+  t1->SetTextColor(1);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+  sprintf(CutName,"Tracks rejected: %f",100.0*fPz/fTotalTracks);
+  t1 = pave->AddText(CutName);
+  t1->SetTextColor(1);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+  
+  sprintf(CutName,"br range: [%f,%f]",fBrMin,fBrMax);
+  t1 = pave->AddText(CutName);
+  t1->SetTextColor(1);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+  sprintf(CutName,"Tracks rejected: %f",100.0*fbr/fTotalTracks);
+  t1 = pave->AddText(CutName);
+  t1->SetTextColor(1);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+  sprintf(CutName,"bz range: [%f,%f]",fBzMin,fBzMax);
+  t1 = pave->AddText(CutName);
+  t1->SetTextColor(1);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+  sprintf(CutName,"Tracks rejected: %f",100.0*fbz/fTotalTracks);
+  t1 = pave->AddText(CutName);
+  t1->SetTextColor(1);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+
+  sprintf(CutName,"eta range: [%f,%f]",fEtaMin,fEtaMax);
+  t1 = pave->AddText(CutName);
+  t1->SetTextColor(1);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+  sprintf(CutName,"Tracks rejected: %f",100.0*fEta/fTotalTracks);
+  t1 = pave->AddText(CutName);
+  t1->SetTextColor(1);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+
+  sprintf(CutName,"y range: [%f,%f]",fRapMin,fRapMax);
+  t1 = pave->AddText(CutName);
+  t1->SetTextColor(1);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+  sprintf(CutName,"Tracks rejected: %f",100.0*fRap/fTotalTracks);
+  t1 = pave->AddText(CutName);
+  t1->SetTextColor(1);
+  t1->SetTextSize(0.04);
+  t1->SetTextAlign(11);
+
+  return pave;
+}
+
+//----------------------------------------//
+void AliAnalysisTrackCuts::PrintTrackCuts()
+{
+  //GetTrackCuts()->Draw();
+
+  AliInfo(Form("**************TRACK CUTS**************"));
+  GetTrackStats();
+  if(fFlagP)
+    GetPStats();
+  if(fFlagPt)
+    GetPtStats();
+  if(fFlagPx)
+    GetPxStats();
+  if(fFlagPy)
+    GetPyStats();
+  if(fFlagPz)
+    GetPzStats();
+  if(fFlagEta)
+    GetEtaStats();
+  if(fFlagRap)
+    GetRapStats();
+  if(fFlagbr)
+    GetBrStats();
+  if(fFlagbz)
+    GetBzStats(); 
+  AliInfo(Form("**************************************"));
+}
diff --git a/ANALYSIS/AliAnalysisTrackCuts.h b/ANALYSIS/AliAnalysisTrackCuts.h
new file mode 100644 (file)
index 0000000..65541a0
--- /dev/null
@@ -0,0 +1,92 @@
+#ifndef ALIANALYSISTRACKCUTS_H
+#define ALIANALYSISTRACKCUTS_H
+/*  See cxx source for full Copyright notice */
+
+
+/* $Id$ */
+
+//-------------------------------------------------------------------------
+//                       Class AliAnalysisTrackCuts
+//   This is the class for the cuts in event & track level
+//
+//    Origin: Panos Christakoglou, UOA-CERN, Panos.Christakoglou@cern.ch
+//-------------------------------------------------------------------------
+
+#include <TPaveText.h>
+#include <TROOT.h>
+#include <TObject.h>
+
+#include "AliESD.h"
+#include "AliESDtrack.h"
+
+class AliAnalysisTrackCuts : public TObject
+{
+ private:
+  Float_t fPMin, fPMax;  //Definition of the range of the P
+  Float_t fPtMin, fPtMax;  //Definition of the range of the Pt
+  Float_t fPxMin, fPxMax;  //Definition of the range of the Px
+  Float_t fPyMin, fPyMax;  //Definition of the range of the Py
+  Float_t fPzMin, fPzMax;  //Definition of the range of the Pz
+  Float_t fEtaMin, fEtaMax;  //Definition of the range of the eta
+  Float_t fRapMin, fRapMax;  //Definition of the range of the y
+  Float_t fBrMin, fBrMax;  //Definition of the range of the br
+  Float_t fBzMin, fBzMax;  //Definition of the range of the bz
+
+  Int_t fP;  //Number of events rejected due to P cut
+  Int_t fPt;  //Number of events rejected due to Pt cut
+  Int_t fPx;  //Number of events rejected due to Px cut
+  Int_t fPy;  //Number of events rejected due to Py cut
+  Int_t fPz;  //Number of events rejected due to Pz cut
+  Int_t fEta;  //Number of events rejected due to eta cut
+  Int_t fRap;  //Number of events rejected due to y cut
+  Int_t fbr;  //Number of events rejected due to br cut
+  Int_t fbz;  //Number of events rejected due to bz cut
+  Int_t fTotalTracks;  //Total number of tracks
+  Int_t fAcceptedTracks;  //Total number of accepted tracks
+   
+  Int_t fFlagP;  //Flag that shows if the P cut was imposed
+  Int_t fFlagPt;  //Flag that shows if the Pt cut was imposed
+  Int_t fFlagPx;  //Flag that shows if the Px cut was imposed
+  Int_t fFlagPy;  //Flag that shows if the Py cut was imposed
+  Int_t fFlagPz;  //Flag that shows if the Pz cut was imposed
+  Int_t fFlagEta;  //Flag that shows if the eta cut was imposed
+  Int_t fFlagRap;  //Flag that shows if the y cut was imposed
+  Int_t fFlagbr;  //Flag that shows if the br cut was imposed
+  Int_t fFlagbz;  //Flag that shows if the bz cut was imposed
+ public:
+  AliAnalysisTrackCuts();
+  
+  ~AliAnalysisTrackCuts();
+
+  void Reset();
+  
+  void SetPRange(Float_t r1, Float_t r2);
+  void SetPtRange(Float_t r1, Float_t r2);
+  void SetPxRange(Float_t r1, Float_t r2);
+  void SetPyRange(Float_t r1, Float_t r2);
+  void SetPzRange(Float_t r1, Float_t r2);
+  void SetBrRange(Float_t r1, Float_t r2);
+  void SetBzRange(Float_t r1, Float_t r2);
+  void SetEtaRange(Float_t r1, Float_t r2);
+  void SetRapRange(Float_t r1, Float_t r2);
+  
+  Bool_t IsAccepted(AliESD *esd,AliESDtrack *esdtrack);
+
+  TPaveText *GetTrackCuts();
+  void PrintTrackCuts();
+  void GetTrackStats();
+  void GetPStats();
+  void GetPxStats();
+  void GetPyStats();
+  void GetPzStats();
+  void GetPtStats();
+  void GetEtaStats();
+  void GetRapStats();
+  void GetBrStats();
+  void GetBzStats();
+  
+  ClassDef(AliAnalysisTrackCuts, 1)
+} ;
+
+#endif
diff --git a/ANALYSIS/AliBalance.cxx b/ANALYSIS/AliBalance.cxx
new file mode 100644 (file)
index 0000000..97da8ee
--- /dev/null
@@ -0,0 +1,393 @@
+/**************************************************************************
+ * Author: Panos Christakoglou.                                           *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+//-----------------------------------------------------------------
+//           Balance Function class
+//   This is the class to deal with the Balance Function analysis
+//   Origin: Panos Christakoglou, UOA-CERN, Panos.Christakoglou@cern.ch
+//-----------------------------------------------------------------
+
+
+#include <stdlib.h>
+
+//ROOT
+#include <Riostream.h>
+#include <TROOT.h>
+#include <TObject.h>
+#include <TSystem.h>
+#include <TObject.h>
+#include <TMath.h>
+#include <TLorentzVector.h>
+
+#include "AliBalance.h"
+
+ClassImp(AliBalance)
+
+//----------------------------------------//
+AliBalance::AliBalance()
+{
+  for(Int_t i = 0; i < MAXIMUM_NUMBER_OF_STEPS; i++)
+    {
+      fNpp[i] = .0;
+      fNnn[i] = .0;
+      fNpn[i] = .0;
+      fB[i] = 0.0;
+      ferror[i] = 0.0;
+    } 
+  fAnalyzedEvents = 0;
+}
+
+//----------------------------------------//
+AliBalance::AliBalance(Double_t P2_Start, Double_t P2_Stop, Int_t P2_bins)
+{
+  this->fP2_Start = P2_Start;
+  this->fP2_Stop = P2_Stop;
+  this->fNumberOfBins = P2_bins;
+  this->fP2_Step = fabs(fP2_Start - fP2_Stop) / (Double_t)fNumberOfBins;
+}
+
+//----------------------------------------//
+AliBalance::~AliBalance()
+{
+}
+
+//----------------------------------------//
+void AliBalance::SetNumberOfBins(Int_t ibins)
+{
+  this->fNumberOfBins = ibins ;
+}
+
+//----------------------------------------//
+void AliBalance::SetInterval(Double_t P2_Start, Double_t P2_Stop)
+{
+  this->fP2_Start = P2_Start;
+  this->fP2_Stop = P2_Stop;
+  this->fP2_Step = fabs(P2_Start - P2_Stop) / (Double_t)fNumberOfBins;
+}
+
+//----------------------------------------//
+void AliBalance::SetAnalysisType(Int_t iType)
+{
+  //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
+  this->fAnalysisType = iType; 
+  if(fAnalysisType==0)
+    {
+      cout<<" ====================== "<<endl;
+      cout<<"||Analysis selected: y||"<<endl;
+      cout<<" ====================== "<<endl;
+    } 
+  else if(fAnalysisType==1)
+    {
+      cout<<" ======================== "<<endl;
+      cout<<"||Analysis selected: eta||"<<endl;
+      cout<<" ======================== "<<endl;
+    }
+  else if(fAnalysisType==2)
+    {
+      cout<<" ========================== "<<endl;
+      cout<<"||Analysis selected: Qlong||"<<endl;
+      cout<<" ========================== "<<endl;
+    }
+  else if(fAnalysisType==3)
+    {
+      cout<<" ========================= "<<endl;
+      cout<<"||Analysis selected: Qout||"<<endl;
+      cout<<" ========================= "<<endl;
+    }
+  else if(fAnalysisType==4)
+    {
+      cout<<" ========================== "<<endl;
+      cout<<"||Analysis selected: Qside||"<<endl;
+      cout<<" ========================== "<<endl;
+    }
+  else if(fAnalysisType==5)
+    {
+      cout<<" ========================= "<<endl;
+      cout<<"||Analysis selected: Qinv||"<<endl;
+      cout<<" ========================= "<<endl;
+    }
+  else if(fAnalysisType==6)
+    {
+      cout<<" ======================== "<<endl;
+      cout<<"||Analysis selected: phi||"<<endl;
+      cout<<" ======================== "<<endl;
+    }
+  else
+    {
+      cout<<"Selection of analysis mode failed!!!"<<endl;
+      cout<<"Choices are: 0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi"<<endl;
+      abort();
+    }
+}
+
+//----------------------------------------//
+void AliBalance::SetParticles(TLorentzVector *P, Double_t *charge, Int_t dim)
+{
+  this->fV = P;
+  this->fCharge = charge;
+  fNtrack = dim;
+}
+
+
+//----------------------------------------//
+void AliBalance::CalculateBalance()
+{
+  fAnalyzedEvents++;
+  Int_t i = 0 , j = 0;
+  Int_t ibin = 0;
+
+  for(i = 0; i < fNtrack; i++)
+    {
+      if(fCharge[i] > 0)
+       fNp += 1.;
+      if(fCharge[i] < 0)
+       fNn += 1.;
+    }
+
+  //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
+   if(fAnalysisType==0)
+    {
+      for(i = 1; i < fNtrack; i++)
+       {
+         for(j = 0; j < i; j++)
+           {
+             Double_t rap1 = 0.5*log((fV[i].E() - fV[i].Pz())/(fV[i].E() + fV[i].Pz())); 
+             Double_t rap2 = 0.5*log((fV[j].E() - fV[j].Pz())/(fV[j].E() + fV[j].Pz())); 
+             Double_t dy = TMath::Abs(rap1 - rap2);
+             ibin = Int_t(dy/fP2_Step);
+             if((fCharge[i] > 0)&&(fCharge[j] > 0))
+               fNpp[ibin] += 1.;
+             if((fCharge[i] < 0)&&(fCharge[j] < 0))
+               fNnn[ibin] += 1.;
+             if((fCharge[i] > 0)&&(fCharge[j] < 0))
+               fNpn[ibin] += 1.;
+             if((fCharge[i] < 0)&&(fCharge[j] > 0))
+               fNpn[ibin] += 1.;
+           }
+       }
+    }//case 0
+   if(fAnalysisType==1)
+    {
+      for(i = 1; i < fNtrack; i++)
+       {
+         for(j = 0; j < i; j++)
+           {
+             Double_t P1 = sqrt(pow(fV[i].Px(),2) + pow(fV[i].Py(),2) + pow(fV[i].Pz(),2)); 
+             Double_t P2 = sqrt(pow(fV[j].Px(),2) + pow(fV[j].Py(),2) + pow(fV[j].Pz(),2));
+             Double_t eta1 = 0.5*log((P1 - fV[i].Pz())/(P1 + fV[i].Pz())); 
+             Double_t eta2 = 0.5*log((P2 - fV[j].Pz())/(P2 + fV[j].Pz())); 
+             Double_t deta = TMath::Abs(eta1 - eta2);
+             ibin = Int_t(deta/fP2_Step);
+             if((fCharge[i] > 0)&&(fCharge[j] > 0))
+               fNpp[ibin] += 1.;
+             if((fCharge[i] < 0)&&(fCharge[j] < 0))
+               fNnn[ibin] += 1.;
+             if((fCharge[i] > 0)&&(fCharge[j] < 0))
+               fNpn[ibin] += 1.;
+             if((fCharge[i] < 0)&&(fCharge[j] > 0))
+               fNpn[ibin] += 1.;
+           }
+       }
+    }//case 1
+   if(fAnalysisType==2)
+    {
+      for(i = 1; i < fNtrack; i++)
+       {
+         for(j = 0; j < i; j++)
+           {
+             Double_t ETot = fV[i].E() + fV[j].E();
+             Double_t PxTot = fV[i].Px() + fV[j].Px();
+             Double_t PyTot = fV[i].Py() + fV[j].Py();
+             Double_t PzTot = fV[i].Pz() + fV[j].Pz();
+             Double_t Q0Tot = fV[i].E() - fV[j].E();
+             Double_t QzTot = fV[i].Pz() - fV[j].Pz();
+             Double_t Snn = pow(ETot,2) - pow(PxTot,2) - pow(PyTot,2) - pow(PzTot,2);
+             Double_t PtTot = sqrt( pow(PxTot,2) + pow(PyTot,2));
+             Double_t Qlong = TMath::Abs(ETot*QzTot - PzTot*Q0Tot)/sqrt(Snn + pow(PtTot,2));
+             ibin = Int_t(Qlong/fP2_Step);
+             //cout<<ibin<<endl;
+             if((fCharge[i] > 0)&&(fCharge[j] > 0))
+               fNpp[ibin] += 1.;
+             if((fCharge[i] < 0)&&(fCharge[j] < 0))
+               fNnn[ibin] += 1.;
+             if((fCharge[i] > 0)&&(fCharge[j] < 0))
+               fNpn[ibin] += 1.;
+             if((fCharge[i] < 0)&&(fCharge[j] > 0))
+               fNpn[ibin] += 1.;
+           }
+       }
+    }//case 2
+  if(fAnalysisType==3)
+    {
+      for(i = 1; i < fNtrack; i++)
+       {
+         for(j = 0; j < i; j++)
+           {
+             Double_t ETot = fV[i].E() + fV[j].E();
+             Double_t PxTot = fV[i].Px() + fV[j].Px();
+             Double_t PyTot = fV[i].Py() + fV[j].Py();
+             Double_t PzTot = fV[i].Pz() + fV[j].Pz();
+             Double_t QxTot = fV[i].Px() - fV[j].Px();
+             Double_t QyTot = fV[i].Py() - fV[j].Py();
+             Double_t Snn = pow(ETot,2) - pow(PxTot,2) - pow(PyTot,2) - pow(PzTot,2);
+             Double_t PtTot = sqrt( pow(PxTot,2) + pow(PyTot,2));
+             Double_t Qout = sqrt(Snn/(Snn + pow(PtTot,2))) * TMath::Abs(PxTot*QxTot + PyTot*QyTot)/PtTot;
+             ibin = Int_t(Qout/fP2_Step);
+             //cout<<ibin<<endl;
+             if((fCharge[i] > 0)&&(fCharge[j] > 0))
+               fNpp[ibin] += 1.;
+             if((fCharge[i] < 0)&&(fCharge[j] < 0))
+               fNnn[ibin] += 1.;
+             if((fCharge[i] > 0)&&(fCharge[j] < 0))
+               fNpn[ibin] += 1.;
+             if((fCharge[i] < 0)&&(fCharge[j] > 0))
+               fNpn[ibin] += 1.;
+           }
+       }
+    }//case 3
+  if(fAnalysisType==4)
+    {
+      for(i = 1; i < fNtrack; i++)
+       {
+         for(j = 0; j < i; j++)
+           {
+             Double_t PxTot = fV[i].Px() + fV[j].Px();
+             Double_t PyTot = fV[i].Py() + fV[j].Py();
+             Double_t QxTot = fV[i].Px() - fV[j].Px();
+             Double_t QyTot = fV[i].Py() - fV[j].Py();
+             Double_t PtTot = sqrt( pow(PxTot,2) + pow(PyTot,2));
+             Double_t Qside = TMath::Abs(PxTot*QyTot - PyTot*QxTot)/PtTot;
+             ibin = Int_t(Qside/fP2_Step);
+             //cout<<ibin<<endl;
+             if((fCharge[i] > 0)&&(fCharge[j] > 0))
+               fNpp[ibin] += 1.;
+             if((fCharge[i] < 0)&&(fCharge[j] < 0))
+               fNnn[ibin] += 1.;
+             if((fCharge[i] > 0)&&(fCharge[j] < 0))
+               fNpn[ibin] += 1.;
+             if((fCharge[i] < 0)&&(fCharge[j] > 0))
+               fNpn[ibin] += 1.;
+           }
+       }
+    }//case 4
+   if(fAnalysisType==5)
+    {
+      for(i = 1; i < fNtrack; i++)
+       {
+         for(j = 0; j < i; j++)
+           {
+             Double_t Q0Tot = fV[i].E() - fV[j].E();
+             Double_t QxTot = fV[i].Px() - fV[j].Px();
+             Double_t QyTot = fV[i].Py() - fV[j].Py();
+             Double_t QzTot = fV[i].Pz() - fV[j].Pz();
+             Double_t Qinv = sqrt(TMath::Abs(-pow(Q0Tot,2) +pow(QxTot,2) +pow(QyTot,2) +pow(QzTot,2)));
+             ibin = Int_t(Qinv/fP2_Step);
+             //cout<<ibin<<endl;
+             if((fCharge[i] > 0)&&(fCharge[j] > 0))
+               fNpp[ibin] += 1.;
+             if((fCharge[i] < 0)&&(fCharge[j] < 0))
+               fNnn[ibin] += 1.;
+             if((fCharge[i] > 0)&&(fCharge[j] < 0))
+               fNpn[ibin] += 1.;
+             if((fCharge[i] < 0)&&(fCharge[j] > 0))
+               fNpn[ibin] += 1.;
+           }
+       }
+    }//case 5  
+  if(fAnalysisType==6)
+    {
+      for(i = 1; i < fNtrack; i++)
+       {
+         for(j = 0; j < i; j++)
+           {
+             Double_t phi1 = TMath::ATan(fV[i].Py()/fV[i].Px())*180.0/TMath::Pi();
+             Double_t phi2 = TMath::ATan(fV[j].Py()/fV[j].Px())*180.0/TMath::Pi();
+             Double_t dphi = TMath::Abs(phi1 - phi2);
+             ibin = Int_t(dphi/fP2_Step);
+             if((fCharge[i] > 0)&&(fCharge[j] > 0))
+               fNpp[ibin] += 1.;
+             if((fCharge[i] < 0)&&(fCharge[j] < 0))
+               fNnn[ibin] += 1.;
+             if((fCharge[i] > 0)&&(fCharge[j] < 0))
+               fNpn[ibin] += 1.;
+             if((fCharge[i] < 0)&&(fCharge[j] > 0))
+               fNpn[ibin] += 1.;
+           }
+       }
+    }//case 6
+}
+
+//----------------------------------------//
+Double_t AliBalance::GetBalance(Int_t p2)
+{
+  fB[p2] = 0.5*(((fNpn[p2] - 2.0*fNnn[p2])/fNn) + ((fNpn[p2] - 2.0*fNpp[p2])/fNp))/fP2_Step;
+
+  return fB[p2];
+}
+    
+//----------------------------------------//
+Double_t AliBalance::GetError(Int_t p2)
+{
+  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;
+
+  return ferror[p2];
+}
+
+//----------------------------------------//
+void AliBalance::PrintResults()
+{
+  Double_t x[MAXIMUM_NUMBER_OF_STEPS];
+  Double_t fSumXi = 0.0, fSumBi = 0.0, fSumBiXi = 0.0;
+  Double_t fSumBiXi2 = 0.0, fSumBi2Xi2 = 0.0;
+  Double_t fSumDeltaBi2 = 0.0, fSumXi2DeltaBi2 = 0.0;
+  Double_t delta_bal_P2 = 0.0, Integral = 0.0;
+  Double_t DeltaErrorNew = 0.0;
+
+  cout<<"=================================================="<<endl;
+  for(Int_t i = 0; i < fNumberOfBins; i++)
+    { 
+      Double_t x = fP2_Start + fP2_Step*i + fP2_Step/2 ;
+      cout<<"B: "<<fB[i]<<"\t Error: "<<ferror[i]<<"\t bin: "<<x<<endl;
+    } 
+  cout<<"=================================================="<<endl;
+  for(Int_t i = 1; i < fNumberOfBins; i++)
+    {
+      fSumXi += x[i];
+      fSumBi += fB[i];
+      fSumBiXi += fB[i]*x[i];
+      fSumBiXi2 += fB[i]*pow(x[i],2);
+      fSumBi2Xi2 += pow(fB[i],2)*pow(x[i],2);
+      fSumDeltaBi2 +=  pow(ferror[i],2) ;
+      fSumXi2DeltaBi2 += pow(x[i],2) * pow(ferror[i],2) ;
+      
+      delta_bal_P2 += fP2_Step*pow(ferror[i],2) ;
+      Integral += fP2_Step*fB[i] ;
+    }
+  for(Int_t i = 1; i < fNumberOfBins; i++)
+    {
+      DeltaErrorNew += ferror[i]*(x[i]*fSumBi - fSumBiXi)/pow(fSumBi,2);
+    }
+  Double_t IntegralError = sqrt(delta_bal_P2) ;
+  
+  Double_t Delta = fSumBiXi / fSumBi ;
+  Double_t DeltaError = (fSumBiXi / fSumBi) * sqrt(pow((sqrt(fSumXi2DeltaBi2)/fSumBiXi),2) + pow((fSumDeltaBi2/fSumBi),2) ) ;
+  cout<<"Analyzed events: "<<fAnalyzedEvents<<endl;
+  cout<<"Width: "<<Delta<<"\t Error: "<<DeltaError<<endl;
+  cout<<"New error: "<<DeltaErrorNew<<endl;
+  cout<<"Interval: "<<Integral<<"\t Error: "<<IntegralError<<endl;
+  cout<<"=================================================="<<endl;
+}
+  
diff --git a/ANALYSIS/AliBalance.h b/ANALYSIS/AliBalance.h
new file mode 100644 (file)
index 0000000..06b4180
--- /dev/null
@@ -0,0 +1,71 @@
+#ifndef ALIBALANCE_H
+#define ALIBALANCE_H
+/*  See cxx source for full Copyright notice */
+
+
+/* $Id$ */
+
+//-------------------------------------------------------------------------
+//                          Class AliBalance
+//   This is the class for the Balance Function analysis
+//
+//    Origin: Panos Christakoglou, UOA-CERN, Panos.Christakoglou@cern.ch
+//-------------------------------------------------------------------------
+
+#include <TROOT.h>
+#include <TObject.h>
+#include <TLorentzVector.h>
+
+#define MAXIMUM_NUMBER_OF_STEPS        1024
+
+class AliBalance : public TObject
+{
+ private:
+  Double_t *fCharge; //matrix with the charge of each track
+  Int_t fNtrack; //number of tracks to compute the BF
+
+  TLorentzVector *fV; //4-momentum vector - info for tracks used to compute the BF
+  
+  Int_t fNumberOfBins; //number of bins of the analyzed interval
+  Int_t fAnalysisType; //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
+  Int_t fAnalyzedEvents; //number of events that have been analyzed
+  Double_t fP2_Start, fP2_Stop, fP2_Step; //inerval info
+       
+  Double_t fNnn[MAXIMUM_NUMBER_OF_STEPS]; //N(--)
+  Double_t fNpp[MAXIMUM_NUMBER_OF_STEPS]; //N(++)
+  Double_t fNpn[MAXIMUM_NUMBER_OF_STEPS]; //N(+-)
+  Double_t fNn, fNp; //number of pos./neg. inside the analyzed interval
+  
+  Double_t fB[MAXIMUM_NUMBER_OF_STEPS]; //BF matrix
+  Double_t ferror[MAXIMUM_NUMBER_OF_STEPS]; //error of the BF
+  
+ public:
+  AliBalance();
+  
+  AliBalance(Double_t P2_Start, Double_t P2_Stop, Int_t P2_Steps);
+  
+  ~AliBalance();
+  
+  void SetParticles(TLorentzVector *P, Double_t *charge, Int_t dim);
+  void SetNumberOfBins(Int_t ibins);
+  void SetAnalysisType(Int_t iType);
+  void SetInterval(Double_t P2_Start, Double_t P2_Stop);
+
+  void CalculateBalance() ;
+  
+  Double_t GetNp() { return 1.0*fNp; }
+  Double_t GetNn() { return 1.0*fNn; }
+  Double_t GetNnn(Int_t p2) { return 1.0*fNnn[p2]; }
+  Double_t GetNpp(Int_t p2) { return 1.0*fNpp[p2]; }
+  Double_t GetNpn(Int_t p2) { return 1.0*fNpn[p2]; }
+  Double_t GetBalance(Int_t p2);
+  Double_t GetError(Int_t p2);
+
+  void PrintResults();
+  ClassDef(AliBalance, 1)
+} ;
+
+#endif
index 656b697..eed8ee2 100644 (file)
@@ -10,7 +10,8 @@ SRCS= TGliteXmlEventlist.cxx\
       AliReaderESD.cxx AliReaderESDTree.cxx \
       AliTrackPoints.cxx AliClusterMap.cxx \
       AliD0toKpi.cxx  AliD0toKpiAnalysis.cxx AliFlowAnalysis.cxx \
-      AliMuonAnalysis.cxx  AliAODv0.cxx
+      AliMuonAnalysis.cxx  AliAODv0.cxx \
+      AliAnalysisEventCuts.cxx AliAnalysisTrackCuts.cxx AliBalance.cxx
 
 HDRS= $(SRCS:.cxx=.h) 
 
@@ -26,6 +27,7 @@ EXPORT:=AliAOD.h AliEventBuffer.h\
       AliReader.h           AliReaderESD.h    \
       AliTrackPoints.h      AliClusterMap.h   \
       AliFlowAnalysis.h     AliReaderESDTree.h \
-      AliMuonAnalysis.h     AliAODv0.h
+      AliMuonAnalysis.h     AliAODv0.h \
+      AliAnalysisEventCuts.h AliAnalysisTrackCuts.h AliBalance.h
 
 EINCLUDE:= TPC ITS