]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
1) correcting the normalization for per-trigger-yield for multidimensional analysis...
authormiweber <miweber@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 28 Nov 2012 19:50:28 +0000 (19:50 +0000)
committermiweber <miweber@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 28 Nov 2012 19:50:28 +0000 (19:50 +0000)
PWGCF/EBYE/BalanceFunctions/AliBalance.cxx
PWGCF/EBYE/BalanceFunctions/AliBalance.h
PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx
PWGCF/EBYE/BalanceFunctions/AliBalancePsi.h
PWGCF/EBYE/macros/drawCorrelationFunctionPsi.C
PWGCF/EBYE/macros/drawCorrelationFunctionPsiChargeIndependent.C
PWGCF/EBYE/macros/readBalanceFunction.C

index ec7d21ca58bf314bcdfc999006b6bf8ddd5844bb..064feaa0a9945952af2daacf32eb10c1b9ad7d44 100644 (file)
@@ -55,7 +55,11 @@ AliBalance::AliBalance() :
   fAnalyzedEvents(0) ,
   fCentralityId(0) ,
   fCentStart(0.),
-  fCentStop(0.)
+  fCentStop(0.),
+  fHistHBTbefore(NULL),
+  fHistHBTafter(NULL),
+  fHistConversionbefore(NULL),
+  fHistConversionafter(NULL)
 {
   // Default constructor
  
@@ -99,13 +103,6 @@ AliBalance::AliBalance() :
     fHistNN[i] = NULL;
 
   }
-
-  //QA histograms
-  fHistHBTbefore = NULL;
-  fHistHBTafter  = NULL;
-  fHistConversionbefore = NULL;
-  fHistConversionafter  = NULL;
-
 }
 
 
@@ -119,7 +116,11 @@ AliBalance::AliBalance(const AliBalance& balance):
   fAnalyzedEvents(balance.fAnalyzedEvents), 
   fCentralityId(balance.fCentralityId),
   fCentStart(balance.fCentStart),
-  fCentStop(balance.fCentStop) {
+  fCentStop(balance.fCentStop),
+  fHistHBTbefore(balance.fHistHBTbefore),
+  fHistHBTafter(balance.fHistHBTafter),
+  fHistConversionbefore(balance.fHistConversionbefore),
+  fHistConversionafter(balance.fHistConversionafter) {
   //copy constructor
   for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
     fNn[i] = balance.fNn[i];
@@ -789,12 +790,12 @@ void AliBalance::PrintResults(Int_t iAnalysisType, TH1D *gHistBalance) {
   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;
+  // 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);
@@ -814,15 +815,15 @@ void AliBalance::PrintResults(Int_t iAnalysisType, TH1D *gHistBalance) {
   
   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: "<<kBFAnalysisType[iAnalysisType].Data()<<endl;
-  cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
-  cout<<"New error: "<<deltaErrorNew<<endl;
-  cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
-  cout<<"=================================================="<<endl;
+  // cout<<"Analysis type: "<<kBFAnalysisType[iAnalysisType].Data()<<endl;
+  // cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
+  // cout<<"New error: "<<deltaErrorNew<<endl;
+  // cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
+  // cout<<"=================================================="<<endl;
 }
  
 //____________________________________________________________________//
-TH1D *AliBalance::GetBalanceFunctionHistogram(Int_t iAnalysisType,Double_t centrMin, Double_t centrMax, Double_t etaWindow,Bool_t correctWithEfficiency, Bool_t correctWithAcceptanceOnly) {
+TH1D *AliBalance::GetBalanceFunctionHistogram(Int_t iAnalysisType,Double_t centrMin, Double_t centrMax, Double_t etaWindow,Bool_t correctWithEfficiency, Bool_t correctWithAcceptanceOnly, Bool_t correctWithMixed, TH1D *hMixed[4]) {
   //Returns the BF histogram, extracted from the 6 TH2D objects 
   //(private members) of the AliBalance class.
   //
@@ -893,7 +894,7 @@ TH1D *AliBalance::GetBalanceFunctionHistogram(Int_t iAnalysisType,Double_t centr
   // withAcceptanceOnly: Data single distributions are normalized to 1 (efficiency not taken into account)
   // else : Data single distributions are normalized to give single particle efficiency of MC
   TFile *fEfficiencyMatrix = NULL;
-  if(correctWithEfficiency){
+  if(correctWithEfficiency && !correctWithMixed){
     if(correctWithAcceptanceOnly) fEfficiencyMatrix = TFile::Open("$ALICE_ROOT/PWGCF/EBYE/macros/accOnlyFromConvolutionAllCent.root");
     else  fEfficiencyMatrix = TFile::Open("$ALICE_ROOT/PWGCF/EBYE/macros/effFromConvolutionAllCent.root");
     if(!fEfficiencyMatrix){
@@ -906,7 +907,7 @@ TH1D *AliBalance::GetBalanceFunctionHistogram(Int_t iAnalysisType,Double_t centr
   // - single particle efficiencies from MC (AliAnalysiTaskEfficiency)
   // - two particle efficiencies from convolution of data single particle distributions 
   //   (normalized to single particle efficiency)
-  if(iAnalysisType == kEta && etaWindow > 0 && correctWithEfficiency){
+  if(iAnalysisType == kEta && etaWindow > 0 && correctWithEfficiency && !correctWithMixed){
 
     TH1F* hEffP  = NULL;
     TH1F* hEffN  = NULL;
@@ -974,7 +975,7 @@ TH1D *AliBalance::GetBalanceFunctionHistogram(Int_t iAnalysisType,Double_t centr
   // - single particle efficiencies from MC (AliAnalysiTaskEfficiency)
   // - two particle efficiencies from convolution of data single particle distributions 
   //   (normalized to single particle efficiency)  
-  if(iAnalysisType == kPhi && correctWithEfficiency){
+  if(iAnalysisType == kPhi && correctWithEfficiency && !correctWithMixed){
 
     TH1F* hEffPhiP  = NULL;
     TH1F* hEffPhiN  = NULL;
@@ -1020,6 +1021,48 @@ TH1D *AliBalance::GetBalanceFunctionHistogram(Int_t iAnalysisType,Double_t centr
     }  
   }
 
+  // do the correction with the event mixing directly!
+  if(correctWithMixed){
+
+    if(hMixed[0] && hMixed[1] && hMixed[2] && hMixed[3]){
+
+      // scale that EM is 1 at 0 for Deta
+      //                    in the region 0-10degree (one 1/2 sector) for Dphi
+      if(iAnalysisType==6){
+       hMixed[0]->Scale(1./(Double_t)hMixed[0]->Integral(hMixed[0]->FindBin(0),hMixed[0]->FindBin(10))*(Double_t)(hMixed[0]->FindBin(10)-hMixed[0]->FindBin(0)+1));
+       hMixed[2]->Scale(1./(Double_t)hMixed[2]->Integral(hMixed[2]->FindBin(0),hMixed[2]->FindBin(10))*(Double_t)(hMixed[0]->FindBin(10)-hMixed[0]->FindBin(0)+1));
+       hMixed[3]->Scale(1./(Double_t)hMixed[3]->Integral(hMixed[3]->FindBin(0),hMixed[3]->FindBin(10))*(Double_t)(hMixed[0]->FindBin(10)-hMixed[0]->FindBin(0)+1));
+      }
+      else{
+       hMixed[0]->Scale(1./(Double_t)hMixed[0]->GetBinContent(1));
+       hMixed[2]->Scale(1./(Double_t)hMixed[2]->GetBinContent(1));
+       hMixed[3]->Scale(1./(Double_t)hMixed[3]->GetBinContent(1));
+      }
+
+      // scale to average efficiency in the pt region (0.3-1.5) and |eta| < 0.8
+      // by multiplying the average single particle efficiencies from HIJING
+      Double_t normPMC = 0.847546;
+      Double_t normNMC = 0.83827;
+      hMixed[0]->Scale(normNMC*normPMC);
+      hMixed[2]->Scale(normNMC*normNMC);
+      hMixed[3]->Scale(normPMC*normPMC);
+
+      // divide by event mixing
+      hTemp1->Divide(hMixed[0]);
+      hTemp2->Divide(hMixed[0]);
+      hTemp3->Divide(hMixed[2]);
+      hTemp4->Divide(hMixed[3]);
+
+      // scale also single histograms with average efficiency
+      hTemp5->Scale(1./normNMC);
+      hTemp6->Scale(1./normPMC);
+
+    }
+    else{
+      AliError("Correction with EventMixing requested, but not all Histograms there!");
+      return NULL;
+    }
+  }
 
 
   if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)) {
@@ -1036,7 +1079,7 @@ TH1D *AliBalance::GetBalanceFunctionHistogram(Int_t iAnalysisType,Double_t centr
   }
 
   // do the acceptance correction (only for Eta and etaWindow > 0)
-  if(iAnalysisType == kEta && etaWindow > 0 && !correctWithEfficiency){
+  if(iAnalysisType == kEta && etaWindow > 0 && !correctWithEfficiency && !correctWithMixed){
     for(Int_t iBin = 0; iBin < gHistBalanceFunctionHistogram->GetNbinsX(); iBin++){
       
       Double_t notCorrected = gHistBalanceFunctionHistogram->GetBinContent(iBin+1);
index fb2e5e29736db7fce96ef3c0add7fa50c8b235c0..acd6e074b8c95ca21c1a408898ddb39246ebe693 100644 (file)
@@ -116,7 +116,7 @@ class AliBalance : public TObject {
   void SetHistNnn(Int_t iAnalysisType, TH2D *gHist) { 
     fHistNN[iAnalysisType] = gHist;}
 
-  TH1D *GetBalanceFunctionHistogram(Int_t iAnalysisType,Double_t centrMin, Double_t centrMax, Double_t etaWindow = -1, Bool_t correctWithEfficiency = kFALSE, Bool_t correctWithAcceptanceOnly = kFALSE);
+  TH1D *GetBalanceFunctionHistogram(Int_t iAnalysisType,Double_t centrMin, Double_t centrMax, Double_t etaWindow = -1, Bool_t correctWithEfficiency = kFALSE, Bool_t correctWithAcceptanceOnly = kFALSE, Bool_t correctWithMixed = kFALSE, TH1D *hMixed[4]=NULL);
   void PrintResults(Int_t iAnalysisType, TH1D *gHist);
 
  private:
index 96dd89b10f800d6003d25f41d01e7d2eb517afdb..02d2516d2aa2d045d23fcdfd388015cf289c292d 100644 (file)
@@ -929,8 +929,7 @@ TH2D *AliBalancePsi::GetCorrelationFunctionPN(Double_t psiMin,
                                              Double_t ptTriggerMin,
                                              Double_t ptTriggerMax,
                                              Double_t ptAssociatedMin,
-                                             Double_t ptAssociatedMax,
-                                             Bool_t   normToTrig) {
+                                             Double_t ptAssociatedMax) {
   //Returns the 2D correlation function for (+-) pairs
   // Psi_2: axis 0
   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax); 
@@ -975,7 +974,7 @@ TH2D *AliBalancePsi::GetCorrelationFunctionPN(Double_t psiMin,
   //c2->cd();
   //fHistPN->Project(0,1,2)->DrawCopy("colz");
 
-  if(normToTrig && (Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
+  if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
     gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
     
   return gHist;
@@ -987,8 +986,7 @@ TH2D *AliBalancePsi::GetCorrelationFunctionNP(Double_t psiMin,
                                              Double_t ptTriggerMin,
                                              Double_t ptTriggerMax,
                                              Double_t ptAssociatedMin,
-                                             Double_t ptAssociatedMax,
-                                             Bool_t   normToTrig) {
+                                             Double_t ptAssociatedMax) {
   //Returns the 2D correlation function for (+-) pairs
   // Psi_2: axis 0
   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax); 
@@ -1013,7 +1011,7 @@ TH2D *AliBalancePsi::GetCorrelationFunctionNP(Double_t psiMin,
 
   //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
   //Printf("Entries (2D): %lf",(Double_t)(fHistNP->Project(0,2,3)->GetEntries()));
-  if(normToTrig && (Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
+  if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
     gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
     
   return gHist;
@@ -1025,8 +1023,7 @@ TH2D *AliBalancePsi::GetCorrelationFunctionPP(Double_t psiMin,
                                              Double_t ptTriggerMin,
                                              Double_t ptTriggerMax,
                                              Double_t ptAssociatedMin,
-                                             Double_t ptAssociatedMax,
-                                             Bool_t   normToTrig) {
+                                             Double_t ptAssociatedMax) {
   //Returns the 2D correlation function for (+-) pairs
   // Psi_2: axis 0
   fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax); 
@@ -1051,7 +1048,7 @@ TH2D *AliBalancePsi::GetCorrelationFunctionPP(Double_t psiMin,
 
   //Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries()));
   //Printf("Entries (2D): %lf",(Double_t)(fHistPP->Project(0,2,3)->GetEntries()));
-  if(normToTrig && (Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
+  if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
     gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
   
   return gHist;
@@ -1063,8 +1060,7 @@ TH2D *AliBalancePsi::GetCorrelationFunctionNN(Double_t psiMin,
                                              Double_t ptTriggerMin,
                                              Double_t ptTriggerMax,
                                              Double_t ptAssociatedMin,
-                                             Double_t ptAssociatedMax,
-                                             Bool_t   normToTrig) {
+                                             Double_t ptAssociatedMax) {
   //Returns the 2D correlation function for (+-) pairs
   // Psi_2: axis 0
   fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax); 
@@ -1089,7 +1085,7 @@ TH2D *AliBalancePsi::GetCorrelationFunctionNN(Double_t psiMin,
 
   //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
   //Printf("Entries (2D): %lf",(Double_t)(fHistNN->Project(0,2,3)->GetEntries()));
-  if(normToTrig && (Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
+  if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
     gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
     
   return gHist;
@@ -1101,8 +1097,7 @@ TH2D *AliBalancePsi::GetCorrelationFunctionChargeIndependent(Double_t psiMin,
                                                             Double_t ptTriggerMin,
                                                             Double_t ptTriggerMax,
                                                             Double_t ptAssociatedMin,
-                                                            Double_t ptAssociatedMax,
-                                                            Bool_t   normToTrig) {
+                                                            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); 
@@ -1156,8 +1151,8 @@ TH2D *AliBalancePsi::GetCorrelationFunctionChargeIndependent(Double_t psiMin,
   gHistNN->Add(gHistNP);
   gHistNN->Add(gHistPN);
 
-  // if normalization to trigger particle, divide by sum of + and - triggers
-  if(normToTrig && (Double_t)(fHistN->Project(0,1)->GetEntries())!=0 && (Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
+  // 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()));
   
   return gHistNN;
index 370a3c28f932d14529cb21ff81832b4798de8983..f206c1893ffc143aeeea960ec9ab49c39ecf57bd 100644 (file)
@@ -71,32 +71,27 @@ class AliBalancePsi : public TObject {
                                   Double_t ptTriggerMin=-1.,
                                   Double_t ptTriggerMax=-1.,
                                   Double_t ptAssociatedMin=-1.,
-                                  Double_t ptAssociatedMax=-1,
-                                  Bool_t   normToTrig=kTRUE);
+                                  Double_t ptAssociatedMax=-1);
   TH2D   *GetCorrelationFunctionNP(Double_t psiMin, Double_t psiMax,
                                   Double_t ptTriggerMin=-1.,
                                   Double_t ptTriggerMax=-1.,
                                   Double_t ptAssociatedMin=-1.,
-                                  Double_t ptAssociatedMax=-1,
-                                  Bool_t   normToTrig=kTRUE);
+                                  Double_t ptAssociatedMax=-1);
   TH2D   *GetCorrelationFunctionPP(Double_t psiMin, Double_t psiMax,
                                   Double_t ptTriggerMin=-1.,
                                   Double_t ptTriggerMax=-1.,
                                   Double_t ptAssociatedMin=-1.,
-                                  Double_t ptAssociatedMax=-1,
-                                  Bool_t   normToTrig=kTRUE);
+                                  Double_t ptAssociatedMax=-1);
   TH2D   *GetCorrelationFunctionNN(Double_t psiMin, Double_t psiMax,
                                   Double_t ptTriggerMin=-1.,
                                   Double_t ptTriggerMax=-1.,
                                   Double_t ptAssociatedMin=-1.,
-                                  Double_t ptAssociatedMax=-1,
-                                  Bool_t   normToTrig=kTRUE);
+                                  Double_t ptAssociatedMax=-1);
   TH2D   *GetCorrelationFunctionChargeIndependent(Double_t psiMin, Double_t psiMax,
                                                  Double_t ptTriggerMin=-1.,
                                                  Double_t ptTriggerMax=-1.,
                                                  Double_t ptAssociatedMin=-1.,
-                                                 Double_t ptAssociatedMax=-1,
-                                                 Bool_t   normToTrig=kTRUE);
+                                                 Double_t ptAssociatedMax=-1);
 
   AliTHn *GetHistNp() {return fHistP;}
   AliTHn *GetHistNn() {return fHistN;}
index b095a2bb38464c7b63432624dc9465fdf687d01a..c93571987130c2ec55cc22e1fb043211608e81c0 100644 (file)
@@ -432,7 +432,7 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed,
       histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
     
     // if normalization to trigger then do not divide Event mixing by number of trigger particles
-    gHistPN[2] = bMixed->GetCorrelationFunctionPN(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,normToTrig);
+    gHistPN[2] = bMixed->GetCorrelationFunctionPN(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
     if(rebinEta > 1 || rebinPhi > 1) gHistPN[2]->Rebin2D(rebinEta,rebinPhi);
     
     // normalization to 1 at (0,0) --> Jan Fietes method
@@ -558,7 +558,7 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed,
       histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
 
     // if normalization to trigger then do not divide Event mixing by number of trigger particles
-    gHistNP[2] = bMixed->GetCorrelationFunctionNP(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,normToTrig);
+    gHistNP[2] = bMixed->GetCorrelationFunctionNP(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
     if(rebinEta > 1 || rebinPhi > 1) gHistNP[2]->Rebin2D(rebinEta,rebinPhi);
 
     // normalization to 1 at (0,0) --> Jan Fietes method
@@ -684,7 +684,7 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed,
       histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
     
     // if normalization to trigger then do not divide Event mixing by number of trigger particles
-    gHistPP[2] = bMixed->GetCorrelationFunctionPP(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,normToTrig);
+    gHistPP[2] = bMixed->GetCorrelationFunctionPP(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
     if(rebinEta > 1 || rebinPhi > 1) gHistPP[2]->Rebin2D(rebinEta,rebinPhi);
 
     // normalization to 1 at (0,0) --> Jan Fietes method
@@ -810,7 +810,7 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed,
       histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
     
     // if normalization to trigger then do not divide Event mixing by number of trigger particles
-    gHistNN[2] = bMixed->GetCorrelationFunctionNN(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,normToTrig);
+    gHistNN[2] = bMixed->GetCorrelationFunctionNN(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
     if(rebinEta > 1 || rebinPhi > 1) gHistNN[2]->Rebin2D(rebinEta,rebinPhi);
 
     // normalization to 1 at (0,0) --> Jan Fietes method
index 034dbc81a13fc970b81b8267e93a3b8c89c1d39f..2d5504f687539ccb4d41266ec743595d66866a6f 100644 (file)
@@ -426,7 +426,7 @@ void draw(TList *list, TList *listBFShuffled, TList *listBFMixed,
       histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})"; 
     
     // if normalization to trigger then do not divide Event mixing by number of trigger particles
-    gHist[2] = bMixed->GetCorrelationFunctionChargeIndependent(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,normToTrig);
+    gHist[2] = bMixed->GetCorrelationFunctionChargeIndependent(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
     if(rebinEta > 1 || rebinPhi > 1) gHist[2]->Rebin2D(rebinEta,rebinPhi);
     
     // normalization to 1 at (0,0) --> Jan Fietes method
index ab2b026eeaf153a1dcec2bbca07cf275def7e4bb..a88e70b9222661043892d8842973f419966b052b 100644 (file)
@@ -6,7 +6,7 @@ const Double_t centralityArray[nrOfCentralities+1] = {0.,5.,10.,20.,30.,40.,50.,
 const Double_t cent[nrOfCentralities]  = {382.8,329.7,260.5,186.4,128.9,85.,52.8,30.,15.8};   // hard coded at the moment for centrality percentiles 
 const Double_t centE[nrOfCentralities] = {3.1,4.6,4.4,3.9,3.3,2.6,2.0,1.3,0.6};               // (0-5,5-10,10-20,20-30,...,70-80)
 
-void readBalanceFunction(Bool_t bHistos = kFALSE, TString inFile = "AnalysisResults.root",Int_t fStartBinBFWidth = 3, Int_t fRebin = 2,Int_t fStartBinBFWidthPhi = 2, Int_t fRebinPhi = 2,TString centEst = "V0M",Double_t etaWindow = -1, Int_t etaBins = -1, Bool_t correctWithEfficiency = kFALSE, Bool_t correctWithAcceptanceOnly = kFALSE) {
+void readBalanceFunction(Bool_t bHistos = kFALSE, TString inFile = "AnalysisResults.root",Int_t fStartBinBFWidth = 3, Int_t fRebin = 2,Int_t fStartBinBFWidthPhi = 2, Int_t fRebinPhi = 2,TString centEst = "V0M",Double_t etaWindow = -1, Int_t etaBins = -1, Bool_t correctWithEfficiency = kFALSE, Bool_t correctWithAcceptanceOnly = kFALSE,  Bool_t correctWithMixed = kFALSE) {
   // Macro to read the output of the BF analysis:  MW: CHANGE THIS!!!!
   //i) Prints and draws the final BF output
   //ii) Plots the QA part of the analysis
@@ -20,7 +20,7 @@ void readBalanceFunction(Bool_t bHistos = kFALSE, TString inFile = "AnalysisResu
   gSystem->Load("libPWGCFebye.so");
 
   //Draw BF       
-  drawBF(bHistos,inFile, fStartBinBFWidth, fRebin,fStartBinBFWidthPhi, fRebinPhi,centEst, "",  etaWindow,etaBins, correctWithEfficiency,correctWithAcceptanceOnly);    
+  drawBF(bHistos,inFile, fStartBinBFWidth, fRebin,fStartBinBFWidthPhi, fRebinPhi,centEst, "",  etaWindow,etaBins, correctWithEfficiency,correctWithAcceptanceOnly,correctWithMixed);    
   
 
   //Merge the output
@@ -28,7 +28,7 @@ void readBalanceFunction(Bool_t bHistos = kFALSE, TString inFile = "AnalysisResu
 }
 
 //___________________________________________________________//
-void drawBF(Bool_t bHistos = kFALSE, TString inFile = "AnalysisResults.root", Int_t fStartBinBFWidth = 1, Int_t fRebin = 1, Int_t fStartBinBFWidthPhi = 1, Int_t fRebinPhi = 1, TString centEst = "V0M",TString extraString = "", Double_t etaWindow = -1, Int_t etaBins = -1, Bool_t correctWithEfficiency = kFALSE, Bool_t correctWithAcceptanceOnly = kFALSE) {
+void drawBF(Bool_t bHistos = kFALSE, TString inFile = "AnalysisResults.root", Int_t fStartBinBFWidth = 1, Int_t fRebin = 1, Int_t fStartBinBFWidthPhi = 1, Int_t fRebinPhi = 1, TString centEst = "V0M",TString extraString = "", Double_t etaWindow = -1, Int_t etaBins = -1, Bool_t correctWithEfficiency = kFALSE, Bool_t correctWithAcceptanceOnly = kFALSE,  Bool_t correctWithMixed = kFALSE) {
   //Function to draw the BF objects and write them into the output file
 
   Int_t maximumCanvases = 13;
@@ -43,6 +43,41 @@ void drawBF(Bool_t bHistos = kFALSE, TString inFile = "AnalysisResults.root", In
   TCanvas *cBF[13][10];
   TCanvas *cBFS[13][10];
 
+  // only one correction method allowed (correctWithEfficiency/Acceptance OR divide by event mixing)
+  if((correctWithEfficiency||correctWithAcceptanceOnly) && correctWithMixed){
+    Printf("only one correction method allowed (correctWithEfficiency/Acceptance OR divide by event mixing)");
+    return;
+  }
+
+  // if divide by event mixing, first get all files for all centralities
+  // here the event mixing file is more or less hard coded (has to be changed if somebody else wants to use this)
+  TH1D *hMixedEta[nrOfCentralities][4];
+  TH1D *hMixedPhi[nrOfCentralities][4];
+  
+  if(correctWithMixed){
+    TFile* fMixed[nrOfCentralities];
+    
+    for (Int_t iCent = 0; iCent < nrOfCentralities ; iCent++){
+      fMixed[iCent] = TFile::Open(Form("/Users/physics/ALICE/balance/BF_EM_new/gridOut/Centrality%.0f_%.0f_oldMethod/Histograms_WMstart%d_rebin%d_WMstartPhi%d_rebinPhi%d_EventMixing_%.0f-%.0f_AnalysisResults.root",centralityArray[iCent],centralityArray[iCent+1],fStartBinBFWidth, fRebin,fStartBinBFWidthPhi, fRebinPhi,centralityArray[iCent],centralityArray[iCent+1]),"READ");
+
+      if(!fMixed[iCent]){
+       Printf("Event Mixing file not available!");
+       return;
+      }
+
+      // at the time of production the event mixing had no centrality info (all in most central bin, centarlity is selected before)
+      hMixedEta[iCent][0] =  (TH1D*)fMixed[iCent]->Get(Form("hPN_eta_0"));
+      hMixedPhi[iCent][0] =  (TH1D*)fMixed[iCent]->Get(Form("hPN_phi_0"));
+      hMixedEta[iCent][1] =  (TH1D*)fMixed[iCent]->Get(Form("hNP_eta_0"));
+      hMixedPhi[iCent][1] =  (TH1D*)fMixed[iCent]->Get(Form("hNP_phi_0"));
+      hMixedEta[iCent][2] =  (TH1D*)fMixed[iCent]->Get(Form("hNN_eta_0"));
+      hMixedPhi[iCent][2] =  (TH1D*)fMixed[iCent]->Get(Form("hNN_phi_0"));
+      hMixedEta[iCent][3] =  (TH1D*)fMixed[iCent]->Get(Form("hPP_eta_0"));
+      hMixedPhi[iCent][3] =  (TH1D*)fMixed[iCent]->Get(Form("hPP_phi_0"));
+      
+    }
+  }
+
   // get the file
   TFile *f = TFile::Open(inFile.Data());
   if(!f) {
@@ -117,7 +152,7 @@ void drawBF(Bool_t bHistos = kFALSE, TString inFile = "AnalysisResults.root", In
     list = (TList*)key->ReadObj();
     listName = TString(list->GetName());
 
-    cout<<"Processing list "<<listName<<endl;
+    //cout<<"Processing list "<<listName<<endl;
 
  
     // ----------------------------------------------------
@@ -281,7 +316,7 @@ void drawBF(Bool_t bHistos = kFALSE, TString inFile = "AnalysisResults.root", In
 
       for(Int_t a = 0; a < 7; a++){
 
-       cout<<"ANALYSE "<<gBFAnalysisType[a]<<endl;
+       //cout<<"ANALYSE "<<gBFAnalysisType[a]<<endl;
 
        // create the BF object
        bf[iList][a]  = new AliBalance();
@@ -327,8 +362,11 @@ void drawBF(Bool_t bHistos = kFALSE, TString inFile = "AnalysisResults.root", In
        bf[iList][a]->SetHistNnn(a, fHistNN[iList][a]);
 
        for(iCanvas = 0; iCanvas < nrOfCentralities; iCanvas++){
+         
+         if(correctWithMixed && a == 1) gbf[iList][iCanvas][a] = bf[iList][a]->GetBalanceFunctionHistogram(a,centralityArray[iCanvas],centralityArray[iCanvas+1],etaWindow,correctWithEfficiency,correctWithAcceptanceOnly,correctWithMixed,hMixedEta[iCanvas]);
+         else if(correctWithMixed && a == 6) gbf[iList][iCanvas][a] = bf[iList][a]->GetBalanceFunctionHistogram(a,centralityArray[iCanvas],centralityArray[iCanvas+1],etaWindow,correctWithEfficiency,correctWithAcceptanceOnly,correctWithMixed,hMixedPhi[iCanvas]);
+         else gbf[iList][iCanvas][a] = bf[iList][a]->GetBalanceFunctionHistogram(a,centralityArray[iCanvas],centralityArray[iCanvas+1],etaWindow,correctWithEfficiency,correctWithAcceptanceOnly);
 
-         gbf[iList][iCanvas][a] = bf[iList][a]->GetBalanceFunctionHistogram(a,centralityArray[iCanvas],centralityArray[iCanvas+1],etaWindow,correctWithEfficiency,correctWithAcceptanceOnly);
          gbf[iList][iCanvas][a]->SetName(Form("BF_%s_Cent_%.0f_%.0f_%d",gBFAnalysisType[a].Data(),centralityArray[iCanvas],centralityArray[iCanvas+1],iList));
 
          cBF[iList][iCanvas]->cd(a+1);
@@ -431,7 +469,9 @@ void drawBF(Bool_t bHistos = kFALSE, TString inFile = "AnalysisResults.root", In
 
        for(iCanvas = 0; iCanvas < nrOfCentralities; iCanvas++){
          
-         gbfs[iList][iCanvas][a] = bfs[iList][a]->GetBalanceFunctionHistogram(a,centralityArray[iCanvas],centralityArray[iCanvas+1],etaWindow, correctWithEfficiency,correctWithAcceptanceOnly);
+         if(correctWithMixed && a == 1) gbfs[iList][iCanvas][a] = bfs[iList][a]->GetBalanceFunctionHistogram(a,centralityArray[iCanvas],centralityArray[iCanvas+1],etaWindow, correctWithEfficiency,correctWithAcceptanceOnly,correctWithMixed,hMixedEta[iCanvas]);
+         else if(correctWithMixed && a == 6) gbfs[iList][iCanvas][a] = bfs[iList][a]->GetBalanceFunctionHistogram(a,centralityArray[iCanvas],centralityArray[iCanvas+1],etaWindow, correctWithEfficiency,correctWithAcceptanceOnly,correctWithMixed,hMixedPhi[iCanvas]);
+         else gbfs[iList][iCanvas][a] = bfs[iList][a]->GetBalanceFunctionHistogram(a,centralityArray[iCanvas],centralityArray[iCanvas+1],etaWindow, correctWithEfficiency,correctWithAcceptanceOnly);
          gbfs[iList][iCanvas][a]->SetName(Form("BFS_%s_Cent_%.0f_%.0f_%d",gBFAnalysisType[a].Data(),centralityArray[iCanvas],centralityArray[iCanvas+1],iList));
          
          cBFS[iList][iCanvas]->cd(a+1);
@@ -541,6 +581,9 @@ void drawBF(Bool_t bHistos = kFALSE, TString inFile = "AnalysisResults.root", In
        fOut = TFile::Open(Form("Histograms_AccCorrWithEfficiency_Window%.1f_Bins%d_WMstart%d_rebin%d_WMstartPhi%d_rebinPhi%d_%s", etaWindow, etaBins, fStartBinBFWidth, fRebin,fStartBinBFWidthPhi, fRebinPhi,inFile.Data()),"RECREATE");
       }
     }
+    else if(correctWithMixed){
+      fOut = TFile::Open(Form("Histograms_CorrEM_Window%.1f_Bins%d_WMstart%d_rebin%d_WMstartPhi%d_rebinPhi%d_%s", etaWindow, etaBins, fStartBinBFWidth, fRebin,fStartBinBFWidthPhi, fRebinPhi,inFile.Data()),"RECREATE");
+    }
     else{
       fOut = TFile::Open(Form("Histograms_AccCorr_Window%.1f_Bins%d_WMstart%d_rebin%d_WMstartPhi%d_rebinPhi%d_%s", etaWindow, etaBins, fStartBinBFWidth, fRebin,fStartBinBFWidthPhi, fRebinPhi,inFile.Data()),"RECREATE");
     }
@@ -550,9 +593,11 @@ void drawBF(Bool_t bHistos = kFALSE, TString inFile = "AnalysisResults.root", In
   }
   fOut->cd();
   for(Int_t i = 0; i < iList+1; i++){
+    if(gbf[i][0][1]){
     cout<<"PROCESS LIST "<<i<<" NOW!"<<endl;
     for(Int_t a = 0; a < 7; a++){
-    cout<<"PROCESS VARIABLE "<<a<<" NOW!"<<endl;
+      if(a==1||a==6){
+      cout<<"PROCESS VARIABLE "<<a<<" NOW!"<<endl;
       
       if(fHistPN[i][a]){
        (fHistPN[i][a]->ProjectionY(Form("hPN_%s_%d",gBFAnalysisType[a].Data(),i)))->Write();
@@ -619,6 +664,7 @@ void drawBF(Bool_t bHistos = kFALSE, TString inFile = "AnalysisResults.root", In
        cout<<"//=========================Centrality "<<centralityArray[j]<<"-"<<centralityArray[j+1]<<"%==================//"<<endl;
        cout<<endl;
       }
+      }
     }
 
     Double_t x,y;
@@ -765,6 +811,7 @@ void drawBF(Bool_t bHistos = kFALSE, TString inFile = "AnalysisResults.root", In
       gintegPS[i]->Write();
 
     }
+    }
   }
   fOut->Close();
   f->Close();
@@ -786,13 +833,13 @@ void GetWeightedMean(TH1D *gHistBalance, Int_t fStartBin = 1, Int_t fStopBin = -
   if(fStopBin > -1) fNumberOfBins = fStopBin;
   Double_t fP2Step    = gHistBalance->GetBinWidth(1); // assume equal binning!
   
-  cout<<"=================================================="<<endl;
-  cout<<"RECALCULATION OF BF WIDTH (StartBin = "<<fStartBin<<")"<<endl;
-  cout<<"HISTOGRAM has "<<fNumberOfBins<<" bins with bin size of "<<fP2Step<<endl;
-  for(Int_t i = fStartBin; i <= fNumberOfBins; i++) { 
-    cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<gHistBalance->GetBinCenter(i)<<endl;
-  } 
-  cout<<"=================================================="<<endl;
+  // cout<<"=================================================="<<endl;
+  // cout<<"RECALCULATION OF BF WIDTH (StartBin = "<<fStartBin<<")"<<endl;
+  // cout<<"HISTOGRAM has "<<fNumberOfBins<<" bins with bin size of "<<fP2Step<<endl;
+  // for(Int_t i = fStartBin; i <= fNumberOfBins; i++) { 
+  //   cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<gHistBalance->GetBinCenter(i)<<endl;
+  // 
+  // cout<<"=================================================="<<endl;
   for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {
     gSumXi += gHistBalance->GetBinCenter(i);
     gSumBi += gHistBalance->GetBinContent(i);
@@ -813,10 +860,10 @@ void GetWeightedMean(TH1D *gHistBalance, Int_t fStartBin = 1, Int_t fStopBin = -
   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<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
-  cout<<"New error: "<<deltaErrorNew<<endl;
-  cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
-  cout<<"=================================================="<<endl;
+  // cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
+  // cout<<"New error: "<<deltaErrorNew<<endl;
+  // cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
+  // cout<<"=================================================="<<endl;
 
   WM  = delta;
   WME = deltaError;
@@ -839,13 +886,13 @@ void GetIntegral(TH1D *gHistBalance, Double_t &integ, Double_t &integE) {
   Int_t fNumberOfBins = gHistBalance->GetNbinsX();
   Double_t fP2Step    = gHistBalance->GetBinWidth(1); // assume equal binning!
   
-  cout<<"=================================================="<<endl;
-  cout<<"RECALCULATION OF BF WIDTH (StartBin = "<<fStartBin<<")"<<endl;
-  cout<<"HISTOGRAM has "<<fNumberOfBins<<" bins with bin size of "<<fP2Step<<endl;
-  for(Int_t i = fStartBin; i <= fNumberOfBins; i++) { 
-    cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<gHistBalance->GetBinCenter(i)<<endl;
-  } 
-  cout<<"=================================================="<<endl;
+  // cout<<"=================================================="<<endl;
+  // cout<<"RECALCULATION OF BF WIDTH (StartBin = "<<fStartBin<<")"<<endl;
+  // cout<<"HISTOGRAM has "<<fNumberOfBins<<" bins with bin size of "<<fP2Step<<endl;
+  // for(Int_t i = fStartBin; i <= fNumberOfBins; i++) { 
+  //   cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<gHistBalance->GetBinCenter(i)<<endl;
+  // 
+  // cout<<"=================================================="<<endl;
   for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {
     gSumXi += gHistBalance->GetBinCenter(i);
     gSumBi += gHistBalance->GetBinContent(i);
@@ -866,10 +913,10 @@ void GetIntegral(TH1D *gHistBalance, Double_t &integ, Double_t &integE) {
   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<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
-  cout<<"New error: "<<deltaErrorNew<<endl;
-  cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
-  cout<<"=================================================="<<endl;
+  // cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
+  // cout<<"New error: "<<deltaErrorNew<<endl;
+  // cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
+  // cout<<"=================================================="<<endl;
 
   integ  = integral;
   integE = integralError;