//ROOT
#include <Riostream.h>
+#include <TCanvas.h>
#include <TMath.h>
#include <TAxis.h>
#include <TH2D.h>
#include "AliMCParticle.h"
#include "AliESDtrack.h"
#include "AliAODTrack.h"
+#include "AliTHn.h"
+#include "AliAnalysisTaskTriggeredBF.h"
#include "AliBalancePsi.h"
//____________________________________________________________________//
AliBalancePsi::AliBalancePsi() :
TObject(),
- bShuffle(kFALSE),
+ fShuffle(kFALSE),
fAnalysisLevel("ESD"),
fAnalyzedEvents(0) ,
fCentralityId(0) ,
fCentStart(0.),
fCentStop(0.),
+ fHistP(0),
+ fHistN(0),
+ fHistPN(0),
+ fHistNP(0),
+ fHistPP(0),
+ fHistNN(0),
+ fHistHBTbefore(0),
+ fHistHBTafter(0),
+ fHistConversionbefore(0),
+ fHistConversionafter(0),
+ fHistPsiMinusPhi(0),
fPsiInterval(15.),
- fPsiNumberOfBins(24) {
+ fDeltaEtaMax(2.0),
+ fHBTCut(kFALSE),
+ fConversionCut(kFALSE),
+ fEventClass("EventPlane"){
// Default constructor
-
- for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
- if(i == 6) {
- fNumberOfBins[i] = 180;
- fP1Start[i] = -360.0;
- fP1Stop[i] = 360.0;
- fP2Start[i] = -360.0;
- fP2Stop[i] = 360.0;
- fP2Step[i] = 0.1;
- }
- else {
- fNumberOfBins[i] = 20;
- fP1Start[i] = -1.0;
- fP1Stop[i] = 1.0;
- fP2Start[i] = 0.0;
- fP2Stop[i] = 2.0;
- }
- fP2Step[i] = TMath::Abs(fP2Start - fP2Stop) / (Double_t)fNumberOfBins[i];
- fCentStart = 0.;
- fCentStop = 0.;
-
- for(Int_t k = 0; k < MAXIMUM_STEPS_IN_PSI; k++) {
- fNn[i][k] = 0.0;
- fNp[i][k] = 0.0;
-
- for(Int_t j = 0; j < MAXIMUM_NUMBER_OF_STEPS; j++) {
- fNpp[i][k][j] = .0;
- fNnn[i][k][j] = .0;
- fNpn[i][k][j] = .0;
- fNnp[i][k][j] = .0;
- fB[i][k][j] = 0.0;
- ferror[i][k][j] = 0.0;
- }
- }
- fHistP[i] = NULL;
- fHistN[i] = NULL;
- fHistPP[i] = NULL;
- fHistPN[i] = NULL;
- fHistNP[i] = NULL;
- fHistNN[i] = NULL;
- }
}
-
//____________________________________________________________________//
AliBalancePsi::AliBalancePsi(const AliBalancePsi& balance):
- TObject(balance), bShuffle(balance.bShuffle),
+ TObject(balance), fShuffle(balance.fShuffle),
fAnalysisLevel(balance.fAnalysisLevel),
fAnalyzedEvents(balance.fAnalyzedEvents),
fCentralityId(balance.fCentralityId),
fCentStart(balance.fCentStart),
fCentStop(balance.fCentStop),
+ fHistP(balance.fHistP),
+ fHistN(balance.fHistN),
+ fHistPN(balance.fHistPN),
+ fHistNP(balance.fHistNP),
+ fHistPP(balance.fHistPP),
+ fHistNN(balance.fHistNN),
+ fHistHBTbefore(balance.fHistHBTbefore),
+ fHistHBTafter(balance.fHistHBTafter),
+ fHistConversionbefore(balance.fHistConversionbefore),
+ fHistConversionafter(balance.fHistConversionafter),
+ fHistPsiMinusPhi(balance.fHistPsiMinusPhi),
fPsiInterval(balance.fPsiInterval),
- fPsiNumberOfBins(balance.fPsiNumberOfBins) {
+ fDeltaEtaMax(balance.fDeltaEtaMax),
+ fHBTCut(balance.fHBTCut),
+ fConversionCut(balance.fConversionCut),
+ fEventClass("EventPlane"){
//copy constructor
- for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
- fP1Start[i] = balance.fP1Start[i];
- fP1Stop[i] = balance.fP1Stop[i];
- fNumberOfBins[i] = balance.fNumberOfBins[i];
- fP2Start[i] = balance.fP2Start[i];
- fP2Stop[i] = balance.fP2Stop[i];
- fP2Step[i] = balance.fP2Step[i];
- fCentStart = balance.fCentStart;
- fCentStop = balance.fCentStop;
-
- fHistP[i] = balance.fHistP[i];
- fHistN[i] = balance.fHistN[i];
- fHistPN[i] = balance.fHistPN[i];
- fHistNP[i] = balance.fHistNP[i];
- fHistPP[i] = balance.fHistPP[i];
- fHistNN[i] = balance.fHistNN[i];
-
- for(Int_t k = 0; k < MAXIMUM_STEPS_IN_PSI; k++) {
- fNn[i][k] = balance.fNn[i][k];
- fNp[i][k] = balance.fNp[i][k];
-
- for(Int_t j = 0; j < MAXIMUM_NUMBER_OF_STEPS; j++) {
- fNpp[i][k][j] = balance.fNpp[i][k][j];
- fNnn[i][k][j] = balance.fNnn[i][k][j];
- fNpn[i][k][j] = balance.fNpn[i][k][j];
- fNnp[i][k][j] = balance.fNnp[i][k][j];
- fB[i][k][j] = balance.fB[i][k][j];
- ferror[i][k][j] = balance.ferror[i][k][j];
- }
- }
- }
}
//____________________________________________________________________//
AliBalancePsi::~AliBalancePsi() {
// Destructor
- for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
- delete fHistP[i];
- delete fHistN[i];
- delete fHistPN[i];
- delete fHistNP[i];
- delete fHistPP[i];
- delete fHistNN[i];
- }
-}
+ delete fHistP;
+ delete fHistN;
+ delete fHistPN;
+ delete fHistNP;
+ delete fHistPP;
+ delete fHistNN;
-//____________________________________________________________________//
-void AliBalancePsi::SetInterval(Int_t iAnalysisType,
- Double_t p1Start, Double_t p1Stop,
- Int_t ibins,
- Double_t p2Start, Double_t p2Stop,
- Double_t psiInterval) {
- // Sets the analyzed interval.
- // Set the same Information for all analyses
- fPsiInterval = psiInterval;
- fPsiNumberOfBins = (Int_t)(360./fPsiInterval);
-
- if(iAnalysisType == -1){
- for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
- fP1Start[i] = p1Start;
- fP1Stop[i] = p1Stop;
- fNumberOfBins[i] = ibins;
- fP2Start[i] = p2Start;
- fP2Stop[i] = p2Stop;
- fP2Step[i] = TMath::Abs(p2Start - p2Stop) / (Double_t)fNumberOfBins[i];
- }
- }
- // Set the Information for one analysis
- else if((iAnalysisType > -1) && (iAnalysisType < ANALYSIS_TYPES)) {
- fP1Start[iAnalysisType] = p1Start;
- fP1Stop[iAnalysisType] = p1Stop;
- fNumberOfBins[iAnalysisType] = ibins;
- fP2Start[iAnalysisType] = p2Start;
- fP2Stop[iAnalysisType] = p2Stop;
- fP2Step[iAnalysisType] = TMath::Abs(p2Start - p2Stop) / (Double_t)fNumberOfBins[iAnalysisType];
- }
- else {
- AliError("Wrong ANALYSIS number!");
- }
+ delete fHistHBTbefore;
+ delete fHistHBTafter;
+ delete fHistConversionbefore;
+ delete fHistConversionafter;
+ delete fHistPsiMinusPhi;
+
}
//____________________________________________________________________//
void AliBalancePsi::InitHistograms() {
- //Initialize the histograms
- TString histName;
- for(Int_t iAnalysisType = 0; iAnalysisType < ANALYSIS_TYPES; iAnalysisType++) {
- histName = "fHistP"; histName += gBFPsiAnalysisType[iAnalysisType];
- if(bShuffle) histName.Append("_shuffle");
- if(fCentralityId) histName += fCentralityId.Data();
- fHistP[iAnalysisType] = new TH3D(histName.Data(),"",(Int_t)(fCentStop-fCentStart),fCentStart,fCentStop,fPsiNumberOfBins,-fPsiInterval/2.,360.-fPsiInterval/2.,100,fP1Start[iAnalysisType],fP1Stop[iAnalysisType]);
-
- histName = "fHistN"; histName += gBFPsiAnalysisType[iAnalysisType];
- if(bShuffle) histName.Append("_shuffle");
- if(fCentralityId) histName += fCentralityId.Data();
- fHistN[iAnalysisType] = new TH3D(histName.Data(),"",(Int_t)(fCentStop-fCentStart),fCentStart,fCentStop,fPsiNumberOfBins,-fPsiInterval/2.,360.-fPsiInterval/2.,100,fP1Start[iAnalysisType],fP1Stop[iAnalysisType]);
+ // single particle histograms
+
+ // global switch disabling the reference
+ // (to avoid "Replacing existing TH1" if several wagons are created in train)
+ Bool_t oldStatus = TH1::AddDirectoryStatus();
+ TH1::AddDirectory(kFALSE);
+
+ Int_t anaSteps = 1; // analysis steps
+ Int_t iBinSingle[kTrackVariablesSingle]; // binning for track variables
+ Double_t* dBinsSingle[kTrackVariablesSingle]; // bins for track variables
+ TString axisTitleSingle[kTrackVariablesSingle]; // axis titles for track variables
- histName = "fHistPN"; histName += gBFPsiAnalysisType[iAnalysisType];
- if(bShuffle) histName.Append("_shuffle");
- if(fCentralityId) histName += fCentralityId.Data();
- fHistPN[iAnalysisType] = new TH3D(histName.Data(),"",(Int_t)(fCentStop-fCentStart),fCentStart,fCentStop,fPsiNumberOfBins,-fPsiInterval/2.,360.-fPsiInterval/2.,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
+ // two particle histograms
+ Int_t iBinPair[kTrackVariablesPair]; // binning for track variables
+ Double_t* dBinsPair[kTrackVariablesPair]; // bins for track variables
+ TString axisTitlePair[kTrackVariablesPair]; // axis titles for track variables
+ /**********************************************************
+
+ ======> Modification: Change Event Classification Scheme
+
+ ---> fEventClass == "EventPlane"
+
+ Default operation with Event Plane
+
+ ---> fEventClass == "Multiplicity"
+
+ Work with kTPCITStracklet multiplicity (from GetReferenceMultiplicity)
+
+ ---> fEventClass == "Centrality"
+
+ Work with Centrality Bins
+
+ ***********************************************************/
+
+ //--- Multiplicity Bins ------------------------------------
+ const Int_t kMultBins = 8;
+ //A first rough attempt at four bins
+ Double_t kMultBinLimits[kMultBins+1]={0,10,20,30,40,50,60,70,80};
+ //----------------------------------------------------------
- histName = "fHistNP"; histName += gBFPsiAnalysisType[iAnalysisType];
- if(bShuffle) histName.Append("_shuffle");
- if(fCentralityId) histName += fCentralityId.Data();
- fHistNP[iAnalysisType] = new TH3D(histName.Data(),"",(Int_t)(fCentStop-fCentStart),fCentStart,fCentStop,fPsiNumberOfBins,-fPsiInterval/2.,360.-fPsiInterval/2.,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
+ //--- Centrality Bins --------------------------------------
+ const Int_t kNCentralityBins = 9;
+ Double_t centralityBins[kNCentralityBins+1] = {0.,5.,10.,20.,30.,40.,50.,60.,70.,80.};
+ //----------------------------------------------------------
- histName = "fHistPP"; histName += gBFPsiAnalysisType[iAnalysisType];
- if(bShuffle) histName.Append("_shuffle");
- if(fCentralityId) histName += fCentralityId.Data();
- fHistPP[iAnalysisType] = new TH3D(histName.Data(),"",(Int_t)(fCentStop-fCentStart),fCentStart,fCentStop,fPsiNumberOfBins,-fPsiInterval/2.,360.-fPsiInterval/2.,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
+ //--- Event Plane Bins -------------------------------------
+ //Psi_2: -0.5->0.5 (in plane), 0.5->1.5 (intermediate), 1.5->2.5 (out of plane), 2.5->3.5 (rest)
+ const Int_t kNPsi2Bins = 4;
+ Double_t psi2Bins[kNPsi2Bins+1] = {-0.5,0.5,1.5,2.5,3.5};
+ //----------------------------------------------------------
- histName = "fHistNN"; histName += gBFPsiAnalysisType[iAnalysisType];
- if(bShuffle) histName.Append("_shuffle");
- if(fCentralityId) histName += fCentralityId.Data();
- fHistNN[iAnalysisType] = new TH3D(histName.Data(),"",(Int_t)(fCentStop-fCentStart),fCentStart,fCentStop,fPsiNumberOfBins,-fPsiInterval/2.,360.-fPsiInterval/2.,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
+ //Depending on fEventClass Variable, do one thing or the other...
+ if(fEventClass == "Multiplicity"){
+ iBinSingle[0] = kMultBins;
+ dBinsSingle[0] = kMultBinLimits;
+ axisTitleSingle[0] = "kTPCITStracklet multiplicity";
+ iBinPair[0] = kMultBins;
+ dBinsPair[0] = kMultBinLimits;
+ axisTitlePair[0] = "kTPCITStracklet multiplicity";
+ }
+ if(fEventClass == "Centrality"){
+ iBinSingle[0] = kNCentralityBins;
+ dBinsSingle[0] = centralityBins;
+ axisTitleSingle[0] = "Centrality percentile [%]";
+ iBinPair[0] = kNCentralityBins;
+ dBinsPair[0] = centralityBins;
+ axisTitlePair[0] = "Centrality percentile [%]";
+ }
+ if(fEventClass == "EventPlane"){
+ iBinSingle[0] = kNPsi2Bins;
+ dBinsSingle[0] = psi2Bins;
+ axisTitleSingle[0] = "#varphi - #Psi_{2} (a.u.)";
+ iBinPair[0] = kNPsi2Bins;
+ dBinsPair[0] = psi2Bins;
+ axisTitlePair[0] = "#varphi - #Psi_{2} (a.u.)";
+ }
+
+ // delta eta
+ const Int_t kNDeltaEtaBins = 80;
+ Double_t deltaEtaBins[kNDeltaEtaBins+1];
+ for(Int_t i = 0; i < kNDeltaEtaBins+1; i++)
+ deltaEtaBins[i] = - fDeltaEtaMax + i * 2 * fDeltaEtaMax / (Double_t)kNDeltaEtaBins;
+ iBinPair[1] = kNDeltaEtaBins;
+ dBinsPair[1] = deltaEtaBins;
+ axisTitlePair[1] = "#Delta#eta";
+
+ // delta phi
+ const Int_t kNDeltaPhiBins = 72;
+ Double_t deltaPhiBins[kNDeltaPhiBins+1];
+ for(Int_t i = 0; i < kNDeltaPhiBins+1; i++){
+ //deltaPhiBins[i] = -180.0 + i * 5.;
+ deltaPhiBins[i] = -TMath::Pi()/2. + i * 5.*TMath::Pi()/180.;
+ }
+ iBinPair[2] = kNDeltaPhiBins;
+ dBinsPair[2] = deltaPhiBins;
+ axisTitlePair[2] = "#Delta#varphi (rad)";
+
+ // pt(trigger-associated)
+ const Int_t kNPtBins = 16;
+ Double_t ptBins[kNPtBins+1] = {0.2,0.6,1.0,1.5,2.0,2.5,3.0,3.5,4.0,5.0,6.0,7.0,8.0,10.,12.,15.,20.};
+ //for(Int_t i = 0; i < kNPtBins+1; i++){
+ //ptBins[i] = 0.2 + i * 0.5;
+ //}
+ iBinSingle[1] = kNPtBins;
+ dBinsSingle[1] = ptBins;
+ axisTitleSingle[1] = "p_{T,trig.} (GeV/c)";
+
+ iBinPair[3] = kNPtBins;
+ dBinsPair[3] = ptBins;
+ axisTitlePair[3] = "p_{T,trig.} (GeV/c)";
+
+ iBinPair[4] = kNPtBins;
+ dBinsPair[4] = ptBins;
+ axisTitlePair[4] = "p_{T,assoc.} (GeV/c)";
+
+ TString histName;
+ //+ triggered particles
+ histName = "fHistP";
+ if(fShuffle) histName.Append("_shuffle");
+ if(fCentralityId) histName += fCentralityId.Data();
+ fHistP = new AliTHn(histName.Data(),histName.Data(),anaSteps,kTrackVariablesSingle,iBinSingle);
+ for (Int_t j=0; j<kTrackVariablesSingle; j++) {
+ fHistP->SetBinLimits(j, dBinsSingle[j]);
+ fHistP->SetVarTitle(j, axisTitleSingle[j]);
}
-}
-//____________________________________________________________________//
-void AliBalancePsi::PrintAnalysisSettings() {
+ //- triggered particles
+ histName = "fHistN";
+ if(fShuffle) histName.Append("_shuffle");
+ if(fCentralityId) histName += fCentralityId.Data();
+ fHistN = new AliTHn(histName.Data(),histName.Data(),anaSteps,kTrackVariablesSingle,iBinSingle);
+ for (Int_t j=0; j<kTrackVariablesSingle; j++) {
+ fHistN->SetBinLimits(j, dBinsSingle[j]);
+ fHistN->SetVarTitle(j, axisTitleSingle[j]);
+ }
- Printf("======================================");
- Printf("Analysis level: %s",fAnalysisLevel.Data());
- Printf("======================================");
- for(Int_t ibin = 0; ibin < ANALYSIS_TYPES; ibin++){
- Printf("Interval info for variable %d",ibin);
- Printf("Analyzed interval (min.): %lf",fP2Start[ibin]);
- Printf("Analyzed interval (max.): %lf",fP2Stop[ibin]);
- Printf("Number of bins: %d",fNumberOfBins[ibin]);
- Printf("Step: %lf",fP2Step[ibin]);
- Printf(" ");
- }
- Printf("======================================");
+ //+- pairs
+ histName = "fHistPN";
+ if(fShuffle) histName.Append("_shuffle");
+ if(fCentralityId) histName += fCentralityId.Data();
+ fHistPN = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
+ for (Int_t j=0; j<kTrackVariablesPair; j++) {
+ fHistPN->SetBinLimits(j, dBinsPair[j]);
+ fHistPN->SetVarTitle(j, axisTitlePair[j]);
+ }
+
+ //-+ pairs
+ histName = "fHistNP";
+ if(fShuffle) histName.Append("_shuffle");
+ if(fCentralityId) histName += fCentralityId.Data();
+ fHistNP = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
+ for (Int_t j=0; j<kTrackVariablesPair; j++) {
+ fHistNP->SetBinLimits(j, dBinsPair[j]);
+ fHistNP->SetVarTitle(j, axisTitlePair[j]);
+ }
+
+ //++ pairs
+ histName = "fHistPP";
+ if(fShuffle) histName.Append("_shuffle");
+ if(fCentralityId) histName += fCentralityId.Data();
+ fHistPP = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
+ for (Int_t j=0; j<kTrackVariablesPair; j++) {
+ fHistPP->SetBinLimits(j, dBinsPair[j]);
+ fHistPP->SetVarTitle(j, axisTitlePair[j]);
+ }
+
+ //-- pairs
+ histName = "fHistNN";
+ if(fShuffle) histName.Append("_shuffle");
+ if(fCentralityId) histName += fCentralityId.Data();
+ fHistNN = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
+ for (Int_t j=0; j<kTrackVariablesPair; j++) {
+ fHistNN->SetBinLimits(j, dBinsPair[j]);
+ fHistNN->SetVarTitle(j, axisTitlePair[j]);
+ }
+ AliInfo("Finished setting up the AliTHn");
+
+ // QA histograms
+ fHistHBTbefore = new TH2D("fHistHBTbefore","before HBT cut",200,0,2,200,0,2.*TMath::Pi());
+ fHistHBTafter = new TH2D("fHistHBTafter","after HBT cut",200,0,2,200,0,2.*TMath::Pi());
+ fHistConversionbefore = new TH2D("fHistConversionbefore","before Conversion cut",200,0,2,200,0,2.*TMath::Pi());
+ fHistConversionafter = new TH2D("fHistConversionafter","after Conversion cut",200,0,2,200,0,2.*TMath::Pi());
+ fHistPsiMinusPhi = new TH2D("fHistPsiMinusPhi","",4,-0.5,3.5,100,0,2.*TMath::Pi());
+
+ TH1::AddDirectory(oldStatus);
+
}
//____________________________________________________________________//
-void AliBalancePsi::CalculateBalance(Float_t fCentrality,
- Double_t gReactionPlane,
- vector<Double_t> **chargeVector) {
+void AliBalancePsi::CalculateBalance(Double_t gReactionPlane,
+ TObjArray *particles,
+ TObjArray *particlesMixed,
+ Float_t bSign,
+ Double_t kMultorCent) {
// Calculates the balance function
fAnalyzedEvents++;
- Int_t i = 0 , j = 0;
- Int_t iBin = 0;
-
+
// Initialize histograms if not done yet
- if(!fHistPN[0]){
+ if(!fHistPN){
AliWarning("Histograms not yet initialized! --> Will be done now");
AliWarning("This works only in local mode --> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
InitHistograms();
}
- Int_t gNtrack = chargeVector[0]->size();
- //Printf("(AliBalancePsi) Number of tracks: %d",gNtrack);
- Double_t gPsiMinusPhi = 0., gPsiMinusPhiPrime = 0.;
- Int_t gPsiBin = -10;
- for(i = 0; i < gNtrack;i++){
- Double_t charge = chargeVector[0]->at(i);
- Double_t rapidity = chargeVector[1]->at(i);
- Double_t pseudorapidity = chargeVector[2]->at(i);
- Double_t phi = chargeVector[3]->at(i);
- gPsiMinusPhi = gReactionPlane - phi;
- gPsiMinusPhiPrime = gPsiMinusPhi;
- if(gPsiMinusPhi < -fPsiInterval/2) gPsiMinusPhiPrime = 360. - gPsiMinusPhiPrime;
- gPsiBin = (Int_t)((gPsiMinusPhiPrime + 360./fPsiNumberOfBins/2.)/(360./fPsiNumberOfBins));
- //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
- for(Int_t iAnalysisType = 0; iAnalysisType < ANALYSIS_TYPES; iAnalysisType++) {
- if(iAnalysisType == kEta) {
- if((pseudorapidity >= fP1Start[iAnalysisType]) && (pseudorapidity <= fP1Stop[iAnalysisType])) {
- if(charge > 0) {
- fNp[iAnalysisType][gPsiBin] += 1.;
- fHistP[iAnalysisType]->Fill(fCentrality,gPsiMinusPhi,pseudorapidity);
- }//charge > 0
- if(charge < 0) {
- fNn[iAnalysisType][gPsiBin] += 1.;
- fHistN[iAnalysisType]->Fill(fCentrality,gPsiMinusPhi,pseudorapidity);
- }//charge < 0
- }//p1 interval check
- }//analysis type: eta
- else if(iAnalysisType == kPhi) {
- if((phi >= fP1Start[iAnalysisType]) && (phi <= fP1Stop[iAnalysisType])) {
- if(charge > 0) {
- fNp[iAnalysisType][gPsiBin] += 1.;
- fHistP[iAnalysisType]->Fill(fCentrality,gPsiMinusPhi,phi);
- }//charge > 0
- if(charge < 0) {
- fNn[iAnalysisType][gPsiBin] += 1.;
- fHistN[iAnalysisType]->Fill(fCentrality,gPsiMinusPhi,phi);
- }//charge < 0
- }//p1 interval check
- }//analysis type: phi
- else {
- if((rapidity >= fP1Start[iAnalysisType]) && (rapidity <= fP1Stop[iAnalysisType])) {
- if(charge > 0) {
- fNp[iAnalysisType][gPsiBin] += 1.;
- fHistP[iAnalysisType]->Fill(fCentrality,gPsiMinusPhi,rapidity);
- }//charge > 0
- if(charge < 0) {
- fNn[iAnalysisType][gPsiBin] += 1.;
- fHistN[iAnalysisType]->Fill(fCentrality,gPsiMinusPhi,rapidity);
- }//charge < 0
- }//p1 interval check
- }//analysis type: y, qside, qout, qlong, qinv
- }//analysis type loop
+ Double_t trackVariablesSingle[kTrackVariablesSingle];
+ Double_t trackVariablesPair[kTrackVariablesPair];
+
+ if (!particles){
+ AliWarning("particles TObjArray is NULL pointer --> return");
+ return;
}
- //Printf("Np: %lf - Nn: %lf",fNp[0],fNn[0]);
+ // define end of particle loops
+ Int_t iMax = particles->GetEntriesFast();
+ Int_t jMax = iMax;
+ if (particlesMixed)
+ jMax = particlesMixed->GetEntriesFast();
+
+ // Eta() is extremely time consuming, therefore cache it for the inner loop here:
+ TObjArray* particlesSecond = (particlesMixed) ? particlesMixed : particles;
+
+ TArrayF secondEta(jMax);
+ TArrayF secondPhi(jMax);
+ TArrayF secondPt(jMax);
+ TArrayS secondCharge(jMax);
+ TArrayD secondCorrection(jMax);
+
+ for (Int_t i=0; i<jMax; i++){
+ secondEta[i] = ((AliVParticle*) particlesSecond->At(i))->Eta();
+ secondPhi[i] = ((AliVParticle*) particlesSecond->At(i))->Phi();
+ secondPt[i] = ((AliVParticle*) particlesSecond->At(i))->Pt();
+ secondCharge[i] = (Short_t)((AliVParticle*) particlesSecond->At(i))->Charge();
+ secondCorrection[i] = (Double_t)((AliBFBasicParticle*) particlesSecond->At(i))->Correction(); //==========================correction
+ }
- Double_t dy = 0., deta = 0.;
- Double_t qLong = 0., qOut = 0., qSide = 0., qInv = 0.;
- Double_t dphi = 0.;
-
- Double_t charge1 = 0;
- Double_t eta1 = 0., rap1 = 0.;
- Double_t px1 = 0., py1 = 0., pz1 = 0.;
- Double_t pt1 = 0.;
- Double_t energy1 = 0.;
- Double_t phi1 = 0.;
-
- Double_t charge2 = 0;
- Double_t eta2 = 0., rap2 = 0.;
- Double_t px2 = 0., py2 = 0., pz2 = 0.;
- Double_t pt2 = 0.;
- Double_t energy2 = 0.;
- Double_t phi2 = 0.;
- //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
- for(i = 1; i < gNtrack; i++) {
-
- charge1 = chargeVector[0]->at(i);
- rap1 = chargeVector[1]->at(i);
- eta1 = chargeVector[2]->at(i);
- phi1 = chargeVector[3]->at(i);
- px1 = chargeVector[4]->at(i);
- py1 = chargeVector[5]->at(i);
- pz1 = chargeVector[6]->at(i);
- pt1 = chargeVector[7]->at(i);
- energy1 = chargeVector[8]->at(i);
- gPsiMinusPhi = gReactionPlane - phi1;
- phi1 -= 360.;
- gPsiMinusPhiPrime = gPsiMinusPhi;
- if(gPsiMinusPhi < -fPsiInterval/2) gPsiMinusPhiPrime = 360. - gPsiMinusPhiPrime;
- gPsiBin = (Int_t)((gPsiMinusPhiPrime + 360./fPsiNumberOfBins/2.)/(360./fPsiNumberOfBins));
-
- for(j = 0; j < i; j++) {
- charge2 = chargeVector[0]->at(j);
- rap2 = chargeVector[1]->at(j);
- eta2 = chargeVector[2]->at(j);
- phi2 = chargeVector[3]->at(j);
- px2 = chargeVector[4]->at(j);
- py2 = chargeVector[5]->at(j);
- pz2 = chargeVector[6]->at(j);
- pt2 = chargeVector[7]->at(i);
- energy2 = chargeVector[8]->at(j);
- phi2 -= 360.;
-
- // filling the arrays
- // RAPIDITY
- dy = rap1 - rap2;
-
- // Eta
- deta = eta1 - eta2;
-
- //qlong
- 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 qxTot = px1 - px2;
- Double_t qyTot = py1 - py2;
- Double_t qzTot = pz1 - pz2;
-
- Double_t eTot2 = eTot*eTot;
- Double_t pTot2 = pxTot*pxTot + pyTot*pyTot + pzTot*pzTot;
- Double_t pzTot2 = pzTot*pzTot;
+ // 1st particle loop
+ for (Int_t i = 0; i < iMax; i++) {
+ //AliVParticle* firstParticle = (AliVParticle*) particles->At(i);
+ AliBFBasicParticle* firstParticle = (AliBFBasicParticle*) particles->At(i); //==========================correction
+
+ // some optimization
+ Float_t firstEta = firstParticle->Eta();
+ Float_t firstPhi = firstParticle->Phi();
+ Float_t firstPt = firstParticle->Pt();
+ Float_t firstCorrection = firstParticle->Correction();//==========================correction
- Double_t q0Tot2 = q0Tot*q0Tot;
- Double_t qTot2 = qxTot*qxTot + qyTot*qyTot + qzTot*qzTot;
+ // Event plane (determine psi bin)
+ Double_t gPsiMinusPhi = 0.;
+ Double_t gPsiMinusPhiBin = -10.;
+ gPsiMinusPhi = TMath::Abs(firstPhi - gReactionPlane);
+ //in-plane
+ if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
+ ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
+ gPsiMinusPhiBin = 0.0;
+ //intermediate
+ else if(((37.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5*TMath::DegToRad()))||
+ ((127.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5*TMath::DegToRad()))||
+ ((217.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5*TMath::DegToRad()))||
+ ((307.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5*TMath::DegToRad())))
+ gPsiMinusPhiBin = 1.0;
+ //out of plane
+ else if(((82.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5*TMath::DegToRad()))||
+ ((262.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5*TMath::DegToRad())))
+ gPsiMinusPhiBin = 2.0;
+ //everything else
+ else
+ gPsiMinusPhiBin = 3.0;
+
+ fHistPsiMinusPhi->Fill(gPsiMinusPhiBin,gPsiMinusPhi);
- Double_t snn = eTot2 - pTot2;
- Double_t ptTot2 = pTot2 - pzTot2 ;
- Double_t ptTot = TMath::Sqrt( ptTot2 );
-
- qLong = TMath::Abs(eTot*qzTot - pzTot*q0Tot)/TMath::Sqrt(snn + ptTot2);
-
- //qout
- qOut = TMath::Sqrt(snn/(snn + ptTot2)) * TMath::Abs(pxTot*qxTot + pyTot*qyTot)/ptTot;
-
- //qside
- qSide = TMath::Abs(pxTot*qyTot - pyTot*qxTot)/ptTot;
-
- //qinv
- qInv = TMath::Sqrt(TMath::Abs(-q0Tot2 + qTot2 ));
-
- //phi
- dphi = phi1 - phi2;
- if(dphi < -180) dphi = 360 - dphi; //dphi should be between -180 and 180!
+ Short_t charge1 = (Short_t) firstParticle->Charge();
+
+ trackVariablesSingle[0] = gPsiMinusPhiBin;
+ trackVariablesSingle[1] = firstPt;
+ if(fEventClass=="Multiplicity" || fEventClass == "Centrality" ) trackVariablesSingle[0] = kMultorCent;
+
+ //fill single particle histograms
+ if(charge1 > 0) fHistP->Fill(trackVariablesSingle,0,firstCorrection); //==========================correction
+ else if(charge1 < 0) fHistN->Fill(trackVariablesSingle,0,firstCorrection); //==========================correction
+
+ // 2nd particle loop
+ for(Int_t j = 0; j < jMax; j++) {
- //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
- if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
+ if(!particlesMixed && j == i) continue; // no auto correlations (only for non mixing)
- // rapidity
- if( dy > fP2Start[kRapidity] && dy < fP2Stop[kRapidity]){
- iBin = Int_t((dy-fP2Start[kRapidity])/fP2Step[kRapidity]);
- if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
-
- if((charge1 > 0)&&(charge2 > 0)) {
- fNpp[kRapidity][gPsiBin][iBin] += 1.;
- fHistPP[kRapidity]->Fill(fCentrality,gPsiMinusPhi,dy);
- }
- else if((charge1 < 0)&&(charge2 < 0)) {
- fNnn[kRapidity][gPsiBin][iBin] += 1.;
- fHistNN[kRapidity]->Fill(fCentrality,gPsiMinusPhi,dy);
- }
- else if((charge1 > 0)&&(charge2 < 0)) {
- fNpn[kRapidity][gPsiBin][iBin] += 1.;
- fHistPN[kRapidity]->Fill(fCentrality,gPsiMinusPhi,dy);
- }
- else if((charge1 < 0)&&(charge2 > 0)) {
- fNpn[kRapidity][gPsiBin][iBin] += 1.;
- fHistPN[kRapidity]->Fill(fCentrality,gPsiMinusPhi,dy);
- }
- }//BF binning check
- }//p2 interval check
- }//p1 interval check
-
- // pseudorapidity
- if((eta1 >= fP1Start[kEta]) && (eta1 <= fP1Stop[kEta]) && (eta2 >= fP1Start[kEta]) && (eta2 <= fP1Stop[kEta])) {
- if( deta > fP2Start[kEta] && deta < fP2Stop[kEta]){
- iBin = Int_t((deta-fP2Start[kEta])/fP2Step[kEta]);
- if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
- if((charge1 > 0)&&(charge2 > 0)) {
- fNpp[kEta][gPsiBin][iBin] += 1.;
- fHistPP[kEta]->Fill(fCentrality,gPsiMinusPhi,deta);
- }
- if((charge1 < 0)&&(charge2 < 0)) {
- fNnn[kEta][gPsiBin][iBin] += 1.;
- fHistNN[kEta]->Fill(fCentrality,gPsiMinusPhi,deta);
- }
- if((charge1 > 0)&&(charge2 < 0)) {
- fNpn[kEta][gPsiBin][iBin] += 1.;
- fHistPN[kEta]->Fill(fCentrality,gPsiMinusPhi,deta);
- }
- if((charge1 < 0)&&(charge2 > 0)) {
- fNpn[kEta][gPsiBin][iBin] += 1.;
- fHistPN[kEta]->Fill(fCentrality,gPsiMinusPhi,deta);
- }
- }//BF binning check
- }//p2 interval check
- }//p1 interval check
-
- // Qlong, out, side, inv
- // Check the p1 intervall for rapidity here (like for single tracks above)
- if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
- if( qLong > fP2Start[kQlong] && qLong < fP2Stop[kQlong]){
- iBin = Int_t((qLong-fP2Start[kQlong])/fP2Step[kQlong]);
- if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
- if((charge1 > 0)&&(charge2 > 0)) {
- fNpp[kQlong][gPsiBin][iBin] += 1.;
- fHistPP[kQlong]->Fill(fCentrality,gPsiMinusPhi,qLong);
- }
- if((charge1 < 0)&&(charge2 < 0)) {
- fNnn[kQlong][gPsiBin][iBin] += 1.;
- fHistNN[kQlong]->Fill(fCentrality,gPsiMinusPhi,qLong);
- }
- if((charge1 > 0)&&(charge2 < 0)) {
- fNpn[kQlong][gPsiBin][iBin] += 1.;
- fHistPN[kQlong]->Fill(fCentrality,gPsiMinusPhi,qLong);
- }
- if((charge1 < 0)&&(charge2 > 0)) {
- fNpn[kQlong][gPsiBin][iBin] += 1.;
- fHistPN[kQlong]->Fill(fCentrality,gPsiMinusPhi,qLong);
- }
- }//BF binning check
- }//p2 interval check
- }//p1 interval check
-
- if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
- if( qOut > fP2Start[kQout] && qOut < fP2Stop[kQout]){
- iBin = Int_t((qOut-fP2Start[kQout])/fP2Step[kQout]);
- if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
- if((charge1 > 0)&&(charge2 > 0)) {
- fNpp[kQout][gPsiBin][iBin] += 1.;
- fHistPP[kQout]->Fill(fCentrality,gPsiMinusPhi,qOut);
- }
- if((charge1 < 0)&&(charge2 < 0)) {
- fNnn[kQout][gPsiBin][iBin] += 1.;
- fHistNN[kQout]->Fill(fCentrality,gPsiMinusPhi,qOut);
- }
- if((charge1 > 0)&&(charge2 < 0)) {
- fNpn[kQout][gPsiBin][iBin] += 1.;
- fHistPN[kQout]->Fill(fCentrality,gPsiMinusPhi,qOut);
- }
- if((charge1 < 0)&&(charge2 > 0)) {
- fNpn[kQout][gPsiBin][iBin] += 1.;
- fHistPN[kQout]->Fill(fCentrality,gPsiMinusPhi,qOut);
- }
- }//BF binning check
- }//p2 interval check
- }//p1 interval check
+ // pT,Assoc < pT,Trig
+ if(firstPt < secondPt[j]) continue;
+
+ Short_t charge2 = secondCharge[j];
+
+ trackVariablesPair[0] = trackVariablesSingle[0];
+ trackVariablesPair[1] = firstEta - secondEta[j]; // delta eta
+ trackVariablesPair[2] = firstPhi - secondPhi[j]; // delta phi
+ //if (trackVariablesPair[2] > 180.) // delta phi between -180 and 180
+ //trackVariablesPair[2] -= 360.;
+ //if (trackVariablesPair[2] < - 180.)
+ //trackVariablesPair[2] += 360.;
+ if (trackVariablesPair[2] > TMath::Pi()) // delta phi between -pi and pi
+ trackVariablesPair[2] -= 2.*TMath::Pi();
+ if (trackVariablesPair[2] < - TMath::Pi())
+ trackVariablesPair[2] += 2.*TMath::Pi();
+ if (trackVariablesPair[2] < - TMath::Pi()/2.)
+ trackVariablesPair[2] += 2.*TMath::Pi();
+
+ trackVariablesPair[3] = firstPt; // pt trigger
+ trackVariablesPair[4] = secondPt[j]; // pt
+ // trackVariablesPair[5] = fCentrality; // centrality
+
+ // HBT like cut
+ if(fHBTCut){ // VERSION 3 (all pairs)
+ //if(fHBTCut && charge1 * charge2 > 0){ // VERSION 2 (only for LS)
+ //if( dphi < 3 || deta < 0.01 ){ // VERSION 1
+ // continue;
- if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
- if( qSide > fP2Start[kQside] && qSide < fP2Stop[kQside]){
- iBin = Int_t((qSide-fP2Start[kQside])/fP2Step[kQside]);
- if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
- if((charge1 > 0)&&(charge2 > 0)) {
- fNpp[kQside][gPsiBin][iBin] += 1.;
- fHistPP[kQside]->Fill(fCentrality,gPsiMinusPhi,qSide);
- }
- if((charge1 < 0)&&(charge2 < 0)) {
- fNnn[kQside][gPsiBin][iBin] += 1.;
- fHistNN[kQside]->Fill(fCentrality,gPsiMinusPhi,qSide);
- }
- if((charge1 > 0)&&(charge2 < 0)) {
- fNpn[kQside][gPsiBin][iBin] += 1.;
- fHistPN[kQside]->Fill(fCentrality,gPsiMinusPhi,qSide);
- }
- if((charge1 < 0)&&(charge2 > 0)) {
- fNpn[kQside][gPsiBin][iBin] += 1.;
- fHistPN[kQside]->Fill(fCentrality,gPsiMinusPhi,qSide);
- }
- }//BF binning check
- }//p2 interval check
- }//p1 interval check
+ Double_t deta = firstEta - secondEta[j];
+ Double_t dphi = firstPhi - secondPhi[j];
+ // VERSION 2 (Taken from DPhiCorrelations)
+ // the variables & cuthave been developed by the HBT group
+ // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
+ fHistHBTbefore->Fill(deta,dphi);
- if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
- if( qInv > fP2Start[kQinv] && qInv < fP2Stop[kQinv]){
- iBin = Int_t((qInv-fP2Start[kQinv])/fP2Step[kQinv]);
- if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
- if((charge1 > 0)&&(charge2 > 0)) {
- fNpp[kQinv][gPsiBin][iBin] += 1.;
- fHistPP[kQinv]->Fill(fCentrality,gPsiMinusPhi,qInv);
+ // optimization
+ if (TMath::Abs(deta) < 0.02 * 2.5 * 3) //twoTrackEfficiencyCutValue = 0.02 [default for dphicorrelations]
+ {
+ // phi in rad
+ //Float_t phi1rad = firstPhi*TMath::DegToRad();
+ //Float_t phi2rad = secondPhi[j]*TMath::DegToRad();
+ Float_t phi1rad = firstPhi;
+ Float_t phi2rad = secondPhi[j];
+
+ // check first boundaries to see if is worth to loop and find the minimum
+ Float_t dphistar1 = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, 0.8, bSign);
+ Float_t dphistar2 = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, 2.5, bSign);
+
+ const Float_t kLimit = 0.02 * 3;
+
+ Float_t dphistarminabs = 1e5;
+ Float_t dphistarmin = 1e5;
+
+ if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0 ) {
+ for (Double_t rad=0.8; rad<2.51; rad+=0.01) {
+ Float_t dphistar = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, rad, bSign);
+ Float_t dphistarabs = TMath::Abs(dphistar);
+
+ if (dphistarabs < dphistarminabs) {
+ dphistarmin = dphistar;
+ dphistarminabs = dphistarabs;
+ }
}
- if((charge1 < 0)&&(charge2 < 0)) {
- fNnn[kQinv][gPsiBin][iBin] += 1.;
- fHistNN[kQinv]->Fill(fCentrality,gPsiMinusPhi,qInv);
- }
- if((charge1 > 0)&&(charge2 < 0)) {
- fNpn[kQinv][gPsiBin][iBin] += 1.;
- fHistPN[kQinv]->Fill(fCentrality,gPsiMinusPhi,qInv);
- }
- if((charge1 < 0)&&(charge2 > 0)) {
- fNpn[kQinv][gPsiBin][iBin] += 1.;
- fHistPN[kQinv]->Fill(fCentrality,gPsiMinusPhi,qInv);
+
+ if (dphistarminabs < 0.02 && TMath::Abs(deta) < 0.02) {
+ //AliInfo(Form("HBT: Removed track pair %d %d with [[%f %f]] %f %f %f | %f %f %d %f %f %d %f", i, j, deta, dphi, dphistarminabs, dphistar1, dphistar2, phi1rad, pt1, charge1, phi2rad, pt2, charge2, bSign));
+ continue;
}
- }//BF binning check
- }//p2 interval check
- }//p1 interval check
+ }
+ }
+ fHistHBTafter->Fill(deta,dphi);
+ }//HBT cut
- // Phi
- if((phi1 >= fP1Start[kPhi]) && (phi1 <= fP1Stop[kPhi]) && (phi2 >= fP1Start[kPhi]) && (phi2 <= fP1Stop[kPhi])) {
- if( dphi > fP2Start[kPhi] && dphi < fP2Stop[kPhi]){
- iBin = Int_t((dphi-fP2Start[kPhi])/fP2Step[kPhi]);
- if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
- if((charge1 > 0)&&(charge2 > 0)) {
- fNpp[kPhi][gPsiBin][iBin] += 1.;
- fHistPP[kPhi]->Fill(fCentrality,gPsiMinusPhi,dphi);
- }
- if((charge1 < 0)&&(charge2 < 0)) {
- fNnn[kPhi][gPsiBin][iBin] += 1.;
- fHistNN[kPhi]->Fill(fCentrality,gPsiMinusPhi,dphi);
- }
- if((charge1 > 0)&&(charge2 < 0)) {
- fNpn[kPhi][gPsiBin][iBin] += 1.;
- fHistPN[kPhi]->Fill(fCentrality,gPsiMinusPhi,dphi);
- }
- if((charge1 < 0)&&(charge2 > 0)) {
- fNpn[kPhi][gPsiBin][iBin] += 1.;
- fHistPN[kPhi]->Fill(fCentrality,gPsiMinusPhi,dphi);
- }
- }//BF binning check
- }//p2 interval check
- }//p1 interval check
+ // conversions
+ if(fConversionCut) {
+ if (charge1 * charge2 < 0) {
+ Double_t deta = firstEta - secondEta[j];
+ Double_t dphi = firstPhi - secondPhi[j];
+ fHistConversionbefore->Fill(deta,dphi);
+
+ Float_t m0 = 0.510e-3;
+ Float_t tantheta1 = 1e10;
+
+ // phi in rad
+ //Float_t phi1rad = firstPhi*TMath::DegToRad();
+ //Float_t phi2rad = secondPhi[j]*TMath::DegToRad();
+ Float_t phi1rad = firstPhi;
+ Float_t phi2rad = secondPhi[j];
+
+ if (firstEta < -1e-10 || firstEta > 1e-10)
+ tantheta1 = 2 * TMath::Exp(-firstEta) / ( 1 - TMath::Exp(-2*firstEta));
+
+ Float_t tantheta2 = 1e10;
+ if (secondEta[j] < -1e-10 || secondEta[j] > 1e-10)
+ tantheta2 = 2 * TMath::Exp(-secondEta[j]) / ( 1 - TMath::Exp(-2*secondEta[j]));
+
+ Float_t e1squ = m0 * m0 + firstPt * firstPt * (1.0 + 1.0 / tantheta1 / tantheta1);
+ Float_t e2squ = m0 * m0 + secondPt[j] * secondPt[j] * (1.0 + 1.0 / tantheta2 / tantheta2);
+
+ Float_t masssqu = 2 * m0 * m0 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( firstPt * secondPt[j] * ( TMath::Cos(phi1rad - phi2rad) + 1.0 / tantheta1 / tantheta2 ) ) );
+
+ if (masssqu < 0.04*0.04){
+ //AliInfo(Form("Conversion: Removed track pair %d %d with [[%f %f] %f %f] %d %d <- %f %f %f %f %f %f ", i, j, deta, dphi, masssqu, charge1, charge2,eta1,eta2,phi1,phi2,pt1,pt2));
+ continue;
+ }
+ fHistConversionafter->Fill(deta,dphi);
+ }
+ }//conversion cut
+
+ if( charge1 > 0 && charge2 < 0) fHistPN->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]); //==========================correction
+ else if( charge1 < 0 && charge2 > 0) fHistNP->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction
+ else if( charge1 > 0 && charge2 > 0) fHistPP->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction
+ else if( charge1 < 0 && charge2 < 0) fHistNN->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction
+ else {
+ //AliWarning(Form("Wrong charge combination: charge1 = %d and charge2 = %d",charge,charge2));
+ continue;
+ }
}//end of 2nd particle loop
}//end of 1st particle loop
- //Printf("Number of analyzed events: %i",fAnalyzedEvents);
- //Printf("DeltaEta NN[0] = %.0f, PP[0] = %.0f, NP[0] = %.0f, PN[0] = %.0f",fNnn[kEta][0],fNpp[kEta][0],fNnp[kEta][0],fNpn[kEta][0]);
}
-
//____________________________________________________________________//
-Double_t AliBalancePsi::GetBalance(Int_t iAnalysisType, Int_t gPsiBin,
- Int_t p2) {
- // Returns the value of the balance function in bin p2
- fB[iAnalysisType][gPsiBin][p2] = 0.5*(((fNpn[iAnalysisType][gPsiBin][p2] - 2.*fNnn[iAnalysisType][gPsiBin][p2])/fNn[iAnalysisType][gPsiBin]) + ((fNpn[iAnalysisType][gPsiBin][p2] - 2.*fNpp[iAnalysisType][gPsiBin][p2])/fNp[iAnalysisType][gPsiBin]))/fP2Step[iAnalysisType];
-
- return fB[iAnalysisType][gPsiBin][p2];
-}
+TH1D *AliBalancePsi::GetBalanceFunctionHistogram(Int_t iVariableSingle,
+ Int_t iVariablePair,
+ Double_t psiMin,
+ Double_t psiMax,
+ Double_t ptTriggerMin,
+ Double_t ptTriggerMax,
+ Double_t ptAssociatedMin,
+ Double_t ptAssociatedMax) {
+ //Returns the BF histogram, extracted from the 6 AliTHn objects
+ //(private members) of the AliBalancePsi class.
+ //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
+ //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
+ // Psi_2
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+
+ // pt trigger
+ if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ }
+
+ // pt associated
+ if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
+ fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ }
+
+ //Printf("P:%lf - N:%lf - PN:%lf - NP:%lf - PP:%lf - NN:%lf",fHistP->GetEntries(0),fHistN->GetEntries(0),fHistPN->GetEntries(0),fHistNP->GetEntries(0),fHistPP->GetEntries(0),fHistNN->GetEntries(0));
+
+ // Project into the wanted space (1st: analysis step, 2nd: axis)
+ TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair);
+ TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair);
+ TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair);
+ TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair);
+ TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle);
+ TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle);
+
+ TH1D *gHistBalanceFunctionHistogram = 0x0;
+ if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
+ gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone();
+ gHistBalanceFunctionHistogram->Reset();
-//____________________________________________________________________//
-Double_t AliBalancePsi::GetError(Int_t iAnalysisType, Int_t gPsiBin,
- Int_t p2) {
- // Returns the error on the BF value for bin p2
- // The errors for fNn and fNp are neglected here (0.1 % of total error)
- ferror[iAnalysisType][gPsiBin][p2] = TMath::Sqrt( Double_t(fNpp[iAnalysisType][gPsiBin][p2])/(Double_t(fNp[iAnalysisType][gPsiBin])*Double_t(fNp[iAnalysisType][gPsiBin])) +
- Double_t(fNnn[iAnalysisType][gPsiBin][p2])/(Double_t(fNn[iAnalysisType][gPsiBin])*Double_t(fNn[iAnalysisType][gPsiBin])) +
- Double_t(fNpn[iAnalysisType][gPsiBin][p2])*TMath::Power((0.5/Double_t(fNp[iAnalysisType][gPsiBin]) + 0.5/Double_t(fNn[iAnalysisType][gPsiBin])),2))/fP2Step[iAnalysisType];
-
- return ferror[iAnalysisType][gPsiBin][p2];
+ switch(iVariablePair) {
+ case 1:
+ gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
+ gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
+ break;
+ case 2:
+ gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
+ gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
+ break;
+ default:
+ break;
+ }
+
+ hTemp1->Sumw2();
+ hTemp2->Sumw2();
+ hTemp3->Sumw2();
+ hTemp4->Sumw2();
+ hTemp1->Add(hTemp3,-1.);
+ hTemp1->Scale(1./hTemp5->GetEntries());
+ hTemp2->Add(hTemp4,-1.);
+ hTemp2->Scale(1./hTemp6->GetEntries());
+ gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
+ gHistBalanceFunctionHistogram->Scale(0.5);
+
+ //normalize to bin width
+ gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
+ }
+
+ return gHistBalanceFunctionHistogram;
}
+
//____________________________________________________________________//
-TGraphErrors *AliBalancePsi::DrawBalance(Int_t iAnalysisType, Int_t gPsiBin) {
- // Draws the BF
- Double_t x[MAXIMUM_NUMBER_OF_STEPS];
- Double_t xer[MAXIMUM_NUMBER_OF_STEPS];
- Double_t b[MAXIMUM_NUMBER_OF_STEPS];
- Double_t ber[MAXIMUM_NUMBER_OF_STEPS];
-
- if((fNp[iAnalysisType][gPsiBin] == 0)||(fNn[iAnalysisType][gPsiBin] == 0)) {
- cerr<<"Couldn't find any particles in the analyzed interval!!!"<<endl;
- return NULL;
+TH1D *AliBalancePsi::GetBalanceFunctionHistogram2pMethod(Int_t iVariableSingle,
+ Int_t iVariablePair,
+ Double_t psiMin,
+ Double_t psiMax,
+ Double_t ptTriggerMin,
+ Double_t ptTriggerMax,
+ Double_t ptAssociatedMin,
+ Double_t ptAssociatedMax,
+ AliBalancePsi *bfMix) {
+ //Returns the BF histogram, extracted from the 6 AliTHn objects
+ //after dividing each correlation function by the Event Mixing one
+ //(private members) of the AliBalancePsi class.
+ //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
+ //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
+ // Psi_2
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+
+ // pt trigger
+ if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
}
-
- for(Int_t i = 0; i < fNumberOfBins[iAnalysisType]; i++) {
- b[i] = GetBalance(iAnalysisType,gPsiBin,i);
- ber[i] = GetError(iAnalysisType,gPsiBin,i);
- x[i] = fP2Start[iAnalysisType] + fP2Step[iAnalysisType]*i + fP2Step[iAnalysisType]/2;
- xer[i] = 0.0;
+
+ // pt associated
+ if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
+ fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
}
-
- TGraphErrors *gr = new TGraphErrors(fNumberOfBins[iAnalysisType],x,b,xer,ber);
- gr->GetXaxis()->SetTitleColor(1);
- if(iAnalysisType==0) {
- gr->SetTitle("Balance function B(#Delta y)");
- gr->GetXaxis()->SetTitle("#Delta y");
- gr->GetYaxis()->SetTitle("B(#Delta y)");
- }
- if(iAnalysisType==1) {
- gr->SetTitle("Balance function B(#Delta #eta)");
- gr->GetXaxis()->SetTitle("#Delta #eta");
- gr->GetYaxis()->SetTitle("B(#Delta #eta)");
- }
- if(iAnalysisType==2) {
- gr->SetTitle("Balance function B(q_{long})");
- gr->GetXaxis()->SetTitle("q_{long} (GeV/c)");
- gr->GetYaxis()->SetTitle("B(q_{long}) ((GeV/c)^{-1})");
- }
- if(iAnalysisType==3) {
- gr->SetTitle("Balance function B(q_{out})");
- gr->GetXaxis()->SetTitle("q_{out} (GeV/c)");
- gr->GetYaxis()->SetTitle("B(q_{out}) ((GeV/c)^{-1})");
- }
- if(iAnalysisType==4) {
- gr->SetTitle("Balance function B(q_{side})");
- gr->GetXaxis()->SetTitle("q_{side} (GeV/c)");
- gr->GetYaxis()->SetTitle("B(q_{side}) ((GeV/c)^{-1})");
- }
- if(iAnalysisType==5) {
- gr->SetTitle("Balance function B(q_{inv})");
- gr->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
- gr->GetYaxis()->SetTitle("B(q_{inv}) ((GeV/c)^{-1})");
- }
- if(iAnalysisType==6) {
- gr->SetTitle("Balance function B(#Delta #phi)");
- gr->GetXaxis()->SetTitle("#Delta #phi");
- gr->GetYaxis()->SetTitle("B(#Delta #phi)");
- }
-
- return gr;
+
+ // ============================================================================================
+ // the same for event mixing
+ AliTHn *fHistPMix = bfMix->GetHistNp();
+ AliTHn *fHistNMix = bfMix->GetHistNn();
+ AliTHn *fHistPNMix = bfMix->GetHistNpn();
+ AliTHn *fHistNPMix = bfMix->GetHistNnp();
+ AliTHn *fHistPPMix = bfMix->GetHistNpp();
+ AliTHn *fHistNNMix = bfMix->GetHistNnn();
+
+ // Psi_2
+ fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+
+ // pt trigger
+ if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
+ fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ }
+
+ // pt associated
+ if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
+ fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ }
+ // ============================================================================================
+
+ //Printf("P:%lf - N:%lf - PN:%lf - NP:%lf - PP:%lf - NN:%lf",fHistP->GetEntries(0),fHistN->GetEntries(0),fHistPN->GetEntries(0),fHistNP->GetEntries(0),fHistPP->GetEntries(0),fHistNN->GetEntries(0));
+
+ // Project into the wanted space (1st: analysis step, 2nd: axis)
+ TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair);
+ TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair);
+ TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair);
+ TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair);
+ TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle);
+ TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle);
+
+ // ============================================================================================
+ // the same for event mixing
+ TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,iVariablePair);
+ TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,iVariablePair);
+ TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,iVariablePair);
+ TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,iVariablePair);
+ TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,iVariableSingle);
+ TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,iVariableSingle);
+ // ============================================================================================
+
+ TH1D *gHistBalanceFunctionHistogram = 0x0;
+ if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
+ gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone();
+ gHistBalanceFunctionHistogram->Reset();
+
+ switch(iVariablePair) {
+ case 1:
+ gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
+ gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
+ break;
+ case 2:
+ gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
+ gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
+ break;
+ default:
+ break;
+ }
+
+ hTemp1->Sumw2();
+ hTemp2->Sumw2();
+ hTemp3->Sumw2();
+ hTemp4->Sumw2();
+ hTemp1Mix->Sumw2();
+ hTemp2Mix->Sumw2();
+ hTemp3Mix->Sumw2();
+ hTemp4Mix->Sumw2();
+
+ hTemp1->Scale(1./hTemp5->GetEntries());
+ hTemp3->Scale(1./hTemp5->GetEntries());
+ hTemp2->Scale(1./hTemp6->GetEntries());
+ hTemp4->Scale(1./hTemp6->GetEntries());
+ hTemp1Mix->Scale(1./hTemp5Mix->GetEntries());
+ hTemp3Mix->Scale(1./hTemp5Mix->GetEntries());
+ hTemp2Mix->Scale(1./hTemp6Mix->GetEntries());
+ hTemp4Mix->Scale(1./hTemp6Mix->GetEntries());
+
+ hTemp1->Divide(hTemp1Mix);
+ hTemp2->Divide(hTemp2Mix);
+ hTemp3->Divide(hTemp3Mix);
+ hTemp4->Divide(hTemp4Mix);
+
+ hTemp1->Add(hTemp3,-1.);
+ hTemp2->Add(hTemp4,-1.);
+
+ gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
+ gHistBalanceFunctionHistogram->Scale(0.5);
+
+ //normalize to bin width
+ gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
+ }
+
+ return gHistBalanceFunctionHistogram;
}
//____________________________________________________________________//
-void AliBalancePsi::PrintResults(Int_t iAnalysisType, TH1D *gHistBalance) {
- //Prints the calculated width of the BF and its error
- Double_t x[MAXIMUM_NUMBER_OF_STEPS];
- Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
- Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
- Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
- Double_t deltaBalP2 = 0.0, integral = 0.0;
- Double_t deltaErrorNew = 0.0;
-
- cout<<"=================================================="<<endl;
- for(Int_t i = 1; i <= fNumberOfBins[iAnalysisType]; i++) {
- x[i-1] = fP2Start[iAnalysisType] + fP2Step[iAnalysisType]*i + fP2Step[iAnalysisType]/2;
- //cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<gHistBalance->GetBinCenter(i)<<endl;
- }
- //cout<<"=================================================="<<endl;
- for(Int_t i = 2; i <= fNumberOfBins[iAnalysisType]; i++) {
- gSumXi += gHistBalance->GetBinCenter(i);
- gSumBi += gHistBalance->GetBinContent(i);
- gSumBiXi += gHistBalance->GetBinContent(i)*gHistBalance->GetBinCenter(i);
- gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(gHistBalance->GetBinCenter(i),2);
- gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(gHistBalance->GetBinCenter(i),2);
- gSumDeltaBi2 += TMath::Power(gHistBalance->GetBinError(i),2);
- gSumXi2DeltaBi2 += TMath::Power(gHistBalance->GetBinCenter(i),2) * TMath::Power(gHistBalance->GetBinError(i),2);
+TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi(Double_t psiMin,
+ Double_t psiMax,
+ Double_t ptTriggerMin,
+ Double_t ptTriggerMax,
+ Double_t ptAssociatedMin,
+ Double_t ptAssociatedMax) {
+ //Returns the BF histogram in Delta eta vs Delta phi,
+ //extracted from the 6 AliTHn objects
+ //(private members) of the AliBalancePsi class.
+ //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
+ //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
+ TString histName = "gHistBalanceFunctionHistogram2D";
+
+ // Psi_2
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+
+ // pt trigger
+ if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ }
+
+ // pt associated
+ if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
+ fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ }
+
+ //AliInfo(Form("P:%lf - N:%lf - PN:%lf - NP:%lf - PP:%lf - NN:%lf",fHistP->GetEntries(0),fHistN->GetEntries(0),fHistPN->GetEntries(0),fHistNP->GetEntries(0),fHistPP->GetEntries(0),fHistNN->GetEntries(0)));
+
+ // Project into the wanted space (1st: analysis step, 2nd: axis)
+ TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
+ TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
+ TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
+ TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
+ TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
+ TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
+
+ TH2D *gHistBalanceFunctionHistogram = 0x0;
+ if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
+ gHistBalanceFunctionHistogram = (TH2D*)hTemp1->Clone();
+ gHistBalanceFunctionHistogram->Reset();
+ gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
+ gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta#varphi (rad)");
+ gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
- deltaBalP2 += fP2Step[iAnalysisType]*TMath::Power(gHistBalance->GetBinError(i),2);
- integral += fP2Step[iAnalysisType]*gHistBalance->GetBinContent(i);
+ hTemp1->Sumw2();
+ hTemp2->Sumw2();
+ hTemp3->Sumw2();
+ hTemp4->Sumw2();
+ hTemp1->Add(hTemp3,-1.);
+ hTemp1->Scale(1./hTemp5->GetEntries());
+ hTemp2->Add(hTemp4,-1.);
+ hTemp2->Scale(1./hTemp6->GetEntries());
+ gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
+ gHistBalanceFunctionHistogram->Scale(0.5);
+
+ //normalize to bin width
+ gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
}
- for(Int_t i = 1; i < fNumberOfBins[iAnalysisType]; i++)
- deltaErrorNew += gHistBalance->GetBinError(i)*(gHistBalance->GetBinCenter(i)*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
-
- Double_t integralError = TMath::Sqrt(deltaBalP2);
-
- Double_t delta = gSumBiXi / gSumBi;
- Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
- cout<<"Analysis type: "<<gBFPsiAnalysisType[iAnalysisType].Data()<<endl;
- cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
- cout<<"New error: "<<deltaErrorNew<<endl;
- cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
- cout<<"=================================================="<<endl;
+
+ return gHistBalanceFunctionHistogram;
}
-
+
//____________________________________________________________________//
-TH1D *AliBalancePsi::GetBalanceFunctionHistogram(Int_t iAnalysisType,Double_t centrMin, Double_t centrMax, Double_t psiMin, Double_t psiMax) {
- //Returns the BF histogram, extracted from the 6 TH2D objects
+TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(Double_t psiMin,
+ Double_t psiMax,
+ Double_t ptTriggerMin,
+ Double_t ptTriggerMax,
+ Double_t ptAssociatedMin,
+ Double_t ptAssociatedMax,
+ AliBalancePsi *bfMix) {
+ //Returns the BF histogram in Delta eta vs Delta phi,
+ //after dividing each correlation function by the Event Mixing one
+ //extracted from the 6 AliTHn objects
//(private members) of the AliBalancePsi class.
- //
- TString gAnalysisType[ANALYSIS_TYPES] = {"y","eta","qlong","qout","qside","qinv","phi"};
- TString histName = "gHistBalanceFunctionHistogram";
- histName += gAnalysisType[iAnalysisType];
-
- SetInterval(iAnalysisType, fHistP[iAnalysisType]->GetZaxis()->GetXmin(),
- fHistP[iAnalysisType]->GetZaxis()->GetXmin(),
- fHistPP[iAnalysisType]->GetNbinsZ(),
- fHistPP[iAnalysisType]->GetZaxis()->GetXmin(),
- fHistPP[iAnalysisType]->GetZaxis()->GetXmax(),
- psiMax-psiMin);
-
- // determine the projection thresholds
- Int_t binMinX, binMinY, binMinZ;
- Int_t binMaxX, binMaxY, binMaxZ;
-
- fHistPP[iAnalysisType]->GetBinXYZ(fHistPP[iAnalysisType]->FindBin(centrMin,psiMin),binMinX,binMinY,binMinZ);
- fHistPP[iAnalysisType]->GetBinXYZ(fHistPP[iAnalysisType]->FindBin(centrMax,psiMax),binMaxX,binMaxY,binMaxZ);
-
- TH1D *gHistBalanceFunctionHistogram = new TH1D(histName.Data(),"",fHistPP[iAnalysisType]->GetNbinsZ(),fHistPP[iAnalysisType]->GetZaxis()->GetXmin(),fHistPP[iAnalysisType]->GetZaxis()->GetXmax());
- switch(iAnalysisType) {
- case kRapidity:
- gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta y");
- gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta y)");
- break;
- case kEta:
- gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta #eta");
- gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta #eta)");
- break;
- case kQlong:
- gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{long} (GeV/c)");
- gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{long})");
- break;
- case kQout:
- gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{out} (GeV/c)");
- gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{out})");
- break;
- case kQside:
- gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{side} (GeV/c)");
- gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{side})");
- break;
- case kQinv:
- gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
- gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{inv})");
- break;
- case kPhi:
- gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta #phi (deg.)");
- gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta #phi)");
- break;
- default:
- break;
- }
-
- TH1D *hTemp1 = dynamic_cast<TH1D *>(fHistPN[iAnalysisType]->ProjectionZ(Form("%s_Cent_%.0f_%.0f_Psi_%.0f_%.0f",fHistPN[iAnalysisType]->GetName(),centrMin,centrMax,psiMin,psiMax),binMinX,binMaxX,binMinY,binMaxY));
- TH1D *hTemp2 = dynamic_cast<TH1D *>(fHistPN[iAnalysisType]->ProjectionZ(Form("%s_Cent_%.0f_%.0f_Psi_%.0f_%.0f_copy",fHistPN[iAnalysisType]->GetName(),centrMin,centrMax,psiMin,psiMax),binMinX,binMaxX,binMinY,binMaxY));
- TH1D *hTemp3 = dynamic_cast<TH1D *>(fHistNN[iAnalysisType]->ProjectionZ(Form("%s_Cent_%.0f_%.0f_Psi_%.0f_%.0f",fHistNN[iAnalysisType]->GetName(),centrMin,centrMax,psiMin,psiMax),binMinX,binMaxX,binMinY,binMaxY));
- TH1D *hTemp4 = dynamic_cast<TH1D *>(fHistPP[iAnalysisType]->ProjectionZ(Form("%s_Cent_%.0f_%.0f_Psi_%.0f_%.0f",fHistPP[iAnalysisType]->GetName(),centrMin,centrMax,psiMin,psiMax),binMinX,binMaxX,binMinY,binMaxY));
- TH1D *hTemp5 = dynamic_cast<TH1D *>(fHistN[iAnalysisType]->ProjectionZ(Form("%s_Cent_%.0f_%.0f_Psi_%.0f_%.0f",fHistN[iAnalysisType]->GetName(),centrMin,centrMax,psiMin,psiMax),binMinX,binMaxX,binMinY,binMaxY));
- TH1D *hTemp6 = dynamic_cast<TH1D *>(fHistP[iAnalysisType]->ProjectionZ(Form("%s_Cent_%.0f_%.0f_Psi_%.0f_%.0f",fHistP[iAnalysisType]->GetName(),centrMin,centrMax,psiMin,psiMax),binMinX,binMaxX,binMinY,binMaxY));
-
- if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)) {
+ //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
+ //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
+ TString histName = "gHistBalanceFunctionHistogram2D";
+
+ if(!bfMix){
+ AliError("balance function object for event mixing not available");
+ return NULL;
+ }
+
+ // Psi_2
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+
+ // pt trigger
+ if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ }
+
+ // pt associated
+ if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
+ fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ }
+
+
+ // ============================================================================================
+ // the same for event mixing
+ AliTHn *fHistPMix = bfMix->GetHistNp();
+ AliTHn *fHistNMix = bfMix->GetHistNn();
+ AliTHn *fHistPNMix = bfMix->GetHistNpn();
+ AliTHn *fHistNPMix = bfMix->GetHistNnp();
+ AliTHn *fHistPPMix = bfMix->GetHistNpp();
+ AliTHn *fHistNNMix = bfMix->GetHistNnn();
+
+ // Psi_2
+ fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+
+ // pt trigger
+ if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
+ fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ }
+
+ // pt associated
+ if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
+ fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ }
+ // ============================================================================================
+
+
+ //AliInfo(Form("P:%lf - N:%lf - PN:%lf - NP:%lf - PP:%lf - NN:%lf",fHistP->GetEntries(0),fHistN->GetEntries(0),fHistPN->GetEntries(0),fHistNP->GetEntries(0),fHistPP->GetEntries(0),fHistNN->GetEntries(0)));
+
+ // Project into the wanted space (1st: analysis step, 2nd: axis)
+ TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
+ TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
+ TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
+ TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
+ TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
+ TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
+
+ // ============================================================================================
+ // the same for event mixing
+ TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,1,2);
+ TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,1,2);
+ TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,1,2);
+ TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,1,2);
+ TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,1);
+ TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,1);
+ // ============================================================================================
+
+ TH2D *gHistBalanceFunctionHistogram = 0x0;
+ if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
+ gHistBalanceFunctionHistogram = (TH2D*)hTemp1->Clone();
+ gHistBalanceFunctionHistogram->Reset();
+ gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
+ gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta#varphi (rad)");
+ gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
+
hTemp1->Sumw2();
hTemp2->Sumw2();
hTemp3->Sumw2();
hTemp4->Sumw2();
- hTemp1->Add(hTemp3,-2.);
+ hTemp1Mix->Sumw2();
+ hTemp2Mix->Sumw2();
+ hTemp3Mix->Sumw2();
+ hTemp4Mix->Sumw2();
+
hTemp1->Scale(1./hTemp5->GetEntries());
- hTemp2->Add(hTemp4,-2.);
+ hTemp3->Scale(1./hTemp5->GetEntries());
hTemp2->Scale(1./hTemp6->GetEntries());
+ hTemp4->Scale(1./hTemp6->GetEntries());
+ hTemp1Mix->Scale(1./hTemp5Mix->GetEntries());
+ hTemp3Mix->Scale(1./hTemp5Mix->GetEntries());
+ hTemp2Mix->Scale(1./hTemp6Mix->GetEntries());
+ hTemp4Mix->Scale(1./hTemp6Mix->GetEntries());
+
+ hTemp1->Divide(hTemp1Mix);
+ hTemp2->Divide(hTemp2Mix);
+ hTemp3->Divide(hTemp3Mix);
+ hTemp4->Divide(hTemp4Mix);
+
+ hTemp1->Add(hTemp3,-1.);
+ hTemp2->Add(hTemp4,-1.);
+
gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
- gHistBalanceFunctionHistogram->Scale(0.5/fP2Step[iAnalysisType]);
+ gHistBalanceFunctionHistogram->Scale(0.5);
+
+ //normalize to bin width
+ gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
}
- PrintResults(iAnalysisType,gHistBalanceFunctionHistogram);
-
return gHistBalanceFunctionHistogram;
}
+
+
+//____________________________________________________________________//
+TH2D *AliBalancePsi::GetCorrelationFunctionPN(Double_t psiMin,
+ Double_t psiMax,
+ Double_t ptTriggerMin,
+ Double_t ptTriggerMax,
+ Double_t ptAssociatedMin,
+ Double_t ptAssociatedMax) {
+ //Returns the 2D correlation function for (+-) pairs
+ // Psi_2: axis 0
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+
+ // pt trigger
+ if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ }
+
+ // pt associated
+ if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
+ fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+
+ //fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
+ //fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
+
+ //TH2D *gHistTest = dynamic_cast<TH2D *>(fHistP->Project(0,0,1));
+ //TCanvas *c1 = new TCanvas("c1","");
+ //c1->cd();
+ //if(!gHistTest){
+ //AliError("Projection of fHistP = NULL");
+ //return gHistTest;
+ //}
+ //else{
+ //gHistTest->DrawCopy("colz");
+ //}
+
+ //0:step, 1: Delta eta, 2: Delta phi
+ TH2D *gHist = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
+ if(!gHist){
+ AliError("Projection of fHistPN = NULL");
+ return gHist;
+ }
+
+ //AliInfo(Form("Entries (test): %lf",(Double_t)(gHistTest->GetEntries())));
+ //AliInfo(Form("Entries (1D): %lf",(Double_t)(fHistP->Project(0,1)->GetEntries())));
+ //AliInfo(Form("Entries (2D): %lf",(Double_t)(fHistPN->Project(0,1,2)->GetEntries())));
+
+ //TCanvas *c2 = new TCanvas("c2","");
+ //c2->cd();
+ //fHistPN->Project(0,1,2)->DrawCopy("colz");
+
+ if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
+ gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
+
+ //normalize to bin width
+ gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
+
+ return gHist;
+}
+
+//____________________________________________________________________//
+TH2D *AliBalancePsi::GetCorrelationFunctionNP(Double_t psiMin,
+ Double_t psiMax,
+ Double_t ptTriggerMin,
+ Double_t ptTriggerMax,
+ Double_t ptAssociatedMin,
+ Double_t ptAssociatedMax) {
+ //Returns the 2D correlation function for (+-) pairs
+ // Psi_2: axis 0
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+
+ // pt trigger
+ if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ }
+
+ // pt associated
+ if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
+ fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+
+ //0:step, 1: Delta eta, 2: Delta phi
+ TH2D *gHist = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
+ if(!gHist){
+ AliError("Projection of fHistPN = NULL");
+ return gHist;
+ }
+
+ //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
+ //Printf("Entries (2D): %lf",(Double_t)(fHistNP->Project(0,2,3)->GetEntries()));
+ if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
+ gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
+
+ //normalize to bin width
+ gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
+
+ return gHist;
+}
+
+//____________________________________________________________________//
+TH2D *AliBalancePsi::GetCorrelationFunctionPP(Double_t psiMin,
+ Double_t psiMax,
+ Double_t ptTriggerMin,
+ Double_t ptTriggerMax,
+ Double_t ptAssociatedMin,
+ Double_t ptAssociatedMax) {
+ //Returns the 2D correlation function for (+-) pairs
+ // Psi_2: axis 0
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+
+ // pt trigger
+ if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ }
+
+ // pt associated
+ if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
+ fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+
+ //0:step, 1: Delta eta, 2: Delta phi
+ TH2D *gHist = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
+ if(!gHist){
+ AliError("Projection of fHistPN = NULL");
+ return gHist;
+ }
+
+ //Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries()));
+ //Printf("Entries (2D): %lf",(Double_t)(fHistPP->Project(0,2,3)->GetEntries()));
+ if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
+ gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
+
+ //normalize to bin width
+ gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
+
+ return gHist;
+}
+
+//____________________________________________________________________//
+TH2D *AliBalancePsi::GetCorrelationFunctionNN(Double_t psiMin,
+ Double_t psiMax,
+ Double_t ptTriggerMin,
+ Double_t ptTriggerMax,
+ Double_t ptAssociatedMin,
+ Double_t ptAssociatedMax) {
+ //Returns the 2D correlation function for (+-) pairs
+ // Psi_2: axis 0
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+
+ // pt trigger
+ if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ }
+
+ // pt associated
+ if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
+ fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+
+ //0:step, 1: Delta eta, 2: Delta phi
+ TH2D *gHist = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
+ if(!gHist){
+ AliError("Projection of fHistPN = NULL");
+ return gHist;
+ }
+
+ //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
+ //Printf("Entries (2D): %lf",(Double_t)(fHistNN->Project(0,2,3)->GetEntries()));
+ if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
+ gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
+
+ //normalize to bin width
+ gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
+
+ return gHist;
+}
+
+//____________________________________________________________________//
+TH2D *AliBalancePsi::GetCorrelationFunctionChargeIndependent(Double_t psiMin,
+ Double_t psiMax,
+ Double_t ptTriggerMin,
+ Double_t ptTriggerMax,
+ Double_t ptAssociatedMin,
+ Double_t ptAssociatedMax) {
+ //Returns the 2D correlation function for the sum of all charge combination pairs
+ // Psi_2: axis 0
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+ fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
+
+ // pt trigger
+ if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
+ }
+
+ // pt associated
+ if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
+ fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+ fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
+
+ //0:step, 1: Delta eta, 2: Delta phi
+ TH2D *gHistNN = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
+ if(!gHistNN){
+ AliError("Projection of fHistNN = NULL");
+ return gHistNN;
+ }
+ TH2D *gHistPP = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
+ if(!gHistPP){
+ AliError("Projection of fHistPP = NULL");
+ return gHistPP;
+ }
+ TH2D *gHistNP = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
+ if(!gHistNP){
+ AliError("Projection of fHistNP = NULL");
+ return gHistNP;
+ }
+ TH2D *gHistPN = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
+ if(!gHistPN){
+ AliError("Projection of fHistPN = NULL");
+ return gHistPN;
+ }
+
+ // sum all 2 particle histograms
+ gHistNN->Add(gHistPP);
+ gHistNN->Add(gHistNP);
+ gHistNN->Add(gHistPN);
+
+ // divide by sum of + and - triggers
+ if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0 && (Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
+ gHistNN->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries() + fHistN->Project(0,1)->GetEntries()));
+
+ //normalize to bin width
+ gHistNN->Scale(1./((Double_t)gHistNN->GetXaxis()->GetBinWidth(1)*(Double_t)gHistNN->GetYaxis()->GetBinWidth(1)));
+
+ return gHistNN;
+}
+
+//____________________________________________________________________//
+Float_t AliBalancePsi::GetDPhiStar(Float_t phi1, Float_t pt1, Float_t charge1, Float_t phi2, Float_t pt2, Float_t charge2, Float_t radius, Float_t bSign) {
+ //
+ // calculates dphistar
+ //
+ Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
+
+ static const Double_t kPi = TMath::Pi();
+
+ // circularity
+// if (dphistar > 2 * kPi)
+// dphistar -= 2 * kPi;
+// if (dphistar < -2 * kPi)
+// dphistar += 2 * kPi;
+
+ if (dphistar > kPi)
+ dphistar = kPi * 2 - dphistar;
+ if (dphistar < -kPi)
+ dphistar = -kPi * 2 - dphistar;
+ if (dphistar > kPi) // might look funny but is needed
+ dphistar = kPi * 2 - dphistar;
+
+ return dphistar;
+}
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff