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",binLimAngle[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){
2099           hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
2100           if(((AliESDEvent*)fInputEvent)->GetPrimaryVertexSPD()) hfetrack.SetMulitplicity(((AliESDEvent*)fInputEvent)->GetPrimaryVertexSPD()->GetNContributors());
2101         }else{
2102           hfetrack.SetAnalysisType(AliHFEpidObject::kAODanalysis);
2103           if(((AliAODEvent*)fInputEvent)->GetPrimaryVertexSPD())  hfetrack.SetMulitplicity(((AliAODEvent*)fInputEvent)->GetPrimaryVertexSPD()->GetNContributors());
2104         }
2105         hfetrack.SetRecTrack(track);
2106         hfetrack.SetCentrality((Int_t)binct);
2107         AliDebug(2,Form("centrality %f and %d",binct,hfetrack.GetCentrality()));
2108         hfetrack.SetPbPb();
2109
2110         // Only TOF PID
2111         if(fMonitorContamination) {
2112           if(fPIDTOFOnly->IsSelected(&hfetrack,0x0,"recTrackCont",0x0)) {
2113             Float_t nsigma = pidResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
2114             valuedeltaphicontamination[3] = nsigma;
2115             fDeltaPhiMapsContamination->Fill(&valuedeltaphicontamination[0]);
2116           }
2117         }
2118
2119         // Complete PID TOF+TPC
2120         if(!fPID->IsSelected(&hfetrack,0x0,"recTrackCont",fPIDqa)) {
2121           continue;
2122         }
2123       }
2124       else {
2125         if(!mcEvent) continue;
2126         if(!(mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel()))))) continue;
2127         AliDebug(2,Form("PdgCode %d",TMath::Abs(mctrack->Particle()->GetPdgCode())));
2128         if(TMath::Abs(mctrack->Particle()->GetPdgCode())!=11) continue;
2129       }
2130     }
2131
2132
2133     /////////////////////////////////////////////////////////////////////////////
2134     // Add candidate to AliFlowEvent for POI and subtract from RP if needed
2135     ////////////////////////////////////////////////////////////////////////////
2136     if(fMonitorQCumulant) {
2137       Int_t idtrack = static_cast<AliVTrack*>(track)->GetID();
2138       Bool_t found = kFALSE;
2139       Int_t numberoffound = 0;
2140       AliDebug(2,Form("A: Number of tracks %d",fflowEvent->NumberOfTracks()));
2141       for(Int_t iRPs=0; iRPs< fflowEvent->NumberOfTracks(); iRPs++) {
2142         AliFlowTrack *iRP = (AliFlowTrack*) (fflowEvent->GetTrack(iRPs));
2143         if(!iRP) continue;
2144         //if(!iRP->InRPSelection()) continue;
2145         if( TMath::Abs(idtrack) == TMath::Abs(iRP->GetID()) ) {
2146           iRP->SetForPOISelection(kTRUE);
2147           found = kTRUE;
2148           numberoffound ++;
2149         }
2150       }
2151       AliDebug(2,Form("Found %d mal",numberoffound));
2152       if(!found) {
2153         AliFlowCandidateTrack *sTrack = (AliFlowCandidateTrack*) MakeTrack(massElectron,track->Pt(),track->Phi(), track->Eta());
2154         sTrack->SetID(idtrack);
2155         fflowEvent->AddTrack(sTrack);
2156         AliDebug(2,"Add the track");
2157       }
2158       AliDebug(2,Form("B: Number of tracks %d",fflowEvent->NumberOfTracks()));
2159     }
2160
2161
2162     /////////////////////
2163     // Fill THnSparseF
2164     /////////////////////
2165
2166     //
2167     valuensparseabis[0] = eventplanesubtracted;
2168     if((fillEventPlane) && (fMonitorEventPlane)) fEventPlaneaftersubtraction->Fill(&valuensparseabis[0]);
2169
2170
2171     if(fMonitorEventPlane)
2172       {
2173         //
2174         valuensparsee[0] = TMath::Cos(2*phitrack);
2175         fCos2phie->Fill(&valuensparsee[0]);
2176         valuensparsee[0] = TMath::Sin(2*phitrack);
2177         fSin2phie->Fill(&valuensparsee[0]);
2178         //
2179         valuensparsee[0] = TMath::Cos(2*eventplanesubtracted);
2180         if(fillEventPlane) fCos2phiep->Fill(&valuensparsee[0]);
2181         valuensparsee[0] = TMath::Sin(2*eventplanesubtracted);
2182         if(fillEventPlane) fSin2phiep->Fill(&valuensparsee[0]);
2183         valuensparsee[0] = TMath::Sin(2*TVector2::Phi_mpi_pi(phitrack-eventplanesubtracted));
2184         if(fillEventPlane) fSin2phiephiep->Fill(&valuensparsee[0]);
2185         //
2186       }
2187
2188     //
2189     valuensparseg[0] = deltaphi;
2190     if(fillEventPlane) fDeltaPhiMaps->Fill(&valuensparseg[0]);
2191
2192     //
2193     valuensparseh[0] = TMath::Cos(2*TVector2::Phi_mpi_pi(phitrack-eventplanesubtracted));
2194     if(fillEventPlane) {
2195       fCosPhiMaps->Fill(&valuensparseh[0]); //TR: fCosPhiQSum+=valuensparseh[0]*TMath:Sqrt(qAna->X()*qAna->X()+qAna->Y()*qAna->Y()); fCosPhiQN++;
2196       if((valuefractioncont[1] >=0) && (valuefractioncont[1] < 11)){
2197         if(fContamination[((Int_t)valuefractioncont[1])]){
2198           Double_t weight = 1.;
2199           if(fAsFunctionOfP) weight = fContamination[((Int_t)valuefractioncont[1])]->Eval(track->P());
2200           else weight = fContamination[((Int_t)valuefractioncont[1])]->Eval(track->Pt());
2201           if(weight<0.0) weight=0.0;
2202           if(weight>1.0) weight=1.0;
2203           fFractionContamination->Fill(&valuefractioncont[0],weight);
2204           if(fv2contamination[((Int_t)valuefractioncont[1])]){
2205             Double_t v2 =  fv2contamination[((Int_t)valuefractioncont[1])]->Eval(track->Pt());
2206             AliDebug(2,Form("value v2 %f, contamination %f and pt %f centrality %d\n",v2,weight,track->Pt(),(Int_t)valuefractioncont[1]));
2207             AliDebug(2,Form("Check for centrality 3: value v2 %f, contamination %f\n",fv2contamination[3]->Eval(track->Pt()),fContamination[3]->Eval(track->P())));
2208             AliDebug(2,Form("Check for centrality 4: value v2 %f, contamination %f\n",fv2contamination[4]->Eval(track->Pt()),fContamination[4]->Eval(track->P())));
2209             AliDebug(2,Form("Check for centrality 5: value v2 %f, contamination %f\n",fv2contamination[5]->Eval(track->Pt()),fContamination[5]->Eval(track->P())));
2210             fContaminationv2->Fill(valuefractioncont[1],valuefractioncont[0],v2,weight);
2211           }
2212           fContaminationmeanpt->Fill(valuefractioncont[1],valuefractioncont[0],TMath::Abs(track->Pt()));
2213         }
2214       }
2215       if(fMonitorEventPlane) {
2216         if(fSP)
2217           valuensparseh[0] *= TMath::Sqrt(qAna->X()*qAna->X()+qAna->Y()*qAna->Y());
2218         fProfileCosPhiMaps->Fill(valuensparsehprofile[1],valuensparsehprofile[2],valuensparseh[0]);  //TR: info
2219       }
2220     }
2221
2222     if(fMonitorPhotonic) {
2223       Int_t indexmother = -1;
2224       Int_t source = 1;
2225       if(mcthere) source = fBackgroundSubtraction->FindMother(mctrack->GetLabel(),indexmother);
2226       fBackgroundSubtraction->LookAtNonHFE(k, track, fInputEvent, 1, binct, deltaphi, source, indexmother);
2227
2228       if((!fAODAnalysis && mcthere) || !mcthere) {
2229         // background
2230         source = 0;
2231         indexmother = -1;
2232         source = FindMother(TMath::Abs(track->GetLabel()),mcEvent, indexmother);
2233         valuensparseMCSourceDeltaPhiMaps[2] = source;
2234         if(mcEvent) fMCSourceDeltaPhiMaps->Fill(&valuensparseMCSourceDeltaPhiMaps[0]);
2235         //LookAtNonHFE(k,track,fInputEvent,mcEvent,binct,deltaphi,source,indexmother);
2236         Int_t taggedvalue = LookAtNonHFE(k,track,fInputEvent,mcEvent,binct,deltaphi,source,indexmother);
2237         if(fMonitorPhotonic) {
2238           // No opposite charge partner found in the invariant mass choosen
2239           if((taggedvalue!=2) && (taggedvalue!=6)) {
2240             //fDeltaPhiMapsTaggedNonPhotonic->Fill(&valuensparseg[0]);
2241             //fCosPhiMapsTaggedNonPhotonic->Fill(&valuensparseh[0]);
2242           }
2243           // One opposite charge partner found in the invariant mass choosen
2244           if((taggedvalue==2) || (taggedvalue==6)) {
2245             fDeltaPhiMapsTaggedPhotonic->Fill(&valuensparseg[0]);
2246             //fCosPhiMapsTaggedPhotonic->Fill(&valuensparseh[0]);
2247           }
2248           // One same charge partner found in the invariant mass choosen
2249           if((taggedvalue==4) || (taggedvalue==6)) {
2250             fDeltaPhiMapsTaggedPhotonicLS->Fill(&valuensparseg[0]);
2251             //fCosPhiMapsTaggedPhotonicLS->Fill(&valuensparseh[0]);
2252           }
2253         }
2254       }
2255     }
2256
2257   }
2258
2259   //////////////////////////////////////////////////////////////////////////////
2260   ///////////////////////////AFTERBURNER
2261   if (fAfterBurnerOn &  fMonitorQCumulant)
2262     {
2263       fflowEvent->AddFlow(fV1,fV2,fV3,fV4,fV5);     //add flow
2264       fflowEvent->CloneTracks(fNonFlowNumberOfTrackClones); //add nonflow by cloning tracks
2265     }
2266   //////////////////////////////////////////////////////////////////////////////
2267
2268
2269
2270   //for(Int_t bincless = 0; bincless < fNbBinsCentralityQCumulant; bincless++) {
2271   //  if((fBinCentralityLess[bincless]< cntr) && (cntr < fBinCentralityLess[bincless+1])) PostData(bincless+2,fflowEvent);
2272   //}
2273
2274   if(fMonitorPhotonic) {
2275     if(fArraytrack) {
2276       delete fArraytrack;
2277       fArraytrack = NULL;
2278     }
2279   }
2280
2281   if(fMonitorPhotonic) fBackgroundSubtraction->CountPoolAssociated(fInputEvent,binct);
2282
2283   PostData(1, fListHist);
2284
2285
2286
2287 }
2288 //______________________________________________________________________________
2289 AliFlowCandidateTrack *AliAnalysisTaskFlowTPCTOFEPSP::MakeTrack( Double_t mass,
2290                           Double_t pt, Double_t phi, Double_t eta) {
2291   //
2292   //  Make Track (Not needed actually)
2293   //
2294
2295   AliFlowCandidateTrack *sTrack = new AliFlowCandidateTrack();
2296   sTrack->SetMass(mass);
2297   sTrack->SetPt(pt);
2298   sTrack->SetPhi(phi);
2299   sTrack->SetEta(eta);
2300   sTrack->SetForPOISelection(kTRUE);
2301   sTrack->SetForRPSelection(kFALSE);
2302   return sTrack;
2303 }
2304 //_________________________________________________________________________________
2305 Double_t AliAnalysisTaskFlowTPCTOFEPSP::GetPhiAfterAddV2(Double_t phi,Double_t reactionPlaneAngle) const
2306 {
2307   //
2308   // Adds v2, uses Newton-Raphson iteration
2309   //
2310   Double_t phiend=phi;
2311   Double_t phi0=phi;
2312   Double_t f=0.;
2313   Double_t fp=0.;
2314   Double_t phiprev=0.;
2315
2316   for (Int_t i=0; i<fMaxNumberOfIterations; i++)
2317   {
2318     phiprev=phiend; //store last value for comparison
2319     f =  phiend-phi0+fV2*TMath::Sin(2.*(phiend-reactionPlaneAngle));
2320     fp = 1.0+2.0*fV2*TMath::Cos(2.*(phiend-reactionPlaneAngle)); //first derivative
2321     phiend -= f/fp;
2322     if (TMath::AreEqualAbs(phiprev,phiend,fPrecisionPhi)) break;
2323   }
2324   return phiend;
2325 }
2326 //_____________________________________________________________________________________________
2327 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)
2328 {
2329   //
2330   // Look At Non HFE
2331   //
2332
2333   // return -1 if nothing
2334   // return 2 if opposite charge within the mass range found
2335   // return 4 if like charge within the mass range found
2336   // return 6 if opposite charge and like charge within the mass range found
2337   //
2338
2339   Int_t taggedphotonic = -1;
2340
2341   Bool_t oppositetaggedphotonic = kFALSE;
2342   Bool_t sametaggedphotonic = kFALSE;
2343
2344   AliDebug(2,Form("fCounterPoolBackground %d in LookAtNonHFE!!!",fCounterPoolBackground));
2345   if(!fArraytrack) return taggedphotonic;
2346   AliDebug(2,Form("process track %d",iTrack1));
2347
2348   TVector3 v3D1;
2349   TVector3 v3D2;
2350
2351   Double_t valuensparseDeltaPhiMaps[5];
2352   Double_t valueangle[3];
2353
2354   valuensparseDeltaPhiMaps[1] = binct;
2355   valuensparseDeltaPhiMaps[2] = track1->Pt();
2356   valuensparseDeltaPhiMaps[0] = deltaphi;
2357   valuensparseDeltaPhiMaps[4] = source;
2358
2359   valueangle[2] = source;
2360   valueangle[1] = binct;
2361
2362   // Pdg code
2363   Int_t pdg1 = CheckPdg(TMath::Abs(track1->GetLabel()),mcEvent);
2364   Int_t numberfound = 0;
2365
2366   //Magnetic Field
2367   Double_t bfield = vEvent->GetMagneticField();
2368
2369   // Get Primary vertex
2370   const AliVVertex *pVtx = vEvent->GetPrimaryVertex();
2371
2372   for(Int_t idex = 0; idex < fCounterPoolBackground; idex++)
2373     {
2374
2375       Int_t iTrack2 = fArraytrack->At(idex);
2376       AliDebug(2,Form("track %d",iTrack2));
2377       AliVTrack* track2 = (AliVTrack *) vEvent->GetTrack(iTrack2);
2378       if (!track2)
2379         {
2380           printf("ERROR: Could not receive track %d", iTrack2);
2381           continue;
2382         }
2383       if(iTrack2==iTrack1) continue;
2384       AliDebug(2,"Different");
2385
2386       // Reset the MC info
2387       valueangle[2] = source;
2388       valuensparseDeltaPhiMaps[4] = source;
2389
2390       // track cuts and PID already done
2391
2392       // if MC look
2393       Int_t pdg2 = -100;
2394       if(mcEvent) {
2395         Int_t source2 = 0;
2396         Int_t indexmother2 = -1;
2397         source2 = FindMother(TMath::Abs(track2->GetLabel()),mcEvent, indexmother2);
2398         pdg2 = CheckPdg(TMath::Abs(track2->GetLabel()),mcEvent);
2399         if(source2 >=0 ) {
2400           if((indexmother2 == indexmother) && (source == source2) && ((pdg1*pdg2)<0.0)) {
2401             if(source == kElectronfromconversion) {
2402               valueangle[2] = kElectronfromconversionboth;
2403               valuensparseDeltaPhiMaps[4] = kElectronfromconversionboth;
2404               numberfound++;
2405             }
2406             if(source == kElectronfrompi0) {
2407               valueangle[2] = kElectronfrompi0both;
2408               valuensparseDeltaPhiMaps[4] = kElectronfrompi0both;
2409             }
2410             if(source == kElectronfrometa) {
2411               valueangle[2] = kElectronfrometaboth;
2412               valuensparseDeltaPhiMaps[4] = kElectronfrometaboth;
2413             }
2414           }
2415         }
2416       }
2417
2418       if(fAlgorithmMA && (!fAODAnalysis))
2419         {
2420           // tracks
2421           AliESDtrack *esdtrack2 = dynamic_cast<AliESDtrack *>(track2);
2422           AliESDtrack *esdtrack1 = dynamic_cast<AliESDtrack *>(track1);
2423           if((!esdtrack2) || (!esdtrack1)) continue;
2424
2425           //Variables
2426           Double_t p1[3];
2427           Double_t p2[3];
2428           Double_t xt1; //radial position track 1 at the DCA point
2429           Double_t xt2; //radial position track 2 at the DCA point
2430           //DCA track1-track2
2431           Double_t dca12 = esdtrack2->GetDCA(esdtrack1,bfield,xt2,xt1);
2432
2433           // Cut dca
2434           if(dca12 > fMaxdca) continue;
2435
2436           //Momento of the track extrapolated to DCA track-track
2437           //Track1
2438           Bool_t hasdcaT1 = esdtrack1->GetPxPyPzAt(xt1,bfield,p1);
2439           //Track2
2440           Bool_t hasdcaT2 = esdtrack2->GetPxPyPzAt(xt2,bfield,p2);
2441
2442           if(!hasdcaT1 || !hasdcaT2) AliWarning("It could be a problem in the extrapolation");
2443
2444           //track1-track2 Invariant Mass
2445           Double_t eMass = 0.000510998910; //Electron mass in GeV
2446           Double_t pP1 = sqrt(p1[0]*p1[0]+p1[1]*p1[1]+p1[2]*p1[2]); //Track 1 momentum
2447           Double_t pP2 = sqrt(p2[0]*p2[0]+p2[1]*p2[1]+p2[2]*p2[2]); //Track 2 momentum
2448           Double_t eE1 = TMath::Sqrt(pP1*pP1+eMass*eMass);
2449           Double_t eE2 = TMath::Sqrt(pP2*pP2+eMass*eMass);
2450
2451           //TLorentzVector v1(p1[0],p1[1],p1[2],sqrt(eMass*eMass+pP1*pP1));
2452           //TLorentzVector v2(p2[0],p2[1],p2[2],sqrt(eMass*eMass+pP2*pP2));
2453           //Double_t imass = (v1+v2).M(); //Invariant Mass
2454           //Double_t angle3D = v1.Angle(v2.Vect()); //Opening Angle (Total Angle)
2455
2456           // daughter
2457           v3D1.SetXYZ(p1[0],p1[1],p1[2]);
2458           v3D2.SetXYZ(p2[0],p2[1],p2[2]);
2459           Double_t openingangle = TVector2::Phi_0_2pi(v3D2.Angle(v3D1));
2460
2461           // mother
2462           TVector3 motherrec = v3D1 + v3D2;
2463           Double_t invmass = TMath::Sqrt((eE1+eE2)*(eE1+eE2)-(motherrec.Px()*motherrec.Px()+motherrec.Py()*motherrec.Py()+motherrec.Pz()*motherrec.Pz()));
2464
2465           // xy
2466           //TVector3 vectordiff = v3D1 - v3D2;
2467           //Double_t diffphi = TVector2::Phi_0_2pi(vectordiff.Phi());
2468           //Double_t massxy = TMath::Sqrt((eE1+eE2)*(eE1+eE2)-(pP1*pP1+pP2*pP2+2*pP1*pP2*TMath::Cos(diffphi)));
2469
2470           // rz
2471           //Double_t difftheta = TVector2::Phi_0_2pi(vectordiff.Eta());
2472           //Double_t massrz = TMath::Sqrt((eE1+eE2)*(eE1+eE2)-(pP1*pP1+pP2*pP2+2*pP1*pP2*TMath::Cos(difftheta)));
2473
2474
2475           Float_t fCharge1 = track1->Charge();
2476           Float_t fCharge2 = track2->Charge();
2477
2478           // Fill Histo
2479           //valueangle[0] = diffphi;
2480           //valueangle[1] = difftheta;
2481           valueangle[0] = openingangle;
2482           if((fCharge1*fCharge2)>0.0) fSameSignAngle->Fill(&valueangle[0]);
2483           else fOppSignAngle->Fill(&valueangle[0]);
2484
2485           // Cut
2486           if(openingangle > fMaxopening3D) continue;
2487           //if(difftheta > fMaxopeningtheta) continue;
2488           //if(diffphi > fMaxopeningphi) continue;
2489
2490           // Invmass
2491           valuensparseDeltaPhiMaps[3] = invmass;
2492           if((fCharge1*fCharge2)>0.0) fSameSignDeltaPhiMaps->Fill(&valuensparseDeltaPhiMaps[0]);
2493           else fOppSignDeltaPhiMaps->Fill(&valuensparseDeltaPhiMaps[0]);
2494
2495           // Cut
2496           if(invmass < fMaxInvmass) {
2497             if((fCharge1*fCharge2)<0.0) oppositetaggedphotonic=kTRUE;
2498             if((fCharge1*fCharge2)>0.0) sametaggedphotonic=kTRUE;
2499           }
2500
2501
2502         }
2503       else
2504         {
2505           Int_t fPDGtrack1 = 11;
2506           Int_t fPDGtrack2 = 11;
2507
2508           Float_t fCharge1 = track1->Charge();
2509           Float_t fCharge2 = track2->Charge();
2510
2511           if(fCharge1>0) fPDGtrack1 = -11;
2512           if(fCharge2>0) fPDGtrack2 = -11;
2513
2514           AliKFParticle ktrack1(*track1, fPDGtrack1);
2515           AliKFParticle ktrack2(*track2, fPDGtrack2);
2516           AliKFParticle recoGamma(ktrack1, ktrack2);
2517
2518           //Reconstruction Cuts
2519           if(recoGamma.GetNDF()<1) continue;
2520           Double_t chi2OverNDF = recoGamma.GetChi2()/recoGamma.GetNDF();
2521           if(TMath::Sqrt(TMath::Abs(chi2OverNDF))>fChi2OverNDFCut) continue;
2522
2523           // DCA
2524           //Double_t dca12 = ktrack1.GetDistanceFromParticle(ktrack2);
2525           //if(dca12 > fMaxdca) continue;
2526
2527           // if set mass constraint
2528           if(fSetMassConstraint && pVtx) {
2529             AliKFVertex primV(*pVtx);
2530             primV += recoGamma;
2531             primV -= ktrack1;
2532             primV -= ktrack2;
2533             recoGamma.SetProductionVertex(primV);
2534             recoGamma.SetMassConstraint(0,0.0001);
2535           }
2536
2537           //Invariant Mass
2538           Double_t imass;
2539           Double_t width;
2540           recoGamma.GetMass(imass,width);
2541
2542           //Opening Angle (Total Angle)
2543           Double_t angle = ktrack1.GetAngle(ktrack2);
2544           valueangle[0] = angle;
2545           if((fCharge1*fCharge2)>0.0) fSameSignAngle->Fill(&valueangle[0]);
2546           else fOppSignAngle->Fill(&valueangle[0]);
2547
2548           // Cut
2549           if(angle > fMaxopening3D) continue;
2550
2551           // Invmass
2552           valuensparseDeltaPhiMaps[3] = imass;
2553           if((fCharge1*fCharge2)>0.0) fSameSignDeltaPhiMaps->Fill(&valuensparseDeltaPhiMaps[0]);
2554           else {
2555             fOppSignDeltaPhiMaps->Fill(&valuensparseDeltaPhiMaps[0]);
2556             /*
2557             if(valueangle[2] == kElectronfromconversionboth) {
2558               printf("Reconstructed charge1 %f, charge 2 %f and invmass %f",fCharge1,fCharge2,imass);
2559               printf("MC charge1 %d, charge 2 %d",pdg1,pdg2);
2560               printf("DCA %f",dca12);
2561               printf("Number of found %d",numberfound);
2562             }
2563             */
2564           }
2565
2566           // Cut
2567           if(imass < fMaxInvmass) {
2568             if((fCharge1*fCharge2)<0.0) oppositetaggedphotonic=kTRUE;
2569             if((fCharge1*fCharge2)>0.0) sametaggedphotonic=kTRUE;
2570           }
2571         }
2572     }
2573
2574   if(oppositetaggedphotonic && sametaggedphotonic){
2575     taggedphotonic = 6;
2576   }
2577
2578   if(!oppositetaggedphotonic && sametaggedphotonic){
2579     taggedphotonic = 4;
2580   }
2581
2582   if(oppositetaggedphotonic && !sametaggedphotonic){
2583     taggedphotonic = 2;
2584   }
2585
2586
2587   return taggedphotonic;
2588 }
2589 //_________________________________________________________________________
2590 Int_t AliAnalysisTaskFlowTPCTOFEPSP::FindMother(Int_t tr, AliMCEvent *mcEvent, Int_t &indexmother){
2591   //
2592   // Find the mother if MC
2593   //
2594
2595   if(!mcEvent) return 0;
2596
2597   Int_t pdg = CheckPdg(tr,mcEvent);
2598   if(TMath::Abs(pdg)!= 11) {
2599     indexmother = -1;
2600     return kNoElectron;
2601   }
2602
2603   indexmother = IsMotherGamma(tr,mcEvent);
2604   if(indexmother > 0) return kElectronfromconversion;
2605   indexmother = IsMotherPi0(tr,mcEvent);
2606   if(indexmother > 0) return kElectronfrompi0;
2607   indexmother = IsMotherC(tr,mcEvent);
2608   if(indexmother > 0) return kElectronfromC;
2609   indexmother = IsMotherB(tr,mcEvent);
2610   if(indexmother > 0) return kElectronfromB;
2611   indexmother = IsMotherEta(tr,mcEvent);
2612   if(indexmother > 0) return kElectronfrometa;
2613
2614   return kElectronfromother;
2615
2616
2617 }
2618 //____________________________________________________________________________________________________________
2619 Int_t AliAnalysisTaskFlowTPCTOFEPSP::CheckPdg(Int_t tr, AliMCEvent* mcEvent) {
2620
2621   //
2622   // Return the pdg of the particle
2623   //
2624
2625
2626   Int_t pdgcode = -1;
2627   if(tr < 0) return pdgcode;
2628
2629   if(!mcEvent) return pdgcode;
2630
2631   AliVParticle *mctrack = mcEvent->GetTrack(tr);
2632
2633
2634   if(mctrack->IsA() == AliMCParticle::Class()) {
2635     AliMCParticle *mctrackesd = NULL;
2636     if(!(mctrackesd = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return pdgcode;
2637     pdgcode = mctrackesd->PdgCode();
2638   }
2639
2640   if(mctrack->IsA() == AliAODMCParticle::Class()) {
2641     AliAODMCParticle *mctrackaod = NULL;
2642     if(!(mctrackaod = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return pdgcode;
2643     pdgcode = mctrackaod->GetPdgCode();
2644   }
2645
2646   return pdgcode;
2647
2648
2649 }
2650 //____________________________________________________________________________________________________________
2651 Int_t AliAnalysisTaskFlowTPCTOFEPSP::IsMotherGamma(Int_t tr, AliMCEvent* mcEvent) {
2652
2653   //
2654   // Return the lab of gamma mother or -1 if not gamma
2655   //
2656
2657   if(tr < 0) return -1;
2658   AliVParticle *mctrack = mcEvent->GetTrack(tr);
2659
2660   if(mctrack->IsA() == AliMCParticle::Class()) {
2661     AliMCParticle *mctrackesd = NULL;
2662     if(!(mctrackesd = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2663     TParticle *particle = 0x0;
2664     particle = mctrackesd->Particle();
2665     // Take mother
2666     if(!particle) return -1;
2667     Int_t imother   = particle->GetFirstMother();
2668     if(imother < 0) return -1;
2669     AliMCParticle *mothertrack = NULL;
2670     if(!(mothertrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2671     TParticle * mother = mothertrack->Particle();
2672     if(!mother) return -1;
2673     // Check gamma
2674     Int_t pdg = mother->GetPdgCode();
2675     if(TMath::Abs(pdg) == 22) return imother;
2676     if(TMath::Abs(pdg) == 11) {
2677       return IsMotherGamma(imother,mcEvent);
2678     }
2679     return -1;
2680   }
2681
2682   if(mctrack->IsA() == AliAODMCParticle::Class()) {
2683     AliAODMCParticle *mctrackaod = NULL;
2684     if(!(mctrackaod = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2685     // Take mother
2686     Int_t imother = mctrackaod->GetMother();
2687     if(imother < 0) return -1;
2688     AliAODMCParticle *mothertrack = NULL;
2689     if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2690     // Check gamma
2691     Int_t pdg = mothertrack->GetPdgCode();
2692     if(TMath::Abs(pdg) == 22) return imother;
2693     if(TMath::Abs(pdg) == 11) {
2694       return IsMotherGamma(imother,mcEvent);
2695     }
2696     return -1;
2697
2698   }
2699
2700   return -1;
2701
2702
2703 }
2704 //
2705 //____________________________________________________________________________________________________________
2706 Int_t AliAnalysisTaskFlowTPCTOFEPSP::IsMotherPi0(Int_t tr, AliMCEvent* mcEvent) {
2707
2708   //
2709   // Return the lab of pi0 mother or -1 if not pi0
2710   //
2711
2712   if(tr < 0) return -1;
2713   AliVParticle *mctrack = mcEvent->GetTrack(tr);
2714
2715   if(mctrack->IsA() == AliMCParticle::Class()) {
2716     AliMCParticle *mctrackesd = NULL;
2717     if(!(mctrackesd = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2718     TParticle *particle = 0x0;
2719     particle = mctrackesd->Particle();
2720     // Take mother
2721     if(!particle) return -1;
2722     Int_t imother   = particle->GetFirstMother();
2723     if(imother < 0) return -1;
2724     AliMCParticle *mothertrack = NULL;
2725     if(!(mothertrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2726     TParticle * mother = mothertrack->Particle();
2727     if(!mother) return -1;
2728     // Check gamma
2729     Int_t pdg = mother->GetPdgCode();
2730     if(TMath::Abs(pdg) == 111) return imother;
2731     if(TMath::Abs(pdg) == 11) {
2732       return IsMotherPi0(imother,mcEvent);
2733     }
2734     return -1;
2735   }
2736
2737   if(mctrack->IsA() == AliAODMCParticle::Class()) {
2738     AliAODMCParticle *mctrackaod = NULL;
2739     if(!(mctrackaod = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2740     // Take mother
2741     Int_t imother = mctrackaod->GetMother();
2742     if(imother < 0) return -1;
2743     AliAODMCParticle *mothertrack = NULL;
2744     if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2745     // Check gamma
2746     Int_t pdg = mothertrack->GetPdgCode();
2747     if(TMath::Abs(pdg) == 111) return imother;
2748     if(TMath::Abs(pdg) == 11) {
2749       return IsMotherPi0(imother,mcEvent);
2750     }
2751     return -1;
2752   }
2753
2754   return -1;
2755
2756 }
2757 //____________________________________________________________________________________________________________
2758 Int_t AliAnalysisTaskFlowTPCTOFEPSP::IsMotherC(Int_t tr, AliMCEvent* mcEvent) {
2759
2760   //
2761   // Return the lab of signal mother or -1 if not signal
2762   //
2763
2764   if(tr < 0) return -1;
2765   AliVParticle *mctrack = mcEvent->GetTrack(tr);
2766
2767   if(mctrack->IsA() == AliMCParticle::Class()) {
2768     AliMCParticle *mctrackesd = NULL;
2769     if(!(mctrackesd = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2770     TParticle *particle = 0x0;
2771     particle = mctrackesd->Particle();
2772     // Take mother
2773     if(!particle) return -1;
2774     Int_t imother   = particle->GetFirstMother();
2775     if(imother < 0) return -1;
2776     AliMCParticle *mothertrack = NULL;
2777     if(!(mothertrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2778     TParticle * mother = mothertrack->Particle();
2779     if(!mother) return -1;
2780     // Check gamma
2781     Int_t pdg = mother->GetPdgCode();
2782     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;
2783     if(TMath::Abs(pdg) == 11) {
2784       return IsMotherC(imother,mcEvent);
2785     }
2786     return -1;
2787   }
2788
2789   if(mctrack->IsA() == AliAODMCParticle::Class()) {
2790     AliAODMCParticle *mctrackaod = NULL;
2791     if(!(mctrackaod = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2792     // Take mother
2793     Int_t imother = mctrackaod->GetMother();
2794     if(imother < 0) return -1;
2795     AliAODMCParticle *mothertrack = NULL;
2796     if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2797     // Check gamma
2798     Int_t pdg = mothertrack->GetPdgCode();
2799     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;
2800     if(TMath::Abs(pdg) == 11) {
2801       return IsMotherC(imother,mcEvent);
2802     }
2803     return -1;
2804   }
2805
2806   return -1;
2807
2808 }
2809 //____________________________________________________________________________________________________________
2810 Int_t AliAnalysisTaskFlowTPCTOFEPSP::IsMotherB(Int_t tr, AliMCEvent* mcEvent) {
2811
2812   //
2813   // Return the lab of signal mother or -1 if not signal
2814   //
2815
2816   if(tr < 0) return -1;
2817   AliVParticle *mctrack = mcEvent->GetTrack(tr);
2818
2819   if(mctrack->IsA() == AliMCParticle::Class()) {
2820     AliMCParticle *mctrackesd = NULL;
2821     if(!(mctrackesd = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2822     TParticle *particle = 0x0;
2823     particle = mctrackesd->Particle();
2824     // Take mother
2825     if(!particle) return -1;
2826     Int_t imother   = particle->GetFirstMother();
2827     if(imother < 0) return -1;
2828     AliMCParticle *mothertrack = NULL;
2829     if(!(mothertrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2830     TParticle * mother = mothertrack->Particle();
2831     if(!mother) return -1;
2832     // Check gamma
2833     Int_t pdg = mother->GetPdgCode();
2834     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;
2835     if(TMath::Abs(pdg) == 11) {
2836       return IsMotherB(imother,mcEvent);
2837     }
2838     return -1;
2839   }
2840
2841   if(mctrack->IsA() == AliAODMCParticle::Class()) {
2842     AliAODMCParticle *mctrackaod = NULL;
2843     if(!(mctrackaod = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2844     // Take mother
2845     Int_t imother = mctrackaod->GetMother();
2846     if(imother < 0) return -1;
2847     AliAODMCParticle *mothertrack = NULL;
2848     if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2849     // Check gamma
2850     Int_t pdg = mothertrack->GetPdgCode();
2851     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;
2852     if(TMath::Abs(pdg) == 11) {
2853       return IsMotherB(imother,mcEvent);
2854     }
2855     return -1;
2856   }
2857
2858   return -1;
2859
2860 }
2861 //____________________________________________________________________________________________________________
2862 Int_t AliAnalysisTaskFlowTPCTOFEPSP::IsMotherEta(Int_t tr, AliMCEvent* mcEvent) {
2863
2864   //
2865   // Return the lab of pi0 mother or -1 if not pi0
2866   //
2867
2868  if(tr < 0) return -1;
2869   AliVParticle *mctrack = mcEvent->GetTrack(tr);
2870
2871   if(mctrack->IsA() == AliMCParticle::Class()) {
2872     AliMCParticle *mctrackesd = NULL;
2873     if(!(mctrackesd = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2874     TParticle *particle = 0x0;
2875     particle = mctrackesd->Particle();
2876     // Take mother
2877     if(!particle) return -1;
2878     Int_t imother   = particle->GetFirstMother();
2879     if(imother < 0) return -1;
2880     AliMCParticle *mothertrack = NULL;
2881     if(!(mothertrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2882     TParticle * mother = mothertrack->Particle();
2883     if(!mother) return -1;
2884     // Check gamma
2885     Int_t pdg = mother->GetPdgCode();
2886     if(TMath::Abs(pdg) == 221) return imother;
2887     if(TMath::Abs(pdg) == 11) {
2888       return IsMotherEta(imother,mcEvent);
2889     }
2890     return -1;
2891   }
2892
2893   if(mctrack->IsA() == AliAODMCParticle::Class()) {
2894     AliAODMCParticle *mctrackaod = NULL;
2895     if(!(mctrackaod = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2896     // Take mother
2897     Int_t imother = mctrackaod->GetMother();
2898     if(imother < 0) return -1;
2899     AliAODMCParticle *mothertrack = NULL;
2900     if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2901     // Check gamma
2902     Int_t pdg = mothertrack->GetPdgCode();
2903     if(TMath::Abs(pdg) == 221) return imother;
2904     if(TMath::Abs(pdg) == 11) {
2905       return IsMotherEta(imother,mcEvent);
2906     }
2907     return -1;
2908   }
2909
2910   return -1;
2911
2912 }