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