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