]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx
HBT like cuts for all pairs (before only like sign)
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliBalancePsi.cxx
index 02d2516d2aa2d045d23fcdfd388015cf289c292d..49ab29f7d6b29715d1fd5c0e7f4d1412e7de175f 100644 (file)
@@ -37,6 +37,7 @@
 #include "AliESDtrack.h"
 #include "AliAODTrack.h"
 #include "AliTHn.h"
+#include "AliAnalysisTaskTriggeredBF.h"
 
 #include "AliBalancePsi.h"
 
@@ -63,8 +64,10 @@ AliBalancePsi::AliBalancePsi() :
   fHistConversionafter(0),
   fHistPsiMinusPhi(0),
   fPsiInterval(15.),
+  fDeltaEtaMax(2.0),
   fHBTCut(kFALSE),
-  fConversionCut(kFALSE) {
+  fConversionCut(kFALSE),
+  fEventClass("EventPlane"){
   // Default constructor
 }
 
@@ -88,8 +91,10 @@ AliBalancePsi::AliBalancePsi(const AliBalancePsi& balance):
   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
 }
 
@@ -108,6 +113,7 @@ AliBalancePsi::~AliBalancePsi() {
   delete fHistConversionbefore;
   delete fHistConversionafter;
   delete fHistPsiMinusPhi;
+    
 }
 
 //____________________________________________________________________//
@@ -128,32 +134,72 @@ void AliBalancePsi::InitHistograms() {
   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"; 
@@ -264,10 +310,11 @@ void AliBalancePsi::InitHistograms() {
 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");
@@ -295,22 +342,28 @@ void AliBalancePsi::CalculateBalance(Double_t gReactionPlane,
   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.;
@@ -338,26 +391,24 @@ void AliBalancePsi::CalculateBalance(Double_t gReactionPlane,
     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 
@@ -376,7 +427,8 @@ void AliBalancePsi::CalculateBalance(Double_t gReactionPlane,
       //       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;
        
@@ -461,10 +513,10 @@ void AliBalancePsi::CalculateBalance(Double_t gReactionPlane,
        }
       }//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;
@@ -550,6 +602,9 @@ TH1D *AliBalancePsi::GetBalanceFunctionHistogram(Int_t iVariableSingle,
     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;
@@ -698,6 +753,9 @@ TH1D *AliBalancePsi::GetBalanceFunctionHistogram2pMethod(Int_t iVariableSingle,
 
     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;
@@ -771,6 +829,9 @@ TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi(Double_t psiMin,
     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;
@@ -917,6 +978,9 @@ TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(Double_t psiMin,
 
     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;
@@ -976,6 +1040,9 @@ TH2D *AliBalancePsi::GetCorrelationFunctionPN(Double_t psiMin,
 
   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;
 }
@@ -1013,6 +1080,9 @@ TH2D *AliBalancePsi::GetCorrelationFunctionNP(Double_t psiMin,
   //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;
 }
@@ -1050,6 +1120,9 @@ TH2D *AliBalancePsi::GetCorrelationFunctionPP(Double_t psiMin,
   //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;
 }
@@ -1087,6 +1160,9 @@ TH2D *AliBalancePsi::GetCorrelationFunctionNN(Double_t psiMin,
   //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;
 }
@@ -1154,6 +1230,9 @@ TH2D *AliBalancePsi::GetCorrelationFunctionChargeIndependent(Double_t psiMin,
   // 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;
 }