1 /**************************************************************************
\r
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
\r
4 * Author: The ALICE Off-line Project. *
\r
5 * Contributors are mentioned in the code where appropriate. *
\r
7 * Permission to use, copy, modify and distribute this software and its *
\r
8 * documentation strictly for non-commercial purposes is hereby granted *
\r
9 * without fee, provided that the above copyright notice appears in all *
\r
10 * copies and that both the copyright notice and this permission notice *
\r
11 * appear in the supporting documentation. The authors make no claims *
\r
12 * about the suitability of this software for any purpose. It is *
\r
13 * provided "as is" without express or implied warranty. *
\r
14 **************************************************************************/
\r
19 // Raphaelle Bailhache <R.Bailhache@gsi.de>
\r
25 #include "TVector2.h"
\r
26 #include "THnSparse.h"
\r
28 #include "TRandom3.h"
\r
29 #include "TProfile.h"
\r
30 #include "TProfile2D.h"
\r
31 #include "TLorentzVector.h"
\r
32 #include "TParticle.h"
\r
34 #include "AliVEventHandler.h"
\r
35 #include "AliAnalysisTaskSE.h"
\r
36 #include "AliAnalysisManager.h"
\r
38 #include "AliVEvent.h"
\r
39 #include "AliESDInputHandler.h"
\r
40 #include "AliMCEvent.h"
\r
42 #include "AliESDEvent.h"
\r
44 #include "AliPIDResponse.h"
\r
45 #include "AliESDVZERO.h"
\r
46 #include "AliESDUtils.h"
\r
47 #include "AliMCParticle.h"
\r
48 #include "AliAODMCParticle.h"
\r
49 #include "AliAODEvent.h"
\r
50 #include "AliAODVertex.h"
\r
51 #include "AliAODTrack.h"
\r
52 #include "AliVTrack.h"
\r
53 #include "AliESDtrack.h"
\r
54 #include "AliESDtrackCuts.h"
\r
55 #include "AliAODTrack.h"
\r
56 #include "AliStack.h"
\r
57 #include "AliMCEvent.h"
\r
59 #include "AliFlowCandidateTrack.h"
\r
60 #include "AliFlowEvent.h"
\r
61 #include "AliFlowTrackCuts.h"
\r
62 #include "AliFlowVector.h"
\r
63 #include "AliFlowCommonConstants.h"
\r
64 #include "AliKFParticle.h"
\r
65 #include "AliKFVertex.h"
\r
67 #include "AliHFEcuts.h"
\r
68 #include "AliHFEpid.h"
\r
69 #include "AliHFEpidQAmanager.h"
\r
70 #include "AliHFEtools.h"
\r
71 #include "AliHFEVZEROEventPlane.h"
\r
73 #include "AliCentrality.h"
\r
74 #include "AliEventplane.h"
\r
75 #include "AliAnalysisTaskHFEFlow.h"
\r
78 //____________________________________________________________________
\r
79 AliAnalysisTaskHFEFlow::AliAnalysisTaskHFEFlow() :
\r
80 AliAnalysisTaskSE(),
\r
82 fAODAnalysis(kFALSE),
\r
83 fUseFlagAOD(kFALSE),
\r
86 fVZEROEventPlane(kFALSE),
\r
87 fVZEROEventPlaneA(kFALSE),
\r
88 fVZEROEventPlaneC(kFALSE),
\r
89 fSubEtaGapTPC(kFALSE),
\r
91 fNbBinsCentralityQCumulant(4),
\r
92 fNbBinsPtQCumulant(12),
\r
93 fMinPtQCumulant(0.2),
\r
94 fMaxPtQCumulant(6.0),
\r
95 fAfterBurnerOn(kFALSE),
\r
96 fNonFlowNumberOfTrackClones(0),
\r
102 fMaxNumberOfIterations(100),
\r
103 fPrecisionPhi(0.001),
\r
104 fUseMCReactionPlane(kFALSE),
\r
107 fChi2OverNDFCut(3.0),
\r
109 fMaxopeningtheta(0.02),
\r
110 fMaxopeningphi(0.1),
\r
111 fMaxopening3D(0.1),
\r
113 fSetMassConstraint(kFALSE),
\r
122 fHFEBackgroundCuts(0),
\r
124 fPIDBackgroundqa(0),
\r
125 fAlgorithmMA(kTRUE),
\r
127 fCounterPoolBackground(0),
\r
128 fHFEVZEROEventPlane(0x0),
\r
131 fEventPlaneaftersubtraction(0x0),
\r
132 fCosSin2phiep(0x0),
\r
137 fSin2phiephiep(0x0),
\r
140 fProfileCosResab(0x0),
\r
141 fProfileCosResac(0x0),
\r
142 fProfileCosResbc(0x0),
\r
145 fProfileCosRes(0x0),
\r
146 fTrackingCuts(0x0),
\r
147 fDeltaPhiMapsBeforePID(0x0),
\r
148 fCosPhiMapsBeforePID(0x0),
\r
149 fDeltaPhiMaps(0x0),
\r
150 fDeltaPhiMapsContamination(0x0),
\r
152 fProfileCosPhiMaps(0x0),
\r
153 fDeltaPhiMapsTaggedPhotonic(0x0),
\r
154 fCosPhiMapsTaggedPhotonic(0x0),
\r
155 fDeltaPhiMapsTaggedNonPhotonic(0x0),
\r
156 fCosPhiMapsTaggedNonPhotonic(0x0),
\r
157 fDeltaPhiMapsTaggedPhotonicLS(0x0),
\r
158 fCosPhiMapsTaggedPhotonicLS(0x0),
\r
159 fMCSourceDeltaPhiMaps(0x0),
\r
160 fOppSignDeltaPhiMaps(0x0),
\r
161 fSameSignDeltaPhiMaps(0x0),
\r
162 fOppSignAngle(0x0),
\r
163 fSameSignAngle(0x0)
\r
167 for(Int_t k = 0; k < 10; k++) {
\r
168 fBinCentralityLess[k] = 0.0;
\r
172 //______________________________________________________________________________
\r
173 AliAnalysisTaskHFEFlow:: AliAnalysisTaskHFEFlow(const char *name) :
\r
174 AliAnalysisTaskSE(name),
\r
176 fAODAnalysis(kFALSE),
\r
177 fUseFlagAOD(kFALSE),
\r
180 fVZEROEventPlane(kFALSE),
\r
181 fVZEROEventPlaneA(kFALSE),
\r
182 fVZEROEventPlaneC(kFALSE),
\r
183 fSubEtaGapTPC(kFALSE),
\r
185 fNbBinsCentralityQCumulant(4),
\r
186 fNbBinsPtQCumulant(15),
\r
187 fMinPtQCumulant(0.0),
\r
188 fMaxPtQCumulant(6.0),
\r
189 fAfterBurnerOn(kFALSE),
\r
190 fNonFlowNumberOfTrackClones(0),
\r
196 fMaxNumberOfIterations(100),
\r
197 fPrecisionPhi(0.001),
\r
198 fUseMCReactionPlane(kFALSE),
\r
201 fChi2OverNDFCut(3.0),
\r
203 fMaxopeningtheta(0.02),
\r
204 fMaxopeningphi(0.1),
\r
205 fMaxopening3D(0.1),
\r
207 fSetMassConstraint(kFALSE),
\r
216 fHFEBackgroundCuts(0),
\r
218 fPIDBackgroundqa(0),
\r
219 fAlgorithmMA(kTRUE),
\r
221 fCounterPoolBackground(0),
\r
222 fHFEVZEROEventPlane(0x0),
\r
225 fEventPlaneaftersubtraction(0x0),
\r
226 fCosSin2phiep(0x0),
\r
231 fSin2phiephiep(0x0),
\r
234 fProfileCosResab(0x0),
\r
235 fProfileCosResac(0x0),
\r
236 fProfileCosResbc(0x0),
\r
239 fProfileCosRes(0x0),
\r
240 fTrackingCuts(0x0),
\r
241 fDeltaPhiMapsBeforePID(0x0),
\r
242 fCosPhiMapsBeforePID(0x0),
\r
243 fDeltaPhiMaps(0x0),
\r
244 fDeltaPhiMapsContamination(0x0),
\r
246 fProfileCosPhiMaps(0x0),
\r
247 fDeltaPhiMapsTaggedPhotonic(0x0),
\r
248 fCosPhiMapsTaggedPhotonic(0x0),
\r
249 fDeltaPhiMapsTaggedNonPhotonic(0x0),
\r
250 fCosPhiMapsTaggedNonPhotonic(0x0),
\r
251 fDeltaPhiMapsTaggedPhotonicLS(0x0),
\r
252 fCosPhiMapsTaggedPhotonicLS(0x0),
\r
253 fMCSourceDeltaPhiMaps(0x0),
\r
254 fOppSignDeltaPhiMaps(0x0),
\r
255 fSameSignDeltaPhiMaps(0x0),
\r
256 fOppSignAngle(0x0),
\r
257 fSameSignAngle(0x0)
\r
263 for(Int_t k = 0; k < 10; k++) {
\r
264 fBinCentralityLess[k] = 0.0;
\r
266 fBinCentralityLess[0] = 0.0;
\r
267 fBinCentralityLess[1] = 20.0;
\r
268 fBinCentralityLess[2] = 40.0;
\r
269 fBinCentralityLess[3] = 60.0;
\r
270 fBinCentralityLess[4] = 80.0;
\r
272 fPID = new AliHFEpid("hfePid");
\r
273 fPIDqa = new AliHFEpidQAmanager;
\r
275 fPIDBackground = new AliHFEpid("hfePidBackground");
\r
276 fPIDBackgroundqa = new AliHFEpidQAmanager;
\r
278 fPIDTOFOnly = new AliHFEpid("hfePidTOFOnly");
\r
280 DefineInput(0,TChain::Class());
\r
281 DefineOutput(1, TList::Class());
\r
282 for(Int_t bincless = 0; bincless < fNbBinsCentralityQCumulant; bincless++) {
\r
283 DefineOutput(bincless+2,AliFlowEventSimple::Class());
\r
287 //____________________________________________________________
\r
288 AliAnalysisTaskHFEFlow::AliAnalysisTaskHFEFlow(const AliAnalysisTaskHFEFlow &ref):
\r
289 AliAnalysisTaskSE(ref),
\r
291 fAODAnalysis(ref.fAODAnalysis),
\r
292 fUseFlagAOD(ref.fUseFlagAOD),
\r
293 fApplyCut(ref.fApplyCut),
\r
294 fFlags(ref.fFlags),
\r
295 fVZEROEventPlane(ref.fVZEROEventPlane),
\r
296 fVZEROEventPlaneA(ref.fVZEROEventPlaneA),
\r
297 fVZEROEventPlaneC(ref.fVZEROEventPlaneC),
\r
298 fSubEtaGapTPC(ref.fSubEtaGapTPC),
\r
299 fEtaGap(ref.fEtaGap),
\r
300 fNbBinsCentralityQCumulant(ref.fNbBinsCentralityQCumulant),
\r
301 fNbBinsPtQCumulant(ref.fNbBinsPtQCumulant),
\r
302 fMinPtQCumulant(ref.fMinPtQCumulant),
\r
303 fMaxPtQCumulant(ref.fMaxPtQCumulant),
\r
304 fAfterBurnerOn(ref.fAfterBurnerOn),
\r
305 fNonFlowNumberOfTrackClones(ref.fNonFlowNumberOfTrackClones),
\r
311 fMaxNumberOfIterations(ref.fMaxNumberOfIterations),
\r
312 fPrecisionPhi(ref.fPrecisionPhi),
\r
313 fUseMCReactionPlane(ref.fUseMCReactionPlane),
\r
314 fMCPID(ref.fMCPID),
\r
315 fNoPID(ref.fNoPID),
\r
316 fChi2OverNDFCut(ref.fChi2OverNDFCut),
\r
317 fMaxdca(ref.fMaxdca),
\r
318 fMaxopeningtheta(ref.fMaxopeningtheta),
\r
319 fMaxopeningphi(ref.fMaxopeningphi),
\r
320 fMaxopening3D(ref.fMaxopening3D),
\r
321 fMaxInvmass(ref.fMaxInvmass),
\r
322 fSetMassConstraint(ref.fSetMassConstraint),
\r
323 fDebugLevel(ref.fDebugLevel),
\r
331 fHFEBackgroundCuts(NULL),
\r
332 fPIDBackground(NULL),
\r
333 fPIDBackgroundqa(NULL),
\r
334 fAlgorithmMA(ref.fAlgorithmMA),
\r
336 fCounterPoolBackground(ref.fCounterPoolBackground),
\r
337 fHFEVZEROEventPlane(NULL),
\r
340 fEventPlaneaftersubtraction(NULL),
\r
341 fCosSin2phiep(NULL),
\r
346 fSin2phiephiep(NULL),
\r
349 fProfileCosResab(NULL),
\r
350 fProfileCosResac(NULL),
\r
351 fProfileCosResbc(NULL),
\r
354 fProfileCosRes(NULL),
\r
355 fTrackingCuts(NULL),
\r
356 fDeltaPhiMapsBeforePID(NULL),
\r
357 fCosPhiMapsBeforePID(NULL),
\r
358 fDeltaPhiMaps(NULL),
\r
359 fDeltaPhiMapsContamination(NULL),
\r
361 fProfileCosPhiMaps(NULL),
\r
362 fDeltaPhiMapsTaggedPhotonic(NULL),
\r
363 fCosPhiMapsTaggedPhotonic(NULL),
\r
364 fDeltaPhiMapsTaggedNonPhotonic(NULL),
\r
365 fCosPhiMapsTaggedNonPhotonic(NULL),
\r
366 fDeltaPhiMapsTaggedPhotonicLS(NULL),
\r
367 fCosPhiMapsTaggedPhotonicLS(NULL),
\r
368 fMCSourceDeltaPhiMaps(NULL),
\r
369 fOppSignDeltaPhiMaps(NULL),
\r
370 fSameSignDeltaPhiMaps(NULL),
\r
371 fOppSignAngle(NULL),
\r
372 fSameSignAngle(NULL)
\r
375 // Copy Constructor
\r
380 //____________________________________________________________
\r
381 AliAnalysisTaskHFEFlow &AliAnalysisTaskHFEFlow::operator=(const AliAnalysisTaskHFEFlow &ref){
\r
383 // Assignment operator
\r
390 //____________________________________________________________
\r
391 void AliAnalysisTaskHFEFlow::Copy(TObject &o) const {
\r
393 // Copy into object o
\r
395 AliAnalysisTaskHFEFlow &target = dynamic_cast<AliAnalysisTaskHFEFlow &>(o);
\r
396 target.fAODAnalysis = fAODAnalysis;
\r
397 target.fUseFlagAOD = fUseFlagAOD;
\r
398 target.fApplyCut = fApplyCut;
\r
399 target.fFlags = fFlags;
\r
400 target.fVZEROEventPlane = fVZEROEventPlane;
\r
401 target.fVZEROEventPlaneA = fVZEROEventPlaneA;
\r
402 target.fVZEROEventPlaneC = fVZEROEventPlaneC;
\r
403 target.fSubEtaGapTPC = fSubEtaGapTPC;
\r
404 target.fEtaGap = fEtaGap;
\r
405 target.fNbBinsCentralityQCumulant = fNbBinsCentralityQCumulant;
\r
406 target.fNbBinsPtQCumulant = fNbBinsPtQCumulant;
\r
407 target.fMinPtQCumulant = fMinPtQCumulant;
\r
408 target.fMaxPtQCumulant = fMaxPtQCumulant;
\r
409 target.fAfterBurnerOn = fAfterBurnerOn;
\r
410 target.fNonFlowNumberOfTrackClones = fNonFlowNumberOfTrackClones;
\r
416 target.fMaxNumberOfIterations = fMaxNumberOfIterations;
\r
417 target.fPrecisionPhi = fPrecisionPhi;
\r
418 target.fUseMCReactionPlane = fUseMCReactionPlane;
\r
419 target.fMCPID = fMCPID;
\r
420 target.fNoPID = fNoPID;
\r
421 target.fChi2OverNDFCut = fChi2OverNDFCut;
\r
422 target.fMaxdca = fMaxdca;
\r
423 target.fMaxopeningtheta = fMaxopeningtheta;
\r
424 target.fMaxopeningphi = fMaxopeningphi;
\r
425 target.fMaxopening3D = fMaxopening3D;
\r
426 target.fMaxInvmass = fMaxInvmass;
\r
427 target.fSetMassConstraint = fSetMassConstraint;
\r
428 target.fAlgorithmMA = fAlgorithmMA;
\r
429 target.fCounterPoolBackground = fCounterPoolBackground;
\r
430 target.fDebugLevel = fDebugLevel;
\r
431 target.fcutsRP = fcutsRP;
\r
432 target.fcutsPOI = fcutsPOI;
\r
433 target.fHFECuts = fHFECuts;
\r
434 target.fPID = fPID;
\r
435 target.fPIDqa = fPIDqa;
\r
436 target.fHFEVZEROEventPlane = fHFEVZEROEventPlane;
\r
439 //____________________________________________________________
\r
440 AliAnalysisTaskHFEFlow::~AliAnalysisTaskHFEFlow(){
\r
444 if(fArraytrack) delete fArraytrack;
\r
445 if(fListHist) delete fListHist;
\r
446 if(fcutsRP) delete fcutsRP;
\r
447 if(fcutsPOI) delete fcutsPOI;
\r
448 if(fHFECuts) delete fHFECuts;
\r
449 if(fPID) delete fPID;
\r
450 if(fPIDTOFOnly) delete fPIDTOFOnly;
\r
451 if(fPIDqa) delete fPIDqa;
\r
452 if(fflowEvent) delete fflowEvent;
\r
453 if(fHFEBackgroundCuts) delete fHFEBackgroundCuts;
\r
454 if(fPIDBackground) delete fPIDBackground;
\r
455 if(fPIDBackgroundqa) delete fPIDBackgroundqa;
\r
456 //if(fHFEVZEROEventPlane) delete fHFEVZEROEventPlane;
\r
460 //________________________________________________________________________
\r
461 void AliAnalysisTaskHFEFlow::UserCreateOutputObjects()
\r
464 //********************
\r
465 // Create histograms
\r
466 //********************
\r
472 //---------Data selection----------
\r
473 //kMC, kGlobal, kTPCstandalone, kSPDtracklet, kPMD
\r
474 //AliFlowTrackCuts::trackParameterType rptype = AliFlowTrackCuts::kGlobal;
\r
475 //AliFlowTrackCuts::trackParameterType poitype = AliFlowTrackCuts::kGlobal;
\r
477 //---------Parameter mixing--------
\r
478 //kPure - no mixing, kTrackWithMCkine, kTrackWithMCPID, kTrackWithMCpt
\r
479 //AliFlowTrackCuts::trackParameterMix rpmix = AliFlowTrackCuts::kPure;
\r
480 //AliFlowTrackCuts::trackParameterMix poimix = AliFlowTrackCuts::kPure;
\r
483 AliVEventHandler *inputHandler = dynamic_cast<AliVEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
\r
484 if(!TString(inputHandler->IsA()->GetName()).CompareTo("AliAODInputHandler")){
\r
485 SetAODAnalysis(kTRUE);
\r
486 //printf("Put AOD analysis on\n");
\r
488 SetAODAnalysis(kFALSE);
\r
493 fcutsRP = AliFlowTrackCuts::GetStandardTPCStandaloneTrackCuts2010();
\r
494 fcutsRP->SetName("StandartTPC");
\r
495 fcutsRP->SetEtaRange(-0.9,0.9);
\r
496 fcutsRP->SetQA(kTRUE);
\r
497 //TList *qaCutsRP = fcutsRP->GetQA();
\r
498 //qaCutsRP->SetName("QA_StandartTPC_RP");
\r
501 fcutsPOI = new AliFlowTrackCuts("dummy");
\r
502 fcutsPOI->SetParamType(AliFlowTrackCuts::kGlobal);
\r
503 fcutsPOI->SetPtRange(+1,-1); // select nothing QUICK
\r
504 fcutsPOI->SetEtaRange(+1,-1); // select nothing VZERO
\r
507 AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
\r
508 cc->SetNbinsMult(10000);
\r
510 cc->SetMultMax(10000.);
\r
511 cc->SetNbinsPt(fNbBinsPtQCumulant);
\r
512 cc->SetPtMin(fMinPtQCumulant);
\r
513 cc->SetPtMax(fMaxPtQCumulant);
\r
514 cc->SetNbinsPhi(180);
\r
515 cc->SetPhiMin(0.0);
\r
516 cc->SetPhiMax(TMath::TwoPi());
\r
517 cc->SetNbinsEta(200);
\r
518 cc->SetEtaMin(-0.9);
\r
519 cc->SetEtaMax(+0.9);
\r
520 cc->SetNbinsQ(500);
\r
528 fHFECuts = new AliHFEcuts;
\r
529 fHFECuts->CreateStandardCuts();
\r
531 fHFECuts->Initialize();
\r
532 if(fAODAnalysis) fHFECuts->SetAOD();
\r
535 //fPID->SetHasMCData(HasMCData());
\r
536 if(!fPID->GetNumberOfPIDdetectors()) fPID->AddDetector("TPC", 0);
\r
537 fPID->InitializePID();
\r
538 fPIDqa->Initialize(fPID);
\r
539 fPID->SortDetectors();
\r
541 if(!fPIDTOFOnly->GetNumberOfPIDdetectors()) fPIDTOFOnly->AddDetector("TPC", 0);
\r
542 fPIDTOFOnly->InitializePID();
\r
543 fPIDTOFOnly->SortDetectors();
\r
545 // HFE Background cuts
\r
547 if(!fHFEBackgroundCuts){
\r
548 fHFEBackgroundCuts = new AliESDtrackCuts();
\r
549 fHFEBackgroundCuts->SetName("nackgroundcuts");
\r
550 //Configure Default Track Cuts
\r
551 fHFEBackgroundCuts->SetAcceptKinkDaughters(kFALSE);
\r
552 fHFEBackgroundCuts->SetRequireTPCRefit(kTRUE);
\r
553 fHFEBackgroundCuts->SetEtaRange(-0.9,0.9);
\r
554 fHFEBackgroundCuts->SetRequireSigmaToVertex(kTRUE);
\r
555 fHFEBackgroundCuts->SetMaxChi2PerClusterTPC(4.0);
\r
556 fHFEBackgroundCuts->SetMinNClustersTPC(50);
\r
557 fHFEBackgroundCuts->SetPtRange(0.3,1e10);
\r
560 // PID background HFE
\r
561 if(!fPIDBackground->GetNumberOfPIDdetectors()) fPIDBackground->AddDetector("TPC", 0);
\r
562 fPIDBackground->InitializePID();
\r
563 fPIDBackgroundqa->Initialize(fPIDBackground);
\r
564 fPIDBackground->SortDetectors();
\r
568 //**************************
\r
569 // Bins for the THnSparse
\r
570 //**************************
\r
573 Int_t nBinsPt = 44;
\r
574 Double_t minPt = 0.1;
\r
575 Double_t maxPt = 20.0;
\r
576 Double_t binLimLogPt[nBinsPt+1];
\r
577 Double_t binLimPt[nBinsPt+1];
\r
578 for(Int_t i=0; i<=nBinsPt; i++) binLimLogPt[i]=(Double_t)TMath::Log10(minPt) + (TMath::Log10(maxPt)-TMath::Log10(minPt))/nBinsPt*(Double_t)i ;
\r
579 for(Int_t i=0; i<=nBinsPt; i++) binLimPt[i]=(Double_t)TMath::Power(10,binLimLogPt[i]);
\r
582 Int_t nBinsPt = 35;
\r
583 Double_t binLimPt[36] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2,
\r
584 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5.,
\r
585 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};
\r
588 Int_t nBinsPtPlus = fNbBinsPtQCumulant;
\r
589 Double_t minPtPlus = fMinPtQCumulant;
\r
590 Double_t maxPtPlus = fMaxPtQCumulant;
\r
591 Double_t binLimPtPlus[nBinsPtPlus+1];
\r
592 for(Int_t i=0; i<=nBinsPtPlus; i++) binLimPtPlus[i]=(Double_t)minPtPlus + (maxPtPlus-minPtPlus)/nBinsPtPlus*(Double_t)i ;
\r
594 Int_t nBinsEta = 8;
\r
595 Double_t minEta = -0.8;
\r
596 Double_t maxEta = 0.8;
\r
597 Double_t binLimEta[nBinsEta+1];
\r
598 for(Int_t i=0; i<=nBinsEta; i++) binLimEta[i]=(Double_t)minEta + (maxEta-minEta)/nBinsEta*(Double_t)i ;
\r
600 Int_t nBinsStep = 7;
\r
601 Double_t minStep = 0.;
\r
602 Double_t maxStep = 7.;
\r
603 Double_t binLimStep[nBinsStep+1];
\r
604 for(Int_t i=0; i<=nBinsStep; i++) binLimStep[i]=(Double_t)minStep + (maxStep-minStep)/nBinsStep*(Double_t)i ;
\r
606 Int_t nBinsEtaLess = 2;
\r
607 Double_t binLimEtaLess[nBinsEtaLess+1];
\r
608 for(Int_t i=0; i<=nBinsEtaLess; i++) binLimEtaLess[i]=(Double_t)minEta + (maxEta-minEta)/nBinsEtaLess*(Double_t)i ;
\r
610 Int_t nBinsCos = 50;
\r
611 Double_t minCos = -1.0;
\r
612 Double_t maxCos = 1.0;
\r
613 Double_t binLimCos[nBinsCos+1];
\r
614 for(Int_t i=0; i<=nBinsCos; i++) binLimCos[i]=(Double_t)minCos + (maxCos-minCos)/nBinsCos*(Double_t)i ;
\r
617 Double_t minC = 0.0;
\r
618 Double_t maxC = 11.0;
\r
619 Double_t binLimC[nBinsC+1];
\r
620 for(Int_t i=0; i<=nBinsC; i++) binLimC[i]=(Double_t)minC + (maxC-minC)/nBinsC*(Double_t)i ;
\r
622 Int_t nBinsCMore = 20;
\r
623 Double_t minCMore = 0.0;
\r
624 Double_t maxCMore = 20.0;
\r
625 Double_t binLimCMore[nBinsCMore+1];
\r
626 for(Int_t i=0; i<=nBinsCMore; i++) binLimCMore[i]=(Double_t)minCMore + (maxCMore-minCMore)/nBinsCMore*(Double_t)i ;
\r
628 Int_t nBinsPhi = 12;
\r
629 Double_t minPhi = 0.0;
\r
630 Double_t maxPhi = TMath::Pi();
\r
631 Double_t binLimPhi[nBinsPhi+1];
\r
632 for(Int_t i=0; i<=nBinsPhi; i++) {
\r
633 binLimPhi[i]=(Double_t)minPhi + (maxPhi-minPhi)/nBinsPhi*(Double_t)i ;
\r
634 //printf("bin phi is %f for %d\n",binLimPhi[i],i);
\r
637 Int_t nBinsPhiLess = 2.0;
\r
638 Double_t minPhiLess = 0.0;
\r
639 Double_t maxPhiLess = 2.0;
\r
640 Double_t binLimPhiLess[nBinsPhiLess+1];
\r
641 for(Int_t i=0; i<=nBinsPhiLess; i++) {
\r
642 binLimPhiLess[i]=(Double_t)minPhiLess + (maxPhiLess-minPhiLess)/nBinsPhiLess*(Double_t)i ;
\r
645 Int_t nBinsTPCdEdx = 140;
\r
646 Double_t minTPCdEdx = -12.0;
\r
647 Double_t maxTPCdEdx = 12.0;
\r
648 Double_t binLimTPCdEdx[nBinsTPCdEdx+1];
\r
649 for(Int_t i=0; i<=nBinsTPCdEdx; i++) {
\r
650 binLimTPCdEdx[i]=(Double_t)minTPCdEdx + (maxTPCdEdx-minTPCdEdx)/nBinsTPCdEdx*(Double_t)i ;
\r
653 Int_t nBinsAngle = 40;
\r
654 Double_t minAngle = 0.0;
\r
655 Double_t maxAngle = 1.0;
\r
656 Double_t binLimAngle[nBinsAngle+1];
\r
657 for(Int_t i=0; i<=nBinsAngle; i++) {
\r
658 binLimAngle[i]=(Double_t)minAngle + (maxAngle-minAngle)/nBinsAngle*(Double_t)i ;
\r
659 //printf("bin phi is %f for %d\n",binLimPhi[i],i);
\r
662 Int_t nBinsCharge = 2;
\r
663 Double_t minCharge = -1.0;
\r
664 Double_t maxCharge = 1.0;
\r
665 Double_t binLimCharge[nBinsCharge+1];
\r
666 for(Int_t i=0; i<=nBinsCharge; i++) binLimCharge[i]=(Double_t)minCharge + (maxCharge-minCharge)/nBinsCharge*(Double_t)i ;
\r
668 Int_t nBinsSource = 10;
\r
669 Double_t minSource = 0.;
\r
670 Double_t maxSource = 10.;
\r
671 Double_t binLimSource[nBinsSource+1];
\r
672 for(Int_t i=0; i<=nBinsSource; i++) binLimSource[i]=(Double_t)minSource + (maxSource-minSource)/nBinsSource*(Double_t)i ;
\r
674 Int_t nBinsInvMass = 50;
\r
675 Double_t minInvMass = 0.;
\r
676 Double_t maxInvMass = 0.3;
\r
677 Double_t binLimInvMass[nBinsInvMass+1];
\r
678 for(Int_t i=0; i<=nBinsInvMass; i++) binLimInvMass[i]=(Double_t)minInvMass + (maxInvMass-minInvMass)/nBinsInvMass*(Double_t)i ;
\r
680 //******************
\r
682 //******************
\r
684 fListHist = new TList();
\r
685 fListHist->SetOwner();
\r
688 fHistEV = new TH2D("fHistEV", "events", 3, 0, 3, 3, 0,3);
\r
690 // Event plane as function of phiep, centrality
\r
691 const Int_t nDima=5;
\r
692 Int_t nBina[nDima] = {nBinsPhi,nBinsPhi,nBinsPhi,nBinsPhi,nBinsC};
\r
693 fEventPlane = new THnSparseF("EventPlane","EventPlane",nDima,nBina);
\r
694 fEventPlane->SetBinEdges(0,binLimPhi);
\r
695 fEventPlane->SetBinEdges(1,binLimPhi);
\r
696 fEventPlane->SetBinEdges(2,binLimPhi);
\r
697 fEventPlane->SetBinEdges(3,binLimPhi);
\r
698 fEventPlane->SetBinEdges(4,binLimC);
\r
699 fEventPlane->Sumw2();
\r
701 // Event Plane after subtraction as function of phiep, centrality, pt, eta
\r
702 const Int_t nDimb=2;
\r
703 Int_t nBinb[nDimb] = {nBinsPhi, nBinsC};
\r
704 fEventPlaneaftersubtraction = new THnSparseF("EventPlane_aftersubtraction","EventPlane_aftersubtraction",nDimb,nBinb);
\r
705 fEventPlaneaftersubtraction->SetBinEdges(0,binLimPhi);
\r
706 fEventPlaneaftersubtraction->SetBinEdges(1,binLimC);
\r
707 fEventPlaneaftersubtraction->Sumw2();
\r
709 // Monitoring of the event Plane cos(2phi) sin(2phi) centrality
\r
710 const Int_t nDimi=3;
\r
711 Int_t nBini[nDimi] = {nBinsCos, nBinsCos, nBinsCMore};
\r
712 fCosSin2phiep = new THnSparseF("CosSin2phiep","CosSin2phiep",nDimi,nBini);
\r
713 fCosSin2phiep->SetBinEdges(0,binLimCos);
\r
714 fCosSin2phiep->SetBinEdges(1,binLimCos);
\r
715 fCosSin2phiep->SetBinEdges(2,binLimCMore);
\r
716 fCosSin2phiep->Sumw2();
\r
718 // Monitoring Event plane after subtraction of the track
\r
719 const Int_t nDime=4;
\r
720 Int_t nBine[nDime] = {nBinsCos, nBinsC, nBinsPt, nBinsEta};
\r
721 fCos2phie = new THnSparseF("cos2phie","cos2phie",nDime,nBine);
\r
722 fCos2phie->SetBinEdges(2,binLimPt);
\r
723 fCos2phie->SetBinEdges(3,binLimEta);
\r
724 fCos2phie->SetBinEdges(0,binLimCos);
\r
725 fCos2phie->SetBinEdges(1,binLimC);
\r
726 fCos2phie->Sumw2();
\r
727 fSin2phie = new THnSparseF("sin2phie","sin2phie",nDime,nBine);
\r
728 fSin2phie->SetBinEdges(2,binLimPt);
\r
729 fSin2phie->SetBinEdges(3,binLimEta);
\r
730 fSin2phie->SetBinEdges(0,binLimCos);
\r
731 fSin2phie->SetBinEdges(1,binLimC);
\r
732 fSin2phie->Sumw2();
\r
733 fCos2phiep = new THnSparseF("cos2phiep","cos2phiep",nDime,nBine);
\r
734 fCos2phiep->SetBinEdges(2,binLimPt);
\r
735 fCos2phiep->SetBinEdges(3,binLimEta);
\r
736 fCos2phiep->SetBinEdges(0,binLimCos);
\r
737 fCos2phiep->SetBinEdges(1,binLimC);
\r
738 fCos2phiep->Sumw2();
\r
739 fSin2phiep = new THnSparseF("sin2phiep","sin2phiep",nDime,nBine);
\r
740 fSin2phiep->SetBinEdges(2,binLimPt);
\r
741 fSin2phiep->SetBinEdges(3,binLimEta);
\r
742 fSin2phiep->SetBinEdges(0,binLimCos);
\r
743 fSin2phiep->SetBinEdges(1,binLimC);
\r
744 fSin2phiep->Sumw2();
\r
745 fSin2phiephiep = new THnSparseF("sin2phie_phiep","sin2phie_phiep",nDime,nBine);
\r
746 fSin2phiephiep->SetBinEdges(2,binLimPt);
\r
747 fSin2phiephiep->SetBinEdges(3,binLimEta);
\r
748 fSin2phiephiep->SetBinEdges(0,binLimCos);
\r
749 fSin2phiephiep->SetBinEdges(1,binLimC);
\r
750 fSin2phiephiep->Sumw2();
\r
752 // Resolution cosres_abc centrality
\r
753 const Int_t nDimfbis=4;
\r
754 Int_t nBinfbis[nDimfbis] = {nBinsCos,nBinsCos,nBinsCos,nBinsCMore};
\r
755 fCosResabc = new THnSparseF("CosRes_abc","CosRes_abc",nDimfbis,nBinfbis);
\r
756 fCosResabc->SetBinEdges(0,binLimCos);
\r
757 fCosResabc->SetBinEdges(1,binLimCos);
\r
758 fCosResabc->SetBinEdges(2,binLimCos);
\r
759 fCosResabc->SetBinEdges(3,binLimCMore);
\r
760 fCosResabc->Sumw2();
\r
762 const Int_t nDimfbiss=4;
\r
763 Int_t nBinfbiss[nDimfbiss] = {nBinsCos,nBinsCos,nBinsCos,nBinsC};
\r
764 fSinResabc = new THnSparseF("SinRes_abc","SinRes_abc",nDimfbiss,nBinfbiss);
\r
765 fSinResabc->SetBinEdges(0,binLimCos);
\r
766 fSinResabc->SetBinEdges(1,binLimCos);
\r
767 fSinResabc->SetBinEdges(2,binLimCos);
\r
768 fSinResabc->SetBinEdges(3,binLimC);
\r
769 fSinResabc->Sumw2();
\r
771 // Profile cosres centrality with 3 subevents
\r
772 fProfileCosResab = new TProfile("ProfileCosRes_a_b","ProfileCosRes_a_b",nBinsCMore,binLimCMore);
\r
773 fProfileCosResab->Sumw2();
\r
774 fProfileCosResac = new TProfile("ProfileCosRes_a_c","ProfileCosRes_a_c",nBinsCMore,binLimCMore);
\r
775 fProfileCosResac->Sumw2();
\r
776 fProfileCosResbc = new TProfile("ProfileCosRes_b_c","ProfileCosRes_b_c",nBinsCMore,binLimCMore);
\r
777 fProfileCosResbc->Sumw2();
\r
779 // Resolution cosres centrality
\r
780 const Int_t nDimf=2;
\r
781 Int_t nBinf[nDimf] = {nBinsCos, nBinsCMore};
\r
782 fCosRes = new THnSparseF("CosRes","CosRes",nDimf,nBinf);
\r
783 fCosRes->SetBinEdges(0,binLimCos);
\r
784 fCosRes->SetBinEdges(1,binLimCMore);
\r
787 const Int_t nDimff=2;
\r
788 Int_t nBinff[nDimff] = {nBinsCos, nBinsC};
\r
789 fSinRes = new THnSparseF("SinRes","SinRes",nDimff,nBinff);
\r
790 fSinRes->SetBinEdges(0,binLimCos);
\r
791 fSinRes->SetBinEdges(1,binLimC);
\r
794 // Profile cosres centrality
\r
795 fProfileCosRes = new TProfile("ProfileCosRes","ProfileCosRes",nBinsCMore,binLimCMore);
\r
796 fProfileCosRes->Sumw2();
\r
798 // Debugging tracking steps
\r
799 const Int_t nDimTrStep=2;
\r
800 Int_t nBinTrStep[nDimTrStep] = {nBinsPt,nBinsStep};
\r
801 fTrackingCuts = new THnSparseF("TrackingCuts","TrackingCuts",nDimTrStep,nBinTrStep);
\r
802 fTrackingCuts->SetBinEdges(0,binLimPt);
\r
803 fTrackingCuts->SetBinEdges(1,binLimStep);
\r
804 fTrackingCuts->Sumw2();
\r
807 const Int_t nDimg=5;
\r
808 Int_t nBing[nDimg] = {nBinsPhi,nBinsC,nBinsPt, nBinsCharge,nBinsEtaLess};
\r
809 fDeltaPhiMaps = new THnSparseF("DeltaPhiMaps","DeltaPhiMaps",nDimg,nBing);
\r
810 fDeltaPhiMaps->SetBinEdges(0,binLimPhi);
\r
811 fDeltaPhiMaps->SetBinEdges(1,binLimC);
\r
812 fDeltaPhiMaps->SetBinEdges(2,binLimPt);
\r
813 fDeltaPhiMaps->SetBinEdges(3,binLimCharge);
\r
814 fDeltaPhiMaps->SetBinEdges(4,binLimEtaLess);
\r
815 fDeltaPhiMaps->Sumw2();
\r
817 // Maps delta phi contamination
\r
818 const Int_t nDimgcont=4;
\r
819 Int_t nBingcont[nDimgcont] = {nBinsPhiLess,nBinsC,nBinsPt, nBinsTPCdEdx};
\r
820 fDeltaPhiMapsContamination = new THnSparseF("DeltaPhiMapsContamination","DeltaPhiMapsContamination",nDimgcont,nBingcont);
\r
821 fDeltaPhiMapsContamination->SetBinEdges(0,binLimPhiLess);
\r
822 fDeltaPhiMapsContamination->SetBinEdges(1,binLimC);
\r
823 fDeltaPhiMapsContamination->SetBinEdges(2,binLimPt);
\r
824 fDeltaPhiMapsContamination->SetBinEdges(3,binLimTPCdEdx);
\r
825 fDeltaPhiMapsContamination->Sumw2();
\r
829 const Int_t nDimgb=3;
\r
830 Int_t nBingb[nDimgb] = {nBinsPhi,nBinsC,nBinsPt};
\r
832 fDeltaPhiMapsBeforePID = new THnSparseF("DeltaPhiMapsBeforePID","DeltaPhiMapsBeforePID",nDimgb,nBingb);
\r
833 fDeltaPhiMapsBeforePID->SetBinEdges(0,binLimPhi);
\r
834 fDeltaPhiMapsBeforePID->SetBinEdges(1,binLimC);
\r
835 fDeltaPhiMapsBeforePID->SetBinEdges(2,binLimPt);
\r
836 fDeltaPhiMapsBeforePID->Sumw2();
\r
839 fDeltaPhiMapsTaggedPhotonic = new THnSparseF("DeltaPhiMapsTaggedPhotonic","DeltaPhiMapsTaggedPhotonic",nDimgb,nBingb);
\r
840 fDeltaPhiMapsTaggedPhotonic->SetBinEdges(0,binLimPhi);
\r
841 fDeltaPhiMapsTaggedPhotonic->SetBinEdges(1,binLimC);
\r
842 fDeltaPhiMapsTaggedPhotonic->SetBinEdges(2,binLimPt);
\r
843 fDeltaPhiMapsTaggedPhotonic->Sumw2();
\r
845 fDeltaPhiMapsTaggedNonPhotonic = new THnSparseF("DeltaPhiMapsTaggedNonPhotonic","DeltaPhiMapsTaggedNonPhotonic",nDimgb,nBingb);
\r
846 fDeltaPhiMapsTaggedNonPhotonic->SetBinEdges(0,binLimPhi);
\r
847 fDeltaPhiMapsTaggedNonPhotonic->SetBinEdges(1,binLimC);
\r
848 fDeltaPhiMapsTaggedNonPhotonic->SetBinEdges(2,binLimPt);
\r
849 fDeltaPhiMapsTaggedNonPhotonic->Sumw2();
\r
851 fDeltaPhiMapsTaggedPhotonicLS = new THnSparseF("DeltaPhiMapsTaggedPhotonicLS","DeltaPhiMapsTaggedPhotonicLS",nDimgb,nBingb);
\r
852 fDeltaPhiMapsTaggedPhotonicLS->SetBinEdges(0,binLimPhi);
\r
853 fDeltaPhiMapsTaggedPhotonicLS->SetBinEdges(1,binLimC);
\r
854 fDeltaPhiMapsTaggedPhotonicLS->SetBinEdges(2,binLimPt);
\r
855 fDeltaPhiMapsTaggedPhotonicLS->Sumw2();
\r
858 const Int_t nDimh=5;
\r
859 Int_t nBinh[nDimh] = {nBinsCos,nBinsC,nBinsPt,nBinsCharge,nBinsEtaLess};
\r
860 fCosPhiMaps = new THnSparseF("CosPhiMaps","CosPhiMaps",nDimh,nBinh);
\r
861 fCosPhiMaps->SetBinEdges(0,binLimCos);
\r
862 fCosPhiMaps->SetBinEdges(1,binLimC);
\r
863 fCosPhiMaps->SetBinEdges(2,binLimPt);
\r
864 fCosPhiMaps->SetBinEdges(3,binLimCharge);
\r
865 fCosPhiMaps->SetBinEdges(4,binLimEtaLess);
\r
866 fCosPhiMaps->Sumw2();
\r
868 const Int_t nDimhb=3;
\r
869 Int_t nBinhb[nDimhb] = {nBinsCos,nBinsC,nBinsPt};
\r
871 fCosPhiMapsBeforePID = new THnSparseF("CosPhiMapsBeforePID","CosPhiMapsBeforePID",nDimhb,nBinhb);
\r
872 fCosPhiMapsBeforePID->SetBinEdges(0,binLimCos);
\r
873 fCosPhiMapsBeforePID->SetBinEdges(1,binLimC);
\r
874 fCosPhiMapsBeforePID->SetBinEdges(2,binLimPt);
\r
875 fCosPhiMapsBeforePID->Sumw2();
\r
877 fCosPhiMapsTaggedPhotonic = new THnSparseF("CosPhiMapsTaggedPhotonic","CosPhiMapsTaggedPhotonic",nDimhb,nBinhb);
\r
878 fCosPhiMapsTaggedPhotonic->SetBinEdges(0,binLimCos);
\r
879 fCosPhiMapsTaggedPhotonic->SetBinEdges(1,binLimC);
\r
880 fCosPhiMapsTaggedPhotonic->SetBinEdges(2,binLimPt);
\r
881 fCosPhiMapsTaggedPhotonic->Sumw2();
\r
883 fCosPhiMapsTaggedNonPhotonic = new THnSparseF("CosPhiMapsTaggedNonPhotonic","CosPhiMapsTaggedNonPhotonic",nDimhb,nBinhb);
\r
884 fCosPhiMapsTaggedNonPhotonic->SetBinEdges(0,binLimCos);
\r
885 fCosPhiMapsTaggedNonPhotonic->SetBinEdges(1,binLimC);
\r
886 fCosPhiMapsTaggedNonPhotonic->SetBinEdges(2,binLimPt);
\r
887 fCosPhiMapsTaggedNonPhotonic->Sumw2();
\r
889 fCosPhiMapsTaggedPhotonicLS = new THnSparseF("CosPhiMapsTaggedPhotonicLS","CosPhiMapsTaggedPhotonicLS",nDimhb,nBinhb);
\r
890 fCosPhiMapsTaggedPhotonicLS->SetBinEdges(0,binLimCos);
\r
891 fCosPhiMapsTaggedPhotonicLS->SetBinEdges(1,binLimC);
\r
892 fCosPhiMapsTaggedPhotonicLS->SetBinEdges(2,binLimPt);
\r
893 fCosPhiMapsTaggedPhotonicLS->Sumw2();
\r
895 // Profile Maps cos phi
\r
896 fProfileCosPhiMaps = new TProfile2D("ProfileCosPhiMaps","ProfileCosPhiMaps",nBinsC,binLimC,nBinsPt,binLimPt);
\r
897 fProfileCosPhiMaps->Sumw2();
\r
899 // Background study
\r
900 const Int_t nDimMCSource=3;
\r
901 Int_t nBinMCSource[nDimMCSource] = {nBinsC,nBinsPt,nBinsSource};
\r
902 fMCSourceDeltaPhiMaps = new THnSparseF("MCSourceDeltaPhiMaps","MCSourceDeltaPhiMaps",nDimMCSource,nBinMCSource);
\r
903 fMCSourceDeltaPhiMaps->SetBinEdges(0,binLimC);
\r
904 fMCSourceDeltaPhiMaps->SetBinEdges(1,binLimPt);
\r
905 fMCSourceDeltaPhiMaps->SetBinEdges(2,binLimSource);
\r
906 fMCSourceDeltaPhiMaps->Sumw2();
\r
908 // Maps invmass opposite
\r
909 const Int_t nDimOppSign=5;
\r
910 Int_t nBinOppSign[nDimOppSign] = {nBinsPhi,nBinsC,nBinsPt,nBinsInvMass,nBinsSource};
\r
911 fOppSignDeltaPhiMaps = new THnSparseF("OppSignDeltaPhiMaps","OppSignDeltaPhiMaps",nDimOppSign,nBinOppSign);
\r
912 fOppSignDeltaPhiMaps->SetBinEdges(0,binLimPhi);
\r
913 fOppSignDeltaPhiMaps->SetBinEdges(1,binLimC);
\r
914 fOppSignDeltaPhiMaps->SetBinEdges(2,binLimPt);
\r
915 fOppSignDeltaPhiMaps->SetBinEdges(3,binLimInvMass);
\r
916 fOppSignDeltaPhiMaps->SetBinEdges(4,binLimSource);
\r
917 fOppSignDeltaPhiMaps->Sumw2();
\r
919 // Maps invmass same sign
\r
920 const Int_t nDimSameSign=5;
\r
921 Int_t nBinSameSign[nDimSameSign] = {nBinsPhi,nBinsC,nBinsPt,nBinsInvMass,nBinsSource};
\r
922 fSameSignDeltaPhiMaps = new THnSparseF("SameSignDeltaPhiMaps","SameSignDeltaPhiMaps",nDimSameSign,nBinSameSign);
\r
923 fSameSignDeltaPhiMaps->SetBinEdges(0,binLimPhi);
\r
924 fSameSignDeltaPhiMaps->SetBinEdges(1,binLimC);
\r
925 fSameSignDeltaPhiMaps->SetBinEdges(2,binLimPt);
\r
926 fSameSignDeltaPhiMaps->SetBinEdges(3,binLimInvMass);
\r
927 fSameSignDeltaPhiMaps->SetBinEdges(4,binLimSource);
\r
928 fSameSignDeltaPhiMaps->Sumw2();
\r
930 // Maps angle same sign
\r
931 const Int_t nDimAngleSameSign=3;
\r
932 Int_t nBinAngleSameSign[nDimAngleSameSign] = {nBinsAngle,nBinsC,nBinsSource};
\r
933 fSameSignAngle = new THnSparseF("SameSignAngleMaps","SameSignAngleMaps",nDimAngleSameSign,nBinAngleSameSign);
\r
934 fSameSignAngle->SetBinEdges(0,binLimAngle);
\r
935 fSameSignAngle->SetBinEdges(1,binLimC);
\r
936 fSameSignAngle->SetBinEdges(2,binLimSource);
\r
937 fSameSignAngle->Sumw2();
\r
939 // Maps angle opp sign
\r
940 const Int_t nDimAngleOppSign=3;
\r
941 Int_t nBinAngleOppSign[nDimAngleOppSign] = {nBinsAngle,nBinsC,nBinsSource};
\r
942 fOppSignAngle = new THnSparseF("OppSignAngleMaps","OppSignAngleMaps",nDimAngleOppSign,nBinAngleOppSign);
\r
943 fOppSignAngle->SetBinEdges(0,binLimAngle);
\r
944 fOppSignAngle->SetBinEdges(1,binLimC);
\r
945 fOppSignAngle->SetBinEdges(2,binLimSource);
\r
946 fOppSignAngle->Sumw2();
\r
949 //**************************
\r
951 //******************************
\r
953 //fListHist->Add(qaCutsRP);
\r
954 fListHist->Add(fPIDqa->MakeList("HFEpidQA"));
\r
955 fListHist->Add(fPIDBackgroundqa->MakeList("HFEpidBackgroundQA"));
\r
956 fListHist->Add(fHistEV);
\r
957 fListHist->Add(fProfileCosRes);
\r
958 fListHist->Add(fProfileCosResab);
\r
959 fListHist->Add(fProfileCosResac);
\r
960 fListHist->Add(fProfileCosResbc);
\r
961 fListHist->Add(fCosSin2phiep);
\r
962 fListHist->Add(fEventPlane);
\r
963 fListHist->Add(fEventPlaneaftersubtraction);
\r
964 fListHist->Add(fCos2phie);
\r
965 fListHist->Add(fSin2phie);
\r
966 fListHist->Add(fCos2phiep);
\r
967 fListHist->Add(fSin2phiep);
\r
968 fListHist->Add(fSin2phiephiep);
\r
969 fListHist->Add(fCosRes);
\r
970 fListHist->Add(fCosResabc);
\r
971 fListHist->Add(fSinRes);
\r
972 fListHist->Add(fSinResabc);
\r
973 fListHist->Add(fTrackingCuts);
\r
974 fListHist->Add(fDeltaPhiMapsBeforePID);
\r
975 fListHist->Add(fCosPhiMapsBeforePID);
\r
976 fListHist->Add(fDeltaPhiMaps);
\r
977 fListHist->Add(fDeltaPhiMapsContamination);
\r
978 fListHist->Add(fCosPhiMaps);
\r
979 fListHist->Add(fProfileCosPhiMaps);
\r
980 fListHist->Add(fDeltaPhiMapsTaggedPhotonic);
\r
981 fListHist->Add(fCosPhiMapsTaggedPhotonic);
\r
982 fListHist->Add(fDeltaPhiMapsTaggedNonPhotonic);
\r
983 fListHist->Add(fCosPhiMapsTaggedNonPhotonic);
\r
984 fListHist->Add(fDeltaPhiMapsTaggedPhotonicLS);
\r
985 fListHist->Add(fCosPhiMapsTaggedPhotonicLS);
\r
986 fListHist->Add(fMCSourceDeltaPhiMaps);
\r
987 fListHist->Add(fOppSignDeltaPhiMaps);
\r
988 fListHist->Add(fSameSignDeltaPhiMaps);
\r
989 fListHist->Add(fSameSignAngle);
\r
990 fListHist->Add(fOppSignAngle);
\r
993 if(fHFEVZEROEventPlane && (fDebugLevel > 2)) fListHist->Add(fHFEVZEROEventPlane->GetOutputList());
\r
996 PostData(1, fListHist);
\r
1001 //________________________________________________________________________
\r
1002 void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/)
\r
1005 // Loop over event
\r
1008 Double_t massElectron = 0.000511;
\r
1009 Double_t mcReactionPlane = 0.0;
\r
1010 Bool_t eventplanedefined = kTRUE;
\r
1012 Float_t cntr = 0.0;
\r
1013 Double_t binct = 11.5;
\r
1014 Double_t binctMore = 20.5;
\r
1015 Double_t binctLess = -0.5;
\r
1016 Float_t binctt = -1.0;
\r
1018 Double_t valuecossinephiep[3];
\r
1019 Double_t valuensparsea[5];
\r
1020 Double_t valuensparseabis[5];
\r
1021 Double_t valuensparsee[4];
\r
1022 Double_t valuensparsef[2];
\r
1023 Double_t valuensparsefsin[2];
\r
1024 Double_t valuensparsefbis[4];
\r
1025 Double_t valuensparsefbissin[4];
\r
1026 Double_t valuensparseg[5];
\r
1027 Double_t valuensparseh[5];
\r
1028 Double_t valuensparsehprofile[3];
\r
1029 Double_t valuensparseMCSourceDeltaPhiMaps[3];
\r
1030 Double_t valuetrackingcuts[2];
\r
1031 Double_t valuedeltaphicontamination[4];
\r
1033 AliMCEvent *mcEvent = MCEvent();
\r
1034 AliMCParticle *mctrack = NULL;
\r
1040 //AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
\r
1041 //if(!esd) return;
\r
1042 AliCentrality *centrality = fInputEvent->GetCentrality();
\r
1043 //printf("Got the centrality\n");
\r
1044 if(!centrality) return;
\r
1045 cntr = centrality->GetCentralityPercentile("V0M");
\r
1046 if((0.0< cntr) && (cntr<5.0)) binct = 0.5;
\r
1047 if((5.0< cntr) && (cntr<10.0)) binct = 1.5;
\r
1048 if((10.0< cntr) && (cntr<20.0)) binct = 2.5;
\r
1049 if((20.0< cntr) && (cntr<30.0)) binct = 3.5;
\r
1050 if((30.0< cntr) && (cntr<40.0)) binct = 4.5;
\r
1051 if((40.0< cntr) && (cntr<50.0)) binct = 5.5;
\r
1052 if((50.0< cntr) && (cntr<60.0)) binct = 6.5;
\r
1053 if((60.0< cntr) && (cntr<70.0)) binct = 7.5;
\r
1054 if((70.0< cntr) && (cntr<80.0)) binct = 8.5;
\r
1055 if((80.0< cntr) && (cntr<90.0)) binct = 9.5;
\r
1056 if((90.0< cntr) && (cntr<100.0)) binct = 10.5;
\r
1058 if((0.< cntr) && (cntr < 20.)) binctt = 0.5;
\r
1059 if((20.< cntr) && (cntr < 40.)) binctt = 1.5;
\r
1060 if((40.< cntr) && (cntr < 80.)) binctt = 2.5;
\r
1062 if((0.0< cntr) && (cntr<5.0)) binctMore = 0.5;
\r
1063 if((5.0< cntr) && (cntr<10.0)) binctMore = 1.5;
\r
1064 if((10.0< cntr) && (cntr<15.0)) binctMore = 2.5;
\r
1065 if((15.0< cntr) && (cntr<20.0)) binctMore = 3.5;
\r
1066 if((20.0< cntr) && (cntr<25.0)) binctMore = 4.5;
\r
1067 if((25.0< cntr) && (cntr<30.0)) binctMore = 5.5;
\r
1068 if((30.0< cntr) && (cntr<35.0)) binctMore = 6.5;
\r
1069 if((35.0< cntr) && (cntr<40.0)) binctMore = 7.5;
\r
1070 if((40.0< cntr) && (cntr<45.0)) binctMore = 8.5;
\r
1071 if((45.0< cntr) && (cntr<50.0)) binctMore = 9.5;
\r
1072 if((50.0< cntr) && (cntr<55.0)) binctMore = 10.5;
\r
1073 if((55.0< cntr) && (cntr<60.0)) binctMore = 11.5;
\r
1074 if((60.0< cntr) && (cntr<65.0)) binctMore = 12.5;
\r
1075 if((65.0< cntr) && (cntr<70.0)) binctMore = 13.5;
\r
1076 if((70.0< cntr) && (cntr<75.0)) binctMore = 14.5;
\r
1077 if((75.0< cntr) && (cntr<80.0)) binctMore = 15.5;
\r
1078 if((80.0< cntr) && (cntr<85.0)) binctMore = 16.5;
\r
1079 if((85.0< cntr) && (cntr<90.0)) binctMore = 17.5;
\r
1080 if((90.0< cntr) && (cntr<95.0)) binctMore = 18.5;
\r
1081 if((95.0< cntr) && (cntr<100.0)) binctMore = 19.5;
\r
1086 if(binct > 11.0) return;
\r
1089 valuensparsea[4] = binct;
\r
1090 valuensparseabis[1] = binct;
\r
1091 valuensparsee[1] = binct;
\r
1092 valuensparsef[1] = binctMore;
\r
1093 valuensparsefsin[1] = binct;
\r
1094 valuensparsefbis[3] = binctMore;
\r
1095 valuensparsefbissin[3] = binct;
\r
1096 valuensparseg[1] = binct;
\r
1097 valuensparseh[1] = binct;
\r
1098 valuensparsehprofile[1] = binct;
\r
1099 valuecossinephiep[2] = binctMore;
\r
1100 valuensparseMCSourceDeltaPhiMaps[0] = binct;
\r
1101 valuedeltaphicontamination[1] = binct;
\r
1103 //////////////////////
\r
1105 //////////////////////
\r
1107 Int_t runnumber = fInputEvent->GetRunNumber();
\r
1108 //printf("Run number %d\n",runnumber);
\r
1110 if(!fPID->IsInitialized()){
\r
1111 // Initialize PID with the given run number
\r
1112 fPID->InitializePID(runnumber);
\r
1114 if(!fPIDTOFOnly->IsInitialized()){
\r
1115 // Initialize PID with the given run number
\r
1116 fPIDTOFOnly->InitializePID(runnumber);
\r
1120 if(!fPIDBackground->IsInitialized()){
\r
1121 // Initialize PID with the given run number
\r
1122 fPIDBackground->InitializePID(runnumber);
\r
1125 fHFECuts->SetRecEvent(fInputEvent);
\r
1126 if(mcEvent) fHFECuts->SetMCEvent(mcEvent);
\r
1133 AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
\r
1135 //printf("No PID response set\n");
\r
1138 fPID->SetPIDResponse(pidResponse);
\r
1139 fPIDTOFOnly->SetPIDResponse(pidResponse);
\r
1140 fPIDBackground->SetPIDResponse(pidResponse);
\r
1142 fHistEV->Fill(binctt,0.0);
\r
1145 //////////////////
\r
1147 //////////////////
\r
1148 if(!fHFECuts->CheckEventCuts("fEvRecCuts", fInputEvent)) {
\r
1149 //printf("Do not pass the event cut\n");
\r
1150 PostData(1, fListHist);
\r
1154 fHistEV->Fill(binctt,1.0);
\r
1156 ////////////////////////////////////
\r
1157 // First method event plane
\r
1158 ////////////////////////////////////
\r
1160 AliEventplane* vEPa = fInputEvent->GetEventplane();
\r
1161 Float_t eventPlanea = 0.0;
\r
1162 Float_t eventPlaneTPC = 0.0;
\r
1163 Float_t eventPlaneV0A = 0.0;
\r
1164 Float_t eventPlaneV0C = 0.0;
\r
1165 Float_t eventPlaneV0 = 0.0;
\r
1166 TVector2 *standardQ = 0x0;
\r
1167 TVector2 *qsub1a = 0x0;
\r
1168 TVector2 *qsub2a = 0x0;
\r
1172 if(fHFEVZEROEventPlane && (!fAODAnalysis)){
\r
1174 //AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
\r
1175 //if(!esd) return;
\r
1177 fHFEVZEROEventPlane->ProcessEvent(fInputEvent);
\r
1179 if(TMath::Abs(fHFEVZEROEventPlane->GetEventPlaneV0A()+100) < 0.0000001) eventPlaneV0A = -100.0;
\r
1181 eventPlaneV0A = TVector2::Phi_0_2pi(fHFEVZEROEventPlane->GetEventPlaneV0A());
\r
1182 if(eventPlaneV0A > TMath::Pi()) eventPlaneV0A = eventPlaneV0A - TMath::Pi();
\r
1185 if(TMath::Abs(fHFEVZEROEventPlane->GetEventPlaneV0C()+100) < 0.0000001) eventPlaneV0C = -100.0;
\r
1187 eventPlaneV0C = TVector2::Phi_0_2pi(fHFEVZEROEventPlane->GetEventPlaneV0C());
\r
1188 if(eventPlaneV0C > TMath::Pi()) eventPlaneV0C = eventPlaneV0C - TMath::Pi();
\r
1191 if(TMath::Abs(fHFEVZEROEventPlane->GetEventPlaneV0()+100) < 0.0000001) eventPlaneV0 = -100.0;
\r
1193 eventPlaneV0 = TVector2::Phi_0_2pi(fHFEVZEROEventPlane->GetEventPlaneV0());
\r
1194 if(eventPlaneV0 > TMath::Pi()) eventPlaneV0 = eventPlaneV0 - TMath::Pi();
\r
1200 eventPlaneV0 = TVector2::Phi_0_2pi(vEPa->GetEventplane("V0", fInputEvent,2));
\r
1201 //printf("eventPlaneV0 %f\n",eventPlaneV0);
\r
1202 if(eventPlaneV0 > TMath::Pi()) eventPlaneV0 = eventPlaneV0 - TMath::Pi();
\r
1203 //printf("eventPlaneV0 %f\n",eventPlaneV0);
\r
1204 eventPlaneV0A = TVector2::Phi_0_2pi(vEPa->GetEventplane("V0A", fInputEvent,2));
\r
1205 if(eventPlaneV0A > TMath::Pi()) eventPlaneV0A = eventPlaneV0A - TMath::Pi();
\r
1206 eventPlaneV0C = TVector2::Phi_0_2pi(vEPa->GetEventplane("V0C", fInputEvent,2));
\r
1207 if(eventPlaneV0C > TMath::Pi()) eventPlaneV0C = eventPlaneV0C - TMath::Pi();
\r
1213 standardQ = vEPa->GetQVector();
\r
1214 Double_t qx = -1.0;
\r
1215 Double_t qy = -1.0;
\r
1217 qx = standardQ->X();
\r
1218 qy = standardQ->Y();
\r
1220 TVector2 qVectorfortrack;
\r
1221 qVectorfortrack.Set(qx,qy);
\r
1222 eventPlaneTPC = TVector2::Phi_0_2pi(qVectorfortrack.Phi())/2.;
\r
1224 // Choose the one used for v2
\r
1226 if(fVZEROEventPlane) eventPlanea = eventPlaneV0;
\r
1227 if(fVZEROEventPlaneA) eventPlanea = eventPlaneV0A;
\r
1228 if(fVZEROEventPlaneC) eventPlanea = eventPlaneV0C;
\r
1229 if(!fVZEROEventPlane) eventPlanea = eventPlaneTPC;
\r
1231 valuecossinephiep[0] = TMath::Cos(2*eventPlanea);
\r
1232 valuecossinephiep[1] = TMath::Sin(2*eventPlanea);
\r
1234 Float_t eventPlanesub1a = -100.0;
\r
1235 Float_t eventPlanesub2a = -100.0;
\r
1236 Double_t diffsub1sub2a = -100.0;
\r
1237 Double_t diffsub1sub2asin = -100.0;
\r
1238 Double_t diffsubasubb = -100.0;
\r
1239 Double_t diffsubasubc = -100.0;
\r
1240 Double_t diffsubbsubc = -100.0;
\r
1241 Double_t diffsubasubbsin = -100.0;
\r
1242 Double_t diffsubasubcsin = -100.0;
\r
1243 Double_t diffsubbsubcsin = -100.0;
\r
1245 //if(fVZEROEventPlane) {
\r
1246 diffsubasubb = TMath::Cos(2.*(eventPlaneV0A - eventPlaneV0C));
\r
1247 diffsubasubc = TMath::Cos(2.*(eventPlaneV0A - eventPlaneTPC));
\r
1248 diffsubbsubc = TMath::Cos(2.*(eventPlaneV0C - eventPlaneTPC));
\r
1250 diffsubasubbsin = TMath::Sin(2.*(eventPlaneV0A - eventPlaneV0C));
\r
1251 diffsubasubcsin = TMath::Sin(2.*(eventPlaneV0A - eventPlaneTPC));
\r
1252 diffsubbsubcsin = TMath::Sin(2.*(eventPlaneV0C - eventPlaneTPC));
\r
1255 qsub1a = vEPa->GetQsub1();
\r
1256 qsub2a = vEPa->GetQsub2();
\r
1257 if(qsub1a) eventPlanesub1a = TVector2::Phi_0_2pi(qsub1a->Phi())/2.;
\r
1258 if(qsub2a) eventPlanesub2a = TVector2::Phi_0_2pi(qsub2a->Phi())/2.;
\r
1259 if(qsub1a && qsub2a) {
\r
1260 diffsub1sub2a = TMath::Cos(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.));
\r
1261 diffsub1sub2asin = TMath::Sin(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.));
\r
1266 /////////////////////////////////////////////////////////
\r
1267 // Cut for event with event plane reconstructed by all
\r
1268 ////////////////////////////////////////////////////////
\r
1270 //if(!fVZEROEventPlane) {
\r
1271 // if(!standardQ) {
\r
1272 // eventplanedefined = kFALSE;
\r
1273 //PostData(1, fListHist);
\r
1278 if((!standardQ) || (!qsub1a) || (!qsub2a)) {
\r
1279 //printf("No event plane\n");
\r
1283 ///////////////////////
\r
1285 //////////////////////
\r
1287 Int_t nbtracks = fInputEvent->GetNumberOfTracks();
\r
1288 //printf("Number of tracks %d\n",nbtracks);
\r
1290 fcutsRP->SetEvent( InputEvent(), MCEvent());
\r
1291 fcutsPOI->SetEvent( InputEvent(), MCEvent());
\r
1292 if( fflowEvent ){
\r
1293 fflowEvent->~AliFlowEvent();
\r
1294 new(fflowEvent) AliFlowEvent(fcutsRP,fcutsPOI);
\r
1295 }else fflowEvent = new AliFlowEvent(fcutsRP,fcutsPOI);
\r
1296 if(mcEvent && mcEvent->GenEventHeader()) {
\r
1297 fflowEvent->SetMCReactionPlaneAngle(mcEvent);
\r
1298 //if reaction plane not set from elsewhere randomize it before adding flow
\r
1299 //if (!fflowEvent->IsSetMCReactionPlaneAngle()) fflowEvent->SetMCReactionPlaneAngle(gRandom->Uniform(0.0,TMath::TwoPi()));
\r
1300 mcReactionPlane = TVector2::Phi_0_2pi(fflowEvent->GetMCReactionPlaneAngle());
\r
1301 if(mcReactionPlane > TMath::Pi()) mcReactionPlane = mcReactionPlane - TMath::Pi();
\r
1302 //printf("MC reaction plane %f\n",mcReactionPlane);
\r
1304 fflowEvent->SetReferenceMultiplicity( nbtracks );
\r
1305 fflowEvent->DefineDeadZone(0,0,0,0);
\r
1306 //fflowEvent.TagSubeventsInEta(-0.8,-0.1,0.1,0.8);
\r
1311 if(fUseMCReactionPlane) {
\r
1312 eventPlanea = mcReactionPlane;
\r
1313 diffsub1sub2a = 0.0;
\r
1317 //////////////////////
\r
1319 //////////////////////
\r
1321 if(eventplanedefined) {
\r
1323 fHistEV->Fill(binctt,2.0);
\r
1326 valuensparsea[0] = eventPlaneV0A;
\r
1327 valuensparsea[1] = eventPlaneV0C;
\r
1328 valuensparsea[2] = eventPlaneTPC;
\r
1329 valuensparsea[3] = eventPlaneV0;
\r
1330 fEventPlane->Fill(&valuensparsea[0]);
\r
1333 if(fDebugLevel > 5) fCosSin2phiep->Fill(&valuecossinephiep[0]);
\r
1335 if(!fVZEROEventPlane) {
\r
1336 valuensparsef[0] = diffsub1sub2a;
\r
1337 fCosRes->Fill(&valuensparsef[0]);
\r
1338 valuensparsefsin[0] = diffsub1sub2asin;
\r
1339 if(fDebugLevel > 5) fSinRes->Fill(&valuensparsefsin[0]);
\r
1340 if(fDebugLevel > 5) {
\r
1341 fProfileCosRes->Fill(valuensparsef[1],valuensparsef[0]);
\r
1345 valuensparsefbis[0] = diffsubasubb;
\r
1346 valuensparsefbis[1] = diffsubbsubc;
\r
1347 valuensparsefbis[2] = diffsubasubc;
\r
1348 fCosResabc->Fill(&valuensparsefbis[0]);
\r
1349 valuensparsefbissin[0] = diffsubasubbsin;
\r
1350 valuensparsefbissin[1] = diffsubbsubcsin;
\r
1351 valuensparsefbissin[2] = diffsubasubcsin;
\r
1352 if(fDebugLevel > 5) fSinResabc->Fill(&valuensparsefbissin[0]);
\r
1353 if(fDebugLevel > 5) {
\r
1354 fProfileCosResab->Fill(valuensparsefbis[3],valuensparsefbis[0]);
\r
1355 fProfileCosResac->Fill(valuensparsefbis[3],valuensparsefbis[1]);
\r
1356 fProfileCosResbc->Fill(valuensparsefbis[3],valuensparsefbis[2]);
\r
1362 ////////////////////////////////////////
\r
1363 // Loop to determine pool background
\r
1364 /////////////////////////////////////////
\r
1365 if( fArraytrack ){
\r
1366 fArraytrack->~TArrayI();
\r
1367 new(fArraytrack) TArrayI(nbtracks);
\r
1370 fArraytrack = new TArrayI(nbtracks);
\r
1372 fCounterPoolBackground = 0;
\r
1373 for(Int_t k = 0; k < nbtracks; k++){
\r
1375 AliVTrack *track = (AliVTrack *) fInputEvent->GetTrack(k);
\r
1376 if(!track) continue;
\r
1379 Bool_t survivedbackground = kTRUE;
\r
1380 if(fAODAnalysis) {
\r
1381 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
\r
1383 AliESDtrack esdTrack(aodtrack);
\r
1384 // set the TPC cluster info
\r
1385 esdTrack.SetTPCClusterMap(aodtrack->GetTPCClusterMap());
\r
1386 esdTrack.SetTPCSharedMap(aodtrack->GetTPCSharedMap());
\r
1387 esdTrack.SetTPCPointsF(aodtrack->GetTPCNclsF());
\r
1388 // needed to calculate the impact parameters
\r
1389 AliAODEvent *aodeventu = dynamic_cast<AliAODEvent *>(fInputEvent);
\r
1391 AliAODVertex *vAOD = aodeventu->GetPrimaryVertex();
\r
1392 Double_t bfield = aodeventu->GetMagneticField();
\r
1393 Double_t pos[3],cov[6];
\r
1394 vAOD->GetXYZ(pos);
\r
1395 vAOD->GetCovarianceMatrix(cov);
\r
1396 const AliESDVertex vESD(pos,cov,100.,100);
\r
1397 esdTrack.RelateToVertex(&vESD,bfield,3.);
\r
1399 if(!fHFEBackgroundCuts->IsSelected(&esdTrack)) {
\r
1400 survivedbackground = kFALSE;
\r
1405 AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
\r
1407 if(!fHFEBackgroundCuts->IsSelected(esdtrack)) survivedbackground = kFALSE;
\r
1411 if(survivedbackground) {
\r
1413 AliHFEpidObject hfetrack2;
\r
1414 if(!fAODAnalysis) hfetrack2.SetAnalysisType(AliHFEpidObject::kESDanalysis);
\r
1415 else hfetrack2.SetAnalysisType(AliHFEpidObject::kAODanalysis);
\r
1416 hfetrack2.SetRecTrack(track);
\r
1417 hfetrack2.SetCentrality((Int_t)binct);
\r
1418 //printf("centrality %f and %d\n",binct,hfetrack.GetCentrality());
\r
1419 hfetrack2.SetPbPb();
\r
1420 if(fPIDBackground->IsSelected(&hfetrack2,0x0,"recTrackCont",fPIDBackgroundqa)) {
\r
1421 fArraytrack->AddAt(k,fCounterPoolBackground);
\r
1422 fCounterPoolBackground++;
\r
1423 //printf("fCounterPoolBackground %d, track %d\n",fCounterPoolBackground,k);
\r
1428 // Look at kink mother in case of AOD
\r
1429 Int_t numberofvertices = 1;
\r
1430 AliAODEvent *aodevent = NULL;
\r
1431 Int_t numberofmotherkink = 0;
\r
1432 if(fAODAnalysis) {
\r
1433 aodevent = dynamic_cast<AliAODEvent *>(fInputEvent);
\r
1435 numberofvertices = aodevent->GetNumberOfVertices();
\r
1438 Double_t listofmotherkink[numberofvertices];
\r
1439 if(fAODAnalysis && aodevent) {
\r
1440 for(Int_t ivertex=0; ivertex < numberofvertices; ivertex++) {
\r
1441 AliAODVertex *aodvertex = aodevent->GetVertex(ivertex);
\r
1442 if(!aodvertex) continue;
\r
1443 if(aodvertex->GetType()==AliAODVertex::kKink) {
\r
1444 AliAODTrack *mother = (AliAODTrack *) aodvertex->GetParent();
\r
1445 if(!mother) continue;
\r
1446 Int_t idmother = mother->GetID();
\r
1447 listofmotherkink[numberofmotherkink] = idmother;
\r
1448 //printf("ID %d\n",idmother);
\r
1449 numberofmotherkink++;
\r
1454 //////////////////////////
\r
1455 // Loop over track
\r
1456 //////////////////////////
\r
1458 for(Int_t k = 0; k < nbtracks; k++){
\r
1460 AliVTrack *track = (AliVTrack *) fInputEvent->GetTrack(k);
\r
1461 if(!track) continue;
\r
1463 if(fAODAnalysis) {
\r
1464 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
\r
1466 AliError("AOD track is not there");
\r
1469 //printf("Find AOD track on\n");
\r
1471 if(aodtrack->GetFlags() != fFlags) continue; // Only process AOD tracks where the HFE is set
\r
1477 valuetrackingcuts[0] = track->Pt();
\r
1478 valuetrackingcuts[1] = 0;
\r
1480 // RecKine: ITSTPC cuts
\r
1481 if(!fHFECuts->CheckParticleCuts(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack, (TObject *)track)) continue;
\r
1482 if(fDebugLevel > 3) fTrackingCuts->Fill(&valuetrackingcuts[0]);
\r
1484 // Reject kink mother
\r
1485 if(fAODAnalysis) {
\r
1486 Bool_t kinkmotherpass = kTRUE;
\r
1487 for(Int_t kinkmother = 0; kinkmother < numberofmotherkink; kinkmother++) {
\r
1488 if(track->GetID() == listofmotherkink[kinkmother]) {
\r
1489 kinkmotherpass = kFALSE;
\r
1493 if(!kinkmotherpass) continue;
\r
1496 AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
\r
1498 if(esdtrack->GetKinkIndex(0) != 0) continue;
\r
1499 } // Quick and dirty fix to reject both kink mothers and daughters
\r
1502 valuetrackingcuts[1] = 1;
\r
1503 if(fDebugLevel > 3) fTrackingCuts->Fill(&valuetrackingcuts[0]);
\r
1505 if(!fHFECuts->CheckParticleCuts(AliHFEcuts::kStepRecPrim + AliHFEcuts::kNcutStepsMCTrack, (TObject *)track)) continue;
\r
1506 valuetrackingcuts[1] = 2;
\r
1507 if(fDebugLevel > 3) fTrackingCuts->Fill(&valuetrackingcuts[0]);
\r
1509 // HFEcuts: ITS layers cuts
\r
1510 if(!fHFECuts->CheckParticleCuts(AliHFEcuts::kStepHFEcutsITS + AliHFEcuts::kNcutStepsMCTrack, (TObject *)track)) continue;
\r
1511 valuetrackingcuts[1] = 3;
\r
1512 if(fDebugLevel > 3) fTrackingCuts->Fill(&valuetrackingcuts[0]);
\r
1514 // HFE cuts: TOF PID and mismatch flag
\r
1515 if(!fHFECuts->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTOF + AliHFEcuts::kNcutStepsMCTrack, (TObject *)track)) continue;
\r
1516 valuetrackingcuts[1] = 4;
\r
1517 if(fDebugLevel > 3) fTrackingCuts->Fill(&valuetrackingcuts[0]);
\r
1519 // HFE cuts: TPC PID cleanup
\r
1520 if(!fHFECuts->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTPC + AliHFEcuts::kNcutStepsMCTrack, (TObject *)track)) continue;
\r
1521 valuetrackingcuts[1] = 5;
\r
1522 if(fDebugLevel > 3) fTrackingCuts->Fill(&valuetrackingcuts[0]);
\r
1524 // HFEcuts: Nb of tracklets TRD0
\r
1525 if(!fHFECuts->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack, (TObject *)track)) continue;
\r
1526 valuetrackingcuts[1] = 6;
\r
1527 if(fDebugLevel > 3) fTrackingCuts->Fill(&valuetrackingcuts[0]);
\r
1531 //printf("Survived\n");
\r
1533 /////////////////////////////////////////////////////////
\r
1534 // Subtract candidate from TPC event plane
\r
1535 ////////////////////////////////////////////////////////
\r
1536 Float_t eventplanesubtracted = 0.0;
\r
1538 //if(eventplanedefined && (!fVZEROEventPlane)) {
\r
1539 if(!fVZEROEventPlane) {
\r
1540 // Subtract the tracks from the event plane
\r
1541 Double_t qX = standardQ->X() - vEPa->GetQContributionX(track); //Modify the components: subtract the track you want to look at with your analysis
\r
1542 Double_t qY = standardQ->Y() - vEPa->GetQContributionY(track); //Modify the components: subtract the track you want to look at with your analysis
\r
1543 TVector2 newQVectorfortrack;
\r
1544 newQVectorfortrack.Set(qX,qY);
\r
1545 eventplanesubtracted = TVector2::Phi_0_2pi(newQVectorfortrack.Phi())/2;
\r
1547 else eventplanesubtracted = eventPlanea;
\r
1549 ///////////////////////////////////////////
\r
1551 //////////////////////////////////////////
\r
1552 Bool_t fillEventPlane = kTRUE;
\r
1553 if(!fVZEROEventPlane){
\r
1554 //if((!qsub1a) || (!qsub2a) || (!eventplanedefined)) fillEventPlane = kFALSE;
\r
1555 if((!qsub1a) || (!qsub2a)) fillEventPlane = kFALSE;
\r
1556 if(fSubEtaGapTPC) {
\r
1557 if(track->Eta() < (- fEtaGap/2.)) eventplanesubtracted = eventPlanesub1a;
\r
1558 else if(track->Eta() > (fEtaGap/2.)) eventplanesubtracted = eventPlanesub2a;
\r
1559 else fillEventPlane = kFALSE;
\r
1566 if(fUseMCReactionPlane) {
\r
1567 eventplanesubtracted = mcReactionPlane;
\r
1568 fillEventPlane = kTRUE;
\r
1571 //////////////////////////////////////////////////////////////////////////////
\r
1572 ///////////////////////////AFTERBURNER
\r
1573 Double_t phitrack = track->Phi();
\r
1574 if (fAfterBurnerOn)
\r
1576 phitrack = GetPhiAfterAddV2(track->Phi(),mcReactionPlane);
\r
1578 //////////////////////////////////////////////////////////////////////////////
\r
1581 ///////////////////////
\r
1582 // Calculate deltaphi
\r
1583 ///////////////////////
\r
1585 // Suppose phi track is between 0.0 and phi
\r
1586 Double_t deltaphi = TVector2::Phi_0_2pi(phitrack - eventplanesubtracted);
\r
1587 if(deltaphi > TMath::Pi()) deltaphi = deltaphi - TMath::Pi();
\r
1589 ////////////////////////////////
\r
1590 // Determine the deltaphi bin
\r
1591 ///////////////////////////////
\r
1594 if(((deltaphi<(TMath::Pi()/4.)) && (deltaphi>0.0)) || ((deltaphi>(3*TMath::Pi()/4.)) && (deltaphi<TMath::Pi()))) valuedeltaphicontamination[0] = 0.5;
\r
1596 if((deltaphi>(TMath::Pi()/4.)) && (deltaphi<(3*TMath::Pi()/4.))) valuedeltaphicontamination[0] = 1.5;
\r
1598 ////////////////////////////////////////
\r
1599 // Define variables
\r
1600 ///////////////////////////////////////
\r
1603 valuedeltaphicontamination[2] = track->Pt();
\r
1604 valuensparsee[2] = track->Pt();
\r
1605 valuensparsee[3] = track->Eta();
\r
1606 valuensparseg[2] = track->Pt();
\r
1607 valuensparseh[2] = track->Pt();
\r
1608 valuensparsehprofile[2] = track->Pt();
\r
1609 valuensparseMCSourceDeltaPhiMaps[1] = track->Pt();
\r
1610 if(track->Charge() > 0.0) {
\r
1611 valuensparseg[3] = 0.2;
\r
1612 valuensparseh[3] = 0.2;
\r
1615 valuensparseg[3] = -0.2;
\r
1616 valuensparseh[3] = -0.2;
\r
1618 valuensparseh[4] = track->Eta();
\r
1619 valuensparseg[4] = track->Eta();
\r
1621 //printf("charge %d\n",(Int_t)track->Charge());
\r
1623 ////////////////////////
\r
1624 // Fill before PID
\r
1625 ///////////////////////
\r
1627 if(fDebugLevel > 5) {
\r
1629 valuensparseg[0] = deltaphi;
\r
1630 if(fillEventPlane) fDeltaPhiMapsBeforePID->Fill(&valuensparseg[0]);
\r
1633 valuensparseh[0] = TMath::Cos(2*TVector2::Phi_mpi_pi(phitrack-eventplanesubtracted));
\r
1634 if(fillEventPlane) {
\r
1635 fCosPhiMapsBeforePID->Fill(&valuensparseh[0]);
\r
1639 ////////////////////////
\r
1641 ////////////////////////
\r
1643 // Apply PID for Data
\r
1646 AliHFEpidObject hfetrack;
\r
1647 if(!fAODAnalysis) hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
\r
1648 else hfetrack.SetAnalysisType(AliHFEpidObject::kAODanalysis);
\r
1649 hfetrack.SetRecTrack(track);
\r
1650 hfetrack.SetCentrality((Int_t)binct);
\r
1651 //printf("centrality %f and %d\n",binct,hfetrack.GetCentrality());
\r
1652 hfetrack.SetPbPb();
\r
1655 if(fPIDTOFOnly->IsSelected(&hfetrack,0x0,"recTrackCont",0x0)) {
\r
1656 Float_t nsigma = pidResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
\r
1657 valuedeltaphicontamination[3] = nsigma;
\r
1658 fDeltaPhiMapsContamination->Fill(&valuedeltaphicontamination[0]);
\r
1661 // Complete PID TOF+TPC
\r
1662 if(!fPID->IsSelected(&hfetrack,0x0,"recTrackCont",fPIDqa)) {
\r
1667 if(!mcEvent) continue;
\r
1668 if(!(mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel()))))) continue;
\r
1669 //printf("PdgCode %d\n",TMath::Abs(mctrack->Particle()->GetPdgCode()));
\r
1670 if(TMath::Abs(mctrack->Particle()->GetPdgCode())!=11) continue;
\r
1675 /////////////////////////////////////////////////////////////////////////////
\r
1676 // Add candidate to AliFlowEvent for POI and subtract from RP if needed
\r
1677 ////////////////////////////////////////////////////////////////////////////
\r
1678 Int_t idtrack = static_cast<AliVTrack*>(track)->GetID();
\r
1679 Bool_t found = kFALSE;
\r
1680 Int_t numberoffound = 0;
\r
1681 //printf("A: Number of tracks %d\n",fflowEvent->NumberOfTracks());
\r
1682 for(Int_t iRPs=0; iRPs< fflowEvent->NumberOfTracks(); iRPs++) {
\r
1683 AliFlowTrack *iRP = (AliFlowTrack*) (fflowEvent->GetTrack(iRPs));
\r
1684 //if(!iRP->InRPSelection()) continue;
\r
1685 if( TMath::Abs(idtrack) == TMath::Abs(iRP->GetID()) ) {
\r
1686 iRP->SetForPOISelection(kTRUE);
\r
1691 //printf("Found %d mal\n",numberoffound);
\r
1693 AliFlowCandidateTrack *sTrack = (AliFlowCandidateTrack*) MakeTrack(massElectron,track->Pt(),track->Phi(), track->Eta());
\r
1694 sTrack->SetID(idtrack);
\r
1695 fflowEvent->AddTrack(sTrack);
\r
1696 //printf("Add the track\n");
\r
1698 //printf("B: Number of tracks %d\n",fflowEvent->NumberOfTracks());
\r
1702 /////////////////////
\r
1703 // Fill THnSparseF
\r
1704 /////////////////////
\r
1707 valuensparseabis[0] = eventplanesubtracted;
\r
1708 if((fillEventPlane) && (fDebugLevel > 5)) fEventPlaneaftersubtraction->Fill(&valuensparseabis[0]);
\r
1711 if(fDebugLevel > 5)
\r
1714 valuensparsee[0] = TMath::Cos(2*phitrack);
\r
1715 fCos2phie->Fill(&valuensparsee[0]);
\r
1716 valuensparsee[0] = TMath::Sin(2*phitrack);
\r
1717 fSin2phie->Fill(&valuensparsee[0]);
\r
1719 valuensparsee[0] = TMath::Cos(2*eventplanesubtracted);
\r
1720 if(fillEventPlane) fCos2phiep->Fill(&valuensparsee[0]);
\r
1721 valuensparsee[0] = TMath::Sin(2*eventplanesubtracted);
\r
1722 if(fillEventPlane) fSin2phiep->Fill(&valuensparsee[0]);
\r
1723 valuensparsee[0] = TMath::Sin(2*TVector2::Phi_mpi_pi(phitrack-eventplanesubtracted));
\r
1724 if(fillEventPlane) fSin2phiephiep->Fill(&valuensparsee[0]);
\r
1729 valuensparseg[0] = deltaphi;
\r
1730 if(fillEventPlane) fDeltaPhiMaps->Fill(&valuensparseg[0]);
\r
1733 valuensparseh[0] = TMath::Cos(2*TVector2::Phi_mpi_pi(phitrack-eventplanesubtracted));
\r
1734 if(fillEventPlane) {
\r
1735 fCosPhiMaps->Fill(&valuensparseh[0]);
\r
1736 if(fDebugLevel > 0) {
\r
1737 fProfileCosPhiMaps->Fill(valuensparsehprofile[1],valuensparsehprofile[2],valuensparseh[0]);
\r
1741 if(fDebugLevel > 1) {
\r
1744 Int_t indexmother = -1;
\r
1745 source = FindMother(TMath::Abs(track->GetLabel()),mcEvent, indexmother);
\r
1746 valuensparseMCSourceDeltaPhiMaps[2] = source;
\r
1747 if(mcEvent) fMCSourceDeltaPhiMaps->Fill(&valuensparseMCSourceDeltaPhiMaps[0]);
\r
1748 Int_t taggedvalue = LookAtNonHFE(k,track,fInputEvent,mcEvent,binct,deltaphi,source,indexmother);
\r
1749 if(fillEventPlane) {
\r
1750 // No opposite charge partner found in the invariant mass choosen
\r
1751 if((taggedvalue!=2) && (taggedvalue!=6)) {
\r
1752 fDeltaPhiMapsTaggedNonPhotonic->Fill(&valuensparseg[0]);
\r
1753 fCosPhiMapsTaggedNonPhotonic->Fill(&valuensparseh[0]);
\r
1755 // One opposite charge partner found in the invariant mass choosen
\r
1756 if((taggedvalue==2) || (taggedvalue==6)) {
\r
1757 fDeltaPhiMapsTaggedPhotonic->Fill(&valuensparseg[0]);
\r
1758 fCosPhiMapsTaggedPhotonic->Fill(&valuensparseh[0]);
\r
1760 // One same charge partner found in the invariant mass choosen
\r
1761 if((taggedvalue==4) || (taggedvalue==6)) {
\r
1762 fDeltaPhiMapsTaggedPhotonicLS->Fill(&valuensparseg[0]);
\r
1763 fCosPhiMapsTaggedPhotonicLS->Fill(&valuensparseh[0]);
\r
1770 //////////////////////////////////////////////////////////////////////////////
\r
1771 ///////////////////////////AFTERBURNER
\r
1772 if (fAfterBurnerOn)
\r
1774 fflowEvent->AddFlow(fV1,fV2,fV3,fV4,fV5); //add flow
\r
1775 fflowEvent->CloneTracks(fNonFlowNumberOfTrackClones); //add nonflow by cloning tracks
\r
1777 //////////////////////////////////////////////////////////////////////////////
\r
1781 for(Int_t bincless = 0; bincless < fNbBinsCentralityQCumulant; bincless++) {
\r
1782 if((fBinCentralityLess[bincless]< cntr) && (cntr < fBinCentralityLess[bincless+1])) PostData(bincless+2,fflowEvent);
\r
1786 delete fArraytrack;
\r
1787 fArraytrack = NULL;
\r
1790 PostData(1, fListHist);
\r
1794 //______________________________________________________________________________
\r
1795 AliFlowCandidateTrack *AliAnalysisTaskHFEFlow::MakeTrack( Double_t mass,
\r
1796 Double_t pt, Double_t phi, Double_t eta) {
\r
1798 // Make Track (Not needed actually)
\r
1801 AliFlowCandidateTrack *sTrack = new AliFlowCandidateTrack();
\r
1802 sTrack->SetMass(mass);
\r
1803 sTrack->SetPt(pt);
\r
1804 sTrack->SetPhi(phi);
\r
1805 sTrack->SetEta(eta);
\r
1806 sTrack->SetForPOISelection(kTRUE);
\r
1807 sTrack->SetForRPSelection(kFALSE);
\r
1810 //_________________________________________________________________________________
\r
1811 Double_t AliAnalysisTaskHFEFlow::GetPhiAfterAddV2(Double_t phi,Double_t reactionPlaneAngle) const
\r
1814 // Adds v2, uses Newton-Raphson iteration
\r
1816 Double_t phiend=phi;
\r
1817 Double_t phi0=phi;
\r
1820 Double_t phiprev=0.;
\r
1822 for (Int_t i=0; i<fMaxNumberOfIterations; i++)
\r
1824 phiprev=phiend; //store last value for comparison
\r
1825 f = phiend-phi0+fV2*TMath::Sin(2.*(phiend-reactionPlaneAngle));
\r
1826 fp = 1.0+2.0*fV2*TMath::Cos(2.*(phiend-reactionPlaneAngle)); //first derivative
\r
1828 if (TMath::AreEqualAbs(phiprev,phiend,fPrecisionPhi)) break;
\r
1832 //_____________________________________________________________________________________________
\r
1833 Int_t AliAnalysisTaskHFEFlow::LookAtNonHFE(Int_t iTrack1, AliVTrack *track1, AliVEvent *vEvent, AliMCEvent *mcEvent,Int_t binct,Double_t deltaphi,Int_t source,Int_t indexmother)
\r
1836 // Look At Non HFE
\r
1839 // return -1 if nothing
\r
1840 // return 2 if opposite charge within the mass range found
\r
1841 // return 4 if like charge within the mass range found
\r
1842 // return 6 if opposite charge and like charge within the mass range found
\r
1845 Int_t taggedphotonic = -1;
\r
1847 Bool_t oppositetaggedphotonic = kFALSE;
\r
1848 Bool_t sametaggedphotonic = kFALSE;
\r
1850 //printf("fCounterPoolBackground %d in LookAtNonHFE!!!\n",fCounterPoolBackground);
\r
1851 if(!fArraytrack) return taggedphotonic;
\r
1852 //printf("process track %d\n",iTrack1);
\r
1857 Double_t valuensparseDeltaPhiMaps[5];
\r
1858 Double_t valueangle[3];
\r
1860 valuensparseDeltaPhiMaps[1] = binct;
\r
1861 valuensparseDeltaPhiMaps[2] = track1->Pt();
\r
1862 valuensparseDeltaPhiMaps[0] = deltaphi;
\r
1863 valuensparseDeltaPhiMaps[4] = source;
\r
1865 valueangle[2] = source;
\r
1866 valueangle[1] = binct;
\r
1869 Int_t pdg1 = CheckPdg(TMath::Abs(track1->GetLabel()),mcEvent);
\r
1870 Int_t numberfound = 0;
\r
1873 Double_t bfield = vEvent->GetMagneticField();
\r
1875 // Get Primary vertex
\r
1876 const AliVVertex *pVtx = vEvent->GetPrimaryVertex();
\r
1878 for(Int_t idex = 0; idex < fCounterPoolBackground; idex++)
\r
1881 Int_t iTrack2 = fArraytrack->At(idex);
\r
1882 //printf("track %d\n",iTrack2);
\r
1883 AliVTrack* track2 = (AliVTrack *) vEvent->GetTrack(iTrack2);
\r
1886 printf("ERROR: Could not receive track %d\n", iTrack2);
\r
1889 if(iTrack2==iTrack1) continue;
\r
1890 //printf("Different\n");
\r
1892 // Reset the MC info
\r
1893 valueangle[2] = source;
\r
1894 valuensparseDeltaPhiMaps[4] = source;
\r
1896 // track cuts and PID already done
\r
1899 Int_t pdg2 = -100;
\r
1901 Int_t source2 = 0;
\r
1902 Int_t indexmother2 = -1;
\r
1903 source2 = FindMother(TMath::Abs(track2->GetLabel()),mcEvent, indexmother2);
\r
1904 pdg2 = CheckPdg(TMath::Abs(track2->GetLabel()),mcEvent);
\r
1905 if(source2 >=0 ) {
\r
1906 if((indexmother2 == indexmother) && (source == source2) && ((pdg1*pdg2)<0.0)) {
\r
1907 if(source == kElectronfromconversion) {
\r
1908 valueangle[2] = kElectronfromconversionboth;
\r
1909 valuensparseDeltaPhiMaps[4] = kElectronfromconversionboth;
\r
1912 if(source == kElectronfrompi0) {
\r
1913 valueangle[2] = kElectronfrompi0both;
\r
1914 valuensparseDeltaPhiMaps[4] = kElectronfrompi0both;
\r
1916 if(source == kElectronfrometa) {
\r
1917 valueangle[2] = kElectronfrometaboth;
\r
1918 valuensparseDeltaPhiMaps[4] = kElectronfrometaboth;
\r
1924 if(fAlgorithmMA && (!fAODAnalysis))
\r
1927 AliESDtrack *esdtrack2 = dynamic_cast<AliESDtrack *>(track2);
\r
1928 AliESDtrack *esdtrack1 = dynamic_cast<AliESDtrack *>(track1);
\r
1929 if((!esdtrack2) || (!esdtrack1)) continue;
\r
1934 Double_t xt1; //radial position track 1 at the DCA point
\r
1935 Double_t xt2; //radial position track 2 at the DCA point
\r
1936 //DCA track1-track2
\r
1937 Double_t dca12 = esdtrack2->GetDCA(esdtrack1,bfield,xt2,xt1);
\r
1940 if(dca12 > fMaxdca) continue;
\r
1942 //Momento of the track extrapolated to DCA track-track
\r
1944 Bool_t hasdcaT1 = esdtrack1->GetPxPyPzAt(xt1,bfield,p1);
\r
1946 Bool_t hasdcaT2 = esdtrack2->GetPxPyPzAt(xt2,bfield,p2);
\r
1948 if(!hasdcaT1 || !hasdcaT2) AliWarning("It could be a problem in the extrapolation");
\r
1950 //track1-track2 Invariant Mass
\r
1951 Double_t eMass = 0.000510998910; //Electron mass in GeV
\r
1952 Double_t pP1 = sqrt(p1[0]*p1[0]+p1[1]*p1[1]+p1[2]*p1[2]); //Track 1 momentum
\r
1953 Double_t pP2 = sqrt(p2[0]*p2[0]+p2[1]*p2[1]+p2[2]*p2[2]); //Track 2 momentum
\r
1954 Double_t eE1 = TMath::Sqrt(pP1*pP1+eMass*eMass);
\r
1955 Double_t eE2 = TMath::Sqrt(pP2*pP2+eMass*eMass);
\r
1957 //TLorentzVector v1(p1[0],p1[1],p1[2],sqrt(eMass*eMass+pP1*pP1));
\r
1958 //TLorentzVector v2(p2[0],p2[1],p2[2],sqrt(eMass*eMass+pP2*pP2));
\r
1959 //Double_t imass = (v1+v2).M(); //Invariant Mass
\r
1960 //Double_t angle3D = v1.Angle(v2.Vect()); //Opening Angle (Total Angle)
\r
1963 v3D1.SetXYZ(p1[0],p1[1],p1[2]);
\r
1964 v3D2.SetXYZ(p2[0],p2[1],p2[2]);
\r
1965 Double_t openingangle = TVector2::Phi_0_2pi(v3D2.Angle(v3D1));
\r
1968 TVector3 motherrec = v3D1 + v3D2;
\r
1969 Double_t invmass = TMath::Sqrt((eE1+eE2)*(eE1+eE2)-(motherrec.Px()*motherrec.Px()+motherrec.Py()*motherrec.Py()+motherrec.Pz()*motherrec.Pz()));
\r
1972 //TVector3 vectordiff = v3D1 - v3D2;
\r
1973 //Double_t diffphi = TVector2::Phi_0_2pi(vectordiff.Phi());
\r
1974 //Double_t massxy = TMath::Sqrt((eE1+eE2)*(eE1+eE2)-(pP1*pP1+pP2*pP2+2*pP1*pP2*TMath::Cos(diffphi)));
\r
1977 //Double_t difftheta = TVector2::Phi_0_2pi(vectordiff.Eta());
\r
1978 //Double_t massrz = TMath::Sqrt((eE1+eE2)*(eE1+eE2)-(pP1*pP1+pP2*pP2+2*pP1*pP2*TMath::Cos(difftheta)));
\r
1981 Float_t fCharge1 = track1->Charge();
\r
1982 Float_t fCharge2 = track2->Charge();
\r
1985 //valueangle[0] = diffphi;
\r
1986 //valueangle[1] = difftheta;
\r
1987 valueangle[0] = openingangle;
\r
1988 if((fCharge1*fCharge2)>0.0) fSameSignAngle->Fill(&valueangle[0]);
\r
1989 else fOppSignAngle->Fill(&valueangle[0]);
\r
1992 if(openingangle > fMaxopening3D) continue;
\r
1993 //if(difftheta > fMaxopeningtheta) continue;
\r
1994 //if(diffphi > fMaxopeningphi) continue;
\r
1997 valuensparseDeltaPhiMaps[3] = invmass;
\r
1998 if((fCharge1*fCharge2)>0.0) fSameSignDeltaPhiMaps->Fill(&valuensparseDeltaPhiMaps[0]);
\r
1999 else fOppSignDeltaPhiMaps->Fill(&valuensparseDeltaPhiMaps[0]);
\r
2002 if(invmass < fMaxInvmass) {
\r
2003 if((fCharge1*fCharge2)<0.0) oppositetaggedphotonic=kTRUE;
\r
2004 if((fCharge1*fCharge2)>0.0) sametaggedphotonic=kTRUE;
\r
2011 Int_t fPDGtrack1 = 11;
\r
2012 Int_t fPDGtrack2 = 11;
\r
2014 Float_t fCharge1 = track1->Charge();
\r
2015 Float_t fCharge2 = track2->Charge();
\r
2017 if(fCharge1>0) fPDGtrack1 = -11;
\r
2018 if(fCharge2>0) fPDGtrack2 = -11;
\r
2020 AliKFParticle ktrack1(*track1, fPDGtrack1);
\r
2021 AliKFParticle ktrack2(*track2, fPDGtrack2);
\r
2022 AliKFParticle recoGamma(ktrack1, ktrack2);
\r
2024 //Reconstruction Cuts
\r
2025 if(recoGamma.GetNDF()<1) continue;
\r
2026 Double_t chi2OverNDF = recoGamma.GetChi2()/recoGamma.GetNDF();
\r
2027 if(TMath::Sqrt(TMath::Abs(chi2OverNDF))>fChi2OverNDFCut) continue;
\r
2030 //Double_t dca12 = ktrack1.GetDistanceFromParticle(ktrack2);
\r
2031 //if(dca12 > fMaxdca) continue;
\r
2033 // if set mass constraint
\r
2034 if(fSetMassConstraint && pVtx) {
\r
2035 AliKFVertex primV(*pVtx);
\r
2036 primV += recoGamma;
\r
2039 recoGamma.SetProductionVertex(primV);
\r
2040 recoGamma.SetMassConstraint(0,0.0001);
\r
2046 recoGamma.GetMass(imass,width);
\r
2048 //Opening Angle (Total Angle)
\r
2049 Double_t angle = ktrack1.GetAngle(ktrack2);
\r
2050 valueangle[0] = angle;
\r
2051 if((fCharge1*fCharge2)>0.0) fSameSignAngle->Fill(&valueangle[0]);
\r
2052 else fOppSignAngle->Fill(&valueangle[0]);
\r
2055 if(angle > fMaxopening3D) continue;
\r
2058 valuensparseDeltaPhiMaps[3] = imass;
\r
2059 if((fCharge1*fCharge2)>0.0) fSameSignDeltaPhiMaps->Fill(&valuensparseDeltaPhiMaps[0]);
\r
2061 fOppSignDeltaPhiMaps->Fill(&valuensparseDeltaPhiMaps[0]);
\r
2063 if(valueangle[2] == kElectronfromconversionboth) {
\r
2064 printf("Reconstructed charge1 %f, charge 2 %f and invmass %f\n",fCharge1,fCharge2,imass);
\r
2065 printf("MC charge1 %d, charge 2 %d\n",pdg1,pdg2);
\r
2066 printf("DCA %f\n",dca12);
\r
2067 printf("Number of found %d\n",numberfound);
\r
2073 if(imass < fMaxInvmass) {
\r
2074 if((fCharge1*fCharge2)<0.0) oppositetaggedphotonic=kTRUE;
\r
2075 if((fCharge1*fCharge2)>0.0) sametaggedphotonic=kTRUE;
\r
2080 if(oppositetaggedphotonic && sametaggedphotonic){
\r
2081 taggedphotonic = 6;
\r
2084 if(!oppositetaggedphotonic && sametaggedphotonic){
\r
2085 taggedphotonic = 4;
\r
2088 if(oppositetaggedphotonic && !sametaggedphotonic){
\r
2089 taggedphotonic = 2;
\r
2093 return taggedphotonic;
\r
2095 //_________________________________________________________________________
\r
2096 Int_t AliAnalysisTaskHFEFlow::FindMother(Int_t tr, AliMCEvent *mcEvent, Int_t &indexmother){
\r
2098 // Find the mother if MC
\r
2101 if(!mcEvent) return 0;
\r
2103 Int_t pdg = CheckPdg(tr,mcEvent);
\r
2104 if(TMath::Abs(pdg)!= 11) {
\r
2106 return kNoElectron;
\r
2109 indexmother = IsMotherGamma(tr,mcEvent);
\r
2110 if(indexmother > 0) return kElectronfromconversion;
\r
2111 indexmother = IsMotherPi0(tr,mcEvent);
\r
2112 if(indexmother > 0) return kElectronfrompi0;
\r
2113 indexmother = IsMotherC(tr,mcEvent);
\r
2114 if(indexmother > 0) return kElectronfromC;
\r
2115 indexmother = IsMotherB(tr,mcEvent);
\r
2116 if(indexmother > 0) return kElectronfromB;
\r
2117 indexmother = IsMotherEta(tr,mcEvent);
\r
2118 if(indexmother > 0) return kElectronfrometa;
\r
2120 return kElectronfromother;
\r
2124 //____________________________________________________________________________________________________________
\r
2125 Int_t AliAnalysisTaskHFEFlow::CheckPdg(Int_t tr, AliMCEvent* mcEvent) {
\r
2128 // Return the pdg of the particle
\r
2132 Int_t pdgcode = -1;
\r
2133 if(tr < 0) return pdgcode;
\r
2135 if(!mcEvent) return pdgcode;
\r
2137 AliVParticle *mctrack = mcEvent->GetTrack(tr);
\r
2140 if(mctrack->IsA() == AliMCParticle::Class()) {
\r
2141 AliMCParticle *mctrackesd = NULL;
\r
2142 if(!(mctrackesd = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return pdgcode;
\r
2143 pdgcode = mctrackesd->PdgCode();
\r
2146 if(mctrack->IsA() == AliAODMCParticle::Class()) {
\r
2147 AliAODMCParticle *mctrackaod = NULL;
\r
2148 if(!(mctrackaod = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return pdgcode;
\r
2149 pdgcode = mctrackaod->GetPdgCode();
\r
2156 //____________________________________________________________________________________________________________
\r
2157 Int_t AliAnalysisTaskHFEFlow::IsMotherGamma(Int_t tr, AliMCEvent* mcEvent) {
\r
2160 // Return the lab of gamma mother or -1 if not gamma
\r
2163 if(tr < 0) return -1;
\r
2164 AliVParticle *mctrack = mcEvent->GetTrack(tr);
\r
2166 if(mctrack->IsA() == AliMCParticle::Class()) {
\r
2167 AliMCParticle *mctrackesd = NULL;
\r
2168 if(!(mctrackesd = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
\r
2169 TParticle *particle = 0x0;
\r
2170 particle = mctrackesd->Particle();
\r
2172 if(!particle) return -1;
\r
2173 Int_t imother = particle->GetFirstMother();
\r
2174 if(imother < 0) return -1;
\r
2175 AliMCParticle *mothertrack = NULL;
\r
2176 if(!(mothertrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
\r
2177 TParticle * mother = mothertrack->Particle();
\r
2178 if(!mother) return -1;
\r
2180 Int_t pdg = mother->GetPdgCode();
\r
2181 if(TMath::Abs(pdg) == 22) return imother;
\r
2182 if(TMath::Abs(pdg) == 11) {
\r
2183 return IsMotherGamma(imother,mcEvent);
\r
2188 if(mctrack->IsA() == AliAODMCParticle::Class()) {
\r
2189 AliAODMCParticle *mctrackaod = NULL;
\r
2190 if(!(mctrackaod = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
\r
2192 Int_t imother = mctrackaod->GetMother();
\r
2193 if(imother < 0) return -1;
\r
2194 AliAODMCParticle *mothertrack = NULL;
\r
2195 if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
\r
2197 Int_t pdg = mothertrack->GetPdgCode();
\r
2198 if(TMath::Abs(pdg) == 22) return imother;
\r
2199 if(TMath::Abs(pdg) == 11) {
\r
2200 return IsMotherGamma(imother,mcEvent);
\r
2211 //____________________________________________________________________________________________________________
\r
2212 Int_t AliAnalysisTaskHFEFlow::IsMotherPi0(Int_t tr, AliMCEvent* mcEvent) {
\r
2215 // Return the lab of pi0 mother or -1 if not pi0
\r
2218 if(tr < 0) return -1;
\r
2219 AliVParticle *mctrack = mcEvent->GetTrack(tr);
\r
2221 if(mctrack->IsA() == AliMCParticle::Class()) {
\r
2222 AliMCParticle *mctrackesd = NULL;
\r
2223 if(!(mctrackesd = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
\r
2224 TParticle *particle = 0x0;
\r
2225 particle = mctrackesd->Particle();
\r
2227 if(!particle) return -1;
\r
2228 Int_t imother = particle->GetFirstMother();
\r
2229 if(imother < 0) return -1;
\r
2230 AliMCParticle *mothertrack = NULL;
\r
2231 if(!(mothertrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
\r
2232 TParticle * mother = mothertrack->Particle();
\r
2233 if(!mother) return -1;
\r
2235 Int_t pdg = mother->GetPdgCode();
\r
2236 if(TMath::Abs(pdg) == 111) return imother;
\r
2237 if(TMath::Abs(pdg) == 11) {
\r
2238 return IsMotherPi0(imother,mcEvent);
\r
2243 if(mctrack->IsA() == AliAODMCParticle::Class()) {
\r
2244 AliAODMCParticle *mctrackaod = NULL;
\r
2245 if(!(mctrackaod = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
\r
2247 Int_t imother = mctrackaod->GetMother();
\r
2248 if(imother < 0) return -1;
\r
2249 AliAODMCParticle *mothertrack = NULL;
\r
2250 if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
\r
2252 Int_t pdg = mothertrack->GetPdgCode();
\r
2253 if(TMath::Abs(pdg) == 111) return imother;
\r
2254 if(TMath::Abs(pdg) == 11) {
\r
2255 return IsMotherPi0(imother,mcEvent);
\r
2263 //____________________________________________________________________________________________________________
\r
2264 Int_t AliAnalysisTaskHFEFlow::IsMotherC(Int_t tr, AliMCEvent* mcEvent) {
\r
2267 // Return the lab of signal mother or -1 if not signal
\r
2270 if(tr < 0) return -1;
\r
2271 AliVParticle *mctrack = mcEvent->GetTrack(tr);
\r
2273 if(mctrack->IsA() == AliMCParticle::Class()) {
\r
2274 AliMCParticle *mctrackesd = NULL;
\r
2275 if(!(mctrackesd = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
\r
2276 TParticle *particle = 0x0;
\r
2277 particle = mctrackesd->Particle();
\r
2279 if(!particle) return -1;
\r
2280 Int_t imother = particle->GetFirstMother();
\r
2281 if(imother < 0) return -1;
\r
2282 AliMCParticle *mothertrack = NULL;
\r
2283 if(!(mothertrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
\r
2284 TParticle * mother = mothertrack->Particle();
\r
2285 if(!mother) return -1;
\r
2287 Int_t pdg = mother->GetPdgCode();
\r
2288 if((TMath::Abs(pdg)==411) || (TMath::Abs(pdg)==421) || (TMath::Abs(pdg)==431) || (TMath::Abs(pdg)==4122) || (TMath::Abs(pdg)==4132) || (TMath::Abs(pdg)==4232) || (TMath::Abs(pdg)==43320)) return imother;
\r
2289 if(TMath::Abs(pdg) == 11) {
\r
2290 return IsMotherC(imother,mcEvent);
\r
2295 if(mctrack->IsA() == AliAODMCParticle::Class()) {
\r
2296 AliAODMCParticle *mctrackaod = NULL;
\r
2297 if(!(mctrackaod = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
\r
2299 Int_t imother = mctrackaod->GetMother();
\r
2300 if(imother < 0) return -1;
\r
2301 AliAODMCParticle *mothertrack = NULL;
\r
2302 if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
\r
2304 Int_t pdg = mothertrack->GetPdgCode();
\r
2305 if((TMath::Abs(pdg)==411) || (TMath::Abs(pdg)==421) || (TMath::Abs(pdg)==431) || (TMath::Abs(pdg)==4122) || (TMath::Abs(pdg)==4132) || (TMath::Abs(pdg)==4232) || (TMath::Abs(pdg)==43320)) return imother;
\r
2306 if(TMath::Abs(pdg) == 11) {
\r
2307 return IsMotherC(imother,mcEvent);
\r
2315 //____________________________________________________________________________________________________________
\r
2316 Int_t AliAnalysisTaskHFEFlow::IsMotherB(Int_t tr, AliMCEvent* mcEvent) {
\r
2319 // Return the lab of signal mother or -1 if not signal
\r
2322 if(tr < 0) return -1;
\r
2323 AliVParticle *mctrack = mcEvent->GetTrack(tr);
\r
2325 if(mctrack->IsA() == AliMCParticle::Class()) {
\r
2326 AliMCParticle *mctrackesd = NULL;
\r
2327 if(!(mctrackesd = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
\r
2328 TParticle *particle = 0x0;
\r
2329 particle = mctrackesd->Particle();
\r
2331 if(!particle) return -1;
\r
2332 Int_t imother = particle->GetFirstMother();
\r
2333 if(imother < 0) return -1;
\r
2334 AliMCParticle *mothertrack = NULL;
\r
2335 if(!(mothertrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
\r
2336 TParticle * mother = mothertrack->Particle();
\r
2337 if(!mother) return -1;
\r
2339 Int_t pdg = mother->GetPdgCode();
\r
2340 if((TMath::Abs(pdg)==511) || (TMath::Abs(pdg)==521) || (TMath::Abs(pdg)==531) || (TMath::Abs(pdg)==5122) || (TMath::Abs(pdg)==5132) || (TMath::Abs(pdg)==5232) || (TMath::Abs(pdg)==53320)) return imother;
\r
2341 if(TMath::Abs(pdg) == 11) {
\r
2342 return IsMotherB(imother,mcEvent);
\r
2347 if(mctrack->IsA() == AliAODMCParticle::Class()) {
\r
2348 AliAODMCParticle *mctrackaod = NULL;
\r
2349 if(!(mctrackaod = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
\r
2351 Int_t imother = mctrackaod->GetMother();
\r
2352 if(imother < 0) return -1;
\r
2353 AliAODMCParticle *mothertrack = NULL;
\r
2354 if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
\r
2356 Int_t pdg = mothertrack->GetPdgCode();
\r
2357 if((TMath::Abs(pdg)==511) || (TMath::Abs(pdg)==521) || (TMath::Abs(pdg)==531) || (TMath::Abs(pdg)==5122) || (TMath::Abs(pdg)==5132) || (TMath::Abs(pdg)==5232) || (TMath::Abs(pdg)==53320)) return imother;
\r
2358 if(TMath::Abs(pdg) == 11) {
\r
2359 return IsMotherB(imother,mcEvent);
\r
2367 //____________________________________________________________________________________________________________
\r
2368 Int_t AliAnalysisTaskHFEFlow::IsMotherEta(Int_t tr, AliMCEvent* mcEvent) {
\r
2371 // Return the lab of pi0 mother or -1 if not pi0
\r
2374 if(tr < 0) return -1;
\r
2375 AliVParticle *mctrack = mcEvent->GetTrack(tr);
\r
2377 if(mctrack->IsA() == AliMCParticle::Class()) {
\r
2378 AliMCParticle *mctrackesd = NULL;
\r
2379 if(!(mctrackesd = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
\r
2380 TParticle *particle = 0x0;
\r
2381 particle = mctrackesd->Particle();
\r
2383 if(!particle) return -1;
\r
2384 Int_t imother = particle->GetFirstMother();
\r
2385 if(imother < 0) return -1;
\r
2386 AliMCParticle *mothertrack = NULL;
\r
2387 if(!(mothertrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
\r
2388 TParticle * mother = mothertrack->Particle();
\r
2389 if(!mother) return -1;
\r
2391 Int_t pdg = mother->GetPdgCode();
\r
2392 if(TMath::Abs(pdg) == 221) return imother;
\r
2393 if(TMath::Abs(pdg) == 11) {
\r
2394 return IsMotherEta(imother,mcEvent);
\r
2399 if(mctrack->IsA() == AliAODMCParticle::Class()) {
\r
2400 AliAODMCParticle *mctrackaod = NULL;
\r
2401 if(!(mctrackaod = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(tr))))) return -1;
\r
2403 Int_t imother = mctrackaod->GetMother();
\r
2404 if(imother < 0) return -1;
\r
2405 AliAODMCParticle *mothertrack = NULL;
\r
2406 if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(TMath::Abs(imother))))) return -1;
\r
2408 Int_t pdg = mothertrack->GetPdgCode();
\r
2409 if(TMath::Abs(pdg) == 221) return imother;
\r
2410 if(TMath::Abs(pdg) == 11) {
\r
2411 return IsMotherEta(imother,mcEvent);
\r