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 \r
  30 #include "AliAnalysisTaskSE.h"\r
  31 \r
  32 #include "AliESDInputHandler.h"\r
  33 #include "AliMCEvent.h"\r
  34 #include "AliESD.h"\r
  35 #include "AliESDEvent.h"\r
  36 #include "AliPIDResponse.h"\r
  37 #include "AliESDVZERO.h"\r
  38 #include "AliESDUtils.h"\r
  39 #include "AliMCParticle.h"\r
  40 \r
  41 #include "AliFlowCandidateTrack.h"\r
  42 #include "AliFlowEvent.h"\r
  43 #include "AliFlowTrackCuts.h"\r
  44 #include "AliFlowVector.h"\r
  45 #include "AliFlowCommonConstants.h"\r
  46 \r
  47 \r
  48 #include "AliHFEcuts.h"\r
  49 #include "AliHFEpid.h"\r
  50 #include "AliHFEpidQAmanager.h"\r
  51 #include "AliHFEtools.h"\r
  52 \r
  53 \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(), \r
  64   fSubEtaGapTPC(kFALSE),\r
  65   fEtaGap(0.0),\r
  66   fNbBinsPtQCumulant(15),\r
  67   fMinPtQCumulant(0.0),\r
  68   fMaxPtQCumulant(6.0),\r
  69   fAfterBurnerOn(kFALSE),\r
  70   fNonFlowNumberOfTrackClones(0),\r
  71   fV1(0.),\r
  72   fV2(0.),\r
  73   fV3(0.),\r
  74   fV4(0.),\r
  75   fV5(0.),\r
  76   fMaxNumberOfIterations(100),\r
  77   fPrecisionPhi(0.001),\r
  78   fUseMCReactionPlane(kFALSE),\r
  79   fMCPID(kFALSE),\r
  80   fcutsRP(0),\r
  81   fcutsPOI(0),\r
  82   fHFECuts(0),\r
  83   fPID(0),\r
  84   fPIDqa(0),\r
  85   fflowEvent(0x0),\r
  86   fHistEV(0),\r
  87   fEventPlane(0x0),\r
  88   fEventPlaneaftersubtraction(0x0),\r
  89   fCos2phie(0x0),\r
  90   fSin2phie(0x0),\r
  91   fCos2phiep(0x0),\r
  92   fSin2phiep(0x0),\r
  93   fSin2phiephiep(0x0),\r
  94   fCosRes(0x0),\r
  95   fDeltaPhiMaps(0x0),\r
  96   fCosPhiMaps(0x0)\r
  97 {\r
  98   // Constructor\r
  99 \r
 100   \r
 101   \r
 102 }\r
 103 //______________________________________________________________________________\r
 104 AliAnalysisTaskHFEFlow:: AliAnalysisTaskHFEFlow(const char *name) :\r
 105   AliAnalysisTaskSE(name),\r
 106   fListHist(), \r
 107   fSubEtaGapTPC(kFALSE),\r
 108   fEtaGap(0.0),\r
 109   fNbBinsPtQCumulant(15),\r
 110   fMinPtQCumulant(0.0),\r
 111   fMaxPtQCumulant(6.0),\r
 112   fAfterBurnerOn(kFALSE),\r
 113   fNonFlowNumberOfTrackClones(0),\r
 114   fV1(0.),\r
 115   fV2(0.),\r
 116   fV3(0.),\r
 117   fV4(0.),\r
 118   fV5(0.),\r
 119   fMaxNumberOfIterations(100),\r
 120   fPrecisionPhi(0.001),\r
 121   fUseMCReactionPlane(kFALSE),\r
 122   fMCPID(kFALSE),\r
 123   fcutsRP(0),\r
 124   fcutsPOI(0),\r
 125   fHFECuts(0),\r
 126   fPID(0),\r
 127   fPIDqa(0),\r
 128   fflowEvent(0x0),\r
 129   fHistEV(0),\r
 130   fEventPlane(0x0),\r
 131   fEventPlaneaftersubtraction(0x0),\r
 132   fCos2phie(0x0),\r
 133   fSin2phie(0x0),\r
 134   fCos2phiep(0x0),\r
 135   fSin2phiep(0x0),\r
 136   fSin2phiephiep(0x0),\r
 137   fCosRes(0x0),\r
 138   fDeltaPhiMaps(0x0),\r
 139   fCosPhiMaps(0x0)\r
 140 {\r
 141   //\r
 142   // named ctor\r
 143   //\r
 144 \r
 145   fPID = new AliHFEpid("hfePid");\r
 146   fPIDqa = new AliHFEpidQAmanager;\r
 147 \r
 148   DefineInput(0,TChain::Class());\r
 149   DefineOutput(1, TList::Class());\r
 150   DefineOutput(2,AliFlowEventSimple::Class()); // first band 0-20%\r
 151   DefineOutput(3,AliFlowEventSimple::Class()); // first band 20-40%\r
 152   DefineOutput(4,AliFlowEventSimple::Class()); // first band 40-60%\r
 153   DefineOutput(5,AliFlowEventSimple::Class()); // first band 60-80%\r
 154   \r
 155 }\r
 156 //________________________________________________________________________\r
 157 void AliAnalysisTaskHFEFlow::UserCreateOutputObjects()\r
 158 {\r
 159 \r
 160   //********************\r
 161   // Create histograms\r
 162   //********************\r
 163 \r
 164   //**************\r
 165   // Cuts\r
 166   //**************\r
 167 \r
 168   //---------Data selection----------\r
 169   //kMC, kGlobal, kTPCstandalone, kSPDtracklet, kPMD\r
 170   //AliFlowTrackCuts::trackParameterType rptype = AliFlowTrackCuts::kGlobal;\r
 171   //AliFlowTrackCuts::trackParameterType poitype = AliFlowTrackCuts::kGlobal;\r
 172 \r
 173   //---------Parameter mixing--------\r
 174   //kPure - no mixing, kTrackWithMCkine, kTrackWithMCPID, kTrackWithMCpt\r
 175   //AliFlowTrackCuts::trackParameterMix rpmix = AliFlowTrackCuts::kPure;\r
 176   //AliFlowTrackCuts::trackParameterMix poimix = AliFlowTrackCuts::kPure;\r
 177 \r
 178   // RP TRACK CUTS:\r
 179   fcutsRP = AliFlowTrackCuts::GetStandardTPCStandaloneTrackCuts2010();\r
 180   fcutsRP->SetName("StandartTPC");\r
 181   fcutsRP->SetEtaRange(-0.9,0.9);\r
 182   fcutsRP->SetQA(kTRUE);\r
 183 \r
 184   //POI TRACK CUTS:\r
 185   fcutsPOI = new AliFlowTrackCuts("dummy");\r
 186   fcutsPOI->SetParamType(AliFlowTrackCuts::kGlobal);\r
 187   fcutsPOI->SetPtRange(+1,-1); // select nothing QUICK\r
 188   fcutsPOI->SetEtaRange(+1,-1); // select nothing VZERO\r
 189 \r
 190   // Flow\r
 191   AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();\r
 192   cc->SetNbinsMult(10000);\r
 193   cc->SetMultMin(0);\r
 194   cc->SetMultMax(10000.);\r
 195   cc->SetNbinsPt(fNbBinsPtQCumulant);\r
 196   cc->SetPtMin(fMinPtQCumulant);\r
 197   cc->SetPtMax(fMaxPtQCumulant);\r
 198   cc->SetNbinsPhi(180);\r
 199   cc->SetPhiMin(0.0);\r
 200   cc->SetPhiMax(TMath::TwoPi());\r
 201   cc->SetNbinsEta(200);\r
 202   cc->SetEtaMin(-0.9);\r
 203   cc->SetEtaMax(+0.9);\r
 204   cc->SetNbinsQ(500);\r
 205   cc->SetQMin(0.0);\r
 206   cc->SetQMax(3.0);\r
 207 \r
 208   \r
 209   // HFE cuts\r
 210 \r
 211   if(!fHFECuts){\r
 212     fHFECuts = new AliHFEcuts;\r
 213     fHFECuts->CreateStandardCuts();\r
 214   }\r
 215   fHFECuts->Initialize();\r
 216   \r
 217   // PID HFE\r
 218   //fPID->SetHasMCData(HasMCData());\r
 219   if(!fPID->GetNumberOfPIDdetectors()) fPID->AddDetector("TPC", 0);\r
 220   fPID->InitializePID();\r
 221   fPIDqa->Initialize(fPID);\r
 222   fPID->SortDetectors();\r
 223   \r
 224   //**************************\r
 225   // Bins for the THnSparse\r
 226   //**************************\r
 227 \r
 228   Int_t nBinsPt = 25;\r
 229   Double_t minPt = 0.001;\r
 230   Double_t maxPt = 10.0;\r
 231   Double_t binLimLogPt[nBinsPt+1];\r
 232   Double_t binLimPt[nBinsPt+1];\r
 233   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
 234   for(Int_t i=0; i<=nBinsPt; i++) binLimPt[i]=(Double_t)TMath::Power(10,binLimLogPt[i]);\r
 235 \r
 236   Int_t nBinsEta = 8;\r
 237   Double_t minEta = -0.8;\r
 238   Double_t maxEta = 0.8;\r
 239   Double_t binLimEta[nBinsEta+1];\r
 240   for(Int_t i=0; i<=nBinsEta; i++) binLimEta[i]=(Double_t)minEta + (maxEta-minEta)/nBinsEta*(Double_t)i ;\r
 241  \r
 242   Int_t nBinsCos = 50;\r
 243   Double_t minCos = -1.0;\r
 244   Double_t maxCos = 1.0;\r
 245   Double_t binLimCos[nBinsCos+1];\r
 246   for(Int_t i=0; i<=nBinsCos; i++) binLimCos[i]=(Double_t)minCos + (maxCos-minCos)/nBinsCos*(Double_t)i ;\r
 247  \r
 248 \r
 249   Int_t nBinsC = 11;\r
 250   Double_t minC = 0.0;\r
 251   Double_t maxC = 11.0;\r
 252   Double_t binLimC[nBinsC+1];\r
 253   for(Int_t i=0; i<=nBinsC; i++) binLimC[i]=(Double_t)minC + (maxC-minC)/nBinsC*(Double_t)i ;\r
 254  \r
 255   Int_t nBinsPhi = 50;\r
 256   Double_t minPhi = 0.0;\r
 257   Double_t maxPhi = TMath::Pi();\r
 258   Double_t binLimPhi[nBinsPhi+1];\r
 259   for(Int_t i=0; i<=nBinsPhi; i++) {\r
 260     binLimPhi[i]=(Double_t)minPhi + (maxPhi-minPhi)/nBinsPhi*(Double_t)i ;\r
 261     //printf("bin phi is %f for %d\n",binLimPhi[i],i);\r
 262   }\r
 263  \r
 264   Int_t nBinsQ = 200;\r
 265   Double_t minQ = -300.0;\r
 266   Double_t maxQ = 300.0;\r
 267   Double_t binLimQ[nBinsQ+1];\r
 268   for(Int_t i=0; i<=nBinsQ; i++) binLimQ[i]=(Double_t)minQ + (maxQ-minQ)/nBinsQ*(Double_t)i ;\r
 269 \r
 270   Int_t nBinsM = 250;\r
 271   Double_t minM = 0.0;\r
 272   Double_t maxM = 800.0;\r
 273   Double_t binLimM[nBinsM+1];\r
 274   for(Int_t i=0; i<=nBinsM; i++) binLimM[i]=(Double_t)minM + (maxM-minM)/nBinsM*(Double_t)i ;\r
 275 \r
 276   Int_t nBinsCELL = 64;\r
 277   Double_t minCELL = 0.0;\r
 278   Double_t maxCELL = 64.0;\r
 279   Double_t binLimCELL[nBinsCELL+1];\r
 280   for(Int_t i=0; i<=nBinsCELL; i++) binLimCELL[i]=(Double_t)minCELL + (maxCELL-minCELL)/nBinsCELL*(Double_t)i ;\r
 281   \r
 282 \r
 283   //******************\r
 284   // Histograms\r
 285   //******************\r
 286     \r
 287   fListHist = new TList();\r
 288 \r
 289   // Histos\r
 290   fHistEV = new TH2D("fHistEV", "events", 3, 0, 3, 3, 0,3);\r
 291   \r
 292   // Event plane as function of phiep, centrality\r
 293   const Int_t nDima=2;\r
 294   Int_t nBina[nDima] = {nBinsPhi, nBinsC};\r
 295   fEventPlane = new THnSparseF("EventPlane","EventPlane",nDima,nBina);\r
 296   fEventPlane->SetBinEdges(0,binLimPhi);\r
 297   fEventPlane->SetBinEdges(1,binLimC);\r
 298   \r
 299   // Event Plane after subtraction as function of phiep, centrality, pt, eta\r
 300   const Int_t nDimb=2;\r
 301   Int_t nBinb[nDimb] = {nBinsPhi, nBinsC};\r
 302   fEventPlaneaftersubtraction = new THnSparseF("EventPlane_aftersubtraction","EventPlane_aftersubtraction",nDimb,nBinb);\r
 303   fEventPlaneaftersubtraction->SetBinEdges(0,binLimPhi);\r
 304   fEventPlaneaftersubtraction->SetBinEdges(1,binLimC);\r
 305 \r
 306   // Monitoring Event plane after subtraction of the track\r
 307   const Int_t nDime=4;\r
 308   Int_t nBine[nDime] = {nBinsCos, nBinsC, nBinsPt, nBinsEta};\r
 309   fCos2phie = new THnSparseF("cos2phie","cos2phie",nDime,nBine);\r
 310   fCos2phie->SetBinEdges(2,binLimPt);\r
 311   fCos2phie->SetBinEdges(3,binLimEta);\r
 312   fCos2phie->SetBinEdges(0,binLimCos);\r
 313   fCos2phie->SetBinEdges(1,binLimC);\r
 314   fSin2phie = new THnSparseF("sin2phie","sin2phie",nDime,nBine);\r
 315   fSin2phie->SetBinEdges(2,binLimPt);\r
 316   fSin2phie->SetBinEdges(3,binLimEta);\r
 317   fSin2phie->SetBinEdges(0,binLimCos);\r
 318   fSin2phie->SetBinEdges(1,binLimC);\r
 319   fCos2phiep = new THnSparseF("cos2phiep","cos2phiep",nDime,nBine);\r
 320   fCos2phiep->SetBinEdges(2,binLimPt);\r
 321   fCos2phiep->SetBinEdges(3,binLimEta);\r
 322   fCos2phiep->SetBinEdges(0,binLimCos);\r
 323   fCos2phiep->SetBinEdges(1,binLimC);\r
 324   fSin2phiep = new THnSparseF("sin2phiep","sin2phiep",nDime,nBine);\r
 325   fSin2phiep->SetBinEdges(2,binLimPt);\r
 326   fSin2phiep->SetBinEdges(3,binLimEta);\r
 327   fSin2phiep->SetBinEdges(0,binLimCos);\r
 328   fSin2phiep->SetBinEdges(1,binLimC);\r
 329   fSin2phiephiep = new THnSparseF("sin2phie_phiep","sin2phie_phiep",nDime,nBine);\r
 330   fSin2phiephiep->SetBinEdges(2,binLimPt);\r
 331   fSin2phiephiep->SetBinEdges(3,binLimEta);\r
 332   fSin2phiephiep->SetBinEdges(0,binLimCos);\r
 333   fSin2phiephiep->SetBinEdges(1,binLimC);\r
 334   \r
 335   // Resolution cosres centrality\r
 336   const Int_t nDimf=2;\r
 337   Int_t nBinf[nDimf] = {nBinsCos, nBinsC};\r
 338   fCosRes = new THnSparseF("CosRes","CosRes",nDimf,nBinf);\r
 339   fCosRes->SetBinEdges(0,binLimCos);\r
 340   fCosRes->SetBinEdges(1,binLimC);\r
 341   \r
 342   // Maps delta phi\r
 343   const Int_t nDimg=3;\r
 344   Int_t nBing[nDimg] = {nBinsPhi,nBinsC,nBinsPt};\r
 345   fDeltaPhiMaps = new THnSparseF("DeltaPhiMaps","DeltaPhiMaps",nDimg,nBing);\r
 346   fDeltaPhiMaps->SetBinEdges(0,binLimPhi);\r
 347   fDeltaPhiMaps->SetBinEdges(1,binLimC);\r
 348   fDeltaPhiMaps->SetBinEdges(2,binLimPt);\r
 349   \r
 350   // Maps cos phi\r
 351   const Int_t nDimh=3;\r
 352   Int_t nBinh[nDimh] = {nBinsCos,nBinsC,nBinsPt};\r
 353   fCosPhiMaps = new THnSparseF("CosPhiMaps","CosPhiMaps",nDimh,nBinh);\r
 354   fCosPhiMaps->SetBinEdges(0,binLimCos);\r
 355   fCosPhiMaps->SetBinEdges(1,binLimC);\r
 356   fCosPhiMaps->SetBinEdges(2,binLimPt);\r
 357   \r
 358   //**************************\r
 359   // Add to the list\r
 360   //******************************\r
 361 \r
 362   fListHist->Add(fPIDqa->MakeList("HFEpidQA"));\r
 363         \r
 364   fListHist->Add(fHistEV);\r
 365 \r
 366   fListHist->Add(fEventPlane);\r
 367   fListHist->Add(fEventPlaneaftersubtraction);\r
 368   fListHist->Add(fCos2phie);\r
 369   fListHist->Add(fSin2phie);\r
 370   fListHist->Add(fCos2phiep);\r
 371   fListHist->Add(fSin2phiep);\r
 372   fListHist->Add(fSin2phiephiep);\r
 373   fListHist->Add(fCosRes);\r
 374   fListHist->Add(fDeltaPhiMaps);\r
 375   fListHist->Add(fCosPhiMaps);\r
 376   \r
 377 \r
 378   PostData(1, fListHist);\r
 379 \r
 380 \r
 381 }\r
 382    \r
 383 //________________________________________________________________________\r
 384 void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/)\r
 385 {\r
 386   //\r
 387   // Loop over event\r
 388   //\r
 389    \r
 390   Double_t massElectron = 0.000511;\r
 391   Double_t mcReactionPlane = 0.0;\r
 392 \r
 393   Float_t cntr = 0.0;\r
 394   Double_t binct = 11.5;\r
 395   Float_t binctt = -1.0;\r
 396   \r
 397   Double_t valuensparsea[2];\r
 398   Double_t valuensparsee[4];\r
 399   Double_t valuensparsef[2];\r
 400   Double_t valuensparseg[3];\r
 401   Double_t valuensparseh[3];\r
 402 \r
 403   AliMCEvent *mcEvent = MCEvent();\r
 404   AliMCParticle *mctrack = NULL;\r
 405   \r
 406   /////////////////\r
 407   // centrality\r
 408   /////////////////\r
 409   \r
 410   AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());\r
 411   if(!esd) return;\r
 412   AliCentrality *centrality = esd->GetCentrality();\r
 413   if(!centrality) return;\r
 414   cntr = centrality->GetCentralityPercentile("V0M");\r
 415   if((0.0< cntr) && (cntr<5.0)) binct = 0.5;\r
 416   if((5.0< cntr) && (cntr<10.0)) binct = 1.5;\r
 417   if((10.0< cntr) && (cntr<20.0)) binct = 2.5;\r
 418   if((20.0< cntr) && (cntr<30.0)) binct = 3.5;\r
 419   if((30.0< cntr) && (cntr<40.0)) binct = 4.5;\r
 420   if((40.0< cntr) && (cntr<50.0)) binct = 5.5;\r
 421   if((50.0< cntr) && (cntr<60.0)) binct = 6.5;\r
 422   if((60.0< cntr) && (cntr<70.0)) binct = 7.5;\r
 423   if((70.0< cntr) && (cntr<80.0)) binct = 8.5;\r
 424   if((80.0< cntr) && (cntr<90.0)) binct = 9.5;\r
 425   if((90.0< cntr) && (cntr<100.0)) binct = 10.5;\r
 426   \r
 427   if((0.< cntr) && (cntr < 20.)) binctt = 0.5;\r
 428   if((20.< cntr) && (cntr < 40.)) binctt = 1.5;\r
 429   if((40.< cntr) && (cntr < 80.)) binctt = 2.5;\r
 430 \r
 431   if(binct > 11.0) return;\r
 432  \r
 433   // centrality\r
 434   valuensparsea[1] = binct;  \r
 435   valuensparsee[1] = binct;    \r
 436   valuensparsef[1] = binct;  \r
 437   valuensparseg[1] = binct;\r
 438   valuensparseh[1] = binct; \r
 439   \r
 440 \r
 441   //////////////////////\r
 442   // run number\r
 443   //////////////////////\r
 444 \r
 445   Int_t runnumber = esd->GetRunNumber();\r
 446    \r
 447   if(!fPID->IsInitialized()){\r
 448     // Initialize PID with the given run number\r
 449     fPID->InitializePID(runnumber);\r
 450   }\r
 451 \r
 452 \r
 453   //////////\r
 454   // PID\r
 455   //////////\r
 456  \r
 457   AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();\r
 458   if(!pidResponse){\r
 459     //printf("No PID response set\n");\r
 460     return;\r
 461   }\r
 462   fPID->SetPIDResponse(pidResponse);\r
 463 \r
 464   fHistEV->Fill(binctt,0.0);\r
 465  \r
 466 \r
 467   //////////////////\r
 468   // Event cut\r
 469   //////////////////\r
 470   if(!fHFECuts->CheckEventCuts("fEvRecCuts", esd)) {\r
 471     PostData(1, fListHist);\r
 472     return;\r
 473   }\r
 474 \r
 475   fHistEV->Fill(binctt,1.0);\r
 476 \r
 477   ////////////////////////////////////  \r
 478   // First method TPC event plane\r
 479   ////////////////////////////////////\r
 480 \r
 481   AliEventplane* esdEPa = esd->GetEventplane();\r
 482   TVector2 *standardQ = esdEPa->GetQVector();  // This is the "standard" Q-Vector\r
 483   Double_t qx = -1.0;\r
 484   Double_t qy = -1.0;\r
 485   if(standardQ) {\r
 486     qx = standardQ->X();\r
 487     qy = standardQ->Y();\r
 488   }  \r
 489   TVector2 qVectorfortrack;\r
 490   qVectorfortrack.Set(qx,qy);\r
 491   Float_t eventPlanea = TVector2::Phi_0_2pi(qVectorfortrack.Phi())/2.;   \r
 492   \r
 493   TVector2 *qsub1a = esdEPa->GetQsub1();\r
 494   TVector2 *qsub2a = esdEPa->GetQsub2();\r
 495   Float_t eventPlanesub1a = -100.0;\r
 496   Float_t eventPlanesub2a = -100.0;\r
 497   if(qsub1a) eventPlanesub1a = TVector2::Phi_0_2pi(qsub1a->Phi())/2.;\r
 498   if(qsub2a) eventPlanesub2a = TVector2::Phi_0_2pi(qsub2a->Phi())/2.;\r
 499   Double_t diffsub1sub2a = -100.0;\r
 500   if(qsub1a && qsub2a) {\r
 501     diffsub1sub2a = TMath::Cos(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.));\r
 502   }\r
 503 \r
 504   //////////////////////////////////////////\r
 505   // Cut for event with TPC event defined\r
 506   /////////////////////////////////////////\r
 507   \r
 508   if(!standardQ) {\r
 509     PostData(1, fListHist);\r
 510     return;\r
 511   }\r
 512 \r
 513   ///////////////////////\r
 514   // AliFlowEvent  \r
 515   //////////////////////\r
 516 \r
 517   Int_t nbtracks = esd->GetNumberOfTracks();\r
 518   //printf("Number of tracks %d\n",nbtracks);\r
 519 \r
 520   fcutsRP->SetEvent( InputEvent(), MCEvent());\r
 521   fcutsPOI->SetEvent( InputEvent(), MCEvent());\r
 522   if( fflowEvent ){ \r
 523     fflowEvent->~AliFlowEvent();\r
 524     new(fflowEvent) AliFlowEvent(fcutsRP,fcutsPOI);\r
 525   }else fflowEvent = new AliFlowEvent(fcutsRP,fcutsPOI);\r
 526   if(mcEvent && mcEvent->GenEventHeader()) {\r
 527     fflowEvent->SetMCReactionPlaneAngle(mcEvent);\r
 528     //if reaction plane not set from elsewhere randomize it before adding flow\r
 529     //if (!fflowEvent->IsSetMCReactionPlaneAngle()) fflowEvent->SetMCReactionPlaneAngle(gRandom->Uniform(0.0,TMath::TwoPi()));\r
 530     mcReactionPlane = TVector2::Phi_0_2pi(fflowEvent->GetMCReactionPlaneAngle());\r
 531     if(mcReactionPlane > TMath::Pi()) mcReactionPlane = mcReactionPlane - TMath::Pi();\r
 532     //printf("MC reaction plane %f\n",mcReactionPlane);\r
 533   }\r
 534   fflowEvent->SetReferenceMultiplicity( nbtracks );\r
 535   fflowEvent->DefineDeadZone(0,0,0,0);\r
 536   //fflowEvent.TagSubeventsInEta(-0.8,-0.1,0.1,0.8);\r
 537 \r
 538   ////////////////\r
 539   // MC\r
 540   ///////////////\r
 541   if(fUseMCReactionPlane) {\r
 542     eventPlanea = mcReactionPlane;\r
 543     diffsub1sub2a = 0.0;\r
 544   }\r
 545 \r
 546   \r
 547   //////////////////////\r
 548   // Fill Histos\r
 549   //////////////////////\r
 550 \r
 551   fHistEV->Fill(binctt,2.0);\r
 552 \r
 553   // Fill\r
 554   valuensparsea[0] = eventPlanea;  \r
 555   fEventPlane->Fill(&valuensparsea[0]);\r
 556 \r
 557   valuensparsef[0] = diffsub1sub2a;\r
 558   fCosRes->Fill(&valuensparsef[0]);\r
 559   \r
 560   //////////////////////////\r
 561   // Loop over ESD track\r
 562   //////////////////////////\r
 563  \r
 564 \r
 565   for(Int_t k = 0; k < nbtracks; k++){\r
 566     \r
 567     AliESDtrack *track = esd->GetTrack(k);\r
 568     if(!track) continue;\r
 569 \r
 570     Bool_t survived = kTRUE;\r
 571     for(Int_t icut = AliHFEcuts::kStepRecKineITSTPC; icut <= AliHFEcuts::kStepHFEcutsTRD; icut++){\r
 572       if(!fHFECuts->CheckParticleCuts(icut + AliHFEcuts::kNcutStepsMCTrack, (TObject *)track)){\r
 573         survived = kFALSE;\r
 574         break;\r
 575       }\r
 576       //printf("Pass the cut %d\n",icut);\r
 577     }\r
 578     \r
 579     if(!survived) continue;\r
 580 \r
 581     // Apply PID for Data\r
 582     if(!fMCPID) {\r
 583       AliHFEpidObject hfetrack;\r
 584       hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);\r
 585       hfetrack.SetRecTrack(track);\r
 586       hfetrack.SetCentrality((Int_t)binct);\r
 587       //printf("centrality %f and %d\n",binct,hfetrack.GetCentrality());\r
 588       hfetrack.SetPbPb();\r
 589       if(!fPID->IsSelected(&hfetrack,0x0,"",fPIDqa)) {\r
 590         continue;\r
 591       }\r
 592     }\r
 593     else {\r
 594       if(!mcEvent) continue;\r
 595       if(!(mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel()))))) continue;\r
 596       //printf("PdgCode %d\n",TMath::Abs(mctrack->Particle()->GetPdgCode()));\r
 597       if(TMath::Abs(mctrack->Particle()->GetPdgCode())!=11) continue;\r
 598     }\r
 599 \r
 600 \r
 601     /////////////////////////////////////////////////////////////////////////////\r
 602     // Add candidate to AliFlowEvent for POI and subtract from RP if needed\r
 603     ////////////////////////////////////////////////////////////////////////////\r
 604     Int_t idtrack = static_cast<AliVTrack*>(track)->GetID();\r
 605     Bool_t found = kFALSE;\r
 606     Int_t numberoffound = 0;\r
 607     //printf("A: Number of tracks %d\n",fflowEvent->NumberOfTracks());\r
 608     for(Int_t iRPs=0; iRPs< fflowEvent->NumberOfTracks(); iRPs++) {\r
 609       AliFlowTrack *iRP = (AliFlowTrack*) (fflowEvent->GetTrack(iRPs));\r
 610       //if(!iRP->InRPSelection()) continue;\r
 611       if( TMath::Abs(idtrack) == TMath::Abs(iRP->GetID()) ) {\r
 612         iRP->SetForPOISelection(kTRUE);\r
 613         found = kTRUE;\r
 614         numberoffound ++;\r
 615       }\r
 616     }\r
 617     //printf("Found %d mal\n",numberoffound);\r
 618     if(!found) {\r
 619       AliFlowCandidateTrack *sTrack = (AliFlowCandidateTrack*) MakeTrack(massElectron,track->Pt(),track->Phi(), track->Eta());\r
 620       sTrack->SetID(idtrack);\r
 621       fflowEvent->AddTrack(sTrack);\r
 622       //printf("Add the track\n");\r
 623     }\r
 624     //printf("B: Number of tracks %d\n",fflowEvent->NumberOfTracks());\r
 625     \r
 626     \r
 627     /////////////////////////////////////////////////////////\r
 628     // Subtract electron candidate from TPC event plane\r
 629     ////////////////////////////////////////////////////////\r
 630 \r
 631     // Subtract the tracks from the event plane\r
 632     Double_t qX = standardQ->X() - esdEPa->GetQContributionX(track);  //Modify the components: subtract the track you want to look at with your analysis\r
 633     Double_t qY = standardQ->Y() - esdEPa->GetQContributionY(track);  //Modify the components: subtract the track you want to look at with your analysis\r
 634     TVector2 newQVectorfortrack;\r
 635     newQVectorfortrack.Set(qX,qY);\r
 636     Float_t eventplanesubtracted = TVector2::Phi_0_2pi(newQVectorfortrack.Phi())/2; \r
 637 \r
 638     ////////////////////////////////////////\r
 639     // Fill pt and eta for the THnSparseF\r
 640     ///////////////////////////////////////\r
 641 \r
 642     valuensparsee[2] = track->Pt();\r
 643     valuensparsee[3] = track->Eta();    \r
 644     valuensparseg[2] = track->Pt();\r
 645     valuensparseh[2] = track->Pt();\r
 646 \r
 647     Bool_t fillTPC = kTRUE;\r
 648     if((!qsub1a) || (!qsub2a)) fillTPC = kFALSE;\r
 649 \r
 650     if(fSubEtaGapTPC) {\r
 651       if(track->Eta() < (- fEtaGap/2.)) eventplanesubtracted = eventPlanesub1a;\r
 652       else if(track->Eta() > (fEtaGap/2.)) eventplanesubtracted = eventPlanesub2a;\r
 653       else fillTPC = kFALSE;\r
 654     }\r
 655 \r
 656     ///////////////\r
 657     // MC\r
 658     //////////////\r
 659     if(fUseMCReactionPlane) {\r
 660       eventplanesubtracted = mcReactionPlane;\r
 661       fillTPC = kTRUE;\r
 662     }\r
 663     \r
 664     //////////////////////////////////////////////////////////////////////////////\r
 665     ///////////////////////////AFTERBURNER\r
 666     Double_t phitrack = track->Phi();    \r
 667     if (fAfterBurnerOn)\r
 668       {\r
 669         phitrack = GetPhiAfterAddV2(track->Phi(),mcReactionPlane);\r
 670       }\r
 671     //////////////////////////////////////////////////////////////////////////////\r
 672 \r
 673 \r
 674     ///////////////////////\r
 675     // Calculate deltaphi\r
 676     ///////////////////////\r
 677     \r
 678     // Suppose phi track is between 0.0 and phi\r
 679     Double_t deltaphi = TVector2::Phi_0_2pi(phitrack - eventplanesubtracted);\r
 680     if(deltaphi > TMath::Pi()) deltaphi = deltaphi - TMath::Pi();\r
 681    \r
 682     /////////////////////\r
 683     // Fill THnSparseF\r
 684     /////////////////////\r
 685 \r
 686     //\r
 687     valuensparsea[0] = eventplanesubtracted;\r
 688     if(fillTPC) fEventPlaneaftersubtraction->Fill(&valuensparsea[0]);\r
 689     \r
 690     //\r
 691     valuensparsee[0] = TMath::Cos(2*phitrack);\r
 692     fCos2phie->Fill(&valuensparsee[0]);\r
 693     valuensparsee[0] = TMath::Sin(2*phitrack);\r
 694     fSin2phie->Fill(&valuensparsee[0]);\r
 695     \r
 696     //\r
 697     valuensparsee[0] = TMath::Cos(2*eventplanesubtracted);\r
 698     if(fillTPC) fCos2phiep->Fill(&valuensparsee[0]);\r
 699     valuensparsee[0] = TMath::Sin(2*eventplanesubtracted);\r
 700     if(fillTPC) fSin2phiep->Fill(&valuensparsee[0]);\r
 701     valuensparsee[0] = TMath::Sin(2*TVector2::Phi_mpi_pi(phitrack-eventplanesubtracted));\r
 702     if(fillTPC) fSin2phiephiep->Fill(&valuensparsee[0]);\r
 703     \r
 704     // \r
 705     valuensparseg[0] = deltaphi;\r
 706     if(fillTPC) fDeltaPhiMaps->Fill(&valuensparseg[0]);\r
 707     \r
 708     //\r
 709     valuensparseh[0] = TMath::Cos(2*TVector2::Phi_mpi_pi(phitrack-eventplanesubtracted));\r
 710     if(fillTPC) fCosPhiMaps->Fill(&valuensparseh[0]);\r
 711     \r
 712   }\r
 713 \r
 714   //////////////////////////////////////////////////////////////////////////////\r
 715   ///////////////////////////AFTERBURNER\r
 716   if (fAfterBurnerOn)\r
 717     {\r
 718       fflowEvent->AddFlow(fV1,fV2,fV3,fV4,fV5);     //add flow\r
 719       fflowEvent->CloneTracks(fNonFlowNumberOfTrackClones); //add nonflow by cloning tracks\r
 720     }\r
 721   //////////////////////////////////////////////////////////////////////////////\r
 722 \r
 723 \r
 724 \r
 725   if((0.0< cntr) && (cntr<20.0))  PostData(2,fflowEvent);\r
 726   if((20.0< cntr) && (cntr<40.0))  PostData(3,fflowEvent);\r
 727   if((40.0< cntr) && (cntr<60.0))  PostData(4,fflowEvent);\r
 728   if((60.0< cntr) && (cntr<80.0))  PostData(5,fflowEvent);\r
 729   \r
 730  \r
 731   PostData(1, fListHist);\r
 732 \r
 733 \r
 734  \r
 735 }\r
 736 //______________________________________________________________________________\r
 737 AliFlowCandidateTrack *AliAnalysisTaskHFEFlow::MakeTrack( Double_t mass, \r
 738                           Double_t pt, Double_t phi, Double_t eta) {\r
 739   //\r
 740   //  Make Track (Not needed actually)\r
 741   //\r
 742 \r
 743   AliFlowCandidateTrack *sTrack = new AliFlowCandidateTrack();\r
 744   sTrack->SetMass(mass);\r
 745   sTrack->SetPt(pt);\r
 746   sTrack->SetPhi(phi);\r
 747   sTrack->SetEta(eta);\r
 748   sTrack->SetForPOISelection(kTRUE);\r
 749   sTrack->SetForRPSelection(kFALSE);\r
 750   return sTrack;\r
 751 }\r
 752 //_________________________________________________________________________________ \r
 753 Double_t AliAnalysisTaskHFEFlow::GetPhiAfterAddV2(Double_t phi,Double_t reactionPlaneAngle) const\r
 754 {\r
 755   //\r
 756   // Adds v2, uses Newton-Raphson iteration\r
 757   //\r
 758   Double_t phiend=phi;\r
 759   Double_t phi0=phi;\r
 760   Double_t f=0.;\r
 761   Double_t fp=0.;\r
 762   Double_t phiprev=0.;\r
 763 \r
 764   for (Int_t i=0; i<fMaxNumberOfIterations; i++)\r
 765   {\r
 766     phiprev=phiend; //store last value for comparison\r
 767     f =  phiend-phi0+fV2*TMath::Sin(2.*(phiend-reactionPlaneAngle));\r
 768     fp = 1.0+2.0*fV2*TMath::Cos(2.*(phiend-reactionPlaneAngle)); //first derivative\r
 769     phiend -= f/fp;\r
 770     if (TMath::AreEqualAbs(phiprev,phiend,fPrecisionPhi)) break;\r
 771   }\r
 772   return phiend;\r
 773 }\r