PWGHF/hfe/AliAnalysisTaskHFEFlow.cxx

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