PWGHF/hfe/AliAnalysisTaskFlowTPCTOFEPSP.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15 //
  16 // Flow task
  17 //
  18 // Authors:
  19 //   Raphaelle Bailhache <R.Bailhache@gsi.de>
  20 //   Theodor Rascanu <trascanu@stud.uni-frankfurt.de>
  21 //
  22 #include "TROOT.h"
  23 #include "TH1D.h"
  24 #include "TH2D.h"
  25 #include "TChain.h"
  26 #include "TVector2.h"
  27 #include "THnSparse.h"
  28 #include "TMath.h"
  29 #include "TRandom3.h"
  30 #include "TProfile.h"
  31 #include "TProfile2D.h"
  32 #include "TLorentzVector.h"
  33 #include "TParticle.h"
  34 #include "TF1.h"
  35 #include <TArrayD.h>
  36
  37 #include <TDirectory.h>
  38 #include <TTreeStream.h>
  39
  40 #include "AliVEventHandler.h"
  41 #include "AliAnalysisTaskSE.h"
  42 #include "AliAnalysisManager.h"
  43
  44 #include "AliVEvent.h"
  45 #include "AliESDInputHandler.h"
  46 #include "AliMCEvent.h"
  47 #include "AliESD.h"
  48 #include "AliESDEvent.h"
  49 #include "AliPID.h"
  50 #include "AliPIDResponse.h"
  51 #include "AliESDVZERO.h"
  52 #include "AliESDUtils.h"
  53 #include "AliMCParticle.h"
  54 #include "AliAODMCParticle.h"
  55 #include "AliAODEvent.h"
  56 #include "AliAODVertex.h"
  57 #include "AliAODTrack.h"
  58 #include "AliVTrack.h"
  59 #include "AliESDtrack.h"
  60 #include "AliESDtrackCuts.h"
  61 #include "AliAODTrack.h"
  62 #include "AliStack.h"
  63 #include "AliMCEvent.h"
  64
  65 #include "AliFlowCandidateTrack.h"
  66 #include "AliFlowEvent.h"
  67 #include "AliFlowTrackCuts.h"
  68 #include "AliFlowVector.h"
  69 #include "AliFlowCommonConstants.h"
  70 #include "AliKFParticle.h"
  71 #include "AliKFVertex.h"
  72
  73 #include "AliHFEcuts.h"
  74 #include "AliHFEpid.h"
  75 #include "AliHFEpidQAmanager.h"
  76 #include "AliHFEtools.h"
  77 #include "AliHFEVZEROEventPlane.h"
  78
  79 #include "AliCentrality.h"
  80 #include "AliEventplane.h"
  81 #include "AliAnalysisTaskFlowTPCTOFEPSP.h"
  82 #include "AliAODMCHeader.h"
  83 #include "TClonesArray.h"
  84 #include "AliHFENonPhotonicElectron.h"
  85
  86 ClassImp(AliAnalysisTaskFlowTPCTOFEPSP)
  87
  88 //____________________________________________________________________
  89 AliAnalysisTaskFlowTPCTOFEPSP::AliAnalysisTaskFlowTPCTOFEPSP() :
  90 AliAnalysisTaskSE(),
  91   fListHist(0x0),
  92   fAODAnalysis(kFALSE),
  93   fFilter(1<<4),
  94   fAODMCHeader(NULL),
  95   fAODArrayMCInfo(NULL),
  96   fBackgroundSubtraction(NULL),
  97   fVZEROEventPlane(kFALSE),
  98   fVZEROEventPlaneA(kFALSE),
  99   fVZEROEventPlaneC(kFALSE),
 100   fSubEtaGapTPC(kFALSE),
 101   fEtaGap(0.0),
 102   fPtBinning(),
 103   fNbBinsCentralityQCumulant(4),
 104   fNbBinsPtQCumulant(12),
 105   fMinPtQCumulant(0.2),
 106   fMaxPtQCumulant(6.0),
 107   fAfterBurnerOn(kFALSE),
 108   fNonFlowNumberOfTrackClones(0),
 109   fV1(0.),
 110   fV2(0.),
 111   fV3(0.),
 112   fV4(0.),
 113   fV5(0.),
 114   fMaxNumberOfIterations(100),
 115   fPrecisionPhi(0.001),
 116   fUseMCReactionPlane(kFALSE),
 117   fSP(kFALSE),
 118   fMCPID(kFALSE),
 119   fNoPID(kFALSE),
 120   fChi2OverNDFCut(3.0),
 121   fMaxdca(3.0),
 122   fMaxopeningtheta(0.02),
 123   fMaxopeningphi(0.1),
 124   fMaxopening3D(0.1),
 125   fMaxInvmass(0.1),
 126   fSetMassConstraint(kFALSE),
 127   fMonitorEventPlane(kFALSE),
 128   fMonitorContamination(kFALSE),
 129   fMonitorPhotonic(kFALSE),
 130   fMonitorWithoutPID(kFALSE),
 131   fMonitorTrackCuts(kFALSE),
 132   fMonitorQCumulant(kFALSE),
 133   fcutsRP(0),
 134   fcutsPOI(0),
 135   fHFECuts(0),
 136   fPID(0),
 137   fPIDTOFOnly(0),
 138   fPIDqa(0),
 139   fflowEvent(NULL),
 140   fAsFunctionOfP(kTRUE),
 141   fHFEBackgroundCuts(0),
 142   fPIDBackground(0),
 143   fPIDBackgroundqa(0),
 144   fAlgorithmMA(kTRUE),
 145   fArraytrack(NULL),
 146   fCounterPoolBackground(0),
 147   fHFEVZEROEventPlane(0x0),
 148   fHistEV(0),
 149   fHistPileUp(0),
 150   fPileUpCut(kFALSE),
 151   fEventPlane(0x0),
 152   fEventPlaneaftersubtraction(0x0),
 153   fFractionContamination(0x0),
 154   fContaminationv2(0x0),
 155   fContaminationmeanpt(0x0),
 156   fCosSin2phiep(0x0),
 157   fCos2phie(0x0),
 158   fSin2phie(0x0),
 159   fCos2phiep(0x0),
 160   fSin2phiep(0x0),
 161   fSin2phiephiep(0x0),
 162   fCosResabc(0x0),
 163   fSinResabc(0x0),
 164   fProfileCosResab(0x0),
 165   fProfileCosResac(0x0),
 166   fProfileCosResbc(0x0),
 167   fCosRes(0x0),
 168   fSinRes(0x0),
 169   fProfileCosRes(0x0),
 170   fTrackingCuts(0x0),
 171   fDeltaPhiMapsBeforePID(0x0),
 172   fCosPhiMapsBeforePID(0x0),
 173   fDeltaPhiMaps(0x0),
 174   fDeltaPhiMapsContamination(0x0),
 175   fCosPhiMaps(0x0),
 176   fProfileCosPhiMaps(0x0),
 177   fDeltaPhiMapsTaggedPhotonic(0x0),
 178   //fCosPhiMapsTaggedPhotonic(0x0),
 179   fDeltaPhiMapsTaggedNonPhotonic(0x0),
 180   //fCosPhiMapsTaggedNonPhotonic(0x0),
 181   fDeltaPhiMapsTaggedPhotonicLS(0x0),
 182   //fCosPhiMapsTaggedPhotonicLS(0x0),
 183   fMCSourceDeltaPhiMaps(0x0),
 184   fOppSignDeltaPhiMaps(0x0),
 185   fSameSignDeltaPhiMaps(0x0),
 186   fOppSignAngle(0x0),
 187   fSameSignAngle(0x0),
 188   fDebugStreamer(0)
 189 {
 190   // Constructor
 191
 192   for(Int_t k = 0; k < 10; k++) {
 193     fBinCentralityLess[k] = 0.0;
 194   }
 195   for(Int_t k = 0; k < 11; k++) {
 196     fContamination[k] = NULL;
 197     fv2contamination[k] = NULL;
 198   }
 199
 200 }
 201 //______________________________________________________________________________
 202 AliAnalysisTaskFlowTPCTOFEPSP:: AliAnalysisTaskFlowTPCTOFEPSP(const char *name) :
 203   AliAnalysisTaskSE(name),
 204   fListHist(0x0),
 205   fAODAnalysis(kFALSE),
 206   fFilter(1<<4),
 207   fAODMCHeader(NULL),
 208   fAODArrayMCInfo(NULL),
 209   fBackgroundSubtraction(NULL),
 210   fVZEROEventPlane(kFALSE),
 211   fVZEROEventPlaneA(kFALSE),
 212   fVZEROEventPlaneC(kFALSE),
 213   fSubEtaGapTPC(kFALSE),
 214   fEtaGap(0.0),
 215   fPtBinning(),
 216   fNbBinsCentralityQCumulant(4),
 217   fNbBinsPtQCumulant(15),
 218   fMinPtQCumulant(0.0),
 219   fMaxPtQCumulant(6.0),
 220   fAfterBurnerOn(kFALSE),
 221   fNonFlowNumberOfTrackClones(0),
 222   fV1(0.),
 223   fV2(0.),
 224   fV3(0.),
 225   fV4(0.),
 226   fV5(0.),
 227   fMaxNumberOfIterations(100),
 228   fPrecisionPhi(0.001),
 229   fUseMCReactionPlane(kFALSE),
 230   fSP(kFALSE),
 231   fMCPID(kFALSE),
 232   fNoPID(kFALSE),
 233   fChi2OverNDFCut(3.0),
 234   fMaxdca(3.0),
 235   fMaxopeningtheta(0.02),
 236   fMaxopeningphi(0.1),
 237   fMaxopening3D(0.1),
 238   fMaxInvmass(0.1),
 239   fSetMassConstraint(kFALSE),
 240   fMonitorEventPlane(kFALSE),
 241   fMonitorContamination(kFALSE),
 242   fMonitorPhotonic(kFALSE),
 243   fMonitorWithoutPID(kFALSE),
 244   fMonitorTrackCuts(kFALSE),
 245   fMonitorQCumulant(kFALSE),
 246   fcutsRP(0),
 247   fcutsPOI(0),
 248   fHFECuts(0),
 249   fPID(0),
 250   fPIDTOFOnly(0),
 251   fPIDqa(0),
 252   fflowEvent(NULL),
 253   fAsFunctionOfP(kTRUE),
 254   fHFEBackgroundCuts(0),
 255   fPIDBackground(0),
 256   fPIDBackgroundqa(0),
 257   fAlgorithmMA(kTRUE),
 258   fArraytrack(NULL),
 259   fCounterPoolBackground(0),
 260   fHFEVZEROEventPlane(0x0),
 261   fHistEV(0),
 262   fHistPileUp(0),
 263   fPileUpCut(kFALSE),
 264   fEventPlane(0x0),
 265   fEventPlaneaftersubtraction(0x0),
 266   fFractionContamination(0x0),
 267   fContaminationv2(0x0),
 268   fContaminationmeanpt(0x0),
 269   fCosSin2phiep(0x0),
 270   fCos2phie(0x0),
 271   fSin2phie(0x0),
 272   fCos2phiep(0x0),
 273   fSin2phiep(0x0),
 274   fSin2phiephiep(0x0),
 275   fCosResabc(0x0),
 276   fSinResabc(0x0),
 277   fProfileCosResab(0x0),
 278   fProfileCosResac(0x0),
 279   fProfileCosResbc(0x0),
 280   fCosRes(0x0),
 281   fSinRes(0x0),
 282   fProfileCosRes(0x0),
 283   fTrackingCuts(0x0),
 284   fDeltaPhiMapsBeforePID(0x0),
 285   fCosPhiMapsBeforePID(0x0),
 286   fDeltaPhiMaps(0x0),
 287   fDeltaPhiMapsContamination(0x0),
 288   fCosPhiMaps(0x0),
 289   fProfileCosPhiMaps(0x0),
 290   fDeltaPhiMapsTaggedPhotonic(0x0),
 291   //fCosPhiMapsTaggedPhotonic(0x0),
 292   fDeltaPhiMapsTaggedNonPhotonic(0x0),
 293   //fCosPhiMapsTaggedNonPhotonic(0x0),
 294   fDeltaPhiMapsTaggedPhotonicLS(0x0),
 295   //fCosPhiMapsTaggedPhotonicLS(0x0),
 296   fMCSourceDeltaPhiMaps(0x0),
 297   fOppSignDeltaPhiMaps(0x0),
 298   fSameSignDeltaPhiMaps(0x0),
 299   fOppSignAngle(0x0),
 300   fSameSignAngle(0x0),
 301   fDebugStreamer(0)
 302 {
 303   //
 304   // named ctor
 305   //
 306
 307   for(Int_t k = 0; k < 10; k++) {
 308     fBinCentralityLess[k] = 0.0;
 309   }
 310   fBinCentralityLess[0] = 0.0;
 311   fBinCentralityLess[1] = 20.0;
 312   fBinCentralityLess[2] = 40.0;
 313   fBinCentralityLess[3] = 60.0;
 314   fBinCentralityLess[4] = 80.0;
 315
 316   for(Int_t k = 0; k < 11; k++) {
 317     fContamination[k] = NULL;
 318     fv2contamination[k] = NULL;
 319   }
 320
 321   fPID = new AliHFEpid("hfePid");
 322   fPIDqa = new AliHFEpidQAmanager;
 323
 324   fPIDBackground = new AliHFEpid("hfePidBackground");
 325   fPIDBackgroundqa = new AliHFEpidQAmanager;
 326
 327   fPIDTOFOnly = new AliHFEpid("hfePidTOFOnly");
 328
 329   DefineInput(0,TChain::Class());
 330   DefineOutput(1, TList::Class());
 331   //for(Int_t bincless = 0; bincless < fNbBinsCentralityQCumulant; bincless++) {
 332   //  DefineOutput(bincless+2,AliFlowEventSimple::Class());
 333   //}
 334
 335 }
 336 //____________________________________________________________
 337 AliAnalysisTaskFlowTPCTOFEPSP::AliAnalysisTaskFlowTPCTOFEPSP(const AliAnalysisTaskFlowTPCTOFEPSP &ref):
 338   AliAnalysisTaskSE(ref),
 339   fListHist(NULL),
 340   fAODAnalysis(ref.fAODAnalysis),
 341   fFilter(ref.fFilter),
 342   fAODMCHeader(ref.fAODMCHeader),
 343   fAODArrayMCInfo(ref.fAODArrayMCInfo),
 344   fBackgroundSubtraction(ref.fBackgroundSubtraction),
 345   fVZEROEventPlane(ref.fVZEROEventPlane),
 346   fVZEROEventPlaneA(ref.fVZEROEventPlaneA),
 347   fVZEROEventPlaneC(ref.fVZEROEventPlaneC),
 348   fSubEtaGapTPC(ref.fSubEtaGapTPC),
 349   fEtaGap(ref.fEtaGap),
 350   fPtBinning(ref.fPtBinning),
 351   fNbBinsCentralityQCumulant(ref.fNbBinsCentralityQCumulant),
 352   fNbBinsPtQCumulant(ref.fNbBinsPtQCumulant),
 353   fMinPtQCumulant(ref.fMinPtQCumulant),
 354   fMaxPtQCumulant(ref.fMaxPtQCumulant),
 355   fAfterBurnerOn(ref.fAfterBurnerOn),
 356   fNonFlowNumberOfTrackClones(ref.fNonFlowNumberOfTrackClones),
 357   fV1(ref.fV1),
 358   fV2(ref.fV2),
 359   fV3(ref.fV3),
 360   fV4(ref.fV4),
 361   fV5(ref.fV5),
 362   fMaxNumberOfIterations(ref.fMaxNumberOfIterations),
 363   fPrecisionPhi(ref.fPrecisionPhi),
 364   fUseMCReactionPlane(ref.fUseMCReactionPlane),
 365   fSP(ref.fSP),
 366   fMCPID(ref.fMCPID),
 367   fNoPID(ref.fNoPID),
 368   fChi2OverNDFCut(ref.fChi2OverNDFCut),
 369   fMaxdca(ref.fMaxdca),
 370   fMaxopeningtheta(ref.fMaxopeningtheta),
 371   fMaxopeningphi(ref.fMaxopeningphi),
 372   fMaxopening3D(ref.fMaxopening3D),
 373   fMaxInvmass(ref.fMaxInvmass),
 374   fSetMassConstraint(ref.fSetMassConstraint),
 375   fMonitorEventPlane(ref.fMonitorEventPlane),
 376   fMonitorContamination(ref.fMonitorContamination),
 377   fMonitorPhotonic(ref.fMonitorPhotonic),
 378   fMonitorWithoutPID(ref.fMonitorWithoutPID),
 379   fMonitorTrackCuts(ref.fMonitorTrackCuts),
 380   fMonitorQCumulant(ref.fMonitorQCumulant),
 381   fcutsRP(NULL),
 382   fcutsPOI(NULL),
 383   fHFECuts(NULL),
 384   fPID(NULL),
 385   fPIDTOFOnly(NULL),
 386   fPIDqa(NULL),
 387   fflowEvent(NULL),
 388   fAsFunctionOfP(ref.fAsFunctionOfP),
 389   fHFEBackgroundCuts(NULL),
 390   fPIDBackground(NULL),
 391   fPIDBackgroundqa(NULL),
 392   fAlgorithmMA(ref.fAlgorithmMA),
 393   fArraytrack(NULL),
 394   fCounterPoolBackground(ref.fCounterPoolBackground),
 395   fHFEVZEROEventPlane(NULL),
 396   fHistEV(NULL),
 397   fHistPileUp(NULL),
 398   fPileUpCut(kFALSE),
 399   fEventPlane(NULL),
 400   fEventPlaneaftersubtraction(NULL),
 401   fFractionContamination(NULL),
 402   fContaminationv2(NULL),
 403   fContaminationmeanpt(0x0),
 404   fCosSin2phiep(NULL),
 405   fCos2phie(NULL),
 406   fSin2phie(NULL),
 407   fCos2phiep(NULL),
 408   fSin2phiep(NULL),
 409   fSin2phiephiep(NULL),
 410   fCosResabc(NULL),
 411   fSinResabc(NULL),
 412   fProfileCosResab(NULL),
 413   fProfileCosResac(NULL),
 414   fProfileCosResbc(NULL),
 415   fCosRes(NULL),
 416   fSinRes(NULL),
 417   fProfileCosRes(NULL),
 418   fTrackingCuts(NULL),
 419   fDeltaPhiMapsBeforePID(NULL),
 420   fCosPhiMapsBeforePID(NULL),
 421   fDeltaPhiMaps(NULL),
 422   fDeltaPhiMapsContamination(NULL),
 423   fCosPhiMaps(NULL),
 424   fProfileCosPhiMaps(NULL),
 425   fDeltaPhiMapsTaggedPhotonic(NULL),
 426   //fCosPhiMapsTaggedPhotonic(NULL),
 427   fDeltaPhiMapsTaggedNonPhotonic(NULL),
 428   //fCosPhiMapsTaggedNonPhotonic(NULL),
 429   fDeltaPhiMapsTaggedPhotonicLS(NULL),
 430   //fCosPhiMapsTaggedPhotonicLS(NULL),
 431   fMCSourceDeltaPhiMaps(NULL),
 432   fOppSignDeltaPhiMaps(NULL),
 433   fSameSignDeltaPhiMaps(NULL),
 434   fOppSignAngle(NULL),
 435   fSameSignAngle(NULL),
 436   fDebugStreamer(0)
 437 {
 438   //
 439   // Copy Constructor
 440   //
 441
 442   for(Int_t k = 0; k < 10; k++) {
 443     fBinCentralityLess[k] = 0.0;
 444   }
 445   for(Int_t k = 0; k < 11; k++) {
 446     fContamination[k] = NULL;
 447     fv2contamination[k] = NULL;
 448   }
 449
 450
 451   ref.Copy(*this);
 452 }
 453
 454 //____________________________________________________________
 455 AliAnalysisTaskFlowTPCTOFEPSP &AliAnalysisTaskFlowTPCTOFEPSP::operator=(const AliAnalysisTaskFlowTPCTOFEPSP &ref){
 456   //
 457   // Assignment operator
 458   //
 459   if(this == &ref)
 460     ref.Copy(*this);
 461   return *this;
 462 }
 463
 464 //____________________________________________________________
 465 void AliAnalysisTaskFlowTPCTOFEPSP::Copy(TObject &o) const {
 466   //
 467   // Copy into object o
 468   //
 469   AliAnalysisTaskFlowTPCTOFEPSP &target = dynamic_cast<AliAnalysisTaskFlowTPCTOFEPSP &>(o);
 470   target.fListHist = fListHist;
 471   target.fAODAnalysis = fAODAnalysis;
 472   target.fFilter = fFilter;
 473   target.fAODMCHeader = fAODMCHeader;
 474   target.fAODArrayMCInfo = fAODArrayMCInfo;
 475   target.fBackgroundSubtraction = fBackgroundSubtraction;
 476   target.fVZEROEventPlane = fVZEROEventPlane;
 477   target.fVZEROEventPlaneA = fVZEROEventPlaneA;
 478   target.fVZEROEventPlaneC = fVZEROEventPlaneC;
 479   target.fSubEtaGapTPC = fSubEtaGapTPC;
 480   target.fEtaGap = fEtaGap;
 481   target.fPtBinning = fPtBinning;
 482   target.fNbBinsCentralityQCumulant = fNbBinsCentralityQCumulant;
 483   target.fNbBinsPtQCumulant = fNbBinsPtQCumulant;
 484   target.fMinPtQCumulant = fMinPtQCumulant;
 485   target.fMaxPtQCumulant = fMaxPtQCumulant;
 486   target.fAfterBurnerOn = fAfterBurnerOn;
 487   target.fNonFlowNumberOfTrackClones = fNonFlowNumberOfTrackClones;
 488   target.fV1 = fV1;
 489   target.fV2 = fV2;
 490   target.fV3 = fV3;
 491   target.fV4 = fV4;
 492   target.fV5 = fV5;
 493   target.fMaxNumberOfIterations = fMaxNumberOfIterations;
 494   target.fPrecisionPhi = fPrecisionPhi;
 495   target.fUseMCReactionPlane = fUseMCReactionPlane;
 496   target.fSP = fSP;
 497   target.fMCPID = fMCPID;
 498   target.fNoPID = fNoPID;
 499   target.fChi2OverNDFCut = fChi2OverNDFCut;
 500   target.fMaxdca = fMaxdca;
 501   target.fMaxopeningtheta = fMaxopeningtheta;
 502   target.fMaxopeningphi = fMaxopeningphi;
 503   target.fMaxopening3D = fMaxopening3D;
 504   target.fMaxInvmass = fMaxInvmass;
 505   target.fSetMassConstraint =  fSetMassConstraint;
 506   target.fMonitorEventPlane = fMonitorEventPlane;
 507   target.fMonitorContamination = fMonitorContamination;
 508   target.fMonitorPhotonic = fMonitorPhotonic;
 509   target.fMonitorWithoutPID = fMonitorWithoutPID;
 510   target.fMonitorTrackCuts = fMonitorTrackCuts;
 511   target.fMonitorQCumulant = fMonitorQCumulant;
 512   target.fcutsRP = fcutsRP;
 513   target.fcutsPOI = fcutsPOI;
 514   target.fHFECuts = fHFECuts;
 515   target.fPID = fPID;
 516   target.fPIDTOFOnly = fPIDTOFOnly;
 517   target.fPIDqa = fPIDqa;
 518   target.fflowEvent = fflowEvent;
 519   target.fAsFunctionOfP = fAsFunctionOfP;
 520   target.fHFEBackgroundCuts = fHFEBackgroundCuts;
 521   target.fPIDBackground = fPIDBackground;
 522   target.fPIDBackgroundqa = fPIDBackgroundqa;
 523   target.fAlgorithmMA = fAlgorithmMA;
 524   target.fArraytrack = fArraytrack;
 525   target.fCounterPoolBackground = fCounterPoolBackground;
 526   target.fHFEVZEROEventPlane = fHFEVZEROEventPlane;
 527   target.fHistEV=fHistEV;
 528   target.fHistPileUp=fHistPileUp;
 529   target.fPileUpCut=fPileUpCut;
 530   target.fEventPlane=fEventPlane;
 531   target.fEventPlaneaftersubtraction=fEventPlaneaftersubtraction;
 532   target.fFractionContamination=fFractionContamination;
 533   target.fContaminationv2=fContaminationv2;
 534   target.fContaminationmeanpt=fContaminationmeanpt;
 535   target.fCosSin2phiep=fCosSin2phiep;
 536   target.fCos2phie=fCos2phie;
 537   target.fSin2phie=fSin2phie;
 538   target.fCos2phiep=fCos2phiep;
 539   target.fSin2phiep=fSin2phiep;
 540   target.fSin2phiephiep=fSin2phiephiep;
 541   target.fCosResabc=fCosResabc;
 542   target.fSinResabc=fSinResabc;
 543   target.fProfileCosResab=fProfileCosResab;
 544   target.fProfileCosResac=fProfileCosResac;
 545   target.fProfileCosResbc=fProfileCosResbc;
 546   target.fCosRes=fCosRes;
 547   target.fSinRes=fSinRes;
 548   target.fProfileCosRes=fProfileCosRes;
 549   target.fTrackingCuts=fTrackingCuts;
 550   target.fDeltaPhiMapsBeforePID=fDeltaPhiMapsBeforePID;
 551   target.fCosPhiMapsBeforePID=fCosPhiMapsBeforePID;
 552   target.fDeltaPhiMaps=fDeltaPhiMaps;
 553   target.fDeltaPhiMapsContamination=fDeltaPhiMapsContamination;
 554   target.fCosPhiMaps=fCosPhiMaps;
 555   target.fProfileCosPhiMaps=fProfileCosPhiMaps;
 556   target.fDeltaPhiMapsTaggedPhotonic=fDeltaPhiMapsTaggedPhotonic;
 557   target.fDeltaPhiMapsTaggedNonPhotonic=fDeltaPhiMapsTaggedNonPhotonic;
 558   target.fDeltaPhiMapsTaggedPhotonicLS=fDeltaPhiMapsTaggedPhotonicLS;
 559   target.fMCSourceDeltaPhiMaps=fMCSourceDeltaPhiMaps;
 560   target.fOppSignDeltaPhiMaps=fOppSignDeltaPhiMaps;
 561   target.fSameSignDeltaPhiMaps=fSameSignDeltaPhiMaps;
 562   target.fOppSignAngle=fOppSignAngle;
 563   target.fSameSignAngle=fSameSignAngle;
 564
 565
 566   for(Int_t k = 0; k < 10; k++) {
 567     target.fBinCentralityLess[k] = fBinCentralityLess[k];
 568   }
 569   for(Int_t k = 0; k < 11; k++) {
 570     target.fContamination[k] = fContamination[k];
 571     target.fv2contamination[k] = fv2contamination[k];
 572   }
 573   target.fDebugStreamer=fDebugStreamer;
 574 }
 575 //____________________________________________________________
 576 AliAnalysisTaskFlowTPCTOFEPSP::~AliAnalysisTaskFlowTPCTOFEPSP(){
 577   //
 578   // Destructor
 579   //
 580
 581
 582   if(fArraytrack) delete fArraytrack;
 583   if(fListHist) delete fListHist;
 584   if(fcutsRP) delete fcutsRP;
 585   if(fcutsPOI) delete fcutsPOI;
 586   if(fHFECuts) delete fHFECuts;
 587   if(fPID) delete fPID;
 588   if(fPIDTOFOnly) delete fPIDTOFOnly;
 589   //if(fPIDqa) delete fPIDqa;
 590   if(fflowEvent) delete fflowEvent;
 591   if(fHFEBackgroundCuts) delete fHFEBackgroundCuts;
 592   if(fPIDBackground) delete fPIDBackground;
 593   //if(fBackgroundSubtraction) delete fBackgroundSubtraction;
 594   //if(fPIDBackgroundqa) delete fPIDBackgroundqa;
 595   //if(fHFEVZEROEventPlane) delete fHFEVZEROEventPlane;
 596   if ( fDebugStreamer ) delete fDebugStreamer;
 597
 598
 599 }
 600 //________________________________________________________________________
 601 void AliAnalysisTaskFlowTPCTOFEPSP::UserCreateOutputObjects()
 602 {
 603
 604   //********************
 605   // Create histograms
 606   //********************
 607
 608   //**************
 609   // Cuts
 610   //**************
 611
 612   //---------Data selection----------
 613   //kMC, kGlobal, kTPCstandalone, kSPDtracklet, kPMD
 614   //AliFlowTrackCuts::trackParameterType rptype = AliFlowTrackCuts::kGlobal;
 615   //AliFlowTrackCuts::trackParameterType poitype = AliFlowTrackCuts::kGlobal;
 616
 617   //---------Parameter mixing--------
 618   //kPure - no mixing, kTrackWithMCkine, kTrackWithMCPID, kTrackWithMCpt
 619   //AliFlowTrackCuts::trackParameterMix rpmix = AliFlowTrackCuts::kPure;
 620   //AliFlowTrackCuts::trackParameterMix poimix = AliFlowTrackCuts::kPure;
 621
 622   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: User create output objects");
 623
 624
 625   // AOD or ESD
 626   AliVEventHandler *inputHandler = dynamic_cast<AliVEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 627   if(!TString(inputHandler->IsA()->GetName()).CompareTo("AliAODInputHandler")){
 628     SetAODAnalysis(kTRUE);
 629     AliDebug(2,"Put AOD analysis on");
 630   } else {
 631     SetAODAnalysis(kFALSE);
 632   }
 633
 634   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: AOD ESD");
 635
 636   // RP TRACK CUTS:
 637   fcutsRP = AliFlowTrackCuts::GetStandardTPCStandaloneTrackCuts2010();
 638   fcutsRP->SetName("StandartTPC");
 639   fcutsRP->SetEtaRange(-0.9,0.9);
 640   fcutsRP->SetQA(kTRUE);
 641   //TList *qaCutsRP = fcutsRP->GetQA();
 642   //qaCutsRP->SetName("QA_StandartTPC_RP");
 643
 644   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: cutsRP");
 645
 646   //POI TRACK CUTS:
 647   fcutsPOI = new AliFlowTrackCuts("dummy");
 648   fcutsPOI->SetParamType(AliFlowTrackCuts::kGlobal);
 649   fcutsPOI->SetPtRange(+1,-1); // select nothing QUICK
 650   fcutsPOI->SetEtaRange(+1,-1); // select nothing VZERO
 651
 652   if( fflowEvent ){
 653     fflowEvent->~AliFlowEvent();
 654     new(fflowEvent) AliFlowEvent(fcutsRP,fcutsPOI);
 655   }
 656   else fflowEvent = new AliFlowEvent(fcutsRP,fcutsPOI);
 657
 658   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: cutsPOI");
 659
 660   // Flow
 661   AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
 662   cc->SetNbinsMult(10000);
 663   cc->SetMultMin(0);
 664   cc->SetMultMax(10000.);
 665   cc->SetNbinsPt(fNbBinsPtQCumulant);
 666   cc->SetPtMin(fMinPtQCumulant);
 667   cc->SetPtMax(fMaxPtQCumulant);
 668   cc->SetNbinsPhi(180);
 669   cc->SetPhiMin(0.0);
 670   cc->SetPhiMax(TMath::TwoPi());
 671   cc->SetNbinsEta(200);
 672   cc->SetEtaMin(-0.9);
 673   cc->SetEtaMax(+0.9);
 674   cc->SetNbinsQ(500);
 675   cc->SetQMin(0.0);
 676   cc->SetQMax(3.0);
 677
 678   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: common constants");
 679
 680   // HFE cuts
 681
 682   if(!fHFECuts){
 683     fHFECuts = new AliHFEcuts;
 684     fHFECuts->CreateStandardCuts();
 685   }
 686   fHFECuts->Initialize();
 687   if(fAODAnalysis) fHFECuts->SetAOD();
 688
 689   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: HFE cuts");
 690
 691
 692   // PID HFE
 693   //fPID->SetHasMCData(HasMCData());
 694   if(!fPID) {
 695     fPID =new AliHFEpid("hfePid");
 696     AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: pid init 0");
 697   }
 698   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: pid init 1");
 699   if(!fPID->GetNumberOfPIDdetectors()) fPID->AddDetector("TPC", 0);
 700   AliDebug(2,Form("AliAnalysisTaskFlowTPCTOFEPSP: GetNumber of PID detectors %d",fPID->GetNumberOfPIDdetectors()));
 701   fPID->InitializePID();
 702   fPIDqa->Initialize(fPID);
 703   fPID->SortDetectors();
 704
 705   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: pid and pidqa");
 706
 707   if(!fPIDTOFOnly->GetNumberOfPIDdetectors()) {
 708     fPIDTOFOnly->AddDetector("TOF", 0);
 709     fPIDTOFOnly->ConfigureTOF(3.);
 710   }
 711   fPIDTOFOnly->InitializePID();
 712   fPIDTOFOnly->SortDetectors();
 713
 714   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: pidtof");
 715
 716   // HFE Background cuts
 717
 718   if(!fHFEBackgroundCuts){
 719      fHFEBackgroundCuts = new AliESDtrackCuts();
 720      fHFEBackgroundCuts->SetName("nackgroundcuts");
 721      //Configure Default Track Cuts
 722      fHFEBackgroundCuts->SetAcceptKinkDaughters(kFALSE);
 723      fHFEBackgroundCuts->SetRequireTPCRefit(kTRUE);
 724      fHFEBackgroundCuts->SetEtaRange(-0.9,0.9);
 725      fHFEBackgroundCuts->SetRequireSigmaToVertex(kTRUE);
 726      fHFEBackgroundCuts->SetMaxChi2PerClusterTPC(4.0);
 727      fHFEBackgroundCuts->SetMinNClustersTPC(50);
 728      fHFEBackgroundCuts->SetPtRange(0.3,1e10);
 729   }
 730
 731   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: hfe background");
 732
 733   // PID background HFE
 734   if(!fPIDBackground->GetNumberOfPIDdetectors()) fPIDBackground->AddDetector("TPC", 0);
 735   fPIDBackground->InitializePID();
 736   fPIDBackgroundqa->Initialize(fPIDBackground);
 737   fPIDBackground->SortDetectors();
 738
 739   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: pid background");
 740
 741   if (fMonitorPhotonic) {
 742     if(!fBackgroundSubtraction) fBackgroundSubtraction = new AliHFENonPhotonicElectron();
 743     if(fAODAnalysis) fBackgroundSubtraction->SetAOD(kTRUE);
 744     fBackgroundSubtraction->Init();
 745   }
 746
 747
 748
 749   //**************************
 750   // Bins for the THnSparse
 751   //**************************
 752
 753   /*
 754   Int_t nBinsPt = 44;
 755   Double_t minPt = 0.1;
 756   Double_t maxPt = 20.0;
 757   Double_t binLimLogPt[nBinsPt+1];
 758   Double_t binLimPt[nBinsPt+1];
 759   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 ;
 760   for(Int_t i=0; i<=nBinsPt; i++) binLimPt[i]=(Double_t)TMath::Power(10,binLimLogPt[i]);
 761   */
 762
 763   Int_t nBinsPt = 20;
 764   Double_t binLimPt[21] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2,
 765                            1.3, 1.4, 1.5, 2., 2.5, 3., 4., 6.};
 766
 767   if(!fPtBinning.GetSize()) fPtBinning.Set(nBinsPt+1, binLimPt);
 768
 769   Int_t nBinsPtPlus = fNbBinsPtQCumulant;
 770   Double_t minPtPlus = fMinPtQCumulant;
 771   Double_t maxPtPlus = fMaxPtQCumulant;
 772   Double_t binLimPtPlus[nBinsPtPlus+1];
 773   for(Int_t i=0; i<=nBinsPtPlus; i++) binLimPtPlus[i]=(Double_t)minPtPlus + (maxPtPlus-minPtPlus)/nBinsPtPlus*(Double_t)i ;
 774
 775   Int_t nBinsEta = 8;
 776   Double_t minEta = -0.8;
 777   Double_t maxEta = 0.8;
 778   Double_t binLimEta[nBinsEta+1];
 779   for(Int_t i=0; i<=nBinsEta; i++) binLimEta[i]=(Double_t)minEta + (maxEta-minEta)/nBinsEta*(Double_t)i ;
 780
 781   Int_t nBinsStep = 6;
 782   Double_t minStep = 0.;
 783   Double_t maxStep = 6.;
 784   Double_t binLimStep[nBinsStep+1];
 785   for(Int_t i=0; i<=nBinsStep; i++) binLimStep[i]=(Double_t)minStep + (maxStep-minStep)/nBinsStep*(Double_t)i ;
 786
 787   Int_t nBinsEtaLess = 2;
 788   Double_t binLimEtaLess[nBinsEtaLess+1];
 789   for(Int_t i=0; i<=nBinsEtaLess; i++) binLimEtaLess[i]=(Double_t)minEta + (maxEta-minEta)/nBinsEtaLess*(Double_t)i ;
 790
 791   Int_t nBinsCos = 50;
 792   Double_t minCos = -1.0;
 793   Double_t maxCos = 1.0;
 794   Double_t binLimCos[nBinsCos+1];
 795   for(Int_t i=0; i<=nBinsCos; i++) binLimCos[i]=(Double_t)minCos + (maxCos-minCos)/nBinsCos*(Double_t)i ;
 796
 797   // Int_t nBinsCosSP = 50;
 798   // Double_t minCosSP = -100.0;
 799   // Double_t maxCosSP = 100.0;
 800   // Double_t binLimCosSP[nBinsCosSP+1];
 801   // for(Int_t i=0; i<=nBinsCosSP; i++) binLimCosSP[i]=(Double_t)minCosSP + (maxCosSP-minCosSP)/nBinsCosSP*(Double_t)i ;
 802
 803   Int_t nBinsC = 11;
 804   Double_t minC = 0.0;
 805   Double_t maxC = 11.0;
 806   Double_t binLimC[nBinsC+1];
 807   for(Int_t i=0; i<=nBinsC; i++) binLimC[i]=(Double_t)minC + (maxC-minC)/nBinsC*(Double_t)i ;
 808
 809   Int_t nBinsCMore = 20;
 810   Double_t minCMore = 0.0;
 811   Double_t maxCMore = 20.0;
 812   Double_t binLimCMore[nBinsCMore+1];
 813   for(Int_t i=0; i<=nBinsCMore; i++) binLimCMore[i]=(Double_t)minCMore + (maxCMore-minCMore)/nBinsCMore*(Double_t)i ;
 814
 815   Int_t nBinsCEvenMore = 100;
 816   Double_t minCEvenMore = 0.0;
 817   Double_t maxCEvenMore = 100.0;
 818   Double_t binLimCEvenMore[nBinsCEvenMore+1];
 819   for(Int_t i=0; i<=nBinsCEvenMore; i++) binLimCEvenMore[i]=(Double_t)minCEvenMore + (maxCEvenMore-minCEvenMore)/nBinsCEvenMore*(Double_t)i ;
 820
 821   Int_t nBinsPhi = 8;
 822   Double_t minPhi = 0.0;
 823   Double_t maxPhi = TMath::Pi();
 824   Double_t binLimPhi[nBinsPhi+1];
 825   for(Int_t i=0; i<=nBinsPhi; i++) {
 826     binLimPhi[i]=(Double_t)minPhi + (maxPhi-minPhi)/nBinsPhi*(Double_t)i ;
 827     AliDebug(2,Form("bin phi is %f for %d",binLimPhi[i],i));
 828   }
 829
 830   Int_t nBinsPhiLess = 2.0;
 831   Double_t minPhiLess = 0.0;
 832   Double_t maxPhiLess = 2.0;
 833   Double_t binLimPhiLess[nBinsPhiLess+1];
 834   for(Int_t i=0; i<=nBinsPhiLess; i++) {
 835     binLimPhiLess[i]=(Double_t)minPhiLess + (maxPhiLess-minPhiLess)/nBinsPhiLess*(Double_t)i ;
 836   }
 837
 838   Int_t nBinsTPCdEdx = 140;
 839   Double_t minTPCdEdx = -12.0;
 840   Double_t maxTPCdEdx = 12.0;
 841   Double_t binLimTPCdEdx[nBinsTPCdEdx+1];
 842   for(Int_t i=0; i<=nBinsTPCdEdx; i++) {
 843     binLimTPCdEdx[i]=(Double_t)minTPCdEdx + (maxTPCdEdx-minTPCdEdx)/nBinsTPCdEdx*(Double_t)i ;
 844   }
 845
 846   Int_t nBinsAngle = 40;
 847   Double_t minAngle = 0.0;
 848   Double_t maxAngle = 1.0;
 849   Double_t binLimAngle[nBinsAngle+1];
 850   for(Int_t i=0; i<=nBinsAngle; i++) {
 851     binLimAngle[i]=(Double_t)minAngle + (maxAngle-minAngle)/nBinsAngle*(Double_t)i ;
 852     AliDebug(2,Form("bin phi is %f for %d",binLimPhi[i],i));
 853   }
 854
 855   Int_t nBinsCharge = 2;
 856   Double_t minCharge = -1.0;
 857   Double_t maxCharge = 1.0;
 858   Double_t binLimCharge[nBinsCharge+1];
 859   for(Int_t i=0; i<=nBinsCharge; i++) binLimCharge[i]=(Double_t)minCharge + (maxCharge-minCharge)/nBinsCharge*(Double_t)i ;
 860
 861   Int_t nBinsSource = 10;
 862   Double_t minSource = 0.;
 863   Double_t maxSource = 10.;
 864   Double_t binLimSource[nBinsSource+1];
 865   for(Int_t i=0; i<=nBinsSource; i++) binLimSource[i]=(Double_t)minSource + (maxSource-minSource)/nBinsSource*(Double_t)i ;
 866
 867   Int_t nBinsInvMass = 50;
 868   Double_t minInvMass = 0.;
 869   Double_t maxInvMass = 0.3;
 870   Double_t binLimInvMass[nBinsInvMass+1];
 871   for(Int_t i=0; i<=nBinsInvMass; i++) binLimInvMass[i]=(Double_t)minInvMass + (maxInvMass-minInvMass)/nBinsInvMass*(Double_t)i ;
 872
 873   Int_t nBinsMult = 100;
 874   Double_t minMult = 0.;
 875   Double_t maxMult = 3000;
 876   Double_t binLimMult[nBinsMult+1];
 877   //for(Int_t i=0; i<=nBinsMult; i++) binLimMult[i]=TMath::Power((Double_t)minMult + (TMath::Sqrt(maxMult)-TMath::Sqrt(minMult))/nBinsMult*(Double_t)i,2);
 878   for(Int_t i=0; i<=nBinsMult; i++) binLimMult[i]=(Double_t)minMult + (maxMult-minMult)/nBinsMult*(Double_t)i;
 879
 880   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: variables");
 881
 882   //******************
 883   // Histograms
 884   //******************
 885
 886   fListHist = new TList();
 887   fListHist->SetOwner();
 888
 889   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: list");
 890
 891   // Minimum histos
 892
 893   // Histos
 894   fHistEV = new TH2D("fHistEV", "events", 3, 0, 3, 3, 0,3);
 895
 896   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: histev");
 897
 898   // V0 multiplicity vs # of tracks vs centraliy
 899   const Int_t nDimPU=4;
 900   Int_t nBinPU[nDimPU] = {nBinsCEvenMore,nBinsCEvenMore,nBinsMult,nBinsMult};
 901   fHistPileUp = new THnSparseF("PileUp","PileUp",nDimPU,nBinPU);
 902   fHistPileUp->SetBinEdges(0,binLimCEvenMore);
 903   fHistPileUp->SetBinEdges(1,binLimCEvenMore);
 904   fHistPileUp->SetBinEdges(2,binLimMult);
 905   fHistPileUp->SetBinEdges(3,binLimMult);
 906   fHistPileUp->Sumw2();
 907
 908   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: eventplane");
 909
 910   // Event plane as function of phiep, centrality
 911   const Int_t nDima=4;
 912   Int_t nBina[nDima] = {nBinsPhi,nBinsPhi,nBinsPhi,nBinsC};
 913   fEventPlane = new THnSparseF("EventPlane","EventPlane",nDima,nBina);
 914   fEventPlane->SetBinEdges(0,binLimPhi);
 915   fEventPlane->SetBinEdges(1,binLimPhi);
 916   fEventPlane->SetBinEdges(2,binLimPhi);
 917   fEventPlane->SetBinEdges(3,binLimC);
 918   fEventPlane->Sumw2();
 919
 920   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: eventplane");
 921
 922   // Fraction of contamination, centrality
 923   const Int_t nDimcont=2;
 924   Int_t nBincont[nDimcont] = {fPtBinning.GetSize()-1,nBinsC};
 925   fFractionContamination = new THnSparseF("Contamination","Contamination",nDimcont,nBincont);
 926   fFractionContamination->SetBinEdges(0,fPtBinning.GetArray());
 927   fFractionContamination->SetBinEdges(1,binLimC);
 928   fFractionContamination->Sumw2();
 929   //
 930   fContaminationv2 = new TProfile2D("Contaminationv2","",nBinsC,binLimC,fPtBinning.GetSize()-1,fPtBinning.GetArray());
 931   fContaminationv2->Sumw2();
 932   //
 933   fContaminationmeanpt = new TProfile2D("Contaminationmeanpt","",nBinsC,binLimC,fPtBinning.GetSize()-1,fPtBinning.GetArray());
 934   fContaminationmeanpt->Sumw2();
 935
 936   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: fraction of contamination");
 937
 938   // Resolution cosres_abc centrality
 939   const Int_t nDimfbis=4;
 940   Int_t nBinfbis[nDimfbis] = {nBinsCos,nBinsCos,nBinsCos,nBinsCMore};
 941   fCosResabc = new THnSparseF("CosRes_abc","CosRes_abc",nDimfbis,nBinfbis);
 942   fCosResabc->SetBinEdges(0,binLimCos);
 943   fCosResabc->SetBinEdges(1,binLimCos);
 944   fCosResabc->SetBinEdges(2,binLimCos);
 945   fCosResabc->SetBinEdges(3,binLimCMore);
 946   fCosResabc->Sumw2();
 947
 948   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: cosresabc");
 949
 950   // Resolution cosres centrality
 951   const Int_t nDimf=2;
 952   Int_t nBinf[nDimf] = {nBinsCos, nBinsCMore};
 953   fCosRes = new THnSparseF("CosRes","CosRes",nDimf,nBinf);
 954   fCosRes->SetBinEdges(0,binLimCos);
 955   fCosRes->SetBinEdges(1,binLimCMore);
 956   fCosRes->Sumw2();
 957
 958   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: cosres");
 959
 960   // Maps delta phi
 961   const Int_t nDimg=5;
 962   Int_t nBing[nDimg] = {nBinsPhi,nBinsC,fPtBinning.GetSize()-1, nBinsCharge,nBinsEtaLess};
 963   fDeltaPhiMaps = new THnSparseF("DeltaPhiMaps","DeltaPhiMaps",nDimg,nBing);
 964   fDeltaPhiMaps->SetBinEdges(0,binLimPhi);
 965   fDeltaPhiMaps->SetBinEdges(1,binLimC);
 966   fDeltaPhiMaps->SetBinEdges(2,fPtBinning.GetArray());
 967   fDeltaPhiMaps->SetBinEdges(3,binLimCharge);
 968   fDeltaPhiMaps->SetBinEdges(4,binLimEtaLess);
 969   fDeltaPhiMaps->Sumw2();
 970
 971   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: deltaphimaps");
 972
 973   // Maps cos phi
 974   const Int_t nDimh=5;
 975   Int_t nBinh[nDimh] = {nBinsCos,nBinsC,fPtBinning.GetSize()-1,nBinsCharge,nBinsEtaLess};
 976   fCosPhiMaps = new THnSparseF("CosPhiMaps","CosPhiMaps",nDimh,nBinh);
 977   fCosPhiMaps->SetBinEdges(0,binLimCos);
 978   fCosPhiMaps->SetBinEdges(1,binLimC);
 979   fCosPhiMaps->SetBinEdges(2,fPtBinning.GetArray());
 980   fCosPhiMaps->SetBinEdges(3,binLimCharge);
 981   fCosPhiMaps->SetBinEdges(4,binLimEtaLess);
 982   fCosPhiMaps->Sumw2();
 983
 984   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: cosphimaps");
 985
 986   //
 987   // fMonitorEventPlane
 988   //
 989   //
 990
 991   if(fMonitorEventPlane) {
 992     // Event Plane after subtraction as function of phiep, centrality, pt, eta
 993     const Int_t nDimb=2;
 994     Int_t nBinb[nDimb] = {nBinsPhi, nBinsC};
 995     fEventPlaneaftersubtraction = new THnSparseF("EventPlane_aftersubtraction","EventPlane_aftersubtraction",nDimb,nBinb);
 996     fEventPlaneaftersubtraction->SetBinEdges(0,binLimPhi);
 997     fEventPlaneaftersubtraction->SetBinEdges(1,binLimC);
 998     fEventPlaneaftersubtraction->Sumw2();
 999
1000     AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: eventplane after sub");
1001
1002     // Monitoring of the event Plane cos(2phi) sin(2phi) centrality
1003     const Int_t nDimi=3;
1004     Int_t nBini[nDimi] = {nBinsCos, nBinsCos, nBinsCMore};
1005     fCosSin2phiep = new THnSparseF("CosSin2phiep","CosSin2phiep",nDimi,nBini);
1006     fCosSin2phiep->SetBinEdges(0,binLimCos);
1007     fCosSin2phiep->SetBinEdges(1,binLimCos);
1008     fCosSin2phiep->SetBinEdges(2,binLimCMore);
1009     fCosSin2phiep->Sumw2();
1010
1011     AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: cossin2phiep");
1012
1013     // Monitoring Event plane after subtraction of the track
1014     const Int_t nDime=4;
1015     Int_t nBine[nDime] = {nBinsCos, nBinsC, fPtBinning.GetSize()-1, nBinsEta};
1016     fCos2phie = new THnSparseF("cos2phie","cos2phie",nDime,nBine);
1017     fCos2phie->SetBinEdges(2,fPtBinning.GetArray());
1018     fCos2phie->SetBinEdges(3,binLimEta);
1019     fCos2phie->SetBinEdges(0,binLimCos);
1020     fCos2phie->SetBinEdges(1,binLimC);
1021     fCos2phie->Sumw2();
1022     AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: cos2phie");
1023     fSin2phie = new THnSparseF("sin2phie","sin2phie",nDime,nBine);
1024     fSin2phie->SetBinEdges(2,fPtBinning.GetArray());
1025     fSin2phie->SetBinEdges(3,binLimEta);
1026     fSin2phie->SetBinEdges(0,binLimCos);
1027     fSin2phie->SetBinEdges(1,binLimC);
1028     fSin2phie->Sumw2();
1029     AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: sin2phie");
1030     fCos2phiep = new THnSparseF("cos2phiep","cos2phiep",nDime,nBine);
1031     fCos2phiep->SetBinEdges(2,fPtBinning.GetArray());
1032     fCos2phiep->SetBinEdges(3,binLimEta);
1033     fCos2phiep->SetBinEdges(0,binLimCos);
1034     fCos2phiep->SetBinEdges(1,binLimC);
1035     fCos2phiep->Sumw2();
1036     AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: cos2phiep");
1037     fSin2phiep = new THnSparseF("sin2phiep","sin2phiep",nDime,nBine);
1038     fSin2phiep->SetBinEdges(2,fPtBinning.GetArray());
1039     fSin2phiep->SetBinEdges(3,binLimEta);
1040     fSin2phiep->SetBinEdges(0,binLimCos);
1041     fSin2phiep->SetBinEdges(1,binLimC);
1042     fSin2phiep->Sumw2();
1043     AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: sin2phiep");
1044     fSin2phiephiep = new THnSparseF("sin2phie_phiep","sin2phie_phiep",nDime,nBine);
1045     fSin2phiephiep->SetBinEdges(2,fPtBinning.GetArray());
1046     fSin2phiephiep->SetBinEdges(3,binLimEta);
1047     fSin2phiephiep->SetBinEdges(0,binLimCos);
1048     fSin2phiephiep->SetBinEdges(1,binLimC);
1049     fSin2phiephiep->Sumw2();
1050     AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: sin2phiephiep");
1051
1052     const Int_t nDimfbiss=4;
1053     Int_t nBinfbiss[nDimfbiss] = {nBinsCos,nBinsCos,nBinsCos,nBinsC};
1054     fSinResabc = new THnSparseF("SinRes_abc","SinRes_abc",nDimfbiss,nBinfbiss);
1055     fSinResabc->SetBinEdges(0,binLimCos);
1056     fSinResabc->SetBinEdges(1,binLimCos);
1057     fSinResabc->SetBinEdges(2,binLimCos);
1058     fSinResabc->SetBinEdges(3,binLimC);
1059     fSinResabc->Sumw2();
1060     AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: sinresabc");
1061
1062     // Profile cosres centrality with 3 subevents
1063     fProfileCosResab = new TProfile("ProfileCosRes_a_b","ProfileCosRes_a_b",nBinsCMore,binLimCMore);
1064     fProfileCosResab->Sumw2();
1065     fProfileCosResac = new TProfile("ProfileCosRes_a_c","ProfileCosRes_a_c",nBinsCMore,binLimCMore);
1066     fProfileCosResac->Sumw2();
1067     fProfileCosResbc = new TProfile("ProfileCosRes_b_c","ProfileCosRes_b_c",nBinsCMore,binLimCMore);
1068     fProfileCosResbc->Sumw2();
1069     AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: profilecosresbc");
1070
1071     //
1072     const Int_t nDimff=2;
1073     Int_t nBinff[nDimff] = {nBinsCos, nBinsC};
1074     fSinRes = new THnSparseF("SinRes","SinRes",nDimff,nBinff);
1075     fSinRes->SetBinEdges(0,binLimCos);
1076     fSinRes->SetBinEdges(1,binLimC);
1077     fSinRes->Sumw2();
1078     AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: sinres");
1079
1080     // Profile cosres centrality
1081     fProfileCosRes = new TProfile("ProfileCosRes","ProfileCosRes",nBinsCMore,binLimCMore);
1082     fProfileCosRes->Sumw2();
1083     AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: profilecosres");
1084
1085     // Profile Maps cos phi
1086     fProfileCosPhiMaps = new TProfile2D("ProfileCosPhiMaps","ProfileCosPhiMaps",nBinsC,binLimC,fPtBinning.GetSize()-1,fPtBinning.GetArray());
1087     fProfileCosPhiMaps->Sumw2();
1088     AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: profilecosphimaps");
1089
1090   }
1091   //
1092   // fMonitorTrackCuts
1093   //
1094
1095   if(fMonitorTrackCuts) {
1096     // Debugging tracking steps
1097     const Int_t nDimTrStep=2;
1098     Int_t nBinTrStep[nDimTrStep] = {fPtBinning.GetSize()-1,nBinsStep};
1099     fTrackingCuts = new THnSparseF("TrackingCuts","TrackingCuts",nDimTrStep,nBinTrStep);
1100     fTrackingCuts->SetBinEdges(0,fPtBinning.GetArray());
1101     fTrackingCuts->SetBinEdges(1,binLimStep);
1102     fTrackingCuts->Sumw2();
1103     AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: trackingcuts");
1104   }
1105
1106   //
1107   // fMonitorContamination
1108   //
1109
1110   if(fMonitorContamination) {
1111     // Maps delta phi contamination
1112     const Int_t nDimgcont=4;
1113     Int_t nBingcont[nDimgcont] = {nBinsPhiLess,nBinsC,fPtBinning.GetSize()-1, nBinsTPCdEdx};
1114     fDeltaPhiMapsContamination = new THnSparseF("DeltaPhiMapsContamination","DeltaPhiMapsContamination",nDimgcont,nBingcont);
1115     fDeltaPhiMapsContamination->SetBinEdges(0,binLimPhiLess);
1116     fDeltaPhiMapsContamination->SetBinEdges(1,binLimC);
1117     fDeltaPhiMapsContamination->SetBinEdges(2,fPtBinning.GetArray());
1118     fDeltaPhiMapsContamination->SetBinEdges(3,binLimTPCdEdx);
1119     fDeltaPhiMapsContamination->Sumw2();
1120     AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: deltaphimapscontamination");
1121
1122   }
1123   //
1124   // fMonitorWithoutPID
1125   //
1126
1127   if(fMonitorWithoutPID) {
1128     //
1129     const Int_t nDimgb=3;
1130     Int_t nBingb[nDimgb] = {nBinsPhi,nBinsC,fPtBinning.GetSize()-1};
1131
1132     fDeltaPhiMapsBeforePID = new THnSparseF("DeltaPhiMapsBeforePID","DeltaPhiMapsBeforePID",nDimgb,nBingb);
1133     fDeltaPhiMapsBeforePID->SetBinEdges(0,binLimPhi);
1134     fDeltaPhiMapsBeforePID->SetBinEdges(1,binLimC);
1135     fDeltaPhiMapsBeforePID->SetBinEdges(2,fPtBinning.GetArray());
1136     fDeltaPhiMapsBeforePID->Sumw2();
1137     AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: deltaphimapsbeforepid");
1138
1139     const Int_t nDimhb=3;
1140     Int_t nBinhb[nDimhb] = {nBinsCos,nBinsC,fPtBinning.GetSize()-1};
1141
1142     fCosPhiMapsBeforePID = new THnSparseF("CosPhiMapsBeforePID","CosPhiMapsBeforePID",nDimhb,nBinhb);
1143     fCosPhiMapsBeforePID->SetBinEdges(0,binLimCos);
1144     fCosPhiMapsBeforePID->SetBinEdges(1,binLimC);
1145     fCosPhiMapsBeforePID->SetBinEdges(2,fPtBinning.GetArray());
1146     fCosPhiMapsBeforePID->Sumw2();
1147     AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: cosphimapsbeforepid");
1148   }
1149   //
1150   // fMonitorPhotonic
1151   //
1152
1153   if(fMonitorPhotonic) {
1154
1155     const Int_t nDimgbp=3;
1156     Int_t nBingbp[nDimgbp] = {nBinsPhi,nBinsC,fPtBinning.GetSize()-1};
1157
1158     fDeltaPhiMapsTaggedPhotonic = new THnSparseF("DeltaPhiMapsTaggedPhotonic","DeltaPhiMapsTaggedPhotonic",nDimgbp,nBingbp);
1159     fDeltaPhiMapsTaggedPhotonic->SetBinEdges(0,binLimPhi);
1160     fDeltaPhiMapsTaggedPhotonic->SetBinEdges(1,binLimC);
1161     fDeltaPhiMapsTaggedPhotonic->SetBinEdges(2,fPtBinning.GetArray());
1162     fDeltaPhiMapsTaggedPhotonic->Sumw2();
1163     AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: deltaphimapstaggedphotonic");
1164
1165     fDeltaPhiMapsTaggedNonPhotonic = new THnSparseF("DeltaPhiMapsTaggedNonPhotonic","DeltaPhiMapsTaggedNonPhotonic",nDimgbp,nBingbp);
1166     fDeltaPhiMapsTaggedNonPhotonic->SetBinEdges(0,binLimPhi);
1167     fDeltaPhiMapsTaggedNonPhotonic->SetBinEdges(1,binLimC);
1168     fDeltaPhiMapsTaggedNonPhotonic->SetBinEdges(2,fPtBinning.GetArray());
1169     fDeltaPhiMapsTaggedNonPhotonic->Sumw2();
1170     AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: deltaphimapstaggednonphotonic");
1171
1172     fDeltaPhiMapsTaggedPhotonicLS = new THnSparseF("DeltaPhiMapsTaggedPhotonicLS","DeltaPhiMapsTaggedPhotonicLS",nDimgbp,nBingbp);
1173     fDeltaPhiMapsTaggedPhotonicLS->SetBinEdges(0,binLimPhi);
1174     fDeltaPhiMapsTaggedPhotonicLS->SetBinEdges(1,binLimC);
1175     fDeltaPhiMapsTaggedPhotonicLS->SetBinEdges(2,fPtBinning.GetArray());
1176     fDeltaPhiMapsTaggedPhotonicLS->Sumw2();
1177     AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: deltaphimapstaggedphotonicls");
1178
1179     const Int_t nDimMCSource=3;
1180     Int_t nBinMCSource[nDimMCSource] = {nBinsC,fPtBinning.GetSize()-1,nBinsSource};
1181     fMCSourceDeltaPhiMaps = new THnSparseF("MCSourceDeltaPhiMaps","MCSourceDeltaPhiMaps",nDimMCSource,nBinMCSource);
1182     fMCSourceDeltaPhiMaps->SetBinEdges(0,binLimC);
1183     fMCSourceDeltaPhiMaps->SetBinEdges(1,fPtBinning.GetArray());
1184     fMCSourceDeltaPhiMaps->SetBinEdges(2,binLimSource);
1185     fMCSourceDeltaPhiMaps->Sumw2();
1186     AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: mcsourcedeltaphimaps");
1187
1188     // Maps invmass opposite
1189     const Int_t nDimOppSign=5;
1190     Int_t nBinOppSign[nDimOppSign] = {nBinsPhi,nBinsC,fPtBinning.GetSize()-1,nBinsInvMass,nBinsSource};
1191     fOppSignDeltaPhiMaps = new THnSparseF("OppSignDeltaPhiMaps","OppSignDeltaPhiMaps",nDimOppSign,nBinOppSign);
1192     fOppSignDeltaPhiMaps->SetBinEdges(0,binLimPhi);
1193     fOppSignDeltaPhiMaps->SetBinEdges(1,binLimC);
1194     fOppSignDeltaPhiMaps->SetBinEdges(2,fPtBinning.GetArray());
1195     fOppSignDeltaPhiMaps->SetBinEdges(3,binLimInvMass);
1196     fOppSignDeltaPhiMaps->SetBinEdges(4,binLimSource);
1197     fOppSignDeltaPhiMaps->Sumw2();
1198     AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: oppsigndeltaphimaps");
1199
1200     // Maps invmass same sign
1201     const Int_t nDimSameSign=5;
1202     Int_t nBinSameSign[nDimSameSign] = {nBinsPhi,nBinsC,fPtBinning.GetSize()-1,nBinsInvMass,nBinsSource};
1203     fSameSignDeltaPhiMaps = new THnSparseF("SameSignDeltaPhiMaps","SameSignDeltaPhiMaps",nDimSameSign,nBinSameSign);
1204     fSameSignDeltaPhiMaps->SetBinEdges(0,binLimPhi);
1205     fSameSignDeltaPhiMaps->SetBinEdges(1,binLimC);
1206     fSameSignDeltaPhiMaps->SetBinEdges(2,fPtBinning.GetArray());
1207     fSameSignDeltaPhiMaps->SetBinEdges(3,binLimInvMass);
1208     fSameSignDeltaPhiMaps->SetBinEdges(4,binLimSource);
1209     fSameSignDeltaPhiMaps->Sumw2();
1210     AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: samesigndeltaphimaps");
1211
1212     // Maps angle same sign
1213     const Int_t nDimAngleSameSign=3;
1214     Int_t nBinAngleSameSign[nDimAngleSameSign] = {nBinsAngle,nBinsC,nBinsSource};
1215     fSameSignAngle = new THnSparseF("SameSignAngleMaps","SameSignAngleMaps",nDimAngleSameSign,nBinAngleSameSign);
1216     fSameSignAngle->SetBinEdges(0,binLimAngle);
1217     fSameSignAngle->SetBinEdges(1,binLimC);
1218     fSameSignAngle->SetBinEdges(2,binLimSource);
1219     fSameSignAngle->Sumw2();
1220     AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: samesignangle");
1221
1222     // Maps angle opp sign
1223     const Int_t nDimAngleOppSign=3;
1224     Int_t nBinAngleOppSign[nDimAngleOppSign] = {nBinsAngle,nBinsC,nBinsSource};
1225     fOppSignAngle = new THnSparseF("OppSignAngleMaps","OppSignAngleMaps",nDimAngleOppSign,nBinAngleOppSign);
1226     fOppSignAngle->SetBinEdges(0,binLimAngle);
1227     fOppSignAngle->SetBinEdges(1,binLimC);
1228     fOppSignAngle->SetBinEdges(2,binLimSource);
1229     fOppSignAngle->Sumw2();
1230     AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: oppsignangle");
1231
1232   }
1233
1234   //**************************
1235   // Add to the list
1236   //******************************
1237
1238   fListHist->Add(fHistEV);
1239   fListHist->Add(fHistPileUp);
1240   fListHist->Add(fEventPlane);
1241   fListHist->Add(fFractionContamination);
1242   fListHist->Add(fCosRes);
1243   fListHist->Add(fCosResabc);
1244   fListHist->Add(fCosPhiMaps);
1245   fListHist->Add(fDeltaPhiMaps);
1246   fListHist->Add(fPIDqa->MakeList("HFEpidQA"));
1247   fListHist->Add(fContaminationv2);
1248   fListHist->Add(fContaminationmeanpt);
1249   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: add default");
1250
1251   if(fMonitorEventPlane) {
1252     fListHist->Add(fProfileCosRes);
1253     fListHist->Add(fProfileCosResab);
1254     fListHist->Add(fProfileCosResac);
1255     fListHist->Add(fProfileCosResbc);
1256     fListHist->Add(fCosSin2phiep);
1257     fListHist->Add(fCos2phie);
1258     fListHist->Add(fSin2phie);
1259     fListHist->Add(fCos2phiep);
1260     fListHist->Add(fSin2phiep);
1261     fListHist->Add(fSin2phiephiep);
1262     fListHist->Add(fSinRes);
1263     fListHist->Add(fSinResabc);
1264     fListHist->Add(fProfileCosPhiMaps);
1265   }
1266   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: add monitor");
1267
1268   if(fMonitorTrackCuts) fListHist->Add(fTrackingCuts);
1269
1270   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: add monitortrackcuts");
1271
1272   if(fMonitorContamination) {
1273     fListHist->Add(fDeltaPhiMapsContamination);
1274   }
1275
1276   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: add deltaphimapscontamination");
1277
1278   if(fMonitorWithoutPID) {
1279     fListHist->Add(fDeltaPhiMapsBeforePID);
1280     fListHist->Add(fCosPhiMapsBeforePID);
1281   }
1282
1283   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: add without pid");
1284
1285   if(fMonitorPhotonic) {
1286   fListHist->Add(fPIDBackgroundqa->MakeList("HFEpidBackgroundQA"));
1287   fListHist->Add(fDeltaPhiMapsTaggedPhotonic);
1288   fListHist->Add(fDeltaPhiMapsTaggedNonPhotonic);
1289   fListHist->Add(fDeltaPhiMapsTaggedPhotonicLS);
1290   fListHist->Add(fMCSourceDeltaPhiMaps);
1291   fListHist->Add(fOppSignDeltaPhiMaps);
1292   fListHist->Add(fSameSignDeltaPhiMaps);
1293   fListHist->Add(fSameSignAngle);
1294   fListHist->Add(fOppSignAngle);
1295   fListHist->Add(fBackgroundSubtraction->GetListOutput());
1296   }
1297
1298   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: add photonic");
1299
1300   if(fHFEVZEROEventPlane && fMonitorEventPlane) fListHist->Add(fHFEVZEROEventPlane->GetOutputList());
1301
1302   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: add event plane");
1303
1304   fListHist->Print();
1305
1306   PostData(1, fListHist);
1307   //for(Int_t bincless = 0; bincless < fNbBinsCentralityQCumulant; bincless++) {
1308   // PostData(bincless+2,fflowEvent);
1309   //}
1310
1311   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: post");
1312
1313
1314 }
1315
1316 //________________________________________________________________________
1317 void AliAnalysisTaskFlowTPCTOFEPSP::UserExec(Option_t */*option*/)
1318 {
1319   //
1320   // Loop over event
1321   //
1322
1323
1324   Double_t massElectron = 0.000511;
1325   Double_t mcReactionPlane = 0.0;
1326
1327   Float_t cntr = 0.0;
1328   Double_t binct = 11.5;
1329   Double_t binctMore = 20.5;
1330   Double_t binctLess = -0.5;
1331   Float_t binctt = -1.0;
1332
1333   Double_t valuecossinephiep[3];
1334   Double_t valuensparsea[4];
1335   Double_t valuensparseabis[5];
1336   Double_t valuensparsee[4];
1337   Double_t valuensparsef[2];
1338   Double_t valuensparsefsin[2];
1339   Double_t valuensparsefbis[4];
1340   Double_t valuensparsefbissin[4];
1341   Double_t valuensparseg[5];
1342   Double_t valuensparseh[5];
1343   Double_t valuensparsehprofile[3];
1344   Double_t valuensparseMCSourceDeltaPhiMaps[3];
1345   Double_t valuetrackingcuts[2];
1346   Double_t valuedeltaphicontamination[4];
1347   Double_t valuefractioncont[2];
1348
1349   AliMCEvent *mcEvent = MCEvent();
1350   AliMCParticle *mctrack = NULL;
1351
1352   // MC info
1353   Bool_t mcthere = kTRUE;
1354   if(fAODAnalysis) {
1355     AliAODEvent *aodE = dynamic_cast<AliAODEvent *>(fInputEvent);
1356     if(!aodE){
1357       //        printf("testd\n");
1358       AliError("No AOD Event");
1359       return;
1360     }
1361     fAODMCHeader = dynamic_cast<AliAODMCHeader *>(fInputEvent->FindListObject(AliAODMCHeader::StdBranchName()));
1362     if(!fAODMCHeader){
1363       mcthere = kFALSE;
1364     }
1365     fAODArrayMCInfo = dynamic_cast<TClonesArray *>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
1366     if(!fAODArrayMCInfo){
1367       mcthere = kFALSE;
1368     }
1369     else {
1370       fHFECuts->SetMCEvent(aodE);
1371       if(fMonitorPhotonic) fBackgroundSubtraction->SetAODArrayMCInfo(fAODArrayMCInfo);
1372     }
1373   }
1374   else {
1375     AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
1376     if(!mcH) mcthere = kFALSE;
1377     else {
1378       if(fMonitorPhotonic) fBackgroundSubtraction->SetMCEvent(fMCEvent);
1379     }
1380   }
1381
1382
1383   /////////////////
1384   // centrality
1385   /////////////////
1386
1387   //AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
1388   //if(!esd) return;
1389   AliCentrality *centrality = fInputEvent->GetCentrality();
1390   //AliDebug(2,"Got the centrality");
1391   if(!centrality) return;
1392   cntr = centrality->GetCentralityPercentile("V0M");
1393   if((0.0< cntr) && (cntr<5.0)) binct = 0.5;
1394   if((5.0< cntr) && (cntr<10.0)) binct = 1.5;
1395   if((10.0< cntr) && (cntr<20.0)) binct = 2.5;
1396   if((20.0< cntr) && (cntr<30.0)) binct = 3.5;
1397   if((30.0< cntr) && (cntr<40.0)) binct = 4.5;
1398   if((40.0< cntr) && (cntr<50.0)) binct = 5.5;
1399   if((50.0< cntr) && (cntr<60.0)) binct = 6.5;
1400   if((60.0< cntr) && (cntr<70.0)) binct = 7.5;
1401   if((70.0< cntr) && (cntr<80.0)) binct = 8.5;
1402   if((80.0< cntr) && (cntr<90.0)) binct = 9.5;
1403   if((90.0< cntr) && (cntr<100.0)) binct = 10.5;
1404
1405   if((0.< cntr) && (cntr < 20.)) binctt = 0.5;
1406   if((20.< cntr) && (cntr < 40.)) binctt = 1.5;
1407   if((40.< cntr) && (cntr < 80.)) binctt = 2.5;
1408
1409   if((0.0< cntr) && (cntr<5.0)) binctMore = 0.5;
1410   if((5.0< cntr) && (cntr<10.0)) binctMore = 1.5;
1411   if((10.0< cntr) && (cntr<15.0)) binctMore = 2.5;
1412   if((15.0< cntr) && (cntr<20.0)) binctMore = 3.5;
1413   if((20.0< cntr) && (cntr<25.0)) binctMore = 4.5;
1414   if((25.0< cntr) && (cntr<30.0)) binctMore = 5.5;
1415   if((30.0< cntr) && (cntr<35.0)) binctMore = 6.5;
1416   if((35.0< cntr) && (cntr<40.0)) binctMore = 7.5;
1417   if((40.0< cntr) && (cntr<45.0)) binctMore = 8.5;
1418   if((45.0< cntr) && (cntr<50.0)) binctMore = 9.5;
1419   if((50.0< cntr) && (cntr<55.0)) binctMore = 10.5;
1420   if((55.0< cntr) && (cntr<60.0)) binctMore = 11.5;
1421   if((60.0< cntr) && (cntr<65.0)) binctMore = 12.5;
1422   if((65.0< cntr) && (cntr<70.0)) binctMore = 13.5;
1423   if((70.0< cntr) && (cntr<75.0)) binctMore = 14.5;
1424   if((75.0< cntr) && (cntr<80.0)) binctMore = 15.5;
1425   if((80.0< cntr) && (cntr<85.0)) binctMore = 16.5;
1426   if((85.0< cntr) && (cntr<90.0)) binctMore = 17.5;
1427   if((90.0< cntr) && (cntr<95.0)) binctMore = 18.5;
1428   if((95.0< cntr) && (cntr<100.0)) binctMore = 19.5;
1429
1430   binctLess = cntr;
1431
1432
1433   if(binct > 11.0) return;
1434
1435   // centrality
1436   valuensparsea[3] = binct;
1437   valuensparseabis[1] = binct;
1438   valuensparsee[1] = binct;
1439   valuensparsef[1] = binctMore;
1440   valuensparsefsin[1] = binct;
1441   valuensparsefbis[3] = binctMore;
1442   valuensparsefbissin[3] = binct;
1443   valuensparseg[1] = binct;
1444   valuensparseh[1] = binct;
1445   valuefractioncont[1] = binct;
1446   valuensparsehprofile[1] = binct;
1447   valuecossinephiep[2] = binctMore;
1448   valuensparseMCSourceDeltaPhiMaps[0] = binct;
1449   valuedeltaphicontamination[1] = binct;
1450
1451   //////////////////////
1452   // run number
1453   //////////////////////
1454
1455   Int_t runnumber = fInputEvent->GetRunNumber();
1456   AliDebug(2,Form("Run number %d",runnumber));
1457
1458   if(!fPID->IsInitialized()){
1459     // Initialize PID with the given run number
1460     fPID->InitializePID(runnumber);
1461   }
1462   if(!fPIDTOFOnly->IsInitialized()){
1463     // Initialize PID with the given run number
1464     fPIDTOFOnly->InitializePID(runnumber);
1465   }
1466
1467   //
1468   if(!fPIDBackground->IsInitialized()){
1469     // Initialize PID with the given run number
1470     fPIDBackground->InitializePID(runnumber);
1471   }
1472
1473   fHFECuts->SetRecEvent(fInputEvent);
1474   if(mcEvent) fHFECuts->SetMCEvent(mcEvent);
1475
1476
1477   //////////
1478   // PID
1479   //////////
1480
1481   AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
1482   if(!pidResponse){
1483     AliDebug(2,"No PID response set");
1484     return;
1485   }
1486   fPID->SetPIDResponse(pidResponse);
1487   fPIDTOFOnly->SetPIDResponse(pidResponse);
1488   fPIDBackground->SetPIDResponse(pidResponse);
1489   if(fMonitorPhotonic) fBackgroundSubtraction->InitRun(fInputEvent,pidResponse);
1490
1491   fHistEV->Fill(binctt,0.0);
1492
1493   //////////////////
1494   // Event cut
1495   //////////////////
1496   if(!fHFECuts->CheckEventCuts("fEvRecCuts", fInputEvent)) {
1497     AliDebug(2,"Does not pass the event cut");
1498     PostData(1, fListHist);
1499     return;
1500   }
1501
1502   fHistEV->Fill(binctt,1.0);
1503
1504
1505   ///////////////////////////////////////////////////////////
1506   // PileUpCut
1507   ///////////////////////////////////////////////////////////
1508
1509   Float_t multTPC(0.); // tpc mult estimate
1510   Float_t multGlob(0.); // global multiplicity
1511   const Int_t nGoodTracks = fInputEvent->GetNumberOfTracks();
1512   for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
1513     AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*>(fInputEvent->GetTrack(iTracks));
1514     if (!trackAOD) continue;
1515     if (!(trackAOD->TestFilterBit(1))) continue;
1516     if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70)  || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.2)) continue;
1517     multTPC++;
1518   }
1519   for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult
1520     AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*>(fInputEvent->GetTrack(iTracks));
1521     if (!trackAOD) continue;
1522     if (!(trackAOD->TestFilterBit(16))) continue;
1523     if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70) || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.1)) continue;
1524     Double_t b[2] = {-99., -99.};
1525     Double_t bCov[3] = {-99., -99., -99.};
1526     if (!(trackAOD->PropagateToDCA(fInputEvent->GetPrimaryVertex(), fInputEvent->GetMagneticField(), 100., b, bCov))) continue;
1527     if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
1528     multGlob++;
1529   } //track loop
1530
1531   Double_t pileup[4];
1532   pileup[0]=fInputEvent->GetCentrality()->GetCentralityPercentile("V0M");
1533   pileup[1]=fInputEvent->GetCentrality()->GetCentralityPercentile("TRK");
1534   pileup[2]=multTPC;
1535   pileup[3]=multGlob;
1536   fHistPileUp->Fill(pileup);
1537
1538   if(fPileUpCut){
1539     if (TMath::Abs(pileup[0]-pileup[1]) > 5) {
1540       AliDebug(2,"Does not pass the centrality correlation cut");
1541       return;
1542     }
1543     if(multTPC < (-36.81+1.48*multGlob) && multTPC > (63.03+1.78*multGlob)){
1544       AliDebug(2,"Does not pass the multiplicity correlation cut");
1545       return;
1546     }
1547   }
1548
1549   // AliVVZERO* vzeroData=fInputEvent->GetVZEROData();
1550   // Double_t mult[3],multV0A(0),multV0C(0);
1551   // for(Int_t i=0; i<32; ++i) {
1552   //   multV0A += vzeroData->GetMultiplicityV0A(i);
1553   //   multV0C += vzeroData->GetMultiplicityV0C(i);
1554   // }
1555
1556   // int ntrk=0;
1557   // for(Int_t k = 0; k < fInputEvent->GetNumberOfTracks(); k++){
1558   //   AliVTrack *track = (AliVTrack *) fInputEvent->GetTrack(k);
1559   //   if(!track) continue;
1560   //   if(!(track->GetStatus()&AliVTrack::kITSrefit)) continue;
1561   //   if(!(track->GetStatus()&AliVTrack::kTPCrefit)) continue;
1562   //   ntrk++;
1563   // }
1564
1565   // mult[0]=fInputEvent->GetNumberOfTracks();
1566   // mult[1]=multV0A+multV0C;
1567   // mult[2]=binctMore;
1568   // fHistPileUp->Fill(mult);
1569
1570   // if(fUpperPileUpCut&&fLowerPileUpCut){
1571   //   if((mult[0]<fLowerPileUpCut->Eval(mult[1])) ||
1572   //      (mult[0]>fUpperPileUpCut->Eval(mult[1]))){
1573   //     AliDebug(2,"Does not pass the pileup cut");
1574   //     PostData(1, fListHist);
1575   //     return;
1576   //   }
1577   // }
1578
1579   ////////////////////////////////////
1580   // First method event plane
1581   ////////////////////////////////////
1582
1583   AliEventplane* vEPa = fInputEvent->GetEventplane();
1584   Float_t eventPlanea = 0.0;
1585   Float_t eventPlaneTPC = 0.0;
1586   Float_t eventPlaneV0A = 0.0;
1587   Float_t eventPlaneV0C = 0.0;
1588   Float_t eventPlaneV0 = 0.0;
1589   TVector2 *qTPC = 0x0;
1590   TVector2 *qsub1a = 0x0;
1591   TVector2 *qsub2a = 0x0;
1592   TVector2 qV0A,qV0C,qV0,*qAna;
1593
1594   // V0
1595
1596   if(fHFEVZEROEventPlane && (!fAODAnalysis)){
1597
1598     //AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
1599     //if(!esd) return;
1600
1601     fHFEVZEROEventPlane->ProcessEvent(fInputEvent);
1602
1603     if(TMath::Abs(fHFEVZEROEventPlane->GetEventPlaneV0A()+100) < 0.0000001) eventPlaneV0A = -100.0;
1604     else {
1605       eventPlaneV0A = TVector2::Phi_0_2pi(fHFEVZEROEventPlane->GetEventPlaneV0A());
1606       if(eventPlaneV0A > TMath::Pi()) eventPlaneV0A = eventPlaneV0A - TMath::Pi();
1607     }
1608
1609     if(TMath::Abs(fHFEVZEROEventPlane->GetEventPlaneV0C()+100) < 0.0000001) eventPlaneV0C = -100.0;
1610     else {
1611       eventPlaneV0C = TVector2::Phi_0_2pi(fHFEVZEROEventPlane->GetEventPlaneV0C());
1612       if(eventPlaneV0C > TMath::Pi()) eventPlaneV0C = eventPlaneV0C - TMath::Pi();
1613     }
1614
1615     if(TMath::Abs(fHFEVZEROEventPlane->GetEventPlaneV0()+100) < 0.0000001) eventPlaneV0 = -100.0;
1616     else {
1617       eventPlaneV0 = TVector2::Phi_0_2pi(fHFEVZEROEventPlane->GetEventPlaneV0());
1618       if(eventPlaneV0 > TMath::Pi()) eventPlaneV0 = eventPlaneV0 - TMath::Pi();
1619     }
1620
1621   }
1622   else {
1623
1624     Double_t qVx, qVy;  //TR: info
1625     eventPlaneV0 = TVector2::Phi_0_2pi(vEPa->CalculateVZEROEventPlane(fInputEvent,10,2,qVx,qVy));
1626     if(eventPlaneV0 > TMath::Pi()) eventPlaneV0 = eventPlaneV0 - TMath::Pi();
1627     qV0.Set(qVx,qVy);
1628     eventPlaneV0A = TVector2::Phi_0_2pi(vEPa->CalculateVZEROEventPlane(fInputEvent,8,2,qVx,qVy));
1629     if(eventPlaneV0A > TMath::Pi()) eventPlaneV0A = eventPlaneV0A - TMath::Pi();
1630     qV0A.Set(qVx,qVy);
1631     eventPlaneV0C = TVector2::Phi_0_2pi(vEPa->CalculateVZEROEventPlane(fInputEvent,9,2,qVx,qVy));
1632     if(eventPlaneV0C > TMath::Pi()) eventPlaneV0C = eventPlaneV0C - TMath::Pi();
1633     qV0C.Set(qVx,qVy);
1634
1635     if(eventPlaneV0<-900) return;
1636     if(eventPlaneV0A<-900) return;
1637     if(eventPlaneV0C<-900) return;
1638
1639
1640     eventPlaneV0=TVector2::Phi_0_2pi(eventPlaneV0);
1641     eventPlaneV0A=TVector2::Phi_0_2pi(eventPlaneV0A);
1642     eventPlaneV0C=TVector2::Phi_0_2pi(eventPlaneV0C);
1643   }
1644
1645
1646   // TPC
1647
1648   qTPC = vEPa->GetQVector();
1649   Double_t qx = -1.0;
1650   Double_t qy = -1.0;
1651   if(qTPC) {
1652     qx = qTPC->X();
1653     qy = qTPC->Y();
1654   }
1655   TVector2 qVectorfortrack;
1656   qVectorfortrack.Set(qx,qy);
1657   eventPlaneTPC = TVector2::Phi_0_2pi(qVectorfortrack.Phi())/2.;
1658
1659   // Choose the one used for v2
1660
1661   if(fVZEROEventPlane){ //TR: info
1662     eventPlanea = eventPlaneV0;
1663     qAna = &qV0;
1664   }
1665   if(fVZEROEventPlaneA){
1666     eventPlanea = eventPlaneV0A;
1667     qAna = &qV0A;
1668   }
1669   if(fVZEROEventPlaneC){
1670     eventPlanea = eventPlaneV0C;
1671     qAna = &qV0C;
1672   }
1673   if(!fVZEROEventPlane){
1674     eventPlanea = eventPlaneTPC;
1675     qAna = &qV0C;
1676   }
1677
1678   valuecossinephiep[0] = TMath::Cos(2*eventPlanea);
1679   valuecossinephiep[1] = TMath::Sin(2*eventPlanea);
1680
1681   Float_t eventPlanesub1a = -100.0;
1682   Float_t eventPlanesub2a = -100.0;
1683   Double_t diffsub1sub2a = -100.0;
1684   Double_t diffsub1sub2asin = -100.0;
1685   Double_t diffsubasubb = -100.0;
1686   Double_t diffsubasubc = -100.0;
1687   Double_t diffsubbsubc = -100.0;
1688   Double_t diffsubasubbsin = -100.0;
1689   Double_t diffsubasubcsin = -100.0;
1690   Double_t diffsubbsubcsin = -100.0;
1691
1692   // two sub event TPC
1693   qsub1a = vEPa->GetQsub1();
1694   qsub2a = vEPa->GetQsub2();
1695
1696   /////////////////////////////////////////////////////////
1697   // Cut for event with event plane reconstructed by all
1698   ////////////////////////////////////////////////////////
1699
1700   if((!qTPC) || (!qsub1a) || (!qsub2a)) {
1701     AliDebug(2,"No event plane");
1702     return;
1703   }
1704
1705   eventPlanesub1a = TVector2::Phi_0_2pi(qsub1a->Phi())/2.;
1706   eventPlanesub2a = TVector2::Phi_0_2pi(qsub2a->Phi())/2.;
1707   diffsub1sub2a = TMath::Cos(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.));
1708   diffsub1sub2asin = TMath::Sin(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.));
1709
1710
1711   // if ( !fDebugStreamer ) {
1712   //   //debug stream
1713   //   TDirectory *backup = gDirectory;
1714   //   fDebugStreamer = new TTreeSRedirector("TaskFlowTPCTOFEPSPdebug.root");
1715   //   if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
1716   // }
1717
1718   // {
1719
1720   //   double v0nrom = TMath::Sqrt(qV0.X()*qV0.X()+qV0.Y()*qV0.Y());
1721   //   double v0Anrom = TMath::Sqrt(qV0A.X()*qV0A.X()+qV0A.Y()*qV0A.Y());
1722   //   double v0Cnrom = TMath::Sqrt(qV0C.X()*qV0C.X()+qV0C.Y()*qV0C.Y());
1723   //   double sub1nrom = TMath::Sqrt(qsub1a->X()*qsub1a->X()+qsub1a->Y()*qsub1a->Y());
1724   //   double sub2nrom = TMath::Sqrt(qsub2a->X()*qsub2a->X()+qsub2a->Y()*qsub2a->Y());
1725
1726   //   (* fDebugStreamer) << "UserExec" <<
1727   //     "binct="<<binct<<
1728   //     "qV0="<<v0nrom<<
1729   //     "qV0A="<<v0Anrom<<
1730   //     "qV0C="<<v0Cnrom<<
1731   //     "qsub1a="<<sub1nrom<<
1732   //     "qsub2a="<<sub2nrom<<
1733   //     "\n";
1734   // }
1735
1736   // three sub events in case of VZEROA and VZEROC
1737   if(!fSP){
1738     diffsubasubb = TMath::Cos(2.*(eventPlaneV0A - eventPlaneV0C));  //TR:
1739     diffsubasubc = TMath::Cos(2.*(eventPlaneV0A - eventPlaneTPC));  //TR:
1740     diffsubbsubc = TMath::Cos(2.*(eventPlaneV0C - eventPlaneTPC));  //TR:
1741   }
1742   else{
1743     if(fVZEROEventPlaneA){
1744       diffsubasubb = qV0A.X()*qV0C.X()+qV0A.Y()*qV0C.Y();
1745       diffsubasubc = qV0A.X()*qTPC->X()+qV0A.Y()*qTPC->Y();
1746       diffsubbsubc = qV0C.X()*qTPC->X()+qV0C.Y()*qTPC->Y();
1747     }
1748     else if(fVZEROEventPlaneC){
1749       diffsubasubb = qV0C.X()*qV0A.X()+qV0C.Y()*qV0A.Y();
1750       diffsubasubc = qV0C.X()*qTPC->X()+qV0C.Y()*qTPC->Y();
1751       diffsubbsubc = qV0A.X()*qTPC->X()+qV0A.Y()*qTPC->Y();
1752     }
1753   }
1754
1755   diffsubasubbsin = TMath::Sin(2.*(eventPlaneV0A - eventPlaneV0C));
1756   diffsubasubcsin = TMath::Sin(2.*(eventPlaneV0A - eventPlaneTPC));
1757   diffsubbsubcsin = TMath::Sin(2.*(eventPlaneV0C - eventPlaneTPC));
1758   // three sub events in case of VZERO all
1759   if(fVZEROEventPlane && (!fVZEROEventPlaneA) && (!fVZEROEventPlaneC)) {
1760     if(!fSP){
1761       diffsubasubb = TMath::Cos(2.*(eventPlaneV0 - eventPlanesub1a));     //TR:
1762       diffsubasubc = TMath::Cos(2.*(eventPlaneV0 - eventPlanesub2a));     //TR:
1763       diffsubbsubc = TMath::Cos(2.*(eventPlanesub1a - eventPlanesub2a));  //TR:
1764     }
1765     else{
1766       diffsubasubb = qV0.X()*qsub1a->X()+qV0.Y()*qsub1a->Y();
1767       diffsubasubc = qV0.X()*qsub2a->X()+qV0.Y()*qsub2a->Y();
1768       diffsubbsubc = qsub1a->X()*qsub2a->X()+qsub1a->Y()*qsub2a->Y();
1769     }
1770
1771     diffsubasubbsin = TMath::Sin(2.*(eventPlaneV0 - eventPlanesub1a));
1772     diffsubasubcsin = TMath::Sin(2.*(eventPlaneV0 - eventPlanesub2a));
1773     diffsubbsubcsin = TMath::Sin(2.*(eventPlanesub1a - eventPlanesub2a));
1774   }
1775
1776   //////////////////////////////////////
1777   // AliFlowEvent  and MC event plane
1778   /////////////////////////////////////
1779
1780   Int_t nbtracks = fInputEvent->GetNumberOfTracks();
1781   AliDebug(2,Form("Number of tracks %d",nbtracks));
1782
1783   if(fMonitorQCumulant) {
1784
1785     fcutsRP->SetEvent( InputEvent(), MCEvent());
1786     fcutsPOI->SetEvent( InputEvent(), MCEvent());
1787     if( fflowEvent ){
1788       fflowEvent->~AliFlowEvent();
1789       new(fflowEvent) AliFlowEvent(fcutsRP,fcutsPOI);
1790     }else fflowEvent = new AliFlowEvent(fcutsRP,fcutsPOI);
1791     if(mcEvent && mcEvent->GenEventHeader()) {
1792       fflowEvent->SetMCReactionPlaneAngle(mcEvent);
1793       //if reaction plane not set from elsewhere randomize it before adding flow
1794       //if (!fflowEvent->IsSetMCReactionPlaneAngle()) fflowEvent->SetMCReactionPlaneAngle(gRandom->Uniform(0.0,TMath::TwoPi()));
1795       mcReactionPlane = TVector2::Phi_0_2pi(fflowEvent->GetMCReactionPlaneAngle());
1796       if(mcReactionPlane > TMath::Pi()) mcReactionPlane = mcReactionPlane - TMath::Pi();
1797       AliDebug(2,Form("MC reaction plane %f",mcReactionPlane));
1798     }
1799     fflowEvent->SetReferenceMultiplicity( nbtracks );
1800     fflowEvent->DefineDeadZone(0,0,0,0);
1801     //fflowEvent.TagSubeventsInEta(-0.8,-0.1,0.1,0.8);
1802
1803     ////////////////
1804     // MC
1805     ///////////////
1806     if(fUseMCReactionPlane) {
1807       eventPlanea = mcReactionPlane;
1808       diffsub1sub2a = 0.0;
1809     }
1810   }
1811
1812
1813   //////////////////////
1814   // Fill Histos
1815   //////////////////////
1816
1817   fHistEV->Fill(binctt,2.0);
1818
1819   // Fill
1820   valuensparsea[0] = eventPlaneV0A;
1821   valuensparsea[1] = eventPlaneV0C;
1822   valuensparsea[2] = eventPlaneTPC;
1823   if(fVZEROEventPlane && (!fVZEROEventPlaneA) && (!fVZEROEventPlaneC)) {
1824     // case VZERO all
1825     valuensparsea[0] = eventPlaneV0;
1826     valuensparsea[1] = eventPlanesub1a;
1827     valuensparsea[2] = eventPlanesub2a;
1828   }
1829   fEventPlane->Fill(&valuensparsea[0]);
1830
1831   // Fill
1832   if(fMonitorEventPlane) fCosSin2phiep->Fill(&valuecossinephiep[0]);
1833
1834   if(!fVZEROEventPlane) {
1835     valuensparsef[0] = diffsub1sub2a;
1836     fCosRes->Fill(&valuensparsef[0]);
1837     valuensparsefsin[0] = diffsub1sub2asin;
1838     if(fMonitorEventPlane) fSinRes->Fill(&valuensparsefsin[0]);
1839     if(fMonitorEventPlane) {
1840       fProfileCosRes->Fill(valuensparsef[1],valuensparsef[0]);
1841     }
1842   }
1843   else {
1844     valuensparsefbis[0] = diffsubasubb;
1845     valuensparsefbis[1] = diffsubasubc;
1846     valuensparsefbis[2] = diffsubbsubc;
1847     fCosResabc->Fill(&valuensparsefbis[0]); //TR: info
1848     valuensparsefbissin[0] = diffsubasubbsin;
1849     valuensparsefbissin[1] = diffsubbsubcsin;
1850     valuensparsefbissin[2] = diffsubasubcsin;
1851     if(fMonitorEventPlane) fSinResabc->Fill(&valuensparsefbissin[0]);
1852     if(fMonitorEventPlane) {
1853       fProfileCosResab->Fill(valuensparsefbis[3],valuensparsefbis[0]);
1854       fProfileCosResac->Fill(valuensparsefbis[3],valuensparsefbis[1]);
1855       fProfileCosResbc->Fill(valuensparsefbis[3],valuensparsefbis[2]);
1856     }
1857   }
1858
1859   ////////////////////////////////////////
1860   // Loop to determine pool background
1861   /////////////////////////////////////////
1862   if(fMonitorPhotonic) {
1863
1864     fBackgroundSubtraction->FillPoolAssociatedTracks(fInputEvent,binct);
1865
1866     if(  fArraytrack ){
1867       fArraytrack->~TArrayI();
1868       new(fArraytrack) TArrayI(nbtracks);
1869     }
1870     else {
1871       fArraytrack = new TArrayI(nbtracks);
1872     }
1873     fCounterPoolBackground = 0;
1874
1875     for(Int_t k = 0; k < nbtracks; k++){
1876
1877       AliVTrack *track = (AliVTrack *) fInputEvent->GetTrack(k);
1878       if(!track) continue;
1879
1880       // Track cuts
1881       Bool_t survivedbackground = kTRUE;
1882       if(fAODAnalysis) {
1883         AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
1884         if(aodtrack) {
1885           AliESDtrack esdTrack(aodtrack);
1886           // set the TPC cluster info
1887           esdTrack.SetTPCClusterMap(aodtrack->GetTPCClusterMap());
1888           esdTrack.SetTPCSharedMap(aodtrack->GetTPCSharedMap());
1889           esdTrack.SetTPCPointsF(aodtrack->GetTPCNclsF());
1890           // needed to calculate the impact parameters
1891           AliAODEvent *aodeventu = dynamic_cast<AliAODEvent *>(fInputEvent);
1892           if(aodeventu) {
1893             AliAODVertex *vAOD = aodeventu->GetPrimaryVertex();
1894             Double_t bfield = aodeventu->GetMagneticField();
1895             Double_t pos[3],cov[6];
1896             vAOD->GetXYZ(pos);
1897             vAOD->GetCovarianceMatrix(cov);
1898             const AliESDVertex vESD(pos,cov,100.,100);
1899             esdTrack.RelateToVertex(&vESD,bfield,3.);
1900           }
1901           if(!fHFEBackgroundCuts->IsSelected(&esdTrack)) {
1902             survivedbackground = kFALSE;
1903           }
1904         }
1905       }
1906       else {
1907         AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
1908         if(esdtrack) {
1909           if(!fHFEBackgroundCuts->IsSelected(esdtrack)) survivedbackground = kFALSE;
1910         }
1911       }
1912       // PID
1913       if(survivedbackground) {
1914         // PID track cuts
1915         AliHFEpidObject hfetrack2;
1916         if(!fAODAnalysis) hfetrack2.SetAnalysisType(AliHFEpidObject::kESDanalysis);
1917         else hfetrack2.SetAnalysisType(AliHFEpidObject::kAODanalysis);
1918         hfetrack2.SetRecTrack(track);
1919         hfetrack2.SetCentrality((Int_t)binct);
1920         AliDebug(2,Form("centrality %f and %d",binct,hfetrack2.GetCentrality()));
1921         hfetrack2.SetPbPb();
1922         if(fPIDBackground->IsSelected(&hfetrack2,0x0,"recTrackCont",fPIDBackgroundqa)) {
1923           fArraytrack->AddAt(k,fCounterPoolBackground);
1924           fCounterPoolBackground++;
1925           AliDebug(2,Form("fCounterPoolBackground %d, track %d",fCounterPoolBackground,k));
1926         }
1927       }
1928     }
1929   }
1930
1931
1932   //////////////////////////
1933   // Loop over track
1934   //////////////////////////
1935   for(Int_t k = 0; k < nbtracks; k++){
1936
1937     AliVTrack *track = (AliVTrack *) fInputEvent->GetTrack(k);
1938     if(!track) continue;
1939
1940     if(fAODAnalysis) {
1941       AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
1942       if(!aodtrack){
1943         AliError("AOD track is not there");
1944         continue;
1945       }
1946       AliDebug(2,"Find AOD track on");
1947       if(!(aodtrack->TestFilterBit(fFilter))) continue;  // Only process AOD tracks where the HFE is set
1948     }
1949
1950
1951     valuetrackingcuts[0] = track->Pt();
1952     valuetrackingcuts[1] = 0;
1953
1954     // RecKine: ITSTPC cuts
1955     if(!fHFECuts->CheckParticleCuts(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack, (TObject *)track)) continue;
1956     if(fMonitorTrackCuts) fTrackingCuts->Fill(&valuetrackingcuts[0]);
1957
1958     // RecPrim
1959     if(!fHFECuts->CheckParticleCuts(AliHFEcuts::kStepRecPrim + AliHFEcuts::kNcutStepsMCTrack, (TObject *)track)) continue;
1960     valuetrackingcuts[1] = 1;
1961     if(fMonitorTrackCuts) fTrackingCuts->Fill(&valuetrackingcuts[0]);
1962
1963     // HFEcuts: ITS layers cuts
1964     if(!fHFECuts->CheckParticleCuts(AliHFEcuts::kStepHFEcutsITS + AliHFEcuts::kNcutStepsMCTrack, (TObject *)track)) continue;
1965     valuetrackingcuts[1] = 2;
1966     if(fMonitorTrackCuts) fTrackingCuts->Fill(&valuetrackingcuts[0]);
1967
1968     // HFE cuts: TOF and mismatch flag
1969     if(!fHFECuts->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTOF + AliHFEcuts::kNcutStepsMCTrack, (TObject *)track)) continue;
1970     valuetrackingcuts[1] = 3;
1971     if(fMonitorTrackCuts) fTrackingCuts->Fill(&valuetrackingcuts[0]);
1972
1973     // HFE cuts: TPC PID cleanup
1974     if(!fHFECuts->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTPC + AliHFEcuts::kNcutStepsMCTrack, (TObject *)track)) continue;
1975     valuetrackingcuts[1] = 4;
1976     if(fMonitorTrackCuts) fTrackingCuts->Fill(&valuetrackingcuts[0]);
1977
1978     // HFEcuts: Nb of tracklets TRD0
1979     if(!fHFECuts->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack, (TObject *)track)) continue;
1980     valuetrackingcuts[1] = 5;
1981     if(fMonitorTrackCuts) fTrackingCuts->Fill(&valuetrackingcuts[0]);
1982
1983     AliDebug(2,"Survived");
1984
1985     /////////////////////////////////////////////////////////
1986     // Subtract candidate from TPC event plane
1987     ////////////////////////////////////////////////////////
1988     Float_t eventplanesubtracted = 0.0;
1989
1990     if(!fVZEROEventPlane) {
1991       // Subtract the tracks from the event plane
1992       Double_t qX = qTPC->X() - vEPa->GetQContributionX(track);  //Modify the components: subtract the track you want to look at with your analysis
1993       Double_t qY = qTPC->Y() - vEPa->GetQContributionY(track);  //Modify the components: subtract the track you want to look at with your analysis
1994       TVector2 newQVectorfortrack;
1995       newQVectorfortrack.Set(qX,qY);
1996       eventplanesubtracted = TVector2::Phi_0_2pi(newQVectorfortrack.Phi())/2;
1997     }
1998     else eventplanesubtracted = eventPlanea;
1999
2000     ///////////////////////////////////////////
2001     // Event plane
2002     //////////////////////////////////////////
2003     Bool_t fillEventPlane = kTRUE;
2004     if(!fVZEROEventPlane){
2005       if((!qsub1a) || (!qsub2a)) fillEventPlane = kFALSE;
2006       if(fSubEtaGapTPC) {
2007         if(track->Eta() < (- fEtaGap/2.)) eventplanesubtracted = eventPlanesub1a;
2008         else if(track->Eta() > (fEtaGap/2.)) eventplanesubtracted = eventPlanesub2a;
2009         else fillEventPlane = kFALSE;
2010       }
2011     }
2012
2013     ///////////////
2014     // MC
2015     //////////////
2016     if(fUseMCReactionPlane) {
2017       eventplanesubtracted = mcReactionPlane;
2018       fillEventPlane = kTRUE;
2019     }
2020
2021     //////////////////////////////////////////////////////////////////////////////
2022     ///////////////////////////AFTERBURNER
2023     Double_t phitrack = track->Phi();
2024     if (fAfterBurnerOn)
2025       {
2026         phitrack = GetPhiAfterAddV2(track->Phi(),mcReactionPlane);
2027       }
2028     //////////////////////////////////////////////////////////////////////////////
2029
2030
2031     ///////////////////////
2032     // Calculate deltaphi
2033     ///////////////////////
2034
2035     // Suppose phi track is between 0.0 and phi
2036     Double_t deltaphi = TVector2::Phi_0_2pi(phitrack - eventplanesubtracted);
2037     if(deltaphi > TMath::Pi()) deltaphi = deltaphi - TMath::Pi();
2038
2039     ////////////////////////////////
2040     // Determine the deltaphi bin
2041     ///////////////////////////////
2042
2043     // in-plane
2044     if(((deltaphi<(TMath::Pi()/4.)) && (deltaphi>0.0)) || ((deltaphi>(3*TMath::Pi()/4.)) && (deltaphi<TMath::Pi()))) valuedeltaphicontamination[0] = 0.5;
2045     // out-of-plane
2046     if((deltaphi>(TMath::Pi()/4.)) && (deltaphi<(3*TMath::Pi()/4.))) valuedeltaphicontamination[0] = 1.5;
2047
2048     ////////////////////////////////////////
2049     // Define variables
2050     ///////////////////////////////////////
2051
2052
2053     valuedeltaphicontamination[2] = track->Pt();
2054     valuensparsee[2] = track->Pt();
2055     valuensparsee[3] = track->Eta();
2056     valuensparseg[2] = track->Pt();
2057     valuensparseh[2] = track->Pt();
2058     valuefractioncont[0] = track->Pt();
2059     valuensparsehprofile[2] = track->Pt();
2060     valuensparseMCSourceDeltaPhiMaps[1] = track->Pt();
2061     if(track->Charge() > 0.0) {
2062       valuensparseg[3] = 0.2;
2063       valuensparseh[3] = 0.2;
2064     }
2065     else {
2066       valuensparseg[3] = -0.2;
2067       valuensparseh[3] = -0.2;
2068     }
2069     valuensparseh[4] = track->Eta();
2070     valuensparseg[4] = track->Eta();
2071
2072     AliDebug(2,Form("charge %d",(Int_t)track->Charge()));
2073
2074     ////////////////////////
2075     // Fill before PID
2076     ///////////////////////
2077
2078     if(fMonitorWithoutPID) {
2079
2080       valuensparseg[0] = deltaphi;
2081       if(fillEventPlane) fDeltaPhiMapsBeforePID->Fill(&valuensparseg[0]);
2082
2083       //
2084       valuensparseh[0] = TMath::Cos(2*TVector2::Phi_mpi_pi(phitrack-eventplanesubtracted));
2085       if(fillEventPlane) {
2086         fCosPhiMapsBeforePID->Fill(&valuensparseh[0]);
2087       }
2088     }
2089
2090     ////////////////////////
2091     // Apply PID
2092     ////////////////////////
2093     if(!fNoPID) {
2094       // Apply PID for Data
2095       if(!fMCPID) {
2096         // pid object
2097         AliHFEpidObject hfetrack;
2098         if(!fAODAnalysis) hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
2099         else hfetrack.SetAnalysisType(AliHFEpidObject::kAODanalysis);
2100         hfetrack.SetRecTrack(track);
2101         hfetrack.SetCentrality((Int_t)binct);
2102         AliDebug(2,Form("centrality %f and %d",binct,hfetrack.GetCentrality()));
2103         hfetrack.SetPbPb();
2104
2105         // Only TOF PID
2106         if(fMonitorContamination) {
2107           if(fPIDTOFOnly->IsSelected(&hfetrack,0x0,"recTrackCont",0x0)) {
2108             Float_t nsigma = pidResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
2109             valuedeltaphicontamination[3] = nsigma;
2110             fDeltaPhiMapsContamination->Fill(&valuedeltaphicontamination[0]);
2111           }
2112         }
2113
2114         // Complete PID TOF+TPC
2115         if(!fPID->IsSelected(&hfetrack,0x0,"recTrackCont",fPIDqa)) {
2116           continue;
2117         }
2118       }
2119       else {
2120         if(!mcEvent) continue;
2121         if(!(mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel()))))) continue;
2122         AliDebug(2,Form("PdgCode %d",TMath::Abs(mctrack->Particle()->GetPdgCode())));
2123         if(TMath::Abs(mctrack->Particle()->GetPdgCode())!=11) continue;
2124       }
2125     }
2126
2127
2128     /////////////////////////////////////////////////////////////////////////////
2129     // Add candidate to AliFlowEvent for POI and subtract from RP if needed
2130     ////////////////////////////////////////////////////////////////////////////
2131     if(fMonitorQCumulant) {
2132       Int_t idtrack = static_cast<AliVTrack*>(track)->GetID();
2133       Bool_t found = kFALSE;
2134       Int_t numberoffound = 0;
2135       AliDebug(2,Form("A: Number of tracks %d",fflowEvent->NumberOfTracks()));
2136       for(Int_t iRPs=0; iRPs< fflowEvent->NumberOfTracks(); iRPs++) {
2137         AliFlowTrack *iRP = (AliFlowTrack*) (fflowEvent->GetTrack(iRPs));
2138         if(!iRP) continue;
2139         //if(!iRP->InRPSelection()) continue;
2140         if( TMath::Abs(idtrack) == TMath::Abs(iRP->GetID()) ) {
2141           iRP->SetForPOISelection(kTRUE);
2142           found = kTRUE;
2143           numberoffound ++;
2144         }
2145       }
2146       AliDebug(2,Form("Found %d mal",numberoffound));
2147       if(!found) {
2148         AliFlowCandidateTrack *sTrack = (AliFlowCandidateTrack*) MakeTrack(massElectron,track->Pt(),track->Phi(), track->Eta());
2149         sTrack->SetID(idtrack);
2150         fflowEvent->AddTrack(sTrack);
2151         AliDebug(2,"Add the track");
2152       }
2153       AliDebug(2,Form("B: Number of tracks %d",fflowEvent->NumberOfTracks()));
2154     }
2155
2156
2157     /////////////////////
2158     // Fill THnSparseF
2159     /////////////////////
2160
2161     //
2162     valuensparseabis[0] = eventplanesubtracted;
2163     if((fillEventPlane) && (fMonitorEventPlane)) fEventPlaneaftersubtraction->Fill(&valuensparseabis[0]);
2164
2165
2166     if(fMonitorEventPlane)
2167       {
2168         //
2169         valuensparsee[0] = TMath::Cos(2*phitrack);
2170         fCos2phie->Fill(&valuensparsee[0]);
2171         valuensparsee[0] = TMath::Sin(2*phitrack);
2172         fSin2phie->Fill(&valuensparsee[0]);
2173         //
2174         valuensparsee[0] = TMath::Cos(2*eventplanesubtracted);
2175         if(fillEventPlane) fCos2phiep->Fill(&valuensparsee[0]);
2176         valuensparsee[0] = TMath::Sin(2*eventplanesubtracted);
2177         if(fillEventPlane) fSin2phiep->Fill(&valuensparsee[0]);
2178         valuensparsee[0] = TMath::Sin(2*TVector2::Phi_mpi_pi(phitrack-eventplanesubtracted));
2179         if(fillEventPlane) fSin2phiephiep->Fill(&valuensparsee[0]);
2180         //
2181       }
2182
2183     //
2184     valuensparseg[0] = deltaphi;
2185     if(fillEventPlane) fDeltaPhiMaps->Fill(&valuensparseg[0]);
2186
2187     //
2188     valuensparseh[0] = TMath::Cos(2*TVector2::Phi_mpi_pi(phitrack-eventplanesubtracted));
2189     if(fillEventPlane) {
2190       fCosPhiMaps->Fill(&valuensparseh[0]); //TR: fCosPhiQSum+=valuensparseh[0]*TMath:Sqrt(qAna->X()*qAna->X()+qAna->Y()*qAna->Y()); fCosPhiQN++;
2191       if((valuefractioncont[1] >=0) && (valuefractioncont[1] < 11)){
2192         if(fContamination[((Int_t)valuefractioncont[1])]){
2193           Double_t weight = 1.;
2194           if(fAsFunctionOfP) weight = fContamination[((Int_t)valuefractioncont[1])]->Eval(track->P());
2195           else weight = fContamination[((Int_t)valuefractioncont[1])]->Eval(track->Pt());
2196           if(weight<0.0) weight=0.0;
2197           if(weight>1.0) weight=1.0;
2198           fFractionContamination->Fill(&valuefractioncont[0],weight);
2199           if(fv2contamination[((Int_t)valuefractioncont[1])]){
2200             Double_t v2 =  fv2contamination[((Int_t)valuefractioncont[1])]->Eval(track->Pt());
2201             AliDebug(2,Form("value v2 %f, contamination %f and pt %f centrality %d\n",v2,weight,track->Pt(),(Int_t)valuefractioncont[1]));
2202             AliDebug(2,Form("Check for centrality 3: value v2 %f, contamination %f\n",fv2contamination[3]->Eval(track->Pt()),fContamination[3]->Eval(track->P())));
2203             AliDebug(2,Form("Check for centrality 4: value v2 %f, contamination %f\n",fv2contamination[4]->Eval(track->Pt()),fContamination[4]->Eval(track->P())));
2204             AliDebug(2,Form("Check for centrality 5: value v2 %f, contamination %f\n",fv2contamination[5]->Eval(track->Pt()),fContamination[5]->Eval(track->P())));
2205             fContaminationv2->Fill(valuefractioncont[1],valuefractioncont[0],v2,weight);
2206           }
2207           fContaminationmeanpt->Fill(valuefractioncont[1],valuefractioncont[0],TMath::Abs(track->Pt()));
2208         }
2209       }
2210       if(fMonitorEventPlane) {
2211         if(fSP)
2212           valuensparseh[0] *= TMath::Sqrt(qAna->X()*qAna->X()+qAna->Y()*qAna->Y());
2213         fProfileCosPhiMaps->Fill(valuensparsehprofile[1],valuensparsehprofile[2],valuensparseh[0]);  //TR: info
2214       }
2215     }
2216
2217     if(fMonitorPhotonic) {
2218       Int_t indexmother = -1;
2219       Int_t source = 1;
2220       if(mcthere) source = fBackgroundSubtraction->FindMother(mctrack->GetLabel(),indexmother);
2221       fBackgroundSubtraction->LookAtNonHFE(k, track, fInputEvent, 1, binct, deltaphi, source, indexmother);
2222
2223       if((!fAODAnalysis && mcthere) || !mcthere) {
2224         // background
2225         source = 0;
2226         indexmother = -1;
2227         source = FindMother(TMath::Abs(track->GetLabel()),mcEvent, indexmother);
2228         valuensparseMCSourceDeltaPhiMaps[2] = source;
2229         if(mcEvent) fMCSourceDeltaPhiMaps->Fill(&valuensparseMCSourceDeltaPhiMaps[0]);
2230         //LookAtNonHFE(k,track,fInputEvent,mcEvent,binct,deltaphi,source,indexmother);
2231         Int_t taggedvalue = LookAtNonHFE(k,track,fInputEvent,mcEvent,binct,deltaphi,source,indexmother);
2232         if(fMonitorPhotonic) {
2233           // No opposite charge partner found in the invariant mass choosen
2234           if((taggedvalue!=2) && (taggedvalue!=6)) {
2235             //fDeltaPhiMapsTaggedNonPhotonic->Fill(&valuensparseg[0]);
2236             //fCosPhiMapsTaggedNonPhotonic->Fill(&valuensparseh[0]);
2237           }
2238           // One opposite charge partner found in the invariant mass choosen
2239           if((taggedvalue==2) || (taggedvalue==6)) {
2240             fDeltaPhiMapsTaggedPhotonic->Fill(&valuensparseg[0]);
2241             //fCosPhiMapsTaggedPhotonic->Fill(&valuensparseh[0]);
2242           }
2243           // One same charge partner found in the invariant mass choosen
2244           if((taggedvalue==4) || (taggedvalue==6)) {
2245             fDeltaPhiMapsTaggedPhotonicLS->Fill(&valuensparseg[0]);
2246             //fCosPhiMapsTaggedPhotonicLS->Fill(&valuensparseh[0]);
2247           }
2248         }
2249       }
2250     }
2251
2252   }
2253
2254   //////////////////////////////////////////////////////////////////////////////
2255   ///////////////////////////AFTERBURNER
2256   if (fAfterBurnerOn &  fMonitorQCumulant)
2257     {
2258       fflowEvent->AddFlow(fV1,fV2,fV3,fV4,fV5);     //add flow
2259       fflowEvent->CloneTracks(fNonFlowNumberOfTrackClones); //add nonflow by cloning tracks
2260     }
2261   //////////////////////////////////////////////////////////////////////////////
2262
2263
2264
2265   //for(Int_t bincless = 0; bincless < fNbBinsCentralityQCumulant; bincless++) {
2266   //  if((fBinCentralityLess[bincless]< cntr) && (cntr < fBinCentralityLess[bincless+1])) PostData(bincless+2,fflowEvent);
2267   //}
2268
2269   if(fMonitorPhotonic) {
2270     if(fArraytrack) {
2271       delete fArraytrack;
2272       fArraytrack = NULL;
2273     }
2274   }
2275
2276   if(fMonitorPhotonic) fBackgroundSubtraction->CountPoolAssociated(fInputEvent,binct);
2277
2278   PostData(1, fListHist);
2279
2280
2281
2282 }
2283 //______________________________________________________________________________
2284 AliFlowCandidateTrack *AliAnalysisTaskFlowTPCTOFEPSP::MakeTrack( Double_t mass,
2285                           Double_t pt, Double_t phi, Double_t eta) {
2286   //
2287   //  Make Track (Not needed actually)
2288   //
2289
2290   AliFlowCandidateTrack *sTrack = new AliFlowCandidateTrack();
2291   sTrack->SetMass(mass);
2292   sTrack->SetPt(pt);
2293   sTrack->SetPhi(phi);
2294   sTrack->SetEta(eta);
2295   sTrack->SetForPOISelection(kTRUE);
2296   sTrack->SetForRPSelection(kFALSE);
2297   return sTrack;
2298 }
2299 //_________________________________________________________________________________
2300 Double_t AliAnalysisTaskFlowTPCTOFEPSP::GetPhiAfterAddV2(Double_t phi,Double_t reactionPlaneAngle) const
2301 {
2302   //
2303   // Adds v2, uses Newton-Raphson iteration
2304   //
2305   Double_t phiend=phi;
2306   Double_t phi0=phi;
2307   Double_t f=0.;
2308   Double_t fp=0.;
2309   Double_t phiprev=0.;
2310
2311   for (Int_t i=0; i<fMaxNumberOfIterations; i++)
2312   {
2313     phiprev=phiend; //store last value for comparison
2314     f =  phiend-phi0+fV2*TMath::Sin(2.*(phiend-reactionPlaneAngle));
2315     fp = 1.0+2.0*fV2*TMath::Cos(2.*(phiend-reactionPlaneAngle)); //first derivative
2316     phiend -= f/fp;
2317     if (TMath::AreEqualAbs(phiprev,phiend,fPrecisionPhi)) break;
2318   }
2319   return phiend;
2320 }
2321 //_____________________________________________________________________________________________
2322 Int_t AliAnalysisTaskFlowTPCTOFEPSP::LookAtNonHFE(Int_t iTrack1, AliVTrack *track1, AliVEvent *vEvent, AliMCEvent *mcEvent,Int_t binct,Double_t deltaphi,Int_t source,Int_t indexmother)
2323 {
2324   //
2325   // Look At Non HFE
2326   //
2327
2328   // return -1 if nothing
2329   // return 2 if opposite charge within the mass range found
2330   // return 4 if like charge within the mass range found
2331   // return 6 if opposite charge and like charge within the mass range found
2332   //
2333
2334   Int_t taggedphotonic = -1;
2335
2336   Bool_t oppositetaggedphotonic = kFALSE;
2337   Bool_t sametaggedphotonic = kFALSE;
2338
2339   AliDebug(2,Form("fCounterPoolBackground %d in LookAtNonHFE!!!",fCounterPoolBackground));
2340   if(!fArraytrack) return taggedphotonic;
2341   AliDebug(2,Form("process track %d",iTrack1));
2342
2343   TVector3 v3D1;
2344   TVector3 v3D2;
2345
2346   Double_t valuensparseDeltaPhiMaps[5];
2347   Double_t valueangle[3];
2348
2349   valuensparseDeltaPhiMaps[1] = binct;
2350   valuensparseDeltaPhiMaps[2] = track1->Pt();
2351   valuensparseDeltaPhiMaps[0] = deltaphi;
2352   valuensparseDeltaPhiMaps[4] = source;
2353
2354   valueangle[2] = source;
2355   valueangle[1] = binct;
2356
2357   // Pdg code
2358   Int_t pdg1 = CheckPdg(TMath::Abs(track1->GetLabel()),mcEvent);
2359   Int_t numberfound = 0;
2360
2361   //Magnetic Field
2362   Double_t bfield = vEvent->GetMagneticField();
2363
2364   // Get Primary vertex
2365   const AliVVertex *pVtx = vEvent->GetPrimaryVertex();
2366
2367   for(Int_t idex = 0; idex < fCounterPoolBackground; idex++)
2368     {
2369
2370       Int_t iTrack2 = fArraytrack->At(idex);
2371       AliDebug(2,Form("track %d",iTrack2));
2372       AliVTrack* track2 = (AliVTrack *) vEvent->GetTrack(iTrack2);
2373       if (!track2)
2374         {
2375           printf("ERROR: Could not receive track %d", iTrack2);
2376           continue;
2377         }
2378       if(iTrack2==iTrack1) continue;
2379       AliDebug(2,"Different");
2380
2381       // Reset the MC info
2382       valueangle[2] = source;
2383       valuensparseDeltaPhiMaps[4] = source;
2384
2385       // track cuts and PID already done
2386
2387       // if MC look
2388       Int_t pdg2 = -100;
2389       if(mcEvent) {
2390         Int_t source2 = 0;
2391         Int_t indexmother2 = -1;
2392         source2 = FindMother(TMath::Abs(track2->GetLabel()),mcEvent, indexmother2);
2393         pdg2 = CheckPdg(TMath::Abs(track2->GetLabel()),mcEvent);
2394         if(source2 >=0 ) {
2395           if((indexmother2 == indexmother) && (source == source2) && ((pdg1*pdg2)<0.0)) {
2396             if(source == kElectronfromconversion) {
2397               valueangle[2] = kElectronfromconversionboth;
2398               valuensparseDeltaPhiMaps[4] = kElectronfromconversionboth;
2399               numberfound++;
2400             }
2401             if(source == kElectronfrompi0) {
2402               valueangle[2] = kElectronfrompi0both;
2403               valuensparseDeltaPhiMaps[4] = kElectronfrompi0both;
2404             }
2405             if(source == kElectronfrometa) {
2406               valueangle[2] = kElectronfrometaboth;
2407               valuensparseDeltaPhiMaps[4] = kElectronfrometaboth;
2408             }
2409           }
2410         }
2411       }
2412
2413       if(fAlgorithmMA && (!fAODAnalysis))
2414         {
2415           // tracks
2416           AliESDtrack *esdtrack2 = dynamic_cast<AliESDtrack *>(track2);
2417           AliESDtrack *esdtrack1 = dynamic_cast<AliESDtrack *>(track1);
2418           if((!esdtrack2) || (!esdtrack1)) continue;
2419
2420           //Variables
2421           Double_t p1[3];
2422           Double_t p2[3];
2423           Double_t xt1; //radial position track 1 at the DCA point
2424           Double_t xt2; //radial position track 2 at the DCA point
2425           //DCA track1-track2
2426           Double_t dca12 = esdtrack2->GetDCA(esdtrack1,bfield,xt2,xt1);
2427
2428           // Cut dca
2429           if(dca12 > fMaxdca) continue;
2430
2431           //Momento of the track extrapolated to DCA track-track
2432           //Track1
2433           Bool_t hasdcaT1 = esdtrack1->GetPxPyPzAt(xt1,bfield,p1);
2434           //Track2
2435           Bool_t hasdcaT2 = esdtrack2->GetPxPyPzAt(xt2,bfield,p2);
2436
2437           if(!hasdcaT1 || !hasdcaT2) AliWarning("It could be a problem in the extrapolation");
2438
2439           //track1-track2 Invariant Mass
2440           Double_t eMass = 0.000510998910; //Electron mass in GeV
2441           Double_t pP1 = sqrt(p1[0]*p1[0]+p1[1]*p1[1]+p1[2]*p1[2]); //Track 1 momentum
2442           Double_t pP2 = sqrt(p2[0]*p2[0]+p2[1]*p2[1]+p2[2]*p2[2]); //Track 2 momentum
2443           Double_t eE1 = TMath::Sqrt(pP1*pP1+eMass*eMass);
2444           Double_t eE2 = TMath::Sqrt(pP2*pP2+eMass*eMass);
2445
2446           //TLorentzVector v1(p1[0],p1[1],p1[2],sqrt(eMass*eMass+pP1*pP1));
2447           //TLorentzVector v2(p2[0],p2[1],p2[2],sqrt(eMass*eMass+pP2*pP2));
2448           //Double_t imass = (v1+v2).M(); //Invariant Mass
2449           //Double_t angle3D = v1.Angle(v2.Vect()); //Opening Angle (Total Angle)
2450
2451           // daughter
2452           v3D1.SetXYZ(p1[0],p1[1],p1[2]);
2453           v3D2.SetXYZ(p2[0],p2[1],p2[2]);
2454           Double_t openingangle = TVector2::Phi_0_2pi(v3D2.Angle(v3D1));
2455
2456           // mother
2457           TVector3 motherrec = v3D1 + v3D2;
2458           Double_t invmass = TMath::Sqrt((eE1+eE2)*(eE1+eE2)-(motherrec.Px()*motherrec.Px()+motherrec.Py()*motherrec.Py()+motherrec.Pz()*motherrec.Pz()));
2459
2460           // xy
2461           //TVector3 vectordiff = v3D1 - v3D2;
2462           //Double_t diffphi = TVector2::Phi_0_2pi(vectordiff.Phi());
2463           //Double_t massxy = TMath::Sqrt((eE1+eE2)*(eE1+eE2)-(pP1*pP1+pP2*pP2+2*pP1*pP2*TMath::Cos(diffphi)));
2464
2465           // rz
2466           //Double_t difftheta = TVector2::Phi_0_2pi(vectordiff.Eta());
2467           //Double_t massrz = TMath::Sqrt((eE1+eE2)*(eE1+eE2)-(pP1*pP1+pP2*pP2+2*pP1*pP2*TMath::Cos(difftheta)));
2468
2469
2470           Float_t fCharge1 = track1->Charge();
2471           Float_t fCharge2 = track2->Charge();
2472
2473           // Fill Histo
2474           //valueangle[0] = diffphi;
2475           //valueangle[1] = difftheta;
2476           valueangle[0] = openingangle;
2477           if((fCharge1*fCharge2)>0.0) fSameSignAngle->Fill(&valueangle[0]);
2478           else fOppSignAngle->Fill(&valueangle[0]);
2479
2480           // Cut
2481           if(openingangle > fMaxopening3D) continue;
2482           //if(difftheta > fMaxopeningtheta) continue;
2483           //if(diffphi > fMaxopeningphi) continue;
2484
2485           // Invmass
2486           valuensparseDeltaPhiMaps[3] = invmass;
2487           if((fCharge1*fCharge2)>0.0) fSameSignDeltaPhiMaps->Fill(&valuensparseDeltaPhiMaps[0]);
2488           else fOppSignDeltaPhiMaps->Fill(&valuensparseDeltaPhiMaps[0]);
2489
2490           // Cut
2491           if(invmass < fMaxInvmass) {
2492             if((fCharge1*fCharge2)<0.0) oppositetaggedphotonic=kTRUE;
2493             if((fCharge1*fCharge2)>0.0) sametaggedphotonic=kTRUE;
2494           }
2495
2496
2497         }
2498       else
2499         {
2500           Int_t fPDGtrack1 = 11;
2501           Int_t fPDGtrack2 = 11;
2502
2503           Float_t fCharge1 = track1->Charge();
2504           Float_t fCharge2 = track2->Charge();
2505
2506           if(fCharge1>0) fPDGtrack1 = -11;
2507           if(fCharge2>0) fPDGtrack2 = -11;
2508
2509           AliKFParticle ktrack1(*track1, fPDGtrack1);
2510           AliKFParticle ktrack2(*track2, fPDGtrack2);
2511           AliKFParticle recoGamma(ktrack1, ktrack2);
2512
2513           //Reconstruction Cuts
2514           if(recoGamma.GetNDF()<1) continue;
2515           Double_t chi2OverNDF = recoGamma.GetChi2()/recoGamma.GetNDF();
2516           if(TMath::Sqrt(TMath::Abs(chi2OverNDF))>fChi2OverNDFCut) continue;
2517
2518           // DCA
2519           //Double_t dca12 = ktrack1.GetDistanceFromParticle(ktrack2);
2520           //if(dca12 > fMaxdca) continue;
2521
2522           // if set mass constraint
2523           if(fSetMassConstraint && pVtx) {
2524             AliKFVertex primV(*pVtx);
2525             primV += recoGamma;
2526             primV -= ktrack1;
2527             primV -= ktrack2;
2528             recoGamma.SetProductionVertex(primV);
2529             recoGamma.SetMassConstraint(0,0.0001);
2530           }
2531
2532           //Invariant Mass
2533           Double_t imass;
2534           Double_t width;
2535           recoGamma.GetMass(imass,width);
2536
2537           //Opening Angle (Total Angle)
2538           Double_t angle = ktrack1.GetAngle(ktrack2);
2539           valueangle[0] = angle;
2540           if((fCharge1*fCharge2)>0.0) fSameSignAngle->Fill(&valueangle[0]);
2541           else fOppSignAngle->Fill(&valueangle[0]);
2542
2543           // Cut
2544           if(angle > fMaxopening3D) continue;
2545
2546           // Invmass
2547           valuensparseDeltaPhiMaps[3] = imass;
2548           if((fCharge1*fCharge2)>0.0) fSameSignDeltaPhiMaps->Fill(&valuensparseDeltaPhiMaps[0]);
2549           else {
2550             fOppSignDeltaPhiMaps->Fill(&valuensparseDeltaPhiMaps[0]);
2551             /*
2552             if(valueangle[2] == kElectronfromconversionboth) {
2553               printf("Reconstructed charge1 %f, charge 2 %f and invmass %f",fCharge1,fCharge2,imass);
2554               printf("MC charge1 %d, charge 2 %d",pdg1,pdg2);
2555               printf("DCA %f",dca12);
2556               printf("Number of found %d",numberfound);
2557             }
2558             */
2559           }
2560
2561           // Cut
2562           if(imass < fMaxInvmass) {
2563             if((fCharge1*fCharge2)<0.0) oppositetaggedphotonic=kTRUE;
2564             if((fCharge1*fCharge2)>0.0) sametaggedphotonic=kTRUE;
2565           }
2566         }
2567     }
2568
2569   if(oppositetaggedphotonic && sametaggedphotonic){
2570     taggedphotonic = 6;
2571   }
2572
2573   if(!oppositetaggedphotonic && sametaggedphotonic){
2574     taggedphotonic = 4;
2575   }
2576
2577   if(oppositetaggedphotonic && !sametaggedphotonic){
2578     taggedphotonic = 2;
2579   }
2580
2581
2582   return taggedphotonic;
2583 }
2584 //_________________________________________________________________________
2585 Int_t AliAnalysisTaskFlowTPCTOFEPSP::FindMother(Int_t tr, AliMCEvent *mcEvent, Int_t &indexmother){
2586   //
2587   // Find the mother if MC
2588   //
2589
2590   if(!mcEvent) return 0;
2591
2592   Int_t pdg = CheckPdg(tr,mcEvent);
2593   if(TMath::Abs(pdg)!= 11) {
2594     indexmother = -1;
2595     return kNoElectron;
2596   }
2597
2598   indexmother = IsMotherGamma(tr,mcEvent);
2599   if(indexmother > 0) return kElectronfromconversion;
2600   indexmother = IsMotherPi0(tr,mcEvent);
2601   if(indexmother > 0) return kElectronfrompi0;
2602   indexmother = IsMotherC(tr,mcEvent);
2603   if(indexmother > 0) return kElectronfromC;
2604   indexmother = IsMotherB(tr,mcEvent);
2605   if(indexmother > 0) return kElectronfromB;
2606   indexmother = IsMotherEta(tr,mcEvent);
2607   if(indexmother > 0) return kElectronfrometa;
2608
2609   return kElectronfromother;
2610
2611
2612 }
2613 //____________________________________________________________________________________________________________
2614 Int_t AliAnalysisTaskFlowTPCTOFEPSP::CheckPdg(Int_t tr, AliMCEvent* mcEvent) {
2615
2616   //
2617   // Return the pdg of the particle
2618   //
2619
2620
2621   Int_t pdgcode = -1;
2622   if(tr < 0) return pdgcode;
2623
2624   if(!mcEvent) return pdgcode;
2625
2626   AliVParticle *mctrack = mcEvent->GetTrack(tr);
2627
2628
2629   if(mctrack->IsA() == AliMCParticle::Class()) {
2630     AliMCParticle *mctrackesd = NULL;
2631     if(!(mctrackesd = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return pdgcode;
2632     pdgcode = mctrackesd->PdgCode();
2633   }
2634
2635   if(mctrack->IsA() == AliAODMCParticle::Class()) {
2636     AliAODMCParticle *mctrackaod = NULL;
2637     if(!(mctrackaod = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return pdgcode;
2638     pdgcode = mctrackaod->GetPdgCode();
2639   }
2640
2641   return pdgcode;
2642
2643
2644 }
2645 //____________________________________________________________________________________________________________
2646 Int_t AliAnalysisTaskFlowTPCTOFEPSP::IsMotherGamma(Int_t tr, AliMCEvent* mcEvent) {
2647
2648   //
2649   // Return the lab of gamma mother or -1 if not gamma
2650   //
2651
2652   if(tr < 0) return -1;
2653   AliVParticle *mctrack = mcEvent->GetTrack(tr);
2654
2655   if(mctrack->IsA() == AliMCParticle::Class()) {
2656     AliMCParticle *mctrackesd = NULL;
2657     if(!(mctrackesd = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2658     TParticle *particle = 0x0;
2659     particle = mctrackesd->Particle();
2660     // Take mother
2661     if(!particle) return -1;
2662     Int_t imother   = particle->GetFirstMother();
2663     if(imother < 0) return -1;
2664     AliMCParticle *mothertrack = NULL;
2665     if(!(mothertrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2666     TParticle * mother = mothertrack->Particle();
2667     if(!mother) return -1;
2668     // Check gamma
2669     Int_t pdg = mother->GetPdgCode();
2670     if(TMath::Abs(pdg) == 22) return imother;
2671     if(TMath::Abs(pdg) == 11) {
2672       return IsMotherGamma(imother,mcEvent);
2673     }
2674     return -1;
2675   }
2676
2677   if(mctrack->IsA() == AliAODMCParticle::Class()) {
2678     AliAODMCParticle *mctrackaod = NULL;
2679     if(!(mctrackaod = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2680     // Take mother
2681     Int_t imother = mctrackaod->GetMother();
2682     if(imother < 0) return -1;
2683     AliAODMCParticle *mothertrack = NULL;
2684     if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2685     // Check gamma
2686     Int_t pdg = mothertrack->GetPdgCode();
2687     if(TMath::Abs(pdg) == 22) return imother;
2688     if(TMath::Abs(pdg) == 11) {
2689       return IsMotherGamma(imother,mcEvent);
2690     }
2691     return -1;
2692
2693   }
2694
2695   return -1;
2696
2697
2698 }
2699 //
2700 //____________________________________________________________________________________________________________
2701 Int_t AliAnalysisTaskFlowTPCTOFEPSP::IsMotherPi0(Int_t tr, AliMCEvent* mcEvent) {
2702
2703   //
2704   // Return the lab of pi0 mother or -1 if not pi0
2705   //
2706
2707   if(tr < 0) return -1;
2708   AliVParticle *mctrack = mcEvent->GetTrack(tr);
2709
2710   if(mctrack->IsA() == AliMCParticle::Class()) {
2711     AliMCParticle *mctrackesd = NULL;
2712     if(!(mctrackesd = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2713     TParticle *particle = 0x0;
2714     particle = mctrackesd->Particle();
2715     // Take mother
2716     if(!particle) return -1;
2717     Int_t imother   = particle->GetFirstMother();
2718     if(imother < 0) return -1;
2719     AliMCParticle *mothertrack = NULL;
2720     if(!(mothertrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2721     TParticle * mother = mothertrack->Particle();
2722     if(!mother) return -1;
2723     // Check gamma
2724     Int_t pdg = mother->GetPdgCode();
2725     if(TMath::Abs(pdg) == 111) return imother;
2726     if(TMath::Abs(pdg) == 11) {
2727       return IsMotherPi0(imother,mcEvent);
2728     }
2729     return -1;
2730   }
2731
2732   if(mctrack->IsA() == AliAODMCParticle::Class()) {
2733     AliAODMCParticle *mctrackaod = NULL;
2734     if(!(mctrackaod = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2735     // Take mother
2736     Int_t imother = mctrackaod->GetMother();
2737     if(imother < 0) return -1;
2738     AliAODMCParticle *mothertrack = NULL;
2739     if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2740     // Check gamma
2741     Int_t pdg = mothertrack->GetPdgCode();
2742     if(TMath::Abs(pdg) == 111) return imother;
2743     if(TMath::Abs(pdg) == 11) {
2744       return IsMotherPi0(imother,mcEvent);
2745     }
2746     return -1;
2747   }
2748
2749   return -1;
2750
2751 }
2752 //____________________________________________________________________________________________________________
2753 Int_t AliAnalysisTaskFlowTPCTOFEPSP::IsMotherC(Int_t tr, AliMCEvent* mcEvent) {
2754
2755   //
2756   // Return the lab of signal mother or -1 if not signal
2757   //
2758
2759   if(tr < 0) return -1;
2760   AliVParticle *mctrack = mcEvent->GetTrack(tr);
2761
2762   if(mctrack->IsA() == AliMCParticle::Class()) {
2763     AliMCParticle *mctrackesd = NULL;
2764     if(!(mctrackesd = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2765     TParticle *particle = 0x0;
2766     particle = mctrackesd->Particle();
2767     // Take mother
2768     if(!particle) return -1;
2769     Int_t imother   = particle->GetFirstMother();
2770     if(imother < 0) return -1;
2771     AliMCParticle *mothertrack = NULL;
2772     if(!(mothertrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2773     TParticle * mother = mothertrack->Particle();
2774     if(!mother) return -1;
2775     // Check gamma
2776     Int_t pdg = mother->GetPdgCode();
2777     if((TMath::Abs(pdg)==411) || (TMath::Abs(pdg)==421) || (TMath::Abs(pdg)==431) || (TMath::Abs(pdg)==4122) || (TMath::Abs(pdg)==4132) || (TMath::Abs(pdg)==4232) || (TMath::Abs(pdg)==43320)) return imother;
2778     if(TMath::Abs(pdg) == 11) {
2779       return IsMotherC(imother,mcEvent);
2780     }
2781     return -1;
2782   }
2783
2784   if(mctrack->IsA() == AliAODMCParticle::Class()) {
2785     AliAODMCParticle *mctrackaod = NULL;
2786     if(!(mctrackaod = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2787     // Take mother
2788     Int_t imother = mctrackaod->GetMother();
2789     if(imother < 0) return -1;
2790     AliAODMCParticle *mothertrack = NULL;
2791     if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2792     // Check gamma
2793     Int_t pdg = mothertrack->GetPdgCode();
2794     if((TMath::Abs(pdg)==411) || (TMath::Abs(pdg)==421) || (TMath::Abs(pdg)==431) || (TMath::Abs(pdg)==4122) || (TMath::Abs(pdg)==4132) || (TMath::Abs(pdg)==4232) || (TMath::Abs(pdg)==43320)) return imother;
2795     if(TMath::Abs(pdg) == 11) {
2796       return IsMotherC(imother,mcEvent);
2797     }
2798     return -1;
2799   }
2800
2801   return -1;
2802
2803 }
2804 //____________________________________________________________________________________________________________
2805 Int_t AliAnalysisTaskFlowTPCTOFEPSP::IsMotherB(Int_t tr, AliMCEvent* mcEvent) {
2806
2807   //
2808   // Return the lab of signal mother or -1 if not signal
2809   //
2810
2811   if(tr < 0) return -1;
2812   AliVParticle *mctrack = mcEvent->GetTrack(tr);
2813
2814   if(mctrack->IsA() == AliMCParticle::Class()) {
2815     AliMCParticle *mctrackesd = NULL;
2816     if(!(mctrackesd = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2817     TParticle *particle = 0x0;
2818     particle = mctrackesd->Particle();
2819     // Take mother
2820     if(!particle) return -1;
2821     Int_t imother   = particle->GetFirstMother();
2822     if(imother < 0) return -1;
2823     AliMCParticle *mothertrack = NULL;
2824     if(!(mothertrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2825     TParticle * mother = mothertrack->Particle();
2826     if(!mother) return -1;
2827     // Check gamma
2828     Int_t pdg = mother->GetPdgCode();
2829     if((TMath::Abs(pdg)==511) || (TMath::Abs(pdg)==521) || (TMath::Abs(pdg)==531) || (TMath::Abs(pdg)==5122) || (TMath::Abs(pdg)==5132) || (TMath::Abs(pdg)==5232) || (TMath::Abs(pdg)==53320)) return imother;
2830     if(TMath::Abs(pdg) == 11) {
2831       return IsMotherB(imother,mcEvent);
2832     }
2833     return -1;
2834   }
2835
2836   if(mctrack->IsA() == AliAODMCParticle::Class()) {
2837     AliAODMCParticle *mctrackaod = NULL;
2838     if(!(mctrackaod = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2839     // Take mother
2840     Int_t imother = mctrackaod->GetMother();
2841     if(imother < 0) return -1;
2842     AliAODMCParticle *mothertrack = NULL;
2843     if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2844     // Check gamma
2845     Int_t pdg = mothertrack->GetPdgCode();
2846     if((TMath::Abs(pdg)==511) || (TMath::Abs(pdg)==521) || (TMath::Abs(pdg)==531) || (TMath::Abs(pdg)==5122) || (TMath::Abs(pdg)==5132) || (TMath::Abs(pdg)==5232) || (TMath::Abs(pdg)==53320)) return imother;
2847     if(TMath::Abs(pdg) == 11) {
2848       return IsMotherB(imother,mcEvent);
2849     }
2850     return -1;
2851   }
2852
2853   return -1;
2854
2855 }
2856 //____________________________________________________________________________________________________________
2857 Int_t AliAnalysisTaskFlowTPCTOFEPSP::IsMotherEta(Int_t tr, AliMCEvent* mcEvent) {
2858
2859   //
2860   // Return the lab of pi0 mother or -1 if not pi0
2861   //
2862
2863  if(tr < 0) return -1;
2864   AliVParticle *mctrack = mcEvent->GetTrack(tr);
2865
2866   if(mctrack->IsA() == AliMCParticle::Class()) {
2867     AliMCParticle *mctrackesd = NULL;
2868     if(!(mctrackesd = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2869     TParticle *particle = 0x0;
2870     particle = mctrackesd->Particle();
2871     // Take mother
2872     if(!particle) return -1;
2873     Int_t imother   = particle->GetFirstMother();
2874     if(imother < 0) return -1;
2875     AliMCParticle *mothertrack = NULL;
2876     if(!(mothertrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2877     TParticle * mother = mothertrack->Particle();
2878     if(!mother) return -1;
2879     // Check gamma
2880     Int_t pdg = mother->GetPdgCode();
2881     if(TMath::Abs(pdg) == 221) return imother;
2882     if(TMath::Abs(pdg) == 11) {
2883       return IsMotherEta(imother,mcEvent);
2884     }
2885     return -1;
2886   }
2887
2888   if(mctrack->IsA() == AliAODMCParticle::Class()) {
2889     AliAODMCParticle *mctrackaod = NULL;
2890     if(!(mctrackaod = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2891     // Take mother
2892     Int_t imother = mctrackaod->GetMother();
2893     if(imother < 0) return -1;
2894     AliAODMCParticle *mothertrack = NULL;
2895     if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2896     // Check gamma
2897     Int_t pdg = mothertrack->GetPdgCode();
2898     if(TMath::Abs(pdg) == 221) return imother;
2899     if(TMath::Abs(pdg) == 11) {
2900       return IsMotherEta(imother,mcEvent);
2901     }
2902     return -1;
2903   }
2904
2905   return -1;
2906
2907 }