PWGHF/hfe/AliAnalysisTaskFlowTPCTOFEPSP.cxx

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