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