#include "AliESDtrack.h"
#include "AliAODTrack.h"
#include "AliTHn.h"
+#include "AliAnalysisTaskTriggeredBF.h"
#include "AliBalancePsi.h"
fHistConversionafter(0),
fHistPsiMinusPhi(0),
fPsiInterval(15.),
+ fDeltaEtaMax(2.0),
fHBTCut(kFALSE),
- fConversionCut(kFALSE) {
+ fConversionCut(kFALSE),
+ fEventClass("EventPlane"){
// Default constructor
}
fHistConversionafter(balance.fHistConversionafter),
fHistPsiMinusPhi(balance.fHistPsiMinusPhi),
fPsiInterval(balance.fPsiInterval),
+ fDeltaEtaMax(balance.fDeltaEtaMax),
fHBTCut(balance.fHBTCut),
- fConversionCut(balance.fConversionCut) {
+ fConversionCut(balance.fConversionCut),
+ fEventClass("EventPlane"){
//copy constructor
}
delete fHistConversionbefore;
delete fHistConversionafter;
delete fHistPsiMinusPhi;
+
}
//____________________________________________________________________//
Int_t iBinPair[kTrackVariablesPair]; // binning for track variables
Double_t* dBinsPair[kTrackVariablesPair]; // bins for track variables
TString axisTitlePair[kTrackVariablesPair]; // axis titles for track variables
-
- //centrality
- /*const Int_t kNCentralityBins = 9;
- Double_t centralityBins[kNCentralityBins+1] = {0.,5.,10.,20.,30.,40.,50.,60.,70.,80.};
- iBinSingle[0] = kNCentralityBins;
- dBinsSingle[0] = centralityBins;
- axisTitleSingle[0] = "Centrality percentile [%]";
- iBinPair[0] = kNCentralityBins;
- dBinsPair[0] = centralityBins;
- axisTitlePair[0] = "Centrality percentile [%]"; */
-
- //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};
- 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.)";
+ /**********************************************************
+
+ ======> 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};
+ //----------------------------------------------------------
+
+ //--- Centrality Bins --------------------------------------
+ const Int_t kNCentralityBins = 9;
+ Double_t centralityBins[kNCentralityBins+1] = {0.,5.,10.,20.,30.,40.,50.,60.,70.,80.};
+ //----------------------------------------------------------
+
+ //--- 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};
+ //----------------------------------------------------------
+
+ //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] = -2.0 + i * 0.05;
+ deltaEtaBins[i] = - fDeltaEtaMax + i * 2 * fDeltaEtaMax / (Double_t)kNDeltaEtaBins;
iBinPair[1] = kNDeltaEtaBins;
dBinsPair[1] = deltaEtaBins;
axisTitlePair[1] = "#Delta#eta";
void AliBalancePsi::CalculateBalance(Double_t gReactionPlane,
TObjArray *particles,
TObjArray *particlesMixed,
- Float_t bSign) {
+ Float_t bSign,
+ Double_t kMultorCent) {
// Calculates the balance function
fAnalyzedEvents++;
-
+
// Initialize histograms if not done yet
if(!fHistPN){
AliWarning("Histograms not yet initialized! --> Will be done now");
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
}
// 1st particle loop
for (Int_t i = 0; i < iMax; i++) {
- AliVParticle* firstParticle = (AliVParticle*) particles->At(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
+
// Event plane (determine psi bin)
Double_t gPsiMinusPhi = 0.;
Double_t gPsiMinusPhiBin = -10.;
Short_t charge1 = (Short_t) firstParticle->Charge();
trackVariablesSingle[0] = gPsiMinusPhiBin;
- trackVariablesSingle[1] = firstPt;
+ trackVariablesSingle[1] = firstPt;
+ if(fEventClass=="Multiplicity" || fEventClass == "Centrality" ) trackVariablesSingle[0] = kMultorCent;
//fill single particle histograms
- if(charge1 > 0) fHistP->Fill(trackVariablesSingle,0,1.);
- else if(charge1 < 0) fHistN->Fill(trackVariablesSingle,0,1.);
+ if(charge1 > 0) fHistP->Fill(trackVariablesSingle,0,firstCorrection); //==========================correction
+ else if(charge1 < 0) fHistN->Fill(trackVariablesSingle,0,firstCorrection); //==========================correction
- // 2nd particle loop (only for j < i for non double counting in the same pT region)
- // --> SAME pT region for trigger and assoc: NO double counting with this
- // --> DIFF pT region for trigger and assoc: Missing assoc. particles with j > i to a trigger i
- // --> can be handled afterwards by using assoc. as trigger as well ?!
- for(Int_t j = 0; j < iMax; j++) {
- if (particlesMixed && j == jMax-1 ) // if the mixed track number is smaller than the main event one
- break;
-
- if(j == i) continue; // no auto correlations
-
- AliVParticle* secondParticle = (AliVParticle*) particlesSecond->At(j);
- Short_t charge2 = (Short_t) secondParticle->Charge();
+ // 2nd particle loop
+ for(Int_t j = 0; j < jMax; j++) {
+
+ if(!particlesMixed && j == i) continue; // no auto correlations (only for non mixing)
+
+ // pT,Assoc < pT,Trig
+ if(firstPt < secondPt[j]) continue;
+
+ Short_t charge2 = secondCharge[j];
- trackVariablesPair[0] = gPsiMinusPhiBin;
+ 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[5] = fCentrality; // centrality
// HBT like cut
- if(fHBTCut && charge1 * charge2 > 0){
+ 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;
}
}//conversion cut
- if( charge1 > 0 && charge2 < 0) fHistPN->Fill(trackVariablesPair,0,1.);
- else if( charge1 < 0 && charge2 > 0) fHistNP->Fill(trackVariablesPair,0,1.);
- else if( charge1 > 0 && charge2 > 0) fHistPP->Fill(trackVariablesPair,0,1.);
- else if( charge1 < 0 && charge2 < 0) fHistNN->Fill(trackVariablesPair,0,1.);
+ 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;
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;
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;
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;
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;
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;
}
//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;
}
//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;
}
//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;
}
// 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;
}