]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
correlation function calculation uses averages of ratios also for centrality
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 25 Apr 2012 08:13:39 +0000 (08:13 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 25 Apr 2012 08:13:39 +0000 (08:13 +0000)
ttr cut value configurable
coverity fix

PWGCF/Correlations/Base/AliUEHist.cxx
PWGCF/Correlations/Base/AliUEHist.h
PWGCF/Correlations/Base/AliUEHistograms.cxx
PWGCF/Correlations/Base/AliUEHistograms.h
PWGCF/Correlations/DPhi/AliAnalysisTaskPhiCorrelations.cxx
PWGCF/Correlations/DPhi/AliAnalysisTaskPhiCorrelations.h
PWGCF/Correlations/macros/dphicorrelations/correct.C
PWGCF/Correlations/macros/dphicorrelations/fit.C

index 496159d5b27385e262639d41795ea820b92ccf00..4f2a145e717f406151784b83eed435929e9a2a24 100644 (file)
@@ -35,6 +35,7 @@
 #include "TCanvas.h"
 #include "TF1.h"
 #include "AliTHn.h"
+#include "THn.h"
 
 ClassImp(AliUEHist)
 
@@ -56,6 +57,8 @@ AliUEHist::AliUEHist(const char* reqHist) :
   fContaminationEnhancement(0),
   fCombineMinMax(0),
   fCache(0),
+  fGetMultCacheOn(kFALSE),
+  fGetMultCache(0),
   fHistogramType(reqHist)
 {
   // Constructor
@@ -324,6 +327,8 @@ AliUEHist::AliUEHist(const AliUEHist &c) :
   fContaminationEnhancement(0),
   fCombineMinMax(0),
   fCache(0),
+  fGetMultCacheOn(kFALSE),
+  fGetMultCache(0),
   fHistogramType()
 {
   //
@@ -474,26 +479,42 @@ Long64_t AliUEHist::Merge(TCollection* list)
 }
 
 //____________________________________________________________________
-void AliUEHist::SetBinLimits(AliCFGridSparse* grid)
+void AliUEHist::SetBinLimits(THnBase* grid)
 {
   // sets the bin limits in eta and pT defined by fEtaMin/Max, fPtMin/Max
   
   if (fEtaMax > fEtaMin)
-    grid->SetRangeUser(0, fEtaMin, fEtaMax);
+    grid->GetAxis(0)->SetRangeUser(fEtaMin, fEtaMax);
   if (fPtMax > fPtMin)
-    grid->SetRangeUser(1, fPtMin, fPtMax);
+    grid->GetAxis(1)->SetRangeUser(fPtMin, fPtMax);
 }  
 
+//____________________________________________________________________
+void AliUEHist::SetBinLimits(AliCFGridSparse* grid)
+{
+  // sets the bin limits in eta and pT defined by fEtaMin/Max, fPtMin/Max
+  
+  SetBinLimits(grid->GetGrid());
+}
+
+//____________________________________________________________________
+void AliUEHist::ResetBinLimits(THnBase* grid)
+{
+  // resets all bin limits 
+
+  for (Int_t i=0; i<grid->GetNdimensions(); i++)
+    if (grid->GetAxis(i)->TestBit(TAxis::kAxisRange))
+      grid->GetAxis(i)->SetRangeUser(0, -1);
+}
+
 //____________________________________________________________________
 void AliUEHist::ResetBinLimits(AliCFGridSparse* grid)
 {
   // resets all bin limits 
   
-  for (Int_t i=0; i<grid->GetNVar(); i++)
-    if (grid->GetGrid()->GetAxis(i)->TestBit(TAxis::kAxisRange))
-      grid->SetRangeUser(i, 0, -1);
+  ResetBinLimits(grid->GetGrid());
 }
-  
+
 //____________________________________________________________________
 void AliUEHist::CountEmptyBins(AliUEHist::CFStep step, Float_t ptLeadMin, Float_t ptLeadMax)
 {
@@ -758,12 +779,57 @@ void AliUEHist::GetHistsZVtx(AliUEHist::CFStep step, AliUEHist::Region region, F
   ResetBinLimits(fEventHist->GetGrid(step));
 }
 
+//____________________________________________________________________
+void AliUEHist::GetHistsZVtxMult(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, THnBase** trackHist, TH2** eventHist)
+{
+  // Calculates a 4d histogram with deltaphi, deltaeta, zvtx, multiplicity on track level and 
+  // a 2d histogram on event level (as fct of zvtx, multiplicity)
+  // Histograms has to be deleted by the caller of the function
+  
+  if (!fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(5))
+    AliFatal("Histogram without vertex axis provided");
+
+  THnBase* sparse = fTrackHist[region]->GetGrid(step)->GetGrid();
+  if (fGetMultCacheOn)
+  {
+    if (!fGetMultCache)
+    {
+      fGetMultCache = ChangeToThn(sparse);
+      // should work but causes SEGV in ProjectionND below 
+//       delete sparse;
+    }
+    sparse = fGetMultCache;
+  }
+  
+  // unzoom all axes
+  ResetBinLimits(sparse);
+  ResetBinLimits(fEventHist->GetGrid(step));
+  
+  SetBinLimits(sparse);
+  
+  Int_t firstBin = sparse->GetAxis(2)->FindBin(ptLeadMin);
+  Int_t lastBin = sparse->GetAxis(2)->FindBin(ptLeadMax);
+  Printf("Using leading pT range %d --> %d", firstBin, lastBin);
+  sparse->GetAxis(2)->SetRange(firstBin, lastBin);
+  fEventHist->GetGrid(step)->GetGrid()->GetAxis(0)->SetRange(firstBin, lastBin);
+    
+  Int_t dimensions[] = { 4, 0, 5, 3 };
+  THnBase* tmpTrackHist = sparse->ProjectionND(4, dimensions, "E");
+  *eventHist = (TH2*) fEventHist->GetGrid(step)->Project(2, 1);
+  
+  ResetBinLimits(sparse);
+  ResetBinLimits(fEventHist->GetGrid(step));
+
+  // convert to THn 
+  *trackHist = ChangeToThn(tmpTrackHist);
+}
+
 //____________________________________________________________________
 TH2* AliUEHist::GetSumOfRatios2(AliUEHist* mixed, AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd)
 {
-  // Calls GetUEHist(...) for *each* vertex bin and performs a sum of ratios:
+  // Calls GetUEHist(...) for *each* vertex bin and multiplicity bin and performs a sum of ratios:
   // 1_N [ (same/mixed)_1 + (same/mixed)_2 + (same/mixed)_3 + ... ]
-  // where N is the total number of events/trigger particles and the subscript is the vertex bin
+  // where N is the total number of events/trigger particles and the subscript is the vertex/multiplicity bin
   // where mixed is normalized such that the information about the number of pairs in same is kept
   //
   // returns a 2D histogram: deltaphi, deltaeta
@@ -772,69 +838,112 @@ TH2* AliUEHist::GetSumOfRatios2(AliUEHist* mixed, AliUEHist::CFStep step, AliUEH
   //   mixed: AliUEHist containing mixed event corresponding to this object
   //   <other parameters> : check documentation of AliUEHist::GetUEHist
   
+  // do not add this hists to the directory
+  Bool_t oldStatus = TH1::AddDirectoryStatus();
+  TH1::AddDirectory(kFALSE);
+
   TH2* totalTracks = 0;
   
-  TH3* trackSameAll = 0;
-  TH3* trackMixedAll = 0;
-  TH1* eventSameAll = 0;
-  TH1* eventMixedAll = 0;
+  THnBase* trackSameAll = 0;
+  THnBase* trackMixedAll = 0;
+  TH2* eventSameAll = 0;
+  TH2* eventMixedAll = 0;
   
   Int_t totalEvents = 0;
   
-  GetHistsZVtx(step, region, ptLeadMin, ptLeadMax, multBinBegin, multBinEnd, &trackSameAll, &eventSameAll);
-  mixed->GetHistsZVtx(step, region, ptLeadMin, ptLeadMax, multBinBegin, multBinEnd, &trackMixedAll, &eventMixedAll);
+  GetHistsZVtxMult(step, region, ptLeadMin, ptLeadMax, &trackSameAll, &eventSameAll);
+  mixed->GetHistsZVtxMult(step, region, ptLeadMin, ptLeadMax, &trackMixedAll, &eventMixedAll);
+  
+//   TH1* normParameters = new TH1F("normParameters", "", 100, 0, 2);
   
-  TAxis* vertexAxis = trackSameAll->GetZaxis();
-  for (Int_t vertexBin = 1; vertexBin <= vertexAxis->GetNbins(); vertexBin++)
+//   trackSameAll->Dump();
+  
+  TAxis* multAxis = trackSameAll->GetAxis(3);
+  for (Int_t multBin = TMath::Max(1, multBinBegin); multBin <= TMath::Min(multAxis->GetNbins(), multBinEnd); multBin++)
   {
-    trackSameAll->GetZaxis()->SetRange(vertexBin, vertexBin);
-    trackMixedAll->GetZaxis()->SetRange(vertexBin, vertexBin);
+    trackSameAll->GetAxis(3)->SetRange(multBin, multBin);
+    trackMixedAll->GetAxis(3)->SetRange(multBin, multBin);
 
-    TH2* tracksSame = (TH2*) trackSameAll->Project3D("yx1");
-    TH2* tracksMixed = (TH2*) trackMixedAll->Project3D("yx2");
+    // get mixed normalization correction factor: is independent of vertex bin if scaled with number of triggers
+    trackMixedAll->GetAxis(2)->SetRange(0, -1);
+    TH2* tracksMixed = trackMixedAll->Projection(1, 0, "E");
     
-    // asssume flat in dphi, gain in statistics
-//     TH1* histMixedproj = mixedTwoD->ProjectionY();
-//     histMixedproj->Scale(1.0 / mixedTwoD->GetNbinsX());
-//     
-//     for (Int_t x=1; x<=mixedTwoD->GetNbinsX(); x++)
-//       for (Int_t y=1; y<=mixedTwoD->GetNbinsY(); y++)
-//     mixedTwoD->SetBinContent(x, y, histMixedproj->GetBinContent(y));
-
-//       delete histMixedproj;
+    // get mixed event normalization by assuming full acceptance at deta at 0 (integrate over dphi)
+    Double_t mixedNormError = 0;
+    Double_t mixedNorm = tracksMixed->IntegralAndError(1, tracksMixed->GetNbinsX(), tracksMixed->GetYaxis()->FindBin(-0.01), tracksMixed->GetYaxis()->FindBin(0.01), mixedNormError);
+    Int_t nBinsMixedNorm = (tracksMixed->GetNbinsX()) * (tracksMixed->GetYaxis()->FindBin(0.01) - tracksMixed->GetYaxis()->FindBin(-0.01) + 1);
+/*    Double_t mixedNorm = tracksMixed->IntegralAndError(tracksMixed->GetXaxis()->FindBin(-0.01), tracksMixed->GetXaxis()->FindBin(0.01), tracksMixed->GetYaxis()->FindBin(-0.01), tracksMixed->GetYaxis()->FindBin(0.01), mixedNormError);
+    Int_t nBinsMixedNorm = (tracksMixed->GetXaxis()->FindBin(0.01) - tracksMixed->GetXaxis()->FindBin(-0.01) + 1) * (tracksMixed->GetYaxis()->FindBin(0.01) - tracksMixed->GetYaxis()->FindBin(-0.01) + 1);*/
+    mixedNorm /= nBinsMixedNorm;
+    mixedNormError /= nBinsMixedNorm;
 
-    // get mixed event normalization by assuming full acceptance at deta of 0
-    Double_t mixedNorm = tracksMixed->Integral(tracksMixed->GetXaxis()->FindBin(-0.01), tracksMixed->GetXaxis()->FindBin(0.01), tracksMixed->GetYaxis()->FindBin(-0.01), tracksMixed->GetYaxis()->FindBin(0.01));
-    mixedNorm /= (tracksMixed->GetXaxis()->FindBin(0.01) - tracksMixed->GetXaxis()->FindBin(-0.01) + 1) * (tracksMixed->GetYaxis()->FindBin(0.01) - tracksMixed->GetYaxis()->FindBin(-0.01) + 1);
+    Float_t triggers = eventMixedAll->Integral(1, eventMixedAll->GetNbinsX(), multBin, multBin);
+//     Printf("%f +- %f | %f | %f", mixedNorm, mixedNormError, triggers, mixedNorm / triggers);
+    
+    mixedNorm /= triggers;
+    mixedNormError /= triggers;
+    
+    delete tracksMixed;
+    
     if (mixedNorm <= 0)
     {
-      Printf("ERROR: Skipping vertex bin %d because mixed event is empty at (0,0)", vertexBin);
+      Printf("ERROR: Skipping multiplicity %d because mixed event is empty at (0,0)", multBin);
+      continue;
     }
-    else
-    {
-      tracksMixed->Scale(1.0 / mixedNorm);
 
-  //     tracksSame->Scale(tracksMixed->Integral() / tracksSame->Integral());
+//     normParameters->Fill(mixedNorm);
       
-      tracksSame->Divide(tracksMixed);
+    TAxis* vertexAxis = trackSameAll->GetAxis(2);
+    for (Int_t vertexBin = 1; vertexBin <= vertexAxis->GetNbins(); vertexBin++)
+    {
+      trackSameAll->GetAxis(2)->SetRange(vertexBin, vertexBin);
+      trackMixedAll->GetAxis(2)->SetRange(vertexBin, vertexBin);
+
+      TH2* tracksSame = trackSameAll->Projection(1, 0, "E");
+      tracksMixed = trackMixedAll->Projection(1, 0, "E");
       
-      // code to draw contributions
-      /*
-      TH1* proj = tracksSame->ProjectionX("projx", tracksSame->GetYaxis()->FindBin(-1.59), tracksSame->GetYaxis()->FindBin(1.59));
-      proj->SetTitle(Form("Bin %d", vertexBin));
-      proj->SetLineColor(vertexBin);
-      proj->DrawCopy((vertexBin > 1) ? "SAME" : "");
-      */
+      // asssume flat in dphi, gain in statistics
+      //     TH1* histMixedproj = mixedTwoD->ProjectionY();
+      //     histMixedproj->Scale(1.0 / mixedTwoD->GetNbinsX());
+      //     
+      //     for (Int_t x=1; x<=mixedTwoD->GetNbinsX(); x++)
+      //       for (Int_t y=1; y<=mixedTwoD->GetNbinsY(); y++)
+      //       mixedTwoD->SetBinContent(x, y, histMixedproj->GetBinContent(y));
+
+      //       delete histMixedproj;
       
-      if (!totalTracks)
-       totalTracks = (TH2*) tracksSame->Clone("totalTracks");
+      Float_t triggers2 = eventMixedAll->Integral(vertexBin, vertexBin, multBin, multBin);
+      if (triggers2 <= 0)
+      {
+       Printf("ERROR: Skipping multiplicity %d vertex bin %d because mixed event is empty", multBin, vertexBin);
+      }
       else
-       totalTracks->Add(tracksSame);
+      {
+       tracksMixed->Scale(1.0 / triggers2 / mixedNorm);
+       // tracksSame->Scale(tracksMixed->Integral() / tracksSame->Integral());
+         
+       tracksSame->Divide(tracksMixed);
+         
+       // code to draw contributions
+       /*
+       TH1* proj = tracksSame->ProjectionX("projx", tracksSame->GetYaxis()->FindBin(-1.59), tracksSame->GetYaxis()->FindBin(1.59));
+       proj->SetTitle(Form("Bin %d", vertexBin));
+       proj->SetLineColor(vertexBin);
+       proj->DrawCopy((vertexBin > 1) ? "SAME" : "");
+       */
+       
+       if (!totalTracks)
+         totalTracks = (TH2*) tracksSame->Clone("totalTracks");
+       else
+         totalTracks->Add(tracksSame);
+
+       totalEvents += eventSameAll->GetBinContent(vertexBin, multBin);
+       
+//     new TCanvas; tracksMixed->DrawCopy("SURF1");
+      }
 
       delete tracksSame;
       delete tracksMixed;
-      
-      totalEvents += eventSameAll->GetBinContent(vertexBin);
     }
   }
 
@@ -853,6 +962,10 @@ TH2* AliUEHist::GetSumOfRatios2(AliUEHist* mixed, AliUEHist::CFStep step, AliUEH
   delete eventSameAll;
   delete eventMixedAll;
   
+//   new TCanvas; normParameters->Draw();
+  
+  TH1::AddDirectory(oldStatus);
+  
   return totalTracks;
 }
 
@@ -2322,3 +2435,25 @@ void AliUEHist::Reset()
   for (Int_t step=0; step<fTrackHistEfficiency->GetNStep(); step++)
     fTrackHistEfficiency->GetGrid(step)->GetGrid()->Reset();
 }
+
+THnBase* AliUEHist::ChangeToThn(THnBase* sparse)
+{
+  // change the object to THn for faster processing
+  
+       // convert to THn (SEGV's for some strange reason...) 
+       // x = THn::CreateHn("a", "a", sparse);
+       
+  // own implementation
+  Int_t nBins[10];
+  for (Int_t i=0; i<sparse->GetNdimensions(); i++)
+    nBins[i] = sparse->GetAxis(i)->GetNbins();
+  THn* tmpTHn = new THnF(Form("%s_thn", sparse->GetName()), sparse->GetTitle(), sparse->GetNdimensions(), nBins, 0, 0);
+  for (Int_t i=0; i<sparse->GetNdimensions(); i++)
+  {
+    tmpTHn->SetBinEdges(i, sparse->GetAxis(i)->GetXbins()->GetArray());
+    tmpTHn->GetAxis(i)->SetTitle(sparse->GetAxis(i)->GetTitle());
+  }
+  tmpTHn->RebinnedAdd(sparse);
+  
+  return tmpTHn;
+}
index ceefd14de53c522fc0845dab53d878d3e4fddb93..0e059edfd7169971c84c91f11a72b0149efe7810 100644 (file)
@@ -21,6 +21,7 @@ class TH2D;
 class TCollection;
 class AliCFGridSparse;
 class THnSparse;
+class THnBase;
 
 class AliUEHist : public TObject
 {
@@ -53,6 +54,8 @@ class AliUEHist : public TObject
   TH2* GetSumOfRatios(AliUEHist* mixed, CFStep step, Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Bool_t etaNorm = kTRUE, Bool_t useVertexBins = kFALSE);
   
   void GetHistsZVtx(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, TH3** trackHist, TH1** eventHist);
+  void GetHistsZVtxMult(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, THnBase** trackHist, TH2** eventHist);
+  
   TH2* GetSumOfRatios2(AliUEHist* mixed, AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd);
 
   TH1* GetTrackEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2 = -1, Int_t source = 1, Int_t axis3 = -1);
@@ -100,7 +103,12 @@ class AliUEHist : public TObject
   void AdditionalDPhiCorrection(Int_t step);
   
   void SetBinLimits(AliCFGridSparse* grid);
+  void SetBinLimits(THnBase* grid);
+
   void ResetBinLimits(AliCFGridSparse* grid);
+  void ResetBinLimits(THnBase* grid);
+  
+  void SetGetMultCache(Bool_t flag = kTRUE) { fGetMultCacheOn = flag; }
   
   AliUEHist(const AliUEHist &c);
   AliUEHist& operator=(const AliUEHist& corr);
@@ -109,6 +117,7 @@ class AliUEHist : public TObject
   virtual Long64_t Merge(TCollection* list);
   void Scale(Double_t factor);
   void Reset();
+  THnBase* ChangeToThn(THnBase* sparse);
   
 protected:
   void SetStepNames(AliCFContainer* container);
@@ -134,9 +143,12 @@ protected:
   
   AliCFContainer* fCache;             //! cache variable for GetTrackEfficiency
   
+  Bool_t fGetMultCacheOn;             //! cache for GetHistsZVtxMult function active
+  THnBase* fGetMultCache;             //! cache for GetHistsZVtxMult function
+  
   TString fHistogramType;             // what is stored in this histogram
   
-  ClassDef(AliUEHist, 8) // underlying event histogram container
+  ClassDef(AliUEHist, 9) // underlying event histogram container
 };
 
 #endif
index 267de5d70cc74efd263cace8eaca91e0e6ff2d39..6d7de7749fb7038aa139546df88cd34e0aa39151 100644 (file)
@@ -432,7 +432,7 @@ void AliUEHistograms::Fill(AliVParticle* leadingMC, AliVParticle* leadingReco)
 }
 
 //____________________________________________________________________
-void AliUEHistograms::FillCorrelations(Double_t centrality, Float_t zVtx, AliUEHist::CFStep step, TObjArray* particles, TObjArray* mixed, Float_t weight, Bool_t firstTime, Bool_t twoTrackEfficiencyCut, Float_t bSign)
+void AliUEHistograms::FillCorrelations(Double_t centrality, Float_t zVtx, AliUEHist::CFStep step, TObjArray* particles, TObjArray* mixed, Float_t weight, Bool_t firstTime, Bool_t twoTrackEfficiencyCut, Float_t bSign, Float_t twoTrackEfficiencyCutValue)
 {
   // fills the fNumberDensityPhi histogram
   //
@@ -561,13 +561,13 @@ void AliUEHistograms::FillCorrelations(Double_t centrality, Float_t zVtx, AliUEH
          Float_t deta = triggerEta - eta[j];
              
          // optimization
-         if (TMath::Abs(deta) < 0.05)
+         if (TMath::Abs(deta) < twoTrackEfficiencyCutValue * 2.5)
          {
            // check first boundaries to see if is worth to loop and find the minimum
            Float_t dphistar1 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 0.8, bSign);
            Float_t dphistar2 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 2.5, bSign);
            
-           const Float_t kLimit = 0.06;
+           const Float_t kLimit = twoTrackEfficiencyCutValue * 3;
 
            Float_t dphistarminabs = 1e5;
            Float_t dphistarmin = 1e5;
@@ -588,7 +588,7 @@ void AliUEHistograms::FillCorrelations(Double_t centrality, Float_t zVtx, AliUEH
              
              fTwoTrackDistancePt[0]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));
              
-             if (dphistarminabs < 0.02 && TMath::Abs(deta) < 0.02)
+             if (dphistarminabs < twoTrackEfficiencyCutValue && TMath::Abs(deta) < twoTrackEfficiencyCutValue)
              {
 //             Printf("Removed track pair %d %d with %f %f %f %f %f %f %f %f %f", i, j, deta, dphistarminabs, phi1, pt1, charge1, phi2, pt2, charge2, bSign);
                continue;
index 3112852a0f299c0ccbc1b98ba2ee2289da796a81..ac2fa581e07d6c2ead3316d798a57929bcade35e 100644 (file)
@@ -28,7 +28,7 @@ class AliUEHistograms : public TNamed
   virtual ~AliUEHistograms();
   
   void Fill(Int_t eventType, Float_t zVtx, AliUEHist::CFStep step, AliVParticle* leading, TList* toward, TList* away, TList* min, TList* max);
-  void FillCorrelations(Double_t centrality, Float_t zVtx, AliUEHist::CFStep step, TObjArray* particles, TObjArray* mixed = 0, Float_t weight = 1, Bool_t firstTime = kTRUE, Bool_t twoTrackEfficiencyCut = kFALSE, Float_t bSign = 0);
+  void FillCorrelations(Double_t centrality, Float_t zVtx, AliUEHist::CFStep step, TObjArray* particles, TObjArray* mixed = 0, Float_t weight = 1, Bool_t firstTime = kTRUE, Bool_t twoTrackEfficiencyCut = kFALSE, Float_t bSign = 0, Float_t twoTrackEfficiencyCutValue = 0.02);
   void Fill(AliVParticle* leadingMC, AliVParticle* leadingReco);
   void FillEvent(Int_t eventType, Int_t step);
   void FillEvent(Double_t centrality, Int_t step);
index 70eb98fad7a1fe930524152c06a63f316799bf79..289dc46448e56cd72b25cdfbf673801104d030fa 100644 (file)
@@ -84,7 +84,7 @@ fFillMixed(kTRUE),
 fMixingTracks(50000),
 fCompareCentralities(kFALSE),
 fTwoTrackEfficiencyStudy(kFALSE),
-fTwoTrackEfficiencyCut(kFALSE),
+fTwoTrackEfficiencyCut(0),
 fUseVtxAxis(kFALSE),
 fSkipTrigger(kFALSE),
 // pointers to UE classes
@@ -244,7 +244,7 @@ void  AliAnalysisTaskPhiCorrelations::CreateOutputObjects()
 
   // event mixing
   Int_t trackDepth = fMixingTracks; 
-  Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemented of AliEventPool
+  Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemented of AliEventPoolManager
    
   Int_t nCentralityBins  = fHistos->GetUEHist(2)->GetEventHist()->GetNBins(1);
   Double_t* centralityBins = (Double_t*) fHistos->GetUEHist(2)->GetEventHist()->GetAxis(1, 0)->GetXbins()->GetArray();
@@ -358,8 +358,9 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseCorrectionMode()
   AliVEvent* inputEvent = fAOD;
   if (!inputEvent)
     inputEvent = fESD;
-  
-  fHistos->SetRunNumber(inputEvent->GetRunNumber());
+
+  if (inputEvent)
+    fHistos->SetRunNumber(inputEvent->GetRunNumber());
 
   TObject* mc = fArrayMC;
   if (!mc)
@@ -661,8 +662,8 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseDataMode()
     
     ((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(1);
     
-    if (fTwoTrackEfficiencyCut)
-      fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracksClone, 0, weight, kTRUE, kTRUE, bSign);
+    if (fTwoTrackEfficiencyCut > 0)
+      fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracksClone, 0, weight, kTRUE, kTRUE, bSign, fTwoTrackEfficiencyCut);
   }
 
   // fill second time with SPD centrality
@@ -719,8 +720,8 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseDataMode()
        if (!fSkipStep6)
          fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracksClone, bgTracks, 1.0 / nMix, (jMix == 0));
 
-       if (fTwoTrackEfficiencyCut)
-         fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracksClone, bgTracks, 1.0 / nMix, (jMix == 0), kTRUE, bSign);
+       if (fTwoTrackEfficiencyCut > 0)
+         fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracksClone, bgTracks, 1.0 / nMix, (jMix == 0), kTRUE, bSign, fTwoTrackEfficiencyCut);
       }
     }
     
index c2078c21c3c581b3dd63d9443135e7366f73bfd8..9f4e46802a375ba5dbac8824f86448d284b58414 100644 (file)
@@ -65,7 +65,7 @@ class  AliAnalysisTaskPhiCorrelations : public AliAnalysisTask
     virtual    void    SetMixingTracks(Int_t tracks) { fMixingTracks = tracks; }
     virtual    void    SetCompareCentralities(Bool_t flag) { fCompareCentralities = flag; }
     virtual    void    SetTwoTrackEfficiencyStudy(Bool_t flag) { fTwoTrackEfficiencyStudy = flag; }
-    virtual    void    SetTwoTrackEfficiencyCut(Bool_t flag) { fTwoTrackEfficiencyCut = flag; }
+    virtual    void    SetTwoTrackEfficiencyCut(Float_t value = 0.02) { fTwoTrackEfficiencyCut = value; }
     virtual    void    SetUseVtxAxis(Bool_t flag) { fUseVtxAxis = flag; }
     virtual     void    SetSkipTrigger(Bool_t flag) { fSkipTrigger = flag; }
     
@@ -111,7 +111,7 @@ class  AliAnalysisTaskPhiCorrelations : public AliAnalysisTask
     Int_t              fMixingTracks;          // size of track buffer for event mixing
     Bool_t             fCompareCentralities;   // use the z vtx axis for a centrality comparison
     Bool_t             fTwoTrackEfficiencyStudy; // two-track efficiency study on
-    Bool_t             fTwoTrackEfficiencyCut;   // enable two-track efficiency cut
+    Float_t            fTwoTrackEfficiencyCut;   // enable two-track efficiency cut
     Bool_t             fUseVtxAxis;              // use z vtx as axis (needs 7 times more memory!)
     Bool_t             fSkipTrigger;             // skip trigger selection
     
@@ -157,7 +157,7 @@ class  AliAnalysisTaskPhiCorrelations : public AliAnalysisTask
     
     Bool_t fFillpT;                // fill sum pT instead of number density
     
-    ClassDef( AliAnalysisTaskPhiCorrelations, 9); // Analysis task for delta phi correlations
+    ClassDef( AliAnalysisTaskPhiCorrelations, 10); // Analysis task for delta phi correlations
   };
 
 class AliDPhiBasicParticle : public AliVParticle
index 7bf3ed24669275d4b7e1140db6229cd19954e854..04b398ed5d34b37128cc84832866d027852cef79 100644 (file)
@@ -6547,24 +6547,24 @@ void ExampleDEtaDPhi(const char* fileNamePbPb, const char* fileNamePbPbMix)
   hist1->Draw("SURF1");
 }
 
-void PlotDeltaPhiEtaGap(const char* fileNamePbPb, const char* fileNamePbPbMix, const char* fileNamepp)
+void PlotDeltaPhiEtaGap(const char* fileNamePbPb, const char* fileNamePbPbMix, const char* fileNamepp, const char* outputFile = "dphi_corr.root")
 {
   loadlibs();
   
   if (!fileNamePbPbMix)
     fileNamePbPbMix = fileNamePbPb;
   
-  file = TFile::Open("dphi_corr.root", "RECREATE");
+  file = TFile::Open(outputFile, "RECREATE");
   file->Close();
   
   Int_t leadingPtOffset = 1;
     
-  Int_t maxLeadingPt = 3;
+  Int_t maxLeadingPt = 4;
   Int_t maxAssocPt = 7;
   if (1)
   {
 //     Float_t leadingPtArr[] = { 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0, 20.0 };
-    Float_t leadingPtArr[] = { 2.0, 4.0, 8.0, 15.0, 20.0 };
+    Float_t leadingPtArr[] = { 2.0, 3.0, 4.0, 8.0, 15.0, 20.0 };
     Float_t assocPtArr[] =     { 0.15, 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 12.0 };
   }
   else
@@ -6580,6 +6580,9 @@ void PlotDeltaPhiEtaGap(const char* fileNamePbPb, const char* fileNamePbPbMix, c
   AliUEHistograms* h2 = (AliUEHistograms*) GetUEHistogram(fileNamepp);
   hMixed2 = (AliUEHistograms*) GetUEHistogram(fileNamepp, 0, kTRUE);
 
+//   h->GetUEHist(2)->SetGetMultCache();
+//   hMixed->GetUEHist(2)->SetGetMultCache();
+  
   if (0)
   {
     h->SetZVtxRange(-0.99, 0.99);
@@ -6589,7 +6592,7 @@ void PlotDeltaPhiEtaGap(const char* fileNamePbPb, const char* fileNamePbPbMix, c
   }
   
   for (Int_t i=0; i<maxLeadingPt; i++)
-    for (Int_t j=0; j<maxAssocPt; j++)
+    for (Int_t j=1; j<maxAssocPt; j++)
     {
       gpTMin = assocPtArr[j] + 0.01;
       gpTMax = assocPtArr[j+1] - 0.01;
@@ -6619,6 +6622,9 @@ void PlotDeltaPhiEtaGap(const char* fileNamePbPb, const char* fileNamePbPbMix, c
       if (1)
       {
        GetSumOfRatios(h, hMixed, &hist1,  step, 0,  10, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, kTRUE); 
+       
+//     new TCanvas; hist1->Draw("SURF1"); return;
+
        GetSumOfRatios(h, hMixed, &hist5,  step, 10,  20, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, kTRUE); 
        GetSumOfRatios(h, hMixed, &hist4,  step, 20,  40, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, kTRUE); 
        GetSumOfRatios(h, hMixed, &hist6,  step, 40,  60, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, kTRUE); 
@@ -6634,11 +6640,11 @@ void PlotDeltaPhiEtaGap(const char* fileNamePbPb, const char* fileNamePbPbMix, c
        GetDistAndFlow(h, hMixed, &hist4,  0, step, 20,  40, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, histType, equivMixedBin, 0, scaleToPairs); 
        GetDistAndFlow(h, hMixed, &hist6,  0, step, 40,  60, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, histType, equivMixedBin, 0, scaleToPairs); 
        GetDistAndFlow(h, hMixed, &hist2,  0, step, 60,  90, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, histType, equivMixedBin, 0, scaleToPairs);
-       step = 6;
+//     step = 6;
        GetDistAndFlow(h2, hMixed2, &hist3,  0, step, 0, -1, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, histType, equivMixedBin, 0, scaleToPairs);
       }
 
-      file = TFile::Open("dphi_corr.root", "UPDATE");
+      file = TFile::Open(outputFile, "UPDATE");
       
       if (hist1)
       {
@@ -10549,6 +10555,8 @@ void PlotQA(const char* fileName)
   
   h->SetPtRange(1.01, 3.99);
   dphi_corr = h->GetUEHist(2)->GetUEHist(AliUEHist::kCFStepReconstructed, 0, 4.01, 14.99, 1, 8);
+  if (dphi_corr->GetEntries() == 0)
+    dphi_corr = h->GetUEHist(2)->GetUEHist(AliUEHist::kCFStepReconstructed+2, 0, 4.01, 14.99, 1, 8);
   
   AliUEHistograms* hMixed = (AliUEHistograms*) GetUEHistogram(fileName, 0, kTRUE);
   if (hMixed->GetUEHist(2)->GetTrackHist(0)->GetGrid(6)->GetGrid()->GetNbins() == 0)
@@ -10559,6 +10567,8 @@ void PlotQA(const char* fileName)
 
   hMixed->SetPtRange(1.01, 3.99);
   dphi_corr_mixed = hMixed->GetUEHist(2)->GetUEHist(AliUEHist::kCFStepReconstructed, 0, 4.01, 14.99, 1, 8);
+  if (dphi_corr_mixed->GetEntries() == 0)
+    dphi_corr_mixed = hMixed->GetUEHist(2)->GetUEHist(AliUEHist::kCFStepReconstructed+2, 0, 4.01, 14.99, 1, 8);
   
   if (runNumber != 0 && runNumber != h->GetRunNumber())
     AliFatal("Run numbers inconsistent");
index 06c11420fdf455a5225b509ecc9f433a2af2e339..d7144027c5e014be74918e048816ca5fd4472ea9 100644 (file)
@@ -117,12 +117,17 @@ Double_t DeltaPhiWidth2DFitFunction(Double_t *x, Double_t *par)
 void SubtractEtaGap(TH2* hist, Float_t etaLimit, Float_t outerLimit, Bool_t scale, Bool_t drawEtaGapDist = kFALSE)
 {
   TString histName(hist->GetName());
+  Int_t etaBins = 0;
 
   TH1D* etaGap = hist->ProjectionX(histName + "_1", TMath::Max(1, hist->GetYaxis()->FindBin(-outerLimit + 0.01)), hist->GetYaxis()->FindBin(-etaLimit - 0.01));
-  Int_t etaBins = hist->GetYaxis()->FindBin(-etaLimit - 0.01) - TMath::Max(1, hist->GetYaxis()->FindBin(-outerLimit + 0.01)) + 1;
+  Printf("%f", etaGap->GetEntries());
+  if (etaGap->GetEntries() > 0)
+    etaBins += hist->GetYaxis()->FindBin(-etaLimit - 0.01) - TMath::Max(1, hist->GetYaxis()->FindBin(-outerLimit + 0.01)) + 1;
 
   TH1D* tracksTmp = hist->ProjectionX(histName + "_2", hist->GetYaxis()->FindBin(etaLimit + 0.01), TMath::Min(hist->GetYaxis()->GetNbins(), hist->GetYaxis()->FindBin(outerLimit - 0.01)));
-  etaBins += TMath::Min(hist->GetYaxis()->GetNbins(), hist->GetYaxis()->FindBin(outerLimit - 0.01)) - hist->GetYaxis()->FindBin(etaLimit + 0.01) + 1;
+  Printf("%f", tracksTmp->GetEntries());
+  if (tracksTmp->GetEntries() > 0)
+    etaBins += TMath::Min(hist->GetYaxis()->GetNbins(), hist->GetYaxis()->FindBin(outerLimit - 0.01)) - hist->GetYaxis()->FindBin(etaLimit + 0.01) + 1;
   
   etaGap->Add(tracksTmp);
 
@@ -1497,12 +1502,15 @@ void DrawExample(const char* histFileName, Int_t i, Int_t j, Int_t histId, TH1**
   c->SaveAs(Form("ex/%s.png", c->GetName()));
   c->SaveAs(Form("ex/%s.eps", c->GetName()));
 
-  c = new TCanvas(Form("ex_%d_%d_%d_a", i, j, histId), Form("ex_%d_%d_%d_a", i, j, histId), 400, 400);
-  hist->SetTitle("");
-  hist->DrawCopy("SURF1");
-  paveText->Draw();
-  c->SaveAs(Form("ex/%s.png", c->GetName()));
-  c->SaveAs(Form("ex/%s.eps", c->GetName()));
+  if (0)
+  {
+    c = new TCanvas(Form("ex_%d_%d_%d_a", i, j, histId), Form("ex_%d_%d_%d_a", i, j, histId), 400, 400);
+    hist->SetTitle("");
+    hist->DrawCopy("SURF1");
+    paveText->Draw();
+    c->SaveAs(Form("ex/%s.png", c->GetName()));
+    c->SaveAs(Form("ex/%s.eps", c->GetName()));
+  }
 }
 
 void DrawExampleAll(const char* histFileName)
@@ -1590,6 +1598,39 @@ void DrawExampleAll(const char* histFileName)
   }
 }
 
+void DrawDoubleHump(const char* histFileName)
+{
+  Float_t exampleI[] = { 0, 0, 0, 1, 1, 1};
+  Float_t exampleJ[] = { 0, 1, 2, 0, 1, 2};
+  
+  TH1* projectionsPhi[9];
+  TH1* projectionsEta[9];
+  
+  TCanvas* c = new TCanvas("DrawDoubleHump", "DrawDoubleHump", 1200, 800);
+  c->Divide(3, 2);
+  
+  TLegend* legend = new TLegend(0.15, 0.7, 0.88, 0.88);
+  legend->SetFillColor(0);
+  
+  for (Int_t i=0; i<6; i++)
+  {
+    DrawExample(histFileName, exampleI[i], exampleJ[i], 0, &projectionsPhi[i], &projectionsEta[i]);
+
+    c->cd(i+1);
+    TH1* clone = projectionsPhi[i]->DrawCopy("");
+    clone->SetLineColor(1);
+    clone->GetYaxis()->SetRangeUser(clone->GetMinimum(), clone->GetMaximum() * 1.2);
+//     clone->GetXaxis()->SetTitle(Form("%s / %s", clone->GetXaxis()->GetTitle(), "#Delta#eta"));
+    clone->GetXaxis()->SetTitle(Form("%s / %s", "#Delta#phi (rad.)", "#Delta#eta"));
+    
+    clone = projectionsEta[i]->DrawCopy("SAME");
+    clone->SetLineColor(2);
+
+    DrawLatex(0.3, 0.85, 1, Form("#Delta#phi projection in |#Delta#eta| < %.1f", 0.8), 0.04);
+    DrawLatex(0.3, 0.80, 2, Form("#Delta#eta projection in |#Delta#phi| < %.1f", 0.8), 0.04);
+  }
+}
+
 void DrawFullCentralityDependence(const char* histFileName)
 {
   Int_t colors[] = { 1, 2, 4, 3, 5, 6 };
@@ -1857,7 +1898,7 @@ void TestTwoGaussian()
 void DrawEtaGapExample(const char* histFileName, Int_t i = 1, Int_t j = 2)
 {
   Float_t etaLimit = 1.0;
-  Float_t outerLimit = 1.79;
+  Float_t outerLimit = 1.59;
   Float_t projectLimit = 0.8;
 
   TFile::Open(histFileName);