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 * 1.5);
 599   fOutputList->Add(hMixedDistTracks);
 600   hMixedDistEvents = new TH2D("hMixedDistEvents", ";centrality;events", nCentBins, centBins, 21, -0.5, 20.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***[nKtBins];
 629   hqmix = new TH3F***[nKtBins];
 630   hqinv = new TH3F***[nKtBins];
 631   for(Int_t k = 0; k < nKtBins; k++)
 632     {
 633       hq[k] = new TH3F**[nEPBins];
 634       hqmix[k] = new TH3F**[nEPBins];
 635       hqinv[k] = new TH3F**[nEPBins];
 636       for(Int_t e = 0; e < nEPBins; e++)
 637         {
 638           hq[k][e] = new TH3F*[nCentBins];
 639           hqmix[k][e] = new TH3F*[nCentBins];
 640           hqinv[k][e] = new TH3F*[nCentBins];
 641           for(Int_t c = 0; c < nCentBins; c++)
 642             {
 643               hq[k][e][c] = new TH3F(Form("hq_%i_%i_%i",k,e,c),Form("hq_%i_%i_%i",k,e,c),20,-0.2,0.2,20,-0.2,0.2,20,-0.2,0.2);
 644               fOutputList->Add(hq[k][e][c]);
 645               hqmix[k][e][c] = new TH3F(Form("hqmix_%i_%i_%i",k,e,c),Form("hqmix_%i_%i_%i",k,e,c),20,-0.2,0.2,20,-0.2,0.2,20,-0.2,0.2);
 646               fOutputList->Add(hqmix[k][e][c]);
 647               hqinv[k][e][c] = new TH3F(Form("hqinv_%i_%i_%i",k,e,c),Form("hqinv_%i_%i_%i",k,e,c),20,-0.2,0.2,20,-0.2,0.2,20,-0.2,0.2);
 648               fOutputList->Add(hqinv[k][e][c]);
 649             }
 650         }
 651     }
 652
 653   // create dummy histograms which just hold the values of the kt, cent, ep bin edges
 654   hktbins = new TH1F("hktbins","kt bins",nKtBins,ktBins);
 655   fOutputList->Add(hktbins);
 656   hcentbins = new TH1F("hcentbins","cent bins",nCentBins,centBins);
 657   fOutputList->Add(hcentbins);
 658   hepbins = new TH1F("hepbins","ep bins",nEPBins,epBins);
 659   fOutputList->Add(hepbins);
 660
 661   Printf("************************");
 662   Printf("using the %s detector for event plane determination!",fEPDet ? "V0C" : "V0A");
 663   Printf("using the %s detector for q-vector determination!",fQPercDet ? "V0C" : "V0A");
 664   Printf("************************");
 665
 666   vertex[0] = vertex[1] = vertex[2] = 0.;
 667
 668   // event mixing pool
 669   fPoolMgr = new AliEventPoolManager*[2];
 670   Int_t poolsize = 1000;
 671   fPoolMgr[0] = new AliEventPoolManager(poolsize, fMixingTracks, nCentBins, centBins, nVzBins, vzBins,nEPBinsMix,epBinsMix);
 672   fPoolMgr[0]->SetTargetValues(fMixingTracks, 0.01, 5); // check these values
 673   fPoolMgr[1] = new AliEventPoolManager(poolsize, fMixingTracks, nCentBins, centBins, nVzBins, vzBins,nEPBinsMix,epBinsMix);
 674   fPoolMgr[1]->SetTargetValues(fMixingTracks, 0.01, 5); // check these values
 675
 676   nCountSamePairs = 0;
 677   nCountMixedPairs = 0;
 678   nCountTracks = 0;
 679
 680   //stopwatch = new TStopwatch();
 681   //stopwatch->Start();
 682
 683   PostData(1, fOutputList);
 684   PostData(2, fHelperPID);
 685   PostData(3, fEventCuts);
 686   PostData(4, fTrackCuts);
 687
 688 }
 689
 690 //________________________________________________________________________
 691 void AliAnalysisTaskFemtoESE::UserExec(Option_t *)
 692 {
 693   // Main loop
 694   // Called for each event
 695
 696   //if(!fAODcase) {cout<<"ESDs not supported"<<endl; return;}
 697
 698   fAOD = dynamic_cast<AliAODEvent*> (InputEvent());
 699   if (!fAOD) {Printf("ERROR: fAOD not available"); return;}
 700
 701   //if(!EventCut(fAOD)) return;
 702   if(!EventCut()) return;
 703
 704   fEventCounter++;
 705   if(fEventCounter%1000==0) Printf("===========  Event # %i  ===========",fEventCounter);
 706
 707   AliCentrality *centrality;// for AODs and ESDs
 708   const AliAODVertex *primaryVertexAOD;
 709
 710   //AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 711   //AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 712   //fPIDResponse = inputHandler->GetPIDResponse();
 713
 714   fBfield = fAOD->GetMagneticField();
 715
 716   /////////////////////////////////////////////////
 717
 718   Float_t centralityPercentile=0;
 719
 720   centrality = fAOD->GetCentrality();
 721   centralityPercentile = centrality->GetCentralityPercentile("V0M");
 722   if(centralityPercentile <= centBins[0]) return;
 723   if(centralityPercentile > centBins[nCentBins]) return;
 724   Double_t centWeight = GetCentralityWeight(centralityPercentile);
 725
 726   primaryVertexAOD = fAOD->GetPrimaryVertex();
 727   vertex[0]=primaryVertexAOD->GetX(); vertex[1]=primaryVertexAOD->GetY(); vertex[2]=primaryVertexAOD->GetZ();
 728   Double_t zvtx = vertex[2];
 729   if(zvtx < vzBins[0] || zvtx > vzBins[nVzBins]) return; // Z-Vertex Cut
 730   //cout<<"Centrality % = " << centralityPercentile << "  z-vertex = " << zvtx << endl;
 731
 732   //stopwatch->Stop();
 733   //Printf("%lf   %lf",stopwatch->RealTime(),stopwatch->CpuTime());
 734   //stopwatch->Start();
 735
 736   // get event plane from V0's
 737   if(!fEventCuts->IsSelected(fAOD,fTrackCuts)) {Printf("Error! Event not accepted by AliAODSpectraEventCuts!"); return;}
 738   Double_t psiV0A = fEventCuts->GetPsiV0A();
 739   Double_t psiV0C = fEventCuts->GetPsiV0C();
 740   Double_t qperc = -999;
 741   qperc = fEventCuts->GetQvecPercentile(fQPercDet);//0: VZERO-A 1: VZERO-C // now works for both LHC10h and LHC11h (for 0-50% centrality only!)
 742   if(!bIsLHC10h)
 743     {
 744       //Printf("%lf   %lf",qperc,GetQPercLHC11h(fEventCuts->GetqV0A()));
 745       hQvecV0A->Fill(centralityPercentile,fEventCuts->GetqV0A(),centWeight);
 746       hQvecV0C->Fill(centralityPercentile,fEventCuts->GetqV0C(),centWeight);
 747       //Printf("q vector = %lf  percentile = %lf",(fQPercDet==0 ? fEventCuts->GetqV0A() : fEventCuts->GetqV0C()),qperc);
 748     }
 749
 750   if(psiV0A == -999) return;
 751   if(psiV0C == -999) return;
 752   if(qperc < fMinQPerc || qperc > fMaxQPerc) return;
 753
 754   //ProcInfo_t procInfo;
 755   //gSystem->GetProcInfo(&procInfo);
 756   //Printf("beginning of event: ResMem %ld VMem %ld", procInfo.fMemResident, procInfo.fMemVirtual);
 757
 758   Double_t psiEP = psiV0A;
 759   if(fEPDet==1) psiEP = psiV0C;
 760
 761   Double_t sin2phi = 0, cos2phi = 0;
 762
 763   TObjArray* tracksPos = new TObjArray();
 764   tracksPos->SetOwner(kTRUE);
 765   TObjArray* tracksNeg = new TObjArray();
 766   tracksNeg->SetOwner(kTRUE);
 767
 768   // Track loop -- select pions
 769   for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) {
 770     AliAODTrack* aodtrack = (AliAODTrack*)fAOD->GetTrack(i);
 771     if (!aodtrack) continue;
 772     if(!TrackCut(aodtrack)) continue;
 773
 774     // filter bit 7 PID method...
 775     Int_t trackPID=999;
 776     for(Int_t m = 0; m < fAOD->GetNumberOfTracks(); m++) {
 777       AliAODTrack* aodtrack2 = (AliAODTrack*)fAOD->GetTrack(m);
 778       if (!aodtrack2) continue;
 779       if(aodtrack->GetID() != (-aodtrack2->GetID() - 1)) continue;// (-aodTrack2->GetID() - 1)
 780       trackPID=fHelperPID->GetParticleSpecies((AliVTrack*)aodtrack2,kTRUE);
 781       hphieta_pid->Fill(aodtrack->Phi()-aodtrack2->Phi(),aodtrack->Eta()-aodtrack2->Eta());
 782       hpt_pid->Fill(aodtrack->Pt()-aodtrack2->Pt());
 783       //cout << aodtrack->Phi() << "   " << aodtrack->Eta() << "   " << aodtrack->Pt() << endl;
 784       //cout << aodtrack2->Phi() << "   " << aodtrack2->Eta() << "   " << aodtrack2->Pt() << "   " << dataID2 << endl;
 785     }
 786
 787     hpidid->Fill((trackPID+1)*aodtrack->Charge());
 788
 789     // select pions
 790     if(trackPID==0)
 791       {
 792         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());
 793         particle->SetPsiEP(psiEP);
 794         particle->SetTPCClusterMap(aodtrack->GetTPCClusterMap());
 795         particle->SetTPCSharedMap(aodtrack->GetTPCSharedMap());
 796
 797         if(particle->Charge()>0)
 798           tracksPos->Add(particle);
 799         if(particle->Charge()<0)
 800           tracksNeg->Add(particle);
 801
 802         // track qa plots
 803         hpx->Fill(particle->Px());
 804         hpy->Fill(particle->Py());
 805         hpz->Fill(particle->Pz());
 806         hpt->Fill(particle->Pt());
 807         hE->Fill(particle->E());
 808         hphieta->Fill(particle->Phi(),particle->Eta());
 809       }
 810     // check event plane angle using tracks in the TPC
 811     if(aodtrack->Pt() < 2 && aodtrack->Pt() > 0.2)
 812       {
 813         sin2phi += (aodtrack->Pt())*sin(2*aodtrack->Phi());
 814         cos2phi += (aodtrack->Pt())*cos(2*aodtrack->Phi());
 815       }
 816
 817   }
 818   // end track loop
 819
 820   Int_t ntracks = tracksPos->GetEntriesFast()+tracksNeg->GetEntriesFast();
 821
 822   // get EP from TPC, just to check
 823   Double_t psiTPC = 0.;
 824   if(ntracks > 0)
 825     psiTPC = 0.5*atan2(sin2phi,cos2phi);
 826   else return;
 827
 828   hPsiTPC->Fill(centralityPercentile,psiTPC,centWeight);
 829   hPsiV0A->Fill(centralityPercentile,psiV0A,centWeight);
 830   hPsiV0C->Fill(centralityPercentile,psiV0C,centWeight);
 831   Double_t dphiEP = psiTPC-psiV0A;
 832   if(dphiEP>TMath::Pi()) dphiEP-=2*TMath::Pi();
 833   if(dphiEP<-TMath::Pi()) dphiEP+=2*TMath::Pi();
 834   hCheckEPA->Fill(centralityPercentile,dphiEP,centWeight);
 835   dphiEP = psiTPC-psiV0C;
 836   if(dphiEP>TMath::Pi()) dphiEP-=2*TMath::Pi();
 837   if(dphiEP<-TMath::Pi()) dphiEP+=2*TMath::Pi();
 838   hCheckEPC->Fill(centralityPercentile,dphiEP,centWeight);
 839
 840   // values for EP shifting method
 841   hShiftTPC->Fill(centralityPercentile,0.,cos(2*psiTPC)*centWeight);
 842   hShiftTPC->Fill(centralityPercentile,1.,sin(2*psiTPC)*centWeight);
 843   hShiftTPC->Fill(centralityPercentile,2.,cos(4*psiTPC)*centWeight);
 844   hShiftTPC->Fill(centralityPercentile,3.,sin(4*psiTPC)*centWeight);
 845   hShiftTPC->Fill(centralityPercentile,4.,cos(6*psiTPC)*centWeight);
 846   hShiftTPC->Fill(centralityPercentile,5.,sin(6*psiTPC)*centWeight);
 847
 848   hShiftV0A->Fill(centralityPercentile,0.,cos(2*psiV0A)*centWeight);
 849   hShiftV0A->Fill(centralityPercentile,1.,sin(2*psiV0A)*centWeight);
 850   hShiftV0A->Fill(centralityPercentile,2.,cos(4*psiV0A)*centWeight);
 851   hShiftV0A->Fill(centralityPercentile,3.,sin(4*psiV0A)*centWeight);
 852   hShiftV0A->Fill(centralityPercentile,4.,cos(6*psiV0A)*centWeight);
 853   hShiftV0A->Fill(centralityPercentile,5.,sin(6*psiV0A)*centWeight);
 854
 855   hShiftV0C->Fill(centralityPercentile,0.,cos(2*psiV0C)*centWeight);
 856   hShiftV0C->Fill(centralityPercentile,1.,sin(2*psiV0C)*centWeight);
 857   hShiftV0C->Fill(centralityPercentile,2.,cos(4*psiV0C)*centWeight);
 858   hShiftV0C->Fill(centralityPercentile,3.,sin(4*psiV0C)*centWeight);
 859   hShiftV0C->Fill(centralityPercentile,4.,cos(6*psiV0C)*centWeight);
 860   hShiftV0C->Fill(centralityPercentile,5.,sin(6*psiV0C)*centWeight);
 861
 862
 863   hcentq->Fill(qperc,centralityPercentile,centWeight);
 864
 865   nCountTracks += ntracks;
 866   //cout << "Found " << ntracks << " pion tracks..." << endl;
 867
 868   hvzcent->Fill(zvtx,centralityPercentile,centWeight);
 869   hcentUnweighted->Fill(centralityPercentile);
 870   hcent->Fill(centralityPercentile,centWeight);
 871   hcentn->Fill(centralityPercentile,ntracks,centWeight);
 872
 873   // resolution histograms
 874   hresV0ATPC->Fill(centralityPercentile,cos(2*(psiV0A-psiTPC))*centWeight);
 875   hresV0CTPC->Fill(centralityPercentile,cos(2*(psiV0C-psiTPC))*centWeight);
 876   hresV0AV0C->Fill(centralityPercentile,cos(2*(psiV0A-psiV0C))*centWeight);
 877
 878   AliEventPool* poolPos = fPoolMgr[0]->GetEventPool(centralityPercentile,zvtx,psiEP);
 879   if (!poolPos) AliFatal(Form("No pool found for centrality = %f, vz = %f, ep = %f", centralityPercentile, zvtx,psiEP));
 880   //Printf("Positive pool found for centrality = %f, vz = %f, ep = %f with %i events in it", centralityPercentile, zvtx,psiEP,poolPos->GetCurrentNEvents());
 881   AliEventPool* poolNeg = fPoolMgr[1]->GetEventPool(centralityPercentile,zvtx,psiEP);
 882   if (!poolNeg) AliFatal(Form("No pool found for centrality = %f, vz = %f, ep = %f", centralityPercentile, zvtx,psiEP));
 883   //Printf("Negative pool found for centrality = %f, vz = %f, ep = %f with %i events in it", centralityPercentile, zvtx,psiEP,poolNeg->GetCurrentNEvents());
 884
 885   TrackLoop(tracksPos,poolPos,psiEP,centralityPercentile);
 886   TrackLoop(tracksNeg,poolNeg,psiEP,centralityPercentile);
 887   //TObjArray* clonedtracks = CloneAndReduceTrackList(tracks,psiEP);
 888   poolPos->UpdatePool(tracksPos);
 889   poolNeg->UpdatePool(tracksNeg);
 890   //cout << "pool contains " << pool->GetCurrentNEvents() << " events and " << pool->NTracksInPool() << " tracks." << endl;
 891   //tracks->Clear();
 892
 893   //delete tracks;
 894
 895   //gSystem->GetProcInfo(&procInfo);
 896   //Printf("end of event: ResMem %ld VMem %ld", procInfo.fMemResident, procInfo.fMemVirtual);
 897
 898   // Post output data.
 899   PostData(1, fOutputList);
 900   PostData(2, fHelperPID);
 901   PostData(3, fEventCuts);
 902   PostData(4, fTrackCuts);
 903 }
 904
 905
 906 void AliAnalysisTaskFemtoESE::TrackLoop(TObjArray *tracks, AliEventPool *pool, Double_t psiEP, Float_t centralityPercentile)
 907 {
 908   Double_t kt = 0;
 909   Double_t qout=0, qside=0, qlong=0;
 910   Double_t pVect1[4] = {0,0,0,0};
 911   Double_t pVect2[4] = {0,0,0,0};
 912   Int_t k, e, c; //bin indices for histograms
 913
 914   Double_t centWeight = GetCentralityWeight(centralityPercentile);
 915
 916   Int_t ntracks = tracks->GetEntriesFast();
 917
 918   Int_t countMixedEvents = 0, countMixedTracks = 0;
 919
 920   // same event loop
 921   for(Int_t j = 0; j < ntracks; j++)
 922     {
 923       //cout << endl << j << "   ";
 924       AliFemtoESEBasicParticle* track1 = (AliFemtoESEBasicParticle*)tracks->At(j);
 925       pVect1[0]=track1->E();
 926       pVect1[1]=track1->Px();
 927       pVect1[2]=track1->Py();
 928       pVect1[3]=track1->Pz();
 929       //cout << pVect1[0] << "   " << pVect1[1] << "   " <<  pVect1[2] << "   " << pVect1[3] << endl;
 930
 931       for(Int_t i = j+1; i < ntracks; i++)
 932         {
 933           AliFemtoESEBasicParticle* track2 = (AliFemtoESEBasicParticle*)tracks->At(i);
 934
 935           Double_t deltaphistar10 = DeltaPhiStar(track1,track2,1.0);
 936           Double_t deltaphistar16 = DeltaPhiStar(track1,track2,1.6);
 937
 938           hphistaretapair10->Fill(deltaphistar10,track1->Eta()-track2->Eta(),kt);
 939           hphistaretapair16->Fill(deltaphistar16,track1->Eta()-track2->Eta(),kt);
 940
 941           if(!PairCut(track1,track2,kFALSE)) continue;
 942
 943           hphistaretapair10a->Fill(deltaphistar10,track1->Eta()-track2->Eta(),kt);
 944           hphistaretapair16a->Fill(deltaphistar16,track1->Eta()-track2->Eta(),kt);
 945           hphistaretapair10b->Fill(deltaphistar10,track1->Eta()-track2->Eta(),kt);
 946           hphistaretapair16b->Fill(deltaphistar16,track1->Eta()-track2->Eta(),kt);
 947
 948           pVect2[0]=track2->E();
 949           pVect2[1]=track2->Px();
 950           pVect2[2]=track2->Py();
 951           pVect2[3]=track2->Pz();
 952
 953           //qinv = GetQinv(pVect1, pVect2); // = qinv**2 = (P1x-P2x)**2 + (P1y-P2y)**2 + (P1z-P2z)**2 - (P1t-P2t)**2
 954           GetQosl(pVect1, pVect2, qout, qside, qlong); // qout, qside, qlong = components of Q=P1-P2 in the P=P1+P2 frame
 955           kt = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.; // = Kt = |pT1+pT2|/2
 956           hkt->Fill(centralityPercentile,kt,centWeight);
 957           hkt3->Fill(kt,track1->Pt(),track2->Pt());
 958           Double_t deltaphi = GetDeltaPhiEP(pVect1[1],pVect1[2],pVect2[1],pVect2[2],psiEP); // angle to event plane in correct range
 959           if(fabs(qout)<0.2 && fabs(qside)<0.2 && fabs(qlong)<0.2) nCountSamePairs++;
 960           if(kt < ktBins[0] || kt > ktBins[nKtBins]) continue;
 961           if(!FindBin(kt,deltaphi,centralityPercentile,k,e,c)) continue;
 962           hktcheck->Fill(kt);
 963           hq[k][e][c]->Fill(qout,qside,qlong,centWeight);
 964           hqinv[k][e][c]->Fill(qout,qside,qlong,GetQinv(pVect1, pVect2)*centWeight);
 965           hNpairs->Fill(deltaphi,centralityPercentile,kt);
 966           hAvDphi->Fill(deltaphi,centralityPercentile,kt,deltaphi);
 967           hPairDphi->Fill(deltaphi);
 968           hqinvcheck->Fill(GetQinv(pVect1, pVect2),kt,centralityPercentile,centWeight);
 969           Double_t dphi = track1->Phi()-track2->Phi();
 970           if(dphi<-TMath::Pi()) dphi += 2*TMath::Pi();
 971           if(dphi>TMath::Pi()) dphi -= 2*TMath::Pi();
 972           hphietapair->Fill(dphi,track1->Eta()-track2->Eta(),kt);
 973           hphietapair2->Fill(dphi,track1->Eta()-track2->Eta(),kt);
 974           //cout << k << "  ";
 975         }
 976     }
 977
 978   // mixed event loop
 979   countMixedTracks = 0;
 980   countMixedEvents = 0;
 981   if (pool->IsReady())
 982     {
 983       Int_t nMix = pool->GetCurrentNEvents();
 984
 985       for (Int_t jMix=0; jMix<nMix; jMix++) // loop over events in pool
 986         {
 987           TObjArray* bgTracks = pool->GetEvent(jMix);
 988           Int_t ntracksmix = bgTracks->GetEntriesFast();
 989           if(ntracksmix <= 0) continue;
 990
 991           // compare event planes of two events
 992           AliFemtoESEBasicParticle* tracktest = (AliFemtoESEBasicParticle*)bgTracks->UncheckedAt(0);
 993           Double_t psiEPmixtemp = tracktest->GetPsiEP();
 994           /*Double_t dphiEPtest = fPsiEPmixtemp-psiEP;
 995           while(dphiEPtest>2*TMath::Pi()) dphiEPtest-=2*TMath::Pi();
 996           while(dphiEPtest<0) dphiEPtest+=2*TMath::Pi();
 997           if(dphiEPtest>TMath::Pi()) dphiEPtest-=TMath::Pi();
 998           if(dphiEPtest>TMath::Pi()/2.) dphiEPtest = TMath::Pi()-dphiEPtest;
 999           if(dphiEPtest > TMath::Pi()/6.) continue; // event planes must be within pi/6
1000           */
1001
1002           Double_t psiEPmix = 0.5*atan2(sin(2*psiEP)+sin(2*psiEPmixtemp),cos(2*psiEP)+cos(2*psiEPmixtemp));
1003           hPsiMix->Fill(centralityPercentile,psiEPmix,centWeight);
1004
1005           Double_t dphimix = psiEP-psiEPmix;
1006           if(dphimix < -TMath::Pi()) dphimix += 2*TMath::Pi();
1007           if(dphimix > TMath::Pi()) dphimix -= 2*TMath::Pi();
1008           Double_t dphi12 = psiEP-psiEPmixtemp;
1009           if(dphi12 < -TMath::Pi()) dphi12 += 2*TMath::Pi();
1010           if(dphi12 > TMath::Pi()) dphi12 -= 2*TMath::Pi();
1011           hCheckEPmix->Fill(dphimix,dphi12);
1012
1013           countMixedTracks += ntracksmix;
1014           countMixedEvents += 1;
1015
1016           for(Int_t j = 0; j < ntracks; j++)
1017             {
1018               //cout << endl << j << "   ";
1019               AliFemtoESEBasicParticle* track1 = (AliFemtoESEBasicParticle*)tracks->At(j);
1020               pVect1[0]=track1->E();
1021               pVect1[1]=track1->Px();
1022               pVect1[2]=track1->Py();
1023               pVect1[3]=track1->Pz();
1024
1025               for(Int_t i = 0; i < ntracksmix; i++)
1026                 {
1027                   AliFemtoESEBasicParticle* track2 = (AliFemtoESEBasicParticle*)bgTracks->UncheckedAt(i);
1028
1029                   if(!PairCut(track1,track2,kTRUE)) continue;
1030
1031                   pVect2[0]=track2->E();
1032                   pVect2[1]=track2->Px();
1033                   pVect2[2]=track2->Py();
1034                   pVect2[3]=track2->Pz();
1035
1036                   if(psiEPmixtemp != track2->GetPsiEP()) AliFatal("Error! Event plane angles are wrong in mixing!!");
1037
1038                   //qinv = GetQinv(pVect1, pVect2); // qinv**2 = (P1x-P2x)**2 + (P1y-P2y)**2 + (P1z-P2z)**2 - (P1t-P2t)**2
1039                   GetQosl(pVect1, pVect2, qout, qside, qlong); // qout, qside, qlong = components of Q=P1-P2 in the P=P1+P2 frame
1040                   kt = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.; // = Kt = |pT1+pT2|/2
1041                   Double_t deltaphi = GetDeltaPhiEP(pVect1[1],pVect1[2],pVect2[1],pVect2[2],psiEPmix); // angle to event plane in correct range
1042
1043                   //Double_t weight = 1./(Double_t)nMix;
1044                   if(fabs(qout)<0.2 && fabs(qside)<0.2 && fabs(qlong)<0.2) nCountMixedPairs++;
1045                   if(kt < ktBins[0] || kt > ktBins[nKtBins]) continue;
1046                   if(!FindBin(kt,deltaphi,centralityPercentile,k,e,c)) continue;
1047                   hqmix[k][e][c]->Fill(qout,qside,qlong,centWeight);
1048
1049                   hPairDphiMix->Fill(deltaphi);
1050
1051                 }
1052             }
1053         }
1054     }
1055
1056   hMixedDistTracks->Fill(centralityPercentile, countMixedTracks);
1057   hMixedDistEvents->Fill(centralityPercentile, countMixedEvents);
1058   //Printf("mixed tracks: %i   mixed events: %i",countMixedTracks,countMixedEvents);
1059 }
1060
1061
1062
1063
1064
1065 //________________________________________________________________________
1066 void AliAnalysisTaskFemtoESE::Terminate(Option_t *)
1067 {
1068
1069   if(ktBins) delete [] ktBins;
1070   if(epBins) delete [] epBins;
1071   if(centBins) delete [] centBins;
1072   if(vzBins) delete [] vzBins;
1073   if(qPercBinsLHC11h) delete [] qPercBinsLHC11h;
1074
1075   // Called once at the end of the query
1076
1077   Printf("Done");
1078
1079 }
1080
1081
1082 /*
1083
1084 //________________________________________________________________________
1085 Bool_t AliAnalysisTaskFemtoESE::AcceptPair(AliChaoticityTrackStruct *first, AliChaoticityTrackStruct *second)
1086 {
1087
1088 if(fabs(first->fEta-second->fEta) > fMinSepPairEta) return kTRUE;
1089
1090 // propagate through B field to r=1m
1091 Float_t phi1 = first->fPhi - asin(first->fCharge*(0.1*fBfield)*0.15/first->fPt);// 0.15 for D=1m
1092 if(phi1 > 2*PI) phi1 -= 2*PI;
1093 if(phi1 < 0) phi1 += 2*PI;
1094 Float_t phi2 = second->fPhi - asin(second->fCharge*(0.1*fBfield)*0.15/second->fPt);// 0.15 for D=1m
1095 if(phi2 > 2*PI) phi2 -= 2*PI;
1096 if(phi2 < 0) phi2 += 2*PI;
1097
1098 Float_t deltaphi = phi1 - phi2;
1099 if(deltaphi > PI) deltaphi -= 2*PI;
1100 if(deltaphi < -PI) deltaphi += 2*PI;
1101 deltaphi = fabs(deltaphi);
1102
1103 if(deltaphi < fMinSepPairPhi) return kFALSE;// Min Separation
1104
1105
1106 // propagate through B field to r=1.6m
1107 phi1 = first->fPhi - asin(first->fCharge*(0.1*fBfield)*0.24/first->fPt);// mine. 0.24 for D=1.6m
1108 if(phi1 > 2*PI) phi1 -= 2*PI;
1109 if(phi1 < 0) phi1 += 2*PI;
1110 phi2 = second->fPhi - asin(second->fCharge*(0.1*fBfield)*0.24/second->fPt);// mine. 0.24 for D=1.6m
1111 if(phi2 > 2*PI) phi2 -= 2*PI;
1112 if(phi2 < 0) phi2 += 2*PI;
1113
1114 deltaphi = phi1 - phi2;
1115 if(deltaphi > PI) deltaphi -= 2*PI;
1116 if(deltaphi < -PI) deltaphi += 2*PI;
1117 deltaphi = fabs(deltaphi);
1118
1119 if(deltaphi < fMinSepPairPhi) return kFALSE;// Min Separation
1120
1121
1122
1123 //
1124
1125 Int_t ncl1 = first->fClusterMap.GetNbits();
1126 Int_t ncl2 = second->fClusterMap.GetNbits();
1127 Int_t sumCls = 0; Int_t sumSha = 0; Int_t sumQ = 0;
1128 Double_t shfrac = 0; Double_t qfactor = 0;
1129 for(Int_t imap = 0; imap < ncl1 && imap < ncl2; imap++) {
1130 if (first->fClusterMap.TestBitNumber(imap) && second->fClusterMap.TestBitNumber(imap)) {// Both clusters
1131 if (first->fSharedMap.TestBitNumber(imap) && second->fSharedMap.TestBitNumber(imap)) { // Shared
1132 sumQ++;
1133 sumCls+=2;
1134 sumSha+=2;}
1135 else {sumQ--; sumCls+=2;}
1136 }
1137 else if (first->fClusterMap.TestBitNumber(imap) || second->fClusterMap.TestBitNumber(imap)) {// Non shared
1138 sumQ++;
1139 sumCls++;}
1140 }
1141 if (sumCls>0) {
1142 qfactor = sumQ*1.0/sumCls;
1143 shfrac = sumSha*1.0/sumCls;
1144 }
1145
1146 if(qfactor > fShareQuality || shfrac > fShareFraction) return kFALSE;
1147
1148
1149 return kTRUE;
1150
1151
1152 }
1153 //________________________________________________________________________
1154 Float_t AliAnalysisTaskFemtoESE::GamovFactor(Int_t chargeBin1, Int_t chargeBin2, Float_t qinv)
1155 {
1156 Float_t arg = G_Coeff/qinv;
1157
1158 if(chargeBin1==chargeBin2) return (exp(arg)-1)/(arg);
1159 else {return (exp(-arg)-1)/(-arg);}
1160
1161 }
1162 //________________________________________________________________________
1163 void AliAnalysisTaskFemtoESE::Shuffle(Int_t *iarr, Int_t i1, Int_t i2)
1164 {
1165 Int_t j, k;
1166 Int_t a = i2 - i1;
1167 for (Int_t i = i1; i < i2+1; i++) {
1168 j = (Int_t) (gRandom->Rndm() * a);
1169 k = iarr[j];
1170 iarr[j] = iarr[i];
1171 iarr[i] = k;
1172 }
1173 }
1174 */
1175 //________________________________________________________________________
1176 Double_t AliAnalysisTaskFemtoESE::GetQinv(Double_t track1[], Double_t track2[]){
1177
1178   Double_t qinv2=1.0;
1179   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.);
1180   if(qinv2 >= 0.) return sqrt(qinv2);
1181   else return -1.*sqrt(-1.*qinv2);
1182 }
1183
1184 //________________________________________________________________________
1185 void AliAnalysisTaskFemtoESE::GetQosl(Double_t track1[], Double_t track2[], Double_t& qout, Double_t& qside, Double_t& qlong){
1186
1187   Double_t p0 = track1[0] + track2[0];
1188   Double_t px = track1[1] + track2[1];
1189   Double_t py = track1[2] + track2[2];
1190   Double_t pz = track1[3] + track2[3];
1191
1192   Double_t mt = sqrt(p0*p0 - pz*pz);
1193   Double_t pt = sqrt(px*px + py*py);
1194
1195   Double_t v0 = track1[0] - track2[0];
1196   Double_t vx = track1[1] - track2[1];
1197   Double_t vy = track1[2] - track2[2];
1198   Double_t vz = track1[3] - track2[3];
1199
1200   if(gRandom->Rndm()<0.5)
1201     {
1202       v0 = -v0;
1203       vx = -vx;
1204       vy = -vy;
1205       vz = -vz;
1206     }
1207
1208   //cout << p0 << "  " << px << "  " << py << "  " << pz << "  " << v0 << "  " << vx << "  " << vy << "  " << vz << "  " << mt << "  " << pt << endl;
1209
1210   qout = (px*vx + py*vy)/pt;
1211   qside = (px*vy - py*vx)/pt;
1212   qlong = (p0*vz - pz*v0)/mt;
1213 }
1214
1215 //________________________________________________________________________
1216 Bool_t AliAnalysisTaskFemtoESE::EventCut(/*AliAODEvent* fevent*/){
1217
1218   // Trigger Cut
1219
1220   Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
1221   Bool_t isSelected2 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral);
1222   Bool_t isSelected3 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kSemiCentral);
1223   if(!isSelected1 && !isSelected2 && !isSelected3) {return kFALSE;}
1224
1225   /*
1226   // Pile-up rejection
1227   AliAnalysisUtils *AnaUtil=new AliAnalysisUtils();
1228   if(!fPbPbcase) AnaUtil->SetUseMVPlpSelection(kTRUE);// use Multi-Vertex tool for pp and pPb
1229   else AnaUtil->SetUseMVPlpSelection(kFALSE);
1230   Bool_t pileUpCase=AnaUtil->IsPileUpEvent(fAOD);
1231   if(pileUpCase) return;
1232
1233   // Vertexing
1234   primaryVertexAOD = fAOD->GetPrimaryVertex();
1235   vertex[0]=primaryVertexAOD->GetX(); vertex[1]=primaryVertexAOD->GetY(); vertex[2]=primaryVertexAOD->GetZ();
1236
1237   if(fabs(vertex[2]) > 10) {cout<<"Zvertex Out of Range. Skip Event"<<endl; return;} // Z-Vertex Cut
1238   ((TH3D*)fOutputList->FindObject("fVertexDist"))->Fill(vertex[0], vertex[1], vertex[2]);
1239
1240   for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) {
1241   AliAODTrack* aodtrack = fAOD->GetTrack(i);
1242   if (!aodtrack) continue;
1243   AODTracks++;
1244   if(!aodtrack->TestFilterBit(BIT(fFilterBit))) continue;// AOD filterBit cut
1245   FBTracks++;
1246   }
1247   ((TH1D*)fOutputList->FindObject("fMultDist2"))->Fill(FBTracks);
1248
1249   //if(fAOD->IsPileupFromSPD()) {cout<<"PileUpEvent. Skip Event"<<endl; return;} // Old Pile-up cut
1250   if(primaryVertexAOD->GetNContributors() < 1) {cout<<"Bad Vertex. Skip Event"<<endl; return;}
1251
1252   ((TH1D*)fOutputList->FindObject("fMultDist3"))->Fill(FBTracks);
1253
1254   fBfield = fAOD->GetMagneticField();
1255
1256   for(Int_t i=0; i<fZvertexBins; i++){
1257   if( (vertex[2] >= zstart+i*zstep) && (vertex[2] < zstart+(i+1)*zstep) ){
1258   zbin=i;
1259   break;
1260   }
1261   }
1262   */
1263
1264   return kTRUE;
1265
1266 }
1267
1268 //________________________________________________________________________
1269 Bool_t AliAnalysisTaskFemtoESE::TrackCut(AliAODTrack* ftrack){
1270
1271   if (!ftrack->TestFilterBit(fFilterBit)) return kFALSE;
1272
1273   if(ftrack->Pt() < 0.14) return kFALSE;
1274   if(ftrack->Pt() > 1.5) return kFALSE;
1275   if(fabs(ftrack->Eta()) > 0.8) return kFALSE;
1276
1277   if(ftrack->GetTPCNcls() < 80) return kFALSE;// TPC nCluster cut
1278
1279   Double_t trackdca[3] = {ftrack->XAtDCA(),ftrack->YAtDCA(),ftrack->ZAtDCA()};
1280   //ftrack->XYZAtDCA(trackdca);
1281   //Double_t dcaxy = sqrt( pow(trackpos[0] - vertex[0],2) + pow(trackpos[1] - vertex[1],2));
1282   //Double_t dcaz = sqrt( pow(trackpos[2] - vertex[2],2));
1283   hdcaxy->Fill(trackdca[0],trackdca[1]);
1284   hdcaz->Fill(trackdca[2]);
1285   //if(dcaxy > 0.2) return kFALSE;
1286   //if(dcaz > 0.15) return kFALSE;
1287
1288
1289   //// FilterBit Overlap Check
1290   //if(fFilterBit != 7){
1291   //    Bool_t goodTrackOtherFB = kFALSE;
1292   //    if(fMCcase && fAOD->GetRunNumber()<=126437) goodTrackOtherFB=kTRUE;// FB7 to FB5 mapping in 10f6a MC does not work
1293   //
1294   //    for (Int_t j = 0; j < fAOD->GetNumberOfTracks(); j++) {
1295   //      AliAODTrack* aodtrack2 = fAOD->GetTrack(randomIndex[j]);
1296   //      if(!aodtrack2) continue;
1297   //      if(!aodtrack2->TestFilterBit(BIT(fFilterBit))) continue;
1298   //
1299   //      if(-(aodtrack->GetID()+1)==aodtrack2->GetID()) {goodTrackOtherFB=kTRUE; break;}
1300   //
1301   //    }
1302   //    if(!goodTrackOtherFB) continue;
1303   //    }
1304
1305   return kTRUE;
1306
1307 }
1308
1309 //________________________________________________________________________
1310 Bool_t AliAnalysisTaskFemtoESE::PairCut(AliFemtoESEBasicParticle* ftrack1, AliFemtoESEBasicParticle* ftrack2, Bool_t mix){
1311
1312   // check for same charge
1313   if(ftrack2->Charge() != ftrack1->Charge()) return kFALSE;
1314
1315   // qinv cut
1316   Double_t trackvec1[4] = {ftrack1->E(),ftrack1->Px(),ftrack1->Py(),ftrack1->Pz()};
1317   Double_t trackvec2[4] = {ftrack2->E(),ftrack2->Px(),ftrack2->Py(),ftrack2->Pz()};
1318   Double_t qinv = GetQinv(trackvec1,trackvec2);
1319   if(qinv < 0.005) return kFALSE; // qinv < 0.005
1320
1321   // deltaEta x deltaPhi* cut
1322   if(fabs(ftrack1->Eta()-ftrack2->Eta()) < fMinSepPairEta)
1323     {
1324       Double_t deltaphistar = DeltaPhiStar(ftrack1,ftrack2,1.0); // angular separation at r=1m
1325       deltaphistar = fabs(deltaphistar);
1326       if(deltaphistar < fMinSepPairPhi) return kFALSE;
1327       deltaphistar = DeltaPhiStar(ftrack1,ftrack2,1.6); // angular separation at r=1.6m
1328       deltaphistar = fabs(deltaphistar);
1329       if(deltaphistar < fMinSepPairPhi) return kFALSE;
1330     }
1331
1332   // share fraction & share quality cut
1333   TBits clusterMap1 = (TBits)(ftrack1->GetTPCClusterMap());
1334   TBits sharedMap1 = (TBits)(ftrack1->GetTPCSharedMap());
1335   TBits clusterMap2 = (TBits)(ftrack2->GetTPCClusterMap());
1336   TBits sharedMap2 = (TBits)(ftrack2->GetTPCSharedMap());
1337
1338   Int_t ncl1 = clusterMap1.GetNbits();
1339   Int_t ncl2 = clusterMap2.GetNbits();
1340   Int_t sumCls = 0; Int_t sumSha = 0; Int_t sumQ = 0;
1341   Double_t shfrac = 0; Double_t qfactor = 0;
1342   for(Int_t imap = 0; imap < ncl1 && imap < ncl2; imap++)
1343     {
1344       if (clusterMap1.TestBitNumber(imap) && clusterMap2.TestBitNumber(imap)) // Both clusters
1345         {
1346           if (sharedMap1.TestBitNumber(imap) && sharedMap2.TestBitNumber(imap)) // Shared
1347             {
1348               sumQ++;
1349               sumCls+=2;
1350               sumSha+=2;}
1351           else {sumQ--; sumCls+=2;}
1352         }
1353       else if (clusterMap1.TestBitNumber(imap) || clusterMap2.TestBitNumber(imap)) // Non shared
1354         {
1355           sumQ++;
1356           sumCls++;
1357         }
1358     }
1359   if (sumCls>0)
1360     {
1361       qfactor = sumQ*1.0/sumCls;
1362       shfrac = sumSha*1.0/sumCls;
1363     }
1364
1365   // sumCls -- number of clusters in track 1 + number of clusters in track 2 (clusters in both tracks counted twice)
1366   // sumSha -- number of shared clusters (counted twice)
1367   // sumQ -- ?
1368
1369   if(!mix)
1370     {
1371       hsharequal->Fill(qfactor);
1372       hsharefrac->Fill(shfrac);
1373     }
1374   else
1375     {
1376       hsharequalmix->Fill(qfactor);
1377       hsharefracmix->Fill(shfrac);
1378     }
1379
1380
1381   if(qfactor > fShareQuality || shfrac > fShareFraction) return kFALSE;
1382
1383   return kTRUE;
1384 }
1385
1386 Double_t AliAnalysisTaskFemtoESE::DeltaPhiStar(AliAODTrack* ftrack1, AliAODTrack* ftrack2, Double_t r)
1387 {
1388   Double_t phi1 = ftrack1->Phi() - asin(0.3*ftrack1->Charge()*(0.1*fBfield)*r/(2.*ftrack1->Pt())); // magnetic field converted from kilogauss to tesla
1389   Double_t phi2 = ftrack2->Phi() - asin(0.3*ftrack2->Charge()*(0.1*fBfield)*r/(2.*ftrack2->Pt()));
1390
1391   Double_t deltaphi = phi1 - phi2;
1392   while(deltaphi > TMath::Pi()) deltaphi -= 2*TMath::Pi();
1393   while(deltaphi < -TMath::Pi()) deltaphi += 2*TMath::Pi();
1394
1395   return deltaphi;
1396 }
1397
1398 /*TObjArray* AliAnalysisTaskFemtoESE::CloneAndReduceTrackList(TObjArray* tracks, Double_t psi)
1399   {
1400   // clones a track list by using AliDPhiBasicParticle which uses much less memory (used for event mixing)
1401
1402   TObjArray* tracksClone = new TObjArray;
1403   tracksClone->SetOwner(kTRUE);
1404
1405   for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
1406   {
1407   AliAODTrack* particle = (AliAODTrack*) tracks->UncheckedAt(i);
1408   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());
1409   copy->SetPsiEP(psi);
1410   copy->SetTPCClusterMap(particle->GetTPCClusterMap());
1411   copy->SetTPCSharedMap(particle->GetTPCSharedMap());
1412
1413   tracksClone->Add(copy);
1414   }
1415
1416   return tracksClone;
1417   }*/
1418
1419 Double_t AliAnalysisTaskFemtoESE::GetDeltaPhiEP(Double_t px1, Double_t py1, Double_t px2, Double_t py2, Double_t psi)
1420 {
1421   // angle of pair
1422   Double_t px = px1+px2;
1423   Double_t py = py1+py2;
1424   Double_t phi = atan2(py,px);
1425
1426   Double_t dphi = phi-psi;
1427   while(dphi < epBins[0]) dphi += 2*TMath::Pi();
1428   while(dphi > epBins[nEPBins]) dphi -= 2*TMath::Pi();
1429
1430   return dphi;
1431 }
1432
1433 Bool_t AliAnalysisTaskFemtoESE::FindBin(Double_t kt, Double_t phi, Double_t cent, Int_t& a, Int_t& b, Int_t&c)
1434 {
1435   a = b = c = -1;
1436   for(Int_t i = 0; i < nKtBins; i++)
1437     {
1438       if(kt >= ktBins[i] && kt < ktBins[i+1])
1439         {
1440           a = i;
1441           break;
1442         }
1443     }
1444   if(a==-1) return kFALSE;
1445
1446   for(Int_t i = 0; i < nEPBins; i++)
1447     {
1448       if(phi >= epBins[i] && phi < epBins[i+1])
1449         {
1450           b = i;
1451           break;
1452         }
1453     }
1454   if(b==-1) return kFALSE;
1455
1456   for(Int_t i = 0; i < nCentBins; i++)
1457     {
1458       if(cent >= centBins[i] && cent < centBins[i+1])
1459         {
1460           c = i;
1461           break;
1462         }
1463     }
1464   if(c==-1) return kFALSE;
1465
1466   return kTRUE;
1467 }
1468
1469
1470
1471 void AliAnalysisTaskFemtoESE::SetKtBins(Int_t n, Double_t* bins)
1472 {
1473   if(ktBins) delete [] ktBins;
1474   nKtBins = n;
1475   nKtBins1 = n+1;
1476   ktBins = new Double_t[nKtBins+1];
1477   for(Int_t i = 0; i < nKtBins+1; i++)
1478     ktBins[i]=bins[i];
1479   Printf("Setting %i kt bins: ",nKtBins);
1480   for(Int_t i = 0; i < nKtBins+1; i++) Printf("%lf",ktBins[i]);
1481 }
1482 void AliAnalysisTaskFemtoESE::SetCentBins(Int_t n, Double_t* bins)
1483 {
1484   if(centBins) delete [] centBins;
1485   nCentBins = n;
1486   nCentBins1 = n+1;
1487   centBins = new Double_t[nCentBins+1];
1488   for(Int_t i = 0; i < nCentBins+1; i++)
1489     centBins[i]=bins[i];
1490   Printf("Setting %i centrality bins: ",nCentBins);
1491   for(Int_t i = 0; i < nCentBins+1; i++) Printf("%lf",centBins[i]);
1492 }
1493 void AliAnalysisTaskFemtoESE::SetVzBins(Int_t n, Double_t* bins)
1494 {
1495   if(vzBins) delete [] vzBins;
1496   nVzBins = n;
1497   nVzBins1 = n+1;
1498   vzBins = new Double_t[nVzBins+1];
1499   for(Int_t i = 0; i < nVzBins+1; i++)
1500     vzBins[i]=bins[i];
1501   Printf("Setting %i vz bins: ",nVzBins);
1502   for(Int_t i = 0; i < nVzBins+1; i++) Printf("%lf",vzBins[i]);
1503 }
1504 void AliAnalysisTaskFemtoESE::SetEPBins(Int_t n, Double_t min, Double_t max)
1505 {
1506   if(epBins) delete [] epBins;
1507   nEPBins = n;
1508   nEPBins1 = n+1;
1509   epBins = new Double_t[nEPBins+1];
1510   for(Int_t y = 0; y < nEPBins+1; y++)
1511     epBins[y] = min+((max-min)/(Double_t)n)*((Double_t)y);
1512   Printf("Setting %i EP bins: ",nEPBins);
1513   for(Int_t i = 0; i < nEPBins+1; i++) Printf("%lf",epBins[i]);
1514 }
1515
1516 Double_t AliAnalysisTaskFemtoESE::GetQPercLHC11h(Double_t qvec)
1517 {
1518   // preliminary attemp at calculating qvector percentile in LHC11h -- still very approximate and only works in 5% bins
1519   if(!qPercBinsLHC11h)
1520     {
1521       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};
1522       nqPercBinsLHC11h = 21;
1523       qPercBinsLHC11h = new Double_t[nqPercBinsLHC11h];
1524       for(Int_t n = 0; n < nqPercBinsLHC11h; n++) qPercBinsLHC11h[n] = tempArray[n];
1525     }
1526
1527   for(Int_t t = 0; t < nqPercBinsLHC11h-1; t++)
1528     if(qvec > qPercBinsLHC11h[t] && qvec < qPercBinsLHC11h[t+1]) return 50.*(2*t+1)/(Double_t)(nqPercBinsLHC11h-1);
1529
1530   if(qvec < qPercBinsLHC11h[0]) return 0.0;
1531   if(qvec > qPercBinsLHC11h[nqPercBinsLHC11h-1]) return 100.0;
1532
1533   return 0.0;
1534
1535 }
1536
1537 Double_t AliAnalysisTaskFemtoESE::GetCentralityWeight(Double_t cent)
1538 {
1539
1540   // use makeCentWeighting.C to fit centrality distribution to obtain this parameterization
1541   Double_t par1 = 2.60629;
1542   Double_t par2 = 0.579333;
1543   Double_t limit1 = 8.71488;
1544   Double_t limit2 = 12.0126;
1545
1546   if(cent < limit1) return 1./par1;
1547   if(cent > limit2) return 1./par2;
1548
1549   Double_t slope = (par2-par1)/(limit2-limit1);
1550   Double_t b = par1-slope*limit1;
1551
1552   return 1./(slope*cent+b);
1553 }