speed up FemtoESE analysis task
authoraohlson <alice.ohlson@cern.ch>
Mon, 29 Sep 2014 09:25:35 +0000 (11:25 +0200)
committerdgangadh <dhevan.raja.gangadharan@cern.ch>
Mon, 29 Sep 2014 09:40:25 +0000 (11:40 +0200)
PWGCF/FEMTOSCOPY/ESE/AliAnalysisTaskFemtoESE.cxx
PWGCF/FEMTOSCOPY/ESE/AliAnalysisTaskFemtoESE.h

index b2ed386..f841f2a 100644 (file)
@@ -21,6 +21,8 @@
 #include "TObjectTable.h"
 #include <vector>
 
+//#include "TStopwatch.h"
+
 #include "AliAnalysisTask.h"
 #include "AliAnalysisManager.h"
 
@@ -51,48 +53,93 @@ ClassImp(AliAnalysisTaskFemtoESE)
 //________________________________________________________________________
 // Default constructor
 AliAnalysisTaskFemtoESE::AliAnalysisTaskFemtoESE() : 
-  AliAnalysisTaskSE(),
-  fAOD(0x0), 
-  fOutputList(0x0),
-  fHelperPID(0x0),
-  fPoolMgr(0x0),
-  fEventCuts(0x0),
-  fTrackCuts(0x0),
-  fFilterBit(128),
-  fSelectBit(AliVEvent::kMB),
-  bIsLHC10h(1),
-  fEventCounter(0),
-  fMixingTracks(10000),
-  fBfield(0.),
-  fMinSepPairEta(0.),
-  fMinSepPairPhi(0.),
-  fShareQuality(0.5),
-  fShareFraction(0.05),
-  nCountSamePairs(0),
-  nCountMixedPairs(0),
-  nCountTracks(0),
-  fMinQPerc(-1000),
-  fMaxQPerc(1000),
-  fQPercDet(0),
-  fEPDet(0),
-  fPsiEPmix(0),
-  fPsiEPmixtemp(0),
-  nKtBins(0),
-  nKtBins1(1),
-  ktBins(0),
-  nEPBins(0),
-  nEPBins1(1),
-  epBins(0),
-  nCentBins(0),
-  nCentBins1(1),
-  centBins(0),
-  nVzBins(0),
-  nVzBins1(1),
-  vzBins(0),
-  hq(0x0),
-  hqmix(0x0), 
-  nqPercBinsLHC11h(1), 
-  qPercBinsLHC11h(0)
+AliAnalysisTaskSE(),
+    fAOD(0x0), 
+    fOutputList(0x0),
+    fHelperPID(0x0),
+    fPoolMgr(0x0),
+    fEventCuts(0x0),
+    fTrackCuts(0x0),
+    fFilterBit(128),
+    fSelectBit(AliVEvent::kMB),
+    bIsLHC10h(1),
+    fEventCounter(0),
+    fMixingTracks(10000),
+    fBfield(0.),
+    fMinSepPairEta(0.),
+    fMinSepPairPhi(0.),
+    fShareQuality(0.5),
+    fShareFraction(0.05),
+    nCountSamePairs(0),
+    nCountMixedPairs(0),
+    nCountTracks(0),
+    fMinQPerc(-1000),
+    fMaxQPerc(1000),
+    fQPercDet(0),
+    fEPDet(0),
+    fPsiEPmix(0),
+    fPsiEPmixtemp(0),
+    nKtBins(0),
+    nKtBins1(1),
+    ktBins(0),
+    nEPBins(0),
+    nEPBins1(1),
+    epBins(0),
+    nCentBins(0),
+    nCentBins1(1),
+    centBins(0),
+    nVzBins(0),
+    nVzBins1(1),
+    vzBins(0),
+    hq(0x0),
+    hqmix(0x0), 
+    nqPercBinsLHC11h(1), 
+    qPercBinsLHC11h(0),
+    hpx(0x0),
+    hpy(0x0),
+    hpz(0x0),
+    hpt(0x0),
+    hE(0x0),
+    hphieta(0x0),
+    hphieta_pid(0x0),
+    hpt_pid(0x0),
+    hvzcent(0x0),
+    hcent(0x0),
+    hcentn(0x0),
+    hphistaretapair10(0x0),
+    hphistaretapair16(0x0),
+    hphistaretapair10a(0x0),
+    hphistaretapair16a(0x0),
+    hphistaretapair10b(0x0),
+    hphistaretapair16b(0x0),
+    hphietapair(0x0),
+    hphietapair2(0x0),
+    hpidid(0x0),
+    hkt(0x0),
+    hktcheck(0x0),
+    hkt3(0x0),
+    hdcaxy(0x0),
+    hdcaz(0x0),
+    hsharequal(0x0),
+    hsharefrac(0x0),
+    hsharequalmix(0x0),
+    hsharefracmix(0x0),
+    hPsiTPC(0x0),
+    hPsiV0A(0x0),
+    hPsiV0C(0x0),
+    hCheckEPA(0x0),
+    hCheckEPC(0x0),
+    hCheckEPmix(0x0),
+    hcentq(0x0),
+    hMixedDist(0x0),
+    hQvecV0A(0x0),
+    hQvecV0C(0x0),
+    hresV0ATPC(0x0),
+    hresV0CTPC(0x0),
+    hresV0AV0C(0x0),
+    hktbins(0x0),
+    hcentbins(0x0),
+    hepbins(0x0)
 {
 }
 
@@ -140,7 +187,52 @@ AliAnalysisTaskFemtoESE::AliAnalysisTaskFemtoESE(const char* name) :
   hq(0x0),
   hqmix(0x0),
   nqPercBinsLHC11h(1), 
-  qPercBinsLHC11h(0)
+  qPercBinsLHC11h(0),
+  hpx(0x0),
+  hpy(0x0),
+  hpz(0x0),
+  hpt(0x0),
+  hE(0x0),
+  hphieta(0x0),
+  hphieta_pid(0x0),
+  hpt_pid(0x0),
+  hvzcent(0x0),
+  hcent(0x0),
+  hcentn(0x0),
+  hphistaretapair10(0x0),
+  hphistaretapair16(0x0),
+  hphistaretapair10a(0x0),
+  hphistaretapair16a(0x0),
+  hphistaretapair10b(0x0),
+  hphistaretapair16b(0x0),
+  hphietapair(0x0),
+  hphietapair2(0x0),
+  hpidid(0x0),
+  hkt(0x0),
+  hktcheck(0x0),
+  hkt3(0x0),
+  hdcaxy(0x0),
+  hdcaz(0x0),
+  hsharequal(0x0),
+  hsharefrac(0x0),
+  hsharequalmix(0x0),
+  hsharefracmix(0x0),
+  hPsiTPC(0x0),
+  hPsiV0A(0x0),
+  hPsiV0C(0x0),
+  hCheckEPA(0x0),
+  hCheckEPC(0x0),
+  hCheckEPmix(0x0),
+  hcentq(0x0),
+  hMixedDist(0x0),
+  hQvecV0A(0x0),
+  hQvecV0C(0x0),
+  hresV0ATPC(0x0),
+  hresV0CTPC(0x0),
+  hresV0AV0C(0x0),
+  hktbins(0x0),
+  hcentbins(0x0),
+  hepbins(0x0)
 {
   
   Printf("*******************************************");
@@ -217,7 +309,52 @@ AliAnalysisTaskFemtoESE::AliAnalysisTaskFemtoESE(const AliAnalysisTaskFemtoESE &
   hq(0x0),
   hqmix(0x0),
   nqPercBinsLHC11h(1), 
-  qPercBinsLHC11h(0)
+  qPercBinsLHC11h(0),
+  hpx(0x0),
+  hpy(0x0),
+  hpz(0x0),
+  hpt(0x0),
+  hE(0x0),
+  hphieta(0x0),
+  hphieta_pid(0x0),
+  hpt_pid(0x0),
+  hvzcent(0x0),
+  hcent(0x0),
+  hcentn(0x0),
+  hphistaretapair10(0x0),
+  hphistaretapair16(0x0),
+  hphistaretapair10a(0x0),
+  hphistaretapair16a(0x0),
+  hphistaretapair10b(0x0),
+  hphistaretapair16b(0x0),
+  hphietapair(0x0),
+  hphietapair2(0x0),
+  hpidid(0x0),
+  hkt(0x0),
+  hktcheck(0x0),
+  hkt3(0x0),
+  hdcaxy(0x0),
+  hdcaz(0x0),
+  hsharequal(0x0),
+  hsharefrac(0x0),
+  hsharequalmix(0x0),
+  hsharefracmix(0x0),
+  hPsiTPC(0x0),
+  hPsiV0A(0x0),
+  hPsiV0C(0x0),
+  hCheckEPA(0x0),
+  hCheckEPC(0x0),
+  hCheckEPmix(0x0),
+  hcentq(0x0),
+  hMixedDist(0x0),
+  hQvecV0A(0x0),
+  hQvecV0C(0x0),
+  hresV0ATPC(0x0),
+  hresV0CTPC(0x0),
+  hresV0AV0C(0x0),
+  hktbins(0x0),
+  hcentbins(0x0),
+  hepbins(0x0)
 {
   /* do nothing yet */  
 } 
@@ -237,163 +374,163 @@ void AliAnalysisTaskFemtoESE::UserCreateOutputObjects()
   if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");
   if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");
   if (!fHelperPID)  AliFatal("HelperPID should be set in the steering macro");
-  
+
   fOutputList = new TList();
   fOutputList->SetOwner();
   fOutputList->SetName("fOutputList");
   
-  TH1D *hpx = new TH1D("hpx","px",200,-2,2);
+  hpx = new TH1D("hpx","px",200,-2,2);
   hpx->GetXaxis()->SetTitle("p_{x}");
   fOutputList->Add(hpx);
-  TH1D *hpy = new TH1D("hpy","py",200,-2,2);
+  hpy = new TH1D("hpy","py",200,-2,2);
   hpy->GetXaxis()->SetTitle("p_{y}");
   fOutputList->Add(hpy);
-  TH1D *hpz = new TH1D("hpz","pz",200,-2,2);
+  hpz = new TH1D("hpz","pz",200,-2,2);
   hpz->GetXaxis()->SetTitle("p_{z}");
   fOutputList->Add(hpz);
-  TH1D *hpt = new TH1D("hpt","pt",100,0,2);
+  hpt = new TH1D("hpt","pt",100,0,2);
   hpt->GetXaxis()->SetTitle("p_{t}");
   fOutputList->Add(hpt);
-  TH1D *hE = new TH1D("hE","E",100,0,2);
+  hE = new TH1D("hE","E",100,0,2);
   hE->GetXaxis()->SetTitle("E");
   fOutputList->Add(hE);
-  TH2D *hphieta = new TH2D("hphieta","track #varphi vs #eta",100,0,2*TMath::Pi(),80,-0.8,0.8);
+  hphieta = new TH2D("hphieta","track #varphi vs #eta",100,0,2*TMath::Pi(),80,-0.8,0.8);
   hphieta->GetXaxis()->SetTitle("#varphi");
   hphieta->GetYaxis()->SetTitle("#eta");
   fOutputList->Add(hphieta);
-  TH2D *hphieta_pid = new TH2D("hphieta_pid","PID check -- #Delta#varphi vs #Delta#eta",100,-0.3,0.3,100,-0.3,0.3);
+  hphieta_pid = new TH2D("hphieta_pid","PID check -- #Delta#varphi vs #Delta#eta",100,-0.3,0.3,100,-0.3,0.3);
   hphieta_pid->GetXaxis()->SetTitle("#Delta#varphi");
   hphieta_pid->GetYaxis()->SetTitle("#Delta#eta");
   fOutputList->Add(hphieta_pid);
-  TH1D *hpt_pid = new TH1D("hpt_pid","PID check -- #Delta p_{t}",100,-0.5,0.5);
+  hpt_pid = new TH1D("hpt_pid","PID check -- #Delta p_{t}",100,-0.5,0.5);
   hpt_pid->GetXaxis()->SetTitle("#Delta p_{t}");
   fOutputList->Add(hpt_pid);
-  TH2D *hvzcent = new TH2D("hvzcent","vz vs cent",nVzBins,vzBins,nCentBins,centBins);
+  hvzcent = new TH2D("hvzcent","vz vs cent",nVzBins,vzBins,nCentBins,centBins);
   hvzcent->GetXaxis()->SetTitle("v_{z}");
   hvzcent->GetYaxis()->SetTitle("centrality");
   fOutputList->Add(hvzcent);
-  TH1D *hcent = new TH1D("hcent","cent",50,0,50);
+  hcent = new TH1D("hcent","cent",50,0,200);
   hcent->GetXaxis()->SetTitle("centrality");
   fOutputList->Add(hcent);
-  TH2D *hcentn = new TH2D("hcentn","cent vs npions",50,0,50,100,0,2000);
+  hcentn = new TH2D("hcentn","cent vs npions",50,0,50,100,0,2000);
   hcentn->GetXaxis()->SetTitle("Centrality");
   hcentn->GetYaxis()->SetTitle("Number of pions");
   fOutputList->Add(hcentn);
-  TH3D *hphistaretapair10 = new TH3D("hphistaretapair10","pair #Delta#varphi* vs #Delta#eta at r=1.0m",100,-0.15,0.15,100,-0.1,0.1,10,0,1);
+  hphistaretapair10 = new TH3D("hphistaretapair10","pair #Delta#varphi* vs #Delta#eta at r=1.0m",100,-0.15,0.15,100,-0.1,0.1,10,0,1);
   hphistaretapair10->GetXaxis()->SetTitle("#Delta#varphi*");
   hphistaretapair10->GetYaxis()->SetTitle("#Delta#eta");
   hphistaretapair10->GetZaxis()->SetTitle("k_{T}");
   fOutputList->Add(hphistaretapair10);
-  TH3D *hphistaretapair16 = new TH3D("hphistaretapair16","pair #Delta#varphi* vs #Delta#eta at r=1.6m",100,-0.15,0.15,100,-0.1,0.1,10,0,1);
+  hphistaretapair16 = new TH3D("hphistaretapair16","pair #Delta#varphi* vs #Delta#eta at r=1.6m",100,-0.15,0.15,100,-0.1,0.1,10,0,1);
   hphistaretapair16->GetXaxis()->SetTitle("#Delta#varphi*");
   hphistaretapair16->GetYaxis()->SetTitle("#Delta#eta");
   hphistaretapair16->GetZaxis()->SetTitle("k_{T}");
   fOutputList->Add(hphistaretapair16);
-  TH3D *hphistaretapair10a = new TH3D("hphistaretapair10a","pair #Delta#varphi* vs #Delta#eta at r=1.0m",100,-0.15,0.15,100,-0.1,0.1,10,0,1);
+  hphistaretapair10a = new TH3D("hphistaretapair10a","pair #Delta#varphi* vs #Delta#eta at r=1.0m",100,-0.15,0.15,100,-0.1,0.1,10,0,1);
   hphistaretapair10a->GetXaxis()->SetTitle("#Delta#varphi*");
   hphistaretapair10a->GetYaxis()->SetTitle("#Delta#eta");
   hphistaretapair10a->GetZaxis()->SetTitle("k_{T}");
   fOutputList->Add(hphistaretapair10a);
-  TH3D *hphistaretapair16a = new TH3D("hphistaretapair16a","pair #Delta#varphi* vs #Delta#eta at r=1.6m",100,-0.15,0.15,100,-0.1,0.1,10,0,1);
+  hphistaretapair16a = new TH3D("hphistaretapair16a","pair #Delta#varphi* vs #Delta#eta at r=1.6m",100,-0.15,0.15,100,-0.1,0.1,10,0,1);
   hphistaretapair16a->GetXaxis()->SetTitle("#Delta#varphi*");
   hphistaretapair16a->GetYaxis()->SetTitle("#Delta#eta");
   hphistaretapair16a->GetZaxis()->SetTitle("k_{T}");
   fOutputList->Add(hphistaretapair16a);
-  TH3D *hphistaretapair10b = new TH3D("hphistaretapair10b","pair #Delta#varphi* vs #Delta#eta at r=1.0m",100,-TMath::Pi(),TMath::Pi(),100,-1.6,1.6,10,0,1);
+  hphistaretapair10b = new TH3D("hphistaretapair10b","pair #Delta#varphi* vs #Delta#eta at r=1.0m",100,-TMath::Pi(),TMath::Pi(),100,-1.6,1.6,10,0,1);
   hphistaretapair10b->GetXaxis()->SetTitle("#Delta#varphi*");
   hphistaretapair10b->GetYaxis()->SetTitle("#Delta#eta");
   hphistaretapair10b->GetZaxis()->SetTitle("k_{T}");
   fOutputList->Add(hphistaretapair10b);
-  TH3D *hphistaretapair16b = new TH3D("hphistaretapair16b","pair #Delta#varphi* vs #Delta#eta at r=1.6m",100,-TMath::Pi(),TMath::Pi(),100,-1.6,1.6,10,0,1);
+  hphistaretapair16b = new TH3D("hphistaretapair16b","pair #Delta#varphi* vs #Delta#eta at r=1.6m",100,-TMath::Pi(),TMath::Pi(),100,-1.6,1.6,10,0,1);
   hphistaretapair16b->GetXaxis()->SetTitle("#Delta#varphi*");
   hphistaretapair16b->GetYaxis()->SetTitle("#Delta#eta");
   hphistaretapair16b->GetZaxis()->SetTitle("k_{T}");
   fOutputList->Add(hphistaretapair16b);
-  TH3D *hphietapair = new TH3D("hphietapair","pair #Delta#varphi vs #Delta#eta",100,-0.1,0.1,100,-0.1,0.1,10,0,1);
+  hphietapair = new TH3D("hphietapair","pair #Delta#varphi vs #Delta#eta",100,-0.1,0.1,100,-0.1,0.1,10,0,1);
   hphietapair->GetXaxis()->SetTitle("#Delta#varphi");
   hphietapair->GetYaxis()->SetTitle("#Delta#eta");
   hphietapair->GetZaxis()->SetTitle("k_{T}");
   fOutputList->Add(hphietapair);
-  TH3D *hphietapair2 = new TH3D("hphietapair2","pair #varphi vs #eta",100,-TMath::Pi(),TMath::Pi(),100,-1.6,1.6,10,0,1);
+  hphietapair2 = new TH3D("hphietapair2","pair #varphi vs #eta",100,-TMath::Pi(),TMath::Pi(),100,-1.6,1.6,10,0,1);
   hphietapair2->GetXaxis()->SetTitle("#Delta#varphi");
   hphietapair2->GetYaxis()->SetTitle("#eta");
   hphietapair2->GetZaxis()->SetTitle("k_{T}");
   fOutputList->Add(hphietapair2);
-  TH1D *hpidid = new TH1D("hpidid","pid id",9,-4.5,4.5);
+  hpidid = new TH1D("hpidid","pid id",9,-4.5,4.5);
   hpidid->GetXaxis()->SetTitle("track PID ID");
   fOutputList->Add(hpidid);
-  TH1D *hkt = new TH1D("hkt","k_{T}",100,0,2);
+  hkt = new TH1D("hkt","k_{T}",100,0,2);
   hkt->GetXaxis()->SetTitle("k_{T}");
   fOutputList->Add(hkt);
-  TH1D *hktcheck = new TH1D("hktcheck","k_{T} check",50,0,1);
+  hktcheck = new TH1D("hktcheck","k_{T} check",50,0,1);
   hktcheck->GetXaxis()->SetTitle("k_{T}");
   fOutputList->Add(hktcheck);
-  TH3D *hkt3 = new TH3D("hkt3","kt vs pt",50,0,1,50,0,5,50,0,5);
+  hkt3 = new TH3D("hkt3","kt vs pt",50,0,1,50,0,5,50,0,5);
   hkt3->GetXaxis()->SetTitle("k_{T}");
   hkt3->GetYaxis()->SetTitle("p_{T,1}");
   hkt3->GetZaxis()->SetTitle("p_{T,2}");
   fOutputList->Add(hkt3);
-  TH2D *hdcaxy = new TH2D("hdcaxy","DCA xy",100,-5,5,100,-5,5);
+  hdcaxy = new TH2D("hdcaxy","DCA xy",100,-5,5,100,-5,5);
   hdcaxy->GetXaxis()->SetTitle("DCA x");
   hdcaxy->GetYaxis()->SetTitle("DCA y");
   fOutputList->Add(hdcaxy);
-  TH1D *hdcaz = new TH1D("hdcaz","DCA z",100,-5,5);
+  hdcaz = new TH1D("hdcaz","DCA z",100,-5,5);
   hdcaz->GetXaxis()->SetTitle("DCA z");
   fOutputList->Add(hdcaz);
-  TH1D *hsharequal = new TH1D("hsharequal","Share Quality",102,-1.02,1.02);
+  hsharequal = new TH1D("hsharequal","Share Quality",102,-1.02,1.02);
   hsharequal->GetXaxis()->SetTitle("Share Quality");
   fOutputList->Add(hsharequal);
-  TH1D *hsharefrac = new TH1D("hsharefrac","Share Fraction",100,0,1);
+  hsharefrac = new TH1D("hsharefrac","Share Fraction",100,0,1);
   hsharefrac->GetXaxis()->SetTitle("Share Fraction");
   fOutputList->Add(hsharefrac);
-  TH1D *hsharequalmix = new TH1D("hsharequalmix","Share Quality -- mixed events",102,-1.02,1.02);
+  hsharequalmix = new TH1D("hsharequalmix","Share Quality -- mixed events",102,-1.02,1.02);
   hsharequalmix->GetXaxis()->SetTitle("Share Quality");
   fOutputList->Add(hsharequalmix);
-  TH1D *hsharefracmix = new TH1D("hsharefracmix","Share Fraction -- mixed events",100,0,1);
+  hsharefracmix = new TH1D("hsharefracmix","Share Fraction -- mixed events",100,0,1);
   hsharefracmix->GetXaxis()->SetTitle("Share Fraction");
   fOutputList->Add(hsharefracmix);
-  TH1D *hPsiTPC = new TH1D("hPsiTPC","TPC EP",100,-1*TMath::Pi(),TMath::Pi());
+  hPsiTPC = new TH1D("hPsiTPC","TPC EP",100,-1*TMath::Pi(),TMath::Pi());
   hPsiTPC->GetXaxis()->SetTitle("#Psi{TPC}");
   fOutputList->Add(hPsiTPC);
-  TH1D *hPsiV0A = new TH1D("hPsiV0A","V0A EP",100,-1*TMath::Pi(),TMath::Pi());
+  hPsiV0A = new TH1D("hPsiV0A","V0A EP",100,-1*TMath::Pi(),TMath::Pi());
   hPsiV0A->GetXaxis()->SetTitle("#Psi{V0A}");
   fOutputList->Add(hPsiV0A);
-  TH1D *hPsiV0C = new TH1D("hPsiV0C","V0C EP",100,-1*TMath::Pi(),TMath::Pi());
+  hPsiV0C = new TH1D("hPsiV0C","V0C EP",100,-1*TMath::Pi(),TMath::Pi());
   hPsiV0C->GetXaxis()->SetTitle("#Psi{V0C}");
   fOutputList->Add(hPsiV0C);
-  TH1D *hCheckEPA = new TH1D("hCheckEPA","Check EP V0A",100,-1*TMath::Pi(),TMath::Pi());
+  hCheckEPA = new TH1D("hCheckEPA","Check EP V0A",100,-1*TMath::Pi(),TMath::Pi());
   hCheckEPA->GetXaxis()->SetTitle("PsiV0A - PsiTPC");
   fOutputList->Add(hCheckEPA);
-  TH1D *hCheckEPC = new TH1D("hCheckEPC","Check EP V0C",100,-1*TMath::Pi(),TMath::Pi());
+  hCheckEPC = new TH1D("hCheckEPC","Check EP V0C",100,-1*TMath::Pi(),TMath::Pi());
   hCheckEPC->GetXaxis()->SetTitle("PsiV0C - PsiTPC");
   fOutputList->Add(hCheckEPC);
-  TH2D* hCheckEPmix = new TH2D("hCheckEPmix","Check EP mixed events",100,-1*TMath::Pi(),TMath::Pi(),100,-1*TMath::Pi(),TMath::Pi());
+  hCheckEPmix = new TH2D("hCheckEPmix","Check EP mixed events",100,-1*TMath::Pi(),TMath::Pi(),100,-1*TMath::Pi(),TMath::Pi());
   hCheckEPmix->GetXaxis()->SetTitle("Psi1 - Psi_mix");
   hCheckEPmix->GetYaxis()->SetTitle("Psi1 - Psi2");
   fOutputList->Add(hCheckEPmix);
-  TH2D *hcentq = new TH2D("hcentq","qvec vs cent",100,0,100,50,0,50);
+  hcentq = new TH2D("hcentq","qvec vs cent",100,0,100,50,0,50);
   hcentq->GetXaxis()->SetTitle("q_{2} percentile");
   hcentq->GetYaxis()->SetTitle("centrality");
   fOutputList->Add(hcentq);
-  TH2D* hMixedDist = new TH2D("hMixedDist", ";centrality;tracks;events", 50, 0, 50, 200, 0, fMixingTracks * 1.5);
+  hMixedDist = new TH2D("hMixedDist", ";centrality;tracks;events", 50, 0, 50, 200, 0, fMixingTracks * 1.5);
   fOutputList->Add(hMixedDist);
-  TH2D *hQvecV0A = new TH2D("hQvecV0A","Qvector in V0A",50,0,50,200,0,5);
+  hQvecV0A = new TH2D("hQvecV0A","Qvector in V0A",50,0,50,200,0,5);
   hQvecV0A->GetXaxis()->SetTitle("Centrality");
   hQvecV0A->GetYaxis()->SetTitle("normalized Qvector");
   fOutputList->Add(hQvecV0A);
-  TH2D *hQvecV0C = new TH2D("hQvecV0C","Qvector in V0C",50,0,50,200,0,5);
+  hQvecV0C = new TH2D("hQvecV0C","Qvector in V0C",50,0,50,200,0,5);
   hQvecV0C->GetXaxis()->SetTitle("Centrality");
   hQvecV0C->GetYaxis()->SetTitle("normalized Qvector");
   fOutputList->Add(hQvecV0C);
 
   // resolution histograms
-  TH1D *hresV0ATPC = new TH1D("hresV0ATPC","cent vs cos(2*(V0A-TPC))",nCentBins,centBins);
+  hresV0ATPC = new TH1D("hresV0ATPC","cent vs cos(2*(V0A-TPC))",nCentBins,centBins);
   hresV0ATPC->GetXaxis()->SetTitle("centrality");
   fOutputList->Add(hresV0ATPC);
-  TH1D *hresV0CTPC = new TH1D("hresV0CTPC","cent vs cos(2*(V0C-TPC))",nCentBins,centBins);
+  hresV0CTPC = new TH1D("hresV0CTPC","cent vs cos(2*(V0C-TPC))",nCentBins,centBins);
   hresV0CTPC->GetXaxis()->SetTitle("centrality");
   fOutputList->Add(hresV0CTPC);
-  TH1D *hresV0AV0C = new TH1D("hresV0AV0C","cent vs cos(2*(V0A-V0C))",nCentBins,centBins);
+  hresV0AV0C = new TH1D("hresV0AV0C","cent vs cos(2*(V0A-V0C))",nCentBins,centBins);
   hresV0AV0C->GetXaxis()->SetTitle("centrality");
   fOutputList->Add(hresV0AV0C);
 
@@ -418,11 +555,11 @@ void AliAnalysisTaskFemtoESE::UserCreateOutputObjects()
     }
 
   // create dummy histograms which just hold the values of the kt, cent, ep bin edges
-  TH1F* hktbins = new TH1F("hktbins","kt bins",nKtBins,ktBins);
+  hktbins = new TH1F("hktbins","kt bins",nKtBins,ktBins);
   fOutputList->Add(hktbins);
-  TH1F* hcentbins = new TH1F("hcentbins","cent bins",nCentBins,centBins);
+  hcentbins = new TH1F("hcentbins","cent bins",nCentBins,centBins);
   fOutputList->Add(hcentbins);
-  TH1F* hepbins = new TH1F("hepbins","ep bins",nEPBins,epBins);
+  hepbins = new TH1F("hepbins","ep bins",nEPBins,epBins);
   fOutputList->Add(hepbins);
 
   Printf("************************");
@@ -433,24 +570,30 @@ void AliAnalysisTaskFemtoESE::UserCreateOutputObjects()
   vertex[0] = vertex[1] = vertex[2] = 0.;
 
   // event mixing pool
+  fPoolMgr = new AliEventPoolManager*[2];
   Int_t poolsize = 1000;
-  fPoolMgr = new AliEventPoolManager(poolsize, fMixingTracks, nCentBins, centBins, nVzBins, vzBins);
-  fPoolMgr->SetTargetValues(fMixingTracks, 0.1, 5); // check these values
-
+  fPoolMgr[0] = new AliEventPoolManager(poolsize, fMixingTracks, nCentBins, centBins, nVzBins, vzBins);
+  fPoolMgr[0]->SetTargetValues(fMixingTracks, 0.1, 5); // check these values
+  fPoolMgr[1] = new AliEventPoolManager(poolsize, fMixingTracks, nCentBins, centBins, nVzBins, vzBins);
+  fPoolMgr[1]->SetTargetValues(fMixingTracks, 0.1, 5); // check these values
 
   nCountSamePairs = 0;
   nCountMixedPairs = 0;
   nCountTracks = 0;
 
+  //stopwatch = new TStopwatch();
+  //stopwatch->Start();
+
   PostData(1, fOutputList);
   PostData(2, fHelperPID);
   PostData(3, fEventCuts);
   PostData(4, fTrackCuts);
+
 }
 
 //________________________________________________________________________
 void AliAnalysisTaskFemtoESE::UserExec(Option_t *) 
-{
+{ 
   // Main loop
   // Called for each event
 
@@ -492,8 +635,11 @@ void AliAnalysisTaskFemtoESE::UserExec(Option_t *)
   //ProcInfo_t procInfo;
   //gSystem->GetProcInfo(&procInfo);
   //Printf("ResMem %ld VMem %ld", procInfo.fMemResident, procInfo.fMemVirtual);
+  //stopwatch->Stop();
+  //Printf("%lf   %lf",stopwatch->RealTime(),stopwatch->CpuTime());
+  //stopwatch->Start();
 
- // get event plane from V0's
+  // get event plane from V0's
   if(!fEventCuts->IsSelected(fAOD,fTrackCuts)) {Printf("Error! Event not accepted by AliAODSpectraEventCuts!"); return;}
   Double_t psiV0A = fEventCuts->GetPsiV0A();
   Double_t psiV0C = fEventCuts->GetPsiV0C();
@@ -504,53 +650,22 @@ void AliAnalysisTaskFemtoESE::UserExec(Option_t *)
       if(fQPercDet == 0) qperc = GetQPercLHC11h(fEventCuts->GetqV0A());
       if(fQPercDet == 1) qperc = GetQPercLHC11h(fEventCuts->GetqV0C());
       //Printf("q vector = %lf  percentile = %lf",fEventCuts->GetqV0A(),qperc);
-      ((TH2D*)fOutputList->FindObject("hQvecV0A"))->Fill(centralityPercentile,fEventCuts->GetqV0A());
-      ((TH2D*)fOutputList->FindObject("hQvecV0C"))->Fill(centralityPercentile,fEventCuts->GetqV0C());
+      hQvecV0A->Fill(centralityPercentile,fEventCuts->GetqV0A());
+      hQvecV0C->Fill(centralityPercentile,fEventCuts->GetqV0C());
     }
   if(psiV0A == -999) return;
   if(psiV0C == -999) return;
   if(qperc < fMinQPerc || qperc > fMaxQPerc) return;
 
-
-  TH1D* hpx = (TH1D*)fOutputList->FindObject("hpx");
-  TH1D* hpy = (TH1D*)fOutputList->FindObject("hpy");
-  TH1D* hpz = (TH1D*)fOutputList->FindObject("hpz");
-  TH1D* hpt = (TH1D*)fOutputList->FindObject("hpt");
-  TH1D* hE = (TH1D*)fOutputList->FindObject("hE");
-  TH2D* hphieta = (TH2D*)fOutputList->FindObject("hphieta");
-  TH2D* hphieta_pid = (TH2D*)fOutputList->FindObject("hphieta_pid");
-  TH1D* hpt_pid = (TH1D*)fOutputList->FindObject("hpt_pid");
-  TH2D* hvzcent = (TH2D*)fOutputList->FindObject("hvzcent");
-  TH1D* hcent = (TH1D*)fOutputList->FindObject("hcent");
-  TH2D* hcentn = (TH2D*)fOutputList->FindObject("hcentn");
-  TH3D* hphistaretapair10 = (TH3D*)fOutputList->FindObject("hphistaretapair10");
-  TH3D* hphistaretapair16 = (TH3D*)fOutputList->FindObject("hphistaretapair16");
-  TH3D* hphistaretapair10a = (TH3D*)fOutputList->FindObject("hphistaretapair10a");
-  TH3D* hphistaretapair16a = (TH3D*)fOutputList->FindObject("hphistaretapair16a");
-  TH3D* hphistaretapair10b = (TH3D*)fOutputList->FindObject("hphistaretapair10b");
-  TH3D* hphistaretapair16b = (TH3D*)fOutputList->FindObject("hphistaretapair16b");
-  TH3D* hphietapair = (TH3D*)fOutputList->FindObject("hphietapair");
-  TH3D* hphietapair2 = (TH3D*)fOutputList->FindObject("hphietapair2");
-  TH1D* hpidid = (TH1D*)fOutputList->FindObject("hpidid");
-  TH1D* hkt = (TH1D*)fOutputList->FindObject("hkt");
-  TH1D* hktcheck = (TH1D*)fOutputList->FindObject("hktcheck");
-  TH3D* hkt3 = (TH3D*)fOutputList->FindObject("hkt3");
-  TH1D* hPsiTPC = (TH1D*)fOutputList->FindObject("hPsiTPC");
-  TH1D* hPsiV0A = (TH1D*)fOutputList->FindObject("hPsiV0A");
-  TH1D* hPsiV0C = (TH1D*)fOutputList->FindObject("hPsiV0C");
-  TH1D* hCheckEPA = (TH1D*)fOutputList->FindObject("hCheckEPA");
-  TH1D* hCheckEPC = (TH1D*)fOutputList->FindObject("hCheckEPC");
-  TH2D* hCheckEPmix = (TH2D*)fOutputList->FindObject("hCheckEPmix");
-  TH2D* hcentq = (TH2D*)fOutputList->FindObject("hcentq");
-  TH2D* hMixedDist = (TH2D*)fOutputList->FindObject("hMixedDist");
-  TH1D* hresV0ATPC = (TH1D*)fOutputList->FindObject("hresV0ATPC");
-  TH1D* hresV0CTPC = (TH1D*)fOutputList->FindObject("hresV0CTPC");
-  TH1D* hresV0AV0C = (TH1D*)fOutputList->FindObject("hresV0AV0C");
+  Double_t psiEP = psiV0A;
+  if(fEPDet==1) psiEP = psiV0C;
 
   Double_t sin2phi = 0, cos2phi = 0;
 
-  TObjArray* tracks = new TObjArray();
-  //tracks->SetOwner(kTRUE);
+  TObjArray* tracksPos = new TObjArray();
+  tracksPos->SetOwner(kTRUE);
+  TObjArray* tracksNeg = new TObjArray();
+  tracksNeg->SetOwner(kTRUE);
 
   // Track loop -- select pions
   for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) {
@@ -574,20 +689,37 @@ void AliAnalysisTaskFemtoESE::UserExec(Option_t *)
     hpidid->Fill((trackPID+1)*aodtrack->Charge());
 
     // select pions
-    if(trackPID==0) 
-      tracks->Add(aodtrack);
-
-      // check event plane angle using tracks in the TPC
-      if(aodtrack->Pt() < 2 && aodtrack->Pt() > 0.2)
-       {
-         sin2phi += (aodtrack->Pt())*sin(2*aodtrack->Phi());
-         cos2phi += (aodtrack->Pt())*cos(2*aodtrack->Phi());
-       }
+    if(trackPID==0)
+      {
+       AliFemtoESEBasicParticle* particle = new AliFemtoESEBasicParticle(sqrt(pow(aodtrack->P(),2)+pow(0.13957, 2)),aodtrack->Px(),aodtrack->Py(),aodtrack->Pz(),aodtrack->Charge(),aodtrack->Phi(),aodtrack->Eta());
+       particle->SetPsiEP(psiEP);
+       particle->SetTPCClusterMap(aodtrack->GetTPCClusterMap());
+       particle->SetTPCSharedMap(aodtrack->GetTPCSharedMap());
+
+       if(particle->Charge()>0)
+         tracksPos->Add(particle);
+       if(particle->Charge()<0)
+         tracksNeg->Add(particle);
+
+       // track qa plots
+       hpx->Fill(particle->Px());
+       hpy->Fill(particle->Py());
+       hpz->Fill(particle->Pz());
+       hpt->Fill(particle->Pt());
+       hE->Fill(particle->E());
+       hphieta->Fill(particle->Phi(),particle->Eta());
+      }
+    // check event plane angle using tracks in the TPC
+    if(aodtrack->Pt() < 2 && aodtrack->Pt() > 0.2)
+      {
+       sin2phi += (aodtrack->Pt())*sin(2*aodtrack->Phi());
+       cos2phi += (aodtrack->Pt())*cos(2*aodtrack->Phi());
+      }
 
   }
   // end track loop
 
-  Int_t ntracks = tracks->GetEntriesFast();
+  Int_t ntracks = tracksPos->GetEntriesFast()+tracksNeg->GetEntriesFast();
 
   // get EP from TPC, just to check
   Double_t psiTPC = 0.;
@@ -595,9 +727,6 @@ void AliAnalysisTaskFemtoESE::UserExec(Option_t *)
     psiTPC = 0.5*atan2(sin2phi,cos2phi);
   else return;
 
-  Double_t psiEP = psiV0A;
-  if(fEPDet==1) psiEP = psiV0C;
-
   hPsiTPC->Fill(psiTPC);
   hPsiV0A->Fill(psiV0A);
   hPsiV0C->Fill(psiV0C);
@@ -613,13 +742,6 @@ void AliAnalysisTaskFemtoESE::UserExec(Option_t *)
 
   hcentq->Fill(qperc,centralityPercentile);
 
-
-  Double_t kt = 0;
-  Double_t qout=0, qside=0, qlong=0;
-  Double_t pVect1[4] = {0,0,0,0};
-  Double_t pVect2[4] = {0,0,0,0};
-  Int_t k, e, c; //bin indices for histograms
-
   nCountTracks += ntracks;
   //cout << "Found " << ntracks << " pion tracks..." << endl;
 
@@ -632,46 +754,70 @@ void AliAnalysisTaskFemtoESE::UserExec(Option_t *)
   hresV0CTPC->Fill(centralityPercentile,cos(2*(psiV0C-psiTPC)));
   hresV0AV0C->Fill(centralityPercentile,cos(2*(psiV0A-psiV0C)));
 
-  AliEventPool* pool = fPoolMgr->GetEventPool(centralityPercentile,zvtx);
-  if (!pool) AliFatal(Form("No pool found for centrality = %f, vz = %f", centralityPercentile, zvtx));
+  AliEventPool* poolPos = fPoolMgr[0]->GetEventPool(centralityPercentile,zvtx);
+  if (!poolPos) AliFatal(Form("No pool found for centrality = %f, vz = %f", centralityPercentile, zvtx));
+  AliEventPool* poolNeg = fPoolMgr[1]->GetEventPool(centralityPercentile,zvtx);
+  if (!poolNeg) AliFatal(Form("No pool found for centrality = %f, vz = %f", centralityPercentile, zvtx));
   //if (!pool) AliFatal(Form("No pool found for centrality = %f, zVtx = %f, Psi_EP = %f", centralityPercentile, zvtx, psiEP));
   //if(pool->IsReady()) hMixedDist->Fill(centralityPercentile, pool->NTracksInPool());
 
+  TrackLoop(tracksPos,poolPos,psiEP,centralityPercentile);
+  TrackLoop(tracksNeg,poolNeg,psiEP,centralityPercentile);
+  //TObjArray* clonedtracks = CloneAndReduceTrackList(tracks,psiEP);
+  poolPos->UpdatePool(tracksPos);
+  poolNeg->UpdatePool(tracksNeg);
+  //cout << "pool contains " << pool->GetCurrentNEvents() << " events and " << pool->NTracksInPool() << " tracks." << endl;
+  //tracks->Clear();
+  
+  //delete tracks;
+
+  // Post output data.
+  PostData(1, fOutputList);
+  PostData(2, fHelperPID);
+  PostData(3, fEventCuts);
+  PostData(4, fTrackCuts);
+}
+
+
+void AliAnalysisTaskFemtoESE::TrackLoop(TObjArray *tracks, AliEventPool *pool, Double_t psiEP, Float_t centralityPercentile)
+{
+  Double_t kt = 0;
+  Double_t qout=0, qside=0, qlong=0;
+  Double_t pVect1[4] = {0,0,0,0};
+  Double_t pVect2[4] = {0,0,0,0};
+  Int_t k, e, c; //bin indices for histograms
+
+  Int_t ntracks = tracks->GetEntriesFast();
 
   for(Int_t j = 0; j < ntracks; j++)
     {
       //cout << endl << j << "   ";
-      AliAODTrack* track1 = (AliAODTrack*)tracks->At(j);
-      pVect1[0]=sqrt(pow(track1->P(),2)+pow(0.13957, 2));
+      AliFemtoESEBasicParticle* track1 = (AliFemtoESEBasicParticle*)tracks->At(j);
+      pVect1[0]=track1->E();
       pVect1[1]=track1->Px();
       pVect1[2]=track1->Py();
       pVect1[3]=track1->Pz();
       //cout << pVect1[0] << "   " << pVect1[1] << "   " <<  pVect1[2] << "   " << pVect1[3] << endl;
 
-      // track qa plots
-      hpx->Fill(pVect1[1]);
-      hpy->Fill(pVect1[2]);
-      hpz->Fill(pVect1[3]);
-      hpt->Fill(track1->Pt());
-      hE->Fill(pVect1[0]);
-      hphieta->Fill(track1->Phi(),track1->Eta());
-
       // same event
       for(Int_t i = j+1; i < ntracks; i++)
        {
-         AliAODTrack* track2 = (AliAODTrack*)tracks->At(i);
+         AliFemtoESEBasicParticle* track2 = (AliFemtoESEBasicParticle*)tracks->At(i);
+
+         Double_t deltaphistar10 = DeltaPhiStar(track1,track2,1.0);
+         Double_t deltaphistar16 = DeltaPhiStar(track1,track2,1.6);
 
-         hphistaretapair10->Fill(DeltaPhiStar(track1,track2,1.0),track1->Eta()-track2->Eta(),kt);
-         hphistaretapair16->Fill(DeltaPhiStar(track1,track2,1.6),track1->Eta()-track2->Eta(),kt);
+         hphistaretapair10->Fill(deltaphistar10,track1->Eta()-track2->Eta(),kt);
+         hphistaretapair16->Fill(deltaphistar16,track1->Eta()-track2->Eta(),kt);
 
          if(!PairCut(track1,track2,kFALSE)) continue;
 
-         hphistaretapair10a->Fill(DeltaPhiStar(track1,track2,1.0),track1->Eta()-track2->Eta(),kt);
-         hphistaretapair16a->Fill(DeltaPhiStar(track1,track2,1.6),track1->Eta()-track2->Eta(),kt);
-         hphistaretapair10b->Fill(DeltaPhiStar(track1,track2,1.0),track1->Eta()-track2->Eta(),kt);
-         hphistaretapair16b->Fill(DeltaPhiStar(track1,track2,1.6),track1->Eta()-track2->Eta(),kt);
+         hphistaretapair10a->Fill(deltaphistar10,track1->Eta()-track2->Eta(),kt);
+         hphistaretapair16a->Fill(deltaphistar16,track1->Eta()-track2->Eta(),kt);
+         hphistaretapair10b->Fill(deltaphistar10,track1->Eta()-track2->Eta(),kt);
+         hphistaretapair16b->Fill(deltaphistar16,track1->Eta()-track2->Eta(),kt);
 
-         pVect2[0]=sqrt(pow(track2->P(),2)+pow(0.13957, 2));
+         pVect2[0]=track2->E();
          pVect2[1]=track2->Px();
          pVect2[2]=track2->Py();
          pVect2[3]=track2->Pz();
@@ -709,6 +855,7 @@ void AliAnalysisTaskFemtoESE::UserExec(Option_t *)
              if(ntracksmix > 0)
                {
                  AliFemtoESEBasicParticle* tracktest = (AliFemtoESEBasicParticle*)bgTracks->UncheckedAt(0);
+                 if(tracktest->Charge() != track1->Charge()){Printf("Charges don't match, there must be a problem..."); continue;}
                  fPsiEPmixtemp = tracktest->GetPsiEP();
                  Double_t dphiEPtest = fPsiEPmixtemp-psiEP;
                  while(dphiEPtest>2*TMath::Pi()) dphiEPtest-=2*TMath::Pi();
@@ -762,19 +909,13 @@ void AliAnalysisTaskFemtoESE::UserExec(Option_t *)
        }
     }
 
-  TObjArray* clonedtracks = CloneAndReduceTrackList(tracks,psiEP);
-  pool->UpdatePool(clonedtracks);
-  //cout << "pool contains " << pool->GetCurrentNEvents() << " events and " << pool->NTracksInPool() << " tracks." << endl;
-  //tracks->Clear();
-  
-  delete tracks;
 
-  // Post output data.
-  PostData(1, fOutputList);
-  PostData(2, fHelperPID);
-  PostData(3, fEventCuts);
-  PostData(4, fTrackCuts);
 }
+
+
+
+
+
 //________________________________________________________________________
 void AliAnalysisTaskFemtoESE::Terminate(Option_t *) 
 {
@@ -798,99 +939,99 @@ void AliAnalysisTaskFemtoESE::Terminate(Option_t *)
 Bool_t AliAnalysisTaskFemtoESE::AcceptPair(AliChaoticityTrackStruct *first, AliChaoticityTrackStruct *second)
 {
   
-  if(fabs(first->fEta-second->fEta) > fMinSepPairEta) return kTRUE;
+if(fabs(first->fEta-second->fEta) > fMinSepPairEta) return kTRUE;
   
-  // propagate through B field to r=1m
-  Float_t phi1 = first->fPhi - asin(first->fCharge*(0.1*fBfield)*0.15/first->fPt);// 0.15 for D=1m
-  if(phi1 > 2*PI) phi1 -= 2*PI;
-  if(phi1 < 0) phi1 += 2*PI;
-  Float_t phi2 = second->fPhi - asin(second->fCharge*(0.1*fBfield)*0.15/second->fPt);// 0.15 for D=1m 
-  if(phi2 > 2*PI) phi2 -= 2*PI;
-  if(phi2 < 0) phi2 += 2*PI;
+// propagate through B field to r=1m
+Float_t phi1 = first->fPhi - asin(first->fCharge*(0.1*fBfield)*0.15/first->fPt);// 0.15 for D=1m
+if(phi1 > 2*PI) phi1 -= 2*PI;
+if(phi1 < 0) phi1 += 2*PI;
+Float_t phi2 = second->fPhi - asin(second->fCharge*(0.1*fBfield)*0.15/second->fPt);// 0.15 for D=1m 
+if(phi2 > 2*PI) phi2 -= 2*PI;
+if(phi2 < 0) phi2 += 2*PI;
   
-  Float_t deltaphi = phi1 - phi2;
-  if(deltaphi > PI) deltaphi -= 2*PI;
-  if(deltaphi < -PI) deltaphi += 2*PI;
-  deltaphi = fabs(deltaphi);
+Float_t deltaphi = phi1 - phi2;
+if(deltaphi > PI) deltaphi -= 2*PI;
+if(deltaphi < -PI) deltaphi += 2*PI;
+deltaphi = fabs(deltaphi);
 
-  if(deltaphi < fMinSepPairPhi) return kFALSE;// Min Separation
+if(deltaphi < fMinSepPairPhi) return kFALSE;// Min Separation
     
   
-  // propagate through B field to r=1.6m
-  phi1 = first->fPhi - asin(first->fCharge*(0.1*fBfield)*0.24/first->fPt);// mine. 0.24 for D=1.6m
-  if(phi1 > 2*PI) phi1 -= 2*PI;
-  if(phi1 < 0) phi1 += 2*PI;
-  phi2 = second->fPhi - asin(second->fCharge*(0.1*fBfield)*0.24/second->fPt);// mine. 0.24 for D=1.6m 
-  if(phi2 > 2*PI) phi2 -= 2*PI;
-  if(phi2 < 0) phi2 += 2*PI;
+// propagate through B field to r=1.6m
+phi1 = first->fPhi - asin(first->fCharge*(0.1*fBfield)*0.24/first->fPt);// mine. 0.24 for D=1.6m
+if(phi1 > 2*PI) phi1 -= 2*PI;
+if(phi1 < 0) phi1 += 2*PI;
+phi2 = second->fPhi - asin(second->fCharge*(0.1*fBfield)*0.24/second->fPt);// mine. 0.24 for D=1.6m 
+if(phi2 > 2*PI) phi2 -= 2*PI;
+if(phi2 < 0) phi2 += 2*PI;
   
-  deltaphi = phi1 - phi2;
-  if(deltaphi > PI) deltaphi -= 2*PI;
-  if(deltaphi < -PI) deltaphi += 2*PI;
-  deltaphi = fabs(deltaphi);
+deltaphi = phi1 - phi2;
+if(deltaphi > PI) deltaphi -= 2*PI;
+if(deltaphi < -PI) deltaphi += 2*PI;
+deltaphi = fabs(deltaphi);
 
-  if(deltaphi < fMinSepPairPhi) return kFALSE;// Min Separation
+if(deltaphi < fMinSepPairPhi) return kFALSE;// Min Separation
   
   
    
-  //
+//
   
-  Int_t ncl1 = first->fClusterMap.GetNbits();
-  Int_t ncl2 = second->fClusterMap.GetNbits();
-  Int_t sumCls = 0; Int_t sumSha = 0; Int_t sumQ = 0;
-  Double_t shfrac = 0; Double_t qfactor = 0;
-  for(Int_t imap = 0; imap < ncl1 && imap < ncl2; imap++) {
-    if (first->fClusterMap.TestBitNumber(imap) && second->fClusterMap.TestBitNumber(imap)) {// Both clusters
-      if (first->fSharedMap.TestBitNumber(imap) && second->fSharedMap.TestBitNumber(imap)) { // Shared
-       sumQ++;
-       sumCls+=2;
-       sumSha+=2;}
-      else {sumQ--; sumCls+=2;}
-    }
-    else if (first->fClusterMap.TestBitNumber(imap) || second->fClusterMap.TestBitNumber(imap)) {// Non shared
-      sumQ++;
-      sumCls++;}
-  }
-  if (sumCls>0) {
-    qfactor = sumQ*1.0/sumCls;
-    shfrac = sumSha*1.0/sumCls;
-  }
+Int_t ncl1 = first->fClusterMap.GetNbits();
+Int_t ncl2 = second->fClusterMap.GetNbits();
+Int_t sumCls = 0; Int_t sumSha = 0; Int_t sumQ = 0;
+Double_t shfrac = 0; Double_t qfactor = 0;
+for(Int_t imap = 0; imap < ncl1 && imap < ncl2; imap++) {
+if (first->fClusterMap.TestBitNumber(imap) && second->fClusterMap.TestBitNumber(imap)) {// Both clusters
+if (first->fSharedMap.TestBitNumber(imap) && second->fSharedMap.TestBitNumber(imap)) { // Shared
+sumQ++;
+sumCls+=2;
+sumSha+=2;}
+else {sumQ--; sumCls+=2;}
+}
+else if (first->fClusterMap.TestBitNumber(imap) || second->fClusterMap.TestBitNumber(imap)) {// Non shared
+sumQ++;
+sumCls++;}
+}
+if (sumCls>0) {
+qfactor = sumQ*1.0/sumCls;
+shfrac = sumSha*1.0/sumCls;
+}
   
-  if(qfactor > fShareQuality || shfrac > fShareFraction) return kFALSE;
+if(qfactor > fShareQuality || shfrac > fShareFraction) return kFALSE;
   
   
-  return kTRUE;
+return kTRUE;
   
 
 }
 //________________________________________________________________________
 Float_t AliAnalysisTaskFemtoESE::GamovFactor(Int_t chargeBin1, Int_t chargeBin2, Float_t qinv)
 {
-  Float_t arg = G_Coeff/qinv;
+Float_t arg = G_Coeff/qinv;
   
-  if(chargeBin1==chargeBin2) return (exp(arg)-1)/(arg);
-  else {return (exp(-arg)-1)/(-arg);}
+if(chargeBin1==chargeBin2) return (exp(arg)-1)/(arg);
+else {return (exp(-arg)-1)/(-arg);}
   
 }
 //________________________________________________________________________
 void AliAnalysisTaskFemtoESE::Shuffle(Int_t *iarr, Int_t i1, Int_t i2)
 {
-  Int_t j, k;
-  Int_t a = i2 - i1;
-  for (Int_t i = i1; i < i2+1; i++) {
-    j = (Int_t) (gRandom->Rndm() * a);
-    k = iarr[j];
-    iarr[j] = iarr[i];
-    iarr[i] = k;
-  }
+Int_t j, k;
+Int_t a = i2 - i1;
+for (Int_t i = i1; i < i2+1; i++) {
+j = (Int_t) (gRandom->Rndm() * a);
+k = iarr[j];
+iarr[j] = iarr[i];
+iarr[i] = k;
+}
 }
 */
 //________________________________________________________________________
-Double_t AliAnalysisTaskFemtoESE::GetQinv(Double_t track1[], Double_t track2[]){
+Double_t AliAnalysisTaskFemtoESE::GetQinv2(Double_t track1[], Double_t track2[]){
   
-  Double_t qinv=1.0;
-  qinv = sqrt( pow(track1[1]-track2[1],2) + pow(track1[2]-track2[2],2) + pow(track1[3]-track2[3],2) - pow(track1[0]-track2[0],2));
-  return qinv;
+  Double_t qinv2=1.0;
+  qinv2 = pow(track1[1]-track2[1],2) + pow(track1[2]-track2[2],2) + pow(track1[3]-track2[3],2) - pow(track1[0]-track2[0],2);
+  return qinv2;
 }
 
 //________________________________________________________________________
@@ -934,46 +1075,46 @@ Bool_t AliAnalysisTaskFemtoESE::EventCut(/*AliAODEvent* fevent*/){
   Bool_t isSelected3 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kSemiCentral);
   if(!isSelected1 && !isSelected2 && !isSelected3) {return kFALSE;}
 
-      /*  
-    // Pile-up rejection
-    AliAnalysisUtils *AnaUtil=new AliAnalysisUtils();
-    if(!fPbPbcase) AnaUtil->SetUseMVPlpSelection(kTRUE);// use Multi-Vertex tool for pp and pPb
-    else AnaUtil->SetUseMVPlpSelection(kFALSE);
-    Bool_t pileUpCase=AnaUtil->IsPileUpEvent(fAOD); 
-    if(pileUpCase) return;
-
-    // Vertexing
-    primaryVertexAOD = fAOD->GetPrimaryVertex();
-    vertex[0]=primaryVertexAOD->GetX(); vertex[1]=primaryVertexAOD->GetY(); vertex[2]=primaryVertexAOD->GetZ();
+  /*  
+  // Pile-up rejection
+  AliAnalysisUtils *AnaUtil=new AliAnalysisUtils();
+  if(!fPbPbcase) AnaUtil->SetUseMVPlpSelection(kTRUE);// use Multi-Vertex tool for pp and pPb
+  else AnaUtil->SetUseMVPlpSelection(kFALSE);
+  Bool_t pileUpCase=AnaUtil->IsPileUpEvent(fAOD); 
+  if(pileUpCase) return;
+
+  // Vertexing
+  primaryVertexAOD = fAOD->GetPrimaryVertex();
+  vertex[0]=primaryVertexAOD->GetX(); vertex[1]=primaryVertexAOD->GetY(); vertex[2]=primaryVertexAOD->GetZ();
     
-    if(fabs(vertex[2]) > 10) {cout<<"Zvertex Out of Range. Skip Event"<<endl; return;} // Z-Vertex Cut 
-    ((TH3D*)fOutputList->FindObject("fVertexDist"))->Fill(vertex[0], vertex[1], vertex[2]);
+  if(fabs(vertex[2]) > 10) {cout<<"Zvertex Out of Range. Skip Event"<<endl; return;} // Z-Vertex Cut 
+  ((TH3D*)fOutputList->FindObject("fVertexDist"))->Fill(vertex[0], vertex[1], vertex[2]);
     
-    for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) {
-      AliAODTrack* aodtrack = fAOD->GetTrack(i);
-      if (!aodtrack) continue;
-      AODTracks++;
-      if(!aodtrack->TestFilterBit(BIT(fFilterBit))) continue;// AOD filterBit cut
-      FBTracks++;
-    }
-    ((TH1D*)fOutputList->FindObject("fMultDist2"))->Fill(FBTracks);
+  for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) {
+  AliAODTrack* aodtrack = fAOD->GetTrack(i);
+  if (!aodtrack) continue;
+  AODTracks++;
+  if(!aodtrack->TestFilterBit(BIT(fFilterBit))) continue;// AOD filterBit cut
+  FBTracks++;
+  }
+  ((TH1D*)fOutputList->FindObject("fMultDist2"))->Fill(FBTracks);
 
-    //if(fAOD->IsPileupFromSPD()) {cout<<"PileUpEvent. Skip Event"<<endl; return;} // Old Pile-up cut
-    if(primaryVertexAOD->GetNContributors() < 1) {cout<<"Bad Vertex. Skip Event"<<endl; return;}
+  //if(fAOD->IsPileupFromSPD()) {cout<<"PileUpEvent. Skip Event"<<endl; return;} // Old Pile-up cut
+  if(primaryVertexAOD->GetNContributors() < 1) {cout<<"Bad Vertex. Skip Event"<<endl; return;}
    
-    ((TH1D*)fOutputList->FindObject("fMultDist3"))->Fill(FBTracks);
+  ((TH1D*)fOutputList->FindObject("fMultDist3"))->Fill(FBTracks);
  
-    fBfield = fAOD->GetMagneticField();
+  fBfield = fAOD->GetMagneticField();
     
-    for(Int_t i=0; i<fZvertexBins; i++){
-      if( (vertex[2] >= zstart+i*zstep) && (vertex[2] < zstart+(i+1)*zstep) ){
-       zbin=i;
-       break;
-      }
-    }
-      */
+  for(Int_t i=0; i<fZvertexBins; i++){
+  if( (vertex[2] >= zstart+i*zstep) && (vertex[2] < zstart+(i+1)*zstep) ){
+  zbin=i;
+  break;
+  }
+  }
+  */
 
-    return kTRUE;
+  return kTRUE;
 
 }
 
@@ -992,26 +1133,26 @@ Bool_t AliAnalysisTaskFemtoESE::TrackCut(AliAODTrack* ftrack){
   //ftrack->XYZAtDCA(trackdca);
   //Double_t dcaxy = sqrt( pow(trackpos[0] - vertex[0],2) + pow(trackpos[1] - vertex[1],2));
   //Double_t dcaz = sqrt( pow(trackpos[2] - vertex[2],2));
-  ((TH2D*)fOutputList->FindObject("hdcaxy"))->Fill(trackdca[0],trackdca[1]);
-  ((TH1D*)fOutputList->FindObject("hdcaz"))->Fill(trackdca[2]);
+  hdcaxy->Fill(trackdca[0],trackdca[1]);
+  hdcaz->Fill(trackdca[2]);
   //if(dcaxy > 0.2) return kFALSE;
   //if(dcaz > 0.15) return kFALSE;
 
 
-    //// FilterBit Overlap Check
-    //if(fFilterBit != 7){
-//     Bool_t goodTrackOtherFB = kFALSE;
-//     if(fMCcase && fAOD->GetRunNumber()<=126437) goodTrackOtherFB=kTRUE;// FB7 to FB5 mapping in 10f6a MC does not work
-//     
-//     for (Int_t j = 0; j < fAOD->GetNumberOfTracks(); j++) {
-//       AliAODTrack* aodtrack2 = fAOD->GetTrack(randomIndex[j]);
-//       if(!aodtrack2) continue;
-//       if(!aodtrack2->TestFilterBit(BIT(fFilterBit))) continue;
-//       
-//       if(-(aodtrack->GetID()+1)==aodtrack2->GetID()) {goodTrackOtherFB=kTRUE; break;}
-//       
-//     }
-//     if(!goodTrackOtherFB) continue;
+  //// FilterBit Overlap Check
+  //if(fFilterBit != 7){
+  //   Bool_t goodTrackOtherFB = kFALSE;
+  //   if(fMCcase && fAOD->GetRunNumber()<=126437) goodTrackOtherFB=kTRUE;// FB7 to FB5 mapping in 10f6a MC does not work
+  //   
+  //   for (Int_t j = 0; j < fAOD->GetNumberOfTracks(); j++) {
+  //     AliAODTrack* aodtrack2 = fAOD->GetTrack(randomIndex[j]);
+  //     if(!aodtrack2) continue;
+  //     if(!aodtrack2->TestFilterBit(BIT(fFilterBit))) continue;
+  //     
+  //     if(-(aodtrack->GetID()+1)==aodtrack2->GetID()) {goodTrackOtherFB=kTRUE; break;}
+  //     
+  //   }
+  //   if(!goodTrackOtherFB) continue;
   //    }
 
   return kTRUE;
@@ -1019,16 +1160,16 @@ Bool_t AliAnalysisTaskFemtoESE::TrackCut(AliAODTrack* ftrack){
 }
 
 //________________________________________________________________________
-Bool_t AliAnalysisTaskFemtoESE::PairCut(AliAODTrack* ftrack1, AliAODTrack* ftrack2, Bool_t mix){
+Bool_t AliAnalysisTaskFemtoESE::PairCut(AliFemtoESEBasicParticle* ftrack1, AliFemtoESEBasicParticle* ftrack2, Bool_t mix){
 
   // check for same charge
   if(ftrack2->Charge() != ftrack1->Charge()) return kFALSE;
   
   // qinv cut
-  Double_t trackvec1[4] = {sqrt(pow(ftrack1->P(),2)+pow(0.13957, 2)),ftrack1->Px(),ftrack1->Py(),ftrack1->Pz()};
-  Double_t trackvec2[4] = {sqrt(pow(ftrack2->P(),2)+pow(0.13957, 2)),ftrack2->Px(),ftrack2->Py(),ftrack2->Pz()};
-  Double_t qinv = GetQinv(trackvec1,trackvec2);
-  if(qinv < 0.005) return kFALSE;
+  Double_t trackvec1[4] = {ftrack1->E(),ftrack1->Px(),ftrack1->Py(),ftrack1->Pz()};
+  Double_t trackvec2[4] = {ftrack2->E(),ftrack2->Px(),ftrack2->Py(),ftrack2->Pz()};
+  Double_t qinv2 = GetQinv2(trackvec1,trackvec2);
+  if(qinv2 < 0.000025) return kFALSE; // qinv < 0.005
 
   // deltaEta x deltaPhi* cut
   if(fabs(ftrack1->Eta()-ftrack2->Eta()) < fMinSepPairEta)
@@ -1077,22 +1218,18 @@ Bool_t AliAnalysisTaskFemtoESE::PairCut(AliAODTrack* ftrack1, AliAODTrack* ftrac
   // sumCls -- number of clusters in track 1 + number of clusters in track 2 (clusters in both tracks counted twice)
   // sumSha -- number of shared clusters (counted twice)
   // sumQ -- ?
-  
-  TH1D* hsharequal;
-  TH1D* hsharefrac;
 
   if(!mix)
     {
-      hsharequal = (TH1D*)fOutputList->FindObject("hsharequal");
-      hsharefrac = (TH1D*)fOutputList->FindObject("hsharefrac");
+      hsharequal->Fill(qfactor);
+      hsharefrac->Fill(shfrac);
     }
   else
     {
-     hsharequal = (TH1D*)fOutputList->FindObject("hsharequalmix");
-     hsharefrac = (TH1D*)fOutputList->FindObject("hsharefracmix");
+      hsharequalmix->Fill(qfactor);
+      hsharefracmix->Fill(shfrac);
     }
-  hsharequal->Fill(qfactor);
-  hsharefrac->Fill(shfrac);
+
   
   if(qfactor > fShareQuality || shfrac > fShareFraction) return kFALSE;
   
@@ -1111,8 +1248,8 @@ Double_t AliAnalysisTaskFemtoESE::DeltaPhiStar(AliAODTrack* ftrack1, AliAODTrack
   return deltaphi;
 }
 
-TObjArray* AliAnalysisTaskFemtoESE::CloneAndReduceTrackList(TObjArray* tracks, Double_t psi)
-{
+/*TObjArray* AliAnalysisTaskFemtoESE::CloneAndReduceTrackList(TObjArray* tracks, Double_t psi)
+  {
   // clones a track list by using AliDPhiBasicParticle which uses much less memory (used for event mixing)
 
   TObjArray* tracksClone = new TObjArray;
@@ -1120,17 +1257,17 @@ TObjArray* AliAnalysisTaskFemtoESE::CloneAndReduceTrackList(TObjArray* tracks, D
 
   for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
   {
-    AliAODTrack* particle = (AliAODTrack*) tracks->UncheckedAt(i);
-    AliFemtoESEBasicParticle* copy = new AliFemtoESEBasicParticle(sqrt(pow(particle->P(),2)+pow(0.13957, 2)),particle->Px(),particle->Py(),particle->Pz(),particle->Charge(),particle->Phi(),particle->Eta());
-    copy->SetPsiEP(psi);
-    copy->SetTPCClusterMap(particle->GetTPCClusterMap());
-    copy->SetTPCSharedMap(particle->GetTPCSharedMap());
+  AliAODTrack* particle = (AliAODTrack*) tracks->UncheckedAt(i);
+  AliFemtoESEBasicParticle* copy = new AliFemtoESEBasicParticle(sqrt(pow(particle->P(),2)+pow(0.13957, 2)),particle->Px(),particle->Py(),particle->Pz(),particle->Charge(),particle->Phi(),particle->Eta());
+  copy->SetPsiEP(psi);
+  copy->SetTPCClusterMap(particle->GetTPCClusterMap());
+  copy->SetTPCSharedMap(particle->GetTPCSharedMap());
 
-    tracksClone->Add(copy);
+  tracksClone->Add(copy);
   }
   
   return tracksClone;
-}
+  }*/
 
 Double_t AliAnalysisTaskFemtoESE::GetDeltaPhiEP(Double_t px1, Double_t py1, Double_t px2, Double_t py2, Double_t psi)
 {
index b03f010..d544061 100644 (file)
@@ -21,6 +21,11 @@ class AliEventPoolManager;
 class AliSpectraAODEventCuts;
 class AliSpectraAODTrackCuts;
 class AliAODTrack;
+class AliEventPool;
+
+//class TStopwatch;
+
+class AliFemtoESEBasicParticle;
 
 #include "AliAnalysisTask.h"
 #include "AliAnalysisTaskSE.h"
@@ -44,6 +49,7 @@ class AliAnalysisTaskFemtoESE : public AliAnalysisTaskSE {
 
   virtual void   UserCreateOutputObjects();
   virtual void   UserExec(Option_t *option);
+  void TrackLoop(TObjArray *tracks, AliEventPool *pool, Double_t psiEP, Float_t centralityPercentile);
   virtual void   Terminate(Option_t *);
 
   AliHelperPID* GetHelperPID() { return fHelperPID; }
@@ -77,11 +83,11 @@ class AliAnalysisTaskFemtoESE : public AliAnalysisTaskSE {
   Double_t GetQPercLHC11h(Double_t qvec);
 
  private:
-  Double_t GetQinv(Double_t[], Double_t[]);
+  Double_t GetQinv2(Double_t[], Double_t[]);
   void GetQosl(Double_t[], Double_t[], Double_t&, Double_t&, Double_t&);
   Bool_t TrackCut(AliAODTrack* ftrack);
   Bool_t EventCut(/*AliAODEvent* fevent*/);
-  Bool_t PairCut(AliAODTrack* ftrack1, AliAODTrack* ftrack2, Bool_t mix);
+  Bool_t PairCut(AliFemtoESEBasicParticle* ftrack1, AliFemtoESEBasicParticle* ftrack2, Bool_t mix);
   Double_t DeltaPhiStar(AliAODTrack* ftrack1, AliAODTrack* ftrack2, Double_t r);
   TObjArray* CloneAndReduceTrackList(TObjArray* tracks, Double_t psi);
   Double_t GetDeltaPhiEP(Double_t px1, Double_t py1, Double_t px2, Double_t py2, Double_t psi);
@@ -91,7 +97,7 @@ class AliAnalysisTaskFemtoESE : public AliAnalysisTaskSE {
   TList                  *fOutputList; //! Compact Output list
   //AliPIDResponse         *fPIDResponse; //! PID response object; equivalent to AliAODpidUtil
   AliHelperPID*     fHelperPID;      // points to class for PID
-  AliEventPoolManager*     fPoolMgr;         //! event pool manager
+  AliEventPoolManager**     fPoolMgr;         //[2] event pool manager
   AliSpectraAODEventCuts* fEventCuts;
   AliSpectraAODTrackCuts* fTrackCuts;
 
@@ -141,6 +147,54 @@ class AliAnalysisTaskFemtoESE : public AliAnalysisTaskSE {
   Int_t nqPercBinsLHC11h;
   Double_t* qPercBinsLHC11h; //[nqPercBinsLHC11h]
 
+  //TStopwatch *stopwatch;
+
+  // histograms
+  TH1D *hpx;
+  TH1D *hpy;
+  TH1D *hpz;
+  TH1D *hpt;
+  TH1D *hE;
+  TH2D *hphieta;
+  TH2D *hphieta_pid;
+  TH1D *hpt_pid;
+  TH2D *hvzcent;
+  TH1D *hcent;
+  TH2D *hcentn;
+  TH3D *hphistaretapair10;
+  TH3D *hphistaretapair16;
+  TH3D *hphistaretapair10a;
+  TH3D *hphistaretapair16a;
+  TH3D *hphistaretapair10b;
+  TH3D *hphistaretapair16b;
+  TH3D *hphietapair;
+  TH3D *hphietapair2;
+  TH1D *hpidid;
+  TH1D *hkt;
+  TH1D *hktcheck;
+  TH3D *hkt3;
+  TH2D* hdcaxy;
+  TH1D* hdcaz;
+  TH1D* hsharequal;
+  TH1D* hsharefrac;
+  TH1D* hsharequalmix;
+  TH1D* hsharefracmix;
+  TH1D *hPsiTPC;
+  TH1D *hPsiV0A;
+  TH1D *hPsiV0C;
+  TH1D *hCheckEPA;
+  TH1D *hCheckEPC;
+  TH2D* hCheckEPmix;
+  TH2D *hcentq;
+  TH2D* hMixedDist;
+  TH2D *hQvecV0A;
+  TH2D *hQvecV0C;
+  TH1D *hresV0ATPC;
+  TH1D *hresV0CTPC;
+  TH1D *hresV0AV0C;
+  TH1F* hktbins;
+  TH1F* hcentbins;
+  TH1F* hepbins;
 
   ClassDef(AliAnalysisTaskFemtoESE, 1);
 };