PWGHF/hfe/AliAnalysisTaskHFEFlow.cxx

   1 /**************************************************************************\r
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
   3 *                                                                        *\r
   4 * Author: The ALICE Off-line Project.                                    *\r
   5 * Contributors are mentioned in the code where appropriate.              *\r
   6 *                                                                        *\r
   7 * Permission to use, copy, modify and distribute this software and its   *\r
   8 * documentation strictly for non-commercial purposes is hereby granted   *\r
   9 * without fee, provided that the above copyright notice appears in all   *\r
  10 * copies and that both the copyright notice and this permission notice   *\r
  11 * appear in the supporting documentation. The authors make no claims     *\r
  12 * about the suitability of this software for any purpose. It is          *\r
  13 * provided "as is" without express or implied warranty.                  *\r
  14 **************************************************************************/\r
  15 //\r
  16 // Flow task\r
  17 // \r
  18 // Authors:\r
  19 //   Raphaelle Bailhache <R.Bailhache@gsi.de>\r
  20 //\r
  21 #include "TROOT.h"\r
  22 #include "TH1D.h"\r
  23 #include "TH2D.h"\r
  24 #include "TChain.h"\r
  25 #include "TVector2.h"\r
  26 #include "THnSparse.h"\r
  27 #include "TMath.h"\r
  28 #include "TRandom3.h"\r
  29 #include "TProfile.h"\r
  30 #include "TProfile2D.h"\r
  31 \r
  32 #include "AliAnalysisTaskSE.h"\r
  33 \r
  34 #include "AliESDInputHandler.h"\r
  35 #include "AliMCEvent.h"\r
  36 #include "AliESD.h"\r
  37 #include "AliESDEvent.h"\r
  38 #include "AliPIDResponse.h"\r
  39 #include "AliESDVZERO.h"\r
  40 #include "AliESDUtils.h"\r
  41 #include "AliMCParticle.h"\r
  42 \r
  43 #include "AliFlowCandidateTrack.h"\r
  44 #include "AliFlowEvent.h"\r
  45 #include "AliFlowTrackCuts.h"\r
  46 #include "AliFlowVector.h"\r
  47 #include "AliFlowCommonConstants.h"\r
  48 \r
  49 #include "AliHFEcuts.h"\r
  50 #include "AliHFEpid.h"\r
  51 #include "AliHFEpidQAmanager.h"\r
  52 #include "AliHFEtools.h"\r
  53 #include "AliHFEVZEROEventPlane.h"\r
  54 \r
  55 #include "AliCentrality.h"\r
  56 #include "AliEventplane.h"\r
  57 #include "AliAnalysisTaskHFEFlow.h"\r
  58 \r
  59 \r
  60 //____________________________________________________________________\r
  61 AliAnalysisTaskHFEFlow::AliAnalysisTaskHFEFlow() :\r
  62   AliAnalysisTaskSE(),\r
  63   fListHist(0x0), \r
  64   fVZEROEventPlane(kFALSE),\r
  65   fVZEROEventPlaneA(kFALSE),\r
  66   fVZEROEventPlaneC(kFALSE),\r
  67   fSubEtaGapTPC(kFALSE),\r
  68   fEtaGap(0.0),\r
  69   fNbBinsCentralityQCumulant(4),\r
  70   fNbBinsPtQCumulant(15),\r
  71   fMinPtQCumulant(0.0),\r
  72   fMaxPtQCumulant(6.0),\r
  73   fAfterBurnerOn(kFALSE),\r
  74   fNonFlowNumberOfTrackClones(0),\r
  75   fV1(0.),\r
  76   fV2(0.),\r
  77   fV3(0.),\r
  78   fV4(0.),\r
  79   fV5(0.),\r
  80   fMaxNumberOfIterations(100),\r
  81   fPrecisionPhi(0.001),\r
  82   fUseMCReactionPlane(kFALSE),\r
  83   fMCPID(kFALSE),\r
  84   fNoPID(kFALSE),\r
  85   fDebugLevel(0),\r
  86   fcutsRP(0),\r
  87   fcutsPOI(0),\r
  88   fHFECuts(0),\r
  89   fPID(0),\r
  90   fPIDqa(0),\r
  91   fflowEvent(0x0),\r
  92   fHFEVZEROEventPlane(0x0),\r
  93   fHistEV(0),\r
  94   fEventPlane(0x0),\r
  95   fEventPlaneaftersubtraction(0x0),\r
  96   fCosSin2phiep(0x0),\r
  97   fCos2phie(0x0),\r
  98   fSin2phie(0x0),\r
  99   fCos2phiep(0x0),\r
 100   fSin2phiep(0x0),\r
 101   fSin2phiephiep(0x0),\r
 102   fCosResabc(0x0),\r
 103   fSinResabc(0x0),\r
 104   fProfileCosResab(0x0),\r
 105   fProfileCosResac(0x0),\r
 106   fProfileCosResbc(0x0),\r
 107   fCosRes(0x0),\r
 108   fSinRes(0x0),\r
 109   fProfileCosRes(0x0),\r
 110   fDeltaPhiMaps(0x0),\r
 111   fCosPhiMaps(0x0),\r
 112   fProfileCosPhiMaps(0x0)\r
 113 {\r
 114   // Constructor\r
 115 \r
 116   for(Int_t k = 0; k < 10; k++) {\r
 117     fBinCentralityLess[k] = 0.0;\r
 118   }\r
 119   \r
 120 }\r
 121 //______________________________________________________________________________\r
 122 AliAnalysisTaskHFEFlow:: AliAnalysisTaskHFEFlow(const char *name) :\r
 123   AliAnalysisTaskSE(name),\r
 124   fListHist(0x0), \r
 125   fVZEROEventPlane(kFALSE),\r
 126   fVZEROEventPlaneA(kFALSE),\r
 127   fVZEROEventPlaneC(kFALSE),\r
 128   fSubEtaGapTPC(kFALSE),\r
 129   fEtaGap(0.0),\r
 130   fNbBinsCentralityQCumulant(4),\r
 131   fNbBinsPtQCumulant(15),\r
 132   fMinPtQCumulant(0.0),\r
 133   fMaxPtQCumulant(6.0),\r
 134   fAfterBurnerOn(kFALSE),\r
 135   fNonFlowNumberOfTrackClones(0),\r
 136   fV1(0.),\r
 137   fV2(0.),\r
 138   fV3(0.),\r
 139   fV4(0.),\r
 140   fV5(0.),\r
 141   fMaxNumberOfIterations(100),\r
 142   fPrecisionPhi(0.001),\r
 143   fUseMCReactionPlane(kFALSE),\r
 144   fMCPID(kFALSE),\r
 145   fNoPID(kFALSE),\r
 146   fDebugLevel(0),\r
 147   fcutsRP(0),\r
 148   fcutsPOI(0),\r
 149   fHFECuts(0),\r
 150   fPID(0),\r
 151   fPIDqa(0),\r
 152   fflowEvent(0x0),\r
 153   fHFEVZEROEventPlane(0x0),\r
 154   fHistEV(0),\r
 155   fEventPlane(0x0),\r
 156   fEventPlaneaftersubtraction(0x0),\r
 157   fCosSin2phiep(0x0),\r
 158   fCos2phie(0x0),\r
 159   fSin2phie(0x0),\r
 160   fCos2phiep(0x0),\r
 161   fSin2phiep(0x0),\r
 162   fSin2phiephiep(0x0),\r
 163   fCosResabc(0x0),\r
 164   fSinResabc(0x0),\r
 165   fProfileCosResab(0x0),\r
 166   fProfileCosResac(0x0),\r
 167   fProfileCosResbc(0x0),\r
 168   fCosRes(0x0),\r
 169   fSinRes(0x0),\r
 170   fProfileCosRes(0x0),\r
 171   fDeltaPhiMaps(0x0),\r
 172   fCosPhiMaps(0x0),\r
 173   fProfileCosPhiMaps(0x0)\r
 174 {\r
 175   //\r
 176   // named ctor\r
 177   //\r
 178   \r
 179   for(Int_t k = 0; k < 10; k++) {\r
 180     fBinCentralityLess[k] = 0.0;\r
 181   }\r
 182   fBinCentralityLess[0] = 0.0;\r
 183   fBinCentralityLess[1] = 20.0;\r
 184   fBinCentralityLess[2] = 40.0;\r
 185   fBinCentralityLess[3] = 60.0;\r
 186   fBinCentralityLess[4] = 80.0;\r
 187   \r
 188   fPID = new AliHFEpid("hfePid");\r
 189   fPIDqa = new AliHFEpidQAmanager;\r
 190 \r
 191   DefineInput(0,TChain::Class());\r
 192   DefineOutput(1, TList::Class());\r
 193   for(Int_t bincless = 0; bincless < fNbBinsCentralityQCumulant; bincless++) {\r
 194     DefineOutput(bincless+2,AliFlowEventSimple::Class()); \r
 195   }\r
 196   \r
 197 }\r
 198 //____________________________________________________________\r
 199 AliAnalysisTaskHFEFlow::AliAnalysisTaskHFEFlow(const AliAnalysisTaskHFEFlow &ref):\r
 200   AliAnalysisTaskSE(ref),\r
 201   fListHist(0x0), \r
 202   fVZEROEventPlane(ref.fVZEROEventPlane),\r
 203   fVZEROEventPlaneA(ref.fVZEROEventPlaneA),\r
 204   fVZEROEventPlaneC(ref.fVZEROEventPlaneC),\r
 205   fSubEtaGapTPC(ref.fSubEtaGapTPC),\r
 206   fEtaGap(ref.fEtaGap),\r
 207   fNbBinsCentralityQCumulant(ref.fNbBinsCentralityQCumulant),\r
 208   fNbBinsPtQCumulant(ref.fNbBinsPtQCumulant),\r
 209   fMinPtQCumulant(ref.fMinPtQCumulant),\r
 210   fMaxPtQCumulant(ref.fMaxPtQCumulant),\r
 211   fAfterBurnerOn(ref.fAfterBurnerOn),\r
 212   fNonFlowNumberOfTrackClones(ref.fNonFlowNumberOfTrackClones),\r
 213   fV1(ref.fV1),\r
 214   fV2(ref.fV2),\r
 215   fV3(ref.fV3),\r
 216   fV4(ref.fV4),\r
 217   fV5(ref.fV5),\r
 218   fMaxNumberOfIterations(ref.fMaxNumberOfIterations),\r
 219   fPrecisionPhi(ref.fPrecisionPhi),\r
 220   fUseMCReactionPlane(ref.fUseMCReactionPlane),\r
 221   fMCPID(ref.fMCPID),\r
 222   fNoPID(ref.fNoPID),\r
 223   fDebugLevel(ref.fDebugLevel),\r
 224   fcutsRP(0),\r
 225   fcutsPOI(0),\r
 226   fHFECuts(0),\r
 227   fPID(0),\r
 228   fPIDqa(0),\r
 229   fflowEvent(0x0),\r
 230   fHFEVZEROEventPlane(0x0),\r
 231   fHistEV(0),\r
 232   fEventPlane(0x0),\r
 233   fEventPlaneaftersubtraction(0x0),\r
 234   fCosSin2phiep(0x0),\r
 235   fCos2phie(0x0),\r
 236   fSin2phie(0x0),\r
 237   fCos2phiep(0x0),\r
 238   fSin2phiep(0x0),\r
 239   fSin2phiephiep(0x0),\r
 240   fCosResabc(0x0),\r
 241   fSinResabc(0x0),\r
 242   fProfileCosResab(0x0),\r
 243   fProfileCosResac(0x0),\r
 244   fProfileCosResbc(0x0),\r
 245   fCosRes(0x0),\r
 246   fSinRes(0x0),\r
 247   fProfileCosRes(0x0),\r
 248   fDeltaPhiMaps(0x0),\r
 249   fCosPhiMaps(0x0),\r
 250   fProfileCosPhiMaps(0x0)\r
 251 {\r
 252   //\r
 253   // Copy Constructor\r
 254   //\r
 255   ref.Copy(*this);\r
 256 }\r
 257 \r
 258 //____________________________________________________________\r
 259 AliAnalysisTaskHFEFlow &AliAnalysisTaskHFEFlow::operator=(const AliAnalysisTaskHFEFlow &ref){\r
 260   //\r
 261   // Assignment operator\r
 262   //\r
 263   if(this == &ref) \r
 264     ref.Copy(*this);\r
 265   return *this;\r
 266 }\r
 267 \r
 268 //____________________________________________________________\r
 269 void AliAnalysisTaskHFEFlow::Copy(TObject &o) const {\r
 270   // \r
 271   // Copy into object o\r
 272   //\r
 273   AliAnalysisTaskHFEFlow &target = dynamic_cast<AliAnalysisTaskHFEFlow &>(o);\r
 274   target.fVZEROEventPlane = fVZEROEventPlane;\r
 275   target.fVZEROEventPlaneA = fVZEROEventPlaneA;\r
 276   target.fVZEROEventPlaneC = fVZEROEventPlaneC;\r
 277   target.fSubEtaGapTPC = fSubEtaGapTPC;\r
 278   target.fEtaGap = fEtaGap;\r
 279   target.fNbBinsCentralityQCumulant = fNbBinsCentralityQCumulant;\r
 280   target.fNbBinsPtQCumulant = fNbBinsPtQCumulant;\r
 281   target.fMinPtQCumulant = fMinPtQCumulant;\r
 282   target.fMaxPtQCumulant = fMaxPtQCumulant;\r
 283   target.fAfterBurnerOn = fAfterBurnerOn;\r
 284   target.fNonFlowNumberOfTrackClones = fNonFlowNumberOfTrackClones;\r
 285   target.fV1 = fV1;\r
 286   target.fV2 = fV2;\r
 287   target.fV3 = fV3;\r
 288   target.fV4 = fV4;\r
 289   target.fV5 = fV5;\r
 290   target.fMaxNumberOfIterations = fMaxNumberOfIterations;\r
 291   target.fPrecisionPhi = fPrecisionPhi;\r
 292   target.fUseMCReactionPlane = fUseMCReactionPlane;\r
 293   target.fMCPID = fMCPID;\r
 294   target.fNoPID = fNoPID;\r
 295   target.fDebugLevel = fDebugLevel;\r
 296   target.fcutsRP = fcutsRP;\r
 297   target.fcutsPOI = fcutsPOI;\r
 298   target.fHFECuts = fHFECuts;\r
 299   target.fPID = fPID;\r
 300   target.fPIDqa = fPIDqa;\r
 301   target.fHFEVZEROEventPlane = fHFEVZEROEventPlane;\r
 302  \r
 303 }\r
 304 //____________________________________________________________\r
 305 AliAnalysisTaskHFEFlow::~AliAnalysisTaskHFEFlow(){\r
 306   //\r
 307   // Destructor\r
 308   //\r
 309   if(fListHist) delete fListHist;\r
 310   if(fcutsRP) delete fcutsRP;\r
 311   if(fcutsPOI) delete fcutsPOI;\r
 312   if(fHFECuts) delete fHFECuts;\r
 313   if(fPID) delete fPID;\r
 314   if(fPIDqa) delete fPIDqa;\r
 315   if(fflowEvent) delete fflowEvent;\r
 316   if(fHFEVZEROEventPlane) delete fHFEVZEROEventPlane;\r
 317   if(fHistEV) delete fHistEV;\r
 318   if(fEventPlane) delete fEventPlane;\r
 319   if(fEventPlaneaftersubtraction) delete fEventPlaneaftersubtraction;\r
 320   if(fCosSin2phiep) delete fCosSin2phiep;\r
 321   if(fCos2phie) delete fCos2phie;\r
 322   if(fSin2phie) delete fSin2phie;\r
 323   if(fCos2phiep) delete fCos2phiep;\r
 324   if(fSin2phiep) delete fSin2phiep;\r
 325   if(fSin2phiephiep) delete fSin2phiephiep;\r
 326   if(fCosResabc) delete fCosResabc;\r
 327   if(fSinResabc) delete fSinResabc;\r
 328   if(fProfileCosResab) delete fProfileCosResab;\r
 329   if(fProfileCosResac) delete fProfileCosResac;\r
 330   if(fProfileCosResbc) delete fProfileCosResbc;\r
 331   if(fCosRes) delete fCosRes;\r
 332   if(fSinRes) delete fSinRes;\r
 333   if(fProfileCosRes) delete fProfileCosRes;\r
 334   if(fDeltaPhiMaps) delete fDeltaPhiMaps;\r
 335   if(fCosPhiMaps) delete fCosPhiMaps;\r
 336   if(fProfileCosPhiMaps) delete fProfileCosPhiMaps;\r
 337 \r
 338 }\r
 339 //________________________________________________________________________\r
 340 void AliAnalysisTaskHFEFlow::UserCreateOutputObjects()\r
 341 {\r
 342 \r
 343   //********************\r
 344   // Create histograms\r
 345   //********************\r
 346 \r
 347   //**************\r
 348   // Cuts\r
 349   //**************\r
 350 \r
 351   //---------Data selection----------\r
 352   //kMC, kGlobal, kTPCstandalone, kSPDtracklet, kPMD\r
 353   //AliFlowTrackCuts::trackParameterType rptype = AliFlowTrackCuts::kGlobal;\r
 354   //AliFlowTrackCuts::trackParameterType poitype = AliFlowTrackCuts::kGlobal;\r
 355 \r
 356   //---------Parameter mixing--------\r
 357   //kPure - no mixing, kTrackWithMCkine, kTrackWithMCPID, kTrackWithMCpt\r
 358   //AliFlowTrackCuts::trackParameterMix rpmix = AliFlowTrackCuts::kPure;\r
 359   //AliFlowTrackCuts::trackParameterMix poimix = AliFlowTrackCuts::kPure;\r
 360 \r
 361   // RP TRACK CUTS:\r
 362   fcutsRP = AliFlowTrackCuts::GetStandardTPCStandaloneTrackCuts2010();\r
 363   fcutsRP->SetName("StandartTPC");\r
 364   fcutsRP->SetEtaRange(-0.9,0.9);\r
 365   fcutsRP->SetQA(kTRUE);\r
 366   TList *qaCutsRP = fcutsRP->GetQA();\r
 367   qaCutsRP->SetName("QA_StandartTPC_RP");\r
 368 \r
 369   //POI TRACK CUTS:\r
 370   fcutsPOI = new AliFlowTrackCuts("dummy");\r
 371   fcutsPOI->SetParamType(AliFlowTrackCuts::kGlobal);\r
 372   fcutsPOI->SetPtRange(+1,-1); // select nothing QUICK\r
 373   fcutsPOI->SetEtaRange(+1,-1); // select nothing VZERO\r
 374 \r
 375   // Flow\r
 376   AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();\r
 377   cc->SetNbinsMult(10000);\r
 378   cc->SetMultMin(0);\r
 379   cc->SetMultMax(10000.);\r
 380   cc->SetNbinsPt(fNbBinsPtQCumulant);\r
 381   cc->SetPtMin(fMinPtQCumulant);\r
 382   cc->SetPtMax(fMaxPtQCumulant);\r
 383   cc->SetNbinsPhi(180);\r
 384   cc->SetPhiMin(0.0);\r
 385   cc->SetPhiMax(TMath::TwoPi());\r
 386   cc->SetNbinsEta(200);\r
 387   cc->SetEtaMin(-0.9);\r
 388   cc->SetEtaMax(+0.9);\r
 389   cc->SetNbinsQ(500);\r
 390   cc->SetQMin(0.0);\r
 391   cc->SetQMax(3.0);\r
 392 \r
 393   \r
 394   // HFE cuts\r
 395 \r
 396   if(!fHFECuts){\r
 397     fHFECuts = new AliHFEcuts;\r
 398     fHFECuts->CreateStandardCuts();\r
 399   }\r
 400   fHFECuts->Initialize();\r
 401   \r
 402   // PID HFE\r
 403   //fPID->SetHasMCData(HasMCData());\r
 404   if(!fPID->GetNumberOfPIDdetectors()) fPID->AddDetector("TPC", 0);\r
 405   fPID->InitializePID();\r
 406   fPIDqa->Initialize(fPID);\r
 407   fPID->SortDetectors();\r
 408   \r
 409   //**************************\r
 410   // Bins for the THnSparse\r
 411   //**************************\r
 412 \r
 413   Int_t nBinsPt = 25;\r
 414   Double_t minPt = 0.001;\r
 415   Double_t maxPt = 10.0;\r
 416   Double_t binLimLogPt[nBinsPt+1];\r
 417   Double_t binLimPt[nBinsPt+1];\r
 418   for(Int_t i=0; i<=nBinsPt; i++) binLimLogPt[i]=(Double_t)TMath::Log10(minPt) + (TMath::Log10(maxPt)-TMath::Log10(minPt))/nBinsPt*(Double_t)i ;\r
 419   for(Int_t i=0; i<=nBinsPt; i++) binLimPt[i]=(Double_t)TMath::Power(10,binLimLogPt[i]);\r
 420 \r
 421   Int_t nBinsPtPlus = 15;\r
 422   Double_t minPtPlus = 0.0;\r
 423   Double_t maxPtPlus = 6.0;\r
 424   Double_t binLimPtPlus[nBinsPtPlus+1];\r
 425   for(Int_t i=0; i<=nBinsPtPlus; i++) binLimPtPlus[i]=(Double_t)minPtPlus + (maxPtPlus-minPtPlus)/nBinsPtPlus*(Double_t)i ;\r
 426 \r
 427   Int_t nBinsEta = 8;\r
 428   Double_t minEta = -0.8;\r
 429   Double_t maxEta = 0.8;\r
 430   Double_t binLimEta[nBinsEta+1];\r
 431   for(Int_t i=0; i<=nBinsEta; i++) binLimEta[i]=(Double_t)minEta + (maxEta-minEta)/nBinsEta*(Double_t)i ;\r
 432  \r
 433   Int_t nBinsCos = 50;\r
 434   Double_t minCos = -1.0;\r
 435   Double_t maxCos = 1.0;\r
 436   Double_t binLimCos[nBinsCos+1];\r
 437   for(Int_t i=0; i<=nBinsCos; i++) binLimCos[i]=(Double_t)minCos + (maxCos-minCos)/nBinsCos*(Double_t)i ;\r
 438  \r
 439   Int_t nBinsC = 11;\r
 440   Double_t minC = 0.0;\r
 441   Double_t maxC = 11.0;\r
 442   Double_t binLimC[nBinsC+1];\r
 443   for(Int_t i=0; i<=nBinsC; i++) binLimC[i]=(Double_t)minC + (maxC-minC)/nBinsC*(Double_t)i ;\r
 444 \r
 445   Int_t nBinsCMore = 20;\r
 446   Double_t minCMore = 0.0;\r
 447   Double_t maxCMore = 20.0;\r
 448   Double_t binLimCMore[nBinsCMore+1];\r
 449   for(Int_t i=0; i<=nBinsCMore; i++) binLimCMore[i]=(Double_t)minCMore + (maxCMore-minCMore)/nBinsCMore*(Double_t)i ;\r
 450 \r
 451   Int_t nBinsPhi = 25;\r
 452   Double_t minPhi = 0.0;\r
 453   Double_t maxPhi = TMath::Pi();\r
 454   Double_t binLimPhi[nBinsPhi+1];\r
 455   for(Int_t i=0; i<=nBinsPhi; i++) {\r
 456     binLimPhi[i]=(Double_t)minPhi + (maxPhi-minPhi)/nBinsPhi*(Double_t)i ;\r
 457     //printf("bin phi is %f for %d\n",binLimPhi[i],i);\r
 458   }\r
 459   \r
 460   //******************\r
 461   // Histograms\r
 462   //******************\r
 463     \r
 464   fListHist = new TList();\r
 465 \r
 466   // Histos\r
 467   fHistEV = new TH2D("fHistEV", "events", 3, 0, 3, 3, 0,3);\r
 468   \r
 469   // Event plane as function of phiep, centrality\r
 470   const Int_t nDima=5;\r
 471   Int_t nBina[nDima] = {nBinsPhi,nBinsPhi,nBinsPhi,nBinsPhi,nBinsC};\r
 472   fEventPlane = new THnSparseF("EventPlane","EventPlane",nDima,nBina);\r
 473   fEventPlane->SetBinEdges(0,binLimPhi);\r
 474   fEventPlane->SetBinEdges(1,binLimPhi);\r
 475   fEventPlane->SetBinEdges(2,binLimPhi);\r
 476   fEventPlane->SetBinEdges(3,binLimPhi);\r
 477   fEventPlane->SetBinEdges(4,binLimC);\r
 478   fEventPlane->Sumw2();\r
 479   \r
 480   // Event Plane after subtraction as function of phiep, centrality, pt, eta\r
 481   const Int_t nDimb=2;\r
 482   Int_t nBinb[nDimb] = {nBinsPhi, nBinsC};\r
 483   fEventPlaneaftersubtraction = new THnSparseF("EventPlane_aftersubtraction","EventPlane_aftersubtraction",nDimb,nBinb);\r
 484   fEventPlaneaftersubtraction->SetBinEdges(0,binLimPhi);\r
 485   fEventPlaneaftersubtraction->SetBinEdges(1,binLimC);\r
 486   fEventPlaneaftersubtraction->Sumw2();\r
 487 \r
 488   // Monitoring of the event Plane cos(2phi) sin(2phi) centrality\r
 489   const Int_t nDimi=3;\r
 490   Int_t nBini[nDimi] = {nBinsCos, nBinsCos, nBinsCMore};\r
 491   fCosSin2phiep = new THnSparseF("CosSin2phiep","CosSin2phiep",nDimi,nBini);\r
 492   fCosSin2phiep->SetBinEdges(0,binLimCos);\r
 493   fCosSin2phiep->SetBinEdges(1,binLimCos);\r
 494   fCosSin2phiep->SetBinEdges(2,binLimCMore);\r
 495   fCosSin2phiep->Sumw2();\r
 496 \r
 497   // Monitoring Event plane after subtraction of the track\r
 498   const Int_t nDime=4;\r
 499   Int_t nBine[nDime] = {nBinsCos, nBinsC, nBinsPt, nBinsEta};\r
 500   fCos2phie = new THnSparseF("cos2phie","cos2phie",nDime,nBine);\r
 501   fCos2phie->SetBinEdges(2,binLimPt);\r
 502   fCos2phie->SetBinEdges(3,binLimEta);\r
 503   fCos2phie->SetBinEdges(0,binLimCos);\r
 504   fCos2phie->SetBinEdges(1,binLimC);\r
 505   fCos2phie->Sumw2();\r
 506   fSin2phie = new THnSparseF("sin2phie","sin2phie",nDime,nBine);\r
 507   fSin2phie->SetBinEdges(2,binLimPt);\r
 508   fSin2phie->SetBinEdges(3,binLimEta);\r
 509   fSin2phie->SetBinEdges(0,binLimCos);\r
 510   fSin2phie->SetBinEdges(1,binLimC);\r
 511   fSin2phie->Sumw2();\r
 512   fCos2phiep = new THnSparseF("cos2phiep","cos2phiep",nDime,nBine);\r
 513   fCos2phiep->SetBinEdges(2,binLimPt);\r
 514   fCos2phiep->SetBinEdges(3,binLimEta);\r
 515   fCos2phiep->SetBinEdges(0,binLimCos);\r
 516   fCos2phiep->SetBinEdges(1,binLimC);\r
 517   fCos2phiep->Sumw2();\r
 518   fSin2phiep = new THnSparseF("sin2phiep","sin2phiep",nDime,nBine);\r
 519   fSin2phiep->SetBinEdges(2,binLimPt);\r
 520   fSin2phiep->SetBinEdges(3,binLimEta);\r
 521   fSin2phiep->SetBinEdges(0,binLimCos);\r
 522   fSin2phiep->SetBinEdges(1,binLimC);\r
 523   fSin2phiep->Sumw2();\r
 524   fSin2phiephiep = new THnSparseF("sin2phie_phiep","sin2phie_phiep",nDime,nBine);\r
 525   fSin2phiephiep->SetBinEdges(2,binLimPt);\r
 526   fSin2phiephiep->SetBinEdges(3,binLimEta);\r
 527   fSin2phiephiep->SetBinEdges(0,binLimCos);\r
 528   fSin2phiephiep->SetBinEdges(1,binLimC);\r
 529   fSin2phiephiep->Sumw2();  \r
 530 \r
 531   // Resolution cosres_abc centrality\r
 532   const Int_t nDimfbis=4;\r
 533   Int_t nBinfbis[nDimfbis] = {nBinsCos,nBinsCos,nBinsCos,nBinsCMore};\r
 534   fCosResabc = new THnSparseF("CosRes_abc","CosRes_abc",nDimfbis,nBinfbis);\r
 535   fCosResabc->SetBinEdges(0,binLimCos);\r
 536   fCosResabc->SetBinEdges(1,binLimCos);\r
 537   fCosResabc->SetBinEdges(2,binLimCos);\r
 538   fCosResabc->SetBinEdges(3,binLimCMore);\r
 539   fCosResabc->Sumw2();\r
 540 \r
 541   const Int_t nDimfbiss=4;\r
 542   Int_t nBinfbiss[nDimfbiss] = {nBinsCos,nBinsCos,nBinsCos,nBinsC};\r
 543   fSinResabc = new THnSparseF("SinRes_abc","SinRes_abc",nDimfbiss,nBinfbiss);\r
 544   fSinResabc->SetBinEdges(0,binLimCos);\r
 545   fSinResabc->SetBinEdges(1,binLimCos);\r
 546   fSinResabc->SetBinEdges(2,binLimCos);\r
 547   fSinResabc->SetBinEdges(3,binLimC);\r
 548   fSinResabc->Sumw2();\r
 549 \r
 550   // Profile cosres centrality with 3 subevents\r
 551   fProfileCosResab = new TProfile("ProfileCosRes_a_b","ProfileCosRes_a_b",nBinsCMore,binLimCMore);\r
 552   fProfileCosResab->Sumw2();\r
 553   fProfileCosResac = new TProfile("ProfileCosRes_a_c","ProfileCosRes_a_c",nBinsCMore,binLimCMore);\r
 554   fProfileCosResac->Sumw2();\r
 555   fProfileCosResbc = new TProfile("ProfileCosRes_b_c","ProfileCosRes_b_c",nBinsCMore,binLimCMore);\r
 556   fProfileCosResbc->Sumw2();\r
 557   \r
 558   // Resolution cosres centrality\r
 559   const Int_t nDimf=2;\r
 560   Int_t nBinf[nDimf] = {nBinsCos, nBinsCMore};\r
 561   fCosRes = new THnSparseF("CosRes","CosRes",nDimf,nBinf);\r
 562   fCosRes->SetBinEdges(0,binLimCos);\r
 563   fCosRes->SetBinEdges(1,binLimCMore);\r
 564   fCosRes->Sumw2();\r
 565 \r
 566   const Int_t nDimff=2;\r
 567   Int_t nBinff[nDimff] = {nBinsCos, nBinsC};\r
 568   fSinRes = new THnSparseF("SinRes","SinRes",nDimff,nBinff);\r
 569   fSinRes->SetBinEdges(0,binLimCos);\r
 570   fSinRes->SetBinEdges(1,binLimC);\r
 571   fSinRes->Sumw2();\r
 572 \r
 573   // Profile cosres centrality\r
 574   fProfileCosRes = new TProfile("ProfileCosRes","ProfileCosRes",nBinsCMore,binLimCMore);\r
 575   fProfileCosRes->Sumw2();\r
 576   \r
 577   // Maps delta phi\r
 578   const Int_t nDimg=3;\r
 579   Int_t nBing[nDimg] = {nBinsPhi,nBinsC,nBinsPt};\r
 580   fDeltaPhiMaps = new THnSparseF("DeltaPhiMaps","DeltaPhiMaps",nDimg,nBing);\r
 581   fDeltaPhiMaps->SetBinEdges(0,binLimPhi);\r
 582   fDeltaPhiMaps->SetBinEdges(1,binLimC);\r
 583   fDeltaPhiMaps->SetBinEdges(2,binLimPt);\r
 584   fDeltaPhiMaps->Sumw2();  \r
 585 \r
 586   // Maps cos phi\r
 587   const Int_t nDimh=3;\r
 588   Int_t nBinh[nDimh] = {nBinsCos,nBinsC,nBinsPt};\r
 589   fCosPhiMaps = new THnSparseF("CosPhiMaps","CosPhiMaps",nDimh,nBinh);\r
 590   fCosPhiMaps->SetBinEdges(0,binLimCos);\r
 591   fCosPhiMaps->SetBinEdges(1,binLimC);\r
 592   fCosPhiMaps->SetBinEdges(2,binLimPt);\r
 593   fCosPhiMaps->Sumw2();\r
 594 \r
 595   // Profile Maps cos phi\r
 596   fProfileCosPhiMaps = new TProfile2D("ProfileCosPhiMaps","ProfileCosPhiMaps",fNbBinsCentralityQCumulant,&fBinCentralityLess[0],nBinsPtPlus,binLimPtPlus);\r
 597   fProfileCosPhiMaps->Sumw2();\r
 598 \r
 599 \r
 600   //**************************\r
 601   // Add to the list\r
 602   //******************************\r
 603 \r
 604   fListHist->Add(qaCutsRP);\r
 605   fListHist->Add(fPIDqa->MakeList("HFEpidQA"));\r
 606   fListHist->Add(fHistEV);\r
 607   fListHist->Add(fProfileCosRes);\r
 608   fListHist->Add(fProfileCosResab);\r
 609   fListHist->Add(fProfileCosResac);\r
 610   fListHist->Add(fProfileCosResbc);\r
 611   fListHist->Add(fCosSin2phiep);\r
 612   fListHist->Add(fEventPlane);\r
 613   fListHist->Add(fEventPlaneaftersubtraction);\r
 614   fListHist->Add(fCos2phie);\r
 615   fListHist->Add(fSin2phie);\r
 616   fListHist->Add(fCos2phiep);\r
 617   fListHist->Add(fSin2phiep);\r
 618   fListHist->Add(fSin2phiephiep);\r
 619   fListHist->Add(fCosRes);\r
 620   fListHist->Add(fCosResabc);\r
 621   fListHist->Add(fSinRes);\r
 622   fListHist->Add(fSinResabc);\r
 623   fListHist->Add(fDeltaPhiMaps);\r
 624   fListHist->Add(fCosPhiMaps);\r
 625   fListHist->Add(fProfileCosPhiMaps);\r
 626 \r
 627   if(fHFEVZEROEventPlane && (fDebugLevel > 2)) fListHist->Add(fHFEVZEROEventPlane->GetOutputList());\r
 628   \r
 629 \r
 630   PostData(1, fListHist);\r
 631 \r
 632 \r
 633 }\r
 634    \r
 635 //________________________________________________________________________\r
 636 void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/)\r
 637 {\r
 638   //\r
 639   // Loop over event\r
 640   //\r
 641    \r
 642   Double_t massElectron = 0.000511;\r
 643   Double_t mcReactionPlane = 0.0;\r
 644   Bool_t   eventplanedefined = kTRUE;\r
 645 \r
 646   Float_t cntr = 0.0;\r
 647   Double_t binct = 11.5;\r
 648   Double_t binctMore = 20.5;\r
 649   Double_t binctLess = -0.5;\r
 650   Float_t binctt = -1.0;\r
 651   \r
 652   Double_t valuecossinephiep[3];\r
 653   Double_t valuensparsea[5];\r
 654   Double_t valuensparseabis[5];\r
 655   Double_t valuensparsee[4];\r
 656   Double_t valuensparsef[2];\r
 657   Double_t valuensparsefsin[2];\r
 658   Double_t valuensparsefbis[4];\r
 659   Double_t valuensparsefbissin[4];\r
 660   Double_t valuensparseg[3];\r
 661   Double_t valuensparseh[3];\r
 662   Double_t valuensparsehprofile[3];\r
 663 \r
 664   AliMCEvent *mcEvent = MCEvent();\r
 665   AliMCParticle *mctrack = NULL;\r
 666   \r
 667   /////////////////\r
 668   // centrality\r
 669   /////////////////\r
 670   \r
 671   AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());\r
 672   if(!esd) return;\r
 673   AliCentrality *centrality = esd->GetCentrality();\r
 674   if(!centrality) return;\r
 675   cntr = centrality->GetCentralityPercentile("V0M");\r
 676   if((0.0< cntr) && (cntr<5.0)) binct = 0.5;\r
 677   if((5.0< cntr) && (cntr<10.0)) binct = 1.5;\r
 678   if((10.0< cntr) && (cntr<20.0)) binct = 2.5;\r
 679   if((20.0< cntr) && (cntr<30.0)) binct = 3.5;\r
 680   if((30.0< cntr) && (cntr<40.0)) binct = 4.5;\r
 681   if((40.0< cntr) && (cntr<50.0)) binct = 5.5;\r
 682   if((50.0< cntr) && (cntr<60.0)) binct = 6.5;\r
 683   if((60.0< cntr) && (cntr<70.0)) binct = 7.5;\r
 684   if((70.0< cntr) && (cntr<80.0)) binct = 8.5;\r
 685   if((80.0< cntr) && (cntr<90.0)) binct = 9.5;\r
 686   if((90.0< cntr) && (cntr<100.0)) binct = 10.5;\r
 687   \r
 688   if((0.< cntr) && (cntr < 20.)) binctt = 0.5;\r
 689   if((20.< cntr) && (cntr < 40.)) binctt = 1.5;\r
 690   if((40.< cntr) && (cntr < 80.)) binctt = 2.5;\r
 691 \r
 692   if((0.0< cntr) && (cntr<5.0)) binctMore = 0.5;\r
 693   if((5.0< cntr) && (cntr<10.0)) binctMore = 1.5;\r
 694   if((10.0< cntr) && (cntr<15.0)) binctMore = 2.5;\r
 695   if((15.0< cntr) && (cntr<20.0)) binctMore = 3.5;\r
 696   if((20.0< cntr) && (cntr<25.0)) binctMore = 4.5;\r
 697   if((25.0< cntr) && (cntr<30.0)) binctMore = 5.5;\r
 698   if((30.0< cntr) && (cntr<35.0)) binctMore = 6.5;\r
 699   if((35.0< cntr) && (cntr<40.0)) binctMore = 7.5;\r
 700   if((40.0< cntr) && (cntr<45.0)) binctMore = 8.5;\r
 701   if((45.0< cntr) && (cntr<50.0)) binctMore = 9.5;\r
 702   if((50.0< cntr) && (cntr<55.0)) binctMore = 10.5;\r
 703   if((55.0< cntr) && (cntr<60.0)) binctMore = 11.5;\r
 704   if((60.0< cntr) && (cntr<65.0)) binctMore = 12.5;\r
 705   if((65.0< cntr) && (cntr<70.0)) binctMore = 13.5;\r
 706   if((70.0< cntr) && (cntr<75.0)) binctMore = 14.5;\r
 707   if((75.0< cntr) && (cntr<80.0)) binctMore = 15.5;\r
 708   if((80.0< cntr) && (cntr<85.0)) binctMore = 16.5;\r
 709   if((85.0< cntr) && (cntr<90.0)) binctMore = 17.5;\r
 710   if((90.0< cntr) && (cntr<95.0)) binctMore = 18.5;\r
 711   if((95.0< cntr) && (cntr<100.0)) binctMore = 19.5;\r
 712 \r
 713   binctLess = cntr;\r
 714    \r
 715   \r
 716   if(binct > 11.0) return;\r
 717  \r
 718   // centrality\r
 719   valuensparsea[4] = binct;  \r
 720   valuensparseabis[1] = binct;  \r
 721   valuensparsee[1] = binct;    \r
 722   valuensparsef[1] = binctMore;  \r
 723   valuensparsefsin[1] = binct;  \r
 724   valuensparsefbis[3] = binctMore;  \r
 725   valuensparsefbissin[3] = binct;  \r
 726   valuensparseg[1] = binct;\r
 727   valuensparseh[1] = binct; \r
 728   valuensparsehprofile[1] = binctLess; \r
 729   valuecossinephiep[2] = binctMore;\r
 730 \r
 731   //////////////////////\r
 732   // run number\r
 733   //////////////////////\r
 734 \r
 735   Int_t runnumber = esd->GetRunNumber();\r
 736    \r
 737   if(!fPID->IsInitialized()){\r
 738     // Initialize PID with the given run number\r
 739     fPID->InitializePID(runnumber);\r
 740   }\r
 741 \r
 742 \r
 743   //////////\r
 744   // PID\r
 745   //////////\r
 746  \r
 747   AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();\r
 748   if(!pidResponse){\r
 749     //printf("No PID response set\n");\r
 750     return;\r
 751   }\r
 752   fPID->SetPIDResponse(pidResponse);\r
 753 \r
 754   fHistEV->Fill(binctt,0.0);\r
 755  \r
 756 \r
 757   //////////////////\r
 758   // Event cut\r
 759   //////////////////\r
 760   if(!fHFECuts->CheckEventCuts("fEvRecCuts", esd)) {\r
 761     PostData(1, fListHist);\r
 762     return;\r
 763   }\r
 764 \r
 765   fHistEV->Fill(binctt,1.0);\r
 766 \r
 767   ////////////////////////////////////  \r
 768   // First method event plane\r
 769   ////////////////////////////////////\r
 770 \r
 771   AliEventplane* esdEPa = esd->GetEventplane();\r
 772   Float_t eventPlanea = 0.0;\r
 773   Float_t eventPlaneTPC = 0.0;\r
 774   Float_t eventPlaneV0A = 0.0;\r
 775   Float_t eventPlaneV0C = 0.0;\r
 776   Float_t eventPlaneV0 = 0.0;\r
 777   TVector2 *standardQ = 0x0;\r
 778   TVector2 *qsub1a = 0x0;\r
 779   TVector2 *qsub2a = 0x0;\r
 780 \r
 781   // V0\r
 782 \r
 783   if(fHFEVZEROEventPlane){\r
 784     \r
 785     fHFEVZEROEventPlane->ProcessEvent(esd);\r
 786     \r
 787     if(TMath::Abs(fHFEVZEROEventPlane->GetEventPlaneV0A()+100) < 0.0000001) eventPlaneV0A = -100.0;\r
 788     else {\r
 789       eventPlaneV0A = TVector2::Phi_0_2pi(fHFEVZEROEventPlane->GetEventPlaneV0A());\r
 790       if(eventPlaneV0A > TMath::Pi()) eventPlaneV0A = eventPlaneV0A - TMath::Pi();\r
 791     }\r
 792     \r
 793     if(TMath::Abs(fHFEVZEROEventPlane->GetEventPlaneV0C()+100) < 0.0000001) eventPlaneV0C = -100.0;\r
 794     else {\r
 795       eventPlaneV0C = TVector2::Phi_0_2pi(fHFEVZEROEventPlane->GetEventPlaneV0C());\r
 796       if(eventPlaneV0C > TMath::Pi()) eventPlaneV0C = eventPlaneV0C - TMath::Pi();\r
 797     }\r
 798 \r
 799     if(TMath::Abs(fHFEVZEROEventPlane->GetEventPlaneV0()+100) < 0.0000001) eventPlaneV0 = -100.0;\r
 800     else {\r
 801       eventPlaneV0 = TVector2::Phi_0_2pi(fHFEVZEROEventPlane->GetEventPlaneV0());\r
 802       if(eventPlaneV0 > TMath::Pi()) eventPlaneV0 = eventPlaneV0 - TMath::Pi();\r
 803     }\r
 804     \r
 805   }\r
 806   else {\r
 807     \r
 808     eventPlaneV0 = TVector2::Phi_0_2pi(esdEPa->GetEventplane("V0", esd,2));\r
 809     if(eventPlaneV0 > TMath::Pi()) eventPlaneV0 = eventPlaneV0 - TMath::Pi();\r
 810     eventPlaneV0A = TVector2::Phi_0_2pi(esdEPa->GetEventplane("V0A", esd,2));\r
 811     if(eventPlaneV0A > TMath::Pi()) eventPlaneV0A = eventPlaneV0A - TMath::Pi();\r
 812     eventPlaneV0C = TVector2::Phi_0_2pi(esdEPa->GetEventplane("V0C", esd,2));\r
 813     if(eventPlaneV0C > TMath::Pi()) eventPlaneV0C = eventPlaneV0C - TMath::Pi();\r
 814   \r
 815   }\r
 816 \r
 817   // TPC\r
 818 \r
 819   standardQ = esdEPa->GetQVector(); \r
 820   Double_t qx = -1.0;\r
 821   Double_t qy = -1.0;\r
 822   if(standardQ) {\r
 823     qx = standardQ->X();\r
 824     qy = standardQ->Y();\r
 825   }  \r
 826   TVector2 qVectorfortrack;\r
 827   qVectorfortrack.Set(qx,qy);\r
 828   eventPlaneTPC = TVector2::Phi_0_2pi(qVectorfortrack.Phi())/2.; \r
 829 \r
 830   // Choose the one used for v2\r
 831 \r
 832   if(fVZEROEventPlane) eventPlanea = eventPlaneV0;\r
 833   if(fVZEROEventPlaneA) eventPlanea = eventPlaneV0A;\r
 834   if(fVZEROEventPlaneC) eventPlanea = eventPlaneV0C;\r
 835   if(!fVZEROEventPlane) eventPlanea = eventPlaneTPC;\r
 836 \r
 837   valuecossinephiep[0] = TMath::Cos(2*eventPlanea);\r
 838   valuecossinephiep[1] = TMath::Sin(2*eventPlanea);\r
 839 \r
 840   Float_t eventPlanesub1a = -100.0;\r
 841   Float_t eventPlanesub2a = -100.0;\r
 842   Double_t diffsub1sub2a = -100.0;\r
 843   Double_t diffsub1sub2asin = -100.0;\r
 844   Double_t diffsubasubb = -100.0;\r
 845   Double_t diffsubasubc = -100.0;\r
 846   Double_t diffsubbsubc = -100.0;\r
 847   Double_t diffsubasubbsin = -100.0;\r
 848   Double_t diffsubasubcsin = -100.0;\r
 849   Double_t diffsubbsubcsin = -100.0;\r
 850 \r
 851   //if(fVZEROEventPlane) {\r
 852   diffsubasubb = TMath::Cos(2.*(eventPlaneV0A - eventPlaneV0C));\r
 853   diffsubasubc = TMath::Cos(2.*(eventPlaneV0A - eventPlaneTPC));\r
 854   diffsubbsubc = TMath::Cos(2.*(eventPlaneV0C - eventPlaneTPC));\r
 855 \r
 856   diffsubasubbsin = TMath::Sin(2.*(eventPlaneV0A - eventPlaneV0C));\r
 857   diffsubasubcsin = TMath::Sin(2.*(eventPlaneV0A - eventPlaneTPC));\r
 858   diffsubbsubcsin = TMath::Sin(2.*(eventPlaneV0C - eventPlaneTPC));\r
 859   //}\r
 860   //else {\r
 861   qsub1a = esdEPa->GetQsub1();\r
 862   qsub2a = esdEPa->GetQsub2();\r
 863   if(qsub1a) eventPlanesub1a = TVector2::Phi_0_2pi(qsub1a->Phi())/2.;\r
 864   if(qsub2a) eventPlanesub2a = TVector2::Phi_0_2pi(qsub2a->Phi())/2.;\r
 865   if(qsub1a && qsub2a) {\r
 866     diffsub1sub2a = TMath::Cos(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.));\r
 867     diffsub1sub2asin = TMath::Sin(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.));\r
 868   }\r
 869   //}\r
 870   \r
 871 \r
 872   /////////////////////////////////////////////////////////\r
 873   // Cut for event with event plane reconstructed by all\r
 874   ////////////////////////////////////////////////////////\r
 875   \r
 876   //if(!fVZEROEventPlane) {\r
 877   // if(!standardQ) {\r
 878   //  eventplanedefined = kFALSE;\r
 879       //PostData(1, fListHist);\r
 880       //return;\r
 881   //}\r
 882   //}\r
 883 \r
 884   if((!standardQ) || (!qsub1a) || (!qsub2a)) return;\r
 885 \r
 886   ///////////////////////\r
 887   // AliFlowEvent  \r
 888   //////////////////////\r
 889 \r
 890   Int_t nbtracks = esd->GetNumberOfTracks();\r
 891   //printf("Number of tracks %d\n",nbtracks);\r
 892 \r
 893   fcutsRP->SetEvent( InputEvent(), MCEvent());\r
 894   fcutsPOI->SetEvent( InputEvent(), MCEvent());\r
 895   if( fflowEvent ){ \r
 896     fflowEvent->~AliFlowEvent();\r
 897     new(fflowEvent) AliFlowEvent(fcutsRP,fcutsPOI);\r
 898   }else fflowEvent = new AliFlowEvent(fcutsRP,fcutsPOI);\r
 899   if(mcEvent && mcEvent->GenEventHeader()) {\r
 900     fflowEvent->SetMCReactionPlaneAngle(mcEvent);\r
 901     //if reaction plane not set from elsewhere randomize it before adding flow\r
 902     //if (!fflowEvent->IsSetMCReactionPlaneAngle()) fflowEvent->SetMCReactionPlaneAngle(gRandom->Uniform(0.0,TMath::TwoPi()));\r
 903     mcReactionPlane = TVector2::Phi_0_2pi(fflowEvent->GetMCReactionPlaneAngle());\r
 904     if(mcReactionPlane > TMath::Pi()) mcReactionPlane = mcReactionPlane - TMath::Pi();\r
 905     //printf("MC reaction plane %f\n",mcReactionPlane);\r
 906   }\r
 907   fflowEvent->SetReferenceMultiplicity( nbtracks );\r
 908   fflowEvent->DefineDeadZone(0,0,0,0);\r
 909   //fflowEvent.TagSubeventsInEta(-0.8,-0.1,0.1,0.8);\r
 910 \r
 911   ////////////////\r
 912   // MC\r
 913   ///////////////\r
 914   if(fUseMCReactionPlane) {\r
 915     eventPlanea = mcReactionPlane;\r
 916     diffsub1sub2a = 0.0;\r
 917   }\r
 918 \r
 919   \r
 920   //////////////////////\r
 921   // Fill Histos\r
 922   //////////////////////\r
 923 \r
 924   if(eventplanedefined) {\r
 925     \r
 926     fHistEV->Fill(binctt,2.0);\r
 927     \r
 928     // Fill\r
 929     valuensparsea[0] = eventPlaneV0A;\r
 930     valuensparsea[1] = eventPlaneV0C;\r
 931     valuensparsea[2] = eventPlaneTPC;\r
 932     valuensparsea[3] = eventPlaneV0;  \r
 933     fEventPlane->Fill(&valuensparsea[0]);\r
 934 \r
 935     // Fill\r
 936     fCosSin2phiep->Fill(&valuecossinephiep[0]);\r
 937     \r
 938     if(!fVZEROEventPlane) {\r
 939       valuensparsef[0] = diffsub1sub2a;\r
 940       fCosRes->Fill(&valuensparsef[0]);\r
 941       valuensparsefsin[0] = diffsub1sub2asin;\r
 942       fSinRes->Fill(&valuensparsefsin[0]);\r
 943       if(fDebugLevel > 1) {\r
 944         fProfileCosRes->Fill(valuensparsef[1],valuensparsef[0]);\r
 945       }\r
 946     }\r
 947     else {\r
 948       valuensparsefbis[0] = diffsubasubb;\r
 949       valuensparsefbis[1] = diffsubbsubc;\r
 950       valuensparsefbis[2] = diffsubasubc;\r
 951       fCosResabc->Fill(&valuensparsefbis[0]);\r
 952       valuensparsefbissin[0] = diffsubasubbsin;\r
 953       valuensparsefbissin[1] = diffsubbsubcsin;\r
 954       valuensparsefbissin[2] = diffsubasubcsin;\r
 955       fSinResabc->Fill(&valuensparsefbissin[0]);\r
 956       if(fDebugLevel > 1) {\r
 957         fProfileCosResab->Fill(valuensparsefbis[3],valuensparsefbis[0]);\r
 958         fProfileCosResac->Fill(valuensparsefbis[3],valuensparsefbis[1]);\r
 959         fProfileCosResbc->Fill(valuensparsefbis[3],valuensparsefbis[2]);\r
 960       }\r
 961     }\r
 962     \r
 963   }\r
 964   \r
 965   //////////////////////////\r
 966   // Loop over ESD track\r
 967   //////////////////////////\r
 968  \r
 969 \r
 970   for(Int_t k = 0; k < nbtracks; k++){\r
 971     \r
 972     AliESDtrack *track = esd->GetTrack(k);\r
 973     if(!track) continue;\r
 974 \r
 975     Bool_t survived = kTRUE;\r
 976     for(Int_t icut = AliHFEcuts::kStepRecKineITSTPC; icut <= AliHFEcuts::kStepHFEcutsTRD; icut++){\r
 977       if(!fHFECuts->CheckParticleCuts(icut + AliHFEcuts::kNcutStepsMCTrack, (TObject *)track)){\r
 978         survived = kFALSE;\r
 979         break;\r
 980       }\r
 981       //printf("Pass the cut %d\n",icut);\r
 982     }\r
 983     \r
 984     if(!survived) continue;\r
 985 \r
 986     // Apply PID\r
 987     if(!fNoPID) {\r
 988       // Apply PID for Data\r
 989       if(!fMCPID) {\r
 990         AliHFEpidObject hfetrack;\r
 991         hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);\r
 992         hfetrack.SetRecTrack(track);\r
 993         hfetrack.SetCentrality((Int_t)binct);\r
 994         //printf("centrality %f and %d\n",binct,hfetrack.GetCentrality());\r
 995         hfetrack.SetPbPb();\r
 996         if(!fPID->IsSelected(&hfetrack,0x0,"",fPIDqa)) {\r
 997           continue;\r
 998         }\r
 999       }\r
1000       else {\r
1001         if(!mcEvent) continue;\r
1002         if(!(mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel()))))) continue;\r
1003         //printf("PdgCode %d\n",TMath::Abs(mctrack->Particle()->GetPdgCode()));\r
1004         if(TMath::Abs(mctrack->Particle()->GetPdgCode())!=11) continue;\r
1005       }\r
1006     }\r
1007 \r
1008 \r
1009     /////////////////////////////////////////////////////////////////////////////\r
1010     // Add candidate to AliFlowEvent for POI and subtract from RP if needed\r
1011     ////////////////////////////////////////////////////////////////////////////\r
1012     Int_t idtrack = static_cast<AliVTrack*>(track)->GetID();\r
1013     Bool_t found = kFALSE;\r
1014     Int_t numberoffound = 0;\r
1015     //printf("A: Number of tracks %d\n",fflowEvent->NumberOfTracks());\r
1016     for(Int_t iRPs=0; iRPs< fflowEvent->NumberOfTracks(); iRPs++) {\r
1017       AliFlowTrack *iRP = (AliFlowTrack*) (fflowEvent->GetTrack(iRPs));\r
1018       //if(!iRP->InRPSelection()) continue;\r
1019       if( TMath::Abs(idtrack) == TMath::Abs(iRP->GetID()) ) {\r
1020         iRP->SetForPOISelection(kTRUE);\r
1021         found = kTRUE;\r
1022         numberoffound ++;\r
1023       }\r
1024     }\r
1025     //printf("Found %d mal\n",numberoffound);\r
1026     if(!found) {\r
1027       AliFlowCandidateTrack *sTrack = (AliFlowCandidateTrack*) MakeTrack(massElectron,track->Pt(),track->Phi(), track->Eta());\r
1028       sTrack->SetID(idtrack);\r
1029       fflowEvent->AddTrack(sTrack);\r
1030       //printf("Add the track\n");\r
1031     }\r
1032     //printf("B: Number of tracks %d\n",fflowEvent->NumberOfTracks());\r
1033     \r
1034     \r
1035     /////////////////////////////////////////////////////////\r
1036     // Subtract electron candidate from TPC event plane\r
1037     ////////////////////////////////////////////////////////\r
1038     Float_t eventplanesubtracted = 0.0;    \r
1039 \r
1040     //if(eventplanedefined && (!fVZEROEventPlane)) {\r
1041     if(!fVZEROEventPlane) {\r
1042       // Subtract the tracks from the event plane\r
1043       Double_t qX = standardQ->X() - esdEPa->GetQContributionX(track);  //Modify the components: subtract the track you want to look at with your analysis\r
1044       Double_t qY = standardQ->Y() - esdEPa->GetQContributionY(track);  //Modify the components: subtract the track you want to look at with your analysis\r
1045       TVector2 newQVectorfortrack;\r
1046       newQVectorfortrack.Set(qX,qY);\r
1047       eventplanesubtracted = TVector2::Phi_0_2pi(newQVectorfortrack.Phi())/2; \r
1048     }\r
1049     else eventplanesubtracted = eventPlanea;\r
1050 \r
1051     ////////////////////////////////////////\r
1052     // Fill pt and eta for the THnSparseF\r
1053     ///////////////////////////////////////\r
1054 \r
1055     valuensparsee[2] = track->Pt();\r
1056     valuensparsee[3] = track->Eta();    \r
1057     valuensparseg[2] = track->Pt();\r
1058     valuensparseh[2] = track->Pt();\r
1059     valuensparsehprofile[2] = track->Pt();\r
1060 \r
1061     Bool_t fillEventPlane = kTRUE;\r
1062     if(!fVZEROEventPlane){\r
1063       //if((!qsub1a) || (!qsub2a) || (!eventplanedefined)) fillEventPlane = kFALSE;\r
1064       if((!qsub1a) || (!qsub2a)) fillEventPlane = kFALSE;\r
1065       if(fSubEtaGapTPC) {\r
1066         if(track->Eta() < (- fEtaGap/2.)) eventplanesubtracted = eventPlanesub1a;\r
1067         else if(track->Eta() > (fEtaGap/2.)) eventplanesubtracted = eventPlanesub2a;\r
1068         else fillEventPlane = kFALSE;\r
1069       }\r
1070     }\r
1071     \r
1072     \r
1073     ///////////////\r
1074     // MC\r
1075     //////////////\r
1076     if(fUseMCReactionPlane) {\r
1077       eventplanesubtracted = mcReactionPlane;\r
1078       fillEventPlane = kTRUE;\r
1079     }\r
1080     \r
1081     //////////////////////////////////////////////////////////////////////////////\r
1082     ///////////////////////////AFTERBURNER\r
1083     Double_t phitrack = track->Phi();    \r
1084     if (fAfterBurnerOn)\r
1085       {\r
1086         phitrack = GetPhiAfterAddV2(track->Phi(),mcReactionPlane);\r
1087       }\r
1088     //////////////////////////////////////////////////////////////////////////////\r
1089 \r
1090 \r
1091     ///////////////////////\r
1092     // Calculate deltaphi\r
1093     ///////////////////////\r
1094     \r
1095     // Suppose phi track is between 0.0 and phi\r
1096     Double_t deltaphi = TVector2::Phi_0_2pi(phitrack - eventplanesubtracted);\r
1097     if(deltaphi > TMath::Pi()) deltaphi = deltaphi - TMath::Pi();\r
1098    \r
1099     /////////////////////\r
1100     // Fill THnSparseF\r
1101     /////////////////////\r
1102 \r
1103     //\r
1104     valuensparseabis[0] = eventplanesubtracted;\r
1105     if(fillEventPlane) fEventPlaneaftersubtraction->Fill(&valuensparseabis[0]);\r
1106     \r
1107 \r
1108     if(fDebugLevel > 3) \r
1109       {\r
1110         //\r
1111         valuensparsee[0] = TMath::Cos(2*phitrack);\r
1112         fCos2phie->Fill(&valuensparsee[0]);\r
1113         valuensparsee[0] = TMath::Sin(2*phitrack);\r
1114         fSin2phie->Fill(&valuensparsee[0]);\r
1115         //\r
1116         valuensparsee[0] = TMath::Cos(2*eventplanesubtracted);\r
1117         if(fillEventPlane) fCos2phiep->Fill(&valuensparsee[0]);\r
1118         valuensparsee[0] = TMath::Sin(2*eventplanesubtracted);\r
1119         if(fillEventPlane) fSin2phiep->Fill(&valuensparsee[0]);\r
1120         valuensparsee[0] = TMath::Sin(2*TVector2::Phi_mpi_pi(phitrack-eventplanesubtracted));\r
1121         if(fillEventPlane) fSin2phiephiep->Fill(&valuensparsee[0]);\r
1122         //\r
1123       }\r
1124 \r
1125     // \r
1126     valuensparseg[0] = deltaphi;\r
1127     if(fillEventPlane) fDeltaPhiMaps->Fill(&valuensparseg[0]);\r
1128     \r
1129     //\r
1130     valuensparseh[0] = TMath::Cos(2*TVector2::Phi_mpi_pi(phitrack-eventplanesubtracted));\r
1131     if(fillEventPlane) {\r
1132       fCosPhiMaps->Fill(&valuensparseh[0]);\r
1133       if(fDebugLevel > 0) {\r
1134         fProfileCosPhiMaps->Fill(valuensparsehprofile[1],valuensparsehprofile[2],valuensparseh[0]);\r
1135       }\r
1136     }\r
1137     \r
1138     \r
1139     \r
1140   }\r
1141 \r
1142   //////////////////////////////////////////////////////////////////////////////\r
1143   ///////////////////////////AFTERBURNER\r
1144   if (fAfterBurnerOn)\r
1145     {\r
1146       fflowEvent->AddFlow(fV1,fV2,fV3,fV4,fV5);     //add flow\r
1147       fflowEvent->CloneTracks(fNonFlowNumberOfTrackClones); //add nonflow by cloning tracks\r
1148     }\r
1149   //////////////////////////////////////////////////////////////////////////////\r
1150 \r
1151 \r
1152 \r
1153   for(Int_t bincless = 0; bincless < fNbBinsCentralityQCumulant; bincless++) {\r
1154     if((fBinCentralityLess[bincless]< cntr) && (cntr < fBinCentralityLess[bincless+1])) PostData(bincless+2,fflowEvent);\r
1155   }\r
1156   \r
1157   PostData(1, fListHist);\r
1158 \r
1159 \r
1160  \r
1161 }\r
1162 //______________________________________________________________________________\r
1163 AliFlowCandidateTrack *AliAnalysisTaskHFEFlow::MakeTrack( Double_t mass, \r
1164                           Double_t pt, Double_t phi, Double_t eta) {\r
1165   //\r
1166   //  Make Track (Not needed actually)\r
1167   //\r
1168 \r
1169   AliFlowCandidateTrack *sTrack = new AliFlowCandidateTrack();\r
1170   sTrack->SetMass(mass);\r
1171   sTrack->SetPt(pt);\r
1172   sTrack->SetPhi(phi);\r
1173   sTrack->SetEta(eta);\r
1174   sTrack->SetForPOISelection(kTRUE);\r
1175   sTrack->SetForRPSelection(kFALSE);\r
1176   return sTrack;\r
1177 }\r
1178 //_________________________________________________________________________________ \r
1179 Double_t AliAnalysisTaskHFEFlow::GetPhiAfterAddV2(Double_t phi,Double_t reactionPlaneAngle) const\r
1180 {\r
1181   //\r
1182   // Adds v2, uses Newton-Raphson iteration\r
1183   //\r
1184   Double_t phiend=phi;\r
1185   Double_t phi0=phi;\r
1186   Double_t f=0.;\r
1187   Double_t fp=0.;\r
1188   Double_t phiprev=0.;\r
1189 \r
1190   for (Int_t i=0; i<fMaxNumberOfIterations; i++)\r
1191   {\r
1192     phiprev=phiend; //store last value for comparison\r
1193     f =  phiend-phi0+fV2*TMath::Sin(2.*(phiend-reactionPlaneAngle));\r
1194     fp = 1.0+2.0*fV2*TMath::Cos(2.*(phiend-reactionPlaneAngle)); //first derivative\r
1195     phiend -= f/fp;\r
1196     if (TMath::AreEqualAbs(phiprev,phiend,fPrecisionPhi)) break;\r
1197   }\r
1198   return phiend;\r
1199 }\r