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