PWGHF/hfe/AliAnalysisTaskFlowTPCTOFEPSP.cxx

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