]> 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
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[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;
2c0951cd 2098 if(!fAODAnalysis){
2099 hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
2100 if(((AliESDEvent*)fInputEvent)->GetPrimaryVertexSPD()) hfetrack.SetMulitplicity(((AliESDEvent*)fInputEvent)->GetPrimaryVertexSPD()->GetNContributors());
2101 }else{
2102 hfetrack.SetAnalysisType(AliHFEpidObject::kAODanalysis);
2103 if(((AliAODEvent*)fInputEvent)->GetPrimaryVertexSPD()) hfetrack.SetMulitplicity(((AliAODEvent*)fInputEvent)->GetPrimaryVertexSPD()->GetNContributors());
2104 }
c683985a 2105 hfetrack.SetRecTrack(track);
2106 hfetrack.SetCentrality((Int_t)binct);
2107 AliDebug(2,Form("centrality %f and %d",binct,hfetrack.GetCentrality()));
2108 hfetrack.SetPbPb();
2109
2110 // Only TOF PID
2111 if(fMonitorContamination) {
2112 if(fPIDTOFOnly->IsSelected(&hfetrack,0x0,"recTrackCont",0x0)) {
2113 Float_t nsigma = pidResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
2114 valuedeltaphicontamination[3] = nsigma;
2115 fDeltaPhiMapsContamination->Fill(&valuedeltaphicontamination[0]);
2116 }
2117 }
2118
2119 // Complete PID TOF+TPC
2120 if(!fPID->IsSelected(&hfetrack,0x0,"recTrackCont",fPIDqa)) {
2121 continue;
2122 }
2123 }
2124 else {
2125 if(!mcEvent) continue;
2126 if(!(mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel()))))) continue;
2127 AliDebug(2,Form("PdgCode %d",TMath::Abs(mctrack->Particle()->GetPdgCode())));
2128 if(TMath::Abs(mctrack->Particle()->GetPdgCode())!=11) continue;
2129 }
2130 }
2131
2132
2133 /////////////////////////////////////////////////////////////////////////////
2134 // Add candidate to AliFlowEvent for POI and subtract from RP if needed
2135 ////////////////////////////////////////////////////////////////////////////
2136 if(fMonitorQCumulant) {
2137 Int_t idtrack = static_cast<AliVTrack*>(track)->GetID();
2138 Bool_t found = kFALSE;
2139 Int_t numberoffound = 0;
2140 AliDebug(2,Form("A: Number of tracks %d",fflowEvent->NumberOfTracks()));
2141 for(Int_t iRPs=0; iRPs< fflowEvent->NumberOfTracks(); iRPs++) {
2142 AliFlowTrack *iRP = (AliFlowTrack*) (fflowEvent->GetTrack(iRPs));
2143 if(!iRP) continue;
2144 //if(!iRP->InRPSelection()) continue;
2145 if( TMath::Abs(idtrack) == TMath::Abs(iRP->GetID()) ) {
2146 iRP->SetForPOISelection(kTRUE);
2147 found = kTRUE;
2148 numberoffound ++;
2149 }
2150 }
2151 AliDebug(2,Form("Found %d mal",numberoffound));
2152 if(!found) {
2153 AliFlowCandidateTrack *sTrack = (AliFlowCandidateTrack*) MakeTrack(massElectron,track->Pt(),track->Phi(), track->Eta());
2154 sTrack->SetID(idtrack);
2155 fflowEvent->AddTrack(sTrack);
2156 AliDebug(2,"Add the track");
2157 }
2158 AliDebug(2,Form("B: Number of tracks %d",fflowEvent->NumberOfTracks()));
2159 }
2160
2161
2162 /////////////////////
2163 // Fill THnSparseF
2164 /////////////////////
2165
2166 //
2167 valuensparseabis[0] = eventplanesubtracted;
2168 if((fillEventPlane) && (fMonitorEventPlane)) fEventPlaneaftersubtraction->Fill(&valuensparseabis[0]);
2169
2170
2171 if(fMonitorEventPlane)
2172 {
2173 //
2174 valuensparsee[0] = TMath::Cos(2*phitrack);
2175 fCos2phie->Fill(&valuensparsee[0]);
2176 valuensparsee[0] = TMath::Sin(2*phitrack);
2177 fSin2phie->Fill(&valuensparsee[0]);
2178 //
2179 valuensparsee[0] = TMath::Cos(2*eventplanesubtracted);
2180 if(fillEventPlane) fCos2phiep->Fill(&valuensparsee[0]);
2181 valuensparsee[0] = TMath::Sin(2*eventplanesubtracted);
2182 if(fillEventPlane) fSin2phiep->Fill(&valuensparsee[0]);
2183 valuensparsee[0] = TMath::Sin(2*TVector2::Phi_mpi_pi(phitrack-eventplanesubtracted));
2184 if(fillEventPlane) fSin2phiephiep->Fill(&valuensparsee[0]);
2185 //
2186 }
2187
2188 //
2189 valuensparseg[0] = deltaphi;
2190 if(fillEventPlane) fDeltaPhiMaps->Fill(&valuensparseg[0]);
2191
2192 //
2193 valuensparseh[0] = TMath::Cos(2*TVector2::Phi_mpi_pi(phitrack-eventplanesubtracted));
2194 if(fillEventPlane) {
2195 fCosPhiMaps->Fill(&valuensparseh[0]); //TR: fCosPhiQSum+=valuensparseh[0]*TMath:Sqrt(qAna->X()*qAna->X()+qAna->Y()*qAna->Y()); fCosPhiQN++;
2196 if((valuefractioncont[1] >=0) && (valuefractioncont[1] < 11)){
2197 if(fContamination[((Int_t)valuefractioncont[1])]){
58a496d1 2198 Double_t weight = 1.;
2199 if(fAsFunctionOfP) weight = fContamination[((Int_t)valuefractioncont[1])]->Eval(track->P());
2200 else weight = fContamination[((Int_t)valuefractioncont[1])]->Eval(track->Pt());
c683985a 2201 if(weight<0.0) weight=0.0;
2202 if(weight>1.0) weight=1.0;
2203 fFractionContamination->Fill(&valuefractioncont[0],weight);
2204 if(fv2contamination[((Int_t)valuefractioncont[1])]){
2205 Double_t v2 = fv2contamination[((Int_t)valuefractioncont[1])]->Eval(track->Pt());
2206 AliDebug(2,Form("value v2 %f, contamination %f and pt %f centrality %d\n",v2,weight,track->Pt(),(Int_t)valuefractioncont[1]));
2207 AliDebug(2,Form("Check for centrality 3: value v2 %f, contamination %f\n",fv2contamination[3]->Eval(track->Pt()),fContamination[3]->Eval(track->P())));
2208 AliDebug(2,Form("Check for centrality 4: value v2 %f, contamination %f\n",fv2contamination[4]->Eval(track->Pt()),fContamination[4]->Eval(track->P())));
2209 AliDebug(2,Form("Check for centrality 5: value v2 %f, contamination %f\n",fv2contamination[5]->Eval(track->Pt()),fContamination[5]->Eval(track->P())));
2210 fContaminationv2->Fill(valuefractioncont[1],valuefractioncont[0],v2,weight);
2211 }
58a496d1 2212 fContaminationmeanpt->Fill(valuefractioncont[1],valuefractioncont[0],TMath::Abs(track->Pt()));
c683985a 2213 }
2214 }
2215 if(fMonitorEventPlane) {
2216 if(fSP)
2217 valuensparseh[0] *= TMath::Sqrt(qAna->X()*qAna->X()+qAna->Y()*qAna->Y());
2218 fProfileCosPhiMaps->Fill(valuensparsehprofile[1],valuensparsehprofile[2],valuensparseh[0]); //TR: info
2219 }
2220 }
2221
2222 if(fMonitorPhotonic) {
2223 Int_t indexmother = -1;
2224 Int_t source = 1;
2225 if(mcthere) source = fBackgroundSubtraction->FindMother(mctrack->GetLabel(),indexmother);
2226 fBackgroundSubtraction->LookAtNonHFE(k, track, fInputEvent, 1, binct, deltaphi, source, indexmother);
2227
2228 if((!fAODAnalysis && mcthere) || !mcthere) {
2229 // background
2230 source = 0;
2231 indexmother = -1;
2232 source = FindMother(TMath::Abs(track->GetLabel()),mcEvent, indexmother);
2233 valuensparseMCSourceDeltaPhiMaps[2] = source;
2234 if(mcEvent) fMCSourceDeltaPhiMaps->Fill(&valuensparseMCSourceDeltaPhiMaps[0]);
2235 //LookAtNonHFE(k,track,fInputEvent,mcEvent,binct,deltaphi,source,indexmother);
2236 Int_t taggedvalue = LookAtNonHFE(k,track,fInputEvent,mcEvent,binct,deltaphi,source,indexmother);
2237 if(fMonitorPhotonic) {
2238 // No opposite charge partner found in the invariant mass choosen
2239 if((taggedvalue!=2) && (taggedvalue!=6)) {
2240 //fDeltaPhiMapsTaggedNonPhotonic->Fill(&valuensparseg[0]);
2241 //fCosPhiMapsTaggedNonPhotonic->Fill(&valuensparseh[0]);
2242 }
2243 // One opposite charge partner found in the invariant mass choosen
2244 if((taggedvalue==2) || (taggedvalue==6)) {
2245 fDeltaPhiMapsTaggedPhotonic->Fill(&valuensparseg[0]);
2246 //fCosPhiMapsTaggedPhotonic->Fill(&valuensparseh[0]);
2247 }
2248 // One same charge partner found in the invariant mass choosen
2249 if((taggedvalue==4) || (taggedvalue==6)) {
2250 fDeltaPhiMapsTaggedPhotonicLS->Fill(&valuensparseg[0]);
2251 //fCosPhiMapsTaggedPhotonicLS->Fill(&valuensparseh[0]);
2252 }
2253 }
2254 }
2255 }
2256
2257 }
2258
2259 //////////////////////////////////////////////////////////////////////////////
2260 ///////////////////////////AFTERBURNER
2261 if (fAfterBurnerOn & fMonitorQCumulant)
2262 {
2263 fflowEvent->AddFlow(fV1,fV2,fV3,fV4,fV5); //add flow
2264 fflowEvent->CloneTracks(fNonFlowNumberOfTrackClones); //add nonflow by cloning tracks
2265 }
2266 //////////////////////////////////////////////////////////////////////////////
2267
2268
2269
2270 //for(Int_t bincless = 0; bincless < fNbBinsCentralityQCumulant; bincless++) {
2271 // if((fBinCentralityLess[bincless]< cntr) && (cntr < fBinCentralityLess[bincless+1])) PostData(bincless+2,fflowEvent);
2272 //}
2273
2274 if(fMonitorPhotonic) {
2275 if(fArraytrack) {
2276 delete fArraytrack;
2277 fArraytrack = NULL;
2278 }
2279 }
2280
2281 if(fMonitorPhotonic) fBackgroundSubtraction->CountPoolAssociated(fInputEvent,binct);
2282
2283 PostData(1, fListHist);
2284
2285
2286
2287}
2288//______________________________________________________________________________
2289AliFlowCandidateTrack *AliAnalysisTaskFlowTPCTOFEPSP::MakeTrack( Double_t mass,
2290 Double_t pt, Double_t phi, Double_t eta) {
2291 //
2292 // Make Track (Not needed actually)
2293 //
2294
2295 AliFlowCandidateTrack *sTrack = new AliFlowCandidateTrack();
2296 sTrack->SetMass(mass);
2297 sTrack->SetPt(pt);
2298 sTrack->SetPhi(phi);
2299 sTrack->SetEta(eta);
2300 sTrack->SetForPOISelection(kTRUE);
2301 sTrack->SetForRPSelection(kFALSE);
2302 return sTrack;
2303}
2304//_________________________________________________________________________________
2305Double_t AliAnalysisTaskFlowTPCTOFEPSP::GetPhiAfterAddV2(Double_t phi,Double_t reactionPlaneAngle) const
2306{
2307 //
2308 // Adds v2, uses Newton-Raphson iteration
2309 //
2310 Double_t phiend=phi;
2311 Double_t phi0=phi;
2312 Double_t f=0.;
2313 Double_t fp=0.;
2314 Double_t phiprev=0.;
2315
2316 for (Int_t i=0; i<fMaxNumberOfIterations; i++)
2317 {
2318 phiprev=phiend; //store last value for comparison
2319 f = phiend-phi0+fV2*TMath::Sin(2.*(phiend-reactionPlaneAngle));
2320 fp = 1.0+2.0*fV2*TMath::Cos(2.*(phiend-reactionPlaneAngle)); //first derivative
2321 phiend -= f/fp;
2322 if (TMath::AreEqualAbs(phiprev,phiend,fPrecisionPhi)) break;
2323 }
2324 return phiend;
2325}
2326//_____________________________________________________________________________________________
2327Int_t AliAnalysisTaskFlowTPCTOFEPSP::LookAtNonHFE(Int_t iTrack1, AliVTrack *track1, AliVEvent *vEvent, AliMCEvent *mcEvent,Int_t binct,Double_t deltaphi,Int_t source,Int_t indexmother)
2328{
2329 //
2330 // Look At Non HFE
2331 //
2332
2333 // return -1 if nothing
2334 // return 2 if opposite charge within the mass range found
2335 // return 4 if like charge within the mass range found
2336 // return 6 if opposite charge and like charge within the mass range found
2337 //
2338
2339 Int_t taggedphotonic = -1;
2340
2341 Bool_t oppositetaggedphotonic = kFALSE;
2342 Bool_t sametaggedphotonic = kFALSE;
2343
2344 AliDebug(2,Form("fCounterPoolBackground %d in LookAtNonHFE!!!",fCounterPoolBackground));
2345 if(!fArraytrack) return taggedphotonic;
2346 AliDebug(2,Form("process track %d",iTrack1));
2347
2348 TVector3 v3D1;
2349 TVector3 v3D2;
2350
2351 Double_t valuensparseDeltaPhiMaps[5];
2352 Double_t valueangle[3];
2353
2354 valuensparseDeltaPhiMaps[1] = binct;
2355 valuensparseDeltaPhiMaps[2] = track1->Pt();
2356 valuensparseDeltaPhiMaps[0] = deltaphi;
2357 valuensparseDeltaPhiMaps[4] = source;
2358
2359 valueangle[2] = source;
2360 valueangle[1] = binct;
2361
2362 // Pdg code
2363 Int_t pdg1 = CheckPdg(TMath::Abs(track1->GetLabel()),mcEvent);
2364 Int_t numberfound = 0;
2365
2366 //Magnetic Field
2367 Double_t bfield = vEvent->GetMagneticField();
2368
2369 // Get Primary vertex
2370 const AliVVertex *pVtx = vEvent->GetPrimaryVertex();
2371
2372 for(Int_t idex = 0; idex < fCounterPoolBackground; idex++)
2373 {
2374
2375 Int_t iTrack2 = fArraytrack->At(idex);
2376 AliDebug(2,Form("track %d",iTrack2));
2377 AliVTrack* track2 = (AliVTrack *) vEvent->GetTrack(iTrack2);
2378 if (!track2)
2379 {
2380 printf("ERROR: Could not receive track %d", iTrack2);
2381 continue;
2382 }
2383 if(iTrack2==iTrack1) continue;
2384 AliDebug(2,"Different");
2385
2386 // Reset the MC info
2387 valueangle[2] = source;
2388 valuensparseDeltaPhiMaps[4] = source;
2389
2390 // track cuts and PID already done
2391
2392 // if MC look
2393 Int_t pdg2 = -100;
2394 if(mcEvent) {
2395 Int_t source2 = 0;
2396 Int_t indexmother2 = -1;
2397 source2 = FindMother(TMath::Abs(track2->GetLabel()),mcEvent, indexmother2);
2398 pdg2 = CheckPdg(TMath::Abs(track2->GetLabel()),mcEvent);
2399 if(source2 >=0 ) {
2400 if((indexmother2 == indexmother) && (source == source2) && ((pdg1*pdg2)<0.0)) {
2401 if(source == kElectronfromconversion) {
2402 valueangle[2] = kElectronfromconversionboth;
2403 valuensparseDeltaPhiMaps[4] = kElectronfromconversionboth;
2404 numberfound++;
2405 }
2406 if(source == kElectronfrompi0) {
2407 valueangle[2] = kElectronfrompi0both;
2408 valuensparseDeltaPhiMaps[4] = kElectronfrompi0both;
2409 }
2410 if(source == kElectronfrometa) {
2411 valueangle[2] = kElectronfrometaboth;
2412 valuensparseDeltaPhiMaps[4] = kElectronfrometaboth;
2413 }
2414 }
2415 }
2416 }
2417
2418 if(fAlgorithmMA && (!fAODAnalysis))
2419 {
2420 // tracks
2421 AliESDtrack *esdtrack2 = dynamic_cast<AliESDtrack *>(track2);
2422 AliESDtrack *esdtrack1 = dynamic_cast<AliESDtrack *>(track1);
2423 if((!esdtrack2) || (!esdtrack1)) continue;
2424
2425 //Variables
2426 Double_t p1[3];
2427 Double_t p2[3];
2428 Double_t xt1; //radial position track 1 at the DCA point
2429 Double_t xt2; //radial position track 2 at the DCA point
2430 //DCA track1-track2
2431 Double_t dca12 = esdtrack2->GetDCA(esdtrack1,bfield,xt2,xt1);
2432
2433 // Cut dca
2434 if(dca12 > fMaxdca) continue;
2435
2436 //Momento of the track extrapolated to DCA track-track
2437 //Track1
2438 Bool_t hasdcaT1 = esdtrack1->GetPxPyPzAt(xt1,bfield,p1);
2439 //Track2
2440 Bool_t hasdcaT2 = esdtrack2->GetPxPyPzAt(xt2,bfield,p2);
2441
2442 if(!hasdcaT1 || !hasdcaT2) AliWarning("It could be a problem in the extrapolation");
2443
2444 //track1-track2 Invariant Mass
2445 Double_t eMass = 0.000510998910; //Electron mass in GeV
2446 Double_t pP1 = sqrt(p1[0]*p1[0]+p1[1]*p1[1]+p1[2]*p1[2]); //Track 1 momentum
2447 Double_t pP2 = sqrt(p2[0]*p2[0]+p2[1]*p2[1]+p2[2]*p2[2]); //Track 2 momentum
2448 Double_t eE1 = TMath::Sqrt(pP1*pP1+eMass*eMass);
2449 Double_t eE2 = TMath::Sqrt(pP2*pP2+eMass*eMass);
2450
2451 //TLorentzVector v1(p1[0],p1[1],p1[2],sqrt(eMass*eMass+pP1*pP1));
2452 //TLorentzVector v2(p2[0],p2[1],p2[2],sqrt(eMass*eMass+pP2*pP2));
2453 //Double_t imass = (v1+v2).M(); //Invariant Mass
2454 //Double_t angle3D = v1.Angle(v2.Vect()); //Opening Angle (Total Angle)
2455
2456 // daughter
2457 v3D1.SetXYZ(p1[0],p1[1],p1[2]);
2458 v3D2.SetXYZ(p2[0],p2[1],p2[2]);
2459 Double_t openingangle = TVector2::Phi_0_2pi(v3D2.Angle(v3D1));
2460
2461 // mother
2462 TVector3 motherrec = v3D1 + v3D2;
2463 Double_t invmass = TMath::Sqrt((eE1+eE2)*(eE1+eE2)-(motherrec.Px()*motherrec.Px()+motherrec.Py()*motherrec.Py()+motherrec.Pz()*motherrec.Pz()));
2464
2465 // xy
2466 //TVector3 vectordiff = v3D1 - v3D2;
2467 //Double_t diffphi = TVector2::Phi_0_2pi(vectordiff.Phi());
2468 //Double_t massxy = TMath::Sqrt((eE1+eE2)*(eE1+eE2)-(pP1*pP1+pP2*pP2+2*pP1*pP2*TMath::Cos(diffphi)));
2469
2470 // rz
2471 //Double_t difftheta = TVector2::Phi_0_2pi(vectordiff.Eta());
2472 //Double_t massrz = TMath::Sqrt((eE1+eE2)*(eE1+eE2)-(pP1*pP1+pP2*pP2+2*pP1*pP2*TMath::Cos(difftheta)));
2473
2474
2475 Float_t fCharge1 = track1->Charge();
2476 Float_t fCharge2 = track2->Charge();
2477
2478 // Fill Histo
2479 //valueangle[0] = diffphi;
2480 //valueangle[1] = difftheta;
2481 valueangle[0] = openingangle;
2482 if((fCharge1*fCharge2)>0.0) fSameSignAngle->Fill(&valueangle[0]);
2483 else fOppSignAngle->Fill(&valueangle[0]);
2484
2485 // Cut
2486 if(openingangle > fMaxopening3D) continue;
2487 //if(difftheta > fMaxopeningtheta) continue;
2488 //if(diffphi > fMaxopeningphi) continue;
2489
2490 // Invmass
2491 valuensparseDeltaPhiMaps[3] = invmass;
2492 if((fCharge1*fCharge2)>0.0) fSameSignDeltaPhiMaps->Fill(&valuensparseDeltaPhiMaps[0]);
2493 else fOppSignDeltaPhiMaps->Fill(&valuensparseDeltaPhiMaps[0]);
2494
2495 // Cut
2496 if(invmass < fMaxInvmass) {
2497 if((fCharge1*fCharge2)<0.0) oppositetaggedphotonic=kTRUE;
2498 if((fCharge1*fCharge2)>0.0) sametaggedphotonic=kTRUE;
2499 }
2500
2501
2502 }
2503 else
2504 {
2505 Int_t fPDGtrack1 = 11;
2506 Int_t fPDGtrack2 = 11;
2507
2508 Float_t fCharge1 = track1->Charge();
2509 Float_t fCharge2 = track2->Charge();
2510
2511 if(fCharge1>0) fPDGtrack1 = -11;
2512 if(fCharge2>0) fPDGtrack2 = -11;
2513
2514 AliKFParticle ktrack1(*track1, fPDGtrack1);
2515 AliKFParticle ktrack2(*track2, fPDGtrack2);
2516 AliKFParticle recoGamma(ktrack1, ktrack2);
2517
2518 //Reconstruction Cuts
2519 if(recoGamma.GetNDF()<1) continue;
2520 Double_t chi2OverNDF = recoGamma.GetChi2()/recoGamma.GetNDF();
2521 if(TMath::Sqrt(TMath::Abs(chi2OverNDF))>fChi2OverNDFCut) continue;
2522
2523 // DCA
2524 //Double_t dca12 = ktrack1.GetDistanceFromParticle(ktrack2);
2525 //if(dca12 > fMaxdca) continue;
2526
2527 // if set mass constraint
2528 if(fSetMassConstraint && pVtx) {
2529 AliKFVertex primV(*pVtx);
2530 primV += recoGamma;
2531 primV -= ktrack1;
2532 primV -= ktrack2;
2533 recoGamma.SetProductionVertex(primV);
2534 recoGamma.SetMassConstraint(0,0.0001);
2535 }
2536
2537 //Invariant Mass
2538 Double_t imass;
2539 Double_t width;
2540 recoGamma.GetMass(imass,width);
2541
2542 //Opening Angle (Total Angle)
2543 Double_t angle = ktrack1.GetAngle(ktrack2);
2544 valueangle[0] = angle;
2545 if((fCharge1*fCharge2)>0.0) fSameSignAngle->Fill(&valueangle[0]);
2546 else fOppSignAngle->Fill(&valueangle[0]);
2547
2548 // Cut
2549 if(angle > fMaxopening3D) continue;
2550
2551 // Invmass
2552 valuensparseDeltaPhiMaps[3] = imass;
2553 if((fCharge1*fCharge2)>0.0) fSameSignDeltaPhiMaps->Fill(&valuensparseDeltaPhiMaps[0]);
2554 else {
2555 fOppSignDeltaPhiMaps->Fill(&valuensparseDeltaPhiMaps[0]);
2556 /*
2557 if(valueangle[2] == kElectronfromconversionboth) {
2558 printf("Reconstructed charge1 %f, charge 2 %f and invmass %f",fCharge1,fCharge2,imass);
2559 printf("MC charge1 %d, charge 2 %d",pdg1,pdg2);
2560 printf("DCA %f",dca12);
2561 printf("Number of found %d",numberfound);
2562 }
2563 */
2564 }
2565
2566 // Cut
2567 if(imass < fMaxInvmass) {
2568 if((fCharge1*fCharge2)<0.0) oppositetaggedphotonic=kTRUE;
2569 if((fCharge1*fCharge2)>0.0) sametaggedphotonic=kTRUE;
2570 }
2571 }
2572 }
2573
2574 if(oppositetaggedphotonic && sametaggedphotonic){
2575 taggedphotonic = 6;
2576 }
2577
2578 if(!oppositetaggedphotonic && sametaggedphotonic){
2579 taggedphotonic = 4;
2580 }
2581
2582 if(oppositetaggedphotonic && !sametaggedphotonic){
2583 taggedphotonic = 2;
2584 }
2585
2586
2587 return taggedphotonic;
2588}
2589//_________________________________________________________________________
2590Int_t AliAnalysisTaskFlowTPCTOFEPSP::FindMother(Int_t tr, AliMCEvent *mcEvent, Int_t &indexmother){
2591 //
2592 // Find the mother if MC
2593 //
2594
2595 if(!mcEvent) return 0;
2596
2597 Int_t pdg = CheckPdg(tr,mcEvent);
2598 if(TMath::Abs(pdg)!= 11) {
2599 indexmother = -1;
2600 return kNoElectron;
2601 }
2602
2603 indexmother = IsMotherGamma(tr,mcEvent);
2604 if(indexmother > 0) return kElectronfromconversion;
2605 indexmother = IsMotherPi0(tr,mcEvent);
2606 if(indexmother > 0) return kElectronfrompi0;
2607 indexmother = IsMotherC(tr,mcEvent);
2608 if(indexmother > 0) return kElectronfromC;
2609 indexmother = IsMotherB(tr,mcEvent);
2610 if(indexmother > 0) return kElectronfromB;
2611 indexmother = IsMotherEta(tr,mcEvent);
2612 if(indexmother > 0) return kElectronfrometa;
2613
2614 return kElectronfromother;
2615
2616
2617}
2618//____________________________________________________________________________________________________________
2619Int_t AliAnalysisTaskFlowTPCTOFEPSP::CheckPdg(Int_t tr, AliMCEvent* mcEvent) {
2620
2621 //
2622 // Return the pdg of the particle
2623 //
2624
2625
2626 Int_t pdgcode = -1;
2627 if(tr < 0) return pdgcode;
2628
2629 if(!mcEvent) return pdgcode;
2630
2631 AliVParticle *mctrack = mcEvent->GetTrack(tr);
2632
2633
2634 if(mctrack->IsA() == AliMCParticle::Class()) {
2635 AliMCParticle *mctrackesd = NULL;
2636 if(!(mctrackesd = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return pdgcode;
2637 pdgcode = mctrackesd->PdgCode();
2638 }
2639
2640 if(mctrack->IsA() == AliAODMCParticle::Class()) {
2641 AliAODMCParticle *mctrackaod = NULL;
2642 if(!(mctrackaod = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return pdgcode;
2643 pdgcode = mctrackaod->GetPdgCode();
2644 }
2645
2646 return pdgcode;
2647
2648
2649}
2650//____________________________________________________________________________________________________________
2651Int_t AliAnalysisTaskFlowTPCTOFEPSP::IsMotherGamma(Int_t tr, AliMCEvent* mcEvent) {
2652
2653 //
2654 // Return the lab of gamma mother or -1 if not gamma
2655 //
2656
2657 if(tr < 0) return -1;
2658 AliVParticle *mctrack = mcEvent->GetTrack(tr);
2659
2660 if(mctrack->IsA() == AliMCParticle::Class()) {
2661 AliMCParticle *mctrackesd = NULL;
2662 if(!(mctrackesd = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2663 TParticle *particle = 0x0;
2664 particle = mctrackesd->Particle();
2665 // Take mother
2666 if(!particle) return -1;
2667 Int_t imother = particle->GetFirstMother();
2668 if(imother < 0) return -1;
2669 AliMCParticle *mothertrack = NULL;
2670 if(!(mothertrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2671 TParticle * mother = mothertrack->Particle();
2672 if(!mother) return -1;
2673 // Check gamma
2674 Int_t pdg = mother->GetPdgCode();
2675 if(TMath::Abs(pdg) == 22) return imother;
2676 if(TMath::Abs(pdg) == 11) {
2677 return IsMotherGamma(imother,mcEvent);
2678 }
2679 return -1;
2680 }
2681
2682 if(mctrack->IsA() == AliAODMCParticle::Class()) {
2683 AliAODMCParticle *mctrackaod = NULL;
2684 if(!(mctrackaod = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2685 // Take mother
2686 Int_t imother = mctrackaod->GetMother();
2687 if(imother < 0) return -1;
2688 AliAODMCParticle *mothertrack = NULL;
2689 if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2690 // Check gamma
2691 Int_t pdg = mothertrack->GetPdgCode();
2692 if(TMath::Abs(pdg) == 22) return imother;
2693 if(TMath::Abs(pdg) == 11) {
2694 return IsMotherGamma(imother,mcEvent);
2695 }
2696 return -1;
2697
2698 }
2699
2700 return -1;
2701
2702
2703}
2704//
2705//____________________________________________________________________________________________________________
2706Int_t AliAnalysisTaskFlowTPCTOFEPSP::IsMotherPi0(Int_t tr, AliMCEvent* mcEvent) {
2707
2708 //
2709 // Return the lab of pi0 mother or -1 if not pi0
2710 //
2711
2712 if(tr < 0) return -1;
2713 AliVParticle *mctrack = mcEvent->GetTrack(tr);
2714
2715 if(mctrack->IsA() == AliMCParticle::Class()) {
2716 AliMCParticle *mctrackesd = NULL;
2717 if(!(mctrackesd = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2718 TParticle *particle = 0x0;
2719 particle = mctrackesd->Particle();
2720 // Take mother
2721 if(!particle) return -1;
2722 Int_t imother = particle->GetFirstMother();
2723 if(imother < 0) return -1;
2724 AliMCParticle *mothertrack = NULL;
2725 if(!(mothertrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2726 TParticle * mother = mothertrack->Particle();
2727 if(!mother) return -1;
2728 // Check gamma
2729 Int_t pdg = mother->GetPdgCode();
2730 if(TMath::Abs(pdg) == 111) return imother;
2731 if(TMath::Abs(pdg) == 11) {
2732 return IsMotherPi0(imother,mcEvent);
2733 }
2734 return -1;
2735 }
2736
2737 if(mctrack->IsA() == AliAODMCParticle::Class()) {
2738 AliAODMCParticle *mctrackaod = NULL;
2739 if(!(mctrackaod = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2740 // Take mother
2741 Int_t imother = mctrackaod->GetMother();
2742 if(imother < 0) return -1;
2743 AliAODMCParticle *mothertrack = NULL;
2744 if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2745 // Check gamma
2746 Int_t pdg = mothertrack->GetPdgCode();
2747 if(TMath::Abs(pdg) == 111) return imother;
2748 if(TMath::Abs(pdg) == 11) {
2749 return IsMotherPi0(imother,mcEvent);
2750 }
2751 return -1;
2752 }
2753
2754 return -1;
2755
2756}
2757//____________________________________________________________________________________________________________
2758Int_t AliAnalysisTaskFlowTPCTOFEPSP::IsMotherC(Int_t tr, AliMCEvent* mcEvent) {
2759
2760 //
2761 // Return the lab of signal mother or -1 if not signal
2762 //
2763
2764 if(tr < 0) return -1;
2765 AliVParticle *mctrack = mcEvent->GetTrack(tr);
2766
2767 if(mctrack->IsA() == AliMCParticle::Class()) {
2768 AliMCParticle *mctrackesd = NULL;
2769 if(!(mctrackesd = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2770 TParticle *particle = 0x0;
2771 particle = mctrackesd->Particle();
2772 // Take mother
2773 if(!particle) return -1;
2774 Int_t imother = particle->GetFirstMother();
2775 if(imother < 0) return -1;
2776 AliMCParticle *mothertrack = NULL;
2777 if(!(mothertrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2778 TParticle * mother = mothertrack->Particle();
2779 if(!mother) return -1;
2780 // Check gamma
2781 Int_t pdg = mother->GetPdgCode();
2782 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;
2783 if(TMath::Abs(pdg) == 11) {
2784 return IsMotherC(imother,mcEvent);
2785 }
2786 return -1;
2787 }
2788
2789 if(mctrack->IsA() == AliAODMCParticle::Class()) {
2790 AliAODMCParticle *mctrackaod = NULL;
2791 if(!(mctrackaod = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2792 // Take mother
2793 Int_t imother = mctrackaod->GetMother();
2794 if(imother < 0) return -1;
2795 AliAODMCParticle *mothertrack = NULL;
2796 if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2797 // Check gamma
2798 Int_t pdg = mothertrack->GetPdgCode();
2799 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;
2800 if(TMath::Abs(pdg) == 11) {
2801 return IsMotherC(imother,mcEvent);
2802 }
2803 return -1;
2804 }
2805
2806 return -1;
2807
2808}
2809//____________________________________________________________________________________________________________
2810Int_t AliAnalysisTaskFlowTPCTOFEPSP::IsMotherB(Int_t tr, AliMCEvent* mcEvent) {
2811
2812 //
2813 // Return the lab of signal mother or -1 if not signal
2814 //
2815
2816 if(tr < 0) return -1;
2817 AliVParticle *mctrack = mcEvent->GetTrack(tr);
2818
2819 if(mctrack->IsA() == AliMCParticle::Class()) {
2820 AliMCParticle *mctrackesd = NULL;
2821 if(!(mctrackesd = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2822 TParticle *particle = 0x0;
2823 particle = mctrackesd->Particle();
2824 // Take mother
2825 if(!particle) return -1;
2826 Int_t imother = particle->GetFirstMother();
2827 if(imother < 0) return -1;
2828 AliMCParticle *mothertrack = NULL;
2829 if(!(mothertrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2830 TParticle * mother = mothertrack->Particle();
2831 if(!mother) return -1;
2832 // Check gamma
2833 Int_t pdg = mother->GetPdgCode();
2834 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;
2835 if(TMath::Abs(pdg) == 11) {
2836 return IsMotherB(imother,mcEvent);
2837 }
2838 return -1;
2839 }
2840
2841 if(mctrack->IsA() == AliAODMCParticle::Class()) {
2842 AliAODMCParticle *mctrackaod = NULL;
2843 if(!(mctrackaod = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2844 // Take mother
2845 Int_t imother = mctrackaod->GetMother();
2846 if(imother < 0) return -1;
2847 AliAODMCParticle *mothertrack = NULL;
2848 if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2849 // Check gamma
2850 Int_t pdg = mothertrack->GetPdgCode();
2851 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;
2852 if(TMath::Abs(pdg) == 11) {
2853 return IsMotherB(imother,mcEvent);
2854 }
2855 return -1;
2856 }
2857
2858 return -1;
2859
2860}
2861//____________________________________________________________________________________________________________
2862Int_t AliAnalysisTaskFlowTPCTOFEPSP::IsMotherEta(Int_t tr, AliMCEvent* mcEvent) {
2863
2864 //
2865 // Return the lab of pi0 mother or -1 if not pi0
2866 //
2867
2868 if(tr < 0) return -1;
2869 AliVParticle *mctrack = mcEvent->GetTrack(tr);
2870
2871 if(mctrack->IsA() == AliMCParticle::Class()) {
2872 AliMCParticle *mctrackesd = NULL;
2873 if(!(mctrackesd = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2874 TParticle *particle = 0x0;
2875 particle = mctrackesd->Particle();
2876 // Take mother
2877 if(!particle) return -1;
2878 Int_t imother = particle->GetFirstMother();
2879 if(imother < 0) return -1;
2880 AliMCParticle *mothertrack = NULL;
2881 if(!(mothertrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2882 TParticle * mother = mothertrack->Particle();
2883 if(!mother) return -1;
2884 // Check gamma
2885 Int_t pdg = mother->GetPdgCode();
2886 if(TMath::Abs(pdg) == 221) return imother;
2887 if(TMath::Abs(pdg) == 11) {
2888 return IsMotherEta(imother,mcEvent);
2889 }
2890 return -1;
2891 }
2892
2893 if(mctrack->IsA() == AliAODMCParticle::Class()) {
2894 AliAODMCParticle *mctrackaod = NULL;
2895 if(!(mctrackaod = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
2896 // Take mother
2897 Int_t imother = mctrackaod->GetMother();
2898 if(imother < 0) return -1;
2899 AliAODMCParticle *mothertrack = NULL;
2900 if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
2901 // Check gamma
2902 Int_t pdg = mothertrack->GetPdgCode();
2903 if(TMath::Abs(pdg) == 221) return imother;
2904 if(TMath::Abs(pdg) == 11) {
2905 return IsMotherEta(imother,mcEvent);
2906 }
2907 return -1;
2908 }
2909
2910 return -1;
2911
2912}