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