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