Initial autocorr checkin.
authorloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 14 Sep 2010 15:45:28 +0000 (15:45 +0000)
committerloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 14 Sep 2010 15:45:28 +0000 (15:45 +0000)
PWG4/TwoPartCorr/AutoCorr.C [new file with mode: 0644]
PWG4/TwoPartCorr/AutoCorr.h [new file with mode: 0644]
PWG4/TwoPartCorr/EventPool.C [new file with mode: 0644]
PWG4/TwoPartCorr/EventPool.h [new file with mode: 0644]
PWG4/TwoPartCorr/TreeClasses.C [new file with mode: 0644]
PWG4/TwoPartCorr/TreeClasses.h [moved from PWG4/TwoPartCorr/MyTreeClasses.C with 98% similarity]
PWG4/TwoPartCorr/rootlogon.C
PWG4/TwoPartCorr/runAutoCorr.C [new file with mode: 0644]

diff --git a/PWG4/TwoPartCorr/AutoCorr.C b/PWG4/TwoPartCorr/AutoCorr.C
new file mode 100644 (file)
index 0000000..a676992
--- /dev/null
@@ -0,0 +1,142 @@
+// $Id:$
+
+#include <Riostream.h>
+#include <TTree.h>
+#include <TChain.h>
+#include <TCanvas.h>
+#include <TClonesArray.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include "Rtypes.h"
+#include "TMath.h"
+#include "TFile.h"
+#include "TSystem.h"
+#include "TROOT.h"
+#include <vector>
+#include <deque>
+#include "TRandom.h"
+#include "TreeClasses.h"
+#include "AutoCorr.h"
+
+Int_t AutoCorr::InitEventPools(Int_t depth, 
+                              Int_t nMultBins, Double_t multbin[], 
+                              Int_t nZvtxBins, Double_t zvtxbin[])
+{
+  // First assign AutoCorr members
+  fNMultBins = nMultBins;
+  fNZvtxBins = nZvtxBins;
+
+  for (int iM=0; iM<nMultBins; iM++) {
+    std::vector<EventPool*> evp;
+    for (int iZ=0; iZ<nZvtxBins; iZ++) {
+      evp.push_back(new EventPool(depth, 
+                                 multbin[iM], multbin[iM+1], 
+                                 zvtxbin[iZ], zvtxbin[iZ+1] ));
+    }
+    fEvPool.push_back(evp);
+  }
+  
+  for (int iM=0; iM<nMultBins; iM++) {
+    for (int iZ=0; iZ<nZvtxBins; iZ++) {
+      fEvPool.at(iM).at(iZ)->SetMultBinIndex(iM);
+      fEvPool.at(iM).at(iZ)->SetZvtxBinIndex(iZ);
+    }
+  }
+  
+  bool print_this = false;
+  if (print_this) {
+    cout << "fEvPool outer size: " << fEvPool.size() << endl;
+    for (int iM=0; iM<nMultBins; iM++) {
+      for (int iZ=0; iZ<nZvtxBins; iZ++) {
+       if(fEvPool.at(iM).at(iZ)) {
+         cout << "multiplicity bin: " << iM;
+         cout << ", z-vertex bin: " << iZ;
+         fEvPool.at(iM).at(iZ)->PrintInfo();
+       }
+      }
+    }
+  }
+  
+  return 0;
+}
+
+EventPool* AutoCorr::GetEventPool(int iMult, int iZvtx)
+{
+  if (iMult < 0 || iMult >= fNMultBins) return 0x0;
+  if (iZvtx < 0 || iZvtx >= fNZvtxBins) return 0x0;
+  return fEvPool.at(iMult).at(iZvtx);
+}
+
+Int_t AutoCorr::UpdatePools(int iEvent, MyHeader* ev, TClonesArray* trk)
+{
+  for (int iM=0; iM<fNMultBins; iM++) {
+    for (int iZ=0; iZ<fNZvtxBins; iZ++) {
+      fEvPool.at(iM).at(iZ)->UpdatePool(iEvent, ev, trk);
+    }
+  }  
+  return 0;
+}
+
+
+Double_t AutoCorr::DeltaPhi(MyPart* t1, MyPart* t2,
+                           double rangeMin, double rangeMax)
+{
+  Double_t dphi = -999;
+  Double_t pi = TMath::Pi();
+  
+  if (!t1 || !t2) return -99.;
+  double phia = t1->Phi();  
+  double phib = t2->Phi();  
+  
+  if (phia < 0)         phia += 2*pi;
+  else if (phia > 2*pi) phia -= 2*pi;
+  if (phib < 0)         phib += 2*pi;
+  else if (phib > 2*pi) phib -= 2*pi;
+  dphi = phib - phia;
+  if (dphi < rangeMin)      dphi += 2*pi;
+  else if (dphi > rangeMax) dphi -= 2*pi;
+  
+  return dphi;
+}
+
+Double_t AutoCorr::DeltaEta(MyPart* t1, MyPart* t2)
+{
+  if (!t1 || !t2) return -99.;
+  return t1->Eta() - t2->Eta();
+}
+
+Bool_t AutoCorr::IsTrackOk(MyPart* t)
+{
+  if (!t) return false;
+  return fabs(t->Eta()) <= 1.2;    
+}
+
+Bool_t AutoCorr::IsPairOk(MyPart* t1, MyPart* t2)
+{
+  if (!t1 || !t2) return false;
+  if (!IsTrackOk(t1) || !IsTrackOk(t2)) return false;
+
+  double deta = DeltaEta(t1, t2);
+  double dphi = DeltaPhi(t1, t2);
+  double dpmax = 0.03;
+  double demax = 0.01;
+
+  double dr = dphi*dphi/(dpmax*dpmax) + deta*deta/(demax*demax);
+  return (dr > 1);
+}
+
+// So far same as IsPairOk()
+Bool_t AutoCorr::IsMixedPairOk(MyPart* t1, MyPart* t2)
+{
+  if (!t1 || !t2) return false;
+  if (!IsTrackOk(t1) || !IsTrackOk(t2)) return false;
+
+  double deta = DeltaEta(t1, t2);
+  double dphi = DeltaPhi(t1, t2);
+  double dpmax = 0.03;
+  double demax = 0.01;
+
+  double dr = dphi*dphi/(dpmax*dpmax) + deta*deta/(demax*demax);
+  return (dr > 1);
+}
+
diff --git a/PWG4/TwoPartCorr/AutoCorr.h b/PWG4/TwoPartCorr/AutoCorr.h
new file mode 100644 (file)
index 0000000..d0308e3
--- /dev/null
@@ -0,0 +1,54 @@
+// $Id:$
+
+#ifndef AutoCorr_h
+#define AutoCorr_h
+
+#include <Riostream.h>
+#include <TTree.h>
+#include <TChain.h>
+#include <TCanvas.h>
+#include <TClonesArray.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include "Rtypes.h"
+#include "TMath.h"
+#include "TFile.h"
+#include "TSystem.h"
+#include "TROOT.h"
+#include <vector>
+#include <deque>
+#include "TRandom.h"
+#include "TreeClasses.h"
+#include "EventPool.h"
+
+class AutoCorr : public TObject
+{
+public:
+  AutoCorr(){;}
+  ~AutoCorr(){;}
+  
+  // RangeMin,Max specify periodic boundaries
+  Double_t DeltaPhi(MyPart* t1, MyPart* t2, 
+                   double rangeMin = -TMath::Pi()/2, 
+                   double rangeMax = 3*TMath::Pi()/2);
+  Double_t DeltaEta(MyPart* t1, MyPart* t2);  
+  Bool_t IsTrackOk(MyPart* t);
+  Bool_t IsPairOk(MyPart* t1, MyPart* t2);
+  Bool_t IsMixedPairOk(MyPart* t1, MyPart* t2);
+//  Bool_t IsEventOk(MyHeader *h);
+
+  Int_t InitEventPools(Int_t depth, 
+                      Int_t nmultbins, Double_t multbins[], 
+                      Int_t nzvtxbins, Double_t zvtxbins[]);
+
+  Int_t UpdatePools(int iEvent, MyHeader* ev, TClonesArray* trk);
+  EventPool* GetEventPool(int iMult, int iZvtx);// { return fEvPool.at(iMult).at(iZvtx); }
+protected:
+  Int_t fNMultBins;
+  Int_t fNZvtxBins;
+
+  std::vector<std::vector<EventPool*> > fEvPool; // [fMultBin][fZvtxBin]
+
+  ClassDef(AutoCorr,1)
+};
+#endif
diff --git a/PWG4/TwoPartCorr/EventPool.C b/PWG4/TwoPartCorr/EventPool.C
new file mode 100644 (file)
index 0000000..20fc2f9
--- /dev/null
@@ -0,0 +1,185 @@
+// $Id:$
+
+#include "EventPool.h"
+
+bool debug = true;
+
+void EventPool::PrintInfo()
+{
+  cout << " --- --- --- " << endl;
+  cout << Form("%20s: %d", "Pool capacity", fMixDepth) << endl;
+  cout << Form("%20s: %d", "Current depth", Depth()) << endl;
+  cout << Form("%20s: %d-%d", "Mult. range", fMultMin, fMultMax) << endl;
+  cout << Form("%20s: %.1f-%.1f", "Z-vtx range", fZvtxMin, fZvtxMax) << endl;
+
+  return;
+}
+
+Bool_t EventPool::IsPoolReady()
+{
+  return (Int_t)fEvents.size()==fMixDepth;
+}
+
+Bool_t EventPool::EventMatchesBin(Int_t mult, Short_t zvtx)
+{
+  // N.B. Lower bin limit included; upper limit excluded.
+  Bool_t multOK = (mult >= fMultMin && mult < fMultMax);
+  Bool_t zvtxOK = (zvtx >= fZvtxMin && zvtx < fZvtxMax);
+  return (multOK && zvtxOK);
+}
+
+Int_t EventPool::Depth()
+{
+  return fEvents.size();
+}
+Int_t EventPool::TracksInPool()
+{
+  int ntrk=0;
+  for (int i=0; i<(int)fEvents.size(); i++) {
+    ntrk += fNTracksInEvent.at(i);
+  }
+  return ntrk;
+}
+
+Int_t EventPool::SetEventMultRange(int multMin, int multMax)
+{
+  fMultMin = multMin;
+  fMultMax = multMax;
+  return 0;
+}
+Int_t EventPool::SetEventZvtxRange(int zvtxMin, int zvtxMax)
+{
+  fZvtxMin = zvtxMin;
+  fZvtxMax = zvtxMax;
+  return 0;
+}
+
+Int_t EventPool::UpdatePool(int iEvent, MyHeader* ev, TClonesArray* trk)
+{
+  // Initialize at any chosen starting event
+  if (!fTracks) fTracks = new TClonesArray("MyPart", 1000);
+
+  fMult = trk->GetEntries();
+  fZvtx = ev->fVz; // Short_t. TODO: check--should fVc be used instead?
+
+  if (!EventMatchesBin(fMult, fZvtx)) {
+    fWasUpdated = false;
+    return -1;
+  }
+
+  // Should see evsize = trsize (= fMixDepth once full).
+  Int_t evsize = fEvents.size();
+  Int_t ntsize = fNTracksInEvent.size();
+
+  if (evsize != ntsize) 
+    cout << "WARNING:  Event array and track counter array sizes do not match."
+        << " evsize = " << evsize
+        << " ntsize = " << ntsize
+        << endl;
+
+  Bool_t firstReachedCapacity = false;
+  if (evsize == fMixDepth - 1) firstReachedCapacity = true;
+
+  // Remove 0th element before appending this event
+  if (evsize >= fMixDepth) {
+    // TODO: find out - does popping delete the elements from memory?
+    fNTracksInEvent.pop_front(); // remove first int
+    fEvents.pop_front();        // remove first track array 
+    fEventIndex.pop_front();
+  }
+
+  // TODO: Clone() is expensive. Maybe a better way available?
+  fTracks = (TClonesArray*) trk->Clone();
+
+  fNTracksInEvent.push_back(ev->fNSelTracks);
+  fEvents.push_back(fTracks);
+  fEventIndex.push_back(iEvent);
+
+  if (firstReachedCapacity) {
+    cout << "Pool " << MultBinIndex() << ", " << ZvtxBinIndex() 
+        << " ready at event "<< iEvent;
+    PrintInfo();
+  }
+
+  fWasUpdated = true;
+
+  bool print = false;
+  if (debug && print) {
+    cout << " Event " << fEventIndex.back();
+    cout << " PoolDepth = " << Depth();
+    cout << " NTracks = " << NTracksInCurrentEvent();
+    // for (int i=1; i<Depth(); i++) {
+    //   cout << " " << NTracksInEvent(iEvent-i);
+    // }
+    cout << " TracksInPool = " << TracksInPool();
+    //    cout << endl;
+  }
+
+  return 0;
+}
+
+MyPart* EventPool::GetRandomTrack()
+{
+  MyPart* trk = 0;
+  TClonesArray* tca = 0;
+  
+ // index of random event in the pool
+  UInt_t ranEvt = gRandom->Integer(fEvents.size()-1);
+  tca = fEvents.at(ranEvt);
+
+ // index of random track in event
+  UInt_t ranTrk = gRandom->Integer(tca->GetEntries());
+  
+  trk = (MyPart*)tca->At(ranTrk);
+  
+  /*
+  if (debug) {
+    // loop over the event buffer
+    for (int i=0; i<(int)fEvents.size(); i++) {
+      TClonesArray* tca = fEvents.at(i);
+      cout << tca->GetEntries() << " ";
+    }
+    cout << endl;
+    // Get first track inside each
+    for (int i=0; i<(int)fEvents.size(); i++) {
+      TClonesArray* tca = fEvents.at(i);
+      trk = (MyPart*)tca->At(0);
+      cout << trk->Pt() << " ";
+    }
+    cout << endl;
+  }
+  */
+
+  return trk;
+}
+
+Int_t EventPool::NTracksInCurrentEvent()
+{
+  return fNTracksInEvent.back();
+}
+
+// Important!! This fn. will break if selective filling is implemented
+// (i.e. pool is not updated for every single event)
+// TODO: Implement internal counters to fix this.
+Int_t EventPool::NTracksInEvent(int iEvent)
+{
+  Int_t n = -1;
+  Int_t curEvent = fEventIndex.back();
+  Int_t offset = curEvent - iEvent;
+  // pos = position of iEvent in rolling buffer
+  int pos = fEventIndex.size() - offset - 1;
+
+  if (offset==0) // (iEvent == curEvent)
+    n = fNTracksInEvent.back();
+  else if (offset < 0 || iEvent < 0) {
+    n = 0;//    cout << " No info for event " << iEvent << " ";
+  }
+  else if (offset > 0 && offset <= (int)fEventIndex.size()) {
+    n = fNTracksInEvent.at(pos);
+  }
+  else
+    cout << "Event info no longer in memory" << endl;
+  return n;
+}
+
+
diff --git a/PWG4/TwoPartCorr/EventPool.h b/PWG4/TwoPartCorr/EventPool.h
new file mode 100644 (file)
index 0000000..a887844
--- /dev/null
@@ -0,0 +1,71 @@
+// $Id:$
+
+#ifndef EventPool_h
+#define EventPool_h
+
+#include <vector>
+#include <deque>
+#include "Rtypes.h"
+#include <Riostream.h>
+#include <TClonesArray.h>
+#include "TFile.h"
+#include "TMath.h"
+#include "TRandom.h"
+#include "TSystem.h"
+#include "TreeClasses.h"
+
+class EventPool : public TObject
+{
+public:
+  EventPool(Int_t d) : fMixDepth(d), fMultMin(-999), fMultMax(-999), 
+    fZvtxMin(-999), fZvtxMax(-999) {;}
+  EventPool(Int_t d, Int_t multMin, Int_t multMax, 
+           Double_t zvtxMin, Double_t zvtxMax) : 
+    fMixDepth(d), fMultMin(multMin), fMultMax(multMax), 
+    fZvtxMin(zvtxMin), fZvtxMax(zvtxMax) {;}
+  ~EventPool(){;}
+
+  // TODO: implement cleanup
+  // Int_t Drain();
+  Int_t SetEventMultRange(int multMin, int multMax);
+  Int_t SetEventZvtxRange(int zvtxMin, int zvtxMax);
+  
+  Int_t UpdatePool(int iEvent, MyHeader* ev, TClonesArray* trk);
+  Bool_t WasUpdated() { return fWasUpdated; }
+  Bool_t EventMatchesBin(Int_t mult, Short_t zvtx);
+  Bool_t IsPoolReady();              // Contains fMixDepth events? 
+  Int_t Depth();                     // Number of events in pool
+  Int_t TracksInPool();              // Total tracks in pool
+  Int_t NTracksInEvent(int iEvent);
+  Int_t NTracksInCurrentEvent();
+  MyPart* GetRandomTrack();
+  void PrintInfo();
+  Int_t MultBinIndex() { return fMultBinIndex; }
+  Int_t ZvtxBinIndex() { return fZvtxBinIndex; }
+  void SetMultBinIndex(Int_t iM) { fMultBinIndex = iM; }
+  void SetZvtxBinIndex(Int_t iZ) { fZvtxBinIndex = iZ; }
+
+protected:
+  TClonesArray* fTracks; // Copy of trk array. Refreshes each event.
+
+  //  Use as ring buffers
+  deque<TClonesArray*> fEvents; // The guts of the class. Holds TObjArrays of MyParts.
+  deque<int> fNTracksInEvent;
+  deque<int> fEventIndex;
+
+  Int_t fMixDepth;             // Number of evts. to mix with
+  Int_t fMultMin, fMultMax;    // Track multiplicity bin range
+  Double_t fZvtxMin, fZvtxMax; // Event z-vertex bin range
+  Int_t fMult;                 // Tracks in current event
+  Short_t fZvtx;               // Current z-vertex
+  Bool_t fWasUpdated;          // Evt. succesfully passed selection?
+
+  Int_t fMultBinIndex;
+  Int_t fZvtxBinIndex;
+
+  ClassDef(EventPool,1)
+};
+#endif
+
+
+
diff --git a/PWG4/TwoPartCorr/TreeClasses.C b/PWG4/TwoPartCorr/TreeClasses.C
new file mode 100644 (file)
index 0000000..1e0ceb1
--- /dev/null
@@ -0,0 +1,3 @@
+// $Id:$
+
+#include "TreeClasses.h"
similarity index 98%
rename from PWG4/TwoPartCorr/MyTreeClasses.C
rename to PWG4/TwoPartCorr/TreeClasses.h
index aecead0..5fb8e53 100644 (file)
@@ -1,7 +1,7 @@
 // $Id:$
 
-#ifndef MyTreeClasses_cxx
-#define MyTreeClasses_cxx
+#ifndef TreeClasses_h
+#define TreeClasses_h
 
 #include <Riostream.h>
 #include <TTree.h>
index c87b8ff..401502a 100644 (file)
@@ -1,21 +1,14 @@
 {
   gSystem->Load("libTree.so");
-  gSystem->Load("libVMC.so");
-  gSystem->Load("libSTEERBase.so");
-  gSystem->Load("libESD.so");
-  gSystem->Load("libAOD.so"); 
-  gSystem->Load("libANALYSIS");
-  gSystem->Load("libANALYSISalice");
-  if (gSystem->Getenv("ALICE_ROOT_ORIG")) 
-    gSystem->AddIncludePath("-I$ALICE_ROOT_ORIG/include");
-  gSystem->AddIncludePath("-I$ALICE_ROOT/include");
   if (gSystem->Getenv("TMPDIR")) 
     gSystem->SetBuildDir(gSystem->Getenv("TMPDIR"));
   if (1) {
       delete gRandom;
       gRandom = new TRandom3(0);
    }
-  gROOT->LoadMacro("MyTreeClasses.C+g");
+  gROOT->LoadMacro("TreeClasses.C+g");
+  gROOT->LoadMacro("EventPool.C+g");
+  gROOT->LoadMacro("AutoCorr.C+g");
 }
 
 
diff --git a/PWG4/TwoPartCorr/runAutoCorr.C b/PWG4/TwoPartCorr/runAutoCorr.C
new file mode 100644 (file)
index 0000000..c3b70a9
--- /dev/null
@@ -0,0 +1,201 @@
+// $Id:$
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include "TreeClasses.h"
+#include "EventPool.h"
+#include "AutoCorr.h"
+#endif
+
+bool debug = true;
+
+void runAutoCorr()
+{
+  const double PI = TMath::Pi();
+
+  //TFile* f = new TFile("output.root", "read");
+  //TList* l = (TList*)f->Get("output");
+  //TTree* tree = (TTree*)l->At(11); 
+  TChain *tree = new TChain("MyTree");
+  tree->Add("res_LHC10e_09132010/mergedruns/merged_*.root");
+  cout << "Entries " << tree->GetEntries() << endl;
+
+  MyHeader* ev = 0;//new MyHeader();
+  TClonesArray* trk = 0;//new TClonesArray("MyPart");
+  AutoCorr* ac = new AutoCorr();
+
+  TBranch* evBranch = tree->GetBranch("header");
+  evBranch->SetAddress(&ev);
+  
+  TBranch* trBranch = tree->GetBranch("parts");
+  trBranch->SetAddress(&trk);
+
+  ////////////////////////
+  // TEST DRIVE
+  ////////////////////////
+  Int_t poolsize = 20;
+  const Int_t nMultBins = 3;
+  Double_t multbin[nMultBins+1] = {180,230, 280, 500}; 
+  const Int_t nZvtxBins = 3;
+  Double_t zvtxbin[nZvtxBins+1] = {-6, -2, 2, 6};
+  ac->InitEventPools(poolsize, nMultBins, multbin, nZvtxBins, zvtxbin);
+
+  // TODO: encapsulate these in AutoCorr (pass arrays above to autocorr ctor)
+  TH1F* hMultBins = new TH1F("hMultBins", "Event multiplicity binning", 
+                            nMultBins, multbin);
+  TH1F* hZvtxBins = new TH1F("hZvtxBins", "Event Z-vertex binning", 
+                            nZvtxBins, zvtxbin);
+
+  if (debug) {
+    cout << "Mult. binning: ";
+  for (int j=1; j<=hMultBins->GetNbinsX(); j++) 
+    cout << hMultBins->GetBinLowEdge(j) << " ";
+  cout << hMultBins->GetBinLowEdge(nMultBins) + 
+    hMultBins->GetBinWidth(nMultBins) << " ";
+  cout << endl;
+    cout << "Z-vtx. binning: ";
+  for (int j=1; j<=hZvtxBins->GetNbinsX(); j++) 
+    cout << hZvtxBins->GetBinLowEdge(j) << " ";
+  cout << hZvtxBins->GetBinLowEdge(nZvtxBins) + 
+    hZvtxBins->GetBinWidth(nZvtxBins) << " ";
+  cout << endl;
+  }
+
+  // More test examples  
+  Short_t zmin = -5, zmax = 5;
+  Int_t multmin = 60, multmax = 100;
+  EventPool* testPool = new EventPool(poolsize); 
+  EventPool* binnedPool = new EventPool(poolsize, multmin, multmax, zmin, zmax); 
+
+
+  // TODO: Put histos in a manager
+  // ---------------------------------
+  Double_t etaMax = 2.2;
+
+  TH2F* hSig[nMultBins];
+  TH2F* hBkg[nMultBins];
+  TH2F* hSB[nMultBins];
+  TH1F* h1 = new TH1F("h1", "h1", 1000, -2, 2);
+  TH1F* h2 = new TH1F("h2", "h2", 1000, -2, 2);
+  TH1F* h3 = new TH1F("h3", "h3", 100, -PI/2, 3*PI/2);
+  
+  for (int iM=0; iM<nMultBins; iM++) {
+    char* name;
+    name = Form("hSig_%i", iM);
+    hSig[iM] = new TH2F(name, name, 100, -PI/2, 3*PI/2, 100, -etaMax, etaMax);
+    name = Form("hBkg_%i", iM);
+    hBkg[iM] = new TH2F(name, name, 100, -PI/2, 3*PI/2, 100, -etaMax, etaMax);
+    name = Form("hSB_%i", iM);
+    hSB[iM] = new TH2F(name, name, 100, -PI/2, 3*PI/2, 100, -etaMax, etaMax);
+  }
+  
+  TH2F* hDR = new TH2F("hDR", "hDR", 100, -0.01, 0.01, 100, -0.005, 0.005);
+  hDR->SetTitle("#Delta#phi-#Delta#eta Distribution (no cut);#Delta#phi;#Delta#eta");
+  TH2F* hDRcut = new TH2F("hDRcut", "hDRcut", 200, -0.2, 0.2, 200, -0.1, 0.1);
+  hDRcut->SetTitle("#Delta#phi-#Delta#eta Distribution;#Delta#phi;#Delta#eta");
+  h2->SetLineColor(kRed);
+  h3->SetLineColor(kBlue);
+  // ---------------------------------
+
+  Int_t nEvents = tree->GetEntries();
+  //nEvents = 10000;
+  for (int n=0; n<nEvents; n++) {
+
+    // GetEntry fills "header" and "parts" branches, reinitializing ev
+    // and trk instances for each event.
+    //tree->GetEntry(n); 
+    evBranch->GetEntry(n);
+    trBranch->GetEntry(n);
+
+//    trk->Print("ls");
+    if (n % 100 == 0) cout << n << endl;
+
+    // tests...don't need
+    testPool->UpdatePool(n, ev, trk);
+    binnedPool->UpdatePool(n, ev, trk);
+
+    ac->UpdatePools(n, ev, trk);
+
+    int SignalEventMultBin = hMultBins->FindBin(ev->fNSelTracks) - 1;
+    int SignalEventZvtxBin = hZvtxBins->FindBin(ev->fVz) - 1;
+    EventPool* pool = ac->GetEventPool(SignalEventMultBin,
+                                      SignalEventZvtxBin);
+
+    if (!pool) continue;     // No mixing pool in range specified by arrays
+    if (!pool->IsPoolReady()) continue; // TODO: Figure out something better
+
+    for (int i=0; i<ev->fNSelTracks-1; i++) {
+      MyPart* t1 = (MyPart*)trk->At(i);
+      if (!ac->IsTrackOk(t1)) continue;
+      
+      h1->Fill(t1->Eta());
+
+      // --- Same-event correlation loop ---
+      for (int j=i+1; j<ev->fNSelTracks; j++) {
+       MyPart* t2 = (MyPart*)trk->At(j);
+       Double_t deta = ac->DeltaEta(t1, t2);
+       Double_t dphi = ac->DeltaPhi(t1, t2);
+       hDR->Fill(dphi, deta); // No cut - contamination pk.
+
+       if (ac->IsPairOk(t1, t2)) {
+         h2->Fill(deta);
+         h3->Fill(dphi);
+         hSig[SignalEventMultBin]->Fill(dphi, deta);
+         hDRcut->Fill(dphi, deta); // Show elliptical hole from cut 
+       }
+      }
+      // ---
+
+      // --- Event mixing loop---
+      // TODO: change 100 to f*<mult> of bkg. event where f is some fixed fraction
+      for (int j=i+1; j<0.2*ev->fNSelTracks; j++) { // TODO: update track selection
+       MyPart* tbg = pool->GetRandomTrack();
+       if (ac->IsMixedPairOk(t1, tbg)) {
+         Double_t deta_bg = ac->DeltaEta(t1, tbg);
+         Double_t dphi_bg = ac->DeltaPhi(t1, tbg);
+         hBkg[SignalEventMultBin]->Fill(dphi_bg, deta_bg);
+       }
+      }
+      // ---      
+
+    }
+    
+  }
+  TCanvas* c1 = new TCanvas("c1", "c1", 1);
+  c1->cd();
+  h2->Draw();
+  h1->Draw("same");
+
+  TCanvas* c2 = new TCanvas("c2", "c2", 1);
+  c2->cd();
+  h3->Draw();
+
+  TCanvas* c4 = new TCanvas("c4", "c4", 1);
+  c4->Divide(2, 1, 0.001, 0.001);
+  c4->cd(1);
+  hDR->Draw("surf2");
+  c4->cd(2);
+  hDRcut->Draw("colz");
+
+  TCanvas* cSig = new TCanvas("cSig", "cSig", 1);
+  cSig->Divide(nMultBins, 1, 0.001, 0.001);
+  for (int iM=0; iM<nMultBins; iM++) {
+    cSig->cd(iM+1);
+    hSig[iM]->Draw("surf1");
+  }
+  TCanvas* cBkg = new TCanvas("cBkg", "cBkg", 1);
+  cBkg->Divide(nMultBins, 1, 0.001, 0.001);
+  for (int iM=0; iM<nMultBins; iM++) {
+    cBkg->cd(iM+1);
+    hBkg[iM]->Draw("surf1");
+  }
+  
+  TCanvas* cSB = new TCanvas("cSB", "cSB", 1);
+  cSB->Divide(nMultBins, 1, 0.001, 0.001);
+  for (int iM=0; iM<nMultBins; iM++) {
+    cSB->cd(iM+1);
+    hSB[iM]->Divide(hSig[iM], hBkg[iM]);
+    hSB[iM]->Draw("surf1");
+  }  
+  return;
+}
+