]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG2/EBYE/AliBalance.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
Adding the SetOwner
[u/mrichter/AliRoot.git] / PWG2 / EBYE / AliBalance.cxx
diff --git a/PWG2/EBYE/AliBalance.cxx b/PWG2/EBYE/AliBalance.cxx
index 3c0882b90e07e3c9999782e72719d04515599800..27789919311221206d7be5f3bb37ee47e95a144a 100644 (file)
--- a/PWG2/EBYE/AliBalance.cxx
+++ b/PWG2/EBYE/AliBalance.cxx
@@ -28,6 +28,7 @@
 #include <TObjArray.h>
 #include <TGraphErrors.h>
 #include <TString.h>
+#include <TH1F.h>
 
 #include "AliVParticle.h"
 #include "AliMCParticle.h"
@@ -43,7 +44,13 @@ AliBalance::AliBalance() :
   TObject(), 
   fAnalysisLevel("ESD"), fNumberOfBins(0),
   fAnalysisType(0), fAnalyzedEvents(0), fP2Start(0),
-  fP2Stop(0), fP2Step(0), fNn(0), fNp(0) {
+  fP2Stop(0), fP2Step(0), fNn(0), fNp(0),
+  fHistfNnn(new TH1F("fHistfNnn","(--) component;;Entries",
+                    fNumberOfBins,fP2Start,fP2Stop)),
+  fHistfNpp(new TH1F("fHistfNpp","(++) component;;Entries",
+                    fNumberOfBins,fP2Start,fP2Stop)),
+  fHistfNpn(new TH1F("fHistfNpn","(+-) component;;Entries",
+                    fNumberOfBins,fP2Start,fP2Stop)) {
   // Default constructor
   for(Int_t i = 0; i < MAXIMUM_NUMBER_OF_STEPS; i++) {
     fNpp[i] = .0;
@@ -52,6 +59,46 @@ AliBalance::AliBalance() :
     fB[i] = 0.0;
     ferror[i] = 0.0;
   } 
+
+  switch(fAnalysisType) {
+  case 0:
+    fHistfNnn->GetXaxis()->SetTitle("#Delta y");
+    fHistfNpp->GetXaxis()->SetTitle("#Delta y");
+    fHistfNpn->GetXaxis()->SetTitle("#Delta y");
+    break;
+  case 1:
+    fHistfNnn->GetXaxis()->SetTitle("#Delta #eta");
+    fHistfNpp->GetXaxis()->SetTitle("#Delta #eta");
+    fHistfNpn->GetXaxis()->SetTitle("#Delta #eta");
+    break;
+  case 2:
+    fHistfNnn->GetXaxis()->SetTitle("q_{long} (GeV/c)");
+    fHistfNpp->GetXaxis()->SetTitle("q_{long} (GeV/c)");
+    fHistfNpn->GetXaxis()->SetTitle("q_{long} (GeV/c)");
+    break;
+  case 3:
+    fHistfNnn->GetXaxis()->SetTitle("q_{out} (GeV/c)");
+    fHistfNpp->GetXaxis()->SetTitle("q_{out} (GeV/c)");
+    fHistfNpn->GetXaxis()->SetTitle("q_{out} (GeV/c)");
+    break;
+  case 4:
+    fHistfNnn->GetXaxis()->SetTitle("q_{side} (GeV/c)");
+    fHistfNpp->GetXaxis()->SetTitle("q_{side} (GeV/c)");
+    fHistfNpn->GetXaxis()->SetTitle("q_{side} (GeV/c)");
+    break;
+  case 5:
+    fHistfNnn->GetXaxis()->SetTitle("q_{inv.} (GeV/c)");
+    fHistfNpp->GetXaxis()->SetTitle("q_{inv.} (GeV/c)");
+    fHistfNpn->GetXaxis()->SetTitle("q_{inv.} (GeV/c)");
+    break;
+  case 6:
+    fHistfNnn->GetXaxis()->SetTitle("#Delta #phi");
+    fHistfNpp->GetXaxis()->SetTitle("#Delta #phi");
+    fHistfNpn->GetXaxis()->SetTitle("#Delta #phi");
+    break;
+  default:
+    break;
+  }
 }
 
 //____________________________________________________________________//
@@ -60,7 +107,13 @@ AliBalance::AliBalance(Double_t p2Start, Double_t p2Stop, Int_t p2Bins) :
   fNumberOfBins(p2Bins), fAnalysisType(0), 
   fAnalyzedEvents(0), fP2Start(p2Start), fP2Stop(p2Stop), 
   fP2Step(TMath::Abs(fP2Start - fP2Stop) / (Double_t)fNumberOfBins), 
-  fNn(0), fNp(0) {
+  fNn(0), fNp(0),
+  fHistfNnn(new TH1F("fHistfNnn","(--) component;;Entries",
+                    fNumberOfBins,fP2Start,fP2Stop)),
+  fHistfNpp(new TH1F("fHistfNpp","(++) component;;Entries",
+                    fNumberOfBins,fP2Start,fP2Stop)),
+  fHistfNpn(new TH1F("fHistfNpn","(+-) component;;Entries",
+                    fNumberOfBins,fP2Start,fP2Stop)) {
   // Constructor
   for(Int_t i = 0; i < MAXIMUM_NUMBER_OF_STEPS; i++) {
     fNpp[i] = .0;
@@ -69,6 +122,46 @@ AliBalance::AliBalance(Double_t p2Start, Double_t p2Stop, Int_t p2Bins) :
     fB[i] = 0.0;
     ferror[i] = 0.0;
   } 
+
+  switch(fAnalysisType) {
+  case 0:
+    fHistfNnn->GetXaxis()->SetTitle("#Delta y");
+    fHistfNpp->GetXaxis()->SetTitle("#Delta y");
+    fHistfNpn->GetXaxis()->SetTitle("#Delta y");
+    break;
+  case 1:
+    fHistfNnn->GetXaxis()->SetTitle("#Delta #eta");
+    fHistfNpp->GetXaxis()->SetTitle("#Delta #eta");
+    fHistfNpn->GetXaxis()->SetTitle("#Delta #eta");
+    break;
+  case 2:
+    fHistfNnn->GetXaxis()->SetTitle("q_{long} (GeV/c)");
+    fHistfNpp->GetXaxis()->SetTitle("q_{long} (GeV/c)");
+    fHistfNpn->GetXaxis()->SetTitle("q_{long} (GeV/c)");
+    break;
+  case 3:
+    fHistfNnn->GetXaxis()->SetTitle("q_{out} (GeV/c)");
+    fHistfNpp->GetXaxis()->SetTitle("q_{out} (GeV/c)");
+    fHistfNpn->GetXaxis()->SetTitle("q_{out} (GeV/c)");
+    break;
+  case 4:
+    fHistfNnn->GetXaxis()->SetTitle("q_{side} (GeV/c)");
+    fHistfNpp->GetXaxis()->SetTitle("q_{side} (GeV/c)");
+    fHistfNpn->GetXaxis()->SetTitle("q_{side} (GeV/c)");
+    break;
+  case 5:
+    fHistfNnn->GetXaxis()->SetTitle("q_{inv.} (GeV/c)");
+    fHistfNpp->GetXaxis()->SetTitle("q_{inv.} (GeV/c)");
+    fHistfNpn->GetXaxis()->SetTitle("q_{inv.} (GeV/c)");
+    break;
+  case 6:
+    fHistfNnn->GetXaxis()->SetTitle("#Delta #phi");
+    fHistfNpp->GetXaxis()->SetTitle("#Delta #phi");
+    fHistfNpn->GetXaxis()->SetTitle("#Delta #phi");
+    break;
+  default:
+    break;
+  }
 }
 
 //____________________________________________________________________//
@@ -81,7 +174,10 @@ AliBalance::AliBalance(const AliBalance& balance):
   fP2Stop(balance.fP2Stop),
   fP2Step(balance.fP2Step),
   fNn(balance.fNn),
-  fNp(balance.fNp) {
+  fNp(balance.fNp),
+  fHistfNnn(balance.fHistfNnn), 
+  fHistfNpp(balance.fHistfNpp), 
+  fHistfNpn(balance.fHistfNpn) {
   //copy constructor
   for(Int_t i = 0; i < MAXIMUM_NUMBER_OF_STEPS; i++) {
     fNpp[i] = .0;
@@ -95,6 +191,9 @@ AliBalance::AliBalance(const AliBalance& balance):
 //____________________________________________________________________//
 AliBalance::~AliBalance() {
   // Destructor
+  if(fHistfNnn) delete fHistfNnn;
+  if(fHistfNpp) delete fHistfNpp;
+  if(fHistfNpn) delete fHistfNpn;
 }
 
 //____________________________________________________________________//
@@ -177,7 +276,7 @@ void AliBalance::SetAnalysisType(Int_t iType) {
 }
 
 //____________________________________________________________________//
-const char* AliBalance::GetAnalysisType() {
+void AliBalance::PrintAnalysisSettings() {
   //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
   TString analysisType;
   switch(fAnalysisType) {
@@ -205,12 +304,15 @@ const char* AliBalance::GetAnalysisType() {
   default:
     break;
   }
-  analysisType += "\nInterval: ";
-  analysisType += fP2Start; analysisType += " - "; analysisType += fP2Stop;
-  analysisType += "\nSteps: "; analysisType += fP2Step;
-  analysisType += "\nBins: "; analysisType += fNumberOfBins;
-
-  return analysisType.Data();
+  
+  Printf("======================================");
+  Printf("Analysis level: %s",fAnalysisLevel.Data());
+  Printf("Analysis type: %s",analysisType.Data());
+  Printf("Analyzed interval (min.): %lf",fP2Start);
+  Printf("Analyzed interval (max.): %lf",fP2Stop);
+  Printf("Number of bins: %d",fNumberOfBins);
+  Printf("Step: %lf",fP2Step);
+  Printf("======================================");
 }
 
 //____________________________________________________________________//
@@ -224,6 +326,7 @@ void AliBalance::CalculateBalance(TObjArray *gTrackArray) {
   AliVParticle* track1 = 0;
   AliVParticle* track2 = 0;
     
+  //Printf("(AliBalance) Number of tracks: %d",gTrackArray->GetEntries());
   Int_t gNtrack = gTrackArray->GetEntries();
   for(i = 0; i < gNtrack; i++) {
     if(fAnalysisLevel == "ESD")
@@ -232,9 +335,13 @@ void AliBalance::CalculateBalance(TObjArray *gTrackArray) {
       track = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
     else if(fAnalysisLevel == "MC")
       track = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
-    Short_t charge = track->Charge();
-    if(charge > 0) fNp += 1.;
-    if(charge < 0) fNn += 1.;
+
+    if(track) {
+      Short_t charge = track->Charge();
+      if(charge > 0) fNp += 1.;
+      if(charge < 0) fNn += 1.;
+    }
+    else continue;
   }
   //Printf("Np: %lf - Nn: %lf",fNp,fNn);
 
@@ -247,10 +354,17 @@ void AliBalance::CalculateBalance(TObjArray *gTrackArray) {
        track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
       else if(fAnalysisLevel == "MC")
        track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
-      Short_t charge1 = track1->Charge();
-      Double_t pZ1 = track1->Pz();
-      Double_t energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
-                                    TMath::Power(track1->M(),2));
+
+      Short_t charge1 = 0;
+      Double_t pZ1 = 0., energy1 = 0.;
+      if(track1) {
+       charge1 = track1->Charge();
+       pZ1 = track1->Pz();
+       energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
+                             TMath::Power(track1->M(),2));
+      }
+      else continue;
+
       for(j = 0; j < i; j++) {
        if(fAnalysisLevel == "ESD")
          track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
@@ -258,19 +372,23 @@ void AliBalance::CalculateBalance(TObjArray *gTrackArray) {
          track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
        else if(fAnalysisLevel == "MC")
          track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
-       Short_t charge2 = track2->Charge();
-       Double_t pZ2 = track2->Pz();
-       Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
-                                      TMath::Power(track2->M(),2));
-
-       Double_t rap1 = 0.5*log((energy1 + pZ1)/(energy1 - pZ1)); 
-       Double_t rap2 = 0.5*log((energy2 + pZ2)/(energy2 - pZ2)); 
-       Double_t dy = TMath::Abs(rap1 - rap2);
-       ibin = Int_t(dy/fP2Step);
-       if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
-       if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
-       if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
-       if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
+
+       if(track2) {
+         Short_t charge2 = track2->Charge();
+         Double_t pZ2 = track2->Pz();
+         Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
+                                        TMath::Power(track2->M(),2));
+         
+         Double_t rap1 = 0.5*log((energy1 + pZ1)/(energy1 - pZ1)); 
+         Double_t rap2 = 0.5*log((energy2 + pZ2)/(energy2 - pZ2)); 
+         Double_t dy = TMath::Abs(rap1 - rap2);
+         ibin = Int_t(dy/fP2Step);
+         if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
+         if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
+         if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
+         if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
+       }
+       else continue;
       }
     }
   }//case 0
@@ -282,9 +400,16 @@ void AliBalance::CalculateBalance(TObjArray *gTrackArray) {
        track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
       else if(fAnalysisLevel == "MC")
        track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
-      Short_t charge1 = track1->Charge();
-      Double_t pZ1 = track1->Pz();
-      Double_t p1 = track1->P();
+
+      Short_t charge1 = 0;
+      Double_t pZ1 = 0., p1 = 0.;
+      if(track1) {
+       charge1 = track1->Charge();
+       pZ1 = track1->Pz();
+       p1 = track1->P();
+      }
+      else continue;
+
       for(j = 0; j < i; j++) {
        if(fAnalysisLevel == "ESD")
          track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
@@ -292,19 +417,23 @@ void AliBalance::CalculateBalance(TObjArray *gTrackArray) {
          track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
        else if(fAnalysisLevel == "MC")
          track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
-       Short_t charge2 = track2->Charge();
-       Double_t pZ2 = track2->Pz();
-       Double_t p2 = track2->P();
-       Double_t eta1 = 0.5*log((p1 + pZ1)/(p1 - pZ1)); 
-       Double_t eta2 = 0.5*log((p2 + pZ2)/(p2 - pZ2)); 
-       Double_t deta = TMath::Abs(eta1 - eta2);
-       ibin = Int_t(deta/fP2Step);
-       
-       if((charge1 > 0.)&&(charge2 > 0.)) fNpp[ibin] += 1.;
-       if((charge1 < 0.)&&(charge2 < 0.)) fNnn[ibin] += 1.;
-       if((charge1 > 0.)&&(charge2 < 0.)) fNpn[ibin] += 1.;
-       if((charge1 < 0.)&&(charge2 > 0.)) fNpn[ibin] += 1.;
-       //Printf("charge1: %d - eta1: %lf - charge2: %d - eta2: %lf - deta: %lf - ibin: %d - fNpp: %lf - fNnn: %lf - fNpn: %lf",charge1,eta1,charge2,eta2,deta,ibin,fNpp[ibin],fNnn[ibin],fNpn[ibin]);      
+
+       if(track2) {
+         Short_t charge2 = track2->Charge();
+         Double_t pZ2 = track2->Pz();
+         Double_t p2 = track2->P();
+         Double_t eta1 = 0.5*log((p1 + pZ1)/(p1 - pZ1)); 
+         Double_t eta2 = 0.5*log((p2 + pZ2)/(p2 - pZ2)); 
+         Double_t deta = TMath::Abs(eta1 - eta2);
+         ibin = Int_t(deta/fP2Step);
+         
+         if((charge1 > 0.)&&(charge2 > 0.)) fNpp[ibin] += 1.;
+         if((charge1 < 0.)&&(charge2 < 0.)) fNnn[ibin] += 1.;
+         if((charge1 > 0.)&&(charge2 < 0.)) fNpn[ibin] += 1.;
+         if((charge1 < 0.)&&(charge2 > 0.)) fNpn[ibin] += 1.;
+         //Printf("charge1: %d - eta1: %lf - charge2: %d - eta2: %lf - deta: %lf - ibin: %d - fNpp: %lf - fNnn: %lf - fNpn: %lf",charge1,eta1,charge2,eta2,deta,ibin,fNpp[ibin],fNnn[ibin],fNpn[ibin]);      
+       }
+       else continue;
       }
     }
   }//case 1
@@ -316,12 +445,19 @@ void AliBalance::CalculateBalance(TObjArray *gTrackArray) {
        track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
       else if(fAnalysisLevel == "MC")
        track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
-      Short_t charge1 = track1->Charge();
-      Double_t pX1 = track1->Px();
-      Double_t pY1 = track1->Py();
-      Double_t pZ1 = track1->Pz();
-      Double_t energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
-                                    TMath::Power(track1->M(),2));
+
+      Short_t charge1 = 0;
+      Double_t pX1 = 0., pY1 = 0., pZ1 = 0., energy1 = 0.;
+      if(track1) {
+       charge1 = track1->Charge();
+       pX1 = track1->Px();
+       pY1 = track1->Py();
+       pZ1 = track1->Pz();
+       energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
+                                      TMath::Power(track1->M(),2));
+      }
+      else continue;
+
       for(j = 0; j < i; j++) {
        if(fAnalysisLevel == "ESD")
          track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
@@ -329,27 +465,31 @@ void AliBalance::CalculateBalance(TObjArray *gTrackArray) {
          track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
        else if(fAnalysisLevel == "MC")
          track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
-       Short_t charge2 = track2->Charge();
-       Double_t pX2 = track2->Px();
-       Double_t pY2 = track2->Py();
-       Double_t pZ2 = track2->Pz();
-       Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
-                                      TMath::Power(track2->M(),2));
-       Double_t eTot = energy1 + energy2;
-       Double_t pxTot = pX1 + pX2;
-       Double_t pyTot = pY1 + pY2;
-       Double_t pzTot = pZ1 + pZ2;
-       Double_t q0Tot = energy1 - energy2;
-       Double_t qzTot = pZ1 - pZ2;
-       Double_t snn = TMath::Power(eTot,2) - TMath::Power(pxTot,2) - TMath::Power(pyTot,2) - TMath::Power(pzTot,2);
-       Double_t ptTot = TMath::Sqrt( TMath::Power(pxTot,2) + TMath::Power(pyTot,2));
-       Double_t qLong = TMath::Abs(eTot*qzTot - pzTot*q0Tot)/TMath::Sqrt(snn + TMath::Power(ptTot,2));
-       ibin = Int_t(qLong/fP2Step);
-       //cout<<ibin<<endl;
-       if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
-       if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
-       if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
-       if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
+
+       if(track2) {
+         Short_t charge2 = track2->Charge();
+         Double_t pX2 = track2->Px();
+         Double_t pY2 = track2->Py();
+         Double_t pZ2 = track2->Pz();
+         Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
+                                        TMath::Power(track2->M(),2));
+         Double_t eTot = energy1 + energy2;
+         Double_t pxTot = pX1 + pX2;
+         Double_t pyTot = pY1 + pY2;
+         Double_t pzTot = pZ1 + pZ2;
+         Double_t q0Tot = energy1 - energy2;
+         Double_t qzTot = pZ1 - pZ2;
+         Double_t snn = TMath::Power(eTot,2) - TMath::Power(pxTot,2) - TMath::Power(pyTot,2) - TMath::Power(pzTot,2);
+         Double_t ptTot = TMath::Sqrt( TMath::Power(pxTot,2) + TMath::Power(pyTot,2));
+         Double_t qLong = TMath::Abs(eTot*qzTot - pzTot*q0Tot)/TMath::Sqrt(snn + TMath::Power(ptTot,2));
+         ibin = Int_t(qLong/fP2Step);
+         //cout<<ibin<<endl;
+         if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
+         if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
+         if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
+         if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
+       }
+       else continue;
       }
     }
   }//case 2
@@ -361,12 +501,19 @@ void AliBalance::CalculateBalance(TObjArray *gTrackArray) {
        track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
       else if(fAnalysisLevel == "MC")
        track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
-      Short_t charge1 = track1->Charge();
-      Double_t pX1 = track1->Px();
-      Double_t pY1 = track1->Py();
-      Double_t pZ1 = track1->Pz();
-      Double_t energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
-                                    TMath::Power(track1->M(),2));
+
+      Short_t charge1 = 0;
+      Double_t pX1 = 0., pY1 = 0., pZ1 = 0., energy1 = 0.;
+      if(track1) {
+       charge1 = track1->Charge();
+       pX1 = track1->Px();
+       pY1 = track1->Py();
+       pZ1 = track1->Pz();
+       energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
+                             TMath::Power(track1->M(),2));
+      }
+      else continue;
+
       for(j = 0; j < i; j++) {
        if(fAnalysisLevel == "ESD")
          track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
@@ -374,27 +521,31 @@ void AliBalance::CalculateBalance(TObjArray *gTrackArray) {
          track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
        else if(fAnalysisLevel == "MC")
          track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
-       Short_t charge2 = track2->Charge();
-       Double_t pX2 = track2->Px();
-       Double_t pY2 = track2->Py();
-       Double_t pZ2 = track2->Pz();
-       Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
-                                      TMath::Power(track2->M(),2));
-       Double_t eTot = energy1 + energy2;
-       Double_t pxTot = pX1 + pX2;
-       Double_t pyTot = pY1 + pY2;
-       Double_t pzTot = pZ1 + pZ2;
-       Double_t qxTot = pX1 - pX2;
-       Double_t qyTot = pY1 - pY2;
-       Double_t snn = TMath::Power(eTot,2) - TMath::Power(pxTot,2) - TMath::Power(pyTot,2) - TMath::Power(pzTot,2);
-       Double_t ptTot = TMath::Sqrt( TMath::Power(pxTot,2) + TMath::Power(pyTot,2));
-       Double_t qOut = TMath::Sqrt(snn/(snn + TMath::Power(ptTot,2))) * TMath::Abs(pxTot*qxTot + pyTot*qyTot)/ptTot;
-       ibin = Int_t(qOut/fP2Step);
-       //cout<<ibin<<endl;
-       if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
-       if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
-       if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
-       if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
+
+       if(track2) {
+         Short_t charge2 = track2->Charge();
+         Double_t pX2 = track2->Px();
+         Double_t pY2 = track2->Py();
+         Double_t pZ2 = track2->Pz();
+         Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
+                                        TMath::Power(track2->M(),2));
+         Double_t eTot = energy1 + energy2;
+         Double_t pxTot = pX1 + pX2;
+         Double_t pyTot = pY1 + pY2;
+         Double_t pzTot = pZ1 + pZ2;
+         Double_t qxTot = pX1 - pX2;
+         Double_t qyTot = pY1 - pY2;
+         Double_t snn = TMath::Power(eTot,2) - TMath::Power(pxTot,2) - TMath::Power(pyTot,2) - TMath::Power(pzTot,2);
+         Double_t ptTot = TMath::Sqrt( TMath::Power(pxTot,2) + TMath::Power(pyTot,2));
+         Double_t qOut = TMath::Sqrt(snn/(snn + TMath::Power(ptTot,2))) * TMath::Abs(pxTot*qxTot + pyTot*qyTot)/ptTot;
+         ibin = Int_t(qOut/fP2Step);
+         //cout<<ibin<<endl;
+         if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
+         if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
+         if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
+         if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
+       }
+       else continue;
       }
     }
   }//case 3
@@ -406,11 +557,18 @@ void AliBalance::CalculateBalance(TObjArray *gTrackArray) {
        track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
       else if(fAnalysisLevel == "MC")
        track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
-      Short_t charge1 = track1->Charge();
-      Double_t pX1 = track1->Px();
-      Double_t pY1 = track1->Py();
-      //Double_t energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
-      //TMath::Power(track1->M(),2));
+
+      Short_t charge1 = 0;
+      Double_t pX1 = 0., pY1 = 0.;
+      if(track1) {
+       charge1 = track1->Charge();
+       pX1 = track1->Px();
+       pY1 = track1->Py();
+       //Double_t energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
+       //TMath::Power(track1->M(),2));
+      }
+      else continue;
+
       for(j = 0; j < i; j++) {
        if(fAnalysisLevel == "ESD")
          track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
@@ -418,24 +576,28 @@ void AliBalance::CalculateBalance(TObjArray *gTrackArray) {
          track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
        else if(fAnalysisLevel == "MC")
          track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
-       Short_t charge2 = track2->Charge();
-       Double_t pX2 = track2->Px();
-       Double_t pY2 = track2->Py();
-       //Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
-       //TMath::Power(track2->M(),2));
-       //Double_t eTot = energy1 + energy2;
-       Double_t pxTot = pX1 + pX2;
-       Double_t pyTot = pY1 + pY2;
-       Double_t qxTot = pX1 - pX2;
-       Double_t qyTot = pY1 - pY2;
-       Double_t ptTot = TMath::Sqrt( TMath::Power(pxTot,2) + TMath::Power(pyTot,2));
-       Double_t qSide = TMath::Abs(pxTot*qyTot - pyTot*qxTot)/ptTot;
-       ibin = Int_t(qSide/fP2Step);
-       //cout<<ibin<<endl;
-       if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
-       if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
-       if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
-       if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
+
+       if(track2) {
+         Short_t charge2 = track2->Charge();
+         Double_t pX2 = track2->Px();
+         Double_t pY2 = track2->Py();
+         //Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
+         //TMath::Power(track2->M(),2));
+         //Double_t eTot = energy1 + energy2;
+         Double_t pxTot = pX1 + pX2;
+         Double_t pyTot = pY1 + pY2;
+         Double_t qxTot = pX1 - pX2;
+         Double_t qyTot = pY1 - pY2;
+         Double_t ptTot = TMath::Sqrt( TMath::Power(pxTot,2) + TMath::Power(pyTot,2));
+         Double_t qSide = TMath::Abs(pxTot*qyTot - pyTot*qxTot)/ptTot;
+         ibin = Int_t(qSide/fP2Step);
+         //cout<<ibin<<endl;
+         if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
+         if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
+         if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
+         if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
+       }
+       else continue;
       }
     }
   }//case 4
@@ -447,12 +609,19 @@ void AliBalance::CalculateBalance(TObjArray *gTrackArray) {
        track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
       else if(fAnalysisLevel == "MC")
        track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
-      Short_t charge1 = track1->Charge();
-      Double_t pX1 = track1->Px();
-      Double_t pY1 = track1->Py();
-      Double_t pZ1 = track1->Pz();
-      Double_t energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
-                                    TMath::Power(track1->M(),2));
+
+      Short_t charge1 = 0;
+      Double_t pX1 = 0., pY1 = 0., pZ1 = 0., energy1 = 0.;
+      if(track1) {
+       charge1 = track1->Charge();
+       pX1 = track1->Px();
+       pY1 = track1->Py();
+       pZ1 = track1->Pz();
+       energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
+                             TMath::Power(track1->M(),2));
+      }
+      else continue;
+
       for(j = 0; j < i; j++) {
        if(fAnalysisLevel == "ESD")
          track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
@@ -460,23 +629,27 @@ void AliBalance::CalculateBalance(TObjArray *gTrackArray) {
          track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
        else if(fAnalysisLevel == "MC")
          track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
-       Short_t charge2 = track2->Charge();
-       Double_t pX2 = track2->Px();
-       Double_t pY2 = track2->Py();
-       Double_t pZ2 = track2->Pz();
-       Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
-                                      TMath::Power(track2->M(),2));
-       Double_t q0Tot = energy1 - energy2;
-       Double_t qxTot = pX1 - pX2;
-       Double_t qyTot = pY1 - pY2;
-       Double_t qzTot = pZ1 - pZ2;
-       Double_t qInv = TMath::Sqrt(TMath::Abs(-TMath::Power(q0Tot,2) +TMath::Power(qxTot,2) +TMath::Power(qyTot,2) +TMath::Power(qzTot,2)));
-       ibin = Int_t(qInv/fP2Step);
-       //cout<<ibin<<endl;
-       if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
-       if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
-       if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
-       if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
+
+       if(track2) {
+         Short_t charge2 = track2->Charge();
+         Double_t pX2 = track2->Px();
+         Double_t pY2 = track2->Py();
+         Double_t pZ2 = track2->Pz();
+         Double_t energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
+                                        TMath::Power(track2->M(),2));
+         Double_t q0Tot = energy1 - energy2;
+         Double_t qxTot = pX1 - pX2;
+         Double_t qyTot = pY1 - pY2;
+         Double_t qzTot = pZ1 - pZ2;
+         Double_t qInv = TMath::Sqrt(TMath::Abs(-TMath::Power(q0Tot,2) +TMath::Power(qxTot,2) +TMath::Power(qyTot,2) +TMath::Power(qzTot,2)));
+         ibin = Int_t(qInv/fP2Step);
+         //cout<<ibin<<endl;
+         if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
+         if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
+         if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
+         if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
+       }
+       else continue;
       }
     }
   }//case 5    
@@ -488,8 +661,15 @@ void AliBalance::CalculateBalance(TObjArray *gTrackArray) {
        track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
       else if(fAnalysisLevel == "MC")
        track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
-      Short_t charge1 = track1->Charge();
-      Double_t phi1 = track1->Phi();
+
+      Short_t charge1 = 0;
+      Double_t phi1 = 0.;
+      if(track1) {
+       charge1 = track1->Charge();
+       phi1 = track1->Phi();
+      }
+      else continue;
+
       for(j = 0; j < i; j++) {
        if(fAnalysisLevel == "ESD")
          track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
@@ -497,14 +677,18 @@ void AliBalance::CalculateBalance(TObjArray *gTrackArray) {
          track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
        else if(fAnalysisLevel == "MC")
          track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
-       Short_t charge2 = track2->Charge();
-       Double_t phi2 = track2->Phi();
-       Double_t dphi = TMath::Abs(phi1 - phi2);
-       ibin = Int_t(dphi/fP2Step);
-       if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
-       if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
-       if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
-       if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
+
+       if(track2) {
+         Short_t charge2 = track2->Charge();
+         Double_t phi2 = track2->Phi();
+         Double_t dphi = TMath::Abs(phi1 - phi2);
+         ibin = Int_t(dphi/fP2Step);
+         if((charge1 > 0)&&(charge2 > 0)) fNpp[ibin] += 1.;
+         if((charge1 < 0)&&(charge2 < 0)) fNnn[ibin] += 1.;
+         if((charge1 > 0)&&(charge2 < 0)) fNpn[ibin] += 1.;
+         if((charge1 < 0)&&(charge2 > 0)) fNpn[ibin] += 1.;
+       }
+       else continue;
       }
     }
   }//case 6
@@ -513,6 +697,33 @@ void AliBalance::CalculateBalance(TObjArray *gTrackArray) {
     Printf("bin: %d - Npp: %lf - Nnn: %lf - Nnp: %lf - Npn: %lf",i,fNpp[i],fNnn[i],fNpn[i],fNpn[i]);*/
 }
 
+//____________________________________________________________________//
+TH1F *AliBalance::GetHistNnn() {
+  //Return the control histogram of the -- component
+  for(Int_t iBin = 1; iBin <= fNumberOfBins; iBin++)
+    fHistfNnn->SetBinContent(iBin,GetNnn(iBin-1));
+
+  return fHistfNnn;
+}
+
+//____________________________________________________________________//
+TH1F *AliBalance::GetHistNpp() {
+  //Return the control histogram of the ++ component
+  for(Int_t iBin = 1; iBin <= fNumberOfBins; iBin++)
+    fHistfNpp->SetBinContent(iBin,GetNpp(iBin-1));
+
+  return fHistfNpp;
+}
+
+//____________________________________________________________________//
+TH1F *AliBalance::GetHistNpn() {
+  //Return the control histogram of the +- component
+  for(Int_t iBin = 1; iBin <= fNumberOfBins; iBin++)
+    fHistfNpn->SetBinContent(iBin,GetNpn(iBin-1));
+
+  return fHistfNpn;
+}
+
 //____________________________________________________________________//
 Double_t AliBalance::GetBalance(Int_t p2) {
   // Returns the value of the balance function in bin p2
@@ -606,20 +817,20 @@ TGraphErrors *AliBalance::DrawBalance() {
     gr->GetYaxis()->SetTitle("B(#Delta #eta)");
   }
   if(fAnalysisType==2) {
-    gr->GetXaxis()->SetTitle("Q_{long} [GeV]");
-    gr->GetYaxis()->SetTitle("B(Q_{long})");
+    gr->GetXaxis()->SetTitle("q_{long} (GeV/c)");
+    gr->GetYaxis()->SetTitle("B(q_{long}) [(GeV/c)^{-1}]");
   }
   if(fAnalysisType==3) {
-    gr->GetXaxis()->SetTitle("Q_{out} [GeV]");
-    gr->GetYaxis()->SetTitle("B(Q_{out})");
+    gr->GetXaxis()->SetTitle("q_{out} (GeV/c)");
+    gr->GetYaxis()->SetTitle("B(q_{out}) [(GeV/c)^{-1}]");
   }
   if(fAnalysisType==4) {
-    gr->GetXaxis()->SetTitle("Q_{side} [GeV]");
-    gr->GetYaxis()->SetTitle("B(Q_{side})");
+    gr->GetXaxis()->SetTitle("q_{side} (GeV/c)");
+    gr->GetYaxis()->SetTitle("B(q_{side}) [(GeV/c)^{-1}]");
   }
   if(fAnalysisType==5) {
-    gr->GetXaxis()->SetTitle("Q_{inv} [GeV]");
-    gr->GetYaxis()->SetTitle("B(Q_{inv})");
+    gr->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
+    gr->GetYaxis()->SetTitle("B(q_{inv}) [(GeV/c)^{-1}]");
   }
   if(fAnalysisType==6) {
     gr->GetXaxis()->SetTitle("#Delta #phi");