PWGCF/FEMTOSCOPY/ESE/AliAnalysisTaskFemtoESE.cxx

   1 #include <iostream>
   2 #include <math.h>
   3 #include "TChain.h"
   4 #include "TFile.h"
   5 #include "TKey.h"
   6 #include "TObject.h"
   7 #include "TObjString.h"
   8 #include "TList.h"
   9 #include "TTree.h"
  10 #include "TH1F.h"
  11 #include "TH1D.h"
  12 #include "TH2D.h"
  13 #include "TH3D.h"
  14 #include "THn.h"
  15 #include "THnSparse.h"
  16 #include "TProfile.h"
  17 #include "TProfile2D.h"
  18 #include "TCanvas.h"
  19 #include "TRandom3.h"
  20 #include "TF1.h"
  21 #include "TObjectTable.h"
  22 #include <vector>
  23
  24 //#include "TStopwatch.h"
  25
  26 #include "AliAnalysisTask.h"
  27 #include "AliAnalysisManager.h"
  28
  29
  30 #include "AliESDEvent.h"
  31 #include "AliESDInputHandler.h"
  32 #include "AliESDtrackCuts.h"
  33
  34 #include "AliAODEvent.h"
  35 #include "AliAODInputHandler.h"
  36 #include "AliAODMCParticle.h"
  37 //#include "AliAnalysisUtils.h"
  38 #include "AliHelperPID.h"
  39 #include "AliEventPoolManager.h"
  40
  41 #include "AliAnalysisTaskFemtoESE.h"
  42
  43 //#include "AliSpectraAODEventCuts.h"
  44 //#include "AliSpectraAODTrackCuts.h"
  45 //#include "/opt/alice/aliroot/master/src/PWGLF/SPECTRA/PiKaPr/TestAOD/AliSpectraAODEventCuts.h"
  46 //#include "/opt/alice/aliroot/master/src/PWGLF/SPECTRA/PiKaPr/TestAOD/AliSpectraAODTrackCuts.h"
  47 //#include "AliSpectraAODEventCuts.h"
  48 //#include "AliSpectraAODTrackCuts.h"
  49
  50 ClassImp(AliAnalysisTaskFemtoESE)
  51
  52
  53 //________________________________________________________________________
  54 // Default constructor
  55 AliAnalysisTaskFemtoESE::AliAnalysisTaskFemtoESE() :
  56 AliAnalysisTaskSE(),
  57     fAOD(0x0),
  58     fOutputList(0x0),
  59     fHelperPID(0x0),
  60     fPoolMgr(0x0),
  61     fEventCuts(0x0),
  62     fTrackCuts(0x0),
  63     fFilterBit(128),
  64     fSelectBit(AliVEvent::kMB),
  65     bIsLHC10h(1),
  66     fEventCounter(0),
  67     fMixingTracks(10000),
  68     fBfield(0.),
  69     fMinSepPairEta(0.),
  70     fMinSepPairPhi(0.),
  71     fShareQuality(0.5),
  72     fShareFraction(0.05),
  73     nCountSamePairs(0),
  74     nCountMixedPairs(0),
  75     nCountTracks(0),
  76     fMinQPerc(-1000),
  77     fMaxQPerc(1000),
  78     fQPercDet(0),
  79     fEPDet(0),
  80     nKtBins(0),
  81     nKtBins1(1),
  82     ktBins(0),
  83     nEPBins(0),
  84     nEPBins1(1),
  85     epBins(0),
  86     nEPBinsMix(0),
  87     nEPBinsMix1(1),
  88     epBinsMix(0),
  89     nCentBins(0),
  90     nCentBins1(1),
  91     centBins(0),
  92     nVzBins(0),
  93     nVzBins1(1),
  94     vzBins(0),
  95     hq(0x0),
  96     hqmix(0x0),
  97     hqinv(0x0),
  98     nqPercBinsLHC11h(1),
  99     qPercBinsLHC11h(0),
 100     hpx(0x0),
 101     hpy(0x0),
 102     hpz(0x0),
 103     hpt(0x0),
 104     hE(0x0),
 105     hphieta(0x0),
 106     hphieta_pid(0x0),
 107     hpt_pid(0x0),
 108     hvzcent(0x0),
 109     hcent(0x0),
 110     hcentUnweighted(0x0),
 111     hcentn(0x0),
 112     hphistaretapair10(0x0),
 113     hphistaretapair16(0x0),
 114     hphistaretapair10a(0x0),
 115     hphistaretapair16a(0x0),
 116     hphistaretapair10b(0x0),
 117     hphistaretapair16b(0x0),
 118     hphietapair(0x0),
 119     hphietapair2(0x0),
 120     hpidid(0x0),
 121     hkt(0x0),
 122     hktcheck(0x0),
 123     hkt3(0x0),
 124     hdcaxy(0x0),
 125     hdcaz(0x0),
 126     hsharequal(0x0),
 127     hsharefrac(0x0),
 128     hsharequalmix(0x0),
 129     hsharefracmix(0x0),
 130     hPsiTPC(0x0),
 131     hPsiV0A(0x0),
 132     hPsiV0C(0x0),
 133     hShiftTPC(0x0),
 134     hShiftV0A(0x0),
 135     hShiftV0C(0x0),
 136     hPsiMix(0x0),
 137     hCheckEPA(0x0),
 138     hCheckEPC(0x0),
 139     hCheckEPmix(0x0),
 140     hAvDphi(0x0),
 141     hNpairs(0x0),
 142     hPairDphi(0x0),
 143     hPairDphiMix(0x0),
 144     hcentq(0x0),
 145     hMixedDistTracks(0x0),
 146     hMixedDistEvents(0x0),
 147     hQvecV0A(0x0),
 148     hQvecV0C(0x0),
 149     hresV0ATPC(0x0),
 150     hresV0CTPC(0x0),
 151     hresV0AV0C(0x0),
 152     hqinvcheck(0x0),
 153     hktbins(0x0),
 154     hcentbins(0x0),
 155     hepbins(0x0)
 156 {
 157 }
 158
 159 //________________________________________________________________________
 160 // Constructor
 161 AliAnalysisTaskFemtoESE::AliAnalysisTaskFemtoESE(const char* name) :
 162   AliAnalysisTaskSE(name),
 163   fAOD(0x0),
 164   fOutputList(0x0),
 165   fHelperPID(0x0),
 166   fPoolMgr(0x0),
 167   fEventCuts(0x0),
 168   fTrackCuts(0x0),
 169   fFilterBit(128),
 170   fSelectBit(AliVEvent::kMB),
 171   bIsLHC10h(1),
 172   fEventCounter(0),
 173   fMixingTracks(10000),
 174   fBfield(0.),
 175   fMinSepPairEta(0.),
 176   fMinSepPairPhi(0.),
 177   fShareQuality(0.5),
 178   fShareFraction(0.05),
 179   nCountSamePairs(0),
 180   nCountMixedPairs(0),
 181   nCountTracks(0),
 182   fMinQPerc(-1000),
 183   fMaxQPerc(1000),
 184   fQPercDet(0),
 185   fEPDet(0),
 186   nKtBins(0),
 187   nKtBins1(1),
 188   ktBins(0),
 189   nEPBins(0),
 190   nEPBins1(1),
 191   epBins(0),
 192   nEPBinsMix(0),
 193   nEPBinsMix1(1),
 194   epBinsMix(0),
 195   nCentBins(0),
 196   nCentBins1(1),
 197   centBins(0),
 198   nVzBins(0),
 199   nVzBins1(1),
 200   vzBins(0),
 201   hq(0x0),
 202   hqmix(0x0),
 203   hqinv(0x0),
 204   nqPercBinsLHC11h(1),
 205   qPercBinsLHC11h(0),
 206   hpx(0x0),
 207   hpy(0x0),
 208   hpz(0x0),
 209   hpt(0x0),
 210   hE(0x0),
 211   hphieta(0x0),
 212   hphieta_pid(0x0),
 213   hpt_pid(0x0),
 214   hvzcent(0x0),
 215   hcent(0x0),
 216   hcentUnweighted(0x0),
 217   hcentn(0x0),
 218   hphistaretapair10(0x0),
 219   hphistaretapair16(0x0),
 220   hphistaretapair10a(0x0),
 221   hphistaretapair16a(0x0),
 222   hphistaretapair10b(0x0),
 223   hphistaretapair16b(0x0),
 224   hphietapair(0x0),
 225   hphietapair2(0x0),
 226   hpidid(0x0),
 227   hkt(0x0),
 228   hktcheck(0x0),
 229   hkt3(0x0),
 230   hdcaxy(0x0),
 231   hdcaz(0x0),
 232   hsharequal(0x0),
 233   hsharefrac(0x0),
 234   hsharequalmix(0x0),
 235   hsharefracmix(0x0),
 236   hPsiTPC(0x0),
 237   hPsiV0A(0x0),
 238   hPsiV0C(0x0),
 239   hShiftTPC(0x0),
 240   hShiftV0A(0x0),
 241   hShiftV0C(0x0),
 242   hPsiMix(0x0),
 243   hCheckEPA(0x0),
 244   hCheckEPC(0x0),
 245   hCheckEPmix(0x0),
 246   hAvDphi(0x0),
 247   hNpairs(0x0),
 248   hPairDphi(0x0),
 249   hPairDphiMix(0x0),
 250   hcentq(0x0),
 251   hMixedDistTracks(0x0),
 252   hMixedDistEvents(0x0),
 253   hQvecV0A(0x0),
 254   hQvecV0C(0x0),
 255   hresV0ATPC(0x0),
 256   hresV0CTPC(0x0),
 257   hresV0AV0C(0x0),
 258   hqinvcheck(0x0),
 259   hktbins(0x0),
 260   hcentbins(0x0),
 261   hepbins(0x0)
 262 {
 263
 264   Printf("*******************************************");
 265   Printf("AliAnalysisTaskFemtoESE named %s",name);
 266   Printf("*******************************************");
 267
 268   // default binning
 269   //SetEPBins(12,-TMath::Pi()/12.,2*TMath::Pi()-TMath::Pi()/12.);
 270   SetEPBins(12,0,2*TMath::Pi());
 271   Double_t ktBinsTemp[5] = {0.2,0.3,0.4,0.5,0.7};
 272   SetKtBins(4,ktBinsTemp);
 273   Double_t centBinsTemp[7] = {0,5,10,20,30,40,50};
 274   SetCentBins(6,centBinsTemp);
 275   Double_t vzBinsTemp[9] = {-8,-6,-4,-2,0,2,4,6,8};
 276   SetVzBins(8,vzBinsTemp);
 277
 278   nEPBinsMix = 6;
 279   nEPBinsMix1 = 7;
 280   epBinsMix = new Double_t[nEPBinsMix1];
 281   for(Int_t ee = 0; ee < nEPBinsMix1; ee++)
 282     {epBinsMix[ee] = (TMath::Pi()/(Double_t)nEPBinsMix)*((Double_t)ee)-TMath::Pi()/2.;}
 283
 284   vertex[0] = vertex[1] = vertex[2] = 0.;
 285
 286   DefineInput(0, TChain::Class());
 287   DefineOutput(1, TList::Class());
 288   DefineOutput(2, AliHelperPID::Class());
 289   DefineOutput(3, AliSpectraAODEventCuts::Class());
 290   DefineOutput(4, AliSpectraAODTrackCuts::Class());
 291
 292 }
 293
 294 //________________________________________________________________________
 295 // destructor
 296 AliAnalysisTaskFemtoESE::~AliAnalysisTaskFemtoESE()
 297 {
 298   /* do nothing yet */
 299 }
 300 //________________________________________________________________________
 301 // copy constructor
 302 AliAnalysisTaskFemtoESE::AliAnalysisTaskFemtoESE(const AliAnalysisTaskFemtoESE &/*obj*/) :
 303   AliAnalysisTaskSE(),
 304   fAOD(0x0),
 305   fOutputList(0x0),
 306   fHelperPID(0x0),
 307   fPoolMgr(0x0),
 308   fEventCuts(0x0),
 309   fTrackCuts(0x0),
 310   fFilterBit(128),
 311   fSelectBit(AliVEvent::kMB),
 312   bIsLHC10h(1),
 313   fEventCounter(0),
 314   fMixingTracks(10000),
 315   fBfield(0.),
 316   fMinSepPairEta(0.),
 317   fMinSepPairPhi(0.),
 318   fShareQuality(0.5),
 319   fShareFraction(0.05),
 320   nCountSamePairs(0),
 321   nCountMixedPairs(0),
 322   nCountTracks(0),
 323   fMinQPerc(-1000),
 324   fMaxQPerc(1000),
 325   fQPercDet(0),
 326   fEPDet(0),
 327   nKtBins(0),
 328   nKtBins1(1),
 329   ktBins(0),
 330   nEPBins(0),
 331   nEPBins1(1),
 332   epBins(0),
 333   nEPBinsMix(0),
 334   nEPBinsMix1(1),
 335   epBinsMix(0),
 336   nCentBins(0),
 337   nCentBins1(1),
 338   centBins(0),
 339   nVzBins(0),
 340   nVzBins1(1),
 341   vzBins(0),
 342   hq(0x0),
 343   hqmix(0x0),
 344   hqinv(0x0),
 345   nqPercBinsLHC11h(1),
 346   qPercBinsLHC11h(0),
 347   hpx(0x0),
 348   hpy(0x0),
 349   hpz(0x0),
 350   hpt(0x0),
 351   hE(0x0),
 352   hphieta(0x0),
 353   hphieta_pid(0x0),
 354   hpt_pid(0x0),
 355   hvzcent(0x0),
 356   hcent(0x0),
 357   hcentUnweighted(0x0),
 358   hcentn(0x0),
 359   hphistaretapair10(0x0),
 360   hphistaretapair16(0x0),
 361   hphistaretapair10a(0x0),
 362   hphistaretapair16a(0x0),
 363   hphistaretapair10b(0x0),
 364   hphistaretapair16b(0x0),
 365   hphietapair(0x0),
 366   hphietapair2(0x0),
 367   hpidid(0x0),
 368   hkt(0x0),
 369   hktcheck(0x0),
 370   hkt3(0x0),
 371   hdcaxy(0x0),
 372   hdcaz(0x0),
 373   hsharequal(0x0),
 374   hsharefrac(0x0),
 375   hsharequalmix(0x0),
 376   hsharefracmix(0x0),
 377   hPsiTPC(0x0),
 378   hPsiV0A(0x0),
 379   hPsiV0C(0x0),
 380   hShiftTPC(0x0),
 381   hShiftV0A(0x0),
 382   hShiftV0C(0x0),
 383   hPsiMix(0x0),
 384   hCheckEPA(0x0),
 385   hCheckEPC(0x0),
 386   hCheckEPmix(0x0),
 387   hAvDphi(0x0),
 388   hNpairs(0x0),
 389   hPairDphi(0x0),
 390   hPairDphiMix(0x0),
 391   hcentq(0x0),
 392   hMixedDistTracks(0x0),
 393   hMixedDistEvents(0x0),
 394   hQvecV0A(0x0),
 395   hQvecV0C(0x0),
 396   hresV0ATPC(0x0),
 397   hresV0CTPC(0x0),
 398   hresV0AV0C(0x0),
 399   hqinvcheck(0x0),
 400   hktbins(0x0),
 401   hcentbins(0x0),
 402   hepbins(0x0)
 403 {
 404   /* do nothing yet */
 405 }
 406 //________________________________________________________________________
 407 // assignment operator
 408 AliAnalysisTaskFemtoESE& AliAnalysisTaskFemtoESE::operator=(const AliAnalysisTaskFemtoESE &/*obj*/)
 409 {
 410   /* do nothing yet */
 411   return *this;
 412 }
 413
 414 //________________________________________________________________________
 415 void AliAnalysisTaskFemtoESE::UserCreateOutputObjects()
 416 {
 417   // Create histograms
 418
 419   if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");
 420   if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");
 421   if (!fHelperPID)  AliFatal("HelperPID should be set in the steering macro");
 422
 423   fOutputList = new TList();
 424   fOutputList->SetOwner();
 425   fOutputList->SetName("fOutputList");
 426
 427   hpx = new TH1D("hpx","px",200,-2,2);
 428   hpx->GetXaxis()->SetTitle("p_{x}");
 429   fOutputList->Add(hpx);
 430   hpy = new TH1D("hpy","py",200,-2,2);
 431   hpy->GetXaxis()->SetTitle("p_{y}");
 432   fOutputList->Add(hpy);
 433   hpz = new TH1D("hpz","pz",200,-2,2);
 434   hpz->GetXaxis()->SetTitle("p_{z}");
 435   fOutputList->Add(hpz);
 436   hpt = new TH1D("hpt","pt",100,0,2);
 437   hpt->GetXaxis()->SetTitle("p_{t}");
 438   fOutputList->Add(hpt);
 439   hE = new TH1D("hE","E",100,0,2);
 440   hE->GetXaxis()->SetTitle("E");
 441   fOutputList->Add(hE);
 442   hphieta = new TH2D("hphieta","track #varphi vs #eta",100,0,2*TMath::Pi(),80,-0.8,0.8);
 443   hphieta->GetXaxis()->SetTitle("#varphi");
 444   hphieta->GetYaxis()->SetTitle("#eta");
 445   fOutputList->Add(hphieta);
 446   hphieta_pid = new TH2D("hphieta_pid","PID check -- #Delta#varphi vs #Delta#eta",100,-0.3,0.3,100,-0.3,0.3);
 447   hphieta_pid->GetXaxis()->SetTitle("#Delta#varphi");
 448   hphieta_pid->GetYaxis()->SetTitle("#Delta#eta");
 449   fOutputList->Add(hphieta_pid);
 450   hpt_pid = new TH1D("hpt_pid","PID check -- #Delta p_{t}",100,-0.5,0.5);
 451   hpt_pid->GetXaxis()->SetTitle("#Delta p_{t}");
 452   fOutputList->Add(hpt_pid);
 453   hvzcent = new TH2D("hvzcent","vz vs cent",nVzBins,vzBins,nCentBins,centBins);
 454   hvzcent->GetXaxis()->SetTitle("v_{z}");
 455   hvzcent->GetYaxis()->SetTitle("centrality");
 456   fOutputList->Add(hvzcent);
 457   hcent = new TH1D("hcent","cent",200,0,50);
 458   hcent->GetXaxis()->SetTitle("centrality");
 459   fOutputList->Add(hcent);
 460   hcentUnweighted = new TH1D("hcentUnweighted","cent - unweighted",200,0,50);
 461   hcentUnweighted->GetXaxis()->SetTitle("centrality");
 462   fOutputList->Add(hcentUnweighted);
 463   hcentn = new TH2D("hcentn","cent vs npions",50,0,50,100,0,2000);
 464   hcentn->GetXaxis()->SetTitle("Centrality");
 465   hcentn->GetYaxis()->SetTitle("Number of pions");
 466   fOutputList->Add(hcentn);
 467   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);
 468   hphistaretapair10->GetXaxis()->SetTitle("#Delta#varphi*");
 469   hphistaretapair10->GetYaxis()->SetTitle("#Delta#eta");
 470   hphistaretapair10->GetZaxis()->SetTitle("k_{T}");
 471   fOutputList->Add(hphistaretapair10);
 472   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);
 473   hphistaretapair16->GetXaxis()->SetTitle("#Delta#varphi*");
 474   hphistaretapair16->GetYaxis()->SetTitle("#Delta#eta");
 475   hphistaretapair16->GetZaxis()->SetTitle("k_{T}");
 476   fOutputList->Add(hphistaretapair16);
 477   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);
 478   hphistaretapair10a->GetXaxis()->SetTitle("#Delta#varphi*");
 479   hphistaretapair10a->GetYaxis()->SetTitle("#Delta#eta");
 480   hphistaretapair10a->GetZaxis()->SetTitle("k_{T}");
 481   fOutputList->Add(hphistaretapair10a);
 482   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);
 483   hphistaretapair16a->GetXaxis()->SetTitle("#Delta#varphi*");
 484   hphistaretapair16a->GetYaxis()->SetTitle("#Delta#eta");
 485   hphistaretapair16a->GetZaxis()->SetTitle("k_{T}");
 486   fOutputList->Add(hphistaretapair16a);
 487   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);
 488   hphistaretapair10b->GetXaxis()->SetTitle("#Delta#varphi*");
 489   hphistaretapair10b->GetYaxis()->SetTitle("#Delta#eta");
 490   hphistaretapair10b->GetZaxis()->SetTitle("k_{T}");
 491   fOutputList->Add(hphistaretapair10b);
 492   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);
 493   hphistaretapair16b->GetXaxis()->SetTitle("#Delta#varphi*");
 494   hphistaretapair16b->GetYaxis()->SetTitle("#Delta#eta");
 495   hphistaretapair16b->GetZaxis()->SetTitle("k_{T}");
 496   fOutputList->Add(hphistaretapair16b);
 497   hphietapair = new TH3D("hphietapair","pair #Delta#varphi vs #Delta#eta",100,-0.1,0.1,100,-0.1,0.1,10,0,1);
 498   hphietapair->GetXaxis()->SetTitle("#Delta#varphi");
 499   hphietapair->GetYaxis()->SetTitle("#Delta#eta");
 500   hphietapair->GetZaxis()->SetTitle("k_{T}");
 501   fOutputList->Add(hphietapair);
 502   hphietapair2 = new TH3D("hphietapair2","pair #varphi vs #eta",100,-TMath::Pi(),TMath::Pi(),100,-1.6,1.6,10,0,1);
 503   hphietapair2->GetXaxis()->SetTitle("#Delta#varphi");
 504   hphietapair2->GetYaxis()->SetTitle("#eta");
 505   hphietapair2->GetZaxis()->SetTitle("k_{T}");
 506   fOutputList->Add(hphietapair2);
 507   hpidid = new TH1D("hpidid","pid id",9,-4.5,4.5);
 508   hpidid->GetXaxis()->SetTitle("track PID ID");
 509   fOutputList->Add(hpidid);
 510   hkt = new TH2D("hkt","k_{T}",nCentBins,centBins,100,0,2);
 511   hkt->GetXaxis()->SetTitle("centrality");
 512   hkt->GetYaxis()->SetTitle("k_{T}");
 513   fOutputList->Add(hkt);
 514   hktcheck = new TH1D("hktcheck","k_{T} check",50,0,1);
 515   hktcheck->GetXaxis()->SetTitle("k_{T}");
 516   fOutputList->Add(hktcheck);
 517   hkt3 = new TH3D("hkt3","kt vs pt",50,0,1,50,0,5,50,0,5);
 518   hkt3->GetXaxis()->SetTitle("k_{T}");
 519   hkt3->GetYaxis()->SetTitle("p_{T,1}");
 520   hkt3->GetZaxis()->SetTitle("p_{T,2}");
 521   fOutputList->Add(hkt3);
 522   hdcaxy = new TH2D("hdcaxy","DCA xy",100,-5,5,100,-5,5);
 523   hdcaxy->GetXaxis()->SetTitle("DCA x");
 524   hdcaxy->GetYaxis()->SetTitle("DCA y");
 525   fOutputList->Add(hdcaxy);
 526   hdcaz = new TH1D("hdcaz","DCA z",100,-5,5);
 527   hdcaz->GetXaxis()->SetTitle("DCA z");
 528   fOutputList->Add(hdcaz);
 529   hsharequal = new TH1D("hsharequal","Share Quality",102,-1.02,1.02);
 530   hsharequal->GetXaxis()->SetTitle("Share Quality");
 531   fOutputList->Add(hsharequal);
 532   hsharefrac = new TH1D("hsharefrac","Share Fraction",100,0,1);
 533   hsharefrac->GetXaxis()->SetTitle("Share Fraction");
 534   fOutputList->Add(hsharefrac);
 535   hsharequalmix = new TH1D("hsharequalmix","Share Quality -- mixed events",102,-1.02,1.02);
 536   hsharequalmix->GetXaxis()->SetTitle("Share Quality");
 537   fOutputList->Add(hsharequalmix);
 538   hsharefracmix = new TH1D("hsharefracmix","Share Fraction -- mixed events",100,0,1);
 539   hsharefracmix->GetXaxis()->SetTitle("Share Fraction");
 540   fOutputList->Add(hsharefracmix);
 541   hPsiTPC = new TH2D("hPsiTPC","TPC EP",nCentBins,centBins,100,-1*TMath::Pi(),TMath::Pi());
 542   hPsiTPC->GetXaxis()->SetTitle("Centrality");
 543   hPsiTPC->GetYaxis()->SetTitle("#Psi{TPC}");
 544   fOutputList->Add(hPsiTPC);
 545   hPsiV0A = new TH2D("hPsiV0A","V0A EP",nCentBins,centBins,100,-1*TMath::Pi(),TMath::Pi());
 546   hPsiV0A->GetXaxis()->SetTitle("Centrality");
 547   hPsiV0A->GetYaxis()->SetTitle("#Psi{V0A}");
 548   fOutputList->Add(hPsiV0A);
 549   hPsiV0C = new TH2D("hPsiV0C","V0C EP",nCentBins,centBins,100,-1*TMath::Pi(),TMath::Pi());
 550   hPsiV0C->GetXaxis()->SetTitle("Centrality");
 551   hPsiV0C->GetYaxis()->SetTitle("#Psi{V0C}");
 552   fOutputList->Add(hPsiV0C);
 553   hShiftTPC = new TH2D("hShiftTPC","TPC shifting values",nCentBins,centBins,6,-0.5,5.5);
 554   hShiftTPC->GetXaxis()->SetTitle("Centrality");
 555   fOutputList->Add(hShiftTPC);
 556   hShiftV0A = new TH2D("hShiftV0A","V0A shifting values",nCentBins,centBins,6,-0.5,5.5);
 557   hShiftV0A->GetXaxis()->SetTitle("Centrality");
 558   fOutputList->Add(hShiftV0A);
 559   hShiftV0C = new TH2D("hShiftV0C","V0C shifting values",nCentBins,centBins,6,-0.5,5.5);
 560   hShiftV0C->GetXaxis()->SetTitle("Centrality");
 561   fOutputList->Add(hShiftV0C);
 562   hPsiMix = new TH2D("hPsiMix","Mix EP",nCentBins,centBins,100,-1*TMath::Pi(),TMath::Pi());
 563   hPsiMix->GetXaxis()->SetTitle("Centrality");
 564   hPsiMix->GetYaxis()->SetTitle("#Psi{Mix}");
 565   fOutputList->Add(hPsiMix);
 566   hCheckEPA = new TH2D("hCheckEPA","Check EP V0A",nCentBins,centBins,100,-1*TMath::Pi(),TMath::Pi());
 567   hCheckEPA->GetXaxis()->SetTitle("Centrality");
 568   hCheckEPA->GetYaxis()->SetTitle("PsiV0A - PsiTPC");
 569   fOutputList->Add(hCheckEPA);
 570   hCheckEPC = new TH2D("hCheckEPC","Check EP V0C",nCentBins,centBins,100,-1*TMath::Pi(),TMath::Pi());
 571   hCheckEPC->GetXaxis()->SetTitle("Centrality");
 572   hCheckEPC->GetYaxis()->SetTitle("PsiV0C - PsiTPC");
 573   fOutputList->Add(hCheckEPC);
 574   hCheckEPmix = new TH2D("hCheckEPmix","Check EP mixed events",100,-1*TMath::Pi(),TMath::Pi(),100,-1*TMath::Pi(),TMath::Pi());
 575   hCheckEPmix->GetXaxis()->SetTitle("Psi1 - Psi_mix");
 576   hCheckEPmix->GetYaxis()->SetTitle("Psi1 - Psi2");
 577   fOutputList->Add(hCheckEPmix);
 578   hAvDphi = new TH3D("hAvDphi","average event plane angle",nEPBins,epBins,nCentBins,centBins,nKtBins,ktBins);
 579   hAvDphi->GetXaxis()->SetTitle("EP bins");
 580   hAvDphi->GetYaxis()->SetTitle("cent bins");
 581   hAvDphi->GetZaxis()->SetTitle("kt bins");
 582   fOutputList->Add(hAvDphi);
 583   hNpairs = new TH3D("hNpairs","number of pairs",nEPBins,epBins,nCentBins,centBins,nKtBins,ktBins);
 584   hNpairs->GetXaxis()->SetTitle("EP bins");
 585   hNpairs->GetYaxis()->SetTitle("cent bins");
 586   hNpairs->GetZaxis()->SetTitle("kt bins");
 587   fOutputList->Add(hNpairs);
 588   hPairDphi = new TH1D("hPairDphi","pairs wrt EP",100,0,2*TMath::Pi());
 589   hPairDphi->GetXaxis()->SetTitle("#Delta#phi");
 590   fOutputList->Add(hPairDphi);
 591   hPairDphiMix = new TH1D("hPairDphiMix","mixed pairs wrt EP",100,0,2*TMath::Pi());
 592   hPairDphiMix->GetXaxis()->SetTitle("#Delta#phi");
 593   fOutputList->Add(hPairDphiMix);
 594   hcentq = new TH2D("hcentq","qvec vs cent",100,0,100,50,0,50);
 595   hcentq->GetXaxis()->SetTitle("q_{2} percentile");
 596   hcentq->GetYaxis()->SetTitle("centrality");
 597   fOutputList->Add(hcentq);
 598   hMixedDistTracks = new TH2D("hMixedDistTracks", ";centrality;tracks", nCentBins, centBins, 200, 0, fMixingTracks * 2.);
 599   fOutputList->Add(hMixedDistTracks);
 600   hMixedDistEvents = new TH2D("hMixedDistEvents", ";centrality;events", nCentBins, centBins, 41, -0.5, 40.5);
 601   fOutputList->Add(hMixedDistEvents);
 602   hQvecV0A = new TH2D("hQvecV0A","Qvector in V0A",50,0,50,200,0,5);
 603   hQvecV0A->GetXaxis()->SetTitle("Centrality");
 604   hQvecV0A->GetYaxis()->SetTitle("normalized Qvector");
 605   fOutputList->Add(hQvecV0A);
 606   hQvecV0C = new TH2D("hQvecV0C","Qvector in V0C",50,0,50,200,0,5);
 607   hQvecV0C->GetXaxis()->SetTitle("Centrality");
 608   hQvecV0C->GetYaxis()->SetTitle("normalized Qvector");
 609   fOutputList->Add(hQvecV0C);
 610
 611   // resolution histograms
 612   hresV0ATPC = new TH1D("hresV0ATPC","cent vs cos(2*(V0A-TPC))",nCentBins,centBins);
 613   hresV0ATPC->GetXaxis()->SetTitle("centrality");
 614   fOutputList->Add(hresV0ATPC);
 615   hresV0CTPC = new TH1D("hresV0CTPC","cent vs cos(2*(V0C-TPC))",nCentBins,centBins);
 616   hresV0CTPC->GetXaxis()->SetTitle("centrality");
 617   fOutputList->Add(hresV0CTPC);
 618   hresV0AV0C = new TH1D("hresV0AV0C","cent vs cos(2*(V0A-V0C))",nCentBins,centBins);
 619   hresV0AV0C->GetXaxis()->SetTitle("centrality");
 620   fOutputList->Add(hresV0AV0C);
 621
 622   hqinvcheck = new TH3F("hqinvcheck","Qinv vs kt vs cent",100,0,1,50,0,1,10,0,50);
 623   hqinvcheck->GetXaxis()->SetTitle("qinv");
 624   hqinvcheck->GetYaxis()->SetTitle("kt");
 625   hqinvcheck->GetZaxis()->SetTitle("centrality");
 626   fOutputList->Add(hqinvcheck);
 627
 628   hq = new TH3F****[2];
 629   hqmix = new TH3F****[2];
 630   hqinv = new TH3F****[2];
 631   for(Int_t z = 0; z < 2; z++) // charge combinations
 632     {
 633       hq[z] = new TH3F***[nKtBins];
 634       hqmix[z] = new TH3F***[nKtBins];
 635       hqinv[z] = new TH3F***[nKtBins];
 636       for(Int_t k = 0; k < nKtBins; k++)
 637         {
 638           hq[z][k] = new TH3F**[nEPBins];
 639           hqmix[z][k] = new TH3F**[nEPBins];
 640           hqinv[z][k] = new TH3F**[nEPBins];
 641           for(Int_t e = 0; e < nEPBins; e++)
 642             {
 643               hq[z][k][e] = new TH3F*[nCentBins];
 644               hqmix[z][k][e] = new TH3F*[nCentBins];
 645               hqinv[z][k][e] = new TH3F*[nCentBins];
 646               for(Int_t c = 0; c < nCentBins; c++)
 647                 {
 648                   hq[z][k][e][c] = new TH3F(Form("hq%i_k%i_e%i_c%i",z,k,e,c),Form("hq%i_k%i_e%i_c%i",z,k,e,c),40,-0.2,0.2,40,-0.2,0.2,40,-0.2,0.2);
 649                   hq[z][k][e][c]->Sumw2();
 650                   fOutputList->Add(hq[z][k][e][c]);
 651                   hqmix[z][k][e][c] = new TH3F(Form("hqmix%i_k%i_e%i_c%i",z,k,e,c),Form("hqmix%i_k%i_e%i_c%i",z,k,e,c),40,-0.2,0.2,40,-0.2,0.2,40,-0.2,0.2);
 652                   hqmix[z][k][e][c]->Sumw2();
 653                   fOutputList->Add(hqmix[z][k][e][c]);
 654                   hqinv[z][k][e][c] = new TH3F(Form("hqinv%i_k%i_e%i_c%i",z,k,e,c),Form("hqinv%i_k%i_e%i_c%i",z,k,e,c),40,-0.2,0.2,40,-0.2,0.2,40,-0.2,0.2);
 655                   //hqinv[z][k][e][c]->Sumw2();
 656                   fOutputList->Add(hqinv[z][k][e][c]);
 657                 }
 658             }
 659         }
 660     }
 661
 662   // create dummy histograms which just hold the values of the kt, cent, ep bin edges
 663   hktbins = new TH1F("hktbins","kt bins",nKtBins,ktBins);
 664   fOutputList->Add(hktbins);
 665   hcentbins = new TH1F("hcentbins","cent bins",nCentBins,centBins);
 666   fOutputList->Add(hcentbins);
 667   hepbins = new TH1F("hepbins","ep bins",nEPBins,epBins);
 668   fOutputList->Add(hepbins);
 669
 670   Printf("************************");
 671   Printf("using the %s detector for event plane determination!",fEPDet ? "V0C" : "V0A");
 672   Printf("using the %s detector for q-vector determination!",fQPercDet ? "V0C" : "V0A");
 673   Printf("************************");
 674
 675   vertex[0] = vertex[1] = vertex[2] = 0.;
 676
 677   // event mixing pool
 678   fPoolMgr = new AliEventPoolManager*[2];
 679   Int_t poolsize = 1000;
 680   fPoolMgr[0] = new AliEventPoolManager(poolsize, fMixingTracks, nCentBins, centBins, nVzBins, vzBins,nEPBinsMix,epBinsMix);
 681   fPoolMgr[0]->SetTargetValues(fMixingTracks, 0.01, 5); // check these values
 682   fPoolMgr[1] = new AliEventPoolManager(poolsize, fMixingTracks, nCentBins, centBins, nVzBins, vzBins,nEPBinsMix,epBinsMix);
 683   fPoolMgr[1]->SetTargetValues(fMixingTracks, 0.01, 5); // check these values
 684
 685   nCountSamePairs = 0;
 686   nCountMixedPairs = 0;
 687   nCountTracks = 0;
 688
 689   //stopwatch = new TStopwatch();
 690   //stopwatch->Start();
 691
 692   PostData(1, fOutputList);
 693   PostData(2, fHelperPID);
 694   PostData(3, fEventCuts);
 695   PostData(4, fTrackCuts);
 696
 697 }
 698
 699 //________________________________________________________________________
 700 void AliAnalysisTaskFemtoESE::UserExec(Option_t *)
 701 {
 702   // Main loop
 703   // Called for each event
 704
 705   //if(!fAODcase) {cout<<"ESDs not supported"<<endl; return;}
 706
 707   fAOD = dynamic_cast<AliAODEvent*> (InputEvent());
 708   if (!fAOD) {Printf("ERROR: fAOD not available"); return;}
 709
 710   //if(!EventCut(fAOD)) return;
 711   if(!EventCut()) return;
 712
 713   fEventCounter++;
 714   if(fEventCounter%1000==0) Printf("===========  Event # %i  ===========",fEventCounter);
 715
 716   AliCentrality *centrality;// for AODs and ESDs
 717   const AliAODVertex *primaryVertexAOD;
 718
 719   //AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 720   //AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 721   //fPIDResponse = inputHandler->GetPIDResponse();
 722
 723   fBfield = fAOD->GetMagneticField();
 724
 725   /////////////////////////////////////////////////
 726
 727   Float_t centralityPercentile=0;
 728
 729   centrality = fAOD->GetCentrality();
 730   centralityPercentile = centrality->GetCentralityPercentile("V0M");
 731   if(centralityPercentile <= centBins[0]) return;
 732   if(centralityPercentile > centBins[nCentBins]) return;
 733   Double_t centWeight = GetCentralityWeight(centralityPercentile);
 734
 735   primaryVertexAOD = fAOD->GetPrimaryVertex();
 736   vertex[0]=primaryVertexAOD->GetX(); vertex[1]=primaryVertexAOD->GetY(); vertex[2]=primaryVertexAOD->GetZ();
 737   Double_t zvtx = vertex[2];
 738   if(zvtx < vzBins[0] || zvtx > vzBins[nVzBins]) return; // Z-Vertex Cut
 739   //cout<<"Centrality % = " << centralityPercentile << "  z-vertex = " << zvtx << endl;
 740
 741   //stopwatch->Stop();
 742   //Printf("%lf   %lf",stopwatch->RealTime(),stopwatch->CpuTime());
 743   //stopwatch->Start();
 744
 745   // get event plane from V0's
 746   if(!fEventCuts->IsSelected(fAOD,fTrackCuts)) {Printf("Error! Event not accepted by AliAODSpectraEventCuts!"); return;}
 747   Double_t psiV0A = fEventCuts->GetPsiV0A();
 748   Double_t psiV0C = fEventCuts->GetPsiV0C();
 749   Double_t qperc = -999;
 750   qperc = fEventCuts->GetQvecPercentile(fQPercDet);//0: VZERO-A 1: VZERO-C // now works for both LHC10h and LHC11h (for 0-50% centrality only!)
 751   if(!bIsLHC10h)
 752     {
 753       //Printf("%lf   %lf",qperc,GetQPercLHC11h(fEventCuts->GetqV0A()));
 754       hQvecV0A->Fill(centralityPercentile,fEventCuts->GetqV0A(),centWeight);
 755       hQvecV0C->Fill(centralityPercentile,fEventCuts->GetqV0C(),centWeight);
 756       //Printf("q vector = %lf  percentile = %lf",(fQPercDet==0 ? fEventCuts->GetqV0A() : fEventCuts->GetqV0C()),qperc);
 757     }
 758
 759   if(psiV0A == -999) return;
 760   if(psiV0C == -999) return;
 761   if(qperc < fMinQPerc || qperc > fMaxQPerc) return;
 762
 763   //ProcInfo_t procInfo;
 764   //gSystem->GetProcInfo(&procInfo);
 765   //Printf("beginning of event: ResMem %ld VMem %ld", procInfo.fMemResident, procInfo.fMemVirtual);
 766
 767   Double_t psiEP = psiV0A;
 768   if(fEPDet==1) psiEP = psiV0C;
 769
 770   Double_t sin2phi = 0, cos2phi = 0;
 771
 772   TObjArray* tracksPos = new TObjArray();
 773   tracksPos->SetOwner(kTRUE);
 774   TObjArray* tracksNeg = new TObjArray();
 775   tracksNeg->SetOwner(kTRUE);
 776
 777   // Track loop -- select pions
 778   for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) {
 779     AliAODTrack* aodtrack = (AliAODTrack*)fAOD->GetTrack(i);
 780     if (!aodtrack) continue;
 781     if(!TrackCut(aodtrack)) continue;
 782
 783     // filter bit 7 PID method...
 784     Int_t trackPID=999;
 785     for(Int_t m = 0; m < fAOD->GetNumberOfTracks(); m++) {
 786       AliAODTrack* aodtrack2 = (AliAODTrack*)fAOD->GetTrack(m);
 787       if (!aodtrack2) continue;
 788       if(aodtrack->GetID() != (-aodtrack2->GetID() - 1)) continue;// (-aodTrack2->GetID() - 1)
 789       trackPID=fHelperPID->GetParticleSpecies((AliVTrack*)aodtrack2,kTRUE);
 790       hphieta_pid->Fill(aodtrack->Phi()-aodtrack2->Phi(),aodtrack->Eta()-aodtrack2->Eta());
 791       hpt_pid->Fill(aodtrack->Pt()-aodtrack2->Pt());
 792       //cout << aodtrack->Phi() << "   " << aodtrack->Eta() << "   " << aodtrack->Pt() << endl;
 793       //cout << aodtrack2->Phi() << "   " << aodtrack2->Eta() << "   " << aodtrack2->Pt() << "   " << dataID2 << endl;
 794     }
 795
 796     hpidid->Fill((trackPID+1)*aodtrack->Charge());
 797
 798     // select pions
 799     if(trackPID==0)
 800       {
 801         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());
 802         particle->SetPsiEP(psiEP);
 803         particle->SetTPCClusterMap(aodtrack->GetTPCClusterMap());
 804         particle->SetTPCSharedMap(aodtrack->GetTPCSharedMap());
 805
 806         if(particle->Charge()>0)
 807           tracksPos->Add(particle);
 808         if(particle->Charge()<0)
 809           tracksNeg->Add(particle);
 810
 811         // track qa plots
 812         hpx->Fill(particle->Px());
 813         hpy->Fill(particle->Py());
 814         hpz->Fill(particle->Pz());
 815         hpt->Fill(particle->Pt());
 816         hE->Fill(particle->E());
 817         hphieta->Fill(particle->Phi(),particle->Eta());
 818       }
 819     // check event plane angle using tracks in the TPC
 820     if(aodtrack->Pt() < 2 && aodtrack->Pt() > 0.2)
 821       {
 822         sin2phi += (aodtrack->Pt())*sin(2*aodtrack->Phi());
 823         cos2phi += (aodtrack->Pt())*cos(2*aodtrack->Phi());
 824       }
 825
 826   }
 827   // end track loop
 828
 829   Int_t ntracks = tracksPos->GetEntriesFast()+tracksNeg->GetEntriesFast();
 830
 831   // get EP from TPC, just to check
 832   Double_t psiTPC = 0.;
 833   if(ntracks > 0)
 834     psiTPC = 0.5*atan2(sin2phi,cos2phi);
 835   else return;
 836
 837   hPsiTPC->Fill(centralityPercentile,psiTPC,centWeight);
 838   hPsiV0A->Fill(centralityPercentile,psiV0A,centWeight);
 839   hPsiV0C->Fill(centralityPercentile,psiV0C,centWeight);
 840   Double_t dphiEP = psiTPC-psiV0A;
 841   if(dphiEP>TMath::Pi()) dphiEP-=2*TMath::Pi();
 842   if(dphiEP<-TMath::Pi()) dphiEP+=2*TMath::Pi();
 843   hCheckEPA->Fill(centralityPercentile,dphiEP,centWeight);
 844   dphiEP = psiTPC-psiV0C;
 845   if(dphiEP>TMath::Pi()) dphiEP-=2*TMath::Pi();
 846   if(dphiEP<-TMath::Pi()) dphiEP+=2*TMath::Pi();
 847   hCheckEPC->Fill(centralityPercentile,dphiEP,centWeight);
 848
 849   // values for EP shifting method
 850   hShiftTPC->Fill(centralityPercentile,0.,cos(2*psiTPC)*centWeight);
 851   hShiftTPC->Fill(centralityPercentile,1.,sin(2*psiTPC)*centWeight);
 852   hShiftTPC->Fill(centralityPercentile,2.,cos(4*psiTPC)*centWeight);
 853   hShiftTPC->Fill(centralityPercentile,3.,sin(4*psiTPC)*centWeight);
 854   hShiftTPC->Fill(centralityPercentile,4.,cos(6*psiTPC)*centWeight);
 855   hShiftTPC->Fill(centralityPercentile,5.,sin(6*psiTPC)*centWeight);
 856
 857   hShiftV0A->Fill(centralityPercentile,0.,cos(2*psiV0A)*centWeight);
 858   hShiftV0A->Fill(centralityPercentile,1.,sin(2*psiV0A)*centWeight);
 859   hShiftV0A->Fill(centralityPercentile,2.,cos(4*psiV0A)*centWeight);
 860   hShiftV0A->Fill(centralityPercentile,3.,sin(4*psiV0A)*centWeight);
 861   hShiftV0A->Fill(centralityPercentile,4.,cos(6*psiV0A)*centWeight);
 862   hShiftV0A->Fill(centralityPercentile,5.,sin(6*psiV0A)*centWeight);
 863
 864   hShiftV0C->Fill(centralityPercentile,0.,cos(2*psiV0C)*centWeight);
 865   hShiftV0C->Fill(centralityPercentile,1.,sin(2*psiV0C)*centWeight);
 866   hShiftV0C->Fill(centralityPercentile,2.,cos(4*psiV0C)*centWeight);
 867   hShiftV0C->Fill(centralityPercentile,3.,sin(4*psiV0C)*centWeight);
 868   hShiftV0C->Fill(centralityPercentile,4.,cos(6*psiV0C)*centWeight);
 869   hShiftV0C->Fill(centralityPercentile,5.,sin(6*psiV0C)*centWeight);
 870
 871
 872   hcentq->Fill(qperc,centralityPercentile,centWeight);
 873
 874   nCountTracks += ntracks;
 875   //cout << "Found " << ntracks << " pion tracks..." << endl;
 876
 877   hvzcent->Fill(zvtx,centralityPercentile,centWeight);
 878   hcentUnweighted->Fill(centralityPercentile);
 879   hcent->Fill(centralityPercentile,centWeight);
 880   hcentn->Fill(centralityPercentile,ntracks,centWeight);
 881
 882   // resolution histograms
 883   hresV0ATPC->Fill(centralityPercentile,cos(2*(psiV0A-psiTPC))*centWeight);
 884   hresV0CTPC->Fill(centralityPercentile,cos(2*(psiV0C-psiTPC))*centWeight);
 885   hresV0AV0C->Fill(centralityPercentile,cos(2*(psiV0A-psiV0C))*centWeight);
 886
 887   AliEventPool* poolPos = fPoolMgr[0]->GetEventPool(centralityPercentile,zvtx,psiEP);
 888   if (!poolPos) AliFatal(Form("No pool found for centrality = %f, vz = %f, ep = %f", centralityPercentile, zvtx,psiEP));
 889   //Printf("Positive pool found for centrality = %f, vz = %f, ep = %f with %i events in it", centralityPercentile, zvtx,psiEP,poolPos->GetCurrentNEvents());
 890   AliEventPool* poolNeg = fPoolMgr[1]->GetEventPool(centralityPercentile,zvtx,psiEP);
 891   if (!poolNeg) AliFatal(Form("No pool found for centrality = %f, vz = %f, ep = %f", centralityPercentile, zvtx,psiEP));
 892   //Printf("Negative pool found for centrality = %f, vz = %f, ep = %f with %i events in it", centralityPercentile, zvtx,psiEP,poolNeg->GetCurrentNEvents());
 893
 894   TrackLoop(tracksPos,poolPos,0,psiEP,centralityPercentile);
 895   TrackLoop(tracksNeg,poolNeg,1,psiEP,centralityPercentile);
 896   //TObjArray* clonedtracks = CloneAndReduceTrackList(tracks,psiEP);
 897   poolPos->UpdatePool(tracksPos);
 898   poolNeg->UpdatePool(tracksNeg);
 899   //cout << "pool contains " << pool->GetCurrentNEvents() << " events and " << pool->NTracksInPool() << " tracks." << endl;
 900   //tracks->Clear();
 901
 902   //delete tracks;
 903
 904   //gSystem->GetProcInfo(&procInfo);
 905   //Printf("end of event: ResMem %ld VMem %ld", procInfo.fMemResident, procInfo.fMemVirtual);
 906
 907   // Post output data.
 908   PostData(1, fOutputList);
 909   PostData(2, fHelperPID);
 910   PostData(3, fEventCuts);
 911   PostData(4, fTrackCuts);
 912 }
 913
 914
 915 void AliAnalysisTaskFemtoESE::TrackLoop(TObjArray *tracks, AliEventPool *pool, Int_t z, Double_t psiEP, Float_t centralityPercentile)
 916 {
 917   // z=0 for positive charges, z=1 for negative charges
 918   Double_t kt = 0;
 919   Double_t qout=0, qside=0, qlong=0;
 920   Double_t pVect1[4] = {0,0,0,0};
 921   Double_t pVect2[4] = {0,0,0,0};
 922   Int_t k, e, c; //bin indices for histograms
 923
 924   Double_t centWeight = GetCentralityWeight(centralityPercentile);
 925
 926   Int_t ntracks = tracks->GetEntriesFast();
 927
 928   Int_t countMixedEvents = 0, countMixedTracks = 0;
 929
 930   // same event loop
 931   for(Int_t j = 0; j < ntracks; j++)
 932     {
 933       //cout << endl << j << "   ";
 934       AliFemtoESEBasicParticle* track1 = (AliFemtoESEBasicParticle*)tracks->At(j);
 935       pVect1[0]=track1->E();
 936       pVect1[1]=track1->Px();
 937       pVect1[2]=track1->Py();
 938       pVect1[3]=track1->Pz();
 939       //cout << pVect1[0] << "   " << pVect1[1] << "   " <<  pVect1[2] << "   " << pVect1[3] << endl;
 940
 941       for(Int_t i = j+1; i < ntracks; i++)
 942         {
 943           AliFemtoESEBasicParticle* track2 = (AliFemtoESEBasicParticle*)tracks->At(i);
 944
 945           Double_t deltaphistar10 = DeltaPhiStar(track1,track2,1.0);
 946           Double_t deltaphistar16 = DeltaPhiStar(track1,track2,1.6);
 947
 948           hphistaretapair10->Fill(deltaphistar10,track1->Eta()-track2->Eta(),kt);
 949           hphistaretapair16->Fill(deltaphistar16,track1->Eta()-track2->Eta(),kt);
 950
 951           if(!PairCut(track1,track2,kFALSE)) continue;
 952
 953           hphistaretapair10a->Fill(deltaphistar10,track1->Eta()-track2->Eta(),kt);
 954           hphistaretapair16a->Fill(deltaphistar16,track1->Eta()-track2->Eta(),kt);
 955           hphistaretapair10b->Fill(deltaphistar10,track1->Eta()-track2->Eta(),kt);
 956           hphistaretapair16b->Fill(deltaphistar16,track1->Eta()-track2->Eta(),kt);
 957
 958           pVect2[0]=track2->E();
 959           pVect2[1]=track2->Px();
 960           pVect2[2]=track2->Py();
 961           pVect2[3]=track2->Pz();
 962
 963           //qinv = GetQinv(pVect1, pVect2); // = qinv**2 = (P1x-P2x)**2 + (P1y-P2y)**2 + (P1z-P2z)**2 - (P1t-P2t)**2
 964           GetQosl(pVect1, pVect2, qout, qside, qlong); // qout, qside, qlong = components of Q=P1-P2 in the P=P1+P2 frame
 965           kt = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.; // = Kt = |pT1+pT2|/2
 966           hkt->Fill(centralityPercentile,kt,centWeight);
 967           hkt3->Fill(kt,track1->Pt(),track2->Pt());
 968           Double_t deltaphi = GetDeltaPhiEP(pVect1[1],pVect1[2],pVect2[1],pVect2[2],psiEP); // angle to event plane in correct range
 969           if(fabs(qout)<0.2 && fabs(qside)<0.2 && fabs(qlong)<0.2) nCountSamePairs++;
 970           if(kt < ktBins[0] || kt > ktBins[nKtBins]) continue;
 971           if(!FindBin(kt,deltaphi,centralityPercentile,k,e,c)) continue;
 972           hktcheck->Fill(kt);
 973           hq[z][k][e][c]->Fill(qout,qside,qlong,centWeight);
 974           hqinv[z][k][e][c]->Fill(qout,qside,qlong,GetQinv(pVect1, pVect2)*centWeight);
 975           hNpairs->Fill(deltaphi,centralityPercentile,kt);
 976           hAvDphi->Fill(deltaphi,centralityPercentile,kt,deltaphi);
 977           hPairDphi->Fill(deltaphi);
 978           hqinvcheck->Fill(GetQinv(pVect1, pVect2),kt,centralityPercentile,centWeight);
 979           Double_t dphi = track1->Phi()-track2->Phi();
 980           if(dphi<-TMath::Pi()) dphi += 2*TMath::Pi();
 981           if(dphi>TMath::Pi()) dphi -= 2*TMath::Pi();
 982           hphietapair->Fill(dphi,track1->Eta()-track2->Eta(),kt);
 983           hphietapair2->Fill(dphi,track1->Eta()-track2->Eta(),kt);
 984           //cout << k << "  ";
 985         }
 986     }
 987
 988   // mixed event loop
 989   countMixedTracks = 0;
 990   countMixedEvents = 0;
 991   if (pool->IsReady())
 992     {
 993       Int_t nMix = pool->GetCurrentNEvents();
 994
 995       for (Int_t jMix=0; jMix<nMix; jMix++) // loop over events in pool
 996         {
 997           TObjArray* bgTracks = pool->GetEvent(jMix);
 998           Int_t ntracksmix = bgTracks->GetEntriesFast();
 999           if(ntracksmix <= 0) continue;
1000
1001           // compare event planes of two events
1002           AliFemtoESEBasicParticle* tracktest = (AliFemtoESEBasicParticle*)bgTracks->UncheckedAt(0);
1003           Double_t psiEPmixtemp = tracktest->GetPsiEP();
1004           /*Double_t dphiEPtest = fPsiEPmixtemp-psiEP;
1005           while(dphiEPtest>2*TMath::Pi()) dphiEPtest-=2*TMath::Pi();
1006           while(dphiEPtest<0) dphiEPtest+=2*TMath::Pi();
1007           if(dphiEPtest>TMath::Pi()) dphiEPtest-=TMath::Pi();
1008           if(dphiEPtest>TMath::Pi()/2.) dphiEPtest = TMath::Pi()-dphiEPtest;
1009           if(dphiEPtest > TMath::Pi()/6.) continue; // event planes must be within pi/6
1010           */
1011
1012           Double_t psiEPmix = 0.5*atan2(sin(2*psiEP)+sin(2*psiEPmixtemp),cos(2*psiEP)+cos(2*psiEPmixtemp));
1013           hPsiMix->Fill(centralityPercentile,psiEPmix,centWeight);
1014
1015           Double_t dphimix = psiEP-psiEPmix;
1016           if(dphimix < -TMath::Pi()) dphimix += 2*TMath::Pi();
1017           if(dphimix > TMath::Pi()) dphimix -= 2*TMath::Pi();
1018           Double_t dphi12 = psiEP-psiEPmixtemp;
1019           if(dphi12 < -TMath::Pi()) dphi12 += 2*TMath::Pi();
1020           if(dphi12 > TMath::Pi()) dphi12 -= 2*TMath::Pi();
1021           hCheckEPmix->Fill(dphimix,dphi12);
1022
1023           countMixedTracks += ntracksmix;
1024           countMixedEvents += 1;
1025
1026           for(Int_t j = 0; j < ntracks; j++)
1027             {
1028               //cout << endl << j << "   ";
1029               AliFemtoESEBasicParticle* track1 = (AliFemtoESEBasicParticle*)tracks->At(j);
1030               pVect1[0]=track1->E();
1031               pVect1[1]=track1->Px();
1032               pVect1[2]=track1->Py();
1033               pVect1[3]=track1->Pz();
1034
1035               for(Int_t i = 0; i < ntracksmix; i++)
1036                 {
1037                   AliFemtoESEBasicParticle* track2 = (AliFemtoESEBasicParticle*)bgTracks->UncheckedAt(i);
1038
1039                   if(!PairCut(track1,track2,kTRUE)) continue;
1040
1041                   pVect2[0]=track2->E();
1042                   pVect2[1]=track2->Px();
1043                   pVect2[2]=track2->Py();
1044                   pVect2[3]=track2->Pz();
1045
1046                   if(psiEPmixtemp != track2->GetPsiEP()) AliFatal("Error! Event plane angles are wrong in mixing!!");
1047
1048                   //qinv = GetQinv(pVect1, pVect2); // qinv**2 = (P1x-P2x)**2 + (P1y-P2y)**2 + (P1z-P2z)**2 - (P1t-P2t)**2
1049                   GetQosl(pVect1, pVect2, qout, qside, qlong); // qout, qside, qlong = components of Q=P1-P2 in the P=P1+P2 frame
1050                   kt = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.; // = Kt = |pT1+pT2|/2
1051                   Double_t deltaphi = GetDeltaPhiEP(pVect1[1],pVect1[2],pVect2[1],pVect2[2],psiEPmix); // angle to event plane in correct range
1052
1053                   //Double_t weight = 1./(Double_t)nMix;
1054                   if(fabs(qout)<0.2 && fabs(qside)<0.2 && fabs(qlong)<0.2) nCountMixedPairs++;
1055                   if(kt < ktBins[0] || kt > ktBins[nKtBins]) continue;
1056                   if(!FindBin(kt,deltaphi,centralityPercentile,k,e,c)) continue;
1057                   hqmix[z][k][e][c]->Fill(qout,qside,qlong,centWeight);
1058
1059                   hPairDphiMix->Fill(deltaphi);
1060
1061                 }
1062             }
1063         }
1064     }
1065
1066   hMixedDistTracks->Fill(centralityPercentile, countMixedTracks);
1067   hMixedDistEvents->Fill(centralityPercentile, countMixedEvents);
1068   //Printf("mixed tracks: %i   mixed events: %i",countMixedTracks,countMixedEvents);
1069 }
1070
1071
1072
1073
1074
1075 //________________________________________________________________________
1076 void AliAnalysisTaskFemtoESE::Terminate(Option_t *)
1077 {
1078
1079   if(ktBins) delete [] ktBins;
1080   if(epBins) delete [] epBins;
1081   if(centBins) delete [] centBins;
1082   if(vzBins) delete [] vzBins;
1083   if(qPercBinsLHC11h) delete [] qPercBinsLHC11h;
1084
1085   // Called once at the end of the query
1086
1087   Printf("Done");
1088
1089 }
1090
1091
1092 /*
1093
1094 //________________________________________________________________________
1095 Bool_t AliAnalysisTaskFemtoESE::AcceptPair(AliChaoticityTrackStruct *first, AliChaoticityTrackStruct *second)
1096 {
1097
1098 if(fabs(first->fEta-second->fEta) > fMinSepPairEta) return kTRUE;
1099
1100 // propagate through B field to r=1m
1101 Float_t phi1 = first->fPhi - asin(first->fCharge*(0.1*fBfield)*0.15/first->fPt);// 0.15 for D=1m
1102 if(phi1 > 2*PI) phi1 -= 2*PI;
1103 if(phi1 < 0) phi1 += 2*PI;
1104 Float_t phi2 = second->fPhi - asin(second->fCharge*(0.1*fBfield)*0.15/second->fPt);// 0.15 for D=1m
1105 if(phi2 > 2*PI) phi2 -= 2*PI;
1106 if(phi2 < 0) phi2 += 2*PI;
1107
1108 Float_t deltaphi = phi1 - phi2;
1109 if(deltaphi > PI) deltaphi -= 2*PI;
1110 if(deltaphi < -PI) deltaphi += 2*PI;
1111 deltaphi = fabs(deltaphi);
1112
1113 if(deltaphi < fMinSepPairPhi) return kFALSE;// Min Separation
1114
1115
1116 // propagate through B field to r=1.6m
1117 phi1 = first->fPhi - asin(first->fCharge*(0.1*fBfield)*0.24/first->fPt);// mine. 0.24 for D=1.6m
1118 if(phi1 > 2*PI) phi1 -= 2*PI;
1119 if(phi1 < 0) phi1 += 2*PI;
1120 phi2 = second->fPhi - asin(second->fCharge*(0.1*fBfield)*0.24/second->fPt);// mine. 0.24 for D=1.6m
1121 if(phi2 > 2*PI) phi2 -= 2*PI;
1122 if(phi2 < 0) phi2 += 2*PI;
1123
1124 deltaphi = phi1 - phi2;
1125 if(deltaphi > PI) deltaphi -= 2*PI;
1126 if(deltaphi < -PI) deltaphi += 2*PI;
1127 deltaphi = fabs(deltaphi);
1128
1129 if(deltaphi < fMinSepPairPhi) return kFALSE;// Min Separation
1130
1131
1132
1133 //
1134
1135 Int_t ncl1 = first->fClusterMap.GetNbits();
1136 Int_t ncl2 = second->fClusterMap.GetNbits();
1137 Int_t sumCls = 0; Int_t sumSha = 0; Int_t sumQ = 0;
1138 Double_t shfrac = 0; Double_t qfactor = 0;
1139 for(Int_t imap = 0; imap < ncl1 && imap < ncl2; imap++) {
1140 if (first->fClusterMap.TestBitNumber(imap) && second->fClusterMap.TestBitNumber(imap)) {// Both clusters
1141 if (first->fSharedMap.TestBitNumber(imap) && second->fSharedMap.TestBitNumber(imap)) { // Shared
1142 sumQ++;
1143 sumCls+=2;
1144 sumSha+=2;}
1145 else {sumQ--; sumCls+=2;}
1146 }
1147 else if (first->fClusterMap.TestBitNumber(imap) || second->fClusterMap.TestBitNumber(imap)) {// Non shared
1148 sumQ++;
1149 sumCls++;}
1150 }
1151 if (sumCls>0) {
1152 qfactor = sumQ*1.0/sumCls;
1153 shfrac = sumSha*1.0/sumCls;
1154 }
1155
1156 if(qfactor > fShareQuality || shfrac > fShareFraction) return kFALSE;
1157
1158
1159 return kTRUE;
1160
1161
1162 }
1163 //________________________________________________________________________
1164 Float_t AliAnalysisTaskFemtoESE::GamovFactor(Int_t chargeBin1, Int_t chargeBin2, Float_t qinv)
1165 {
1166 Float_t arg = G_Coeff/qinv;
1167
1168 if(chargeBin1==chargeBin2) return (exp(arg)-1)/(arg);
1169 else {return (exp(-arg)-1)/(-arg);}
1170
1171 }
1172 //________________________________________________________________________
1173 void AliAnalysisTaskFemtoESE::Shuffle(Int_t *iarr, Int_t i1, Int_t i2)
1174 {
1175 Int_t j, k;
1176 Int_t a = i2 - i1;
1177 for (Int_t i = i1; i < i2+1; i++) {
1178 j = (Int_t) (gRandom->Rndm() * a);
1179 k = iarr[j];
1180 iarr[j] = iarr[i];
1181 iarr[i] = k;
1182 }
1183 }
1184 */
1185 //________________________________________________________________________
1186 Double_t AliAnalysisTaskFemtoESE::GetQinv(Double_t track1[], Double_t track2[]){
1187
1188   Double_t qinv2=1.0;
1189   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.);
1190   if(qinv2 >= 0.) return sqrt(qinv2);
1191   else return -1.*sqrt(-1.*qinv2);
1192 }
1193
1194 //________________________________________________________________________
1195 void AliAnalysisTaskFemtoESE::GetQosl(Double_t track1[], Double_t track2[], Double_t& qout, Double_t& qside, Double_t& qlong){
1196
1197   Double_t p0 = track1[0] + track2[0];
1198   Double_t px = track1[1] + track2[1];
1199   Double_t py = track1[2] + track2[2];
1200   Double_t pz = track1[3] + track2[3];
1201
1202   Double_t mt = sqrt(p0*p0 - pz*pz);
1203   Double_t pt = sqrt(px*px + py*py);
1204
1205   Double_t v0 = track1[0] - track2[0];
1206   Double_t vx = track1[1] - track2[1];
1207   Double_t vy = track1[2] - track2[2];
1208   Double_t vz = track1[3] - track2[3];
1209
1210   if(gRandom->Rndm()<0.5)
1211     {
1212       v0 = -v0;
1213       vx = -vx;
1214       vy = -vy;
1215       vz = -vz;
1216     }
1217
1218   //cout << p0 << "  " << px << "  " << py << "  " << pz << "  " << v0 << "  " << vx << "  " << vy << "  " << vz << "  " << mt << "  " << pt << endl;
1219
1220   qout = (px*vx + py*vy)/pt;
1221   qside = (px*vy - py*vx)/pt;
1222   qlong = (p0*vz - pz*v0)/mt;
1223 }
1224
1225 //________________________________________________________________________
1226 Bool_t AliAnalysisTaskFemtoESE::EventCut(/*AliAODEvent* fevent*/){
1227
1228   // Trigger Cut
1229
1230   Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
1231   Bool_t isSelected2 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral);
1232   Bool_t isSelected3 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kSemiCentral);
1233   if(!isSelected1 && !isSelected2 && !isSelected3) {return kFALSE;}
1234
1235   /*
1236   // Pile-up rejection
1237   AliAnalysisUtils *AnaUtil=new AliAnalysisUtils();
1238   if(!fPbPbcase) AnaUtil->SetUseMVPlpSelection(kTRUE);// use Multi-Vertex tool for pp and pPb
1239   else AnaUtil->SetUseMVPlpSelection(kFALSE);
1240   Bool_t pileUpCase=AnaUtil->IsPileUpEvent(fAOD);
1241   if(pileUpCase) return;
1242
1243   // Vertexing
1244   primaryVertexAOD = fAOD->GetPrimaryVertex();
1245   vertex[0]=primaryVertexAOD->GetX(); vertex[1]=primaryVertexAOD->GetY(); vertex[2]=primaryVertexAOD->GetZ();
1246
1247   if(fabs(vertex[2]) > 10) {cout<<"Zvertex Out of Range. Skip Event"<<endl; return;} // Z-Vertex Cut
1248   ((TH3D*)fOutputList->FindObject("fVertexDist"))->Fill(vertex[0], vertex[1], vertex[2]);
1249
1250   for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) {
1251   AliAODTrack* aodtrack = fAOD->GetTrack(i);
1252   if (!aodtrack) continue;
1253   AODTracks++;
1254   if(!aodtrack->TestFilterBit(BIT(fFilterBit))) continue;// AOD filterBit cut
1255   FBTracks++;
1256   }
1257   ((TH1D*)fOutputList->FindObject("fMultDist2"))->Fill(FBTracks);
1258
1259   //if(fAOD->IsPileupFromSPD()) {cout<<"PileUpEvent. Skip Event"<<endl; return;} // Old Pile-up cut
1260   if(primaryVertexAOD->GetNContributors() < 1) {cout<<"Bad Vertex. Skip Event"<<endl; return;}
1261
1262   ((TH1D*)fOutputList->FindObject("fMultDist3"))->Fill(FBTracks);
1263
1264   fBfield = fAOD->GetMagneticField();
1265
1266   for(Int_t i=0; i<fZvertexBins; i++){
1267   if( (vertex[2] >= zstart+i*zstep) && (vertex[2] < zstart+(i+1)*zstep) ){
1268   zbin=i;
1269   break;
1270   }
1271   }
1272   */
1273
1274   return kTRUE;
1275
1276 }
1277
1278 //________________________________________________________________________
1279 Bool_t AliAnalysisTaskFemtoESE::TrackCut(AliAODTrack* ftrack){
1280
1281   if (!ftrack->TestFilterBit(fFilterBit)) return kFALSE;
1282
1283   if(ftrack->Pt() < 0.14) return kFALSE;
1284   if(ftrack->Pt() > 1.5) return kFALSE;
1285   if(fabs(ftrack->Eta()) > 0.8) return kFALSE;
1286
1287   if(ftrack->GetTPCNcls() < 80) return kFALSE;// TPC nCluster cut
1288
1289   Double_t trackdca[3] = {ftrack->XAtDCA(),ftrack->YAtDCA(),ftrack->ZAtDCA()};
1290   //ftrack->XYZAtDCA(trackdca);
1291   //Double_t dcaxy = sqrt( pow(trackpos[0] - vertex[0],2) + pow(trackpos[1] - vertex[1],2));
1292   //Double_t dcaz = sqrt( pow(trackpos[2] - vertex[2],2));
1293   hdcaxy->Fill(trackdca[0],trackdca[1]);
1294   hdcaz->Fill(trackdca[2]);
1295   //if(dcaxy > 0.2) return kFALSE;
1296   //if(dcaz > 0.15) return kFALSE;
1297
1298
1299   //// FilterBit Overlap Check
1300   //if(fFilterBit != 7){
1301   //    Bool_t goodTrackOtherFB = kFALSE;
1302   //    if(fMCcase && fAOD->GetRunNumber()<=126437) goodTrackOtherFB=kTRUE;// FB7 to FB5 mapping in 10f6a MC does not work
1303   //
1304   //    for (Int_t j = 0; j < fAOD->GetNumberOfTracks(); j++) {
1305   //      AliAODTrack* aodtrack2 = fAOD->GetTrack(randomIndex[j]);
1306   //      if(!aodtrack2) continue;
1307   //      if(!aodtrack2->TestFilterBit(BIT(fFilterBit))) continue;
1308   //
1309   //      if(-(aodtrack->GetID()+1)==aodtrack2->GetID()) {goodTrackOtherFB=kTRUE; break;}
1310   //
1311   //    }
1312   //    if(!goodTrackOtherFB) continue;
1313   //    }
1314
1315   return kTRUE;
1316
1317 }
1318
1319 //________________________________________________________________________
1320 Bool_t AliAnalysisTaskFemtoESE::PairCut(AliFemtoESEBasicParticle* ftrack1, AliFemtoESEBasicParticle* ftrack2, Bool_t mix){
1321
1322   // check for same charge
1323   if(ftrack2->Charge() != ftrack1->Charge()) return kFALSE;
1324
1325   // qinv cut
1326   Double_t trackvec1[4] = {ftrack1->E(),ftrack1->Px(),ftrack1->Py(),ftrack1->Pz()};
1327   Double_t trackvec2[4] = {ftrack2->E(),ftrack2->Px(),ftrack2->Py(),ftrack2->Pz()};
1328   Double_t qinv = GetQinv(trackvec1,trackvec2);
1329   if(qinv < 0.005) return kFALSE; // qinv < 0.005
1330
1331   // deltaEta x deltaPhi* cut
1332   if(fabs(ftrack1->Eta()-ftrack2->Eta()) < fMinSepPairEta)
1333     {
1334       Double_t deltaphistar = DeltaPhiStar(ftrack1,ftrack2,1.0); // angular separation at r=1m
1335       deltaphistar = fabs(deltaphistar);
1336       if(deltaphistar < fMinSepPairPhi) return kFALSE;
1337       deltaphistar = DeltaPhiStar(ftrack1,ftrack2,1.6); // angular separation at r=1.6m
1338       deltaphistar = fabs(deltaphistar);
1339       if(deltaphistar < fMinSepPairPhi) return kFALSE;
1340     }
1341
1342   // share fraction & share quality cut
1343   TBits clusterMap1 = (TBits)(ftrack1->GetTPCClusterMap());
1344   TBits sharedMap1 = (TBits)(ftrack1->GetTPCSharedMap());
1345   TBits clusterMap2 = (TBits)(ftrack2->GetTPCClusterMap());
1346   TBits sharedMap2 = (TBits)(ftrack2->GetTPCSharedMap());
1347
1348   Int_t ncl1 = clusterMap1.GetNbits();
1349   Int_t ncl2 = clusterMap2.GetNbits();
1350   Int_t sumCls = 0; Int_t sumSha = 0; Int_t sumQ = 0;
1351   Double_t shfrac = 0; Double_t qfactor = 0;
1352   for(Int_t imap = 0; imap < ncl1 && imap < ncl2; imap++)
1353     {
1354       if (clusterMap1.TestBitNumber(imap) && clusterMap2.TestBitNumber(imap)) // Both clusters
1355         {
1356           if (sharedMap1.TestBitNumber(imap) && sharedMap2.TestBitNumber(imap)) // Shared
1357             {
1358               sumQ++;
1359               sumCls+=2;
1360               sumSha+=2;}
1361           else {sumQ--; sumCls+=2;}
1362         }
1363       else if (clusterMap1.TestBitNumber(imap) || clusterMap2.TestBitNumber(imap)) // Non shared
1364         {
1365           sumQ++;
1366           sumCls++;
1367         }
1368     }
1369   if (sumCls>0)
1370     {
1371       qfactor = sumQ*1.0/sumCls;
1372       shfrac = sumSha*1.0/sumCls;
1373     }
1374
1375   // sumCls -- number of clusters in track 1 + number of clusters in track 2 (clusters in both tracks counted twice)
1376   // sumSha -- number of shared clusters (counted twice)
1377   // sumQ -- ?
1378
1379   if(!mix)
1380     {
1381       hsharequal->Fill(qfactor);
1382       hsharefrac->Fill(shfrac);
1383     }
1384   else
1385     {
1386       hsharequalmix->Fill(qfactor);
1387       hsharefracmix->Fill(shfrac);
1388     }
1389
1390
1391   if(qfactor > fShareQuality || shfrac > fShareFraction) return kFALSE;
1392
1393   return kTRUE;
1394 }
1395
1396 Double_t AliAnalysisTaskFemtoESE::DeltaPhiStar(AliAODTrack* ftrack1, AliAODTrack* ftrack2, Double_t r)
1397 {
1398   Double_t phi1 = ftrack1->Phi() - asin(0.3*ftrack1->Charge()*(0.1*fBfield)*r/(2.*ftrack1->Pt())); // magnetic field converted from kilogauss to tesla
1399   Double_t phi2 = ftrack2->Phi() - asin(0.3*ftrack2->Charge()*(0.1*fBfield)*r/(2.*ftrack2->Pt()));
1400
1401   Double_t deltaphi = phi1 - phi2;
1402   while(deltaphi > TMath::Pi()) deltaphi -= 2*TMath::Pi();
1403   while(deltaphi < -TMath::Pi()) deltaphi += 2*TMath::Pi();
1404
1405   return deltaphi;
1406 }
1407
1408 /*TObjArray* AliAnalysisTaskFemtoESE::CloneAndReduceTrackList(TObjArray* tracks, Double_t psi)
1409   {
1410   // clones a track list by using AliDPhiBasicParticle which uses much less memory (used for event mixing)
1411
1412   TObjArray* tracksClone = new TObjArray;
1413   tracksClone->SetOwner(kTRUE);
1414
1415   for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
1416   {
1417   AliAODTrack* particle = (AliAODTrack*) tracks->UncheckedAt(i);
1418   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());
1419   copy->SetPsiEP(psi);
1420   copy->SetTPCClusterMap(particle->GetTPCClusterMap());
1421   copy->SetTPCSharedMap(particle->GetTPCSharedMap());
1422
1423   tracksClone->Add(copy);
1424   }
1425
1426   return tracksClone;
1427   }*/
1428
1429 Double_t AliAnalysisTaskFemtoESE::GetDeltaPhiEP(Double_t px1, Double_t py1, Double_t px2, Double_t py2, Double_t psi)
1430 {
1431   // angle of pair
1432   Double_t px = px1+px2;
1433   Double_t py = py1+py2;
1434   Double_t phi = atan2(py,px);
1435
1436   Double_t dphi = phi-psi;
1437   while(dphi < epBins[0]) dphi += 2*TMath::Pi();
1438   while(dphi > epBins[nEPBins]) dphi -= 2*TMath::Pi();
1439
1440   return dphi;
1441 }
1442
1443 Bool_t AliAnalysisTaskFemtoESE::FindBin(Double_t kt, Double_t phi, Double_t cent, Int_t& a, Int_t& b, Int_t&c)
1444 {
1445   a = b = c = -1;
1446   for(Int_t i = 0; i < nKtBins; i++)
1447     {
1448       if(kt >= ktBins[i] && kt < ktBins[i+1])
1449         {
1450           a = i;
1451           break;
1452         }
1453     }
1454   if(a==-1) return kFALSE;
1455
1456   for(Int_t i = 0; i < nEPBins; i++)
1457     {
1458       if(phi >= epBins[i] && phi < epBins[i+1])
1459         {
1460           b = i;
1461           break;
1462         }
1463     }
1464   if(b==-1) return kFALSE;
1465
1466   for(Int_t i = 0; i < nCentBins; i++)
1467     {
1468       if(cent >= centBins[i] && cent < centBins[i+1])
1469         {
1470           c = i;
1471           break;
1472         }
1473     }
1474   if(c==-1) return kFALSE;
1475
1476   return kTRUE;
1477 }
1478
1479
1480
1481 void AliAnalysisTaskFemtoESE::SetKtBins(Int_t n, Double_t* bins)
1482 {
1483   if(ktBins) delete [] ktBins;
1484   nKtBins = n;
1485   nKtBins1 = n+1;
1486   ktBins = new Double_t[nKtBins+1];
1487   for(Int_t i = 0; i < nKtBins+1; i++)
1488     ktBins[i]=bins[i];
1489   Printf("Setting %i kt bins: ",nKtBins);
1490   for(Int_t i = 0; i < nKtBins+1; i++) Printf("%lf",ktBins[i]);
1491 }
1492 void AliAnalysisTaskFemtoESE::SetCentBins(Int_t n, Double_t* bins)
1493 {
1494   if(centBins) delete [] centBins;
1495   nCentBins = n;
1496   nCentBins1 = n+1;
1497   centBins = new Double_t[nCentBins+1];
1498   for(Int_t i = 0; i < nCentBins+1; i++)
1499     centBins[i]=bins[i];
1500   Printf("Setting %i centrality bins: ",nCentBins);
1501   for(Int_t i = 0; i < nCentBins+1; i++) Printf("%lf",centBins[i]);
1502 }
1503 void AliAnalysisTaskFemtoESE::SetVzBins(Int_t n, Double_t* bins)
1504 {
1505   if(vzBins) delete [] vzBins;
1506   nVzBins = n;
1507   nVzBins1 = n+1;
1508   vzBins = new Double_t[nVzBins+1];
1509   for(Int_t i = 0; i < nVzBins+1; i++)
1510     vzBins[i]=bins[i];
1511   Printf("Setting %i vz bins: ",nVzBins);
1512   for(Int_t i = 0; i < nVzBins+1; i++) Printf("%lf",vzBins[i]);
1513 }
1514 void AliAnalysisTaskFemtoESE::SetEPBins(Int_t n, Double_t min, Double_t max)
1515 {
1516   if(epBins) delete [] epBins;
1517   nEPBins = n;
1518   nEPBins1 = n+1;
1519   epBins = new Double_t[nEPBins+1];
1520   for(Int_t y = 0; y < nEPBins+1; y++)
1521     epBins[y] = min+((max-min)/(Double_t)n)*((Double_t)y);
1522   Printf("Setting %i EP bins: ",nEPBins);
1523   for(Int_t i = 0; i < nEPBins+1; i++) Printf("%lf",epBins[i]);
1524 }
1525
1526 Double_t AliAnalysisTaskFemtoESE::GetQPercLHC11h(Double_t qvec)
1527 {
1528   // preliminary attemp at calculating qvector percentile in LHC11h -- still very approximate and only works in 5% bins
1529   if(!qPercBinsLHC11h)
1530     {
1531       Double_t tempArray[21] = {0.0, 0.22995, 0.33047, 0.410831, 0.480728, 0.545566, 0.606841, 0.66634, 0.725193, 0.783813, 0.843311, 0.904185, 0.96796, 1.03522, 1.10768, 1.18774, 1.27808, 1.3857, 1.52438, 1.73633, 4.95};
1532       nqPercBinsLHC11h = 21;
1533       qPercBinsLHC11h = new Double_t[nqPercBinsLHC11h];
1534       for(Int_t n = 0; n < nqPercBinsLHC11h; n++) qPercBinsLHC11h[n] = tempArray[n];
1535     }
1536
1537   for(Int_t t = 0; t < nqPercBinsLHC11h-1; t++)
1538     if(qvec > qPercBinsLHC11h[t] && qvec < qPercBinsLHC11h[t+1]) return 50.*(2*t+1)/(Double_t)(nqPercBinsLHC11h-1);
1539
1540   if(qvec < qPercBinsLHC11h[0]) return 0.0;
1541   if(qvec > qPercBinsLHC11h[nqPercBinsLHC11h-1]) return 100.0;
1542
1543   return 0.0;
1544
1545 }
1546
1547 /*Double_t AliAnalysisTaskFemtoESE::GetCentralityWeight(Double_t cent)
1548 {
1549
1550   // use makeCentWeighting.C to fit centrality distribution to obtain this parameterization
1551   Double_t par1 = 2.60629;
1552   Double_t par2 = 0.579333;
1553   Double_t limit1 = 8.71488;
1554   Double_t limit2 = 12.0126;
1555
1556   if(cent < limit1) return 1./par1;
1557   if(cent > limit2) return 1./par2;
1558
1559   Double_t slope = (par2-par1)/(limit2-limit1);
1560   Double_t b = par1-slope*limit1;
1561
1562   return 1./(slope*cent+b);
1563   }*/
1564
1565 Double_t AliAnalysisTaskFemtoESE::GetCentralityWeight(Double_t cent)
1566 {
1567
1568   // use makeCentWeighting.C to fit centrality distribution to obtain this parameterization
1569   Double_t par1 = 2.60629;
1570   Double_t par2 = 0.579333;
1571   Double_t limit1 = 8.71488;
1572   Double_t limit2 = 12.0126;
1573
1574   if(cent < limit1) return 1.;
1575   if(cent > limit2) return 1.;
1576
1577   Double_t slope = (par2-par1)/(limit2-limit1);
1578   Double_t b = par1-slope*limit1;
1579
1580   if(cent < 10.)
1581     return par1/(slope*cent+b);
1582   else
1583     return par2/(slope*cent+b);
1584 }