PWGHF/hfe/AliAnalysisTaskHFEFlowData.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 "AliVEventHandler.h"\r
  33 #include "AliAnalysisTaskSE.h"\r
  34 #include "AliAnalysisManager.h"\r
  35 \r
  36 #include "AliVEvent.h"\r
  37 #include "AliESDInputHandler.h"\r
  38 #include "AliMCEvent.h"\r
  39 #include "AliESD.h"\r
  40 #include "AliESDEvent.h"\r
  41 #include "AliPIDResponse.h"\r
  42 #include "AliESDVZERO.h"\r
  43 #include "AliESDUtils.h"\r
  44 #include "AliMCParticle.h"\r
  45 #include "AliVTrack.h"\r
  46 #include "AliAODTrack.h"\r
  47 \r
  48 #include "AliFlowCandidateTrack.h"\r
  49 #include "AliFlowEvent.h"\r
  50 #include "AliFlowTrackCuts.h"\r
  51 #include "AliFlowVector.h"\r
  52 #include "AliFlowCommonConstants.h"\r
  53 \r
  54 #include "AliHFEcuts.h"\r
  55 #include "AliHFEpid.h"\r
  56 #include "AliHFEpidQAmanager.h"\r
  57 #include "AliHFEtools.h"\r
  58 #include "AliHFEVZEROEventPlane.h"\r
  59 \r
  60 #include "AliCentrality.h"\r
  61 #include "AliEventplane.h"\r
  62 #include "AliAnalysisTaskHFEFlowData.h"\r
  63 \r
  64 \r
  65 //____________________________________________________________________\r
  66 AliAnalysisTaskHFEFlowData::AliAnalysisTaskHFEFlowData() :\r
  67   AliAnalysisTaskSE(),\r
  68   fListHist(0x0), \r
  69   fAODAnalysis(kFALSE),\r
  70   fUseFlagAOD(kFALSE),\r
  71   fApplyCut(kTRUE),\r
  72   fFlags(1<<4),\r
  73   fVZEROEventPlane(kFALSE),\r
  74   fVZEROEventPlaneA(kFALSE),\r
  75   fVZEROEventPlaneC(kFALSE),\r
  76   fSubEtaGapTPC(kFALSE),\r
  77   fEtaGap(0.0),\r
  78   fDebugLevel(0),\r
  79   fHFECuts(0),\r
  80   fPID(0),\r
  81   fPIDqa(0),\r
  82   fHistEV(0),\r
  83   fEventPlane(0x0),\r
  84   fCosResabc(0x0),\r
  85   fCosRes(0x0),\r
  86   fDeltaPhiMaps(0x0),\r
  87   fCosPhiMaps(0x0)\r
  88 {\r
  89   // Constructor\r
  90 \r
  91   \r
  92 }\r
  93 //______________________________________________________________________________\r
  94 AliAnalysisTaskHFEFlowData:: AliAnalysisTaskHFEFlowData(const char *name) :\r
  95   AliAnalysisTaskSE(name),\r
  96   fListHist(0x0),\r
  97   fAODAnalysis(kFALSE),\r
  98   fUseFlagAOD(kFALSE),\r
  99   fApplyCut(kTRUE),\r
 100   fFlags(1<<4), \r
 101   fVZEROEventPlane(kFALSE),\r
 102   fVZEROEventPlaneA(kFALSE),\r
 103   fVZEROEventPlaneC(kFALSE),\r
 104   fSubEtaGapTPC(kFALSE),\r
 105   fEtaGap(0.0),\r
 106   fDebugLevel(0),\r
 107   fHFECuts(0),\r
 108   fPID(0),\r
 109   fPIDqa(0),\r
 110   fHistEV(0),\r
 111   fEventPlane(0x0),\r
 112   fCosResabc(0x0),\r
 113   fCosRes(0x0),\r
 114   fDeltaPhiMaps(0x0),\r
 115   fCosPhiMaps(0x0)\r
 116 {\r
 117   //\r
 118   // named ctor\r
 119   //\r
 120   \r
 121   fPID = new AliHFEpid("hfePid");\r
 122   fPIDqa = new AliHFEpidQAmanager;\r
 123 \r
 124   DefineInput(0,TChain::Class());\r
 125   DefineOutput(1, TList::Class());\r
 126    \r
 127 }\r
 128 //____________________________________________________________\r
 129 AliAnalysisTaskHFEFlowData::AliAnalysisTaskHFEFlowData(const AliAnalysisTaskHFEFlowData &ref):\r
 130   AliAnalysisTaskSE(ref),\r
 131   fListHist(0x0),\r
 132   fAODAnalysis(ref.fAODAnalysis), \r
 133   fUseFlagAOD(ref.fUseFlagAOD),\r
 134   fApplyCut(ref.fApplyCut),\r
 135   fFlags(ref.fFlags),\r
 136   fVZEROEventPlane(ref.fVZEROEventPlane),\r
 137   fVZEROEventPlaneA(ref.fVZEROEventPlaneA),\r
 138   fVZEROEventPlaneC(ref.fVZEROEventPlaneC),\r
 139   fSubEtaGapTPC(ref.fSubEtaGapTPC),\r
 140   fEtaGap(ref.fEtaGap),\r
 141   fDebugLevel(ref.fDebugLevel),\r
 142   fHFECuts(0),\r
 143   fPID(0),\r
 144   fPIDqa(0),\r
 145   fHistEV(0),\r
 146   fEventPlane(0x0),\r
 147   fCosResabc(0x0),\r
 148   fCosRes(0x0),\r
 149   fDeltaPhiMaps(0x0),\r
 150   fCosPhiMaps(0x0)\r
 151 {\r
 152   //\r
 153   // Copy Constructor\r
 154   //\r
 155   ref.Copy(*this);\r
 156 }\r
 157 \r
 158 //____________________________________________________________\r
 159 AliAnalysisTaskHFEFlowData &AliAnalysisTaskHFEFlowData::operator=(const AliAnalysisTaskHFEFlowData &ref){\r
 160   //\r
 161   // Assignment operator\r
 162   //\r
 163   if(this == &ref) \r
 164     ref.Copy(*this);\r
 165   return *this;\r
 166 }\r
 167 \r
 168 //____________________________________________________________\r
 169 void AliAnalysisTaskHFEFlowData::Copy(TObject &o) const {\r
 170   // \r
 171   // Copy into object o\r
 172   //\r
 173   AliAnalysisTaskHFEFlowData &target = dynamic_cast<AliAnalysisTaskHFEFlowData &>(o);\r
 174   target.fAODAnalysis = fAODAnalysis;\r
 175   target.fUseFlagAOD = fUseFlagAOD;\r
 176   target.fApplyCut = fApplyCut;\r
 177   target.fFlags = fFlags;\r
 178   target.fVZEROEventPlane = fVZEROEventPlane;\r
 179   target.fVZEROEventPlaneA = fVZEROEventPlaneA;\r
 180   target.fVZEROEventPlaneC = fVZEROEventPlaneC;\r
 181   target.fSubEtaGapTPC = fSubEtaGapTPC;\r
 182   target.fEtaGap = fEtaGap;\r
 183   target.fDebugLevel = fDebugLevel;\r
 184   target.fHFECuts = fHFECuts;\r
 185   target.fPID = fPID;\r
 186   target.fPIDqa = fPIDqa;\r
 187   \r
 188 }\r
 189 //____________________________________________________________\r
 190 AliAnalysisTaskHFEFlowData::~AliAnalysisTaskHFEFlowData(){\r
 191   //\r
 192   // Destructor\r
 193   //\r
 194   if(fListHist) delete fListHist;\r
 195   if(fHFECuts) delete fHFECuts;\r
 196   if(fPID) delete fPID;\r
 197   if(fPIDqa) delete fPIDqa;\r
 198  \r
 199 \r
 200 }\r
 201 //________________________________________________________________________\r
 202 void AliAnalysisTaskHFEFlowData::UserCreateOutputObjects()\r
 203 {\r
 204 \r
 205   //********************\r
 206   // Create histograms\r
 207   //********************\r
 208 \r
 209   //**************\r
 210   // Cuts\r
 211   //**************\r
 212 \r
 213   AliDebug(1,"AliAnalysisTaskHFEFlowData: create output objects");\r
 214 \r
 215   // HFE cuts\r
 216 \r
 217   if(!fHFECuts){\r
 218     fHFECuts = new AliHFEcuts;\r
 219     fHFECuts->CreateStandardCuts();\r
 220   }\r
 221   fHFECuts->Initialize();\r
 222   if(fAODAnalysis) fHFECuts->SetAOD();  \r
 223 \r
 224   AliDebug(1,"AliAnalysisTaskHFEFlowData: HFE cuts initialize");\r
 225 \r
 226   // PID HFE\r
 227   //fPID->SetHasMCData(HasMCData());\r
 228   if(!fPID->GetNumberOfPIDdetectors()) fPID->AddDetector("TPC", 0);\r
 229   fPID->InitializePID();\r
 230   fPIDqa->Initialize(fPID);\r
 231   fPID->SortDetectors();\r
 232 \r
 233   AliDebug(1,"AliAnalysisTaskHFEFlowData: pid and pidqa");\r
 234   \r
 235   //**************************\r
 236   // Bins for the THnSparse\r
 237   //**************************\r
 238 \r
 239   Int_t nBinsPt = 24;\r
 240   Double_t binLimPt[25] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2,\r
 241                            1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 3., 3.5, 4., 5.,\r
 242                            6.};\r
 243   \r
 244   Int_t nBinsEtaLess = 2;\r
 245   Double_t minEta = -0.8;\r
 246   Double_t maxEta = 0.8;\r
 247   Double_t binLimEtaLess[nBinsEtaLess+1];\r
 248   for(Int_t i=0; i<=nBinsEtaLess; i++) binLimEtaLess[i]=(Double_t)minEta + (maxEta-minEta)/nBinsEtaLess*(Double_t)i ;\r
 249  \r
 250   Int_t nBinsCos = 50;\r
 251   Double_t minCos = -1.0;\r
 252   Double_t maxCos = 1.0;\r
 253   Double_t binLimCos[nBinsCos+1];\r
 254   for(Int_t i=0; i<=nBinsCos; i++) binLimCos[i]=(Double_t)minCos + (maxCos-minCos)/nBinsCos*(Double_t)i ;\r
 255  \r
 256   Int_t nBinsC = 11;\r
 257   Double_t minC = 0.0;\r
 258   Double_t maxC = 11.0;\r
 259   Double_t binLimC[nBinsC+1];\r
 260   for(Int_t i=0; i<=nBinsC; i++) binLimC[i]=(Double_t)minC + (maxC-minC)/nBinsC*(Double_t)i ;\r
 261 \r
 262   Int_t nBinsCMore = 20;\r
 263   Double_t minCMore = 0.0;\r
 264   Double_t maxCMore = 20.0;\r
 265   Double_t binLimCMore[nBinsCMore+1];\r
 266   for(Int_t i=0; i<=nBinsCMore; i++) binLimCMore[i]=(Double_t)minCMore + (maxCMore-minCMore)/nBinsCMore*(Double_t)i ;\r
 267 \r
 268   Int_t nBinsPhi = 25;\r
 269   Double_t minPhi = 0.0;\r
 270   Double_t maxPhi = TMath::Pi();\r
 271   Double_t binLimPhi[nBinsPhi+1];\r
 272   for(Int_t i=0; i<=nBinsPhi; i++) {\r
 273     binLimPhi[i]=(Double_t)minPhi + (maxPhi-minPhi)/nBinsPhi*(Double_t)i ;\r
 274   }\r
 275 \r
 276   Int_t nBinsCharge = 2;\r
 277   Double_t minCharge = -1.0;\r
 278   Double_t maxCharge = 1.0;\r
 279   Double_t binLimCharge[nBinsCharge+1];\r
 280   for(Int_t i=0; i<=nBinsCharge; i++) binLimCharge[i]=(Double_t)minCharge + (maxCharge-minCharge)/nBinsCharge*(Double_t)i ;\r
 281   \r
 282   AliDebug(1,"AliAnalysisTaskHFEFlowData: array of bins");\r
 283 \r
 284   //******************\r
 285   // Histograms\r
 286   //******************\r
 287     \r
 288   fListHist = new TList();\r
 289   fListHist->SetOwner();\r
 290 \r
 291   AliDebug(1,"AliAnalysisTaskHFEFlowData: list created");\r
 292   \r
 293   // Histos\r
 294   fHistEV = new TH2D("fHistEV", "events", 3, 0, 3, 3, 0,3);\r
 295   \r
 296   AliDebug(1,"AliAnalysisTaskHFEFlowData: fHistEv");\r
 297 \r
 298   // Event plane as function of phiep, centrality\r
 299   const Int_t nDima=5;\r
 300   Int_t nBina[nDima] = {nBinsPhi,nBinsPhi,nBinsPhi,nBinsPhi,nBinsC};\r
 301   fEventPlane = new THnSparseF("EventPlane","EventPlane",nDima,nBina);\r
 302   fEventPlane->SetBinEdges(0,binLimPhi);\r
 303   fEventPlane->SetBinEdges(1,binLimPhi);\r
 304   fEventPlane->SetBinEdges(2,binLimPhi);\r
 305   fEventPlane->SetBinEdges(3,binLimPhi);\r
 306   fEventPlane->SetBinEdges(4,binLimC);\r
 307   fEventPlane->Sumw2();\r
 308 \r
 309   AliDebug(1,"AliAnalysisTaskHFEFlowData: fEventPlane");\r
 310   \r
 311   // Resolution cosres_abc centrality\r
 312   const Int_t nDimfbis=4;\r
 313   Int_t nBinfbis[nDimfbis] = {nBinsCos,nBinsCos,nBinsCos,nBinsCMore};\r
 314   fCosResabc = new THnSparseF("CosRes_abc","CosRes_abc",nDimfbis,nBinfbis);\r
 315   fCosResabc->SetBinEdges(0,binLimCos);\r
 316   fCosResabc->SetBinEdges(1,binLimCos);\r
 317   fCosResabc->SetBinEdges(2,binLimCos);\r
 318   fCosResabc->SetBinEdges(3,binLimCMore);\r
 319   fCosResabc->Sumw2();\r
 320 \r
 321   AliDebug(1,"AliAnalysisTaskHFEFlowData: fCosResabc");\r
 322 \r
 323   // Resolution cosres centrality\r
 324   const Int_t nDimf=2;\r
 325   Int_t nBinf[nDimf] = {nBinsCos, nBinsCMore};\r
 326   fCosRes = new THnSparseF("CosRes","CosRes",nDimf,nBinf);\r
 327   fCosRes->SetBinEdges(0,binLimCos);\r
 328   fCosRes->SetBinEdges(1,binLimCMore);\r
 329   fCosRes->Sumw2();\r
 330 \r
 331   AliDebug(1,"AliAnalysisTaskHFEFlowData: fCosRes");\r
 332   \r
 333   // Maps delta phi\r
 334   const Int_t nDimg=5;\r
 335   Int_t nBing[nDimg] = {nBinsPhi,nBinsC,nBinsPt, nBinsCharge,nBinsEtaLess};\r
 336   fDeltaPhiMaps = new THnSparseF("DeltaPhiMaps","DeltaPhiMaps",nDimg,nBing);\r
 337   fDeltaPhiMaps->SetBinEdges(0,binLimPhi);\r
 338   fDeltaPhiMaps->SetBinEdges(1,binLimC);\r
 339   fDeltaPhiMaps->SetBinEdges(2,binLimPt);\r
 340   fDeltaPhiMaps->SetBinEdges(3,binLimCharge);\r
 341   fDeltaPhiMaps->SetBinEdges(4,binLimEtaLess);\r
 342   fDeltaPhiMaps->Sumw2();  \r
 343 \r
 344   AliDebug(1,"AliAnalysisTaskHFEFlowData: fDeltaPhiMaps");\r
 345 \r
 346   // Maps cos phi\r
 347   const Int_t nDimh=5;\r
 348   Int_t nBinh[nDimh] = {nBinsCos,nBinsC,nBinsPt,nBinsCharge,nBinsEtaLess};\r
 349   fCosPhiMaps = new THnSparseF("CosPhiMaps","CosPhiMaps",nDimh,nBinh);\r
 350   fCosPhiMaps->SetBinEdges(0,binLimCos);\r
 351   fCosPhiMaps->SetBinEdges(1,binLimC);\r
 352   fCosPhiMaps->SetBinEdges(2,binLimPt);\r
 353   fCosPhiMaps->SetBinEdges(3,binLimCharge);\r
 354   fCosPhiMaps->SetBinEdges(4,binLimEtaLess);\r
 355   fCosPhiMaps->Sumw2();\r
 356 \r
 357   AliDebug(1,"AliAnalysisTaskHFEFlowData: fCosPhiMaps");\r
 358 \r
 359   //**************************\r
 360   // Add to the list\r
 361   //******************************\r
 362 \r
 363   //fListHist->Add(qaCutsRP);\r
 364   fListHist->Add(fPIDqa->MakeList("HFEpidQA"));\r
 365   fListHist->Add(fHistEV);\r
 366   fListHist->Add(fEventPlane);\r
 367   fListHist->Add(fCosRes);\r
 368   fListHist->Add(fCosResabc);\r
 369   fListHist->Add(fDeltaPhiMaps);\r
 370   fListHist->Add(fCosPhiMaps);\r
 371 \r
 372   AliDebug(1,"AliAnalysisTaskHFEFlowData: added to the list");\r
 373   \r
 374   \r
 375   PostData(1, fListHist);\r
 376 \r
 377   AliDebug(1,"AliAnalysisTaskHFEFlowData: Post Data");\r
 378 \r
 379 }\r
 380    \r
 381 //________________________________________________________________________\r
 382 void AliAnalysisTaskHFEFlowData::UserExec(Option_t */*option*/)\r
 383 {\r
 384   //\r
 385   // Loop over event\r
 386   //\r
 387 \r
 388   AliDebug(1,"AliAnalysisTaskHFEFlowData: UserExec");\r
 389    \r
 390   Float_t cntr = 0.0;\r
 391   Double_t binct = 11.5;\r
 392   Double_t binctMore = 20.5;\r
 393   Float_t binctt = -1.0;\r
 394   \r
 395   Double_t valuensparsea[5];\r
 396   Double_t valuensparsefbis[4];\r
 397   Double_t valuensparsef[2];\r
 398   Double_t valuensparseg[5];\r
 399   Double_t valuensparseh[5];\r
 400 \r
 401   AliDebug(1, "Variable initialized");\r
 402 \r
 403   \r
 404   /////////////////\r
 405   // centrality\r
 406   /////////////////\r
 407   \r
 408   //AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());\r
 409   //if(!esd) return;\r
 410   AliCentrality *centrality = fInputEvent->GetCentrality();\r
 411   //printf("Got the centrality\n");\r
 412   if(!centrality) {\r
 413     PostData(1, fListHist);\r
 414     return;\r
 415   }\r
 416   cntr = centrality->GetCentralityPercentile("V0M");\r
 417   if((0.0< cntr) && (cntr<5.0)) binct = 0.5;\r
 418   if((5.0< cntr) && (cntr<10.0)) binct = 1.5;\r
 419   if((10.0< cntr) && (cntr<20.0)) binct = 2.5;\r
 420   if((20.0< cntr) && (cntr<30.0)) binct = 3.5;\r
 421   if((30.0< cntr) && (cntr<40.0)) binct = 4.5;\r
 422   if((40.0< cntr) && (cntr<50.0)) binct = 5.5;\r
 423   if((50.0< cntr) && (cntr<60.0)) binct = 6.5;\r
 424   if((60.0< cntr) && (cntr<70.0)) binct = 7.5;\r
 425   if((70.0< cntr) && (cntr<80.0)) binct = 8.5;\r
 426   if((80.0< cntr) && (cntr<90.0)) binct = 9.5;\r
 427   if((90.0< cntr) && (cntr<100.0)) binct = 10.5;\r
 428   \r
 429   if((0.< cntr) && (cntr < 20.)) binctt = 0.5;\r
 430   if((20.< cntr) && (cntr < 40.)) binctt = 1.5;\r
 431   if((40.< cntr) && (cntr < 80.)) binctt = 2.5;\r
 432 \r
 433   if((0.0< cntr) && (cntr<5.0)) binctMore = 0.5;\r
 434   if((5.0< cntr) && (cntr<10.0)) binctMore = 1.5;\r
 435   if((10.0< cntr) && (cntr<15.0)) binctMore = 2.5;\r
 436   if((15.0< cntr) && (cntr<20.0)) binctMore = 3.5;\r
 437   if((20.0< cntr) && (cntr<25.0)) binctMore = 4.5;\r
 438   if((25.0< cntr) && (cntr<30.0)) binctMore = 5.5;\r
 439   if((30.0< cntr) && (cntr<35.0)) binctMore = 6.5;\r
 440   if((35.0< cntr) && (cntr<40.0)) binctMore = 7.5;\r
 441   if((40.0< cntr) && (cntr<45.0)) binctMore = 8.5;\r
 442   if((45.0< cntr) && (cntr<50.0)) binctMore = 9.5;\r
 443   if((50.0< cntr) && (cntr<55.0)) binctMore = 10.5;\r
 444   if((55.0< cntr) && (cntr<60.0)) binctMore = 11.5;\r
 445   if((60.0< cntr) && (cntr<65.0)) binctMore = 12.5;\r
 446   if((65.0< cntr) && (cntr<70.0)) binctMore = 13.5;\r
 447   if((70.0< cntr) && (cntr<75.0)) binctMore = 14.5;\r
 448   if((75.0< cntr) && (cntr<80.0)) binctMore = 15.5;\r
 449   if((80.0< cntr) && (cntr<85.0)) binctMore = 16.5;\r
 450   if((85.0< cntr) && (cntr<90.0)) binctMore = 17.5;\r
 451   if((90.0< cntr) && (cntr<95.0)) binctMore = 18.5;\r
 452   if((95.0< cntr) && (cntr<100.0)) binctMore = 19.5;\r
 453 \r
 454   \r
 455   if(binct > 11.0) {\r
 456     PostData(1, fListHist);\r
 457     return;\r
 458   }\r
 459  \r
 460   AliDebug(1, "Centrality");\r
 461 \r
 462   // centrality\r
 463   valuensparsea[4] = binct;  \r
 464   valuensparsef[1] = binctMore;  \r
 465   valuensparsefbis[3] = binctMore;  \r
 466   valuensparseg[1] = binct;\r
 467   valuensparseh[1] = binct; \r
 468   \r
 469   //////////////////////\r
 470   // run number\r
 471   //////////////////////\r
 472 \r
 473   Int_t runnumber = fInputEvent->GetRunNumber();\r
 474   \r
 475   if(!fPID->IsInitialized()){\r
 476     fPID->InitializePID(runnumber);\r
 477   }\r
 478 \r
 479   AliDebug(1, "Run number");\r
 480 \r
 481   //////////\r
 482   // PID\r
 483   //////////\r
 484  \r
 485   AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();\r
 486   if(!pidResponse){\r
 487     PostData(1, fListHist);\r
 488     return;\r
 489   }\r
 490   fPID->SetPIDResponse(pidResponse);\r
 491 \r
 492   AliDebug(1, "PID");\r
 493 \r
 494   fHistEV->Fill(binctt,0.0);\r
 495  \r
 496   AliDebug(1, "fHistEv");\r
 497 \r
 498   //////////////////\r
 499   // Event cut\r
 500   //////////////////\r
 501   if(!fHFECuts->CheckEventCuts("fEvRecCuts", fInputEvent)) {\r
 502     PostData(1, fListHist);\r
 503     return;\r
 504   }\r
 505 \r
 506   AliDebug(1, "Event cut");\r
 507 \r
 508   fHistEV->Fill(binctt,1.0);\r
 509 \r
 510   ////////////////////////////////////  \r
 511   // First method event plane\r
 512   ////////////////////////////////////\r
 513 \r
 514   AliEventplane* vEPa = fInputEvent->GetEventplane();\r
 515   Float_t eventPlanea = 0.0;\r
 516   Float_t eventPlaneTPC = 0.0;\r
 517   Float_t eventPlaneV0A = 0.0;\r
 518   Float_t eventPlaneV0C = 0.0;\r
 519   Float_t eventPlaneV0 = 0.0;\r
 520   TVector2 *standardQ = 0x0;\r
 521   TVector2 *qsub1a = 0x0;\r
 522   TVector2 *qsub2a = 0x0;\r
 523 \r
 524   // V0\r
 525 \r
 526   eventPlaneV0 = TVector2::Phi_0_2pi(vEPa->GetEventplane("V0", fInputEvent,2));\r
 527   if(eventPlaneV0 > TMath::Pi()) eventPlaneV0 = eventPlaneV0 - TMath::Pi();\r
 528   eventPlaneV0A = TVector2::Phi_0_2pi(vEPa->GetEventplane("V0A", fInputEvent,2));\r
 529   if(eventPlaneV0A > TMath::Pi()) eventPlaneV0A = eventPlaneV0A - TMath::Pi();\r
 530   eventPlaneV0C = TVector2::Phi_0_2pi(vEPa->GetEventplane("V0C", fInputEvent,2));\r
 531   if(eventPlaneV0C > TMath::Pi()) eventPlaneV0C = eventPlaneV0C - TMath::Pi();\r
 532   \r
 533   AliDebug(1, "V0 event plane");\r
 534   \r
 535   // TPC\r
 536 \r
 537   standardQ = vEPa->GetQVector(); \r
 538   Double_t qx = -1.0;\r
 539   Double_t qy = -1.0;\r
 540   if(standardQ) {\r
 541     qx = standardQ->X();\r
 542     qy = standardQ->Y();\r
 543   }  \r
 544   TVector2 qVectorfortrack;\r
 545   qVectorfortrack.Set(qx,qy);\r
 546   eventPlaneTPC = TVector2::Phi_0_2pi(qVectorfortrack.Phi())/2.; \r
 547 \r
 548   AliDebug(1, "TPC event plane");\r
 549 \r
 550   // Choose the one used for v2\r
 551 \r
 552   if(fVZEROEventPlane) eventPlanea = eventPlaneV0;\r
 553   if(fVZEROEventPlaneA) eventPlanea = eventPlaneV0A;\r
 554   if(fVZEROEventPlaneC) eventPlanea = eventPlaneV0C;\r
 555   if(!fVZEROEventPlane) eventPlanea = eventPlaneTPC;\r
 556 \r
 557   Float_t eventPlanesub1a = -100.0;\r
 558   Float_t eventPlanesub2a = -100.0;\r
 559   Double_t diffsub1sub2a = -100.0;\r
 560   Double_t diffsubasubb = -100.0;\r
 561   Double_t diffsubasubc = -100.0;\r
 562   Double_t diffsubbsubc = -100.0;\r
 563   \r
 564   diffsubasubb = TMath::Cos(2.*(eventPlaneV0A - eventPlaneV0C));\r
 565   diffsubasubc = TMath::Cos(2.*(eventPlaneV0A - eventPlaneTPC));\r
 566   diffsubbsubc = TMath::Cos(2.*(eventPlaneV0C - eventPlaneTPC));\r
 567   \r
 568   qsub1a = vEPa->GetQsub1();\r
 569   qsub2a = vEPa->GetQsub2();\r
 570   if(qsub1a) eventPlanesub1a = TVector2::Phi_0_2pi(qsub1a->Phi())/2.;\r
 571   if(qsub2a) eventPlanesub2a = TVector2::Phi_0_2pi(qsub2a->Phi())/2.;\r
 572   if(qsub1a && qsub2a) {\r
 573     diffsub1sub2a = TMath::Cos(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.));\r
 574   }\r
 575 \r
 576   AliDebug(1, "Diff");\r
 577   \r
 578   /////////////////////////////////////////////////////////\r
 579   // Cut for event with event plane reconstructed by all\r
 580   ////////////////////////////////////////////////////////\r
 581   \r
 582   if((!standardQ) || (!qsub1a) || (!qsub2a)) {\r
 583     PostData(1, fListHist);\r
 584     return;\r
 585   }\r
 586 \r
 587   AliDebug(1, "Number of tracks");\r
 588   \r
 589   Int_t nbtracks = fInputEvent->GetNumberOfTracks();\r
 590   \r
 591   //////////////////////\r
 592   // Fill Histos\r
 593   //////////////////////\r
 594 \r
 595   fHistEV->Fill(binctt,2.0);\r
 596   \r
 597   // Fill\r
 598   valuensparsea[0] = eventPlaneV0A;\r
 599   valuensparsea[1] = eventPlaneV0C;\r
 600   valuensparsea[2] = eventPlaneTPC;\r
 601   valuensparsea[3] = eventPlaneV0;  \r
 602   fEventPlane->Fill(&valuensparsea[0]);\r
 603   \r
 604   if(!fVZEROEventPlane) {\r
 605     valuensparsef[0] = diffsub1sub2a;\r
 606     fCosRes->Fill(&valuensparsef[0]);\r
 607   }\r
 608   else {\r
 609     valuensparsefbis[0] = diffsubasubb;\r
 610     valuensparsefbis[1] = diffsubbsubc;\r
 611     valuensparsefbis[2] = diffsubasubc;\r
 612     fCosResabc->Fill(&valuensparsefbis[0]);\r
 613   }\r
 614     \r
 615   \r
 616   //////////////////////////\r
 617   // Loop over ESD track\r
 618   //////////////////////////\r
 619  \r
 620   AliDebug(1, "Loop tracks");\r
 621 \r
 622   for(Int_t k = 0; k < nbtracks; k++){\r
 623     \r
 624     AliVTrack *track = (AliVTrack *) fInputEvent->GetTrack(k);\r
 625     if(!track) continue;\r
 626 \r
 627     if(fAODAnalysis) {\r
 628       AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);\r
 629       if(!aodtrack){\r
 630         AliError("AOD track is not there");\r
 631         PostData(1, fListHist);\r
 632         return;\r
 633       }  \r
 634       //printf("Find AOD track on\n");\r
 635       if(fUseFlagAOD){\r
 636         if(aodtrack->GetFlags() != fFlags) continue;  // Only process AOD tracks where the HFE is set\r
 637       }\r
 638     }\r
 639     \r
 640     if(fApplyCut) {\r
 641       Bool_t survived = kTRUE;\r
 642       for(Int_t icut = AliHFEcuts::kStepRecKineITSTPC; icut <= AliHFEcuts::kStepHFEcutsTRD; icut++){\r
 643         if(!fHFECuts->CheckParticleCuts(icut + AliHFEcuts::kNcutStepsMCTrack, (TObject *)track)){\r
 644           survived = kFALSE;\r
 645           break;\r
 646         }\r
 647       }\r
 648       if(!survived) continue;\r
 649     }\r
 650     \r
 651     // Apply PID\r
 652     AliHFEpidObject hfetrack;\r
 653     if(!fAODAnalysis) hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);\r
 654     else hfetrack.SetAnalysisType(AliHFEpidObject::kAODanalysis);\r
 655     hfetrack.SetRecTrack(track);\r
 656     hfetrack.SetCentrality((Int_t)binct);\r
 657     hfetrack.SetPbPb();\r
 658     if(!fPID->IsSelected(&hfetrack,0x0,"recTrackCont",fPIDqa)) {\r
 659       continue;\r
 660     }\r
 661     \r
 662     /////////////////////////////////////////////////////////\r
 663     // Subtract electron candidate from TPC event plane\r
 664     ////////////////////////////////////////////////////////\r
 665     Float_t eventplanesubtracted = 0.0;    \r
 666 \r
 667     if(!fVZEROEventPlane) {\r
 668       // Subtract the tracks from the event plane\r
 669       Double_t qX = standardQ->X() - vEPa->GetQContributionX(track);  //Modify the components: subtract the track you want to look at with your analysis\r
 670       Double_t qY = standardQ->Y() - vEPa->GetQContributionY(track);  //Modify the components: subtract the track you want to look at with your analysis\r
 671       TVector2 newQVectorfortrack;\r
 672       newQVectorfortrack.Set(qX,qY);\r
 673       eventplanesubtracted = TVector2::Phi_0_2pi(newQVectorfortrack.Phi())/2; \r
 674     }\r
 675     else eventplanesubtracted = eventPlanea;\r
 676 \r
 677     ////////////////////////////////////////\r
 678     // Fill pt and eta for the THnSparseF\r
 679     ///////////////////////////////////////\r
 680 \r
 681     valuensparseg[2] = track->Pt();\r
 682     valuensparseh[2] = track->Pt();\r
 683     if(track->Charge() > 0.0) {\r
 684       valuensparseg[3] = 0.2;\r
 685       valuensparseh[3] = 0.2;\r
 686     }\r
 687     else {\r
 688       valuensparseg[3] = -0.2;\r
 689       valuensparseh[3] = -0.2;\r
 690     }\r
 691     valuensparseh[4] = track->Eta();\r
 692     valuensparseg[4] = track->Eta();\r
 693 \r
 694     ///////////////////////////////\r
 695     // Event plane without track\r
 696     /////////////////////////////\r
 697     Bool_t fillEventPlane = kTRUE;\r
 698     if(!fVZEROEventPlane){\r
 699       if((!qsub1a) || (!qsub2a)) fillEventPlane = kFALSE;\r
 700       if(fSubEtaGapTPC) {\r
 701         if(track->Eta() < (- fEtaGap/2.)) eventplanesubtracted = eventPlanesub1a;\r
 702         else if(track->Eta() > (fEtaGap/2.)) eventplanesubtracted = eventPlanesub2a;\r
 703         else fillEventPlane = kFALSE;\r
 704       }\r
 705     }\r
 706     \r
 707     ///////////////////////\r
 708     // Calculate deltaphi\r
 709     ///////////////////////\r
 710     Double_t phitrack = track->Phi();  \r
 711     Double_t deltaphi = TVector2::Phi_0_2pi(phitrack - eventplanesubtracted);\r
 712     if(deltaphi > TMath::Pi()) deltaphi = deltaphi - TMath::Pi();\r
 713    \r
 714     /////////////////////\r
 715     // Fill THnSparseF\r
 716     /////////////////////\r
 717 \r
 718     valuensparseg[0] = deltaphi;\r
 719     if(fillEventPlane) fDeltaPhiMaps->Fill(&valuensparseg[0]);\r
 720     \r
 721     //\r
 722     valuensparseh[0] = TMath::Cos(2*TVector2::Phi_mpi_pi(phitrack-eventplanesubtracted));\r
 723     if(fillEventPlane) {\r
 724       fCosPhiMaps->Fill(&valuensparseh[0]);\r
 725     }\r
 726     \r
 727   }\r
 728 \r
 729 \r
 730   AliDebug(1, "Post data");\r
 731   \r
 732   PostData(1, fListHist);\r
 733 \r
 734 \r
 735  \r
 736 }\r