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